A DYNAMIC RIVER BASIN
WATER QUALITY
MODEL
EPA" 910/9-91-019"
John Yearsley
EPA Region 10
September 1991
-------
TABLE OF CONTENTS
Chapter Page Number
1 Model Development
1.1 Introduction 1
1.2 Mathematical Development 4
1.3 Method of Solution 19
1.4 References 22
2 Computer Program Description
2.1 Introduction 23
2.2 Data Preparation 23
2.3 Software Organization 24
2.4 References 34
Appendices
I Data Input Preparation
11 Fortran Source Code Listing
-------
LIST OF FIGURES
Figure Number Page Number
1.1 Schematic diagram of typical river basin
showing major hydrologic features
included in the model 3
1.2 Segmentation schemd for river or river-run
reservoir showing typical reach and typical
computational element 5
1.3 Segmentation scheme for stratified reservoir
showing typical reach and typical
computational element 6
1.4 Flow diagram for ecologic state variables
in the river basin model 8
2.1 Schematic diagram of branching river system 26
n
-------
LIST OF TABLES
Table Number Page Number
2.1 Options available for computing the
rate of reaeration 28
111
-------
CHAPTER 1
MODEL DEVELOPMENT
1.1 Introduction
Knowledge of the physical, chemical and biological processes,
acquired from laboratory and field studies, has greatly increased our
understanding of aquatic ecosystems. In the case of some state variables
which characterize these ecosystems, the knowledge has reached a level at
which it is possible to describe important features of the ecosystem in terms
of mathematical relationships. For those processes which can be described
in these terms, simulation of water quality with mathematical models can
be a useful tool for water resource planning. These relationships are
developed by combining the empirical results and concepts from the field
and laboratory with the fundamental laws of conservation of energy,
momentum and mass. Formulation of these relationships leads to the
equations for a mathematical model.
In their most general form, the equations developed to describe
aquatic ecosystems are extremely difficult to solve. However, by limiting the
range of the analysis to certain time and length scales and by using
numerical methods, it is possible to obtain approximate solutions to such
equations . The limitations in the model are invoked by making certain
assumptions. While invoking assumptions may make the solution of the
equations more tractable, it will also limit the range of applications of the
model. It is, therefore, important for those using the model to understand
both the capabilities and the limitations of an ecosystem model.
The model described in this report makes use of ecosystem concepts
which have been used in other modelling efforts (e.g., Thomann et al., 1975;
-------
Patten et al., 1975., DiToro et al., 1975; Chen and Orlob, 1975; Scavia, 1980).
Based upon these concepts, the mathematical model described in this
report has been developed to simulate the following state variables:
(1) Carbonaceous biological oxygendemand
(2) Dissolved oxygen
(3) Algal biomass
(4) Organic nitrogen
(5) Ammonia nitrogen
(6) Nitrite + nitrate nitrogen
(7) Organic phosphorus
(8) Orthophosphorus
(9) Temperature
(10) Coliform bacteria
(11) Conservative constituent #1
(12) Conservative constituent #2
The hydrologic setting is that of a river basin within which there can
be free-flowing river segments, river-run reservoirs and stratified
reservoirs (Figure 1.1). The model simulates the time history of the ecologic
state variables dynamically for time scales of hours and greater. Length
scales of computational segments can be of the order of 100's of feet
longitudinally and 5-10 feet vertically. The ability of the model to resolve
changes at these time and length scales depends, of course, on the
availability of appropriate data, as well as on the structural accuracy of the
model.
-------
RESERVOIR
POINT SOURCE
NON-POINT SOURCE/
GROUNDWATER/
WITHDRAWAL
Figure 1.1. Schematic diagram of typical river basin shown major
hydrologic features included in model.
-------
1.2 Mathematical Development
The proper application of a mathematical model requires a
knowledge of the model's capabilities and limitations. These limitations
are determined by the assumptions upon which the model has been based.
The general assumptions associated with this model are:
Horizontal and vertical advection and vertical eddy diffusion are
the primary physical processes for water and mass transport
The vertical eddy diffusivity is the same for all state variables
The lateral variations of properties in the waterbodies are
negligible compared to longitudinal and vertical variations of the
properties
Rate constants for the various reactions do not change over a given
length segment
Hydrodynamic characteristics are a function of the stream , river,
or reservoir geometry, only
The river system can be divided into a finite number of segments
within which hydrodynamic characteristics are constant
Hydrodynamic characteristics of free-flowing river segments and
river-run reservoirs can be expressed as a simple function of the flow
in any segment
Hydrodynamic characteristics of stratified reservoir segments are
a function of the density structure of the reservoir
The time required for flow in a reach to adjust to changes in
elevation is small compared to the travel time of some constituent.
Another way of viewing this is in terms of the speed of the gravity
wave carrying elevation information compared to the average river
velocity.
Simulated state variables of the ecosystem are averages over a
given computational element (Figure 1.2 for free-flowing rivers or
river-run reservoirs and Figure 1.3 for stratified reservoirs) and a
finite time interval
-------
-------
Figure 1.3. Segmentation scheme for stratified reservoir
showing typical reach and typical computational element.
-------
For some state variable, C, which is time- and length-averaged over a
computational element, the general conservation equation in the ijt
flowing river or river-run reservoir segment can be written as:
d(CV) A
_ = A(QC)X + (Q C)n + O.- F.
Similarly, for the ij^ stratified reservoir segment, the general
conservation equation is
d(CV).. NS
^-1= A(QC)X + A(QC)Z + A(KAC)Z + 2, (QpCp)n + % - rij d-lb)
n=l
where,
V = the volume of the ijth computational element where i refers
to the segment number and j to the computational element
number within the segment,
C = the length- and time-averaged value of some state variable
over the ij**1 computational element,
A(QC)X = the advective transfer in the longitudinal (x-) direction,
A(QC)Z = the advective transfer in the vertical (z-) direction,
A(QC)n = the transfer of flows from the computational element
due to inputs or outputs such as point source discharges,
non-point source return flows and withdrawals for
drinking water or irrigation,
A(KAC)Z= the eddy diffusion in the vertical (z-) direction,
-------
A = the surface area of the ijth element,
Ojj = the source term for the state variable, Cij,
FJJ = the sink term for the state variable, GJJ.
For all state variables, gains and losses due to the physical processes
of advection and eddy diffusion are treated in the same manner. The
source, ij, and sink, Fij, for each of the state variables are determined
from existing knowledge of physical, chemical and biological processes. A
flow diagram of the interaction among state variables is shown in Figure
1.4. A complete description of the source and sink terms in the mass
balance of each state variable is given below.
Carbonaceous Biological Oxygen Demand (CBQD). Ci
The major sink term for (CBOD) included in this model is the
following:
Stabilization of CBOD by microorganisms
The stabilization of CBOD is represented by a first-order,
temperature-dependent process. The differential equation for the mass
balance, excluding the physical processes of advection and diffusion, is as
follows:
d(CV)
--=-Kcv (L2)
where,
8
-------
SOLAR ENERGY 1 | NUTRIENTS
r
...] L
ORGANIC
N
PHYTO-
PLANKTON
ATMOSPHERIC
(^ DO
NBOD
DO
^
CBOD
ORGANIC
P
Figure 1.4. Flow diagram for ecologic state variables in the river basin model.
-------
20
K, =K, 9.T
20 o i
K = the deoxygenation rate at 20 C, days"
01T = the temperature correction factor for deoxygenation.
Dissolved Oxygen (DO). 9
The major source and sink terms for DO include the following:
Stabilization of CBOD by microorganisms
Nitrogenous biological oxygen demand (NBOD)
Respiration of phytoplankton
Photosynthesis by phytoplankton
The mass balance is
d(C V)
where,
K.2 = the reaeration rate, days'1,
20
= K2 °2T
20 o
K = the reaeration rate at 20 C
(C9-20.0)
62T = 1.024
Csat = saturation level of dissolved oxygen, mg/1
10
-------
= the stoichiometric relationships between oxygen and
nitrogen for the oxidation of ammonia to nitrate,
= the nitrification rate for converting ammonia to nitrate,
days"1.
Algal Biomass. 03
The major elements characterizing the dynamics of algal biomass
are
growth driven by energy from sunlight in the presence of the
macronutrients, nitrogen and phosphorus
respiration of organic carbon stores
settling of plankton due to gravitational influences
Mathematical formulation of the mass balance for these processes is
given by
d(C V) w
_L_=(G.R._±)C3V (1.4)
where,
G = Gmax fT(C9) fid) fN(C5,C6) fp(C6)
= the maximum growth rate for alglal biomass, days'1,
= the function describing the temperature-dependency of the
algal growth rate
T r
-e
11
-------
T -
opt
.3( _
whenC9>Topt
fl(I) = the function describing the dependency of the growth rate on
solar radiation
*o ~^Z2 0 ~YZi
fp = the photo period, fraction of days,
y = the extinction coefficient, meters'1,
IQ = net solar radiation at the water surface, kcal/meter2/second,
Is = the optimal radiation for algal growth, kcal/meters2/second,
zi = depth below surface of top of the ijth element, meters,
Z2 = depth below surface of botttom of the ijth element, meters.
12
-------
growth limiting factor for nitrogen,
C. + CL
(K, + C. + CJ
A o 6
= the half-saturation constant for nitrogen, mg/1 N,
C$ = the concentration of ammonia nitrogen, mg/1 N,
Cg = the concentration of nitrate nitrogen, mg/1 N,
fp = growth limiting factor for phosphorus,
C
Kp = the half-saturation constant for phosphorus, mg/1 P,
Cg = the concentration of orthophosphate, mg/1 P.
Organic Nitrogen. 04
The major sources and sinks for organic nitrogen are
waste products due to respiration of algae
mineralization to ammonia nitrogen due to bacterial action.
The corresponding equation for mass balance is
d(C V)
13
-------
where
K44 = the mineralization rate of organic nitrogen, days'1,
= K20e
44 4T
20 o
K = the mineralization rate at 20 C
44
(C -20.0)
04T = 1.084
ctNC = the nitrogen/carbon ratio in algae,
Ammonia Nitrogen.
The major source and sink terms for ammonia nitrogen are
mineralization of organic nitrogen to ammonia
nitrification of ammonia to nitrate
uptake of ammonia by algal growth.
The mass balance equation is written
d(C V)
GC3 + K44 C4 - K55 ^ V (1.6)
3
f = algal preference factor for ammonia uptake
3
= the nitrification rate for converting ammonia to nitrate,
days'1
14
-------
(Cg - 20.0)
65T = 1.084
Nitrate Nitrogen. Cg
The major sources and sinks for nitrate nitrogen are
uptake of ammonia by algal growth
nitrification of ammonia to nitrate.
The mass balance equation is
d(C V)
where,
K56
Organic Phosphorus. C?
The major sources and sinks for organic nitrogen are
Waste products due to respiration of algae
Mineralization to organic phosphorus due to bacterial action.
The mass balance equation is
d(C7V)
_L_=(aRC-KC)V (1.8)
where,
K56
15
-------
ape = the phosphorus/carbon ratio in algae,
= the mineralization rate of organic phosphorus, days'1.
Qrthophosphate. Cg
The major sources and sinks for orthophosphate are
mineralization of organic phosphorus
uptake of orthophosphate by algal growth.
The mass balance equation is
d(C V)
(1.9)
Temperature. 9
The heat budget method is used to simulate changes in water
temperature in the river basin. The elements of the heat budget include
Net short wave radiation, qsn,
Net atmospheric (long wave) radiation, qat,
Water surface (long wave) radiation, qw,
Evaporative heat flux, qe,
Convective heat flux, qc.
The mass balance equation is
d(C V)
16
-------
where
p = the water density, kg/meters3,
Cp = the specific heat capacity of water, kcal/kg/°C,
= Qsn + Qat -Qw - Qe
= the net heat flux across the air- water interface,
kcal/meters2/second,
Az = the surface area of the computational element, meters2.
The methodology used to estimate the individual components of the
heat budget is similar to that described by Water Resources Engineers, Inc.
(1967).
Coliform Bacteria. Cm
The major source and sink terms for coliform bacteria are
mortality due to hostile environmental conditions
The mass balance equation is written
d'ciov'
dt B 10
where,
KB = the mortality rate for bacteria, seconds'1,
= 2.31xlO-6(l + .01111 Cio)
17
-------
Conservative Constituents. Cn and CT?
For the conservative constituents, the only processes affecting changes
in concentration are those of dilution, advection and diffusion. These
processes are treated in the same way for all constituents.
18
-------
1.3 Method of Solution
The state-space formulation for the ecologic model described above is
based upon the conservation equation for a well-mixed control volume of
finite dimensions. The differential equation for each state variaable can be
written in the following form:
After rearranging, eq. (1.12) can be written as
(1.13)
The conservation equations for all twelve state variables lead to the
following system of first-order, nonlinear differential equations:
dC
(1.14)
19
-------
The two-step Runge-Kutta method (Press et al, 1986) is used to solve
this system of equations for each computational element, beginning at the
first headwater reach and working downstream through the entire system in
the sequence specified in the problem description. The accuracy and the
stability of the solution will depend upon the time and space increments used
to characterize the problem. Because this is basically an explicit formulation
in time, stability criteria are associated with NI, the ratio of the residence
time of a computational element and the computational interval, where
N =^
iN
At = the computational interval,
Q = the net volume associated with the computational element,
V = the volume of the element.
and N2, the ratio of the diffusion time between elements and the
computational interval
N2= At
Az = the thickness of the computational element,
K = the vertical coefficient of eddy diffusivity.
Stability criteria associated with NI and N2 apply to the stratified
reservoirs. Advection and dilution are the only hydrodynamic processes
affecting concentration in the free-flowing and river -run reservoir
20
-------
segments. Therefore, only the criterion associated with NI applies to these
segments. Since the system of equations is nonlinear, exact criteria cannot
be derived. Results from the analysis of the linearized difference equations,
as well as common practice, suggests NI should be approximately one (1.0)
for best results and that instabilities will occur for values larger than one. N2
should be less than 0.50 (Bella and Dobbins, 1967) to prevent instabilities from
occurring in the solution.
21
-------
1.4 References
Bella, D.A. and W.E. Dobbins. 1967. Finite difference modelling of river and
estuary pollution. In: Proceedings of national symposium on
estuarine pollution - Stanford University. Stanford University,
Stanford, CA., 612-645.
Chen, C.W. and G.T. Orlob. 1975. Ecologic simulation for aquatic
environments. In: Systems Analysis and Simulation in Ecology, Vol
III. Academic Press, New York, NY. 354 pp.
DiToro, D.M., D.J. O'Connor and R.V. Thomann. 1975. Phytoplankton-
zooplankton-nutrient interaction model for Western Lake Erie. In:
Systems Analysis and Simulation in Ecology, Vol. III. Academic
Press, New York, NY. 354 pp.
Patten, B.C., D.A. Egloff and T.H. Richardson. 1975. Total ecosystem model
for a cove in Lake Texacoma. In: B.C. Patten (Editor), Systems
Analysis and Simulation in Ecology, Vol. IV Academic Press, New
York, NY., 329-371.
Press, W.H., B.P. Flannery, S.A. Teukolsky and W.T. Vetterling. 1986.
Numerical recipes, the art of scientific computing. Cambridge
University Press. 818 pp.
Scavia, D. 1980. An ecological model of Lake Ontario. Ecologic Modelling,
8, 49-78.
Water Resources Engineers. 1967. Prediction of thermal energy
distribution in streams and reservoirs. Prepared for the State of
California Department of Fish and Game. 88 pp.
22
-------
CHAPTER 2
COMPUTER PROGRAM DESCRIPTION
2.1 Introduction
The river basin model described in Chapter 1 has been designed to
analyze the impact of
point source wastes from industries and municipalities
non-point sources
water diversions
upon the aquatic ecosystems of freely-flowing rivers, river-run reservoirs
and stratified reservoirs. The model design is based upon the river basin
concept and the software provides the capability of analyzing branching
systems. This volume describes way in which the software implementing
the model is structured.
2.2 Data Preparation
For the purposes of preparing input data, the river basin being
examined must be first divided into reaches or segments. Within a given
reach, the hydraulic characteristics and base reaction rates are constant.
The quality and quantity of non-point source, or distributed inputs are also
assumed to be constant throughout each reach. In turn, each reach is
divided into a number of computational elements. The size of the
computational element is specified by the user according to requirements
for achieving a stable solution, as well as for resolving changes at a scale
consistent with the needs of river basin planners.
Water quality and quantity for all headwaters of the river system,
point sources, non-point sources and diversions must be determined. For
-------
those reaches which are characterized as stratified reservoirs, the
operating schedule of the reservoir must be available. Rate constants for
chemical and biological reactions in each reach must be determined. In
freely-flowing river segments or river-run resei'voirs, the relationship
between river flow and river depth and river flow and river velocity must be
specified. For stratified reservoirs, it is necessary to know the relationship
between volume and stage and surface area and stage. Sediment oxygen
demand in each reach must also be estimated.
The heat budget method is used to simulate water temperature. The
necessary components of the heat budget can be determined using the
methods described by Water Resources Engineers, Inc. (1967). The
software is designed so there can be a number of meteorological provinces.
This makes it possible to simulate temperatures in those river basins for
which the geographical extent is such that meteorological conditions may
vary substantially from one region to another.
2.3 Software Organization
Subroutine BEGIN
Data describing the basic structure of the river basin are read in this
subroutine. This includes system parameters such as an alphanumeric
description of the problem, computation time interval, simulation period
parameters and number of meteorological observations per day. Next,
parameters describing the topology of the river basin network are read. The
sequence in which these data are entered is important because the program
logic computes water quality in the same sequence. The first reach should
be the headwaters reach of the main stem, followed in sequence by
downstream reaches until a confluence is encountered. The reach
24
-------
containing the confluence must have a unique, though not necessarily
sequential number, NJUNC(N). The reach following must be a headwater
reach of one of the tributaries framing the confluence. If there is more than
one tributary in the confluence, the order in which tributary headwaters
are considered is not important. From the tributary headwater the
sequence is downstream until another confluence is encountered. When
all the headwaters which form the confluence have been included, the
sequence continues downstream on the main stem.
As an example, consider the network of reaches A, B, C, D, E, F, G,
H, I, J, K, L, and M, shown in Figure 2.1. Correct sequences include the
following:
(1) A, B, C, D, E, F, G, H, I, J, K, L, M
(2) B, A, D, C, E, F, G, H, K, I, J, L, M
(3) C, D, E, A, B, F, G, H, J, K, I, L, M
(4) D, C, E, B, A, F, G, H, I, K, J, L, M,
as well as a number of others.
Within each reach the user must specify a unique number for the
headwater, NHEAD(N), if the reach is a headwater; the number of point
sources, NPOINT(N), within the reach; the number of diversions in the
reach, NDIV; a unique, but not necessarily sequential number, NRESRV, if
the reach is a reservoir reach; and the number of computational elements,
NCELM, in the reach.
The program then reads a number of important parameters
characterizing the dynamics of the ecosystem. These parameters include
25
-------
H
Figure 2.1. Schematic diagram of branching river system.
26
-------
the oxygen/carbon ratio in photosynthesis, OCRAT; the carbon/chlorophyll
a ratio, CCLRAT; the nitrogen/carbon ratio in algae, NCRAT; the
phosphorus/carbon ratio in algae, PCRAT; half-saturation constants for
nitrogen and phosphorus, KMN and KMP; the fraction of ammonia in
nitrogen uptake by algae, NH3PRF; optimal solar radiation for algal
growth, QOPT; sinking speed of algae, WSINK; optimal, minimum and
maximum water temperatures for algal growth, TOPT, TLO, TUP; and
maximum algal growth rate and respiration rate, PGO and PRESO.
The rates for chemical and biological parameters include: the
reaeration rate, XKDO; the deoxygenation, rate, XKBODL; the rate of decay
of organic nitrogen, XKN44; the rate of decay of ammonia nitrogen, XKN55;
the deposition rate of phosphorous, XKP77; the sediment oxygen demand,
SOD; and the background light extinction coefficient, EXCO.
Several options are available for specifying the reaeration rate,
XKDO. The user may specify the value by inputting the desired positive
number in units of days'1 (base e). Input of a negative number for XKDO
will result in the computation of the reaeration rate according to one of the
formulae given in Table 2.1. The values of all rates are adjusted for
temperature as described in Volume I.
27
-------
Table 2.1 Options available for compute the reaeration rate, K2
XKDO* Reaeration Rate, K2
(days'1, base e)
Reference
>0.0
= XKDO
-1.0
-2.0
-3.0
K _ 2.3*5.026 U
2 1.673
.0.969
0.5
D
"
2.3*9.4 U
0.67
D
1.85
Churchill et al (1962)
O'Connor and Dobbins (1958)
0wens et al(l964)
See Card Group II, Type 6 in Appendix I
28
-------
The geometric characteristics of the reach are also input in this
routine. If the segment is a free-flowing river reach or a river-run
reservoir, coefficients relating velocity and depth to flow are required. The
coefficients, DEPTH1 and DEPTH2, relate the depth of the reach, D, to the
flow, Q, in the following way:
D=DEPTH1*QDEPTH2
The coefficients, VEL1 and VEL2 relate the velocity of the reach U, to
the flow, Q as
U = VEL1 * QVEL2
If the segment is a stratified reservoir, the lowest level of active
storage, Z(NV,1), the thickness of each computational element, ZLAYER,
and volume coefficients, AVC(NV) and BVC(NV) are input. The volume
coefficients are used to compute reservoir volume, V, at any depth N
according to
V = AVC(NV) + (Z(NV,N)-Z(NV,1))*BVC(NV)
The quantity and quality of non-point sources are specified by the
parameters, QRET and CRET(N). Headwater quantity and quality are
specified by QHEAD and CHEAD(N). Point source quantity and quality are
QPOINT and CPOINT(N), respectively, and the location of the point sources
29
-------
in i*iver miles is RMP. Quantity of diversion water is QDIV and the location
of the diversion in river miles is RMDIV.
The last card for the segment is text 'END1. The software keys on this
delimiter to indicate to the user whether the amount of information
provided is consistent with the amount required. If it is not, the software
issues a diagnostic indicating the reach in which the inconsistency occurs.
Subroutine SYSTEM
Subroutine SYSTEM maintains control of the simulation and output
portions of the program. Checking to determine whether the reach to be
simulated is a river segment or a stratified reservoir segment is performed
as are summations of water and water quality for junctions. Calls to
RESMOD and RIVMOD are made from this routine depending upon
whether the reach has been described as a river-run reservoir or a
stratified reservoir.
Subroutine QUALTY
Numerical solutions to the first-order differential equations for
ecosystem state variables are obtained in this subroutine using the two-step
Runge-Kutta method (Press et al, 1986). The subroutine software
implements the mathematical development described in Chapter 1.
Subroutine QUALTY is called from either RIVMOD or RESMOD to advance
the state variables of a computational element ahead one time step. The
results obtained in QUALTY are stored for use as initial conditions for the
next time step and for output, if desired.
30
-------
Subroutine WRITE 1
Subroutine WRITE 1 provides the output function for the ecosystem
model, software. Entry at WRITE 1 occurs at the beginning of each reach,
after all reach parameters, waste loads and diversions have been included.
Entry at WRITE2 occurs at the end of each computational element to print
the predicted values of the water quality constituents and the status of the
water budget. Entry at WRITES occurs at the end of the reach to print
surface elevation, water depth and velocity, dissolved oxygen saturation and
algal dissolved oxygen production and respiration.
Subroutine RESMOD
The physical processes of dilution, advection and turbulent diffusion
are developed for each computational element in the stratified reservoir
segment. Inflows are assigned to computational elements (layers) within
the reservoir segment based upon their density. This determination is
made by a call to Subroutine LAYRIN with the temperature of the inflow as
an argument. A call is also made to Subroutine MIX to see if the reservoir
density profile is stable. This is of importance during the fall as the surface
layers begin to cool, initiating overturn in the reservoir.
Subroutine RESMOD accounts for surface phenomena associated
with transfer of thermal energy and solar radiation. After performing
these calculations, a call is made to Subroutine QUALTY to advance the
estimate of the state variables one computational time increment. Upon
returning from Subroutine QUALTY, the current time level is compared
with time level for which output has been defined. If the two time levels
31
-------
match, The appropriate entry point in Subroutine WRITE 1 is called with
output to the printer and/or the plotter file.
Subroutine RIVMOD
Subroutine RIVMOD performs functions for the free-flowing river
and river-run reservoir segments similar to those performed in Subroutine
RESMOD for stratified reservoirs. There is, however, no vertical diffusion,
nor is there any account kept of vertical density structure. This is because
each segment in the free-flowing river and river-run reservoir segments
are assumed to be well-mixed both laterally and vertical.
Subroutine LAYRIN
The density, p, of inflowing water to a stratified reservoir is estimated
from the relationship
P =
(C9+283.)(C9-3.98)2
(503.57(C9+67.26)
The density of inflow water is compared with density of each
computational segment, beginning at the surface and proceeding
downward. The segment into which the inflow is placed is the first one
encountered for which the segment density exceeds the density of the
inflow.
Subroutine MIX
After the state variables have been projected ahead one time step in a
stratified reservoir, a call to this subroutine is made from Subroutine
32
-------
RESMOD to determine if there are instabilities in the density profiles. An
instability is defined as condition for which the density decreases with
depth. If such an instability is found, the reservoir is mixed uniformly
from the surface to the level of the instability. The resulting concentrations
are returned to Subroutine RESMOD as the updated state variables for the
next time step.
Subroutine ENERGY
Meteorologic data including net solar radiation, QNS, net
atmospheric radiation, QNA, dry bulb temperature, DBT, wind speed,
WIND, and vapor pressure EA are used to calculate thermal exchange
between the computational element at the water surface and the
atmosphere. The resulting heat budget is a source term in simulating the
water temperature.
2.4 Description of Input Data
Instructions for preparing the input data are given in Appendix I.
2.5 Source Code
The software is written in FORTRAN 77. An effort has been made to
simplify the code so as to be easily transported to other FORTRAN
compilers. A listing of the source code is given in Appendix II.
33
-------
2.6 References
Churchill, M.A., H.L. Elmore, and R.A. Buckingham. 1962. The
prediction of stream reaeration rates. ASCE Journ of Sanitary
Engineering, SA-4, 88, 1-46
OjConnor, D.J. and W.E. Dobbins. 1958. Mechanism of reaeration in
natural streams. ASCE Trans., 123. 641-684.
Owens, M., R.W. Edwards and J.W. Gibbs. 1964. Some reaeration studies
in streams. Int. Jour. Air and Water Pollution. 8_, 469-486.
Press, W.H., B.P. Flannery, S.A. Teukolsky and W.T. Vetterling. 1986.
Numerical i-ecipes, the art of scientific computing. Cambridge
University Press. 818pp.
Water Resources Engineers. 1967. Prediction of thermal energy
distribution in streams and reservoirs. Prepared for the State of
California Department of Fish and Game. 88 pp.
34
-------
Appendix I
Data input formats for the
river basin ecosystem model
-------
Card Group I
Card Variable
Type Columns Format Name
Description
1 1-80
2 1-5
6-10
11-15
16-20
20A4
F5.0
F5.0
F5.0
F5.0
XTITLE
DT
DAY1
DAY2
WOBSPD
Alphanumeric description of
analysis being performed.
Simulation time interval, days
First day of simulation
Last day in simulation
Meteorological observations
4
(If
ZPLOT>0)
21-25 F5.0 DAYPRT
1-5
16-20
1-5
6-10
F5.0 REACHX
6-10 F5.0 LAT
11-15 F5.0 ZLOW
F5.0 ZPLOT
15
15
NPLOT(l)
NPLOTC2)
per day
Time interval for printed
output, days
Number of reaches to be
simulated
Average latitude of river basin
Minimum thickness in the
surface layer of a stratified
reservoir
Switch indicating number of
reaches for which a plot is
desired.
Number of first reach for
which plot is desired
Number of second reach for
which plot is desired
15 NPLOTC.) Number of ZPLOTth reach for
which plot is desired
1-1
-------
Card Group I (continued)
Card Variable
Type Columns Format Name
Description
1-5
F5.0 OCRAT
6-10 F5.0 CCLRAT
11-15 F5.0 NCRAT
16-20 PCRAT
21-25 F5.0 KMN
26-30 F5.0 KMP
31-35 F5.0 NH3PRF
36-40 F5.0 QOPT
41-45 F5.0 WSINK
46-50 F5.0 TOPT
51-55 F5.0 TLO
56-60 F5.0 THI
61-65 F5.0 PGO
66-70 F5.0 PRESO
Oxygen:carbon ratio in
photosynthesis
Carbon:chlorophyll a ratio for
algae
Nitrogen:carbon ratio for algae
Phosphorus:carbon ratio for
algae
Half-saturation constant for
nitrogen, mg/1 N
Half-saturation constant for
phosphorus, mg/1 P
Fraction of NH3 in inorganic
nitrogen utilized by algae
Optimal radiation level for
algae, kcal/m2/second
Algal sinking rate, meters/day
Optimal temperature for algal
growth, °C
Minimum temperature for
algal growth, °C
Maximum temperature for
algal growth, °C
Maximum growth rate for
algae, days'1
Maximum respiration rate for
algae, days'1
1-2
-------
Card Group II (continued)
Card Variable
Type Columns Format Name
Description
1-20 A20 RNAME(N)
21-25 F5.0 RMILEl(N)
26-30 F5.0 RMILE2CN)
31-35 F5.0 ELEV(N)
1-5
15 NHEAD
6-10 15 NPOINT(N)
11-15 15 NDIV(N)
16-20 15 NJUNC
21-25 15 NRESRV(N)
26-30 15 NCELM(N)
Alphanumeric description of
reach
Beginning (upstream) river
mile of reach
Ending (downstream) river
mile of reach
Elevation of surface of reach
above Mean Sea Level (MSL),
feet
Identification number for
headwaters reach. Must be
unique, but does not have to be
sequential.
Number of point sources in
reach.
Number of diversions in reach
Identification number if
downstream boundary has a
confluence with other reaches.
Does not have to be sequential,
but all reaches with a common
confluence must have the
same identification number
Identification number if reach
is a stratified reservoir. Must
be unique, but does not have to
be sequential.
Number of computational
elements in reach. If
NRESRV(N)>0, the maximum
number of computational
elements in the stratified
reservoir
1-3
-------
Card Group II (continued)
Card Variable
Type Columns Format Name
Description
5 31-35 15 NWPROV(N)
6 1-10 F10.0 XKDO(N)
11-20 F10.0 XKBODL(N)
21-30 F10.0 XKN44CN)
31-40 F10.0 XKN55(N)
41-50 F10.0 XKP77(N)
51-60 F10.0 XKBACT(N)
61-70 F10.0 SOD(N)
71-80 F10.0 EXCO(N)
****************************
* Use Card Type 7a if NRESRV=0
7a 1-10 F10.0 DEPTHl(N)
Identification number for
meteorological province in
which reach is located.
Reaeration rate, days*1
Deoxygenation rate, days'1
Rate of mineralization of
organic nitrogen, days'1
-1
11-20 F10.0 DEPTH2CN)
Nitrification rate, days
Rate of mineralization of
organic phosphorus, days'1
Rate of dieoff for coliform
bacteria, days'1
Sediment oxygen demand,
mg/l(O2)/day
Background light extinction
coefficient, meters'1
*********************
*
*********************
Depth coefficient, Dl, in the
formula:
Depth=Dl*FlowD2
Depth coefficient, D2, in the
formula given above. Depth is
in feet, Flow is in cubic
feet/second
1-4
-------
Card Group II (continued)
Card Variable
Tvpe Columns Format Name
Description
7a 21-30 F10.0 VELl(N)
(continued)
31-40 F10.0 VEL2(N)
41-50 F10.0 DEPTHO
41-50 F10.0 WIDTHO
Velocity coefficient, VI, in the
formula:
Velocity=Vl*FlowV2
Velocity coefficient, V2, in the
formula given above. Velocity
is in feet/ second
Initial depth of reach, feet
Initial width of reach, feet
*************************************************
* Use Card Type 7b if NRESRV*0 *
*************************************************
7b
MO F10.0 Z(NV,1)
11-20 F10.0 ZLAYER
21-30 F10.0 AVC(NV)
31-40 F10.0 BVC(N)
41-50 F10.0 ZOUT
51-60 F10.0 QQRES
Elevation, feet above MSL, of
lowest, active portion of
reservoir
Thickness of computational
elements in stratified
reservoir, feet
Volume coefficient, Al, in the
formula:
Volume=Al+Bl*Depth
Volume coefficient, Bl, in the
above formula. Volume is in
cubic feet and Depth is feet
above Z(NV,1)
Elevation of reservoir
discharge, feet above MSL
Reservoir discharge, cfs
1-5
-------
Card Group II (continued)
Card Variable
Type Columns Format Name
Description
8
9
9
1-10 F10.0 QRET(N)
Quantity of ground return flow
to the reach, cubic feet/second
1-5 F5.0 CRET(1,N) CBOD of return flow, mg/1
6-10 F5.0 CRET(2,N)
11-15 F5.0 CRET(3,N)
16-20 F5.0 CRET(4,N)
21-25 F5.0 CRET(5,N)
26-30 F5.0 CRET(6,N)
31-35 F5.0 CRET(7,N)
36-40 F5.0 CRET(8,N)
41-45 F5.0 CRETC9.N)
46-50 F5.0 CRET(10,N)
51-55 F5.0 CRET(11,N)
56-60 F5.0 CRET(12,N)
DO of return flow, mg/1
Algal biomass of return flow,
mg/1 C
Organic nitrogen in return
flow, mg/1
Ammonia-nitrogen in return
flow, mg/1
Nitrate-nitrogen in return
flow, mg/1
Organic phosphorus in return
flow, mg/1
Orthophosphate-phosphorus
in return flow, mg/1
Temperature of return flow, °C
Coliform bacteria in return
flow, MPN
Conservative constituent #1 in
return flow
Conservative constituent #2 in
return flow
1-6
-------
******************* **
*
*
Card Types 10 and 11 should be omitted if the reach is not a
headwater reach (NHEAD=0)
*
*
*
*
Card Group II (continued)
Card Variable
Type Columns Format Name
Description
10 1-10
11 1-5
6-10
11 11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
F10.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
QHEAD(N)
CHEADU.N)
CHEAD(2,N)
CHEAD(3,N)
CHEAD(4,N)
CHEAD(5,N)
CHEAD(6,N)
CHEAD(7,N)
CHEAD(8,N)
CHEAD(9,N)
CHEAD(10,N)
Headwater flow, cfs
CBOD of headwater flow, mg/1
DO of headwater flow, mg/1
Algal biomass of headwater
flow, mg/1 C
Organic nitrogen in headwater
flow, mg/1
Ammonia-nitrogen in
headwater flow, mg/1
Nitrate-nitrogen in headwater
flow, mg/1
Organic phosphorus in
headwaters flow, mg/1
Orthophosphate-phosphorus
in headwaters flow, mg/1
Temperature of headwaters
flow, °C
Coliform bacteria in
51-55 F5.0 CHEADdl.N)
headwaters flow, # of
coliforms/100 ml
Conservative constituent #1 in
headwaters flow
1-7
-------
Card Group II (continued)
Card Variable
Type Columns Format Name
Description
11 56-60 F5.0 CHEAD(12,N) Conservative constituent #2 in
(continued) headwaters flow
*
*
Card Types 12 and 13 should be omitted if there are no point
sources in the reach (NPOINT=0)
*
*
*
* *
*************************************************
12 1-10 F10.0 QPOINT(N)
11-20 F10.0 RMP(N)
21-30 F10.0 XNAME(N)
Point source flow, cfs
River mile of point source, N
Alphanumeric description of
point source
1-8
-------
Card Group II (continued)
Card Variable
Type Columns Format Name
Description
13 1-5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
F5.0
CPOINTd.N) CBOD of point source flow,
mg/1
CPOINT(2,N) DO of point source flow, mg/1
CPOINT(3,N) Algal biomass of point source
flow, mg/1 C
CPOINT(4,N) Organic nitrogen in point
source flow, mg/1
CPOINT(5,N) Ammonia-nitrogen in point
source flow, mg/1
CPOINT(6,N) Nitrate-nitrogen in point
source flow, mg/1
CPOINT(7,N) Organic phosphorus in point
source flow, mg/1
CPOINT(8,N) Orthophosphate-phosphorus
in point source flow, mg/1
CPOINT(9,N) Temperature of point source
flow, <>C
CPOINT(10,N) Coliform bacteria in point
source flow, # of coliforms/100
ml
CPOINTdl.N) Conservative constituent #1 in
point source flow
CPOINT(12,N) Conservative constituent #2 in
point source flow
*
*
*
*
Repeat Card Types 12 and 13 NPOINT times in each reach.
NPOINT is read as Card Type 4
*
*
*
*
*************************************************
1-9
-------
Card Group II (continued)
Card Variable
Type Columns Format Name Description
****:*:::******:(:*#***:: *********:£******:»:*:»:******:(:****
* *
* Card Types 14 should be omitted if there are no diversions *
* in the reach (NDIV=0) *
* *
*************************************#*:<;:(:*:£*:(:*#*#
14 1-10 F10.0 QDIV(N) Quantity of ground return flow
to the reach, cubic feet/second
11-21 F10.0 RMDIV(N) River mile of diversion
**************************:*:*#********:$:*:(::(<*:£*:£:(:%:$::<:
* *
* Repeat Card Type 14 NDIV times in each reach. *
* NDIV is defined on Card Type 4 *
* *
* Repeat Card Types 4- 14 REACHX times. *
* REACHX is defined on Card Type 3 *
* *
****************#****#*****#******#**************
Meteorological Data
The heat budget method is used to simulate water temperature. The
data needed for the heat budget of each meteorological province must be
stored in a binary file. The required data and the order which they should
occur are given in Table A.I. A set of data is required for each of the
WOBSPD (defined on Card Type 2, above) periods per day, beginning on
DAY1 (defined on Card Type 2, above) and continuing for a total number of
days as defined by the input variable, DAY2 (defined on Card Type 2, above).
I-10
-------
Table A.I Meteorological variables for heat budget estimates and order
in which they must occur on the binary storage file
QNS Net solar radiation
kcal/meter2/second
QNA Net atmospheric radiation
kcal/meter2/second
DBT Dry bulb temperature
°C
WIND Wind speed
meters/second
PF 6.41xlO'4 * (Air pressure)
CO-1
EA Water vapor pressure
mb
PHOTO Photo period
Fraction of days
Ml
-------
Appendix II
Listing of FORTRAN 77 Source Code
-------
PROGRAM RBM10
C
C Dynamic river basin model for simulating water quality in
C branching river systems with freely-flowing river segments,
C river-run reservoirs and stratified reservoirs. Documentation
C is given in EPA 910/9-91-019.
C
C John Yearsley
C EPA Region 10 ES-098
C 1200 Sixth Ave
C Seattle, WA 98101
C (206)553-1532
C
CHARACTER*30 NAMEI
INCLUDE :RBM10.COM
C
C Open file containing reach data
C
WRITE(*,2600)
CALL FNAME(NAMEI)
OPEN(UNIT=4,FILE=NAMEI,STATUS='OLD1)
C
C Open file for output
C
WRITE(*,2700)
CALL FNAME(NAMEI)
OPEN(UNIT=7,FILE=NAMEI,STATUS='NEW)
C
C Call systems programs to get started
C
CALL BEGIN
CALL SYSTEM
C
C Close file after simulation is complete
C
CLOSE(UNIT=4)
CLOSE(UNIT=7)
1500 FORMAT(SOAl)
1600 FORMAT(SFIO.O)
2600 FORMATC Name of file containing river reach data: ')
2700 FORMATC Name of output data file: ')
STOP
END
II-l
-------
PROGRAM RBM10
C
C Dynamic river basin model for simulating water quality in
C branching river systems with freely-flowing river segments,
C river-run reservoirs and stratified reservoirs. Documentation
C is given in EPA 910/9-91-019.
C Modified for Macintosh Classic II on October 1, 1992
C For additional information contact:
C
C John Yearsley
C EPA Region 10 ES-098
C 1200 Sixth Ave
C Seattle, WA 98101
C (206)553-1532
C
CHARACTER*30 NAMEI
INCLUDE :RBM10.COM
C
C Open file containing reach data
C
WRITE(*,2600)
CALL FNAME(NAMEI)
OPEN(UNIT=4,FILE=NAMEI,STATUS='OLD')
C
C Open file for output
C
WRITE(*,2700)
CALL FNAME(NAMEI)
OPEN(UNIT=7,FILE=NAMEI,STATUS='NEW')
C
C Call systems programs to get started
C
CALL BEGIN
CALL SYSTEM
C
C Close file after simulation is complete
C .
CLOSE(UNIT=4)
CLOSE(UNIT=7)
1500 FORMAT(SOAl)
1600 FORMAT(SFIO.O)
2600 FORMATC Name of file containing river reach data: ')
2700 FORMATC Name of output data file: ')
STOP
END
II-1
-------
SUBROUTINE BEGIN
CHARACTER END*3,NAMEI*30,DLIM*3
REAL*4 LAT,DDATA(7)
C
INCLUDE :RBM10.COM
CHARACTER* 1 EXT(IO)
CHARACTER*!! PFILE
CHARACTER* 12 PPFILE
CHARACTER*20 BLANK
DATA DLIM/'END7,PFILE/'RIVPLOT.DAT7
DATA EXT/'O'.'l'.^'.'S'.'^.'S'.'G'.'T'.'S'.'g1/
DATA BLANK/1 7
C
C Initialize arrays of dimension 10
C
DO 9 N= 1,10
PNAME(N)=BLANK
HDNAME(N)=BLANK
RSNAME(N)=BLANK
WPNAME(N)=BLANK
9 CONTINUE
C
C Initialize arrays and constants
C
DO19N=1,10
NINJ(N)=0
NPLOT(N)=0
QHEAD(N)=0.
QPOINT(N)=0.
QDIV(N)=0.0
RMP(N)=-100.
RMDIV(N)=-100.
19 CONTINUE
C
C Initialize rate constants and reach name
C
DO 29 N= 1,15
HEAD(N)=.FALSE.
NPOINT(N)=0
NDPNT(N)=.FALSE.
RNAME(N)=BLANK
NDIV(N)=0
QRET(N)=0.
NRESRV(N)=0
XKBACT(N)=1.0E-10
XKDO(N)=1.0E-10
XKBODL(N)=1.0E-10
XKN44(N)=1.0E-10
XKN55(N)=1.0E-10
XKN66(N)=1.0E-10
XKP77(N)=1.0E-10
29 CONTINUE
II-2
-------
D039N=1,12
DO 39 NN= 1,100
DO39NNN=1,2
CONC(N,NN,NNN)=0.0
39 CONTINUE
NPONT=0
NDIVRS=0
NRES=0
IHEAD=0
IRES=0
IWPROV=0
C
C Card Group I
C
C Card Type 1. Alphanumeric information for title
C
READ(4,1020)XTITLE
C
C Card Type 2. Simulation time interval, starting day, number of days
C to be simulated,number of meteorological observations
C per day.
C
READ(4,1040) DT,DAY1,DAY2,WOBSPD,DAYPRT
C
C Read number of reaches, average latitude,
C maximum number of elements
C in any stratified reservoir and minimum thickness for the the surface
C element of a stratified reservoir.
C
READ(4,1040)REACHX,LAT,ZLOW,ZPLOT
C
C Change floating point constants to integers
C
IPLOT=ZPLOT
NWPD=WOBSPD
NCONST=12
NDAYS=DAYSX
NDPRNT=DAYPRT
LDAY1=DAY1
LDAY2=DAY2
PD=1./DT
DT=86400.*DT
DT2=DT/2.
NPD=PD
C
C Determine period of weather observations in terms of number of
C simulations per day
C
NWMOD=NPD/NWPD
NREACH=REACHX
II-3
-------
c
C Convert DT from fraction of days to seconds
C
C
C Check to see if plot output has been requested. If so, read
C reach numbers for which there will be plotter output and
C open RIVPLOT.DAT for output
C
IF(IPLOT.EQ.O) GO TO 55
READ(4,1044) (NPLOT(I),I=1,IPLOT)
DO49IP=1,IPLOT
NFILE=19+IP
PPFILE=PFILE//EXT(IP)
OPEN(UNIT=NFILE,FILE=PPFILE,STATUS='NEW)
49 CONTINUE
55 CONTINUE
C
C Card Group Ha. Oxygen :carbon ratio, carbon:chlorophyll a ratio,
C nitrogen :caron ratio, phosphophorus ratio,
C Michaelis-Menton term for N and P, algal preference
C for ammonia, optimal light, plankton settling rate,
C optimal, upper and lower temperatures for algal growth,
C maximum algal growth rate, maximum algal respiration rate
C
READ(4,1040) OCRAT,CCLRAT,NCRAT,PCRAT,KMN,KMP,NH3PRF
,QOPT,WSINK,TOPT,TLO,TUP,PGO,PRESO
C
C Convert meters/day to feet/second
C
VSINK=3.2808*WSINK/86400.
C
C Card Group lib. Reach characteristics
C
DO499N=1,NREACH
C
C Card Type 3. Reach description, begin and end river mile, elevation
C
READ(4,1050) RNAME(N),RMILEl(N),RMILE2(N),ELEV(N)
C
C Card Type 4. Headwater ID #, # of point sources, # of diversions,
C junction ID #, reservoir ID #, # of computational elements,
C weather province ID #.
C
READ(4,1044) NHEAD,NPOINT(N),NDIV(N))NJUNC,NRESRV(N),NCELM(N)
,NWPROV(N)
IF(NWPROV(N).GT.IWPROV) IWPROV=NWPROV(N)
C
C Card Type 5. Rate constants - XKDO,XKBODL,XKN44,XKN55,XKP77,SOD,EXCO
C
READ(4,1048) XKDO(N),XKBODL(N),XKN44(N),XKN55(N),XKP77(N),
XKBACT(N),SOD(N),EXCO(N)
II-4
-------
c
c
C Card Type 7a. River reaches: Depth and velocity coefficients,
C initialize segment volume.
C
IF(NRESRV(N).EQ.O) THEN
READ(4,1048) DEPTH1(N),DEPTH2(N),VEL1(N),VEL2(N),DEPTHO,WIDTHO
VOL(N)=DEPTHO*WIDTHO*(RMILElfN)-RMILE2(N))*5280.
C
C ***CARD TYPE 7b. Reservoir reaches. Layer thickness, bottom depth
C and coefficients for estimating reservoir geometry
C
ELSE
XKDO(N)=-10.0
NV=NRESRV(N)
RSNAME(NV)=RNAME(N)
IRES=IRES+1
READ(4,1048) Z(NV,1),ZLAYER,AVC(NV),BVC(NV),ZOUT,QQRES
ZSURFCNV, 1)=ELEV(N)-Z(NV, 1)
IOUT=((ZOUT-Z(NV, 1))/ZLAYER)+1
QRES(IOUT,NV)=QQRES
AAC(NV)=BVC(NV)
BAC(NV)=BVC(NV)
C
C Establish initial reservoir volume
C
VRES(NV,1)=AVC(NV)+ZSURF(NV,1)*BVC(NV)
ZHIGH(NV)=ZLOW+ZLAYER
NFIX=NCELM(N)
D079NF=1,NFIX+1
F=NF
Z(NV,NF)=(F-1. )*ZLAYER
ASURF(NV,NF)=AAC(NV)
79 CONTINUE
VOLO=AVC(NV)
VSEG(NV,1)=VOLO+ZLAYER*BVC(NV)
JSURF(NV)=1
DO 89 NF=2,NFIX
IF(Z(NV,NF).LT.ZSURF(NV,1)) JSURF(NV)=NF
VSEG(NV,NF)=ZLAYER*BVC(NV)
89 CONTINUE
END IF
C
C Card Types 8 and 9. Groundwater return quantity and quality.
C
READ(4,1065)QRET(N),(CRET(I,N),I=1,12)
II-5
-------
c
C Check to see if this is a headwater reach
C
IF(NHEAD.EQ.O) GO TO 100
HEAD(N)=.TRUE.
IHEAD=IHEAD+1
NMHEAD(IHEADj=NHEAD
C
C Card Types 10 and 11. Headwater quantity and quality.
READ(4,1063) QHEAD(IHEAD),HDNAME(IHEAD),(CHEAD(I,IHEAD),I=1,12)
C
100 CONTINUE
C
C Check to see if there are point sources in the reach.
C
IF(NPOINT(N).EQ.O) GO TO 150
NCYCLE=NPOINT(N)
DO139NN=1,NCYCLE
NPONT=NPONT+1
NRCH(NPONT)=N
C
C Card Types 12 and 13. Point source quantity and quality.
C
READ(4,1060) QPOINT(NPONT),RMP(NPONT),PNAME(NPONT),
+ (CPOINT(I,NPONT),I=1,12)
139 CONTINUE
150 CONTINUE
NCYCLE=NDIV(N)
C
C Check for diversions
C
IF (NDIV(N).EQ.O) GO TO 180
DO 159NN=1,NCYCLE
NDIVRS=NDIVRS+1
C
C Card Type 14. Diversion quantity and river mile of diversion.
READ(4,1048) QDIV(NDIVRS),RMDIV(NDIVRS)
C
159 CONTINUE
180 CONTINUE
C
C Check for stream junction. If NJUNC.NE.O set junction traps
C
IF (NJUNC.EQ.O) GO TO 250
NDPNT(N)=.TRUE.
NJNCTN(N)=NJUNC
NINJ(NJUNC)=NINJ(NJUNC)+1
250 CONTINUE
C
C Card Type 15. Delimiter card.
C
READ(4,1080) END
II-6
-------
c
C Checking for card sequence error. If there is, terminate program
C with diagnostic identifying reach # with error
C
IF(END.EQ.DLIM) GO TO 499
WRITE(*,3000) N
499 CONTINUE
800 CONTINUE
DO899I=1,IWPROV
NWTAPE=50+I
WRITE(*,2500) I
CALL FNAME(NAMEI)
OPEN(UNIT=NWTAPE,FILE=NAMEI)
READ(NWTAPE,1400) WPNAME(I)
DO899II=1,LDAY1-1
READ(NWT APE, 1500) LL,(DDATA(JJ,J=1,7)
899 CONTINUE
C
C Call to output routine to write system information
C
CALL WRITEO
1020 FORMAT(ASO)
1040 FORMATU6F5.0)
1042 FORMAT(8I10)
1044 FORMATU6I5)
1048 FORMAT(SFIO.O)
1050 FORMAT((A20,12F5.0))
1060 FORMAT(2F10.0,A20/16F5.0)
1063 FORMAT(F10.0,A20/16F5.0)
1065 FORMAT(F10.0/16F5.0)
1080 FORMAT(AS)
1145 FORMAT(8F10.2)
1152 FORMAT(6I3)
1400 FORMAT(A20)
1500 FORMAT(I5,7F10.0)
2500 FORMATC Energy budget file for meteorologic province - ',15)
3000 FORMATdHO,' Card sequence error for data in Reach # - ',15)
C
C Return to RMAIN
/< ******************************************************
C
RETURN
900 END
II-7
-------
SUBROUTINE SYSTEM
DIMENSION CONCJ(12,10),NINJA(10),QNJ(10),WDATA(5,7),EDATA(7)
EQUIVALENCE (EDATA(1),QNS)
C
INCLUDE :RBM10.COM
C
nl=l 5/13/92
n2=2 5/13/92
DO 999 ND=LDAY1,LDAY2
WRITE(*,*) ' DAY =',ND
DO999NDD=1,NPD
C
C Read weather data from files if time period is correct
C
IMOD=MOD(NDD,NWMOD)
IFdMOD.EQ.DTHEN
DO9I=1,IWPROV
NWR=50+I
READ(NWR,1500) LDUMM,(WDATA(I,II),II=1,7)
1700 FORMAT(I5,7E11.3)
9 CONTINUE
END IF
C
C Begin reach computations
C
IND=ND
IPD=NDD
DAY=ND
20IH1=1
QSUM=0.0
IHEAD=1
NPONT=1
NDIVRS=1
DO 39 1=1,12
C1(I)=CHEAD(I,IHEAD)
39 CONTINUE
DO 49 11= 1,10
NINJA(II)=0
QNJ(II)=0.0
DO 49 1=1,12
CONCJ(I,II)=0.0
49 CONTINUE
DAY=ND
QSUM=QHEAD(IHEAD)
C
C Begin cycling through the reaches.
C
NRR=0
II-8
-------
DO899N=1,NREACH
C
C Read meteorological data from the appropriate file
C
IWR=NWPROV(N)
C QNS=WDATA(IWR,1)
C QNA=WDATA(IWR,2)
C DBT=WDATA(IWR,3)
C WIND=WDATA(IWR,4)
C PF=WDATA(IWR,5)
C EA=WDATA(IWR,6)
C PHOTO=WDATA(IWR,7)
DO 59 1=1,7
EDATA(I)=WDATA(IWR,I)
59 CONTINUE
NR=N
RM1=RMILE1(N)
RM2=RMILE2(N)
C
C Check for reservoir. If NRESRV(N).NE.O set reservoir traps
C
IF (NRESRV(N).NE.O) THEN
CALL RESMOD
ELSE
CALL RIVMOD(RM1,RM2)
END IF
260 CONTINUE
C
C Check for a junction
C
IF (.NOT.NDPNT(N)) GO TO 300
DO 279 NJ= 1,10
IF (NJNCTN(N).NE.NJ) GO TO 279
QNJ(NJ)=QSUM+QNJ(NJ)
NINJA(NJ)=NINJA(NJ)+1
DO 269 U=l,12
CONCJ(IJ,NJ)=CONCJ(IJ,NJ)+C1(IJ)*QSUM
IF (NINJA(NJ).EQ.NINJ(NJ)) THEN
C 1(IJ)=CONCJ(IJ,NJ)/QNJ(NJ)
END IF
269 CONTINUE
QSUM=0.0
IF (NINJA(NJ).EQ.NINJ(NJ)) QSUM=QNJ(NJ)
279 CONTINUE
300 CONTINUE
II-9
-------
c
C Check for new headwaters
C
IF(HEAD(N+1))THEN
IHEAD=IHEAD+1
QSUM=QHEAD(IHEAD)
DO 399 1=1, 12
Cl(I)=CHEAD(I,IHEAD)
399 CONTINUE
END IF
450 CONTINUE
899 CONTINUE
ntmp=nl 5/13/92
nl=n2 5/13/92
n2=ntmp 5/13/92
999 CONTINUE
1500 FORMAT(I5,7F10.0)
2600 FORMATC 1615)
C
Return to RMAIN
****** ********** *********** +c^**H<^
C
950 RETURN
END
11-10
-------
SUBROUTINE QUALTY(TIME,CT,NLL,NSURF)
REAL*4IAJIS,KB,KT,K1,K2,K44,K55,K66,K77,KPI
DIMENSION CT(12),C(12),DCDT( 12)
INCLUDE :RBM10.COM
DATA ONRAT/3.42857/
C
C Light limitation function
C
FLIGHT(F,Z1,Z2,IA,IS,GAMMA)=
. (2.718*F/(GAMMA*(Z2-Z1)))*
. (EXP((-IA*EXP(-GAMMA*Z2))/IS)
. -EXP((-IA*EXP(-GAMMA*Z1))/IS))
C
C Nutrient limitation function
C
FNUTR(Y,HALF)=Y/(Y+HALF)
C
C Initialize concentrations
C
D049N=1,NCONST
C(N)=CT(N)
DCDT(N)=0.0
49 CONTINUE
DOFAC=C(1)/(0.5+C(D)
C
DO899NRNG=1,2
C
C Increment volume in second step of Runge-Kutta method
C
FCTR=0.5*(NRNG-1)
V=VELM+FCTR* DVOL
DVDT=DVOL/DT
QFCTR=QNOUT+DVDT
D=DEPTH/3.2808
C
C Compute typical temperature factors for various biological
C processes
C
T=C(9)
TM20=T-20.0
TF45=1.045**TM20
TF84=1.084**TM20
C
C Calculate rate constants which are in feedback loop
C from nutrients to algae to nutrients
C
DZQ=Z2-Z1
DIN=C(5)+C(6)
DIP=C(8)
IA=VLIGHT
IS=QOPT
11-11
-------
IF(T.GT.TOPT) GO TO 130
TLIM=EXP(-2.3*((TOPT-T)/(TOPT-TLO))**2)
GO TO 140
130 CONTINUE
TLIM=EXP(-2.3*((TOPT-T)/(TOPT-TUP))**2)
140 CONTINUE
NLIM=FNUTR(DIN,KMN)
PLIM=FNUTR(DIP,KMP)
QLIM=FLIGHT(PHOTO,Z1,Z2,IA,IS,GAMMA)
XLIM=QLIM
KLIM=1
IF(XLIM.LT.NLIM.AND.XLIM.LT.PLIM) GO TO 144
KLIM=2
XLIM=NLIM
144 IF(XLIM.LT.PLIM) GO TO 148
KLIM=3
148 CONTINUE
PG=TLIM*QLIM*NLIM*PLIM*PGO/86400.
C
C
PRES=PRESO*TF45/86400.
KPI=PG-PRES
C
C Calculate algal concentration - mg/1 of carbon
C
DCDT(3)=KPI*C(3)+(PLOAD(3)-C(3)*QFCTR)/V
C
C Calculate organic nitrogen
C
150 CONTINUE
C
C Temperature factor for Organic-N mineralization
C
TF4=TF84
IF(TF4.LT.1.0E-5) TF4=1.0E-5
K44=TF4*XKN44(NR)/86400.
DCDT(4)=NCRAT*PRES*C(3)-K44*C(4)+(PLOAD(4)-C(4)*QFCTR)/V
C
C Calculate ammonia-nitrogen
C
200 CONTINUE
C
C Temperature factor for NH4-N nitrification
C
TF5=TF45
IFOT5.LT.1.0E-5) TF5=1.0E-5
K55=TF5*XKN55(NR)/86400.
DCDT(5)=-NCRAT*NH3PRF*PG*C(3)+K44*C(4)-K55*C(5)
+(PLOAD(5)-C(5)*QFCTR)/V
II-12
-------
c
C Calculate nitrate-nitrogen
C
250 CONTINUE
C
C Temperature factor for NO3-N denitrification
C
TF6=TF45
IFCTF6.LT.1.0E-5) TF6=1.0E-5
K66=TF6*XKN66(NR)/86400.
DCDT(6)=-NCRAT*(1.-NH3PRF)*PG*C(3)+K55*C(5)
-K66*C(6)+(PLOAD(6)-C(6)*QFCTR)/V
C
C Calculate organic phosphorus
C
300 CONTINUE
C
C Temperature factor for Organic-P mineralization
C
TF7=TF84
IFCTF7.LT.1.0E-5) TF7=1.0E-5
K77=TF7*XKP77(NR)/86400.
DCDT(7)=PCRAT*PRES*C(3)-K77*C(7)+(PLOAD(7)-C(7)*QFCTR)/V
C
C Calculate inorganic phosphorus
C
350 CONTINUE
C
DCDT(8)=-PCRAT*PG*C(3)+K77*C(7)+(PLOAD(8)-C(8)*QFCTR)/V
C
C Calculate carbonaceous BOD
C
400 CONTINUE
C
C BOD
C
C Temperature factor for BOD deoxygenation
C
TF1=TF45
IF(TFl.LT.l.OE-S) TFl=1.0E-5
K1=TF1*XKBODL(NR)/86400.
DCDT(1)=-K1*C(1)+(PLOAD(1)-C(1)*QFCTR)/V
C
C Calculate dissolved oxygen
C
500 CONTINUE
C
C Temperature factor for DO rearation
C
TF2=1.024**TM20
IFCTF2.LT.1.0E-5) TF2=1.0E-5
11-13
-------
c
C Saturation level
C
CSAT=(14.62-0.3898*T+0.006969*T**2-5.897E-5*T**3)
CSAT=CSAT*((1.-(6.97E-6*ELEV(NR)))**5.167)
SDMND=TF45*SOD(NR)/(DEPTH*86400.)
C
IF (XKDO(NR).GT.O.O) GO TO 519
ZPOINT=ABS(XKDO(NR)j
IPOINT=ZPOINT
REARC=0.0
IF(IPOINT.EQ.O) GO TO 520
505 CONTINUE
GOTO (511,513,515,518),IPOINT
C
C Wind-driven effects on reaeartion in lakes and reservoirs
C
511 CONTINUE
REARC=(0.64+0.128*WIND*WIND)*3.2808/DZQ
GO TO 520
C
C Churchill-Elmore-Buckingham equation for reaeration
C
513 REARC=11.6*(U**0.969)/(DEPTH**1.673)
GO TO 520
C
C O'Connor-Dobbins equation for reaeration
C
515 REARC=12.9*U**0.5/DEPTH**1.5
GO TO 520
C
C Owens-Edwards-Gibbs equation for reaeration
C
518 REARC=21.6*((U**0.67)/(DEPTH**1.85))
GO TO 520
C
C User-defined reaeration rate
C
519 CONTINUE
REARC=XKDO(NR)
C
520 K2=TF2*REARC/86400.
C
C OCRAT and ONRAT are stoichiometric ratios for oxygen produced by
C carbon fixation and nitrate uptake. Defined in DATA statement in
C RBM10.COM
C
O2PROD=PG*C(3)*(OCRAT+ONRAT*NCRAT*(1.-NH3PRF))
02LOSS=OCRAT*PRES*C(3)
11-14
-------
c
DCDT(2)=-K2*(C(2)-CSAT)-K1*C(1)-4.57*K55*C(5)
. +02PROD-02LOSS-SDMND+(PLOAD(2)-C(2)*QFCTR)/V
C
C Temperature calculation
C R= Rho * Cp * Conversion Factor= 1000. * 1.0 / 3.2808
C (Converts energy budget from MKS units to English units)
C Initialized in DATA statement above
C
DCDT(9)=(PLOAD(9)-C(9)*QFCTR)/V
600 CONTINUE
C
C
C Calculate coliform concentrations
C
C
C Temperature-dependent rate constant
C
KB=XKBACT(NR)*2.3E-6*(1.+0.111*T)
DCDT(10)=-KB*C(10)
. +(PLOAD(10)-C(10)*QFCTR)/V
IF(NRNG.EQ.1)THEN
DO849N=1,NCONST
C(N)=C(N)+DT2*DCDT(N)
849 CONTINUE
END IF
899 CONTINUE
DO 949 N=1,NCONST
CT(N)=CT(N)+DT*DCDT(N)
IF(CT(N).LT.O.O) CT(N)=0.0
949 CONTINUE
999 CONTINUE
C
(_/
C Return to Subroutine RESMOD/RIVMOD
Q ********************************** ********************
c
RETURN
END
C
11-15
-------
SUBROUTINE WRITEO(XARG)
C
CHARACTER*! CMMA,LIM(3)
INCLUDE :RBM10.COM
DATA CMMA/Y/.LIM/'LYN'.'P1/
C
C Print general information regarding river system
C
WRITE(7,2010)XTITLE
C
C General systems parameters
C
WRITE(7,2015) NREACH,DT,LDAYl,LDAY2,NWPD,NDPRNT
C
C Headwaters
C
WRITE(7,2020) (NMHEAD(I),HDNAME(I),I=1,IHEAD)
C
C Point sources
C
WRITE(7,2025) (I,PNAME(I),NRCH(I),I=1,NPONT)
C
C Reservoirs
C
WRITE(7,2030)(I,RSNAME(I),I=1,IRES)
C
C Meteorologic provinces
C
WRITE(7,2035) (I,WPNAME(I),I=1,IWPROV)
C
C Parameters for phytoplankton dynamics
C
WRITE(7,2040) OCRAT,CCLRAT,NCRAT,PCRAT
,KMN)KMP,NH3PRF
,QOPT,WSINK
,TOPT,TLO,TUP
,PGO,PRESO
RETURN
C
C First entry point from river/reservoir modules
C
ENTRY WRITE l(XARG)
JPLOT=0
DO 19 I=1,IPLOT
IF(NR.NE.NPLOTd)) GO TO 19
IPFILE=19+I
JPLOT=I
19 CONTINUE
RM1=RMILE1(NR)
RM2=RMILE2(NR)
WRITE(7,2045) XTITLE
WRITE(7,2050) NR,RNAME(NR),RM1,RM2,ND
11-16
-------
IFONOT.HEAD(NR)) GO TO 30
WRITE(7,2060) QHEAD(IHEAD),(CHEAD(I,IHEAD),I=1,10)
30 CONTINUE
WRITE(7,2080) QRET(NR),(CRET(I,XR),I=1,10)
N'CYCLE=NPOINT(XR)
NPNTUO
. NDV1=0
NR1=NR-1
IF(NR1.EQ.O)GOTO50
D049I=1,NR1
NPNT1=NPNT1+NPOINT(I)
NDV1=NDV1+NDIV(I)
49 CONTINUE
50 CONTINUE
C
C Write titles for point sources
C
IF(NCYCLE.EQ.O) GO TO 100
WRITE(7,2052)
DO99I=1,NCYCLE
NPNT1=NPNT1+1
XBOD=CPOINT(1,NPNT1)*5.4*QPOINT(NPNT1)
XNH3=CPOINT(5,NPNT1)*5.4*QPOINT(NPNT1)
XN23=CPOINT(6,NPNT1)*5.4*QPOINT(NPNT1)
XPO4=CPOINT(8,NPNT1)*5.4*QPOINT(NPNT1)
WRITE(7,1053) PNAME(NPNT1),RMP(NPNT1),QPOINT(NPNT1)
,CPOINT(2,NPNT1),XBOD,XNH3,XN23,XPO4
,CPOINT(9,NPNT1),CPOINT( 10.NPNT1)
99 CONTINUE
100 CONTINUE
C
C Write titles for diversions
C
NCYCLE=NDIV(NR)
IF(NCYCLE.EQ.O) GO TO 200
WRITE(7,2054)
DO 199 I=1,NCYCLE
NDV1=NDV1+1
QDV1=QDIV(NDV1)
WRITE(7,1055) RMDIV(NDV1),QDIV(NDV1)
199 CONTINUE
200 CONTINUE
11-17
-------
IF(NRESRV(NR).EQ.O) THEN
WRITE(7,2056) DEPTHl(NR),DEPTH2(NR),REARC,XKBODL(NR),QNS
. ,VEL1(NR),VEL2(NR),XKN44(NR),PHOTO,XKN55(NR),XKN66(NR)
. ,XKP77(NR)
GO TO 250
ELSE
NV=NRESRV(NR)
WRITE(7,2057)REARC,XKBODL(NR),QNS
. ,XKN44(NR),PHOTO,XKN55(NR),XKN66(NR),XKP77(NR)
WRITE(7,2058) NV,QOUT,ZSURF(NV, 1),VRES(NV,2),ZSURF(NV,2)
END IF
250 CONTINUE
WRITE(7,2045) XTITLE
WRITE(7,2050) NR,RNAME(NR),RM1,RM2,ND
C
C Call first entry point in output routine. Pass argument with
C information about hydrologic regime.
C
IF(XARG.GE.O.O) THEN
WRITE(7,2100)
ELSE
WRITE(7,2110)
END IF
C
C First return point
C
C
r< ****************** ************************************
C Return to Subroutine RESMOD/RIVMOD
£J ******************************************************
C
RETURN
C
C Entry point to write simulated values of water quality and water budget
C
ENTRY WRITE2CXARG)
D0269N=1,12
COUT(N)=CONC(N,NRR,nl) 5/13/92
269 CONTINUE
C
XSAT=100.*COUT(2)/CSAT
XALGAE= 1000.*COUT(3)/CCLRAT
PKA=0.0902+(2730./(COUT(9)+273.2))
PERN=1./(1.+10.**(PKA-COUT(11)))
UNNH3=1000.*COUT(5)*PERN
XPORG=1000.*COUT(7)
XPO4=1000.*COUT(8)
XPTOT=XPO4+XPORG
ISAT=XSAT
11-18
-------
IF(XARG.GE.O.O) THEN
XAVE=XARG
WRITE(7,2400) XAVE,COUT(2)IISAT1COUT(1),XALGAE,LIM(KLIM),COUT(5)
. .COUT(6),XP04,COUT(9),COUT(10),Q1.QTPNT,QRETRN,QTDIV,DVDT,Q2
ELSE
XAVE=-XARG
QV1= QVRT1
QV2=-QVRT2
WRITE(7,2410)XAVE,COUT(2)1ISAT,COUT(1),XALGAE,LIM(KLIM)
. 1COUT(5),COUT(6),XPO4,COUT(9),COUT(10),QRINN,QV1,QV2,DVDT,QROUTN
350 CONTINUE
END IF
IF(JPLOT.EQ.O) GO TO 380
IFdPD.EQ.l.AND.NR.EQ.NPLOT(JPLOT)) then
BOD5=COUT(2)*(1-EXP(-5*XKBODL(NR)))
C
C WRITE PLOTTER OUTPUT TO RIVPLOT.DAT. RECORD CONTAINS BOD(S-DAY),
C DO, Chlorophyll a,AMMONIA NITROGEN.NITRITE + NITRATE NITROGEN,
C DISSOLVED ORTHOPHORUS, AND TEMPERATURE.
C
WRITE(IPFILE,1010) ND,CMMA,COUT(2),CMMA,BOD5,CMMA,XALGAE
. ,CMMA,COUT(5),CMMA,COUT(6),CMMA,COUT(7),CMMA,COUT(9)
end if
380 CONTINUE
1010 FORMAT(I5,3(A1,F6.1),3(A1,F6.3))A1,F6.1)
1025 FORMAT(16I5)
1050FORMAT(7X,I3,A20,7X,F7.0,8X,F7.0,8X,F7.0)
1052 FORMAT(40I2)
1053 FORMAT(T12,A20,T34,F6.1,1X,F7.1,4X,F4.1,1X,F6.0,1X,F6.0,
. 2X,F6.0,1X,F6.0,1X,F6.0,1X,F6.0,3X,F4.1)
1054 FORMAT( 16X.F6.1.8X.F6.1.6X.F6.1.6X.F6.1)
1055 FORMAT(T12,F6.1,T25,F6.1)
1056 FORMAT(14X,F6.3,8X,F6.3,8X,F6.3,8X,F6.3,8X,F6.3,4X)
1070 FORMAT(34X,F6.3,7X,F6.3,A20/
.(14X,F6.3,6X,F6.3,6X,F6.3,6X>F6.3,6X,F6.3,12X))
1080FORMAT(14X,F6.1,7X,F6.1,8X,I6,7X,F6.1)
2010 FORMAT(1H1//T12
.,' ',33X,'RIVER BASIN MODEL'
,32X,' '//
T12,' ',A80,' ')
2015 FORMATC//T12,' NUMBER OF REACHES - ',157
T12,1 TIME INCREMENT - ',F5.0,' seconds'/
T12,' STARTING DAY - ',157
T12,1 ENDING DAY - ',157
T12,1 WEATHER DATA INCREMENT - ',15,' per day'/
T12,1 PRINTOUT INTERVAL - ',15,' days')
2020 FORMAT(//T12,' HEADWATERS1/
T12,' # NAME'/
T12,' 7
T12.I5.A20)
11-19
-------
2025 FORMAT(//T12,' POINT SOURCES'/
T12,1 # NAME REACH NO.7
T12,' 7
T12,I5,A20,I5)
2030 FORMAT(//T12,' RESERVOIRS7
T12/ # NAME 7
T12,1 'I
T12,I5,A20)
2035 FORMAT(//T12,' METEOROLOGIC PROVINCES7
T12,1 # NAME '/
T12,1 7
T12.I5A20)
2040 FORMAT(//T12,' Parameters for Algal Dynamics 7
T12,' - - 7
. T12,10:C Ratio - ',F5.2/
T12,' Carbon:Chlorophyll Ratio -',F5.1/
. T12,' N:C Ratio - ',F5.3/
. T12,1 P:C Ratio - ',F5.3/
T12,' Nitrogen Half-Saturation - ',F5.3,' mg/17
T12,' Phosphorus Half-Saturation - ',F5.3,' mg/17
T12,' Algal Preference for NH4 - ',F5.1/
T12,' Optimal Light Intensity - ',F5.3,' kcal/m**2/sec7
T12,' Sinking Rate - ',F5.3,' meters/day7
T12,'Optimal Temperature - \F5.1,'Deg. C7
T12,' Minimum Temperature - ',F5.1,' Deg. C7
T12,' Maximum Temperature - ',F5.1,' Deg. C7
T12,1 Maximum Growth Rate -',F5.1,'l/days/7
T12,' Maximum Respiration Rate - ',F5.2,' I/days')
2045 FORMAT(1H1////1X,A80)
2050 FORMAT(1X,'REACH NUMBER ',I2/1X,A20,/1X,
. 'R.M. '.F5.1,1 TO R.M. '.F5.1/1X,
. 'DAY -',15)
2052 FORMAT(////T12,'POINT SOURCES1/
T12,1 7
. T12,1NAME',T35,1RIVER FLOW DO BODL NH3 N02+N03',
. ' PHOS TEMP BACT7
. T36/MILE (CFS) (MG/L) (LB/D) (LB/D) (LB/D) (LB/D) (CENT)'
. ,' (/100ML)'/)
2054 FORMAT(////T12,'DIVERSIONSV
T12,' 7
. T12,'RIVER',T25,1FLOW7T12,'MILE',T25,'(CFS)7)
2056 FORMAT(////T12,'HYDAULIC COEFFICIENTS',18X,'RATE CONSTANTS(BAS',
.'E E), DAYS**-1',10X,'HEAT BUDGET PARAMETERS'/
.iix,1 '.isx,1 -'
.,10X,' - 7
. T12/DEPTH = ',F7.4,'*FLOW**',F6.4,T51,'K2(DO) = ',F6.3,
. T68/KKBOD) = 1,F6.4,T92,'Q(NET SOLAR) = '.F6.4,1 KCAL7
. T12,'VELOCITY = 1,F7.4,'*FLOW**I1F6.4,T51,' KN44 = ',F6.3
,T92,'PHOTO PERIOD = ,F5.2,' DAYS/DAY/
T5i; KN55 = ',F6.3/
T51,' KN66 = ',F6.4/
T5i; KP77 = ',F6.4)
2057 FORMAT(////T51,'RATE CONSTANTS(BAS',
11-20
-------
:E E), DAYS**-r,iox,'HEAT BUDGET PARAMETERS'/
.nx; '.isx,' '
.,iox; 7
T51,'K2(DO) = ',F6.3,
. T68/KKBOD) = '.F6A,T92,'QtNET SOLAR) = ,F6.4.' KCAL7
T5i; KN44 = ',F6.3
,T92,'PHOTO PERIOD = ',F5.2.' DAYS/DAY'/
T5i; KN55 = ',F6.3/
T51,' KN66 = ',F6.4/
T51,1 KP77 = ',F6.4)
2058 FORMAT(////T12,'RESERVOIR NUMBER ',127
. T12,' 7/T12,'OUTFLOW',9X,' = ',F10.0,
. ' FT**3/SEC.',T54,'INITIAL DEPTH = ',F5.0,' FEET7
. T12,'RESERVOIR VOLUME = '.1PE10.2,' FT**3',T54,
. 'FINAL DEPTH = ',OPF5.0,' FEET')
2060 FORMAT(////T12,'HEADWATERS', 10X/FLOW BODL DO ALGAE1,
' ORG-N NH3-N N03-N ORG-P P04-P TEMP BACT 7
T12,1- ',
T32,'(CFS) (MG/L) (MG/L) (MG/L) (MG/L) (MG/L) (MG/L)',
1 (MG/L) (MG/L)(CENT) (/100ML)1//
28X)F7.0,3X,8(F5.1,2X),1X,F4.1,2X,F8.0)
2080 FORMAT(////T12,1GROUNDWATER',9X,'FLOW DO BODL ALGAE',
' ORG-N NH3-N N03-N ORG-P PO4-P TEMP BACT 7
T12,'RETURN',14X,'(CFS) (MG/L) (MG/L) (MG/L) (MG/L)',
' (MG/L) (MG/L) (MG/L) (MG/LXCENT) (/100ML)'/
T12,'- '/
28X,F7.0,3X,8(F5.1,2X),1X,F4.1,2X,F8.0//)
C
2100 FORMAT(////1X,'RIVER DO BODL ALGAE NH3-N N02+NO3 ',
. ' PHOS',
. ' TEMP BACT INFLOW + POINT + SEEPAGE - DIVERSIONS',
. ' - DVDT = OUTFLOW71X,'MILE (MG/L)(%SAT)(MG/L) (uG/L) ',
. '(MG/L)',
. ' (MG/L) (uG/L) (CENT) (/100ML) (CFS) (CFS) (CFS)',
(CFS) (CFS) (CFS)'//)
2110 FORMAT(////1X,'ELEV DO BODL ALGAE NH3-N N02+N03 ',
. ' PHOS',
.' TEMP BACT INFLOW + VERT(N-l) + VERT(N)',
. ' - DVDT = OUTFLOW71X,'(FT) (MG/L)(%SAT)(MG/L) (uG/L) ',
. '(MG/L)',
. ' (MG/L) (uG/L) (CENT) (/100ML) (CFS) (CFS) (CFS)',
(CFS) (CFS)'//)
2400 FORMAT(1X,F5.1,1X,F4.1,2X,I3,2H% ,F6.1,1X,F6.1,A1
. ,2X,F5.2,2X,F5.2,2X,F6.1,4X,F4.1,1X,F8.0,3X
. ,F7.1,3X,F7.1,3X,F7.1,4X,F7.1,4X,F7.1,3X,F7.1)
2410 FORMAT(1X,F5.1,1X,F4.1,2X,I3,2H% ,F6.1,1X,F6.1,A1
. ,2X,F5.2,2X,F5.2,2X,F6.1,4X,F4.1,1X,F8.0,3X
. ,F7.1,3X,F7.1,3X,F7.1,4X,F7.1,4X,F7.1,3X,F7.1)
2450 FORMAT(8F10.2)
11-21
-------
C Return to Subroutine RESMOD/RIVMOD
/~* a*:******;*:;*::*:** :*::*::*::*: A*******:*:***:*:-.*:;*: ^r^***:*:****:*:;**;*:
C
RETURN
C
C Entry point to write end-of-reach values
C
ENTRY WRITES
O2L=O2LOSS*86400.
O2P=O2PROD*86400.
C
WRITE(7,2500)
WRITE(7,2700) ELEV(NR), DEPTH, U,CSAT,02P,02L,SOD(NR)
2500 FORMAT(1H1,20X,'END-OF-REACH VALUES FOR VARIOUS PARAMETERS'//)
C
2700 FORMAT(21X,'SURFACE ELEVATION - ',F8.1,' FEET7
21X,'WATER DEPTH - '.F8.1,1 FEET7
21X,'WATER VELOCITY - '.F8.4,1 FEET/SEC'/
21X/DISSOLVED OXYGEN SATURATION - '.F8.1,1 MG/L7
2 1X/OXYGEN PRODUCTION - '.F8.4,' MG/L/DAYV
21X/RESPIRATION RATE - ',F8.4,' MG/I7DAY/
21X,'SEDIMENT OXYGEN DEMAND - ',F8.2,' GM/M**2/DAY/)
C Return to Subroutine RESMOD/RIVMOD
O **!|t**##*j|c****#J)t*** #*****#*******###***#***# ***********
C
RETURN
END
11-22
-------
SUBROUTINE RESMOD
REAL*4C(12),CRES(12,10),DCDT(12),INFRED(10),QRIN(10),QROUT(10)
,QVERT(0:10),RSLOAD( 12,10)
INCLUDE .-RBM10.COM
DATA DIFF/l.OE-4/
GFUNCfA,B,EL)=A+B*EL
HFUNC(A.B,EL)=A*EXP(B*EL)
C
C Initialize important counters, constants and variables
C
NFIX=NCELM(NR)
NDVRN=NDIV(NR)
NPNN=NPOINT(NR)
NPN=NPONT
NV=NRESRV(NR)
NSURF=JSURF(NV)
NSRFM1=NSURF-1
IF(NSRFMl.LE.O) NSRFMl=l
QIN=QSUM
QOUT=0.0
D019N=1,NFIX
AXY(N)=BVC(NV)
IR=NRR+N
DO 19NN=1,NCONST
CRES(NN,N)=CONC(NN,IR,nl) 5/13/92
19 CONTINUE
D029N=1,NFIX
QRIN(N)=0.0
QROUT(N)=0.0
D029NN=1,12
RSLOAD(NN,N)=0.0
29 CONTINUE
C
C Determine level of inflow for upstream flow and update loading for
C appropriate element
C
TMPIN=CONC(9,NRR,nl) 5/13/92
CALL LAYRIN(CRES,TMPIN,NSURF,NLYR)
QRIN(NLYR)=QRIN(NLYR)+QSUM
DO 39 NNN=1,NCONST
RSLOAD(NNN,NLYR)=RSLOAD(NNN,NLYR)+QSUM*C1(NNN)
39 CONTINUE
C
C Determine level of inflow for point sources and update loading
C for appropriate element
C
IF(NPNN.EQ.O) GO TO 95
D089NN=1,NPNN
QIN=QIN+QPOINT(NPN)
TMPIN=CPOINT(9,NPN)
CALL LAYRIN(CRES,TMPIN,NSURF,NLYR)
11-23
-------
QRIN(NLYR)=QRIN(NLYR)+QPOINT(NPN)
DO49NNN=1,NCONST
RSLOAD(NNN,NLYR)=RSLOAD(NNN,NLYR)+QPOINT(NPN)*CPOINT(NNN,NPN)
49 CONTINUE
NPN=NPN+1
89 CONTINUE
95 CONTINUE
IF(NDVRN.EQ.O) GO TO 105
DO 99 NN=1,NDVRN
QOUT=QOUT+QDIV(NDIVRS)
QROUT(NSURF)=QROUT(NSURF)+QDIV(NDVR)
NDIVRS=NDIVRS
99 CONTINUE
105 CONTINUE
C
C Determine inflow layer for return flow and update loading for
C appropriate element
C
QIN=QIN+QRET(NR)
TMPIN=CRET(9,NR)
CALL LAYRIN(CRES,TMPIN,NSURF,NLYR)
QRIN(NLYR)=QRIN(NLYR)+QRET(NR)
DO 149NNN=1,NCONST
RSLOAD(NNN,NLYR)=RSLOAD(NNN,NLYR)+QRET(NR)*CRET(NNN,NR)
149 CONTINUE
QOUT=QOUT+QRES(NSURF,NV)
QDWN=QRES(NSURF,NV)
QROUT(NSURF)=QROUT(NSURF)+QRES(NSURF,NV)
DO 189 NN=1,NSRFM1
QOUT=QOUT+QRES(NN,NV)
QDWN=QDWN+QRES(NN,NV)
QROUT(NN)=QROUT(NN)+QRES(NN,NV)
189 CONTINUE
C
C Mass balance, estimate new volum and surface elevation
C
DVOL=(QIN-QOUT)*DT
VRES(NV,2)=VRES(NV, D+DVOL
ZSURF(NV,2)=(VRES(NV,2)-AVC(NV))/BVC(NV)
C
C Calculate vertical velocities from continuity
C
QVERT(1)=QRIN(1)-QROUT(1)
IF(NSRFMl.LE.l) GO TO 300
DO249NN=2,NSRFM1
QVERT(NN)=QRIN(NN)-QROUT(NN)+QVERT(NN-1)
249 CONTINUE
300 CONTINUE
QVERT(NSURF)=0.0
11-24
-------
c
C Calculate advective load using upstream weighting
C
DO349NN=1,NSRFM1
XR1=NN+1
XR2=XN
FQ=1.0
IFCQVERT(NN).LE.O.O) THEN
NR1=NN
NR2=NN+1
FQ=-1.0
END IF
DO 347 NNN=1,NCONST
RSLOAD(NNN,NR1)=RSLOAD(NNN,NR1)+FQ*QVERT(NN)*CRES(NNN,NR2)
347 CONTINUE
349 CONTINUE
C
C Calculate diffusion load
C
DO 399 N=1,NSRFM1
Nzl=N
Nz2=N+l 5/13/92
dz=z(nv,nz2)-z(nv,nzl) 5/13/92
DO 399 NN=1,NCONST 5/13/92
DFFUSN
. =DIFF*AXY(N+l)*(CRES(NN,Nz2)-CRES(NN,Nzl))/dz 5/13/92
RSLOAD(NN,N)=RSLOAD(NN,N)+DFFUSN
RSLOAD(NN,N+ 1)=RSLO AD(NN,N+ D-DFFUSN
399 CONTINUE
C
C Call Subroutine QUALTY to calculate concentrations in each element
C
INFRED(NSURF+ 1)=0.5*QNS
DZZ=ZSURF(NV, l)-Z(NV.NSURF)
AREA2=ASURF(NV,NSURF+1)
XKDO(NR)=-1.0
DO419N=NSURF,1,-1
AREA1=ASURF(NV,N)
C
C Account for sinking rates of plankton
C
RSLOAD(3,N)=RSLOAD(3,N)+VSINK*(AREA2*CRES(3>N+1)-AREA1*CRES(3,N))
C
C Light shading due to plankton
C
CHLR=1000.*CRES(3,N)/CCLRAT
GAMMA=(EXCO(NR)+0.0088*CHLR+0.054*CHLR**0.67)/3.2808
INFRED(N)=HFUNC(INFRED(N+1),GAMMA,-DZZ)
HEAT=INFRED(N+1)*AREA2-INFRED(N)*AREA1
11-25
-------
c
C Update thermal load. RFAC=Rho*Cp*Conversion=1000.*l./3.2808
C =304.8
C
RSLOAD(9,N)=RSLOAD(9,N)+HEAT/304.8
IF(N.GT.l) DZZ=(Z(NV,N)-Z(NV,N-1))
C
C' Initialize constituent loading factors
C
DO 409 NN=1,NCONST
PLOAD(NN)=RSLOAD(NN,N)
409 CONTINUE
IF(N.EQ.NSURF.OR.NSURF.EQ.UTHEN
TSURF=CRES(9,N)
CALL ENERGY(TSURF,QSURF)
VELM=GFUNC(AVC(NV),BVC(NV),ZSURF(NV,1).Z(NV,NSURF))
VSURF=VELM
IF(NSURF.EQ.l) VSURF=VSURF+VOLO
C
C Account for heat flux at surface
C
PLOAD(9)=PLOAD(9)+QSURF*AREA2/304.8
ELSE
VELM=VSEG(NV,N)
DVOL=0.0
QSURF=0.0
END IF
Z 1=ZSURF(NV, 1)-Z(NV,N+1)
Z2=ZSURF(NV, 1)-Z(NV,N)
QNOUT=QROUT(N)+AMAX1(QVERT(N),0.0)-AMIN1(QVERT(N-1),0.0)
VLIGHT=INFRED(N+1)
CALL QUALTY(TIME,CRES( 1,N),N,NSURF)
NLIMT(N)=KLIM
XKDO(NR)=0.0
ARE A2=AREA 1
419 CONTINUE
C
C Call Subroutine MIX if there is a temperature instability
C
IF(NSURF.GT.1)THEN
CALL MIX(CRES,VSURF,NSURF,NV)
END IF
DO429N=NSURF,1,-1
DO 429 NN=1,NCONST
nrrn=nrr+n
CONC(NN,NRRN,n2)=CRES(NN,N)
429 CONTINUE
NRR=NRR+NSURF
DO439N=NSURF,1,-1
KLIM=NLIMT(N)
QRINN=QRIN(N)
QROUTN=QROUT(N)
QVRT1=0.0
11-26
-------
IF(N.NE.l) QVRT1=QVERT(N-1)
QVRT2=QVERT(N)
XARG=-(Z(NV,N+1)+Z(NV,N))*0.5
NPMOD=MOD(ND,NDPRNT)
IF(NPMOD.NE.O) GO TO 435
IFdPD.EQ.DTHEN
' IF(N.EQ.NSURF) CALL WRITE IfXARG)
CALL WRITE2(XARG)
END IF
435 CONTINUE
NRR=NRR-1
439 CONTINUE
C
C Update downstream loading by computing loading released from
C reservoir
C
D0459NN=1,NCONST
CSUM=0.0
DO 449 N=1,NSURF
CSUM=CSUM+QRES(N,N"V)*conc(NN,N+nrr,n 1) 5/13/92
449 CONTINUE
C1(NN)=CSUM/QDWN
459 CONTINUE
C
C Combine surface element with next lower one if thickness is less than
C minimum
C
ZDIFF=ZSURF(NV,2)-Z(NV,NSURF)
IF(ZDIFF.LE.ZLOW) THEN
NR1=NSURF
NR2=NSURF+1
NRR1=NRR+NSURF-1
NRR2=NRR+NSURF
DVOL1=GFUNC(AVC(NV),BVC(NV),Z(NV,NR2))
-GFUNC(AVC(NV),BVC(NV),Z(NV,NR1))
DVOL2=GFUNC(AVC(NV),BVC(NV),ZSURF(NV,1))
. -GFUNC(AVC(NV),BVC(NV),Z(NV,NR2))
DO 469 NN=1,NCONST
CONC(NN,NRRl,n2)
.=(DVOLl*CONC(NN,NRRl,n2)+DVOL2*CONC(NN,NRR2,n2)) 5/13/92
./(DVOL1+DVOL2)
C
C Zero out concentration in the element which was eliminated so
C there will be no artificial fluxes from the top down
C
CONC(NN,NRR2,n2)=0.0
469 CONTINUE
NSURF=NSURF-1
END IF
11-27
-------
c
C If layer thickness exceeds tolerance, create a new element with
C concentration equal to value in the old surface layer
C
IF(ZDIFF.GE.ZHIGHfNV)) THEN
NR1=NSURF
NR2=NSURF+1
DO 479 NN=1,NCONST
CONC(NN,NR2,n2)=CONC(NN,NRl,n2) 5/13/92
479 CONTINUE
NSURF=NSURF+1
END IF
C
C Increment sub-reach index (NRR) by NFIX since that is the
C number of places reserved for vertical layers in each reservoir.
C NFIX is initialized in the input data
C
NRR=NRR+NFIX
C
C Update reservoir surface elevation and volume
C
ZSURF(NV,1)=ZSURF(NV,2)
VRESCNV, 1)=VRES(NV,2)
C
o ******************************************************
C Return to Subroutine SYSTEM
Q ******************************************************
C
RETURN
END
11-28
-------
SUBROUTINE RIVMOD(RM1,RM2)
REAL*4C(12),DCDT(12)
INCLUDE :RBM10.COM
EQUIVALENCE (NCONST,NC),(AXY(1),AREA1),(AXY(2)IAREA2)
DATANone'l/ 5/13/92
C
C Specify some constants and determine the number of
C computational elements in the reach.
C
NCYCLE=NCELM(NR)
CYCLE=NCYCLE
NSURF=1
DXM=(RMl-RM2yCYCLE
DX=5280.*DXM
QRETRN=QRET(NR)/CYCLE
XKDO(NR)=-2.0
C
C Calculate properties in each of the computational elements
C
DO249NN=1,NCYCLE
C
C Update loading rates with return flow
C
DO99NNN=1,12
PLOAD(NNN)=QRETRN*CRET(NNN,NR)
99 CONTINUE
C
C Index segment number (NRR)
C
NRR=NRR+1
C
C Initialize reach flows
C
QTDIV=0.0
QTPNT=0.0
C
C Location of upstream and downstream boundaries
C
X1=RM1-DX*(NN-1.)/5280.
X2=RM1-DX*NN/5280.
XAVE=(Xl+X2)/2.
C
C Check for diversions
C
IF (NDIV(NR).EQ.O) GO TO 150
135 XMDIV=RMDIV(NDIVRS)
IF (XMDIV.GT.X1.0R.XMDIV.LT.X2) GO TO 150
QTDIV=QDIV(NDIVRS)+QTDIV
NDIVRS=NDIVRS+1
GO TO 135
150 CONTINUE
11-29
-------
c
C Check for point sources
C
IF (NPOINT(NR).EQ.O) GO TO 165
155 XMP=RMP(NPONT)
C
C if a point source is within reach boundaries, combine
C the waste flow with the river flow.
C
IF (XMP.GT.X1.OR.XMP.LT.X2) GO TO 165
DO159NNN=1,12
PLOAD(NNN)=PLOAD(NNN)+CPOINT(NNN,NPONT)*QPOINT(NPONT)
159 CONTINUE
QTPNT=QTPNT+QPOINT(NPONT)
NPONT=NPONT+1
GO TO 155
165 CONTINUE
C
C Determine the water balance
C
Q1=QSUM
Q2=QSUM+QRETRN+QTPNT-QTDIV
QIN=QSUM+QRETRN+QTPNT
QNOUT=QIN
C
C Update loading rates with upstream input
C
DO169NNN=1,12
PLOAD(NNN)=PLOAD(NNN)+Q1*C1(NNN)
169 CONTINUE
C
C Perform mass balance and estimate volume change
C
DVOL=(Q1-Q2)*DT
VELM=VOL(NR)/CYCLE
QAVE=0.5*(Q1+Q2)
U=VEL1(NR)*QAVE**VEL2(NR)
DEPTH=DEPTH1(NR)*QAVE**DEPTH2(NR)
DO179NNN=1,2
AXY( NNN)=VE LM/D EPTH
179 CONTINUE
Z 1=0.0
Z2=DEPTH
C
C Initialize concentrations prior to calling the
C Runge-Kutta differential equation solver
C
DO 189 NNNN= 1,12
C(NNNN)=CONC(NNNN,NRR,n 1)
189 CONTINUE
11-30
-------
c
C Account for sinking rates of plankton
C
PLOAD(3)=PLOAD(3)-VSINK*AREA1*C(3)
C
C Surface exchange of thermal energy
C
TSURF=C(9)
CALL ENERGYCTSURF.QSURF)
PLOAD(9)=PLOAD(9)+AREA2*(0.5*QNS+QSURF}/304.8
VLIGHT=0.5*QNS
C
C Calculate light extinction coefficient
C
CHLR=1000.*C(3)/CCLRAT
GAMMA=EXCO(NR)+0.0088*CHLR+0.054*CHLR**0.67
C
C Call SUBROUTINE QUALTY, the Runge-Kutta
C differential equation solver
C
CALL QUALTY(TIME,C,None,None)
NPMOD=MOD(ND,NDPRNT)
IF(NPMOD.NE.O) GO TO 195
IF(IPD.EQ.1)THEN
IF(NN.EQ.l) CALL WRITEl(XAVE)
CALL WRITE2(XAVE)
END IF
195 CONTINUE
DO 199 NNNN= 1,12
CONC(NNNN,NRR,n2)=C(NNNN) 5/13/92
Cl(NNNN)=CONC(NNNN,NRR,nl) 5/13/92
199 CONTINUE
QSUM=Q2
249 CONTINUE
IF(IPD.EQ.l.AND.NPMOD.EQ.O) CALL WRITES
C
Q
C Return to Subroutine SYSTEM
Q
C
RETURN
END
11-31
-------
FUNCTION RHO(T)
RHO=1000.-(((T-3.98)**2)*(T+283./y'(503.57*(T+67.26))
RETURN
END
11-32
-------
SUBROUTINE LAYRINfCRES.TMPIX.NSURF.NLYR)
REAL*4CRES(12,10.)
INCLUDE :RBM10.COM
IFfXSURF.EQ.l; THEN
NLYR=1
RETURN
END IF
N=NSURF
10 CONTINUE
RHO1=RHO(TMPIN)
RHO2=RHO(CRES(9,N);
IF(RH01.GT.RH02) THEN
N=N-1
IF(N.EQ.1)GOT050
GO TO 10
50 CONTINUE
END IF
NLYR=N
C
Return to Subroutine RESMOD
*******#**********************************************
C
RETURN
END
11-33
-------
SUBROUTINE MIX(CRES,VSURF,NSURF,NV)
DIMENSION CRESf 12,lOXCSUMf 12)
INCLUDE :RBM10.COM
NMIX=XSURF
NSRFM1=NSURF-1
VSUM=VSURF
TMIX=CRES(9,NSURF)
RHO1=RHO(TMIX)
DO49NN=1,NCONST
CSUM(NN)=CRES(NN,NSURF)*VSURF
49 CONTINUE
DO99N=NSRFM1,1,-1
RHO2=RHO(CRES(9,N))
IF(RH01.GT.RH02) THEN
VSUM=VSUM+VSEG(NV,N)
DO 79 NN=1,NCONST
CSUM(NN)=CSUM(NN)+CRES(NN,N)*VSEG(NV,N)
79 CONTINUE
TMIX=CSUM(9)/VSUM
RHO1=RHO(TMIX)
NMIX=N
END IF
99 CONTINUE
DO 199N=NSURF,NMIX,-1
DO 199 NN=1,NCONST
CRES(NN,N)=CSUM(NN)/VSUM
199 CONTINUE
C
/-< ************************************* *****************
C Return to Subroutine RESMOD
f~< ******************************************************
C
RETURN
END
11-34
-------
SUBROUTINE ENERGY(TSURF.QSURF)
REAL*4 LVP
INCLUDE :RBM10.COM
DATA PI/3.14159/.EVRATE/1.5E-9/
EO=2.1718E8*EXP(-4157.0/(TSURF+239.09))
RB=PF*(DBT-TSURF)/(EO-EA)
LVP=597.0-0.57*TSURF
QEVAP=1000.*LVP*EVRATE*WIND*(EO-EA)
QCONV=RB*QEVAP
BOWEN=2.-SIN(2.*PI*DAY/365.+PI/6.)
QWS=6.693E-2+1.471E-3*TSURF
QSURF=0.5*QNS+QNA-QWS-QEVAP+QCONV
C
/~1 ^^^^A^^***^*^*****^*****^^*^^**^*^^^^^*^*^^^^^^*^^*^^^:
C Return to Subroutine RESMOD/RIVMOD
Q il:**^*********)!:**;****^*********** *********************!*=
C
RETURN
END
11-35
-------
SUBROUTINE FNAME(INFILE)
CHARACTER*30 INFILE
READ(*,1000) INFILE
1000 FORMAT(ASO)
C
/"I M'Jf-M*.-*-**,**^. #.***** -J£ ***%
C Return to RJVIAIN
. %^%*.^
900 RETURN
END
11-36
------- |