The ASAD atmospheric chemistry integration package and chemical reaction database

Computer Physics Communications ELSEVIER

Computer Physics Communications 105 (1997) 197-215

The ASAD atmospheric chemistry integration package and chemical reaction database Glenn D. Carver l, Paul D. Brown, Oliver Wild 2 Centre for Atmospheric Science, Cambridge University, Chemistry Department, Lensfield Road, Cambridge, CB2 I EIV, UK

Received 18 April 1997; revised 12 May 1997

Abstract The ASAD software package frees the atmospheric chemist from writing code to solve chemical time dependent equations and can be used for a wide range of problems including stratospheric chemistry, tropospheric chemistry or pollution studies, and for ID, 2D and 3D models. Since little programming is required, the effects of programming errors are greatly reduced and consistent results between different models ensured. Through the use of input files to list the desired chemical species and their properties and the chemical reactions, different chemical scenarios can be rapidly developed and tested without any additional programming. ASAD comprises a plug-in subroutine library to solve atmospheric chemistry time-dependent problems, a chemical reaction database and utility programs. ASAD has been designed to be coupled to any type of atmospheric model and be used as far as possible as a 'black box'. The code has been carefully optimized for performance and it will vectorize well on machines with this capability. It includes a choice of integrators based on published methods and supports approximations often used in atmospheric chemistry modelling, such as chemical families, fractional products and the steady state approximation. ASAD also supports the inclusion of user-supplied schemes for wet and dry deposition, source emissions, heterogeneous chemistry and photolysis. @ 1997 Elsevier Science B.V. PACS: 82.40.We; 02.60.Cb; 47.70.Fw Ke)words: Atmospheric chemistry; ODEs; Time integration; Chemical reactions; Stratospheric chemistry; Tropospheric chemistry; Kinetics

PROGRAM SUMMARY

1~tle of program: ASAD

Computer for which the program is designed and others on which it is operable: Cray YMP, Cmy J90, Sun SPARCstation, Silicon

Graphics workstation, IBM RS6000, Apple Macintosh, Pentium and Pentium Pro PC (any computer with a Fortran77 compiler). Operating systems under which the program has been tested: Unix

Catalogue identifier: ADGC

(Solaris 1 and 2, IRIX, AIX, UNICOS, Linux), MacOS System 7

Program obtainable from: CPC Program Library, Queen's Univer-

sity of Belfast, N. Ireland

Programming language used: FORTRAN77

t Corresponding author; e-mail: [email protected]. 2 Now at: Earth System Science, University of California, lrvine, CA 92697, USA. 0010-4655/97/$17.00 (~) 1997 Elsevier Science B.V. All fights reserved. PII S0010-4655(97)00056-8

198

G.D. Carver et aL/Computer Physics Communications 105 (1997) 197-215

Memory required to exec,tte w#h typical data: Highly dependent on problem; number of species, number of reactions and whether the code is being used to solve a ID, 2D or 3D problem. The number of words required for static (COMMON) storage is given approximately by the formula NIVORDS = npts(3ntr + 9nsp + 2nr + 5) + t,r( ]ntr + ]nsp + 16) + 21 (nsp + ntr), where npts is the number of gridpoints that ASAD is passed in each call, ntr is the number of chemical tracers transported by the calling model, nsp is the total number of chemical species and nr is the total number of chemical reactions used. Memory required will be slightly higher than this due to temporary arrays declared in subroutines. Number of bits h~ a word: No explicit sizing has been used. Therefore the code will run in 64 bit precision (REAL*8) on Cray (or other 64 bit word) computers and 32 bit precision (REAL*4) on other computers or workstations. However, on 32 bit computers, the use of compiler options to force real variables to 64 bits is strongly recommended for anything other than the simplest of problems. No. of processors used: 1 Has the code been vectorized: Yes, the vector length is over the spatial domain of the problem

Keywords: Atmospheric chemistry, ODEs, time integration, chemical reactions, stratospheric chemistry, tropospheric chemistry, kinetics Nature of physical problem and its method of solution ASAD has been designed to be the chemical part of an atmospheric chemistry-transport model. It is not a complete atmospheric chemistry model. It solves the time-dependent chemical equations. The problem is specified by inputting a list of chemical species, their properties and a list of chemical reactions. The user can also include some of the common approximations (such as chemical families) used in atmospheric chemistry models. ASAD has been designed to be treated as a black box and therefore has a clearly defined interface with its calling transport model. No code for the chemistry needs to be developed or changed and the user can alter the chemical reaction system quickly and easily to enable rapid development and testing. A number of published time integration methods for the chemical equations have also been built into the package. The user can specify which method to use, dependent on the desired accuracy and cost. Interfaces are provided in the code for the user to supply source and sink terms from emissions and wet and dry deposition. Although designed for atmospheric chemistry problems, with a suitable choice of integrator, ASAD could also be used for other chemical problems such as those found in combustion studies.

Peripherals ,tsed: None No. of bytes bz distributed program, hlcluding test data, etc.: 674048

Restrictions on the complexity of the problem The size of the problem is controlled by parameters in the code. Therefore, the available computing resources essentially restrict the complexity of the problem.

Distribution format: uuencoded compressed tar file Subprograms used: SVODE (from NETLIB) and D02EAF (from the NAG library) may be chosen as stiff ordinary differential equation integrators. All the necessary subroutines to use SVODE are included with ASAD. The use of the D02EAF subroutine requires the user to have access to the NAG library. Separate documentation available: A detailed user guide is included with the software or available from the corresponding author. There is also a World Wide Web page for this software at http://www, arm. ch. cam. ac. uk/acmsu/anad/index, html.

LONG

Typical running time Highly dependent on the problem; number of species, reactions, size of the spatial domain, choice of integrator and attached user code for photolysis or heterogeneous reactions. As an example, however, a box (single point) model of a typical stratospheric chemistry scheme with 30 species, 150 reactions, detailed photolysis and heterogeneous reaction schemes takes 4 CPU seconds for 10 days on a Sun HyperSPARC (150Mhz) workstation (rated at 4.7 SPECfp95) using a 10 minute chemical timestep and the IMPACT integrator.

WRITE-UP

1. I n t r o d u c t i o n A t m o s p h e r i c c h e m i s t r y m o d e l s a r e u s e d in a variety o f r e s e a r c h studies, c o v e r i n g a w i d e r a n g e o f spatial a n d t e m p o r a l scales; f r o m s t u d i e s o f t h e d i u r n a l v a r i a t i o n o f u r b a n air p o l l u t i o n to s e a s o n a l g l o b a l s t u d i e s o f s t r a t o s p h e r i c o z o n e o r c l i m a t e s i m u l a t i o n s o f t r o p o s p h e r i c m e t h a n e . T h e y vary f r o m s i m p l e m o d e l s that r e p r e s e n t a s i n g l e p o i n t in t h e a t m o s p h e r e ( s o - c a l l e d b o x m o d e l s ) , t h r o u g h to m u l t i - d i m e n s i o n a l L a g r a n g i a n

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

199

(trajectory) and Eulerian (gridpoint) models. Similarly, the complexity of the chemistry in such models can vary enormously. Typically, the more computationally expensive 2D and 3D models use fewer species and reactions than cheaper box or trajectory models. A common practice is to use box models with very detailed chemistry and accurate numerical integration methods to develop and test simplified chemistry schemes for use in the more expensive 3D models. It would clearly be an advantage to have a chemistry package that could be added to any model of the atmosphere and used essentially as a 'black box', but which was flexible enough that it could support very detailed and complex chemistry schemes and also support the use of the approximations often used in simplified schemes. It was for this purpose that we developed ASAD and its accompanying chemical reaction database and utility programs. The use of a single chemistry code ensures consistency between models and greatly reduces the development time since no programming is required to solve the chemical rate equations. Another very important benefit is that programming errors that might have arisen if the user had to write the code are avoided. Different chemical scenarios can therefore be rapidly developed and tested. The idea of a FORTRAN program to solve the ordinary differential equations (ODEs) arising from a chemical system is not new. The inspiration for ASAD and its reaction database was the DELOAD program and UMIST database of interstellar reactions presentcd in the Ph.D. thesis [I] and Ref. [2]. The FACSIMILE [3] and AutoChem [4] packages are similar in that they also allow atmospheric chemistry schemes to be constructed. However, unlike ASAD, DELOAD is a 'code-writer' in that the user supplies files listing chemical species and reactions and DELOAD generates the FORTRAN77 code to compute the chemical tendencies. The user must then form their own subroutines based on this code. On the other hand, FACSIMILE is a commercial package which does not use FORTRAN77 but a proprietary language and cannot therefore be coupled to existing atmospheric transport models. AutoChem (also inspired by DELOAD) is written in FORTRAN77, but like FACSIMILE only uses stiff integration methods to integrate species individually and does not support chemical families [5]. The high computational cost limits its use to single point (box) or ID models of the atmosphere. In contrast, ASAD takes these ideas further. It was designed to be a comprehensive package of 'plug-in' subroutines that could be used in all types of atmospheric models. Whilst using similar input files to DELOAD, it does not output code but solves the equations for the chemical part of the atmospheric chemistry model and passes the solution back to the calling model. It also supports the use of the common approximations, such as chemical families (see Section 2.1), steady state species (see Section 3) and fractional products (see Section 4.4), used in atmospheric chemistry models. In addition, ASAD has been written in FORTRAN77 so that it can be modified or customised if necessary. The following sections first introduce some basic ideas behind the ASAD package. Then the process of choosing a set of chemical species is discussed, followed by a description of the ASAD chemical reaction database and how to use it to construct a chemical scheme. After an overview of the ASAD program structure, we present some comments and mention likely future developments. It is not intended that this long write-up provide all the details necessary for using ASAD. A detailed user guide which takes the reader step-by-step in setting up an atmospheric chemistry model with ASAD is available with the package.

2. Basic principles An atmospheric chemistry model can be generally considercd to consist of two main sections, a dynamical/transport component and a chemistry component, as shown schematically in Fig. 1. The dynamical part is responsible for the transport of the chemical species. The chemical part of the model is responsible for computing and integrating the chemical rates of change. To do this, it may use information from a photolysis scheme, possibly also a heterogeneous scheme (for example, reactions on cloud droplets), and perhaps also a surface scheme for emissions of source gases and wet and dry deposition schemes. Information can flow both

200

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

I Chemistry scheme (ASAD)

Atmospheric transport or dynamical model

Photolysis scheme

Heterogeneous chemistry scheme

Surface emissions & deposition model

[

Fig. I. Schematicillustration of an atmosphericchemistrymodelshowingthe two main components,the transportmodel and the chemistry model.The transportmodelmay be a box, 1D, 2D or 3D. The chemistrymodelmay optionallyincludeuser-suppliedschemesfor photolysis. heterogeneouschemistry,emissions and deposition. ways between the packages. As indicated in the figure, ASAD is designed to take the place of the chemical part of the model and makes use of user-supplied photolysis, heterogeneous, deposition and emission schemes. Consider the total rate of change of a species, y, dy

dy

dy

The transport component of the atmospheric chemistry model will have the responsibility for the first term on the r.h.s, of Eq. (1), whilst ASAD takes the responsibility for computing the remaining tendencies. ASAD computes the second term from tabulated reaction rates stored in input files and photolysis rates from the user-supplied photolysis scheme, whilst the third term is computed by calling the other user-supplied packages. The user may of course choose to use ASAD just for the second term and use their own method to integrate the deposition and emissions tendencies outside of A S A D . ASAD therefore computes the total chemical rate of change of each species (the sum of the second and third terms on the r.h.s, of Eq. ( I ) ) and integrates it forward in time. This rate of change of a species, y, is given by dy = p _ L y _ ly 2 + E dt

Dy,

(2)

where P is the total chemical production of the species, L is the chemical loss rate, I is the loss.rate when the species self-reacts, E is a production rate from emissions from the surface and D is the loss rate from dry and wet deposition. ASAD will compute the rate of production ( P ) and the total chemical loss rate ( L y - ly 2) from the information in the reaction rate input files. If the user wishes to include a source of emissions and loss through deposition, they must supply replacements for dummy subroutines in ASAD, otherwise these terms are zero and ignored by ASAD. At this stage it is necessary to distinguish between the chemistry variables that are known and stored in the calling atmospheric transport model, which may be either chemical species or families, and the chemical species that ASAD uses internally. We will refer to the former as 'tracers' or 'families' since they are transported and the latter as 'species'. A S A D only integrates the tendencies for those tracers and families known to the calling model. For example, some species may be specified as constant or in steady state, in which case they are not integrated in time or known to the calling transport model, being defined within the ASAD subroutines only.

G_D. Carver et aL/Computer Physics Communications 105 (1997) 197-215

201

2.1. Chemical families A unique feature of ASAD is its ability to support the use of chemical families (see [5] ). A chemical family involves grouping related species so that the short lived chemical species are parametrized in terms of a longer lived species and the total family concentration is transported. By short-lived, we mean that the species lifetime is much smaller than the model timestep, so that the concentration of the species is purely determined by its chemical relationship to the long lived species, rather than by transport. For example, an odd oxygen family, Ox, is often used in atmospheric chemistry modelling which contains the species Ox = 03 + O(ID) + O ( 3 p ) ,

(3)

where ozone, 03, is the long lived species and O(ID) and O(3p) are short lived and assumed to be in steady state. The calling atmospheric model will transport the variable Ox whilst ASAD will separate it into its members to work out the rate of change of Ox (and other species that depend on the members of the family) and integrate it before returning the updated values to the model. The use of chemical families is computationally advantageous for several reasons. First, they reduce the stiffness of the chemical equations, Eq. (2), and allow the use of longer timesteps or less costly integration methods. Second, it reduces the number of tracers that need to be transported in the model. In general a chemical family, f , is defined by

f = nmYm + E niYi' i

(4)

where n denotes the number of 'odd atoms' of the family member and y denotes the concentration of the family member. The subscript m denotes the main family member whilst the subscript i denotes the members assumed to be in steady state. For example, for the chlorine family ClOx with the species Cl, ClO and Cl202 as members, Cl202 contributes two odd atoms of chlorine (Cl) to the ClOx family whereas there is only one in ClO (because C1202 photolyses to give Cl + Cl + 02). Note that this is not the same as counting the number of atoms in the molecule. For another example, in the odd oxygen family, Ox, ozone (03) has 3 oxygen atoms but only contributes one odd oxygen atom to the Ox family. The calling transport model will supply the value of the family, f . The subroutine FI'OY in ASAD is responsible for computing the concentrations of the family members, y, before the computation of the chemical tendencies. In order to do this, an iterative Jacobi-like procedure is followed in which we write the species as ratios, by dividing Eq. (4) by the main family member and inverting R,,I =

1

lira -]- E i niRim

,

(5)

where

Rmf = Ym/ f

(6)

and

Rim = YilYm.

(7)

The concentrations of the species Yi and hence the ratios of the family members, Rim, are calculated from the assumption of steady state. From Eq. (2), assuming no quadratic" loss term,

Yi=

P+E L+D"

If l in Eq. (2) is nonzero then a quadratic expression is solved for y.

(8)

202

G.D. Carver et aL/Computer Physics Communications 105 (1997) 197-215

The iterative procedure initially partitions the family so that y,, = f and all the Yi = 0. ASAD computes the production and loss of each species using the supplied input reaction rate files. Then the ratios Rim are computed via Eq. (8) followed by the ratio Rmf using Eq. (5). This is used to give a new value of Ym using Eq. (6) and then new values of Yi using Eq. (7). The procedure repeats until the values of all the y concentrations converge or the number of iterations specified by the user have been done. Convergence is achieved when the difference between two iterations gives a value less than RTOL x y or when y is less than ATOL, where RTOL is the relative tolerance and ATOL is the absolute tolerance. The user is free to chose the values of RTOL and ATOL. Since the values of the y are bounded by the use of ratios, the family members' concentrations converge within a small number of iterations. However, during this iteration, the concentrations of steady state species which are not a member of a family must also be computed using Eq. (8). The concentrations of these species are usually small and often strongly dependent on each other as well as the family members. Consequently, it is generally the convergence of the steady-state species which determines the number of iterations required. The user must experiment to determine the optimum number required for the desired accuracy. ASAD will use all the reactions that a specie s takes part in to compute the ratios. However, some reactions may have an insignificant contribution and in 'hand-coded' chemistry schemes these reactions are sometimes omitted to reduce computational costs with little effect on the results. At the moment, ASAD does not provide any facility for this since the reactions to be omitted are highly dependent on the chemistry scheme, the choice of families, reactions and topic under study. ASAD also has the facility to treat species as tracers when their chemical lifetimes are long enough, but as members of a family when their lifetime becomes short compared to the timestep. An example is the species NO3 in the lower troposphere. By day, its lifetime is of the order of seconds but by night it is roughly 12 hours or more. ASAD can include the species in the family when the lifetime is short and integrate it as a separate tracer when its lifetime increases. Note that in this case, the calling model is transporting this species as a tracer, even though at times it is a family member. However, the use of this option can introduce a computational overhead. This is because, for some of the time integration subroutines, it is possible to precompute coefficients once only at the start of the run. If a species moves in and out of a family during the run, these coefficients must be recomputed at each timestep. 2.2. Process splitthzg of chemistry and transport The chemical lifetime of a species or family often requires a smaller timestep for the chemical integration than is needed by the transport model alone. Therefore the integration of the transport tendency (the first term on the r.h.s, of Eq. (1)) is usually separate from the chemical integration of that tracer. For example, the transport model may use a 'dynamical' timestep of 30 minutes whilst the desired accuracy of the chemical solution may require a 10 minute 'chemical' timestep. ASAD would then perform 3 internal chemical sub-steps before returning to the calling model. During these sub-steps, the physical variables passed by the model, temperature and pressure, do not change. ASAD passes back the updated values of the tracers and the time averaged chemical tendencies. It can therefore accommodate models which processes split the chemical integration and simply require the updated values, and also models that add the chemical tendencies to the transport tendencies in a nonprocess split fashion. 2.3. Vectorization and performance ASAD has been designed for use in any atmospheric model with any number of spatial dimensions. To accomplish this, it has arrays with only one spatial dimension which is used for the innermost loop in the code wherever possible. This enables the code to veetorize on computers with this facility. The spatial dimension

G.D. Carver et aL /Computer Physics Communications 105 (1997) 197-215

203

could be longitude, latitude, height or a combination of these. For instance, in our 3D atmospheric chemistry transport model we pass a longitude-height slice of the model grid to ASAD to achieve good vector length (approximately 1000 elements), ttowever, some of the integrator subroutines (particularly the stiff integrators) built into ASAD to solve the chemical rate of change equations, will only compute a single gridpoint at a time (see Section 5.3 for more details). Considerable effort has been taken to ensure that the FORTRAN code executes efficiently on vector computer architectures, in order to make ASAD attractive to users with expensive 3D atmospheric chemistry models. The code has been carefully written to avoid unnecessary calculations wherever possible. The PERFVIEW facility on a Cray YMP was used extensively to analyze memory references to avoid memory conflicts and improve pipeline performance. Optimization techniques such as loop unrolling have also been used to improve the algorithmic efficiency. Whilst all the code optimization for ASAD was performed on a vector computer, ASAD also performs well on cache based microprocessors (e.g. UNIX workstations). However, long vectors are not appropriate for these types of processors and we encourage the user to pass ASAD a much smaller set of points. As microprocessor cache designs vary considerably, it is impossible to give advice about the ideal number of points to pass. Some experiments with an atmospheric model on microprocessors were carried out in [6] where they noted that good performance was obtained with vector (inner loop) lengths of about 16 on a Cray T3D node (DEC alpha processor).

3. Choosing a species set

Two types of files are read by ASAD. The first is the chch.d file which contains a list of the chemical species the user wishes to use and their properties (recall that we distinguish between 'tracers' transported by the model and passed to ASAD, and 'species' that ASAD uses internally). The other type of file is the set of reaction ratefiles which we describe in the following section. Together, these files define the chemistry scheme. The user must run a program called SELECT which takes the user's chch.d file and, using a given reaction database, scans that database and extracts only those reactions involving the species listed in the chch.d file. This process is also described in more detail in the following section. When creating a chch.d file, it is very important that the species must be identified by the names given to them in the A S A D reaction database, otherwise the species will not be recognised. Species names are restricted to 10 characters. We have devised a nomenclature system for species whose commonly used names or chemical formulae are longer than this imposed limit (see Appendix). A list of species currently featured in the ASAD reaction database is given in the file species.d, which is included with the package. An example of a chch.d file is shown in Table 1 and the meaning of each column is given here. FORTRAN list directed I / O is used to read this file so precise spacing of the columns is not important. Column 1. An integer is used to identify each species for the user. This is ignored by ASAD but is used by some of the utility programs. Column 2. This column contains the species names, exactly as used in the ASAD reaction database and listed in the species.d file. Column 3. This indicates the number of "odd atoms" the species contributes to the family it belongs to (see Section 2.1). This is ignored if a species does not belong to a family. Column 4. This is a two character code that is used to specify how the species should be treated in ASAD. These must be in upper case. The following options ,are available: TR

The species is integrated as an independent chemical tracer (not part of a family, not constant or in steady state e.g. N 2 0 ) . This tracer must be passed to ASAD by the calling model.

FM

The species is to be treated as a member of the family named in column 5. This family is one of the

204

G.D. Car~'er et aL/Computer Physics Communications 105 (1997) 197-215

Table 1 Example chch.d file format Species number

1

Species name

No. of Odd atoms

O(ID)

2 3 4 5 6 7 8 9 I0 lI 12

O(3P) 03 CI C1202 OCIO CIO H202 OH CH4 N20 MeOOH CO2 N2

13 14

Treatment

Family

'FM' 'FM' 'FM' 'FM" 'FM' 'Fr' 'FM' 'SS' !SS" 'TR'

'Ox' 'Ox' 'Ox' 'ClOx' 'CIOx' 'CIOx' 'CIOx' ' ' ' ' ' '

9T R '

9 '

'TR'

' '

"CT'

~ '

~CT'

9 9

tracers defined in the atmospheric transport model. F'r

The species is to be integrated as an independent tracer unless its lifetime during any one particular model timestep is less than half the chemistry timestep (the two timesteps may be different - see Section 2.2), in which case it is treated as a member of the family specified in column 5- Both the species and its associated family will be tracers in the calling transport model (note, however, that the concentration of the family stored in the transport model should not include this species)

SS

The species is neither stored in the calling model nor part of a family but, for the purposes of the chemistry, can be assumed to be in steady-state. This species will only be known internally to ASAD.

CT

The species is not stored in the model but, for the purposes of the chemistry, is assumed to have a fixed volume mixing ratio. The species for which this may be used are: C02 (3.5 x 10-6), //2 (5.0 x 10-7), 02 (0.20945) and N2 (0.78084). These concentrations are 'hard-wired' into the code. Other species and/or concentrations may be added to this list by modifying ASAD source code.

CF

Like the 'CT' species type, this species is also treated as a constant, but one which can have a different value at each point in the spatial domain. In order to set this constant field, a subroutine which the user must supply, INICNT, is called at the appropriate place in ASAD.

Note that only the variables which represent the families and species integrated separately (types FM, FF and TR) are passed between the calling model and ASAD. Family members, steady state species (SS) and constant species (CT and CF) are not passed. Column 5. This column contains the name of the family (for types FM and FT), if required, to which the species belongs. A blank space indicates the species is not a member of a family. Family names must be 10 characters or less, and must not be the same as any of the species names. ASAD will always use the last named species of a particular family as the major species in that family. Therefore, it is important to order the species which are family members correctly in the chch.d file.

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

205

Table 2 Example bimolecularratefile Index

Reactants

8 9 10 11

Br Br BrO BrO

Products 03 OCIO BrO BrO

BrO BrO Br Br2

02 CIO Br 02

02

ko

a

1.70E- l I 2.60E- l I 9.20E-13 1.80E-13

0.00 0.00 0.00 0.00

fl 800.0 1300.0 250.0 250.0

Notes

C C

4. Chemical reaction database 4.1. Ratefiles

There are four master ratefiles in the ASAD chemical reaction database, one for each category of reactions: bimol.d Bimolecular reactions; trimol.d Termolecular and unimolecular reactions; photol.d Photolysis reactions; hetero.d Heterogeneous reactions. The reactions have been compiled mainly from the IUPAC Supplement IV assessment [7] and the JPL reaction handbook [8]. Reaction data in refereed papers have also been used on occasion for reactions that have not previously been measured or for new reactions. These are indicated in the ratefile. In general, however, we felt it was appropriate for the master ratefiles to adhere to carefully reviewed kinetic data. However, the user is free to create their own master ratefiles and adopt their own policy with regard to usage of published kinetic data. It is our intention to update the ASAD reaction database when new assessments are published. All ratefiles are in ASCII and use an easy to read format. In any ratefile all ASAD programs will treat a line beginning with a ' # ' as a comment and ignore it. 4.1.1. Bimolecular ratefile

A few lines from the bimolecular ratefile are shown in Table 2. The ratefiles are laid out in columns; the first is the reaction number in the file, the next 2 contain the reactants, the next 3 columns contain the products and the subsequent 3 columns contain the parameters needed to calculate the rate coefficient for the reaction. The last column is a character string used for comments, or to refer to notes at the bottom of the ratefile. For example, to list references for those reactions taken from sources other than [7] and [8] or to denote kinetic data which is only known to an upper limit. Blank entries are acceptable where a reaction does not have enough products. If it is necessary to have more than 3 products, the user must alter a parameter in the ASAD code and reformat the ratefiles accordingly. A utility tool is provided for this purpose. ASAD uses a fixed format FORTRAN read and write for these ratefiles so positioning of the columns is important. Species names must be no IOhger than 10 characters and must be a recognised species name in the ASAD nomenclature system (see Appendix). For the bimolecular ratefile, columns 7 (k0), 8 ( a ) and 9 (/3) are used to compute the rate coefficient, k, from the Arrhenius expression T

k=k0(3--6-6)

exp(-~Tfl )

(9)

in units of cm-3s - I , where T is the temperature in Kelvin. The subroutine BIMOL is used to compute all the bimolecular rate coefficients.

G.D. Carveret aLIComputer Physics Communications 105 (1997) 197-215

206 Table 3 Example termolecularratefile Index

Reactants

1 2

BrO C1202

NO2 rn

Products

f

kl

BrONO2 m CIO CIO

327.0 0.6

4.70E-3 i 1.35E-05

ai -3.10 1.00

/~l

k2

0.0 1.70E-I I 8720.0 1.80E-I-15

a2 -0.60 0.00

/72 0.0 8450.0

There are exceptions to this general formula. For example, the reaction between CO and O H is pressure dependent. However, there are only a handful of reactions that do not use Eq. (9). Rather than alter the ratefile format, these reactions are flagged in column 10 and notes on how to compute the reaction rate are given at the end of the ratefile. For such reactions listed in the A S A D reaction database, code specific to that reaction has been included in BIMOL to give the correct rate coefficient. A question mark '?' in any of the product columns indicates that the products are unknown even though a measurement for the reaction rate exists (this applies to all ratefiles). We will discuss the significance of these reactions in choosing a reaction subset in the following section. There are reactions which have two or more possible branches. Each branch must have a separate entry in the ratefile. Where the branching ratio is unknown, we have assumed that the branches occur at equal rates i.e. "50 : 50 ratio in the case of two branches, 33 9 33 : 34 for three branches and so on. This is accomplished by dividing the total measured ko parameter by the number of branches. 4.1.2. Termolecular ratefile Although called the termolecular ratefile we also include unimolecular reactions in this file as unimolecular and termolecular reactions can be thought of as consisting of the same elementary processes and can exhibit varying order depending on pressure (see [9]).The termolecular file follows a similar format to the bimolecular ratefile and an example is shown in Table 3. For all reactions a third body is assumed. For unimolecular reactions, the character ' m ' appears in place of the second reactant. The number of parameters needed to calculate the rate coefficients for terrnolecular reactions is greater than in the bimolecular case. The temperature and pressure dependent rate coefficient, k, is found from the rate constants for the low pressure, ko, and high pressure, ko~, limits by

1

k0tM]

"~ _ (, + (log to[ ~,l/k~ ) 2) - '

(lO)

in units of cm-6s - I where [M] is the concentration of the third body (molecules cm -3) (using the notation of [7]). Typically, Fc is reported as a constant value, however it sometimes has a significant temperature dependence. To include this we use the following approach. If f is the parameter in column 6, when this value is less than 1, Fc otherwise, Fc = e x p ( - T / f ) where T is the temperature (in Kelvin) (see [7] for more discussion). The low pressure rate constant, which is also first order in pressure, is found from k0=ki \300)

exp

,

(11)

where hi, tel and fll are the parameters given in columns 7-9 of the ratefile (Table 3). Similarly, the high pressure limit, k, is found from

k~=kzt-~-~j

exp

,

where k2, te2 and/32 are the parameters given in columns 10-12 of the ratefile.

(12)

G.D. Car~'eret al./Computer Physics Communications 105 (1997) 197-215

207

One complication is that the low pressure limit, k0, can sometimes be reported as dependent on the nature of [M], i.e. a different rate is given for N2 and 02. If only the kl parameter is different, then we take the weighted mean of the specified values for kl. However, if the exponential factors, a l , are different, we list the reaction as having two branches in the ratefile and scale the parameter kl by the proportion of [M] in air. In both cases, notes are added to the ratefile for such reactions.

4. !.3. Photolysis and heterogeneous ratefiles The photolysis and heterogeneous ratefiles are similar to the bimolecular ratefile. The columns used to specify the reaction rate are unused in these master ratefiles. This is because it is impossible to define a format which suits all photolysis and heterogeneous schemes. These columns are available to the user for specifying parameters for each reaction as part of their scheme. ASAD will read all the columns in the ratefile and the parameters are available in COMMON blocks in the code. 4.2. Selecting reactions Once the choice of species for a chemical scheme has been made, the reaction set must be defined. The utility program SELECT is designed to scan the supplied master ratefiles (or alternate ones) and select only those reactions involving the species listed in the chch.d file. This subset of reactions is written to new ratefiles in the same format as the master ratefiles. Two types of output from SELECT are possible; either an 'open' reaction set or a 'closed' reaction set. By 'open' we mean that reactions are selected if reactants are listed in the chch.d file. In a 'closed' set both reactants and products must be listed in the chch.d file. There are advantages and disadvantages with both types of scheme. An open reaction set may lead to nonconservation since products which are not listed in the chch.d file will not be integrated. On the other hand, a closed set of reactions will conserve but may exclude reactions that are known to be important loss mechanisms for particular species, even though the products of those reactions are unknown or irrelevant to the overall aim of the chemistry scheme. We feel that the best course of action is tO take the open set of ratefiles and modify them as necessary. If the open reaction set reveals species that ought to be included in the chemistry, these can be added to the chch.d file and SELECT rerun with the revised list of species. However, the user is of course free to make the final choice.

4.3. Modifying the reaction ratefiles As well as the need to study the differences between the open and closed set of ratefiles, there are several more reasons why they may have to be modified before being used by ASAD. For instance, if a reaction is included which has unknown products but is generally considered to be an important loss mechanism for one of its reactants, the user can make an assumption about the products formed. The user may also wish to exclude those reactions where the measured rate is listed as an upper limit only (indicated in the comment column in the ratefile) since their inclusion may introduce unrealistic effects. Or it may be desirable to retain them in the final ratd'files in case new kinetic data subsequently becomes available. In this case, the reactions can effectively be turned off by setting all the rate coefficient parameters (columns 7 to 9) to zero. The user may also wish to attribute a fraction of the reaction to one of the specified products (see next section). Another reason for editing the ratefiles before using them in a chemistry model is to reduce the number of species to decrease the computational cost. This is distinct from including species in families. Species that might initially be included in the chch.d file but might be subsequently removed are those that react with either 02 or N2 almost instantly upon formation, and are of relatively little interest on their own. For example, CH3 formed in the oxidation of Ctt4 reacts very quickly with 02 to form C H 3 0 0 (written as M e O 0 in the ASAD database; see Appendix). The species Ctt3 could then be removed from the chemistry scheme completely by deleting it from the chch.d file and carefully editing the ratefiles. For this example, all reactions in which CH3 is

208

G.D. Carveret aL/Conzputer Physics Communications 105 (1997) 197-215

one of the reactants should be deleted and in those reactions where it is formed as a product, the species M e O 0 should replace CH3. If a species deleted in this way forms more than one product, they must all be included as reaction products. For example, HCO is another species frequently ignored for stratospheric chemistry schemes. It reacts to form CO and H02, both of which must therefore be listed as products. It should be clear that this type of manipulation of the reactions in the ratefiles should be done with extreme care and always well commented. It is important to note that the ratefiles contain information about the reactants and products of a reaction and, in the case of the bimolecular and termolecular ratefiles, the parameters necessary to calculate the reaction rate coefficient. It is therefore straightforward to modify the chemistry scheme by altering the parameters in the ratefiles, without altering the ASAD code or recompiling it. For example, reactions may be turned off just by setting the reaction rate to zero. Although the reactions in the master ratefiles will be bona fide ones from published references, the user is entirely free to add pseudo reactions to represent or parametrize certain chemical processes or for testing sensitivity to newly proposed reactions. 4.4. Fractional products Long and complex reaction chains can make ASAD expensive to run. Such long chains are frequently found in atmospheric hydrocarbon modelling and often approximated with reactions in which the final products are written as a fraction of the production rate of the reaction. There are no fractional products used in the master ratefiles, but the user may include them in the ratefiles produced by SELECT. For example, with a reaction A+B~C+D

witharatek

(13)

the rate of change of C from this reaction is dC/dt = f k A B , where f is the fraction that contributes to the production of C. At present, f must be a constant. A S A D does not directly support any variation of f due to temperature or pressure. If this is necessary, the user must edit the ASAD source code. Whilst strictly nonconservative, the use of fractional products can greatly reduce the number of reactions and therefore reduce the computational cost of the chemistry scheme. On a separate line under each of the products, a fraction expressed in FORTRAN ' F ' format is placed. Instead of a number in the first (index)column of the ratefile, a '$' character is used as a continuation marker. Positioning of these entries is important as A S A D uses fixed format I / O to read these ratefiles. Examples are given in the A S A D userguide that accompanies the package.

5. Overview of p r o g r a m structure 5.1. Model-ASAD interface Fig. 2 shows the relationship between ASAD and the calling atmospheric model and tile files that ASAD will read. There are two ASAD subroutines that the atmospheric model must call, CINIT and CDRIVE. The subroutine CINIT is called once only at the start of the run. It is responsible for initialising the chemistry, and will read the chch.d file and the reaction ratefiles previously set by the user as described above. It will also initialise those species designated as constant in the chch.d file. All the options provided by ASAD, such as the chemistry timestep, choice of integrator, whether emissions are included and so on, are controlled by variables passed through the argument list to CINIT. Full details of the arguments to CINIT are provided in the ASAD userguide that accompanies the package and as comments in the source code. The subroutine CDRIVE is the subroutine that is called by the atmospheric model to carry out the integration of the chemistry (and return the chemical tendencies). CDRIVE will be called at least once every model timestep. The actual number of calls depends on the length of the vector passed to ASAD. For instance, in a

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

209

ASAD

Atmospheric model

g m

i

~

read files: chch.d and ratefiles

9 .

Initialise common blocks Call user supplied initialization code

m m

Loop over model timesteps

~

C

D

R

I

V

E

~

I

~

"

Compie reaction rates

Compute species from input tracers

m m

-IJ

-" .~ -"

loop over chemistry timesteps

Compute product and loss Compute tendencies Integrate species I

Fig. 2. Flow diagram of ASAD and interface with the calling atmospheric model.

3D model in which ASAD is passed a longitude-height slice, CDRIVE will be called once for each latitude (see Section 2.3). CDRIVE should be passed the tracer concentrations and the temperature and pressure of the gridpoints. If the chemistry is not being integrated, the concentrations are unchanged and ASAD returns the chemical tendencies for these tracers. If the chemistry is being integrated (normal usage), then the updated chemical tracers are returned together with the chemical tendency averaged over the integration period (i.e. the model timestep). CDRIVE will accept tracers in either volume mixing ratio (vmr) or number density (molecules cm-3). Internally, ASAD subroutines always use number density as this removes the density dependence from the kinetic terms and hence reduces the number of multiplications required. CDRIVE will convert vmr to number density if required. CDRIVE will then compute the bimolecular and termolecular reaction rates from the kinetic data read by CINIT. If the user has specified the use of emission rates or deposition rates, then user-supplied subroutines will be called to set these. Note that all these rates are dependent only on the model variables, such as temperature and pressure and are therefore only computed once each time CDRIVE is called. The chemistry time loop in CDRIVE begins with the computation of the photolysis and heterogeneous reaction ~ites. These must come from user-supplied subroutines since such schemes are dependent on the problem under study. Note that both are called inside the chemistry time loop. The photolysis reaction rates depend on the zenith angle which in turn depends on the local time at each gridpoint which changes during the chemistry timestepping. Similarly, the heterogeneous reaction rates may depend on the concentrations of the chemical species. The next stage is to set the species concentrations, as listed in the chch.d file, from the tracers and families passed from the calling model. This is done by the subroutine FTOY. For tracers which are simply species, the values are copied. For families, the family members are separated as described previously in Section 2.1. During the calculation of the family members, the values of the steady state species are also computed. Once all the concentrations of the species are known, the production and loss terms are computed and hence the chemical tendencies of each species. The tracer/family tendencies are then formed and the user's

210

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

choice of integrator is called to advance the chemical concentrations of the tracers and families. More details of the subroutine structure called by CDRIVE are given in the ASAD userguide that accompanies the package. 5.2. User-supplied subroutbtes

ASAD includes dummy subroutines for the initialisation and calculation of photolysis rates, heterogeneous reaction rates, emission and dry and wet deposition and to set constant species of type 'CF'. If the user wishes to include these effects, replacement subroutines must be supplied. There are two types of user-supplied subroutines for the photolysis, heterogeneous chemistry, deposition and emissions; an initialisation routine which CINIT will call, and a routine to compute the actual rates during the timestepping. The user need not use the initialisation routine but must supply a replacement subroutine for the ASAD routine which computes the rates. For example, CINIT will call the routine INPHOT to initialise the photolysis scheme. The ASAD version is a dummy subroutine but the user may use it to do any setting up of their photolysis scheme. The subroutine PHOTOL will be used by ASAD to set and return the photolysis rates. If the user wishes to use a photolysis scheme, PHOTOL should be replaced by a subroutine to interface ASAD to their own code. The user must read the ASAD userguide in order to find out how to correctly place these rates in the A S A D COMMON blocks. 5.3. bltegrators

ASAD has a choice of built-in integrators which can be used to solve the ODEs in Eq. (2). Integrators designed for stiff and nonstiff problems are available. The choice of method to use will depend on the required accuracy of the solution and the computational cost of the chemistry. For instance, if families are not used, the problem is generally stiff. A chemistry based on families used in a 3D atmospheric transport model could use one of the computationally efficient nonstiff solvers. The first nonstiff integrator is the Quasi-Steady State Approximation (QSSA) as described in [ 10]. This is a cheap and explicit time integration method which is well known in the literature. It is included in ASAD as it is a popular method. However, it can suffer from poor conservation and we do not recommend its regular use without careful testing and possibly modifications to the source code to improve the conservation of the scheme. We have included the IMPACT algorithm as described in [ 11 ]. This is an accurate implicit time integration method with good conservation properties, whilst efficient and highly suited to modelling stratospheric and tropospheric chemistry in 2D and 3D models when families are used to reduce the overall stiffness of the chemical system. Two highly accurate but expensive integrators suitable for stiff problems may also be chosen by the user. The first is the NAG library routine D02EAF [ 12] and the second is the SVODE integrator [ 13] from the NETLIB subroutine library. The NAG integrator can only be used if the NAG library is available; it is not provided as part of A S A D . The call to this routine is commented out by default to avoid compilation problems. The SVOlgE subroutine is included with ASAD. The SVODE integrator can solve both stiff problems, using backward differentiation formulae, or nonstiff problems using an Adams method. Both of these options are provided in ASAD. It is very important to note that if either the SVODE or NAG integrators are selected, they will only solve the equations for a gridpoint at the time. Since the inner loop is over gridpoints in the ASAD code, this means that the ASAD subroutines which represent most of the computational work will not vectorize and the performance will suffer on a vector machine. We therefore do not recommend the use of these integrators for multi-dimensional problems on vector processor architectures. A typical scenario for the development of an atmospheric chemistry scheme would be to use either the SVODE or N A G integrator initially for testing and revising the chemical species and reactions in a box model

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

211

of the atmosphere. It could also be used to compare the results when chemical families are used. Then, the computationally cheaper integrators could be used after development when the chemistry scheme is used in a 2D or 3D model.

6. Test runs There are several tests to be conducted using the ASAD software to ensure that it is working correctly. The first test checks that the ratefiles can be read correctly. The second and third tests involve running ASAD on a sample problem using both stiff and nonstiff integrators. Our sample problem is a stripped-down stratospheric chemistry scheme run at a single point in the atmosphere for 10 days. The scheme includes daylight averaged photolysis rates (computed from our detailed stratospheric photolysis scheme) but does not include any stratospheric heterogeneous processes. Note that the stratospheric scheme is supplied for testing only and should not be used for realistic stratospheric simulations, although it might form a starting point for building a realistic stratospheric scheme. The first test, uses the SELECT program to duplicate the supplied rate files for the box model test runs. The Unix job script uses these ratefiles as master ratefiles and then runs SELECT given the same list of species that created the ratefiles initially. The output ratefiles should be identical to the input ratefiles. In the second and third tests (for which output is provided), a sample box model code has been supplied. A Unix job script runs this model twice, the same set of species are used in both runs but their treatment in ASAD is altered. In the first run, the c h c h . d file is configured to use families for the chemical species and the IMPACT integrator is used. In the second run, no families are used, each species (apart from those specified as constant) is integrated using the SVODE integrator. Differences to the 3rd significant figure may be expected because of the tolerances used in the test runs and the differences between hardware and compilers.

7. Summary and future development In this paper we have described the ASAD atmospheric chemistry package and reaction database. ASAD is a collection of subroutines which can be incorporated into any model of the atmosphere and used as a 'black box' to solve the chemical rate equations. The chemical scheme is defined by a set of input files listing the required species and their properties and reaction rate files. A number of utility programs are provided for use with the ratefiles. No programming is required by the user to compute and solve the time dependent equations. This enables much faster development and testing of chemistry schemes than when coding chemistry schemes for atmospheric models from scratch. It reduces the risk of programming errors and ensures consistent results across different models. There is a simple interface between ASAD and the calling transport model with dummy subroutine 'placeholders' so that the user may include their photolysis, heterogeneous and surface emission and deposition schemes. A choice of ODE integrators is built into the package. As the kinetic reaction rate data is stored in input files no recompilation of the code is necessary when these data are changed. Rate data can be added and updated quickly and sensitivity experiments can be easily performed. ASAD supports methods commonly used in atmospheric chemistry, such as chemical families, steady state species and the use of fractional products. Further developments in ASAD are largely user driven. At the t i r e of writing, we plan to investigate new numerical algorithms, such as additional ODE integrators to improve the performance and choice for the user and, perhaps more significantly, to provide a massively parallel processing (MPP) framework for ASAD. Atmospheric chemistry is an application that is highly suited to MPP computers and this is currently a high priority.

212

G.D. Carver et aL/Computer Physics Communications 105 (1997) 197-215

Within our research group ASAD has proved to be a valuable tool and a number of stratospheric and tropospheric chemistry schemes have been developed (see [14] for example). It has greatly reduced the development time for our atmospheric chemistry models because of the decreased programming and debugging required. When developing both the stratospheric and the tropospheric chemistry schemes for our 3D atmospheric chemistry transport model, we were able to rapidly develop and test a number of different schemes using a box model with the NAG integrator, then run our 3D model by simply using the same species and reaction rate files for ASAD in the 3D model, with no further code development. We believe that ASAD may prove useful in other types of studies such as pollution studies and possibly also computational fluid dynamics studies such as detonations, if coupled with the appropriate integrators.

Acknowledgements This work has been carried out as part of the UK Universities Global Atmospheric Modelling Programme (UGAMP) within the Centre for Atmospheric Science. UGAMP is a NERC funded community research programme. The Centre for Atmospheric Science is a joint initiative between the departments of Chemistry and Applied Mathematics in the University of Cambridge. P. Brown is grateful to the UK Department of the Environment for financial support. The authors would like to thank all those members of the Centre for Atmospheric Science who contributed to the development of ASAD and its testing. In particular we thank; Martyn Chipperfield, Sarah Goodchild, Jamie Kettleborough, David Lary and Dudley Shallcross. We would also like to thank Peter Stott for discussions on the IMPACT timescheme, Tom Millar and Lida Nejad for copies of the UMIST reaction database and the DELOAD program, and our scientific colleagues using ASAD whose questions served to improve our code and documentation.

Appendix A. Species naming convention A complete list of the species known to ASAD can be found in the file species.d provided with the package. To identify a species in a chemical scheme uniquely, it may be necessary not only to supply the full molecular formula, but also isomeric information. This is clearly impractical in code with restrictions on the length of species names, and hence a number of simplifications have been adopted. The emphasis has been to stress the connectivity and functional groups present in a compound, so as to avoid obscuring its chemical behaviour. Inorganics

-ONO is used for most NO2 -ONO2 for most NO3 -O2NO2 for most NO4

Organics

Me Et Pr Bu -CHO

for for for for for

methyl CH3ethyl CH3CH2propyi C3H7- (i- iso; n- normal) butyl C4H9 (i, n, t) aldehydes

Halocarbons are given an identifying letter by species type, followed by a number representing the number of each type of atom in the compound (the standard notation for CFCs), omitting any finalzeros. A single letter suffix, based on the number of F atoms, may then be added to identify a particular isomer:

G.D. Carver et al./Computer Physics Communications 105 (1997) 197-215

213

F freon ( h a l o n with 1 c a r b o n atom assumed) H other h a l o n s R h a l o c a r b o n radical, c a r b o n atom activated first digit: n u m b e r o f C atoms less 1 second: n u m b e r o f H atoms p l u s 1 third: n u m b e r o f F atoms fourth: CI atoms fifth: B r atoms t i f the terminal group is - C F 3 d if the terminal group is - C F 2 A m if the terminal group is - C F A B For example, C2HF4Ct is often referred to as H C F C I 2 4 , The n a m e we use is H 1 2 4 1 t and identifies the structure more clearly whilst r e m a i n i n g sufficiently close to the c o m m o n term to be instantly recognisable,

References [ i ] L.A.M. Nejad, Ph.D. thesis (University of Manchester, UK, 1986). [2l T.J. Millar, J.M.C. Rawlings, A. Bennett, ED. Brown and S.B. Chamley, Gas phase reactions and rate coefficients for use in Astrochemistry. The UMIST ratefile, Astron. Astrophys. Suppl. Set. 87 (1991) 585. [3] A.R. Curtis and W.P. Sweetenham, FACSIMILE Release H Users Manual, AERE Report R-11771 (H.M. Stationery Office, London, 1985). [4l D_J. Lary, M.P. Chipperfield, R. Toumi, The potential impact of the reaction OH + ClO .--* HCI + 02 on polar ozone photochemistry, J. Atm. Chem. 21 (1995) 61. [5] J. Austin, On the explicit versus family solution of the fully diurnal photochemical equations of the stratosphere, J. Geophys. Res. 96 (1991) 12941. [6] D. Dent, L, Isaksen, G. Mozdzynski, G. Robinson, E Wollenweber and M. O'Keefe, The 1FS model performance measurements, Coming of age: Proceedings of the sixth ECMWF workshop on the use of parallel processors in meteorology, G.-R. Hoffman and N. Kreitz, eds. (World Scientific, London, 1995). [ 7 ] R. Atkinson, D.L. Baulch, R.A. Cox, R.E Hampson Jr, J.A. Kerr and J. Troe, Evaluated kinetic and photochemical data for atmospheric chemistry: supplement IV, Atm. Env. 26A (1992) ! 187. [8] W.B. DeMote et al., Chemical kinetics and photochemical data for use in stratospheric modelling, Evaluation number 10, Jet Propulsion Laboratory, Pasedena, Calif., Publication 92-20 (1992). 191 R.E Wayne, Chemistry of Atmospheres, 2nd ed. (Clarendon Press, Oxford, 1991) p. 103. ! 10] E. Hesstwedt,. Hov and I.S.A. lsaksen, Quasi-steady-state approximations in air pollution modelling, Int. J. Chem. Kinet. 10 (1978) 971. [ 11 ] G.D. Carver and P.A. Stott, IMPACT: an implicit time integration scheme for chemical species and families, Annales Geophysicae (1997), submitted. 112l NAG, The NAG Fortran Library Manual, D02 chapter, NAG Ltd (Wilkinson House, Jordan Hill Road, Oxford, UK, i996). [13] P.N. Brown, G.D. Byme, and A.C. Hindmarsh, VODE: A Variable Coefficient ODE Solver, SIAM J. Sci. Star. Comput. 10 (1989) 1038. [ 14] D.Z. Stockwell, J.A. Pyle, G.D. Carver, M.P. Chipperfield, K.S. Law and D.E. Shallcross, Modelling NOx emissions from lightening, J. Geophys. Res. (1997), submitted.

G.D. Can'er et aL/Computer Physics Communications 105 (1997) 197-215

214

TEST RUN OUTPUT Both tcst runs use a stripped down version o f our stratospheric chemistry scheme run at a single point in the atmosphere. The run is for 10 days and the same species are used in both runs, only their treatment varies to illustrate the flexibility in A S A D . In both cases the output to F O R T R A N channel 6 is the same and is shown below. ************************************************************ * *

ASAD ATMOSPHERIC CHEMISTRY INTEGRATION PACKAGE

* *

* *

COPYRIGHT (c) CENTRE FOR ATMOSPHERIC SCIENCE CHEMISTRY DEPARTMENT, CAMBRIDGE UNIVERSITY, UK

* *

* *

AUTHORS: GLENN CARVER ([email protected]) PAUL BROWN WEB: http://www.atm.eh.cam.ac.uk/acmsu/asad/

* * *

*

**************************~*************~******************

***

CHEMISTRY INFORMATION

***

CHEMISTRY IS TREATING ADVECTED TRACERS IN THE ORDER:

i: 0(3P) 6: NO 11: CIO 16: HCI 21:H202 26: MeO

2: O(ID) 7:NO3 12:N20 17: HOCI 22:H20

3:03 8:NO2 13:N205 18:CION02 23: HCHO

4: OH 9: CI 14:HON02 19:C0 24:MeO0

5:H02 i0:C1202 15:HO2N02 20:CH4 2 5 : MoOOH

IF THE TRACERS WERE NOT INITIALISED IN THIS ORDER THEN THE MODEL RESULTS ARE WORTHLESS STEP NO. =

1

The model will continue to write the model step number until the end o f the run. The first run uses a chch.d file which includes families. We use the IMPACT integrator in this example. The model writes two files, one containing those species in families (fams.dat), the other containing the rest o f the species (species.dat). The model outputs the species values to these files as the run proceeds. The user should check their model run against the test output files included in the package to verify the A S A D software. We show only the last line here for brevity:

G.D. Carver et aL /Computer Physics Communications 105 (1997) 197-215

fams.dat

(last line) - Family run

DAY 9.937500 H02N02 1.09E-10 H20 4.50E-06

0X 2.46E-06 HCL 2.09E-09 HCHO 6.26E-12

spec.dat

(last line) - Family run

DAY 9.937500 N02 2.00E-09

03P 4.59E-13 NO3 4.34E-13

N0X 2.56E-09 HOCL 1.87E-13 MEOOH 1.95E-14

OID 1.37E-21 CL 3.11E-15

CLOX 4.15E-12 CLON02 1.10E-09

03 2.46E-06 CLO 4.13E-12

N20 1.50E-07 C0 1.O1E-08

N205 6.15E-10 CH4 1.50E-06

HN03 8.00E-09 H202 5.27E-13

OH 9.17E-15 CL202 7.14E-15

H02 5.13E-14 MEO0 1.43E-14

NO 5.61E-I0 MEO 9.02E-19

215

The second run uses a c h c h . d file in which all the nonconstant species are to be integrated separately as if tracers in a model. For this test we must use the stiff integrator SVODE. Although we are not using families in this run, for convenience we write the same species to the same files as given for the nonstiff case: fams.dat

(last line) - Stiff run

DAY 9.937500 HO2N02 1.10E-IO H20 4.50E-06

OX 2.46E-06 HCL 2.09E-09 HCHO 6.26E-12

NOX 2.56E-09 HOCL 1.87E-13 MEOOH 1.95E-14

CLOX 4.15E-12 CLON02 1.fOE-09

N20 1.50E-07 CO I.OIE-08

N205 6.15E-I0 CH4 1.50E-06

HN03 8.00E-09 H202 5.28E-13

OH 9.18E-15 CL202 7.14E-15

H02 5.14E-14 MEO0 1.43E-14

NO 5.61E-10 MEO 9.02E-19

spec.dat (last line) - Stiff run DAY 9.937500 NO2 2.00E-09

03P 4.59E-13 NO3 4.33E-13

OlD 1.37E-21 CL 3.11E-15

03 2.46E-06 CLO 4.13E-12