United States
Environmental Protection
Agency
Corvallis Environmental
Research Laboratory
Corvallis, Oregon 97330
PREDICTION OF INITIAL MIXING FOR
MUNICIPAL OCEAN DISCHARGES
by A. M. Teeter and D. J. Baumgartner
CERL - 043

-------
PREDICTION OF INITIAL MIXING FOR
MUNICIPAL OCEAN DISCHARGES
by A. M. Teeter and D. J. Baumgartner
CERL - 043
Marine Systems Division
Corvallis Environmental Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
200 S.W. 35th Street
Corvallis, Oregon 97330

-------
PREDICTION OF INITIAL MIXING FOR
MUNICIPAL OCEAN DISCHARGES
by
A. M. Teeter and D. J. Baumgartner
Marine Systems Division
Corvallis Environmental Research Laboratory
Corvallis, Oregon 97330
CERL PUBLICATION 043
May, 1979
CORVALLIS ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
CORVALLIS, OREGON 97330

-------
CONTENTS
Page
Abbreviations, Symbols, etc	 i"11
Acknowledgments	.		v
1.	Introduction
Mixing Zone Concepts . . 			1
Ocean Discharges			2
Buoyant Plume Models 		4
2.	Recommendations
Methods		7
Appropriate Conditions 			7
Mixing Zone Specification. 		
3.	Summary and Conclusions
General. . .		13
Simplified Results for Plumes in Stratified Ambients ....	15
Plumes in a Stagnant Environment		16
Plumes in a Flowing Environment . 			17
Merging Plumes in a Stagnant Environment		19
Merging Plumes in a Flowing Environment 		20
Non-Uniform Stratification			22
Plumes in Uniform Density Ambients 	 ....	22
4.	PLUME Model
Theoretical Development			26
Model Description		29
Example Input, Output and Model Listing. . . 		30
5.	OUTPLM Model
Theoretical Development and Model Description		41
Example Input, Output and Model Listing	42
6.	DKHPLM Model
Theoretical Development. ....... 			51
Model Description. . 					54
Example Input, Output and Model Listing		 . .	55
References				85
ii

-------
LIST OF ABBREVIATIONS, SYMBOLS, ETC.
B - Radius or characteristic half-width of the plume.
C - Concentration by volume, = 1/Sa.
Ca - Concentration in ambient.
Ce - Concentration in effluent.
Cf - Concentration after initial dilution.
D - Diameter of discharge port.
E - Entrainment.
1/2
Fn - Froude number, = U./(g D A0/po)
J
g - Gravity.
Q - Volume rate of discharge.
q - Volume rate of discharge per unit length of diffuser.
r - Radial distance from the plume's centerline.
s - Distance along plume axis.
Sa - Flux averaged dilution or the total plume efflux divided by the volume
of the discharge.
sg - Length of the zone of flow establishment.
T - Temperature.
T' - Fluctuating component of T.
U - Ambient current speed.
U.	-	Initial velocity of the jet.
3
u'	-	Fluctuating component of V.
V	-	Velocity at some point in the plume in the s direction,
v	-	Velocity normal to the plume's axis.
i i i

-------
v' - Fluctuating component of v.
x - Horizontal distance in the direction of a jet discharge and ambient
currents or the horizontal distance parallel to a diffuser.
y - Horizontal distance 90 to x.
z - Vertical distance from the discharge depth.
2 - Depth of water.
A - Normalized density disparity, = (p^ - p)/(pQ " Pd).
Ao ~ Initial density disparity of the waste, = p0 - p^.
Ah - Maximum height of rise or height of rise to the trapping level.
p - Density of the plume at some point.
p0 - Ambient density at the level of the discharge.
Pd - Density of the discharge.
Pw - Ambient density at some level.
dp/dz - Absolute value of the ambient density gradient.
1 m3 = 35.31 ft3 = 264 gallons.
iv

-------
ACKNOWLEDGMENTS
The assistance and cooperation of Drs. Lorin Davis, Larry Winiarski, and
Safa Shirazi, and Mr. Walter Frick of the Ecosystem Modeling and Assessment
Branch of the U.S. Environmental Protection Agency's Corvallis Environmental
Research Laboratory are gratefully acknowledged. A similar acknowledgment is
due Bill McDougal of the Marine and Freshwater Branch.
v

-------
SECTION. 1
INTRODUCTION
Initial dilution is the rapid turbulent mixing between wastewater dis-
charged at depth in the ocean and seawater. It takes place close to the point
of release, resulting from physical interactions between these components. The
understanding of this process can be used to predict the initial dilution
attained under given conditions. Such information enables comparison of waste-
water constituent concentrations, subsequent to initial dilution, to water
quality standards and thus can assist in evaluating the discharge's effect on
the marine environment.
The purpose of this repprt is to present procedures for calculating the
initial dilution and describing the zone of initial dilution surrounding a
discharge site. These calculations would be required under regulations pro-
posed by the U.S. Environmental Protection Agency on April 25, 1978 (1) imple-
menting section 301(h) of the Clean Water Act (2), draft regulation implement-
ing section 403(c), and provisions of the State of California's Ocean Plan (3).
MIXING ZONE CONCEPTS
When wastewater is injected into the ocean, it has buoyancy and momentum.
This energy is used in an initial mixing or dilution process. The wastewater
discharge forms a plume which rises and entrains ambient water and is thus
diluted. The mixing process slows markedly when plumes reach a position of
equal buoyancy and momentum with the ambient (no relative energy) or the
surface. Then the dilution process becomes more dependent on ambient oceanic
1

-------
processes. Thus, the mixing zone is a region where rapid mixing takes place
between waste and ambient waters. It can be thought of either as a concentra-
tion isopleth or as the zone where initial mixing takes place. It is not a
"treatment for wastes but a transition region (4). By insuring that water
quality standards are not violated outside this zone, impacts on the marine
environment from the waste discharge should be controlled. While the dilution
of some pollutants may have little effect on their impacts, biological uptake,
concentration, and deposition, and physical processes, like aggregation and
sedimentation, are concentration dependent.
OCEAN DISCHARGES
Ocean discharges of municipal wastes have a wide variety of physical
characteristics which can affect initial dilution. Although single port out-
falls are still used, many outfalls now are multi-port designs to optimize the
dilution and vertical rise attained by their plumes. The numbers of ports and
their orientation can form complex patterns reflecting a great deal of design
effort. The discharged material is almost uniform, with physical characteris-
tics similar to those of fresh water. The marine environment, however, varies
greatly along the miles of our country's shoreline and the physical conditions
at the discharge site have a great effect on initial dilution. The ambient
receiving waters are of a higher density than the waste discharge and are
typically vertically stratified. Currents are usually present but may have
periodicity caused by tidal and other forces.
Discharges initially have momentum and buoyancy as they are injected into
1 /2
the environment. Densimetric Froude numbers (Fn = U./(g D A0/po) ) describe
J
the ratio between these two forces. The higher the initial Froude number, the
more closely the resulting plume resembles a momentum or forced jet. The
2

-------
smaller the number, the more the plume resembles a purely buoyant element.
Because the density difference between the waste and the ambient varies only
slightly and jet velocity is bounded by practical considerations, Froude numb-
ers for ocean discharges generally vary between 0 and about 30. The ends of
this range are characteristic of single-port and multi-port diffusers, respec-
tively. In this range of Froude numbers, buoyancy is likely to dominate the
initial mixing process, making volume rate of discharge, currents, ambient
stratification, and possibly water depth, important. Especially in a current,
buoyancy dominates since the plume's momentum becomes insignificant compared to
the momentum entrained from the ambient.
Existing coastal discharge sites vary in depth up to about 75 meters. The
water columns at these sites are typically stratified vertically by differences
in temperature and/or salinity. This leads to density stratification. The
amount or degree of stratification varies widely in time and space. Proximity
of large sources of fresh water, surface heating, advection, and wind mixing
all can influence density stratification. Salt wedges can form near the mouth
of estuaries with large fresh-water flows. Advection and nonlinear interac-
tions between wind-induced surface turbulence and buoyancy gradients also can
produce pycnoclines, areas of strong density gradients or discontinuities.
Strong pycnoclines can arrest the vertical motion of buoyant plumes, ending the
initial mixing process. Winds of sufficient duration can erode the density
gradient by the turbulence they produce, leaving the ambient well mixed and
allowing buoyant plumes to rise to the surface.
Currents are very complex in coastal areas. These are areas of high
energy received from a number of processes. Currents at discharge sites have
tidal components which can be strong. Oceanic currents impinge on the coast
and intensify on the western edges of the oceans (the Florida Current for
3

-------
instance). Local forces such as wind shear and waves also generate currents.
Instantaneous current values (time scale of minutes), rather than long-term
(time scale of hours, days, etc.) net currents, are important in the initial
dilution process.
BUOYANT PLUME MODELS
Ocean outfalls can be modeled by properly scaled hydraulic models, or by
mathematical models. The methods described here include only the latter,
although in cases of complex geometry or special conditions, physical or hy-
draulic modeling may be appropriate. The characteristics and behavior of
plumes and jets in other aquatic settings, the atmosphere, as well as for
ocean outfalls, have been studied by means of mathematical models. Many
excellent reviews have been published (5, 6, 7).
Early work in this field was concerned with purely buoyant convection (no
initial momentum). Rouse et aL (8) studied the buoyant plume above a contin-
uous heat source in a uniform, stagnant environment. Equations of mass contin-
uity, momentum, and energy were integrated by assuming Gaussian approximations
for the lateral distributions of velocity and temperature across the plume and
with experimentally determined spreading constants. The results were used to
describe the distribution of velocity and temperature as functions of mass
flow and height of rise. Priestly and Ball (9) studied thermal plumes issuing
into a stratified, as well as uniform, environment. They used essentially the
same equations as Rouse et aL for mass, momentum, and energy, but their
lateral Gaussian profiles of temperature and velocity assumed a lateral length
scale.
Morton et a2. (10) studied buoyant plume convection. They introduced an
entrainment function based on local characteristic plume velocity and length
4

-------
scale to replace the spreading ratio. This approach has been adapted by many
later investigators (5). Later, Morton (11) assumed other lateral profiles of
velocity and temperature to account for the difference between eddy transport
of heat and momentum. Another coefficient was introduced and estimated from
experimental data.
Abraham (12) investigated neutrally-buoyant forced plumes, purely-buoyant
plumes, and plumes with combinations of both properties. Equations of mass,
momentum, and energy were integrated using lateral Gaussian profiles recogniz-
ing the difference in transfer rates between momentum and heat (or concentra-
tion). Entrainment was considered by integrating the continuity of mass equa-
tion with appropriate spreading coefficients. Abraham's results included
analytic expressions for plume centerline velocity and concentration as func-
tions of initial plume characteristics, height of rise, and entrainment coef-
ficients for momentum and heat. The limiting cases of neutrally-buoyant
forced jets and purely-buoyant plumes were considered.
Fan (13) reported on plumes issuing at arbitrary angles into a stagnant,
stratified environment. Velocity and buoyancy profiles were used similar to
those used by Morton (11) and included a lateral characteristic length. Solu-
tions were found by simultaneously solving six equations (continuity of mass,
vertical momentum, horizontal momentum, conservation of buoyancy, vertical
geometry, and horizontal geometry).
Fan (13) also investigated buoyant plumes in a flowing, uniform-density
ambient. A theoretical solution was given as a function of plume drag forces
and an entrainment function. The entrainment function was such that drag and
entrainment coefficients varied with plume conditions, a disadvantage in
applications beyond the range of experimental validation (14). Abraham (14)
introduced an entrainment function based on subdomains influenced by initial
5

-------
momentum and buoyancy. With entrainment coefficients defined for these two
subdomains and a single coefficient describing the drag force, the theoretical
solution was in good agreement with the experimental results of Fan.
Plumes from multi-port discharges have been studied less extensively.
Cederwall (15) studied the flow regimes of line sources for discharges in
confined and unconfined environments. Roberts (16) studied the flow regimes
of line sources in unstratified and unconfined environments. Soti1 (17)
proposed a mathematical model for slot jets or continuous line sources in
stagnant, stratified ambients. Koh and Fan (18) proposed a mathematical model
of multi-port diffusers by interfacing single jet and slot jet solutions at a
transition point. A good review of this subject is given by Davis (7).
Three computer models are described in this paper. The PLUME model by
Baumgartner et aK (19) simulates a solitary plume in a stagnant environment
of arbitrary stratification. OUTPLM by Winiarski and Frick (20) models a
single plume in a stagnant or flowing environment of arbitrary stratification.
DKHPLM by Davis (21) describes single plumes in flowing, stratified ambients
which are allowed to merge with identical adjacent plumes to simulate closely
spaced, multiple-port diffusers. The results of these models have been used
to develop analytic solutions to describe the initial dilution of ocean dis-
charges in terms of the principal quantities involved. These three computer
models and the analytic solutions resulting from them will be presented along
with suggestions and example applications.
6

-------
SECTION 2
RECOMMENDATIONS
METHODS
The buoyant-plume phase of waste dispersion can be described by a variety
of published methods. However, to evaluate the water quality impacts of
municipal ocean discharges, the methods presented here are well suited and
recommended for consistency. The analytic methods or computer models can be
used as the situation warrants or requires. When especially high resolution
is necessary, the computer models should be used. For unusual situations or
conditions, other methods will be considered.
Although these methods are considered reliable, they should be continu-
ally evaluated relative to theoretical developments and measurements, espec-
ially quality field measurements.
APPROPRIATE CONDITIONS
The dilution achieved during the initial mixing process is dependent on
ambient and discharge conditions and is, therefore, highly variable. To avoid
the inception of biological effects, the occurrence of pollutant concentra-
tions greater than water quality criteria should be small. In evaluating a
discharge's effect on water quality, therefore, the appropriate conditions to
consider are those which result in the "lowest" dilution and those which occur
at times when the environment is most sensitive. For example, minimum dilu-
tion will be predicted using a combination of maximum vertical density strati-
fication, maximum waste flow rate, and minimum currents for a particular site.
7

-------
Other situations may be more critical to water quality depending on the ambi-
ent water quality, and applicable water quality standards.
Predicting "lowest" dilution relies on the availability of quality meas-
urements or data with which to estimate absolute worst case conditions, and is
therefore not without problems. The statistical uncertainty in estimates of
absolute worst case conditions is generally great. Also there are inherent
biases to some oceanographic measurements. For example, current measuring
instruments have finite thresholds. It therefore becomes difficult to dis-
tinguish low values (which may be as high as 0.05 m/sec) from zeroes with
these data. In estimating environmental conditions, a more reliable estima-
tion can be made at the lowest 10 percentile on a cumulative frequency distri-
bution. Data on ambient density structure are not collected by remote means
so that the frequency of observations is generally low. To increase the
reliability of "worst-case" estimates, data should be evaluated not only from
the discharge site but from nearby coastal areas of similar environmental
setting.
Defining "worst-case" conditions as a combination of affecting conditions
each taken at the.worst 10 percentile on cumulative frequency distributions is
recommended. This approach allows a reliable estimation of these conditions
to be made and prevents the unlikely occurrence of more extreme conditions
from biasing the predictions. The probability of these conditions occurring
simultaneously is much less than 10 percent insuring that the predicted dilu-
tion will be exceeded most of the time. Given the sensitivity of dilution to
external conditions and the ranges of these conditions, the evaluation of data
and selection of model input conditions are critical to the proper evaluation
of water quality impacts from municipal ocean discharges.
The characteristics of the discharge needed for evaluation are:
8

-------
a)	port sizes, spacings, depths and orientation,
b)	flow rates through ports,
c)	wastewater density, temperature and apparent salinity at the point
of discharge.
The hydraulic characteristics of diffusers should be checked and the
diffuser segment with the highest flowrate per unit diffuser length used for
subsequent initial dilution predictions. Apparent salinity refers to the
major inorganic ions of the wastewater which contribute to its density. It
has been found that a salinity, based on conductivity ratios for seawater,
gives an acceptable value. Alternately, an apparent salinity which corres-
ponds to a measured density can be used. Apparent salinities of municipal
effluent generally range from 1 to 5 /oo.
As mentioned, the principal quantities to consider in dilution prediction
are volume rate of discharge (Q) or volume rate of discharge per unit diffuser
length (q), ambient stratification (dp/dz), initial density difference between
waste and ambient (A0), and ambient currents (U). These quantities should be
considered at:
a)	times of seasonal maximum waste outflow,
b)	times of seasonal maximum and minimum stratification,
c)	times of low ambient water quality,
d)	times of low net circulation or flushing, and
e)	times of exceptional biological activity.
These periods have durations generally of 1 to 3 months. The quantities se-
lected to represent these periods should reflect conditions of the worst 10
percent of occurrences during these periods; highest volume rate of dis-
charge, extreme ambient stratification, lowest initial density difference, and
9

-------
lowest current speed. If a reliable estimate can't be made from the available
data, a conservative estimate should be substituted, e.g. zero currents.
To select conditions at the worst 10 percentile, they must be ranked in
some way. Current speed data usually consist of discrete values and can be
ranked into cumulative frequency distributions by computer or by hand. It is
more difficult to rank existing data on ambient stratification and discharge
rate. A suitable value to use for the volume rate of discharge is that value
exceeded 2 to 3 hours during the average day for the period. The average over
the period of daily values exceeded 2 to 3 hours per day or the value exceeded
2 to 3 hours on a day which has average total flow over the period could be
used, depending on the form of the flow data. Data from the last 2 years or
longer should be used to insure that the flows are representative.
The worst stratification is that which is greatest over the height of
rise of the plume. This condition will result in the lowest dilution if other
conditions are constant. If the density gradients are uniform with depth they
can be ranked numerically. If not, measured profiles can be input to the
computer models PLUME or OUTPLM and the results used to rank the profiles.
Alternatively, the average or overall density gradient over the height of rise
of the plume can be estimated by the method given in the next Section (Non-
Uniform Stratification) and these average gradients ranked. The density dif-
ference between the ambient and the effluent will probably vary only a small
amount and will probably be due to changes in the ambient density. The dens-
ity difference can, therefore, be selected by ranking the ambient density at
the level of the discharge.
In unusual cases, the discharge ports may resemble thin (< 1/10 D),
sharp-edged orifices. In these cases, the port diameter may be taken as 0.62
times the actual diameter of the orifice and the discharge velocity adjusted
10

-------
accordingly. This will increase the computed Froude number of the discharge.
Where the orifice is rounded or thick, as is generally the case, no correction
should be made.
MIXING ZONE SPECIFICATION
After initial dilution, the concentration of waste constituents (Cf) is a
function of the average dilution achieved (Sa) and their concentration in the
ambient (Ca) and the effluent (Ce):
Cf = Ca + (Ce-Ca)/Sa
If the effluent has been adequately treated and disposed of in an environment-
ally appropriate area, Cf's for its various constituents should be at or below
water quality standards.
The zone surrounding the discharge site geometrically bounding the criti-
cal initial dilutions is termed the zone of initial dilution (ZID) to disting-
uish it from other mixing zone definitions. It is a concentration isopleth.
The ZID describes an area in which inhabitants, including the benthos, may be
chronically exposed to concentrations of pollutants in excess of water quality
standards or at least to concentrations greater than those predicted for the
critical conditions described above. The ZID does not attempt to describe the
area bounding the entire initial mixing process for all conditions or the area
impacted by the sedimentation of organic material.
The purpose for defining a zone around the discharge point is to disting-
uish biological impacts caused by exposure of organisms to concentrations
higher than Cf from those impacts caused by the re-concentration of pollutants
by physical, chemical, and biological means. The former is unavoidable with-
out pre-dilution or treatment of the waste to meet water quality standards.
The latter indicates that dilution alone is not sufficient to ameliorate the
11

-------
environmental and ecological impacts of the discharge and additional treat-
ment, relocation of the discharge, or influent source control may be neces-
sary.
The ZID can be specified by analyzing prediction results for a range of
conditions. However, it can be more simply specified using the maximum height
of rise predicted for the critical conditions as a radial distance measured
horizontally from the outfall diffuser or port. This distance will often
equal the depth of water at the discharge site. During periods of higher
currents, the plume trajectory will be more horizontal and initial dilutions
higher than predicted for the critical conditions. The dilution achieved over
that portion of the trajectory within the ZID, however, will be approximately
equal to the initial dilution predicted for the critical conditions.
12

-------
SECTION 3
SUMMARY AND CONCLUSIONS
GENERAL
The plume attributes of interest here are the dilution attained at the
end of convective ascent through ambient waters, and the plume's position.
For our purposes, the poorest or lowest dilutions are important. Dilution is
defined as the total volume of a parcel or sample divided by the volume of
waste contained within it. This is equivalent to the inverse of the volume
fraction (or concentration) of waste that a parcel or sample contains.
The behavior of buoyant plumes can be mathematically modeled by properly
considering mass, momentum, and energy. Some form of entrainment function
must be assumed and fitted to experimental data. Other assumptions generally
made are that all flows are steady and incompressible, pressure is hydrostatic
throughout, the plume iis fully turbulent and axisymmetric, and that turbulent
diffusion dominates and is significant only in the radial direction. Distri-
butions of velocity and concentration may also be assumed. The physics of
this problem can be solved in various ways. Systems of differential equations
can be integrated across the plume, reducing the equations to one independent
variable, length along the axis of the plume. Difference equations can be
used with Lagrangian coordinants. Simultaneous equations or three-dimensional
numerical schemes also have been used.
Mathematical models of buoyant plumes are systems with internal variables
(mass, momentum, and energy), external variables (discharge characteristics,
ambient vertical density, and currents) and internal mechanisms (entrainment).
13

-------
These models have behavior which is affected by principal quantities. Compu-
ter model studies can be used to describe plume behavior of interest in terms
of effecting principal quantities. For stagnant conditions, the principal
quantities are initial density difference between the waste and ambient (A0),
ambient stratification (dp/dz), and flow rate or flow rate per unit length for
multi-port diffusers (Q or q). Current speed (U), if considered, is also a
principal quantity. Froude number is not important to the dilution attained
in the neighborhood considered (0 < Fn < 30), especially in currents. The
Froude number and orientation of the discharge affect the plume's trajectory
in stagnant conditions but have only a small effect in a flowing ambient. The
spacing of multiple ports, for a given discharge rate per diffuser length,
does not significantly effect the dilution attained by merging plumes in a
flowing ambient. In stagnant conditions, port spacing has a moderate (perhaps
20 percent) effect on initial dilution.
Closed-form solutions for the behaviors or attributes of interest (dilu-
tion at maximum height of rise and height of maximum rise) have been found in
semi-theoretical, state-transition modeling. These solutions have been mapped
or tuned to the behavior of the computer models. Plume behavior can, there-
fore, be solved analytically in terms of the principal quantities. These
solutions also allow the importance of the principal quantities to be seen
directly. The computer models are still useful for detailed analysis or spec-
ial conditions. In many cases, however, the analytic models are quite ade-
quate to calculate dilution and maximum height of rise.
The discharge configuations which will be considered include single
ports, and multiple ports spaced such that they interact at some point during
convective ascent. These configuations are termed, respectively, plumes, and
merging plumes. Some other restrictions and assumptions must be made. The
14

-------
discharge site is assumed to be at least ten port diameters	deep and uncon-
fined. Flow in the ambient is assumed to be uniform. Each	of the three
models used in this study is suited to certain combinations	of conditions.
Recommendations are given in Table 1. Application criteria	for the analytic
or simplified models are summarized on Figures la and lb at the end of this
section.
TABLE 1. RECOMMENDED MODEL APPLICATION.
Discharge Ambient Ambient Surfacing Analytic 	Computer Models	
Configuration Density Currents Plumes? Models PLUME OUTPLM DKHPLM
Single Plumes Stratified Stagnant Yes	X	X
"	"No	XX
11 Flowing Yes	X	X
No	X	X
Uniform Stagnant Yes	X	X
11 Flowing Yes	X	X
Merging Plumes Stratified Stagnant Yes	X	X
No	X	X
" Flowing Yes	X	X
11	No	X	X
SIMPLIFIED RESULTS FOR PLUMES IN A STRATIFIED AMBIENT
Analytic solutions are given in this section for dilution and height of
rise of plumes and merging plumes. For stagnant conditions, the height of
rise is taken to the trapping level. The trapping level is where the plume
first loses its positive buoyancy. Above this level the plume is expected to
"pool" and spread horizontally. For flowing conditions, height of rise is
taken to the centerline of the plume where it first becomes horizontal.
Dilutions predictions are made for these locations.
The results described cover ambient conditions including vertical density
stratification, both stagnant and flowing, and the previously described dis-
15

-------
charge configurations. Vertical density stratification is normally found in
the marine environment and, since it adversely affects dilution, should be
considered. A summary of the analytic solutions covering the conditions and
configurations listed in Table 1 will be given. Following these, sections on
non-uniform stratification and single plumes in uniform-density ambients are
presented. SI units (kilograms, meters, seconds) are used.
Plumes in a Stagnant Environment
The dilution for a single port discharge into a stably-stratified and
quiescent ambient is given by::
Sa = 0.46 Q-1/4 A03/4 (dp/dz)"5/8	(1)
where Sa is the flux averaged dilution at the trapping level. The height of
the trapping level can be determined by:
, Ah = 4.6 (A0 Q)1/4 (dp/dz)"3/8	(2)
The dilution of the plume is most affected by initial density difference
(A0) and ambient density gradient (dp/dz), and secondarily by the volume rate
of discharge (Q). Flux of buoyancy and ambient density gradient both affect
height of rise. The height of rise (Ah) given by equation (2) is a maximum
value for the trapping level. This expression does not consider the angle of
discharge from the vertical or the Froude number which can reduce the actual
trapping level height. If, in application, it is found that the height of
rise given by equation (2) exceeds the water depth (Z) at the discharge site,
dilution can be calculated by:
Sa = 0.10 Z (A0/Q)1/2 (dp/dz)"1/4	(3)
This relationship assumes linear dilution with height of rise which is not
strictly observed in stagnant conditions. As the water depth becomes smaller
than Ah, equation (3) predicts dilution which becomes less than those pre-
16.

-------
dieted by the PLUME program. If greater accuracy is needed the PLUME program
should be used.
Equation (1) is in good agreement with model results for jets with Froude
numbers of up to about 50. For jets with Froude numbers above this value the
PLUME program should be used.
Example
A discharge has a maximum flow rate of 0.10 m3/sec and a density of
1000 kg/m3. The ambient at the discharge site has a density at the surface of
1020 kg/m3 and, at 30 m depth, a density of 1025 kg/m3 during periods of maxi-
mum stratification. From equations (1) and (2):
Sa = 0.46 (0.10)"1/4 (25)3/4 (0.166)~5/8
= 28
Ah = 4.6 ((0.10)(25))1/4 (0.166)~3/8
= 11.3 m
If the depth over the discharge were only 10 m and other conditions the same,
then from equation (3):
Sa = 0.10 (10.) (25/0.10)1/2 (25)1/2 (0.166)"1/4
- 25
Plumes in a Flowing Environment
The dilution and height of rise of a single plume discharged into a
stably-stratified, flowing ambient are:
Sa = 3.0 (U/Q)1/3 (A0/(dp/dz))2/3	(4)
Ah = 2.3 (A0 Q/U(dp/dz))1/3	(5)
17

-------
The dilution of the plume is most strongly dependent on the ratio of initial
density difference to ambient density gradient. Equation (4) will predict
increased dilution (greater than equation (1)) when:
U > 0.0036 (AoQ)1/4 (dp/dz)1/8
Neither the dilution nor the height of rise are sensitive to discharge angle
or typical Froude numbers. The maximum height of rise is for the centerline
so that the plume actually extends one radius up from this level. The plume
1/2
radius can be estimated by: B = (SaQ/nU) If Ah + B exceeds the depth of
1/3
water at the discharge site (when Z < 3.3 (A0Q/U(dp/dz)) ), the computed
dilution can be corrected by:
Sa (corrected) = Sa  Z/(Ah + B)	(6)
where B is the plume's radius or half width. By substitution into equation
(6), the dilution of a plume which reaches the surface is:
Sa = 0.92 Z (U/Q)2/3 (Ac/(dp/dz))1/3	(7)
Equation (7) is applicable when U > 0.036 (A0Q)^ (dp/dz)^, otherwise
use equation (3).
Example
In the previous example, A0 = 25 kg/m3, Q = 0.10 m3/sec, dp/dz =
0.166 kg/m3/m. If the lowest 10 percentile of current speed observations is
0.10 m/sec, what is the dilution and height of rise?
Sa = 3(0.10/0.10)1/3 (25/0.166)2/3
= 85
Ah = 2.3((25)(0.10)/(0.10)(0.166))1/3
= 12.2 m.
1 /?
The radius can be estimated as: B = (56(0.10)/n(0.10))
= 5.2 m
18

-------
If the water at the discharge site is only 10 m deep, dilution as cor-
rected by equation (6) is:
Sa (corrected) = 85(10)/(12.4 + 5.2)
= 49
Alternately, dilution can be re-computed using equation (7):
Sa = 0.92 (10) (0.10/0.10)2/3 (25/0.166)1/3
= 49
Merging Plumes in a Stagnant Environment
Brooks (22) gives expressions for the dilution and height of rise of slot
plumes in stagnant, stratified ambients. Slightly re-arranged, re-dimensioned,
and tuned to DKHPLM for consistency with the other analyses presented here,
dilution can be predicted by:
Sa = 0.32 q~1/3A02/3(dp/dzf1/2	(8)
The height of rise of the merging plume is:
Ah = 5.9 (A0q)1/3 (dp/dz)"1/2	(9)
The height of rise is predicted from equation (9) is the maximum trapping level
height.
Widely spaced plumes which merge after a significant rise will attain a
higher dilution (perhaps 20 percent) than predicted by equation (8).
If the height of rise of the plumes exceeds the depth of water, dilution
can be calculated by:
Sa = 0.054 Z q"2/3 A01/3	(10)
This dilution will be lower than predicted by simulation model DKHPLM because
of the assumed linearity of dilution with height of rise and because the trap-
ping level is assumed to be the maximum. If more accuracy is needed, DKHPLM
should be used.
19

-------
Example
A diffuser discharges 0.02 m3/sec/m in 50 m of water. If the efflu-
ent has a density of 1000 kg/m3, A0 = 25 kg/m3, and the ambient density gradi-
ent is 0.05 kg/m3/m, the dilution by equation (8) will be:
Sa = 0.32(0.02)"1/3 (25)2/3 (0.05)~1/2
= 45
The maximum height of rise to the trapping level, from equation (9), will be:
Ah = 5.9 ((0.02)(25))1/3 (0.05)"1/2
= 20.9 m
If the water depth were 10 m, the dilution, from equation (10), would be:
Sa = 0.054 (10) (0.02)~2/3 (25)1/3
= 21
Merging Plumes in Flowing Environment
For multi-port diffusers where plumes interact or for slot plumes, the
principal quantities which affect dilution and height of rise of the plume are
initial density difference (A0), discharge rate per unit length of diffuser
(q), current speed (U) and ambient density gradient (dp/dz). The dilution from
these discharges can be found by:
Sa = 2.7 (A0U/q(dp/dz))1/2	(11)
and the height of rise by:
Ah = 1.3 (A0q/U(dp/dz))1/2	(12)
Like single plumes, the dilution and height of rise of multiple plumes in a
current are not sensitive to Froude number and discharge orientation. Appar-
ently the port spacing or level of interaction also have little effect on the
values ultimately attained by the plume.
Currents will increase the dilution attained by merging plumes when:
20

-------
U > 0.014 (A0q)1/3
If currents are not perpendicular to the diffuser, the discharge rate (q) can
be adjusted to consider the length of the diffuser projected by the current.
Below a relative angle of about 45 degrees, however, prediction should be
based on stagnant conditions.
The height of rise of merging plumes is about equal to the half-width.
1 /2
If Z < 2.6 (A0q/U(dp/dz)) the plume will surface and its dilution, cor-
rected by equation (6), calculated by:
Sa = 1.0 Z U/q	(13)
The dilution of surfacing, merging plumes is affected by currents when:
U > 0.054 (A0q)1/3
These analyses, as well as DKHPLM, assume that the diffuser length is very
long compared to the height of rise. Otherwise edge effects complicate the
problem.
Example
A discharge of 12.5 m3/sec 300 MGD) is made in 50 m of water.
The ambient density during maximum stratification varies from 1024 kg/m3 at
the surface to 1025 kg/m3 at the bottom. The diffuser is 1250 m long so that
the discharge rate per unit length of diffuser (q) is 0.01 m3/sec/m. Density
of the effluent is 1000 kg/m3. If the lowest 10 percentile of current speed
occurrences is 0.05 m/sec, the dilution will be about:
Sa = 2.7((0.05)(25)/(0.01)(0.02))1/2
= 213
And the expected height of rise will be:
Ah = 1.3((25)(0.01)/(0.05)(0.02))1/2
= 21 m
21

-------
If the discharge site were only 20 m deep, equation (13) predicts the dilution
to be:
Sa = 1.0 (20)(0.05)/0.01
= 100
Non-Uniform Stratification
The preceding analyses are for uniform, or monotonical, density gradi-
ents. Often density gradients deduced from field measurements are not
strictly monotonical but can be approximated as such. For cases where the
vertical density distribution is interrupted by a sharp discontinuity or
pycnocline, other approximations are necessary. Brooks (22) gives one method.
The height of rise can first be estimated by using the density gradient of the
lower layer. If this height of rise extends to the upper layer, the density
gradient is taken as the average gradient to that point. It may be necessary
to re-estimate maximum height of rise if the first adjustment resulted in a
large change in the average density gradient. The density gradient can then
be re-estimated. By this method, an average density gradient over the ex-
pected height of rise of the plume can be estimated and dilution calculated.
This procedure is useful for evaluating stratification data and identifying
the most adverse stratification. If the density structure of the ambient is
very complex, or if more accuracy is needed, then PLUME or OUTPLM models
should be used.
PLUMES IN UNIFORM DENSITY AMBIENTS
Conditions may exist at a discharge site which mix the ambient.to a uni-
form state. Normally these are not the conditions of interest since ambient
stratification is more critical to the initial dilution process. The analytic

-------
solutions for single plumes in stratified ambients, because of their form,
approach infinity as the ambient density gradient approaches zero. Thus, they
may over-estimate dilution when stratification is slight. Different analytic
solutions are therefore necessary for single plumes in uniform ambients. The
analytic solutions for merging plumes in stagnant and flowing conditions,
equations (10) and (13), are applicable to uniform density ambients.
Analytic solutions will be given here for cases of single plumes in stag-
nant and flowing ambients with little or no stratification. For a single
plume in a stagnant, uniform ambient:
Sa = 0.074 Z5/3 Q~2/3	(14)
1 /4
For a single plume in a uniform ambient flowing at U > 0.1 (Q/Z) , dilution
can be calculated by:
Sa = 1,6 Z2 U4/3 Q_1	(15)
Equation (15) should be applied when:
dp/dz < 0.18 Z"3 if2 Q A0	(16)
The criteria presented for the application of equation (14) in Figure la is of
limited value because of the uncertainty in predicting dilutions in stagnant,
stratified ambients when the height of rise is limited and much less than the
trapping level height. In very siight stratification, this is the case.
However, in either the flowing or the stagnant case, results from equation
(14) or (15) should be compared with equation (3) or (6) when stratification
is slight. The lesser result is correct since stratification should reduce
the dilution of single plumes.
23

-------
PLUMES
Stagnant
F1owi ng
U < 0.0036 (A0Q)1/4 (dp/dz)1/8 < U
Stratified Sa =
0.46 Q"1/4 A03/4 (dp/dz)"5/8
Sa =
3.0 (U/Q)1/3(A0/dp/dz)2/3
Z < 4.6 (A0Q)1/4 (dp/dz)"3/8 Z < 3.3 (A0Q/U(dp/dz))1/3
U < 0.036 (A0Q)1/4 (dp/dz)1/8 < U
Surfacing
Sa =


Sa =


0.10 Z (A0/$)1/2 (dp/dz)"1/4

0.92 Z (U/Q)2/3 (A0/dp/dz)1/3

dp/dz < 3.
3 z"8/3Z2/3 Ao2
dp/dz < 0.18 Z"3 I
l"2 A0Q


U < 0.1 (Q/Z)1/4 < U





1
Uniform
Sa = 0.074
z5/3 q-2/3

Sa = 1.6 Z2 Q"1 U4/3
Figure la. Application criteria summary for the simplified analysis of plumes
in various ambients.
24

-------
MERGING PLUMES
Stagnant
Flowing
1
U < 0.014 (A0q) < U
Stratified
Surfacing/
Uniform
1/3 (dp/dz)"1/2 Z < 2.6 (A0q/U(dp/dz))1/2
1/3
U < 0.054 (A0q)
Sa =
2.7 (A0U/q(dp/dz))1/2
Sa = 1.0 Z (U/q)
-2/3 A 1/3
Sa =
-1/3 A 2/3
(dp/dz)
-1/2
Figure lb. Application criteria summary for the simplified analysis of merging
plumes in various ambients.
25

-------
SECTION 4
"PLUME" MODEL
THEORETICAL DEVELOPMENT
The computer model PLUME (5, 19) considers a buoyant plume issuing at an
arbitrary angle into a stagnant, stratified environment. Two zones of plume
behavior are considered. Profiles of velocity and concentration are considered
to be uniform or "top hat" at the port orifice. The region required to develop
fully established profiles is the zone of flow establishment. Beyond this
point, similarity profiles are assumed for velocity and concentration. This
allows the equations for mass continuity, momentum, and density disparity to be
integrated across the plume, reducing them to one independent variable, the
length along the axis of the plume.
Zone of Flow Establishment
The length of flow establishment (sg) has been found to depend on the
initial Froude number (12). As initial Froude numbers approach infinity,
s = 5.6 D. At low Froude numbers, this distance can be much shorter. The
e
centerline value of concentration is assumed equal to the initial discharge
conditions. The inclination of the plume's axis at sg is determined by linear
interpolation using 0, sg, and Fn. The centerline velocity at sg depends on
the Froude number and, at low Froude numbers, may be greater than the initial
velocity. To avoid estimation, the velocity at s (V ) is determined by s^
e me	e
and empirical spreading coefficients.

-------
Zone of Established Flow
The following assumptions are made:
a)	flow is steady and incompressible,
b)	the plume is fully turbulent,
c)	turbulent diffusion is significant only in the radial direction,
d)	pressure is hydrostatic throughout,
e)	plume flow is axisymmetric.
It is also assumed that the distributions of velocity and concentration are
Gaussian across the plume, so that:
V/Vm = exp (-k (r/s)2)	(17)
C/C^ = exp (-k|j(r/s)2)	(18)
where V and C are the velocity and concentration, respectively, at some radial
distance, r, and the subscript m refers to the centerline. The distance along
the plume's trajectory is s, and k and n are the empirical coefficients.
The integrated governing equations are:
conservation of mass,
d/ds / V r dr = -r. u.	(19)
o	b b
conservation of momentum,
d/ds f V2 r dr = sine / (p0o_p)g r/pd dr	(20)
conservation of density disparity,
d/ds f V A r dr = /" V r dr  dp^/ds	(21)
O	O
(Po-Pd)
conservation of pollutants,
d/ds f V C r dr = 0	(22)
O
Substituting in the velocity profiles (17) and (18), and simplifying
yields:
27

-------
momentum,
d/ds (Vm3 s3/k3/2) = 3 g((p0-pd) Vm Am s pk3/2 sine	(23)
density disparity,
d/ds (Vm Am s2/k(p + 1)) = (l/(Po-Pd)) Vm s2/k-dPo6/ds	(24)
conservation of pollutant,
Vm Cm s2/k(|J + ]) = 2 Vo Co/4	(25)
The equations can be cast dimensionlessly as:
momentum,
dE*/ds* = 3/Fn2-((1 + M)/pk1/2)s* R* sine	(26)
density disparity,
dR*/ds* = sVk1/2  E*1/3 (dp/ds*)	(27)
conservation of pollutant,
Cm = (k1/2 (M + 1 ))/(4E"tl/3 s*)	(28)
where s* = s/D, V* = Vm/V0, E* = (VJ s*/k1/2)3, R* = V* Am s*2/k(|j + 1) =
e*1/3	sVkl/2 (m + 1)( and	= Poo/(Po-pd)
The initial conditions for the zone of established flow are obtained by evalu-
ating E* and R* at s* = sg*. It is assumed that ambient stratification is
negligible in the zone of flow established so A = C and R * = 1/4. Since
33	m m	e
1/2	3
Cme = C0, Eg* = (k Cm + l)/4 se*) The values of p and k are mean values
of the limiting cases of purely-buoyant and purely-momentuni jets found by
Abraham (12).
The angle of the plume (6) can be evaluated by considering conservation
of horizontal momentum,
d/ds f V2 cos0 r dr = 0	(29)
O
which combined with equation (17) yields
6 = cos'1 ((Ee*/E*)2/3 cos0e)	(30)
28

-------
MODEL DESCRIPTION
The zone of flow established is first solved depending on the Froude
number:
s* = 2.8 Fn2/3	Fn < 2
e
s* = 0.113 Fn2 +4.	2 < Fn < 3.1
e
s* = (5.6 Fn2)/(Fn4 + 18)1/2 Fn > 3.1
Then other initial conditions, 9e, Eg*, and Re*, are solved for. Equations
(26) and (27) are then solved at steps along the plume's axis by Runge-Kutta
approximation. E* is used to calculate the inclination 6 and concentration C.
Then the position of the plume's axis can be incremented by dx/ds = cos0 and
dz/ds = sin0, where x and z are horizontal and vertical coordinates, respec-
tively.
The centerline concentration is inverted to describe the centerline dilu-
tion. To convert the centerline dilution to a flux averaged dilution, the
distribution of concentration must be weighted by distribution of velocity.
With the assumed distributions, the flux average dilution can be found by:
Sa = 1.77 (1/CJ or = 1.77 Sm	(31)
m	m
where S is the centerline dilution from PLUME. Calculation of dilution is
m
terminated when the density of the plume equals the density of the ambient.
This is normally somewhat below the maximum height of rise but is where simi-
larity ends and above which the plume will spread and become passive, possibly
interfering with further dilution.
Computation of time and plume diameter have been added to the original
program. Those changes are set off by asterisks in the program listing given
here. The centerline velocity, averaged over each step length, is divided
into the step length to obtain time. The plume diameter is found by D = 0.308
s.
29

-------
EXAMPLE INPUT, OUTPUT, AND MODEL LISTING
Input
The example problem is that of a single port with a diameter of 0.25 m
discharging 0.1 m3/sec of effluent with a density of 1000 kg/m3. The dis-
charge is horizontal. The discharge site varies in density from 1020 kg/m3 at
the surface to 1025 kg/m3 at the discharge depth of 30 m.
A schematic of the computer-punch input cards for this problem follows.
30

-------
I 1
t
1
5
Header Information



)
8
1
Plume Test. Case



r
s
tr
1
L-.,
~ 1 1
A
; T i i ii i ; 1 r i ( i i i i
-n	!rr-
TTTT"f
r riTTTit  11 i		
i
ii " f i i i i 1 i :  ! i i
' I
tf
I 1
^rii?rL
ir
_L_!
JSj
 4

I
P o P G 
I 8 ! f g
q CJ pn 3
n
f S g s
 5 c
3* I ?
 S.2
 3,8
 0
tfl 3
5- 5 3
(9
X
CU
3 O
O O) T3
--s r
Q Q.C
3
T3

-------
CO
ro
DEPTH
1 13 F5.0
F10.0
F10.0
Fin.n

i 1'
METERS
ANGLE
DEPTH
Units, 1=MKS, blank^FPS
Number of ambient density points to be specified
Port angle measured from the horizontal
Port diameter, m.
Depth of port
Density of effluent, gm/cm
Blank
Total volume rate of discharge
Number of ports
Print out interval
filMIT WTTIT1
hh ro m

-------
H..HL 
DP(r)
o.o
30.0
RHO(I) |
1.0200 |
1.0250
F10.0
F10.0
Im
lC>!
SAL
Card 3
F10.0
Mi
Card 4
rm i
icj
in
nnnmi
1 !
DP (I) - Depth measured from the surface
RH0(I)- Density (gm/cm3), or temperature (C) if Sal specified
SAL - If density not specified in RH0(I) then salinity, o/oo.

-------
Output
The first part of the model output repeats the input and provides other
computed initial conditions. The Froude number given is the square of the
definition used here. The discharge velocity and characteristics of the zone
of flow establishment are also given in this section.
The next section of the output traces the plume and its characteristics
along its trajectory. T is a time based on the centerline velocity. S is the
distance along the centerline. X and Z are distances measured horizontally in
the direction of the discharge and vertically from the water's surface to the
center of the plume, respectively. Elevation is the vertical distance to the
centerline of the plume from the level of the discharge. D is the diameter of
the plume and DILN is the centerline dilution. The dilution and height of
rise to the trapping level are also given. As mentioned, equation (31) is
used to convert this to a flux averaged dilution. In this case, Sa = 1.77 Sm
= 27.3. The example output and program listing follow.
34

-------
S3 PLUME VERSION 2.3 9/12/77
GJ
cn
BUOYANT plume in a density stratified mediae
PLUME TEST CASE
CASE NO. 1 INITIAL CONDITIONS ....
UNITS: mks
PORT ANGLF ...... 	 . .	0
FROUOE NUMRER 		67.7
LENGTH FOR FLOW ESTABLISHMENf . . . 1.40
INTEGRATION STEP LENGTH 		.076
PRINTOUT INTERVAL . 		3.00
XO		 1.39
ZO 	 ..........	29.84
DISCHARGE DENSITY . 	 ... 1.00000
PORT DEPTH		30.00
FLOWRATE	 1.OOOOOE-O1
NUMBER OF PORTS 		 1
DISCHARGE VELOCITY 		2.04
PORT DIAMETER 		2.50000E-01
DENSITY STRATIFICATION DEPTH RHO
0 1.02000
30.00 1.02500
T
S
X
Z
D
ELEV
THETA
2.35
4.29
4.06
28.85
1.32
1.15
36.2
6.43
7.25
5.99
26.62
2.23
3.38
59.1
11.51
10.22
7.27
23.95
3. 15
6.05
60.1
17.64
13.18
8.30
21.17
4.06
8.63
70.3
25.61
16. 14
9.41
18.42
4.97
11.58
63.1
33.91
18.12
10.74
17.07
5.58
12.93
0
LAST LINE ABOVE IS FOR MAXIMUM HEIGHT OF RISE.
TRAPPING LEVEL IS	21.5 WITH DILUTION OF 15.4321
HEIGHT OF RISE= 28.3 PERCENT OF DEPTH
OILN
3.3354
7.0703
11.6926

-------
0001	PROGRAM PLUME
0002	COMMON GFK.FM,COSTH.SINTH.COSTHE.DS.Cl.C2E13.FLAG.GRAV
0003	DIMENSION ZD(SO).00(50).RHO(50).DP(50).CMT(5)
0004	INTEGER NDC,TRAPPD,FLAG.CHGOEN.METERS
0005	WRITEI6.110)
0006	110 FORMAT("0S3 PLUME VERSION 2.3 9/12/77")
0007	NOCASE=0
0008	CUBRT = 1,/3.
0009	TW0THD=2./3.
0010	20 READC5.222)(CMT(I),1=1.5)
0011	222 FORMAT(5A8)
0012	IF (EOF(5)) CALL EXIT
0013	READ(5.73) NDC.METERS,NPTS.ANGLEtDIA.DEPTH,RHOJ.RFD.Q.FN.PS
0014	IF(EOF(5)) GOTO 7*
0015	73 FORMAT(2I1I3F5.07F10.0)
0016	FD=A9S(RFD)
0017	DO 102 1=1.NPTS
0018	RE AD(5.75) OP(I).RHO(I).SAL
0019	IF(EOF(5)) GOTO 999
0020	IFISAL.EQ.O.) GOTO 102
0021	RHO(I)=1.*.0 01SIGMAT(SAL.RHO(I)>
0022	102 CONTINUE
0023	CALL ASORT(OP,RHO,.NPTS)
10024	999 IF(NPTS.NE.1) GOTO 76
10025	NPTS=2
10026	DP(2)=DEPTH
0027	RH0(2>=RH0<1)
0028	76 N0CASE=N0CASE*1
0029	75 FORMAT(8F10.0)
0030	6RAV=32.172
0031	IF(METERS) GRAV=9.80665
0032	DO 55 1 = 1.NPTS
0033	IF(DP(I).GE.DEPTH) GOTO 56
0034	55 CONTINUE
0035	WRITE(6.59) NOCASE
0036	59 FORMATC'-NO DENSITY INFORMATION FOR JET LEVEL. EXECUTION FOR".
0037	  CASE NO.".12." DELETED.")
0038	T0 20
0039	56 NP=I
0040	NM=I-1
0041	RHOB =(OEPTH-DP(NM))(RHO(NP)-RhO(NM))/(DP(NP)-DP(NM))*RHO(NM)
0042	niSP=RHOJ-RHOH
0043	00 54 I=1.NM
0044	J=NP-I
0045	ZD(I)=(OEPTH-OP(J))/DIA
 0046	54 r>G(I) = (RHO(J.l)-RHO(J> ) D I A/(DI SP (OP (J* 1) -DP (J)) )
:0047	IF(NDC) GOTO 58
'0048	U0=Q/(FN.78539620IA*DIA)
>0049	RFD=UOUORHOJ/<-DISPDIAGHAV)
0050	FO=A0S(RFD)
10051	58 .IF (FD.LE.4.01.0R.FD.GE.9.99) GOTO 61
'0052	S=.113FU*4.
>0053	GOTO 62 j *-o
10054	61 IF (FD.LE.4.01) S = 2.8FDCUBRT	/
10055	lF if/fO.Q^ <1.3 ^ S = 5". 0/SV(?.f IF0
10056	62	TEMP=ATANU.416667S/FD)
) 00 5 7	THETA0=.Oil1ANGLE(1.5708-TEMP).TEMP
10058	C0STHE=C0S(THETA0)
10059	SINTHE=SIN(THETAO)
10060	DSI=DEPTH/(177.DIA)
10061	IF(OGU) .EO.O.) GOTO 77
10062	DGTEMP=.01/DG(1)
J0063	IF(DGTEMP.LE.O.) DGTEMP=-DGTEMP
36

-------
00064'	"hS I =. 12 1.*  f#2^ ( ACtT5lrb'*l DGTEMPj )
00065	77 NPO=IFIX(PS/(DSIDIA)~.5)
00066	N=0
00067	Z=SSINTHE
00068	X=SCOSTHE
00069	F=<4./S)<3
00070	R=.25
00071	C1=E**TW0THD
00072	C2=75/RFD
00073	IPTS=1
00074	6=PG1IPTS)
00075	7LIM = ZDUPTS>
00076	CHGDEN=0
00077	FLAG=0
00078	TRAPPD=0
00079	XP-X*DI A
00080	7P=DEPTH-ZDIA
00031	nSIP=DSIDIA
00062	SP=SDIA
00083	C
00084	T=2.SP/(U0tl.336.67(SS)))
00085	C ooooooo(H>ooooooooooooooooo4ooeoooooooooeoi>
00086	C
00037	C WHITE INITIAL CONDITIONS
00088	C
00089	WRITE(6.23) (CMT(I>* I = 1,5).NOCASE
00090	23 FORMAT C  0A BUOYANT PLUME IN A DENSITY STRATIFIED "
00091	MEDi A".I = 1.NPTS)
00117	24 FORMAT(DENSITY STRATIFICATION DEPTH RHO"//
00118	10<24X.F6.2,F11.5/>)
00119	C
00120	C SET INITIAL CONDITIONS
00121	C
00122	DS=DSI
00123	WRITE(6.109)
00124	C
00125	109 FORMAT < lH-,6X"T".9X,S".9X.,,X',9XZ"t9X.D"7X"ELEV'S
00126	6X."THETA"f4X."DILN")
00127	C 									
00128	GOTO 16
00129	11 DS=DSI
37

-------
urso
45
DElX^COS fHDS
0131

PELZ=SINThOS
'0132

DELE=FK
'0133

OELT=FM
0134

SP0S2=S*DS*.5
0135

no 10 1=1,2
0136

CALL SDEWIV(SPD52,E.5FK,R.5*FM)
0137

0ELX = DELX*2.C0,iTHDS
>0138

0ELZsUELZ*2.#SINTHDS
'01 39

0ELE=DELE*2.#FK
10140

0ELT=DELT*2.FM
>0141
10
CONTINUE
'0142

CALL SOERIV(S*OSE*FkR*FM)
'0143

ZLAST=Z
>0144

7INCR=(DELZ*SINTHDS)/6.
>0145

Z=Z*ZINCW
>0146

IF(CHGOEN) GOTO 41
0147

IF(Z.GT.ZLIM) GOTO 40
<0148
4 3
X = X*(DELX*COSTHOS)/6.
;0149
C
  o o    o  o o e   o   o   o      ft o tt
>0150

E0LD=E
>0151
C
fteftftftoft
>01 52

E=E*(DELE*FK)/6.
J0153

R=R*(0ELTFM)/6.
10154
C

10155

0T=.218DIA0S/(U0(E0LD,l(l./3.)/S*E0159

IFIE.LE.O.) FL AG=1
>0160
c

10161
c
THIS STOPPING CRITERIA IS BASED ON VELOCITY GOING TO ZERO
>0162
c

>0163

IF(FLAG) GOTO 13
>0164
16
CALL SDERIV(S,ER)
>0165

IF(TRAPPD) GOTO 14
>0166

IFIR.GT.O.) GOTO 15
10167

RAT =R/(RO-R)
>0168

F13TRP=E13*(E13-E130)
>0169

7TRAP=DEPTH-(Z*ZINCRRAT)OIA
>0170

SMTRAP=.245(S*DSRAT)E13TRP
'0171

TRAPPD=1
>0172

GOTO 14
>0173
15
RO=R
>0174

F130=E 1 3
>0175
14
N=N* 1
>0176

IF(N^(N/NPO)*NPO.NE.O) GOTO 11
>0177
13
XP=XD1A
>0178

ELEV=ZOIA .
>0179

?P=OEPTH-ELEV
>0160

SP-S01 A
>0181
C

>0162

R=.308SP
>0183
C
o u  o  ft  ft # ft ft o  a     o  ft  a  a o  ft ft ft  ft o a ft ft ft ft ft ft ft ft  o ft  ft ft ft
10184

nlLN=.245S0185

THETA=ACOSF(COSTH)57.2958
>0186

IF(.NOT,TRAPP01 GOTO 7?
10187
C
ftaftaofta0aaaafte
JO 188

WRITE(6,101) TSPXPZPB,ELEV,THETA
10189
C
aeftftftftattftftefteaafteafteo
10190

GOTO 71
10191
C

JO 192
72
WRITE(6,101) T,SP,XP,ZP,8,ELEV,THETA,0ILN
>0193
101
FORMAT(6F10.2,F10.1,F10.4)
10194
C
- e  a a e    a   a                  a      o    a a       
50195
71
IF I .NOT.FLAG) GOTO 11
38

-------
J0196
30197
00198
30199
30200
30201
30202
J 020 3
30204
30205
30206
30207
30208
30209
30210
30211
90212
90213
90214
0O21S
00216
00217
00218
00219
00220
00221
00222
00223.
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00236
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
PlJfcf'l M=100(2TRaP/0ePTm 100.)
WRITE(6.65)
65 FORMAT("-LAST LINE ABOVE IS FOR MAXIMUM HEIGHT OF RISE.")
67	FORMAT{'STRAPPING LEVEL IS'SFlO.lt
 ' WITH DILUTION OF".F10.4)
TF(TRAPPO) WRITEI6.67) ZTRAP.SMTRAP
IF(.NOT.TRAPPO) WRITE(6.6B)
IF(TRAPPD)WRITE(669) PI1EPTH
69 FOHMAT (" HEIGHT OF RISE = ".F5. 1." PERCENT OF DEPTH")
68	FORMAT("ATRAPPING LEVEL NOT REACHED")
GOTO 20
C
C FIND NEXT STRATIFICATION AND RECOMPUTE LAST STEP IF
C NECESSARr
C
40	OS=DS(ZLlrt-ZLAST)/(Z-ZLAST)
CALL SDERmS.E.R)
CHG0EN=1
Z=ZLAST
GOTO 45
41	CHGDEN=0
IPTSIPTS*1
IF(IPTS.GT.NM) GOTO 42
G=DG(IPTS)
ZLIM=ZD(IPTS)
GOTO 43
42	WRITE(644)
44 FORMAT("OPLUME HITS SURFACE"/
 HEI6HT OF RISE3100.0 PERCENT OF DEPTH")
flag=i
ROTO 43
74 CALL EXIT
END
SUBROUTINE SDERIV(StE.R>
COMMON G.FK FM COSTH.SINTHtCOSTHEDSC1C2E13FLAGGRAV
INTEGER FLAG
F13=0.
TF(E.LE.O.) GOTO 3
E13=E" < 1./3.)
COSTH=COSTHE*C1/(E13E13)
IF(C05TH.LT.1) GOTO 1
C
C THIS STOPPING CRITERIA IS BASED ON PLUME BECOMING HORIZONTAL
C AGAIN
C
FLAG=1
SINTH=0.
COSTH* 1 .
GOTO 3
1	SINTH=5QRT(1.-C0STHC0STH)
3 FK=C2SRSINTHDS
FM=.109SE13*GDS
RETURN
END
FUNCTION S1GMAT(SAL.T)
SIGO=<((6.8E-6#SAL)-4.B2E-4)SAL*.8149)SAL".093
R=1.E-6T(<.01667T-.fll64)T*18.03)
A*.001T*((.0010843*T-.09818)T*4.7867)
SUMT=(T-3.9fl)(T-3.9fl)o(T*283.)/(503.57*(T*67.26))
SIGMAT=(SIGO*.l324){l.-A+B'MSIGb-.1324))-SUMT
RETURN
END
SUBROUTINE ASORT (DPRMO 'Kj PTS)
DIMENSION DP(1).RHO(l)
NESTED=NPTS
L=NESTED-1
39

-------
00262
00263
00264
00265
00266
00267
00268
00269
00270
-0021L
0027 S
00276
00277
no,10 M=1L
nested=nested-i
no 10 I=l.NESTE0
IF(OP 11.LE.DP 
DP 11 * 11=DUMHY
nUHMY=KHO(I)
RHO( Ii =RHO(I ~ 1)
JPHO(1*1)=DUMMY
GOTO 10
""CONTINUE"
RETURN
END
40

-------
SECTION 5
"OUTPLM" MODEL
THEORETICAL DEVELOPMENT AND MODEL DESCRIPTION
The computer model OUTPLM (18, 23) considers a plume element or puff. By
following this element, the characteristics of a continuous single plume in a
flowing ambient are described. The original model has been adapted for ocean
discharges. Density, velocity, temperature, and salinity are assumed to be
average properties of the element. The quantities of mass, momentum, and
energy are conserved. Two explicit functions for entrainment are used. An
equation of state is used to calculate the density of the ambient and the
plume element at each time step.
Initially, the plume segment has a horizontal velocity component (compu-
ter variable U) and a vertical component (V). It has temperature (T) and
salinity (S) and is injected into a uniformly flowing (UW) ambient. The
segment has a radius (B) and thickness (H) and is tracked in two dimensions,
depth (Z) and horizontal distance (X).
Entrainment brings ambient mass (plus momentum, temperature, and salin-
ity) into the plume element. Entrainment is assumed to consist of two mechan-
isms. One mechanism (EINS) is due to the impingement of current on the plume.
That mass which crosses the exposed area of the plume element projected by
currents is entrained. Since the plume element is assumed to be cylindrical,
two corrections to the projected area are included in EINS. They are a radius
correction and a correction for the curvature of the plume trajectory (see
(23)). The other mechanism (ZWEI) is aspiration entrainment which captures

-------
0.1 times the product of the external area of the plume element and its shear
velocity. The assumption that the total entrainment is the larger of these
two mechanisms gives good agreement with data.
In the computer program, the entrained mass (DM) is added to the ele-
ment's mass (PM) to become the new mass (SUM). The new temperature and salin-
ity of the element are the averages of the old values and the entrained ambi-
ent values (TA and SA) weighted by their relative masses. The new density
(DEN), and thus buoyancy, creates a vertical acceleration on the plume seg-
ment. Since the element is considered to be one of a train, each following
the preceding element, drag is assumed to be negligible and the vertical
velocity is incremented (by DV). The horizontal velocity is found by conser-
vation of horizontal momentum with the entrainment (DM) having the velocity of
the current. The segment length (H) is changed in proportion to the total
velocity (VEL) to conserve mass and pollutant. The radius (B) is changed to
correspond to the new mass and density. Dilution (DIL) is calculated by
comparing the initial volume to that of the element. The program terminates
execution when the vertical velocity reaches zero, the surface is reached, or
length scales or execution step limits are reached.
EXAMPLE INPUT, OUTPUT, AND MODEL LISTING
Input
Consider, for example, a horizontal, single port discharge into a flow-
ing, stratified ambient. The maximum discharge velocity is 2.0 m/sec. The
discharge port radius is 0.125 m. The temperature and salinity of the waste
are 15.0C and 1.09 /oo, respectively. The ambient at the surface is 15.0C
and 27.2 /oo. At 30 m, the level of the discharge, the ambient is 15.0C and
42

-------
33.71 /oo. Currents are assumed to be 0.10 m/sec. The input punch cards for
this case follow.
43

-------
jul


uw
T
15.00;
F10.3
F10.3 i
Si
1 .09!
F10.3
e :
125
F10.3 i
1.0
F10.3
A ;
o.io
F10.3
125
F10.3
minmmi
riTTTTTTTTT
V - Vertical velocity of the plume, m/sec.
UW- Current speed, m/sec.
T - Temperature of the plume, C.
S - Salinity of the plume, o/oo.
H - Thickness of the plume element C^B), m.
E - Impingement entrai.nment coefficient, (recommended =
A - Aspiration entrainment coefficient (recommended = 0
B - Plume radius, m.
1.0)
.1)
5	moo
p	x co cr
I	n> -s I
*	3 CL -a
i	"a r
1
13
"a
*
beamu
J.

-------
u
KSCALE
YSCALE
NPTS
HERB
] |
tR| |
2.0
500
50C
2
j
1000 i 25
1
! F10.3 1 F5.0
F5.C
15
15 i 15
\
\
!
{
i




*
1
?
vi
i! j
TTCT
1
i _

I
i

1
!
.-ULillL
| |

Mill
1 | i 1 ! 1 M ; i! ;
1 1 : i 1 i ! I i i ! f
Hi
-i I ;
j ! !	
Ml

O
. i  rj;=
1 !oi
H
J . i

i is
1 j i - ^
	 j...	 !..	i j 	;c. 4
cn
U - Horizontal Velocity of the plume.
XSCALE - Number of diameters in the X direction at which run will be cut-off.
YSCALE - Number of diameters in the Z direction at which run will be cut-off.
NPTS - Number of points at which ambient is specified.
ITERB - Number of steps allowed.
IR - Interval of printout.
mOO
x o> cr
o> I
3 Q- -o
r~
-T03
CD
3
-o
C
c+

-------
DP (I)
SACI) :
TA r :
27.20 <
15.00
Card 3
33.71 i
5.00
Card 4
F10.3
F10.3
F10.3 i
DP(I) - Depth where ambient is specified, m.
SA(I) - Salinity, o/od.
TA(I) - Temperature, C.
m
C
O
X
a>
c=
0J
-s
-H
3
a.
"U
TD

i

OJ
3
a>


ti
a>

3
r+

XJ
O

c


c+



-------
Output
Model output and program listing follow. The line above the header
echoes the first two input cards. Ambient stratification is then given with
input depth, salinity, temperature, and computed sigma-t densities (density,
kg/m3, 1000). Other discharge conditions are given. K is a velocity
ratio, = U/UW. Q is the volume rate of discharge. The first data line in the
listing are the initial conditions. The horizontal distance in the direction
of the jet discharge is the X dimension. Z is depth measured from the water's
surface. B is the radius, THICK is the plume element thickness along the
centerline. MASS is the mass of the plume element being considered. EINS and
ZWEI are the previously mentioned entrainments. DENDIF is the plume's density
difference with the ambient. DILUTION is an average volumetric dilution, the
inverse of CONC. TEMPDIF is the temperature difference between the plume
element and the ambient at that level.
47

-------
0	.100 15.000	1.090	.125	1.000	.100	.125	2.000 500 500 2 1000 25
      
OUTFALL BUOYANT (JET) PLUME IN FLOWING. STRATIFIED AMBIENT
0S3 VERSION 2.3 5-16-79
AMBIENT STRATIFICATION- DEPTH,M SALIN TEMP.C SIGMaT
0 27.20 15.00 20.01
30.00 33.71 15.00 25.00
K FROUOE	(i CURRENT
2.00E 01 8.O0E 00	9.82E-02 1.00E-01
	MODEL INPUT (LINE 1) AND MODEL OUTPUT	
x Z B THICK MASS EINS ZWEI DILUTION DENDIFF HOR VEL	VER VEL TOT VEL TEMPO IF
OE 00 3.00E 01	1.25E-01 1.25E-01 6.14E 00 1.12E-02 1.91E Ul 1.00E 00 2.50E 01 2.00E 00	OE 00 2.00E 00 OE 00
4.44E-03 3.00E 01	1.25E-01 1.24E-01 6.18E 00 1.12E-02 1.91E 01 1.01E 00 2.48E 01 1.99E 00	5.43E-04 1.99E 00-2.33E-10
1.21E-01 3.00E 01	1.47E-01 1.06E-01 7.30E 00 2.97E-03 4.97E-02 1.18E 00 2.10E 01 1.70E 00	1.36E-02 1.70E 00 2.33E-10
2.67E-01 3.00E 01	1.74E-01 9.02E-02 8.68E 00 4.27E-03 5.91E-02 1.40E 00 1.77E 01 1.44E 00	2.76E-02 1.44E 00 OE 00
4.41E-01 3.00E 01	2.05E-01 7.69E-02 1.03E 01 6.21E-03 7.03E-02 1.66E 00 1.49E 01 1.23E 00	4.22E-02 1.23E 00 1.40E-09
CO 6.50E-01 3.00E 01	2.42E-01 6.57E-02 1.23E 01 9.13E-03 8.36E-02 1.98E 00 1.25E 01 1.05E 00	5.7SE-02 1.05E 00 9.31E-10
8.99E-01 3.00E 01	2.84E-01 5.64E-02 1.46E 01 1.36E-02 9.95E-02 2.34E 00 1.05E 01 8.99E-01	7.49E-02 9.02E-01 2.33E-10
1.20E 00 2.99E 01	3.34E-01 4.86E-02 1.74E 01 2.05E-02 1.18E-01 2.78E 00 8.83E 00 7.72E-01	9.38E-02 7.77E-01 2.33E-10
1.55E 00 2.99E 01	3.90E-01 4.22E-02 2.06E 01 3.13E-02 1.41E-01 3.31E 00 7.42E 00 6.65E-01	1.15E-01 6.75E-01 6.98E-10
1.97E 00 2.98E 01	4.54E-01 3.70E-02 2.45E 01 4.78E-02 1.68E-01 3.93E 00 6.23E 00 5.75E-01	1.37E-01 5.91E-01 9.31E-10
2.47E 00 2.97E 01	5.26E-01 3.28E-02 2.92E 01 7.25E-02 1.99E-01 4.66E 00 5.22E 00 4.99E-01	1.62E-01 5.25E-01 6.98E-10
3.03E 00 2.94E 01	6.03E-01 2.96E-02 3.47E 01 1.08E-01 2.37E-01 5.54E 00 4.36E 00 4.36E-01	1.87E-01 4.74E-01 9.31E-10
3.65E 00 2.91E 01	6.86E-01 2.73E-02 4.13E 01 1.56E-01 2.83E-01 6.59E 00 3.62E 00 3.82E-01	2.11E-01 4.37E-01 1.16E-09
4.32E 00 2.87E 01	7.72E-01 2.56E-02 4.91E 01 2.18E-01 3.37E-01 7.83E 00 2.98E 00 3.38E-01	2.31E-01 4.09E-01 1.40E-09
5.02E 00 2.82E 01	8.64E-01 2.43E-02 5.84E 01 2.96E-01 4.01E-01 9.31E 00 2.43E 00 3.00E-01	2.47E-01 3.88E-01 1.40E-09
5.76E 00 2.75E 01	9.64E-01 2.32E-02 6.94E 01 3.92E-01 4.77E-01 1.11E 01 1.94E 00 2.68E-01	2.57E-01 3.71E-01 4.66E-10
6.52E 00 2.68E 01	1.07E 00 2.22E-02 8.26E 01 5.11E-01 5.67E-01 1.32E 01 1.51E 00 2.41E-01	2.62E-01 3.56E-0X 6.98E-10
7.32E 00 2.58E 01	1.20E 00 2.12E-02 9.82E 01 6.62E-01 6.74E-01 1.56E 01 1.13E 00 2.19E-01	2.60E-01 3.40E-01 6.98E-10
8.15E 00 2.48E 01	1.34E 00 2.01E-02 1.17E 02 8.06E-01 7.42E-01 1 .86E 01 7.98E-01 2.00E-01	2.52E-01 3.22E-01 6.98E-10
8.97E 00 2.38E 01	1.52E 00 1.87E-02 1.39E 02 9.59E-01 7.75E-01 2.2IE 01 5.11E-01 1.84E-01	2.36E-01 2.99E-01 1.16E-09
9.79E 00 2.27E 01	1.73E 00 1.71E-02 1.65E 02 1.14E 00 7.84E-0) 2.63E 01 2.69E-01 1.71E-01	2.14E-01 2.74E-01 9.31E-10
1.06E 01 2.17E 01	1.99E 00 1.54E-02 1.96E 02 1.36E 00 7.57E-01 3.13E 01 7.07E-02 1.59E-01	1.87E-01 2.46E-01 2.33E-10
1.15E 01 2.08E 01	2.31E 00 1.36E-02 2.33E 02 1.61E 00 6.82E-01 3.72E 01-8.51E-02 1.50E-01	1.57E-01 2.17E-01 2.33E-10
1.23E 01 1.99E 01	2.70E 00 1.18E-02 2.78E 02 1.92E 00 5.61E-01 4.42E 01-1.97E-01 1.42E-01	1.24E-01 1.88E-01-b.98E-10
I.32E 01 1.93E 01	3.I6E 00 1 .02E-02 3.30E 02 2.28E 00 4.24E-01 5.26E 01-2.66E-01 1.35E-01	9.14E-02 1.63E-01-1.40E-09
1.41E 01 1.88E 01	3.69E 00 8.93E-03 3.93E 02 2.71E 00 3.19E-01 6.25E 01-3.00E-01 1.30E-01	6.01E-02 1.43E-01 4.66E-10
1.52E 01 1.84E 01	4.26E 00 7.98E-03 4.67E 02 3.23E 00 3.07E-01 7.44E 01-3.09E-U1 1.25E-01	2.64E-02 1.28E-01 9.31E-10
1.61E 01 1.83E 01	4.49E 00 7.71E-03 5.04E 02 3.48E 00 6.63E-01 7.97E 01-3.04E-01 1.23E-01-2.07E-03 1.23E-01 2.33E-10
NUMBER OF STEPS= 636

-------
00001
00002
00003
00004
00005
0000b
U0007
000UB
oooov
00010
00011
00012
00013
00014
00015
0001b
0001 7
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
0003S
00036
00037
00038
00034
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
0005H
00059
00060
00061
00 062
00063
PROGRAM OUTPLM
DI MENS I ON DP(30) ,SA(30) ,TA(30).DENP(30)
DATA (G=9.) , < I W0=2. I ,
6	FORMAT(BF10.3*/,IF]0,3*2F5.0,3151
9 FORMAT(3F10,3)
14 FORMAT 126X,4F6.2/)
READ (10,61 V,UW,T,SH,E,A,B,U,XSCALEZSCALE,NPTS,ITERB, IR
WRITE(11*51 V*UW)
8 CONTINUE
WRITE (11,14) ((OP(I)iSAtl)TA(I),DNP(1)> ,1=1,NPTS)
DT =OD=DB=1.
X=BSAVE=ZWEI	= ZERO
Z=DPINPTSI
OEL = T-T A(NPTS)
AK = 99<59.
DIL = 1 
BA=B
BO=BTWO
VI =VEL=SUR7(VV*UU)
OENA=1000.SIGMAT(SA(NPTS),TA(NPTS))
DEN=1000.*SIGMAT(S*T)
RHO=OEN
0EN01FF=DENA-DEN
PM=PMO=PIHRHOEN
(J=HM/DENYEL/H
IF(UW.NE.ZERO) AK=VEL/UW
FR=\/EL/SURT (DENl)lFF/OEN*TWOBG)
7	FORMAT!/"'	"/
1"	K FROUDE	0 CURRENT"/1X4E9.2//
2"	MODEL INPUT (LINF II ANO MODEL OUTPUT	X
3	Z	B	THICK	MASS	EINS	ZWEI DILUTION DENDIFF
4	HOrt VEL VEH VEL TOT VEL TEMPOIF")
WRITE(11,7) AK
ZWEI=TW0P1DENAHBABSIVEL-UW*U/VEL)A0T
DM=E1NS
IF(ZWEI.GT.EINS) DM=ZWEI
GAMMA=OM/.0069555/PM
DT=DT/GAMMA
EINS=EINS/GAMMA
49

-------
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
0007b
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
oooaa
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00103
00109
00110
00111
00112
00113
OUl 14
00115
dm=dm/gamha
13 FORMAT (1X.14E9.2)
IFtJ.EO.1)WHITE(11t13)X.Z.Pt.H,PM,EIN5tZWEI
1OIL tDENOIFF.U.VVELDEL
SUH=PM*DM
UO=U
U= /IH.NE.1)G0 TO 99
IF (RATIOZ.GT.ZSCALE.OR.RATIOX.GT.XSCALE)GO TO 999
98	RATIOX=X/BO
WATI0Z=Z/dO
DIL = (RHOPM)/
-------
SECTION 6
"DKHPLM" MODEL
THEORETICAL DEVELOPMENT
The computer model DKHPLM (7, 21, 24) is an approach to the problem of
merging plumes which considers the detailed dynamics of the process instead of
simplifying the problem to a idealized slot plume or to a combination of plume
and slot plume. DKHPLM considers three zones of plume behavior; zones of flow
establishment, established flow, and merging. The first two zones are based
on the analyses of Hirst (25, 26) for a plume in a stratified, flowing en-
vironment. In the zone of merging, neighboring plumes are superimposed. This
allows a smooth transition as single plumes begin to compete for dilution
water then gradually merge with their neighbors. Equations for the conserva-
tion of mass, pollutant, energy, and momentum are developed for all three
zones. Entrainment is an explicit function dependent on the local Froude
number, plume spacing, excess velocity, and ambient velocity. Similar lateral
profiles, a 3/2 power approximation of a Gaussian, are assumed for velocity,
concentration, and temperature. These profiles are superimposed in the merg-
ing zone. A complete theoretical development of this model is beyond the
scope of this paper but can be found in (21, 24, 25, and 26). The following
is a brief summary.
Zone of Flow Establishment
All quantities are assumed uniformly distributed in the plume at the
point of discharge. In the zone of flow establishment, these uniform profiles

-------
change to similar profiles as the boundary layer diffuses inward to the cen-
terline of the jet. The rate at which the profiles of velocity, concentration
and temperature develop may vary. The integrated forms of the governing
equations are:
conservation of mass,
d/ds / V r dr = -Tim r (r v) = E	(32)
o
conservation of energy,
d/ds J V (T-T ) r dr = - d T /ds J V r dr - lim r -  (r v'T') (33)
O .	00	O
conservation of pollutant,
d/ds f V (C-C ) r dr = -d C /ds J V r dr - lim r ->  (r v' C') (34)
O	00	00 O
conservation of momentum in the s equation,
d/ds f r dr = U E sinBj cos02 + J g (p - p)/pH) r dr sin02
o	o	00	Q
- lim r a (r u' v')	(35)
where 0X is the horizontal angle between the centerline and the x axis and 02
is the angle between the centerline and the horizontal. Two.additional inte-
gral equations have been developed from equation (35) to describe momentum in
two additional plume coordinants. These "natural" coordinants of the plume,
described in (25), are converted to conventional three-dimensional Cartesian
coordinates for model output. Implicit in the derivation of.these equations
are the assumptions that:
a)	flow is steady in the mean,
i		>
b)	flow is fully turbulent,
c)	fluid is incompressible and density variations are included only in
the buoyancy terms,
d)	all other fluid properties are constant,
e)	no frictional heating,
f)	pressure variations are purely hydrostatic,
52

-------
g)	changes in density are small enough to be approximated by a linear
equation of state,
h)	flow within the jet is axisymmetric,
i)	flow within the jet can be approximated as boundary layer flow,
j) the ambient is infinite in extent.
Several of the assumptions are compensated for in the solution. The zone of
flow establishment uses a special entrainment function (see equation (128) of
(24)) which is a function of local Froude number, velocities, port diameter,
spacing, and thickness of the developed flow region.
Zone of Established Flow
The equations of mass, energy, pollutant, and momentum (equations 32-35),
and the additional momentum equations mentioned above are also solved in the
zone of established flow, but are cast into a slightly different form. The
governing equations are written in a cylindrical coordinant system where  is
the circumferential angle around the plume and cross section and the independ-
ent variables are r and s. These are evaluated using the assumed 3/2 power
approximation to Gaussian lateral profiles. The angle is the angle between
the centerline projected to the xy plane and the x axis, 02 is the angle be-
tween the centerline and the xy plane. These angles relate the two coordinate
systems. Another entrainment function is used in the zone of established flow
which is a function of the local Froude number, velocities, plume diameter,
and spacing.
Zone of Merging
When adjacent plumes begin to overlap, the discharge is no longer consid-
ered axisymmetric. The distributions of plume properties are superimposed.

-------
Another entrainment function is used which also considers the variable en-
trapment surface during merging. A drag term is also introduced to account
for the additional bending of the plumes after merging.
MODEL DESCRIPTION
An equation of state for sea water has been added to the original program
(subroutine SIGMAT) replacing the linear equation of state. In the zone of
flow establishment, the system of six governing equations are solved simultan-
eously (subroutine SIMQ) and stepped forward in space by a Hamming's modified
predictor-corrector method (subroutine HPCG). This procedure continues until
velocity, temperature and concentration become fully developed. Subroutine
0UTP1 contains the results which are stored as initial conditions for the zone
of established flow.
In the zone of established flow, similar profiles and the integral method
allow solution of the six governing equations for the six unknowns initially
by Runge-Kutta integral approximation and then by the Hamming's modified pre-
dictor-corrector method. At the point where the plumes overlap, the assumed
similarity no longer applies. The merging plumes have axes of symmetry along
the discharge line and normal to it. Only one quadrant of a plume, taken to a
midpoint of the overlap area, is evaluated. The profiles are superimposed in
this region using integral similarity coefficients. Above the point where the
plume and the ambient have equal density, results are obtained by extrapola-
tion.
The program itself contains many comments and explanations which serve as
further documentation.
54

-------
EXAMPLE INPUT, OUTPUT, AND MODEL LISTING
Input
As an example case, suppose a diffuser 1250 m long discharges 12.5 m3/sec.
The ports are 0.178 m in diameter, oriented vertically, and spaced 5 m apart.
The ambient currents are assumed to be 0.05 m/sec and normal to the discharge
line. (Note that the angle of currents to the diffuser should be within about
20 of normal.) The velocity through the discharge ports is 2.0 m/sec. The
effluent is 15.01C and 1.09 /oo. (Note that the effluent should have some
finite temperature difference with the ambient.) The ambient at the level of
discharge is 15.0C and 33.71 /oo (1025 kg/m3). The ambient temperature
gradient is negligible but must be finite. The ambient salinity gradient is
-0.026 /oo per m (a density gradient of 0.020 kg/m3/m). Depth of the dis-
charge is 50 m. Computer punch card format for this input follows.
55

-------
HEADER INFORMATION (Nil)
0 DKHPLM EXAMPLE RUN
J3 - Number of runs to follow C0=one run).

-------
-.02600 !
15.01
1.09
15.00
33.71
.00001
F10.6
F10.6
F10.6 i
F10.6
F10.6
F10.6
F10.6
F10.6i
UO - Initial jet velocity, m/sec.
UA - Current speed, m/sec.
TO - Initial jet temperature, C.
CO - Initial jet salinity, o/oo.
TAO- Ambient temperature at the level of discharge, C,
CAO- Ambient salinity at the level of discahrge, C.
TIZ- Temperature gradient, C/m.
CIZ- Salinity gradient, o/oo/m.
m
 o
X
Q>
Q)
-s ^
3
Q. "O
-o


ro 2
(D

HH

3

"O

c

c+


-------
en
CO
r *- 
i
konecnr^
j
DO j
i
!'
TH1 |
]
i
TH2
i
1
SF !
SPACE
H
I
0.178
i
i
90.0 j
90.0 L
!
500.0 \
5.0
50.0
F10.6
1
l
F10.6 j
F10.6 !
i .
I
F10.6 j
F10.6
F10.6
5
SPACE
Port diameter, m.
Horizontal discharge angle (DEG) relative to the current	(70
-------
Output
Example output and program listing follow. The initial input conditions
are printed along with dimensionless conditions. The first section of data
listing concerns the zone of flow establishment. The columns list: length
along the plume axis (S), horizontal distance parallel with discharge line
(X), horizontal distance normal to the discharge line (Y), vertical distance
from the level discharge (Z), the horizontal angle from the plume's axis to
the discharge line (TH1), the angle the plume's axis makes from the horizontal
(TH2), the radius (B), potential core widths for velocity, temperature and
concentration (RU, RT, and RC), normalized centerline disparities of velocity,
temperature, and concentration with the ambient (DUCL, DTCL, and DCCL), the
ambient density normalized by the density of the discharge (PI), and non-
dimensional time (TIME). In the zone of established flow, width and average
dilution (Q/Qo) are also given. Time is given in seconds. The centerline
dilution before zone of merging is:
Sm = 0.52 Q/Qo	(36)
m
after complete merging (which never actually occurs):
Sm = 0.70 Q/Qo	(37)
based on the assumed distributions of concentration and velocity.
If stagnant conditions are specified, the dilution should be taken at the
point of buoyant equilibrium with the ambient. The location is announced in
the model output.
59

-------
SOLUTION TO MULTIPLE BUOYANT DISCHARGE PROBLEM WITH AMBIENT CURRENTS AND VERTICAL GRADIENTS
DKHPLM EXAMPLE CASE	I.D. = 0
DISCHARGE VELOCITY = 2.00-M/S o TEMP. = 1^,.01-DEG C	 SALINITY = 1.09-PPT
niAMETEW = . lfl-M  SPACING = 5.00-M  DEPTH =	50.00-M
AMBIENT CONDITIONS AT DISCHARGE ELEVATION. VELOCITY =	,05-M/S  TEMP. = 15.O0-DEG C  SALINITY = 33.71 PPT
AMBIENT STRATIFICATION GRADIENTS  TEMPERATURE	.000010-DEG C/M  SALINITY  -0.026000-PPT/M
AMHIENT CONDITIONS AT DISCHARGE - NONDIMENS IONAL
TFMP= .99933 SAL INITY = 30.92661 DENS IT Y= 1.02501 T EMPGR = .000000 SALGR = -0.002123 DENGR=-0.0000018
PORT SPACING L/D = 28.09
FROUDE NO =	9.58 VELOCITY RATIO = .025 STRATIFICATION NO = 14U96.6

ZONE OF
FLOW ESTABLISHMENT
ALL LENGTHS
ARE IN
METERS




S X
Y
Z
THl
TH2
B
RU
RT
RC
DUCL DTCL DCCL
PI
TIME
n o
0
n
90.00
90.00
0
.089
.089
.089
1.000 1.000 1.000
1.02501
0
.09 .00
.00
.09
90.00
99.91
.018
.083
.083
.083
1.000 1.000 1.000
1 .02501
1.0
in .oo
.00
. l 8
90.00
89.81
.038
.077
.077
.077
1.000 1.000 1.000
1.02501
2.0
.27 .00
.00
.27
90.00
89.70
.059
.070
.070
.070
1.000 1.000 1.000
1.02501
3.0
.3ft .00
.00
. 36
90.00
89.57
.082
.063
.062
.062
1.000 1.000 1.000
1.02500
4.0
.44 .00
.00
. 44
90. 00
89.43
.107
.054
.053
.053
1.000 1.000 1.000
1.02500
5.0
.53 .00
. 0 0
.53
90.00
89.28
.133
.045
.043
.043
1.000 .999 1.000
1.02500
6.0
.6? .00
.00
.62
90.00
89. 1 1
.162
.035
.033
.033
1.000 .999 1.000
1.02500
7.0
.7 1 -0.00
.01
. n
90.00
88.92
.193
.024
.021
.021
1.000 .999 .999
1.02500
8.0
.80 -0.00
.01
.80
90.00
88.72
.226
.011
.00 7
.007
.999 .999 .999
1.02499
9.0
STARTING LFNGTH,
T =
.844









.84 -0.00
.01
.84
90.00
88.62
.243
.004


.999 .999 .999
1.02499
.4
STARTING LFNGTH,
VELOCITY
=
.066








.87 -n.oo
.0 1
.87
90.00
88.57
.252 -
0.000


.999 .930 .930
1 .02499
9.7

ZONE OF
ESTABLISHED FLOW








S X
Y
Z
THl
TH2
WIDTH
DUCL
OTCL
DCCL
PI TIME
Q/UO

.87 -0.00
.01
.87
90.00
88.57
.51
.997
.930
.930
1.02499 .43
2.07

1.5ft .00
.05
1 .58
90.00
P4.59
1 .24
.429
.354
.354
1.02498 1.00
5.48

2.29 .00
. 14
2.2a
90.00
80.81
2.03
.284
. 196
.196
1.02497 2.03
9.99

3.00 .00
.2*
2.98
90.00
77.58
2.82
.223
.126
.127
1.02495 3.43
15.67

3.71 .00
.45
3.67
90.00
74.86.
3.59
.190
.088
.089
1.02494 5.12
22.48

4.43 .00
.65
4.36
90.00
72.50
4.35
.169
.065
 066
1.02492 7.04
30.37


ZONE OF MERGING
FLOW








5.14 .00
.87
5. 03
90.00
69.97
5.09
.154
.050
.051
1.02491 9.14
39.29

5.65 .00
1.14
5.69
90.00
66.86
5.71
.144
.040
.041
1.02490 11.40
47.81

6.56 .00
1.43
6.34
90.00
64 . 78
6.27
.138
.033
.034
1.02488 13.76
55.90

7.27 .00
1 .74
6.98
90.00
63.27
6.81
.132
.027
.028
1.02487 16.20
63.73

8.01 .00
2.08
7.63
90.00
62.08
7.35
.126
.022
.024
1.02486 18.81
71.60

9.43 .00
2. 77
8. 88
90.00
60.32
8.44
.113
.015
.017
1.02483 24.21
86.36

10.86 .00
3.49
10.11
90. 00
58.89
9.64
.099
.010
.011
1.02481 30.19
100.48


-------
DkmPLM EXAMPLE CASE
1.0. a
0
s
X
Y
2
TH1
TH2
WIDTH
DUCL
DTCL
DCCL
PI
TIME
Q/O
12.2*
.00
4.24
11.3?
90.00
57.55
11.18
.087
.006
.008
1.02479
36.87
114.06
13.75
.00
5.04
12.55
90.00
56.06
12.63
.084
.004
.005
1.02476
44.28
127.61
15.17
.00
5.85
13.72
90.00
54.44
14.13
.080
.002
.004
1.02474
51.69
140.44
16.60
.00
6.70
14.86
90.00
52.58
15.73
.076
-0.000
.002
1.02471
59.38
152.93
in.02
.00
7.58
15.98
90.00
50.38
17.43
.071
-0.002
.000
1.024ts9
67.41
165.00

PLUMF HAS
REACHED
ITS EQUILIBRIUM
HEIGHT -
STRATIFIED
ENVIRONMENT




ia.2
.00
8.10
16.59
90.00
48.94
18.45
.068
-0.003
-0.000
1.02468
72.11
171.58
19.00
.00
8.22
16.72
90.00
48.60
18.68
.067
-0.003
-0.001
1.02468
73. 17
173.02
19. lfl
.00
8.34
16.86
'90.00
48.25
18.91
.066
-0.003
-0.0 0 1
1.02467
74.24
174.45
19.36
.00
8.46
16.99
90.00
47.89
19.14
.065
-0.004
-0.001
1.02467
75.32
175.88
19.53
.00
8.58
17.12
90.00
47.52
19.38
.065
-0.004
-0.001
1.02467
76.40
177.29
19.71
.00
8.70
17.25
90.00
47.15
19.62
.064
-0.004
-0.001
1.02467
77.50
178.69
19.89
.00
8.82
17.38
90.00
46. 7b
19.86
.063
-0.004
-0.001
1.02466
78.60
180.08
20.07
.00
8.94
17.51
90.00
46 . 36
20.11
.062
-0.004
-0.002
1.02466
79.72
181.46
20.25
.00
9.06
17.64
90.00
45.95
20.35
. 062
-0.005
-0.002
1.02466
80.84
182.84
20.4 2
.00
9.19
17.77
90.00
45.53
20.60
.061
rO.005
-0.002
1.02466
61.97
184.20
20.6 0
.00
9.31
17.89
90.00
45.09
20.85
.060
-0.005
-0.002
1.02465
83.11
185.55
20.7B
.00
9.44
18.0?
90.00
44.65
21.10
.059
-0.005
-0.002
1.02465
84.26
186.89
20 . <36
.00
9.57
18.14
90.00
44. 19
21.36
.058
-0.005
-0.002
1.02465
85.42
188.21
21.14
.00
9.70
18.27
90.00
43.72
21.62
.058
-0.005
-0.003
1.02465
86.59
189.53
21.31
iOO
9.H 2
18.39
90.00
4 3.23
21.87
.057
-0.006
-0.003
1.02464
87.77
190.83
21 .49
.00
9.95
18.51
90.00
42.73
22.14
.056
-0.006
-0.003
1.02464
88.96
192.12
21 .67
.00
10.09
IB.63
90.00
42.21
22.40
.055
-0.006
-0.003
1.02464
90. 17
193.39
21.05
.00
10.22
18.75
90.00
41 .68
22.67
.054
-0.006
-0.003
1 .02464
91.38
194.65
22.03
.00
10.35
18.87
90.00
41.13
22.93
.053
0.006
-0.003
1.02463
92.61
195.90
22.20
.00
10.49
18.98
90.00
40.57
23.20
.052
-0.007
-0.003
1 .02463
93.85
197.13
22.47
.00
1 0.69
19.16
90.00
39.68
23.60
.051
-0.007
-0.004
1.02463
95.73
198.95
22.83
.00
10.97
19.38
90.00
38. 44
24.15
.049
-0.007
-0.004
1.02462
98.29
201.32
23. 1R
.00
11.25
19.60
90.00
37. 12
24.70
.04 7
-0.007
-0.004
1.02462
100.91
203.62
23.54
.00
11.53
19.81
90.00
35.70
25.27
.045
-0.008
-0.004
1.02462
103.58
205.85
23.89
.00
11.83
20.01
90.00
34. 19
25.84
. 044
-0.008
-0.005
1.02461
106.32
20S.00
24.25
.00
12.12
20.21
90.00
32.58
26.41
.042
-0.008
-0.005
1.02461
109.13
210.06
24,61
.00
12.43
20.40
90.00
30.86
26.98
. 040
-0.009
-0.005
1.02460
112.01
212.02
24.96
.00
12.73
20.57
90.00
29.02
27.55
. 038
-0.009
-0*006
1.02460
114.96
213.89
25.32
.00
13.05
20.74
90.00
27.05
28.12
.036
-0.009
-0.006
1.02460
118.00
215.66
25.67
.00
13.37
20.90
90.00
24.96
28.67
.034
-0.010
-0.006
1.02459
121.11
217.30
26.03
.00
13.69
21 .04
90.00
22.73
29.21
.032
-0.010
-0.006
1.02459
124.30
218.83
26.39
.00
14.03
21.17
90.00
20.36
29.72
.030
-0.010
-0.006
1.02459
127.58
220.22
26.74
.00
14.36
21.29
90.00
17.85
30.20
.029
-0.010
-0.007
1.02459
130.95
221.47
27. 10
.00
14.70
21.39
90.00
15.21
30.64
.027
-0.011
-0.007
1.02458
134.39
222.58
27.45
.00
15.05
21.48
90.00
12.45
31.04
.026
-0.011
-0.007
1.02458
137.91
223.52

DKHPLM example case





I.D.
= 0


S
X
r
Z
THl
TH2
HIDTH .
DUCL
DTCL
DCCL
PI
TIME
Q/O
27.81
.00
15.40
21.54
90.00
9.57
31 .36
. 025
-0.011
-0.007
1.02453
141.49
224.31
28.34
.00
15.93
21.61
90.00
5.09
31.71
.023
-0.011
-0.007
1 .02458
146.96
225.16
28.88
.00
16.46
21.64
90.00
.50
31 .88
. 023
-0.011
-0.007
1 .02458
152.51
225.63
29.06
.00
16.64
21.64
90.00
-1.03
31.91
.023
-0.011
-0.007
1.02458
154.36
225.72
PLOMF HAS REACHED ITS MAXIMUM HEIfiHT - STRATIFIED environment
MAXIMUM Pl.lJMF HEIGHT. EXTRAPOLATED	LINEAR = 44.52	PARACOLIC = 43.41
NO OF INTEGRATION STEPS= 898 NO OF HALVINGS = 11 ABSERR= .00100 PRMT<5>= 1.000 FINAL SPACE = 4.000

-------
00001
00002
00003
00004
00005
onoob
D0007
00008
00009
00010
000 11
0001?
00013
00014
II001S
000 If.
00017
000 1H
000 1*
00020
00021
00022
00023
00024
00026
00026
00027
(>0 0?ft
000 29
1)0030
90031
00032
00033
00034
00035
00036
J 0 0 3 7
>0038
00 0 39
0 004 0
10041
00042
.10043
1)0044
11004S
00046
00047
.') 0 0 4 R
;>0049
00 0 SO
00 051
DOO'52
0 0 0 b 3
i)OOS4
0 0 055
00056
00057
U0058
00059
J 0 0 6 0
)OOt>l
00062
00063
PROGRAM DKHPLM
THIS PROGRAM CALCULATES THE CHARACTFRISTICS OF A LINE UF EOUALLY
SPACFD BUOYANT DISCHARGES INTO A FLOWING STRATIFIED AMBIENT
IT IS ASSUMED THAT BUOYANCY IS DUE TO TEMPERATURE ANL) SALINITY
THE METHOD OF SOLN INVOLVES 7 ORDINARY DIFFERENTIAL EQUATIONS
WHICH ARE 1. CONSERVATION OF MASS! ?. CONSERVATIONS OF ENERGY?
3. CONSERVATION OF CONCENTRATION (SALINITY)! 4. DENSITY DEFICIENCY
(DEPTVFn FROM 2)3)5 5.6.7. MOMENTUM F'JNS IN THE Z (AXIAL).
K (VERTICAL) ) J(HORIZONTAL AND PARALLEL TO AMBIENT CURRENT)'01MECT1 ON
EON 4 IS USED ONLY IN ALGEBRAIC FORM.
The SIX EQUATIONS are WRITTEN FOR A CONTROL VOLUME WHICH IS FINITE
IN A DIRECTION PERPENDICULAR TO THE JET AND INFINITESIMAL IN THE
DIRECTION OF THE JET AXIS.
 INPUT VARIABLES 09o
INPUT CONSISTS OF THREE CARDS FOR EACH CASE RUN - SETS OF .THREE
CARDS CAN BE STACKED FOR MULTIPLE RUNS
 CARD NO. 1
The FIRST INPUT IS AN INDICATOR AS TO THE NUMBER OF RUNS
REMAINING AFTER PRESENT PUN. ZERO OR BLANK INDICATES LAST
RUN OR A SINGLE RUN. FORMAT (14)
THE REST OF THIS CAPO CONTAINS ANY INFORMATION THE USER WISHES
TO HAVE PRINTED OUT AT Thf TOP OF EACH PAGE.
1.

1.
2.
3.
4.
5.
6.
7.
H.
1.
2.
CARD NO.
FORMAT (8F10.0)
DISCHARGE VELOCITY (M/S)
AMBIENT VELOCITY - 
AMBIENT SALINITY AT DISCHARGE LEVEL IPPT)
LINEAR AMBIENT TEMPERATURE GRADIENT (DEG C/M)
LINEAR AMBIENT SALINITY GRADIENT (PPT/M)
~ CARD NO. 3 
FORMAT
(6F10.0)
DISCHARGE DIAMETER (M)
HORIZONTAL 01SCHAriGF ANGLE (DEG) RELATIVE TO THE AMBIENT CURRENT
90 DEG IS FOR CURRENT NORMAL TO LINE OF DISCHARGE PORTS
MODEL IS NOT GOOD FOR LARGE VARIATIONS FROM THIS DIRECTION
DISCHARGE ANGLE (DEG) RELATIVE TO THE HORIZONTAL.
90 DEG IS VERTICAL - ZF.RO IS HORIZONTAL DOWN STREAM -
DISTANCE ALONG PLUMP CENTFRLINF SIMULATION IS TO STOP (M)
DISTANCE BETWEEN PORTS (M) .
USE A LARGE VALUE TO SIMULATE SINGLE PORT DISCHARGES
DEPTH UF DISCHARGE (M)
EXTERNAL DF.RIV.OUTP
PFAL I I 1.112.I 13.I 14
DIMENSION PRMT(S). Al IX(16.6) . F(6). FP(6) .Nll(lH)
COMMON /AMH/ TIZ.CIZ.CIO.TIO.PIO.GA.AL.RE.FR.Rl.FRL
COMMON /OUT/ L.TH1.Th2.H,Z.SNEW.OS.R.K.N11.J3.J2.J5.JJ.SPACE.S1GMA
XDO.CD.A7.XXX
COMMON /PLTE/ Jb.J7.J8.J9.J10.XlM.YlM.X2M.Y2M.Y3M
62

-------
000& E>KovaJc)
00073

T T 3 = 0.4S1.07034  OO O ^"?
I!4 = 0.31558*1.094 1
00074

00075
1
oFAD (5.102) J.I,Nil
0007b

TF(EOF(5).NE.O) GO TO 221
00077

.J5 = 0
00078

RE = 1.4253E-4
00079

GA=-7.6669E-4
00080

cn=2-
00081

ns=4.
00082

Ml = l
00083

jj = l
00084

.11=0
00085

J2 = l
00086

,J6=0
00087

11 = 1
00088

M = 16
00084
C Nil CONTAINS ALPHAMERIC INFO WHICH IS PRINTED IN OUTPUT. J3 IS ID NO.
00090
102
FORMAT (14,16A4.2A2)
00091
220
PRINT 200. N11.J3
0009?
200
FORMAT (1H1, SOLUTION TO MULTIPLE BUOYANT DISCHARGE",
00093

1 PROHLEM WITH AMBIENT CURRENTS AND VERTICAL GRADIENTS",//,
00094

? P0X.18A4," 1.0. = "14.//)
00095
42
CONTINUE
00096
37
RFAD(5.lOSHJO.UA.TO.CO.TAO.CAO,TIZ,CIZ
00097
105
FORMAT (8F10.6)
00098

RFAll (5,100) DO, TH1,TH?,SF,SPACE,H
00099

PRINT 2000,UO,TO.CO.r>0.SPACE,H
00100

PRINT 200l,UA,TAO,CAO
00101

PRINT 2002,TIZ.C1Z
00102
2000
FORMAT (1M ,"DISCHARGE VELOCITY = "F5.?."-M/S  TEMP. = ",F5.2
00103

x."-OEG C  SALINITY = ,F5.2,"-PPT",/,1H ," DIAMETER = ",F5.2,
00104

K.I-M SPACING = "F5.2 M oo DEPTH =  , F 7 . 2 , "-M")
00105
2001
FORMAT (1H ,"AMBIENT CONDITIONS AT DISCHARGE ELEVATION, VELOCITY
00106

X = "F 7.2 "-M/S  TEMP. = ",F 7.2,"-DEG C  SALINITY = ",F7.2,"
00107

XPPT")
00 I 08
2002
FORMAT (]H ."AMblENT STRATIFICATION GRADIENTS TEMPERATURE
00109

X  ",F10.6,"-0EG C/M  SALINITY  ",F 1 0 .6,"-RPT/M"//>
00110

h = H/D02.
0011 1

CALL SIGMAT(PIO,TO,CO)
001 1?

CALL SIGMAT(PAO,TAO,CAO)
001 13

POVO= 05D0/U0
00114

P=UA/UO
00115

BF=flETO/PIO
00116

GA=GA CO/P10
00117

PI0=ABS(PA0-PI0)
00118

FR = SQRT (PII)/PI0D09.8>
00119

FR=U0/FR
00120

TIZ=TIZDO/<2.T0)
00121

CIZ=CIZD0/(2.CO)
00122

TIO=TAO/TO
00123

CIO=CAO/CO
00124

SF = SF/DO2.
00125

SPACE = SRACE/DO
00126
100
FORMAT (aFl0.6)
00127
9876
R1 = R
00128

RE = 1.0 E*6
00129

PTO = 1.0 ~(8E(1.0-T10) *GA(I.fl-CIO))
63

-------
U.O 1 JU

= -Ht11/ - GACI/
oom

PRINT 207. TI0.CI0.PI0,TIZ,C1Z,PIZ
0 0 132
207
FORMAT (1H AMBIENT CONDITIONS AT DISCHARGE
00133

1" NONI) JMENSIONAL"./, " TEMP=", F H . 5 . 2X , "S AL INI T V = " , F8 . 52X,
00 1 34

? 'nENSTTY = ".F8.b3X.HTFMPGR = ,,,F9.6.2X."SALGR = ",F9.6,2X,
00135

3 '.'DFNGR = ", F 1 0 . 7 ,/)
00136
40
PRINT 202. SPACE
00137
202
FORMAT ( 1H ,> PORT SPACING L/D = ",F6.2)
00138

IF (ABS(PIZ) - 1.0E-10) 90.90.91
00139
90
T = 1.0 E*8
00 140

GO TO 92
oom
91
T = (PIO-l.0)/(-P1Z)
001 <.2
92
PRINT 206. FR.R.T
00 1 o
206
FORMAT (1HU." FROUDE NO ='.F9.24X."VELOC1TY RATIO =".
00144

1 F7.3.4*."STRATIFICATION NO = ",F7.1./)
00 145

FR = 2.FRFR
on 146

SPACE = 2.SPACE
00147

nsn = OS
00148

TTT1 = THl
00149

TTT2 = TH2
0(1150

R.3 = RSIN(0.01745329TH1)COS(0.01745329,,TH2)
00151

IF (R3-1.0) 60.61.60
0015?
61
93 = 0.0
00153
60
CALL ZFE
00154
C ZFF TAKES CARE OF COMPUTING IC AT END OF ZFE.
00155

OS = DS11
00156

THl = TTT10.01745329
00157

TH2 = TTT20.01745329
00158

TIE = T1E0.01745329
00159

T2F = T2E0.01745329
001 60
C CONVFRT TH FROM DEGREES TO RADIANS
00161
C TH
ANGLE AT S=0. TE = ANGLE AT SE
00162

SI = SIN(TIE)
00163

C? = COS(T2E)
00164

PS 12 = R'Sl^C?
00165

IF (RE - 1.0H4) 20,21.21
00166
20
PRINT 214, RE
00167
214
FORMAT (1H FLOW MAY NOT HE FULLY TURBULENT, RE = ",F8.0,/I
00168
21
PRMT(1) = SET
00169

DRMT<2> = SF
00170

RMT<3> = LIS
00171

PRMT(4)= 0.001
001 72
36
1 = -1
00173

N = 6 .
00174
C L=-
1 USED TO INITIALIZE OUTP. N = NUMBER OF DFFERENT1AL EONS.
00175

= 0
00176
C K COUNTS NO OF INTEGRATIONS.
00177

no 2 I=l6
001 78
2
FP(I) = 0.1666667
00179

STOMA = HUE
001 fiO

PRE? = BHd'RBE
001 ft)

F(l) = RRE2CI I 1~(.5-1 I 1)RS12)
001R2

F(2)=OTEt3HE3&( I IP* fl Il-I
001*3

F<3> = l)CEortRE2 (112* ( 1 Il-I I2MRS12)
00184

F(4)=A8E2II2(1.-RSI 2)O02*2.<> (1.-RSI2)RS1?8HE2111
00185

1 ~RS12RS12RHE2.5
00186

F(5) = TIE
00187

F (6) .= T2E
00 188
C SET INITIAL CONDITIONS AT END OF POTENTIAL CORE.
00 189

PRINT 208
00 1 90
208
FORMAT (1H0.20X,"Z0NF OF ESTAbLISHEO FLOW',/
00191

1 1H  5x , "S", 8X,"X",7X,'
-------
uoiso


Ih 1 IMLt- -111 1 Ut 1 1 , I 1
00197

1 1
CALL EXTR {S.X,Y,Z.B.PU,DT,5.D0)
00198

10
PR INT 212, K,IhLF.PMMT(4),PMT(5),PRMT(3)
00199

212
FORMAT (1H NO OF INTEGRATION STEPS=,14.3X."NO OF HALVINGS"
00200


1 I3.3X,'>ABSEHR = ",F9.5,3X."0PMT (5) =",F9.3." FINAL SPACE ="
30201


2 .F9.3)
00202


GO TO 1
J0203

221
STOP
00204


FND
0020S


SUBROUTINE SIGMAtIrho.T.S)
90206


SIGO=-4.fl2K-4>oS.Rl49>S-.O93
'10207


C = 1 .E-6T( (.01667T-.8164)
00226
C
PRFLF COMPUTES INTEGRAL OUANITITIES WITHIN THE JET
00227
C
F8
IS THE INTEGRAL OF DENSITY DIFFERENCE, AND 8 THE WIDTH OF TE JET
00228


FP<1) = E
00229


FP(2)=-F(1)TIZS2 ~ VTI
00230


FP<3) =-F(1)CIZS2 ~ VCI
00231


FP (4) = E*RS1C2 FPS2 ~ VUI
00232


F9 = F(4) - 0.25E0E ~ VVI
00233


TF < BBS(C2) - 0.01) 1.1.2
00234

1
F (5) = F(6)
00235


FP(5) = 0.0
00236


GO TO 3
00237

2
FDY=0.
00238


UNX=-(C22C1S1 I
00239


UNY=S22C202*C12
00240


I)NZ = -(S2C2S1)
00241


UN=SORT
00247


FP(5)= (ERC1)/(C2F9)*FDY
00248

3
FP(6) = (F8C2 - (EM*FD)SlS2)/F9
00249


K =K ~ 1
00250


RFTURN
0 0 251


FND
00252


SUBROUTINE XINT(I,A.F1INT.F21NT.F3INT)
00253


DIMENSION XM25) , X 11 (2b) ,X 12(25)
00254


DATA (XA=1.0E10,1.,1.1.1.2,1.3.1.4,1.5.1.6,I.7,1.6.1.9,
00255


X2.0,100.),(X I 1 = .4314,.4314..42^5,.4259,.4238,.4221,.4209,
00256


X.4201,.4196,.4194,.4193..4193), (X I 2=.436..3905..3576,
00257


X.3348,.3201,.31 14,.3068..304 7,.304..3038,.30377, .303?7)
0025H

1
TF ( A-XA.I I ) ) 2,3,4
00259

2
I = 1-1
00260


GO TO, 1
00261

3
F1 I NT = X 11 (I )
65

-------
.>0262
10 263
00264
D0 26S
00266
00267
0026B
00269
J0270
J0271
0027,?
:)02 7 3
(10274
00 275
0O27t>
00277
00278
00279
0 028 0
002B1
00282
002H3
00284
00285
00286
00287
0028H
00289
0029 0
00291
0029?
00293
00294
00295
00296
00297
00298
00299
00 300
0030]
00302
00303
00304
00305
00306
00307"
00308
00309
00310
00311
00312
0 0313
00314
00315
00316
0031 7
00318
00319
00320
003? 1
00322
00323
00324
00325
00326
00327
F 2 1 N I = A I 2 ( 1 )
r,0 TO 8
4	IF(A-XA(I*1)>5.6.7
6	FlINT=XI1{I*1)
F2INT = XI2(IM>
T = I ~ 1
r,0 TO 8
7	T = i * 1
RO TO 4
5	clINT = XIl(I)*(XIl(l*l)-XIl(I> )/(XA(1*1)-X A(I) ) 
-------
00328
00329
1)0330
00331
0033?
00333
00334
00335
00336
00337
II033H
00334
00340
00341
00342
00343
0034 4
00345
00346
00347
0034B
00349
00350
003S1
0035?
00353
00354
00355
00356
00357
00358
00359
00360
00361
0036?
00363
00364
00365
00366
00367
0036H
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
003HO
00381
0038?
00383
00384
00385
00386
00387
00388
00389
00390
00391
0039?
00393
f4=4.1 1 J* A I
C5=U12U12AI3/PI
C6=2.II4AI2/P1
C7=2.II3U12AI1/PI
r,0 TO 504
C oooo mfrgeo PLUME
503	CALL XINT(I.A.Al1,Al?,AI 3)
C1=.45AI3/PI
r?=U12Al3/Pl
C3=.315SBAI3/P1
CW.#U12C1
C5=U12>C2
C6=C 3
C7=0.5C4
504	IF(AHSIU12)-.00]>112,112.113
112	RSUR=F(1)F(1)C3//Cl-C3Z.BC2"F(1)/C1/C1-F(4)
CC=C3F(1> F(1)/Cl/Cl
RSQR=-BR-SURT 
nSQR=BSQR/(2.AA>
12 =SORTIBSQP)
MGMA = B
DLUM=F(1)/(BSQRC1)-C2/C1
IJCL = DLUM*U12
DTCL = F(2)/(HSQH(01 UMC6*C7))
DCCL = F (3)/(HSiiR (Ol UMC6*C7) )
PICL = REDTCL ~ GADCCL
IF I 11/ ((PI 0-1.)FR>
IF(A.LT.2.)FH=2.BS0PIJ3P1CLAI1/IPIFR(PIO-1.))
IF(A.LE..95)FH=FaAI3/AI1/2.14068
C	COMPUTE F8 - INTEGRAL OF DENSITY DIFFERENCE
FRLI = PICLB/(FR(PIO-1.0)UCLUCL)
GO TO 3
2 F8 = 0.0
FRLI = 0.0
3 IF < AB5(FRLI) - 1 0.0ARS(1.O/FRL)) 10,13,13
13 FRLI = 1. O/FRL
10 CONTINUE
Al = .05
A?2 = 0.0
AA = Al ~ A22FKLISIN(F (6) ).
CC = A7WSQHT <1.0-S12S1?>
RH=RABS(UCL-RS12)
28 TF (A.GE.2. )F = AA (BB (1 .-XXX/A) ~CCH)
IF(A.LT.2.)E=AA
RF AL 111.112,113,114
DIMENSION F(6).FP(to),PRMT(5),SO(4),SN(4).N11(18)
COMMON /OUT/ L,TH1,Th2,H,Z,S1,DSR,KN11,J3J2,J5JJSPACE~SIGMA
X,00,CO,A 7,XXX
COMMON /AMB/ T IZ , CIZ ,C10,T 10P10,GA,ALBEtFR,R1,FRL
COMMON /PLTE/ Jb,J7,J8.J9,J 10,X1M,Y1M,X2M,Y2M,Y3M
COMMON /ZF9/ SEr,XEYEZE,BBEDCE,DTE,LCNTRUE,RTEWCE, T1ET2E
COMMON /STR95/ II.M.R3
COMMON /CONS/ 111,112,113,114,111
COMMON/WGM/MERGCK,ROVO
I = I I I
PII=3.14159
67

-------
00394
0039s
00396
00397
0039H
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
004 1 0
U04 1 1
00412
0041 3
00414
00415
004 16
0041 7
0041ft
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
oo43a
00439
00440
0044 1
00442
00443
00444
00445
00446
00447
00448
00449
004SO
0 04b 1
0 0452
00453
00454
00455
00456
00457
00458
00459
IF (L) 6.8.7
6 L= 0
L 19 = 0
Ml =-1
M2 = 0
71 = -OS
SNEW = PWMI(1)
SN (1) = SIN(TIF)
SN<2) = 51N ( T 2 E)
SN(3) = COS(Tit)
SN(4) = COS(T2E)
TIME = SET
x = XE
Y	= YE
7. = ZE
urc = l.o
S 1 = PRMT <1>
T 3 I = T 10
C31 = CIO
IF(CIO-1.O.LT..000001)C3I=0.
IF (TIO- 1.0) 63.64.63
64 T 31 = 0.0
63 GO TO 10
8 . I. = l
SNEW = OS ~ SI
SO(1) = SN(1)
SO(2) = SN(2)
SO(3) = SN(3)
SO(4) = SN(4)
7 TF (M2) 18.18.19
19	ns = S-SNEw
TF (S-Sl - 1.0) 83.1,1
18 IF (S - SNEW) 20.1.1
20	IF (46S(S-SNEW) - O.OOl'SNEW) 1.1.83
1 SN(1)=S[N(F(5))
SN(?)=SIN(F(6))
SN(3)=COS(F(5))
SN(4)=COS(F(6))
X= X ~0.5 S0 < 4 ) )  ( S-S 1)
Y	= Y ~0.5 ~ S0(1)S0U)->
71 = Z
7 = Z~05 (SN ( 2) ~ S0(?) IMS-51)
no 9 1=1.4
9	SO(I)=SN(I)
10	PIZ = -f)ET 17. - GAC I Z
TI = TIO ~ TIZZ
CI = CIO ~ CIZZ
PI = Plo ~ PIZ*Z
IF (F ( 4 ) ) 21.21.55
55 Rl4=RSN(1)SN(4)
&=SPACF/SIGMA
IF(A-2.I 501.500.500
501	IF(A-.95)503.503.502
C SINGLE PLUME 
500 r1=1 I 1
C?=14/2.
C3=I 12
C4 = 2.1411 1
r5=R14R14/2.
C6=I I 2
C7 = II 1R14
0,0 TO 504
C	MERGING PLUMES 00
502	CONTINUE
519 FORMAT(/20X."ZONE OF MERGING FLOW'/).
IF(MEHGCK.EO.O)PRINT 519
68

-------
00460


MFRGCK=1
00461


CALL X1NT(I.A.A 11A IA I 3)
00462


rl=2.II3Al 1/PII
00463


C2=R14A13/PI 1
00464


C3 = 2.MI4AI2/PII
00465


C4=4.II3AI1R14/PII
00466


r5=Rl4Rl4Al3/PII
00467


C6=2.II4AI2/PII
00466


r7=2.II3Rl4Al1/PIT
00469


r,0 TO 504
00470
C
oo MERGED PLUMES 
0047 1

503
CALL XINT(1,AAI1,AI?,AI3)
00472


ri=.4SAl3/PII
00473


C2=W14A13/P11
00474


r3=.3155HAI3/PII
00475


C4=2.R14C1
00476


C5=R14C2
00477


C6=C3
00478


07=0.5C4
00479

504
IF(A0S(RU)-.OO1) 12.12.99
00480

12
BSOR=F<1> F < 1 >C3/(F(4)C1C1>
00481


GO TO 15
00482

99
AA=C3C2C2/C1/C1-C4C2/C1*C5
00483


RR=C4F( I)/C1-C32.C2F(1>/C1/C1-F<4>
00484


CC=C3F(1)F(1)/C1/C1
00485


BSOR = -BB-SURT(H8RB-4.AACC)
U0486


RSOR = BSQW/(2.AA)
00487

15
B = SQKT(BSQR)
00488


SIGMA = B
00489


S1GMA0 =SI(iMADO
00490


DLUM = F(1)/AOELC I 1 .O-CIO)
00502
c
DELU= (UCL - UI*51C2)/UJO - UI> DELT = ( TCL - T 1) ~ ( TO - TI)
00503


F55 = F(5)57.29578
00504


F 66 = F(6)57.29578
00505
c
CONVFRT ANGLES BACK TO DFGRFES
00506


DMA = 0.5(UCL ~ UCCl
00507


TIME = TIME ~ (S-Sll/UMA
00508


TF (FR -0.98E*6) 110.111,111
00509

110
IF (ABS((PI-PCL)0> -1.0E-6ABS(UCL2#FR(PIO-1.0)1) 112,113,
00510

113
FRL = FR(PIO-l.b>UCLUCL/
00511


fiO TO 112
00512

111
FRL = FR
00513

112
IF (L19) 70.70,71
00514

70
SOO = S/2.D0
00515


XO = X/2.D0
00516


YO = Y/2.n0
00517


ZO = //2.00
00518


P0=2.F(1)
00519


TTMER=TIMEWOVO
00520


PRINT 104,SOO,XO,YO,70.F55.F66,SIGMAO,DELU,DELT,DELC
00521

PI,TIMER,HO
00522

104
FORMAT (1H , F9.2.1X.5F8.2, 1X,F8.2,3F9.3,1X,F9.5, 2F9.2)
00523


Ml = Ml ~ 1
00524


M = M 1
00525


1.19 = L19 ~ 1
69

-------
00526


go iu i d
00S27

71
1.19 = 0
00528

72
PICL = PI-PCL
005?')


PICLO = PIO-l.O
00530


IF C M 21 40*40*4 2
00531

40
IF (SIGN(l.O.PICL) - SIGN<1.O.PICLO)) 41,42.41
00532

41
ns = 0.0
00533


Ml = 0
00534 .


M2 = 1
00535


MFXT=1
0 05 36


PRINT 213
00537

213
FORMAT (1H . 10X."PLUME HAS REACHED ITS EQUILIBRIUM HEIGHT
00538


1  STRATIFIED ENVIRONMENT")
00539

42
TF (Ml - 10) 36.32.3?
00540-

32
Ml = 0
00541


OS = 2.0DS
00542

36
TF (Z) 3.3.31
00543

31
IF (M2) 90,90.91
00544

91
TF (Z-Zl) 21.90.90
00545

21
PRMTI5) = 1.0
00546


TH - 11
00547


PRINT 212
00548

212
format <1H .9X." PLUME HAS REACHED ITS MAXIMUM HEIGHT - " ,
00549


1 "STRATIFIED ENVIRONMENT"./)
00550


IF  83.83.81
00551

90
IF (7.-H) 3.3,4
00552

4
PRMTI5) =1.0
00553
C
thf
CARDS BELOW ARE FOR. USE WHEN THE PROGRAM IS USED FOR
00554
C
DISCHARGE into WATER OF DEPTH H. This MUST HE READ INTO PROGRAM
00555


PRINT 214
00556

214
FORMAT (1H ."PLUME HAS REACHED.WATER SURFACE"./)
00557
C
IF
WF HAVE REACHED THE WATER SURFACE - STOP.
00558


TF (J 6 > 83.83.81
00559

3
SNEW = SNEW ~ DS
00560


SI = s
00561


UCC = UCL
00562


TF 
-------
OOSId
00593
00594
00595
00596
00S97
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
Sl(l) = bl t 1 *1>
Xl(I) = X1 (I~1)
rlili = v l (I~l)
71(1) = 21 (1*1)
R1(I) = 81(1*1)
mil (I) = ou 1(1*1)
DT1I1I = 011(1*1)
5 CONTINUE
M=3
1 SI = Z
B1 (M> = B
nu1(m > = ou
DTI(M) = OT
MR = M
RFTURN
4 IF (MR-3) 202.6.202
6 MR2 = MR-2
MR1=MR-1
7MAX = Zl(MRl) * DU1(MR1)/(2.0C))D0
207	PRINT 102.ZMAX. ZM
102 FORMAT 1 1 HO." MAXIMUM PLUME HEIGHT. EXTRAPOLATED
1 "LINEAR ="F8.2" PARABOLIC =".F8.2./>
RFTURN
202	PRINT 203
203	FORMAT <1H .ERROR IN HPCG. IHLF = 11''./)
RFTURN
FND
SUHHOUTINE ZFE
THIS PROGRAM COMPUTES ThE SOLUTION TO BUOYANT JET PROBLEMS WITHIN
THF ZONE OF FLOW ESTABLISHMENT. I.E. UP TO XE.
THF SOLUTION CALCULATES RU.HT.RCB.TH1.TH2 AS FUNCTIONS OF S.
DIMENSION PRMT(5).AUX(16.6).F(6).FP(6).N11(18)
COMMON /AMB/ TIZ.CIZ.CI0.TI0.PI0.GA.AL.BE.FR.R1.FRL
COMMON /OUT/ L.TH1.TH2*H,Z,SNfw,DSHKNl1.J3J2J5.JJ.SPACESIGMA
X.DO.CD.A7.XXX
COMMON /ZFsi/ SET.XE.YE.ZE.0HEr>CE.DTE.LCNT,RUERTE.HCE,TlE.T2E
EXTERNAL DERIV1OUTP1
PIZ = -BETIZ - GA'CIZ
TH1 = 0.01745329TH1
TH2 = 0.01745329TH2
PRMT(1) = 0.0
PRMT(2) = 100.0
PRMT(3) = 0.01
PRMT(4) = 0.0001
TF (ABS(COS(TH2)) - 0.01)1.1.2
1	TH1 = TH2
2	N=6-
L =0
LCNT = 0
71

-------


no 41=1,6
00659
3
FP(I> = 0.16666667
00660

F(l) = 1.0
0066]

F(2t = 1.0
00662

F(3I = 1.0
00663

F 14) = 0.0
00664

F(5) = TH1
(10665

F (6) = IH2
00666
C Ffl) = HU. F12) = Hi. FI3> = HC  F<4) =H, F<5> = Ihl. F<6) = TH2
U0667

rm_F = 0
00668

PRINT 210
U0669
210
POftMAT 11H0.20Xi"ZONF OF FLOw HSTARLISHMENT  ALL LENGTHS "
00670

p'ARE IN MEIERS ',/
00671

1 1H .5X.''S",9X.''X",7X. Y''.7X."Z".6X,"TH1''.5X."TH2',.9X
00672

1  "R I,SX."BU"6)U"WT".6X."RC,,.5X "OUCL" .4*"DT CL" 4X t "ijCCL" 
00673

2 6X."Pl"7<"TIME"./l
00674

CAUL 0UTP1 (PHM r(l)F.FPlHLf"  N PKMT )
01)675

TALL HPCG (PHMT,F,FP,N,IHLFil)ERlV1.0UTPl.AUX>
00676

TF  . FA =B/AO
006B4

TH1 = F 15)
(10685

TH2 = F(6)
00636

no 9 1=1.6
006H7
9
FPU) = 0.1666667
OO6R0

PRMT11) = SET
00689

PRMT12) = SET ~ 100.0
00690

PRMT 13) = 0.01
0069 1

PHMT(4) = 0.000 i
00692

1HLF = n
00693

("ALL 0UTP1 (PHMT 11) .F.FP.lHLF.N.PRMD
00694

CALL HPCG )
00718

R12 = RS1*C2
00719

SI 2- = S1C2
00720

nTO=1.0-TJO - TIZZ
00721

OCO=1.0-C10 - CIZaZ
00722

F = El (0.02040.0144,1F (4) 1 o (1.0*E2iiS2/FRI o (A8SC1 ,0-H12> >t"3aH
00723

1 'SORT(1.0-S122>)(1-E4/SPACE)
72

-------
00724
C
E
IS THE ENTRAPMENT FUNCTION.
007?S


01=.45*.55R12
00726


011=2.*<.12857*.37143*R12)
00727


02 = .31558*.13442*R1?
00728


n3=2.* t.06b76..06181R]2)
00729


n4=.31558*.26fl04R12. .4155R12R12
00730


05=2. (.066 76*. 12362*R12*.3nS)62R 12H12)
00731


flf> a DTOBEMF2)*F ~ .4SF (2> F U) >
00732


nsaOe* OCbGA(Fl3FI3)O.S *0.i257F(4)*F(4 ~45*Ft3>fUH
00733


0 = 0.5F (1) F (1) ~ F(])F(4)[)ft *F(4)F <4)*0.5i>5 -0.25EE
00 734


A A ] =0.5F (1 ) F<] ) ~nlF<] FU) ~ D11'*Q.5*F I 4) u F (4)
00735


A (1.1 > = F (1 ) ~ 01F(4)
00736


A(1.2>= 0.0
00737


A(1.3) = 0.0
00738


All.4) = D1F < 1) ~ D1]FUI
007 3")
C
FIRST EQUATION IS CONTINUITY.
00740


A(21) = 0.0
00741


AI2.2) = DT0F(2) ~ n2DT0F(4)
00742


A<2.3> = 0.0
00743


A(2.4) = PT0*(D2F(2> ~ D3*F14))
00744
c
SECOND EQU&TION IS ENERGY.
00745


A(3 ~ 1> = 0.0
00746


A (3.2) = 0.0
00747


A ( 31 3 > =" DCO* (F ( 3) ~ .D?F(4)
00748


A (3.4) = DC0(D2*F(3> ~ D3F < 4))
00749
c
THIRD EGUATON IS CONCENTRATION.
00/50


A(4 * 1) = F < 1) '~ D4*F(4)
00751


A ( 4.2) = 0.0
00752


A(4 * 3 > = 0.0
00753


A I 4.4) = 04F (1 > ~ DSF(4)
00754
c
FOURTH EQUATION IS MOMENTUM
00/55


FP(1) = E
00 756


TIZ7 = TIZ
00757


FP ( 2) = TIZZ*S21-AAI0.5*F(2)? ~ F(2>FI4)D2 ~0.5*F(4)*2D3)
00758


FP(2)=FP(2)~VTE
00759


FP(3) = CIZS2*(-AAI ~0.5*F(3)2 *T(3)F(4)D2 ~0.5*F(4)*2*D3>
00760


FP(3) = FP< 3> +VCE .
00761


IF (L-2) 3.4.3
00762

4
A(2.2 > = O.S03F(412
00763


A<2.4) = 03F(2)F(4)
00764


A(3.3) = A(2.2)
00 765


A(3.4) = U3*F(3)*F(4)
00766


FP(2J=-TIZZ*S2*AAI*VTEnTE
00767


FPO) = -CIZS2AAI*VCE"DCE
00768


06 = 
-------
00790	A?(2.2)=A(4.4)
00791	FP(2)=FP(4)
00792	CALL SIMQ(A2FP.JK)
00793	GO TO 7
00794	20 .1 = 3
00795	A3(1.1)=A(1.1)
00796	A3(l2)=0.
00797	A3I1.3)=A(1.4)
00798	A3(2.1>=0.
00799	A3(2.2)=A(2.2)
00800	A3I2.3)= A(2.4)
OOflOl	A3(3.1)=A(4,1)
1)0802	A3 (3.2) =0.
00*03	A3(3.3)=A(4.4)
00804	FP(3)=FP(4)
00805	CALL SIMQ(A3.FP.J.K)
00806	GO TO 7
00807	6 .1 = 3
00808	A3(1.1)=A(1.1)
00809	A3(1.2>=0.
00810	A3(1.3)=A(1,4)
00811	A3(2.1)=0.
00812	A3(22)=A(33)
00813	A 3 (2. 3) = A (3 . 4 )
00814	A3(3.1)=A(4,1)
00815	A3(3.2)=0.
00816	A3(3.3)=A(4.4)
00817	FP(2)=FP(3)
00818	FP(3)= FP(4)
00819	CALL SIMO(A3.FP.J.K)
00820	7 IF (K) 1.2.1
00821	2 TF(J-3)131.132.136
00822	131 rP(4)=FP(2)
00823	FP(3)=0.
00824	FP(2)=0.
00825	GO TO 136
00826	132 TF(DCO.NE.O.)GO TO 134
00827	FP(4)=FP(3)
00828	FP(3)=0.
00329	GO TO 136
00830	134 FP(4)= FP(3)
00831	FP(3)=FP(2)
00832	FP(2)=0.
00833	136 UNX = -(C22C1S1)
00834	l.iNY = S22*C22,tCl?
00835	UNZ=- = ERC 1 /(0C2)+FDY
00843	FP (6) = (F8C2 - (ER+Fn>S1S2>/Q
00844	RFTURN .
00845	1 WRITE(6.100)
00846	100 FORMAT (1HO."SINGULAR MATRIX-EXECUTION TERMINATED AT THIS POINT")
00847	STOP
00848	FNO
00849	SUBROUTINE 0UTP1 (SF,FP.IH.N,PRmT)
00850	COMMON./AMB/ T 11. C IZ . C 10. T 11),PI0 .GA . AL . BE. FW. R 1. FWL
00851	COMMON /OUT/ L . IH1 , TH2.M. Z . SNEhl. OS. R .K , N1 1 . J3. J2. J5. JJ. SPACE. SIGMA
00852	X.00.CO.A 7.X XX
00853	COMMON /ZF9/ SET.XEYE.ZE.HBE.DCE.OTE.LCNT.HUEKTE.RCE.TIE.T2E
00854	COMMON /STR95/ II.M.R3
00855	COMMON/WGM/MERGCKROVO
74

-------
00857
0085B
00659
00660
00661
0066?
00863
00864
00865
00866
00867
00666
00669
00870
00871
00672
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00688
00889
00690
00691
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
DIMENSION f (6)lb).PHMI lb).bU14 IbN(tI.Ml I IItt)
IF(L) 1.1.2
1	SNtl) = SIN(THl)
SN<2) = SIN < TH2)
SNI3) = COS(TH1)
SN(4> = COS(TH2)
TJ = PRMT(l)
r)S = PRMTI3)
BUN s 0.0
RTN=0.0
RCN = 0  0
RN=1.0
SSN=0.0
T1N=TH157.29576
T?N=TH257.29578
XN = 0.0
YN = 0.0
7N = 0.0
X =0 . 0
v = 0 . 0
7 = 0.0
JP = 100
IF  30.31 .30
31 L = 1
GO TO 6
30 L = 2
X = XE
Y	= YE
7 = Z E
GO TO 6
2	IF (ABSIS-TJ >-0.01) 3.3.4
3	SN <1> = S1N(F < 5))
SN(2> = SIN(F(6))
SN(3) = COS IF(5))
SN(4) r COS(F(6))
X = X ~ 0  5* (SN ( 3) SN < 4 ) ~ SO I 3>S0(4))DS
V	= Y ~ 0.5MSN(1)*SN(4) ~ SOI] )S0(4) )r)S
7=2* 0.5*(SN(2) ~ SO(2))DS
6 no 5 1=1.4
5 SOI I) = SN 11)
TIME = S
PTZ = -BETIZ - GACIZ
CI = CIO *CIZZ
TI = TIO*TIZZ
PI = PIO ~ PIZZ
F55 = F(5)57.29576
F66 = F (6) <>57.29578
TF (JP-100) 50.51.51
51 DUU = (1.0-RSN(1)SN<4)I/I1.0-R3)
DELT = (l.O-Tll/ll.O-TIO)
DELC=0.
nFNOM=l.-CIO
IFIOENOM.NE.O.)DELC=(1,0-CIl/OENOM
C Ft?) = RT OK (TM-TI) /TO. F(3) SAME F"OP C
C DUU = IUO-U12)/(UO - U120 >. DELT = ITM-II)/ITO-T10)
IF (L-2) 22.23.22
23 I1ELT = F I 2)/ 11  0-T 10)
nLC=0.
nENOM=i.o-cio
IF IDEN0M.NE.0.)DELC=FI 3)/DENOH
100 FORMAT (1H .F9.2lX,SFfl.2.2X2F8.3.16X.3F8.32X F9.5.F9.1)
22 SOO = S/2.00
xoo = x/2.no
YOO = Y/2.l)0
ZOO = Z/2.D0
RRO = F(4)/2.D0
75

-------
00922

RUOO = F( 1)/2."00
00923

riTOO = F 121/2.*00
00924

ncoo = F<3)/2.DO
00925

TTMER=TJMER0V0
00 926

TF (1-2)1001,1002,1001
00927
1002
PRTNT 10n.SOO,XOO,YOO.ZOO.F55Fb6HRO.RUOO,r>Ull.DELT,
00928

IDFLC.PNTIHfR
009?9

C,0 TO 1004
00930
1001
PRINT 100 31SOO,XOO. YOO, 700,F55,F66,BRO,RUOO.OT 00.OCOO,UUU
00931

lOELT.OELC.PI.TIME
00932
11)04
RPI_OT=F(4)2.S0RTIP.)
00933
1 003
FORMAT (1H .F9.2,lX.SFfl.2.2X.7K6.3,2XF9.5,F9.1>
00934
C
WRITE(7.1000)SY,X,Z.BFLOT,F 11>F<2),F(3>OUUOELT tOELC.
00935
C
TIMEW
00936
C1000
FORMAT(5F8.3.6F5.3.FQ.3)
009 37

M = M ~ 1
00 93ft

TF MYO-YN)/I T 200. SOO
00984
2on
FORMAT (1H STARTING LENGTH, T =",F9.3,)
00985

(. = -2
00986

PRMTI5) = 1.0
0048 7

pfturn
76

-------
0098B

13
TJ = TJ ~ PRMT < 3)

009B9


RETURN

00990

40
TJ = TJ ~ PRMT(3)

00991


RETURN

00992

21
IF (Fill) 10.10.40

00993

10
PRMT < 5) = 1.0

00994


SFT = SSN - RUN"(SSN-SS>/(RUN-RUOI

00995


RUE = RUO ~ (SE TSS1(RUO-HUN)/(SS-SSN)

00996


ORE = BO ~ (SE T-SS)(BORN)/ISS-SSN)

00997


me. = WTO ~ (SEr-SS)(RTO-RTN>/(SS-SSN)

0099ft


nCE = RCO ~ (SE1-5SI

00999


TIE = T10 ~ (SEr-SSIMTlOTlN)/lYO-YN)/(SS-SSN>

01003


7F = ZO ~ (SFT-SS)  (70-7N) / (i,5-SSN

0 1 004


OUT = <1.0-RSIN<0.0]745oTlE)C0S(0.01745FNOM=l .O-CIO

01008


IFIDENOM.NE.O.)OCT =OCE/DENOM

01009


PI = PIO ~ PIZ'ZE

01010


SOO = SET/2.DO

01011


PRINT 205.500

01012

205
FORMAT <1H STARTING LENGTH, VELOCITY = ". F9.31

01013


*00 = XE/2.00

01014


YOO = YE/2.nO

0101S


700 = ZF_/2.DO

01016


HRO = H8E/2.*D0

01017


RUOO = RU/?.*D0

01018


DTOO = DTE/2.

01013


OCOO = OCE/2.

01020


PRINT 1'00SOO.XOO.YOOZOOT1ET2FRBO.RUOUDUTiDTT

01021


lOCT.PI.TlME

01022

4
RFTURN

01023


FND

01 024


SUBROUTINE SIMQ(A.B.N.KS)
SIMO
01025
C


SIMQ
01026
c


SIMO
01027


0 I MENS I ON A (1) ,d( 1)
SIMQ
01028
c


SIMQ
01029
c

FORWARD solution
SIMO
01 030
c


SIMO
01031


TOL=0.0
SIMO
01032


K S = 0
SIMO
01033


JJ = -N
SIMQ
01034


no 65 J=1,N
SIMO
01035


>(Y = J*1
SIMO
01036


JJ=JJ*N.1
SIMQ
01037


fUGA = 0
SIMO
01038


IT = JJ-J
SIMQ
01039


DO 30 I=J.N
SIMQ
01040
c


SIMO
01041
c

SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN
SIMQ
01042
c


SIMQ
01043


T J=IT 1
SIMU
01044


IF(ABS(BIGA)-ABS(A(Ij)>) 20.30.30
SIMQ
01045

20
niGA=A(IJ)
SIMQ
01046


TMAX=I 
SIMQ
01047

30
CONTINUE
SIMQ
01048
c


SIMQ
01049
c

TEST FOR PIVOT Lt'SS THAN TOLERANCE (SINGULAR MATRIX*
SIMQ
01050
c


SIMQ
01051


IF 
-------
01 os4
C


010SS
C

interchange hows tf necessary
01 056
C


010S7

40
T1=J*N(J-2)
010SH


tt=imax-j
010 59


no SO K = J.N
01060


T1 = 11*N
01 061


12= 11~ I T
0106?


SAVE=A( I 1 )
01 063


A(I 1)=A(12)
01 064


A(12)= S A V E
01 06S
C


01066
c

DIVIDE EQUATION Hy LEADING COEFFICIENT
01 067
c


01068

so
A CI I)=A(I 1)/BIGA
01 069


<;ave = h ( I MA A )
01070


rt(I MAX)< J)
01071


R(J)=SAVE/BIGA
01072
c


01 073
c

ELIMINATE NEXT V4RIARLE
01 074
c


01 07S


TF(J-N) SS.70.5S
01076

5S
TQS=N(J-l)
01077


no 65 IX=JY.N
01078


TXJ=IQS* IX
01 079


T T = J - I X
01080


no 60 JX=JY.N
01081


TXJX=N)a(IXJ)>
01085
c


01 086
c

BACK SOLUTION
01 087
c


01088

70
NY=N-1
01 089


T T = N 0 N
01090


no 80 J=1.NY
01091


I A =IT-J
01092


TR=N-J
01 093


TC = N
010 94


no 80 K=1.J
01095


R(IB)=fl(IBl-A(IA)H(IC)
01096


T A=IA-N
01097

80
TC=IC-1
0109H


RFTURN
01 099


FNO
01 100
c


01101
01102
c




SURROUTINE HPCG(PRMT.Y.DEWY.NDIM.IHLF.KCT.OUTP,AUX)
01103
c

SUBROUTINE HPCG
0110 4
c


01 10S
c

PURPOSE
01106
c

TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL
01107
c

DIFFERENTIAL EOUATIONS WITH GIVEN INITIAL VALUES.
0110S
c


01 109
c

USAGE
01 1 10
c

CALL HPCG (PRMT.Y.DERY.NDIM,IHLF.FCI.OUTR.AUX)
01111
c

PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.
01112
c


01113
c

DESCRIPTION OF PARAMETERS
01114
c

PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER
01115
c

OR EQiial TO 5. WHICH SPECIFIES 1 HE PARAMETERS OF
01116
c

THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR
01117
c

COMMUNICATION BETWEEN OUTPUT SUrtROUTINE (FURNISHED
01 llfi
c

HY THE USER) AND SURROUTINE HPCG. EXCEPT PRMT(S) 1
01119
c

THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE 1
MMU
SIMU
S I HQ
SIMQ
SIMO
SIMU
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMO
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMO
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMU
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMU
SIMQ
SIMQ
SIMU
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
78

-------
01
I lit)
C

rIPCG ANO 1 HEY ARf
HPCG
01
121
C
PRMT( 1 )-
LOWER BOUND OF THE INTERVAL (INPUT 1,
HPCG
01
122
C
PHMT(2) -
UPPER ROUND OF THE INTERVAL (INPUT),
HPCG
01
123
C
PRMT(3)-
INITIAL INCREMENT OF THE INDEPENDENT VARIABLE
HPCG
01
124
C

(INPUT).
HPCG
01
125
C
PHMT (4)-
UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS
HPCG
0 1
126
C

GREAfFR Than PRMT(4). INCREMENT .GETS HALVED.
HPCG
01
127
C

if increment is less than prmtoi and absolute
HPCG
01
128
C

ERROR LESS THAN PRMT(41/50. INCREMENT GETS DOUBLED
.HPCG
01
129
C

THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS
HPCG
01
130
C

OUTPUT SUBROUTINE.
HPCG
01
131
C
PRMT<5)-
NO INPUT PARAMETER. SUBROUTINE HPCG INITIALIZES
HPCG
01
132
C

PRMT(5)=0. IF The USER WANTS TO TERMINATE
HPCG
01
133
C

SUBROUTINE HPCG AT ANY OUTPUT POINT, HE HAS TO
HPCG
01
1 34
C

CHANGF PMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE
HPCG
01
135
C

OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE
HPCG
01
1 3*
C

FEASIBLE IF ITS DIMENSION IS DEFINED GREATER
HPCG
0 I
137
C

THAN 5. HOWEVER SUBROUTINE HPCG DOES NOT REQUIRE
HPCG
01
138
C

ANO CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
HPCG
0]
139
c

FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
HPCG
01
140
c

(CALLING HPCG) WHICH ARE OUT AINEO BY:SPECIAL
HPCG
01
141
c

MANIPULATIONS WITH OUTPUT OATA IN SUBROUTINE OUTP.
HPCG
01
142
c
Y
INPUT VECTOR OF INITIAL VALUES. (DESTROYED)
HPCG
01
143
c

LATEWON Y IS THE RESULTING VECTOR OF DEPENDENT
HPCG
01
144
c

VARIARLES COMPUTED AT INTERMEDIATE POINTS X.
HPCG
01
145
c
DEMY
INPUT VECTOR OF ERROR WEIGHTS. (OESTROYEO)
HPCG
01
14b
c

THE SUM OF ITS COMPONENTS MUST HE EQUAL TO 1.
HPCG
0 1
147
c

LATERON OERY IS THE VECTOR OF DERIVAT 1 l/ES. WHICH
HPCG
01
148
c

BELONG TO FUNCTION VALUES Y AT A POINT X.
HPCG
01
149
c
NDI M
AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
HPCG
01
150
c

EQUATIONS IN THE SYSTEM.
HPCG
01
151
c
IhLF
AN OUTPUT VALUE. WHICH SPECIFIES THE NUMBER OF
HPCG
01
152
c

BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS
HPCG
01
153
c

GREATFR THAN 10. SURPOUTINE HPCG RETURNS WITH
HPCG
01
154
c

F.HROR MESSAGE I HLF = 1 1 INTO MAIN PROGRAM.
HPCG
01
155
c

ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE
HPCG
0 1
156
c

PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)
-HPCG
01
157
c

PRMT(1)) RESPECTIVELY.
HPCG
01
158
c
FCT
THE NAME OF AN EXTERNAL SUBROUTINE USED. IT
HPCG
01
159
c

COMPUTES THE RIGHT HANI) SIDES OERY OF THE SYSTEM
HPCG
01
160
c

TO GIVEN VALUES OF X AND Y. ITS PARAMETER LIST
HPCG
01
161
c

MUST RE X.Y.OERY. THE SUBROUTINE SHOULD NOT
HPCG
01
162
c

DESTROY X AND Y.
HPCG
01
163
c
OUTP
THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.
HPCG
01
164
c

its parameter list must be x.y.oery,ihlf.ndim.prmt
.HPCG
01
165
c

none OF THESE PARAMETERS (EXCEPT, IF NECESSARY,
HPCG
01
166
c

PRMT(4).PRMT(5)....) SHOULD BE CHANGED BY
HPCG
01
167
c

SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO
.HPCG
01
16B
c

SUBROUTINE HPCG IS TERMINATEO.
HPCG
01
169
c
AUX
AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM
HPCG
01
170
c

COLUMNS.
HPCG
01
171
c


Hf?CG
01
172
c
REMARKS

HPCG
01
173
c
THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM. IF
HPCG
0 1
174
c
(1) MORE
THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE
HPCG
01
175
c
NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE
HPCG
m
176
c
IHLF
= 11) .
HPCG
01
177
c
(21 INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN
HPCG
01
178
c
(ERROR MESSAGES THLF=12 OR IHLF=13),
HPCG
01
179
c
(31 THE
WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,
HPCG
01
180
c
(4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.
HPCG
0]
181
c


HPCG
01
182
c
SUBROUTINES
AND FUNCTION SUHPROGRAMS REQUIRED
HPCG
01
183
c
THE EXTERNAL SUBROUTINES FCT(X,Y,DERY ANO
HPCG
01
184
c
OUTP(X.Y
.DERY.IHLF.NDIM.PRMT) MUST HE FURNISHED BY THE USER
.HPCG
01
185
c


HPCG
79

-------
U] 18b
C

METHOD
HPCG
01 107
C

EVALUATION IS DONE RY MEANS OF H AMM IN6S MODIFIED PREDICTOR-
HPCG
01181
C

CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD* USING 4
HPCG
01 19
C

PRECEEDI N(> POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE
HPCG
01 190
C

DEPENDENT VARIABLES.
HPCG
011^1
C

FOURTH ORDER RuNGE-KUTTA METHOD SUGGESTED BY RALSTON IS
HPCG
Ul 192
C

USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR
HPCG
01 193
C

COMPUTATION OF STARTING VALUES.
HPCG
01 194
C

SUHROUTINE HPCG AUTOMATICALLY.ADJUSTS THE INCREMENT DURING
HPCG
01 19b
C

THE WHOLE COMPUTATION BY HALVING OR DOUBLING.
HPCG
01 196
C

TO GET FULL FLEXIBILITY IN OUTPUT. AN OUTPUT SUBROUTINE
HPCG
01 197
C

MUST HE CODEO RY THE USER.
HPCG
01 19B
C

FOR REFERENCE. SEF
HPCG
01 199
C

(1) RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL
HPCG
01 200
C

COMPUTERS, WILEY. NFW YORK/LONDON, I960. PP.95-109.
HPCG
01201
C

(2) RALSTON, OUNGE .Y < 1 )f)ERY(1>,AUX(16.1)
HPCG
01209
c


HPCG
01210
c


HPCG
01211


N=1
HPCG
01212


THI.F = 0
HPCG
01213


X=PRMT(1)
HPCG
01214


h=prm r13)
HPCG
01215


DPMT(5)=0.
HPCG
0121b


no 1 I=1.NUIM.
HPCG
0 1217


AUX(16.I)=0.
HPCG
01 218


All* (15,1) =OERY ( I )
HPCG
01219

1
AUX(1I)=Y(I)
HPCG
0 1220


IF
-------
01252


N=1
HPCG
0 1 253


T5W = 2
HPCG
01254


r,oro loo
HPCG
01255
C


HPCG
01256

1.3
X = XH
HPCG
01257


CALL FCT (X.Y.DEKY)
HPCG
01 2S8


N = 2
HPCG
01259


no 14 I = 1t NO IM
HPCG
0 J 260


MIX (2. 1 ) =Y ( I J
HPCG
01 261

14
MJX(9I)=DEY(1)
HPCG
"1262


TSW = 3
HPCG
01263


noro loo
HPCG
01264
C


HPCG
01265
C

COMPUTATION OF TEST VALUE DELT
HPCG
01266

15
OFL T = 0.
HPCG
<> 1 26?


no 16 i=i.NnjM
HPCG
01268

16
DFL T =DELT~ AUX(15.1)ABS1Y(I)-AUX(4. I))
HPCG
01 269


DFL T = 06666667DEL T
HPCG
01270


TFH(.375AUX (8.I)~.791666 7AUX (9.I)
HPCG
0 1 295


I-.208 3333"AIJX(10*1)~.041666670EPY(I))
HPCG
01296

23
X = X.H
HPCG
01297


N=N*l
HPCG
01298


CALL FCT(X.Y.DERY)
HPCG
01299


CALL OUTP(X.Y.DERY.IHLF.NDIM.PRMT)
HPCG
01300


T F1PRMT(b) J6.24.6
HPCG
01301

24
TFtN-4)25.200.200
HPCG
01302

25
no 26 1=1.NOIM
HPCG
01303


AUX(N.I)=Y(I)
HPCG
0 1 304

26
4UX(N*7.I)=DERY(I)
HPCG
01 305


IF(N-3)27.29.20U
HPCG
01306
C


HPCG
01307

27
no 28 1=1.NDIM
HPCG
01308


nFLT = AUX(9.1)~ AUX(9.1)
HPCG
01309


dflt=oelt*uelt
HPCG
01310

28
Y(l)= AUX(1.I>~.J333333H(AUX(8.I>~DELT*AUX(10.I))
HPCG
01311


ROTO 23
HPCG
01312
C


HPCG
01313

29
no 30 1=1.NDIM
HPCG
01314


nELT = AUH(9.1)~AUX(10.I).
HPCG
01315


nFLT = DELT*DELT ~L)ELT
HPCG
0131b

30
Y(I)=AUX(1.I).375H0(AUX(8.I)OELT.AUX(1 1 . I))
HPCG
01317


ROTO 23
HPCG
81

-------
0
31 ft
C


HPCG
0
314
C

THp FOLLOWING PART OF SUHROUTINF HPCG COMMUTES HY MEANS OF
HPC6
0
3?0
C

RUNGE-KUTTA METHOD STARTING VALUES FOR THE" NOT SELK-START 1NG
HPCG
0
321
C

PREDICTOR-CORRECTOR METHOD.
hpcg
0
322

100
no 101 I = l.Nf)lM .
hpcg
0
323


7=HAUX =Z
hpcg
u
325

101
Y(I)= AUXIN.T)~.4Z
hpcg
0
326
C

7 IS AN AUXILIARY STORAGE LOCATION
hpcg
n
327
c


hpcg
0
328


7 = X.4 H
hpcg
u
329


CALL FCT (Z . Y.DEHY)
hpcg
0
330


no 102 i = ino Ih
hpcg
0
331


7=HOERY(l>
hpcg
0
332


AUX (6*1) -Z
hpcg
0
333

102
Y ( 1 ) =AUX (N. I) ~ .2969776AUX (5, 1 ) ~ . I587596Z
HPCG
0
334
c


HPCG
a
335


7=x ~ .4557372H
HPCG
a
336


CALL FCT(Z.Y.OERY)
HPCG
0
337


no 103 I = 1NDIM
HPCG
0
33fl


7=HDERY(I)
HPCG
0
339


AUX(7 * I ) =Z
HPCG
0
340

103
Y(I)= AUX(N.I)*.2181004AUX(5,I)-3.050965"AUX(6.1)*3.a32Hb5Z
HPCG
a
34 1
c


HPCG
0
34 2


7 = X H
HPCG
i)
34 3


CALL FCT(Z.Y.DERY)
HPCG
a
344


no 104 I=l,NOIM
HPCG
0
345

1040Y ( I ) = AUX (N. I ) ~. 1747603AUX IS. I ) -.S51 4807AUX 16. I )
HPCG
0
34 6


1 ~1 .205536AUX< 7.1)~.171184flHnKRY(I>
HPCG
0
347


r,010(9. 13. 15.21 ) .ISW
HPCG
0
34R
c


HPCG
0
349
c

POSSIBLE HREAK-POINT FOR LINKAGE
HPCG
0
350
c


HPCG
0
351
c

STARTING VALUES ARE COMPUTED.
HPCG
0
35?
c

NOW START HAMMINGS MODIFIED PEDICTOR-COHRfCTOR METHOD.
HPCG
0
353

200
TSTEP=3
HPCG
0
354

201
T F (N-8) 204. 202.204
HPCG
0
355
c


HPCG
0
356
c

n= causes the mows of aux to change their storage locations
HPCG
0
357

20?
no 203 N=2.7
HPCG
6
35B


no 203 1=1.ND1M
HPCG
0
359


AUX (N-l . I ) =AUX =AUX(N*7i I)
HPCG
0
361


N = 7
HPCG
0
362
c


HPCG
0
363
c

N LESS THAN 8 CAUSES N*1 TO (SET N
HPCG
0
364

204
N=N*1
HPCG
0
365
c


HPCG
0
366
c

COMPUTATION OF NEXT VECTOR Y
HPCG
0
36 7


DO 205 [=1.NDIM
HPCG
0
368


AUX(N-l.I)= Y(I)
HPCG
0
369

205
aux in6. i i =r>ERy < I >
HPCG
0
370


X = X + H
HPCG
0
371

206
TSTEP=ISTEP+1
HPCG
0
372


DO 20 7 I=1.NDIM
HPCG
0
373


OnFL T = AUX(N-4.I) ~1.3J3333*H6(AUXIN*6.1)*AIH IN*6.I)-AUXIN*5.1)~
HPCG
0
374


1AUX(N*4,It*AUX(N*4.I 1 )
HPCG
0
3 75


Y(I)=DELT-.925619flAUX(16,If
HPCG
0
376

207
AUX(16.1)=OELT
HPCG
0
377
c

PREDICTOR IS NOni GENERATED IN WOW 16 OK AUX. MODIFIED PREDICTOR
HPCG
0
37M
c

IS GENERATED IN Y. DFLT MEANS AN AUXILIARY STORAGE.
HPCG
0
379
c


HPCG
0
3*0


CALL FCT(X.Y.DERY)
HPCG
0
3fll
c

DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY
HPCG
0
3fl2
c


HPCG
0
3H3


110 208 1 = 1. ND IM
HPCG
82

-------
U1 384	ODFL I = . 12b(9.AUX (N-l,l)-AUX*J.H(UEkY*	HPCG
01305	1AUX(N*6.I)-AUX	HPCG
01388	C	HPCG
01389	C	TEST WHETHEP H MUST RE HALVED OR DOUBLED	HPCG
01390	OFLT=0.	HPCG
01391	no 209 1 = 1.NDIM	HPCG
01392	209 DFLT=DELT*AUX(15.1)ARS(AUX(16,I>)	HPCG
01393	TF(0LT-PRMT<4)1210.222.222	HPCG
01 39<	C	HPCG
01395	C H MUST NOT BE HALVED. THAT MEANS Yd) AWE GOOD.	HPCG
0J396	210 CALL FCT(X.Y.DEHY)	HPCG
01397	CALL OUTP(X.Y.DERY.IHLF.NDIM.PRMT)	HPCG
01398	TF(PRMT(5)>212.211.212	HPCG
01399	211 IF(IHLF-11)213.212.212	HPCG
01400	212 RETURN	HPCG
01401	213 TFIH(X-PRMT(2)))214.212.212	HPCG
01402	214 IF(ABS(X-PRMT(2))-.1ARS(H>)212,215.215	HPCG
01403	215 IF(DELT-.02PHMT(4))216.216.201	HPCG
01404	C	HPCG
01405	C	HPCG
01406	C	H COULD HE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE	HPCG
01407	C	AVAILABLE	HPCG
01408	216 TF(IHLF)201201217	HPCG
01409	217 fF=AUX(N-2, I >	HPCG
01418	AUX(N 2 ,1) = AUX(N4. I )	HPCG
01419	Al)X(N-3.I)=AUX(N-6.I)	HPCG
01420	AUX(N*6.I)=AUX(N*5.I)	HPCG
01421	AUX(N*5.I)=AUX(N*3,1 )	HPCG
01422	AtJX (N* 4  I ) = AUX (N* 1 . I )	HPCG
01423	nFLT=AUX(N*b.I)*AUX(N*5.I)	HPCG
01424	DFLT=OELTDELT*OELT	HPCG
01425	22 inAUX(16. I)=8.962963*(Y(1)-AUX(N-3.I)>-3.361111H(DERY(I) ~DELT	HPCG
01426	1*AUX(N*4.1)>	HPCG
01427	ROTO'201	HPCG
01428	C	HPCG
01429	C	HPCG
01430	C	H MUST 8E HALVED	HPCG
01431	222 IHLF = IHLF~ 1	HPCG
01432	IF(IHLF-10)223.223.210	HPCG
01433	223 H=.5H	HPCG
01434	TSTEP=0	HPCG
01435	DO 224 I = 1.NO IM	HPCG
01436	OY(I)=.00390625(80.AUX(N-l.I)~135.AUX(N-2.I)*40.AUX(N-3. I) ~ HPCG
01437	1AUX(N-4.I)>-. 1 171875(AUX(N*b.I)-6.AUX(N*S.I)-AUX(N+4,I))H	HPCG
01438	r>AUX (N4. I) =.00390625 (12.AUX (N-l. 1 ) *135.MUX (N-2. I) ~ ,	HPCG
01439	1108.AUX(N-3.I)~AUX(N-4,I) >-.0234375(AUX ~18.AUX(N*5.I>- HPCG
01440	29. AUX(N* 4.I))"H	HPCG
01441	AUX(N-3.I)=AUX(N2.1)	HPCG
01442	224 AUX(N4.I)=AUX(N.b.I)	HPCG
01443	X=X-H	HPCG
01444	OELT=X-(H*H)	HPCG
01445	CALL FCT(DELT.Y.DERY)	HPCG
01446	HO 225 I = 1NDIM	HPCG
01447	AUX(N2.I)=Y(I)	HPCG
01448	AUX(N*5.I)=DEWY(I>	HPCG
01449	225 Y(I>=AUX(N-4,I)	HPCG
83

-------
014b0 DF-'L r =UtL I - (rlH I	nrtu
01451 CALL FC T (DEL T  Y DEH Y)	HPCG
0145? 00 226 I = 1NO IM	HPCG
01453	DFLT=AUX	HPCG
01454	r>FLT = OELT*DELT*OtLT	HPCG
01 455	niUJX ( 16. I > = 0.962*63 UUX (N- 1 ,	I) -Y ( I) ) -3. 361 111H ( AUX 	HPCG
01457	226 AIJX (N j, 1 ) =DEKY ( I )	HPCG
01458	r,oT0 206	HPCG
0)45
-------
REFERENCES
1.	U.S. Environmental Protection Agency. April 25, 1978. Modification of
Secondary Treatment Requirements for Discharges into Marine Waters.
Federal Register, Vol. 43, No. 80, p. 17484-17508. Washington, D.C.
2.	33 USC 1251, et sec). PL 92-500 as amended December 28, 1977.
3.	California State Water Resources Control Board. 1978. Water Quality
Control Plan for Ocean Waters of California. Sacramento, California.
4.	EPA. 1976. Quality Criteria for Water. U.S. Environmental Protection
Agency, Washington, D.C. EPA-440/9-76-023.
5.	Baumgartner, D. J. and D. S. Trent. 1970. Ocean Outfall Design: Part
I. Literature Review and Theoretical Development. U.S. Dept. of Int.,
Federal Water Quality Adm. (NTIS No. PB-203 749).
6.	Briggs, G. A. 1969. Plume Rise. U.S. Atomic Energy Commission, Oak
Ridge Tennessee. (NTIS No. TID-25075).
7.	Davis, L. R. and M. A. Shirazi. 1978. A Review of Thermal Plume Model-
ing. Keynote address. In: Proceedings of the Sixth International Heat
Transfer Conf., ASME, Aug. 6-11, 1978, Toronto, Canada.
8.	Rouse, H., C. S. Yih and H. W. Humphreys. 1952. Gravitational Convec-
tion from a Boundary Source. Tellus, Vol. 4, pp. 201-210.
9.	Priestley, C. H:. B. and F. K. Ball. 1955. Continuous Convection from an
Isolated Source of Heat. Quarterly Journal of the Royal Meteor. Soc.,
Vol., 81,"p. 144-157.
85

-------
10.	Morton, B. R., G. I. Taylor and J. S. Turner. 1956. Turbulent Gravita-
tional Convection from Maintained and Instantaneous Sources. Proc. of
the Royal Soc. of London, Vol. A234, pp. 1-23.
11.	Morton, B. R. 1959. Forced Plumes. Journal of Fluid Mechanics, Vol. 5,
pp. 151-163.
12.	Abraham, G. 1963. Jet Diffusion in Stagnant Ambient Fluid. Delft Hy-
draulics Publication No. 29. Delft, Netherlands.
13.	Fan, L. H. 1967. Turbulent Buoyant Jets into Stratified and Flowing
Ambient Fluids. Cal. Inst. Tech., Keck Hydraulics Lab., Rep. No. KH-R-15.
Pasadena, California.
14.	Abraham, G. 1971. The Flow of Round Buoyant Jets Issuing Vertically
into Ambient Fluid Flowing in a Horizontal Direction. Delft Hydraulics
Publication No. 81, Delft, Netherlands.
15.	Cederwall, K. 1971. Buoyant Slot Jets into Stagnant or Flowing Environ-
ment. Cal, Inst, of Techn., Keck Lab. Rep. No. KH-R-25. Pasadena,
California.
16.	Roberts, P. J. W. 1977. Dispersion of Buoyant Waste Water Discharged
from Outfall Diffusers of Finite Length. Cal. Inst, of Techn., Keck Lab.
Report No. KH-R-35. Pasadena, California.
17.	Sotil, C. A. 1971. Computer Program for Slot Buoyant Jets into Strati-
fied Ambient Environments. Cal. Inst, of Techn., Keck Lab. Tech. Memo
71-2. Pasadena, California.
18.	Koh, R. C. and L. N. Fan. 1970. Mathematical Models for the Prediction
of Temperature Distribution Resulting from the Discharge of Heated Water
in Large Bodies of Water. U.S. Environmental Protection Agency, Water
Poll. Cont. Res. Series Rep. 1613 0DW0/70.

-------
19.	Baumgartner, D. J., D. S. Trent, and K. V. Byram. 1971. User's Guide
and Documentation for Outfall Plume Model. U.S. Environmental Protection
Agency, Pac. N.W. Water Lab. Working Paper No. 80. Corvallis, Oregon.
(NTIS No. PB 204-577/BA).
20.;	Winiarski, L. D. and W. E. Frick. 1976. Cooling Tower Plume Model.
U.S. Environmental Protection Agency, Corvallis Environ. Res. Lab.
EPA-600/3-76-100. Corvallis, Oregon.
21\ Davis, L. R. 1975. Analysis of Multiple Cell Mechanical Draft Cooling
Towers. U.S. Environmental Protection Agency, Corvallis Environ. Res.
Lab. EPA-660/3-75-039. Corvallis, Oregon.
22. Brooks, N. H. 1973. Dispersion in Hydrologic and Coastal Environments.
U.S. Environmental Protection Agency, Corvallis Environ. Res. Lab.
EPA-660/3-73-010. Corvallis, Oregon.
23: Winiarski, L. D. and W. E. Frick. 1978. Methods of improving Plume
Models. Presentation at Cooling Tower Environ.1978, Univ. of Maryland,
College Park, Maryland.
24! Kannberg, L. D. and L. R. Davis. 1976. An Experimental/Analytical In-
vestigation of Deep Submerged Multiple Buoyant Jets. U.S. Environmental
Protection Agency, Corvallis Environ. Res. Lab. EPA-600/3-76-001. Cor-
vallis, Oregon.
25.	Hirst, E. A. 1971. Analysis of Round, Turbulent, Buoyant Jets Dis-
charged into Flowing Stratified Ambients. U.S. Atomic Energy Commission,
Oak Ridge Nat. Lab., Rep. ORNL-4685. Oak Ridge, Tennessee.
26.	Hirst, E. A. 1971. Analysis of Buoyant Jets within the Zone of Flow
Establishment. U.S. Atomic Energy Commission, Oak Ridge Nat. Lab., Rep.
N. 0RNL-TM-3470. Oak Ridge, Tennessee.

-------
Guthrie, D. L. 1975. Documentation for OUTFALL: A Computer Program for
the Calculation of Outfall Length Based on Dilution Requirement. USEPA,
Region II. New York, New York.
88

-------