EPA-600/3-80-028b
February 1980
MODELING OF SIMULATED PHOTOCHEMICAL
SMOG WITH KINETIC_MECHANISMS
Volume 2. CHEMK: A Computer Modeling
Scheme for Chemical Kinetics
by
G. Z. Whitten
H. Hogo
Systems Applications, Incorporated
950 Northgate Drive
San Rafael, California 94903
Contract No. 68-02-2428
Project Officer
Marcia C. Dodge
Atmospheric Chemistry and Physics Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
This report has been reviewed by the Environmental Sciences Research
Laboratory, U.S. Environmental Protection Agency, and approved for publica-
tion. Approval does not signify that the contents necessarily reflect the
views and policies of the U.S. Environmental Protection Agency, nor does
mention of trade names or commercial products constitute endorsement or
recommendation for use.
ii
-------
ABSTRACT
Mechanisms that describe the formation of photochemical smog are devel-
oped using a computer modeling technique directed toward the simulation of
data collected in two smog chambers: an indoor chamber and a dual outdoor
chamber. The results of simulating 164 different experiemnts are presented
in Vol. I. Individual compounds for which specific experiments were simu-
lated and mechanisms developed include the following: formaldehyde, acet-
aldehyde, ethylene, propylene, butane, and toluene. Experiments in both
chambers were simulated for all these compounds. The mechanisms reported
describe the decay of the precursor organic compound, formation and decay of
secondary organics, conversion of nitrogen oxides, formation of nitrates,
and the appearance and decay of ozone. Special emphasis is given to the
chemistry of toluene. Also included is a study of a generalized smog-based
or carbon-bond mechanism developed in a previous study. Vol. II contains
the user's manual and coding for a chemical kinetics computer program, .CHEMK.
This report was submitted to the U.S. Environmental Protection Agency
in fulfillment of Contract No. 68-02-2428 by Systems Applications, Incorpor-
ated. This report covers the period 23 August 1978 to 23 August 1979, and
work was completed as of 23 August 1979.
iii
-------
CONTENTS
Abstract iii
Figures vi
Tables vi
1. Introducti on 1
2. CHEMK Input 3
Basic input for a new problem 3
Program options 11
References 22
Appendices
A. Program listing A-l
B. Input data—CHEMK format B-l
C. Input data—EPA format C-l
D. Output data D-l
E. Example use of variable photolysis and print options E-l
-------
FIGURES
Number Page
1 Data input flow chart 4
TABLES
Number Page
1 Input Card Description 16
-------
SECTION 1
INTRODUCTION
Calculation of the time-dependent concentration profiles of a number
of reacting species in a complex reaction mechanism is central to the
science of chemical kinetics. Such calculations allow the validation of
kinetic schemes and reaction rate constants through the comparison of
predicted concentrations with experimental data. The need to rely on
approximation methods such as steady-state assumptions, quantum yields,
and a "rate-determining step" in making these calculations adds uncertainty
to predicted concentrations and could possibly conceal inaccuracies
in the kinetic mechanism if the experimental data were highly scattered.
These approximation methods are also time-consuming in that the user must
formulate them and evaluate their impact on the predictions.
Work by Whitten (1974) and Overton (1976) led to the development of
fast and efficient computer programs that make complex calculations common-
place. These computer programs integrate the system of ordinary differen-
tial equations (DDEs) associated with a kinetic mechanism, given a set of
initial conditions. The method used to integrate the system of ODEs was
developed by Gear (1971) and modified by Hindmarsh (1974). Initial work
by Whitten (1974) led to the development of a computer program that could
handle 50 chemical species and 200 chemical reactions. That program,
CHEMK, used the modified Gear routines documented by Hindmarsh (1974).
Since that time Spellmann and Hindmarsh (1975) have modified the Gear
routines to handle sparse matrices. The Jacobian, or the matrix of partial
differential equations, associated with a given chemical kinetic mechanism
is usually sparse, meaning that it contains more zero elements than non-
zero elements. Thus, the use of a routine designed specifically for sparse
matrices would increase the efficiency of the computer program CHEMK in
terms of central processing (CP) time. The work of Spellman and Hindmarsh
1
-------
and our experience with the limitations of the first versions of CHEMK led
us to modify CHEMK.
In this report, we present a new version of CHEMK that can handle
chemical kinetic mechanisms containing up to 89 species and up to 200
reactions. The new version can also handle photolysis rate constants and
temperatures that vary with time, and reaction schemes with stoichiometric
coefficients. The Gear routine used is that of Spellmann and Hindmarsh
referred to above.
The version of CHEMK documented in this report requires approximately
44000 decimal words to load on a Control Data Corporation (CDC) 7600 com-
puter. The efficiency of the new program has been increased about a factor
of 3 from the original versions of CHEMK. Rather than present a detailed
description of the program and its various routines, the emphasis in this
report has been placed on providing a description of program inputs and
options available to the user. A complete program listing is provided in
Appendix A. Appendices B and C give the input data for an example run in
the CHEMK and the EPA format, respectively. Appendix D gives the output
of that example run in EPA format, and Appendix E provides an example of
variable photolysis and print options.
-------
SECTION 2
CHEMK INPUT
A brief description of each data input card and its argument list is
given below. Both initial program specification and modified operation
are treated in this presentation. The data input flow is illustrated
schematically in Figure 1.
BASIC INPUT FOR A NEW PROBLEM
General Data Card
The arguments that can be entered on this card are:
> The last number of the reaction sequence.
> Choice of reaction format.
> Print option to print only instantaneous concentrations
at specified intervals.
> Dilution factor.
> Number of time intervals in which the photolysis rate
constants are changing.
> The time interval of each photolysis number.
> The number of photolysis reactions.
> The number of temperature values.
> The time interval of each temperature value.
> Option to plot the photolysis values as a function of time.
The last number of the reaction sequence should correspond to the number
on the last card in the reaction scheme, even though it is not the highest
reaction number. The next parameter (choice of reaction format) is to
enable users of EPASIM (Overton, 1976) and users of the first versions of
-------
*«x
i
i £
is*
< o w
IP!
i «•= =
T
i
t.
in
°
-
*-
Is5
a
5
u
S
£ I
g
i/i
512
U U
—• hU
IE G»
M W»
Is
sili
5
3
-------
CHEMK to use the current version of CHEMK. A more detailed discussion of
this option is presented in the section on kinetic mechanism cards.
An option has been added to the CHEMK program that will control the
printing of the species' net rates of formation and the net rates of reac-
tion. When this option is not used, only species concentrations are printed
at specified time intervals.
To activate this option so that species net rates of formation and
net rates of reaction will be printed at the same specified time intervals,
the user enters any alphanumeric character anywhere in columns 11 through
14 of the General Data Card.
The dilution factor represents a constant dilution rate for every
species such that:
dC.
where
C. = concentration of species i
D = dilution factor.
In many experiments all species are being diluted by a common first-order
factor. If this factor is entered in the third field of the General Data
Card (the one that identifies the last reaction and on all subsequent
appearances of such a card), then this option will be activated. If one
species (e.g., 02) is part of the dilution gas, a first order "flow out"
(e.g., 02> the only reactant with no products) reaction with a negative
rate constant equal to the dilution factor will cancel the dilution for
such a species.
-------
The next three parameters (number of time intervals in which the
photolysis rate constants are changing, the time interval of each photoly-
sis number, and number of photolysis reactions) are used to input the con-
ditions of the time-varying photolysis rate constants. If a blank or zero
is entered for the number of time intervals, then the time interval and
number of photolysis reactions should not be entered. If the number of
time intervals is nonzero then the first card after the General Data Card
should contain the vector of reaction numbers that identify the reactions
with varying photolysis rate constants. Up to seven numbers may be placed
on each card.
The next card(s) will contain the initial photolysis value and the
values at the end of each time interval. If, for example, the user inputs
5 as the number of time intervals and 60 (minutes) as the time interval,
then six photolysis values should be entered on the card (one for the
initial value and one for every 60 minutes thereafter). Up to seven
photolysis values may be entered on one card. The total time in which the
photolysis rate constants are changing in the above example is 300 minutes.
If the simulation is carried for a longer period of time (e.g., 600 minutes)
then all photolysis rate constants are set to zero after 300 minutes. The
maximum number of time intervals is 97 for each simulation. If more values
are needed the MORE option may be used (see section on MORE option).
If the number of temperature values is nonzero, the next set of cards
is the vector of temperature values (initial plus the end of each time inter-
val specified on the General Data Card). The format for the temperature
values is similar to that for the photolysis values. All temperature values
must be in units of degrees Kelvin (K). All activation energies entered on
each reaction card must also be in units of degrees Kelvin. The temperature
at each specified printout time is printed with the rest of the species
under the column TEMP.
-------
If the user wishes to obtain a plot of the photolysis values as a
function of time, then a nonzero value is declared in columns 56 to 60.
A line printer plot is produced at the end of the species concentration
versus time profiles. This option should be activated only when the plot
option is declared.
2. Kinetic Mechanism Cards
Each chemical reaction in the mechanism must be entered individually.
Entry requires (1) that the reaction be numbered, (2) that its components
be specified in alphanumeric form, and (3) that the rate constant at 298K be spec-
ified. If the reaction is temperature-dependent, then the activation energy
(in degrees Kelvin) can also be specified. A maximum of 200 reactions is
allowed by the program.
The user has a choice of two different input formats for the reactions.
To accommodate users of EPASIM (Overton, 1976) CHEMK allows nonunity stoichi-
ometric coefficients for the products of the reactions. The reaction input
format described in EPASIM and summarized in Table 1 is available for users
of EPASIM and the current version of CHEMK. A maximum of three reactants and
three products (with stoichiometric coefficients) are allowed under the
EPASIM format. To use the EPASIM format, the user must enter a nonzero value
in the second field of the General Data Card (see preceding subsection).
Reaction identification numbers are not used under the EPASIM format
as described by Overton (1976). Therefore, we introduced this parameter in
columns 17, 18, and 24 of the reaction card. The first field represents the
reaction number from 1 through 99. If the user has more than 99 reactions,
a modulus factor must be specified in column 24 to convert the reaction num-
ber to a number greater than 100. For example, a reaction identification
number of 13 is given on the reaction card in columns 17 and 18 while column
24 is left blank. A reaction identification number of 123 is given on the
-------
reaction card first by input of 23 in columns 17 and 18 and 1 in column 24.
For reaction 200, columns 17 and 18 are left blank (representing zero) and
a 2 is entered in column 24.
For users of earlier versions of CHEMK and new users a different reaction
input format is available. A maximum of three reactants and four products
(without stoichiometric coefficients) is allowed under this format.
For the reacting species, any combinations of letters and numbers can
be used to identify the individual entries with the exception that the letter
M is exclusively reserved to designate the total pressure and TEMP is reserved
for the temperature.
To specify a "flow in" reaction, that is, the continuous addition of a
particular species, a reaction with only products can be used where the reac-
tion rate k is equal to the flow rate. A reaction with a reactant but no
products implies a first-order "flow out" process.
When preparing the data, one must ensure that the rate constant is in
units consistent with the order of the reaction as specified in the reaction
mechanism. Instead of input of a rate constant for each photolysis reaction,
a photolysis ratio is entered when variable photolysis is used (see the section
on the General Data Card). This ratio is the ratio of the actual photolysis
rate constant to the photolysis values entered on the cardfs) following
the General Data Card. For example, if the user inputs photolysis values,
then each of the photolysis reactions will have a photolysis constant
determined at each instant in time from the product of the interpolated
photolysis value and the ratio entered on the reaction card. All photol-
ysis values must be declared on the card after the General Data Card,
as discussed earlier. Therefore, for example, if the photolysis values
represent.N02 photolysis rate constants, then the N02 photolysis reaction
identification number must be declared and a value of 1.0 entered for
the rate constant. It is also wise to keep in mind the following points:
8
-------
> H02 and HOO will be treated as different species.
> 2ND will be treated as a new species and not as two NOs
unless the EPASIM format is used.
> The maximum number of different species is 89.
Title Card
Each simulation is titled by entering the desired notation in columns 1
through 28 of this card. Anything may be entered in the first four columns
of this card with the following exceptions:
> All blank: This card will terminate the program.
> CONT: This is an option indicating that concentrations and
species from a previous run will be used as initial con-
ditions for a new run.
> MORE: This choice indicates that new reactions or new photo-
lysis and temperature values are to be added.
> PLOT: This entry activates the plotting routine.
Output Specification Card
The first two items on this card specify the number of reactants with
nonzero initial concentrations and the number of integration steps between
print cycles. Normally, the program suppresses the printing of the concen-
trations of all species as a function of the integration steps. However,
the user may key the output to the frequency of the integration steps by
entering a negative number of steps between print-outs in the appropriate
slot on the data card. Thus, if the user enters a value of -5, the pro-
gram will print the instantaneous concentrations every 5 iteration steps.
If the user only wants concentrations printed at specified intervals, then
this option should not be activated (i.e. leave this slot blank). To
obtain output at specific times, the user enters the initial print time
and then either a fixed time increment or a fixed multiple to be used
between subsequent printings. For example, a user may specify an initial
-------
print time of 10 minutes and set the time increment at 10 minutes. In
this case, output will be printed at 10 minutes, 20 minutes, 30 min-
utes, and so on. Or the user may specify an initial print time of 10 min-
utes and then set the print multiple at 2, in which case output is pro-
vided at 10 minutes, 20 minutes, 40 minutes, 80 minutes, and so on.
Initial Reactant Identification Card(s)
The vector containing the alphabetic designation of all reactants with
a nonzero initial concentration is given on this card(s).
Initial Reactant Concentration Card(s)
The initial concentrations of all reactants, in appropriate units, are
provided on this card in the same sequence as species are entered on the
Initial Reactant Identification Card(s).
Time Specification Card
The simulation starting time and ending time are the first two arguments
on this card. Their time unit should be identical to that used to specify
rate constants. In addition, the absolute temperature and the error toler-
ance for the integration step may be entered on this card. Default values
_2
of these parameters are 298K and 10 , respectively.
After each simulation the user may choose the following options:
> Start a new simulation with the same mechanism by first
inserting a new title card.
> Terminate the program by inserting a blank card.
> Provide a new or changed reaction sequence and start a
new calculation using the MORE option.
10
-------
- Continue the calculation using the current results but
with some parameter changed using the CONT option.
Plot the computed output using the PLOT option.
All of these possibilities are discussed in the following section.
PROGRAM OPTIONS
More Than One Computation
Rather than enter a terminating blank card, insert a new title card
followed by a new set of initial condition data. The new data will provide
computational results using the previously specified reaction scheme.
To Continue with Previously Calculated Species
If the new title card has CONT in its first four columns, the species
from the previous calculation and their concentrations will be used as initial
conditions for a problem in which either new reactions are entered or the
temperature or error tolerance is changed. The next card must in this case
be the time specification card, having the starting and ending times, error
tolerance, and temperature.
New or Changed Reactions
To introduce either new or changed reactions into the kinetic mechanism,
enter a card with the word MORE punched in columes 1 through 4 followed by
i
the set of cards used for specifying the changes in the kinetic mechanism.
This set includes a first card similar to the General Data Card. Any
parameter controlled on the original General Data Card can be changed
on this new version. The difference between this card and the original
is that all blank fields will default to the parameters controlled on
the original card or to the present values changed by previous uses of
11
-------
the MORE option. For example the dilution factor may be the only change
desired. To accomplish this a MORE card followed by a card with only the
new dilution factor in columns 21 through 30 would be used. The next
card in this example might then be a CONT to continue the previous cal-
culation or a new title to start an entirely new calculation with the
same mechanism modified only by the new dilution factor.
The most frequent utilization of this option is typically for addi-
tions or alterations to the original mechanism. To alter any reaction a
new reaction card is entered with the same identification number as the
original. The new card must contain the complete reaction in its new
form since it will replace the original. New reactions must have identi-
fication numbers not previously used in the original mechanism. In any
case the identification number of the last reaction card to be read in at
this time must appear on the new General Data Card in columns 1 to 3.
A new variable photolysis vector, a new variable temperature vector,
or a new set of photolysis reactions controlled by the variable photolysis
vector can also be introduced using this option. If a new photolysis reac-
tion is to be added to the original set that varies with the time dependent
photolysis vector, then the entire new set, including both the originals
and the new additions, must be read in.
Plot Option Card
At the end of the computational cycle, a card with PLOT in its first
four columns must be entered. This signals the program that the plot option
is to be used and plot data are to follow.
Plot Title Card
If a plot has been selected by entering the word PLOT on the above card,
then the user must supply a plot title card, which may have storage option
12
-------
entry and an option to print the species concentrations at certain time steps.
The storage option entry determines whether the plotted output will be stored
on a tape file for later access. If the output is to be stored, the user need
only insert any nonzero number in the option argument space.
The data is stored on a file titled TAPE 7 in unformatted form. The
following parameters are written to the file TAPE 7:
> NTIT—The plot title.
> NAME(L)--The Lth species being plotted.
> CLOW—The lowest concentration entered on the plot
specification card.
> CHIGH—The highest concentration value entered on the
plot specification card.
> TLOW—The lowest time value entered on the plot
specification card.
> THIGH--The highest time value entered on the plot
specification card.
> NDAT--Number of observed data points.
> NPNT--Number of calculated data points.
> (TIME(J),J=1,NDAT)—The vector of time values of the
observed data.
> (DATA(J),J=1,NDAT)—The vector of observed
concentration values.
> (SAVTIM(J),J=l,NPNT)--The vector of time values of
the calculated data.
> (SAVCON(L,J),J=1,NPNT)~The vector of calculated concentrations.
The twelve parameters are written with a single WRITE statement and should
be read in the same manner.
If the user is interested in the actual concentrations of the species
being plotted, then a nonzero value entered for this option will activate the
13
-------
printing of the actual concentrations at the times closest to every one-
eightieth of the total simulation time. Thus, if the total simulation time
is 400 minutes then calculated concentrations at 5-minute intervals are
printed.
Plot Specification Card
The plotting routine can handle up to three species on one plot. This
Is done by specifying the number of species to be plotted on the plot speci-
fication card. The information entered on this card is: The number of spe-
cies to be plotted, the concentration labels for the vertical axis (five labels
must be entered), the lowest and highest values of the concentrations, and
the lowest and highest values of the integration time.
Species Information Card
This card contains the name of the species to be plotted (the species
identification code), the number of data points to be plotted from experi-
mental data, and the plot output symbol. If the number of experimental data
points is zero, the next card is not entered. A species identification card
must be entered for each species to be plotted.
Experimental Plot Data Cards
A unique feature of the PLOT program is its ability to overlay
experimental data over a computer concentration profile. To do this, pairs
of time and concentration values must be entered for each species. Four
pairs of these values appear on each card, and the number of entries corres-
ponds to the value specified on the previous data card.
End of Plot Data Card
A blank card must be entered to signal the end of the plot data. The
next card can start a new simulation.
14
-------
End of Program Card
A blank card entered instead of a new title or option card ends the
program.
Table 1 provides an individual description of each data card, its
argument list, format code, and column entries.
Calculation of the Round-Off Error
The variable UROUND in block data ALPHA! should be set to the round-
off error associated with each computer system. Currently, UROUND is set
to the round-off error of 7.5 x 10 associated with a CDC 7600 computer
system. To reset UROUND, the card (W.I6) in block data ALPHA1 should be
changed:
DATA UROUND/user's round-off error/ (card W.16)
UROUND is calculated from the number of significant digits (N) used
for the mantissa of a floating point constant:
UROUND = 2~N
For example, on the UNIVAC 1110 computer, each word contains 36-BITs of
-27
which 27 are used for the mantissa. Thus, 2 is equal to approximately
7.5 x 10"9. This is the value required for UROUND in CHEMK.
15
-------
TABLE 1. INPUT CARD DESCRIPTION
Card
Column Format
Description
General Data Card
Next card(s)
Next card(s)
1-3 13 Last number of reaction sequence
6-10 15 Choice of reaction input format:
Nonzero=EPASIM, default=CHEMK
format
11-14 A4 Print option for net rates of
reaction and species rate of
formation: default=printing is
suppressed
15-20 — Not read
21-30 F10.0 Dilution factor
31-35 15 Number of photolysis time units
36-40 15 Time interval of each photolysis
time unit
Number of photolysis reactions
Number of temperature time units
Time interval of each temperature
time unit
Option to plot photolysis values
as a function of time
61-80 — Not read
1-10 7110 If there are variable photolysis
11-20 rate constants, the reaction
21-30 identification vector identify-
ing number of each photolysis
reaction is read here
•
61-70
71-80 — Not read
1-10 7F10.0 The initial photolysis and the
11-20 value at the end of each time
21-30 unit. If photolysis occurs for
31-40 N units of time, then there
should be N+l photolysis values
41-45 15
46-50 15
51-55 15
56-60
61-70
71-80
Not read
16
-------
Table 1 (Continued)
Card
Column Format
Description
Next card(s)
Kinetic Mechanism Card
(EPASIM format)
1-10
11-20
21-30
31-40
61-70
71-80
1-4,
7-11,
13-16,
17-18
19-23
24
25-28
32-36
37-40
41-45
7F10.0 If variable temperature is
used then the temperature
values are read here. If
there are N units of time
in which the temperature
changes, then there should
N+l temperature values
A4
A4
A4
12
Not read
The names of the reactants
(up to three are allowed)
The identification number
of the reaction (number
may be between 0 and 99; 0
only if next field is
nonzero)
F5.0 The stoichiometric coeffi-
cient of the first product;
a blank implies a default
value of 1
II Modulus factor for reaction
identification numbers
greater than 100 (values are
0 or blank, 1, and 2); if a
value of 1 is entered, the
reaction identification num-
ber will be between 100 and
199 depending on the value
in columns 17 and 18
A4 The name of the first
product
F5.0 The stoichiometric coeffi-
cient of the second product
A4 The name of the second
product
F5.0 The stoichiometric coeffi-
cient of the third product
17
-------
Table 1 (Continued)
Card
Column Format
Description
Kinetic Mechanism Card
(CHEMK format)*
Title Card
Output Specification
Card
49-52 A4 The name of the third
product
55-64 F10.0 The rate constant at 298K
or, for photolysis reac-
tions, the ratio of the
photolysis value to reac-
tion number 1
66-72 F7.0 The activation energy in
degrees Kelvin (K)
13 The reaction identifica-
tion number
A4 The names of the reactants
(up to three allowed)
A4 The names of the products
(up to four allowed)
1- 3
6- 9
11-14
16-19
21-24
26-29
31-34
36-39
41-50
51-60
1-72
F10.0 The rate constant at 298K
or the photolysis ratio
for photolysis reactions
F10.0 The activation energy in
degrees Kelvin (K)
18A4 The first four columns must
not all be blank or contain
one of the following: MORE,
CONT, PLOT
1-3 13 Number of reactants with
initial nonzero concentra-
tions (maximum is 89)
21-30 110 Number of integration steps
between print.cycles;
default is 104 to suppress
printing keyed to step size.
31-40 F10.0 Specific time of initial
print (independent of step
size)
18
-------
TABLE 1 (Continued)
Card
Column Format
Description
Initial Reactant Card(s)
41-50 F10.0 Time interval
51-60 F10.0 Time multiple (note that
at least one of the above
two fields must be blank)
1-4, A4 Alphanumeric identification
11-14 of the reactants with ini-
21-24 tial nonzero concentrations
seven species per card
Initial Reactant Con-
centration Card(s)
Time Specification Card
61-64
65-80
1-10
11-20
21-30
61-70
71-80
1-10
11-20
21-30
Not read
F10.0 Initial concentrations of
reactants. Concentrations
must be given in the same
sequence as the species are
given on Initial Reactant
Card(s)
Not read
F10.0 Simulation starting time.
Time must be in the same
units as the rate constants
F10.0 Simulation ending time.
Time must be in the same
units as the rate constants
F10.0 Temperature (constant
throughout the simulation
if variable temperature
profile is not used;
default=298K)
31-40 F10.0 Error tolerance; default=10
41-80 -- Not read
-2
19
-------
TABLE 1 (Continued)
Card
Col limn Format
Description
Options Card
Plot Option Cards:
Plot Title
Plot Specification
Card
1- 4 A4 Plot initiated if PLOT is
entered. Additional kin-
etic steps are processed
if MORE is entered. A
new starting time is used
if CONT is entered. Pro-
gram terminates if blank
card is inserted.
5-80 — Not read
1-12 3A4 Title of simulation
13-54 — Not read
55-60 Option to print actual con-
centrations of all plotted
species
61-65 15 Data storage option. If
nonzero, plot data is saved
on a file named TAPE?
66-80 — Not read
1-5 — Not read
6-7 12 Number of species to be
plotted; up to three species
may be placed on one plot
8-15 — Not read
16-19 A4 Concentration labels:
21-24 Five arguments must be
26-29 entered
31-34
36-39
41-50 F10.0 Lowest concentration value
51-60 F10.0 Highest concentration value
61-70 F10.0 Lowest time value
71-80 F10.0 Highest time value
20
-------
TABLE 1 (Concluded)
Card
Column Format
Description
Species Identification
Card (one card for each
species)
1-4 A4 Name of species to be
plotted
6- 7 12 Number of data points to
be plotted from experimen-
tal data. Value must be
between 1 and 80. If set
at zero, no experimental data
are read
Experimental Plot Data
End of Plot Data
End of Program
8
9
10-80
1-10
11-20
21-30
31-40
61-70,
71-80
1- 4
5-80
1- 4
5-80
—
Al
—
F10.0
F10.0
F10.0
A4
—
A4
—
Not read
Symbol for plotted output
Not read
Time
Concentration
(A total of four sets
time and concentration
placed on each card)
Blank card
Not read
Blank card
Not read
of
are
* Default reaction format
21
-------
REFERENCES
Gear, C. W. (1971), Numerical Initial Value Problems In Ordinary Differen-
tial Equations (Prentice-Hall, Englewood Cliffs, New Jersey).
Hindmarsh, A. C. (1974), "GEAR: Ordinary Differential Equation System
Solver," Report UCID-30001, Rev. 3, Lawrence Livermore Laboratory,
Livermore, California.
Overton, J. H. (1976), "Users Guide to EPASIM--A Chemical Kinetics Simu-
lation Program," TN-262-1643, Northrop Services, Incorporated,
Huntsville, Alabama.
Spellmann, 0. W., and A. C. Hindmarsh (1975), "GEARS: Solution of Ordinary
Differential Equations Having a Sparse Jacobian Matrix," Report UCID-
30116, Lawrence Livermore Laboratory, Livermore, California.
Whitten, G. Z. (1974), "Rate Constant Evaluations Using a New Computer
Modeling Scheme," 167th National Meeting, American Chemical Society.
22
-------
APPENDIX A
PROGRAM LISTING
A-l
-------
A 3
C****AAA******A*************AAA******************AAAAAAAAAAA***********A /\ 4
C A 5
C CHEM< IS A FORTRAN PROGRAM WHICH, WHEN GIVEN A PREDETERMINED A 6
C KINETIC MECHANISM, COMPUTES THE CONCENTRATIONS OF THE VARIOUS A 7
C REACTANTS IN TIME. IT WAS WRITTEN BY GARY Z. WHITTEN OF SYSTEMS A 8
C APPLICATIONS INCORPORATED IN SAN RAFAEL, CA AND INCORPORATES A 9
C THE GEAR INTEGRATION PACKAGE OF A. HINDMARSH OF LAWRENCE LIVERMORE A 10
C LABORATORIES IN LIVERMORE, CA. THE PRINTER-PLOT ROUTINE WAS A 11
C PROVIDED BY DAVID C. WHITNEY OF SAI AND THE PROGRAM WAS CONVERTED A 12
C ANSI FORTRAN BY JIM MEYER OF SAI. THE COMPUTER CODE THAT FOLLOWS A 13
C IS EXECUTABLE ON ANY IBM OR CDC MACHINE WITH A FORTRAN IV COMPILER A 14
C A 15
C TWO SEPARATE SETS OF DATA ARE NEEDED TO EXECUTE THE PROGRAM. ONE A 16
C SET CONTROLS THE COMPUTATIONAL PART OF THE PROGRAM AND IS A 17
C NECESSARY TO ESTABLISH THE INTEGRATION ROUTINE. THE SECOND SET OF A 18
C CONTROLS THE PLOTTER OUTPUT AND PROVIDES THE PARAMETERS NECESSARY A 19
C SPECIFY THE OUTPUT FORMAT. THE READER IS REFERED TO THE PROGRAM A 20
C FOR FURTHER INFORMATION. A 21
C A 22
C GARY Z. WHITTEN A 23
C SYSTEMS APPLICATIONS INCORPORATED A 24
C 950 NORTHGATE DRIVE A 25
C SAN RAFAEL, CALIFORNIA 94903 A 26
C (415) 472-4011 A 27
C A 28
C**'********A*****A************AAAAA***'****A***AAAA**AAA**A*******'**l**AAA ^ £9
C A 30
COMMON /DAW NR,KR(200,7),A(200),S(200),ITITLE(7),TEM?,ERR,START, A 31
lSTOPP,PC,NP,SIG(91)fIP(91),ITYPE(200),R(200),BK,SG,DILUT A 32
COMMON /NAMES/ SPECIS(91),REACT(91),NS A 33
COMMON /FRPLOT/ NIT(3),SAVCON(90,80),SAVTIM(80),JGRID(89,40),NT A 34
COMMON /ALPHA/ IGO(4),IBLANK,MBLANK,OINTER A 35
COMMON /APLOT/ JVERT(52,2),JBLANK,JSTAR,JPLUS,JBAR A 36
COMMON /GEAR!/ T,GUESS,HMIN,HMAX,EPS1,UROUNO,NC,MF1,KFLAG1,JSTART A 37
COMMON /GEAR2/ YMAX(IOO) A 38
COMMON /INOUT/ IN,IOUT,ITAPE A 39
COMMON /HEAT/ CV,Q(200),SC(200,7),ISC(200,3),ITEM3 A 40
COMMON /SPARSE/ IA(91),JA(1000) A 41
COMMON /STORE/ AST(35),IPL(7),TEMEND,NTEMP,TMI,NPHOT,PHI,IL,NFRST, A 42
1IPH(30),QM(100),PM(100),PSTOP A 43
COMMON /PHOTR/ IPP,RDAT(80),RTIM(80),RR1,IN10 A 44
DIMENSION C(91), IRS(200,7), KRS(7), SC1(7) A 45
INTEGER SPECIS,REACT,PHI,TMI,AST A 46
C A 47
C ALPHAMERIC DATA ARE PREASSIGNED IN BLOCK DATA ALPHA! A 48
C A 49
DATA NPLOT/4HSAVE/,IBLK/1H / A 50
DATA IN2,IN3,IN5,IN6,IN7,IN8,IN9/7*1/ A 51
DATA XN4,XN6,XN9/3*1.0/ A 52
A-2
-------
C A 53
C INITIAL PARAMETERS A 54
C A 55
IN=5 A 56
IOUT=6 A 57
ITAPE=7 A 58
C A 59
C CLEAR PATTERN MATRIX AND SET THE FIRST ELEMENTS A 60
C ALSO SET STOICHIOMETRIC COEFFICIENTS EQUAL TO 1 A 61
C A 62
DO 5 0=1,200 A 63
DO 5 K=l,7 A 64
IRS(J,K)=IBLANK A 65
SC(J,K)=1. A 66
5 KR(J,K)=0 A 67
DO 10 J=l,200 A 68
A(J)=0. A 69
R(J)=0. A 70
10 S(J)=0. A 71
NR=0 A 72
NS=0 A 73
C A 74
C NX = LAST NUMBER OF THE REACTION SEQUENCE A 75
C NPLOT = PLOT OPTION A 76
C DILUT= DILUTION FACTOR A 77
C A 78
READ (IN.185) NX,NEPA,NPLIT,DILUT,NPHOT,PHI,IL,NTEMP,TMI,IPP A 79
GO TO 20 A 80
15 NS=NS-1 A 81
READ (IN.185) NX,IN2,IN3,XN4,IN5,IN6,IN7,IN8,IN9,IN10 A 82
IF (IABS(IN2).NE.O) NEPA=IN2 A 83
NPLIT=IN3 A 84
IF (ABS(XN4).NE.O) DILUT=XN4 A 85
IF (IABS(IN5).NE.O) NPHOT=IN5 A 86
IF (IABS(IN6).NE.O) PHI=IN6 A 87
IF (IABS(IN7).NE.O) IL=IN7 A 88
IF (IABS(IN8).NE.O) NTEMP=IN8 A 89
IF (IABS(IN9).NE.O) TMI=IN9 A 90
IF (lABS(INlO).NE.O) IPP=IN10 A 91
20 IF (IL.LE.O.OR.IN7.LE.O) 60 TO 25 A 92
READ (IN.200) (IPH(I),I=1,IL) A 93
25 IF (NPHOT.LE.O.OR.IN5.LE.O) 60 TO 30 A 94
NPHOT=NPHOT+1 A 95
READ (IN.190) (PM(I),I=1,NPHOT) A 96
PM(NPHOT+1)=2.0*PM(NPHOT)-PM(NPHOT-1) A 97
PM(NPHOT+2)=3.0*PM(NPHOT)-2.0*PM(NPHOT-1) A 98
PSTOP=FLOAT((NPHOT-1)*PHI) A 99
30 IF (NTEMP.LE.O.OR.IN8.LE.O) 60 TO 35 A 100
NTEMP=NTEM?+1 A 101
READ (IN,190) (QM(I),I=1,NTEMP) A 102
A-3
-------
qM(NTEMs-H)*2.*QM(NTEMJ)-qM(NTef>-l) A 103
QM(NTEMP+2)=3.0*QM(NTEMP)-2.0*qM(NTEMP-l) A 104
TEMENO=FLOAT((NTEW>-1)*TMI) A 105
35 IF (NX.GT.O) WRITE (IOUT.130) A 106
IF (NX.LE.O) NS-NS+1 A 107
IF (NX.LE.O) 60 TO 80 A 108
C A 109
C REACTION INPUT DATA A 110
C A 111
40 IF (NEPA.LE.O) READ (IN.135) J,(IRS(J,I),I=1,7),A(J),S(J) A 112
IF (NEPA.LE.O) 60 TO 55 A 113
READ (IN.220) (KRS(I),I=1,3),J,SC1(1),JJ,KRS(4),(SC1(LL-3),KRS(LL) A 114
1,LL=5,6),RTE,ENERGY A 115
IF (JJ.GT.O) J=JJ*100+J A 116
DO 45 11-1,6 A 117
45 IRS(J,II)=KRS(II) A 118
A(J)=RTE A 119
S(J)=£NERGY A 120
DO 50 11=4,6 A 121
50 SC(J,II)=SCl(II-3) A 122
55 DO 60 1=1,7 A 123
60 IF (SC(J,I),LE.O.) SC(J,I)-1. A 124
DO 75 K=l,7 A 125
NL=K*5 A 126
NF=NL-4 A 127
DO 65 LK=NF,NL A 128
65 AST(LK)=IBLK A 129
IPL(K)=IBLK A 130
IF (K.EQ.1.0R.K.EQ.4) 60 TO 70 A 131
IF (IRS(J.K).NE.IBLANK) IPL(K-1)=JPLUS A 132
70 IF (IRS(J,K).NE.IBLANK) CALL VALU (SC(J,K),K,NF,NL) A 133
75 CONTINUE A 134
WRITE (IOUT.140) J,(AST(I),I=1,5),IRS(,J,1),IPL(1),(AST(I),I=6,10), A 135
1IRS(J,2),IPL(2),(AST(I),1*11,15),IRS(J,3),IPL(3),(AST(I),1=16,20), A 136
2IRS(J,4).IPL(4)f(AST(I),1*21,25),IRS(J,5),IPL(5),(AST(I),1=26,30), A 137
3IRS(J,6),IPL(6),(AST(I),1=31,35),IRS(J,7),IPL(7),A(J),S(J) A 138
KR(J,1)=100 A 139
IF (J.GT.NR) NR
-------
IF (ITITLE(l).EQ.IGO(l)) GO TO 15 A 153
IF (ITITLE(l).EQ.160(2)) GO TO 115 A 154
IF (ITITLE(l).EQ.160(3)) GO TO 85 A 155
IF (ITITLE(l).EQ.IBLANK) STOP A 156
GO TO 90 A 157
C A 158
C CALL PLOT A 159
C A 160
85 READ (IN,180) (NIT(I),I»1,3),IDT,KALCM> A 161
NT=NT-2 A 162
CALL PLOT (NIT,NT,NS,SPECIS,SAVTIM,SAVCON,IDT,KALCM?,JGRID) A 163
GO TO 80 A 164
C A 165
C OPTIONS CARD A 166
C A 167
90 READ (IN,155) N,NPRNT,TPRNT,TSTEP,TFACT A 168
C A 169
C TIME STEP SKIP OPTION A 170
C A 171
IF (NPRNT*NPRNT.EQ.O) NPRNT=100000 A 172
C A 173
C TIME STEP LENGTH OPTION A 174
C A 175
IF (TPRNT*TPRNT.EQ.O.) TPRNT=1.E10 A 176
C A 177
C CONCENTRATION OF SPECIES INITIALLY PRESENT A 178
C A 179
READ (IN.125) (REACT(I),I=1,N) A 180
READ (IN,190) (C(I),I=1,N) A 181
C A 182
C STARTING AND ENDING INTEGRATION TIMES A 183
C A 184
READ (IN.175) START,STOPP.TEW,ERR A 185
IF (START*START.EQ.O.) START=0. A 186
C A 187
C SPECIFICATION OF THE TEMPERATURE IF UNSPECIFIED IN INPUT A 188
C A 189
IF (TEMP.LE.O.) TEMP=298. A 190
C A 191
C SPECIFICATION OF THE ERROR BOUND IF UNSPECIFIED IN INPUT A 192
C A 193
IF (ERR*ERR.EQ.O.) ERR=l.E-2 A 194
C A 195
C OUTPUT OF THE INITIAL CONDITIONS A 196
C A 197
WRITE (IOUT.145) (ITITIE(I),1=1,7),(REACT(I),I=1,N) A 198
95 WRITE (IOUT,150) (C(I),I=1,N) A 199
C A 200
C COMPUTE THE NET RATES OF REACTION A 201
C A 202
A-5
-------
IF (NTEMP.6T.O) TEM?=QM(1) A 203
IF (ITITLE(1).EQ.IGO(2).AND.NTEMP.GT.O) TEMP=TEMOLD A 204
CALL RATES (C,N) A 205
C A 206
C OUTPUT OF THE TEMPERATURE AND ERROR BOUND A 207
C A 208
WRITE (IOUT.160) TEM>,ERR A 209
C A 210
C OUTPUT OF THE DILUTION FACTOR A 211
C A 212
IF (DILUT*DILUT.EQ.O.) DILUTX). A 213
IF (DILUT.NE.O.) WRITE (IOUT.170) DILUT A 214
C A 215
C SET LIMITS FOR TIMED OUTPUTS A 216
C A 217
IF (TSTEP.NE.O.) YMAX(1)=TSTEP A 218
IF (TSTEP.NE.O.) PC=TSTEP A 219
IF (TSTEP*TSTEP.EQ.O.) YMAX(1)=1.E10 A 220
IF (TSTEP*TSTEP.EQ.O.) PC*1.E10 A 221
IF (TFACT.NE.O.) YMAX(2)»TFACT A 222
IF (TFACT*TFACT.EQ.O.) YMAX(2)=1. A 223
C A 224
C RATE CONSTANTS A 225
C A 226
WRITE (IOUT,165) (R(I),M,NR) A 227
C A 228
C WRITE INITIAL CONDITIONS OF THE CELL A 229
C A 230
IF (NPHOT.LE.O) 60 TO 105 A 231
WRITE (IOUT.215) (IPH(I),I=1,IL) A 232
DO 100 1=1,IL A 233
J=IPH(I) A 234
YMAX(10*I)=R(J) A 235
100 CONTINUE A 236
WRITE (IOUT,225) (YMAX(10+I),M,IL) A 237
NPT=NPHOT A 238
WRITE (IOUT.205) PHI,(PM(I),I«1,NPT) A 239
105 IF (NTEMP.LE.O) 60 TO 110 A 240
NPT=NTEM? A 241
WRITE (IOUT,210) TMI,(OM(I),M,NPT) A 242
110 6UESS=1.E-10 A 243
T=START A 244
IF (NPLOT.EQ.IBLANK) YMAX(3)=0. A 245
IF (NPLOT.NE.IBLANK) YMAX(3)=1. A 246
IF (NPLIT.NE.IBLANK) YMAX(4)=1. A 247
IF (NPLIT.EQ.IBLANK) VMAX(4)=0. A 248
C A 249
C INITIALIZE PARAMETERS A 250
C A 251
CALL YFIX (NS,TPRNT,C,NPRNT,INX) A 252
A-6
-------
INX=4 A 253
TEMOLD=TEMP A 254
GO TO 80 A 255
C A 256
C CONTINUATION OF DATA A 257
C A 258
115 READ (IN.175) START,STOPP,TEMP,ERR A 259
IF (TEMP.LE.O.) TEMP=298. A 260
IF (TEMOLD.NE.298..AND.TEMOLD.NE.O.) TENP=TEMOLD A 261
IF (ERR*ERR.EQ.O.) ERR=l.E-2 A 262
DO 120 1=1,NS A 263
120 REACT(I)=SPECIS(I) A 264
C A 265
C INITIAL CONDITIONS A 266
C A 267
WRITE (IOUT.145) (ITITLE(I),I=1,7),(REACT(I),I=1,NS) A 268
N=NS A 269
GO TO 95 A 270
C A 271
C A 272
125 FORMAT (7(A4,6X)) A 273
130 FORMAT (1H1.14H THE REACTIONS,86X,13HRATE CONSTANT.2X.15HACT. ENER A 274
1GY (K)) A 275
135 FORMAT (I3,2X,7(A4,1X),2F10.0) A 276
140 FORMAT (/,2X,I3,2X,3(5A1,1X,A4,2X,A1),1H=,2X,4(5A1,1X,A4,2X,A1),1P A 277
1E11.3.2X.E13.3) A 278
145 FORMAT (1H1,30X,7A4,//,23H INITIAL CONCENTRATION ,//(10X,10(4X,A4, A 279
14X))) A 280
150 FORMAT (/(8X.1P10E12.3)) A 281
155 FORMAT (I3,17X,I6,4X,3F10.0) A 282
160 FORMAT (/,34H THE TEMPERATURE OF THE CELL WAS =,1PE9.2,26H AND THE A 283
1 ERROR TOLERANCE =,E9.2) < A 284
165 FORMAT (/,29H THE RATE CONSTANTS USED WERE,/,(/,8X,lP10E12.3)) A 285
170 FORMAT (/,30H THE OVERALL DILUTION RATE WAS.1PE9.2) A 286
175 FORMAT (8F10.0) A 287
180 FORMAT (3A4,43X,2I5) A 288
185 FORMAT (I3,2X,I5,A4,6X,F10.2,6I5) A 289
190 FORMAT (7F10.3) A 290
195 FORMAT (7A4) A 291
200 FORMAT (7110) A 292
205 FORMAT (/44HOTHE N02 PHOTOLYSIS NUMBERS IN INTERVALS OF ,I3,15H TI A. 293
1ME UNITS ARE,/,(1HO,1P10E13.3)) A 294
210 FORMAT (/34HOTHE TEMPERATURES IN INTERVALS OF ,I3,15H TIME UNITS A A 295
1RE/(1HO,1P10E13.3)) A 296
215 FORMAT (/29HOTHE PHOTOLYSIS REACTIONS ARE,/(1HO,9I13)) A 297
220 FORMAT (2(A4,2X)>A4,I2,F5.1,I1,A4,2X,2(F5.0,1X,A4,2X),F10.2,1X,F7. A 298
12) A 299
225 FORMAT (57HOTHE PHOTOLYSIS RATIOS OF EACH OF THE REACTIONS ABOVE A A 300
1RE/(1HO,4X,1P9E13.3)) A 301
END A 302-
A-7
-------
SUBROUTINE YFIX ( NO, TLAST,C,NQ, INDEX)
C
C
C THIS IS THE NORMAL OUTPUT ROUTINE
C
C
C
C
C
COMMON /SPARSE/ IA(91),JA(1000)
COMMON /DATA/ NR,KR(200,7),A(200),S(200),ITITLE(7),TEMP,ERR,START,
1STOPP,PC,NP,SIG(91) ,IP(91) ,ITYPE(200) ,R(200) ,BK,SG,DILUT
COMMON /NAMES/ SPECIS(91),REACT(91),NS
COMMON /FRPLOT/ NIT(3),SAVCON(90,80),SAVTIM(80),JGRID(89,40),NT
COMMON /6EAR1/ D,H,DUM(4),IDUM(4)
COMMON /6EAR2/ YMAX(IOO)
COMMON /INOUT/ INPP,IOUT,ITAPE
COMMON /HEAT/ CV,Q(200),SC(200,7),ISC(200,3),ITEMP
COMMON /STORE/ AST(35),IPL(7),TEMEND,NTEMP,TMI,NPHOT,PHI,IL,NFRST,
1IPH(30) ,QM(100) .PMdOO) .PSTOP
COMMON /PHOTR/ IPP,RDAT(80),RTIM(80),RR1,IN10
DIMENSION C(91), RT(200)
INTEGER SPECIS,REACT,TMI,PHI,AST
DATA IGO/4HCONT/
DATA NT1/1/
DATA START1/0.0/
N=NS-1
NS2*NS-2
NFRST=1
TIME INTERVALS
IF (ITITLE(l).NE.IGO) GO TO 30
IF (NT1.EQ.O) GO TO 30
TDC=(STOPP-START1)/80.
TD=START1+TDC
J=0
NTSV=NT
DO 25 1=1,80
5 J*J+1
IF (J.GT.NTSV) GO TO 20
IF (SAVTIM( J) .GE.TD) GO TO 10
GO TO 5
10 SAVTIM(I)=SAVTIM(J)
DO 15 K-1.NS
15 SAVCON(K,I)=SAVCON(K,J)
NT=I
GOTO 25
20 IF (I.EQ.l) NT=1
SAVTIM(I)=TD
25 TD=TD+TDC
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
A-8
-------
60 TO 40 8 51
30 TDC=(STOPP-START)/80. B 52
TD=START B 53
DO 35 1=1,80 B 54
TD=TD+TDC B 55
35 SAVTIMdHD B 56
START1=START B 57
NT=1 B 58
40 NC=0 B 59
C B 60
C INITIALIZE PARAMETERS B 61
C B 62
NHS=0 B 63
MS=IABS(NQ) B 64
TPRNT=TLAST+START B 65
TSTEP=PC B 66
TFACT=YMAX(2) B 67
IF (TFACT.GT.l.) TSTEP=0. B 68
MT=50-((NS-1)/10)*3-(NR-1)/10 B 69
NN=NS+1-MINO(NS/10,1) B 70
IF (YMAX(3).EQ.O.) NPLOT=0 B 71
IF (YMAX(3).NE.O.) NPLOT=1 B 72
IF (YMAX(4).EQ.O.) NOPT=0 B 73
IF (YMAX(4).NE.O.) NOPT=1 B 74
IF (NOPT.EQ.O) MT=50-((NS-1)/10)*2 B 75
NTP=0 B 76
NT1=NPLOT B 77
IF (NPLOT.EQ.O) NT=-1 B 78
IN=INDEX B 79
T=START B 80
TNEXT-SWM. B 81
CALL DIFFUN (NS2,TNEXT,C,RT) B 82
CALL SAVLIN (T.C.NS2) B 83
NFRST=2 B 84
IF (T-START) 60,60,50 B 85
45 TNEXT=AMIN1(TPRNT,STOPP) B 86
IF (T.EQ.START) GO TO 55 B 87
50 IF (NQ.LT.O.AND.IABS(NQ).NE.O) IN=3 B 88
55 CALL DRIVES (NS2,T,H,C,TNEXT,ERR,21,IN,IA,JA) B 89
T=TNEXT B 90
IF (T.GE.STOPP) TPRNT=STOPP B 91
60 IF (T.EQ.START) TIMNW=START B 92
IF (T.NE.START) TIMNW=TPRNT B 93
IF (T.EQ.START) GO TO 65 B 94
IF (NQ.LT.O.AND.IABS(NQ).NE.O) TIMNW=T B 95
NHS=NHS+1 B 96
IF (NHS.GE.MS.OR.T.GE.TPRNT) GO TO 65 B 97
GO TO 45 B 98
65 NHSO B 99
CALL DIFFUN (NS2,TIMNW,C,RT) B 100
A-9
-------
IF (NOPT.LE.O) NTP=NTP+(N/10)+4 B 101
IF (NOPT.GT.O) lfTP*tfl?+(N/10)*3+(NR-l)/lCH-12 B 102
IF (NTP.GT.MT.OR.T.EQ.START) WRITE (IOUT,120) B 103
IF (T.EQ.START) 60 TO 70 B 104
IF (NOPT.LE.O.AND.NTP.LE.MT) 60 TO 75 B 105
70 CONTINUE B 106
IF (NTP.6T.MT) NTP=0 B 107
IF (NN.6E.11) WRITE (IOUT,H5) (SPECIS(I),I=1,NN) B 108
IF (NN.LE.10) WRITE (IOUT,125) (SPECIS(I),I=1,NN) B 109
75 CONTINUE B 110
IF (NS.6E.11) WRITE (IOUT,130) TIMNW,(C(I),I=1,10),H,(C(I),I*11,NS B 111
12)JEMP,SG B 112
IF (NS.LE.10) WRITE (IOUT,130) TIMNW,(C(I),M,NS2),TEM?,SG,H B 113
IF (NOPT.EQ.O) WRITE (IOUT.135) B 114
IF (NOPT.EQ.O) GO TO 95 B 115
RT(NS)=0. B 116
RT(N)=0. B 117
DO 80 1=1,NS2 B 118
80 RT(NS)=RT(NS)+RT(I) B 119
WRITE (IOUT.105) (RT(I),I=1,NS) B 120
DO 90 M.NR B 121
J=KR(I,1) B 122
IF (J.EQ.O) RT(I)=0. B 123
IF (J.EQ.O) GO TO 90 B 124
JT*ITYPE(I) B 125
XT»1. B 126
DO 85 L=1,JT B 127
J=KR(I,L) B 128
X>1. B 129
IF (J.6T.O.AND.J.LE.NS2) XJ=C(J)**ISC(I,L) B 130
IF (ISC(ItL).EQ.-l) XJ=C(J)**SC(I,L) B 131
IF (J.EQ.NS) X>SG B 132
XT=XT*XJ B 133
85 CONTINUE B 134
RT(I)-XT*R(I) B 135
90 CONTINUE B 136
WRITE (IOUT.110) (RT(I),I=1,NR) B 137
95 IF (TIMNW.EQ.START) 60 TO 55 B 138
IF (T.6E.STOPP) RETURN B 139
IF (IN.EQ.3.AND.T.LT.TPRNT) GO TO 45 B 140
TPRNT=TPRNT*TFACT^TSTEP B 141
GO TO 45 B 142
C B 143
C B 144
105 FORMAT (/,10H NET RATES,2X,1P10E12.3,/,(12X,1P10E12.3)) B 145
110 FORMAT (//,2X,22HTHE REACTION RATES ARE,/,(1H .1P10E13.2)) B 146
115 FORMAT (/,4X,5HTIME ,4X,10(4X,A4,4X),/,2X,8HINTERVAL,3X,10(4X,A4,4 B 147
1X),/,(13X,10(4X,A4,4X))) B 148
120 FORMAT (1H1) B 149
125 FORMAT (/,4X,5HTIME ,4X,10(4X,A4,4X),/) B 150
A-10
-------
130 FORMAT (/,lPllE12.3/tllE12.3/,(12X,10E12.3)) B 151
135 FORMAT (1H ) B 152
END B 153-
A-ll
-------
SUBROUTINE RATES (C,N)
C
C
C THIS SUBROUTINE SETS THE INITIAL CONCENTRATIONS AND CALCULATES THE
C RATE CONSTANTS
C
C
C
C
C
C
C
C
C
C
INTEGER SPECIS,REACT
COMMON /DATA/ NR,KR(200,7) ,A( 200) ,S(200),ITITLE(7), TEMP, ERR,START,
1STOPP,PC,NP,SIG(91),IP(91),ITYPE(200),R(200),BK5SG,OILUT
COMMON /NAMES/ SPECIS(91),REACT(91),NS
COMMON /HEAT/ CV,Q(200),SC(200,7),ISC(200,3),ITEMP
COMMON /ALPHA/ IGO(4),IBLANK,MBLANK>JINTER
DIMENSION C(91)
FCT=l./298.-l./TEMP
SIGG=0.
DO 5 1=1,91
5 SIG(I)=0.
DO 15 I=1,N
DO 10 J=1,NS
IF (SPECIS(J).EQ.REACT(I)) SIG(J)=C(I)
IF (SPECIS(O).EQ.REACT(D) GO TO 15
10 CONTINUE
SIGG=SIGG+C(I)
15 CONTINUE
N=NS
M=N-1
C(N)=0.
DO 20 1=1 ,M
C(N)=C(N)+SIG(I)
20 C(I)=SIG(I)
BK=0.
IF SIG(N) DOES NOT EQUAL ZERO, IT IMPLIES THAT THAT THE CONCENTRAT
SPECIES N HAS BEEN READ. OTHERWISE M IS THE SUM OF THE INITIAL
CONCENTRATIONS.
C(N)=C(N)+SIGG
IF (SIG(N).NE.O.) BK=SIG(N)-C(N)
IF (SIG(N).NE.O.) C(N)=SIG(N)
NP=0
DO 25 1=1, NR
IF (KR(I,1).EQ.O) GOTO 25
CALCULATE THE RATE CONSTANTS
IF (ABS(S(D).EQ.O.) R(I)=A(I)
IF (ABS(S(D).NE.O.) R(I)=A(I)*EXP(S(I)*FCT)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
A-12
-------
25 CONTINUE C 51
IF (NS.LE.9) SPECIS(NS+1)=JINTER C 52
RETURN C 53
END C 54-
A-13
-------
SUBROUTINE DIFFUN (N,T,X,XT) D 1
C D 2
C******AA*************************************************************** Q 3
C 04
C THIS SUBROUTINE CALCULATES THE DERIVATIVE VECTOR OF THE ODEIS D 5
C D 6
***AAAAAAA*********************************************************** D 7
D 8
COMMON /DATA/ NR,KR(200,7),A(200),S(200),ITITLE(7),TEM5,ERR,START, D 9
1STOPP,PC,NP,SI6(91),IP(91),ITYPE(200),R(200),BK,SG,DILUT 0 10
CONON /NAMES/ SPECIS(91) ,REACT(91) ,NS 0 11
COWDN /HEAT/ CV,Q(200),SC(200,7),ISC(200,3),ITEW> D 12
COWON /STORE/ AST(35),IPL(7) JEteND,NTE^TMI,NPHOT,PHI,IL,NFRST, D 13
1IPH(30),QM(100),PM(100),PSTOP D 14
COMMON /PHOTR/ IPP,RDAT(80),RTIM(80)fRl,IN10 0 15
DIMENSION XT(N), X(N) D 16
INTEGER PHI.TMI D 17
INTEGER SPECIS D 18
P=BK D 19
DO 5 I=1,N D 20
XT(I)=-OILUT*X(I) D 21
5 P=P+X(I) D 22
SG=P D 23
IF (NFRST.EQ.l) FCT=((3355.7046E-6)*TEMP-1.)/TEW> D 24
IF (NFRST.EQ.l) GO TO 10 D 25
IF (T.EQ.TOLD) GO TO 30 D 26
IF (T.LE.l.) GOTO 30 D 27
10 IF (NPHOT.LE.O) GO TO 30 D 28
IF (NFRST.NE.l) GO TO 15 D 29
PINT=FLOAT(PHI) D 30
PINV=1./PINT D 31
15 CONTINUE D 32
IF (T.GT.PSTOP) GO TO 20 D 33
IZ«IFIX(T*PINV)+1 D 34
Z=T*PINV-FLOAT(IZ-1) D 35
IF (T.LE.PINT) Rl=PM(l)+(0.5*PM(3)*(Z-l.h0.5*PM(l)*(Z-3.)-PM(2)*( D 36
1Z-2.))*Z D 37
IF (T.LE.PINT) GO TO 20 D 38
Rl=PM(IZ)-K).25*Z*(5.*PM(IZ+l)-3.0*PM(IZ)-PM(IZ-l)-PM(IZ+2)+(PM(IZ- D 39
11)-PM(IZ)-PM(IZ+1)+PM(IZ+2))*Z) D 40
20 IF (Rl.LT.O.) Rl*0. D 41
DO 25 IK=i,IL D 42
IR=IPH(IK) D 43
R(IR)=R1*A(IR) D 44
25 IF (R(IR).LT.O.) R(IR)=0. D 45
30 TNOW=TEf*> D 46
IF (NFRST.EQ.l) GO TO 35 D 47
IF (T.EQ.TOLD) GO TO 50 D 48
IF (T.LE.l.) GO TO 50 D 49
35 CONTINUE D 50
A-14
-------
IF (T.GT.TEMEND) 60 TO 50 D 51
IF (NFRST.NE.l) GO TO 40 D 52
TINT=FLOAT(TMI) D 53
TINV-l./TINT D 54
40 CONTINUE D 55
IZ=IFIX(T*TINV)+1 D 56
Z=T*TINV-FLOAT(IZ-1) D 57
IF (T.LE.TINT) TNOW=OM(l)-H(0.5*QM(3)*(Z-l.)-»-0.5*QM(l)*(2-3.)-QM(2) D 58
1*(Z-2.))*Z D 59
IF (T.LE.TINT) GO TO 45 D 60
TNOW=QM(IZ)+0.25*Z*((5.*QM(IZ+l)-3.0*QM(IZ)-QM(IZ-l)-QM(IZ+2))+(QM D 61
1(IZ-1)-QM(IZ)-QM(IZ+1)+QM(IZ+2))*Z) 0 62
45 CONTINUE 0 63
IF (TNOW.NE.TEMP.AND.TNOW.GT.O.) FCT=((3355.7046E-6)*TNOW-1.)/TNOW D 64
50 DO 70 IR=1,NR D 65
I-KR(IR,1) D 66
IF (I.EQ.O.OR.R(IR).EQ.O.) GO TO 70 D 67
C D 68
C CHECK FOR A ZEROTH ORDER REACTION D 69
C D 70
RT=1. D 71
JT=ITYPE(IR) D 72
DO 55 L»1,JT D 73
I"KR(IR,L) D 74
IF (I.EQ.99.0R.I.EQ.O) GO TO 55 D 75
IF (I.NE.NS) RT=RT*X(I) D 76
IF (I.EQ.NS) RT=RT*P D 77
55 CONTINUE D 78
IF (ABS(S(IR)).NE.O..AND.T.GT.l.) R(IR)=A(IR)*EXP(S(IR)*FCT) D 79
RT=RT*R(IR) D 80
DO 60 1=1,JT D 81
I=KR(IR,L) D 82
IF (I.EQ.O.OR.I.GT.N) GOTO 60 D 83
XT(I)=XT(I)-RT D 84
60 CONTINUE D 85
DO 65 K=4,7 D 86
I=KR(IR,K) D 87
IF (I.EQ.O) GO TO 70 D 88
C D 89
C I WILL BE NEGATIVE IF CLEAN HAS BEEN CALLED D 90
C D 91
IF (I.LT.O) GO TO 65 D 92
XT(I)=XT(I)+SC(IR,K)*RT D 93
65 CONTINUE D 94
70 CONTINUE D 95
TEMP=TNOW D 96
TOLD=T D 97
RETURN D 98
END D 99-
A-15
-------
SUBROUTINE MATRX (C,KR)
C
C
C
C
C
C
MATRX CREATES A NR X 7 MATRIX OF THE REACTION SCHEME WHERE EACH
ROW REPRESENTS A REACTION. THE FIRST THREE COLUMNS ARE REACTANTS
AND THE LAST FOUR ARE THE PRODUCTS. THE ELEMENTS CORRESPOND
TO THE INDIVIDUAL SPECIES AND WILL BE USED AS SUBSCRIPTS.
C
C
C
C
C
C
C
C
C
C
INTEGER SPECIS.REACT
COMMON /DATA/ NR,IR(200,7),A(200),S(200),ITITLE(7)JEW>fERR,START,
1STOPP,PC,NP,SI6(91) ,IP(91) ,ITYPE(200) ,R(200) ,BK,SG,DILUT
COMMON /NAMES/ SPECIS(91),REACT(91),NS
COMMON /GEAR1/ T,H,HMIN,HMAX,EPS1,UR(XJND,NC,MF1,KFLAG1,JSTART
COMMON /ALPHA/ 160(4), IBLANK.MBLANK.JINTER
COMMON /HEAT/ CV,Q(200),SC(200,7),ISC(200,3),ITEMP
DIMENSION C(91), KR(200,7)
NOLD=NS+1
DO 90 I=1,NR
SKIP REACTIONS ALREADY PROCESSED
IF (IR(I,2).EQ.NOLD.OR.IR(I,3).EQ.NOLD) IR(I,3)=99
IF (IR(I,2).EQ.NOLD) IR(I,2)*0
IF (IR(I,1).EQ.NOLD) IR(I,1)=99
IF (IR(I,1).EQ.NOLD) IR(I,3)*99
IF (IABS(IR(I,1)).NE.100) GO TO 90
IF LESS THAN THREE REACTANTS, FILL FIRST SLOTS.
IF (KR(I.l).NE.IBLANK) GO TO 5
IF (KR(I,3).NE.IBLANK) KR(I,1)=KR(I,3)
IF (KR(I,3).NE.IBLANK) SC(I,1)=SC(I,3)
SC(I,3)=1.
KR(I,3)=IBLANK
IF (KR(I.l).NE.IBLANK) GO TO 5
KR(I,1)=KR(I,2)
KR(I,2)=IBLANK
SC(I,1)=SC(I,2)
SC(I,2)=1.
5 IF (KR(I,2).NE.IBLANK) GO TO 10
IF (KR(I,3).NE.IBLANK) SC(I,2)=SC(I,3)
SC(I,3)=1.
IF (KR(I,3).NE.IBLANK) KR(I,2)=KR(I,3)
KR(I,3)=IBLANK
10 DO 15 K=4,7
GET RID OF M AS A PRODUCT
E
E
E
E
E
E
E
E
E
E
£
E
E
E
E
E
E
E
E
E
E
I
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
1
2
4
5
6
7
8
9
« f\
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
A-16
-------
C E 51
IF (KR(I,K).EQ.MBLANK) KR(I,K)=IBLANK E 52
15 CONTINUE E 53
C E 54
C IF LESS THAN FOUR PRODUCTS, FILL FIRST SLOTS. E 55
C E 56
IF (KR(I,4).NE.IBLANK) GOTO 30 E 57
DO 20 K»l,3 E 58
INDEX=8-K E 59
IF (KR(I,INDEX).NE.IBLANK) GO TO 25 E 60
20 CONTINUE E 61
INDEX=5 E 62
25 KR(I,4)=KR(I,INDEX) E 63
KR(I,INDEX)=IBLANK E 64
SC(I,4)=SC(I,INDEX) E 65
SC(I,INDEX)=1. E 66
IF (KR(I,1).NE.IBLANK.OR.KR(I,4).NE.IBLANK) 60 TO 30 E 67
IR(I,1)=0 E 68
GO TO 90 E 69
30 CONTINUE E 70
IF (KR(I,5).NE.IBLANK) GO TO 35 E 71
IF (KR(I,7).NE.IBLANK) KR(I,5)=KR(I,7) E 72
IF (KR(I,7).NE.IBLANK) SC(I,5)=SC(I,7) E 73
SC(I,7)»1. E 74
KR(I,7HBLANK E 75
IF (KR(I,5).NE.IBLANK) GO TO 35 E 76
KR(I,5)=KR(I,6) E 77
KR(I,6)=IBLANK E 78
SC(I,5)=SC(I,6) E 79
SC(I,6)=1. E 80
35 K=KR(I,6) E 81
IF (K.NE.IBLANK) GO TO 40 E 82
KR(I,6)=KR(I,7) E 83
KR(I,7)=K E 84
SC(I,6)=SC(I,7) E 85
SC(I,7)=1. E 86
40 DO 85 J=l,7 E 87
K=KR(I,J) E 88
C E 89
C PROCESS REACTANTS HERE E 90
C E 91
IF (J.GT.3) GO TO 65 E 92
C E 93
C ALL M DEPENDENT REACTIONS ARE TO HAVE A 99 IN THE THIRD SLOT E 94
C E 95
IF (K.NE.MBLANK) GO TO 60 E 96
GO TO (45,50,55),J E 97
45 KR(I,1)=KR(I,2) E 98
SC(I,1)=SC(I,2) E 99
50 KR(I,2)=KR(I,3) E 100
A-17
-------
SC(I,2)»SC(I,3) B 101
SC(I,3)-1. E 102
KR(I,3)sMBLANK E 103
55 IR(I,3)=99 E 104
60 K=KR(I,J) E 105
C E 106
C ZERO ORDER REACTIONS HAVE 99 IN FIRST SLOT E 107
C E 108
IF (KR(U).EQ.IBLANK) IR(I,1)=99 E 109
IF (KR(I.l).EQ.IBLANK) K*99 E 110
C E 111
C ALL BLANKS ARE SET EQUAL TO ZERO E 112
C E 113
IF (J.EQ.3.AND.K.EQ.MBLANK) K=99 E 114
65 IF (K.EQ.MBLANK.OR.K.EQ.IBLANK) IR(I,J)*0 E 115
IF (K.EQ.IBLANK.OR.K.EQ.99) GO TO 85 E 116
IF (NS.NE.O) GO TO 70 E 117
NS*1 E 118
GO TO 80 E 119
70 DO 75 L»1,NS E 120
IF (K.NE.SPECIS(L)) GOTO 75 E 121
C E 122
C SLOT SET TO SPECIES NUMBER E 123
C E 124
IR(I,J)=L E 125
GO TO 85 E 126
75 CONTINUE E 127
C E 128
C IF NO SPECIES ARE FOUND, ADD ONE TO THE LIST E 129
C E 130
IF (SPECIS(NS).NE.ITEM>) NS=NS+1 E 131
80 SPECIS(NS)=K E 132
C(NS+2)-C(NS+l) E 133
C(NS+1)=C(NS) E 134
C(NS)=0. E 135
IR(I,J)-NS E 136
85 CONTINUE E 137
90 CONTINUE E 138
IF (SPECIS(NS).NE.ITENP) NS«NS+1 E 139
SPECIS(NS)=ITEff E 140
SPECIS(NS+1)=MBLANK E 141
NS=NS+1 E 142
DO 95 IK=1,NR E 143
DO 95 MT-1,3 E 144
J=IFIX(SC(IK,MT)-KJROUND) E 145
ISC(IK,MT)=J E 146
95 IF (SC(IK,WT)-FLOAT(J).GT.4.*UROUND) ISC(IK,MT)=-1 E 147
DO 100 1=1,NR E 148
IF (IR(I.l).EO.O) GOTO 100 E 149
ITYPE(I)»2 E 150
A-18
-------
C E 151
C CHECK FOR 99 CODE AND SUBSTITUTE M WHICH IS THE LAST SPECIES E 152
C E 153
IF (IR(I,1).EQ.99.AND.IR(I,3).EQ.99) IR(I,1)=NS E 154
IF (IR(I,1).EQ.NS.AND.IR(I,3).EQ.99) IR(I,3)=0 E 155
IF (IR(I,2).EQ.O.AND.IR(I,3).EQ.99) IR(I,2)=NS E 156
IF (IR(I,2).EQ.NS.AND.IR(I,3).EQ.99) IR(I,3)=0 E 157
IF (IR(I,2).EQ.O) ITYPEdH E 158
IF (IR(I,3).NE.O) ITYPE(I)-3 E 159
IF (IR(I,3).EQ.99) IR(I,3)=NS E 160
100 CONTINUE E 161
RETURN E 162
END E 163-
A-19
-------
SUBROUTINE PLOT (NTIT,NPNT,NTOT,NA«,SAVTIM,SAVCON,IOT,KCP,JGRID)
SUBROUTINE ****** PLOT
C
C
C THIS SUBROUTINE READS THE PLOT CARDS AND PLOTS THE RESULTS AS PART
C OF THE PRINTED OUTPUT ~ IT DOES NOT DRIVE A PLOTTER.
C
C SYMBOL DESCRIPTIONS -
C
C CGRID
C CHIGH
C CLOW
C CSPAN
C DATA
C IDT
C J
C JBLANK
C XONC
C JGRID
C JN
C JPLUS
C JSTAR
C JSYMB
C JVERT
C JX
THE LENGTH OF THE VERTICAL AXIS, PPM
HIGHEST CONCENTRATION VALUE, PPM
LOWEST CONCENTRATION VALUE, PPM
CONCENTRATION NORMALIZATION FACTOR
CONCENTRATION DATA POINTS, PPM, UP TO 80
IF NOT ZERO THEN PRINT RAW DATA USED FOR PLOTS
DO-LOOP INDICES OR LOCAL POINTERS
A HOLLERITH WORD OF FOUR BLANK CHARACTERS
CONCENTRATION LABELS
THE PLOTTING GRID
HOLLERITH SYMBOL FOR NO DATA POINTS
THE CHARACTER >+>
THE CHARACTER >*>
SYMBOL TO BE USED FOR PLOTTING SAVED POINTS
VERTICAL LEGEND
THE CHARACTER >X>
C K DO-LOOP INDICES OR LOCAL POINTERS
C KCP IF NOT EQUAL TO ZERO SAVE DATA FOR OFFLINE PLOTTING
C KCON CONCENTRATION COORDINATE ON GRID
C KTIM TIME COORDINATE ON GRID
C L DO-LOOP INDICES OR LOCAL POINTERS
C M DO-LOOP INDICES OR LOCAL POINTERS
C MAXCON LIMIT ON NUMBER OF VERTICAL POINTS
C MAXPNT MAXIMUM NUMBER OF SAVED TIME AND CONCENTRATION POINTS
C MAXTIM LIMIT ON NUMBER OF HORIZONTAL POINTS
C N DO-LOOP INDICES OR LOCAL POINTERS
C NAME SPECIES NAMES, ONE PER SPECIES
C NDAT NUMBER OF CONCENTRATION DATA POINTS
C NIN THE FORTRAN INPUT UNIT (NORMALLY 5)
C NOUT THE FORTRAN OUTPUT UNIT NUMBER (NORMALLY 6)
C NPLT THE NUMBER OF SPECIES TO BE PLOTTED
C NPNT NUMBER OF SAVED TIMES AND CONCENTRATIONS
C NTEST SPECIES NAME FOR TESTING
C NTIT USER-INPUT TITLE FOR PRINTOUT, 3 FOUR-CHARACTER WORDS
C NTOT TOTAL NUMBER OF SPECIES
C NX FLAG FOR CORRECT SPECIES TEST
C SAVCON SPECIES CONCENTRATIONS, PPM, ONE PER SPECIES AT 80 TIMES
C SAVTIM TIMES THAT CONCENTRATIONS ARE SAVED, MIN, UP TO 80 VALUES
C TGRID THE LENGTH OF THE HORIZONTAL AXIS, MIN
C THIGH HIGHEST TIME VALUE, MIN
C TIME TIMES AT WHICH CONCENTRATIONS ARE INPUT, MIN, UP TO 80
C TLOW LOWEST TIME VALUE, MIN
F
F
F
F
F
F
F '
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
A-20
-------
C TPRINT TIMES FOR PRINTOUT ON HORIZONTAL AXIS, MIN F 51
C TSPAN TIME NORMALIZATION FACTOR F 52
C F 53
C BEGINNING OF PROGRAM. F 54
C F 55
C ENTRY POINT F 56
C F 57
C SET DIMENSIONS OF INCOMING ARRAYS F 58
C F 59
DIMENSION SAVCON(90,80), SAVTIM(SO), NTIT(3), NAME(91) F 60
C F 61
C SET DIMENSIONS OF LOCAL ARRAYS F 62
C F 63
DIMENSION XONC(5), TIME(80), DATA(80), JN(3) F 64
DIMENSION JGRID(89,40), TPRINT(9), KVERT(52) F 65
C F 66
C DEFINE THE VERTICAL LABEL AND ESTABLISH ALPHAMERIC DATA F 67
C F 68
COMMON /APLOT/ JVERT(52,2),JBLANK,JSTAR,JPLUS,JBAR F 69
COMMON /PHOTR/ IPP,RDAT(80),RTIM(80),RR1,IN10 F 70
C F 71
C DEFINE MISCELLANEOUS DATA VALUES F 72
C F 73
COMMON /INOUT/ NIN,NOUT,ITAPE F 74
DATA MAXTIM/89/,MAXCON/40/,MAXPNT/80/ F 75
DATA TGRID/88./,CGRID/40./,JX/lHX/,JNO/lH / F 76
DATA KVERT/11*4H ,4H P ,4H H ,4H 0 ,4H T ,4H 0 ,4H L ,4H F 77
1Y ,4H S ,4H I ,4H S ,4H ,4H C ,4H 0 ,4H N ,4H S ,4H T F 78
2,4H A ,4H N ,4H T ,22*4H / F 79
C F 80
C DEFINE VERTICAL AXIS VIA ASSIGNMENT STATEMENTS F 81
C F 82
IF (IDT*IDT.EQ.O) GO TO 10 F 83
WRITE (NOUT.150) F 84
WRITE (NOUT,145) (NAME(I),I=1,NTOT) F 85
DO 5 J=1,NPNT F 86
5 WRITE (NOUT.155) SAVTIM(J),(SAVCON(I,J),I=1,NTOT) F 87
10 DO 15 J=l,40 F 88
JGRID(1,J)=JBAR F 89
15 CONTINUE F 90
JGRID(1,1)=JPLUS F 91
JGRID(1,11)=JPLUS F 92
JGRID(1,21)=JPLUS F 93
JGRID(1,31)=JPLUS F 94
JN(1)=JSTAR F 95
JN(2)=OPLUS F 96
JN(3)=JX F 97
TSPAN=-10. F 98
C F 99
C READ PLOT CONTROL CARD F 100
A-21
-------
C F 101
20 READ (NIN,130) NPLT,XONC,CLOW,CHIGH,TLOW,THIGH F 102
C F 103
C TEST FOR END PLOTTING F 104
C F 105
IF (NPLT.EQ.O) GO TO 105 F 106
WRITE (NOUT.160) NTIT F 107
C F 108
C SET NORMALIZATION FACTORS AND VERTICAL LABELS F 109
C F 110
CSPAN*CGRID/(CHIGH-CLOW) F 111
TSPAN*TGRID/(THIGH-TLOW) F 112
JVERT(1,2)=JCONC(5) F 113
JVERT(11,2)=OCONC(4) F 114
JVERT(21,2)=JCONC(3) F 115
JVERT(31,2)=JCONC(2) F 116
C F 117
C SET HORIZONTAL TIME LABELS F 118
C F 119
DO 25 J=l,9 F 120
TPRINT(J)=FLOAT(J-1)/8.*(THIGH-TLOW)+TLOW . F 121
25 CONTINUE F 122
C F 123
C CLEAR GRID F 124
C F 125
DO 35 K-1.MAXCON F 126
DO 30 J*2,MAXTIM F 127
JGRID(J,K)*OBLANK F 128
30 CONTINUE F 129
35 CONTINUE F 130
NX=0 F 131
DO 95 LK'l.NPLT F 132
READ (NIN.135) NTEST.NDAT.JSYMB F 133
C F 134
C TEST NUMBER OF DATA POINTS AND READ DATA F 135
C F 136
IF (NDAT.LE.O) GO TO 45 F 137
IF (NDAT.LE.MAXPNT) GO TO 40 F 138
WRITE (NOUT,185) MAXPNT F 139
GO TO 125 F 140
C F 141
C READ DATA POINTS F 142
C F 143
40 READ (NIN.140) (TIME(J),DATA(J),J=1,NDAT) F 144
C F 145
C TEST FOR CORRECT SPECIES NAME F 146
C F 147
45 DO 50 L-l.NTOT F 148
IF (NTEST.EQ.NAME(L)) GO TO 55 F 149
50 CONTINUE F 150
A-22
-------
WRITE (NOUT,190) NTEST F 151
NX=NX+1 F 152
IF (NPLT.EQ.1.0R.NX.EQ.3) GO TO 20 F 153
IF (NPLT.EQ.2.AND.NX.EQ.2) GO TO 20 F 154
GO TO 95 F 155
C F 156
C IF THERE ARE CALCULATED POINTS, GET THEIR COORDINATES F 157
C F 158
55 IF (NPNT.LE.O) GO TO 65 F 159
DO 60 J=1,NPNT F 160
KTIM=IFIX((SAVTIM(J)-TLOW)*TSPAN+1.5) F 161
Y=(SAVCON(L,J)-CLOW)*CSPAN-0.5 F 162
IF (Y.LT.0.0) GO TO 60 F 163
KCON=IFIX(Y) F 164
KCON=MAXCON-KCON F 165
C F 166
C CHECK FOR BEING WITHIN GRID, THEN PLACE ON GRID F 167
C F 168
IF (KTIM.LT.2) GO TO 60 F 169
IF (KCON.LT.l) GO TO 60 F 170
IF (KTIM.GT.MAXTIM) GO TO 60 F 171
IF (KCON.GT.MAXCON) GO TO 60 F 172
JGRID(KTIM,KCON)=JSYMB F 173
60 CONTINUE F 174
C F 175
C IF THERE ARE DATA POINTS, GET THEIR COORDINATES F 176
C F 177
65 IF (NDAT.LE.O) GO TO 75 F 178
DO 70 J=1,NDAT F 179
KTIM=IFIX((TIME(J)-TLOW)*TSPAN+1.5) F 180
Y=(DATA(J)-CLOW)*CSPAN-0.5 F 181
IF (Y.LT.0.0) GO TO 70 F 182
KCON=IFIX(Y) F 183
KCON=MAXCON-KCON F 184
C F 185
C CHECK FOR BEING WITHIN GRID, THEN PLACE ON GRID F 186
C F 187
IF (KTIM.LT.2) GO TO 70 F 188
IF (KCON.LT.l) GO TO 70 F 189
IF (KTIM.GT.MAXTIM) GO TO 70 F 190
IF (KCON.GT.MAXCON) GO TO 70 F 191
IF (LK.EQ.l) JGRID(KTIM,KCONHSTAR F 192
IF (LK.EQ.2) JGRID(KTIM,KCON)=JPLUS F 193
IF (LK.EQ.3) JGRID(KTIM,KCON)=JX F 194
70 CONTINUE F 195
75 IF (NDAT.GT.O) 60 TO 85 F 196
WRITE (NOUT,195) NTEST,JNO.JSYMB F 197
GO TO 90 F 198
85 WRITE (NOUT.195) NTEST,JN(LK),JSYMB F 199
C F 200
A-23
-------
C SAVE DATA FOR OFFLINE PLOTS IF DESIRED F 201
C F 202
90 IF (KCP*KCP.EQ.O) GO TO 95 F 203
WRITE (ITAPE) NTIT,NAME(L),aOW,CHIGH,TLOW,THI6H,NDAT,NPNT,(TIME(J F 204
1),>1,NDAT),(DATA(J),>1,NDAT),(SAVTIM(J),>1,NPMT),(SAVCON(L,J),J F 205
2=1,NPNT) F 206
95 CONTINUE F 207
DO 100 K=1,MAXCON F 208
WRITE (NOUT.165) JVERT(K,1)$JVERT(K,2),(JGRID(J,K),J»1,MAXTIM) F 209
100 CONTINUE F 210
C F 211
C PRINT THE HORIZONTAL AXIS AND LABELS F 212
C F 213
WRITE (NOUT.170) JCONC(l) F 214
IF (THIGH.GT.80.) WRITE (NOUT.175) TPRINT F 215
IF (THIGH.LE.80.) WRITE (NOUT.180) TPRINT F 216
GO TO 20 F 217
C F 218
C END OF SUBROUTINE ~ RETURN TO CALLER F 219
C F 220
105 IF (IPP.LE.O) RETURN F 221
IF (IN10.LE.O) RETURN F 222
IF (TSPAN.LE.O.) RETURN F 223
OVERT(1,2HH0.80 F 224
JVERT(11,2)=4H0.60 F 225
JVERT(21,2)*4H0.40 F 226
JVERT(31,2)=4H0.20 F 227
JCONCUHHO.OO F 228
CSPAN=CGRID/0.8 F 229
WRITE (NOUT.200) NTIT F 230
DO 110 I=1,MAXCON F 231
DO 110 J=2,MAXTIM F 232
110 JGRID(J,I)=JBLANK F 233
DO 115 I=1,NPNT F 234
KTIM=IFIX((RTIM(I)-TLOW)*TSPAN+1.5) F 235
Y=(RDAT(I)-aOW)*CSPAN-0.5 F 236
IF (Y.LT.Q.O) GO TO 115 F 237
KCON'IFIX(Y) F 238
KCON=MAXCON-KCON F 239
IF (KCON.LT.l) GO TO 115 F 240
IF (KTIM.LT.2) GO TO 115 F 241
IF (KTIM.GT.MAXTIM) GO TO 115 F 242
IF (KCON.GT.MAXCON) GO TO 115 F 243
JGRID(KTIM,KCON)=JSTAR F 244
115 CONTINUE F 245
DO 120 K=1,MAXCON F 246
WRITE (NOUT,165) KVERT(K),JVERT(K,2),(JGRID(J,K),>1,MAXTIM) F 247
120 CONTINUE F 248
WRITE (NOUT.170) XONC(l) F 249
IF (TPRINT(9).GT.80.) WRITE (NOUT,175) TPRINT F 250
A-24
-------
IF (TPRINT(9).LE.80.) WRITE (NOUT,180) TPRINT F 251
RETURN F 252
125 CONTINUE F 253
STOP F 254
C F 255
C LIST OF FORMAT STATEMENTS F 256
C F 257
C F 258
C F 259
130 FORMAT (5X,I2,8X,5(A4,1X),4F10.0) F 260
135 FORMAT (A4,1X,I2,1X,A1) F 261
140 FORMAT (8F10.0) F 262
145 FORMAT (/,4X,5HTIME ,4X,10(4X,A4,4X),/,(13X,10(4X,A4,4X))) F 263
150 FORMAT (1H1) F 264
155 FORMAT (1P11E12.4,/,(12X,10E12.4)) F 265
160 FORMAT (1H1,//////////,62X,3A4,/,32X,7HSPECIES,2X,5HEXPT.,2X,4HSIM F 266
1.) F 267
165 FORMAT (16X,2A4,89A1) F 268
170 FORMAT (20X,A4,1H+,8(11H +)) F 269
175 FORMAT (I25,2X,8F11.0,/,62X,14HTIME (MINUTES),/) F 270
180 FORMAT (I25,2X,8F11.3,/,62X,14HTIME (MINUTES),/) F 271
185 FORMAT (33H PROGRAM CANNOT HANDLE MORE THAN ,I4,28H PLOT POINTS - F 272
1 JOB ABORTED.) F 273
190 FORMAT (33X,13HSPECIES NAME ,A4,26H NOT FOUND IN SPECIES LIST) F 274
195 FORMAT (33X,A4,6X,A1,5X,A1) F 275
200 FORMAT (1H1,//////////,62X,3M,//) F 276
END F 277-
A-25
-------
SUBROUTINE VALU (VAL.N.NF.NL) G 1
COWON /STORE/ AST(35),IPL(7),TEI€ND,NTEMP,TMI,NPHOT,PHI,IL,NFRST, G 2
1IPH(30),OM(100),PM(100),PSTOP G 3
CCWON /GEAR1/ T,H,FWIN,H1AX,EPS1,UROUND,NC,MF1,KFLAG1,JSTART G 4
DIMENSION AC(5), AD(5), BC(5) G 5
INTEGER AC,AD,BC,AST G 6
DATA IBLK/1H /,IPER/1H./,IZRO/1HO/ G 7
AVAL'ALOGIO(VAL) G 8
IF (AVAL.EQ.O.) RETURN G 9
DO 5 1=1,5 G 10
ACUHBtK G 11
AD(I)=IZRO G 12
5 BCUHZRO G 13
IAD=IFIX(AVAL+UROUND) G 14
IF (IAD) 30,10,10 G 15
10 IREM=3-IAD G 16
IAD=IAD+1 G 17
JREM=IFIX(VAL-HJROUND) . G 18
REM-VAL-FLOAT(JREM) G 19
IF (REM.GT.UROUND.AND.IREM.GT.O) GO TO 15 G 20
CALL CONVT (JREM,AC,5) G 21
GO TO 50 G 22
15 CALL CONVT (JREM,AD,IAD) G 23
JREM=IFIX(REM*(10.**IREM)+0.1) G 24
CALL CONVT (JREM,BC,IREM) G 25
DO 20 J=1,IAD G 26
20 AC(J)=AD(J) G 27
AC(IAD+1)=IPER G 28
DO 25 K=1,IREM G 29
25 AC(K+IAD+1)=8C(K) G 30
GO TO 40 G 31
30 IAD=IABS(IFIX(AVAL-0.1))-1 G 32
IVAL=IFIX(VAL*10000.+0.1) ' G 33
CALL CONVT (IVAL,AC,5) G 34
AC(1)=IPER G 35
DO 35 1=1,IAD G 36
35 ACU+1HZRO G 37
40 IF (AC(S).NE.IZRO) GO TO 50 G 38
DO 45 K-1,4 G 39
L=6-K G 40
45 AC(L)=AC(L-1) G 41
AC(1}=IBLK G 42
60 TO 40 G 43
50 DO 55 I=NF,NL G 44
K-I-NF+1 G 45
55 AST(I)=AC(K) G 46
RETURN • G 47
END G 48-
A-26
-------
SUBROUTINE CONVT (NUM.L.N) H 1
C H 2
C SUBROUTINE CONVT CONVERTS INTEGERS TO ALPHANUMERICS H 3
C FOR PRINTING H 4
C H 5
C ASSUMES VALUE OF INTEGER IS POSITIVE H 6
C H 7
DIMENSION L(5), JDIGIT(IO) H 8
C H 9
DATA JDIGIT/1HO,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ H 10
DATA JBLANK/1H / H 11
C H 12
NI=NUM H 13
DO 5 M,N - H 14
L(I)=JBLANK H 15
5 CONTINUE H 16
C H 17
DO 10 K=1,N H 18
I=N-K+1 H 19
NEXT=NI/10 H 20
NDX=(NI-NEXT*10)+1 H 21
l(I)=JDIGIT(NDX) H 22
IF (NEXT.LE.O) GO TO 15 H 23
NI=NEXT H 24
10 CONTINUE H 25
C H 26
15 RETURN H 27
END H 28-
A-27
-------
SUBROUTINE SAVLIN (T,C,N) I 1
COWON/FRPLOT/ NIT(3),SAVCON(90,80),SA\rriM(80),JGRID(89,40),NT I 2
CCW3N /PHOTR/ IPP,RDAT(80),RTIM(80),RR1,IN10 I 3
DDCNSION C(91) I 4
DATA NFRST/1/ I 5
IF (NFRST.EQ.l) 60 TO 5 16
IF (T.EQ.TOLD) RETURN I 7
5 IF (NT.LT.O) RETURN I 8
NFRST*2 I 9
TOLD=T I 10
IF (NT.6T.80) RETURN I 11
IF (T.LT.SAVTIM(NT)) RETURN I 12
SAVTIM(NT)-T - I 13
RTIM(NT)=T I 14
DO 10 M,N I 15
10 SAVCON(I,NTK(I) I 16
RDAT(NT)=RR1 .1 17
NT*NT+1 I 18
RETURN I 19
END I 20-
A-28
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
SUBROUTINE SPLNA (N,X,Y,J,D,C,W)
DIMENSION X(10), Y(10), 0(2), C(30), W(30)
OVER THE INTERVAL X(I) TO X(I+1), THE INTERPOLATING
POLYNOMIAL
Y=Y(I)+A(I)*Z+B(I)*Z**2+E(I)*Z**3
WHERE Z=(X-X(I))/(X(W)-X(I))
IS USED. THE COEFFICIENTS A(I),B(I) AND E{I) ARE CONFUTED
BY SPLNA AND STORED IN LOCATIONS C(3*I-2),C(3*I-1) AND
C(3*I) RESPECTIVELY.
WHILE WORKING IN THE ITH INTERVALJHE VARIABLE Q WILL
REPRESENT Q=X(I+1) - X(I), AND Y(I) WILL REPRESENT
Y(H-1)-Y(I)
Q=X(2)-X(1)
YI=Y(2)-Y(1)
IF (J.EQ.2) GO TO 5
IF THE FIRST DERIVATIVE AT THE END POINTS IS GIVEN,
A(l) IS KNOWN, AND THE SECOND EQUATION BECOMES
MERELY B(1)+E(1)=YI - 0*D(1).
C(1)=0*D(D
C(2)=1.0
W(2)=YI-C(1)
GO TO 10
IF THE SECOND DERIVATIVE AT THE END POINTS IS GIVEN
B(l) IS KNOWN, THE SECOND EQUATION BECOMES
A(1)+E(1)=YI-0.5*0*Q*D(1). DURING THE SOLUTION OF
THE 3N-4 EQUATIONS,A1 WILL BE KEPT IN CELL C(2)
INSTEAD OF C(l) TO RETAIN THE TRIDIAGONAL FORM OF THE
COEFFICIENT MATRIX.
5 C(2)=0.0
W(2) =0.5*Q*Q*D(1)
10 M=N-2
IF (M.LE.O) GO TO 20
UPPER TRIANGULARIZATION OF THE TRIDIAGONAL SYSTEM OF
EQUATIONS FOR THE COEFFICIENT MATRIX FOLLOWS—
DO 15 I=1,M
AI=Q
Q=X(I+2)-X(I+l)
H=AI/Q
C(3*I)=-H/(2.0-C(3*I-1))
W(3*I)=(-YI-W(3*I-1))/(2.0-C(3*I-1))
C(3*I+1)=-H*H/(H-C(3*I))
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
• J
J
J
J
J
J
J
J
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
A-29
-------
r
i«
C
C
C
r
r
C
C
r
v»
r
L.
C
C
C
C
r
i*
W(3*I+1MYI-W(3*I))/(H-C(3*I))
YI»Y(I+2)-Y(I+l)
C(3*I+2)»1.0/(1.0-C(3*I+1))
15 W(3*I+2MYI-W(3*I+1))/(1.0-C(3*H-1))
E(N-l) IS DETERMINED DIRECTLY FROM THE LAST EQUATION
OBTAINED ABOVE, AND THE FIRST OR SECOND DERIVATIVE
VALUE GIVEN AT THE END POINT.
20 IF (J.EQ.l) 60 TO 25
C(3*N-3MQ*Q*0(2)/2.0-W(3*N-4))/(3.0-C(3*N-4))
GO TO 30
25 C(3*N-3MQ*D(2)-YI-W(3*N-4))/(2.0-C(3*N-4))
30 M-3*N-6
IF (M.LE.O) GO TO 40
BACK SOLUTION FOR ALL COEFFICENTS EXCEPT
A(l) AND B(l) FOLLOWS--
DO 35 II=1,M
T=M-I 1+3
35 C(I)-M(I)-C(I)*C(W)
40 IF (J.EQ.l) GO TO 45
IF THE SECOND DERIVATIVE IS GIVEN AT THE END POINTS,
A(l) CAN NOW BE COMPUTED FROM THE KNOWN VALUES OF
8(1) AND E(l). THEN A(l) AND B(l) ARE PUT INTO THEIR
PROPER PLACES IN THE C ARRAY.
C(1H(2)-Y(1)-W(2)-C(3)
C(2)-W(2)
RETURN
45 C(2)=W(2)-C(3)
RETURN
END
J
J
J
J
.— .1
J
J
J
... J
J
J
J
J
J
J
.— .1
J
J
J
J
J
J
... 1
J
J
J
J
... J
J
J
J
J
J
J
51
52
53
54
EC
99
56
57
58
CO
•jy
60
61
62
63
64
65
fifi
oo
67
68
cq
UJ
70
71
72
73
74
/H
75
76
77
78
7Q
fy
80
81
82
83
84
85-
A-30
-------
SUBROUTINE DRIVES (N,TO,HO,YO,TOUT,EPS,MF,INDEX,IA.JA) K 1
DIMENSION IA(1), JA(1), YO(N) K 2
DIMENSION Y(91,6) K 3
COMMON /6EAR1/ T,H,HMIN,WAX,EPSC,UROUND,NC,MFC,KFLAG,JSTART K 4
COMMON /GEAR2/ YMAX(100)/GEAR3/ERROR(100)/GEAR4/W1(90,3) K 5
COMMON /GEARS/ IW1(91,9)/GEAR6/W2(2000)/GEAR7/IW2(2000) K 6
COVMON /GEARS/ EPSJ,IPTI2,IPTI3,IPTI4,IPTR2,IPTR3,NGRP K 7
COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE,NZA,NPL,NPU,NZL,NZU,NZRO K 8
COMMON /INOUT/ INP,LOUT,ITAPE K 9
DATA NMX/90/,LENW2/2000/,LENIW2/2000/ K 10
NGP=0 K 11
IF (INDEX.EQ.4) GO TO 15 K 12
IF (INDEX.EQ.O) GO TO 30 K 13
IF (INDEX.EQ.2) GO TO 35 K 14
IF (INDEX.EQ.-l) GO TO 40 K 15
IF (INDEX.EQ.3) GO TO 45 K 16
IF (INDEX.NE.l) GO TO 135 K 17
IF (EPS.LE.O.) GO TO 120 K 18
IF (N.LE.O) GO TO 125 K 19
IF ((TO-TOUT)*HO.GE.O.) GO TO 130 K 20
MITER=MF-10*(MF/10) K 21
IF ((MITER.NE.1).AND.(MITER.NE.2)) GO TO 15 K 22
NP1=N+1 K 23
NZA=IA(NP1)-1 K 24
MAX=LENIW2/2 K 25
IPTI2=MAX+1 K 26
CALL SORDER (N,IA,JA,IW1,IW1(1,5),MAX,IW2,IW2(IPTI2),IER) K 27
IPTI2=NZA+1 K 28
IF (IPTI2+NZA-1.GT.LENIW2) GOTO 145 K 29
DO 5 I=1,NP1 K 30
5 IW1(I,2)=IA(I) K 31
DO 10 1=1,NZA K 32
10 IW2(I)=JA(I) K 33
CALL NSCORD (N,IW1(1,2),IW2,IW1(1,3),IW2(IPTI2),IW1,IW1(1,5),IW1(1 K 34
1,8)) K 35
MAXPL=(LENIW2-NZA)/2 K 36
IPTI3=IPTI2+MAXPL K 37
MAXPU=LENIW2-IPTI3+1 K 38
CALL NSSFAC (N,IW1(1,2),IW2,MAXPL,IW1(1,3),IW2(IPTI2),IW1(1,4),MAX K 39
lPU,IWl(l,5),IW2(IPTI3),IWl(l,6),Y(lt6),IWl(l,9),Y,Y(l,2),Y(l,3),IW K 40
21(l,7),IWl(l,8)fY(l,4),Y(l,5),IER) K 41
NPL«M(N,4) K 42
NPU=IW1(N,6) K 43
NZL=M(N+1,3) K 44
NZU=IW1(N+1,5) K 45
IPTR2=NZA+1 K 46
IPTR3=IPTR2+MAXO(NZA,NZL) K 47
IF (IPTR3+MAXO(NZA,NZU)-1.6T.LENW2) GO TO 145 K 48
15 DO 20 1=1,N K 49
YMAX(I)=ABS(YO(I)) K 50
A-31
-------
IF (YMAX(I).EQ.O.) YMAX(I)=1.E-10 K 51
20 Y(I,1)=YO(I) K 52
NC»N K 53
T=TO K 54
H-HO K 55
NZROO K 56
TST=EPS*1.E-10 K 57
00 25 1=1,N K 58
25 IF (Y(I.l).GT.TST) NZRO-NZROH K 59
NZRO=MAXO(NZRO,1) K 60
NOLD=NZRO K 61
HMIN=ABS(HO) K 62
HMAX=ABS(TO-TOUT)*10. K 63
HWX=AMIN1(WAX,20.) K 64
EPSC=£PS K 65
J^OMF •• K 66
JSTARTX) K 67
NO=N K 68
NMX1=NMX+1 K 69
EPSJ=SOJ^T(UROUND) K 70
NHCUT=0 K 71
GO TO 50 K 72
30 H1AX=ABS(TOUT-TOUTP)*10. K 73
FWAX*AMIN1(WAX,20.) K 74
GO TO 80 K 75
35 HMAX=ABS(TOUT-TOUTP)*10. K 76
HMAX=AMIN1(HMAX,20.) K 77
IF ((T-TOUT)*H.GE.O.) GO TO 150 K 78
GO TO 85 K 79
40 IF ((T-TOUT)*H.GE.O.) GO TO 140 K 80
JSTART*-! K 81
NC*N K 82
EPSC=EPS K 83
^FC=MF K 84
45 CONTINUE K 85
50 CALL STIFFS (Y,NO,IA,JA,W1,NMX,IW1,NMX1) K 86
KGO=1-KFLAG K 87
GO TO (55,95,110,100),KGO K 88
55 CONTINUE K 89
D=0. K 90
NZRO=0 K 91
DO 70 1*1,NC K 92
IF (Y(I.l).GE.O.) GOTO 65 K 93
NGP=NGP+1 K 94
DO 60 J=l,6 K 95
K=(J-1)*N+I K 96
60 Y(K,1)=0. K 97
65 CONTINUE K 98
IF (Y(I,1).GT.TST) NZRO=NZRO+1 K 99
AYI=ABS(Y(I,1)) K 100
A-32
-------
YMAX(I)=AMAX1(1.E-10,AYI) K 101
70 D=D+(AYI/YMAX(I))**2 K 102
NZRO=MAXO(NZRO,1) K 103
IF (NZRO.NE.NOLD) JSTART=-1 K 104
D=D*(UROUND/EPS)**2 K 105
DO 75 II-l.NC K 106
75 YO(II)=Y(II,1) K 107
CALL SAVLIN (T.YO.N) K 108
IF (D.GT.FLOAT(N)) GO TO 115 K 109
IF (INDEX.EQ.3) GO TO 150 K 110
IF (INDEX.EQ.2) GO TO 85 K 111
80 IF ((T-TOUT)*H.LT.O.) GO TO 45 K 112
CALL INTERP (TOUT,Y,NO,YO) K 113
GO TO 160 K 114
85 IF (T.GE.TOUT) GO TO 90 K 115
IF (((T+H)-TOUT).LE.O.) GO TO 45 K 116
H=(TOUT-T)*(1.+4.*UROUND) K 117
JSTART=-1 K 118
GO TO 45 K 119
90 JSTART=-1 K 120
H=AMIN1(H,1.) K 121
GO TO 150 K 122
95 CONTINUE K 123
100 IF (NHCUT.EQ.10) GO TO 105 K 124
NHCUT=NHCUT+1 K 125
HMIN=.1*HMIN K 126
H=.1*H K 127
JSTART=-1 K 128
GO TO 45 K 129
105 WRITE (LOUT.165) K 130
IF (KGO.EQ.4) WRITE (LOUT,180) T K 131
STOP K 132
110 WRITE (LOUT.170) T,H K 133
STOP K 134
115 WRITE (LOUT.175) T K 135
KFLAG=-2 K 136
STOP K 137
120 WRITE (LOUT,185) K 138
STOP K 139
125 WRITE (LOUT,190) K 140
STOP K 141
130 WRITE (LOUT,195) K 142
STOP K 143
135 WRITE (LOUT,200) INDEX K 144
STOP K 145
140 WRITE (LOUT.205) T,TOUT,H K 146
STOP K 147
145 WRITE (LOUT.210) K 148
STOP K 149
150 TOUT=T K 150
A-33
-------
DO 155 M,N K 151
155 YO(I)=Y(I,1) K 152
CALL SAVLIN (TOUT,YO,N) K 153
160 INDEX=KFLAG K 154
TOUTP=TOUT K 155
HO=HUSED K 156
IF (KFLAG.NE.O) HO=H . K 157
RETURN K 158
C K 159
C K 160
165 FORMAT (//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//) K 161
170 FORMAT (//35H KFLAG « -2 FROM INTEGRATOR AT T » .E16.8.5H H =,E16 K 162
1.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) K 163
175 FORMAT (//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/56H EPS K 164
1TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/) K 165
180 FORMAT (//35H KFLAG = -3 FROM INTEGRATOR AT T » .E16.8/45H CORREC K 166
1TOR CONVERGENCE COULD NOT BE ACHIEVED/) K 167
185 FORMAT (//28H ILLEGAL INPUT.. EPS .LE. O.//) K 168
190 FORMAT (//25H ILLEGAL INPUT.. N .LE. O//) K 169
195 FORMAT (//36H ILLEGAL INPUT.. (TO-TOUT)*H .GE. O.//) K 170
200 FORMAT (//24H ILLEGAL INPUT.. INDEX =,I5//) K 171
205 FORMAT (//44H INDEX = -1 ON INPUT WITH (T-TOUT)*H .GE. 0./4H T =,E K 172
116.8,9H TOUT =,E16.8,6H H *,E16.8) K 173
210 FORMAT (//42H INSUFFICIENT WORKING STORAGE IN IW2 OR W2//) K 174
END K 175-
A-34
-------
SUBROUTINE STIFFS (Y,NO,IA,JA,W1,NMX,IW1,NMX1) L 1
COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART L 2
COMMON /GEAR2/ YMAX(100)/GEAR3/ERROR(100) L 3
COMMON /GEARS/ W2(2000)/GEAR7/IW2(2000) L 4
COMMON /GEARS/ EPSO,IPTI2,IPTI3,IPTI4,IPTR2,IPTR3,NGRP L 5
COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE,IDUMMY(5),NZRO L 6
COMMON /DATA/ NR,KR(200,7),A(200),S(200),ITITLE(7),TEM),ERR,START, L 7
1STOPP,PC,NP,SIG(91),IP(91),ITYPE(200),R(200),BK,SG,DILUT L 8
COMMON /HEAT/ CV,Q(200),SC(200,7),ISC(200,3),ITEMP L 9
DIMENSION Y(NO,1), IA(1), JA(1), W1(NMX,1), IWl(NMXl.l) L 10
DIMENSION EL(13), TQ(4), RT(3) L 11
DATA EL(2)/1./,OLDLO/1./ L 12
KFLAG=0 L 13
TOLD=T L 14
IF (JSTART.GT.O) GO TO 50 L 15
IF (JSTART.NE.O) GO TO 10 L 16
CALL DIFFUN (N,T,Y,W1) L 17
DO 5 I-l.N L 18
5 Y(I,2)=H*W1(I,1) L 19
METH=MF/10 L 20
MITER=MF-10*METH L 21
NQ=1 L 22
L=2 L 23
IDOUB=3 . L 24
RMAX=1.E4 L 25
RC=0. L 26
CRATE=1. L 27
HOLD=H L 28
MFOLD=MF L 29
NSTEP=0 L 30
NSTEPJ=0 L 31
NFE=1 L 32
NJE=0 L 33
IRET=3 L 34
GO TO 15 L 35
10 IF (MF.EQ.MFOLD) GO TO 25 L 36
MEO=METH L 37
MIO=MITER L 38
METH=MF/10 L 39
MITER=MF-10*METH L 40
MFOLD=MF L 41
IF (MITER.NE.MIO) IWEVAL=MITER L 42
IF (METH.EQ.MEO) GO TO 25 L 43
IDOUB=L+1 L 44
IRET-1 L 45
15 CALL COSET (METH,NQ,EL,TQ,MAXDER) L 46
LMAX=MAXDER+1 L 47
RC=RC*EL(1)/OLDLO L 48
OLDLO=EL(1) L 49
20 FN=FLOAT(NZRO) L 50
A-35
-------
EDN=FN*(TQ(1)*EPS)**2 L 51
E=FN*(TQ(2)*EPS)**2 L 52
EUP=FN*(TQ(3)*EPS)**2 L 53
BND=FN*(TQ(4)*EPS)**2 L 54
EPSOLD-EPS L 55
NOLD=NZRO L 56
GOTO (30,35,50),IRET L 57
25 IF (EPS.EQ.EPSOLD.AND.NZRO.EQ.NOLD) GO TO 30 L 58
IRET=1 L 59
GO TO 20 L 60
30 IF (H.EQ.HOLD) GO TO 50 L 61
RH=H/HOLD L 62
H=HOLD L 63
IREDO=3 L 64
GO TO 40 L 65
35 RH=AMAX1(RH,HMIN/ABS(H)) L 66
40 RH=AMIN1(RH,HMAX/ABS(H),RMAX) L 67
Rl-1. L 68
DO 45 J=2,L L 69
R1=R1*RH L 70
DO 45 I-l.N L 71
45 Y(I,J)=Y(I,J)*R1 L 72
H=H*RH L 73
RC=RC*RH . L 74
IDOUB-L+1 L 75
IF (IREDO.EQ.O) GO TO 290 L 76
50 IF (ABS(RC-1.).GT.0.3) IWEVAL41ITER L 77
IF (NSTEP.GE.NSTEPJ*20) IWEVAL*MITER . L 78
T=T+H L 79
DO 55 J1-1.NQ L 80
DO 55 J2=J1,NQ L 81
J»(NO*J1)-J2 L 82
DO 55 I-l.N L 83
55 Y(1,J)-Y(I,J)+Y(I,M) L 84
60 DO 65 1=1,N L 85
65 ERROR(I)=0. L 86
M=0 L 87
CALL DIFFUN (N,T,Y,W1(1,2)) L 88
NFE=NFE+1 L 89
IF (IWEVAL.LE.O) GO TO 140 L 90
IWEVAW) L 91
RC=1. L 92
NJE-NJE+1 L 93
NSTEPJ=NSTEP L 94
CON=-H*EL(1) L 95
ISV=M L 96
LSV=L L 97
NZ=IA(Nfl)-l L 98
DO 70 I=1,NZ L 99
70 W2(I)=0. L 100
A-36
-------
DO 125 IR=1,NR L 101
IF (KR(IR,1).EQ.O.OR.KR(IR,1).EQ.99) GO TO 125 L 102
MT=ITYPE(IR) L 103
DO 85 1=1,W L 104
IF (KR(IR,I).EQ.N+2) 60 TO 85 L 105
JX=I+l-I/3*3 L 106
KX=I+2-I/2*3 L 107
J=KR(IR,JX) L 108
K=KR(IR,KX) L 109
XJ=1. L 110
IF (J.EQ.O) GO TO 75 L 111
XJ=Y(J,1) L 112
IF (J.EQ.N+2) XJ=SG L 113
75 XK=1. L 114
IF (K.EQ.O) GO TO 80 L 115
XK=Y(K,1) L 116
IF (K.EQ.N+2) XK=SG L 117
80 RT(I)=R(IR)*XJ*XK L 118
85 CONTINUE L 119
DO 120 K=1,MT L 120
I=KR(IR,K) L 121
IF (I.EQ.N+2) GO TO 120 L 122
DO 100 L«1,MT L 123
J=KR(IR,L) L 124
IF (J.EQ.N+2) GO TO 100 L 125
M=IA(J)-1 L 126
90 M=M+1 L 127
IF (I-JA(M)) 90,95,90 L 128
95 W2(M)=W2(M)-RT(L) " L 129
100 CONTINUE L 130
DO 115 L=4,7 L 131
J=KR(IR,L) L 132
M=IA(I)-1 L 133
IF (J) 115,120,105 L 134
105 M=M+1 L 135
IF (J-JA(M)) 105,110,105 L 136
110 W2(M)=W2(M)+RT(K)*SC(IR,L) L 137
115 CONTINUE L 138
120 CONTINUE L 139
125 CONTINUE L 140
DO 135 J-l.N L 141
KMIN=IA(J) L 142
KMAX=IA(J+1)-1 L 143
DO 130 K=KMIN,KMAX L 144
W2(K)=W2(K)*CON L 145
IF (JA(K).EQ.J) W2(K)=W2(K)+1.-CON*DIIUT L 146
130 CONTINUE L 147
135 CONTINUE L 148
CALL NSCORA (N,IA,JA,W2,IW1(1,2),W2(IPTR3),W2(IPTR2),IW1,IW1(1,7), L 149
1IW1(1,8)) L 150
A-37
-------
CALL NSNFAC (N,IW1(1,2),IW2,W2,IW1(1,3),IW2(IPTI2),IW1(1,4),W2(IPT L 151
lR2),Wl(l,3),IWl(lt5),IW2(IPTI3),IWl(l,6),W2(IPTR3),Wl,IWl(l,7),M L 152
2(1,8),IER) L 153
M*ISV L 154
L=LSV L 155
IF (IER.NE.O) 60 TO 160 L 156
140 DO 145 I-l.N L 157
IF (M.LE.O) 60 TO 145 L 158
IF (-H*W1(I,2)*10..6T.Y(I,1)) 60 TO 155 L 159
145 W1(I,1)»H*W1(I,2)-(Y(I,2)+ERROR(I)) L 160
CALL NSBSLV (N,IWl,IWl,M(l,3),IW2(IPTI2),IWl(li4),W2(IPTR2),Wl(l L 161
1,3),IW1(1,5),IW2(IPTI3),IW1(1,6),W2(IPTR3),W1(1,2),W1,W2) L 162
D*0. L 163
DO 150 1=1,N L 164
ERROR(I)=ERROR(I)+W1(I,2) L 165
D=D+(W1(I,2)/YMAX(I))**2 L 166
150 W1(I,1H(I,1)+EL(1)*ERROR(I) L 167
IF (M.NE.O) CRATE=AMAX1(.9*CRATE,D/D1) L 168
IF ((D*AMIN1(1.,2.*CRATE)).LE.BND) 60 TO 175 L 169
Dl=0 L 170
M=MH L 171
IF (M.EQ.3) 60 TO 155 L 172
CALL DIFRJN (N,T,W1,W1(1,2)) L 173
GO TO 140 L 174
155 NFE=NFE+2 L 175
IF (IWEVAL.EQ.-l) GO TO 170 L 176
160T*TOLD L 177
RMAX=2. L 178
DO 165 J1=1,NQ L 179
DO 165 J2=J1,NQ . L 180
J=(NQfJl)-J2 L 181
DO 165 I-1.N L 182
165 Y(I,J)«Y(I,J)-Y(I,J+1) L 183
IF (ABS(H).LE.HMIN*1.00001) GO TO 285 L 184
RH».25 ' L 185
IREDO=1 L 186
GO TO 35 L 187
170 IWEVAL=MITER L 188
GO TO 60 L 189
175 IF (MITER.NE.O) IWEVAL=-1 L 190
NFE=NFE4« L 191
D«0. L 192
DO 180 I=1,N L 193
180 WHERROR(I)/YMAX(I))**2 L 194
IF (D.GT.E) GO TO 195 L 195
KFLAG=0 L 196
IREDOO L 197
NSTEP=NSTEP+1 L 198
HUSED=H L 199
NQUSED=NQ L 200
A-38
-------
DO 185 J=1,L L 201
DO 185 1=1,N L 202
185 Y(I,J)=Y(I,J)+EL(J)*ERROR(I) L 203
IF (IDOUB.EQ.l) 60 TO 205 L 204
IDOUB=IDOUB-1 L 205
IF (IDOUB.GT.l) GO TO 295 L 206
IF (L.EQ.LMAX) GO TO 295 L 207
DO 190 1=1,N L 208
190 Y(I,LMAX)=ERROR(I) L 209
GO TO 295 L 210
195 KFLAG=KFLAG-1 L 211
T=TOLD L 212
DO 200 J1=1,NQ L 213
DO 200 J2=J1,NQ L 214
J=(NCH-J1)-J2 L 215
DO 200 I=1,N L 216
200 Y(I,J)=Y(I,J)-Y(I,J+1) L 217
RMAX=2. L 218
IF (ABS(H).LE.HMIN*!.00001) GO TO 275 L 219
IF (KFLAG.LE.-3) GO TO 265 L 220
IREDO=2 L 221
PR3=1.E+20 L 222
GO TO 215 L 223
205 PR3=1.E+20 L 224
IF (L.EQ.LMAX) GO TO 215 L 225
Dl=0. L 226
DO 210 1=1,N L 227
210 D1=D1+((ERROR(I)-Y(I,LMAX))/YMAX(I))**2 L 228
ENQ3=.5/FLOAT(L+1) L 229
PR3=((Dl/EUP)**ENQ3)*1.4+1.4E-6 L 230
215 ENQ2=.5/FLOAT(L) L 231
PR2=((D/E)**ENQ2)*1.2+1.2E-6 L 232
PR1=1.E+20 L 233
IF (NQ.EQ.l) GO TO 225 L 234
D=0. L 235
DO 220 1=1,N L 236
220 D=D+(Y(I,L)/YMAX(I))**2 L 237
ENQ1=.5/FLOAT(NQ) L 238
PRH(D/EDN)**ENQl)*1.3+1.3E-6 L 239
225 IF (PR2.LE.PR3) GO TO 230 L 240
IF (PR3.LT.PR1) GO TO 240 L 241
GO TO 235 L 242
230 IF (PR2.GT.PR1) GO TO 235 L 243
NEWQ=NQ L 244
RH=1./PR2 L 245
GO TO 255 L 246
235 NEWQ=NQ-1 L 247
RH=1./PR1 L 248
GO TO 255 L 249
240 NEWQ=L L 250
A-39
-------
RH=1./PR3 L 251
IF (RH.LT.1.1) GO TO 250 L 252
DO 245 I-l.N L 253
245 Y(I,NEWO+1KRROR(I)*EL(L)/FLOAT(L) L 254
GO TO 260 L 255
250 IDOUB=10 L 256
GO TO 295 L 257
255 IF ((KFLAG.EQ.O).AND.(RH.LT.1.1)) GO TO 250 L 258
IF (NEWQ.EQ.NQ) GO TO 35 L 259
260 NQ=NEWQ L 260
L-NQfl L 261
IRET=2 L 262
GO TO 15 L 263
265 IF (KFLA6.EQ.-9) GO TO 280 L 264
RH=10.**KFLAG L 265
RH=AMAX1(HMIN/ABS(H),RH) L 266
H=H*RH L 267
CALL DIFFUN (N,T,Y,W1) L 268
NFE=NFE+1 L 269
00 270 1=1,N L 270
270 Y(I,2)=H*W1(I,1) L 271
IWEVAL*MITER L 272
IDOUB=10 L 273
IF (NQ.EQ.l) GO TO 50 L 274
NQ=1 L 275
L*2 L 276
IRET=3 L 277
GO TO 15 L 278
275 KFLAG=-1 L 279
GO TO 295 L 280
280 KFLAG=-2 L 281
GO TO 295 L 282
285 KFLAG=-3 L 283
GO TO 295 L 284
290 RMAX=100. L 285
295 HOLD=H L 286
JSTART=NQ L 287
RETURN L 288
END L 289-
A-40
-------
SUBROUTINE COSET (METH,NQ,EL,TQ,MAXDER) M 1
DIMENSION PERTST(12,2,3), EL(13), TQ(4) M 2
DATA PERTST/1.,1.,2.,1...3158,.07407,.01391,.002182,.0002945,.0000 M 3
13492,.000003692,.0000003524,1.,1.,.5,.1667,.04167,1.,1.,1.,1.,1.,1 M 4
2.,1.,2.,12.,24.,37.89,53.33,70.08,87.97,106.9,126.7,147.4,168.8,19 M 5
31.0,2.0,4.5,7.333,10.42,13.7,1.,!.,!.,1.,1.,1.,1.,12.0,24.0,37.89, M 6
453.33,70.08,87.97,106.9,126.7,147.4,168.8,191.0,1.,3.0,6.0,9.167,1 M 7
52.5,1.,!.,!.,!.,!.,!.,!.,I./ M 8
5 MAXDER=5 M 9
GOTO (10,15,20,25,30),NQ M 10
10 ELUH.O M 11
GO TO 35 M 12
15 EL(1)=6.6666666666667E-01 M 13
EL(3)=3.3333333333333E-01 M 14
GO TO 35 M 15
20 EL(1)=5.4545454545455E-01 M 16
EL(3HL(1) M 17
EL(4)=9.0909090909091E-02 M 18
GO TO 35 M 19
25 EL(1)=0.48 M 20
EL(3)=0.7 M 21
EL(4)=0.2 M 22
EL(5)=0.02 M 23
GO TO 35 M 24
30 EL(1)=4.3795620437956E-01 M 25
EL(3)=8.2116788321168E-01 M 26
EL(4)=3.1021897810219E-01 M 27
EL(5)=5.4744525547445E-02 M 28
EL(6)=3.6496350364964E-03 M 29
35 DO 40 K*l,3 M 30
40 TQ(K)=PERTST(NQ,METH,K) M 31
TQ(4)=.5*TQ(2)/FLOAT(NO+2) M 32
RETURN M 33
END M 34-
A-41
-------
SUBROUTINE NSCORA (N,IA,JA,A,IAP,JAWORK,AWORK,C,IR,ICT) N 1
INTEGER IA(1),JA(1),IAP(1),JAWORK(1),C(1),IR(1),ICT(1) N 2
REAL A(1),AWORK(1) N 3
00 5 K»1,N N 4
I(X<(K) N 5
5 IR(ICK)=K N 6
JMIN-1 N 7
DO 15 K=1,N N 8
ICK-C(K) N 9
JW\X=JMIN+IA(ICK+1)-IA(ICK)-1 N 10
IF (JMIN.GT.JMAX) SO TO 15 N 11
IAINK=IA(ICK)-1 N 12
DO 10 J=JMIN,JMAX N 13
IAINK»IAINK+1 N 14
JAOUTJ=JA(IAINK) N 15
JAOUTJ=IR(JAOUTJ) N 16
JAWORK(J)=JAOUTJ N 17
10 AWORK(J)=A(IAINK) N 18
15 JMIKWMAX-H N 19
DO 20 I-l.N N 20
20 ICT(I)«IAP(I) N 21
JMIN-1 N 22
DO 30 I-l.N N 23
ICK-C(I) N 24
JMflX«JMIM*IA(IOC+l)-IA(ICK)-l N 25
IF (JMIN.GT.JMAX) GO TO 30 N 26
DO 25 J*JMIN,JMAX N 27
JAOUTJ=JAWORK(J) N 28
ICTJ=ICT(JAOUTJ) N 29
A(ICTJ)=AWORK(J) N 30
ICT(JAOUTJ)=ICTJ+1 N 31
25 CONTINUE N 32
30 JMIN=JMAX+1 N 33
RETURN N 34
END N 35-
A-42
-------
SUBROUTINE NSNFAC (N,IA,JA,A,IL,JL,ISL,L,D,IU,JU,ISU,U,X,IRL,JRL,I 0 1
1ER) 0 2
INTEGER IA(1),JA(1),IL(1),JL(1),ISL(1) 0 3
INTEGER IU(1),JU(1),ISU(1),IRL(1),JRL(1) 0 4
REAL A(1),L(1),D(1),U(1),X(1) 0 5
REAL LKI 06
IER=0 0 7
00 5 K=1,N 0 8
IRL(K)=IL(K) 0 9
5 JRL(K)=0 0 10
DO 90 K=1,N 0 11
X(K)=0. 0 12
11=0 0 13
IF (JRL(K).EQ.O) GO TO 15 0 14
I=JRL(K) 0 15
10 I2=JRL(I) 0 16
JRL(I)=I1 0 17
11=1 0 18
X(I)=0. 0 19
1=12 0 20
IF (I.NE.O) GO TO 10 0 21
15 JMIN=ISU(K) 0 22
JMAX=JMIN+IU(K+1)-IU(K)-1 0 23
IF (JMIN.GT.JMAX) GO TO 25 0 24
DO 20 J=JMIN,JMAX 0 25
JUJ=JU(J) 0 26
20 X(JUJ)=0. 0 27
25 JMIN=IA(K) 0 28
JMAX=IA(K+1)-1 0 29
DO 30 J=JMIN,JMAX 0 30
JAJ=JA(J) 0 31
30 X(JAJ)=A(J) 0 32
1=11 0 33
IF (I.EQ.O) GO TO 50 0 34
35 IRLI=IRL(I) 0 35
LKI=-X(I) 0 36
L(IRLI)=-LKI 0 37
JMIN=IU(I) 0 38
JMAX=IU(I+1)-1 0 39
IF (JMIN.GT.JMAX) GO TO 45 0 40
ISUB=ISU(I)-1 0 41
DO 40 J=JMIN,JMAX 0 42
ISUB=ISUB+1 0 43
JUJ=JU(ISUB) 0 44
40 X(JUJ)=X(JUJ)+LKI*U(J) 0 45
45 I=JRL(I) 0 46
IF (I.NE.O) GO TO 35 0 47
50 IF (X(K).EQ.O.) GO TO 95 0 48
DK=1./X(K) 0 49
D(K)=DK • 0 50
A-43
-------
IF (K.EQ.N) GO TO 90 0 51
JMIN-IU(K) 0 52
JMAX=IU(K+1)-1 0 53
IF (JMIN.GT.JMAX) GO TO 60 0 54
ISUB=ISU(K)-1 0 55
DO 55 J=JMIN,JMAX 0 56
ISUB=ISUB+1 0 57
JUJ-JU(ISUB) 0 58
55 U(JH(JUO)*DK 0 59
60 CONTINUE 0 60
1*11 0 61
IF (I.EQ.O) GO TO 85 0 62
65 IRL(I)-IRL(IH 0 63
I1=JRL(I) 0 64
IF (IRL(I).GE.IL(I+1)) GOTO 80 0 65
ISLB=IRL(I)-IL(I)+ISL(I) 0 66
J-JL(ISLB) 0 67
70 IF (I.GT.JRL(J)) GO TO 75 0 68
J=JRL(J) 0 69
GO TO 70 0 70
75 JRL(I)=JRL(J) 0 71
JRL(J)=I . 0 72
80 1=11 0 73
IF (I.NE.O) GO TO 65 0 74
85 ISLK=ISL(K) 0 75
IF (IRL(K).GE.IL(K-H)) GO TO 90 0 76
>JL(ISLK) 0 77
JRL(K)=JRL(J) 0 78
JRL(J)=K 0 79
90 CONTINUE 0 80
RETURN 0 81
95 IER=K 0 82
RETURN 0 83
END 0 84-
A-44
-------
SUBROUTINE NSBSLV (N,R,C,IL,JL,ISL,L,D,IU,JU,ISU,U,X,B,Y) P 1
DIMENSION R(l), IL(1), JL(1), IU(1), JU(1), C(l), ISL(l), ISU(l) P 2
DIMENSION L(l), X(l), B(l), U(l), Y(l), D(l) P 3
INTEGER R,RK,C,CK P 4
REAL L P 5
DO 5 K=1,N P 6
RK=R(K) P 7
5 Y(K)=B(RK) P 8
DO 15 K=1,N P 9
JMIN=IL(K) P 10
JMAX=IL(K+1)-1 P 11
YK=-D(K)*Y(K) P 12
Y(K)=-YK P 13
IF (JMIN.GT.J1AX) GO TO 15 P 14
ISLB=ISL(K)-1 P 15
DO 10 J=JMIN,JMAX P 16
ISLB=ISLB+1 P 17
JLJ=JL(ISLB) P 18
10 Y(JLJ)=Y(JLJ)+YK*L(J) P 19
15 CONTINUE P 20
K=N P 21
DO 30 1=1,N P 22
SUM=-Y(K) P 23
JMIN=IU(K) P 24
JMAX=IU(K+1)-1 P 25
IF (JMIN.GT.JMAX) GO TO 25 P 26
ISUB=ISU(K)-1 P 27
DO 20 J=JMIN,JMAX P 28
ISUB=ISUB+1 P 29
JUJ=JU(ISUB) P 30
20 SUM=SUM+U(J)*Y(JUJ) P 31
25 Y(K)=-SUM P 32
CK=C(K) P 33
X(CK)=-SUM P 34
30 K=K-1 P 35
RETURN P 36
END P 37-
A-45
-------
SUBROUTINE YSMER (A,K,A1) Q 1
INTEGER A,A1(2) Q 2
COWON /INOUT/ INP.LOUT.ITAPE Q 3
WRITE (LOUT.5) A,K,A1(1),A1(2) Q 4
RETURN Q 5
C Q 6
C Q 7
C 08
5 FORMAT (1X,A10,I6,2A10) Q 9
END Q 10-
A-46
-------
SUBROUTINE INTERP (TOUT,Y,NO,YO) R 1
COVWON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART R 2
DIMENSION YO(NO), Y(NO,1) R 3
DO 5 I=1,N R 4
5 YO(I)=Y(I,1) R 5
L=JSTART+1 R 6
S=(TOUT-T)/H R 7
Sl=l. R 8
DO 15 J=2,L R 9
S1=S1*S R 10
DO 10 1=1,N R 11
10 YO(I)=YO(I)+S1*Y(I,J) R 12
15 CONTINUE R 13
RETURN R 14
END R 15-
A-47
-------
SUBROUTINE SPARS (IA,JA,N) S 1
COMMON /DATA/ NRtKR(200,7),A(200),S(200),ITITLE(7),TEMP,ERR,START, S 2
1STOPP,PC,NP,SI6(91),IP(91),ITYPE(200),R(200),BK,S6,DILUT S 3
DIMENSION IA(N), JA(N) S 4
DO 5 I«1,M S 5
5 IA(IH S 6
KT=0 S 7
IA(N+1)=1 S 8
JA(1)=0 S 9
DO 70 IR»1,NR S 10
IF (KR(IR,l).EQ.N+2) GO TO 70 S 11
IF (KR(IR,1).EQ.O.OR.KR(IR,1).EQ.99) GO TO 70 S 12
MT=ITYPE(IR) S 13
DO 65 K=1,MT S 14
I*KR(IR,K) S 15
IF (KR(IR,K).EQ.N+2) GO TO 65 S 16
DO 30 L=1,MT S 17
J=KR(IR,L) S 18
IF (KR(IR,L).EQ.N+2) GO TO 30 S 19
K1=IA(J) S 20
K2-IA(J*1)-1 S 21
IF (K1.GT.K2) GO TO 15 S 22
DO 10 MHC1.K2 S 23
10 IF (I.EQ.JA(M)) GO TO 30 S 24
15 DO 20 M=J,N S 25
20 IA(M*1)«IA(MH)+1 S 26
KT=KT+1 S 27
KD=KT-K2 S 28
K2*K2+1 S 29
DO 25 M«1,KD S 30
25 JA(KT+2-M)=JA(KT+l-M) S 31
JA(K2)=I S 32
30 CONTINUE S 33
Kl-IA(I) S 34
DO 60 1=4,7 S 35
K2-IA(W)-1 S 36
J=KR(IR,L) S 37
IF (KR(IR,L).EQ.N+2) GO TO 60 S 38
IF (J) 60,65,35 S 39
35 IF (K1.GT.K2) GO TO 45 S 40
DO 40 M=K1,K2 S 41
IF (J.EQ.JA(M)) GO TO 60 S 42
40 CONTINUE S 43
45 DO 50 M=I,N S 44
50 IA(M*1)-IA(M*1)+1 S 45
KT=KT+1 S 46
KD=KT-K2 S 47
K2=K2+1 S 48
DO 55 M-l.KD S 49
55 JA(KT+2-M)=JA(KT+l-M) S 50
A-48
-------
JA(K2)=J S 51
60 CONTINUE S 52
65 CONTINUE S 53
70 CONTINUE S 54
DO 80 I-l.N S 55
K1=IA(I)+1 S 56
K2=IA(I+1)-1 S 57
IF (K1.GT.K2) GO TO 80 S 58
MT=K2-K1+1 S 59
DO 75 K=1,MT S 60
DO 75 M=K1,K2 S 61
IF (JA(M).GT.JA(M-l)) GO TO 75 S 62
J=JA(M-1) S 63
JA(M-1)=JA(M) S 64
JA(M)=J S 65
75 CONTINUE S 66
80 CONTINUE S 67
M=N S 68
DO 95 I=1,M S 69
IF (lA(I-H).GT.IA(I)) GO TO 95 S 70
NM=I+1 S 71
NN=N+1 S 72
KMIN=IA(NM) S 73
KMAX=IA(NN) S 74
DO 85 J=KMIN,KMAX S 75
KM=KMAX+KMIN-J S 76
85 JA(KM)=JA(KM-1) S 77
KNOW=IA(I) S 78
JA(KNOW)=I S 79
DO 90 LL=NM,NN S 80
90 IA(LL)=IA(LL)+1 S 81
95 CONTINUE S 82
RETURN S 83
END S 84-
A-49
-------
SUBROUTINE SORDER (N,IA,JA,P,Q,MAX,V,L,IER) T 1
INTEGER IA(1),JA(1),P(1),Q(1),V(1),L(1) T 2
INTEGER S,SFS,PI,PJ,PK,VI,VJ,VK,QVK,DTHR,DMIN T 3
IER*0 T 4
DO 5 S=1,MAX T 5
5 L(S)-S*1 T 6
SFS=1 T 7
L(MAX)0 T 8
DO 10 K-l.N T 9
P(K)=K T 10
Q(K}* T 11
V(K)=1 T 12
10 L(K)=0 T 13
SFS=SFS+N T 14
DO 50 K=1,N T 15
JMIN-IA(K) T 16
JMAX=IA(K+1)-1 T 17
IF (JMIN.GT.JMAX+1) GO TO 145 T 18
KDIAG=0 T 19
DO 45 J=JMIN,JMAX T 20
VJ=JA(J) T 21
IF (VJ.NE.K) GO TO 15 - T 22
KDIAG=1 T 23
GO TO 45 T 24
15 LLK=K T 25
20 LK=LLK T 26
LLX=L(LK) T 27
IF (LLK.EQ.O) GO TO 25 T 28
IF (V(LLK)-VJ) 20,30,25 T 29
25 LLK-SFS T 30
IF (LLK.EQ.O) GO TO 150 T 31
SFS<(SFS) T 32
V(K)=V(KH T 33
V(LLK)=VJ T 34
L(LLKK(LK) T 35
L(LK)=LLK T 36
30 UK=VJ T 37
35 LK=LLK T 38
LLK=L(LK) T 39
IF (LLK.EQ.O) GO TO 40 T 40
IF (V(LLK)-K) 35,45,40 T 41
40 LLK=SFS T 42
IF (LLK.EQ.O) GO TO 150 T 43
SFS=L(SFS) T 44
V(VJH(VJ)+1 T 45
V(LLK)=K T 46
L(LLK)=L(LK) T 47
L(LK)=LLK T 48
45 CONTINUE T 49
IF (KDIAG.EQ.O) GO TO 160 T 50
A-50
-------
50 CONTINUE T 51
J=0 T 52
DTHR=0 T 53
DMIN=N T 54
1=0 T 55
55 1=1+1 T 56
IF (I.GT.N) GO TO 140 T 57
JMIN=MAXO(J+1,I) T 58
IF (JMIN.GT.N) GO TO 70 T 59
60 DO 65 J=JMIN,N T 60
VI=P(J) T 61
IF (V(VI).LE.DTHR) GO TO 75 T 62
IF (V(VI).LT.DMIN) DMIN=V(VI) T 63
65 CONTINUE T 64
70 DTHR=DMIN T 65
DMIN=N T 66
JMIN=I T 67
GO TO 60 T 68
75 PJ=P(I) T 69
P(J)=PJ T 70
Q(PJ)=J T 71
PI=VI T 72
P(I)=PI T 73
Q(PI)=I T 74
LI=VI T 75
80 LI=L(LI) T 76
IF (LI.EQ.O) GO TO 105 T 77
VK=V(LI) T 78
LLK=VK T 79
LJ=VI T 80
85 LJ=L(LJ) T 81
IF (LJ.EQ.O) GO TO 100 T 82
VJ=V(LJ) T 83
IF (VJ.EQ.VK) GO TO 85 T 84
90 LK=LLK T 85
LLK=L(LK) T 86
IF (LLK.EQ.O) GO TO 95 T 87
IF (V(LLK)-VJ) 90,85,95 T 88
95 LLK=SFS T 89
IF (LLK.EQ.O) GO TO 155 T 90
SFS=L(SFS) T 91
V(VK)=V(VK)+1 T 92
V(LLK)=VJ T 93
L(LLK)=L(LK) T 94
L(LK)=LLK T 95
GO TO 85 T 96
100 IF (V(VK).GT.V(VI)) GO TO 80 T 97
1=1+1 T 98
QVK=Q(VK) T 99
PI=P(I) T 100
A-51
-------
P(QVK)*PI T 101
Q(PI)*QVK T 102
P(I)=VK T 103
Q(VK)=I T 104
60 TO 80 T 105
105 LI-VI T 106
110 IF (L(LI).EQ.O) GO TO 135 T 107
LI-L(LI) T 108
VK=V(LI) T 109
LLK=VK T 110
QVK=MINO(Q(VK),I) T 111
115 LK=LLK T 112
LLK=L(LK) T 113
IF (LLK.EQ.O) GO TO 120 T 114
VJ=V(LLK) • T 115
IF (Q(VO).GT.QVK) GO TO 115 T 116
V(VK)=V(VK)-1 T 117
L(LKK(LLK) T 118
L(LLK)=SFS T 119
SFS=LLK T 120
LLK=LK T 121
GO TO 115 T 122
120 IF (Q(VK).LE.I) GO TO 130 T 123
IF (V(VK).LE.DTHR) GO TO 125 T 124
IF ((DTHR.LT.V(VK)).AND.(V(VK).LT.DMIN)) DMIN=V(VK) T 125
GO TO 110 T 126
125 J=MINO(Q(VK)-1,J) T 127
DTHR=V(VK) T 128
GO TO 110 T 129
130 L(LK)=SFS T 130
SFS=L(VK) ' T 131
GO TO 110 T 132
135 L(LI)=SFS T 133
SFS=L(VI) T 134
GO TO 55 T 135
140 RETURN T 136
145 CALL YSMER (3HROW,K,13H OF A IS NULL) T 137
GO TO 165 T 138
150 CALL YSMER (3HROW,K,16H EXCEEDS STORAGE) T 139
GO TO 165 T 140
155 CALL YSMER (6HVERTEX.VI.16H EXCEEDS STORAGE) T 141
GO TO 165 T 142
160 CALL YSMER (6HCOLUMN,K,19H.. DIAGONAL MISSING) T 143
165 IER=1 T 144
RETURN T 145
END T 146-
A-52
-------
SUBROUTINE NSSFAC (N,IA,JA,MAXPL,IL,JL,ISL,MAXPU,IU,JU,ISU,P,V,IRA U 1
1,JRA,IRAC,IRL,JRL,IRU,JRU,IER) U 2
INTEGER IA(1),JA(1),IRA(1),JRA(1),IL(1),JL(1),ISL(1) U 3
INTEGER IU(1),JU(1),ISU(1),IRL(1),JRL(1),IRU(1),JRU(1) U 4
INTEGER P(1),V(1),IRAC(1) U 5
INTEGER VI,VJ,VK,PK,PPK,PI,CEND,REND U 6
IER=0 U 7
DO 5 K=1,N U 8
IRAC(K)=0 U 9
JRA(K)=0 U 10
JRL(K)=0 U 11
5 JRU(K)=0 U 12
DO 10 K=1,N U 13
IAK=IA(K) U 14
IF (IAK.GT.IA(K+1)) GO TO 200 U 15
IF (JA(IAK).GT.K) GO TO 10 U 16
JAIAK=JA(IAK) • U 17
JRA(K)=IRAC(JAIAK) U 18
IRAC(JAIAK)=K U 19
10 IRA(K)=IAK U 20
JLPTR=0 U 21
IL(1)=1 U 22
JUPTR=0 U 23
IU(1)=1 U 24
DO 195 K=1,N U 25
P(l)=l U 26
V(1)=N+1 U 27
LSFS=2 U 28
VJ=IRAC(K) U 29
IF (VJ.EQ.O) GO TO 30 U 30
15 PPK=1 U 31
20 PK=PPK U 32
PPK=P(PK) U 33
IF (V(PPK)-VJ) 20,220,25 U 34
25 P(PK)=LSFS U 35
V(LSFS)=VJ U 36
P(LSFS)=PPK U 37
LSFS=LSFS+1 U 38
VJ=JRA(VJ) U 39
IF (VJ.NE.O) GO TO 15 U 40
30 LASTI=0 U 41
I=K U 42
35 I=JRU(I) U 43
IF (I.EQ.O) GO TO 60 U 44
PPK=1 U 45
JMIN=IRL(I) U 46
JMAX=ISL(I)+IL(I+1)-IL(I)-1 U 47
IF (LASTI.GT.I) GO TO 40 U 48
LASTI=I U 49
LASTID=JMAX-JMIN U 50
A-53
-------
IF (JL(JMIN).NE.K) LASTID=tASTICH-l U 51
40 IF (JMIN.6T.JW) SO TO 35 U 52
DO 55 J=JMIN, J1AX U 53
VJ=OL(J) U 54
45 PK=PPK U 55
PPK=P(PK) U 56
IF (V(PPK)-VJ) 45,55,50 U 57
50 P(PK)=LSFS U 58
V(LSFSHJ U 59
P(LSFS)=PPK U 60
PPK=LSFS U 61
LSFS-LSFS+1 U 62
55 CONTINUE U 63
60 TO 35 U 64
60 PI-P(l) U 65
IF (V(PI).NE.K) 60 TO 225 U 66
IF (LASTI.EQ.O) 60 TO 65 U 67
IF (LASTID.NE.LSFS-3) 60 TO 65 U 68
IRLL=IRL(LASTI) U 69
ISL(K)»IRLL+1 U 70
IF (JL(IRLL).NE.K) ISL(K)-ISL(K)-l U 71
IL(K-H)«IL(K)-HASnD U 72
IRL(KHSL(K) U 73
60 TO 80 U 74
65 ISL(K)=JLPTR+1 . U 75
PI-P(l) U 76
PI-P(PI) U 77
VI-V(PI) U 78
70 IF (VI.GT.N) 60 TO 75 U 79
JLPTR=JLPTR+1 U 80
IF (JLPTR.6T.MAXPL) 60 TO 230 U 81
JL(JLPTR)*VI U 82
PI-P(PI) U 83
VI-V(PI) U 84
60 TO 70 U 85
75 IRL(K)«ISL(K) U 86
IL(K+1)-IL(K)+OLPTR-ISL(1C)+1 U 87
80 P(l)=l U 88
V(1)=N+1 U 89
LSFS=2 U 90
JMIN*IRA(K) U 91
JMAX»IA(K+1)-1 U 92
IF (JMIN.6T.JMAX) 60 TO 100 U 93
DO 95 J=OMIN,JMAX U 94
VJ=JA(J) U 95
PPK=1 U 96
85 PK=PPK U 97
PPK=P(PK) U 98
IF (V(PPK)-VJ) 85,205,90 U 99
90 P(PK)=LSFS U 100
A-54
-------
V(LSFS)=VJ U 101
P(LSFS)=PPK U 102
95 LSFS=LSFS+1 U 103
100 LASTI=0 U 104
I-K U 105
105 I»JRL(I) U 106
IF (I.EQ.O) 60 TO 130 U 107
PPK=1 U 108
JMIN=IRU(I) U 109
JMAX-ISU(I)+IU(I+1)-IU(I)-1 U 110
IF (LASTI.GT.I) GO TO 110 U 111
LASTI=I U 112
LASTID=JMAX-JMIN U 113
IF (JU(JMIN).NE.K) LASTID=LASTI[>1 U 114
110 IF (JMIN.GT.JMAX) GO TO 105 U 115
DO 125 J=JMIN,JMAX U 116
VJ=JU(J) U 117
115 PK=PPK U 118
PPK=P(PK) U 119
IF (V(PPK)-VJ) 115,125,120 U 120
120 P(PK)=LSFS U 121
V(LSFS)=VJ U 122
P(LSFS)=PPK U 123
PPK=LSFS U 124
125 LSFS=LSFS+1 U 125
GO TO 105 U 126
130 PI=P(1) U 127
IF (V(PI).NE.K) GO TO 210 U 128
IF (LASTI.EQ.O) GO TO 135 U 129
IF (LASTID.NE.LSFS-3) GO TO 135 U 130
IRUL=IRU(LASTI) U 131
ISU(K)=IRUL+1 U 132
IF (JU(IRUL).NE.K) ISU(K)=ISU(K)-1 U 133
IU(K+1)=IU(K)+LASTID U 134
IRU(K)=ISU(K) U 135
GO TO 150 U 136
135 ISU(K)=JUPTR+1 U 137
PI=P(1) U 138
PI=P(PI) U 139
VI=V(PI) U 140
140 IF (VI.GT.N) GO TO 145 U 141
JUPTR=JUPTR+1 U 142
IF (JUPTR.GT.MAXPU) GO TO 215 U 143
JU(JUPTR)=VI U 144
PI=P(PI) U 145
VI=V(PI) U 146
GO TO 140 U 147
145 IRU(K)=ISU(K) U 148
IU(K+1)»IU(K)+JUPTR-ISU(KH U 149
150 I=K U 150
A-55
-------
155 Il-JRL(I) U 151
CEND-ISL(I)+IL(I+1)-IL(I) U 152
IF (IRL(I).GE.CEND) 60 TO 160 U 153
IRLI«IRL(I) U 154
J-JL(IRLI) U 155
JRL(I)=JRL(J) U 156
ORL(J)-I U 157
160 1=11 U 158
IF (I.EQ.O) GOTO 165 U 159
IRL(I)-IRL(IH U 160
GO TO 155 U 161
165 I * U 162
170 Il-JRU(I) U 163
REND*ISU(I)+IU(I+1)-IU(I) U 164
IF (IRU(I).GE.REND) GO TO 175 U 165
IRUI-IRU(I) U 166
.WUflRUI) U 167
JRU(I)=JRU(J) U 168
JRU(J)=I U 169
175 1=11 U 170
IF (I.EQ.O) GO TO 180 U 171
IRU(I)«IRU(IH U 172
GO TO 170 U 173
180 I=IRAC(K) U 174
IF (I.EQ.O) GO TO 195 U 175
185 I1=JRA(I) U 176
IRA(I)=IRA(IH U 177
IF (IRA(I).GE.IA(I+1)) GO TO 190 U 178
IRAI«IRA(I) U 179
IF (JA(IRAI).GT.I) GO TO 190 U 180
JAIRAI=OA(IRAI) U 181
JRA(I)=IRAC(JAIRAI) U 182
IRAC(JAIRAI)=I U 183
190 1=11 U 184
IF (I.NE.O) GO TO 185 U 185
195 CONTINUE U 186
ISL(N)=OLPTR U 187
ISU(N)=JUPTR U 188
RETURN U 189
200 CALL YSMER (3HROW,K,13H OF A IS NULL) U 190
GO TO 235 U 191
205 CALL YSMER (3HROW,K,20H HAS DUPLICATE ENTRY) U 192
GO TO 235 U 193
210 CALL YSMER (3HROW,K,17H HAS A NULL PIVOT) U 194
GO TO 235 U 195
215 CALL YSMER (3HROW,K,19H EXCEEDS JU STORAGE) U 196
GO TO 235 U 197
220 CALL YSMER (3HCOL,K,20H HAS DUPLICATE ENTRY) U 198
GO TO 235 U 199
225 CALL YSMER (3HCOL,K,17H HAS A NULL PIVOT) U 200
A-56
-------
GO TO 235 U 201
230 CALL YSMER (3HCOL,K,19H EXCEEDS JL STORAGE) U 202
235 IER=1 U 203
RETURN U 204
END - U 205-
A-57
-------
SUBROUTINE NSCORD (N,IA,JA,IAWORK,JAWORK,C,IR,ICT) V 1
INTEGER IA(1),JA(1),IAWORK(1),JAWDRK(1),C(1),IR(1),ICT(1) V 2
DO 5 I«1,N V 3
5 ICT(I)=0 V 4
IAWDRK(1)=1 V 5
DO 15 K-l.N V 6
ICK-C(K) V 7
JMIN=IAWORK(K) V 8
JMAX=JMIN+IA(ICK+1)-IA(ICK)-1 V 9
IAWORK(K+1)=JMAX+1 V 10
IF (JMIN.GT.JMAX) GO TO 15 V 11
IAINK=IA(ICK)-1 V 12
DO 10 J»JMIN,JMAX V 13
IAINK-IAIMM V 14
JAOUTJ=OA(IAINK) V 15
JAOUTJ=IR(JAOUTJ) V 16
JAWDRK(J)=JAOUTJ V 17
10 ICT(JAOUTJHCT(JAOUTJH V 18
15 CONTINUE ' V 19
IA(1H V 20
DO 20 I=1,N V 21
IA(M)=IA(IHCT(I) V 22
20 ICT(I)=IA(I) V 23
DO 30 I=1,N V 24
JMIN=IAWORK(I) V 25
JMAX=IAWORK(I+1)-1 V 26
!F (JMIN.GT.JMAX) GO TO 30 V 27
DO 25 J=JMIN,JMAX V 28
JAOUTJ=JAWORK(J) V 29
ICTJ*ICT(JAOUTJ) V 30
JA(ICTJ)«I V 31
25 ICT(JAOUTJ)=ICTJ+1 V 32
30 CONTINUE V 33
RETURN V 34
END V 35-
A-58
-------
BLOCK DATA ALPHA! W 1
COMMON /ALPHA/ IGO(4),IBLANK,MBLANK,JINTER W 2
COMMON /APLOT/ JVERT(52,2),JBLANK,JSTAR,JPLUS,JBAR W 3
COMMON /GEAR1/ T.GUESS,HMIN,HMAX,EPS1,UROUND,NC,MF1,KFLAG1,JSTART W 4
COMMON /HEAT/ CV,Q(200)>SC(200,7),ISC(200,3),ITEMP W 5
COMMON /STORE/ AST(35),IPL(7),TEMEND,NTEMP,TMI,NPHOT,PHI,IL,NFRST, W 6
1IPH(30),QM(100),PM(100),PSTOP W 7
COMMON/PHOTR/ IPP,RDAT(80),RTIM(80),RR1,MO W 8
DATA IGO(1)/4HMORE/,IGO(2)/4HCONT/,IGO(3)/4HPLOT/,IGO(4)/4H /,! W 9
1BLANK/4H /,MBLANK/4HM /.JINTER/4HINTV/ W 10
DATA JVERT/12*4H ,4H C ,4H 0 ,4H N ,4H C ,4H E ,4H N ,4H W 11
IT ,4H R ,4H A ,4H T ,4H I ,4H 0 ,4H N ,4H ,4H P ,4H P W 12
2,4H M ,75*4H / W 13
DATA JBLANK/4H /,JSTAR/1H*/,JPLUS/1H+/,JBAR/1HI/ W 14
DATA ITEMP/4HTEMP/,TEMEND/0.0/,PSTOP/0.0/,IN10/1/ W 15
DATA UROUND/7.5E-15/ W 16
END W 17-
A-59
-------
APPENDIX B
INPUT DATA--CHEMK FORMAT
B-l
-------
APPENDIX B
INPUT DATA--CHEMK FORMAT
69 SAVE
1 R03
2 0
3 RO 03
4 R02 03
J R02 0
6 OH 03
7 1102 03
o on no2
9 Oil CO
10 no HO
1 1 no H03
12 H02 (103 ICO
is ro 1102
14 n02 1102
IS 0 PAR
16 Oil PAH
17 0 OLE
'fl 0 OLE
. ? Oil OLE
20 03 OLE
S 1 03 OLE
S2 0 ETII
23 0 nil
24 Oil ETII
23 03 ETII
26 RO AC03
27 NO nnn2
2(1 no RA02
24 RO HE02
TO RO ME02
31 HO ME02
32 03 HD02
S3 03 RA02
34 03 HE02
33 ou cum
36 OU CARD
37 CARD
30 aa
39 aa
-0 CARD
* PAH x
42 1(02 AC03
•S3 PAfl
•;•> 1102 AC03
->3 1102 ME02
46 MO CRIC
47 .102 CHIC
4.1 CARO CRIC
•»? no nciic
a> .102 ncnc
a i CAIUI note
13 >:iiu:
nil I:MIU
H4 Clllli
15 item;
oo IICIMI
97 III:IK;
ia tn:uc
vt uu ARO
60 UU AIIO
6 1 OU A1\O
. 000324
RO 0
03
H02
R03
RO
1102
on
I1N03
H02 CO2
N02 H02
R02 H02
R02 OU
KE02 Oil
HE02
HEU2 AC03 X
CAIUJ PAR
RA02
CARO CRIC
CARD ncnc x
HC02 1102 CO
CAIIO I'AR
RD02
CARO CRIC
R02 HE02 C02
R02 CARD CARO 1102
N03 CARD CARO 1103
H02 CARO neoa x
R02 CARD 1102
RRAT
CARD CARB H02
CARB CARB H02
CARB H02
11UU CO
AC03 X
00
1102 1102 CO
NE02 X U02 CO
CO
PAN
AC03 H02
R02 CARD
.103 CARD
!»02 CAtUJ P,\I\
.•O3
2. 700 E* 03
2.?oor»on
4.200L+04
O.OOOE-03
8.000E-03
6.000002
6.0IIOO02
1.200I1+04
2.400E-03
3.UOOE+03
1.200E+04
1 . 200F>04
4.000E+03
8 ••••£* 93
S.000E+02
B.«OOE+««
2.000E+92
S.OOOE^OO
9.SOOOU3
9.300E+03
I.OOOE*00
I.OOOH+02
I.OOOE»02
3.SOOE-01
I.OOOE*03
2.000C+03
2.QOOE-02
4.000E+03
4.000C+03
l.200t+04
O. 000l>03
2.000C*03
1 . 20Or.« 04
n.oo(in»oo
2.ouoi>o:i
ft .Tnoi'to'j
U. 4IIIII lli'l
V. (MIDI. Hl|
1 i &OOI1*U«
II. linn uu
4. U&1IO03
a.»o«K»oi
«.000t>03
I.60OE«03
1,SOOE»04
1400.*
2430.0
1000.0
1S23.0
-10600.0
12300.0
62 CARD W
B-2
-------
*• HO Y RO CAM PAH 3.«0«E+01
04 HO Y 103 AERO 1.S00E+OI
•• HO3 Y CAM CAM 3.S00E+04
M 03 Y AERO 6.009E-01
47 OH GLY B02 OC CO 1.000E+04
•a cc Y Y Y 1.099E+92
M CLY KE02 B02 CC 2.928E-01
•DM
TO SAVE .••6324
I 1102 HO 0 .3
•7 CAM OA .00064
40 CAHB CO .90932
M HO HE02 R02 CAHB HE02 X 2779.
30 HO JJE02 H02 CAHB 802 923ft.
«« CLY KE02 B02 CC .»72
Tl RX OB .»3
72 co .003
73 03 .001
TO H20 -.••6324
BC-237
10 69. 60.
•O H02 OLE PAR AHO CARS B20
RX CO ETH .. .
.377 .1*6 .09 7.38 .177 .1*12 24000.
.994 .9 .87S
36«. 3*3.2
PLOT
EC-237
3 •••• «.2» «.49 «.6« •.86 9. .8 •. 4«9.
03 23 0
•. .992 IB. .992 39. .91 48. .»42
M. .11 78. .IBS 90. .264 198. .338
12». .391 138. .482 180. .811 168. .887
1M. .894 198. .62 2i«. .64 228. .68
24*. .683 288. .63 270. .648 288. .638
3M. .628 318. .613 330. .693 348. .894
36*. .884
HO 23 »
• . .377 IS. .272 30. .169 48. .*83
M. . »4 78. .024 90. .«1B 108. .913
12*. .911 138. .01 ISO. .097 168. .007
IB*. .006 198. .006 210. .006 228. .006
24«. .••« 288. .006 270. .008 288. .008
30*. .993 318. .008 330. .098 348. .90S
360. .008
•02 28 2
0. .106 18. .198 39. .29 48. .383
M. .368 78. .38 99. .324 108. .297
120. .268 138. .234 180. .208 168. .171
180. .142 198. .117 210. .096 22S. .07
240. .074 288. .068 270. .063 285. .08
300. .039 318. .038 330. .037 348. .037
369. .036
BLANK CARD
BLAHK CARD
B-3
-------
APPENDIX C
INPUT DATA-ERA FORMAT
c-i
-------
APPENDIX C
INPUT DATA-ERA FORMAT
n
KB
o
HO
no3
«os
01
003
01
08
NO
NO
NO2
NO
B03
0
OH
O
O
oa
03
03
0
0
OH
03
NO
NO
NO
NO
rt
NO
NO
03
00
03
OU
OH
CAM
CAM
FUR
.102
PAN
BOS
U02
NO
N02
CAM
no
NO2
CAM
CHIC
CRIC
CHIC
ncnc
ncnc
HUM:
note
OH
OH
on
CAM
no
NO
Noa
03
OH
CC
CLY
SAVE
O3 3
03 4
O 3
03 4
03 7
N02 8
CO *
NO 10
Noa ii
N03 020 12
D02 13
U02 14
P.XR IS
PAR 16
OLE 17
OLE IB
OLE l«
OLE 20
OLE 21
ETO 22
ETtl 23
ETO 24
CHI 23
ACO3 26
RI102 27
RA02 28
HE02 3»
30
HC02 31
HE02 33
ROO2 33
RAOS 34
NEO2 33
CAM 36
CAM 37
38
34
40
41
X 43
AC03 43
.000334
HO
03
1(03
H03
HO
BO2
OB
BH03
003
3. N02
3. H02
R02
HE02
HE03
HE03
CAM
RA02
CAM
CAM
HE03
CAM
MOS
CAM
H03
KM
N02
R02
HEOS
N02
NRAT
2. CAM
3. CAM
CAM
003
AC03
OA
3. 803
AC03
HE02
CRIC
CRIC
CRIC
ncnc
ncnc
new
CO
PAH
AC03
ARO
ARO
ARO
V .
T
T
T
T
CLT
43
46
47
48
49
SB
31
32
S3
34
as
36
S?
SB
S«
60
41
62
63
44
49
66
47
68
6« S. T
rt
NO2
HOS
N02
N03
CO
3. 003
HC02
HC03
CAM
•02
BOS
OH
M
M3
3. CAM
AERO
COS
OB
00
AC03
PAR
CRIC
ncnc
002
PAR
CRIC
HE03
3. CAM
3. CAM
CAM
X
CAM
B02
802
HOS
CO
X
CO
H02
CAM
CARD
CARD
CAM
C02
00
BOS
2. 002
a. v
CLY
CLY
CAM
AIM
CC
CO
COS
CO
X
v
PAH
00
06
1.000E*00
4.400E*06
S.300OOI
4.M6C-62
1.340e«04
r.r0oe*0i
1480.
«4«0.
isaa.
X
X
CO
COS
003
DOS
pp
DOS
B03
PAR
PAR
1.S60E-03 -10*00.
1.200E+04
l.SOOE+04
2.000E+41
2I700E+03
2.700E»03
4.2OOE+04
B.OOOE-03
B.OOOE-03
6.000E*02
6.000E*02
1.200E»04
2.400E-03
3.800E«03
1.200E+O4
1.200E+04
4.000E+03
1 . OOOE+03
a.000£*93
S.OOOE+02
3.000E«>«0
3. OOOE+02
S.OOOE«00
9.300E+03
l!oO0E»00
1. OOOE+02
l.OOOE+62
3.300E-01
2l OOOE+03
a.aooE-oa tssoe,
4.000E+03
4.000E+03
1.200E+04
a.OOOE+«3
2.000E+03
1.200E+O4
a.OOOE+03
2.00UE+03
2.400E»02
3.400E»OS
4.230E»0S
a.SOOE+6!
6.AOOE+03
1.600E+63
I.S00E*«4
l.000E»0a
3.M6E+41
1.000E«04
3!023C-«t
C-2
-------
•DM
74 t SAVE .88*324
Mtt t HO O .8
CARB 38 OB CO ,MM4
CARB 41 CO .66632
HO HEM 29 1102 CARB PP 2778.
RO KE03 31 B02 CARB H02 9230.
CLY 76 HE02 102 CC .972
RX 71 OH .63
72 CO .663
03 73 .Ml
B20 74 -.•••324
EC-237
!• 60. 66.
HO HOa OLE PAR ARO GARB H30
RX CO tTB
.377 .1*6 .99 7.38 .177 .1*13 24666.
.••4 .9 .878
36*. 363.2
run
CC-337
3 •••• «.2« 6.4* •.«• «.M ». .S t. 4M.
O3 29 O
• . .M2 18. .t«2 3«. .•! 48. .»42
«0. .11 78. .108 90. ,264 1*8. .938
126. .391 138. .463 IB*. .811 168. .887
IBO. .894 198. .62 210. .64 228. .68
24*. .689 238. .68 37». .648 263. .638
3««. .620 318. .413 939. .6t3 348. .894
3*6. .884
HO 23 H
• . .37? 18. .272 36. .169 48. .683
6«. .64 78. .424 90. .fl!8 105. •••13
126. .611 138. .81 180. .697 168. .8*7
180. .0*6 198. .666 210. .006 229. .8*6
246. .866 388. .806 276. .60S 283. .8*8
3*6. .668 318. .888 336. .668 348. .888
366. .888
H02 29 2
6. .166 18. .198 36. .29 48. .383
66. .368 78. .98 98. .384 108. .297
136. .268 138. .234 186. .268 168. .171
186. .142 198. .117 21*. .696 228. .67
246. .674 288. .668 270. .663 288. .68
3*6. .639 318. .638 336. .637 348. .637
368. .836
BLAR CAM)
CARD
C-3
-------
APPENDIX D
OUTPUT DATA
(using EPA input format for the chemical reactions)
D-l
-------
APPENDIX D
OUTPUT DATA
33 S S
5 u u u u u
S O 9 9 a 9
U a is 9 tt o
» •» 9 a 9
§
t»
U
2
• • — M«->-999O-9OO9O9999999999O999
II I I I I I I I I I I I I I I I I I I I I I I I I I
U U U U
x u
8 S ? =
o u S o < £ u
n tt el n a n a
rinfiontfifi ft ooococ.
o
I]
S « s
s § !
i i i f
- 8
* 3
y r
Hi
.§
2
n
n
e
9B
W U X U
n
« Cl S 0
2 1 S I
n ct n
fi n ti
2 = 2 2 * 2 « t 2 * ?, « S 8 S S S S S S S
0-2
-------
e
i
o
o o
i i
o e
i i
o o> e o
i i i i
e o
i i
• tt e n n e
e e e e e e
* -* + * * *
u u u u u
e o e o
e e e o
- e e «a
n
s ;
n n
e e
o n »
e e o
n n ri
e e e
ct -
e e
N n
e o
o •
+ *
u u
n » B
22!
sis
e a ci
u u u
e e o
o e e
o e o
P. — — ~
i ?
a. a.
ri
o o
u u
ri n ci
222°
x S
§ g
.
u o
CI
o
u
fl CI >- ».
= 00 J J
o = e >• o 3
H g q
< S 5
u o u
N N
n ri
o ri o
U0OU
-------
?'••••••
I I I I I I
4 2 a 2 2 2 2
*
3 3
3 I
siiii.
>• »i
• * » »
2 S 8 S 8
: s s
D-4
-------
2
>•
• «SSSe«SSS
u'uuuuuuuwu
^
a
i
S §
o < <
e u u u
fl b J
Sa o e o
o u E B
a o
o u
N
§
* s
B E W U
O
2 g d S S3
— B — <* ~ » ~ tt m
-------
; 2888 Xv 88
ES i lit |H
I. ••••«!••
a .; a - - ~ « -
«• c« « a ct a ~ M
* »*?*???
M U Hi 11
8 22^:22:
9! !!!•!!!!
I 22:2:2;
M Is IIHI
;M : 2 2 2 2 2 :
is! iiiiiii
• « -< a a o • o •
J * ei -•<>•«•- a
— " « _ n • n •> -•
aki IB w u ill hi ia M u
S I ? S S 2 S S S
• £ & ••«•••»
* 5 i » » ? i * * *
^MU feiUMUliiblMlil
K8:898S9
*Btl»H«5N
2 - ? ? • S ? ? ! •
Mi: 1 I I I i I! i
•>•» ei — *•«!• — —
• «l
-;Sa»*---*-"
s .EgS,^
S4 &J «4 7 *P^*Jfll^^JO
!*• » Z!r C ^
-------
52,
fj O O O O O O
O
I
u III
-> O
C— -J O
-
rj
3
O
O
a
OOCO3OO
oo
ii
-f n
so
ii
§§
1 1 N
9 - «>
I I
U
I
*
O
nci
00
UU
NO
oor o
-------
tiff
ess.
* — c -
» C r-
c to
' '.
. .
~
»««»
— «ri
e-c:--z
W M H* iL w. —
»r. :is tut ;
erac- — e -r ti
c r: .-3 ci c j
c r. s c c f
-»*> I. e>»
e c £ i. - i c
ci :i I- £ II i
fan B a o
ee eee
« i i iii
UUid VklU
SiJ §«s
c * - ti e n
I I I I I I I
uuuuuuu
« e> eo
*—»
eee
ii
-
» c =-
ens
» v » e B e
te-o-»
at.cn
eeee
r»ene
ecs-r?
eeee eee
o —«M
oecn
Oefr>
eeee eee
i i i * t t i
Hg. uuuu uuu
u < £ «N;iri i»no
— BU — e«s>B aee
O«(- o«-r.o -v^
uuuuuuuu
aaca i»-r
eeee
iiit
uuuu
B —«t>
oeee
tiii
uuuu
esoo
n —«•— — o — <
slis
eooe
i i t i
uuuu
FiSHC'
nee«»^
< — aesn — o 11 o-
SB
ne
ee
4- +
uu
ec>
en
tin
»«»*»««»
eeeeeee
i i i i t i i
uuuuuuu
« ti f n « ss »
t— «» — —
oeeeeee
i i i i i i i
M U U U M U (io>«— o —
i i i i i i i
uu uuu L;
n n » — n c
— o—
e» — « N r « n u
\»
ecoeeeeo
i i i i i i i i
uuuuuuuu
•ariGtic-etti-
o —r»oe-o —a
— *
eeeeeeee
i i i i i i i i
uuuuuuuuu
e — rso
-------
S3»
OS
f 1 1
aau
1
eoaaoeo
i i i i i i t
U U U U bl Cd U
r3rse>o o oci
f j o — o «i o a
8e»
oo
i i <
«uu
• 00
eot-
SNN
fl —
I
f I I I f 1 I
uuuuuuu
.«e>oo — oo
o o a o o e- o
- * ei ci - cj —
S J
uu
o r- ci n N a
ooo ooo
+ 11 iii
uuu uuu
ecio
— ei —
OOOOOOO
I I I I I
! — — ONCIO
N a » ei n - ei
u o
_ U
e-i
UU
or-ci
000
nan
oeo
UU U
OOO
n
•» ^ s- o cira t
i3O«or»r-o
••a ci ci si * c> —
n§
g=§=
sS
ooo
uui
823
ooo
iii
a*W
f. — a
one
-Q-
NOTI
ooo
i i i
uuu
n a o *
oooo
I I I -f
uuuu
?§C10
r« — n ci
— ci n ci
0930
n f- t- ci
n t. -.• n
rjr.no
o ooo
1 1 1 1
-UUU
oa-aa
-an GO
nci — o
oano
o ooo
•TO CIO
oeeio
anoo
UblU
O>— f.
030CI
-00*
ooo
N»C1
C1S-
ON*
eeo
— ri o — — v
enan
ooae
uuuu
— f -ci
e ooo
•' •* •'• u
n 3 NO
S-T3-
ci ci n ci
Q
0-00
ncir.0
i
floor*
ooeo
i i i i
uuuu
?,??*-
OOOQ
till
hlUUbl
^ rt C3 rt
Nar-a-
uuuu
cin cie
NO nci
cm — a>
i i r
I
I I I I I f I
UUUUUUU
n o ci o c = n
o u t- ci o e- a
o-a —nod —
•c r» * f o o n
o oooooo
I I 1 I I I I
II
sis
000
I I I
uuu
ooo
•f-CI
ooo
iii
NNCI
oo —
O 13 f
— O-&
ooo
I I I
uuu
OO 3
ci N o r* N r- o
C O — •!• CJ — «•
NO*
OOO
I I I
uuu
nno
ont-
uuu
»« o
o n tf>
— — cs
i i
• - — T -. s •-. is a
^o«j-'-cio-«> 2s
o N o « « o e>
ooooooo
i i i i i i i
uuuuuuu
ci f- a ci
o - o o c e> o
o 3 c a ci o
r. a n to r> ci ti
n«3cccc
.
uuuuuuuu
o r. e> •* T ct ci a
ci n o CT n •*« »
sooeoooo
i i i i i t i i
u u u u u u u u
Q —
ci •» o n « •» n n
oni3« cici «• o
'/> ei --•»-
1
2 ei c t a r r- —
D-9
-------
fll»C
ee —
i i i
». —r. *|»'«
an* —«««•
•eeeeeo
I I I i I i i
UkJUUUUU
o » e c »
»ee
e ——
i i i
ri • n « » n n
B e -
- us;
•HE enc
c :i» ? si is f
t» * — —!• — f
8'eS SiS
II II
rtctn
SSS KKt-
•>eo
nS-c*-o-
or-n ne>n
eae eee
+ i i itt
UUU UUU
can o — e-
isi.0 nea
eta» atio
i i i i i i i
uuuuuuu
2§S 8SS
ss:
•>Q
ee
«B«S- -O-
eeeee
uuuuwuu
• — * —
UUU
non oeo
>— cio a
BBC
SS
poo one
:eo -•»«•
e»— aeo
oeo eee
tit iii
uuu uuu
no— »s:i
—f I-
»£§!
nan
eeee eea
i i i * iii
uuuu uuu
e> ci t» ri » n ci
n-ee Sea
»ri —«j -ti-e
i i
eeee eeee
iiii i i i i
uuuu uuuu
n — tt e e c ci n
VNOO onNO
BO — e oocio
1
siis
e
i •* i i
uuuu
§t7»e>
-cio
O — !• —
eeee
i i t i
uuuu
— OO
no
pe
man—
CO-J1C6
IO-3CCI
l-0t300»
CO-3V9C
oooo
uuuu
goae
t
n —
eee
wuu
eeeeeee
i i t t i i i
uuuuuuu
——-CIO
ooooooo
I I I I 1 I t
uuuuuuu
eeeeeee
i i i i i i i
uuuuuuu
i i i i i i i
uuuuuuu
- »i. n o n o
U-«O
eeeeeeee
i i i i i i i *
uuuuuuuu
eeeeeeee
i i i i i i i »
uuuuwuu
S
oeoeooo
HUUUUUUU
uuuuuuu
c
e-oneoeov
-n — ««— co
eeeooeeo
UUUUUUUU
— laoco
——tsee
oeeeeeee
— i i i i i i i i
uuuiuuuuu
— fjo»n— cs
D-10
-------
N
Cl
1
u
-OBCI
tfi
LJ
2
••
o ci
B
O
M
Cl
i
5
if
s
H
o
n
0
o
o
I
>
O.O.C
D-ll
-------
APPENDIX E
EXAMPLE USE OF VARIABLE
PHOTOLYSIS AND PRINT OPTIONS
E-l
-------
APPENDIX E
EXAMPLE USE OF VARIABLE
PHOTOLYSIS AND PRINT OPTIONS
Ml 14 *O T 1
I |2 IT 18 I* 24 29
.•16 .13 .29 .39 .46 .8 .82
.93 .83 .9 .46 .3* .29 .13
01*
mn i no o i.
O 2 03, 4400000.
.3 HO 3 H02 23.8 I4SO.
hO2 03 4 H03 .0474 2430.
H«J2 0 8 HO 13000.
M3 HO 63. H02 20000.
»02 M3 7 H20S 3910. -861.
CSOS 8 H02 H03 6.03 10320.
k^OS 92. BK03 .02
hi) R02 102. BOHO .0000003
HOMO UOHO 11 HO BOS .000018
noao 12 on MO . 10
OB B02 13 IBI03 16000.
OH HO 14 BONO 160OO.
002 HO 13 R02 OU 12000.
DOS 1102 16 BOOB 7400.
UOOB 172. OU* .00079
03 IB 010 .002?
03 I* 0 .069
OlD 2O 0 4.4SE10 -97.4
010 212. On 6.8E09
OB 03 22 BOS 76.9 1000.
03 802 23 OB 2.03 SBO.
BCUO 24 .0046
UCaO 2S2. U02 .003
1KUO OB 26 U02 14000.
HE02 HO 2O HEO H02 12000.
MEO 29 UCIIO 002 200000.
HCO H02 3fl 18000.
HEO «02 31 BCUO BONO 4400.
HC02 B02 32 74OO.
REACTIVITY 1000 fPMC HETBAHE
3 60. 60.
HO M2 C04
.1 .1 1000.
840. 303. .01
fLOT
IUOCTIVITY 1*00 erne KETILXHE
I 0.00 0.40 0.80 1.20 1.40 0. 1.6 0. 880.
03 0
BLANK CAW)
BLAUK CARD
E-2
-------
^,
i
u
s
n
0
u
e
n
»»»«l. OS'S — fl
n«f.O~--fio*'^»f"!ST-*— f»
— — — • ri ci fi ei N ct N ri M :i n n n
E-3
-------
B
£
•.
S
0
f
»
I
S
•
^
•
i
» u
§ i
"
o
§u
1
-
e
»
n
f
•
M
Eg
5
1
o
n
e
i
S
C9
•
9
y
U rf
U ul
c S
!
0
I
ei
••
E
a
S
«•
n
i
•
•
;
i
ei
i
t>
8
u
IS
•
e
ui
e>
t»
0
ei
•
e
A
&
o
»
»
•
•V
«
u
«
N
V
«
I
»
«
8
u
e
ci
:
g
•
K
n
E
e
e
»
»•
f
S
*
u
o
«
»
e
+
u
e
e
<9
"
o
••
S
e
e
0
|
e*
I
•»
i
»
e
i
e
»
rt
I
«
CJ
e
ill
e
«
»
S
•^
u
s3
fl
O
«
5 S
4> 4>
2 u
S 0
ri e
— »
a N
» s
»
N
»
i
2 S
|
i
E
?
^ ^
u
g
Ik
e
S3
3 S
< W U
Sfc.
e e
= S
8
i
e
ca
8
J>
e
o
«
n
e
u
e
o
N
ei
*
e
^
1
r»
_
o
%
••
P
i
H
O
*
1
X
_l
<
>
1
s
2
3
S
E
u
e
e
e
o
i
s
0
i
0
o
1
hi
e
e
M
0
O
i
o
I
»
e
u
i
V
C)
o
u
0>
ti
e
i
u
e
S
rt
|
0
i
u
^J
e
u
§
ri
e
u
e
e
fe g
n
o
§n
e
ul
5
e
3
3
^*
S
e>
£
•j
S
a
3
V
e
2
u
5
u u u
o s e
* » «
u
T
•T
2
»
w)
i*»
?~
£
55
-J
P
d
_ *
•ri
• '
3
H
^
r™
u
p*
o
»
u
o
9
O
v
ER
5
:—
i^
rt
i
y
P
Cl
a
i .
u
e
0
^
~
o
i
a
«
T
•o
0)
E-4
-------
• <8 >3 K C- •C
O O O O O O
I I I f I I
w u u u u u
ft - «• t- 0> fl fl
o » »• * N » o
_: « ei « o « n
o o « a a e r-
— oooooo
U U hi U U U U
O -»««> — N
e * — — — o » »
O « O * * -a * o
§or> nn en *n »n »n »n »»
oo> oe oo oo oe oo eo oo
4- I + 14- I * 14- 14- 14- I4> 14-
O u uu uu uu uu uu uu uu uu
e o fo »o no »o —o no —o o —
O_ o -to -so oo — o no fie fia - ei- •*— o— »— n-
8 14^ 14- 14- (4> 14- I4> 14- 14-
no, u uu uu uu uu uu uu uu uu
o2 o cio ao oo o>o oo cia -a ao
EU n an «n «>n «n vn an en nn
Bt- o ns a>o ao eo DO oo t-o »o
on nn en cin — n no —a fin an
N— »•— t-o i>-o* >e» w» »o> csa
o— a*- a— a— oe oa oa ee
ii it ii it ii it ti ii ^~^
n uu uu uu uu uu uu uu uu __
oo fti ao no Nn e> r- no «o> w "-^
ci y vri «n o a- nn —»- «r- «a »» — — C
i"
4->
«e> CO aia or- at- f-o oo nn c
oo oo oo oo oo oo oo ao Q
uu uu uu uu uu uu uu uu O
on no> TO
as -a -ra -oo
ou (if s> f i aa —a e>a — e>
EC o « * f i o f i — •» no fio
oe a— —— -» n— o'» w-
n nn on tin «n fin tin -ft -fi ?
o oo oa oa oe oo 30 oo ee v*
4 14- 14- 14- t 4- 14- 14- 14- 14- GJ
u u u uU Uu uu uu uu uu uiu f&
T o — o •* o no 90 oo c^^ f 19 ^ t*
n= a no ee TO oo ao ae> Oe> *«
ou o no fie fie oo QO no —e> «a> 4J
» - - * » » ^
e>« e* »e> an 5 ft Sfi a— a— 3
oo oe oa oo oo oo oo oe O
it ii it it ii ii ii ii
o uu uU Uu uU uu uu Uu f.i ^j
a n-r nc> «o «\o »•» oo on fio
o = oc^ on on oe cl— ao ao oo
oo fin *» a'fi -N -fi — r*n aa «o on
z= o nci se sr «-r — ci *N o-r rs —
— o' e« c i » fl * n » * — n — N — f I — —
E-5
-------
t
u
e
n
n
i-
, i
s
M
S
s
I
i
to
e
1.
ss
1 +
uu
ss
uu
II
S3
I +
is
?!
uu
51
88
i ^>
uu
si
88
i +
uu
• a
S3
sfe si
S3
ss
uu
oo
t +
uu
ao
»o
ao
-e
as
-e
uu
S3
ss
n«
00
I. t
oo
I t
a*
oo
si
2 S-
uu
t»e<
uu
TJ
a>
ss
uu
ss
uu
oo
ii
*<•
So
ii
ee
SS
§
—n
oo
uu
oa
a»
88
uu
a v
on
eo
uu
—o
*•
(I«
ss
uu
S"
ei«
SS
MM
S?
88
uu
gs
(A
O UU
S Sc~
>= ar>
10-3606
00-3610
So
uu
SS
uu
an
s;
• i
uu
uu
e»o
aei
act
uu
25 i?
-w
uu
00
uZ
an
o-
iu
S3
n-
82
i i
uu
oo
§3
uu
oo --
PIFI
eo
= uu
net
00
uu
s
-------
o
s
e
§
S
o
§
o
o
o
f.
o>
c
0
in
&
— o
ii
O
o
fl
N
O
un
o
ei
E-7
-------
* *
*
•o
•o
o
o
O.
4J
O
E-8
-------
TECHNICAL REPORT DATA
(I'lcasc read Instructions on th.L reverse before completing)
1 REPORT NO.
EPA-600/3-30-028b
2.
TITLE AND SUBTITLE
MODELING OF SIMULATED PHOTOCHEMICAL SMOG WITH KINETIC
MECHANISMS Volume 2. CHEMK: A Computer Modeling
Scheme for Chemical Kinetics
5 REPORT DATE
February 1980
6. PERFORMING ORGANIZATION CODE
3. RECIPIENT'S ACCESSION-NO.
7. AUTHOR(S)
G. Z. Whitten and H. Hogo
8. PERFORMING ORGANIZATION REPORT NO.
EF78-107R
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Systems Applications,
950 Northgate Drive
San rafael, California
Incorporated
94903
10. PROGRAM ELEMENT NO.
1AA603 AC-054 (FY-79)
11. CONTRACT/GRANT NO.
Contract No. 68-02-2428
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory-RTF, NC
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, N.C. 27711
13. TYPE OF REPORT AND PERIOD COVERED
Final 7/78-9/79
14. SPONSORING AGENCY CODE
EPA/600/09
15. SUPPLEMENTARY NOTES
16. ABSTRACT
Mechanisms that describe the formation of photochemical smog are developed
using a computer modeling technique directed toward the simulation of data collected
in two smog chambers: an indoor chamber and a dual outdoor chamber. The results
of simulating 164 different experiments are presented in Vol. 1. Individual com-
pounds for which specific experiments were simulated and mechanisms developed
include the following: formaldehyde, acetaldehyde, ethylene, propylene, butane,
and toluene. Experiments in both chambers were simulated for all these compounds.
The mechanisms reported describe the decay of the precursor organic compound,
formation and decay of secondary organic compounds, conversion of nitrogen oxides,
formation of nitrates, and the appearance and decay of ozone. Special emphasis is
given to the chemistry of toluene. Also included is a study of a generalized
smog-based or carbon-bond mechanism developed in a previous study. Volume 2
contains the user's manual and coding for a chemical kinetics computer program,
CHEMK.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
c. COSATl Field/Group
*
*
*
*
*
*
Air pollution
Reaction kinetics
Photochemical reactions
Test chambers
Mathematical models
Computerized simulation
13B
07D
07E
14B
12A
09B
13. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (ThisReport!
UNCLASSIFIED
21. NO. OF PAGES
113
JS (This page)
22. PRICE
EPA Form 2220-1 (9-73)
------- |