U.S. DEPARTMENT OF COMMERCE
National Technical Information Service
PB-270 778
STRAM - An Air Pollution Model
Incorporating Nonlinear Chemistry,
Variable Trajectories and Plume
Segment Diffusion
Battelle Pacific Northwest Labs, Richland, Wash
Prepared for
Environmental Protection Agency, Research Triangle Park, N C Office of Air
Quality Planning and Standards
Apr 77
-------
PB 270 778
EPA-450/3-77-012
April 1977
STRAM
AN AIR POLLUTION MODE
INCORPORATING NONLINEA]
CHEMISTRY, VARIABL
TRAJECTORIES, AND PLUM]
SEGMENT DIFFUSIO1
U.S. ENVIRONMENTAL PROTECTION AGENCY
Office of Air and Waste Management
Office of Air -Quality Planning and Standards
Research Triangle Park, North Carolina 27711
l TECHNICAL
INFORMATION SERVICE
U.S. DEPARTMENT OF COMMERCE
SPRINGFIELD, VA. 22161
-------
NOTICE
THIS DOCUMENT HAS BEEN REPRODUCED
FROM THE BEST COPY FURNISHED US BY
THE SPONSORING AGENCY. ALTHOUGH IT
IS RECOGNIZED THAT CERTAIN PORTIONS
ARE ILLEGIBLE, IT IS BEING RELEASED
IN THE INTEREST OF MAKING AVAILABLE
AS MUCH INFORMATION AS POSSIBLE.
-------
STRAM - AN AIR POLLUTION MODEL
INCORPORATING NONLINEAR CHEMISTRY,
VARIABLE TRAJECTORIES,
AND PLUME SEGMENT DIFFUSION
by
J.M. Hales, D.C. Powell, and T.D. Fox
Battelle Pacific Northwest Laboratories
P.O. Box 99
Richland, Washington 99352
Contract No. 68-02-1982
Program Element No. 2AC129
EPA Project Officer: Joseph A. Tikvart
Prepared for
ENVIRONMENTAL PROTECTION AGENCY
Office of Air and Waste Management
Office of Air Quality Planning and Standards
Research Triangle Park, North Carolina 27711
April 1977
-------
This report is issued by the Environmental Protection Agency to report
technical data of interest to a limited number of readers. Copies are
available free of charge to Federal employees, current contractors and
grantees, and nonprofit organizations - in limited quantities - from the
Library Services Office (MD-35), Research Triangle Park, North Carolina
27711; or, for a fee, from the National Technical Information Service,
5285 Port Royal Road, Springfield, Virginia 22161.
This report was furnished to the Environmental Protection Agency by
Battelle Pacific Northwest Laboratories, P.O. Box 99, Richland, Washington
99352, in fulfillment of Contract No. 68-02-1982. The contents of this report
are reproduced herein as received from Battelle Pacific Northwest Labora-
tories . The opinions, findings, and conclusions expressed are those of
the author and not necessarily those of the Environmental Protection Agency.
Mention of company or product names is not to be considered as an endorse-
ment by the Environmental Protection Agency.
Publication No. EPA-4fft)/3-77-012
11
-------
ABSTRACT
This document provides a technical description, user's guide and program
listing for (1) STRAM - a variable trajectory, reactive plume-segment model
for ground level air pollution assessments resulting from multi-source
emissions on a multi-state scale, and (2) a supporting program, Random-to-
Grid, which generates gridded wind data for STRAM from synoptic wind data
at arbitrarily located observing stations.
The advection is defined over a Cartesian grid of which the sampling grid,
also Cartesian, is a subset, usually of greater density. There is also
provision for further sampling at particularly specified receptor locations.
The reactive plume chemistry is calculated by a Subroutine STRAC and
related subroutines, which calculate the diffusion, the wet and dry
depletion, and the reactive chemistry within each plume segment.
The formulations of reactive chemistry, of wet and dry plume depletion, of
plume rise, and of diffusion will differ for separate applications and can
be expected to change as the states of knowledge in each of these areas
evolve. Therefore, STRAM isolates each of these calculations in individual
subroutines so that the user can easily effect timely modifications of the
code as is deemed appropriate.
STRAM also has provision for time-varying definitions of stability
and mixing height. However, the manner of specification and of entry
of the data is left to the discretion of the user, notwithstanding the
included example.
The principal output of STRAM is concentrations on the sampling grid and
at each particularly specified sampling point for each of the analyzed
chemical components. These are available for three averaging periods
(1) once for the entire running time, (2) serially for the basic sampling
interval, and (3) serially for an arbitrarily specified intermediate time.
Matrices of maximum values over all matrices of this last type ara also
printed out.
Output is made to magnetic tape of all matrices averaged over the basic
sampling interval and of such input parameterization as is essential to
interpret the results of the run.
iii
-------
CONTENTS
ABSTRACT iii
FIGURES vi
TABLES vii
ACKNOWLEDGMENTS viii
I. INTRODUCTION 1
II. TECHNICAL DISCUSSION 3
A. Plume Chemistry and Dispersion 3
1. Theory of Basic Equations 3
2. Finite Difference Approximation 8
3. Examples of Applicability 11
B. Trajectory Analysis 12
1. Grid Definitions 12
2. Theory of Basic Equations 15
3. Approximation of the Wind Integral by Finite
Summation 16
4. Examples of Application 19
C. Calculation of Input Parameters 19
1. Preparation of Data for Calculation of Trajectories. . 19
2, Calculations of Plume Chemistry and Dispersion
Data 19
III. USER'S GUIDE 25
A. INTRODUCTION 25
1. General Data Processing Requirements 25
2. STRAM System Flow Chart 25
IV
-------
CONTENTS (Cont.)
B. RANDOM-TO-GRID 25
1. Description, Function, Constraints 25
2. Flow Diagram - Major Functions 27
3. Executive Control Language and Deck Setup 29
4. Input Description and Formats 30
5. Output Description and Formats 35
6. Diagnostic Messages 36
C. STRAM PROGRAM 36
1. Description, Functions and Constraints 36
2. Flow Diagram - Major Functions m. . 52
3. Executive Control Language and Deck Setup 52
4. Input Description and Formats 54
5. Output Description and Formats 60
6. Diagnostics 62
IV. INSTALLATION - TEST RUN 65
A. Random-to-Grid Test Run 65
1. Input Data 65
2. Output Data 71
B. STRAM Test Run 91
1. Input Data 91
2. Output Data 95
V. SOURCE PROGRAM LISTINGS 106
VI. REFERENCES 140
VII. GLOSSARY 141
VIII. ADDENDA 147
-------
FIGURES
1 Plume Increments and Segments From Two Sources 1
2 Schematic Illustration of Mass Balance Principle 3
3a Portion of Truncated 47 x 47 NMC Grid Used in
Establishing Coordinate System for STRAM. The Origin
in This System Has Been Shifted to the North Pole 14
3b Expanded View of the Box in Figure 3a 14
4 Theoretical Trajectories Step. ... 15
5 Two-step Approximation of Trajectory Step , 17
6 Washout Coefficients for Size Codistributed Aercsols 21
7 System Flow Chart , 26
8 Random-to-Grid Flow Diagram - Major Functions 28
9 Truncated NMC Grid, Used as a Basis for Interpolation
Grid Positioning on Input Card #5 to RNGRD 33
10 Flow Diagram for Subroutine STRAC 43
11 Flow Diagram for Subroutine INTGRT 46
12 Schematic Illustration of Plume as Seen in Subroutine
SAMPLE 50
13 Flow Diagram of the Major Functions of STRAM 53
VI
-------
TABLES
Forms of the Washout and Dry Deposition "Generation"
Integrals in Equation 9
2 Base Points and Weighting Factors for Ten-Point
Gauss-Laguerre Integration 8
3 Suggested Reaction Rate Parameters 20
4 Coefficients for Dispersion Parameters Formulae 24
5 Card Input Descriptions 30
6 Wind Data Tape Input Description 34
7 Card Input Description for STRAM 55
8 Card Input for Test Run of Random-to-Grid 65
9 Card Input for Test Run of STRAM 91
vii
-------
ACKNOWLEDGMENTS
The wind-field portion of the plume model described in this report is an
expansion of an original scheme developed by Dr. L. L. Wendell of Battelle-
Northwest. The authors express their sincere appreciation to Dr. Wendell
for his helpful guidance and suggestions given during the course of this
project.
Vlll
-------
SECTION I
INTRODUCTION
STKAM (Source-Transport-Receptor Analysis Model) is a user-oriented computer
program that provides computed estimates of concentrations and removal rates
in airborne, multicomponent, reactive plumes. A multiple source capability
is provided, thus allowing calculations, at points on a regional-scale map,
of ground-level concentrations which have arisen from a superimposition of
a number of plumes originating at different geographical locations.
At the outset it is important to observe that STRAM provides calculated plume
distributions that pertain to a given instant of time. They would correspond,
for example, to a snapshot taken from a satellite of plumes on the earth's
surface below. A visual portrayal of such computed distributions is shown
in Figure 1, which (for this example) is a view of the plumes divided into
plume segments.
O PLUME INCREMENT
PLUME SEGMENT
Figure 1. Plume increments and segments from two sources
-------
Such distributions can be computed for successive increments of time, and a
time series of these distributions can be averaged to obtain climatological
pollutant concentration distributions.
The workings of STRAM consist of three major routines; 1) a trajectory
routine that calculates trajectories of plume increments over the regional-
scale area using historical gridded wind data, 2) a dispersion routine that
calculates dispersion parameters, wet and dry removal, and transformation
of mass within each plume by chemical processes, 3) a sampling routine that
interpolates plume concentrations to grid intersections and other speci-
fied locations.
The output of STRAM to line printer and magnetic tape consists of concentra-
tions averaged over specified time intervals and certain maxima over these
time intervals.
-------
SECTION II
TECHNICAL DISCUSSION
A. PLUME CHEMISTRY AND DISPERSION
1. Theory of Basic Equations
The fundamental basis for the STRAM plume model is an integral material
balance over a finite volume element of plume, as indicated in Figure 2.
Denoting pollutant concentration by c and the wind speed in the downwind
direction by u, one can represent the fluxes of pollutant into and from the
represented volume element by
uc dy dz
and
00 00
uc dy dz
x+<5x
respectively.
X+6X
REMOVAL
Figure 2. Schematic illustration of mass balance principle
-------
The rate of gain* of pollutant in the volume element by reaction and deposi-
tion processes will be denoted here by the term G**. Although a number of
mechanisms can act to deplete or generate material, it is appropriate at
this point simply to lump them into a single term, which can be disaggregated
later as desired. The rate of gain of material by reaction, washout, and
dry deposition can be described by the term
ff
•^ ''—CO
G(x,y,z) dy dz 6x
Since the time rate of change of pollutant abundance in the volume element
is equal to the sum of contributions by the above mechanisms, one can perform.
a fundamental material balance to obtain the following expression:
c dy dz 6x =» time rate of change of pollutant
mass in macroscopic volume element
i-f f
dtJ0 L
• L L -**,-/ /
o -°°
n
•* r\ "£_OO
—oo /• oo
uc dy dz
G(x,y,z) dy dz 6x . (1)
Equation (1) is directly extendable to multicomponent systems simply by sub-
scripting c and G (c., i=A,B,C—N; G., i»A,B,C...N where i is a pollutant
species index), and writing N simultaneous equations of the form (1) for
the n components existing in the plume.
One should note that equation (1) is highly general, its primary precepts
being those of mass conservation and differential calculus.
We now shall proceed to limit this generality by imposing a number of stipula-
tions / which are itemized as follows:
• A quasi-steady state exists, where the left-hand side of (1) can
be set equal to zero, at least for short intervals of time.
*Note that loss is a negative "gain." The sign convention will be used
that denotes gain of material as a positive entity.
**G has units of mass/length -time.
-------
• The shapes of the concentration profiles are known. We will
consider two cases at present; these are where the concentration
profiles follow (1) a bivariate-nonnal form, given by
exp
(z+h)
20 2
z -i
(2)
and (2) a limited mixing height form, given by
f(x)
1/2TT U HO
exp -
(3)
20
• The wind profile is uniform vertically, and does not change appre-
ciably in the interval (x, x+6x).
Applying the above assumptions and taking limits in x in the conventional
manner, one obtains the simpler form
dx
• ~
«oo * oo «oo p t
u/ / c dy dz • = / /
JQ J-ao Jo J-°o
G(x,y,z) dy dz.
(4)
Now, since c (or at least its profile shape) is known and integrable, one
can express the bracketted term in equation (4) as an algebraic quantity.
Thus,
dx
G(x,y,z) dy dz
(5)
or, in the multicomponent case,
«i r
="JL L
fo
Q "A
in th
dy dz, i=A,B,C...N
(5a)
Q =/ /o,, u c dy dz may be interpreted physically as the amount of material
in the plume passing a downwind point per unit time. At x=0, or, for plumes
where G=»0, this is simply the source strength. Under other conditions Q may
be visualized as a "virtual" source strength, which may increase or decrease
with distance downwind from the source depending upon whether the pollutant
-------
specie is being added to or depleted from the plume with distance. Under
the assumptions stipulated above, Q is also equal to the f (x) term appearing
in equations (2) and (3). Thus, one can solve for concentration values in
any plume if appropriate values of Q are known. Q in turn may be generated
by solving the integral equations (5a)—a task that provides the basis for
a major portion of the STRAC algorithm.
Solution of (5a) requires mathematical description of G, and a number of
additional modeling assumptions are introduced with this step. These will
be subdivided into those regarding chemical transformation, precipitation
scavenging, and dry deposition in the following text.
Chemical Reaction; The basis for calculating the chemical-reaction component
of G. is simply the conventional rate expression utilized in chemical kinetics,
r (x,y,z) = r.(T,P,hv...c ,c ,c .. .c ) . (6)
1 -L A a (_ W
Here r. is interpreted as the rate of gain of pollutant i by chemical reaction
in a differential volume element of atmosphere. Since the reaction of com-
ponent i may depend upon the concentrations of other species, this term serves
as a coupling entity between equations (5a).
One should note here that concentrations of real plumes fluctuate in time,
and equations (2) and (3) represent averages over these fluctuations. If
equation (6) is nonlinear in concentration, then a positive or negative error
(depending upon the direction of deviation from linearity) will result in
calculation of the transformation term. A fundamental assumption of this
modeling effort is that such errors are negligibly small.
Precipitation Scavenging; The basis for calculating the precipitation scaveng-
ing component of G. is the washout-coefficient concept. Defining the amount
of material removed from a differential volume element of plume per unit time
by
- W^X^Z) =» A^tx^Z) , (7)
where A, is the washout coefficient, provides a direct means for evaluating
this component. One should note here that equation (7) implies the basic
modeling assumption of irreversible scavenging.
Dry Deposition; The basis for calculating dry deposition in the present model
is the deposition-velocity concept. Here the flux of material delivered to
an incremental area of ground surface is
F (x,y,0) = V, c.(x,y,0) , (8)
a. d. i
i i
where V (with units of length/time) is the deposition velocity for plume
constituent i.
-------
The mathematical approach utilized by the STRAC algorithm is to evaluate the
integral terms in (5a) prior to solving the differential equation in x. From
6, 7, and 8, one obtains
r r r r
I I Gi(x,y,z) dy dz - I I r±(x,ytz)
J0 •'-oo X J0 '-co
« r f
/ I " Ai ci(x'y'2) d* dz + I ~Vd Ci
•'n •'-oo J-ta *
dy dz +
(x,y,0) dy
(9)
Since r,(x,y,z) depends upon concentration in a manner that may be rather
complicated, this model utilizes a numerical technique for evaluation of its
integral. The corresponding terms for washout and dry deposition may be
evaluated directly by formal integration in conjunction with equations (2)
and (3). The results of these integrations are given in Table 1.
TABLE 1. FORMS OF THE WASHOUT AND DRY DEPOSITION "GENERATION" INTEGRALS
IN EQUATION 9.
Bivariate
normal
Limited
mixing
height
Washout
Dry Deposition
u
v
exp
A. Q.
i *i
u
u H
-------
2. Finite Difference Approximations
As indicated above, finite differencing techniques are employed by the STRAC
algorithm to approximate the reaction-rate integral in equation (9), and also
to evaluate the x dependencies of Q. in equations (5a). The first of these
applications, i.e., that of evaluating the integral
Cr. a f I ri(x'*'
Jn •'—oo
z) dy dz
(10)
is accomplished in a straightforward manner using a ten-point Gauss-Laguerre
quadrature:*
-5
10
i
(11)
Here the £. are selected values of the independent variable and the w. are
weighting factors, as shown in Table 2. The expected error from this approxi-
mation technique is
5.4 x 10
*6
in
(12)
TABLE 2. BASE POINTS AND WEIGHTING FACTORS FOR TEN-POINT GAUSS-
LAGUERRE INTEGRATION
Independent variable
base points
0
0
1
3
5
8
11
16
21
29
.13779
.72945
.80834
.40143
.55249
.33015
.84378
.27925
.99658
.92069
34705
45494
29017
36978
61400
27467
58379
78313
58119
70122
40
03
40
55
64
64
00
78
81
74
Weighting
w.
i
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
30844
40111
21806
62087
95015
75300
28259
42493
18395
99118
factors
11157
99291
82.876
45609
16975
83885
23349
13984
64823
27219
65
55
12
87
18
88
60
96
98
61
x
x
X
X
X
X
X
— T
10_2
10
10~4
10~£
10
10~?0
— 1 2
10 L
-------
Equation (10) is evaluated using (11) as follows: first, the y and z extents
of the plume are scaled so that the primary portions of the plume profile are
well covered by the independent variable points £. given in Table 2. r.(x,y/2)
is then multiplied by the factor e^ to compensate for its reciprocal term in
(11), and the resulting function is integrated over y at each chosen level of
z (dictated by £.). The ten values of the y integral resulting from this
operation are then doubled (to reflect the fact that the limits for the y
integration are (-00, °°) rather than (o,°°), and are used as base points for
a second integration in the z-direction. In applying (11) for this final
integration the base functions once again are premultiplied by e^ to compen-
sate for its reciprocal in the integral.
Values for y and z are scaled for the numerical quadrature in a somewhat
arbitrary manner to provide good coverage of the plume by the independent
integration variable £, whose base points (note Table 2) range approxi-
mately between zero and thirty. In the y-direction we let
5 - ft , (13,
y
thus,
I (z) * / e * F(y) dy - j*- / e s F(£) d£ . (14)
(p(y) - ey ri(x,y,z)
For the z-direction we let
£ » — + 15 - —
U U
Z Z
unless h < 5a , in which case we let
. (16)
z
These transformations provide the integral representations
r r
r - I e~Z Hy(z) dz = ~ I e"C Hy(?) d?
•' •'
("y ' '* V")
-------
Scaling y and z using (13), (15), and (16), integrating using (11), and then
descaling using (14) and (17) allows the integral (10) to be approximated in
a relatively straight-forward manner.
Once the integral Ir. is evaluated it may be added to the appropriate entries
of Table 1 to obtain the total integral in equation (5a). For a bivariate-
normal plume this is
A. Q
(18)
For the limited mixing-height case, it is given by
ri
u
u H
(19)
This effectively reduces the integro-differential equations (5a) into differen-
tial equations of the form
d x
(x)
(20)
The STRAC algorithm generates finite-difference approximations to solutions
of equation (20) using a second-order Runge-Kutta technique. This technique
utilizes an initial condition
q. at x =* 0
(21)
plus a derivative evaluation at x = 0, to estimate a derivative
dQi
f 3 —
Ax dx
Ax/2
at the half-interval step Ax/2. This value is then utilized with condition
(21) to provide an estimate of Q. at x = Ax:
Ax
+ Ax f
Ax
(22)
Values of 0. at larger x are generated by repeating this procedure. The initial
condition (21) siznply is replaced by the derived value of Q, at the last point
of calculation, and the process is continued until the desired downwind distance
10
-------
is attained. The fractional truncation error* associated with this method
is of the order of
where X is the total downwind distance covered by the integration.
tot
In summary, the numerical procedure utilized to generate solutions to equations
(5a) follows the following sequence:
1. Evaluate the double integral in equation (10) using a double
application of the Gauss-Laguerre quadrature. Repeat for
all components in the system.
2. Integrate differential equation (20) for one step in the
x-direction using the second-order Runge-Kutta algorithm.
Repeat for all components in the system.
3. Repeat steps 1 and 2 until the desired downwind distance has
been attained. Generate ground-level concentration values as
desired using equations (2) and (3) in conjunction with the
computed values of Q..
3. Examples of Applicability
•
The plume model described in the previous text is applicable to a wide
variety of situations. In addition to being useful as a single-source model
for predictive and diagnostic purposes, it can be utilized as a multiple source
descriptor. When combined with descriptions of wind fields and other mete-
orological parameters, this model can be utilized to provide concentrations,
at points on a map, arising from multicomponent, reactive plumes originating
at a variety of source locations. The following sections describe how this
combination is accomplished within the framework of this overall modeling
scheme.
*Expressed as the fraction of the maximum value of Q.
11
-------
B. TRAJECTORY ANALYSIS
1. Grid Definitions
Horizontal trajectories are computed in time segments of equal length At
using a series of gridded wind fields. The various grids to be defined are
all derived from the NMC 47 by 47 octagonal grid. A portion of this grid
is in Figure 3a.
The NMC grid map is a polar stereographic projection of an octagonal area
centered about the north pole and extended sufficiently that all the northern
hemisphere north of the 18th parallel is included.2 Some of the equations
pertaining to the map are as follows. Let R be the radius in NMC grid units
from the north pole to a point on the map at latitude $ (radians) . The
equations relating R and are:
••"•""T-raH • U3)
* - 1 - 2 tan'1 ^S . (24)
The polar stereo projection is such that a basic grid interval
A = 3.81 x 107 cm (25)
is defined as the interval at 60 °N between consecutive grid intersections.
At any latitude the grid interval is given by
D » A/ (26)
where
-------
Let F be an arbitrary floating point number selected as the ratio of spacing
on thl advection grid to that on the sampling grid. Using 0 defined by (26)
we can write the expression for the grid intervals as functions of latitude
on the advection grid and on the sampling grid. We have
Da - D/I , (28)
at a,
D= = D /F , (29)
s as
for grid intervals on the advection and sampling grids, respectively. The
latitude dependence comes from the map factor in D. Let the X and Y direction
within the advection grid be defined as the horizontal and vertical directions,
respectively, on the NMC grid. In this case, the X direction is approximately
west-east over the U.S., exactly so, at 80°W longitude. There are three right-
handed coordinate systems to be considered within the original grid—one for
each of the previously named grids. For the coordinate system of the NMC
grid, we will consider the origin to be at the pole.* The negative Y axis
extends downward across North America at Longitude 80°W. Over North America,
the X coordinates are negative west of 80°W, and positive east of 80°W.
Therefore, all grids selected over the United States will have negative Y
coordinates in this system, and most would have either zero or negative x
coordinates for the left margin. This coordinate system is used only in one
read statement to identify the advection grid on the NMC grid. No other
input, and no gridded output is in this coordinate system. The advection
grid assigns the origin to the lower left corner of the advection grid.
The numbering of the X and Y axis is in increments of D .
a
The sampling grid is a subset of the advection grid. The origin may be at any
point on the advection grid, other than those on the upper or right boundaries.
The numbering of the X and Y axes is in increments of D .
An enlargement of the three grids is shown in Figure 3b. The heavy lines
are lines connecting the NMC grid points. They mark off the unit distance D
on the NMC grid and on the original grid. The lighter solid lines together
with the heavy lines mark off the unit distance on the advection grid for
which I =2, i.e., 2 units D on the advection grid equal one unit D on the
original grid. The dashed lines together with the solid lines mark off the
unit distance on the sampling grid for which F =4.
S
The coordinate designation (-3,-16) is the coordinates on the NMC grid of
the origin of the advection grid. Similarly the coordinate designation
(1,-12) is the coordinates on the NMC grid of the upper right corner of the
advection grid. The numbers 0 through 8 along the lower and left edges are
the scale in the advection grid system.
The sampling grid has its lower left and upper right corners at (1,1) and
(4,4), respectively, on the advection grid.
*Note difference between this and coordinate definitions for the same grid
in Random-to-Grid, see Figs. 3a, 9.
13
-------
-15
-20
-5
Figure 3a. Portion of truncated 47 x 47 NMC grid used in establishing
coordinate system for STRAM. The origin in this system
has been shifted to the North Pole
8
7
6
5
4
3
2
1
1-3, -16IQ
C
(1
*
-H-l-
-GJ
~U-[
1
-M-l-
4.+-1-
-t-t-t-
- t— 1— . -
~ ^™"l™n '
4-44.
-H-j-
T T'
-ttj:
°s
3 4
0
03
'.
5678
Figure 3b. Expanded view of the box in Figure 3a
14
-------
The stars are source locations. The dots in the sampling grid are particu-
larly specified sampling locations. The numbers in parentheses are coordinates
on the NMC grid.
2. Theory of Basic Equations
Let the horizontal projection of a portion of a trajectory over one square
of the advection grid be as in Figure 4. Assume that the curved line between
the two dots is the horizontal projection of the path of advection during
the interval At. The change of coordinates due to advection over the directed
distance As is expressed by resolving it into the components in the X and Y
directions. Specifically over the interval At, the components of the advec-
tion are:
t+At
Ax - r-
U(t') cos a(t')dt'
(30)
AY - =-
t+At
U(t') sin a(t')dt'
(31)
where U(t') is the horizontal wind speed. D is used to convert the advection
components from the length units of U to length units on the advection grid.
t
Y
X
Figure 4. Theoretical trajectories step
15
-------
Let the components of the wind in the X and Y direction along the trajectory
be designated by
U[t;X(t),Y(t)]
V[t;X(t),Y(t)], respectively.
The capital letters are used to distinguish these wind directions from those
used in Subroutine STRAC, which are parallel and lateral to the mean path of
advection rather than being aligned with the grid. Then the increments of
advection over a time At are given by
t+At
AX - — I U[t';X(t'),Y(t')]dt' , (32)
D I
a J
t+At
AY » — I V[t';X(t'),Y(t')]dt' . (33)
After the new coordinates X(t+At) and Y(t+At) are calculated for each incre-
ment of a given plume, a check is made to see if any of these values are off
the grid. If so, not only the affected plume increment, but the entire portion
of the plume emitted before the emission time of the affected increment is
deleted from further consideration.
3. Approximation of the Wind Integral by Finite Summation
The integration of the wind components over the trajectory is approximated
by a two-step iteration method involving bilinear interpolation in space and
in time (see Figure 5) .
The various positions in Figure 5 are best defined by equations:
^ = X(t) + U[t;X(t) ,Y(t)]At , (34)
Y = Y(t) + V[t;X(t) ,Y(t)]At , (35)
X2 = Xl
Y2 = Yl + V[t4-At;X1,Y1]At , (37)
X(t+At) = 0.5[X(t) + X] , (38)
Y(t+At) = 0.5[Y(t) + Y ] . (39)
16
-------
(X(t) ,
[X(t+At) ,Y(t+At))
2'Y2
Figure 5. Two-step approximation of trajectory step
Thus, the new position (X(t+At) ,Y(t+At) ) is found by adding the vector average
of two advection increments to the former position (X(t),Y(t)). The first
advection increment is calculated by advecting the plume increment for the
entire interval At using the wind effective at X(t),Y(t) at time t. The
addition of this increment to the position (X(t),Y(t)) yields the position
(X_,Y.). However, (X ,Y ) is not assumed to be the actual end of the advec-
tion increment because tne wind may change along the trajectory. To correct,
at least to some degree, for this situation, a second advection increment
is calculated using (X.. ,Y.. ) as the starting point and the wind at that point
effective at the end or tne time increment. The addition of this advection
increment to position (X-,Y.) yields position (X ,Y_) . Then the new plume
increment position (X(t+At) ,Y(t+At) ) is taken to Tse halfway along the line
from (X(t),Y(t)) to (X2,Y2).
The bilinear interpolation by which the effective wind components U(t) , and
V(T) are calculated works as follows. Let t and t^^ be the effective times
and
Let t and t
of the two wind maps closest to time t. Time interpolation weights t
t_ are defined by
17
-------
t - t
m
:2 " t , - t
m+1 m
:, * 1 - t_
*» - t< V1 '
(40)
(41)
Next, we note the location of (X(t),Y(t)) on Figure 5. The coordinate values
are between 3 and 4 for X(t) and between 2 and 3 for Y(t). Accordingly, we
define:
This yields
0(t)
X - X(t) - 3
q
x - i - x
p q
Y = Y(t) - 2
Y - 1 - Y
P
U(tm;3,2)
U(4'2) + '
qp
VpYq U(V3'3) + VpYq U(tm+l'3'3) '
U(t;4'3) + T U('4'3)*
m
qq l'
Similar equations hold for V(t), U(t+At), and V(t+At) .
Therefore the advection equations 32 and 33 are calculated from
1
Ax
AY
2D
2D
U(t) + U(t+At)
V(t) + V(t+At)
At
At
The subsequent positions are of course given by
X(t+At) * X(t) + AX
Y(t+At) = Y(T) + AY
(42)
(43)
(44)
(45)
(46)
(47)
(48)
(49)
(50)
18
-------
4. Examples of Application
An example of individual trajectories generated by application of a gridded
wind field can be found in the paper by L. L. Wendell (1972).3 The gridded
winds are derived by application of a random-to-grid subroutine (see III.B)
to the tower data.
C. CALCULATION OF INPUT PARAMETERS
1. Preparation of Data for Calculation of Trajectories
The calculation of trajectories does not depend on the selection of experi-
mentally determined values analogous to those enumerated in II.C.2 for plume
chemistry and dispersion calculations. However, it is necessary to specify
the location of the advection grid on the NMC grid, and the location of the
sampling grid on the advection grid, as well as the scales of the advection
and sampling grids. Also, it is necessary to specify the initial hour for
the run using the time zone that is used on the wind tape, the duration in
hours of the basic advection step, the time interval between printed con-
centration maps, and the total number of hours in the run. Meteorological
data is in Greenwich time (GMT).
Since all these relations are governed by input data that are read in from
punch cards prepared by the user, the detailed description necessary for
the user is found in III.C.4 on input parameters.
Furthermore, a matrix of map factors must be used to establish the scale on
the original grid (and consequently the scales on the other grids) as a
function of latitude. However, this is done automatically once the proper
punch card input data are prepared.
2. Calculations of Plume Chemistry and Dispersion Data
a. Chemical Reaction Rates
The rates of SO transformation in plumes have been discussed in detail in
a previous report.1* From this review, it is apparent that the wide variety
of conversion rates observed to date cannot be correlated in any simple or
straightforward manner. Moreover, our level of understanding at the present
time does not justify incorporation of anything more complex than linearized
kinetics parameterizations, at least for routine calculations.*
The previous review indicates that it is expedient to segregate oxidation
rates into rather broad condition classes, which are based upon meteorological
circumstances and ambient pollutant levels. Table 3 presents a summary of
*It should be noted here that research currently in progress is expected
to change this situation markedly within the next two years, by providing
more reliable parameterizations which are likely to be nonlinear in
character.
19
-------
TABLE 3. SUGGESTED REACTION RATE PARAMETERS
Condition Conversion rate Rate constant
class %/hr
Dry nighttime or daylight- 1 2.78 x 10
low hydrocarbon background
Nighttime or daylight 8 2.22 x 10~5
RH > 80%
Sunny, with high hydro- 5 1.39 x 10~5
carbon levels (imbedded
in urban plume)
Recommended overall value 2 5 . 56 x 10
for simple modeling purposes
selected conversion rates, based upon the ranges given in the review, for
each condition class. The rate recommended for trial use with the STRAM
code is 2%/hr. First-order rate constants appearing in Table 3 are defined
by the equation
i / moles
- rso2 a - rso ' - k c so2 (units are s
4
In concordance with the code, the units of k are reciprocal seconds. This
rate constant is entered into the code directly in subroutine form, as will
be described in Section III.
b. Washout Rates
Although wet removal of SO and SO can occur through a number of mechanisms,
it is appropriate for present purposes to lump these into parameterized
expressions based simply upon direct, nonreactive uptake of the species by
raindrops. For SO. aerosol one can proceed directly by assessing the
particle jize distribution and applying existing washout theory to determine
representative washout coefficients for a number of aerosol and rain spectra,
which car. be utilized directly for this purpose.5 One such set of curves,
pertaining to washout under general frontal rain conditions , is shown in
Figure 6. Using this curve in conjunction with a log-normally distributed
aerosol with a mass-median diameter of 0.05 y and logarithmic standard
deviation of 2 (roughly typical of observed sulfate distributions in
plumes) , one obtains a washout coefficient of
A _ = 0.0045 J , (52)
oCJ ~
4
where J is the rainfall rate in mm/sec , and A is in sec
20
-------
OS
a
o.
i
LL.
o
c*
< E
ee. o
&
o
en
B
-a
0)
4J
3
.-3
W
•H
•a
o
o
0)
N
•H
cn
0
14-1
-P
c
0)
•H
o
-H
u-i
M-l
0)
0
o
cn
•r(
fa
o
e5
8
cJ
I-
21
-------
Although SO washout does not follow the irreversible behavior typical of
aerosol removal, an exact treatment of this phenomenon lies beyond the
present capabilities of the STRAM code . Because of this , the recommended
procedure for estimating SO removal is to apply an effective washout coeffi-
cient based upon relative rates of removal observed in past studies of plume
washout behavior. These investigations have indicated that, as a rough
approximation at points more than a few kilometers downwind from the source ,
Aso2 - °-01 Aso4= • (53)
Values of the washout coefficients given above can be incorporated by chang-
ing Subroutine WSHOUT of the STRAM code, as described in Section III. To
utilize these effectively for long-range calculations, however, one must
incorporate a spatial and temporal description of rain intensity or climato-
logical averages of the same. Consideration of these spatial and temporal
variations is beyond the present capability of the STEAM code. In addition,
the current STRAM code assumes that washout rates are equal to zero.
c. Deposition Velocities
Deposition velocities for SO_ have been observed by a number of investigators
to be in the range of 1 cm/sec. The only known measurements of 30= deposition
velocities in the real atmosphere have fallen in the range of 0.1 cm/sec.6
Accordingly, the values
V, =1 cm/sec (54)
0.1 cm/sec (55)
are recommended for present test calculations. As with previous parameters
these are entered directly into STRAM in subroutine form as described in
Section III.
d. Plume Rise
For long-range calculations, it is usually most expedient to presume the
plume has attained its total effective rise and, rather than burden the user
with costly repetitive calculations, STRAM uses the virtual stack height
throughout to represent the effective plume height. This height may be
changed at regular intervals, typically hourly. The user may wish to rewrite
the Subroutine PLMRIZ to use Briggs' modified formulae or any of a number of
alternative techniques.
e. Stability and Mixing Height
Like virtual stack height, stability and mixing height data may be read in
at regular intervals, typically hourly.
22
-------
f. Calculation of Dispersion Parameters
The dispersion parameters a and a are calculated for distances out to
100 km using analytical formulae dlrived from the curves in the Workbook
of Atmospheric Dispersion Estimates by D. B. Turner.7 This was a quite
straightforward procedure for ay since the curves given by Turner on log-log
paper are approximately straight lines. For a the Turner curves are not
straight lines. However, the curvature is slight enough that good approxima-
tions for downstream distances of 1 to 100 km can be made by assuming the
curves to be straight lines through the values for those points. Under
these conditions extrapolation of the az formulae to distances outside this
range is not to be considered. For distances greater than 100 km formulae
by Heffter are used.8
Since stability changes during real-time dispersion, the formulae to be given
below cannot be used as such. If this were done, for example over a time
interval when the stability is increasing, the dispersion parameter values
would be represented as decreasing functions of tine. Therefore, instead
of computing from a formulae such as av(s) (where s is total travel distance)
Y
or Oy(t), the derivative is used as follows:
a (s+As)
da
a (s) + As ,
y ds
s+1/2 As
(56)
da
a (t+At) =« a (t) + At
y y dt
t+l/2 At
(57)
Similar equations apply- for a .
z *
The integral formulae for ay and a for travel distances less than 100 km
are of the following forms
Vs'00
Vs'00 = Zas
102m 5 s < 105m
103m < s £ 105m
(58)
(59)
Here a is a stability index designated as A, B, C, D, E, or F for the Pasquill
stability categories. The formulae are not valid if the stability has changed
over the travel distance s. The coefficients Ja and Z yield a and az in
meters when s is also measured in meters. They are given in Table 4.
The derivative forms used in carrying out the actual computations are, of
course,
da
ds
= 0.9 Y s
a
-0.1
102m 5 s 5 105m
(60)
23
-------
TABLE 4. COEFFICIENTS FOR DISPERSION PARAMETERS FORMULAE
a
A
B
C
D
E
F
Y
ot
0.36
0.25
0.19
0.13
0.096
0.063
ClO _ M
2 b •• 1
— - - b Z s a
as a a
Z
a
0.00023
0.058
0.11
0.57
0.85
0.77
3 5
10 m < s S 10 m
b
a
2.10
1.09
0.91
0.58
0.47
0.42
^
(61)
Note that at s = 0, STRAM sets a = 0 and a =0.
The integral formulae for ay and a given by Heffter for travel distances
greater than 100 km are
a (t) = 0.5t , t in sec , a in m (62)
a (t) = (2K t)1/2 , K = 5 m2 sec"1 (63)
z z z
The derivatives used in computation when s is greater than 10 m are
da
-rr^- = 0.5 , (64)
dt
- 1/2 (2K )1/2t-1/2
at z
(65)
24
-------
SECTION III
USER'S GUIDE
A. INTRODUCTION
1. General Data Processing Requirements
The data processing requirements for the STRAM system include the operation
of two computer programs: Random-to-Grid and STRAM. The input to the STRAM
program contains gridded data for wind components generated by the Random-to-
Grid program. Because the STRAM program incorporates the map factor for the
NMC grid (see Section II.B.I), the gridded wind data must be defined at the
points of a rectangular subset or some integral subdivision of this grid.
This grid is called the advection grid. An example of how the gridded wind
component data may be prepared from randomly spaced station data using the
Random-to-Grid program is discussed in Section III.B.
The program is designed to accept from cards at evenly spaced intervals
DTIME new values of stability, mixing height, and source strength for
each chemical component. The program is equipped with one option that
reduces the requirement for stability and mixing height data to one 24-
hour cycle and another option that reduces the source strength requirement
to one set of constant values. Details are given in Section III.C.4.
The STRAM output data tape is discussed in Section III.C.5.
2. STRAM System Flow Chart
The STRAM system flow chart is given in Figure 7.
B. RANDOM-TO-GRID
1. Description, Function, Constraints
The function of this program is to use upper air multi-level wind speed
and direction observations from a randomly distributed network of reporting
stations to compute gridded fields of the average u- and v-component winds
in a layer of a specified thickness and height above the ground.
The major inputs to the program are station wind direction and speed data
from rawinsonde observations. These data are on magnetic tapes written at
the Air Resources Laboratories in Silver Springs, Maryland. The output con-
sists of a series of u- and v-component wind fields at 12-hour intervals.
These data are written to tape in a form which may then be used by the
STRAM program.
25
-------
READ GOVERNING
PARAMETERS
FROM CARD
PRINT GOVERNING
PARAMETER
VALUES
PARAMETER
VALUES CREDIBLE
PIBALAND
RAW1N SONDE
DATA
RANDOM-TO-GRID
PROGRAM
BASIC ADVECTION AND SAMPLING LOOP
L READ WIND STABILITY, MIXING HEIGHT,
SOURCE STRENGTHS
2. CALCULATE VALUES FOR ALL MODaED
PHYSICAL PROCESSES
3. ENTER NEW SAMPLING CALCULATIONS ON
TAPE
4. MAKE SCHEDULED OUTPUTS TO LINE PRINTER
FALSE - DO NEXT
ITERATION
TRUE
PRINT FINAL CONCENTRATION
MATRICES
END FILE OUTPUT TAPE
STOP
\ /
CONCENTRATION
MATRICES
TAPE
Figure 7. System Flow Chart
26
-------
The program consists of a main program and four subroutines, plus six
w subroutines which are used in transmittal of information from the data
tapes.
The computation of the gridded wind fields takes place in the main program.
In calculating the values at each point, the values at stations nearest
the grid point are averaged, weighing the station values according to the
reciprocal of the square of the distance between each station and the grid
point. To avoid having to search for the nearest stations and calculate
their distances from each grid point for every observation time, this opera-
tion is performed only once, at the beginning of the program. This requires
that the latitudes, longitudes, and identification numbers of all stations
which may be used later in the program must be input on cards (see input
cards No. 6, 7, 8, and 9 in Section III.E.4, "Input Description and Formats")
before entering the wind field calculation loop (see the second and third
boxes in the "Flow Diagram - Major Functions," Section III.B.2).
The following subroutines perform specific functions in the program.
w
FASCND orders the elements of a vector according to their magnitude from
the least to the greatest. FASCND is called in the section of the main
program which orders stations according to their distances from each grid
point.
v MAPTIME is used to advance the given time by a specified interval in hours,
in order to select data from the desired observation times, and skip informa-
tion from other times.
RESON2 reads rawinsonde and pibal station data from a specified file on the
input tape.* It uses LYRWND to compute the average wind for the specified
^ layer at each station. It returns to the main program u- and v-component
layer averages for all available stations at a given observation time.
Six subroutines are employed to read the 32-bit word data on the input
tapes, and convert to 36-bit Univac words. These subroutines, READ1, READ2,
READ3, BUMP, NEXT, and SPLIT, are discussed in "Input Description and Formats,"
w Section III.B.4.
2. glow Diagram - Major Functions
The execution sequence of the major operations performed in PJ4GRD is sche-
matically represented in Figure 8.
*RESON2 searches for the desired observation time on the input tape. The
tape must be prepositioned however, using appropriate job control language
(see Section III.B.3) before the execution of the program.
27
-------
Input of control parameters for
windfield generation
Read in station identification
and location information
-*-No-
At each grid point generate two vec-
tors, one containing the distances in
ascending order to the nearest ten
stations, and the other containing
the subscripts of these stations
Specify starting and ending times in
the series, and layer boundary heights
,
1
\
Call RESON2 to, obtain u- and v-com-
ponent layer averages at the stations
J
1
\
Compute gridded u- and v-wind com-
ponents from layer wind averages
at the randomly spaced stations
* f
Read upper wind data for
a single observation time
Compute the layer averaged
wind at each station
I J
Write the gridded u- and v-windfields
to tape or disk for further use
\
,
_/ Is this the final map N.
"N. in the series? /
Y$
s
res—<
Is another series of
maps to follow?
Figure 8. Random-to-Grid Flow Diagram - Major Functions
28
-------
3. Executive Control Language and Deck Setup
The following set of control cards is required to execute the Random- to-
Grid program on the UNIVAC 1110.
RNGRD
EXECUTIVE CONTROL LANGUAGE (ECL)
@RUN , PRI/OPTS RUNID/CORE/TAPE/DISC , ACCT , PROJ/CLASS , RUNTIME , PAGES/CARDS , STIME
@ASG,T O.UAL*WINDTAPE . , 16N ,999999 (ASSIGN HEFTER WIND TAPE)
@USE 10 , OUAL*WINDTAPE . (ACCESS HEFTER FILE ON FORTRAN UNIT 10)
@MOVE QUAL*WINDTAPE. ,N (POSITION HEFTER TAPE TO FIRST DESIRED
FILE, EACH FILE HAS 1/2 MONTHS DATA,
MOVE HEFTER FILE N NUMBER OF FILES)
@ASG,T/W QUAL*RNGRDTAPE . , 16N, 888888 (ASSIGN RNGRD OUTPUT TAPE)
@USE 3 , QUAL*RNGRDTAPE . (WRITE RNGRD TAPE ON FORTRAN UNIT 3)
@ASG,A QUAL*PROGFILE. (ASSIGN RNGRD PROGRAM FILE)
@XQT QUAL*PROGFILE . RNGRDPGM (EXECUTE RNGRD ABSOLUTE PROGRAM)
(INSERT THE DESIRED RNGRD CONTROL CARDS HERE)
@FIN
@@
NOTE: ALL UNDERLINED PARAMETERS TO BE SUPPLIED BY THE USER,
RUN CARD PARAMETERS, QUALIFIERS, FILE NAMES, REEL
NUMBERS AND RNGRD CONTROL CARD DATA.
29
-------
4. Input Description and Formats
Card input to the KNGRD program is detailed in Table 5.
TABLE 5. CARD INPUT DESCRIPTION
Card No. Columns
1 1-64
2 1-5
2 6-10
2 11-15
2 16-20
2 21-25
2 26-30
2 31-35
2 36-40
Format Name
16A4 ITLE
15 NP
15 NSX
15 NSY
15 NFX
15 NFY
15 KTAPE
15 NSTL
15 INT
Definition
Title of calculation run
Number of stations for which data
may be available
Starting x-direction grid point
value. x=0 for NSX=1
Starting y-direction grid point
value. y=0 for NSY=1
Final x-direction grid point value
Final y-direction grid point value
Write wind fields to tape : l=yes ,
0=no
Minimum number of stations to be
used in any grid point calculation
Print interval, print out one of
every "INT" u- and v- fields. If
INT=0, no fields will be printed
1-80 1615 KSX(I),
1=1,NFY
1-80 1615 MFX(I),
1=1,NFY
1-10 F10.0 RCH
11-20 F10.0 XNMC
21-30
F10.0 YNMC
For variable left boundary. Speci-
fies the lowest x-value in each row,
from 1 to NFY
For variable right boundary, Speci-
fies the lowest x-value in each row,
from 1 to NFY
Radius in grid units within which
all stations are to be used in an
interpolation for a given point
With respect to the octagonal 47 x 47
NMC grid (see Figure 9) the x-value
of grid point (1,1) on the inter-
polation grid
With respect to the octagonal 47 x 47
NMC grid (see Figure 9) the y-value
of grid point (1,1) on the inter-
polation grid
30
-------
TABLE 5 (continued). CARD INPUT DESCRIPTIONS
Card No. Columns Format Name
6*
Definition
31-40
41-50
1-80
F10.0
RFC
F10.0
1615
THET
KSTA(I),
Resolution factor. The ratio of
the interpolation grid spacing to
the NMC grid spacing. RFC=1.0 yields
grid with the resolution of the NMC
grid (Ax V380 km at 60°N)
Number of degree of clockwise rota-
tion of the x-axis on the interpola-
tion grid from the x-axis on the
NMC grid
Last four digits of each station
number for stations whose observa-
tions might be used
Latitudes of stations whose observa-
tions might be used
Longitude of stations whose observa-
tions might be used
Names or abbreviations of stations
whose observations might be used .
Three letter abbreviations for the
months of the year
Last two digits in starting year
of the series
Starting month of the series
Starting day of the series
Starting hour of the observation
series
Year of the last observation time
in the series
Month of the last observation time
in the series
Day of the last observation time
in the series
Hour of the last observation time
in the series
Perform random-to-grid for another
series of maps after this one:
l=yes, 0=no
*There might be more than one card necessary to read in this array.
7*
8*
9*
10
11
11
11
11
11
11
11
11
11
1-80
1-80
1-80
1-48
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
16F5 . 2
16F5.2
16 (IX,
A4)
12A4
15
15
15
15
15
15
15
15
15
1=1, NP
1=1, NP
NAMSI(I)
KMON(I) ,
1=1,12
IYS
MS
IDS
IHS
IYL
ML
IDL
IHL
KNTHR
31
-------
TABLE 5 (continued). CARD INPUT DESCRIPTIONS
Card No. Columns Format
11
12
12
13
13
13
13
13
Name
Definition
46-50
1-5
6-10
1-10
11-20
21-30
31-40
41-45
15
P5.0
F5.0
F10.1
F10.1
F10.1
F10.1
15
KREW
ZB
ZT
FLATMN
FLATMX
FLONMN
FLONMX
IDT
Rewind input tape before starting
this series : 1 • yes, 0 = no
Height above ground of the bottom
of the desired layer
Height above ground of the top of
the desired layer
Minimum latitude of station search
area
Maximum latitude of station search
area
Minimum longitude of station search
area
Maximum longitude of station search
area
Selected map interval for analysis.
IDT is an integer multiple of the
smallest time interval (in hours)
between observed data, and less than
or equal to 24 hours
Magnetic tape input to the RNGRD program consists of pibal and rawinsonde
data which has been reduced from ETAC Upper Air data tapes by the NOAA Air
Resources Laboratories and converted for use on EPA's Univac 1110 computer.
The data are in 80-column card-image format.
Each month's data are divided into two files on the tape, the first covering
days 1-15, and the second covering days 16 through the end of the month.
Data for each observation time are headed by a descriptor record (record
type #1 in Table 6) containing the data and time, number of reporting sta-
tions, and number of records to be skipped to get to the next observation
time. Data are reported four times daily (0000, 0600, 1200, and 1800 GMT).
Because of the relatively sparse data coverage at 0600 and 1800, RNGRD uses
only the 0000 and 1200 GMT observations.
Within each observation time, data are grouped according to station number,
and each station's data are organized in ascending order by level of observa-
tion. Data for each station are headed by a descriptor record (record type
#2 in Table 6), giving the station number, latitude and longitude (in
degrees and hundredths of degrees), station elevation (in meters), elevation
of the surrounding terrain (in meters, averaged over a 1° by 1° square),
and the number of observation levels reported for the station.
32
-------
10
15
NMC OCTAGONAL GRID
20 25 30
35
40
45 47
-
rw/?. a...
; t I *^ rf ' • .
.1
45 47
Truncated 47x47 grid
The pole point is (I,J) = (24,24)
Figure 9. Truncated NMC grid, used as a basis for interpolation grid
positioning specified on input card #5 to RNGRD
33
-------
The data for each observation level resides on a single logical record
(record type #3 in Table 6), and includes the height of observation (in
meters), wind direction (in degrees), and wind speed (in meters per second).
TABLE 6. WIND DATA TAPE INPUT DESCRIPTION
Record
Type
#1
#1
#1
#1
#1
#1
#2
#2
#2
#2
#2
f2
#3
#3
#3
Character No.
(column)
1-3
4-7
8-9
10-11
12-15
16-20
1-5
6-10
11-17
18-22
23-27
28-29
1-5
6-8
9-12
Format
A3
14
12
12
14
15
15
F5.2
F7.2
F5.0
F5.0
12
F5.0
F3.0
F4.1
Name
IM
J2
J4
J5
NS
NREC
KSTAT
FLAT
FLON
KHT
KTZ
NLP
KZ
DIR
SPD
Definition
Month of observation
Year of observation
Day of observation
Hour of observation (GMT)
Number of stations reported for
this observation time
Number of records to skip to reach
next observation time
Station number
Latitude (degrees and decimal fraction)
Longitude (degrees and decimal fraction)
Station height (meters)
Terrain height (meters)
Number of observation levels
Observation height (meters)
Wind direction (degrees)
Wind speed (meters per second)
The data on these tapes were written in the 32-bit word configuration of
the IBM equipment used by ARL. As such, they must be converted to the
36-bit word Univac configuration. This is accomplished in RNGRD by use
of the six input and conversion subroutines, READ1, READ2, KEAD3, NEXT,
SPLIT, and BUMP.
Call sta-cements in RESON2 to the first three of these subroutines correspond
to and perform the tasks of Fortran READ statements which would read and
return data from record types #1, #2, and #3 described above.
34
-------
Subroutine NEXT is called by the three READ subroutines to obtain informa
tion from the tape. Subroutines SPLIT and BUMP are used by the READ sub-
routines to rearrange the IBM bit structure into the Univac configuration.
5. Output Description and Formats
There are two forms of output from this program. The most important of
these is the tape output. This contains gridded u- and v-component wind
fields appropriate for use with the STRAM program. The second form is
the printed output which gives some of the input information, and, if
desired, station layer-averaged u- and v-component winds followed by
gridded u- and v-component wind field values at each observation time.
Tape Output. Information is written to a tape or disk using the three
unformatted binary write statements below.
WRITE(3) IYR, MON, IDA, IHR, ZB, ZT, NSX, NFX, NSY, NFY
WRITEC3) ((UG(I,J), I=NSX,NFX), J=NSY,NFY)
WRITE(3) {(VG(I,J), I-NSX,NFX), J=NSY,NFY)
The first statement writes the year (last two digits) , month (integer) ,
day, and hour of the observation time, the height in meters above ground
of the bottom, and of the top of the averaging layer, the starting and
final x-coordinates , and the starting and final y-coordinates on the wind
field grid. The second statement writes the entire gridded u-component
array, and the third writes the entire gridded v-component array.
This sequence is repeated at each of the observation times for which wind
field data is computed. There are no end-of-file markers between records
or observation times.
Printed Output. The printed output may be divided into two types, standard
and optional. The standard output is always present. The presence of the
optional output depends on the value of the input control parameter INT (see
"Input Description and Formats," card number 2).
Standard output includes the following:
• Page 1: listing of station numbers, latitudes, and longitudes which
were read from input cards.
• Page 2: title of calculation run; station names; x-coordinates of
stations; y-coordinates of stations; starting and final x-grid coordinates
for each u; the height above ground of the bottom and the top of the
averaging layer; and the minimum and maximum latitudes and longitudes
of the station search area.
35
-------
The optional output applies to information from each observation time.
This information is printed out once in every "INT" observation time.
If INT-0, no optional information is printed. For each observation time
for which data are to be printed, the following information is included
in the output:
• Page 1 (Optional Printout): heading containing time of observation
and boundaries of the layer; a table giving station number, latitude,
longitude, station elevation, and the layer-averaged u- and v-component
winds, for each station for which data were read from the tape.
• Page 2 (Optional Printout): u- and v-component station layer winds,
printed to check their transfer from RESON2 into the main program; the
number of stations used in the interpolation for each grid point,* gridded
layer u-components (for the observation time, and the gridded layer v-
components* for the observation time.
6. Diagnostic Messages
Subroutine RESON2. If there is no observational data available on the
source tape for the designated time, the following message will be printed,
with the number of the month designated by mm, the day by dd, the year by
yy, and the hour by hh:
"RECORD NOT AVAILABLE FOR mm/dd/yy, hhZ"
Subroutine LYRWND. There are two diagnostic messages which may appear in
the listing of layer-averaged u- and v-component winds (see Section III.B.5,
"Page 1 (Optional Printout)"). Since no winds can be interpolated for a
layer which lies entirely below, or entirely above all of the wind observa-
tions, no layer averages are calculated in such cases. The diagnostics
associated with this situation are:
" - ALL OBSERVATION BELOW LAYER," and "ALL OBSERVATION ABOVE LAYER."
C. STRAM PROGRAM
1. Description, Functions, and Constraints
The purpose of -his program is to model the distribution of ground-level
concentrations over a rectangular area approximately 1000 km on each side
resulting from emission of multicomponent plumes from several elevated
sources within the area. Broadly, the program can be divided into three
functions—a. trajectory function, a dispersion function, and a sampling
function.
*The values appear in the printout as they would appear on the grid, i.e.,
with the value in the lower left corner corresponding to the lower left
grid point.
36
-------
The executive control and trajectory function are in the main program.
The main program reads in all input parameters, prints these out, defines
the advection and sampling grids on the 47 by 47 NMC grid, calls the major
subroutines, computes plume trajectories, and governs all output, of which
there are two kinds, line printer and magnetic tape.
The input and output parameters are described in detail in III.C.4 and
III.C.5, respectively. The basic equations used are given in Section II.
Therefore, this section will be primarily of a descriptive nature, explaining
how these things are used, how the program functions, and what constraints
the user must be aware of.
a. Grid Definition and Spacing
There are three grids to be considered by the user, the 47 by 47 NMC grid,
the advection grid, and the sampling grid. Average concentrations are
computed for all intersections of the sampling grid (as well as for up
to ten particularly specified other locations within the sampling grid),
which is a rectangular area within (or equal to) the advection grid.
Each of the three grids has its own coordinate system. The definitions
of these and the relationships among them are described in II.B.I.
The reasons for providing the capability to make the sampling grid smaller
than the advection grid involve both dispersion physics and computer logis-
tics. The number of intersections that can feasibly be specified for the
sampling grid are limited primarily by computer core size in view of all
the other core requirements in the code. The spacing between consecutive
intersections is limited by the fact that if the spacing is too large,
much of the concentrations from the plumes will go undetected. Therefore,
it may be necessary to settle for a sampling grid that covers only part
of the advection grid. Also, in this way, contributions from sources
outside the sampling grid may be considered. The present dimensions of
the code allow for a sampling grid of not more than 13 by 13 equally spaced
intersections.
The location and subdivision of the sampling grid are defined by the input
parameters specifying the locations on the advection grid of the lower
left and upper right corners of the sampling grid, along with the ratio
of the scale of the sampling grid to that of the advection grid (see
III.C.4). The specifications can be any positive floating point numbers
subject to the constraints of dimensions of the sampling and advection
grids.
The advection grid is limited to an array of not more than 17 intersections
in the horizontal direction of the NMC grid and not more than 13 inter-
sections in the vertical direction. The location and subdivision of the
advection grid are defined by integer input parameters specifying the loca-
tions on the NMC grid of the lower left and upper right corners of the
advection grid, along with the scale ratio of the advection grid to the
37
-------
NMC grid (see III.C.4). Since these are all integer specifications, the
boundaries of the advection grid (also called the rectangle) are along
lines connecting the intersections of the NMC grid, and all intersections
of the NMC grid within and on the boundaries of the advection grid are
also intersections on the advection grid. However, the total grid within
the rectangle may be a finer grid than the NMC grid, formed by integral
subdivision. That is, the spacing on the advection grid may be two or
three times as dense as that on the NMC grid.
Usually the spacing of the intersections on the advection grid is either
identical to that of the NMC grid or twice as dense. In the first case,
IDIV is 1; in the second, IDIV is 2. The choice depends on the spacing
of gridded wind observations provided. If gridded winds are created
especially for the NMC grid, IDIV must be 1. The NMC wind tapes available
through the National Center for Atmospheric Research (NCAR) are examples.
But this spacing is approximately 350 km over the United States, and is
hardly optimum on a grid that is approximately 1000 km on a side. If the
user can make his own gridded wind tapes using synoptic weather station
data over some portion of the United States for input data to a random-
to-grid type program, he may conclude from inspection of a map of these
stations that a denser grid would reflect the detail of the wind field
more fully. Using this type of data, our choice for spacing in the
advection grid is to double the grid density, i.e., to make IDIV equal
to 2. This gives an average distance of about 175 km between successive
grid intersections over the United States.
An estimate of the grid spacing on the sample grid can be made by con-
sidering the effective width of a Gaussian plume, which we will consider
to be 6cr . For a plume increment that has traveled 200 km, this figure
would be about 50 km. Therefore, grid spacing on the sampling grid of
this order would be appropriate. For example, if DIV is 4.0, the grid
spacing on the sampling grid is one-fourth that on the advection grid,
and a typical grid spacing on the sampling grid would be about 45 km
(if IDXV is 2) .
b. Other Initial Functions
After the grids have been defined and all input parameters have been read
in, there is a printout of the essential information (see III.C.5) . Next,
certain of the input parameter values are checked for acceptability. If
an unacceptable value is found, it is printed out and the run is stopped
(see III.C.6).
If the run survives the diagnostic check, the next computations are of
map factors for the intersection points of the advection grid (see II.B.I).
These factors are necessary because the scale on the NMC grid is a function
of latitude, and the basic grid spacing holds only at 60eN. The next step
consists of computing the latitudes and longitudes of all sampling grid
intersections and printing these out.
38
-------
The last step before beginning the major loop, the advection loop, is to
initialize appropriate values. The average concentration matrices CHISAV
and CHILAV are zeroed out since their values are computed by summation
with the same matrix element on both .sides of the equal sign in the Fortran
statement. The matrices containing a and a for each plume element,
SIGMY and SIGMZ are zeroed out for the same reason. The variable TIME,
which is incremented once in each iteration of the advection loop by
TIME = TIME + DTIME ,
is also set to zero. This variable gives the current running time in hours
at any time during the run. The variable TMAP, which is incremented also,
once, during each iteration of the advection loop
TMAP - TMAP + DTIME
is also set at zero. TMAP is used to call in new gridded wind maps at the
appropriate times. Since the interval between maps is DMAP hours, (see
III.C.4), TMAP is allowed to increase until it equals DMAP, at which time
a new map for each of U and V, the wind components in the horizontal and
vertical directions, respectively (on the NMC map), is read in. At this
time, TMAP is reset at zero.
Also, the vector NINCV containing the number of increments in each of the
NSOURC plumes, is initialized at 1. This is because there is just one
plume increment to be advected from each source as the run begins. Each
element of the vector is incremented by 1 at the end of each major
iteration, thus keeping track of the total number of increments in each
plume.
c. Beginning the Major Loop
With all these values initialized, the major iterations are ready to begin.
The first time through, two sets of wind maps are read in. The first set
consists of maps for U and V at the initial time of the run, and the second
set is for U and V DMAP hours later than the initial time. These two sets
of wind maps are used for the bilinear interpolation in time and space
(see II.B.3) that is used to calculate the wind for any point on the grid
at any time between two map times. The fraction TMAP/DMAP is used to
calculate the time interpolation. On subsequent iterations, if TMAP is
less than DMAP, the wind map reading part is skipped. When TMAP is not
less than DMAP, the last of the two sets of wind maps in core is trans-
ferred to the core locations of the first set of maps and a new second
set is read in.
Following the wind map section, a Subroutine KALEND is called that com-
putes the hour (0 to 23), the day, the month, and the year, effective at
the end of the current major iteration. Next, the time interpolation
values, with respect to wind maps CA and CB, are calculated. Then the
39
-------
matrices CHI and CH, the average concentrations over 'the interval DTIME
for both the gridded and particular sampling locations are set at back-
ground levels.
The stability, mixing heights and source strengths are read in as functions
of time according to two options available to the user, see III.C.4.
d. The Advection Routine
The advection process is accomplished by three nested Do loops. The outer
one is for sources, the next is for increments from each source, and the
inner one consists of two steps to carry out the two-part iteration for
each advection step described in II.B.3.
The X and Y positions on the advection grid of the plume increments are
kept track of in the matrices XP(K,NP) and YP(K,NP), respectively, where
the first index is for source and the second for increment. The value of
NP is in order of release from the source, i.e., NP » 1 for the first
increment released. This increment will usually be the furthest from
the source. As the run progresses from one major advection step to another,
the individual plume increments retain their index values while the position
values are altered by addition or subtraction, depending on the direction
of the advection.
Within the two inner Do-loops, the advection calculations are made for
each increment of each plume. The calculated values include the new X
and Y positions of the increments on the advection grid, the average wind
speed in cm sec (UU(K,NP)) for each increment, and the total distance
advanced in cm (DS(K,NP)) for each increment. The wind speed and travel
distance matrices are entered later into Subroutine STRAC for use in
diffusion computations.*
Just before the end of the loop for plume increments is completed a Sub-
routine PLMRIZ is called that computes the effective stack height for
the emission if OPTION » 1.
Immediately after the Do-loop for the increments has been completed for
all increments, but still within the Do-loop for the sources, another
Do-loop for the increments is called to ascertain whether or not any
plume increments have wandered off the advection grid. If it is found
that the ith increment is off the grid, this increment and all those emitted
before the ith increment are removed from further consideration. Since this
removes increments 1 through i, the data pertaining to the remaining incre-
ments, positions, wind speeds, sigma values, total distance and time
traveled, are all reindexed with the data that was for increment i + n,
now being for increment n.
*See addenda
40
-------
e. Subroutine STRAC and INTGRT
Following the reindexing of the increment data, Subroutine STRAC is called
to determine ground-level concentrations at positions relative to the plume
centerline and source location. Subroutine STRAC fulfills this function
by providing numerical approximations to the solutions of equations (5a),
and utilizing the resulting Q. values in conjunction with (2) and (3) to
determine corresponding concentration values.
The subroutine call necessary to execute subroutine STRAC is given as
follows:
CALL STRAC (M,SUBMIN,SUBMAX,FACTOR,MODE,
TIME,OPTION,NCOMP).
The arguments of the call list are defined in the glossary of computer
terms. Additional information is exchanged between the subroutine and
the main program via the COMMON statement which reads as follows in the
main program:
COMMON NINCV(IO),QQ(6,100,10),CBKG(6),DS(10,100),UU(10,100),H(10),
SIGMY(10,100),CGRND(7,100,6),ISTAB(10,100),HH(10,100),XP(10,100),Y
P(10,100),XM(10,100),SIGMZ(10,100),PTIME(10,100),XRES(10),YRES(10)
,VIRTH(10,100)
In the calling sequence given above M is the number of the particular
source, and MODE dictates the calculation sequence in STRAC as indicated
in the glossary. It should be noted that for computational efficiency
MODE should be set as low as possible. Unnecessarily high values of
MODE (e.g., MODE set equal to 4 and a nonreactive plume) will result in
essentially the same output, but will direct STRAC to perform nonessential
calculations. Setting OPTION to 1 or 2 dictates whether calculations in
STRAC utilize the conventional Guassian or limited mixing-height distri-
butions. Under option 2 STRAC assumes that emissions are within the
mixed layer at all times regardless of stack height.
SUBMIN and SUBMAX are minimum and maximum allowable downwind increment
lengths for numerical integration in the x direction. Choice of these
values is dictated by accuracy requirements, and also by the rates of
plume transformation and removal. A reasonable rule of thumb is to try
to maintain the increment lengths such that no more than 5% of the plume
is transformed or removed in any given downwind step. With MODE = 1
(inert, nondepositing plume) no numerical integrations are performed,
and thus downwind increment lengths will be as large as allowed by plume
drift during the current time step. Under this condition, no values for
SUBMIN and SUBMAX are required.
Since proportionately more transformation often occurs early in a plume's
history, it is desirable to relax the downwind grid spacing requirement
with increasing distance. This is accomplished by FACTOR, which allows
the x grid spacing to be expanded for each successive increment as follows:
41
-------
Ax
Ax
step n+1
* FACTOR (66)
step n
STRAC limits expansion by Equation (66) to Ax intervals less than or equal
to SOBMAX.
The remainder of the input variables are self-explanatory from Table 7.
As noted above, much of the essential information is exchanged between
STRAC and the main program via the COMMON statement. In particular, the
computed ground-level concentrations
CGRND(I,N,K)
I » crosswind (y) position index
N = downwind increment index
K = component index,
as well as a crosswind spread parameter
SIGMY (M,N)
M « source index
are conveyed to the main program in this manner.
The layout of Subroutine STRAC is shown in Figure 10. Upon being called
from the main program, this subroutine first initializes pertinent varia-
bles. Downwind distance and plume depletion factors, as well as the grid
control variable TOTAL are set equal to zero. Following this the initial
grid spacing in the x direction is set. If MODE equals 1, no numerical
integrations are performed, and thus no grid subdivisions are required.
The subroutine simply computes the concentrations at the given increment
DS(M,N) and proceeds. If MODE does not equal 1, then the initial grid
spacing is set to the value of SUBMIN or DS(M,N), whichever is smaller.
x is then advanced one-half step, plume dispersion parameters are calcu-
lated, and values of the virtual source terms appropriate to that distance
are obtained by adding the computed values for generation of material in
the half-interval by reaction and deposition (taken to be zero for the
first half-step). This completes the first step of the two-step Runge-
Kutta integration.
Following this initial sequence, an integration over the full interval
DX is performed. This is accomplished by calculating generation rates
over the full interval using concentration estimates computed at the half
interval on the basis of source terms calculated from step 1. Virtual
source terms for that plume increment are then adjusted to reflect the
gain of material during the distance step DX, arriving at a new set of
terms QQF(K) = QQ(K,N,M) for each component. Various steps of this calcu-
lation are bypassed, according to the value of MODE.*
*See Addenda
42
-------
i
zf
§
•H
4J
3
s
w
fi
tn
«
•H
-o
43
-------
Once the (QQ(K,N,M) K=1,NCOMP) array has been calculated at the full inter-
val DX, ground-level concentrations may be computed. This, however, is
performed only when the end of the DX interval happens to fall on the end
of the major plume interval DS(M,N). If this is not the case, an additional
downwind step is taken and the integration procedure repeats until this
condition is satisfied. Under conditions when a major plume interval has
been traversed, ground-level concentrations are calculated at seven points
across the y-axis of the plume, starting at y=0 and proceeding to one side
at intervals of l/2a . Equation (2) or (3) is used for this purpose,
depending upon whether a bivariate-nonnal (OPTION=1) or limited mixing-
height (OPTION*2) model is being employed.
Regardless of whether or not a major plume interval is completed on a given
step, a grid-space adjustment is always attempted before the next integration
step is perforiued. If a major plume increment has been completed on the last
interval, the next step size, DOS, is set equal to that of the preceding step
(DDSOLD) or to that of the next major interval DS(M,N), whichever is smaller.
If a major interval has not been completed on the current step, then DOS is
set equal either to its current value times the grid expansion factor FACTOR,
or to SUBMAX, whichever is smaller. Furthermore, if DOS exceeds the distance
to the end of the major interval DS(M,N), then it is set equal to this dis-
tance. This provision guarantees that calculations ultimately will coincide
with points at the ends of the major intervals.
Upon reaching the end point of a major interval and printing, the increment
index N is reduced by 1, and computations for the next increment are
initiated. N equal to zero signifies that calculations for the total
plume have been completed; under this condition control is returned to
the main program.
After an x grid-space adjustment has been made, the subroutine re-enters
the first step of the two-step Runge-Kutta integration. In contrast to
the initial entry which assumed all generation processes were negligible,
each subsequent entry calculates values of the generation of material in
the half-step. This is performed in the same manner as used previously
for the full-step integration, calling Subroutine INTGRT to approximate
equation (10) for the chemical reaction term, and using the formulae in
Table 1 to calculate corresponding terms for dry-deposition and washout.
Following calculation of the generation terms the subroutine advances
the value of x one-half increment, adjusts the source terms as appropriate
to the .iew half interval position and proceeds.
As indicated above, the STRAC subroutine calculates through NINCV(J) major
increments of plume, storing ground-level concentrations CGND(I,N,K)
(1*1,7; N*l, NINCV(J); K=1,NCOMP) for subsequent use by the main program.
Here, the index J denotes source. For any given time increment, STRAC
must be called M times to compute concentrations of plumes from M
different source locations.
44
-------
Subroutine INTGRT. Subroutine INTGRT serves the important function of
numerically approximating solutions to the double integral equation (5).
It is interrogated by Subroutine STRAC by the call
CALL INTGRT (OPTION,QQX,DX,TIME,SIGY,SIGZ,M,N,NCOMP),
where QQX is the array of current virtual source strengths and the remain-
ing symbols have been defined previously. The primary variable returned
to the calling program is the generation of material by chemical reaction
GAIN(K), K»1,NCOMP, which is communicated via the blank common route.
Two major paths exist in INTGRT/ corresponding to the bivariate-normal
and limited mixing-height models; these are followed with OPTION set
equal to 1 or 2, respectively. The paths are shown in the flow diagram
in Figure 11, which, because of the involved nature of the subroutine,
employs a somewhat more detailed method for portraying the nested loops.
With OPTION»2, the program enters into a nested loop which starts by
setting the first value of z where integration over y is to be performed.
This is accomplished by choosing the first value of £^ shown in Table 2,
and scaling according to equation (15) or (16). The program proceeds to
the next loop which scales the y values along the path of integration
in a similar manner. This is done sequentially, using equation (2) to
calculate concentrations C(K) of the NCOMP components at the given (x,y)
location, and then calculating the corresponding reaction rates by the
subroutine call: CALL RATE (C,FUNCT,M,N,TIME). The two-dimensional array
FUNCT returned by RATE corresponds to the reaction rate (moles/cm3 sec)
of each component at each y position at fixed z. When all of these
values have been computed for a given z-position they are integrated over
y (cf equation (5)) sequentially for each component. This is accomplished
in a DO loop, which indexes K, stores the values of FUNCT in a one-
dimensional dummy array FUN, and then calls the Subroutine LAGURE, which
performs the actual Gauss-Laguerre numerical integration of the function
with base values FUN. The integrated values ZINT are descaled according
to equation (13) and stored in the array ZI(I,K). This array in inte-
grated by a second call of subroutine LAGURE, the results are descaled
according to equation (16), and the desired values GAIN(K) are calculated.
In the event that limited mixing-height calculations are specified
(OPTION=2) the subroutine performs integrations in a manner similar to
that described above. The differences between these two paths are two-
fold. First, the limited mixing-height equation (3) is utilized for
OPTION*2, rather than its counterpart equation (2). Second, the numerical
integration process is conducted in the y-direction only, since the assumed
uniformity of concentration with height allows a formal integration to be
performed in the z-direction.
45
------- An error occurred while trying to OCR this image.
-------
f. Subroutines RATE, WSHOUT, and DRYDEP
The subroutines RATE, WSHOUT, and DRYDEP are written, in concordance
with the discussion on pages 19 to 22, to return respective values of
conversion rate, washout coefficient, and deposition velocity when
interrogated by Subroutines STRAC and INTGRT. The three subroutines
are considered to be modular, and are expected to be replaced as appro-
priate by the user for investigation of systems other than the example
given here.
As indicated in equation (6) Subroutine RATE calculates the microscopic
rate of chemical reaction at a given point in the plume where the con-
centrations are denoted by the vector C(K). One call of the subroutine
results in generation of reaction rates for all components in the system,
which are denoted by FUNCT (J,K). Here K is the component index and J
pertains to crosswind location of the point of current calculation.
The Subroutine RATE shown in Table 3, as discussed previously, is limited
to two components (ostensibly SO. and sulfate). Their conversion rates
are assumed to be first order in SO. concentration, and have the rate
constant
k - 5.56 x 10~6 sec"1
as indicated in Table 3. Thus the conversion rate for SO (K=l) is given
by
FUNCT(J,1) = - 5.56 E - 6 * C(l) .
Since decay of SO results in formation of sulfate (K=Z), the correspond-
ing conversion rate for sulfate is
FUNCT(J,2) = - FUNCT(J,1) .
One should note that modification of this subroutine to include more complex
reaction phenomena is trivial. For example, a reaction whose decay is
second order with respect to its own concentration, first order in com-
ponent 2, and whose generation is 3/2 order in component 3 can be expressed
readily by the form
FUNCT(J,1) = k1*C(l)**2 + k2*C(2) - k3*(C(3)**1.5 .
Subroutine WSHOUT operates in a manner similar to RATE, in that it produces
washout coefficients (sec~^) for each of the plume components when interro-
gated by the calling subroutine. These coefficients are denoted by WOC(K),
47
-------
where K is the species index. At present all values of WOC are set equal
to zero, presuming a no-rain situation. Other values may be set by the
user as desired; this option may include time-variant values of the washout
coefficient, if rain intensity data are available as a function of time
and location.
Subroutine DRYDEP provides deposition velocities (cm/sec) of each plume
component, which may be time variant at the user's option. They are
denoted by VD(K) where K is the component index. For the present case the
deposition velocity for S0_ is set equal to 1.0 cm/sec while that for sul-
fate assumes the value 0.1 cm/sec, in accordance with the previous discussion.
g. Subroutine SIGMA
Subroutine SIGMA computes values of a and a in cm and stores them in
matrices SIGMY(M,N) and SIGMZ(M,N) respectively, where the index M is for
sources and the index N is for plume increments. The computations are done
in accordance with the equations in II.C.2.f.
As stated in this section, the 0 values are computed as functions of travel
distance out to 100 km and as functions of travel time thereafter. In both
cases the a values are added onto at each computation by computing a deriva-
tive at the current value of travel distance or travel time, multiplying
the derivative by an interval of distance or time, and adding to the previous
value. The derivatives with respect to travel distance are functions of
stability, indicated by a variable KSTAB, which may take on values of
1,2,3,4,5,6. The initial values for a are set in the main program as zero.
In order to compute the derivatives at the proper travel distance or travel
time, the travel distance of each plume increment, XM(M,N) and the travel
time of each plume increment PTIME(M,N) are kept track of in SIGMA.
If the user should try to modify the subroutine, some confusion may result
from the fact that the spatial derivatives, DSYDX and DSZDX are computed
from empirical equations 60 and 61, with coefficients from Table 4, are
for travel distance in meters, while the final a values are in cm. Thus
XC, the current travel distance is computed from cm quantities X and DX
by multiplying 0.01 to convert to m. However, the final
-------
i. Subroutine SAMPLE
After STRAC returns control to the main program and still within the source
loop, the Subroutine SAMPLE is called to interpolate to the intersections
of the sampling grid the concentrations calculated by STRAC for each
increment of the plume. SAMPLE also has an option to interpolate the con-
centrations to particularly specified locations on the sampling grid. The
situation facing SAMPLE is illustrated in Figure 12. Here a plume with
six increments is shown. The segmented line extending from the source,
out to the first increment of the plume, is the axis of each increment of
the plume. The two external lines are at a distance of 3a from the axis
on either side. Each crossplume line is perpendicular to the sourceward
plume axis segment. The dots along the crossplume lines indicate the
distances from the plume axis (0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0) a at which
the concentrations computed in STRAC apply. The seven positions referred
to in the description of STRAC are the axis position and the six positions
on one side. The concentrations at the six positions on the other side
are found by reflection. Concentrations have been computed by STRAC for
13 positions (including reflection) on all the crossplume lines shown.
The dots were left off the last two for clarity.
It will be noticed that six of the intersections of the sampling grid fall
within the bounds of the plume, so defined. SAMPLE locates these points
and computes applicable concentrations for each, using concentrations for
two of the points on each of the two crossplume lines defining the four-
sided polygon the intersections lie in. The concentrations assigned to
the grid intersections are stored in a matrix CHI(KP,I,J), where KP refers
to the particular chemical component, and I and J are the subscripts for
the particular grid locations.
Positions P]_ and P2 are positions of particular sampling stations at which
concentrations are to be computed. SAMPLE locates point P^ as farther than
3d from the plume axis and assigns no contribution from this source.
SAMPLE locates point P2 as nearer the plume axis of the second increment
than 3ay and computes a concentration for each chemical component using the
same interpolation method that is used for the grid intersections. These
concentrations are stored as CH(KP,L), where KP refers to the chemical
component and L to the particular sampling location.
When SAMPLE returns control to the main program, the DO-loop for the sources
has been completed for one source. During each major advection step the
loop is traversed NSOURC times during which SAMPLE superimposes the contri-
butions from plumes from different sources on the grid and at the particular
sampling locations additively. Thus, the average concentration for each
chemical component at each grid intersection over the interval DTIME that
is conveyed to output is the sum of contributions from several sources.
The same is true for each particular sampling site.
49
-------
u
•H
O
(N
H
0)
»J
cn
•H
50
-------
j. Completion of the Major Loop
After the source loop has been completed for all sources, contributions are
made for each chemical component and for each grid intersection from the
matrix CHI to the longer time averaging matrices CHISAV and CHILAV. The
contributions are added as exposures, i.e., time integrated concentrations,
by the statements
CHISAV(KP,I,J) - CHISAV(KP,I,J) + CHI(KP,I,J)*DTIME
CHILAV(KP,I,J) - CHILAV(KP,I,J) + CHI(KP,I,J)*DTIME ,
where both matrices are calculated in mole-hr/cm , KP indicates the chemical
component, I and J are the subscripts for the grid intersection location. As
indicated in the discussion of program output (section III.c.5) CHISAV is
averaged over some short period of time (DPRINT) such as 24 hours and is
routinely printed out. Just before the print-out of CHISAV, the average con-
centrations are obtained in ygm/m in Subroutine PRCHI by use of an appropriate
factor with DPRINT the denominator.
Following the print-out CHISAV is zeroed out.
CHILAV is added onto, in the same way. At the end of the run, CHILAV is
divided by the time for the entire run TLIMIT, and simultaneously converted
to ygm/m3.
Also, just before the printing of CHISAV, its elements are compared for
storage in a maximum value array CHIMAX. The statement is
CHIMAX (KP*, I, J) * AMAX1 (CHIMAX (KP, I, J), CHISAV (KP, I, J)) .
Thus, CHIMAX stores the largest value from CHISAV for each chemical component
at each grid intersection.
The concentration data, CH, for the particular sampling locations are averaged
into longer periods in just the same way, using matrices called CHSAV and
CHLAV. The maximum values over all matrices called CHSAV are stored in the
matrix CHMAX.
Following these printing and storing steps, the number of increments in
each plume is incremented by one in preparation for the next major iteration.
The last statement in the major iteration loop checks the value of TIME
against that of TLIMIT. If TIME is less than TLIMIT control skips back to
the beginning of the major loop and starts the next iteration.
k. Termination of the Run
The program steps associated with the termination of the run are executed
after TIME is checked against TLIMIT, if TIME is equal to or greater than
TLIMIT. At this point all values on the matrices CHILAV, CHIMAX, CHLAV,
51
-------
CHMAX are printed out. An end-of-file is written on the tape that has been
accumulating print-out of the matrices CHI and CH at the interval DTIME.
The above-mentioned steps are also executed if the run is aborted either
by encountering an end-of-file on the wind tape or by a plume having too
many increments on the map. These terminations are discussed in the
section on diagnostics (see III.C.6).
2. Flow Diagram (Major Functions)
A flow diagram of the major function of STRAM is given in Figure 13. The
long box on the left is the main program containing the executive and
trajectory functions. Inside the long box are several nested boxes indi-
cating the structure of the principal set of DO-loops. The major loop is
traversed at intervals DTIME from the beginning to the completion of the
run. The source loop is traversed once for each plume source. The plume
increment loops are traversed once for each increment of the plume from
a particular source. The two-step advection loop is a two-step iteration
used to find the length in grid units of the interval traversed by one
plume increment during the interval DTIME. All the subroutines are shown
as smaller rectangles with probes extending into the loop from which they
are called. For example, the principal subroutines STRAC and SAMPLE are
called from within the source loop in the main program.
The minor subroutines are shown by the smallest size of rectangle. The
functions of these subroutines are written immediately adjacent the probe.
3. Executive Control Language and Deck Setup
The following set of control cards is required to execute the STRAM program
on the UNIVAC 1110.
STRAM
EXECUTIVE CONTROL LANGUAGE(ECL)
@RUN,PRI/OPTS RUNID/CORE/TAPE/DISC,ACCT,PROJ/CLASS,RUNTIME,PAGES/CARDS,STIME
@ASG,T QUAL*RNGRDTAPE.,16N,888888 (ASSIGN RNGRD TAPE)
@USE 3,QUMj*RNGRDTAPE. (ACCESS RNGRD TAPE ON FORTRAN UNIT 3)
@ASG,T/W QUAL*STRAMTAPE.,16N,777777 (ASSIGN STRAM OUTPUT TAPE)
@USE 2,QUAL*STRAMTAPE. (WRITE STRAM TAPE ON FORTRAN UNIT2)
@ASG,A QUAL*PROGFILE. (ASSIGN STRAM PROGRAM FILE)
52
-------
START
: STOP If \
MO PARAMETER
VAWESAftE I
ENCOUNTERED/
EXECUTIVE CONTROL AND
TRAJECTORY CALCULATION
KAO INPUT PARAMETERS
PRINT INPUT PARAMETERS
WRITE INPUT PARAMETERS TO TAPE
INITIALIZE VALUES
COMPUTE AND PRINT
L MAP FACTORS FOR ADVECTION GRID
i LATITUDE AND LONGITUDE FOR SAMPLING GRID
MAJOR LOOP
, COMPUTE HOUR.
DAY, MONTH. YEAR
READ STABILITIES
AND MIXING HEIGHTS
READ SOURCE
STRENGTH
COMPUTE EFFECTIVE
STACK HEIGHT
RETURN
BEGINNING
OF LOOP
SOURCE LOOP
PLUME INCREMENT
LOOP-I
COMPUTE FOR THE
ADVECTION STEP OVER
TIME INTERVAL DTIME
LOCATION OF TRAJECTORY
END POINT
MEAN WIND SPEED
. DISTANCE TRAVELED
PLUME INCREMENT
LOOP-II
REMOVE PORTIONS OF PLUME
OFT- ADVECTION GRID
TIME • TIME * DTIME
WRITE CONCENTRATION
MATRICES TO TAPE
ADO CONTRIBUTIONS INTO
LONGER TIME-AVERAGING
MATRICES
EXCERCISE OPTIONS FOR
' MATRIX PRINTOUT
TIME-LT-TLIMIT
,ENOflLE OUTPUT TAPE
PRINT FINAL CONCENTRATION MATRICES
STRAC
CALCULATE PLUME
KyANDdj
CONCENTRATIONS
SCAVENGING
DRY DEPOSITION
fOR A GIVEN POINT IN TIME
SAMPLE
DEFINE ON THE SAMPLING GRID
A QUADRILATERAL IN WHICH
SAMPLINT. Of A GIVEN PLUME •—
SEGMENT WILL BE CALCULATED
INTERPOLATE CONCENTRATIONS
FROM STRAC TO SAMPLING GRID
AND TO PARTICULAR POINTS
WITHIN QUADRILATERAL
SUPERIMPOSE PLUMES FROM THE
SEVERAL SOURCES
LEIN
Figure 13. Flow diagram of the major function of STRAM
53
-------
@XQT QUAL*PROGFILE.STRAMPGM (EXECUTE STRAM ABSOLUTE PROGRAM)
(INSERT THE DESIRED STRAM CONTROL CARDS HERE)
@FIN
@@
@@
NOTE: ALL UNDERLINED PARAMETERS TO BE SUPPLIED BY THE USER,
RUN CARD PARAMETERS, QUALIFIERS, FILE NAMES, REEL
NUMBERS AND STRAM CONTROL CARD DATA.
4. Input Description and Formats
The input variables are read in either from tape or from cards. In the
present version the only data input from a tape is the gridded wind maps.
Each information unit consists of a header record, a U-component array and
a V-component array that are identical to the output from the Random-to-
Grid subroutine described in Section III.B.5.
The read statements in STRAM corresponding to the write statements in
Random-to-Grid (see III.B.5) are
READ(3)IYR,MON,IDA,IHR,ZB,ZT,NXA,NX,NYA,NY
READ(3)((U(I,J),I-NXA,NX),J=NYA,NY)
READ(3)((V(I,J),I=NXA,NX),J»NYA,NY)
In these statements IYR, MON, IDA, IHR are year, month, day and hour of
tape record, ZB and ZT are heights (see Table 5), NXA and NX are the I
subscripts of the initial and final values in the X-direction, and NYA
and NY are the J subscripts of the initial and final values in the
Y-direction. Actually, NXA and NYA should both be 1. NX and NY should
be the dimensions of the array.
54
-------
TABLE 7. CARD INPUT DESCRIPTION FOR STRAM
Card No.
Columns Format
Name
Definition
1
2
3
4
5
6
6
6
6
7
7
7
7
7
7
1-78
1-78
1-78
1-78
1-30
1-10
11-20
21-30
31-40
1-5
6-10
11-15
16-20
21-25
26-30
13A6
13A6
13A6
13A6
615
P10.0
P10.0
F10.0
F10.0
15
15
15
15
15
15
KXI
KXISAV
KXIMAX
KXILAV
KCOMP
TLIMIT
DTIME
DMAP
DPRINT
NSOURCE
NCOMP
LLX
LLY
MRX
MRY
Title for matrix CHI
Title for matrix CHISAV
Title for matrix CHIMAX
Title for matrix CHILAV
Vector of 6 integers, 0 or 1 to
govern output to tape for up to 6
chemical components
Total number of hours for the run
Number of hours in each basic
advection step
Number of hours between successive
wind maps
Averaging interval (hours) for
routine printout of concentrations
Number of source locations
Number of chemical components under
analysis
These four integers are the X and Y
coordinates on the 47 by 47 NMC grid
(origin at the north pole) of the
lower left and upper right corners
of the advection grid
7
7
7
7
41-45
46-50
51-55
56-60
15
15
15
15
IHR
IDAY
IMON
HER
7
7
31-35 15 IDIV Number of advection grid linear units
per NMC grid linear unit
36-40 15 IMAX Maximum number of increments that can
be considered in one plume
These four integers are the initial
time of the run—hour, day, month,
year
61-65 15 KPRINT Option, equal 0 or 1, to govern print-
out of concentrations at interval
DTIME
66-70 15 NRE Number of nongridded sampling points
71-75 15 JSTAB Option equal 0 to 1: 0 - read in
one 24-hour cycle of stabilities and
mixing heights; 1 - read in new
stabilities and mixing heights
throughout run
55
-------
TABLE 7 (continued). CARD INPUT DESCRIPTION FOR STRAW
Card No. Columns Format Name Definition
8
8
8
8
1-10
11-20
21-30
31-40
F10.0
F10.0
F10.0
F10.0
XLL
YLL
XUR
YUR
76-80 15 JSOORC Option equal 0 or 1: 0 - read in
one set of source strengths; 1 -
read in new source strengths every
interval DTIME
These four numbers are the X and Y
coordinates on the advection grid of
the lower left and upper right corners
of the sampling grid
41-50 F10.0 DIV Number of sampling grid linear units
per advection grid linear unit
1-5 15 OPTION Integer, equal 1 or 2, opting for
vertical distribution of pollutants
that is (1) Gaussian or (2) even
6-10 15 MODE Integer, equal 1,2,3,4, indicating
mode of analysis pursued in Subroutine
STRAC (see glossary)
11-20 F10.0 FACTOR Ratio by which the interval of inte-
gration along the plume axis (STRAC)
is expanded for each successive step
9
9
10
11
No. of
NCOMP
NSOURC
NRE
varies
21-30
31-40
1-60
1-60
Cards
1-78
1-10
11-20
21-30
1-10
11-2C
varies
F10.0
F10.0
6E10.0
6F10.0
INPUT
13A6
F10.0
F10.0
F10.0
F10.0
F10.0
15,
F10.0
SUBMIN
SUBMAX
CBKG
WMOLE
These two numbers are the minimum
and maximum intervals in cm used in
the integrations carried out in STRAC
Background levels of concentration
for the NCOMP components under analysis
Vigm/m3
Molecular weights for each of the
same - gin/mole
CARDS THAT VARY IN NUMBER
KPOLUT
(KP,D,
1=1,13
SORLAT(K)
SORLON(K)
H(K)
RELAT(L)
RELON(L)
MSTABV(I) ,
HHMV(I)
Names of the individual pollutants,
KP=1, NCOMP
Source latitude in degrees ,
Source longitude in degrees,
Source height in cm,
(K=l, NSOURC)
Receptor latitude in degrees,
Receptor longitude in degrees ,
(L=l, NRE)
Stability (1,2,3,4,5,6) and mixing
height (cm) prevailing for the
interval DTIME
varies
varies 8E10.3 QV(KP,K)
Source strengths for NCOMP components
at NSOURC locations for the interval
DTIME (gm/sec)
56
-------
The first four cards give the alphameric title vectors, KXI, KXIMAX, KXISAV,
KXILAV, one vector on each card. They are the titles to be printed on the
output pages for the matrices CHI, CHIMAX, CHISAV, and CHIIAV, respectively.
These matrices (each) are dimensioned to output concentrations for six
chemical components on a 13 by 13 sampling grid. They are defined as follows:
CHI is the matrix of concentrations averaged over the basic advec-
tion interval DTIME, which is the smallest time inteval considered
in the analysis.
CHISAV is the matrix of concentrations averaged over a longer
interval DPRINT. The matrix is zeroed out after each print-out.
Each print-out of this matrix is calculated by averaging DPRINT/
DTIME contributions from matrix CHI.
CHIMAX is the maximum number occurring in each element of the
CHISAV matrix when all the matrices of this kind calculated
during a run are considered.
CHILAV is the matrix of concentrations averaged over the entire
run.
KCOMP is a vector of 6 integers, either 0 or 1. The purpose of
this vector is to restrict the number of chemical components,
for which the concentrations, CHI are written on output tape.
For each nonzero value read into the vector, a corresponding
chemical component will be written on the output tape. For
example, if the first and fourth members are the only nonzero
members in the vector, the first and fourth chemical components
will be written to output tape.
The next card is a. series of floating point time parameters. All are given
in hours.
TLIMIT is the total number of hours for the run.
DTIME is the number of hours in each basic advection step. A
typical value is 1. This is also the minimum value compatible
with the dimensions and use of MSTABV.
DMAP is the number of hours between successive wind maps. A
typical value is 12.
DPRINT is the arbitrary number of hours chosen for averaging
concentrations for the routine print-out. A typical value is 24.
If the suggested values are read in, the program will run through the loop
that calculates the trajectories and all the other data TLIMIT times, cal-
culating hourly concentrations each time. There is an option to suppress
the print-out of all these hourly concentrations, but not the 24 hour
57
-------
concentrations, which are routinely printed out. Also, every 12 hours,
new maps of the gridded horizontal wind components are read in.
The next card read in consists entirely of integer parameters. NSOURC is
the number of source locations. It must not be greater than 10. NCOMP
is the number of chemical components under analysis. It must not be greater
than 6. The next 4 numbers are X and Y coordinates in the NMC grid
system of the lower left and upper right corners of the rectangular portion
of the NMC grid chosen for the advection grid. LLX and LLY are the X and
Y coordinates, respectively, of the lower left corner. MRX and MRY are
the X and Y coordinates, respectively, of the upper right corner. If the
rectangle is over the United States, LLY and MRY will both be negative,
and LLX will most likely be negative. All 4 of these numbers are
found by counting grid intersections from the pole on the NMC map (see
III.C.I). IDIV is the number of divisions made of each NMC grid linear
unit to create the final advection grid unit. If the advection grid is
to have squares with sides half the length of the original NMC squares,
IDIV is 2. If no further subdivision is desired, IDIV is 1.
The advection grid is limited to 17 intersections in the horizontal direc-
tion and 13 intersections in the vertical direction.
IMAX is the maximum number of increments that can be considered in one plume.
As the program is now dimensioned, this number must not be greater than
100.
IHR, IDAY, MON, AND HER are the initial time for the run—the hour, day,
month, and year (minus 1900) as given on the wind tape.
KPRINT is an option governing the printout of hourly average concentration
values. If KPRINT is 0, these are omitted. If KPRINT is 1, they are printed
out.
NRE indicates the number of particular sampling points that will be included.
As the program is now dimensioned, this number cannot be greater than 10.
If NRE is zero, the sampling is restricted to grid intersections.
JSTAB and JSOURC are options to curtail the reading of stabilities and mixing
heights (JSTAB) and source strengths (JSOURC). If JSTAB = 0, reading one
24-hour cycle of stabilities and mixing heights suffices. This cycle is
repeated for every 24-hour period of the analysis. There are 24/DTIME
cards, each with one stability value (the integers 1-6 stand for the six
Pasquill stability categories A-F) and one mixing height value in cm. If
JSTAB is not equal to zero, new values must be read in at interval DTIME
throughout the run. If OPTION = 1, the mixing heights may be omitted.
If JSOURC is zero, one set of source strengths for the NCOMP components at
the NSOURC locations will suffice. Otherwise, a new set must be read in at
each interval DTIME.
58
-------
These options are made available to curb excessive requirements for card
reading inherent in the present version of the program. Ideally these
data would be available from tapes on a real time basis throughout* the
run, but the provision of such data is beyond the scope of the present
effort.
The next card read in consists of floating point numbers that locate the
sampling grid on the advection grid. XLL and YLL are the X and Y coordinates
in the system of the advection grid of the lower left corner of the sampling
grid. XUR and YUR are the X and Y coordinates in the system of the advec-
tion grid of the upper right corner of the sampling grid. Since the origin
of the advection grid is at the lower left corner, all four of these numbers
are necessarily non-negative. DIV is the number of divisions made of each
advection grid linear unit to create the sampling grid unit. If the sampling
grid is to have squares with sides half the length of the advection grid
squares, DIV is 2.0. If no further subdivision is desired, DIV is 1.0. The
sampling grid is limited to 13 intersections in each direction.
The next card read in contains parameters used by Subroutine STRAC. OPTION
is an integer equal to 1 or 2—1 if STRAC is to calculate Gaussian vertical
distributions, and 2 if a limited mixing-height is imposed instead.
MODE dictates the calculation sequence in STRAC as indicated in the glossary.
FACTOR is the ratio by which the interval of integration along the plume
axis is expanded for each successive step. SUBMIN and SUBMAX are the minimum
and maximum intervals in cm used in the integrations carried out in STRAC.
The initial interval of integration is SUBMIN, following which FACTOR is
used to graduate the interval to SUBMAX (or DS(M,N), whichever is smaller,
see Ill.C.l.e). The relation between FACTOR, SUBMIN and SUBMAX just described
is reiterated for each increment of each plume every time STRAC is called.
The next card gives the background concentrations for the NCOMP components
in ygm/m .
From this point on the number of cards read in varies in accordance with pre-
viously entered input parameters. The first set of such cards is the alpha-
meric title matrix for the pollutants, KPOLUT. On each card is read
KPOLUT(KP,I),1 = 1,13, which is the name of one pollutant. Since the
format is 13A6, the name may occupy as many as 78 places. The other index
KP varies from 1 to NCOMP, there being NCOMP of these cards.
The next set of cards read in gives the latitudes and longitudes and stack
heights of the NSOURC source locations. There is one card for each source
location. Each card has for one source latitude, longitude, and the stack
height in cm.
The next set of cards gives the latitudes and longitudes of the NRE particu-
lar sampling sites. Each card has the latitude and longitude of one such
site.
59
-------
The remainder of the input cards are for stabilities, MSTABV(I), mixing
heights HHMV(I), and source strengths QV(KP,K). As explained above, if
JSTAB is 0, one 24-hour cycle of stabilities and mixing heights will
suffice for a run of arbitrary length since the cycle is repeated for
each 24-hour period. Also, if JSOURC is 0, one set of source strength
values will suffice for the entire run. Detailed discussion of these
inputs has already been given in the discussion of the options JSTAB and
JSOURC. The necessary orders of the cards when JSTAB and JSOURC are both
zero is illustrated in Section IV.
5. Output Description and Formats
There are two output files, one to line printer and one to tape. The out-
puts will be discussed in the order of their occurrence within the program.
The write statements to tape are unformatted. The number of successively
written output tapes per run will depend on the length of the run and will
have to be determined by the user. Examples of output to line printer are
given in Section IV.
The first printing consists of most of the input parameters just described.
The formats are arranged so that the numbers are in place within descrip-
tive sentences.
Immediately following this output, a Subroutine DIAGNO is called, wherein
several of the input parameter values are tested for reasonableness. If
one of these values is out of range, it is printed out along with its name,
and the program is stopped. The output of this subroutine is described
in detail in III.C.6.
If the Subroutine DIAGNO is executed without stopping the run, the input
parameters are written to tape in three unformatted write statements.
They are
First write statement to tape:
NSOURC, NCOMP, LLX, LLY, MRX, MRY, IDIV, IHR,
IDAY, IMON, HER, OPTION, MODE
Second write statement to tape:
TLIMIT, DTIME, DMAP, XLL, YLL, XUR, YUR, DIV,
DM, DMS, (XSOURC(K), YSOURC(K), H(K), K= 1,10)
Third write statement to tape:
NRE, RELAT, RELON.
All these items except DM, DMS, XSOURC and YSOURC are read-in parameters
and are defined in the previous section on the input parameters. DM is the
basic distance in centimeters between consecutive grid intersections on the
advection grid. DMS is the same for the sampling grid. The derivation of
the corresponding actual scales on the grids as functions of latitude is
60
-------
discussed in the earlier descriptive section on the entire code. XSOURC
and YSOURC are X and Y coordinates of source positions on the advection
grid.
The latitudes and longitudes of the intersections of the sampling grid are
also printed out.
The remainder of the output to line printer is the information on concentra-
tions, both gridded and at particular sites.
The routine print-out of concentration matrices is one set, containing one
matrix for each chemical component, printed out for each DPRINT hours. The
concentrations in each matrix are averaged over DPRINT hours. Typically
DPRINT may equal 24. After each print-out, the matrix is zeroed out so
that averaging can begin for the next period of DPRINT hours.
Since each concentration matrix is limited for each chemical component to
a grid of 13 by 13 points, they are easily printed out on one page. Such
a page has
(1) a heading indicating the type of matrix (averaged over DPRINT
hours;
(2) name of the chemical component average concentrations;
(3) the time and date at the end of the averaging interval;
(4) the average concentration matrix in E format. Concentrations are
in micrograms per m at ground level;
(5) additional print-out of the largest number on the matrix.
If in addition there is sampling at particular locations, the page contains
(1) a heading so indicating;
(2) X- and Y-coordinates on sampling grid of each particular location;
(3) the average concentration at each such location.
Pages of exactly the same format are printed out for three other types of
concentration matrices. These pages differ in appearance only by the page
headings that describe what type of concentration matrix is being printed
out. At the end of the run two of these types of concentration matrices
61
-------
are printed out. One gives the maximum concentration for each chemical
component averaged over a period of DPRINT hours that occurred at each
sampling location. The other type is the average concentration for the
entire period of the run for each chemical component and at each sampling
location.
In addition there is an option to print out all the concentration matrices
for the basic advection interval DTIME (typically 1 hour). If this option
is desired, the number of pages of associated printed output should be
anticipated before running the program this way.
The concentrations averaged over DTIME hours are also written unformatted
to tape, both for the grid and for the particular locations. The chemical
component, indicated by N and the time and date, indicated by JHR, JDAY,
JMON, and JIER, are also included. The two write statements are:
WRITE(2)TIME,KP,WHOLE(KP),((CHI(KP,I,J),1=1,NSX),J=1,NSY);
IF(NRE.NE.O)WRITE(2)(CH(KP,L),L-1,NRE)
The CHI matrix is for the grid, and the CH matrix is for the particular
sampling locations. Concentrations are in gram moles/m . At the completion
of the run an end-of-file is written on this tape.
6. Diagnostics
There are a number of diagnostic messages printed out from a Subroutine
DIAGNO. This subroutine checks a number of the input parameters for
reasonable values. If one of them is found to be outside of a prescribed
range, a message is printed out and the run stops before any trajectory
or dispersion calculation is made.
The subroutine does not print descriptive statements but simply prints out
the unacceptable parameter's name, the parameter value and the acceptable
range. The user must make his own diagnosis from his knowledge of the
definitions of the input parameters, and from the brief descriptive write-
up to follow here. Diagnostics on NSOURC and NCOMP indicate that either
(1) a value has been entered that exceeds the dimensioned capacity or
(2) inadvertantly a zero or negative value has been entered. The maximum
allowable values for NSOURC and NCOMP are 10 and 6, respectively.
NX and NY are the number of intersections in the X and Y directions,
respectively, of the advection grid. These are limited by the current
dimension statements to 17 and 13, respectively. They are calculated
from the values read in of the lower left (LLX,LLY) and upper right
(MRX.MRY) corner coordinates on the NMC grid of the advection grid and
from IDTV, which indicates the subdivision of the NMC grid made to form
the advection grid. If the values exceed this capacity, or are less than
two, they are printed out, and the program is stopped.
62
-------
NSX and NSY cure the number of intersections in the horizontal and vertical
dimensions on the sampling grid. They are limited by current dimension
statements to 13 each. They are calculated by the program from the
dimensions of the advection grid and the further subdivision made of it.
If they exceed 13 or are less than 2, the values are printed out and the
program stops.
LLX and LLY are the coordinates in the NMC grid system of the lower left
corner of the advection grid. Since we use a right-handed system with
the origin at the pole, the positive Y axis pointing upward and the posi-
tive X axis point rightward, both coordinates of the lower left corner
of a rectangle over the United States are usually negative—always so for
the Y-coordinate, LLY. The diagnostics for LLX and LLY are designed to
stop a run if the specified location of the rectangle is obviously not
over the United States or Canada. Such an error could be made simply by
neglecting the negative signs. Thus, if the user wants the rectangle
over some other portion of the northern hemisphere, these diagnostic
statements must be removed from the program.
The source coordinates (XSOURC(I), YSOURC(D) are computed from inputs of
longitude and latitude, respectively. The coordinates of the particular
receptor sites (XRES(L),YRES(L)) are computed from the same type input
data. Subroutine DIAGNO checks all these. If any are found to be off
the advection grid, the index (the I or L above) of any such source and/or
receptor location is printed out and the run stopped.
While these diagnostic statements do not make the program foolproof, they
do prevent execution in case of some of the most obvious errors.
There are other cases in which a diagnostic message may be printed out
during the course of the trajectory and dispersion calculations. This
occurs if for some reason the run cannot be continued. In this case .
the run terminates by calculating average concentrations for the portion
of the run executed. An end-of-file mark is also put on the output tape.
There are two cases in which this may occur. One is if an end-of-file
mark is encountered on the wind tape. Another is if the number of incre-
ments remaining on the map for one of the plumes exceeds the dimension
(currently set at 100). If premature termination occurs for either reason,
the effective hour and date at the termination as well as the total number
of hours in the run are printed out by the line printer.
The user should be aware that for prolonged periods of light wind con-
ditions the number of increments per plume on the advection grid may
exceed 100. Possible solutions to this problem should be considered in
advance.
In its present form the program occupies about 44000 core locations in the
CDC-66000 version. If the allotment of increments per plume is increased
by redimensioning, 222 more locations will be required for each additional
63
-------
increment allowed per plume. Thus, if the allotment is raised from 100 to
200, 22200 more core locations will be required, bringing the total to more
than 65000.
If this is prohibitive and no smaller allotment seems safe, one of the
following alternatives may be considered:
(1) The only way to effect a major reduction in core requirements if the
allotment of increments per plume is increased is to reduce the number
of sources.
(2) The core demands of the problem may be reduced by making the advec-
tion grid representative of a smaller area.
64
-------
SECTION IV
INSTALLATION - TEST RUN
A. RANDOM-TO-GRID TEST RUN
The following test run of the Random-to-Grid program is included for the
purpose of providing an example against which the user can make a com-
parison with a run on the computer to which he has access.
1. Input Data
a. From Tape
The data tape used to produce this example consists of pibal and rawinsonde
data reduced from ETAC Upper Air data tapes by the NOAA Air Resources
Laboratory and converted for use on a Univac 1110 computer.
b. From Cards
The content of each card of input data is tabulated below in Table 8.
TABLE 8.
Card No.
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
3
3
3
3
CARD INPUT FOR TEST RUN
Columns Format Name
1-64
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
16A4
15
15
15
15
15
15
15
15
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
15
15
15
15
15
15
15
15
15
15
15
15
ITLE
NP
NSX
NSY
NFX
NFY
KTAPE
NSTL
INT
Value (or alphameric specification)
TEST OF NAMER TAPE FORMAT
45
1
1
13
9
1
3
1
KSX(l)
KSX(2)
KSX(3)
KSX(4)
KSX(5)
KSX(6)
KSX(7)
KSX(8)
KSX(9)
KSX(IO)
KSX(ll)
KSXU2)
1
1
1
1
1
1
1
1
1
1
1
1
65
-------
TABLE 8 (continued). CARD INPUT FOR TEST RUN
Card No. Columns Format Name Value (or alphameric specification)
3
3
3
3
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
5
5
5
5
5
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
61-65
66-70
71-75
76-80
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
1-10
11-20
21-30
31-40
41-50
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
F10.0
F10.0
F10.0
F10.0
F10.0
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
KSX(13)
KSX(14)
KSX(15)
KSX(16)
MFX(l)
MFX(2)
MFX(3)
MFX(4)
MFX(5)
MFX(6)
MFX(7)
MFX(8)
MFX(9)
MFX(IO)
MFX(ll)
MFX(12)
MFX(13)
MFX(14)
MFX(15)
MFX(16)
RCH
XNMC
YNMC
RFC
THET
KSTA(l)
KSTA(2)
KSTA(3)
KSTA(4)
KSTA(5)
KSTA(6)
KSTA(7)
KSTA(8)
KSTA(9)
KSTA(IO)
KSTA(ll)
KSTA(12)
KSTA(13)
KSTA(14)
KSTA(15)
KSTAU6)
1
1
1
1
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
2.0
21.0
8.0
.5
0.0
2228
2304
2311
2317
2327
2340
2349
2355
2357
2402
2403
2408
2414
2425
2429
2433
66
-------
TABLE 8 (continued). CARD INPUT FOR TEST RUN
Card No.
Columns Format
Name
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
KSTA(17)
KSTA(18)
KSTA(19)
KSTA(20)
KSTAC21)
KSTA(22)
KSTA(23)
KSTA(24)
KSTA(25)
KSTA(26)
KSTA(27)
KSTA(28)
KSTA(29)
KSTA(30)
KSTA(31)
KSTA(32)
KSTA(33)
KSTA(34)
KSTA(35)
KSTA(36)
KSTA(37)
KSTA(38)
KSTA(39)
KSTA(40)
KSTA(41)
KSTA(42)
KSTA(43)
KSTA(44)
KSTA(45)
X(l)
X(2)
X(3)
X(4)
X(5)
X(6)
X(7)
X(8)
X(9)
X(10)
X(ll)
X(12)
X(13)
X(14)
X(15)
X(16)
2451
2456
2518
2520
2528
2532
2534
2553
2606
2627
2637
2645
2654
2655
2701
2712
2722
2734
2747
2811
2826
2836
2848
2852
2853
4399
4486
4494
4695
33.57
35.27
33.95
36.08
36.25
34.73
36.88
34.60
35.23
37.85
38.98
39.88
38.37
38.37
39.87
38.65
Value (or alphameric specification)
67
-------
TABLE 8 (continued). CARD INPUT FOR TEST RUN
Card No.
Columns Format
Name
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
8
8
8
8
8
8
8
8
8
8
8
8
3
8
8
8
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
X(17)
X(18)
X(19)
X(20)
X(21)
X(22)
X(23)
X(24)
X(25)
X(26)
X(27)
X(28)
X(28)
X(30)
X(31)
X(32)
X(33)
X(34)
X(35)
X(36)
X(37)
X(38)
X(39)
X(40)
X(41)
X(42)
X(43)
X(44)
X(45)
Y(l)
Y(2)
Y(3)
Y(4)
Y(5)
Y(6)
Y(7)
Y(8)
Y(9)
Y(10)
Y(ll)
Y(12)
Y(13)
Y(14)
Y(15)
Y(16)
37.77
39.07
42.75
40.53
42.93
40.67
41.78
41.37
43.65
45.47
42.97
44.48
44.38
45.58
45.83
46.87
46.37
46.47
48.57
50.22
53.20
51.27
53.83
49.90
49.82
43.72
40.65
41.67
35.68
86.75
75.55
83.32
79.95
86.57
92.23
93.90
98.40
97.47
75.48
77.47
75.25
81.60
82.55
84.12
88.97
Value (or alphameric specification)
68
-------
TABLE 8 (continued)
Card No. Columns
CARD INPUT FOR TEST RUN
Format Name Value (or alphameric specification)^
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
F5.2
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
Y(17)
Y(18)
Y<19)
Y(20)
Y(21)
Y(22)
Y(23)
Y(24)
Y(25)
Y(26)
Y(27)
Y(28)
Y(29)
Y(30)
Y(31)
Y(32)
Y(33)
Y(34)
Y(35)
Y(36)
Y(37)
Y(38)
Y(39)
Y(40)
Y(41)
Y(42)
Y(43)
Y(44)
Y(45)
NAMSI(l)
NAMSK2)
NAMSK3)
NAMSK4)
NAMSK5)
NAMSI(6)
NAMSK7)
NAMSK8)
NAMSK9)
NAMSI(IO)
NAMSI(ll)
NAMSK12)
NAMSK13)
NAMSK14)
NAMSK15)
NAMSK16)
99.97
95.63
73.80
80.23
78.73
89.68
87.75
96.02
70.32
73.75
83.73
88.13
98.22
94.07
66.43
68.02
75.98
84.37
93.38
66.27
70.90
80.65
89.87
97.23
99.65
65.25
73.78
69.97
75.90
2228
2304
2311
2317
2327
2340
2349
2355
2357
2402
2403
2408
2414
2425
2429
2433
69
-------
TABLE 8 (continued) . CARD INPUT FOR TEST RUN
Card No.
Columns Format
Name
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
10
10
10
10
10
10
10
10
10
10
10
10
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
1-4
5-8
9-12
13-16
17-20
21-24
25-28
29-32
33-36
37-40
41-44
45-48
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
IX, A4
A4
A4
A4
A4
A4
A4
A4
A4
A4
A4
A4
A4
NAMSI(17)
NAMSK18)
NAMSK19)
NAMSI(20)
NAMSK21)
NAMSK22)
NAMSI(23)
NAMSI(24)
NAMSI(25)
NAMSK26)
NAMSK27)
NAMSK28)
NAMSK29)
NAMSK30)
NAMSK31)
NAMSK32)
NAMSK33)
NAMSK34)
NAMSI(35)
NAMSK36)
NAMSK37)
NAMSK38)
NAMSK39)
NAMSI(40)
NAHSK41)
NAMSK42)
NAMSI(43)
NAMSI(44)
NAMSK45)
KMON(l)
KMON(2)
KMON(3)
KMON(4)
KMON(5)
KMON(6)
KMON(7)
KMON(8)
KMON(9)
KMON(IO)
KMON(ll)
KMON(12)
2451
2456
2518
2520
2528
2532
2534
2553
2606
2627
2637
2645
2654
2655
2701
2712
2722
2734
2747
2811
2826
2836
2848
2852
2853
4399
4486
4494
4695
JAN
FEE
MAR
APR
MAY
JUN
JUL
AUG
SEP
OCT
NOV
DEC
Value (or alphameric specification)
70
-------
TABLE 8 (continued). CARD INPUT FOR TEST RUN
Card No. Columns Format Name Value (or alphameric specification)
11
11
11
11
11
11
11
11
11
11
12
12
13
13
13
13
13
1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
1-5
6-10
1-10
11-20
21-30
31-40
41-45
15
15
15
15
15
15
15
15
15
15
F5.0
F5.0
F10.1
F10.1
F10.1
F10.0
15
IYS
MS
IDS
IMS
IYL
ML
IDL
IHL
KNTHR
KREN
ZB
ZT
FLATMN
FLATMX
FLONMN
FLONMX
IDT
74
4
5
0
74
4
9
0
0
0
100.
1000.
33.0
54.0
65.0
100.0
12
2. Output Data
The resulting output is shown on the following pages.
71
-------
1 o
ia
IN
la.
i"
IB
' in
i on
IN
! *
; N
'in
2
• N
, CO
o
*
' an
ju?
i
' r- co
M (M (M tM
"1 in f- ao ^
N M N (SI »
» ,\! LI U^ 31
. z> fo in (si* en
m in. to ca a*
M N N Njar
:sss^s
N U1 U2 CO' »
N N M f* f
in r~ 01
Si^R
MM,*
en ay
« •
I 00'
a- tn:
• «,
N 3
a u
Mj I
11
(M
Ni H » .10 X
II .— ~
a
•
in
. • a
a ••
.a M
a fo
s
72
-------
a*
S
J (M
8T
o \a a -> _) u
i r» o a a
O rt U> IfJ ,1 UJ
J j> » ?> vj P.
91 iH 91 O (t .6J
I i-* *H 91 91
I ' 9> 91 91 I I
O
'7
38
QM
a
j
.91 U) I
91
91
91 91
•-I £ i-l N » M J
rt 91 -4 rt I ,-(
91 I
GO 3'/> M i
y> n -i »
91 «4 f-4
Sfj »2 1-1 u> is-.ui « «-• t-t ji N.C y c^j
f-i >H IH ;T i t4 i 3 1-4 t*
• 01 i i i i ' i
N 3 tfl IO
| >4 M |
II
3 ^ . . • • ••••• •
'Z ZO^-flN 39 •* 31 N •• » -O
i* S 2 Q
r >. <
3>* w
m _i
J J 71 ul » N J J>
UJU1C4IA P»»oC3^-SOf» '-O Nl f*r**O SI'OSJMPJ ^33).
Q^oro9ixi«TM9»y*»aEuiaj-^x?i9i-j>a:®>J'^^ic3i*itli-iN»-*'
.n r» »*i i/> r»
i JaJN^J 71
>-in'O9i>-j):j-op«i/>^->"Mv»v»9'in>'vioao9ijc3>-30»»ja)io»-ojJjinr*'.')oi
a ** "-» .'~r*j j 'j j"r>" ^J. "°P". "* ^
> > >
o o o
3> S CD
I
3 • • • —
^ U1 ^O ^ O ^
1-1 ^ "II 1 -^ .1
o
3j _
<* « « « « , « ) , II
^!QS5ns93B^e?^S^i3S:5yi«2^jl!5lo-3 .i?3?3"iS5 =>'» !S si ^ =5 3 4 & .3 w
• Z • ••«.• *Z«Z«Z**«Z •• ••••Z«f •• •« ••.•• cc^** • z
^1 r^ A O CO O 91 Cf 30 P* C^ O «^J —J ^ O ^4 ^0 C' ^ ^ ^1 ^1 *O 'J' ^ 09 C "1 ^ ^ ' 9^ *O 1-4 c
« < «
2 2 -3
f- 3 91 I
(M » » I
^ *4 P4 iW ^ (M M M Pa N CJ .J (M ^ - «,r»
i I l 'i l I
I/I
h.
cc
73
------- An error occurred while trying to OCR this image.
-------
3 I
33
_ o;
K«
UJ »
>•
a
O.M
'l
c
.
3
3
a
" 3 *.
u4inu>i:
5 i^ in'f
H>U)in'
w»f~a
CD
§
3 - » . a -3
• M M •» ,•
'•N «
ui
,
>
io-Hin o
JT a- -o •-<
te
^
'^a-a> j
»ir» 3
75
-------
1 1
a a \
m ) V
• • • •
CM U9 '91
*4
j °3 =J ^ =5
I ^s9s
! ' 91 I 01
: «v » i-» CM
i _d • _» •
KIIOS-MMMnPOM
Ol fO ui .N
. _. _i <-ifler>»y»in
I »* "J4 • ,j a «• 3- a u
I^^T w»JiM'0"-vi-T»» rlarijoooi
ul on 0*. u)
ar-r-o: i 'in^aaoioj
I M M M » 10 » U1 .
30 N Ni
I I -I
• • • • • t
A CM in • » A J -J
i .; ; •-Irti-trHHHH rt H' H H «M N
rt 01 U) J> .
I ocMuicM • enM^oHa-tsoiucj i 5'-5c?*Jf2*r1'-?
I WHU3J1 ' r»iflii=0'^J 01 91 -a
! ' '
iiflCMMaicjjf-iuicM nmio. oiaHH9'ior'CMooo>C3C3
ii-l! ' i-trtr^H •-t;HI-(CM'J
ui » (si r-. u> oo J
r4 r-t '4
^ -I ». *. ' v> M o o ^ ^ ^ ^ o ^ -1 * •". "Z •« ^ *. i
J1 I 4 J> • l|T»J13131^q 31»»»-4^-0-.0»
7^°' n-.v,v,w«^» ^ HHa*l^*CO:.T-4'^l.v> O aogjJl0^t/t^»'30CSt
• • • • 2Pn*'}vlv^lV)'v>*-'O>o 3JU)*Tt3T>Xtr» ^3jt30.'31ui37>JJ
!tf«'»^» - 2 H-! >S * J.J Jri-: £ J .-si -J,d J .-oi -?»
Z W Olll D. ^M^r^-CIIII
•u *>' f~ r- J1 » j» v-j i» .31 J » -( J M PO ..1 l>- -O >« l~
a. rt-n'i- (T'-rtNcn'^airjT (oMaiOL3
r ii < a ..!.•.<>• a €••<>.•..
O i- . j'l.NNtNUlujf^U). J^- U. -I -I
o c. •H'^?y^'0t5'»f^ P tS5]S;r*aC:ri?;
Q. MM Mr, -o.ro «M v, z pJ^^'^^J^^ 5 ^^o^
2 ^ I I I fc: I I I
i T"'l I i I
76
-------
' i
> *> «o
o «o » o i-i
a • »* *
N:
T 3?3
>• M »«
3 :
3 5^J ,
-1 Ul
3
t-
^j
X
77
-------
s;
*i°.».°. i
»"
en
p- O1 » O
la « 4 rJ
S9I
01 O
J • « •
41 91 Of 01
i 01 I, O1
i j4» ' m
i eg ui * u)
i |; '
i a- a fi »
j ertw 1«
i a « P »
i • • t •
I 01 N 01 CM
i a i ffl
9« in oj M
ov 91 I
M rs» r- 01
-< I I
ca m
•:«S J Jl
i
a a a a
.. > • •
fC 9 9
»c54
i i
f) KI M ro i
I M » 3> 1 T "1 •>
1 co ^ in ^o H
J» 3 (N N J
M N N. ao r»j m ti «i co
(M U3 9< H i-ti «B » r- '"1
in >•> i u> >JJ» H 'J >-i
So a H N! » ^. f M
C3 a- M m, a a- f~ m
il il M: I H' ,1 M r-
II I
P*>^U>NLO P*N1
•^ao^jp^N'^n iflp*
!9F3SSW|S5 Sifl
iii i i i i i
Jl-^ 3- n M MI
'f
o.
•z
s
a tn = = i =
• i • « • .' O
Jl »O Jt .M -4
CH C^ | ;h-
cn o^ <
>
o
*-* 31 a- u m p fH'
^ -4 .** 30 J1 O> 31
(-H ffl 3) ® 3*' N M iH M
j> j) in a 9- N «A r* •&
I I I I I t I II
^ CO CO
IN » 3-
N
i1*- i*- p u3 i, j. ui a- ^ .*o
I i r T i iT l i
* ;*• 3i * J> f* 5 S O3
• •***»• ••
** j> ii :»• ffi » ao ^*
I I I I I I 1 II
a-^couj-^co^ v>r*«
^ ffl M 03 PI N Irt P* *H
N M 'O P* 3 P GO -p ^3
:i
M' PI PI) .*O >*O M PI Gi
£1-4 J) 05 C1 CO '^ ''O
irf) PO PI ^ ^J 'Pi CO
'^j _t ,H rj ui eu c^ ai r* r-
te I l l 1 I I l i I
78
-------
X
V
3
» in oo
en 3\ a
77
M J LJ in CO kl in I
•O q CJ r» 03 r- r-4
°i "2 *?'
R3fTTT?
OJ^Jlt-t^,— »»j13»~-^>.N(MJ1a-3»T> O C3U)OOTIr>)>r) O
_„..,.-,*. 3Hr»ineB4rnr*^U9i'^'^#ki#Hi49* -3,o -) J -j ns ci us o j>
oe*»««» ••••••••••••••••••••• • •••••« •
oj » f. ifi m (N 9ir-*i»ti/iMvrOf-imiar~ i/>^-^-^»»>-t'o n ji3>?iHi»-N 7^
X 91iH H -4, t-tJIJ'l 0>
a
3
Z Z3">J1Na-H
* s ^ >
O O
3 5
2 5
-------
a a
« •
on an
» s i
4 o 01 q
« • « •
1«TR ;
n8«G •
«•••!
cn m en '
,a|!8
en u> a in
«•••'
r-iN 01 » i
i
vq to «o M !
•s i
' 01 en
09 a n o ;
4* • • t
at «i en i
HS '« ;
r '1 r^ cj m i^ a
a r» en 01 in, » ij, 01 o
i i
i i
HN*UI««W»» _ IfVfCJI I 7-7 Iff
o- o a a
-8. "J ^ 'B *! "3 -^
u>i/>ip. ifld->rtLA
71 -< Ol I
4.1 -( U1 ul
I I
» J Jl
* t t
t M M m ••*•) r-O' •»<*•> *O i
.a' a ;g
3 in in r^
• t • •
i^ -n o a
•ui o 3- a -~-
a. 01 <-* 01
si en 71
o
o
o
«t
N -t -t » fl '
3 J ~t IX) I
1 » 3 *> '
• J 'N! r-i 'n a*
"I JJ 71 3
rt i-l 30'%
I I I I I I I I I
3"-3T»-Oin.nJ!-13>
*OP»JJ'J1'O3'3'.^J
jj? 3*.nr»Qor*J3 ^
I I I I I , I I I I
Tl -4 rsl T> 3. *-
^i » f~ » H -n
2.S
." j r^ *o **i ^ GO r^ r^ o
i i i i i * i i i
3- *1 tn ."1. H J} T
-M 3 J J>' •< t-( 3
I I r J ij5 f» r» -0
r^» r* "0
i i i
^ j • jj ^ ,•,; ^ 4
O, I I I ' t -I I I
OiMfjilJlNNr**1 > O J> ^ J *O (Jl •-* N ^J
(siujuioarjNui'o^o r^ooDr-r-ouia^
fO , fO .*O ro
«J Ol50'J>uiU>'iT'^l^M jj_ «-*(f-4
r.f 711
Q
I
• fO N> fO M fOt'*l M m ^1
o
z
£CJr-013-rHC03-0)O O'
oiuir-ajr»(Mr^HO a.
a: ......... cc
in cj
fj rj
t- i I
80
-------
5r?-f(
3
? s ? ? T ?'? ? r ? ~
8
z z o
* 2-
r i- «
jj < w
3 »
x o
a J
.1
j
53
3
81
-------
a a !
« ^ '
in a *j a {
HSrtS i
_i "" M S3
H OB 01 01
to
I
01 ' 01
u>j CM cv| O
VI 11 «1> W I
« • * •
I r» oo fl 01 !
j : I . I !
. . _ . !
n •? j
a> tsi ~t oo j
00) (M 1
I
I I
| a 3 a 09
< t • • •
I 9> cvj 91 09
! » 8M
0v ' Si T
01 01
«-4 01 a to
,N in o i
t rl m PJ M M M ro
t
i *-» 1,1
i a- oo r- oo | ',
' 3 t3 rt £ i i&a•**>*!
j H H I I 1
; C3.CJC3C3 ! ^^^n.rt^^vi^v,
! _ t
I Or Ol Si 1H
P S S W
•. ••• f^^jr^3"''TM'inro
3- -T ^ rH ,
•••••••••»
OOuO (Jl U1 ^j 00 f"» V) OT
J> H M N (N M fs*
^ ^ iH r-» r-4 *-*
i i^ ^ aj 3J n» JJ ^
^ H ,-4 r-l.
fl' • • • •{
till ; I I ! I
W'CJh-NCliOlO'ltNi1
i ? ? 7* ! I* ? "?
I
ifl-O 30 Q OH f^ fip (M irt
I I I I li I I I I
N:rr»ooH'UJa>aif-i
f 7 TTfif 7-CT
i i i i i 11,11
4 I I N rf tf u5 L0
III I I
r» N •••*>» 3Jjyi
»yiJJJlJ>r~J,-->-
N N ^*> '^1
I I I I
"i-o.Tif'O •^*=rr* "OuyvM
Tia-H^sr 3 » 31 o
•n^Tffl^ao'Oji.')-^
r^ yi oo N ^ -H ft .N N
W rH J3 UT CO -^ M ?J r-
jiT»^,rr* .^^-rnar
II >'•
^50 -T TI
• . . • ' 2
Tf
•r
o
ui en
o
O l/l 30
j "i .n
00 U3 .N 5T vfl tfl
ajmr^OI J»r- -4
• ••• ••••
>
a
^•> ^j r 4 ci co >
a
i (^ H 0*1 CO
£U1
^D
» ro ro n ro f*> rn
u) rHrOPOPOfsiH
Q.
X
te I I
32
-------
d
H.
^j
a
»<>MMM«M»»»a.»»»u>tntni.-iiou>iJu)k3Kr-r-r-r-e>«eoa- oj
CJ M '(M CJ M M M (M M rj N (M (SJ N (SI PJ (SI M (VJ r^ M M -N fJ (SI (~J ,"J M CJ fJ CJ !M «• »• ffl
83
-------
91
01
a
•
01
»>
SffiSS '
01 O1 •
r. 5) PI H ,
in o co «o
38*8
01 SI
a M a f-
(31 N 01
01 <-« 71
01 on
a » a
01 (S| O1
01 1 J1
rl 9 T
I I I I I- I I I I
m in co
* I
CO J-> sj 0)
r. ul N 1*1 r-« 03 uj o in
» i i i i • T T r i
a> N Q »
• • • •
o» m 01 «-i
t»t««>««-««
I I i I I I I I I
0* 01
CO rH 1-4 lH
» • • • PI PI ,
ui co rH o
t-i •
i
r- oo Ji r-
• • • • PI f k
» in H
»• rH a- »
a ^ o r-
• • • • •O"OP>W1P1P1P1P1P)
Sh f o tfi '
(-4 91 |
it IM ffl .N .
. . • • MP»p>~1P1P)P1MrO
Ol rj » Ifl
H -4 I
I
G
.•>r33-i*1^3-3
* • !••••••
'f
i f J
i 1 I
.n o m co r- *•» (st fst
H t |
i^-J-O^Nul^f*1-
• ••• m>o'nfO'io»*'*i'o
i -* i
g
^i » a-^ijif^r.rsjy>
-4 3 J> T» 50 H •-* i^ •-«
-3 J _i in r- al ' J
-1OV11O'OJV» O"O
* a <7> J
• * • • •tOM*-»foro^'*tr*i)i
lH Jt r-t 0> ^^.
C 23 *-
Jl 31 i-(
1 -^' r-t
i— . a» 03 f :ri o oo TO
yi1^ rj ui f*--7i co^i
I I l
H-oa^-m'tf-p-^y
77 f
^ a » 13 3- r
2 • • • • O
O 0^ uO OT ^ —4
n ^1 I 0^ I V—-
r 01 01 <
CJ ^
'> 'u.
Q
S o'
o z
uj 1 -O .3 "O '* ^ -1 3 0^ J ^ >- -"* -1 %>I u"> J "O PI 2i •
a.
X
fc- III
84
-------
> » m i-t m M o i-i 0» r» »> o N »* ^ *> u> !•• ^ M *« o • «i> •"» r »• H la H la >t> ^'V> **> r a a
9
r<
i *S
§N- Ul
U1 _l
M UJ
a 3'3» iO»NAJ»3'O»'
I /> 3 "» l>- J5
VI UiUiMuir-"iur-«op-uif^r«r-
w ' S 10 fo 01 Ji M TI' » » » m H (ji ji
r^tsi tf> "> (
-» 3 ^ wi ,
a.ji r- ui •>! jn
4 "?'
'J §
* .-I
X
CO
6jnr»^w wr"i
a *H *H (W * * *n i
M M M 'PO (*1 PO fOM
CM (M N M M M N N f J
M CNJ (-4 (Si
-------
I
:5
fS ''8
«i
ma N a
• • « «
|i e
I ft I f
on OB M «r
•I ft * m
i-i ft -i-t
i
01 . 91 I
Mm CM M '
• • « •
M N 01 ft I
i 'i T • i
I r-i^in
01 CO N 00
* • • •
U CO i-l M
^1°.
Hi*""' T
o u> a 31
« • • «
i Oil H 91 N
: u» «
I Hi I
i -*!•• O> 5T
: -Jlrf U* ?J
a! a o a
.....
rofO»Mi»lM'^''ii^
••^ w> » » » 1 "» "t l
1 vD. l^ I CSJ
i ^'1^=1
Of PO ^" y
I i
j 00' M vfl. "1
I "1 J> M .«•
U1 ^ 1^ If)
rol a-
I
i u^aoTiai'^immfca
»> 10 M CM
ulHtviiuJi-uJNri
i-
oo
.•••<««•«
91 ^ Q n O CO «i
1fl
H « r» a j a
in s OK r^ in ^
III ^ HH
Q M il1 » fO H H («• .•»>
•^ » -1 I » 3- "> fO
I I I
>
3
I I I I -(
anN .oa'.njO'OD
'^r-c^in 'Jloi 31-4 .j
i i i i 3 .4
^cj^^in'so Miiji
^i/i NI f» r- H . •<
'Z ' t/1
u a ») o u> 12
'j; 91. ff (^ H. N r^ H oo as
I A1 .f. .«.••*• LA ••••.••.•
i ii I » Ji JJ !0 | H •! W| u HI Ul ^ » ' '•• 3 N M
3' I I I I 3 I I I I r* H i-l
conjoin o.a.N
-T -1
jjj. ^J f f ^) vp -H f-i a:
i I i
a:
tr
86
-------
II
o
c:
>
Mronton^nMtffvffvvjrv
f-
ct
' te
I *N N fi f-J N CM N C*4 f-J ftrf (SJ
-------
•as
01
SrJlS
eUo?
SitS
r~1 r~i •
J,.JU
"I"
f;':0.
i I
« . • • * PI PI « PI PI P» PI PI i
A «r ui Pi • ' i
HJ c\j 9i a
H H
Jd
Or- C3 ai N
NisnsS* i
H; I H I
rt ^ ^1 M
1 ! I
PI PI p» PI pi »
i I
PI PI PI PI P*PI
; j
i :
piipiPiPip^piioia-.-oi
3i?3f
a; 3 H 3
p>' cniDi 01
I i
I !
• MipiPiPiPtlmpya-cii
Hi H
uva
i\*i
3/
UI N 9*i091l0**l0 ffVUiiAU9
aswasasss
i ••••••i«»j4
1111:1 H
jiTa-»3|30.-^r-n
I a* f^' P" *JP d* PI
I lil III I
57 911 in PV 3 T* ui og
•N N! M N. f» oo a uS
r«l'»P*»9l3IHrf>HJ1
III I ! I H1 H I . I '
i ' '! ' I
HilONKMrHIIMUI-Cn r»
oonpi'HNIN Min o>
•,•••••{)«••
?, ' -! --";'-'•
•It »i • • • r « •
HiN N|H I » in.f» «
I, I I I j
ni 01 oo u) a-
HJM r-. a> oo
or CM aho PI g in PI en
tn. Q in v t0i CO (M< f oo
•<••••!• •< • •>
tn a r- a- i PI r«i oo oo
UllPI N a f-'H H'uJ Cn
PDI<«S'9-1|3J*^ -i
ui,3
U>:91
D|Q M..-
IM a
H ^)'
•O|H oa'a- H;91 PI,U1 H
p-,O » a- 3'ifl cviW H
5
1
is
10
15
b
•i • * • z
;
P» P» *P» »
1""
01 !M 3. ST j Z ' i
•i t « • ' O P> M Pll M Ml » Pl'«l P):
Oli UT ff\ '.M ; f-1 • '
01' | Cf I ' i— ' ' ; '
' • w PI PI. PI PI! PI PI PI PI'
°l
Ui
3;
«,
21
£'
Si
fe-
S'SSisslsas^i
ol'r* *rfd!4-JJ J
I I I
r r v
u> H uv r» in) P> O' a «r
l
Ol
al
•
v^'on "n in miPi oi'H INI
P>:H 1 09 ial* »>in :M
•It •'• *!• «.•_!
H.P) » m 'j>|ni
IS 'PI hft'
0} U3 3
r fs
^ d* tOivd
H. og u) op o q rsi |
I T ~\Tt H; I
l 11 I l
£'
a:i
HI
M'CSI fOCO-HMHf*-
88
-------
I I
2
a >-
31 I , : .
•-I ' I iH
3 iii
U : ' i
M
1 ' ' C5
' ii
*•• *••
* *
• fc
K- f1* CNJ,fJ C^J N N f^ "4 ^J (Sj ,N CM »N ^ -N ">l f J N f^ ?5 N C4 N N ^ »M TJ CVJ C4 N (M N N N » ff JT 2
A Z
89
-------
ffl
oi
1 l| 1
a
t 4 •
N co at N
i |i
u«f*i rt co
•i • « •
H rll*>
f, t I
137'
dm a >
01| *0 01 I
S ft
SCM£N
«i • • «
I i
01 in at to
I I
MMPIMMJMf'IMP'l
PT| »•> «f n n
lags j : =
i • • • i mt ^ 3" r
i PO eo in J
I I' I • '
n ! ' ' ,
j TT'j I*"
I ','!',
»1 (M OJ O 1
f-|N 3 (>J ! i ' '
!• v i . !
M r» UJ 91 I ! I
P»l U) -* If) ' I
,;'-*' i : I
I ^j«) r«-J T» i , !
\ *i • « * ! • i*n ^ *o *•> ro » r*v w *TI
i Ml -H r-* H ; H- , I
"!' •• 15 : : !
MI a P» a- I
tf^JTjff I
u a: H a in
jo oil * ov
•CL 9>i a> !*•
| x 9*: art , <
lo ; ' v
,u iR
O |V1 1*1 "1 1 '11 "1 "» f*> "V
3!
<;
$ !
•i./
Q cn i
mi a .
• a) a 01 r a
i M p. m u) CM
W N M ^
'I
! lA UJ
in n M
Si ia » i
rfU A\*
I | i i ' *
»;N nla
HiH rt"
n!in
TIT f'f f
u» ca
H
1
•i • «i i
NIH n r
I I I'l
;i
i!
i en f» M 01 !<• KI
ION CM,i
i I »» ui 14 N i •*» 9»
li I I; I
N-Nr>iv0ul^a>eOHi
•« l.tTSl'SIt
I CM Hi U> '
I II
ux f-i H
r?r?
i .a » >* ~t x a >* -i 3
» J» ^ f»< H *> SO •»>'
O»lkOJH:r«O.HN.
J'd,
.0 rg OM \n co HI
Sfflifl^wotaas
i i i i * i : •+
» » (SK
1 1 li
ioi^ o •*> Si H 5>; H H-
, I |f*1uir^T<
i I I I l~« I I I '
i I '
1 '
i ui oo o 'Q «ii n mi co o
i »'.n i r» r»]f- H' oo CM1
r-i w H ;
(•i i'
i
J] ^j ^ /M C3 J> (M f*»i
*'?? cj H1 -4 4 j? r~'
111,1 I '
I
, ^ r-f 01 O41 CN H CO kd
O'
to'
<'
a:
H ro r» M co'u j UN'
"i CM H« m'uJ ") « sr
ro
NNirf T.THITT
la ;
• M'M PTn MlM fO ro PT
10 :
iz { !
O
CL.
>'LD iH
•OH
m
-------
B. STRAM TEST RUN
The following test run of the STRAM program is included for the purpose of
providing an example against which the user can make a comparison with a
run on the computer to which he has access. The run will be 24 hours in
length, the basic advection and sampling interval will be 1 hour, the
short term averaging matrix CHISAV will be averaged over 12 hours, and
wind maps will be read in every 12 hours. Core requirement is 45000.
1. Input Data
a. From Tape
The data tape used to produce this example was produced by the Random-to-
Grid program using an original tape supplied by the Air Resources Labora-
tories at Silver Springs, Maryland (see III.B.I). The tape so-produced is
a standard label 9-track tape written at 1600 bpi.
b. From Cards
The contents of each card of input data is tabulated below in Table 9.
TABLE 9. CARD INPUT FOR TEST RUN
Card No. Columns Format Name
6
6
6
6
7
7
7
7
1-78
1-78
1-78
1-78
1-30
1-10
11-20
21-30
31-40
1-5
6-10
11-15
16-20
13A6
13A6
13A6
13A6
615
F10.0
F10.0
F10.0
F10.0
15
15
15
15
KXI
KXISAV
KXIMAX
KXILAV
KCOMP
TLIMIT
DTIME
DMAP
DPRINT '
NSOURC
NCOMP
LLX
LLY
Value (or alphameric specification)
CONCENTRATIONS AVERAGED OVER BASIC
ADVECTION STEP DTIME
CONCENTRATIONS AVERAGED OVER
STANDARD INTERVAL DPRINT
MAXIMUM VALUES OVER THE CHISAV
MATRICES
CONCENTRATIONS AVERAGED OVER THE
ENTIRE RUNNING TIME TLIMIT
24.
1.
12.
12.
2
2
-3
-16
91
-------
TABLE 9 (continued). CARD INPUT FOR TEST RUN
Card No.
Columns Format
Name
7
7
7
7
7
7
7
7
7
7
7
7
8
8
8
8
8
9
9
9
9
9
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
1-10
11-20
21-30
31-40
41-50
1-5
6-10
11-20
21-30
31-40
15
15
15
15
15
15
15
15
15
15
15
15
F10.0
F10.0
F10.0
F10.0
F10.0
15
15
F10.0
F10.0
F10.0
MRX
MRY
IDIV
IMAX
IHR
IDAY
IMON
HER
KPRINT
NRE
JSTAB
JSOURC
XLL
YLL
XUR
YUR
DIV
OPTION
MODE
FACTOR
SUBMIN
SUBMAX
3
-12
2
100
0
5
4
74
0
2
0
0
1.
1.
4.
4.
4.
2
4
2.
250000.
5000000
Value (or alphameric specification)
10
11
1-60 6F10.0 CBKG 5. 1.
1-20 2F10.0 WMOLE 64. 96.
12
13
14
14
14
15
15
15
16
16
17
17
18
18
19
19
19
19
1-78
1-78
1-10
11-20
21-30
1-10
11-20
21-30
1-10
11-20
1-10
11-20
1-5
6-15
1-10
11-20
21-30
31-40
13A6
13A6
F10.0
F10.0
F10.0
F10.0
F10.0
F10.0
F10.0
F10.0
F10.0
F10.0
15
F10.0
E10.0
E10.0
E10.0
E10.0
KPOLUT SULFUR DIOXIDE
SULFATE
SORLAT(l)
SORLON(l)
HS(l)
SORIAT(2)
SORLON(2)
HS(2)
RELAT(l)
RELON(l)
RELAT(2)
RELON(2)
MSTABV(l)
HHMV(l)
QV(1,1)
QV(2,1)
QV(1,2)
QV(2,2)
41.00
89.50
0.
37.00
86.00
0.
40.00
89.00
38.00
85.00
6
70000.
1.0E4
1.0E2
1.0E4
1.0E2
92
-------
TABLE 9 (continued). CARD INPUT FOR TEST RUN
Card Ho. Columns format
Name
20
20
21
21
22
22
23
23
24
24
25
25
26
26
27
27
28
28
29
29
30
30
31
31
32
32
33
33
34
34
35
35
36
36
37
37
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
15
P10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
-F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
15
F10.0
NSTABV(2)
HHMV(2)
MSTABV(3)
HHMV(3)
MSTABV(4)
HHMV(4)
MSTABV(5)
HHMV(5)
MSTABV(6)
HHMV(6)
MSTABV(7)
HHMV(7)
MSTABV(8)
HHMV(8)
MSTABV(9)
HHMV(9)
MSTABV(IO)
HHMV(IO)
NSTABV(ll)
HHMV(ll)
MSTABVU2)
HHMV(12)
MSTABV(13)
HHMVU3)
MSTABV(14)
HHMV(14)
MSTABV(15)
HHMV(15)
MSTABVU6)
HHMV(16)
MSTABV(17)
HHMV(17)
MSTABV(IS)
HHMV(18)
MSTABV(19)
HHMV(19)
6
60000.
6
50000.
6
40000.
6
30000.
6
20000.
4
40000.
4
70000.
4
90000.
3
120000.
3
140000.
3
150000.
3
160000.
2
170000.
2
180000.
3
180000.
3
170000.
4
160000.
5
140000.
Value (or alphameric specification)
93
-------
TABLE 9 (continued). CARD INPUT FOR TEST RUN
Card No. Columns Format Name Value (or alphameric specification)
38
38
39
39
40
40
41
41
42
42
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
1-5
6-15
IS
F10.0
15
F10.0
IS
F10.0
IS
F10.0
IS
F10.0
MSTABV(20)
HHMV(20)
MSTABV(21)
HHMV(21)
MSTABV(22)
HHMV(22)
MSTABV(23)
HHMVC23)
MSTABV(24)
HHMV(24)
5
120000.
5
100000.
5
90000.
5
80000.
5
70000.
Cards 17 through 41 simulate with card input the more desirable situation
in which there would be hourly input from one tape of stability MSTABV
(ITIME) and mixing height HHMV(ITIME) and hourly input from another tape
of source strength QV(KP,K) for each component KP at each source K. The
particular numbers in MSTABV and HHMV in this example are arbitrary but
typical of conditions in the temperate latitudes under clear skies, if
one assumes a sinusoidal cycle of mixing height as is commonly advocated.9
Other parameter values are input by the code. They are
1. The rate constants:
FUNCT(J,1) = -5.56E-6*C(1)
FUNCT(J,2) - -FUNCT(J,1)
2. The deposition velocities are given by:
VD(1) » 1.0
VD(2) = 0.1
in cm sec
94
-------
3. The washout coefficients are given by:
WOC(l) - 0.
WOC(2) - 0.
2. Output Data
The resulting output is shown in the following pages.
95
-------
£ S
| o
•» •
«
o
o 3 a o
— 1/1 a. at
IU •*
Z OC •" -J
iu Z a.
$
a ••
z o
u. « u ae
u ac ae •- at ac
O 0.* Z 8
Z UJ h- U UJ O
« at j — . ac —
< v» »l- •- «c ac
a 3 « xui o o o
vi ac Zufr- ac at • • ac z
3 xx on»- ac ^ uj ac -< z _j
O O »Z 0«U|0«I v» o t-
. < — ac ac ee « t^i—o«
i ^- -j 3 at a. i- >o ,r 0*0
a *— x z •o»r**ac
ac w>«oa a -.-»
[ V) * * (^ A UL O X JC Q X O * * t* O ^ 4C
i iu vi <*> P^ «^ uj Z O u» Q • O O ^ ^ *^ 4 O *"* *^ O
i UJ 4 • oC OC >- X *- *M J U. ^ _j m — O 3
10 z o uj t- z x Q. — o ac — a
— jtxo^»-'>-xwz ac o ac r
EX o or^aoz^uj>a9 •ur^Oacdoacw ac >• ac Q o us a.
_ o u. uj C > < a. c Z -j o a. XCC4
v> 3 x H* z (*>*_jvj>t-z z vi ^ a 3 z o«v>
>• a z < ^ o a x o u 9 >^ o uj z •• w < uj ou>
-i x*** • -J • x Z — "J o a> — ac — *- — X Z f\j uj uj
O uj 4 O ^ Z X Z X UJ —1 ^» *t 9 ^— *N O t/i ^> o
^ ^ UJ ^ O 3 ^ 3 •• ** *• 3
ac XXO uj « o O < ac z O uj H-
uj ujv-»oooujz< o * «n ui o vi — as *** — u>
OC X iU ^ wi hO<-tO< X
UJ O-JtJUJUJUJuJUJ^-UJ UJ«£ UJ ^ UJ UJ UJ Q,
x atxxxxxxr x — ^ xis ri r <
96
-------
.*> » -o m •» —« o OB AI m *r«« r- « o m AI *« m • OB m »•• AI m
m o — o r-o •» » o » « » AI OB » • m » •• r- r-r- •» r-^ o
» «of "•*f O O" Ol*t Of d*i
-« .* * * oo * « *<» <*
•»» mm m .» f- r- »•*
mm «* m r-m m »r o» «f .o ** Aim OB m m AI ••• AJ r- *>• m ^ o <••
»« .«.*o • « «M •« *r
•» •• eO f* O «
•«m — rf> o & o
•r " •a •o » — o« m «m -• o» -r * r- » oin « •< -or-
«« om « u> m« o>« mm «m • »m o>m «m •» vm r>m r>m -o in
« •« o» mm «« «m »•• «m o>«
m1*) o*m ^m « r> « «m
o o
•« Aim mo «m AI^ m^ + «
r> m^ ^« mm AI* to^uukw
.« oo
« m o>« «o mm «r- » •• AI^ m« am — r- * — r>m
mo -^^
•«r- Or- or- or- Or- o*r- o> r- «r- •r*> r- r- r-r>
^•9 *f m -ra •* & <* in mj» -n« mat m« m« mti
00
«-l •—»
• »
*^m m r*Ai oo m9k ^m otr- m-« i-«> mr. o>r- «« Aim «m ,».» *• -t r-m
•• •• •• •• •• •« . • . • • • •• •• •• • • O Al
— • O • or- or- »r- cpr- «D f- «r- • r- r-r- r-r- r- r- «r- -.
a «m «« o » m — J -» c^r* AIO mm « •« -• o «m «r- — —
_ _• « r-m •<• -t o>r 4>m IMAI o- ^ m«« ~> O • o ^ o- o « r-ca
r>ttt r-r- r-r- J3 r-
m » m gD m« m«
y O *** •">•> «A tn f1* ^» O O *V *^ >A >«O f"" ff» O -^ *H -oO -O tf» 9* >*^ ^ "«O 4
•4 •«••-» r* o mo 9> * ^- • 4> ^«o *^ u> r*- -.f en c^ Of\ -o AJ
•/I •• •• •• »• •• •• •• •• •» •• •• • • ••
-•a* O 9- o « 9>« Ch* <^ « » « « « * « r»* r* « r- « o «
^v O «O 1,0 liA (^ -tj1 ^ n lift f^ ••* *V • «•* --.f O O O ^* 0^ W O.D O ^ <-0 I
a
^
o
z
w •• •• *• •• •• •• •• •• •• •• •• •• •
97
------- An error occurred while trying to OCR this image.
-------
ae ooaaooooooom
Ul OOOOOOOOOOO«
t- ooooooooooom
ooooooooor-
OOOOOOOOOO
ooaooooQQ*
o
o
oooooooaaooo
oaaooooooooa
oooooooooooo
ooooommoooooo
ooooommooeoao
o
OOOOO — U1OOOOOO #
ooooommoooooo o
.J
<
>
or
UJ
Z
o
at
o
z
5
oc
o
UJ
o
^
at
IU
>
< IU
* u.
a -j
— 3
^ *rt
^
at
z
IU
a
•»
p>
*
m
*»
>
oc
z
u.
a
o
z
u*
,_
«t
•^
O
O
o
•M
o
o
"
0
o
a
o
o
0
**
0
o
o
"
_
«
•*
o
u
o
0
0
o
o
0
*
.0
•*
o
CJ
0
g
o
o
o
o
o
o
o
a
a
o
o
0
0
0
o
o
0
0
o
o
0
a
o
o
o
o
o
o
o
0
o
o
0
o
0
o
0
0
o
o
0
o
o
o
o
n
o
o
o
9
O
O
O
O
O
O
O
0
000
000
090
o o o
U U 0
o o o
0 O 0
000
0 0 O
o o o
0 O O
•a
*
•"
X
fc-
X
z
h-
z
UJ
X
IU
.J
UJ
X
0 IX
ct OO O
O -a-H o
>- • • o
OC «
O o
O. 0.
X X
Q a
u. Z z <
O Q UJ
a a u
U l_l Z
a
99
-------
OOOOO — •>J--0»^-» —
ooooooooooo-ir»
-OoOOO-«-» •* «• O O 3
«
•• a
| «
a
ct u>u\u%iAiAinin^>o a.«^o^ •
< rsi r> lu u>
o u
x -a iu
•< oe
»—
u-> OO^^JOOOOOOOOO JC 3 O
oouoOo>->uuo < — —
-J ac oe
3OO
a r» x — o o
— •- z z
a u -J of ae — —
uj a < *- < j j
o — > •< a. a. a.
v)
>az* ....... •••••_ u.
<>^iAtAU>u>i/tiAirtu%tAmirvu>in zz
cf' t- f< a a
- ...... .......X OCQ'Jl-
z^u»i/%^iAir»ti%j^i*"\i/\iA^irtint2 ^-^^^
ui r -> a a uj
«j — ui _ — a a i_j
100
-------
•* in 11 ix n
«••»«••
* in o » »
« • • • •
"V «• -« «• »*
o
ooojmr-^_r-a>o
oooo»oZtf\!nnio
3
II O p O O O
3
X OOOOOO'M — «ir<
o o o o o •» >n « o o o
* <« •• -* m
•or
•* •>« ••« "• «•
O-»«»"»-««».»m«
— o
Ooogr~uri'MuiOOOoO
— — ooor-o^iXooooo
osoooo^mp-oooooo
^••••••••••••*
^T^4<^*«*^«4tn'V-4"4«4*««4*^
oo30...o
"*^ ****m^«*^«^".*^«^-**^f^ ab*^o^ .
a
z «
oomoooooooooo
— -- cioow-— — —
~ _ v^jo<-*uucjoooo
o o o ««oouoooOoo>n
c
•" * »- « -l-l
° * « Q. Q. Q.
* u,
«UJ- — — -• — -,-, — -t-.«« — _ Zz
"•Su, Z u. 0° u,
*"•=> i « «-, Z
= =; Z3UJUJQ
-30 UJ |J t— ^. M
<-<"* OoaC003OOOOOO_) — <« t-
4UIOOOOOOOOOOOOOUJ U.ZZ«
ae 0000300000000 —— ct
~ ac o o
3 < «ot z
= •>•£ - - §§ o
* Z X Z U u Z
« — c
101
------- An error occurred while trying to OCR this image.
-------
s §
o o
«r o
» —
•» o
ac ooooo«
vu ooooo»
•- OOOOO —
m -t ni ** ~* m
ooooom*Mr-r»ooo
ooooom^-mni-exoo
ooooo — -omorjooo
».«f-4.«.«>**«4
roooo — tn-J'-'-oooo
h-OOOOO-<-*-"-'OOOO
ooow>o-«>noeooo
oo
eo
ooo^^oooooooo***
oooo«v^ooooooo«> ac oo in
ooo»«niooooooo-r a «-j o
o — m»ooooooOoo ocoa
o^-'o « •" ••
O«^->uoUoooO<->Ovn _iucac
J^r^----------x i0o
acoooooooooooooz 55
oeiuooooooooooooo ae««
!u™oooooooooooooz a ^, ^
Su^Z ^^^----------"" * ZZ
>_ !_ •« a a
•n « u. Z ui tn
uiu.a ui«>^>nz
= _i Z3UJU4Q
iS0 u, o K- »- —
«>nzoooaooooooooo-J — * S *~
>UlOOOOOOOOOOQOOUI U.ZZ«
00«3<30oOOC3e30C5c3 ^ —— «
§5^J-J-J-J-----3 « ac at z
z i _i a a uj
-•Muj --CC«J
103
-------
<_> ooooom*r»ino«o
>» OOOOO«NniO
Z OOOOOWf~»o — — O
OOOOOO«-«P>INOO
OOOOOO-* — — -OOOO
ooooooa>*»3OOOO
oooorinOooooa
inmininOo^ininmininin
"
UJ
ac *
X
I&J
IU
I
H-
ac <-»
>
a
O UJ -I
IU O <
< x ac
ac a uj
o
o
o
in
o
o
o
m
0
a
0
o
o
o
in
~
-*
r^
•"«
o
O
o
o
o
3
in
o
01
in
•o
•*
0
a
o
in
IN
*
in
to
a)
in
O
a
o
a.
fB
in
o
0
o
o
o
o
m
o
o
o
in
O
o
o
o
o
o
*
o
o
o
m
0
0
o
o
o
o
m
o
o
o
in
o
o
0
o
o
3
m
3
O
O
in
O
8
0
O
0
3
in
o
0
o
in
o
o
o
o
in
•a
wl
X
ac
^
X
*
r-
-»
o
«
in
O O
o o
a<
a
ZU-C3 luac^i^Z
O — t X 3 uj uj O
•«z>a uj o ^- »— •»
c-vtZOOOOOOOOOOOOO-J _<<»~
4IUOOOOOOOOOOOOOOI a.ZZinminiAminininmiOino ••lacac^
iu z _i a a u
u»*!u —«aa^j
z z x E <-» 1-1 z
o — «« — a
W *~ JC ^1 X >- w
104
-------
in
o
o
at ooooor-<*r»m9»infMa*
if ooooo^minr-Xor-...
f- ooooooa>iM>M«p«a
IM
o
o
ooooo^>->m^o>«o«
§OOOOinr-tno>r-mor-
oopp^r^MiM^ppp
S"« O> IM O «
m IM m o >M
*» m m o o >M
* r-
• •
0> M
z •*•
•o
OOOJ-OOOO« 3 11 -. o
• •••«*•••••••• ^*««o
^**^^*I^I**«*I>^*^*^^*»^»N«-«^ Okf^a^ •
P"«moo0oooaoo
- - uoooupouw
_ OOOUpo
aeoo *«**«*-*»*»rt**«».4**M»*w* O
a x — o o
_ — »- z z
o -j ae ae •• «•
si " " -i_i
o » « a. o. o. u.
«ocoooooooopooooz zx o
atiuoOOOOOOOOOOOO oc«« iu
ui»-oooooooooooooz o <^ m
>Z>-<.... •....._ It X
«U — « — -.-.« — — «^_l — _4_ ZZ >-
H* ^ irt Q a •«
5<^ Z UJ « X
»"•= luacuiMZ
Bio XSiuuiOZ
— 3O UjO^^MlU
^-^zoqaooopooooop-i — < < >- t.
ooooooooooooo
-i a a lu
•« a a
-------
SECTION V
SOURCE PROGRAM LISTINGS
The complete listing of the Random-to-Grid and STRAM programs as used to
produce the installation test run example of Section IV is given below.
106
-------
O U O l_ ri_,UUL>CjU|jl-iULjou<-'>->UL CJ U Lj LJ I
|UUOU^UUUUUUOUUUCJU|_|UU<->U|_II_|<->U<->LJ
o a o 2
,uuuuu uuu
aurluU tiua
107
-------
s-t- • »w * L..UI
(-.T! > -• rf« II V — •—•=!'. t-l
C'C'«. "t- I' II* >-«»U.^*^*-«»*V-
*MM MK >-> MM •— Z«Z'rti--X — ">•• »"«<--1«
r^CnO «• M M» — X M U H tl • » II
K— Zli»l'i-'«»»»»-l»-:— f.
t., JT o _ < 1- i-' — . —• — — u- ».
i« 51- w*'1 -• " 2" x x - >-
««"3C-jl"l*'**» tlU.b' t>
i:^ —> -- f x x >• >• r >: t ^
oaoouooo ucjouoaucjcicj
au
o o
108
-------
ououuoauoouaouuouaoouuuciocuuuuouucuuuuuucjuuuciuauouciuacioucj
109
-------
n
z
55
x >- >-
4,1 yi Irt
i
i l' i
X
I
•3
>.
k.
X
I
I
II
— S
X _l <
*2-
iCi
, u
K _i -i — o
r> —
z -,
S|^ ?i
— r. "-i M w.
f-i . — in
. c o o
u >- s > —
Z l-t — •" K
I ^ 3
M ^ . M 5 M
h M H — f- t-
O o
t) -•;
^^
ao a>
§S
: 0 -* >-t *^ *c
u o: o: c: o L-
M s a * t> i_i
. fj
i. _i '' -
K -> =
II I ~ T.
_l >• * •-!
•>££-
41 * ^->
t- «
to u z x
z1-1 c
o M o: o
o T a-1_
W " uj
M " ro
•» ->«
>- • >• • u
t: — L M ^
* t3 *• ** *-
• *" • U> ^4
v_i y > j> t,
I- T" i« •> — I
Z I Z I UT
s^hK-^,;a
v • f.
•- CJ »->->- 1,'i V >_ <
z ij» L «: I* it* i" z: ii
H z M M Z n a.
KCM ccft;o"0£O
c.n-)t.urj->c.u.
Ofl
c. •
esr-
J^ K<
iO CJ
U"V
, o
°g
t *•
'^v
. MX
>• -I,
§i
oo
y tr.
tH >
>- W
Z «•
O L.
O HI
«r t.
U3.
MM
U. k. l>-
o 1.1 >•• <" i
U "~
LJ II U c.
•_, r, I
5Soz
»• M >• X
z tr —
o a. o L.
o x o M
» K.
(.• u
CO
r-
o
u X
I
_ l-l
-. 5
.° fc
fi X
u r-
Lt •—
h-
* <
i: r t_
"O-Q C.
?.0>- E
0.1- n u
r: M
LJ CJ
1>1 k)
ur-oj(T>u ^t
- oj 01 LJ r-i r,, c> »
110
------- An error occurred while trying to OCR this image.
------- An error occurred while trying to OCR this image.
-------
2 — .-i
- f I
r I. '>
Z cj t-
• *"
^-M o *"-»M^^»-*n
3t L1 M l-i » « ,, « II
oz ^«j: M**«-eriii-»M—
cu. v»:i »<«'r— >-<— j
i — ni'M— -J—
r> n u u u. u M u u u. ^ « _ u — D
t^oj:j;nxTi-'_iM«n <^«z
» *i rj
(j u u
U (J u
cjo0uuU<-'oucjucjtJuiJuuui_iC;ou
OOUCJOQCJUOOUOOUOUCJUUOOO
113
------- An error occurred while trying to OCR this image.
------- An error occurred while trying to OCR this image.
-------
u
o
L
i» £._ CD JT »
OCw V V u
U t- 0: _ X 0
-I 3 u. = « ,.
o o 5 r r c
a a u LJ u cj u t- oi
116
-------
oc
a
x
a
U
0 U
U Cl I
»« ^4 •
L.
(u
^ •
•• »
U
f*
-<
u
0
•4
•H
L.
O
U
U
M
U
c,
^4
U
CI
u a-
p« »
a
n
i)
Cl
*^
U
U
L.
n^rt — '4u'4ajiH<-ii-ia'4i-i
cjor1;»*;"f;«ei'---f''"»'';u>^cuuiu-<«>i(n»i.iu)r.wf~<0
^ay^^SSot'S?,'?'''H2r'^"'o'>i"'-''"MMM
«
117
-------
— — a-
._
U l^ U C- U CJ l_i
r- CO (T> (_) r-t
un m u^ u> LO
118
-------
119
------- An error occurred while trying to OCR this image.
-------
t r-
»c.;
>•• P*
_ »" X
?: s
01'
c- o
C'
y
o
«
X C'r^^uj --^J f, ^j o*»fl>-<
<~i «-»f'(^LjocjQ:uCi.tttjc_i«uJi-<
i. u -J .J i.*- «=.i;^nxl ( »- s _i
i: n _i — «•• cj ii " >• ii > >-- ii >~ c — >-, c: _j c..
CC'*>--l_jOCJCiL'(^l_l.'cil.jl_:L.c:;0-»^
< i s; i> M i- u. «; c- ^ re x c- « s; c- >•• M r: L. u L..
' O LJ U U U LJ C) tJ U U (^ U t. L- (-» t
121
-------
PROGRAM HA IN II HPUT .OU TPliT ( TAPES' INPUT , TAP Eb'OUTPUT .TAP E2 ,T APE 3 t
DIMENSION CHI (6. 13. 13 I.CHI MAX (6, 1 3 . 1 3 ) ,CH 1 SA V ( 6 , 13 , 13 ) .CHILAVU.13
l,13>.XSaURCI10l,YSOURCI10) .KCOMPI 6 ) .KPOLUT 16 , 1 3) .MSTAB V 1 24 )
DIMENSION FACMAP(17.13),UAU7.13),U8U7,13).VA<17,13) ,V8< 17.13)
DIMENSION U(2).vmtOX(2).OY(2),G(2>,XG<2> . YG I 2 ) . WHOLE 16 >
UMENS10N MI <13),KX1SAV( 11) .KX1MAX(13),KXILAVU3) .HHMVI24)
01 HENS ION XRE(IO), YRE ( 10 I . RELAT (10 > .RELQNI 10). CHI 6, 10) .CHM AX I 6 . 10)
DIMENSION CHSAVI6.10) , CHLAVI 6 , 10) , OVI 6 . 10 ) .SORLAT 1 10 ) , SORLONI 10 )
COMMON NINCV(10).00«>,100.!0).CBKG<6).OS< 10.100) .UU( 10.100) .HI 10) .
1SICMYI 10,100) .C&RNDI7. lOO.b). ISTABI10, laO I ,HH 1 10. 100 ) , XP ( 10 , 100 ) , Y
2P
-------
XSUR»(XUR-XLL)»OIV
YSUR«(YUR-YLL)»OIV
NSX«XSUR»1.001
NSY»YSUR+1.001
NX1»NWX1«IOIV
NY1>NUV1*IOIV
FNX-NX1
FNY«NY1
NX«NX1»1
CFAC'3600.*OTIME/DH
RcAD(5.1005)OPTION.W)Oe,FACTOR,SUBMlN.SUeHAX
RtAOIS.lOMHCBKGUPI ,KP«1.NCOMP)
FURHATIdElO.OI
R6AO<5.1001MHMULEIKP),KP»1,NCOMPI
FMOLE«10.0»»12
C3 1G45 KP«1.NCCMP
1C45 CBKG(KP)«CBKG/< 1.0+SINIA) )
XNMC»-R«SIN«B-EIGHTYI
YNHC'-R*COS(B-EICHTY1
XRE(L)«IXNMC-FLCAT(LLXI I»FLOAT( 10IVI
YRE(L)*(YNNC-FLOAT(LLY))«FLOAT(IOIV)
XRES(L)«(XRE(L )-XLL)»OIV
115 YRES(LJ«(YRE(L)-YLL»*OIV
12 WRITEI6.120U
U01 FORMAT (1H1 ."SOURCE TRANSPORT RECEPTOR ANALYSIS MODEL*)
URITE(6.1203)TL1HIT.OTIME.OHAP.OPRINT
1203 FORNATI1HO. "TOTAL TINE FOR RUN IS-F7.0," HOURS. "/1H ."BASIC TIME 1
INCREMENT IS-F7.0." HOURS. HAP INTERVAL IS"F7.0." HOURS. -/1H , "AVE
2RACINC INTERVAL BETWEEN PRINTED MATRIX SERIES IS-F6.1," HCURS")
WRITE(6.120SUSCURC
U05 FOSMATdHO. "THERE ARE"I3." SOURCE LOCATIONS.")
WRITEI6.1207)
1207 FORMATI1H ." LATITUDES. LONGITUDES. X AND Y COORDINATES ON THE
1 AOVECT10N GRID. AND EFFECTIVE STACK HEIGHTS ARE")
CC 121 K'l.NSOURC
121 WRITE 16.1 311 ISORLATIK) .SORLUNIK) . XSOURC(K) .YSOURC(K) ,H(K)
1211 FQRMATdH ,Zf 10.2 , F20 . 2 .F10.2 , F12.0. " CM.")
HRITEI6.1221 )NCOHP
1221 F3RMATI1H .11," CHEMICAL COMPONENTS ARE ANALYZED.")
123
-------
WRITEI6.1222)I WHOLE(KP ).KP»1,NCOMP )
1222 FORMAT!1H ."MOLECULAR WEIGHTS ARE"6F10.1)
WRITE(6,1223MCBKG7 GRID BY SELEC
1TINO A RECTANGULAR PORTION AND FURTHER SUBDIVIDING.")
WRITE(6.1232)LLX.LLY
1232 FORNATI1H ."THE COORDINATES ON THE NHC GRID OF THE LOWER LEFT CORN
1ER OF THE RECTANGLE ARE"215,".">
WRITEI6.1233)
1233 FCjRMATdH ."THESE ARE SPECIFIED WITH RESPECT TO A RIGHT HANDED COO
10ROINATE SYSTEM WITH THE ORIGIN AT THE POLE.")
WR1TE(6.123<>>
1234 FORMAT(1H ,«THE NEGATIVE Y AXIS EXTENDS DOWNWARD ALONG THE HERIOI
IAN BO DEGREES WEST.")
WRITE(6.12<>1)NWX1,NWY1
12
-------
13*1 FORMAT.YS
1UURC(KI.H(K).K«1.10)
HR!TE(2)NRE.RELAT,RELON
CO IS KP'l.NCOMP
CO IS JM.NSY
DC 15 1*1,NSX
CHISAVIKP.I,J)«0.0
IS CHILAV0.0
XM(K,NP)'0.0
V1RTH(K.NP)«0.0
16 PTIME(K.NP>«0.0
CALL MAPFACIPI.LLX,LLY,101 V.FACMAP.NX,NY)
CALL LATLQN(LLX,LLY,IOtV.XLL.YLL.DlV.NSX.NSY,CHi.PI)
19 TIME'0.0
TMAP»0.0
TPRINT'0.0
ITIME«0
KZFLAG'O
00 20 I'l.NSOURC
20 NINCVI11*1
21 CALL KALENDIIHR.IOAY.IMON.11ER.OTIMt.JHR.JOAY.JMON.JIER>
ITIME»ITIME»1
IF(ITIME.&T.KYCLE)ITIME«1
IFUFST.NE.OIGO TO 22
CALL W'ND(UA.VA,NXA,NXB.NYA,NYB,NREAO,KFST.LFLAG.TIME.CALMbO)
NXTEST«NXB-NXA»1
NYTEST»NYB-«YA»1
IFINXTEST.EO.KX.ANO.NYTEST.EO.NYIGO TO 2<>
MRITEI6.2101)
2101 FOR1ATI1HO."DIMENSIONS ON UINu TAPE DISAGREE WITH THOSE CALCULATED
1 FROM INPUT PARAMETERS")
«RITE(6.2102)NXTEST.NX.NYTEST,NY
;102 FORMATdH . "NXTEST .NX .NYTEST ,NY"<.I 51
STOP
22 TMAP-TMAP+OTIME
IF(TMAP.LE.OMAP)GO TO 25
THAP.0.0
00 23 1*1.NX
CO 23 J«1.NY
UA( 1. J)«UB( I. J)
23 VAIi.J)'VB(I.J)
2
-------
1F
XP(K.N1NC)'X50URC(K)
YP(K,NIMC)*Y50URC(K)
DO 35 NP'l.NINC
XGI1 I'XPI K.NP)
Yd 1)*YP( K.NPI
C8'TMAP/OMAP
CA'l.O-CB
DO 33 KT«1.2
1«XG(KT)+1 .
J»YG«KT)*1.
ISTAB(K.NP)*NSTAB( I.J)
HHIK.NPJ'AMAXl (riHM ( I . J ) ,HH ( K , NP ) )
FI *1-1
JA=J+1
IFUA.GT.NXtlA'NX
IF(JA.GT.NY)JA*NY
RX »XG(KT)-FI
RY»YC(KT»-FJ
RXl'l.O-RX
RY1=1.0-»Y
A'RX1»RY1
8 »RX1»RY
C»RX»RY
D«RX»RY1
CST»A«(CA*UA( I . J)»C8*UBI I.J))
CST*CST*B»(CA»UA( I .JA)»CB»UB( I .JA»
CST«CST+C"(CA«UA( IA, JA )+C8*UB( IA.JA) )
CST'CST+0«
-------
CST«CST»C*CCA«VA«1 A.JA)+C8»VB(1A.JA))
CST»CST»0»(CA»VA<1A,J)+CB»VBI JA.JI)
V(KT)«C5T
GUT»«FACMAP»OX(1I
YG(2)«YG(l)»DY(l)
IF(XG(2».LE.O.O.OR.XG<2».GE.FNX)GO TO 32
IF(YG(2).GT.O.O.AND.YG(2).l.T.FNY)GO TO 33
32 XPIK.NPI«X&(2)
YP(K.NPI«Y&(2)
GO TO 35
33 CONTINUE
CXM«0.5«(OX(1I»OX(2) I
CYH»0.5»(OY(11+OYI2))
XP(K,NP|"XP(K.NP1+OXM
YPIK,HP)»YP(K.NP1+OYM
UU(KiNP)'50RT((U(1)+U(2I)»»2«(V(1)+V(2
UU(K.NP)»0.5»UU(K,NPI
DSH»SCRT(OXM»«2»DYM»»2)
OSIK.NP)«2 .0»OSM»OM/(i(!)»&(2>)
IF(DPT I ON.EJ.1 (CALL PLMRIZ»SIGMZ
05U.I )>OS(K.I«NPR)
UUIK.1 )«UUIK,I»NPR )
XMIK.II-XMIK.1+NPR)
PTIME(K,I)»PTIHE(K.I»NPR)
ISTABIK,I)«ISTAB(K.1«NPRI
HHIK,I)»HH(K.I«NPR)
37 VIRTH1K.I)-VIRTHIK.l+NPR)
38 N1NCV(KI=NP-1
GO TO 392
39 CONTINUE
392 CALL STRACIK.SUBMIN.SUBMAX.FACTOR,MODE.TIME.OPTION,NCOMP)
IF(NINCVIK).LT.21GO TO 40
CALL SAMPLE(XP.YP,XSOURC.YSOURC,NINCV,CHI.CGRND.FACMAP.SIGMY.OHS.K
l.NX.NY.NCOHP.XLL.YLL.NSX.NSY,01V,NRE.XRES,YRf$,CH
<>Q CONTINUE
TIHE«TIME+DTIME
TPRINT-TPRINT+OTIME
127
-------
DO 42 KP'l.NCOMP
IFtKCOMPIKP).EO.OIGU TO 42
WRlTE(2)T|HE.KP,WMOLEiKP),IICHllKP.I.J).I«l.NSX).J«l.NSY)
IFNINCV(KI»1
IF(N1NCVIK).GT.IMAX1CQ TO 65
61 CONTINUE
IFITIME.LT.TLIMITIGC TO 21
GO TO 80
65 WRITE(6.6501)N1NCV
6501 FORMAT(IHO.-CNE PLUME HAS TOO HANY INCREMENTS"1014)
NR1TEI6,65021JHR.JOAY,JHON.JIER
6502 FORMATUH ."TIME OF TERMINATION IS"415." I HOUR.DAY.MONTH,YEAR)"I
URITEI6,6503)TIME
6503 FORMAT(1H ."RUNNING TIME IS-F7.1." HOURS.")
WRITE 16.6504)IMAX
6504 FORHATtlH ."MAXIMUM ALLOWED NUMBER OF INCREMENTS PER PIUME IS-I4)
TLIM1T«TIME
GO TO 80
70 HRITEI6.7001)
7C01 FQRHATI1HO.-ENO OF FILE ENfQUNTEREO ON WIND TAPE")
WRITE 16,6502)JHR.JOAY,JMON.JIER
WRITE!6.65031T IME
TLIMIT»TIME
aO CALL PRCHIIKXlMAX.CHIMAX.NSX.NSY.DPRINT.NCaHP.NRE.XRES.YRES.CHHAX.
1JHR.JOAY.JMQN.JIER.KPQLUT.FMOLE,WHOLE)
CALL PRCHI(KXILAV.CHILAV.NSX.NSY.TLINlT.NCOMP.NRE.XIiES.YRES.CHLAV.
1JHR.JDAY.JHQN.JIER.KPOLUT.FHOLE,WHOLE)
ENOflLE 2
WR1TEI6.8001)
128
-------
8001 FORMAT UNO."TAPE WRITTEN WITH EOF")
STOP
END
SUBROUT INE 0 IAGNOINSOURC,NfONP.NX,NY.NSX.NSY.LLX.LLY.XSOURC,YSOURC
l.FNX.FNY.NRE.XRES.YRES.XSUR.VSUR)
DIMENSION XSOURCI10).YSOURft10),XRES(10).YRES(10)
IFINSOURC.GT.O.AND.NSOURC.LE.10IGO TO 11
WRITEI6.1001INSCURC
1001 FORMATI1M ."NSOURC —IS,'. 10 15 UPPER LIMIT."}
STOP
11 IFINCOMP .GT.O.AND.NO1MP.LE .6100 TO 12
WRITEI6.1101INCOMP
1101 FQRMATI1H ."NCOMP «"15.". 6 IS UPPER LIMIT")
STOP
12 1FINX.GE.2.ANO.NX.LE.17IGO TO 13
WRITE! 6.120DNX
1201 FORMATI1H ."NX '"15.". LIMITS ARE 2 AND 17")
STOP
13 IFINY.CE.2.ANO.NY.LE.13ICQ TO !<•
WRITEI6.1301IKY
1301 FORNATI1H ,»NY «"I5.". LIMITS ARE 2 AND 13")
STOP
14 IF(NSX.GE.2.ANO.NSX.LE.13)GQ TO IS
wRlTE(b,l<.01)N5X
1*01 FORMATIIH ."NSX '"IS.". LIMITS ARE 2 AND 13")
STOP
IS IF(NSY.GE.2.ANO.NSY.LE.13)GO TO 16
WRITEI6.1501INSY
1S01 FGRMATI1M ."NSY «"!5.". LIMITS ARE 2 AND 13")
STOP
16 IF (LLX.GE.-2't.ANO.LLX.LE.l IGO TO 17
WR ITEIb.16011LLX
1601 FQRMATIlrt ."LLX «»I6.". LIMITS ARE -2* AND +1")
STOP
17 IF(LLY.GE.-26.AND.LLY.LE.-5)GO TO 175
WRITEI6, 1701ILLY
1701 FORMATI1H ,"LLY. «"16.". LIMITS ARE -26 AND -5")
STOP
175 KXFLAG'O
KYFLAG«0
LXFLAG'O
LYFLAG'O
00 18 1=1,NSOURC
IFUSOURCin.LT.O.O.OR.XSOURCt 1>.GT.FNX)KXFLAG*1
IF(YSOURC(I).LT.O.O.OR.YSOURCID.GT.FNYIKYFLAG'l
IF1KXFLAG.NE.OWRITEI6.1751 (KXFLAG
IF(KYFLAG.NE.O)WRITE(6.1752)KYFLAG
1751 FORMATIIH ."X COORDINATE OF SOURCE"!*.." OFF ADVECTION GRID")
1752 FORMATIIH ."Y COORDINATE OF SOURCE"!*." OFF ADVECTION GRID").
18 CONTINUE
IFINRE.EO.OIGQ TO 20
CO 19 1*1,NRE
I FIX RESt I) .LT.O.O.OR.XRESII>.GT.XSUR)LXFLAG«I
!F(YRES(D.LT.O.O.OR.YRESII>.GT.YSUR>LYFLAG-I
IF(LXFLAC.NE.O) WRITE(6.1851ILXFLAG
IFILYFLAG.NE.OIHRITEI6. 185i)LYFLAG
1?S1 FORMATIIH ."X CCOR01NATE Of RECEPTOR "!<.." OFF SAMPLING GRID")
129
-------
1652 FORMAT (1H ,"Y COORDINATE OF RECEPTOR "!<.." OFF SAMPLING GRID")
19 CONTINUE
20 IFtKXFLAG.NE.O.OR.KYFLAG.NF.OJSTOP
!F(LXFLAG.NE.O.OR.LYFLAG.NE.O)STOP
99 RFTURk
END
SUBROUT INE KALENOI I MR, I DAY . I HON. 1 1 ER.DTIHE . JHR , JDAY, JMON. J IER )
DIMENSION KDAYI12)
DATA KDAY/ 31. 28, 31. 30, 3 1.30, 31. 31, 30. 3 1,30. 31 /
CATA KF5T/0/
IF(KFST.NE.O»GO TO 2
KFST«1
IF( 1.EQ.I1ER)KDAY(2)«29
JHR>I(-R
JDAYMDAY
JMON'IMON
Jl ER'IIER
2 IT!ME«DTIME+0.01
JHR*JHR+ITIME
IF ( JHR.i.T.24 IRETURN
JHR*JHR-24
JDAY«JDAY+1
IF ( JDAY.LE.KOAYIJMQNI (RETURN
JDAY'l
JHDN'JMON+1
IFUMON.LE .1? IRtTURN
JMON»1
JIER'J1ER+1
IF(KOAY(2).E0.28)GO TO 3
KDAY(2)«28
RETURN
3 I"»»(JIER/<>>
IF(I .EQ.JIER)KOAY(2)>29
RETURN
END
SUBROUTINE PRCH1(LTITLE.X,NSX,NSY,DVT1ME.NCOMP,NRE.XRES,YRES.Y.JHR
1, JDAY. JMON, J YE R.KPOLUT.FHOLE. WHOLE )
DIMENSION LTITLE(13).X(6,13,13I , XRES ( 1 0 ) , YRES ( 10 ) , Y ( 6 , 10 ) , KPOLUT 1 6
1.13I.HMULE (6)
30 9 N=1.NCDKP
MRITEI&.1001 ILTITLE
1C01 FORMAT UH1.13A6)
URITE(6,1002IN.(KPOLUT(N.I),I>1,13)
1002 FORMATdH . U , 5X . 13A6 , " Ml fROGRAMS/CUB 1C METER"!
WRITE! 6,1003 ) JHR, JDAY, JMON. JYER
1003 FORMATdH ."TIME AT END OF INTERVAL"*! 5." (HOUR .DAY .MONTH, YEAR )n )
HC'0.0
ZIKFAC=FMOLE«WMOLE(N)/DVTIME
CO 5 J'l.NSY
IFCDVTIME.EO.l.OIGO TO 3
CO 2 1=1. NSX
2 X(N , I , JR)'X(N. i, JR) »Z IKFAC
3 HR:TE(6.3001)(X(N,I .JR) ,1'l.NSX)
3001 FORMAT11H0.13F10.3)
DO <• I'l.NSX
130
-------
4
MJ01
5
5C01
6
6C01
7001
SCO I
85
9
1
1001
2
2C01
HC'AMAXl (MC.XIN.l ,JR) )
X(N. I , JR»«0.0
FORMAT (1HO."HAX|MUM ELEMENT IN MATRIX IS"F12.3)
CONTINUE
HR1TE(6.4001)HC
IFINRE. £0.0)00 TO 9
WRITEtb, 50011
FORMATUH .//1HO. "SIMILAR FIGURES FOR PARTICULAR RECEPTOR SITES")
00 6 L'l.NRE
V(N.L)«Y
CO SS L'l.NRE
YIN.D'0.0
CONTINUE
RETURN
END
SUBROUTINE
CIHENSION
CIV'ICIV
SlN60>SIN(PI/3.)
DC 1 J'l.NY
YC*FLOATILLV)»«FLOAT(J)-1.0)/DIV
DO 1 Ml. NX
xc = FLOAT (LLx)»t FLOAT (i >-I.OI/DIV
R«SORT(XC«XC*YC»YC)
PHI»0.5»PI-2.0»ATAN(R/31.18<.)
FMI1 ,J )«(! .0 + 51 N60)/( 1.3+SINI PHI ) I
»RITE( 6.1001)
FORMAT UHO. "HAP FACTORS")
00 2 J'l.NY
K=NY+1-J
HRITE(6,2001)IFM(I ,K) ,I»1,NX)
FORNATUH .20F6.3I
RETURN
END
SUBROUTINE LATLONILLX.LLY. 101 V .XLL .YLL .01 V ,NSX .NSY .CHI .P I )
DIMENSION CHK6.13.13l
MAPFACIPI.LLX.LLY.IDIV.FM.NX.NY)
FM!17,13)
1101
XL 'FLOAT! LLXl+XLl/ FLOAT! IOIV)
YL»FLOAT(LLY)+YLL/ FLOAT! ID IV)
DEL=1.0/(FLOAT( 10IV)»OIV)
DO 1 J'l.NSY
Y»YL*DEL«(FLOAT(J)-1.0)
CO 1 I'l.NSX
X»XL+OEL«(FLOAT( D-1.0)
R'SORT(X«X*Y»Y)
CHI! 1,I,J)=90.0-CONV»2.0»ATAN»R/31
CHI (2, I. ,)«aO.O+CONV«ATAN2<-X,-Y)
WRITEI6. 1101 )
FORMATI1H1. "LATITUDE AND LONGITUDE
1")
DO 2 J*1.NSY
FOR SAMPLING GRID INTERSECTIONS
131
-------
JR«NSY+1-J
Wft!TE(6.1201)(CHIU.ItJR).IM,NSX)
1201 FORMAT(1H0.13F10.21
WRITE(6.1301MCHI(2.1,JR).I*1,NSX)
1301 FQRMAT11H .13F10.2)
i CONTINUE
DO 3 K*1.2
00 3 J'l.NSY
CO 3 1*1.NSX
3 CHI(K.I,J)-0.0
RETURN
END
SUBROUTINE *INDIP.O.NXA.NX.NYA.NY,NREAO.KFST.L FLAG,TIME.CALMKC)
DIMENSION 10(121.U<13,9),V«13,91,P(17,131,0117,131
10 REAOOHYR.lOm.lDU) ,IDISI.ZD.ZT.NXA.NX.NYA.NY
IFIEQF(3) 120.1 1
11 *RlTE{6,noniYR.lDC3I.IO«.),IDt5),ZO,ZT.NXA.NX,NYA,NY
1101 FQRMATI1H .',112,2612.3.411?)
REAOI3H(U(I,J),I«NXA,NX),J*SYA,NY )
READ(3)1»100»IO«.)*IDI5)
IFIISREAD.GT.NTAPE )GO TO 10
KFST=1
15 CO 16 I'NXA.NX
00 16 J'NYA.NY
PI I.J)»1CO.O»U
IF (SPEED.GE.CALMMOIGO TO 19
IFISPEED.GT.O.OIGO TO 17
TMETA»TIME/H20.0«3.1<.1592h5)
GO TO 18
17 THETA«ATAN2(Q( I,J).Pi I .J) )
18 PM.J) »CAIM«0*COSITHETA)
0( I. J) *CALMMO»SIN(THFTA)
19 CONTINUE
RETURN
20 LFLAG*!
RETURN
END
SUBROUTINE PLHR I Z(UU.ri,VIRTH.K,NP)
DI KENS ION UU I 10,100) .HI 10) .V I RTHl 10.100 )
VlRTHIK.NPl'HIK)
RETURN
END
SUBRCUT INE STRAC IM .SUBMIN.SUBMAX.FACTOR.MODE,TI ME.OPT ION,KKK)
COMMON SNUO). 00(6. IOC. 10) .CBKGI6) .US ( 10,100 ) . UI 10 . IOC ) ,H( 101 ,
1S1GMYI 10,100),CGRND(7,100,61,1 STAB 110,100),HH(10.100),XP(10,100)
2YP(10.100).XH(10.100I,SIGMZ(10.100 ).PTI«EI10,100).XRES(IO).
3YRESI10)-VIRTHI10,100)
CIMENSION VO(6),UOC(6)
DIMENSION AD(6 ) , Ab 16>,GAINI 61,GOT{6).COP I 6 I,YDUMI 10,100)
INTEGER OPTION
1002 FORMATI1H .12F10.0)
132
-------
x*o.
N«NN(M)
00 10 K»l,KKK
GA1N GO TO 145
CALL INTGRTIOPTlON.OOT.DX.THE.SIGY.SIGZtN.N.KKK.GAINI
IF(MOOE.E0.3) G3 TO 210
C CALCULATION OF WASHOUT AND DRY DEPOSITION FLUXES
145 CALL MSHGUTINtUCCtN.KKKtTIHE)
CALL ORYOEPIN,VD.M.KKK.TIME )
D2*VIRTri(M.N)**2/(2.»SIGZ*SIGZ)
CO 367 K'l.KKK
IHOPTION.EQ.2) GO TO 366
AD(K)«0.0
IF(C2.LE.30.0IAO(K)«0.7979»QOTlK>«VOtK)»EXPJ-D2)«OX/ISIGZ«U»HH(M,N))
367 AU(K)'WOC(K)«QUTIK)»OX/UIM.N)
C CALCULATION OF VIRTUAL SOURCE STRENGTHS AT INTERVAL END
210 DO 211 K=1.KKK
QQF
C6»VIRTH(N.N)»»2./(2.»SIGZ«SIGZ»
DC 2<>0 K'l.KKK
C5*UO(K,N.M)/(3.1<>159*$1GY*SIGZ»U1.7
YDUM(ItNI'FLOAT I!-l)»OY
133
-------
CGRNDll.N,K)*Q.O
ARG«D6+0.5»IYOUM( 1 .N)/S 1&Y)*»2
IFUR&.&T.40.01GO TO 240
CGRNDlI.N,K)*D5»EXPI-ARG)
24C CONTINUE
GO TO 445
444 D1*0.3989/(U
525 SIGY*5IGMY(M.,M)
51GZ'5IGMZ(H,N>
IF (MCOE.E0.2) GO TO 527
CALL I NT CRT(OPT ION.OOF.DX.TIME.SI&Y.SIGZ.M.N.KKK.&AIN )
IFIHODE.EQ.3) Gb TO 100
527 CALL WSHOUTIN.UGC.H.KKK.TIME)
CALL DRYDEP1N.VU.M.KKK.T1ME)
CD 867 K-l.KKK
IFIQPTION.EQ.2 ) GO TO 866
l)2»G.5«(VIRTH
-------
2 *.330152746764, 11 .8437*5 837900 . 16.279257831378.
3 21.99658581198, 29.92069701227W
IF(OPTION.E0.2I CO TO 666
D7».1592/(U{M.N»»SIGY»SIGZ»
CO 130 1*1.10
SCALING FOR QUADRATURE IN I
Z«SIGZ»
DO 120 JM.10
SCALING FOR QUADRATURE IN V
Y«SIGY»ZLAGR(JJ/6.
D9»EXPt-V«Y/(2.»SIGY»SIGYI>
CC 110 K'l.KKK
110 CIK)'QextK)«07*D8»D9+CBKGtH>
120 CALL RATE(C,FUNCT.J.M.N.TIHEI
DO 130 K'l.KKK
00 135 11-1,10
135 FUNIIII'FUNCT(II.K)
INTEGRATION OVER Y
CALL LAGURE(ZINT.FUN.ZLAGR.U)
DESCALING OF Y INTEGRAL
133 ZlNtI,K>«ZlHT«Sl&Y/6.
DO 140 K'l.KKK
CO 143 I'l.lO
143 I II D'ZINl I.K)
INTEGRATION OVER Z
CALL LAGURE(ZINT.Zl.ZLAGR.W)
DESCALING OF Z INTEGRAL AND CALCULATION OF GAIN BY REACTION
140 GAINIK)*2.»OX«SlGZ*ZINT/3.
GO rC 999
666 01'.3989/(U(H.N)«HH(M,N)*S1GY)
C2»SIGY»SIGY«2.
SCALING FOR QUADRATURE IN Y
CO 778 IM.10
Y«S1GY«ZLAGR( I 1/6.
09 «£XP(-Y»Y/D2I
30 779 K*1.KKK
779 C(K)*01»09»QQX(K)»C8KG(K)
778 CALL RATE(C.FUNCT.1.H.N.TIME)
CO 780 K»1.KKK
CO 781 I I'l.lO
781 FUNl I I 1'FUNCTlI I ,K)
INTEGRATION IN Y
CALL LAGURE(ZINT.FUN.ZLAGR.U)
CESCALING OF Y INTEGRAL
ZIN(l.K>>ZINT*SIGY/6
780 &AIN(Kl»HH(M,N)*ZIN{1,KI«OX«2.
999 RETURN
END
SUBROUTINE SIGHAI X .S IGY .SK.Z .H .N)
COMMON NINCVdOl.&0l&,100.101,CBK&<6),DSt10.100) .U( 10,100) ,H(10) ,
135
-------
1SIGHYI 10,100) .CoRNDI 7, 100,6) .1 STAB! 10,100 ) ,HH( 10.1 001, XPI1 0,100).
2YPU0.100>.XNU0.100),SIGHni0.100).PTINEUO,100>.XRESUC) .YRESI 10
3) .VIRTHl 10.100)
CIHENS10N YA(6),ZA«6) ,ZB<6 )
3 ATA YA/0.33.0.22,0.1 7,0.12.0 .086, 0.05 7 /
DATA Z A/ 0. 00048, 0.06 3, 0.100,0. 33. 0 .<• 1,0. 32 1
CATA ZB/1.10.0.G9.-0.09.-0.42.-0.53.-0.58/
OX *X-XM(M,N)
CT'OX/UIM.NI
PTIHE(H,N)»PT1HE
OSYOX«YAtKSTAB)»XC»»(-0.1)
CSZDX«ZA( fSTAB)»XC»»(ZaiKSTABM
5IGY»S1&MY(H.M)+DSYDX»CX
S1GZ«5IGMZ(M.N)+OSZDX»OX
GO TO 5
<. TC *PTIME(M.N)-0.5»OT
SIGZ«51GKZI».N)+158.»OT/SCISTITC>
5 SIGHY(N.N)0.
WQC(2)'0.
RETURN
END
SUBROUTINE DRYDEP ( N ,VD ,M,Nf OHP , T I IE )
CIHENSION NINCV(IO) ,VO (6)
VD( ll«l.
VDIil- .1
RETURN
END
SUBROUT INE SAHPLE ( XP .YP .XSOURC .YSOUAC , NINC V.CHI .CGRNO , FACHAP, .MGHY
1 .DHS.K.NX.NY.NCUMP.XLL .YLL .NSX.NSY ,DI V .NRE ,XRES , YRES .CH)
DIMENSION XP(IO.IOO) .YPU0.100) .XSDURC(IO) .rSQURCI 10>,NINCV(10)
CIHENSION CHI (6,13.13) .CGRNDI 7. 100.6 I .FACH API 17.13 ) .S1GHY I 10.100).
IXRES(10).YRES(10),CHI6,10)
136
-------
N»MNCV(K)
EPSIL»0.1««6
CO <* HP«Z,N
M-N+2-NP
IFIM.LT.NJGO TO 1
PX«(XSOURC(Kl-XLL)*OIV
PY*IYSQURC(K)-YLL)»OlV
ax«(XP(K.M)-XLL1«01V
IF(AB5IQX-PXI.LT.EPSIL)QX»PX»EPS1L
IF(ABSUY-PYI.LT.EPSILIQY*PY»EPSIL
KX«XPCK,H»»1.5
LY »YP(K.M) + 1.5
1F(KX.LT.1)KX«1
IF(KX.GT.NSXIKX«NSX
IFILY.LT .1)LY«1
IF(LY.GT.NSY»LY«NSY
F*FACHAP3.0*SIGMY(K,H)
CALL LE IN(PX.PY.OX.OY,S30.0HS,AX.AY.BX .BY .S.F)
PX>OX
PY.QY
S3P'S30
1 3X 'UP(K.M-n-XLL) >01V
QY«IYP(K,M-1)-YLL)»OIV
A8P3X«ABS(OX-PX)
ABPQY«ABS(OY-PY)
IF (ABPOX.LT.EPS ID CX»PX»tPSU
IF(ABPOY.LT.EPSILIOY«PY»EPSIL
KX >XP(K.N-1)*1 .5
LY*YP(K.N-1l+l.5
IFUX.LT.l IKX'l
IF(KX.GT.NSXIKX'NSX
IFILY.LT.1ILY»1
IF(LY.GT.NSYILY«NSY
S3Q«3.0»SIG«Y»K,H-1)
CALL LElNtPX.PY,OX.QV,S3Q.DHS.CX.CY.OX.OY.S.F)
XMAX«AMAX1(AX.BX.CX ,DX)
XM1N.&M1NKAX.BX.CX.OX)
YHAX'AHAXH AY.BY.CY.OY)
YHIN'ANINUAY.BY.CY.OY)
IMN«XNIN»1.999
1MAX«XNAX*1.001
JMAX«YMAX+ 1.001
IF(IMlN.LT.l)IH1N«1
1F(1MAX.GT.NSX)IMAX>N5X
IF(JMIN.LT .1)JM1N«1
IFIJMAX.&T.NSY)JMAX«NSY
IF! 1HAX.LT.1M1NICO TO 31
IF(JMAX.LT.JMINI&Q TO 31
CO 3 J'JMIN.JMAX
Y'J-1
DO 3 1'IHIN.IHAX
X-I-1
0 =SQRT((PX-QXI«»2*IP
YQI$T»«(OY-PV1»X»(PX-OX)*Y*QX»PY-PX*OY)/D
137
-------
YDIST'ABSIYDIST)
SYGRID«S30»F/DMS
IFIYUIST.GT.SYGRIDIGO TU 3
XL «(X*S*(S»PX+Y-PY)»/( 1.0»S*SI
TEST»(CX-XL)*(PX-XU
IF(TEST.GT.O.O)GC TO 3
hC»(XL-PX)/(OX-PX)
WP*1.0-«IQ
SIGY3«53P»HP+S30»«P
S1GY3«SIGY3*F/DMS
FN1«1.0»6.0»YCIST/SIGY3
Nl'FNl
W1*1.0-M2
1F( FN1.&T.7.0)&C TO 3
IF(N2.GT.7IN2»7
CC 2 KP'l.NCOHP
A'C&RNDINl .H.KP)*M1*HP
B'CGRNO(N2tM.KP)«M2*MP
C'C(,RND(Nl.M-l.KP)*Wl»We
0H.KP>aM2«MP
C=CGRND(N1.M-1 ,KP)«W1«WO
0-CGRNOIN2.M-1.KP) »H2»U6
3<. CH (KP.L ) =CHIKP ,L l+A+8 + C+O
35 CONTINUE
138
-------
S3P»S3U
AX'CX
AY»CY
sX *OX
BY*DY
CONTINUE
RETURN
END
SUBROUTINE LEIMPX.PY,CX.OY,S3C,OMS.CX,CY.OX,OY.S,F)
S« OY-PY l/(gX-PX)
SI >-l .0/S
CP'SQRTI 1 .0/1 1 .0«S 1»S I ) I
5P=SORT(1.C-CP»CPI
IFISl.LT .O.OIS?'-SP
OEX»Oh*CP
OEY«DH*SP
CX«;X-OEX
CY«3Y-DEY
ex »JX + CEX
CY«JY+DEY
Rt TURN
END
139
-------
SECTION VI
REFERENCES
1. Carnaham, B., Luther, H. A., and Wilkes, J. O., Applied Numerical
Methods, Wiley, New York, (1969).
2. Jenne, R. L., "The NMC Octagonal Grid," National Center for Atmo-
spheric Research, Boulder, CO, (1970) unpublished.
3. Wendell, L. L., "Mesoscale Wind Fields and Transport Estimates
Determined for a Network of Wind Towers," Monthly Weather Review,
No. 100, pp. 565-578 (1960).
4. Levy, A., Drewes, D. R., and Hales, J. M., "SO Oxidation in Plumes:
A Review and Assessment of Relevant Mechanistic and Rate Studies,"
Final Report to EPA, Contract 68-02-1982 (1976) in press.
5. Dana, M. T. and Hales, J. M., "Statistical Aspects of the Washout
of Polydisperse Aerosol," Atm. Environ., No. 10, pp. 45-50 (1976)
6. Engelmann, R. G. and Sehmel, G. A., "Atmosphere-Surface Exchange of
Particulate and Gaseous Pollutants," ERDA Document CONF-740921, U.S.
Energy Research and Development Administration, Oak Ridge, TN (1976).
7. Turner, D. B., Workbook of Atmospheric Dispersion Estimates, Public
Health Service Publication No. 99-AP-26, U.S. Department of Health,
Education and Welfare, Cincinnati, OH (1969).
8. Heffter, J. L. and Ferber, G. J., A Regional-Continental Scale Trans-
port, Diffusion, and Deposition Model - Part II; Diffusion-Deposition
Models, NOAA Technical Memorandum ERL ARL-50, Air Resources Laboratories,
Silver Springs, MD, pp. 17-28 (1975).
9. Norton, C. and Hoidale, G., "The Diurnal Variation of Mixing Height By
Month Over White Sands Missile Range, NM," Research and Development
Technical Report ECOM-5579, Atmospheric Sciences Laboratory, U.S. Army
Electronics Command, White Sands Missile Range, NM 88002, November 1975.
140
-------
SECTION VII
STRAM GLOSSARY
AD - Vector of loss rates of material from a plume segment due to dry
deposition, mole/sec
AW - Vector of loss rates of material from a plume segment due to wet
deposition, mole/sec
CA - 1.0 - CB.
CB - Fraction of map interval DMAP between time when last map was read
in and end of current advection interval DTIME.
CGRND - Ground level concentration matrix computed in subroutine STRAC
in coordinate system along the plume increment axis (crosswind position,
increment, chemical component). moles/cm3*
CH - Matrix of concentrations for particularly specified sampling sites
averaged over the basic advection interval DTIME. printed as ygm/nr
CHI - Matrix of gridded concentrations averaged over the basic advection
interval DTIME, typically one hour, printed as ygm/m , written to tape
as gram mole/cia3
CHILAV - Matrix of gridded concentrations averaged over the total running
time TLIMIT. printed as ygm/m3
CHIMAX - Matrix of maximum values over all matrices CHISAV. printed as ygm/m3
CHISAV - Matrix of gridded concentrations averaged over the intermediate
time interval DPRINT, which is typically 24 hours, printed as ygm/m3
CHLAV - Matrix of particular point concentrations averaged over the total
running time TLIMIT. printed as ygm/ra3
o
CHMAX - Matrix of maximum values over all vectors CHSAV. printed as ygm/mj
CHSAV - Matrix of particular point concentrations averaged over the inter-
mediate time interval DPRINT. printed as ygm/m3
DPS - Subinterval of major downwind grid increment (DS) used for Runge-
Kutta integration procedure, cm
DDSOLD - Previous value of DOS. cm
DIV - Linear subdivision factor between the sampling grid and the advec-
tion grid.
*Read in as ygm/m
141
-------
DM - Interval between consecutive grid intersection on the advection grid
at latitude 60°N. cm
DMAP - Time between successive gridded wind maps, hr
QMS - Interval between consecutive grid intersection on the sampling grid
at latitude 60°N. cm
DPRINT - Time of the intermediate interval over which concentrations are
averaged—the standard interval between concentration matrix printouts. hr
DS_ - Interval traversed by a plume increment during a basic advection step—
matrix (source, increment), cm
DTIME - Time of the basic interval used for trajectory calculation—also
the shortest interval over which concentrations are averaged, hr
DSYDX, DSZDY - Derivatives with travel distance of SIGMY and SIGMZ,
respectively.
DX - Downwind distance increment for Runge-Kutta integration, cm
FACTOR - Ratio by which the interval of the integration along the plume
axis in Subroutine STRAC is increased from one step to the next within
the same basic advection step.
fi «-6 "} "}
FMOLE - (10 vigm/gm)/(10 m /cm ) - conversion factor.
FUN - Rate of chemical reaction (dummy variable). moles/cm sec
FUNCT - Rate of chemical reaction, moles/cm sec
GAIN - Gain of material in plume segment by chemical reaction, moles/sec
H_- Height vector for the plume sources. cm
HH - Nondecreasing mixing height applying to a particular plume segment,
drawn from the gridded mixing height matrix HHM according to plume location.
HHM - Matrix of gridded mixing heights, cm
HHMV - Vector of mixing heights sufficient for one 24-hour cycle—one for
each interval DTIME. cm
IDAY - Initial day of the run, in time zone of the wind tape.
IDIV - Linear subdivision factor between the advection grid and the 47 x
47 NMC grid.
IHR - Initial hour of the run in time zone of the wind tape.
HER - Initial year of the run minus 1900.
IMAX - maximum number of increments allowed for each plume.
142
-------
IMON - Initial month of the run.
JDAY - Current day within the run in the time zone of the wind tape.
JHR - Current hour within the run in the time zone of the wind tape.
JIER - Current year minus 1900 within the run.
JMON- Current month within the run.
JSOURC - Option allowing one set of source strengths to serve as constants
for an entire run. 0 - read in one set of source strengths; 1 - read in
new source strengths at every interval DTIME.
JSTAB - Option allowing one 24-hour cycle of stabilities and mixing heights
to be repeated diurnally throughout the entire run. 0 - read in one 24-hour
cycle of stabilities and mixing heights, one set of two numbers at each
interval DTIME; 1 - read in new values at interval DTIME throughout the
entire run.
JYER - Current year within the run minus 1900.
KCOMP - A vector for six integers either 0 or 1, determining which of the
six chemicals will be represented by average concentrations over the
interval DTIME on the output tape.
KP - Chemical component index with a value of 1, 2, 3, 4, 5, or 6.
KPOLUT - Alphameric matrix of names for the NCOMP chemical components
under analysis. Each pollutant name is a vector.
KPRINT - Option governing printout of concentration values averaged over
interval DTIME. 0 - omit; 1 - print.
KSTAB - The stability applying to a particular plume segment, drawn from
the gridded stability matrix MSTAB according to the location of the plume.
KXI - Alphameric title vector for the matrix CHI.
KXILAV - Alphameric title vector for the matrix CHILAV.
KXIMAX - Alphameric title vector for the matrix CHIMAX.
KXISAV - Alphameric title vector for the matrix CHISAV.
LLX, LLY - X and Y coordinates on the NMC 47 x 47 grid of the lower left
corner of the advection grid, counted from an origin at the pole. X is
positive to the right, and Y is positive toward the top of the grid.
143
-------
MODE - Dictates calculation sequence in Subroutine STRAC:
MODE » 1 - neither kinetics or deposition are calculated
MODE » 2 - deposition calculated
MODE » 3 - kinetics calculated
MODE * 4 - kinetics and deposition calculated
MRX, MRY - X and Y coordinates on the NMC 47 x 47 grid of the upper right
corner of the advection grid, counted from an origin at the pole. X is
positive to the right, and Y is positive toward the top of the grid.
MSTAB - Matrix of gridded stability (1, 2, 3, 4, 5, 6 according to Pasquill
categories).
MSTABV - Vector of stabilities (integers 1-6 for Pasquill A-F) sufficient
for one 24-hour cycle - one number for each interval DTIME.
NCOMP - Number of chemical components involved in the chemical analysis.
NINCV - The vector indicating the number of increments of each plume
currently on the advection grid.
NP_ - Index for plume increments.
NRE - Option indicating whether or not sampling will be done at particular
receptor sites on the grid, in addition to grid intersections.
NSOURC - Number of plume sources.
NSX, NSY - Numbers of intersection in the horizontal and vertical direc-
tions, respectively, on the sampling grid.
NX, NY - Numbers of intersections in the X and Y directions, respectively,
on the advection grid.
OPTION - An inter-option indicating whether vertical distribution will be
considered (1) Gaussian with no lid limitations or (2) even between the
surface and a lid.
PTIME - Total travel time for plume increments—matrix (source, increment).
hr
Q2 - Source strength matrix (chemical component, source, increment).
moles/sec
QQF - Dummy array for source strengths, moles/cm
QQX - Dummy array for source strengths, moles/cm
QV - Matrix of source strengths (component, source) of emission at one
interval DTIME, or by option for the entire run. gm/sec
KELAT, RELON - Vectors of latitude and longitude, respectively, for the
particularly specified sampling locations.
144
-------
SIGMY, SIGMZ - Lateral and vertical standard deviation (a , a ), respec-
tively, at a point along the plume axis—matrices (source, increment), cm
SORIAT, SORLON - Latitude and longitude vectors, respectively, for the
plume sources.
SUBMAX - Maximum interval of integration along plume axis in Subroutine
STRAC. cm
SUBMIN - Minimum interval of integration along plume axis in Subroutine
STRAC. cm
TIME - Time between the initial time of the run and the initial time of
the current basic advection interval, hr
TLIMIT - Total time of a run. hr
TMAP - Time counter incremented between values of zero and DMAP at which
value new wind maps are read in. hr
TOTAL - Distance from last emitted end of plume segment centerline to
end of current integration interval, cm
TPRINT - Time counter incremented between values of zero and DPRINT, at
which value the average concentration matrices CHISAV and CHSAV are printed
out. hr
UU - Average wind speed for a plume increment during a basic advection
step—matrix (source, increment). cm/sec
VD - Deposition velocity—vector (chemical component). cm/sec
W - Vector of weighting factors in 'Gauss-Laguerre integration.
WMOLE - Vector of molecular weights for the six chemical components, gm/mole
WQC - Washout rate—vector (chemical component). sec
X - (in Subroutine PRCHI) Printed matrix of gridded concentrations, ugm/m
XC_ - Total travel distance for midpoint of a given plume increment, m
XLL, YLL - X and Y coordinates on the advection grid of the lower left
corner of the sampling grid. The origin of the advection grid is at the
lower left corner of the advection grid.
XM - Total travel distance for plume increments—matrix (source, incre-
ment) . cm
XP, YP - X and Y coordinates on the advection grid of the plume incre-
ments—matrices (source, increment).
145
-------
XRES, YRES - Vectors of X and Y coordinates on the advection grid of par-
ticular receptor locations.
XSOURC, YSOURC - Vectors of X and Y coordinates on the advection grid for
the plume sources.
XUR, YUR - X and Y coordinates on the advection grid of the upper right
corner of the sampling grid. The origin of the advection grid is at the
lower left corner of the advection grid.
Y_ - (in Subroutine PRCHI) Printed concentrations of particular receptor
locations. ygm/m
ZIKFAC - Conversion factor from mole-hr/cm to ygm/m .
ZINT - Integral returned by Gauss-Laguerre integration subroutine (dummy
variable).
ZIAGR - Transformed independent variable for Gauss-Laguerre quadrature.
146
-------
SECTION VIII
ADDENDA
A. CALM WINDS AND RELATION TO PARAMETER SUBMIN (page 40)
Each set of wind component matrices is edited so that wind speeds are never
less than an arbitrary value defined in the code by CALMWD. Presently its
value is set at 100 cm/sec. Such a provision is necessary because Sub-
routine STRAC executes an infinite loop if the advection over interval
DTIME is less than the value of SUBMIN. If a wind speed is found to be
less than CALMWD but greater than zero, the prevailing wind direction is
used. Otherwise a direction is assigned.
B. CONSTRAINT ON DRY DEPOSITION CALCULATION (page 42)
The present algorithm for computation of dry deposition is not suited to
computations for ground level releases if Gaussian vertical distribution
is assumed (OPTION » 1) under very stable conditions.
\
147
-------
LEGAL NOTICE
This report was prepared by Battelle as an account of sponsored
research activities. Neither Sponsor nor Battelle nor any person acting
on behalf of either:
Makes any warranty or representation, express or implied, with
respect to the accuracy, completeness, or usefulness of the information
contained in this report, or that the use of any information, apparatus,
process, or composition disclosed in this report may not infringe
privately owned rights; or
Assumes any liabilities with respect to the use of, or for damages
resulting from the use of, any information, apparatus, process, or
composition disclosed in this report.
------- |