Mathematical Model of the Columbia River from the Pacific Ocean to Bonneville Dam Part I


FEDERAL WATER POLLUTION CONTROL ADMINISTRATION
i
NORTHWEST REGION. PACIFIC NORTHWEST WATER LABORATORY
MATHEMATICAL MODEL
OF THE COLUMBIA RIVER
FROM THE PACIFIC OCEAN
TO BONNEVILLE DAM
PART 1
CLEA
November 1969
DlSAlVOniTMEIT
r-	NOTE A	[
Columbia River entrance
Tb« proi«d depth it 48 loot Controlling depth! aia
publuhad monthly in th* lor.il Nolle* Ic Mannari by
Ihe (J S Colli Guild and monthly ill tha Mydrogriphit
011k. Notua lo M4>ir>«>» Additional information m«,
h» »bt*in«d from ih» Corpi ol Cn|in«an. IJ S. Army,
Portland Oregon
tWirklon Ml
UNITED STATES-WEST COAST
OREGON-WASHINGTON
COLUMBIA RIVER

-------
FEDERAL WATER POLLUTION CONTROL ADMINISTRATION
NORTHWEST REGION, PORTLAND, OREGON
James L. Agee, Regional Director
PACIFIC NORTHWEST WATER LABORATORY
CORVALLIS, OREGON
A. F. Bartsch, Director
NATIONAL THERMAL
POLLUTION RESEARCH
Frank H. Rainwater
NATIONAL COASTAL
POLLUTION RESEARCH
D. J. Baumgartner
BIOLOGICAL EFFECTS
Gerald R. Bouck
MANPOWER AND TRAINING
Lyman J. Nielson
NATIONAL EUTROPHICATION
RESEARCH
A. F. Bartsch
WASTE TREATMENT RESEARCH
AND TECHNOLOGY: Pulp &
Paper; Food Processing;
Wood Products & Logging;
Special Studies
James R. Boydston
CONSOLIDATED LABORATORY
SERVICES
Daniel F. Krawczyk
NATIONAL COASTAL POLLUTION
RESEARCH PROGRAM
D. J. Baumgartner, Chief
R. J. Callaway
M. H. Feldman
B. D. Clark
G. R. Ditsworth
W. A. DeBen
L. C. Bentsen
D. S. Trent
D. L. Cutchin
E. M. Gruchalla

-------
MATHEMATICAL MODEL OF THE COLUMBIA
RIVER FROM THE PACIFIC OCEAN TO
BONNEVILLE DAM
PART I
Theory, Program Notes and Programs
by
R. J. Callaway
K. V. Byram
G. R. Dltsworth
United States Department of, the Interior
Federal Water Pollution Control Administration, Northwest Region
Pacific Northwest Water Laboratory
200 South Th1rty-f1fth Street
Corvallis, Oregon 97330

-------
CONTENTS
Em
INTRODUCTION		1
Models of Pollution Problems 		5
Other (Large Scale) Models 		9
General Model Features 	 .........	13
General Flow Diagram 		16
MATHEMATICAL METHODS		21
Differential Equations - Terminology
and Assumptions			21
Finite Differences and Explicit Solutions. .....	27
The Explicit Solution of 3u/3t = 32u/ax2. ...	31
Runge-Kutta Solution of Hydraulic Equations		33
Diffusion and Dispersion			37
Stability, Numerical Mixing and Other Errors ....	41
Hydraulic Equations 		42
Dispersion Equation 		44
HEAT BUDGET TERMS		47
Heat Flux Through the Water Surface		49
Temperature Dependent Terms - Computation		50
Back Radiation				50
Evaporation Heat Exchange. . 			51
Convection Heat Exchange 		52
Summary of Heat Budget Step		53
Equilibrium Temperature		53
SOMATIZATION		55
General		55
Base Charts and Data Sources		55
Datum Planes		60
River Boundaries		60
Junction Point Layout 		61
Criteria for Selecting Junction Points		61

-------
CONTENTS (CONT.)
2m
Junction Data		63
Junction Numbers 		63
Junction Surface Area		64
Initial Head		64
Number of Channels at a Junction 		64
Junction Depths		65
Tributary Stream Flows 		68
Channel Data		69
Channel Numbers		69
Channel Length 		69
Channel Depths 			69
Channel Widths 		70
Cross-sectional Area		70
Channel Flow		70
Channel Velocity 		72
DISCUSSION 			73
ACKNOWLEDGEMENTS 		77
REFERENCES		79
APPENDIX I			83
Program and Notes for HYDRA 		83
APPENDIX II		103
Program and Notes for HYDEX 	 . 		103
APPENDIX III		117

-------
LIST OF FIGURES
Figure	Page
1	Columbia River - Oregon and Washington -
Pacific Ocean to Bonneville Dam 	 2
2	Schematic of Junction-Channel Network and
Possible Transfer Processes 	 14
3	General Flow Diagram: Constituent
Concentration and Temperature 	 17
4	Schematic of Currents, Columbia River 		26
5	Definition Sketch for Difference Equations		28
6	Explicit Integration Scheme 		31
7	Energy Balance Terms		48
8	Columbia River Schematization - River Mile
0.0 (Pacific Ocean) to River Mile 28.3
(Jim Crow Point)	 56
9	Columbia River Schematization - River Mile
28.3 (Jim Crow Point) to River Mile 87.0
(Lewis River)		 57
10	Columbia River Schematization - River Mile
87.0 (Lewis River) to River Mile 146.1
(Bonneville Dam)			 58
11	Columbia River Oregon-Washington
Schematization - River Mile 28 to River
Mile 35 (Exploded View of Figure 9 Top)	 59
12	Integer Terms Employed in Scheduling Print

-------
LIST OF TABLES
Table Page

-------
ABSTRACT
The Columbia River from the Pacific Ocean to Bonneville Dam
is treated as a series of two-dimensional finite elements in the
formulation of a mathematical model of the system.
Currents and stages are simulated along the river via an
explicit solution of the one-dimensional equations of motion and
continuity; two-dimensional conditions in the horizontal are
approached by means of a branched network of connecting channels and.
i
junctions. Computed net velocities and stages are used as inout to
the advection-diffusion equation and solutions are obtained for
any coupled (e.g., BOD-DO) or uncoupled, first order reaction,
conservative and/or non-conservative substance.
Emphasis is placed on obtaining a solution for temperature as
the dependent variable. Allowance is made for input of meteorological
variables and a stepwise heat budget computation is made in order
to predict temperature conditions on an hourly basis.
A discussion of some existing pollution models, numerical
methods and error sources is given; computer programs and program

-------
INTRODUCTION
This report is one result of a 1968 FWPCA decision to model
the Columbia River system from the Canadian border to the Pacific
Ocean for the purpose of evaluating existing and/or potential
thermal pollution problems. Described here are the mathematical
procedures» elementary theory, and documentation of computer pro-
grams employed in the lower Columbia study.
Part II of this report describes input procedures, provides
a test program and gives examples of actual output. Verification
procedures will also be given.
This work considers that portion of the Columbia from the
Pacific Ocean to Bonneville Dam (Figure 1). The system above
Bonneville has been treated as comprised of unstratified reservoirs
(Morse, 1969) and stratified reservoirs (WRE, 1969).
At low flow, tidal effects in the form of a small diurnal
tidal rise and fall are observable at the dam; by some definitions,
the system up to the dam could be considered an estuary. However,
the estuarine portion is usually restricted to that semi-enclosed
part of the lower river where salt water is present. The fresh-
water portion of the river, where ocean generated tidal effects
occur, is called the tidal river.
In order to model the entire 146 miles of the Columbia to the

-------
Jim Crow Pt.
A»toriq
u
OREGON
AREA SHOWN
Lonqview
Cowlitz R.
WASHINGTON
Lewis R.
yoncouver
Washougal R.
i
Bonneville Dam
Portland
:si
Willamette:
i;R.
, Sandy R.
0	8 10 20 30
STATUTE MILES
FIGURE I. COLUMBIA RIVER - PACIFIC OCEAN TO

-------
3
the system includes estuarine to river-run circulation patterns,
it was felt that a time record of longitudinal tidal flows and
stage elevations would be required. Because there are many islands
and tributaries on the main stream, it was also felt that these
features should be incorporated in the model.
One method of solving this problem would have been to start
from scratch and develop an in-house model of the system. However,
the existence of a rather ingenious model developed by Water
Resource Engineers of California (WRE) for the San Francisco Bay-
Delta system tended to discourage this approach especially as the
model had proved to be quite versatile in handling a number of
situations. Foremost in these considerations was the ability of
the model to approach two-dimensional (horizontally) conditions;
not the least was the fact that several years of running experience
were built into it.
Accordingly, copies of the program and decks were obtained
from the Southwest Region of the FWPCA - the original contracting
agency* - through the courtesy of Dr. Howard Harris and Mr. Ken
Feigner**, who explained the basic workings of the programs, made
suggestions on schematization and provided us with many helpful
suggestions and comments.
*Actually, Public Health Service before the creation of the FWPCA.
**Note added 1n press. Mr. Feigner 1s currently completing a
documentation of FWPCA experience with the model 1n San Francisco
and San Diego waters. It will serve as a valuable accompaniment

-------
4
The computer program developed by WRE did not treat
temperature as such, hence, it was modified by us to accept
temperature as a variable. Because of the general nature of the
water quality portion of the program, i.e., its ability to handle
conservative and non-conservative substances and to couple them if
required, it was decided to retain those features rather than strip
the model to handle only temperature. Thus, the description that
follows emphasizes the methods employed in the temperature compu-
tations, but not at the exclusion of DO or BOD or any other sub-
stance.
It should be borne in mind that this is in actuality a one-
dimensional model although provision is made to branch flows at
junctions. Any substance discharged into a junction, is by the
one-dimensional assumption, assumed to completely mix throughout
the junction at each time step before being advected or diffused
to another junction through a channel. Any numerical model will
have similar artificialities; unfortunately, it is usually left to
the user to uncover them for himself. Because the program
received was undocumented, it was felt necessary to document it
for those who might want to employ it for production runs. Rather
than give only a description of card input requirements, a full
documentation was developed because of the rather extensive knowl-

-------
5
of the subjects involved: open channel flow, diffusion and
dispersion processes, numerical methods, sanitary engineering and
heat budget methods in addition to estuarine flow processes. If
the user is to be other than a knob turner, he should develop capa-
bility in these fields. If the program notes and literature cited
are studied carefully, they should provide an independent start
for getting on the estuary bandwagon complete with thermal pollu-
tion weaponry.
Models of Pollution Problems
Based on the premise that many marine pollution problems can
be solved via computer methods, this section sets forth assump-
tions and limitations of some models currently in use. The "prob-
lem" is stipulated to be relatable to the physical environment,
i.e., whether or not a bad situation will result is predictable on
the basis of the pollutant's reaction rate and the hydrodynamic
situation in the effluent discharge area. The condition is thus
restricted to the prediction of the concentration of a specific
pollutant at a given time and place given certain information on
discharge rates, concentrations, and flow and diffusion 1n the
estuary. How these predicted concentrations will affect the biota
or whether or not they will lead to synergistic or antagonistic

-------
6
Deterministic (as opposed to stochastic) models of the
environment are either steady-state or time varying. The steady-
state assumption simply means that there is no concentration
change of a substance or property with time. The effluent is
discharged at a constant rate, and has been discharged for a long
enough period to come into equilibrium with the receiving waters;
any fresh water flow to the environment is constant, diffusion
rates and other characteristics are also steady. The topography
of the estuary can be modeled quite closely, i.e., any tide level,
cross-sectional areas can be incorporated to show the irregular
nature of the geographic setting. However, the effect of tidal
height variations on cross-sectional areas (hence, water volume
changes) and tidal current fluctuations cannot be modeled here
except by repeated application of the steady-state case, in which
case there would evolve a process of simulation. Simulation of
various reaction rates, river flows, diffusion and reaeration
rates is a logical extension of the steady-state assumption and
perhaps the best justification for its use. For, by simulation,
the expected range of concentration of a given pollutant can be
easily explored by use of a steady-state digital model. Input
information to a complex area can be obtained from existing hydro-
graphic charts, flows can usually be extracted from federal or

-------
7
developed steady-state model, as opposed to the judgement
necessary to carry out a realistic simulation, is elementary.
(Interpretation of results is, as always, the ultimate hangup;
however, this does not relate to the present discussion.) The
steady-state model, then, is useful in a situation where a rapid,
first-cut approximation to a situation will suffice. In a highly
complex industrialized setting such as the Delaware Estuary, the
steady-state case has been used as the foundation of a linear pro-
gramming method of meeting certain water quality standards. For
instance, if wastes of known volume, concentration, and reaction
rates are discharged at various locations along some miles of an
estuary and a dissolved oxygen standard of, say, 5 ppm is to be
obtained, the linear programming concept can be used in conjunction
with the steady-state case to ensure that this goal will be met
most of the time and at the least expense to the parties involved.
Various external constraints are, of course, involved here, but
the tools are available for the exercise of logical and unarbi-
trary decision making. Progress in extending this concept to the
dynamic situation is underway. It is safe to predict that tool-
maklng will precede the implementation of these devices. The
reason for this will be obvious to any manager who is or has been
involved with a decision that has crossed political boundaries

-------
8
While the steady-state model has its uses, it also has its
drawbacks. The fact is steady-state situations in nature don't
really exist; hence, the absolute verification of such a condition
is impossible. Most such problems have escape hatches; with the
environmental scientist or engineer, the size of the hatch opening
depends on how loose a definition of steady-state he is willing to
accept. The purist will not be satisfied that steady-state veri-
fication has, in fact, been accomplished; nagging doubts will
remain until he has gained: 1) experience with such models,
2) judgement on how critical a condition of, say, flow variation
with time really is, 3) the realization that one is not usually
concerned with precision in, e.g., the second decimal of the D.O.
concentration.
Thus far, mention has not been made of the dimensionality of
the problem. Here is meant the variation of water quality con-
ditions with depth across stream and along the axis of the stream.
The first stream model, proposed by Streeter and Phelps (1925),
dealt with a freshwater condition and no variation of density was
assumed with depth. Lateral (cross-stream) variations were also
neglected, hence, the only gradient in concentration allowed was
longitudinal (along the stream axis). Vertical variations in
density occur In fresh water bodies, but unless the stream is
deep, turbulent mixing ensures that such gradients are minimized.

-------
9
headwaters of a reservoir. The reservoir may be markedly stratified
during summer; use of a one-dimensional model obviously doesn't
make sense in such a case although it could be implemented to
grind out neat rows of numbers.
Proceeding from the freshwater to the seawater environment
also usually means leaving the quasi-one-dimensional state and
entering at least a periodically stratified water body. In the salt
water portion of an estuary, one-dimensionality has in the past
been inferred from a vertical profile of salinity showing little or
no variation. The steady-state velocity distribution was also
assumed to be invariant from top to bottom. Recent theoretical
investigations (Hansen and Rattray, 1965) have shown that the verti-
cal current profile need not be exactly related to the salinity
distribution, although one's intuition would probably argue other-
wise.
Other (Large Scale) Models
Presented here is a brief discussion of the basic philosophy
and assumptions underlying models such as used by Thomann,
O'Connor, and others on the East Coast and the modified Water
Resources Engineers model used here, Then a description of the
general flow diagram of the entire system is given in order that
the functional interrelationships of the different parts of the

-------
10
It is noted that the primary difference between the WRE
model and that of Thomann (1963) is that the former representation
of estuarine flow computes intratidal velocities, while the latter
doesn't. There is, then a difference in viewpoint on how big a
time average one is justified in taking. The original Thomann
model used a time average of one day (numerical step size is
smaller). One reason for this large time average is a matter of
philosophy, namely that pollution control measures (measures that
the model output indicates should be taken) on the order of a day
are feasible, but those on a scale of hours generally are not. A
recent paper by O'Connor, et al., (1968), indicates that the
"...flux due to the tidal velocity, however, is too complex to be
explicitly included in the mass balance." O'Connor's model
integrates from slack tide to slack tide "...when the tidal
velocity is zero."
One may argue that Thomann's original model took too large a
time bite; but it must be remembered that his verification period
consisted in simulating the dissolved oxygen profile at various
points in the Delaware for one year. Shorter time periods, on the
order of the WRE model, could have been included but the input-
output problems would have been horrendous, to put it mildly.
Accepting the idea that control measures in the Delaware need not

-------
11
been gained by reducing the time step significantly, rf it is
assumed that the short time hydraulic effects do not affect the
overall waste distribution computed. In any event, the time
average employed is quite an important consideration and must be
carefully spelled out.
The hydrography of the Columbia River is quite different from
the Delaware. Discharge at the mouth is some 40 to 100 times as
great; saltwater penetration is at most 25 miles upstream (Hansen,
1965), while in the Delaware, it is about three to four times that;
tides are mixed, etc. The Columbia contains many islands and
several channels may cut through small areas of the river. Tidal
current reversal in the Columbia occurs some 75 to 100 miles
upstream during low flow, although tide effects (vertical motion)
can be seen at Bonneville Dam (Mile 146).
The steady-state assumption is an attractive one if for no
other reason than that programming and computational effort necessary
to achieve it is slight compared to the transient cases.
The use of a one-dimensional model is another questionable
assumption, even though the model discussed is a "quasi" two-
dimensional system. Obviously, in a stream as large as the
Columbia, cross-channel velocity variations will be quite large;
simulating a point source outfall on one bank of the river and
then insisting that the effluent will be immediately and completely

-------
12
devout simulation enthusiast to swallow a bit much. This is
particularly true in the light of recent evidence that Taylor-
type mixing probably won't occur for some distance (large
diffusion time) downstream of a point source (Fischer, 1968). It
is also true that the downstream distribution of waste discharged
close to a bank in a large river system (width/depth ratio >>1)
will usually be constrained to remain near that bank for some
distance downstream. The utility of the model being used, then,
is not in the simulation of small scale events, but as an indicator
of the meteorological effects on a very large system. It is not
unlikely that a small-scale model of waste heat discharge to the
Columbia will encompass a single junction of the large model.
Recent work by Leendertse (1967) and W. Hansen (1966) on two-
dimensional modeling will certainly provide a major step forward
in solving pollution problems in embayments and coastal regions.
In treating the nonsteady dispersion equation, a great deal
of computational effort is devoted to computing velocities in a
(finite-grid) network of channels. If one were able to specify
the velocities functionally then the largest part of the problem
would be solved since the dispersion equation could be solved
directly. In addition, if the diffusion coefficients were known
functionally or could be assumed constant, another saving in
labor could be effected. Such is not the case, unfortunately, and

-------
13
continuity equations in such a fashion that 1) the tide wave
amplitude and phase are verified with distance from the input wave
and 2) flows and direction of flows are in reasonable agreement
with known input lateral and mainstream flows. The constitution
of a "reasonable" agreement between observed and computed flow is
not easy to discern since there will always be discrepancies in,
among other things, input conditions assumed and those actually
occurring. This is, then, a problem of verification, which is
discussed in Part II.
General Model Features
Several widely scattered papers have been published on the
water quality aspects of the WRE model, for instance, Shubinski
et al. (1965), and Orlob et al. (1967). A recent paper by Orlob
(1968) discusses the various processes involved in modifying con-
centrations, particular their relation to the model's channels and
junctions (or nodes).
In the model, physical characteristics of a real setting
(Figure 2) are represented by junctions which occur at physical
branches or at somewhat arbitrary spacings between branches in a
network of channels and junctions.
Junctions connect short straight segments of regular cross-
section; these segments are termed channels. Inflow-outflow

-------
NODE OR
JUNCTION
Short a Long
Wave Radiation
Sensible Heat
CHANNEL
Growth
(e.g. Bacteria) Import
Surface
Mass Transfer
sv*:
Back
Radiation
Evaporation
Sensible Heat
Transfer
Decay Export
r • « i
MAIN
STREAM

TRIBUTARY
FIGURE 2. SCHEMATIC OF JUNCTION-CHANNEL NETWORK
AND

-------
surface area, and head; in addition, constituent mass, decay, and
growth rates are junction properties.
A channel is characterized by length, width, cross-sectional
area, and hydraulic radius; in addition, net flows, velocities
and friction are channel properties.
In essence, storage is provided at the junctions as well as
potential and input-output; the channels provide conveyance
between junctions.
Figure 2 summarizes the processes occuring in a schematic
junction of a channel network.
At Junction 2, net change in heat or mass, AM, during a time
step is brought about by the following:
AM = Advection ± Diffusion ± Heat transfer process
+ Surface mass transfer + Growth + Import
- Evaporation - Decay - Export.
Of course, change may only occur by processes of advection and
diffusion or in combination with the remaining terms. If only
temperature is being considered then only the first three terms
on the right are used, since evaporation is computed separately.
The above is an expression of the advection-diffuslon
equation with source and sink terms and is solved numerically

-------
16
in the advectiori term. If the solution is in terms of temperature
(as a constituent) then it is an expression of the energy equation.
General Flow Diagram
In summary, the programs solve for current velocity and tide
stage in one program; net velocities and heads are averaged over a
suitable time period in a subroutine which is used as input to a
stepwise solution of the dispersion equation.
Referring to Figure 3, the following step-by-step description
of the general computer program can be used to define what is
happening from the initial step of obtaining depth information to
the final one of printing out predicted temperature or concentra-
tions. As indicated above, most of the work occurs in step 2

-------
^ Initial
Conditions
Junction and
Channel Data
Print Input
Data,
Y,V,Q _
Print Min/Max
and Ave. flows
and heads
f Store
Wet Flows
V Heads
Compute
Y,V,Q
HYDRA
Condense
HYDRA for
Input to
QUAL
HYDEX
FIGURE 3. General Flow Diagram -
Constituent Cone. & Temperatures
(See text)
Print Heads
vs Time
Scratch
Print Min/Max
Flows & Cone.
Print Meteor
Input Data
Input Meteor.
Terms.
METDTA
Interpol ate
Input Data
EXQUA
Condense
QUAL Output
QUAL
Compute
Concentration
Heat Budget

-------
18
General Flow Diagram
For Use with Figure 3
1.	Initial and boundary conditions, junction and channel
data, such as cross-sectional areas and channel lengths, are read
into Program HYDRA.
2.	HYDRA computes heads at each junction and velocities
and flows in each channel. These data are printed versus time
(3) and/or summarized in subroutine HYDEX (4). Program will
terminate after (3) if HYDEX is not called.
3.	Print routine can be scheduled to list all or portions
of the output.
4.	HYDEX averages the data over certain time intervals
(e.g., 15-30 minutes) for input to binary tape (6) and/or averages
heads, flows and velocities for tidal cycles, days, etc., for
printout.
5.	Printout from HYDEX; execution will terminate here if
QUAL is not used.
6.	Net flows computed in HYDEX are stored in binary for
use by QUAL.
7.	Program QUAL needs average net channel flows calculated
in HYDEX to run. From input initial and boundary conditions,
QUAL computes concentrations of substances released at any net-

-------
19
8.	Local climatological data (net radiation - computed Qr
observed, air temperature, cloud cover, wind speed, etc.) are read
into subroutine METDTA if temperature is a constituent (9).
9.	METDTA interpolates incoming radiation and other terms
to conform to the selected quality time step,
10.	Printout of meteorological data.
11.	Printout of temperature and/or concentrations (up to
five constituents are allowed) occurs here. Program may termi-
nate here or pass to 13.
12.	If subroutine EXQUA is called, data is stored on binary
for execution by EXQUA. (Subroutine EXQUA is not discussed 1n this
report but will be made available on request.)
13.	EXQUA can be reprogrammed to summarize data 1n a manner
similar to HYDEX (4).
14.	Printout of a computation using EXQUA would be the final

-------
MATHEMATICAL METHODS
Differential Equations - Terminology and Assumptions
The programs discussed present numerical solutions to one-
dimensional linear or nonlinear partial differential equations
that are coupled or uncoupled for substances that are conservative
or nonconservative. The foregoing jargon is helpful in seeing
through the bramble bush of the leapfrog solutions and other
manipulations which are conceptually simple, but sometimes hard
to follow. When all is said and done, we are faced with solving
the "fundamental equation of linear sanitary engineering" which,
in operational form in one-dimension, is:

2	— (D —) + u —
at ax * L 3x 3x
(U C, T,...) = ES,
where BOD (L), D.O. (C), Temperature (T), etc., can be expressed
as a sum of sources and sinks (sS).
Expressing equation (1) in the simplest form for all three
variables:
2)	ik . D 92L . .. 3L u ,
3t " DL	u ax * -KiL
3)	3^C . |, 3C _ |> I . |> /£ - C)
sT °L 3X2-+ u w- K-jL + K2 (Cs I)
4)	II . n 32T x „ 3T » /T T)

-------
22
L	=	BOD concentration
C	=	DO concentration
Cs	= DO saturation concentration
Dl	=	Coefficient of long, dispersion
Kj	=	Decay rate
K2	=	Reaeration rate
K3	=	A thermal exchange rate
I = Water temperature
T = Equilibrium temperature
C
u « Mean velocity.
The equations differ only in the source and sink terms which
are peculiar to the particular substance. If there were no
reaction terms (ICj = K2 = Kg = 0), then the solution for one sub-
stance would be a simple multiple of another if, and only if,
the diffusion rate, D^, for each were equal and constant or varied
alike with distance. (Such a condition is known as the Reynolds
analogy, I.e., assuming that turbulent transfer rate of, say, heat
is the same as that of oxygen.)
It can be seen that (2) must be solved for L before (3) can
be computed (if the reaction rates are nonzero). The two equations
are thus said to be coupled through L. If the reaction rates are
nonzero, the substance (e.g., BOD) is said to be nonconservative;

-------
23
The one-dimensional assumption is inherent in equation (1)
as change is allowed only in the x (longitudinal) direction.
Equation (1) is the simplified form of the local time change. In
full bloom, the operator is written as:
where the y, z diffusion terms are allowed to vary. This equation
is merely a statement; it says nothing about what processes are
affecting the distribution (see Sverdrup, et al., Chapter V, 1942,
for an excellent discussion of the distribution of variables).
Sometimes the longitudinal and/or vertical terms are neglected
because the velocities involved are assumed to be very small. It
may be that that is so, but it can also be true that the y and z
gradients are very large so that the products (v |y, w |j) may not
be negligible. If the cross-stream and vertical advection terms
are to be neglected then they must either effectively cancel each
other or be very small. When calling on the one-dimensional
assumption, the gradients involved must be assumed to be negligible',
this is the condition that obtains when an estuary is "wel1-mixed."
The remaining terms to be discussed concern linearity. The
so-called non-linear terms, if not neglected, cause dreadfuls to

-------
is found for the system, then additional solutions can be obtained
by multiplying the answers {which might be the longitudinal BOD
concentrations) by any given number. This number might correspond
to, say, an increase or decrease in waste treatment. At any rate,
the solutions are said to be superposable. If the system is non-
linear, then multiplying by a number in one position will not
necessarily give a proportional output as the answer somewhere else.
As a result, many, many analytical solutions may be required to
determine the output in a nonlinear system, where a single solution
may suffice in a linear one.
In dealing with the hydraulic equation (for a complete dis-
cussion, see Dronkers (1964), Baltzer and Lai (1968), and
j	» 2
Leendertse (1967)), retention of the nonlinear term (u jT
is usually required since it may be at least equal in magnitude
to the linear terms. A tide wave becomes distorted with distance
upstream because of changes in channel configuration and roughness
through this nonlinear term and the nonlinear frictional term
(ku2). This is implied from the characteristic of linear systems
by noting that the output generated by a sinusoidal input 1s
also strictly sinusoidal even though the phase may be shifted
*Also called the convective-inertial term or the advection of

-------
25
and its amplitude modified. A problem in the prediction of water
height in estuaries and tidal rivers concerns the nonlinearity of
the system as the wave is distorted with its passage upstream. A
wave describable by a single harmonic (for a short period) at the
estuary mouth may require many harmonics for its description
further upstream.
There are some problems in the practical use of equations 1,
2 in estuaries. First of all, there are irregular boundaries; u
as used here is the net freshwater velocity (Q/A) and we really
should consider u = Q/A + utcos(wt), where ut is a tidal term and
w = 2tt/T» where T is a tidal period. In practice, the cosine term
can be replaced by a Fourier series to represent any degree of tide
complexity required. The one-dimensional pitfalls are fairly
obvious, but one should bear in mind that this implies a uniform
velocity from top to bottom and side to side (no shear). If there
isn't any shear, then the primary turbulence generating mechanism
is lost. We can overcome this (ignore 1t) by simply assigning a
certain value to the diffusion term. (D^ In this case.) The sur-
prising thing about all this, considering the assumptions, etc.,
is that with a finite difference model 1t can all be made to work,
I.e., serve as a pollution planning and management tool.
We'll need to know, or might want to calculate, the velocity

-------
26
This varies from strictly seaward directed river flow in the upper
reaches (with a bit of a sine wave thrown in) to a diurnally
reversing current in the estuary as shown in Figure 4 for the
Columbia.
Ebb
Ebb
Ebb
			
"O

Flood
F1 ood
F1 ood


ASTORIA
BEAVER
BONNEVILLE

FIGURE 4. Schematic of Currents, Columbia River
Current velocity is obtained by solving the equations of
motion (5) and continuity (6):
5)
6) 3h . 1 da _ n _ _ *
it * F ax _ °' q _ Au"
Suffice it to say that equations 5, 6 are written in finite
difference form, the estuary is schematized, i.e., depths,
areas (A), widths (b), roughness coefficients (k), are determined,
initial conditions are specified and u (and h) are solved for by
the "leapfrog" method which is employed in solving the coupled

-------
27
initial conditions of velocity and stage are read into the computer
along with the boundary conditions. Velocity and flow are computed
from the momentum equation; the computed flow is substituted into
the continuity equation to obtain a new stage elevation which is
then used in place of the initial condition to obtain new velocity
and flow values. The new flow obtained is again substituted into
the continuity equation and the process leapfrogs until the cycle
is complete.
Finite Differences and Explicit Solutions
In dealing with non-analytical solutions to differential
equations, it is necessary to express derivatives in a form
that the computer can handle, namely, finite difference approxi-
mations of infinitesimal quantities. What one really wants is to
make the infinitesimally small derivative as big (finite) as
possible while still satisfying the equation of motion or any
other equation.
The usual drill is to start with the Taylor series expansion
about x of a function, say u(x), which doesn't contain any sudden

-------
28
and
8)	+
Equation (7) could be used to predict the value of u(x) a dis-
tance ax ahead of it if the function and its derivatives were
known at x. How good the approximation is will depend on how
large h is.
Difference approximations are classified as either forward,
backward, or central. A particular computing scheme may make use
of one or more of these approximations, and the proper formula-
tion must be employed to ensure a stable and convergent solution.
Referring to Figure 5, it can be seen that the first deriva-
tive of u(x) centered about the point P can be written by
inspection as:
9) u'(x) =	= ^L. |u(x+ax) - u(x-ax)|, where the chord AB
is tangent to u(x) at P.
u(x)f
0
X-AX
X
X+AX

-------
29
The same expression can be obtained by subtracting equation {8)
from (7) and neglecting terms greater than or equal to ax3. The
error is then said to be of order 3, is written as 0(ax3), and is
the result of chopping off (truncating) the higher order terms.
Truncation error is inevitable, but can be made insignificant.
The second derivative of u(x) at P can be written by adding
equations (7) and (8) and neglecting terms of 0(ax^) and higher:
U(x+AX) + u(x-AX) ¦ 2u(x) + AX2	+ OfAX4)
or
a ju(x+Ax) + u(x-AX) - 2u(x)j.
The forward difference approximation of the slope at P
(^) 1s:
1°) u'(x) =-L |u(x+ax) - u(x)}-,
hence, values only at P and forward of 1t are used. Similarly,
the backward difference 1s:
11> U'(x) - jl{u(x) - U(X-4X)}.
Since nonsteady-state problems must be dealt with, provision
must be made to move the solutions ahead 1n time as well as along

-------
30
It is often desirable to solve a class of problems in such a
manner that recomputing needn't be done every time there is a
change in scale of a particular geometric or physical property, i.e.,
it shouldn't be required to compute the temperature distribution in
a rod for every length of rod imaginable. Such a process occurs
when the equations are expressed in terms of nondimensional
variables. For instance, the parabolic heat equation describing
the transient temperature distribution* in a rod can be written as:
12)	aU _ r 32U
sT a5P" '
where c is a constant; U is temperature; X is the distance from
one end of the (thin, uniform) rod; and T is time. By making
suitable transformations, this equation can be expressed in non-
dimensional form as:
13)	au „ a2u
at ~ ax2"
A finite difference grid can be used for the numerical solution
of this equation. The "explicit" method is Illustrated because
it is used herein to solve the hydraulic and dispersion equations.
Advantages and disadvantages of the explicit scheme are discussed
later.
*This is also one form of the F1ck diffusion equation which 1s

-------
31
The Explicit Solution of.|£ = |^T
The finite difference form of the nondimensional heat equation
is:
14) U-! -!X1 ' U-5 ¦} U-!x1 i	i + ui 1 1
t.J+i	= i+l .J 1»J	iJ
At	(AX)*
A forward difference is used for the time step and a central
difference is used for the second (space) derivative; the sub-
scripts are shown in Figure 6.
- j+1
¦ j

Ui J+1

I c
!+
L

A
s
s'
~
X '
X-'
X
X ^
ui-l,j
"1.J
ui+i .j
' u
1+2. j

1-1
¦+¦1 AX
n
h
H
FIGURE 6. Explicit Integration Scheme
Rearranging equation (14) for the value of u(x,t) after one time
step:

-------
32
The value of ^ (located at the j+1 row by an x) can next
be computed from the values u. u,x1 and u.., .. The whole
l ,J IT I , J	1+t , J
scheme can be repeated until values are known for row j+1; these
can then be used to obtain new (i.e., j+2) values at the next time
step. This explicit formulation requires that the initial and
boundary conditions be given.
Of critical importance in the numerical solution of the
parabolic heat equation is that the ratio of time step, At, to the
square of the space step, (Ax)2, must lie between 0 and 1/2. This
relates to the "stability" of the solution, a subject which will
be treated later on in the treatment of the wave equation which
has a somewhat different stability criterion. The use of a central
difference formulation can create problems; these are also discussed
in the section on stability.
It is possible to solve the system of equations simultaneously
by matrix inversion or some other "implicit" method which has the
advantage of being unconditionally stable for large time steps.
Even here, however, short time steps may be required to obtain the
necessary accuracy and to minimize numerical violations of water
mass and constituent concentration conservation. It is not clear
if an implicit solution for the Columbia would have justified the
considerable reprogramming effort that would have had to be

-------
33
Runge-Kutta Solution of Hydraulic Equations
Although any method of forward integration could be used, a
two (rather than the usual four) step Runge-Kutta (R-K) procedure
is employed in the solution of the equations of motion and con-
tinuity. Other methods are known to be more efficient but have
not yet been considered. The principal advantages in using R-K
methods lies in their independence of past computing stages, i.e.,
the method is self-starting. The R-K method is also stable when
grid spacings are uneven or change during computation. It is
difficult, however, to estimate the truncation error at a given
point in the computation although estimates can be obtained (see,
e.g., Hildebrand, p. 238, 1956).
For a channel with constant width and employing a slightly
different notation than before, the continuity equation 1s:
16) kWBM.O.
where
V a Average channel velocity during a time step (it)
A « Cross-sectional area of channel
H * Height of water surface above (arbitrary) horizontal
datum
B * Channel width.
The equation of motion In the x-direction 1s:

-------
34
where
g = Acceleration due to gravity
K = Friction coefficient
and the absolute value of V ensures the proper direction of the
frictional force, namely opposite to the direction of V.
The assumptions on which 16and 17are based are as follows:
1.	Acceleration and momentum transfer normal to the x-axis
is negligible. Thus, tributary inflows contribute to a change in
junction head, but impart no momentum during the contribution.
2.	Wave length is at least twice channel depth. If not,
the shallow water assumption utilized here (in which wave celerity
t/gfi") would not hold. The "wave" referred to here is of tidal
period, not a wind wave.
3.	Coriolis and wind forces are negligible.
4.	Each channel is straight (hence, no centrifugal effects)
and the cross-sectional area does not vary over its length.
The steps outlined below are contained in sequence numbers
147 - 207 in the listing for program HYDRA.
Following the notation of Shubinski and Sheffey (1966), and
Shubinski, et al. (1965, 1967) equations can be written for

-------
35
where
V. = ith channel velocity
At = time step (for R-K integration)
aV
•— = velocity gradient evaluated as suggested by Lai (1966)
n
as follows:
equation 16 is rewritten as
91 = _ B 9H _ V ciA
3X J 3t " J 3X *
expressed in finite difference form and substituted as
"n
K. = frictional resistance coefficient
g = gravitational acceleration
aH a head (potential) difference between junctions at ends
of channel
L.,. = channel length.
Similarly, the continuity equation is
19) aH. Q,
aF" 18 IT 9 w^ere ^ 1s now a Junction Index and
J
Qj = net flow at j during a time step, At
Aj ¦ junction surface area (constant)

-------
36
The solution of 18 and 19 using a two-step R-K (leapfrog)
procedure is as follows:
1. Initial and boundary conditions are specified so that
the system state is known at time t. Predictions are required at
time (t + At) and multiples thereof. Superscripts t, t+1/4, t+1/2,
t+1, imply values at time t, t+At, t+At, t+At, respectively. The
4 2
superscript t+1/4 indicates a term using mixed time steps.
2. Compute half-interval velocities and "quarter"-interval
channel flow	t
l1 /9	+ a+ +	U^)
wt+1/2 _ yt . At /wt A n |/t|wt|,,t _ AH*
vt - vi+ r (vi a^- Kilvilvi - s q-
gt+l/4 _ yt+l/2 At
1	1
3.	Compute half-interval heads and quarter-interval
channel areas	t
*Tm - Hj+ f <*j->
At+V4 = At + ^1. (AHt+l/2 _ AHtj
4.	Compute full-interval velocities and three-quarter
interval channel flow	t+1/2
V?+1 = v*+1/2 + At(V?+1^2
1	1	v 1 AXn
„t+l/2,vt+l/2,vt+l/2 „ 4Ht+1/21
i ivi ivi s rj- )
nt+3/4 _ wt+1 *t+1/2

-------
37
5.	Compute full-interval heads and three quarter interval
areas ,+1/2
H*+1 = hV*"1/2 + At(^- )
J	J	Aj
At+3/t At+l/2 + h (aHW . iHt*l/2)
6.	Upgrade system parameters, K, Q, A, which can be com-
puted from geometric considerations,etc.
7.	Continue at step 2 until cycle is complete.
Diffusion and Dispersion
The spreading out of material from a point source is easy to
visualize in terms of an instantaneous release, but real life
effluent discharges are more likely to be continuous or periodic.
A continuous release, however, can be synthesized analytically
from a sum of instantaneous releases so that a discussion of the
longitudinal diffusion of material properly starts with Instan-
taneous releases, (See Okubo and Karwelt, 1969, for a discussion
of the above as well as on the effect of shear on diffusion.)
Estuarine pollution models derive from the advectlon-
diffusion equation (as does a river model) which can be written
as:
20)
— + u — + v — + w — - f—(d —) + —(D —) + —(D —* eS

-------
38
By comparison with equation (1), it can be seen that the
cross-stream and vertical terms were neglected in getting to the
one-dimensional equation.
Confusion as to the meaning of dispersion as opposed to
diffusion is easily rectified if equation (20) is referred to as
the "dispersion equation." Dispersion will then include the
advective transport of material as well as its diffusion due to
turbulent flux.
If a one-dimensional coordinate system moves with the center
of mass of material, equation (20) degenerates into the Fick
equation originally developed to describe molecular scale phenomena
in which local concentration changes are due to diffusion only
(and diffusion is constant):
21)	3c . n 32c
9t UX 3X7" *
In large scale problems, the eddy diffusion analog to the
Fick Equation is often used.
The instantaneous point source solution of (21) is:
22)	C = —-— e-x2/4Dt
(4-rrDt) V2
which describes a Gaussian distribution curve about x^O, and
where M = Cdx. Since the Fick Equation is statistical in nature,

-------
39
high to low concentrations. It can be seen from (22) that as time
increases, the normal curve will flatten; the area under the curve
remains finite and equal to the total mass of marked particles
released at the source. It should be recognized that the solution
for predicting the concentration of a diffusing substance is also
a probabilistic equation, i.e., the mean concentration of marked
particles at x is directly proportional to the probability of
finding marked particles at that point. This idea is reinforced
by noting that a so-called Monte Carlo simulation of diffusion can
be obtained quite easily with a table of random numbers. The
longitudinal distribution of a substance introduced into a pipe can
be obtained by this method without employing the diffusion equation
at all (See Crank, 1955, p. 216, for an example).
Because of the feeling of uneasiness generated in scaling
up molecular analogs of diffusion to geophysical size, considerable
research has been devoted to more satisfactory descriptions.
Employed in the quality program is a form of the Kolmogoroff
hypothesis (Orlob, 1959) for computing the diffusion coefficient:
I = Empirical constant (dlmenslonless)
E = v.j •g-j^, an energy dissipation term*
*Note that for E constant the dispersion term becomes the "4/3 law,"

-------
= Channel velocity
g = Gravitational constant
AH/L^ = Potential (head) difference at ends of channel i
I = Scale of phenomenon; written in terms of the hydraulic
radius, R.
The dimensions (M, L, T) of the energy dissipation term are:
E = j(LT"1)(LT~2)J1^3 = L^V1
and of the scale terms are
I = {L}^3, so that
Dd = L2/3T_1.L4/3 = L^T"1
and (I is dimension!ess.
The diffusion of mass per time step is then
Ac.
Diff/At = C|Q|R -A
1
with dimension
MT"1 = (•)(LaT"1)(L)(ML"s)(L"1)
Ac
and where 	d is mass concentration gradient and Q is flow.
Li
As stated elsewhere, numerical errors can contribute to the
spreading out of material. When this is not accounted for, the
errors will be hidden in leading to erroneous conclusions as
to the relative magnitude of the advection and diffusion terms.

-------
41
term in the one-dimensional equation is river flow * cross-
sectional area. Diurnal tidal variations are then not implicit,
but are in reality responsible for producing the spread of
material with time such that peak concentrations occur both up-
stream and downstream of the source. This tidal displacement has
to be accounted for even in the tidaily-averaged equation and is
generally dumped into some form of the diffusion coefficient. This
coefficient is not a pure turbulent diffusion term, but is, rather,
a catch-all.
Stability, Numerical Mixing and Other Errors
A finite difference representation of differential equations
means that one will obtain solutions at discrete points at certain
time steps. Because of this and the fact that computing machines
carry only a finite number of decimal places, problems of trun-
cation, roundoff error, convergence and stability will always arise.
Certain of these concepts are stated concisely by O'Brien, et al.
(1951):
"Let D represent the exact solution of the partial differen-
tial equation, a represent the exact solution of the partial
difference equation, and N represent the numerical solution of the
partial difference equation. We call (D - a) the truncation
error; 1t arises because of the finite distance between points of

-------
42
the problem of convergence. We call (a - N) the numerical error.
If a faultless computer working to an infinite number of decimal
places were employed, the numerical error would be zero. Although
(a - N) may consist of several kinds of errors, we usually consider
it limited to round-off errors. To find the conditions under which
(a - N) is small throughout the entire region of the integration is
the problem of stability.
"The principal problem in the numerical solution of partial
differential equations is to determine N such that (D - N) is
smaller than some preassigned allowable error throughout the whole
region considered. We can assert that
{D - N) e (D - A) + (A - N)
is small for a numerical calculation over a fine mesh using a
stable, convergent difference scheme."
It should be noted that other definitions exist for truncation
and convergence. If a Taylor series expansion is used to approxi-
mate derivatives, only a few terms are carried; the higher ordered
terms are dropped and the series is said to be truncated. Like-
wise, a particular computing scheme may converge to a proper solu-
tion at a relatively fast or slow rate depending on the scheme
employed and the choice of initial conditions.
Hydraulic Equations
Two types of errors can occur in the programs under dis-

-------
43
program (HYDRA), stability is generally inferred from the so-
called "Courant Condition" for explicit finite representations of
the hydraulic (open-channel) equations. The Courant criterion can
be written as:
24> h »ivi ^ >4t •
where
Hmax 88 Max'imum channel depth
At = Integration step
The term /gHmax is the speed of a shallow water wave and holds
where the wave length is greater than twice the channel depth. The
approximation L- ^ v/giT At usually suffices in schematization as
is discussed later. "Wave length" refers to the length of a tidal
wave with a period that is on the order of 12.4 hours.
It has been found (See, e.g., Perkins, 1968) that even though
the Courant Condition is met, instability may occur and that this
instability is due to the presence of the non-11 near frictlonal
resistance term, K^V. (Vj, in the equation of motion. This term
is written:
n2|Vi
25) K^IV, =
where
(1.49)2 R1^5
Vi * K1V1

-------
44
R.j = Hydraulic radius of ith channel
= {(n/1.49)2 R,~4/3)
The modified Courant Condition is then written:
26) L. 5. IV. ±	g.K'lAt ,
which says that for a given integration step and channel depth,
the channel length must be at least of a certain length if stabi-
lity is to be maintained.
During the process of verifying current and stage, the
Manning coefficient can be adjusted in various reaches. This may
result in instabilities if n becomes too large, however, and a
shorter time step may become necessary or the schematization re-
examined.
Checks are available in the program to determine the serious-
ness of violation of water mass conservation resulting from
numerical procedures. These are discussed in part II.
Dispersion Equation
Two types of instability occur in the quality program.
Recently, attention has been directed to these aspects by various
authors (Orlob, et al., 1967; Bella and Dobbins, 1968; Prych, 1969).
Briefly, the problem occurs in the form of numerical errors

-------
45
is not conserved and a pseudo-dispersion of substance occurs. If
diffusion is included in the dispersion equation, the error is
masked as a longitudinal spreading of material in a manner that
appears to be a turbulent diffusion of the substance. If the
diffusion term is not included in the equation, then an initial
load distributed evenly throughout a given channel should move as
a self-contained parcel, i.e., it should not spread out with time.
Because the channel lengths and integration time steps are
fixed (in the analysis under discussion), the velocity in a given
channel times the time step (with resultant dimension as length)
may not exactly equal the particular channel length. If it were
exactly equal, there would be no problem, hence,a condition simi-
lar to the Courant Conditions for maximum stability would hold. In
essence, then, more material may be transported into a junction
than the junction can hold or more may be withdrawn than actually
exists. Program statements are written to prevent negative con-
centrations or this type of supersaturation for dissolved oxygen
concentrations. For other substances the statements are an
indication of instability and a determination of the seriousness of
the condition must be made.
The transport term AMps due to the "pseudo-dispersion"
phenomenon can be expressed as follows:

-------
46
where
K = Pseudo-dispersion coefficient (L2T_1)
ps
Ai = Channel cross-section area (L2)
aC/L^ = Concentration gradient (ML"4)
At = Time step (T).
The term KpS will depend on the particular difference scheme
employed and (Bella, 1969) can be roughly computed from:
vi
27) Kps = ji- [(1 - 2Y)L1 - Vjit]
y = 0 for a backward difference solution
= .25 for a quarter-point difference solution
=0.5 for a central difference solution
=1.0 for a forward difference solution.
While some choices of y will minimize AMpg, they may prove to con-
tribue to instability in the form of oscillations about the
solution points.
The pseudo-dispersion transport term is minimized in the
quality program by employing the quarter-point method which yields
"reasonably" accurate and stable solutions. Further testing of
the model with the diffusion term omitted and various difference

-------
HEAT BUDGET TERMS
Only a brief discussion of heat budgetry will be given here
since a certain familiarity with the subject 1s assumed and there
are any number of excellent texts and papers readily available.
Figure 7 illustrates the heat exchange processes at an earth
boundary during the day and at night. Similar magnitudes of
energy transfer components hold for water surfaces, except for the
evaporation term which may be relatively larger. It can be seen
that some processes considerably outweigh others; in the heat
budget formulation used here, this is reflected in the neglect of
radiation pseudo-convection and heat conduction. Also neglected
are terms for conduction of heat through the earth-water interface
and advection by rain.
In any given time period, the temperature at a fixed point
(Eulerlan analysis) will be raised or lowered or remain constant
depending on the heat balance of the heat budget terms and the
advection and diffusion rates during that time period. The time
period used 1n this discussion corresponds to the time step
employed 1n the temperature simulation (program QUALTEMP). The net
rate can be computed from empirical formulas; the formula summary
prepared by TVA (1967) was used 1n the following resume as a basic
reference. The system of units employed are: length in meters (M),
mass in kilograms (KG), time in seconds (SEC), pressure in milli-
bars (MB), temperatures in degrees centigrade (°C) and degrees

-------
UNIVERSAL SPACE
EXTRATERRESTRIAL
SOLAR RADIATION
lREFLECTIOI
\ FROM ,
CLOUDS/
ABSORPTION
DIFFUSE SCATTERING
BACK RADIATION
EFFECTIVE OUTGOING
S®	RADIATION
RADIATION FROM SKY
EVAPORATION
CONVECTION
RADIATIVE PSEUDO CONDUCTION
— HEAT CONDUCTION
SURFACE
TO GROUND
SUPPLIED
Energy balance at noon on a sunny day. The
width of the arrows is proportional to the
amount of energy transferred.
EFFECTIVE OUTGOING
RADIATION
EVAPORATION
SURFACE
BACK RADIATION
CONVECTION
RADIATIVE PSEUDO CONDUCTION
THERMAL CONDUCTION
FORMATION OF DEW
SUPPLIED FROM
GROUND
Energy balance at night drawn to the same
scale as above.
FIGURE 7. Energy Balance Terms (From R. Gieger, "The Climate near
the Ground," pp. 7, Harvard University Press, Cambridge,

-------
49
Heat Flux Through the Water Surface, q^, (KCAL M 2SEC 1)
«H =
qsn + qatn + qw + qe + V where
<"sn "
Net solar radiation flux, +:incoming
''atn
Net atmospheric radiation flux, +:incoming
% =
Water surface radiation flux, -:outgoing
qe =
Evaporative heat flux, outgoing
qc =
Convective heat flux, ±neither way, depending on

the difference in air and water temperature

(+. 1f Ta>Ts) .
The first two terns on the right of (28) are discussed only
briefly. They are quite complicated functions, but easily comput
able. It is assumed that the available meteorological programs
have been made use of to obtain qsn + Qatn *or sPec^° t-®mes an<*
geographic locations under discussion, or that direct measurements
independent of the water sur-
are available. These two terms are inaeyci
* x. ^ , ,i,	n 1 and can be computed by an
face temperature (unlike qfi, qw> ana
external program not necessarily linked to that under discussion.
Since a one-d1mens1onal model 1s employed, all net Incoming
radiation 1s absorbed at the surface; U 1s distributed evenly
throughout the water column via the one-dimensional assumption.
The temperature dependent heat budget terms are computed at

-------
50
is used in the formulas for q , q » q , to obtain new values which
c w w
are in turn summed with the independent terms to compute a new net
flux.
Temperature Dependent Terms, - Computation
The dependence of the surface temperature is direct for q
W
and q and somewhat indirect for q as can be seen in the following
v	C
approximations:
29)	qw - a-(Ts + 273.16)4
30)	qc = b-(Ts - Ta)
31)	qe = c-(es - ea) ,
where
T = Surface water temperature (°C)
d
T, = Air temperature (°C)
a
e& = Pressure of saturated water vapor at temperature T
ea = Pressure of water vapor in ambient air
a,b,c ® Empirical coefficients.
Back Radiation, qw, (KCAL M^SEC"1)
All bodies emit radiation at a rate proportional to the
fourth power of the absolute temperature (TQ) of their surface.

-------
51
back radiation term, qw- The surface radiation formula is:
32)	qw = e-a«T S where
e = 0.97, the emissivity
a = 1.36 x 10"11 KCAL M"2SEC'1 °)f\ the Stefan-fiol tzman
constant.
Evaporation Heat Exchange, qe* (KCAL M~2S£C_1)
Heat loss by the vaporization of water is expressed by:
33)	qe « PW*E'HV, where
E * Rate of water loss due to evaporation, M SEC"1
HV » Latent heat of vaporization, KCAL KG"1
pw a Water density, 1000 KG M"3,
E is computed by means of the formula:
34)	E ¦ N*U*(e_ - ea), where
S 
-------
52
The heat vaporization term (KCAL KG ) is written:
HV = 597. - 0.57-Ts.
The vapor pressure terms (MB's) are computed through
exponential approximation formulae first employed by Lamoreaux
(1962) through the Clausius-Clapeyron equation. The coefficients
used are taken from a WRE report (1969):
es = 2.1718 x 108 exp(-4157.0/(239.09 + Tg))
ea = 2.1718 x 108 exp(-4157.0/(239.09 + Twb))
"AP(Ta " Twb)(6>6 x 10~* + 7"59 * 10"' (Twb"'
where
Twb = Wet bulb temPerature (°c)
AP = Air pressure (MB).
Convection Heat Exchange, qc, (KCAL M"2SEC_1)
Convective exchanges, as sensible heat transfer, far outweigh
conduction heat exchanges (which are neglected).
Although direct measurements of both q and q are possible,
B	v
their measurement is quite complex due in part to instrumentation
difficulties and the necessity to somehow measure turbulent flux
terms (which are masked in transfer coefficients).
The method used here is to employ the Bowen ratio:
-------
53
since qg is easily computed (but not necessarily an accurate
estimate) qc can be evaluated through:
35) qc = BR-qe .
The Bowen ratio is computed as follows:
BR = 6.1 x lO'^-AP-
Ta"Ts
es " ea
Summary of Heat Budget Step
Initial conditions are used to compute the dependent heat
budget terms; these are summed algebraically and added to the
independent terms. The net flux (q^, which will be zero, positive
or negative, depending on the relative magnitude of the terms)
during a computation step (1 hour, here) 1s multiplied by the time
step (At) divided by density, depth (d) and specific heat (Cp) and
added to the most recent temperature term:
36) Tnew = Told +	'
During the next computational Interval, Tnew becomes T^t
advection and diffusion steps and time changes in depth are applied
1n the program just prior to the net heat flux step.
Equll1br1um Temperature, Te, °C
For a check on the temperature as computed above or as a
-------
54
temperature" approach. The most recent work on this subject has
been conducted by Edinger, Geyer and associates whose publications
(1965, 1967, 1968, e.g.) should be examined for a complete des-
cription of the subject. Briefly, the equilibrium temperature
method is a shortened approximation to the net heat transfer
method outlined above in that linear approximations to the vapor
pressure and back radiation terms are employed. Temperature
estimates can be made rather rapidly using a desk calculator if
single water parcels (Lagrangian analysis) are dealt with.
An option is provided in the program to compute the exchange
coefficient, equilibrium temperature, and the water parcel
temperature according to the equation:
37> Ts = Te + (Told-Te» exP(-^T3-).
where
T = Equilibrium temperature, °C
V
K = Thermal exchange coefficient, KCAL M"2SEC"1 °C~1,
Equation 37 could be used in itself for an analysis where the
coordinate system moved with the water parcel; it is known as the
-------
SCHEMATIZATION
General
Details of the schematizatlon of the Columbia River under tidal
influence are described 1n detail below to exhibit the geographic
and hydrologic input data for the model. The total schematizatlon
(Figures 8-11) consists of 396 finite elements called "junctions,"
each of which is an arbitrarily-shaped area centered about a junction
point; the junctions are connected by 432 "channels." The large
scale work charts of the schematization shown in Figures 8-10
which include the numbering system and other detail are on file at
our laboratory. Part II explains how to select boundary conditions
such that only portions of the schematization need be used.
Base Charts and.Data Sources
The schematizatlon was prepared primarily from U.S. Coast
and Geodetic Survey navigation charts numbers 6151, 6152, 6153,
6154, and 6156, scale 1:40,000. Coast and Geodetic Survey Chart
number 6155, scale 1:20,000, and U.S. Army Corps of Engineers
dredge sheets, scale approximately 1:5040, were also used to obtain
geographic data for selected areas of the river system.
Flow data were obtained from records of the U. S. Artqy Corps
of Engineers, Portland District, and summarized by FWPCA personnel
-------
STATUTE MILES
RM 28.
Jim Crow Point
RM 0
.CHANNELS
'JUNCTIONS
FIGURE a COLUMBIA RIVER SCHEMATIZATION
RIVER MILE 0.0 (PACIFIC OCEAN) TO
-------
Jim Crow Point ^See figure II for exploded view
STATUTE MILES
?o
RM 28
Island
CHANNELS
FLOOD

FUNCTIONS
' Cowlitz Rtwr
Longview
l I I
m
RM 50
STATUTE MILES
Cowlitz River
-Lewis
wer
¦I 1 I iw > •
STATUTE MILES
FIGURE 9. COLUMBIA RIVER SCHEMATIZATION
RIVER MILE 28.3 (JIM CROW POINT)
-------
Portland
.< Willamette
River
River
FLOOD
Sauvies Island
St. Helens ¦ ~
TUTE MILES
CHANNELS
UNCTIONS
Vancouver
RM 112,
WashougaI River

rtland
Willamette
River

*xSon
120.5
STATUTE MILES
Bonneville Dom
STATUTE MILES
Washouga I River
ANNEL NUMBERS
(SEE TEXT)
Sandy
River
RIVER CURRENT
FIGURE 10 COLUMBIA RIVER SCHEMATIZATION
RIVER MILE 87.0 (LEWIS RIVER) TO
-------
WASHINGTON
Au	y.-..',-.-




'-1V k* .\V|-U <«
%< tf/K/J* *
XXXII

OREGON
FLOOD
FIGURE II. COLUMBIA RIVER OREGON-WASHINGTON SCHEMATIZATION
RIVER MILE 28 TO RIVER MILE 35 (EXPLODED VIEW OF
-------
60
Tidal information was taken from Coast and Geodetic Survey
tide tables and navigation charts.
Datum Planes
Stream depths on the navigation charts are referenced to mean
lower low water at lowest river stages, Columbia River Datum; mean
sea level is used as the datum plane.
River Boundaries
Heavy black lines on the base charts mark the river bed
boundaries. These lines were traced on overlay paper to provide
an outline of the river. All islands, lower reaches of major
tributaries, and major sloughs were outlined. Minor tributaries
and sloughs were not included.
An exception to these boundaries is from river mile 125 to
river mile 146. In this segment, the boundaries of the river were
taken at the M.L.L.W. level; that is, the first dotted line within
the riverbed proper. On the charts, these Tines separate inter-
tidal areas (green colored) from the water (blue colored). The
decision to use this line as a boundary in this reach of the river
was based on the assumption that tidal influence was minimal in
this reach and that the primary use of the model would be to
-------
61
Junction Point Layout
The overlay was placed on the base chart and the apparent
main flow routes sketched. The depth and width of the riverbed,
islands, etc. were considered in this preliminary sketch. After
it was felt the flow pattern was reasonably represented, the
junction points were plotted on the overlay. Establishing these
points (as well as sketching in the flow) is somewhat of an art;
however, certain criteria must be met. These are discussed in the
following section.
Boundaries were then established around the junction points.
That part of a boundary between junction points connected by a
channel was drawn across the channel, near its midpoint, to the
edge of the estuary schematization or to the boundary separating
unconnected junction points. Nominally these bounds are perpen-
dicular to the channels. In the narrower parts of the schematiza-
tion and in those areas where junction points lie near the side of
the schematization, the schematization boundary forms part of the
junction boundary. In the wider parts of the schematization
bounds between unconnected junction points were somewhat arbitrary
and generally were drawn along intertldal areas and shoal areas,
Criteria for Selecting Junction Points
The selection of junction points and the distance between
-------
62
an "average" channel depth between junctions. As stated earlier
the maximum average depth, Hmax, determines the speed of a shallow
water wave according to /gHmax. For a given integration time, At,
the channel length, L.., is limited according to the quantity-^
Li 2 4t •/9f& *
As mentioned before, the first "free hand" schematization is
based on a pre-selected At. The actual channel depths, the areas
of interest, the divergency and convergence of channels, and the
detail one wishes to go into enter into the selection of junction
points. Once this first selection has been made, it is possible
to compute the required topographic information from a knowledge
of the mean depths in each junction and the surface areas of the
junctions.
Data Obtained from the Schematization
Upon completion of the schematization, pertinent input data
for the model program were obtained from each junction and each
channel.
Each junction has as input data: a number, from one to five
channels connected to it, a surface area and an initial head. In
U The relation is more complicated as was discussed earlier;
however, this simple formula was used to estimate the successive
-------
63
addition, those junctions located where a tributary enters the
river has as input data the flow of that tributary.
Each channel has as input data: a number, two junctions
connected to it, a width, a depth, an initial streamflow velocity,
and a Manning coefficient.
Derivation of these data from the schematization and other
pertinent records are discussed in the following sections.
Junction Data
Junction Numbers
There are 396 junctions or junction points in the entire
scheme numbered from 1 through 396 inclusive. The schematization
was prepared in two sections. Section I extends from river mile
28.3 (Jim Crow Point) to river mile 146.1 (Bonneville Dam) and the
junctions are numbered consecutively from 1 though 260.
Section II of the scheme extends from river mile 0 (Pacific
Ocean) to river mile 28.3. To facilitate the location of starting
points 1n this part of the scheme, the two junctions at the sea-
ward end were numbered 1 and 2, respectively. Those junctions near
river mile 28.3 which had been numbered 1, 2, and 3 in Section I
were renumbered 261, 262, and 263, respectively.*
*It would have been easier to rewrite the program to handle any
numbering system at the ocean end. This has been done by D. Fitz-
gerald (Northeast FWPCA Regional Office) for Boston Harbor;
-------
64
Additionally, junctions in Section II were renumbered 1
through 142 inclusive in order that Section II of the scheme could
be operated independently of Section I.
Junction Surface Area
The surface area of each junction, at mean tide level datum, was
measured with a planimeter and recorded in square feet.
Those junctions which have major sloughs entering them have
included in their surface area the surface area of those sloughs.
Initial Head
The initial heads for each junction (the approximate height of
the water surface at a flow of 147,200 C.F.S. at Bonneville Dam)
were obtained from Corps of Engineers records. A graph of the
heads at selected river mile intervals was prepared and the appro-
priate datum taken from the graph for each junction. Initial heads
could have been taken as 0.0 throughout at the expense of a delay
in convergence in the iteration.
Number of Channels at a Junction
From one to five channels may enter each junction. The
number of each channel entering is listed as input data. The
lowest numbered channel entering is listed first, the highest
numbered channel last. More than five channels may be accommodated
-------
65
Junction Depths
After the schematization was prepared, the mean depth of each
junction was determined.* The technique for doing this 1s described
below:
A transparent grid overlay consisting of 225-600 foot square
squares, scale 1:40,000, was prepared. The steps given below out-
line the procedure for finding the depth of a junction.
1.	The grid was placed over a junction outline on the base
chart.
2.	The depth at the center of each 600-foot square was read
and recorded. A detailed explanation of this procedure is given in
the paragraphs following these steps.
3.	The square was marked with a grease pencil and counted.
4.	In each junction, there were always some grid squares
that fell on the junction boundaries, putting only parts of the
grid squares within the junction boundary. These parts of grid
squares were summed mentally to make a whole grid square and the
depth estimated and recorded.
5.	The sum of the squares read for each junction divided
into the sum of the depths gave the mean depth for the junction.
*The depths thus calculated are not used directly 1n the program
but were made to provide an independent deck of depths computed
-------
66
The procedure and the data entered on each card are described
below.
1.	The junction number was read and entered on the punch
card in columns 1-3.
2.	The card number was entered in column 4. Most of the
junctions required that more than one card be used to record all
the depth readings. The cards required for each junction were
numbered sequentially from 1 through the number required.
3.	The size of the grid squares being used was entered in
columns 5-7.
4.	The stage correction was entered in columns 8-10.
The stage correction was applied to depths obtained from the
chart, referenced to mean lower low water to obtain a depth
referenced to mean tide level. The stage correction varied for
different reaches of the river. Near the seaward end, it was
4 feet; in the reach immediately below Bonneville Dam, it was
taken as 0. The stage correction was made to the nearest whole
-------
67
TABLE 1
STAGE CORRECTIONS VS. RIVER MILES
Stage Correction	River Miles
4 feet	0-28
3 feet	28 - 50
2 feet	50 - 76
1 foot	76 - 122
No correction	122 - 146
5.	The number of channels entering a junction was entered
in column 11.
6.	All of the data in columns 1-11 were entered on each
card being used for a particular junction.
7.	Depth data, read directly from the base charts, for each
grid square was entered in columns 12 through 80. Three columns
were used for each reading.
8.	As noted in 4 above, all depths read directly from the
chart were referenced to mean lower low water. This situation
caused intertidal areas, shown In green on the charts, to be above
the datum from which depths were read. In order to accommodate
-------
68
correction applicable to that particular junction, was read by the
reader when such an area occurred under a grid square. The depth
was entered on the card with a 90 preceding it. For example, the
intertidal area of a junction in that reach of the river having
a stage correction of 3 feet would be entered on the card as 903.
9.	After all the grid squares had been read and accounted
for, a 999 was entered on the card to indicate the end of data for
that junction.
10.	Frequently, the reader and the keypunch operator would
change roles and the junctions would be read a second time.
11.	The two independent mean depths were compared; if the
difference between them was less than two feet, the mean of the
two readings was taken as the junction depth. If the two depths
varied by two feet or more, the junction was read one or more
times to obtain a usable junction depth.
Tributary Stream Flows
Flow data for several tributaries which enter the Columbia
River were available from Corps of Engineers records. These data
were entered as input data for the junctions in which the tribu-
-------
69
Channel Data
Channel Numbers
Channels in Section I of the scheme were numbered from 1
through 276 inclusive.
In Section II of the scheme, the two channels near the seaward
end of the scheme were numbered 1 and 2, respectively. Those
channels near river mile 28.3, which had been numbered 1, 2f and
3, were renumbered 277, 278, and 279, respectively, Additionally,
the channels in Section II were renumbered from 1 through 159
inclusive, in order that it could be operated independently of
Section I.
Channel Length
The length of each channel between two connected junctions
was measured in feet and ranged from about 2,000 feet to about
12,000 feet.
Channel Depths
The depth of each channel is taken as the mean of the two
depths of the junctions which that channel connects. It was felt
that the preliminary smoothing effected by this averaging would
compensate for channels lying partly in deep water and partly 1n
-------
70
Channel Widths
The widths of each channel were measured along the junction
boundary which crossed a channel near its midpoint. Widths were
measured in feet. Widths were measured at both M.L.L.W, and M.T.L.
Cross-sectional Area
In Section I, cross-section areas were constructed along each
junction boundary crossing a channel, planimetered, and reported
in square feet.
Cross-sectional areas in Section II were obtained by multi-
plying the M.T.L. width of a channel by its mean depth referenced
to the appropriate datum.
Channel Flow
Streamflow in each channel was used to calculate initial
velocities. Arbitrary initial velocities could also have been
used; the extra work involved here was felt worthwhile in order to
reduce the possibility of instability due to a bad choice of
initial conditions which might have been difficult to correct.
The total flow in the river (at the mouth) was taken as the
sum of the flow at Bonneville Dam plus the flow from tributaries,
for which data were available, during the modeling period.
Flow in each channel was derived from the flow in the
channel immediately upstream from it plus any flow entering from
-------
71
For example, in channel number 274, immediately below
Bonneville Dam (see Figure 10), the flow is 147,200 C.F.S. (meas-
ured at Bonneville).
The flow remains constant for all channels downstream through
number 265. This channel branches into channel numbers 264 and
270, respectively.
To find the flow in each of these channels, a straight line
partitioning was done in the following fashion.
The sum of the cross-sectional areas of channel numbers 264
and 270 was found and the percentage each channel contributed to
this total was calculated. This percentage was then multiplied by
the flow in channel number 265 (the branching channel) to give the
flow in 264 and 270, respectively.
Flow in channel #265 = 147,200 ft3/sec.
Cross-section area channel #270 * 21,714 ft2 * 36% of total
Cross-section area channel #264 * 38,164 ft2 s 64% of total
Total « 59,878 ft2
Flow in channel #264 = (0.64)(147,200 ft3/sec) » 93,800 ft3/sec.
Flow in channel #264 = (0.36)047,200 ft3/sec) ¦ 53,360 ft3/sec.
The flow in channel numbers 267, 268, and 269 remain the same
as the flow in channel 270. Similarly, flow in channel numbers
261, 262, and 263 remain the same as the flow in channel number
-------
72
The flow in channel number 260 is the sum of the flow in
channel numbers 261 and 267.
Similar partitioning and summing of flows was done throughout
the scheme; about four hours were required to complete the entire
channel initialization.
Channel Velocity
The initial water velocity (owing to streamflow) in each
channel was found by dividing the flow in that channel by its
cross-sectional area. If a channel ended in a slough or in a
tributary with no recorded flow data, the velocity was set to
zero.
Computing the velocities in this manner resulted in what
were to be unrealistic velocities in some channels, on the order of
10 feet per second. When such velocities occurred, the width or
depth of the channel was arbitrarily reduced an appropriate amount
to make the flow realistic. Such changes were made in channel
numbers 69, 70, 71, 72, 85, 206, 222, 234, 235, 236, 237, 240,
-------
DISCUSSION
Deterministic pollution models of the environment usually
involve analytical or numerical solutions of the dispersion (or
advection-diffusion) equation. The numerical methods employed
are likely to be identical regardless of the constituent involved,
hence a generalized model should obviate the necessity of deriving
a new model for different topographical settings and constituents.
This being the case an existing model was modified to handle
temperature; the solutions obtained are thus forms of the energy
equation.
Where Coriolis terms are unimportant and stratification
(either vertical or horizontal) is slight, the finite element
representation of two-dimensional environments may be quite satis-
factory. If the Coriolis force is not negligible, then the methods
used will not suffice since the velocity term in the y-direction
is required in a solution of the x-direction equation of motion.
Since many open estuarine areas consist of numerous scoured
channels, the junction-channel representation of these areas may
not be as forced as it may first appear. Certainly, the use of a
one-dimensional approach to sections of a tidal river may be ques-
tioned, but it is also questionable if a fine grid model incorporat-
ing horizontal shear terms would add a great deal to our present
-------
74
verification procedure which would be required for such a model,
especially under different wind, tide and runoff conditions. It
is well known that the only real limitation we have in the complete
solution of the Navier-Stokes equations is one of computer hard-
ware. Doubtless we will have mind-bogglingly fast machines with
almost unlimited storage capacity sometime in the near future, but
the big question is likely to remain on how to handle and verify
the rather simple models we have even now. Such problems will
always face the model user; if the uses of modeling are to be
well served, verification will go hand-in-hand with modeling use.
In flood routing problems and for high-accuracy displays of
periodically exposed tidal flats, procedures must be employed to
allow for time-varying channel widths. Where very accurate repre-
sentations of tidal flow are required, it may not be enough to
vary only the cross-section; in this model, however, rectilinear
channels were assumed and cross-section variation occurs only by
a change in stage elevation. The two-dimensional model employed
by Leendertse (personal communication) on Jamaica Bay, New York,
apparently accounts for tide flat exposure every five time steps.
Such models are highly desirable if not absolute necessities in
shallow estuaries such as Tillamook Bay and other areas where
only a stream cuts through extensive tide flat areas at low water.
-------
75
the modeler or manager and perhaps less rigorous approaches may
suffice for certain aspects of a particular pollution model. It
should be noted that a model such as Leendertse's or the one
described here requires several years of continuous development,
and pollution agencies usually operate on far more demanding time
scales.
In the matter of the heat budget calculations (where Part II
provides examples) it is known that the heat exchange process at
the surface is much more involved than would be implied by the
equations employed. Air-sea interaction occupies a large area of
research in the oceanographic and meteorological community and
involves studies of the flux of heat and momentum to and from the
atmosphere. Turbulent processes at the surface are still rather
mysterious, so ultimate solutions of heat budget processes are not
likely in the immediate future; the approximations employed,
however, have given reasonably satisfactory answers where verifi-
cation has been possible.
Of more concern in certain areas, particularly the stratified
marine environment, is the role that vertical velocity and density
variations play in the overall pollution dispersion problem. Here
again, field studies with an end to verification are major under-
takings. As an indication of the upper part of the scale, field
Programs conducted by the U. S. Arnjy Corps of Engineers for the
-------
76
approximately $250,000. Aside from the usually back-breaking pro-
cess of field collection is the general inadequacy (in terms of
ease of use and reliability) of water quality measuring devices.
(If the data has to be collected over one or more tidal periods,
it is usually a toss-up as to whether the electronic gear will
give out before the field personnel do.)
Finally, if it is assumed that the well-planned survey goes
off without a hitch and the measuring devices do not balk, data
reduction and analysis will surely manage to contain unplanned for
and/or uncorrectable situations.
So much for executing the faultless survey; what is one to
do? Short of designating the problem of verification and field
collection as someone else's business, it behooves the model user
and builder to be aware of what goes into the various terms and
coefficients in order that they may be properly sampled or esti-
mated at the appropriate time. He should also be aware of the
realities and limitations of field techniques and existing instru-
mentation, as well as being aware of possible alternative solutions
-------
ACKNOWLEDGEMENTS
What was thought to be a fairly polished draft of the manu-
script was reviewed by Dr. D. J. Baumgartner, Mr. K. Feigner,
Dr. B. A. Tichenor and Mr. J. Yearsley, all of the FWPCA. They
supplied many constructive, well-taken and, hopefully, well-
received comments which vastly improved the overall effort. Any
technical errors or grammar violations that survived remain with
-------
REFERENCES
Bella, D. E. and W. E. Dobbins, 1968, "Difference modeling of
stream pollution,'1 J. San. Enq. Div. ASCE, 94, No. SA5,
pp. 995-1016.
Bella, D.E., 1969, Tidal flats in estuarine water quality
analysis, Progress Report for FWPCA Research Grant
WP-01385-01, Dept. of Civ. Eng., Oregon State Univ.
9/30/69. 45pp processed.
Baltzer, R. A. and C. Lai, 1968, "Computer simulation of unsteady
flows in waterways," J. Hyd. Div. ASCE, 94, No. HY4,
pp. 1083-1118.
Crank, J., 1956, The mathematics of diffusion, Oxford Univ. Press,
London, 347 pp.
Dronkers, J. J., 1964, Tidal computations in riyers and coastal
waters, John Wiley and Sons, Inc., New York, 518 pp.
Dronkers, J. J., 1969, "Tidal computations for rivers, coastal
areas, and seas," J. Hyd. Div., Proc. Amer. Soc. of Civ.
Eng., 95, No. HY1, January 1969, pp. 29-77.
Edinger, J. E., D. W. Duttweiler, and J. C. Geyer, 1968, "The
response of water temperatures to meteorological conditions,"
Water Res. Research, 4, 5, pp. 1137-1143.
Edinger, J. E. and J. C. Geyer, 1965, Heat exchange in the
environment, Edison Elect. Inst., Pub. 65-902, New York.
Edinger, J. E. and J. C. Geyer, 1967, Analyzing stream electric
power plant discharges. Proc. Nat. Symp. on Est. Poll.,
Stanford university, August 1967, pp. 462-485.
Fischer, H. B., 1967, "The mechanics of dispersion in natural
streams, J. Hyd. Div., Proc. Amer. Soc. of Civ. Eng., 93,
No. HY6, November 1967, pp. 187-216.
Hansen, Donald V., 1965, Currents and mixing in the Columbia River
estuary, Ocean Sci. and Ocean Eng. Trans, of the Joint Conf.
Mar. Tech. Soc. and Amer. Soc. of Limn, and Ocean., Washington.
-------
80
Hansen, D. V. and M. Rattray, Jr., 1965, "Gravitational
circulation in straits and estuaries," J. Mar. Res.,
pp. 104-122.
Hansen, W., 1966,"The reproduction of the motion in the sea by
means of hydrodynamical-numerical methods," Mitteil. Inst.
Meeresk, Hamburg, 5i, 57pp.
Hildebrand, 1956, Introduction to numerical analysis, McGraw-Hill
Book Co., New York, 511 pp.
Lai, C., 1966, "Discussion of 'Computer simulation of estuarial
networks'". J. Hyd. Pi v., ASCE, No. 3, pp. 96-99.
leendertse, Jan J., 1967, Aspects of a computational model for
long-period water-wave propagation, Memo, RM-5Z94-PR,
The Rand Corporation, Santa Monica, California, 165pp.
Morse, W. E., 1969, Stream temperature prediction model,
Presented at AGU Annual Regional Conference, Portland, Oregon,
October 16-17, 1969.
O'Brien, G. G., M. A. Hyman, and S. Kaplan, 1951, "A study of the
numerical solution of partial differential equations,"
J. Math, and Phys., 29, pp. 223-251.
O'Connor, D. J., J. P. St. John, and D. M. DiToro, 1968, "Water
quality analysis of the Delaware River estuary," J. San. Enq.
Div. ASCE, 94, No. SA6, pp. 1225-1252.
Okubo, Akira and M. J. Karweit, 1969, "Diffusion from a continuous
source in a uniform shear flow," Limn. & Ocean., 14, 4,
pp. 514-520.
Orlob, G, T.. 1959, "Eddy diffusion in homogeneous turbulence,"
J. Hyd. Div., ASCE, 85, No. HY9, pp. 75-101.
Orlob, Gerald T., 1968, Estuarial system analysis quantity and
Duality considerations, Proc. Nat. Symp. on the Anal, of
ater-Resource System, July 1-3, 1968, Denver, Colorado,
pp. 341-358.
Orlob, G. T., R. P. Shubinski, and K. D. Feigner, 1967,
Mathematical modeling of water quality in estuarial systems,
P**oc. Nat. Symp. on Est. Poll., Stanford University,
-------
81
Perkins, F. E., 1968, The role of damping on the stability of
finite difference schemes, ASCE Envir. Eng. Conf,,
Chattanooga, Tennessee, April 1968, 12 pp.
Prych, Edmund A., 1969, "Discussion," J. San. Eng. Div., Proc.
Amer. Soc. of Civ. Eng., 95, No. SA5, October 1969,
pp. 959-964.
Shubinski, R. P., J. C. McCarty, and M. R. Lindorf, 1965,
Computer simulation of estuarial networks, ASCE Water Res.
Eng. Conf., Mobile, Alabama, March 8-12, 1965, Con. Preprint
168, pp. 1-28.
Shubinski, R. P., J. C. McCarty and M. R. Lindorf, 1967, "Closure
to 'Computer simulation of estuarial networks'", J. Hyd. Div..
ASCE, No. 1, pp. 68-69.
Shubinski, R. P. and C. F. Scheffey, 1966, Wave propagation in
estuarial networks, Proc. Sec. Aust. Conf. on Hyd. and F1.
Mech., Univ. of Auckland, N. Z., pp. A81-A96.
Streeter, H. W. and E. B. Phelps, 1925, A study of the pollution
and natural purification of the Ohio River, Public Health
Bull. 146, U.S. Public Health Service, Washington, D.C.S 75 pp.
Sverdrup, H. U., M. W. Johnson, and R. H. Fleming, 1942,
The Oceans, their physics, chemistry, and general biology,
Prentice-Hall, Inc., New York, 1087pp.
TVA, Division of Water Control Planning, Eng. lab., 1967, Heat and
mass transfer between a water surface and the atmosphere,
Revised, May 1968, 98pp. processed.
Thomann, R. V., 1963, "Mathematical model for dissolved oxygen,"
J. San. Eng. Div. ASCE, 89, No. SA5, 30pp.
U. S. Geological Survey, 1952, Water-loss investigations. Vol. 1 -
Lake Hefner Studies, U. S. Geological Survey Circ. 229, ~
153 pp.
Wada, A., 1967, Study on recirculation of cooling water of power
station sited on a bay, Japan Soc. of civ. Enar. Ifl. lQfi7.
pp. 143-170. 		~ '
Water Resources Engineers, Inc., 1965, A water Quality model of
the Sacramento-San Joaguin D6lta, keport of an investigation
-------
82
Water Resources Engineers, Inc., 1966, A hydraulic-watfer quality
model of Siusun and San Pablo Bays% Report to the FMPCA.
Southwest Region, March 1966, 34pp.
Water Resources Engineers, Inc., 1969, Mathematical models for the
prediction of thermal energy changes in impoundments, MSS
-------
APPENDIX I
NOTES FOR PROGRAM HYDRA
"...free from bugs,...
if possible,
If you know any such."
-------
85
APPENDIX I
NOTES FOR PROGRAM HYDRA
Logical unit (or data set reference) numbers:
A time sharing conversational computer system was used for
much of the work on the programs done in connection with this
paper. Several separate data files were created for different
portions of the input data; the program then referenced each file
with a different unit number.
Most operating systems have facilities for equating different
unit numbers to the same source. If card input is used, units
listed below as card image input should be equated to the card
reader.
Unit Number	Use
5
images)
Standard Output (printer)
Junction Input (card images)
Channel Input (card images)
Restart Output (card images)
Output for HYDEX (binary tape images)
Control and System Input (card
6
7
8
9
-------
86
PROGRAM NOTES
Program HYDRA
Sequence No. FORTRAN Name
Dimensioned variables
AREA
AREAT
V
VT
Y
YT
Dimensioned variables
AK
ALPHA
AREAS
6
CIEN
CN
Comments
In the Runge-Kutta integration
scheme, the variable quantities
head (Y), velocity (V), and
cross-sectional area (AREA), are
extrapolated forward a half-
cycle. The extrapolated values
are held in arrays YT (for Y
temporary), VT, and AREAT. The
temporary values are then used
to compute the new values at
the time one cycle forward.
Heads are junction properties,
velocities and areas are
channel properties.
Constant used in obtaining
frictional resistance term =
(G*CN**2/2*21), see line 146.
Alphanumeric information used
to identify printout.
Surface area of junction
determined by, e.g., plani-
metering.
Channel width.
Channel length.
Friction coefficient.
JPRT
Junction numbers where data
-------
87
Sequence No. FORTRAN Name
NCHAN
NJUNC
Q
QIN
R
Variables in COMMON
DELT
NC
NCYC
NCYCC
NJ
NOPRT
NPRT
PERIOD
Comments
Channels (maximum of 5)
connected to junction J., e. g.,
NCHAN (J,l) » channel number
of first channel connected to
Junction 0, etc.
Junctions (maximum of 2)
attached to either end of channel
N, e.g	NJUNC (N.I) =
junction number at one end of
channel N, etc.
Channel flow.
Tributary inflow.
Hydraulic radius («AREA/B).
Time increment (seconds) used
as integration step.
Total number of channels.
Number of cycles hydraulic
program will run.
Used to hold the current cycle
for use outside the main loop.
Total number of junctions.
Number of junctions for which
data is to be printed.
Printing output interval, i.e.,
results are printed every
NPRTth cycle,
Period of input wave at ocean
-------
88
Sequence No. FORTRAN Name
Comments
19-24
20
20
20
Variable names listed as they occur in the program:
Read input information and
write a general heading.
20
23
23
23
25-30
27
34-40
TZERO
NETFLW
ISTART
INPSUP
IPRT
IWRTE
KPNCHI
I STOP
Initial time of this run.
Switch to call HYDEX: ?0.calls.
If the run is a restart, this
is the next cycle to be pro-
cessed. (= to 1 + the cycle
number written in the run to
be restarted.)
If j*0, suppresses printing
junction and channel data.
Start printing on cycle	.
Write binary output on 10
beginning at cycle	.
The interval at which the
restart records will be
written.
Position tape 10 if this is a
restart.
The number of the last cycle
processed in the run to be
restarted.
Read junction initial conditions;
unless this is a restart, set
YT=Y; determine if the input
data cards are in the correct
order.
41-44
Write junction input data
-------
89
Sequence No. FORTRAN Name
48-53
54-58
62
63-64
68	NEXIT
69-88
92-101
105	DELT2
106
107
108	W
Comments
Read channel initial conditions,
calculate channel cross-
sectional area (see text),
determine if the input data
cards are in the correct order.
Write channel input data unless
INPSUP *0.
Read the numbers of the junctions
for which data will be printed.
Read and write the amplitude,
phase and period of the input
tide wave.
A flag which is non-zero if
there is a compatibility or
restart problem.
Determine that the junction
and channels are connected to
each other properly.
Output the control and system
data to unit 10. If IWRTE is 0,
(hydraulic output desired for
every cycle) calculate channel
flows and output Initial head,
speed and flow to unit 10.
Half a time step,
Convert the Initial time to
seconds.
Convert the period in hours to
seconds.
Used in the Fourier series
representation of the tidal
-------
90
Sequence No. FORTRAN Name
109-119
Comments
The restart provision has two
purposes: if the program
terminates abnormally, (over
time estimate, system failure,
etc.) a restart record can be
used to start the program at
some mid-point without wasting
all the computer time used in
getting to the point. A restart
record is also made at normal
termination, so that the run
can be extended if desired.
Writing the restart record
itself uses extra time. If
abnormal termination is not a
problem, set KPNCHI=0, and the
restart record will be made
only at the end. To protect
against abnormal termination,
set KPNCHI>200 or so. A
restart record will be made
every KPNCHIth cycle. To use
the restart unit, equate it to
units 7 and 8, and run the pro-
gram using cards only for the
unit 5 read statements, get
TZERO from the last line of
printed output saying, "TZERO
F0R RESTARTING = ..." Set
ISTART = 1 + last cycle for
which restart "deck" was made.
120	If there has been a compatibi'
lity problem, or an error in
KPNCHI, stop.
121	G	Gravitational constant.
125-131	Calculate the frictional con-
stant in each channel; reorder
the junction numbers if
necessary.
132	T	Set time equal to the initial
-------
91
Sequence No. FORTRAN Name
140
140-143
147-156
150-151
152-155	DVDX
160
161-170
171
175
Comments
Loop through sequence number
289 (statement 285) for each
hydraulic cycle, through NCYC
cycles.
Replace current values of
NCYCC, T2 and T each time
through the loop. See section
on Runge-Kutta integration,
explicit solution, leapfrog
methods.
Compute channel speeds and
f1 ows.
Divide the frictional constant
by the hydraulic radius raised
to the 4/3 power after comput-
ing the current value of the
hydraulic radius.
Compute the velocity gradient
(dv/dx) and channel velocity
from initial or last computed
values of velocity and head.
Compute the Fourier series
representation of the input tide
wave for as many junctions as
are at the ocean.
Determine the sum of tributary
and channel flows into each
junction.
Find a new junction head based
on the amount of water added to
or subtracted from each
junction.
Perform second step of Runge-
Kutta Integration substituting
values previously calculated as
-------
92
Sequence No. FORTRAN Name
204-207
211-215
216-237
Comments
Compute new values of the cross-
sectional area based on the
increase or decrease of heads
previously computed.
If this is a binary output
cycle, write the cycle number,
junction numbers, heads, the
speeds, and flows in the
channels on unit 10.
If this is a print cycle or the
last cycle, enter the "selective
print routine."
Print Routine
221-222
223-237
241-245
250-258
Convert time to hours and
print a general heading.
For the junctions that are to
be printed, print the junction
heads. For each channel con-
nected to said junctions, print
the channel flows and speeds
after determining the correct
sign.
If the channel speeds during
any cycle exceed a predeter-
mined value print the cycle
and channel number and EXIT.
If this is a restart record
cycle or the last cycle, write
on unit 9. These data can be
used for initial conditions
to another case or used to
start up again in case the run
was interrupted.
-------
Sequence No. FORTRAN Name
267-274
278
93
Comments
Make a copy of the restart
Information on the printer.
Call hydraulic extract
-------
PROGRAM HYDRA	OOOOI
£	##•#~•#»~•~•~»••••»••~•«~~*»••~«*•••##•##~###»•#«~*•*#•	00002
C	~	EXPLICIT SOLUTION FOR	•	00003
c	~ DYNAMIC FLO* In a TWO-DIMFNSICNAL SYSTEM	•	00004
c	•	PACIFIC NORTHWEST WATER LABORATORY	#	00005
C	•	FfOERAL WATER POLLUTION CONTROL ADMINISTRATION •	00006
C	00007
DIMENSION ALPHA(72),Y<5)»YT (5)»AREAS(5).QIN<5) •	00008
•	NCHAN(*5»5> *CLEN(5) tB(5) «ARfA(5) *AREAT(5) •	00009
•	CN<5)•V(5)»VT<5) »Q<5),R{5)• AK f5) ~	00010
 »JPRt <75) »DFEP<5)	00011
COMMON ALPHA»Y»YT.AREA.Q.AREAS.QTN»V,B»CLENtR,CN,OELT»	00012
•	NCHANfNjUNC#JPRTtNj.NC*NCYCtNPRT,MOPRT.PERIOOtNCYCC	00013
REWjNO 10	00014
REW|NO 9		00015
£•••»•*•••••«••••*•••••«•»••«*«•**»«»•«•••#«•*»•»«••*•*•«•»•»«*«*•»*«•	00016
$•	Read* print* and check data	#	oooit
(;#••••~•*•*«»•****•#••*»•»•*«»«**•••*•*«««»##«•«*••*••«•#»••»«««»«»«««	00018
RE AD 5 f 100) (ALPHA < H«I*i«36)	00019
READ(5»105)NJ»NCtNCYCfNPRT#N0PRT,DELT»TZER0.NETFLW»I5TART»lNPSUP	00020
V#RftE(6tU0) (ALPHA(I) ,i«l,36)	00021
READ(Sf530) IPRT»IWRTEtKPNCHI	00022
WnCJ^lliS) NJtNC»NCYC»NPRT.DELT»TZEROtIWRTEtNCYCtKPNcHl»lPRT	00023
IF. .	00026
IST0P»|START-1	00027
IPRT«(ISTyP/NPRT~1)»NPRT	00028
DO 116 j«IWRTEfISTCP	00029
116 »EAD'dO)	00030
C	•~•##•#•###~*~#~»###•#	00031
C	• JUNCTION DATA •	00032
C	00033
118 DC 119 J«1»NJ	00034
*EAD(7»120>JJ»AWEAS(J),(NCHAN(J*K)»K»lt5)»Y(j)»OlN(J),yt(J)	00035
IEjtiSTART.EQ.O)YT(J)«Yf j)	00036
IF(JJ.EQ.J>60 TO 119	00037
W?iTE<6tIl7) JJfJ	00038
STOP	00039
-------
ir{INPSUP.NE,O)GO TO 12I	000*1
**l7pi**l2+>	00042
WTjfi«»I25J	, AREAS (J> tGIN(J) ,  *P(N)«CN(N)»B(m),V(N)	00049
A5EA'
-------
155 v«NCHAN^J,K>	00081
00 160 I«l#2	00082
iFIJ.EO.NJUNCrNtDJeO TO 165	00003
160 CONTINUE	00084
HCXIT»1	00085
*RlTE«6fH5) N,J	00086
|65 CONTINUE	0008 T
. 170 ^ .CCNfi^MC . 			00088
C ###~••••#••##«*•~•#••#•#*»#•~#•#•#«##•«~•#•*»••»##»»#~•#«##*«###*#•	00089
c * wRiTe in?tiav»,0Eometric, ano descriptive data	o* unit 10 » oooqo
C	00091
iriIST*»fvNE«0)9C TO 301	00092
WRlTfllO)fALPHAtl)»I«i»36)»NJ»NC»DELT»(CNCN),R(N)»B tM)*	00093
•	.CUCN»AREAS(J> «QlN(J) ~ (NCHAN(JtKI *K« 1,51 »J*1«NJ> «	00095
•	(AREA(N)»VCN)f(NJUNC(N»I)•Ist«2)tN®ltNC)	00096
XFOC TO 301	0009T
OS 300 N»I»NC	00098
0INJ ¦ AREA(N) • V(N)	00099
3do .continue « - .	00100
**|T£ ,jWRTEt0C TO 51	00112
48 IFfkwRiTE^LE^PNCHI^ISTARTJGO TO 52	00113
KWRlfE ¦ KWRITE - KPNCHI	00114
INK ¦ INK ~ 1	00115
IF ClNI(»t.t«I0)®0 TO 48	00116
WRIT^(6#406>KPNCHI»NCYC	OOllT
NEXIT-I,	00118
51 KWRITE ¦ NCYC	00119
-------
6 m 32.1739	00121
c	••**«#«•»»#**•••»#••»••••»•«*•*•••	00122
£	# COMPUTE CHANNEL CONSTANTS •	00123
C	••••»•••*#**»•*»••»»»•«•»»•**»*»»«	00124
00 19Q N«1#NC	00125
AK • 6 * (CN(H)**2/2»208J96)	00126
|F(nJUNCtNtJ)»LE«NJUNCtN*2)>00 TO 190	00127
K£EP»NJUNC(N»i>	00128
NJUNC(N,1)«NJUNC(N,2)	00129
NJUNC(N*2)*KEEP	00130
196 CCNflNUE	00131
T a TtE^O	00132
!FU$TART.EG.O> JSTAPT*!	00133
£**•*•**••*••*•#*«••«**•«*«*»«**•*«******•**»*•***»»******»***»****»»*» 00134
C* MAfN LCCP * 00135
00136
C	»*•~****•«#*»•»*»*»*	00137
C	• COMPUTATION *	00138
C	**•«#«*•••*«••**«***	00139
00 28§ JCYC-lSTARttNCYC	00140
NCYCC ¦ ICYC	00141
|2 ft* DEUT2	00142
T • Tf.DE^T.	00143
^	•*•*••••••*••••«•*••«••••••••*••«•••••	00144
C	« COMPUTE yA^F CYCLE VELOCITIES •	00145
C	~*••*»••#»*»#~~•••#•####•#***•*•#*#~»•	00146
00 m N*lfNC	00147
ML*NJUNC(N#J)	00148
MH«NJUNC{N«2)	00149
*|N) • ArEA(N) i B(N)	00150
AKT ¦ AK(N) ~ (R(N)**1«333333)	00151
D¥DX « (l«d/R(N)l*('j-AKT *V(N)*ABSF (V )	00154
•	--Y{NL>>>	00155
204 0{N>*VT|N)*AREAjN)	00156
C	•##•*•**•#•#~*•»~*»*~**#*»»##•*»»	00157
C	* COMPUTE HALF,CYCLE HEADS •	00158
C	•#••*****#»##•#•••#•»•****••*••»•	00159
-------
ii;
OC 2£5 J«2»NJ
SUMQaQlN(J)_
00 220 K«lf5
IF(NCHAN(J,K).EQ#0)GC TC 225
N-NCHAN(J.K)
IF0	CONTINUE
225 	Yf(vl>	MPEUT/ARFASCJ))«0.5)»SUM0
c	*•»•#*•»•»•••••••••*•*»••••••«•••»••«•••••••«*•»»••»«•••
C	* . . . Cv*fl?UTE HALF CYCI^E ArEaS, full C^CLE VeLOCITIES •
c
00 230 N«1*NC
NU»NJUNC(N,1)
MH"NJUNp(N»2)
ARE AT (N}»AR£A (N) ~0»5«B f N» * C Yf (NH) -Y(NH) *YT )
ft IN) ¦ AREAT(N) / B (N)
AKT2 » AK(N) / ,(R (N) o*l ,333333)
OVDX ¦ (1*0/R(N))*i({YT(NH)-Y(NH)~YT(NL)-Y(NL))/OELT)
~	(VT(N)/CLEN(N)) * (YT(NH)-YT(Nt)))
V(N)"V(N)*DELT*J(VTJN)#PVDX) -AKT2 *VT(N)«ABSF
-------
c
c
c
* COMPUTE FULL CYCLE AREAS
C
C
C
C
c
c
00 256 N*1»NC
ML«NJUNC(N,U
NH»NJUNC(N,?)
256 	AREAfN) ¦ AREATjN>*0.5«B(N)*(Y(NM)-yt(NH)*Y(NL)-YT(NL))
C*			MAIN LOOP (CpNTiNUED) OUTPUT	»
IF(ICYC#LT»IWRTE)6C t? ?59
* BINARY tape cut •
.	ICYC» (YJ J) » J«1 »NJ) « (V (N> »Q (N) «M»1 »NC)
259 iF(lCYC.NE*IPRT.ANOtlCYC,NE#NCYC>$C TO 263
f. -f.PRINTER.OUT •
IPRT*IPRT*NPRT
f* j/3*6o.O
WRIT516.302) ICYC.TIME
00 3*0 i«l»fiOPRT
J-JPRTCI)
WRITE<6t305l J.YtJ)
00 335 K»i»5 . .
IFIKCHAN(J,K>,EQ*O>0C TO 335
N»NCH*fMJ«JO
1FI^NE,NJUNC(N,I))GC TO 320
VELbVINI
FLCW«Q(N)
00 to 325
320	VEL—VIN)
FLCW«-0(N)
325	WRltE(6t33o') NtVEUFLOW
335	C°NTINUE
340	CONTINUE
C	•' CHECK FOR REASONABLE VELOCITIES *
00201
00202
00203
0020*
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
0021B
00219
00220
00221
00222
00223
00224
00225
00226
0022?
0022B
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
to
-------
263 00 275 N«1»NC	00241
IF(ABSFiV(N)).uf.20#)6C TO 275	00242
WRITE(61,270) TCYCtN	00243
STOP	00244
275	CONTINUE	00245 £
fF(iCYC.NE»NCYC.AND,iqYC.LT#KWRltE)GC TO 285	00246 °
C	•••»•#~»••*#•#•#•#•####«#*	00247
C	* ,. make Restart tape„ •	00248
c	00249
WRIt?(9>120) (J,AREAS(J)', fQlKj(j),	0025o
*	YT(J)00251
WRlff (9tl30)(N,CLENCN)» (NJUNC(N»k) .K«1.2j,R(N)»CN(m)~ R(IM).	00252
~	„ V(Nt»N«1»NC)	00253
KWRiTE«^WRlTE*KPNCHt	00254
ENDFILE 9	00255
REMIMO 9	00256
T2ERC2»T/36d0t	00257
WRITE(6t281) ICYC,TZER02	00258
285 CPNTjNUf	00259
£••••»••••••**•#••««•#••*#•##•*•»••*•••»«*«•••«»»••»«*»*»•••«•«*#«»•»•	00260
C»	END MAin LOOP	•	00261
00262
ENOFILE 10	00263
£	•••»•••*•»•«••«••••••••••«•	00264
c	* PRINT RESTART INFO *	00265
C	•••••••••••*••••••»•••••*••	00266
WR|Tc(6t432)	00267
WRITEI6f402)	00268
WRITE <6#404) (J»Y(J)»AREAS(J)»QIW(J),(NCHAN(J,K),K = 1,5)00269
W?iTEjl6f4l01	00270
WRlTEJ6t4l2)	(NtCLEN(N)*B(N)tAREA(N)«rN(N)*v(N)tR(N)*	00271
•(NJUNC(NfK)»K«1»2)pNml,NC)	00272
WRfTeC6t?99) IWRTEtNCYCC	00273
WRITEt6t422) NCyCC	00274
C	00275
C	• CALL HYDRAULIC EXTRACT PROGRAM •	00276
C	00277
IECNETFLW.NE.O)CALL HYDEX	00278
STOP	00279
-------
c
c
• END-MATN. PROGRAM *	002B1
•••»*•»•*•*#•*#••»•*•••••	00282
100 FORMAT(IBA4)	00283
105 FORMAT<5i5»2F10.0»3l5)	00284
110 FORMAT (1H1///	00285
•	1H l8A4t5X,47H FEDERAL WATER POLLUTION CONTRql ADMINiSTRaT 00286
•ION/	00287
•	1H 18A4»5X«35H PACIFIC NORTHWEST WATER LABORATORY////)	00288
115 FORMAT<132H JUNCTIONS CHANNELS CYCLES OUTPUT INTERVAL TIME 00289
•	INTERVAL INITIAL TIME WRITE BINARY TAPE RESTART INTERVAL	00290
•START PRINT//	00291
•	1H CYCLES»F11.0»5H SEC.tF12.3.14H MRS. CYCLES I4.4H T 00292
~0 I*jIftl9H CYCLES CYCLE 14////)	00293
117 FORMATC4QH0JUNCTION DATA CARD CUT OF SEQUENCE. JJ® I4,4H,J« 14)	00294
120 FbRMATti5»F!o*Of5Xt5I3tF10#5.F10.0»Fl0.5)	00295
124	FORMAT|fH *25X«21H*« JUNCTION DATA •*///)	00296
125	FORMAT <86H JUNCTION INITIAL HEAD SURFACE AREA InPUTI OUTPUT	00297
•	CHANNELS ENTERING JUNCTION//<1H • I6tF15.4*Fl7.0»Fll.2,Ii2,	00298
•	00299
127	FpRMAT<39H0CHANNEL OaTA CARD OUT OF SEQUENCE. NN* I4,4h,n» 14)	00300
128	fORMATcJin/Z/lH ?25X»20H#* CHANNEL OAT A ••///)	00301
150 FORMATjimF8.0»2l3tF6#i;F5»3tF5.0f5XfF10.3)	00302
135 FORMAT* 97H CHANNEL LENGTH WIDTH AREA MANNING VfLOCfT 00303
~V 0Yp RADIUS JUNCTIONS AT ENDS//	00304
C?H l5tFll.0.F8.0»F10.lfF9.3.Fl0.5,Fi3,l.I23,I6))	00305
137 FORMAT mlsi	00306
145 FORMAT(?OHOCOMPATIBILITy CHECK, CHANNEL I4.11H# JUNCTInN 14)	00307
177 FORMAtiftFiO.b*	00308
179 FORMAT(1H///	00309
» lHfl5*«32H##C0EFFlClENTS FOR TIDAL INPUT»»///	00310
•	6X*2HAl»8Xt2HA2;8X,2HA3*8Xt4HPHI2*8Xt4HPHI3«8X«6HPERl0D//	00311
•	5FiOt6,FiO,2///	00312
•	6H.WHE&E//	00313
•	41H YU)« Al*A2#SIN(WT*PHl2)~A3#SIN(WT*PHI3))	00314
270 FQRMAT(34H0VEL0CITY EXCEEDS 20 FPS IN CYCLE l3,10Hf CHANNEL 13,	00315
*23Ht EXECUTION TERMINATED.)	00316
281 FORMAT(4«H0RESTART DECK TAPE WAS LAST WRITTEN AFTER CYcLE 14	00317
••26H T2ER0 FOR RESTARTING « F7.4)	00318
299 FORMAT132H0TAPE 10 WAS WRITTEN FROM CYCLE I6.10H TO CYCLE I*//)	00319
-------
•	27H SYSTEM STATUS AFTER CYCLE„I4tFl2.2#6H HOURS//	003?l
*	54h junction head channel velocity	flow/ 00322
•	§*H N^MbE* (FT) NUMBER tFPS)	(CF$)) 00323
305 F0RMATUH0I5fFl3.4>	003?4
330 FORMAT (1H I28,Fl4.5,Fi2.n	0032*
402 FORMAT (1H1///	00326
•	32H JUNCTION DATA FOR RFSTART DECK///)	0032?
404 FORMAT (B6H JUNCTION INITIAL HEAD SURFACE AREA InPUT-OuTPUT 00328
•	CHANNELS ENTERING JUNCTI0N//(1H tI6.F15.4.F!7.0tFU.2»Il2,	00329
~	4J$)>	00330
406 FORMAT<#0WJTH KpNCHIM.16,# AND NCYC»#»I6#	00331
*	* RESTART RECORDS WILL BE MADE TOO MANY TIMES*)	00332
410 FORMAT (1H1///	00333
*	11H CHANNEL DATA FOR RESTART DECK///)	00334
412 FORMAfi 97H CHANNEL LENGTH WIDTH AREA MANNING	VFLCCfT 00335
•Y HYD RADIUS JUNCTIONS At ENDS//	00336
•UH i5tFil.dtFB#0fF10.1tF9.3»Fl0.5fFl3.2. I23tl6>)	00337
422 FCRMftT142HQEN0 OF TWC-DiMENSlONAt EXPLICIT PROGRAM. I4*8h	CYCLES.) 00338
432 FORMAT(36H0END OF FILF WAS WRITTEN ON TAPE *0.)	00339
530 FORMAT(315)	00340
-------
APPENDIX II
-------
105
APPENDIX II
NOTES FOR SUBROUTINE HYDEX
Unit	Use
3	Output to QUALTEMP (binary
records)
5	Control and System Input
(card images)
6	Standard Output (printer)
10	Input from HYDRA (binary
-------
106
PROGRAM NOTES
Subroutine HYDEX
Sequence No. FORTRAN Name
Comments
Dimensioned variables in addition to those discussed in HYDRA.
7-13	ARMIN, ARMAX
NMIN, NMAX
QEXMIN, QEXMAX
QEXT
QNET
RANGE
VEXT
VMIN, VMAX
YAVE
Minimum and maximum channel
cross-sectional areas over the
entire run.
Cycle when the minimum and
maximum head occurred.
Minimum and maximum quality
cycle average flows over the
entire run; i.e., minimum and
maximum QEXT's.
Accumulates channel flows over
a quality cycle; becomes the
average over the cycle.
Accumulates channel flows over
the entire run; becomes the
average over the run.
(YMAX-YMIN)
Accumulates channel velocities
over quality cycle; becomes the
average over the cycle.
Minimum and maximum channel
velocities over the entire run.
(Note: entire run means
hydraulic cycles used by
HYDEX-NSTART to NSTOP)
Used to accumulate junction
heads over the entire run; be-
comes an average for the entire
-------
107
Sequence No. FORTRAN Name
YMIN, YMAX
YNEW
Variable names as they
19
24
25	NODYN
26	FNODYN
30-33
34	NSTOP
35	NSTART
Comments
Minimum and maximum junction
heads over the entire run.
Updated value for junction
head.
occur in the program.
Rewind unit 10 which contains
information from HYDRA.
Read heading information to be
printed later.
Read the number of hydraulic
cycles per dynamic (water
quality) cycle. Example: If a
water quality cycle of one
hour is to be used and the
integration period in HYDRA
is DELT =120-0, then NODYN =
3600/120 =30.
Floating Point NODYN.
Read system information com-
puted by HYDRA and stored on
unit 10.
Set NSTOP equal to the total
number of cycles in HYDRA
(NCYCC is passed through
COMMON.)
Start HYDEX a specified
number of tidal cycles from
NSTOP.
Example:
PERIOD = 12 hours
¦ 12*3600 seconds
DELT » 120 seconds
NCYCC ¦ 961
NSTART = 961 - 3600 X 12 - 601
-------
108
Sequence No. FORTRAN Name
36-39
37	DELTQ
38
39	JRITE
43 ICYCTF
43
45-58
59-65
69
Comments
This allows convergence to be
achieved in HYDRA before
extracting in HYDEX. This
(converged) cycle can then be
run repeatedly in QUAL for any
number of dynamic steady-state
cycles.
Write the alphanumeric infor-
mation from HYDRA with a
general heading.
Find the quality cycle in hours.
Print information from the
hydraulics program as well as
starting, stopping and
interval cycles used.
The hydraulic cycle number
when the next quality cycle
begins.
Read from unit 10. The
hydraulic cycle number which is
currently being processed.
Read and ignore the hydraulic
output from HYDRA on unit 10
until the hydraulic cycle read
is the same as the starting
cycle in HYDEX.
When the starting hydraulic
cycle is read from the tape,
initialize several variables.
As each hydraulic cycle is read
from the tape, update minima
and maxima, and add to
accumulator variables.
After initializing, branch to
write initial conditions on
-------
109
Sequence No. FORTRAN Name
71-80
81-85
86-94
95-105
110-116
117-122
123-126
127
128-131
132
133
Comments
If the velocity in a channel is
zero, compute area from
junction heads; otherwise, from
Q/V.
Initialize the area variables
if this is the (NSTART + 1)
hydraulic cycle (KFLAG - 1),
If this is the first hydraulic
cycle in the next dynamic
cycle, summarize.
When one quality cycle is
through, complete averages
over the cycle, update minima,
maxima and add to accumulator
variables for entire run.
Adjust the flow and velocity
accumulators to include only
1/2 of the current hydraulic
cycle.
If this is the first quality
cycle (KFLAG2 = 1) initialize
minima and maxima variables.
Otherwise, update maxima and
minima.
Output the flow and velocity
average to unit 3.
Reinitialize the flow and
velocity accumulators.
Skip to the summary portion if
this is the last cycle.
Output the cycle number, and the
initial heads for the next
-------
110
Sequence No. FORTRAN Name
134
Summary
142	FNSMNS
143-153
155-165
166-171
174	K
176-180
Comments
Update JRITE (hydraulic cycle
number at beginning of next
quality cycle).
Section
Floating point representation
of (NSTOP-NSTART).
Compute average flow, area, and
hydraulic radius, range, and
average heads.
Output average flows,
descriptive and geometric
information on unit 3.
Print summary results.
Number of dynamic cycles
processed.
Print a few values from the
ocean end of the estuary for
each quality cycle, to check
-------
e
c
c
c
c
SUBROUTINE. HYOEX
####~####•#~•#•«
•	NET FLOW PROGRAM	•
•	pacific northwest water laboratory	•
•	FfOfRAL WATER PORTION CONTROL ADMINISTRATION •
OIMENSjCn VHIN(5)«VmAX(5),ARMIN(5>»ARmAX(51 *
•	QEXMIN C5)tQEXMAX(5)»YMlN(5)• VMAX(5),RANGE(5)•
•	ARAV£(5)»NMIN{5)f NMAX(5)
DIMENSION AtPHA(72),Y(5>•AREAS(5)«0lN(5),NCHftN<5,5>~
•	V(51*0(5)»AREA(5)#B15)•CLfN(S),R<5>•
•	CN(5),MJUNC(5,2),jPRT(75)»YNEW(5),ONET(5> ,
•	QEXT<51.#VEXT(5),YT<5)*YAVE<5>
COMMON AlPHA»Y»YTfAREA,Q.AREAS.QTN»V,B,CLEN,P,CN,DELT*
•	NCHAN t NJUNC9 «^PRT ? N./t NC t NC YC •NPRT,N$PR T ~ PERIOD»NC YCC
C*	JN?UT ANO INITIALIZATION	*
REWIND 10
*EWIN0 3 ...
c
G	*. . . CWRO^^ARD$ .UNfQUE fp,HYOEX •
C	••••»••••••••••*•••««•••••»»»•*•»•••••
READ(5*103)(ALPHA(I),1-37,72)
READ(5*BO) nooyn
FN5DYN«FLpAfFJNOOYN)
c
£	~ _SYSTEM. INFORMATION FROM UNIT lo •
C
REAOJX01	CALPHA(1)«I«1,36)«NJ*NC«DELT•(CN(NI»R(N)*®(N)•
•	C^EN(N),N«l,NC)
RE AO j( 10)	CYC J) »AR£gS(J) tQIN (J) , (NCHANC J»K) »K«l »5) , J»1 »Nj> ,
•	(AREA CN)»V(N) t(NJUNC(Ntl),I«1,2)tN«l»NC)
N$TOP.».NCYCC
NSTART ¦ NCYCC - (PERIOD / OELT)
WRlTE<6»f05j| CALPHA(p tf«l»72)
OEUT P«OELT*fNCOYN/3$00•0
WRITEC6»351) NSTART,NSTCP#OELT.NCOYN,OELTQ
JRfTf ¦ NSTART
00001
00002
00003
00004
00005
0000*
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
OOO^o
00031
00032
00033
00034
00035
00036
00037
00038
00039
-------
C* 	 MAJN LOOP	* 000*1
£•••~#~#••#*#••*#••#•••#•*••••#•**••*•*••*##•»•*#•*#»•»##»»•»~#~••»••»	00042
202 R^AOdOl ICYCTF, (YNEW( J) , (V(N) fO tNl fN«l,NCl	000*3
IF(ICYCTF- NSTaRT)20?;264»z68	0004*
204 05 206 N*1 »NC . ,	00045
C	00046 ^
C * PROCESS FIRST HYDRAULIC CYCLE *	00047
C ~ ?FOR INITIALIZATION)	• 00048
C	00049
QNEf(N) ¦ 0,5*Q(N)	00050
QEXT(N) « 0.5*Q ¦ 0.0	00060
YMiN
-------
156
1ST
00
21q
IF(KEL*&,VE*1)9C TO 1ST
OS 156 N*lpNC
¦ A*EA/N)
§PHAX(H> » AREA(N)
continue
2JO N»1»NC
ONET(N) ¦ QNET(N) ~
OfXfjN) » Q£XT(N) ~
V£XT * YNEW(J)
NHA*00 TO 180
YM|NtJ> • YHEH(j>
yMlNJJJ ¦ ICYCTF
C3Hf{NUE
IF ¦
V£XT|N> ¦
continue
IF(KFLA6?.NE.1)6C TO
00 la* N**»NC
OEXMiN(N) « °EXT(N1
OEXMAX(N) ¦ OEXTIN)
OEXT(N) • 0.5«0
-------
181 CONTINUE	00121
60 TO Is$	00122
183 DC 1|7 M»1»NC	00123
IF(OEXT (N)•GT«Q£XMAX(N)JQEXMaX(N> «GEXT 	0012*
!F(OfXT(N)tLT.QEXMiN(N))QEXMtN(N)«QEXT(N)	00125
187	CONTINUE	00126
188	WRitf(3> *J»l.NJ>	00133
J*lTp ¦ JRITE ~ NOOYN	00134
	OO.Ty . ?02		00135
00U6
		END MAIN tOOP	* 00137
C	••••»•*•«»*•»••»•••»••«•*»**•	00139
C	* §UMMARfZE ALU CYCLES *	00140
C	••••••••»••»••***•*»*•••••«••	00141
220 FNSMNS«fl,CATF (NSTGPwNSTART)	00142
00 222 N»1»NC	00143
QNEfJN> ¦ QNET(N) - 0.5»Q(N)	00144
QNETCN) « QNET(N)/FNSMNS	00145
ARAVE ¦ YAVE(J) / FNSMNS	00151
Y(J) ¦ 0.0	00152
260 CONIINUE	00153
REWIND 10	00154
WR|fE(3)(ONET(N)»N«IfNC)	00155
WRITE(3? 
-------
*Re4(HfmAfrE4(Nt + (B(N)/2.0r*((Y + (Y(J)-YNFW(NLi)1	00161
2*6 CONTINUE ..	0026?
W*lT£ C3) (Y(J),AK£A5(J),QlNtJ),(NCHANfJiK)pK*l,5),	00163
•	«	00166
•	VMAg(N)•ARMIN(N),ARMAX(N)•ARAVE*M«1»NC)	00167
REWIND 3	00168
WRIT£'<6>262>  ICYCTF# (YNEW(J)t JaltNJ)	00177
REA0I3) (QEXT(N)tVEXTtN)•N*1*NC>	00178
^Plfr^t232) ICYCTF, YNEW(l) ,OEXf (1) tQEXf (2)	00179
234 .CONTINUE	00180
RfWlND 3	00181
WTE<6t240)	00182
RETURN	00183
C	00184
C •	fND fNfIRE,SUBROUTINE •	00185
C •••••*•»»•»•*••••#•••••«•••«••	00186
80 FORNiffSfSl	00187
J03 FORMAT<%8A4)	00188
105 FORMAT 4}M1///	00189
» 1H 18A4,5X»47h FEDERAL WATER POLLUTION CONTROL AoMlNiSTRAT 00190
•ION/	00191
•	I% 18A4.§X#37H PACIFIC NORTHWEST WATER LABORATORY /	00192
V. 1H 18A4/JH 18A4////J	00193
224 FORMAT J119H •<•••• FLOW	»•••<• 00194
•	• • VELOCITY • ~	* • • CROSS-SECTIONAL AREA * # •/ 00195
•	II8H CHANNEL NET FLOW MIN.	MAX. 00196
•	MIN. MAX. MIN. MAX*	AvE./ 00197
•	, _U9H NUMBER  (CFS)	(CFS) 00198
•	_.CCFS> (CFS)„ 
-------
?32 FORMAT(I7»5X»F10.2»6X.F11#2*FI2.2>
240 FORMAT(25H0EN0 OF NET FLOW PROGRAM,I
2*? FORMAT(lHl///
•	53H	OUTPUT FOR CHECKING DATA ON EXTRACTED TAPE •#*«///
•	AW hydraulic head AT	»flcw in channel#/
•	49H CYCLE JUNCTION NO*1	NO*1	NO.2//)
262 FORMAT CIHi////
•	98H JUNCTION MINIMUM HEAD OCCURS AT MAXIMUM HEAD CCCU
»RS AT AVERAGE HEAD TIDAL RANGE/
•	94H NUMBER	(FT)	CYCLE		CY
•CLE (FT)//
•	<1H i^tF15.2»n3fF16,2tI13.FI6f2.Fl5.2))
351 FORMAT(B8H •~•*•••• FROM HYDRAULICS PROGRAM #«~»«•*»
•	CYCLES.0ER TIME INTERVAL in/
•87h start cycle stop cycle time interval
•	. quALItY P*OGrAM'/ -
•1H I7»!l4*Fll#0f9H SECONDS.10X«I6»12X,F9»2.7H HOURS/////)
END
HYDRAULIC
Quality cycle
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
-------
APPENDIX III
-------
APPENDIX III
119
NOTES FOR PROGRAM QUAITEMP
Unit	Use
2	Control Input (except unit 11)
(card images)
3	Input from HYDEX (binary
records)
9	Restart output (card images)
10	Output for extracting program
(binary records)
11	Waste flow input (card images)
16	Standard Output (printer)
61	Standard Output (printer or
-------
120
PROGRAM NOTES
QUALTEMP
Sequence No. FORTRAN Name
Comments
Dimensioned variables in addition to those discussed in HYDRA
and HYDEX
7-18
ALPH
ASUR
BETA
C
CIN
CLIMIT
CMASS
CONST
Intercept (millibars) used in
temperature, vapor-pressure
approximation.
Junction surface area.
Proportionality coefficient
(MB° C"1) used in the linear
approximation of the tempera-
ture, vapor-pressure relation.
(See Edinger, et al., 1965).
Initial (or present) concentra-
tion at a junction.
Concentration of the ocean
input water. CIN(M,K), for
example, is the concentration
in the ocean for constituent
M,K quality cycles into a tidal
cycle.
Upper concentration limit for
a constituent. If exceeded
during computation, execution
is terminated.
Mass of a constituent in a
junction.
Constant mass of pure constit-
uent added to a diverted flow,
which appears at the junction
to which the diverted flow is
-------
121
Sequence No. FORTRAN Name
CSAT
CSPEC
DECAY
DIFFK
EQTEM
FACTR
JDIV1
JDIV2
JRET1
Comments
Saturation concentration of a
(dissolved gas) constituent. If
the concentration ever exceeds
this value, the concentration
is forced to the saturation
value and an error message
printed.
Concentration of a waste dis-
charge. This differs from
CONST in that it is in MG/L
and dependent on the diverted
flow rate.
Decay coefficient (K,) in SEC"1
for BOD or other substance with
a reaction rate. Base e.
Diffusion coefficient, computed
from CDIFFK and channel
dimensions.
The equilibrium temperature at
a junction.
Multiplication factor applied
to the concentrations to
accelerate convergency. (See
NJSTOP, for example)
Junction number where a
diversion is to occur.
Junction number where a
diversion is to occur.
Junction number to which the
diversion from JDIV1 is re-
-------
122
Sequence No. FORTRAN Name
JRET2
MARK
NCQNDK
NCONOX
NGROUP
NJSTRT
NJSTOP
Comments
Junction number to which the
diversion from JDIV2 is
returned. (There may be up to
NUNITS of the JDIV1-JRET1,
JDIV2-JRET2 pairs.)
Contains the quality cycle
numbers which bound a series
of KOUNTT quality cycles;
used in keeping track of the
binary output for the quality
extraction program.
Contains the nonconservative
numbers. For example, NCONDK(I)
is the constituent number for
the Ith nonconservative con-
stituent; a decay rate is
associated with each such
constituent.
Contains constituent numbers
with an associated reoxygena-
tion rate. The value of
NCONOX(K) (^0) is the constit-
uent number which is paired with
constituent number NCONDK(K).
Number of groups of junctions
for each constituent. For
example, NGROUP (K) is the
number of groups of junctions
to which multiplication factors
are applied for the Kth con-
stituent. Each group consists
of junctions in the series of
to
, inclusive. For example,
suppose that for constituent
number 3, a convergency factor
of 1.5 is to be applied to
junctions 1-4, and that a con-
vergence factor of 2.5 is to be
-------
123
Sequence No. FORTRAN Name
ODECAY
QC
QE
QINWQ
QNET
QTOT
QW
Comments
N6R0UPS (3) = 2,
NJSTRT (3,1) =	1, and NJSTOP
(3.1)	= 4, and	FACTR (3,1) =
1.5.
NJSTRT (3,2) =	7, and NJSTOP
(3.2)	¦ 8, and	FACTR (3,1) =
2.5.
1.0 - DECAY (SEC"1)
Convective heat exchange at a
junction.
Evaporative heat exchange at a
junction.
The flow into a waste producing
entity. There is a QINWQ for
each junction. It is zero if
no water is being removed from
the junction, or if water is
returned to the junction only
after having been removed from
another junction. If it is
negative, waste is being added
to the junction from an
external source which obtained
the water from outside the
system. If it is positive, water
is being removed from the
junction, and may or may not
be returned to another junction.
The net channel flow during a
quality cycle.
The total heat exchanged at a
junction.
The back radiation from a
-------
124
Sequence No. FORTRAN Name
REOXK
RETRNF
VOL
VOLQIN
Undimensioned variable names
discussed in HYDRA, HYDEX.
20-21	A
AP
BB
QRNET
TA
Variable names listed as
36	NOJ
Coircnents
Reoxygenation rate (K2) in
SEC"1 for dissolved oxygen.
Base e.
Proportion of constituent that
is returned to a junction after
a diversion.
Volume of a junction.
Volume of wastewater removed
during a quality cycle
(QINWQO * DELTQ).
in COMMON in addition to those
Coefficient in evaporation
equation.
Air pressure, millibars. Inter-
polated as necessary in sub-
routine METDTA.
Coefficient in evaporation
equation.
Net incoming radiation (KCAL
M"2SEC"1) calculated for each
cycle in subroutine METDTA.
Dry bulb temperature, °C,
calculated in METDTA.
occur in the program.
The number of junctions that
are at the ocean end of the
model. If NOJ is 2, then
junctions 1 and 2 are assumed
-------
125
Sequence No. FORTRAN Name
ITEMP
IEQTEM
37	NRSTRT
INCYC
NQCYC
NOEXT
CDIFFK
NT AG
38	IPRT
NQPRT
Comments
The number of the constituent
that is temperature; i.e., the
ITEMPth constituent is tempera-
ture. (=0 if temperature is not
a constituent).
Switch which, when non-zero,
indicates that the equilibrium
temperature should be computed.
The first hydraulic cycle to be
used from the extract type.
The first quality cycle which
is to be processed by QUAL.
The last quality cycle which is
to be processed by QUAL.
Switch which can be used to
indicate that EXQUA should be
called.
Constant used in computing the
eddy diffusion coefficient.
Counter to indicate how many
quality cycles have passed in
one tidal cycle. Initially,
it describes where in the tidal
cycle the program will start.
NTAG runs up to NSPEC, and is
then reset to zero.
Holds the next quality cycle
which will generate printed
output. (See Figure 12)
Print output every NQPRTth
cycle.
NEXTPR
Output is printed every
NQPRTth cycle for one tidal
-------
126
Sequence No. FORTRAN Name
INTBIG
IWRITE
NEXTWR
IWRINT
39	N0JP1
40-53	K
54-55
Comments
quality cycles are skipped, and
output is then printed again at
the NEXTPRth cycle, every
NQPRTth cycle for one tidal
cycle, etc. See Figure 12)
Analagous to IPRT for the
binary output.
Analagous to NEXTPR for the
binary output.
Analagous to INTBIG for the
binary output.
Number of ocean junctions plus
one.
The number of quality cycles on
the extracted tape, usually a
tidal cycle. Bypass all of the
extracted hydraulic infor-
mation on tape to obtain the
net flow and system information
from HYDEX.
Read additional alphameric
information from cards, and
print the aggregate alphameric
information.
56	DELTQ	Length of a quality cycle,
seconds.
57	DELTQ1	Length of a quality cycle, hours.
58	DELTQ2	Length of the printing interval,
in hours.
59-61	Print constants, counters,
flags, for the run.
62
NUMCON
The number of constituents in
-------
10
20
40
100 quality cycle number
30
50
60
80
90
70
1.	IPRT = 10. Begin printing at the 10th quality cycle.
2.	NQPRT = 2. Print every 2nd cycle (thus, the 10th, 12th, 14th, etc.), until...
3.	a tidal cycle has elapsed, then
4.	start printing again at NEXTPR after
5.	INTBIG cycles have elapsed since beginning of last print cycle.
-------
128
Sequence No. FORTRAN Name
64-66	NALPHA
67	CLIMIT
68-80
84-94	NUNITS
98-111
98	NFIRST
116
118
119
Comments
The number of alphameric
variables to be read which
depend on the number of con-
stituents. Read and print a
card descriptive of each con-
stituent.
Read the limiting concentrations
for each constituent.
Read and print the reoxygenation
and decay coefficients. If
there are none (exit from loop
with k=l) print a message.
Read the number of diversion
return combinations. Then read
and print the diversion-return
information. If no diversion-
return information, so print.
Input waste flow information.
Since all of it won't fit on a
card for five constituents,
read two sets of cards if there
are more than three constituents.
Check each card for proper
sequencing (JJ=J). If a
sequence error is noted, stop.
Is used in reading the two sets
of cards.
Print a "multiplication factor"
heading. Then, for each con-
stituent,
1)	read a card containing the
number of groups for that
constituent.
2)	if there are no groups for
that constituent, go to the
-------
129
Sequence No. FORTRAN Name
121
123-129
131
136-139
143	NSPEC
148	NOPRT
149	OPRT
150-163
164	METDTA
168	NCOUNT
Coironents
3)	for each group, read infor-
mation describing the group.
4)	with that information, apply
the return factors to the
starting concentrations.
Print "no concentration factors
applied" for constituents for
which NGROUP is zero.
Print the wasteflows, and
adjusted concentrations.
The number of times concentra-
tions are input at the ocean
end. There should be enough
for one tidal cycle at the
appropriate interval. For
example, if the tidal cycle is
25 hours and the quality cycle
is 20 minutes, NSPEC=75. Read
cards containing NSPEC con-
centrations.
The number of junctions for
which printout is desired.
Contains the junction numbers
of NOPRT junctions.
Print the channel and junction
geometry.
The entry point in the
meteorological data subroutine
which inputs the meteorological
data.
Counts the number of times data
is printed in a tidal cycle, so
that when the tidal cycle is
over, NEXTPR may be used to
-------
130
Sequence No. FORTRAN Name
169	KOUNTT
170	NOB
171	NEOT
172-177
181-194	KVOL
196-198
200-210
211-214
218-222
Comments
Used in controlling the binary
output notes in MARK.
The number of binary output
tidal cycles elapsed.
ICYCTF at end of extracted
tape.
Reorder the junctions connected
to a channel, if necessary.
A flag to indicate that the
loop following has been gone
through twice. The junction
volumes are computed twice.
They are computed from CLEN, B,
and R, where the first time, R
is an average from all of the
quality cycles. The volume
thus calculated is an average
volume.
The average volume, and other
initial and descriptive para-
meters are output as initial
information to the quality
binary output tape.
A new R is computed, which will
make the volume correct for a
specific junction head, Y.
The new junction volume is then
computed from the R's when this
is done, skip to 774.
The heads which started the
hydraulics extract tape are
input, and the volumes
corrected for these heads.
Initial mass concentrations
are computed from the initial
-------
131
Sequence No. FORTRAN Name
226-227
228-230
231-239
244	NQCYCC
248
249
250
255	VOLFLW
254-272	FACTOR
Comments
Eddy diffusion constants are
computed from channel geometry.
Waste water volumes which
change during a quality cycle
are computed.
If binary output is to be made
from the first cycle, write the
initial concentrations.
Used to retain the value of
ICYC when the main loop is
complete.
Channel flows and volumes are
obtained from the extracted
information.
If all of the information from
the hydraulics program has
been read, rewind it, so that
it can be read again to con-
tinue the quality computations
for an indefinite time with
the same basic hydraulic infor-
mation.
Read junction heads from
extracted information.
The volume of flow during a
quality cycle.
Depending on whether the
channel is connected to the
ocean, and the direction of
flow in the channel, the
quarter point solution tech-
nique is applied to the channel
concentration gradient. For a
discussion, see Orlob, et al.
-------
132
Sequence No. FORTRAN Name
273-274
280-291
Comments
Adjust the mass concentrations
for advection and diffusion.
Adjust the mass concentrations
for decay and reoxygenation.
If NC0NDK(I)=0, the Ith con-
stituent is conservative, and
no correction is made for that
constituent. If it is non-zero,
a correction is made, and
NCONOX(K) checked. If non-
zero, a reoxygenation correction
is made based on the constit-
uent.
295-296	NTAG	The ocean concentrations are
input for one tidal cycle, at
each quality cycle. If the
tidal cycle is complete, reset
NTAG, so that the ocean con-
centration information can be
used again.
297-300	Set the concentrations at the
ocean junctions to the ocean
concentration.
305-310	Adjust for waste flows. If the
waste flow is negative, an in-
flowing waste from outside the
system is assumed, and the mass
concentration at the junction
is adjusted using the volume
flow and CSPEC, the concentra-
tion of the waste.
311-314	If the waste flow is positive,
an outflow is assumed, and the
mass concentration is adjusted
using the concentration at the
junction.
318	If NUNITS is zero, bypass the
loop for adjusting concentrations
-------
133
Sequence No. FORTRAN Name
319-330
332
335	FMD
336
337	RHOW
338	TWC
339	HV
340	ROXDR
341	EA
Comments
For each diversion, adjust the
mass concentration of the
receiving junction on the basis
of the volume flow from the
contributing junction, the
original concentration in that
junction, the return factor
(RETRNF), and CONST, (which
allow for pollutant which enters
during the diversion.)
If temperature isn't being com-
puted, skip the next loop.
The subroutine which Fetches
the Meteorological Data
according to ICYC.
DO for each junction where CIN
isn't an input, i.e., where
the concentration isn't fixed
by the junction being in the
ocean.
Water density in KGM"3.
Initialize surface temperature
(or update it).
Compute latent heat of
vaporization.
The reciprocal of RHOW*DEPTH
with a conversion factor for
ft to meters.
Saturation vapor pressure (MB)
at the wet bulb temperature of
the air.
343
ES
Saturation vapor pressure (MB)
at the temperature of the
-------
134
Sequence No.	FORTRAN Name
344	DELVAP
345
346	T1
349	DELTMP
350	B0WM0D
353	QDEP
354	QR
356
357
358-359
360
361, 362 T2, T3
Comments
Difference of the above two.
Since evaporation can still go
on at (reported) wind speeds
of 0 MSEC'1 wind is set to
0.05 MSEC-1 even if it is a
calm day.
Temporary variable used in
computing QE, QC.
Difference in air-water
temperature.
Modified Bowen ratio.
Sum of the terms dependent
on surface water temperature.
KCAL M-2SEC"1.
Atmospheric radiation terms
(measured or computed) which
are independent of surface
temperature.
Computed temperature gain or
loss since initial condition,
or updated from the last com-
puted value.
Bypass, 1f not Interested 1n
the equilibrium temperature.
The limits of the table used 1n
computing EEKTEMP are 0-30°C.
Find out where the initial (or
last computed) temperature
lies within the table.
Temporary variables, using
ALPH, BETA to compute the heat
exchange coefficient and the
-------
135
Sequence No. FORTRAN Name
363	XCHCF
366	DTEM
367
368-371
TEM2
376-381
385-392
394-408
Comments
The exchange coefficient.
Difference between last
temperature and equilibrium
temperature.
Temperature calculated by the
equilibrium temperature
equation.
See if TEM2 (above) uses values
of ALPH, BETA originally
assumed. If not, use the newly
computed value of temperature
to obtain new values and go
through the loop again.
Compute new junction volumes,
and from them, new specific
concentrations.
If negative concentrations
occur, set the concentrations
to zero. This condition can
occur if the time step is too
large. This is one form of
instability and this corrective
procedure doesn't really cure
the instability, although it
may be partly justified if the
concentrations are low. If this
is a print cycle, an error
message is printed.
Two checks are made on the high
side of concentration. If the
constituent is nonconservative,
a check is made that the paired
constituent does not exceed its
saturation value. If it does,
it is set to the saturation
-------
136
Sequence No. FORTRAN Name
409-418
419-438
423	KOUNTT
439-458
426-477
Comments
All constituents are checked
against an upper limit (CLIMIT).
If they are higher than that
limit, execution is terminated.
Binary output is made for the
quality extracting program
according to a rather confusing
sequence of counters and flags.
Basically, these are arranged
to give output every quality
cycle for one tidal cycle,
and then to skip several
tidal cycles, and output again
later in the program.
Accumulates the number of times
in a tidal cycle that output
has occurred. When KOUNTT =
NSPEC, the last output is
written, and output is
suppressed until NEXTWR, which
is computed from when the last
series of output began.
A restart card deck containing
the non-constant variables is
made if this is the last cycle.
These could be used to continue
the program if the values have
not stabilized in the length of
run selected initially. Because
all of the information will not
fit on a card, two loops are
necessary if there are more than
three constituents.
Printed output is now made
according to a series of
counters and flags described
-------
137
Sequence No. FORTRAN Name
478-489
Comments
If the equilibrium temperature
is to be calculated, it is
printed in addition to the heat
budget terms which are listed
in both KCAL/M2-sec, and
BTU/hour.
497-499
620-END	Meteorological
A subroutine to extract the
quality data in a form somewhat
similar to HYDEX has been used
elsewhere. However, it is
rather specific to a given
locale and not particularly
useful in the case of the
Columbia River and is not
described in this report.
Data Subroutine
629
632
635
INT
NPTS
NQCSM
QRNETAO
UWINDAO
TAA()
TAW()
APA()
IDQ
The interval, in seconds,
between data points on input.
The number of data points input
(should be enough for one day).
Time, expressed as the number
of quality cycles since mid-
night at the start of the
quality program.
An array of QRNET's (net
radiation).
An array of UWIND's (wind
speed).
An array of TA's (air
temperature).
An array of TAW's (wet bulb
temperature).
An array of AP's (air
pressure).
The integer representation of
the length of the quality
-------
138
Sequence No. FORTRAN Name
636	FINT
637	LOT
639	FMD
646	ITIM
647
648	I
Comments
The floating representation of
the interval in seconds
between meteorological data
points.
The length of the meteorological
table in seconds.
"Fetch meteorological data."
Called when the table is to be
referenced at a particular
time. The values are inter-
polated from the arrays, and
stored in COMMON in variables
QRNET, UWIND, TA, TAW, and AP.
The seconds of elapsed time in
the simulation since the start
of the program.
The entry in the meteorological
table which immediately pre-
cedes the present time.
ITIT	The seconds of elapsed time
since the start of the last
meteorological data set.
649
FACT
A factor for interpolating
between the Ith and (I+l)th
-------
. .... P9O0BAM. OUALTEMP		0000J
(;*#*•*#•**#**••*•«***•**••#»*»•##•*••«##*»•*•«*•#»•*»••••*#**••«••#*«•	00002
C*	QUALITY PROGRAM QUARTS POINT VERSION	*	00003
C»	PACIFIC NORTHWEST WATfR LABORATORY	~	00004
C* . .Ff^pERAL . WATER PORTION CONTROL ADMINISTRATION	#	00005
OJMENSiON 0ECAY(5> ,RE0xk<5>•NCCn6k<5),NC0n0X(5),CSaT<5i »CDEC&Y(5i	OOOOT
OI*£N§iQN N6RCUP(10)»FACTRC5«I0).NJSTRT(5»10)*NJSTOP(5,l0>•	00008
•	MARK(10.2)	00009
OIMENSjpN jpiVl <|20) tJ0iV2(20> • JftfTHSO) »JRET2<20> »R£TRhF (20»5) ,	00010
•	„ CONST(20*51	00011
0IMENSJ3N YNEW(5) »V0LQiN(5),Q(5,5)»CSPEC<5.5),	00012
•	, CIN(5flOO) '|VOL(5) «ASUR(5I •QfNWQ(5> .QNET'( 5) ~	00013
•	CMASSI5*5)»DIFFK( 5) t. ALPHA<198)»CLIMIT(5)	00014
OTMEN$I$N Y(5)•AREAS(5)»QIN(5)«NCHAN(5,5)•	00015
•	V( 5>«0( 5)»AREA( 5)»B( 5>.CLEN( 5>tR< 5).	00016
•	_ CNJ 5>»NJUNC( 5f2)~JPRT( 75),DEEP(5>	0001T
DIMENSION QC(5)fQW(5),o^(5)#eqTem(5).QTcTtS)	00018
COMMON ALPHA• C»CSPEC»VOL• QINWQtNSPEC• DELTQ* NllMCON»NALPHA,NJ11CYC~	00019
•	NOOYNfNSTART»NSTOP,ASUR#MARK»NOB#lTEMPfIEQTEM»QRNET»	00020
•	UW|NDfTAtTAW«AP»A«BB	00021
CCMMON/OATA/ALPh(6).BETA(6)	00022
0AfA}ALPH»5#7»4,0».75T,.5.4i,-15.29,.30.43)•	00023
•	(BETA«0.62« 0» 842•U107•I•459.1.898•2.449)	00024
EQUIVALENCE (AREAS,ASUR)t(QlNtQTNWQ),(CNtOIFFK)*(VOLQlNtONET)	00025
•»(CMASStNCHAN)•(YN£w»AREA)	00026
REWIND |	00027
REWIND 10	00028
Dp 38 NQB-1.10	00029
0Q-3* K«l,2	00030
36 MA«KiNOB,K) « 0	00031
3B CONTINUE 		00032
£	•*••*••*••«••••••••••••«*•••••••	00033
?	• . iNPuf.CpNfROL CONSTANTS *	00034
C	-	00035
*EA0(2.80) NJtNC*NSTARf;NSTOP#NOnYNtNOJ#ITEMP,IEQTEM	00036
RfAD(2tBi) NRSTRT«INCYCfNQCYCtNOFXT•COIFFKtNTA6	00037
-------
N8JP1»N9/*1-	00040
c	•«••••«••••••••#*••«•••«#«•«•*«#«•«#•••••«•••#••	00041
c	* FiNO SYSTEM INFORMATION FROM HYOEX UNIT •	00042
C	~#*••#*••*•••~•~•~•~••••••#•~~••••••••••••••••••	00043
K » |NSTpP-NSTART)/NODYN	00039
0p^8$ J * 1#K	00044
RCA0(3J	00045
R£APi3?	00046
96 CONTINUE	00047
RE*D<3> (QN^T(NtfN*IfNC)	00048
READ(3>(ALPHA(I)»I«lt36>tNJiNCtOELTt(CN(N)#R(N>.B(N),	00049
#	, CLEWN>.»N*1»NC) _	00050
REA0C3)IY f J)»AREAS IJ)»0IN » (NJUNC(N»i) »I*Tt2>	00052
*E*W 3	00053
RCAp/?»103)fALPHA(I)fI«37»7?l	00054
WRITE<«!«i05> (ALPHA(D,1*1,72)	00055
6eL™«DELT#FLCATF	00056
DELTgJ*6eLTQ/3660.0	00057
DCtTp2«§EtTQl»FL0ATF(N0PRT)	00058
WRjTE.(ftlf 106) N5TARTfN5TCP»DfLT	00059
WRITE(l6»107)NRsTRTfINCYC»N0CYC,INTBl8»OELT02,DELTQl#CDlFFK	00060
WRff^IALPHA(I;,i«X09?NALPHA)	00066
REA0(2»110) {CLIMIT(KJ,K«1#NUMCCN)	00067
C •••••••*••••••«««••••«•«•«•*«*•*»••«***•##«*»•*«•#•••»•••«»*	000^8
C • INPUT/OUTPUT—REOXyGENaTICn AND decay COEFFICIENTS *	00069
C «*••*«••«••*«•«•»*•••««*••*•*•«*•«•»*•«#*•«*•*•••«*•«•«*••»*	00070
REA0(2»*d* {NCONOK(KKNCONOxtK) ,K«1	00075
IF(NfONOX(K).EQ,0)WRITFC16»58>NCCN0K(K)•OECAY(K)	00076
*F #DECAy ~	00077
~ NCCNCX(K)«REQxMK)»C5AT(K)	00078
44 C^TIHUE	00079
-------
C	• ,jNpUT/QUTPuf«-pxvERSiON-BETUBN FACtC*S	•
ReAD(2fU2) NUNITS
ir
OC 11T I«l.NUNITS
RCAbf2«il6)JDXyf(I)9JDfv2(Z>«JREti(li»JRET2t
~	(RETRNF fCONST(f#M) »M«1 •NUMCON)
WRfTEfl«»350) ItJOiyl (T) »J0JV2(I) tJRET1 (I) »JRET2(I) t
•	(RETRNF(I.H)tCCNST(I»M),M«1,NUMCON)
117 CCNTfNMt	
c	*•••*••»••*••••••••••*•••«*•••**•«•••••••»••»••»•*»•»
C	• ..fMfUT«-WASTE FLOW CONCENTRATIONS AND AMOUNTS *
C	»•~•»#••»»»•••••••••#~»#••##•*••*#»•»#•~#######•#••#»
UB NFJR$T?3
IF(NUMCSN.LT.3)NFIRST«NUMCCN
00	. . 		_
READ(11.200) JJtOlNWQ(JI»CC(J»K),CSPEC(J»K)»K-1*NFIRST)
- - IF|^.C0«J»60 TO 206
202 WRITE(16»204)JJ,J
STOP
206 CONTfNUjr
if INU(JC0«.LE.3)60 TO 212
NFIRST«NF1RST*1
Ofi	_
R^AD^t1,200) JJt (C(J»K).CSPEC(J,K)•K«NF1RST,NUMCCN>
TO 202
210 CONTjNUf
j INPUT—MULTiPL|CAT|ONHFACTCRS TC BE APPLIED j
21? WRftEl16*224)
00 Z%Z ItltNUHCCN
READ(2.112) NGROUP(i)
ir(NGRO|)P(l).EO.O)GC TO 222
N9 « N8R0l,P(I>
000*1
000B2
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
001 OB
00109
OOUO
00111
00112
00113
00114
00115
00116
00117
ooiie
00119
-------
c
c
c
READ(?f 220) (FACTRC11K),NJSTRT fJ#K)*NJSTOPIX«K>«K«1tM6)
WRlTEU^f?28) I# (K»FACtP(I»K) •NJ^TRT(I»K) ~NJSTCP(!»K) »KaltNR)
OC 234 K«J.N6
Njf m htjSTRf(ltK)
¦ M^STCP(ItK)
Dp 23* j«NJlfNJ?
C(J.I) • C(J»I) • FACTR(I»K)
234 CONTINUE
222 CCNT|NUf
DO 232 J«1»NUMCCN
TF.EQ.0)tfRITE(l6,?i6)i
232 CONTINUE 	
# CUTPUTr-ApUUSfED CONCENTRATION? *
WRIT? {iftt 241)
OC 283-J»l»NJ
WRlTF(l6t2fl2> J»0INWQ(J)t(C(J«K)•CSPFC NSPEC
oc i86_M*i*NUMCcN _ _
RCA0(U,l84) CCIN(M,I) ~I«l,NSPfC)
WRfTg<6Jtl88)M» »K«1«2)»Nt Y (N)•(NCHAN(N»I)•Tsl• *>«N»l«Nl)
Nj * Ml ~ 1
tFtNj - MC)76»79»78
78 WRITEt\95)tJtY(J)»(MCHANlJ»K).K*l,5>,J*N1.N2>
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
4*
-------
GO Tjj 79	00161
76 WRITE(61t194)(N»CLEN(N)«B(N)•AREA(N),CN(N)«QNET(N)•	00162
•	RCNi#(yjUNCCALL METOTA	00164
00165
C*	INITIALIZATION	• 00166
00167
NCCUNT ¦ 0	00168
KCUNff * 0	00169
NCS » 0	00170
NEST ¦ N5TCP - NODYN	00171
DC 358 N«1»NC	00172
IF(N^UN£(Nf1>•LE»NJUNC(N*2))60 TO 358	00173
KEEP»NjjjNC»NjUNCIN»2>	00175
N,jU#C$UM/ASUM „	00192
VpL(J) ¦ ASUR(J)«OBAR	00193
373 CONTINUE	00194
IFacVPLtNE.djGC TC 774	00195
WRjTE(fg) (ALPHA(I),I«1.NALPHA)	00196
WRltE(lO) N^tN06YNtNSPEC*0EUTQ«<0!NW0{J)#VCL(J>»ASljR{J)•	00197
•	CCSPECfJ»K)tK«l,NUHCON)00198
KVCt » f	00199
00 710 NaltNC	00200
4k
-------
c
c
c
R(N) • AREA(NI / 8
710 C5NTINUC
60 TC 359
• CORRECT VOLUMES FOR STARTING HEADS *
774 READI3)	ICYCTF#*YNFW*J)»J«1*NJ>
ircicrCTF.eC.NRSTRTjeo T0_776
READ(3)	(0(N>tV(N),N«1»NC)
as to
776 r>0 780 j*l»NJ
VOMJ) « VOt>30»32»34
30 IWRITE ¦ INCYC
90 To 3$
3? wpitfiip)	IWRlte»ttCU«K>tK«itNUMCCN)
NOB * NOB + 1
MARK(NCR.1) « IWRITE
WplTf!6t»*93>N0B»lWRlTE
KOUNTT s KOUNTT ~ 1
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
002! 3
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
-------
34 CCNTfWUF	00239
£•#»•#••#*##•~#•*••*•*••»•#~*#«•#»•~»•*••~»•••••»•»«»»*»•*»»»•»*»«»»»#	00240
r» MAIN LOOP #	00241
OC 536 iCYC ¦ INCYC»NQCYC	00243
NQCYCC • ICYC	00244
C •««»•«••*•«•««••«**•«««•«•«••«••«•»*•	00245
c * Input hydraulics information »	00246
C	00247
R£A$<3! tO(N)« V (N) fN*1 »NlC)	0024ft
IF(ICYCTF.6E.NECT)REWIND 3	00249
READ{3) ICYCTF, jfYNfWtJ)	00250
C	00251
C f ADJUST FOR ADVEctlCN ANO DIFFUSION •	00252
C	00253
Op NslfNC	00254
V3LFi> « «(N) • DEL*0	00255
Nl ¦ NJUNC(N»1)	00256
NH ¦ NJUNC(Nf2)	00257
IF(N90T«D|OJi GC TO 406	00258
tF'(QCN)tQE*0«)80 TO 404	00259
FACtb»»Q«	00260
60 TO M2	00261
404 FACTOR S.1*00	00262
00 Tp 412	00263
406 IF(OfN),0E.O.)SO TO 410	00264
FACtOR « 0.25	00265
©0 TO 4*2	00266
*10 FACTQB ¦ 0.75	00267
41? 00 icyltNUHCON	00268
QGRAO « C - C * 06RAD	00272
CMASS(NHtK) « CMASS(NH,K) ~ ADMASS ~ DIMASS	00273
CMASS « CMASS(Nt.K) - ADMASS - DIMASS	00274
414 CONTINUE	00275
416 CONTINUE	00276
C •«•••»••••»••*«•••*•»»•»»	00277
r *	ApjusT FOR DECAy •	00278
C	00279
-------
IF INCQNDKCK) vLE,0)G0 TO 424	00281
NCDN. ¦ NCCNDKtKJ	00282
NCCNQ rn NCCNOX(K)	00283
00 4?0 J»NCJPlfNJ	00204
CMAS$(J,NCGN)»CMASS * OECAY	00285
iF(NCCNC»LC#0!SQ TO 4£0	00286
CMAS$(J»MCCNC) ¦ CMASSirJ•NCCNC) « C(j**CGN) # VOL  * COECAYffO	00287
• ~ RECXK••••••••••••••••••••«•«••*•••*••••«•••«••••••	00294
NfAG « NTAG ~ 1	00295
IF (Nf AG # (?E •NSPEC)NTAG«6	00296
OS 429 K«i*NUMCCN	00397
00 4|9 J«lfNOJ	00298
C(JfK> ¦ CIN(K»NTAG*l)	00299
429	CONTINUE		00300
C ••••«•*••*•«•«•••••»*••»*«•••••»«»**••*•«••••**••«•»••«•«	00301
C • ADJUST FOR WASTE OUTFLOWS• AND FOR WASTE INFLOWS #	00302
c * ,frqm external sources	* 00303
DC 434 j«NOjPltNJ	00305
IFt0Xfc|WQ(J))430f434#432	00306
430	OC 4^ K«t»NUHC0N	0030?
CMAS*(J«lO«CMASS(J»K) - CSPEC(J*K) • VCLQIN(j)	00308
431	CCNtlNUE	00309
GO TO *34	00310
432	DO 433 K«1»NUMCCN	00311
CMA$S«CMASS - CtJ,K> ~ VCLQlNtJ*	00312
433	CONTINUE	00313
434	CONTINUE		00314
C ••••••••••••••ft**********************	00315
c • adjust #or diversion returns •	oo3i6
C »••••••*••»•••••*••*•«•*«***««•«*••••	00317
IF00 TO 442	00316
DO 440 1*1tNUNlTS	00319
-------
JD2 » JDIV2(I)	003?1
¦ JRETKI)	00322
JR? a JRET2tI>	00323
DC 438 w|»ltNUMCCN _	00324
CMASS UR1#M)»CMASS(JR1,M>~(C tJDl,M)•VOLQTN  *	003*5
~	CCNSfUtM)	00326
CMASS{JR2tH)-CHASS(jR?,m>~IC(J02•M)•VCLQIN(JD2>«RET*NF(I,M))~	00327
» CCNSTlltH)	00328
43* CONTINUE	00329
440 CONTINUE	00330
442 IF(XTEMR.EQ.O)GC TO 443	00311
C	00312
C * ADJUST TEMPERATURE FOR NON-AHVECTIVF HEAT transfers #	00333
C •»••«••••#»••»•••••••••*•••»•*»*•»••••••«•«»•••*«•**•*«»•«»•	00314
CALL FMD	00315
DC 1500 J«NCJP1»NJ	00336
RHCWaiddO..~ B	00337
TWC ¦ C(JtlTEMP)	00338
HV«597*-*S7«TWC	00339
R0*DR«l40/t(VOL(J)*304.80061)/ASURUWjND«6.05	00345
T1«RHCW«HV»|A^BB»UWIN0)	00346
OC«J'»*fi*OELVAP	00347
IF(DELVAP«LT.0.6>0E(J)«d»6	00348
DELTMP«f fcfC-T A	00349
B0WHp0«p,61*Tl	00350
OC C J)¦80WM0D#0ELTHP	00351
QW(J)»7.36E-2*1.17E~3»TwC	00352
QOePBQEij}«QW«0.0)T#C»0.0	00359
-------
T?«BEfA(MNI~6»1E-4»aP	00361
T3«*iPH(MN)-EA-6•IE-4»AP»TA	00362
XCHCP«f.l7E-3*Tl*T2	00363
DNUM«0R-7*36E-2-Tl*T3	00364
eQTEMlJfONUM/XCHCF	00365 _
0TEM»f«C-E0fEMcj)	00366 £
TEM2«£0TEM(J)*0TEM«EXP(-C(XCHCF»0ELT0»*PCXDR)>	00367
ii-ifix(TEM?)/5*1 _	00368
IF(II.EQ.NN>GC TO 1500	00369
TWC-TEM?	00370
60 to 4	00371
1500 C9NTWE			00372
C •»••••••*•••*•«**•••«•«•~•*•»•••*«»*»»*•#*«»*«*•»•*«*#*»•#«•*«»««	00373
C # C?MPytf SPECIFIC CONCENTRATIONS FROM MASS CONCENTRATIONS •	00374
C #•«~«~•#•~•#•~•••~•»~•#««•«•««**~«•««~»«*»«»«»~**«*••*•»**«~*«*•«¦	00375
443 HO 446 ^«NCJPltNJ	00376
VCLiJ) « VGl(J> ~ ASUR(J)»(YNEW(J)-Y(J))	00377
00 444 K«1(NUMC0N	00378
C(JtK)"CMASS{J,K)/VOL(J)	00379
444 pONtlNllE	00380
446 CONTINUE	00381
C	••••••*•••«•*••*»••*•*«•*•*•«**•*••••«	00382
C	• CHECK NEGATIVE CONCENTRATIONS ~	00383
C	••••••••••«•*•»•»••«»*•«•••»«»•»«•»•••	00384
Op 466 j«l«NJ	00385
YUJ « iH€UWRITE<61*466>J»ICYC«K»C(J*K)	00389
CCJtK) * 0.0	00390
CMASSl J»iO" 0*0	00391
464 CONTINUE	00392
466 CONTINUE	00393
IP(NCCNDK tl)«EQ,Q)GO T$ *7*	00394
C	00395
C	^ CH5C* SUPERSATURATI9N •	00396
C	•••*•«*•«*•••••••«**•«•«*»••«•	00397
05 476 S«1#NUMCCN	00398
tr(NCONOK(K),EQ.O.CR.NCONCX(K).EO.O)6C TC 476	00399
-------
00 475 JmlfNJ	0040!
IF00 TO *75	00402
WR!te{Af»*74>NCCNtJ#lCYCtC(JtNCCN>	00*03
C(J»NCCN) » CSA7(K)	00404
CMAS§ ¦ C(J,NCCN) * VOL (J)	00405
475	CONTINUE	00*0*
476	CCNTjNUE	00407
479	CONTINUE	00408
DC 482 J»1»NJ	00409
C •••••••••••«•*«*•««««»**«•«•«««•«	00410
C • . CHECK OVER MAXIMUM LIMIT «	00411
q «#«•«•••*•»«*•#»•*•••••»*••»*»•»»	00412
DC *80 K»ltNUMCON	00*13
IF(C{J«K)•LE.CLIMITjK))60 To 480	00414
WRltE(6l»478)K»CLIMlT(K)tJtfCYC	00415
STOP	00*1&
480	CCNfjNUf	00417
48? CONTINUE	00418
IFt CICYC^NSPEC)-NQCYC)486«484•490	00419
C #*#**•••••##•••~##~#*•#»#~#»~#~~*~•»~~#~#~~•~•«•#«	004?0
C ~, ,mAKE.BINARY OUTPUT FOR EXTRACTING PROGRAM •	00421
C	00422
484 KOUNTT « 0	00*23
AA AAA	AA424
486 IFIICYCjUT.iwRlTE)©© to 500	00425
490 K0UNTT*15yUNtT~ 1	00426
IF(KOUNTT -J)492t492*494	00427
492 NOB • N08,*l	00428
MARK(N08»II • ICYC	00429
WITPnS#*93) NO® 11 CYC	00*30
494 IF«KOUNTT,tT.(NSPEC*l)>60 To 498	00411
MARK'(N§B92>?ICYC	00432
Wfimn6#*97) NoB*ICYC	00433
KCUNTT.q	00434
IWRITE « NEXTWR	00435
NEXTWR ¦ NEXTWR ~ IWRINT _	00436
498 WRITEnd) I CYC» ((C(J»K)»K«1 •NUMCON) *Jal«NJ)	00437
500 CONTINUE	00*38
IF(jtCYC»L.T«NQCYC)GC TO 520	00*39
C	00**0
c * make Restart deck #	00**1
-------
NFIRST-3	00443
IF(>HIHCSM.LT.3)NFXRST»NUMCCN	00444
WRffE«9ifd3i{ALPHA III«T«1«7?)	00445
00 556	00446
WRtTf«9#555> J#OIWWOCJ>»(C«J.K>,cSPEC	00447
556 CONTINUE	00448
IEiNUMCCN.LE«3>8C 10 517	00449
NFIRST « NFIRST ~ 1	00450
08 558 JaI*NJ	00451
WRITEI9»557) J»(C(J»K),CSPEC	0045?
558 CbwTINUe	00453
517 CONTINUE - . ~	00454
WRItf(l6»5l8)ICYCfICYCTF«NTAG	00455
END FILE 9	OO4**
REWIND 9	00457
520 CONTINUE . _	0045B
C ••*•••**•»••••«•«•••••••»•	00459
C ,CUfpUT TO PRINTER •	00*60
C *•••*«*•••**••*•••#•••*«»«	00461
IE(icYC*NSPEC*l.GE»NOCYC)GC TO 5?8	00462
IF(ICYCftT.lPRTjGO TO 536	00463
IPRT ¦ fpRT *NQpRT	00464
MC5UNT ¦ NCOUNT ~ 1	00465
IF(NCCUNf.LT.NSPEC/NQPRT*2)G0 TO 528	00466
NCOUNf » 0	00467
IPRT » NE*TPR	00468
. NfcxtpR • NEXTPR ~ INT8j$	00469
528 HOURS ¦ DELTO # FLOATF / 3600,0	00470
KP4YS ¦ HOURS / 23.99999	00471
HOURS ¦ HOURS - FL0ATF(f4 • KOAy*)	00472
WRitEil6»53d> ICVC,KOAYS»HOURS	00473
DC 534 |»1tNOPRf	00474
j«jPRT{I)	0O475
„ WRiTE<61t532)J»V(J)»CCrj»K)»K*l»NUMC5N)	00476
534 CONTINUE	00477
IF
-------
DO 535 j«l»NGPRT	00481
J»„JPRT{J)	004Q2
8TUl»0TptfJ)M327»29	00483
Rfu3«QW|j)•[327.29	00**4
BTU4aQE(J)*1327,29	00485
RTUSpQCiJ)*1327,29	00486
WRIT?(61«533)QTCT(J)•BTUl»QR«BTU?»QW(J)•BTU3,oE(Ji »BTU4f	00487
•	QC(J)i.8TU5»EQTEM(J)	00488
535 CONffNUE	00489
53ft CONTjNyi?	00490
00491
c* 	END MaIN.LOOP	•	00*92
00493
REWINP 3	00494
WRjf?(6i#546) ((MARK(J.KMK«l,2)»J«ltN0B)	00495
W®CYCC 		00496
C	00497
C	•	CHL.TP QUALITY EXTRACTION PRCQRAM GOES HERE #	00498
C	»•••*•••»••*••**••••••••»••»••••«*«••*•*••••«•«****«•	00499
STOP .	00500
40 FORMAT	00509
58 FORMAT(JHO/	005l0
•	lTHOCONSTITUfNT NO. |J«59H IS TREATED AS A N0N-CCN5ERV ATIVE	00511
•	WITH OCCAY COEFFICIENT « Fl0.7t45H BUT IS NOT PAIRED WITH Any CTH	00512
•ER C5^STITUENT//J	00513
80 FORMAT<715)	00514
B1 FORMftT(38H0N0 WASTE WATER RETURN FACTORS APPLIED//)	00515
84 FpRHAT(4|5»F10«6tl5)	005l6
103 FORMAT CIBA4)	00517
105 FORMATUHI////	00518
•	1H i*8A4#5Xt47H FEDERAL WATER POLLUTION CONTROL ArjMlNfSTRAT	00519
-------
•	1H 18A4»5X»37H PACIFIC NORTHWEST WATER LABORATORY /	00521
•	jfH 18A*/JN }8A4////)	00522
106 FORMAT I42H *•*••*•* FROM HYORAULfCS PROGRAM •••••»•*/	00573
•	4gH START CYCLE . STOP CYCLE TIME INTERVAL//	00524
•	I« I7fl!*fF!2.6t9H SECONOS/////)	00525
ldt FpRMATUI™ STARTING CYCLE IWltlAL QUALITY TOTAL OUALfTy *	00526
#•# OUTPUT INTERVALS *** TIME INTERVAL IN	CONSTANT FOR/	0052T
•	122H ON HYO, EXTRACT TAPE	CYCLE	CvClFS	00528
•	CYCLES	HOURS	QUALITY PROGRAM DIFFUSION COEFFICIENT	00529
•S//	.	. .	00530
•	113.118.116.113.F14.2.F17.3.6H HOURS.F17.3////)	00531
109	FORMAf C31H PRlNfOUT IS TO BEGIN AT CYCLF 14//	00532
•	49H quality tape for extracting is to begin at cyClet5////>	00533
110	FORMAT(5F10.0)	00534
112 FORMAf(j5)	00535
116 FORMAT(|3t3l4»5(F5,6.E«.2))	00536
120 FORMAT(1HOI5»42H CONSTITUENTS BEING CONSIDERED IN THIS RUN//)	00537
12? FORMAT<1H018A4)	00538
184 FORMAT(7F10.0)	00539
188 FORMAT(55H0SRECIFIED C-FACtORS aT jUncTiOn 1 FOR CONSTITUENT nO. I	00540
«!//	00541
•	UH TF12.2))	00542
192 FORMAT(fiiS)	00543
194	FORMATCl5#2F8.0»F9.6»F8.3.Fl2.2.Fl0.1tI9.l6)	00544
195	FCRM^f(S5Xtl9tF8.2tl8t4i5)	-	00545
196	FORMAT(|HI////	00546
•	..13?H **************************** CHANNEL DATA »*•»*•*•	00547
••*•••«•*««•*«•«•«•»«« ••*••«««•• JUNCTION DATA *	00548
••••••••••/	00549
•	132H CHAN. LENGTH WIDTH AREA MANNING NET FLOW HYD.	00550
•Radius Junc. at ends	junc. head channels entering	00551
•	JUNCTION//	_	00552
•	
-------
•NSTlfUENT/	00561
•	132HUNIT NO, 1 NO. 2 NO. 1 NO. 2 COEFF. CONST.	00562
•	CC£FFt CONST. COEFF. CONST, COEFF. CONST. COEFF.	00563
•	CONST?//)	00564
200 FpRMAf(XStTFlO.O)	00565
204 FORMAT(3fHODATA CARD OUT OF SEQUENCE, JJ« I4(3HVJ« U)	00566
216 FC*MATt52H0NC MULTIPLICATION FACTOR APPLIED TO CONSTITUENT Nc.12/)	00567
22o FORMAT	005*«
224 F0PHATC70Hl»»»»»MULTIPLiCATl0N FACTORS APPLIED To OBTAtN STARTING	00569
•CONCENTRATIONS//	00570
•	5fH CONSTITUENT ©ROUP FACTOR JUNCTION NUMBFRS)	00571
228 FORMAT OH //l8»Ill#FU.?til2.2H -,14/	00572
•	(119,FU.2*112.2H -,I4)>	00573
24i FORMAT(|HI////	00574
•	.UOH WATEP	00575
•QUALITY DATA	00576
•	12QH * FIrST cONStiTUENt • SECOND CONSTjTIJEnT	00577
•	* THIRD CONSTITUENT • FOURTH CONSTITUENT • FIFTH CONSTITUENT #/	00578
•	. U«H INITIAL INFLOW INITIAL INFLOW	00579
•	INITIAL INFLOW INITIAL INFLOW INITIAL INFLOW/	00580
•	U8H JUNC. INFLOW CONC. CONC. CONC. CCNC.	00581
•	. CvNC# CONC. CONC. CONC. CONC# CONC.//)	00582
282 FORMAT(I4«F|0*4fFl2»l«2Fid*l* Fli,lt3Fl0.1»Fn.l,2F10.1>	00583
350 F0RMAT(l3tI8*I79U0»I7vF9.2,El2«2t4(F7.2,E12$2))	00584
460 FORMATC39H DEPLETION CORRECTION MADE AT JUNCXION I3t7H CYCLE I4»	00585
•	21H FOR CONSTITUENT N0S U»12H# CONC. WAS FlO.2)	00586
474 FORMATC^HOSIlPERSAtllRATjON OF CONSTITUENT NO. I1.23H PoEvENTfD AT	00587
•JUNCTION 14,7H CYCLE I4.10H CONC, WAS F10.2//)	00588
478 F0PMAT(34H0C0NCENTRATlcN OF CONSTITUENT NO. Ilt8H EXCEE0S#F7.1.	00589
•	13H IN JUNCTION I3tl4H DURING CYCLE I5.25H. EXECUTION TERMINATE	00590
•0#I -	00591
493 FC*MAT	00592
497 FORMAt(£//6H HArK(I2»5H,2) «I5///)	00593
518 FORMAT(JH1////46HRESTART DECK TAPE WAS LAST WRITTEN AFTER CYCLEI*/	00594
•	50H HYDRAULIC CYCLE ON EXTRACT TAPE FOR PESTARTtNQ * 15/	00595
•	fH NTA6 * 13///)	00596
530 FORMAT
-------
*	CONCENTRATION FACTORS #*»*«*#»*»»¦»##»#•»•»•*####/	00601
~	J09H JUNCTION HEAD 1 ST. CCNSTIT. 2mD, CO^STI	00602
*T» ¦ _3?% CC'^TIT. 4TH. CONSTlT• 5TH. CONSTIT,/ 		00603
#	" 1 ISM NT IRE R (FT) (M~GL> (MGI.)	00604
•		00605
531 FORMAT(#1 #»»»«»«#»#»«»#» paDIATTON TERMS ANO CQU?L*»	00606 £
«MI«RTlW TFMPEHAtUHKS #»**»»»«»#«~»*»*/	00607
»*0fKCAL JMPI.IES KH.ogpam-CALORIES PEP SQUARE METER pER *,	00608
**SFCfINO. hTh IMPLIES otij p£r .SQUAREFCGT_P£P.HCURi */ „ 	_	QMQS
«*. NET 3«niATI0N 1MCGMIMG SCLAR BACK RADIATION FvAPCRA*,	OO610
•	*T!0!" CCNIM1CTICI FQUlL TEMP#/	00611
#1X«5(# KCAL HTM *)«* CFNTIGRADF*)	00612
533 FCRMAT.3.1X) ,fq.2)	00613
53? F0RMATn-ni5fFl;>.2fF?0.2»4Fl7,2)	00614
54Q F^MATn?'H00,'ALlTY TAPE. WAS WRTTtEN FROM CYCtF»T*»9H Tf! rYCLF« T6/)	00615
54? FORMAT(P^H-ENO OF QiiALTTY RUN.«T«;.9H CYCLES,)	00616
555 FCRMATU^ne.^bFlO.?)	00617
557 FCPMAT(I=..^F10.?)	00618
FNn	00619
susrcuttm^ metota	006.20
C »¦»«»»»»»»»»»»«»»<*»•»»»»»»»»»»»»«»»»»»»»»»»»«»»»»»»»»¦»»»»»»«»»»»»	Q06? 1
C • St 10RCUTINF TO TNpUT METFORCLOGICAL DATA ~ "• ~~	00622
C *	FC^ QUALITY PROGRAM				 »		00623
C •»•»~«»*»»»»***#*»##»**#»»»*»#•«•»»»#»~«»»»~»#«»»#»»~»•#»»»#»~»»	00624
HTME^SIOM ORNETa i?A) *tiwTMDA(24) ,tAA(?4> ,TAWA(?4) ,APA<£0	00625
COMMCM ALPHAU9«) »C (S,S) tCSPFC f5«5) »V0L<«5) rQiNWQTAWACf1 tAPA(I)		00633
lOO CCMTTNtJF	00634
jnosiFlx(OELr^)	00635
FINTaFLfATFiirlT)	00636
LCT=JNT»MPTs	00617
-------
,?^TKY FyO . . ... 		 			:			00639
£»#•#*«»*###»» »«•¦»*«#«¦««•##»««¦««#«»«»»#»»#»#»#«»#*»#»*#»#«~»*»*»**»»»~	00640
C« ENTRY POINT re FFTCM CURRENT (INTERPOLATED) VALUE	» _ 00641
C« ZF mETEORDLCPJTCAL VARIABLES •	00642
nnnnnn»"»»*»»*	00643
NOCSMsmoc^+1	00645
	JTlM=Knc5'**Im:?			 __		Q06$£
ITTTalfH-CITiM/LCTinisT	00647
T=TTIT/jMT _	_ 00648
FACTsFLCATF(ITTT-I«IMT)/FINT	00649
1 = 1 *1 			_ 00650
Isf+T	00651
TF < J.fiT^vPTS) J=1	___			00652
QRNETe(T A(J)-OWMET A(I))«FACT+QRNFTAt1)	00653
UWlMns	00659
11	F0RMATCF4.4,3F3.1,F4) 		00660
1? FCRMAT(*1 »»»»»***##»*#»•»»•»#»*»»«#•« TABLE CF Met^CrC#*	00661
• *LOGICAL DATA »*»»»»»#»#»»»»»«¦»«»¦»»»»»•»»~*/	00662
**0 MET DRY WET*/	00663
•* INCOMES irflNO RllLB BULR ATMOSPHERIC*/ 			_ 00664
•* RADIATION SPEED tEHP TEMP PRESSURE*/	00665
	•* 
-------
As the Nation's principal conservation agency, the Depart-
ment of the Interior has basic responsibilities for water,
fish, wildlife, mineral, land, park, and recreational re-
sources. Indian and Territorial affairs are other major
concerns of America's "Department of Natural Resources."
The Department works to assure the wisest choice in manag-
ing all our resources so each will make its full contribu-
-------