EPA-R4-73-013b
CONTROLLED EVALUATION OF THE
REACTIVE ENVIRONMENTAL
SIMULATION MODEL (REM)
Volume II
User's Guide
Pacific Environmental Services, INC
2932 Wilshire Boulevard, Suite 202, Santa Monica, California 90403
-------
CONTROLLED EVALUATION OF THE REACTIVE
ENVIRONMENTAL SIMULATION MODEL (REM)
VOLUME II: USER'S GUIDE
EPA R4-73-013b
by
Allan Kokin, Lowell G. Wayne and Melvin Weisburd
Pacific Environmental Services, Inc.
2932 Wilshire Boulevard (Suite 202)
Santa Monica, California 90403
(213) 829-3658
Contract No. 68-02-0345
Project Officer: Ralph C. Sklarew
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND MONITORING
ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
February 1973
-------
This report has been reviewed by the Environmental
Protection Agency and approved for publication.
Approval does not signify that the contents
necessarily reflect the views and policies of the
Agency, nor does mention of trade names or commer-
cial products constitute endorsement or recommen-
dation for use.
-------
TABLE OF CONTENTS
SECTION I. INTRODUCTION I.I
SECTION II. DESCRIPTION OF THE REACTIVE ENVIRONMENTAL
SIMULATION MODEL (REM) II.1
A. Introduction II. 1
B. Grid System II. 4
C. Program Code and Machine Use II.4
D. Operation of the REM II. 4
E. Control Strategies II.6
F. REM Modules II. 7
1. Initialization Program II. 7
2. Numerical Integration II. 8
3. Chemical Kinetics II. 8
4. Meteorology II. 9
5. Barrier Check 11.10
6. Ultraviolet Module 11.11
7. Source Emissions 11.13
8. Output Routine 11.13
9. Plotting Routine 11.22
10. Reverse Trajectory Program 11.22
SECTION II. MATHEMATICAL DESCRIPTION OF THE REACTIVE
POLLUTION SIMULATION MODEL (REM) III.l
A. Meteorological Module III.l
B. Ultraviolet Irradiance Module III. 4
C. Kinetics Module Equations III. 7
D. Emission Module III.11
E. Numerical Integration III.12
SECTION IV. FLOWCHARTS IV. 1
SECTION V. MODEL INPUTS V.I
SECTION VI. SAMPLE INPUT DECK FOR REM VI. 1
SECTION VII * GLOSSARY VII. 1
A. Common Regions VII. 1
-------
ii
B. MAIN VII. 3
C. SMGOUT VII. 3
D. DIFFUN VII. 4
E. METEOR VII. 4
F. ULTRAV VII. 5
G. SOURCE VII. 6
SECTION VIII. PROGRAM LISTING OF REACTIVE POLLUTANT
ENVIRONMENTAL SIMULATION MODEL (REM) VIII.1
SECTION IX. OPERATION OF THE REVERSE TRAJECTORY PROGRAM IX.1
SECTION X. REVERSE TRAJECTORY INPUTS X.I
SECTION XI. REVERSE TRAJECTORY SAMPLE INPUT DECK XI. 1
SECTION XII. REVERSE TRAJECTORY PROGRAM LISTING XII. 1
-------
I.I
Section I
Introduction
This document describes the design, operation and use of the
Reactive Pollution Environmental Simulation Model (REM). Its main
purpose is to explain to the user how REM functions and from experience
gained in operation, how it can best be utilized.
The latter includes imparting to the user the following:
• Detailed instructions for setting up input decks;
• How REM may be utilized in performing control strategies;
• An understanding of the methodology, techniques and
calculational procedures employed in REM; and
• Sufficient information to enable program modifications
to be accomplished.
The above tasks are facilitated by the incorporation of algorithm
descriptions, flow charts, complete data input instructions, sample
input decks, a glossary of program terms and program listings in this
volume.
-------
II.1
Section II
Description of the Reactive
Environmental Simulation Model (REM)
A. INTRODUCTION
The REM effort has been directed towards the construction of a
sophisticated photochemical mechanism which is capable of closeJv
replicating smog chamber data and then extending this capability to
the urban atmosphere. REM is based on a Lagrangian (moving coordinate)
system. This system enables the numerical computation of the chemical
reactions to take place in a moving parcel or column of air. The
parcel is bounded by the inversion base above and the ground below.
The locations of the base of the column at successive points in timt
describe the path or trajectory that the air parcel traverses across
the region. The trajectory is computed by special routines contained
in the REM program from wind velocity information contained in the data
base as a function of time of day and location. An overview of the
model's operation is presented in Figure II-1.
The computer program has been designed as a detailed episode mode;
for chemical kinetics, with additional modules linked up to form a
complete atmospheric system. The modules presently in the system deter
mine the necessary meteorological parameters, the rate of absorption
of actinic light by N02, and emissions due to traffic and area sources.
A schematic diagram of REM appears in Figure II-2. The facility exists
to incorporate modules to calculate the emissions due to large point
sources and the effect of vertical diffusion should practical and
effectual methods be found to perform these tasks.
-------
II.2
PHOTOCHEMICAL MODELING
For pollution modeling, the air space within each
geographic grid is assumed to emit an approximate
mean value of pollutants for a given time. The
model includes nitrogen oxides, high and low re-
active hydrocarbons and carbon monoxide. Values
are provided by local environmental control
agencies.
Meteorological stations provide
current data on winds, tempera-
ture, pressure, dew point, etc.,
to the model for predictive
movement of air parcels through
the various grids.
Pollutant monitoring stations
measure types and amounts of
pollutants, i.e., oxidants, hydro-
carbons, nitrogen oxides, car-
bon monoxide.
TERRAIN
INVERSION LAYER
The column of air
is one square meter
in cross-sectional
area. Its height is
determined by ter-
rain and inversion
layer conditions.
The latter are de-
rived directly from
radiosonde meas-
urements.
The model predicts the amount and
type of photochemicals present in a
column of air as it traverses a com-
puted trajectory path to a pollutant
monitoring station. At this point
predicted and observed values are
compared.
This model simulates a region's photochemical pollutant content
by representing in a computer the chemical reactions occurring
in a column of air as it traverses a dynamic wind trajectory. The
model is tested by acquiring appropriate meteorological and
pollutant data from a dense network of monitoring stations.
The path of a moving parcel of air is computed on the basis of
wind station data. As the parcel passes over the urban area it
picks up pollutants from vehicular traffic, power plants, factories
and other sources. Pollutants, or emission rates, for each grid are
contained in the data base. With the sun's energy acting as the
catalyst, the pollutants are transformed into a variety of photo-
chemical pollutants, producing the well known effects of smog.
Predicted concentrations for these pollutants are printed out for
points along the air parcel path and compared with values meas-
ured at actual air monitoring stations.
COMPUTER
In a given area this photo-
chemical model could predict the effects of
new smog control devices for any given time.
period after introduction, determine a pro-
posed freeway's effect on pollutant concen-
trations, or the contributory effects of new
industries and proposed rapid transit sys-
tems. Decreases in pollutants that would'
result if modified or advanced propulsion
systems are used in vehicles or if nuclear
power plants are installed also could be
predicted by the model.
Figure II-l. PES Reactive Environmental Model (REM)
-------
II. 3
WIND. TEMPERATURE. HUMIDITY
EMISSIONS INVENTORY, RADIO-
SONDE CURVED SOLAR
CURVES, EMiSSiON~F ACTORS
DETERMINATION
OF EMISSIONS
INPUT INTO COLUMN
METEOROLOGICAL
EVALUATION
SOLAR
RADIATION
MIXING
EPTH, TEM
HUMIDITY
WIND
ABSORPTION
RATE
EMISSIONS
RATES
RATES OF
CHANGE OF
CONCENTRA-
TIONS AND
LOCATION
NUMERICAL
INTEGRATION
(ADAMS)
NEW TIME
AND
LOCATIO
NEW
CONCEN
TRATIONS
SMGOUT&SMGPLT
OUTPUT
TABLES & GRAPHS
Figure II.2.
Schematic of the PES Reactive
Environmental Model (REM)
-------
II.4
B. GRID SYSTEM
The geometric setting for REM is a standard two-dimensionsl coor-
dinate system. REM utilizes the mile as the standard unit of measure
with all urban grid squares chosen as 4-square miles. The kilometer
can also be used as the standard.
To set up the grid system, the origin is placed at a point such
that the region of interest resides entirely in the positive quadrant.
The coordinate positions of all air monitoring and meteorological sta-
tions are established in either miles or kilometers. These values are
then fed into the model at execution time. In this manner, REM can be
applied in any region.
Figure II-3 illustrates the grid system used by REM for the Los
Angeles Basin. The locations of the air monitoring stations are desig-
nated as circles and the coordinate points are given in miles.
C. PROGRAM CODE AND MACHINE USE
REM is coded in the FORTRAN IV computer language and has operated
with both the Level G and H compilers. It is thus operational on any
computer system with 150K in core and a standard FORTRAN compiler.
The Los Angeles validation version of the model has been success-
fully run on the IBM System 360, Models 50 and 65, and the IBM System
370, Models 155 and 165. The ratio of simulation time to computer time
used was approximately 150:1 on the 370/155 and 550:1 on the 370/165.
In regions where much less data is available than in Los Angeles, REM
will operate at even greater speeds.
D. OPERATION OF REM
REM has been designed to operate as a multi-receptor point model.
The user is required to provide a minimum number of inputs after the
data base has been prepared. In its execution mode, the model trans-
ports a parcel of air, bounded by the ground and the inversion layer,
-------
II.5
53
SI
49
47
45
»7
33
33
31
39
27
15
23
21
19
17
13
13
11
>LA
EUI
o
*:
\
17 19 21 23 23 27 29 31 33 33 37 39 41 43 43 47 49 31 53 55 57 59 61 63 65 67
Figure II.3. REM Grid System for the Los Angeles Basin
-------
II.6
along a dynamic trajectory that is generated utilizing wind data from
numerous meteorological stations. At each point along the path, the
speed and direction of the trajectory are recalculated.
The user may select the starting point or the terminus for each
simulation run. In either case, he must indicate the length and starting
time (of day) of the simulation, and initial values for the major pollu-
tants.
If the user selects the terminus, he must first operate the reverse
trajectory program (see Section F.10) to determine the correct starting
position and major pollutant values. The starting location is necessary
because this point must be input into REM in order for the trajectory
to reach the selected terminus at the end of the full simulation run.
The initial pollutant concentrations are determined by interpolation of
observed values at air monitoring stations to the trajectory starting
point. Once the location and starting values have been calculated for
a simulation run, the same trajectory may be run many times over with-
out utilizing the reverse trajectory program again.
If the user executes the model by selecting a starting point, he
may wish to choose a monitoring station where the initial pollutant
values will be known. As an alternative, he may select another point of
interest, or pick a completely random location. For these latter cases,
some assumptions may be made regarding the starting values, or minimum
or standard background values may be input.
CONTROL STRATEGIES
After a given full simulation run has been made, it may be consi-
dered as a standard run for the purposes of evaluating particular con-
trol strategies. Control strategies are implemented to reduce pollu-
tant concentrations and to attain required air quality standards.
-------
II.7
By varying any of the numerous emissions inputs to the REM, either
individually or in combination, the possible effects of control strategies
can be determined before they are actually implemented. These strategies
may include for example: reducing total traffic, reducing traffic in
specified areas, reducing emission rates, varying traffic patterns, or
altering operating times of area sources.
F. REM MODULES
The modular construction of REM has facilitated the model's
development and alteration, when necessary. Any given module may be
unplugged and replaced as desired; or special modules may be created
and added to handle special situations.
1. Initialization Program
The initialization program reads in the data base and the infor-
mation needed for a complete simulation run. This program converts
all of the data, except station locations, to the metric system for
easy utilization by the other modules, and calls the numerical integra-
tion, output and plotting routines as required by monitoring the simu-
lation time. At the end of the run, it checks to see if another simu-
lation is to take place. To make additional runs, only control informa-
tion must be input thereby facilitating the operation of REM for the
user.
The initialization routine outputs the data base that has been
read in to provide the user with a permanent record of the information
employed in the simulation. A complete list of the data categories
required for input follows:
• Trajectory start position and time of day;
• Length of simulation;
• Starting pollutant concentration values;
-------
II. 8
• Control parameters for numerical integration routine;
• Meteorological stations - location and elevation;
• Meteorological variables - hourly data (wind speed and
direction, temperature, dew point);
• Emission factors for freeway and street traffic in grams
per vehicle mile;
• Grid distribution of hourly area emissions;
• Radiosonde data - temperature vs. elevation profile; and
• Solar zenith angle as a function of time.
2. Numerical Integration
The numerical integration module is the heart of REM. This routine
accepts as input derivatives calculated by the kinetics and meteoro-
logical modules and computes pollutant concentrations and the parcel
position forward in simulated time. The kinetics module supplies the
derivative of each of the photochemical pollutants computed as a
function of time. The meteorological module provides wind velocity
outputs as the time derivatives of the position vectors.
The economic use of REM is dependent upon the numerical integration
routine. By attaining larger step sizes, the program will execute
faster, use less computer time, and be less expensive to operate.
REM utilizes the Adams method of integrating ordinary differential
equations of the Gear numerical integration package. A sensitivity
analysis has been conducted on the Adams routine in REM. As a result
of this examination, it was found that the Adams method executes most
economically with the maximum step size at 0.1 minutes and the error
limit at .05 in pollutant concentration units.
3. Chemical Kinetics
The chemical kinetics module has operated on reaction mechanisms
C. W. Gear, DIFSUB for Solution of Ordinary Differential Equations,
Communications of the ACM, Volume 14, Number 3, March 1971.
-------
II.9
containing 13-41 reactions and is based on a set of stoichiometrically
valid elementary reactions. The current version contains a 32-reaction
mechanism that preserves the flexibility of chemical detail necessary
for adequate simulation of both primary and secondary contaminants.
The set of elementary reactions, rate constants and chemical species
is easily modified and expanded as required.
This routine calculates the elementary reaction rates and the
derivatives of each of the accumulating chemical species as they change
with time. It calls the meteorological, ultraviolet and emissions
modules into operation, as necessary. The kinetics program adds the
perturbations from the traffic and area sources to the derivatives
computed for the chemistry for input into the numerical integration
routine.
4. Meteorology
This module consists of a series of related routines that employ
wind speed and velocity, temperature, relative humidity, elevation,
and radiosonde data to calculate the meteorological values, mixing
depth and rates of change in the wind at the pollutant column position.
a. Spatial Interpolation
Wind speed and direction, temperature, and relative humidity
undergo both spatial and temporal interpolations from numerous
meteorological stations to the trajectory location of the air
parcel. The spatial interpolation is inversely proportional to
the square of the distance from station to parcel. This procedure
has two advantages over simply using a ratio inversely propor-
tional to the distance involved. Firstly, it ensures that the
meteorological values assigned to the trajectory position are
more closely in keeping with the data at the nearest stations; and
secondly, it permits the computer program to execute faster by
eliminating the necessity of computing square roots.
-------
11.10
b. Temporal Interpolation
The temporal interpolation is linear. The meteorological
items for the current hour and the next hour are utilized in this
process. The portion of each item contributing to the result is
directly proportional to the number of minutes in the hour that
have passed. For example, at 10:15 A.M. one-quarter of the value
for hour 10 would be added to three-quarters of the value for hour 11.
c. Mixing Depth
The model calculates mixing depth on the basis of a temperature
versus height profile derived from temperature soundings at various
altitudes aloft recorded by radiosonde measurement or aircraft.
The mixing height is taken as a function of temperature at the
ground and the most recent vertical temperature soundings. The
mixing depth is then obtained by subtracting ground elevation at
the location of interest. The ground temperature is interpolated
at the trajectory point from meteorological data at each station.
In this technique the temperature versus height profile is assumed
fixed during the day unless more frequent soundings are provided.
The mixing depth changes with variations in ground temperature and
elevation.
d. Trajectory Derivatives
The rates of change of the wind components at the air parcel
position are utilized by the numerical integration routine to
calculate the next trajectory position. Two differential equations
are reserved for the trajectory calculation. They employ these
derivatives in the same manner in which the integration of the
differential equations for the chemistry occurs.
-------
11.11
5. Barrier Check
A "barrier" is an imaginary line or set of connected line segments
superimposed on the region of interest with the coordinate system for
that area. Its purpose is to represent a topographical feature such as
a mountain range which may serve to separate sections of an air quality
region. The mountain range may tend to influence differing climato-
logical characteristics and wind regimes. For example, on a given day
in the Los Angeles Basin, the meteorological characteristics of the
Westwood area are likely to differ from those in Van Nuys due, in part,
to the intervention of the Santa Monica Mountains. Experience in the
Los Angeles area has shown that significant temperature and wind regime
differences between the San Fernando Valley and the central Basin tend
to occur on any given day.
REM contains two such barriers. One bisects the Santa Monica
Mountains and separates the San Fernando Valley from the main section
of Los Angeles. The other barrier passes through the Puente and San
Jose Hills at the east side of the Basin.
The main function of the barrier routine is to determine the
"permissible" meteorological stations. Permissible stations are simply
stations which are located in the same general section of the Basin as
the trajectory point of interest, as defined by the barriers. That
is, both the trajectory point and all stations to be used in the cal-
culation of the trajectory path and the interpolated meteorological
values at any trajectory point are located on the same side of each
barrier. This is graphically shown in Figure II.4.
The barrier check is intended to ensure that a station which is
very close in distance to a trajectory point of interest, but separated
from the trajectory point by a mountain range does not exert any
influence on the calculation of the trajectory path, temperature, and
relative humidity at the trajectory point.
-------
M
2
•
(Wind Direction)
M
Wind direction
Line T - T. is the trajectory path
a D
T, - T trajectory point of interest
' a
T9 - T downwind trajectory point along the same trajectory
c. d
P-, Permissible meteorological stations
for T, (M4, M5, M6> M?, Mg, Mg)
Permissible meteorological stations
for T2 (M2, M3,
, Mg,
, Mg)
M. meteorological stations
Figure II.4 Demonstration of Barrier Check Technique
-------
11.13
6. Ultraviolet Module
This module calculates a diurnal ultraviolet irradiance function
based on measurement of sky cover, latitude and local calendar time.
The irradiance determines the specific ultraviolet absorption rate of
NO-, which is used in the kinetics module.
7. Source Emissions
This routine computes emissions of NO, CO, reactive hydrocarbons
and less reactive hydrocarbons for freeway and street traffic, and
calculates emissions from area sources of NO, reactive hydrocarbons,
and less reactive hydrocarbons.
The inputs to this program for vehicular sources include two
2
diurnal traffic curves (see Figures II.5 and II.6), emission factors
for reactive hydrocarbons, less reactive hydrocarbons, CO and NO in
grams per vehicle mile, and distributions of daily vehicle traffic on
streets and freeways given in a system of 650 4-square mile grids that
2
cover the Los Angeles Basin (see Figures II.7 and II.8). The emissions
of the above pollutants are computed as a function of time by applying
their respective emission factors to the traffic levels at the air
parcel location qualified by the daily traffic usage curve.
Emissions from area sources are also calculated as a function of
time. A constant diurnal usage curve is applied to grid system dis-
3
tributions of hourly pollutant emissions of NO, reactive hydrocarbons,
2
and less reactive hydrocarbons (see Figures II.9, 11.10 and 11.11).
2
Appendix A, Development of a Simulation Model for Estimating Ground
Level Concentrations of Photochemical Pollutants by P. J. W. Roberts,
P. M. Roth and C. L. Nelson; Report 71SAI-6, March 1971, prepared by
Systems Applications, Inc.
3
Oxides of nitrogen taken as NO
-------
2
O
U
O
O
U
.09
.08
.07
.06
.05
.04
.03
.02
.01
M
\ 1
1 1 1 1
42468
1
10
I I
12 14
1 1 1 1
16 18 20 22 2
DAILY HOUR
Figure II.5. Temporal Distribution of Daily Street Traffic
(Interpreted from Data in Reference 2, pg. 11.13)
-------
z
3
O
o
O
Q
LL
O
z
o
H
O
<
cc
u.
.10
.09
.08
.07
.06
.05
.04
.03
.02
.01
Ul
24
10
12
14
16
18
DAILY HOUR
20
Figure II.6. Temporal Distribution of Daily Freeway Traffic
(Interpreted from Data in Reference 2, pg. 11.13)
22 24
-------
11.16
35
60
137
239
139
25
9
0
47
12
155
255
250
118
30
25
18
64
.-— -^.H
34
202
238
249
107
17
0
0
99
~*
149
223
238
250
194
12
2
7
117
N
60
146
252
266.
264
182
59
30
98
231
476
X
125
172
221
295
414
350
63
70
283
512
254
305
y
2N
50
135
152
251
328
268
121
32
403
346
338
255
266
261
U89
V
•
43
85
233
276
280
155
260
549
459
261
255
413
296
255
215
9\j 360
1 H239
IK 1 121
1 1
V
49
•~*
28
66
54
212
247
116
478
512
478
412
311
392
397
328
412
324
304
188
28
-------
11.17
59
169
271
.....
311
383
\
.__:
166
94
155
242
222
890
352
356
242
101
\
119
68
456
137
473
235
\
\
144
156
126
426
141
511
295
\
\
^
,
/
V
k
158
17
494
.140
333
51
160
404
348
351
I
)
•'•v_
103
104
178
275
215
426
32
391
^
150
395
194
205
466
326
^X
190
423
290
436
795
531
391
366
305
257
217
625
131
143
F
/
J
51
335
620
452
449
343
45
13
193
257
585
623
159
241
236
545
107
V
37
90
109
257
518
576
2B3
281
302
131
140
-•-.
28
297
224
332
129
56
323
306
79
444
527
314
383
278
195
529
X
311
234
205
265
33
30
289
102
349
\
41
254
285
263
104
243
59
213
181
\
74
362
147
156
253
83
144
112
153
\
66
291
70
12
330
183
78
115
X
65
214
59
251
145
25
150
160
80
V18
43
203
42
15
138
68
398
163
97
6
^
11
170
21
108
156
257
67
148
95
N
173
32
3
207
137
161
131
169
61
159
.!<
1!
16
85
.....
77
91
59
Figure II.8. Distribution of Freeway Traffic by Grid Area
(thousands of vehicle miles per day)
[From Reference 2, pg. 11.13]
-------
11.18
0
0
8
e
e
0
0
0
8
0
i
9
9
8
0
0
0
8
0
1
9
13
12
0
0
0
8
0
1
18
18
14
0
1
1
9
\J
9
10
1?
15
12
4
3
6
19
9
X
9
10
15
15
8
4
7
10
19
18
19
S'
\
9
12
15
17
8
4
9
10
19
20
21
25
29
V8
\
5\
0
/
V
^
4
9
.12
15
15
8
e
18
10
. 20
20
21
25
29
29
21
21
19
18
18
>-^v_
6
1
11
12
12
9
9
19
11
20
20
20
22
25
25
27
27
24
21
21
^
1
9
1
7
10
10
19
12
21
19
19
20
22
22
34
34
29
23
23
20
^
0
7
1
4
10
10
19
12
21
19
20
20
22
22
34
34
29
23
&
ft
0
2
8
8
9
9
19
12
21
37
38
31
24
24
25
25
30
35
^
0
0
4
8
1
9
19
12
21
37
38
31
24
24
25
25
30
35
^
0
0
1
8
9
18
18
11
20
29
29
26
24
24
24
24
26
28
-JD-
0
0
0
8
10
16
18
10
19
21
21
22
23
23
23
23
22
21
^
0
0
0
3
10
8
18
10
19
21
21
22
23
23
23
23
22
21
13
"Is
0
0
0
0
16
16
19
13
22
9
19
22
27
27
19
19
18
11
11
11
\
0
0
0
0
19
18
19
13
22
1
5
18
27
27
19
11
11
11
11
11
11
\
0
0
0
0
16
18
20
12
21
6
3
7
23
23
11
11
11
11
11
11
11
11
X
•
0
0
0
0
14
14
20
11
20
15
9
0
11
11
11
11
11
11
11
11
11
11
11
V
0
0
0
0
13
13
20
11
20
19
2
0
11
11
11
11
11
11
11
11
11
11
11
11
\^
*-^
0
0
0
0
9
9
17
9
1
18
1
0
5
11
11
11
11
11
11
11
0
11
11
11
^
0
0
0
0
9
9
17
9
1
18
1
0
0
8
11
11
11
11
11
11
0
0
0
0
0
V.
0
0
0
0
9
7
9
1
1
9
0
0
0
0
2
5
8
11
11
0
0
0
.0
0
0
0
0
0
0
8
0
9
10
10
8
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Figure II.9. Distribution of Reactive Hydrocarbon Emissions
(average kilograms/hour) by Grid Area
[From Reference 2, pg. 11.13]
-------
11.19
0
0
30
30
30
0
0
0
30
0
0
31
31
31
1
0
0
30
0
0
31
31
31
1
0
0
30
1
3
35
35
30
0
0
5
35
X
'
32
36
39
39
30
15
5
24
39
40
N
IT -"•
32
36
39
39
30
15
15
39
39
40
40
V
K
-.-
31
42
53
53
34
19
23
33
33
45
45
51
56
^
6\
1
/
Vs
10
31
42
53
53
34
34
33
33
33
45
45
51
56
56
31
31
'31
30
30
15
1
40
48
48
38
38
36
35
35
47
47
51
56
56
46
46
39
32
32
^
0
37
13
33
42
42
39
37
37
49
49
51
56
56
61
61
47
33
33
31
>
0
27
13
23
42
42
39
37
37
49
49
51
56
56
61
61
47
33
?f
LS
r
0
10
31
31
36
36
48
60
60
110
110
90
70
70
36
36
34
33
^
0
0
16
31
6
36
48
60
60
110
110
90
70
70
36
36
34
33
fyrf-
0
0
5
30
35
35
44
53
53
82
82
69
57
57
36
36
34
33
"~»
0
0
0
25
35
35
40
45
45
56
56
50
45
45
36
36
34
33
""A
0
0
0
15
35
35
40
45
45
56
56
50
45
45
36
36
34
33
14
X
0
0
0
0
37
37
49
62
62
34
34
35
36
36
34
34
32
11
11
11
\
0
.. -
0
1
0
0
37
37
49
62
62
4
19
35
36
36
34
IS
13
11
11
11
11
\
0
0
0
0
39
39
43
49
49
22
12
20
33
33
13
13
11
11
11
11
11
11
V
0
0
0
0
41
4i
38
36
36
39
34
4
11
11
11
11
11
11
11
11
11
11
11
\
0
0
0
0
41
41
38
36
36
39
9
4
11
11
11
11
11
11
11
11
11
11
11
10
^
*~—
0
0
0
0
30
30
31
31
1
32
2
0
5
11
11
11
11
11
11
11
0
11
11
11
"">
0
0
0
0
15
30
31
31
1
32
2
0
0
6
11
11
11
11
11
11
0
0
0
0
\°
0
0
0
0
0
20
32
4
4
1
1
0
0
0
2
5
10
11
11
0
0
0
. 0
0
0
0
0
0
0
0
10
30
30
30
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Figure 11.10. Distribution of Oxides of Nitrogen Emissions Taken
as Nitric Oxide (average kilograms/hour) by Grid Area
[From Reference 2, pg. 11.13]
-------
11.20
0
0
80
80
79
0
0
0
79
C=l
0
3
84
84
81
2
1
0
79
t=*=
0
3
84
84
81
2
1
0
79
ia_
6
8
105
105
81
2
11
20
99
X
91
108
125
125
80
41
21
79
118
120
\
10
'
91
108
12S
125
80
41
41
118
118
121
120
V90
'\
85
133
180
180
96
57
74
90
90
145
145
189
233
\214
\<
^
3
>
y
11
85
133
180
100
96
96
93
90
90
145
145
189
233
233
96
96
78
79
79
*•— •— -
26
3
101
158
158
113
113
105
98
93
150
150
179
214
214
173
173
130
96
96
^
0
107
56
116
129
129
118
106
106
156
156
175
193
193
251
251
182
111
111
95
10^
0
83
56
76
129
129
118
106
106
156
156
175
193
193
251
251
182
111
f?
/<•/
/26
0
22
81
81
103
103
153
204
204
426
426
333
245
245
131
131
144
156
Wj
•
0
2
52
81
24
103
153
204
204
426
426
333
245
245
131
131
144
156
J^T
0
1
11
80
102
102
136
174
174
307
307
267
197
197
120
120
122
125
"3r-
0
0
0
70
101
101
122
143
143
189
189
170
150
150
108
108
101
94
^
0
0
0
30
101
101
122
143
143
189
189
170
150
150
108
108
101
94
55
X
0
0
0
0
108
108
159
212
212
90
98
118
137
137
94
94
86
44
44
44
X
0
0
0
0
108
108
159
212
212
19
59
118
137
137
94
59
51
44
44
44
44
\
0
0
0
0
118
118
137
158
158
69
99
68
109
109
51
51
44
44
44
44
44
44
X
0
0
0
0
129
120
117
104
104
117
108
19
•44
44
44
44
44
44
44
44
44
44
44
\
0
0
0
0
129
120
117
104
104
117
38
19
44
44
44
44
44
44
44
44
44
44
44
V40
0
0
0
0
80
80
B2
85
6
86
7
1
23
45
44
44
44
44
44
44
0
44
44
44
^-
0
0
0
0
41
80
82
85
6
86
7
1
1
31
44
44
44
44
44
44
0
0
0
0
. 0
V
0
0
0
0
2
62
90
21
21
4
4
1
0
0
5
20
40
44
44
0
0
0
0
0
0
0
0
0
0
1
21
97
115
115
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Figure 11.11. Distribution of Less Reactive Hydrocarbon
Emissions (average kilograms/hour) by Grid Area
[From Reference 2, pg. 11.13]
-------
TIHE» 0.380199370 03 MINUTES _.STE_Pi..66.
_POSI TIO_N« 5.8.4039.. 40_. 3 573 STEP SIZE' 0.100000000 00
PAGE 13
_LXPEC.!EO_AAT_i_OF JLHANGJL_Of CONCENTRATIONS DUE TO EMISSIONS FROM REGIONAL SOURCES
PPHM/H1N CRAH-MOLES/MIN
PROPYLENE
NO
CO
OUHHC
PROPYLENE
NO
CO
OUHHC
STREET TRAFFIC
FREEHAY TRAFFIC
0.36600-02 0.10250-01 0.2196D 00 0.49180-01
0.14460-02 0.40490-02 0.86770-01 0.19440-01
AREA SOURCES
0.17520-02 0.15470-02
0.72520-01
0_._5_2 74E-06 0.1477E-05 0.3165E-04 0.7089E-05
0.2084E-06 Q.5836E-06 0.1251E-04 0.2801E-05
0.2525E-06 0.2230E-06 0.1045E-04
TEMPERATURE
•CENT.)
27.2137
conpouNo
1 N02
2 NO
3 03
4 C3H6
5 HN03
6 CH3CHO
7 HCHO
8 CO (PPMI
9 CH30N02
10 C2H3N05
11 FIXED 02
12 DUMHC
13
14 H
15 OH
16 H02
17 :H3
18 CH30
19 CH302
20 C2H30
21 C2H302
22 C2H4Q2
23 0
24 N03
25'C2H303
RELATIVE NIXING HIKING
HUMIDITY ELEVATION HEIGHT DEPTH
(METERS) (METERS* (METERSI
0.2871 181.3918 757.9106 576.5190
CONC PPHM CHANGE. PPHM/MIN
0.10570770 02
0.74260920 00
0.13083660 02
0.5829044D 00
0.29533170 Cl
0.71383390 01
0.34705030 01
0.4641978D 01
0.39845050 00
0.55473800 00
0.0
0.11756080 03
0.0
0.38162460-07
0.12446960-04
0.52582540-02
0.68931480-08
0.40259910-04
0.38551340-02
0.14202760-06
0.79538800-03
0.10959920-05
0.28038470-05
0.1H61o70-05
6.57995750-03
-0.15016660-01
-0.28831680-02
0.27850620-01
-0.48061620-02
0.11190490-01
-0.55822400-04
-0.70338840-02
-0. 38063320-02
-0.15815400-03
0.41342320-03
0.0
-0. 33399040-01
.. .
Figure 11.12. Sample REM Output
-------
11.22
8. Output Routine
This program prints output data at various times along the tra-
jectory route. The user may specify the times at which the printouts
are to occur and the routine checks to see if a requested output time
has been reached. A sample output is shown in Figure 11.12. The
primary data categories printed include:
• Time of output and trajectory position;
• Step size;
• Rate of change of concentrations from source routine
in pphm/minute and in gram-moles/minute;
• Meteorological values at the trajectory point including
mixing depth;
• Pollutant concentrations;
• Rates of change of accumulating (primary) pollutants; and
• Rates of the elementary reactions.
9. Plotting Routine
At the end of the simulation run, this program outputs three digi-
tal plots containing the history of each accumulating pollutant repre-
sented in time vs. concentration curves. An example of this output is
shown in Figure 11.13.
10. Reverse Trajectory Program
The reverse trajectory program calculates wind trajectories and
interpolates initial concentration values for use by the simulation
model. Its primary inputs consist of hourly average wind data moni-
tored at a network of meteorological stations, and hourly average
pollutant data recorded at a group of pollutant monitoring stations.
The coordinate positions of these stations are also required.
In addition, the user must specify a coordinate position and the
length of time for which the simulated trajectory is to be calculated.
The program then considers the given coordinate position to be the
-------
}.3901^_03_J1IN_SIMULAT_10N BS^O.90000 02 JfUflUI'O. 10000_OQ PRHTI4I-0.10000-01
»«CIF|C ENVIRONMENTAL SERVICES, INC. — PHOTOCHEMICAL MODEL IREMI — VERSION C8
F"» OEnCHSTRlTlq'i OCTOBER, 19T2 .
RIIK BFGINS AT 0700 HOURS ON SEPTEMBER 11, 1469
TRAJECTORY TO DOWNTOWN LOS ANGELS AT 1130 HOURS
6.2
7.9
IS."
23.2
30.9
19.5
46.2
53.9
61.6
69.2
76.9
84.6
92.2
99.T
- 115.2
122.9
130.6
138.2
145.9
193.6
161.3
163.9
176.6
184.3
Mil. 9
-' 199.6
207.3
214.9
22?. 6
230.3
237.9
245.6
253.3
261.0 "
268.6
276.1
284.^
291.6
299.1
307. ^
314.6
322.3
330.0
337.6
345.3
353.0
360.6
3<>9.J
176.0
• *•*.*.••*.*.••*.••
i 3 — - .. X4 " --.-.- .._.
> 3 241
3 241
> ' 741
> 3 2 4 1
.__,
532 4 1
532 4 1
5234 I
5234 1
! 2 4 3 1
5243 1
5 t 43 1
524 3 1
'
t 4 31
25 4 31
25 4 31
' 4 s 13
24 5 1 3
24 5 13
1°, 5 13
24 5 13
245 31
245 31
245 31
245 "3 I
2 4 5 3 '
0.20
a. oo
16.85
25.55
34.25
__
60.95
"" 7«I»0~
87.JO '
114.80
123.80
1".80
141.80
168.80
186. 60
222.80
231.80
240.80
L _ ..
276.35
285.13
293.90
302.90
329.90
338.49
34T.4J
356.45
M
•
N5
Figure 11.13. Digital Plot of Primary Pollutant Time vs. Concentration Histories
Vertical axis represents time in minutes; horizontal axis depicts pollutant concentrations in
pphm. The species plotted are as follows: 1 - NCK, 2 - NO, 3 - 0 , 4 - C H , and 5 - HNO .
-------
11.24
terminus of the trajectory, and computes consecutive points along the
trajectory in a direction opposite to that of the wind. At the con-
clusion of the trajectory (the determination of the starting point), an
interpolation routine calculates initial concentrations for nitrogen
dioxide, nitric oxide, ozone and carbon monoxide from the pollutant
data.
Although its primary use has been for calculating reverse trajec-
tories, the program may also be used for computing forward trajectories.
In this case, the pollutant interpolation routine would yield concen-
tration values which could be considered as validation points.
Theoretically, the program could operate with one meteorological
and one pollutant station. However, the greater the number of stations
utilized, the more realistic and accurate will be the results.
-------
III.l
Section III
Mathematical Description of the Reactive
Pollution Simulation Model (REM)
A. METEOROLOGICAL MODULE
1. Introduction
The meteorological module calculates weather parameters at the
air parcel location after considering the effect of natural barriers
such as mountain ranges. This routine also determines the mixing
depth and the rates of change of the wind components at this point.
2. Derivations
a. The Barrier Routine and Definition of Permissible Stations
A barrier is a set of adjoining line segments utilized by the
REM to represent a natural obstacle such as a mountain range. The
barrier program's function is to determine the set of "permissible"
meteorological stations employable by the meteorological routine
for interpolating wind speed and direction, temperature, relative
humidity and ground elevation. A station is permissible if and
only if it is located on the same side of a barrier as is the air
parcel. This is expressed mathematically below.
Each barrier is defined by an integer ND, which is one greater
O
than the number of barrier sides, and a sequence of points (x., y.),
i=l, 2, . . , ND, with each consecutive pair of points defining a
D
side of the barrier. Let (x , y ) and (x , y ) be, respectively,
o o s s
the coordinate of the trajectory point and a given station. Let
m = (y -y )/(x -x ), b = y -mx be the slope and y the intercept
SO SO O O
of the line joining the trajectory point to the station. For each
i = 1, 2, . . . , N_-l let
D
-------
III. 2
b. = y.-m.xi
y -m.x -b.
sill
ix o s' ^o Js'
y±+rmwb
-
A station located at (x , y ) is said to be permissible if,
s s
and only if, for each i=l, 2, . . . , N -1 either
U
\±*Q, X±21, M^SO or /i^sl.
b. Interpolation of Meteorological Variables
A spatial interpolation is performed automatically for the
data values from two successive hours. This is followed by a linear
temporal interpolation for the two spatially interpolated values
for each of four meteorological variables. The spatial interpola-
tion is given as follows: Let. (x., y.) be the coordinates of each
J J
permissible station (j=l,2,...N ). Let d.. .d-.d , . . . ,dN be the dis-
tance from these stations to the trajectory point and let
x..,x , ...,XJT be, for example, the corresponding set of observed x
p
components of wind velocity. The x component of wind velocity, x ,
at the' trajectory point is estimated according to:
N
(III.l)
Each of the three remaining trajectory point variables is estimated
by the corresponding analogue to equation (III.l).
-------
III. 3
For example, if the calculation is being conducted for a tra-
jectory point at 0715 hours, then two interpolations are conducted
for each of four variables (x and y wind components, temperature,
and relative humidity). One interpolation is calculated for 0700
and the other for 0800.
This is followed by a temporal interpolation for each set of
spatially interpolated variables given by:
x = x + (x2-x1)-(t-t1) (III.2)
where: x is the interpolated value;
x1 is the spatially interpolated value at hour 1;
x_ is the spatially interpolated value at hour 2;
tj is the time of hour 1;
t is the time of day where t>t .
In the example above, equation (III.2) would cause 25% of
the 0700 value to be added to 75% of the 0800 value.
c. Rates of Change of the Wind Components
The wind component rates of change are easily computed by
dividing the hourly interpolations of these variables by 60.
These values are then passed to the numerical integration routine
where they are utilized in the integration process by the ordinary
differential equations reserved for the trajectory calculation.
d. Computation of Mixing Height
The mixing height is computed by passing a straight line of
slope -.00986 (the dry adiabatic lapse rate) through the trajectory
height-temperature point and determining the abscissa of its inter-
section with the radiosonde curve. The mixing depth is calculated
by subtracting interpolated height-above-sea-level from the mixing
height.
Given the interpolated values (z ,T ) of height above sea level
and temperature at our trajectory point, the module selects the
-------
III.4
appropriate A.M. or P.M. radiosonde curve depending of the value
t . Each of these curves is defined by an integer N or N
o ' 6 a.m. p.m.
and a sequence of points (z^.T,) k=l,2,...N ( .. The mixing
height, h , is computed according to:
T -T, +m. z, -az
, o k K k o
n =
m
where k is the first integer such that l
-------
III. 5
responding atmospheric pressure p , the mixing height, h , the weather
o m
cover, WQ, and the average concentration of NO-, C , in the column.
The module's output is the specific absorption rate of NO , k .
2. Derivations
As a function of the time of day, t (Pacific Standard Time), a
value of the solar zenith angle 6, is interpolated from a table derived
from the equinox curve of Figure 2 in Chapter II of [1]. The air mass,
m, is interpolated as a function of 6 in Table 2 of [1].
Depending of an input cloud-type parameter, a value of average
fraction of radiation transmitted, (3, is interpolated in Table 10 of
o
[1], For each wavelength band of 100A in a wavelength region of
3000A
-------
7.2
III. 7
4000A
, ,273+T
in"16'
10
A=3000A
In passing we note that in our formulation equation (III. 3) reduces
to (11-29) of [1] when a=0 and w =0 (\=w ), but equation (11-38) does not.
3. Units
All lengths are in meters .
Time is in minutes .
Temperature is in °C.
All pressures are in millibars.
-2 -1 7 -1
I . are in photons meters minutes (100A)
0 A.
C is in pphm.
2 -1
k.Trt is in minutes
NO 2
w is dimensionless .
o
P is dimensionless .
4. References
[1] P. A. Leighton, Photochemistry of Air Pollution, Academic
Press, New York, 1961.
[2] Handbook of Physics and Chemistry
C. KINETICS MODULE EQUATIONS
1. Introduction
We first characterize an abstract k step -N^ species constant mass,
constant volume chamber photochemical mechsnism. As a consequence of
this characterization, using standard principles of chemical kinetics,
we can write, in compact notation, the differential equations of species
concentration for the constant volume, constant mass situation. Next
we expand our conceptual framework to the chemical kinetics of the con-
-------
III. 8
tents of a pollutant column whose motion over a time-varying field of
urban emission sources induces variations in average column temperature,
colume mass, etc. This is accomplished by adding a series of appro-
priate first order perturbation terms to the species concentration time
derivatives. It is to be emphasized that the number of reactions and
number of species characterizing our imbedded kinetic mechanism are
limited only by computer storage and computation time considerations.
For this reason we are utilizing the IBM 370/155, a high-speed, large
storage computer, in conjunction with an auxiliary FORTRAN system whose
function is to derive a set of differential equations from a corres-
ponding set of photochemical equations. Thus, we have complete flexi-
bility in reaction number, species number and, more importantly, the
photochemical equations themselves.
The inputs to the module are the outputs of the weather, diffusion,
and ultraviolet irradiance modules.
2. Derivations
An abstract k-reaction, N(j--species, constant-volume, constant-mass
photochemical mechanism is characterized by (i) a binary valued kxNg-
matrix B whose (i,j) element, b , is 1 if, and 0 unless, the jth
th
species appears as a reactant in the i (elementary) step; (ii) an
integer valued kxNg- species conservation matrix, P, whose (i,j) element,
P.., is the number of molecules of species j produced (p..>0) or con-
1J th r
sumed (p. .<0) in the i step; (iii) a K-element vector F of reaction
rate factors k., j=l, 2....K; and (iv) a N^-element vector, C , of
initial concentrations, C,g(0), £=!, 2,...,NCT. From the above defini-
tions and elementary considerations of chemical kinetics, we may formu-
late the following set of differential equations:
-------
III. 9
J-i
«>)
(III.
Ci(0) ' Cio
(i = l,2,...,No)
or, in matrix-vector form:
C(t) = p'R/t)
C(0)
where
C(t) =
C0(t)
c3(t)
-------
III.10
R(t)
n
C£(t)
ct0
m
3. Units
All concentrations are in parts per hundred million except CO, which
is in ppm. Rates of change of concentration are in pphm • min and h
m
-------
III.11
(mixing depth) is in meters.
D. EMISSION MODULE
1. Introduction
The emission module computes the contributions to the pollutant
column of street vehicular traffic, freeway traffic and stationary
area sources. In characterizing the module, the Los Angeles Basin is
subdivided into a grid system of 650 equal 4-square-mile squares corres-
ponding to a 50-by-52-mile square total area. To each of these grid
areas is assigned the total number of vehicular freeway and non-freeway
miles travelled daily; and the rate of emission in kilograms per hour
of reactive hydrocarbons, nitric oxide*, and less reactive hydrocarbons
from area sources. The traffic section of the routine requires
emission factors in grams /vehicle mile for reactive hydrocarbons, nit-
ric oxide, carbon monoxide, and less reactive hydrocarbons.
The inputs to the module are the time(t), the position of the
pollutant column (x,y), the mixing depth (h ), the emission distri-
m
butions described above, emission factors described above, and diurnal
usage curves. The diurnal curves are illustrated in Figures II.5 and
II.6. The grid distributions are shown in Figures II.7 - 11.11.
The module's outputs include the emissions from the five sources
street traffic
rr. vehicular sources
freeway traffic
reactive hydrocarbons (propylene)
nitric oxide
less reactive hydrocarbons
area source emissions
in propylene, nitric oxide, carbon monoxide, and less reactive hydro-
carbons. These results are given in both pphm/min and gram-moles/min.
*0xides of nitrogen taken as nitric oxide
-------
Ill. 12
2. Derivations
Given the air parcel or column location (x,y) (in miles) such that
11 < x < 67 and 1 < y < 53, the projection of the column centroid falls
in grid sub-area i, where is is computed according to:
+ 2S[ttl -!].-„
where [z] = integer part of z.
Once the primary grid square is selected, the program determines
the three next closest grids to the pollutant column; that is, the three
grids whose midpoints are nearest the column. At this point, the
emissions attributable to the column are computed utilizing values
from four grids. The contribution to each of the sources for each of
the pollutants is calculated on the basis of the distance of the grid
midpoint to the pollutant column and is inversely proportional to the
square of this distance. These contributions are then modified by the
appropriate diurnal usage curves.
3. Units
All computations are calculated in metric units with the exception
of distances which remain in miles.
E. NUMERICAL INTEGRATION
DIFSUB is the numerical integration program. It is more commonly
known as the Gear routine named after its developer, Professor C. W.
Gear of the University of Illinois.
This routine has three numerical integration routines imbedded in
it. REM utilizes the option known as the Adams method.
No flowchart or other descriptive material is provided for DIFSUB
since this routine was not programmed by PES staff. For additional
information contact Professor Gear at the Department of Computer Science,
University of Illinois; or see the Communications of the Association
for Computing Machinery, March, 1971.
-------
IV. 1
Section IV
Flowcharts
-------
START
MAIN
I */ DIFSUB
DIFFUN
SHGOUT W-
SMGPLT W-
METEOR
BARIER
OLTRAV
PIF2
to
SOURCE
OVERVIEW OF REM
-------
1 OF 15
IV. 3
START
MAIN
N
J
MODEL
INITIALIZATION
INPUT
COMPLETE
DAILY
DATA
BASE
CONVERT
DATA FOR
COMPUTER RUN
OUTPUT
STARTING
VALUES
FLOW OF MAIN PROGRAM
-------
2 OF 15
IV. 4
-------
3 OF 15
IV. 5
DIFSUB*
1. Called by MAIN
2. Calls DIFFUN
DIFSUB always returns to the calling
program after completing its tasks.
*Flowchart of DIFSUB not provided since
no programming was done on this rou-
tine by PES.
-------
4 OF 15
IV. 6
(• START "N
V DIFFUN J
TRANSFER
ACCUMULATING
SPECIES
SET
NON-
ACCUMULATING
SPECIES
RESET
SPECIES
TO ZERO
COMPUTE
CHAIN
YIELD
SET NEW
TRAJECTORY
POSITION
COMPUTE
RATES OF
REACTION
COMPUTE
RATES OF
CHANGE
5ft
FLOW OF DIFFUN SUBROUTINE
-------
IV. 7
5 OF 15
TRW.
POSITION OFF
MAP
NO
ADD CO
PERTURBATIONS
TO RATES OF
CHANGE
ADD
PERTURBATIONS
TO RATES OF
CHANGE
MODIFY
RATES OF
CHANGE
RESTORE
ACCUMULATING
SPECIES
RETURN
DIFSUB
v
J
-------
IV. 8
6 OF 15
SET TIME AND
CONCENTRATIONS
INTO PLOT
ARRAY
SAVE PEAK
VALUES, TIMES
AND
HALF TIMES
UTP
ONLY AT
NTERVAL
TIME TO
OUTPUT.
YES
\OUTPUT
\CONCENTRATIONS
VlME, LOCATION/
V EMISSION 7
\ RATES /
SET NEGATIVE
CONCENTRATIONS
TO ZERO
RETURN
MAIN
A
J
FLOW OF SM60UT SUBROUTINE
-------
IV.9
7 OF 15
f START "\
V METEOR J
INITIALIZATION
SET AM
RADIOSONDE
CURVE
CALL
BARIER
P. 14
SET BARRIER
INDICATORS
FOR TRAJECTORY
POINT
NEW
RADIOSONDE
C^RVE
NO
1
r
SET PM
RADIOSONDE
CURVE
^
\
t
SET DISTANCES
FROM TRAJECTORY
POINT TO
MET. STATIONS
TRAJ.
POINT ON
TATI
NO
MET. STATIONS
UALI
DETERMINE
INTERPOLATED
MET.VALUES
FOR NEXT HOUR
DETERMINE
INTERPOLATED
MET. VALUES FOR
CURRENT HOUR
FLOW OF METEOR SUBROUTINE
-------
IV. 10
8 OF 15
TIME
INTERPOLATE
METEOROLOGICAL
VALUES
SET
INTERPOLATED
GROUND ELEVATION
VALUE
SET TIME
INTERPOLATED
METEOROLOGICAL
VALUES
SET
INTERPOLATED
GROUND ELEVATION
VALUE
CALCULATE
PRELIMINARY
MIXING HEIGHT
SET MIXING
DEPTH FROM
MIX. HT. AND
GROUND ELEV.
BELO
MINIMUM MIX.
DEP
RESET
MIXING HT.
AND
MIXING DEPTH
SET WIND
COMPONENTS
XDOT AND YDOT
RETURN
DIFFUN
C8(0
-------
IV. 11
SET UP ARRAYS
FOR SUN ANGLE,
AIR MASS, CLOUD
TYPE AND WATER
VAPOR PRESSURE
CALL
PIF2
SET SUN
ANGLE OF
INCIDENCE
CALL
PIF2
P. 13
SET AIR MASS
AT THIS ANGLE
9 OF 15
CALL
PIF2
\ PHa
SET BETA,
AVERAGE FRACTION
OF INCIDENT
RADIATION
TRANSMITTED
/ CALL \
( PIF2 )
13 /
SET PARTIAL
PRESSURE OF
WATER VAPOR IN
THE ATMOSPHERE
LAMBDA=
3000
ANGSTROMS
10A
FLOW OF ULTRAV SUBROUTINE
-------
OF 15
IV. 12
DETERMINE
DIRECT AHD
SKY RATES OF
IRRADIANCE
\
1
LAMDA =>
LAMDA +
1000A
CALCULATE
SPECIFIC
ABSORPTION
RATE OF UV
IRRADIANCE
(RETURN SMGFCT)
-------
IV. 13
11 OF 15
SET TIME
TO DAYLIGHT
TIME
EVALUATE
TRAFFIC
DISTRIBUTION
CURVES
FIND GRID
IN WHICH
TRAJECTORY
POINT LIES
FIND THREE
CLOSEST
ADJACENT
GRIDS
SET
MIDPOINTS OF
THE FOUR
GRIDS
FIND DISTANCES
OF TRAJECTORY
POINT TO
GRID MIDPOINTS
TRAJ.
POINT ON
MI DPT
NO
CALCULATE
CONTRIBUTIONS
OF TRAFFIC GRIDS
(1/DIST2)
FLOW OF SOURCE SUBROUTINE
-------
IV. 14
12 OF 15
ADD TRAFFIC
CONTRIBUTIONS
OF TRAJECTORY
GRID ONLY
CALCULATE
CONTRIBUTIONS OF
SOURCE GRIDS
(1/DIST2)
ADD AREA SOURCE
CONTRIBUTIONS
OF TRAJECTORY
GRID ONLY
YLIGH
SAVINGS
TIME
YES
RESET TIME
TO STANDARD
TIME
-------
C START PIF2 J
YES
VAR.
PAST END OF
VARIABLE
,LIST
NO
IND. YES
VAR. BEFORE
START OF IND.
VARIABLE
DETERMINE
INTERVAL OF IND,
VARIABLE IN
INDEPENDENT
VAR. LIST
NO
THIRD
POINT ON
HIGH
SIDE
THIRD
POINT ON
LOW
SIDE
SELECT THIRD
POINT CLOSEST
TO THE
INDEPENDENT
VARIABLE
IV. 15
INTERPOLATE
INTO
DEPEN. VAR.!
LIST. USING TWO
HIGHEST POINTS.
INTERPOLATE
INTO DEPEN.
VAR. LIST
USING TWO
LOWEST POINTS
INTERPOLATE
INTO DEPEN.
VAR.LIST USING
TWO'INTERVAL
PTS. AND ONE
ON LOW SIDE
INTERPOLATE
INTO DEPEN.
VAR.LIST USING
TWO INTERVAL
PTS. AND ONE
ON HIGH SIDE
INTERPOLATE
INTO DEP. VAR.
LIST USING TWO
INTERVAL
POINTS AND
THIRD POINT
SELECTED
13 OF 15
(RETURN ULTRA?)
-(RETURN ULTRA?)
FLOW OF PIF2 SUBROUTINE.
-------
START BAR IERJ
DETERMINE
SLOPE AND
INTERCEPT FROM
TRAJ. PT. TO EACH,
METEOR. .STATION
INITIALIZE
STATION,
BARRIER AND
SEGMENT INDICES
IV. 16
DETERMINE
SLOPE AND
INTERCEPT OF EACH
BARRIER SEGMENT
ALL
BARRIERS
CHECKED
14 Of* 15
STEP BARRIER
INDEX
INITIALIZE
SEGMENT INDEX.
INCLUDE
METEOROLOGY
STATION
ALL
STATIONS
CHECKED
STEP STATION
INDEX
INITIALIZE
SEGMENT, BARRIER
INDICES
DOES\ YES
BARRIER
BLOCK
STATION
STEP
BARRIER
SEGMENT
INDEX
ALL
SEGMENTS
CHECKED
EXCLUDE
STATION
FLOW OF BARRIER SUBROUTINE.
-------
IV. 17
15 OF 15
START SMGPLT
INITIALIZATION
OUTPUT TIME
CONCENTRATION
DATA
SET
PLOTTING
LIMITS
SET
PLOT SCALE
OUTPUT
PLOT
HEADINGS
DETERMINE
OUTPUT
CHARACTERS
^OUTPUT ONE LINE
OF CHARACTERS,
AND VERTICAL,
AXES
OUTPUT
HORIZONTAL
AXIS
ALL PLOTS
COMPLETED
(RETURN MAIN
FLOW OF SMGPLT SUBROUTINE.
-------
V.I
Section V
Model Inputs
Model inputs must be prepared exactly as described. A complete sample
input deck is shown in Section VI.
1. Card Type 1
Heading data 1 array 4 cards.
a. Alphanumeric data used as heading information. May be blank
cards.
FORMAT (20A4)
0 8
1 0
2. Card Type 2
Starting data 4 items 1 card.
XPOS - X coordinate of starting position in miles. Range: 17.0
to 67.0. Cols 1-10.
YPOS - Y coordinate of starting position in miles. Range: 1.0
to 53.0. Cols 11-20.
GOTIME - Simulation start time (local standard time) in hours and
decimal fraction hours. Range: 7.000 to 18.999. Cols 21-
30.
XMIN - Minimum mixing depth in meters. Current value used is
generally 50 meters. Cols. 31-40.
-------
FORMAT (4F10.3)
V.2
0 11 22 33 44 i
1 01 01 01 01 C
XPOS
YPOS
GOTIME
XMIN
3. Card Type 3
Starting data 6 items 1 card.
PRMT (2) - Length of simulation in minutes and decimal fraction of
minutes. Cols. 1-10.
PRMT (3) - Minimum step size for numerical integration in decimal
fraction of minutes. Normal range is from .000001 to .001.
Cols. 11-20.
PRMT (4) - Maximum step size for numerical integration in minutes
and decimal fraction of minutes. Normal range is from 0.1 to
20.0 minutes. Cols. 21-30.
PRMT (5) - Initial value of step size for numerical integration in
minutes and decimal fraction of minutes. Normal range is from
.0001 to .01 minutes. Cols. 31-40.
PRMT (6) - Error criteria for numerical integration in decimal
fraction of minutes. Normal range is from .0001 to .10 minutes.
Cols. 41-50.
PRMT (7) - Indicator for plot routine. Determines how often data
is to be saved for output in the plot. Reasonable value is
usually .25*PRMT (2) in minutes. For example, if PRMT (2) is
300.0 minutes, then PRMT (7) would be 75.0. Cols. 51-60.
-------
V.3
FORMAT (6F10.3)
0 11 22 33 44 55 66 {
1 01 01 01 01 01 01 C
PRMT(2)
PRMT(3)
PRMT(4)
PRMT(5)
PRMT(6)
PRMT(7)
4. Card Type 4
Starting data 14 items 1 card (all values must be right justi-
fied.
NDIM - Number of accumulating species. For CO simulation only,
value must be 1. Cols. 1-5.
NC - Total number of species. For CO simulation only, value must
be 3. Cols. 6-10.
NEQ - Number of elementary steps in mechanism. For CO simulation
only, value must be 1. Cols. 11-15.
NLL - Number of lines to be used in plots, 0 gives 50 lines. Cols 16-20
NSTA - Number of meteorological stations. Cols. 21-25.
MF - Indicator for numerical integration method. Range: 0 to 2;
values are: 0 - Adam's method
1 - Gear (matrix inversion) for less than 5
equations
2 - Gear (Gaussian elimination and back substitu-
tion) for 5 or more equations. Cols. 26-30.
MAXDER - Maximum number of derivatives calculated. Use MAXDER=12
if MF=0, or use MAXDER=6 if MF=1 or 2. Cols. 31-35.
ICLOUD - Index which defines seven different cloud types for ultra-
violet radiation module. Range: 2 to 8, values are:
Leighton, Phillip A., Photochemistry of Air Pollution, Academic Press,
New York, 1961, Page 40, Table 10.
-------
V.4
2 - cirrus
3 - Cirrostratus
4 - Altocumulus
5 - Altrostratus
6 - Stratocumulus
7 - Stratus
8 - Fog
Cols. 36-40.
NRAD - Number of radiosonde points (maximum = 15) in the curve
describing the mixing height. Cols. 41-45.
ITRAF - Traffic module operation indicator. 0 - module does not
operate, 1 - module operates. Cols. 46-50.
ISOUR - Area source module operation indicator. 0 - module does not
operate, 1 - module operates. Cols. 51-55.
IHT - Indicates if mixing height is entered by way of radiosonde
data, or if it is input directly in meters. 0 - radiosonde
data, 1 - input mixing height. Cols. 56-60.
IN - Indicates if model is to operate again with data from the
previous run or a new set of data is necessary. 0 - a complete
set of data is necessary for the run, 1 - return with the
previous data. Must input only card types 1, 2, 3, 4, 5, 6 and 7.
Need card type 6 only if IHT = 1. Cols. 61-65.
IPST - Time indicator. 0 - traffic module operates on daylight time,
1 - traffic module operates on standard time, Cols. 66-70.
-------
V.5
FORMAT (1415)
!0 00 11 11 22 22 33 33 44 44 55 55 66 66 77 £
1 56 01 56 01 56 01 56 01 56 01 56 01 56 01 C
\
"c
\
\
\
*F
\
\
\
\
\
x
Hr
\
5. Card Type 5
Output data 1 item 1 array 1 card.
INTR - Integer. Positive value causes outputs at approximate
times specified in array STEP. Zero or negative value
causes outputs each time through. Cols. 1-2.
STEP - Array of up to 15 times for output in minuts. Output
usually occurs within 1/2 minute of time specified. Cols.
6-80.
FORMAT (12, 3X, 15F5.2)
.000 00 11 11 22 22 33 33 44 44 55 55 66 66 77 77 8
123 56 01 56 01 56 01 56 01 56 0J. 56 01 56 01 56 Q
I
N
T
R
STEP
(1)
STEP
(2)
STEP
(3)
STEP
(4)
STEP
(5)
STEP
(6)
STEP
(7)
STEP
(8)
STEP
(9)
STEP
(10)
STEP
(ID
STEP
(12)
STEP
(13)
STEP
(14)
STEP
(15)
6. Card Type 6
Mixing height data 1 array 1 card
XHITE - Array of up to 12 values of mixing height in meters.
Cols. 1-72, optional (used only if IHT is set to 1; see IHT,
card type 4). Represent hours 0700 to 1800, respectively.
-------
FORMAT (12F6.0)
V.6
0 00 11 11 22 33 33 44 44 55 66 66 77 8
1 67 23 89 45 01 67 23 89 45 01 67 23 „
XHITE
(1)
XHITE
(2)
XHITE
(3)
XHITE
(4)
XHITE
(5)
r, (
XHITE
,(6)
XHITE
XHITE
(8) ,
XHITE
(9).
ICHITE
(10)
ICHITE
(11)
ICHITE
(12)
7. Card Type 7
Species data—2 arrays—j^—j + 1 cards, 4 species per card
**
LABEL - Alphanumeric array of NC species names. Cols. 1-8, 21-28,
41-48, 61-68.
Y - NC double precision values of initial species concentrations
in pphm, except CO is in ppm. Cols. 9-20, 29-40, 49-60, 69-80.
0 00 22 22 44 44 66 66 8
1 8Q 01 89 01 89 01 89 P
LABEL
Y
LABEL
Y
LABEL
Y
LABEL
Y
8. Card Type 8
Sun Angle Data 1 array 2 cards
ZENITH - Twelve pairs of an hour in standard time coupled with the
sun angle at that time. Hour in cols. 1-5, 11-15, 21-25,
31-35, 41-45, 51-55, 61-65, 71-75. Sun angle in cols. 6-10,
16-20, 26-30, 36-40, 46-50, 56-60, 66-70, 76-80. Pairs carry
*[ ] means largest integer.
**See card type 4.
C_H, is set to .015 x CO in pphm
j o
DUMHC is set to .15 x CO in pphm
-------
V.7
over onto second card.
0 00 11 11 22 22 33 33 44 44 55 55 66 66 77 77 8
1 56 01 56 01 56 01 56 01 56 01 56 01 56 01 56 0
Z
(1,1)
Hour
Z
(1,2)
kngle
Z
(2,1)
lour
Z
(2,2)
ingle
Z
(3,1)
lour
Z
(3,2)
Angle
Z
(4,1)
Hour
Z
(4,2)
ingle
Z
(5,1)
lour
Z
(5,2)
ingle
Z
(6,1)
tour
Z
(6,2)
ingle
Z
(7,1)
Hour
Z
(7,2)
Angle
Z
(8,1)
Hour
Z
(8,2)
Angle
Z - array ZENITH
0 00 11 11 22 22 33 33 44 8
.1 56 01 56 01 56 01 56 01 o
Z
(9,1)
Hour
Z
(9,2)
ingle
Z
0-0, D
Hour
Z
(M, 2)
kngle
Z
(H,l)
Hour
Z
01,2)
Angle
Z
O2, W
lour
Z
(12,2)
Angle
Z = array ZENITH
9. Card Type 9
Meteorological data—2 arrays—maximum of 32 stations.
AIR - 3 cards per station, 4 hours per card, 12 hours per station.
Cards are as follows:
IDNUM - column 1-3, station ID code
AIR - column 4-5, start hour of data, values are 7, 11 or 15
6, number of readings on card, max = 4
7-23
24-40
41-57
58-74
6 items in 17 columns in the following order:
3 columns - wind direction in degrees
3 columns - wind speed in miles/hour, last column is
in tenths of a mile/hour
3 columns - temperature in °F
2 columns - dew point in °F
6 columns - blank
hourly data—17 columns
-------
V.8
FORMAT (A3, 3X, 4(F3.0, F3.1, F3.0, F2.0, 6X), 6X)
Conversion; Wind speed and direction are converted to XDOT and
YDOT, the X and Y components of the wind velocity at the column.
Temperature is converted to °C. Relative humidity is determined
from temperature and dew point. The station order on card types
8 and 9 must be exactly coincidental.
0 00000 01 11 nil
1 34567 90 23 5678
22 22 23 3333 .44 44 44 4555
34 67 90 2345
01 34 67 9012
55566 66 6666
78901 34 6789
I
D
H
o
U
R
W
D
I
R
W
S
P
D
T
E
M
P
H
U
M
I
D
W
D
I
R
W
S
P
D
T
E
M
P
H
U
M
I
D
W
D
I
R
W
S
P
D
T
E
M
P
H
U
M
I
D
W
D
I
R
W
S
P
D
T
E
M
P
H
U
M
I
D
77
45
10. Card Type 10
Meteorological Stations—2 arrays—maximum of .8 cards, 4 stations
per card
WPOS - Array of triplets, 4 stations per card, maximum of 32
stations. Cards are as follows:
IDMET - column 1-3, alphanumeric station ID code.
WPOS - column 4-8, X coordinate of station position in miles, last
2 columns in hundredths of miles. Column 9-13, Y coordinate of
station position in miles, last 2 columns in hundredths of
miles. Column 14-20, ground elevation of station in feet, last
2 columns in hundredths of feet. This 20 column grouping is
repeated in columns 21-40, 41-60, 61-80 as required, and on
additional cards if necessary.
-------
V.9
FORMAT (4(A3, 2F5.2, F7.2))
The station order on card types 9 and 10 must be exactly coincidental.
0 0
1 3
T
D
0 0
4 8
X
0 1
9 3
Y
1 2
4 0
ELEV
2 2
1 3
T
LI
2 2
4 8
X
2 3
9 3
Y
3 4
4 ' 0
ELEV
4 4
1 3
I
L)
4 4
4 8
X
4 5
9 3
Y
5 6
4 0
ELEV
6 6
1 3
1
LI
6 6
4 8
X
6 7
9 3
Y
7 8
* 0
ELEV
11. Card Type 11
Traffic data 1 array 1 card
EMIT - Array of emission factors in grams per vehicle mile for
chemical compounds in the following order:
1 - Propylene (C_H,). Cols. 1-5
J o
2 - Nitric oxide (NO). Cols. 6-10
3 - Carbon monoxide (CO). Cols. 11-15
4 - Dummy hydrocarbons (DUMHC). Cols. 16-20
FORMAT (10F5.2)
0 00 11 11 22 8
1 56 01 56 01 0
EMIT
(D
C3H6
JMIT
(2)
NO
EMIT
(3)
CO
EMIT
(4)
DIMHC
-------
V.10
12. Card Type 12
Omit if ITRAF is not set to 1 (see card type 4)
Vehicular mileage by 4-square mile grid 2 arrays 52 cards
FGRID - 650 values representing the freeway traffic mileage per
grid in thousands of miles. 25 values at 3 columns each per
card.
SGRID - 650 values representing the street traffic mileage per
grid in thousands of vehicular miles. 25 values at 3 columns
each per card.
FORMAT (25F3.0)
0 00 00 01 11 11 11 22 22 22 33 33 33 34 44 44 44 55 55 55 66 66 66 67 77 77 £
1 34 67 90 23 56 89 12 45 73 01 34 57 90 23 5,6 89 12 45 78 01 3,4 67 90 23 56 C
G
R
I
D
6
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
C
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
G
R
I
D
C
R
I
D
G
R
I
D
G
R
I
D
FGRID and SGRID have the same formats.
13. Card Type 13
Omit if ISOUR is not set to 1 (see card type 4)
Area source emissions by 4-square mile grid 3 arrays 78 cards
AGRID1 - 650 values representing distribution of reactive hydro-
carbon (propylene) emissions by grid in units of kilograms per
hour. 25 values at 3 columns each per card.
AGRID2 - 650 values representing distribution of nitric oxide (NO)
emissions by grid in units of kilograms per hour. 25 values
at 3 columns each per card.
AGRID4 - 650 values representing distribution of less reactive
-------
V.ll
hydrocarbons (DUMHC) emissions by grid in units of kilograms
per hour. 25 values at 3 columns each per card.
FORMAT (25F3.0)
0 00 00 01 11 11 11 22 22 22 33 33 33 34 44 44 44 55 55 55 66 66 66 67 77 77 8
1 34 67 90 23 56 89 12 45 78 01 34 67 90 23 56 89 12 45 78 01 34 67 90 23 56 0
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
A
G
R
D
AGRID1, AGRID2, AGRID4 have the same formats.
14. Card Type 14
Chemical mechanism data 3 arrays NEQ*cards (see card type 3)
RCT - Alphanumeric array of NEQ chemical reactions. Cols. 1-24.
IR - 5 integer values per card of species designations. These
values reflect the mechanism being used. Columns 25-39,
3 columns per number, right justified.
C - array of NEQ double precision reaction rate constants.
Cols. 40-54.
FORMAT (3A8, 513, D15.8, 26X)
22 22 33 33 33 34 55 8
45 78 01 34 67 90 45 0
RCT
I
R
0)
I
R
(2)
I
R
0)
I
R
(4)
I
R
(5)
C
*See card type 4.
-------
V.12
15. Card Type 15
Chemistry parameters 1 array 1 card
Z - Chain termination rate factors
Z(l) Combination of free radicats - current value used is 100.0
Z(2) Combination of radicals with NO - current value used is
0.05
FORMAT (2D20.8)
0 22
1 01
Z(l)
44 1
01 C
2(2)
16. Card Type 16
Omit if IHT is 1 (see card type 4)
Radiosonde data 2 arrays maximum of 8 cards
ARADIO - AM radiosonde curve of NRAD (see card type 4) elevations
in meters followed by NRAD temperatures in degrees centigrade.
The first temperature is the temperature at the lowest
elevation, the second temperature is the one at the next
elevation etc. The temperatures must begin on a separate card.
If NRAD is greater than 8, two cards must be used for elevation,
followed by two for temperature. Ten columns are used for
each elevation and temperature. The last column in each field
is in tenths of a meter or degree.
PRADIO - PM radiosonde curve of NRAD elevations and temperatures.
Data is identical to that for the array ARADIO. If only one
radiosonde curve is available, the same one is used for both
-------
V.13
curves. If two curves are available, the model assumes they
have the same number of points.
FORMAT (8F10.1)
0 11 22 33 44 55 66 77 e
1 01 01 01 01 01 01 01 o
ELBV
ELBV
ELEV
ELBV
ELXV
ELEV
ELEV
ELEV
0 11 22 33 44 55 66 77 . 8
1 01 01 01 01 01 01 01 0
TEMP
TEMP
TEMP
TEMP
TEMP
IHff
TEMP
TCMP
APRADIO
PRADIO
-------
VI. 1
Section VI
Sample Input Deck for REM
Card Type 1
TRAJECTORY ftUN ON SEPTEMBER 26, 1972
SEPTEMBER 30, 1969
390 MINUTE HUN
25 METEOROLOGICAL STATIONS
RUN TO POMONA AT 13:30
Card Type 2
57.5261 10.9335
7.01
50.0
Card Type 3
390.0 .0001
0. 1
.001
.01
100.0
Card Type 4
12 25 33
0 25
12
2-
1
0
Card Type 5
IS 0. 0. 30. 60. 100. 140. 180. 220. 260. ISOO. 330. 360. 380. 386. 339.
Card Type 6*
260
Card Type 7
220
2G0T 300 , 315. 450 500
400
N02 0.106302D 02NJ 0.3649160 0203 0.16H827D 01C3H6 0.25 0 02
HN03 0.0 CH3CHO 0.0 ' HCHO 0.0 CO < PPM )0.1666381,' 02
D 03
CH30N02
CH3
C2H302
C2H303
0.0
0.0
0.0
0.0
0.0
C2HJN05
H '
CH30
C2H402
0.0
0.0
0.0
0.0
FIXED 020.0
OH
CH302
0
0.0
0.0
0.0
D'JMHC
HO 2
C2H30
NO 3
0.25
0.0
0.0
0.0
Card Type 8
87
7
15
57
9
16
63
68
9
17
52
80
10
U
36
12 34
13
38
*Card type is shown only for illustrative purposes. As a result of inputs
shown on Card Type 3, this data should be omitted.
-------
VI. 2
Card Type 9
LAX07418 5
LAX11422 8
LAX15423 13
BUR07428 3
eU°11415 7
BUS15414 9
IGB074 90 1
LGB11417 8
LGB15428 12
NTB07400 00
NT81 1400 00
NTB15420. 13
CAP07418 3
CAPU422 7
CAP1S4.?7 6
HOL074225 2
HOL 114225 2
HOH54225 2
PIC07427 1
PIC. 114 IB 2
PIC15422 1
RLA074 4'3 I
R.LAU4180 1
RL A 15422 5 1
VEK07413Q 2
VER114180 9
VfcfU5422514
HNT074 7 3
HMT11420 7
WNT15420 8
WLA0742S 1
WLA 11420 2
WLA15422 5
WHT074KiO 2
WHT1H203 3
WHT154158 5
PAS07434 1
PAS11416 3
PAS15420 5
VEN'07413 1
VEN11420 5
VSN15420 10
RVA074 ~f3 2
PVA114156 5
P. VA 15422 5 9
RES07427 1
RHS11422 3
P.ES15427 8
SB 0741S 5
RB 114?2 5
RB 15420 7
CPK074315 1
CPK114725 ?
CPK154315 4
KF1074180 1
KF1114225 5
KFI154203 8
ENC07422 3
ENC114.16 4
66621114
76641111
76601084
7052
9150
9245
65571122
805911 19
R7531038
64621115
0000
78621083
18 6
24 9
25 10
12 3
15 9
9 9
11 3
13 12
31 8
13 2
22 9
20 7
13 3
22 8
22 6
225 'i
225 2
225 6
22 2
18 2
18 1
180 1
225 1
225 1
130 3
18010
22.5 9
11 3
22 9
19 6
• 1>J 1
:?o 3
22 3
158 2
203 6
ISC 4
13 1
16 3
20 3
16 2
25 5
22 5
90 3
15n 6
225 6
20 1
27 3
29 6
20 4
22 5
18 7
248 1
315 3
338 4
135 2
225 7
225 6
24 2
12 8
69641115
78631103
75581074
7753
915iJ
944>i
70601125
82571108
85521081
72621116
81601104
73611075
24 8
25 12
24 12
9 9
ID 9
10 9
17 7
17 14
10 7
19 6
20 12
00 00
18 6
22 13
22 6
13:3 2
225 1
225 6
18 2
22 1
18 1
45 1
180 1
2? 5 1
135 2
IbO 8
225 7
27 3
20 9
16 7
20 2
22 4
22 3
l.'JO 3
203 6
158 5
9 3
20 4
16 2
18 2
22 5
22 3
90 5
203 8
203 5
20 1
25 6
31 2
20 2
22 7
20 9
22b 1
293 6
338 3
203 3
225 9
225 7
34 2
13 8
71651119
79621095
73611074
8251
9455
8651
74591125
82581098
80611078
762011.17
82591095
0000
113 3
25 13
21 8
90 7
15 14
10 9
17 6
17 10
10 6
20 5
20 13
17 6
18 5
27 11
22 6
225 2
180 2
225 2
18 2
18 1
22 1
180 1
225 2
225 1
135 7
18010
225 7
24 5
20 3
16 5
20 1
22 5
22 1
1HO 3
180 7
15 S 3
11 2
22 7
16 1
18 5
25 6
13 1
108 ti
20310
203 5
22 3
25 8
11 1
20 4
20 7
20 6
248 1
293 6
158 1
225 2
225 8
225 5
2 3
27 8
76611089
70631070
9345
8153
68561122
84591091
73611081
30621113
80531087
70621073
093019
093020
093021
093025
093026
093027
093028
093029
093030
093037
093038
093039
093046
093047
093048
093055
093056
093057
093064
093C65
093066
-------
VI. 3
Card Type 9 (continued)
ENC15427
ELM074315
ELMU4225
ELM154180
COM07416
CUM11422
CUM15429
POM074135
POM114270
POM154293
AZU07413
AZU11422
AZ.U15420
ONT07436
UNT11423
ONT15426
9
1
4
8
2
3
14
2
5
5
1
6
13
3
9
16
6949
8553
9161
27 6
158 1
203 6
13012
29 2
27 7
29 10
ISO 2
248 6
270 5
13
22
20
17
28 12
26 16
7753
9354
8955
29 5
180 1
18010
180 6
22 3
29 9
29 7
203 2
270 6
270 4
18 3
20 9
20 7
20 6
27 14
00 00
7952
9344
0000
6 4
180 3
22511
180 8
18 4
29 12
31 8
270 2
248 7
293 3
18 3
20 12
20 3
27 6
00 00
27 9
093069
093070
093071
093082
093083
0930B4
8659
0000
8156
Card Type 10
LAX 32
CAP 40
VER 41
PAS 4V
RB 31
ELM 51
ONT 72
27
33
31
42
21
35
35
99
250
180
870
10
300
916
BUR 33
HOL 35
WNT 63
VEN 26
MSH 28
COM 41
45
37
30
30
48
23
699 LGB 45
300 PIC 33
950 WLA 29
5 RVA 48
850 KFI 53
100 PCM 69
18
34
3'V
29
23
34
34 NTS 51
100 RLA 45
200 WHT 51
180 RES 22
100 ENC 24
350 AZU 59
16
26
28
46
43
41
27
100
250
800
800
600
Card Type 11
2 4 80 16
Card Type 12
11 22
143 45107 5293'.
13134354534032.3195
625
391326217
351 32 257
34h) 305
404 366
295160 391
278
236
241 56383102243253
1-59131 3142(19104
302 527 ?.C
281129444 ? 3-
283332 265
18
80 6
97 95 61
115160163148
153 7B 67169
131112183150 257131
;>213 144 25398156161
145 68 137
59 83330251138108207
15
59
91
77
85
16
12
70
59 42 21
15
16
101295511 51 531449623576 205263156
473141333426466795452585513224 79234285147
242137 205436620257257297306311254362291214203170173159
356 215194290335193109
352 140275 423 13 90 28 41 74 66 65 43 11
169271311383390456426494178395190 51
222 126 17104150 37
242 156158109
155 68144
59 94119
F01
F02
F03
F04
FO'j
F06
F07
FOB
f-09
F10
Fll
F12
F13
F14
F15
F16
F17
F18
F19
F20
F21
F22
F23
F2-V
F25
-------
VI. 4
Card Type 12 (continued)
166
5119 18
7218196 52
3+104277 94114 10
3612'j!.07248215100 29
241331.35116 74 80 11
33 48114195116166203190 98
75157126153175 6012521 2.:! 002 1 92492661 82
15
26
19
32
751162C9:!202b5149168 163 28
55 36
4 5 25176 83
30 49 28142145
4612L188266254157413306317162
1239304249175 '94283 1^2282 1 65 155 1 73 2 59 i!2 8304 1691 23 109 11
9360324231124 79204268430204 T.J81 8630'5 IJ4232 3 24111 3 31 6
124215412322210121270312301 93 96104270:. 79290197 82 48 30
lcl925532335!U76265293257298241230135H4 62129114 95 57 4
2226129b3972B53142I>429;>25539626322923413i;!40156l03 21 9 7
106266413392331393273322105265253190295119 42 22 1 37 1 5
88301-2^525531 1'.734773642922732991 981 05 69 21 62 42 90 33 12
1254762543382614124l351232?324293264195127 43 39:.'27121 45 53 27 5
47 64 9911723151234<>45947l}468642579290199186 75 1081 73 1 3 1 1 80142 18 43 50 73
18 7 98283403549512572473226217229285243250242114190280244 71 10 99
9 25 2 30 70 322604783752591641161902.36227236185160130155232143 66 77
25 30 17 12 59 6J121155H6 25?45242191162318308241218183 77126149151 76 32
139118107174182 }f-02b82B02471.';l0204204 21173241212156 91 5 7 1 6 6 2
239250249250264414323276212215 4131 58 91113 21 1 2 12
137255238238266295251233 54 20 27157 50 4 1 '2
60155202223252221152 85 66 94 73 7 2 2 2
35 12 3414914t17213? 43 28 17 1 421221 2
60125 50
F26
SOI
S02
S03
S04
SOS
S06
SO 7
S08
S09
S10
Oil
S12
S13
S14
S15
S16
S17
S18
S19
S?0
521
S22
S23
S24
S25
S26
Card Type 13
30
30
30
30
30
1
31
31
31
30
1
31
31
31
20
35
t>
3U
35
35
3
1
25
40
39
24
5
15
30
39
39
36
32
10
30
40
40
39
3.9
15
15
30
39
39
36
32
10
25
15
1
6
6
94
56
51
45
45
33
33
23
19
34
53
53
42
31
15
30
30
31
31
31
56
56
51
45
45
33
33
33
34
34
53
53
42
31
21
32
32
39
72
46
56
56
51
47
47
35
35
36
38
38
48
48
40
1
5
31
33
33
47
61
61
56
56
51
'•9
49
37
37
39
42
42
33
13
37
11
18 18 18
33 67 33
50 78 34
63 45 36
61 36 36
56 70 70
56 70 70
51 90' 90
491 101 10
49110110
37 oO 60
37 60 60
39 48 48
42 36 36
42 36 6
23 31 31
13 31 16
27 10
18
35
34
36
39
57
57
70
82
82
53
53
44
35
35
30
5
3
33
34
36
36
50
45
50
56
56
45
45
40
35
35
25
5
14
33
34
36
36
45
45
50
56
56
45
45
40
35
35
15
5
11
11
11
32
34
34
46
42
35
34
34
62
62
49
37
37
5
11
11
11
11
13
15
34
36
36
35
19
4
62
62
43
37
37
5
11
11
11
11
11
11
13
13
33
33
20
12
22
49
49
'-3
39
39
11
11
11
11
11
11
11
11
11
11
11
4
34
39
36
36
38
41
41
3
10
11
11
11
11
11
11
11
11
11
11
1 1
4
9
39
36
36
3d
41
41
1
11
5
11
11
11
11
11
11
11
11
11
11
5
2
32
1
31
31
30
30
3
11
11
11
11
11
11
11
6
2
32
1
31
31
30
15
11
11
10
5
2
1
1
4 30
4 30
32 30
20 10
HOI
H02
H03
H04
H05
H06
H07
H08
H09
H10
Hll
H12
H13
H14
H15
H16
H17
H18
H19
H20
H21
H22
H23
H24
H25
H26
NO I
NO 2
-------
VI. 5
Card Type 13 (continued)
8 8
8 8
8 9
8 9
1
4
8 9
1
1
12 14
13 18
9 18
1 1
10
9
19
6
3
4
12
15
17
10
9
4
14
10
5
6
584
6 29
13 25
19 21
18
19
10
7
4
8
15
15
10
9
6
20
19
10
9
4
a
17
15
12
9
11
60
40
-\
24
54
5
5 11
6 15 20 2 5 11 11
18 21 23 25 37 37 30 22 13 11 11
18 21 23 23480 35 55 21 21 11 11
19 24 29 63603 30 26 22 22 13 11
21371 34 61141 25 24 23 23 1 9 11
21 27 34 34 25 25 24 23 23 19 19
29 25 22 22 24 24 24 92 2316? 27
29 25 22 22 24 2* 24 23 23106 27
25 22 20 20 31 31 26 22 22 22 18
21 20 19 20 30 3P 29 21 21 19 5
20 20 19 19 37 37 29 21 21
20 20 21 21 21 21 20 19 19
10 11 12 12 12 12 11 10 10
13 19 19 19 19 19 18 18 18
8 9 10 10 9 9 18 16 8
8 9 10 10 9 1 9 10 10
15 12 7488383
15 12 1 1 3 4 1
12 11 9 7 2
9 1 1
10
26 72 95 26 2.0
79 9611 1 6391 91 91131 55
79 96111111491156145 94 94
78130182208576144122101101
96433251271218131120108108
96173251251131131152103108
33 2332141931932452 45 19720.? 150
9 1
12 22
13 13
19 19
16 18
16 19
20
20 44
44 44
44 44
44 44
66 51
94 59
94 94
5 11 11
11 11 11
11 11 11
11 11 11
11 11 11
11 11 11
11 11 11
11 11 11
11 11 11
23 11 11
23 11 11
7
392
6 15 19
21 20 20
12 11 11
20 20 20
18 14 13
16 14 13
40
20 44 44
44 44 44
44 44 44
44 44 44
44 44 44
44 44 44
44 44 44
51 44 44
51 44 44
241137109 44 44
3023323321419.U93245245197150150196137109 44 44
79 79
1
2
79 81
80 84
80 84
3
411
79 991
20
1 11
2 2
81 81
841051
341051
10
71
20
901
891
1201451
1211451
13118
79113
2-1
41
80
41
41
80
251251
25
1251
3 81C81081
6
91
91
90
90
74
57
96
301
801
89179175175333333277170170118113
451501561564264263071891U9
451501561564 26 ^2630 71 89 189
90 98106106204204174143143
98 59
98 19
68 19 19
99108 38
69117117
<':l 2212158104104
90 93 106 1062 042 04 174 14 3 143 21 22121
9310511dll8153153136l22122
9611312-1129103103102101 101
96113129129103 24102101101
80153116 76 81 01 80 70 30
H01SS 56 56 81 5?. 11
58104104
159159137117117
1081061
.'. 081081
18128128
18129129
11
11
11
11
11
11
11
11
11
5
1
18
1
9
17
9
9
15
44
44
44
44
44
44
44
44
'44
45
23
1
7
86
6
85
82
RO
80
11
11
11
11
11
11
8
1
18
I
9
17
9
9
44
44
44
44
44
44
31
1
1
7
86
6
85
82
80
41
11
11
8
5
2
9 8
1 10
1 10
9 9
7
9 8
44
44
40
20
5
1
4
4
21115
21115
90 97
62 21
2 1
?3133101107 83 22 2 1
85
85 3
N03
N04
NO 5
N06
N07
N08
NO 9
NlO
Nil
N12
N13
N14
N15
N16
N17
N18
N19
N20
N21
N22
N23
N24
N25
N26
L01
L02
L03
L04
L05
L06
L07
LOS
L09
L10
Lll
L12
L13
L14
L15
L16
L17
LIB
L19
L20
L21
L22
L23
L24
L25
L26
-------
Card Type 14
VI. 6
NU2+KW
0»02+M
Cm NO
N02+03
C2H303+U2
N03+NO
C3H6*OH =
NO+H02
H+02+M
.C2H30 + 02
N02+OH
OH»03
OH + CO
CH3«-02*M
CH302»Nl>
CH3CH02 + NO
C2H303+NO
C2H302+NO
C2H402+NO
CH30+02
C3H6+0
C3H6+03
C3H6+CH-02
C2H402+02
C3H6«-H02
C3H6»ME02=
C2H30+M
CH30+N02
= MOO
^03 + M
=N02*02
=N03*U2
= C2H3U2+Q3
= 2N02
CH3CHCI4-CH3
=N02+OH
=H02+M
= C2H303
= HN03
=H02*02
=H+C02
= CH3r)2*M
=CH30+N02
=CH3G2*N02
= C2H3U2+NU2
=C2H30+N02
= C.H3CHQ*N02
= HCHO + H(J2
=CH3+C2H30
=HCMU*C2H402
= HCHO-»-C2HM02
= C2H303+nH
=CH30+CH3CHO
•CH3»MEO»C2H30
= CH3+C(!+p''l
=CH30N02
C2H303 + N02 = C2H3C13N02
DUMHC*0
HOMO » HV
=CH3*C2H30
= H + h * CO
N02 +N03 + H20 * HNOJ*- "
1
23
3
1
25
24
4
2
14
20
I
1?
15
17
19
18
2
21
22
18
4
^
4
22
4
4
20
18
1
12
7
1
0
0
2
3
0
2
15
16
0
0
15
3
3
0
2
2
25
2
2
0
23
3
23
0
16
19
0
1
25
23
0
24
2
3
1
24
21
1
6
1
16
25
5
16
14
19
18
19
1
20
6
7
17
7
7
15
18
17
17
9
10
17
14
5
23
0
0
0
3
1
17
15
0
0
0
0
0
0
1
I
21
1
1
16
20
22
22
25
6
18
8
0
0
20
14
5
0
0
6
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
• o
0
20
0
0
0
0
a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
.0
0
0
0
0
0
0
.40
.14
.4
.1
.2
.245
,1
.1
.1
.1
.1
.1
.3
.14
.2
.1
.25
.2
. 1
.2
.350
.5
.175
.3
.1
.1
.1
.1
.2
.70
.3
.1
0
D
0
00
07
00
D-04
D
0
D
0
D
D
0
D
0
0
D
0
D
n
D
0
D
02
03
04
02
07
06
03
02
01
07
01
02
01
02
03
03
02
D-04
1)
0
D
02
03
01
D-Ol
D
D
D
0
02
01
00
01
0-02
0
03
3
4
5
6
7
8
'I
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
Card Type 15
.1 D 03
.5
0-01
Card Type 16
33.5
17.7
50.0
20.0
97.5
18.9 .
100.6
23.0
289.6
23.0
l'^9. 4
22.7
591.3 11
27.7 ;
359.7 ii|
23.9 !
0.4
5.5
'O.I '
8.5
1517.9
22.4
1527.1
22.7
-------
VII. 1
Section VII
Glossary
A. COMMON REGIONS
1. KINCOM
A Heading label
BB Numerical integration counter
C Pollutant reaction rates
CC Plot data
D Rates of reaction
DERY Rates of change of species
INTR Output list indicator
ITRAF—Execute traffic sources function indicator
KDIF NDIF+2
LABEL—Species names
M Page number of output
N Temporary storage
NC No. of species
NDIF No. of accumulating species
NEQ No. of mechanism steps
NLL No. of plot lines
NN No . of printouts
NflC NC-NDIF
NNDIF—NDIF+1
PEAK1—Peak value of concentration 1
PEAK3—Peak value of concentration 3
PRMT Parameter input list
RCT Reaction list
STEP Output list
TIM1 Half time of concentration 1
TIM2 Half time of concentration 2
TIM4 Half time of concentration 4
TIME1—Peak time of concentration 1
TIME3—Peak time of concentration 3
XPOS X coordinate (miles)
YPOS Y coordinate (miles)
YY Non-accumulating species
ZZ Chain yield
-------
VII. 2
2. MYCOM
AIR Meteorological data
AGRIDl-Area source grid distribution of reactive hydrocarbons
AGRID2-Area source grid distribution of nitric oxide
AGRIDA-Area source grid distribution of less reactive hydrocarbons
ARADIO-A.M. radiosonde curve
EMIT Emission factors for vehicular sources
FGRID—Grid distribution of freeway traffic
ICLOUD-Type of clouds in sky
IDM First time through METEOR module indicator
IDS First time through SOURCE module indicator
IDU First time through ULTRAV module indicator
IHT Mixing height values input indicator
IPST Pacific standard time indicator
NRAD-—No. of radiosonde points
NSTA No. of meteorological stations
POS Stationary source position in meters
PRADIO-P.M. radiosonde curve
SGRID—Grid distribution of street traffic
TIME Time of day
WPOS Meteorological stations position and elevation
XDOT X component of wind
XHITE—Array of input mixing height values
XMIN Minimum mixing depth
YDOT Y component of wind
ZENITH-12 hourly values of sun angle at Los Angeles
3. ALLCOM
ANS Array of calculated meteorological values for the trajectory point
GOTIME-Simulation start time
GPOLL—Source inputs in gram-moles/minute
ISOUR—Execute source module indicator
POLL Source inputs in pphm/minute
XDEPTH-Mixing depth
XDEPRV-Previous mixing depth
4. BARER
IDB First time through barrier check routine indicator
-------
VII. 3
B. MAIN
1. Arrays
CSAVE—Holds correction terms and derivatives for DIFSUB
ERROR—Contains estimated one step error
IDMET—Pollutant station names
IDNUM—Meteorological station names
IP Temporary storage space for pivots
PW Temporary storage space
SAVE Saves matrix Y
Y Dependent variables and their scaled derivatives
YMAX Maximum values of Y
YSAVE—Temporary storage for Y
S calars
ATEMP—Temporary temperature
BAROM—Barometric pressure
DEWPT—Dew point
IN Another simulation run indicator
KFLAG—DIFSUB completion code indicator
MAXDER-Maximum derivative number
MF Method indicator for Gear
NDIM Order of set of ordinary differential equations
C. SMGOUT
1. Arrays
Y Dependent variables of differential equations (accumulating
and non-accumulating species)
2. Scalars
H Step size
HAL1 Half the current value of concentration 1
HAL2 Half the current value of concentration 2
HAL4 Half the current value of concentration 4
PREVYl-Previous value of concentration 1
PREVY2-Previous value of concentration 2
PREVY4-Previous value of concentration 4
PTIM1—Time when PREVY1 is saved
PTIM2—Time when PREVY2 is saved
PTIM4—Time when PREVY4 is saved
X Time into simulation
-------
VII. 4
D. DIFFUN
1. Arrays
DY First derivatives of accumulating species
Y Set of accumulating and non-accumulating species
YSAVE—Temporary storage for Y
2. S calars
STORE—Temporary rate of change of mixing depth
X Time into simulation run
XPREV—Previous simulation time last step
E. METEOR
1. Arrays
DISINV-Sum of the square of the inverse distances from the air parcel
location to the meteorological stations
DIST Square of the distance from the air parcel point to a meteoro-
logical station
JBAR Permissible station indicators
RADIO—Radiosonde curve in use
RSLOPE-Slopes of the line segments of RADIO
VAL1 Spatially interpolated meteorological values at first hour
VAL2 Spatially interpolated meteorological values at second hour
2. Scalars
HMIN Minimum elevation point on the radiosonde curve
HOUR Hour of the day
ITIME—Hour of the day minus 6
SLOPE—Adiabatic lapse rate
TIMEHR-Number of minutes into current hour
TMIN Minimum temperature point on the radiosonde curve
X X-coordinate position of the pollutant column
XPREV—Time previous cycle through routine
XTIME—Time into simulation run in minutes
XTIMEl-Temporary storage of current time
XXDOT—X component of rate of change of wind
Y Y-coordinate position of the pollutant column
YYDOT—Y-component of rate of change of wind
-------
VII. 5
F. ULTRAV
1. Arrays
ABSORP-Absorption coefficients of N02 and N20^
AIRMS—Air mass at various zenith angles
CLOUD—Average fractions of incident radiation transmitted by 7
cloud types
FLIST_-Dependent variable list corresponding to independent values
in XLIST_
OZONE—Absorption coefficients of ozone in the ultraviolet region
by wavelength, 3000A-4000A
SOLAR—Mean solar spectral irradiance outside, the atmosphere for
spectral bandwidths of 100A centered at wavelengths 3000A-4000A
VAPOR—Values of water vapor pressure for different temperatures
XLIST_-Independent variable interpolation list
XMOLE—Indices of refraction and molecular scattering coefficients
for a homogeneous standard atmosphere at 3000A-4000A
2. Scalars
AIRMAS-Values of the air mass at various zenith angles
ALBEDO-Reflectively of the ground
ANGLE—Zenith angle of the sun
BETA Value of air mass for given zenith angle
CANGLE-COS of zenith angle
Cl Temporary storage
C2 Temporary storage
DIRECT-Direct radiation from the sun in the smog layer
DUST Dust content of the atmosphere; 1 near surface of urban areas,
0 in upper atmosphere
G Parameter between 0 or 1 which is a function of the directional
distribution of the scattered radiation and of the amount of
multiple scattering
LAMDA—Wavelength in Angstroms
OZ3 Quantity of ozone in the atmosphere at standard temperature
and pressure
PI 3.14159
PRESS—Water vapor pressure
PZ Partial pressure of atmospheric water vapor
RATE Specific ultraviolet absorpiton rate of N02
RATE1—Temporary storage
RATE2—Temporary storage
SCATOZ-Transmissivity with respect of ozone absorption
SCTMOL-Transmissivity relative to molecular scattering
SCTPAR-Transmissivity due to particle scattering
SKY Radiation reflected from the sky into the smog layer
-------
VII. 6
G. SOURCE
1. Arrays
DIST Square of the distance from the air parcel position to the
midpoint of the three closest grid areas
EMITWT-Temporary storage
EPOS Midpoints of three closest grid areas
FCURVE-Temporal distribution curve of freeway traffic
IBOX Number of grid containing the air parcel
SCURVE-Temporal distribution curve of street traffic
WEIGHT-Molecular weight of major pollutants
2. S calars
AREA!—Reactive hydrocarbon emissions in grid area containing parcel
(kg./hr.)
AREA2—NO emissions in grid area containing parcel (kg./hr.)
AREA4—Less reactive hydrocarbon emission in grid area containing
parcel (kg./hr.)
AX Value of diurnal area source usage curve at any time during
the day
DISINV-Sum of the square of the inverse distance from the air parcel
location to the midpoint of an adjoining grid square
FACTOR-Conversion factor inversely proportional to mixing depth
FWY Daily freeway mileage in grid square containing air parcel
FX Value of diurnal curve for freeway traffic at any time during
the day
GMIN Minutes into current hour of the day
IHR Current hour of the day
IPT Grid square quadrant indicator
IXM X-coordinate of midpoint of grid square containing air parcel
IYM Y-coordinate of midpoint of grid square containing air parcel
STREET-Daily street traffic mileage in grid square containing air parcel
SX Value of diurnal curve for street traffic at any time during
the day
X X-coordinate of air parcel location
Y Y-coordinate of air parcel location
-------
Section VIII
Program Listing of Reactive Pollutant
Environmental Simulation Model (REM)
-------
VIII. 1
COMPILER
C
C
OS/360 FORTRAN H
OPTIONS - NAME' MAIN* OPT »_01»yNEJNT»55,S_IZE»OppOJ^»
SOURC E, BC Of NOL I Sf »"NOOEC K» LQIAD* MA?, NOEDIT, 10* NOXREF"
REACTIVE ENVIRONMENTAL MODEL _. VERSION CB
DIMENSION SAVEJ 13*40) »C SAVE U0» 3), YMAXUO), ERROR UO),PW< 1600),
1 IP(40),Y(13*40)WSAVE(40), IDNUM(32V»lbMET(32)*IR(50»5),
2 D A I R (12» 32 » 2 ) » RS OR IDJ 650) » RFGR I0(650). RGR I DJ (650) »
3 RGRID2(650)VRGR'iD4(650J ".. " "
DOUBLE PRECISION C,D»Y,YY,ZZ»DERY,SAVE*CSAVE,YMAX,ERROR*YSAVE»
I POLL»PRMT»BB*Z
LABEL»BLANK/« '/,RCT
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
t
C
C
C
C
C
C
C
C
C
C
C
C
C
KINCOM KINCOM KINCOM KINCOM KINCOM
A HEADING LABEL
8B SAVE PLOT DATA COUNTER
C POLLUTANT REACTION RATES
CC PLOT DATA
"D— —-RATES OF REACTION
RATES OF CHANGE OF SPECIES
-TYPE OF DISPERSION CURVES
OUTPUT LIST INDICATOR
SOURCE EMISSIONS MODULE INDICATOR
TRAFFIC EMISSIONS MODULE INDICATOR
IDISP
INTR-
ISOUR
ITRAF
LABEL-SPECIES NAMES
FOR NUM. INT
-NO. OF COLUMNS
-NO. OF SPECIES
-NO. OF ACCUMULATING SPECIES
-NO. OF MECHANISM STEPS
•NJ. OF PLOT LINES
-NO. OF PRINTOUTS
-NC-NDIF
•NDIF+1
-PEAK VALUE OF CONCENTRATION 1
-PEAK VALUE OF CONCENTRATION 3
-PARAMETER INPUT LIST
-REACTION" LIST
-OUTPUT LIST
-HALF "TIME OF CONCENTRATION 1
•HALF TIME OF CONCENTRATION 2
•HALF TIME OF CONCENTRATION 4
TI ME 1-PEAK TIME OF C£NCE^NTR_^T_IOJN 1
T I ME3-PE AK " f I ME"' OF ~C ONCENTRAT I ON 3 "
XPOS—X COORDINATE (MILES)
YPOS--Y COORDINATE (MILES)
YY NON-ACCUMULATING SPECIES
ZZ — -CHAIN YIELD
~ND IM- -NO. OF AC C lJM"uI A"fTNG~ S ?£ CIIS ""*' 2 ~
f-, CONCENTRATIONS OF SPECIES
COMMON /KINCOM/ Z2(4)»YY(40),LABEL(40)*D(50),RCT(50»3).PRMT(8)
NC
NOIF-
NEQ —
.NJL.L™
NN
NNC —
NNDIF
PEAK1
PEAK3
PRMT-
"RCT--
STEP-
TIM1-
TIM2-
TIM4-
00001000
00002000
00003000
00004000
00005000
00005100
00006000
00007000
00008000
00009000
00010000
00011000
00012000
00013000
00014000
00015000
00016000
00017000
ooouooo
00019000
00020000
00021000
00022000
00023000
00024000
00025000
00026000
00027000
00028000
66029000
00030000
00031000
00032000
00033000
00034000
00035000
00036000
00037000
00038000
00039000
00040000
00041000
00042000
00043000
00044000
00045000
00046000
6oo47b"oo
00048000
00049000
00050000
-------
VIII.2
1 OERY < 40 1 >C < 50 ) » 88 ? ZJ 2J » N» NC_»NEQ» N0.1f ^KOIF » M, NN, INTR, NNC »
2 NNO'lF, I UV> IM» IS»ITRAF ,NLL"icC< 16, 100 ) » A("«61 *STEP (15 ) ,T~IM~1,
3 TIM2,TIM4,PEAKl,TIMEi,PEAK3*TIME3*XPOS,YPOS
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
MYCOM MYCOM MYCOM-- -MYCOM. MYCOM—
EMISSIONS
EMISSIONS
EMISSIONS
AGRID1-AREA SOURCE PROPYLENE
AGRID1-AREA SOURCE NO
AGRID1-AREA SOURCE OUMHC
AIR METEOROLOGICAL DATA
EMIT TRAFFIC MODULE EMISSION RATES
ICLOUD-TYPE OF CLOUDS IN SKY
I DM FIRST TIME THROUGH METEOR MODULE
IDS FIRST TIME THROUGH SOURCE MODULE
IDU FIRST TIME. THROUGH ULTR_AV MODULE
IHT MIXING HEIGHT"VALUES INPUT IWHCATOR
IPST PACIFIC STANDARD TjME INDICATOR
NRAD NO. OF RADIOSONDE POINTS
NSTA NO. OF M_ETEOR_QLC)G.ICAL STATIONS
NSTACK-NO. OF STATIONARY SOURCES
TINJL TLMf OF PAT
"STATIONS"
INDICATOR
INDICATOR
INDICATOR
WPOS METEOROLOGICAL
•X COMPONENT OF
•ARRAY OF INPUT
•MIXING HEIGHT
•Y 'COMPONENT OF
POSITION AND ELEVATION
MIXING HEIGHT VALUES
XHITE —
KHT
YDOT-—Y 'COMPONENT OF WIND
ZENITH-12 HOURLY VALUES OF SUN ANGLE AT LOS ANGELES
COMMON /MYCOM/ TIME»ICLOUD,IDM»IDU,IDD*IDS^IPST.NSTA,NSTACK»NRAD»
1 AIR ( 12» 32. 6 ) » WPOS ( 32. 3) » POS ( 32» 3 )»S>OS ( 32» 3 jit HS (12 ) »
2 DS(12),VS(12)»TS<12)»QS( 12) ,S8AROM» STEM_P» XDOT» YDOT»
3 AR AD I 0 (15. 2 ) » PR AD 10"( 15,2)» ZENi TH ( 12"» 2 )» SGR ID ( b50 J ,
* _f.9R15 < 6-5P' » A?R1PJ (&50 ) » A_GR^p2 (65plf_AGRJJA (650 )* EJIj T (10) »
5 """ XHTTE(12)VVH"f»XMIN "~ " " "~ ~ "
COMMON" /BARRR/ IDB
IDB FIRST TIME THROUGH BARRIER CHECK INDICATOR
ALL COM ALLCOM ALLCOM ALLCOM ALLCOM ALLCOM
ANS ARRAY OF CALCULATED METEOR." VA'CuTsTTOR THE TRAJ. PT.
GOTIME-SIMU_LATION START TIME
ISOUR —SOURCES MODULE INDICATOR
POLL ARRAY OF CALCULATED TRAFFIC AND ARE_A SOURCES
XDEPThf-MIXING^bEPTH AT TRAJECTORY POINT
XOEPRV-LAST CALCULAIED ..MIXING DEPTH
COMMON /ALLCOM/ POLL(4,3),GPOLL(4,3)»ANSC8),GOTIME,XDEPTH,ISOUR»
OOJD51000
00052000
00053000
00054000
00055000
00056000
00057000
00058000
00059000
00060000
00061000
00062000
00063000
00064000
00065000
00066000
00067000
00068000
00069000
00070000
00071000
00072000
00073000
00074000
00075000
00676000
00077000
00078000
00079000
oooeoboo
oopjaipoo
"o ob a 2000
00083000
00084000
00085000
00066000
00087000
6o6~88000
00089000
00090000
00091000
00092000
_0 009 3000
00094000
00095000
00096000
00097000
OOOT98000
00099000
boTdooibo
00101000
00162600
00103000
-------
VIII.3
C
C
C
C
C
C
C
C
C
C
C
EXTERNAL OIFFUN
?H(ATEMP.DEWPT)
1
XDEPRV
10.*M7.5*(DEwPTMDEWPT>.237.3)
-ATEMPMATEMP + 237.3) ) )
20
30
INITIALIZATION
INPUT HEADING INFORMATION
IN = 0
READ <5,1000»SND=900) A
JSTART = 0
00 10 1=1,40
00 5 J=l,13
Y( Jf I ) = O.DO
YMAX(I) = 1.000
YSAVEU ) = 0.00
LABEUI) = BLANK
D< I ) = O.ODO
DERY(I) = O.ODO
30 20 1=1,4
ZZ< I) = 0.000
^ = 0
M = 1
NN = 1
88 = 1.0E08
108 = 0
100 = 0
ion - o
OS = 0
DU = 0
PST = 0
UV = 0
M = 0
S = 0
00 30 1=1,4
00 30 J=l,3
3POLLU »J) = 0.
POLL( I» J) = 0.DO
BAROM
:OVER
1013.
0.
INPUT START POSITION. START TlflE, MINIMUM XDEPTM
INPUT LENGTH OF SIMULATION, HMIN, HMAX, INITIAL STEP SIZE, ERROR
CRITERION, PLOT VALUE JNOEX"
READ (5,1001) XPQS.YP'DS»GOTIME»XMIN
?EAD (5,1006) (PRMT( I)»I»2,7)
00103200
00104000
00105000
00106000
00107000
00103000
00109000
00110000
00111000
00112000
00113000
00114000
00115000
00116000
00117000
00118000
00119000
00120000
00121000
00122000
00123000
00124000
00125000
00126000
00127000
00128000
00129000
00129100
00130000
00131000
00132000
00133000
00134000
00135000
00136000
00137000
00133000
00133100
00139200
00139000
OOl'-OOOO
00141000
00141200
00142000
OOU3000
00144000
00145000
00146000
00147000
00143000
00149000
00150000
00151000
-------
VIII.4
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
IF (PRNT(7>) 70,70,75
70 PRMTm = 2.DO*PRMT(2)
75 PRMT( I) = 0.
PPMT(8) = PRMT(l)
INPUT NO. OF ACCUMULATING SPECIES, TOTAL NO. OF
MECHANISM STEPS, PLOT LINE PARAMETER, NO,
00152000
00153000
00154000
00155000
00156000
SPECIES, NO. OF 00157000
OF METEOROLOGICAL 00158000
STATIONS, NUMERICAL INTEGRATION METHOD, NO. OF DERIVS, CLOUDOO 1 59000
TYPES, NO. OF RADIOSONDE POINTS, TRAFFIC
MODULE, MIXING HEIGHT, SAME DAY FOR NEXT
STANDARD TIME INDICATORS
?EAD (5,1002) NDIM,NC»NEQ,NLL»NSTA,MF*MAXDER» 1C
1 ISOUR.IHT, IN, IPST
INPUT OUTPUT INDICATOR AND LIST OF APPROXIMATE
READ (5,1003) INTR,STEP
SOIF = NOIM
>4DIM = NDIM+2
KDIF = NDIM
NNDIF = NDIF+l
INPUT MIXING HEIGHT
IF (IHT . NE. 1 ) GO TO 80
READ (5,1020) XHITE
INPUT SPECIES LIST AND INITIAL CONCENTRATIONS
CONCENTRATIONS ARE IN PPHM EXCEPT FOR CD
60 READ (5,1004) ( L AB EL < I ) , Y ( 1 , I ) , I - 1 . NC )
SAVE NON-ACCUMULATING SPECIES
DO 85 I=NNDIF,40
YY( I ) = Yd, I )
85 Y( 1,1 ) = O.DO
CHECK FOR SECOND RUN
IF (IN .EQ. 1) GO TO 300
INPUT SUN_ANGLE AT LOS ANGELES
READ (5,1032) ( ( ZEN 1 TH ( I , J ) , J = l , 2 ) , I = 1, 12 >
INPUT METEOROLOGICAL DATA.
MODULE., SOURCE 00160000
RUN, AND PACIFIC 00161000
00162000
00163000
LOUD»NRAD, ITRAF, 00164000
00165000
00166000
OUTPUT TIMES 00167000
00168000
00169000
00170000
00171000
00171200
00172000
00173000
00174000
00175000
00176000
00177000
00178000
00179000
IN PPM 00179100
00180000
00181000
00162000
00183000
00134000
00165000
00166000
00187000
00168000
00189000
00190000
00191000
00192000
00193000
00194000
0019500J
00196000
00197000
00198000
00199000
00200000
00201000
00202000
-------
VIII.5
100
c
c
c
c
c
c
c
c
c
105
10?
120
130
DO 100 J=1,NSTA
READ (5,1005) I DNUM ( J ) , ( ( A I R ( I » J, K ) , K = 1,4 ) , I =1 , 4 )
READ (5,1005) IONUM( J) »( (AIRU, J,K),K*1,4) ,1=5,8)
READ (5.1005) I DNUM (J ) , ( ( A IR ( I'., J,K ) * K= 1, 4 ) , I = 9, 1 2 )
INPUT METEOROLOGICAL STATION OATA
WRITE (6,1013)
READ (5.1014) (IDMET(I),(WPOS(I,J).J=i,3),I=1,NSTA)
CONVERT OATA TO METRIC UNITS
DO 120 I»1,NSTA
*POS(I..3) = WPOS(I,3)*.3048
DO 120 J=l,12
IF (AIR(J,I,3) .EQ. 0) GO TO 105
AIRU.I.3) = (AIRU. I.3)-32. )*. 55555
IF (AIRU.I.4) .EQ. 0) GO TO 120
AIR(J,I,4» = (AIRU, I.41-32. )*. 55555
CONTINUE
DETERMINE RELATIVE HUMIDITY.
150
160
165
170
ISO
0.) GO TO 140
00 150 1*1, 12
DO 150 J=1»NSTA
IF (AIR(I.J»3) .NE,
AIR(I»J»4) = 0.
GO TO 150
IF (AIR(I»J*4) .EO. 0.) GO TO 150
ATEMP = AIRU» J»3)
OEWPT = AIR(I.J.4)
AIR(I,J,4) = RH(ATEMP.DEWPT)
:D^4TINUE
OUTPUT METEOROLOGICAL DATA
C
C
WRITE (6*1015)
IF (IHT .NE. 1)
WRITE (6,1021)
KOUT = 0
WRITE (6,1007)
DO 180 1=1,12
IH = (I+6)*100
DO 180 J = 1.,NSTA
IF (KOUT .LE. 50)
KOUT = 0
WRITE (6,1007)
GO TO 180
KOUT = KOUT+1
WRITE (6,1008)
SAVE WIND DATA
( IDMETU ), (WPOS( I » J)» J*l,3),
GO TO 165
XHITE
l.NSTA)
GO TO 170
IDNUM(J)»IH,(AIR(I,J»K),K=l,4)
00203000
00204000
00205000
00206000
00207000
00208000
00209000
00210000
00211000
00212000
00213000
00214000
00215000
00216000
00217000
00218000
00219000
00220000
00221000
00222000
00223000
00224000
00225000
00226000
00227000
00228000
00229000
00230000
00231000
00232000
00233000
00234000
00235000
00236000
00237000
00233000
00239000
002^0000
00241000
00242000
00243000
00244000
00245000
00246000
00247000
00248000
00249000
00250000
00251000
00252000
00253000
00254000
00255000
-------
VIII.6
200
C
C
C
232
235
240
233
250
C
C
C
00 200 J=l,12
DO 200 I=l,NSTA
DAIR
-------
VIII. 7
280
300
320
READ (5,1038) Z
INPUT AM AND PM RADIOSONDE CURVES
;F (IHT .EQ. 1) GO TO 280
READ (5.1022) I ARAO I 0(I,1),I = 1.NRAD)
READ (5,1022) (ARAD I 0(I,2),I = I,NRAD)
(PRADIO'J I,"l), 1 = 1,NRAD)
(PRADIO(1,2),I»l,NRAD)
(ARADIO( I,I),ARAD 10!1.2).1 = 1,NRAD)
(PRADn(I,l),PRADIO(I,2),I=l»NRAD)
( (ZENITHd, J), J=1,2),I = 1,12)
READ
.READ
WRITE
WRITE
WRITE
IN *
WPITE
WRITE
DO 320
WRITE
WRITE
WRITE
WRITE
(5,1022)
(5,1022)
(6,1023
(6, 1024
(6, 1033
I
6, 1025
6, 1010
(PRMT(
I=2»7),GOf IME,MF,MAXDER»A
I=1»NEQ
6,1012) I,(RCT(I,J),J=1,3)»C(I),(IR(I,K),K=1,5)
6,1039) Z
6, 1031 >
6,1037) (I,LA8EL(I),Y(1,I),DERY(I),1=1,NC)
C
c
C
SET UP POSITION FOR INTEGRATION
Y(l»NDIF+l) = XPQS
Y(l,NOIF+2> = YPOS
400 CALL DIFSU8 (NDIM,PRMT(1),Y,SAVE,CSAVE,PRMT(5),PRMT(3 ),PRMT(4 ) ,
I PRMT(6),1F,YMAX»ERROR,KFLAG,JSTART*MAXDER»PW»IP)
IF (KFLAG .NE. 1) GO TO POO
DO 420 J=l,NC
420 YSAVE(J ) = Y( 1,J)
CALL SMGOUT (PRMT<1),YSAVE,°RMT<5 ) )
JSTART =1
IF (PRMT(l) .LT. PRMT(2» GO TO 400
IF (NEQ .EQ. I) GO TO I
CALL SMGPLT
30 TO 1
900 WRITE (6,1026) KFLAG
STOP
1000 FORMAT(20A4)
1001 FOR1AT(4F10.3)
1002 FORMATl 1415)
1003 f=ORMAT( 12,3X, 15F5.2 )
1004 FORMAT(4(A8,012.3))
1005 FORf1AT(A3,3X.4(F3.0,F3.1.F3.0,C2.o,6X),6X)
1006 FORMAT(6F10.3)
1007 FORMAT( 1H1, 'STAT I ON•,6X,•0 AT A•,6X,•WINO1»6X,•WIND'»16X,•REL.',
1 /' ','NUMBER',6X,'HOUR'»6X,'OIR. ',5X,'SPEED1 .6X, 'TEMP' ,
2SX,'HUM.',
/• ',22X,MDE<
•4X,
3' (CENT)', / /)
1008 FORMAT(« •,A7,6X,14,6F 10. 3)
1009 FORMAT(25F3.0)
1010 FORMAT!'OK-VALUES FGR REACTIONS1//)
1011 FORMAT (3A8,5 13,015.8.26X)
00295100
00296000
00297000
00298000
00299000
00300000
00301000
00302000
00303000
00304000
00305000
00306000
00307000
00308000
00309000
00310000
00311000
00311100
00312000
00313000
00313200
00313220
00313240
00313260
00313280
00314000
00315000
00316000
00317000
00318000
00319000
00320000
00321000
00321200
00322000
00323000
00324000
00325000
00326000
00327000
00328000
00329000
00330000
00331000
00332000
00333000
00334000
00335000
00336000
00337000
00338000
00339000
00340000
-------
VIII.8
1012 FORMAT!' •.I4,5X,3A8,020.S»5X»5 I 3) 00341000
1013 FORMAT! IH1.30X,'METEOROLOGICAL STA~T I ONS • , 1IHO, 2 1 X, ' X-COORO • > 11 X* 00342000
I1Y-COORD',10X,«ELEVATION'./' ',3X,'STATI ON',11X>•(MILES)'>1IX, 00343000
2' !MILES)',HX, '(METERS) '/) 00344000
1014 FORMAT!4( A3,2F5.2,F7.2) ) 00345000
1015 FORMAT!' •,A 10,2F13.3.F19.3) 00346000
1020 FORMAT! 12F6.0) 00347000
1021 FORMAT!//////' INPUT MIXING HE IGHTS---0700 'HOURS TO 1800 HOURS',//00348000
110X,12F9.1) 00349000
1022 FORMAT!8F10.I) " 00350000
1023 FORMAT!1H1.20X,'AM RAOIOSDNOE CURVE'»///1H ,23X,'HEIGHT',6X,•TEMP'00351000
1./1H ,21X,'{METERS) (CENT.)«,//(IH , 19X,If 10.1) I 00352000
1024 FORMAT! 1HO»20X.'PM RADIOSONDE CURVE',///1H ,23X,'HE IGHT',6X»'TEMP'00353000
1,/1H »21X, ' JMETERS) CCENT. )'..//! 1H , 19X, 2F 10. 1 ) ) 00354000
1025 FORMAT!'I1,010.4,' MIN SIMULATI ON'»4X,•HMIM=•,010.4,4X,«HMAX=•, 00355000
1 010.4,AX,'INITIAL H=',D10.4/30X,'ERROR CRI TER I A* '., 010. 4, I 3X, 00356000
2 '88=',010.4//1 START TIME IS'»F8.4»' HOURS',21X,'MF=•,I 2,14X» 00357000
3 'MAXOER='» I3//4! ' ',20A4,/» 00357100
1026 FORMAT!//////• KFLAG= »,I3) 00358000
1027 FORMAT! 10F5.2) " 00359000
1028 FORMAT!//////• TRAFFIC EMISSION RATES',//1 PPQPYLENE=',F5.2,10X, 00360000
l'NO='»F5.2,IOX,'CO=',F5.2»10X»'DUMHC='»F5,2) 00361000
1029 FORMAT!1H1,'GRID DISTRIBUTION OF FREEWAY TRAFFIC IN THOUSANDS OF V00362000
1EHICLF MILES PER DAY'» / /(25F5,0) > 00363000
1030 FORMAT11H1,'GRID DISTRIBUTION OF STREET TRAFFIC IN THOUSANDS OF VE00364000
HICLE MILES PER DAY« / /(25F 5.0) ) " " 00365000
1031 FORMAT! IH1,'STARTING VALUES•///,8X,' COMPOUND CONCN PPHM CHA003&6000
INGE. PPHM/MIN'/) 00367000
1032 FORMAT! 3!2F5. l» 00368000
1033 FORMAT!'0',20X,'ANGLE OF THE SUN AT LOS ANGELES',//'0',25X,•HOJR',00369000
I 5X,'ANGLE',//(1H . 1 9<,2F 10.1)) 00370000
1034 FORMAT I1H1,'GRIO 0 ISTRI BUT I ON "OF REACTIVE HYDROCARBONS (PROPYL00371000
IENE) IN KILOGRAMS PER HOUR',//!25F5.0)) 00372000
1035 FORMAT! 1HI,'GRID DISTRIBUTION OF NO IN KILOGRAMS PER HOUR'»//(25F500373009
1.0)) 00374000
1036 FORMAT!IH1,'GRID DISTRIBUTION OF LESS -RE ACT IVE HYDROCARBONS (DUMHC00375000
1) IN KILOGRAMS PER HOUR',//!25F5.0)) 00376000
1037 FORMAT!' ',5X, I 3, 1 X,A8,20 16.7 ) 00377000
103S FORMAT (2D20.8) 00377100
1039 FORMAT
-------
VIII.9
OS/360 FORTRAN H
COMPILER OPTIONS - NAME= MA IN,OPT = Ol,LINECNT=55,SIZE = OOOOK*
SOURCE,BCD,NOLIST,NaDECK,LOAD,MAP,NOEDIT,'rD,NOXREF
SUBROUTINE DIFSUB(N*T,Y»SAVE,CSAVE*H,HMIN,HMAX*EPS*MF,YMAX,
1 ERROR,KFLAG,JSTART,MAXDER,PW,IP)
IMPLICIT REAL*8 (A-H.Q-Z)
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$
THIS SUBROUTINE INTEGRATES A SET OF N ORDINARY DIFFERENTIAL FIRST
ORDER EQUATIONS OVER ONE STEP OF LENGTH H AT EACH CALL. H CAN BS
SPECIFIED BY THE USER FOR EACH STEP, BUT IT MAY BE INCREASED OR
DECREASED BY DIFSUB WITHIN THE RANGE HMIN TO HMAX IN ORDER TO
ACHIEVE AS LARGE A STEP AS POSSIBLE WHILE NOT COMMITTING A SINGLE
STEP ERROR WHICH IS LARGER THAN EPS IN THE L-2 NORM, WHERE EACH
COMPONENT OF THE ERROR IS DIVIDED BY THE COMPONENTS OF YMAX.
THE PROGRAM REQUIRES FOUR SUBROUTINES NAMED
DIFFUN(T,Y,DY)
DECQMP(N,M,PW,IP)
SOLVE(N,M,PW,CSAVEd, 1)* IP)
PEDERV(T,Y,PW,M)
THE FIRST, DIFFUN* EVALUATES THE DERIVATIVES OF THE DEPENDENT
VARIABLES STORED IN Yd,I) FOR I = 1 TO N, AND STORES THE
DERIVATIVES IN THE ARRAY DY. THE NEXT TWO ARE CALLED ONLY IF
METHOD FLAG MF IS SET TO 1 OR 2 FOR STIFF METHODS. DECOMP IS
DECOMPOSITION WITH PIVOTING THAT DECOMPOSES THE MATRIX
THE PIVOTS IN THE INTEGER ARRAY IP. M IS THC DECLARED
IP(N) IS SET TO 0 IF PW IS SINGULAR SOLVE PERFORMS
BACK SUBSTITUTION ON THE CONTENTS OF CSAVEU,1), LEAVING THE
RESULTS THERE. THE MATRIX IN BOTH SHOULD BE SINGLE PRECISION,
WHILE ALL FLOATING POINT VECTORS SHOULD 3E DOUBLE PRECISION.
PEDERV IS USED ONLY IF MF IS I, AND COMPUTES THE PARTIAL
DERIVATIVES OF THE DIFFERENTIAL EQUATIONS AS
MF PARAMETER.
THS
A
STANDARD LU
PW, LEAVING
SIZE OF PW.
DESCRIBED UNDER THE
THE PROGRAM USES DOUBLE PRECISION ARITHMETIC
POINT VARIABLES EXCEPT THOSE STARTING WITH P
SINGLE PRECISION TO SAVE TIMF AND SPACE.
FOR ALL FLOATING
THE FORMER ARE
SEVERAL ARRAYS.
CSAVE U,l)
CSAVE(I,2)
0037900Q
00380000
00361000
*00382000
*00363000
*00334000
*003fl5000
*00386000
*00387000
*00388000
*00369000
100390000
*00391000
*00392000
*00393000
$00394000
*00395000
*00396000
*00397000
*00398000
$00399000
$00400000
*00401000
''00402000
$00403000
$00404000
*00405000
$00406000
$00407000
$00408000
$00410000
$00412000
13000
TH? TEMPORARY STORAGE SPACE IS PROVIDED BY THE CALLER IN THE
INTEGER ARRAY IP, THE SINGLE PRECISION ARRAY PW, AND THE DOUBLE
PRECISION ARRAYS SAVE AND CSAVE. THE ARRAY PW IS USED ONLY TO HOLD
THE MATRIX OF THE SAME NAME* AND SAVE IS USED TO SAVE THE VALUES
OF Y IN CASE A STEP HAS TO BE REPEATED, BUT CSAVE IS USED TO HOLD
THE REGIONS USED ARE
IS USED MAINLY TO HOLD THE CORRECTION TERMS IN THE
CORRECTOR LOOP, AND HOLDS THE DERIVATIVES DURING
JACOBIAN EVALUATIONS.
IS USED TO SAVE THE VALUES OF THE SUMS OF ALL OF THE
CORRECTION TERMS IN THE PREVIOUS STEP AFTER THE_Y
HAVE BEEN ACCUMULATED IN THE ARRAY ERROR IN~THE
CURRENT STEP. THIS ENABLES THE BACKWARDS DIFFERENCE
OF ERROR TO" BE FORMED. IT IS USED TO ESTIMATE THE'
STEP SIZE FOR ONE ORDER HIGHER THAN CURRENT.
*00417000
*004,?0000
*00422000
*OOA23000
$00426000
*00427000
*00428000
*00429000
-------
VIII.10
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*
c*
c*
c*
CSAVEU>3) IS USED TO STORE THE DERIVATIVES WHEN THEY ARE
COMPUTED BY DIFFUN.
THE
THE
N
PARAMETERS TO THE SUBROUTINE DIFSUB HAVE
FOLLOWING MEANINGS..
N
SAVE
CSAVE
H
HMIN
HMAX
EPS
MF
THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS,
MAY BE DECREASED ON LATER CALLS"IF THE NUMBER OF
ACTIVE EQUATIONS REDUCES., BUT IT MUST NOT BE
INCREASED WITHOUT CALLING WITH JSTART = 0.
THE INDEPENDENT VARIABLE.
A 13 BY N ARRAY CONTAINING THE DEPENDENT VARIABLES AND
THEIR SCALED DERIVATIVES. Y(J-H,I) CONTAINS
THE J-TH DERIVATIVE OF Yd) SCALED BY
H$$J/FACTORIAL(J) WHERE H IS THE CURRENT
STEP SIZE. ONLY Y(l,[) NEED BE PROVIDED BY
THE CALLING PROGRAM ON THE FIRST ENTRY.
IF IT IS DESIRED TO INTERPOLATE TO NON MESH POINTS
THESE VALUES CAN BE USED. IF THE CURREt4T STEP SIZE
IS H AND THE VALUE AT T * E IS NEEDED, FORM
S = E/H, AND THEN COMPUTE
NQ
YCIMT + E) = SUM Y( J+l,I)»S$$J
J = 0
A BLOCK OF AT LEAST 13$N FLOATING POINT LOCATIONS.
N*3 FLOATING POINT LOCATIONS USED BY THE SUBROUTINES.
THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP
H MAY BE ADJUSTED UP OR DOWN BY THE PROGRAM
IN ORDER TO ACHEIVE AN ECONOMICAL INTEGRATION.
HOWEVER, IF THE H PROVIDED 3Y THE USER DOES
NOT CAUSE A LARGER ERROR THAN REQUESTED, IT
WILL BE USED. TO SAVE COMPUTER TIME, THE USER IS
ADVISED TO USE A FAIRLY SMALL STEP FOR THE FIRST
CALL. IT WILL BE AUTOMATICALLY INCREASED LATER.
THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE
INTEGRATION. NOTE THAT ON STARTING THIS MUST BE
MUCH SMALLER THAN THE AVERAGE H EXPECTED SINCE
A FIRST ORDER METHOD 'is USED INITIALLY.
THE MAXIMUM SIZE TO WHICH THE STEP WILL BE INCREASED
THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES
DIVIDED BY YMAX(I) MUST BE LESS THAN THIS
IN THE EUCLIDEAN NORM. THE STEP AND/OR ORDER IS
ADJUSTED TO ACHEIVE THIS.
INDICATOR. THE FOLLOWING ARE ALLOWED..
AN ADAMS PREDICTOR CORRECTOR IS USED.
A MULTI-STEP METHOD SUITABLE FOR STIFF
SYSTEMS IS USED. IT WILL ALSO WORK FOR
NON STIFF SYSTEMS. HOWEVER THE USER
MUST PRpVIOE_A SUBROUTINE PEDERV WHICH
EVALUATES THE. PARTIAL DERIVATIVES OF
THE DIFFERENTIAL EQUATIONS WITH RESPECT
TO THE Y'S. THIS IS DONE BY CALL
PEDERV(T,Y,PW,M). PW IS AN N BY N ARRAY
THE
METHOD
0
1
$00430000
$00431000
$00432000
$00433000
$00434000
$00435000
$00436000
$00437000
$00438000
+00439000
$00440000
$00441000
$00442000
$00443000
$00444000
$00445000
$00446000
$00447000
$00448000
$00449000
$00450000
$00451000
$00452000
$00453000
*00454000
$00455000
$00456000
$00457000
$00458000
$00459000
$00460000
$00461000
$00462000
$00463000
$00464000
$00465000
$00466000
$00467000
$00463000
$00469000
$00470000
$00471000
$00472000
$00473000
$00474000
$00475000
$00476000
$00477000
$00478000
$00479000
$00480000
$00481000
$00482000
-------
VIII.11
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*
YMAX
ERROR
KFLAG
JSTART AN INPUT
WHICH MUST BE SET TO THE PARTIAL OF
THE I-TH EQUATION WITH RESPECT
TO THE J DEPENDENT VARIABLE IN PW(I,J).
PW IS ACTUALLY STORED IN AN M BY M
ARRAY WHERE M IS THE VALUE OF N USED ON
THE FIRST CALL TO THIS PROGRAM.
THE SAME AS CASE 1, EXCEPT THAT THIS
SUBROUTINE COMPUTES THE PARTIAL
DERIVATIVES BY NUMERICAL DIFFERENCING
OF THE DERIVATIVES. HENCE PEDERV IS
NOT CALLED.
N LOCATIONS WHICH CONTAINS THE MAXIMUM
SEEN SO FAR. IT SHOULD NORMALLY BE SET TO
COMPONENT BEFORE THE FIRST ENTRY. (SEE THE
DESCRIPTION OF EPS.)
AN ARRAY OF N ELEMENTS WHICH CONTAINS THE ESTIMATED
ONE STEP ERROR IN EACH COMPONENT.
A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
•M THE STEP WAS SUCCESFUL,
-I THE STEP WAS TAKEN WITH H = HMIN, BUT THE
REQUESTED ERROR WAS NOT ACHIEVED.
-2 THE MAXIMUM ORDER SPECIFIED WAS FL'UND TO
BE TOO LARGE.
-3 CORRECTOR CONVERGENCE COULD NOT BF
ACHIEVED FOR H .GT. HMIN.
-4 THE REQUESTED ERROR IS SMALLER THAN CAN
BE HANDLED FOR THIS PROBLEM.
INDICATOR WITH THE FOLLOWING MEANINGS..
AN ARRAY OF
OF EACH Y
1 IN EACH
-1 REPEAT THE LAST STEP WITH A NEW H
0 PERFORM THE FIRST STEP. THE FIRST STEP
MUST BE DONE WITH THIS VALUE OF JSTART
SO THAT THE SUBROUTINE CAN INITIALIZE
ITSELF.
+1 TAKE A NEW STEP CONTINUING FROM THE LAST.
JSTART IS SET TO NO, THE CURRENT ORDER OF THE METHOD
AT EXIT. NO IS ALSO THE ORDER OF THE MAXIMUM
DERIVATIVE AVAILABLE.
MAXDER
+00483000
+00484000
+00485000
+00486000
+00487000
+00488000
+00489000
+00490000
+00491000
+00492000
+00493000
$00494000
+00495000
+00496000
+00497000
+00493000
+00499000
+00500000
+00501000
+00502000
+00503000
+00504000
+00505000
+005C600D
*00507000
+00508000
+00509000
+00510000
+00511000
4-00512000
+00513000
+00514000
+00515000
+00516000
+00517000
*0051flOOO
*00519000
*00520000
•*00521000
+00522000
*00523000
+00524000
THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN THE
METHOD. SINCE THE OROER IS EQUAL TO THE HIGHEST
DERIVATIVE USED, THIS RESTRICTS THE ORDER. IT MUST
BE LESS THAN 13 FOR ADAMS AND 7 FOR STIFF METHODS.
A BLOCK OF AT LEAST N++2 FLOATING POINT LOCATIONS.
«**************«**:**** ***:(<**«*************« 00525000
DIMENSION Y( 13,1), YMAX(l), SAVEU3.1), ERROR ( 1 ), PWT I I , 00526000
I A(13),PERTST( 12,2,3),CSAVE(N,3),IP( I) 00527000
«**********<.****** + ** + ***** + ********* ** + **00528000
THE COEFFICIENTS IN PERTST ARE USED IN SELECTING THE STEP AND +00529000
ORDER, THEREFORE ONLY ABOUT ONE PERCENT ACCURACY IS NEEDED. +00530000
***************** *****##*$********** + *******')'***** ****$**************00531000
DATA PERfST/2.,4.5, 7.333,10.42,13.7,17.15,1., !.,!.,T. ,1. »!.» ' 00532000
1 12.,12.,24..37.89,53.33,70.08,87.97,100.9,126.7,147.3,00533000
2 2168.8, 191. 4,3.,6. , 9. 167, 12 . 5, 15 . 98, 1. , l'. , 1. , I."',!., 1., 005 34 000
3 1.,312.,24.,37.89,53.33,70.08,87.97,136.9,126.7,147.3,00535000
-------
VIII.12
4 168. 9, 4191. 4, 211. ,!.,!.,. 5,. 1667, . 04 1 67, . 003333 ,
6 .0021813, .00029, 6. 000035, .0000037, ,00000035 /
DATA A(2) / -1.0 /
IRET = 1
KFLAG =1
IF (JSTART.LE.O) GO TO 140 ___
C * + **«*****<<***«**««*«*** «;***^*^*«* *«****««** ***** + ^« ************** **^
C* BEGIN BY SAVING INFORMATION FOR POSSIBLE RESTARTS AND CHANGING
C* H BY THE FACTOR R IF THE CALLER HAS CHANGED~H. ALL VARIABLES
C* DEPENDENT ON H MUST ALSO BF CHANGED.
C* E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NO. EUP ~IS
C* TO TEST FOR INCREASING THE ORDER, EDWN FOR DECREASING THE ORDER.
C* HNEW IS THE STEP SIZE THAT WAS USED ON THE LAST CALL.
C*****:******************************************************* ********.
100 00 110 I = 1,N
DO 110 J = l,K
110 SAVE( J, I) = Y(J, I )
MOLD = HNEW
IF IH.EQ.HOLO) GO T3 130
120 ?ACUtf = H/HOLD
I '. E T 1 = 1
GO TO 750
125 HOLD = H
130 SQOLD = NQ
TOLD = T
3ACUM = 1.0
IF ( JSTART.GT.O) GO TO 250
GO TO 170
140 IF ( JSTA.RT.EO.-l) GO TO 160
C* ON THE FIRST CALL, THE ORDER IS SET TO I AND THE INITIAL
C* DERIVATIVES ARE CALCULATED.
BR = 1.0
MO = I
N3 = N
ytf = N**2
CALL DIFFUN(T,Y,CSAVE( 1,3) )
DO 150 I = 1,N
150 Y(2,I ) = CSAVE(I,3)*H
HNEW = H
K = 2
30 TO 100
C* REPEAT LAST STEP BY RESTORING SAVED INFORMATION.
160 IF (NQ. EQ.NQOLD) JSTART = 1
T = TOLD
NO = NQOLD
< = NQ +• I
GO TO 120
00536000
00537000
00536000
00539000
00540000
00541000
00542000
***6~0543000
*00544000
*00545000
*00546000
*00547000
*00548000
*00549000
»**00550000
00551000
00552000
00553000
00554000
00555000
00556000
00557000
00558000
00559000
00560000
00561000
00562000
00563000
00564000
00565000
x**nnciAAnno
r *• *r \J \J J
-------
VIII.13
C* SET THE COEFFICIENTS THAT DETERMINE THE ORDER AND THE METHOD *00589000
.C* TYPE. CHECK FOR EXCESSIVE ORDER.' THE LAST TWO STATEMENTS OF " *00590000
C* THIS SECTION SET IWEVAL .GT.O IF PW IS TO BE RE-EVALUATED *00591000
C* BECAUSE OF THE ORDER CHANGE, AND THEN REPEAT THE INTEGRATION *00592000
C* STEP IF IT HAS NOT YET BEEN DONE (IRET = 1) OR SKIP TO A FINAL *0.0593000
C* SCALING BEFORE EXIT IF IT HAS BEEN COMPLETED (IRET = 2). *00594000
170 IF (MF.EQ.O) GO TO 180 00596000
IF (NQ.GT.6) GO TO 190 00597000
GO TG (221,222,223,224,225,226),NO 00598000
180 IF(NO.GT.12) GO TO 190 00599000
GO TD (205,206,207,208,209,210,211,212,213,214,215,216), NQ 00600000
190 KFLAG = -2 00601000
RETURN 00602000
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
C*
THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO THE MAXIMUM
ACCURACY PERMITTED BY THE MACHINE. THEY ARE, IN THE ORDER USED.
-1
-1/2,-1/2
-5/12,-3/4,-1/6
-3/3,-11/12,-1/3,-1/24
-251/720,-25/24,-35/72,-5/48.-1/120
-95/288,-l37/120,-5/9,-l7/96,-l/40,-l/720
-19087/60460,-49/40,-203/270,-49/I 92,-7/144,-7/1440,-1/5040
IN
METHODS TO CALCULATE FURTHER COEFFICIENTS CAN BE FOUND
HENRICI, DISCRETE VARIABLE METHODS IN ORDINARY D. E.
-1
-2/3,-1/3
-12/25,-7/10,-1/5,-1/50
-120/274,-225/274,-85/274,-15/274,-1/274
-180/441,-58/63,-15/36,-25/252,-3/252,-I/1764
205
206
207
208
209
Ad) '
GO TO
Ad ) =
A( 3) =
GO TO
Ad) =
A(3) =
A(4) --
GO TO
Ad) =
A(3) =
A(4) >
U5) -
GO TO
Ad) =
A(3) =
A(4) =
A(5) =
-1
230
" •
-0
230
-0
-0
-0
230
-0
• -o
-0
-0
230
-0
•• -1
-0
-0
.0
50000000000
.500000000
.4166666666666667
.750000000
.1666666666666667
.3750000000000
.9166666666666667
.3333333333333333
.04166666666666667
.3486111111111111
.0416666666666667
.4861111111111111
.1041666666666667
*00604000
*00605000
*00606000
*00607000
*00608000
*00609000
*00610000
*00611000
*00612000
*00613000
*00614000
*00615000
*00616000
*00617000
*00618000
*00619000
*00620000
*00621000
*00622000
00624000
00625000
00626000
00627000
00628000
00629000
00630000
00631000
00632000
00633000
00634000
00635000
00636000
00637000
00638000
00639000
00640000
00641000
-------
VIII.14
008333333333333333
3298611111111111
1416666666666667
,625000000
1770833333333333
0250000000
,001388368888888889
3155919312169312
225000000
7518518518518519
2552083333333333
04861111111111111
004861111111111111
0001984126984126984
A<6) = -0
GO TO 230"
210 A(l) = -0
A(3) = -1
A(4) - -0
A(5) = -0
A(6) = -0
A(71 = -0
30 TO 230
211 A(l ) = -0
M 3) = -1
A(4) = -0
A ( 5) = -0
A(6) = -0
A(7) = -0
A ( 8) = -0
GO TO 230
212 A(l) = -0.3042245370370370
A(3> = -1.296428571428485
A(4) = -0.8685185185185185
A(5) = -0.3357638883388883
A(6) = -0.07777777777777777
A<7) = -0.01064814814814815
A(3) = -0.0007936507936507937
A(9) = -0.00002480158730157189
30 TO 230
213 A(1I = -0.2948680004409171
A(3) = -1.358928571423505
A(4) = -0.9765542323041257
A(5) = -0.41718750000
A(6) = -0.1113541666666667
A(7) = -0.0187450000
A(8) - -0.001934523309523809
A(9) = -0.0001116071428570911
A(10) = -0.00000275573192239706
GO TO 230
214 A(l) = -0.2869754464285714
A<3) = -1.414484126933951
A(4) = -1.077215603465301
A(5) = -0.4985670194001107
A(6) = -0.148437500000
A(7) = -0.0290605709876226
A(3) = -0.0037202380952323
A<9) = -0.0002996858465601941
A(10) = -0.00001377365961195214
A(ll) = -0.00000027557319223878
GO TO 230
215 Ad = -0.230189564439367
A(3 = -1.464484126986034
A(4 = -1.171514550268255
A(5 = -0.5793581900385553
A<6 = -0.1883228615536993
A(7 = -0.04143036265485243
00642000
00643000
00644000
00645000
00646000
00647000
00643000
00649000
00650000
00651000
00652000
00653000
00654000
00655000
00656000
00657000
00658000
00659000
00660000
00661000
00662000
00663000
00664000
00665000
00666000
00667000
00663000
00669000
00670000
00671000
00672000
00673000
00674000
00675000
00676000
00677000
00678000
00679000
00680000
00681000
006fl2000
00683000
00684000
00665000
00686000
00687000
00688000
00689000
00690000
00691000
00692000
00693000
00694000
-------
VIII.15
A(8) = -0.00621114418,0002914
= -0.0006252066799684720 "
10) = -0.00004041740152967666
-0.000001515652557372595
-.00000000250521083864932
274265540031599
509938672424507
260271163994607
659234182076101
2304580026342793
05569724610169335
009439484126243074
-0.001119274966826774
00009093915342935554
000004822530863610993
0000001503126502923833
000000002087675698479380
-0.6666666666666667
-0.3333333333333333
30
= - 0.5454545454545455
,09090909090909091
,480000000
,700000000
200000000
020000000
437956204379562
8211678^32116788
,3102189781021898
05474452554744526
0036496350364963504
4081632653061225
,9206349206349206
,4166666666666667
0992063492063492
,0119047619047619
000566893424036282
216
221
222
223
224
225
226
230
A(9) = -0,
A(10) = -(
A(ll) - -(
A(12) = -,
GO TO 230
A(l) =.-0,
A(3) = -I,
A(4) = -1,
A (5) = -0,
A (6) = -0.
A(7) = -0,
A(8) = -0,
A (9) = -0,
A<10) = -(
A( 11) = -(
A( 12) - -(
A( 13) = -(
GO TO 230
A ( 1 ) = - 1 ,
GO TO 230
All) = -0,
A(3) = -0,
GO TO 230
Ad) = - (
A ( 3 ) = AC
A (4) = -0,
GO TO 230
A( 1 ) = -0,
A(3) = -0,
A(4 ) = -0,
A(5) = -0,
30 TO 230
A( 1) = -0,
A(3) = -0,
A(4) = -0,
A(5) = -0,
A(6> = -0.
GO TO 230
A(l) = -0,
A (3) = -0,
A(4) - -0,
A (5) = -3,
AI6) = -0,
A<7) = -0,
< = NQ+1
IDOUB = K
MTYP * (4
ENQ2 = ,5i
ENQ3 = .5,
ENQ1 = 0. !
PEPSH = El
- MF1/2
5/FLOATINO *
5/FLOAT(NQ +
0.5/FLOATINQ)
I)
2)
EUP
PERTST(NQ,MTYP,2)*PEPSH)
00695000
00696000
00697000
00693000
00699000
00700000
00701000
00702000
00703000
00704000
00705000
00706000
00707000
00708000
00709000
00710000
00711000
00712000
00713000
'00714000
00715000
00716000
00717000
00713000
00719000
00720000
00721000
00722000
00723000
00724000
00725000
00726000
00727000
0072SOOO
00729000
00730000
00731000
"00732000
00733000
00734000
00735000
00736000
00737000
00738000
00739000
00740000
00741000
00742000
00743000
"00744000
00745000
00746000
00747000
-------
VIII.16
E » MTYP,3T*PEPSH)**2 00749000
IP (EDWN.EQ.O) GO TO 780 00750000
8ND = **2 00751000
240 IWEVAL = MF 00752000
GO TO { 250 , 680 ).IRET 00753000
c* THIS SECTION'COMPUTES THE "PREDICTED VALUES BY EFFECTIVELY $00755000
C* MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE *00756000
C* MATRIX. *00757000
250 T = T + H 00759000
00 260 J * 2»K 00760000
DO 260 Jl =~J»K " "" "" " 00761000
J2=K-Jl+J-l 00762000
DO 260 I = 1*N 00763000
260 Y(J2»I) = YU2.I) + Y(J2M,M 00764000
C***********************************************************************00765000
C* J? TO 2 CORRECTOR ITERATIONS ARE TAKEN._CONVERGENCE IS TESTED BY *00766000
C* REQUIRING THE L2 NORM OF CHANGES TO BE LESS THAN BND WHICH IS *00767000
C*' DEPENDENT ON THE ERROR TEST CONSTANT. *00768000
C* THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE ARRAY *00769000
C* ERRORU). IT IS EQUAL TO THE K-TH DERIVATIVE OF Y MULTIPLIED *00770000
C* BY H**K/(FACTORIAL(K-1)*A(K)). AND IS THEREFORE PROPORTIONAL *00771000
C* TO THE ACTUAL .ERRORS TO THE LOWEST _POWER_ OF H PRESENT. <.H**K) *00772009
00 270 I = 1,N 00774000
270 ERRORU) = 0.0 00775000
DO 430 L=l»2 00776000
CALL DIFFUN (T,1,CSAVE(1»3>) 00777000
C***********************************************************************g0778000
C* IF THERE HAS BEEN A CHANGE OF" ORDER OR THERE HAS BEEN TROUBLE *00779000
C* WITH CONVERGENCE* PW IS RE-EVALUATED PRIOR TO STARTING THE *00780000
C* CORRECTOR ITERATION IN THE CASE OF STIFF METHODS. IrtEVAL IS *00781000
C* THEN SET TO -1 AS AN INDICATOR THAT IT HAS BEEN DONE. *00782000
C***********************************************************************00783000
IF (IWEVAL.LT.1) GO TO 350 0078400J
IF (MF.EQ.2) GO TO 310 " 00785000
CALL PEDERV(T,Y,PW,N3I 00786003
R = A(l)*H 00787000
DO 230 I = 1.N4 00788000
280 PW(I) = PW(I)*R 00789000
C **«*#****!•<*************************************************************Q 0790000
C* ADD THE IDENTITY MATRIX TO THE JACOBIAN~AND DECOMPOSE INTOlu * PW *00791000
C***********************************************************************Q0792000
290 DO 300 I = 1»N 00793000
300 PWt I*(N3 + 1)-N3) = 1.0 * PW( I *(N3+1 )-N3) 00794000
IWEVAL = -1 00795000
CALL DECOMP(N»N3*PW>IP) 00796000
IF (IP(N) .NE. 0) GO TO 350 ' 00797000
GO TO 440 00798000
C************************************************************#**********00799000
C* EVALUATE THE JACOBIAN INTO PW BY NUMERICAL DIFFERENCING. R IS THE *00800000
-------
VIII. 17
C* CHANGE MADE TO THE ELEMENT OF Y. IT IS EPS RELATIVE TO Y WITH *00801000
C* A MINIMUM OF EPS**2. F STORES THE UNCHANGED VALUE OF Y. *00802000
310 DO 340 J = l,N 00804000
F = Y(1»J) 00805000
P. = EPS*DMAX1S»DABS(F) ) 00806000
Y(l.J) = Y(l.J) + R 00807000
D = Ad)*H/R OOSOSOOO
CALL DIFFUN(T, Y,CSAVEd,l ) ) 00809000
DO 330 I = 1»N 00810000
330 PW( I+< J-1)*N3) = (CSAVEdhl) -CSAVEd,3» * D 00811000
340 Y(1,J) = F 00812000
GO TO 290 00313000
350 DO 360 I = UN ' 00814000
360 CSAVEtl.l) = Y(2,I) - CSAVE(I.3)+H 00815000
IF(MF.EQ.O) GO TO 410 00816000
CALL SQLVE(N,N3,PW.CSAVE( 1, 1),IP) 00817000
C**************************************************x--********************00813000
C* CORRECT AND COMPARE DEL* THE L2 NORM OF CHANGE/YMAX, WITH BND. *00319000
C* ESTIMATE THE VALUE OF THE L2 NORM OF THE NEXT CORRECTION BY *00320000
C* BR*2*DEL AND COMPARE WITH BNO. IF EITHER IS LESS, THE CORRECTOR *00821000
C* IS SAID TO HAVE CONVERGED. *00822000
Q *********************************************************** ************OQ823000
410 DEL = 0.000 00824000
DO 420 I = 1..N 00825000
Y(l,I) = Yd, I) +• Ad )*CSAVF ( I ,1 ) 00826000
Y(2>I) = Y(2,I) - CSAVE(I.l) 00827000
ERROR(I) = ERRORU) + CSA«/E(I,1) 00828000
DEL = DEL * (CSAVEl I, 1) /YMAXI I) l**2 00829000
420 CONTINUE 00830000
426 IFCL.GE.2) BR = DMA XI ( . 9*BR, DEL/DELI) 00831000
DELI = DEL " 00832000
IF(OMINl(DEL, BR*DEL*2.0) .LE. BND) GO TO 490 00333000
430 CONTINUE 00534000
Q**«******* **************************************** ************Q0842000
440 T = TOLD 00843000
IF { (H.LE. «HMIN*l. 00001 )). AND. (< IWEVAL - MT Y° ) . L T . - 1 ) ) GO TO 460 008^4000
IF ( (MF.EQ.O) .OR. ( IWEVAL. NE.O) ) RACUM = R AC UM*0. 25DO 005^.5000
IWEVAL = MF 00846000
IRET1 = 2 00847000
30 TO 750 00843000
460 KFLAG = -3 00349000
MQ = NQOLD 00850000
470 DO 480 I = l,N 00851000
DO 480 J = 1,K 00852000
490 Y(J,I) = SAVE(J,I) 00853000
-------
VIII. 18
H = HOLD 00354000
JSTART = NO 00855000
RETURN . 00856000
********** *************************************** *+*********00857000
C* THE CORRECTOR CONVERGED AND CONTROL IS PASSED TO STATEMENT 520 *00858000
C* IF THE FRROR TEST IS O.K.. AND TO 5*0 OTHERWISE. *00859000
C* IF THE STEP IS O.K. IT IS ACCEPTED, IF IDOUB HAS SEEN REDUCED *00860000
C* TO ONE, A TEST IS MADE TO SEE IF THE STEP CAN BE INCREASED *00861000
C* AT THE CURRENT ORDER OR 3Y GOING TO ONE HIGHER OR ONE LOWER. *00862000
C* SUCH A CHANGE IS ONLY MADE IF THE STEP CAN BE INCREASED BY AT *00863000
C* LEAST 1.19 IF NO CHANGE IS POSSIBLE IDQUB IS SET TO 3 TO *00864000
C* PREVENT FUTHER TESTING FOR 6 STEPS. *00865000
C* IF A CHANGE IS POSSIBLE, IT IS MADE AND IOOUB IS SET TO *00866000
C* NQ * 1 TO PREVENT FURTHER TESING FOR THAT NUMRER OF STEPS. *00867000
C* IF THE ERROR WAS TOG LARGE, THE OPTIMUM STEP SIZE FOR THIS OR *00868000
C* LOWER ORDER IS COMPUTED, AND THE STEP RETRIED. IF IT SHOULD *00869000
C* FAIL TWICE MORE IT IS AN INDICATION THAT THE ' DER I VAT I VES THAT *00870000
C* HAVE ACCUMULATED IN THE Y ARO.AY HAVE ERRORS OF THE WRONG ORDER *0087LOOO
C* SO THE FIRST DERIVATIVES ARE RECOMPUTED AND THE CRDER IS SEJ *00872000
C* TO 1. *00873000
D = 0.0 00875000
DO 500 I = l,N ; 00876000
500 D = D + (ERRORU )/YMAX( I I )**2 00877000
IWEVAL = 0 00878000
IF (D.GT.E) GO TO 540 " 00879000
IF (K.LT.3) GO TO 520 00880000
C*#*#*$$**$*******$******************************** ******** *************00861000
C* COMPLETE THE CORRECTION OF THE HIGHER ORDER DERIVATIVES AFTER A *00882000
C* SUCCESFUL STEP. *00883000
DO 510 J = 3,K " ' 00885000
DO 510 I = 1»N 0088600.")
510 Y(J,I> * Y(J.I) + A( J)*ERROR( I) 00887000
520 KFLAG = +1 00888000
HNEW = H 00889000
IF IIDaU8.LE.ll GO TO 550 00890000
IDOUB = IOOUB - I 00891000
IF (IDOUB. GT.l) GO TO 700 00892000
DO 530 I = 1,N 00893000
530 CSAVE(I,2) = ERPOR(I) 00894000
30 TO 700 00895000
C* REDUCE THE FAILURE FLAG COUNT fo CHECK FOR MULTIPLE FAILURE'S. *00897000
C* RESTORE T TO ITS ORIGINAL VALUE AND TRY AGAIN UNLESS THERE HAVE *00898000
C* THREE FAILURES. IN THAT CASE THE DERIVATIVES ARE ASSUMED TO HAVE *00899000
C* ACCUMULATED ERRORS SO A RESTART FROM THE CURRENT VALUES OF Y IS *00900000
C* TRIED. THIS IS CONTINUED UNTIL SUCCESS OR H = HMIN. *00901000
540 KFLAG = KFLAG - 2 ~"~ 00903000
IF (H.LE. (HMIN*l. 00001) ) GO TO 740 00904000
T = TOLD 00905000
IF (KFLAG. LE. -5) GO TO 720 00906000
-------
VIII.19
C* PRl. PR2, AND PR3 WILL CONTAIN THE AMOUNTS 3Y WHICH THE STEP SIZE
C* CHOULD BE DIVIDED AT ORDER ONE LOWER. AT THIS ORDER* AND AT ORDER
C* ONE HIGHER RESPECTIVELY.
550
PR2 = tD/E)**ENQ2*l. 2
PR3 » l.E+20
IP ((NQ.GE.MAXDER).OR. (KFLAG.LE.-l) )" GO TO 570
D =* 0.0
DO 560 I = l.N
0 = D + ((ERROP(I) - CS AVE ( I , 2 ) ) / YMAX ( I) ) **2
PR3 = (D/EUP)**ENQ3*1.4
PRl = l.E+20
IF (NQ.LE.l \ GO TO 590
D = 0.0
DO 580 I = 1»N
0=0* ( Y(K, I) /YMAXC I ) )**2
PRl = CO/EDWN)**ENQ1*1.3
CONTINUE
IF (PR2.LE.PR3) GO TO 650
IF (PR3.LT.PR1) GO TO 660
R = 1.0/AMAXKPRl.l. E-4)
\4EWQ = NO - 1
610 IDOUB = 8
IF ((KFLAG.EQ. II. AND. (R.LJ. (1. 1 ) ) ) GO TO 700
IF (NEWQ.LE.NO) GO TO 630'
560
570
580
590
600
C* COMPUTE ONE ADDITIONAL SCALED DERIVATIVE IF ORDER IS INCREASED.
DO 620 I * l.N
620 Y(NEWO+UI) = ERROR (I ) * A(K)/DFLOAT ( K)
630 K = NEWQ "V 1
IF (KFLAG.EQ.1) GO TO 670
RACUM = RACUM*R
IRET1 = 3
GO TC 750
640 IF (NEWQ.EQ.NQ) GO TO 250
NO = NEWQ
GO TC 170
650 IF (PR2.GT.PRl) GO TO 600
NEWQ = NO
* = 1.0/AMAXKPR2* l.E-4)
GO TO 610
660 R * l.O/AMAXl(PR3.1. E-41
NEWQ » NQ * 1
GO TO 610
670 IRET = 2
R = DMINl(R»H«AX/DABS(HM
= H "
IF (NQ.EQ.NEWQ) GO TO 680
NQ = NEWQ
GO TD 170
*00908000
$00909000
*00910000
11000
00912000
00913000
00914000
00915000
00916000
00917000
00918000
00919000
00920000
00921000
00922000
00923000
00924000
00925000
00926000
00927000
OOQ28000
00929000
00930000
00931000
00932000
*00934000
00936000
00937000
00938000
00939000
00940000
00941000
00942000
00943000
00944000
00945000
00946000
009^7000
00948000
00949000
00950000
00951000
00952000
00953000
00954000
00955000
00956000
00957000
00958000
00959000
-------
VIII. 20
680 Rl = 1.0
00 690 J = 2,K
RI •= RI*R
00 690 I « 1»N
690 Y(J»I) « Y(J,I)*R1
IDOUB = K .
700 DO 710 I = 1»N
710 YMAX(I) = DPAXl(YttAXm,O.ABS( Yd, I ) ) )*
JSTART * NO
RETURN
720 IF(NQ.GT.l) GO TO 725
IF (H .LE. 2.DO*HMIN) GO TO 780
GO TO 550
725 DO 730 1= UN
Yd, I ) = SAVEd,! )
SAVE(2»I> = HOLD*CSAVE< 1.3)
730 Y(2,I) * SAVE12»I)*R
R = H/HOLD
CALL OIFFUN ( T ,Y »CS A VE ( I, 3 ) )
NQ * I
KFLAG = I
GO TO 170
740 KFLAG = -I
HNEW = H
JSTART = NQ
RETURN
00960000
00961000
00962000
00963000
00964000
00965000
00966000
00967000
00968000
00969000
00970000
00971000
00972000
00973000
00974000
00975000
00976000
00977000
00976000
00979000
00980000
00981000
00982000
00983000
00984000
00985000
C* THIS SECTION SCALES ALL VARIABLES CONNECTED WITH H AND RETURNS *00987000
C* TO THE ENTERING SECTION. *00988000
C* *********************************************************** ***********Q0969000
750
760
770
760
RACUM = OMAX1(DABS(HMIN/HOLO),RACUM)
RACUM = DMIN1 '(RACUM,DABS(HMAX/HOLD) )
Rl = l.o
00 760 J » 2,K
Rl = Rl*RACUM
DO 760 I » 1>N
Y(J,I) = SAVE(J,I)*R1
H = HOLD*RACUM
00 770 I = 1,N
Yd, I) = SAVEd»I)
IDOUB
30 TO
KFLAG
GO TO
END
* K
( 125
» -4
470"
250
640 l.IRETl
00990000
00991000
00992000
00993000
0099AOOO
00995000
00996000
00997000
0099SOOO
00999000
01000000
01001000
01002000
01003000
01004000
-------
VIII.21
1AY 71)
OS/360 FORTRAN H
1PILER OPTIONS - NAME = MAIN,OPT=01.LINECNT=55,SIZE=OOOOK,
SQURCE,BCO,NOLIST,NOOEClOLOAO/MAP,NOEDIT,ID>NaXREF
SUBROUTINE DECOMP(N»NO IM,A. IP)
IMPLICIT REAL*8(B-H,Q-Z)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A IS SINGLE
PRECISION.)
MATRIX TRIANGULARIZATIQN BY GAUSSIAN ELIMINATION.
INPUT...
N = ORDER OF MATRIX
NDIM = DECLARED DIMENSION OF ARRAY A.
H = MATRIX TO BE TRIANGULARI ZED. (FOR STIFF METHODS*
PRECISION; ALL OTHER VARIABLE ARE DOUBLE
OUTPUT...
AU»J), "l.LE.J = UPPER TRIANGULAR FACTOR. U.
A(I.J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I-L
IP(K)* K.LT.N = INDEX OF K-TH PIVOT ROW.
IP(N) = (-1)**(NUMBE<* OF INTERCHANGES) OR 0.
USE 'SOLVE1 TO OBTAIN SOLUTION OF LINEAR SYSTEM.
DETERM(A) = IP(N ) *A(1,1)* A(2,2)*...*A(N,N).
IF IP(N) « 0, A IS SINGULAR, SOLVE WILL DIVIDE BY ZERO.
DIMENSION A(NDIM,NOIM(»IP(NDIM)
DIMENSION B(4,4)
COMMON NFNS»NW,B
Nw = NW * 1
IP(N) = 1
00 6 K=1,N
IF(K.EQ.N) GO TO 5
KP1 = K+l
1 = K
DO 1 I=KP1»N
IF(ABS(A(I,K) ).GT.ABS(A(M,K)1 I M=I
I CONTINUE
IP IK) * M
IF(M.NE.K) IP(N) = -IP(N)
T = A(M,K)
A(M*K) = A(K,K)
A(K,K) = T
IF(T.EQ.O) GO TO 5
DO 2 I=KP1»N
2 A( I,K) = -A(I.K) / T
DO 4 J=KP1,N
T = A(,1,J)
A(M»J)=A(K,J)
A(K,J) = T
IF(T.EQ.O.) GO TO 4
DO 3 I=KP1,N
3 A(I»J) = A( I,J) * A(I,K)*T
4 CONTINUE
5 IF(A(K,KJ.EQ.O.) IP«N) = 0
6 CONTINUE
RETURN
END
01005000
01006000
01007000
01008000
01009000
01010000
01011000
01012000
01013000
01014000
01015000
01016000
01017000
OlOlflOOO
01019000
01020000
01021000
01022000
01023000
01024000
01025000
01026000
01027000
01028000
01029000
01030000
01031000
01032000
01033000
0103AOOO
01035000
01036000
01037000
01038000
01039000
01040000
01041000
01042000
01043000
0 1044000
01045000
010^6000
01047000
01048000
01049000
01050000
01051000
01052000
01053000
01054000
01055000
-------
VIII.22
05/360 FORTRAN H
COMPILER OPTIONS - NAME = MAI N,OPT=01,LINECNT=55»SIZE = OOOOK»
SaURCE»8CO»NOLISt»NOOECK,LOAD»MAP^NOEDIT»lO^NOXREF
SUBROUTINE SQLVE(N»NO IM»A.8»IP)
IMPLICIT REAL*fl <8-H,Q-Z)
C
c
C
c
c
c
c
c
c
c
c
SOLUTION OF LINEAR SYSTEM/ A*X = B.
INPUT...
N = ORDER OF MATRIX.
NDIM = DECLARED DIMENSION OF ARRAY A.
A = TRIANGULARIZED MATRIX OBTAINED FROM •
3 = RIGHT HAND SIDE VECTOR.
IP = PIVOT VECTOR OBTAINED FROM 'DECOMP'.
OUTPUT...
3 = SOLUTION VECTOR, X.
DIMENSION A(NDIM.NDIM), B(NDIM), IP(NDIM)
IF(N.EQ.1» GO TO 9
NM1 = N-l
DO 7 K*1,NM1
KP1 = K + 1
M = IP(K)
T = B(M)
B(M) = B(K)
3
-------
VIII.23
OS/360 FORTRAN H
COMPILER OPTIONS - NAME* MA I N, OPT=01»IINECNT*55»SIZE=OOOOK,
SOURCE.BCD.NOLIST.,NODECK»IOAD,MAP.NOEDIT, ID.NOXREF
SUBROUTINE PEDERV (T.r*PW»-N3) 01089300
DIMENSION Y(13»40).PW(1600) 01089320
DOUBLE PRECISION T.Y 01089340
RETURN 01089360
END 01089380
-------
VIII.24
OS/360 FORTRAN H
COMPILER OPTIONS - NAME = MA I N, a<>T=0 1* L I NECNT= 55» S I Z£ =OOOOK,
SOURCE»BCQ»NOLIST,NOD?CK.LaAD»MAP»NOEDIT»ID»NOXREF
SUBROUTINE SMGOUT
-------
VIII.25
c
c
c
c
c
30
40
70
71
73
80
61
63
90
91
93
110
120
130
140
145
SAVE PEAK
IF (Yd ).
PEAK1- Y(
TIME1= X
IF (Y(3).
PEAK3= Y(
TIME3
IF (Y
IF (Y
T IM1 =
PREVY
PTIM1
IF (Y
IF (Y
TI12 =
PREVY
PTIM?
IF (Y
IF (Y
TIM4 =
-
d
(1
PT
1 =
= X
(2
(2
PT
2=
= X
(4
(4
PT
OREVY4=
PTIM4
CHECK
IF ( I
IF (X
WRITE
IF (I
WRITE
IF (I
WRITE
WRITE
IF (I
WRITE
1
WRITE
WRITE
WRITE
WRITE
= X
I
X
)-
) .
VALUES, P
LT.PEAK1)
1)
LT.PEAK3)
3)
PREVY1) 71
GT.HAL1.0R
IMI-(PTIM1-X
Y(
) -
) .
1)
PREVY2) 91
GT.HAL2.0R
IM2-IPTIM2-X
Y(
)-
).
2)
PREVY4) 91
GT.HAL4.QR
IM4-(PTIM4-X
Y(
F
NTR)
•
(
TR
(
TR
(
(
LT
6,
AF
6.
AF
6,
6,
SOUR
(
6,
4)
OUTPUT THI
130, 130.
. (STEP(fl)
EAK TIMES.
GO TO 40
GO TO 70
, 71
.PRE
)*((
. 3.1
.PRE
)*( (
, 91
.PRE
)*( (
S TI
120
-0.5
, 60
VYl.LT
PREVYl
, 90
VY2.LT
_ AND_ HALF-TIMES.
.MALI) GO TO 73
-HALl)MPREVYl-Y(l) ) )
.HAL2) GO TO 83
PREVY2-HAL2)/(PREVY2-Y(2)»
,110
VY4.LT
PREVY4
ME.
) ) GO
1000) X»NN, XPOS, YPQS
.EQ. 0 .AND.
1012)
.EQ. 0)
1013) (POL
ISOUR
.HAL4) GO TO 93
-HAL4)/(PREVY4-Y(4) ) )
TO 170
,H,M
.EQ. 0) GO TO 145
GO TO 140
L( I,
1 ) . I = 1
1007) (POLLt I, 2), 1=1
.EQ. 0)
1011) POLL
,4), (GPOLU I.I). 1 = 1,4)
,4), (GPOLL( 1,2) , 1=1,4)
GO TO 145
(1.3
,GPOLL (4,
(
(
(
(
DO 150
150
WRITE
CONTI
(
6.
6,
6.
6.
J
6,
no9 )
1010) TEMP
1001 )
1002) (I,L
=NNDIF,NC
1003) J.L
)»POLL
3)
. HUMID. HT,
ABEL
ABEL
(.1 ) » Y(
( J),YY
(2, 3), POLL (4, 3), GPOLU lr 3)
XHT,XDEPTH
I)*OERY( I) . I = 1,NDIF)
(J )
NUE
IF (NEQ .
160
170
WRITE
WRITE
M = M
IF (B
(
(
+ 1
B
6,
6,
EQ. 1)
1005)
1006) ( I. (
.EO. l.OE+08
GO T
RCT(
0 160
I, J) . J
» GO TO I
= i,3),0( I ), 1=1, NEQ)
30
NN=NN+l
C
C
180
88 =
CHECK
0.
0
FOR
NEGATIVE
CONC
ENTRAT
IONS OF SPECIES.
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,GPOLL(2.3)0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
1
I
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
138000
139000
140000
141000
142000
143000
144000
145000
146000
1^7000
148000
149000
150000
151000
152000
153000
154000
155000
1156000
1
1
1
1
1
1
I
I
1
I
1
I
1
1
1
I
1
I
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
157000
158000
159000
160000
161000
162000
163000
164000
165000
166000
166200
166400
167000
168000
168200
169000
170000
170100
171000
172000
173000
174000
175000
176000
177000
177200
173000
179000
130000
181000
182000
183000
154000
165000
-------
VIII.26
400 00 405 I*1»NDIF
IF (Yd)) 401*405.405
401 Y( I ) = 0.
405 CONTINUE . .
DO 415 I=NNDIF.NC
IF mm.) 411,415.415. . _
411 YY(I) » 6.
415 CONTINUE
3B=BB>1.0
RETURN
1000 FORMAT('1','TIME=',015.8,• MINJTES ' • 5X»'STEP='•\ 3*5X»•POS ITION=
1F8.4,.','.F8.4.6X,'STEP S I.ZE» • »DX§. a» 6X_» • P AGE ' »I 4,///)
1001 FORMATl//' COMPOUND CQNC PPHM
1002 = ORMAT(' ',I3.1X.A8.2016.7)
1003 FORMAT!' •,I 3»IX»A8,016.7)
1005 PORMAT!«i',16X,'RATES OF REACTION'/)
1006 FQRMAT( I 4,4X,3A8»Dl6.5 )
1007 FORMAT!'OFREEWAY TR AFF I C ' . 4X> 40 1 2. 4 . 4X~, 4E12~. 4 )
CHANGE. PPHM/MINM
1008 FORMAT CODIFFUSION EMISSIONS
I 6X,'NO='.E15.8)
1009 FORMAT(//1HO»16X,'RELATIVE'
16X,'MIXING',/ ' TEMPERATURE
IN GRAMS/MIN'.6X,
HUMIDITY1
,E15. 8.
.18X.'MIXING'.
2L E VAT I ON',6 X, 'HEIGHT' ,7X, 'DEPTH «,/6X, MCE NT. ) '.
3 ' • (METERS) (METERS) (METERS)')
1010 FORMAT( 1H0.5F 12.4)
1011 FOR MAT( 'OAREA SOURCES•.7X»2012.4*12X*D12.4,4X»2E12.4*12X*E12.4/
1012 FORMAT(30X,'EXPECTED RATE OF CHANGE OF CONCENTS! AT IONS DUE TO EM
1IONS FROM REGIONAL SOURCES'//40X,'PPHM/MIN'.42X.'GRAM-MOLES/MIN
222X,'PROPYLENE'.7X, ' NG' . 10X» ' CD « . 8X, '.DUMHC ' . 9X , ' PROP YLENE ' . 7X . •
3, 10X. 'CO',8X, 'DUMHC ' )
1013 FORMAT!'OSTREET TRAFFIC'.5X»4012.4.4X,4E12.4)
END
01186000
01187000
01188000
01189000
01190000
01191000
01192000
Ofl93000
01194000
01195000
01196000
01197000
01198000
01199000
01200000
01201000
01202000
01203000
01204000
01205000
01206000
01207000
01208000
01209000
01210000
01211000
E01212000
01213000
01214000
01215000
// )01216000
ISS01216110
'//01216120
N0'01216130
01216140
01217000
01218000
-------
VIII.27
OS/360 FORTRAN H
COMPILER OPTJONS - NAME = MAIN . OPT = 01.»L I NcCN.I»_55.»S I2E»OOOOK*
SOURCE»BCO.NOL"lST»N"obECK.LbAD»MAP»NaEDIT, IO,NOXREF
SUBROUTINE DIFFUN (X,YSAVE,OY) 01219000
DIMENSION YSAVE(13. 40),Y<40)»DY(40) 01220000
DOUBLE PRECISION YSAVE,Y,OY,C»0»Z,X»YY,ZZ,3B»DERY,POLL.PRMT,XPREV 01221000
REAL*8 LABEL,RCT " " 01222000
COMMON /KINCOM/ ZZ ( 4 ) , Y Y ( 40 ) , L A BELJ 40) ? D ( 50 ) ,RCT ( 50, 3 ) , PRMT ( 8) , 01223000
"l OERY«40)»C(50")j3B»Z(2)»"N»NC*NEQ,NOIF»KDIFfMlNN. INTR»NNC. " 01224000
2 NNDIF,IUV,IM» IS.ITRAF,NLL»CC(16,100),A(80)VSTEP(15)»TIM1. 01225000
3 TIM2»TIM4,PEAKl.T~IMElF>EA'K3»TIME3*XPOS,YPOS' * 01226000
C 01227000
COMMON /ALLCOP/ POLL ( 4» 3 ) » GPOLL < 4, 3 ) » ANS ( 3 ) » GOT I ME > XDEPTH* I SOUR • 01228000
1 XOEPRV _ 01228200
DO 2 1 = 1,NC " "~ " 01229000
2 Yd) = YSAVEdrl) 01230000
IF (NEQ.EO. 1) GO TO 18 01230200
C 01231000
C COMPUTE NON-ACCUMULAING SPECIES. 01232000
C 01233000
IF (Y(2M 6»6»8 01234000
6 Y(2) = (0(1) )MC(3)*Y(3)+C(6)*Y<24) +C(8)*Y(16)+ 01235000
1C(16)*Y(19) * C(17)*Y(18) * C(18)*Y(25) + C(l9)*Y(2l) * 01236000
2C(20)*YC22)) 01237000
DERY(2) * 0. 01238000
C _ _ _ 01238500
8 YY(14)"« (0(14)*100".00*2.00*0(32 )>/(C(9) + ZZ(3») 01239000
C 01240000
YY(15) * (0(8»*0(11)+D(25) )/ 012*1000
I (C(7)*YU)+C(12)*Y(l)+C(13)*Y(3)«-C(l'+ZZ(3)) 01253000
C 01254000
YY(21) = (0(5)+D(18) )/(C(19)*Y(2)+ZZ(3)) 01255000
C 01256000
YY(22) = (0(23)*D(24) >/"(C ( 20 ) *Y ( 2 ) + C ( 25 ) ) " 01257000
C 01253000
YY(23) = 0(1)/(C(2)*C(22)*Y(4)*C('24)*Y(4)*C(31)*Y(12) ) 01259000
C 01260000
YY(24) = 0(4)/(C(6)*Y(2)+C(33)*Y(1) ) 01261000
C _ 01262000
YY(25) «= (D( 10)+0(25) )/ (C( 18) *Y( 2")*C ('30 )*Y ( 1 ) *C ( 5 ) *ZZ ( 3 ) ) " 01263000
C 01264000
C TEST FOR NEGATIVE S3LID STATE SPECIES 01265000
C 01266000
-------
VIII.28
DO 10 I=NNDIF,NC
IF (YY(I) .IT. 0.0)
10 :ONTINUE
YY( I )
0.0
c
c
c
COMPUTE CHAIN YIELD.
ZZ(2) - YY(14)+YY(15)+YY( 16)+YY< 17)+YY(18)+YY(19)+YY(20)+YY(21>+
1 YY(25)
ZZ(3) * ZZ(2) *Z( 1)*Y(2)*Z(2)
SET TRAJECTORY POSITION
13 XPOS = Y(NOIF+1)
YPOS = Y(ND~IF + 2)
) ,OER Y ( NO IF+2 ) )
CALL METEOR ( XPOS, YPOS» 0* X. XPRE V/DERY (NDI
IF (NEQ .EQ. 1) GO TO 40
IF ( IUV .NE. 0) GO TO 20
CALL ULTRAV (C<1»)
20 IUV = IUV-H
IF ( IUV .GT. 100) IUV=0
COMPUTE RATES OF REACTION.
0(1 ) = C( 1)*Y( I)
0( 2) = C( 2)*YY(23)
IF (Y(?) .LE. 0.0)
1 Y(2) = (0( 1) )/(C(3)*Y< 3)+C(6)*YY(24)+Cm*YY(16)+C( 16)*YY( 19) +
2 C(17)*YY(13)+C(18)*YY(25)*C(19)*YY(21)+CJ26>*YY(22))
3
D
D
0
D
0
D
D
0
0
0
D
D
0
0
0
0
0
0
D
D
D
( 3)
( 4)
( 5)
( 6)
( 7)
( 8)
( 9)
(10)
(11)
(12)
(13)
(14)
(15)
( 16)
( 17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
_ f*
= C
= C
= C
= c
= c
= c
= c
= 0
= c
= c
= c
= c
= c
= c
= c
= c
= c
» c
= c
= c
= c
( 3)
( 4)
( 5)
( 6)
( 7)
( 8)
( 9)
(10)
.00
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
*Y( 3)
*Y( 1)
*YY(25
*YY
*YY
*YY
*YY
*YY
*YY
*YY
*YY
*YY
*YY
*YY
24
15
16
14
20
15
15
15
17
19
18
*YY(25
*YY(21
*YY(22
*YY(18
*YY(23
* Y( 4
*YY(23
*Y(2
*Y(3
)
) * Y (
) *Y(
)*Y(
)
)
)*Y(
)*Y(
)*Y(
)
) *Y (
)*Y(
)*Y(
)*Y(
)*Y(
)
)*Y(
)*Y(
)*Y(
)
)
2
4
2
1
)
)
)
)
3)
8
2
2
2
2
2
4
3
4
)
)
)
)
)
)
)
)
)
01267000
01268000
01269000
01270000
01271000
01272000
01273000
01273100
01274000
01275000
01276000
01277000
01273000
01279000
01279500
01260000
012»0200
01280400
01281000
01282000
01283000
01284000
01285000
01266000
01287000
01288000
01289000
0 1290000
01291000
01292000
01293000
01294000
01295000
01296000
0129700')
01293000
0129900-)
01300000
01301000
01302000
01303000
01304000
01305000
01306000
01307000
01303000
01309000
01310000
01311000
01312000
01313000
01314000
01315000
-------
VIII.29
0(25) = C(25)*YY(22)
0(26) = C(26)*YY(16)*Y(4)
0(27) = C(27)*YY(19)4'Y(4)
D(28) = C(28)*YY(20)
0(29) = C(29)*YY(18)*Y(1)
0(30) = C(30)*YY(25)*Y(1)
0(31) = CJ31)*YY(23^*Y( 12)
0(32) = C(32)* Y( 7)
0(33) * C(33)*YY(24)*Y(II
COMPUTE RATES OF CHANGE.
DERY( 1)=+0( 3)+0( 6)+D( 6)
1+D(19)+D(20)-D( l)-0( 4)
DERY( 2)=+D( 1) -D( 3)-0(
1-D(18)-D(19)-0<20)
OERY(3) = 0(2 )-D(3)-D(4)+D(5)-D(13)-D(23)
DERY( 4)=-D(22)-D(23)-D(24)-D(26)-D(27)-D(7)
DERY( 5)= + D(12) + 0(33) + D(33)
DERY( 6)=*D(20)+D(26) *• D(7)
OERY(7)=+0(21)+D(23)+D(24)-D(32)
DERY (8) = -0(14) * (D(28)-t-D(32) )/100.DO
DERY( 9)=+D(29)
DERY(10)=+0(30)
OERY(ll) = 0.00
*D( 8) *0(16)+D(17)+D(18)
-D(12)-D(29)-D(30) -0(33)
6) -0( 8) -D(16)-D(17)
OERY( 12) = -0(31)
'.O
55
60
65
C
C
C
70
80
131
200
IF (
IF (
IF (
CALL
I SOUR
XPOS
YPOS
•
.L
.L
SOURCE
IF (NEQ .
OESY(2) =
OERY(4) =
EQ.
E.
E.
0
18.
2.
(XPOS
EQ. 1
OERY
)
(2)
.AND.
.OR.
.OR.
»YPO.S)
ITPAF
XPOS .
YPOS .
GO TO 60
* PGLLJ2,
+ POLL( 1.
I)
1)
.EQ.
GE.
GE.
0
66
52
+PQLL(
+POLLI
)
«
*
2
I
GO TO 65
) GO TO 65
) GO TO 65
,2)+PQLL(2»3
»2)+POLL(l»3
(8) = DERY(8) > (POLL13.D+POLL (3.2) )/100.DO
DERY(12) = DERY(12) + POLL(4,1)+POLL(4,2)+POLL(4,3 )
GO TO 65
OERY(l) = (PQLL(3.l)*PQLL"(3.2))/100.DO
?F !.X-DEP.TH .LE. XDEPRV)_ GO TO 80
"I"F~"(X"".E"Q'. XP'REV")' 'GO'T'O "so" "
STORE = ((XDEPRV-XDEPTH)*XDEPRV)/(CX-XPREV)*XDEPTH*XDEPTH)
MODIFY CONCENTRATIONS DUE TO CHANGE IN MIXING DEPTH
DO 70 I»1,NDIF
DERY(I)"= OERYd )
00 131 I-l.KDIF
.2X11 > .!L_DIRTJJJ
YSAVEfl.I) * Yd)
RETURN
END
01316000
01317000
01318000
01319000
01320000
01321000
01322000
01323000
01324000
01325000
01326000
01327000
01328000
01329000
01330000
01331000
01332000
01333000
01334000
01335000
01336000
01337000
01338000
01339000
01340000
01341000
01342000
01343000
01344000
01345000
01345100
01345110
01345200
01346000
01347000
01347100
01348000
01349000
01350000
01351000
01351200
01351800
01352000
01353000
01354000
01354100
01354200
01354300
01355000
"01356000
01357000
01358000
01359000
OT360000
01361000
01362000
-------
VIII.30
OS/360 FORTRAN H
COMPILER OPTIONS - NAME = MAIN,QPT=01,LINECNT=55*SIZE=OOOOK,
SaURC~E,BCD»NOLIST,NaDECk»LaAD»"MAP,NbEblT, ID,NOX~REF
SUBROUTINE flETEOR (X,Y,I I DO,XTIME,XPREV,XXOOT,YYDOT)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
ARRAYS
MR—12 HOURLY AVERAGES OF SIX VARIABLES FOR EACH STATION.
WIND_SPEED AND DIRECTION, TEMPERATURE. HUMIDITY*, BAROMETRIC
PRESSURE» AND CLOUD COVER".
ANS —METEOROLOGICAL DATA FOR AN INTERPOLATED POINT.
ARADIO--AM RADIOSONDE CURVE
DISINV--DISTANCE INVERSES, INTERPOLATED POINT TO STATIONS.
GIST—DISTANCES SQUARED FROM INTERPOLATED POINT TO STATIONS.
JBAR —BARRIER CHECK INDICATORS, 1 MEANS ELIGIBLE STATION.
PRADIO--PM RADIOSONDE CURVE
RADIO —RADIOSONDE CURVE
RSLOPE—SLOPES OF RADIOSONDE CURVE SEGMENTS
WPOS--NSTA STATION LOCATIONS, THREE DIMENSIONS
1
DIMENSION DIST(32),J3AR(32)»RADIO(15,2),RSLOPE(14),JFUSE(32»32),
DIS INV(8) . VALK8) »V
-------
VIII.31
C
c
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
35 IF (IHT .EQ. 1)
ISONDE = 5
NNRAD = NRAD-1
IF (ITIME-ISONDE)
GO TO 50
40,50,50
SET AM RADIOSONDE CURVE.
40 00 45 L=1*NRAD
RADIO(L»1) ~ ARAOIQ(L,1)
45 RAOIO(L»2) = ARADICXL.2)
GO TO 60
SET PM RADIOSONDE CURVE.
50 DO 55 L=1*NRAD
?ADIO(L»1) = PRADIO(L.l)
55 RADIOIL»2) =• PRADIO(L»2)
ISONOE = 100
60 00 65 L*2,NRAD
65 RSLOPE(L-l) = )
GO TO 75
BARRIER CHECK POINT AGAINST ALL STATIONS.
75 CALL BARIER (X,Y,WPQS,JBAR,NSTA)
IF (ITIME-ISONDE) 100,80,80
SET PM RADIOSONDE CURVE.
80 DO 85 K=1,NR*D
*ADIO(K»1) = PRADIO(K,1)
85 RADIOIK,2) = PRADIO(K.2)
ISONDE = 100
00 90 L=2,NRAD
<)0 RSLOPE(L-l) =• (RADIO
-------
VIII.32
C
C
C
C
C
C
C
C
C
C
C
INTERPOLATION POINT ON A STATION.
I ONE = 1
GO TO 250
120 CONTINUE
NO STATIONS'QUALIFY. ' " ' " ~"
IF (NPTS .EQ. 0) RETURN
SUM CONTRIBUTIONS OF EACH STATION DIVIDED 3Y ITS DISTANCE TO THE
INTERPOLATED POINT SQUARED. DIVIOE JJESULTS BY 1 OVER THE
SU1 OF THE INVERSE DISTANCE'S. SET~RESUL TS INTO 'ANS~'~.
DO 140 J*l,4
00 130 I=1,NSTA
IF (JBAR(I) .NE, 1) GO TO 130
IF CAIR(ITIME+1,I.J) .EQ. O.J GO TO 130
DISINV(J) » DISINV(J) + l./OISTd")
VAL2(J) = VAL2(J)+AIRUTIME+1»I» J»/DIST( I)
CONTINUE
IF (DISINV(J) .NE. 0) GO TO 135
VAL2U) = 0.
GO TO 140
VAL2U) = VAL2(J)/DiSINV( J)""
CONTINUE
DO 1*5 1=1,8
DISINV(I) = 0.
DO 160 J = l,4
DO 150 I=1»NSTA
IF (JBAR(I) .NE. I) GO TO 150
IF (AIR
-------
VIII.33
250
270
ANS(l) = ANS( D/60.
ANS<2) * ANS(2)/60.
GO TO 400
TRAJECTORY POINT ON STATION. SET ANS TO STATION VALUES,
DO 270 J=l,6
VAL2U) = AIR(ITIME+l.IQNE,J)
VAU(J) = AIR( ITIflE. IONE, J)
ANS(J) * VALH J) + (VAL2< JJ-VALH J) )*TIMEHR
HT = WPOSUONE-,3)
C
C
C
C
C
C
C
400
DETERMINE MIXING HEIGHT ANO MIXING HEIGHT DERIVATIVE.
XDEPRV = XDEPTH
XPREV = XTI1E1
IF (IHT .EQ. I) GO TO 430
DO 420 I=l,NNRAD
XHT = (ANS(3»-RADIQ(I»2)+RSLOPE( I )*RADIO(I,1)-SLOPE*HT)/
1 (RSLOPE ( D-SLOPE)
IF UXHT .GT. RAQIOUrll) /AND.
1 (XHT .LE. RADIO! I-t-1.1 )) .AND.
2 (XHT .GE. HT)) GO TO 450
420 CONTINUE
TMIN = 100.
DO 425 I=1»NRAD
IF ((HT .GT. RADIO(I,1M .OR. (RADIO(I»2) .GE. TM1JN"))
TMIN = RADIO( N2)
HMIN - RADIO(1,1)
CONTINUE
IF ((TEMP+SLOPE*(HMIN-HT)) .GE. TMIN) GO TO 450
XDEPTH = A1AX 1 (XMIN»XOEPR V)
XHT = XDEPTH+HT
GO TO 500
-------
VIII.34
OS/360 FORTRAN H
COMPILER OPTIONS - NAME= MAIN*0°T = 01,LINECNT=55»SIZE=OOOOK,
SaURCE,BCD,NaLIST,NODECK,'LOAD,MAP,NOEClT,ID,
SUBROUTINE ULTRAV (RATE)
NOXREF
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
ULTRAVIOLET RADIATION MODULE
COMMON /MYCOM/ TIME. I CLOUD, I OM, I DU» I DD» I OS, I PST, NST A, NST ACK, NRAD,
1 AIR(12,32,6)»WPOS<32,3>,PaS(32,3),SPaS(32,3),HS( 12)»
2 DS(12),VS<12)»TS(12)»OS<12>,SBAROM,STEMP,XOOT»YDOT,
3 ARADIO(15.»2)»PR AD 10(15»2)»ZENITH(12»2).SGR 10(650),
4 F GRID (650), AGRIDK650), AGRID2(650),AGRID4(650), EMIT ( 10) ,
5 XHITE(12)*IHT,XMIN
DOU3LE PRECISION POLL
COMMON' /ALLCOM/ POLL u,3),G°OLL(4,3),ANS(8) ,GOTIME> XOEP'TH, ISOUR,
1 XOEPRV
EQUIVALENCE
1
(TEMP,ANS(3)),(HUMIO,ANS(4)),(BAROM,ANS(5) ),
(COVER,ANS(6)),(HT,ANS(7) ) ,(XHT.ANS(8 ) )
BY
DECLARE ARRAYS AND PRESET DATA
ABSORP--11 VALUES OF ABSORPTION COEFFICIENTS OF N02 AND N204
A IRIS 12 VALUES OF THE AIR MASS AT VARIOUS ZENITH ANGLES
CLOUD 4 AVERAGE FRACTIONS OP INCIDENT RADIATION TRANSMITTED
7 CLOUD TYPES
3ZONE 11 ABSORPTION COEFFICIENTS OF OZONE IN THE ULTRAVIOLET
REGION BY WAVELENGTH, 3000A-4000A
SOLAR U MEAN SOLAR SPECTRAL IRRADIANCE OUTSIDE THE ATMOSPHERE
FOR SPECTRAL 8ANDWIDTHS OF 100A CENTERED AT WAVELENGTHS
3000A-4000A
XMOLE 11 INDICES TF REFRACTION AND MOLECULAR SCATTERING COEFFI-
CIENTS FOR A HOMOGENEOUS STANDARD ATMOSPHERE AT 3000A-4000A
ZENITH—13 VALUES OF TIME OF DAY AND ANGLE OF SUN TO VERTICAL
VAPOR VALUES OF WATER VAPOR PRESSURE FOR DIFFERENT TEMPERATURES
XLIST INDE°ENDENT VARIABLE INTERPOLATION LIST
FLIST DEPENDENT VARIABLE LIST CORRESPONDING TO INDEPENDENT
VALUES IN FLIST
DIMENSION AIRMS(9,2).CLOUD(4,8)•OZONE(11),SOLAR!ll),XMOLE<11),
1 VAPOR(9,2),ABSORP( 11),
2 XLISTZ(50),XLISTA(50),XLISTC(50),XLISTV(50),
3 FLISTZ(50),FLISTA(50)»FLISTC(50)»FLISTV(5C)
DECLARE PRESET CONSTANTS
DUST .==" OUST CONTENT OF ATMOSPHERE
DZ3 =* OZONE CONTENT OF ATMOSPHERE
C- =* REFLECTION COEFFICIENT
ALBEDO == FRACTION OF
PI == 3.14159
TIME TIME OF DAY
TEMP HOURLY TEMPERATURE OF THREE CLOSEST WEATHER STATIONS
BAROM HOURLY BAROMETRIC PRESSURE OF 3 CLOSEST WEATHER STAS
HUMID HOURLY HUMIDITY OF THREE CLOSEST WEATHER STAS (PERCT)
01571000
01572000
01573000
01574000
01575000
01576000
01577000
01578000
01579000
01580000
015B0200
01581000
01561200
01582000
01583000
01584000
Ol585000
01586000
01567000
01588000
01589000
01590000
01591000
01592000
01593000
01594000
01595000
01596000
01597000
01598000
01599000
01600000
01601000
01602000
01603000
01604000
01605000
01606000
01607000
01608000
01609000
01610000
01611000
01612000
01613000
01614000
01615000
O16'l6000
01617000
01618000
01619000
-------
VIII.35
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
COVER HOURLY CLOUD COVER OF THE 3 CLOSEST WEATHER STATIONS
XHT CALCULATED MIXING HEIGHT (M)
ZZ INTERPOLATED ELEVATION AT BASE OF COLUMN (M)
TEMPO RATE OF CHANGE OF TEMPERATURE
BAROMD BAROMETER
XHTD MIXING HEIGHT
DOUBLE PRECISION RATE
DATA ALBEDO /O.O/
DATA DUST /l.O/
DATA G tO. 51
DATA OZ3 /2.2/
DATA PI /3. 14159/
DATA
1 A I RMS/ 15 . ,20. ,30., 40., 4 5., 60., 70. ,30. ,85. ,1,03, 1.06, 1. 15,
2 1.30,1.41,2.00,2.90,
3 CLOUD/1.1,2.0,3.0,4.0,.85,
4 5047454141
5 2424171718
1 OZONE/. 44,. 12, .032, .0085,.
£ SOLAR/. 92,1. 19, 1.37, 1.91,1
£ 3.10/,
5 X,10LE/.530».461,.402». 353,
1 .158/,
£ VAPOR/0. ,5. ,10. ,15. ,20. ,25
3 1.2788,1.7535,2.3756
'CHECK FOR FIRST TIME THROUGH "UL
IF ( IDU .EQ. 1 ) GO TO 250
5.60, 10. 39/,
.89».82».80».34».73,.7l».65».52»
41 41 35 34 32 31 25 25
18/
0020, .0005,5*0 /,
.90,2.07,2. 10,2.48,2.36,2.20,
. 3 H, .275,. 24 5,. 2 18,. 195,. 17 5,
.,30. ,35., 40., . 4579,. 6 54 3, .9209,
,3. 1824,4.2175,5.5324/
TRAV
SET UP ARRAY ZENITH FOR INTERPOLATION ( XL I S TZ , FL I STZ )
SET UP ARRAY AIRMS FOR INTERPOL
ATION (XLISTA,FL IS'TA)
SET UP ARRAY CLOUD FOR INTERPOLATION i XL I STC » FL I S TC )
SET UP ARRAY VAPOR FOR "iNTERPOL
DO 50 1-1,12'
XLISTZU = ZENITHU,!)
50 FLISTZl I = ZENITH( I ,2)
DO 100 =1,9
XLISTAU = AIRMS(I,l)
100 FLISTAl I = AIRMSd, 2)
DO 150 I 1,4
XLISTCU = CLOUDU,!)
150 FLISTCU = CLOUDU, ICLOUD)
DO 200 =1,9
XL~IStV(I = VAPORU,!)
200 FLISTV( I » VAPORU, 2)
IDU = 1
ATIQN (XLISTV.FLISTV)
01620000
01621000
01622000
01623000
01624000
01625000
01626000
01627000
01628000
01629000
01630000
01631000
01632000
01633000
01634000
01635000
01636000
01637000
01638000
01639000
01640000
01641000
01642000
01643000
01644000
01645000
01646000
01647000
01643000
01649000
01650000
01651000
01652000
01653000
01654000
01655000
01656000
0 1657000
01658000
01659000
01660000
01661000
01662000
01663000
01664000
01665000
01666000
01667000
01668000
01669000
01670000
01671000
01672000
-------
VIII.36
C
C
C
C
C
250
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
300
DETERMINE ANGLE OF SUN
AIR MASS AT THIS ZENITH ANGLE
BETA (INCIDENT RADIATION) FOR THIS
PARTIAL PRESSURE OF WATER VAPOR IN
ANGLE = P I F 2 ( T IME, XL I STZ, 1 2* FL I STZ )
CANGLE = COS( ANGLE*PI/180. )
AIRMAS = PIF2(ANGLE,XLISTA,9,FLISTA)
SETA = PIF2( AIRMAS. XL I S TC » 4, FL I STC )
PRESS = PIF2( TEMP»XLISTV»9»FLISTV)
PZ = PRESS*HUMID
W = 2. 3*PZ*10.**(-HT/22000. )
Cl = (l.-CaVER»/CANGLE+2.0*BETA*COVER
C2 = 2.0*( l.-COVER+BETA*COVER)
MASS AND A CLOUD
THE ATMOSPHERE
TYP
(BY HUNDREDS)*
FOR EACH WAVELENGTH IN ANGSTROMS FROM 3000 TO 4000
DETERMINE THE FOLOWING VALUES
A MOLECULAR SCATTERING COEFFICIENT (TA3LE XMOLE)
B TRANSMISSION DUE TO MOLECULAR SCATTERING
C TRANSMISSION DUE TO PARTICLE SCATTERING
0 ABSORPTION COEFFICIENT OF OZONE (TABLE OZONE)
E TRANSMISSION DUE TO OZONE
F MEAN SOLAR SPECTRAL IRRAOIANCE OUTSIDE THE ATMOSPHERE
(TABLE SOLAR)
G DIRECT AND SKY RADIATION
H CONTRIBUTION TO ULT3A VIOLET REACTAMT RATE PARAMETER
FOR ONE WAVELENGTH
THEN SUM THE CONTRIBUTION TO GET THE SPECIFIC ABSORPTION RATE
L = 1
RATE1 = 0.
RATE2 = 0.
RATE * 0.
DO 300 LAMDA=3000»3300> 100
A,8
SCTMOL = -XMOLE(L)*BAROM*AIRMAS/1013.25
XLAMD4 = 1.E+04/LAMDA
SCTPAR = -(3.75E-03*W*XLAMO
-------
VIII.37
OS/360 FORTRAN H
COMPILER OPTIONS - NAME = MAIN, OPT»01» L INECNT«55.» S I ZE=OOOOK,
SOURCE,BCD,NOLIST,NaDECK,La"AD,MAP,NOEDrr/iD,NOXREF
SUBROUTINE SOURCE(X,Y)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
SOURCE EMISSIONS MODULE
COMMON /MYCOM/ TIME,ICLOUD,I DM, IDU,I DO,IDS,IPST,NSTA,NSTACK»NRAO,
I AIR(12,32,6),WPOS(32,3),POS(32,3)»SPOS(32,3),HS( 12),
2 DS(12>,VS<12),TSl12 ) ,QS(12),SBAROM,STEMP,XDOT,YDOT,
3 ARADIO(15,2),PRADIO(15,2),ZENITH(12,2),SGRIO(650),
4 FGRID(650).AGRID1(650),AGR102(650),AGR104(650) »EMIT( 10) ,
5 XHITE<12)»IHT,XMIN
DOUBLE PRECISION POLL
COMMON /ALLCOM/ POLL(A,3)•GPDLL(4,3)»ANS(8),GOT I ME,XDEPTH»I SOUR,
I XDEPRV
DECLARE ARRAYS AND PRESET DATA
POLL ARRAY OF 4 POLLUTANTS, 3 SOURCES
EMIT -EMISSION RATES OF POLLUTANTS
WEIGHT — GRAM MOLECULAR WEIGHTS OF POLLUTANTS
FCURVE--TEMPORAL DISTRIBUTION CURVE OF FREEWAY TRAFFIC
SCURVE--TEMPORAL DISTRIBUTION CURVE OF STREET TRAFFIC
IBOX INDEX TO GRID POPULATION ARRAYS
EPOS MIDPOINTS OF ADJACENT GRIDS TO TRAJECTORY POINT
FGRID FREEWAY MILEAGE DISTRIBUTION BY GRID
SGRID---STREET MILEAGE DISTRIBUTION BY GRID
ACRID AREA SOURCES USAGE DISTRIBUTION BY GRID
DIMENSION WEIGHT(10)»FC.URVE(26),SC'JRVE<26), ACURVE(26) >IBOX(4),
I EPOS<4,2) ,DIST( A),EMITWT( 10)
ORDER OF MOLECULAR WE IGHTS--PROPYLENE,N0»CO,OUMHC
DATA WEIGHT /42.,30..28.,25.,6*0./
DATA FCURVE /.03022*5*.00776,.017B».0591,.0768,.06^d».0536,
1 " A*. 049425,.0569,3*.0746,.0598,5*.03022..00776/
DATA SCURVE /.030775,6*.006783,.0293,2*.065 I,2*.0502,5*.06088,
1 2*.08195,2*.05405.A*.030775,.006783/
IF STANDARD TIME, DON'T RESET
IF (IPST .EQ. 1) GG TO 20
TIME = TIME+1.0
CHECK FOR FIRST TIME THRU ROUTINE
20 IF (IDS .EQ. 1) GO TO 50
AX = .0633
POLLI3.3) = 0.
DO 40 1=1,4
01731000
01732000
01733000
01734000
01735000
01736000
01737000
01738000
01739000
01740000
01741000
01741200
01742000
01742200
01743000
01744000
01745000
01746000
01747000
01748000
017
-------
VIII.38
C
C
C
C
C
C
C
C
C
40 EMITWTU) = EMIT( I) / I WE IGHT( I)*60.)
IDS = 1
50 IHR = TIME+2.0
GMIN = TIME+2.0-IHR
EVALUATE TEMPORAL DISTRIBUTION CURVES
100
IF (GMIN .GT/0.5) GO TO 100
FX = FCURVE( IHR-1)*( ,5-GMIN)
SX * SCURVE( l'HR-l)*( ,5-GMIN)
30 TO 150
FX = FCURVEt IHR)*(1.5-GMIN) *
SX = SCURVEC IHR)*(1.5-GMIN> +
FCURVEJIHR)*(
SCURVE
-------
VIII.39
GO TO 400
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
X POSITIVE,
300 IBOX(2) * I
IBOXO) > I
IBOX(4) * I
IPT = 4
Y POSITIVE
BOX( D+25
BOX(1)>26
BOX( 1)+1
SET MIDPOINTS OF ADJACENT GRIDS
400 GO TO (420,
420 EPOS(2.2) *
EPOS (3,1) =
GO TO 500
440 EPOS(2,2) =
EPOSO,!) =
GO TO 500
460 EPOS(2,2) =
EPOSO,!) =
GO TO 500
480 EPOS(2,2) =
"EPOS (3, f) "=
500 EPOS(2,1) =
EPOSO, 2) =
EPOS(4,l) «
EPOS(4,2» =
440,460,480), IPT
IYM-2
IXM-2
IYM-2
IXM + 2
IYM4-2
IXM-2
IYMt-2
IXM>2
IXM
EPOS(2,2)
EPOSO, 1)
I YM
CALCULATE EMISSION CONTRIBUTIONS TO THE COLUMN OF THE FOUR
ADJACENT GRIDS INVERSELY PROPORTIONAL TO THE SQUARE OF
THE DISTANCE FROM THE. TRAJECTORY POINT TO THE Gfifl 0 MIDPTS,
(Y-IYM)**2
G0 TO 600
540
DIST(l) » (X-IXM)**2
IF (DISjm . EQ. 0.) _ _
DISINV = 1/DISt ( 1) """ " ~ ...... '
STREET = SGRID( IBOX(l) ) /DISTI 1)
FWY = FGRIDC IBOX( 1) )/DIST( I)
IF (ISOUR .NE. 1) GO TO 540
AREA! = AGRIOK IBQX(l) )/DIST(l)
AREA2 = AGRID2( IBOX(l) )/OIST(l)
AREA4 » AGRID41 IBOX(l) i'/DIST (1 ) .....
DO 550 1=2,4
IF HBOXU) .EQ. 0) GO TO 550
aiSTdJ = (X-EPOSU, l» )**2 * (Y-EPOS( 1,2) )**2
01833000
01834000
01835000
01836000
01837000
01838000
01839000
01840000
01841000
01842000
01843000
01844000
01845000
01846000
01847000
01843000
01849000
01850000
01851000
01852000
01853000
01854000
01855000
01856000
01857000
01858000
01859000
01860000
01861000
01862000
01863000
01864000
01865000
01866000
01867000
01868000
01869000
01870000
01871000
01872000
01873000
01874000
01875000
01876000
01877000
01878000
01879000
01880000
01681000
01882000
01883000
01884000
01885000
-------
VIII.40
550
DISINV =
STREET =
FWY
IF (I SOUR
A R E AI =
AREA2 =
A&EA4 =
:ONTINUE
STREET *
FWY
IF (I SOUR
AREA1 =
AR5A2 =
AREA4 =
30 TO 620
620
640
650
680
700
DI_SINV-H/DIST< I) _
STREET + SGRIDUBOXt I) )/OIST(I )
FWY + FGRIDJ IBOXU) )/DIST< I )
.NE. 1) GO TO 550
AREA! + AGRIDK I30XU ) )/DIST( I )
AREA2 + AGRID2( IBOXt IM/DIST< I)
AREA4 + AGRID4J IBaxU ) )/OISTUJL
STREET/OISINV
FWY/DISINV
.NE. 1) GO TO 620
AREA1/DISINY
AREA2/DISINV
AREA4/DISINV
600
TRAJECTORY ON MIDPOINT OF GRID
STREET = SGRIDl I BOX(1))
FWY = FGRID(I30X(1))
IF (ISOUR .NE. i) GO TO 620
AREA1 = AGRIDU I8QX( 1) )
AREA2 = AGRID21IBOX(I))
AREA* = AGRIDAtI80X(1))
CALCULATE EMISSIONS FROM SOURCES
STTEMP = SX*STREET/10359386.
FTEMF = FX*FWY/10359396.
IF (ISOUR .NE. 1) GO TO 640
ATEMP1
ATEMP2
ATEMP4
FACTOR
DO 650 1=
GPOLLU »1)
GPOLLi 1.2)
POLL(I.I)
POLL( 1.2)
IF (ISOUR
GPOLL(1>3)
5POLL(2*3)
GfOLL(4.3>
POLL(I.3)
POLL(2»3)
POLL(4,3)
IF (IPST .
AX*ARJA1/10359386.
AX*AREA2/10359386.
AX*AREA4/10359386.
4.0E6/XOEPTH
= STTEMP*E'1ITWT( I )
= FTEMP*EMITWT(I)
• GP'OLU I » 1)*FACTOP
« GPOLLII»2)*FACTOR
, NE." 1) GO TO 680
« ATEMP1*EMITWT(1)
= ATEMP2*E1ITWT(2)
EQ
_....
GPOLL"("i.3)*FACTQR"
GPOLL(2,3)*FACTOR
GPOLL<4»3)*FACTOR
1) GO TO 700
TIME =
RETURN
END "~
TIME-1.0
01886000
01887000
ouaaooo
01889000
01890000
01891000
01892000
01893000
01894000
01895000
01896000
01897000
01898000
01899000
01900000
01901000
01902000
01903000
01904000
01905000
01906000
01907000
01908000
01909000
01910000
01911000
01912000
01913000
01914000
01915000
01916000
01917000
01918000
01918500
01919000
01920000
01920200
01921000
01921200
01922000
01923000
01923200
01924000
01924200
01925000
01925200
01926000
01927000
01928000
01929000
-------
COMPILER
C
C
C
C
C
C
C
C
C
VIII.41
OS/360 FORTRAN H
OPTIONS - NAME3 MAJJ1,OPT*01,LLNECNI>55,S!ZE = P
SaURCE,BCD*NOLiST,NaDECK»LOAD»MAP,NbEOIT,ID,NOXREF
FUNCTION PIF2 (X,XL IST,N,FLIST»
SECOND ORDER ONE VARIABLE INTERPOLATION
DIMENSION XLISTUOO), FLIST(IOO)
BLI(FLO,FL1,FL2,FL3.FL4") = (( FL 1-FLO )*(FL 3-FL 4) / ( FL~2-FL 1) + FL3 )
IF (X-XLISTtN) ) 2,1,1
PAST END OF INDEPENDENT VARIABLE LIST.
1 I - N-l
GO TO 5
2 IF (X-XLISTUM 4, 4,6
BEFORE START OF INDENOENT VARIABLE LIST.
4 I « 1
5 < = 1
GO TO 30
6 X = 2
DETERMINE INDEPENDENT VARIABLE INTERVAL.
7 DO 3 1*1,N
IF (X-XLISTlI)) 9»9,8
8 CONTINUE
I = N
91=1-1
30 BLIF1 = BLKX.XLISTl I),XLIST( 1*1 ),FLIST(I ),FLIST(I*l) )
10 IF (K-I) 11,11,12
11 PIF2 = BLIF1
RETURN
FINO BEST THIRD POINT.
12 IF (( I*2)-N> 13,13,16
13 IF< (I-l)-l) 15,14,14
14 IF(ABS(XLIST( I-1)-X) - A8S(XL IST(I+2 ) -X) ) 16,15,15
15 L = 1+2
30 TO 17
16 L = 1-1
3LIF2 = (BLI(X»XLIST(i),XLIST(L),FLIST(I),FLIST(L»))
PIF2 = BLI(X,XLIST(I+l),XLIST(L),BLIFl,BLIF2)
18 RETURN
END
01930000
01931000
01932000
01933000
01934000
01935*000
01936000
01937000
01938000
01939000
01940000
01941000
01942000
01943000
01944000
01945000
01946000
01947000
01948000
01949000
01950000
01951000
01952000
01953000
01954000
01955000
01956000
01957000
01958000
01959000
019t>0000
01961000
01962000
01963000
01964000
01965000
01966000
01967000
01968000
01969000
01970000
01971000
01972000
01973000
01974000
-------
VIII.42
OS/360 FORTRAN H
COMPILES OPTIONS - NAME = MA IN, OPT=01, LJNECNT*55, S I ZE = OOOOK,
SOURC E~,BCD, NOL I Sf",~NODECK» LOAD /MAP,NOE DIT, ID, NOXRE F
SUBROUTINE BARIER (X,Y,WPCS,JBAR,NPOS) 01975000
C 01976000
C DETERMINES IF INTERSECTION OF LINE_ SEGMENT CONNECTING TRAJECTORY 01977000
C POINT AND STATION AND THE BARRIER' LINE SEGMENT LIES ON THE BARRI ECO 1978000
C ITSELF. EVALUATE FUNCTIONS ALPHA AND BETA, IF THEY ARE BOTH LESS01979000
c" "
r
\*
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
50
100
c
c
c
200
c
c
'THAN'! ANO'GTRE'ATER THAN o, STATION Fs ELIMINATED.
COMMGN /BARRR/ IDS
DIMENSION WPOS(32,3), JB AR < 32 )» B AR < 10, 10 ) , 8SLOPE < 10, 10 ) »
I TSLOPE(32,2)
SBAR IS NUMBER OF BARRIERS
ARRAY BAR ODD COLUMNS CONTAIN X COORDINATES
EVEN COLUMNS CONTAIN Y COORDINATES
ROW I IN HDD COLUMNS CONTAIN NO. OF POINTS IN BARRIER
ARRAY BSLOPE--ODD COLUMNS CONTAIN M VALUES OF BARRIER LINE SEGMTS
01980000
01981000
01982000
01983000
01984000
01985000
019~86000
01987000
01988000
01989000
01990000
EVEN COLU1NS CONTAIN Y VALUES OF BARRIER LINE S EGMENTSO 1991 000
ARRAY TSLOPE--ODO COLUMN CONTAINS ~M VALUES OF TRAJECTORY POINT
WEATHER STATION LINE SEGMENT
EVEN COLUMN CONTAINS ~B VALUES OF TRAJECTORY POINT
WEATHER STATION LINE_SEGMENT
DATA N3AR/ 2/» B AR/ 2. » 16. » 38. > 8*0. , 2* 40. , 7*0. , 4 . * 53. .» 62 . , 66 . » 71 . ,
1 " 6*0. .31 ,,27.,32./23. ,65*0. /
DETERMINE SLOPE AND INTERCEPT FROM INTERPOLATION POINT TO EACH
STATION.
DO 50 1=1, NPOS
TSLOPE( 1,1) = 0.
TSLOPEt 1,2) = 0.
DO 100 1=1, NPOS
IF( (WPOSI I»1)-X) .EQ. 0.) GO TO 100
TSLOPE(I.L) = (WPOS ( I ,2)-Y)/(WPOS( I.D-X)
TSLOPE(I,2) = Y-TSLOPE( I,1)*X
CONTINUE
DETERMINE SLOPE AND INTERCEPT OF EACH SEGMENT OF EACH BARRIER
IF ( IOB .EQ. 1) GG TO 210
00 200 J=l,NBAR
NN3AR = BAR(l,2*J-l)
DO 200 K=2,NNBAR
IF ( (BAR(K<-i,2*J-l)-BAR(K,2*J-l) ) .EO. 0.) GO TO 200
3SLOPE(K,2*J-1) = (3A3(K+l,2*J)-BAR(K,2*J) )/
1 (8AR(K+1,2*J-1)-BAR(K»2*J-1) )
BSLOPE(K,2*J) = BAR(K,2*J)-BSLOPE(K,2*J-1)*BAP(K,2*J-1)
CONTINUE
IDS = I
DETERMINE INTERSECTION FUNCTIONS
01992000
01993000
01994003
01995000
01996000
01997000
01998000
01999000
02000000
02001000
02002000
02003000
02004000
02005000
02006000
02007000
02008000
02009000
02010000
02011000
02012000
02013000
02014000
02015000
02016000
02017000
02018000
02019000
02020000
02021000
"02022000
02023000
02024000
02025000
-------
VIII.43
IF VALUE OF INSECTJON IS LE
INCLUDE STATION.
210 00 400 I»1.NPOS
DO 375 J=l,NBAR
NNBAR = 8AR(1»2*J-1)
00 375 K=2. NNBAR
TO ZERO OR GE 1*
GO
TO 220
1,2) )
IF (2*-m EQ. 0.)
DENOM1 = BSLOPE(K,2*J-l)*(X-WPOSU»l) >-(Y
IF (DENOM1 .EO. 0.) GO' TO 350
ALPHA = (WPOS( I,2)-BSLOPE
-------
VIII. 44
OS/ 360 FORTRAN H
COMPILER
R OPTIONS - NAME= MAIN, OPT=Ol, L I NECNT*55, S I ZE=OOOOK,
SOURCE>BCD»NOLIST»NaDECK,LOAD»MAP»NOEDIT*ID*NOXREF
SUBROUTINE SMGPLT
DIMENSION OUT( 101), YPR( 11)
DOUBLE PRECISION C, D, Z, DERY., ZZ, PRMT, BB, POLL , YY
COMMON /KINCOM/ ZZ ( 4 ) , YY( 40 ) » L ABEL < 40 ) » D ( 50 ) » RC T ( 50* 3 ) » PRMT ( 8 ) ,
1 DERY(40),C( 50),BB,Z(2),N,NC,NEQ,NDIF,KDIF,M,NN, INTR,NNC,
2 NNDIF, IUV, IM, I S, I TR AF/NLL , CC fl&, 100 ) , A i'( 80)"* STEfM 1 5 ) » T I Ml,
3 TIM2,TIM4,PEAK1,TIME1,PEAK3,TIME3»XPOS*YPOS
COMMON /ALLCOM/ POLL ( 4, 3 ) ,GPOLL ( 4, 3 ) ,ANS(8) , GOT I ME ,XOEPTH, I SOUR,
1 XOEPRV
REAL*8 LABEL, RCT
REAL BLANK/ • ' / , ANG ( 16 ) / • X' * • 1 • * • 2 • » ' 3 ' .» • 4 • * • 5 • , * 6 • , • 7 • , • 8 ' , f 9 •
2'At,'B',lC',lDl,'E'»lF'/,IB/Mt/
YMIN = 0.
DO 5 1=1,101
5 OUT( I ) = 0.
DO 6 1=1,11
6 YPRU ) = 0.
'1M = NDIF+1
NN=NN-l
WRITE(6. 9) NDIF,NN
IF(CC(2,NN).NE.PEAK1 . AND. CC ( 3 . NN ) . NE . PE AK 3) GO TO 10
WRITE (6*115)
10 WRITE (6,110) PEAK1,TIME1,PEAK3,TIME3
IF(TIMl.NE.O.O) WRITE (6,120) TIM1
IF(TI*2.NE.O.O) WRITE (6,130) TIM2
IF(TIM4.NE.O.O) WRITE (6,140) TIM4
11 NGRAPH=1
IF_(MM-9) 15, 15> 12
12 NGRAPH=2
IF (MM-14) 14,14,13
13 NGRAPH=3
14 JMAX»6
SO TO 16
15 JMAX=MM
16 JMIN=2
IF(NLL) 18, 18,20
Ifl ^)LL = 51
20 IF(NLL.LT.NN)NLL=NN+l
AA=NLL-1
XS.CAL = (PRMT(2 )-PRMT( 3) ) /A A _
IF" "(PRMT(2) .GT. CC(i,NN)*-5".0) XSC AL= (CC ( 1 , NN ) -PRMT ( 8 ) ) /AA
XSCLHF=XSCAL/2.0
DO 95 K=1,NGRAPH
IF(NGRAPH-K) 21*21*22
21 JAX=MM
GO TO 23
22JAXSK*5*1
23 JIN=(K-l)*5+2
125 YMAX=6.0
DO 40 1=1, NN
02065000
02066000
02067000
02068000
02069000
02070000
02071000
02072000
02073000
02073200
02074000
,0"2075000
02076000
02077000
02078000
02079000
02080000
02081000
02082000
02083000
02084000
02085000
02086000
02087000
02088000
02089000
02090000
02091000
02092000
02093000
02094000
02095000
020°6000
02097000
02098000
02099000
02100000
02101000
02102000
02103000
02104000
02105000
02106000
02107000
02108000
02109000
02110000
02111000
02112000
02"! 13 000
02114000
-------
VIII.45
30
40
PO 40 J.= JMIN, JMAX
IF(CC(J,I)-YMAX) 40,40,30
YMAX * CC(J»I )
CONTINUE
YSCAL*YMAX/100.
XB-CCC1,!)
L- 1
1=1
WRITE (6,1)
WRITE (6,7)
45
46
47
48
50
55
56
57
59
60
70
80
66
90
95
99
i
2
3
PRMT(l).PRMTm,PRMT(5),PRMT(6)»A
XPR=XB + F* XSCAL
IF (XPR-CCU.L )-l
0) 43,48,47
GO
IF
IF
00
50, 50» 70
TO 46
,Eio.4,' MIN SIMULATION'
L» PR^T(4)=«,E10.4,/,4( • «,20A4,M)
FORMATC1H ", F 8 . 1, 2X, 101 Al, F 7. 2 )
FORMAT (1H »Ffl.l,2X, •IV,99X,« I')
BB=',E10.4,
PRMT(3)=',EIO
02115000
02116000
02117000
02118000
02119000
02120000
02121000
02122000
02123000
02124000
02125000
02126000
02127000
02128000
02129000
02130000
02131000
02132000
02133000
02134000
02135000
02136000
02137000
02135000
02139000
02140000
02141000
02142000
02143000
02144000
021*5000
02146000
02147000
02148000
02149000
02150000
02151000
02152000
02153000
02154000
02155000
02156000
02157000
02158000
02159000
02160000
02161000
02162000
__ 02163000
,4,02164000
02165000
02166000
02167000
-------
VIII.46
7 FORMAL*
8 FOR MAT ('
9 FORMAT!
IF DATA.
100 FORMAT!
110 FORMAT!
I/ ,
115 FORMAT!
120 FORMAT!
130 FORMAT!
140 FORMAT!
150 FORMAT!
END
1H _»J: 9X» 1.0 U OH* . l>. IH*
IH ,3"xVllF"l6".2) ~
'OTHERS ARE'»I3,' DIF. CONCN.
•,'MM WAS CHANGED TO 10, NN WAS CHANGED TO 50.')
N02 PEAK CONCN PPHM WAS "' »E 13". 7, ' AT TIME "'*E9~. 3
•»E9.3,' MIN.')0217AOOO
0^168000
02169000
STEPS IN THIS SET 002170000
02171000
02172000
',02173000
MIN*
1 03 PEAK CONCN PPHM WAS t,E13.7,1 AT TIME
IH ,'BOTH PEAK VALUES NOT
•OHALF-TIME FOR N02 IS
•OHALF-TIME FOR NO IS
•OHALF-TIME FOR OLEFIN
1H1)
REACHED, HIGHEST
MINUTES.1)
«,F7.3, ' MINUTES. • )
IS ',F7.3,' MINUTES.')
VALUES WERE')
02175000
02176000
02177000
02173000
02179000
02180000
-------
IX. 1
Section IX
Operation of the Reverse Trajectory Program
The reverse trajectory program calculates air parcel paths beginning with
the terminal position. It uses hourly average wind information, and inter-
polates initial pollutant concentrations from hourly pollutant monitoring data.
Starting values are computed for nitrogen dioxide, nitric oxide, ozone, and
carbon monoxide. Trajectories may also be calculated in a forward direction
with the pollutant interpolation occurring at the end point. In this case,
the values computed for N0_, NO, 0. and CO amy be considered for validation
purposes.
The program's operation, subroutine by subroutine follows:
1. MAIN
a. Reads in all input data to the program.
b. Converts all data to the metric system.
c. Initialization function.
d. Prints meteorological and pollutant data.
e. Calls SMGDRV.
f. Calls POLLUT.
g. Outputs pollutant interpolation values.
h. Checks for multiple runs.
2. SMGDRV
a. Serves numerical integration function.
b. Controls processing end time.
c. Calls RWIND.
d. Calls PUTOUT.
3. RWIND
a. Calls BARTER.
b. Determines wind velocity inputs for the particular trajectory
point.
-------
IX. 2
c. Determines wind velocity rates of change (derivatives) of wind
velocity data.
4. BARIER
a. Determines if a natural barrier blocks a meteorological station
or a pollutant monitoring station from the trajectory point in
consideration. If so, the station is discounted.
5. PUTOUT
a. Prints trajectory positions, time into trajectory in minutes.
6. POLLUT
a. Calls BARIER
_b. Determines interpolated concentrations of NO-* NO, 0_, and CO at
final trajectory point.
-------
X.I
Section X
Reverse Trajectory Inputs
Inputs must be prepared exactly as described. A sample input deck is shown
in Section X.
1. Card Type 1
Heading data array A 4 cards
Alphanumeric data used as heading information. May be blank cards.
o 8
l o
2. Card Type 2
Starting data—7 Items—1 Card
XPOS - X coordinate of starting position in miles. Range: 17.0 to
67.0; cols. 1-10.
YPOS - Y coordinate of starting position in miles. Range: 1.0 to
53.0; cols. 11-20.
GOTIME - Simulation start time (local standard time) in hours and
decimal fraction of hours. Range: 7.000 to 18.999; cols. 21-30.
PRMT(2) - Length of simulation in minutes and decimal fraction of
minutes. Cols 31-40. For backward trajectory, value must be
negative. For forward trajectory, value must be positive.
PRMT(3) - Initial value of step size for numerical integration in
minutes and decimal fraction of minutes. For backward trajec-
tory, value must be negative. For forward trajectory, value
-------
X.2
must be positive. Normal range in magnitude is from .001 to
.5; cols. 41-50.
PRMT(4) - Upper bound of error in numerical integration. Determines
allowable difference between predictor and corrector. Normal
range is from .5 to .01; cols. 51-60.
TINT - Output interval. Specifies in minutes how often position print
outs are to occur. Normal range is from 5.0 to 30.0 minutes;
cols. 61-70.
FORMAT (4F10.3, 3F10.5)
0 11 22 33 44 55 66 77 8
1 01 01 01 01 01 01 01 0
XPOS
YPOS
GOTIME
PRMT(2)
PRMT(3)
PRMT(4)
TINT
3. Card Type 3
Starting data—4 Items—1 Card (all values must be right justified)
NSTA - Number of meteorological stations; cols. 1-10.
NDIF - Number of accumulating species, usually 0; cols. 11-20.
NPOL - Number of pollutant monitoring stations; cols. 21-30.
IN - Indicates if program is to operate again with data from
previous run or a new set of data is necessary.
0 - A complete set of data is necessary for the run.
1 - Rerun with previous data. Must input only card types
1, 2, and 3; cols. 31-40.
-------
X.3
0 11 22 33 J44& 8
1 01 " 01 01 '^OIP 0
NSTA
NDIF
NPOL
IN
4. Card Type 4
Wind velocity data—1 array—maximum of 32 stations.
Use card type 9 from model, or set up cards as follows:
AIR - 3 cards per station, 4 hours per card, 12 hours per station.
IDNUM - Column 1-3, alphanumeric station ID
AIR - Column 4-5, start hour of data, values are 7, 11 or 15,
Column 6, number of readings on card, max = 4
Column 7-12
Column 24-29
Column 41-46
Column 58-63
2 items in the 6 columns, 3 columns per item wind direction in degrees,
wind speed in miles/hour, last column is in tenths of a mile/hour.
Conversion: wind speed and direction are converted into XDOT and
YDOT, the X and Y components of wind velocity at the trajectory point.
FORMAT (A3, 3X, 4(F3.0, F3.1, 11X))
The station order on card types 4 and 5 must be exactly coincidental.
hourly data—6 columns
-------
X.4
0 00000 01 11 22 22 23 44 44 44 55 66 66 8
1 34567 90 23 34 67 90 01 34 67 78 01 34 0
I
D
H
0
u
p
4
W
D
I
R
W
S
P
D
W
D
I
R
W
S
p
D
W
D
I
R
W
S
P
D
W
D
I
R
W
S
P
D
5. Card Type 5—(See Figure 4)
Meteorological Stations—2 arrays—maximum of 8 cards, 4 stations per card.
WPOS - Array of quadruplets, 4 stations per card, maximum of 32 stations.
Cards are as follows:
IDMET - Column 1-3, alphanumeric station ID code.
WPOS - Column 4-8, X coordinate of station position in miles, last
2 columns in hundreths of miles. Column 9-13, Y coordinate of
station position in miles, last 2 columns in hundreths of miles.
Column 14-20, ground elevation of station in feet, last
2 columns in hundreths of feet. This 20 column grouping is
repeated in columns 21-40, 41-60, 61-80 as required, and on
additional cards if necessary.
FORMAT (4(A3, 2F5.2, F7.2))
The station order on card types 4 and 5 must be exactly coincidental.
0 00 00 11 22 22 22 33 44 44 44 55 66 66 66 77 8
1 34 89 34 01 34 89 34 ' 01 34 89 34 01 34 89 34 o
I
D
X
Y
ELEV
I
D
X
Y
ELEV
I
D
X
Y
ELEV
I
D
X
Y
ELEV
-------
X.5
6. Card Type 6
Pollutant monitoring data—2 arrays—maximum of 5 cards, 4 stations
per card.
PPOS - Array of quadruplets, 4 stations per card, maximum of 20 stations.
Cards are as follows:
IDPOL - Column 1-3, alphanumeric station ID code.
PPOS - Column 4-8, X coordinate of station position in miles, last
2 columns in hundreths of miles.
Column 9 - 13, Y coordinate of station position in miles, last
2 columns in hundreths of miles.
Columns 14-20, blank.
This 20 column grouping is repeated in columns 21-40, 41-60, 61-80 as
required, and on additional cards if necessary.
FORMAT (4(A3, 2F5.2, F7.2))
The station order card types 6 and 7 must be exactly coincidental.
0 K) 00 11 22 22 22 33 44 44 44 55 66 66 66 77 8
1 34 89 34 01 34 89 34 01 34 89 34 01 34 89 34 0
I
D
X
Y
D
X
Y
D
X
Y
I
D
X
Y
7. Card Type 7
Pollutant monitoring data—2 arrays—maximum of 20 stations.
CONG - 4 cards per station, 3 hours per card, 12 hours per station.
Set up card as follows:
-------
X".6
IDPOL - Column 1-3, alphanumeric station ID code.
CONG - Column 4-5, start hour of data, values are 7, 10, 13 or 16,
6, number of readings on card, max =3.
7-27
28-48 hourly data—21 columns, all values are in pphm
49-69
5 items in the 21 columns as follows:
4 columns - N0», last column is in parts/billion.
4 columns - NO, last column is in parts/billion.
4 columns - 0_, last column is in parts/billion.
4 columns - blank, field not used but has been set aside
for input of hydrocarbon data or another pollutant
of interest.
5 columns - CO, last column is in parts/billion.
0 00000 11 11 11 22 22 33 33 34 44 44 55 55 66 66 , 67 8
1 34567 01 45 89 23 78 12 56 90 34 89 23 67 01 45 90 Q
I
D
H
0
U
R
3
N02
NO
°3
CO
N02
NO
°3
CO
N02
NO
°3
CO
FORMAT (A3, 3X, 3(3F4.1, 4X, F5.1), 11X)
The station order or card types 6 and 7 must be exactly coincidental.
-------
XI. 1
Section XI
Reverse Trajectory Sample Input Deck
Card Type 1
SEPTEMBER 29, 1969 11 30 AM
TERMINATE AT PASADENA
9 KCTfOKOlOGlCAL. STATIOf**r~*
Z POtLUfANT STATIONS I
Card Type 2
42.0
'13.5
15.0
Card Type 3
\z
Card Type 4
W8114
NT8154
NZJ074
NZJ114
NZJ154
CAP074 90 2
CAP114270 8
CAP154270 8
ONT074
ONT114
ONT154
IAX074350 5
LAX114240 9
LAX15426016
BUR074 7
BUR 114 1
BUR154315 4
L6B074315 2
LCB 114270 5
L6B15427013
ENC074113 2
ENC114293 3
ENC154293 6
EtflOTAUO 3
EU* 11*220 7
Eli §115*2*01 2
67611114
86621111
89561094
63631101
88561107
90581091
70
87
87
6854
9542
9952
65621113
80591116
72621100
7151
9246
10047
67611122
87571119
86571102
7459
9453
9539
225 3
270 9
270 8
350 5
26014
26016
1
2
315 4
270 3
270 8
27010
158 2
293 5
270 6
100 1
220 9
25Q10
90611105
84541093
92571101
88591093
75
88
86
7054
9745
9760
69631118
74621114
71601099
7356
9839
9642
71611125
91561112
81581102
8055
97S2
93L»*
270 5
270 9
270 5
240 9
26014
26016
1
180 3
338 2
270 4
27010
270 7
293 3
293 7
270 5
140 2
22010
Zfi9.lt,
78611117
88601099
82521096
77621113
92611096
86601093
79
86
82
10048
9248
72641121
73621109
71601096
7756
9950
9143
77611125
89581108
77571102
8557
98*6
86 J* ,.„.
270 6
270 9
270 6
240 9
26014
230 8
1
270 4
270 1
293 7
27010
270 3
293 3
293 6
270 3
230 4
22011
,.,210 8
80621114
92581095
76521099
84601113
92581091
77591093
BO
86
78
9047
10056
8638
74631119
73621103
70581096
8849
9950
9038
84601125
86591105
72551105
9055
94*3
•JSO
-------
XI. 2
Card Type 4 (Continued)
COW074 45 3
COW114293 9
COM15429315
CPK074203 1
CPK114135 1
CPK154338 6
AZU074315 2
AZU114203 5
AZU15420310
WHT074338 4
WHTU4248 5
WHT154270 6
WLA074270 1
ULA114225 4
WLA154225 4
WNT074 23
MNT114248
WNT154203
VER074270 (>
VERU4270 3
VER15422511
VEN074293 1
ijfEN 114225 6
'\TEN154248ll"
RVA074 23 3
«VAll*203 5
RVA154248 7
RSD074113 1
RSD114203 1
RSD154270 7
R8 074293 3
RB 114203 5
SB 154225 8
RLA074225 1
RLA114225 1
RLA154225 1
POM074U3 2
POM114315 4
POM154293 6
PIC074315 1
PICU4225 1
PIC154180 1
MAL074225 2
HAL114225 5
MAL15422510
KFI074 23 1
KF1114293 4
KFI15431510
HOL074270 1
HOL114225 1
HOL154225 2
315 3
29311
27D17
135 1
203 2
338 6
90 2
203 7
20311
293 2
270 6
270 6
203 1
225 5
203 4
158
248
248 8
270 2
270 5
270 5
248 2
24810
248 8
338 3
225 6
225 7
248 1
68
293
248
225
203
360
225 1
225 1
135 2
315 4
315 6
360 1
225 1
1«0 1
225 3
225 7
225 8
270 1
315 6
315 9
135 6
180 1
225 2
293 4
29314
29316
180 1
158 2
338 4
113 2
225 7
203 9
180 2
270 6
270 4
225 I
225 4
203 2
293
203
248 7
270 2
225 3
270 9
270 5
,24810,
225 6
23 3
248 8
248 6
180 1
135 2
248 6
248 4
225 7
203 5
270 1
270 1
225 2
248 2
293 4
293 4
225 I
225 1
158 2
225 5
90 2
225 I
293 8
293 7
270 I
225 9
,225 6
293 9
293 IS
29312
135 1
315 4
315 2
68 1
203 9
203 5
248 2
270 6
248 3
225 2
225 3
203 1
158
203
180 4
270 3
225 5
27010
248 3
248JU.
203 7
180 4
248 9
225 5
68 1
90 2
270 3
225 3
225 8
203 6
180 1
225 1
225 2
270 2
315 6
293 5
225 1
180 1
203 2
225 8
90 3
293 2
31510
293 5
180 1
225 5
180 2.
-------
XI. 3
Card Type 4 (Continued)
BKT074 90 5
BKT11427Q 6
BKT15427011
PAS074338 2
PAS114180 2
PAS154225 6
CMC074 60 3 7159
CMC114240 7 9355
CMC 15426011 9852
90 3
270 7
27012
45 1
203 3
225 3
150 2 7957
26010 9450
27010 8354
90 4
270 8
270 9
225 I
180 3
225 2
190 5 8558
26012 9152
250 8 7854
270 4
27010
270 7
158 2
225 4
45 1
220 8
26012
230 5
8955
9148
7455
Card Type 5
MTB 51
LAX 32
ELM. 52
MHT 51
VEN 26
RLA 45
KPI 53
,E*C 44
16
27
:^5
3C
30
26
23
31
27
99
300
250
5
100
100
180
NZJ 70
BUR 33
COM 39
WLA 30
RVA 48
POM 70
HOL 35
8
45
22
34
30
33
37
380
530
100
200
180
850
300
CAP 40
LGB 45
CPK 19
WNT 63
RSD 22
PIC 33
BKT 66
33
18
46
30
46
34
36 .
300 ONT
34 ENC
BOO AZU
950 VER
333 RB
100 MAL
950 PAS
76 35
23 44
59 41
41 31
31 21
18 34
47 42
916
800
600
180
10
20
870
Card Type 6
" .CH.A 40
LGB 45
PAS 47
35
18
42
AZU 59
RSD. 22
WHT 51
41
46
30
BUR 33
- POM TO
CMC 44
45
33
31
WLA
LAX
ELM
30 34
32 27
52 35
Card Type 7
"OLA073
OLA103
OLA133
OLA163
AZU073
AZU103
AZU133
A2U163
BUR073
BUR 103
BUR133
BUR163
MLA073
MLA103
MLA133
MLA163
LGB073
LG8103
LGB133
L 68 163
3 23
18 11
5
7 5
13 4
0 1
7 1
17 34
28 7
10 2
12 35
14 3
8 1
9 3
9 26
20 9
19 16
11 21
2
12
7
2
3
22
38
7
2
11
14
3
1
10
10
2
1
2
1
1
17
3
6
10
8
9
, * •-,
Ifl
11
7
5
17
4
5
13
6
5
5
11
10
7
16
5
.V 7
25
19
11
21
8
10
17
22
6
:<»
36
5
2
6
3
1
1
27
3
3
38
1
5
18
8
18
21
5
9
5
1
9
24
25
2
4
14
11
2
3
12
6
1
1
6
1
1
18
6
2
5
11
6
6
•_>, , 5-
17
10
7
6
17
5
4
4
11
6
5
5
19
6
5
11
5
11
29
11
17
19
8
9
28
13
8
30 8
8
3 3
8 1
1 16
28
1 14
1 .A .
14 7
21
1 5
12 2
6 5
10
2 4
1 3
10 2
I
21 1
20 I.
*
3
6
4
10
7
6
7
9
11
5
12
8
^
4
9
5
4
-------
XI. 4
Card Type 7 (Continued)
RS0073
RS0103
KSD133
RSD163
POH073
POM 103
POM133
POM163
LAX073
LAX103
LAX133
LAX163
PAS073
PAS103
PAS133
PAS163
WHT073
MHT103
WHT133
WHT163
CMC073
CMC103
CMC133
CMC163
ELM073
ELM103
ELM133
ELN163
13
10
3
18
10
13
10
11
7
10
6
5
5
19
41
10
12
16
25
9
4
11
1
1
6
1
2
20
4
3
16
1
1
1
34
2
1
2
32
7
9
8
4
14
9
4
3
19
21
6
2
4
8
4
2
26
20
6
2
21
19
4
2
13
14
2
16
6
3
5
9
6
5
6
9
6
3
6
14
13
1
1
17
9
3
4
6
3
21
4
3
21
10
5
13
13
5
8
9
12
8
6
30
21
9
9
26
14
7
3
7
1
1
4
1
1
2
9
1
3
2
1
2
2
9
1
1
1
22
6
6
22
9
9
10
6
3
9
14
21
6
2
10
5
4
9
41
12
3
6
34
12
5
9
2
2
2
11
4
3
5
9
5
6
5
5
4
5
11
1
1
15
6
3
3
8
9
2
4
17
2
4
15
8
14
10
6
9
6
6
9
37
10
11
32
10
7
5
1
1
1
1
1
3
6
2
5
1
1
5
3
1
1
15
3
9
12
17
7
2
16
15
16
3
3
7
4
1
17
32
9
2
11
30
7
3
5
1
6
23
5
3
8
5
5
6
/
4
6
6
5
3
6
3
11
2
1
1
15
4
4
3
6
5
2
5
-------
Section XII
Reverse Trajectory Program Listing
-------
XII. 1
G LEVEL
C
20
MAIN
DATE - 72342
21/45/45
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
£
C
C
C
C
C
C
C
C
BD0001
DI M E N S ION" Y ( 4 0 ) , D E R Y (40)VI 0 N U M~( 32), I~R ( 5 0, 5)"» 6 A I fiI (12!, 3 2» 2!) , P R M T ( 6 )., 8 0 6 0 0 2
1 A(80),IDMET(32),IDPOL(32)
DOUBLE PRECISION AUX,C,D.DELT,DERY,H,X,Y,YY,Z*ZZ
COMMON /RECOM/ TIME*ANS(fl).OERIV(8)* NSTA,NPOL»AIRC12»32»2 )»
1 WPOS(32,3)»GOTIME
COMMON /BARRR/ IDB
COMMON /OUTCOM/ TOUT>TINT
COMMON /MILES/ COORD(2).DXY(2)
REAL*4 N02»NO
DIMENSION CONC( 12»20»5)»PPOS(32»3»»ID(20»
GOTIME —SIMULATION START TItfE
EXTERNAL RWIND»PUTOUT
INITIALIZATION
INPUT HEADING INFORMATION
1 READ <5»1000»END«900) A
DO 10 I=l»40
Y(I) = Q.ODO
10 DERY(J) = 0.ODO
TOUT = 0.0
N = 0
IDB * 0
SIZE* UPPER BOUND ERROR LIMIT* AND OUTPUT INTERVAL
BE0003
BD0004
800005
3D0006
BD0007
BD0003
800009
B00010
BD0011
BE0012
800013
800014
B00015
BD0016
BD0017
BD0018
8D0019
B00020
BO0021
BD0022
BD0024
BD0025
BO0026
~BD002>
800028
800029
800030
BD0031
BD0033
B00034
BD0035
READ (5.1001) XPOS»YPOS»GOTIME»PRMT(2)»PRMT(3).PRMT(4)»TINT
INPUT START POS IT I ON.START TIME, LENGTH OF SIMULATION, INTIAL STEPBD0032
75 PRMT(l) * 0. BD0036
READ NO. METEOR. STATIONS, NO. Of EOUAT., NO. OF POLLUTANT STA. 8E003620
8E003640
800037
READ (5,1002) NSTA,NDIF,NPOL»IN BE0038
BD0039
B00040
BDOO^l
SET UP POSITION FOR LNTEGRAT ION 8D0.042
"" " " 800043
Y{NDIF*l) * XPQS+1609.3 800044
Y(NDIF+2) » YPOS*1609.3 800045
BD0046
CHIECK FOR SECOND RUN BDOOA7
800048
IF (IN .EQ. 1) GO TO 300 8000^9
B00051
IMPUT METEOROLOGICAL DATA. BD0052
800053
-------
XXI. 2
G LEVEL 20
MAIN
DATE - 72342 21/45/45
C
C
C
C
C
C
00 100 J-ifNSTA . .. .
READ (5*1005) ID-NUIK J)» ((A1R< I* J»K)» K«l»2)» I«l» 4)
READ (5»IOQ1) IDNUM(J),((AIR(I,J^) »K«lj_2) ,J-5* 8)
100 READ (5»lb05~) 1DNUM( J)» UAIRU* J*k),K»i>2)» I«9,12)
INPUT METEOROLOGICAL STATION DATA
WRIT!
READ ( 5» 1014JL.UD.MfT.U > * (WPOJJ I*J
CONVERT DATA TO METRIC UNITS
-.UJLLtl2JUN.ST * >
80.0054
"BE 0055
BE0056
BE0057
80.005.8
800*0*5-9
800060
"800061"
BE0062
"800063
B00064
B00065
BD0066
_DO 120 I«1>NSTA .
120 WPOS(I»3) - WPOS<1*3)*.3048 BD0067
C _____ _ _•_ •_ BD0068
C READ POLLUTANT STATIONi'~POSITIONS" " "" : " BD0069~"~
C _ _ _BJ>0070
READ C5»1014) ( IDPOL( r)»(>POS( I» J>»J«i»"3)"*I»r»NPOL) BE0071
C BD0072
C" READ POLLUTANT DATA ~" BD0073
C _ _ _ _ BD0074
00 140 J*l»~NPbL " "" ~ " ' BD0075
READ <5»_1?0^)_ I_DJ J»»(tCONC(L» J» K»_»K"i»4l_,»i-l.?.3J ______ _BE_00_76_
READ (5»1006) " ID(J)»(tCONC(I,J,K)»K*1*4)»I»4,6) BE0077
RE AiD_l 5» 1006) ID(J JjrJ J CONC ( I» J> K )li. „_ _BD0093_
C """" ' •" ' ~ ~ " •" " BD0094
C VRITE POLLUT^NT__pATA BD0095
.„... ... ._ - B00096
WRITE (6.100J) __ __ _ BD0097
WRITE (6»"10'16) (10(1 J»PPOS(I*l)»PPOSriV2)»I-iVNPOL) "- " " CB0098
K3UT « 0 _ __ _BP_9 \P_P_
WRITE (6,"loTbl ~ ": ' ' " - ~ - BOO'lO'f
00 15_0__ IjLijt_l? BD0102
IH « ( I*6)*100 "" ~ " " " 800103
00 150 J»1>NPOL BD0104
' IF'lKOUf .LE. '501" GO TO 145 " '"."" " BD0105
KOUT = 0 BD0106
-------
XII. 3
G LEVEL 20
MAIN
DATE = 72342
21/45/45
WRITE (6,1010)
G3 TO 150
145 K3UT = KOUT + 1
150 WRITE (6.1011) I0,IH,(CONC(I,J»K),K=1,4)
: SAVE WIND DATA
0!) 200 J=l,12
DO 200 I=1.NSTA
DAIRU. I.I ) = AIR(J.I.2)
200 DAIRU, I ,2) = AIPU, I ,1)
I CONVERT WIND DATA TO RADIANS
DD 230 J=l,12
DO 230 I=1,NSTA
AIRU.I.1) = -OAIRU, I,l)*SIN(DAIR(J, I,2)*3. 14159/180. )
230 AIR(J,I,2) = -DAIRU, 1,1 )*COS(OAIR(J, I,21*3. 14159/130. >
300 WRITE (6,1025) PRMU2). PRMT(3),PRMT(4).GOTIME,A
CALL SMGDRV (PRMT.Y»DERY»2,2»I MLF,RWIND,PUTOUT)
CALL POLLUT (T[ME.COORD(1),COORD<2>,CONC,NPOL,N02,NO*03,CO,PPOS)
WRITE (6,1012) N32.NO.a3,CO
WRITE (6,1020) IHLF
G!3 TO 1
POO STOP
1000 FORMAT(20A4)
1001 F1RMAT(4F10.3.3F10.5)
1002 FORMAT(4HO)
1005 F3RMAT(A3.3X>4(F3.0»F3.1,11X))
1006 FDRMAT(A3»3X,3(3F
-------
XII. 4
G LEVEL 20
SMGDRV
DATE = 72342
21/45/45
C
C
C
C
C
C
C
C
SUBROUTINE SMGORV (PRMT,Y,DERY,NO I!% NSPEC,IHLF.SMGFCT,SMGOUT)
DIMENSION PRMT(1),Y(1),DERY(1),AUX<16,50)
DOUBLE PRECISION AUX»C»0»DELT,QERY,H,X.Y,YY>It II
N»l
IHLF=0
X=PRMT<1)
OELT = 0.000
H=PRMT<3)
ISTE" = 0
DO 1 I=1»NDIN
AJX(16,1 )=O.DO
AUX(15»I ) = DERY( I )
1 AUXd, I)=Y( I )
IF(H*(PRMT(2)-X) )3»2»4
E^ROR RETURNS
2 IHLF=12
GOTO 4
3 IHLF=13
COMPUTATION OF DERY FDR STARTING VALUES
4 CALL SMGFCT (XtY»DERY»-1 )
RECORDING OF STARTING VALUES
CALL SMGOUT (X,Y,OERY,OELT,H»IHLF)
5 IF ( IHLF)7,7,6
6 RETURN
7 DO 8 I=1»NDIM
8 AUX(8, I )=DERY( I)
COMPUTATION OF AUXJ2. I)
ISW=1
GOTO 100
9 X=X+H
03 10 I=1.NOIM
10 AUX(2, I)=Y(I )
INCREMENT H IS TESTED BY MEANS OF BISECTION
11 IHLF=IHLF*1
X=X-H
DD 12 I=1.NDIM
12 AUX(4, I)»AUX(2»I)
H=.500*H
N=l
ISW = 2
GOTO 100
13 X=X*H
CALL SMGFCT (X,Y,DERY,l)
N=2
8D0156
800159
BD0160
BD0161
BD0162
300163
3D0164
3D0165
BD0166
BD0167
BD0168
8D0169
800170
BD0171
BD0172
8D0173
BD0174
BD0175
BOO 176
BD0177
BD0178
BD0179
BDOlflO
BD0161
BD0162
BD0183
8E0184
3D0165
800136
BD0167
BD0183
B00199
BD0190
3D0191
BD0192
BD0193
B00194
BD0195
B00196
BD0197
BD0198
BD0199
B00200
800201
300202
8D0203
3D0204
800205
800206
BD0207
800208
800209
BD0210
-------
XII. 5
G LEVEL 20
SMGDRV
DATE
72342
C
C
.03 14 I»_1»NDIM
AUX (2, I)=Y(I )
14 AUX(9,I)=DERY( I )
ISW=3
GOTO 100
COMPUTATION OF TEST VALUE DELT
SET DELT TO MAX PREDICTOR-CORRECTOR DIFFERENCE
15 OELT = DABS(Y< 1)-AUX(4,1)>
03 16 I=2,NSPEC
IF (DA8S(Y( I )-AUX(4, I ) ) .GT. DELT) DELT = DABS ( Y( I ) -AUX ( 4 , I ) )
16 CONTINUE
IF (DELT-PRNT14) ) 19, 19,17
17 IF ( IHLF-10) 11* lfl, Id
NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.
18 IHLF=11
21
GOTO 4
THERE IS SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS.
19 X=X+H
DO 20 I=l,NDtM
AJX(3f I )=Y( I )
20 AJX( 10,1 ) = DERY( I )
ISW =
G3TO
100
N=l
X=X4H
CALL SMGFCT (X.Y,OERY,1)
X=PRMT(l )
DD 22 IM.NDIM
AUX( 11,1 )=DERY( I)
220Y(I)'»AUX(l»I)+H*(.37500*AUX('8,I
1-. 20833333 33 333 333 3DO* AUX ( 10, I )
23 X=X+H
•»•. 79 16666666666667DO*AUX ( 9» I )
. 04 1 66666666666666700* DER Y ( I ) )
CALL SMGFCT ( X , Y,DER Y , -I )
CALL SMGOUT ( X, Y, DER Y, DE LT, H» I HL F )
24 IF (N-4) 25.201*201
25 03 26 IM.NDIM
AUX(N, I)=Y( I )
26 AUXCN+7, I I=DERY( I )
IF (N-3) 27,29,201
27 03 28 I=1»NDIM
DELT=AUX(9,I) *AUX(9, I )
DELT=DELT+DELT
28 Y( I )=AUX(1,I )•»-. 33333333333333333DO*H*< AUX(d, I ) *DEL T+ AUX( 10 > I ) )
GOTO 23
BD0211
B00212
B00213
BO 0214
BD0215
BD0216
BD0217
3D0218
3D0219
3D0220
B00221
8D0222
3D0223
300224
BD0225
BD0226
BD0227
800228
B00229
3D0230
BD0231
800232
3D0233
300234
BD0235
BD0236
BD0237
3D0238
BD0239
BD0240
B00241
300242
B00243
8D0244
800245
8D0246
BD0247
BD0243
300249
BD0250
BE0251
300252
B00253
BD0254
8002.55
800256
B00257
BD0258
BD0259
BD0260
800261
B00262
BD0263
-------
XII. 6
G LEVEL 20
SMGDRV
DATE = 72342
21/45/45
C
C
C
C
C
C
C
29 DD 30 1 = 1, NO IN . .... .
OELT*AUX(9»I)+AUXMO» I)
DELT=DELT+DELT+DELT
30 Y< I)*AUX(1,I )+.375DO*H*< AUX(8>I)+DELT+AUX( 11,1 M
GOTO 23
THE FOLLOWING PART OF SUBROUTINE OHPCG COMPUTES BY MEANS OF
RJNGE-KUTTA METHOD START I NG VALUE S FOR THE NOT SELF-STARTING
PREDICTOR-CORRECTOR METHOD.
100 00 101 I=1,NDIM
Z=H*AUX(N+7, I >
AUX(5»II«Z
101 Y< I )=AUX(N,I )+.4DO*Z
Z IS AN AUXILIARY STORAGE LOCATION
Z=X+.4DO*H
CALL SMGFCT < Z, Y, DER Y, 1 )
DO 102 I»1.NDIM
Z=H*DERY(I)
AUX(6,I)=Z
102 'Y( I )=AUX(N, I )«-. 2969776092477 5360 00* AUX( 5, I
. 1 5875964497 1035 83DO*
Z=X*. 4557372 54 21 8789 4300*H
CALL SMGFCT ( Z > Y »DE3 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS
202 DO 203 N=2»7
03 203 I=1,NDIM
AUX(N-1» I)»AUX(N.»I)
800304
800305
B00306
800307
800308
8D0309
800310
800311
8D0312
B00313
BD0314
800315
800316
-------
XII. 7
G LEVEL 20
SMGDRV
DATE = 72342
.21/45/45
203 AUXIN+6, I)=AUX(N+7,I )
C
C
C
C
204
N LESS
N=N+1
THAN 8 CAUSES NG L TO GET N
COMPUTATION OF NEXT
DO 205 I=1,NDIM
AUXIN-l, I)=Y( I )
AUX(N*6, I)=DERY(I)
VECTOR Y
205
206 DO 207 I=1,NDIM
ODELT=AUX(N-4,I)+ 1.3333333333333333DO*H*(AUX(N+6,1)+AUX(N+6,I)-
IAUXCN + 5.. Il+AUX (N + 4, I ) +AUX(N+4» I) )
Y( I>«DELT-.9256 19834710743600*AUX(16,I )
207 AJX( 16,1 )=OELT
PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR
IS GENERATED IN Y. DELT MEANS AN AUXILIARY STOPsAGE.
C
C
C
C
C
CALL SMGFCT -AUX(N-3,I)+ 3.DO*H*(DERf(I )+AUX(N+6,
UAUX(N+6,I )-AUX(N+5, I ) ) )
AUX< 16,1 >=AUX< 16. D-DELT
208 Y( I ) =DELT+.074380 1652892562000*AUX(16,1)
TEST WHETHER H MUST BE HALVED OR DOUBLED
SET DELT TO MAX PREDICTOR-CORRECTOR DIFFERENCE
DELT = DABS(AUX(16,1>)
DO 209 I=2,NSPEC
1= (DABS(AUX(16,1)) .GT. DELT)
209 CONTINUE
IF (DELT .GE. PRMT(4)) GO TO 222
IHLF = 0
H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD.
210 CALL SMGOUT (X,Y»DERY,DELT,H,IHLF )
211 IF(IHLF-11)213,212,212
212 RETURN
DELT = DABS(AUX(16, I)J
CHECK IF END TIME REACHED.
213 IF(H*'(X-PRMT(2 )-0. ID D) 214.250,250
214 IF(DABS(X-PRMT(2)-0.ID 1)-.100*0ABS(H))
215 IF (DELT .LT. . 01*PRMT (4)") GO TO 217
216 ISTEP = 0
GO TO 201
250,215,215
H COULD BE DOUBLED IF ALL NECESSARY
AVAILABLE. CHECK FOR a CONSECUTIVE
217 IF USTEP-7) 213,219,219
218 ISTEP = ISTEP+l
PRECEEDING VALUES ARE
DELTS WITHIN LIMIT.
8D0317
BD0318
300319
BD0320
BD0321
BD0322
800323
' 8D0324
BD0325
BD0326
' 800327
8D0328
300329
BD0330
800331
8D0332
B00333
800334
BD0335
BD0336
300337
8D0338
300339
I JBD0340
B00341
300342
BD0343
BD0344
BD0345
SD0346
3D0347
BD0348
600349
B00350
300351
BD0352
3D0353
B00354
BE0355
ED0356
B00357
3D0358
.300359
BD0360
BD0361
BD0362
BD0363
3D0364
BD0365
BD0366
BD0367
B00366
BD0369
-------
XII. 8
G LEVEL 20 SMGDRV DATE = 723.42 21/45/45
G3 TO 201 BD0370
219 IF (DABS(H)-17.0) 220.216*216 ' ~BD0371
220 H=H+H BD0372
ISTEP=0 BD0373
00 221 I=1*NDIM B00374
AUXIN-l,I
AUXCN-2,I
AUXtN-3.I
AUX(N+6» I
AUX(N+5»I
AUX(N+4, I
=AUX(N-2»I) BD0375
=AUXJN-4,I) BD0376
=AUX(N-6,I) " B00377
*AUX(N+5»I) B00378
=AUX(N+3,I) B00379
*AUX(N+1»I) B00380
DELT=AUX(N+6*I)+AUX(N+5»I) BD0381
DELT=DELTfDELT+OELT 800382
2210AUXI16»I >=8.962962962962963DO*(Y(I)-AUXIN-3* I ) ) BD0383
1-3.361llllll1111111DO*H*(DERY(I)+DELT+AUXI)* BD0396
1108.00*AUX(N-3> I)+AUX(N-4.. I ) )-. 023437 500* ( AUXOELT*DELT BD0412
OAUX( 16,1 ) = 8.9629629629&2963DO*UUX
-------
XII. 9
G LEVEL 20
PUTOUT
DATE = 72342
21/45/45
SUBROUTINE PUTOUT .( T I M5 » Y» DER Y .D ELTA, 3 I ZE, I HL F >
DIMENSION Y(I)»DERY( 1)
DOUBLE PRECISION DELTA.DERY.SIZE*TINE,Y
CDMMON /OUTCOM/ TOUT,TINT
COMMON /MILES/ COORD(2 ) ,OXY(2 )
IF ( IHLF .NE. 30) GO TO 50
GD TO 100
50 IF (OABS(TIME) .GE. TOUT) GO TO 100
GD TO 2000
100 WRITE ' (6.1000) TI ME,COORDl1)»COORD(2»,SIZE
TOUT = TOUT+TINT
1000 FORMAT(«OTIME=t»Dl5.8»l5X.'POSITI ON*•,2PI 5,
2000 RETURN
END
15X,«SIZE=',015.a)
BE0421
BD0422
BD0423
800424
300425
BE0426
BD0428
BD0429
BD0430
CB0431
B00432
CB0433
BD0436
BD0437
-------
XII.10
G LEVEL 20
RWINO
DATE = 72342
21/45/45
SUBROUTINE Rl«l I N0( XT I ME. Y» OER Y, NOT )
C
c
C
10
15
C
c
c
c
c
c
DIMENSION Y(1),DERY(1),0 I STl32)»JBAR(32)»VAL1(9)»VAL2(8).
1 DISINV(32)
DOUBLE PRECISION DERY,XTIME»Y
COMMON /MILES/ COOR 0 ( ?.) , DXY ( 2 )
COMMON /RECOM/ TI ME, ANS<8),DERI V<8 ) , NSTA,NPOL>AIR(12,32*2),
1 WPQS<32,3) .GOTIME
SET TIME AND POSITION
TIME = GQTIME+XTIME/60.
ITIME = TIME-6
CDOROt 1) * Yd )/1609. 3
COORD(2) = Y(2)/1609.3
DO 10 1=1,32
DISINV(I) = 0.
DIST»COORD(2)»WPOSfJBAR,NSTA)
DO 30 1=1,NSTA
IF (JBARU ) .NE. 1) GO TO 30
OIST(I) = (COORD(1)-WPOS(I,l))**2 * (COORD(2)-WPOS(I,2))**2
IF (OIST(I) .NE. 0) GO TO 30
IDNE = I
GO TO 100
CONTINUE
SET XDOT AND YDOT INTERPOLATED BETWEEN TrfO SUCtSSIVE HOURS
DO 50 J=l»2
OD ^0 1=1,NSTA
IF (JBAR(I> .NE. 1) GO TO 40
IF (AIR( ITIME + 1, I,J) .EO. 0.) GO TO 40
DISINVU) = DISINVt J)+l. /DISTC I >
VAL2(J) = VAL2(J)+AIR(ITIME+1,I»J)/DIST(I)
CDNTINUE
IF (DlSINV(J) .NE. 01 GO TO 45" " " "
VAL2(J) = 0.
GO TO 50
45 VAL2U) = VAL2( J)/DISINV( J)
50 CONTINUE
00 60 1=1,2
60 DISINVtI) = 0.
DO 80 J*l,2
DO 70 1=1,NSTA
IF (JBARU ) .NE. 1) GO TO 70
30
40
BD0438
BD0439
CA0440
CA044010
BD0441
BD0442
8D0443
BD0444
B00445
3D0446
BD0447
BD0448
800449
BD0450
BD0451
BD0452
B00453
BD0454
BD0455
CA045520
CA045540
BD0456
BD0457
BD0458
BD0459
BD0460
BD0461
BD0462
8D0463
5D0464
BD0465
CA0466
BD0467
3D0468
CA0469
B0047Q
BD0471
300472
BD0473
CA0474
BD0475
CA0476
BD0477
BD0473
CA0479
3D0480
CA0481
BD0482
CA0433
CA048320
CA046340
CA048360
CA048380
-------
XII.11
G LEVEL 20 RWIND DATE » 72342 21/45/45
IF (AIRUTIME, I.J) .EQ. 0.) GO TO170 CA0484
OISINV(J) = DI~S7NV< j'F+~l."/GIST! I f~" " " ~ ~ CA048420
VAL1U) = VALK J)+AIR( ITIME» I. J)/DIS_T( I ) CA0^8440
70 CONTINUE " """ " CA048460
IF (OISINV(J) .NE. 0.) GO TO 75 __ CAQASAflO
VALKJ) =0. CAO<»85
GO TO 80 CA048520
75 VALKJ)"- VALK J)/Ol"SINV(J) . """ " CA048540
80 CONTINUE CAOA8560
TIMEHR » TIHE-( ITIME + 6) - CA048580
00 90 J=l»2 CA0436
90 AMSU) " VALK J )-M VAL 2 ( J >-VAL K J ) )*TIMEHR CA048620
GO TO 120 CA04S640
C CA048660
C OMLY ONE STATION CA048670
C CA048680
100 OD 110 J*l,2 • CA0487
110 ANS(J) = AIR(ITIME»IONE»J) CA0438
120 00 130 I=l»2 CA0489
OXY(I) » ANS(I)/60. - 300A90
130 OERY(I) = DXY« I )*1609.3 CAOA91
RETURN 300492
END 300493
-------
XII.12
G LEVEL 20
POLLUT
DATE a 72342
21/45/45
C
C
C
C
C
C
SUBROUTINE POLLUT (T I ME, X, Y,CO_NC*NPaL, A l»_A2» A3, A4, PPOS ) (
CALCULATE POLLUTANT STARTING VALUES
DIMENSION CONC<12*20,5),VAL<5)»DISINV<32),DIST<32),JBAR(32),
1 PPOS<32,3)
ITIME = TIME-6
00 10 I-l.NPOL
DISINVlI) = 0.
10 DISTd ) = 0.
DO 15 1=1,4
15 VAL(I) = 0.
DETERMINE USABLE STATIONS AND SET STATION DISTANCES
30
(Y-PPOSJ 1,2) )**2
30
CALL BARIER (X,Y,PPOS,JBAR,NPOL)
DO 30 I=1,NPOL
IF (JBAR(I) .NE. 1) GO TO
DISTd) = +CONC( ITIME,!, J»/DIST( I )
40 CONTINUE
IF (OISINVU) .NE. 0) GO TO 45
VAL(J) = 0.
GD TO 50
45 VAL( J) » VAL( J I/DISINVU)
50 CONTINUE
GO TO 80
60 DO 70 J=l,4
70 VAL(J) = CONC( ITIME,IONE,J)
80 Al = VALd )
A2 = VAL«2)
A3 = VAL(3)
A4 = VAL(4)
RETURN
EMD
BE0494
BD0495
8D0496
3D0497
BE0498
B00499
3D0500
BD0501
BD0502
BOO503
8E0504
BD0505
BD0506
BD0507
BD0508
BD0509
BD0510
BD0511
BD0512
BD0513
BD0514
BD0515
BD0516
BD0517
BD0518
B00519
BE0520
BD0521
BD0522
BD0523
3D0524
BD0525
BD0526
B00527
B00528
BD0529
B00530
BD0531
BD0532
BE 0533
BD0534
BD0535
BOO536
3D 0537
BD0533
BD0539
BD0540
-------
XII.13
G LEVEL 20
3ARIER
DATE = 72342
21M5/45
C
C
C
C
SUBROUTINE 8ARIER (X»Y,WPOS,JBAR,NPOS)
_ BD0541
300542
DETERMINES IF INTERSECTION OF LINE SEGMENT CONNECTING TRAJECTORY 8005A3
PDINT AND STATION AND THE BARRIER LINE SEGMENT LIES ON THE BARRIERBD0544
ITSELF. EVALUATE FUNCTIONS ALPHA AND BETA. IF THEY ARE BOTH LESSBD0545
C
C.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
50
100
C
C
C
200
C
C
C
C
THAN 1 AND GREATER THAN 0. STATION IS ELIMINATED.
n I ME NS~ION WPOS ( 32, 3 ) , JBAR < 32 ) /BAR ( 10, 10 ) , BSLOPE ( F6", 1 0) ,"
I TSLOPE(32,2)
NBAR IS NUMBER OF BARRIERS
A*RAY BAR ODD COLUMNS CONTAIN X COORDINATES
EVEN COLUMNS CONTAIN Y COORDINATES
ROW 1 IN ODD COLUMNS"!: ONTA IN NO. "OF "POINTS IN BARRIER
ARRAY BSLOPE--ODD COLUMNS CONTAIN M VALUES OF BARRIER LINE SEGMTS
EVEN COLUMNS CONTAIN Y VALUES OF BARRIER LINE SEGMENTS
ARRAY TSLOPE--ODD COLUMN CONTAINS M VALUES OF TRAJECTORY POINT
WEATHER STATION LINE SEGMENT
EVEN COLUMN CONTAINS B VALUES OF TRAJECTORY POINT
WEATHER STATI-ON L INE SEGMENT
COMMON /BARRR/ IDB
DATA NBAR/ If, BAR/2 . » 16. » 38. , 3*0. ,2*40. » 7*0. , 4. » 5.3 . ,60. ,66. »71. ,
G 6*0. , 31. ,23. ,32. ,23., 65*0. /
"'DETERMINE SLOPE AND INTERCEPT FROM INTERPOLATION' POINT TO EACH"
STATION.
DO 50 I=1,NPOS
TSLOPEU,!) = 0.
TSLOPECI,2) = 0.
DD 100 I=1.NPOS
IFUWPOS(I,l)-X) .EO. 0.) GOTO 100
TSLOPEd,!) =« (WPOS( I.2)-Y)/(WPOS( I,l)-X)
TSLOPECI,2) = Y-TSLOPEII,1)*X
CONTINUE
DETERMINE SLOPE AND~ INTERCEPT OF~ ~E AC H SEGMENT OF EACH BARRIER
IF ( IDB .EQ. 1 ) GO TO 210
00 200 J=1,NBAR
NSBAR = 8AR( 1,2*J-1)
DO 200 K=2,NNBAR
IF ( (BAR(K + 1»2*J-1>-BAR(K,2*J-1) ) .EQ. 0.) GO TO 200
8SLOPE(K,2*J-1 ) = (BAR(K + 1,2*J)-BAR(K,2*J) )/
1 (BAR(K+1*2*J-1)-BAR(K,2*J-1) )
BSLOPE(K,2*J) = BAR(K*2*J)-BSLOPE(K,2*J-1)*8AR(K,2*J-1)
CONTINUE
IOB = 1
DETERMINE INTERSECTION FUNCTIONS
IF VALUE OF INSECTION IS LE TO ZERO OR GE 1,
INCLUDE STATION.
BD05A6
B00547
BD05A8
3005^9
BD0550
BD0551
BD0552
BD0553
3D0554
BD0555
BD0556
BD0557
BD0558
BD0559
BD0560
BD0561
8D0562
8D0563
BD0564
BD0565
BOO 566
8D0567
BD0568
BD0569
BD0570
600 571
BOO 572
BD0573
BD0574
BD0575
B00576
BD0577
BD0578
BD0579
BO 05 80
8D0581
BD0582
BD0583
BD058A
BD0585
BD0586
BD0587
BD0588
800589
BD0590
BD0591
BD0592
BD0593
-------
XII.14
G LEVEL 20
BARIER
DATE
72342
21/45/45
210
220
240
260
00 400
00 375
NNBAR '
DO 375
I=l,NPOS
J=l,NBAR
BAR!1»2*J-1)
K«2,NNBAR
IF (
------- |