How to output photolysis rates?

[PaulHamer]

I am currently running the 4.10 open source version of EMEP. In order to help diagnose the chemical behaviour in the model I would like to look at the photolysis rates in the model. I am specifically interested in JNO2. Ideally I would like to be able to plot the photolysis rates on an hourly basis. I also understand it can be difficult to output variables from the model if no specific diagnostic already exists in the output routines.

So far I just tried to crudely print the JNO2 photolysis rates. I have figured out where the rates are read in from the look up tables and how they get modified to consider clouds in the subroutine DefPhotolysis_ml.f90. I have added a print statement at the end of DefPhotolysis_ml.f90 (see below) and I have been able to analyse the rates that get printed.

          if(MasterProc) write(*,*) "JNO2 KCHEMTOP",rcphot(IDNO2,KCHEMTOP)
          if(MasterProc) write(*,*) "JNO2 KMAX_MID",rcphot(IDNO2,KMAX_MID)

The problem I have is that I have no way of figuring out which time the outputted rates correspond to. It seems that within runchem, that setup_phot gets called many times each hour. I have tried some crude analysis techniques on the outputted JNO2 values I printed to try and separate out the hourly variations by filtering out JNO2 values that repeat, but the results are not satisfying. I can see the dynamic range of JNO2 now, which is helpful, but I was looking for something a little more detailed, and perhaps have the chance to plot the variability of JNO2 for a specific day in the year, e.g., June 21st or December 21st.

Another way of doing this might be to "hack" the lookup tables for the information I need. The problem I have here is that I am not sure how to figure out the indexing within the tables to be able to extract a diurnal profile of JNO2 for a specific day.

Any help you could give to solve this issue I have would be greatly appreciated. Thanks in advance.

[avaldebe]

@JanEiofJonson please answer this

[JanEiofJonson]

Regarding hourly output of photolysis rates

As you have already found out photolysis rates are are pre-calculated for solar zenith angles 0 - 90. In the model we then read in the photolysis rates from the solar zenith and weighting between clear sky and cloudiness.

In order to print out the hourly photolysis rates you have to add it in Derived/My_Derived_ml

I have been on vacation + a meeting, and I need to do some extra work to give more advice on this.

[mifads]

Hopefully the "Derived" scheme will get easier in future, though we keep finding that all flexible solutions become complex quite quickly

If you just need data for individual points, a quick hack might be to print out from the Solver_ml, since there the model knows about rcphot(IDNO2,k). To activate the debug_flag, set DEBUG%RUNCHEM to be T in config_emep.nml, and set the DEBUG%IJ coordinates as you want. Then in Solver_ml, one could add something like:

  if ( debug_flag ) then
   call datewrite('JNO2', IDNO2, [ rcphot(IDNO2,20), rcphot(IDNO2,KCHEM_TOP) ] )  ! or : for k
  end if

towards the end, maybe even just before 'end subroutine chemistry'. The output will be to your standard output, but use 'grep JNO2' on that to get just these values. I haven't tested, but something like this should work.

(The photolysis scheme is getting pretty old, so we would be very interested to hear about any comparisons with measurements or other schemes.)

Dave

[gitpeterwind]

If you want values defined each hour, section 3.2 of the user guide shows how to output any value defined in the model. The values are then written in the standard netCDF output files. This is a rv4_15 feature though.

[JanEiofJonson]

I am using 4_14, and not 4_15. (mainly because I have been/still are running for a project testing the effects of different emissions and don't want to change version right now.

Will it still work for this version?

[JanEiofJonson]

To print out I made the following changes (as suggested by Peter Wind)

In DefPhotolysis_ml.f90

Add special3d

  use MetFields_ml       , only : cc3d,cc3dmax,z_bnd, special3d

Summing up at the end of the subroutine

           special3d(i,j,KCHEMTOP:KMAX_MID,1)= 0. 
           special3d(i,j,KCHEMTOP:KMAX_MID,1)=rcphot(IDNO2,KCHEMTOP:KMAX_MID) 

          end if   !  end izen <  90 (daytime)  test

    end subroutine setup_phot

In MetFields_ml.f90

Nspecial3d has to be set to 1 (originally usually 0)

  integer, parameter ::Nspecial3d = 1

In config_Outputs**.nml

Include in the list of output in the nml file

 'PhotNO2       ',  'MET3D', 'special3d1 ', '-', 's-1  ', -99,  -99,      F,   1.0, T, 'YMH',

[PaulHamer]

Apologies for the long silence. I have been busy on other projects until now.

I have quickly tried out your suggested implementation. Unfortunately, I run into a problem in the compilation. I have traced this to the fact that I am using the rv4.10 opensource version (see my original query), and I see that Nspecial3d and special3d were not yet implemented in this version of the code. I guess I should move to using rv4.14 in that case? At this stage we are using the photolysis rates for diagnostic purposes for interpreting the model output, so any changes in the resulting concentrations in the newer version won't affect this analysis.

[PaulHamer]

After a bit of work I have figured out that you can extract the essential pieces of code from the rv4.15 version of MetFields.f90 and implement them in the rv4.10 version of MetFields.f90. Please find attached the version of the MetFields.f90 and DefPhotolysis_ml.f90 file that I have successfully been able to use to output the photolysis fields. Although, note, in addition to the code changes recommended by Jon, it is also necessary to add the following lines of code in DefPhotolysis_ml.f90 in order to output photolysis rates during the night.

    if ( Grid%izen > 90 ) then     ! Photolysis rates zero when the sun is
                                                 ! below the horizon
         rcphot(:,:) = 0.0
         special3d(i,j,KCHEMTOP:KMAX_MID,1)= 0.d0
         special3d(i,j,KCHEMTOP:KMAX_MID,1)=rcphot(IDNO2,KCHEMTOP:KMAX_MID)

    else !! (izen < 90)  -- sun above horizon

With the code as recommended by Jon, you get photolysis rates >0 during the night since the code lines to creates special3d do not get called during the night time.

ps. I cannot upload the code files without changing the files to .txt files. DefPhotolysis_ml.f90.txt MetFields_ml.f90.txt

[avaldebe]

@PaulHamer It seems that your question was satisfactory answered, so I'll close this issue. If it was not, please re-open it.

[gitpeterwind]

In the version rv4_17, the NO2 photolysis rate can simply be obtained with a line in the config file using:

OutputMisc=
  'J(NO2)'         ,'USET','D3_J(NO2)'  ,'photorate','1/s' ,-99,-99,F,1.0,T,'YM',