EPA-600/4-75-005-a
May 1975
Environmental Mnnitoi mi)
DEVELOPMENT
OF AN URBAN
AIR QUALITY SIMULATION MODEL
WITH COMPATIBLE RAPS DATA
VOLUME I
:'f,
IW.
*V „-?
U.S. Environmental Protection Agency
Office of Research and Development
Washington, D. C. 20460
-------
EPA-600/4-75-005-3
DEVELOPMENT
OF AN URBAN
AIR QUALITY SIMULATION MODEL
WITH COMPATIBLE RAPS DATA
VOLUME I
by
C.C. Shir and L.J. Shieh
IBM Research Laboratory
San Jose, California 95193
Contract No. 68-02-1833
ROAP No. 26AAI23
Program Element No. 1AA003
EPA Project Officer: Robert E. Eskridge
Chemistry and Physics Laboratory
Office of Research and Development
Research Triangle Park, N. C. 27711
Prepared for
U.S. ENVIRONMENTAL PROTECTION AGENCY
Office of Research and Development
Washington, D. C. 20460
May 1975
-------
EPA REVIEW NOTICE
This report has been reviewed by the National Environmental Research
Center - Research Triangle Park, Office of Research and Development,
EPA, and approved for publication. Approval does not signify tlrat the
contents necessarily reflect the views and policies of the Environmental
Protection Agency, nor does mention of trade names or commercial
products constitute endorsement or recommendation for use.
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S . Environ-
mental Protection Agency, have been grouped into series. These broad
categories were established to facilitate further development and applica-
tion of environmental technology. Elimination of traditional grouping was
consciously planned to foster technology transfer and maximum interface
in related fields. These series are:
1. ENVIRONMENTAL HEALTH EFFECTS RESEARCH
2 . ENVIRONMENTAL PROTECTION TECHNOLOGY
3. ECOLOGICAL RESEARCH
4. ENVIRONMENTAL MONITORING
5. SOCIOECONOMIC ENVIRONMENTAL STUDIES
6. SCIENTIFIC AND TECHNICAL ASSESSMENT REPORTS
9. MISCELLANEOUS
This report has been assigned to the ENVIRONMENTAL MONITORING
series. This series describes research conducted to develop new or
improved methods and instrumentation for the identification and quanti-
fication of environmental pollutants at the lowest conceivably significant
concentrations. It also includes studies to determine the ambient concen-
I rations of pollutants in the environment and/or the variance of pollutants
as a function of time or meteorological factors.
This document is available to the public for sale through the National
Technical Information Service, Springfield, Virginia 22161.
Publication No. EPA~600/4-75-005-a
11
-------
DEVELOPMENT OF URBAN AIR QUALITY SIMULATION
MODEL WITH COMPATIBLE RAPS DATA
by
C. C. Shir and L. J. Shieh
IBM Reseach Laboratory, San Jose, California 95193
IBM Scientific Center, Palo Alto, California 94304
Final Report
Contract No, 68-02-1833
Prepared for
Meteorology Laboratory
0. S. Environmental Protection Agency
Research Triangle Park, N.C, 27711
May 1975
-------
PAGE 2
Abstract
An advanced generalized urban air quality model (IBMAQ-2) is
developed based on the theory utilized in an existing model
(IBMAQ-1) as prescribed in Ref, 1. The model, based on
numerical integration of the concentration equation,
computes temporal and three-dimensional spatial
concentration distributions resulting from specified urban
point and area sources by using NEDS {National Emission Data
System) and simulated RAMS (Regional Air Monitoring System)
data. The urban surface roughness values, estimated from
measured mean building heights and calculated building
density, are treated explicity as input parameters. The UTM
(Universal Transverse Metric) coordinates are used in all
geographical, source emission, and monitoring data. A new
method to incorporate point sources into the grid
computation is developed by using a Lagrange trajectory
method. Many model options are provided which enable users
to study conveniently the significant effects which these
options have on the final concentration distributions.
These options offer the user not only flexibility in the
adoption of data, but also offer one timely optimization of
the model.
The program description is included to provide a guide for
users. The program is constructed in a modular form which
allows users to change or improve each component
conveniently. The input auxiliary model, which processes
geographical, source emission, and monitoring data, is also
included.
The model was tested using simulated meteorological
conditions for a three-day period. The significant effects
of urban surface roughness en concentration distribution
were also investigated.
-------
PAGE 3
Acknowledgements
The authors wish to acknowledge the helpful discussions of
Dr. H. Eskridge, Dr. K. Demerjian, and Mr, K. L. Calder of
the Meteorology Laboratory, U. S. Environmental Protection
Agency, The appreciation is extended to Dr. R, H. Kay, Dr.
B. H. Armstrong, and Dr. E.. Clementi for their encouragement
and assistance during the course of this study. Special
thanks are due to Dr. D. Fox, who was the first government
project officer, for his guidance in the early part of the
study, and Ms. G. Smiley of the Meteorology Laboratory,
Environmental Protection Agency, for her assistance in
testing the program on the EPA computer facility. Our
gratitude is also due Dr. B. Armstrong and Dr. P. K.
Halpern for their proofreading of the manuscript.
This study was supported by the 0. S. Environmental
Protection Agency under Contract No. 68-02-1833.
-------
PAGE 4
°f Contents
Part I: Model Description
1. Introduction 8
2. Model Formulation 10
A. Simulation Region and Grid System.................... 10
Bi Equations, Boundary Conditions and
Initial Conditions...................................11
C. Parameterization-.................................... 15
D. Method of Analysis 35
E. Modeling Options.. 48
3. Discussion................................................51
Part II: Program Description
1. Introduction..............................................54
2. Program IBMAQ-2 — Diffusion Computation.................. 56
A. System fieguirements.................................. 56
B. Program Structure.............................. ...58
C. Suaaary of Subroutines and Their Functions........... 63
D. Variables Dsed in the Program. ....98
3. Computational Procedure...................................109
A. Parameter Specification.............................. 111
B. Input Specification 6 Requirement.. 116
C. Output Specification & Requirement..... ..,,....,121
D. Results Storage... 124
E. JCL Specifications... ..................,,~.......127
-------
PAGE 5
4. Auxiliary Program — Input Data Preparation 130
A. Area Source Data..................................... 136
B« Point Source Data... 137
C, Surface Parameter Dat.a.. ............................. 139
D. BAHS Station Location Data 140
E. Validity Check of All the Data , 141
F. Function of Each Auxiliary Program.,.. 141
References
-------
PAGE 6
of Figures
Figure 1-1. St. Louis metropolitan area and horizontal
numerical grid specification
Figure II-1. General Flow Chart of IBMAQ-2
Figure II-2- Detailed Flow Diagram of BAIN program
Figure II-3. Detailed Flow Diagram of SUBROUTINE AACOHP
Figure Il-4. Schematic Diagram of Processing NEDS
and Geographical Data
-------
PAGE 7
Table II-1. Subroutines Included in the Program IBMAQ-2
Table II-2. Arguments Passed in Each Subroutine
Table II-3.
DIMENSION, COMMON and EQDIVALENCE
statements of Main Program.
Table II-H I/O Units Used in the Program
Table II-5. Example of NAMELIST/INLIST/Input Data
Table II-6. Output Operations Controlled by LWRITE (N)
Table II-7. Sample of JCL, Example 1
Table II-8. Sample of JCL, Example 2,
Table II-9. Sample of Area Source Emission Data
Table 11-10. Sample of NEDS Point Source Data
-------
PAGE 8
I. HODEL DESCRIPTION
1. Introduction
An urban air quality simulation model can be used to study the
complex relation between the source emissions and ambient
concentrations, and the dependence of this relation on the
meteorological and urban surface conditions. The development of
such a model depends on current knowledge of urban air pollutant
transport, on the diffusion mechanism, and on the availability
and quality of source emission, meteorological, and surface
condition data. In addition, the purpose and computation
requirements of the model should be taken into consideration.
Therefore, the degree of sophistication of a model should be
commensurate with data availability, computational facility, and
application purposes. A brief summary of the development of
various modeling approaches and the objectives of developing a
grid model are discussed in Ref, 1, which is included in the
appendix.
He are concerned herein with a research-oriented model. Its
purpose is not only to understand air pollutant transport
phenomena, but also to study various significant parameters which
affect the urban air quality. Although the model utilizes
currently available data, through study of its results one can
infer additional characteristics or types of data that nay be
useful and therefore should be collected in the future. In
-------
PAGE 9
addition, the model can be used to study and verify various
assumptions used in simpler model formulations. Therefore, it
can serve as a guide for the future development of an improved
simpler model.
In this work, an advanced urban air quality simulation model is
developed based on the numerical integration of the concentration
equation over a specified set of grid elements. The model
offers a large degree of flexibility by taking into account the
temporal and spatial variation of meteorological variables and
surface conditions. These parameters are not normally
incorporated in most simplified models. Moreover, the model is
designed to be able to accept both extensive and limited amounts
of data. It is especially developed to adopt BAPS and RAMS data
collected in the St. Louis metropolitan area for evaluation and
validation. However, during the course of model development, the
data acquisition system of these programs had not yet been
completed. Therefore, the model development could not benefit
from the proposed extensive data collection, and further
improvement or modification may be needed when this data becomes
available. The model is applied to determine the SO
distribution in St. Louis. It also can be extended to study
multiple-species distributions without extensive modification.
-------
PAGE 10
2. Model Fornulation
2-A. Simulation Region and Grid System
The region of interest is the St. Louis metropolitan area as
defined in the RAPS project. The factors which are
considered, are as follows: spatial distributions of point
and area sources, area sources inventory resolution,
monitoring station distribution, computational grid
resolution and computational requirements. Most of the area
sources are concentrated in a region of 20 by 20 km in the
St. Louis area. The point sources are distributed along the
Mississippi river, extending from the Alton area in the
north to 60 km south downstream. A region of 60 km
north-south by 40 kin east-west centered at HAMS station 101
is sufficient to cover the major point and area sources as
well as most of the RAMS stations. The individual area
source inventory grid elements in the St. Louis urban area
are 1 k& square. Using a uniform 1 km grid system, a 40 x
60 network with a total of 2400 grid elements is required
for the region. This number is too large for our
computation. In view of this, a "stretch" grid system is
designed to obtain optimum resolution and economical
computations. This system is composed of 30 grid elements
east-west and 40 grid elements north-south. In additions,
there are finer grid elements 1 km square at the center
urban region, and coarser grid elements of 1 km by 2 km
-------
PAGE 11
rectangular size and 2 km x 2 km size on the outer suburban
region. The total volume in the region and under the mixing
height is to be simulated. There are four remote RAMS
stations outside the region which may provide some
information about the background concentrations. There are
also three significant point sources outside the region, and
their effects on the concentrations inside the region
require a further investigation.
2-B. Equations, Boundary Conditions and Initial Conditions
Since this area has a reasonably flat terrain, the effects
of topography are neglected. The surface roughness
parameter z is used to represent the effects of urban
buildings, forests, and rolling hills, S02 emission is
considered passive and therefore does not alter the
meteorological conditions. The turbulent diffusion of SO
is assumed to be of the gradient diffusion type. The
governing equation of SO^ in the atmosphere based on the
mass conservation law is
where C is the mean concentration of SO , V = (U, V, W) is
4-r —V
the mean wind vector, Q is the source strength rate, R is
the chemical reaction rate, K is the horizontal eddy
-------
PAGE 12
2
diffusivity, VTT is tne horizontal Laplacian operator, and K
H
is the vertical eddy diffusivity. The upper and lower
boundary conditions are:
K ~ = 0 , z = 0, H (2)
d Z ;
where H = H{t), the mixing depth of the planetary boundary
layer, is assumed constant in space due to the fact that the
data for a spatially varying H are difficult to obtain and
not likely to be available in the near future. The
absorption of the SO by the ground surface is neglected
because it requires data describing the surface properties
and their absorption rate. The significance of the
spatially varying mixing height and surface absorption rate
must be investigated, so that their priority in a data base
can be established. As these data become available, their
incorporation into the model is a simple matter. The inflow
lateral boundary conditions are:
c = cb , if v • 5 < o (3)
-------
PAGE 13
where C, is the background concentration, and ^ is the unit
outward vector nornal to the lateral surface. The outflow
lateral boundary conditions are:
0=° . »-<>. • W
and
Where x and y are the east and north boundries of
max J max
the area respectively. The data required to specify the
exact lateral boundary conditions are not available. These
conditions in (4), (5), which extrapolate the concentration
outside the region, serve as a reasonable approximation when
the region of conputation is large enough. In practice, and
as we have found, the lack of well-posed boundary
conditions (due to the absence of data required for the
proper conditions) does not cause serious problems. This is
because the horizontal advection terms, which dominate the
horizontal diffusion terms, are only first order in the
-------
PAGE 14
space derivative. The computation may be affected by the
inflow boundary conditions. However, if there is no high
concentration outside the inflow boundary this influence is
minimal, and the linear extrapolation could offer a fair
approximation. In this study, the computational region is
about ten times larger than the urban area in which major
sources are located. Moreover, there is no nearby major
urban area to influence the inflow boundary conditions,
although there are three major point sources located to the
west and south of this region. Their influences may not be
negligible under certain meteorological conditions. The
effects of these remote point sources and return flow cannot
be treated by this model alone. One approach would be an
experimental study to determine the background concentration
in order to evaluate the necessity of developing a larger
regional model to estimate the background concentration.
Since the background concentrations may be from elevated
sources, the surface concentration measurement alone may not
be sufficient to describe their three-dimensional spatial
distributions.
The initial conditions of concentration are arbitrarily set
to zero. The reguired three-dimensional initial
concentration distributions cannot be generated easily from
the surface measurements alone. The model reguires the
initialization of the surface concentration distributions
based on the observed data from the monitoring stations as
-------
PAGE 15
well as the vertical distributions of concentration. These
input data are not usually available. However, the
computations show that the concentrations reach observed
levels approximately within two hours under average wind
speed conditions. Hence, the initial conditions are not
iaportant for the concentration computations after this
initial time. This required tiae period depends on the wind
speed and the region of interest. This tiae interval can be
estiaated by T^L/u, here L is the distance downwind froa
the sources and u is average wind speed. However, this may
not be valid when the wind changes its direction drastically
during this time period. In this case, the initial period
aay need to be longer. However, it is better, if possible,
to select an initial period in which the wind directions do
not vary too much. The aodel siaulation can begin at an
arbitrary time. However, it is simplest to begin at a time
which is consistent with the time period used for
acquisition of the input data (e.g., if the air quality
measurements are on a 24 hourly basis beginning 1400 LST,
then the model computations should be started at 1400 LST).
In addition, in order to compare the computed and observed
results, one should have compatible average time periods of
computations and field measurments.
2-C. Parameterization
The parameters required for the integration of the equation
-------
PAGE 16
(2) are V, K , K, Q, R, and H, which, except for Q, are not
provided explicitly or sufficiently. The following
t
discussion outlines the methods which were used to obtain
the required input data for the model.
It must be pointed out that the choice of parameterization
is a very delicate matter and may be somewhat controversial.
The parameterization of a model depends on the availability
and the accuracy of input data, the best available theory,
and computer resources available. Since none of these
factors is static, but all are continuously changing and
improving, one must leave some leeway in the model for
future improvement. In addition, the model can also be used
to study the significance and ramifications of various
atmospheric variables available from future measurements,
Thus, more parameterization options than are essential
should be included in order to allow users to study their
potential effect, to optimize the model, and to point the
way for future data acquisition. The following methods of
parameterization represent the best efforts of the authors
based on available data, current theory, and computer
capability. We emphasize that many options of the
parameterization have been included in the model for the
purposes of optimization, sensitivity testing, directing
future data acquisition, and future model improvement. When
all SAPS and RAMS data become available, one may have to
re-examine these paraaeterizations.
-------
PAGE 17
2-C-a) lind .field
U) Hori^ggtal wind distribution
The wind vector V = V (x, y, z, t) is required at every grid
element for each time step of integration. The
time-averaged surface wind field for the total region was
obtained by using a weighted interpolation scheme. Data
from the measurement stations were interpolated to a square
grid. The choice of the grid size for this analysis grid
system requires an empirical study. In general, it depends
on the spatial distribution and density of the observation
network of wind vector and the spatial variability of the
wind field within the region of the study. Currently, the
grid size of U KM is used. Several schemes have been tested
in this analysis. It was found that reasonable results can
be obtained from the equations^
Z (V/rmn2) Z (
«± . = ^ and 7^=5^
m,n
where u.. and v- • are components of the wind vector at
analysis grid points in x and y directions, respectively;
-------
PAGE 18
^* *w
u and v are initial quess field values at analysis grid
points and r is the distance from grid point point (i,j)
inn
to grid point {m,n). The initial-guess field is obtained by
assuming
u = u,
mn k
v = v,
mn k
for minimum of r.
where u, and v are wind vector components measured at
station k and r is the distance from a grid point to
station k.
From this analyzed wind field a linear interpolation is
employed to obtain a wind vector at each grid point used in
the numerical scheme. The interpolation formula for the u
component of the wind vector is
where 0 < a = 6x/D < 1; 0 < b = 6y/D < 1. 6x and 6y are
-------
PAGE 19
the distances between numerical grid point (k,l) and wind
analysis grid point (irj) in * an(* 7 direction
respectively. D is the grid size of wind analysis grid.
Similar expression is used for V component of the wind
vector.
In summary, the procedure for computing the wind vector at
each grid point used in the numerical scheme is as follows:
» A square grid system (wind field analysis grid) is
selected which is U km in size.
» An initial-guess wind vector {u,v) is obtained at each
grid point (m,n) from nearest available observation
station.
• The weighted interpolation formula is applied to the
initial-guess field to obtain the wind vector at each
analysis grid point.
• Simple linear interpolation is used to obtain the wind
vector at each grid point used in the numerical scheme
from the wind vector values at the four nearest
adjacent analysis grid points.
The analysis of a mesoscale wind field is a very complex
problem. In reality, we do not expect the simple method
-------
PAGE 20
employed here will be able to handle all possible
meteorological conditions. Due to the unavailability of
RAMS data, whether a more complex method is necessary has
yet to be deterained.
The model also includes an option to use a subjective wind
field analysis. This is included in order to compare the
reliability of various wind-analysis schemes, It also
includes two other options; using a uniform wind field, and
generating the wind field directly from RAMS data. The
former could be used to study the significance of
non-uniform distribution of wind, and the latter is included
for future daily operation of the model.
Vertical distribution of wind
The spatial distribution of upper layer wind data was not
available, and the present knowledge of urban meteorology
cannot offer much information about this distribution.
Fortunately, the concentration at those locations where the
local sources dominated are not sensitive to the upper wind.
The vertical wind profiles at each grid location are assumed
to be of power-law form:
= |V I (z/z )p , (6)
'+S1 ' S
-------
PAGE 21
where V and V are the upper and surface wind velocities at
->• ->s
the height z and z respectively. The power constant p is
s
determined by
p = £n (|¥l/IXl>/ An(z/Z) , (7)
where V and V are the winds at the height z_ and z,. The
-*•-*•! 31
values of p are restricted between 0.15 and 0.65 which are
the usually observed values (Ref. 2). The direction of the
upper wind is also unknown. It is understood that the upper
wind usually has a direction to the right of the surface
wind in the north4rn hemisphere. However, this angle cannot
be determined quantitatively from any theory under general
conditions. Two methods have been previously tested. One
method assumes no directional change and the other assumes
equal angle change with height from z to z . Results from
5 -J
both methods differ very little. Since the upper air
sounding data for BAPS were not available, it did not seem
profitable to develop any new scheme to use these
unavailable data.
-------
PAGE 22
(iii) Yertica. 1 wind components
The vertical winds are important under certain
meteorological conditions such as strong convergence due to
the urban heat island effect or to frontal passage.
However, these data are usually not available. Thus, the
vertical winds were calculated from horizontal winds through
the continuity equation. In the previous study, the
interpolated horizontal winds are so smooth that the
vertical winds are too siaall to influence the concentration
distributions significantly. In addition, the incomplete
knowledge of the vertical winds resulting from the urban
heat island effects make the simulation of the vertical wind
difficult. Moreover, when the upper air sounding data are
not recorded at a fixed location, substantial pre-analysis
is required before they can be used for the model. The
effects of vertical wind on the urban atmospheric transport
remain to be investigated.
2-C-b) Atmospheric stability
It is well known that turbulent diffusion of air pollutants
is the second most important mechanism next to the wind
transport in affecting air pollutant dispersion. Its
magnitude is directly related to the atmospheric turbulence
intensity. However, due to the complex and heteorogeneous
distribution of the atmospheric turbulence, it is not simple
-------
PAGE 23
to describe or categorize the turbulence intensity. The
aost significant factors influencing the turbulence are
wind, atmospheric stability, and surface conditions. The
surface condition, parameterized here as surface roughness,
has been especially neglected in most of the previous air
pollution diffusion theories and air pollution models.
However, it is the fundamental mechanical property
influencing the atmospheric flows. A practical method to
categorize the turbulence by considering the three
components, namely, wind, stability, and surface roughness,
has been developed and described in detail in Ref. 1. The
method first entails the calculation of the gross modified
Pasguill's stability class, which is a continuous function,
and then estimates the Monin-Obukhov length by the following
expressions:
1/L = ± [d ' Hn (1.2 + 10/ZQ)]2 ' 10f(S) , (8)
and
f(S) = -a/(l+b|s|) , (9)
-------
PAGE 24
where S is the stability class; a=4, b=1.3, c=.85,
d=0.216586, and z = ZQ (x,y) is the surface roughness
parameter in meters. The stability parameter S=0 denotes
neutral conditions, and -negative and positive values of S
denote unstable and stable conditions respectively- Thus,
the sign of L in equation (8) must be the same as that of
S, Moreover, the continuous values of S as discussed
previously allow smooth changes between stability classes.
Since this method which requires minimal data of wind,
insolation, and surface roughness, is much easier to use
than measuring vertical temperature distributions, its
merits should not be overlooked. Certainly, it may require
further improvements such as continuous insolation
classifications, non-uniform spatial stability
distributions, and effects of precipitation. The accuracy
of this method has to be evaluated using available data,
Future RAMS data will provide the vertical temperature
difference at 12 stations. It is not yet clear how to more
usefully utilize these data to describe the spatial
distribution of turbulence than the method used here because
these data are not yet available.
From the previous study, we did not find any systematic
errors (the difference between the computed and observed
concentration) with respect to the stability parameter.
Hence, no obvious errors from this method of determining
-------
PAGE 25
stability can be determined.
2-C-c) Eddj diffusivity
The treataent of turbulent diffusion as gradient diffusion
is a classical method. In view of the scarcity of the
relevent meteorological data needed to describe the
heteorogeneous turbulence distributions over an urban area
and test any turbulence theory, it is doubtful that a more
sophisticated method is necessary. The method to determine
the eddy diffusivity developed here is based on experimental
data and computed results from the turbulent transport
model which is described in Bef. 1, and can be expressed as
follows:
-b z/H
K = u*A/h , a = kQze ° , (10)
where u^, § , k and H are friction velocity,
noa-dimensional temperature gradient, von Karman constant,
and the height of the boundary layer, respectively. The
factor b , which should be a function of stability, has a
value of U determined by the turbulent transport model (Ref.
3). Its dependence on the stability requires further
investigation. From given L and Vg (wind velocity at z. -
10m), we can calculate Ks from the following procedures:
-------
Calculate u. from V , L and z by
* ->-s o
PAGE 26
»* = kol¥sl/V
(11)
where
z
z
(12)
$ is the non-dimensional wind shear. According to
Businger, et al (Bef.U) , 4> can be expressed as follows:
rm
1 +
(1 -
-1/4
where
-------
PAGE 27
= z/L, a = 4.7, and 3 = 15.
By use of last equations, equation
becomes
(12) after integration
(1 + z/z ) +
for
Z/ZQ) -
[(1
-1
+ 2 tan w - ir/2 , for C<0, (13)
where
= l/4>
2. Calculate K from u and A, by
7f T f->
K =
(14) '.
-------
PAGE 28
where
Y + a£
yd - 3"?)
and
Y = 0.74, 3' = 9
Formulas (12) to (14) are valid in the surface layer. Se
calculate K = K(z=10m) and then extrapolate K to higher
s
altitude by K = Ks £/£ . This approach by no means implies
s
that the turbulent diffusion of air pollutants over an urban
area is really so simple. These formulas, derived from
experimental data, are based on the assumption of
equilibrium for the turbulence. This assumption may not
hold over an urban area where the horizontal inhomogeneity
plays aa important role, because turbulence is not in
equilibrium in the vicinity of a change in surface roughness
and temperature. However, little information is available
about the urban effects on the eddy diffusivity. The
formulas used here are subject to improvement when more
knowledge is available. However, we feel that this method
is a significant improvement over existing methods because
the effects of the surface roughness are taken into account
-------
PAGE 29
in the calculation of I and K. The effects of the surface
roughness which are neglected in most other approaches, can
influence the concentration distribution significantly.
The horizontal eddy diffusivity, K^, has less effect on the
n
surface concentration distributions than K, the vertical
eddy diffusivity. However, it does significantly affect the
concentration from point sources and may not be neglected.
As the numerical scheme for the advection is significantly
improved, the artificial diffusion from finite difference
approximation error is reduced. This will require the
horizontal eddy diffusivity to be estimated aore accurately.
Both the horizontal and the vertical eddy diffusivity depend
on the height, wind, stability, and surface roughness.
Moreover, K also depends on the meso-scale wind motion
which may be prominent under stable conditions. The
meso-scale wind influences are too complex to describe in a
simple formula. As a preliminary formulation, we have found
that the horizontal eddy diffusivity can be related to the
vertical eddy diffusivity K by the relation:
KTT = f. (z,S) K
where to a first approximation, f is dependent on height
-------
PAGE 30
and atmosphere stability conditions. Although it is not
possible to determine the form of f analytically, its
behavior and values as an empirical representation can be
infered in principle from experimental data or from
calculation with an appropriate turbulent transport model.
Recently our study indicated that the value of f varies
slowly with height and stability conditions, and is about 10
near neutral conditions toward the top of the surface layer.
However, further study is required to obtain more precise
results and a more adequate functional representation of
f (Z,5). Currently. f = 10 is assumed in the model
k K
computations.
2-C-d) Surface roughness
The surface roughness, which is the fundamental mechanical
surface property affecting the atmospheric flow, is
neglected in most air guality models. There nay be several
reasons for this. In the past, both Gaussian Plume
diffusion theory and the empirical determination of the
diffusion parameters (a , a ) do not explicitly include the
z y
effects of surface roughness. Recently, Pasguill (Ref, 5)
included the effects of surface roughness on the a and a •
x Y
Another reason for omitting the effect of surface roughness
is related to the lack of such data over an urban area.
Thus, surface roughness effects may be treated implicitly
-------
PAGE 31
through a "tune up" process rather than calculated
explicitly. The most widely used "tune up" method for
including the surface roughness in a Gaussian plume model is
to modify the initial values of 0- according to the average
z
building height cf the area source grid elements.
The following simple formula was used to estimate the
roughness length of the urban area (Ref. 6):
ZQ = 0.5 rh , (15)
where r is the silhouette area ratio and h is the effective
mean height of the roughness elements. Data for the mean
building heights and mean surface altitudes are available
for the St. Louis area. Those of natural roughness
elements, such as forest, and the silhouette ratio are not
available. The silhouette ratio, which is the density of
the roughness elements, has values ranging from about 0.05
to 0,5 from rural to urban area, respectively. One crude
method to estimate the building density is to assume that
the building density is proportional to the annual area
source emission from space heating. This enables one to
estimate the urban surface roughness approximately without
extensive survey data. Although this method is not totally
-------
PAGE 32
satisfactory, it constitutes a plausible effort to include
the effects of surface roughness when sufficient data are
not available.
The formulas used to estimate the silhoutte ratio and the
effective mean height of the roughness elements are the
following:
r = 0.04 (1 + Q3(x,y)/Q ) , (16)
ct d
= h, +0.2h, +5
b t
where Q and Q are the annual mean areas source strength
a a J
and its spatial average value respectively; h, and b. are
the mean building height and terrain height respectively.
Inclusion of the terrain height is based on the assumption
that the rolling hills in the south-west area also
contribute to the turbulent mixing. The validity of
equation (16) requires further investigation. The constant
5m term in the formula for h represents our assumption of
the contribution of the natural roughness elements to the
effective mean height. This term leads to a minimum value
of surface roughness of 10 cm from the formula for z . This
o
method will be used until the extensive survey data become
-------
PAGE 33
available. The effects of surface roughness on the
distribution of the surface concentration are significant.
The magnitude of influence depends on the atmospheric
stability condition and heights of emitting sources. In St.
Louis area, our model computations indicate that a
reasonable results of the surface concentration values can
still be achieved by the use of a less complicated method
for estimating surface roughness. However, it is important
that the estimated values of surface roughness should
represent the characteristics of the surface obstacles in
the simulation region.
2-C-e) T.he mixing height
The mixing height is an important parameter to determine the
volume in which air pollutants can be diluted by turbulent
mixing. This mixing height usually varies temporally and
spatially, and is difficult to predict accurately by either
surface measurements alone or planetary boundary layer
theory. However, when the source-receptor distance is not
too large such as in St. Louis, and the mixing height is
greater than 1000 m, it does not have profound effects on
the surface concentration in the central urban area because
the time of advection is shorter than the time of diffusion.
In general, the determination of the mixing height requires
at least two measurements: one each for the daily maximum
and minimum height. In this study, since no data were
-------
PAGE 34
available, the average values of previously measured data
are used. In modeling the mixing height, the air pollutants
can be totally confined within the mixing height during its
diurnal variation. This approach seems to be reasonable in
the morning during the growth period of mixing height. On
the other hand, in the evening period, the temperature lapse
rate usually approaches a neutral profile and then a stable
lapse rate starts to grow near the surface. This causes a
certain portion of the air pollutants to be left above the
nixing height. when this type of discontinuous mixing
height occurs, the air pollutants can no longer be assumed
totally confined underneath the mixing height.
2-C-f) Chemical Jtsaction rate
The chemical reaction rates of SO is expressed by R = -k c,
^ cl
where k is the reaction rate constant. The chemical
cl
reaction parameterization in the model could be considered
controversial due to its complexity. Its dependencies on
various factors such as insolation, water vapor, and other
pollutants are not yet completely clear. These chemical
reaction effects could be significant under low wind-speed
condition, when the time scale of a chemical reation is
comparable to the time scale of advection. A high value for
the chemical reaction rate, namely a half-life about two
hours, was used in the previous study. The results
-------
PAGE 35
nonetheless showed that the model constantly over-predicted
values for the lov wind-speed cases. This is not easily
explained by the choice of the chemical reaction rate alone.
It involves more complex factors such as the assumed
emission rate and emission height of the area sources. With
the current magnitude of the uncertainty in the area source
emission, this discrepancy requires further investigation of
both the chemical reactions and area source emission
inventory.
2-D. Methods of Analysis
2-D-a) JiEiii specifications
The three dimensional 30x40x14 grid system consists of
16800 grid elements (Fig. 1-1). The x, y, z axes are
oriented E-W, N-S and vertically, respectively.
The horizontal grid sizes are specified as follows:
Ax. = 1000 m , g <_ i <_ 25 ,
= 2000 m , otherwise ,
Ay. = 1000 m , 11 <. j ± 30 ,
= 2000 m , otherwise
-------
PAGE 36
XUTM=725
XUTM=765
YUTM
=4312
YUTM k
=4252
Ax=
1km
i i i l*~*l l i i i l i
Fig. 1-1 The St. Louis metropolitan area and horizontal
numerical grid specification.
-------
PAGE 37
This system, which may be termed a "stretched" grid system,
can cover a U0x60 square ka area, and encompass smaller grid
sub^elements compatible with the area source inventory grid
size at the central urban area. It also has other
advantages such as the program is invariant to grid size,
and excessive index computation can be avoided. It must be
noted that the horizontal coordinate system used here is the
Universal Transverse Mercator (UTM) systeii (Ref. 7). The
vertical grid sizes are specified as follows:
=
f 20m , 1 <_ k <_ 5 ;
25m , 6 <_ k <_ 9 ;
(H-200)/4 , 10 <_ k <_ 13, for H >_ 300m;
125m , 10 < k < 13, for H < 300m,
-------
PAGE 38
where fl is the nixing height when H > 300m. The lower nine
grid points under 200m have a fixed spatial size because the
effective heights of point and area sources are within this
layer. The levels of the upper four grid points are
determined by the hourly varying mixing height. The grid
spacing for these four grid points is set to be larger than
or equal to 25m. Thus, the minimum height of the grid
system is 300m. When the mixing height is lower than the
top of the grid system, the grid spacing remains constant.
However, small values of eddy diffusivity were forced at
those levels located in the inversion layer the mixing
height, where the mixing is small.
-------
PAGE 39
2-D-b) numerical Methods
& second-order central finite-difference scheme was used to
integrate the advection and horizontal diffusion term,
Since the concentration fields usually have large
variations, phase errors and damping resulting from the
finite difference method are large, and careful treatment is
required. Our finite difference approximation to the
concentration equation is
W[C*ijk]
VCijkJ + At Q.jk/AV - At
where AV = AxAyAz is the grid volume, and (1, V, W are
finite difference operators for horizontal (U,V) and
vertical (W) advection, and V and fl are analogous
A, £ 2
operators for horizontal and vertical diffusion. The time
step index is denoted by n, and i, j, Jt are grid indices
-------
PAGE 40
along the x, y, z directions respectively. The asterisk.
denotes the interia values when a time-splitting operation
which will be explained later has been performed. The
operators U, I/, and W are defined by
U[C. ., ] = F [C ., ] - F ' [C. ., ] ,
ijk x ijk x
Fy CCijk] - V[Cijk]'
[CijkJ -
where the F-factors, representing the amount of the
concentration flux on the boundary of the grid element, are
defined as follows:
= min (Ci_1/jkf Fx) , if axi > 0 ; (18a)
= min (Cijk , Fx) , if axi < 0 ,
Fx = min {Cijk ' Fx} ' if ax*i > ° ;
= min (Ci+1^jk , Fx) , if ax*± < 0 ;
-------
PAGE
-------
PAGE 42
Cljk = C
and
**
The derivation of finite difference scheme for horizontal
advection is given in Appendix 6, This method preserves the
amount of concentration when it is transported from one grid
element to the next. Because the method restricts the flux
advected from any grid elements to be larger than or equal
to the mass in that elements. In addition, it prevents the
concentration values from becoming negative due to phase
errors. For x direction advection, by comparing a Taylor
series expansion and the finite-difference formula used in
the model, it can be shown that the truncation error e for
an uniform grid is
= (1 - oT)Ax2 | UAt 9^C _ , 3N
8X~
-------
PAGE
where
• n 2C 3C \
N = min °' (l+a)Ax " ^T/
a = UAt/Ax
The detailed characteristics of this finite difference
operator are not easily analyzed. The reason for this lies
in the non-linearity which causes interaction between nodes
(
}
of different wavelengths. The first tern of the error is
that of the second-order central difference error. The
second tern is due to the flux correction described in
equation (18a) and (18b). One always can use a method with
higher accuracy at the expense of computing time. However,
this extra expense is meaningful only when it produces a
commensurate improvement. The sophistication of the
difference scheme should be compatible with the
sophistication of the remainder of the model and the
accuracy of the available data. He have attempted to
achieve an appropriate balance among these factors in this
model.
The operators P and V are defined by
Xy z
-------
PAGE 44
and
where
= At K[r . C. . .. - (1 + r . )C. ..
H V 1 "I •- 1 ~lK VT T~lV
J.1 A. J. _L_L^_JJV J^.J- -LJjS.
Ci+lf j
- (1+ryk)CiJk + Ci)/A] ' (20)
(2D
Yk = At/(Azk + Azk_1))
rk = Azk/AZk-l
and 0 is a parametric constant.
The stability criteria for the terms of advection,
horizontal diffusion, and chemical reaction rate are well
known (Ref. 8) as follows:
-------
PAGE 45
a , a < 1, where subscript a denotes x, y, or z;
a. a.
At Ku/(Ax2 + Ay2) '< 1/4, and At k < 1.
H ~~ <*
and
4(1-6) ' (22)
Condition (22) is used to determine the value of 9 in
equation (21). The value of 6 must be close to 1, if a
large value of y. is used. When the grid spacing is a
function of tine, a correction due to the changing volume is
needed, because equation (17) only applies to a constant
grid spacing system. It can be shown that -C. .. a£nAz, /at
should be added on the right hand side of equation (17).
The implementation of this correction is simplified by
changing the height H(t) for each hourly period. With this
arrangement, the concentrations at the upper five grid
points are adjusted hourly according to the height, H as
follows:
[c* "to./**. ' k - kn + 1? (23)
Cijk = Wl + r/s\ rn . _ .
J 8/ n n\ C._.v , k - kn ,
-------
PAGE
where
Sn = (H " zk )/(R " Zk }' and kn
n n
When a sudden decrease of the mixing height occurs, one
should use an alternative method, One such approach
provided as an option in the model avoids a compression of
the concentration by allowing some concentration to be left
above the nixing height. The same concentration profile is
therefore kept at the same height- This method is as
follows:
(24)
where
k' = k
n - sn
-------
PAGE U7
2-D-c) Point source formulation
One of the prevailing criticisms of the grid model is
related to the accuracy of the incorporation of the point
sources into the grid system. These comments are usually
qualitative rather than quantitative in nature. One
particular problem is related to initial diffusion. When
the concentration due to point source is injected into the
grid system, the concentration at a grid element represents
only the mean value over the grid element volume. The
sub-grid variation can only be accommodated fay a separate
sub-grid system which has not yet been developed. In this
study, an experimental method utilizing a Lagrange plume
trajectory is introduced. That is, one keeps track of the
location of the plume center and the plume width initially,
and incorporates the plume into the grid system when the
plume width is comparable with the grid size (see Section
2-E-e). In other words, the effective point source location
is calculated according to the wind field and the diffusion
speed. The method can be described as follows:
T
T = mir , ,
' 16 Jin (10) K.
H
-------
PAGE 46
where r , r , and v are the initial, and new point source
->o ->• .>
location, and the wind vector at the initial point source
location, respectively. The effective height of the point
source is h and K is the average vertical eddy diffusivity
below the stack height. The time scale of relocation is the
minimum of the time which the surface concentration requires
to reach its maximum, (first terra on the right-hand side of
equation 25) or the time required for the horizontal plume
area to become comparable with the grid element area {Second
term on the right-hand side of equation 25).
2-E, Modeling Options
The model developed does not consist of a single method, and
perhaps should be thought of as a family of methods, which
aim to describe the phenomena of air pollutant transport in
as realistic a manner as possible. This is achieved by
introducing a variety of options which offer the user the
ability to study the significance of the components as
parameters. It also makes the model more general in dealing
with input data, and offers a convenient method of
optimization. The optimization of the model here means to
select the most realistic and effective method to describe
the nature of air pollutant transport with available input
data. It is the most essential and time-consuming task
during the development of an air quality model. The
existence of these options makes the model a useful and
-------
PAGE U9
general tool to understand and describe the physical system
under investigation. The options included are as follows:
2-E-a) Choice of surface wind field
The surface wind field is the most important input to the
model other than the sources. Four options are available,
namely, a subjectively analyzed wind field, an objectively
analyzed wind field, data from RAMS stations, and a uniform
wind. The subjective wind field can provide a reference
with which to compare various objective wind analysis
schemes. The wind-field generating program which generates
the surface wind field over the simulated region by using
the observed wind data from the monitoring stations as
described in section (2-C-a i), is included, thus, it saves
one step of preanalysis of the wind field. The option of
uniform wind field can be used to test the effects of
non-uniform distributions of the wind field.
2-E-b) Choice of vertical variation of wijod field
One can construct the upper wind aligning it with the
surface wind, or can specify an angle change from the
surface wind, or uniform distribution of the upper wind.
The upper wind affects the high point sources more than
surface area sources, especially when the plume travels a
long distance.
-------
PAGE 50
2-E-c) Choice of vertical Hind components
The surface concentration could be significantly sensitive
to the vertical wind components. However, such data are
usually not available. Three options are included: no
vertical wind components, calcuating raw vertical wind
coaponents from the continuity equation, and smoothing the
calculated vertical wind components.
2-E-d) Choice of concentration adjustment .w_it]i mixing
varia.tj.on
The hourly variation of the mixing height influences the
concentration distributions. One of the options included is
to Jceep the total mass constant below the mixing height.
Thus, the concentrations below the mixing height are
compressed or diluted accordingly. Another option is to
keep the concentration distribution unchanged with respect
to altitude, while the grid size is changed,
2-E-e) Choice o_f £O_iiit sou_rce modeling
The options included are: to calculate the relocation of a
plume by using the local wind speed and appropriate time
scale of plume diffusion according to the specified formula,
equation (25), or simply specify the relocation of all
-------
PAGE 51
planes by a single tiae scale which is Just equal to one
tiae step used for numerical tiae integration. It is noted
that this tiae scale is used for calculating relocation—not
for integration. The latter is used when hypothetical eddy
diffusivity is eaployed. For instance, if one assumes the
urban area is extremely saooth with roughness length less
than 1 ca, then one vill have very saall values of eddy
diffusivity which, according to equation (25), aay lead to
an unusually large tiae scale. Other options included are:
to incorporate the relocated pluae into the grid element
nearest to the center of the pluae or distribute the plume
into the nearby four grid eleaents.
3. Discussions and Recommendations
k research-oriented urban air quality siaulation model such
as described herein is a powerful tool to understand the
phenomena of transport and diffusion of air pollutants over
an urban area. It is an atteapt to put together in a
practical balance as much as we know about air pollution and
urban aeteorology, ataospheric turbulence, source emission,
and aathematical technique, to describe the phenoaena, along
with a little ataospheric cheaistry. Moreover, the
availability and quality of the input data, and computer
capability should be considered siaultaneously. He hope the
aodel has been designed such that it is easily used and
-------
PAGE 52
modified by other users for their purposes.
Among all of these considerations, the most important item
is still the quality of the input data. A poor quality of
data introduces uncertainty or confusion, and makes
validation or evaluation of the model difficult.
In particular, the quality of the emission data, either
directly measured or derived, is the most important factor
in the process of development and validation of an urban air
pollution model. No sensible validation can be made without
good quality emission data. In addition, the quality of the
meteorological data, surface condition data, and air
pollution observation data are also important. However, the
required detail and accuracy of these data is not easily
determined. This must be investigated through use of the
nodel on these data.
On the other hand, some basic questions concerning the grid
•odel are: the adequacy of the eddy coefficient to describe
the turbulent flux, the representation of a concentration
field by a finite number of grid elements, the sub-grid
variation of concentration, especially near a point source
or monitoring station, and the numerical errors from the
difference scheme. These questions cannot be resolved
quantitatively without careful evaluation of the model with
data of good quality. The priority for improvement of these
-------
PAGE 53
limitations in the model oust also be deterained following a
thorough investigation with good data. The grid aodel is a
promising approach because of its inherent physical
validity, its generality and flexibility. We believe its
advantages as a research tool will becoae more and more
apparent as it is systematically iaproved, feature by
feature, in step with is proved data quality, and improved
knowledge of atmospheric structure.
-------
PAGE 5U
II. PROGRAM DESCBIPTION
1. Introduction
In part I of this report we described the physical
representation of our air quality simulation model (IBMAQ-2)
and presented the modeling raethologies, input requirements
and computational procedures. In this part we shall
consider the operational problem of using the computer
program based on these equations and assumptions to simulate
a concentration field in the St. Louis metropolation area.
The notation used in part I is retained in the following
discussion unless specified otherwise.
The computer program codes developed under this contract can
be divided into two parts. They are all written in FORTRAN
IV. The first part is the program for the main diffusion
computations (program IBMAQ-2), which consists of one main
program and 39 subroutines. The second part is the group of
auxiliary programs, which process and prepare input data for
the main diffusion model.
The program codes of IBMAQ-2 constitute a new version of our
earlier program IBMAQ-1 (Reference 1), The major changes in
-------
PAGE 55
this program are the following: a) implementation of
improved physical assumptions in the model formulation; b)
adaption of the program to operate with data to be obtained
from the RAPS project and RAMS stations in St. Loais; c)
reconstruction of i/o specifications; and d) streamlining
the general structure of the program so that it is easy to
modify for future work. It should be pointed out here, due
to the delays of the HAPS project, the program code of
IBBAQ-2 was completed without testing on the new data.
Thus, the I/O specification and certain methodologies
employed in this program may need to be revised when such a
data set is available.
In section II-2 and II-3, we shall discuss the program
IBMAQ-2. The auxiliary programs are briefly described in
section II-4. The program codes of IBMAQ-2 include quite
comprehensive comment statements, which are necessary for
understanding the program logic. The program listing, which
is given in Appendix 1, is a self-explanatory document. We
believe the user will have no difficulty in understanding
it. Therefore, in the following discussion we shall
concentrate on an over-view of the program and how to use
it.
-------
PAGE 56
Program IBMAQ-2. Diffusion Computations
2-A. System requirements
The model program is written in FORTRAN IV. The system
requirements depend upon the usage of the model and its
operational modes. The following three operational
configurations are used in our developaent and operation of
this model in a time sharing system.
2-A-a) Debugging runs:
64k words {one word equal to four bytes) of core are
required in this mode along with at least three sequential
files or data sets (disk space). The CPU (Central
Processing Unit) time is about half of the amount required
by a production run. This particular operational aode is to
check the execution of I/O statement in the program and to
detect any obvious programming errors (e.g. compatibility of
passing arguments from one routine to another). The
reduction of CPU core requirement is achieved by computing
concentration field only at five vertical layers and
assuming no vertical variation of wind field.
-------
PAGE 57
2-A-b) Testing runs:
75K words of core are required and at least three sequential
files or data sets. The CPU time is equivalent to
production runs. The concentration fields are computed at
14 vertical layers in the testing run. However, there is no
vertical variation of wind field. This operational mode is
to check the implementations of modeling methodologies in
the program codes.
2-A-c) Production runs;
96K words of core are required for production runs and at
least five sequential files or data sets. The CPU time is
2-5 minutes for each 24 hour real tine simulations on the
IBM 360/195. The precise running time depends upon the
meteorological conditions of the simulated days.
The line printer is required for all the above-listed
operation configurations. If the job is to be submitted in
a batch mode, a card reader is necessary.
The code was developed on an IBM System 360/195 running
OS/MVT and TSO. It has been tested on the EPA Research
-------
PAGE 58
Triangle Park computer facilities running on a time sharing
system. The code runs equally well in either environment.
2-B. General structure of the program
In this section, we shall briefly describe the overall
structure of the IBMAQ-2 program.
This model is based upon the numerical solution of a
concentration equation, and is, of course, quite different
from a Gaussian plume-type aodel. The pollutant
concentration is computed in three spatial dimensions in
successive time steps. The procedure is continued to obtain
the tia»e evolution of the concentration distribution until
the new values of source emissions and/or meteorological
variables are read in. in some instances, the steady state
solution may be reached before the next input time; in this
event it is useless to repeat the computation for the
subsequent time steps. In these cases, it is convenient to
pass directly to the new values of meteorological and/or
source data. Thus, the computational procedure has two
nested computational loops—a main loop and a time-step
loop. The former one reads in meteorological and source
emission data and computes the necessary parameters. The
-------
PAGE 59
latter one is nested in the Bain loop. this time-step loop
computes an instantaneous spatially varying concentration
field for each small time interval. This time interval is
dictated by the numerical stability of the finite-
difference scheme.
Fig. T-1,-1 shows the general flow diagram of program IBMAQ-2.
The program can be divided into three parts. They are as
follows:
2-B-a) Init ializatign:
This part of the program is executed only once for each run
of the program. The CPU storage allocation is first
executed, then the model parameters and options are read in,
constants are defined and variables and time indices are
initialized in this part of the program. It also inputs
geographical data of the modeled region and annual source
emission data. The names of the subroutines called in this
part of the program are: CDTOT_P, CONSIS, 6EQIN, XIUTHS.
lYOUrt. PRINTS, MITES, and HBITEX. These will be defined
in the sequel.
-------
PAGE 60
Specification statements
(DIMENSION; EQUIVALENCE;
NAME LIST; COMMON; DATA)
Input:
NAME LIST (model parameters)
Initializations:
Constants & variables
Input geographic & annual emission
data & edit data for output
Print out input data
L^
Initialize numerical grid system
i
r
Set up I/O devices for storage of
computed result & output
'
^
Input time varying meteorological
& source emission data
1
Compute time varying parameters
required by the model & output
'
r
Compute time step of numerical
integration
i
r
Compute concentration field for each
time step
I
Edit computed results for output
if required
^ Store model parameters & computed
"" results on disk or tape
~ Print out time varvinq input data &
*" model parameters
Figure 11-1. General Flow Chart of IBMAQ-2
-------
PAGE 61
2-B-b) Main loop:
This part of the program is executed every "LPBIHT" time
interval. This tiae interval is the smallest time interval
over which the hourly average source emission and
meteorological data are available and are to be input to the
program.
This main loop accomplishes the input of these data and
edits them for model use. The meteorological parameters
required by the model are also computed in this part of the
program. In additon, the variable vertical grid resolution
and the numerical time step are determined. The computed
average concentration field values at each grid point and
BANS stations are also transmitted as output to a peripheral
storage device for future usage. The names of the
subroutines used in this part of the program which will be
defined later are: OjJIA££f SOUSIN, 8IHPIH, HINDER.
VIM* MIF, STABIT, HHCJLLC. DIHEN1, SHIFTH, CAD JDS.
AKZCAJ,, PBI8TA, PBINfB, HRITEI. and DTTEST.
2-B-c) Time stejj loop:
This loop is nested in the main loop. The execution of this
part of the program is managed by subroutine AACOBP. In
-------
PAGE 62
this loop, the instantaneous concentration field is computed
every tiae step according to the finite-difference
approximation of the concentration equation. The hourly
average concentration field is also computed. Except for
subroutines PHI.NTC* STNCOH, SHIFTN, CCHECK, .WBITEX, HfllTEZ,
ZZGBID and TIMEJC, t&e subroutines used in this part of the
program are for the purpose of computing each term of
concentration equation. They are: OV£yiI» UliU, XjfDIFF,
^IDIFF, SQDBCE, and CHEHIC, and along with the previously
mentioned subroutines will be described in the next
section.
We believe that a program as complex as this model requires
should be written in a form so that it is easy for the user
to modify any part of it. This ease of modification can
extend to the improvement of physical parameterization of
the model, to the computational methodologies used, and the
I/O specification of the program. Hence, the current
prograa has been constructed in modular form to achieve this
objective. The execution of the program is controlled by
two routines, namely the MAIN program and the subroutine
AACOMP- The flow diagrams of these two routines are given
in section II-2-c.
The dimensional arrays used in the computation are all
-------
PAGE 63
passed through call arguments and their dimensions are
dynamically allocated in the subroutines. If one wishes to
change any subroutines or dimension statements, one only has
to deal with the calling routine and the referenced
subroutine.
Finally, all I/O operations are executed in only a few
subroutines. It is very easy to change I/O specifications
through the control parameters included in the HAHELIST
input. (See section II-3)
2-C. Summary of subroutines and their function
The computer program of this model consists of a main
program and a total of 39 subroutines. There are multiple
entries among five of these subroutines. For quick
reference. Table II-1 lists a brief description of each
subroutine. This listing is part of the comment statements
appearing in the main program (see Appendix 1). Table II-2
lists the variable arguments passing from call statements to
called routine or vice versa. Variables included in the
*.
common statements are not listed. (see Table II-3 for
variables in common.) Following are the detailed
descriptions of each subroutine:
-------
PAGE 64
MAIN PROGRAM: CALLED FROM: None
ROUTINE CALLED: AACOMP, AKZCAL, CADJUS, CONSIN
CDTOTP, DIHENS, DIMEN1, HHCALC,
DTTESX, GEOIN, OOTAPE, PRINTS
PBINTA, PRINTS, PRINTC, SHIFTN,
SOUSIN, TIMEX, WINDIN.
I/O PEHFOBBED: Yes
The main program calls various computational and I/O
subroutines. It sets up main computational loops to coapute
concentration fields. It is a program controlling a general
flow of execution of the IBMAQ-2 model. The detailed flow
diagram of this routine is shown in Fi£. XJ[^2.
-------
PAGE 65
NAME
AACOMP
AKZCAL
CADJUS
CCHECK
CHEM1C
COUSIN
*CDTOTP
DIMENS
*DIMEN1
*HHCALC
DTTEST
GEOIN
OUTAPE
POSITV
PRINTS
*PRINTA
*PRINTB
*PBINTC
SHIFTS
SOURCE
SOUSIN
STABIT
STNCON
TIMEX
UVZF
UVFLUX
UVINTP
WWZF
WINDER
iflNDGR
WINDIN
WRITES
*HRITEX
*WRITEZ
WWFLUX
XYDIFF
XYUTMS
*XYUTN1
ZZDIFF
ZZGRID
CALLED FROM
MAIN
HAIN
MAIN
AACOMP
AACONP
MAIN
RAIN
MAIN
MAIN
WAIN
HAIN
MAIN
MAIN
(NOT USED)
MAIN
HAIN
HAIN
MAIN
MAIN,AACOMP
AACOMP
MAIN
WINDIN
PHINTC
MAIN
HINDIN
AACOMP
WINDIN,WINDGR
WINDIN
WINDER,UVZF
WINDIN
MAIN
PRINTS
PRINTS
PRINTC
AACOMP
AACOMP
GEOIN
GEOIN,SOURCE
AACOMP
SOURCE
Z5JC1IOMS
COMPUTE THE CONCEHTBATION FIELDS
COMPUTE EDDI DIFFUSIVITY
ADJUST C VALUES DUE TO CHANGE IN GRID DIHEN,
CHECK FOR STEADY STATE CONDITION
COMPUTE CHEMICAL DECAY
SPECIFY MODEL PARAMETERS
PRINT CARD IMAGE OF NAMELIST INPUT
INITIALIZE GRID SYSTEM
SET VERTICAL GRID SYSTEM
COMPUTE TIME VARYING MIXING HEIGHT
TIME STEP FOR NUMERICAL METHOD
INPUT GEOGHAPH. AND ANNUAL EMISSION DATA
OUTPUT RESULTS TO TAPE OR DISK
SET C=0 IF IT IS LESS THAN ZERO
PRINT GEOGRAPH. AND ANNUAL EMISSION DATA
PRINT TIME VARYING EMISSION RATES
PRINT TIME VARYING METEOROLOGICAL DATA
PRINT CONCENTRATION FIELDS
SHIFT ARRAY A TO ARRAY fl
ADD NEW SOURCE EMISSION INTO THE SYSTEM
INPUT TIHE VARYING SOURCE EMISSION RATE
ESTIMATE CONTINUOUS STABILITY CLASSES
COMPUTE CONC. VALUES AT RAMS STATIONS
FIX TIHE INDICES
COMPUTE VERTICAL WIND PROFILE OF (U,V)
COMPUTE HORIZONTAL ADVECTION
INTERPOLATE ANALYZED U,V TO NUMERICAL GRID
COMPUTE i COMPONENT OF WIND
CONVERT WIND VECTOR TO COHP. OR VICE VERSA
GENERATE SFC. WIND FIELD FROM RAMS DATA
READ IN SURFACE WIND FIELD AND HAMS DATA
PRINT DATA ARRAY
PRINT HORIZONTAL FIELD OF ARRAY
PRINT VERTICAL CROSS-SECTION OF ARRAY
COMPUTE VERTICAL ADVECTION
COMPUTE HORIZONTAL DIFFUSION
COMPUTE UTH COORDINATE OF NUMERICAL GRID
CONVERT (X,Y) FROM UTM TO NUMERICAL GRID
COMPUTE VERTICAL DIFFUSION
CONVERT Z FROM METER TO NUMERICAL GRID UNIT
(NOTE: * DENOTE ENTRY POINT IN LAST STATED SUBROUTINE)
TABLE II-1 SUBROUTINES INCLUDED IN THE PROGRAM IBHAQ-2
-------
P^-E 66
SUBROUTINE AACOMP (CPl.C.CC.COLn.U.V.M.ZS.OA.P'V ,AKZ,AKN ./W.W
SUBROUTINE AKZCAL ( AKZ.U.V.ZO.AKF.AKH.Z'MMJ'ijJKM.KM)
SUBROUTINE. CADJUS (CPl.C.RDZ, IM\JM,K'^IMJ" J JKM.KTI.RH.
SUBROUTINE CCHECK (CPl,COLD,IJKM,IMJM,ICHECK,nCMri)
SUBROUTINE CHEMIC (CPl,C,IM,JM,KM,IM,ri,MKM,AKA)
SUBROUTINE CONSIN (I'1t,ritKM,Kr!tr\]'1fMKM,MK!l)
ENTRY CDTOTP ( JUNITl.IUNITl ,IUNIT2)
SUBROUTINE DIMENS (DX.DY.DZ .ROX.ROY.RnZ ,DXS .D^S .DZS ,Z ,Z'1,IMJ..n,KM)
ENTRY DIMEN1 ( nZ,RDZ,OZS,Z,ZM,K'1,KTI ,RH)
ENTRY MHCALC (HH.HMIN.HMAX)
SUBROUTINE DTTEST (U.V.AKH .AKZ.AKF.AKAHR.nZF.VZF.nX.HY ,Z,r%JM,K'1,TMn,
IJKN.KN)
SUBROUTINE OEOPI (OB ,ZS ,ZO,POD, ,XPMTM,YP!ITM,XP,YP,Zn ,Z^,Nn ,T!TI i,]MT"tnX
SUBROUTINE OUTAPE (CC.ICAL JOBS /ISX.PVJf.IWYTP JHnTP ,L''An")
SUBROUTINE POSITV (AtIMf.]f!tKMtimMtI,lK'l)
SUBROUTINE PRINTS ( CPl,CC,Cl,COLn,IORS ,IC.AL,ICVtt,Kr"L,rrv?S ,!J,V,",!'l
*,Vl,UU,W,AKZ,UZF,VZF,WZF,A!V
*,UZF,V7.F,AKF,AKH,AKZ,nX,nY,nZ,Z,EK,FK,TM,,l'\KM,Kf!J".rt,L")
SUBROUTINE SOUSIN (OB ,ZS ,ZO ,POB ,XP ,YP ,ZP,ZR,ri,,]r\LM)
Table II-2 ARHUMENTS USED IN E^CII SMBnOMTr!E
-------
.«nE 67
SUBROUTINE STABIT (V,ESKY,IHR,STAB,ITRAT)
SUBROUTINE STNCON (CC.IOBS .ICAL.KOBS.KCAL.ITOBS.NS .NSX.I'I.J'I.XS ,YS ,JXS ,
SUBROUTINE TIMEX (NMONDY)
SUBROUTINE IJVZF (U,V,IJ1,V1,UZF,VZF,WZF,Z,IM,J'1,H,!CI,LX)
SUBROUTINE UVFLUX (CPlfCtU,UZFtDX,nY,RDX,RDY,ZtEK,FK,I",JMtK'i,ri,ri
*,IJKM,LX,LY,MKN)
SUBROUTINE UVINTP (UU,VV,IN,JN,U,V,IMtJM,DX,DY,DELX)
SUBROUTINE IflZF (UfV,W,DXtDYtDZ,IM,JMfKM,KN,Uf'f,JIINIT)
SUBROUTINE MINDER (U.V.UV.THETA.I)
SUBROUTINE VJINDGR (Ul.Vl.lJU.VV.MEAR.USTH.VSTN.XSUT'l.YSUT^.IMn.JMT"
*,nX,DY,IM,JM,IN,JM,NS)
SUBROUTINE WINDIN (U,V, ''1,1)1, VI, UU,W, NEAR, USTN,VSTN,UZF,VZF,WZF,XS!!TM,YSI!TM,
*,IUTM,,]UTM,Z,DX,DY,DZ,ri,JM,K'1,IMJM,IJKN,IM,JN,KN
*,IS,NS)
SUBROUTINE '-/RITES (0,IXf1AX,IYf1AX,ISIZE,JSIZE,IBF^,JBEn,IUT'\JMTM,IFORM
-*,J UN IT .RATIO, TITLE)
ENTRY MRITEX (Q.IXMAX.IYMAX .IS IZE.JSIZE.IFORM.JtniT, RATIO .TITLE)
ENTRY HRITEZ ( A,Z,RATIO,IM,jn,KM,iriC,JMC,TITLA)
i
SUBROUTINE 1-MFLUX (CPl,C,M,MZF,Z,nX,DY,nZ,EK,FK,RnX,RnY,RnZ,IM,J"1,Kri,
*,IfWMtIJKH,I,JKN)
SUBROUTINE XYDIFF (CPl.C.AKH.AKZ.AKF.RDX.RDY.DXS.DYSJM.JM.n.IM.n.UKM)
SUBROUTINE XYIJTMS (IUTn,JUTM,IXBER,IYBEn,IBEn,JBEG,nX,DY,DXA,nY«,n,JM)
ENTRY XYUTM1 (XUP1,YIJTM,IDIMAX,Xn,YD,IDIM,L''!AY)
SUBROUTINE ZZDIFF (CPl,C,AKZ,AKF,RDZ,DZS,Z'1,EK,FK,n,JM,KM,IM,]M,IJK'l)
SUBROUTINE ZZRRID (ZA,ZB,Z,IMA,I^,KM)
TABLE II-2 (continued)
-------
P."^E 63
DIMENSION
* CP1(30,40,14),C(30,40,14),I)(30,40,7),V(30,40,7),''I(30,40,7)
* CP1(30,40,14),C(30,40,14),M(30,40,1),V(30,40,1),"(30,40,1)
* CP1(30,40,5),C(30,40,5),U(30,40,1),V(30,40,1),M(30,40,1)
* U(30,40,7),V(30,40,7),M(30,40,7)
* U(30,40,1),V(30,40,1), '-1(30,40,1)
* ,COLD(30,40), CC(30,40), Cl(30,40), 1)1(30,40), Vl(30,40)
* , OA(30,40), OB(30,40), ZS(30,40), ZO(3Q,40), AKZ(30,40)
* , UU 9,13), VV( 9,13),NEAR( 9,13)
* , DX(30), RDX(30), DXS(30), HY(40), RDV(40)
* , DZ(14), RDZ(14), DZS(14), Z(13), ZM(14),
* , AKF(14), UZF(14), VZF(14), MZF(14)
* , NP(150), XP(150), YP(150), ZP(150), ZR(150),
* ,PQA(150), PQB(150), EK(150), FK(150)
* , IS(25), XS(25), YS(25), JXS(25)
* ,1065(26), ICAL(26), KOBS(26), KCAL(26)
* ,IUTM(30),JUTM(40),XPUTM(150),YP!ITM(150)
* , AKA(24), NMONDY(12)
COMMON /AADATA/
* If1l,JMl,Kni,JUNIT,KUNITC,KU'1ITn,KUNIPT,K!!NITS,K!INIT'f
* ,IYR,IMO,IDAY,IHR,ITM,ITMNR,ITSEC,ITOTHR,ITSTEP,DT,TM,TSEC
* ,LPRINT,LTSTOP,LTSOUS,LT'fIND
* .LHRITE(IO) ,LSOIJS(2) ,LTOP,UITOP,LPO,U'H,LWIND,nINO,LCR!!N,LCHEM
* ,RAMS(6,25),PARM(10),Al(4),AK,HG,HP,HS,OLMi:j,DCMIN
* ,Pf1AX,PMIN,RIR,ZMAX,ZRPO,ZRISE,OBTOT,POBTOT,!)0,PHTFHZ,!!ZF
COMMON /CBLOCK/ CP1(30,40 ,14) ,C(30,40,14)
COMMON /CBLOCK/ CP1(30,40,5),C(30,40,5)
EQUIVALENCE
* ,( U(I,'I,'D! ''-1(1,1,1))
! (QA(I,'D;OB(I',I)); (pnA(i)l
, (UZF(.1),VZF(1)),
*
* ', (USTN(1),EK(120)), (VSTN(1),FK(120))
TARLE II-3 DIMENSION, COMMON, 1ND EOMprALEMCE STATEMENTS
OF MAIN PROGRAM
-------
Operation Sequence Subroutine Called Variables Defined and/or Comments
P1GE 69
'DIMENSION'; 'EQUIVALENCE';
'NAMELIST'; 'COMMON'
(See Table 11-2 for variable
values initialized)
(Copy NAMELIST input data to
JUNIT2 and print data in card
image format)
(See Table 11-3 for variable inputs
in this statement)
QA,ZS,Zo,PQA,XPUTM,YPUTM,
XP,YP,ZP,ZR,NP,IUTM,JUTM,
LM,XSUTM,YSUTM,XS,YS,IS,
JXS,JYS,NS,QBSUM,PQBSUM,
QBTOT.PQBTOT
CALL
DIMENS
DXS,RDX,DYS,RDY,DZS,RDZ,
Z,ZM,HH,PARM(5)
CALL
PRINTS
(Print geographic and annual
emission data)
***Main Loop Begin***
Figure 11-2. Detailed Flow Diagram of MAIN Program
-------
PAGE 70
Operation Sequence Subroutine Called Variables Defined and/or Comments
CALL
OUTAPE
Yes
(*Set up I/O UN
"Output compu
*Set CC=0
UNIT=KUNITP & KUNITC\
ted result I
(Title for print output)
CALL
SOUSIN
QB,PQB,ZR,QBTOT,PQBTOT,
PARM(6),PARM(7)
^-^
— ^(100)
(Print time varying source emission
data)
Figure 11-2. (continued)
-------
PAGE 71
Operation Sequence Subroutine Called Variables Defined and/or Comments
HH
CP1
CALL
1
WINDIN
UU,VV,U,V,W,RAMS,NEAR,
USTN,VSTN,UZF,VZF,WZF
PARM(1),PARM(2),PARM(3),
PARM(4),PARM(8),PARM(9)
Figure 11-2. (continued)
-------
Operation Sequence Subroutine Called Variables Defined and/or Comments
CALL
AKZCAL
RIB,UO,PHIFHZ,PARM(10)
AKZ.AKF
PAGE 72
CALL
DTTEST
DT
(Print time varying meteor, data
and model parameters)
""Time Step Loop Begin***
Figure 11-2. (continued)
-------
PAGE 73
SOBBOUTIMS AiCOMP: CALLED FROH: HAIg
BOOTIME CALLED: CCHECK; CHEHIC, SHIFTS;
SOURCE; UVFLDX; BHFLUX;
XYDIFF; ZZDIF
I/O PEBFORHED: None
This routine is the control program, which supervises the
proper execution of the steps for obtaining the numerical
solution of the concentration equation. It calls various
computational subroutines to compute the concentration field
contribution of each term of the equation. It is executed
each time step in the numerical integration of the equation.
The hourly average concentration field is also computed here.
The general flow diagram of this control program is shown in
Tig. II-3,
-------
PAGE 74
Operation Sequence Subroutine Called Variables Defined and/or Comments
No
Variables passed
from MAIN
1
Specification
Statements
1
CALL
1
CALL
1
CALL
1
CALL
1
CALL
1
CALL
1
CALL
1
CALL
1
\. ? ;S
j^es
CALL
1
CALL
» 1
*+
CALL
SHIFTN
SHIFTN
SOURCE
SHIFTN
UVFLUX
SHIFTN
(See Table 1 1-2 for argument listing)
'DIMENSION'; 'COMMON'
COLD
C
ZS,ZRR,XP,YP,ZP,CP1
(Add new sources)
C
CP1
(Computes U-flux)
C
UVFLUX
SHIFTN
WWFLUX
SHIFTN
XYDIFF
CP1
(Computes V-flux)
C
CP1
(Computes W-flux)
C
CP1
(Computes horizontal diffusion)
Figure 11-3. Detailed Flow Diagram of SUBROUTINE AACOMP
-------
P1GE 75
Operation Sequence Subroutine Called Variables Defined and/or Comments
CALL
SHIFTN
C
CALL
ZZDIFF
CP1
(Compute vertical diffusion)
CALL
SHIFTN
C
CP1
(Compute chemical decay)
CCHECK
ICHECK
(Check if steady state being reached?)
(Next time step will be remainder of hour)
Figure 11-3. (continued)
-------
PAGE 76
SUBROUTINE AKZCAL: CALLED FROM: HACK
ROUTINE CALLED: None
I/O PERFORMED: None
This subroutine computes the vertical eddy diffusivity for
the surface layer AKZ(I,J) and its vertical variation in a
function form AKF(K), This routine is called whenever new
meteorological data is input into the model. The routine can
be divided into three parts. The first part computes the
bulk Richardson number RIB, if upper air measurements are
available (only for reference purpose). The second part
computes eddy diffusivities for the surface layer from the
continuous atmospheric stability class S, and surface
routhness ZO(I,J), The equations used in this computation
are given in Section I-2-B. Finally, a coaputation is
performed of the vertical variation of eddy diffusivity
according to equation (10) given in Section I-2-B.
SDBROUTINE C&DJDS: CALLED FROM: MAIN
ROUTINE CALLED: None
I/O PERFORMED: None
-------
PAGE 77
This subroutine adjusts the concentration values of upper
levels (froa grid KHH to KH) after a change in the volume of
grid elements due to the variation of mixing height. This
routine is executed every one hour of real-tiae computation,
and when the ratio of old ti»e step and new tiae step grid
voluae BH is not equal to one. The coaputation method is
given in Section I-2-C, equations 23 and 24.
SUBROUTINE CCHECK: CALLED FHOH: AACOMP
ROUTINE CALLED: None
I/O PEBFOBMBD: None
Subroutine CCHECK checks whether the computed concentration
field in each tine step has reached a steady state solution
of the concentration equation. Sets ICHECK=0 if steady state
is reached, otherwise sets ICHECK=1.
ENTRY CDTOTP: CALLED FBOH: MAIN
ROUTINE CALLED: None
I/O PEBFOBHED: Yes
-------
PAGE 78
This is an entry in subroutine CQNSIN. It only executes once
in the beginning of the program. It is used to copy card
input data to I/O unit = IUNIT2 and to print out input data
in card image format.
SUBROUTINE CHEHIC: CALLED FROM: AACOMP
ROUTINE CALLED: None
I/O PERFORMED: None
This subroutine computes chemical decay of SO, It is
executed every time step, if ICHEH^0 froo NAHELIST input.
The equation of coaputation is given in Section I-2-B.
SDBRODTINE CONSIN: CALLED FROM: HAIN
ROUTINE CALLED: None
I/O PERFORMED: None
-------
PAGE 79
In this routine, values of model constants are specified. It
only executes once at the beginning of the program.
SOBBOOTINE DIHENS: CALLED FBOH: MAIN
ROUTINE CALLED: None
I/O PERFORMED: None
MULTIPLE EHTBY: DIHEN1; HHCALC
This subroutine is called once in the beginning of the
program. It initializes the necessary parameters for the
grid system of the model. It also initializes the mixing
height HH at height of planetary-boundary layer.
ESTRY DIHEN1: CALLED FBOH: BAIN
ROUTINE CALLED: None
I/O PEBFOBKED: None
This routine is executed every hour of real time when mixing
height is recomputed. It is an entry of subroutine DIHENS
-------
PAGE 80
and coaputes vertical grid syste® doe to changes in the
mixing height. The computation of the vertical grid systea
as function of mixing height is discussed in Section I-2-C.
SUBROUTINE DTTESI: CALLED FBOM: MAIN
ROUTINE CALLED: None
I/O PERFORMED: None
This subroutine computes and checks the tine step for the
numerical integration of the concentration equation. The
numerical stability criteria for horizontal advection,
horizontal diffusion and chemical reaction are all included,
See Section I-2-C for detailed stability criteria.
SOBBOailNE GEOIN: CALLED FROM: MAIN
BOOTINE CALLED: XYUTMS; XYUTH1
I/O PERFORMED: Yes
This subroutine reads in geographical and annual emission
-------
PAGE 81
data. It executes once at the beginning of the program. The
input variables are A) Origin of nuaerical grid in UTM
coordinate; B) UTM coordinates of BAMS stations
XSUTM(L),ISOTM(L); C) Surface roughness parameter ZO(I,J); D)
Effective emission height of area sources ZS(I,J); E) UTM
coordinate of point sources and their stack heights,
XPUTM(L), YPUTM(L) and ZP(L); F) Annual averaged area source
esission rates QB(I,J) and their sum QBTOT; G) Annual average
point source emission rates PQB(L) and their sum PQBTOT; and
H) Normalized plume rise on each point source ZH
-------
PAGE 82
programmed aethods are discussed ID Section I-2-B, This
routine is subject to modification or revision when store
information about mixing height or vertical stratification of
boundary layer atsosphere in the St. Louis area is
available.
SUBROUTINE OUTAPE: CALLED FBO«: HAIN
ROUTINE CALLED: None
I/O PEBFOBMED: Yes
This subroutine outputs the computed results on two I/O units
= KDNITP and KUNITC. The detail output specification is
given in Section II-3-D. This routine is executed in two
ways. During the initialization of the program (ITH^O), it
sets up I/O units KUNITP and KUSITC and checks whether there
are previously computed results stored in these I/O units.
If the answer is yes, it sets LHARH=1, and gets these I/O
units ready for additional writes. Otherwise, it sets
LHA8H=0 and returns to main prograa. In the subsequent
execution, it is only called every LPBINT interval of real
time to output the averaged concentration values into these
two I/O units according to the output specification.
-------
PAGE 83
SUBROUTINE POSIT?
This subroutine sets variable A(I,J,K) - 0, if it is a
negative value. Currently, it is not used in the model. It
is included in the program for the possibility of using other
finite difference methods which Bay result in negative values
of concentration due to nuaerical or truncation error*
SUBROUTINE PRINTS: CALLED FROM:
ROUTINE CALLED:
I/O PERFORMED:
MULTIPLE ENTRY:
MAIN
WHITEX, iRITEZ
Yes
PRINTA, PRINTS, PRINTC
This subroutine prints out geographical and annual emission
data according to output specification (Section II-3-C), It
also calls subroutine IRITES to initialize the dummy
arguments needed in entry IRITEX and iRITEZ. It only
executes one tiae in the program.
-------
PAGE 8U
ENTRY PHINTA: CALLED FBOM:. MAIN
ROUTINE CALLED: WRITEX
I/O PEBFOBMED: Yes
This subroutine is an entry in subroutine PRINTS. It prints
out time varying emission rates of area and point sources.
It is executed every LTSOUS interval. See Section II-3-C for
output specification.
ENTRY PHINTB: CALLED FBOM: MAIN
ROUTINE CALLED: HRITEX
I/O PERFORMED: Yes
This subroutine is an entry in subroutine PRINTS. It prints
out time varying meteorological data, RAMS stations data, and
model parameters. It is executed every LTWIND interval. See
Section II-3-C for output specification.
-------
PAGE 85
ENTRY PBINTC: CALLED FBOM: MAIN
BOUTINS CALL: WBITEX, WBITEZ, STNCON
I/O PEBFOFMED: Yes
This subroutine is also an entry in the subroutine PBINTS.
It prints out computed results according to output
specification (see Section II-3-C). This routine is called
every time step. However, it is executed in two ways. For
every time step, it prints out instantaneous concentration
values at BAMS stations and returns to the main prograa at
the beginning of the time step loop. For every LPSINT
interval and ICBON^O, it prints out: A) Instantaneous surface
concentration field, C1,(I,J), B) Two vertical cross sections
of instantaneous concentration for I=IMC and J=JMC; C)
Averaged surface concentration field CC{I,J) for LPBINT
interval, D) Averaged concentration values at BAMS stations,
and C) 24-hour average concentration values at BAMS
stations.
-------
PAGE 86
SOBHODTINE SHIFTN: CALLED FROM:
ROUTINE CALL:
I/O PERFORMED:
MAIN; AACOHP
None
None
This subroutine shifts array A(I,J,K) to B(I,JfK). Its
execution is determined by the numerical method employed in
the model. See Fig. II-2 and II-3 for current usage of this
routine.
SUBROUTINE SOURCE: CALLED FROM: AACOMP
ROUTINE CALL: ZZGRID, XYUTM1
I/O PERFORMED: None
This subroutine adds new source emission into the system. It
is executed every tiae step. The program can be divided into
two parts: add area sources and add point sources. The
logical path which determines the execution of each part is
specified by LSOUS{1) and LSOUS(2). The method by which the
new sources are added into the system is discussed in Section
I-2-F. The effective stack: heights of point sources are
coaputed in this routine for each LPRINT interval.
-------
PAGE 87
SOBBQOTIME SOOSIN: CALLED PHOH: MAIH
HOUTIBE CALL: Hone
I/O PEBFORHED: Yes
This subroutine reads in tiie-varying source emission data
QB(I,J) and PQB(L). It is executed every LTSQOS interval.
See Section II-3-B for input specification.
SOBBOOTINB STABIT: CALLED FBOH: WINDIN
BOUTIME CALL: None
I/O PERFOfiSED: Hone
This subroutine estimates a continuous atmospheric stability
class as function of average wind speed, solar insolation,
and sky condition (cloud cover). It is executed every LTHIND
interval or upon the input of new neteorological data. More
data and further study are necessary to generalize this
subroutine. A method of analysis applied in this routine is
discussed in Section I-2-B.
-------
PAGE 88
SOBfiOUTINE STNCON: CALLED FROM: PEINTC
ROUTINE CALL: None
I/O PERFORMED: None
This subroutine stores concentration values measured at RAMS
stations on array IOBS(L) and computes simulated
concentration values at stations on array ICAL(L) by bilinear
interpolation of CC{I,J). It also computes 24-hour-average
concentration at stations KOBS (L) and KCAL (L) for observed
and computed values, respectively. In addition, the average
values of each array CAL (L) , IOBS (L), KCAL(L) and KOBS(L) for
L=1, NS are calculated. These values represent the spatial
average concentration and they are stored on same array with
index of NSX. NS is the total aumber of BAMS station and
NSX=NS+1. This routine is executed every LPRINT interval.
SDBBODTINE TIHEX: CALLED FBQM: MAIN
BOUTINE CALL: None
I/O PERFOBMED: None
-------
PAGE 89
This subroutine fixes the tiae indices for both siaulation
tiae and real time. It is executed every tiae step.
SUBROUTINE UVZF: CALLED FHOH: «INDIN
HOUTIBE CALL: ilNDBB
I/O PERFORMED: Done
This subroutine computes vertical variation of horizontal
wind coaponents U and V. Three aethods are included in this
routine. The choice of aethod of coaputation is controlled
by LUIND, which is defined in NAHELIST input. These three
aethods are discussed in Section I-2-D. This routine is
executed upon the input of new aeteorological data.
SUBROUTINE UVFLUX: CALL HOOTIME: AACOMP
ROUTINE CALLED: None
I/O PERFOBHBD: None
-------
PAGE 90
This subroutine computes horizontal advection terms in the
concentration equation. It is executed every tiae step. The
current version of the program is the second-order central
finite difference scheme with flux correction. It has been
discussed in detail in Section I-2-C. The program works for
both X-direction and Y-direction advection Cy interchanging
the actual arguments in call statements. The actual
arguments listed in Table II-2 for this subroutine is to
compute x-direction advection. For computing y-direction
advection, the actual arguments in the CALL statement is
(CP1, C, U, UZF, DY, DX, BDY, 8DX, Z, EK, FK, JM, IM, KM,
IMJH, IJKM, LY, LX, IJO),
SUBROUTINE UVINTP: CALLED FROM: BINDIN, WIHDGB
ROUTINE CALL: None
I/O PERFORMED: Hone
This subroutine obtains the surface wind field at all
numerical grid points by bilinearly interpolating analyzed
wind fields over a wind analysis grid system, which has a
grid size of DELX. {It is assumed that DELX>DX, DY). This
routine is executed whenever there is an input of new wind
-------
PAGE 91
data.
SUBROUTINE SWZF: CALLED FROM: WINCIM
ROUTINE CALL: None
I/O PEBFOBHED: None
This subroutine computes the vertical conponent of wind
W(I,J,K). The control paraaeter LWi supervises the execution
of this subroutine. For LWW=0, it assumes W (I,J,K)=0.0. For
L»H=1, a continuity equation is used to coapute W(I,J,K),
For LW8=2, after obtaining W(I,J,K), a nine-point averaging
scheae is used to smooth the W field. This routine is
executed after each new input of Meteorological data.
SUBROUTINE WINDER:
CALLED FHOH: WINDER, WWZF
ROUTINE CALLED: None
I/O PERFORMED: None
This subroutine converts wind vector (given wind direction
-------
PAGE 92
and speed) to wind components in X and Y directions or vice
versa.
SUBROOTINE HINDER: CALLED FROM: WINDIN
ROUTINE CALL: WINDER, OVINTP
I/O PERFORMED: Yes
This subroutine generates surface wind field from observed
wind data at RAMS stations. Two methods are included in this
program. The first oethod assumes a spatially uniform wind
field which is assigned a value from a wind observation at a
station nearest to the center of the city. The second method
employs a weighted interpolation scheme to obtain a wind
field on the wind analysis grid, UU (1,3), VV (I,J) which has a
grid size of DELX. Then subroutine UVINTP is called to
obtain U(I,J), V(I,J) at all numerical grid points. The
interpolation scheme is discussed in Section I-2-D.
-------
SOBBOUTIHE HIHDIH: C1LLED FHOH:
ROUTINE CALL:
I/O PEBFOBHED:
PAGE 93
HAIH
MINDBB, ilHDEB, UVIHTP,
STABIT, UVZF, HWZF
Yes
This subroutine reads in Meteorological data and data
observed at BAHS stations. Thereafter, various subroutines
are called to generate a three-dimensional wind field at
numerical grid points. It also calls subroutine STABIT to
obtain a continuous atmospheric stability class. It is
executed every LTUIHD interval* The input specification is
discussed in Section II-3-B.
SUBBOUTIHE iBITES:
CALLED FBOH:
BOUTIHE CALL:
I/O PEBFOBMED:
HOLTIPLE EMTBY:
PBIMTS
None
Hone
HBITEX, tfBITEZ
This duiay subroutine is executed once in the beginning of
the program to initialize the dimension and data statements
needed in the entries HBITEX and HBITEZ.
-------
PAGE 94
ENTRY WHITEX: CALLED FROM:. PBINTS, PSINTA, PRINTC
ROUTINE CALL: None
I/O PERFOBMED: Yes
This subroutine is an entry in subroutine WRITES. It prints
out array Q(I,J) under the title 'TITLE'. A variable format
is used in the WRITE statements. Parameter IFQRM determines
the number of columns to be printed on each line and controls
a page layout.
ENTRY WRITEZ: CALLED FBOW: PBINTC
ROUTINE CALL: None
I/O PERF08MED: Yes
This subroutine also is an entry in subroutine WRITES. It
prints out two vertical cross sections of array A at J=JMC
and I=IMC.
-------
PAGE 95
SOBBOUTIME HHFLUX: CALLED FROfl: &ACOHF
ROUTINE CALL: None
I/O PERFOBHED: None
This subroutine computes the vertical advection term in the
concentration equation. The sane numerical scheme as used in
subroutine OVFLOX is applied in this routine. It is executed
every time step. See Section I-2-C for a detailed
discussion.
SOBBOQTINE XYDIFF: CALLED FBOH: AACOHP
BOUTINE CALL: None
I/O PERFORMED: None
This subroutine computes the horizontal diffusion terms in
the concentration equation. It is executed every time step.
The second-order central finite difference scheme is used.
i
See Section I-2-C for detailed discussion.
-------
PAGE 96
SUBROUTINE XYUTMS:
CALLED FBOM: GEOIN
ROUTINE CALL: None
I/O PEfiFOflFSED: None
MULTIPLE ENTRY: XYUTM1
This subroutine computes UTM coordinates for all numerical
grid points. It is executed once at the beginning of the
program.
ENTRY XIUTM1: CALLED FROM: GEOIN, SOURCE
ROUTINE CALL: None
I/O PEHFOBMED: None
This subroutine converts point (XUTM,YUTM) in UTM coordinates
to point (XD,YD) in coordinates based on numerical grid
system which has an origin of (0,0). It is an entry in
subroutine XYUTMS.
-------
PAGE 97
SUBBOUTINE ZZDIFF: CALLED FROM: AACOHP
ROUTINE CALLED: Hone
I/O PEBPOHMED: Hone
This subroutine computes the vertical diffusion term in the
concentration equation. It is executed every time step. An
implicit second-order central finite difference scheme with
variable grid sizes and variable eddy diffusivities is used.
See Section II-2-C for detailed discussion.
SUBROUTINE ZZGRID: CALLED FROM: SOURCE
ROUTINE CALLED: None
I/O PERFORMED: None
This subroutine converts height ZA in units of meters into
numerical grid unit ZB.
-------
PAGE 98
2-D.
VARIABLES OSED IN THE PROGRAM
The most important variables and constants used in the
program are listed in this section. They are arranged in
alphabetical order. This listing includes the array size,
at which part of the program that array been defined, the
unit used and definition of each variable.
NAME ABRAY HHERE UNIT DEFINITION
SIZE DEFINED
AK CONSIN
AKA(N) 24 INLIST
AKF(K) KM AKZCAL
AKH(K) KM INLIST
AKZ (I, J) IFUM AKZCAL
A1 (N)
CONSIN
C (I, J, K) IJKM AACOMP
CC (I, J) IMJM AACOMP
COLD {I, J) IMJM AACOMP
CP1 (I, J,K) IJKM AACOMP
Von Karman constant,
sec~l Chemical reaction rate constants
for hour N of a day,
Function determines the vertical
variation of AKZ.
ia2/sec Normalized horizontal eddy diffusi-
vity over layer K.
m2/sec Surface vertical eddy diffusivity
at grid point (I,J),
Four constants in in Businger's
foraula.
ug/ffi3 Last computed concentration values
at grid point (I,J,K).
ug/sa3 Average concentration at grid point
(I,J) over LPRINT time interval.
ug/is3 Surface concentration at grid point
(I,J) for previous time step.
ug/tn3 New computed concentration value at
grid point (I,J,K).
C1(I,J) IMJM AACOMP ug/m3 Surface concentration layer of CP1.
-------
PAGE 99
DCHIN
DT
DX(I)
DXS(I)
DY(J)
DYS{J)
DZ (K)
DZS(K)
EK,FK(N)
RFZ
HG
HH
HHAX
HMIN
HP
HS
ICAL(N)
IDAY
IDAYTP
IHR
IHRTP
IJKM
DTTEST
IM INLIST
IK DIMENS
JH IMLIST
Jfl DIMENS
KM DIMEN1
KH DIMEN1
LH
AKZCAL
INLIST
HHCALC
ISLIST
INLIST
INLIST
INLIST
NSX STNCON
TIMEX
INLIST
TIMEX
INLIST
CONSIN
sec
m
.«
m
m*
m
n*
m
m
m
m
m
m
ug/i
IHLIst ug/a3 Value of criterion for checking con-
vergence of concentration field.
Time step in numerical integration.
Grid interval size in X direction.
DX(I) squared.
Grid size in Y direction.
DY(J) squared.
Vertical grid interval size.
DZ(K) squared.
Temporary storage array.
Integral of non-dimensional wind
shear.
Thickness of planetary boundary
layer-
Height of nixing layer.
Maximum height of mixing layer
during a day.
Minimum height of mixing layer
during a day
Height of upper wind measurement.
Effective height of surface wind.
ug/m3 Computed average concentration value
at monitoring station N; N=NSX is
spatial average of ICAL{N), N=1,NS.
Day of the month.
Starting day for which the computed
and observed concentration value are
stored on I/O UNIT=KUNITP and KUNITC.
Time of the day.
Starting hour of IDAYTP.
IM*JH*KM.
-------
PAGE 100
IMJM
IHO
IBM
IN
IOBS (N)
IS{N)
ITM
ITOTHB
ITSEC
ITSTEP
ITOBS(N)
IDNIT
IUNIT1
IDNIT2
IUTH (I)
CONSIN
TIflEX
CONSIN
INLIST
NSX STNCON ug/m3
&S GEOIN
TIMEX
TIMEX Hour
TIMEX 9ec
TIMEX
NSX STNCON
MAIN
MAIN
MAIN
IM XYDTMS ffl
IYR
JM
JHC
TIMEX
INLIST
INLIST
Number of grid points in X direction.
Grid point in x direction where a
¥-Z cross section of concentration
field to be printed.
IM*JM.
Month of the year.
IM-1.
Number of wind grid points in X
direction for which the analyzed
wind field is stored,
Observed average concentration at
monitoring station N; N=NSC is spatial
average of IOBS (N) , N=1,NS.
Identification number of monitoring
station N.
Equals TM, simulated time in seconds
(integer) .
Total real time being simulated in a
run.
Simulation time starting from each one
hour interval (integer, ITSEC^TSEC).
Number of time steps being simulated.
Number of observed datum values at
monitoring station N during a 24-hour
period.
An I/O unit for card input (=IUNIT2).
I/O unit for card reader.
I/O unit for temporary storage unit,
UTM coordinate of numerical grid
point in X direction.
Year of the centry.
Number of grid points in Y direction,
Grid point in y direction where a X-Z
-------
PAGE 101
JH1
JM
JONIT
JXS(N)
JYS(H)
ON
KOBS (N)
KUHITC
KONITG
KDHI5TP
KUNITS
KONITW
KHIND
CONSIN
INLIST
INLIST
HS GEOlN
NS GEOIN
KCAL(N)
KB
KH1
KM
NSX STNCON
INLIST
CONSIN
INLIST
INLIST
NSX STNCON
INLIST
INLIST
INLIST
INLIST
INLIST
INLIST
cross section of concentration
field to be printed.
JM-1.
Sane as IN, except in Y direction.
I/O unit for line printer.
X coordinate of monitoring station
N in grid units.
Y coordinate of monitoring station
H in grid units.
2U-hour average of ICAL(N).
Number of grid points i& vertical
direction.
KM-1.
Number of levels at which 3-dimen-
sional wind field is computed.
Number of fixed grid intervals in
vertical direction.
ug/a3 24-hour average of IOBS (N) .
I/O unit for storing computed and
observed average concentration
(ICAL,IOBS) .
Input unit for geographical and
annual emission data.
I/O unit for storing computed surface
concentration field.
Input unit for the varying source
emission data.
Input unit for time-varying meteoro-
logical data (RAMS) and surface wind
field (UU,VV).
Control flag for the choice of
obtaining surface wind field.
=1, Input the preprocessed objective
analyzed wind field {01,?1) ;
-------
PAGE 102
LCHEM
INLIST
LCRON
INLIST
LHJUS
INLIST
LM
LPRINT
GEOIN
CONSIN sec
LPQ
INLIST
= 2, Input the subjective analyzed
wind field (UU,VV) ;
= 3, Input BAMS data, spatially
uniform wind field is assumed;
=4, Input HAMS data, variable wind
field is generated.
Control flag for computing chemical
reaction.
=0, Chemical decay of pollutant is
not computed;
-1, Compute chemical decay.
Control flag for a run.
-0, I/O test run;
=1, actual run.
Control flag for adjusting concen-
tration values due to the change in
the voluiae of grid cells,
=0, keep total mass constant;
= 1, keep mass constant, if RH <1,
otherwise, mass is not constant;
=2, keep mass constant, if HH >1,
otherwise mass is not constant,
= 3, mass is not constant,
Total number of point sources.
Minimum value of LTSoUs and LTHIND,
Time average concentration field CC
is computed over this time interval.
This parameter is also used for
controlling the logic path of execut-
ing the main progran and certain I/O
operations.
Control flag for modeling point
sources.
=1, source is added at one grid
point only with downwind
distance of U*min (T1,T2);
-------
PAGE 103
LSODS(M) 2
INLIST
LTOP
INLIST
LTSTOP
LTSOUS
LTWIND
LWARM
INLIST Hour
INLIST Sec
INLIST Sec
OUTAPE
LVIND
INLIST
=2, source is added at four neigh-
boring grid points;
= 3, same as =1 but »ita downwind
distance of V/DT;
=.4, same as =2 but with downwind
distance of V/DT.
Control flags for adding area (M=1)
and point (M=2) source into the system.
=0, *M' source is not considered;
=1, 'H1 source is injected.
Control flag for choice of upper
boundary condition for vertical
diffusion computations,
=0, sets C=0 at the boundary.
=1, reflecting boundary.
Maximum hours of real time to be
simulated in a run.
Time interval of the time-varying
source emission data that is to be
input.
Time interval of the time varying
meteorological data is to be input.
A parameter indicating the status of
I/O UNIT = KUNITC and KUNITP.
=0 I/O units are initiated in this
run;
=1 I/O units contain results from
previous runs.
Control flag for choice of computing
vertical variation of horizontal wind
components
-------
PAGL 104
LWBITE(N) 10 INLIST
LWTOP
INLIST
LWW
INLIST
NEAR (Ir J) IN, JN WINDES
NMONDY(N) 12 MAIN
NP(L) LM GEOIN
NS
NSX
OLMIN
P ABM(N)
INLIST
MAIN
INLIST
10
N=1
HINDIN
WINDIN
=3, vertical change in wind direction
is assumed linearly proportional
to the height.
Control flag for output options (see
Section II-B-C).
Control flag for choice of upper
boundary condition for vertical
advection computations.
=0, no vertical flux crosses the
boundary;
=1, flux is allowed to cross the
boundary,
control flag for choice of computing
vertical wind component W(I,J,K).
= 0, set H (I,J,K) ^Q;
= 1, W(I,J,K) is computed by conti-
nuity equation;
= 2, after computing W(I,J,K) by con-
tinuity equation, it is smoothed
again by 9-p»int average.
The nearest RAMS station to wind grid
point (I,J) has observed wind data.
Number of days in Nth month of a year.
Point source identification number
(refer to the NEDS data set).
Total number of monitoring stations.
NS+1.
Minimum value of Monin-Obukhov length..
Array to store various meteorological
and source parameters.
m/sec Mean wind speed.
deg. Wind direction measured at center
of the city.
N=3
WINDIN °C
1st-level temperature measured at
center of the city,
-------
PAGE 105
»N=5«
H=7«
• H= 10'
PHIFHZ
PHAX
PHIN
PQA(L)
PQB(L)
PQBSOH
PQBTOT
QA(I,J)
QB(I,J)
QBSUH
LH
LH
HINDIS
flAIN
SOOSIN
SOUSIN
ilNDIN
WINDIN
AKZCAL
AKZCAL
INLIST
INLIST
SOUSIN
SOUSIN
Continuous atmospheric stability
class (-S) .
a =HH, ailing height.
g/sec =• QBTOT, total area source eaission
rate.
g/sec =PQBTOT, total point source emission
rate.
°C
2nd-level teiperature leasurei at
center of the city.
Radiation or sky condition.
=OL, Hanin-ObukhoT length.
Non-di«ensional temperature gradient.
NaxiiUB ralue of exponent in wind
profile pover law.
SOUSIN
SOQSIN
IBJH SOOSIN
IHJM SOUSIN
GEOIN
value of exponent in wind
profile power law.
g/sec Emission rate of point source L.
g/sec Eiission rate of point source L at
current time plus LTSQOS tiie
interval (Not used in the current
program).
g/sec Total point source emission rate
used in the model.
g/sec Total point source emission rate.
g/sec/ Area source strength at grid point
km* (I,J) at current time plus LTSOUS
time interval (not used in the
current program) .
g/sec/ Area source strength at grid point
km*
-------
PAGE 106
RAMS (N)
• N=1 '
•N=2'
•N-J'
1 N= 5 '
•N=6'
PDX (I)
RDY (J)
RDZ {K)
RH
RIB
TM
6 HINDIN
-
-
-
-
-
IM DIMENS
JM DIMENS
KM DIMEN1
DIMEN1
AKZCAL
TlMEX
TSEC
TIMEX
m/sec
deg.
°C
°C
ug/m3
U(I,J,K) IJKN HINDIN
USTN (N) NS HINDER
UU (I, J) IN, JN WINDIN
UO AKZCAL
01 (I,J) IMJM HINDIN
UZF(K) TM UVZF
V(I,J,K) IJKM WINDIN m/sec
sec
sec
m/sec
tn/sec
m/sec
lli/SCC
m/sec
RAMS station data.
Surface wind speed.
Surface wind direction.
1st-level temperature,
2nd-level temperature.
S02 concentration.
Radiation or sky condition.
DX (1) /DX (I- 1) , 1=2, IM,
DY (J) /DK (J-1) , J = 2, JM,
DZ {K)/DZ (K-1) , K=2, KM-
Ratio of old and new vertical grid
size at upper level.
Bulk Richardson number.
Simulated time in seconds {floating
point number).
Simulated time starting from each
one-hour interval (floating point
n umber).
wind component in X direction at
grid point (I,J,K).
Surface wind component in X direction
at monitoring station.
Surface wind component in X direction
at wind grid point (I,J).
Surface friction velocity.
Surface wind component in X direction
at grid point (I, J) .
Function determining the vertical
variation of wind component in X
direction.
tfind component in Y direction at
grid point (I,J,K).
-------
PAGE 107
TSTH(M) NS WINDGR a/sec
?7(I,J) IN,JN HINDIN m/sec
VI (I,J) IHJH HINDIN B/sec
TZP(K) KB OVZF
I(I,J,K) IJO HWZF a/sec
flZF(K) KB HWZF
(D
LH GEOIN
XPUTM(L) LH GEOIN KB
85 GEOIN
XSUTB(N) US GEOIN fca
LH GEOIN
YPUTB(L) Lfl GEOIN
YS{»)
NS GEOIN
YSOTB(H) NS GEOIN lea
Z(K) KB DIBEN1 to
ZM(K) KM DIMEN1 •
ZBAX INLIST n
ZP(N) NS GEOIN B
Surface wind coaponent in T direction
at Bonitoring station N.
Surface wind coaponent in Y direction
at wind grid point (I,J).
Surface wind coaponent in Y direction
at grid point (I,J).
Function deteraining the vertical
variation of wind component in Y
direction.
Wind coaponent in Z direction at
grid point (I,J,K).
Function determining the vertical
variation of wind component in Z
direction.
X coordinate of point source L in
nuaerical grid units.
UTB coordinate of point source L in
X direction.
X coordinate of HiflS station N in
nuaerical grid unit.
UTB coordinate of HAMS station H in
X direction.
Y coordinate of point source L in
nuaerical grid units.
UTM coordinate of point source L in
Y direction.
Y coordinate of BABS station N in
numerical grid units.
UTB coordinate of BABS station N in
Y direction.
Height of vertical grid point K.
Height of points in the aiddle of
vertical grid eleaent.
Baxiaua aixing height used in the
aodel.
Physical stack height of point
source N.
-------
PAGE 108
ZPR (N) NS SOU.RCE
ZR
ZRPQ
ZOMEAN
NS SOUSIN
INLIST
ZRTSE INLIST
ZS (I,J) IMJM GEOIN m
ZO(I,J) IMJM GEOIN m
INLIST m
Effective stack, height of point
source N in numerical grid uuits,
Normalize plume rise of point source
N.
Constant used to get plume rise from
its emission rate (if plume rise is
not supplied).
Parameter for adjusting plume rise.
Effective emission height of area
source at qrid point (I,J).
Surface roughness length at grid
point (I,J).
The uniform value of surface rough-
ness for entire region, If ZOMEAN^O. ,
then original ZO(I,J) is used.
-------
PAGE 109
3. Computational Procedures
The model requires three input data files and two output
data files. The current program provides various options in
I/O operations. The users can exercise these options
through the paraneters in included in the NAMELIST input and
jJob Control Language (JCL) specifications. In this section,
we shall describe the specifications and requirements of the
I/O devices for the model and the JCL set-up to run the
program. Table II-4 lists the I/O units used in this
program. The third column of the table indicates that I/O
unit is for input (I) or for output (0), The fourth column
indicates in which subroutine that I/O operations are
executed.
-------
PAGE 110
OMIT
DS8AHE
IZP.
IUNIT
IUNIT1
IUNIT2
JUNIT
JUNIT1
KUNITG EPAGE02
KUNITS EPASOUS
KUNITH 1INDDATA
KUNITH HINDDATB
KUNITH EPARAHS
KUNITC EPASTN01
KUNITP EPACOSC1
I
I
I/O
0
0
I
I
I
I
I
o
0
SAIN
CDTOTP
CDTOTP
(ALL)
CDTOTP
GEOIN
SOUSIN
HINDIN
¥INDIN
HINDIN
OUTAPE
OUTAPE
'NAHELIST
(UNIT
(SCRA:
{UNIT
(UNIT
XRAMS,ZS(
KHR,KMO,[
KYR,KMO, I
KYR,KHO,!
KYR,KMO,i
IYR,IMO, J
IYR,IHO,;
(UNIT FOB CARD READER)
(SCRATCH STORAGE UNIT)
{UNIT FOR LINE PRINTER)
(OMIT FOB LINE PRINTER)
S,ZS,ZO,QE,PQB,
,KMO,KDAY,KYB,QB,PQB
,KMO,KDAY,KHR,U1,V1,RAH
,KHO,KDAY,KHR,UU,VV,RAM
,KMO,KDAY,KHH/RAHS
,IMO,IDAY,IHR,PARM,ICAL
,IHO,IDAY,IHH,PAHM,CC,
ICAL,IOBS
Table II-U I/O Units Used in the Prograa
-------
PAGE 111
3-A, Parameter specification
All the necessary paraaeters for the model are initialized
bv means of input NABELIST/IHLIST/- Some of the secondary
parameters and values of constants are specified in
subroutine COUSIN. Initialization of variable arrays and
certain constants are specified in DATA statements in
various subroutines vhen it is necessary. The numerical
grid system is initialized vhen subroutine DIHENS is called
and executed.
Except for NAHELISI input, all other initializations of the
model parameters are default processes and shall be
consistent with the aodel formulation and programming logic.
Therefore, in the following discussion we shall only deal
with NAfiELIST input.
The parameters included in the BAHELIST/INLIST/ and their
sample values are given in Table II-5. The definition of
each variable is given in Section II-2-D.
For better perspective, let us discuss these parameters in
various groups according to their functions in the program.
-------
PAGE 112
SINLIST
T M
«~* A
IK-3Q, JH=40,
LM=150,
i, n - i 3 u ,
IMC=15, JMC=19,
DX=5*2000.,20*1000.,5*2000.,DY=1
DZ=5*20.,9*25.,TM=0.0, DT=120,,A
HS=10., HP=140., HG-1000., ZMAX=
ZRPQ=0.6.ZBISE=1.0,Pa&X=0.5,PMIN
.,DY=10*2000.,20*1
120,,AKA=24*1.0E-4
Table II-5. Example of NAMELIST/IHLIST/to INPDT DATA
-------
PAGE 113
3-A-a) Paraaeters specifying grid system;
IH,JH,KH,KHN,IN,JB,KN,LM,NS,DX,DY,DZ.
DX, DY and DZ specify non-uniform grid size for a numerical
method in the x,y and z directions respectively. The choice
of DX and DY shall correspond to the area source inventory
grid. The dimensions of DX, DY, DZ are Ifl,JH,KM,
respectively. Therefore, they shall have IB values for DX
and so forth. IM,JM, KH,IN,JH,KN,LH,SS are also ased for
dynamic allocation of storage when various subroutines are
called. Hence, the values specified for these parameters
shall not exceed the number used in the DIMENSION statement
in the HAIH program. As mentioned in section II-2-A, this
program can run on three operational configurations, namely,
debug, test, and production run. To choose one particular
operational mode, the following rules shall be followed.
For the debug run, set KM=5, and KN=1; for the test run, set
KH*1»» and, 0=1; and for the production run, set KM=14 and
KM«7, In addition, the DIMENSION specification in SAIN
program for three dimensional arrays, namely, CP1, C, 0, V,
1, shall be changed accordingly. (The program provides three
cards for this purpose. The user only needs to select one
card and keep the other two as comment cards.)
3-A-b) £a£4JBeter_s fo£ tine indices;
TB,DT,IHB,IDAY,IMO,IYH,IDAYTP,IHBTP,LTSTOP.
TH shall b« set equal to zero. DT is the initial guess of
-------
PAGE 114
the time step for numerical integration. It will be
recomputed in the execution of the program,
IHB,IDAY,IHO,IYB specify the starting real time for which
the model computation is to be performed. LTSTOP determines
the number of hours to be simulated in a run. It is assumed
that the necessary input data for time-varying source
emission and meteorological data are available through the
tiae span of simulation. (See section II-3-B). IDAYTP and
IHBTP specify the starting day and the hour for which the
computed concentrations have been stored on I/O units KONITC
and/or KUNITP- (see section II-3-D).
3-A-c) Physical parameters:
AKH,AKA,HS,HG,ZMAX,flMIN,HMAXrZHPQ,ZBISB,PMAX,PHIH,
DCMIN, OLMIN, ZOMEAN
This group of the parameters deal with physical assumptions
currently employed in the program. The specification of
these values is confined to the context of the current model
formulation, and within this formulation they have
reasonable physical meaning. ie suggest that unless one is
completely familiar with the model and its programming
logic, caution should be observed in changing these values
from those given in Table II-4.
-------
PAGE 115
3-A-d) Parameters ££r I/O operation;
LTSOUS,LTiIMD,IHC,JHC,lIiHITE.
LWBITE control the output options, and will be described in
section II-3-C and section II-3-D. IHC,JMC shall have
values less than IH,JH respectively. These two numbers
decide where vertical cross-sections of concentration field
are to be printed out. LTSOUS and LTSIHD specify at what
time interval the time-varying source emission and
meteorological data were prepared. The units for these two
parameters are in seconds. However, they shall be chosen in
the multiple of 3600 seconds with minimum values of 3600
seconds™
3-A-e) Number for I/O unit ;
JUSIT,JUNITC,KUNI?G,KUNITP,KOHITS, KOHITM.
These numbers specify various I/O units to be used in the
program. The values given to these parameters shall be
consistent with JCL specifications. More detailed
description is given in section II-3-E.
3-A-f) Parameters for selection of jsodelj.iiHj oj>tions:
LCROH, LHJUS, LCHEM, LHW, LTOP, LWTOP, LSOOS, LPQ, KWIND,
LWIND.
The values given to these parameters shall be in accordance
with the users need and shall be confined to the context of
the logic used in the program. Each value given to a
-------
PAGE 116
parameter has a special leaning in the computational method.
Hence, they Bust be defined carefully. The variables
described in section II-2-D offer a quick reference.
3-B. Input specifications and requirements
In addition to the NAHELIS? input described in the previous
section, the model requires three sets of input data. They
are: a) geographical and annual emission data, b) time-
varying source emission data (if available), and c) time-
varying meteorological data.
3-B-a) Geographical and annual emission data.
These data are input through subroutine GEOIN. The I/O unit
is KUNITG, They are prepared by auxiliary programs (see
section II-4). Since this data set is prepared in
card-input format, the following discussion is based on an
80-column card format. The listing of this input data set
is given in Appendix 3.
This data set shall consist of eight (8) data groups. They
should be arranged in the following sequence.
-------
PAGE 117
(i) Region specification
IXBEG, IYBEG,IBEG, JBEG: Format = (415).
These four numbers appear in one card. IXBEG and IYBBG
specify the UTM coordinates of the southwestern corner of
the area source inventory map. IBEG and JBEG are the x and
y distances (in kilometers) from point (IXBEG, IIBEG).
Point (IXBEG*IBEG, IYBEG+JBEG) is the origin of the
numerical grid systea.
(ii) RAMS st at ion locations
HS Format = (15)
(FHTDS(H) ,H=1,10) Format = (10A4)
(IS(L), XSUTa(L), ISUTH (L) , L=1,NS) Format = (FHTDS)
Note that variable format is used to read in data of BAMS
station identification and locations. FWTDS specify the
format statement by which the following data are to be read
in. The same method is applied to input the remainder of
the data in KDNITG.
(iii) Surface roughness parameters
-------
PAGE 118
(FMTDS(N) ,N=1,10) Foraat = (10AU)
( (ZO (I,J) ,1=1,111) , J=1 ,JM) Format = (FMTDS)
(iv) Effective emission heights of error
(FMTDS (N) ,N=1, 10) Format = {10A4)
{ (ZS(I,J) ,1=1, IM) , J=1,JH) Format = (FMTDS)
(v) Pfiilii source locations and
LMAX Format = (15)
(FMTDS (N) ,N=1, 10) Format = (10A4)
(NP(L) ,XPUTM(L) ,YPUTM{1) ,L=1,LMAX) Forfflat= (FMTDS)
(vi) ASJiM^jL average emission rates of area source
(FMTDS (N) ,N=1,10) Format = (10AU)
((QB(I,J) ,1=1,IM) , J=1,JH) Format = (FMTDS)
QBTOT Format = {F10. 1)
(vii) Annual average emissioii rajte of £oin^ sources
(FMTDS(N),N=1,10) Fornat= (10AU)
-------
PAGE 119
(PQB(tf),L=1,LHAX) Format = (FMTDS)
PQBTOT Format = (F10. 1)
(viii) Normalized glume r^ge for point sources
FMTDS(N) ,N=1,10 Format = (10A4)
-------
PAGE 120
L=1,LR), (ZB (0),L=1,LB), QBTOT,PQBTOT.
Unformatted read statement is used to input this data set.
Therefore, it is necessary to prepare these data in KUHITS
based on unformatted write statement. KYR,KMO,KDAT,KHB
define the ending time of each averaged emission rate data.
3-B-c) Meteorological
Meteorological data are input through subroutine WINDIN
every LTilHD second of real time. Siailarily to time
varying source emission data, the unformatted write and read
is used in the I/O processes. Data shall also be prepared
in advance of the model computation and stored on I/O unit =
KOSITU. It shall at a minimum cover the entire time period
for which model computation is to be performed, i.e.
(ITOTHH*3600/LT«IHD) sets of these data.
There are four options in the model for constructing surface
wind fields. These options are controlled by parameter
K¥IMD. Thus, the meteorological data should be prepared
consistent with the choice of options.
K¥IND=1: KYB,KBO, KDAY,KHB, ((^1 (I,J) ,1= 1 ,IK) , J= 1, JM)
-------
PAGE 121
({¥1 (I,J) ,1=1, IB) ,J=1,JH) , ((BAMS(H,N) ,H=1,6) ,N=1,BS)
KiIHD=2: KYB,KHO,KDAY,KHB, < (OU{I,J) ,1=1, IN) ,J=1,JN),
((VV (I,J) ,1=1, IN) ,J=1,JM) , <{RAfiS(H,H) ,M= 1 ,6) , N= 1 , NS)
: KYfi,KHO,KDAY,KHB, { (BAMS(H,N) ,H=1,6) , N=1,HS)
3-C. Output specifications
For each computation run, the NAMELIST input is
automatically printed out in its original image. This is
executed in subroutine CDfOTP, For e?ery pass through the
•ain conputational loop in MAIN program, the real time of
the simulated period is printed out. It serves as a title
heading for other output.
All other print outputs of the aodel computations are
executed in subroutine PRINTS (including entry point
PBIHTA,PBIHTB and PBIHIC). The storage of computational
results on tape or disk is performed in subroutine QUTAPE.
ie shall discuss the storage of computed results in the
next section. In the following discussion, we shall only
deal with print output, its I/O unit is JONIT.
-------
PAGE 122
Depending on the aode of each computational run, this
program provides to the users the selection of which
variables are to be printed out. When the surface field and
vertical cross-section are needed, subroutine WRITES and
WRITEZ are automatically called to perform such duties. The
choice of print output is controlled by parameter LWRITE (N)
where the naximui value of H is ten. These ten integer
values of LWRITE shall be specified by the user in the
NAMELIST input. In table II-6 the options of printing each
variable is listed. The subroutine names listed on the
first column of table II-6 indicate that where the output
operations are executed. The 3rd coluan of the table
specifies the condition that the value of LWHITE is needed
in order to output or print out variables listed in the Qth
column. In general, the highlights of Table II-6 are as
follows:
a) If LH8ITE(H)=0, for all N, there will be no output
performed.
b) Each element of LHRITE control a group of
variables to be printed out. The function of each
-------
PAGE 123
LBRITE(B)
ROOTI8B
PBIHTS
COHTBQL FLAG FOB OOTPOT
H
B*(to
1
COHD.
.GE. 1
fARIABLES TO BE OOTPOT
IS,XSaTH,ISOTa,XS,7S,JXS,JIS(L) rL=1,
(ZO (I,J) ,ZS(I,J) ,QB(I,J) rQBTOT,QBSD«
PBIHTA
PRINTS
PBIHTC
OOTAPE
2
3
4
5
6
7
8
9
10
.GE.
= 2,
= 3 ,
.GE.
s
.GE.
35
.HE.
.GE.
.GE.
~
s
.GE.
s
* 2,
.GE.
.GE.
. GE,
=
.GE.
.GE.
1
4
4
1
2
2
3
0
1
1
2
3
1
2
4
1
3
1
2
1
1
PQBTOT,PQBSDM
QBTOT,PQBTOT,QBSOH,PQBSOH
QB(I,J)
PQB(L)
HAHS(H,L)
UU(I,J) ,fV(I,J)
01 (I,J) ,¥1 (I,J)
*8(I,J,1), (IF L8»=1)
U(I,J,K)rY(I,J,K)
4-8(1, J,K) , (IF LBW=1)
Z(K),BH
OZF(K) ,TZF(R)
«8ZF(K), (IF LMW=1)
AKF(K)
AKZ(I,J), FOB J=JH/2
AKZ(I,J)
PARS (L) ,OT
aO,PRIFHZ,BFZ,BIB
CP1 (JIS,JIS)
CP1 (IHC,J,K) ,CP1 (I, JMC,K)
CP1 (I, J, 1)
CP1 (I,J,K) , (IF IHR=0)
ICAL(L) ,IOBS(L)
* KCAL(L) ,KOBS(L) , (IF IHB=0)
CC(I,J)
ITBrIHO,IDAY,IHB,PABH(fl) ,ICAL(L) ,IOBS(L)
— OM I/O UHIT=KOHITC
IYB,IHO,IDAY,IBRfPABH{H) ,CC(I,J) ,ICAL(L)
S IOBS(L) — OM I/O OUT = KOHITP
Table II-6 Output Operations Controlled by LiHITE(H)
-------
PAGE 124
element of LHHITE can be summarized as:
LHEITE{1) : Geographical and annual average emission data;
LWBITE(2) : Time- varying source emission data
LHHITE(3) : RAMS station data and subjective analyzed
wind field (if KHIND=2)
L»BITE(U) : 3-D wind field at grid points and vertical
grid system (vertical coaponent only print
when LSB^O)
LWBITE{5)
LIHITE{6)
LWHITE{7)
LMHITE(8)
LHRITE(9)
LHBITE<10)
Eddy diffusivities
Tiae step and meteorological paraaeters
Instantaneous concentration fields
Average concentration fields
Storing computed results on I/O unit=KUNITC
Storing coaiputed results on I/O unit=KONITP
c) The largest value which can be assigned to each
element of LHKITE varies from one to four. When these
largest values are used, it will result in a large
volume of print output.
d) The optimum print output for a production run is:
L»RITE=1,1,1,1,2,2,2,2.1,1.
This will result in three pages of output per each
LPBINT interval of simulation.
A sample of print outputs is given in appendix 4.
3-D. Storage of gesglts
Frequently the user will wish to store the computational
-------
PAGE 125
results on a disk or a tape. Such storage will provide a
data base for farther analysis of the model performance or
for other usages, e.g., graphic display and cliaatological
study.
The storage of computed results is executed by subroutine
OUTAPE. The control parameters are LWaiTE(9) and
LWRITE(IO). They Bust be set greater than or equal to one,
if the user wishes to exercise this option. The I/O units
are KOHITC and KUNITP.
The variables to be stored on KUHITC are:
ITB,IMO,IDAI,IHB,(FARM(N),H=1,10),
(ICAL(H) ,B=1,NSX) , (IOBS (N) ,N~1,HSX) .
On KUHITP, they are:
ITB,IHO,IDir,IHB,
-------
PAGE 126
index shall be in chronological order. In addition, the
starting time (day and hour) that data were stored on these
I/O unit should be specified in the NiHELIST input. They
are IDAITP and IHfiTP. The following example illustrates how
it works:
Example 1i
LSRITE(9)=1, LHBITE{10}=0, IHB=0, IDiY=3, IHBTP=1,
IDAITP=2.
The program will write its results on I/O unit = KUNITC
following data for the 24th hour of the second day. If such
data cannot be found, then the run will be terminated due to
EOF or Error on I/O unit =KONITC.
Example 2:
LWHITE(9) = 1,1WRITE (10)=0,IHR=0,IDAY=1,IHHTP=1,IDAYTP=1
The program will assume that the I/O unit=KUNITC does not
contain any previously computed results. It will start its
output execution froa the beginning of the allocated space
on KONITC.
-------
PAGE 127
3-B. JCL Specification
To run the program, the user is required to prepare the
input data according to the above described input
specifications. He is also required to specify the JCL
cards. One of the aost important considerations in
preparing the job deck for running the program is that the
parameters specified in the NAMELIST input must be
consistent with JCL specification, with the operational aode
of the run, and with the intended modeling options. In this
section, »e shall briefly illustrate these factors. He
shall assume that all the necessary input data described in
section II-3-B are created, catalogued and stored on the
peripheral storage device. The DSHAKE for each data set is:
a) Geographical and Annual emission data
DSN=JCSU165.EPAGE02.DATA
b) Time-varying emission data
DSN=JCS4165.EPASOUS.DATA
c) Tine-varying meteorological data
DSM=JCS4165.EPAHAMS.DATA
There are various ways of running the program. The simplest
way is to use all the source decks and go through the
compiling, linkediting and executing steps. An example of
th« job set-up is given in Table II-7. In this example, it
is assumed that the source program of IBHAQ-2 has been
previously copied to a data set named JCS4165.IBaAQ2.POBT,
stored on a direct access device and catalogued. If the
source program is to be submitted on the card reader, the
-------
PAGE 128
//EPA35A JOB
-------
PAGE 129
//EPA35A JOB
-------
P1GE 130
JCL set-up is given in Table II-8- It is noted that in both
exaaples, JCL specification of the I/O unit for geographical
and annual emission data is 13; time-varying source emission
data is 14, time-varying meteorological data is 15, and
storage of coaputed results is specified by 16 and 17. All
these I/O units in JCL specifications are consistent with
NAMELIST input parameters of KUNiTG=13, KONITS=14,
KUNITI=15, KOHITC=16 and KONITP=17 respectively. I/O unit
=12 is for scratch storage.
Auxiliary Prograa—Input Data Preparation
In this section we shall briefly describe the auxiliary
program which processes and analyzes geographical and annual
average eaission data for the aodel input. The input
specifications were described in section II-3-B, above.
For this study, we obtained the following data sets:
1) Annual average area source emission data of S0_ for the
St. Louis area: This data set consists of 20<49 punched
cards. Each card represents an area source. The area
source grid is non-uniforo and consists of square
elements ranging from one to ten kiloaeters in size. It
-------
PAGE 131
covers an area of 200x140 kiloaeter squares. This data
set was provided by EPA.
2) NEDS point source data for the St. Louis area: This set
consists of about 4700 cards. There are 175 plants and
factories included in this source inventory. Many of
then have multiple stacks vithin a plant. The structure
and format of this data set is described in "Guide for
compiling a Comprehensive Emission Inventory"- published
by EPA (Heference 9).
3) UTH conversion table and program: This data set and
program is provided by EPA. It is to be used for
converting a location in longitude and latitude units to
a UTS coordinate system. It also can be used for
converting one zone of OTM as an extension of the
adjacent zone.
4} RAHS station locations: Twenty-five RAMS station
locations in longitude and latitude units are provided by
EPA in printed form.
5) Topographical and building information: These data were
obtained by us as a few pages of handwritten notes.
(Reference 10), It consists of information related to
average buildiag height, height of ground surface above
the bottom of Hiss. Biverr and the estimated surface
-------
PAGE 132
roughness. A nested grid system is used for these data.
The origin of tbe grid is placed at the intersection of
highways 67 and 40. These data is obtained from Dr. Fred
H. ?ukovich of Research Triangle Institute, Research
Triangle Park, H.C.
6) St. Louis Maps: They are published by U.S. Geological
Survey and by various oil companies.
It is obvious that the above-described data are not
coapatible with model input specifications. To convert
these data sets into the form required for model usage is
guite a time-consuming task. Our experience has indicated
that the systematic analysis of input data is as important
as the model computation. Therefore, the analysis
procedures described in this section may be useful in the
further model development.
The first objective of processing these data is to select a
unified coordinate systea, the computational region of the
model, and its grid systea. Although, we try to be as
objective as possible, some of the work still involves
extensive subjectivity. In this respect, the computer
programs can only systematically screen the data, tabulate
it and display it in map forms, and the investigator must
-------
PAGE 133
make the decisions.
The selection of a coordinate system is straightforward.
All the data should be based on the UTM coordinate system.
Consequentially, the numerical grid systea in the model
computation should also be related to this coordinate
system. However, the selection of the computational region
of the model and its detailed grid system requires careful
consideration. The accuracy of model computations and the
requirements of computer capability for performing these
computations are closely related to the characteristics of
the modeled region and its computational mesh. The factors
that need to be considered are as follows:
1) The computational region shall be sufficiently large so
that we can minimize the numerical error in the central
region (where we are most concerned) due to the specified
inflow boundary condition.
2) The computational region shall include the major
area-source and point-source emission in the area. This
will make the computed concentration field consistent
with the total pollutant mass emitted.
3) The computational region shall encompass all the RAMS
stations (if possible). This will permit utilization of
-------
PAGE 134
the data obtained from the BiWS stations.
4) The numerical grid resolution shall be compatible with
the area source emission inventory grid. It is desirable
that numerical confutations shall be based on the finest
mesh size available in the area source emission
inventory. Thus, the accuracy of numerical computations
will be compatible with the given input data.
5) In addition to the above-mentioned criteria, we should
keep in aind the computer capability for performing the
computation involving a large region to be covered, and
utilizing a fine grid resolution. The CPU core storage
and CPO time required to run a job may outweight other
factors. Our most important concern in this work was
that the model be operable on EPA*s computer facilities
in a reasonable amount of CPO time.
It is obvious that the selection of the computational region
of the model and its grid system shall be performed
concurrently with other data analysis. Based on the
criteria described above, we arrived at the following design
features:
1) The computational region employed is 40x60 km centered
-------
PAGE 135
at RAMS station 101 (center of St. Louis city). This
total area is about the sane as in our earlier model and
includes all the major area sources from St. Louis, E.
St. Louis, Alton city and Granite City and the major
point sources located at Alton and Merames. A total of
21 RAMS stations are vithin this region.
2) The smallest grid element of the area source inventory
is 1x1 km . It is clear that the model computational
grid must have the same resolution.
3) With 14 levels in the vertical grid adopted , the 30x40
element horizontal grid is the maximum resolution
acceptable by most computer facilities. {Note that
increasing the number of grid cells also means increasing
the required CPU time for computation) .
4) A non-uniform numerical grid is adopted in order to
cover the 40x60 km area with a 30x40 element grid
system. He call this non-uniform system a "stretched"
grid. The center of the computational region is a 20x20
grid of one kilometer squares. This is compatible with
the finest area-source resolution. In the outer area,
the grid elements are stretched, and become 2x1, 1x2,
and 2x2 km in size.
In addition to the above analysis, the St. Louis map was
-------
PAGE 136
digitized. This is for graphic display purposes. Now, let
us briefly discuss our analysis procedure and the function
of each auxiliary program.
4-A. Area Source data
1) Objective
o To make the best use of NEDS data for the model and
obtain the information for determining the computa-
tional region and the numerical grid system.
2) Pertinent factors:
o Geographical distribution of area sources and their
emission strengths.
o Grid resolution used in the area source inventory.
o Area source distribution as related to location of
BAMS stations.
o Computational requirements of the model.
3) Analysis procedures:
o Screen the original data set for error and
consistency check.
o Convert the original data to a uniform grid
-------
PAGE 137
of 1x1 ka size and compute total emission rate
from all sources.
o Print data in aap fora for checking their
distribution.
o Choose 40x60 grid elements and check the sub-total
emission rates. (This is iterated several tines).
o Convert 40x60 grid system to 30x40 stretched
system which is coapatible to the computational
(numerical) grid.
o Convert eaission rate to density units.
o Estimate approximate effective emission height
or each area source from current data and from
1964-1965 data.
o Check emission data with 1964-1965 monthly
emission data.
4-B. Point source data
1) Objective
o Select significant point sources from NEDS
data.
o Combine snail stack emission.
o Estimate plume rise from available stack
parameters.
o Subjectively fill in the missing information
-------
PAGE 138
in the original data.
2) Pertinent factors
o Geographical distribution of point sources.
o Relative location to RAMS stations.
o Computational requirements of the model,
o Retain the stack identification parameter
for reference.
3) Analysis procedure
o Tabulate and screen NEDS data (we found
many error and missing cards).
o Subjectively edit original data set. (Mostly
fill in nissing cards or place cards in proper
order).
o Select stacks with SC^ emission and compute
their total emission.
o Estimate DIM coordinate for selected source.
o (lap source locations on area source grid and
coaputational grid for checking their
distribution.
o Select necessary stack parameters from original
data set.
o Compute pluae rise for each stack based on
available stack parameters (various formulas
-------
PAGE 139
were used) .
o Combine saall stacks together to reduce total
number of point sources.
o Estimate values for all the missing data.
o Sort and arrange point-source sequence within
the computational region and compute the
sub-total emission rate.
o Check emission data with 1964-1965 data.
4) Comment
This is most difficult data analysis the authors
have ever been called upon to perform. It is
obvious that a better source inventory is required,
if a model validation is to be considered.
U-C. Surface parameter data
a) Objective
o Estimate surface roughness parameters.
o Estimate effective emission height of
area sources.
b) Pertinent Factor
o The effect of results on the final
-------
PAGE 140
diffusion computation,
c) Analysis procedure
o Screen and tabulate P.M. Vukovich data,
o Convert F.M. Vukovich data froa highway
coordinate to UTM coordinate.
o Convert F.M. Vukovich data from nested
grid to numerical grid.
o Estimate values for vacant region (when
nested grid is converted to stretched grid,
a few strips of the region are without data),
o Estimate silhouette area ratio of each grid element.
o Estimate surface roughness parameters based
on Lettau*s forsula.
o Estimate effective emission height of area
sources,
o Compare with 1964-1965 data.
4-D. RIMS Station location data
The task required for this data set is to convert station
location in longitude and latitude to UTM coordinates and to
compute station locations in terms of the grid unit used by
-------
PAGE 141
the aodel.
4-E- Validity check of all the data
The aboye-analyzed data are coibined together in order to
•eet the input specifications of the aodel. The aodel
computation was performed based on this geographical and
annual eaission data. The 1965 Meteorological data was used
in these coaputations. The computed results were checked
with our previous aodel results and observed data at ten
aonitorforg stations in Feb. 1965. This was done to check
the consistency of the current data set with respect to
aodel foraulation. No validation analysis was atteapted.
4-F. Function of each auxiliary prograa
As aentioned previously, the objective of these auxiliary
prograas is to prepare the data for aodel input. The final
result of data analysis is to create a data set
"EPAGE02.DATA." This data set is read into the aodel by I/O
uait^KUHITG. The scheaatic diagraa shown in Figure II-4
describes how this data set is created by these auxiliary
prograas. In this diagraa, IXIZX.F designates a FORTRAN
prograa and XXXXX.D is a input or output data. The listing
of these auxiliary prograas are given ia Appendix 2.
-------
PAGE 142
(St. Louis Map)
'I I 1
EPAPT.D EPACLKTB.D EPASTA1.D
| j
f EPAPT.F j f EPACLKTB.F j
1 1
F
EPAPT1.D EPACLK.D -+•( EPASTA.F J
I I
EPAPT1.F Y« ! EPASTA2.D
|
^m.n-r«r, _^ ^nAr,T, ,- "N
1 1
EPAPT3.F ~N
EPAPT4.F 7 Prnt 1
1 ^ 1
/- X 1
Prnt >l EPAPT5.F V- »• EPAPT3.D |
I !
EPAPT4.D 1
1
( Ed t ) '
EPAPT5.D |
j
^ EP APT6. F ").« 1
1
EPAPT6.D
^
| |
EPAAREAO.D EPASFCZ1.D
1 EPASFCZ2.D -,
rEPAAREAl Fj \
T 1 EPASFCZ F J
EPAAREA1.D |
i
EPASFCZ3.D
1
EPAAREA2.F J
|
EPAAREA2.D
EPASFCZ4.D
EPASFCZ5.D
*Note:
T XXXXX.F J
Fortran Program
XXXXX.D
Data
Figure 11-4. Schematic Diagram of Processing NEDS & Geographical Data
-------
PAGE 143
1)
Program
Input
Output
Function
EPAPTBL.FORT
NEDS point source data
Print only
a) Tabulate NEDS point source data.
b) Identify missing or misplaced cards in
data set.
c) The result of tbis information is used
to edit (subjectively) original data.
2)
Program
Input
Output
Function
EPAPT.FOBT
EPAPT.DiTJ (edited NEDS point source data)
EPAPT1.DATA
a) Tabulate data.
b) Create new data set EPAPT1 DATA, which
eliminate stacks with zero 302 emission
and only contains selected stack parameters.
3)
Program
Input
Output
Function
EPACLKTB.FOBT (EPA-provided program)
EPACLKTB.DATA
EPACLK.DATA
Copy UTM conversion table from sequential
data set to a direct-access device.
Program
Input
Output
Function
EPAPT1.FOBT
EPAPT1.DATA
EPACLK.DATA
EPAPT2.DATA
a) Convert UfH coordinate of stacks in
zone 16 as extension of zone 15.
b) compute total emission rates from
all point source in tons/year.
5)
Program
Input
Output
Function
Comments
EPAPT2.FOBT
EPAPT2.DATA
Print only
Compute plume rise from each stack by
various formulae.
Current program included five formulae.
They are:
a) Hoses and Carson (Beference 11).
b) Hodified Csanady (Beference 12),
c) CONCAHE (Beference 13),
d) Modified COHCAWEtl (Beference 12).
e) Modified COMCAWEI2 (Reference 12).
6)
Program
Input
Output
EPAPT3.FOBT
EPAPT4.FOBT
BPAPT2.DATA
Print only
-------
PAGE
Function
Comment
a) Hap point source onto 1x1 km2 grid
for whole region in units of tons/year
and gs/sec,.
b) Map point source location onto numerical
grid.
c) Coapute point source emission per each
numerical grid element.
d) List point sources located outside the
region.
e) Compute total point source emission in the
map region and its ratio to original total
emission rate.
f) Computes total number of point source in
the map region.
The information resulting from this program is
used to:
a) select computational region.
b) help in writing program "EPAPTS.FORT.
7)
Program
Input
Output
Function
Comment
EPAPT5.FOHT
EPAPT2.DATA
EPAPT3.DATA
EPAPT4.DATA
a) Combines small-source emission stacks within
a plant to a large stack.
b) Conpute normalized plume rise by modified
CONCAWEI1 Formulae.
c) Create new data EPAPT3,DATA, which has same
format as input data set EPAPT2.DATA,
d) Create new data EPAPT4.DATA, which contains
only x,y,z location of the stack, S02
emission rate and normalized plume rise,
However, the stacks for which data are
missing are identified.
a) EPAPT4.DATA is used to create EPAPT5.DATA,
by terminal input to fill in all the
missing data.
b) The criteria for combining stacks together
are:
o all stacks must be within a plant.
o stack emission rate shall be less than
150 tons/year.
o if combined emission rate of small stacks
is greater than 150 tons/year, then treat
as an independent stack.
o if it is less than 150 tons/year, then
it is combined to other stack which has
emission greater than 150 tons/year.
o weighted average method is used to
computed average physical stack height
and plume rise for combined stack.
The weighting factor is the
emission rate.
-------
PAGE 145
8)
Program
Input
Output
Function
Consent
EPAPT6.FOHT
EPAPT5.DATA
EPAPT6.DATA
a) Arranges and prints point source sequence in
a particular order for analysis purpose.
This program provides four methods of
ordering point source sequence. These
methods are according to source strength;
x,y locations; regions, and the distance
to the center of computational region.
b) Creates EPAPT6.DATA, which eliminates
sources outside the computational region
and arranges point source sequence in
a particular order.
Format for EPAPT6.DATA is put on the first card
of the data set.
Program
Input
Output
Function
Comment
EPASTA.FOBT
EPACLK.DATA
EPASTA1.DATA
-------
PAGE 146
Output
Function
EPAAREA2.DATA
EPASFCZ3,DATA
EPASFC25.DATA
a) Reads in area sources in 200x140 uniform
grid elements from EPAA8EA1, and chooses
40x60 desired area sources, and converts
40x60 to 30x40 area sources corresponding
to computational non-uniform grid
elements.
b) Beads in terrain height and building height
data in Bukhovich's highway coordinates,
and converts them to UTM coordinates,
c) Fills up the undefined area by interpolation
and extrapolation, and converts data from
his nested grid system to uniform grid.
d) Convert the data from uniform grid to
stretched grid for desired region.
e) Combines the terrain height, building
height, and area sources to estimate the
surface roughness.
13}
Program
Input
Output
Function
EPAGE01.FOBT
EPAPT6.DATA
EPASTA2.DATA
EPAABEA2.DATA
EPASFCZ4. DATA
EPASFCZ5.DATA
EPAGE02.DATA
Creates EPAGE02. DATA, which meets the input
specification of model requirement for
geographical and annual average emission
data (see section II-3-B)
14)
Program
Input
Output
Function
EPAGEOIN.FORT
EPAGE02.DATA
print only
Prints out EPAGE02.DATA in map form for
checking purposes.
15}
Program
Input
Output
Function
Consent
EPAMAP.FOBT
EPAMAP2.DATA
print only
prints out St. Louis map on line printer.
EPAMAP2.DATA is a data set which contains a
digitized information of St. Louis map.
-------
PAGE
REFERENCES
1. Shir, C.C. and L. J, Shieh, 1974: A Generalized Urban Air
Pollution Model and Its Application to the Study of S02
Distributions in the St. Louis Metropolitan Area. J.Appl.
Meteo, 13, pp. 185-204.
2. Panofsky, H-A., A.K. Blackadar. and G.E. HcVehil, 1960: The
diabatic wind profile. Quant. J. Roy. Meteorol. Soc. Vol.
86, pp 390-398.
3. Shir, C.C. 1973: A preliminary numerical study of
atmospheric turbulent flows in the idealized planetary
boundary layer, J. Atmos. Sci., Vol. 30, pp 1327-1339.
4. Businger, J.A., et al., 1973: Flux-profile relationships in
the atmospheric surface layer. J. Atmos. Sci., vol. 28, pp
181-189.
5. Pasguill, F., 1974: Atmosgheric. Pit:fjig.JQ.a. London, Van
Nostrand.
6. Lettan, H.H., 1970: Physical and meteorological basis for
mathematical models of urban diffusion process. Proc. symp.
Multiple-Source Urban Diffusion Models, A.C. Stern, Ed.,
APCO Publ. No. AP-86, Research Triangle Park, N.C.
7. U.S. Department of Army, 1957: Universal Transverse Hercator
(»rid« U.S. Dept. of Army Washington, D.C., Publication No.
TMS-247-8.
8. Hitchmyer, H.D., and K.W. Morton, 1967: Difference Methods
for Initial Vaj-ue Problems. New York, Interscience, 405
pp.
9. U.S. Environmental Protection Agency, 1973: GajLde for
Compiling a Comprehensive Emission Inventory (Revised).
Office of Air Quality planning and Standard, EPA,
Publication No. APTD-1135, Research Triangle Park, N.C.
10. Vukovich, P.M., 1974: private communications.
11, Moses, H and 3.E. Carson, 1968: Stack Design Parameters
Influencing Plume Rise. J. Air Poll. Control Assoc,, 18 pp
454-457.
12. Thomas, F.H., S.B. Carpenter and H.C. Colbaugh, 1970: Plume
Sise Estimates for Electric Generating Stations. J, Air
Poll- Control Assoc., 20, pp 170-177.
13. Brummage, K.G. et. al. 1966: Tfeg Calculation of Atmospheric
Dispersion from a S_tack, Stitching COSCAHE, The Hague, The
Netherlands, 57 pp.
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/4-75-005a
3. RECIPIENT'S ACCESSION-NO.
4. TITLE AND SUBTITLE
DEVELOPMENT OF AN URBAN AIR OUALITY SIMULATION MODEL
WITH COMPATIBLE RAPS DATA: VOLUME I
5. REPORT DATE
May 1975
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
C. C. Shir and L. J. Shieh
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
IBM Research Laboratory
San Jose, California 95193
IBM Scientific Center Palo Alto, Calif. 94304
10. PROGRAM ELEMENT NO.
1AA003
11. CONTRACT/GRANT NO.
68-02-1833
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Protection Agency
Environmental Research Center
Research Triangle Park, N.C. 27711
13. TYPE OF REPORT AND PERIOD COVERED
Final Report 7/1/74-5/30/75
14. SPONSORING AGENCY CODE
15. SUPPLEMENTARY NOTES
Issued as Volume I of 2 Volumes
16. ABSTRACT
An advanced generalized urban air quality model (IBMAQ-2) is developed based on the
theory utilized in an existing model (IBMAQ-1) as prescribed in Ref. 1. The model,
'" based on numerical integration of the concentration equation, computes temporal and
three-dimensional spatial concentration distributions resulting from specified
urban point and area sources by using NEDS (National Emission Data System) and
simulated RAMS (Regional Air Monitoring System) data. The UTM (Universal Trans-
verse Metric) coordinates are used in all geographical, source emission, and
monitoring data. A new method to incorporate point sources 'into the grid computa-
tion is developed by using a Lagrange trajectory method. Many model options are
provided which enable users to study conveniently the significant effects, which
these options have on the final concentration distributions.
The program description is included to provide a guide for users. The program is
constructed in a modular form which allows users to change or improve each compo-
nent conveniently. The input auxiliary model,.which processes geographical, source
emission, and monitoring data, is also included.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
Air Pollution
Boundary Layer Modeling
Mathematical Model ins
8. DISTRIBUTION STATEMENT
UNLIMITED
19. SECURITY CLASS (ThisReport)
UNCLASSIFIED
21. NO. OF PAGES
20. SECURITY CLASS (This page)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (9-73)
-------
INSTRUCTIONS
1. REPORT NUMBER
Insert the EPA report number as it appears on the cover of the publication.
2. LEAVE BLANK
3. RECIPIENTS ACCESSION NUMBER
Reserved for use by each report recipient.
4. TITLE AND SUBTITLE
Title should indicate clearly and briefly the subject coverage of the report, and be displayed prominently. Set subtitle, if used, in smaller
type or otherwise subordinate it to main title. When a report is prepared in more than one volume, repeat the primary title, add volume
number and include subtitle for the specific title.
5. REPORT DATE
Each report shall carry a date indicating at least month and year. Indicate the basis on which it was selected (e.g., date of issue, date of
approve!, date of preparation, etc.).
6. PERFORMING ORGANIZATION CODE
Leave blank.
7. AUTHOR(S)
Give name(s) in conventional order (John R. Doe, J. Robert Doe, etc.). List author's affiliation if it differs from the performing organi-
zation.
8. PERFpRMING ORGANIZATION REPORT NUMBER
Insert if performing organization wishes to assign this number.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Give name, street, city, state, and ZIP code. List no more than two levels of an organizational hirearchy.
10. PROGRAM ELEMENT NUMBER
Use the program element number under which the report was prepared. Subordinate numbers may be included in parentheses.
11. CONTRACT/GRANT NUMBER
Insert contract or grant number under which report was prepared.
12. SPONSORING AGENCY NAME AND ADDRESS
Include ZIP code.
13. TYPE OF REPORT AND PERIOD COVERED
Indicate interim final, etc., and if applicable, dates covered.
14. SPONSORING AGENCY CODE
Leave blank.
15. SUPPLEMENTARY NOTES
Enter information not included elsewhere but useful, such as: Prepared in cooperation with, Translation of, Presented at conference of,
To be published in, Supersedes, Supplements, etc.
16. ABSTRACT
Include a brief ^200 words or less) factual summary 'of the most significant information contained in the report. If the report contains a
significant bibliography or literature survey, mention it here.
17. KEY WORDS AND DOCUMENT ANALYSIS
(a) DESCRIPTORS - Select from the Thesaurus of Engineering and Scientific Terms the proper authorized terms that identify the major
concept of the research and are sufficiently specific and precise to be used as index entries for cataloging.
*••
(b) IDENTIFIERS AND OPEN-ENDED TERMS - Use identifiers for project names, code names, equipment designators, etc. Use open-
ended terms written in descriptor form for those subjects for which no descriptor exists.
(c) COSATI FIELD GROUP - Field and group assignments are to be taken from the 1965 COSATI Subject Category List. Since the ma-
jority of documents are multidisciplinary in nature, the Primary Field/Group assignment(s) will be specific discipline, area of human
endeavor, or type of physical object. The application(s) will be cross-referenced with secondary Field/Group assignments that will follow
the primary posting(s).
18. DISTRIBUTION STATEMENT
Denote releasability to the public or limitation for reasons other than security for example "Release Unlimited." Cite any availability to
the public, with address and price.
19. & 20. SECURITY CLASSIFICATION
DO NOT submit classified reports to the National Technical Information service.
21. NUMBER OF PAGES
Insert the total number of pages, including this one and unnumbered pages, but exclude distribution list, if any.
22. PRICE
Insert the price set by the National Technical Information Service or the Government Printing Office, if known.
EPA Form 2220-1 (9-73) (Reverse)
------- |