Photochemistry of the Generic PCM - Revision history

Milcareck at 10:38, 14 January 2026

2026-01-14T10:38:45Z

Maxime Maurice: /* Boundary conditions */

2025-03-05T09:57:08Z

‎Boundary conditions

Maxime Maurice: Created page with "= 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..."

2025-03-05T08:55:30Z

Created page with "= 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..."

New page

= 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.
Other flags related to photochemistry are:
* '''photoheat''': turn on or off heating by UV absorption.
* '''jonline''': turn on or off online photolysis.
* '''stellarflux''': name of the stellar flux file.
* '''depos''': turn on or off deposition.

== traceur.def ==
All chemical species must be listed in [[The_traceur.def_Input_File|'''traceurs.def''']] ('''which needs to follow the modern layout''').

== Reaction network ==
Besides the species- and reaction-specific data files ([[#Photolysis reaction|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 [[The_traceur.def_Input_File|'''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 ==
=== Online photolysis ===
In the baseline use of the photochemistry model, photolysis are computed at the runtime (online), by solving the UV radiative transfer using the TUV solver, considering species-specific UV cross section. o do so, first make sure that the apropriate flag is set to true in the '''callphys.def''':
jonline = .true.
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: [https://www.uv-vis-spectral-atlas-mainz.org/uvvis/ Mainz spectral atlas], [https://phidrates.space.swri.edu/ phidrates], [https://home.strw.leidenuniv.nl/~moldata/ 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.
=== Offline photolysis ===
=== Stellar flux ===
In order to compute photolysis (is it also the case for offline photolysis?), you need to provide a UV stellar spectrum. This file must be located in datadir/stellar_spectra/photochem_stellar_spectra, an is composed of a 2-lines header and 2 columns, the left one indicating wavelength (in nm) and the right one indicating the flux (in photons/cm²/s/nm).
=== Constant photodissociation rates ===
Sometimes photodissociation rates are simply given as a constant throughout the atmosphere. One way to reproduce that is to treat these reactions as type 1 reaction (see below).

== Reaction rate formulae ==
Many reaction rate parametrizations can be found in [https://jpldataeval.jpl.nasa.gov/pdf/NASA-JPL%20Evaluation%2019-5.pdf 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 [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Other_GCM_Configurations_worth_knowing_about#Compiling_a_test_case_.28TRAPPIST-1e.29 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 :<math>k=a\left(\frac{T}{T_0}\right)^ce^{\frac{-b}{T} }[M]^d</math>Its parameters (to give in the '''reactfile''' in the same order, starting from column 106) are <math>a, b, c, T_0, d</math>.
# Termolecular reactions have low- and high-pressure-limiting rate constant (<math>k_0</math> and <math>k_\infty</math>, respectively). Their effective rate constant is then given by :<math>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 } }</math>Its parameters (to give in the '''reactfile''' in the same order, starting from column 106) are <math>k_0, n, a_0, k_\infty, m, b_0, T_0, f_c, g, h, d_{\text up}, d_{\text down}</math>.
# Finally, a last form of the reaction rate is :<math>\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,</math>with <math>F_{\text cent}=1-ae^{-\frac{T}{T^{***} } }+ae^{-\frac{T}{T^* } }+ae^{-\frac{T^{**}}{T } }</math>, and :<math>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},</math>with <math>c=-0.4-0.67\log_{10}(F_{\text cent})</math>, <math>N=0.75-1.27\log_{10}(F_{\text cent})</math>, and <math>d=0.14</math>.
::Its parameters (to give in the '''reactfile''' in the same order, starting from column 106) are <math>k_0, n, a_0, k_\infty, m, b_0, T_0, a, T^{***}, T^*, T^{**}, d_{\text up}, d_{\text down}</math>.
For example, the following line describes the 3-body reaction CO + O + M = CO2 + M, whose reaction rate is <math>k=2.2\times10^{-33}e^{-\frac{1780}{T} }[M]</math> (be careful, here <math>[M]</math> 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

= Boundary conditions =
== Deposition velocity ==
== Production rate ==
== Escape ==

= 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 reactfile 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.

== Hard-coding reactions with 3 products ==
The model does not natively handle reactions with 3 products. The trick is to instead split it into 2 "semi-reactions", for instance instead of HNO3 + OH -> H2O + O + NO2, add 0.5 HNO3 + 0.5 OH -> H2O + 0.5 NO2 and 0.5 HNO3 + 0.5 OH -> O + 0.5 NO2 (or any stoichiometrically-equivalent couple). Each semi-reaction then has the same constant rate as the total reaction. Don't forget to count them as two distinct reactions.

= Python post-processing tools for photochemistry =
A mini python library is available in the '''util''' directory of the Generic model (LMDZ.GENERIC/utils), along with a demonstration jupyter notebook.

@@ Line 17: / Line 17: @@
 == Reaction network ==
-Besides the species- and reaction-specific data files ([[#Photolysis reaction|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:
+Besides the species- and reaction-specific data files ([[#Photolysis reaction|see below]]), you need one file: '''network.def''' (previously '''reactfile''' where it was 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 '''network.def''' is as follows:
 * Columns 1-50 are reserved for reactants. Their names must match those given [[The_traceur.def_Input_File|'''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.
@@ Line 28: / Line 28: @@
 In the baseline use of the photochemistry model, photolysis are computed at the runtime (online), by solving the UV radiative transfer using the TUV solver, considering species-specific UV cross section. o do so, first make sure that the apropriate flag is set to true in the '''callphys.def''':
   jonline = .true.
-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/'''):
+Photolysis reactions are identified by a 0 at column 101 of the '''network.def'''. 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 cm<sup>2</sup>, 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: [https://www.uv-vis-spectral-atlas-mainz.org/uvvis/ Mainz spectral atlas], [https://phidrates.space.swri.edu/ phidrates], [https://home.strw.leidenuniv.nl/~moldata/ 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 O<sub>2</sub> into two O atoms in the ground state:
+The parameters for photolysis reaction in the '''network.def''' 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 O<sub>2</sub> 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.
@@ Line 42: / Line 42: @@
 == Reaction rate formulae ==
-Many reaction rate parametrizations can be found in [https://jpldataeval.jpl.nasa.gov/pdf/NASA-JPL%20Evaluation%2019-5.pdf 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 [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Other_GCM_Configurations_worth_knowing_about#Compiling_a_test_case_.28TRAPPIST-1e.29 the TRAPPIST-1e example]). The chosen formula is by 1, 2, or 3 at column 51 or the '''reactfile'''.They read respectively:
+Many reaction rate parametrizations can be found in [https://jpldataeval.jpl.nasa.gov/pdf/NASA-JPL%20Evaluation%2019-5.pdf 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 [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Other_GCM_Configurations_worth_knowing_about#Compiling_a_test_case_.28TRAPPIST-1e.29 the TRAPPIST-1e example]). The chosen formula is by 1, 2, or 3 at column 51 or the '''network.def'''.They read respectively:
-# The simple form of the reaction rate is <p>:<math>k=a\left(\frac{T}{T_0}\right)^ce^{\frac{-b}{T} }[M]^d</math>Its parameters (to give in the '''reactfile''' in the same order, starting from column 106) are <math>a, b, c, T_0, d</math>.</p>
+# The simple form of the reaction rate is <p>:<math>k=a\left(\frac{T}{T_0}\right)^ce^{\frac{-b}{T} }[M]^d</math>Its parameters (to give in the '''network.def''' in the same order, starting from column 106) are <math>a, b, c, T_0, d</math>.</p>
-# Termolecular reactions have low- and high-pressure-limiting rate constant (<math>k_0</math> and <math>k_\infty</math>, respectively). Their effective rate constant is then given by <p>:<math>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 } }</math>Its parameters (to give in the '''reactfile''' in the same order, starting from column 106) are <math>k_0, n, a_0, k_\infty, m, b_0, T_0, f_c, g, h, d_{\text up}, d_{\text down}</math>.</p>
+# Termolecular reactions have low- and high-pressure-limiting rate constant (<math>k_0</math> and <math>k_\infty</math>, respectively). Their effective rate constant is then given by <p>:<math>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 } }</math>Its parameters (to give in the '''network.def''' in the same order, starting from column 106) are <math>k_0, n, a_0, k_\infty, m, b_0, T_0, f_c, g, h, d_{\text up}, d_{\text down}</math>.</p>
 # Finally, a last form of the reaction rate is <p>:<math>\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,</math>with <math>F_{\text cent}=1-ae^{-\frac{T}{T^{***} } }+ae^{-\frac{T}{T^* } }+ae^{-\frac{T^{**}}{T } }</math>, and :<math>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},</math>with <math>c=-0.4-0.67\log_{10}(F_{\text cent})</math>, <math>N=0.75-1.27\log_{10}(F_{\text cent})</math>, and <math>d=0.14</math>.
-::Its parameters (to give in the '''reactfile''' in the same order, starting from column 106) are <math>k_0, n, a_0, k_\infty, m, b_0, T_0, a, T^{***}, T^*, T^{**}, d_{\text up}, d_{\text down}</math>.</p>
+::Its parameters (to give in the '''network.def''' in the same order, starting from column 106) are <math>k_0, n, a_0, k_\infty, m, b_0, T_0, a, T^{***}, T^*, T^{**}, d_{\text up}, d_{\text down}</math>.</p>
 For example, the following line describes the 3-body reaction CO + O + M = CO<sub>2</sub> + M, whose reaction rate is <math>k=2.2\times10^{-33}e^{-\frac{1780}{T} }[M]</math> (be careful, here <math>[M]</math> 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
@@ Line 72: / Line 72: @@
 = 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 reactfile 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:
+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 network.def 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 network.def). 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).

@@ Line 51: / Line 51: @@
 = Boundary conditions =
-== Deposition velocity ==
+All processes described in this section are implemented in the '''deposition_source.F90''' file in LMDZ.GENERIC/libf/aeronostd.
-== Production rate ==
+For each chemical species, two types of surface boundary conditions can be enforced in the chemistry, specified by the '''SF_mode''' in the [[The_traceur.def_Input_File|traceur.def]] file. If none of the flags mentioned below is given, the default boundary conditions are no flux.
 == Escape ==
 = Using hard-coded rate constants =