• chem

The CLaMS chemistry module chem

1.1 The program chem

The program chem is a box-trajectory model that is capable of running on a multiple trajectories at the same time. It calculates the change of chemical composition along each given trajectory. The integration is achieved using the ASAD package provided by Glenn Carver at Cambridge. The photolysis code is basically the David Lary code with a few additions. The heterogeneous chemistry is based on the code of Ken Carslaw and used analytical expressions to determine the condensation of sulphate, nitrate and water vapor. The ASAD routines are written in Fortran 77 whereas all other routines are now updated to Fortran 90. The program can run on many systems like AIX and Linux systems including JUMP.

1.2 How to get the source code

The source code of chem in Jülich is under the control of the CVS (Concurrent Version System) utility. To get the source code first make sure that the environment variable CVSROOT is set to /usr/local/icg/icg1/archive:

cvs co chem

This command creates the subdirectories chem and its subdirectories. (Alternatively, the command cvs co clams creates all CLaMS programs. The chemistry program is then located in the directory clams/chem). The files are placed into the following directory structure:

2 Compiling the model

Using the Makefile in the directory chem the program chem can be compiled using the Makefile:

make chem

or for the MPI version for parallel computing:

make useMPI=true chem

The actual version contains the het. chemistry parametrisation of Shi et al. (2001).

The makefile contains all dependencies between the different source and object files. Especially after a change of dimensions (i.e. in file paramcdef.h) the program needs to be recompiled. There are a lot of options in which this Makefile can be used. Type make for a list of the options.

3 Running the model

3.1 Input and control parameters

To run the program chem an input control file chem.inp must be provided. Typically it is created in a script file scripts/*.exec that afterwards starts the program. It contains model control parameters as the filenames of trajectory files, switches for the different use of the model as explained below. The informations are provided via a namelist, i.e. parameters that are left out will be set to default values. The different namelists in the appearing order in chem.inp are summarized in the overleaf table and explained below. The order of variables within one namelist is irrelevant. The trajectory data contained in file name dsn in (namelist CNTL) This data which should conform to the standard NetCDF trajectory format and should include at least pressure temperature data reside in the directory testData for the given example scripts. Any other directory may be chosen. The solar zenith angle (SZA) will be calculated by the model at each chemistry timestep (if chemical substeps are chosen, nchemsteps > 1, at each substep) The meaning of the namelists and its parameters that need to be in the file chem.inp is summarized in the table below.

3.2 Initialization

There are two possibilities for providing the initial mixing ratios of the chemical species:

1. Initialization via namelist SPECIES in chem.inp: The number of species to be initialized with nonzero value needs to be specified by the variable mspecies. This needs to be followed by this number of namelists containing species name and initial mixing ratio in arbitrary order. In the case of multiple trajectories, each individual trajectory is then initialized with the same values.

2. Initialization via NetCDF input file (default initial.nc). This method is chosen if the variable mspecies in the namelist SPECIES is set to zero. The filename of the initialization file can be set by the variable file init in the namelin the namelist SPECIES.

4 ASAD Parameters

4.1 Chemical Species

The names chemical species (and potentially chemical families) are provided in the array chch_t. defined in the source code under asad/clams/clams_chem_data_<type>.f90. The dimensions (number of chemical species and number reactions) in these files is set in clams_chem_data.f90 for the chemistry reaction scheme <type> (see below).

The second column of this file contains the chemical species name after the ASAD naming conventions (see ASAD manual for details). The fourth column contains the kind of chemical species (TR: chemical tracer, FM: chemical family member, SS: species in steady state, CT: constant mixing ratio). For family members, the fifth column contains the corresponding family name and the third column is the number of times, how often a species molecule counts within the family.

4.2 Chemical Reactions

The chemical reactions and their reaction rates are set in the arrays ratb_t (bimolecular reactions), ratt_t (termolecular reactions), ratj_t (photolysis reactions) and rath_t (heterogeneous reactions) defined in the source code under asad/clams/clams_chem_data_<type>.f90.

• The dimensions (number of chemical species and number reactions) in these

files is set in clams_chem_data.f90 for the chemistry reaction scheme <type> (see below).

5 Chemistry Schemes

There are the following chemistry schemes implemented in chem. They can be chosen through the variable chemdata_type in the namelisnt CNTL

6 Solver Algorithms

Within the ASAD package the user can choose between different algorithms that solve the underlying differential equations. The choice of the solver is controlled by the variable method in the namelist CNTL.

IMPACT (method=1): The solver IMPACT is based on the Family concept. Several species must be combined into families, preferably in a way that chemical changes within the family may be fast, whereas the mixing ratio of the total family changes relatively slowly, e.g.

[Ox] := [O3] + [O(1D)] + [O(3P)]

The definition of the families is given in the file chch.d, in which the entry ’FM’ in the third column tells that a species is a family member of the family denoted in column 4.

NEWTON RAPHSON (method=3) Integration using the Newtaon Raphson integrator. In case of no convergence with the timestep dt half timestep ist used (iterative, max_redo=16 compared to dt/128) from M. Prather and O. Wild. Some comparison between IMPACT (method=1) and NEWTON RAPHSON (method=3) have been performed, see here.

NAG (method=10) and SVODE (method=11): The solver name NAG stands for a routine of the NAG library (d02eaf) that solves a stiff system of ordinary differential equations. SVODE is a similar solver provided with the ASAD code. The advantage of a ”stiff solver” is that no approximations about families have to be made. However these solvers are numerically costly. There is no principal difference between NAG and SVODE, besides that SVODE is faster. Depending if family solver or a stiff solver is used, different files for the definition of chemical species/families (data/chch.d) and model dimensions (include/paramcdef.h) are needed. Examples for standard (family) and stiff files are located in data/datatypes/chch.* and include/incltypes/paramcdef.*. After a change of paramcdef.h, the model needs to be re-compiled. This procedure can be done in one step using the makefile command make stiff or make standard. These commands use also a program called check dims that checks the dimensions of the file paramcdef.h and corrects it automatically.

7 Timesteps

According to the flexible use of the chemistry model, there are different timesteps:

1. The dynamical timestep: It corresponds to the frequency at which dynamic data are available in the trajectory file.

2. The chemical input timestep: This corresponds to the frequency at which the chemistry routine is called. The variable nchemsteps tells, how many chemistry calls are conducted per dynamical timestep. This is needed, if the dynamical timestep catches not the fast variations of the solar zenith angle.

Be careful when using nchemsteps > 1 together with a trajectory with unequal timesteps, e.g. a 7h long trajectory with a 2h timestep and nchemsteps=2. This would result in two final steps of 0.5 h

3. The regular internal chemistry timestep at which the solver is called: ncdt (seconds). It is given in the namelist CNTL.

4. The output timestep: This is the timestep at which data are stored in the output file. The variable iout_skip controls that: The results are written out at every iout_skip chemistry time step.

5. The timestep sedidt is given for writout into the files to be read by the program sedi. It should be divisable by the chemical input timestep (2.).

8 Heterogeneous Chemistry

The aerosol particle microphysics and the heterogeneous reaction rates are calculated by version 6 of the code developed by Ken Carslaw and updated using the parametrisation by Shi et al. (2001). Four different phases of the stratospheric aerosols are possible:

1. liquid aerosol as ternary H2SO4/HNO3/H2O solution

2. solid sulfuric acid tetrahydrate (SAT)

3. solid nitric acid trihydrate (NAT)

4. solid water ice (ice)

Also combinations of the 4 phases may exist. All input parameters that controls the particle microphysics heterogeneous reaction rates need to be put into the file hetinit6.dat. It contains the variables in the following order:

For the reactions that can occur on more than one substrate, the program chem uses the sum of the reaction rates on each substrate. The number density of the single aerosol types (N(aero), N(SAT), N(NAT), N(ice)) and the amount of sulfuric acid in the aerosol aer H2SO4 for each trajectory are written to the output file. If N(aero) and aer H2SO4 are present and positive in the initialization file, these values are used instead of N0 and ppbh2so4 in the file hetinit6.dat. By default the initial particle phase is liquid. However, if N(SAT), N(NAT), or N(ice) are present in the initialization file, it will be accounted for.

9 Photolysis

9.1 Lary-Photolysis

The Lary photolysis code solves the radiative transfer equation on the basis of a spherical geometry scheme. It is chosen if in the namelist PHOTO the variable iphot is set to 1 (default). The control parameters are read in via the namelist PHOTO. The albedo needs to be specified. The ozone profile (ozdsn) that is needed for the calculation can either be a single ozone profile provided as ASCII file or multiple ozone profiles given as NetCDF-file. In the latter case ozone profiles for different latitudes can be used that is important for global, multiple trajectory calculations. Also a set of ozone profiles for different times is possible. Here, the integer variable ntime o3 controls which timestep is chosen. The form of the NetCDF file needs to be as the output file of the Mainz 2-D model. If the Mainz 2-D model output file is used here, ntime o3 corresponds to the number of month.

Since the initialization of this calculation is numerically costly, there is the possibility to saving the parameters (e.g. table of enhancement factors as a function of wavelength altitude, SZA, and latitude) for the next similar calculation. To do this, the logical variable quickread must be set to true. Then the table is written to dissoc.save at the first run and read in from this file at the following runs.

The parameters for the Lary-photolysis code that need to be defined prior to compilation are set in the module photodim M (in file source/pmodules.f90): The number of model levels is given through jpslevall. The variable ptop defines the top pressure level in the photolysis calculation. The model uses either the top of the given ozone profile or ptop, whatever is the lower pressure level. The parameter maxlats contains the maximum number of different latitudes for which the photolysis table is calculated. maxlats=18 is needed for global 2-D model ozone data, whereas one needs only maxlats=1 for a box model or a calculation in a narrow latitude domain. Be aware that large values of maxlats cause the need of high runtime memory.

9.2 Diurnal average photolysis rates

If iphot=0 is set, the model calculates latitude and altitude dependent diurnally averaged photolysis rates. This can be used for longterm integrations in which the diurnal cycle of chemical species is not important.

9.3 Röth-Approximation Formula

The approximation formula by E.P. Röth is chosen if in the namelist PHOTO the variable iphot is set to 2. In that case the remaining information in the namelist PHOTO is irrelevant and all information is read in from the file roeth.dat. It must contain the 3 parameters A, B and C of each photolysis reaction for the conditions needed in the calculation (altitude, latitude, albedo, overhead ozone column . . . ). The photolysis rates are calculates using the formula J = A * exp{B*(1-1/cos(C*SZA))}

10 Outputs

The model creates a NetCDF file containing an almost total description of the the trajectory integration contained in dataset dsn out (namelist CNTL). Information of important switches are written to global attributes starting with exp CHEM.

11 Examples

11.1 Single Box Model Calculation

The example script file boxtest_svo.exec is a single box model calculation using the SVODE solver. The initialization is given through the namelist input. The timesteps of trajectory, chemistry, and output are 1h, also the solver is called every hour (stiff differential equation solvers adjust the timestep into variable steps internally). The elementary reaction rates and the rate constants are written onto the output file.

11.2 Multiple Box Model Calculation

The example script chem_nh.exec is a 2-day calculation for a trajectory set covering the northern Hemisphere. The trajectory timestep is 24h, which is divided into nchemsteps=24 1-hour chemistry steps. The output is written every 12h. The internal solver timestep (ncdt) is 10 min. The chemical initialization is read from the file testData/initial_nh_97022712.nc. The ozone profile for the photolysis is the second month (ntime_o3) taken from the 2-D model output file testData/o3_nh_97T04.ref.nc.

12 Convergence problems with IMPACT

Due to the extension of CLaMS to the lower altitudes in the troposphere, there are increased problems that the solver IMPACT does not converge. This is likely due to the used family concept for the NOx species, for details see the master thesis of David Scheinmann (ASAD: Eine kritische Analyse der mathematischen Modelierung und numerischen Simulation der Reaktionskinetik in der Atmosphärenchemie, Uni Wuppertal, 2012). In some cases the divergence created impossible high concentrations and NaNs. These values will be filtered out and replaced by missing values. Unlike higher altitudes, for zeta < 300K that happens significanly often (status 25.04.2013). This problem needs to be solved. The use of the Newton Raphson Integation (method=3) is more stable