-------
For stable conditions or unlimited mixing:
* (A12)
where g~ was defined at the beginning of this appendix.
In unstable or neutral conditions and if a is greater than 1.6 times
the mixing height, L, the distribution below the mixing height is uniform
with height provided that both the effective height, H, and the receptor
height, z, are below the mixing height:
f = 1/L (A13)
(If H or z is above the mixing height, f = 0.)
In all other unstable or neutral conditions, that is, if a is less
than 1.6 times the mixing height:
f = g3/ [ az(2v}h ] (A14)
LINE SOURCE COMPUTATIONS
The line source and receptor relationships are shown in Figure 4. Points
A and B are the beginning and ending points of the line, x and y are the
upwind and crosswind distances from a receptor to a point on the line source.
x and y are given by equations (AT) and(A2).
66
-------
R and S are the coordinates of any point on the line and are functions
P p
of length, A, along the line and total length of the line, D:
RP 4 (RB - V + RA (A15)
SB -V + SA
Since R and S are functions of a (see Figure 2), x and y are also functions
of a.
The upwind distances, x, of the receptor to points A and B are calculated.
If both of these x's are negative, then the line source is downwind of the
receptor and the concentration at the receptor due to the line source is
zero. If either or both x's are positive, it is determined if there is a
point on the line source directly upwind of the receptor. The equation for
the line source is:
(R - RA) / (S - SA) = (RB - RA) / (SB - SA) (A17)
The equation of the line through the receptor in the upwind direction is:
R - R = (S - S ) tan 0 (A18)
Solving these two equations simultaneously for R and S by letting:
m= (RB-RA)/ (SB-SA).
m SA - Sr tan 0 + Rp - RA (A19)
m - tan 0
67
-------
and R = (S - S ) tan 0 + R (A20)
This resulting intersection R, S (see point I, Figure 4) must be tested
to see if it is upwind of the receptor by testing for positive x using
equation (Al), and if it lies on the line source between points A and B.
If both of these conditions are met, the line source is considered to be
made up of two segments: The first segment from point A to the intersection
point, and the second segment from the intersection point to point B. The
point directly upwind has a significant contribution to the concentration
at the receptor. The line is broken into two segments so that this point
will be considered in the integration estimate that follows. If there are
two segments to the line source, the total concentration is given by the
sum of the two concentrations from the segments.
The concentration from a line source is then given by the following
equation.
q rD
X = TT- / f. d* (A21)
where:
q = line source emission rate, g sec m
X/
u = wind speed, m sec
D = line source length, meters
and f = point source dispersion function
68
-------
The point source dispersion function is given by one of the three
following equations where a and a are evaluated for the upwind distance,
x, and the stability class.
For stable conditions, or unlimited mixing:
f =
9192 (A22)
2-7T a a
y z
/—••
where g, and g? have been previously defined.
In unstable or neutral conditions, if a is greater than 1.6 times the
mixing height, L (meters), the distribution below the mixing height is
uniform with height regardless of source or receptor height, provided both
are less than the mixing height.
f = - 91 (A23)
P a ify^\%
In all other unstable or neutral conditions:
99
z
193 (A24)
°
It should be noted that in PAL, initial a 's and oz's, such as to account
for initial dispersion in the turbulent wake behind vehicles, are input values.
The virtual distance for the given stability class is calculated within the
algorithm. These virtual distances are added to the physical distances prior
to determining a and a for each concentration computation.
69
-------
The integral, A21, is evaluated using the trapezoidal rule by making <
first estimate, C,, given by:
c =V^Jl C fp(o) + fp(lOA£) ] + I f (TA£)l (A25)
1 U \J i = 1 J
where AS, = D/10, and f is defined above. For each of the 11 evaluations of
f , the upwind distance, x, of the point on the line from the receptor is
determined as a function of Ad,and x is used to determine a and a in the
x y
function f . Portions of the line where f equals zero are eliminated and
the line is redefined. Then equation (A25) is again evaluated.
A second estimate, C-, is determined from:
2u
10
z M-Ai/2 + JAA) 1 (A26)
. = ,
J *
P
A calculation is made of (C? - C,)/C«. If this value is within a set
criteria, x = C?' ^ not w""^hin the criteria, a third estimate is made by
evaluating f for additional points midway between those points previously
calculated and assuming linear concentration change between adjacent points
(trapezoidal rule). This is continued until the approximation to the integral
converges, that is, until two successive estimates are within the criteria.
The last estimate is the concentration, x» fi"°m this line source.
70
-------
For line sources with an intersection point upwind from the receptor,
the contribution from the second segment must also be determined. If the
point directly upwind was not used to break the line into two segments,
and the entire line source is considered as one segment, it is possible
for the entire portion of the line that contributes significantly to the
concentration to be between two of the eleven points, and thus an erroneous
concentration of zero would result.
For line sources with multiple lanes, the integral would be evaluated
for each lane and the concentration summed to represent the total concentration
at that receptor from the line source.
CURVED SOURCE COMPUTATIONS
Figure 5 shows the curved source and receptor relationships. Points A
and C are the end points of the curved path. Point B is an arbitrary point
on the curve between A and C. The curve is assumed to be of constant radius.
The radius of the circle, of which the curve is a part, is determined
in the following way: The coordinates of point D midway between A and B
on a chord are found by averaging the R and S coordinates of A and B. Similarly
the coordinates of point E midway between B and C on a chord are found. The
perpendicular bisector of a chord of a circle will pass through the center
of a circle. Also the slope of the perpendicular bisector is the negative
reciprocal of the slope of the chord. The slope of chord AB is:
71
-------
(RB - RA) / (SB - SA)
and the slope of chord BC is:
(Rc - RB) / (Sc - SB)
The slope of the bisector through D, call this m,, is:
- (SR - S.) / (RR - R.) = m,
DM DM I
The slope of the bisector through E, call this m?, is:
- (Sr - SD) /(Rr - R0) = nu
U D U D L
The resulting equation of the bisector through D is:
R - R= m, (S - SJ
0 D 1 0 D
and the equation of the bisector through E is:
R0 - RE = m2 (SQ - SE) (A28)
R and S can be determined from these two equations, for example:
So = (mlSD ' m2 SE " RD + RE) ' (ml " IT12) (A29)
and RQ = m] (SQ - SD) + RD (A30)
72
-------
The radius of the circle can then be determined from the coordinates of, the
center (R , S0) and any one of the three points on the curve, for example:
S\ L. , t r\ n \ •—
. ... M J + (RA - RJ
f\ \j r\ \J
In order to provide a proper estimate-of the concentrations from the
curved source at a receptor, it is desirable to determine if a line in the
upwind direction from the receptor intersects the curve. The equation of
the line through the receptor in the direction of the upwind azimuth is:
R - Rr = (S - Sr) tan 0 (A32)
The equation of the circle of which the curve is a part is:
p2 = (S - So)2 + (R - RQ)2 (A33)
The two possible intersections of the line and the circle are found by solving
a quadratic equation for one of the two variables, for example, of the form:
a R2 +b R + c = 0 (A34)
2
where a = 1 + cot 0 (A35)
b = -2 Rr cot20 + 2 Sr cote -2 SQ cote - 2 RQ (A36)
and c = Rr cot 0 - 2 Rr $r cote - 2 R S cot© - 2 SQ S^
o o o 9
J_C^-LC^-_LI">'- £• tl\?~J\
+ $r + SQ + RQ -p (A37;
73
-------
Two possible values of R are found from:
(A38)
2
If b -4ac is negative, R has imaginary roots and there is no intersection
2
of the circle and the line. If b -4ac is zero, the line is tangent to the
2
circle and there is one intersection. If b -4ac is positive there are two
intersections and both roots of the equation are found. If one or two
values of R are found, then one or two values of S are found from:
S = (R-Rr) cote + Sr (A39)
If one or more intersections were found for the line and circle, the
upwind distance of this (these) point(s) from the receptor can be determined
using equation (Al). If this value is negative, indicating the point is
downwind from the receptor, it need not be considered further. However,
for a point with a positive x, it must be determined if this intersection
on the circle is also on the curve, that is, between points A and C.
The direction from the center of the circle (R , S ) to a point (R ,
S ) on the circle is:
* - tan"1 [ (Rp-RQ) / (Sp - SQ) ] (A40)
-------
The directed azimuth from the circle's center to any intersection and to
points A, B, and C can then be found. With some logic to consider the
crossover point (from 360° to 0°), it can be determined if the for the
intersection is within the range ot $ swept out by the curve.
After determining if any intersections are on the curve, the curve is
considered in three pieces if there are two intersections, two if there is
one intersection, and as a single curve if there are no intersections. The
computation is done very similarly to that for line sources. Whereas the
length along the line is used to make calculations at intervals for the
line source, the variation of the azimuth from the center of the circle,
<|>, is used for the curved source. The coordinates of a point on the curve
for which the concentration contribution is being evaluated are found
from:
R = RQ + p sin <|> (A41)
S = SQ + p sin (A42)
The upwind and crosswind distances, x and y, of a given point on the
curve with coordinates (R , S ) from the receptor (R , S ) are calculated
as previously using equations (Al) and (A2). The same point source dispersion
function, f , as given in the "Line Source Computations" section is used.
75
-------
Portions of the curve that contribute no concentration to the receptor
are eliminated, and the curve is redefined. An initial estimate of the
concentration from the remaining curve is made using eleven equally spaced
points along the curve. These and the ten pointss in between are then used
for a second estimate. Comparison is made of the two estimates and, as with
the area and line source estimates, if these two are within the given criteria,
the concentration is that of the second estimate. If additional estimates
are made, concentration contributions for additional points, midway between
these points already considered, are determined assuming the concentration
varies linearly between adjacent points, that is, trapezodial integration.
If the curve was broken into pieces because of intersections, the
contribution from the additional pieces must be determined as above.
76
-------
« I- E
M Q UJ
a H i~
I- O I" fsl
O M UJ
H 3 i:
I ^
to x:
«< ISI —
i: I
u
H
I/I
3 1
u
a
U
z
.J
CO
O
a.
K
o
u.
a
o
o
: z
u
(9
t- s:
i/i a *
UJ C —
C I- —
O I/I z;
o
u
l/l
J _f
O U.
a: —
sc u t/i
1/1
UJ
a
o
o
o
o
T2
<
UJ
«T CJ t/>
O O X
1/1
<
UJ
.
or or -5:
uj o *•
a a
,f~ a
a N
(o a
u> r--
,i-l C
o
a
c
^ h-
•* X
U U
u LJ
z x
n X
z t-
t/1 -* • • » •
< M 3C r-l r-t i-H rH
x: — i
O C O C.
C- O C O
x *- o o o a'
i H Z • * t *
M ' a a o o
o u* to u? a>
1/1 a
o;
u o
O'
CL H ~
O C.
tjl UJ
C, C
IB U)
o
z
.
Ul < I- fl
U S cr
- UJ L .
I = E
30. -
: ^J M H-
O UJ iJ
M »- CC.CC
UJ »- t/l • » • •
u>u u i a
z tr o uj
o
a
Q.
o
c
c
<
2
o
Ui Ul
u o
tr 3
= -J
O L-
v r
Z <
M U
o a
a. *'
uJ cc
> 3
CO
o; z
o IH
X -J
U LT
IM
U. '_J
w'* Z
1/1 3C —
t^ uJ t— O
Z O o UJ
IH a z LI
V- O 3 Ljt X
& o c o
z w »- -,
IT
o
o
ul (j o r:
1/1 or a; i: I
« 3 L. o
o: o
z
o c
o o
n •=:
t- a o o o
I i «• o o c o
z' ujh-s: c;oco
u.1 u o i c c "3 n^
azarzu c u. tr. o
MI3U-UJ f .H.-IOJ
M'_JOIXH C.OC.C;
, t- x o o a c
M. i/l U • . . •
o o
z
77
-------
1
1
I 1 1 ' ! 1
i 1
1 i 1 1
! 1
1 1
;
i i
i
— |
i
i '
!
I
1
I
1 i
1
I
i
i
i
i
o o a a
i \ft in u> m,
) . • * .
1 rH r4 *4 tH
^— *
O
LU
^
Z
i i
i
o o o o
C 0 O O
i l 1
j ;
1 -T
1
I
s ! J
a a a d a
< 3: ui a
H K h- •
Q a uj
uj H X
| r a —
i i
1 !
1 1
t
(
1
\A U
a: a
' ' U T U •
f < h- H- a
K Q UJ CM
•gas
1
-l t-l >-
t ' ^ ^
' 1 U)
! M
i M
' !_i
' ' <
; i (_(
o
c.
i ^— *• •
M >- x rn
21 w
i , H
|
1
s
1 r-<
1 X *» •
; '. *-• x
, • • • : cc «
W M M Kl
1 i
1 j
' co ad**
i1" *1
It
•
H
llli1!
, 1
i ' i i
! l : ; i i
i ' : i i 1 ; ' ' t
t ' i
i ! 1
f
i
i i
i i
• i ; '
! 1 ' i
iii '
!i i >
i
|
, ; i
,i i -j
1
1
1
1
• ! ' !< * !
i ; .MHO
: : o < a
1
1
1
1
i
i
i i
1
1 I
1
|
O O O i !
000
• . .
to to to
« Jrt
o o o
ff-l tH *-«
t-4 r4 «-l
• * •
1
z a o o o
M r- r»- h» **•
I i o ' a
o o a a
i-l H H .4
P ! H H M H
™ ' • f t *
g } M H H H
O
i
^ | 1 f
•—
M «0 M r»
1 u> u) r»- r*
I 1 1 ' Q O Q O
j ; t t " t
rf\ I i
S i :
a. H-'— •
; vil i
; i
O C O O
t t • f
' 1
1 1
, 1
! |
COCO
V V
f
}
i
1 [
i z'
1 W H
, U
! a.1
i T O
> o
*
**
555
. . .
oca
H H *4
•H H H
• •/ •
, !
O (_> a u
K) M K» ro
't- U I'- 0
H X •
u> ie
& ^ *™ »^ « r^ r^
I Z CVX • • • •
:0 ; i3!;S
1 1 1
!N 1 L
H Ul H'H O
Ul O X CC O
0: z ce o u o
' W 3 M t- »
o a a D ) Q-ioujuj
a o o o ' u1 xiac
o c o o ; .x
j O O O O '
4 V. a3 UJ O
1 O C O r4
1 O C O C
| !
—
o ' ' c
1 x -c
! aooo w uiK-2:c3
• ••• ' LJCJOlQ
i > Z K Z O O
' ! H 3 UJ LJ r
cz _j o o: n a
; to »~ s. a
'in^r-* j 3 i^ot
1 • ' ' '***
0 '
i '
1 1
a o a
o o o
a o a
• • •
1 1
J |
; !
i
1 j
i
i
:
i
I
11
t
1 '
i
i
i
i
1 :
!
i
i LJ CU O
a. a ,
1 • VI 1 •
! H i
i
i ' J
<
H ui o
— O Z Q
! o ul H U
j ' '0. _l 0
i ZJ"! I • 1
i i * r i
l ' i-* i i
C3 Q '
i UI X O
, J U, >HO
i ' o cc «; a ,
| ! 1=3 ft-0 !
!
1 1
1
1 1
I 1 I/»,U j • i
J»: 1 H
d 1
1 • w,_l !
I h- «t
1 O Z O
! < o ul o
1 Ct M Z C
' U. M M O
; - tt J •
, O H
i "1 * , 1
1 z; '
1 O
! 1 l HO
1 i 1 ! t— t a
i ' i «t ui a
M' K. O
1
(
!
I
|
I
1
1
i x
(
•
'
i
i
!
j
1
-J
<
K-
H
< W
M I
o u t-
I ,
i ;
.
.
.
;
l
'
i
]
I
',
I I
[ j
I |
] |
i ;
ro m to f
Q o a c
'lilt
h- IB CO M
r- co u L
I
M rn »o M
o d a a
i I1 1 1
a- r* m CM
ui U1 >> CO Cl * M
M Jf CNJ fH
j
i i
M «H H
•
ooocaooa
c c o c
a c o a
ccacoooa
U.Q.O. 'CDocopoa
1
i
l X
l/l
,
_J
** W
H UJ
i
i
1
c o o c
C O 0 C
i
c a a c
o o o c
t ouz ooooooca
' ' KL « • u. u r->
' <, I-i L. Q. J
1 i
i !
!
(
|
1
|
i
•
i
1
I
:
>
1 ' 1 1
1 -J; 1
«tl 'O
z! t- O
CC Z O '
' 3, H O
M, O «
i Oj Q.H •
i i
* ** •
co. a •
f! •£ i
1 i J '
1 1 i 1
1 '
1 i
i si i . ;
[ 1 1 'n
1 , o! 01
ui: fj
! i °i i
• •-;
i i
i
• - ' * ;
j
i
i
1
1
: 2: a
! ' - io
; o
-J, H
• H '
i i >
• J
1 — a; a Cv a o o a a > , \tr
i M
1
i
i
i
I i in!
1 V
!
1
I
1
!
! X
H
1
'Q
UJ 11
> x
cr cccccc
' * a
o c
1 1
o a oo f-
C C « r-
a a o a
c c c c
1 occ»- COa>r^OCOC
: u, u a.
1
,-» i
<
1 Z
1 o
1 M to
a X »-< u
o u: z
* or
, u
.
1
VI
a
lo
[V—
jo.
UJ
o
u f-i
x -i
Ul
X <
O Ul
K K
U. <
I
n
UJ a: t-
! f Z
K tt tt t— '
UJ U_ 0
h-
^
. (L
lt
i *.
,K t
' O M V-
:«t *-
LJ X
1 . i ton o
Si . >
i '» ui
! i A'
t c o c. M c a c'o i • s!i o
! ,*: rj^^-r-ioi^ioH1 ' x •
*— HOOQOOOO , * UI
a:
i * "i
o a o a
l>3 . .
>C9 ! l
!o I :
i
a lj 3| •
! • i *: h-M*»u>a-*a-tfi UJi in
l
0 — r-lrfCC". NMMN OQ1
K
a. a
lUI
1
| :
j
' . ' i
•a<
' ' ^ i i
I Oiu ! '
x!
! UK;
Q H H i O H N l»% * Ifl.U) h- 00 h- • ;
ooo • i . ; i 1 • H
o a o ' u o
t • •
i
1 ' Z1
: uj o
1 : : z
' ,« i , j ! } x I • .
1 '. ! ! ! i ' !
o c r- *
r* f\
c c o a
(M n
H H »-* C
a c r'i d a o H «i
c— ^ *^ i: c; t-< cs
o a a K
a a M »O iO K
O O O C
1 t 1 1
f^ f^J *O C
r— ?— co .-
00 Ul LD (N
,tO 3- H r-
1
co in tn 4
.3 13 3 £
1 1 1 1
tn a tsj r
r-l LO CM »-
M t" U) f
N I"! If) t-
I
1
o c a c
o o M a
Q a » H
i
1
i »o io1 ro ^
> a o a c
till
* r- M (\
* \jft m r- «H
j to a * M
to * «^ K
p1
i
a o a c
A G c c: c
t a c a c
i a a a c
i
i
l
1
" C 0,0 C
oouaoaoc:
| U Z U1 |
i^ct
! 0
;z u
i-i x,
CC
! n. o
0 OX
I/I U >-
'MX CC
< C 0
i— a
O Z
1 (3 (-
< 0-
I
1
j
COOMOOOC
f\lf*-*HOih-inr<
H 0 O O O O O C
•
,
Z Ul ' '
a M u aaraaooac
Ul t-
f«- ^' CO U
J * * * 51
Hl/lttl/l HHOOfMCMNfV
,Z 0 ul.
(-
L. I-
•a
'o cr
*•
Z Z
u
o c
'
u c
i *-
I1 w*i I™11 ii iwi' iw
! . . o l t-i N HI * ! i ' . !
! ! i * z, 1
i
1 1
! i < ! ! : ;
l •
1 I ,
i ' 1 ' •
* 1
j , i ' i
! ' ' : : : I
;»
1
I
X
o
H
0. •
Ul O
> u z ,
Ul
a
3
O
x
I
1
1
j ;
I
1
H N fO d
i
j
1
1
|
i
!
1
f « u) r- e:
1
i
78
-------
CO
<
€>
LO
cc
O
O
uu
cc
©
O
Q.
CC
O
a.
LU
O
LU
CC
©
(O
cc
O
a.
LU
O
CC
O
a.
LU
a
o
o
CO
00
cc
O
O
LU
cc-
s
CM
O
O
CNJ
CO
cc
LU
o
to
o
o
o
in
cc
o
a.
LU
O
LU
CC
o
LO
CM
O
O
CM'
O
LO
O
O
O
in
SH3131AI
79
-------
NORTH
WIND
SOURCE 1
SOURCE 2
RECEPTOR
EAST
Figure 2. Upwind and crosswind distances of point sources
from a receptor.
80
-------
NORTH
WIND
AREA
SOURCE
RECEPTOR
-*1 EAST
Figure 3. Minimum, XL and maximum, X2, upwind distances of an area source
from a receptor and crosswind line sources.
81
-------
NORTH
WIND
Figure 4. Line source and receptor relationships.
82
-------
NORTH
WIND
(Rp.Sp)
RECEPTOR
(Rr.Sr)
EAST
Figure 5. Curved source and receptor relationships.
83
-------
PART 2
The slant (or special line) and special curved path source calculations
are made in a similar way as the line source and curved source calculations,
respectively. Unlike Equation A21 , the concentration from special sources is
* • i ft «,
-------
The vehicle acceleration can also be expressed as
2 2
\li - r
where Xf = total distance of travel.
The distance of travel to any point "P" can be expressed as
where X = distance of travel to P.
Solving for t in Equation 5 and substituting in Equation 3,
The emission rate q ( •? ) is now given by
X» S 6C"" m
,t - Tr i (7,
3600
where TV = traffic volume.
The emission rate is inversely proportional to V . In order to ensure that
q will not approach infinity as V goes to zero, a simple technique is used
to set a minimum vehicle speed. The minimum speed (V ) is calculated using
the traffic volume and a gross estimate of average vehicle length (VL).
.« in\ vtii \ O /
S 3600
V is then in (~-)> and physically it is the slowest speed the vehicles
could be going and still maintain the traffic volume. If V is less than Vs,
then V is set equal to Vg.
85
-------
For most applications this change in vehicle speed will have negligible
effect on the concentration estimates.
86
-------
APPENDIX B
FORTRAN STATEMENTS
87
-------
PAL
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
10
C
PAL VERSION 78001
PAL CALCULATES HOURLY CONCENTRATIONS FROM POINT,AREA, AND
LINE SOURCES.
PAL READS ALL INPUT DATA CARDS.SOURCE INPUT DATA CARDS
INCLUDE POINT,AREA,HORIZONAL LINE,SPECIAL LINE.HORIZONAL
CURVED PATH,AND SPECIAL CURVED PATH DATA CARDS.
ALL SOURCE INPUT DATA CARDS ARE OPTIONAL.
METEOROLOGY AND RECEPTOR INPUT DATA CARDS ARE
NOT OPTIONAL.
THE FOLLOWING SUBROUTINES ARE CALLED BY PAL.
POINT—POINT SOURCE SUB-MODEL.
AREA--AREA SOURCE SUB-MODEL.
HRZLN--HORIZONAL LINE SOURCE SUB-MODEL.
CRVLN--CURVED PATH SOURCE SUB-MODEL.
SPCLN--SPECIAL LINE SOURCE SUB-MODEL.
SPCCR--SPECIAL CURVED PATH SOURCE SUB-MODEL.
COMMON /SOP/ QP(31),HPP(31),TSP(31),VSP(3D,DP(31) ,VFP(3D,RQP(31)
* \ niri"vTi/^\^> e~irrr\ri\/'^'*\ *~i /~\ JIT T-\ / o ^ o. 11 \ TNTTTI/^^I— \ ~r TI rr r\
(3D,CO.NA(31
V^wliiiv/ll /k^^J./ *at ••• \.J I /}****• \«J'/JAfc-'A \ .J • / J * fc-/* V _J ' / ) •L/ •*• V ,J I / ) » *• A
1 ,SQP(3D ,SYOP(3D ,SZOP(31) ,CONP(31,24) ,DVP(25) ,IUZP
COMMON /SOA/ QA(3D,HQ(31),RQ(31),SQ(31),DEST(3D ,DNOR
,24),DVA(25),IUZA
COMMON /SOL/ QLN(3D,HLN(31),RAQ(3D ,SAQ(31),RBQ(31),SBQ(3D,SYO(3
1 1),SZO(3D,CONLH(31,24),DVH(25),IUZH
COMMON /SOCP/ QLNS(3D,HLNS(3D ,RBQS(3D ,SBQS(3D , RMQS( 3 1) , SMQS( 3 1
1) ,REQS(3D ,SEQS(3D ,SIYO(3D ,SIZO(3D , CONCLN( 3 1 , 24) , DDVH( 31 ) , IUZC ,
2RADIUS(3D ,NLANE(3D
COMMON /SOLS/ QLS( 31) , HAS( 31) ,HBS(3D ,RAS(3D ,SAS(3D , RBS( 3-1 ) ,SBS(
131) ,SYOS(3D ,SZOS(3D ,CONLS(31 ,24) ,DVS(25) , IUZS,SPDI( 3 1 ) ,SPDF(3D ,
2TVSL(3D ,VSSL(3D
COMMON /SOCS/ QLNA(3D ,HCL(3D ,RBQA(3D ,SBQA( 3 1) , RMQA( 3 1) ,SMQA(3D
1 ,REQA(3D ,SEQA(3D ,SIYA(3D ,SIZA(3D ,CONCLA(31 ,24) ,DVHA(3D ,IUZE,S
2PEK3D ,SPEF(3D ,TVCL(3D,VSCL(3D
COMMON /EEC/ RR(31),SR(31),ZR(31)
COMMON /WEA/ WTHET(25),WU(25),MKST(25),WHL(25),WTA(25),UHGT
NO','YES'/
, ACAC3D,
. . QLTS(3D,
.IVER/78001/
DIMENSION ALP(20), CONT(30,24), ACP(31)
11), ACT(3D, QLT(20), YAN(3), ACLC(31)
DATA YAN /«
IRD = 5
IWRI=6
MAXQ=31
MAXR=31
READ (IRD,790) ALP
WRITE (IWRI,800) ALP
ALP IS 80 COLS OF
ACLH(3D,
ACAC(3D
ACLS(3
,IVER
.__ _. ALPHANUMERIC INFORMATION TO HELP IDENTITY OUT
READ (IRD,810) PINA,PINL,IQP,IUZP,IQA,IUZA,IQLH,IUZH,IQLC,IUZC,IQL
1S,IUZS,IQAC,IUZE,IAVG,IDRNL,UHGT
C
C
C
C
C
C
C
C
PINA
PINL
IQP
IUZP
IQA
IUZA
IQLH
IUZH
X.
X.
X
X
X
X
X
X
XXXX
XXXX
TEST FOR AREA INTEGRATION ACCURACY.
TEST FOR LINE INTEGRATION ACCURACY,
CONTROL FOR POINT SOURCES.
WIND INCREASE WITH HEIGHT,POINT
CONTROL FOR AREA SOURCES.
WIND INCREASE WITH HEIGHT,AREA.
CONTROL FOR HORIZ. LINE SOURCES.
WIND INCREASE WITH HOT.,HORIZ.LINE,
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
0010
0020
0030
0040
0050
0060
0070
0080
0090
0100
0110
0120
0130
0140
0150
0160
0170
0180
0190
0200
0210
0220
0230
0240
0250
0260
0270
0280
0290
0300
0310
0320
0330
0340
0350
0360
0370
0380
0390
0400
0410
0420
0430
0440
0450
0460
0470
0480
0490
0500
0510
0520
0530
88
-------
PAL
c
c
c
c
c
c
c
c
c
c
IQLC
IUZC
IQLS
IUZS
IQAC
IUZE
IAVG
IDRNL
UHGT
CO!
X
X
X
X
X
X
X
X
x:
JT:
c
20
30
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
40
50
60
C
70
80
C
CONTROL FOR CURVED PATH SOURCES.
WIND INCREASE WITH HGT.,CURVED PATH.
CONTROL FOR SLANT LINE SOURCES.
WIND INCREASE WITH HGT.,SLANT LINE.
CONTROL FOR SPECIAL PATH SOURCES.
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
WIND INCREASE WITH HGT. ,SPECIAL PATH. (DIMENSIONLESS '.
(DIMENSIONLESS)
(DIMENSIONLESS)
(M)
CONTROL FOR AVG.CONCENTRATIONS.
CONTROL FOR DIURNAL FACTORS.
XXX.X HEIGHT APPLICABLE TO WIND SPEED.
VALUES: IQ ARE 1 FOR NON-INCLUSION, 2 FOR INCLUSION.
WRITE (IWRI,820) PINA,PINL,YAN(IQP+1),YAN(IUZP+1),YAN(IQA+1),YAN(I
1UZA+1),YAN(IQLH+1),YAN(IUZH+1),YAN(IQLC+1),YAN(IUZC+1),YAN(IQLS+1)
2,YAN(IUZS+1),YAN(IQAC+1),YAN(IUZE+1),YAN(IAVG+1),YAN(IDRNL+1),UHGT
GO TO (bO,20), IQP
WRITE HEADING FOR POINT SOURCE INFORMATION.
WRITE (IWRI,830)
1=1
READ (IRD,810) ICARD,QP(I),HPP(I),TSP(I),VSP(I),DP(I),VFP(I),RQP(I
1) ,SQP(I) ,SYOP(I),SZOP(I)
DECIMALS MUST BE INCLUDED IN REAL VARIABLES.
A COMMA MUST BE PLACED BETWEEN EACH VARIABLE.
IF A VARIABLE IS NOT USED SET IT = TO 0.
ICARD=1 ON ALL SOURCE INPUT CARDS EXCEPT THE LAST THEN ICARD=2.
*
I CARD
QP
HPP
TSP
VSP
DP
VFP
RQP
SQP
SYOP
SZOP
* * P 0 I
X
XXXXX.XX
XXX. X
XXX. X
XX. X
XX. XX
xxxx.x
xxx.xxxx
xxx.xxxx
x.x
x.x
NT SOURCE CARD***
EQUALS 1 OR 2 SEE ABOVE.
POINT SOURCE STRENGTH.
PHYSICAL STACK HEIGHT.
STACK GAS TEMPERATURE.
STACK GAS VELOCITY.
STACK INSIDE DIAMETER.
STACK GAS VOLUME FLOW.
EAST COORDINATE OF STACK.
NORTH COORDINATE OF STACK.
INITIAL SIGMA Y.
INITIAL SIGMA Z.
(DIMENSIONLESS;
(G/SEC)
(M)
(DEC K)
(M/SEC)
(M)
(M**3/SEC)
(KM)
(KM)
(M)
(M)
IF (I-MAXQ) 40,40,50
WRITE (IWRI,840) I,QP(I),HPP(I),TSP(I),VSP(I),DP(I),VFP(I),RQP(I),
1SQP(I),SYOP(I),SZOP(I)
IF (ICARD.GT.l) GO TO 60
1 = 1 + 1
GO TO 30
WRITE (IWRI,850) MAXQ
CALL EXIT
GO TO (110,70), IQA
tvRITE HEADING FOR AREA SOURCE INFORMATION.
WRITE (!WRI,b60)
J=l
READ (IRD,810) ICARD,QA(J),HQ(J) ,RQ(J) ,SQ(J) ,DEST(J),DNOR(J)
* * * A R E A SOURCE CARD***
ICARD X EQUALS 1 OR 2 SEE ABOVE. (DIflKNSIONLESS)
QA X.XXXXX AREA SOURCE STRENGTH. (G/SEC-M**2)
HQ XX.XX AREA SOURCE HEIGHT. (M)
0540
0550
0560
0570
0580
0590
0600
0610
0620
0&30
0640
0650
06bO
0670
0680
0690
0700
0710
0720
0730
0740
0750
0760
0770
0780
0790
0800
0810
0820
0830
Ob40
0850
0860
0870
0880
0890
0900
0910
0^30
0940
0^50
Oy60
Oy70
1000
1010
1020
1030
1040
1050
1060
89
-------
PAL
C
C
C
C
90
100
110
C
120
C
130
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
140
150
C
C
C
C
160
RQ XXX.XXX
SQ XXX.XXX
DEST XX.XX
DNOR XX.XX
EAST COORDINATE OF S.W. CORNER.
NORTH COORDINATE OF S.W. CORNER.
EAST-WEST SIZE.
NORTH-SOUTH SIZE.
IF (J-MAXQ) 90,yO,100
WRITE (IWRI,870) J,QA(J),HQ(J),RQ(J) ,SQ(J) ,DEST(J),DNOR(J)
IF (ICARD.GT.l) GO TO 110
J=J + 1
GO TO 80
WRITE (IWRI,880) MAXQ
CALL EXIT
GO TO (210,120) , IQLH
WRITE HEADING FOR HORIZONTAL LINE SOURCE INFORMATION.
WRITE (IWRI,890)
READ A SOURCE CARD.
K=l
READ (IRD,810) ICARD,HLN(K),NL,RAQ(K),SAQ(K),RBQ(K),SBQ{K),SYO(K),
1SZO(K),WT,WM,(QLT(J),J=1,NL)
THE PROGRAM WILL HANDEL 1 LANE OR AN EVEN NUMBER OF LANES.
AN ODD NUMBER OF LANES MUST BE SEPERATED INTO 2 I
AN EVEN NUMBER PLUS 1 SOURCE.
* *
ICARD
HLN
NL
RAQ
SAQ
RBQ
SBQ
SYO
SZO
WT
WM
QLT(l)
QLT(2)
•
•
•
QLT(NL;
* H 0 R I Z 0
CARD
X
XX. X
XX.
xxx.xxxx
xxx.xxxx
xxx.xxxx
xxx.xxxx
x.x
x.x
XX. X
XX. X
x.xxxxx
x.xxxxx
x.xxxxx
x.xxxxx
x.xxxxx
1 X.XXXXX
NTAL LINE SOURl
1
EQUALS 1 OR 2 SEE ABOVE.
LINE SOURCE HEIGHT
NO. OF LINES IF MULTI-LINE
EAST COORDINATE, POINT A
NORTH COORDINATE, POINT A
EAST COORDINATE, POINT B
NORTH COORDINATE, POINT B
INITIAL SIGMA Y
INITIAL SIGMA Z
TOTAL WIDTH
WIDTH OF MEDIAN
LINE SOURCE STRENGTH
ONE FOR EACH OF ANL SOURCES
ORDERING IS FIRST SOURCE
STRENGTH IS FOR LEFT-MOST
LINE WHEN LOOKING FROM
POINT A TO POINT B
IF (K-MAXQ) 140,140,200
IF (NL-1.) 190,190,150
ANL=NL
IF (K+NL-1-MAXQ) 160,160,200
WT IS THE TOTAL WIDTH OF THE LINE SOURCE INCLUDING THE MEDIAN.
WM IS THE WIDTH OF THE MEDIAN (METERS)
LANES ARE ORDERED FROM LEFT TO RIGHT WHEN LOOKING FROM POINT
TO POINT B.
WRITE (IWRI,900) ANL,WT,WM
WL=(WT-WM)/ANL
RA=RAQ(K)
SA=SAQ(K)
KM)
KM)
KM)
KM)
DNOR(J)
iTION.
,SBQ{K),SYO(K) ,
' LANES.
IRCES
E * * *
(DIMENSIONLESS)
(M)
(DIMENSIONLESS)
(KM)
(KM)
(KM)
(KM)
(M)
(M)
(M)
(M)
(G/SEC-M)
> >
i t
i >
i t
i t
:NG THE MEDIAN.
1G FROM POINT A
1070
1080
1090
1100
1110
1120
1130
1140
1150
1160
1170
1180
1190
1200
1210
1220
1230
1240
1250
1260
1270
1280
1290
1300
1310
1320
1330
1340
1350
1360
1370
1380
1390
1400
1410
1420
1430
1440
1450
1460
1470
1480
1490
1500
1510
1520
1530
1540
1550
1560
1570
1580
1590
90
-------
PAL
RB=RBQ(K) 1600
SB=5BQ(K) 1610
SYON=3YO(K) 1620
SZON=SZO(K) 1630
HLNN=HLN(K) 1640
DELR=RB-RA 1650
DELS=SB-SA 1660
DIST=SQRT(DELS*DELS+DELR*DELR) 1670
NLIM=NL/2 1660
ALIM=NLIM 1690
DO 170 KN=1,NLIM 1700
A=KN 1710
DL=(-0.5)*WM+((-1)*ALIM-0.5+A)*WL 1720
DUM=DL*0.001/DIST 1730
ID=K+KN-1 1740
RAQ(ID)=RA+DELS*DUM 1750
RBQ(ID)=RB+DELS*DUM 1760
SAQ(ID)=SA-DELR*DUM 1770
SBQ(ID)=SB-DELR*DUM 1780
SYO(ID)=SYON 1790
SZO(ID)=SZON 1800
QLN(ID)=2LT(KN) 1810
HLN(ID)=HLNN 1820
170 WRITE (IWRI,910) ID,QLN{ID),HLN(ID),RAQ(ID),SAQ(ID),RBQ(ID),SBQ(ID 1830
1),SYO(ID),SZO(ID) 1840
NS=NLIM+1 1850
AS=NS 1860
DO IdO KN=NS,NL 1870
A=KN 1880
DL=0.5*WM+(0.b+A-AS)*WL 1890
DUM=DL*0.001/DIST 1900
ID=K+KN-1 1910
RAQ(ID)=RA+DELS*DUM 1920
RSQ(ID)=RB+DELS*DUM 1930
SAQ(ID)=SA-DELR*DUM 1940
SBQ(ID)=SB-DELR*DUM 1950
SYO(ID)=3YON 1960
SZO(ID)=SZON 1970
QLN(ID)=QLT(KN) 1980
HLN(ID)=HLNN 1990
180 WRITE (IwRI,910) ID,QLN(ID),HLN(ID),RAQ(ID),SAQ(ID),RBQ(ID),SBQ(ID 2000
1),SYO(ID),SZO(ID) 2010
K=K+NL 2020
IP (ICARD.GT.l) GO TO 210 2030
GO TO 130 2040
C WRITE LINE SOURCE INFORMATION. 2050
190 QLN(K)=QLT(1) 2060
WRITE (IWRI,910) K,QLN(K),HLN(K),RAQ(K),SAQ(K),RBQ(K),SBQ(K),SYO(K 2070
1),SZO(K) 2080
IF (ICARD.GT.l) GO TO 210 2090
K=K+1 2100
GO TO 130 2110
200 WRITE (IWRI,920) MAXQ 2120
91
-------
PAL
210
C
220
C
230
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
240
CALL EXIT
CONTINUE
GO TO (290,220) , IQLC
WRITE HEADING FOR CURVED HORIZONAL LINE SOURCE INFORMATION
WRITE (IWRI,930)
READ A SOURCE CARD
KK=0
LL=1
READ (IRD,810) ICARD,HLNS(LL),NL,RBQS(LL),SBQS(LL),RMQS(LL),SMQS(L
1L),REQS(LL),SEQS(LL),SIYO(LL),SIZO(LL),WT,WM,(QLTS(JJ),JJ=1,NL)
THE PROGRAM VvILL HANDEL 1 LANE OR AN EVEN NUMBER OF LANES.
AN ODD NUMBER OF LANES MUST BE SEPERATED INTO 2 SOURCES
AN EVEN NUMBER PLUS 1 SOURCE.
A COMMA MUST BE PLACED BETWEEN EACH VARIABLE
*** CURVED PATH SOURCES ***
EQUALS 1 OR 2 SEE ABOVE.
PATH SOURCE HEIGHT.
NUMBER OF LANES.
EAST COORDINATE,PT.A
NORTH COORDINATE,PT.A.
EAST COORDINATE,PT.B.
NORTH COORDINATE,PT.B.
EAST COORDINATE,PT.C.
NORTH COORDINATE,PT.C.
INITIAL SIGMA Y.
INITIAL SIGMA Z.
TOTAL WIDTH OF PATH.
WIDTH OF MIDIAN.
PATH SOURCE STRENGTH.
ONE FOR EACH LANE.
ORDERING OF EMISSION RATES
ARE FROM OUTSIDE LANES TO
INSIDE.
ANL=NL
NLANE(LL)=NL
A=l.
ALIM=ANL/2.
WL=(WT-WM)/ANL
DO 270 KL=1,NL
KK=KK+1
IF (KK.GT.31) GO TO 280
QLNS(KK)=QLTS(KL)
RADIUS IS A CORRECTION FACTOR TO ADJUST THE RADIUS
FOR DIFFERENT LANES.
RADIUS(KK)=(A*.5*WM+ALIM*WL-.5*WL*A)*.001
ALIM=ALIM-1.
IF (ALIM.LT.O.l.AND.ALIM.GT.-.l) ALIM=-1.
IF (ALIM.LT.O.) A=-l.
IF (KL.EQ.l) GO TO 240
GO TO 250
WRITE (IWRI,940) KK,QLNS(KK),HLNS(LL),RBQS(LL),SBQS(LL),RMQS(LL),S
IMQS(LL),REQS(LL),SEQS(LL),SIYO(LL),SIZO(LL),NL,WT,WM
ICARD X
HLNS XX.X
XX
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX.XXXX
x.x
x.x
XX.X
XX. X
x.xxxxx
NL
RBQS
SBQS
RMQS
SMQS
REQS
SEQS
SIYO
SIZO
WT
QLTS
(DIMENSIONLESS)
(M)
(DIMENSIONLESS)
(KM)
(KM)
(KM)
(KM)
(KM)
(KM)
(M)
(M)
(M)
(M)
(G/SEC-M)
2130
2140
2150
2160
2170
2180
2190
2200
2210
2220
2230
2240
2250
2260
2270
2280
2290
2300
2310
2320
2330
2340
2350
2360
2370
2360
2390
2400
2410
2420
2430
2440
2450
2460
2470
2480
2490
2500
2510
2520
2530
2540
2550
2560
2570
2580
2590
2600
2610
2620
2630
2640
2650
92
-------
PAL
250
260
270
280
290
C
300
C
310
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
320
330
340
C
350
C
360
C
C
GO TO 270
IF (KL.GT.2) GO TO 260
WRITE (IWRI,770)
WRITE (IWRI,780) KK,QLNS(KK)
CONTINUE
IF (ICARD.GT.l) GO TO 290
LL=LL+1
GO TO 230
WRITE (IWRI,950) MAXQ
GO TO (340,300), IQLS
WRITE HEADING FOR SLANT OR VERTICAL LINE SOURCE INFORMATION.
WRITE (IWRI,960)
READ A SOURCE CARD.
L=l
READ (IRD,810) ICARD,QLS(L),HAS(L),HBS(L),RAS(L),SAS(L),RBS(L),SBS
(L),SPDI(L),SPDF(L),SYOS(L),SZOS(L),TVSL(L),VSSL(L)
A COMMA MUST BE PLACED BETWEEN EACH VARIABLE
***SLANT OR VERTICAL LI
ICARD X EQUALS 1 OR 2 SEE ABOVE.
QLS X.XXXXX LINE SOURCE STRENGTH.
HAS XX.X HEIGHT OF POINT A.
HBS XX.X HEIGHT OF POINT B.
RAS XXX.XXXX EAST COORDINATE,POINT A.
SAS XXX.XXXX NORTH COORDINATE,POINT A.
RBS XXX.XXXX EAST COORDINATE,POINT B.
SBS XXX.XXXX NORTH COORDINATE,POINT B.
SPDI XXX.X SPEED AT POINT A.
SPDF XXX.X SPEED AT POINT B.
SYOS X.X INITIAL SIGMA Y.
SZOS X.X INITIAL SIGMA Z.
TVSL XXXXX. VEHICLE VOLUME.
VSSL XX.X GROSS ESTIMATE OF VEH. SIZE.
N E
(DIMENSIONLESS!
(G/SEC)
(M)
(M)
(KM)
(KM)
(KM)
(KM)
(M/SEC)
(M/SEC)
(M)
(M)
(VEH/HR)
(M)
IF (L-MAXQ) 320,320,330
WRITE (IWRI,970) L,QLS(L),HAS(L),HBS(L),RAS(L),SAS(L),RBS(L),SBS(L
1) ,SPDI(L),SPDF(L),SYOS(L),SZOS(L),TVSL(L),VSSL(L)
IF (ICARD.GT.l) GO TO 340
L=L+1
GO TO 310
WRITE (IWRI,980) MAXQ
CALL EXIT
GO TO (390,350) , IQAC
WRITE HEADING FOR SPECIAL PATH SOURCE INFORMATION
WRITE (IWRI,990)
READ A SOURCE CARD
LC=1
READ (IRD,dlO) ICARD,QLNA(LC),HCL(LC),RBQA(LC),SBQA(LC),RMQA(LC),S
IMQA(LC),REQA(LC),SEQA(LC),SPEI(LC),SPEF(LC),SIYA(LC),SIZA(LC),TVCL
2(LC),VSCL(LC)
IF (LC.GT.30) GO TO 370
A COMMA MUST BE PLACED BETWEEN EACH VARIABLE
**SPECIAL PATH SOURCES ***
2660
2670
2680
2690
2700
2710
2720
2730
2740
2750
2760
2770
2780
2790
2800
2810
2820
2830
2840
2850
2860
2870
2880
2890
2900
2910
2920
2930
2940
2950
2960
2970
2980
2990
3000
3010
3020
3030
3040
3050
3060
3070
3080
3090
3100
3110
3120
3130
3140
3150
3160
3170
3180
93
-------
PAL
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
ICARD
QLNA
HCL
RBQA
SBQA
RMQA
SMQA
REQA
SEQA
SPEI
SPEF
SIYA
SIZA
TVCL
VSCL
X
X.XXXXX
XX. X
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX.XXXX
XXX. X
XXX. X
x.x
x.x
xxxxx.
XX. X
370
380
C
390
C
400
C
C
C
C
C
C
C
C
410
420
430
C
C
440
450
460
EQUALS 1 OR 2 SEE ABOVE.
PATH SOURCE STRENGTH.
PATH SOURCE HEIGHT.
EAST COORDINATE,POINT A.
NORTH COORDINATE,POINT A.
EAST COORDINATE,POINT 8.
NORTH COORDINATR,POINT 6.
EAST COORDINATE,POINT C.
NORTH COORDINATE,POINT C.
SPEED AT POINT A.
SPEED AT POINT C.
INITIAL SIGMA Y.
INITIAL SIGMA Z.
VEHICLE VOLUME.
GROSS ESTIMATE OF VEH. SIZE.
(DlrtENSIONLESS)
(G/SEC)
(M)
(KM)
(KM)
(KM)
(KM)
(KM)
(KM)
(M/SEC)
(M/SEC)
(M)
(M)
(VEH/HR)
(M)
) LC,QLNA(LC) ,HCL(LC) ,RBQA(LC) ,S6QA(LC) ,RMQA(LC) ,S
,SEQA(LC) ,SPEI(LC) ,SPEF(LC) ,SIYA(LC) ,SIZA(LC) ,TVCL
WRITE
IMQA(LC) ,REQA(LC
2(LC) ,VSCL(LC)
IF (ICARD.GT.l) GO TO 380
LC=LC+1
GO TO 360
WRITE (IWRI,1010) MAXQ
CONTINUE
WRITE HEADING FOR PRINT-OUT OF RECEPTOR INFORMATION.
WRITE (IWRI,1020)
READ A RECEPTOR CARD.
H=l
READ (IRD,810) ICARD , RR(N ) ,SR(N ), ZR (N )
A COMMA MUST BE PLACED BETWEEN EACH VARIABLE
***RECEPTOR CARD***
ICARD X EQUALS 1 OR 2 SEE ABOVE.
RR XXXX.XXX R COORDINATE
SR XXXX.XXX S COORDINATE
ZR XXX. X HEIGHT ABOVE GROUND
A COMMA MUST BE PLACED BETWEEN EACH VARIABLE
IF UN-MAXR) 4io,4iu,420
WRITE RECEPTOR INFORMATION.
WRITE (IWRI,1030) N ,RRU ) ,SR(N ) , ZR(N )
IF (ICARD.GT.l) GO TO 430
N=N+1
GO TO 400
WRITE (IWRI,1040) MAXR
CALL EXIT
M=l
WRITE MET HEADING.
WRITE (IWRI,1050)
READ A METEOROLOGY CARD.
GO TO (450,460), IDRNL
READ (IRD,810) ICARD ,WTHET(M) ,WU (M) ,MKST(M) ,
GO TO 470
READ (IRD,810) ICARD ,WTHET(M) ,WU (M) ,MKST(M) ,WHL( M) ,WTA(M) ,DVP (M) ,D
1VA(M) ,DVH(M) ,DDVH(M) ,DVS(M) ,DVHA(M)
(DIMENSIONLESS )
(KM)
(KM)
(M)
,WTA(M)
3190
3200
3210
3220
3230
3240
3250
3260
3270
32bO
3290
3300
3310
3320
3330
3340
3350
3360
3370
3380
3390
3400
3410
3420
3430
3440
3450
3460
3470
34bO
34*0
3500
3510
3520
3530
3540
3550
3560
3570
3580
3590
3600
3610
3620
3630
3640
3650
3660
3670
3680
3690
3700
3710
94
-------
PAL
(DIMENSIONLESS)
(DEC AZIMUTH)
(M/SEC)
(DIMENSIONLESS)
(M)
(DEC K)
(DIMENSIONLESS)
(DIMENSIONLESS)
(DIMENSIONLESS)
C A COMMA MUST BE PLACED BETWEEN EACH VARIABLE
C ***METEOROLOGY CARD***
C ICARD X EQUALS 1 OR 2 SEE ABOVE.
C WTHET XXX. WIND DIRECTION
C WU XX. X WIND SPEED
C MKST X STABILITY CLASS
C WHL XXXX. MIXING HEIGHT
C WTA XXX. X AMBIENT AIR SURFACE TEMP.
C DVP X.XXX DIURNAL VARIATION, POINT SOURCES (DIMENSIONLESS )
C DVA X.XXX DIURNAL VARIATION, AREA SOURCES (DIMENSIONLESS)
C DVH X.XXX DIUR. VAR. , HORIZ. LINE SOURCES (DIMENSIONLESS)
C DVS X.XXX DIUR. VAR., SLANT LINE SOURCES
C DDVH X.XXX DIUR. VAR., CURVE PATH SOURCES
C DVHA X.XXX DIUR. VAR. , SPECIAL PATH SOURCES
C THE DIURNAL VARIATION VALUES ARE THE RATIO OF THAT EMITTED
C DURING THIS HOUR TO THE INPUT Q. IF THE SAME AS THE INPUT,
C THE VALUE IS 1.000
470 CONTINUE
IF (M-25) 490,4«0,480
480 WRITE (IWRI,1060)
CALL EXIT
490 GO TO (500,510), IDRNL
500 DVP(M)=1.
DVA(M)=1.
DVH(M)=1.
DVS(M)=1.
DDVH(M)=1.
DVHA(M)=1.
C WRITE METEOROLOGY INFORMATION.
510 WRITE (IWRI,1070) M ,WTHET(M) ,WU (M) , MKST (M) ,WHL( M) ,WTA( M) , DVP (M) ,DV
1A ( M ) , DVH ( M ) , DDVH ( M ) , DVS ( M ) , DVHA ( M )
IF (ICARD. GT.l) GO TO 520
M=M+1
GO TO 440
520 CONTINUE
DO 530 IM=1,M
DO 530 IN=1,N
CONP(IN,IM)=0.
CONA(IN,IM)=0.
CONLH(IN,IM)=0.
CONCLN(IN,IM)=0.
CONLS(IN,IM)=0.
CONCLA(IN,IM)=0.
530 CONT(IN,IM)=0.
GO TO (550,540) , IQP
S40 CALL POINT (M,I,N)
550 GO TO (570,560) , IQA
5bO CALL AREA U,J,N,PINA)
570 GO TO (590, 5bO), IQLH
580 CALL HRZLN (M,K,N,PINL)
590 GO TO (olO,600), IQLC
bOO CALL CRVLN (M,LL,N, PINL )
610 GO TO (630,620), IQLS
3720
3730
3740
3750
3760
3770
3780
3790
3800
3810
3820
3830
3d40
3850
3860
3870
3880
3890
3900
3910
3920
3930
3940
3950
3960
3970
3980
3990
4000
4010
4020
4030
4040
4050
4060
4070
4080
4090
4100
4110
4120
4130
4140
4150
4160
4170
4180
4190
4200
4210
4220
4230
4240
95
-------
PAL
620 CALL SPCLN (M,L,N,PINL) 4250
630 GO TO (550,640), IQAC 4260
fc.40 CALL SPCCR (M,LC,N,PINL) 4270
650 CONTINUE 4280
DO 660 IM=1,M 4290
DO 660 IN=1,N 4300
660 CONT(IN,IM)=CONP(IN,IM)+CONA(IN,IM)+CONLH(IN,IM)+CONLS(IN,IM)+CONC 4310
1LN(IN,IM)+CONCLA(IN,IM) 4320
DO 670 IiH=l,w 4330
WRITE (IWRI,1080) 4340
DO b70 IN=1,N 4350
670 WRITE (IWRI,1090) IM,IN,RR(IN),SR(IN),ZR(IN),CONP(IN,IM),CONA(IN,I 4360
1M),CONLH(IN,IM),CONCLN(IN,IM),CONLS{IN,IM),CONCLA(IN,IM),CONT(IN,I 4370
2M) 4380
IF (M-l) 740,740,680 4390
680 GO TO (740,690), IAVG 4400
690 DO 700 IN=1,N 4410
ACP(IN)=0. 4420
ACA(IN)=0. 4430
ACLH(IN)=0. 4440
ACLC(IN)=0. 4450
ACLS(IN)=0. 4460
ACAC(IN)=0. 4470
700 ACT(IN)=0. 4480
DO 710 IM=1,M 4490
DO 710 IN=1,N 4500
ACP(IN)=ACP(IN)+CONP(IN,IM) 4510
ACA(IN)=ACA(IN)+CONA(IN,IM) 4520
ACLH(IN)=ACLH(IN)+CONLH{IN,IM) 4530
ACLC(IN)=ACLC(IN)+CONCLN(IN,IM) 4540
ACLS(IN)=ACLS(IN)+CONLS(IN,IM) 4 550
ACAC(IN)=ACAC(IN)+CONCLA(IN,IM) 4560
710 ACT(IN)=ACT(IN)+CONT(IN,IM) 4570
DO 720 IM=1,N 4580
ACP(IN)=ACP(IN)/M 4590
ACA(IN)=ACA(IN)/M 4600
ACLH(IN)=ACLH(IN)/M 4610
ACLC(IN)=ACLC(IN)/M 4620
ACLS(IN)=ACLS(IN)/M 4630
ACAC(IN)=ACAC(IN)/M 4640
720 ACT(IN)=ACT(IN)/M 4650
WRITE (IWRI,1100) M 4660
NDUt4=0 4670
WRITE (IWRI,1080) 4680
DO 730 IN=1,N 4690
730 WRITE (IrtRI,1090) NDUM,IN,RR(IN),SR(IN),ZR(IN),,ACP(IN),ACA(IN),AC 4700
ILH(IN),ACLC(IN),ACLS(IN),ACAC(IN),ACT(IN) 4710
740 READ (IRD,810,END=760) KTL 4720
IF (KTL) 760,760,750 4730
750 GO TO (10,390,430,760), KTL 4740
760 CALL EXIT 4750
C 4760
770 FORMAT (/,10X,'SOURCE STRENGTH FOR THE REMAINING LANES',/,15X,'LAN 4770
96
-------
PAL
780
790
800
810
820
830
840
850
860
870
880
890
•900
910
920
930
940
950
1E',5X,'SOURCE STRENGTH',/)
FORMAT (17X,I2,T28,F10.8)
FORMAT (20A4)
FORMAT (1H1,20A4,5X,'VERSION ',15)
FORMAT ()
FORMAT (1HO,'PINA =',F10.5,' PINL =',
1D INCREASE',/,13X,'INCLUDED WITH
F10.5,//,14X,'SOURCE',5X,'WIN
HEIGHT',//,4X,'POINT',6X,A3,1
22X,A3,//,4X,'AREA',7X,A3,12X,A3,//,1X,'HORIZONTAL',/,1X,'LINE SOUR
3CE' , 3X,A3,12X,A3,//,4X, '
CURVE' ,/,
1X, 'PATH
SOURCE' , 3X, A3 , 1 2X , A3 , //
4,3X,'SPECIAL',/,1X,'LINE SOURCE',3X,A3,12X,A3,//,3X,'SPECIAL',/,1X
5,'PATH SOURCE',3X,A3,12X,A3,1 OX,'AVERAGE',5X,A3,' DIURNAL',3X,A3
6,5X,'HEIGHT AT WIND SPEED',1X,F6.1,
FORMAT (1HO,'* **POINT SO
10INT PHYSICAL STACK STACK
2DINATES INITIAL SIGMAS'
3 GAS DIAMETER FLOW
1X,'METERS',//)
URGES***',/,7X,'NO. P
STACK VOLUME COOR
,/13X,'SOURCE HEIGHT TEMP
EAST NORTH Y
4',/12X,'STRENGTH (METERS) (DEG-K) VELOCITY (METERS) (CU
(KM)
(M)
Z
M/S
(M)V13X,'(G/SEC)',24X,'(M/SE
THE STORAGE A
AR
5EC) (KM)
6C)',/)
FORMAT (1H ,I7,3X,F10.2,2F10.1,3F10.2,2F10.3,2F10.1)
FORMAT (1HO,'THE NUMBER OF POINT SOURCE CARDS EXCEED
1LLOTTED FOR THEM IN THE DIMENSION STATEMENT:',13)
FORMAT (///,'* **AREA SOURCES** *',/,7X,'NO,
1EA AREA COORDINATES AREA SIZE',/,13X,'SOURCE
2 SOURCE SW-CORNER',/,12X,'STRENGTH HEIGHT EAST
3 NORTH EAST-WEST NORTH-SOUTH',/,10X,'(G/SEC-M**2) (METERS) (KM
4) (KM) (KM) (KM)',/)
(1H ,I7,3X,F10.8,F10.1,4F10.3)
(1HO,'THE NUMBER OF AREA SOURCE CARDS EXCEED THE STORAGE AL
FOR THEM IN THE DIMENSION STATEMENT:',13)
(///,'* **HORIZONTAL LINE SOURCES*
LINE LINE',17X,'COORDINATES',/,13X,'SOURC
FORMAT
FORMAT
1LOTTED
FORMAT
1 * *',/,7X,'NO.
2E SOURCE
3X,'NO. TOTAL
4 NORTH EAST
POINT A'
13X,
,12X,
NORTH',8X,'Y'
MEDIAN',/,
'POINT B',10X,'INITIAL
'STRENGTH HEIGHT
SIGMAS',11
W
5IDTH',/,11X,
(G/SEC-M) (METERS)
EXCEED
13)
ALP
THE STORAGE AL
A T H
'S 0
EAST
9X,'Z',19X,'OF WIDTH
(KM) (KM) (KM)
6 (KM) (M) (M)',17X,'LANES (METERS) (METERS) ')
FORMAT (1H ,100X,F10.0,2F10.1)
FORMAT (1H ,I7,3X,F10.8,F10.1,4F10.3,2F10.2)
FORMAT (1HO,'THE NUMBER OF LINE SOURCE CARDS
1LOTTED FOR THEM IN THE DIMENSION STATEMENT:',
FORMAT (///,'* **CURVED HORIZON
1URCES *** ',/,2X,'NO.',T10,'PATH',T21
2ES',/,T9,'SOURCE',T20,'SOURCE',T32,'POINT A'
3INT C',T87,'INITIAL SIGMAS',T108,'NO.',4X,'TOTAL MEDIAN',/,T8,
4'STRENGTH',T20,'HEIGHT',T29,'EAST',T38,'NORTH',T50,'EAST',T59,'NOR
5TH',T68,'EAST',T77,'NORTH',T89,'Y',T99,'Z',T108,'OF WIDTH
6WIDTH1,/,T7,'(G/SEC-M)',T19,'(METERS)',T29,'(KM)',T39,'(KM)',T50,'
7(KM)',T60,'(KM)',T68,'(KM)',T78,'(KM)',T88,'(M)',T98,'(M)',T108,'L
8ANES (METERS) (METERS)')
FORMAT (/,2X,I2,T7,F10.8,T17,F9.3,T26,F9.3,T36,F9.3,T47,F9-3,T57,F
19.3,T66,F9.3,T76,F9.3,T86,F7.2,T96,F7.2,T108,I2,T115,F7.2,2X,F7.2)
FORMAT (///,3X,'THE NUMBER OF CURVED PATH SOURCES EXCEEDS',' THE N
'PATH',T42,'COORDINAT
B',T71,'PO
4780
4790
4800
4810
4820
4830
4840
4850
4860
4870
4880
4890
4900
4910
4920
4930
4940
4950
4960
4970
4980
4990
5000
5010
5020
5030
5040
5050
5060
5070
5080
5090
5100
5110
5120
5130
5140
5150
5160
5170
5180
5190
5200
5210
5220
5230
5240
5250
5260
5270
5280
5290
5300
97
-------
PAL
1UMBER ALLOTTED FOR THEM IN THE DIMENSION STATEMENT:',13) 5310
960 FORMAT (///,'* **SLANT OR VERTICAL LINE S 5320
10URCES** *',/,7X,'NO. SOURCE',7X,'SOURCE HEIGHT',10X,'PO 5330
2INT A1,13X,'POINT B',9X,'INITIAL FINAL',3X,'INITIAL SIGMAS',2X,'T 5340
3RAFFIC1,1X,'VEH.',/,12X,'STRENGTH',8X,'(METERS)',4X,2(5X,'EAST 5350
4 NORTH'),5X,'SPEED1,4X,'SPEED',7X,'Y',6X,'Z',4X,'VOLUME1,2X,'SIZE 5360
51,/,12X,1 (G/SEC) POINT A',3X,'POINT B',4X,4('(KM)',6X),2('(M/S 5370
6EC)',2X),2X,'(M)',4X,'(M)',2X,'(VEH/HR)',1X,'(M)',/) 5380
970 FORMAT (1H ,17,3X,F10.2,2F10.1,6F10.3,F8.2,F7.2,2X,F5.0,IX,F5.1) 5390
980 FORMAT (1HO,'THE NUMBER OF SLANT LINE SOURCE CARDS EXCEED THE STOR 5400
1AGE ALLOTTED FOR THEM IN THE DIMENSION STATEMENT:',13) 5410
990 FORMAT (///,'* **SPECIAL PATH SOURCES ***', 5420
i/, ix, •NO. ', ix, •SOURCE' ,4x, •SOURCE ' ,9x, •POINT A' ,nx, 'POINT B',IIX 5430
2,fPOINT C1,7X,•INITIAL1 ,3X,•FINAL1,4X,'INITIAL SIGMAS' ,2X,'TRAFFIC 5440
3 VEH.',/,4X,'STRENGTH',3X,'HEIGHT',7X,'EAST1,4X,'NORTH1,5X,'EAST1, 5450
44X,1NORTH',5X,'EAST',4X,'NORTH1,4X,'SPEED1,5X,'SPEED1,7X,'Y',6X,'Z 5460
51 ,5X,'VOLUME SIZE',/,5X,' (G/SEC) ' ,2X,( (METERS) ' ,6X,' (KM)' ,5X, ' (KM 5470
6) ' ,5X,'(KM)',5X,' (KM)',5X,' (KM)',5X,' (KM)',4X,'(M/SEC)',2X,'(M/SEC 5480
7)',5X,'(M)',4X,'(M)',3X,'(VEH/HR)',1X,'(M)',/) 5490
1000 FORMAT (IX,12,IX,F8.2,2X,F8.1,2X,6(F8.3,IX),3X,2(F7.3,2X),1X,2F7.3 5500
1,2X,F5.0,1X,F5.1) 5510
1010 FORMAT (///,' THE NUMBER OF SPECIAL PATH SOURCES EXCEEDS THE1,1 ST 5520
10RAGE ALLOTTED FOR THEM IN THE DIMENSION STATEMENT:',13) 5530
1020 FORMAT (///,'* **RECEPTORS** *',/,7X,'NO. RREC(KM) 5540
1 SREC(KM) Z (M)'/) 5550
1030 FORMAT (1H , 110,2F10.3,FlO.1) 5560
1040 FORMAT (1HO,'THE NUMBER OF RECEPTOR CARDS EXCEED THE STORAGE1,' AL 5570
1LOTTED FOR THEM IN THE DIMENSION STATEMENT: ', 13 ) 5580
1050 FORMAT (///,'* **METEOROLOGY** *1,/,7X,1NO. THETA1, 5590
I1 (DEC) U (M/SEC) KST HL (M) T (DEG-K) DIURNAL V 5600
2ARIATIONS (FRACTIONS OF GIVEN Q)',/,T90,'HORIZONAL',2X,'CURVED',3X 5610
3, 'SPECIAL1 ,2X,'SPECIAL' ,/,T74,'POINT',5X,'AREA1 ,5X,1LINE1 ,5X,'PATH 5620
4',5X,'LINE1,5X,'PATH1) 5630
lObO FORMAT (1HO,'THE NUMBER OF MET CARDS EXCEED 24.') 5640
1070 FORMAT (1H ,19 , ' . ' ,FlO.0,FlO.1,110,2F10.0,10X,6(F8.4,IX)) 5650
1080 FORMAT (///,'* **CONCENTRATIONS AT RECEP 5660
IT 0 R S * * *',/,5X,'CONCENTRATIONS IN GRAMS PER CUBIC METER1,T63, 5670
21FROM',6X,'FROM',6X,'FROM',6X,'FROM',/,lX,'HOUR',1X,'RECEPTOR ',3X 5680
3,'RECEPTOR COORDINATES',4X,'FROM1,6X,'FROM',6X,'HORIZONAL',1X,'CUR 5690
4VED',3X,'SPECIAL',3X,'SPECIAL1,/,8X,'NO.',5X,'EAST',4X,'NORTH',3X, 5700
5'HEIGHT1,3X,'POINTS',bX,'AREAS',5X,'LINES',5X,'PATHS',5X,'LINES1,5 5710
6X,'PATHS',5X,'TOTAL',//) 5720
1090 FORMAT (1H ,2X,12,T9,12,T14,2(F8.3,IX),F6.2,T43,1P7(E8.3,2X)) 5730
1100- FORMAT (1H1,'AVERAGE CONCENTRATIONS FOR1,14,' HOURS.') 5740
C 5750
END 5760
98
-------
POINT
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
10'
20
30
40
SUBROUTINE POINT (M,I,N)
POINT CALCULATES HOURLY CONTRATIONS FROM POINT SOURCES,
THE FOLLOWING SUBROUTINES ARE CALLED BY POINT.
FPLUME--HEIGHT OF FINAL PLUME HISE(BRIGGS).
A DISTANCE X(BRIGGS)
CONCENTRATION AT
AT A GIVEN UPWIND
XPLUME--EFFECTIVE PLUME RISE AT
RCONCP—DETERMINES THE RELATIVE
A RECEPTOR FROM A POINT
AND CROSSWIND DISTANCE.
THE FOLLOWING FUNCTIONS ARE CALLED BY POINT.
XVY—CALCULATES THE VIRTUAL DISTANCE NECESSARY
FOR THE INITIAL CROSSWIND DISPERTION.
THE VIRTUAL DISTANCE NECESSARY
FOR THE INITIAL VERTICAL DISPERTION.
COMMON /SOP/ QP(3D ,HPP(3D ,TSP(31) ,VSP(31) ,DP(3D ,VFP(3D ,RQP(3D
,SQP(3O ,SYOP(3D ,SZOP(3D,CONP(31 ,24) ,DVP(25) , IUZP
COMMON /EEC/ RR(31),SR(31),ZR(31)
COMMON /WEA/ WTHET(25),WU(25),MKST(25),WHL(25),WTA(25),UHGT
DIMENSION PUZ(6)
DATA PUZ /O.15,0.15,0.20,0.25,0.40,0.60/
TO ACCOUNT
XVZ—CALCULATES
TO ACCOUNT
STACK
STACK
STACK
EAST
QP
HPP
TSP
VSP
DP
VFP
RQP
SQP
SYOP
SZOP
IWRI=6
DO 200 NW=1 ,M
DO FOR EACH HOUR.
THETA=WTHET(NW)
TR=THETA/57.2958
SINT=SIN(TR)
COST=COS(TR)
U = WU(Nwf)
UZ = U
KST=MKST(NW)
DTHDZ=0.
HL=WHL(NW)
T=WTA(NW)
PWR=PUZ(KST)
CALCULATE FOR ALL
DO 190 N3=1,1
K£H=1
IF (SYOP(NS))
XVYP=0.
GO TO 30
SYON=SYOP(NS)
XVYP=XVY(SYON,KST)
IF (SZOP(NS)) 40,40,50
XVZP=0.
POTNT SOURCE STRENGTH.
PHYSICAL STACK HEIGHT.
STACK GAS TEMPERATURE.
GAS VELOCITY.
INSIDE DIAMETER.
GAS VOLUME FLOW.
COORDINATE OF STACK
NORTH COORDINATE OF STACK
INITIAL SIGMA Y.
INITIAL SIGMA Z.
(G/SEC)
(M)
(DEG K)
(M/3EC)
(M)
(M**3/SEC)
(KM)
(KM)
(M)
(M)
SOURCES,
10,
10,20
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00180
00190
00200
00210
00220
00230
00240
00250
00260
00270
00230
00290
00300
00310
00320
00330
00340
00350
00360
00370
00380
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
99
-------
POI JT
uu
C
C
7U
C
oO
100
no
120
130
140
loO
C
1/0
GO TO oO
3i,ON=S£OP I.MS;
AV'£P=AV2 (S'iDN , v\Sr )
CONTINUE
CALCULaT'J FOR ALL kF.CCP rOtv3 .
IF (X) lou,loU,"/0
DETERi'tl.lC CKOSGfl'JD OI3T,AiiCC
Y=o"SI.MT-K*COST
GO TO (30,l4U), KSil
E3TI laTf] PLtJ.^F '0
CALL XPLliun
,UZ ,
ESTInAT1;: UCLATI\/K CONCBN'f RATION (CHI*[j/j).
CALL RCONCP (iK(-iC) ,H, flL,XZ,Xy,Y,KST,AN,.'iaUi''.,3Y,Si,RC)
liO
200 CONTINUE
RETURN
C
210 PORi'lAT (1HO,'WO COhPGTftTIONS ARE rtTTKi'.PTBO AG SOUSCS ' , 13 , ' HAS A CT
1ACK GAS TEi-lPERATURE LESS T4AN THK AilSI-JNT AIR TGMPE MATURE . ' )
C
END
COboU
J C 00 U
UOolJ
UUu2U
UOb4U
OOooO
OOo7U
OUboO
00 '/JU
00710
00720
00730
00/40
OU/jG
007oO
U07 70
007 !iO
00 /yJ
OOuOO
00o20
00o40
OUODU
OOdoU
Ooou
uuy20
00lj30
00^40
OO^oO
00»70
01000
100
-------
AREA
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
SUBROUTINE AREA ( NHR , NQ , NR , PIN )
AREA CALCULATES HOURLY CONCENTRATIONS FROM
AREA SOURCES.
THE FOLLOWING SUBROUTINES ARE CALLED BY AREA.
PGSIG--CALCULATES SIGMA if, AND SIGMA Z.FROM
THE STABILITY AND DOWNWIND DISTANCE.
RCONCA — CALCULATES CHI*U/Q .RELATIVE CONCENTRATION
NORMALIZED FOR WIND SPEED FOR A RECEPTOR
DOWNWIND OF A CROSSWI.'JD INFINITE LINE SOURCE.
COMMON /SOA/ QA(31 ) ,HQ(3D ,RQ(3D ,SQ(3D ,DEST(3D ,DNOR(3D ,COrJA(31
1 ,24) ,DVA(25) ,IUZA
COMMON /REC/ RR( 3 1 ) ,SR( 3 1 ) , ZR(3 1 )
COMMON /WEA/ WTHET(25) ,WU(25) ,MKST(25) , WHL(25) ,VJTA(25) ,UHGT
DIMENSION GA(78), RA(4), SA(4), XL(4), PUZ(6)
GA ARE VALUES OF THE CUMULATIVE NORMAL DISTRIBUTION IN INCREMEN
OF 0.1 S.
DATA GA /.O, .0001 , .000 1 , .0002 , .0002, .0003, .0005, .0007,. 0010,. 0013,
1 .0019, .0026, .0035, .0047, .0062, .0092, .0107, -0139, .0179, .0228, .0287,
2.0359, .0446, .0548, .0668, .0808, .0968, . 1 151 , . 1 357, .1587, . 1841 , .21 19 ,
3.2420, .2743, .3085, .3446, .3821, .4207, .4602, .5000, .5398 , .5793 , .6 179 ,
4.6554, .6915, .7257, .7580, .7881,. 8159,. 8413, .8643, .8849, .9032, .9192,
5.9332, .9452, .9554, .9641, .9713, .9772, .9821, .9861,. 9893, .9918, .9938,
6.9953, .9965, .9974 , .9981 , .9987 , .9990, .9993, -9995 , .9997, .9998, .9998,
7.9999, -9999/
DATA PUZ /O. 15, 0.15, 0.20, 0.25,0. 40, 0.60/
OA AREA SOURCE STRENGTH. (G/SEC-M**2)
HQ AREA SOURCE HEIGHT. (M)
RQ EAST COORDINATE OF S.W. CORNER. (KM)
SQ NORTH COORDINATE OF S.W. CORNER. (KM)
DEST EAST-WEST SIZE. (KM)
DNOR NORTH-SOUTH SIZE. (KM)
SL IS THE Y VALUE (SOURCE UNITS) OF THE INTERSECTION OF LINES
XUSE AND R
RL IS THE X VALUE (SOURCE UNITS) OF THE INTERSECTION OF LINES
XUSE AND S
SL(R)=( (XUSE- ((R-RREC)*SINT)) /COST )+SREC
HL(S)=((XUSE-((S-SREC)«COST))/SINT)+RREC
X IS THE UPWIND DISTANCE (SOURCE UNITS) OF SOURCE AT (R,S)
FROM RECEPTOR AT (RPS.SPS)
Y IS THE CROSSWIMD DISTANCE (SOURCE UNITS) OF SOURCE AT (R,S)
FROM RECEPTOR AT (RPS.SPS)
X(R,S)=(R-RREC)*SIMT+(S-SREC)*COST
Y(R,S)=(S-SREC)*SINT-(R-RREC)*COST
IWRI=6
DN IS NUMBER USED TO DETERMINE X (UPWIND FROM RECEPTOR) STEP
SIZE.
IDN=10
DN=IDN
DO 680 NW=1 ,NHR
THETA=WTHET(NW)
T=THETA/57.2958
T IS WIND DIRECTION (RADIANS)
SINT=SIN(T)
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00180
00190
00200
00210
00220
00230
00240
00250
00260
00270
00280
00290
00300
00310
00320
00330
00340
00350
00360
00370
00380
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
101
-------
AREA
COST=COS(T)
U=WU(NW)
UZ=U
KST=MKST(NW)
HL=WHL(NW)
PWR=PUZ(KST)
DO 670 NS=1,NQ
C DETERMINE FOUR BOUNDARIES OF AREA SOURCE. (IN SOURCE UNITS)
RB1=RQ(NS)
RB2=RBl+DEST(NS)
SB3=SQ(NS)
SB4=S63+DNOR(NS)
Q=QA(NS)
H=HQ(NS)
GO TO (20,10), IUZA
10 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR
IF (UZ.LT.1.0) UZ=1.0
20 DO 660 NC=1,NR
RREC=RR(NC)
SREC=SR(NC)
Z=ZR(NC)
C INITIALIZE VARIABLES FOR INTEGRATION PROCEDURE.
KDN=IDN
KNT=0
TEST=1.0
TOT1=0.
TOT2=0.
DMLX=0.
C KROL IS USED FOR PROGRAM FLOW CONTROL (SEE STEP 590).
KROL=1
C DETERMINE UPWIND DISTANCE OF EACH SOURCE CORNER FROM RECEPTOR.
XL(1)=X(RB1,SB3)
XL(2)=X(RB1,SB4)
XL(3)=X(RB2,SB3)
XL(4)=X(RB2,SB4)
C NEXT 8 STATEMENTS DETERMINE MAX AND MIN X (UPWIND) DISTANCES
C OF SOURCE CORNERS.
XMAX=XL(1)
XMIN=XL(1)
DO 60 L=2,4
IF (XMAX-XL(L)) 30,40,40
30 XMAX=XL(L)
40 IF (XMIN-XL(L)) 60,60,50
50 XMIN=XL(L)
60 CONTINUE
IF (XMAX) 660,660,70
C NO PART OF AREA SOURCE
C DETERMINE LIMITS FOR X
70 IF (XMIN) 30,90,90
80 XMIN=0.
C DELX IS INTEGRATION X WIDTH(SOURCE
90 DELX=(XMAX-XMIN)/DN
IS UPWIND OF RECEPTOR (660)
ITERATION. (IN SOURCE UNITS.)
UNITS) DMLX IS SAME(METERS)
00540
00550
00560
00570
00580
00590
00600
00610
00620
00630
00640
00650
00660
00670
00680
00690
00700
00710
00720
00730
00740
00750
00760
00770
00780
00790
00800
00810
00820
00830
00840
00850
00860
00870
00880
00890
00900
00910
00920
00930
00940
00950
00960
00970
00980
00990
01000
01010
01020
01030
01040
01050
01060
102
-------
AREA
DMLX=DELX*1000. 01070
C XUSE IS THE CURRENT VALUE OF UPWIND DISTANCE FROM THE RECEPTOR 01080
C FOR THIS CALC. OF CONG. (A NUMBER OF ITERATIONS ARE DONE FOR 01090
C THE TRAPEZOIDAL INTEGRATION.) 01100
XUSE=XMIN OHIO
C JDN IS CONTROL FOR PROPER INCLUSION OF FIRST AND LAST TERMS OF 01120
C TRAPEZOIDAL INTEGRATION. 01130
JDN=0 01140
100 TOT=0. 01150
C DETERMINE TWO LOCI.(IN SOURCE UNITS) TO 215 + 1 01160
IF (XUSE) 510,560,110 01170
110 1=0 01180
IF (COST+0.0001) 130,120,120 01190
120 IF (COST-0.0001) 160,160,130 01200
130 Sl=SL(RBl) 01210
IF (S1-SB3) 160,140,140 01220
140 IF (S1-SB4) 150,150,160 01230
150 1=1+1 01240
RA(I)=RBl 01250
SA(I)=S1 01260
160 IF (SINT+0.0001) 180,170,170 01270
170 IF (SINT-0.0001) 210,210,180 01280
IbO R3=RL(SB3) 01290
IF (R3-RB1) 210,190,190 01300
190 IF (R3-RB2) 200,200,210 01310
200 1=1+1 01320
RA(I)=R3 01330
SA(I)=SB3 01340
IF (1-2) 210,330,210 01350
210 IF (COST+0.0001) 230,220,220 01360
220 IF (COST-0.0001) 260,260,230 01370
230 S2=SL(RB2) 01380
IF (S2-SB3) 2bO,240,240 01390
240 IF (S2-SB4) 250,250,260 01400
250 1=1+1 01410
RA(I)=RB2 01420
SA(I)=S2 01430
IF (1-2) 260,330,260 01440
260 IF (SINT+0.0001) 280,270,270 01450
270 IF (SINT-0.0001) 310,310,280 01460
280 R4=RL(SB4) 01470
IF (R4-RB1) 310,290,290 01480
290 IF (R4-RB2) 300,300,310 01490
300 1=1+1 01500
RA(I)=R4 01510
SA(I)=SB4 01520
310 CONTINUE 01530
IF (1-2) 550,330,320 01540
320 IERR=215 01550
GO TO 52U 01560
C DETERMINE CROSSWIND DISTANCE OF THE TWO LOCI. 01570
330 RDUM=RA(1) 01580
SDUM=SA(1) 015yO
103
-------
AREA
YA=Y(RDUM,SDUM) 01600
RDUM=RA(2) 01610
SDUM=SA(2) 01620
YB=Y(RDUM,SDUM) 01630
C YA AND YB IN SOURCE UNITS. 01640
XDUM=XUSE 01650
C XDUM IS UPWIND DISTANCE (IN KM.). 01660
C FIND DISPERSION PARAMETER VALUES FOR THIS DISTANCE. 01670
CALL PGSIG (XDUM,XDUM,KST,SY,SZ) 01680
C EXPRESS THE CROSSWIND DISTANCE (Y) IN NORMAL DISTRIBUTION 01690
C DEVIATES. (SMAX AND SMIN) 01700
IF (YA-YB) 350,550,340 01710
340 SMAX=YA*1000./SY 01720
SMIN=YB*1000./SY 01730
GO TO 360 01740
350 SMAX=YB*1000./SY 01750
SMIN=YA*1000./SY 01760
C DETERMINE THE CUMULATIVE VALUE OF THE NORMAL DIST. FOR EACH S. 01770
3bO IF (SMAX+3.8) 550,550,370 01780
370 IF (SMIN-3.8) 380,550,550 01790
380 IF (3MAX-3.8) 400,400,390 01800
390 CFMAX=1.0 01810
GO TO 410 01820
400 SSD=SMAX*10.+40. 01830
ISD=SSD 01840
SSO=ISD 01850
CFMAX=GA(ISD)+(SSD-SSO)*(GA(ISD+1)-GA(ISD)) 01860
410 IF (SMIN+3.8) 420,420,430 01870
420 CFMIN=0.0 01880
GO TO 440 01890
430 SSD=SMIN*10.+40. 01900
ISD=SSD 01910
SSO=ISD 01920
CFMIN=GA(ISD)+(SSD-SSO)*(GA(ISD+1)-GA(ISD)) 01930
440 CONTINUE 01940
C DETERMINE THE AREA (CFM) OF THE NORMAL DIST. BETWEEN SMIN AND Oly50
C SMAX. 01960
CFM=CFMAX-CFMIN 01970
IF (CFM) 460,450,470 01980
450 CWI=0. 01990
GO TO 480 02000
460 IERR=405 02010
GO TO 520 02020
C DETERMINE THE CROSS-WIND INTEGRATED CONCENTRATION FOR THIS X. 02030
470 CALL RCONCA (Z,H,HL,XDUM,KST,AN,M,SZ,CWI) 02040
C CWI IS THE CONCENTRATION FROM AN INFINITE LINE SOURCE 02050
C IN SECONDS PER SQUARE METER. 02060
480 CONTINUE 02070
GO TO (490,540), KROL 02080
C STATEMENTS 490-530 PROVIDE FOR ADDITION OF THE TERMS IN THE 02090
C FORMULA FOR TRAPEZOIDAL INTEGRATION THAT CORRESPOND TO THE 02100
C ENDPOINTS OF THE INTEGRATION. 02110
490 IF (XUSE-XMIN) 510,530,500 02120
104
-------
AREA
500
510
520
530
C
C
540
C
550
560
570
580
590
C
600
610
C
C
620
630
640
C
C
650
660
670
680
C
690
700
C
TERMS FOR CURRENT INTEGRATION
FINITE LENGTH OF SOURCE.
APPROXIMATION,
IF (JDN-KDN) 540,530,510
IERR=430
GO TO 520
WRITE (IWRI.690) IERR
CALL EXIT
CWI=0.5*CWI
TOT IS A SUBTOTAL OF INTEGRATION
(SEC/METER**2) ADJUSTED FOR THE
TOT=TOT+CWI*CFM
JDN = KDN TERMINATES THIS INTEGRATION
IF (JDN-KDN) 560,570,580
XUSE=XUSE+DELX
JDN=JDN+1
GO TO 110
IF (TOT) 580,660,590
IERR=475
GO TO 520
GO TO (600,610), KROL
PROGRAM FLOW REACHES 600 ONLY ON FIRST INTEGRATION ATTEMPT,
TOT1=TOT*DMLX
JDN=1
KNT=2
XUSE=XMIN+DELX/2.
KROL=2
GO TO 100
TOTZ=TOT1/2.+(TOT*DMLX)/2.
TOT1 IS PREVIOUS ESTIMATE OF INTEGRATION (SEC/METER).
TOT2 IS'CURRENT ESTIMATE OF INTEGRATION (SEC/METEH).
TEST=A8S(TOT2-TOT1)/TOT2
CHECK PRECISION OF INTEGRATIONS AGAINST DESIRED VALUE.
IF (TEST-PIN) 650,650,620
IF (KNT-10) 640,640,630
WRITE (IWRI.700) IOT1,TOT2
IERR=500
GO TO 520
DELX=DELX/2.
XUSE=XMIN+DELX/2.
DMLX=DMLX/2.
TOT1=TOT2
KNT=KNT+1
JDN=1
KDN=2*KDN
RETURN TO CALCULATE ADDITIONAL VALUES FOR X DISTANCE
MIDWAY BETWEEN X DISTANCES USED FOK LAST ITERATION.
GO TO 100
CONTINUE
CONA(NC,NW)=CONA(NC,NW)+TOT2*Q*DVA(NW)/UZ
CONTINUE
CONTINUE
CONTINUE
RETURN
FORMAT
FORMAT
END
(1HO,'ERROR AT',15)
(1HO,'INTERVAL HAS BEEN
HALVED 10 TIMES',2E15.8)
105
02130
02140
02150
02160
02170
02180
02190
02200
02210
02220
02230
02240
02250
02260
02270
02230
02290
02300
02310
02320
02330
02340
02350
02360
02370
02380
02390
02UOO
02410
02420
02430
02440
02450
02460
02470
02480
02490
02500
02510
02520
02530
02540
02550
02560
02570
02580
02590
02600
02610
02620
02630
02640
02650
02660
02670
02630
-------
HRZLN
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
r
\s
c
c
c
c
c
c
r1
-------
HRZLN
QL=QLN(NS) 00b40
H=HLN(NS) 00550
GO TO (20,10) , IUZH 00560
10 UZ=U*(H/UHGT)**PWR 00b70
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR 00580
IF (UZ.LT.1.0) UZ=1.0 OOSyO
20 IF (SYO(NS)) 40,30,40 OObOO
30 KVYL=0. 00610
GO TO 50 00620
40 SYON=SYO(NS) 00630
XVYL=XVY(SYON,KST) 00640
50 IF (SZO(NS)) 70,60,70 006^0
60 XVZL=U. 00660
GO TO 80 00670
70 SZON=SZO(i\IS) 00680
XVZL=XVZ(SZON,KST) 00690
80 CONTINUE 00700
DO 920 NC=1,NR 00710
RREC=RR(NC) 00720
SREC=SR(NC) 00730
Z=ZR(NC) 00740
X1=X(R1,S1) 00750
X2=X(R2,S2) 00760
IF (XI) yO,110,110 00770
90 IF (X2) 100,110,110 00780
100 RC=0. 00790
GO TO 920 00800
110 DELR=R2-R1 00810
DELS=S2-S1 00820
Y1=Y(R1,S1) 00830
Y2=Y(R2,S2) 00840
IF (Y1-Y2) 120,420,120 00650
C IF Yl = Y2, LINE SOURCE IS PARALLEL TO UPWIND AZIMUTH FROM RECE 00860
120 IF (COST+0.0001) 190,130,130 00870
130 IF (COST-0.0001) 140,140,190 00880
140 IF (DELR+0.0001) 170,150,150 00a90
150 IF (DELR-0.0001) 160,160,170 00900
160 SLOC=SREC OOS10
RLOC=R1 00920
GO TO 310 00930
170 SLP=DELS/DELR 00940
IF (SLP) 180,420,180 00950
180 SLOC=SREC 00960
RLOC=(SLOC-S1)/SLP+R1 00970
GO TO 310 00980
190 IF (SINT+0.0001) 240,200,200 00990
200 IF (SINT-0.0001) 210,210,240 01000
210 IF (DELR+0.0001) 230,220,220 01010
220 IF (DELR-0.0001) 420,420,230 01020
230 SLP=DELS/DELR 01030
RLOC=RREC 01040
SLOC=SLP*(RLOC-R1)+S1 01050
GO TO 310 01060
107
-------
HRZLN
240 IF (DELR+0.0001) 270,250,250 01070
250 IF (DELR-0.0001) 260,260,270 01080
260 RLOC=R1 01090
SLOC=(RLOC-RREC)*COST/SINT+SREC 01100
GO TO 310 OHIO
270 IF (DELS+0.0001) 300,280,280 01120
280 IF (DELS-0.0001) 290,300,300 01130
290 SLOC=S1 01140
RLOC=(SLOC-SREC)*SINT/COST+RREC 01150
GO TO 310 01160
300 TATH=SINT/COST 01170
C TATH IS TANGENT (THETA) 01180
SLP=DELS/DELR 01190
C SLP IS SLOPE OF LINE SOURCE. 01200
RLOC=(RREC/TATH+S1-SLP*R1-SREC)/(1/TATH-SLP) 01210
SLOC=(RLOC-RREC)/TATH+SREC 01220
C RLOC, SLOC IS LOCUS OF LINE THROUGH LINE SOURCE AND UPWIND VECT 01230
310 XLOC=X(RLOC,SLOC) 01240
IF (XLOC) 420,420,320 01250
C XLOC IS POSITIVE IF LOCUS IS UPWIND. 01260
320 IF (S2-S1) 330,330,340 01270
330 SMAX=S1 01280
SMIN=S2 01290
GO TO 350 01300
340 SMAX=S2 01310
SMIN=S1 01320
350 IF (R2-R1) 360,360,370 01330
360 RMAX=R1 01340
RMIN=R2 01350
GO TO 3bO 01360
370 RMAX=R2 01370
RMIN=Rl 01380
C SEE IF UPWIND LOCUS IS ON LINE SOURCE. 01390
380 IF (RLOC-RMIN) 420,390,390 01400
390 IF (RMAX-RLOC) 420,400,400 01410
400 IF (SLOC-SMIiO 420,410,410 01420
410 IF (SMAX-SLOC) 420,430,430 01430
420 INDIC=1 01440
C INDIC =1 FOR NO LOCUS ON LINE SOURCE. 01450
XA=X1 01460
YA=Y1 01470
XB=X2 01480
YB=Y2 01490
GO TO 440 01500
430 INDIC=2 01510
C INDIC =2 FOR LOCUS ON LINE SOURCE. 01520
XA=X1 01530
YA=Y1 01540
XB=XLOC 01550
YB=0. 01560
440 DISX=XB-XA 01570
DISY=YB-YA 01580
DISI=3QRT(DISX*DI3X+DISY*DISY) 01590
108
-------
HRZLN
C DISI IS LENGTH(KM) OF LINE CONSIDERED.
IF (DISI) 460,450,460
450 CURR=0.
GO TO 660
460 DDI=DI3I*1000./20.
C ONE-HALF IS INCLUDED IN THE 20.
C DDI IS ONE-HALF TIMES 1/10 OF DISI (M).
DX=DISX/10.
DY=DISY/10.
PRE\7=0.
KNTRL=1
XI=XA
YI = YA
KNT=0
DO 580 1=1,11
C STORE EACH XI,YI.
XST(I)=XI
YST(I)=YI
IF (XST(I)) 470,470,480
470 RC=0.
GO TO 490
480 XZ=XI+XVZL
XY=XI+XVYL
CALL RCONCP (Z,H,HL,XZ,XY,YI,KST,AN,M,SY,SZ,RC)
490 GO TO (500,540), KNTRL
C IF RC IS ZERO, CONTINUE UNTIL RC IS POSITIVE,
500 IF (RC) 570,570,510
510 IF (1-1) 520,520,530
520 KNTRL=2
GO TO 560
C RESET POINT A TO LAST ONE PREVIOUS.
530 XA=XST(I-1)
YA=YST(I-1)
KNTRL=2
GO TO 560
540 IF (RC) 550,550,560
C RESET POINT B IF REACH ZERO CONCENTRATION.
550 XB=XI
YB=YI
GO TO 590
560 KNT=KNT+1
570 XI=XI+DX
YI=YI+DY
580 CONTINUE
590 IF (KNT) 610,610,600
bOO IF (KNT-6) 440,440,650
C IF GET TO 610, CONC. FROM THIS SEGMENT IS 0.
610 GO TO (620,630,640), INDIC
620 RC=0.
GO TO 920
630 FIRST=0.
GO TO 890
640 RC=FIRST
01600
01610
01620
01630
01640
01650
01660
01&70
01680
01690
01700
01710
01720
0173U
01740
01750
01760
01770
01780
01790
01800
01810
01820
01830
01840
01850
01860
01870
01880
01890
01900
01910
01920
01930
01940
01950
01960
01970
01980
01990
02000
02010
02020
02030
02040
02050
02060
02070
02080
02090
02100
02110
02120
109
-------
HRZLN
t> 50
C
C
C
C
660
670
680
690
700
710
720
730
740
C
750
C
GO TO
CONTINUE
DO AiMOTHER TRAPEZOIDAL INTEGRATION PROM A TO B IN TEN STEPS,
IT IS LIKELY THAT A TO B HAVE BEEN REDEFINED.
DISX=XB-XA
DISY=YB-YA
DISI=S2RT(DI3X*DISX+DISY*DISY)
DISI IS DISTANCE (KH) FROM A TO B
DELD=DISI*10U.
DELD IS 1/10 DISI IN METERS.
DX=DISX/10.
DY=DI3Y/10.
760
XDIM=XA
YDUM=YA
IF (XDUM) 660,660,670
RC=0.
GO TO 680
XZ=XDUM+XVZL
XY=XDUM+XVYL
CALL RCONCP ( Z , H , HL,XZ ,XY , YDUM ,KST, AN , M,SY ,SZ ,RC )
SUM=SUM+RC/2.
DO 710 1=1,9
XDUM=XDUM+DX
YDUM=YDUM+DY
IF (XDUM) 6yO,690,700
RC=0.
GO TO 710
XZ=XDUM+XVZL
XY=XDUM+XVYL
CALL RCONCP ( Z ,H , HL,XZ ,XY ,YDUM ,KST, AN, M,SY ,SZ ,RC )
SUM=SUM+RC
XDUM=XDUM+DX
YDUM=YDUM+DY
IF (XDUM) 720,720,730
RC=0.
GO TO 740
XZ=XDUM+XVZL
XY=XDUM+XVYL
CALL RCONCP ( Z ,H,HL,XZ ,XY, YDUM ,KST, AN, M,SY ,SZ ,RC )
SUM=SUM+RC/2.
INTEGRATED VALUE IS CURR.
CURR=SUM*DELD
ILIM=10
PREV=CURR
EVALUATE FOR POINTS IN BETWEEN THOSE ALREADY EVALUATED.
DELD=DELD/2.
XDUM=XA+DX/2.
YDUM=YA+DY/2.
DO 790 I=1,ILIM
IF (XDUM) 760,760,770
RC = 0.
GO TO 780
02130
02140
02150
02160
02170
02160
02190
02200
02210
02220
02230
02240
02250
02260
02270
02280
02290
02300
02310
02320
02330
02340
02350
02360
02370
02380
02390
02400
02410
02420
02430
02440
02450
02460
02470
02480
02490
02500
02510
02520
02530
02540
02550
02560
02570
02580
02590
02600
02610
02620
02630
02640
02650
110
-------
HRZLN
770
C
780
790
C
800
C
810
820
830
840
850
C
350
870
880
390
900
9 1 0
920
930
940
XZ=XDUM+XVZL
XY=XDUM+XVYL
CALL RCOHCP ( Z ,H, HL ,XZ ,XY , YDUM ,KST, AN , M,SY ,SZ ,RC)
NOTE ADD THESE TO fiC'S FOUND ABOVE.
SUM=SUM+HC
XDUM=XDUM+DX
YDUM=YDUM+DY
CURR=SUM*DELD
TEST=ABS( (CURR-PREV)/CURR)
IF WITHIN PIN OF LAST i/ALUE (PREV), CONSIDER THIS AS FINAL VALU
IF (TEST-PIN) 860,800,800
ILIM=ILIM*2
PREV=CURR
EVALUATE POINTS IN BETWEEN.
DELD=DELD/2.
DX=DX/2.
DY=DY/2.
XDUM=XA+DX/2.
YDUM=YA+DY/2.
DO 840 1=1 ,ILIM
IF (XDUH) 810,810,820
RC=0.
GO TO 830
XZ=XDUM+XVZL
XY=XDUM+XVYL
CALL RCONCP ( Z , H , HL , XZ ,XY , YDUM , KST, AN , M,SY ,SZ , RC)
XDUM=XDUM+DX
YDUM=YDUM+DY
CURR=SUM*DSLD
TEST=ABS((CUrtR-PREV)/CURR)
IF (TEST-PIN) 860,850,850
ILIM=ILIM*2
DX=DX/2.
DY=DY/2.
GO TO 750
AT 860 rlAVt FINAL VALUE
GO TO (370,380,900), INDIC
RC=CURH
GO TO 910
FIRST=CURH
OF INTEGRATION IN CURR
XA=XLOC
YA=0.
Xb=X2
YB=Y2
GO TO 440
RC=FIRST+CURR
CONLH ( NC , NW ) =CONLH( HC , M ) +RC*QL*DVH (KU) /UZ
CONTINUE
CONTINUE
CONTINUE
RETURN
in
02660
02670
02680
02690
02700
02710
02720
02730
02740
02750
02760
02770
02780
02790
02800
32810
02820
02830
02840
02850
0286G
02870
02880
02690
02900
02910
02920
02930
02940
02950
02960
02970
02980
02990
03000
03010
03020
03030
03040
03050
03060
03070
03080
03090
03100
03110
03120
03130
03 1 40
03150
03160
03170
03130
03190
03200
-------
CRVLN
SUBROUTINE CRVLN (NHR , NQ , NR , PIN)
C CRVLN CALCULATES HOURLY CONCENTRATIONS FROM
C MULTI-LANE CURVED PATH SOJRCES.
c THE: FOLLOWING SUBROUTINES ARE CALLED BY CRVLN.
C CURLIi^ —CALCULATES THE INTERSECTION OF THE
C LINE FORMED BY THE WIND DIRECTION AND
C RECEPTOR COORDINATES AND THE CIRCLE
C DETEhMINED FROM POINTS A,B,C.
C RCONCP--DSTERMINES THE RELATIVE CONCENTRATION AT
C A RECEPTOR FROM A POINT AT A GIVEN UPWIND
C AND CROSSWIMD DISTANCE.
C THE FOLLOWING FUNCTIONS ARE CALLED BY CRVLN.
C XVY--CALCULATliS THE VIRTUAL DISTANCE NECESSARY
C TO ACCOUNT FOR THE INITIAL CROSSWIMD DISPERTION.
C XVZ---CALCULATES THE VIRTUAL DISTANCE NECESSARY
C TO ACCOUNT FOR THE INITIAL VERTICAL DISPERTION.
C ANGARA-DETERMINES THE APPROPRIATE ARCTAM(M/N)
C WITH THE RESULTING ANGLE BETWEEN 0 AND
C 360 DEGREES.
C DIFAHG—DETERMINES THE DIFFERENCE BETWEEN TWO ANGLES.
C THE RESULTING ANGLE IS BETWEEN 0 AND 360 DEGREES,
COMMON /SOCP/ QLNS(31),HLNS(31),RBQS(31),SBQS(31),RMQS(31),SMQS(31
1) ,REQS(31 ) ,SEQS(3D ,SIYO(31) ,SIZO(3D ,CONCLN(31 ,24) ,DDVH(3O ,IUZC,
2RADIU3(31) ,NLANE(3O
COMMON /REG/ RR(31),SR(31) , ZR(31)
COMMON /WLA/ WTHET(25),WU(25),MK3T(25),WHL(25),WTA(25),UHGT
DIMENSION PUZ(6), RRR(11), SSS(11), DIR(11), ARC(5), XANG(5)
C
c
c
c
c
c
c
c
c
ft
u
c
c
c
c
c
c
c
c
c
c
c
c
DATA P
QLHS
HLNS
RBQS
SBQS
RMQS
SMQS
REQS
SEOS
SIYO
SIZO
CONCLN
DDVH
IUZC
RR EA
SR NO
ZR
WTHET
WU WI
MKST
WHL M
WTA A
WP A
/O.15,0.15,0.20,0.25,0.40,0.607
PATH SOURCE STRENGTH FOR EACH LANE
HEIGHT OF THE PATH SOURCE
EAST COORDINATE,POINT A
NORTH COORDINATE,POINT A
EAST COORDINATE,POINT B
NORTH COORDINATE,POINT B
EAST COORDINATE,POINT C
NORTH COORDINATE,POINT C
INITIAL SIGMA Y
INITIAL SIGMA Z
CONCENTRATION
VARIATION,CURVED
PATH SOURCE
DIURNAL
WIND INCREASE WITH HEIGHT
COORDINATE OF RECEPTOR
! COORDINATE OF RECEPTOR
HEIGHT OF RECEPTOR ABOVE GROUND
WIND DIRECTION
D SPEED
STABILITY CLASS
NG HEIGHT
SNT AIR SURFACE TEMP.
AMBIENT AIR SURFACE PRESSURE
X(R,S)=(R-RREC)*SINT+(S-SREC)*COST
X IS UPWIND DISTANCE OF R,S FROM RREC.SREC
Y(R,S)=(S-SREC)*SINT-(R-RREC)*COST
(G/SEC-M)
(M)
(KM)
(KM)
(KM)
(KM)
(KM)
(KM)
(M)
(M)
(G/M**3)
(DIMENSIONLESS)
(DIMENSIONLESS)
(KM)
(KM)
(M)
(DEC AZIMUTH)
(M/SEC)
(DIMENSIONLESS)
(M)
(DEC K)
(MB)
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
001 10
00120
00130
00140
00150
00160
00170
00180
00190
00200
00210
00220
00230
00240
00250
00260
00270
00280
00290
00300
00310
00320
00330
00340
00350
00360
00370
00380
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
112
-------
CRVLN
C Y IS CROSSWTND DISTANCE OF R,S FROM RREC,SREC
C START DO LOOP FOR EACH HOUR.
DO 460 NW=1,NHR
KUMC=0
THETA=WTHET(NW)
TR=THETA/57.2958
SINT=SIN(TR)
COST=COS(TR)
U=WU(NW)
UZ=U
KST=MKST(Nto)
HL=WHL(NVv)
PWR=PUZ(KST)
C START DO LOOP FOR EACH CURVED PATH.
DO 460 NS=1,NQ
RB=RBQS(NS)
RM=RMQS(NS)
RE=REQS(NS)
SB=SBQS(NS)
SW=SMQS(NS)
SE=SEQS(iSIS)
H=HLNS(NS)
GO TO (20,10), IUZC
10 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(20G./UHGT)**PrtR
IF (UZ.LT.1.0) UZ=1.0
20 IF (SIYO(NS)) 40,30,40
30 XVYL=0.
GO TO 50
40 SYON=SIYO(NS)
XVYL=XVY(SYON,KST)
50 IF (SIZO(NS)) 70,60,70
bO XVZL=0.
GO TO 80
70 SZON=SIZO(NS)
XVZL=XVZ(SZON,KST)
dO CONTINUE
C THE COORDINATES OF THE CENTER OF THE CIRCLE,
C DETERMINED FROM THE INTERSECTION OF THE PERPENDICULAR
C BISECTORS OF CHORDS AS AND BC,AND THE RADIUS OF THE
C CIRCLE ARE CALCULATED.
RMID=(RB+RM)/2.
SMID=(3B+SM)/2.
RNID=(Rh+RE)/2.
SNID=(3M+SE)/2.
IF (RM-RB.EQ.O) GO TO yO
IF (RE-RM.Es2.U) GO TO 100
C SL01 IS THE SLOPE OF THE PERPENDICULAR BISECTOR
C OF CHORD AB.
C SL02 IS THE SLOPE OF THE PERPENDICULAR BISECTOR
C OF CHORD BC.
SL01—( (RM-RB)/(SM-SB) )
SL02=-((RE-RM)/(SE-SM))
00540
00550
00560
00570
00580
00590
00600
00610
00620
00630
00640
00650
00660
00670
00680
00690
00700
00710
00720
00730
00740
00750
00760
00770
007HO
00790
00800
00810
00820
00830
00840
00850
00860
00b70
00880
00890
00900
00910
00920
00930
00940
00950
00960
00y70
00980
00990
01000
0101J
01020
01030
01040
01050
01060
113
-------
CRVLN
C RC AND SC ARE THE COORDINATES OF THE CENTER OF THE
C CIRCLE.
RC=(SL01*RMID+SNID-SMID-SL02*RNID)/(SL01-SL02)
SC=SL01*(RC-RMID)+SMID
GO TO 110
90 RC=RMID
SL02=-((RE-RN)/(SE-SM))
SC=SL02*(RC-RNID)+SNID
GO TO 110
100 RC=RNID
SL01=-((RM-RB)/(SM-SB))
SC=SL01*(RC-RMID)+SMID
GO TO 110
C RAD IS THE RADIUS OF THE CIRCLE.
110 RAD=SQRT((RC-RB)**2+(SC-SB)**2)
DELM=RB-RC
DELN=SB-SC
C ANGB IS THE DIRECTED ANGLE FROM THE CENTER TO POINT A.
ANGB=ANGARC(DELM,DELN)
DELM=RM-RC
DELN=SM-SC
C ANGM IS THE DIRECTED ANGLE FROM THE CENTER TO POINT 8,
ANGM=ANGARC(OELM,DELN)
DELrf=RE-RC
DELN=SE-SC
C ANGE IS THE DIRECTED ANGLE FROM THE CENTER TO POINT C,
ANGE=ANGARC(DELM,DELN)
ANGBE=ABS(DIFANG(ANGB,ANGE))
ANGBM=A8S(DIFANG(ANGB,ANGM))
ANGME=ABS (DIFANG(ANGM, ANGE) )
IF (ANGBM.LT.ANGBE.AND.ANGME.LT.ANGBE) GO TO 120
GAMDEG=360.-ANGBE
GO TO 130
120 GAWDEG=ANGBM+ANGME
130 A=l.
IF (GAMDEG.GE.180.) GO TO 140
BEDEG=DIFANG(ANGE,ANGB)
IF (BEDEG.LT.O.) A=-l.
GO TO 150
140 RDEG=ANGB+180.
RMDEG=DIFANG(ANGM,RDEG)
REDEG=DIFANG(ANGE,RDEG)
IF (RMDEG.GE.O..AND.REDEG.LE.0.) A=-l.
150 CONTINUE
NUMLNE=NLANE(NS)
NUMC=NUMLNE+NUMC
NLL=NUMC-NUMLNE+1
RDI=RAD
DO 460 NL=NLL,NUMC
RAD=RDI+RADIUS(NL)
QL=QLNS(NL)
DO 460 NC=1,NR
RREC=RR(NC)
01070
01080
01090
01100
OHIO
01120
01130
01140
01150
01160
01170
01180
01190
01200
01210
01220
01230
01240
01250
01260
01270
01280
01290
01300
01310
01320
01330
01340
01350
01360
01370
01380
01390
01400
01410
01420
01430
01440
01450
01460
01470
01480
01490
01500
01510
01520
01530
01540
01550
01560
01570
01580
01590
114
-------
CRVLN
SREC=SR(NC)
Z=ZR(NC)
C CURLIN CALCULATES THE COORDINATES OF ANY UPWIND LOCI ON ARC AB
C FROM HERE TO STATEMENT 310 THE BEGINNING AND ENDING
C POINTS,WITH THEIR CORRESPONDING DIRECTIONS,FOR EACH SEGMENT
C OF THE ARC ARE CALCULATED
CALL CURLIN (RC,SC,RREC,SREC,Z,RAD,THETA,RL,SL,RLl,SLl,NLOCI)
IF (NLOCI) 160,170,160
C NLOCI HAS VALUES OF -1,0,OR 1 FOR:
C NO UPWIND LOCI ON ARC AB: 1 LOCI ON ARC AB:
C OR 2 LOCI ON ARC AB
160 XANG(1)=ANGB
ARC(1)=GAMDEG
J=l
GO TO 310
170 K=l
GO TO 190
IdO K=2
190 N=l
DO 270 KK=1,K
IF (KK.EQ.2) GO TO 200
DELM=RL-RC
DELN=SL-SC
ANGRL=ANGARC(DELM,DELN)
RANG=ANGRL
GO TO 210
200 DELM=RL1-RC
DELN=SL1-SC
ANGRL=ANGARC(DELM,DELN)
RLANG=ANGRL
210 IF (GAMDEG.GT.180..AND.A.GT.O. ) GO TO 220
IF (GAMDEG.GE.180..AND.A.LT.O.) GO TO 230
IF (GAMDEG.LT.ldO..AND.A.LT.O. ) GO TO 240
BIN=DIFANG(ANGRL,ANGB)
IF (BIN.LE.O..OR.BIN.GE.GAMDEG) GO TO 260
ARC(N)=BIN
GO TO 250
220 BIN=DIFANG(ANGRL,ANGB)
IF (6IN.GT.O.) ARC(N)=BIN
IF (BIN.GT.O.) GO TO 250
IF (BIN.LE.O.) DIF=(GAMDEG-360.)-BIN
IF (DIF.LT.O.) GO TO 2bO
ARC(N)=BIN+360.
GO TO 250
230 BIN=DIFANG(ANGRL,ANGB)
IF (BIN.LT.O.) ARC(N)=BIN
IF (BIN.LT.O.) GO TO 250
IF (BIN.GE.O..AND.BIN.LE.360.-GAMDEG) GO TO 260
ARC(N)=360.-BIN
GO TO 250
240 BIN=DIFANG(ANGRL,ANGB)
IF (BIN.GE.O..OR.BIN.LE.-GAMDEG) GO TO 260
ARC(N)=BIN
01600
01610
01620
01630
01640
01650
01660
01670
01680
01690
01700
01710
01720
01730
01740
01750
01760
01770
01780
01790
OldOO
01&10
01620
Old30
01840
01850
01660
01870
OlBdO
01390
01900
01910
01920
01930
01940
01950
01960
01970
01980
01990
02000
02010
02020
02030
02040
02050
02060
02070
02080
02090
02100
02110
02120
115
-------
CRVLN
GO TO 250
250 N=N+1
INDEX=-1
IF (KK.EQ.2) INDEX=1
GO TO 270
260 IF (K.NE.l) GO TO 270
XANG ( 1 ) =ANGB
ARC ( 1 ) =GAMDEG
270 CONTINUE
XANG ( 1 ) =ANGB
IF (N-l.EQ.O) GO TO 300
IF (N-1.EQ.2) GO TO 280
ARC ( 2 ) =GAMDEG-ABS ( ARC ( 1 ) )
XANG(2)=RANG
IF (INDEX.EQ.l) XANG(2)=RLANG
J=2
GO TO 310
280 ARC(1)=ABS(ARC(1) )
ARC(2)=ABS(ARC(2) )
IF (ARC(l) .GT.ARC(2) ) GO TO 290
ARC(2)=ARC(2)-ARC(1)
ARC ( 3 ) =GAMDEG-ARC ( 1 ) -ARC ( 2 )
XANG(2)=RANG
XANG(3)=RLANG
J=3
GO TO 310
290 AA=ARC(1)
ARC(1)=ARC(2)
ARC(2)=AA-ARC(2)
ARC ( 3 ) =GAMDEG-ARC ( 1 ) - ARC ( 2 )
XANG(2)=RLANG
XANG(3)=RANG
J = 3
GO TO 310
300 J=l
ARC(1)=GAMDEG
310 TOTAL=0.
DO 450 M=1,J
KNT=0
NM=0
LCOUNT=0
C CHI VALUES ARE CALCULATED
C THE SEGMENT IS REDUCED IN
C VALUES ARE NON-ZERO.
320 DO 350 L=l,ll
DIR(L)=XANG(M)+(A*(L-1)/10. )*ABS(ARC(M)
IF (DIR(L).GT.360.) DIR(L )=DIR(L)-360.
IF (DIR(L) .LT.O.) DIR(L)=DIR(L)+360.
DIRRAD=DIR(L)/57. 29578
RRR(L)=RC+RAD*SIN(DIRRAD)
SSS(L)=SC+RAD*COS(DIRRAD)
XZ=X(RRR(L) ,SSS(L) )+XVZL
IF (XZ-XVZL) 350,350,330
FOR 11 POINTS
SIZE UNTIL AT
ON THE ARC
LEAST 7 OF
SEGMENT.
THE 11 CHI
02130
02140
02150
021oO
02170
02160
02190
02200
02210
02220
02230
02240
02250
02260
02270
02280
02290
02300
02310
02320
02330
02340
02350
02360
02370
02380
02390
02400
02410
02420
02430
02440
02450
02460
02470
02480
02490
02500
02510
02520
02530
02540
02550
02560
02570
02580
02590
02600
02610
02620
02630
02640
02650
116
-------
CRVLN
330
340
C
C
350
360
C
C
370
C
380
390
400
XY=X(RRR(L) ,SSS(L) )+XVYL
YI=Y(RRR(L) ,SSS(L) )
CALL RCONCP ( Z ,H-,HL,XZ ,XY, YI ,KST, AN ,MI ,SY,SZ ,CHI )
IF (CHI) 350,350,340
NM=NM+1
THE NEXT TWO STEPS DETERMINE THE DIRECTION FROM
THE CENTER TO THE NEW STARTING POINT
IF (L.EQ.l) DIRINT=DIR(L)
IF (NM.EQ.l.AND.L.GT.l) DIRINT=DIR(L-1 )
IF (L.GT.l.AND.NM.EQ.l) NM=2
LCOUNT=0
IF (L.NE.ll) LCOUNT=1
CONTINUE
NM=NM+LCOUNT
IF (NM) 430,430,360
IF (NM.GT.7) GO TO 370
VAR=NM-1
ARC(M)=ABS(ARC(M) )*(VAR/10. )
XANG(M)=DIRINT
NM=0
GO TO 320
DO A TRAPEZOIDAL INTEGRATION FROM A TO B IN TEN STEPS,
IT IS LIKELY THAT A OR B HAVE BEEN REDEFINED.
11=11
L=l
KK=1
SUM=0 .
VAR=NM-1
ARR=ABS ( ARC ( M ) ) * ( VAR/1 0 . )
DELL=1/10 OF THE ARC LENGTH IN METERS
DELL=( (ARR/10. )/57 . 29578 ) *RAD*1000 .
DO 400 I=L,II,KK
DIRR=DIRINT+ ( A* ( 1-1 ) /DIV ) *ARR
DIRRAD=DIRR/57.29578
RRR1=RC+RAD*SIN(DIRRAD)
SSS1=SC+RAD*COS(DIRRAD)
XZ=X (RRR1 ,SSS1 )+XVZL
IF (XZ-XVZL) 400,400,390
XY=X(RRR1,SSS1)-I-XVYL
YI=Y(RRR1,SSS1)
CALL RCONCP ( Z ,H,HL,XZ ,XY , YI ,KST, AN, MI ,SY,SZ ,CHI )
IF (I.EQ.l) Ql=CHI/2.
IF (I.EQ.II) Q2=CHI/2.
SUM=SUM+CHI
CONTINUE
KNT=KNT+1
IF (KNT.GT.l) GO TO 410
TOTQ=SUM-Q1-Q2
TOTQ2=TOTQ*DELL
SUM=0 .
DELL=DELL/2.
11 = 21
02660
02670
02680
02690
02700
02710
02720
02730
02740
02750
02760
02770
02780
02790
02800
02810
02820
02830
02840
02850
02860
02870
02880
02890
02900
02910
02920
02930
02940
02950
02960
02970
02980
02990
03000
03010
03020
03030
03040
03050
03060
03070
03080
03090
03100
03110
03120
03130
03140
03150
03160
03170
03180
117
-------
410
420
430
440
450
460
C
CRVLN
L=2
KK=2
GO TO 380
CONTINUE
TOTQ=TOTQ+SUM
TOTQ1=TOTQ*DELL
TEST=ABS((TOTQ1-TOTQ2)/TOTQl)
IF (TEST-PIN) 440,420,420
TOTQ2=TOTQ1
DELL=DELL/2.
11=2*11-1
SUM=0.
GO TO 380
TOTQ1=0.
CONTINUE
TOTAL=TOTAL+TOTQ1
CONTINUE
CONCLN(NC,NW)=CONCLN(NC,NW)+TOTAL*QL*DDVH(NW)/UZ
CONTINUE
RETURN
END
03190
03200
03210
03220
03230
03240
03250
03260
03270
03280
03290
03300
03310
03320
03330
03340
03350
03360
03370
03380
03390
03400
118
-------
SPCLN
SUBROUTINE SPCLN (NHR,NQ,NR,PIN) 00010
C SPECIAL LINE SOURCE 00020
C THIS SUBROUTINE DETERMINES CONCENTRATIONS FROM 00030
C SPECIAL LINE SOURCES.CASES CONSIDERED INCLUDE 00040
C DIFFERENT HEIGHTS*ABOVE GROUND)OF END POINTS; 00050
C CONSTANT VELOCITY OF SOURCES ALONG LINE,OR 00060
C CONSTANT ACCELERATION(LINEAR VELOCITY CHANGE); 00070
C DECELERATION IS NEGATIVE ACCELERATION,AND 00080
C VARIATION OF WIND WITH HEIGHT. 00090
C THE FOLLOWING SUBROUTINES ARE CALLED BY SPCLN. 00100
C RCONCP--DETERMINES THE RELATIVE CONCENTRATION AT 00110
C A RECEPTOR FROM A POINT AT A GIVEN UPWIND 00120
C AND CROSSWIND DISTANCE. 00130
C THE FOLLOWING FUNCTIONS ARE CALLED BY SPCLN. 00140
C XVY—CALCULATES THE VIRTUAL DISTANCE NECESSARY 00150
C TO ACCOUNT FOR THE INITIAL CROSSWIND DISPERTION. 00160
C XVZ--CALCULATES THE VIRTUAL DISTANCE NECESSARY 00170
C TO ACCOUNT FOR THE INITIAL VERTICAL DISPERTION. 00180
COMMON /SOLS/ QLS(31),HAS(31),HBS(31),RAS(31),SAS(31),RBS(31),SBS( 00190
13D ,SYOS(3D ,SZOS(3D ,CONLS(31 ,24) ,DVS(25) ,IUZS ,SPDI(31) ,SPDF(3D , 00200
2TVSL(3D ,VSSL(3D 00210
COMMON /REG/ RR(31),SR(31),ZR(31) 00220
COMMON /WEA/ WTHET(25),WU(25),MKST(25),WHL(25),WTA(25),UHGT 00230
DIMENSION XST(11), Y3T(11), HST(11), PUZ(6) 00240
DATA PUZ /O.15,0.15,0.20,0.25,0.40,0.60/ 00250
C QLS LINE SOURCE STRENGTH. (G/SEC) 00260
C HAS HEIGHT OF POINT A. (M) 00270
C HBS HEIGHT OF POINT 8 (M) 00280
C RAS EAST COORDINATE,POINT A. (KM) 00290
C SAS NORTH COORDINATE,POINT A. (KM) 00300
C RBS EAST COORDINATE,POINT B. (KM) 00310
C SBS NORTH COORDINATE,POIMT B. (KM) 00320
C SPDI SPEED AT POINT A. (M/SEC) 00330
C SPDF SPEED AT POINT 3. (M/SEC) 00340
C SYOS INITIAL SIGMA Y. (M) 00350
C SZOS INITIAL SIGMA Z. (M) 00360
C TVSL VEHICLE VOLUME. (VEH/HR) 00370
C VSSL GROSS ESTIMATE OF VEH. SIZE. (M) 00330
X(R,S)=(R-RREC)*SINT-t-(S-SREC)*COST 00390
C X IS UPWIND DISTANCE OF R,S FROM RREC.SREC 00400
Y(R,S)=(S-SRfiC)*SIrlT-(R-RREC)*COST 00410
C Y IS CROSSWIND DISTANCE OF R,S FROM RREC.SREC 00420
DISTA(RM,SM,HM)=SQRT((RM-XS)«*2.+(SM-YS)«*2.+(HM-HGT)»*2.) 00430
DO 1080 NW=1,NHR 00440
C DO FOR EACH HOUR. 00450
THETA=WTriET(NW) 00460
TR=THETA/57.2958 00470
SINT=SIM(TR) 00480
COST=COS(TR) 00490
U=WU(NW) 00500
UZ=U 00510
KST=MKST(NW) 00520
00530
119
-------
SPCLN
P^R=PUZ(KST) 00540
DO 1070 NS=1,NQ 00550
C DO FOR EACH LINE SOURCE. 00560
R1=RAS(NS) 00570
S1=SAS(NS) 005dO
R2=RBS(NS) 00590
S2=3BS(NS) 00600
QL=QLS(NS)*(TVSL(NS)/3t>00.) 00610
SCALV=VSSL(NS)*(TVSL(NS)/3600.) 00620
Hl=HA3(NS) 00630
HGT=H1*0.001 00640
H2=HBS(NS) 00650
DELR=R2-R1 00660
DELS=S2-S1 00670
DIST=SQRT(DELR**2+DELS**2) 00680
DELH=H2-H1 00690
VI=SPDI(NS)**2 00700
VF=SPDF(NS)**2 00710
DELSPD=VF-VI 00720
DISTOL=SQRT((R2-R1)**2.+(S2-S1)**2.+(H2*0.001-Hl*.001)**2.) 00730
IF (SYOS(NS)) 10,10,20 00740
10 XVYL=0. 00750
GO TO 30 00760
20 SYON=SYOS(NS) 00770
XVYL=XVY(SYON,KST) 00780
30 IF (SZOS(NS)) 40,40,50 00790
40 XVZL=0. 00800
GO TO 60 00610
50 SZON=SZ03(NS) 00820
XVZL=XVZ(SZON,KST) 00830
60 CONTINUE 00640
DO 1060 NC=1,NR 00d50
C DO FOR EACH RECEPTOR. 00860
RREC=RR(NC) 00670
SREC=SR(NC) 00880
Z=ZR(NC) 00690
X1=X(R1,S1) 00900
XS=X1 00910
X2=X(R2,S2) 00920
IF (XI) 70,70,80 00930
70 IF (X2) 1060,1060,80 00940
60 Y1=Y(R1,S1) 00950
YS=Y1 00960
Y2=Y(R2,S2) 00970
C IF Yl AND Y2 ARE BOTH POSITIVE OR BOTH NEGATIVE NO LOCUS 00980
C ON THE LINE SOURCE ,GO TO 390. 00990
IF (Y1.GE.O.O.AND.Y2.GE.O.O) GO TO 390 01000
IF (Y1.LE.O.O.AND.Y2.LE.O.O) GO TO 390 01010
C IF Yl = Y2, LINE SOURCE IS PARALLEL TO UPWIND AZIMUTH FROM RECE 01020
IF (DIST) 90,390,90 01030
90 IF (COST+0.0001) 160,100,100 01040
100 IF (COST-0.0001) 110,110,160 01050
110 IF (DELR+0.0001) 140,120,120 01060
120
-------
SPCLN
C AT (110) COST APPROX. ZERO, WIND 90 DEC. OR 271). 01070
120 IF (DELR-0.0001) 130,130,140 01080
130 SLOC=SREC 01090
C AT(130) LINE SOURCE IS NORTH-SOUTH. 01100
RLOC=R1 OHIO
GO TO 260 01120
140 SLP=DELS/DELR 01130
C AT(140) WIND STILL 90 DEC. OR 270. 01140
IF (SLP) 150,390,150 01150
C IF SLP ZERO LINE PARALLEL TO WIND. OllbO
150 SLOC=SREC 01170
RLOC=(SLOC-S1)/SLP+R1 01180
GO TO 280 01190
IbO IF (SINT+0.0001) 210,170,170 01200
C AT(160) WIND NOT yO DEG. OR 270. 01210
170 IF (SINT-0.0001) 180,180,210 01220
IbO IF (DELR+0.0001) 200,l90,lyO 01230
C AT(180) wIND 0 DEG. OR 160. 01240
190 IF (DELR-O.U001) 390,390,200 01250
C GO TO 390 IF LINE SOURCE PARALLEL TO WIND. 012bO
200 SLP=DELS/DELR 01270
C AT(200) WIND 0 DEG. OR 160 DEG. AND DELR IS NOT ZERO. 01280
RLOC=RREC 01290
SLOC=SLP*(RLOC-R1)+31 01300
GO TO 280 01310
210 IF (DELR+0.0001) 240,220,220 01320
C AT(210) WIND NOT FROM 4 CARDINAL DIRECTIONS. 01330
220 IF (DELR-0.0001) 230,230,240 01340
230 RLOC=R1 01350
C AT(230) LINE SOURCE IS NORTH-SOUTH. 01360
SLOC=(RLOC-RREC)*COST/SINT+SREC 01370
GO TO 280 013«0
240 IF (DELS+0.0001) 270,250,250 01390
250 IF (DELS-0.0001) 260,270,270 01400
260 SLOC=S1 01410
C AT(260) LINE SOURCE IS EAST-WEST. 01420
RLOC=(3LOC-SREC)*SINT/COST+RREC 01430
GO TO 280 01440
270 TATH=SINT/COST 01450
C AT(270) GENERAL CASE,WIND NOT FROM CARDINAL DIRECTION 01460
C AND LINE SOURCE NOT EAST-WEST NOR NORTH-SOUTH. 01470
C TATH IS TANGENT (THETA) 01480
SLP=DELS/DELR 01490
C SLP IS SLOPE OF LINE SOURCE. 01500
RLOC=(RREC/TATH+S1-SLP*R1-SREC)/(1/TATH-SLP) 01510
SLOC=(RLOC-RREC)/TATH+SREC 01520
C RLOC, SLOC IS LOCUS OF LINE THROUGH LINE SOURCE AND UPWIND VECT 01530
280 XLOC=X(RLOC,SLOC) 01540
IF (XLOC) 390,390,290 01550
C XLOC IS POSITIVE IF LOCUS IS UPWIND. 01560
290 IF (S2-S1) 300,300,310 01570
300 SMAX=S1 01580
SMIN=S2 01590
121
-------
SPCLN
GO TO 320
310 SMAX=S2
SMIN=S1
320 IP (R2-R1) 330,330,340
330 RMAX=R1
RMIN=R2
GO TO 350
340 RMAX=R2
RMIN=R1
C SEE IF UPWIND LOCUS IS ON LINE SOURCE.
350 IF (RLOC-RMIN) 390,360,360
360 IF (RMAX-RLOC) 390,370,370
370 IF (SLOC-SMIN) 390,380,380
380 IF (SMAX-SLOC) 390,400,400
390 INDIC=1
C INDIC =1 FOR NO LOCUS ON LINE SOURCE.
XA=X1
YA=Y1
HA=H1
XB=X2
YB=Y2
HB=H2
GO TO 410
400 INDIC=2
C INDIC =2 FOR LOCUS ON LINE SOURCE.
XA=X1
YA=Y1
HA=H1
XB=XLOC
YB=0.
RDUM=RLOC-Rl
SDUM=SLOC-S1
DDUM=SQRT(RDUM*RDUM+SDUM*SDUM)
HB=(DDUM/DIST)*(H2-H1)+H1
HLOC=HB
410 DISX=XB-XA
DISY=YB-YA
DISH=0.001*(HB-HA)
DISI=SQRT(DISX**2+DISY**2+DI3H**2)
C DISI IS LENGTH(KM) OF LINE CONSIDERED.
IF (DISI) 420,420,430
420 CURR=0.
GO TO 1000
430 DDI=DISI*1000./20.
C ONE-HALF IS INCLUDED IN THE 20.
C DDI IS ONE-HALF TIMES 1/10 OF DISI (M)
DX=DISX/10.
DY=DISY/10.
DH=(HB-HA)/10.
C DH IS IN METERS.
PREV=0.
KNTRL=1
XI=XA
01600
01610
01620
01630
01640
01650
01660
01670
01680
01690
01700
01710
01720
01730
01740
01750
01760
01770
01780
01790
01800
01310
01820
01830
01840
01850
01860
01870
01880
01890
01900
01910
01920
01930
01940
01950
01960
01970
01980
01990
02000
0201C
0202C
0203C
0204C
0205C
0206C
02071
02081
0209(
02101
0211
0212
122
-------
SPCLN
YI=YA
HI=HA
KNT=0
DO 570 1=1,11
C STORE EACH XI,YI.
XST(I)=XI
YST(I)=YI
HST(I)=HI
IF (XST(I)) 440,440,450
440 RC=0.
GO TO 480
450 XZ=XI+XVZL
XY=XI+XVYL
H=HI
GO TO (470,460) , IUZS
460 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR
IF (UZ.LT.1.0) UZ=1.0
470 CALL RCONCP (Z,H,HL,XZ,XY,YI,KST,AN,M,SY,SZ,RC)
RC=RC/UZ
HH=H*0.001
DIS=DISTA(XI,YI,HH)/DISTOL
SPEED=SQRT(VI+DELSPD*DIS)
IF (SPEED.LT.SCALV) SPEED=SCALV
RC=RC/SPEED
4dO GO TO (4yO,530), KNTRL
C IF RC IS ZERO, CONTINUE UNTIL RC IS POSITIVE.
490 IF (RC) 560,560,500
500 IF (1-1) 510,510,520
510 KNTRL=2
GO TO 550
C RESET POINT A TO LAST ONE PREVIOUS.
520 XA=XST(I-1)
YA=YST(I-1)
HA=HST(I-1)
KNTRL=2
GO TO 550
530 IF (RC) 540,540,550
C RESET POINT B IF REACH ZERO CONCENTRATION.
540 XB=XI
YB=YI
HB=HI
GO TO 580
550 KNT=KNT+1
560 XI=XI+DX
YI=YI+DY
H1=HI+DH
570 CONTINUE
580 IF (KNT) 600,600,5yO
590 IF (KNT-6) 410,410,640
C IF GET TO 370, CONC. FROM THIS SEGMENT IS 0.
600 GO TO (610,620,b30), INDIC
610 RC=0.
02130
02140
02150
02160
02170
02180
02190
02200
02210
02220
02230
02240
02250
02260
02270
02280
02290
02300
02310
02320
02330
02340
02350
02360
02370
023oO
02390
02400
02410
02420
0243U
02440
02450
02460
02470
024oO
02490
02bOO
02510
02520
02530
02540
02550
02560
02570
025bO
0259U
U2600
U2610
02620
02630
02o40
02o50
123
-------
SPCLN
GO TO 1060
620 FIRST=0.
GO TO 1030
630 RC=FIRST
GO TO 1050
640 CONTINUE
C DO ANOTHER TRAPEZOIDAL INTEGRATION FROM A TO 3 IN TEN STEPS.
C IT IS LIKELY THAT A TO B HAVE BEEN REDEFINED.
DISX=X8-XA
DISY=YB-YA
DISH=0.001*(HB-HA)
DISI=SQRT(DI3X**2+DISY**2+DISH**2)
C DI3I IS DISTANCE(KM) FROM A TO 6
DELD=DISI*100.
C DELD IS 1/10 DISI IN METERS.
DX=DISX/10.
DY=DISY/10.
DH=(HB-HA)/1U.
SUM=0.
XDUM=XA
YDUM=YA
IF (XDUM) 650,650,660
b50 RC=0.
GO TO 690
660 XZ=XDUM+XVZL
XY=XDUM+XVYL
H=HA
GO TO (680,670), IUZS
b70 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR
IF (UZ.LT.1.0) UZ=1.0
680 CALL RCONCP (Z,H,HL,XZ,XY,YDUM,KST,AN,M,SY,SZ,RC)
RC=RC/UZ
690 CONTINUE
HH=H*0.001
DIS=DISTA(XDUM,YDUM,HH)/DISTOL
SPEED=SQRT(VI+DELSPD*DIS)
IF (SPEED.LT.SCALV) SPEED=SCALV
RC=RC/SPEED
SUM=SUM+RC/2.
DO 750 1=1,9
XDUM=XDUM+DX
YDUM=YDUM+DY
H=H+DH
IF (XDUM) 700,700,710
700 RC=0.
GO TO 740
710 GO TO (730,720), IUZS
720 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR
IF (UZ.LT.1.0) UZ=1.0
730 XZ=XDUM+XVZL
XY=XDUM+XVYL
02660
02670
02660
02690
02700
02710
02720
02730
02740
02750
02760
02770
02780
02790
02800
02810
02820
02830
02840
02850
02860
02870
02880
02890
02900
02910
02920
02930
02940
02950
02960
02970
02980
02990
03000
03010
03020
03030
03040
03050
03060
03070
03080
03090
03100
03110
03120
03130
03140
03150
03160
03170
03180
124
-------
SPCLN
CALL RCONCP U,H,HL,XZ,XY,YDUM,KST,AN,i«l,SY,SZ,RC)
RC=RC/UZ
740 CONTINUE
HH=H*0.001
DIS=DISTA(XDUM,YDUM,HH)/DISTOL
SPEED=SQRT(VI+DELSPD*DIS)
IF (SPEED.LT.SCALV) SPEED=SCALV
RC=RC/SPEED
750 SUM=SUM+RC
XDUM=XDUM+DX
YDUM=YDUM+DY
H=H+DH
IF (XDUM) 760,760,770
760 RC=0.
GO TO bOO
770 GO TO (790,7bO), IUZS
780 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(200./UHGT)**PwR
IF (UZ.LT.1.0) UZ=1.0
7yO XZ=XDUM+XVZL
XY=XDUi*l+XVYL
CALL RCONCP (Z,H,HL,XZ,XY,YDIM ,KST,AN,M,SY,SZ,RC)
RC=RC/UZ
BOO CONTINUE
HH=H*0.001
DIS=DISTA(XDUM,YDUM,HH)/DI3TOL
SPEED=SQRT(VI+DELSPD*DIS)
IF (SPEED.LT.SCALV) SPEED=SCALV
RC=RC/SPEED
SUM=SUM+RC/2.
C INTEGRATED VALUE IS CURR.
CURR=SUh*DELD
ILIM=10
810 PREV=CURR
C EVALUATE FOR POINTS IN BETWEEN THOSE ALREADY EVALUATED.
DELD=DELD/2.
XDUM=XA+DX/2.
YDUM=YA+DY/2.
H=HA+DH/2.
GO TO (830,820), IUZS
820 UZ=U*(H/UHGT)**PWR
IF (H.GT.200.) UZ=U*(200./UHGT)**PVvR
IF (UZ.LT.1.0) UZ=1.0
830 DO 690 I=1,ILIM
IF (XDUM) 840,840,850
840 RC=0.
GO TO 860
850 XZ=XDUM+XVZL
XY=XDUM+XVYL
CALL RCONCP (Z,H,HL,XZ,XY,YDUM,KST,AN,rt,SY,SZ,RC)
RC=RC/UZ
360 CONTINUE
HH=H*0.001
03200
03210
03220
03230
03240
03250
03260
03270
03280
03290
03300
03310
03320
03330
03340
03350
03360
03370
03380
033^0
03400
03410
03420
03430
03440
03450
03460
03470
03480
03490
03500
03510
03520
03530
03540
03550
03560
03570
03580
03590
03600
03610
03620
03630
03640
03650
03660
03670
03680
03690
03700
03710
125
-------
SPCLN
DIS=DISTA(XDUMfYDUMrHH)/DISTOL 03720
SPEED=SQRT(VI+DELSPD*DIS) 03730
IF (SPEED.LT.SCALV) SPEED=SCALV 03740
RC=RC/SPEED 03750
C NOTE ADD THESE TO EC'S FOUND ABOVE. 03760
SUM=SUM+RC 03770
H=H+DH 03780
GO TO (680,870), IUZS 03790
870 UZ=U*(H/UHGT)**PfcR 03800
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR 03810
IF (UZ.LT.1.0) UZ=1.0 03820
880 XDUM=XDUM+DX 03830
890 YDUM=YDUM+DY 03840
CURR=SUM*DELD 03850
TEST=ABS((CURR-PREV)/CURR) 03860
C IF WITHIN PIN OF LAST VALUE (PREV), CONSIDER THIS AS FINAL VALU 03870
IF (TEST-PIN) 1000,900,900 03880
900 ILIM=ILIM*2 03890
PREV=CURR 03900
C EVALUATE POINTS IN BETWEEN. 03910
DELD=DELD/2. 03920
DX=DX/2. 03930
DY=DY/2. 03940
DH=DH/2. 03950
XDUM=XA+DX/2. 03960
YDUM=YA+DY/2. 03970
H=HA+DH/2. 03980
GO TO (920,910), IUZS 03990
910 UZ=U*(H/UHGT)**PWR 04000
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR 04010
IF (UZ.LT.1.0) UZ=1.0 04020
920 DO 980 I=1,ILIM 04030
IF (XDUM) 930,930,940 04040
930 RC=0. 04050
GO TO 950 04060
940 XZ=XDUM+XVZL 04070
XY=XDUM+XVYL 04080
CALL RCONCP (Z,H,HL,XZ,XY,YDUM,KST,AN,M,SY,SZ,RC) 04090
RC=RC/UZ 04100
950 CONTINUE 04110
HH=H*0.001 04120
DIS=DISTA(XDUM,YDUM,HH)/DISTOL 04130
SPEED=SQRT(VI+DELSPD*DIS) 04140
IF (SPEED.LT.SCALV) SPEED=SCALV 04150
RC=RC/SPEED 04160
SUM=SUM+RC 04170
H=H+DH 04180
GO TO (970,960), IUZS 04190
960 UZ=U*(H/UHGT)**PWR 04200
IF (H.GT.200.) UZ=U*(200./UHGT)**PWR 04210
IF (UZ.LT.1.0) UZ=1.0 04220
970 XDUM=XDUM+DX 04230
980 YDUM=YDUM+DY 04240
126
-------
SPCLN
CURR=SUM*DELD
TEST=ABS((CURR-PREV)/CURR)
IF (TEST-PIN) 1000,990,990
990 ILIM=ILIM*2
DX=DX/2.
DY=DY/2.
DH=DH/2.
GO TO 810
C AT 1000 HAVE FINAL VALUE OF INTEGRATION IN CURR.
1000 GO TO (1010,1020,1040)", INDIC
1010 RC=CURR
GO TO 1050
1020 FIRST=CURR
1030 INDIC=3
XA=XLOC
YA=0.
HA=HLOC
XB=X2
YB=Y2
HB=H2
GO TO 410
1040 RC=FIRST+CURR
1050 CONLS(NC,NW)=CONLS(NC , NW)+RC*QL*DVS(NW,
1060 CONTINUE
1070 CONTINUE
1060 CONTINUE
RETURN
C
END
04250
04260
04270
04280
04290
04300
04310
04320
04330
04340
04350
04360
04370
04380
04390
04400
04410
04420
04430
04440
04450
04460
04470
044BO
04490
04500
04510
04520
04530
127
-------
SPCCR
C
C
C
C
C
C
C
o
u
C
C
C
C
C
C
C
r^
\^
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
SUBROUTINE SPCCR (MHR,NO,NR,PIM)
SPCCR CALCULATES HOURLY CONCENTRATIONS FROM A
SINGLE LANE CURVED PATH SOURCE. VEHICLES ARE
ALLOWED TO HAVE A CONSTANT ACCELERATION ALONG
THE PATH SOURCE.
THE FOLLOWING SUBROUTINES ARE CALLED BY SPCCR.
CURLIN--CALCULATES THE INTERSECTION OF THE
LINE FORMED BY THE WlMD DIRECTION AND
RECEPTOR COORDINATES AwD THE CIRCLE
DETERMINED FROM POINTS A,B,C.
RCONCP—DETERMINES THE RELATIVE CONCENTRATION AT
A RECEPTOR FROM A POINT AT A GIVEN UPWIND
AND CROS3WIND DISTANCE.
THE FOLLOWING FUNCTIONS ARE CALLED BY SPCCR.
XVY--CALCULATES THE VIRTUAL DISTANCE NECESSARY
TO ACCOUNT FOR THE INITIAL CROSSWIND DISPERTION.
XVZ—CALCULATES THE VIRTUAL DISTANCE NECESSARY
TO ACCOUNT FOR THE INITIAL VERTICAL DISPERTION.
ANGARA-DETERMINES THE APPROPRIATE ARCTAw(M/N)
WITH THE RESULTING ANGLE BE.TWEEN 0 AND
360 DEGREES.
DIFANG—DETERMINES THE DIFFERENCE BETWEEN TWO ANGLES.
THE RESULTING ANGLE IS BETWEEN 0 AND 360 DEGREES.
COMMON /SOCS/ QLNA(3D ,HCL(3D , RBQA ( 3 1 ) ,SBQA( 3 1) ,RMQA ( 3 1) ,SMQA(3D
1 ,REQA(3D ,SEQA(3D ,SIYA(31) ,SIZA(3D ,CONCLA(31 ,24) ,DVHA(3O ,IUZE,S
2PEK3O ,SPEF(31) ,TVCL(31),VSCL(31)
COMMON /REC/ RR(31),SR(31),ZR(31)
COMMON /WEA/ WTHET(25),WU(25),MKST(25),WHL(25),WTA(25),UHGT
DIMENSION PUZ(6), RRH(11), SSS(11), DIR(11), ARC(5), XANG(5)
DATA PUZ /0.15,0.15,0.20,0.25,0.40,0.607
QLNA PATH SOURCE STRENGTH FOR EACH LANE
HCL HEIGHT OF THE PATH SOURCE
RBQA EAST COORDINATE,POINT A
SBQA NORTH COORDINATE,POINT A
RMQA EAST COORDINATE,POINT B
SMQA NORTH COORDINATE,POINT B
REQA EAST COORDINATE,POINT C
SEQA NORTH COORDINATE,POINT C
SIYA INITIAL SIGMA Y
SIZA INITIAL SIGMA Z
CONCLA CONCENTRATION
DVHA DIURNAL VARIATION,CURVED PATH SOURCE
IUZE WIND INCREASE WITH HEIGHT
SPEI INITIAL VEHICLE SPEED
SPEF FINAL VEHICLE SPEED
TVCL VEHICLE VOLUME
VSCL GROSS ESTIMATE OF VEHICLE SIZE
RR EAST COORDINATE OF RECEPTOR
SR NORTH COORDINATE OF RECEPTOR
ZR HEIGHT OF RECEPTOR ABOVE GROUND
WTHET WIND DIRECTION
WU WIND SPEED
MKST STABILITY CLASS
(G/SEC)
(M)
(KM)
(KM)
(KM)
(KM)
(KM)
(KM)
(M)
(M)
(G/M**3)
(DIMENSIONLESS)
(DIMENSIOMLESS)
(M/SEC)
(M/SEC)
(VEH/HR)
(M)
(KM)
(KM)
(M)
(DEG AZIMUTH)
(M/SEC)
(DIMENSIQNLESS)
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00130
00190
00200
00210
00220
00230
00240
00250
00260
00270
00280
00290
00300
00310
00320
00330
00340
00350
00360
00370
00380
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
128
-------
SPCCR
C WHL MIXING HEIGHT (M)
C WTA AMBIENT AIR SURFACE TEMP. (DEC
C Vi/P AMBIENT AIR SURFACE PRESSURE (MB)
X(R,S)=(R-RREC)*SINT+(S-S REC)*COST
C X IS UPWIND DISTANCE OF R,S FROM RREC,SREC
Y(R,S)=(3-SREC)*SINT-(R-RREC)*COSr
C Y IS CROSSWINO OISTA^CS OP R,S FROM RREC,SREC
C START 00 LOOP FOR SACH HOUR
DO 500 NW=1,NHR
THETA=WTHST(iM)
TR=THETA/57.295ti
SINT=SIN(IR)
COST-COS(TR)
U=VvU(NW)
UZ=U
KST=MKST(fW)
HL=WHL(NW)
PWR=PUZ{KST)
C START DO LOOP FOR EACH CURVED PATH
DO 500 NS=1,NQ
QL=QLNA(N3)*('IVCL(HS)/3600. )
SCALV=VSCL(NS)*(TVCL(WS)/3600.)
RB=RBQA(NS)
RM=RMQA(NS)
RE=REQA(NS)
SB=SBQA(NS)
SM=SMQA(NS)
SE=SEQA(NS)
VI=SPEI(NS)**2
VF=SPEF(NS)**2
DELSPD=VF-VI
H=HCL(NS)
GO TO (20,10), IUZE
10 UZ=U*(H/UHGT)**P^R
IF (H.GT.200.) UZ=U*(200./UHGT)**t>^R
IF (UZ.LT.1.0) UZ=1.0
20 IF (SIYA(NS)) 40,30,40
30 XVYL=0.
GO TO 50
40 SYON=SIYA(NS)
XVYL=XVY(SYON,KST)
50 IF (SIZA(NS)) 70,60,70
60 XVZL=0.
GO TO 80
70 SZON=SIZA(NS)
XVZL=XVZ(SZOr^,KST)
80 CONTINUE
C FROM HERE THROUGH STATSMENT 150 THE CENTER AND
C CIRCLE,THROUGH POINTS A,B AND C,ARE CALCULATED
RMID=(RB+RM)/2.
SMID=(SB+SM)/2.
RNID=(RM+RE)/2.
SNID=(SM+SE)/2.
RADIUS OF THE
00540
00550
00560
00570
00580
00590
00600
00610
00620
00630
00640
00650
00660
00670
00680
00690
00700
00710
00720
00730
00740
00750
00760
00770
00780
00790
00800
00810
00820
00830
00840
00850
00860
00870
00880
00890
00900
00910
00920
00930
00940
00950
00960
00970
00980
00990
01000
01010
01020
01030
01040
01050
01060
129
-------
SPCCR
90
100
110
120
130
C
C
C
140
150
C
C
C
IF (RM-R13.SQ.O.) GO TO 90
IF (RE-RM.EQ.O, ) GO TO 100
SL01=-( (RM-RB)/(3'4-SlJ) )
SL02=-( (RE-Rd)/(3E-SM) )
RC= ( SL01 *RMI D-K9N I D-SMI D-S L02 *RN I 0 ) / ( SL01-S 002
SC=SLOl*(RC-RrtIO)-K<3inID
GO TO 110
RC=RMID
SL02=-( (RE-RM)/(3B-SM) )
SC=SL02*(RC-RNID)-t-SNID
GO TO 110
RC=RNIO
SL01=-( (RM-RI3)/(3>l-Srj) )
SC=SL01*(RC-RC4IO)+3'<1IO
GO TO 110
RAD=SJ.RT( U3-RB)**2+(3C-S!3)**2)
DELM=RB-RC
DELN=SB-3C
ANGB=ANGARC ( OWfjrt , DELN )
DELM=RM-RC
DELN=SM-SC
ANGi«l=ANGARC ( Ofiij 4 , OfiM )
DELM=RE-RC
DEL_N=SE-SC
ANGE=ANGARC ( O
ANGBE=A3S ( 0 1 H1
L-4 , OELM )
^ 3 ( AMGB , ANGE ) )
ANGBM=A33 ( 01 PANG ( ANG3 , AMG'1 ) )
ANGME=ABS ( 0 CFA'.IG ( ANG«i , ANGS ) )
IF (ANGBM.LT.AMGBR.AMD.ANGME.ivr.ANJGBE)
GAMDEG=360.-AttGBE
GO TO 130
GO TO 120
CONTINUE
GAMDEG 13 THE ANGLR SUBTENDED B^
IF A=l THE ROTATION FROM A TO 6
IF A=-l THE ROTATION FROM A TO ',
A=l.
IF (GAMDEG.Gfi. 180. ) GO TO 140
BEDEG=DIFANG(ANGe,ANGB)
IF (3EDEG.LT.O.) A=-l.
GO TO 150
RDEG=ANG8+li3U.
'IMK AUC A.HC
IS CLOCKS [Si:
C3 COUNTER-CLOCKWISE
A=-l .
REDEG=DIFASG(ANGE,RO EG)
IF (RMDEG.GE.O..AND.REDEG.LE. 0.)
CONTINUE
DO 500 NC=1,NR
RREC=RR(NC)
SREC=SR(NC)
Z=ZR(NC)
CURLIN CALCULATES THE COORDINATES De1 aN\T UPWIND L.OCI ON ARC A3
FROM HERE TO .STA'CE^ENT 31i) [MS BEGINNING AND ENDING
POINTS,WITH TUcJtH CORRESPONDING OIREC PIONS , ^)R EACd SEGMENT
01070
01080
01090
01100
OHIO
01120
01130
01140
01150
01160
01170
01180
011*0
01200
01210
01220
01230
01240
01250
01260
01270
01280
01290
01300
01310
01320
01330
01340
01350
01360
01370
01380
01390
01400
01410
01420
01430
01440
01450
01460
01470
014bO
01490
01500
01510
01520
01530
01540
01550
01560
01570
01580
01590
130
-------
SPCCR
C OF THE ARC ARE CALCULATED
CALL CURLIN (RC,SC,RREC,SRBC,2,RAD,THETA,RL,SL,RLl,SL1,NLOCI)
IF (NLOCI) 160,170,180
C NLOCI HAS VALUES OF -1,0,OR 1 FOR:
C NO UPWIND LOCI OM ARC AS: 1 LOCI ON ARC AB:
C OR 2 LOCI OM ARC AB
160 XANG(l)=AflG8
ARC(1)=GAMDEG
. J=l
GO TO 310
170 K=l
GO TO 190
180 K=2
190 N=l
DO 270 KK=1,K
IF (KK.EQ.2) GO TO 200
DELM=RL-RC
DELN=SL-SC
ANGRL=ANGARC(OSLM,OELN)
RANG=ANGRL
GO TO 220
200 DELM=RL1-RC
210 IF (K.NE.l) GO TO 270
XANG(1)=ANGB
ARC(1)=GAMDEG
DELN=SL1-SC
ANGRL=ANGARC(DSLM,DELN)
RLANG=ANGRL
220 IF (GAi4DBG.GT.l80. . ANfD. A.GT. 0. ) GO TO 230
IF (GAMDEG.GR.180. .AND. A.LT. 0 .) GO TO 240
IF (GAMDEG.LT.180..AND.A.LT.0,) GO TO 250
BIN=DIFANG(ANGRL,ANGB)
IF (BIN.LS.O..OR.BIN.GE.GAMDEG) GO TO 210
ARC(N)=BIN
GO TO 260
230 BIN=DIFANG(ANGRL,ANGB)
IF (BIN.GT.O.) ARC(N)=BIN
IF (BIN.GT.O.) GO TO 260
IF (BIN.LE.O.) DIF=(GAMDEG-360.)-BIN
IF (DIF.LT.O.) GO TO 210
ARC(N)=BIN+360.
GO TO 260
240 BIN=DIFAtf3(AN3RL,ANGB)
IF (BIN.LT.O.) ARC(N)=8IN
IF (BIN.LT.O.) GO TO 260
IF (BIN.GE.O..AND.BIN.LE.360.-GAMDEG) GO TO 210
ARC(N)=360.-BIN
GO TO 260
250 BIN=DIFANG(ANGRLfANG3)
IF (BIN.GE.O..OR.BIN.LE.-GAMDEG) GO TO 210
ARC(N)=BIM
GO TO 260
260 N=N+1
01600
01610
01620
01630
01640
01650
01660
01670
01680
01690
01700
01710
01720
01730
01740
01750
01760
01770
01780
01790
01800
01810
01820
01830
01840
01850
01860
01870
01880
01890
01900
01910
01920
01930
01940
01950
01960
01970
01980
01990
02000
02010
02020
02030
02040
02050
02060
02070
02080
02090
02100
02110
02120
131
-------
SPCCR
INDEX=-1
IF (KK.EQ.2) INDEX=1
GO TO 270
270 CONTINUE
XANG(1)=ANGB
IF (N-l.EQ.O) GO TO 300
IF (N-1.EQ.2) GO TO 280
ARC ( 2 ) =GAMDEG-ABS ( ARC ( 1 ) )
XANG ( 2 ) =RANG
IF (INDEX. EQ. I) XA4G(2)=RLA43
J=2
GO TO 310
280 ARC(1)=AB3(ARC(1) )
ARC(2)=A3S(ARC(2) )
IF (ARC(D.GT.ARC(2) ) GO TO 290
ARC ( 2 ) =ARC ( 2 ) -ARC ( 1 )
ARC ( 3 ) =GAMOS3-ARC ( 1 ) - ARC ( 2 )
XANG(2)=RA*I3
XANG(3)=RLANG
J = 3
GO TO 310
2yO AA=ARC(1)
ARC(1)=ARC(2)
ARC ( 2 ) =AA-ARC ( 2 )
ARC ( 3 ) =G AMD EG -ARC ( 1 ) - ARC ( 2 )
XANG ( 2 ) =RLAtfo
XANG(3)=RAMG
J = 3
GO TO 310
300 J=l
ARC(l)=GA:J10f-r,
310 lOTAL-O.
DISTOT=RAO* ( 3
-------
SPCCR
340
3bO
360
C
C
370
C
380
C
C
390
C
400
GO TO 3Si!
CALL RCOMCP U,^,HL,XZ,XY,YI,KS'C,AN,rtl,SY,SZ,CHI)
IF (CHI) 370,370,360
NM=NM+1
NCHI=NCHI+1
THE NEXT TwO STEPS DETERMINE THc] DCKECTCON «'R.DM
THE CENTER TO THE MEW START ['•IG P-1CNT
IF (L.EQ.l) DIRI*JT=UIR(L)
IF (NM.EQ.l.ANO.L.GT.l) 01R [N't^-H «( ',-1 )
IF (L.GT.l.ANO.NM.EQ.l) N>1=2
LCOUNT=0
IF (L.NE.ll) LCXINI.^1
CONTINUE
NM=NM+LCOUNT
NM-1 IS T:iE NUMBER OP NON-ZERO '3ECT CONS Or' ['HE ARC SEGMENT
IF (NM) 470,470,380
IF (NCHI.GT.6) GO TO 390
VAR=NM-1
ARC(M)=ABS(ARC(M))*(7AR/10.)
XANG(M)=DIRINT
NM=0
NCHI=0
GO TO 320
DO A TRAPEZOIDAL IN UEGKA'L" £ON t-'ROM A TO 8 IN TEN STEPS.
IT IS LIKELY THAT A OR B HAVE BEEN REDEFINED.
11=11
L=l
KK=1
SUM=0.
VAR=NM-1
ARR=A3S(ARC(M))*(7AR/10.)
DELL=((ARR/10.)/57.29578)*RAD*1000.
DELL=1/10 OF THE ARC LENGTH IN METERS
DO 440 I=L,II,KK
410
420
430
DIRR=DIRINT+(A*(I-1)/DIV)*ARR
DIRRAD=DIRR/57.29b78
RRR1=RC+RAD*SIN(DIRRAD)
SSS1=SC+RAD*COS(DIRRAD)
XZ=X(RRR1,SSS1)+XVZL
XY=X(RRR1,SSS1)+XVYL
YI=Y(RRR1,SSS1)
DTHETA=DIFANG(DIRR,ANGB)
IF (A.LT.O..AND.DTHETA.GT.O.)
IF (A.GT.O..AND.DTHETA.LT.O.)
DIS=(RAD*(DTHETA/57.29578))/DISTOT
SPEED=SQRT(\/I+DELSPO*DIS)
IF (SPEED.LT.SCALV) SPEED=SCALV
IF (XZ-XVZL) 410,410,420
CHI=0.
GO TO 430
CALL RCONCP (Z,H,HL,XZ,XY,YI,KST,AN,MI,3Y,SZ,CHI)
CHI=CHI/SPEEO
OTBETA=3 60.-ABS(DTHETA)
DTHETA=360.-ABS(DTHETA)
02660
02670
02680
02690
02700
02710
02720
02730
02740
02750
02760
02770
02780
02790
02800
02810
02820
02830
02840
02850
02860
02870
02880
02890
02900
02910
02920
02930
02940
02950
02960
02970
02980
02990
03000
03010
03020
03030
03040
03050
03060
03070
03080
03090
03100
03110
03120
03130
03140
03150
03160
03170
03180
133
-------
SPCCR
440
450
C
C
C
C
C
460
470
460
490
500
C
IF (I.EQ.l) Ql=CHI/2.
IF (I.EO..TT.) Q2=C4I/2.
SUM=SUM+Cfl [
CONTINUE
KNT=KNT+i
IF (KNT.GT.l)
-------
XPLUME
SUBROUTINE XPLUME (HX,F,DELHF,HP,U,X)
C BRIGGS EFFECTIVE HEIGHT AT THE DISTANCE X.
C OUTPUT VARIABLES ARE...
C HX EFFECTIVE PLUME HEIGHT FOR DISTANCE X (METERS)
C INPUT VARIABLES ARE...
C F BUOYANCY FLUX (M**4/SEC**3)
C DELHF FINAL PLUME RISE (METERS)
C HP PHYSICAL STACK HEIGHT (METERS)
C U WIND SPEED (M/SEC)
C X DOWNWIND DISTANCE (KM)
XM=1000.*X
C XM IS X IN METERS.
C STATEMENT 10 IS EQUATION (2), BRIGGS(1971).
10 DELHX=1.6*F**0.333333*XM**0.666667/U
IF (DELHX.GT.DELHF) DELHX=DELHF
HX=HP+DELHX
RETURN
C
END
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00180
00190
135
-------
FPLUME
SUBROUTINE FPLUME (HF,DISTF,F,DELHF,HP,TS,VS,D,VF,KST,U,DTHDZ ,T) 00010
C BRIGGS FINAL HEIGHT 00020
C D. B. TURNER, ENVIRONMENTAL APPLICATIONS BRANCH 00030
C METEOROLOGY LABORATORY, ENVIRONMENTAL PROTECTION AGENCY 00040
C RESEARCH TRIANGLE PARK, N C 27711 00050
C (919) 549 - 8411, EXTENSION 4565 00060
IF (T) 10,10,20 00070
C T = 0. MEANS NO AMBIENT TEMPERATURE GIVEN. USE T = 293. 00080
10 T=293. 00090
C OUTPUT VARIABLES ARE... 00100
C HF FINAL EFFECTIVE PLUME HEIGHT (METERS) 00110
C F BUOYANCY FLUX (M**4/SEC**3) 00120
C DELHF FINAL PLUME RISE (METERS) 00130
C DISTF DISTANCE OF FINAL PLUME RISE FROM SOURCE (KM) 00140
C INPUT VARIABLES ARE... 00150
C HP PHYSICAL STACK HEIGHT (METERS) 00160
C TS STACK GAS TEMPERATURE (DEC K) 00170
C VS STACK GAS EXIT VELOCITY (M/SEC) 00180
C D INSIDE STACK DIAMETER (METERS) 00190
C VF STACK GAS VOLUMETRIC FLOW RATE (i«I**3/SEC) 00200
C KST STABILITY (CLASS), SEE PAGE 209 OF PASQUILL, 00210
C ATMOSPHERIC DISPERSION. CLASSES DEFINED BY... 00220
C 1 IS PASQUILL STABILITY CLASS A 00230
C 2 IS PASQUILL STABILITY CLASS B 00240
C 3 IS PASQUILL STABILITY CLASS C 00250
C 4 IS PASQUILL STABILITY CLASS D 00260
C 5 IS PASQUILL STABILITY CLASS E 00270
C b IS PASQUILL STABILITY CLASS F 00260
C U WIND SPEED (M/SEC) 00290
C DTHDZ POTENTIAL TEMPERATURE LAPSE RATE (DEG K/METER) 00300
C T AMBIENT AIR TEMPERATURE (DEG K) 00310
C IF VF IS NOT GIVEN, CALCULATE IT FROM STACK DATA. 00320
20 IF (VF) 30,30,40 00330
30 VF=0.785398*VS*D*D 00340
C THE CONSTANT 0.785398 = PI/4 00350
40 F=3.12139*VF*(TS-T)/TS 00360
C THE CONSTANT 3.12139 IS THE ACCELERATION DUE TO GRAVITY / PI. 00370
C GO TO APPROPRIATE BRANCH FOR STABILITY CONDITION GIVEN. 00380
C IF UNSTABLE OR NEUTRAL GO TO 50, IF STABLE GO TO 90. 00390
GO TO (50,50,50,50,90,90,90), KST 00400
C DETERMINE APPROPRIATE FORMULA FOR CALCULATING XST, DISTANCE AT 00410
C WHICH TURBULENCE BEGINS TO DOMINATE. THE FORMULA USED DEPENDS 00420
C UPON BUOYANCY FLUX. STATEMENTS 60 AND 70 ARE EQUATION (7), 00430
C BRIGGS(1971). 00440
50 IF (F-55.) 60,70,70 00450
60 XST=14.*F**.625 00460
GO TO 80 00470
70 XST=34.*F**.4 00480
80 DISTF=3.5*XST 00490
DELHF=1.6*F**0.333333*DISTF**0.666667/U 00500
GO TO 160 00510
90 IF (DTHDZ) 100,100,130 00520
C IF DTHDZ IS NEGATIVE OR ZERO ASSIGN TO IT A VALUE OF 0.02 OR 00530
136
-------
FPLUME
C 0.035 IF STABILITY IS SLIGHTLY STABLE OR STABLE, RESPECTIVELY. 00540
100 GO TO (50,50,50,50,110,120,120), KST 00550
110 DTHDZ=0.02 00560
GO TO 130 00570
120 DTHDZ=0.035 00580
130 S=9.8061b*DTHDZ/T 00590
C THE CONSTANT 9.80616 IS THE ACCELERATION DUE TO GRAVITY. 00600
C S IS A STABILITY PARAMETER. 00610
C CALCULATE PLUME RISE ACCORDING TO EQUATION (2), BRIGGS(1972). 00620
DHA=2.4*(F/(U*S))**0.333333 00630
C CALCULATE PLUME RISE BY EQUATION (5), BRIGGS(1971) FOR LIGHT 00640
C WIND CONDITIONS ACCORDING TO MORTON, TAYLOR, AND TURNER. 00650
DELHF=5.0*F**0.25/S**0.375 00660
IF (DHA-DELHF) 140,140,150 00670
140 DELHF=DHA 00680
C DISTANCE TO FINAL PLUME RISE IS GIVEN BY THE FOLLOWING 00690
150 DISTF=3.14159*U/S**0.5 00700
C CALCULATE FINAL EFFECTIVE HEIGHT. 00710
160 HF=HP+DELHF 00720
DISTF=DISTF/1000. 00730
RETURN 00740
C 00750
END 00760
137
-------
PGSIG
SUBROUTINE PGSIG (X,XY,KST,SY,SZ) 00010
C D. B. TURNER, ENVIRONMENTAL APPLICATIONS BRANCH 00020
C METEOROLOGY LABORATORY, ENVIRONMENTAL PROTECTION AGENCY 00030
C RESEARCH TRIANGLE PARK, N C 27711 00040
C (919) 549 - 8411, EXTENSION 4565 00050
C VERTICAL DISPERSION PARAMETER VALUE, SZ DETERMINED BY 00060
C SZ = A * X ** B WHERE A AND B ARE FUNCTIONS OF BOTH STABILITY 00070
C AND RANGE OF X. 00080
C HORIZONTAL DISPERSION PARAMETER VALUE, SY DETERMINED BY 00090
C LOGARITHMIC INTERPOLATION OF PLUME HALF-ANGLE ACCORDING TO 00100
C DISTANCE AND CALCULATION OF 1/2.15 TIMES HALF-ARC LENGTH. 00110
DIMENSION XA(7), XB(2), XD(5), XE(8), XF(9), AA(8), BA(8), AB(3), 00120
1BB(3), AD(6), BD(6), AE(9), BE(9), AF(10), BF(10) 00130
DATA XA /.5,.4,.3,.25,.2,.15,.I/ 00140
DATA XB /.4,.2/ 00150
DATA XD /30.,10.,3.,1.,.3/ 00160
DATA XE /40. ,20. ,10. ,4. ,2. ,1.,.3,.l/ 00170
DATA XF /60.,30. ,15.,7. ,3.,2.,1.,.7,.2/ 00180
DATA AA /453.85,346.75,258.89,217.41,179.52,170.22,158.08,122.8/ 00190
DATA BA /2.1166,1.7283,1.4094,1.2644,1.1262,1.0932,1.0542,.94477 00200
DATA AB /109.3j51d7md"51j7gx"/ 00210
DATA BB 71.0971,0.98332,0.931987 00220
DATA AD /44.053,36.650,33.504,32.093,32.093,34.459/ 00230
DATA BD /0.51179,0.56589,0.60486,0.64403,0.81066,0.86974/ 00240
DATA AE 747.618,35.420,26.970,24.703,22.534,21.628,21.628,23.331,2 00250
14.26/ 00260
DATA BE 70.29592,0.37615,0.46713,0.50527,0.57154,0.63077,0.75660,0 00270
1.81956,0.83667 00280
DATA AF /34.219,27.074,22..651,17.836,16.187,14.823,13.953,13.953,1 00290
14.457,15.2097 00300
DATA BF 70.21716,0.27436,0.32681,0.41507,0.46490,0.54503,0.63227,0 00310
1.68465,0.78407,0.815587 00320
GO TO (10,40,70,80,110,140), KST 00330
C STABILITY A (10) 00340
10 TH=(24.167-2.5334*ALOG(XY))/57.2958 00350
IF (X.GT.3.11) GO TO 170 00360
DO 20 ID=1,7 00370
IF (X.GE.XA(ID)) GO TO 30 00380
20 CONTINUE 00390
ID=8 00400
30 SZ=AA(ID)*X**BA(ID) 00410
GO TO 190 00420
C STABILITY B (40) 00430
40 TH=(18.333-1.8096*ALOG(XY))/57.2958 00440
IF (X.GT.35.) GO TO 170 00450
DO 50 ID=1,2 00460
IF (X.GE.XB(ID)) GO TO 60 00470
50 CONTINUE 00480
ID=3 00490
60 SZ=AB(ID)*X**B3(ID) 00500
GO TO 180 00510
C STABILITY C (70) 00520
70 TH=(12.5-1.0857*ALOG(XY))/57.2958 00530
138
-------
PGSIG
SZ=61.141*X**0.91465
GO TO 160
C STABILITY D (faO)
80 TH=(d.3333-0.72382*ALOG(XY))/57.2958
DO yo 10=1,5
IF (X.GE.XD(ID)) GO TO 100
yO CONTINUE
ID=6
100 SZ=AD(ID)*X**6D(ID)
GO TO 130
C STABILITY E (110)
110 TH=(6.25-0.54287*ALOG(XY))/57.2958
DO 120 ID=1,8
IF (X.GE.XE(ID)) GO TO 130
120 CONTINUE
ID=9
130 SZ=AE(ID)*X**BE(ID)
GO TO 180
C STABILITY F (140)
140 TH=(4.1667-0.36191*ALOG(XY))/57.2958
DO 150 ID=l,y
IF (X.GE.XF(ID)) GO TO 160
150 CONTINUE
ID=10
loO SZ=AF(ID)*X**BF(ID)
GO TO 180
17U SZ=5000.
GO TO 190
180 IF (SZ.GT.5000.) SZ=5000.
190 SY=465.11b*XY*SIN(TH)/COS(TH)
C 465.116 = 1000. (M/KM) / 2.15
RETURN
END
00540
00550
00560
00570
00580
00590
00600
00610
00620
00630
00640
00650
00660
00670
00680
00690
00700
00710
00720
00730
00740
00750
00760
00770
007aO
00790
00800
00810
00820
00830
00840
00850
00860
00870
139
-------
RCONCP
c
c
c
c
c
c
c
SUBROUTINE RCONCP (Z,H,HL,X,XY,Y,KST,AN,M,3Y,SZ,RC)
D. 6. TURNER, RESEARCH METEOROLOGIST* MODEL DEVELOPMENT
DIVISION OF METEOROLOGY, ENVIRONMENTAL PROTECTION
ROOM 314B, NCriS
MAILING ADDRESS-
' ON
BUILDING, RTF. PHONE (919) 549-8411
DM, EPA, RESEARCH TRIANGLE PARK, NC
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
10
20
30
40
c
c
50
c
60
C
C
C
c
ASSIGNMENT FROM NATIONAL OCEANIC AND ATMOSPHERIC
ADMINISTRATION, DEPARTMENT OF COMMERCE.
SUBROUTINE RCONCP CALCULATES CHI*U/Q CONCENTRATION VALU
RCONCP CALLS SUBROUTINE PGSIG TO OBTAIN STANDARD DEV
THE INPUT VARIABLES ARE....
Z RECEPTOR HEIGHT (M)
H EFFECTIVE STACK HEIGHT (M)
HL=L HEIGHT OF LIMITING LID (M)
X DISTANCE RECEPTOR IS DOWNWIND OF SOURCE (KM)
XY X+VIRTUAL DISTANCE USED FOR AREA SOURCE APPROX.
Y DISTANCE RECEPTOR IS CROSSvJIND FROM SOURCE (KM)
KST STABILITY CLASS
THE OUTPUT VARIABLES ARE
AU THE NUMBER OF TIMES THE SUMMATION TERM IS E-VALUAT
AND ADDED IN.
RC RELATIVE CONCENTRATION (SEC/M**3)
THE FOLLOWING EQUATION IS SOLVED —
RC = (1/(2*PI*SIGMA Y*SIGMA Z))* (EXP(-0.5*(Y/SIGMA
(EXP(-0.5*((Z-H)/SIGMA Z)**2) + EXF(-0.5*((Z+H)/SIGMA
PLUS THE SUM OF THE FOLLOWING 4 TERMS K TIMES (N=1
TERM 1- EXP(-0.5*((Z-H-2NL)/SIGMA Z)**2)
TERM 2- EXP(-0.5*((Z+H-2NL)/SIGMA Z)**2)
TERM 3- SXP(-0.5*((Z-H+2NL)/SIGMA Z)**2)
TERM 4- EXP(-0.5*((Z+H+2NL)/SIGMA Z)**2)
THE ABOVE EQUATION IS SIMILAR TO EQUATION (5-8) P 36 IN
WORKBOOK OF ATMOSPHERIC DISPERSION ESTIMATES WITH THE
OF THE EXPONENTIAL INVOLVING Y,TURNER(1970).
IWRI IS CONTROL CODE FOR OUTPUT
IWRI=6
IF THE SOURCE IS ABOVE THE LID, SET RC = 0., AND RETURN
IF(KST.GE.5) GO TO 50
IF (H-HL) 10,10,20
IF (Z-HL) 50,50,40
IF (Z-HL) 40,30,30
WRITE (IWRI,460)
RC = 0.
RETURN
IF X IS LESS THAN 1 METER, SET RC=0. AND RETURN. THIS
PROBLEMS OF INCORRECT VALUES NEAR THE SOURCE.
IF (X-0.001) 40,60,60
CALL PGSIG TO OBTAIN VALUES FOR SY AND SZ
CALL PGSIG (X,XY,KST,SY,SZ)
SY = SIGMA Y, THE STANDARD DEVIATION OF CONCENTRATION
Y-DIRECTION (M)
SZ = SIGMA Z, THE STANDARD DEVIATION OF CONCENTRATION
Z-DIRECTION (M)
C1 = 1 .
IF (Y) 70,90,70
1 BRANCH,
CY.
EXT 4564
27711
S3
TION3.
KM)
'ED
)**2))
;A Z)**2)
,K) --
ADDITION
.
AVOIDS
IN THE
IN THE
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00130
00190
00200
00210
00220
00230
00240
00250
00250
00270
00280
00290
00300
00310
00320
00330
00340
00350
00360
00370
00380
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
140
-------
RCONCP
70 YD=1000.*Y 00540
C YD IS CROSSWIND DISTANCE IN METEHS. 00550
DUM=YD/3Y 00560
TEMP=0.5*DUM*DUM 00570
IF (TEMP-50.) 80,40,40 00580
80 C1=EXP(TSMP) 00590
90 IF (KST-4) 100,100,110 00600
100 IF (HL-5000.) 190,110,110 00610
C IF STABLE CONDITION OR UNLIMITED MIXING HEIGHT, 00620
C USE EQUATION 3.2 IF Z = 0, OR EQ 3.1 FOR NON-ZERO Z, 00630
C TURNER(1970). 00640
110 C2=2.*SZ*SZ 00650
IF (Z) 40,120,140 00660
120 C3=H*H/C2 00670
IF (C3-50.) 130,40,40 00680
130 A2=1./EXP(C3) 00690
C WADE EQUATION 3.2,TURNER(1970). 00700
RC=A2/(3.14159*SY*SZ*C1) 00710
M=1 00720
RETURN 00730
140 A2=0. 00740
A3=0. 00750
CA=Z-H 00760
CB=Z+H 00770
C3=CA*CA/C2 00730
C4=CB*CB/C2 00790
IF (C3-50.) 150,160,160 00800
150 A2=1./EXP(C3) 00810
160 IF (C4-50.) 170,180,180 00820
170 A3=1./EXP(C4) 00830
C WADE EQUATION 3 • 1 ,TURNEH( 1970) . 00840
180 RC=(A24.A3)/(6.2&313*SY*SZ*C1) 00850
M=2 00860
RETURN 00870
C IF SIGMA-Z IS GREATER THAN 1.6 TIMES THE MIXING HEIGHT, 00880
C THE DISTRIBUTION BELOW THE MIXING HEIGHT IS UNIFORM WITH 00890
C HEIGHT REGARDLESS OF SOURCE HEIGHT. 00900
190 IF (SZ/HL-1.6) 210,210,200 00910
C WADE EQUATION 3.5,TURMER(1970). 00920
200 RC=1./(2.5066*SY*HL*C1) 00930
M=3 00940
RETURN 00950
C INITIAL VALUE OF AN SET = 0. 00960
210 AN=0. 00970
IF (Z) 40,370,220 00980
C STATEMENTS 220 TO 360 CALCULATE RC, THE RELATIVE CONCENTRATION, 00990
C USING THE EQUATION DISCUSSED ABOVE. SEVERAL INTERMEDIATE 01000
C VARIABLES ARE USED TO AVOID REPEATING CALCULATIONS. 01010
C CHECKS ARE HADE TO BE SURE THAT THE ARGUMENT OF THE 01020
C EXPONENTIAL FUNCTION IS NEVER GREATER THAN 50 (OR LESS THAN 01030
C -50). IF 'AN' BECOMES GREATER THAN 45, A LINE OF OUTPUT IS 01040
C PRINTED INFORMING OF THIS. 01050
C CALCULATE MULTIPLE EDDY REFLECTIONS FOR RECEPTOR HEIGHT Z. 01060
141
-------
RCONCP
220 A1=1./(6.28318*SY*SZ*C1) 01070
C2=2.*SZ*SZ 01080
A2=0. 01090
A3=0. 01100
CA=Z-h 01110
CB=Z+H 01120
C3=CA*CA/C2 01130
C4=CB*CB/C2 01140
IF (C3-50.) 230,240,240 01150
230 A2=1./£XP(C3) 01160
240 IF (04-50.) 250,260,260 01170
250 A3=1./EXP(C4) 01180
260 SUM=0. 01190
THL=2.*HL 01200
270 Ari=Atf+1. 01210
A4 = 0. 01220
A5=0. 01230
A6=0. 01240
A7=0. 01250
C5=AW*THL 01260
CC=CA-C5 01270
CD=CB-C5 01280
CE=CA+C5 01290
CF=CB+C5 01300
C6=CC*CC/C2 01310
C7=CD*CD/C2 01320
C8=CE*CE/C2 01330
C9=CF*CF/C2 01340
IF (C6-50.) 280,290,290 01350
280 A4=1./EXP(C6) 01360
290 I? (C7-50.) 300,310,310 01370
300 A5=1./EXP(C7) 01380
310 IF (C8-50.) 320,330,330 01390
320 A6=1./EXP(C8) 01400
330 IF (C9-50.) 340,350,350 01410
340 A7=1./EXP(C9) 01420
350 T=A4+A5+A6+A7 01430
SUM=SUM+T 01440
IF (T-0.01) 360,270,270 01450
360 ftC=A1*(A2+A3+SUM) 01460
tf=5 01470
RETUhN 01480
C CALCULATE MULTIPLE EDDY REFLECTIONS FOR GROUND LEVEL RECEPTOR H 01490
370 A1=1./(6.28318*SY*SZ*C1) 01500
A2=0. 01510
C2=2.*SZ*SZ 01520
C3=H*H/C2 01530
IF (C3-50.) 380,390,390 01540
380 A2=2./EXP(C3) 01550
390 SUH=0. 01560
THL=2.*HL 01570
400 AN=AN+1. 01580
A4=0. 01590
142
-------
RCONCP
A6=0. 01600
C5=AN*THL 01610
CC=H-C5 01620
CE=H+C5 01630
C6=CC*CC/C2 01640
C8=CE*CE/C2 01650
IF (C6-50.) 410,420,420 01660
410 A4=2./EXP(C6) 01670
420 IF (C8-50.) 430,440,440 01680
430 A6=2./EXP(C8) 01690
440 T=A4+A6 01700
SUM=SUM+T 01710
IF (T-0.01) 450,400,400 01720
450 RC=A1*(A2+SUM) 01730
M=4 01740
RETURN 01750
C 01760
460 FORMAT (1HO,'BOTH H AND Z ARE ABOVE THE MIXING HEIGHT SO A RELIABL 01770
1E COMPUTATION CAN NOT BE MADE.') 01780
C 01790
END 01800
143
-------
RCONCA
C
C
C
C
C
C
C
C
C
C
o
u
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
10
20
30
C
C
40
50
60
SUBROUTINE RCONCA (Z,H,HL,X,KST,AN,M,SZ,RCZ)
SUBROUTINE RCONCA CALCULATES CHI*U/Q, RELATIVE CONCENTRATION
NORMALIZED FOR WIND SPEED FOR A RECEPTOR DOWNWIND OF A
CROSSWIND INFINATE LINE SOURCE.
D. B. TURNER, RESEARCH METEOROLOGIST* MODEL APPLICATIONS BRANCH
METEOROLOGY LABORATORY, ENVIRONMENTAL PROTECTION AGENCY.
ROOM 316B, NCHS BUILDING, RTP. PHONE (919) 549-8411 EXT 4564
MAILING ADDRESS: MTL.EPA, RESEARCH TRIANGLE PARK, NC 27711.
* ON ASSIGNMENT FROM NATIONAL OCEANIC AND ATMOSPHERIC
ADMINISTRATION, DEPARTMENT OF COMMERCE.
THE INPUT VARIABLES ARE
Z RECEPTOR HEIGHT (M)
H EFFECTIVE STACK HEIGHT (M)
HL= HEIGHT OF LIMITING LID (M)
X DISTANCE RECEPTOR IS DOWNWIND OF SOURCE (KM)
KST STABILITY CLASS
SZ = SIGMA Z, THE STANDARD DEVIATION OF CONCENTRATION IN THE
Z-DIRECTION (M)
THE OUTPUT VARIABLES ARE
AN THE NUMBER OF TIMES THE SUMMATION TERM IS EVALUATED
AND ADDED IN.
RCZ RELATIVE CONCENTRATION (DIMENSIONLESS)
IWRI 13 CONTROL CODE FOR OUTPUT
THE FOLLOWING EQUATION IS SOLVED --
RC = (1/2.5066 *SIGMA Z))* (EXP(-0.5*(Y/SIGMA Y)**2)) *
(EXP(-0.5*((Z-H)/SIGMA Z)**2) + EXP(-0.5*((Z+H)/SIGMA Z)**2)
PLUS THE SUM OF THE FOLLOWING 4 TERMS K TIMES (N=1,K) --
EXP(-0.5*((Z-H-2NL)/SIGMA Z)**2)
EXP(-0.5*((Z+H-2NL)/SIGMA Z)**2)
EXP(-0.5*((Z-H+2NL)/SIGMA Z)**2)
EXP(-0.5*((Z+ri+2ML)/SIGMA Z)**2)
SQUARE HOOT OF 2 * PI
TO 330 CALCULATE RC, THE RELATIVE
TERM 1-
TERM 2-
TEfiM 3-
TERM 4-
2.5066 IS THE
STATEMENTS 190
CONCENTRATION,
USING THE EQUATION DISCUSSED ABOVE. SEVERAL INTERMEDIATE
VARIABLES ARE USED TO AVOID REPEATING CALCULATIONS.
CHECKS ARE MADE TO BE SURE THAT THE ARGUMENT OF THE
EXPONENTIAL FUNCTION IS NEVER GREATER THAN 50 (OR LESS THAN
-50). IF 'AN' BECOMES GREATER THAN 45, A LINE OF OUTPUT IS
PRINTED INFORMING OF THIS.
CALCULATE MULTIPLE EDDY REFLECTIONS FOR RECEPTOR HEIGHT Z.
SOURCE IS ABOVE
S) GO TO 40
10,10,20
40,40,180
180,30,30
IF THE SOURCE IS ABOVE THE LID, SET HC = 0
IF(KST.GE
IF (ri-HL)
IF (Z-HL)
IF (Z-HL)
WRITE (IivRI,430)
RETURN
IF X IS LESS THAN 1 METER, SET HC=0. AND
PROBLEMS OF INCORRECT VALUES NEAR THE SOURCE
(X-0.001) 160,50,50
(KST-4) 60,60,70
(HL-5000.) 150,70,70
AND RETURN
THIS AVOIDS
IF
IF
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00180
00190
00200
00210
00220
00230
00240
00250
00260
00270
00280
00290
00300
00310
00320
00330
00340
00350
00360
00370
00330
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
144
-------
RCONCA
C IF STABLE CONDITION OR UNLIMITED MIXING HEIGHT, 005*10
C USE EQUATION 3-2 IF Z = 0, OR EQ 3-1 FOR NON-ZERO Z, 00550
C TURN-ER(1970) . 00560
70 C2=2.*SZ*SZ 00570
IF (Z) 180,80,100 00580
80 C3=H*H/C2 00590
IF (C3-50.) 90,180,180 00600
90 A2=2./EXP(C3) 00610
C WADE EQUATION 3 .2 ,TURNER<1970). 00620
RCZ=A2/(2.5066*SZ) 00630
M=1 00640
RETURN 00650
100 A2=0. 00660
A3=0. 00670
CA=Z-H 00680
CB=Z+H 00690
C3=CA*CA/C2 00700
C4=CB*CB/C2 00710
IF (C3-50.) 110,120,120 00720
110 A2=1./EXP(C3) 00730
120 IF (C4-50.) 130,140,140 00740
130 A3=1./EXP(C4) 00750
C WADE EQUATION 3.1,TURNER(1970). 00750
140 RCZ=(A2+A3)/(2.5066*SZ) 00770
M=2 00780
RETURN 00790
C IF SIGMA-Z IS GREATER THAN 1.6 TIMES THE MIKING HEIGHT, 00800
C THE DISTRIBUTION BELOW THE MIXING HEIGHT IS UNIFORM WITH 00810
C HEIGHT REGARDLESS OF SOURCE HEIGHT. 00820
150 IF (SZ/HL-1.6) 170,170,160 00830
C WADE EQUATION 3.5,TURMER(1970). 00840
160 RCZ=1,/HL 00850
M=3 00860
RETURN 00870
C INITIAL VALUE OF AN SET = 0. 00880
170 AN = 0. 00690
IF (Z) 180,340,190 00900
1BO RCZ=0. 00910
RETURN 00920
190 A1=1./(2.506&*SZ) 00930
C2=2.*SZ*3Z 00940
A2=0. 00950
A3=0. 00960
CA=Z-H 00970
CB=Z+H 00980
C3=CA*CA/C2 00990
C4 = CB*Ct3/C2 01000
IF (C3-50.) 200,210,210 Q1Q10
200 A2=1./EXP(C3) 01020
210 IF (C4-5U.) 220,230,230 01030
220 A3=1./EXP(C4) 01Q40
230 3UM=0. 01050
THL=2.*HL 01060
145
-------
RCONCA
240
250
260
270
260
290
300
310
320
330
C
340
350
360
370
380
390
400
410
420
A4 = 0.
A5=0.
A6 = 0.
A7=0.
C5=AM*THL
CC=CA-C5
CD=CB-C5
CE=CA+C5
CF=CB+C5
C6=CC*CC/C2
C7=CD*CD/C2
C3=CE*CE/C2
C9=CF*CF/C2
IF (C6-50.) 250,260,260
A4=1 ./EXP(C6)
IF (C7-50.) 270,280,280
A5=1 ./EXP(C7)
IF (C8-50.) 290,300,300
A6=1 ./EXP(C8)
If (C9-50.) 310,320,320
A7=1 ./EXP(C9)
T=A4+A5+A6+A7
SUM=SUH+T
IF (T-0.01) 330,240,240
M=5
RETUKN
CALCULATE MULTIPLE EDDY REFLECTIONS FOR GROUND LEVEL RECEPTOR H
A1=1 ./(2.5066*SZ)
A2=0.
C2=2.*SZ*SZ
C3=H*H/C2
IF (C3-50.) 350,360,360
A2=2./EXP(C3)
SUM=0.
THL=2.*HL
A4=0.
A6=0.
C5=AN*THL
CC=H-C5
CE=H+C5
C6=CC*CC/C2
C8=CE*CE/C2
IF (C6-50.) 380,390,390
A4=2./EXP(C6)
IF (C8-50.) 400,410,410
A6=2./EXP(C8)
T=A4+A6
SUM=SUM+T
IF (T-0.01) 420,370,370
RCZ=A1*(A2+SUM)
01070
01080
01090
01100
01110
01120
01130
01UO
01150
01160
01170
01180
01190
01200
01210
01220
01230
01240
01250
01260
01270
01280
01290
01300
01310
01320
01330
01340
01350
013&0
01370
01330
01390
01400
01410
01420
01430
01440
01450
01460
01470
01480
01490
01500
01510
01520
01530
01540
01550
01560
01570
01580
01590
146
-------
RCONCA
M=4 01600
RETURN 01610
C 01620
430 FORMAT (1HO,'BOTH H AND Z ARE ABOVE THE MIXING HEIGHT SO A RELIABL 01630
1E COMPUTATION CAN NOT BE MADE.') 01640
C 01650
END 01660
147
-------
XVY
FUNCTION XVY (SYO,KST)
C XVY CALCULATES THE VIRTUAL DISTANCE NECESSARY TO
C ACCOUNT FOR THE INITIAL CROSSWIND DISPERTION.
GO TO (10,20,30,40,50,60), KST
10 XVY=(SYO/213.)**1.1148
RETURN
20 XVI=(SYO/155.)**1.097
RETURN
30 XVY=(SYO/103.)**1.092
RETURN
40 XVY=(SYO/68.)**1.07fa
RETURN
50 XVY=(SYO/50.)**1.086
RETURN
60 XVY=(SYO/33.5)**1.083
RETURN
C
END
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
00180
148
-------
XV Z
FUNCTION XVZ (SZO,KST) 00010
C XVZ CALCULATES THE VIRTUAL DISTANCE NECESSARY 00020
C TO ACCOUNT FOR THE INITIAL VERTICAL DISPERTION. 00030
DIMENSION SA(7), SB(2), SD(5), SE(8), SF(9), AA(d), AB(3), AD(6), 00040
1AE(9), AF(10), CA(8), CP(3)r CD(6), CE(9), CF(10) 00050
DATA SA 713.95,21.40,29.3,37.67,47.44,71.16,104.657 00060
DATA SB /20.23,40./ 00070
DATA SD /12.09,32.09,65.12,134.9,251.2/ 00080
DATA SE 73.534,8.698,21.626,33.469,49.767,79.07,109.3,141.867 00090
DATA SF 74.093,10.93,13.953,21.627,26.976,40.,54.89,68.84,83.257 00100
DATA AA /122.8,158.08,170.22,179.52,217.41,258.89,346.75,453.85/ 00110
DATA AB /90.673,98.483,109.3/ 00120
DATA AD /34.459,32.093,32.093,33.504,3b.650,44.053/ 00130
DATA AE 724.26,23.331,21.628,21.628,22.534,24.703,26.97,35.42,47.6 00140
118/ 00150
DATA AF 715.209,14.457,13.953,13.953,14.823,16.187,17.836,22.651,2 00160
17.074,34.2197 00170
DATA CA /1.0585,.9486,.9147,.8879,.7909,.7095,.5786,.4725/ 00180
DATA CB /I.073,1.017,.91157 00190
DATA CD 71.1498,1.2336,1.5527,1.6533,1.7671,1.95397 00200
DATA CE 71.1953,1.2202,1.3217,1.5854,1.7497,1.9791,2.1407,2.6585,3 00210
1.3793/ 00220
DATA CF /I.2261,1.2754,1.4606,1.5816,1.8348,2.151,2.4092,3.0599,3. 00230
16448,4.60497 00240
GO TO (10,40,70,80,110,140), KST 00250
C STABILITY A(10) 00260
10 DO 20 ID=1,7 00270
IF (SZO.LE.SA(ID)) GO TO 30 00280
20 CONTINUE 00290
10=8 00300
30 XVZ=(SZO/AA(ID))**CA(ID) 00310
RETURN 00320
C STABILITY B (40) 00330
40 DO 50 ID=1,2 00340
IF (SZO.LE.SB(ID)) GO TO 60 00350
50 CONTINUE 00360
ID=3 00370
60 XVZ=(SZO/AB(ID))**CB(ID) 00380
RETURN 00390
C STABILITY C (70) 00400
70 XVZ=(SZO/61.141)**!.0933 00410
RETURN 00420
C STABILITY D (80) 00430
80 DO 90 ID=1,5 00440
IF (SZO.LE.SD(ID)) GO TO 100 00450
90 CONTINUE 00460
ID=6 00470
100 XVZ=(SZO/AD(ID))**CD(ID) 00480
RETURN 00490
C STABILITY E (110) 00500
110 DO 120 ID=1,8 00510
IF (SZO.LE.SE(ID)) GO TO 130 00520
120 CONTINUE 00530
149
-------
xvz
ID=9 00540
130 XVZ=(SZO/AE(ID))**CE(ID) 00550
RETURN 00560
C STABILITY F(140) 00570
140 DO 150 ID=1,9 00580
IF (SZO.LE.SF(ID)) GO TO 160 00590
150 CONTINUE 00600
ID=10 00610
160 XVZ=(SZO/AF(ID))**CF(ID) 00620
RETURN 00630
C 00640
END 00650
150
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
10
20
30
c
c
c
c
CURLIN
SUBROUTINE CURLIN (RC,SC,R,S,Z,RAD,THETA,RL,SL,RLl,SL1,NLOCI)
THE INPUT DATA FOR CURLIN ARE THE RECEPTOR AND CENTER
COORDINATES,THE RADIUS OF THE CIRCLE AND THE WIND DIRECTION.
THE FUNCTION OF CURLIN IS TO DETERMINE WHERE THE LINE FORMED BY THE
DIRECTION AND RECEPTOR COORDINATES INTERSECTS THE CIRCLE.
R,S ARE RECEPTOR COORDINATES
RC,SC ARE CENTER COORDINATES
RAD 13 THE RADIUS
THETA IS THE WIND DIRECTION
RL,SL AND RL1,SL1 ARE LOCI ON THE CIRCLE.
RA=3.141b92bS4/160.
IF (THETA.EQ.O..OR.THETA.EQ.130.) GO TO 10
IF THE WIND DIRECTION IS 0 OR loO DEGREES THAN THE RECEPTOR
COORDINATE CALCULATIONS FOLLOW.IF THE WIND IS FROM ANY OTHER
DIRECTION THAN THE RECEPTOR COORDINATE CALCULATIONS BEGIN AT
STATEMENT 10.
GO TO 20
RL=R
RL1=R
A=l.
B=-2.*SC
C=SC**2-RAD**2+(R-RC)**2
DISC=B**2-4.*A*C
IF (DISC.LT.O.) GO TO 50
SL=(-B+SQRT(DISC))/(2.*A)
SLl=(-B-SQRT(DISC))/(2.*A)
GO TO 30
IF (THETA.GT.O..AND.THETA.LE.90.) SLOPE=1 ./TAN (THETA*RA)
IF (THETA.GT.90..AND.THETA.LT.160.) SLOPE=-1./TAN(IdO.*RA-RA*THETA
1)
IF (THETA.GT.180..AND.THETA.LE.270.) SLOPE=1./TAN(THETA*RA-RA*lfaO.
1)
IF (THETA.GT.270..AND.THETA.LT.360.) SLOPE = -TAN(THETA*RA-RA*270.)
DEL=-SLOPE*R+S-SC
A=1.+SLOPE**2
B=-2.*RC+2.*SLOPE*DEL
C=-(RAD**2-RC**2-DEL**2)
DISC=B**2-4.*A*C
IF (DISC.LT.O.) GO TO 50
RL=(-B+SQRT(DISC))/(2.*A)
RL1=(-B-SQRT(DISC))/(2.*A)
SL=S+SLOPE*(RL-R)
SL1=S+SLOPE*(RL1-R)
X=(RL-R)*SIN(THETA*RA)+(SL-S)*COS(THETA*RA)
X AND XI ARE USED TO CHECK IF THE INTERSECTION VALUES ARE UPWIND.
IF BOTH X AND XI ARE POSITIVE THEN THERE ARE TWO UPWIND
INTERSECTIONS. IF BOTH ARE NEGITIVE THEN THERE ARE NO UPWIND
00010
00020
00030
00040
00050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
OOltiO
00190
00200
00210
0022C
00230
00240
00250
002t>0
00270
00280
00290
00300
00310
00320
00330
00340
00350
003&0
00370
00380
00390
00400
00410
00420
00430
00440
00450
00460
00470
00480
00490
00500
00510
00520
00530
151
-------
c
c
c
40
50
60
C
CURLIN
INTERSECTIONS. IF EITHER ONE IS POSITIVE AND THE OTHER IS
NEGATIVE THAN THERE IS ONE UPWIND INTERSECTION.
X1=(RL1-R)*SIN(THETA*RA)+(SL1-S)*COS(THETA*RA)
NLOCI=1
IF (X.LT.O..AND.X1.LT.O.) GO TO 50
IF (X.LT.O.) GO TO 40
IF (Xl.LT.O.) NLOCI=0
GO TO 60
NLOCI=0
RL=RL1
SL=SL1
GO TO 60
NLOCI=-1
RETURN
END
00540
00550
00560
00570
00580
00590
00600
00610
00620
00630
00640
00650
00660
00670
00680
006*0
00700
152
-------
DIFANG
FUNCTION DIFANG (ANGA,ANGB) 00010
C DETERMINES THE DIFFERENCE BETWEEN TWO ANGLES. 00020
C THE RESULTING ANGLE IS BETWEEN 0 AND 360 DEGREES. 00030
DIFANG=ANGA-ANGB 00040
IF (DIFANG) 10,40,50 00050
10 IF (DIFANG+laO.) 20,30,40 00060
20 DIFANG=DIFANG+360. 00070
GO TO 10 OOObO
30 DIFANG=ldO. OOOyO
40 RETURN 00100
50 IF (DIFANG-130.) 40,30,60 00110
60 DIFANG=DIFANG-360. 00120
GO TO 50 00130
C 00140
END 00150
153
-------
ANGARC
FUNCTION AUGARC (DELM,DELN)
C ANGARC DETERMINES THE APPROPRIATE ARCTArJ (M/N) WITB
C THE RESULTING ANGLE BETWEEN 0 AND 360 DEGREES.
IF (DELN) 10,40,bO
10 IF (DELM) 20,30,20
20 ANGARC=57.2y57d*AIAN(DELM/DELN)+ldO.
RETURN
30 ANGARC=l80.
RETURN
40 IF (DELM) 50,60,70
50 ANGARC=27U.
RETURN
60 ANGARC=0.
RETURN
70 ANGARC=090.
RETURN
00 IF (DELM) 90,100,110
90 ANGARC=57.2957b*ATAN(DELM/DELN)+360.
RETURN
100 ANGARC=3bO.
RETURN
110 AHGARC=57.2957b*ATAN(OELM/DELN)
RETURN
C
END
00010
00020
00030
00040
0-0050
00060
00070
00080
00090
00100
00110
00120
00130
00140
00150
00160
00170
OOlbO
00190
00200
00210
00220
00230
00240
002t>0
154
* U.S. GOVERNMENT PRINTING OFFICE: 1978—740-261/335 Region No. 4
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/4-78-013
2.
3. RECIPIENT'S ACCESSION-NO.
4. TITLE AND SUBTITLE
USER'S GUIDE FOR PAL
A Gaussian-Plume Algorithm for Point, Area, and
Line Sources
5. REPORT DATE
February 1978
6. PERFORMING ORGANIZATION CODE
. AUTHOR(S)
William B. Petersen
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
10. PROGRAM ELEMENT NO.
1AA603 AB-25(FY-78)
11. CONTRACT/GRANT NO.
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory - RTP, NC
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, North Carolina 27711
13. TYPE OF REPORT AND PERIOD COVERED
In-house
14. SPONSORING AGENCY CODE
EPA/600/09
15. SUPPLEMENTARY NOTES
16. ABSTRACT
PAL is an acronym for this point, area, and line source algorithm. PAL is a
method of estimating short-term dispersion using Gaussian-plume steady-state
assumptions. The algorithm can be used for estimating concentrations of non-
reactive pollutants at 30 receptors for averaging times of from 1 to 24 hours,
and for a limited number of point, area, and line sources (30 of each type).
Calculations are performed for each hour. The hourly meteorological data
required are wind direction, wind speed, stability class, and mixing height.
Single values of each of these four parameters are assumed representative for
the area modeled.
This algorithm is not intended for application to entire urban areas but is
intended, rather, to assess the impact on air quality, on scales of tens to
hundreds of meters, of portions of urban areas such as shopping centers, lar
parking areas, and airports. Level terrain is assumed.
large
7.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS C. COS AT I Field/Group
* Air pollution
* Atmospheric models
* Algorithms
* Dispersions
13B
14A
12A
8. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (ThisReport)
UNCLASSIFIED
21. NO. OF PAGES
163
20. SECURITY CLASS (This page)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (9-73)
155
-------
03 ^
8-8
m
"O
o
o
-ti.
i
00
I
o
o o
CD 3. C O C
s*1 &s*
I ° 3 * I
rr. Fr\ •-v ^* CD
tf, O -. > <0
5^ g- CD (o
a: 5.1 5-
rr' m "^ •*••
Q^J O
o 5 5
>. C5 CD CD
^ a r»
^ I' §~T3
O ^ fb CU
"^ 5 g> Sj
>. o S- CD
8-a.g-o
.
ft •<
3
CD 3-
to S
f 5
CC3
CD
s-
O
S-
U)
-c
(D
O
-n
CD o
o c
o i
O
m
(A
J»
i?
31 O rn
If I
3 CD jo
O CD 2, ^
D' 2- -jj ^
^- — CD m
3 C? CD ,
oj w CD r?
i" CD -< f*
- ' fa o r~
O o =r -o
01
K>
(D
cx>
g O
3 CD
CD ^
— CD
° o
o ?
CD 5
3 CD
m
m
q
°
Z
>
CD
m
z
o
C
en
rn
Z
o o
I H
— ^
m
7 u)
H ni
> >
CO T) D
™ xO -«
01 r*v —i
u rn
H m
m (/)
n
H
O
Z
Q
m
Z
n
TJ
>
-------