Difference between revisions of "FlexChem"

From Geos-chem
Jump to: navigation, search
(Add production & loss families (if desired))
 
(13 intermediate revisions by the same user not shown)
Line 1: Line 1:
On this page we provide information about FlexChem in GEOS-Chem.
+
This content has been migrated to:
  
== Overview ==
+
* [https://kpp.readthedocs.io <tt>kpp.readthedocs.io</tt>]: User manual for The Kinetic Preprocessor (KPP)
 
+
* [https://geos-chem.readthedocs.io/en/latest/geos-chem-shared-docs/supplemental-guides/using-kpp-with-gc.html '''Update chemical mechanisms with KPP''' at <tt>geos-chem.readthedocs.io</tt>]
The '''[https://github.com/geoschem/kpp Kinetic Pre Processor (KPP) package]''' (A. Sandu et al) creates optimized chemical solver code in Fortran-90 from mechanism defined in user-editable text files. The resulting Fortran-90 code can be added into chemical models such as GEOS-Chem.  Adding changes to a mechanism can be simply done by editing the mechanism's configuration files and then re-running KPP to generate new Fortran output files.
+
 
+
FlexChem is '''[https://github.com/geoschem/kpp/tree/GC_updates a clean implementation of KPP'''] that has been customized for GEOS-Chem v11-01 and later versions.  Within FlexChem, there is a single supported chemical mechanism (named [https://github.com/geoschem/geos-chem/blob/main/KPP/fullchem/fullchem.eqn fullchem]), but users may also define their own custom mechanisms.
+
 
+
The source code in <tt>flexchem_mod.F90</tt> serves as the connection between the chemical mechanism solver files created by KPP and GEOS-Chem. It passes initial species concentrations, photolysis rates, meteorology fields, etc. to the chemical mechanism solver source code files created by KPP.  These source code files compute reaction rates, perform the forward integration, and pass the updated species concentrations back to GEOS-Chem.
+
 
+
The main benefits of FlexChem are:
+
 
+
#Better documentation of chemical mechanisms;
+
#Easier to drop in other chemical mechanisms;
+
#Optimized chemistry computations; and
+
#Removal of the old SMVGEAR solver (used prior to GEOS-Chem v11-01).
+
 
+
== Requirements ==
+
 
+
FlexChem/KPP requires the following:
+
 
+
#A C language compiler (such as gcc, from the GNU Compiler Collection)
+
#The <tt>flex</tt> library.  This is often installed on many computer systems, or can be easily installed with a package manager such as [https://github.com/spack/spack Spack].
+
 
+
Depending on your setup, you might have to load these packages with the <tt>module load</tt> or <tt>spack load</tt> commands.  Ask your sysadmin for more information.
+
 
+
== Quick start ==
+
 
+
=== Install KPP locally ===
+
 
+
'''Current KPP version for use with GEOS-Chem: <span style="color:red">2.3.0_gc</span>'''
+
 
+
Download KPP to your home directory with this command:
+
 
+
<nowiki>git clone -b GC_updates https://github.com/geoschem/KPP.git</nowiki>
+
 
+
* The <tt>-b GC_updates</tt> will check out the <tt>GC_updates</tt> branch instead of the <tt>main</tt> branch.
+
* IMPORTANT! Do not download the KPP source code into your GEOS-Chem source code directory. This will avoid confusion with the KPP folder that is already within GEOS-Chem, which contains Fortran files created by KPP that define the chemical mechanism.
+
 
+
=== Build the KPP executable ===
+
 
+
Once you have downloaded KPP and have made sure that all the [[#Requirements|required libraries] exist on your system, you may use these commands to build the KPP executable file:
+
 
+
cd KPP/kpp-code
+
make distclean
+
make all
+
 
+
If the build completes successfully, you will see an executable file named <tt>kpp in the <tt>KPP/kpp-code/bin/</tt> folder.
+
 
+
=== Set environment variables for KPP ===
+
 
+
Once have built KPP, you must add the path to the KPP executable to your Unix <tt>PATH</tt> variable.
+
 
+
If you use the bash Unix shell, add these lines to your <tt>~/.bash_aliases</tt> file.  (If you don't have a <tt>~/.bash_aliases</tt> file, you can add these lines to your <tt>~/.bashrc</tt> file instead.)
+
 
+
export PATH=$PATH:/PATH_TO_KPP/KPP/kpp-code/bin/
+
export KPP_HOME=PATH_TO_KPP/KPP/kpp-code
+
 
+
If you use the csh or tcsh Unix shell, add these lines to your <tt>~/.cshrc</tt> file:
+
 
+
setenv PATH $PATH:/PATH_TO_KPP/KPP/kpp-code/bin/
+
setenv KPP_HOME=PATH_TO_KPP/KPP/kpp-code
+
 
+
NOTES:
+
*In the examples above, <tt>PATH_TO_KPP</tt> is the path to the top-level KPP folder.
+
*For example, if you installed KPP into your home directory, then <tt>PATH_TO_KPP</tt> would be <tt>~/KPP</tt>, etc.
+
 
+
=== Use KPP to create Fortran-90 source code for GEOS-Chem ===
+
 
+
At this point you can now use KPP to generate Fortran-90 source code files that will solve your chemical mechanism in an efficient manner.  Navigate to this folder in your GEOS-Chem source code:
+
 
+
* GEOS-Chem 12.9.3 and prior versions: <tt>KPP</tt>
+
* GEOS-Chem 13.0.0 and later versions: <tt>src/GEOS-Chem/KPP</tt>
+
* GCHP 13.0.0 and later versions: <tt>src/GCHP_GridComp/GEOSChem_GridComp/geos-chem/KPP</tt>
+
 
+
Here you will find two sub-folders: <tt>fullchem</tt> and <tt>custom</tt>, and a script named <tt>build_mechanism.sh</tt>.
+
 
+
The <tt>fullchem</tt> folder contains chemical mechanism specification files (<tt>fullchem.eqn</tt> and <tt>gckpp.kpp</tt>) and Fortran-90 solver files for the default "out-of-the-box" GEOS-Chem chemical mechanism.  You should leave these files untouched.  This will allow you to revert to "out-of-the-box" "fullchem" mechanism if the need should arise. 
+
 
+
The <tt>custom</tt> folder contains sample chemical mechanism specification files (<tt>custom.eqn</tt> and <tt>gckpp.kpp</tt>) which have been copied from fullchem.  You can edit these files to [[#Specifying a custom chemical mechanism|define your own custom mechanism (see subsequent sections for detailed instructions)]].
+
 
+
Once you are satisfied with your custom mechanism specification you may now use KPP to build the source code files for GEOS-Chem.  Return to the the KPP folder containing <tt>build_mechanism.sh</tt> and then type:
+
 
+
./build_mechanism.sh custom
+
 
+
You will see output similar to this:
+
 
+
This is KPP-2.3.0_gc.
+
+
KPP is parsing the equation file.
+
KPP is computing Jacobian sparsity structure.
+
KPP is starting the code generation.
+
KPP is initializing the code generation.
+
KPP is generating the monitor data:
+
    - gckpp_Monitor
+
KPP is generating the utility data:
+
    - gckpp_Util
+
KPP is generating the global declarations:
+
    - gckpp_Main
+
KPP is generating the ODE function:
+
    - gckpp_Function
+
KPP is generating the ODE Jacobian:
+
    - gckpp_Jacobian
+
    - gckpp_JacobianSP
+
KPP is generating the linear algebra routines:
+
    - gckpp_LinearAlgebra
+
KPP is generating the utility functions:
+
    - gckpp_Util
+
KPP is generating the rate laws:
+
    - gckpp_Rates
+
KPP is generating the parameters:
+
    - gckpp_Parameters
+
KPP is generating the global data:
+
    - gckpp_Global
+
KPP is generating the driver from none.f90:
+
    - gckpp_Main
+
KPP is starting the code post-processing.
+
+
KPP has succesfully created the model "gckpp".
+
+
+
Reactivity consists of 172 reactions
+
Written to gckpp_Util.F90
+
 
+
If this process is successful, the <tt>custom</tt> folder should now be populated with several <tt>.F90</tt> source code files:
+
 
+
CMakeLists.txt*      gckpp_Initialize.F90  gckpp_LinearAlgebra.F90  gckpp_Precision.F90
+
custom.eqn          gckpp_Integrator.F90  gckpp.map                gckpp_Rates.F90
+
gckpp_Function.F90  gckpp_Jacobian.F90    gckpp_Model.F90          gckpp_Util.F90
+
gckpp_Global.F90    gckpp_JacobianSP.F90  gckpp_Monitor.F90        Makefile_gckpp
+
gckpp_HetRates.F90@  gckpp.kpp            gckpp_Parameters.F90
+
 
+
These files contain optimized instructions for solving the chemical mechanism that you just defined.
+
 
+
=== Tell GEOS-Chem to use your custom mechanism ===
+
 
+
You must explicitly tell GEOS-Chem that it should use the custom mechanism that you just built rather than the default "out-of-the-box" fullchem mechanism.  To do this, you must pass the <tt>-DCUSTOMMECH=y</tt> flag to CMake at configuration time. 
+
 
+
For example, to build GEOS-Chem "Classic" with your custom mechanism, navigate to your run directory, and type:
+
 
+
cd build
+
cmake ../CodeDir -DCUSTOMMECH=y
+
 
+
You should see output such as this written to the screen:
+
 
+
  ... etc before ...
+
-- General settings:
+
  * CUSTOMMECH:  <span style="color:darkgreen">'''ON'''</span>  OFF
+
  ... etc after ...
+
+
This lets you know that GEOS-Chem will use the custom mechanism instead of the "out-of-the-box" fullchem mechanism. Once GEOS-Chem has been configured properly, you may proceed to build the GEOS_Chem executable:
+
 
+
make -j
+
make -j install
+
 
+
and this will build the <tt>gcclassic</tt> executable in the run directory.
+
 
+
The process is the same when building GCHP, make sure to pass <tt>-DCUSTOMMECH=y</tt> at configuration time.
+
 
+
== Specifying a custom chemical mechanism ==
+
 
+
To create a custom mechanism, you will need to edit the following files:
+
 
+
In GEOS-Chem 12.9.3 and prior versions:
+
 
+
# <tt>KPP/custom/custom.eqn</tt>
+
# <tt>KPP/custom/gckpp.kpp</tt>
+
 
+
In GEOS-Chem Classic 13.0.0 and later versions:
+
 
+
# <tt>src/GEOS-Chem/KPP/custom/custom.eqn</tt>
+
# <tt>src/GEOS-Chem/KPP/custom/gckpp.kpp</tt>
+
 
+
In GCHP 13.0.0 and later versions:
+
 
+
# <tt>src/GCHP_GridComp/GEOSChem_GridComp/geos-chem/KPP/custom/custom.eqn</tt>
+
# <tt>src/GCHP_GridComp/GEOSChem_GridComp/geos-chem/KPP/custom/gckpp.kpp</tt>
+
 
+
The <tt>custom.eqn</tt> file specifies the chemical species and reaction list. The <tt>gckpp.kpp</tt> file specifies the choice of solver, language, list of chemical families, and rate-law functions.  The "out-of-the-box" <tt>custom.eqn</tt> and <tt>gckpp.kpp</tt> are copies of the the default GEOS-Chem fullchem mechanism configuration files (<tt>fullchem.eqn</tt> and <tt>gckpp.kpp</tt>).  This will easily let you create your own mechanism using the fullchem mechanism as your starting point.
+
 
+
=== Add species ===
+
 
+
List chemically-active (aka variable) species in the <tt>#DEFVAR</tt> section of <tt>custom.eqn</tt>, as shown below:
+
 
+
#DEFVAR
+
+
A3O2      = IGNORE; {CH3CH2CH2OO; Primary RO2 from C3H8}
+
ACET      = IGNORE; {CH3C(O)CH3; Acetone}
+
ACTA      = IGNORE; {CH3C(O)OH; Acetic acid}
+
...etc ...
+
 
+
List species whose concentrations do not change in the <tt>#DEFFIX</tt> of <tt>custom.eqn</tt>, as shown below:
+
 
+
#DEFFIX
+
+
H2        = IGNORE; {H2; Molecular hydrogen}
+
N2        = IGNORE; {N2; Molecular nitrogen}
+
O2        = IGNORE; {O2; Molecular oxygen}
+
... etc ...
+
 
+
Species may be listed in any order, but we have found it convenient to list them alphabetically.
+
 
+
=== Add gas-phase reactions ===
+
 
+
List gas-phase reactions first in the <tt>#EQUATIONS</tt> section of <tt>custom.eqn</tt>.
+
 
+
#EQUATIONS
+
//
+
// Gas-phase reactions
+
//
+
// NOTES:
+
// ------
+
// (1) Be sure to use "D" exponents to force double precision values!
+
//    (i.e. write 1.70d-12 instead of 1.70e-12, etc.).
+
//        -- Bob Yantosca (16 Dec 2020)
+
//
+
// (2) This file might not render properly if the right hand side of the
+
//    equation is longer than ~100 characters.  This seems to be an issue
+
//    with the KPP code itself.  See this Github issue at geoschem/KPP:
+
//    https://github.com/geoschem/KPP/issues/1
+
//        -- Bob Yantosca (16 Dec 2020)
+
//
+
// (3) To avoid useless CPU cycles, we have introduced new rate law functions
+
//    that skip computing Arrhenius terms (and other terms) that would
+
//    evaluate to 1.  The Arrhenius terms that are passed to the function
+
//    are in most cases now noted in the function name (e.g. GCARR_abc takes
+
//    Arrhenius A, B, C parameters but GCARR_ac only passes A and C
+
//    parameters because B=0 and the (300/T)*B would evaluate to 1).
+
//    This should be much more computationally efficient, as these functions
+
//    are called (sometimes multiple times) for each grid box where we
+
//    perform chemistry.
+
//        -- Bob Yantosca (25 Jan 2020)
+
//
+
//
+
O3 + NO = NO2 + O2 :                        GCARR(3.00E-12, 0.0, -1500.0);
+
O3 + OH = HO2 + O2 :                        GCARR(1.70E-12, 0.0, -940.0);
+
O3 + HO2 = OH + O2 + O2 :                    GCARR(1.00E-14, 0.0, -490.0);
+
O3 + NO2 = O2 + NO3 :                        GCARR(1.20E-13, 0.0, -2450.0);
+
... etc ...
+
 
+
==== General form ====
+
 
+
No matter what reaction is being added, the general procedure is the same. A new line must be added to <tt>custom.eqn</tt> of the following form:
+
 
+
A + B = C + 2.000D : RATE_LAW_FUNCTION(ARG_A, ARG_B ...);
+
 
+
The denotes the reactants (A and B) as well as the products (C and D) of the reaction. If exactly one molecule is consumed or produced, then the factor can be omitted; otherwise the number of molecules consumed or produced should be specified with at least 1 decimal place of accuracy. The final section, between the colon and semi-colon, specifies the function (<tt>RATE_LAW_FUNCTION</tt>) and its arguments which will be used to calculate the reaction rate constant k.  Rate-law functions are specified in the <tt>gckpp.kpp</tt> file.
+
 
+
For an equation such as the one above, the overall rate at which the reaction will proceed is determined by <tt>k[A][B]</tt>.  However, if the reaction rate does not depend on the concentration of A or B, you may write it with a constant value, such as:
+
 
+
A + B = C + 2.000D : 8.95d-17
+
 
+
This will save the overhead of a function call.
+
 
+
==== Rates for two-body reactions according to the Arrhenius law ====
+
 
+
For many reactions, the calculation of k follows the Arrhenius law:
+
 
+
k = a0 + ( 300 / TEMP )**b0 + EXP( c0 / TEMP )
+
 
+
For example, the JPL chemical data evaluation (Feb 2017) specifies that the reaction O3 + NO produces NO2 and O2, and its Arrhenius parameters are A = 3.0x10^-12 and E/R = c0 = 1500.
+
 
+
To specify a two-body reaction whose rate follows the Arrhenius law, you can use the <tt>GCARR</tt> rate-law function, which is defined in <tt>gckpp.kpp</tt>.  For example, the entry for the <tt>O3 + NO = NO2 + O2</tt> reaction can be written as in <tt>custom.eqn</tt> as:
+
 
+
O3 + NO = NO2 + O2 : GCARR(3.00E12, 0.0, -1500.0);
+
 
+
==== Other rate-law functions ====
+
 
+
The <tt>gckpp.kpp</tt> file contains other rate law functions, such as those required for three-body, pressure-dependent reactions. Any rate function which is to be referenced in the <tt>KPP/custom/custom.eqn</tt> file must be available in <tt>gckpp.kpp</tt> prior to building the reaction mechanism.
+
 
+
==== NEAR-FUTURE UPDATE: UPDATING RATE-LAW FUNCTIONS TO AVOID COMPUTING TERMS THAT EVALUATE TO 1 ====
+
 
+
In an early version 13 release, we will use optimized rate-law functions that to avoid computing terms that evaluate to 1.  This is more computationally efficient and will eliminate wasted CPU clock cycles. 
+
 
+
For example, the rate for the <tt>O3 + NO =  NO2 + O2</tt> reaction mentioned above can be computed as:
+
 
+
k = 3.0x10^-12 + EXP( 1500 / TEMP )
+
 
+
Note that the term <tt>(300/TEMP)**b0</tt> evaluates to 1 and thus can be omitted. Because the <tt>EXP</tt> and <tt>**</tt> mathematical operations are costly in terms of CPU clock cycles, computing terms that would only evaluate to 1 should be avoided as much as possible
+
 
+
To optimize rate-law functions, we have created a set of parallel functions that can be called depending on which arguments are nonzero.  For example, the Arrhenius law function <tt>GCARR</tt> will be split into multiple functions:
+
 
+
#<tt>GCARR_abc(a0, b0, c0)</tt>: Use when a0 > 0 and b0 > 0 and c0 > 0
+
#<tt>GCARR_ab(a0, b0)</tt>: Use when a0 > 0 and b0 > 0
+
#<tt>GCARR_ac(a0, c0)</tt>: Use when a0 > 0 and c0 > 0
+
 
+
Thus we can write the O3 + NO reaction in <tt>custom.eqn</tt> as:</blockquote>
+
 
+
      O3 + NO = NO2 + O2 : GCARR_ac(3.00d12, -1500.0d0);
+
 
+
using the rate law function for when both a0 and c0 are nonzero.
+
 
+
Tests have shown that using these optimized rate functions results in a speedup of about 20 minutes for our 4&deg; x 5&deg; GEOS-Chem Classic benchmark simulation.
+
 
+
=== Add heterogeneous reactions ===
+
 
+
List heterogeneous reactions after all of the gas-phase reactions in <tt>custom.eqn</tt>, according to the format below:
+
 
+
//
+
// Heterogeneous reactions
+
//
+
HO2 = O2 :                                  HET(ind_HO2,1);                      {2013/03/22; Paulot2009; FP,EAM,JMAO,MJE}
+
NO2 = 0.500HNO3 + 0.500HNO2 :                HET(ind_NO2,1);
+
NO3 = HNO3 :                                HET(ind_NO3,1);
+
NO3 = NIT :                                  HET(ind_NO3,2);                      {2018/03/16; XW}
+
... etc ...
+
 
+
Implementing new heterogeneous chemistry requires an additional step. For the reaction in question, a reaction should be added as usual, but this time the rate function should be given as an entry in the <tt>HET</tt> array. A simple example is uptake of HO2, specified as
+
 
+
HO2 = O2 : HET(ind_HO2,1);
+
 
+
Note that the product in this case, O2, is actually a fixed species, so no O2 will actually be produced. O2 is used in this case only as a dummy product to satisfy the KPP requirement that all reactions have at least one product.  Here, <tt>HET</tt> is simply an array of pre-calculated rate constants. The rate constants in <tt>HET</tt> are actually calculated in <tt>gckpp_HetRates.F90</tt>.
+
 
+
To implement an additional heterogeneous reaction, the rate calculation must be added to this file. The following example illustrates a (fictional) heterogeneous mechanism which converts the species XYZ into CH2O. This reaction is assumed to take place on the surface of all aerosols, but not cloud droplets (this requires additional steps not shown here). Three steps would be required:
+
 
+
#Add a new line to the <tt>custom.eqn</tt> file, such as <tt>XYZ = CH2O : HET(ind_XYZ,1);</tt><br><br>
+
#Add a new function to <tt>gckpp_HetRates.F90</tt> designed to calculate the heterogeneous reaction rate. As a simple example, we can copy the function <tt>HETNO3</tt> and rename it <tt>HETXYZ</tt>. This function accepts two arguments: molecular mass of the impinging gas-phase species, in this case XYZ, and the reaction's "sticking coefficient" - the probability that an incoming molecule will stick to the surface and undergo the reaction in question. In the case of <tt>HETNO3</tt>, it is assumed that all aerosols will have the same sticking coefficient, and the function returns a first-order rate constant based on the total available aerosol surface area and the frequency of collisions.<br><br>
+
#Add a new line to the function <tt>SET_HET</tt> in <tt>gckpp_HetRates.F90</tt> which calls the new function with the appropriate arguments and passes the calculated constant to <tt>HET</tt>. Example: assuming a molar mass of 93 g/mol, and a sticking coefficient of 0.2, we would write <tt>HET(ind_XYZ,  1) = HETXYZ(9.30E1_fp, 2E-1_fp)</tt>
+
 
+
The function <tt>HETXYZ</tt> can then be specialized to distinguish between aerosol types, or extended to provide a second-order reaction rate, or whatever the user desires.
+
 
+
=== Add photolysis reactions ===
+
 
+
List photolysis reactions after the heterogeneous reactions, as shown below.
+
 
+
//
+
// Photolysis reactions
+
//
+
O3 + hv = O + O2 :                          PHOTOL(2);      {2014/02/03; Eastham2014; SDE}
+
O3 + hv = O1D + O2 :                        PHOTOL(3);      {2014/02/03; Eastham2014; SDE}
+
O2 + hv = 2.000O :                          PHOTOL(1);      {2014/02/03; Eastham2014; SDE}
+
... etc ...
+
NO3 + hv = NO2 + O :                        PHOTOL(12);    {2014/02/03; Eastham2014; SDE}
+
... etc ...
+
 
+
A photolysis reaction can be specified by giving the correct index of the <tt>PHOTOL</tt> array. This index can be determined by inspecting the file [[FAST-JX_v7.0_photolysis_mechanism#FJX_j2j.dat|<tt>FJX_j2j.dat</tt>]] file. 
+
 
+
* NOTE: See the [[The_input.geos_file#Photolysis|<tt>PHOTOLYSIS MENU</tt> section of the <tt>input.geos</tt>]] file in your run directory for the folder in which <tt>FJX_j2j.dat</tt> is located).
+
 
+
For example, one branch of the NO3 photolysis reaction is specified in the <tt>KPP/custom/custom.eqn</tt> file as
+
 
+
NO3 + hv = NO2 + O : PHOTOL(<span style="color:red">12</span>)
+
 
+
Referring back to <tt>FJX_j2j.dat</tt> shows that reaction 12, as specified by the left-most index, is indeed NO3 = NO2 + O:
+
 
+
  <span style="color:red">12</span> NO3      PHOTON    NO2      O                      0.886 /NO3  /
+
 
+
If your reaction is not already in <tt>FJX_j2j.dat</tt>, you may add it there. You may also need to modify <tt>FJX_spec.dat</tt> (in the same folder ast <tt>FJX_j2j.dat</tt>) to include cross-sections for your species. Note that if you add new reactions to <tt>FJX_j2j.dat</tt> you will also need to set the parameter <tt>JVN_</tt> in module <tt>Headers/CMN_FJX_MOD.F</tt> to match the total number of entries.
+
 
+
If your reaction involves new cross section data, you will need to follow an additional set of steps. Specifically, you will need to:
+
 
+
#Estimate the cross section of each wavelength bin (using the correlated-k method), and
+
#Add this data to the <tt>FJX_spec.dat</tt> file.
+
 
+
For the first step, you can use tools already available on the Prather research group website. To generate the cross-sections used by Fast-JX, download the file "UCI_fastJ_addX_73cx.zip" from: [ftp://128.200.14.8/public/prather/Fast-J_&_Cloud-J/]. You can then simply add your data to <tt>FJX_spec.dat</tt> and refer to it in <tt>FJX_j2j.dat</tt> as specified above. The following then describes how to generate a new set of cross-section data for the example of some new species MEKR:
+
 
+
To generate the photolysis cross sections of a new species, come up with some unique name which you will use to refer to it in the <tt>FJX_j2j.dat</tt> and <tt>FJX_spec.dat</tt> files - e.g. MEKR. You will need to copy one of the addX_*.f routines and make your own (say, addX_MEKR.f). Your edited version will need to read in whatever cross section data you have available, and you'll need to decide how to handle out-of-range information - this is particularly crucial if your cross section data is not defined in the visible wavelengths, as there have been some nasty problems in the past caused by implicitly assuming that the XS can be extrapolated (I would recommend buffering your data with zero values at the exact limits of your data as a conservative first guess). Then you need to compile that as a standalone code and run it; this will spit out a file fragment containing the aggregated 18-bin cross sections, based on a combination of your measured/calculated XS data and the non-contiguous bin subranges used by Fast-JX. Once that data has been generated, just add it to <tt>FJX_spec.dat</tt> and refer to it as above. There are examples in the addX files of how to deal with variations of cross section with temperature or pressure, but the main takeaway is that you will generate multiple cross section entries to be added to <tt>FJX_spec.dat</tt> with the same name.
+
 
+
An important complication: if your cross section data varies as a function of temperature AND pressure, you need to do something a little different. The acetone XS documentation shows one possible way to handle this; Fast-JX currently interpolates over either T or P, but not both, so if your data varies over both simultaneously then this will take some thought. The general idea seems to be that one determines which dependence is more important and uses that to generate a set of 3 cross sections (for interpolation), assuming values for the unused variable based on the standard atmosphere.
+
 
+
=== Modifying the gckpp.kpp file ===
+
 
+
(2) Modify the <tt>gckpp.kpp</tt> file if needed. This file defines the KPP integrator and tells KPP to use the specified <tt>.eqn</tt> file to build the mechanism. This file also defines the prod/loss families as [[FlexChem#Production_.26_loss_diagnostic|described below]].
+
#Rebuild the mechanism with KPP by navigating to the top-level KPP directory and typing <tt>./build_mechanism.sh Custom</tt>. The <tt>build_mechanism.sh</tt> script will call on KPP to build the <tt>gckpp_*.F90</tt> files. You may choose to rebuild the other supported chemical mechanisms by using <tt>./build_mechanism.sh [NAME]</tt>, but we recommend testing with the Custom mechanism first to avoid breaking the other mechanisms.
+
#Configure GEOS-Chem for the new mechanism with <
+
 
+
=== Add production & loss families (if desired) ===
+
 
+
Certain common families (e.g. POx, LOx) have been pre-defined for you.  You will find the family definitions near the top of the <tt>KPP/custom/gckpp.kpp</tt> file: 
+
 
+
#FAMILIES
+
POx : O3 + NO2 + 2NO3 + PAN + PPN + MPAN + HNO4 + 3N2O5 + HNO3 + BrO + HOBr + BrNO2 + 2BrNO3 + MPN + ETHLN + MVKN + MCRHN + MCRHNB + PROPNN + R4N2 + PRN1 + PRPN + R4N1 + HONIT + MONITS + MONITU + OLND + OLNN + IHN1 + IHN2 + IHN3 + IHN4 + INPB + INPD + ICN + 2IDN + ITCN + ITHN + ISOPNOO1 + ISOPNOO2 + INO2B + INO2D + INA + IDHNBOO + IDHNDOO1 + IDHNDOO2 + IHPNBOO + IHPNDOO + ICNOO + 2IDNOO + MACRNO2 + ClO + HOCl + ClNO2 + 2ClNO3 + 2Cl2O2 + 2OClO + O + O1D + IO + HOI + IONO + 2IONO2 + 2OIO + 2I2O2 + 3I2O3 + 4I2O4;
+
LOx : O3 + NO2 + 2NO3 + PAN + PPN + MPAN + HNO4 + 3N2O5 + HNO3 + BrO + HOBr + BrNO2 + 2BrNO3 + MPN + ETHLN + MVKN + MCRHN + MCRHNB + PROPNN + R4N2 + PRN1 + PRPN + R4N1 + HONIT + MONITS + MONITU + OLND + OLNN + IHN1 + IHN2 + IHN3 + IHN4 + INPB + INPD + ICN + 2IDN + ITCN + ITHN + ISOPNOO1 + ISOPNOO2 + INO2B + INO2D + INA + IDHNBOO + IDHNDOO1 + IDHNDOO2 + IHPNBOO + IHPNDOO + ICNOO + 2IDNOO + MACRNO2 + ClO + HOCl + ClNO2 + 2ClNO3 + 2Cl2O2 + 2OClO + O + O1D + IO + HOI + IONO + 2IONO2 + 2OIO + 2I2O2 + 3I2O3 + 4I2O4;
+
PCO : CO;
+
LCO : CO;
+
PSO4 : SO4;
+
LCH4 : CH4;
+
PH2O2 : H2O2;
+
 
+
NOTES:
+
#The Ox and CO rates are used in GEOS-Chem for computing budgets in the [http://acmg.seas.harvard.edu/geos/geos_benchmark.html 1-month benchmark simulations].
+
#PSO4 is required for simulations using [[TOMAS_aerosol_microphysics|TOMAS aerosol microphysics]].
+
 
+
To add a new prod/loss family, add a new line to the <tt>#FAMILIES</tt> section with the format
+
 
+
FAM_NAME : MEMBER_1 + MEMBER_2 + ... + MEMBER_N;
+
 
+
The family name must start with <tt>P</tt> or <tt>L</tt> to indicate whether KPP should calculate a production or a loss rate.
+
 
+
The maximum number of families allowed by KPP is currently set to 300. Depending on how many prod/loss families you add, you may need to increase that to a larger number to avoid errors in KPP. You can change the number for <tt>MAX_FAMILIES</tt> in <tt>KPP/kpp-code/src/gdata.h</tt> and then [[FlexChem#KPP_source_code|rebuild KPP]].
+
 
+
#define MAX_EQN        1500    /* KPP 2.3.0_gc, Bob Yantosca (11 Feb 2021)  */
+
#define MAX_SPECIES    1000    /* KPP 2.3.0_gc, Bob Yantosca (11 Feb 2021)  */
+
#define MAX_SPNAME      30
+
#define MAX_IVAL        40
+
#define MAX_EQNTAG      12    /* Max length of equation ID in eqn file    */
+
#define MAX_K          150    /* Max length of rate expression in eqn file */
+
#define MAX_ATOMS        10
+
#define MAX_ATNAME      10
+
#define MAX_ATNR        250
+
#define MAX_PATH        120
+
#define MAX_FILES        20
+
#define MAX_FAMILIES    300
+
#define MAX_MEMBERS    150
+
#define MAX_EQNLEN      200
+
 
+
<span style="color:red">'''IMPORTANT:''' When adding a prod/loss family or changing any of the other settings in <tt>gckpp.kpp</tt>, the chemistry mechanism will need to be rebuilt with KPP as [[#Building a custom chemical mechanism|described above]].</span>
+
 
+
--[[User:Melissa Payer|Melissa Sulprizio]] ([[User talk:Melissa Payer|talk]]) 17:13, 9 November 2016 (UTC)
+
 
+
=== Saving out production & loss rates ===
+
 
+
==== NetCDF format ====
+
 
+
[[GEOS-Chem v11-02]] introduces the option to save GEOS-Chem diagnostics to netCDF format by compiling with <tt>NC_DIAG=y</tt>. To save out the prod/loss diagnostics to netCDF format, you can add the following field names to your defined collection(s) in <tt>HISTORY.rc</tt>:
+
 
+
  'Prod_?PRD?',                  'GIGCchem', 
+
  'Loss_?LOS?',                  'GIGCchem', 
+
 
+
where <tt>?PRD?</tt> and <tt>?LOS?</tt> are wildcards that will be placed with all production and loss species as defined in GEOS-Chem. '''''NOTE: GCHP doesn't accept wildcards at this time, so each prod/loss species will need to be listed separately.'''''
+
 
+
--[[User:Melissa Payer|Melissa Sulprizio]] ([[User talk:Melissa Payer|talk]]) 17:21, 6 February 2018 (UTC)
+
 
+
== Previous issues that are now resolved ==
+
 
+
=== FlexChem bug fix: do not zero ACTA, EOH, HCOOH ===
+
 
+
<span style="color:green">'''''This fix was included in [[GEOS-Chem 12#12.0.0|GEOS-Chem 12.0.0]].'''''</span>
+
 
+
'''''[[User:Katherine Travis|Katie Travis]] wrote:'''''
+
 
+
<blockquote>I am working on a VOC simulation, and noticed that in my copy of v11-02f, the following species are set to zero in two places:</blockquote>
+
 
+
      ! Zero certain species
+
      C(ind_ACTA) = 0.e0_dp
+
      C(ind_EOH) = 0.e0_dp
+
      C(ind_HCOOH) = 0.e0_dp
+
 
+
<blockquote>And</blockquote>
+
 
+
  C(ind_ACTA)  = 0.0_dp
+
  C(ind_HCOOH) = 0.0_dp
+
 
+
<blockquote>Since none of these species are fixed in Tropchem.eqn, shouldn’t they NOT be set to zero?</blockquote>
+
 
+
'''''Mike Long wrote:'''''
+
 
+
<blockquote>I think the code should be removed.  This must have been a patch added to maintain parity with SMVGEAR w/o anticipating that the species would become active.</blockquote>
+
 
+
--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 16:19, 17 May 2018 (UTC)
+
 
+
=== Remove calls to UPDATE_SUN, UPDATE_RCONST from gckpp_Integrator.F90 ===
+
 
+
<span style="color:green">'''''This update was included in [[GEOS-Chem v11-01 benchmark history#v11-01 provisional release|GEOS-Chem v11-01 provisional release]]'''''</span>
+
 
+
A slow down in GEOS-Chem run time was observed following the implementation of FlexChem in [[GEOS-Chem v11-01|v11-01]]. To resolve this, a temporary workaround was implemented. This fix may be replaced with a more permanent solution in [[GEOS-Chem v11-01#v11-01 public release|GEOS-Chem v11-01 public release]].
+
 
+
'''''[[User:Bmy|Bob Yantosca]] wrote:'''''
+
 
+
<blockquote>KPP automatically places calls to <tt>UPDATE_SUN</tt> and <tt>UPDATE_RCONST</tt> in the routines <tt>FunTemplate</tt> and <tt>JacTemplate</tt> (both in the <tt>gckpp_Integrator.F90</tt> module).  This assumes that you are not interfacing KPP into any other model, and that you will use KPP to compute the sun angles at each timestep.
+
   
+
We now call <tt>UPDATE_RCONST</tt> once per grid box before calling the KPP integrator.  Also, because we use FAST-JX to get the photo rates, we no longer need to call <tt>UPDATE_SUN</tt>. These duplicate calls were causing a performance bottleneck, as <tt>UPDATE_RCONST</tt> was being called more than 7 million times per day.
+
   
+
We have removed the duplicate calls from the <tt>gckpp_Integrator.F90</tt> modules in each of the chemistry mechanisms.  But we will also need to make sure that when building KPP fresh from an equation file, that this step gets automatically added to the build sequence.</blockquote>
+
 
+
--[[User:Melissa Payer|Melissa Sulprizio]] ([[User talk:Melissa Payer|talk]]) 19:53, 9 January 2017 (UTC)
+
 
+
== Known issues ==
+
 
+
None at this time.
+

Latest revision as of 20:24, 4 August 2022

This content has been migrated to: