Input for KPP

KPP basically handles two types of input files: Kinetic description files and auxiliary files. Kinetic description files are in KPP syntax and described in the following sections. Auxiliary files are described in the section entitled Auxiliary files and the substitution preprocessor.

KPP kinetic description files specify the chemical equations, the initial values of each of the species involved, the integration parameters, and many other options. The KPP preprocessor parses the kinetic description files and generates several output files. Files that are written in KPP syntax have one of the suffixes .kpp, .spc, .eqn, or def.

The following general rules define the structure of a kinetic description file:

A KPP program is composed of KPP sections, KPP commands, and Inlined Code. Their syntax is presented in BNF description of the KPP language.
Comments are either enclosed between the curly braces { and }, or written in a line starting with two slashes //.
Any name given by the user to denote an atom or a species is restricted to be less than 32 character in length and can only contain letters, numbers, or the underscore character. The first character cannot be a number. All names are case insensitive.

The kinetic description files contain a detailed specification of the chemical model, information about the integration method and the desired type of results. KPP accepts only one of these files as input, but using the #INCLUDE command, code from separate files can be combined. The include files can be nested up to 10 levels. KPP will parse these files as if they were a single big file. By carefully splitting the chemical description, KPP can be configured for a broad range of users. In this way the users can have direct access to that part of the model that they are interested in, and all the other details can be hidden inside several include files. Often, the atom definitions (atoms.kpp) are included first, then species definitions (*.spc), and finally the equations of the chemical mechanism (*.eqn).

KPP sections

A # sign at the beginning of a line followed by a section name starts a new KPP section. Then a list of items separated by semicolons follows. A section ends when another KPP section or command occurs, i.e. when another # sign occurs at the beginning of a line. The syntax of an item definition is different for each particular section.

#ATOMS

The atoms that will be further used to specify the components of a species must be declared in an #ATOMS section, e.g.:

#ATOMS N; O; Na; Br;

Usually, the names of the atoms are the ones specified in the periodic table of elements. For this table there is a predefined file containing all definitions that can be used by the command:

#INCLUDE atoms.kpp

This should be the first line in a KPP input file, because it allows to use any atom in the periodic table of elements throughout the kinetic description file.

#CHECK

KPP is able to do mass balance checks for all equations. Some chemical equations are not balanced for all atoms, and this might still be correct from a chemical point of view. To accommodate for this, KPP can perform mass balance checking only for the list of atoms specified in the #CHECK section, e.g.:

#CHECK N; C; O;

The balance checking for all atoms can be enabled by using the #CHECKALL command. Without #CHECK or #CHECKALL, no checking is performed. The IGNORE atom can also be used to control mass balance checking.

#DEFVAR and #DEFFIX

There are two ways to declare new species together with their atom composition: #DEFVAR and #DEFFIX. These sections define all the species that will be used in the chemical mechanism. Species can be variable or fixed. The type is implicitly specified by defining the species in the appropriate sections. A fixed species does not vary through chemical reactions.

For each species the user has to declare the atom composition. This information is used for mass balance checking. To ignore mass balance checking for a given species, one can declare the predefined atom IGNORE as being part of the species composition. Examples for these sections are:

#DEFVAR
  NO2 = N + 2O;
  CH3OOH = C + 4H + 2O;
  HSO4m = IGNORE;
  RCHO = IGNORE;
#DEFFIX
  CO2 = C + 2O;

#EQUATIONS

The chemical mechanism is specified in the #EQUATIONS section. Each equation is written in the natural way in which a chemist would write it:

#EQUATIONS

<R1> NO2 + hv = NO + O3P :  6.69e-1*(SUN/60.0);
<R2> O3P + O2 + AIR = O3 :  ARR_ac(5.68e-34,  -2.80);
<R3> O3P + O3 = 2O2 :       ARR_ab(8.00e-12, 2060.0);
<R4> O3P + NO + AIR = NO2 : ARR_ac(1.00e-31,  -1.60);
//... etc ...

Note

The above example is taken from the saprc99 mechanism (see models/saprc99.eqn), with some whitespace deleted for clarity. Optional equation tags are specified by text within < > angle brackets. Functions that compute saprc99 equation rates (e.g. ARR_ac, ARR_ab) are defined in util/UserRateLaws.f90 and util/UserRateLawsInterfaces.f90.

Only the names of already defined species can be used. The rate coefficient has to be placed at the end of each equation, separated by a colon. The rate coefficient does not necessarily need to be a numerical value. Instead, it can be a valid expression (or a call to an inlined rate law function) in the target language. If there are several #EQUATIONS sections in the input, their contents will be concatenated.

A minus sign in an equation shows that a species is consumed in a reaction but it does not affect the reaction rate. For example, the oxidation of methane can be written as:

CH4 + OH = CH3OO + H2O - O2 : k_CH4_OH;

However, it should be noted that using negative products may lead to numerical instabilities.

Often, the stoichiometric factors are integers. However, it is also possible to have non-integer yields, which is very useful to parameterize organic reactions that branch into several side reactions:

CH4 + O1D = .75 CH3O2 + .75 OH + .25 HCHO + 0.4 H + .05 H2 : k_CH4_O1D;

KPP provides two pre-defined dummy species: hv and PROD. Using dummy species does not affect the numerics of the integrators. It only serves to improve the readability of the equations. For photolysis reactions, hv can be specified as one of the reagents to indicate that light ($h\nu$) is needed for this reaction, e.g.:

NO2 + hv = NO + O : J_NO2;

When the products of a reaction are not known or not important, the dummy species PROD should be used as a product. This is necessary because the KPP syntax does not allow an empty list of products. For example, the dry deposition of atmospheric ozone to the surface can be written as:

O3 = PROD : v_d_O3;

The same equation must not occur twice in the #EQUATIONS section. For example, you may have both the gas-phase reaction of N2O5 with water in your mechanism and also the heterogeneous reaction on aerosols:

N2O5 + H2O = 2 HNO3 : k_gas;
N2O5 + H2O = 2 HNO3 : k_aerosol;

These reactions must be merged by adding the rate coefficients:

N2O5 + H2O = 2 HNO3 : k_gas + k_aerosol;

#FAMILIES

Chemical families (for diagnostic purposes) may be specified in the #FAMILIES section as shown below. Family names beginning with a P denote production, and those beginning with an L denote loss.

#FAMILIES
  POx : O3 + NO2 + 2NO3 + HNO3 + ... etc. add more species as needed ...
  LOx : O3 + NO2 + 2NO3 + HNO3 + ... etc. add more species as needed ...
  PCO : CO;
  LCO : CO;
  PSO4 : SO4;
  LCH4 : CH4;
  PH2O2 : H2O2;

KPP will examine the chemical mechanism and create a dummy species for each defined family. Each dummy species will archive the production and loss for the family. For example, each molecule of CO that is produced will be added to the PCO dummy species. Likewise, each molecule of CO that is consumed will be added to the LCO dummy species. This will allow the PCO and LCO species to be later archived for diagnostic purposes. Dummy species for chemical families will not be included as active species in the mechanism.

#INITVALUES

The initial concentration values for all species can be defined in the #INITVALUES section, e.g.:

#INITVALUES
  CFACTOR = 2.5E19;
  NO2 = 1.4E-9;
  CO2 = MyCO2Func();
  ALL_SPEC = 0.0;

If no value is specified for a particular species, the default value zero is used. One can set the default values using the generic species names: VAR_SPEC, FIX_SPEC, and ALL_SPEC. In order to use coherent units for concentration and rate coefficients, it is sometimes necessary to multiply each value by a constant factor. This factor can be set by using the generic name CFACTOR. Each of the initial values will be multiplied by this factor before being used. If CFACTOR is omitted, it defaults to one.

The information gathered in this section is used to generate the Initialize subroutine (cf ROOT_Initialize). In more complex 3D models, the initial values are usually taken from some input files or some global data structures. In this case, #INITVALUES may not be needed.

#LOOKAT and #MONITOR

There are two sections in this category: #LOOKAT and #MONITOR.

The section instructs the preprocessor what are the species for which the evolution of the concentration, should be saved in a data file. By default, if no #LOOKAT section is present, all the species are saved. If an atom is specified in the #LOOKAT list then the total mass of the particular atom is reported. This allows to check how the mass of a specific atom was conserved by the integration method. The #LOOKATALL command can be used to specify all the species. Output of #LOOKAT can be directed to the file ROOT.dat using the utility subroutines described in the section entitled ROOT_Util.

The #MONITOR section defines a different list of species and atoms. This list is used by the driver to display the concentration of the elements in the list during the integration. This may give us a feedback of the evolution in time of the selected species during the integration. The syntax is similar to the #LOOKAT section. With the driver general, output of #MONITOR goes to the screen (STDOUT). The order of the output is: first variable species, then fixed species, finally atoms. It is not the order in the MONITOR command.

Examples for these sections are:

#LOOKAT NO2; CO2; O3; N;
#MONITOR O3; N;

#SETVAR and #SETFIX

The commands #SETVAR and #SETFIX change the type of an already defined species. Then, depending on the integration method, one may or may not use the initial classification, or can easily move one species from one category to another. The use of the generic species VAR_SPEC, FIX_SPEC, and ALL_SPEC is also allowed. Examples for these sections are:

#SETVAR ALL_SPEC;
#SETFIX H2O; CO2;

KPP commands

A KPP command begins on a new line with a # sign, followed by a command name and one or more parameters. Details about each command are given in the following subsections.

Default values for KPP commands
KPP command	default value
#AUTOREDUCE	`OFF`
#CHECKALL
#DECLARE	`SYMBOL`
#DOUBLE	`ON`
#DRIVER	`none`
#DUMMYINDEX	`OFF`
#EQNTAGS	`OFF`
#FUNCTION	`AGGREGATE`
#HESSIAN	`ON`
#INCLUDE
#INTEGRATOR
#INTFILE
#JACOBIAN	`SPARSE_LU_ROW`
#LANGUAGE
#LOOKATALL
#MEX	`ON`
#MINVERSION
#MODEL
#REORDER	`ON`
#STOCHASTIC	`OFF`
#STOICMAT	`ON`
#UPPERCASEF90	`OFF`

#AUTOREDUCE

The #AUTOREDUCE ON command can be used with #INTEGRATOR rosenbrock to enable automatic mechanism reduction as described in Lin et al. [2023]. Automatic mechanism reduction is disabled by default.

#DECLARE

The #DECLARE command determines how constants like dp, NSPEC, NVAR, NFIX, and NREACT are inserted into the KPP-generated code. #DECLARE SYMBOL (the default) will declare array variables using parameters from the ROOT_Parameters file. #DECLARE VALUE will replace each parameter with its value.

For example, the global array variable C is declared in the ROOT_Global file generated by KPP. In the small_strato example (described in Running KPP with an example stratospheric mechanism), C has dimension NSPEC=7. Using #DECLARE SYMBOL will generate the following code in ROOT_Global:

! C - Concentration of all species
  REAL(kind=dp), TARGET :: C(NSPEC)

Whereas #DECLARE VALUE will generate this code instead:

! C - Concentration of all species
  REAL(kind=dp), TARGET :: C(7)

We recommend using #DECLARE SYMBOL, as most modern compilers will automatically replace each parameter (e.g. NSPEC) with its value (e.g 7). However, if you are using a very old compiler that is not as sophisticated, #DECLARE VALUE might result in better-optmized code.

#DOUBLE

The #DOUBLE command selects single or double precision arithmetic. ON (the default) means use double precision, OFF means use single precision (see the section entitled ROOT_Precision).

Important

We recommend using double precision whenever possible. Using single precision may lead to integration non-convergence errors caused by roundoff and/or underflow.

#DRIVER

The #DRIVER command selects the driver, i.e., the file from which the main function is to be taken. The parameter is a file name, without suffix. The appropriate suffix (.f90, .F90, .c, or .m) is automatically appended.

Normally, KPP tries to find the selected driver file in the directory $KPP_HOME/drv/. However, if the supplied file name contains a slash, it is assumed to be absolute. To access a driver in the current directory, the prefix ./ can be used, e.g.:

#DRIVER ./mydriver

It is possible to choose the empty dummy driver none, if the user wants to include the KPP generated modules into a larger model (e.g. a general circulation or a chemical transport model) instead of creating a stand-alone version of the chemical integrator. The driver none is also selected when the #DRIVER command is missing. If the command occurs twice, the second replaces the first.

#DUMMYINDEX

It is possible to declare species in the #DEFVAR and #DEFFIX sections that are not used in the #EQUATIONS section. If your model needs to check at run-time if a certain species is included in the current mechanism, you can set to #DUMMYINDEX ON. Then, KPP will set the indices to zero for all species that do not occur in any reaction. With #DUMMYINDEX OFF (the default), those are undefined variables. For example, if you frequently switch between mechanisms with and without sulfuric acid, you can use this code:

IF (ind_H2SO4=0) THEN
  PRINT *, 'no H2SO4 in current mechanism'
ELSE
  PRINT *, 'c(H2SO4) =', C(ind_H2SO4)
ENDIF

#EQNTAGS

Each reaction in the #EQUATIONS section may start with an equation tag which is enclosed in angle brackets, e.g.:

<R1> NO2 + hv = NO + O3P :  6.69e-1*(SUN/60.0);

With #EQNTAGS set to ON, this equation tag can be used to refer to a specific equation (cf. ROOT_Monitor). The default for #EQNTAGS is OFF.

#FUNCTION

The #FUNCTION command controls which functions are generated to compute the production/destruction terms for variable species. AGGREGATE generates one function that computes the normal derivatives. SPLIT generates two functions for the derivatives in production and destruction forms.

#HESSIAN

The option ON (the default) of the #HESSIAN command turns the Hessian generation on (see section ROOT_Hessian and ROOT_HessianSP). With OFF it is switched off.

#INCLUDE

The #INCLUDE command instructs KPP to look for the file specified as a parameter and parse the content of this file before proceeding to the next line. This allows the atoms definition, the species definition and the equation definition to be shared between several models. Moreover this allows for custom configuration of KPP to accommodate various classes of users. Include files can be either in one of the KPP directories or in the current directory.

#INTEGRATOR

The #INTEGRATOR command selects the integrator definition file. The parameter is the file name of an integrator, without suffix. The effect of

#INTEGRATOR integrator_name

is similar to:

#INCLUDE $KPP_HOME/int/integrator_name.def

The #INTEGRATOR command allows the use of different integration techniques on the same model. If it occurs twice, the second replaces the first. Normally, KPP tries to find the selected integrator files in the directory $KPP_HOME/int/. However, if the supplied file name contains a slash, it is assumed to be absolute. To access an integrator in the current directory, the prefix ./ can be used, e.g.:

#INTEGRATOR ./mydeffile

#INTFILE

Attention

#INTFILE is used internally by KPP but should not be used by the KPP user. Using #INTEGRATOR alone suffices to specify an integrator.

The integrator definition file selects an integrator file with #INTFILE and also defines some suitable options for it. The #INTFILE command selects the file that contains the integrator routine. The parameter of the command is a file name, without suffix. The appropriate suffix (.f90, .F90, .c, or .m is appended and the result selects the file from which the integrator is taken. This file will be copied into the code file in the appropriate place.

#JACOBIAN

The #JACOBIAN command controls which functions are generated to compute the Jacobian. The option OFF inhibits the generation of the Jacobian routine. The option FULL generates the Jacobian as a square NVAR x NVAR matrix. It should only be used if the integrator needs the whole Jacobians. The options SPARSE_ROW and SPARSE_LU_ROW (the default) both generate the Jacobian in sparse (compressed on rows) format. They should be used if the integrator needs the whole Jacobian, but in a sparse form. The format used is compressed on rows. With SPARSE_LU_ROW, KPP extends the number of nonzeros to account for the fill-in due to the LU decomposition.

#LANGUAGE

Attention

The Fortran77 language option is deprecated in KPP 2.5.0 and later versions. All further KPP development will only support Fortran90.

The #LANGUAGE command selects the target language in which the code file is to be generated. Available options are Fortran90, C, or matlab.

You can select the suffix (.F90 or .f90) to use for Fortran90 source code generated by KPP (cf. #UPPERCASEF90).

#MEX

Mex is a Matlab extension that allows to call functions written in Fortran and C directly from within the Matlab environment. KPP generates the mex interface routines for the ODE function, Jacobian, and Hessian, for the target languages C, Fortran77, and Fortran90. The default is #MEX ON. With #MEX OFF, no Mex files are generated.

#MINVERSION

You may restrict a chemical mechanism to use a given version of KPP or later. To do this, add

#MINVERSION X.Y.Z

to the definition file.

The version number (X.Y.Z) adheres to the Semantic Versioning style (https://semver.org), where X is the major version number, Y is the minor version number, and Z is the bugfix (aka “patch”) version number.

For example, if #MINVERSION 2.4.0 is specified, then KPP will quit with an error message unless you are using KPP 2.4.0 or later.

#MODEL

The chemical model contains the description of the atoms, species, and chemical equations. It also contains default initial values for the species and default options including a suitable integrator for the model. In the simplest case, the main kinetic description file, i.e. the one passed as parameter to KPP, can contain just a single line selecting the model. KPP tries to find a file with the name of the model and the suffix .def in the $KPP_HOME/models subdirectory. This file is then parsed. The content of the model definition file is written in the KPP language. The model definition file points to a species file and an equation file. The species file includes further the atom definition file. All default values regarding the model are automatically selected. For convenience, the best integrator and driver for the given model are also automatically selected.

The #MODEL command is optional, and intended for using a predefined model. Users who supply their own reaction mechanism do not need it.

#REORDER

Reordering of the species is performed in order to minimize the fill-in during the LU factorization, and therefore preserve the sparsity structure and increase efficiency. The reordering is done using a diagonal Markowitz algorithm. The details are explained in Sandu et al. [1996]. The default is ON. OFF means that KPP does not reorder the species. The order of the variables is the order in which the species are declared in the #DEFVAR section.

#STOCHASTIC

The option ON of the #STOCHASTIC command turns on the generation of code for stochastic kinetic simulations (see the section entitled ROOT_Stochastic. The default option is OFF.

#STOICMAT

Unless the #STOICMAT command is set to OFF, KPP generates code for the stoichiometric matrix, the vector of reactant products in each reaction, and the partial derivative of the time derivative function with respect to rate coefficients (cf. ROOT_Stoichiom and ROOT_StoichiomSP).

#CHECKALL, #LOOKATALL

The shorthand commands #CHECKALL and #LOOKATALL apply #CHECK and #LOOKAT, respectively, to all species in the mechanism.

#UPPERCASEF90

If you have selected #LANGUAGE Fortran90 option, KPP will generate source code ending in .f90 by default. Setting #UPPERCASEF90 ON will tell KPP to generate Fortran90 code ending in .F90 instead.

Inlined Code

In order to offer maximum flexibility, KPP allows the user to include pieces of code in the kinetic description file. Inlined code begins on a new line with #INLINE and the inline_type. Next, one or more lines of code follow, written in the target language (Fortran90, C, or Matlab) as specified by the inline_type. The inlined code ends with #ENDINLINE. The code is inserted into the KPP output at a position which is also determined by inline_type as shown in KPP inlined types. If two inline commands with the same inline type are declared, then the contents of the second is appended to the first one.

List of inlined types

In this manual, we show the inline types for Fortran90. The inline types for the other languages are produced by replacing F90 by C, or matlab, respectively.

KPP inlined types
Inline_type	File	Placement	Usage
F90_DATA	ROOT_Monitor	specification section	(obsolete)
F90_GLOBAL	ROOT_Global	specification section	global variables
F90_INIT	ROOT_Initialize	subroutine	integration parameters
F90_RATES	ROOT_Rates	executable section	rate law functions
F90_RCONST	ROOT_Rates	subroutine	statements and definitions of rate coefficients
F90_UTIL	ROOT_Util	executable section	utility functions

F90_DATA

This inline type was introduced in a previous version of KPP to initialize variables. It is now obsolete but kept for compatibility. For Fortran90, F90_GLOBAL should be used instead.

F90_GLOBAL

This inline type can be used to declare global variables, e.g. for a special rate coefficient:

#INLINE F90_GLOBAL
  REAL(dp) :: k_DMS_OH
#ENDINLINE

Inlining code can be useful to introduce additional state variables (such as temperature, humidity, etc.) for use by KPP routines, such as for calculating rate coefficients.

If a large number of state variables needs to be held in inline code, or require intermediate computation that may be repeated for many rate coefficients, a derived type object should be used for efficiency, e.g.:

#INLINE F90_GLOBAL
  TYPE, PUBLIC :: ObjGlobal_t
     ! ... add variable fields to this type ...
  END TYPE ObjGlobal_t
  TYPE(ObjGlobal_t), TARGET, PUBLIC :: ObjGlobal
#ENDINLINE

This global variable ObjGlobal can then be used globally in KPP.

Another way to avoid cluttering up the KPP input file is to #include a header file with global variables:

#INLINE F90_GLOBAL
! Inline common variables into KPP_ROOT_Global.f90
#include "commonIncludeVars.f90"
#ENDINLINE

In future versions of KPP, the global state will be reorganized into derived type objects as well.

F90_INIT

This inline type can be used to define initial values before the start of the integration, e.g.:

#INLINE F90_INIT
  TSTART = (12.*3600.)
  TEND = TSTART + (3.*24.*3600.)
  DT = 0.25*3600.
  TEMP = 270.
#ENDINLINE

F90_RATES

This inline type can be used to add new subroutines to calculate rate coefficients, e.g.:

#INLINE F90_RATES
  REAL FUNCTION k_SIV_H2O2(k_298,tdep,cHp,temp)
    ! special rate function for S(IV) + H2O2
    REAL, INTENT(IN) :: k_298, tdep, cHp, temp
    k_SIV_H2O2 = k_298 &
      * EXP(tdep*(1./temp-3.3540E-3)) &
      * cHp / (cHp+0.1)
  END FUNCTION k_SIV_H2O2
#ENDINLINE

F90_RCONST

This inline type can be used to define time-dependent values of rate coefficients. You may inline USE statements that reference modules where rate coefficients are computed, e.g.:

#INLINE F90_RCONST
  USE MyRateFunctionModule
#ENDINLINE

or define variables directly, e.g.:

#INLINE F90_RCONST
  k_DMS_OH = 1.E-9*EXP(5820./temp)*C(ind_O2)/ &
    (1.E30+5.*EXP(6280./temp)*C(ind_O2))
#ENDINLINE

Note that the USE statements must precede any variable definitions.

The inlined code will be placed directly into the subroutines UPDATE_RCONST and UPDATE_PHOTO in the ROOT_Rates file.

F90_UTIL

This inline type can be used to define utility subroutines.

Auxiliary files and the substitution preprocessor

The auxiliary files in the $KPP_HOME/util subdirectory are templates for integrators, drivers, and utilities. They are inserted into the KPP output after being run through the substitution preprocessor. This preprocessor replaces several placeholder symbols in the template files with their particular values in the model at hand. Usually, only KPP_ROOT and KPP_REAL are needed because the other values can also be obtained via the variables listed in KPP inlined types.

KPP_REAL is replaced by the appropriate single or double precision declaration type. Depending on the target language KPP will select the correct declaration type. For example if one needs to declare an array BIG of size 1000, a declaration like the following must be used:

KPP_REAL :: BIG(1000)

When used with the command #DOUBLE ON, the above line will be automatically translated into:

REAL(kind=dp) :: BIG(1000)

and when used with the command #DOUBLE OFF, the same line will become:

REAL(kind=sp) :: BIG(1000)

in the resulting Fortran90 output file.

KPP_ROOT is replaced by the root file name of the main kinetic description file. In our example where we are processing small_strato.kpp, a line in an auxiliary Fortran90 file like

USE KPP_ROOT_Monitor

will be translated into

USE small_strato_Monitor

in the generated Fortran90 output file.

List of auxiliary files for Fortran90

Auxiliary files for Fortran90
File	Contents
`dFun_dRcoeff.f90`	Derivatives with respect to reaction rates.
`dJac_dRcoeff.f90`	Derivatives with respect to reaction rates.
`Makefile_f90` and `Makefile_upper_F90`	Makefiles to build Fortran-90 code.
`Mex_Fun.f90`	Mex files.
`Mex_Jac_SP.f90`	Mex files.
`Mex_Hessian.f90`	Mex files.
`sutil.f90`	Sparse utility functions.
`tag2num.f90`	Function related to equation tags.
`UpdateSun.f90`	Function related to solar zenith angle.
`UserRateLaws.f90` and `UserRateLawsInterfaces.f90`	User-defined rate-law functions.
`util.f90`	Input/output utilities.

List of symbols replaced by the substitution preprocessor

Symbols and their replacements
Symbol	Replacement	Example
KPP_ROOT	The `ROOT` name	`small_strato`
KPP_REAL	The real data type	`REAL(kind=dp)`
KPP_NSPEC	Number of species	7
KPP_NVAR	Number of variable species	5
KPP_NFIX	Number of fixed species	2
KPP_NREACT	Number of chemical reactions	10
KPP_NONZERO	Number of Jacobian nonzero elements	18
KPP_LU_NONZERO	Number of Jacobian nonzero elements, with LU fill-in	19
KPP_LU_NHESS	Number of Hessian nonzero elements	10
KPP_FUN_OR_FUN_SPLIT	Name of the function to be called	`FUN(Y,FIX,RCONST,Ydot)`

Controlling the Integrator with `ICNTRL` and `RCNTRL`

In order to offer more control over the integrator, KPP provides the arrays ICNTRL (integer) and RCNTRL (real). Each of them is an array of 20 elements that allow the fine-tuning of the integrator. All integrators (except for tau_leap and gillespie) use ICNTRL and RCNTRL. Details can be found in the comment lines of the individual integrator files in $KPP_HOME/int/.

ICNTRL

Summary of ICNTRL usage in the f90 integrators. Here, Y = used, and s = solver-specific usage.
ICNTRL	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17
beuler		Y	Y	Y	Y	Y									Y
dvode															Y
exponential
feuler															Y	Y	Y
gillespie
lsode		Y		Y						s					Y
radau5		Y		Y	Y	Y					Y				Y
rosenbrock_adj	Y	Y	Y	Y		s	s	s							Y
rosenbrock	Y	Y	Y	Y											Y	Y
rosenbrock_tlm	Y	Y	Y	Y								s			Y
rosenbrock_autoreduce	Y	Y	Y	Y								s	s	s	Y	Y
runge_kutta_adj		Y	Y	Y	Y	s	s	s	s	s	Y				Y
runge_kutta		Y	Y	Y	Y	Y				s	Y				Y
runge_kutta_tlm		Y	Y		Y	Y	s		s	s	Y	s			Y
sdirk4		Y		Y											Y
sdirk_adj		Y	Y	Y	Y	Y	s	s							Y
sdirk		Y	Y	Y	Y	Y									Y
sdirk_tlm		Y	Y	Y	Y	Y	s		s			s			Y
seulex	Y	Y		Y						s	s	s	s	s	Y
tau_leap

ICNTRL(1)

= 1: $F = F(y)$, i.e. independent of t (autonomous)

= 0: $F = F(t,y)$, i.e. depends on t (non-autonomous)

ICNTRL(2)

The absolute (ATOL) and relative (RTOL) tolerances can be expressed by either a scalar or individually for each species in a vector:

= 0 : NVAR -dimensional vector

= 1 : scalar

ICNTRL(3): Selection of a specific method.

ICNTRL(4): Maximum number of integration steps.

ICNTRL(5): Maximum number of Newton iterations.

ICNTRL(6)

Starting values of Newton iterations (only avaialble for some of the integrators).

= 0 : Interpolated

= 1 : Zero

ICNTRL(11): Gustafsson step size controller

ICNTRL(12): (Solver-specific for rosenbrock_autoreduce) Controls whether auto-reduction of the mechanism is performed. If set to = 0, then the integrator behaves the same as rosenbrock.

ICNTRL(13): (Solver-specific for rosenbrock_autoreduce) Controls whether in auto-reduction species production and loss rates are scanned throughout the internal time steps of the integrator for repartitioning.

ICNTRL(14): (Solver-specific for rosenbrock_autoreduce) If set to > 0, then the threshold is calculated based on the max of production and loss rate of the species ID specified in ICNTRL(14) multiplied by RCNTRL(14).

ICNTRL(15)

This determines which Update_* subroutines are called within the integrator.

= -1 : Do not call any Update_* subroutines

= 0 : Use the integrator-specific default values

> 1 : A number between 1 and 7, derived by adding up bits with values 4, 2, and 1. The first digit (4) activates Update_SUN. The second digit (2) activates Update_PHOTO. The third digit (1) activates Update_RCONST. |

For example ICNTRL(15)=6) (4+2) will activate the calls to Update_SUN and Update_PHOTO, but not to Update_RCONST.

ICNTRL(16)

Treatment of negative concentrations:

= 0 : Leave negative values unchanged

= 1 : Set negative values to zero

= 2 : Print warning and continue

= 3 : Print error message and stop

ICNTRL(17)

Verbosity:

= 0 : Only return error number

= 1 : Verbose error output

ICNTRL(18) ... ICNTRL(20): currently not used

RCNTRL

Summary of RCNTRL usage in the f90 integrators. Here, Y = used, and s = solver-specific usage.
RCNTRL	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19
beuler	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
dvode
exponential
feuler
gillespie
lsode	Y	Y	Y
radau5		Y		Y	Y	Y	Y	Y	Y	Y	Y
rosenbrock_adj	Y	Y	Y	Y	Y	Y	Y
rosenbrock	Y	Y	Y	Y	Y	Y	Y
rosenbrock_tlm	Y	Y	Y	Y	Y	Y	Y
rosenbrock_autoreduce	Y	Y	Y	Y	Y	Y	Y					s		s
runge_kutta_adj	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
runge_kutta	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
runge_kutta_tlm	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
sdirk4	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
sdirk_adj	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
sdirk	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
sdirk_tlm	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y
seulex	Y	Y	Y	Y	Y	Y	Y	Y		s	s	s	s	s	s	s	s	s	s
tau_leap

RCNTRL(1): Hmin, the lower bound of the integration step size. It is not recommended to change the default value of zero.

RCNTRL(2): Hmax, the upper bound of the integration step size.

RCNTRL(3): Hstart, the starting value of the integration step size.

RCNTRL(4): FacMin, lower bound on step decrease factor.

RCNTRL(5): FacMax, upper bound on step increase factor.

RCNTRL(6): FacRej, step decrease factor after multiple rejections.

RCNTRL(7): FacSafe, the factor by which the new step is slightly smaller than the predicted value.

RCNTRL(8): ThetaMin. If the Newton convergence rate is smaller than ThetaMin, the Jacobian is not recomputed.

RCNTRL(9): NewtonTol, the stopping criterion for Newton’s method.

RCNTRL(10): Qmin

RCNTRL(11): Qmax. If Qmin < Hnew/Hold < Qmax, then the step size is kept constant and the LU factorization is reused.

RCNTRL(12): (Solver-specific for rosenbrock_autoreduce) Used to specify the threshold for auto-reduction partitioning, if ICNTRL(12) = 1, and ICNTRL(14) = 0. Will be ignored if ICNTRL(14) > 0.

RCNTRL(14): (Solver-specific for rosenbrock_autoreduce) Used to specify the multiplier for threshold for auto-reduction partitioning, if ICNTRL(12) = 1, and ICNTRL(14) > 0, RCNTRL(14) is multiplied against max of production and loss rates of species ICNTRL(14) to produce the partitioning threshold, ignoring RCNTRL(12).

RCNTRL(10) ... RCNTRL(19): (Solver-specific for seulex)

RCNTRL(20): currently not used