EPA R4-73-012b
October 1972 Environmental Monitoring Series
USER GUIDE
TO DIFFUSION/KINETICS
: :.. : : :: .
. .
. .
. .
S .
: : ::: : : . : . : : . k ; .::
U S>
: : : . : : : :

-------
EPA-R4-73-012b
USER’S GUIDE
TO DIFFUSION/KINETICS
(DIFKIN) CODE
by
J.R. Martinez
General Research Corporation
P.O. Box 3587
Santa Barbara, California 93105
Contract No. 68-02-0336
Program Element No. 1A1009
EPA Project Officer. Ralph C. Skiarew
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND MONITORING
U.S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
October 1972

-------
This report has been reviewed by the Environmental Protection Agency and
approved for publication. Approval does not signify that the contents
necessarily reflect the views and policies of the Agency, nor does
mention of trade names or commercial products constitute endorsement
or recommendation for use.
ii

-------
CONTENTS
SECTION PAGE
1 INTRODUCTION 1
2 INPUTS 4
2.1 Reading the Data Cards 4
2.2 Problem—Control Inputs 5
2.3 Additional Inputs 8
3 OUTPUTS 17
3.1 Summary of Input Data 17
3.2 Output of Species Concentrations 23
3.3 Utility Outputs 26
4 TECHNIQUES FOR USING DIFKIN 27
4.1 Changing the Chemical Model 27
4.2 Multiple Runs 28
4.3 Date and Time Subroutines 29
5 PROGRAM MODIFICATIONS 30
6 CONTROL OF INTEGRATION STEP SIZE 34
APPENDIX A LIST OF SUBROUTINES AND THEIR USE 35
APPENDIX B LISTING OF INPUT DECK PRODUCED BY SUBROUTINE PREDAT 47
APPENDIX C COMMON BLOCKS 59
REFERENCES 61
iii

-------
ILLUSTRATIONS
NO. PAGE
1.1 Schematic of Diffusion Model for Air Pollution Simulation 1
1.2 Schematic Flow of DIFKIN Processing 2
2.1 Process for Reading Data Used by DIFKIN 4
3.la Summary of Input Data 19
3.lb Summary of Input Data: Stoichiometric Matrix of Product
Species 20
3.lc Summary of Input Data: Description of Chemical Reactions
Used in Model 21
3.ld Summary of Input Data: Boundary Conditions and Molecular
Weights 21
3.le Summary of Input Data: Initial Diffusivity Profile 22
3.2a DIFKIN Output of Initial Concentrations 24
3.2b DIFKIN Output of Computed Concentrations 25
iv

-------
TABLES
NO. PAGE
2.1 Problem—Control Inputs 6
2.2 Additional Inputs 9
5.1 Dimension Changes Required to Add More Species 31
5.2 Dimension Changes Required to Add More Reactions 32
5.3 Dimension Changes Required to Add More Space Mesh Points 33
C.l Subroutine Location of Common Blocks 59
V

-------
1 INTRODUCTION
OEis manual contains instructions for using the computer code deve-
loped by General Research Corporation to model photochemical smog. The
mathematical details of the code have been described in previous reports 1
and will not be repeated here. These instructions are thus limited to
describing the features of the code which a user must know in order to
successfully utilize the program. In addition, instructions are provided
regarding certain modifications of the code.
DIFKIN is written in FORTRAN IV compatible with the IBM 360 system.
OEe code is composed of a main program and twenty—five subprograms.
Appendix A contains a list of the subprograms and a brief description of
the function of each. OEe code operates in single—precision mode and
requires approximately 16,000 words of core storage for execution.
Figure 1.1 is a diagram of the concepts used in the model. An air
parcel moves along a path determined by the local prevailing wind speed
.— LUTANT INFLUXES AT ANY
ElEVATION (INCLUDING THE
GROUND) ARE IMPOSED BY THE
EMISSION SOURCE FUNCTIONS
Figure 1.1. Schematic of Diffusion Model for Air Pollution Simulation
SUNLIGHT IS GIVEN AS
A FUNCTION OF TIME
SPACE/TIME TRACI
_-TH ROUGH THE SOURCE
GRID IS DERIVED
FROM WIND DAT
——
TIME.DEPENDENT MIXING
AND REACTION IS COMPUTED
FOR AIR PARCEL UP TO THE
MIXING HEIGHT h
1

-------
and direction at various times during the day. The air parcel is divided
into a mesh of equally—spaced points in the vertical direction. As the
air parcel moves over a region it receives pollutant emissions (fluxes)
from the ground. The pollutants within the air parcel then undergo
chemical reactions and mixing, with the incident sunlight driving several
of the chemical reactions that occur. At each mesh point, the concentra-
tions of several species of pollutants are computed; the output of the
program consists of these concentrations as functions of time and height.
Figure 1.2 is a simplified flow chart showing DIFKIN input, process-
ing, and output. The user provides meteorological and emission data,
along with chemical rate constants and ultraviolet radiation intensity.
OEe program computes the concentration of each of the species at each
mesh point at specified time intervals along the path of the moving mass
of air.
DIFFUS
MIXING
IVITY
HEIGHT
INITIAL CONCENTRATI
AND BOUNDARY CONDITIONS
(EMISSION DATA)
ULTRAVIOLET
INTENS ITY
I SPECIES CONCENTRATION
-..- AS A FUNCTION OF HEIGHT
[ AND TIME
Figure 1.2. Schematic Flow of DIFKIN Processing
METEOROLOG I CAL
MODULE
CHEMICAL
MODULE
2

-------
It should be noted at this time that several subroutines used in
the code are specialized for the chemical model and in general must be
modified if the chemical model is changed. These subroutines are dis-
cussed in Sec. 4.1. However, changes in the numerical values of the
chemical rate constants require no program modifications. Hence the
rate constants can be changed at will simply by changing the appropriate
input cards.
OEe meteorological inputs used in the program consist of inversion
height data and the value of the diffusion coefficients. Wind speed and
direction are not inputs to DIFKIN, but are used in an external auxiliary
program which generates the proper emissions history for the specified
wind trajectory. The auxiliary program produces as output a deck of
fluxes along the trajectory, and these fluxes are the inputs to DIFKIN.
A forthcoming version of DIFKIN will combine the current auxiliary pro-
grams with the version of DIFKIN described in this manual.
3

-------
2 INPUTS
Before discussing the input parameters, we shall provide a descrip-
tion of the input process. OEis is necessary because the proper logical
unit numbers must be assigned by the job control cards in order to run
the program.
2.1 READING THE DATA CARDS
The process of reading the data cards is illustrated in Fig. 2.1.
The data cards are read by subroutine PREDAT from the card reader, which
PREDAT refers to as unit 3, e.g., READ (3,N) LIST. PREDAT reads the
data until it detects an end—of—file mark, and then produces a listing
of the data deck on the line printer, which is designated as unit 6, e.g.,
WRITE (6,N) LIST. Appendix B shows an example of the output produced
by PREDAT. The subroutine also writes the data on a scratch peripheral
memory unit such as a disk. This unit is designated by the logical unit
number 5. After writing the data on unit 5, the subroutine rewinds unit
READS FROM CARD READER
TO WHICH IT ASSIGNS
LOGICAL UNIT NUtIBER 3
PREDAT WRITES THE
DATA DECK ON UNIT 5.
UNIT 5 IS REWOUND
AFTER WRITING
Figure 2.1. Process for Reading Data Used by DIFKIN
PREDAT PRODUCES A
LISTING OF IMAGES
OF DATA CARDS ON
LINE PRINTER
4

-------
5 and the data are then ready to be read by DIFKIN, which reads from
unit 5,
2.2 PROBLEM-CONTROL INPUTS
OEe input parameters of the program will be discussed in the order
in which they are read by the code. OEe reader should refer to Appendix
B, which contains a copy of the listing produced by PREDAT, in order to
obtain a better picture of the actual structure of the data deck. Table
2.1 shows a list of the problem—control inputs along with a description
of the function of each input variable, the name of the variable used in
the code, and the format.
5

-------
TABLE 2.1
PROBLEM-CONTROL INPUTS
Card Variable Name
Variable Function Format Remarks
Number Used in Code
I T I tLE Label for each page of the (20A4) May contain any coenta user wishea
output to make. All 80 columns of the card
may be used.
2 TIMIN Initial time of run in (40X,FlO.O)
minu tea
3 TSTOP Final time of run In (40X,FlO.O)
minutes
4 DELT Initial value of integra— (40X.F1O.O) OEe value of DELT is changed by the
tion step size in minutes code during execution. OEe initial
value should be small as a general
rule. e.g., DELT < 0 01 minutes.
5 TSTMAX Step size control parameter (40X,FlO.0) See Sec. 6 for explanation of use.
6 TSTMIN Step size control parameter (40X,FlO.0) See Sec. 6 for explanation of use.
7 PRNTIN Print interval (40X,FlO.0) Concentrations are printed every
PRI4TIN minutes.
8 UPLIM Upper limit on step size (40X,Fl0.0) OEe step size cannot be greater than
UPLIM. A typical value is 0.5 mm.
9 LO%JLIM Lower limit on step size (40X,Fl0 0) OEe step size cannot be smaller than
LOWLIM. A typical value is 0.002 mm.
10 DEL l Height of one cell of (40X,F10.0) DELl . H/(N—l) where H is the mix—
space mesh in meters ing depth and N is the number of
vertical mash points including the
ground and the upper edge, where
N — NOSTAT (see card 25).
11 SUNTIN Time interval in minutes (40X,FlO.0) OEe updating is done periodically
for updating NO 2 photo— since k 1 is a function of solar
dissociation rate constant, zenith angle. May have any dummy
value if k 1 is held constant (see
card 17). A typical value is 10
minutes.
12 FIXSTP Integration step size (40X,A4) OEis is a Hollerith variable which is
fixed (YES) or variable either YES or NO. In either case the
(NO) word must start in column 41 and must
contain no blanks.
13 KOUT Print a list of all the (40X,A4) See comment for card 12. OEe list
time steps (YES) or not contains the number of the time step,
(NO) the time at which the integration
step was performed, and the step sice
used. OEe list is written on logical
unit 9, for output on the line printer
via control cards after program exe-
cution is ended.
6

-------
TABLE 2.1 (Cant.)
Number Variable Function Format Remarks
14 TI ICtJT Print at every print inter— (40X,A.4) Sea comment for card 12. See Fig. 3.2b
val (sea card 7) the number for an example of the output.
of integration steps car-
ried out between print
intervals, and tha central
processor tia spent in the
calculation (YES) or not (NO)
15 KEVERY Print the concentrations (40 1, 1.4) See comment far card 12.
for every integration step
(YES) or not (NO)
16 SPUN Output concentrations on (401,44) See comment for card 12. OEe format used
punched cards (YES) or not for punching is compatible with the fQr—
(NO) mat used by DIFKIN for reading initial
concentrations.
1 1 QRLTE NO photodiesocietion rste (40 1,44) See comment for card 12. If gRAOE NO
2 there is no updating of k . See also
constant variable (YES) or 1
not (NO) comment for card 11.
18 IFLUX Fluxes time—variable (YES) (401,44) See comment for card 12.
or constant (NO)
19 QHTINV Inversion height time— (401,44) See comment for card 12. See Sec. 2.3.7
varying (YES) or constant for expLanation.
(NO)
20 NW4STP Maximum number of integrs— (401,110) Integer variable. Program execution is
tion steps terminated if NLJMSTP is exceeded. A
message is printed to this effect to
advise user.
21 NOREAC Number of chemical (40X,IlO) Integer variable. Cannot exceed 20 with
reactions present array dimensions.
22 NOSP1 Total number of species, (401,110) Integer variable. Cannot exceed 20 with
including resctive and present array dimensions.
non—reactive species
23 NSTDY Number of reactive species (401,210) Integer variable.
in steady state which are
to be explicitly calculated
24 NTRACE Number of non—reactive (401,110) Integer variable. OEis number ia either
(tracer) species zero or one. A zero means no tracer
species is input.
25 NOSTAT Number of vertical mesh (401,110) Integer variable. Cannot exceed 10 with
points to be used present array dimensions.
7

-------
2.3 ADDITIONAL INPUTS
The remaining inputs required by the program are discussed below.
All the inputs are presented in Table 2.2, following the same format as
Table 2.1 except that the card number is now changed to a group number.
This is because the number of cards in these inputs varies, depending on
the number of species and the number of mesh points. Each of the nine
groups of inputs is discussed in the following sections.
2.3.1 Species Name and Molecular Weight
After card number 25, the name and molecular weight of each species
are input. These are input on the same card using format (40X, A4, 6X,
FlO.2). The name and molecular weight are stored in arrays SPEC and
WTNOLE, respectively. The number of cards input is NOSP1 (see Table 2.1,
card 22).
OEe order in which the species are input is important for program
execution. The program must have
. All active species first (in any order)
. NSTDY steady—state species next (in any order)
The tracer species last (if any declared by NTRACE, card 20)
The so—called active species are those whose concentrations are computed
by integration of the rate equations. This is in contrast to the steady—
state species, which are obtained algebraically.
This input step determines the order of appearance of all the spe-
cies in any subsequent inputs which are tied to the species.
2.3.2 Initial Concentrations for Each Species
The initial concentration of each species at each mesh point in
parts per hundred million (pphm) is input next. OEis is stored in array
PPHM, and PPI-LM(i,j) denotes the concentration of the ith species at the
jth mesh point. The format used is (20X, 5ElO.4) and the input for each
8

-------
TABLE 2.2
ADDITIONAL INPUTS
Group Variable Name
Variable Function Format Remarks
Number Used in Code
26 SPEC 1 WT?K)LE Species name and molecular (40X,A4,6X,FlO.2) Each card contains the name and
weight molecular weight of one species
OEus the number of cards required
is NOSP1 (see card 22). See
Sec. 2.3.1.
27 PPHN Initial concentration of (20X,5E10.4) See Sec. 2.3.2.
each species at each mesh
point in pphm
28 RATRON Rate constants for chemical (40X,F10.0) One card per reaction makes a
reactions in units consis— total of NOREAC cards in this
tent with species concen— group. See Sec. 2.3.3.
tretion, e.g., pphc 1 min’
29 REACT Description of chemical (8A4) See Sec. 2.3.4.
reactione in the model.
30 STOYK R, STOYKP Stoichiometric coefficients (2F 10.O) See Sec. 2 3.5.
for reactants and products
31 RATK1 Table of rate constants of (2 0X,E15.5) See Sec. 2.3.6 Required only
NO 2 photodissociation reac— if QRATE YES (card 17)
tion, k 1 . in minute 1
32 HT INV Table of diffusion coef— (20X, SE 1O.4) See Sec. 2.3.7 Required only
ficients as functions of if QHTINV YES (card 19).
time and altitude, in meters
squared per minute
33 FLXTIN Table of times when fluxes (40X,Fl0.0) See Sec. 2.3.8. Required only if
are updsted, in minutes IFLUX a YES (card 18).
34 FLXWAL, VALWAL, Boundary conditions for each (30X,3E15.S) See Sec. 2.3.9.
FLXDCE species: flux at ground,
concentration at ground, and
flux at upper edge,
respectively.
9

-------
species is the concentration at each mesh point from bottom to top, with
up to five mesh points to a card as indicated by the diagram below.
CONCENTRATION AT MESH POINT 1 CONCENTRATION AT MESH POINT 5
— —
1O2 P ’ O. 1’ • P I
112 • I 11011 IZI3 l4fl TJ fl2O11fl ?110173,2 1303 137 343S3 3)3III404l42 e41404434t3 SI 5 .505 373IS IOIi I 01414 6. ,
,; 343614 .414
The number of cards per species is equal to [ (NOSTAT — 1)15] + 1 , where
NOSTAT is the number of space mesh points (see Table 2.1, card 25), and
the notation [ ] denotes the integral part of the expression within
the brackets. Cards for the tracer species must not be included if
NTRACE = 0 (card 24).
OEe initial concentrations of the steady—state species need not be
input because the code automatically computes them. However, the user
must include cards in the input deck for each steady—state species. In
our sample run, four species are designated as being in steady—state: NO 3 ,
N 2 0 5 , OH , and R0 2 . Referring to Appendix B, cards 40—43 show that the
initial concentrations of these species have not been specified. OEe
initial concentrations of these species will be computed by the code from
the given initial concentrations of the five active species, NO , HC ,
NO 2 , 03 , and HNO 2 . Figure 3.2a in Sec. 3.2 shows the computed initial
concentrations of the steady—state species.
2.3.3 Rate Constants for Chemical Reactions
These are input one to a card, i.e., one
format (40X, Fl0.O). OEey are stored in
of NOREAC rate—constant data cards. OEe
is the first one In the kinetic model as
it may be assigned any position provided
consistent with the assigned order.
card per reaction, using
array RATKON. OEere is a total
NO 2 photodissociation reaction
currently structured. However,
that any subsequent inputs are
10

-------
2.3.4 Description of Chemical Reactions in the Model
OEis input is for information only and plays no role in the compu-
tation. OEe input is in Hollerith form using format (8A4), one reaction
to a card. OEe user may make any notations in the allotted space, the
first 32 columns of the card. Appendix B, cards 60 to 75, shows an example
of these data cards. OEere is a total of NOREAC cards and the symbols
are stored in array REACT. This input need not be repeated when multiple
runs are made (see Sec. 4.2 for multiple—run procedures).
2.3.5 Specification of the Chemical System (Stoichiometric Coefficients )
A set of chemical reactions is specified by a set of NOREAC stoichio—
metric equations of the form:
v .M. ki ,M , = 1, 2, ... , NOREAC
ij j ijJ
j=l j=l
where M. = chemical formula of jth species
. = stoichiometric coefficient of reactant species
13
= stoichiometric coefficient of product species
k = rate constant of ith reaction [ RATKON(1), see Sec.
2.3.3]
$ = number of chemically active species, i.e.,
s = NOSP1 — NTRACE
OEe input of the DIFKIN code consists of the stoichiometric coef-
ficients and V j . OEese coefficients are input two per card,
v 1 and , using format (2FlO.0). The ordering scheme is best ex-
plained by an example. Consider the reaction set:
11

-------
hv+NO 2 - NO+0 (1)
O+O 2 +M- O 3 +M (2)
N0+0 3 -*N0 2 +0 2 (3)
The stoichiometry input would be as follows:
Input Using Format (2F10.O )
Card Number Species Reaction Number v
1 NO 2 1 1. 0
2 2 0 0
3 3 0 1.
4 NO 1 0 1.
5 2 0 0
6 3 1. 0
7 03 1 0 0
8 2 0 1.
9 3 1. 0
And so on for the other species. OEus it should be clear that the stoi—
chiometric coefficients are grouped by species and that for each species
the stoichiometry must be specified for all reactions. Hence, for
NOREAC reactions and (NOSP1 — NTRACE) species we have
(NOREAC) x (NOSP1 — NTRACE) cards in the stoichioinetric matrix input.
We note that the order of the species—reaction groups in this input must
be consistent with the name ordering, the concentration ordering, and so
forth. Also, it can be seen that this kind of input is very general and
affords great flexibility in specifying various chemical systems. OEe
stoichiometric matrix is input only once when making multiple runs (see Sec.
4.2 for an explanation of multiple—run procedures).
12

-------
2.3.6 Input for Variable NO 2 Photodissociation Rate Constant
These data are read only if QRATE = YES (see Table 2.1, card 17);
otherwise they must not be included in the deck.
OEe format used for these cards is (20X, El5.5) and the rate con-
stants are input one to a card, one card for each time. They are stored
in array RATK1 and the maximum number of cards is 100. OEe time interval
between updates is that specified on card 11 of Table 2.1. The last card
of the table must contain a negative number to stop the code from reading.
If QRATE = YES , then the first card on the input table must con-
tain the value of the rate constant which corresponds to the initial time
of the run. OEe program updates RATKON(1) (see Sec. 2.3.3) automatically
to reflect this fact. Hence, the value of the rate constant input in
Sec. 2.3.3 need not be the correct value inasmuch as it will be updated
when RATK1 is read. OEe updating is done by a call to the subroutine
UP RATE .
2.3.7 Table of Diffusion Coefficients as Functions of Time and Altitude
These data are read only if QHTINV = YES (see Table 2.1, card 19);
otherwise they must not be included in the deck.
The format used for these cards is (20X, 5E10.4). Each card con-
tains up to five numbers arranged in columns 21—30, 31—40, 41—50, 51—60,
61—70. The diffusivities are updated as functions of time and each entry
of the input table must contain the following: the update time, the
inversion height, and the diffusivities at the various points in the mesh.
(OEe inversion height is informational only and is not used in the com-
putation, so it really need not be input. It has been included for the
sake of completeness.) The number of diffusivities is equal to NOSTAT + 1
since each diffusivity must be specified at the midpoint of each cell as
well as at the ground and at the upper edge of the space mesh. The table
is arranged as follows:
13

-------
Card 1——Columns 21—30 Update time in minutes from start of
run
Columns 31-40 Inversion height in meters
Columns 41—50, 51—60, Diffusivity at ground and at the next
61—70 two points in the mesh, in meters
squared per minute
Card 2——Columns 21—30, etc., Diffusivities at other points in the
and thereafter
mesh, 5 to a card
Thus, for each update time there must be a set of diffusivities.
If the number of mesh points is 5, for example, two cards per update time
will be needed. OEe last group of cards in the array must contain a
negative number in the position allocated for the time in order to stop
*
the code from reading. The maximum number of update times allowed is
20 and the table is stored in array HTINV.
2.3.8 Table of Flux Update Times
OEis table contains the times when the fluxes are to be updated,
in minutes from the start of the run. OEe first entry of the table must
match the initial time of the run, TIMIN (Table 2.1, card 2). OEe last
entry of the table must be equal to or greater than TSTOP (Table 2.1,
card 3). The times are stored in array FLXTIM, and the format used is
(40X,Fl0.0). OEe very last card of this group must contain a negative
number in Columns 41-50 to stop the code from reading. OEe maximum num-
ber of cards allowed is 100, under present dimensions.
2.3.9 Boundary Conditions
Two kinds of boundary conditions are allowed in the code, a flux
boundary condition and a wall boundary condition. In the former, fluxes
*
Note that this last group must have the same number of cards as the pre-
ceding ones, even though only the first carries any information.
14

-------
are specified at the ground and at the upper edge of the mesh. In the
latter, a particular value of concentration is specified at the ground.
For reactive species, only the flux condition is allowed. For the tracer
species, the user may specify either of the two kinds of boundary
conditions.
OEese inputs are stored in arrays FLXWAL, FLXDGE, and VALWAL, where
the first one denotes the flux at the ground, the second the flux at the
upper edge of the air parcel, and the third the concentration at the
ground. Normally, FLXDGE is set to zero since no mass is allowed to flow
through the upper boundary of the air parcel; this is equivalent to say-
ing that the pollutants are trapped below the inversion base.
The boundary conditions are read as needed when the updating time
has arrived. Hence, there is no need to provide storage space for the
boundary conditions at various times. The order in which the boundary
conditions are input must be consistent with the order of the update
times previously read in (see Sec. 2.3.8). OEe boundary conditions are
read and updated by subroutine IJPFLUX, the main program being designed
to call TJPFLUX at the appropriate time. For each time, one card specifies
the boundary conditions for one species in the format (30X, 3El5.5). OEe
order of the boundary conditions on each card is FLXWAL, VALWAL, FLXDGE
in Columns 31—45, 46—60, 61—75, respectively. OEe units of the flux
depend on the units used for the species concentration. In our case the
program performs the computations in mass—fraction units, and the fluxes
are thus specified as volumetric fluxes with units of meters/minute. Any
scaling that needs to be done on the boundary conditions must be done by
IJPFLUX since there is no provision for scaling in the main program. In
our case the fluxes are scaled by a factor of i0 8 in order to preserve
accuracy in the computation.
For the tracer species, the value of VALWAL is significant if it
is other than zero. A zero value of VALt.IAL indicates to DIFKIN that no
wall boundary condition is being used (see subroutine CRANK).
15

-------
Appendix B, lines 339 to 416, shows a typical set of inputs for
flux boundary conditions. In this case only NO and HC have fluxes
at the ground and no flux penetrates the upper edge of the air parcel.
Note that the fluxes on lines 339 and 340 appear also in Fig. 3.ld, scaled
by a factor of io8. OEe program “knows” that only two fluxes need to be
read because subroutine UPFLUX, which reads the fluxes, is written to
read only two fluxes. If more or fewer fluxes are input, UPFLUX must
be modified to reflect this fact.
16

-------
3 OUTPUTS
DIFKIN produces two kinds of output: line printer and punched
card. The latter is optional, being produced at the command of the user
(Table 2.1, Card 16). The printed output is shown in the accompanying
figures. The interpretation of the output is discussed below.
The main line printer output consists of the following:
1. Images of the data cards in the input deck (Appendix B)
2. Summary of input data (Figs. 3.la through 3.le)
3. Concentrations of the various species as a function of
time and height (Figs. 3.2a and 3.2b).
In addition to these three items, several utility outputs are available
for the convenience of the user. These are described in Sec. 3.3.
The punched—card output consists of the species concentrations as
functions of time and height. The punched cards are generated at the
same time the species concentrations are printed. The format of the
cards matches the format used for input, group 27 of Table 2.2.
3.1 SUMMARY OF INPUT DATA
The line printer output shown in Fig. 3.1 consists of a summary of
the input data used in the problem. At the top of the first page, Fig.
3.la, the run label appears first. This appears on every page of the
output and corresponds to the data on the first card of the input deck.
Also appearing on the first line are the central processor time (CP TIME)
in seconds, the date of the run, and the page number. Proceeding down
the page, several input quantities are listed, beginning with the ini-
tial time and ending with the lower limit fractional change used to
control the step size. On Fig. 3.la these entries are cross—referenced
to Table 2.1.
The next listing of inputs consists of the reactant and product
stoichiometric matrices, together with the rate constant of each reaction.
17

-------
The stoichiometric matrix listings shown in Figs. 3.la and 3.lb are
interpreted as follows. The reactions correspond to the rows of the
matrix and the species to the columns. Thus consider reaction 2: the
stoichiometrjc matrix of reactants, Fig. 3.la, says that only species
NO and OZON appear in reaction 2, each with unit stoichiometric coeffi-
cient. The rate constant of reaction 2 is given as .267 pphm min 1 . The
stoichiometrjx matrix of products, Fig. 3.lb, is read in the same way
as the matrix for reactants. Thus for reaction 2 the product is one
molecule of N02. The schematic reaction obtained from reading the
stoichiometrjc matrices is (l)NO + (l)OZON - (l)N02, where the quantities
in parentheses are the stoichiometric coefficients. (Note that this
schematic reaction does not balance. The true reaction is NO + 03 - NO 2
+ O2 but 02 is not carried explicitly in the computation, hence the
omission.)
Figure 3.lc shows a list of all the reactions used in this sample
run. These reactions correspond to the schematic reactions defined by
the stoichiometric matrices, of course. Thus the listing in Fig. 3.lc
can be used to interpret the stoichiometric matrices and vice versa.
Figure 3.ld lists the species in the prescribed order, together
with the boundary conditions at the initial time and the molecular weight
of each species. In the example shown, only NO and HC (hydrocarbon) have
fluxes at the ground (FLXWAL) and no flux is permitted to cross the upper
edge of the air parcel (FLXDGE = 0). (See Sec. 2.3.9 for information
regarding VALWAL.) Note that the value of FLXWAL shown for NO and HC
is equal to that shown on cards 339 and 340 of Appendix B scaled by a
8
factor of 10 .
The last listing of inputs, shown in Fig. 3.le, is a list of the
initial values of the diffusivities. These values correspond to the
quantities shown in cards 285 and 286 of Appendix B.
18

-------
OTFO IN 01111 .
INPUT noTo 100 VERTICA l DIEFtiSION PROWLER WITH PH000CHERISTRT
I IIT I AL TIME IS 0. MIN(JTE S (INPUT CARD 2)
FINAL TI l ‘ S i.S00000.l12 MINIITIS (INPUT CARD 3)
TIMF STEP 5720 IS 1.0 00 0 00—0 2 MI II IITES (INPUT CARD 4)
1 1000 10 * 1 “1511 IN0 0 0MENT I S I.I5000E.07 METEOS (INPUT CARD 10)
S l I PpfD 00 11 fl(flflN5 IS is (INPUT CARD 21)
1U000 01 50001 15 15 q (INPUT CARD 22)
u I fO Of SPECIES I STEATIT STATE IS 4 (INPUT CARD 23)
5110510 Or 01100fW 5 00(1 15 IS U irowsoy EXCEED II(INPSIT CARD 24)
5 ,100 0, n O o000ICAL ME o POISTS Is S INCLUDINO THE DR0, 1 HD AN ’) THE 0000(INPUT CARD 25)
STfP I ?E 00 010 015
1 100 0 ( li i !” 000(0 105*5 (0*1100 IS 3.40 5000—n p (INPUT CARD 5)
0W P 1 IW IT EWACTIONAL ( ‘* 1 100 IS I.0A101 0—02 (INPUT CARD 6)
C TIME • I7.I o s l o 07/75 7I P ors
or so p
——‘FACT
0. 0 1.0
0. 0 Din
1. 0 0.0
1.0 0,0
0.0 U.n 0.0
0.0 1.0 0.1
0.0 0.0 0 .0
1 .0 1.0 0.0
1.0 0.0 1.0
0. 0 0.0
Din 1.0 1.1
0.0 1.0 0.0
0.0 0.0 0.0
0.0 0 , 0 0 .0
0.0 l. A 0.0
0.0 1.0 0.0
STOICMI
0 .0
0.0
0.0
0
0.0
0.0
0,0
0 .0
0.0
1.0
flMf T oT . . .
0.0 0 .0
0.0 0. 4
0.0 0 .0 0,0
0 ,0 0 ,0
0.0 0.0 r,n
0.0 0.0 0.0
0,0 0.0
0.0 0.0 l,fl
0.0 0.0 0,0
0.0 0.0 n,n
In.. PolE rn Sin’ I
IPPM• P I’
S • 7 l iar—fl
.05700
I .n000n —n
so o o
loon .os
? 0040
i . non
in • o On
.,0 0 0 00r—n
0.01 1211 _ n
Figure 3.la. Summary of Input Data
H OOp ‘lO i ‘120 5
0 2 0 W
ANT
0.0
1.0
0.0
0.0
0 , 0
‘In
‘IT ’ ?
or or n m , ,
1 0.0
0.0
? 1.0
(1.0
1 0 .0
1,, 0
4 0.0
0 ,0
1.0
I’ l l
1 0.0
1.0
7 l i i’
0.0
, 0. 0
0.0
o
0.0
in n.o
0 , 1 1
I i n,o
0.0
I , 0. 0
0.0
I i 0 , 0
( ‘.0
IN 0 . 0
0.0
IS 3.0
‘.0
1& 0 ,0
(‘, O
0.0 0.0 0.0 0 .0 s.0000nr—n
0.’ 1 .0 0.0 0,0 A 5 ,nfln
4.0 0.0 1.0 . f l 10. 100
0. 0 0.0 1.0 0,0 10 . 000
0.0 0•0 0,0 0 , 0 I .0 n nD nr n’
0 .n 0 ,0 0 , 0 n n 3 , O O o 0 n _ n- I
19

-------
(ltS0lN 0&0L0 ‘HI,, _ CR tTN S (7.17 *tF o,,,si7 0*00 0
Pr o nt O” — — — P 0 0 II ( ‘0 S T (‘ 1 ( I 0 0 0 0 0 — . .
I (.00 0’ O 1.00 O. 1 ’0 0.00 0.00 0.00
0• 0 0
7 (,.00 0.I1O 1.1 1 0 0.00 0.00 0.00 0.0 0 ‘1.00 ‘0
0.00
1 0.00 0.00 0.00 0.00 0.00 0 00 0.00 0.00
0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00
0 0.00 0.00 1.00 0.00 0.00 0. 00 0.00 .13
0.00
A 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
O • 00
7 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00
0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00
0 0.00 0.01’ 0.00 0.00 0.00 0.00 0.00 0.00
1.00
10 1 .00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
0.00
I I 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00
0 • 00
I ? 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00
0.00
13 0 .00 0.00 1.00 0.00 0.00 1.00
0.00
I A 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0 0 0
0 • 00
(0 0.00 0.00 0.00 0.00 2.00 0.00 0 .00 0.00
0 • 0 0
( 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00
Figure 3.lb. Summary of Input Data: Stoichiometric Matrix of Product Species
20

-------
01’pI N SA’PLC O l I N — T!ME • (7.40 flAy 07/75/77 PI pE
I. PPIOTOM.N02.MO .03
7, NO ‘03 ‘$02.02
3. 0 • Me • PROP
4. ON • MC PR O?
3. *02 • NO . MOP • .Ip1RflN
4. *02 • ND ? • PAM
?. ON ‘NO •MN07
A. O W •MOZSMNO3
0. 03 .Me •R0
tO. PHOTON • MONO . ON • NO
11. NO? ‘03 •N0 IO2
I ?. ‘i03 • $02 • $205
U. Np03 • $03 • NO ?
I’. M705 • $70 •
IS. NO • $02 • $70 I PMNO2
14. MO? • PAPTICLES • PRODUCTS
Figure 3.lc. Summary of Input Data: Description of Chemical Reactions
Used in Model
AT RIN SAMPLE PUN — CR TINt • 17.04 flAyF 07/20/71 Riot
SPEe r ts PI OSAL V’L 4L rL’not Mfli MT
NO 3.07054(’OI 0. 0. 1.O0100 ,0I
‘C 0.4437309) 1 0. o.
NO? ‘. 0. 0. 4. O1OOE’0 (
020$ 0. ‘. 0. •. 00000t.O l
•N0P 0. •. 0 ,
NO) 0. 3. 0. 4. 2 0 00 0p. O I
?05 I I. 0. 0. 1. 0 1 000 0’OP
I i , 0. 0, 1 . l 00 00t. 0 (
P02 0. 0. 0 ,
Figure 3.ld. Summary of Input Data: Boundary Conditions and Molecular
Weights
21

-------
CTIFkIII SAHPL! RUN CP TIMF • 7 1 fllTF 0 , / P S / i l Djfl(
P4171*1 DIFFUSION COEFFICIENTS
0 * 1 U E 5HflWN j O E FOR ROINYS HaLFWAY REYWEEN SY*tIOUS
57 4 11DM WEIGHT IH (YFRSI Y IFFIJSI0IYY
! 0.00 30.000 0 D E
30 .0 000 00
p I IR.flT
30.000 100
230.00
10.000000
4 1AC •Qfl
10.000000
5 460.00 30.000000
Figure 3.le. Summary of Input Data: Initial Diffusivity Profile
22

-------
3.2 OUTPUT OF SPECIES CONCENTRATIONS
Figures 3.2a and 3.2b show the output of the species concentrations
produced by DIFKIN. Figure 3.2a contains the initial concentrations
which have been input by the user. These correspond to cards 35—39 of
Appendix B. Each column of the output corresponds to one species, and
each row corresponds to a given height. The first row shows a height
of zero meters and is, of course, the ground. Subsequent rows correspond
to increasing altitudes. Note that the first column contains two species,
NO and R02; these are read in the order indicated by the order of the
name of the species, i.e., first NO and then RO2. The species are printed
folded together in multiples of 8. Hence, for 17 species the first column
would contain 3 species, and so on.
Figure 3.2b shows the output of the computation after 30 minutes
of real time have elapsed. Because the print interval, PRNTIN (card 7
in Table 2.1 and Appendix B), is equal to 30 minutes the program will
print the species concentrations and the time every 30 minutes. Some
additional information is printed on the page: the value of the diffu-
sion coefficient at the ground, the value of k 1 , and the fluxes of
NO, HC, and CO, the last being the tracer species. Since we have no
tracer species in this example, the CO flux is set to zero. Note that
the value of k 1 shown in Fig. 3.2b matches the quantity shown on card
220 of Appendix B. Similarly, the fluxes correspond to those shown on
lines 349 and 350 of Appendix B scaled by io8. Clearly, the value of
and the flux are not the values of those quantities at 0700, which
corresponds to 30 minutes. This is simply because the printing is done
before updating k 1 and the fluxes. After printing, the program will
update these quantities and proceed to the next integration step.
The last line of Fig. 3.2b shows that it took 121 integration steps
to get from 0 to 30 minutes of real time (cycle No. 1), and that this
took 19.439 seconds of central processor time. This output is optional
and is controlled by TIMOUT (Table 2.1, Card 14).
23

-------
flT kTN AMP1 RuIP C YTUF I7.A ATF ( ;/; /7 P 4F
IN!TTAL PPfl;ILFS OF C0 4CENTP*T!ONS PPNM
M(TGWt —M. NO NC N02 OZOPJ NNO2 N7O
02
o.o ,.Q0000F•OO 1.9QOOo .O I. oooo •oo 4.ll000F—o7 4.n0000F.oo S.A’ 74F—n6 .4 QR r—nP 4.0 ,9?f_
4 .94404F—06
115.00 5.ooooo .oo I.lI000F.0 ).30000r.oo 5.?8non —o2 4.onooor•oo 7.p S gr—nA 5.474fl2r.nO 7.T4lA?F—o
?3O.O ?.P000CIF.0O S.OOAAOF.OC A.500O0F OI 7 5o0fl( 2 4.00000E.0O 1.O737R —n7 S.51 n2F 0M
I .?4I6R —fl5
3 45.on ,.Ioooo .oo 4. 8 0000r.OO .4000or—oI 8.12OOo( Op 4,1 OOfl0(.OO I.III.Pr l7 T.e4na5 — A I.43Iq4F—n
1 .30899 .05
460.On ?.Ifl0C( •OO 4.dOflOOE.flA P.40 000F—0I 8.L2Ooo —o2 4.000ooEoo 1.III4AF—,17 5.A4fl4 r—, I.43 1o4F—
I .30 q9 -
Figure 3.2a. DIFKIN Output of Initial Concentrations

-------
OTFKTN SAMPLE PaiN —
(P yrs, • 37 flP n AT E fl;/?S/T PanE y
Figure 3.2b. DIFKIN Output of Computed Concentrations
CONCFNTPATIOPJ P ROFILES PA RTS EP MUMDP(fl MIILTflN
TIMEc 3.00000E .0I MINUTES
1(1 • g.74432 ont.o2
NC FIrmA 33.A73300
NO? O7ON
GPO(PMO fl1FFu’ 1VTTV • 30.000000
NITRIC o10E FLuX • i 0 7tq;on
NEIONT.M. NO NC
80?
fl.nfl 1.Q1464 1•O1
?.471?I —0’
115.Ofl 4.438371*00
1.1251 6F05
230.00 1 .A340*r.00
P .5933SE—0 .
3 4 •00 1.414flF.0fl
?.QOSAI F
4A0.0fl I.4fl5?AE.flO
?.Q l8Q?F0 .
NI’MNFP
C
0
. ,
. )
: 1 :
‘ C
2 .9 3O f .01
I .14.2PE.0I
5.25 5S6E.0fl
A •bQAQ0 .00
4. O7 OQ TF •n0
OF TIME STEPS
2.e4382F .00
2. 2626SF .00
I .r.8444F.00
1 .50325r.00
1 .499A8Fs00
IN CYCLE MO.
S.O?IQ GF—0 ?
I .81537E—0l
3.PS Q7AE —01
3.5202AF .OI
3.532 10E—0I
I IS 121
N u n?
4.01 O1OF.00
3. QQ OISE •00
1. Q8215F.flfl
3 . 9 14 • 00
1. Q 8143F .ofl
F) IFEFRI NT I At
CO FLLIX C fl
Nfl 1
6.R71 tOE— OP
2.4 8 1 0 SF— f l Y
4.4A 0 mr— cal
4. 8IA32F. 07
4 •P1 ?8flF—fl l
CP TIME •
*J205
I .09727C.O7
1.3Q466r n7
4 .26PS?r o7
4 .37 142r—0
4.37 772r—117
ON
4 .O2P7AE-OA
I .08 flfl 05
7,A If 12F.n 5
2.1 PN .ii
2.1 9A 2SF—nt
N,

-------
3.3 UTILITY OUTPUTS
Several convenient utility outputs may be produced on an optional
basis. These are:
I A list of all the time steps used in the integration
(Table 2.1, card 13). This is written on logical unit 9.
A note of the number of integration steps taken during a
print interval and the central processor time used (Table
2.1, card 14; Fig. 3.2b).
A list of the species concentrations (Fig. 3.2b) after
every integration step rather than at every print interval
(Table 2.1, card 15).
In addition to these we use peripheral units such as disks to
print certain messages which assist the user in ascertaining that the
run is proceeding properly. These messages are:
. The update time and updated value of the fluxes, k 1 ,
and diffusivities. The last are written on logical unit
number 11, and the other two on logical unit number 10.
The occurrence of negative concentrations. This is
written on logical unit 8.
By means of control cards the user can cause these various messages to
output to the line printer after program execution is terminated or, if
desired, suppress the output.
26

-------
4 TECHNIQUES FOR USING DIFKIN
4.1 CHANGING THE CHEMICAL MODEL
DIFKIN provides the user with great flexibility in changing the
chemical model. OEe remark8 that follow are intended to assist the user
in changing models efficiently.
Changing the chemical model requires that new stoichiometric matrices
be input (see Sec. 2.3.5). In this regard, it is important to keep in
mind the need to be consistent in ordering the species and the reactions.
For example, suppose that the reaction NO + NO 3 - 2NO 2 were to be added
to the reaction set shown in Fig. 3.lc, and let this be reaction number
17. Note that no new species have been added to the system. The first
step that must be taken consists of changing the number of reactions from
16 to 17 (card 21, Appendix B). OEen the appropriate rate constant and
reaction description are inserted after cards 60 and 75 in Appendix B.
OEe next step is to modify the stoichiometric matrices. This requires
that a card be inserted for each species describing the stoichiometry of
the species in the new reaction. Thus for NO we have to insert a card
after card 91 in Appendix B which has a 1 in Column 1 and a zero in
Column 11, since NO is a reactant but not a product. For hydrocarbon
we have to insert a card after card 107 in Appendix B, containing zeros
since hydrocarbon appears neither as a reactant nor a product in this
reaction. And so on for the remaining species.
OEe procedure for adding a species is described below by means of
an example. For simplicity, we assume that no new reactions are added.
OEus suppose that in the present chemical model it is desired to compute
explicitly the concentration of PAN in reaction 6, R0 2 + NO 2 + PAN
(card 65, Appendix B; Fig. 3.lc). OEis adds one new species but no extra
reactions. Suppose further that PAN is assigned position number 6 on
the list of species. The first step is to increase the number of species
in the input deck from 9 to 10 (card 22, Appendix B and Table 2.1). We
27

-------
shall also assume that PAN is not in steady state, hence NSTDY (card 23,
Appendix B and Table 2.1) does not change. The name and molecular weight
of PAN must then be inserted after card 30 of Appendix B, and the initial
concentration must be input after card 39 of Appendix B. OEe next step
is to describe the stoichiometry of PAN. OEis is done by inserting 16
cards (one for each reaction) with the stoichiometric coefficients of PAN
after card 155 of Appendix B. In this example PAN is a product of reac-
tion 6 and does not appear in any other reaction. Hence the stoichiometry
of PAN consists of zeroes on all cards except for reaction 6 which has a
zero in Column 1 and a one In Column 11. Finally, the subroutines men-
tioned in the next paragraph must be modified since a new species has
been added and the position of NO 3 , N 2 0 5 , OH , and RO 2 in the species
vector has been altered.
In addition to the changes in the inputs mentioned above, certain
subroutines must be modified or changed altogether when the chemistry
is changed. OEese are: JACOB, NUBETA, NUR.ATE, STATE, and UPRATE. OEese
subroutines must also be modified if the order of the species or the
order of the reactions is changed. OEe reader is referred to Appendix A
for an explanation of the contents of these subroutines.
4.2 MULTIPLE RUNS
DIFKIN is designed to allow the user to stack several cases using
different inputs. It is presupposed, however, that the chemical model
does not change from case to case, except that the values of the rate
constants may be changed at will.
To run multiple cases it is required that a card which has the word
“MORE” punched in Columns 1—4 be inserted between data decks. OEe pro-
gram continues to read cards until either the word “END” or “MORE” is
detected. OEe former causes execution to terminate; the latter causes
the prcgram to recycle and begin reading data.
28

-------
When making multiple runs it is not necessary to input more than
once the schematic description of the reaction set (see Sec. 2.3.4) and
the stoichiometric matrices (see Sec. 2.3.5), since the basic chemical
model does not change. All the other cards in the input deck must be
input, however. For example, to run multiple cases using the sample deck
shown in Appendix B, the user must omit cards 60 through 219 in any sub-
sequent data decks, and include everything else.
4.3 DATE AND TIlIE SUBROUTINES
OEe code utilizes system routines for obtaining the date and the
machine clock time. OEese system routines are peculiar to the computer
center where the code is to be run and thus will require special handling
at each installation. We have provided dummy routines which call the sys-
tem routines in order to avoid the need to modify the main program. OEese
dummy routines are MDATE and SECOND and are described in Appendix A.
29

-------
5 PROGRAM MODIFICATIONS
OEis section contains information needed by the user in case cer-
tain basic modifications need to be carried out. The modifications con-
sist of increasing the capacity of the program to allow for more species,
more reactions, and more mesh points.
It may be recalled that the current array dimensions of the program
limit the number of species and the number of reactions to 20 and the
number of mesh points to 10. Tables 5.1, 5.2, and 5.3 show the arrays
whose dimensions must be modified if more capacity is needed for species,
reactions, and mesh points, respectively.
Additional changes that must be made in the main program are:
1. Change the 20 in the DATA statement DATA NROW/20/ to the new
maximum number of species.
2. Change the 20 in the DATA statement DATA MAXR/20/ to the new
maximum number of reactions.
3. Change the 10 in the DATA statement DATA MA)C’ISH/lO/ to the
new maximum number of space mesh points.
30

-------
TABLE 5.1
DIMENSION CHANCES REQUIRED TO ADD MORE SPECIES
*
Array Location Modification Instructions
ISTOYR COMMON/GEM! C = MNS
RATE COMMON/CHEM/ R = t INS
BETA COMMON/GEM! C = tINS
CONIN COMMON/CHEM/ R = MNS
WTMOLE COMMON/CHEM! Length = MNS
II COMNON/CHEM! Length = tINS x MNR
JJ COMMON/CHEM ! Length = MNS + 1
VALWAL COMMON/TRACER! Length = MNS
FLXWAL COMMON/TRACER! Length =
FLXDGE COMMON/TRACER! Length = 1 1145
STOYKR Main Program C = 111 15
STOYKP Main Program C = MNS
SPEC Main Program Length = 11145
CONUP Main Program R = tINS
CONEW Main Program R = 11145
PPHN Main Program R = MNS
WORK Main Program R = MNS
CZERO Main Program R = MNS
AJACOB Main Program R = MNS, C = tINS x M W
Q MainProgram RMNS,C=MNSxMNP
*
R = number of rows of a two—dimensional array
C = number of columns of a two—dimensional array
tINS = maximum number of species
MNR = maximum number of reactions
NNP = maximum number of mesh points
31

-------
TABLE 5.2
DIMENSION CHANGES REQUIRED TO ADD MORE REACTIONS
*
Array Location Modification Instructions
ISTOYR COMMON/CHEM/ R = MNR
BETA COMMON/CHEM/ R = MNR
RATEFF COMMON/CHEM/ Length = MNR
RATKON COMMON/CHEM/ Length = MNR
I I COMMON/C l iENt Length = MNR x MNS
STOYKR Main Program R = MNR
STOYKP Main Program R = MNR
REACT Main Program R = MNR
*
R = number of rows of a two—dimensional array
C = number of columns of a two—dimensional array
MNS = maximum number of species
MNR = maximum number of reactions
MNP = maximum number of mesh points
32

-------
TABLE 5.3
DIMENSION CHANGES REQUIRED TO ADD MORE SPACE MESH POINTS
*
Array Location Modification Instructions
RATE COMMON/CHEM/ C = MNP
CONIN COMNON/CHEM/ C = MNP
DFINIT COMMON/SKY! Length = MNP + 1
ZEE COMMON/SKY/ Length = MNP
TRACIN COMMON/TRACER! Length = MNP
SCALDF COMMON/TRACER! Length = MNP + 1
CONUP Main Program C = MNP
CONEW Main Program C = MNP
PPHN Main Program C = MNP
WORK Main Program C = MNP
GZERO Main Program C = MNP
AJACOB Main Program C = MNP x MNS
Q Main Program C = MNP x MNS
SCALE Main Program Length = MNP
DSCALE Main Program Length = MNP
DCOF Main Program Length = MNP
HTINV Main Program C = MNP + 3
TRIDAT Subroutine CRANK Length = (4 x MNP) — 2
* R = number of rows of a two—dimensional array
C = number of columns of a two—dimensional array
MNS = maximum number of species
MNR = maximum number of reactions
MNP maximum number of mesh points
33

-------
6 CONTROL OF INTEGRATION STEP SIZE
The parameters TSTMAX and TSTMIN control the size of the integra-
tion step (see Table 2.1, cards 5 and 6). OEe step size is controlled by
a fractional—change criterion as follows. Let the concentration of the ith
species at the jth mesh point at the nth step be denoted by c . . The test
performed uses the fractional change parameter
n+1 n n
a.. = C . — c. Ic
1J iJ ij ii
If > TSTMAX for any (i,j) combination, the step size is halved for
the next integration cycle. If TSTMIN for all i and i the
step size is doubled for the next integration cycle. For any particular
problem some experimentation is required to determine the values of TSTMAX
and TSTMIN that yield a good combination of speed and accuracy. Our exper-
ience has been that values of TSTMAX = .03 and TSTMIN = .01 work very
well for the current model.
34

-------
APPENDIX A
LIST OF SUBROUTINES AND THEIR USE
A.1 SUBROUTINE COEQ(N,A,X,FL)
Use. Solves a set of simultaneous linear equations whose matrix
is tridiagonal. In DIFKIN it is used to solve the equations which yield
the concentration of the tracer species. It is called by subroutine CRANK.
Arguments . N = number of rows in the matrix
A = matrix elements stored in a linear vector of
4N — 2 elements
X = solution vector
FL = failure index: FL = 1 if matrix is singular, and
FL = 0 otherwise
Remarks . The nonzero elements of the vector A look like
Row 1 A(l) A(2) A(3)
2 A(4) A(5) A(6) A(7)
3 A(8) A(9) A(1O) A(ll)
N—2 A(4N—12) A(4N—ll) A(4N—iO) A(4N—9)
N—i A(4N—8) A(4N—7) A(4N—6) A(4N—5)
Row N A(4N—4) A(4N—3) A(4N—2)
A.2 SUBROUTINE CRANK(N,DELT)
Use. Solves the linear diffusion equation Bc/st = fDz [ D(z c/3z] ,
which has a diffusion coefficient which is a function of z . Uses the
Crank—Nicolson method of differencing. Used with inert species only.
Arguments . N = number of mesh points
DELT = integration step size
35

-------
Remarks . Calls subroutine COEQ for solving the resulting system
of equations. The concentrations and other variables are passed through
COMMON/TRACER!.
A.3 SUBROUTINE DFFUSE(NOSTAT, DELZ)
Use. Sets initial values of diffusion coefficients when diffusivi—
ties are not time—varying. The diffusivities are stored In the vector
DFINIT, which is transferred to the subroutine through COMMON/SKY!.
Arguments . NOSTAT number of mesh points
DELZ = mesh interval
A.4 SUBROUTINE DIFCOF(SCALDF, NOSTAT, DCOF, SCALE, DSCALE)
Use. Generates scaling vectors DCOF, SCALE, DSCALE, for use in
the integration of the rate equations of the reactive species.
Arguments . SCALDF = Input vector of scaled diffusivitles, defined
in main program by SCALDF(I) = D(I)/DELZ 2 ,
where D(I) = iII diffusivity and DELZ = mesh
interval.
NOSTAT = number of mesh points
DCOF = a vector used for scaling the matrix equation
which is the result of the finite—difference
method. DCOF(I) = SCALDF(I) + SCALDF(I+l)
SCALE = vector of scale factors
SCALE(I) = SCALDF(I+l)/SCALDF(I)
DSCALE = vector of scale factors
DSCALE(I) = l/SCALDF(I)
36

-------
A.5 SUBROUTINE JACOB(A,B,R,N,W,NN,C)
Use. Computes the Jacobian matrix of the chemical rate equations
having the form dYi/dt f 1 (y 1 ,y 2 , ‘ ‘n where y is the con-
centration of the ith species, i = 1,2, . . . , N . Each element of the
Jacobian is defined by A 1 = Df /3Yj i,j = 1,2, . . . N.
Arguments . A = Jacobian matrix of size N x N
B = vector of species concentrations
R = vector of chemical rate constants
N = order of the Jacobian matrix, equal to the number
of rate equations to be integrated
= vector of molecular weights of all the chemical
species
NN = total number of species in the system, including
species in steady—state but not the tracer species
C work vector
A.6 SUBROUTINE MATADD(A,B,C,NRA,NCA)
ENTRY MATSUB(A,B,C,NRA,NCA)
Use. NATADD computes the sum of two matrices C = A + B
MATSUB computes the difference of two matrices C = A — B
Arguments . A, B input matrices
C = resultant matrix
NRA = number of rows of ‘A, B, C
NCA = number of columns of A, B, C
A.7 SUBROUTINE MATINV(A,N,NNAX,B,M,PIVOT,IPIVOT,INDEX,DETER.t)
Use. Computes the inverse of the matrix A or the solution of the
matrix equation Ax = B . Called by MATNVT.
37

-------
Arguments . A = matrix to be inverted of order N x N
N = order of matrix A
NHAX = dimension of A as declared in calling program,
NMAX>N
B = right—hand side of matrix equation Ax = 3 ,
dimensioned B(NNAX,M)
M = number of columns of B; if M < 0 only the
matrix inverse is found, without solving
matrix equations
PIVOT = work vector of size N
IPIVOT = work vector of size N
INDEX = work array dimensioned INDEX(N,2)
DETERN = determinant of A
A.8 SUBROUTINE MATMUL(A,B,C,NRA,NCARB,NCB)
Use. Computes the matrix product C = A x B.
Arguments . A = matrix of order NRA x NCARB
B = matrix of order NCARB x NCB
C=Ax B
NRA = number of rows of A
NCARB = number of columns of A and rows of B
NCR = number of columns of B
A.9 SUBROUTINE MATNVT(A,N,D)
Use. Computes the inverse of matrix A by calling MATINV.
Arguments . A = matrix of order N x N; destroyed in the inversion
process as its inverse is substituted for it
38

-------
N = order of matrix A, limited to 24
D = determinant of A
A.lO SUBROUTINE MATSCL(X,A,NR,NC)
Use. Scales the matrix A by factor X.
Arguments . X = scaling factor
A = matrix to be scaled
NR = number of rows of A
NC = number of columns of A
A.ll SUBROUTINE MDATE(IDATE)
Use. This is a dummy subroutine which calls other system routines
to obtain the date of the run. The date is printed on every page of the
output. User must modify MDATE in order to make it compatible with his
installation.
ArZument . IDATE vector which contains the date of the run
A. 12 SUBROUTINE NEWPAG(TITLE ,LSKIP ,IDATE)
Use. This is a utility subroutine used for labeling each page of
the output. The subroutine labels each page with the run label provided
by the user, the date of the run, the internal machine time, and the page
number.
Arguments . TITLE = vector of Hollerith characters which contains
the run label which user prescribed in the first
card of the input deck (see Table 2.1, card 1)
LSKIP = number of lines skipped before printing heading;
if LSKIP = 0 the page heading is printed at the
top of the page
IDATE = vector which contains the date of the run
39

-------
A.13 SUBROUTINE N1JBETA(BETA,NR)
Use. Used to modify the net stoichiometric matrix when required
by stationarity assumptions.
Arguments . BETA = net stoichiometric matrix;
BETA(i,j) = V j Vjj i = 1,2, . . . , r, j = 1,2, . . . , 8,
where r = number of reactions and s number of
species (see Sec. 2.3.5 for definition of
and v )
NR = number of rows of BETA as declared in DIMENSION
statement
A.14 SUBROUTINE NURATE(RATEFF ,RATKON ,CONUP)
Use. This routine computes any pseudo—rate constants which result
from stationarity assumptions.
Arguments . RATEFF = vector of effective rate constants which con-
tains all the pseudo—rate constants
RATKON = vector of basic rate constants which are used
to compute the pseudo—rate constants
CONIJP = vector of species concentrations
A.15 SUBROUTINE PREDAT
Use. Reads input data from card reader and produces a listing of
the image of the input deck. See Sec. 2.1 for a full explanation.
A.16 SUBROUTINE QGEN(A,B,C,N)
Use. This subroutine generates the Q matrix which enters into
the matrix equation Q(C — C 0 ) = —2G 0 . The Q matrix is stored in the
same location as the matrix A and is defined by
ll = — C A 1 = B —
40

-------
Arguments . A = matrix of order N x N, modified by the subroutine
B = constant factor added to the diagonal elements of A
C = constant factor which scales all the elements of A
N order of matrix A
A.l7 SUBROUTINE SECOND(A)
Use. This is a dummy subroutine which calls other system routines
to obtain the internal machine time. It is used by the main program for
timing the computation. It is also called by NEWPAG.
Arguments . A = internal machine time in any desired units
A.18 SUBROUTINE SOURCE(CONIJP)
Use. Computes the rate function of the chemical rate equations.
Arguments . CONUP = work vector
Remarks . Consider a set of chemical equations defined by:
v M - v M i 1,2, . . . , r
j=l j=l
where = chemical formula of the jth species
k 1 = rate constant of ith reaction
vll = stoichiometric coefficient of jth reactant in the ith
reaction
= stoichiometric coefficient of the th product in the
ith reaction
r = number of reactions
s = number of species
The rate equations for such a system are defined by:
= R 1 = (v — ujj)ki 71’ cim ; j 1,2,.. .,s
i1 m1
41

-------
where Cj = concentration of the jth species
R rate function of the jth species
Subroutine SOURCE computes R. as shown above. The quantities needed by
SOURCE are provided through COMMON/CFEEM/.
A.19 SUBROUTINE STATE(CONEW,CONUP ,NROW)
Use. Computes algebraically the concentrations of those species
assumed to be in steady state.
Arguments . CONEW = array of updated species concentrations; CONEW(i,j)
is the concentration of the ith species at the jth
mesh point
CONUP = work vector
NROW = number of rows of CONEW as declared in a
DIMENSION statement
Remarks . Other quantities needed by STATE are provided through
COMMON/CHEM/. This subroutine must be modified when the chemical model
is changed or the order of the species and reactions is altered.
A. 20 SUBROUTINE TIMEX(KSTEP ,TITLE )IDATE)
Use. This subroutine is used for timing various phases of the
computational process. It also counts the number of integration steps
during a print interval.
Arguments . KSTEP = total number of integration steps taken to date
TITLE = Hollerith vector which contains the label of
output page
IDATE = vector which contains the date of the run
Remarks . This subroutine contains an entry called ENDRUN which is
called at the end of the computation to compute and print the net time
42

-------
taken by the computation. Subroutine TIMEX calls subroutines SECOND and
NEWPAG.
A.21 SUBROUTINE TRIVRT(A,D,SCALE,WORX,NRA,NCA,NMAT,NNR,U)
Use. Inverts a tridiagonal matrix, A, whose elements also are
matrices. The lower—diagonal entries are identity matrices. The upper—
diagonal elements are scaled identity matrices.
Arguments . A = doubly—dimensioned array of size NRA x NCA
D = array of size NRA x NMAT; right—hand side of matrix
equation
SCALE = vector of scale factors of the upper—diagonal
elements of A
WORK = work array of size NRA x NRA
NRA = number of rows of A as declared in a DIMENSION
statement in the calling program
NCA = number of columns of A, NCA = NRA x NHAT
NNAT = number of columns of D
NNR = actual order of the diagonal elements of A
U = solution array of size NRA x NNAT
Remarks . Given a matrix equation of the form
AU = D
where A is block tridiagonal, this subroutine computes the solution U.
The form of A is
A 11 a 1 1 0 . 0
I A 22 2I 0
O AN...lN...l N—l 1
o 0 I
43

-------
where N = NMAT, and the c ’s correspond to the scale factors stored in
SCALE. It can be seen that there is no need to store A in full since only
the diagonal elements, Au, are important. The zeros can be excluded, of
course, and the identity matrices are known and need not be stored. Thus
the matrix A that is input into TRIVRT contains only the diagonal elements,
Au, of the full matrix. Each is of order NNR, and the first entry
of each A 1 is stored in A 1 ., where j = (i—l)NRA + 1, i 1, . . . , NNAT.
Thus the actual A looks like
A (A 11 ,A 22 AN_i , N-i’ ANN } NRA
NMAT x NRA
Similarly, U and D are stored in the form
U = [ U]!U2,••• NRA
NMAT
D = [ D 1 ,D 2 , . . . , D 1}NRA
— -- —
NNAT
and each U 1 and D is a vector of NRA components.
A.22 SUBROUTINE UPFLUX(TIME,NTRACE)
Use. This subroutine is called by the main program at time TIME
to update the fluxes.
Arguments . TIME = time at which fluxes are updated
NTRACE = flag which tells the subroutine whether or not
a tracer species is present (see Table 2.1,
card 24)
44

-------
Remarks . The fluxes are scaled by a factor of io8 for the purpose
of improving computational accuracy in the integration process. Any
required additional scaling of the fluxes must be done by this subroutine.
The subroutine prints a message on logical unit 10 which lets the user
know the update time and the value of the new fluxes.
A.23 SUBROUTINE IJPRATE(TIME ,RATNEW,RATOLD,RATEFF)
Use. This subroutine is used to update the magnitude of rate
constants which are functions of time.
Arguments . TIME = update time
RATNEW = factor used for updating the rate constant
RATOLD = vector of old rate constants to be updated
RATEFF = work vector of old rate constants
Remarks . In our photochemical smog model there are two reactions
which are directly triggered by sunlight. They are:
k 1
hv + N0 2 + NO + 0
k 2
hv + HN.O 2 - OH + NO
As sunlight intensity changes with time of day both k 1 and k 2 must be
updated. In our model k 2 is a function of k 1 and thus the variable RATNEW
is actually the new value of k 1 . Then k 2 is computed from k 1 and stored
in the appropriate location in both RATOLD and RATEFF. OEe new value of
k 1 is also stored in these two vectors.
A.24 SUBROUTINE WXMIT(NWDS,NCLS ,SF,FROM,TO)
Use. This subroutine is used to transfer scaled values of array
FROM to vector TO.
45

-------
Arguments . NWDS = number of words to be transferred; if NWDS < 0,
then the quantity SF x FROM(l,l) is transferred
to NWDS locations in TO.
NCLS = number of rows of FROM
SF = scale factor
FROM two—dimensional array
TO = vector that receives scaled elements of FROM
A.25 SUBROUTINE XMIT(N,A,B)
Use. Transmits N words from A to B.
Arguments . N = number of words to be transmitted; if N < 0,
subroutine transmits A to N locations in B
A = value of vector to be transmitted
B = vector receiving new values
46

-------
APPENDIX B
LISTING OF INPUT DECK PRODUCED BY SUBROUTINE PREDAT
47

-------
V I IAGES OF DATA—CARPS
I i i 21 31 41 51 71
CARD .. •.. .•. • • • , • • • • .... . ....
1 DIFKIM SAMPLE RUN -
P INITIAL TIME 0.
3 FINAL TIME 360.
A INITIAL TIME STEP .01
‘. UPPER LIMIT FRACTIONAL CHANGE .03
t, LOWER LIMIT FRACTIONAL CHANGE 0.01
7 PRINT INTERVAL 30.
P UPPER LIMIt AN nELT 0.5
0 LOWER LIMIT ON DELT 0.002
10 VERTICAL MESH INTERVAL IMETERSI 115.
1 11 21 31 41 51 61 71
C ARD .... ..... •... .•... .
I i TIME INTERVAL FOR UPDATING i 10.
I ? 00 WE WANT A FIXED STEP SIZE NO
13 SHALL WE OUTPUT THE TIME STEPS NO
1 4 SHALL WE OUTPUT THE EXECUTION TIMES YES
15 SMALL WE OUTPUT EVERY CYCLE NO
16 00 WE WANT P(INCHED OUTPUT NO
17 IS K) VARIABLE Y(S
I R ARE SOURCES TIME—VARIABLE YES
19 IS INVERSION HEIGHT VARIABLE YES
?0 NUMBER OF INTEGRATION STEPS 23000
1 11 21 31 41 31 61 71
CARD •.........t.s.......+.s.....s............s...... 1........t..... . . . . . .t .t1 . . .
21 NUMBED OF REACTIONS 16
22 NUMBER OF SPECIES
23 NUMRER OF SPECIES IN STEADY STATE 4
24 NUMBER OF TRACER SPECIES 10 OR P 0
25 NUMBER OF VERTICAL STATIONS 5
26 SPECIES NAME AND MOLE WEIGHT NO 30.01
27 MC 72.2 0
2$ MOp 46.01
29 OZON 4$.
30 HNO2 47.00
1 11 21 31 41 51 A l 71
CARD . .. .... .. ......... • ....•. .... • •s• s s•i • • •.•.•• ••.•• •..•• ..... .......
31 P 403 62.
32 P 4205 10$.
33 ON 17.
34 RO ? 48,1
35 Nfl PPNM 7.9 5.0 ?. ? 2.1 ?. )
36 Nt PPMM 19.9 11.1 5.0 4.8 4.B
37 NO? PPMM 1.6 1.1 0.83 O.A4 0.84
38 03 PPMM .0411 .fl2$ .0785 .0R12 .0A12
39 MM I I ? PPMM 4 , 4. 4. 4. 4.
40 NO3 PPMM - COMPUTFO 9Y CODE
48

-------
TWAPIS OF DAT*—CARflS
1 11 21 31 41 31 61 7)
CARO •
41 P4205 PPHM. COPPUTfO OY COOl
42 ON PPHM— CONPUTID BY COOl
43 P02 PPI4M. COMPUTID BY COOl
44 P4 ,1 CONSTaNTS — PIACTION NO. 1 IKI) f BtA1Nln FPON —UPPATI—
45 RIACTION NO. 2 .267
46 QIACTION NO. 3 1.0001—6
47 P *CTION NO, 4 50.
40 PIACTION NO. 5 1.01.3
49 PCACTION NO. 0 2.0
50 RIACTION NO. 7 35.
I 31 21 31 43 5 1 71
CARD
51 Ptac’iOw NO. e 36. -
52 R aCT1ON NO, 9 4.001—5
33 PIACTION NO. 10 O0?AINIO FPO$’ •UP AT(—
54 RIACTION NO. 11 5.01—5
55 REACTION NO. 12 45.
55 REACTION NO. 13 34.
67 REACTION NO, 14 60.3
5$ REACTION NO. iS 1.01—5
39 REACtION NO. 16 .001
SO PWOYON.NOp.NO .03
1 11 21 31 41 51 63 71
CAb
*1 NO .03 .NO?.O?
62 0 .NC •$P02
63 OW •HC •1R02
SO? • NO • NO? • .1250N
$3 P C? • NO? • PAN
$6 OW •NO •NNO2
67 ON • P402 S P4 1 103
6603 •HC .PO2
66 PHOTON • PIONO. OH . P40
70 NO? .03 •N03.O2
1 11 21 31 41 51 63 P1
CARD •.........•.........•.........•
71 N03 • P402 • P4205
72 N203 • P403 • NO?
73 N205 • P420 2HN03
14 NO • P402 • P420 • 2HN0?
75 NO? • PAPTICLIS . PRODUcTS
76 0. 1. BIOIN NO STOICWTONETRIC COIFFTC?ENTS
77 I. 0
7$ 0 0
79 0 0
60 1. 0
49

-------
TP A fS O DITA—CARr S
1 I I 21 31 41 SI 61 71
CARD ...... •.... ••.... .•.•••. •...•..•..... ...
81 0 0
82 1. 0
83 0 0
84 0 0
85 0 1.
R6 0. 0.
87 0. 0.
88 0. fl.
89 0. 0.
90 1. 0.
) 11 2 ) 31 41 51 61 71
CARD •.........•.........• ..• .•.•••••.....••.....,...•....•...,•...
9) 0. 0. (ND MO
92 0. 0. AtOIM HYOROCAROOM
93 0 0
94 ). 0.
95 1. 0
96 0 0
97 0 0
98 0 0
99 0 0
100 1. 0
1 11 21 31 4) 5) 61 7)
CARD
101 0 0
102 0 0
103 0. 0.
104 0. 0.
105 0. 0.
106 0. 0.
107 0. 0. ( MI) HYDROCARRON
108 1. 0. 8(0TH P402
109 0 1,
110 0 0
I II 21 31 41 5 1 61 7)
CARD
111 0 0
112 0 1.
113 1. 0
114 0 0
)1S 1. 0
116 0 0
117 0 0
118 1. 0
319 3. 0.
120 0. 1.
50

-------
1MA( !S or DATA—CARDS
1 11 21 31 41 51 61 71
CARD
121 0 0
122 1, 0•
123 1. 0. E n NO
124 0. 1. BEGIN OZONE
125 1. 0
126 0 0
127 0 0
126 0 0
129 0 0
130 0 0
I 11 21 31 41 31 61 71
CARD •.........•.........+......,..•.......,.•.........•.......,.•.........•,....,...
131 0 0
132 1. 0
133 0 0
134 1. 0.
133 0. 0.
136 0. 0.
137 0. 0,
136 0. 0.
139 0. 0. E dD OZONE
160 0. 0. BEGIN HNO2
1 11 21 31 41 51 61 71
CARD ...... . . ..•.. . ..... ..... .. . ..••. . ...... ...... • •.•••. . . .... . ..... . . ...... ..... . . .
161 0 0
142 0 0
143 0 0
144 0 0
145 0 0
146 0 1.
147 0 0
146 0 0
149 1. 0
150 0 0
1 11 21 31 41 51 61 7 !
CARD •.........•. ••.•..•••••••.......••..•.....••..,......•..
151 0. 0.
152 0. 0.
153 0. 0.
154 0. 2.
155 0. 0, (Nfl HNO?
156 0. 0. BEGIN ND3
157 0 0
158 0
159 0 0
lAO 0 0
51

-------
TMAflFS or DATA-CA nS
1 I I 21 31 41 31 61 71
CARD *1....... ..... • •••• .. ..., . •••••• •••••• ••,•,•,•••• ••••••• ,••••• ••••••
IAI 0 0
162 0 0
163 0 0
264 0 0
163 0 0
166 0 1.
167 1. 0
168 0. 1.
1A9 0. 0
170 0. 0.
1 11 21 31 41 3 2 61 71
CARD •.........•.........•,,..,..,.•.........•.......,.•.........•.........e.....,. ..
171 0. 0. (ND N03
172 0. 0. BEO!N N205
273 0. 0.
174 0. 0.
175 0. 0.
176 0. 0.
177 0. 0.
178 0. 0.
179 0.
180 0. o.
1 11 21 31 41 31 61 71
CARD •.........,.................s[[[
181 0. 0.
18? 0. 0.
183 0. 1.
184 1. 0.
183 ), 0.
186 0. 0,
187 0. 0, END N203
186 0. 0. BEGIN OH
189 0 0
190 0 0
1 1) 21 31 41 31 61 71
CARD •s..••e.i••.••.•i.s••..•ies, ,.•s.s.....,•........i•.........•.........•.........
191 1. 0
19? 0. 0.123
193 0 0
194 1. 0
193 1. 0

-------
fwjr, (5 0? DATA—CARnS
1 11 21 31 41 51 61 7 %
CARD
201 0. 0.
202 0. 0.
203 0. 0. END Oil
204 0. 0. PEA lS P02
205 0 0
206 0. 8.
20? 0. 8 ,
20 1 1. 0
289 I. 0
210 0 0
1 11 2% 31 41 51 61 71
CARD . ... ..... •4 .. .
211 0 0
212 0. 1.
213 0 0
214 0 0
215 0. 0.
216 0. 0.
217 0. 0.
211 0. 0.
219 0. 0. (Nfl PO P (ND STA TCHTDMFTPY
220 630 K ) 5.423391—02
1 I I 21 31 4 % 51 61 7%
ciQo
221 640 7.623011 .02
222 650 9.7443 2 1— 02
223 700 1. 178671—01
224 710 1.374981—01
225 720 1.563301—01
226 730 1.743621—01
227 740 1.915911—01
228 750 2.0801S 1— 01
229 800 2.236391—01
230 810 2.384591—01
1 I I 2% 31 41 51 6) 7)
CARD . •
23) 020 2.524791—01
232 830 2.657031—01
233 840 2.78136 1.0 )
234 SS 2.89783 1—01
235 900 1.00 8511—01
234 9 )0 1.107491—ni
237 920 3.193401—01
238 930 3.279621—01
239 940 3.35082 1—01
240 050 3.431111—01
53

-------
TMAr,FS .,, flAT —CAROS
ii 21 31 41 SI 61 71
C ARO ...... • •••• •••••••• •••••
241 1 (100 3,4966 1 1— 01
242 lOin 3.5 5 543 1—01
243 1020 1.6077 1 1—01
244 l )30 3.653 8 1flI
245 104(1 3.693201—01
246 IflSO 1.726701—01
247 1100 3.754221—01
246 1130 3.775901—01
249 1120 3.79183 1— ( 1I
250 1130 3.802111—01
1 11 21 31 41 51 61 71
CARD •........• ••••••••••••••...•....
251 1140 1.806791—01
P52 1150 3.805911-03
253 1200 3.799461—0%
254 1210 3.787401—03
255 1220 1.769661—01
256 130 3.746141—01
257 1240 3.736741—0%
258 1250 3.6 61311—01
259 1300 3.639711—01
260 1310 3.591811—01
1 11 21 33 41 51 63 71
CARD ••.......••. .•.........• •‘••••.••• ••.•. •. •. • ,••s.••• ,• ,••. , , . , ,•. •
261 1320 3.537461—01
262 1330 3.47653 1—01
263 1340 3.408871—01
264 1350 3.334371—03
265 1400 3.252921—01
266 3410 3,171931—01
267 1420 3.076l0(— 01
268 1430 2.972631—0)
269 1440 2.661431 .01
270 1450 2.762401—0)
I ii 21 31 41 51 61 71
CARD •.........‘....s....•.......................................... ..... ,...., ..,...
271 1500 2.615491—01
272 1510 2.48004 101
273 3520 2.337901—03
274 3530 2.386941—01
27S 1S40 2.028031—01
276 1550 1.861081—01
277 1600 1 ,696081—01
278 IALO 1.503041—01
279 1620 1.111981—01
280 16313 1. 31293 1— 01
54

-------
TMA ( 1ES OF DATA—CARDS
I I I 21 31 41 31 61 7 1
CARD ••••••••s••• ••• •s••••s•••••••,•••ses••s••••s••••s•. .. ..•s . .. .. . .. .. .... .....
281 1640 9. 05 928E. 02
28? 1650 6.91033(—02
253 1700 4.6 0302E—02
284 LAST CARD X I —10.
283 01FF. UPDATE 0. 0. 30. 30. 10.
286 30. 30. 30.
peT oi . uPDATE 130. 57.5 30. 10.
288 30. 30. 30.
59 01FF. UPDATE 225. 172.3 60. 800. 800.
290 30. 30. 30.
1 11 21 31 41 51 61 71
CARD •.
291 01FF. UPDATE 300. 257.3 600. boo. 600.
292 600. 30. 30.
293 01FF. UPDATE 373. 402.3 600. 600. 60(1.
294 600. 600. 30.
295 01FF. UPDAtE 3 90. 460. 600. 800. 600.
294 400. 600. 600.
297 LAST CARD —10.
298 2ND LAST CARD
299 FLUX TINES 430 0. Y!MF MATCHES INITIAL RUN TTMF
300 5.
1 11 21 31 41 31 61 71
CARD •.,.......•... •....•....•.....,...•........••. .. ...•.
301 640 10.
302 441 11.
303 650 20.
304 853 23.
305 700 30.
306 703 35.
307 710 40.
308 715 45.
109 710 48.
310 720 30.
1 11 21 31 41 51 61 71
CARD •.........•. .. .•. ...... ..•..... ....•. ........•. I• •I S •••• ...
311 751 81.
112 752 82.
313 000 90.
314 020 110.
315 039 129.
116 040 130.
317 855 145.
318 60 150.
319 900 150.
320 900 150.
55

-------
IMAGES OF flaTA.CARr S
1 11 21 31 41 51 61 71
CARD
321 Q10 160.
3U 920 170.
323 921 171.
324 939 Ipq.
1000 210.
326 1003 213.
327 1o2 6 236.
328 1043 2 3,
329 1100 270.
330 1112 282.
1 11 21 31 41 51 61 71
CARD •,, ••••,, ,.,.., ,. •IIt•• ...... ••... .. ...•• ...... .••.. .... •..•.. •• •. •ISSSSI
331 1113 283.
332 1130 300.
333 1143 313.
334 1200 330.
335 1202 332.
336 1212 342.
337 1230 360. TIME MATCHES FINAl. PUN TINE
33$ LAST CARD OF VECTOR OF FLUX TIMES — . HAS DUMMY NEGATIVE FLUX TIME
339 NO 630 FLUXES 3.07 086E—07
340 HC 630 3.44323E07
1 11 21 31 41 51 61 71
CARD .. . . . . ... ..... , , , . . ... . .... .. ,•,. . . . ...... . ... . . . ..... .. ... . •.... . . .. ... . .. ... . .
341 NO 35 3. 070 $61 .0T
342 MC A9 3 5.46873 1—07
343 NO 640 3.070461—07
344 MC 640 5. 4978? —07
344 NO 641 3.4146Th—a?
344 MC 441 6.446931—07
347 NO 650 3.4146 8 1o1
34$ M C 630 6.3144 5 107
349 NO 633 1 .074951—07
358 MC 433 3.347331—0?
1 11 21 31 41 31 61 71
CARD
351 NO 700 1.724691—07
332 MC 700 4.540371—01
353 NO 703 1.726691—0?
334 MC 705 4.532121—07
353 NO 710 1.72669!07
356 NC 710 4.52 15 3107
357 P lO 7 15 1.726691—07
354 MC 715 4.413091—07
339 MO 714 2.014911—07
340 MC 718 5.030241—01
56

-------
YMAGrS O DATA—CARFIS
1 11 21 31 41 51 61 71
CARD •.. ...... .•.... .....•..... • ...e.. ...••• ,•...• ....S... •••••• I••. •• I•I•I ••IS•I ••
361 NO 720 2.oLe9 lE—07
362 MC 720 S.0279i!—07
363 N f l 751 l.7?664 !07
364 ‘ IC 75) 4.50 665E— 07
345 NO 752 4.82533E— 07
366 IIC 752 6.7511fl 07
36? NO MOO 4.359?4 —O1
364 MC 600 6.16496E07
344 NO $20 4.35924E—07
370 MC $20 8,o9$7B —07
) i i 2) 3 ) 4) 51 61 71
dM0 •........ .•... • ••.•••. ...... ..•.. ...... .•...• .....•.. ...... .•.. .. .....••• ......
371 NO 639 3.92 453t—01
372 MC $39 6.6O554 —O7
373 NO 640 3.9!953E—07
374 MC $ 60 6.80012C—0T
375 NO 65 6 3.Q2953 —07
374 MC 655 6.7370 6 —07
377 NO $60 1.37077F07
376 MC 660 3.1)9062—07
379 NO 900 1.125532—07
360 MC 900 2.617232—07
1 11 2) 31 41 5) 61 71
CAMO •.........•.........•.........•.........•.........•...................•.........
341 NO 900 1.60145207
342 MC 900 2 ,762192—07
343 NO 910 1.6014 5 1—07
344 M C 910 2.767962—07
345 NO 920 1.6014 5107
366 MC 920 2.769942—07
347 NO 921 7.69501 208
364 MC 921 1.663462—07
369 NO 939 1.9256)2—07
390 MC 939 3.372122—07
1 1 ) 2) 31 4) 51 63 71
CARO •t..I .I.It+SS .. .eSSS•I .SS•S .II• .. .SI.gSI•StS •ItI•I•I.S •S. . . •• .SII . .Se .•IISSII . .I
191 NO 1000 1.945492—07
392 MC 1000 3. 779662—01
393 NO 1003 9.669762—08
394 MC 1003 2.0 6196F—07
395 NO 1026 1.017972—07
196 MC 1026 2.4037)2—07
397 NO 1043 1.75 130F—07
396 MC 1043 2.566722—07
399 NO 1100 1,q 630$F—Ol
400 MC 1100 2.q0 972207
57

-------
THMWS (W DATA—CAQflS
11 21 31 41 51 61 7 1
CARD ...... ...• • ..•..... • ..••...... ...•....
MU NO I I I ? 1.05725E—07
402 NC 1112 2.32002 —07
403 NO 1113 7.84901E—08
404 MC 1113 1.88121E—07
405 NO 1130 7.A4q0 1r—o8
406 NC 1130 1.R8576 07
407 1143 9.3qS9o —oq
408 MC 1143 2.67180 —O8
NC) 1200
410 MC 1200 2.A7180 —08
1 11 21 31 41 51 Al 71
CARD ... . .......s..................................•.... ........
411 NO 1202 3.34153(—08
412 1202 1.0 5322E07
413 NO 1212 6. 11234E 08
414 MC 1212 1.S2622 07
415 MO 1230 6.11234E—08
$16 MC 1?30 ND FLUX 1.5262 V—07
$17 EPID NOTE TO USER. FOR STACKING RUNS WRITE MORE INSTEAD OF END IN COLS. 1.4
58

-------
APPENDIX C
COMMON BLOCKS
Three COMMON blocks are used in DIFKIN. The locations of these
blocks are shown in Table C.l. All three blocks are located in the main
program also.
TABLE C.l
*
SUBROUTINE LOCATION OF COMMON BLOCKS
COMMON Block Name Subroutine Location
CHEM SOURCE, STATE
SKY DFFUSE
TRACER CRANK, UPFLUX
*
All COMMON blocks are also found in the main program.
59

-------
REFERENCES
1. A. Q. Eschenroeder and J. R. Martinez, Further Development of the
Photochemica]. Smog Model for the Los Anzeles Basin , General Research
Corporation CR—l—191, March 1971.
61

-------
7. Author(s Martinez
9.
5. Rep. rt Dote
October 1972
6.
8.
Performing Organization
No.CR2 273
Rept.
10.
Project/Task/Work
Unit
No
11. Contract/Grant No.
68-02-0336
BIBLIOGRAPHIC DATA r . Report No 2.
SHEET I EPA-R4-73-012b
3. Recipient’s Accession No.
4. title and ubtLtie
User’s Guide to Diffusion/Kinetics (DIFKIN) Code
Performing Organization Name and Address
General Research Corporation
‘ . o. Box 3587
Santa Barbara, California 93105
12. Sponsoring Organization Name and Address
ENVIRONMENTAL PROTECTION AGENCY
Research Triangle Park, North Carolina 27711
13. Type of Report & Period
Covered
14.
15. Supplementary Notes
16. Abstracts
This manual is intended for users of the GRC Diffusion/Kinetics (DIFKIN) code
for simulating photochemical smog. The general structure and operational
capabilities of the code are described. Detailed instructions for program use
are provided to assist prospective users in operating the code.
17. Key Words and Document Analysis 17o l)cscripcors
Air pollution
Simulation
Photochemi cal reacti ons
Smog
Diffusion
Mathemati cal models
Computer programs
Manuals
Input output routines
17b. Identifiers/Open-Ended Terms
Diffusion/Kinetics (DI FKIN) code
Atmospheric modeling
DIFKIN
17c. COSATI Field/Group
18. Availability Statement 19. Security CIa,s (This 21. No. of Pages
Report)
Unlimited UNCLASSIFIED 71
20. Security Class (This 22. Price
Page
UNCLASSIFIED
ORM NT )s-35 IRCV. 3-72)
USCOMM-OC 14052 .P72

-------
INSTRUCTIONS FOR COMPLETING FORM NTIS-35 (10-70) (Bibliographic Data Sheet based on COSATI
Guidelines to Format Standards for Scientific and Technical Reports Prepared by or for IIe Federal Government,
P 8-180 600).
1. Report Number Each individually bound report shall carry a unique alphanumeri’e dcsignation sclected by the performing
organization or providcd by the sponsoring organization. Use uppercase letters and Arabic numerals only. Examples
FASE13-NS-87 and FAA-RD-68-09.
2 Leave blank
3. Recipient’s Accession Number. Restrvcd for use by each rc port recipient
4. Title ond Subtitle I irk should indit iti . Ic irly and brolly ihc suhjtci covi ragc of h rcport, and be displayed promi-
nently Set subtitlt , if us, d, in smalL r typt or otherwist subordinate it to main titli . When a rt port is prepared in more
than one olumt , repeat t he ptirnars t it It ad! olumt number md inc lode subtitle for t In spec i t o . volume
5. Report Dote. I ach it port shall carry a dan indicating at It ast month and year. Indicate the basis on which it was seiected
(c g , date of issue, dait. of approval, d itt of prepazaiioii
6 Performing Organi lotion Code. Least. blink
7. Author(s) Give name(s) in conventional orII r (e g , john Ft Dot, or J Robert l)oi ) l.ist author’s affiliation if it differs
from the performing organizaiion
8. Performing Organi zotion Report Number liiseri if ptrforiniiig org.initac ion wishes to issign this number.
9 Performing Organi zation Home and Address. (‘let name, sircet, i ii y, state, and tip itch Lisi no morc than two levels of
an org ini.’at ional hit ran hy Display ihi name of the organization txaetly as it should ap it in Government indexes such
as USGRDR-I
10. Project Task/Work Unit Number I st thi project task an,! sork otto numbers undtr wliii h th report was prepared.
I Contract/Grant Number Ins ri c itnir itt or grant numbe r undt r wli h report w is pre pire i i
12. Sponsoring Agency Name and Address Inc hid 1 tip odi
13 Type of Report and Period Covered Itidic in inn rim, Ictial, i ic , and, if applicable, dati s overcd
14- Sponsoring Agency Code Lease blink
15 Supplementary Notes I nr c r infiitttiar ion not tic Iuded , lst when but useful, such ,i Pri p ired in coupe ration with
1 r.inslation of Pri sented at c ,,nfc 0 nit’ of I i i Ii, published iii 1 ,up r—i It Supplements . .
16 Abstract Inc Iudi_ a brief (201) words or Ii ss) actual sutnu ltiry of the most significant information cont iincd in the report.
If the report contains a significant bihliugr iph ) or literaturt survty, mention it here
17. Key Words and Document Anolysis (a) Descriptors. ‘ n Itt t from th u. I hcsaurus of I tigittet ring and St it ntifiu. I crms the
proper authorited tcrms that idtnt if> the major concept of the research and are suffie ii. nily spec fit and precise to be used
as indes entries for cataloging
(b). Identifiers and Open- Ended Terms Ls ud niifit-rs for project names, code name—. equipment designators, etc Use
open-ended terms written in des’. ripior form for those subje is for which no descriptor x ists
(c). COSATI Field/Group. Field md Group issignmt nts art to be taken from the 1965 CUSA 1 I Subject Category List
Since tht majority of dot uments art mu lt id use i pit nary in nature, the primary Field /Group ass ignme nt(s) will be the spec ific
discipline, area of human endeavor, or t>pe of physical object The application(s) will he cross-referenced with secondary
Field/Group ass ignmt nts that will follow the primary posting(s)
18. Distribution Statement Denote re Icasahulit > to the public or limitation for reasons other than security for example “Re-
lease unlimited’’. Cite any as lilabilit > iii ilte public, vath address ano price.
19 & 20. Security Classificotion. I)o not submit classified reports to the National feehnical
21. Number of Pages. Insert tht total numhc r of pages, including this one and unnumbered pages, but excluding distribution
list, if any.
22. Price. Insert the price set hy the Nat iotial 1 eehnieal Information Service or the (‘overnmcnt Printing Office, if known.
FORM NTiS3S RE’d 3 721 uscbMM-oc , 4 95a-P12

-------