The CLaMS chmistry 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 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:

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 GNUMake:

gmake chem

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 gmake 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 Dimensions

The dimensions of the arrays on the program are defined in the include file include/paramc.h and include/paramcdef.h. A change in any of these dimensions requires a re-compilation of the program. The following table summarizes the dimension variables:

jpnl Maximum number of airparcels.

jpctr Number of model tracers/families in the model.

jpspec Number of species in the chemistry.

jpbk Number of bimolecular reactions.

jptk Number of termolecular reactions.

jppj Number of photolysis reactions.

jphk Number of heterogeneous reactions.

The dimensions in the include file include/paramcdef.h may be checked and eventually corrected using the make file command gmake check. It is also possible to define the maximum number of airparcels directly through a make command by typing e.g.

gmake jpnl=100 check.

4.2 Chemical Species

The names chemical species (and potentially chemical families) are provided in the file chch.d. 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.3 Chemical Reactions == The chemical reactions and their reaction rates are set in the files ratb.d (bimolecular reactions), ratt.d (termolecular reactions), ratj.d (photolysis reactions) and rath.d (heterogeneous reactions) located in the directory data. The number of reactions in these files must be consistent with those in include/paramcdef.h. This can be checked and eventually corrected by typing gmake check (using the program check dims).

5 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.

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 gmake stiff or gmake standard. These commands use also a program called check dims that checks the dimensions of the file paramcdef.h and corrects it automatically.

6 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.

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.).

7 Heterogeneous Chemistry

The aerosol particle microphysics and the heterogeneous reaction rates are calculated by version 6 of the code developed by Ken Carslaw. 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: sigma: standard deviation of the particle size log normal distribution

N0: initial liquid aerosol density [cm−3] (see also below)

ppbh2so4: Gas phase equivalent of the amount of sulfuric acid in the aerosol [ppbv]

n_ice: ice number density [cm−3] (if formed during the calculation)

n_nat: NAT number density [cm−3] (if formed during the calculation)

SAT melting allowed? 1=yes, 0=no (see Koop and Carslaw, Science 272, 1638-1641, 1996) Parameterization for reactions on NAT/SAT: Hanson and Ravi (1) or Abbat and Molina/Zhang (0)

Transformation matrix: 4×4-matrix containing the allowed (1) and forbidden (0, 99) phase transitions between the 4 phases. The right upper triangle stands for decreasing temperatures and the lower left triangle for increasing temperatures.

Saturation matrix: similar 4×4-matrix containing possible super-saturation needed for the occurrence of freezing.

gamma values: Reaction probabilities for the individual reactions on the different substrates. Gamma is calculated by the routine for those values set to 1.0 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.

8 Photolysis

8.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.

8.2 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))}

9 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.

10 Examples

10.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.

10.2 Multiple Box Model Calculation

The following example (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.