Area Source Radiological Emission Analysis Code
(AREAC)
Technical Note
ORP-EAD-76-6
U. S. ENVIRONMENTAL PROTECTION AGENCY
Office of Radiation Programs
-------
Area Source Radiological Emission Analysis Code
(AREAC)
by
David Michlewicz
October 1976
U.S. ENVIRONMENTAL PROTECTION AGENCY
Office of Radiation Programs
Environmental Analysis Division
Protective Action Planning and Investigation Branch
Washington, B.C. 20460
-------
FOREWORD
The Office of Radiation Programs carries out a National Program
designed to evaluate the exposure of man to ionizing and nonionizing
radiation and to promote development of controls necessary to protect
the public health and safety and assure environmental quality.
In order to supplement our field measurement capabilities, we
rely on computer models to estimate, doses to man from releases of
radioactive materials. Existing air pathway models estimate doses
from point source releases of atmospheric radioactivity, i.e., where
the size of the emitting source is small and does not perturb the
resultant calculations. However, these models are not accurate at
estimating doses at close-in distances when the area of the source
is large as in the case of a uranium mill tailings pile, for example.
Therefore, this model has been designed to more accurately assess
close-in doses from large area sources.
Comments on this analysis would be appreciated. These should be
sent to the Director, Environmental Analysis Division, of the Office
of Radiation Programs.
Floyd L. Galpin
Director
Environmental Analysis Division
iii
-------
ABSTRACT
A computer code designed to calculate potential radiological
impact of atmospheric releases of radionuclides from area sources
is presented and discussed. The code is written in FORTRAN IV,
requires 48 K storage, and runs about 12 seconds on an IBM 370
system. The code can calculate radionuclide concentrations and
individual inhalation doses at up to six specific receptor loca-
tions and at up to 192 general locations around an area source.
Population doses can also be calculated.
The code accounts for area source shape, cloud diffusion,
ground and inversion-lid reflections, and radionuclide decay by
time of flight. It is dose model independent and requires a dose
conversion factor as part of input data to calculate doses propor-
tional to radionuclide concentrations.
The code is extensively annotated and simply written, hopefully
facilitating its use, and an effort was made to provide the user
with a high degree of flexibility in the utilization of this code.
v
-------
CONTENTS
Page
FOREWORD ill
ABSTRACT v
1. Introduction. . , 1
2. Mathematical Models 2
2.1 Diffusion of Airborne Emissions from Area Sources 2
2. 2 Diffusion Function , , . . . 4
2.3 AREAC Coordinate System 6
2.4 Diffusion Geometry , 7
2. 5 Receptor Location Geometry , , 7
2.6 Area Source Geometry, , 11
2.7 Diffusion at Large Distances from Area Sources 13
3. AREAC Input 15
4. AREAC Output 21
5. AREAC Structure 22
6. Sample Problem Description , 24
Appendix A - Point Source Approximation of an Area Source at
Large Receptor Distances 26
Appendix B - Glossary of Terms in AREAC , . . . 34
Appendix C - AREAC Sample Problem Input, Output, and Program
Listing , 37
References 55
FIGURES
Figure 1. AREAC Coordinate System 6
Figure 2. Geometry of AREAC - Plan View 8
Figure 3. Vertical Dispersion Coefficient as a Function of
Downwind Distance from the Source 10
Figure 4. Simplified AREAC Flow Chart 23
vii
-------
CONTENTS
Page
Figure 1A. Geometry of Long Term Sector Averaged Diffusion
from A Line Source Normal to Wind Direction 27
Figure 2A. Geometry of Diffusion from a Line Source
Tangential to Wind Direction. 30
TABLES
Table 1. Classification of Atmospheric Stability.. 9
Table 2. AREAC Input Card Sequence and Data Format 16
viii
-------
Area Source Radiological Emissions Analysis Code
(AREAC) ~ ~
1. Introduction
The Office of Radiation Programs in the Environmental Protection
Agency is developing a comprehensive dose computational system (CDCS)
for the analysis of emissions and effluents from the nuclear fuel
cycle. To assess the potential radiological effects of airborne
radionuclide releases from point sources, a computer code AIREM
has been developed by the Office of Radiation Programs in the past (JL) .
The code described in this program manual has been based to a large
extent on the AIREM code and represents an initial attempt at
developing a quantitative model for analyzing the potential radio-
logical impact of airborne constant, continuous releases of gaseous
radionuclides from area sources; principally inactive uranium tailings piles,
As presently written, the code consists of two main parts.
In the first part, a sector-averaged Gaussian diffusion equation
is utilized to calculate radionuclide dispersion coefficients
(X/Q's), radionuclide concentrations, and resultant inhalation
yearly doses at up to six specific receptor locations in the vicinity
of an area source. These, in practice, would correspond to the
maximally exposed individuals identified to be living near the pile.
In the second part, the same diffusion equation is employed to
calculate radionuclide dispersion coefficients, concentrations, and
yearly doses and population doses for 16 wind sectors and up to
12 downwind distances (192 sector-segments) around the area source.
-------
The manner in which the diffusion equation is utilized to account
for the distributed nature of an area source is described in
detail in the program manual.
The code was designed to be flexible by providing the user
with various options. The user can specify the shape of the area
source (rectangular or circular) and the degree of accuracy in
accounting for the distributed nature of the source. The code
provides an option for calculating inhalation doses only at
specific receptor locations or at locations in the population
wheel, or both. Since the code calculates radionuclide dispersion
coefficients and radionuclide concentrations at various desired
locations, by appropriate redefinitions of the dose conversion
factor, the user can adapt the code to calculate any desired quantity
which is proportional to radionuclide concentration. Thus, the code
can be adopted to calculate semi-infinite cloud gamma doses, and,
by setting the decay constant to zero, airborne concentrations of
non-radioactive gaseous emissions from area sources.
The code, utilized to its maximum capacity, requires 48 k bytes
of storage space, and runs about 12 seconds on an IBM 370. It is
written in FORTRAN IV.
2. Mathematical Models
2.1 Diffusion of Airborne Emissions from Area Sources
o
In general, the steady State concentration, X (Ci/M ), of
airborne radionuclides emitted from a point source can be specified
-------
at an arbitrary receptor location by the product of source strength,
Q (radionuclide rate of emission - Ci/sec), and a function X/Q,
which depends only on the mechanism by which the radionuclides
are transported from the point of emission to the receptor location.
Expressed mathematically:
X = Q • X/Q . (1)
Likewise, the concentration, dX, of airborne radionuclides,
which are emitted at the rate of q Ci/sec per unit of area from
a differential area element dA of a distributed source is:
dX = X/Q qdA (2)
where both X/Q and q are functions of dA for a fixed receptor location.
The concentration of airborne radionuclides emitted from the
entire area source may be obtained by integrating equation 2 over
the surface of the source. Performing the integration:
X = fx/Q • qdA
A
or, if q is independent of position,
X = qfX/Q dA (3)
A
Solution of equation 3 requires a two or, if the height of
the area source is not uniform, three dimensional integration of
a complicated diffusion function. However, an approximate solution
may be obtained by approximating the integral with a summation of
the integrand over the source area. Thus, if the source area, A,
-------
is divided into n small area elements, then
X - q E X/Q AA
n k k
where:
AA, is the area of the k area element and X/Q^ is the
value of the diffusion function X/Q at source position AA, .
If all the area elements are of the same size, AA, they can
be taken outside the summation sign, and the radionuclide concen-
tration at the fixed receptor location is:
X - qAA EX/Q
k k
Furthermore, if Q is the rate of radionuclide emission from the
entire source, it follows that
X = £ £ X/Q (4)
n n k
If one defines the function X/Q = _! E X/Q, to be an area source
A n n *
diffusion function, then
X * Q • X/QA (5)
and the similarity of equation 5 to equation 1 becomes apparent.
The definition of X/Q. is used in AREAC to calculate area
A
source dispersion coefficients (X/Q.'s) and corresponding
J\
radionuclide concentrations and doses at receptor locations
specified by the code user.
2.2 Diffusion Function
The point source diffusion function used in AREAC to
-------
calculate the area source diffusion function is a standard sector-
averaged Gaussian diffusion equation (_2,3) modified to include
radionuclide decay by time of flight.
-------
2.3 AREAC Coordinate System
The basic coordinate system used in AREAC is a two-dimensional
rectangular coordinate system in which the positive x-axis points
in the eastern direction, the positive y-axis points in the
northern direction, and the origin is in the center of the area
source. In instances where the polar coordinate system is used,
its coordinates are defined in terms of the rectangular coordinate
system as shown in the following figure:
N
(x,y)
(r,6)
0
Figure. 1. AREAC Coordinate System
-------
2.4 Diffusion Geometry
The geometrical scheme of AREAC is similar to a cartwheel.
It is illustrated in figure 2. Winds blow the emissions from
each source area element down pie-shaped sectors centered on
the area element. The concentration of radionuclides in
the sectors is assumed to have a Gaussian distribution in the
vertical direction, centered at the effective release height,
and to be uniform in the horizontal direction across each
sector. As illustrated in figure 2, there are 16 sectors
corresponding to 16 compass directions (N, NNE, NE...NNW) toward
which the winds blow.
Atmospheric stability, which determines the standard deviation
(a ) of the vertical distribution of radionuclide concentration,
Z
is characterized by six stability classes which' are described in
table 1.
The concentration standard deviation in the vertical direction
increases monotonically with distance from the source area element.
according to figure 3 (3^), until it is set equal to 0.8 of the
mixing layer height (3). This value determines the maximum
extent of cloud diffusion in the vertical direction. Values of
mixing layer heights may be obtained from references (1) and (4_).
2.5 Receptor Location Geometry
AREAC consists of two main parts. In the first part, area
source dispersion coefficients (X/Q 's), radionuclide concentrations,
-------
N
NNE
NE
UP TO 12
ANNULAR RINGS
STABILITY
CLASSES
PER SECTOR
SEGMENT
ENE
Figure 2. Geometry of AREAC - plan view
-------
Table 1. Classification of Atmospheric Stability
Stability
classification
Extremely unstable
Moderately unstable
Slightly unstable
Neutral
Slightly stable
Moderately stable
o
Pasquill 09 Temperature change
categories (degrees) with height
(°C/100 m)
A 25.0 <-1.9
B 20.0 -1.9 to -1.7
C 15.0 -1.7 to -1.5
D 10.0 -1.5 to -0.5
E 5.0 -0.5 to 1.5
F 2.5 1.5 to >4.0
Standard deviation of horizontal wind direction fluctuation over
a period of 15 minutes to 1 hour. The values shown are averages for
each stability classification.
Reference: AEC Safety Guide 23.
-------
1,000
100
I-J
b
DISTANCE DOWNWIND, kr
Figure 3. Vertical Dispersion Coefficient as a Function of
Downwind Distance from the Source (3)
10
-------
and yearly individual inhalation doses are calculated at up to six
(6) receptor positions specified by the code user. These locations
have to be specified in polar coordinates defined in section 2.3.
In the second part, area source dispersion coefficients,
radionuclide concentrations, individual and, if desired, population
doses are calculated at up to 192 locations in a population wheel
around the source. The population wheel, which is illustrated in
figure 2, has its center in the center of the area source and
is divided into sixteen (16) 22.5 degree sectors and up to twelve
(12) annular rings. The sixteen sectors and up to twelve annular
rings form a circular grid of up to 192 sector segments. Individual
doses within this grid are calculated at the centers of the sector
segments. The contribution to the total population dose from each
sector segment is the product of the dose and population within the
sector segment.
2.6 Area Source Geometry
The code user has the option of approximating the shape of
the area source by either a rectangle or a circle. In either case,
as discussed in section 2.1, the source area is divided into a
number of equal area elements to account for the distributed nature
of the source.
A. Rectangular Area Source
If the shape of the source is approximated by a rectangle,
11
-------
it is assumed to have two sides' of length XL parallel to the
x-axis and two sides of length YL parallel to the y-axis (see
section 2.3). It is divided into a grid of NX by NY area elements,
where NX is the number of divisions along the x-axis and NY is the
number of divisions along the y-axis, and radionuclides are assumed
to be emitted from the centers of the area elements.
NX and NY can range from 1 to 10. If both NX and NY are
equal to 1, the radionuclides are assumed to be emitted from the
center of the area source. This is equivalent to approximating an
area source by a point source.
By choosing appropriate values of XL and YL, distributed
sources ranging in shape from squares to finite line sources can
be approximated.
B. Circular Area Source
If the shape of the source is approximated by a circle of
radius RL, then NX is the number of sectors and NY is the number of
annular rings into which it is subdivided, and radionuclides are
assumed to be emitted from the centers of the NX • NY sector segments
which are formed. The division of the circle into sectors is,
except for a minor difference, similar to the division of the
population wheel.
In order to assure that the areas of all the sector segments
are of the same size, the following procedure was chosen for
dividing the circle into the NY annular rings.
12
-------
If AA is the area of a sector segment, then NXAA is the
area of an annular ring, and since AA has to be constant, the
areas of all annular rings must be the same.
If RL is the radius of the circular area source, then the
area of the source
2
A = TT RL
and, since there are NY annular rings of equal area, the area of
a ring
irRL2
R NY
Thus, if R_ is set equal to zero and R. is the radial
distance to the outer boundary of the i annular ring, then
p2 D2 TTRL2
TrR. = irR. . H --
1 i-1 NY
or, dropping the constant TT ,
. =./R2 , + RL /NY
i V i-l
This procedure for determining annular boundaries will assure
that all area elements in the circular source grid are of the same
size.
2.7 Diffusion at Large Distances from Area Sources
In a problem in which the area source is divided into 10 by
10 area elements and dose calculations are performed at 6 specific
and 192 (16 sectors x 12 annular rings) general receptor locations,
there are 10 x 10 x (192 + 6) = 19,800 source to receptor distances
for which the point source diffusion equation has to be solved.
13
-------
Since at large receptor distances an area source approaches
a point source, this amount of calculation is most often
unnecessary, and to save computation time (and money), AREAC
provides the code user with an option for treating the area
source as a point source at large receptor distances.
The manner in which this is accomplished is outlined in
Appendix A.
14
-------
3. AREAC Input
Up to 80 lines of input data are required to run AREAC.
The input card sequence and data format are summarized in
table 2 and illustrated in the sample program in Appendix C.
First Card
This is a title card which identifies the input problem.
No calculations are made with data on this card.
Second Card
KGEOM is an area source geometry number. If KGEOM = 0, the
shape of the source is approximated by a rectangle. If KGEOM = 1,
the shape of the source is approximated by a circle.
Third Card
If the source has the shape of a rectangle (KGEOM = 0), the
card should contain the X-length, Y-length, and height of the source
(meters). If the source has the shape of a circle (KGEOM = 1), the
card should contain the radius and height of the source (meters).
Fourth Card
Source strength (Ci/s), radionuclide decay constant (s ), and
3
dose conversion factor (mrem/yr per Ci/m ) are required on this
card.
Fifth Card
Area Source Division Numbers
If source has the shape of a rectangle, NX and NY are,
respectively, the number of divisions in the X and Y directions
into which the source is subdivided.
-------
Table 2. AREAC Input Card Sequence and Data Format
Card Sequence
I card
1 card
1 card
1 card
1 card
1 card
16 cards
16 cards
1 card
1 card
Columns
1-20
21-25
27-30
31-35
37-40
41-45
1
1-10
11-20
21-30
or
1-10
11-20
1-10
11-20
21-30
1- 5
6-10
1
1-60
1-60
1-10
1
Description3
Source name
Number of months of data
Beginning month
Beginning year
Ending month
Ending year
KGEOMa (0 or 1)
X-length of area source (m)J
Y-length of area source (m)JKGEOM=0
Height of source (m) \
Radius of area source (m)| „„_„,, ,
T, . , . r / \ f iMjtUrl— _L
Height of source (m) I
Source strength (Ci/s)
Radionuclide decay constant (s )
Dose conversion factor (mrem/yr
per Ci/m3)
NXa (10 max.)
NYa (10 max.)
KDIFa
Wind frequencies by stability class
(%) , 10 columns each
Wind speeds by stability class
(m/s) , 10 columns each
SIGMAX3 (m)
KPROBa
Format
5A4
F5.0
A4
15
A4
15
11
F10.0
F10.0
F10.0
F10.0
F10.0
E10.3
E10.3
E10.3
15
15
11
6F10.0
6F10.0
F10.0
11
See glossary in Appendix B.
16
-------
Table 2. AREAC Input Card Sequence and Data Format (continued)
Card Sequence
Columns
Description5
Format
1 card
1 card
1 card
1 card
1-60
1-60
1- 2
MRAD cards
2 MRAD cards
1-30
1-80
Number of specific receptors H
(IND) (6 max.)
Receptor radial positions (m), 6F10.0
10 columns each
Receptor angular positions (degrees) 6F10.0
10 columns each
Number of annuli in poulation wheel 12
(MRAD) (12 max.)
Midpoint, lower and upper radii (m), 3F10.0
10 columns each
Annular ring population, 8 sectors 8F10.0
to a card, 10 columns each
17
-------
If source has the shape of a circle, NX and NY are, respectively,
the number of annular rings and sectors into which the source is
subdivided.
Sixth Card
KDIF is an area source diffusion number which provides the user
with an option for treating the area source as a point source at
large source to receptor distances. If KDIF = 1, the area source
is treated as a distributed source at all receptor locations.
If KDIF ^ 1, the area source is treated as a point source at
large receptor distances.
Wind Rose Data
Wind frequencies (%) and wind speeds (m/s) are entered in a
clockwise direction beginning with the north sector. Sixteen
frequency and sixteen wind speed cards are required (1 card per
sector). Each card contains six entries which correspond to
stability classes A through F. Wind direction is the direction
toward which the wind blows.
SIGMAX Card
SIGMAX sets the maximum value of the standard deviation in
the vertical direction. It should be equal to 0.8 times the
height of the inversion layer (1,3).
Receptor Set Identification Card
KPROB is a number which identifies the set of receptors for
which dose calculations will be performed. If KPROB = 1, only
18
-------
dose calculations for specific receptors will be performed. If
KPROB = 2, only-dose calculations for the population wheel will be
performed. If KPROB is any other number, calculations for both
specific receptors and the population wheel will be performed.
Specific Receptor Number Card
If KPROB = 2, this card should be omitted.
This card should contain the number of specific receptors (IND)
for which dose calculations are to be performed (1 < IND £ 6).
Specific Receptor Radial Positions
If KPROB = 2, this card should be omitted.
This card should contain the distances in meters of receptors
1 through IND as measured from the center of the area source.
Specific Receptor Angular Positions
If KPROB = 2, this card should be omitted.
This card should contain the angular positions in degrees of
receptors 1 through IND as measured in the clockwise direction from
the y-axis (see section 2.3). If KPROB = 1, this is the last data
card.
Population Wheel Annuli Card
This card should contain the number of annuli (MRAD) in the
population wheel (1 £ MRAD < 12). The number of sectors is always 16
Annular Boundary Cards
These MRAD cards should contain the radial distances from the
center of -the area source to the midpoints, inner, and outer
19
-------
boundaries of the MRAD population wheel annuli. The minimum midpoint
radius which can be used is 100 m + DTEST/2 .
Population Data Cards
Population data are entered clockwise starting in the
north sector and within inner annulus. Two MRAD cards are
required. If KPROB ^ 1, this is the last set of data cards.
See glossary in Appendix B for definition of term.
20
-------
4. AREAC Output
The output of AREAC consists of four main parts: source input
data, meteorological input data, specific receptor input data and
calculational results, and population wheel input data and calcula-
tional results.
If no calculations are desired for either a set of specific
receptors or a population wheel, then that part is omitted. Error
messages are printed out, and program is terminated if the source
area shape is not specified or the number of specific receptors is
outside of the allowed range. If a specific receptor is less than
100 m from the area source, calculations for that receptor are
omitted. However, if this occurs at a location in the population
wheel, the program is terminated. The sample problem in Appendix C
illustrates the output sequence and format.
21
-------
5. AREAC Structure
AREAC consists of a main program, subroutines DIFUSN, GEOMO
and GEOM1, and function SIGZ. Figure 4 presents a simplified
program flow chart.
All input and output operations and many calculations are
performed in the main program. If the source has the shape of a
rectangle, subroutine GEOMO is called by the main program to calculate
the centers of area elements into which the area is subdivided. If
the source has the shape of a circle, this function is performed by
subroutine GEOM1. Subroutine DIFUSN is called by the main program
to calculate an area source X/Q at a receptor location specified
by the main program. DIFUSN also informs the main program when a
receptor is too close to the area source. Function SIGZ (_!) is
called by DIFUSN to provide the concentration standard deviation
in the vertical direction (a ) for a receptor location and meteoro-
z
logical stability class specified by DIFUSN. Communication between
the main program and subroutines is accomplished through the use
of a blank COMMON.
22
-------
Read/Write
'Source Input Data,
Determine Coordinates of Centers
of Source Area Elements
Subroutine
GEOMO
GEOM1
Read/Write
'Meteorological Input Data
True
/ Read/Write
Specific Receptor Input Data
Calculate
X/Q's, Concentrations, Doses
Write
'x/Q's, Concentrations, Doses
True
False
Read/Write
'Population Wheel Input Data
Calculate
X/Q's, Doses, Population Doses
Write
'X/Q's, Doses, Population Doses/
END
Figure 4. Simplified AREAC Flow Chart
23
-------
6. Sample Problem Description
AREAC was employed in a sample problem to assess the potential
radiological impact of airborne emissions of radon-222 from a
radium bearing phosphate gypsum pile. The shape of the pile
was approximated by a circular disc with a radius of 590 m and an
average height of 29 m. The pile area was divided into a grid of
10 by 10 area elements and, at large receptor to source distances,
the pile was treated as a point source. Based on measurements
of radon exhalation rate from the pile, the pile emission rate
of radon-222 was taken to be approximately 4.3 x 10 Ci/s. A
12
concentration to dose conversion factor of 4 x 10 mrem/yr per
3 3
Ci/m (4 mrem/yr per pCi/m ) to the bronchial epithelium region
of the lung was assumed (5).
Meteorological data for the pile area, covering a period of
60 months from January 1969 to December 1973, was obtained from the
Environmental Data Service, National Oceanic and Atmospheric
Administration, United States Department of Commerce, and a
SIGMAX value of 1000 meters was assumed. The code was used to
calculate yearly doses at six specific locations west of the
pile and to calculate yearly doses, population doses, and the total
population dose in a population wheel comprised of 16 sectors and
12 annuli extending to a distance of 80 kilometers from the pile
The population distribution was 'Obtained from the Census Bureau's
Master Enumeration District List with Coordinates as edited by the
24
-------
Office of Telecommunications of the U.S. Department of Commerce,
The input data for this problem is shown in Appendix C together
with the output and the program listing.
25
-------
APPENDIX A
Point Source Approximation of an Area Source
at Large Receptor Distances
This appendix presents the derivation of a criterion for
treating an area source as a point source at large sources to
receptor distances. The criterion is derived by showing that
in calculating the yearly average radionuclide concentrations
at receptor positions beyond a certain distance from an area
source, both the crosswind and the downwind extent of the source
can be neglected.
Consider a ground level line source of airborne emissions
of length L, which is bisected by line AB which is perpendicular
to the source (figure 1A). Let point A be the point of intersection,
and let line AB bisect a wind rose sector whose apex is located at
point A.
By method analogous to that which is discussed in section 2.1,
the yearly average radionuclide concentration at a point located on
line AB may be calculated, approximately, by dividing the line
source into a number of point sources and summing the contribution
of each source to the radionuclide concentration at the receptor
location.
If the line source is assumed to be a point source located at
the center of the line source (point A), then the radionuclide
concentration at a receptor location on line AB is obtained by
26
-------
FIGURE 1 A. GEOMETRY OF LONG TERM SECTOR AVERAGED DIFFUSION FROM A LINE SOURCE
NORMAL TO WIND DIRECTION
27
-------
considering the contribution of only point A.
In both cases, the procedure for calculating the dispersion
coefficients from an individual point is identical. The point to
receptor distance and wind rose sector orientation are determined,
and the sector averaged Gaussian diffusion equation is solved
for the appropriate distance and joint wind-speed and direction-
frequency distribution.
If the line source is assumed to be a point source, a receptor
will always be located in the same wind rose sector. If, for
example, line AB points north, then a receptor will always be in
the northern sector regardless of its distance from the source.
For a point located on a line source, however, the sector
orientation of a receptor on line AB will depend on its distance
from the source. For example, a receptor located immediately to
the north of point A will be in the eastern sector, as viewed from
point C, and will approach the northern sector with increasing source to
receptor distance. At the distance where the receptor is in the
same sector with respect to points A and C, it is in the same sector
for all points between A and C and, by symmetry, all points on the
line source.
If, in figure 1A angle FCE is the angular width of a wind
rose sector centered on line CG, then point E will be in the same
sector with respect to all points on the line source. If 9 is the
width of a sector, then according to elementary geometry, angle CEA
28
-------
is equal to 9/2, and the length of line segment AE is:
— _ L
2 tan (8/2) .
The distance from a point on the line source to point E
ranges from the length of line segment AE to the length of
line segment CE. If R is the ratio of CE~ to AEf, then
R =
cos(0/2) .
In a wind rose consisting of 16 sectors, 6 = 22.5 degrees and
AE = 2.51 L
and R = 1.02.
Therefore, since beyond point E the joint wind speed and
direction-frequency distribution is identical for all points
on the line source, and the source to receptor distance is
constant to within 2 percent, distance AE may be taken as
the distance beyond which the line source can be treated as a
point source.
Since a distributed area source can be thought of as an
infinite number of line sources, the above discussion indicates
that if L is the width of the source in the crosswind direction,
the crosswind extent of the source can be neglected beyond a
distance of 2.51 L from the source. This will be true of receptor
Indeed, the joint wind speed and direction-frequency distribu-
tion becomes identical at point H, but at that distance, the source
to receptor distance can no longer be assumed to be constant.
29
-------
locations anywhere within a sector since the joint wind speed
and direction-frequency distribution is constant within a sector.
It can also be shown that, at the distance of 2.51L, the extent of
the area source in the prevailing wind direction can also be neglected.
Consider a ground level line source of length L extending
from -L to +L along an axis pointing in the prevailing wind direction.
1 2
The origin of the axis coincides with the center of the line source,
and a receptor is located at a distance D from the center of the
source. This situation is illustrated in figure 2A.
Figure 2A. Geometry of Diffusion from a
Line Source Tangential to Wind Direction
wind
-L
2
0 +1
e
I
>
direction
)
The concentration of radionuclides at point D can be obtained
by integrating the point source difussion function along the line
source from D -L^ to D +L. It should also be possible to obtain the
2 2
concentration at point D by assuming that all radionuclides are
emitted from some point located between -L_ and +L and at a distance
2 2" I
X from point D.
2
Since the point source diffusion function is continuous on
F-L, +_LJ , then, according to the Mean Value Theorem for integrals,
such a point must exist (6).
30
-------
If f(x) is the diffusion function which gives the radionuclide concen-
tration at a distance x from the point of emission, and q is the radio-
nuclide release rate per unit of source length, then xo must satisfy
the following equation:
Lqf(x0) =\ q f(x)dx . (2)
Assuming a uniform release rate one can see that:
(3)
or that x is that distance from the receptor at which f(x) is equal to
its average value over the source length.
For low source height, the long term sector averaged diffusion
function is proportional to (xaz) , where az is the concentration
standard deviation in the crosswind direction. If az is considered to
be proportional to x°, then f(x) is proportional to x , where c is a
constant and r = 1+c. Using figure 3-3, reference 3, and a source to
receptor distance of 1 km, the values of c, the scope of the curve,
and r, as a function of stability class, are approximately:
31
-------
Stability Class
A 2.1 3.1
B 1.1 2.1
C .92 1.92
D .67 1.67
E .61 1.61
F .57 1.57
When the x r expression for f(x) is substituted into equation 3,
then, after carrying out the integration and simplifying the result,
one obtains:
r in1-"
£ L_±
L
2Dl - 1 - 2D< . (4)
1-r
At the source to receptor distance at which the crosswind extent of the
area source can be neglected L_ = Tan (11.25°), and substituting in the
2D
appropriate values for r one finds that the value of — ranges from
0.97 for stability class A to 0.98 for stability class F. At larger
source to receptor distances the value of L approaches 0 and conse-
Xo 2D 3
quently, the value of -— quickly approaches unity.
Since D is the distance from the receptor to the center of the line
source, it has thus been shown that a distance approximately equal to
2.51 times the greatest dimension of an area source, the extent or the
source in both the crosswind and downwind directions can be neglected,
and all radionuclides can be assumed to be emitted from the center of
the area source.
o
This can easily be shown by expanding terms on the right side
of equation 4 in Taylor series (6_) .
32
-------
In AREAC, the results of the previous discussion have been
utilized to provide the code user with an option for treating
an area source as a point source at large receptor distances.
This calculational "shortcut" is based on a modified form of
equation 1 in which L was chosen to be the largest projection of
the area source in the crosswind direction. If the source has
the shape of a rectangle, then L is its diagonal, and if the
source is a circle, then L is its diameter. Also, since the
smallest distance between the edge of an area source and a
receptor located at a distance D from the center of the source
can be equal to D-L^, the distance from the center of an area
"2
source at which the point source assumption can be made has been
conservatively extended by I,. Thus, the area source is treated
2
as a point source at receptor distance x at which
x > L \l + 1 I
2 tan(11.25°) I.
33
-------
APPENDIX B
Glossary of Terms in AREAC
TERM
APOS(I)
CHIOQ
CHIQ(I,J)
CONG
DCF
DEC
DIR(I)
DOS(I,J)
DOSE
DTEST
FACIL
FREQ(I,J)
IND
DESCRIPTION
Angular position of I specific receptor (degrees)
Ratio of concentration to release rate at a
specific receptor location
Ratio of concentration to release rate at the
center of population wheel sector segement (I,J)
3
Radionuclide concentration in air (Ci/m )
3
Dose conversion factor (mrem/yr per Ci/m )
Radionuclide decay constant (s )
Compass direction of the I sector (N, NNE...NNW)
Dose at center of population wheel sector
segment (I,J) (rem/yr)
Dose at a specific receptor location
(mrem/yr)
The longest possible projection of the area source
in a crosswind direction (m). If source has the
shape of a circle, DTEST is its diameter. If
source has the shape of a rectangle, DTEST is its
diagonal. It is used in approximating the area
source by a point source at large receptor to
source distances (see Appendix A).
Area source description (name, etc.)
Frequency that the wind blows toward sector I
in stability class J (%)
Number of specific receptors
34
-------
KDIF
KGEOM
KK
KPROB
MONTHS
MONTH1
MONTH2
MRAD
MSEC
NX
NY
NYR1
NYR2
POP(I,J)
POPDOS(I,J)
Area source diffusion number. If KDIF = 1,
the area source is treated as a distributed
source at all receptor locations. If KDIF / 1,
the area source is treated as a point source
at large source to receptor distances.
Area source geometry number (0 or 1). If
KGEOM = 0 or 1, the shape of the source is
assumed to be a rectangle or a circle respectively.
Number by which the main program is informed
whether a receptor is located too close to the
area source
Receptor set identification number. If KPROB = 1,
only calculations for specific receptors are
performed. If KPROB = 2, only locations in the
population wheel are considered. If KPROB ^ 1 and
KPROB ^ 2, calculations for both specific receptors
and the population wheel are performed.
Number of months represented by input data
Beginning month represented by input data
Ending month represented by input data
Number of annuli in the population wheel (12 max.)
Number of sectors in the population wheel (16 always)
Area source division numbers. If area source
has the shape a rectangle, NX and NY are,
respectively, the number of divisions in the X
and Y directions (see section 2.6) into which
the source is subdivided.
If area source has the shape of a circle, NX and
NY are, respectively, the number of annular rings
and sectors into which the source is subdivided.
Beginning year represented by input data
Ending year represented by* input data
Population in population wheel sector segment (I,J)
Population dose received by population
in sector segment (I,J) (1000 person-rem/yr)
35
-------
RINN(I)
RL
RMID(I)
ROUT(I)
RPOS(I)
SIGMAX
SOURCE
SPEEB(1,J)
SX(I,J)
SY(I,J)
SZ
WSA
XL
XPOS(I)
YL
YPOS(I)
ZH
Radial distance to inner boundary of I annular
ring from center of area source (m)
Radius of circular area source (m)
Radial distance to middle of I annular ring
from center of area source (m)
Radial distance to outer boundary of I annular
ring from center of area source
Radial position of I specific receptor (m)
Maximum value of SZ
Area source radionuclide emission rate (Ci/s)
Average wind speed in sector I for stability
class J (m/s)
X-coordinate of center of source area element
Y-coordinate of center of source area element
(I,J) (m)
Standard deviation of concentration in the
vertical direction, a
Angular width of a sector (22.5 degrees)
Width of a rectangular area source in the
X-direction (m^
X-coordinate of I ; specific receptor (m)
Width of a rectangular area source in the
Y-direction (m)
Y-coordinate of the I specific receptor (m)
Height (average) of the area source (m)
36
-------
APPENDIX C
AREAC sample problem input, output, and program listing
37
-------
AREAC
AREA SOURCE RADIOLOGICAL EMISSIONS ANALYSIS CODE
ENVIRONMENTAL PROTECTION AGENCY
OFFICE OF RADIATION PROGRAMS
ENVIRONMENTAL ANALYSIS DIVISION
401 M STREET, S.W.
WASHINGTON, D.C. 20460
SOURCE INPUT DATA
FACILITY, NO. MONTHS OF DATA, PERIOD
C SAMPLE PILE 60. JAN 1969 DEC 1973
SOURCE AREA HAS THE SHAPE OF A CIRCLE
PILE RADIUS = 590. M. , HEIGHT = 29. M.
SOURCE STRENGTH,CI/SEC DECAY CONSTANT,I/SEC DOSE CONVERSION FACTOR,MREM/YR PER CI/CU.M.
4.280E-06 2.100E-06 4.000E+12
THE SOURCE AREA HAS BEEN SUBDIVIDED INTO 10 ANNULAR RINGS AND 10 SECTORS
AT LARGE DISTANCES THE SOURCE AREA WILL BE TREATED AS A POINT SOURCE
-------
METEOROLOGICAL INPUT DATA
WliJD
DIK
N
NNt
Nt
ENE
E
ESE
SE
SSE
S
SSw
sw
wsrf
w
WNM
NW
NNvJ
WINO
Dirt
N
NNt
INlE
ENc
E
ESE
iE
SSE
S
ss*
sw
MS*
w
WNX
NM
NNw
FREQUENCY
A
0.02
0.01
0.01
u.03
0.02
0.0
0.0
U.02
0.02
0.02
0.02
0.01
0.02
0.0
0.01
0.0
SPEEDb IN
A
1.06
1.69
1.69
1.62
1.59
0.0
0.0
1.59
1.59
1.59
1.59
1.69
1.59
0.0
1.69
0.0
IN PERCENT
B
0.3T
0.35
0.40
0.58
1.02
0.19
0.16
0.17
0.31
0.31
0.28
0.29
0.46
0.42
0.38
0.40
METERS PER
B
2.18
2.56
2.86
3.08
2.88
2.16
1.89
2.14
2.21
2.31
2.26
2.12
2.15
1.80
2.19
2.25
BY STABILITY CLASS
C
0.85
0.86
0.71
1.12
2.61
0.52
0.41
0.41
0.51
0.50
0.73
1.11
1.95
1.02
0.81
0.98
SECOND BY
C
3.78
3.81
3.75
3.58
4.45
3.60
4.40
3.64
3.43
3.60
3.58
3.55
3.46
3.70
3.19
3.64
D
2.50
1.43
0.93
1.18
3.62
1.94
2.04
1 .86
1.98
1.69
2.64
3.43
=,.68
2.91
2.38
1.95
STABILITY
D
4.18
4.44
3.82
4.39
4.70
4.75
5.14
5.09
4.38
4.26
4.38
4.17
4.18
4.21
4.27
4.21
FOR EACH SECTOR
E
0.80
0.35
0.47
0.36
0.82
0.77
1.13
0.86
1.18
0.92
1.38
2.36
3.25
1.69
1.05
0.85
CLASS FOR EACH
E
2.98
2.92
2.88
3.25
3.25
3.17
3.29
3.25
3.32
3.29
3.26
3.14
3.09
3.06
3.20
3.00
F
0.74
0.34
0.44
0.42
0.66
0.72
1.29
1.01
1.45
1.04
1.95
3.04
4.67
2.13
1.23
1.03
SECTOR
F
1.30
1.24
1.48
1.28
1.35
1.45
1.49
1.37
1.44
1.39
1.44
1.41
1.42
1.26
1.21
1.35
SIuMAA= 10UO. METERS
-------
FORIRAN IV G LEVEL 21
MAIN
DATE = 76012
14/57/19
PAGE 0001
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
COMMON NX,NY,SIGMAX,SX(10»10),SY<10»10),FREQ<16,6),SPEED(16,6) ,ZH»
1 PI.DEC,DTEST,KDIF
DIMENSION RPOS(6),APOS(6>,XPOb<6>,YPOS(6>,DIR(16).R1NN(12>,ROUT(12
1).005(12.16).POPDOS(12,16).CH10(12.16).POP(12,16),RMID(12),FACIL<5
1)
REAL MONTHS
Pl=3.141593
DATA DIH/'N '.'NN.E '.'NE <»iENE •»'£ '.'ESE '»'SE '»«SSE ','
IS ','SSrt '»'SW ','WSW ','W «,'WNW ','NW '»«NNW •/
WRITE(6,10>
10 FORMAT(1H1,T32»'AREAC'//
1 1X,T12.'AREA SOURCE RADIOLOGICAL EMISSIONS ANALYSIS CODE'/
2 lH*,T12,'_',Tl7,'_',Ta4,'_',T37,'_',T47,'_',T56.'_'//
3 IX,T20.'ENVIRONMENTAL PROTECTION AGENCY1/
4 IX,T20,'OFFICE OF RADIATION PROGRAMS'/
5 IX,T20,-ENVIRONMENTAL ANALYSIS DIVISION'/
6 1X,T20,'401 M STREET,S.W.'/
7 IX,T20.'WASHINGTON O.C. 20460'//)
READ/WRITE SOURCE DATA
WRITE16.20)
20 FORMAT(1HO,'SOURCE INPUT DATA'/)
WHITE(6,30)
30 FOKMATUHO,'FACILITY, NO.MONTHS OF DATA, PERIOD')
RE AD (5, 40) FACIL. MONTHS, MONTH1 ,NYR1 .MONTH2.NYR2
40 FOf
-------
FORTRAN IV G LEVEL 21
MAIN
DATE = 76069
16/53/04
PAGE 0002
-p-
Ul
C READ/WRITE SOURCE STRENGTH.RADIONUCLIDE DECAY
C CONSTANT(DEC)'AND DOSE CONVERSION FACTOR(DCF)
0035 1010 READ(5,130)SOURCE«OrC,DCF
0036 130 FORMATOE10.3)
0037 _ WRITE<6«140>
0038 140 FORMAT)'OSOURCE STRENGTH,CI/SEC'.T30t•DECAY CONSTANT,1/SEC',T60,'D
10SE CONVERSION ?ACTOR.MREM/YR PER CI/CU.M.')
0039 WRITE(6t150>SOURCE.DEC,DCF
0040 150 FORMAT) T5 , 1PE 1 0 . 3 , T35, 1PE.1 0 .3, T65, 1PE1 0 . 3/>
C READ/WRITE NX AND NY.
C IF THE PILE IS A RECTANGLE.NX IS THE NUMBER OF DIVISIONS
C ALONG THE X AXIS.ANf) NY IS THE NUMBER OF DIVISIONS ALONG
C THE Y AXIS INTO WHICH THE PILE HAS BEEN SUBDIVIDED.
C THE POSITIVE X-AXIS IS ORIENTED TOWARD THE EAST
C THE POSITIVE Y-AXIS IS ORIENTED TOWARD THE NORTH
C THE ORIGIN OF THE X-Y COORDINATE SYSTEM IS IN THE
C CENTER OF THE AREA ^OURCE
c IF THE PILE is A CIRCLE.NX is THE NUMBER OF ANNULAR RINGS,
C AND NY IS THE NUMBER OF SECTORS INTO WHICH THE PILE HAS
C BEEN SUBDIVIDED.
0041 _ REAO(5,160)NX.NY
0042 160 FORMATI2I5)
0043 ... IF (KGEOM.EQ.O) WRITF(6,I70) NX,NY
0044 170 FORMAT(1HO,'THE SOURCE AREA HAS BEEN SUBDIVIDED INTO ',13,' BY',13
i«' RECTANGULAR AREA ELEMENTS'/)
0045 IF(KGEOM.EO.l) WftlTF(6,175) NX,NY
0046 175..FORMATUHO, 'THE SOURCE AREA HAS BEEN SUBDIVIDED INTO ',I3,« ANNULA
1R RINGS AND',13,' SECTORS'/)
C CALCULATE THE X-Y COORDINATES OF SOURCE AREA ELEMENTS
0047 IF(KGEOM.EQ.O) CALL GEOMO(XL.YL)
0048 IF(KGEOM.EO.1) CALL GEOMl(RL)
c READ KDIF.IF KDIF is NOT EQUAL TO I.THE SOURCE AREA WILL BE
C TREATED AS A POINT SOURCE AT LARGE DISTANCES FROM THE PILE.
C IF KDIF IS EQUAL TO It THE SOURCE AREA WILL BE TREATED AS A
C DISTRIBUTED SOURCE EVEN AT LA^QE DISTANCES FROM THE PILE.
0049 ' READ(5,"l80) KDIF
0050 180 FORMAT(Il)
0051 IF(KOIF.EQ.1) GO TO 1020
0052 . WRITE(6,190>
0053 190 FORMAT(1HO,'AT LARGE DISTANCES THE SOURCE AREA WILL BE TREATED AS
1A POINT SOURCE')
0054 f>Q TO 1030
0055 1020 WRITE(6,200)
0056 200 FORMAT11HO,'SOURCE AREA WILL BE TREATED AS A DISTRIBUTED SOURCE IN
1 ALL CALCULATIONS')
C READ/WRITE METEOROLOGICAL DATA
C THERE ARE ALWAYS 16 COMPASS DIRECTIONS
c nNo 6 METEOROLOOICAI STABILITY CLASSES
00-t-7 1030 WOITE(iS,?.10)
00^3 210 FORMAT(]H1,'METEOROLOGICAL INPUT.DATA')
C READ/WRITE WIND FREQUENCIES BY 'STABILITY CLASS
0059 00 10401=1.16
0060 READ(5,220)(FREQ(i,J)»J=l,6)
0061 220 FORMAT(6F10.0)
0062 1040 CONTINUE
0063 WRITE16.230)
0064 230^FORMAT(1HO,'WIND FREQUENCY IN PERCENT BY STABILITY CLASS FOR',
1' EACH SECTOR' »/
-------
FO^TRftN IV G LEVEL 21
MAIN
OATF = 76069
16/52/04
PAGE 0003
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
007S
0076
0077
0078
0079
0080
0081
0032
0063
0084
0085
0086
0087
0088
00fl9
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
2 IX,-DIR',8X,'A',9X,,9X,'C',9X,•D•,9X,'E',9X,'F«>
DO 10501=1,16 ' .
WRITE (6.2.40) DIR( I) , (FREQC I,J),J=1,6)
240 FORMAT(1X,A4,6F10.2)
1050 CONTINUE
c READ/WRITE WIND SPEFOS RY STABILITY CLASS
00 10601=1,16
READ(5,220)(SPEEO(I.J),J=1,6>
1060 CONTINUE
WRITE(6,250)
250 FORMAT(1HO,'WIND SPEEDS IN METERS PER SECOND BY STABILITY CLASS FO
1R EACH SECTOR',/ '
1 IX,'DIR'»8X,'A',9X,'B',9X,'C',9X,«D',9X,«E',9X,'FM
00 10701=1,16
WRITE(6,240) DIR(I),(SPEED(I,J),J=1,6)
1070 CONTINUE
READ(5,260)SIGMAX
260 FORMAT(F10.0)
WRITE(6,?7o)SIGMAX
270 FORMATf1HO.'SIGMAX='.F10.0,• METERS')
C READ KPROB.IF KPR.OB =1,ONLY DOSE CALCULATIONS FOR SPECIFIC
C RECEPTORS WILL BE PERFORMED. IF KPROB =2,ONLY DOSE CALCULATIONS
C FOR THE POPULATION WHEEL WILL BE PERFORMED. IF KPROB IS ANY
C OTHER NUMBER,BOTH CALCULATIONS WILL BE PERFORMED.
READ15.280) KPROB
280 FORMAT(Il)
!F(KPR03.EQ.2) GO TO 1150 '
C READ/WRITE. NUMBER' OF SPECIFIC RECEPTORS
READ(5,290)INO
290 FORMATII3)
IF( INO.LT.1.0R.IND.GT.6) GO TO'l26S
WRITF(6,300>INO
300 FORMAT(1H1,'THERE ARE »,I3,' SPECIFIC RECEPTORS')
C READ/WRITE RADIAL'- AND ANGULAR POSITIONS OF
c SPECIFIC RECEPTORS
C RADIAL POSITION IS THE RADIAL DISTANCE,IN METERS,
C FROM THE CENTER OF THE AREA SOURCE TO THE RECEPTOR POSITION
C ANGULAR POSITION IS THE CLOCKWISE ANGULAR DISPLACEMENT OF THE
c RECEPTOR POSITION,IN DEGREES,FROM THE POSITIVE Y-AXIS
READ(5,310)(RPOS(I),1=1,IND)
READ(5,310)(APOS(I).1=1,IND)
310 FORMAT(6F10.0)
WRITE(6,320)
320 FORMAT(1HO,'RECEPTOR NUMBER',T30,'RADIAL POSITION.METERS',T60,'ANG
1ULAR POSITION,DEGREES')
DO 10801=1,IND
WRITE(6,330)I»RPOS(I)«APOS .
330 FORMAT < IX »I3fT35».FlO.O,T65,FlO.'0)
1080 CONTINUE
C CONVFRT RFCEPTOP .POSITIONS TO X-Y COORDINATES
no 11001=1.IND
THFTA = PIltAPOS( I ) /ISO.
XPOS(I)=RPOS(I)°SIN(THETA)
YPOS(I)=RPOS(I)"COS(THETA)
1100 CONTINUE
C CALCULATE AND OUTPUT DOSERATES.ETC..
C AT SPECIFIC RECEPTOR LOCATIONS
WRITE<6,380)
-------
FORTRAN IV G LEVEL 21
MAIN
DATE = 76069
16/53/04
PAGE 0004
0104
0105
0106
0107
oioa
0109
0110
on i
0112
01 13
0114
0115
0116
0! 17
0118.
01 1Q
0120
0121
0122
0123
0124
0125
0126
0127
012P
0129
0130
0131
0132
0133
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
380 FORMAT(1HO,'RECEPTOR NUMBER',T35.'CHI/Q',T55,'RAOIONUCLIDE CONCENT
1 RATION',T93,-DOSE RATE',/»T33,(SEC/CU.M.•,T63,'CI/CU.M.',T94,'MREM
2/YRt)
DO 11401=1,IND
CALL DIFUSNIXPO'S ( I) ,YPOS ( I ) ,CHIOQ.KK)
IF(KK.EQ.O) GO 'TO 1130
CONC=CHIOQ»SOURCE
•005E = CONC<>DCF :
WRITE(ft.390)I,CHIOO.CONC.DOSE
390 FORMAT <1X,I3,T30,1PE12.3.'"T60,1PE12.3,T90,OPF10.2)
60 TO 1140
1130 WRITFI6.400) 1
400 FORMAT(IX,13.T20»'RrCEPTOR IS TOO CLOSE TO THE SOURCE AREA'1
1140 CONTINUE
WP.ITF (6,405)
405 FORMAT(////////////////////)
TF(KPROR.EQ.1) GO TO 1270
0 RElOA'SITE THE -NUMRER OF ANNULI(MRAD) IN THE POPULATION WHEEL
c THF NUMBER OF SECTORS is ALWAYS 16
IF (KOROH.NE.2) GO TO 1155
1150 UKITF(6,406)
406 FORMAT(1 HI)
1155 MSEC=)6
READ(5,410) MW10
410 FORM AT(12)
WP I TE (6> '-2C-) '".Stc.Mp.iD
420 FOPM'.T ( I'-JO , ' THERE ARE ',12,' 'SECTORS AMD ',12,' ANNULI IN THE POPU
! I \~ ION r.'HEEL ' )
c REin/wRiiE RADIAL DISTANCES,IN METERS,FROM THE CENTER OF THE
C SOURCE AREA TO THE MIDDLE,AND INNER AND OUTER BOUNDARIES OF THE
C MRftO POPULATION WHEFL ANNULI
WR'TTF (6,430)
430 FORMAT(1HO,'DISTANCES-IN METERS,TO THE MIDDLE,AND INNER AND OUTER
iFOUNcupiES OF POPULATION WHEEL ANNULI•//• ANNULUS NUMBER',T2o»»RMI
20' ,T40, 'RINN' , T60, '^ROUT' )
DO 11601=1,MRAD
REAO(5,440) RMIO(I),RINN(!),ROUT(I)
440'FORMAT(3F10.0)
WRlTE(6v450) I,RMID(D«RINN(t).ROUT(I)
450 FORM&T(1X,T5,I2.T15.F1Q.O,T35«F10.0,T55,F10.C)
1160 CONTINUE
C REoD/rfRITE POPULATT'N IN SECTOR SEGMENTS BOUNDED 8Y THE POPULATION
C WHEEL ANNULI. POP(I.J) IS THE POPULATION IN ANNULUS I AND SECTOR J
00 1170I = 1«MRAD-
RFAD(5,460) (PQPd.J) »J=i »8> .
PFAD(5,460) (POP(I.J) »J = 9«16)'
460 FORMAT(BF10.0)
1170 CONTINUE
WRITE(6.470)
470 FORMAT(1H1,'POPULATIONS IN SECTOR SEGMENTS.TWO LINES PER RADIUS',
1'IREAD CLOCKWISE). AND TOTAL POPULATION IN ANNULI'/
21X,< I•,IX,'N/S/TOTAL'«3X»
3'NNE/SSW',5X,'NE/SW',3X,'ENE/WSW',7X,
4'E/W',3X,'ESE/WSW',5X, 'SE/NW'. ,3X, 'SSE/NNW' )'
DO 11901=1,MRAO-
W RITE (6,480) Iv-(POP(I,J) ,J = 1,8)
480 FORMAT (IX,12,T5.,-8F10.0)
WRIT£(6,490) (POPd.J) tJ = 9»16J
-------
FORTRAN IV G LEVEL 21
MAIN
DATE = 76069
16/52/04
PAGE 0005
00
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0150
0160
0161
016?
0163
0164
0165
0166
0167
0168
016Q
0170
0171
017?
0173
017^
0175
0176
0177
0173
0179
0180
0131
0182
0183
OU',4
01H5
0186
0137
o i .-;A
Olq.<3
0191
0191
0192
0194
490 FORMAT(1X,T5.8F10.0>
SUM=0.
00 1180J=1»16
SUM=SUM+POP(ItJ)
1180 CONTINUE
WRITE(6,500> SUM.
500 FORM4T(lX,T5tF10.t>>
1190 CONTINUE '. •_
C NOW INDIVIDUAL AND COPULATION DOSE CALCULATIONS
C WILL BE PERFORMED
DO -1200 J=l,16
ANGLE = 22.5<>FLOATtJ-l>*PI/180.
SINUS=SlN(ANGLei
COSIN=COS(ANGLE>
00 1200 I=l,MrtAD
C CALCULATE THE X-Y COORDINATES•OF THE CENTER OF
C POPULATION WHEEL SECTOR SEGMENT
X=RMIO(I)°SINUS
Y=RMID(I)°COSIN
CALL DIFUSN(X,Y,CHIOQ»KK)
IF(KK.FO.O) GO TO 1260
CHIO'?ITE(6.530) OIRfj).(CHIO(I,J).I = l.MRAO)
530 FORMAT(1X,A4,1P12tl0.2)
lain CONTINUE
WHITE (6.540)
-------
FORTRAN IV G LEVEL ?1
MAIN
PATE-= 76ofe9
16/52/04
PAGE'0006
0195
0196
0197
0198
0199
0200
0201
0202
0203
020*
1250 WRITE(6,580)
5flO FORMAT (1H1 »
59".FORMAT(1H1.'ERROR.ONE SECTOR SEGMENT TOO CLOSE TO SOURCE AREA.PROG
1RAM TERMINATED.')
GO TO 1270
1265 WRITE(6.600)
600 FORM4K1H1,'ERROR.N'JMBER OF SPECIFIC RECEPTORS OUTSICE OF ALLOWED
1RANGE (1-6).PROGRAM TERMINATED.•i
1270 STOP
END
-------
FORTRAN IV G LEVEL 21
DIFUSN
DATE = 76012
14/57/19
PAGE 0001
0001
0002
Ln
O
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
SUBROUTINE DIFUSN(XX,YY,TCHIOQ*KK>
COMMON NX,NY,bIGMAX,SX<10,10),SY<10»10>,FREQ<16,6).SPEED(16.6),ZH,
1 PI.DEC.'OTEST.KDIF
C THIS SUBROUTINE PERFORMS DIFFUSION CALCULATIONS
C THE CHIOQ CONTRIBUTION FROM EACH SOURCE AREA ELEMENT IS
C ADDED TO OBTAIN THE TOTAL CHIOQ AT THE RECEPTOR POSITION
C THE DIFFUSION CALCULATION UTILIZES EQUATION 3.144 IN METEOROLOGY
C AND ATOMIC ENERGY,A:.,D FUNCTION SIGZ. THE VALUE OF
C SIGZ IS SET TO A MAXIMUM VALUE OF SIGMAX
C THE CHIOQ CONTRIBUTION FROM EACH AREA ELEMENT IS DECAYED
C BY TIME OF FLIGHT
C THE TOTAL RADIONUCLIDE CONCENI RAT ION DISPERSION FACTOR,TCHIOQt
c AT A RECEPTOR LOCATION,is NORMALIZED WITH RESPECT TO THE TOTAL
C SOURCE AREA BY DIVIDING IT BY THE NUMBER OF AREA ELEMENTS,(NX«NY)t
C AND PRINTED OUT. THEN,IT IS MULTIPLIED BY THE TOTAL SOURCE AREA
C RADIONUCLIDE RELEASE RATE AND DOSE CONVERSION FACTOR TO
C YIELD A YEARLY DOSE .
C THIS PROCEDURE IS VALID ONLY /<2.«PI>
CONST=0.0203
C THE ANGULAR WIDTH OF A SECTOR,WSA,IS 22.5 DEGREES.
WSA=22.5
II=NY
JJ = NX
IF (KDIF.EQ.l) GO TO 500
C THE FOLLOWING STATEMENTS DETERMINE WHETHER THE AREA
C CAN BE TREATED AS A POINT SOURCE.
RTEST=SQRT(XX°XX»YY»YY)
WSR=PI«WSA/(lbO.«2.)
D=DTEST«(1.*1./TAN(WSR>)/2.
IF(RTEST.GT.D) 11=1
IF(RTEST.GT.D) JJ=1
C INITIALIZE THE TOTAL CHIOQ TO ZERO
500 TCHIOQ=0.
DO 6000 1=1,11
DO 6000 J=1,JJ
C DETERMINE THE DISTANCE FROM TH£ (I,J) AREA ELEMENT TO THE
C RECEPTOR LOCATION
DIFX=XX-SX(I,J)
DIFY=YY-SY(I,J)
IF(KOIF.EQ.I) GO TO 600
IF(RTEST.GT.D) DIFX=XX
IF(RTEST.GT.D) DIFY=YY
600 R=SURT
-------
FORTRAN IV G LEVEL 21
DIFUSIM
DATE * 76012
14/57/19
PAbt 0003
0029
0030
0031
0032
0033
0034
0035
0036
OOJ7
0038
0039
0040
00*1
0042
0043
0044
0045
0046
0047
0048
0049
0050
OOal
0052
0053
0054
0055
C NOW THE-T HUNS FROM 0 TO 360 DEGREES
C DETERMINE IN WHICH SECTOR THET FALLS
C KSEC IS THE SECTOR NUMBER
THET=THET*WSA/2.
KSEC=THET/WbA«l.
C SUM CHIOQ OVER STABILITY KLASSES
C CONVERT RECEPTOR DISTANCE FROM METERS TO KILOMETERS
RKM=0.001aR
CHIOQ=0.
C SUM OVER STABILITY CLASSES
DO 5000 KLASS=1,6
C TO AVOID UNDERFLOWS DUE TO RAPID DECAY OF ShORTLIVED RADIONUCLIDES
C THE CHIOQ WILL BE SET TO ZERO IF THE RADIONUCLIDE HAS DECAYED
C BY MORE THAN JO MEAN LIVES.
C THE CHIOQ WILL ALSO BE SET EQUAL TO ZERO IF FREQ OR SPEED ARE LESS
C THAN 0.01
IF(FREQ 4700.4700.4600
4600 XLIVES=DEC°R/SPEED
IF(30.-XLIVES) 4700,4700,4800
4700 B=0.
GO TO 4900
4800 IF(SIGMAX.LT.l.) SIGMAX=1,E«04
SZ=SIGZ(KLASS.HKM)
IF(SZ.GT.5IGMAX) SZ=SIGMAX
DECAY=EXP(-XLIVES>
ARG = ZH*ZH/(2.°SZ<1SZ)
IF(ARG.GT.100.) ARG=100.
DENOM=SZ°SPEED°R
B=CONST°FREQ(KSEC,KLASS)'DECAY*EXP(-ARG)/DENOM
4900 CHIOQ=CHIOQ*B
5000 CONTINUE
TCHIOQ=TCHIOQ»CHIOQ
6000 CONTINUE
TCHIOQ=TCHIOQ/FLOAT(II*JJ)
KK = 1
7000 RETURN
END
-------
FORTRAN IV G LEVEL 21
GEOMO
DATE = 76013
PAGE 0001
0001
0002
0003
000^
0005
OOU6
0007
0008
0009
0010
0011
SUBROUTINE GEUMO(XL.YL)
C THIS SUBROUTINE DETERMINES THE X-Y COORDINATES OF THE CENTERS OF
C AREA ELEMENTS.THERE A*E NX BY NY ELEMENTS ARRANGED IN A MATRIX OF
C NY ROWS AND NX COLUMNS.ELEMENT (1,1) IS IN THE NORTHWESTERN CORNER
C OF THE MATRIX, AND ELEMENT (NY,NX) IS IN THE SOUTHEASTERN CORNER.
C SX(I.J) AND SY(I,J) ARE,RESPECTIVELY,THE X AND Y COORDINATES
C OF THE CENTER OF AREA ELEMENT (I,J)
C ALL AREA ELEMENTS HAVE THE SAME AREA
COMMON NX,NY,SIGMAX,SX(10»10).SY(10»10>,FREQ<16«6).SPEED(16,6),ZH,
1 PI.DEC.DTEST.KDIF
DELTX=XL/FLOAT(NX)
DELTY=YL/FLOAT(NY>
DO 600 1=1,NY
DO 600 J=1,NX
SY(I,J)=YL/2.-OELTY»(2»I-l)/2.
SX(I,J)=-XL/2.»DELTX° <2»
600 CONTINUE
RETURN
END
-------
FORTRAN IV 6 LEVEL 21
GEOM1
DATE = 76012
14/57/19
PAGE 0001
0001
Ui
LO
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
00^2
0023
0024
0025
0026
0027
0028
0029
0030
SUBROUTINE GEOM1(RL)
C THIS SUBROUTINE DETERMINES THE X-Y COORDINATES OF CENTERS OF AREA
C ELEMENTS INTO WHICH THE CIRCULAR PILE HAS BEEN DIVIDED.RL IS THE
c RADIUS OF THE PILE.NX is THE NUMBER OF ANNULI AND NY is THE NUMBER
C OF SECTORS.SXII.J) AND SY(I.J) ARE.RESPECT IVELY,THE X AND Y
C COORDINATES OF THE AREA ELEMENT(I,J).I DENOTES THE ANNULUS NUMBER.
C STARTING WITH 1 NEAR THE CENTER AND INCREASING OUTWARD,AND J
C DENOTES THE SECTOR WHICH SPANS THE ANGLE FROM
R(I)=SQRT(ARG)
10 CONTINUE
DO 20 1=1.NX
RM(I) = (R(I) *H(I.+ 1) )/2.
20 CONTINUE
DELTR=2.°PI/YNY
DO 30 1=1,NY
THET(I)=(2»I-l)°DELTR/2.
30 CONTINUE
DO 35 1=1.NX
DO 35 J=1,NY
SX(I,J)=RM(I)«SIN(THET(J))
SY(I,J)=RM(I)*COS(THET(J>)
35 CONTINUE
GO TO 45
40 SXU ,1)=0.
SY(1,11=0.
45 RETURN
END
-------
FORTRAN IV 0 LEVEL 21
SIGZ
DATE = 76012
14/57/19
PAGE 0001
0001
0002
OOU3
0004
0005
0006
0007
0008
000.9
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
OOiO
0021
00^2
0023
0024
0025
0026
100
110
120
130
FUNCTION SlbZ IKLASS,X>
FUNCTIONS Ant OAP/EJA FITS TO SIGMA ZEES
IN EPA OFFICE OF AIR PROGRAMS DOCUMENT NO. AP-26. (USGPO » WASH. ,DC.
STOCK NO. 5503-0015, PRICE ONE DOLLAR)
FUNCTION CALCULATES THE STANDARD DEVIATION OF PLUME CONCENTRATION
IN THE Z DIRECTION FOR STABILITY KLASS ( 1 = A, 2 = B, 3 = C,
4 = D, 5 = Et 6 = F) AND THE DISTANCE ALONG THE CENTERLINE OF
THE PLUME IN KILOMETERS, x.
THE STANDARD DEVIATION IS IN DETERS.
AS WRITTEN, A CAPPING LAYER ISIGMAX) UPPER BOUNDS SIGZ
THE RUN IS STOPPED IF X (MM) IS LESS THAN 100. METERS
IF(X-O.I) 110,130,130
140
150
160
170
180
190
200
210
230
220
FORMAT (IX, "ONE DISTANCE R IS TOO SMALL')
STOP
XX=ALOG(X>
SIGMAX=10000.
GO TO (140, 170, 180, 190, 20 0,210) , KLASS
IFIX-1, 5) 150,150,160
SIGZ=EXP(6.126786+XXO(2.214445.XX«(-0.041129*XX»
1 (-0.379863+ XX* (-0.099597) ) ) ) )
IF (SIGZ-SIGMAX) 220,220,230
IF (SIGZ-SIGMAX) 220,220,230
S I GZ = EXP ( 4. 686302*X X«<1. 062550 »XX»( 0.0 18771) ) )
IF (SIGZ-SIGMAX) 220 .220,230
SIGZ=61.141032*X*«. 91*651
IF (SIGZ-SIGMAX 1220,220,230
SIGZ = EXP<3.416367*Xxf>(0.729577*XX*(-0.031207) ) »
IF (SIGZ-SIGMAX) 220, 220,230
S I GZ = EXP( 3. 057629*XX«(0.67908y»XX<*< -0.044892) ) )
IF (SIGZ-SIGMAX 1220,220,230
SIGZ=EXP(2.62i>488»Xx«<0.658866«XX*<-0.054137>-M
IF (SIGZ-SIGMAX) 220.220,230
SIGZ=SIGMAX
RETURN
END
-------
REFERENCES
(I) MARTIN, J. A., C. B. NELSON, and P. A. CUNY. AIREM Program
Manual - A computer code for calculating doses, population
doses, and ground depositions due to atmospheric emissions of
radionuclides. Office of Radiation Programs, Environmental
Protection Agency, Washington, D.C. 20460 (May 1974).
(_2) SLADE, D. H., editor, Meteorology and Atomic Energy, 1968,
TID-24190. National Technical Information Service, U.S.
Department of Commerce, Springfield, Virginia 22151 (July 1968),
(_3) TURNER, D. BRUCE. Workbook of atmospheric dispersion estimates,
Publication No. AP-26. U.S. Environmental Protection Agency,
Office of Air Programs, Research Triangle Park, North
Garolina 27711 (July 1971).
(_4) HOLZWORTH, G. C. Mixing heights, wind speeds and potential for
urban air pollution throughout the contiguous United States,
Publication No. AP-101. U.S. Environmental Protection Agency,
Office of Air Programs, Research Triangle Park, North
Carolina 27711.
(5) ENVIRONMENTAL PROTECTION AGENCY. Environmental Analysis of the
Uranium Fuel Cycle, Part I - Fuel Supply. Office of Radiation
Programs, Washington, D.C. 20460 (October 1973).
(6.) THOMAS, G.B., Jr. Calculus and Analytic Geometry, Fourth
Edition, Addison-Wesley Publishing Company, p. 183 (1968).
55