Interfacing GEOS-Chem with KPP
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 wiki page.
Note: The interface between KPP and GEOS-Chem was changed in version 8-02-04. The differences implied in how to interface GEOS-Chem with KPP between (1) version 8-02-03 and (2) versions 8-02-04 and higher are explicitly given in the adequate paragraphs.
Contents
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 (*)
(*) 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 improved HO2 uptake by aerosol update that went into GEOS-Chem v9-02 introduced the following reaction in globchem.dat:
A 0 3.30E+01 2.0E-01 0 0 K 0.00 0. 0.
HO2 +
= + + +
+ + + +
+ + + +
+ + + +
The 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 globchem.eqn 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 globchem.spc file:
DUMMY = IGNORE;
Once you have made those modifications, you should be able to run KPP by using kpp gckpp.kpp.
--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 kpp-2.2.1/src/gdata.h from:
#define MAX_EQN 400
to:
#define MAX_EQN 1000
You will need to recompile KPP and rerun using kpp gckpp.kpp. Note that the warnings for duplicate equations are normal and should not be a problem for KPP.
--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 *.f90 to *.F90 (required for GEOS-Chem v9-01-02 and later versions.
- Modify the generated kpp files
- 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 [**] 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 Initialize 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 WLAMCH
!$OMP THREADPRIVATE( Eps, First )
gckpp_Rates.F90
For v8-02-04 and higher
- 1. Comment all rate law functions and update_sun
- 2. 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
- 3. Reduce UPDATE_PHOTO to a stub:
SUBROUTINE Update_PHOTO ( )
USE gckpp_Global
END SUBROUTINE Update_PHOTO
For v8-02-03
- 1. Comment all rate law functions and update_sun
- 2. 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
- 3. Reduce UPDATE_PHOTO to a stub:
SUBROUTINE Update_PHOTO ( )
USE gckpp_Global
END SUBROUTINE Update_PHOTO
gckpp_Util.F90
NOTE: This modifications is not always needed.
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 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 [Talk:Interfacing_GEOS-Chem_with_KPP|this discussion]].
To resolve this problem:
- Check which elements of V() are not initialized in shuffle_user2kpp
- Identify the species name for that element in the gckpp.map file generated by KPP
- Locate the species in globchem.dat and identify which of the missing indices are associated with dead species and active species. Alternatively, check gckpp_Initialize.f90 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>
--Melissa Sulprizio (talk) 22:16, 20 October 2015 (UTC)
