United States
Environmental Protection
Agency
Environmental Sciences Research
Laboratory
Research Triangle Park NC 2771 1
EPA-600/8-84-027
October 1984
Research and Development
x>EPA
INPUFF—A Single
Source Gaussian
Puff Dispersion
Algorithm
User's Guide
-------
EPA-600/8-84-027
October 1984
INPUFF - A SINGLE SOURCE GAUSSIAN
PUFF DISPERSION ALGORITHM
User ' s Guide
by
William B. Petersen
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, NC 27711
and
Joseph A. Catalano, Thomas Chico, and Tsanying S. Yuen
Aerocomp, Inc.
3303 Harbor Boulevard
Costa Mesa, CA 92626
Contract No. EPA 68-02-3750
Project Officer
D. Bruce Turner
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, NC 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NC
-------
DISCLAIMER
This document has been reviewed in accordance with U.S. Environmental
Agency policy and approved for publication. Mention of trade names or
commercial products does not constitute endorsement or recommendation for
use.
AFFILIATION
Mr. William B. Petersen is a meteorologist in the Meteorology
and Assessment Division, Environmental Protection Agency,
Research Triangle Park, NC. He is on assignment from the
National Oceanic and Atmospheric Administration, U.S. Department
of Commerce. Mr. Joseph A. Catalano, Mr. Thomas Chico, and Mr.
Tsanying S. Yuen are employed by Aerocomp, Inc., Costa Mesa, CA.
11
-------
PREFACE
One area of research within the Meteorology and Assessment
Division is development, evaluation, validation, and application
of models for air quality simulation, photochemistry, and
meteorology. The models must be able to describe air quality and
atmospheric processes affecting the dispersion of airborne
pollutants on scales ranging from local to global. Within the
Division, the Environmental Operations Branch adapts and
evaluates new and existing meteorological dispersion models and
statistical technique models, tailors effective models for
recurring user application, and makes these models available
through EPA's UNAMAP system.
INPUFF is an integrated puff model with a wide range of
applications and flexibility. It is designed to model
semi-instantaneous or continuous point sources over a spatially
and temporally variable wind field. A software plotting package
is also provided to display concentration versus time plots for
each receptor and the puff trajectories after each simulation
time.
Although attempts are made to thoroughly check computer
programs with a wide variety of input data, errors are
occasionally found. Revisions may be obtained as they are issued
by completing and returning the form on the last page of this
guide.
The first four sections of this document are directed to
managers and project directors who wish to evaluate the
applicability of the model to their needs. Sections 5, 6, 9, an!
11 are directed to engineers, meteorologists, and other
scientists who are required to become familiar with the details
of the model. Finally, Sections 7 through 11 are directed to
persons responsible for implementing and executing the program.
Comments and suggestions regarding this publication should be
directed to:
Chief, Environmental Operations Branch
Meteorology and Assessment Division (MD-80)
Environmental Protection Agency
Research Triangle Park, NC 27711.
iii 6-84
-------
Technical questions regarding use of the model may be asked
by calling (919) 541-4564. Users within the Federal Government
may call FTS 629-4564. Copies of the user's guide are available
from the National Technical Information Service (NTIS),
Springfield, VA 22161.
The magnetic tape containing FORTRAN source code for INPUFF
will be contained (along with other dispersion models) in future
versions of the UNAMAP library, which may be ordered from
Computer Products, NTIS, Springfield, VA 22161 (phone number:
(703) 487-4763).
This user's guide is intended to be a living document that is
updated as changes are required. Each page of the User's Guide
to INPUFF has a month and year typed in the lower right hand
corner. Future revisions to this document will be indicated on
this page, and every page that is changed due to the revision
will have a new date printed in the lower right hand corner. The
current version number of INPUFF is 84107. The current version
number of INPUFF and the date associated with it will be given in
the preface of the user's guide. The version number is also
maintained in the source code allowing the user to confirm that
his user's guide and source code are current.
iv 6-84
-------
ABSTRACT
INPUFF is a Gaussian Integrated PUFF model. The Gaussian
puff diffusion equation is used to compute the contribution to
the concentration at each receptor from each puff every time
step. Computations in INPUFF can be made for a single point
source at up to 25 receptor locations. In practice, however, the
number of receptors should be kept to a minimum. In the default
mode, the model assumes a homogeneous wind field. However, the
user has the option of specifying the wind field for each
meteorological period at up to 100 user-defined grid locations.
Three dispersion algorithms are utilized within INPUFF for
dispersion downwind of the source. These include Pasquill's
scheme as discussed by Turner (1970) and a dispersion algorithm
discussed by Irwin (1983), which is a synthesis of Draxler's
(1976) and Cramer's (1976) ideas. The third dispersion scheme is
used for long travel times in which the growth of the puff
becomes proportional to the square root of travel time. A
software plotting package is provided to display concentration
versus time for a given receptor and the puff trajectories after
each simulation time.
6-84
-------
CONTENTS
Preface iii
Abstract v
Figures ix
Tables *
Symbols and Abbreviations xi
Acknowledgments xiii
Executive Summary 1
1. Introduction 2
2. Data-Requirements Checklist 4
3. Features and Limitations 6
4. Basis for INPUFF 8
Gaussian puff methodology 8
P1 ume rise 8
Dispersion algorithms 10
5. Technical Description 11
Gaussian puff equations 11
P1 ume rise 13
Dispersion algorithms 14
Mixing height 18
Atmospheric stability 20
Gridding schemes 20
6. Example Problems 22
Example 1 -- Moving source 22
Example 2 -- Low level source with low
wind speed conditions 24
Example 3 -- Variable wind field 26
7. Computer Aspects of the Model 33
INPUFF 33
Program modules 35
Plot postprocessor 36
8. Input Data Preparation 38
Record input sequence for INPUFF 38
Input data for plot postprocessor 45
9. Sensitivity Analysis 48
Puff combination -- ALPHA 48
Frequency of puff combination -- ISTTIM ... 50
Puff size — STANUM 50
Size of modeling region 52
vii 6-84
-------
CONTENTS (continued)
10. Execution of the Model and Sample Test
11. Interpretation of Output
Example 1 -- Moving source . . . .
Example 2 -- Low level source with
wind speed conditions
low
References
Append ices
A. Plume Rise
B. Listing of
FORTRAN Source Code
54
61
61
69
79
81
85
VI 1 1
6-84
-------
FIGURES
Number Page
1 Gaussian puff model 9
2 Effect of variable mixing height on puff dispersion . 19
3 A possible arrangement of modeling and
receptor grids 21
4 Source path for example 1 23
5 Source-receptor geometry for example 2 25
6 Emission rate versus time plot for example 2 .... 2,5
7 Topographic relief of the modeling region
in example 3 27
8 Wind fields and puff trajectories of each hour
for example 3 29
9 Pollutant concentration field of example 3 32
10 Structure of INPUFF 34
11 Sens t vity of CPU time to ALPHA 49
12 Sens t vity of CPU time to ISTTIM 51
13 Sens t vity of CPU time to STANUM 51
14 Sens t vity of concentrations to STANUM 53
15 Sens t vity of CPU time to size of modeling region . 53
16 Sample job stream for INPUFF 54
17 Output for the sample test 57
18 Annotated output of example 1 64
19 Concentration versus time plots for example 1 .... 68
20 Annotated output of example 2 71
21 Puff locations at the end of each simulation time
for example 2 77
ix 6-84
-------
TABLES
Number Page
1 A Comparison of INPUFF to Other Commonly Used
2
3
4
5
6
7
8
9
10
11
1?,
Definition of Variables Used in Plume Rise Equations
Comparison of INPUFF and PTPLU
Record Input Sequence for INPUFF
Record Input Sequence for Plot Postprocessor ....
Percent Change in Concentrations Using Different
ALPHA Values
Input Data for Example 1
14
15
24
26
28
36
38
46
SO
62
70
6-84
-------
SYMBOLS AND ABBREVIATIONS
Dimensions are abbreviated as follows:
m = mass, 1 = length, t = time, K = temperature
C -- pollutant concentration (m/13)
d -- stack inside diameter (1)
F -- buoyancy flux parameter (l*/t3)
fy -- nondimensional function of travel time for
horizontal dispersion
fz -- nondimensional function of travel time for vertical
dispersion
g -- acceleration due to gravity (1/t )
H -- effective height of plume (1)
h -- stack height above ground (1)
h' -- stack height adjusted for stack downwash (1)
L -- mixing layer depth (1)
Q -- emission rate (m/t)
r -- radial distance from center of puff (1)
s -- stability parameter (t~2)
t -- travel t ime (t)
T — ambient air temperature (K)
Ts -- stack gas temperature (K)
u -- wind speed at stack top (1/t)
vs -- stack gas exit velocity (1/t)
x -- downwind distance (1)
xf -- distance to final rise (1)
x* — distance at which atmospheric turbulence begins
to dominate entrainment (l)
y -- crosswind distance (1)
z -- height above ground (1)
Ah -- plume rise (1)
AT -- temperature difference between ambient air and
stack gas (K)
(AT)C -- temperature difference for crossover from momentum
to buoyancy-dominated plume (K)
36/3z -- vertical potential temperature gradient of a layer
of air (K/l)
TT -- pi , 3.14159
oa -- standard deviation of the horizontal wind angle
(radians)
ae -- standard deviation of the vertical wind angle
(radians)
ar -- horizontal dispersion parameter (1)
xi 6-84
-------
SYMBOLS AND ABBREVIATIONS (continued)
aro
av
aw
°x
°y
°z
azo
initial horizontal
standard deviation
component of the
standard deviation
dispersion (1)
of the horizontal crosswind
wind (1/t)
of the vertical component
of the wind (1/t)
dispersion parameter in the downwind direction (1)
lateral dispersion parameter (1)
vertical dispersion parameter (1)
initial vertical dispersion (1)
XI 1
6-84
-------
ACKNOWLEDGMENTS
The authors wish to express their appreciation to Mr. D.
Bruce Turner and Mr. John S. Irwin for helpful conrments regarding
aspects of the work presented here. Portions of this text were
excerpted from the CRSTER, MPTER, and PTPLU user's guides.
In addition to the listed authors of the user's guide, other
Aerocomp personnel made important contributions to the document.
Source code review and comment optimization were made by Mr.
Frank V. Hale III, who also designed the sensitivity analysis.
Ms. Sarah Cunningham prepared all technical illustrations and
produced the final document. Editorial review was provided by
Ms. Vicki Catalano.
xi i i 6-84
-------
EXECUTIVE SUMMARY
The INPUFF UNtegrated PUFF) computer code is designed to
simulate dispersion from semi - instantaneous or continuous point
sources over a spatially and temporally variable wind field. The
algorithm is based upon Gaussian puff assumptions including a
vertically uniform wind direction field and no pollutant removal
or chemical reactions. INPUFF can estimate concentrations from a
single point source at up to 25 receptors.
INPUFF utilizes three distinct dispersion algorithms. For
short travel time dispersion, the user has the option of using
either the Pasqui11-Gifford (P-G) scheme (Turner, 1970) or the
on-site scheme (Irwin, 1983). The third dispersion algorithm was
designed for use in conjunction with the P-G or on-site schemes.
It is used for long travel times where the growth of the puff is
assumed proportional to the square root of travel time.
Features of the INPUFF computer code include:
Optional stack downwash;
Wind speed extrapolated to release height;
Temporally variable source characteristics;
• Temporally and spatially variable wind field
(user-supplied);
Consideration of terrain effects through
user-supplied wind field; and
• Consideration of moving source.
In addition, a software plotting package has been provided to
display concentration versus time for a given receptor and the
puff trajectories after each simulation time.
In the interest of stimulating research, a simple sensitivity
analysis of several user options is provided in Section 9. Tips
on minimizing computer costs without sacrificing accuracy are
also suggested.
6-84
-------
SECTION 1
INTRODUCTION
INPUFF is a Gaussian plume, integrated puff model with a wide
range of applications. The implied modeling scale is from tens
of meters to tens of kilometers. The model is capable of
addressing the accidental release of a substance over several
minutes, or of modeling the more typical continuous plume from a
stack. Several requests to the Environmental Operations Branch
for assistance in modeling the air quality downwind of
incineration ships prompted the development of an integrated puff
model. INPUFF is, therefore, capable of simulating moving point
sources as well as stationary sources.
Computations in INPUFF can be made for a single point source
at up to 25 receptor locations. In practice, however, the number
of receptor locations should be kept to a minimum. INPUFF is
primarily designed to model a single event during which one
meteorological transition period may occur, such as, going from
afternoon to evening conditions. Up to 144 separate
meteorological periods of the same length may be used to
characterize the meteorology during the event; this provides a
time resolution that ranges from minutes to an hour. The user
has the option of specifying the wind field for each
meteorological period at up to 100 grid locations or allowing the
model to default to a homogeneous wind field. Three dispersion
algorithms are used within INPUFF for dispersion downwind of the
source. The user may use the Pasqui11-Gifford (P-G) scheme
(Turner, 1970) or the on-site scheme (Irwin, 1983) for short
travel time dispersion. The on-site scheme, so named because it
requires specification of the variances of the vertical and
lateral wind direction, is a synthesis of work performed by
Draxler (1976) and Cramer (1976). The third dispersion scheme is
for long travel times in which the growth of the puff becomes
proportional to the square root of time. A software plotting
package has also been provided to display concentrations versus
time for a given receptor and the puff trajectories after each
simulation period.
This document is divided into three parts, each directed to a
different audience: managers, dispersion meteorologists, and
computer specialists. The first four sections are aimed at
managers and project directors who wish to evaluate the
applicability of the model to their needs. Sections 5, 6, 9, and
2 6-84
-------
11 are directed toward dispersion meteorologists or engineers who
are required to become familiar with the details of the model.
Finally, Sections 7 through 11 are directed toward persons
responsible for implementing and executing the program. A
listing of the FORTRAN source statements and a detailed
description of the plume rise algorithm are included in the two
appendices.
6-84
-------
SECTION 2
DATA-REQUIREMENTS CHECKLIST
INPUFF requires data on user options, grid dimensions,
sources, meteorology, receptors, and plotter control. The user
must indicate whether the following options are to be employed:
• Stack-tip downwash,
Source update,
• User-supplied wind field,
Intermediate concentration output, and
Puff information output.
The dimension of the modeling grid must be specified. If the
user-supplied wind field option is implemented, then the
dimension of the meteorological grid along with the size of each
grid rectangle must also be indicated. It is recommended that
both grids be given a common origin. If a puff travels outside
the modeling region, it is deleted from further consideration.
If it travels outside the meteorological grid, but is still
within the modeling region, the last wind experienced by the puff
is used to advect it further.
Information required on the source includes the following:
Locat ion (km),
• Emission rate (g/sec),
Physical stack height (m),
• Stack gas temperature (K),
Stack diameter (m),
• Stack gas velocity (m/sec),
Stack gas volume flow (m3/sec), and
• Initial ar and az (m).
4 6-84
-------
Also the direction and speed of the source, if it is moving, must
be provided as input.
The meteorological data needed for the computations are as
follows:
Wind direction (deg),
Wind speed (m/sec),
Mixing height (m),
• Stability class (dimensionless),
Standard deviation of elevation angle (radians),
• Standard deviation of azimuth angle (radians),
• Ambient air temperature (K), and
• Anemometer height (m).
The user has the option of updating the meteorological
information after each meteorological time period. The location
and height of each receptor must be indicated. If dispersion is
characterized by the on-site scheme, then the standard deviations
of the azimuth and elevation angles are required.
The following information is required by the plot routines:
• Type of plot desired,
• Location of concentration versus time plots, and
Plott ing grid.
The plot routines were developed on a UNIVAC 1110 and use CALCOMP
plotting software.
6-84
-------
SECTION 3
FEATURES AND LIMITATIONS
Several requests to the Environmental Operations Branch for
assistance in modeling the air quality downwind of incineration
ships stimulated the development of INPUFF, a model capable of
simulating a moving point source in a spatially variable wind
field. The model also possesses the following features which
increase its flexibility and range of application:
Optional stack-tip downwash,
Wind speed extrapolated to release height,
• Temporally variable source characteristics,
• Temporally and spatially variable wind field,
Up to 25 receptors,
Some consideration of terrain effects through
the wind field, and
Optional graphics display.
The implied modeling scale is from tens of meters to tens of
kilometers. INPUFF is capable of addressing the accidental
release of a substance over a short time period, or of modeling
the more typical continuous plume from a stack.
Although INPUFF has several advantages over its continuous
plume counterpart, it still retains several limitations,
including:
Wind direction constant with height,
• No consideration of pollutant removal or
chemical reactions, and
No explicit treatment of complex terrain.
Also, INPUFF has no provision for calculating the effect of a
large area source or multiple point sources.
6 6-84
-------
The table below presents a comparison of
those of other air quality models.
INPUFF features with
TABLE 1. A COMPARISON OF INPUFF TO OTHER COMMONLY USED
AIR QUALITY MODELS
X -- used by model
O -- opt i ona 1
MODEL TYPE
Gaussian plume
Gaussian puff
GRID SIZE
AVERAGING PERIOD
Hour
3-hour
24-hour
Annua 1
TYPE OF SOURCES
Single stack
Mu 1 1 ip 1 e s tacks
Area sources
Line sources
RECEPTORS
Number of
Cartesian coordinates
Cartesian coordinates w/ elevations
Polar coordinates
Polar coordinates w/ elevations
Program generated grid
METEOROLOGICAL DATA
RAMMET preprocessor
\1ESOPAC preprocessor
STAR file
User spec i f led
Program generated
POLLUTANT
Non-reac t i ve
Hal f-1 i fe
PLUME RISE
Stack tip down wash
Gradual plume rise
Buoyancy- induced dispersion
TERRAIN ADJUSTMENTS
I
N
P
U
F
F
X
100
0
O
O
X
25
X
X
X
O
0
P P P
T T T
P D M
L I T
U S P
XXX
R
A
M
X
M
P
T
E
R
X
c
R
S
T
E
R
X
V
A
L
L
E
Y
X
P
A
L
X
M
E
S
O
P
U
F
F
vl
X
1600
XXX
X
XXX
25
50 30
X5
X
X X
X
XXX
O
0 X
O
O
X
O
0
X
250
100
180
X
X
O
X
0
X
O
O
O
O
0
O
0
O
X
250
180
X
X
X
X
X
0
X
O
O
O
O
0
X
X
X
X
X
19
180
X
X
X
O
O
0
O
0
X
50
50
112
X
X
X
X
0
O
O
0
X
O
0
X
99
99
99
99
X
X
X
X
O
O
0
X
10
X"
X
X
X'
0
(1) Specially suited for long range
t ranspor t .
(2) Co-located stacks.
(3) Total of 50 point and/or area
sources .
(4) Concentrations computed at each
grid point and at up to 50
user-specified receptors.
(5) Concentrations computed at
fixed downwind distances.
(6) Linear conversion of SO2 to
SO4 and linear deposition of
SO2 and SO,.
6-84
-------
SECTION 4
BASIS FOR INPUFF
GAUSSIAN PUFF METHODOLOGY
A graphical representation of the INPUFF model is given in
Figure 1. Here the first puff (the puff with the longest
trajectory) was first exposed to east-southeast winds, followed
by slightly stronger winds from the south and the
south-southeast. The second puff was released at the time the
winds shifted from east-southeast to south. The third puff was
released when winds were from the south-southeast. The stability
conditions need not be equal for the various time steps, though
in the figure, stability is shown to be fairly constant with time
(i.e., the rate of puff growth is constant over the time frame).
INPUFF assumes ax = av, thus puffs remain circular throughout
their lifetime. Puffs A, B, and C represent the location
three emitted puffs at time
of the
In Gaussian-puff
a series of puffs
algorithms, source emissions are treated as
emitted into the atmosphere. Constant
conditions of wind and atmospheric stability are assumed during a
time interval. The diffusion parameters are functions of travel
time. During each time step, the puff centers are determined by
the trajectory and the in-puff distributions are assumed to be
Gaussian. Thus, each puff has a center and a volume which are
determined separately by the mean wind, atmospheric stability,
and travel time.
PLUME RISE
Plume rise is calculated using the methods of Briggs (see
Section 5). Although plume rise from point sources is usually
dominated by buoyancy, plume rise due to momentum is also
considered. Building downwash, buoyancy-induced dispersion, and
gradual plume rise are not presently treated by INPUFF.
Stack-tip downwash (optional) can be considered using the
methods of Briggs. In such an analysis, a height increment is
deducted from the physical stack height before momentum or
buoyancy rise is determined. Use of this option primarily
affects computations from stacks having small ratios of exit
velocity to wind speed.
6-84
-------
B
SOURCE
Figure 1. Gaussian puff model.
6-84
-------
DISPERSION ALGORITHMS
Three dispersion algorithms are used within INPUFF for
dispersion downwind of the source:
P-G scheme as discussed by Turner (1970),
On-site scheme formulated by Irwin (1983), and
Long travel time scheme.
The user has the option of choosing either the P-G or the on-site
algorithm (for short travel time dispersion) and specifying when
the long travel time dispersion parameters are to be implemented.
Dispersion downwind of a source, as characterized by the P-G
scheme, is a function of stability class and downwind distance.
Stability categories are commonly specified in terms of wind
speed and solar radiation. The on-site dispersion algorithm is a
synthesis of Draxler's (1976) and Cramer's (1976) ideas and
requires specification of the variances of the vertical and
lateral wind directions. The third dispersion scheme is used in
conjunction with the other two and is for long travel times in
which the growth of the puff is proportional to the square root
of time.
10 6-84
-------
SECTION 5
TECHNICAL DESCRIPTION
This section presents the mathematical formulation of the
Gaussian-puff model.
GAUSSIAN PUFF EQUATIONS
The concentration, C, of a pollutant at x, y, z from an
instantaneous puff release with an effective emission height, H,
is given by the following equation:
C(x,y,z,H) = Q expr-j/x-utX2] expr-!/_y\21
(2TT)3'2axayoz [ 2\ a* / J L 2 \°y / J
jexp["-i/z+H\2l + exPr-i/z-H\21) (1)
—
Since each puff is free to move in response to changing wind
speed, u, and is not constrained to a single centerline, the
diffusion parameters are given as functions of travel time, t,
rather than of downwind distance.
Following the puff and assuming ax equals ay, expressed as
ar, where r = /(x-ut)2 + y2, the puff equation can be rewritten
as follows:
C(r,z,H)= Q expr-j./^Y1iexpr-l/^+H\2"|+expi"-jL/z-jl\2'|) (2)
When az becomes larger than eight tenths of the mixed depth
layer, L, the puff is assumed to be well mixed and the
concentration equation is expressed as,
C(r,z,H)= Q expf-j./_r_\2"| for az > 0.8L. (3)
j./_r_\2]
2\CTr/ J
The total contribution from all the puffs is summed at each
receptor after each time step.
Although a Gaussian-puff model, such as INPUFF, is useful in
estimating pollution dispersion under unsteady and nonuniform
flow, it has several limitations:
11 6-84
-------
(1) Pollution dispersion within the puff is assumed to be
Gaussian and meteorological conditions within a time step
are assumed to be spatially and temporally uniform. These
assumptions may cause significant error in estimating
concentrations, especially at long travel distances.
(2) The diffused material is assumed to be stable and to remain
suspended in the atmosphere over a long period of time.
Chemical reactions, atmospheric removal, and other nonlinear
processes are not handled by INPUFF.
(3) Data for puff diffusion is sparse and there is no ordering
of the a curves by stability; therefore, many Gaussian-puff
models use plume a's. However, similarity theory for puff
diffusion (Batchelor, 1952) suggests that there is a region
in which puff growth is greater than plume growth. For
downwind distances where travel time is larger than sampling
time, the use of plume a's in a puff model may be
inappropriate. However, as long as the variations in
meteorological conditions are not simulated to any finer
resolution than 3 to 10 minute periods, the use of plume
characterizations of dispersion may still be reasonable.
(4) As mentioned, the primary purpose of the integrated puff
model is to simulate a continuous plume. Plume diffusion
formulas apply to continuous plumes, where the sampling time
is long compared to the travel time from source to
receptor. Since INPUFF uses the plume characteristics of ay
and az, one would expect that the concentration estimates
from INPUFF would yield the best agreement with observations
if the travel time was short compared to the sample
duration of the concentration estimates. Since this
assumption is violated, the model estimates relate more to
the average of many realizations of the same experiment,
recognizing that the correspondence of any one experiment
may differ greatly in comparison to the average obtained
from many experiments.
(5) Given the complex nature of the wind field, sampling the
flow so that it can be completely defined from a
mathematical point of view is impossible. There can always
be any number of solutions which could stem from one initial
state, while satisfying all other requirements.
The most important difference between Gaussian-plume models
and INPUFF is that INPUFF can handle changing meteorological
conditions, whereas typical Gaussian-plume models assume spatial
and temporal uniformity in the meteorology.
12 6-84
-------
INPUFF is similar to grid models in that both require the
construction of a spatial and temporal wind field. Grid models,
which solve the diffusion equation via finite-differencing
schemes, are of an Eulerian type, whereas INPUFF is of a
Lagrangian type. Grid models can incorporate nonlinear
processes, such as chemical reactions. However, they require
more computer processing time. Another problem with grid models
is that because of computer memory limitations, the grid spacing
is usually too large to incorporate physical processes such as
plume rise and stack downwash.
PLUME RISE
Plume rise from point sources is calculated using the methods
of Briggs (1969, 1971, 1973, and 1975). These equations are
based on the assumption that plume rise depends on the inverse of
the mean wind speed and is directly proportional to the
two-thirds power of the downwind distance from the source, with
different equations specified for neutral or unstable conditions
and for stable conditions. Only the final rise equations are
summarized below. The reader is referred to Appendix A for the
details of the formulation.
For unstable or neutral atmospheric conditions, the downwind
distance of final plume rise is
Xf = 3.5x*,
where
x* = 14F5/8 for F < 55 mVsec3
and
x* = 34F2/5 for F > 55 m*/sec3.
The final plume rise under these conditions is
H = h' + [1.6FI/3(3.5x*)2/3/u(hH. (4)
For stable atmospheric conditions, the downwind distance of
final plume rise is
Xf = 0.0020715u(h)s-1/z
where
s = gO6/3z)/T.
P1 ume rise is
H = h' + 2.6{F/[u(h)s j}1'3 for windy conditions (5)
13 6-84
-------
and
H = h
4F
1/"-3/8
for near-calm conditions. (6)
The lower
(5 and 6)
and units
in Table 2.
of the two values obtained from the above two equations
is taken as the final effective height. Definitions
of variables mentioned in this section are summarized
TABLE 2.
Symbol
F
g
H
h'
s
T
u(h)
Xf
X*
DEFINITION OF VARIABLES USED IN PLUME RISE
Def ini t i on
Buoyancy flux parameter
Acceleration due to gravity
Effective height of plume
Stack height adjusted for stack downwash
Stability parameter
Ambient air temperature
Wind speed at stack top
Distance to final rise
Distance at which atmospheric turbulence
begins to dominate entrainment
EQUATIONS
Units
mH/sec3
m/sec2
m
m
sec"2
K
m/sec
m
m
DISPERSION ALGORITHMS
The primary purpose of the- integrated puff model is to
simulate a continuous or semi-continuous plume for varying
meteorological conditions. The vertical and lateral dispersion
parameters for continuous plume dispersion models are used in
INPUFF. Under steady meteorological conditions, the output
concentrations of INPUFF should, all other factors such as plume
rise being equal, approximate the results calculated by a
Gaussian-plume model such as PTPLU. To demonstrate this,
concentration estimates of INPUFF and a modified version of
PTPLUT are compared. The meteorology used in this comparison is
as follows:
• Wind speed -- 3 m/sec,
Wind direction -- 270°,
Mixing height -- 1500 m, and
Stabi1ity class -- C.
t PTPLU modified to output concentrations at various downwind
distances instead of maximum concentration and distance to the
maximum.
14
6-84
-------
INPUFF was executed for
steady-state conditions.
a
2-hour simulation to bring about
Table 3 sumnarizes the results. The last column shows the
percent difference in the computed concentrations for the two
models. Although they differ by a factor of five at receptors
close to the source, the percent difference decreases to less
than 5% near the maximum concentrations. The results show that
INPUFF can indeed simulate a continuous plume.
TABLE 3. COMPARISON OF INPUFF AND PTPLU
Downwi nd
distance
(km)
2.154
2.448
2.783
3.162
3.594
4.084
4.642
5.275
5.995
6.813
7.743
8.799
10.000
11.365
12.915
Concentra
(ug/m3
INPUFF
0.20
1.06
4.27
13.19
32.41
65.07
108.50
158.00
202.90
233.30
252.10
254.80
245.30
224.80
197.30
t ion
)
PTPLU
0.04
0.44
2.60
10.12
28.24
60.62
105.75
156.52
203.42
238.55
257.88
261.42
251.87
233.15
209.25
Difference t
(%)
400.00
140.91
64.23
30.34
14.77
7.34
2.60
0.95
- 0.26
- 2.20
- 2.24
- 2.53
- 2.61
- 3.58
- 5.71
t Difference(%) = [(INPUFF - PTPLU)/PTPLU] 100
Three dispersion algorithms are incorporated within the model
to account for initial dispersion, short travel time dispersion,
and long travel time dispersion. The initial dispersion
algorithm handles the finite size of the release through the use
of initial dispersion parameters. Once the puff leaves the
source its growth is determined by the short travel time
dispersion algorithm. This algorithm has two schemes: the
Pasqui11-Gifford scheme which characterizes dispersion as a
function of downwind distance and the on-site scheme which
characterizes dispersion as a function of travel time. For long
travel time, a dispersion algorithm that allows the puff to grow
as a function of the square root of time is used.
15
6-84
-------
Initial Dispersion
The initial dispersion of the plume at the source is modeled
by specifying the initial horizontal and vertical dispersion
parameters, aro and ozo. For tall stacks these parameters,
generally, have little influence on downwind concentrations.
However, if the source is large enough or close enough to the
ground, then initial size is important in determining ground
level concentrations near the source. For a source near the
ground, the initial horizontal dispersion can be calculated by
dividing the initial horizontal dimension of the source by 4.3,
and the initial vertical dispersion parameter is derived by
dividing the initial height of the source by 2.15. This method
of accounting for the initial size of near ground level release
gives reasonable concentration estimates at downwind distances
greater than about five times the initial horizontal dimension of
the source.
Short Travel Time Dispersion
Dispersion downwind of the source can be characterized by the
P-G scheme, which is a function of stability class and downwind
distance, or by the on-site scheme, which is a function of travel
time.
Pasqui11-Gifford Scheme
The P-G values that appear as graphs in Turner (1970) are
used in the model. However, for neutral atmospheric conditions
two dispersion curves as suggested by Pasquill (1961) are
incorporated into the model. Dispersion curves Dl and D2 are
appropriate for adiabatic and subadiabatic conditions,
respectively. The D2 curve is used in Turner (1970) for neutral
practical point of view, since temperature
be available we refer to the Dl and D2 curves
P-G stability classes are numerical inputs
Stability classes A through D-day are
classes D-night through F are specified by
conditions. From a
soundings may not
as D-day and D-night.
to the puff model .
specified by 1-4, and
5-7, respectively.
On-site Meteorology Scheme
The o-curves of the P-G scheme above are based on data of
near-ground level releases and short-range dispersion studies.
These data are used to extrapolate the P-G curves to high release
heights and far receptor distances. In view of this, INPUFF has
an option of using on-site meteorological data to estimate
dispersion. This scheme is a result of the recommendations of
the American Meteorological Society's workshop on stability
classification schemes and sigma curves (Hanna et
Irwin (1983) proposed characterizing Oy and a.
similar to Cramer (1976) and Draxler tl976).
deviation of the crosswind concentration
al.f 1977).
in a manner
The standard
distr ibuti on
1 S
16
6-84
-------
= avtf,
(7)
where av is the standard deviation of the horizontal crosswind
component of the wind, t is the downstream travel time of the
pollutant, and fy is a nondimensional function of travel time.
The standard deviation of the vertical concentration
distribution, az, for an elevated source, when oz is less than
the source height, is
z, (8)
where aw is the standard
the wind, and fz is
dependent upon travel time. The nondimensional
fz were characterized by Irwin (1983) as
deviation of the vertical component of
a nondimensional function, primarily
functions f., and
fy = !.
fz = 1,
0.9(t/1000)1/2], (9)
for unstable conditions
and
fz = !./[! + 0.9(t/50)1/2] for stable conditions. (10)
Besides the P-G stability class,
which are assumed to be typical of
the scheme
condi t ions
and
requires i
,>w, wiiiun are ussuineu to ue typical 01 conditions at final plume
height. For small angles, av = aau and aw = aeu where u is the
wind speed at measurement height and aa and ae are the standard
deviations of the horizontal and vertical wind angle,
respectively. The puff model requires aa and ae as data input
and computes av and aw.
Long Travel Time Dispersion
That the dispersion parameters used in INPUFF satisfy the
diffusion theory developed by Taylor (1921) is desirable. Taylor
showed that for an ensemble average of particle displacements
during stationary and homogeneous conditions, the dispersion
parameters can be written as,
Td t
= 2(vw)'2J J
00
R(r)dTdt
(11)
where R(T) is
component of
variances of the lateral
velocity, respectively;
the
horizontal and vertical
respectively instead of
the Lagrangian autocorrelation of
the wind velocity fluctuation;
or vertical components of
and Td is the diffusion
appropriate
(vw)'2 are the
vertical components of the
Td is the di ffusion time.
diitusions, v'2 and w12 are
The autocorrection starts
(v'w1 ) z.
and approaches 0 for large diffusion time. Therefore
wind
For
used
at 1
from
17
6-84
-------
Eq. 11, while the growth of the puff is linear with time near the
source, the growth becomes proportional to the square root of
time at large distances. In the model, after the puff has
attained a specified horizontal dimension, the algorithm
automatically goes to a long travel time growth rate proportional
to the square root of time. The size of the puff at that time is
specified by the user. For example, the user may decide that
when or for the puff is greater than 1000 meters the long travel
time dispersion parameters should by utilized. A very large
SYMAX value (see page 45) results in the long travel time code
not being executed.
MIXING HEIGHT
Depending on the stack height, plume rise, and height of the
mixing layer, the puffs can be above or below the mixed depth
layer, L. If the puffs are above L then there are two cases that
govern their growth. Initially the puffs are allowed to grow
according to the P-G, F curve, or if the on-site scheme is used,
the puffs are restricted to a vertical growth rate characterized
by aw = 0.01 m/sec. After the puffs attain a given size of ar
(not actual puff size) specified by the user, the growth rate is
specified by the /I.
When the puffs are below L, then there are four cases that
must be considered. Cases one and two are puffs which are not
well mixed vertically and whose growth rates are characterized by
the short travel time sigmas or by /t. Cases three and four are
puffs that are well mixed vertically and whose growth for or is
for short travel times or according to /"t. During the modeling
simulation, every puff is given a key to indicate whether it is
above or below L and whether its growth rate is characterized by
the short travel time sigmas or by /T.
In the modeling design, puffs are allowed to change their
dispersion keys. When the height of L becomes greater than the
puff height, the puffs are allowed to grow at the rate
characterized by surface measurements. Normally this is a
neutral or unstable situation. This transition period is likely
to occur in the morning hours. In the afternoon, despite the
decay of active mixing, a puff remains well mixed through the
maximum mixing lid as shown in Figure 2. The maximum height of L
is stored for each puff and is never allowed to decrease. This
method assures that concentration does not increase with downwind
distance or travel time, so as to violate the second law of
thermodynami cs.
18 6-84
-------
1000-1
TIME OF DAY (hours)
PUFF ELEMENTS
LEGEND
MIXING LID —
STORED
— MAXIMUM
MIXING LID
Figure 2. Effect of variable mixing height on puff dispersion.
19
6-84
-------
ATMOSPHERIC STABILITY
As discussed earlier, short travel time dispersion can be
characterized by two schemes, the P-G scheme and the on-site
scheme. The P-G scheme uses the empirical P-G curves and
stability classification to estimate dispersion coefficients
(Turner, 1970), whereas the on-site scheme relates diffusion
directly to turbulence. If on-site meteorological data are not
available, only the widely used P-G scheme can be adopted. If
on-site meteorological data are available, either scheme can be
used.
INPUFF's on-site scheme adopts Irwin's algorithm (1983) in
characterizing ay and az. This scheme essentially requires
information on the standard deviations of horizontal (aa) and
vertical (ae) wind fluctuations and wind speed at measurement
height. Stability is classified as stable or unstable from the
near-surface data for temperature difference, Richardson Number,
or stability parameter.
GRIDDING SCHEMES
INPUFF requires a meteorological preprocessor to compute wind
speed and direction at each grid square. The user is required to
specify the format of the meteorological data file. The
coordinate and size of each grid square, as well as the extent of
the meteorological region, must be defined in the input. The
modeling region need not be the same as the meteorological
region, but the southwest corner of both should coincide. If the
meteorological region is smaller than the modeling region and the
puffs travel outside of the meteorological region, then they are
advected according to their last wind speed and direction. If
the meteorological region is larger than the modeling region and
the puffs travel outside the modeling region, they are eliminated
from further consideration. The source must stay within the
modeling region; otherwise, all puffs are eliminated.
To improve the spatial resolution of the concentration
pattern, receptors in INPUFF are specified by the user. The
resolution of the receptors can be more detailed than that of the
meteorological grid. The receptors may be placed independent of
the meteorological grid. Figure 3 illustrates a possible
arrangement of the modeling region, meteorological grid, and
receptor locations. In this example the receptors are
concentrated along part of the puff trajectory with a spatial
resolution two times finer than the meteorological grid.
20 6-84
-------
LU
z
Q
o:
O
O
o
O
CO
i
oc
o
MODELING
REGION
METEOROLOGICAL
GRID
RECEPTOR
GRID
X (km)
EAST-WEST COORDINATE
Figure 3. A possible arrangement of modeling and receptor grids.
21
6-84
-------
SECTION 6
EXAMPLE PROBLEMS
In this section, problems are provided to illustrate
different modeling scenarios and to demonstrate several unique
features of INPUFF. Details concerning input and output of the
first two example problems are discussed in Section 11 after the
reader has become familiar with INPUFF input data preparation.
EXAMPLE 1 -- MOVING SOURCE
This example uses a unique feature of INPUFF that allows the
source to move at a constant speed and direction over a specified
time. Figure 4 shows the source path and receptor locations.
The source is initially southwest of the receptors and travels
due east remaining south of all receptors. Southerly winds at
3.5 m/sec are observed and the atmosphere is slightly unstable.
Twenty minutes into the simulation the source assumes a northeast
heading. Atmospheric conditions become neutral, wind speed
increases to 4 m/sec, and wind direction changes slightly from
180° to 170°. The stack parameters of the source are as follows:
• Emission rate -- 600 g/sec,
Stack height -- 30 m,
Stack gas temperature — 390 K,
Stack gas velocity -- 15 m/sec, and
Stack diameter -- 2 m.
The impacts at the receptors are outlined in Table 4. As
shown in the table, INPUFF provides average concentrations for
each meteorological time period and for the total simulation
time. As expected, impacts are greatest at the western receptors
(1, 2, 5, and 6) during the first meteorological period and to
the eastern receptors (3, 4, 7, and 8) during the second
meteorological period.
22 6-84
-------
3-
N
t
E 9-
Receptor locations
.7
.3
.2 .4
Source direction 90°
-End of first
meteorological period
km
Figure 4. Source path for example 1
23
6-84
-------
TABLE 4. COMPUTED CONCENTRATIONS FOR EXAMPLE 1
Receptor
number
1
2
3
4
5
6
7
8
0-20 min.
135
167
22
<1
181
221
4
<1
Ambient concentrations (ug/m3)
ave . 20-40 min. ave. 40
<1
8
123
13
<1
2
180
13
min ave.
68
87
72
7
90
111
92
6
The input stream and
provided in Section 11. The
also demonstrated in Section
output listing for this problem are
plotting features of the model are
11.
EXAMPLE 2 -- LOW LEVEL SOURCE WITH LOW WIND SPEED CONDITIONS
This problem illustrates the model simulation of a low level
release during conditions of light and variable winds. Another
feature highlighted in the problem is that of temporally variable
source characteristics.
Twelve periods of 10-minute duration are used to simulate a
2-hour release. Both meteorology and source characteristics are
updated every 10 minutes. The wind speeds are light at 0.5
m/sec, and wind direction fluctuates from 145 to 210
dispersion measurements of aa and ae
in the simulation. Values of other
parameters are listed below:
On-s i te
are available and are used
pertinent meteorological
Mixing height -- 5000 m,
aa -- 0.393 radians,
ae -- 0.035 radians, and
• Temperature -- 290 K.
The source-receptor geometry shown in Figure 5 was
based on the observed southeast to south-southwest
Receptors are located along two radial arcs approximately
and 1.0 km from the source. Figure 6 shows how the
strength decays with time. Initially the emission rate
g/sec, but by the 12th period it has dropped to 12 g/sec.
chosen
winds .
0.5 km
source
is 825
24
6-84
-------
en
e
*t
A
O>
V>
to
o
3
A
A
*»
0)
c.
to
I
EMISSION RATE (g/s)
O
•»
A
X
O
to
OQ
C
A
O
C
O
A
I
A
O
A
-0
TO
A
I
o
•t
A
X
P
A
to
NORTH-SOUTH COORDINATE. Y (km)
O -t M U «
m
v>
m
(0
o
g
z
m
x
Sr
W-
oo
-------
Average concentrations at each receptor for the simulation
time are listed in Table 5. As expected, impacts are greatest at
receptors (3 and 8) due north of the source.
TABLE 5. COMPUTED CONCENTRATIONS FOR EXAMPLE 2
Receptor
number
1
2
3
4
5
6
7
8
9
10
2-hour average concentrations
(ug/m3)
13
536
4846
323
4
<1
180
16300
44
<1
The input stream and abridged output listing for this problem
are provided in Section 11. The plotting features of the model
are also shown there.
EXAMPLE 3 -- VARIABLE WIND FIELD
The user-specified horizontal wind field option is exercised
in this example. The topographic relief of the modeling region
is illustrated in Figure 7, along with the location of the point
source and the receptors. Characteristics of the stack are as
follows:
Source strength -- 700 g/sec,
Stack height -- 75 m,
Stack gas temperature -- 455 K,
Stack gas velocity -- 16 m/sec, and
• Stack diameter -- 3 m.
The boundaries of the region in Figure 7 are also the limits of
the modeling and meteorological regions for the INPUFF
s imulat ion.
26 6-84
-------
20-
18-
16-
14-
12-
J,oJ
8-
6-
4-
2-
0
© SOURCE
• RECEPTOR
ELEVATION CONTOUR (m)
• 100
©
I
2
6 8
T
10
12
km
i
14
16 18 20 22
I
24
Figure 7. Topographic relief of the modeling r.egion in example 3
27
6-84
-------
The hourly meteorological data observed at the source are
given below.
TABLE 6. HOURLY METEOROLOGICAL DATA FOR EXAMPLE 3
Hour Direction Speed Stability Mixing Height
(degrees) (m/sec) (m)
1
2
3
4
5
6
272
265
272
282
298
315
2.4
2.4
2.6
4.0
4.9
6.0
F
F
E
D, day
C
C
200
200
200
400
1000
1500
Wind fields in this example were generated using a nondivergent
wind model. Wind fields for the 6 hours modeled, along with the
corresponding puff trajectories, are illustrated in Figure 8. In
the first 3 hours of the simulation, the terrain feature causes
distortion of the wind field and results in meandering of the
puffs. As the atmosphere becomes less stable (hours 4, 5, and
6), flow over the terrain occurs and the field is more uniform.
The wind fields of Figure 8 corresponding to the hours in
question illustrate this point. The puffs released in the last 3
hours follow a straight path.
The pollutant concentration field is provided in Figure 9,
with isopleths representing 6-hour averages. Notice how the
contours are "bent" by the terrain feature. In general, the flow
carried the puffs south of the peak. Thus peak concentrations
occurred along the southward facing slope.
28
6-84
-------
6Z
oq
c
T
(0
oo
Q.
C/l
D
3
a
•o
c
(B
u^.
(B
a
tn
n
a>
a
o
c
o
-5
(0
X
0
•o
0)
cc
c
* -
01 —
s-
3 M •
s-
s-
M _
8-
km
i i i i i i i i i i
* * * 4 * * * \ \ \
4 * 4 I * * * \ \ \
44f*f**\\\
4 4 f x > » A x \ \
4 * * / ^ # \ \ I \
4 4 1 M / \ M \
4 1 1 i \ \ 1 M \
4 4 4 \ \ \ | \ \ \
i i i \ \ \ * ; i i
i ; ; i ^ ^ i i ; i
i i i ; ; ; » i ; *
4 4 * 1 1 i * * \ \
i
o
c
Km
c
M-
»-
o~
3 M~
s-
o~
M_
M
»J
'?r???rr?f?
i
JSi
*&
HI
/K^^
K^=^
<
M «
'' ,. .
•
.-
S -
3 **
s-
8-
M _
M
km
i i i i i i i i i i
* / t 4 * * * I 4 *
* / * / * * * 4 I I
* 4 * / > ^ * * * *
f / / x ' * \ \ \ *
I 4 I / ' / t \ 4 »
4 4 1 M / \ 1 1 »
4 4 4 i 1 I i \ 4 4
4 4 4 * \ 1 / 1 4 4
4 4 4 \ ^ 4 / / 4 4
4414444444
4444*44444
4444444444
o
c
Km
<
M-
.-
s~
x- -*_
3 M
S~
o~
M
M
M_
sYr??frr?f?
J
w
-------
(panuijuoo) -g
ZZ OZ
91
Zl 01 8
I ' '
9 » Z °
I I I
Z 0
-Z
->
-8
-0
-«H
-Zl
•«
-91
-81
-OZ
»Z ZZ OZ 81 9t W Z» 0» B 9
I - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1-0
-Z
-»
-9
-8
^^^^^^_^^^^^-"
-Zl
-»l
-91
-Bl
v nnoH _ l°*
Ult\
n ZZ OZ 81 91 H Zl 01 8 8 f Z 0
I - 1 - 1 - 1 - 1 _ I _ I _ I _ I _ I _ I I I 0
-I
-B
-<»
-Zl
-w
-81
-81
e HHOH oz
oo
i
o
CO
»Z ZZ OZ 81 91 »l Zl 01 8 9 » Z 0
I i I i I I I ] I I 1 I
-Z
-8
-8
-Ol
-Zl
-W
-91
-81
-OZ
-------
o g g » « o
m
cc
15
O
I
. O
«
-S
I I I
LU>|
5 E
- S
tr
O
i
xxxxx
xxxxx
xxxxx
X XX XX
xxxxx
xxxxx
XX XXX
xxxxx
xxxxx
xxxxx
XXX XX
xxxxx
xxx
xxx
xxx
xxx
xxx
XXX
xxx
XXX
xxx
xxx
xxx
xxx
XX
XX
XX
XX
XX
XX
XX
X X
XX
X X
XX
X X
3
C
c
o
t)
oo
3
bo
r- S
_ o
S 2
31
6-84
-------
20-
18-
16
14-
12-
8-
6-
4-
2-
0
T
2
© SOURCE
• RECEPTOR
— CONCENTRATION ISOPLETH (jitg/m3)
T
4
T
6
"T
8
~r
10
12
km
~1 1 1 1 1 1
14 16 18 20 22 24
Figure 9. Pollutant concentration field of example 3.
32
6-84
-------
SECTION 7
COMPUTER ASPECTS OF THE MODEL
INPUFF
This section discusses the general framework of INPUFF. The
section is intended to give the reader a general knowledge of the
computer program, rather than a detailed description of each
subroutine. The general flow of INPUFF, the structure of the
computer subroutines and functions, and a brief description of
each subroutine and function are included.
The main routine reads the following types of information:
• Options to be exercised during program execution,
Simulation information and puff characteristics,
Specifications of the modeling region,
• Source characteristics,
• Receptor coordinates, and
• Meteorological data.
INPUFF is a single source model that permits source
characteristics to be updated every meteorological period. The
meteorology during the modeling exercise can be specified by up
to 144 equal length meteorological periods. Concentration
estimates can be made for 25 locations. If the user does not
apply the wind field option, all input data are read by the main
routine, which also produces all printed output.
Figure 10 shows the structure of the
functions. PUFF is the main routine that reads
stores the appropriate data in common
subroutines. A brief description of the
subroutines, and functions follows.
subroutines and
input data and
with the other
main program,
33
6-84
-------
PUFF«
PLMRS
PGSIG
VTIME — JSISIG
ADVECT*
PROCES
-XVY
•XVZ
-PGSIG
-LTSIG
•VTIME — JSISIG
•JSISIG
CONCEN
Indicates optional call
Figure 10. Structure of INPUFF,
6-84
-------
PROGRAM MODULES
PUFF
inputs except
exerci sed.
PUFF: PLMRS,
The main
the initial
puffs, and
-- PUFF is the main program that reads all
wind field data, if that option is to be
The following subroutines are called by
PGSIG, VTIME, ADVECT, PROCES, and CONCEN,
program generates the puffs, computes
dispersion and plume rise, advects the
combines the puffs where appropriate. PUFF prints out
the input data and the concentration estimates at each
receptor for each time period.
-- This subroutine is called by PUFF and computes the
plume rise based on Briggs1 equations.
-- This subroutine is called by PUFF and PROCES to compute
the P-G dispersion parameters.
-- This subroutine is called by PUFF and computes the
virtual time necessary to account for the current size
of the puff using the on-site scheme. The process is
analogous to the virtual distance concept. The only
subroutine called by VTIME is JSISIG.
ADVECT — This subroutine is called by PUFF if the user-supplied
wind field option is exercised (i.e., IADT = 1).
ADVECT reads the gridded wind field data from unit 21,
and computes the appropriate wind speed and direction
for each puff.
PROCES — Called directly by PUFF, the major functions of PROCES
are to: determine which dispersion routine is called
for each puff, assign dispersion keys (KEYP) for each
puff, account for the effect of the mixed depth layer
PLMRS
PGSIG
VTIME
for each puff, and compute ar and
PROCES calls subroutines PGSIG,
JSISIG, and functions XVY and XVZ.
oz for
LTSIG,
each puff.
VTIME, and
CONCEN --
This subroutine
concent rat i on
locat i on. Th i s
are within
each other
is called by PUFF and
from each puff for
subroutine also combines
computes the
each receptor
puffs if they
a user supplied distance (ALPHA times ar) of
XVY
XVZ
This function is called by PROCES and calculates the
virtual distance necessary to account for the initial
crosswind dispersion using the P-G scheme.
This function is called by PROCES and calculates the
virtual distance necessary to account for the initial
vertical dispersion using the P-G scheme.
35
6-84
-------
LTSIG -- This subroutine is called by PROCES and computes the
dispersion parameters for long range transport.
JSISIG -- This subroutine is called by VTIME and PROCES and
computes the dispersion parameters for short travel
time transport using the on-site scheme.
Subroutines PLMRS, PGSIG, ADVECT, CONCEN, XVY, XVZ, LTSIG, and
JSISIG are lowest-level subroutines, i.e., they do not call any
other routine.
The table below gives the input/output units used by the
model.
TABLE 7. INPUT/OUTPUT UNITS USED BY THE MODEL
Unit number Mode Contents
5 Input Program control and input data
6 Output Output listing
21 Input User-supplied wind field
22* Output/input Input data for plotting software
* Output from the main routine and input for plotting routine.
PLOT POSTPROCESSOR
The plot routine reads the following types of information:
Type of plots desired,
Location of concentration versus time plots, and
Plott ing grid.
The above information is read from unit 5. The following
information, generated by the main routine if IP22 = 1 (see Table
8), is read from unit 22:
• Number of meteorological periods,
• Length of each meteorological period,
• Total simulation time,
• Location of each receptor,
• Computed concentrations at each receptor, and
Location of each puff and its ar and az value.
36 6-84
-------
The plot routines were developed on a UNIVAC 1110 and call
CALCOMP plotting software. The main program calls two
subroutines which actually do the plotting. These are PLTOON,
which generates concentration versus time plots at specified
receptor locations, and PLTTRJ, which plots puff trajectories and
receptor locations. The input data for the plot routines are
shown in Table 9 and are described in the next section.
37 6-84
-------
SECTION 8
INPUT DATA PREPARATION
RECORD INPUT SEQUENCE FOR INPUFF
There are eleven record types that are read by INPUFF. Seven
of these are free format input, three are alphanumeric, and one
is a user-specified format. While the free format is very easy
to use, care should be taken to ensure that every variable is
given a value in the correct order. Each variable should be
separated by a comma and should conform to the variable name
type. Three of the eleven input records are optional, depending
on the options exercised on record 2. A brief description of
each input parameter is given in Table 8 with the appropriate
units; the metric system of units is used throughout the model.
Thus horizontal coordinates of source and receptor locations are
in kilometers, temperatures in degrees kelvin, and emission rates
in grams per second. Under the "Formal" column of Table 8, AN
refers to alphanumeric, FF represents free format, and US
indicates user-specified format.
TABLE 8. RECORD INPUT SEQUENCE FOR INPUFF
Record type <5c
Variable Format Variable description Units
Record 1
ALP
Record 2
I OPT
IOPS
AN 80-character title to describe
output
FF Stack downwash option
(1 = use option, 0 = do not use)
FF Update source characteristics
option (1 = source information
update every meteorological
period, 0 = no update)
(continued)
38
6-84
-------
TABLE 8. (continued)
Record type &
Variable Format
Variable description
IADT
IP22
IPCC
FF
FF
IPIC
FF User-supplied wind field option
(1 = use option, 0 = do not use)
Unit 22 output option
(1 = use option, 0 = do not use)
Print out option (1 = puff
information to be printed each
ITIME, 0 = no print out)
FF Print intermediate concentration
option (1 = print out, 0 = no
print out)
Record 3 -- Optional
XSWC FF
YSWC
NUMX
NUMY
DX
DY
Record 4
FRMAT
Record 5
NTIME
ITIME
FF
FF
FF
FF
FF
Optional
AN
Units
East-west coordinate of SW corner km
of wind field grid
North-south coordinate of SW km
corner of wind field grid
Number of grid squares in
east-west direction
Number of grid squares in
north-south direction
East-west width of grid squares km
North-south width of grid squares km
Format of unit 21. Subroutine
ADVECT reads unit 21 according to
the format specified in FRMAT
FF Number of meteorological periods
for simulation
FF Time span of each meteorological sec
per iod
(cont inued)
39
6-84
-------
TABLE 8. (continued)
Record type &
Var iable
Format Variable description
Units
I STEP
INC
I SIM
ISTTIM
IW
DISKEY
NREC
ANHGT
Record 6
XGRID
YGRID
XSIZE
YSIZE
FF Time between puff release (If sec
ISTEP < 0 a value for ISTEP will
be computed based on the stability
class, wind speed, and minimum
distance from the source to the
receptor - GDIS).
FF Time interval at which sec
concentrations are printed out.
INC must be evenly divisible
into ITIME.
FF Time to start concentration sec
calculat ions
FF A multiple of ISTEP for puff
combi nat i on
FF Unit number for write statements
FF Dispersion parameter options
DISKEY = 1 for P-G scheme
DISKEY = 2 for on-site scheme
FF Number of receptors
FF Anemometer height m
FF East-west coordinate of SW km
corner of modeling region
FF North-south coordinate of SW km
corner of modeling region
FF East-west size of modeling region km
FF North-south size of modeling km
region
(cont inued)
40
6-84
-------
TABLE 8. (continued)
Record type &
Var i able
STANUM
ALPHA
SYMAX
Record 7
XSORC
YSORC
QP
HPP
TSP
DP
VSP
VFP
SYOP
SZOP
SDIR
SSPD
Record 8
XREC
Format
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
FF
Variable description
If the puff center is at a greater
distance than the product of
STANUM and ar from a receptor, no
concentration calculation is made
at that receptor for that puff.
Fraction of ar for puff
combi nat i on
Maximum size of ar before going to
LTSIG routine to calculate the
dispersion parameters
East-west coordinate of the source
North-south coordinate of the
source
Emission rate
Height of release
Stack gas temperature
Stack diameter
Stack gas velocity
Stack gas volume flow
Initial ar
Initial az
Source direction, if moving
Source speed
East-west coordinate of receptor
Units
—
m
km
km
g/sec
m
K
m
m/sec
m3/sec
m
m
deg
m/sec
km
YREC
FF North-south coordinate of receptor km
(cont i nued)
41
6-84
-------
TABLE 8. (continued)
Record type &
Var iable
ZREC
Record 9
FRMT
Record 10
WDIR
WSPD
HL
KST
SGPH
SGTH
TEMP
GDIS
Format
FF
AN
US
US
US
US
US
US
US
US
Variable description
Height of receptor
Format of the meteorological
data to be read on record 10
(see page 55)
Wind direction
Wind speed
Mixing height
Stability class (see page 45)
Standard deviation of elevation
angle, cre
Standard deviation of azimuth
angle, aa
Air temperature
Minimum distance from source to
Units
m
deg
m/sec
m
--
radians
radians
K
km
receptor
Record 11 -- Optional
FF Record 11 is optional and is
read if IOPS = 1 (record 2).
If IOPS = 1 then there must be
NTIME-1 number of records in this
group to update source data. The
description of pertinent variables
is given in record type 7 beginning
with QP.
US = User-specified format
Most of the input data are straightforward and typical of the
kind of information required for Gaussian models. However, there
are some input variables which are unique to this code and
require additional explanation to ensure proper assignment of
values.
42
6-84
-------
Records 3 and 4
Records 3 and 4 are read if IADT = 1. The information on
record 3 defines the coordinates of the SW corner of the gridded
region and the size of each grid square. The user must define
the wind speed and direction in each grid square. There are a
few caveats associated with using gridded meteorological data.
The source must stay within the defined region. The
meteorological region defined on record 3 need not be the same as
the modeling region defined on record 6, but the southwest corner
of both should have the same coordinates. If the meteorological
region is smaller than the modeling region and the puffs travel
outside of the meteorological region, then they will be advected
according to their last wind speed and direction. If the
meteorological region is larger than the modeling region and the
puffs travel outside the modeling region, they will be eliminated
from further consideration. Record 4 requires the user to input
the format of his meteorological data file. This file has to be
assigned to unit 21, and is read by subroutine ADVECT according
to the format specified on record 4. This is the only statement
read outside of the main routine. If the option to specify the
wind field is exercised, then the meteorological data read on
record 10 must be appropriate for the grid square that contains
the source. Record 10 must be supplied whether or not the wind
field option is exercised.
Record 5
The data requested on record 5 give the program information
regarding the modeling design. NTIME is the number of
meteorological periods simulated in a run. ITIME is the time
period associated with the meteorological data. For example, if
the meteorological data are recorded in 20-minute averages and
the user wants to make a 3-hour simulation, then NTIME = 9 and
ITIME = 1200 seconds. ISTEP is the time interval between puff
releases. If ISTEP is assigned a negative value the model
computes ISTEP based on the stability class, wind speed, and
minimum distance from source to receptor. The minimum value that
can be assigned to ISTEP is 1 second. When assigning ISTEP for a
moving source, be sure to take into account the path of the
source when computing the minimum distance between source and
receptor (GDIS), specified on record 10. ISTEP should always be
divisible into ITIME and no greater than one-half the value of
INC. INC is a time interval for which intermediate concentration
values are printed out. For example, if ITIME = 1200 and
INC = 300, then four 5-minute average concentration tables are
printed (if IPIC = 1) as well as the 20-minute average
concentration table.
43 6-84
-------
The next two input parameters, ISIM and ISTTIM, are used to
reduce computing time. ISIM is the time when concentration
calculations are to begin. For many cases ISIM is assigned a
value of zero. However, if the minimum source-receptor distance
is large and requires a substantial amount of travel time for the
puffs to reach the receptor, a value for ISIM can be assigned
which would advect the puffs downwind but would delay the
concentration calculations until the current time equaled ISIM.
ISTTIM has to be an integer multiple of ISTEP. For example, if
it is desired to test for puff combination every third puff
release, then ISTTIM = 3. The user should not assign a value of
ISTTIM such that ISTTIM * ISTEP is greater than ITIME.
Record 6
XGRID, YGRID, XSIZE, and YSIZE on record 6 are used to define
the modeling region. STANUM and ALPHA are used to reduce
computation time, and SYMAX is used to specify when the model
should go to long travel time dispersion.
The parameter ALPHA controls when puff combinations take
place. Combinations occur only for adjacent puffs in the release
sequence which have the same dispersion key. A puff can have one
of six possible dispersion keys: (1) puff is below the mixing
height and using short travel time dispersion; (2) puff is using
long travel time dispersion; (3) puff is above the mixing height;
(4) puff is well mixed and using either P-G or on-site
dispersion; (5) puff is above the mixing height and using long
travel time dispersion; and (6) puff is well mixed and using long
travel time dispersion. For instance, suppose two puffs are
adjacent in time and have identical dispersion keys. If ALPHA is
1 then the puffs combine when their centers are within ar of each
other (ar of the younger puff is used for the test). If ALPHA
equals 2, then the puffs combine when their centers are within
2ar of each other. A value of ALPHA equal to 0 results in no
puff combinations. ALPHA can be assigned any value; however, in
practice, ALPHA equal to 1 is a reasonable value for puff
combinat ion.
Upon combining puffs, the position, displacement, and travel
time are averaged between the two puffs. The mass is summed and
the dispersion parameters and virtual position are that of the
younger puff.
The parameter STANUM controls whether or not a concentration
estimate is made for a given puff at a given receptor. STANUM is
the number of crr' s in the horizontal dimension that the puff has
to be from the receptor before a concentration estimate is made
for that puff. For example, if STANUM equals 5, no concentration
estimate is made for any puff at a given receptor if the distance
between the receptor and the puff center is greater than 5ar.
44 6-84
-------
SYMAX is the maximum size of ar for any puff before the
program calls LTSIG to compute the dispersion parameters. SYMAX
can be assigned any size (in meters) depending on how soon the
user wants the model to compute the dispersion parameters as a
function of the square root of time. If it is desired not to
call LTSIG, then a very large value of SYMAX should be assigned.
Record 7
Source information data are specified on this record. The
user should ensure that the coordinates of the source are within
the modeling region. If the source is outside of the modeling
region or moves outside of the modeling region, all puffs emitted
after that point are eliminated from further consideration.
Source direction is specified in degrees from north. A source
direction of 90° means the source is heading east (180° = south;
270° = west; and 360° or 0° = north).
Record 9
The user specified format for record 10 is provided on this
record. The format can be either free or formatted.
Record 10
With the exception of stability class (KST) the variables on
this record are typical of many air quality models. As mentioned
in Section 5, INPUFF considers seven stability categories with
the inclusion of D-day and D-night. Thus stability classes A
through D-day are specified by 1-4, and classes D-night through F
are specified by 5-7, respectively.
INPUT DATA FOR PLOT POSTPROCESSOR
The input data for the plot postprocessor, assigned on four
input records, are read using free format (indicated by an FF in
Table 9). Table 9 shows the input parameters for each record
with the appropriate units. The main routine of the plotting
package reads the input data and the information generated on
unit 22 by the main routine of the puff model. There are two
plots which are optional output in the execution of the plotting
routine. One is a plot of concentration versus time and the
other is a plot of the puff trajectory at the end of each
meteorological period. Either one or both of the plots may be
requested during a given simulation.
45 6-84
-------
TABLE 9. RECORD INPUT SEQUENCE FOR PLOT POSTPROCESSOR
Record type &
Variable Format Variable description Units
Record 1
IPLT FF Plotting options:
1 = plot concentration versus
t ime
2 = plot puff trajectory
3 = plot both
Record 2
IYR
NUMR
ITPT
FF Order of magnitude of
concentration to be plotted on
the y-axis. (Default value is 6)
FF Number of receptors for which
concentration versus time is
plotted
FF Number of periods for which
concentration versus time is
plotted. ITPT must be evenly
divisible into NTIME. (If
ITPT > 999, all periods are
plotted together.)
XSI
YSI
Record 3
IREC
Record 4
XMIN
YMIN
XSIZE
YSIZE
AXL
AYL
FF
FF
FF
FF
FF
FF
FF
FF
FF
Length of x-axis
Length of y-axis
Receptor number for concentration
versus time plots. (NUMR integers
are read on this record.)
East-west coordinate of SW
corner of plotting grid
North-south coordinate of SW
corner of plotting grid
East-west size of plotting grid
North-south size of plotting grid
Length of x-axis
Length of y-axis
in
in
--
km
km
km
km
in
in
46
6-84
-------
On record 2, NUMR is the number of receptor locations that a
plot of concentration versus time is generated. The actual
receptor numbers are read on record type 3. For example, if the
user has made concentration estimates at ten locations and wishes
to see the concentration versus time plots for receptors 1, 3,
and 8, then NUMR = 3 and the array on record 3 is assigned the
values 1, 3, and 8. The third parameter on record 2 is ITPT.
This parameter allows the user to combine meteorological periods
for the concentration versus time plots. If ITPT = 1, then a
concentration versus time plot is generated each ITIME for all
receptors specified on record type 3. However, for ease in
observing the time variations in concentrations, the periods can
be combined. For example, if NTIME = 3 and ITIME = 3600 (i.e., a
3-hour simulation) and a plot of concentration versus time is
desired for the entire 3 hours, ITPT should be set to greater
than 999. ITPT must be evenly divisible into NTIME, or be
greater than 999.
47 6-84
-------
SECTION 9
SENSITIVITY ANALYSIS
This section presents a simple analysis designed to acquaint
the user with the magnitude of changes expected in pollutant
concentrations and CPU time when certain model inputs are
varied. The verification run presented in Section 10 was used as
a basis for this analysis.
PUFF COMBINATION -- ALPHA
Integrated puff models are by their nature computationally
time consuming. To minimize computational time required in the
model, the puffs are combined or deleted, or in certain
situations no computation is made. For instance, if a puff is
not close to a receptor no computations may take place. The
parameter ALPHA controls the rate of puff combinations. If the
value of ALPHA is 1, then the puffs combine when their centers
are within one lateral standard deviation of each other.
As noted in Figure 11, CPU time increases rapidly as ALPHA
approaches zero due to increased number of puffs. Execution time
for ALPHA equal to 0.2 is more than four times longer than for an
ALPHA of 1. CPU time levels off for ALPHA greater than 1.
Increasing ALPHA from 1 to 3 results in only a 50% reduction in
execution time.
The sensitivity of concentrations to ALPHA is shown in Table
10. Varying ALPHA from 0 to 1 has little effect on
concentrations, but for values of ALPHA greater than 1,
computation errors greater than 10% may occur. This result, in
conjunction with decreased computer costs with increasing ALPHA
(see Figure 11), suggests that ALPHA equal to 1 is a reasonable
value for puff combination.
48 6-84
-------
4.0 r
3.0
LU
2
2.0
Q.
o
0.0
1.0 2.0
a
3.0
Figure H. Sensitivity of CPU time to ALPHA.
49
6-84
-------
TABLE 10. PERCENT CHANGE IN CONCENTRATIONS USING DIFFERENT
ALPHA VALUES*
Downwi nd
di c t" o n f» P
(km)
0.5
1.0
2.0
3.0
5.0
10.0
20.0
30.0
50.0
0.4
0
0
0
0
0
0
0
0
0
0.6
0
+ 1
0
+ 1
+ 1
+ 1
0
0
0
ALPHA
1.0
0
+ 4
+ 1
+ 1
+ 1
-1
+ 1
0
0
2.0
0
+ 10
+ 6
+ 8
+ 4
_ i
-2
-5
+ 2
3.0
0
+ 15
+ 9
+ 3
+ 15
+ 4
+ 2
+ 3
-6
Concentrations were compared with those
ALPHA equal to 0.2.
computed with
FREQUENCY OF PUFF COMBINATION — ISTTIM
The parameter ISTTIM controls the frequency at which the
model investigates the puffs to determine if they are to be
combined. Interrogating the puffs every time step is time
consuming in itself; therefore, ISTTIM is specified in multiples
of the time step. For example, if one desires to test for puff
combination every third puff release, then ISTTIM equals 3.
ISTTIM also controls the frequency puffs are interrogated to
determine whether they are outside the modeling region.
Figure 12 shows the relationship between execution time and
ISTTIM. It is apparent from the figure that computer costs are a
minimum when ISTTIM is between 3 and 8. Execution time increases
as ISTTIM approaches 1 because the code interrogating the puffs
is executed more often. It also increases when ISTTIM is greater
than 8 because the program must keep track of more puffs.
PUFF SIZE -- STANUM
The parameter STANUM determines whether or not a
concentration estimate is made for a given puff at a given
receptor. STANUM is the number of lateral standard deviations
defining the maximum puff-receptor distance for concentration
estimates. For example, if STANUM is equal to 5, then
concentration calculations are made only if the puff center is
within 5ar of the receptor.
50
6-84
-------
CPU TIME [STANUM]
CPU TIME [STANUM =10]
P
m
7*
o
OT
(D
3
w
o
w>
O
I
c
OS
I
oo
C
•i
(D
en
(D
3
w
v-
O
•-•,
O
cj
o
»•*
OJ
CPU TIME [ISTTIM]
CPU TIME [ISTTIM=5]
o P
'-» M
O
U
9
o>
z e
2 S
-------
Figure 13 shows the sensitivity of execution time to STANUM.
Execution time is normalized by dividing the CPU time of each run
with the CPU time for STANUM equal to 10. As expected, computer
costs increase as STANUM increases since concentration
calculations are made more often. However, choosing STANUM equal
to 1 or 2 may result in unacceptable underestimates as
demonstrated by Figure 14.
Precision and accuracy in the concentration estimates improve
as STANUM is increased, but the precision gained quickly becomes
unnecessary. In Figure 14, the concentrations at various
downwind receptors are divided by the concentrations computed
when STANUM equals 10. At STANUM equal to 3, the computed
concentrations are almost identical to those computed when STANUM
equals 10. Increasing STANUM from 5 to 10 increases the computer
costs but yields no change in concentrations. Thus choosing
STANUM equal to 3, 4, or 5 should be sufficient for any
simulation case.
SIZE OF MODELING REGION
By defining the modeling region carefully, the user may save
substantial computer costs as illustrated in Figure 15. For
example, it makes little sense to extend the modeling region 50
kilometers downstream of the source when all the receptors are
within 5 kilometers. INPUFF keeps track of all puffs in the
modeling region regardless of their distance from a particular
receptor. It might, nevertheless, be useful to have a large
modeling region under some circumstances, such as in a dramatic
wind shift situation that blows puffs back over the receptors.
52 6-84
-------
CPU TIME [SIZE]
CPU TIME [50 km]
co
c
>-i
(D
CO
(D
9
CO
O
T3
C!
3
5^
*
S
o
z 8
1
5
53
>
o
? 8
3T
o
£
CONCENTRATION [STANUM]
CONCENTRATION [STANUM=10]
9 o p o p
o. a ~J o» «o
b
. . .^
.
•
.
•
-
'
CO
^
z
V.. 2
"X
z
n
•^
1
CO
z
c
2
I ii
I CA
s
5
z
c
S
n
' V
O
-------
SECTION 10
EXECUTION OF THE MODEL AND SAMPLE TEST
INPUFF produces an error-free compile on IBM MVS and UNIVAC
EXEC 8 computers with comparable execution results. A sample job
stream is presented below.
UNIT 5 = DATA
UNIT 22 = OUTPUT
PLOTTER DATA
UNIT 21 = INPUT
WIND FIELD
UNIT 6= PRINTER
EXECUTE INPUFF
JOB CARD
IF IP22 =1
IF IADT = 1
Figure 16. Sample job stream for INPUFF,
54
6-84
-------
Sample test data for model verification are as follows
INPUFF VERIFICATION RUN
1,0,0,0,0,0
2,3600,-1,3600,0,5,6,1,7,10.
0. ,0. ,25. ,40. ,5. ,1. , 1000.
O.,20.,2750.,165.,425.,4.5,38.,0.,1.5,1.5,0.,0.
0.5
1.0
2.0
3.0
5.0
10.0
20.0
20
20
20
20
20
20
,0.
,0.
,0.
,0.
,0.
,0.
20. ,0.
(F4.0,F3.0,F6.0,I2,3F5.0,F3.0)
270. 3. 1500. 4 .112 .175 290. .5
270. 3. 1500. 4 .112 .175 290. .5
A job stream
following form:
for a UNIVAC EXEC 8 system might have the
@RUN,R/R JOB-ID,ETC
@ASG,A MODELS*LOAD.
@ASG,A WINDS
©USE 21,WINDS
@ASG,R PLOT
@USE 22,PLOT
@XQT MODELS*LOAD.INPUFF
(input records shown above)
©FIN
Not needed for
verification run
The following is a
OS or MVS. Units
preallocated.
sample job stream
21 and 22 are
for an IBM system under
assumed to have been
//JOBID JOB
//XINPUFF EXEC
//STEPLIB DD
//FT21F001 DD
//FT22F001 DD
//FT06F001 DD
//FT05F001 DD
(input
/*
II
records
(PROJ,ACCT,OTHER),CLASS=A,TIME=1
PGM=INPUFF,TIME=(,30)
DSN=USER.MODELS.LOAD,DISP=SHR
DSN=USER.WINDS.DATA,DISP=SHR\ Not needed for
DSN=USER.PLOT.DATA,DISP=SHR / verification run
SYSOUT=A
*
shown above)
55
6-84
-------
A sample job stream for a CDC system under Scope 3.14 may
look as follows:
XX,T05,P4.
USER,HALE,EPA.
PROJECT,*PRJ*XX.
ATTACH,LIB,MODELSLIB,ID=XX.
ATTACH,TAPE21,WINDS,ID=XX.\ Not needed for
ATTACH, TAPE 2 2 , PLOT, I D=XX. J verification run
LIBARY,LIB.
INPUFF.
*
(input records shown above)
*
Figure 17 provides the output for the sample test. Users may
verify the proper execution of the program by comparing their
results with those given in the figure.
56 6-84
-------
INPUFF VERIFICATION RUN
tn
VERSION 84107
OPTIONS
A "1" INDICATES THAT THE OPTION
HAS BEEN EXERCISED
STACK DOWNWASH 1
UPDATE SOURCE DATA 0
USER SUPPLIED WIND FIELD 0
UNIT 22 OUTPUT OPTION 0
PRINT PUFF INFORMATION 0
INTERMEDIATE CONCENTRATIONS 0
INPUT PARAMETERS
ISIM= 0
ISTT1M= 5
STANUM= 5.00
ALPHA= 1.00
SYMAX= 1000.0
ANHGT= 10.0
I
00
Figure 17. Output for the sample test.
-------
«•• SOUKCL INFORMATION •••
SOURCE STACK STACK
STRENGTH HEIGHT TEMP.
(G/SEC) (M) (DEG-K)
.275E+04 165.00 425.000
*** METEOROLOGY
WIND Dili. WIND Sl'D. MIXING
(DEC) (M/SEC) (M)
270.0 3.000 1500
STACK GAS STACK VOLUME COORDINATES AT TIME 0 SOURCE SOURCE PLUME
VELOCITY DIAMETER FLOW EAST NORTH SPEED DIRECTION HEIGHT
(M/SEC) (M) (M**3/SEC) (KM) (KM) (M/SEC) (DEC) (M)
38.000 4.500 0.000 0.000 20.000 0.000 0.0 558.22
* » •
INITIAL SIGMAS
HGT. PROF.EP STABILITY U PLUME TEMP (R) (Z) SIGMA TH. SIGMA PH.
(DIMEN) (CLASS) (M/SEC) (K) (M) (RAD.) (RAD.)
0.150 4 4.702 290.0 1.5 1.5 0.1750 0.1120
SIMULATION PERIOD SIMULATION TIME PUFF RELEASE RATE SOURCE RECEPTOR DISTANCE DISPERSION
START (SEC) STOP (SEC) (SEC) (SEC) (KM) TYPE
0 3600 3600 15 O.SO 1
cn
00
3600 SEC AVG. CONCENTRATION AT RECEPTORS FOR SIMULATION PERIOD 0 TO 3600 SECONDS
RECEPTORS
X Y Z
0.500 20.000 0.000
1.000 20.000 0.000
2.000 20.000 0.000
3.000 20.000 0.000
5.000 20.000 0.000
10.000 20.000 0.000
20.000 20.000 0.000
Oi
1
00
4*.
CONCENTRATION (G/M»'3)
O.OOOE+00
O.OOOE+00
7.799E-20
8.515E-13
7.678E-08
1.872E-OS
2.919E-1-1
Figure 17. (continued)
-------
Cn
<0
••* SOURCE INFORMATION ••*
SOURCE STACK STACK
STRENGTH HEIGHT TEMP.
(G/SEC) (M) (DEG-K)
.275E+04 165.00 425.000
» « • METEOROLOGY
WIND DIR. WIND SPD. MIXING
(DEC) (M/SEC) (M)
270.0 3.000 1500
STACK GAS STACK VOLUME COORDINATES AT TIME 3600 SOURCE SOURCE PLUME
VELOCITY DIAMETER FLOW EAST NORTH SPEED DIRECTION HEIGHT
(M/SEC) (M) (M**3/SEC) (KM) (KM) (M/SEC) (DEC) (M)
38.000 4.500 0.000 0.000 20.000 0.000 0.0 558.22
• * *
INITIAL SIGMAS
IIGT. PROF.EP STABILITY U PLUME TEMP (R) (Z) SIGMA Til. SIGMA PH.
(DIMEN) (CLASS) (M/SEC) (K) (M) (RAD.) (RAD.)
0.150 4 4.702 290.0 1.5 1.5 0.1750 0.1120
SIMULATION PERIOD SIMULATION TIME PUFF RELEASE RATE SOURCE RECEPTOR DISTANCE DISPERSION
START (SEC) STOP (SEC)
3600 7200
(SEC) (SEC) (KM) TYPE
3600 15 0.50 1
3600 SEC AVG. CONCENTRATION AT RECEPTORS FOR SIMULATION PERIOD 3600 TO 7200 SECONDS
RECEPTORS
X Y Z
0.500 20.000 0.000
1.000 20.000 0.000
2.000 20.000 0.000
3.000 20.000 0.000
5.000 20.000 0.000
10.000 20.000 0.000
20.000 20.000 0.000
CONCENTRATION (G/M»«3)
O.OOOE+00
O.OOUE+00
1.268E-19
1.285E-12
1.227E-07
4.906E-05
1.360E-04
Figure 17. (continued)
oo
-------
2.00 HR AVG. CONCENTRATION AT RECEPTORS FOR ALL SIMULATION PERIODS
RECEPTORS
X Y Z CONCENTRATION (G/M»*3)
0.500
,000
000
3.000
5.000
10.000
20.000
20.000
20.000
20.000
20.000
20.000
20.000
20.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
O.OOOE+00
O.OOOE+00
1.024E-19
1.068E-12
9.974E-08
3.389E-05
6.801E-05
I
00
Figure 17. (continued)
-------
SECTION 11
INTERPRETATION OF OUTPUT
The output of INPUFF has ten parts, three of which are
optional. The output begins with printing the title of the run,
which can be up to 80 characters in length. The next printed
information is a list of options, followed by a list of the input
parameters. Next are the source data followed by a printout of
meteorological conditions used in the execution of the model for
the current simulation period. These are followed by five pieces
of information regarding how INPUFF simulates the release,
including: simulation period, simulation time, puff release rate,
minimum source-receptor distance, and dispersion type. The next
two output sections are optional. If IPIC = 1, then intermediate
concentrations are written every INC seconds. The time period for
which the averages are appropriate is printed in the first line
of the intermediate concentration output. If IPCC = 1, then
information on each puff is printed each ITIME in additi.on to
average concentrations at each receptor. Finally, a table of
average concentrations is output giving averages for each
receptor for all meteorological periods.
There is one other optional output available to the user. If
IP22 = 1, then information is written to unit 22, which can be
used later for plotting purposes.
The input stream and output listing of example problems 1 and
2 of Section 6 are presented in the next two sections. The
reader is referred to the earlier section for the physical
description of each problem. Intricacies of the input data are
discussed and the output listing is annotated for ease of
interpretat ion.
EXAMPLE 1 -- MOVING SOURCE
This example demonstrates an unique feature of INPUFF that
allows the source to move at a constant speed and direction over
a specified time (ITIME). If the source is changing speed and
direction more frequently than the meteorology, then ITIME can be
adjusted to match the frequency of source changes. The
meteorology is also changed accordingly. Table 11 lists the
input data; outputs of the example problem are given in Figure
18. Since IPCC = 1, the output includes puff information printed
for each ITIME.
61 6-84
-------
TABLE 11. INPUT DATA FOR EXAMPLE 1
Record Record Type
EXAMPLE 1 MOVING SOURCE 1
1,1,0,1,1,0 2
2,1200,-1,60,0,5,6,1,8,10. 5
0.,0.,25.,15.,10.,.75,1000. 6
O.,.2,600.,30.,390.,2.,15.,0.,l.,l.,90.,2. 7
1.54,1.19,0. 8
1.65,1.35,0. 8
2.,1.5,0. 8
2.35,1.35,0. 8
1.08,1.38,0. 8
1.3,1.7,0. 8
2.,2.,0 . 8
27170 8
(F4!o,F4.0,F6.0,I2,3F5.0,F3.0) 9
180. 3.5 3000. 3 .074 .105 290. .5 10
170. 4.0 3000. 4 .047 .067 288. .5 10
600.,30.,390.,2.,15.,0.,!.,!.,45.,2. 11
Note that the source information is updated every 20 minutes
for two periods. If, however, the source speed and direction
were changing every 5 minutes, NTIME could be changed to 8 and
ITIME to 300 seconds. There would now be eight record type 10's
and seven record type ll's. The first four meteorology records
would be the same, and the second four would also be the same.
In this example, the source speed and direction are assumed to
change at 20-minute intervals.
The information printed for each puff includes: puff number
and coordinates, time of puff release, total mass of the puff, ar
and oz for the puff, and its dispersion key. Because the puffs
combine as they travel downwind, the puff's characteristics are
adjusted each time it combines with another puff. For example,
puff 1 has a total mass of 60,000 grams. Since the source
strength is 600 g/sec and the puff release rate is 20 seconds,
this represents the combination of five puffs. All the
parameters are affected by puff combinations except the
dispersion key (KEYP). Puffs with different KEYP values do not
combine.
Plots of concentration versus time for each of the eight
receptors are shown in Figure 19. The coordinates of each
receptor are printed at the top of each plot. The input data
used in the execution of the plot programs are very short and are
shown below.
62 6-84
-------
Input Data Records Data
1 1
2t ~lyOy*yi)«jO«
3 1, 2, 3, 4, 5, 6, 7, 8
63 6-84
-------
o>
*>.
EXAMPLE 1
OPTIONS
MOVING SOURCE Run title
A "1" INDICATES THAT THE OPTION
HAS BEEN EXERCISED
VERSION 84107
STACK DOWNWASH 1
UPDATE SOURCE DATA 1
USER SUPPLIED WIND FIELD 0
UNIT 22 OUTPUT OPTION 1
PRINT PUFF INFORMATION 1
INTERMEDIATE CONCENTRATIONS 0
INPUT PARAMETERS
1SIM= 0
ISTT1M= 5
STANUM=10.00
ALPHA= 0.75
SYMAX= 1000.0
ANHGT= 10.0
Options and input parameters
exercised by the user
05
i
00
Figure 18. Annotated output of example 1
-------
Cn
1
oo
• • • SOURCE INFORMATION • • •
SOURCE STACK STACK STACK GAS STACK VOLUME COORD 1 NATI.b AT 1'IMI. U SOURCE SOURCE
STRENGTH HEIGHT TEMP. VELOCITY DIAMETER FLOW EAST NORTH SPEED DIRECTION
(G/SEC) (M) (DEG-K) (M/SEC) (M) (M»*3/SEC) (KM) (KM) (M/SEC) (DEG)
.600E+03 30.00 390.000 15.000 2.000 0.000 0.000 0.200 2.000 90.0
•** METEOROLOGY «••
INITIAL S1GMAS
WIND DIR. WIND SPD. MIXING HGT. PROF . EP STABILITY U PLUME TEMP (It) (Z) SIGMA TH . SIGMA PH.
(DEG) (M/SEC) (M) (DIMEN) (CLASS) (M/SEC) (K) (M) (RAD.) (RAD.)
180.0 3.500 3000. 0.100 3 4.462. 290.0 1.0 1.0 0.1050 0.0740
PLUME
HEIGHT
(M)
113.47
SIMULATION PERIOD SIMULATION TIME PUFF RELEASE RATE SOURCE RECEPTOR DISTANCE DISPERSION \ D ff • i ,.• „
START (SEC) STOP (SEC) (SEC) (SEC) (KM) TYPE \ ^UJJ Simulation
o 1200 i2oo-«— ITIME 20+—ISTEP o-so« — CDIS i / information
re lease puff |
puff location time mass ag P-G dispersion
(m) (sec) (gj scheme
PUFF* X Y Z TIME TOTAL Q SY SZ KEYP
1 60.002 5420.866 113.471 30 48000.00 459.583 277.755 1
2 210.002 5086.195 113.471 105 36000.00 433.114 261.459 1 Since IPCC = 1 the location
3 340.002 4796.146 113.471 170 48000.00 410.007 247.259 1 , J J,
• release time, mass, size, and
* • dispersion fay are printed for
31 2319.999 378.491 113.471 1160 uooo.oo 21.965 13.302 i each puff within the modeling
32 2359.999 289.246 113.471 1180 12000.00 11.953 7.468 1 -r>f>njr>n
33 2399.999 200.000 113.471 1200 12000.00 1.000 1.000 1 J-K'dl'u"-
ITIME + 1200 SEC AVG. CONCENTRATION AT RECEPTORS FOR SIMULATION PERIOD 0 TO 1200 SECONDS
RECEPTORS
X Y Z CONCENTRATION (G/M»»3)
1.540 .190 0.000 1.353E-04
1.650 .350 o.ooo 1.069E-04 Average concentrations at each
2.000 .500 0.000 2.168E-05
I'OBO 'ago o'ooo }'«o8E-o! feceptor are printed at the end
1.300 .700 0.000 2.206E-04 ,, , , , . , • ,
2.000 2.000 o.ooo 3.675E-06 °J &acn meteoroiogicai period.
2.700 1.700 0.000 2.921E-14
Figure 18. (continued)
-------
O5
O5
CD
I
00
Next meteorological period. Source parameters and
meteorology are different from the previous period.
s o u it c E
N FORM AT ION
SOURCE STACK
STRENGTH HEIGHT
(G/SEC) (M)
.600E+03 30.00
STACK
TliMP.
(UEG-K)
390.000
STACK GAS STACK VOLUME
VELOCITY DIAMETER FLOW
(M/SEC) (M) (M*«J/SEC)
COORDINATES AT TIMI. 12(10
EAST NOKI'll
(KM) (KM)
15.000
2.000
0.000
2.400
0. 200
SOURCE
SPEED
(M/SEC)
2. 000
SOURCE
DIRECTION
(DEG)
45.0
PLUML
HI.1GIIT
(M)
100.17
METEOROLOGY
WIND DIH.
(DEG)
170.0
WIND SPD. MIXING IIGT.
(M/SEC) (M)
PROF.EP
(DIMEN)
STABILITY U PLUME TEMP
(CLASS) (M/SEC) (K)
INITIAL SIGMAS
(R) (/,) SIGMA Til. SIGMA PH.
(M) (RAD. ) (RAD. )
4.000
3000.
0. 150
5.652
288 .0
1 .0
1 .0
0.0670
0.0470
SIMULATION PERIOD SIMULATION TIME
START (SEC) STOP (SEC) (SEC)
1200 2400 1200
PUFF RELEASE RATE
(SLC)
12
SOURCE RECEPTOR DISTANCE
(KM)
0. 50
DISPERSION
TYPE
1
PUFF*
1
2
3
45
46
47
7.340 9589.681
207.339 9143.456
407.339 8697.231
4039.568 1996.690
4068.315 1946.872
4097.062 1897.054
113.471
113.471
113.471
TIME
592
692
792
TOTAL Q
SY
KEYP
100.167 2376
100.167 2388
100.167 2400
60000.00
60000.00
60000.00
7200.00
7200.00
7200.00
595.576
563.776
531 . 498
11 . 662
6.552
1 .000
273.751
253.405
233. 123
6.484
3.967
1 .000
1
1
1
•
1
1
1200 SEC AVG. CONCENTRATION AT RECEPTORS FOH SIMULATION PERIOD 1200 TO 2400 SECONDS
RECEPTORS
X
1.540
1.650
2.000
2.350
1.080
1.300
2.000
2.700
Y
1. 190
1.350
1.500
1.350
1.380
1.700
2.000
1.700
Z
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
CONCENTRATION
1. 154E-07
7 .611E-06
1.233E-04
1.310E-05
4.781E-10
2.136E-06
1.803E-04
1 .285E-05
(G/M»*3)
Average concentrations for the
second meteorological period.
Figure 18. (continued)
-------
o>
-3
Length of simulation time
0.67 HR AVG. CONCENTRATION AT RECEPTORS FOR ALL SIMULATION PERIODS
X
1.540
1.650
2.000
2.350
1.080
1.300
2.000
2.700
RECEPTORS
Y
1.190
1.350
1.500
1.350
' 1.380
1.700
2.000
1.700
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
CONCENTRATION (G/M**3)
6.770E-05
8.725E-05
7.247E-05
6.555E-06
9.032E-05
1.114E-04
9.198E-05
Average concentrations at each
receptor over the modeling period.
6.423E-06
I
00
Figure 18. (continued)
-------
I.14. t.ltl
It If M 31 *• 4f M
It-'-
-i.-1
-I.'1
Tint 111 mmrrtl
o* i I.M. i.Jli
• • it tt
X «t •» M
Tine I" Htmt
If1-
-It-'
•ICfPVO* t I.M. t.311
TIM IH NlMirU
MCCTM < a.««. *.*ti
f-
-I."'
• I !•
t I It
-It*'
TUC III HtNWTCf
Figure 19. Concentration versus time plots for example 1
68
6-84
-------
EXAMPLE 2 -- LOW LEVEL SOURCE WITH LOW WIND SPEED CONDITIONS
This problem illustrates the model simulation for a low level
release during conditions of light and variable winds. The input
data stream is shown in Table 12 and the abridged output in
Figure 20. A very important difference between this example and
the previous example is that for this example DISKEY on record 5
has been assigned a value of 2. Dispersion downwind of the
source is no longer characterized by travel distance but by
travel time using the on-site dispersion scheme. The values
assigned to aa and ae are not used in the P-G characterization of
dispersion. However, in the on-site scheme, CTV and az are
functions of oa and ae. If measures of these fluctuation
statistics are available, they should be used. If none are
available, then typical values can be used based on the state of
the atmosphere.
In this example, twelve simulation periods of 10-minute
duration are used to simulate the 2-hour release. The atmosphere
is stable with large fluctuations in the wind direction. aa has
been assigned a large value, typical of low wind speed
conditions. The strength of the source is decaying with time;
initially the source strength is 825 g/sec, but by the 12th
period it has dropped to 12 g/sec.
In this simulation, average concentrations at each receptor
are printed every 10 minutes. The puff locations at the end of
each 10 minute period are plotted in Figure 21. The input data
for the plot program are shown below:
Input Data Record Data
1 2
2 O.,0.,5.,5.,5.,5.
Circles drawn around the centers of the puff positions have a
radius equal to ar.
69 6-84
-------
TABLE 12. INPUT DATA FOR EXAMPLE 2
Record
EXAMPLE 2 LOW LEVEL
0,1,0,1,0,0
12, 600, -1,300, 0,5, 6, 2, 10,
0. ,0. ,25. ,15. ,10. ,1. ,1000
2. ,1. ,825. ,3. ,290. , .5,10.
1.54,1.19,0.
1.65,1.35,0.
2. ,1.5,0.
2.35,1.35,0.
2.46,1.19,0.
1.08,1.38,0.
1.3,1.7,0.
2. ,2. , 0.
2.7,1.7,0.
2.92,1.38,0.
(F4
180
210
175
145
155
210
200
182
170
195
185
195
562
383
261
178
121
83.
56.
38.
26.
18.
12.
.0
•
•
•
•
•
•
•
•
•
•
•
•
• »
• >
• »
• >
• >
,3
,3
,3
,3
,3
,3
,F3.0,F6.0
•
•
•
•
•
•
•
•
•
•
•
•
3
3
3
3
3
•
•
•
•
•
•
5
5
5
5
5
5
5
5
5
5
5
5
.
.
•
.
•
>
>
>
>
»
>
5000
5000
5000
5000
5000
5000
5000
5000
5000
5000
5000
5000
,290.
,290.
,290.
,290.
,290.
290. ,
290. ,
290.,
290. ,
290. ,
290. ,
. 6
. 6
. 6
. 6
. 6
. 6
. 6
. 6
. 6
. 6
. 6
. 6
,.5
,.5
,.5
,.5
,.5
.5,
.5,
.5,
.5,
.5,
.5,
,1
•
•
•
•
•
•
•
•
•
•
•
•
,1
Record type
SOURCE LOW WIND SPEED 1
2
10. 5
6
,0.,1.,1.,0.,0. 7
8
8
8
8
8
8
8
8
8
8
2,3F5.0,F3.0)
035
035
035
035
035
035
035
035
035
035
035
035
0. ,0
,10. ,0
,10. ,0
,10. ,0
,10. ,0
10
10
10
10
10
10
. ,0.
. ,0.
. ,0.
. ,0.
. ,0.
. ,0.
.393
.393
.393
.393
.393
.393
.393
.393
.393
.393
.393
.393
.,1.
.,1.
.,1.
.,1.
.,1.
,1.,
,1-,
,1.,
,1.,
,1.,
,1.,
290.
290.
290.
290.
290.
290.
290.
290.
290.
290.
290.
290.
,1. ,0
,1. ,0
,1. ,0
,1. ,0
,1. ,0
1. ,0.
1 . , 0 .
1 . , 0 .
1. ,0.
1. ,0.
1. ,0.
.5
.5
.5
.5
.5
.5
.5
.5
.5
.5
.5
.5
.,0.
. ,0.
. ,0.
. ,0.,
. ,0.
,0.
,0.
,0.
,0.
,0.
,0.
9
10
10
10
10
10
10
10
10
10
10
10
10
11
11
11
11
11
11
11
11
11
11
11
70
6-84
-------
05
I
00
EXAMPLE 2
OPTIONS
LOW LEVEL SOURCE LOW WIND SPEED
A "1" INDICATES THAT THE OPTION
HAS BEEN EXERCISED
VERSION 84107
STACK DOWNWASH 0
UPDATE SOURCE DATA 1
USER SUPPLIED WIND FIELD 0
UNIT 22 OUTPUT OPTION 1
PRINT PUFF INFORMATION 0
INTERMEDIATE CONCENTRATIONS 0
INPUT PARAMETERS
ISIM= 0
ISTTIM= 5
STANUM=10.00
ALPHA= 1.00
SYMAX= 1000.0
ANHGT= 10.0
Figure 20. Annotated output of example 2.
-------
SOURCE INFORMATION
SOURCE STACK
STRENGTH HEIGHT
(G/SEC) (M)
.825E+03 3.00
STACK
TEMP.
(DEG-K)
290.000
STACK GAS STACK VOLUME COORDINATES AT TIME
VELOCITY DIAMETER FLOW EAST NORTH
(M/SEC) (M) (M*»3/SEC) (KM) (KM)
10.000 0.500
0.000
2.000
1.000
SOURCE
SPEED
(M/SEC)
0.000
SOURCE
DIRECTION
(DEG)
0.0
PLUME
HEIGHT
(M)
11.50
to
METEOROLOGY
WIND DIR.
(DEG)
180.0
WIND SPD. MIXING HGT.
(M/SEC) (M)
PROF.EP
(DIMEN)
0.500
5000.
0.350
INITIAL SIQV1AS
STABILITY U PLUME TEMP (R) (Z) SIGMA TH. SIGMA PH.
(CLASS) (M/SEC) (K) (M) (RAD.) (RAD.)
6
0.525
290.0
1.0
1.0
SIMULATION PERIOD SIMULATION TIME
START (SEC) STOP (SEC) (SEC)
0 600 600
PUFF RELEASE RATE
(SEC)
100
SOURCE RECEPTOR DISTANCE
(KM)
0.50
0.3930
0.0350
DISPERSION
TYPE
2
600 SEC AVG. CONCENTRATION AT RECEPTORS FOR SIMULATION PERIOD
0 TO
600 SECONDS
1.540
1.650
2.000
2.350
2.460
1.080
1.300
2.000
2.700
2.920
RECEPTORS
Y
1.190
1.350
1.500
1.350
1.190
1.380
1.700
2.000
1.700
1.380
CONCENTRATION (G/M*»3)
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
6.153E-15
1.533E-10
9.132E-07
1.533E-10
6.153E-15
O.OOOE+00
O.OOOE+00
O.OOOE+00
O.OOOE+00
O.OOOE+00
O>
I
oo
*>.
Figure 20. (continued)
-------
OURCE INFORMATION
05
I
oo
SOURCE STACK STACK
STRENGTH HEIGHT TEMP.
(G/SEC) (M) (DEG-K)
.562E+03
3.00
290.000
STACK GAS STACK VOLUME COORDINATES AT TIME 600 SOURCE SOURCE
VELOCITY DIAMETER FLOW EAST NORTH SPEED DIRECTION
(M/SEC) (M) (M»*3/SEC) (KM) (KM) (M/SEC) (DEC)
10.000
0.500
0.000
2.000
1.000
0.000
0.0
PLUME
HEIGHT
(M)
11.50
Co
METEOROLOGY
WIND DIR.
(DEG)
210.0
WIND SPD. MIXING HGT.
(M/SEC) (M)
PROF . EP
(DIMEN)
0.500
5000.
0.350
STABILITY U PLUME TEMP
(CLASS) (M/SEC) (K)
INITIAL SIGMAS
(R) (Z) SIGMA TH. SIGMA PH.
(M) (RAD.) (RAD.)
6
0.525
290.0
1.0
1 .0
SIMULATION PERIOD SIMULATION TIME
START (SEC) STOP (SEC) (SEC)
600 1200 600
PUFF RELEASE RATE
(SEC)
100
SOURCE RECEPTOR DISTANCE
(KM)
0.50
0.3930
0.0350
DISPERSION
TYPE
2
600 SEC AVG. CONCENTRATION AT RECEPTORS FOR SIMULATION PERIOD
600 TO
1200 SECONDS
RECEPTORS
Y
1.540
1.650
2.000
2.350
2.460
1.080
1.300
2.000
2.700
2.920
1.190
1.350
1.500
1.350
1.190
1.380
1.700
2.000
1.700
1.380
I CONCENTRATION (G/M«*3)
0.000 1.S75E-10
0.000 2.521E-07
0.000 3.030E-03
0.000 5.136E-04
0.000 3.8S9E-06
0.000 7.277E-21
0.000 1.028E-14
0.000 4.525E-07
0.000 3.338E-08
0.000 2.713E-12
Figure 20. (continued)
-------
£
0>
1
oo
••* SOURCE INFORMATION •
SOURCE STACK STACK STACK GAS STACK
STRENGTH HEIGHT TEMP. VELOCITY DIAMETER
(G/SEC) (M) (DEG-K) (M/SEC) (M)
.383E+03 3.00 290.000 10.000 0.500
*•* METEOROLOGY • * •
WIND D1R. WIND SPD. MIXING HGT. PROF.EP
(DEC) (M/SEC) (M) (DIMEN)
175.0 0.500 5000. 0.350
SIMULATION PERIOD SIMULATION TIME PUFF
START (SEC) STOP (SEC) (SEC)
1200 1800 600
600 SEC AVG. CONCENTRATION AT RECEPTORS
RECEPTORS
X Y Z CONCENTRATION
1.540 .190 0.000 2.960E-08
1.650 .350 0.000 8.590E-06
2.000 .500 0.000 1.287E-02
2.350 .350 0.000 1.617E-03
2.460 .190 0.000 1.456E-05
1.080 .380 0.000 5.749E-13
1.300 .700 0.000 2.048E-08
2.000 .000 0.000 3.338E-03
2.700 .700 0.000 2.508E-05
2.920 .380 0.000 4.712E-09
* *
VOLUME COORDINATES AT TIME 1200 SOURCE SOURCE PLUME
FLOW EAST NORTH SPEED DIRECTION HEIGHT
(M**3/SEC) (KM) (KM) (M/SEC) (DEG) (M)
0.000 2.000 1.000 0.000 0.0 11.50
INITIAL SIGMAS
STABILITY U PLUME TEMP (R) (Z) SIGMA TH. SIGMA PH.
(CLASS) (M/SEC) (K) (M) (RAD.) (RAD.)
6 0.525 290.0 1.0 1.0 0.3930 0.0350
RELEASE RATE SOURCE RECEPTOR DISTANCE DISPERSION
(SEC) (KM) TYPE
100 0.50 2
FOR SIMULATION PERIOD 1200 TO 1800 SECONDS
(G/M**3) Output is abridged. The following
meteorological periods are missing
from the sample output:
1800 to 2400 sec,
2400 to 3000 sec,
3000 to 3600 sec,
3600 to 4200 sec,
4200 to 4800 sec,
4800 to 5400 sec,
5400 to 6000 sec, and
6000 to 6600 sec.
Figure 20. (continued)
-------
SOURCE INFORMATION
OO
SOURCE STACK STACK
STRENGTH HEIGHT TEMP.
(G/SEC) (M) (DEG-K)
STACK GAS STACK VOLUME
VELOCITY DIAMETER FLOW
(M/SEC) (M) (M**3/SEC)
.120E+02 3.00 290.000 10.000 0.500 0.000
Source strength has decayed to 12 g/sec
from an original value of 825 g/sec.
*** METEOROLOGY »•*
COORDINATES AT TIME 6600
EAST NORTH
(KM) (KM)
2.000
1.000
SOURCE
SPEED
(M/SEC)
0.000
SOURCE
DIRECTION
(DEG)
0.0
PLUME
HEIGHT
(M)
11.50
WIND DIR.
(DEG)
195.0
WIND SPD. MIXING HOT.
(M/SEC) (M)
PROF.EP
(DIMEN)
0.500
5000.
0.350
STABILITY U PLUME TEMP
(CLASS) (M/SEC) (K)
INITIAL SIGMAS
(R) (Z) SIGMA TH.
(RAD.)
0.525
290.0
(M)
1.0 1.0
SIMULATION PERIOD SIMULATION TIME
START (SEC) STOP (SEC) (SEC)
6600 7200 600
PUFF RELEASE RATE
(SEC)
100
SOURCE RECEPTOR DISTANCE
(KM)
0.50
0.3930
SIGMA PH.
(RAD.)
0.0350
DISPERSION
TYPE
2
600 SEC AVG. CONCENTRATION AT RECEPTORS FOR SIMULATION PERIOD 6600 TO 7200 SECONDS
RECEPTORS
Y
1.540
1.650
2.000
2.350
2.460
1.080
1.300
2.000
2.700
.190
.350
.500
.350
.190
.380
.700
.000
.700
2.920
1.380
Z CONCENTRATION (G/M*«3)
0.000 4.081E-09
0.000 6.713E-07
0.000 5.585E-04
0.000 2.517E-05
0.000 2.147E-07
0.000 1.294E-10
0.000 3.769E-07
0.000 3.437E-03
0.000 1.465E-05
0.000 1.026E-08
Figure 20. (continued)
-------
Oi
1
oo
2.00 HR AVG. CONCENTRATION
RECEPTORS
X Y Z
1.540 1.190 0.000
1.650 1.350 0.000
2.000 1.500 0.000
2.350 1.350 0.000
2.460 1.190 0.000
1.080 1.380 0.000
1.300 1.700 0.000
2.000 2.000 0.000
2.700 1.700 0.000
2.920 1.380 0.000
AT RECEPTORS FOR ALL SIMULATION PERIODS
CONCENTRATION (G/M**3)
1.346E-05
5.363E-04
4.846E-03
3.234E-04
3.827E-06
1.945E-07
1.803E-04
1.630E-02
4.375E-05
3.052E-08
Figure 20. (continued)
-------
A. TRAJECTORV OF PUFF(SI
B. TRAJECTORV OF PUFF(S)
Q TRAJECTORY OF PUFF(S)
4
•1
2
*- i —
8 —
1
.7
•6 .^
c
R
.9
•
9
1
1 •
b —
•
-
.
; i
1
. 7
.6 .2
E
R
, 9
P'.\ '
•
e
•
!
5-
1 :
g •
S t
0 3
z
3
O
vt
X
1 S
C
(
I
.7
•6 .1B<
ITT-I-
Z
-a
^ .9
Till
e
-rr-T-r-
•
s
x, EAST-UEST COORDINATE x(k«J, EAST-UEST COORDINATE x(k«), EAST-UEST COORDINATE
Q TRAJECTORV OF PUFF(S) £ TRAJECTORY OF PUFF (5 ) P TRAJECTORV OF PUFF! S )
5— 1 .
•1 -
-
z—
-
*" 1 —
e-
e
t
r
.7 J^
16 ft
* ^
a
}
/ .9
3
, S
3
0
4
5
1 :
0
X
1 :
r -
X
9
1
I
ff
•fl
.6 m
. 1^
z
R
.9
3.< J
.5
:
e
•
S
(
1
pA.
.7 %!
.6 &
.1
a
k .9
1 * '
r >s
•
0
.
5
0» x(k«i>, EAST-UEST COORDINATE x, EAST-UEST COORDINATE x(kn). EAST-UEST COORDINATE
00
Figure 21. Puff locations at the end of each simulation period.
A - L represents 10-minute intervals
-------
Ill I I I 1
JIBS 1011003 Mjnos-HivoN -IUHU
3JBX10HOCD
78
6-84
-------
REFERENCES
Batchelor, O. G. 1952. The Theory of Homogeneous Turbulence.
Cambridge University Press, London.
Briggs, G. A. 1969. Plume Rise. USAEC Critical Review Series.
TID-25075, National Technical Information Service,
Springfield, VA. 81 pp.
Briggs, G. A. 1971. Some Recent Analyses of Plume Rise
Observation. In: Proceedings of the Second International
Clean Air Congress, H. M. Englund and W. T. Beery, eds.,
Academic Press, New York. pp. 1029-1032.
Briggs, G. A. 1973. Diffusion Estimation for Small Emissions.
NOAA Atmospheric Turbulence and Diffusion Laboratory,
Contribution File No. (Draft) 79. Oak Ridge, TN. 59 pp.
Briggs, G. A. 1975. Plume Rise Predictions. In: Lectures on
Air Pollution and Environmental Impact Analysis, D. A.
Haugen, ed. Am. Meteorol. Soc., Boston, MA. pp. 59-111.
Cramer, H. E. 1976. Improved Techniques for Modeling the
Dispersion of Tall Stack Plume. In: Proceedings of the
Seventh International Technical Meeting on Air Pollution
Modeling and its Application. No. 51, NATO/CCMS, pp. 731-780
(NTIS PB 270 799).
Draxler, R. R. 1976. Determination of Atmospheric Diffusion
Parameters. Atmos. Environ., 10: 99-105.
Hanna, S. R., G. A. Briggs, J. Deardorff, B. A. Egan, F. A.
Gifford, and F. Pasquill. 1977. AMS-Workshop on Stability
Classification Schemes and Sigma Curves—Summary of
Recommendations. Bull. Am. Meteorol. Soc., 58: 1305-1309.
Irwin, J. S. 1983. Estimating Plume Dispersion - A Comparison
of Several Sigma Schemes. J. Climate Applied Meteorol.,
22: 92-114.
Pasquill, F. 1961. The Estimation of the Dispersion of
Windborne Material. Meteorol. Magazine, 90: 33-49.
79 6-84
-------
Pasquill, F. 1976. Atmospheric Dispersion Parameters in
Gaussian Plume Modeling. Part II. Possible Requirements for
Change in the Turner Workbook Values. EPA-600/4-76-030b,
U.S. Environmental Protection Agency, Research Triangle
Park, NC. 44 pp.
Taylor, G. I. 1921.
Proceedings of
20: 196.
Diffusion by Continuous Movements. In:
the London Mathematical Society, Series 2,
Turner, D. B. 1970.
Estimates. Office
(NTIS PB 191 482).
Research Triangle Park,
Workbook of Atmospheric Dispersion
of Air Programs Publication No. AP-26
U.S. Environmental Protection Agency,
NC. 84 pp.
80
6-84
-------
APPENDIX A
PLUME RISE
The use of the methods of Briggs to estimate plume rise and
effective height of emission are discussed below. In all
calculations, it is assumed that actual or estimated wind speed
at stack top, u(h), is available.
STACK DOWNWASH
To consider stack downwash, the physical stack height is
modified following Briggs (1973, p. 4). The h' is found from
h' = h + 2{[vs/u(h)] - 1.5}d for vs < 1.5u(h), (A-l)
h' = h for vs > 1.5u(h),
where h is physical stack height (meters), vs is stack gas
velocity (meters per second), and d is inside stack-top diameter
(meters). The h' is used throughout the plume height
computation. If stack downwash is not considered, h' = h in the
equat i ons.
BUOYANCY FLUX
For most plume rise calculations, the value of the Briggs
buoyancy flux parameter, F (mVs3), is needed. The following
equation is equivalent to Briggs' Eq. 12 (1975, p. 63):
F = (gvsd2AT)/(4Ts), (A-2)
where AT = Ts - T, Ts is stack gas temperature (degrees kelvin),
and T is ambient air temperature (degrees kelvin).
UNSTABLE OR NEUTRAL: CROSSOVER BETWEEN MOMENTUM AND BUOYANCY
For cases with stack gas temperature greater than or equal to
ambient air temperature, it must be determined whether the plume
rise is dominated by momentum or buoyancy. The crossover
temperature difference (AT)C is determined for (1) F less than 55
and (2) F greater than or equal to 55. If the difference between
stack gas temperature and ambient air temperature, AT, exceeds or
equals the (AT)C, plume rise is assumed to be buoyancy dominated;
if the difference is less than (AT)C, plume rise is assumed to be
81 6-84
-------
momentum dominated (see below).
The crossover temperature difference is found by setting
Briggs' Eq. 5.2 (1969, p. 59) equal to the combination of Briggs1
Eqs. 6 and 7 (1971, p. 1031) and solving for AT. For F less than
55,
(AT)C = 0.0297v^/3Ts/d2'3.
For F equal to or greater than 55,
(AT)C = 0.00575v|/3Ts/d1/3
UNSTABLE OR NEUTRAL: BUOYANCY RISE
(A-3)
(A-4)
For situations where AT exceeds or is equal to (AT)C as
determined above, buoyancy is assumed to dominate. The distance
to final rise xf (in kilometers) is determined from the
equivalent of Briggs1 Eq. 7 (1971, p. 1031), and the distance to
final rise is assumed to be 3.5x*, where x* is the distance at
which atmospheric turbulence begins to dominate entrainment. For
F less than 55,
Xf = 0.049F5'8.
For F equal to or greater than 55,
xf = 0.119F2'5.
(A-5)
(A-6)
The plume height, H (in meters), is determined from the
equivalent of the combination of Briggs' Eqs. 6 and 7 (1971,
p. 1031). For F less than 55,
H = h' + 21.425F3/lf/u(h),
and for F equal to or greater than 55,
H = h' + 38.71F3/5/u(h).
UNSTABLE OR NEUTRAL: MOMENTUM RISE
(A-7)
(A-8)
For situations where the stack gas temperature is less than
the ambient air temperature, it is assumed that the plume rise is
dominated by momentum. Also, if AT is less than (AT)C from Eq.
A-3 or A-4, it is assumed that the plume rise is dominated by
momentum. The plume height is calculated from Briggs' Eq. 5.2
(1969, p. 59):
H = h' + 3dvs/u(h).
(A-9)
Briggs (1969) suggests that this equation is most applicable when
vs/u is greater than 4. Since momentum rise occurs quite close
82
6-84
-------
to the point of release, the distance to final rise is set equal
to zero.
STABILITY PARAMETER
For stable situations, the stability parameter s is
calculated from the following equation (Briggs, 1971, p. 1031):
s = g(36/3z)/T. (A-10)
As an approximation, for stability class E, 36/3z is taken as
0.02 K/m, and for stability class F, 36/3z is taken as 0.035 K/m.
STABLE: CROSSOVER BETWEEN MOMENTUM AND BUOYANCY
For cases with stack gas temperature greater than or equal to
ambient air temperature, it must be determined whether the plume
rise is dominated by momentum or buoyancy. The crossover
temperature difference (AT)C is found by setting Briggs' Eq. 59
(1975, p. 96) equal to Briggs1 Eq . 4.28 (1969, p. 59), and
solving for AT. The result is
(AT)C = 0.019582vsT s 1/2 . (A-ll)
If the difference between stack gas temperature and ambient air
temperature (AT) exceeds or equals (AT)C, the plume rise is
assumed to be buoyancy dominated; if AT is less than (AT)C, the
plume rise is assumed to be momentum dominated.
STABLE: BUOYANCY RISE
For situations where AT is greater than or equal to (AT)C,
buoyancy is assumed to dominate. The distance to final rise (in
kilometers) is determined by the equivalent of a combination of
Briggs' Eqs. 48 and 59 (1975, p. 96):
xf = 0.0020715u(h)s
-1/2
The plume height is determined by the equivalent of Briggs'
Eq. 59 (1975, p. 96):
H = h' + 2.6{F/[u(h)s]}1/3 . (A-13)
The stable buoyancy rise for calm conditions (Briggs, 1975,
pp. 81-82) is also evaluated:
H = h' + 4F1/*s-3'8. (A-14)
The lower of the two values obtained from Eqs. A-13 and A-14 is
taken as the final effective height.
83 6-84
-------
By setting Eqs. A-13 and A-14 equal to each other and solving
for u(h), one can determine the wind speed that yields the same
plume rise for the wind conditions (A-13) as does the equation
for calm conditions (A-14). This wind speed is
u(h) = (2.6/4)3F1/'is1/8
= 0.2746F1/*s1/8. (A-15)
For wind speed less than or equal to this value, Eq. A-14
should be used for plume rise; for wind speeds greater than this
value, Eq. A-13 should be used.
STABLE: MOMENTUM RISE
When the stack gas temperature is less than the ambient air
temperature, it is assumed that the plume rise is dominated by
momentum. If AT is less than (AT)C as determined by Eq. A-ll, it
is also assumed that the plume rise is dominated by momentum.
The plume height is calculated from Briggs' Eq. 4.28 (1969,
p. 59):
H = h' + 1.5{(v|d2T)/[4Tsu(h)]}1/3s-1/6. (A-16)
The equation for unstable or neutral momentum rise (A-9) is
also evaluated. The lower result of these two equations is used
as the resulting plume height.
REFERENCES
Briggs, G. A. 1969. Plume Rise. USAEC Critical Review Series.
TID-25075, National Technical Information Service,
Springfield, VA. 81 pp.
Briggs, G. A. 1971. Some Recent Analyses of Plume Rise
Observation. In: Proceedings of the Second International
Clean Air Congress, H. M. Englund and W. T. Beery, eds.,
Academic Press, New York. pp. 1029-1032.
Briggs, G. A. 1973. Diffusion Estimation for Small Emissions.
NOAA Atmos. Turb. and Diff. Lab., Contribution File No.
(Draft) 79. Oak Ridge, TN. 59 pp.
Briggs, G. A. 1975. Plume Rise Predictions. In: Lectures on
Air Pollution and Environmental Impact Analysis, D. A.
Haugen, ed., Am. Meteorol. Soc., Boston, MA. pp. 59-111.
84 6-84
-------
APPENDIX B
LISTING OF FORTRAN SOURCE CODE
INPUFF
The source code listing of INPUFF follows. The program
consists of a main module, 8 subroutines, and 2 functions.
85 6-84
-------
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42. .
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.
57.
58.
59.
60.
61.
62.
63.
64.
65.
66.
67.
68.
69.
70.
71.
72.
73.
74.
75.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
INPUFF VERSION 84107
PROGRAM ABSTRACT: INPUFF
INPUFF IS A GAUSSIAN INTEGRATED PUFF MODEL. THE
GAUSSIAN PUFF DIFFUSION EQUATION IS USED TO COMPUTE
THE CONTRIBUTION TO THE CONCENTRATION AT EACH RECEPTOR
FROM EACH PUFF EVERY TIME STEP. COMPUTATIONS IN INPUFF
CAN BE MADE FOR A SINGLE POINT SOURCE AT UP TO 25
RECEPTOR LOCATIONS. IN PRACTICE HOWEVER, THE NUMBER OF
RECEPTORS SHOULD BE KEPT TO A MINIMUM IN THE DEFAULT
MODE THE MODEL ASSUMES A HOMOGENEOUS WIND FIELD.
HOWEVER, THE USER HAS THE OPTION OF SPECIFYING THE
WIND FIELD FOR EACH METEOROLOGICAL PERIOD AT UP TO
100 USER DEFINED GRID LOCATIONS. THREE DISPERSION
ALGORITHMS ARE UTILIZED WITHIN INPUFF FOR DISPERSION
DOWNWIND OF THE SOURCE. THE FIRST TWO ARE PASQUILL'S
SCHEME AS DISCUSSED BY TURNER(1970) AND A DISPERSION
ALGORITHM DISCUSSED BY IRWIN(1983). THIS LATTER SCHEME
IS A SYNTHESIS OF DRAXLER'S(1976) AND CRAMER1S(1976)
IDEAS. THE THIRD DISPERSION SCHEME IS FOR LONG TRAVEL
TIMES IN WHICH THE GROWTH OF THE PUFF BECOMES PROPORTIONAL
TO THE SQUARE ROOT OF TIME. A SOFTWARE PLOTTING PACKAGE
IS PROVIDED TO DISPLAY CONCENTRATION VERSUS TIME FOR A
GIVEN RECEPTOR AND THE PUFF TRAJECTORIES AFTER EACH
SIMULATION TIME.
REFERANCES:
TURNER, D. BRUCE, 1970: WORKBOOK OF ATMOSPHERIC DISPERSION
ESTIMATES, OFFICE OF AIR PROGRAMS PUBLICATION NO. AP-26.
U. S. ENVIRONMENTAL PROTECTION AGENCY, RESEARCH TRIANGLE
PARK, NC. 84 PP. [NTIS PB 191 482].
IRWIN, J. S., 1983: ESTIMATING PLUME DISPERSION—A COMPARISON
OF SEVERAL SIGMA SCHEMES. JOURNAL OF CLIMATE AND APPLIED
METEOROLOGY, 22, 92-114.
DRAXLER, R. R., 1976: DETERMINATION OF ATMOSPHERIC
DIFFUSION PARAMETERS, ATMOS. ENVIRON., 10, 99-105.
CRAMER, H. E., 1976: IMPROVED TECHNIQUES FOR MODELING THE
DISPERSION OF TALL STACK PLUME. PROCEEDINGS OF THE 7TH
INTERNATIONAL TECHNICAL MEETING ON AIR POLLUTION MODELING
AND ITS APPLICATION. N. 51, NATO/CCMS, 731-780.
[NTIS PB 270 799]
AUTHOR OF MODEL CODE FOR INPUFF:
W. B. PETERSEN
ON ASSIGNMENT TO THE ENVIRONMENTAL PROTECTION AGENCY FROM
THE NATIONAL OCEANIC AND ATMOSPHERIC ADMINSTRATION,
DEPT. OF COMMERCE.
PROGRAM SUPPORTED BY:
ENVIRONMENTAL OPERATIONS BRANCH
MAIL DROP 80, EPA.
RTP, NC. 27711
PHONE (919) 541-4564, FTS 629-4564
EXTERNAL FILES:
THERE ARE TWO OPTIONAL EXTERNAL FILES WHICH MAY BE
UTILIZED IN THE EXECUTION OF INPUFF. IF THE USER SUPPLIES
HIS OWN GRIDDED WIND FIELD DATA, THE WIND SPEED AND
DIRECTION ARE READ BY INPUFF FROM UNIT 21. THE FORMAT FOR
UNIT 21 IS USER DEFINED. UNIT 22 IS AN OUTPUT FILE GENERATED
BY INPUFF IF IP22=1. THIS FILE IS USED BY THE PLOTTING
ROUTINES OR COULD BE USED BY THE USER FOR ADDITIONAL
ANALYSIS OF THE CONCENTRATION ESTIMATES.
FLOW DIAGRAM-
INPUFF
PUF0010
PUF0020
PUF0030
PUF0040
PUF0050
PUF0060
PUF0070
PUF0080
PUF0090
PUF0100
PUF0110
PUF0120
PUF0130
PUF0140
PUF0150
PUF0160
PUF0170
PUF0180
PUF0190
PUF0200
PUF0210
PUF0220
PUF0230
PUF0240
PUF0250
PUF0260
PUF0270
PUF0280
PUF0290
PUF0300
PUF0310
PUF0320
PUF0330
PUF0340
PUF0350
PUF0360
POT0370
PUF0380
PUF0390
PUF0400
PUF0410
PUF0420
PUF0430
PUF0440
PUF0450
PUF0460
PUF0470
PUF04SO
PUF0490
PUF0500
PUF0510
PUF0520
PUF0530
PUF0540
PUF0550
PUF0560
POT0570
PUF0580
PUFOS90
PUF0600
PUF0610
PUF0620
PUF0630
PUF0640
PUF0650
PUF0660
PUF0670
PUF0680
PUF0690
POT0700
PUF0710
PUF0720
PUF0730
PUF0740
PUF0750
86
6-84
-------
76.
77.
78.
79.
80.
81.
82.
83.
84.
85.
86.
87.
88.
89.
90.
91.
92.
93.
94.
95.
96.
97.
98.
99.
100.
101.
102.
103.
104.
105.
106.
107.
108.
109.
110.
111.
112.
113.
114.
115.
116.
117.
118.
119.
120.
121.
122.
123.
124.
125.
126.
127.
128.
129.
130.
131.
132.
133.
134.
135.
136.
137.
138.
139.
140.
141.
142.
143.
144.
145.
146.
147.
148.
149.
150.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
1
READ OPTIONS AND INPUT DATA
1---
Ljwwr wit MH . r L.I\. i (jus
- PLMRS (COMPUTE PLUME RISE)
- PGSIG (COMPUTE SY, SZ FOR PG STABILITY' CLASSES)
- XVY (COMPUTE VIRTUAL DISTANCE FOR SY)
- XVZ (COMPUTE VIRTUAL DISTANCE FOR SZ)
- VTIME (COMPUTE VIRTUAL TRAVEL TIMES)
1
1- JSIS1G (COMPUTE SY, SZ BASED ON TRAVEL TIME)
LOOP ON EACH ! STEP
- XVY
- XVZ
- VTIME
- ADVECT (IF IADT NE 0 READ MET. FILE AND COMPUTE
GRID COORDINATES OF PUFFS)
- PLMRS
1
*
* *
* •
* PG •
NO-* CURVES*
* «
• *
* *
*
i
YES
1
1- XVY
1 - XVZ
1
BASED ON 6 POSSIBLE VALUES OF KEYP CALL
1
1- PGSIG
1
1 AND/OR
1
1- LTSIG (COMPUTE SY AS SQRT(TIME) )
1 — VT 1 \TP
I ~ VI 1 Ivui
1
1
•BASED ON 6 POSSIBLE VALUES OF KEYP CALL
1
1- JSISIG
1
I AND/OR
1
1- LTSIG
CJJN(_£.N
i
1 1 LUOr ON HhChrTORS
1 1
1 I---I PUFF LOOP
1 1 1
PUF0760
PUF0770
PUF0780
PUF0790
PUF0800
PUF0810
PUF0820
PUF0830
PUF0840
PUF0850
PUF0860
PUF087U
PUF0880
PUF0890
PUF0900
PUF0910
PUF0920
PUF0930
PUF0940
PUF0950
PUF0960
PUF0970
PUF0980
PUF0990
PUF 1000
PUF1010
PUF 1020
PUF1030
PUF1040
PUF1050
PUF 1060
PUF 1070
PUF 1080
PUF 1090
PUF1100
PUF 1110
PUF1 12U
PUF 1130
PUF 1140
PUF1 150
PUF1 160
PUF1 170
PUF1180
PUF1 190
PUF1200
PUF 1210
PUF 12 20
PUF 12 30
PUF 1240
PUF 12 50
PUF1260
PUF1270
PUF 12 80
PUF 1290
PUF 1300
PUF1310
PUF 13 20
PUF1330
PUF1340
PUF 13 50
PUF1360
PUF1370
PUF1380
PUF1390
PUF1400
PUF 1410
PUF 1420
PUF1430
PUF1440
PUF 1450
PUF1460
PUF 14 70
PUF1480
PUF1490
PUF1500
87
6-84
-------
151.
152.
153.
154.
155.
156.
157.
158.
159.
160.
161.
162.
163.
164.
165.
166.
167.
168.
169.
170.
171.
172.
173.
174.
175.
176.
177.
178.
179.
180.
181.
182.
183.
184.
185.
186.
187.
188.
189.
190.
191.
192.
193.
194.
195.
196.
197.
198.
199.
200.
201.
202.
203.
204.
205.
206.
207.
208.
209.
210.
211.
212.
213.
214.
215.
216.
217.
218.
219.
220.
221.
222.
223.
224.
225.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
I
I I COMPUTE CONCENTRATIONS EACH 1STEP AND SUM.
I I
I---END OF PUFF LOOP
I
END OF RECEPTOR LOOP
I COMBINE PUFFS
| END OF [STEP LOOP
-END OF MET. LOOP
STOP
INTEGER DISKEY
DIMENSION CONCT(25),PP(7),ACON(25),CONC1(25),CONCIA(25),
1TEMP(144),CDIS(144),FRMT(18),ALP(18),TACON(25),FRMAT(18)
COMMON /XP/ XPUFF(600) ,YPUFF(600) ,ZPUFF(600), ITPUFF(600),
1XDSP(600),QPUFF(600),SY(600),SZ(600),CONC(25),NREC,KEYP(600)
COMMON AVEA/ WDIR( 144 ) ,WSPD( 144 ) , HL( 144 ) , SGTH( 144 ) , SGPH( 1 44 ) ,
1PREPU44),XYOP(600) ,XZOP(600>,I TIME,KST(144),STANUM,SW( 144),
2SV(144),UPLM(144)
COMMON /STA/ QP.HPP,TSP,VSP,DP,VFP,SYOP.SZOP.DH
COMMON /REC/ XREC(25),YREC(25),ZREC(25).ALPHA,ANHGT
COMMON /ADVCT/ PDIR( 600 ), PSPD( 600 ), XSWC, YSWC, NUMX, NUMY, DX,DY, FRMAT
DATA PP/0.07,0.07,0.10,0.15,0.15,0.35,0.55/.MAXPUF/600/
PI = 3.141593/180.
NPUFF = 1
SDIR = 0.
SSPD = 0.
READ (5,390) ALP
READ (5,») IOPT,IOPS,IADT,IP22,IPCC,IPIC
I OPT
IOPS
IADT
IP22
IPCC
IPIC
XSWC
YSWC
NUMX
NUMT
DX
DY
(DIMENSIONLESS)
(D1MENS1ONLESS)
STACK DOWNWASH OPTION
IF IOPT=0 NO DOWNWASH
UPDATE SOURCE CHARACTERISTICS
EACH MET. PERIOD. IF IOPS=0 NO UPDATE
USER SUPPLIED WIND FIELD. IF IADT=0 (DIMENSIONLESS)
NO WIND FIELD DATA.
UNIT 22 OUTPUT OPTION. IF 1P22=0 (DIMENS IONLESS)
NO OUTPUT TO UNIT 22.
OPTION TO PRINT OUT PUFF INFORMATION;DIMENSIONLESS)
EACH ITIME. IF 1PCC=0 NO PRINT OUT.
OPTION TO PRINT INTERMEDIATE (DIMENSIONLESS)
CONCENTRATIONS. IF IPIC=0 NO PRINT OUT.
IF (IADT.NE.O) READ (5,») XSWC,YSWC,NUMX,NUMY,DX,DY
EAST-WEST COORDINATE OF THE SW ( DIMENS IONLESS )
CORNER OF METEOROLOGICAL REGION
NORTH-SOUTH COORDINATE OF THE SW
CORNER OF METEOROLOGICAL REGION
NUMBER OF GRID SQUARES IN EAST-WEST (DIMENS IONLESS)
DIRECTION
NUMBER OF GRID SQUARES IN NORTH-
SOUTH DIRECTION
EAST-WEST WIDTH OF GRID SQUARE
NORTH-SOUTH WIDTH OF GRID SQUARE
WIND SPEED AND DIRECTION ARE READ IN FOR EACH GRID
SQUARE, A ROW AT A TIME, FROM WEST TO EAST (LEFT TO
RIGHT). ROWS ARE READ FROM SOUTH TO NORTH (BOTTOM TO
TOP) .
IF (IADT.NE.O) READ (5,390) FRMAT
FRMAT IS THE FORMAT OF FILE 21 (MET. DATA).
READ (5,*) NTIME,ITIME,I STEP,INC,ISIM,ISTTIM,IW,DISKEY,NREC,ANHGT
(DIMENSIONLESS)
(DIMENSIONLESS)
(KM)
(KM)
NTIME NUMBER OF PERIODS OF SIMULATION (DIMENS IONLESS )
ITIME SIMULATION TIME (SECONDS)
I STEP TIME BETWEEN PUFF RELEASES (SECONDS)
IF ISTEP IS NEGITIVE A VALUE FOR ISTEP WILL BE
COMPUTED BASED ON THE STABILITY CLASS, WIND SPEED
AND MINIMUM DISTANCE FROM SOURCE TO RECEPTOR(GDIS).
INC SAMPLING TIME FOR CONCENTRATIONS (SECONDS)
ISIM TIME TO START CONCENTRATION (SECONDS)
CALCULATIONS
ISTTIM A MULTIPLE OF ISTEP FOR PUFF (DIMENSIONLESS)
COMBINATIONS.
PUF1510
PUF1520
PUF1530
PUF1540
PUF1550
PUF1560
PUF1570
PUF1580
PUF1590
PUF1600
PUF1610
PUF1620
PUF1630
PUF1640
PUF1650
PUF1660
PUF1670
PUF1680
PUF1690
PUF1700
PUF1710
PUF1720
PUF1730
PUF1740
PUF1750
PUF1760
PUF1770
PUF1780
PUF1790
PUF1800
PUF1810
PUF1820
PUF1830
PUF1840
PUF1850
PUF1860
PUF1870
PUF1880
PUF1890
PUF1900
PUF1910
PUF1920
PUF1930
PUF1940
PUF1950
PUF1960
PUF1970
PUF1980
PUF1990
PUF2000
PUF2010
PUF2020
PUF2030
PUF2040
PUF2050
PUF2060
PUF 2070
PUF2080
PUF2090
PUF2100
PUF 21 It)
PUF2120
PUF2130
PUF2140
PUF2150
PUF2160
PUF2170
PUF2180
PUF2190
PUF2200
PUF2210
PUF2220
PUF2230
PUF2240
PUF2250
88
6-84
-------
226.
227 .
228 .
229.
230.
231 .
232.
233.
234.
235.
236.
237.
238.
239.
240.
241.
242.
243.
244.
245.
246.
247.
248.
249 .
250.
251.
252.
253.
254.
255.
256.
257 .
258.
259 .
260.
261 .
262.
263.
264.
265.
266.
267.
268.
269.
270.
271.
272.
273.
274.
275.
276.
277 .
278.
279.
280.
281 .
282.
283.
284.
285.
286.
287.
288.
289 .
290.
291.
292.
293.
294.
295.
296 .
297.
298.
299 .
300.
£
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
5
c
c
c
c
c
£
c
IW UNIT NUMBER FOR WRITE STATEMENTS (DIMENSIONLESS)
DISKEY DISPERSION OPTION (DIMENSIONLESS)
DISKEY=1 FOR PG CURVES
DISKEY=2 FOR TRAVEL TIME CURVES
NREC NUMBER OF RECEPTORS (DIMENSIONLESS)
ANHGT ANEMOMETER HEIGHT (M)
READ (5,») XGRID, YGRID, XSIZE, YSIZE, STANUM, ALPHA, SYMAX
XGRID EAST-WEST COORDINATE OF (KM)
S.W. CORNER OF MODEL REGION
YGRID NORTH-SOUTH COORDINATE OF S.W. (KM)
CORNER OF MODEL REGION
XSIZE EAST -WEST SIZE OF MODEL REGION (KM)
YSIZE NORTH -SOUTH SIZE OF MODEL REGION (KM)
STANUM THE NUMBER OF SIGMA Y'S THE PUFF IS (DIMENSIONLESS)
FROM THE RECEPTOR SUCH THAT NO
CONCENTRATION ESTIMATE IS MADE
ALPHA FRACTION OF CROSSWIND DISPERSION (DIMENSIONLESS)
FOR PUFF COMBINATION
SYMAX MAXIMUM SIZE OF SIGMA Y BEFOR GOING (M)
TO LTSIG ROUTINE.
READ (5,«) XSORC, YSORC, QP, HPP, TSP, DP, VSP, VFP, SYOP, SZOP, SDIR.SSPD
XSORC X COORDINATE OF SOURCE (KM)
YSORC Y COORDINATE OF SOURCE (KM)
QP EMISSION RATE (G/SEC)
HPP HEIGHT OF RELEASE (M)
TSP STACK GAS TEMPERATURE (KELVIN)
DP STACK DIAMETER (M)
VSP STACK GAS VELOCITY (M/SEC)
VFP STACK GAS VOLUME FLOW (\I**3/SEC)
SYOP INITIAL SIGMA Y (M)
SZOP INITIAL SIGMA Z (M)
SDIR SOURCE DIRECTION (DEGREES)
SSPD SOURCE SPEED (M/SEC)
READ (5,*) (XREC( I ) ,YREC( I ) ,ZREC( I ) , 1=1, NREC)
XREC X COORDINATE OF RECCEPTOR (KM)
YREC Y COORDINATE OF RECEPTOR (KM)
ZREC Z COORDINATE OF RECEPTOR (M)
READ (5,390) FRMT
FORMAT FOR METEOROLOGICAL DATA
READ (S.FRMT) (WDIRt I ) , WSPD( I ) , HL( I ) , KST( I ) , SGPIK I ) , SGTH( I ) ,
1TEMP( I ) ,CDIS( I ) , 1=1, NTIME)
WDIR WIND DIRECTION (DEGREES)
WSPD WIND SPEED (M/SEC)
HL MIXING HEIGHT (M)
KST STABILITY CLASS (DIMENSIONLESS)
SGPH SIGMA PHI (RADIANS)
SGTH SIGMA THETA (RADIANS)
TEMP AIR TEMPERTURE (K)
GDIS MINIMUM DISTANCE SOURCE TO RECEPTOR (KM)
1F( ALPHA. LT. 0. ) ALPHA=1.0
IF(STANUM.LT.O. ) STANUM=3.0
IF( ISTTIM.LT.O) ISTTIM=3
ISTSA = I STEP
SYTT = SYOP
SZTT = SZOP
DO 5 J=l,25
TACON ( J ) = 0 .
CONT I NUE
WRITE ( IW.370) ALP
WRITE ( IW.380)
WRITE (IW.400) IOPT, IOPS, IADT, IP22, IPCC, IPIC
WRITE (IW.410) ISIM, ISTTIM, STANUM, ALPHA, SYMAX, ANHGT
XGRID2 = XGRID + XSIZE
YGRID2 = YGRID + YSIZE
TD = 0.
IF (IP22.NE.O) WRITE (22) NTIME , ITIME, NREC
THE FOLLOWING INFORMATION IS WRITTEN TO UNIT 22 IF
IP22 .NE. 0.
WRITE NTIME, ITIME, NREC
WRITE XREC, YREC, ZREC
( 1 WRITE FOR EACH OF NREC NUMBER OF RECEPTORS)
WR I TF TT I STPP 1
nrvi IE, iifioidr ——— — — |
WRITE ITM,(CONC( I ) , 1 = 1, NREC) NTIME NUMBER
PUF2260
PUF2270
PUF2280
PUF2290
PUF2300
PUF2310
PUF2320
PUF2330
PUF2340
PUF2350
PUF2360
PUF2370
PUF2380
PUF2390
PUF2400
PUF2410
PUF2420
PUF2430
PUF2440
PUF2450
PUF2460
PUF2470
PUF2480
PUF2490
PUF2500
PUF2510
PUF2520
PUF2530
PUF 2 540
PUF2550
PUF2560
PUF2570
PUF2580
PUF2590
PUF2600
PUF2610
PUF2620
PUF2630
PUF2640
PUF2650
PUF2660
PUF 2670
PUF2680
PUF2690
PUF2700
PUF2710
PUF2720
PUF2730
PUF2740
PUF2750
PUF2760
PUF2770
PUF2780
PUF279U
PUF2800
PUF2810
PUF2820
PUF 2830
PUF2840
PUF2850
PUF2860
PUF2870
PUF2880
PUF2890
PUF2900
PUF2910
PUF2920
PUF2930
PUF2940
PUF2950
PUF2960
PUF2970
PUF2980
Of TC 9 Q Q fl
rUr L y y u
PUF3000
89
6-84
-------
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
301.
302.
303.
304.
305.
306.
307.
308.
309 .
310.
311 .
312.
313.
314.
315.
316.
317.
318.
319.
320.
321.
322.
323.
324.
325.
326.
327 .
328.
329.
330.
331.
332.
333.
334.
335.
336.
337.
338.
339.
340.
341.
342.
343.
344.
345.
346.
347.
348.
349.
350.
351.
352.
353.
354.
355.
356.
357.
358.
359.
360.
361.
362.
363.
364.
365.
366.
367.
368.
369.
370.
371.
372.
373.
374.
375.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
10
20
C
C
30
C
C
40
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
(1 WRITE EACH ISTEP) OF WRITES PUF3010
WRITE NPUFF (it OF PUFFS) I PUF3020
WRITE XPUFF,YPUFF,ZPUFF,SY.SZ I PUF3030
(NPUFF NUMBER OF WRITES) I PUF3040
PUF3050
THE FOLLOWING CODE CAN BE USED TO READ UNIT 22. PUF3060
KEAD(22) NTIME,ITIME.NREC PUF3070
DO 10 I=1,NREC PUF3080
10 READ(22) XREC(I),YREC(I),ZREC(I) PUF3090
20 READ(22) IT,ISTEP PUF3100
30 READ(22) ITM,(CONCtK),K=1,NREC) PUF3110
IF (ITM.LT.IT'ITIME) GO TO 30 PUF3120
READ(22) NPUFF PUF3130
DO 40 NP=1,NPUFF PUF3140
40 READ(22) XPUFF(NP),YPUFF(NP),ZPUFF(NP),SY(NP),SZ(NP) PUF3150
IF (IT.LT.NTIME) GO TO 20 PUF3160
IF (IP22.EQ.O) GO TO 20 PUF317U
DO 10 J=1,NREC PUF318U
WRITE (22) XREC(J),YREC(J),ZREC(J) PUF3190
CONTINUE PUF3200
CONTINUE PUF3210
BEGIN DO LOOP FOR MET. PERIODS PUF3220
DO 360 IT=1,NTIME PUF3230
IF (IT.EQ.l) GO TO 30 PUF3240
IF (1OPS.NE.O) READ (5,*) QP.HPP,TSP,DP,VSP,VFP,SYOP, PUF3250
1 SZOP.SDIR.SSPD PUF3260
IF IOPS .NE. 0 READ NEW SOURCE DATA. PUF3270
CONTINUE PUF3280
IF (IOPS.EQ.O) SYOP = SYTT PUF3290
IF (IOPS.EQ.O) SZOP = SZTT PUF3300
ISTEP = ISTSA PUF3310
IF (IT.EQ.l) GO TO 40 PUF3320
IF (HL(IT).LT.HLtIT-1).AND.KSTtIT).LE.4) HL(IT) = HL(IT-l) PUF3330
UNDER NEUTRAL OR UNSTABLE CONDITIONS THE MIXED LAYER PUF3340
CANNOT DECREASE WITH TIME. PUF3350
CONTINUE PUF3360
KS = KST(IT) PUF3370
PREP(IT) = PP(KS) PUF3380
TD = ITIME PUF3390
CALL PLMRS (KS,WSPD(IT),TEMP(IT),PREP(IT),IOPT,ANHGT,DH,HE) PUF3400
DH=PLUME RISE PUF3410
HE=EFFECTIVE STACK HEIGHT PUF34'0
PHGT = HE PUF3430
EXTRAPOLATE WIND SPEED UP TO PLUME HEIGHT OR 200 \I. PUF3440
WHICH EVER IS LOWER PUF3450
IF (HE.GT.200.) PHGT = 200. PUF3460
USCAL = (PHGT/ANHGTC'PREPt IT) PUF3-170
IF (ANHGT.GT.PHGT) USCAL = 1. PUF3480
ANEMOMETER HEIGHT SHOULD BE NEAR THE SURFACE. IF PUF349U
THE WIND SPEEDS GIVEN ARE APPROPRIATE FOR HEIGHTS PUF3500
ABOVE THE PLUME HEIGHT NO EXTRAPOLATIONS ARE PUF3510
PERFORMED. PUF 3 5 2 0
UPLM(IT) = WSPD(1T)'USCAL PUF3530
SV(IT) = SGTH(IT)*UPLM(IT) PUF3540
SW(IT) = SGPH(IT)*UPLM(IT) PUF3550
USPD = -UPLM(IT)*SIN(WDIR(IT)'PI) PUF356U
VSPD = -UPLM(IT)»COS(WDIR(IT)«PI) PUF3570
XSPD = SSPD«SIN(SDIR'PI) PUF3580
YSPD = SSPD*COS(SDIR*PI) PUF3590
USPD U COMPONENT OF WIND. PUF3600
VSPD V COMPONENT OF WIND. PUF3B10
XSPD X COMPONENT OF SOURCE SPEED. PUF3620
YSPD Y COMPONENT OF SOURCE SPEED. PUF3630
CALL PGSIG (GDIS,GDIS,KS.SYT.SZT) PUF3640
DPLM = SQRT((USPD-XSPD)«'2 + (VSPD-YSPD)••2) PUF3650
IF (DPLM.LT.UPLMtIT)) DPLM = UPLMtIT) PUF3660
SRT = (2.*SYT)/DPLM PUF3670
ISRT = SRT PUF3680
IF (ISTEP.GT.O) GO TO 100 PUF3690
IF ISTEP IS NEGITIVE A VALUE FOR ISTEP WILL BE PUF3700
COMPUTED BASED ON THE STABILITY CLASS, WIND SPEED PUF3710
AND MINIMUM DISTANCE FROM SOURCE TO RECEPTOR (CD I S ) . PUF3720
IF (ISRT.GE.INC) GO TO 90 PUF3730
IF (ISRT.LE.l) GO TO 80 PUF3740
DO 60 1=1,ISRT PUF3750
90
6-84
-------
2
2
2
2
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
1
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
376.
377.
378.
379.
380.
381.
382.
383.
384.
385.
386.
387.
388.
389.
390.
391.
392.
393.
394.
395.
396.
397.
398 .
399.
400.
401.
402.
403.
404.
405.
406.
407 .
408.
409 .
410.
411 .
412.
413.
414.
415.
416.
417.
418.
419.
420.
421.
422.
423.
424.
425.
426.
427.
428.
429.
430.
431.
432.
433.
434.
435.
436.
437.
438.
439.
440.
441 .
442.
443.
444.
445.
446.
447 .
448.
449.
450.
60
70
80
90
100
C
C
C
C
C
C
C
C
C
C
110
C
120
C
C
C
C
IRR = ISRT - I + 1
IR = MODI INC, IRR)
IF ( IR.EQ.O) GO TO 70
CONT I NUE
[STEP = IRR
GO TO 100
I STEP = 1
GO TO 100
1STEP = INC/ 2
CONTINUE
IF ( IT.NE. 1 ) GO TO 110
ITPUFF(l) = 0
XPUFF( 1 ) = XSORC'1000.
YPUFF(l) = YSORC* 1000.
ZPUFFd ) = HE
KEYP(l) = 1
KEYP( NPUFF )=1 PUFF BELOW MIXED LAYER AND SHORT RANGE
KEYP(NPUFF)=2 USE SIGMAS APPROPRIATE FOR LONG TRAVEL TIMES.
KEYP(NPUFF)=3 PUFF ABOVE MIXED LAYER
KEYP(NPUFF)=4 PUFF WELL MIXED AND SHORT RANGE
KEYP(NPUFF)=5 PUFF ABOVE MIXED LAYER AND LONG RANGE
KEYP(NPUFF)=6 LONG TRAVEL TIME AND WELL MIXED.
SY(1) = SYOP
SZ(l) = SZOP
IF (ZPUFF(l) .GT.HU IT)) KEYP(l) = 3
QPUFF EMISSIONS DURING PERIOD I STEP.
XPUFF.YPUFF.ZPUFF PUFF COOR.
QPUFF ( 1 ) = QP'ISTEP
KS = KST( IT)
IF (KEYP(1).EQ.3) KS = 7
COMPUTE VIRTUAL DISTANCE
XYOP(l) = XVY(SYOP.KS)
XZOP(l) = XVZ(SZOP.KS)
XDSP(l) = 0.
IF (DISKEY.EQ. 1) GO TO 1 1 0
IF (KEYP( 1) .EQ.3.0R.KEYP( D.EQ.5) KS = 7
COMPUTE VIRTUAL TIME
CALL VTIME ( KEYP( 1 ) , IT , KS , SY( 1 ) , SZ( 1 ) ,XYOP( 1 ) , XZOP( 1 ) )
CONTINUE
ITMS = 0
WRITE (IW.420) I TPUFF ( NPUFF ) ,QP,HPP,TSP, VSP, DP, VFP, XSORC, YSORC,
ISSPD.SDIR.HE
WRITE ( IW.430) WD1R( IT) ,WSPD( IT) ,HL( IT) ,PREP( IT) ,KST( IT) ,UPLM( IT) ,
1TEMPI IT) ,SYOP,SZOP,SGTH( IT) ,SGPH( IT)
ISTRP = IT'ITIME - IT1ME
ISTOP = IT'ITIME
WRITE (IW.440) ISTRP, ISTOP, ITIME, ISTEP,CDIS( IT) ,DISKEY
DO 120 1=1, NREC
ZERO CONCENTRATION ARRAYS.
CONCTl I ) - 0.
CONC( I) =0.
CONCI(I) = 0.
CONTINUE
IST10 = ISTTIM'ISTEP
IF (IP22.NE.O) WRITE (22) IT.ISTEP
BEGIN DO LOOP FOR I STEP WITHIN LOOP FOR MET. PERIOD
DO 350 I = ISTEP, ITIME, I STEP
ITRAV = I
NPUFF = NPUFF + 1
IF ( NPUFF. GT.MAXPUF) STOP 100
XSORC = XSORC + (XSPD'ISTEP) / 1000.
YSORC = YSORC + ( YSPD* I STEP ) / 1 000 .
ADVECT SOURCE COOR.
XPUFFl NPUFF) = XSORC'1000.
YPUFF( NPUFF) = YSORC* 1000.
XPUFF.YPUFF NEW PUFF COOR.
ZPUFFt NPUFF) = HE
QPUFF (NPUFF) = QP'ISTEP
KEYP( NPUFF) = 1
SY( NPUFF) = SYOP
SZ(NPUFF) = SZOP
IF (ZPUFF(NPUFF) .GT.HL( IT) ) KEYP(NPUFF) = 3
QPUFF IS EMISSION RATE OF CURRENT PUFF.
KS = KST( IT)
IF (KEYP(NPUFF) .EQ.3) KS = 7
PUF3760
PUF3770
PUF3780
PUF3790
PUF3800
PUF3810
PUF3820
PUF3830
PUF3840
PUF3850
PUF3860
PUF3870
PUF3880
PUF3890
PUF3900
PUF3910
PUF3920
PUF3930
PUF3940
PUF3950
PUF3960
PUF3970
PUF3980
PUF3990
PUF4000
PUF4010
PUF4020
PUF4030
PUF4040
PUF4050
PUF4060
PUF4070
PUF4080
PUF4090
PUF4100
PUF41 10
PUF4120
PUF 4 130
PUF4140
PUF4150
PUF4160
PUF 4170
PUF4180
PUF4190
PUF4200
PUF4210
PUF4220
PUF 4230
PUF4240
PUF4250
PUF4260
PUF427U
PUF4280
PUF4290
PUF4300
PUF4310
PUF4320
PUF4330
PUF4340
PUF4350
PUF4360
PUF4370
PUF4380
PUF4390
PUF4400
PUF4410
PUF4420
PUF 4430
PUF4440
PUF4450
PUF4460
PUF4470
PUF4480
PUF4490
PUF4500
91
6-84
-------
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
2
2
2
2
2
2
2
2
2
2
3
3
3
3
3
3
2
2
2
3
3
3
3
3
3
3
3
3
3
3
3
3
3
2
2
2
2
3
3
3
2
451.
452.
453.
454.
455.
456.
457.
458.
459.
460.
461.
462.
463.
464.
465.
466.
467.
468.
469.
470.
471.
472.
473.
474.
475.
476.
477 .
478.
479.
480.
481.
482.
483.
484.
485.
486.
487.
488.
489.
490.
491.
492.
493.
494.
495.
495.5
496.
497.
498.
499.
500.
501.
502.
503.
504.
505.
506.
507 .
508.
509.
510.
511.
512.
513.
514.
515.
516.
517.
518.
519.
520.
521.
522.
523.
524.
525.
130
C
C
140
150
160
C
170
180
C
190
210
C
220
C
XYOP(NPUFF) = XVY(SYOP.KS)
XZOP(NPUFF) = XVZ(SZOP.KS)
IF (DISKEY.EQ.1) GO TO 130
IF (KEYP(NPUFF).EQ.S.OR.KEYP(NPUFF).EQ.5) KS = 7
CALL VTIME (KEYP(NPUFF),IT.KS,SY(NPUFF).SZfNPUFF),
1 XYOP(NPUFF),XZOP(NPUFF))
XDSP(NPUFF) = 0.
ITM = ITM + I STEP
ITPUFF(NPUFF) = ITM
ITPUFF TIME OF RELEASE OF CURRENT PUFF.
NN = NPUFF - 1
ADVECT ALL PREVIOUS PUFFS.
IF (IADT.EQ.O) GO TO 140
CALL ADVECT (NPUFF,ITRAV,I STEP)
IF (XSPD.LT.0.01) GO TO 140
UZP = PSPD(NPUFF)
IF (ABS(WSPD(IT)-UZP).LE.0.1) GOTO 140
CALL PLMRS(KS,UZP,TEMPI IT),PREP(IT), IOPT,ANHGT,DH,HE)
ZPUFF(NPUFF) = HE
WSPD(IT) = UZP
PHGT = HE
IF (HE.GT.200.) PHGT = 200.
USCAL = (PHGT/ANHGT)»«PREP(IT)
IF (ANHGT.GT.PHGT) USCAL = 1.
UPLM(IT) = WSPD(IT)'USCAL
CONTINUE
DO 160 J=1,NN
L = NN * 1 - J
IF (IADT.EQ.O) GO TO 150
USPD = -SIN(PDIR(L)»PI)«PSPD(L)»USCAL
VSPD = -COS(PDIR(L)«PI)«PSPD(L)»USCAL
CONTINUE
XPUFF(L) = XPUFF(L) + USPD'ISTEP
YPUFF(L) = YPUFF(L) + VSPD*I STEP
CONTINUE
CALL PROCES (NPUFF,ISTEP,DISKEY,IT,ITRAV.SYMAX,IAVT)
IF (ITM.LT.ISIM) GO TO 170
CONC. FROM THE PUFFS AT TIME=ITPUFF.
CALL CONCEN (NPUFF,IST10,I)
IF (IP22.NE.O) WRITE (22) ITM,(CONCtUP),IJP=1,NREC)
CONTINUE
IRM = MOD(I,IST10)
IF (IRM.NE.O) GO TO 210
NPU = NPUFF
DO 180 KPF=1,NPU
IF (QPUFF(KPF).LE.O.) KEYP(KPF) = -1
XTG = XPUFF(KPF)/1000.
YTG = YPUFF(KPF)/1000.
IF (XTG.LT.XGRID.OR.XTG.GT.XGRID2.OR.
1 YTG.LT.YGRID.OR.YTG.GT.YGRID2) KEYP(KPF) = -1
CONTINUE
L = 0
DO 190 KPF = l.NPU
COMPUTE NUMBER OF PUFFS AFTER PUFF ELIMINATION.
IF (KEYP(KPF).LT.O) GO TO 190
L + 1
XPUFF(KPF)
YPUFF(KPF)
ZPUFF(KPF)
QPUFF(KPF)
ITPUFF(KPF)
SY(KPF)
SZ(KPF)
XYOP(KPF)
XZOP(KPF)
KEYP(KPF)
XDSP(KPF)
XPUFF(L)
YPUFF(L)
ZPUFF(L)
QPUFF(L)
ITPUFF(L)
SY(L)
SZ(L)
XYOP(L)
XZOP(L)
KEYP(L)
XDSP(L)
CONTINUE
NPUFF = L
CONTINUE
SUMMING CONC. FOR SAMPLING PERIOD, I.E. 1 HR.
DO 220 L=1,NREC
CONCT(L) = CONCT(L) + CONC(L)
CONTINUE
SUMMING CONC. FOR INTER. SAMPLING PERIODS,I.E. 5 MIN.
DO 230 L=1,NREC
PUF4510
PUF4520
PUF4530
PUF4540
PUF4550
PUF4560
PUF4570
PUF4580
PUF4590
PUF4600
PUF4610
PUF4620
PUF4630
PUF4640
PUF4650
PUF4660
PUF4670
PUF4680
PUF4690
PUF 4700
PUF4710
PUF4720
PUF4730
PUF4740
PUF4750
PUF4760
PUF4770
PUF4780
PUF4790
PUF4800
PUF 4810
PUF4820
PUF4830
PUF4840
PUF4850
PUF4860
PUF4870
PUF4880
PUF4890
PUF4900
PUF4910
PUF4920
PUF4930
PUF4940
PUF4950
PUF4955
PUF4960
PUF4970
PUF4980
PUF4990
PUF5000
PUF5010
PUF5020
PUF5030
PUF5040
PUF5050
PUF5060
PUF5070
PUF5080
PUF5090
PUF5100
PUF5110
PUF5120
PUF5130
PUF5140
PUF5150
PUF5160
PUF5170
PUF5180
PUF5190
PUF 5200
PUF5210
PUF522D
PUF5230
PUF5240
PUF5250
92
6-84
-------
3 526. CONCI(L) = CONCI(L) + CONC(L) PUF5260
3 .527. 230 CONTINUE POT5270
2 528. IF (I.NE.ITIME) GO TO 290 PUF5280
2 529. IF (IP22.NE.O) WRITE (22) NPUFF PUF5290
2 530. IF (IP22.EQ.O) GO TO 250 PUF5300
2 531. DO 240 INPF=1,NPUFF PUF5310
3 532. WRITE (22) XPUFF(INPF),YPUFF(INPF).ZPUFF(INPF), PUF5320
3 533. 1 SY(INPF),SZ(INPF) PUF5330
3 534. 240 CONTINUE PUF5340
2 535. 250 CONTINUE PUF5350
2 536. IF (IPCC.EQ.O) GO TO 270 POT5360
2 537. WRITE (IW.450) PUF5370
2 538. DO 260 IPR=1,NPUFF PUF5380
3 539. WRITE (IW.460) IPR,XPUFF(I PR),YPUFFlI PR), PUF5390
3 540. 1 ZPUFF(IPR),ITPUFF(IPR),QPUFF(IPR),SY(IPR), PUF5400
3 541. 2 SZ(IPR),KEYP(IPR) PUF5410
3 542. 260 CONTINUE PUF5420
2 543. 270 CONTINUE PUF5430
2 544. C IF I=ITIME END OF SIMULATION PERIOD. PUF5440
2 545. DIV = FLOAT( ITIME)/FLOAT! ISTEP) PUF5450
2 546. C AVG. CONC. FOR SAMPLING TIME. PUF5460
2 547. DO 280 L=1,NREC PUF5470
3 548. ACON(L) = CONCT(L)/DIV PUF5480
3 549. TACON(L) = TACON(L) + ACON(L) PUF5490
3 550. 280 CONTINUE PUF5500
2 551. 290 CONTINUE PUF5510
2 552. IRM = MOD(ITM,INC) PUF5520
2 553. IF (IRM.NE.O) GO TO 350 PUF5530
2 554. C IF IRM=0 END OF INTER. SAMPLING PERIOD, I.E. 5 MIN. PUF5540
2 555. DIV = FLOAT (I NO/FLOAT (I STEP) PUF5550
2 556. e AVG. CONC. AT END OF INTER. SAMPLING PERIOD. PUF5560
2 557. DO 300 J=1,NREC PUF5570
3 558. CONCIA(J) = CONCI(J)/DIV PUF5580
3 559. 300 CONTINUE PUF5590
3 560. C ZERO CONC. ARRAYS FOR INTER. SAMPLING TIMES. PUF5600
2 561. DO 310 L=1,NREC POT5610
3 562. CONCI(L) = 0. PUF5620
3 563. 310 CONTINUE PUF5630
2 564. IF (IPIC.EQ.O) GO TO 330 PUF5640
2 565. ITMV1 = ITM - INC PUF5650
2 566. WRITE (IW.470) INC, ITMV1, ITM PUF5660
2 567. DO 320 J=1,NREC PUF5670
3 568. WRITE (IW.480) XREC(J),YRECtJ),ZREC(J),CONCIAtJ) PUF5680
3 569. 320 CONTINUE PUF5690
2 570. 330 CONTINUE PUF5700
2 571. IF (I.NE.ITIME) GOTO 350 PUF5710
2 572. WRITE (IW.490) ITIME,ISTRP,I STOP PUF5720
2 573. DO 340 J=1,NREC PUF5730
3 574. WRITE (IW.480) XREC(J),YREC(J),ZREC(J),ACON(J) PUF5740
3 575. 340 CONTINUE PUF5750
2 576. 350 CONTINUE POT5760
2 577. C END LOOP OVER ISTEP PUF5770
1 578. 360 CONTINUE PUF5780
1 579. C END LOOP OVER MET. PERIOD POT5790
580. DIV = FLOAT(NTIME) PUF5800
581. DO 365 J=1,NREC PUF5810
1 582. TACON(J) = TACON(J)/DIV PUF5820
1 583. 365 CONTINUE POT5830
584. TOTT = (NTIME«ITIME)/3600. PUF5840
585. WRITEt IW.500) TOTT PUF5850
586. DO 367 J=1,NREC PUF5860
1 587. WRITEtIW,480) XREC(J),YREC(J),ZREC(J),TACON(J) PUF5870
1 588. 367 CONTINUE PUF5880
589. STOP PUF5890
590. 370 FORMAT ('1',3X,18A4) PUF5900
591. 380 FORMAT (60X, ' VERS ION 84107') PUF5910
592. 390 FORMAT ( 18A4) PUF5920
593. 400 FORMAT (3X,'O P T I O N S ',10X,'A "1" INDICATES THAT THE OPTION'/ PUF5930
594. 1T34,'HAS BEEN EXERCISED'//3X,'STACK DOWNWASH1,T33,I3/3X, PUF5940
595. 2'UPDATE SOURCE DATA',T33,I3/3X, PUF5950
596. 3'USER SUPPLIED WIND FIELD',T33,I3/3X,'UNIT 22 OUTPUT OPTION1,T33, PUF5960
597. 4I3/3X,'PRINT PUFF INFORMATION',T33,I3/3X,'INTERMEDIATE CONCENTRA1, PUF5970
598.- 5'TIONS' ,T33, I3//) PUF5980
599. 410 FORMAT (3X,'I NPUT PARAMETERS '//3X,1ISIM=',I6/3X, PUF5990
600. 1'ISTTIM=',I6/3X,'STANUM='.F5.2/3X,'ALPHA='.F5.2/3X,'SYMAX=',F9.1/ PUF6000
93
6-84
-------
601. 23X,'ANHGT=',F5.1//) PUF6010
602. 420 FORMAT ('1 * • * SOURCE INFORMATION ***'// PUF6020
603. 1T2,'SOURCE',T11,'STACK',T2i,'STACK',T31,'STACK GAS',T42,'STACK', PUF6030
604. 2T51,'VOLUME',T62,'COORDINATES AT TIME ',I 5,T90, 'SOURCE1 ,T100,'SO', PUF6040
605. 3'URGE' ,T115, 'PLUME' /T2, 'STRENGTH' ,T11, 'HEIGHT' ,T21, 'TEMP. ' ,T31, PUF6050
606. 4'VELOCITY',T41,'DIAMETER1,T52,'FLOW1,T64,'EAST',T77,'NORTH',T90, PUF6060
607. 5'SPEED',T100,'DIRECTION',T115,'HEIGHT'/T3,'(G/SEC)',T12,'(M)',T20, PUF6070
608. 6'(DEG-K)',T32,'(M/SEC)',T43,'(M)',T50,'(M»*3/SEC)',T64,'(KM)',T78, PUF6080
609. 7'(KM) ' ,T90,'(M/SEC)',T100,'(DEC)',T116,'(M)'//T2,E8.3,T11,F6.2, PUF6090
610. 8T21,F7.3,T32,F7.3,T41,F6.3,T51,F8.3,T61,F9.3,T75,F9.3,T91,F7.3, PUF6100
611. 9T101,F5.1,T115,F7.2/) PUF6110
612. 430 FORMAT (///'••• METEOROLOGY •••'//T76,'INITIA1, PUF6120
613. 1'L SIGMAS'/T2,'WIND DIR.',T14,'WIND SPD.',T25,•MIXING HGT.',T39, PUF6130
614. 2'PROF.EP',T50,'STABILITY',T60,' U PLUME',T70,'TEMP1,T79,'(R) (', PUF6140
615. 3'Z)',T91,'SIGMA TH.',T102,'SIGMA PH.'/T4,'(DEC)',T15,'(M/SEC)', PUF6150
616. 4T28,'(M) ',T39,'(DIMEN)' ,T50,'(CLASS)',T61,'(M/SEC)',T71,'(K)',T82, PUF6160
617. 5'(M)',T92,'(RAD.)',T104,'(RAD.)'//T3.F5.1,T16,F6.3,T27,F5.0,T39, PUF6170
618. 6F7.3.T50,I4.T60,F7.3,T70,F5.1,T77,F5.1,T83,F5.1,T91,F6.4,T102, PUF6180
619. 7F6.4/) PUF6190
620. 440 FORMAT (T6,'SIMULATION PERIOD',T26,'SIMULATION TINE',T45,'PUFF R', PUF6200
621. 1'ELEASE RATE1,T65,'SOURCE RECEPTOR DISTANCE',T95,'DISPERSION1/T3, PUF6210
622. 2'START (SEC) STOP (SEC)',T31,'(SEC)',T51, '(SEC)',T75,'(KM) ' ,T98, PUF6220
623. 3'TYPE'/T3,16,T16,I 6,T32,I 4,T5I,I 4,T72,F8.2.T100,I III) PUF6230
624. 450 FORMAT (///' PUFF* ' ,T15 , 'X' ,T26 ,'Y.',T37 ,'Z1 ,T44 ,'TIME' ,T54 , 'TOTA' , PUF6240
625. 1'L Q',T67,'SY',T77,'SZ',T85,'KEYP'/) PUF6250
626. 460 FORMAT (2X,I 4,3X,31F10.3,IX),I 6,3X.F10.2,2F10.3,T85,11) PUF6260
627. 470 FORMAT (///2X.I6,1 SEC',' AVG. CONCENTRATION AT RECEPTORS DURING', PUF6270
628. I1 INTERMEDIATE PERIOD ',16,' TO ',16,' SECONDS'//IOX,'RECEPTORS'/ PUF6280
629. 2T5,'X',T15,'Y1,T25,'Z',T35,'CONCENTRATION
-------
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
666.
667 .
668 .
669.
670 .
671 .
672.
673.
674.
675.
676.
677 .
678.
679.
680.
681.
682.
683 .
684.
685 .
686 .
687.
688 .
689 .
690.
691.
692.
693.
694.
695.
696.
697 .
698 .
699.
700.
701.
702.
703.
704.
705.
706.
707 .
708.
709 .
710.
711.
712.
713.
714.
715.
716.
717.
718.
719.
720.
721.
722.
723.
724.
725.
726.
727.
728.
729 .
730.
731.
732.
733.
734.
735.
736.
737 .
738.
739.
740.
C
C
C
C
C
C
C
10
C
20
30
40
50
C
C
C
C
C
C
60
70
C
SUBROUTINE CONCEN (NPUFF,IST10,ITP)
SUBROUTINE CONCEN COMPUTES THE CONCENTRATION AT
EACH RECEPTOR.
COMMON /WEA/ WDIRI144),WSPD(144),HL(144),SGTH(144),SGPH(144),
1PREP(144),XYOP(600),XZOP(600).ITIME,KST(144),STANUM,SW(144),
2SV(144),UPLM(144)
COMMON fX.fl XPUFFUOO) ,YPUFF(600) ,ZPUFF(600) , ITPUFFUOO ) ,
IXDSPUOO),QPUFF(600),SY(600),SZ(600),CONC(25),NREC,KEYP(600)
COMMON /STA/ QP,HPP,TSP,VSP,DP,VFP,SYOP,SZOP,DH
COMMON /REC/ XRECl25),YRECl25),ZREC(25),ALPHA,ANUGT
PI = 3.141593
PHI2 = 15.7496
RECEPTOR LOOP.
DO 50 IREC=1-NREC
XR = XREC XPUFF(I))/2.
YPUFF(I -1) = (YPUFF(I - 1) * YPUFF(I))/2.
ZPUFF(I -1) = (ZPUFF(I -1) + ZPUFF(I))/2.
QPUFF(I - 1) = QPUFF(I -1) + QPUFF(I)
ITPUFF(I-l) = (ITPUFF(I-l) * ITPUFF(I))/2.
XDSP(I-l) = (XDSP(I-l) + XDSP(I))/2.
QPUFF(I) = - 1.
CONTINUE
N = 0
DO 90 1=1,NPUFF
REMOVE TAGGED PUFFS AND COMPUTE NUMBER OF PUFFS
CON0010
CON0020
CON0030
CON0040
CON0050
CON0060
CON0070
CON0080
CON0090
CON 0100
CON 0110
CON0120
CON0130
CON0140
CON01SO
CON0160
CON0170
CON0180
CON0190
CON0200
CON0210
CON0220
CON0230
CON0240
CON0250
CON0260
CON0270
CON028U
CON0290
CON0300
CON0310
CON0320
CON0330
CON0340
CON0350
CON0360
CON0370
CON0380
CON0390
CON0400
CON0410
CON0420
CON0430
CON0440
CON0450
CON0460
CON0470
CON0480
CON0490
CON0500
CON0510
CON0520
CON0530
CON0540
CON0550
CON0560
CON0570
CON0580
CON0590
CON0600
CON0610
CON0620
CON0630
CON0640
CON0650
CON0660
CON0670
CON0680
CON0690
CON0700
CON0710
CON0720
CON0730
CON0740
CON0750
95
6-84
-------
741.
742.
743.
744.
745.
746.
747.
748.
749.
750.
751.
752.
753.
754.
755.
756.
757 .
758.
C AFTER PUFF COMBINATION
IF (QPUFFU:
80 N
XPUFF(N) =
YPUFF(N) =
ZPUFF(N) =
QPUFF(N) =
ITPUFF(N) =
SY(N) =
SZ(N)
XYOP ( N )
XZOP(N)
KEYP(N)
XDSP(N)
90 CONTINUE
NPUFF = N
100 RETURN
END
)) 90,80
N + 1
XPUFF( 1
YPUFFl 1
ZPUFFi 1
QPUFF ( I
ITPUFF(
SY( I )
SZ( I )
XYOP ( I )
XZOP( I)
KEYP( 1)
XDSP( 1 )
,80
)
)
)
)
1 )
CON0760
OON0770
OON0780
OON0790
CON0800
CON 0810
CON0820
CON0830
CON0840
CON0850
CON0860
CON0870
OON0880
CON0890
CON0900
CON0910
CON0920
CON0930
759.
760.
761.
762.
763.
764.
"765.
766.
767 .
768.
769.
770.
771.
772.
773.
774.
775.
SUBROUTINE JS I S IG(KEY, IT , T , KS , SY, SZ )
C SUBROUTINE JS1SIG COMPUTES SIGMAS BASED ON TRAVEL TIME
C BASED ON IRWIN (198J) WHICH CULMINATES WORK BY DRAXLER
C (1976) AND CRAMER (1976).
COITION /WEA/ WDIR( 144 ) , WSPD( 1 44 ) , HL( 1 44 ) , SGTH( 144 ) , SGPH( 1 44 ) ,
1PREP( 144) ,XYOP(600) ,XZOP(600) , ITIME.KST( 144) ,STANUM,
2SW( 144) ,SV(144) ,UPLM(144)
SIGV = SV( IT)
SIGW = SW( IT)
IF (KEY.EQ.3.OR.KEY.EQ.5) SIGW = 0.01
= 1.0 + 0.9*SQRT(T/1000. )
= SIGV*T«(1. /FY)
FY
SY
FZ = 1.
IF (KS.GT.4) FZ = 1.0
SZ = SIGW*T»(1./FZ)
RETURN
END
0 . 9»SQRT(T/ 50 . )
JSI0010
JSI0020
JSI0030
JSI0040
JSI0050
JSI0060
JSI0070
JSI0080
JSI0090
JSI0100
JSI0110
JSI0120
JSI0130
JSI0140
JSI0150
JSI0160
JSI0170
776. SUBROUT1NL LTS 1G( SYY , UT, SY )
777. C SUBROUTINE LTSIG COMPUTES THL SIGMAS FOR LONG RANGH
778. C TRANSPORT, PROPORTIONAL TO Tilt SQRT(TIME).
779. C BASED ON DRAXLER (1976).
780. TSTR = 9166.6667
781. A = 10.4446
782. DELT = ((SYY/A)*"2 - TSTR) + DT
783. SY = A*SQRT(TSTR " DELT)
784. RETURN
785. END
LFSOOIO
LTS0020
LTS0030
LTS0040
LTS0050
LTS0060
LTS0070
LTS0080
LTS0090
LTS 0100
96
6-84
-------
1
1
1
1
1
1
1
1
786 .
787.
788 .
789.
790.
791 .
792.
793.
794.
795.
796.
797.
798.
799.
800.
801.
802.
803.
804.
805.
806.
807 .
808 .
809.
810.
811.
812.
813.
814.
815.
816.
817.
818.
819 .
820.
821 .
822.
823.
824.
825.
826.
827 .
828.
829 .
830.
831.
832.
833.
834.
835.
836.
837.
838.
839.
840.
841.
842.
843.
844.
845.
846.
847 .
848.
849.
850.
851.
852.
853.
854.
855.
856.
857.
858.
859.
860.
C
C
C
C
C
C
C
C
C
C
C
10
20
30
C
40
50
60
C
70
C
75
C
80
90
100
C
110
120
130
SUBROUTINE PGSIG (X,XY,KST,SY,SZ) PGS0010
D. B. TURNER, ENVIRONMENTAL APPLICATIONS BRANCH PGS0020
METEOROLOGY LABORATORY, ENVIRONMENTAL PROTECTION AGENCY PGS0030
RESEARCH TRIANGLE PARK, N C 27711 PGS0040
(919) 549 - 8411, EXTENSION 4565 PCS0050
VERTICAL DISPERSION PARAMETER VALUE, SZ DETERMINED BY PCS0060
S'l = A * X *• B WHERE A AND B ARE FUNCTIONS OF BOTH STABILITY PGS0070
AND RANGE OF X. PGS0080
HORIZONTAL DISPERSION PARAMETER VALUE, SY DETERMINED BY PGS0090
LOGARITHMIC INTERPOLATION OF PLUME HALF-ANGLE ACCORDING TO PGS0100
DISTANCE AND CALCULATION OF 1/2.15 TIMES HALF-ARC LENGTH. PGS0110
DIMENSION XA(7), XB( 2 ) , XD(5), XE(8), XF ( 9 ) , AA(8), BA( 8 ) , ABO), PGS0120
1BB13), AD16), BD(6), AE(9), BE19), AF(10), BF(10) PGS0130
DATA XA /.5,.4, .3, .25,.2, .15,.I/ PGS0140
DATA XB /.4,.27 PGS0150
DATA XD /30. , 10. ,3. , 1 . , .3/ PGS01SO
DATA XE /4 0. ,20. ,10. ,4. , 2 . , 1 . , .3,. 1 / PGS 0170
DATA XF 760.,30.,15.,7.,3.,2.,1.,.7,.2/ PGS0180
DATA AA 7453.85,346.75,258.89,217.41,179.52,170.22, 158.08,122.8, PGS0190
DATA BA II.1166,1.7283,1.4094,1.2644,1.1262,1.0932,1.0542, .9447 / PGS0200
DATA AB 7109.30,98.483,90.6737 PGS0210
DATA BB 71.0971,0.98332,0.931987 PGS0220
DATA AD 744.053,36.650,33.504,32.093,32.093,34.4597 PGS0230
DATA BD 70.51179,0.56589,0.60486,0.64403,0.81066,0.86974.' PGSOJ40
DATA AE 747.618,35.420,26.970,24.703,22.534,21.628,21.628,23.331, PGS0250
124.267 PGS0260
DATA BE 70.29592,0.37615,0.46713,0.50527,0.57154,0.63077,0.75660, PGS0270
10.81956,0.83667 PGS 0280
DATA AF 734.219,27.074,22.651,17.836,16.187,14.823,13.953,13.953, PGS0290
114.457,15.2097 PGS030U
DATA BF 70.21716,0.27436,0.32681,0.41507,0.46490,0.54503,0.63227, PGS031U
10.68465,0.78407,0.81558/ PGS0320
GOTO (10,40,70,75,80,110,140), KST PGS0330
STABILITY A (10) PGS0340
TH = (24.167 - 2.5334*ALOG(XY))757.2958 PGS035U
IF (X.GT.3.11) GO TO 170 PGS 0 3 6 0
DO 20 ID=1,7 PGS037U
IF (X.GE.XAtID)) GO TO 30 PGS038U
CONTINUE PGS0390
ID = 8 PGS0400
SZ = AA(ID)*X«*BA(ID) PGS0410
GO TO 190 PGS0420
STABILITY B (40) PGS04JO
TH = (18.333 - 1.8096*ALOG(XY))/57.2958 PGS0440
IF (X.GT.35.) GO TO 170 PGS0450
DO 50 ID=1,2 PGS0460
IF (X.GE.XBlID)) GO TO 60 PGS0470
CONTINUE PGS0480
ID = 3 PGS0490
SZ = AB(ID)»X»0BB(ID) PGS0500
GO TO 180 PGS0510
STABILITY C (70) PGS0520
TH = (12.5 - 1.0857»ALOG(XY))/57.2958 PGS0530
SZ = 61.141»X*»0.91465 PGS0540
GO TO 180 PGS0550
D DAY TIME PGS0560
TH = (8.3333 - 0.72382»ALOG(XY))757.2958 PGS0570
SZ = 30.9057*X»*0.8273 PGS0580
GO TO 180 PGS0590
STABILITY D (80) PGS0600
TH = (8.3333 - 0.72382»ALOG(XY))757.2958 PGS0610
DO 90 1D=1,5 PGS0620
IF (X.GE.XDtID)) GO TO 100 PGS0630
CONTINUE PGS0640
ID = 6 PGS0650
SZ = AD(ID)*X«»BD(ID) PGS0660
GO TO 180 PGS0670
STABILITY E (110) PGS0680
TH = (6.25 - 0.54287»ALOG(XY))757.2958 PGS0690
DO 120 ID=1,8 PGS0700
IF (X.GE.XE(ID)) GOTO 130 PGS0710
CONTINUE PGS0720
ID = 9 PGS0730
SZ = AE(ID)«X»*BE(ID) PGS0740
GO TO 180 PGS0750
97
6-84
-------
861. C STABILITY F (140) PGS076G
862. 140 TH = (4.1667 - 0.36191»ALOG(XY))/57.2958 PGS0770
863. DO 150 ID=1,9 PGS0780
1 864. IF (X.GE.XF(ID)) GOTO 160 PGS0790
1 865. 150 CONTINUE PGS0800
866. ID = 10 PGS0810
867. 160 SZ = AF(ID)«X**BF(ID) PGS0820
868. GO TO 180 PGS0830
869. 170 SZ = 5000. PGS0840
870. GO TO 190 PGS0850
871. 180 IF (SZ.GT.5000.) SZ = 5000. PGS0860
872. 190 SY = 465.116«XY»SIN(TH)/COS(TH) PGS0870
873. C 465.116 = 1000. (M/KM) / 2.15 PGS0880
874. RETURN PGS0890
875. END PGS0900
98 6-84
-------
876.
877 .
878.
879.
880.
881.
882.
883.
884.
885.
886.
887.
888.
889 .
890.
891.
892.
893.
894.
895.
896.
897.
898.
899 .
900.
901 .
902.
903.
904.
905.
906.
907 .
908 .
909 .
910.
911.
912.
913.
914.
915.
916.
917.
918.
919 .
920.
921 .
922.
923.
924.
925.
926.
927 .
928.
929.
930.
931.
932.
933.
934.
935.
936.
937.
938.
939.
940.
941.
942.
943.
944.
945.
946.
947.
948.
949.
950.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
40
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
50
C
C
C
C
C
C
60
C
C
C
SUBROUTINE PLMRS(KST,U,TEMP,PL,IOPT,HANE,DELH,H) PLR0010
COMMON /STA/ QP,HPP,TSP,VSP,DP,VFP,SYOP,SZOP,DH PLR0020
VS = VSP PLR0030
TS = TSP PLR0040
D = DP PLR0050
MODIFY WIND SPEED BY POWER LAW PROFILE IN ORDER TO TAKE INTO PLR0060
ACCOUNT THE INCREASE OF WIND SPEED WITH HEIGHT. PLR0070
ASSUME WIND MEASUREMENTS ARE REPRESENTATIVE FOR HEIGHT=HANE. PLR0080
THT IS THE PHYSICAL STACK HEIGHT PLR0090
THT = HPP PLR0100
POINT SOURCE HEIGHT NOT ALLOWED TO BE LESS THAN 1 METER. PLR0110
IF (THT.LT.l.) THT = 1. PLR0120
U - WIND SPEED AT HEIGHT 'HANE' PLR0130
PL - POWER FOR THE WIND PROFILE PLR0140
UPL - WIND AT THE PHYSICAL STACK HEIGHT PLR0150
UPL = U*(THT/HANE)»*PL PLR0160
WIND SPEED NOT ALLOWED TO BE LESS THAN 1 METER/SEC. PLR0170
IF (UPL.LT.l.) UPL = 1. PLR0180
BUOY = 2.45153«VS*D**2 PLR0190
TEMP- THE AMBIENT AIR TEMPERATURE FOR THIS HOUR PLR0200
DELT = TS - TEMP PLR0210
F = BUOY'DELT/TS PLR0220
CALCULATE H PRIME WHICH TAKES INTO ACCOUNT STACK DOWNWASH PLR0230
BRIGGS(1973) PAGE 4 PLR0240
HPRM = THT PLR0250
IF 1OPT=0, THEN NO STACK DOWNWASH COMPUTATION PLR0260
IF (IOPT.EQ.O) GO TO 40 PLR0270
DUM = VS/UPL PLR0280
IF (DUM.LT.1.5) HPRM = THT + 2.«D*(DUM-1.5) PLR0290
'HPRM' IS-BRIGGS' H-PRIME PLR0300
IF (HPRM.LT.O.) HPRM = 0. PLR0310
CONTINUE PLR0320
CALCULATE PLUME RISE AND ADD H PRIME TO OBTAIN EFFECTIVE PLR0330
STACK HEIGHT. PLR0340
PLUME RISE CALCULATION PLR0350
IF (KST.GT.5) GO TO 60 PLR0360
PLUME RISE FOR UNSTABLE CONDITIONS PLR0370
IF (TS.LT.TEMP) GO TO 70 PLR0380
IF (F.GE.55.) GO TO 50 PLR0390
DETERMINE DELTA-T FOR BUOYANCY-MOMENTUM CROSSOVER(F<55) PLR0400
FOUND BY EQUATING BRIGGS(1969) EQ 5.2, P 59 WITH PLR0410
COMBINATION OF BRIGGS(1971) EQUATIONS 6 AND 7, P 1031 PLR0420
FOR F<55. PLR0430
DTMB = 0.0297»TS»VS*»0.33333/D»«0.66667 PLR0440
IF (DELT.LT.DTMB) GO TO 70 PLR0450
DISTANCE OF FINAL BUOYANT RISE(0.049 IS 14*3.5/1000) PLR0460
BRIGGS(1971) EQUATION 7,F<55, AND DIST TO FINAL RISE IS PLR0470
3.5 XSTAR DISTF IN KILOMETERS PLR0480
DISTF = 0.049»F»*0.625 PLR0490
COMBINATION OF BRIGGSU971) EQUATIONS 6 AND 7, P 1031 FOR PLR0500
F<55. PLR0510
DELH = 21.425»F»»0.75/UPL PLR0520
GO TO 90 PLR0530
DETERMINE DELTA-T FOR BUOYANCY-MOMENTUM CROSSOVER(F>55) PLR0540
FOUND BY EQUATING BRIGGS(1969) EQ 5.2, P 59 WITH PLR0550
COMBINATION OF BRIGGSU971) EQUATIONS 6 AND 7, P 1031 PLROS60
FOR F>55. PLR0570
DTMB = 0.00575*TS»VS**0.66667/D«*0.33333 PUR0580
IF (DELT.LT.DTMB) GO TO 70 PLR0590
DISTANCE OF FINAL BUOYANT RISE (0.119 IS 34*3.5/1000) PLR0600
BR1GGSU971) EQUATION 7, F>55, AND DIST TO FINAL RISE PLR0610
IS 3.5 XSTAR. DISTF IN KILOMETERS PLR0620
DISTF = 0.119»F»«0.4 PLR0630
COMBINATION OF BRIGGS(1971) EQUATIONS 6 AND 7, P 1031 PLR0640
FOR F>55. PLR0650
DELH = 38.71«F»»0.6/UPL PLR0660
GO TO 90 PLR0670
PLUME RISE FOR STABLE CONDITIONS. PLR0680
DTHDZ = 0.02 PLR0690
IF (KST.GT.5) DTHDZ = 0.035 PLR0700
S = 9.80616*DTHDZ/TEMP PLR0710
IF (TS.LT.TEMP) GO TO 80 PLR0720
DETERMINE DELTA-T FOR BUOYANCY-MOMENTUM CROSSOVER(STABLE) PLR0730
FOUND BY EQUATING BRIGGSU975) EQ 59, PAGE 96 FOR STABLE PLR0740
BUOYANCY RISE WITH BRIGGSU969) EQ 4.28, PAGE 59 FOR PLR0750
99
6-84
-------
951. C STABLE MOMENTUM RISL. PLR0760
952. DTMB = 0.019582*TEMP*VS*SQRT(S) PLR0770
953. IF (DELT. LT. DTMB) GO TO 80 PLR0780
954. C STABLE BUOYANT RISE FOR WIND CONDITIONS.(Wl ND NOT ALLOWED PLR0790
955. C LOW ENOUGH TO REQUIRE STABLE RISE IN CALM CONDITIONS.) PLR0800
956. C BRIGGSU975) EQ 59, PAGE 96. PLR0810
957. DELH = 2.6*(F/(UPL»S))»*0.333333 PLR0820
958. C COMBINATION OF BRIGGS(1975) EQ 48 AND EQ 59. NOTE DISTF PLR0830
959. C IN KM. PLR0840
960. DISTF = 0.0020715«UPL/SQRT(S) PLR0850
961. GO TO 90 PLR0860
962. C UNSTABLE-NEUTRAL MOMENTUM RISE PLR0870
963. C BRIOGS(1969) EQUATION 5.2, PAGE 59 NOTE: MOST ACCURATE PLR0880
964. C WHEN VS/U>4; TENDS TO OVERESTIMATE RISE WHEN VS/U<4 PLR0890
965. C (SEE BR1GGS(1975) P 78, FIG 4.) PLR0900
966. 70 DELH = 3.*VS*D/UPL PLR0910
967. DISTF = 0. PLR0920
968. GO TO 90 PLR0930
969. C STABLE MOMENTUM RISE PLR0940
970. 80 DHA = 3.*VS*D/UPL PLR0950
971. C BRIGGSU969) EQUATION 4.28, PAGE 59 PLR0960
972. DELH = 1.5«(VS»VS*D*D*TEMP/(4.*TS*UPL))»»0.333333/S»»0.166667 PLR097C
973. IF (DHA. LT. DELH) DELH = DHA PLR0980
974. 90 H = HPRM + DELH PLR0990
975. RETURN PLR1000
976. END PLR1010
100
6-84
-------
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
977 .
978.
979 .
980.
981 .
982.
983.
984.
985.
986.
987 .
988 .
989 .
990.
991.
992.
993.
994.
995.
996.
997 .
998 .
999.
1000.
1001 .
1002.
1003 .
1004.
1005.
1006 .
1007 .
1008.
1009.
1010.
1011 .
1012.
1013.
1014.
1015.
1016.
1017.
1018.
1019.
1020.
1021.
1022.
1023.
1024.
1025.
1026.
1027.
1028.
1029.
1030.
1031.
1032.
1033.
1034.
1035.
1036.
1037.
1038.
1039.
1040.
1041 .
1042.
1043.
1044.
1045.
1046.
1047 .
1048.
1049.
1050.
1051.
C
C
C
C
10
20
C
C
30
C
C
C
C
C
C
40
50
C
C
C
C
C
C
C
C
C
C
60
C
SUBROUTINE PHOCES (NPUFF,I STEP.DISKEY,IT,1TKAV,SYMAX,IADT)
SUBROUTINE PROCES ASSIGNS KEYP VALUES TO EACH PUFF AND
DETERMINES SIGMA VALUES FOR EACH PUFF.
INTEGER DISKEY
DIMENSION FRMATU8)
COMMON /WEA/ WDIIU 1 44 ) ,VVSPD( 1 44 ) , HL( 1 44 ) , SGTH( 1 44 ) , SGPHt 1 44 ) ,
1PREP(144),XYOP(600),XZOP(600),IT I ME,KST(144),STANUM,SW(144),
2SV(144),UPLM(144)
COMMON /XP/ XPUFF(600),YPUFF(600),ZPUFF(600),ITPUFF1600),
1XDSP(600),QPUFF(600),SY(600),SZ(600),CONC(25),NREC,KEYP(600)
COMMON /ADVCT/ PDIR( 600 ) , PSPD( 600 ) , XSWC, YSWC, NUMX.NUMY , DX, DY, FRMAT
DT = I STEP
HL8 = 0.8«HL(IT)
IHL = HL(IT)
IF (IHLT.EQ.IHL) GO TO 20
IF CURRENT MIXED LAYER IS DIFFERENT THAN PREVIOUS
ONE, UPDATE KEYP(J)'S.
DO 10 1=1,NPUFF
IF (KEYP(I).EQ.4.OR.KEYP(1).EQ.6) GOTO 10
IF (KEYP(I).EQ.3.AND.ZPUFF(I).LT.IHL) KEYP(I) = 1
IF (KEYP(I).EQ.5.AND.ZPUFF(I).LT.IHL) KKYP(I) = 2
IF (KST(IT).GE.5) GO TO 10
IF (KEYP(I).EQ.3.OR.KEYP(I).EQ.5) GOTO 10
IF (SZ(I).GE.HL8) KEYP(I) = 4
CONTINUE
IHLT = IHL
CONTINUE
GO TO (30,130), DISKEY
DISKEY=1 FOR PG CURVES.
DISKEY=2 FOR TRAVEL TIME CURVES
CONTINUE
IF (IT.EQ.1) GO TO 50
IF (ITRAV.NE.ISTEP) GO TO 50
AT THE FIRST TIME STEP OF EACH NEW MET. PERIOD FIND THE
VIRTUAL DISTANCES FOR EACH PUFF.
DO 40 1=1,NPUFF
XDSP(I) = 0.
SET TRAVEL DISTANCED. AT THE BEGINNING OF NEW
MET. PERIOD.
IF (KEYP(I).EQ.6) GO TO 40
NO VIRTUAL DISTANCE COMPUTATION IS NECESSARY IF KEYP=6
KS = KST( IT)
IF (KEYP(I).EQ.3.OR.KEYP
-------
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1052.
1053 .
1054.
1055.
1056.
1057 .
1058.
1059.
1060.
1061.
1062.
1063.
1064.
1065.
1066.
1067.
1068.
1069.
1070.
1071.
1072.
1073.
1074.
1075.
1076.
1077.
1078.
1079.
1080.
1081.
1082.
1083.
1084.
1085.
1086.
1087.
1088.
1089.
1090.
1091.
1092.
1093.
1094.
1095.
1096.
1097.
1098.
1099.
1100.
1101.
1102.
1103.
1104.
1105.
1106.
1107.
1108.
1109.
1110.
1111.
1112.
1113.
1114.
1115.
1116.
1117.
1118.
1119.
1120.
1121.
1122.
1123.
1124.
1125.
1126.
70
C
80
C
90
C
100
C
110
C
120
C
130
C
140
150
C
GO TO 120
HL8) KEYP(J) = -1
HL8) SZ(J) = HL8
SYMAX.AND.KEYPl J) .EQ.4) KEYP(J) = 6
IF (KS.GE.5)
IF (SZ( J) .GE.
IF (SZ( J) .GE
IF (SY(J) .GE
GO TO 120
CONTINUE
KEYP=2
XDSP(J) = XDSPU) + X
XVZL - XDSP(J) + XZOP(J)
SYY = SY!J)
CALL LTSIG (SYY.DT.SYD)
SY(J) = SYD
CALL PGSIG ( XVZL, XVYL.KS, SYD ,SZD)
SZ(J) = SZD
IF (KS.GE.5) GO TO 120
IF (SZ( J) .GE.HL8) KEYP(J) = 6
IF (SZ(J) .GE.HL8) SZ(J) = HL8
GO TO 120
CONTINUE
KEYP=3
XDSP(J) = XDSP(J) + X
XVYL = XDSP(J) + XYOP(J)
XVZL = XDSP(J) * XZOP(J)
CALL PGSIG (XVZL, XVYL, 7, SYD, SZD)
SY(J) = SYD
SZ(J) = SZD
IF (SY(J).GE.SYMAX) KEYP(J) = 5
IF (SY(J).GE.SYMAX) SY(J) = SYMAX
GO TO 120
CONTINUE
KEYP=4
XDSP(J) = XDSP(J) + X
XVYL = XDSP(J) + XYOP(J)
CALL PGSIG (XVZL, XVYL.KS, SYD, SZD)
SY(J) = SYD
IF (SZ(J).LT.HLS) SZ(J) = HL8
IF (SY(J).GE. SYMAX) KEYP(J) = 6
IF (SY( J) .GE. SYMAX) SY(J) = SYMAX
GO TO 120
CONTINUE
KEYP=5
XDSP(J) = XDSP(J) + X
XVZL = XDSP(J) * XZOP(J)
CALL PGSIG (XVZL, XVYL, 7, SYD, SZD)
SZ(J) = SZD
SYY = SY(J)
CALL LTSIG (SYY, DT, SYD)
SY(J) = SYD
GO TO 120
CONTINUE
KEYP=6
XDSP(J) = XDSP(J) + X
SYY = SY(J)
CALL LTSIG (SYY, DT, SYD)
SY(J) - SYD
IF (SZU) .LT.HL8) SZ(J) = HL8
CONTINUE
END LOOP OVER PUFFS FOR PG CURVES
GO TO 230
CONTINUE
DISPERSION AS A FUNCTION OF TIME.
IF (IT.EQ.l) GO TO 150
IF (ITRAV.NE.ISTEP) GO TO ISO
DO 140 1=1,NPUFF
XDSP(I) = 0.
IF (KEYP( I) .EQ.6) GOTO 140
KS = KST( IT)
IF (KEYP( I) .EQ.3.0R.KEYP(I) .EQ.5) KS = 7
CALL VT1ME (KEYP( 1 ) , IT ,KS ,SY( I ) , SZ( I ) ,XYOP( I )
CONTINUE
CONTINUE
BEGIN LOOP OVER PUFFS FOR TRAVEL TIME CURVES
DO 220 J=1,NPUFF
KS = KST(IT)
TI = ITPUFF(NPUFF) - ITPUFF(J)
,XZOP( I ) )
PRO0760
PRO0770
PRO0780
PRO0790
PRO0800
PRO0810
PHO0820
PRO0830
PRO0840
PR00850
PRO0860
PRO0870
PRO0880
PRO0890
PRO0900
PRO0910
PRO0920
PRO0930
PRD0940
PRO0950
PRO0960
PRO0970
PRO0980
PRO0990
PRD1000
PRO1010
PRO1020
PRO1030
PRO1040
PRO1050
PRO1060
PRO1070
PR01080
PRO1090
PB01100
PRO1110
PRO1120
PRO1130
PRO1140
PRO1150
PR01160
PR01170
PRO1180
PRO1190
PRO1200
PR01210
PRO1220
PR01230
PRO12-40
PRD1250
PRO1260
PRO1270
PR01280
PRO1290
PRO1300
PRO1310
PRO1320
PRO1330
PRO1340
PR01350
PRO1360
PR01370
PRO1380
PRO1390
PR01400
PRO1410
PRO1420
PR01430
PRO1440
PRO1450
PRO1460
PRO1470
PRO1480
PR01490
PRO1500
102
6-84
-------
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1127 .
1128.
1129.
1130.
1131 .
1132.
1133.
1134.
1135.
1136.
1137 .
1138.
1139.
1140.
1141.
1142.
1143.
1144.
1145.
1146.
1147.
1148.
1149 .
1150.
1151 .
1152.
1153.
1154.
1155.
1156.
1157.
1158.
1159 .
1160.
1161.
1162.
1163.
1164.
1165.
1166.
1167 .
1168.
1169.
1170.
1171.
1172.
1173.
1174.
1175.
1176.
1177.
1178.
1179.
1180.
1181.
1182.
1183.
1184.
1185.
1186.
1187.
1188.
1189.
1190.
1191.
1192.
1193.
1194.
1195.
1196.
1197.
1198.
1199.
1200.
1201.
1202.
1203.
1204.
160
C
170
C
180
C
190
C
200
C
210
C
220
C
230
IF ( TI . LT . 0 . 5 ) GO TO 220
GOTO (160,170,180,190,200,210), KEYP(J)
CONTINUE
KEYP=1
XDSP(J) = XDSP(J) + I STEP
TVY = XDSP(J) + XYOP(J)
TVZ = XDSP(J) + XZOP(J)
CALL JSISIG (KEYP(J),IT,TVY,KS,SYD,SZD)
SY(J) = SYD
IF (SY(J).GE.SYMAX) KEYP(J) = 2
IF (SY(J).GE.SYMAX} SY(J) - SYMAX
CALL JSISIG (KEYPtJ),IT,TVZ,KS,SYD,SZD)
SZ(J) = SZD
IF (KS.GE. 5 ) GO TO 220
IF (SZ(J).GE.HL8) KEYP(J) = 4
IF (SZ(J).GE.HLS) SZ(J) = HL8
IF (SY(J).GE.SYMAX.AND.KEYP(J).EQ.4) KEYP(J) = 6
GO TO 220
CONTINUE
KEYP=2
XDSP(J) = XDSP(J) + I STEP
SYY = SY(J)
SZZ - SZ(J)
CALL LTSIG (SYY.DT.SYD)
SY(J) = SYD
TVZ = XDSP(J) + XZOP(J)
CALL JSISIG (KEYP(J),IT,TVZ,KS,SYD,SZD)
SZ(J) = SZD
IF (KS.GE.5) GO TO 220
IF (SZ(J).GE.HLB) KEYP(J) = 6
IF (SZ(J).GE.HL8) SZ(J) = HL8
GO TO 220
CONTINUE
KEYP=3
KS = 7
XDSP(J) = XDSP(J) + 1STEP
TVY = XDSP(J) + XYOP(J)
TVZ = XDSP(J) + XZOP(J)
CALL JSISIG (KEYP(J),IT,TVY,KS,SYD,SZD)
SY(J) = SYD
IF (SY(J).GE.SYMAX) KEYP(J) = 5
IF (SY(J).GE.SYMAX) SY(J) = SYMAX
CALL JSISIG (KEYPtJ),IT,TVZ,KS,SYD,SZD)
SZ(J) = SZD
GO TO 220
CONTINUE
KEYP=4
XDSP(J) = XDSP(J) + I STEP
TVY = XDSP(J) * XYOP(J)
CALL JSISIG (KEYP(J),IT,TVY,KS,SYD,SZD)
SY(J) = SYD
IF (SY(J).GE.SYMAX) KEYP(J) = 6
IF (SY(J).GE.SYMAX) SY(J) = SYMAX
IF (SZ(J).LT.HL8) SZ(J) = HL8
GO TO 220
CONTINUE
KEYP=5
KS = 7
XDSP(J) = XDSP(J) + ISTEP
SYY = SY(J)
TVZ = XDSP(J) + XZOP(J)
CALL LTSIG (SYY,DT,SYD)
SY(J) = SYD
CALL JSISIG (KEYP(J),IT,TVZ,KS,SYD,SZD)
SZ(J) = SZD
GO TO 220
CONTINUE
KEYP=6
XDSP(J) = XDSP(J) + ISTEP
SYY = SY(J)
CALL LTSIG (SYY,DT,SYD)
SY(J) = SYD
IF (SZ(J).LT.HL8) SZ(J) = HL8
CONTINUE
END LOOP OVER PUFFS FOR TRAVEL TIME CURVES
CONTINUE
RETURN
END
PRO1510
PRO1520
PRO1530
PRO1540
PRO1550
PRO1560
PRO1570
PRO1580
PRO1590
PRO1600
PRO1610
PRO1620
PRO1630
PRO1640
PRO1650
PRO1660
PRO1670
PRO1680
PRO1690
PRO1700
PRO1710
PR01720
PRO1730
PRO1740
PRO1750
PRO1760
PRO1770
PRO1780
PRO1790
PRO1800
PRO1810
PRO1820
PRO1830
PRO1840
PRO1850
PRO1860
PRO1870
PRO1880
PRO1890
PRO1900
PRO1910
PRO1920
PRO1930
PRO1940
PRO1950
PRO1960
PRO1970
PRO1980
PH01990
PRO2000
PRO2010
PRO2020
PRO2030
PRO2040
PRO2050
PRO2060
PRO2070
PRO2080
PRO2090
PRO2100
PRO2HO
PRO2120
PRO2130
PRO2140
PR02150
PRO2160
PRO2170
PR02180
PRO2190
PRO2200
PRO2210
PRO2220
PR02230
PH02240
PR02250
PRO2260
PRO2270
PRO2280
103
6-84
-------
1
1
1
2
2
3
3
3
3
3
3
3
2
2
1
1
1
1
1205.
1206.
1207.
1208.
1209.
1210.
1211.
1212.
1213.
1214.
1215.
1216.
1217.
1218.
1219.
1220.
1221.
1222.
1223.
1224.
1225.
1226.
1227.
1228.
C
C
140
145
150
155
200
SUBROUTINE VTI ME(KEY,IT,KS,SYO,SZO.XYO.XZO)
SUBROUTINE VTIME COMPUTES VIRTUAL TIMES. THIS PROCESS IS
ANALOGOUS TO COMPUTING THE VIRTUAL DISTANCE.
DO 200 IC=1,2
X1NC = 0.
KD
DO
= 10
150 K=1,KD
DX = 10.*«(KD-K)
DO 140 J=1,1U
XINC = XINC + DX
XINCR = XINC/1000.
CALL JSISIG(KEY,IT,XINCR,KS,SYD,SZD)
DIFF = SYD - SYO
IF (IC.EQ.2) DIFF = SZD - SZO
IF (DIFF) 140,155,145
CONTINUE
XINC = XINC - DX
CONTINUE
CONTINUE
IF (IC.EQ.1) XYO = XINCR
IF (IC.EQ.2) XZO = XINCR
CONTINUE
RETURN
END
VT I 0 010
VTI0020
VTI0030
VTI0040
VTI0050
VTI006U
VT10 0 7 0
VTI0080
VT10090
VTI 0 1 0 0
VTI 011 0
Vri0120
VT I 0 1 3 U
VTI0140
VT 1 0 1 5 0
VTI 016 0
VTI 017 U
VTI 018 U
v r i o 19 o
VTI0200
VTI0210
VTI 0 2 2 0
VTI 0 2 3 0
VTI0240
1229.
1230.
1231.
1232.
1233.
1234.
1235.
1236.
1237.
1238.
1239 .
1240.
1241.
1242.
1243.
1244.
1245.
1246.
1247.
1248.
C
C
10
20
30
40
C
50
60
70
FUNCTION XVY (SYO.KST)
XVY CALCULATES THE VIRTUAL DISTANCE NECESSARY TO
ACCOUNT FOR THE INITIAL CROSSWIND DISPERTION.
GO TO (10,20,30,40,50,60,70), KST
XVY = (SYO/213.)*•!. 1148
RETURN
XVY = (SYO/155.)••!.097
RETURN
XVY = (SYO/103.)**!.092
RETURN
XVY = (SYO/68.)»*!.076
RETURN
D(DAY) AND D(NIGHT) .ARE THE SAME FOR SIGMA Y.
XVY = (SYO/68.)*•!.076
RETURN
XVY = (SYO/50.)»•!.086
RETURN
XVY = (SYO/33.5)**1.083
RETURN
END
XVY0010
XVY0020
XVY0030
XVY 0040
XVY0050
XVY0060
XVY0070
XVY0080
XVY0090
XVY0100
XVY0110
XVY0120
XVY0130
XVY0140
XVY0150
XVY0160
XVY017U
\\YO180
XVY0190
XVY0200
104
6-84
-------
1
1
1
1
1
1
1
1
1
1
1249.
1250.
1251.
1252.
1253.
1254.
1255.
1256.
1257.
1258.
1259.
1260.
1261.
1262.
1263.
1264.
1265.
1266.
1267 .
1268.
1269.
1270.
1271.
1272 .
1273.
1274.
1275.
1276.
1277.
1278 .
1279.
1280.
1281.
1282.
1283.
1284.
1285.
1286.
1287.
1288.
1289.
1290.
1291.
1292.
1293.
1294.
1295.
1296.
1297 .
1298.
1299.
1300.
1301.
1302.
1303.
1304.
1305.
1306.
1307 .
1308 .
1309 .
1310.
1311.
1312.
1313.
1314.
1315.
C
C
C
10
20
30
C
40
50
60
C
70
C
75
C
80
90
100
C
110
120
130
C
140
150
160
FUNCTION XVZ (SZO.KST) XVZ0010
XVZ CALCULATES THE VIRTUAL DISTANCE NECESSARY XVZ0020
TO ACCOUNT FOR THE INITIAL VERTICAL DISPERTION. XVZ0030
DIMENSION SA(7), SB(2), SD(5), SE(8), SF(9), AA(8), AB(3), AD(6), XVZ0040
IAE(9), AFUO), CA(8), CB(3), CD(6), CEO), CF(10) XVZ0050
DATA SA 713.95,21.40,29.3,37.67,47.44,71.16,104.65/ XVZ0060
DATA SB 720.23,40.7 XVZ0070
DATA SD 712.09,32.09,65.12,134.9,251.27 XVZ0080
DATA SE 73.534,8.698,21.628,33.489,49.767,79.07,109.3,141.867 XVZ0090
DATA SF 74.093,10.93,13.953,21.627,26.976,40.,54.89,68.84,83.257 XVZ0100
DATA AA 7122.8,158.08,170.22,179.52,217.41,258.89,346.75,453.857 XVZ0110
DATA AB 790.673,98.483,109.3/ XVZ0120
DATA AD 734.459,32.093,32.093,33.504,36.650,44.0537 XVZ0130
DATA AE 724.26,23.331,21.628,21.628,22.534,24.703,26.97,35.42, XVZ0140
147.6187 XVZ0150
DATA AF 715.209,14.457,13.953,13.953,14.823,16.187,17.836,22.651, XVZ0160
127.074,34.2197 XVZ0170
DATA CA 71 . 0585 , .9486 ,. 9147 ,. 8879 ,. 7909, . 7095, . 5786, .47257 XVZ0180
DATA CB 71.073, 1.017,.91157 XVZ0190
DATA CD 71.1498,1.2336,1.5527,1.6533,1.7671,1.95397 XVZ0200
DATA CE 71.1953,1.2202,1.3217,1.5854,1.7497,1.9791,2.1407,2.6585, XVZ0210
13.37937 XVZ0220
DATA CF 71.2261,1.2754,1.4606,1.5816,1.8348,2.151,2.4092,3.0599, XVZ0230
13.6448,4.60497 XVZ0240
GOTO (10,40,70,75,80,110,140), KST XVZ0250
STABILITY A(10) XVZ0260
DO 20 ID=1,7 XVZ0270
IF (SZO.LE.SAtID)) GO TO 30 XVZ0280
CONTINUE XVZ0290
ID = 8 XVZ0300
XVZ = (SZO/AA(ID))*«CA(ID) XVZ0310
RETURN XVZ0320
STABILITY B (40) XVZ0330
DO 50 ID=1,2 XVZ0340
IF (SZO.LE.SBtID)) GO TO 60 XVZ0350
CONTINUE XVZ0360
ID = 3 XVZ0370
XVZ = (SZO/AB(ID))'»CB(ID) XVZ0380
RETURN XVZ0390
STABILITY C (70) XVZ0400
XVZ = CSZO/61.141)'»1.0933 XVZ0410
RETURN XVZ0420
STABILITY D(DAY) XVZ0430
XVZ = (SZO730.9057)»»1.2088 XVZ0440
RETURN XVZ0450
STABILITY D (80) XVZ0460
DO 90 ID=1,5 XVZ0470
IF (SZO.LE.SD(ID)) GO TO 100 XVZ0480
CONTINUE XVZ0490
ID = 6 XVZ0500
XVZ = (SZO/AD(ID))«»CD(ID) XVZ0510
RETURN XVZ0520
STABILITY E (110) XVZ0530
DO 120 ID=1,8 XVZ0540
IF (SZO.LE.SEtID)) GOTO 130 XVZ0550
CONTINUE XVZ0560
ID = 9 XVZ0570
XVZ = (SZO/AE(ID))«»CE(ID) XVZ0580
RETURN XVZ0590
STABILITY F(140) XVZ0600
DO 150 ID=1,9 XVZ0610
IF (SZO.LE.SF(ID)) GOTO 160 XVZ0620
CONTINUE XVZ0630
ID = 10 XVZ0640
XVZ = (SZO7AF(ID))'»CF(ID) XVZ0650
RETURN XVZ0660
END XVZ0670
105
6-84
-------
PLOT POSTPROCESSOR
The source code listing of the plotting routine is presented
in subsequent pages. The program consists of a main module and
two subroutines.
106
6-84
-------
1
1
1
1
2
1
1
1
1
1
1
1
1
1
I
1
I
1
1
1
2
1
1
1
1
1
1
1
1
1
1
1
34
35
3b
37
38
39
40
41
42
43
44
45
4b
47
48
44
50
5 1
52
53
54
5 5
56
5 7
58
59
60
bl
o2
63
64
b5
66
67
b8
69
70
7 1
72
73
74
75
3 .
4 .
5 .
6 .
7
8 .
9 .
10.
11 .
12.
13.
14.
15.
16 .
17.
18 .
19 .
20.
21 ,
22.
23.
24.
25.
2b .
2 7 .
28 .
29 .
JO .
3 1 .
J 2 .
J 3 .
34 .
35 .
Jb .
J7 .
38 .
39 .
40 .
41 .
42.
43 .
44 .
45 .
4b .
47 .
48 .
44 .
50 .
5 1 .
52 .
5 J .
54 .
5 5 .
56 .
5 7 .
58 .
59 .
60 .
bl .
o2.
63.
b4.
b5 .
66 .
67 .
b8 .
69.
70.
7 1 .
72.
73.
74.
75 .
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
10
20
30
40
50
C
bO
COMMON /PUFTRJ/ XPAGE , YPAGE ,XPUF( 600 ) , YPUF( 600 ) ,SRR( 600 ) ,XMIN,XMAX
1,YMIN,YMAX,AXL,AYL,XREC(25),YREC(25),NREC
COMMON /PUFCON/ XPG , YPG ,XS I , YS I .CONTU60U ) ,TIME( 3600 ) ,TIMIN
DIMENSION CONC(25), IREC125)
XSI=0.
READ (5,«) IPLT
IF
IF
IF
IPLT = 1 PLOT CONCENTRATION VS. TIME
IPLT = 2 PLOT PUFF TRAJECTORY
IPLT = 3 PLOT BOTH
IF (IPLT.EQ.l.OR.IPLT.EQ.3) READ (5,») IYR,NUMR,ITPT,XSI,YSI
IYR IS THE ORDERS OF MAGNITUDE TO BE PLOTTED ON THE
Y AXIS. IF IYR = 0 THE DEFAULT VALUE WILL BE SET AT
IYR = 6.
NUMR IS THE NUMBER OF RECEPTORS THAT
CONCENTRATION VS. TIME IS PLOTTED FOR.
ITPT IS THE NUMBER OF PERIODS YOU WANT PLOTTED TOGETHER.
ITPT MUST BE EVENLY DIVISABLE INTO NTIME. IF ITPT > 999
ALL PERIODS ARE PLOTTED TOGETHER.
XSI LENGTH OF X AXIS IN INCHES
YSI LENGTH OF Y AXIS IN INCHES
IF (IPLT.EQ.l.OR.IPLT.EQ.3) READ (5
IF (IPLT.EQ.2.OR.IPLT.EQ.3) READ (5
(IREC(I),1=1,NUMR)
XMIN,YMIN,XSIZE,YSIZE,AXL,A
1YL
XMIN EAST-WEST COORDINATE OF SW CORNER OF
PLOTTING GRID.
YMIN NORTH-SOUTH COORDINATE OF SW CORNER OF
PLOTTING GRID.
XSIZE EAST-WEST SIZE OF PLOTTING GRID (KM)
YSIZE NORTH-SOUTH SIZE OF PLOTTING GRID (KM)
\XL LENGTH OF X AXIS IN INCHES
\YL LENGTH OF Y AXIS IN INCHES
XMAX=XMIN+XSIZE
YMAX=YMIN*YSIZE
INITIALIZE PLOT AND DEFINE ORIGIN.
CALL PLOTS (I,I,14)
CALL PLOT (3. ,2. ,-3)
XPAGE=0.
YPAGE=0.
XPG=0.
YPG=0.
IF ( IPLT. NE. LAND. IPLT. NE . 3) GO TO 90
DO 80 IRP=1,NUMR
CMAX=0.
IVAL=0
READ (22) NTIME,ITIME.NREC
DO 10 1=1,NREC
READ (22) XREC(I),YREC(I),Z
XRR=XREC(IREC(IRP))
YRR=YREC(IREC(IRP))
READ (22) IT,I STEP
TIMIN=((IT-1)'ITIME)/60.
IR=MOD(IT,ITPT)
IF (IR.EQ.O) TIMIN=((IT-ITPT)•ITIME)/60.
IF (ITPT.GE.999) TIMIN=0.
READ (22) ITM,(CONC(K),K=1,NREC)
IVAL=IVAL+1
TIME(IVAL)=ITM
CONT(IVAL)=CONC(IREC(IRP))
IF (CMAX.LT.CONT(IVAL)) CMAX=CONT(IVAL)
IF ( ITM.LT.IT*I TIME) GO TO 30
READ (22) NPUFF
DO 40 NP=1,NPUFF
READ (22) D1,D2,D3,D4,D5
IF IITPT.GT.999) GO TO SO
IF (IR.NE.O) GO TO 60
CONTINUE
XM=XSI + 2.
DEFINE NEW ORIGIN FOR NEXT PLOT.
CALL PLOT (XM.O.,-3)
CALL PLTCON (IVAL,IYR.XRR,YRR.CMAX)
IVAL=0
IF (IT.LT.NTIME) GO TO 20
GO TO 70
CONTINUE
(KM)
(KM)
PLT0010
PLT0020
PLT0030
PLT0040
PLT0050
PLT0060
PLT0070
PLT0080
PLT0090
PLT0100
PLT0110
PLT0120
PLT0130
PLT0140
PLT0150
PLT0160
PLT0170
PLT0180
PLT0190
PLT0200
PLT0210
PLT0220
PLT0230
PLT0240
PLT0250
PLT0260
PLT0270
PLT0280
PLT0290
PLT0300
PLT0310
PLT0320
PLT0330
PLT0340
PLT0350
PLT0360
PLT0370
PLT0380
PLT0390
PLT0400
PLT0410
PLT0420
PLT0430
PLT0440
PLT0450
PLT0460
PLT0470
PLT0480
PLT0490
PLT0500
PLT0510
PLT0520
PLT0530
PLT0540
PLT0550
PLT0560
PLT0570
PLT0580
PLT0590
PLT0600
PLT0610
PLT0620
PLT0630
PLT0640
PLT0650
PLT0660
PLT0670
PLT0680
PLT0690
PLT0700
PLT0710
PLT0720
PLT0730
PLT0740
PLT0750
107
6-84
-------
1
1
1
1
1
1
1
1
1
76 .
77 .
78 .
79.
80.
81.
82.
83.
84.
85.
86.
37.
88.
89.
90 .
91.
92.
93.
94.
95.
96.
97 .
98.
99.
100.
101 .
102.
103 .
104.
105 .
106 .
107 .
108 .
109 .
110.
111.
112.
70
80
90
100
110
120
130
MO
C
150
160
C
170
180
C
IF { IT.LT.NTIME) GO TO 20
IF (IT.EQ.NTIME.AND.ITPT.GT.999) GO TO 50
WRITE (6,170)
CONTINUE
REWIND 22
CONTINUE
CONTINUE
IF (IPLT.EQ.l) GO TO 160
IF IIPLT.EQ.2.OR.IPLT.EQ.3) GO TO 100
GO TO 150
CONTINUE
READ (22) NTIME,ITIME.NREC
DO 110 I=1,NREC
READ (22) XREC(I),YREC(I),Z
READ (22) IT.ISTEP
READ (22) ITM,(CONC(K),K=1,NREC)
IF (ITM.LT.ITMTIME) GO TO 130
READ (22) NPUFF
DO 140 NP=1,NPUFF
READ (22) XPUF(NP),YPUF(NP),Z,SY,SZ
SRR(NP)=SY/1000.
XM=AXL+2
IF (XSI.GT.AXL) XM=XSI+2.
DEFINE NEW ORIGIN FOR NEXT PLOT.
CALL PLOT (XM,0.,-3)
CALL PLTTRJ (NPUFF)
IF (IT.LT.NTIME) GO TO 120
GO TO 160
WRITE (6, 180) IPLT
CONTINUE
CALL PLOT (10.,0.,3)
CALL PLOT (10. ,0. ,999 )
STOP
FORMAT (5X,'NTIME IS NOT AN EVEN MULTIPLE OF ITPT',/)
FORMAT (5X,'IPLT= ',12,' IPLT MUST HAVE A VALUE OF 1,2,OR 3',/)
END
PLT0760
PLT0770
PLT0780
PLT0790
PLT0800
PLT08IO
PLT0820
PLT0830
PLT0840
PLT0850
PLT0860
PLT0870
PLT0880
PLT0890
PLT0900
PLT0910
PLT0920
PLT0930
PLT0940
PLT0950
PLT0960
PLT0970
PLT0980
PLT0990
PLT1000
PLT1010
PLT1020
PLT1030
PLT1040
PLT1050
PLT1060
PLT1070
PLT1080
PLT1090
PLT1100
PLT1UO
PLT1120
PLT1130
108
6-84
-------
114. SUBROUTINE PLTCONf IVAL , IYR.XR, YR.CMAX)
115. COMMON /PUFCON/ XPG,YPG,XSI,YSI,CONT(3600),TIME(3600),TIMIN
116. DIMENSION XRAYt3600),CON(3600)
117. IFUYR.LE.O) IYR=6
118. DO 10 I=1,IVAL
1 119. CT=CONT(I)
1 120. IF(CT.LT.1.OE-35) CT=1.0E-35
1 121. CON(I)=CT
1 122. 10 XRAY(I)=TIME(I)/60.
1 123. C XC IS SCALE FOR X AXIS.
124. XC=(XRAY(IVAL)-TIMIN)/XSl
125. DO 20 1=1,70
1 126. IP=-(35-I)
1 127. UPVAL=CMAX/10.»»IP
1 128. IF(UPVAL.LT.L) GO TO 30
1 129. 20 CONTINUE
130. 30 BMAX=10.«*IP
131. BMIN=BMAX/10.«»IYR
132. DO 40 1=1,IVAL
1 133. 40 IF(CON(I).LT.BMIN) CON(I)=BMIN
134. YC=IYR/YSI
135. C YC IS SCALE FOR Y AXIS.
136. C YSI LENGTH OF Y AXIS
137. C PLOT Y AXIS AT XPG AND XSI.
138. CALL LGAXS(XPG,YPG,22HCONCENTRATION (G/M»»3),22,YSI,90.,
139. 1BMIN.YC)
!40. CALL LGAXStXSI , YPG , 22HCONCENTRATION (G/M"3 ) , - 22 , YSI , 90 . ,
141. 1BMIN.YC)
142. C PLOT X AXIS AT YPG AND YSI.
14-5- CALL AXIS(XPG,YPG, 14HTIME (MINUTES ) , - 1 4 ,XS I , 0 . ,TIMIN ,XC)
'44. CALL AXIS(XPG,YSI,1H ,1,XSI,0.,TIMIN,XC)
145. XRAY(IVAL+l)=TtMIN
146. XRAY(IVAL+2)=XC
147. CON(IVAL+1)=BMIN
148. CON(IVAL+2)=YC
149. C PLOT DATA.
15U. CALL LGLINCXRAY,CON,IVAL,1,0,0,1)
151. XP=XPG+1
152. YP=YPG+YSI*0.5
153. C WRITE TITLE AT TOP OF PLOT.
154. CALL SYMBOL(XP.YP,0.14,9HRECEPTOR ,0.,9)
155. CALL NUMBER(999.0,YP,0.14,XR,0.,3)
156. CALL SYMBOL(999.0,YP,0.14,3H , ,0.,3)
157. CALL NUMBER(999.0,YP,0.14,YR,0.,3)
158. RETURN
159. END
PLC0010
PLC0020
PLC0030
PLC0040
PLCD050
PLC0060
PLC0070
PLC0080
PLC0090
PLC0100
PLC0110
PLC0120
P LCD 130
PLC0140
PLCD 150
P LCD 160
PLC0170
PLC0180
PLCD 190
PLC0200
PLC0210
PLC0220
PLC0230
PLC0240
PLC0250
PLC0260
PLC0270
PLC0280
PLC0290
PLC0300
PLC0310
PLC0320
PLC0330
PLC0340
PLC0350
PLC0360
PLC0370
PLC0380
PLC0390
PLC0400
PLC0410
PLC0420
PLC0430
PLC0440
PLC0450
PLC0460
109
6-84
-------
160. SUBROUTINE PLTTRJ(NPUFF) PLJ0010
161. COMMON /PUFTRJ/ XPAGE,YPAGE,XPUF(600),YPUF(600),SRR(600), PLJ0020
162. IXMIN,XMAX,YMIN,YMAX,AXL,AYL,XREC(25),YREC(25),NREC PLJ0030
163. DIMENSION XP(602),YP(602) PLJ0040
164. DO 10 1=1,NPUFF PLJ0050
1 165. XP(I)=XPUF(I)/1000. PLJ0060
1 166. 10 YP('I )=YPUF( D/1000. PLJ0070
167. XC=(XMAX-XMIN)/AXL PLJ0080
168. YC=(YMAX-YMIN)/AYL PLJ0090
169. C XC IS SCALE FOR X AXIS. PLJ0100
170. C YC IS SCALE FOR Y AXIS. PLJ0110
171. XP(NPUFF+1)=XMIN PLJ0120
172. XP(NPUFF+2)=XC PLJ0130
173. YP(NPUFF+1)=YMIN PLJ0140
174. YP(NPUFF+2)=YC PLJ0150
175. C PLOT X AXIS AT YPAGE AND YPP. PLJ0160
176. CALL AXIS(XPAGE,YPAGE,26HX(KM) EAST-WEST COORDINATE,-26,AXL, PLJ0170
177. 10.,XMIN,XC) PLJ0180
178. YPP=YPAGE+AYL PLJ0190
179. XPP=XPAGE+AXL PLJ0200
180. CALL AXIS(XPAGE,YPP,26HX(KM) EAST-WEST COORDINATE,26,AXL, PLJ0210
181. 10.,XMIN,XC) PLJ0220
182. C PLOT Y AXIS AT XPAGE AND XPP. PLJ0230
183. CALL AXIS(XPAGE,YPAGE,28HY(KM) NORTH-SOUTH COORDINATE,28,AYL, PLJ0240
184. 190.,YMIN,YC) PLJ0250
185. CALL AX1SUPP,YPAGE,28HY(KM) NORTH-SOUTH COORDINATE,-28,AYL, PLJ0260
186. 190.,YMIN,YC) PLJ0270
187. DO 20 1=1,NPUFF PLJ0280
1 188. X=(XP(I)+SRR(I)-XMIN)/XC PLJ0290
1 189. Y=(YP(I)-YMIN)/YC PLJ0300
1 190. SR=SRR(I)/XC PLJ0310
1 191. C DRAW CIRCLE AROUND PUFF CENTER, RADIUS= 1 SIGMA. PLJ0320
L 192. CALL CIRCL(X,Y,0.,360.,SR,SR,0.) PLJ0330
i 19J. 20 CONTINUE PLJ0340
194. DO 30 I=1,NREC PLJ0350
1 195. XI=I PLJ0360
1 196. XR=(XREC(I)-XMIN)/XC PLJ0370
1 197. YR=(YREC(I)-YMIN)/YC PLJ0380
1 198. C PLOT RECEPTOR LOCATIONS. PLJ0390
1 199. CALL NUMBER(XR,YR,0.14,XI,0.,-1) PLJ0400
1 200. 30 CONTINUE PLJ0410
1 201. C PLOT LINE THROUGH PUFF CENTERS. PLJ0420
202. CALL LINEUP,YP,NPUFF, 1,0,0) PLJ0430
203. RETURN PLJ0440
204. END PLJ0450
110
6-84
-------
Date
Chief, Environmental Operations Branch
Meteorology and Assessment Division (MD-80)
U. S. Environmental Protection Agency
Research Triangle Paark, NC 27711
I would like to receive future revisions to the
"INPUFF User's Guide."
Name
Organizat ion
Address
City
State Zip Code
Phone (Optional) ( )
6-84
------- |