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

-------