Difference between revisions of "Interfacing GEOS-Chem with KPP"

From Geos-chem
Jump to: navigation, search
 
(2 intermediate revisions by 2 users not shown)
Line 1: Line 1:
[[Image:obsolete.jpg]]
+
This content has been migrated to [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>]
 
+
<span style="color:red">'''''This page describes how to install KPP and how to interface it with GEOS-Chem v10-01 and earlier versions. For [[GEOS-Chem v11-01]] and later, instructions can be found on the [[FlexChem|''FlexChem'' wiki page]].'''''</span>
+
 
+
'''Note:''' The interface between KPP and GEOS-Chem was changed in  [[GEOS-Chem v8-02-04|version 8-02-04]]. The differences implied in how to interface GEOS-Chem with KPP between (1) [[GEOS-Chem v8-02-03|version 8-02-03]] and (2) versions 8-02-04 and higher are explicitly given in the adequate paragraphs.
+
 
+
== Installing KPP and Running the examples ==
+
The source code, examples and the documentation to install and run KPP are available directly from [http://people.cs.vt.edu/~asandu/Software/Kpp/ A. Sandu web page].
+
 
+
For version 2.2.1 of KPP to work with GEOS-Chem, you must modify the source code so duplicate and proportional reactions in the input files will trigger a warning instead of an error. (Thanks to Paul Eller.)
+
 
+
In the file kpp/src/scanner.c , search for the following line of code around line 665 (may be slightly different):
+
ScanError( "Duplicate equation: "
+
            " (eqn<%d> = eqn<%d> )", i+1, EqnNr+1 );
+
and change it to the following:
+
ScanWarning( "Duplicate equation: "
+
              " (eqn<%d> = eqn<%d> )", i+1, EqnNr+1 );
+
Then replace :
+
        ScanError( "Linearly dependent equations: "
+
                  "( %.0f eqn<%d> = %.0f eqn<%d> )",
+
                  r1, i+1, r2, EqnNr+1 );
+
with
+
        ScanWarning( "Linearly dependent equations: "
+
                    "( %.0f eqn<%d> = %.0f eqn<%d> )",
+
                    r1, i+1, r2, EqnNr+1 );
+
 
+
This should allow you to create the KPP model. Having duplicate and proportional equations should not be a problem for KPP.
+
 
+
 
+
(1) Installation of KPP is simple. You need to double check the location of the flex library. Ask your system administrator, or just use the "find" command.
+
 
+
(2) Running examples should not be a problem either. The kpp executable being on the $path, you just need to be in the directory were the model file (*.kpp) is to run kpp on it. (Some examples names are different between the distribution and the manual, but it is straightforward to figure it out.)
+
 
+
 
+
== Generating KPP input files from GEOS-Chem globchem.dat ==
+
"geos2kpp_parser.pl" is a Perl script from Paul Eller that automatically generates the three input files for KPP from globchem.dat:
+
> geos2kpp_parser.pl globchem.dat
+
 
+
the script is available in two versions, one for the original globchem.dat format ([[geos2kpp_parser.pl|geos2kpp_parser.pl]]), one for the new format used by Fabien Paulot for the alternate isoprene chemistry ([[paulot_geos2kpp_parser.pl|paulot_geos2kpp_parser.pl]]).
+
 
+
this will create 3 files (globchem.def, globchem.eqn, globchem.spc). Then you must use "globchem" for the MODEL flag in the *.kpp file. So create a "gckpp.kpp", with a rosenbrock solver, like that:
+
 
+
#MODEL      globchem
+
#INTEGRATOR rosenbrock
+
#LANGUAGE  Fortran90
+
#DRIVER    none
+
#HESSIAN    on
+
#STOICMAT  on
+
 
+
The "none" driver is enough to get the files needed for GEOS-Chem. We are not generating a standalone code.
+
 
+
 
+
 
+
== Running KPP ==
+
 
+
Put the gckpp.kpp file and the globchem* files in the same directory, and call kpp there:
+
 
+
> kpp gckpp.kpp
+
 
+
this will generate all the necessary files:
+
 
+
gckpp_Function.f90 (*)
+
gckpp_Global.f90  (*)
+
gckpp_Hessian.f90
+
gckpp_HessianSP.f90
+
gckpp_Initialize.f90  (*)
+
gckpp_Integrator.f90
+
gckpp_Jacobian.f90
+
gckpp_JacobianSP.f90
+
gckpp_LinearAlgebra.f90  (*)
+
gckpp_Model.f90
+
gckpp_Monitor.f90
+
gckpp_Parameters.f90
+
gckpp_Precision.f90
+
gckpp_Rates.f90  (*)
+
gckpp_Stoichiom.f90
+
gckpp_StoichiomSP.f90
+
gckpp_Util.f90 (*)
+
 
+
(*) are modified in the next step. The following are also generated but not needed for GEOS-Chem (but keep gckpp.map):
+
 
+
Makefile_gckpp
+
gckpp.map
+
gckpp_mex_Fun.f90
+
gckpp_mex_Hessian.f90
+
gckpp_mex_Jac_SP.f90
+
 
+
To the 17 files above, you must add two more, '''gckpp_comode_mod.f90''' and '''Makefile''', that you can copy from an existing one in a distribution of GEOS-Chem.
+
 
+
'''For v8-02-03:''' Copy the files from the v8-02-03 GEOS-Chem distribution, KPP/43t/gckpp_comode_mod.f90, and KPP/43t/Makefile. They do not need to be modified.
+
 
+
'''For v8-02-04 and higher:''' Copy the files from the current GEOS-Chem distribution, KPP/standard/gckpp_comode_mod.f90, and KPP/standard/Makefile. They do not need to be modified.
+
 
+
=== For GEOS-Chem v9-02 and higher ===
+
 
+
==== Misplaced ';' warning ====
+
 
+
The [[NOx-Ox-HC-aerosol#Improved_HO2_uptake|improved HO2 uptake by aerosol update]] that went into [[GEOS-Chem v9-02]] introduced the following reaction in <tt>globchem.dat</tt>:
+
 
+
A    0 3.30E+01  2.0E-01      0 0 K  0.00    0.    0.       
+
      HO2          +                                                   
+
=                  +                  +                  +           
+
+                  +                  +                  +           
+
+                  +                  +                  +           
+
+                  +                  +                  +
+
 
+
The [http://www.cs.vt.edu/~asandu/Software/Kpp/Download/kpp-2.1_UsersManual.pdf KPP User's Manual] states, "One restriction is that the list of products must not be empty. If you have such a reaction (e.g. the dry deposition of atmospheric species to the surface), you can define a DUMMY species as the product." To get around this issue, find the appropriate HO2 reaction in <tt>globchem.eqn</tt> and change it from:
+
 
+
{272} HO2 =  :                                        EXTARR(3.30E+01, 2.0E-01, 0.0);
+
 
+
to:
+
 
+
{272} HO2 = DUMMY :                                    EXTARR(3.30E+01, 2.0E-01, 0.0);
+
 
+
You will also need to add the following line to the end of the <tt>globchem.spc</tt> file:
+
 
+
    DUMMY          =          IGNORE;
+
 
+
Once you have made those modifications, you should be able to run KPP by using <tt>kpp gckpp.kpp</tt>.
+
 
+
--[[User:Melissa Payer|Melissa Sulprizio]] 17:58, 3 June 2014 (EDT)
+
 
+
==== "Too many equations" error ====
+
 
+
Starting with [[GEOS-Chem v9-02]], we exceed the number of equations that KPP allows. This may lead to the following errors when running KPP:
+
 
+
This is KPP-2.2.1.
+
+
KPP is parsing the equation file.Warning :globchem.eqn:132: Duplicate equation:  (eqn<93> = eqn<94> )
+
Warning :globchem.eqn:447: Duplicate equation:  (eqn<302> = eqn<304> )
+
Warning :globchem.eqn:448: Duplicate equation:  (eqn<303> = eqn<305> )
+
Error : Too many equations
+
Error : Too many equations
+
Error : Too many equations
+
+
Fatal error : 3 errors and 3 warnings encountered
+
 
+
To get past the error, change the following line in <tt>kpp-2.2.1/src/gdata.h</tt> from:
+
 
+
#define MAX_EQN        400
+
 
+
to:
+
 
+
#define MAX_EQN        1000
+
 
+
You will need to recompile KPP and rerun using <tt>kpp gckpp.kpp</tt>. Note that the warnings for duplicate equations are normal and should not be a problem for KPP.
+
 
+
--[[User:Melissa Payer|Melissa Sulprizio]] 11:24, 4 June 2014 (EDT)
+
 
+
== Interfacing the generated Code into GEOS-Chem ==
+
 
+
There are three steps:
+
*Rename all generated kpp files from <tt>*.f90</tt> to <tt>*.F90</tt> (required for [[GEOS-Chem v9-01-02]] and later versions.
+
*[[#Modify kpp-generated files|Modify the generated kpp files]]
+
*[[#Modify GEOS-Chem|Modify the GEOS-Chem code]]
+
 
+
=== Modify kpp-generated files ===
+
 
+
We have two categories of modification: one specifically to interface with GEOS-Chem and the other to work with OpenMP. The modifications are described below.
+
 
+
==== gckpp_Function.F90 ====
+
   
+
:1. After declaration of variable A, add:
+
 
+
    !$OMP THREADPRIVATE( A )
+
 
+
==== gckpp_Global.F90 ====
+
 
+
:1. Add following declarations:
+
 
+
    ! VAR_ADJ - Concentrations of variable species (global) [**]
+
      REAL(kind=dp) :: VAR_ADJ(NVAR)
+
    ! V_CSPEC - Concentrations of variable species (global)
+
      REAL(kind=dp) :: V_CSPEC(NVAR)
+
    ! V_CSPEC_ADJ - Concentrations of variable species (global) [**]
+
      REAL(kind=dp) :: V_CSPEC_ADJ(NVAR)
+
   
+
    ! NJ - Number of cost function being evaluated [**]
+
      INTEGER, PARAMETER :: NJ = 1
+
    ! NTT - Total number of tropospheric grid cells
+
      INTEGER :: NTT
+
    ! JLOOP - Total number of tropospheric grid cells
+
      INTEGER :: JLOOP
+
    ! SMAL2 - Parameter for insuring positive tracer values, same as in reader.f
+
      REAL(kind=dp), PARAMETER :: SMAL2 = 1.0d-99
+
    ! NCOEFF - Number of reaction rate coeff adjoints [**]
+
      INTEGER, PARAMETER :: NCOEFF = 24
+
    ! VAR_R_ADJ - Concentrations of reaction rate adjoint (global) [**]
+
      REAL(kind=dp) :: VAR_R_ADJ(NCOEFF)
+
    ! JCOEFF - Reaction numbers for each (define in INIT_KPP) [**]
+
      INTEGER :: JCOEFF(NCOEFF)
+
    ! IND - Reaction numbers for each (define in INIT_KPP)
+
      INTEGER :: IND(NREACT)
+
 
+
:Note: variables indicated with <tt>[**]</tt> are for the adjoint model only, and can be deleted... or kept if the adjoint integrator may be used.
+
   
+
:2. Add the following lines to work with OpenMP:
+
 
+
    !$OMP THREADPRIVATE(VAR,VAR_ADJ,VAR_R_ADJ,V_CSPEC,V_CSPEC_ADJ, C )
+
    !$OMP THREADPRIVATE(FIX,JLOOP,RCONST,TIME)
+
   
+
    ! Move stack_ptr here and make THREADPRIVATE for OMP parallelization (dkh, 07/28/09)
+
      INTEGER :: stack_ptr = 0 ! last written entry
+
    !$OMP THREADPRIVATE( stack_ptr )
+
 
+
:3. Comment out the equivalence for OpenMP:
+
 
+
    ! VAR, FIX are chunks of array C
+
          ! EQUIVALENCE( C(1),VAR(1) )
+
          ! EQUIVALENCE( C(88),FIX(1) )
+
 
+
==== gckpp_Initialize.F90 ====
+
 
+
:1. Rewrite subroutine <tt>Initialize</tt> as:
+
 
+
    SUBROUTINE Initialize ( )
+
     
+
      USE gckpp_Global
+
      USE gckpp_Util,    ONLY : Shuffle_user2kpp
+
      USE gckpp_Monitor
+
     
+
      INTEGER :: i
+
   
+
      CALL Shuffle_user2kpp(V_CSPEC,VAR)
+
   
+
      DO i = 1, NFIX
+
        FIX(i) = 1.d0
+
      END DO
+
   
+
      ! these two loops are for the adjoint only
+
      DO I =1, NVAR
+
        C(I)=VAR(I)
+
      ENDDO
+
   
+
      DO I = 1, NFIX
+
        C(NVAR+I) = FIX(I)
+
      END DO
+
           
+
    END SUBROUTINE Initialize   
+
 
+
==== gckpp_LinearAlgebra.F90 ====
+
 
+
:1. Add THREADPRIVATE for the two SAVE variables of the function <tt>WLAMCH</tt>
+
 
+
    !$OMP THREADPRIVATE( Eps, First )
+
 
+
==== gckpp_Rates.F90 ====
+
 
+
===== For v8-02-04 and higher =====
+
   
+
:1. Comment all rate law functions and <tt>update_sun</tt>
+
 
+
:2. Reduce <tt>UPDATE_RCONST</tt> to:
+
 
+
          SUBROUTINE Update_RCONST (R_KPP)
+
         
+
            USE gckpp_Monitor
+
         
+
            ! INPUT ARGUMENT:
+
            REAL*8, INTENT(IN) :: R_KPP(:,:)
+
         
+
            INTEGER :: N
+
         
+
            DO N = 1, NREACT
+
              RCONST(N) = R_KPP(JLOOP,IND(N))
+
            END DO
+
               
+
          END SUBROUTINE Update_RCONST
+
 
+
:3. Reduce <tt>UPDATE_PHOTO</tt> to a stub:
+
 
+
              SUBROUTINE Update_PHOTO ( )
+
             
+
                USE gckpp_Global
+
                   
+
              END SUBROUTINE Update_PHOTO
+
 
+
===== For v8-02-03 =====
+
 
+
:1. Comment all rate law functions and <tt>update_sun</tt>
+
   
+
:2. Reduce <tt>UPDATE_RCONST</tt> to:
+
 
+
          SUBROUTINE Update_RCONST ( )
+
         
+
            USE gckpp_COMODE_MOD, ONLY : R_KPP
+
            USE gckpp_Monitor
+
         
+
            INTEGER :: N
+
         
+
            DO N = 1, NREACT
+
              RCONST(N) = R_KPP(JLOOP,IND(N))
+
            END DO
+
               
+
          END SUBROUTINE Update_RCONST
+
 
+
:3. Reduce <tt>UPDATE_PHOTO</tt> to a stub:
+
 
+
              SUBROUTINE Update_PHOTO ( )
+
             
+
                USE gckpp_Global
+
                   
+
              END SUBROUTINE Update_PHOTO
+
 
+
==== gckpp_Util.F90 ====
+
 
+
'''NOTE: This modifications is not always needed.'''
+
 
+
'''''[[User:Plesager|Philippe Le Sager]] wrote:'''''
+
 
+
:[This modification] was needed for the standard 43 tracers and for the Fabien Paulot alternative isoprene chemistry, but not for the 54 tracers simulation .
+
 
+
:The problem was an initialization issue. Some tracers are not initialized in subroutine <tt>shuffle_user2kpp</tt> of <tt>gckpp_Util.f90</tt>. The non-initialization was leading to different results with 1 and 4 processors. The fix consists in forcing the initialization of dead species to 1d-99, and of active species to their original value in the <tt>gckpp_Initialize.f90</tt>. See [Talk:Interfacing_GEOS-Chem_with_KPP|this discussion]].
+
+
To resolve this problem:
+
* Check which elements of <tt>V()</tt> are not initialized in <tt>shuffle_user2kpp</tt>
+
* Identify the species name for that element in the <tt>gckpp.map</tt> file generated by KPP
+
* Locate the species in <tt>globchem.dat</tt> and identify which of the missing indices are associated with dead species and active species. Alternatively, check <tt>gckpp_Initialize.f90</tt> to see which indices are not initialized (dead).
+
+
In [[GEOS-Chem v8-03-02]] the following solutions were applied:
+
+
:1. For the standard 43 tracers simulation:
+
 
+
    ! Tracer #15 is LISOPOH and it's initialization is missing.
+
    ! It's a dead tracer for this scheme so initialized at 1d-99
+
    ! (ccc, 9/15/10)
+
    V(15)=1d-99
+
+
:2. For Fabien Paulot's alternate isoprene chemistry:
+
 
+
    V(41, 43) = 1d-99
+
    V(50, 93, 98) = 1d-20
+
 
+
=== Modify GEOS-Chem ===
+
 
+
Two scenario are possible:
+
 
+
(1) '''Modify an existing simulation:''' Overwrite the 19 files in the corresponding KPP subdirectory (e.g., KPP/standard/ or KPP/soa/).
+
 
+
(2) '''Add a new simulation or option that uses a new globchem.dat:''' Create a new subdirectory in the KPP directory, and put the 19 KPP files there. You can use the new simulation/option by calling :
+
 
+
    make CHEM=<subdirectory name>
+
 
+
--[[User:Melissa Payer|Melissa Sulprizio]] ([[User talk:Melissa Payer|talk]]) 22:16, 20 October 2015 (UTC)
+

Latest revision as of 20:25, 4 August 2022

This content has been migrated to Update chemical mechanisms with KPP at geos-chem.readthedocs.io