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

From Geos-chem
Jump to: navigation, search
(For GEOS-Chem v9-02 and higher)
 
(11 intermediate revisions by 2 users not shown)
Line 1: Line 1:
This page describes how to install and run KPP, and how to interface it with GEOS-Chem.
+
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>]
 
+
'''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.
+
 
+
 
+
== Interfacing the generated Code into GEOS-Chem ==
+
 
+
There is two steps: modifying the generated kpp files, and modifying 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 list of modifications is:
+
 
+
    -------------------------------------------------
+
    gckpp_Function.f90
+
    -------------------------------------------------
+
   
+
    after declaration of variable A, add:
+
   
+
    !$OMP THREADPRIVATE( A )
+
     
+
 
+
    -------------------------------------------------
+
    gckpp_Global.f90
+
    -------------------------------------------------
+
   
+
    (i) 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)
+
   
+
   
+
    (ii) and 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 )
+
   
+
   
+
   
+
    (iii) and comment the equivalence for OpenMP:
+
   
+
    ! VAR, FIX are chunks of array C
+
          ! EQUIVALENCE( C(1),VAR(1) )
+
          ! EQUIVALENCE( C(88),FIX(1) )
+
   
+
   
+
    Note: variables indicated with [**] are for the adjoint only, and
+
    can be deleted... or kept if the adjoint integrator may be used.   
+
 
+
    -------------------------------------------------
+
    gckpp_Initialize.f90
+
    -------------------------------------------------
+
    rewrote subroutine 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
+
    -------------------------------------------------
+
   
+
    add THREADPRIVATE for the two SAVE variables of the
+
    function WLAMCH
+
   
+
    !$OMP THREADPRIVATE( Eps, First )
+
 
+
    -------------------------------------------------
+
    gckpp_Rates.f90 '''for v8-02-03'''
+
    -------------------------------------------------
+
   
+
    (i) comment all rate law functions and update_sun
+
   
+
    (ii) reduce UPDATE_RCONST 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
+
   
+
   
+
   
+
    (iii) reduce UPDATE_PHOTO to a stub:
+
   
+
   
+
              SUBROUTINE Update_PHOTO ( )
+
             
+
                USE gckpp_Global
+
                   
+
              END SUBROUTINE Update_PHOTO
+
 
+
    -------------------------------------------------
+
    gckpp_Rates.f90 '''for v8-02-04 and higher'''
+
    -------------------------------------------------
+
   
+
    (i) comment all rate law functions and update_sun
+
   
+
    (ii) reduce UPDATE_RCONST 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
+
   
+
   
+
   
+
    (iii) reduce UPDATE_PHOTO to a stub:
+
   
+
   
+
              SUBROUTINE Update_PHOTO ( )
+
             
+
                USE gckpp_Global
+
                   
+
              END SUBROUTINE Update_PHOTO
+
 
+
    -------------------------------------------------
+
    gckpp_Util.f90
+
    -------------------------------------------------
+
+
    The modifications here is not always needed. It 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 shuffle_user2kpp of
+
    gckpp_Util.f90. 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 gckpp_Initialize.f90. See this [[Talk:Interfacing_GEOS-Chem_with_KPP|discussion]].
+
+
    To proceed: (i) check which element of V() is not
+
    initialized in shuffle_user2kpp, then (ii) by cross
+
    referencing globchme.dat and the kpp-generated gckpp.map
+
    file, check which ones of the missing index are dead and
+
    which ones are active, or (ii) check gckpp_Initialize.f90
+
    to see which index is not initialized (dead).
+
+
    Here are the solutions we used so far:
+
+
    For the standard 43 tracers simulation, we need:
+
+
        V(15)=1d-99
+
+
    For Fabien Paulot alternate isoprene chemistry,
+
 
+
        V(41, 43) = 1d-99
+
        V(50, 93, 98) = 1d-20
+
 
+
=== Modify GEOS-Chem ===
+
Most of it is already done. Now, depending on what kind of globchem.dat you are modifying, several scenario are possible:
+
 
+
(1) you are ''modifying'' an existing simulation. Then overwrite the 19 files in the corresponding subdirectory in the KPP directory of the GEOS-Chem distribution (e.g., KPP/standard/ or KPP/soa/). That's it.
+
 
+
(2) you are ''adding'' a new simulation or option that uses a new globchem.dat. You need to create a new subdirectory in the KPP directory, and put the 19 kpp files there. That's it, you can use the new simulation/option by calling :
+
 
+
    make CHEM=<subdirectory name>
+
 
+
 
+
 
+
--[[User:Plesager|phs]] 14:23, 16 October 2009 (EDT)
+

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