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

From Geos-chem
Jump to: navigation, search
m (Generating KPP input files from GEOS-Chem globchem.dat)
(Interfacing the generated Code into GEOS-Chem: list of required modifications)
Line 74: Line 74:
  
 
== Interfacing the generated Code into GEOS-Chem ==
 
== Interfacing the generated Code into GEOS-Chem ==
(coming soon)
+
 
 +
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
 +
    -------------------------------------------------
 +
   
 +
    (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_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 tracerssimulation .
 +
 
 +
    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.
 +
 +
    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)
 
--[[User:Plesager|phs]] 14:23, 16 October 2009 (EDT)

Revision as of 17:50, 20 November 2009

This page describes how to install and run KPP, and how to interface it with GEOS-Chem.

Installing KPP and Running the examples

The source code, examples and the documentation to install and run KPP are available directly from 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), one for the new format used by Fabien Paulot for the alternate isoprene chemistry (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 


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 
   -------------------------------------------------
   
   (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_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 tracerssimulation .
 
   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. 

   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>


--phs 14:23, 16 October 2009 (EDT)