PB83-200626
COMPLEX/PPM Air Quality Model, User's Guide
Environmental Research and Technology, inc.
Concord, MA
Prepared for
Environmental Sciences Research Lab,
Research Triangle Park, NC
May 83
^BH» •WMBBHPWPRPRPHI ^^w ^IWPPPHHWHi w™
Todniuf Mormttion Senrfce
-------
IPA-600/8-83-015
May 1983
PB83-200626
COMPLEX/PFM AIR QUALITY MODEL
User's Guide
by
D.C. Striroaitis
J.S. Scire
A. Base
Environmental Research & Technology, Inc.
Concord, Massachusetts 01742
Contract No. 68-02-2759
Project Officer
John F. Clarke
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, NC 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NC 27711
HFRGDUCED BT
NATIONAL TECHNICAL
INFORMATION SERVICE
US OEPMINEII! OF COHHEtCE
O. M. inti
-------
TECHNICAL REPORT DATA
(Please read Inunctions on the reverie lit fore campieiin/f)
REPORT NO.
HPA-600/ 8-83-015
3. RECIPIENT'S ACCESSION»NO.
REPORT SATE
Hay 1983
00626
TITLE ANQ SUBTITLE
COHPLEX/PFM AIR QUALITY MODEL
User's Guide
i, PERFORMING ORGANIZATION CODE
AUTHORIS!
D. G. Strimaftis, J. S. Sci're and A. Bass
8. PERFORMING ORGANIZATION REPORT NO.
PERFORMING ORGANIZATION NAME AND ADDRESS
Environmental Research and Technology, Inc.
Concord, Massachusetts 01742
10. Pr
-------
NOTICE
This deuusciit has been reviewed in accordance with
U.S. Environmental Protection Agency policy and
approved for publication. Mention of trade names
or commercial products does not constitute endorse-
ment or recommendation, for use.
11
-------
PREFACE
This publication contains a user's guide for the COMPLEX/PFM disper-
sion model which applies potential flow theory for calculating the impact
of non-reactive pollutants on the windward aide of the first significant
terrain feature downwind of a point source. The potential flow model
offers a realistic physical description of the interaction of plumes
with terrain features; plume path and dispersion coefficients are mod-
ified based on the height of the plume relative to the height and shape
of the terrain feature.
The PFM option within COHPLEX/PFM requires twice-daily representative
radiosonde data as input and incorporates an algorithm to interpolate a
temperature sounding for each hour of the day. The interpolated sounding
is used to calculate the mixing height, plume rise, and critical dividing
streamline for that hour. Plumes above the critical streamline are
assumed to pass over the terrain feature and, for neutral and stable
stratification, potential flow theory is the basis for specifying the in-
teraction of the plume with the terrain feature. Plumes below the cri-
tical streamline are assumed to Impact on the terrain. Since potential
flow theory does not apply for plume impaction, concentrations are calcu-
lated as in the COMPLEX 1 dispersion model.
Potential flow theory has not been extensively applied to dispersion
problems. The COHPLEX/PFM model is being made available through this
document and the UNAMAP system to encourage its application for the
purpose of accruing evaluation statistics for diversified complex terrain
applications. No regulatory standing should be ascribed to COMPLEX/PFM.
Its use for any regulatory application should be considered in light of
EPA's Guideline on Air Quality Models.
The model represents a first step towards a refined complex terrain
model and may be modified in the future based on feedback from users and
on the results of the EPA multi-year complex terrain model development
project. Corrections and modifications of this model may be obtained as
they are issued by completing and sending the form on the last page of
this guide. In addition, the mailing form may also be used to obtain a
sensitivity analysis comparing concentration, mixing height, and plume
rise calculations by COMPLEX/PFM wluh those for COMPLEX I and COMPLEX II.
The report should be available by October 1983.
The computer code for COMPLEX/PFM will be available on version 5 of
the User's Network for Applied Modeling of Air Pollution (UNAMAP) system
by spring 1983. The UNAMAP system may be purchased on magnetic tape from
NTIS for use on the user's computer system. For information on UNAMAP
contact: Chief, Environmental Operations Branch, MD-80, U.S. Environmental
i-rotection Agency, Research Triangle Park, NC 27711.
ill
-------
ABSTRACT
A user's guide has been assembled to describe the purpose, design,
and operation of the COHPLEX/PFM air quality modeling system. The system
combines the features of the Potential Flow Model (PFM) with those of the
EPA COMPLEX I and COMPLEX II models to produce a potential flow complex
terrain model for routine application.
Potential flow dispersion calculations using the theory cf Hunt and
Mulhearn (1973) and Hunt and Snyder (1978) may be selected as an option
within COM?LEX/?FM. When this option is selected, the model requires
hourly wind speed and temperature profiles in order to calculate hourly
mixing heights, hourly plume rise (using a layered plume rise equation),
and hourly values of the critical dividing streamline height. A
preprocessor is provided to interpolate hourly profile data from twice
daily (morning and evening) radiosonde data. Potential flow calculations
are performed whenever the plume lies above the dividing streamline
height, providing that the surface stability class Is D, E, or F.
COMPLEX I (22.5° sector averaging) calculations are performed whenever the
plume lies below the dividing streamline height, and COMPLEX II
calculations are performed whenever the surface stability class is A, B,
or C.
iv
-------
CONTENTS
Preface ill
Abstract iv
Figures . vii
Tables ix
Abbreviations and Symbols xi
Acknowledgement xiii
1. Introduction ......... 1
1.1 Objective ...... ............ 1
1.2 Major Features. ......... ... 2
1.3 COMPLEX/PFM Modeling Package 4
1.4 Summary of DatA Required for PFM Options » 6
1.4.1 Meteorological Data. ......... 6
1.4.2 Source Data. ................... 6
1.4.3 Receptor Data. ....... . 6
1.4.4 Program Control Parameters 7
2. Technical Description. ............ 8
2.1 Introduction. ..., ....... 8
2.2 The Potential Flow Model 8
2.2.1 Diffusion Theory 8
2.2.2 PFM Factors 21
2.2.3 Plume Trajectory Analysis - Neutral Flow ..... 22
2.2.4 Plume Trajectory Analysis - Stratified Flow. ... 28
2.3 The COMPLEX/PFM Option ............. 31
2.3.1 PFM Option: Long and Short. 32
2,3.2 Plume Rise ........ ....... 36
2.3.3 Froude Number and HC Analysis. 41
2.4 Temperature and Velocity Profile Interpolation. ..... 44
2.4.1 Required Data and Defaults 44
2.4.2 Mixing Heights and Profile Interpolation ..... 44
-------
2.5 Terrain Description Guidance. 48
2.5.1 PMF-Short Obstacle Selection ... 48
2.5.2 PHF-Long Obstacle Selection 50
3. User's Instructions. .............. .. 55
3.1 SETUP Instructions. . 55
3.2 READ56 Instructions , . 55
3.3 PROFIL Instructions ........ 59
3.4 COMPLEX/PFM Instructions 61
3.4.1 COMPLEX/PFM Options. 61
3.4.2 Input Format Specifications. . . 63
References., 84
Appendices
A. PFM Bluff Cross-Sections 86
B. Program Description - PROFIL. ..... 89
C. Program Description - COMPLEX/PFM ....... 91
D. Subroutine Description - SPFM . , 95
E. Test Case 101
vi
-------
FIGURES
Number Page
1 The COMPLEX/PFM modeling system ..... 5
2 Horizontal dispersion coefficients for neutral flow
over a circular ridge ..». 16
3 Vertical dispersion coefficients for neutral flow
ever a circular ridge . . . 17
4 Horizontal dispersion coefficients for neutral flow
over a hemispherical hill 18
5 Vertical dispersion coefficients for neutral flow
over a hemispherical hill ............. 19
6a Velocity speed up factors for neutral flow over a
circular ridge. 20
6b Velocity speed up factors for neutral flow over a
hemispherical hill. ........ 20
7 Plume path and surface streamlines derived from
flow over a step. 25
8 Height of the source streamline above the hill
crest for various stack heights and Froude Numbers. 29
9 Relationship between PFM and actual sources and
receptors in PFM-Long 35
10 Relationship between PFM receptors and COMPLEX
receptor in PFM-Short 37
11 Calculation of neutral plume rise . . 40
12 Relationship between PFM geometry and plume height
when HC, the dividing streamline height, is
not zero 43
13 Sample construction of temperature profile for fourth
hour after A.M. sounding 46
vii
-------
FIGURES
Page
Relationship between north and the orientation of
the bluff coordinate system ............ 49
15 Definition sketch for selecting cross-wind and
a long-wind aspect ratios 51
16 Relationship between XO, YO, and OBSANG 52
17 The COMPLEX/PFM modeling system 56
18 Input data deck setup for COMPLEX/PFM ........ 64
viii
-------
TABLES
Number
1 Obstacle shapes available far applications of
PFM-Long 53
2 PARAMS namelist - READ56 57
3 Data identification record - READ56 58
4 FROFIL program control card variables. ....... 60
5 Options available in OOMPLEX/PFM .......... 62
6 COMPLEX/PFM card type 1, 2, and 3 - title. 65
7 COMPLEX/PFM card 4 - control and constants ..... 66
8 COMPLEX/PFM card 5 - options 67
9 COMPLEX/PFM card 6 - wind and terrain . 68
10 COMPLEX/PFM card type 7 - point source ....... 69
11 COMPLEX/PFM card type 8 - specified significant
sources 70
12 COMPLEX/PFM card type 9 - met. data identifiers. . . 71
13 COMPLEX/PFM card type 10 - polar coordinate
receptors. .................... 72
14 COMPLEX/PFM card type 11 - polar coordinate
receptor elevations 73
15 COMPLEX/PFM card type Iz - receptors ........ 74
16 COMPLEX/PFM card type 13 - obstacle information. . . 75
17 COMPLEX/PFM card type 14 - grid center 76
18 COMPLEX/PFM card type 15 - receptors ........ 77
ix
-------
TABLES
Number Page
19 COMPLEX/PFM card type 16 - segmented run ...... 78
20 COMPLEX/PFM card type 17 - meteorology 79
21 COMPLEX/PFM optional input file - meteorological
data ......... 81
22 COMPLEX/PFM optional input file - emission data. . . 82
23 COMPLEX/PFM optional input file - hourly sounding
data 83
-------
LIST OF SYMBOLS AND ABBREVIATIONS
SYMBOL
a height of obstacle
a,b,c half the length of the primary axes of an ellipsoid
C concentration
C1 dimensionless concentration
D depression of a neutral streamline to approximate stratified
flow
Di,2 normal and crosswind dimensionless diffusivities
Dl,2 normal and crosswind "effective" diffusivities
fl,2 trial solution functions
F buoyancy flux
FR Froude number
g(s) trial solution function
h step height (bluff problem)
H plume height above flat terrain
HC critical dividing streamline of the flow
HE effective source height (above HC)
HH pet plume height above the surface
normal dimensional diffusivity
crosswind dimensional diffusivity
H local neutral streamline height above the obstacle
L streamline lift near an obstacle
LJI streamline lift near an obstacle in neutral flow
n dimensionless vertical coordinate normal to plume trajectory
ns distance along normal from plume centerline to surface
N Brunt-Vaisala frequency
Q emission rate of pollutant
R(s) distance from axis of symmetry to plume centerline for
axisymmetric obstacle
s path length along plume centerline
s,n,y curvilinear coordinates for system that follows the plume
centerline
S stability parameter
S speed-up factor
t advection tiina along plume centerline
t parametric complex variable in bluff streamline equations
T(s) line integral for crosswind flow deformation
un local velocity along the plume in neutral flow
us local velocity along the plume in stratified flow
xi
-------
-u(s) scaled along-streamline velocity
U(s) along-streamline velocity |
H» uniform velocity upwind of source (u-component)
VCD uniform velocity upwind of ;source (v-componont)
v(s) scaled velocity normal to streamline
JV(s) velocity normal to streamline
x distance along X axis from source
Xg downwind distance from soutce to obstacle center
ZG height of streamline above crest for neutral flow
ZCF height of streamline above crest for stratified flow
fl' entrainment parameter
r> plume height perturbation factor
A ellipsoidal coordinate
p density
Gy,z crosswind and vertical dispersion coefficients
CTyf,zf crosswind and vertical dispersion coefficients in flat
terrain
<(>(s) line integral for normal (vertical) to streamline flow
deformation
4 velocity potential
i|J0 stream function along surface
\})8 stream function through source
ij»p streamline equation of plume centerline
ABBREVIATIONS
cm centimeter
EPA U.S. Environmental Protection Agency
ERT Environmental Reasearch & Technology, Inc.
FMF Fluid Modeling facility
K Kelvin
km kilometer
m meter
NCC National Climate Center
NWS National Weather Service
PGT Pasquill-Gifford-Turner
s seconds
xil
-------
ACKNOWLEDGEMENT
The authors wish to acknowledge a.my helpful discussions with Dr. A.
Venkatram as well as the technical contributions of Dawna Paton, John
Beebe, Anthony Curreri, and Dr. John Nuber in the construction and testing
of many of the algorithms contained in this model. Special thanks are
extended to John F. Clarke for his timely suggestions and overall
direction on the project and also for his perseverance in identifying
many of the final bugs in the program and preparing the test case
computations.
XLll
-------
SECTION 1
INTRODUCTION
I.I Objective
The focus of the development work that produced the modeling system
described in this manual has been the application of potential flow theory
to the problem of predicting pollutant concentrations on significant
terrain features. Two classes of meteorological conditions are often
associated with the likelihood of large ground-level concentrations: (1)
low wind speed and stable cases and (2) moderate or high wind speed and
neutral or slightly stable cases. The first class of conditions generally
leads to high concentrations through direct plume impingement or terrain
blocking. The second class of situations promotes high concentrations
because, as the plume is transported over terrain features, it is forced to
pass close to the terrain surface. Physical mechanisms relevant to this
class of conditions include terrain-induced alteration of the plume
centerline trajectory and kinematic constraints on horizontal and vertical
dispersion. Potential flow theory is particularly useful in modeling
concentrations produced by the second class of conditions.
Other methods have been used to account for plume trajectory
deformations over terrain features. For example, those EPA models
suggested for use in complex terrain on a site-specific basis (COMPLEX I,
COMPLEX II, and MPTER, all available on UNAMAP-Version 4) modify the height
of the plume centerline over the terrain. This is done by using plums path
coefficients that vary according to the dispersion stability class of the
surface layer. These coefficients do not vary with obstacle shape or
position along the plume path and are only weakly responsive to the
magnitude of the speed of the flow relative to the density stratification.
Consequently, a better representation of plume dynamics near terrain is
sought through the application of potential flow theory.
The first phase in the study of the potential flow theory approach
culminated in the construction and testing of what has become known as the
Potential Flow Model, PFM (Bass et al. 1981). That version of PFM
contained potential flow streamline solutions for flows over the center of
an isolated hemisphere and an isolated circular cylinder. Flow over
obstacles of intermediate shape between these two was approximated by a
streamline weighting scheme. The model was designed to estimate pollutant
concentrations from a single plume over a single terrain feature only when
the atmosphere's density stratification is neutral to slightly Ltable (that
is, no blocking or impingement).
-------
A more general model is needed, however, to apply potential flow
theory calculations in regulatory permitting situations. Concentration
estimates are needed when the flow is either very stable or well-mixed
(unstable) as well as when it is neutral to slightly stable. In addition,
several sources often need to be modeled simultaneously for sequential
hourly periods totaling as much as one year so that three hour, 24- ho'ir,
and annual pollutant concentration averages at many receptors can be
estimated.
Potential flow calculations themselves must also be generalized. Hill
shapes other than the sphere-cylinder variants are needed to better
approximate single terrain features. Moreover, plume trajectories that do
not pass over the center of a terrain feature are necessary, as the proper
simulation of wind direction variability is an essential feature of a
successful sequential model.
The objective of the second phase of this study was the generalization
of Che PPM code, and the design of a complete modeling system that would
satisfy each of these requirements. That modeling system is COMPLEX/PFM.
The PFM calculations performed by this model are applicable to more general
terrain shapes and wind flow angle. And the host model (EPA COMPLEX)
contains the structure for handling many point sources, many hourly
simulation periods, and any averaging time for reporting air polluti»n.t
concentration statistics. The host model also provides pollutant
concentration algorithms for impingement situations and strongly convective
situations.
1.2 Major Features
COMPLEX/PFM is a modified version of COMPLEX I/COMPLEX II which
contains the PFM algorithm as an option. If the PFM option is not invoked,
the model performs COMPLEX I (22.5° crosswind sector averaging and Gaussian
vertical spread based on az) computations for stability classes five
and six (E, F), and COMPLEX II (bivariate Gaussian spread based on cry
and oz) computations for stability classes one through four (A through
D). Both COMPLEX versions are based on the MPTER model (Pierce et al.
1980) with an expanded list of terrain algorithms. Major features of
COMPLEX include:
o Averaging periods of longer than one hour, if selected by user.
o Hourly meteorological data that may be read off punched cards for
each hour, or from a tape (or disk) containing a year's data
(same data as used for RAM or CRSTER).
o Optional terrain adjustments as a function of stability class.
o Inclusion or omission of stack downwash.
o Inclusion of gradual plume rise, or final rise only.
o Inclusion or omission of buoyancy-induced dispersion of pollutant
during plume rise using the Paaquill method.
o Input of anemometer height.
o Input of wind profile power law exponents as functions of
stability.
-------
o Concentration contributions that are available per hour and/or
for the selected averaging period at each receptor from up to 25
sources.
o Concentrations available hourly and/or for the selected averaging
period at each receptor.
o Optional output of the following information.* average
concentration over length of record, plus highest five
concentrations for each receptor for four end-to-end averaging
times (1-, 3-, 8-, and 24 hour), and an additional averaging time
selected by the user.
o Optional output files for further processing of concentrations
that are available per hour and f>r each averaging period.
o Up to 50 spatially separated point sources.
o Up to 180 receptors with no restrictions on location.
Several new features are introduced when the PFM option is invoked.
These include the following:
o Refined hourly mixing height calculation based upon twice daily
radiosonde profiles of wind and temperature.
o Layered plume rise calculations that use hourly wind and
temperature profiles interpolated from twice-daily radiosonde
observations.
o Assessment of wind flow characteristics by means of a critical
dividing streamline (HC) and Froude number (PR) analysis based on
the hourly wind and temperature profiles.
o Selection of the COMPLEX I, COMPLEX II, or PFM algorithm based on
the surface stability class, FR, HC, and initial plume height
instead of on the stability class alone.
o Inclusion of explicit plume size and trajectory deformations
computed from potential flow theory when PFM is applicable.
o An economical long-term version (PFM-Long) for computing annual
average concentrations, for assessing regions of frequent high
concentration and for identifying critical periods of meteorology
that produce the highest expected pollutant concentrations.
o A short-term, worst-case or critical period version (PFM-Short)
for a refined analysis of maximum pollutant concentrations
expected.
Along with these new features come several added restrictions to the
original COMPLEX features:
o all sources (maximum = 50) are at the same point,
o the receptor pattern is radial; and
o there is no gradual plume rise option.
These restrictions arise from the necessity of maintaining a simple, fixed
source-terrain relationship.
The new features should be viewed as a natural extension of the
COMPLEX modeling system when representative profiles of wind and
temperature are used. Characterization of the relationship between HC, FR,
-------
plume height, and terrain height is essential to assessing plume behavior
and to choosing the most appropriate algorithm to simulate it. When plume
impingement is likely, based upon these quantities and not simply upon the
dispersion stability class of the surface layer, a VALLEY-like (Burt 1977)
(COMPLEX I with the recommended parameter choices) computation is
performed. When potential flow theory applies, a PFM computation is
performed. If the boundary layer is convectively unstable, a COMPLEX II
computation is performed with a suitable plume path coefficient.
1.3 COMPLEX/PFM Modeling Package
The COMPLEX/PFM modeling package is schematically illustrated in
Figure 1. The meteorological preprocessing steps on the left side of the
figure are required by most sequential air quality dispersion models in use
today, especially those available through UNAMAP. This preprocessor
requires hourly surface data and twice-daily mixing height data. For
further information about the preprocessor program, consult the User's
Guide for the Single Source (CRSTER) Model (EPA 1977).
The preprocessing steps on the right side of the figure pertain only
to the use of the PFM option within COMPLEX/PFM. For short-term runs,
user-specified meteorological data may be input on cards. Hourly
temperature and wind profiles are required when the PFM option is
invoked. These can be constructed by the PROFIT, preprocessor from
twice-daily profile data and hourly surface data. When on-site tower data
are available, tower wind and temperature data should also be included by
the user although the code for doing this is not yet included in the
preprocessor.
Two sources of twice-daily profile data may be available:
representative profiles obtained at a National Weather Service (NWS)
Upper-Air Network Station, or on-site profiles obtained for aite-specific
model application. The NWS data are obtained from the National Climate
Center (NCC) on their TDF5600 series data tapes. When the NWS data are
used, program READ56 reads the TDF5600 tape, checks for missing data, and
reformats the wind, temperature, and pressure data. At this point, the
user must ensure that missing data it. properly flagged. When on-site data
are used, the user is responsible for formatting the data and flagging
missing or bad data.
All source, receptor, and program control data are entered by card
deck. The control information determines which algorithms are to be used
in the computations and what sort of program output is desired.
Prior to execution, an unformatted look-up table of intermediate PFM
computations must be prepared for use with the PFM-Long option. A
formatted file of this table accompanies the COMPLEX/PFM source code.
Program SETUP (lower right on Figure 1) converts this file to the form
needed at execution time. SETUP must be run only once, when the model is
first compiled on a new system.
-------
Figure 1. The COMPLEX/PFM modeling system.
-------
1.4 Summary of Data Required for PFM Options
Four types of input data are required for the COMPLEX/PFM diapersion
model:
o meteorological data,
o source data,
o receptor data, and
o program control parameters.
Moat of this data is needed by the host model and therefore has been
described at length in the MPTER user's guide (Pierce et al. 1980). Only
those data required to operate the PFM option will be discussed here.
1.4.1 Meteorological Data
The PFM option requires hourly temperature and wind profiles. These
are interpolated from twice daily observations in program PROFIL. When
short periods of time are simulated, the user should analyze these profiles
by hand and create the hourly sequence file directly. The hourly data
should include:
o the number of data levels,
o the mixing height (m),
o the height of each data level (m),
o the wind speed (tn/s) and wind direction from which the wind blows
(degrees from North) at each data level, and
o the temperature at each data level (°K).
1.4.2 Source Data
The PFM option requires the same source data as COMPLEX I and II, but
two restrictions are imposed:
o no more than 50 point sources are allowed and
o all sources must be placed at the same location (as in CRSTER).
1.4.3 Receptor Data
A total of 180 receptors may be used with the PFM option. These must
be aligned along rays from the source* The position and number of
receptors along any ray is independent of the position and spacing of
receptors on other rays. Also, the angular spacing between rays need not
be regular.
The receptor information required for PFM-Long includes:
o coordinates of the origin (source location)(user units),
o receptor name,
o receptor bearing and distance (user units),
o receptor height above terrain (user units),
o the height of terrain above sea level at receptor position (user
unite),
o relief height of controlling terrain feature (user units),
-------
o the distance between the source and the controlling terrain
feature (user units), and
o the obstacle shape code that best describes the controlling
terrain shape.
When PFM-Short is invoked, additional obstacle information is needed:
o source coordinates relative to the obstacle coordinate system
(user units),
o the relief height of the terrain feature (user units),
o the angle of rotation from North of the obstacle coordinate
system (degrees clockwise),
o cross-wind and along-wind obstacle aspect ratios, if the obstacle
is a hill, and
o the blu^f shape code, if the obstacle is a bluff.
1.4.4 Program Control Parameters
COMPLEX/PFM requires one more program control parameter than needed to
run COMPLEX I or COMPLEX II. The PFM option is set by option number 26. A
value of one (1) will invoke the PFM-Long algorithm, whereas a value of two
(2) will invoke the PFM-Short option. If zero is specified, no PFM
computations will be done.
-------
SECTION 2
TECHNICAL DESCRIPTION
2.1 Introduction
COMPLEX/PFM is a modified version of COMPLEX (I and II) which contains
an option for potential flow computations. The potential flow computations
are derived from PFM. Section 2 presents the theoretical background of PFM
and the method of implementing the theory within the framework of the
COMPLEX model system.
Section 2.2 contains a description of the diffusion theory appropriate
for deforming potential flow streamlines and the methods of solution for
potential flows over isolated simple bluffs and hills. It also describes
how computations based on the theory are used to modify the flat terrain
Gaussian diffusion equation. Section 2.3 presents a description of the PFM
option within COMPLEX. Both the long-term and short-term versions are
discussed.
Sections 2.4 and 2.5 deal with data requirements for COMPLEX/PFM. The
need for hourly temperature and velocity profiles requires a new
meteorological preprocessor to interpolate twice-daily profile observations
when more frequent soundings are unavailable. Section 2.4 describes how the
preprocessor interpolates between soundings to produce hourly Bounding
information. Section 2.5 provides guidance in specifying bluff shapes and
ellipsoid aspect ratios to represent real terrain features.
2.2 The Potential Flow Model
2.2.1 Diffusion Theory
The potential flow complex terrain model incorporates the theory
developed by Hunt and Mulhearn (1973) and Hunt and Snyder (1978) for
turbulent plumes embedded within potential flow fields. They developed
solutions to the diffusion equations describing flow fields over
two-dimensional and three-dimensional axisymmetric terrain obstacles.
Qualitatively, these solutions are of Gaussian form, with crosswind and
vertical dispersion coefficients evaluated as line integrals of the velocity
field along the plume centerline trajectory. Because the mathematics used
are not immediately familiar, some details of the derivation are presented;
the reader is urged to consult the original references.
-------
The. modeling technique applies potential flow theory to a Gaussian
point source model. It permits explicit evaluation of the terrain-induced
kinematic constraints and trajectory variations. In adopting the present
technique, an attempt has been made to broaden its applicability.
Approximations have been developed to extend the neutral formulation to
slightly stable cases and to allow for three-dimensional terrain obstacles
that are not axisymmetric. This approach is suggested for the following
topographic and meteorological situations:
o isolated, simple bluffs of arbitrary height and slope and simple
hills of arbitrary height and arbitrary aspect ratio in the
cLOSSwind and downwind directions, and
o neutral to slightly stable density stratifications.
At the outset, the limitations of the potential flow model should be
stressed. These limitations arise mainly from physical effects that are not
described by the model:
o The presence of a realistic surface boundary layer is ignored.
o Relevant physical phenomena such as velocity shear, radiative
heating, and flow separation are omitted.
o The theoretical model requires that tbe plume remain "thin"
compared to its height above the terrain. This criterion is often
violated near the hill crest.
[Strictly speaking, the first two of these limitations also restrict the
theoretical validity of a conventional Gaussian solution in flat terrain
situations (Pasquill 1978), unless the set of dispersion coefficients
specifically accounts for these effects.] It is prudent to apply the
potential flow theory approach only to the windward side of obstacles and
not to situations dominated by lee wake or separation effects not treated by
potential flow theory. In many cases, the computed plume dimension will
violate the restriction to "thin plumes", but it is not unreasonable to
"stretch" the theoretical formulation in these cases (Hunt, Puttock, and
Snyder 1979, Bass et al. 1981).
Two~Diroensional Obstacles
The starting point for the calculation is the two-dimensional advection
diffusion equation in the streamline coordinate system (s,n,y):
f \ OC OV 0\> ~ / \O V »»/\(7*-* /n-iV
u(s) •=— + n •=- T— = D,(s)—— + D,{s)—— (2-1.)
Here s, n, and y are dimensionless directional coordinates in the along
plume, normal (nearly vertical), and crosswind directions relative to the
plume trajectory. When the terrain is flet and the trajectory is level, the
3 and n directions correspond to the downwind distance, x, and the vertical
distance from plume centerline, z. In Equation 2-1, u and v are velocities
along and normal to the trajectory, and DI(S) and D2
-------
dif fusivities in the n and y d_;-ections that vary along the trajectory. C*
is the dimension less concentration. The variables in Equation 2-1 are
nondimensionalized by the height of the obstacle, a, and the uniform
velocity upwind of the source, Ifeo, so that:
C' C UOQ a2/Q (2-2a)
u(s) = U(s)/UU (2-2b)
v(s) - V(g)/Uko (2-2c)
D^s) = K1(s)/IaUco] (2-2d)
D2(s) » K2(s)/[aUtoJ (2-2e)
where Q is the source strength (mass/unit time), U(s) and V(s) are the
dimensional velocities (length/unit time), and Kj(s) and K2(s) are the
dimensional diffusivities (length squared/unit time). The relationship
between these- dif f us ivities and the PGT dispersion cc fficients (Turner
1970) will be discussed later.
Substituting the continuity equation:
fe * fe ' ° (2-3>
into Equation 2-1 yields:
t \ 9C' 3" 3C' r> I llfil ,. n /• 1 J2c' {-, A)
u(s) -5 — - n -5- -5 — - D (a) — — + D (s) — . (2-4;
3s 3s 3n 1 2 2 2
Solutions to Equation 2-4 are sought that are Gaussian in form:
C'(S,n,y) = ~ exp [ - f^s) n2 - £2 (B) y2] (2-5)
and also satisfy the additional constraint that the mass flux across the
plume is constant, i.e.:
+00 -KO
J J u(s) C'dn dy - 1.0 (2-6)
10
-------
This constraint requires that:
1/2
g(s) = TT u(s) / [fL(s) £2(s)] (2-7)
Substituting Equation 2-5 into 2-4 yields:
lu(s)]2/ [4$(s) = J D1(s1) u (s')ds'
and
T(s) = J ^(s'VuCs'MdB1 (2-103
o
By Equation 2-7:
g(s) = 4ll[(})(s) T(s)]1/2 (2-11)
In more familiar terms, the dispersion coefficients in the Gaussian form of
the diffusion equation solution, crn and Cfy, are given by:
f.(s) = 1/[20 2J (2-12a)
1 n
f2(s) - l/[2ay2] (2-12b)
11
-------
or
an2 - l/[2 f^s)] = 2]2 (2-13)
ay2 = l/[2 f2(s)J = 2T(s> (2-14)
Three-Dimensional Axisymmetric Obstacles
The arguments for three dimensions proceed much the same as in the
previous section. The appropriate non-dimensional diffusion equation is:
u(s) ££l + n &• 2£1 = Di(s) JL£1 + D2(s) -&£. (
3s 8n In 3n2 y2
where the coordinate system is now with respect to (s,n,y). Here, y ^-s
the angular, azimuthal coordinate from the axis of symmetry along the flow
with:
y * Y^(s) (2-16)
where R(s) is the distance from the axis of symmetry of the obstacle to the
plume centerline, and y is the crosswind distance from the plume
centerline. The continuity equation in this coordinate system becomes:
] . 0 (
n s ds
Tne trial soluticn sought is of the form:
C'(s,n,Y) = ~r exp[-f1(s) n2 -f2(s)y2] (2-18)
subject to the constraint:
2TT •*»
' _l u(s)C'(s,n,Y)R(s) dndy = 1.0 (2-19)
Substitution of Equation 2-18 into 2-15 yields:
f(s) = Iu(8)R(s)]2 /[45>(s)] (2-20a)
12
-------
£2(s) = l/[4l(s)] (2-20b)
where
<|>(B> = J D1(s') [R(8')]2u(8') ds1 (2-21)
s „
T(s) = J [D^(s')/[u(s') R (s^JJds1 (2-22)
2
o
and substitution of Equations 2-18 and 2-20 into 2-19 gives:
g(s) - 4f fo(s) T(s)5/2 (2-23)
The forms of the familiar Gaussian dispersion coefficients are obtained by
substituting Equation 2-20 into Equation 2-12a and b:
a n = 2 <|>(s) / [u(s) R(s)J2 (2-24)
ffy = 2 [R(s)]2 T(s) (2-25)
PGT Scaling of Dispersion Coefficients
Evaluating the terrain-influenced dispersion coefficients (Equations
2-13, 2-14 and 2-24, 2-25) with this formulation requires specifying the
crosswind and the normal diffusivities, D£ and D^. To compare model
calculations with analogous flat terrain situations, an approximation scheme
was implemented, using the PGT dispersion coefficients (Turner 1970) as a
calibrating scale. Other dispersion coefficient systems can be easily
incorporated within the basic methodology used here.
Qualitatively, the diffusivity at a given distance from the source
along the plume centerline streamline is taken as that for the same transit
time in flat terrain. The consequence of this assumption is that model
calculations of dispersion coefficients reduce to flat terrain values in the
limits of large downwind distances or small obstacles.
13
-------
Substituting Equations 2~2l and 2-22 into Equations 2-24 and 2-25 gives
explicit expressions for az and 0y. (From here on, the dispersion
coefficient in the normal direction is denoted as az.) At a distance,
a, from the source along the streamline, the dispersion coefficients are
given by:
2 2 D (s) 8 2
a (s) - A - J R (s») uCs^ds1 (2-26)
S
a2 (s) s 2R2(s) D(s) J ds'/[R2(sf) u(s')] (2-27)
2
o
where D^, D2 are typical mean diffusivity values at distance s.
These values are assumed to be given by the PGT values (ozf ,
Oyf):
dj(s)/2t (2-28)
= o (a)/2t (2-29)
where s is the integrated path length along the plume trajectory and t is
the advection time:
s * J ds1 (2-30)
o
s , i
f — (2-31)
u(s') J '
-------
Substituting Equations 2-28 through 2-31 into Equations 2-26 and 2-27 yields
the final form of the dispersion coefficients:
a
zf
I R2(s')
t u2(s) R2(S)
(2-32)
= a
J
ds'
'(s1) u(s'
(2-33)
The bracketed quantities approach 1 as streamline deformation becomes
negligible, and oz, Oy reduce to their flat terrain values.
Figures 2 and 3 illustrate the plume spread statistics oy and
az, respectively, for neutral flow over a circular ridge. The ridge is
1.0 km high; the point source, 200 m high, is 4.5 km from the ridge center.
Calculations are shown for diffusivities scaled to reproduce the PGT neutral
stability class values for Oy, and oz in the flat terrain case.
This illustrates that only the vertical dimensions of the plume (as
characterized by Oz) differ from flat terrain values in flow over
two-dimensional obstacles. (In each of Figures 2 through 5, the center of
the ridge is indicated by a short vertical bar at the downwind distance 4.5
km; the horizontal bar extending from 3.5 km to 5.5 km denotes the
windward-leeward extent of the ridge.)
As streamlines diverge to the leeward and windward side of the hill,
<7Z increases. Over the ridge, vertical compression of the streamlines
leads to a drastic decrease in oz. Over the ridge crest the vertical
plume thickness is about equal to the height of the plume centerline above
the surface, ns. This marginally obeys the thin plume approximetion.
With the same source-obstacle geometry retained, but the terrain
obstacle type changed from a circular ridge to a hemispherical hill, Figures
4 and 5 display the three-dimensional neutral flow results developed from
Equations 2-24 and 2-25. In this case, both cry and az are affected
by the potential flow field, as evidenced by a marked increase in CTV
over the hill top. Comparing
-------
1000
> 100
o
I
Crest
Hid
••Extent*
I 1
1 (1) Computed Stability Class D (ridge)
—.— 12) PGT Stability Class 0 (flat)
_L
I
1 2 3 4 5 7 10 20
Downwind Distance from Source (km)
Figure 2. Horizontal dispersion coefficients for neutral flow
over a circular ridge.
16
-------
1000
— 100
N
o
10
i
Crest
—U
Hill
••Extent*
_| L
—"— (1) Computed Stability Class D I ridge)
—..— (2) PGT Stability Class 0 (flat)
I
I 2 3 4 S 7 10 20
Downwind Dinancs from Source (km)
Figure 3. Vertical dispersion coefficients for neutral flow
over a circular ridge.
17
-------
1000
(1) Computed Stability Class D (hill)
— (2J PQT Stability Class D (flat)
3 4 5 7 10 20
Downwind Distance from Source (km)
Figure 4. Horizontal dispersion coefficients for neutral flow
over a hemispherical hill.
18
-------
1000
— 100
M
o
I
10
1 •• "" ' (1) Computed Stability Class D thill)
— (2) PGT Stability Class D (flat)
Crest
I
I 2 3 4 5 7 10 20
Downwind Distance from Source (km)
Figure 5. Vertical dispersion coefficients for neutral flow
over a hemispherical hill.
19
-------
(fl
o
•rf
I
D.
3
•a
s
2.0
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1.0
.9
.8
.7
.6
.5
Cylinder
(2D)
1.0 2.0 3.0
Downwind Distance from Source (km)
4.0
Figure 6a, Velocity speed up factors for neutral flow over a
circular ridge.
1.5
1.4
1.3
1.2
1.1
1.0
.9
.8
.7
.6
.5
to
o
8
LL
a
3
•Q
I
(0
Sphere
(3D)
1.0 2.0 3.0
Downwind Distance from Source (km)
4.0
Figure 6b. Velocity speed up factors for neutral flow over a
hemispherical hill.
20
-------
the hemispherical hill. Both streamline patterns upwind and downwind of the
obstacle crests are symmetric; therefore, the speed-up factor has been
plotted only from source to crest. The theoretical speed-up value at the
surface of each of the crests is 2.0 for the cylindrical ridge, and 1.5 for
the hemispherical hill. Because the source streamline does not touch the
surface at the crest, the surface speed-up values are approached, but not
attained, in Figure 6. Note that the influence of the cylindrical ridge
extends farther upstream (and downstream) than does the influence of the
hemispherical hill.
2.2.2 PFM Factors
Th«? discussion of the previous section presented the theory used in PFM
for calculating concentrations on the surface of a simple obstacle, when a
potential flow streamline trajectory had been computed. The method for
scaling diffusivities using sigma-curves (e.g. PGT) was also presented.
This scaling approach provides a simple way of quantifying the plume
deformations predicted by PFM so that "PFM-like" computations can be
simulated within other Gaussian plume models.
The Gaussian plume equation for surface concentrations under unlimited
mixing conditions for level terrain is:
C - 2 C0 exp (-0.5(H/ozf)2) exp (-0.5 (y/ayf)2) (2-34)
where
C0 = Q/(2ir 0yf azf um) (2-35)
and
Q - pollutant emission rate
Urn = wind speed
H =* plume height above the surface
y = lateral distance to plume centerline
O2f = standard deviation of concentration in the vertical (flat
terrain)
Oyf = standard deviation of concentration in the horizontal (flat
terrain)
The corresponding equation in the PFM for a plume near an obstacle is
Equation 2-5 (or 2-18). However, when the functions g(s), fjCs), and
f2(s) are rewritten in terms of u(s), Oz(s), and Oy(s) the form of
Equation 2-5 (or 2-18) is the same as Equation 2-34:
Ch = 2 Coh exp(-0.5(n/az(a))2)exp (-0.5(y/oy(s))2) (2-36)
21
-------
; where "'
C0h = Q/(2ir are computed from potential flow theory within PFM. Ratios of
the sigmas are readily found according to Equations 2-32 and 2-33 and their
properties have already been discussed. For convenience, four PFM "factors"
are defined in terms of these quotients:
o (s) a (s) u{s)
(2-40a)
(2-40b)
(2-40c)
(2~40d)
These factors are independent of wind speed and dispersion parameters. They
depend only on the source-terrain geometry, receptor distance, and the
density (((ratification of the atmosphere. Hence, more complicated
computations such as plume trapping and buoyant plume enhancement may
readily extend the usefulness of PFM computations. In particular, these PFM
factors are central to the PFM option within COMPLEX, as discussed in
Section 2.3.2.
2.2.3 Plume Trajectory Analysis - Neutral Flow
Streamlines are computed for two basic obstacle shapes in PFM - a
two-dimensional bluff and a three-dimensional ellipsoid. In each case, the
shape is kept simple to minimize computer costs. Although few real terrain
features resemble pure ellipsoids or simple bluffs, it should be remembered
22
-------
that potential flow itself ia an approximation to the wind field over the
earth's surface. The intent of the model ia to compute first-order changes
in the plume trajectory and shape as it passes over major terrain obstacles'
Flow Over a Bluff
Potential flow is non-divergent flow of a constant density fluid
without vorticity. The only condition placed on the velocity field of such
a flow near a solid surface is that there should be no velocity component
into the surface. This means that the velocity component normal to the
surface must be zero. The component parallel to the surface may take on any
value. Therefore, the surface must represent a streamline of the flow and
conversely, any streamline of the flow may be considered a solid surface.
This property of potential flow is used to generate streamlines for
flow over a family of bluff shapes. The basic solution contained in PFM is
for flow over a step (Milne-Thomson, 1960). However, by considering any
streamline of this solution as a possible solid surface, the PFM actually
contains an infinite nu^er of bluff shapes. Of the twenty shapes
explicitly incorporated, each streamline has its own characteristic alope
(they range from 0.15 to 3.2). A user only has to select the shape most
like the terrain feature to be modeled. Section 2.5 contains instructions
for selecting the most appropriate bluff shape,
Caution must be used in applying potential flow theory to very steep
bluffs. As the bluff face becomes steeper, the probability of streamline
separation and secondary circulation increases. The potential flow solution
completely neglects these effects and therefore becomes inappropriate. In
fact, the theoretical velocity field close to the crest of the step feature
tends to infinity. Consequently, the maximum bluff slope represented in PFM
is about 3:1.
The solution to the step problem requires complex variables and the
final equations are transcendental functions of a complex parametric
variable, t. Let RT be the real part of t, and MT be the imaginary part of
t; then:
nx/h
iry/h
(sinh RTXsin MT)
* RT + (sinh RT)
-------
In these equations, x and y are the space variables in the horizontal
and vertical directions, respectively; ij> is the streamfunction, a constant
along a particular streamline; u and v are the velocity components along the
streamline. Their values depend upon the height of the step (h), and the
upstream velocity (Um) (see Figure 7).
An iterative technique (Newton-Raphson) is used to evaluate a
streamline trajectory and the along-streamline velocity. A value for the
atrearafunction (Y) is defined when a bluff shape is specified, and
numerous points along the x-axis are specified between the source position
and the position of the last downwind receptor. For each x, RT and MT are
evaluated through iteration uging Equations 2-41a and b. Then the
streamline height (y) and the velocity (u,v) are computed.
This process works only when Uo, and h are properly specified,
because the streamfunctions used in PFM correspond to the particular choir.es
Uoo = 2 and h = 100. Therefore, user's dimensions are rescaled into
"step space" for computation of the streamline coordinates and velocities
and then scaled back again before the PFM line integrations are performed.
As shown in Figure 7, the scaling process depends on which bluff shape
is chosen. Equations 2-4]a through 2-41g are based on a coordinate system
with its origin at the base of the step. The heights of all streamlines are
measured from this reference point. The operational bluff height in any
application, however, is measured from the upwind height of what is labeled
the surface streamline. The net height change (DH) of this streamline
between a position 100 step heights upwind of the step and a position 2 step
heights downwind of the step is defined as the scaling measure for the bluff
height and plume height ia»the re^. setting. This model bluff height (DH)
decreases as surface streamlines farther from the origin are chosen (smaller
slopes).
In practice, PFM scales all space dimensions by the hill (bluff)
heights, and it scales all velocities by the mean flow speed. This is done
at the start of the program. When the plume path over a particular bluff is
later computed, all length scales are multiplied by DH, and the new velocity
is multiplied by Uoo (=2). This transforms the flow speed, source
location, plume height, and bluff height into the dimensional units of "step
space." The scaled plume height at the source is then added to the height
of the surface streamline at the source position, and the streamline through
this point (the plume streamline) is obtained. Once the series of points
along this streamline and the corresponding velocity components are
calculated, the height of the bluff base (indicated by the dashed line in
Figure 7) is subtracted from the streamline height, and old lengths and
velocities are made nondimensional once again by dividing by DH and Uo,.
This analysis has been based on the two-dimensional problem of flow
normal to an infinitely wide terrain obstacle. Because flows of oblique
incidence are just as likely to occur, a simple geometric approximation has
been incorporated in PFM.
-------
\\\\\\\\\\\^^^
Source Position
Figure 7. Plume path and surface streamlines derived from flow over a step.
-------
As the angle between the wind direction and the direction normal to the
crest line becomes nonzero, trajectory distances to the crest increase. The
trajectory no doubt tends to curve along the face of the bluff somewhat
before sweeping back over the crest. Also, the effective slope of the bluff
face along the trajectory changes. If the streamline curvature is
neglected, then the problem remains two-dimensional and can be handled in
the model by changing the source to obstacle distance and by changing the
effective bluff shape. Each of these changes depends on the cosine of the
wind angle. When the geometry is altered and a new bluff shape is computed,
the plume trajectory analysis along the wind direction proceeds as before.
Flow Over an Ellipsoid
The general case of flow over an ellipsoid of arbitrary shape is
neither two-dimensional nor axisymmetric. Consequently, a unique
streamfunction cannot be defined for each streamline as it was in the bluff
analysis. This means that streamline trajectories must be generated by
selecting a point in space, solving for the velocity at that point and
moving a short distance along ths indicated trajectory to select another
point, and start the process over again.
Outside an ellipsoid moving along its positive x—axis at a speed Ua,
in a fluid at rest, the velocity at any point is derived from the velocity
potential $ (Lamb, 1945).
,m H1 i
$ - K x J ~=-^ (2-42a)
X (a + A')A
A= <(a2 + A'Xb2 + A'Xc2 + A))1/2 (2-42b)
when the equation of the ellipsoid surface is:
222
% * L + i. =1 (2-43)
a b c
The integral is a function of the elliptic coordinate A, defined in
terms of x, y, z, a, b, and c by :
(2-44)
26
-------
Hence, A = constant corresponds to one of a family of confocal ellipsoids,
with X =0 being the surface of the obstacle. The constant K also
an integral over X :
K = a b c IW(2 - AQ) (2-45a)
A = abc I —p - (2-45b)
0 o (az + X'M
Having obtained the velocity potential, the velocity of the fluid
surrounding this translating ellipsoid is:
(u, v, w) = (3§/0x, a#/ayt 3*/az) (2-46)
This relationship changes when the ellipsoid is stationary and when the
fluid is moving along the positive x-axis at a speed Oo>:
(u, v, w) = Ot/ax + Uoo, 3f/3y, 8f /8z) (2-47)
This is the form used in PFM to address air flow over ellipsoidal obstacles.
When the wind is not directed along the x-axis, the incident flow has a
v-component as well as a u— component . In that case, velocities due to each
conponent are calculated separately, then added together. The solution for
flow along the y-axis may be written from symmetry with Equations 2-42 and
2-45:
§ = K y J -~ (2-48a)
X (b2-:- X ' )A
K = a b c Vro/(2 - !„) (2~48b)
The velocity is obtained again from derivatives of the velocity
potential §:
(u, v, w) - Of/ax, 3*/3y * Vco, 3f/32) (2-49)
27
-------
It is important to emphasize that computed plume trajectories over
these arbitrary ellipsoids are approximate because a streamfunction cannot
be specified. The local velocity-increment scheme will either undershoot or
overshoot the actual trajectory each time the tangential velocity is used to
advance to the next point. The accuracy attained is directly related to the
density of sampling points and the curvature of the flow.
2.2.4 Plume Trajectory Analysis - Stratified Flow
Potential flow trajectory solutions presented in Section 2.2.3 are
appropriate for neutral density stratification only. Ambient stratification
can significantly affect plume behavior and resultant atmospheric pollutant
concentrations. For moderate plume heights relative to hill size and for
moderately stable stratification, laboratory experiments suggest that the
plume will go over the top of the hill, but the path of the plume will be
closer to the surface and the flow over the crest will be faster than in
neutral flow. For small plume heights under strong stratification, the
plume will tend to go around the hill rather than over it. If the hill is a
long ridge, the flow may be "blocked" and stagnate upwind. To approximate
first-order effects of moderate stratification, results of EPA Fluid
Modeling Facility (FMF) experiments on plume behavior in a stratified water
channel flow over a polynomial hill (Hunt et al. 1978) have guided the
development of a streamline modification algorithm.
Stratification is assessed in terms of the Froude number (FR) based
upon the hill height (H):
»
FR = (2-50a)
N = I
-------
0.50
0.45
0.40
0.35
0.30
OJ5
0.20
0.15
0.10
0.05
O LARGE TOW TANK
a WIND TUNNEL, Fr = °°
A SMALL TOW TANK, Fr = 1.7
0SMALL TOW TANK. Fr = 1.0
7 SMALL TOW TANK, Fr = »
OF POTENTIAL
IflTHEOHY (NEUTRAL FLOW)
8
Source: Hunt et al. 1978
Figure 8. Height of the source streamline (ns) above the hill crest
(height = a) for various stack heights (Hs) and Froude
numbers (Fr).
29
-------
Call the streamline lift at the crest in neutral flow Ln. Then an
approximate form with the right asymptotic limits (Ln and U«./N) for
the streamline lift in stratified flow (D is:
L - FR H/(l + FR H/Ln) (2-51)
Because the neutral lift is also known away from the hill crest from the
potential flow solutions, Equation 2—51 is assumed to hold everywhere. Note
that as LH becomes vanishingly small far from the hill, L goes to zero as
well.
PFM keeps track of the actual streamline trajectory, not the streamline
lift near the hill. Therefore, Equation 2-51 is transformed from a lift
expression to a depression expression:
D0 = Lj, - L • !„/(! + FRL> (2-52a)
FRL = FR (H/I^) = UV(NLn) <2-52b)
In the neutral limit where Ucc/N » Lj, Oo approaches zero, and
no stratification effects are simulated. In the strongly stratified limit
where U»/N « L,j, Do approaches 1^ - u^/H, the difference
between the neutral streamline deflection, and the length scale for vertical
deflections in stratified flow. When Ln is zero, so is Do,
Application of Do to potential flow streamlines tends to "relax"
vertical streamline deflections around a hill with increasing stratifica-
tion, without regard for the consequences of allowing streamlines to pass
through the hill. Although some of the flow should not pass over the crest
for Froude numbers less the one, comparison of Equation 2-52 with Figure 8
shows that the computed relaxation is much too severe. Streamlines are
straightened out too rapidly with increasing stratificaton. The expression
for Do must include some adjustment for the presence of the hill.
An adjustment region of order Uoo/N close to the hill surface is
missing in Equation 2-52. This is the length scale in which significant
vertical motion is allowed. All streamlines that rise over the hill must be
compressed into this region. This adjustment region is incorporated in the
depression equation (2-52a) by adding an exponential adjustment factor. The
factor is ad hoc and merely requires the computed depression to tend toward
zero when the neutral streamline passes very close to the obstacle (compared
to the length scale Uto/N). If the neutral streamline is several scale
heights above the local surface, the correction factor tends to 1 and the
depression is equal to Do. The corrected streamline depression is called
D;
D - D0 (1 - exp(-S,/FR H)) (2-53)
where I is the local neutral streamline height above the obstacle.
30
-------
When Che streamline height is adjusted using Equation (2-53), the along
streamline velocity must be recomputed. Conservation of mass indicates that
the ratio of the stratified flow velocity (us) to that in neutral
flow (un) depends upon the vertical rate of change of the streamline
depression:
~) <2-54)
un dz
Furthermore, a change in the trajectory requires a redistribution of the
along-streamline wind speed among the u, v, w, components so that the
velocity still points along the streamline. PFM adjusts the streamline wind
speed according to Equation 2-54 and then computes the individual components
based upon the slopes of the depressed streamline trajectory.
2.3 The COMPLEX/PFM Option
The COMPLEX models (COMPLEX I and COMPLEX II) are multiple point source
sequential terrain models formulated by the Complex Terrain Team at the EPA
Workshop on Air Quality Models held in Chicago in February, 1980. COMPLEX I
is a univariate Gaussian horizontal sector-averaging model (sector width ~
22,5°), while COMPLEX II computes off-plume-centerline concentrations
according to a bivariate Gaussian distribution function. Both models are
very closely related to the MPTER model in both structure and operation.
Anyone who is not familiar with either COMPLEX or MPTER should consult the
MPTER user's manual (Pierce, et al. 1980).
Terrain treatment in COMPLEX varies with stability class. Neutral and
unstable classes use a 0.5 terrain adjustment, while stable classes use no
terrain adjustment when the recommended options are selected. With 22.5°
sector averaging, COMPLEX I performs sequential VALLEY plume impingement
calculations for stable cases. COMPLEX II plume impingement calculations
are similar, with the exception that sector averaging is not used.
COMPLEX/PFM has the ability to utilize PFM calculations for neutral to
moderately stable flows. The PFM option invokes either COMPLEX I, COMPLEX
II, or PFM computations depending upon the stability class and the Froude
number. Unlike previous versions, however, all sources must be located at
the same point (as in CRSTER).
COMPLEX II is invoked whenever the stability class is either 1, 2, or 3
(A, B, or C), regardless of the Froade number. In these cases plume growth
is rapid and the details of terrain adjustment are not so important. A 0.5
terrain adjustment is an adequate representation of average plume behavior.
31
-------
PFM is invoked for stability classes 4, 5, and 6 (D, E, and F) whenever
the flow along the plume streamline has enough kinetic energy to rise
against the stable density gradient and surmount the highest terrain
elevation along the wind direction. Such a plume is said to be above the
critical dividing streamline of the flow (Section 2.3.3).
COMPLEX I is invoked whenever the plume is found to be beneath the
critical dividing streamline of the flow. Flumes beneath the dividing
streamline no longer pass over the terrain peak and therefore nay impinge on
the face of the the hill somewhere. Thus, the PFM option defaults to
VALLEY-like computations for impingement cases. This can potentially occur
in conjunction with stability classes 4, 5, and 6; but, class 4 occurrences
may be rare.
The PFM option enhances the ability of COMPLEX to perform complex
terrain Gaussian plume dispersion computations in two important areas.
Firstly, it incorporates plume deflections and distortions derived from
potential flow theory. This enhancement approximates at least first-order
terrain effects on plume geometry. And, because the streamline computations
vary with obstacle shape, plume height and Froude number, plume distortions
are coupled directly to meteorological variations and the approximate
terrain geometry in a way that no single terrain adjustment could be.
Secondly, the use of the PFM option requires vertical temperature and
velocity information to characterize the Froude number, the critical
dividing streamline, and stable plume rise. Availability of the Froude
number and the dividing streamline removes the assumption of coupling
between the surface dispersion stability class and the dynamics of the flow
aloft at plume elevation under stable conditions. It is not necessary to
identify plume impingement with class E or F dispersion conditions.
Elements of the PFM option are presented in the following sections.
Section 2.3.1 explains how PFM computations are made in both a long- and
short-term version of COMPLEX/PFM; Section 2.3.2 summarizes the plume rise
computations; Section 2.3.3 addresses the computation of dividing
streamlines and Froude numbers and their use in the model; Section 2.3.4
describes incidental changes to COMPLEX required for the PFM option.
2.3.1 PFM Option: Long and Short
Two PFM option choices are provided in COMPLEX/PFM. The short-term
option calls PFM as a subroutine and is designed for critical period
analyses in which 3-hour or 24-hour concentrations near one particular
terrain feature are to be evaluated. The long-term option searches a
look—up table of intermediate PFM results and is designed to be run for a
full year of sequential hourly meteorological data, producing concentrations
and concentration statistics for up to 180 receptors. It may have as many
terrain features as there are receptors. The long-term version of the PFM
option should be used as a sequential analysis model for calculating annual
average concentrations, for siting air quality monitors, and for identifying
32
-------
those periods within the meteorological data which produce the highest
predicted concentrations.
PFM-Long
PFM-Long is constructed around a series of PFM calculations for 5 plume
heights, 6 Froude numbers, and 20 discrete obstacles shapes. It is not
designed to reproduce actual PFM results for the specific plume heights,
meteorological conditions, and obstacle shapes of a particular application.
Instead, it is designed to provide credible, conservative approximations.
All of the PFM computations in the look-up table assume that the plume
travels directly over the crest of the obstacle, perpendicular to the crest
line. The user must supply the model with those obstacle shapes that best
represent each of the terrain features. The user is given 4 bluff choices
and 16 ellipsoid choices. These are summarized in Section 2.5.
One controlling obstacle must be identified for each receptor.
Usually, there will be one obstacle specified along each receptor ray and
all receptors along that radial will have the same controlling obstacle.
However, if it is not clear which of a series of terrain features along a
radial is most important, several obstacles may be specified. In this case,
receptors more than one hill height downwind of the nearest obstacle should
be identified with the next closest obstacle downwind of the receptor.
Note that PFM (or COMPLEX) does not include wake physics. If a plume
is strongly deformed by the leading obstacle along a receptor ray, then
model predictions downwind of this obstacle are invalid. Predicted
concentrations downwind of the first terrain feature should be considered
adequate only if upwind terrain features are no larger than one half of the
expected plume heights.
The look-up table of PFM calculations contains 8 variables at 49
non-dimensional downwind distances for each possible plume height, Froude
number, and obstacle shape combination (440 in all). Of these eight
variables, six are presently used to calculate PFM factors. These six are:
PTA, the elapsed time from the source along the plume trajectory; ANS, the
distance from the plume centerline to the surface; US, the wind speed along
the plume centerline; R2S, the square of the local effective radius of
curvature of the obstacle; PHI, the line integral of UR^ along the plume
trajectory; and TEE, the line integral of 1/{UR2) along the plume
trajectory.
The source position end the position of the last PFM receptor in the
tabulated computations are selected so that the source is located upwind of
the region of substantial influence around the obstacle and the final
receptor is about one obstacle length downwind of the obstacle. This allows
PFM factors (Equation 2-40) to be calculated for arbitrary source
positions. If the source is loceted closer to the obstacle than the source
33
-------
in the PFH computations, the Case 1 version of Equation 2-40 is used. If
the source is located farther from the obstacle, the Case 2 version is used
(see Figure 9).
Case 1 assumes that the source is located at the one downwind PFM
receptor point (S) c" sest to the actual source position and all actual
receptors are located at the nearest PFM receptor point downwind of the
source. Then, approximate PFM factors for each of these receptors (R) are
computed using the following equations:
CFAC(R)
SZFAC(R) -
f
f
- PHI(S))(TEE(R) - TEE(S}}
(US(S)(PTA(R) - PTA(S)))2
PHI(R) - PHI(S)
(— i— .)
R2S(R)(PTA(R) - PTA(S)) US(R)
(2-55a)
(2-55b)
SYFAC(R)
HFAC(R)
f
R2S(R)(TEE(R) - TEE(S))
PTA(R) - PTAfS)
ANS(R) / ANS(S)
(2-55c)
(2-55d)
Here PHI, TEE, and PTA are the line integrals <{>(s), T(s), and t(s);
therefore, only the change in these values from S to R are needed. All
other quantities are evaluated at the receptor point R. US(s) is u (s),
R2S(s) is R2(s), and ANS(s) is the distance from the plume centerline to
the surface (see Equations 2-32, 2-33, and 2-40).
Case 2 assumes the nondimensional distance from the actual source to
the FFM source is d. Because plume deformation is negligible over this
distance, the line integrals along the plume streamline are constant and all
PFM factors are unity. Consequently, PFM factors only need to be calculated
for receptors that fall within the 49 array entries. Again, the PFM
receptor points closest to actual receptors «re selected, and PFM factors
are calculated with the following equations:
CFAC(R)
(R2S(l)US(l)d + PHI(R))(d/(R2S(l)US(l)) + TEE(R))
+ PTA(R)))2
(2-56a)
SZFAC(R)
SYFAC(R)
HFAC(R)
f
R2S(1) US(1) d + PHlQU , 1
R2S(R)(d/US(l) + PTA(R)) ^US(R)
)
R2S(R)(d/(R2S(l)US(l))
TEE(R))
d/usd) H
ANS(R) / ANS(l)
PTA(R)
(2~56b>
(2-56c)
(2-56d>
34
-------
Casa I: Actual source is located downwind of the PFM source
O PFM Source
• PFM Receptor
Actual Source Location
Actual Receptor Location
Case II: Actual source is located upwind of the PFM source
O PFM Source
• PFM Receptor
1
Actual Source Location
Actual Receptor Location
s
8
Figure 9. Relationship between PFM and actual sources and
receptors in PFM-Long,
Note: PFM receptors labeled R,S, or 1 are used in Equation 2-55 and
Equation 2-56-
35
-------
Again, the line integrals PHI, TEE, and PTA are adjusted so that they
encompass the entire interval between the source and the receptor. For
example, R2S(l)US(l)d is the "PHI" integral from the actual source to the
PFM source and PHI (R) is its value from the PFM source to the receptor
evaluation point (R).
PFM-Long starts out by setting up factor arrays for each receptor.
Because each receptor has a fixed relationship to its own controlling
obstacle and the location of the sources, factors are calculated for all
plume height-Froude number combinations (22). This is done once at the
start of the program.
Receptor-specific PFM factors modify the COMPLEX model's concentration
equations in much the same way as they modified Equations 2-34 and 2-35.
Within a particular hour loop, one of the six Froude number classes is
selected on the basis of the hourly Froude number, and the PFM factors at
each receptor are interpolated for each source. Each source can have a
different plume height, so the PFM factors are interpolated in height only.
These factors do not affect computed downwind and^ off-axia distances within
the host model.
PFM-Short
PFM-Short contains the PFM code as a subroutine. Instead of
approximating the source-terrain geometry, plume height, and meteorology as
in PFM-Long, PFM-Short performs the appropriate PFM calculation. This means
that more obstacle shapes are at the user's command, the obstacle
orientation to the wind trajectory is preserved, and actual plume heights
and Froude numbers are incorporated in the computation.
The user may specify only one terrain obstacle. This is consistent
with the refined purpose of the short-term model. Additional obstacle
descriptors required for PFM-Short are described in Section 2.5.
PFM factors are computed at 49 points along the wind trajectory for
each source. COMPLEX determines the downwind distance and crosswind
distance to each model receptor and selects the PFM factors from the FFM
receptor point closest to the downwind distance. The geometry is summarized
in Figure 10. When the particular factors are selected, concentrations are
computed as in PFM-Long.
2.3.2 Plume Rise
Wind and temperature information available from upper dir soundings is
used in the plume rise calculation when the PFM option is selected. A
stepwise integration is performed through each vertical layer, allowing
different vertical potential temperature gradients and wind speeds in each
layer. Thus, the effects of elevated inversion and vertical wind shear can
be accounted for in the plume rise calculations.
36
-------
'stance
10.
37
-------
The technique described by Brigga (1975) Is used to integrate the
stable plume rise equation. The reduction of buoyancy flux is calculated
layer-by-layer, until all of the plume's buoyancy is consumed. The height
at which the buoyancy flux goes to zero is the equilibrium or final plume
rise height.
For a vertical plume,
Fz = Fj - 0.027 Sj Fj/3 (Z8/3 - Zj8/3) (2-57)
where Fz is the buoyancy flux at height 1 (height above stack height);
Fj is the Iruoyancy flux at the bottom of the jth layer;
Fo is the initial buoyancy flux of the plume;
Sj is a stability parameter [(g/Ta)80/3Z] in thF jth layer.
For a bent-over plume,
F - F. - (1/3)S'2 S u. (Z3- Z,3) (2-58)
•J J J . J
where B" is an entrainment parameter, and
uj is the average wind speed in layer j .
The formulation (either bent-over or vertical plume) that gives the largest
decrease in buoyancy flux in a given layer is used in that layer.
The entraikiment parameter, $', is assumed equal to 0.41314. Thus, with a
single layer of constant potential temperature gradient and wind speed, Equation
2-58 reduces to standard Briggs final stable plume rise equation.
Z - 2.6 (F0/u Si)l3 (2-59)
To calculate neutral plume rise, the Briggs (1969) neutral plume rise
equation is applied in a stepwise fashion with a constant wind speed assumed
in each layer. Plume rise is calculated in each layer by means of the
average wind speed in that layer. A virtual source technique is usea to
match the plume rise at the interface of layers.
38
-------
Figure 11 illustrates the procedure for calculating neutral plume
rise. The virtual distance (Xv); at which the plume enters layer j is:
Hj UJ
(X). = ( ~
v 3 1.6 F 1/3
3/2
(2-60)
The actual downwind distance, Xa, at which the plume enters layer j
j = 2
Xj_ + I [(Xv)i - (Xv)i_il j >_ 3
i=3
(2-61)
The layer, J, in which final plume height is reached is determined when
Xa _> Xf, where
3.5 X*
(2-62)
X*
14F5/8
34F2/5
55
F > 55
(2-63)
Thus, the virtual final rise distance, (Xv)f,
(Xv)f = (Xv)j + (Xf -
(2-64)
The final neutral plume rise, Zf, is:
1.6
(Xv)f
2/3
UJ
(2-65)
39
-------
Layer
Number
(Xv)4-
Figure 11. Calculation of neutral plume rise.
-------
If the wind speed does not vary with height, (Xv)f equals Xf
and Equation (2-65) reduces to single layer Briggs final plume rise equation.
The stable plume rise equation allows unrealistically large plume rise
when the vertical potential temperature gradient obtained from the sounding
data is negative or near zero. Likewise, the layered neutral plume rise
equations may overpredict plume rise under conditions when an elevated
inversion exists above a neutral layer. For these reasons, the plume rise
used by the model when the PFM option is requested is always the minimum of
the layered stable and neutral plume rise equations, regardless of stability
classification at the surface.
If sounding data are missing for a particular hour (or if the PFM
option is not requested), the single layer Briggs stable and neutral plume
rise equations are used to determine plume height.
2.3.3 Froude Number and HC Analysis
The Froude number characterizes the balance of inertial and buoyancy
forces over the scale of the obstacle height. If it is less than one,
streamlines close to the surface upwind of the hill are more likely to go
around the obstacle than over it. PFM computations do not address these
streamlines, so it is important to define the height of the critical
dividing streamline (HC), and apply PFM calculations only to the flow above
this level.
The height HC is found from a simple energy argument suggested by
Snyder, et al. (1981). The kinetic energy of the dividing streamline upwind
of the obstacle must just balance the energy required to lift the streamline
from its initial height (HC), to the top of the obstacle (H) through the
density gradient:
H
0.5 p u2 = g I (H-Z)(-dp/dz)dz (2-66)
HC
If the potential energy gained in lifting a streamline exceeds the
initial kinetic energy, then that streamline must lie below HC. Conversely,
if the potential energy is less theu the initial kinetic energy, the
streamline must lie above HC. When the kinetic energy of the surface
streamline exceeds the accrued potential energy, there is no critical
dividing streamline (HC=0).
COMPLEX/PFM evaluates HC through Equation 2-66 by integrating over the
temperature soundings. The temperature gradient is constant between
observation levels, and the wind speed is obtained by linear interpolation
between observation levels.
41
-------
A bulk Froude number for the layer above HC is calculated from:
HT
FR =
N(H-HC)
(2-67)
HC
where HT is a height above the obstacle equal to the maximum of either 1-5
the obstacle height, or the sum of plume height and obstacle height. The
wind speed is the mean for the layer between HG and HT, and N is based upon
the temperature difference across the layer.
Both the HC and FR computations are treated as local analyses for the
plume streamline. Each must be based upon terrain heights along the plume
trajectory. In PFM-long, each receptor radial ie assumed to pass over the
crest of the controlling obstacles. PFM-short incorporates details of the
plume trajectory over a single three-dimensional obstacle; the maximum
terrain height along the wind direction in this case is computed from the
equation of the ellipsoidal hill surface.
After the dividing streamline height and Froude number are computed,
effective plume and obstacle heights are defined. The layer beneath HC has
been assumed to either stagnate, or flow around rather than over the
obstacle. Plumes in this layer are treated in the COMPLEX I branch of the
model. Plumes above HC still pass over the obstacle, but now they pass
closer to the surface since the lower layer no longer competes for space
above the obstacle. Consequently, the depth of the lower layer is
subtracted from both the plume height and the obstacle height before PFM
factors are calculated.
The resulting factors for plume deformation and centerline
concentration can be used directly as before, but a correction must be made
for the plume centerline height (Figure 12). The net plume height above the
surface is not solely determined by the effective plume height times HFAC.
Away from the obstacle, the net height should be equal to the full plume
height before the HC correction. Closer to the obstacle, the surface of the
terrain rises toward HC, and the net height decreases.
This process is accounted for by setting the net plume height (HH)
equal to the sum of the effective plume height (HE) times HFAC and the
difference between HC and the local terrain height (ZH, less than or equal
to HC):
HH = (HEHHFAC) + (HC - ZH) (2-68)
42
-------
HH(B) = HE*HFAC
HH(A)=HE*HFAC + HC-
Plume Path
Figure T2. Relationship between PFM geometry and plume height (HH)
when HC, the dividing streamline height, is not zero.
-------
This is done only when 2H is less than HC, When ZH exceeds HC, the
second term in Equation 2-68 is set equal to zero.
2.4 Temperature and Velocity Profile Interpolation
Adequate meteorological profile information is essential to modeling
elevated plumes in complex terrain. The behavior of streamlines at plume
height depends on wind speeds and temperature gradients from the surface to
the top of the important terrain features and beyond. If such profiles of
temperature and wind speed are not available, then the utility of
fabricating Proude numbers and dividing streamline heights from standard
surface observations alone is questionable. COMPLEX/PFM predictions using
such Froude numbers offer no advantage over predictions of either COMPLEX I
or COMPLEX II.
2.4.1 Required Data and Defaults
The most desirable meteorological data include hourly wind velocity and
temperature soundings from the modeling site. The next most desirable
on-site data include at least twice-daily temperature and velocity soundings
(early morning and evening) and an on-site instrumented tower at least as
tall as the source stack. The least desirable (but still modelable) data
set includes representative off-site upper-air data and surface data. This
last choice is essentially equivalent to using standard available airport
data.
Only the first of the above possible data sets contains hourly profiles
for calculating hourly values of plume rise, FR, and HC. Consequently, a
profile preprocessor has been constructed to interpolate between twice-daily
soundings and thereby create an hourly data set. It is described in the
next section.
A full year of hourly modeling data requires a complete set of
twice-daily profiles, including an evening sounding from the last day of the
preceeding year and a morning sounding from the first day of the following
year. If the complete set is not available, then gaps will exist in the
interpolation set. These gaps will be flagged in the program, and default
calculations will be performed in COMPLEX/PFM using the host COMPLEX I or
COMPLEX II algorithms, depending on the surface stability class.
2.4.2 Mixing Heights and Profile Interpolation
The sounding preprocessor is designed to be used with the data obtained
from the National Climate Center (NCC). Sounding data from the national
network of upper-air stations are contained on "TDF 5600" tapes.
The first part of the preprocessor, READ56, reads pressure,
temperature, wind speed, and wind direction data from the sounding tape.
These are formatted for further analysis and written to a file for editing.
Only pressure levels with good temperatures are placed into the file, and
-------
missing winds are flagged. The user must elect to either remove pressure
levels with bad winds, or replace the missing data with "reasonable
values". The user must also specify the uppermost standard pressure level
needed for profile interpolations. The default is the 700 mb level.
Temperature data must be available at this level on all soundings. Any
missing entries are flagged by the preprocessor. When data other than an
NCC TDF 5600 tape is used, the user must prepare and format the data to look
like the output from READ56.
The second part of the preprocessor, PROFIL, computes hourly mixing
heights and interpolates the temperature and wind data to construct hourly
sounding profiles between the actual sounding times. Mixing height
computations are based in part upon the Benkley-Schulraan scheme (1979).
Temperature and wind velocity interpolation methods are consistent with the
hourly mixing height determinations.
Temperature Interpolation ~ Morning to Evening
Interpolated soundings between the early morning (1200 GUI} and the
evening (0000 GMT) observations must account for the growth and partial
decay of the daytime convective layer when this layer depth exceeds the
mechanically mixed layer. The interpolation begins with a calculation of
the maximum height of the convectively mixed layer (HLMAX) for the day.
First, a temperature profile for the time of the maximum relative surface
temperature is constructed by interpolating all vertical profile levels in
time. Then a dry adiabat is drawn from the maximum surface temperature to
the constructed temperature profile. The point of intersection defines
HLMAX. This method differs from the Benkley-Schulman method in that each
data point in the profiles determines the local temperature shift in time
due to advection, rather than just the 700 mb level.
Once HLMAX is determined, the profile interpolation above and below
HLMAX is handled differently. All profile points above HLMAX are unaffected
by surface heating; therefore, hourly temperatures are found by simple
linear interpolation in time at all levels greater than or equal to HLMAX
(Figure 13). Below HLMAX, the surface temperature progression must be
incorporated. The hourly adjustment is as follows:
o adjust profile below HLMAX by shifting temperature an amount equal
to the shift at HLMAX;
o follow the adiabat from the surface temperature to the point of
intersection with the shifted profile;
o compare the height of the intersection with the mechanical mixed
layer height (see Benkley et al. 1979);
o create a new profile point at the maximum of these two heights and
neglect all profile points below it; and
o call this height the mixing height for this hour (HL).
-------
AM
PM
> Linear Interpolation in Time Along
i Constant Pressure Surfaces
HLMAX
Figure 13. Sample construction of temperature profile for fourth hour
after A.M. sounding.
46
-------
If HL is less than the previous value of the mixing height, and if the
surface temperature has not reached its maximum value for the day, HL is set
to its preceeding value, and the temperature at HL is incremented an amount
equal to the change at HLMAX.
Once the maximum surface temperature has been reached, the profile
below HLMAX will slowly relax, approaching the shape of the evening
sounding. The mixing height persists at HLMAX until the night-time
mechanical shear layer becomes dominant, and this transition is determined
as in the Benkley-Schulraan method. Consequently, temperatures below HLMAX
are linearly interpolated from the dry adiabat to the actual evening
sounding.
Throughout this discussion, reference has been made to a "surface"
temperature. If on-site meteorological tower data are available, the
"surface" temperatures referred to above should actually be the uppermost
measurement height. In this way, site-specific surface layer changes can be
included directly. The user may then merge the lower tower measurements
into the profile data set with a suitable data manipulation routine.
Temperature _Inter]^o_l_ajti_qn •^Evening to Morning
The remnants of the daytime mixed layer continue to decay in the
evening. Cooling proceeds more rapidly at the surface than at the middle
levels of the profile. Therefore, linear interpolation of temperatures at
all sounding levels is appropriate. Regions of greatest temperature change
will have the greatest hourly rate of change, so the lower profile will
stabilize most rapidly.
Velocity Interpolation
Wind velocity interpolation is likely to be far less accurate than
temperature interpolation. For this reason, the results of the wind speed
interpolation are used for FR-HC analysis and plume rise calculations only,
and the wind directions are not used at all. Linear interpolation is used
throughout after profile height levels are computed for the temperatures.
The interpolation proceeds as follows:
o interpolate all speeds and directions between soundings;
o interpolate in height to compute winds at new pressure levels
created in the temperature profile interpolation (i.e., HLMAX); and
o compute wind speed and direction at the hourly HL heights using
the surface stability class, the recommended power law exponents,
and the "surface" wind speed and direction from the surface
meteorology data file.
-------
This interpolation method separatee the surface layer flow from the upper
level flow. Winds measured on-aite or at the nearest NWS station control
the wind speed and direction throughout the mixed layer. Radiosonde winds
observed only twice a day are not allowed to influence low-level winds at
all.
2.5 Terrain Description Guidance
Two classes of obstacle shapes are incorporated in PFM: bluffs and
ellipsoids. Bluffs are strictly two-dimensional features with one ascending
face, but no descending face* Ellipsoids are truly isolated hills, and may
be either two- or three-dimensional. The t^n-dimensional ellipsoid is
applicable to long ridges, while the three-dimensional ellipsoid is
applicable to mounds, peaks, or abort ridge elements.
Twenty bluff profiles and an infinite number of ellipsoid shapes are
available within PFM-Short. A subset of four bluff profiles and sixteen
ellipsoid shapes is contained in PFM-Long.
2.5.1 PMF-Short Obstacle Selection
Profiles of twenty bluff sections are presented in Appendix A. Each
section has a characteristic mean slope, and a characteristic crest location.
To pick the bluff shape most like an actual terrain feature, construct
a terrain section through a representative part of the feature. Make sure
that the section is along a cut perpendicular to the crest line. Scale the
length units by the crest height in such a way that the profiles are the
same size as the bluff profiles contained in Appendix A.
Now match the general shape of the terrain profile to one of the twenty
bluff shapes. The most important features to match are terrain slopes near
the top, and the overall scale of the terrain shape. The number of the
bluff shape most like the terrain section should be specified as parameter
"BLFSH." The actual bluff height should be entered as "HAO."
The source location (XO, YO) relative to the center of the obstacle,
and the orientation (OBSANG) of the obstacle must also be specified. Each
bluff shape has a zero point identified on the horizontal axis. The actual
perpendicular distance of your source from the corresponding point on your
terrain section is XO. YO if zero. Note that the value entered for XO must
be negative, since the positive X-direction points from the source to the
obstacle, and the source is positioned on the negative side of the origin
(unless your source lies between the zero point and the crest).
The direction of the positive X-axis also defines the bluff
orientation. OBSANG is the angle between North and the vector direction X.
This angle is measured in a clockwise direction (Figure 14),
-------
North
Wind
Figure 14. Relationship between north and the orientation of the
bluff coordinate system.
49
-------
Ellipsoid shapes are based on along-wind (AWARD) and cross-wind (CWARO)
aspect ratios (Figure 15). An aspect ratio is defined as one-half the
obstacle length (or breadth) divided by the obstacle height. The along-wind
aspect ratio is calculated using the major ellipse axis that penetrates the
windward face. The cross-wind aspect ratio is calculated using the major
ellipse axis perpendicular to this axis. The "windward" face assumes that a
wind is blowing from the source toward the center of the obstacle. When the
windward face is correctly chosen, the source should lie within a 45° cone
about the -x axis.
As with the bluff, the user should try to pick aspect ratios that
produce an ellipsoid most like the shape of the actual terrain feature. A
good starting point is to calculate aspect ratios based on the shape of the
terrain feature along its one-half-height contour. Sketches of the
resulting ellipsoids should then be compared to terrain cross-sections*
After selecting the ellipse axes, the obstacle origin is defined
(center of ellipse), and the source position (XO, YO) and obstacle
orientation can be determined. The positive X-axis points away from the
source along the "along-wind" ellipsoid axis. Therefore, XO and YO are
readily obtained and OBSANG is once again the angle measured clockwise from
North to the vector pointing along the positive X-axis (Figure 16). The
terrain height (KAO) is the maximum height of the terrain feature that is
consistent with the specification of the ellipsoid.
2.5.2 PFM-Long Obstacle Selection
Specification of obstacle shapes and position parameters for PFM-Long
is a simplified version of the PFM-Short process. The user may select a
bluff shape from four possible profiles, or an ellipse shape from sixteen
possible combinations of along-wind and cross-wind aspect ratios. Once this
is done, only the source-to-obstacle distance and the maximum obstacle
height need be obtained since obstacle orientation to the wind is neglected.
All available obstacle shapes are numbered one through twenty. The
first four are bluff shapes, the rest are ellipsoids. Table 1 summarizes
which bluff shapes and ellipsoids are included. The obstacle shape is
specified as parameter IOBSH. The distance between the source and the
obstacle center is specified as parameter ODIST. This quantity should not
be negative. It is a true d5stance and not a source coordinate relative to
the obstacle center.
Selection of a controlling obstacle for each receptor depends on the
complexity of the terrain. FFM computations aia only valid for plumes
upwind of and over obstacles. No wake effects in the lee of obstacles are
included. Therefore, only the dominant terrain features surrounding the
source should be considered.
50
-------
Cross-wind aspect ratio
CWARO = B/C
Wind
Along-wind aspect ratio
AWARD =A/C
Figure 15. Definition sketch for selecting cross-wind and along-wind
aspect ratios.
51
-------
North
(XO.YO) Source
Figure 16. Relationship between XO, YO, and OBSANG.
-------
TABLE 1. OBSTACLE SHAPES AVAILABLE
FOR APPLICATION,1- OF PFM-LONG
IOBEH
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
BLFSH
1
8
12
16
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
CWARO
0
0
0
0
1
2
5
10
1
2
5
10
1
2
5
10
1
2
5
10
AWARO
0
0
0
0
1
1
1
1
2
2
2
2
5
5
5
5
10
10
10
10
53
-------
However, when major plumes rise above nearby obstacles, the most
important terrain features may lie farther away. In cases where the
dominant features are not readily apparent, more than one obstacle may be
specified along a single receptor radial. In this case, receptors more than
one obstacle length downwind of one obstacle should be associated with the
next obstacle along the radial. This method effectively treats succeeding
features in isolation from one another.
54
-------
SECTION 3
USER'S INSTRUCTIONS
The COMPLEX/PFM modeling package is schematically illustrated in
Figure 17. User's instructions for the SETUP file restructuring program,
the READ56 and PROFIL preprocessing programs, and the COMPLEX/PFM model
program are presented below.
3.1 SETUP Instructions
The SETUP program rearranges the data contained in the formatted PFM
look-up table supplied with COMPLEX/PFM. No program control variables are
required to run SETUP because the structure of both the input and the output
file is predetermined. System channel unit 8 is assigned to the formatted
input file, and unit 20 is assigned to the random access output file. This
program should be run only when the COMPLEX/PFM modeling system is first
installed on the user's system.
3.2 BEAD56 Instructions
The READS6 program reads an NCC TDF 5600 data tape of upper air data,
selects those sounding levels within a pressure window specified by the
user, checks each sounding time for missing or multiple entries, and flags
entries that contain missing data. Pressures, temperatures, calculated
heights, and wind velocities are written to a formatted file for editingj
and for later use in the program PROFIL (see Section 3.3).
One program control card is needed to reset nameliet variables when
other than default values are desired. The namelist variables are described
in Table 2.
Two system channel units are needed for data file input and output.
The file containing TDF 5600 data is assigned unit 10, and the output file is
assigned unit 7.
The first record in the output file created by READ56 repeats the
namelist variables in PARAMS. This record is followed by one data
identification record and several data records for each sounding. Table 3
describes the contents of the data identification record. Data from each
sounding level follows, with four strings of pressure (mb), height (m),
temperature (°K), wind direction (degrees), and wind speed (ra/s) per
record. If missing data were found in the sounding, or if a sounding is
55
-------
Figure 17. The COMPLEX/PFM modeling system.
56
-------
TABLE 2. PARAMS NAMELIST - READ 56
Variable
Type
Deacription
Default
START Integer Year and Julian day number
of the first complete day of
profiles (e.g., 81001)
STOP Integer Year and Julian day number of —
the last comp'.ete day of
profiles (e.g., 81365)
SPTOOZ Real Starting pressure level for 2000.
data to be extracted from the
OOZ sounding (mb)
EPTOOZ Real Ending pressure level for data 700.
to be extracted from the OOZ
sounding (mb)
SPT12 Real Starting pressure level for data . 2000.
to be extracted from the 12Z
sounding (mb)
EPT12 Real Ending pressure level for data 700.
to be extracted from the 122
sounding (mb)
57
-------
TABLE 3. DATA IDENTIFICATION RECORD - BEAD 56
Variable Type
DeBcription
Default
ITPDK
NOSTA
LVL
NLVL
3X.I4
5X.I5
IOBTM(4) 5X.4I2
5X.I2
80X,I3
Label identifying the tape
as aeries "5600"
Identification number for the
upper air station
Year4 month, day, and hour
of the sounding (GMT)
Total number of sounding
levels in the profile
Number of sounding levels
extracted from the profile
58
-------
missing, or entered more than once, additional records containing warning
information are written before the next sounding is obtained.
3.3 PROFIL Instructions
The PROFIL program requires one input card to specify program control
variables, one input disk file containing the processed surface
i meteorological data, and one input disk file containing edited twice daily
temperature and wind velocity soundings. It produces a list file, and a
file of hourly wind and temperature profiles. The default system unit
number for the sounding file is 7; the default system unit number for the
surface data is 2.
Namelist data entry format is used on the program control card. The
control variables that make up namelist "INPUT" and their default values are
described in Table A. The user only needs to change variables whose default
values are unacceptable.
The surface meteorology file is unformatted and is normally obtained as
output from RAMMET or the CRSTER preprocessor. Its content is described in
Section 3.4.2 (Table 21).
The file of twice-daily soundings (OOZ and 12Z) is a formatted data
set. The first record is a namelist statement (PARAMS), - h sounding that
follows contains one parameter record and several recorab - jounding data.
Each record type is described in Sectic.i 3.2. -
Before the soundings file can be used by PROFIL, the user must add the
soundings for the day after the last day to be processed. This is necessary
because the 12Z profile on the final day really corresponds to a morning
sounding. Both the evening sounding plus the morning sounding of the
following day are needed to complete the mixing heights and profile
interpolations through hour 23 of the final day.
The user must also edit the soundings file to identify missing data
levels or missing soundings. If the wind speed is missing at one level, for
example, the user might either interpolate a value based on the values at
adjacent levels or remove the sounding level entirely. If the latter action
ia taken, the number of good levels in the sounding (NLEV) must be reduced
accordingly. If an entire sounding is missing, then only the sounding
information record should be present, with NLEV set to zero.
PROFILE will interpret this zero as a gap in the sounding data, and no
hourly profile data will be calculated. Instead, the output file for this
day will contain 24 data entries reporting a mixing height of 9999, and a
maximum number of profile points of zero. CCMPLEX/PFM will then ignore the
sounding data for this day, and perform either COMPLEX I or COMPLEX II
computations, depending upon the stability class.
59 x*
-------
TABLE 4. PROFIL PROGRAM CONTROL CARD VARIABLES
Variable
NDAYS
JTIME
Type
Integer
Integer
Description
Number of days to be run
Local time of the 002
Default
19
sounding (1900 EST)
ZMEAS Real Anemometer height above the
surface (meters)
ZQ Real Roughness length scale
representative of the area
(meters)
LMETOT Logical Control for writing program
output to disk
IMET Integer System unit number for the
program output data file
ISURF Integer System unit number for the
surface meteorology file
IUP Integer System unit number for the
upper air meteorology file
ICARD Integer -System unit number for card input
ILIST Integer System unit number for the program
output list file
PEXP Real (7) Wind profile power law exponents,
one for each stability
class (1-7)
10.0
0.05
.TRUE.
5
6
0.10
0.15
0.20
0.25
0.30
0.30
0.30
60
-------
't-3.4 COMPLEX/PFM Instructions
Most program control options available in COMPLEX/PFM are unchanged
from MPTER. These will not be restated here to the same depth. The user
should read and understand the MPTER user's guide prior to running
COMPLEX/PFM.
Portions of the model which are new, or portions which have been
changed, are discussed in the following sections. Section 3.4.1 provides a
list of all options available in the model, and describes those which are
new and those old ones which are not compatible with the PFM options.
Section 3.4.2 contains a description of input card and tape format.
3.4.1 COMPLEX/PFM Options
All COMPLEX/PFM options are listed in Table 5. Two of these,
options 25 and 26, are not described in the MPTER guide.
Option 25 was added when COMPUEX I and COMPLEX II were designed from
the MPTER code. This option controls which of several terrain adjustment
algorithms is to be used. In brief, these are:
» OPT 25 = 0 Terrain adjustment follows MPTER.
• OPT 25 = 1 Partial height correction made, but plume may not come
closer to the surface than ZMIN (meters).
• OPT 25 = 2 One calculation is made with the receptor at plume
height over level terrain. A second calculation is
made as in OPT 25 = 1, and the lesser of these two
calculations is used.
• OPT 25 = 3 Concentrations are calculated as if there is no
terrain, except that the receptor is placed at the
actual terrain elevation above sea level.
• OPT 25 = 4 Concentrations are calculated as in OPT 25 = 1 unless
the plume path correction is zero and the terrain
exceeds the plume height. In that case, the solid
ground is placed a distance ZMIN (meters) below the
plume, and the receptor is placed at the actual
terrain height above the plume.
• OPT 25 = 5 One calculation is made with the receptor at plume
height over level terrain. A second calculation is
made as in OPT 25 - 4. The lesser of these two
calculations is used.
61
-------
TABLE 5. OPTIONS AVAILABLE IN COMPLEX/PFM
Option
Description
Technica1 Opt ions
lOPT(l) Use Terrain Adjustments
IOPT(2) No Stack Downwash
IOPT(3) No Gradual Plume Rise
IOPT(4) Include Buoyancy Induced Dispersion
Input Options
IOPT(5) Met, Data on Cards
IOPT(6) Read Hourly Emissions
IOPT(7) Specify Significant Sources
IOPT(8) Input Radial Distances and Generate Polar
Coordinate Receptors
Printed Output Options
IOPT(9) Delete Emissions With Height Table
lOPT(lO) Delete Resultant Met. Date Summary for Averaging
Period
lOPT(ll) Delete Hourly Contributions
IOPT(12) Delete Met. Data on Hourly Contributions
IOPT(13) Delete Final Plume Height and Distance to Final
Rise on Hourly Contributions
IOPT(14) Delete Hourly Summary
IOPT(15) Delete Met. Data on Hourly Summary
IOPT(16) Delete Final Plume Height and Distance to Final
Rise on Hourly Summary
IOPT(17) Delete Averaging-Period Contributions
IOPT(18) Delete Averaging-Period Summary
IOPT(19) Delete Average Concentrations and High-Five Table
Other Control and Output Options
IOFT(20) . Run is Part of a Segmented Run
IOPT(2l) Write Partial Concentrations to Disk or Tape
IOPT(22) Write Hourly Concentrations to Disk or Tape
IOPTC23) Write Averaging Period Concentrations to Disk or
Tape
IOPT(24) Punch Averaging Period Concentrations on Cards
Terrain Adjustment Options
IOPT(25) Complex Terrain Option
IOPT(26) Potential Flow Option
62
-------
Option 26iwas added to control the PFM option selection. If zero is
specified, standard COMPLEX calculations are made, with the exception that
COMPLEX I is called for stability classes 5 and 6 (E and F) and COMPLEX II
is called for stability classes 1-4 (A-D). In this sense, the COMPLEX
"host" model for PFM computations is a mix of the two previously available
COMPLEX models.
When option 26 is set to 1 or 2, PFM calculations are triggered.
Selection 1 is the PFM-Long option which is designed to approximate PFM
results in a long, sequential calculation mode. Selection 2 is the
PFM-Short option which calls PFM as a subroutine. PFM-Short is designed for
short runs where refined, critical-period calculations are needed.
Option 26 (non-zero) also triggers the layered plume rise
calculations. The user has access to this feature only in combination with
the entire PFM option.
A few other options should be specifically set when option 26 equals 1
or 2 to be compatible. Option 3, the gradual plume rise option, should be
set for no gradual plume rise because the layered plume rise equations
produce only a final rise height. Option 1 should be used because the PFM
option requires terrain correction. Option 25 should be set to 1. And
option 8, the MPTER polar receptor grid option, should not be used because a
more versatile radial receptor package is contained in the PFM option
sequence.
3.4.2 Input Format Specifications
Card Input
Tables 6 through 20 list the card input needed to run COMPLEX/PFM.
Individual values on a card must be separated by either a comma or a space
when the data format is specified as "free."
Table 6 describes card type 1, which applies to the first three cards
(see Figure 18). These identify the output run for the user's records.
Table 7 describes card type 4, which initializes several program control
variables and conversion constants. Table 8 describes card type 5, which
controls option branch points within the model. Table 9 describes card type
6, which sets the values for the surface wind measurement height, the wind
power law exponents, the terrain adjustment factors, and the minimum plume
standoff distance. Table 10 describes card type 7, the source specification
cards (up to 50).
Most of the remaining cards depend upon option selections. Table 11
describes card type 8, which is used only when significant source
contributions are specified (OFT 7=1). Table 12 describes card type 9,
which is used when standard meteorological data is read in from a file (OPT '
*= 0). Tables 13 and 14 describe card types 10 and 11, which are used when
the CRSTER type of radial receptor pattern ia wanted (OPT 8=1). Card
63
-------
Mataorology
Segmented
Run
Up to
Receptors
Grid Center
Obstacle
Information card
TVP"
13
Receptors
Polar Coord.
Rocoptora
Polnr Coord,
Rocdptor Etav. Cnrd
Type
Met Data
Identifiers
ispeciflad
Slgnlf, Sources
C.idl
jurces
Card
/
Card
Type
B
Wind ond
Terrain
Options
Control end •
ConnWnlD
Throo
Heading
Cards
Card
4
Card
6
Cord
6
Card
«»
9
Type
10
11
Card
Typo
?2
Cord
Typo
14
Card
Type
15
Card
Type
IB
If Option 5=1
If Option 20=1
/ II Option 26 = 1.2
II Option 26 =1,2
If Option 26 = 2
If Option 26 =
If Options 1 and 8 = T and Option 26=0
If Option 8 -1 and Option 26 • 0
If Opt Ion 7=1
If Option 5=0
Figure 18. Input data deck setup for COMPLEX/PFM.
-------
TABLE 6. COMPLEX/PFM CARD TYPE 1, 2, AND 3 - TITLE (3 CARDS)
Variable Format Description Units
LINE! 20A4 80 alphanumeric characters for heading
LINE2 20A4 80 alphanumeric characters for heading -
LINE3 20A4 80 alphanumeric characters for heading
65
-------
TABLE 7. COMPLEX/PFM CARD 4 - CONTROL AND CONSTANTS (1. CARD)
Variable Format
Description
Units
IDATE(l)
IDATE(2)
IHSTRT
NAVG
IPOL
NSIGP
NAV5
CONONE
CELM
HAFL
Free Two digit year
Format Starting Julian day for this run
Starting hour for this run
Number of averaging periods to be run
Number of hours in an averaging period
Pollutant indicator:
3 = S02,
4 = suspended participates
Number of significant point sources,
max =25.
Number of hours in the user specified
per?od for which a high-five
concentration table is generated,
(most likely equal to 2, 4, 6, or 12)
Multiplier constant, user distance units to
km example multipliers:
feet to km: 3.Q48E-04
miles to km: 1.609344
meters to km? l.OE-03
Multiplier constant, user ht units to m
Pollutant half-life
(an entry of zero will cause skipping
of pollutant loss calculations)
seconds
66
-------
TABLE 8. COMPLEX/PFM CARD 5 - OPTIONS (1 card)
Variable Format
Description
Technical Optiona
lOPT(l) F-'ee Use Terrain Adjustments
IOPT(2) Format No Stack Downwash
IOPT(3) No Gradual Plume Rise
IOFr(4) Include Buoyancy Induced Dispersion
Input Optiona
IOPT(5) Met. Data on Cards
IOPT(6) Read Hourly Emissions
IOPT(7) Specify Significant Sources
IOPT(8) Input Radial Distances and Generate Polar
Coordinate Receptors
Printed Output Options
IOPT(9) Delete Emissions With Height Table
IOPT(10) Delete Resultant Met. Data Summary for Averaging
Period
lOPT(ll) Delete Hourly Contributions
IOPT(12) Delete Met. Data on Hourly Contributions
IOPT(13) Delete Final Plume Height and Distance to Final
Rise on Hourly Contributions
IOPT(14) Delete Hourly Summary
IOPT(15) Delete Met. Data on Hourly Summary
IOPT(16) Delete Final Plume Hieght and Distance to Final
Rise on Hourly Summary
IOPT(17) Delete Averaging-Period Contributions
IOPT(18) Delete Averaging-Period Summary
IOPT(19) Delete Average Concentrations and High-Five Table
Other Control and Output Options
IOPT(20) Run is Part of a Segmented Run
IOPT(21) Write Partial Concentrations to Disk or Tape
IOPT(22) Write Hourly Concentrations to Disk or Tape
IOPT(23) Write Averaging Period Concentrations to Disk or
Tape
IOPT(24) Punch Averagir.^ Period Concentrations on Cards
Terrain Adjustment Options
IOPT<25) Complex Terrain Option
IOPT(26) Potential Flow Option
Note:
Integer values: 0 or 1 (0 means don't use option) on options
1-24.
67
-------
TABLE 9, COMPLEX/PFM CARD 6 - WIND AMD TERRAIN (1 card)
Variable
Format
Description
Units
HANE
PLd), 1=1,6
CONTER(l), 1=1,6
ZMIH
Free Anemometer height meters
Format Wind increase with height
exponents for each stability
class
Terrain adjustment factors
for each stability class
(real numbers from 0 to 1)
Minimum approach distance of meters
plume centerlins above the
terrain surface
68
-------
TABLE 10. COMPLEX/PFM CARD TYPE 7 - POIHT SOURCE* (up to 50 cards)
Variable
RNAME
SOURCE (1, HPT)
SQURCE<25NPT)
SODRCE(3,NPT)
SOURCE(4,NPT)
SOURCE ( 5, NPT)
SOURCE(6,NPT)
. SOURCE ( 7 ,KPT)
SOURCE(8,NPT)
ELP(NPT)
Format
2A6
F8.2
F8.2
F8.2
F8.2
F8.2
F8.2
F8.2
F8.2
F4.0
Deacription
12 alphanumeric characters for source
identification
East coordinate of point source
North coordinate of point source
Sulfur dioxide emission rate
Particulate emission rate
Physical stack height
Stack gas temperature
Stack inside diameter
Stack gas exit velocity
Source ground-level elevation
Units
.
user units
user units
g s~
-i
g s
meters
Kelvin
meters
m B
user ht unit!
*Card with ENDPOINTS in columns 1-9 is read in after the last point source card.
-------
TABLE 11. COMPLEX/PFM CARD TYPE 8 - SPECIFIED S'v.jTFICANT
SOURCES (1 card)(used with option 7=1)
Variable Format Description Units
NPT 13 Number of user specified significant
point sources
MPS 2513 Point source numbers user wants to be
considered significant
70
-------
TABLE 12. COMPLEX/PFM CARD TYPE 9 - MET. DATA IDENTIFIERS
(used with option 5=0)
Variab le Format Description Units
ISFCD Free Surface met station identifier 5 digits
ISFCYR Format Year of surface met data 2 digits
IMXD Upper-air station identifier 5 digits
IMXXR Year of mixing-height data 2 digits
71
-------
TABLE 13. COMPLEX/PFM CARD TYPE 10 - POLAR COORDINATE RECEPTORS
(1 card) (used with option 8=1)
Variable Format Description : TJnita
RAUIL(I) Free Up to five radial distances at user units
1=1,5 format which 36 receptors are generated
around points CENTX, CENTY on
azimuths 10 to 360 degrees
CEHTX East coordinate about which user units
radials are centered
GENTY North coordinate about which user units
radials are centered
72
-------
TABLE 14^ COMPLEX/PFM CARD TYPE 11 - POLAR COORDINATE RECEPTOR
ELEVATIONS (36 cards) (used if options 1 and 8
are both 1)
Variable Format Description : _ Units
IDUM 12 Azimuth indicator (l to 36)
8X (8 blank columns)
ELRDUM 5F10.0 Receptor ground-level user ht units
elevations for this
azimuth for up to five
distances
73
-------
TABLE 15. COMPLEX/PPM CARD TYPE 12 - RECEPTORS* (up to 180 cards)
(use If option 26 = 0)
Variable Format
Description
Units
RNAME 2A4 8 alphanumeric characters for
station identification
RREC F10.3 East coordinate of receptor
SREC F10.3 North coordinate of receptor
ZR F1G.O Receptor height above local
ground-level
ELR F10.0 Receptor ground—level elevation
user units
user units
meters
user ht units
*Card with ENDREC incolurans 1-6 should follow the last receptor card
(a maximum of 180 receptors are allowed).
74
-------
TABLE 16. COMPLEX/PFH CARD TYPE 13; - OBSTACLE INFORMATION (used if
option 26 =* 2)
Variable Format
Description
Units
XO Free Source coordinates as represented
YO format in hill coordinate system
HAD Obstacle relief height
OBSANG Angle from North to hill x-axis
CWARO Obstacle cross-wind aspect ratio
(ellipsoid)
AWARO Obstacle along-wind aspect ratio
(ellipsoid)
BLFSH Obstacle bluff shape number (bluff)
user units
user units
user ht units
degrees
75
-------
TABLE 17. COMPLEX/PFM CARD; TYPE 14 - GRID CENTER
(Used if optionJ26 « 1, 2)
Variable Format Description i ; Units
CENTX Free User coordinates for the center user units
CENTY Format of the polar receptor grid (must user units
coincide with source position)
76
-------
TABLE 18. COMPLKX/PFM CARD TYPE 15 - RECEPTORS* (up to 180 cards)
(used if option 26 = 1, 2)
Variable
Format
Description
Unita
RNAME 2A4 8 alphanumeric characters for
station identification
2X
RTHETA F10.0 Angle to receptor (clockwise from
north)
RADIUS F10.0 Distance tovreceptor
ZR F10.0 Receptor height above local
ground-level
ELR F10.0 Receptor ground-level elevation
HTERAN F10.0 Obstacle relief height
GDIS F1Q.O Source, distance from obstacle
center
TOBSH 12 Obstacle shape number (1-20)
degrees
user units
meters
user ht units
user ht units
user units
*Card with ENDREC in columns 1-6 should follow the last receptor card (a
maximum of 180 receptors are allowed).
77
-------
TABLE 19. COMPLEX/PFM CARD TYPE 16 - SEGMENTED RUN (1 card)
(used with option 20 = 1)_
Variable Format Descrlptioo -, , r ,-^
IDAY Free Number of days already processed
LDRUN Format Last day to be processed in this run
78
-------
TABLE 20. COMPLEX/PM CARD TYPE 17 - METEOROLOGY
(used with option 5=1)
Variable
Format
Description
Units
JYR Free Year of met data
DAY1 Format Julian day of met data
JHR Hour of met data
IKST Stability class for this hour
QU Wind speed for this hour
QTEMP Ambient air temperature for
this hour
QTHETA Wind direction for this hour
QHL Mixing height for this hour
2 digits
3 digits
2 digits
~1
m s
Kelvin
degrees azimuth
meters
79
-------
type 11 is reads only when a terrain correction is also specified (OPT
1 = 1). Table 15 describes card type 12, the £«•_? receptor card. It
contains the name, coordinates, and heights of receptors placed anywhere in
the field. If less than 180 receptors were specified by the radial grid
(OPT 8=1), additional receptors may be read in here. If the PFM option is
used, card type 12 is replaced by card type 14.
Tables 16 through 18 describe the receptor and terrain information
cards (types 13, 14, and 15) needed when the PFM option is used (OPT 26 =
1,2). Card type 13 is used only with PFM-Short (OPT 26 = 2). It contains
the location, size, shape, and orientation of the single terrain feature
under study. Note that the obstacle relief height is the height of the hill
above stack base. Card type 14 is used with both PFM options. It contains
the user coordinates for the center of the radial receptor grid. These
coordinates must be the same as those for each of the sources (card type 7).
Card type 15 is also used with both PFM options. It takes the place of the
regular receptor card (type 12) and contains receptor coordinates and
heights as well as obstacle information for the one terrain feature (if any)
closest to the receptor along the receptor radial.
Table 19 describes card type 16, which sets parameters to allow a
segmented run to be made (OPT 20 " 1). Table 20 describes card type 17,
which contains the surface meteorological data on cards (OPT 5=1).
FileInput
Table 21 shows the file structure for the meteorological data that are
read from unit 11 if option 5=0. This file is unformatted and is normally
obtained as output from MMMET or the CRSTER preprocessor.
Table 22 shows the file structure for the emission data that is read
from unit 15 if option 6 • 1. This file is unformatted and must be
generated by the user.
Table 23 shows the file structure for the hourly profile data that is
read from unit 16 if option 26 = 1, or 2. This file is unformatted and is
usually generated by the PROFIL preprocessor. If the soundings are prepared
and analysed by hand (e.g., in the case of a PFM-Short analysis), the user
must generate this file.
80
-------
TABLE 21. COMPLEX/PFM OPTIONAL Il:?UT FIJLi - METEOROLOGICAL DATA
(unit 11) (Input if option 5=0)
Variable
Diraensiona
Description
Units
Record 1
ID
IYEAR
IBM
IYR
Surface station identifier
Year of surface data
Mix ht station identifier
Year of mix ht data
Record Type 2 (one for each day of year)
JYR Year
IMO Month
DAY1 Julian day
IKST 24 Stability class
QU 24 Wind speed
QTEMP 24 Ambient air temperature
DUMR 24 Flow vector to 10°
QTHETA 24 Randomized flow vector
HLH 2,24 Mixing ht
5 digits
2 digits
5 digits
2 digits
m
Kelvin
deg— azimuth
deg-azimuth
meters
81
-------
TABLE 22. CokPLEX/PFM OPTIONAL INPUT FILE - EMISSION DATA
(tinit 15) (input if optidn 6 = 1)
Variable Dimensions Description Units
Record Type 1 (one for each hour of simulation)
IDATP -- Date-time indicator
consisting of:
Year 2 digits
Julian day 3 digits
Hour 2 digits
SOURCE
(IPOL, I), 1=1,
NPT
Emission rate for the
pollutant IPOL for
each source g e~l
82
-------
TABLE 23. COMPLEX/PEM OPTIONAL INPUT FILE - HOURLY SOUNDING DATA
(unit 16) (input if option 26 = 1, or 2)
Variable Dimensions Description Units
NLEVM - Number of levels in sounding -
for thi.s hour
HKLX - Mixing height for this hour meters
HT NLEVM Height of observation meters
WSP NLEVM Wind speed m s~l
WDIR NLEVM Wind flow vector degrees-azimuth
IMPURE NLEVM Air temperature Kelvin
83
-------
REFERENCES
Bass, A., D.G. Striraaitis, B.A. Egan 1981. Potential FlowModel for
Gaussian Plume Interaction with Simple Terrain Features, U.S. EPA
Office of Research and Development, EPA-600/54-81-008. Research
Triangle Park, NC.
Benkley, C.W., L.L. Schulman 1979. Estimating Hourly Mixing Depths
from Historical Meteorological Data. Journal of Applied
Meteorology, 18(6); 772-780.
Briggs, G.A. 1975. Plume Rise Prediction. Lectures on Air
Pollution and Environmental Impact Analysis, American
Meteorolgical Society, pp. 80-84.
Briggs, G.A. 1969. Flume Rise. Critical Review (TID-25075). Atomic
Energy Commission, Division of Technical Information,
Oak Ridge, TN.
Burt, E.W. 1977. Valley Model User's Guide. EPA-450/2-77-018. U.S.
Environmental Protection Agency, Research Triangle Park, NC.
Hunt, J.C.R. and R.J. Hulliearn 1973. Turbulence Dispersion from
Sources Near Two-Dimensional Obstacles. Journal of Fluid
MjBchanica. 61; 245-274.
Hunt, J.C.R., J.S. Puttock, and W.H. Snyder 1979. Turbulent Diffusion
from a Point Source in Stratified and Neutral Flows Around a
Three-Dimensional Hill. Part I - Diffusion Equation Anlysis.
Atmospheric Environment 13; 1227-1239.
Hunt, J.C.R. and W.H. Snyder 1978. Flow Structure and Turbulent
Diffusion Around a Three-Dimensional Hill. Part II - Surface
Concentrations due to Upstream Sources (unpublished manuscript).
Lamb, H. 1945. Hydrodynamics. New York: Dover Publications.
Milne-Thomson, L.M. 1960. Theoretical Hydrodynamics. New York:
HacMillan.
Pasquill, F. 1978. Atmospheric Dispersion Parameters in Plume
Modeling. EPA-600/4-78-021.U.S. Environmental Protection
Agency, Research Triangle Park, NC.
-------
Fierce, T.E., D.B. Turner 1980: Uaerjs Guide for MPTER - A Multiple
Point Gaussian Dispersion Algorithm with Optional Terrain
Adjustment, U.S. EPA Office of Research and Development,
EPA-600/8-80-016. Research triangle Park, NC.
Turner, D.B. 1970. Workshop of Atmospheric Dispersion Estimates.
U.S. Department of Health, Education and Welfare. Public Health
Service Publication 999-AP-26, 88 pp.
U.S. Environmental Protection Agency 1977: Uaer's Manual for Single
Source (CRSTER) Model. Monitoring and Data Analysis Division,
EPA-450/2-77-Q13.Research Triangle Park, NC.
85
-------
APPENDIX A
PFM BLUFF CROSS-SECTIONS
A.I Preparation of Terrain Profiles
Bluff obstacle shapes should be selected whenever a wide terrain
feature has a sharp ascending face and no significant descending
face. The condition that the feature be wide helps to differentiate
between a two-dimensional bluff-like shape, and a three-dimensional
ellipsoid shape which has a large along-wind aspect ratio.
If the feature has a cross-wind aspect ratio greater than ten,
then it should be considered two-dimensions! (i.e., the feature looks
as though it has an infinite extent in the cross-wind direction.) If
such a feature then has no apparent descending face, a bluff profile
should be chosen for use with the PFM option. Had the feature
possessed an apparent descending face, then an ellipsoid with a large
along-wind aspect ratio would have been preferable.
Once a decision is made to model with a bluff profile, a
representative cross-section of the terrain feature should be
prepared. Follow these five steps:
1. Select a representative portion of the terrain feature.
2. Find the height of the plateau beyond the rising face. Use
this hill height as a scaling length for horizontal and
vertical measurements.
3. Prepare a table of terrain height versus distance starting
at least eight hill heights upwind of the crest, and
continuing at least two hill heights beyond the crest.
4. Divide these measurements by the hill height to form
non-dimensional heights and distances.
5. Plot the tabulated values. Use a scale of 2 cm equals one
length unit for the horizontal distances, and a scale of
10 cm equals one length unit for the height?.
You should now have a cross-section of the terrain feature drawn to
the same scale as the model bluff profiles presented in the next
section.
86
-------
*f A,2 Selection r^pf a Model Bluff
Twenty bluff shapes are contained in the PFM model. The user
' selects one of these by specifying a bluff shape number (parameter
i BLFSH) between one and twenty.
.,-j.... Ten of these shapes are presented in Figure A-l. They correspond
to bluff shapes 1-19, odd. Each is drawn to the same scale: 10 cm
equals one hill height on the vertical axis; 2 cm equals one hill
height on the horizontal axis.
To select one of the shapes, place a tracing of the actual
terrain shape prepared according to the instructions in Section A.I
over Figure A-l. Slide the tracing back and forth along the
horizontal axis to obtain the best general fit to one of the
underlying curvea. The portion of the bluff above one-half the hill
height is most important in the matching process. Now interpolate
between the nearest curves to obtain the best bluff shape number.
Once the shape has been chosen, note where the zero line of
Figure A-l cuts across your terrain tracing. This is the "center" or
reference point that must be used in specifying distances between your
source and the terrain feature. Be sure to convert length units back
to your user units before referencing this number.
If PFM-Long is the option being used, then only four of the
twenty bluff shapes are available to you. These four are implicity
contained in Figure A-l and must be referenced by parameter "IOBSH".
IOBSH 1 is shape 1, IOBSH 2 is shape 8, IOBSH 3 is shape 12, and
IOBSH 4 is shape 16.
87
-------
00
00
1.2-
1.0-
-6.0
Figure A-l. BXuff shapes available in PFM-Short (20 in all)
Note: Vertical scale is 5 times the horizontal scale.
-------
APPENDIX B
PROGRAM DESCRIPTION - PROFIL
B.I Program Structure
Program PROFIL computes hourly mixing heights and temperature and
wind profiles from twice daily temperature and wind velocity
soundings. These soundings must either be formatted by READS6 if the
data come from NCC 5600 tapes, or they must be formatted by the user
i£ the data comes from some other source.
Major loops and all subroutines are charted in Figure, B-l. Within
-------
PROHL
! "Read Control Data
•Initialize First Day's Data
—-Loopon Days
•Read Surface and Profile Data
"REPLAC
r~' "" Loop on Hours From Midnight to Morning Sounding
"INTERP
* Calculate Mechanical Mixing Height
" INSERT
" REMOV
* Find Maximum Mixing Height (HMAX) Between Morning and Evening Sounding
" TEXTR
*" INTERP
" CONV
* Calculate Mechanical Mixing Height
" DTINT
—-Loop on Hours From Morning to Evening Sounding
•• INTERP
* Adjust Profile Below HMAX For T Change at HMAX
" CONV
* Calculate Mechanical Mixing Height
" INSERT
" REMOV
I
\ .
Loop On Hours From Evening Sounding to Midnight
" INTERP
" CONV
* Calculate Mechanical Mixing Height
" INSERT
" REMOV
3
s
End
Figure B-l. Flow chart for program PROFIL.
90
-------
APPENDIX C
PROGRAM DESCRIPTION - COHPLEX/PFM
C.I Program Structure
The COMPLEX/PFM code structure is nearly the same as that for
COMPLEX (I, II) and MPTER. Changes required to include the PFM option
have been kept to a minimum; most of the new code appears as new
subroutines.
Major loops and all subroutines are charted in Figure C-l. When
this diagram is compared with the corresponding charts for MPTER or
COMPLEX (I, II), the major loops are identical. New branches within
these loops are activated when the PFM option is turned on (OPT 26 =
1,2). A brief description of the function of each subroutine follows
(^denotes subroutines unique to PFM options).
COMPLCX/PFM -
Main Program: reads and checks input parameters; calls
POLAR for ?FM receptor information; initializes
variables; writes out initial information; determines
significant sources; controls loops on days, averaging
time, and hours; calls PTR for point source
concentrations; reads and writes tapes/discs; calls
OUTHR to print concentration contributions and
summaries.
POLAR *- Subroutine called from Main, POLAR handles receptor and
terrain input when PFM option is selected. Standard
terrain information as well as PFM terrain description
is obtained. If PFM-Long (QPT26=1) is being executed,
PFMFAC is called to calculate PFM factors.
PFMFAC *- Subroutine called by POLAR, PFMFAC returns PFM factors
for five plume heights and six Froude number classes at
each receptor.
ANGARC - Function called from Main, ANGARC calculates the arctan
of east resultant wind component over the north
resultant wind component and returns the angle between
0" and 360°.
PTR - Subroutine called by Main, PTR calculates the point
source contribution at each receptor. It controls the
receptor loop and the source loop.
91
-------
COMPLiX/PFM
1
L_
Read Input Data
"POLAR (OPT 26 = 1,2)
1 •" PFMFAC (OPT 26 = 1)
J—- Loop for Calendar Days
- Loop for Averaging Time
* Read Standard Met Data
"ANGARC
Loop on Hours
• Read Hourly Sounding Data (OPT 26 = 1,2)
»*«** PTR
* Calc Max Hill Height Along Wind Direction (OPT 26 = 2)
- Loop on Receptors
««HCCALC(OPT26 = 1,2)
- Loop on Sources
•*NRISE(OPT26 = 1,2)
«»SRISE(OPT26=1,2)
»«FRNUM (OPT 26 = 1,2; KST > 3; H > HC)
"CMPLX (OPT 26 = 0) or (KST < 4) or (H < HC)
i
. *• RCP 2 (KST < 4} or (KST = 4; OPT 26 = 0)
!
i
PGYZ
1 " RCP (KST > 4) or (KST = 4; OPT 26 = 1.2)
' ** PQZ
"» RCP 2 (OPT 26 = 1.2; KST > 3; H > HC)
"CSTFAC (OPT 26 = 1)
*«WCFAC(OPT26 = 2)
SPFM
PGYZ
" Rank
** OUTHR
•" OUTAVG (Entry Point in OUTHR)
Exit
Figure C-l. Flow chart for prograa COMPLBX/PFM.
92
-------
RANK - Subroutine called by Main, RANK ordered concentrations
for four or five averaging times so that the five
highest concentrations are presented for each receptor.
OUTHR - Subroutine called by Main, OUTHR arranges and prints
tables of pollutant concentration.
HGCALC *- Subroutine called by PTR when PPM option is used,
HCCALC computes the critical dividing streamline of the
flow.
NRISE *- Subroutine called by PTR when PFM option is used, NRISE
calculates the neutral plume rise in the presence of
wind shear.
SRISE
PRNUM
CMPLX
RCP
RCP2
CSTFAC
WCFAC
*- Subroutine called by PTR when the Ptii option is used,
SRISE calculates the stable plume rise in the presence
of wind shear and density stratification.
*- Subroutine called by PTR when the PFM option is used,
FRNUM returns the Froude number controlling the flow
above the critical dividing streamline.
- Subroutine called by PTR when a PFM calculation is
inappropriate, CMPLX adjusts the plume heights
according to option 25 and calls either RCP2 or RCP to
perform a COMPLEX II or COMPLEX I computation,
respectively.
- Subroutine called by CMPLX, RCP returns values of az
and relative concentrations computed with sector
averaging (COMPLEX I) RCP calls PGZ for az.
- Subroutine called by CMPLX or PTR, RCP2 returns values
of relative concentrations computed without sector
averaging (COMPLEX II) and ay and oz. When
PFM-Long is used, RCP2 calls CSTFAC for PFM factors.
When PFM-Short is used, RCP2 calls WCFAC for PFM
factors.
RCP2 calls PGYZ for 0y and az.
*- Subroutine called by RCP2, CSTFAC obtains approximate
PFM factors for each source-receptor combination by
interpolating the factors computed by PFMFAC. CSTFAC
is used only when the Plftl-Lo-ig option is selected (OPT
26=1)
- Subroutine called by RCP2, WCFAC controls the call fo
subroutine SPFM and selects PFM factors for each
receptor. WCFAC is used only when the PFM-Short (or
"worst-case") option is selected (OPT 26-2).
93
-------
SPFM *- Subroutine called by WCFAC, SPFM performs the potential
flow computations for each specific source-terrain
geometry (one obstacle) for this hour. SFFM factors
are returned at 49 fixed points along the wind
trajectory. (More details are presented in Appendix D.)
PGZ *- Subroutine called by RCP, PGZ calculates oz for a
given downwind distance arid stability class using the
Pasquill-Gifford curves.
PGY2 - Subroutine called by RCP2, PGYZ calculates CTy and
Oz for a given downwind distance and stability
class using the Pasquill Gifford curves.
94
-------
APPENDIX D
SUBROUTINE DESCRIPTION - SPFM
D.I Program Structure
The SPFM computer code has six primary stations. Section. 1
performs data entry and scaling functions; Section 2 performs the
appropriate streamline trajectory computations; Section 3 alters the
computed streamline if the Frouds number is sufficiently small;
Section 4 performs the line integrals described in Section 2.2.1;
Section 5 computes the distance from the plume centerline to the
terrain surface; and Section 6 computes the PFM factors at 49 points
along the plume trajectory. Each section is briefly described below.
Data Entry
SPFM version 4.0 is designed to be called as a subroutine within
the COMPLEX model described in Section 2.3. Therefore, some of the
option switches and variables are set internally and the user no
longer has access to them. Their function is still described by muans
of the comment cards in the listing so readers who are interested
should consult the code.
Most source and obstacle information is passed to SPFM through
common PFMTER. The source location (XO, YO) is measured from an
origin centered on the terrain obstacle of height HAD. The x-axis
points along the obstacle axis closest to the wind direction if it is
an ellipsoid or if it is perpendicular to the bluff face. The bluff
shape is specified through parameter BLFSH, while the ellipse shape is
specified by the crosswind and alongwind aspect ratios (CWARO,
AWARD). WNDANG is the angle the wind makes with the x-axis.
The remaining source and meteorological information is passed
directly through the SPFM subroutine call. NSRCE is the source
identification number (no larger than 50), HO is the plume height, and
FROUDE is the Froude number.
All length scales are divided by the obstacle height because SPFM
computations are nondimensional. Default constant diffusivities are
specified, but their actual values are unimportant because the PFM
factors are independent of them. And the 49 PFM receptor points are
selected along the wind direction. Their dimensional spacing (DELINT)
is passed back to the calling program through common PFMTER,
95
-------
Plume Path and Streamline Velocities
The plume centerline trajectory is defined numerically by placing
up to 600 streamline points in an array PATH. Along-streamline
velocity components for each of these points is placed in an array
VEL. If the bluff shape parameter BLFSH is zero, then subroutine
ELLIPS is called and the appropriate streamline path over an ellipsoid
characterized by CWARO and AWARO is computed. When BLFSH is
non-zero.subroutine BLUFF is called and the appropriate streamline
over the specified bluff is computed.
Stratification Adjustment
PATH and VEL array elements are adjusted to approximate
first-order effects of stable density stratification when the Froude
number is less than or equal to ten. Again, either a bluff subroutine
or an ellipse subroutine is called depending on the type of obstacle.
Separate routines are needed because the local streamline height
required in Equation 2—53 depends on the shape of the ground surface.
Subroutine BLFSTR is called for adjustment of bluff streamlines, and
ELLSTR is called for the ellipsoid.
Plume Path Line Integrals
Subroutine LININT is called to evaluate the line integrals from
the source to each PFM receptor location. These line integrals
determine the value of the $ and T integrals described in Section
2.2.1 and the elapsed time. LININT als-J obtains the effective radius
of curvature of the ellipsoidal obstacles and substitutes this for the
distance to the axis of symmetry. (Recall that the theory of Section
2.2.1 was developed for axisymmetric three-dimensional obstacles.)
Also, because the PFM receptor points do not coincide with all
elements in the PATH and VEL arrays, LlillNT interpolates the y and z
coordinate for each receptor point x, and it also interpolates the
along-streamline velocity.
Plume Distance From the Surface
The effective plume centerline distance from the surface of the
obstacle is not necessarily the vertical coordinate of the plume
streamline minus the local terrain elevation. Rather, it must be the
distance to the surface along the normal to the streamline because
this is the direction in which Oz is defined.
Two separate subroutines are required for this calculation
because the surface is defined differently for the bluff than for the
ellipsoid. The value of BLFSH once again determines which is called.
The bluff subroutine is BLFDIS and the ellipsoid subroutine is ELLDIS.
96
-------
PFM Factors
Stored array elements from the line integrations and the
centerline displacement above the surface at each PFM receptor are
combined to form CFAC, SYFAC, SZFAC, and HFAC arrays (see Section
2.2.2). Each factor array has an element for each receptor point*
allowing either exact concentration computations at these points, or
interpolated calculations for points in between.
The arrays are also dimensioned for up to 50 individual sources
because they are passed back to the calling program through common
FACWC. This means that SPFM as a subroutine must be called for each
source for each hour that is simulated. The advantage to keeping SPFM
a subroutine is that the calling model that has been selected is
already established with a well-defined input and output structure
(see Section 2.3). Details of averaging times for reporting
concentrations and stability class determinations are external to
SPFM, so its structure can remain compact and uncluttered. In this
form new users may readily become familiar with it.
Figure D-l also shows the inter-relationship of all subroutines
within SPFM. The function of each is summarized below.
SPFM - Subroutine called by WCFAC, SPFM is the main calling program
for the PFM branch within COMPLEX/PFM. It controls the
transformation and scaling of the input data, provides the
structure for the remaining five major sections of the code,
and performs the SFFM factor computations.
BLUFF - Subroutine called by SPFM, BLUFF is responsible for
computing the plume path and velocities along the path for
flow over a two-dimensional bluff. It scales the problem
into "step space", computes the plume streamfunction through
a combination of calls to BLFHT and GETPSI, and then fills
up the "path" and "velocity" arrays through BLFPTH to define
the plume trajectory, and the along-streamline velocities.
BLFHT - .Subroutine called by BLUFF, BLFDIS and GETPSI, BLFHT returns
the height of a streamline over a step given the value of
the streamfunction, and the horizontal distance from the
step.
IMRAF - Subroutine called by BLFHT, IMRAF uses the Newton-Raphson
iterative technique to solve for the real part of the
parametric variable "t" (see Section 2.2.3).
GETPSI - Subroutine called by BLUFF, GETPSI returns the value of the
streamfunction (psi) that passes through a given point.
BLFPTH - Subroutine called by BLUFF, BLFPTH computes the plume path
and the velocities along the plume once the streamfunction
for the plume streamline is found. It solves the
transcendental equations at up to 600 points along the
streamline through repeated calls to IMRAF.
97
-------
SPFM
Section 1
Section 2
Section 3
Section 4
SactionS
Section 6
•Read Input Data
• Compute PFM Receptor Spacing
* Transform User Units to Nondimenslonal Units
• Obtain Plume Path Coordinates and Wind Speed Along Trajectory
• BLUFF (1BLFSH>0)
* Adjust Distances end BIuH Shape for Wind Angle
* Find OH. Surface Streamfunction, Plume Streamfunction
; " SLFHT " GET PSI
'"IMRAF ''"BLFHT
"" IMRAF
•BLFPTH
• - Loop Fo>- Downwind Distance
"IMRAF
* Increment Distance
ELUPS0)
• - - Loop over PFM Rsceotora
• XINT5R
•BLFHT
' " IMBAF
ELLOIS (IBLFSH * 0)
Loop over PFM Receptors
XINTER
L_
• Calculate CFAC. SZFAC, SYFAC, HFAC
Exit
Figure D-l. Flow chart for subroutine SPFM.
? 98
-------
I'fELLIPS - Subroutine called by SPFM, ELLIPS is responsible for
r computing the plume path and velocity along the plume path
for flow over an arbitrary ellipsoid. It controls the
specification of up to 600 points along the streamline, and
solves for the net velocity at these points due to mean flow
components along each of the two horizontal ellipsoid axes.
XRAPH - Subroutine called by ELLIPS, XRAPH returns the local
velocities at a point along the streamline due to the mean
flow component along the "x-axis" of the ellipsoid.
YRAPH - Subroutine called by ELLIPS, YKAPH returns the local
velocities at a point along the streamline due to the mean
flow component along the "y-axis" of the ellipsoid.
GETLAM - Subroutine called by ELLIPS and XRAPH, GETLAM solves a cubic
equation for the ellipsoidal coordinate lamda.
LAMNEW - Subroutine called by XRAPH and YRAPH, LAMNEW returns the
ellipsoidal coordinate lamda by using thu Newton-Raphson
iterative technique, given an initial value for lamda*
POTINX - Subroutine called by XRAPH, POTINX sets up the integration
needed to find the velocity potential due to the mean flow
along the "x-axis."
POTINY - Subroutine called by YRAPH, POTINY sets up the integration
needed to find the velocity potential due to the mean flow
along the "y-axis".
DQATR - Subroutine called by POTINX and POTINY, DQATR performs a
numerical integration of either the functions Fl, F2
(POTINX), or the functions F3, F4 (POTINY).
BLFSTR - Subroutine called by SPFM, BLFSTR alters the plume path and
the along-plume velocities for flow over a bluff to
approximate the effect of density stratification.
ELLSTR - Subroutine called by SPFM, ELLSTR alters the plume path and
the along-plume velocities for flow over an ellipsoid to
approximate the effect of density stratification.
LININT - Subroutine called by SPFM, UNINT evaluates the "phi" and
"tee" line integrals along the plume. These are evaluated
at 49 downwind distances.
KERML -
Subroutine called by LININT, KERNL returns the value of the
"phi" and "tee" integral kernels at a point and also
computes the advective time and distance along the plume
streamline.
99
-------
DIFF - Subroutine called by LININT if PGT-scaled diffusivities are
requested, DIFF returns dimensionless diffusivities scaled
from PGT sigraa curves.
PGT - Subroutine called by DIFF, PGT returns aigma-y and sigma-z
given a stability class and a downwind distance.
BLFDIS - Subroutine called by SPFH, BLFDIS returns the distance from
the plume centerline to the surface along the normal to the
streamlines at each of 49 receptor points.
ELLDIS - Subroutine called by SPFH, ELLDIS returns the distance from
the plume centerline to the surface along the normal to the
streamline at each of 49 receptor points.
XINTER - Subroutine called by KERNL, BLFDIS, and ELLDIS, XINTER
interpolates between points along a streamline to return the
y,z coordinates, and the corresponding velocity components
at an arbitrary downwind distance, x.
100
-------
APPENDIX E
TEST CASE
The computer code for COMPLEX/PFM is available on version 5 of the
User's Network for Applied Modeling of Air Pollution (UNAMAP). It con-
sists of two files. The first file, COMPLEX/PFM file, contains the main
program and subroutines, preprocessing programs and test data, and input
data streams for test cases. The second file, COMPLEX/PFMFACTORS, contains
a formated table of PFM factors. Preprocessor program SETUP converts
this formatted file to an unformatted random access file for use with
COMPLEX/PFM-Long.
Test cases are provided for the meteorological preprocessors and for
both PFM-Long and PFM-Short options of the main program. There are two
meteorological preprocessor programs - READ56 and PRQFIL,. READ56 reads an
unformatted TDF 5600 upper air file and produces a formatted file of
twice-daily soundings as specified by the program control variables in
Table 2. A segment of a formatted TDF 5600 file has been included in
COMPLEX/PFM along with program FTOUF-UP which converts these data to an
unformatted file for input to READ56. This is the starting point for the
test case and the user's output should be compared with the test output of
READ56 listed in the test output file of UNAMAP before proceeding.
Program PROFIL requires the formatted output of READ56 and the unfor-
mattted surface meteorology output of the preprocessors RAMMST or CRSMET.
The surface meteorological preprocessing steps have been omitted here
because the code for performing this task (RAMMET or CRSMET) is contained
elsewhere in UNAMAP. However, a segment of the processed formatted surface
data file for use in the test case has been provided in the COMPLEX/PFM
UNAMAP file along with program FTOUF-SFC which converts the data to an
unformatted file as required by PROFIL. The output of PROFIL is also
unformatted. A formatted version of PROFIL output for the test case is
included in the COMPLEX/PFM file along with program SND. The user should
format his version of the PROFIL output file using SND to verify that the
numbers match.
Four individual test case executions of the main program are Included.
They represent computations for flow over a bluff and flow over an ellipsoid,
using both PFM-Long and PFM-Short options. The cases are needed to verify
that the PFM look-up table Is constructed correctly (PFM-Long), and to
verify that SPFM computations (PFM-Short) for bluffs and ellipsoids are
being done correctly. Input data streams are contained in the COMPLEX/
PFM file and output from running the tests are contained In the test out-
put file of UNAMAP.
101
-------
Date
Chief, Environmental Applications Branch
Meteorology and Assessment Division (MD-80)
U.S. Environmental Protection Agency
Research Triangle Park, NC 27711
I would like to receive future revisions
to the User's Guide for COMPLEX/PFM
Name
Organization
Address
City State Zip
Phone (Optional) ( )
102
------- |