Difference between revisions of "UKCA Chemistry and Aerosol Tutorial 6"

From UKCA
Line 1: Line 1:
 
[[UKCA Chemistry and Aerosol Tutorials | Back to UKCA Chemistry and Aerosol Tutorials]]
  
[[UKCA Chemistry and Aerosol Tutorials | Back to UKCA Chemistry and Aerosol Tutorials]]
   
−
==What you will learn in this Tutorial==
 +
==What you will learn in this tutorial==
   
−
During this tutorial you will learn how to make new UM ancillary files. Then you will learn how to add new emissions into UKCA so that they emit into one of your new tracers.
 +
During this tutorial you will learn how UKCA specifies different chemical reactions. You will then add a new reaction involving the new tracers that you have added.
   
  +
==Task 6.1: Add a bimolecular reaction==
−
At the end of the [[UKCA & UMUI Tutorial 4 | previous tutorial]] you will now know how to create new tracers for use by UKCA. However, after completing the tasks, your tracers will still be empty, as nothing has been put into them. This tutorial will teach you how to create an emissions ancillary file that the UM will read, and that you can then tell UKCA to use and emit into your tracer(s).
  
   
  +
<span style="color:green">'''TASK 6.1:''' You should now add in the bimolecular reaction of '''ALICE''' with '''OH''' to form '''BOB''' and a secondary organic compound (labelled in UKCA as '''SEC_ORG'''). This reaction is given by:</span>f
−
This tutorial will go through the steps needed to make an emission into a tracer which UKCA does not currently emit into. The steps in making the ancillary file will be the same for a species which is currently emitted into, although in this simplier case you would not need to make any code changes.
  
  +
<math>
  +
\textrm{ALICE} + \textrm{OH} \longrightarrow \textrm{BOB} + \textrm{SEC\_ORG}
  +
</math>
   
  +
{| border="1"
−
During this tutorial you will be tasked with making a new emissions ancillary file, and adding it in to one of your new tracers.
  
  +
! Parameter || Value
  +
|-
  +
| <math>k_{0}</math> || 2.70E-11
  +
|-
  +
| <math>\alpha</math> || 0.00
  +
|-
  +
| <math>\beta</math> || -390.00
  +
|}
   
  +
'''Note:''' If you were unable to successfully complete [[UKCA Chemistry and Aerosol Tutorial 5#Task 5.2: make the required code changes to add your emission into UKCA| Task 5.2]], then please take a copy of the '''f''' job from the Tutorial experiment (''Tutorial: solution to Task 5.2 - adding new chemical emissions in UKCA'') and work from there, as this will allow you to only make the required changes. Please also make a new branch and merge-in branch '''branch''' at revision number '''rev''' to allow you to proceed.
−
==Task 6.1: Create a new emissions file and use it in your job==
  
   
  +
==Adding new Chemical Reactions==
−
<span style="color:green">'''Task 6.1:''' In the <pre>/work/n02/n02/ukca/Tutorial/vn8.4/Task6.1</pre> directory on ARCHER, or the <pre>/projects/ukca/Tutorial/vn8.4/Task6.1</pre> directory on MONSooN, there is the file '''Emissions_of_ALICE.nc''' which is a 0.5x0.5 degree resolution surface emission field. You should regrid this file to '''N96''', and then make a new surface emissions ancillary file with this as slot '''316'''. You should then use this new file, and the corresponding user STASHmaster file, in your UMUI job.</span>
  
   
  +
UKCA currently uses two different methods of defining the chemical reactions solved in the model. The first is a backward Euler solver, and is used for the ''RAQ'' and ''StdTrop'' chemistry schemes where the solver itself is created by a code-writer. The second makes use of the [http://www.atm.ch.cam.ac.uk/acmsu/asad/ ASAD chemical integration software package], and is used for the ''CheT/TropIsop'', ''CheS/Strat'', and ''CheST/StratTrop'' chemistry schemes. ASAD can use many different solvers, although currently it uses symbolic Newton-Raphson solver. In this tutorial we will only consider the ASAD framework, as this is easily extended by a user.
−
'''Note:''' If you were unable to successfully complete [[UKCA Chemistry and Aerosol Tutorial 5#Task_5.2: Add these two new tracers to UKCA | Task 5.2]], then please take a copy of the '''e''' job from the Tutorial experiment (''Tutorial: solution to Task 5.2 - adding new chemical tracers to UKCA'') and work from there, as this will allow you to make only the changes required for this task. Please also make a new branch and merge-in branch '''branch''' at revision number '''rev''' to allow you to proceed.
  
   
  +
ASAD considers four different types of chemical reactions: bimolecular reactions, termolecular reactions, heterogeneous reactions, and photolysis reactions. To make changes and add reactions you will need to make changes to the UKCA source code which can be found in
−
==Part 1. Making a new Emissions Ancillary File==
  
   
  +
vn8.4_<span style="color:blue">your_branch_name</span>/src/atmosphere/UKCA
−
The UM uses its own format to read-in initial and emissions data, the ''ancillary file''. UKCA makes use of these files for the surface and aircraft emissions, and these files can easily be made up from netCDF data using the [http://cms.ncas.ac.uk/documents/xancil/ Xancil] program. However, before we can use Xancil to create our emissions ancillary, we may first need to use [http://badc.nerc.ac.uk/help/software/xconv/ Xconv] to regrid the raw data to the correct resolution of the UM configuration that you are using.
  
   
  +
During this tutorial you will be tasked with adding a new reaction into your branch.
−
===Regridding data with Xconv===
  
   
  +
==Biomolecular Reactions==
−
We are using a model at N96L85 resolution. In the horizontal this is 1.875 degrees by 1.25 degrees. There are 85 vertical levels.
 
   
  +
For most bimolecular reactions, it is sufficient to provide the <math>k_{0}</math>, <math>\alpha</math>, and <math>\beta</math> coefficients that are used to compute the rate coefficient <math>k</math> from the Arrhenius expression
−
The '''N96''' UM grid
 
   
  +
<math>
−
# Has 192 points in longitude (x) and 145 points in latitude (y)
  
  +
k = k_{0} \left(\frac{T}{300}\right)^{\alpha} \textrm{exp} \left(\frac{-\beta}{T}\right)
−
# Starts at 0.0 longitude with a spacing of 1.875 degrees
  
  +
</math>
−
# Starts at -90.0 latitude (i.e. the South Pole) with a spacing of 1.25 degrees
  
   
  +
===Bimolecular Reaction Definition===
−
====Horizontal Regridding====
  
   
  +
The bimolecular reactions are defined in the '''ukca_chem_<span style="color:blue">scheme</span>.F90''' routines using the '''ratb_t''' Fortran type specification, and are held in arrays. At the end of this routine the '''ratb_defs_<span style="color:blue">scheme</span>''' array is created from these, and if that scheme is selected in UKCA these reactions are copied across into the master '''ratb_defs''' array.
−
[[Image:Xconv_Trans_N96.png|thumb|right|Figure 1: The Xconv ''Trans'' window.]]
  
−
Horizontal regridding in Xconv is straight-forward. First, open your dataset in Xconv by
  
   
  +
The format of this '''ratb_t''' type is
−
xconv -i <span style="color:blue">file</span>.nc
  
   
  +
ratb_t('Reactant 1','Reactant 2','<span style="color:blue">Product 1 </span>','<span style="color:red">Product 2 </span>','<span style="color:green">Product 3 </span>',&
−
and double-click on the field that you want to regrid and then click the '''Trans''' button at the far top right. This will show a window similar to Figure 1 (which in fact shows the default settings for a N96 field).
  
  +
'<span style="color:purple">Product 4 </span>', <math>k_{0}</math>, <math>\alpha</math>, <math>\beta</math>, <span style="color:blue">Fraction of Product 1 produced</span>, <span style="color:red">Fraction of Product 2 produced</span>, &
  +
<span style="color:green">Fraction of Product 3 produced</span>, <span style="color:purple">Fraction of Product 4 produced</span>), &
   
  +
If fractional products are not required for a reaction, then the ''fraction of each product'' formed should be set to 0.000. If fractional products are required for any one of the products then the fraction of each product formed should be set to its correct value.
−
As we are regridding an emission we need to select '''area weighted interpolation''' as we need to conserve the total amount of quantity emitted. Then scroll-down to the boxes at the end and enter the following values:
  
   
  +
The specifications of the individual reactions are done as, e.g.
−
* Number of columns = 192
  
−
* First longitude = 0.000000
  
−
* Column spacing = 1.875000
  
−
* Number of rows = 145
  
−
* First latitude = -90.000000
  
−
* Row spacing = 1.250000
  
   
  +
ratb_t('O3 ','C5H8 ','HO2 ','OH ',' ',& ! B133
−
which match up with the grid definition above, and ensure that
  
  +
' ', 3.33E-15, 0.00, 1995.00, 0.750, 0.750, 0.000, 0.000), & ! B133 IUPAC2007*
  +
...
  +
ratb_t('OH ','C5H8 ','ISO2 ',' ',' ',& ! B144
  +
' ', 2.70E-11, 0.00, -390.00, 0.000, 0.000, 0.000, 0.000), & ! B144 IUPAC2009
  +
...
  +
ratb_t('OH ','HCl ','H2O ','Cl ',' ',& ! B159
  +
' ', 1.80E-12, 0.00, 250.00, 0.000, 0.000, 0.000, 0.000), & ! B159 JPL2011
   
  +
The first reaction in these examples takes its kinetic data from [http://www.iupac-kinetic.ch.cam.ac.uk/ IUPAC]. Going to this website, this reaction is defined [http://www.iupac-kinetic.ch.cam.ac.uk/datasheets/xhtml/HOx_VOC8_HO_CH2C%28CH3%29CHCH2%28isoprene%29.xhtml_mathml.xml here]. The second reaction above takes its kinetic data from [http://jpldataeval.jpl.nasa.gov/ NASA's Jet Propulsion Laboratory]. The rate for this can be found on page 1-19 of the [http://jpldataeval.jpl.nasa.gov/pdf/JPL%2010-6%20Final%2015June2011.pdf JPL2011 document]. When adding new reactions you will need to increment the size of the array holding the <tt>ratb_t</tt> type.
−
* Pole longitude = 0.000000
  
−
* Pole latitude = 90.000000
  
   
  +
To add new bimolecular reactions you will need to append equivalent lines for the new reactions to the end of the '''ratb_defs_<span style="color:blue">scheme</span>''' array (increasing the array sizes accordingly). If there is a reaction that is an exception to the general Arrhenius equation then special code needs to be placed in the '''asad_bimol.F90''' routine, which is held in the <tt>UKCA/</tt> source-code directory.
−
Now click '''Apply'''. You should be given a message similar to
  
   
  +
===Increase the size of JPBK (and JPNR)===
−
Area weighted interpolation from 720x360 Regular grid to 192x145 Regular grid
  
   
  +
As well as adding these reactions to the ''ukca_chem_<span style="color:blue">scheme</span>.F90'' routine (and incrementing the size of the arrays in that routine accordingly, you will also need to increase the values of two parameters that UKCA needs. These are
−
in the dialogue window. You can now extract this data as a new netCDF file (you cannot re-save a file in Xconv) by putting a name in the '''output file name''' box and clicking '''convert'''.
  
   
  +
* '''JPBK''' is the number of bimolecular reactions
−
====Vertical Regridding====
  
  +
* '''JPNR''' is the total number of reactions
   
  +
These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the <code>&RUN_UKCA</code> namelist in '''CNTLATM'''. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the ''CNTLATM'' file in your <tt>$HOME/umui_jobs/<span style="color:blue">jobid</span></tt> directory.
−
It is not possible to vertically regrid data using Xconv. You will need to do this in another way. If you need to vertically regrid data, and am unsure of the best way, please contact [[User:Nla27 | Luke Abraham]] for advice.
 
   
  +
==Termolecular Reactions==
−
Remember that
  
   
  +
As well as defining reactions involving a third body, the termolecular rate definition can also be used to define unimolecular reactions.
−
* Emissions data must be re-gridded in a mass conserving way, so you will probably need to integrate the field on one grid and then decompose it again on the new grid.
  
−
* Tracer data can be fitted to the profile of the field on the old grid.
  
   
  +
The pressure and temperature dependent rate, <math>k</math>, of a termolecular reaction is given by
−
===Choosing a STASH slot for your new emission(s)===
  
   
  +
<math>
−
To make a new ancillary file for your new emission(s), you should first decide on a STASH item for it/them. Currently UKCA makes use of the '''user single-level ancillary file''' and '''user multi-level ancillary file''' which uses STASH section 0 items 301-320 (single-level) and 321-340 (multi-level). What these numbers correspond to is set in the file '''ukca_setd1defs.F90''', as well as in the user STASHmaster file associated with the job you are using (which can be found in '''Model Selection &rarr; Atmosphere &rarr; STASH &rarr; User-STASHmaster files. Diags, Progs & Ancills''').
 
  +
k = \left(\frac{k_{0}\left[M\right]}{1+k_{0}\left[M\right]/k_{\infty}}\right)F_{c}^{\left(1+\left[\textrm{log}_{10}\left(\frac{k_{0}\left[M\right]}{k_{\infty}}\right)\right]^{2}\right)^{-1}}
  +
</math>
   
  +
where the low pressure rate constant <math>k_{0}</math> is given by
−
====Listing of emissions from STASH====
  
   
  +
<math>
−
{| border="1"
  
  +
k_{0} = k_{1} \left(\frac{T}{300}\right)^{{\alpha}_{1}} \textrm{exp} \left(\frac{-{\beta}_{1}}{T}\right)
−
! Stash code || Emission
  
  +
</math>
−
|-
  
−
| 301 || NOx surf emissions
 
−
|-
 
−
| 302 || CH4 surf emissions
 
−
|-
 
−
| 303 || CO surf emissions
 
−
|-
 
−
| 304 || HCHO surf emissions
 
−
|-
 
−
| 305 || C2H6 surf emissions
 
−
|-
 
−
| 306 || C3H8 surf emissions
 
−
|-
 
−
| 307 || ME2CO surf emissions
 
−
|-
 
−
| 308 || MECHO surf emissions
 
−
|-
 
−
| 309 || C5H8 surf emissions
 
−
|-
 
−
| 310 || BC fossil fuel surf emissions
 
−
|-
 
−
| 311 || BC biofuel surf emissions
 
−
|-
 
−
| 312 || OC fossil fuel surf emissions
 
−
|-
 
−
| 313 || OC biofuel surf emissions
 
−
|-
 
−
| 314 || Monoterpene surf emissions
 
−
|-
 
−
| 315 || NVOC surf emissions
 
−
|-
 
−
| 322 || BC BIOMASS 3D EMISSION
 
−
|-
 
−
| 323 || OC BIOMASS 3D EMISSION
 
−
|-
 
−
| 340 || NOX AIRCRAFT EMS IN KG/S/GRIDCELL
 
−
|}
  
   
  +
and the high pressure rate constant <math>k_{\infty}</math> is given by
−
====Code in ukca_setd1defs.F90====
  
   
  +
<math>
−
The species emitted are set in two places, firstly in the definition of an array called '''em_chem_spec''' which is scheme specific, and secondly in a block of code which searches through the ''em_chem_spec'' array and assigns a STASH number to it (as defined by the list above).
  
  +
k_{\infty} = k_{2} \left(\frac{T}{300}\right)^{{\alpha}_{2}} \textrm{exp} \left(\frac{-{\beta}_{2}}{T}\right)
  +
</math>
   
  +
===Termolecular Reaction Definition===
−
For example, for the ''CheST/StratTrop'' chemistry (not using aerosol chemistry), em_chem_spec is set to
  
   
  +
The termolecular reactions are defined in the '''ukca_chem_<span style="color:blue">scheme</span>.F90''' routines using the '''ratt_t''' Fortran type specification, and are usually held in one single array (there are not usually enough reactions to require splitting the reactions over several arrays).
−
em_chem_spec = &
  
−
(/'NO ','CH4 ','CO ','HCHO ', &
  
−
'C2H6 ','C3H8 ','Me2CO ','MeCHO ', &
  
−
'C5H8 ','NO_aircrft'/)
  
   
  +
To format of this '''ratt_t''' type is
−
This can be found in the '''Stratospheric Chemistry''' section, controlled by the IF block where <code>(L_ukca_strattrop .AND. .NOT. L_ukca_achem)</code>. If you are using aerosol chemistry then the number of emissions are increased accordingly.
  
   
  +
ratt_t('Reactant 1','Reactant 2','<span style="color:blue">Product 1 </span>','<span style="color:red">Product 2 </span>', <math>f</math>, &
−
Further down the code there is this block of code:
  
  +
<math>k_{1}</math>, <math>{\alpha}_{1}</math>, <math>{\beta}_{1}</math>, <math>k_{1}</math>, <math>{\alpha}_{1}</math>, <math>{\beta}_{1}</math>, <span style="color:blue">Fraction of Product 1 produced</span>, <span style="color:red">Fraction of Product 2 produced</span>), &
   
  +
and as in <tt>ratb_t</tt>, where the fraction of a product should be set to 0.000 if this functionality does not need to be used.
−
J = n_use_tracers
  
−
IF (n_chem_emissions+n_3d_emissions+n_mode_emissions > 0) THEN
  
−
DO i=1,n_chem_emissions + n_3d_emissions
  
−
UkcaD1Codes(J+i)%section = 0
  
−
UkcaD1Codes(J+i)%item = n_emiss_first+i-1 ! trop chemistry
  
−
UkcaD1Codes(J+i)%len_dim1 = row_length ! uses stash codes
  
−
UkcaD1Codes(J+i)%len_dim2 = rows ! 301-309 for
  
−
UkcaD1Codes(J+i)%required = .true. ! surface emissions
  
−
UkcaD1Codes(J+i)%prognostic = .true. ! from Section 0
  
−
! Special cases, emissions already available in UM
  
−
IF (em_chem_spec(i)(1:7) == 'SO2_low') THEN
  
−
UkcaD1Codes(J+i)%item = 58
  
−
IF (.NOT. L_SO2_SURFEM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='SO2 surface emissions from UM are not flagged'
  
−
errcode=58
  
−
 
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:7) == 'SO2_nat') THEN
  
−
UkcaD1Codes(J+i)%item = 121
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
IF (.NOT. L_SO2_NATEM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='SO2 natural emissions from UM are not flagged'
  
−
errcode=121
  
−
 
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:8) == 'SO2_high') THEN
  
−
UkcaD1Codes(J+i)%item = 126
  
−
IF (.NOT. L_SO2_HILEM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='SO2 high-level emissions are not flagged'
  
−
errcode = UkcaD1Codes(J+i)%item
  
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:3) == 'NH3') THEN
  
−
UkcaD1Codes(J+i)%item = 127
  
−
IF (.NOT. L_NH3_EM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='NH3 surface emissions from UM are not flagged'
  
−
errcode = UkcaD1Codes(J+i)%item
  
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i) == 'BC_fossil ') THEN
  
−
UkcaD1Codes(J+i)%item = 310
  
−
ELSEIF (em_chem_spec(i) == 'BC_biofuel') THEN
  
−
UkcaD1Codes(J+i)%item = 311
  
−
ELSEIF (em_chem_spec(i) == 'OC_fossil ') THEN
  
−
UkcaD1Codes(J+i)%item = 312
  
−
ELSEIF (em_chem_spec(i) == 'OC_biofuel') THEN
  
−
UkcaD1Codes(J+i)%item = 313
  
−
ELSEIF (em_chem_spec(i) == 'Monoterp ') THEN
  
−
UkcaD1Codes(J+i)%item = 314
  
−
ELSEIF (em_chem_spec(i) == 'NVOC ') THEN
  
−
UkcaD1Codes(J+i)%item = 315
  
−
ELSEIF (em_chem_spec(i) == 'BC_biomass') THEN
  
−
UkcaD1Codes(J+i)%item = 322
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ELSEIF (em_chem_spec(i) == 'OC_biomass') THEN
  
−
UkcaD1Codes(J+i)%item = 323
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ELSEIF (em_chem_spec(i) == 'SO2_biomas') THEN
  
−
UkcaD1Codes(J+i)%item = 324
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ELSEIF (em_chem_spec(i)(1:3) == 'DMS') THEN
  
−
UkcaD1Codes(J+i)%section = 17
  
−
UkcaD1Codes(J+i)%item = 205
  
−
UkcaD1Codes(J+i)%prognostic = .false.
  
−
IF (.NOT. L_DMS_EM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='DMS surface emissions from UM are not flagged'
  
−
errcode = UkcaD1Codes(J+i)%section*1000 + &
  
−
UkcaD1Codes(J+i)%item
  
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:7) == 'NO_airc') THEN
  
−
UkcaD1Codes(J+i)%item = 340
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ENDIF
  
−
ENDDO
  
−
ENDIF
  
   
  +
The <math>f</math> value is used to define the <math>F_{c}</math> value by
−
This block of code is rather complicated, but what it essentially means is that for the STASH codes 301-309, the emissions are assumed to be in the order of the species in ''em_chem_spec'', but for the other emissions the STASH numbers are explicitly defined. As you can see from the table above, for 2D (surface) emissions the numbers 301-315 are reserved, and for 3D emissions the numbers 322, 323, and 340 are reserved.
 
   
  +
<blockquote>
−
This means that if you are adding in a new surface emission(s) it is best to use the slots 316-320, unless you need more than 5 slots. For 3D emissions you have more leeway.
  
  +
If <math>f < 1.0</math> then <math>F_{c} = f</math><br/>
  +
else <math>F_{c} = \textrm{exp}\left(-T/f\right)</math>
  +
</blockquote>
   
  +
as <math>F_{c}</math> may or may not be highly temperature dependent.
−
====Emissions STASHmaster File====
  
   
  +
Examples of these reactions are
−
Now that you have selected your slot(s), you need to create a new STASH specification for it/them. The easiest way to do this is to copy the existing user STASHmaster file that defines your current (possible) emissions, and extend that. This is found in '''Model Selection &rarr; Atmosphere &rarr; STASH &rarr; User-STASHmaster files. Diags, Progs & Ancills'''. This will be a different STASHmaster file to the one that contains the UKCA tracers, and e.g. in the UKCA Tutorial job, is called '''emiss_TCMIM_Aero.presm'''.
  
   
  +
ratt_t('N2O5 ','m ','NO2 ','NO3 ', 0.3, & ! T023
−
This contains entries like
  
  +
1.30E-03, -3.50, 11000.00, 9.70E+14, 0.10, 11080.00, 0.000, 0.000), & ! T023 IUPAC 2002
  +
ratt_t('NO ','NO ','NO2 ','NO2 ', 0.0, & ! T024
  +
3.30E-39, 0.00, -530.00, 0.00E+00, 0.00, 0.00, 0.000, 0.000) & ! T024 IUPAC 2001
   
  +
To add new termolecular reactions you will need to append equivalent lines for the new reactions to the end of the '''ratt_defs_<span style="color:blue">scheme</span>''' array (increasing the array sizes accordingly).
−
#
  
−
1| 1 | 0 | 301 |NOx surf emissions |
  
−
2| 2 | 0 | 1 | 1 | 5 | -1 | -1 | 0 | 0 | 0 | 0 |
  
−
3| 000000000000000000000000000000 | 00000000000000000001 | 3 |
  
−
4| 1 | 0 | -99 -99 -99 -99 -99 -99 -99 -99 -99 -99 |
  
−
5| 0 | 531 | 0 | 129 | 0 | 0 | 0 | 0 | 0 |
  
−
#
 
−
 
−
#
  
−
1| 1 | 0 | 340 |NOX AIRCRAFT EMS IN KG/S/GRIDCELL |
  
−
2| 2 | 0 | 1 | 1 | 2 | 10 | 11 | 0 | 0 | 0 | 0 |
  
−
3| 000000000000000000000000000000 | 00000000000000000001 | 3 |
  
−
4| 1 | 0 | -99 -99 -99 -99 -99 -99 -99 -99 -99 -99 |
  
−
5| 0 | 520 | 20 | 65 | 0 | 0 | 0 | 9999 | 0 |
  
−
#
  
   
  +
===Increase the size of JPTK (and JPNR)===
−
The top entry (''NOx surf emissions'') defines a 2D field, and the other (''NOX AIRCRAFT EMS IN KG/S/GRIDCELL'') defines a 3D field. You can see that there are differences in the numbers (other than the 301/340 item number) in various places in these specifications, which effectively (in this instance) tell STASH if the field is 2D or 3D.
 
   
  +
As with the bimolecular reactions, you will also need to increase the values of two parameters that UKCA needs. These are
−
Full details on what each of these numbers mean can be found in appendix 3 of ''Unified Model Documentation Paper C4'' which can be found on the [http://collab.metoffice.gov.uk/twiki/pub/Support/Umdp Met Office Collaboration Twiki (password required)].
  
   
  +
* '''JPTK''' is the number of termolecular reactions
−
You should copy either the 2D or 3D specification, depending on what type of emission you are adding in, and edit only the '''STASH item number''', '''name of field''', and the '''field code'''. These can be found here:
  
  +
* '''JPNR''' is the total number of reactions
   
  +
These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the <code>&RUN_UKCA</code> namelist in '''CNTLATM'''. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the ''CNTLATM'' file in your <tt>$HOME/umui_jobs/<span style="color:blue">jobid</span></tt> directory.
−
#
  
−
1| 1 | 0 | <span style="color:red">301</span> |<span style="color:red">NOx surf emissions</span> |
  
−
2| 2 | 0 | 1 | 1 | 5 | -1 | -1 | 0 | 0 | 0 | 0 |
  
−
3| 000000000000000000000000000000 | 00000000000000000001 | 3 |
  
−
4| 1 | 0 | -99 -99 -99 -99 -99 -99 -99 -99 -99 -99 |
  
−
5| 0 | <span style="color:red">531</span> | 0 | 129 | 0 | 0 | 0 | 0 | 0 |
  
−
#
 
   
  +
==Heterogeneous Reactions==
−
For the field code (531 above), if you were making a new emission at 316, it is advisable that you increase the current code above by 15 as well, e.g. ''546''.
  
   
  +
Heterogeneous reactions are those that occur on aerosol surfaces. There is no functional form defined for these reactions, with special code needed to be added for each case.
−
You should make these changes to your copy of the original emissions user STASHmaster file, save this, and then replace the original file in the UMUI's ''Model Selection &rarr; Atmosphere &rarr; STASH &rarr; User-STASHmaster files. Diags, Progs & Ancills'' table with your new file. Now go to '''Model Selection &rarr; Atmosphere &rarr; STASH &rarr; Initialisation of User Prognostics''' and scroll down the table until you find your new emission. Set the value in the '''Option''' column to '''2''' (''Initialise to User Ancillary File'').
  
   
  +
===Heterogeneous Reaction Definition===
−
As you need to make up a new ancillary file, you should copy this user STASHmaster file onto the supercomputer, as it will be needed by Xancil when you make the new ancillary file. It is advisable to put it in the same directory as the one containing your new emission(s) file.
  
   
  +
The heterogeneous reactions are defined in the '''ukca_chem_<span style="color:blue">scheme</span>.F90''' routines using the '''rath_t''' Fortran type specification, usually in one array.
−
===Using Xancil===
  
  +
To format of this '''rath_t''' type is
   
  +
rath_t('Reactant 1','Reactant 2','<span style="color:blue">Product 1 </span>','<span style="color:red">Product 2 </span>','<span style="color:green">Product 3 </span>',&
−
====Extract your current emissions====
  
  +
'<span style="color:purple">Product 4 </span>', <span style="color:blue">Fraction of Product 1 produced</span>, <span style="color:red">Fraction of Product 2 produced</span>, &
  +
<span style="color:green">Fraction of Product 3 produced</span>, <span style="color:purple">Fraction of Product 4 produced</span>), &
   
  +
i.e. there is no rate information provided. For reactions on PSCs special code has been added to the routines in '''ukca_hetero_mod.F90''', and for other reactions there is code in '''asad_hetero.F90'''. Examples of this type are
−
Due to current limitations of the UM, you can only have one user single-level and one multi-level file. If you want to add a new emissions field, you must include the existing emissions in your new file along with it.
  
   
  +
rath_t('ClONO2 ','H2O ','HOCl ','HONO2 ',' ', &
−
You can use Xconv to extract these fields. You must first find the location of the current file(s). To do this, go to '''Model Selection &rarr; Atmosphere &rarr; Ancillary data and input data files &rarr; Climatologies and potential climatologies''' and either open the '''User multi-level ancillary file & fields''' or '''User single-level ancillary file & fields''' panel. This will give two boxes giving the '''directory name or environment variable''' and the '''file name'''. This first box will contain an environment variable which sets the directory location. You can find the value of this environment variable in '''Model Selection &rarr; Input/Output Control and Resources &rarr; Time Convention and SCRIPT Environment Variables'''.
  
  +
' ', 0.000, 0.000, 0.000, 0.000), &
  +
...
  +
rath_t('SO2 ','H2O2 ','NULL0 ',' ',' ', & !HSO3+H2O2(aq)
  +
' ', 0.000, 0.000, 0.000, 0.000), &
   
  +
To add new heterogeneous reactions you will need to append equivalent lines for the new reactions to the end of the '''ratt_defs_<span style="color:blue">scheme</span>''' array (increasing the array sizes accordingly), before adding code to either '''ukca_hetero_mod.F90''' or '''asad_hetero.F90'''.
−
Once you have found the required file, extract all the fields within it to one of your own directories (ideally the one containing the netCDF file of your new emission(s)).
  
   
−
====Make a new emissions ancillary file====
 +
===Increase the size of JPHK (and JPNR)===
   
  +
As with the bimolecular and termolecular reactions, you will also need to increase the values of two parameters that UKCA needs. These are
−
Xancil is installed on both HECToR at
  
   
  +
* '''JPHK''' is the number of heterogeneous reactions
−
/work/n02/n02/hum/bin/xancil
  
  +
* '''JPNR''' is the total number of reactions
   
  +
These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the <code>&RUN_UKCA</code> namelist in '''CNTLATM'''. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the ''CNTLATM'' file in your <tt>$HOME/umui_jobs/<span style="color:blue">jobid</span></tt> directory.
−
and on MONSooN (the postproc03 machine) at
  
   
  +
==Photolysis Reactions==
−
/projects/um1/linux/bin/xancil
  
   
  +
These define a reaction where a chemical compound is broken down by photons. There is no functional form defined for this type of reaction. Instead, either (in the troposphere) input files are used to define the reaction rates for each species, while (in the stratosphere) on-line look-up tables are generated for the rates for each species, or separate photolysis codes, '''Fast-J''' or '''Fast-JX''', are used to interactively calculate the rate of reaction throughout the troposphere (for Fast-J) or the whole atmosphere (for Fast-JX). These interactive schemes are preferred as they take the effect of aerosols or clouds into account at each timestep, allowing for more feedbacks to be investigated. In the upper stratosphere there are some wavelength regions that Fast-JX does not consider, and so the 3D on-line look-up tables are also used for these regions.
−
You may already have this location in your PATH and so can just launch Xancil from the command line by typing <code>xancil</code>. When you do this it will load up the Xancil window, which is initially rather empty. You should click on the '''Xancil''' text in the top-left corner. This will give 4 options:
 
   
  +
===Tropospheric Off-Line Photolysis===
−
* Configuration
  
−
* Atmosphere Ancillary Files
  
−
* Ocean Ancillary Files
  
−
* Generalised Ancillary Files
  
   
  +
If Fast-JX is not being used, then the off-line two-dimensional (zonally average) tropospheric photolysis is used (for all schemes). It is based on the work of Hough (1988)[1] and Law ''et al'' (1998)[2].
−
When using UKCA you will need to make use of all of these options, with the exception of the ''Ocean Ancillary Files''.
 
   
  +
This scheme makes use of datafiles which define the reaction rate for a particular species (e.g. H2O2), or if no rate is known, a '''nil''' rate can be used. For UM 8.2 these files (in ASCII format) can be found in
−
'''You should view the ''[[Using Xancil]]'' page for more details on each of these sub-menus.'''
  
−  
−
Load your netCDF and user STASHmaster files into the '''Xancil &rarr; Configuration''' panel, and [[Using Xancil#Grid Configuration | define the vertical levels]] if you are making a 3D ancillary file.
 
 
 
  +
/work/n02/n02/hum/vn8.2/ctldata/UKCA/tropdata/photol
−
Click on the ''Atmosphere Ancillary Files'' text and either open the '''Multi-level User Fields''' or '''Single-level User Fields''' panels. You should now
  
   
  +
on HECToR, and in
−
* Set the output file name
  
−
* Define the dates. If the ones in the netCDf file are fine to use then you can use them, or you can specify them. For this latter option you can either use the dates from the [[Using Xancil#Grid Configuration | grid configuration panel]], or you can define them again here
  
−
* Enter the number of ancillary fields. This will need to be the total of the number of fields in the original ancillary file, plus the number of new emissions you are adding
  
−
* For each individual field you should select it with the mouse, then
  
−
** Define which netCDF file to use (it will remember the preference from the previous field)
  
−
** Enter the STASH code (the PP code should be filled in automatically)
  
−
** Either enter or confirm the variable name. Xancil takes the variable name from the ''field code'' defined in the STASHmaster file specification for each field. If there are multiple fields in the netCDF file with the same field code then Xancil may select the wrong one. If the name does not match a field code you will need to select it manually.
  
   
  +
/projects/um1/vn8.2/ctldata/UKCA/tropdata/photol
−
Once you have entered all the data for all fields you should
  
   
  +
on MONSooN. To use this scheme, in the UMUI go to '''Model Selection &rarr; Atmosphere &rarr; Model Configuration &rarr; UKCA Chemistry and Aerosols &rarr; PHOTO''' and click '''2D Photolysis Scheme'''. You will then need to give the location of the files (above). The code controlling this scheme is held in '''ukca_phot2d.F90'''.
−
* Use the '''Save/Save As''' button to save the job, as it may be needed later
  
−
* Click the '''Output Anc. Files''' button to create the new ancillary file. Any errors will appear in the '''Output messages''' window, or to the terminal.
  
   
  +
It is advised that this scheme is no longer used, and interactive photolysis should be used instead. For the ''CheS/Strat'' or ''CheST/StratTrop'' schemes, Fast-JX should be used as this covers the stratosphere as well as the troposphere.
−
===Use your new emissions file in the UMUI===
  
   
  +
'''References'''
−
Now that you have created your new emissions file you can use this in the UMUI. Go to '''Model Selection &rarr; Atmosphere &rarr; Ancillary data and input data files &rarr; Climatologies and potential climatologies''' and either open the '''User multi-level ancillary file & fields''' or '''User single-level ancillary file & fields''' panel, and set the '''Directory name or Environment Variable''' to the directory containing your new emissions file, and the '''file name''' to the name of your new file.
  
  +
# Hough, A. M.: The calculation of photolysis rates for use in global modelling studies, Tech. rep., UK Atomic Energy Authority, Harwell, Oxon., UK, 1988
  +
# Law, K., Plantevin, P., Shallcross, D., Rogers, H., Pyle, J., Grouhel, C., Thouret, V., and Marenco, A.: Evaluation of modeled O3 using Measurement of Ozone by Airbus In-Service Aircraft (MOZAIC) data, J. Geophys. Res., 103, 25721–25737, 1998
   
  +
===Stratospheric Look-Up Table Photolysis===
−
'''Note:''' On HECToR this directory must be located on '''/work''' as the ''/home'' directory cannot be read at run time. This includes any symbolic links from /home to /work.
  
   
  +
In a chemistry scheme which has stratospheric chemistry, such as ''CheS/Strat'' and ''CheST/StratTrop'', if interactive photolysis is not used, then above 300hPa the look-up table approach of Lary and Pyle (1991)[1] is used (below 300hPa the tropospheric scheme described above is used). To use this scheme, in the UMUI go to '''Model Selection &rarr; Atmosphere &rarr; Model Configuration &rarr; UKCA Chemistry and Aerosols &rarr; PHOTO''' and click '''2D Photolysis Scheme'''. The code for this scheme is held in '''ukca_photolib.F90'''.
−
==Solution to Task 6.1: Create a new emissions file and use it in your job==
  
   
  +
'''References'''
−
Please see [[Solution to UKCA Chemistry and Aerosol Tutorial 5 Task 5.1 |this page]] for a solution to [[#Task 6.1: Create a new emissions file and use it in your job|Task 6.1]]
  
  +
# Lary, D. and Pyle, J.: Diffuse-radiation, twilight, and photochemistry, J. Atmos. Chem., 13, 393–406, 1991.
   
  +
===Interactive Photolysis===
−
==Task 6.2: make the required code changes to add your emission into UKCA==
  
   
  +
The Fast-J scheme (Wild ''et al'', 2000)[1] uses 7 different wavelength bins appropriate for the troposphere, and the Fast-JX scheme (Neu et al, 2007)[2] adds up to an extra 11 bins allowing use in the stratosphere.
−
<span style="color:green">'''TASK 6.2:''' You should now make the UKCA code changes to add your emission into the '''ALICE''' tracer. No run-time processing of this surface emission is required. You will also need to add-in the molar-mass of '''ALICE'''.</span>
  
   
  +
To use these schemes, in the UMUI go to '''Model Selection &rarr; Atmosphere &rarr; Model Configuration &rarr; UKCA Chemistry and Aerosols &rarr; PHOTO''' and click either '''FASTJ Photolysis Scheme''' or '''FASTJX Photolysis Scheme'''. You will then need to give the location of several input data files used by these schemes. The code for Fast-J is in the <tt>UKCA/</tt> directory in the '''fastj_*.F90''' files (controlled by '''ukca_fastj.F90'''), and the code for Fast-JX is in the '''fastjx_*.F90''' files (controlled by '''ukca_fastjx.F90''').
−
{| class="collapsible collapsed wikitable"
  
−
|-
  
−
! Hint
 
−
|-
  
−
| You can calculate the molar mass from the mass of air and the conversion factor defined in [[UKCA Chemistry and Aerosol Tutorial 5#Task 5.2: Add these two new tracers to UKCA|Task 5.2]].
  
−
|}
  
   
  +
Further details on the Fast-JX scheme, and how it is used in UKCA, can be found in [http://www.geosci-model-dev.net/6/161/2013/gmd-6-161-2013.html Telford ''et al'' (2013)][3].
−
'''Note:''' If you were unable to successfully complete Task 6.1 above, then please take a copy of the '''f''' job from the Tutorial experiment (''Tutorial: solution to Task 6.1 - adding a new chemical emissions ancillary file'') and work from there, as this will allow you to make only the changes required for this task. Please also make a new branch and merge-in branch '''branch''' at revision number '''rev''' to allow you to proceed.
  
   
  +
The Fast-J/Fast-JX data files are held in
−
You can also find a copy of an emissions ancillary file, with the required emissions, at
  
   
−
/work/n02/n02/ukca/Tutorial/vn8.4/Task6.1/solution/Task6.1_AR5_aero_2000.anc
 +
/work/n02/n02/hum/vn8.2/ctldata/UKCA/fastj
   
−
on ARCHER, and at
 +
on HECToR, and
   
−
/projects/ukca/Tutorial/vn8.4/Task6.1/solution/Task6.1_AR5_aero_2000.anc
 +
/projects/um1/vn8.2/ctldata/UKCA/fastj
   
 
on MONSooN.
  
on MONSooN.
   
  +
'''References'''
−
==Part 2. UKCA Code Changes==
  
  +
# Wild, O., Zhu, X., and Prather, M.: Fast-J: accurate simulation of in- and below-cloud photolysis in tropospheric chemical models, J. Atmos. Chem., 37, 245–282, doi:10.1023/A:1006415919030, 2000
−  
  +
# Neu, J., Prather, M., and Penner, J.: Global atmospheric chemistry: integrating over fractional cloud cover, J. Geophys. Res., 112, D11306, 12 pp., doi:10.1029/2006JD008007, 2007
−
===Changes to ukca_setd1defs.F90===
  
  +
# Telford, P. J., Abraham, N. L., Archibald, A. T., Braesicke, P., Dalvi, M., Morgenstern, O., O'Connor, F. M., Richards, N. A. D., and Pyle, J. A.: Implementation of the Fast-JX Photolysis scheme (v6.4) into the UKCA component of the MetUM chemistry-climate model (v7.3), Geosci. Model Dev., 6, 161-177, doi:10.5194/gmd-6-161-2013, 2013.
−  
−
The ukca_setd1defs.F90 tells UKCA what fields it should expect to find from the UM to allow it to run. Previously you edited this routine so that UKCA knew about your new tracers, now you must edited to tell it about your new emissions. This is done in two places
  
−  
−
====Add the species to em_chem_spec====
  
−  
−
You will need to find the '''em_chem_spec''' definition for the scheme that you are using. For example, the ''CheST/StratTrop'' chemistry is located in the ''Stratospheric Chemistry'' section and contained in the IF block controlled by <code>(L_ukca_strattrop .AND. .NOT. L_ukca_achem)</code> if you are not using aerosol chemistry, and by <code>(L_ukca_strattrop .AND. L_ukca_achem)</code> if you are using aerosol chemistry. For the former case this is defined by
  
−  
−
n_chem_emissions = 9
  
−
n_3d_emissions = 1 ! aircraft NOX
  
−
 
−
ALLOCATE(em_chem_spec(n_chem_emissions+n_3d_emissions))
  
−
em_chem_spec = &
  
−
(/'NO ','CH4 ','CO ','HCHO ', &
  
−
'C2H6 ','C3H8 ','Me2CO ','MeCHO ', &
  
−
'C5H8 ','NO_aircrft'/)
  
−  
−
You should edit the equivalent block. You should first increase the value of '''n_chem_emissions''' (for surface emissions) or '''n_3d_emissions''' (for 3D emissions) by the number of emissions that you are adding, and you should then add the names of the species (as they appear in the '''nm_spec''' array in ''ukca_setd1defs.F90'', and how they appear in the ''CHCH_DEFS'' specification in the ''ukca_chem_<span style="color:blue">scheme</span>.F90'' routine) that you are emitting into. Convention is that these are in ascending order by STASH code. The first 9 species in the above list are STASH codes 301-309, and <code>'NO_aircrft'</code> is STASH code 340, so all new species should be placed before <code>'NO_aircrft'</code>. You should make sure that the size of ''em_chem_spec'' is correct for the number of species within it.
  
−  
−
====Tell UKCA the STASH code associated with your new emission====
  
−  
−
In the [[UKCA & UMUI Tutorial 5#Code in ukca_setd1defs.F90 | previous ukca_setd1defs.F90 discussion]] above we saw the following block of code
  
−  
−
J = n_use_tracers
  
−
IF (n_chem_emissions+n_3d_emissions+n_mode_emissions > 0) THEN
  
−
DO i=1,n_chem_emissions + n_3d_emissions
  
−
UkcaD1Codes(J+i)%section = 0
  
−
UkcaD1Codes(J+i)%item = n_emiss_first+i-1 ! trop chemistry
  
−
UkcaD1Codes(J+i)%len_dim1 = row_length ! uses stash codes
  
−
UkcaD1Codes(J+i)%len_dim2 = rows ! 301-309 for
  
−
UkcaD1Codes(J+i)%required = .true. ! surface emissions
  
−
UkcaD1Codes(J+i)%prognostic = .true. ! from Section 0
  
−
! Special cases, emissions already available in UM
  
−
IF (em_chem_spec(i)(1:7) == 'SO2_low') THEN
  
−
UkcaD1Codes(J+i)%item = 58
  
−
IF (.NOT. L_SO2_SURFEM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='SO2 surface emissions from UM are not flagged'
  
−
errcode=58
  
−
 
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:7) == 'SO2_nat') THEN
  
−
UkcaD1Codes(J+i)%item = 121
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
IF (.NOT. L_SO2_NATEM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='SO2 natural emissions from UM are not flagged'
  
−
errcode=121
  
−
 
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:8) == 'SO2_high') THEN
  
−
UkcaD1Codes(J+i)%item = 126
  
−
IF (.NOT. L_SO2_HILEM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='SO2 high-level emissions are not flagged'
  
−
errcode = UkcaD1Codes(J+i)%item
  
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:3) == 'NH3') THEN
  
−
UkcaD1Codes(J+i)%item = 127
  
−
IF (.NOT. L_NH3_EM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='NH3 surface emissions from UM are not flagged'
  
−
errcode = UkcaD1Codes(J+i)%item
  
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i) == 'BC_fossil ') THEN
  
−
UkcaD1Codes(J+i)%item = 310
  
−
ELSEIF (em_chem_spec(i) == 'BC_biofuel') THEN
  
−
UkcaD1Codes(J+i)%item = 311
  
−
ELSEIF (em_chem_spec(i) == 'OC_fossil ') THEN
  
−
UkcaD1Codes(J+i)%item = 312
  
−
ELSEIF (em_chem_spec(i) == 'OC_biofuel') THEN
  
−
UkcaD1Codes(J+i)%item = 313
  
−
ELSEIF (em_chem_spec(i) == 'Monoterp ') THEN
  
−
UkcaD1Codes(J+i)%item = 314
  
−
ELSEIF (em_chem_spec(i) == 'NVOC ') THEN
  
−
UkcaD1Codes(J+i)%item = 315
  
−
ELSEIF (em_chem_spec(i) == 'BC_biomass') THEN
  
−
UkcaD1Codes(J+i)%item = 322
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ELSEIF (em_chem_spec(i) == 'OC_biomass') THEN
  
−
UkcaD1Codes(J+i)%item = 323
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ELSEIF (em_chem_spec(i) == 'SO2_biomas') THEN
  
−
UkcaD1Codes(J+i)%item = 324
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ELSEIF (em_chem_spec(i)(1:3) == 'DMS') THEN
  
−
UkcaD1Codes(J+i)%section = 17
  
−
UkcaD1Codes(J+i)%item = 205
  
−
UkcaD1Codes(J+i)%prognostic = .false.
  
−
IF (.NOT. L_DMS_EM .AND. (L_ukca_aerchem .OR. &
  
−
L_ukca_achem)) THEN
  
−
cmessage='DMS surface emissions from UM are not flagged'
  
−
errcode = UkcaD1Codes(J+i)%section*1000 + &
  
−
UkcaD1Codes(J+i)%item
  
−
CALL EREPORT('UKCA_SETD1DEFS',errcode,cmessage)
  
−
ENDIF
  
−
ELSEIF (em_chem_spec(i)(1:7) == 'NO_airc') THEN
  
−
UkcaD1Codes(J+i)%item = 340
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−
ENDIF
  
−
ENDDO
  
−
ENDIF
  
−  
−
You will need to add code at the end of this block to tell UKCA what STASH code is associated with which species. The best thing to do is to copy one of the explict blocks at the end, and adjust accordingly, e.g.
  
−  
−
For surface (2D) emissions you should add in the following
  
−  
−
ELSEIF (em_chem_spec(i) == '<span style="color:blue">Your Species</span>') THEN
  
−
UkcaD1Codes(J+i)%item = <span style="color:blue">Your Species STASH code (301-320)</span>
  
−  
−
and for 3D emissions you should add in the following
  
−  
−
ELSEIF (em_chem_spec(i) == '<span style="color:blue">Your Species</span>') THEN
  
−
UkcaD1Codes(J+i)%item = <span style="color:blue">Your Species STASH code (321-340)</span>
  
−
UkcaD1Codes(J+i)%len_dim3 = tr_levels
  
−  
−
Remember that the character string length for a UKCA species is '''10 characters'''.
  
−  
−
===Changes to ukca_constants.F90===
  
−  
−
As you are adding in an emission you may need to define the molar mass (in g/mol) of the species that you are emitting in to, if it is not already defined. Add this definition into '''ukca_constants.F90''', e.g.
  
−  
−
...
  
−
REAL, PARAMETER :: m_co = 28.
  
−
REAL, PARAMETER :: m_hcho = 30.
 
−
...
  
−  
−
===Changes to ukca_emission_ctl.F90===
  
   
  +
===Photolysis Reaction Definition===
−
You will need to ensure that the '''molmass''' array is filled for the tracer that you are emitting in to. This is filled in the '''WHERE''' block which has entries like
  
   
  +
The photolysis reactions are defined in the '''ukca_chem_<span style="color:blue">scheme</span>.F90''' routines using the '''ratj_t''' Fortran type specification, usually in several arrays.
−
...
  
  +
To format of this '''ratj_t''' type is
−
WHERE (em_chem_spec == 'CO ') molmass = m_co
  
−
WHERE (em_chem_spec == 'HCHO ') molmass = m_hcho
  
−
...
  
   
  +
ratj_t('Reactant 1','Reactant 2','<span style="color:blue">Product 1 </span>','<span style="color:red">Product 2 </span>','<span style="color:green">Product 3 </span>',&
−
You should add an additional line for each new emission.
 
  +
'<span style="color:purple">Product 4 </span>', <span style="color:blue">Fraction of Product 1 produced</span>, <span style="color:red">Fraction of Product 2 produced</span>, &
  +
<span style="color:green">Fraction of Product 3 produced</span>, <span style="color:purple">Fraction of Product 4 produced</span>, Quantum Yield, Look-up Label), &
   
  +
The '''Look-Up Label''' is used to define the file used for the 2D photolysis, and is used by Fast-J/Fast-JX to find the correct values for each species in the input data files. This is a 10-character string, although only the first '''7''' characters are read by Fast-JX.
−
For surface emissions, unless your new emissions data requires further run-time processing, such as adding a diurnal cycle (see how isoprene (C5H8) is treated below), then you will not need to make too many additional changes to the '''ukca_emission_ctl.F90''' routine. If you do need to add functionality such as this then you will need edit the following IF block:
  
   
  +
Examples of this type are
−
DO l=1,n_emissions
  
−
IF (advt(k) == em_chem_spec(l) .AND. &
  
−
em_chem_spec(l) == 'NO ' ) THEN
  
−
! Convert from kg NO2/m2/s to kg NO/m2/s
  
−
em_field(:,:,k) = emissions(:,:,l)*m_no/m_no2
  
−
ELSE IF (advt(k) == em_chem_spec(l)(1:3) .AND. &
  
−
em_chem_spec(l) == 'SO2_low ' ) THEN
  
−
! Convert from kg S/m2/s to kg SO2/m2/s and take off sulphate fraction
  
−
em_field(:,:,k) = emissions(:,:,l)* &
  
−
(1.0 - mode_parfrac/100.0)*m_so2/m_s
  
−
ELSE IF (advt(k) == em_chem_spec(l) .AND. &
  
−
em_chem_spec(l) == 'DMS ' ) THEN
  
−
! Convert from kg S/m2/s to kg DMS/m2/s
  
−
em_field(:,:,k) = emissions(:,:,l)*m_dms/m_s
  
−
ELSE IF (advt(k) == 'MeOH ' .AND. &
  
−
em_chem_spec(l) == 'NVOC ' ) THEN
  
−
! Convert from kg C/m2/s to kg CH3OH/m2/s
  
−
em_field(:,:,k) = emissions(:,:,l)*meoh_factor* &
  
−
m_meoh/(m_c*3.0)
  
−
ELSE IF (advt(k) == em_chem_spec(l) .AND. &
  
−
em_chem_spec(l) == 'Monoterp ' ) THEN
  
−
! Convert from kg C/m2/s to kg C10H16/m2/s
  
−
em_field(:,:,k) = emissions(:,:,l)*m_monoterp/(m_c*10.0)
  
−
! === biogenic emissions ===
  
−
ELSE IF (advt(k) == em_chem_spec(l) .AND. &
  
−
em_chem_spec(l) == 'C5H8 ') THEN
  
−
IF (L_ukca_diurnal_isopems) THEN
  
−
tmp_in_em_field(:,:) = emissions(:,:,l)*(m_isop/(5.0*m_c))
  
−
! DEPENDS ON: ukca_diurnal_isop_ems
  
−
! testdcycl = .TRUE.
  
−
CALL UKCA_DIURNAL_ISOP_EMS(row_length, rows, &
  
−
tmp_in_em_field, cos_zenith_angle, &
  
−
int_zenith_angle, &
  
−
sin_theta_latitude, FV_cos_theta_latitude, &
  
−
tan_theta_latitude, timestep, tmp_out_em_field,&
  
−
testdcycl)
  
−
em_field(:,:,k) = tmp_out_em_field(:,:)
  
−
ELSE
  
−
em_field(:,:,k) = emissions(:,:,l)*(m_isop/(5.0*m_c))
  
−
END IF
  
−
ELSE IF (advt(k) == em_chem_spec(l) ) THEN
  
−
em_field(:,:,k) = emissions(:,:,l)
  
−
ENDIF ! end advt(k)
  
   
  +
ratj_t('H2O2 ','PHOTON ','OH ','OH ',' ', &
−
If nothing needs to be done to the emissions field, then the final section
  
  +
' ', 0.0, 0.0, 0.0, 0.0, 100.000,'jh2o2 ') , &
  +
ratj_t('HCHO ','PHOTON ','HO2 ','HO2 ','CO ', &
  +
' ', 0.0, 0.0, 0.0, 0.0, 100.000,'jhchoa ') , &
   
  +
===Increase the size of JPPJ (and JPNR)===
−
ELSE IF (advt(k) == em_chem_spec(l) ) THEN
  
−
em_field(:,:,k) = emissions(:,:,l)
  
   
  +
As with the bimolecular, termolecular, and heterogeneous reactions, you will also need to increase the values of two parameters that UKCA needs. These are
−
just adds the emissions field from the ancillary file into the correct place to be emitted to the tracer.
 
   
  +
* '''JPPJ''' is the number of photolysis reactions
−
After this is done the UM boundary layer mixing scheme is called (''TR_MIX'') for each tracer, and this is also where the emissions are added to the tracer field.
  
  +
* '''JPNR''' is the total number of reactions
   
  +
These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the <code>&RUN_UKCA</code> namelist in '''CNTLATM'''. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the ''CNTLATM'' file in your <tt>$HOME/umui_jobs/<span style="color:blue">jobid</span></tt> directory.
−
For 3D emissions you will need to explicitly add this into the required tracer using the ''TRSRCE'' subroutine. This is done for lightning emissions (calculated on-line) and aircraft emissions (read-in from the multi-level user ancillary file). You should copy what is done for e.g. aircraft emissions and adapt accordingly.
  
   
−
==Solution to Task 6.2: make the required code changes to add your emission into UKCA==
 +
==Solution to Task 6.1: Add a bimolecular reaction==
   
−
Please see [[Solution to UKCA & UMUI Tutorial 5 Task 5.2 |this page]] for a solution to [[#Task 6.2: make the required code changes to add your emission into UKCA|Task 6.2]]
 +
Please see [[Solution to UKCA Chemistry and Aerosol Tutorial 6 Task 6.1 |this page]] for a solution to [[#Task 6.1: Add a bimolecular reaction|Task 6.1]].
   
 
----
  
----

Revision as of 11:15, 25 February 2014

Back to UKCA Chemistry and Aerosol Tutorials

What you will learn in this tutorial

During this tutorial you will learn how UKCA specifies different chemical reactions. You will then add a new reaction involving the new tracers that you have added.

Task 6.1: Add a bimolecular reaction

TASK 6.1: You should now add in the bimolecular reaction of ALICE with OH to form BOB and a secondary organic compound (labelled in UKCA as SEC_ORG). This reaction is given by:f

Parameter Value
2.70E-11
0.00
-390.00

Note: If you were unable to successfully complete Task 5.2, then please take a copy of the f job from the Tutorial experiment (Tutorial: solution to Task 5.2 - adding new chemical emissions in UKCA) and work from there, as this will allow you to only make the required changes. Please also make a new branch and merge-in branch branch at revision number rev to allow you to proceed.

Adding new Chemical Reactions

UKCA currently uses two different methods of defining the chemical reactions solved in the model. The first is a backward Euler solver, and is used for the RAQ and StdTrop chemistry schemes where the solver itself is created by a code-writer. The second makes use of the ASAD chemical integration software package, and is used for the CheT/TropIsop, CheS/Strat, and CheST/StratTrop chemistry schemes. ASAD can use many different solvers, although currently it uses symbolic Newton-Raphson solver. In this tutorial we will only consider the ASAD framework, as this is easily extended by a user.

ASAD considers four different types of chemical reactions: bimolecular reactions, termolecular reactions, heterogeneous reactions, and photolysis reactions. To make changes and add reactions you will need to make changes to the UKCA source code which can be found in 

vn8.4_your_branch_name/src/atmosphere/UKCA

During this tutorial you will be tasked with adding a new reaction into your branch.

Biomolecular Reactions

For most bimolecular reactions, it is sufficient to provide the , , and coefficients that are used to compute the rate coefficient from the Arrhenius expression

Bimolecular Reaction Definition

The bimolecular reactions are defined in the ukca_chem_scheme.F90 routines using the ratb_t Fortran type specification, and are held in arrays. At the end of this routine the ratb_defs_scheme array is created from these, and if that scheme is selected in UKCA these reactions are copied across into the master ratb_defs array.

The format of this ratb_t type is 

ratb_t('Reactant 1','Reactant 2','Product 1 ','Product 2 ','Product 3 ',&
'Product 4 ',  ,  ,  , Fraction of Product 1 produced, Fraction of Product 2 produced, &
Fraction of Product 3 produced, Fraction of Product 4 produced), &

If fractional products are not required for a reaction, then the fraction of each product formed should be set to 0.000. If fractional products are required for any one of the products then the fraction of each product formed should be set to its correct value.

The specifications of the individual reactions are done as, e.g. 

ratb_t('O3        ','C5H8      ','HO2       ','OH        ','          ',& ! B133 
'          ',  3.33E-15,  0.00,   1995.00, 0.750, 0.750, 0.000, 0.000), & ! B133 IUPAC2007*   
...
ratb_t('OH        ','C5H8      ','ISO2      ','          ','          ',& ! B144 
'          ',  2.70E-11,  0.00,   -390.00, 0.000, 0.000, 0.000, 0.000), & ! B144 IUPAC2009   
...
ratb_t('OH        ','HCl       ','H2O       ','Cl        ','          ',& ! B159 
'          ',  1.80E-12,  0.00,    250.00, 0.000, 0.000, 0.000, 0.000), & ! B159 JPL2011

The first reaction in these examples takes its kinetic data from IUPAC. Going to this website, this reaction is defined here. The second reaction above takes its kinetic data from NASA's Jet Propulsion Laboratory. The rate for this can be found on page 1-19 of the JPL2011 document. When adding new reactions you will need to increment the size of the array holding the ratb_t type.

To add new bimolecular reactions you will need to append equivalent lines for the new reactions to the end of the ratb_defs_scheme array (increasing the array sizes accordingly). If there is a reaction that is an exception to the general Arrhenius equation then special code needs to be placed in the asad_bimol.F90 routine, which is held in the UKCA/ source-code directory.

Increase the size of JPBK (and JPNR)

As well as adding these reactions to the ukca_chem_scheme.F90 routine (and incrementing the size of the arrays in that routine accordingly, you will also need to increase the values of two parameters that UKCA needs. These are

  • JPBK is the number of bimolecular reactions
  • JPNR is the total number of reactions

These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the &RUN_UKCA namelist in CNTLATM. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the CNTLATM file in your $HOME/umui_jobs/jobid directory.

Termolecular Reactions

As well as defining reactions involving a third body, the termolecular rate definition can also be used to define unimolecular reactions.

The pressure and temperature dependent rate, , of a termolecular reaction is given by

where the low pressure rate constant is given by

and the high pressure rate constant is given by

Termolecular Reaction Definition

The termolecular reactions are defined in the ukca_chem_scheme.F90 routines using the ratt_t Fortran type specification, and are usually held in one single array (there are not usually enough reactions to require splitting the reactions over several arrays).

To format of this ratt_t type is 

ratt_t('Reactant 1','Reactant 2','Product 1 ','Product 2 ', , &
,  ,  , ,  ,  , Fraction of Product 1 produced, Fraction of Product 2 produced), &

and as in ratb_t, where the fraction of a product should be set to 0.000 if this functionality does not need to be used.

The value is used to define the value by

If then
else

as may or may not be highly temperature dependent.

Examples of these reactions are 

ratt_t('N2O5      ','m         ','NO2       ','NO3       ',     0.3,    & ! T023  
  1.30E-03, -3.50, 11000.00,  9.70E+14,  0.10, 11080.00, 0.000, 0.000), & ! T023 IUPAC 2002   
ratt_t('NO        ','NO        ','NO2       ','NO2       ',     0.0,    & ! T024  
  3.30E-39,  0.00,  -530.00,  0.00E+00,  0.00,     0.00, 0.000, 0.000)  & ! T024 IUPAC 2001

To add new termolecular reactions you will need to append equivalent lines for the new reactions to the end of the ratt_defs_scheme array (increasing the array sizes accordingly).

Increase the size of JPTK (and JPNR)

As with the bimolecular reactions, you will also need to increase the values of two parameters that UKCA needs. These are

  • JPTK is the number of termolecular reactions
  • JPNR is the total number of reactions

These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the &RUN_UKCA namelist in CNTLATM. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the CNTLATM file in your $HOME/umui_jobs/jobid directory.

Heterogeneous Reactions

Heterogeneous reactions are those that occur on aerosol surfaces. There is no functional form defined for these reactions, with special code needed to be added for each case.

Heterogeneous Reaction Definition

The heterogeneous reactions are defined in the ukca_chem_scheme.F90 routines using the rath_t Fortran type specification, usually in one array. To format of this rath_t type is 

rath_t('Reactant 1','Reactant 2','Product 1 ','Product 2 ','Product 3 ',&
'Product 4 ', Fraction of Product 1 produced, Fraction of Product 2 produced, &
Fraction of Product 3 produced, Fraction of Product 4 produced), &

i.e. there is no rate information provided. For reactions on PSCs special code has been added to the routines in ukca_hetero_mod.F90, and for other reactions there is code in asad_hetero.F90. Examples of this type are 

rath_t('ClONO2    ','H2O       ','HOCl      ','HONO2     ','          ', &
'          ', 0.000, 0.000, 0.000, 0.000), &
...
rath_t('SO2       ','H2O2      ','NULL0     ','          ','          ', & !HSO3+H2O2(aq)
'          ', 0.000, 0.000, 0.000, 0.000),                               &

To add new heterogeneous reactions you will need to append equivalent lines for the new reactions to the end of the ratt_defs_scheme array (increasing the array sizes accordingly), before adding code to either ukca_hetero_mod.F90 or asad_hetero.F90.

Increase the size of JPHK (and JPNR)

As with the bimolecular and termolecular reactions, you will also need to increase the values of two parameters that UKCA needs. These are

  • JPHK is the number of heterogeneous reactions
  • JPNR is the total number of reactions

These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the &RUN_UKCA namelist in CNTLATM. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the CNTLATM file in your $HOME/umui_jobs/jobid directory.

Photolysis Reactions

These define a reaction where a chemical compound is broken down by photons. There is no functional form defined for this type of reaction. Instead, either (in the troposphere) input files are used to define the reaction rates for each species, while (in the stratosphere) on-line look-up tables are generated for the rates for each species, or separate photolysis codes, Fast-J or Fast-JX, are used to interactively calculate the rate of reaction throughout the troposphere (for Fast-J) or the whole atmosphere (for Fast-JX). These interactive schemes are preferred as they take the effect of aerosols or clouds into account at each timestep, allowing for more feedbacks to be investigated. In the upper stratosphere there are some wavelength regions that Fast-JX does not consider, and so the 3D on-line look-up tables are also used for these regions.

Tropospheric Off-Line Photolysis

If Fast-JX is not being used, then the off-line two-dimensional (zonally average) tropospheric photolysis is used (for all schemes). It is based on the work of Hough (1988)[1] and Law et al (1998)[2].

This scheme makes use of datafiles which define the reaction rate for a particular species (e.g. H2O2), or if no rate is known, a nil rate can be used. For UM 8.2 these files (in ASCII format) can be found in 

/work/n02/n02/hum/vn8.2/ctldata/UKCA/tropdata/photol

on HECToR, and in 

/projects/um1/vn8.2/ctldata/UKCA/tropdata/photol

on MONSooN. To use this scheme, in the UMUI go to Model Selection → Atmosphere → Model Configuration → UKCA Chemistry and Aerosols → PHOTO and click 2D Photolysis Scheme. You will then need to give the location of the files (above). The code controlling this scheme is held in ukca_phot2d.F90.

It is advised that this scheme is no longer used, and interactive photolysis should be used instead. For the CheS/Strat or CheST/StratTrop schemes, Fast-JX should be used as this covers the stratosphere as well as the troposphere.

References

  1. Hough, A. M.: The calculation of photolysis rates for use in global modelling studies, Tech. rep., UK Atomic Energy Authority, Harwell, Oxon., UK, 1988
  2. Law, K., Plantevin, P., Shallcross, D., Rogers, H., Pyle, J., Grouhel, C., Thouret, V., and Marenco, A.: Evaluation of modeled O3 using Measurement of Ozone by Airbus In-Service Aircraft (MOZAIC) data, J. Geophys. Res., 103, 25721–25737, 1998

Stratospheric Look-Up Table Photolysis

In a chemistry scheme which has stratospheric chemistry, such as CheS/Strat and CheST/StratTrop, if interactive photolysis is not used, then above 300hPa the look-up table approach of Lary and Pyle (1991)[1] is used (below 300hPa the tropospheric scheme described above is used). To use this scheme, in the UMUI go to Model Selection → Atmosphere → Model Configuration → UKCA Chemistry and Aerosols → PHOTO and click 2D Photolysis Scheme. The code for this scheme is held in ukca_photolib.F90.

References

  1. Lary, D. and Pyle, J.: Diffuse-radiation, twilight, and photochemistry, J. Atmos. Chem., 13, 393–406, 1991.

Interactive Photolysis

The Fast-J scheme (Wild et al, 2000)[1] uses 7 different wavelength bins appropriate for the troposphere, and the Fast-JX scheme (Neu et al, 2007)[2] adds up to an extra 11 bins allowing use in the stratosphere.

To use these schemes, in the UMUI go to Model Selection → Atmosphere → Model Configuration → UKCA Chemistry and Aerosols → PHOTO and click either FASTJ Photolysis Scheme or FASTJX Photolysis Scheme. You will then need to give the location of several input data files used by these schemes. The code for Fast-J is in the UKCA/ directory in the fastj_*.F90 files (controlled by ukca_fastj.F90), and the code for Fast-JX is in the fastjx_*.F90 files (controlled by ukca_fastjx.F90).

Further details on the Fast-JX scheme, and how it is used in UKCA, can be found in Telford et al (2013)[3].

The Fast-J/Fast-JX data files are held in 

/work/n02/n02/hum/vn8.2/ctldata/UKCA/fastj

on HECToR, and 

/projects/um1/vn8.2/ctldata/UKCA/fastj

on MONSooN.

References

  1. Wild, O., Zhu, X., and Prather, M.: Fast-J: accurate simulation of in- and below-cloud photolysis in tropospheric chemical models, J. Atmos. Chem., 37, 245–282, doi:10.1023/A:1006415919030, 2000
  2. Neu, J., Prather, M., and Penner, J.: Global atmospheric chemistry: integrating over fractional cloud cover, J. Geophys. Res., 112, D11306, 12 pp., doi:10.1029/2006JD008007, 2007
  3. Telford, P. J., Abraham, N. L., Archibald, A. T., Braesicke, P., Dalvi, M., Morgenstern, O., O'Connor, F. M., Richards, N. A. D., and Pyle, J. A.: Implementation of the Fast-JX Photolysis scheme (v6.4) into the UKCA component of the MetUM chemistry-climate model (v7.3), Geosci. Model Dev., 6, 161-177, doi:10.5194/gmd-6-161-2013, 2013.

Photolysis Reaction Definition

The photolysis reactions are defined in the ukca_chem_scheme.F90 routines using the ratj_t Fortran type specification, usually in several arrays. To format of this ratj_t type is 

ratj_t('Reactant 1','Reactant 2','Product 1 ','Product 2 ','Product 3 ',&
'Product 4 ', Fraction of Product 1 produced, Fraction of Product 2 produced, &
Fraction of Product 3 produced, Fraction of Product 4 produced, Quantum Yield, Look-up Label), &

The Look-Up Label is used to define the file used for the 2D photolysis, and is used by Fast-J/Fast-JX to find the correct values for each species in the input data files. This is a 10-character string, although only the first 7 characters are read by Fast-JX.

Examples of this type are 

ratj_t('H2O2      ','PHOTON    ','OH        ','OH        ','          ', &
     '          ',    0.0,   0.0,   0.0,   0.0, 100.000,'jh2o2     ') ,  &
ratj_t('HCHO      ','PHOTON    ','HO2       ','HO2       ','CO        ', &
     '          ',    0.0,   0.0,   0.0,   0.0, 100.000,'jhchoa    ') ,  &

Increase the size of JPPJ (and JPNR)

As with the bimolecular, termolecular, and heterogeneous reactions, you will also need to increase the values of two parameters that UKCA needs. These are

  • JPPJ is the number of photolysis reactions
  • JPNR is the total number of reactions

These are set automatically in the UMUI (depending on what scheme is chosen), and are placed in the &RUN_UKCA namelist in CNTLATM. You will need to make a hand-edit to change these accordingly. The current values can be found by saving and processing the job, and then viewing the CNTLATM file in your $HOME/umui_jobs/jobid directory.

Solution to Task 6.1: Add a bimolecular reaction

Please see this page for a solution to Task 6.1.

Written by Luke Abraham 2014

Retrieved from "http://www.ukca.ac.uk/wiki/index.php?title=UKCA_Chemistry_and_Aerosol_Tutorial_6&oldid=2258"
Powered by MediaWiki