Photochemistry
Contents
Generic PCM
General Use
Compilation
Photochemistry in the generic model is largely independent of the and does not require to recompile it (just make sure that you have linked the LAPACK and BLAS libraries).
callphys.def
In order to use photochemistry, you simply need to have the following line in callphys.def:
photochem = .true.
traceur.def
All chemical species must be listed in traceurs.def (which needs to follow the modern layout).
Reaction network
Besides the species- and reaction-specific data files (see below), you need one file: reactfile (located in the chemnetwork folder) for the reaction network. The first line of the file is ignored, then for each line, the syntax of the reactfile is as follows:
- Columns 1-50 are reserved for reactants. Their names must match those given traceurs.def, in particular, if you have both water vapour and water ice as tracers, named h2o_vap and h2o_ice respectively, the name to use for water vapour in chemical reaction is h2o_vap. Reactant species must be separated by any number of blank spaces.
- Columns 51-100 are reserved for products, with the same rules applying. Notice that in the case of a third body reaction, the third body appears only with the reactants and is omitted on the products side.
- Column 101 is an integer designating the type of reaction: 0 for photolysis, 1, 2 or 3 for other reactions, depending on their reaction rate formula. The four next spaces are blank.
- Column 106 and after are dedicated to parameters specific to each reaction type.
Note that lines can be commented out by starting them with the character "!".
Photolysis
Photolysis reactions are identified by a 0 at column 101 of the reactfile. They need two types of files (which need to be placed in datadir/cross_sections/):
- cross sections (syntax: first line: wavelength in nm, second line cross section in cm2, first line is ignored)
- yield rates (syntax: first line: wavelength in nm, second line: photodissociation yield, first line is ignored)
These data can be found in benchmark cases for common species (e.g. water), or on the following pages for more exotic ones: Mainz spectral atlas, phidrates, Leiden University. The parameters for photolysis reaction in the reactfile are given as follows, from column 106: an identifier character string starting with "j", an integer corresponding to the number of cross section files (one for each temperature range), the lower bound of each temperature range, the name of each cross section file, the name of the yield rate file. For example, the following line describes the photolysis of O2 into two O atoms in the ground state:
o2 hv o o 0 jo2_o 4 150 200 250 300 O2_150.abs O2_200.abs O2_250.abs O2_300.abs o2_o.yld
In this case, 4 different cross section files exist, corresponding to the temperature range 150 K-200 K, 200 K-250 K, 250 K-300 K, >300 K.
Reaction rate formulae
Many reaction rate parametrizations can be found in here. The model currently accounts for three different formula for the reaction rate (more exotic ones can be hard-coded in the chimiedata_h.F90 file, as indicated in the TRAPPIST-1e example). The chosen formula is by 1, 2, or 3 at column 51 or the reactfile.They read respectively:
- The simple form of the reaction rate is
\[k=a\left(\frac{T}{T_0}\right)^ce^{\frac{-b}{T} }[M]^d\]Its parameters (to give in the reactfile in the same order, starting from column 106) are \(a, b, c, T_0, d\).
- Termolecular reactions have low- and high-pressure-limiting rate constant (\(k_0\) and \(k_\infty\), respectively). Their effective rate constant is then given by
\[k=ge^{-\frac{h}{T} }+\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} }[M]^{d_{\text up} } }{1+\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M]^{d_{\text down} } }fc^{\frac{1}{1+\left[\log_{10}\left(\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M] \right) \right]^2 } }\]Its parameters (to give in the reactfile in the same order, starting from column 106) are \(k_0, n, a_0, k_\infty, m, b_0, T_0, f_c, g, h, d_{\text up}, d_{\text down}\).
- Finally, a last form of the reaction rate is
\[\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} }[M]^{d_{\text up} } }{1+\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M]^{d_{\text down} } }[F_{\text cent}]^X,\]with \(F_{\text cent}=1-ae^{-\frac{T}{T^{***} } }+ae^{-\frac{T}{T^* } }+ae^{-\frac{T^{**}}{T } }\), and \[X=\left[1+\left[\frac{\log_{10}\left(\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M] \right)+c }{N-d\left(\log_{10}\left(\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M] \right) \right)+c } \right]^2 \right]^{-1},\]with \(c=-0.4-0.67\log_{10}(F_{\text cent})\), \(N=0.75-1.27\log_{10}(F_{\text cent})\), and \(d=0.14\).
- Its parameters (to give in the reactfile in the same order, starting from column 106) are \(k_0, n, a_0, k_\infty, m, b_0, T_0, a, T^{***}, T^*, T^{**}, d_{\text up}, d_{\text down}\).
For example, the following line describes the 3-body reaction CO + O + M = CO2 + M, whose reaction rate is \(k=2.2\times10^{-33}e^{-\frac{1780}{T} }[M]\) (be careful, here \([M]\) corresponds to the molecular density of the atmosphere, not that of the third body):
co o M co2 1 2.2e-33 -1780.0 0.0 1 1.0
Using hard-coded rate constants
You might sometimes want to track a reaction whose rate constant is not capture by either of the three above parametrizations. In this case you'll need to get your hands dirty and hard-code the rate constant parametrization in the code. Notice that the reaction should not be included in the reactfil in this case (it is however a good practice to indicate it as a commented line and write "hard-coded', both to keep track of your network and because the post-processing tools will notify it by reading the reactfile). Here is how to proceed:
First, identify the type of your reaction: currently, the model distinguishes between 1) bimolecular reactions (A + B -> C (+ D)), 2) quadratic reactions (A + A -> B (+C)), 3) quenching reactions (A + C -> B + C), and 4) heterogeneous reactions (e.g. A + ice -> B + C). Bimolecular reactions are given type v4, quadratic reactions are given type v3, and both quenching an heterogeneous reactions are given type vphot (the same as photodissociation reactions).
Then, open the file LMDZ.GENERIC/libf/aeronostd/chimiedata_h.F90. This file has a dedicated section for hard-coded reactions. Some of them are already there, commented out. We'll use one of them (the reaction CO + OH -> CO2 + H) as a example on how to proceed. This is a bimolecular reaction, so we first need to increment the counter of hard-coded reactions of type v_4:
integer, parameter :: n4_hard_coding = 1
In the master branch version of the code, the variable n4_hard_coding is set to 0. If you're version of the code already has hard-coded reactions, then simply add one to the current value of n4_hard_coding.
In the same file, you will find a subroutine called indice_HC for registering hard-coded reactions. You can create a space for your reaction, taking example of what's already there. For instance our example reaction is already there:
!===========================================================
! e001 : CO + OH -> CO2 + H
!===========================================================
! nb_reaction_4 = nb_reaction_4 + 1
! indice_4(nb_reaction_4) = z4spec(1.0, indexchim('co'), 1.0, indexchim('oh'), 1.0, indexchim('co2'), 1.0, indexchim('h'))
To use it, you must uncomment the two last lines. The top one increments the number of reactions of type v4, and the second one registers the reacting species for your new reaction.
Then go to the subroutine reactionrates_HC. This is where you will write down the parametrization formula into the reaction rates table (the variable v_4 for bimolecular reactions, v_3 for quadratic reactions, or v_phot for all others). Again, you can take example of what's already there. If you want to use our example reaction, find the dedicated space, opening with the line:
!--- e001: oh + co -> co2 + h
and simply uncomment the following lines:
!nb_reaction_4 = nb_reaction_4 + 1
to update the index at which the rate is written in table v_4, and then (for instance):
! jpl 2015 !do ilev = 1,nlayer !!oh + co -> h + co2 ! rate1 = 1.5e-13*(t(ilev)/300.)**(0.0) !!oh + co + m -> hoco ! ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0) ! ak1 = 1.1e-12*(t(ilev)/300.)**(1.3) ! rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1) ! xpo2 = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2) ! ! v_4(ilev,nb_reaction_4) = rate1 + rate2*0.6**xpo2 !end do
(except the first one) to use the parametrization from the JPL 2015 collection (linked above). Notice that several parametrizations are implemented; to use another one simply uncomment another block, or write your own.
