UNITED STATES ENVIRONMENTAL PROTECTION AGENCY
REGION II
26 FEDERAL PLAZA
NEW YORK. NEW YORK 1OOO7
Documentation for
ES001
"A Steady-state, One Dimensional, Estuarine
Water Quality Model"
Steven C. Chapra, Sanitary Engineer
Systems Analysis Section
Seymour Gordimer, Chief
Computer Systems Section
Data Systems Branch
September, 1973
-------
UNITED STATES ENVIRONMENTAL PROTECTION AGENCY
REGION II
26 FEDERAL PLAZA
NEW YORK. NEW YORK 1OOO7
Documentation for
ES001
'A Steady-state, One Dimensional, Estuarine
Water Quality Model"
Steven C. Chapra, Sanitary Engineer
Systems Analysis Section
Seymour Gordimer, Chief
Computer Systems Section
Data Systems Branch
September, 1973
-------
Table of Contents
Introduction 2
Theory 5
Basic Theory 7
Development of Water Quality Equations
for an Estuary 12
Other Applications 28
Segmentation of Estuary Systems 29
Evaluation of Boundary Conditions 31
Sample Problem 41
Method for Evaluation and Presentation
of Boundary Conditions for
Complex Estuaries 47
Nomenclature 64
Sensitivity 67
Description of the Computer Program 80
Description of Main Program with
Flowcharts 81
Description of Subroutines 88
Description of Functions 100
User's Manual 104
Initial Steps in Formulating ES001 104
Discussion of System Parameters 106
Program Run Preparation 109
-------
Data Description and Input Requirements 113
Output 124
Restrictions 128
Example Runs 130
References 159
Appendix (ES001 listing)
1A
-------
INTRODUCTION
- 2 -
-------
In 1968, Hydroscience, Inc. prepared a report for the FWPCA entitled,
"Mathematical Models for Water Quality for the Hudson-Champlain and Metropolitan
Coastal Water Pollution Control Project."-'- The report presented the develop-
ment procedures of water quality models for various waterways including the
New York Harbor Complex, and parts of the Raritan and Hudson Rivers. As part
of the project, several computer programs were developed to model the various
water systems. Unfortunately, inadequate documentation was available for these
programs and an effort has been made by the Environmental Protection Agency to
publish such documentation as well as to modify the programs for efficiency
and ease of use. It is the purpose of this report to document one of the
programs - ES001/2.
ES001 is presently programmed for use on the IBM 370 while ES002 is
compatible with an IBM 1130 system. Thus, the designation ES001/2 is used
when referring to both programs. In general, the two programs have identical
input, output and theoretical basis. The remainder of this documentation,
unless explicitly stated, will deal with ES001 in particular with most of the
information relevant to ES002.*
In its present form, ES001 is capable of modelling general one-
dimensional estuarine systems for a variety of substance concentrations.
The program requires segmentation of the system being modelled into sec-
tions within which the various geomorphological, physical and hydro-
logical parameters of the estuary are constant. The sections are then
connected mathematically, the junction points of these segments being
boundary points where these parameters can change. Several types of
junctions are allowed including triple junctions, dams, etc., which in
combination allow the modelling of numerous types of configurations.
The model is assumed to be at steady state and to be tidally averaged.
Steady state means that the variables and parameters describing a system
do not vary with time. Investigation of steady—state concentrations is a
useful planning strategy for many situations. For example, in summer months,
usually the most critical period as far as some types of pollution are concerned,
the hydrograph of fresh water flow entering an estuary may be low and
relatively constant creating an excellent opportunity for application
of this model. Tidally averaged means that fluctuations in response caused
by tidal ascillations are not treated explicitly in the model. The effects
of tidal mixing are manifested indirectly through inclusion in a bulk
dispersion term.
Finally, reaction kinetics are assumed to be of the first-order for
the reasons that many substances reacting in an estuary approximate this
type of behavior.
*An addendum to this report includes the source deck for ES002 and a
description of its specifications.
- 3 -
-------
This report assumes that the user has a basic knowledge of stream
and estuarine analysis. The basic approach used in ES001 was developed
by O'Connor and several references may be consulted to supplement this
development (1,2,3,4).
ES001 is meant to be used to either furnish insight into a particular
phenomena, or as a predictive device for use in water quality planning.
Care should be taken at all times to consider all the assumptions
underlying its formulation and by no means could it ever be construed
to apply to any and every estuarine situation. With this in mind it is
an excellent tool for the use of those interested in applying rational
approaches to the problems of the deterioration of the environment.
- 4 -
-------
THEORY
- 5 -
-------
The Law of Conservation of Mass is the basis from which the ES001
model has been developed. Using this law and making certain assumptions,
differential equations may be derived which describe the behavior of a
mass in a one-dimensional estuary. Since the assumptions used to derive
the equations state that the body of water have constant hydraulic and
geometric characteristics, the estuary must be divided into segments in
which parameters, such as cross-sectional area, depth, width, etc., are
constant. Segments are also formed at points of waste input. Individual
differential equations are then written for each segment and solved in
terms of unknown coefficients. The segments are then joined in a mathe-
matical fashion, by means of boundary conditions relating concentration
and flux of mass at each face of the junctions between the segments, and
the unknown coefficients are determined mathematically. Finally, the
coefficients are placed in the original solutions which are then used
to generate continuous values of the parameter under study throughout
each section.
The following sections are designed to introduce the user to the
theory upon which program ES001 is based, as well as the analytic forms
in which this theory is applied in the program. The chapter is divided
into two areas: a theoretical discussion and a section on the specific
way in which the theory is applied in ES001.
The former includes a section on basic theory, a development of
water quality equations for an estuary, sections on segmentation of the
system and evaluation of boundary conditions and a sample problem. The
latter is basically a review of the evaluation and presentation of bound-
ary conditions as particular to program ES001.
In general, discussion will focus on the BOD-DO deficit system. However,
the program is capable of modeling analogous systems of what might be called
sequential reactions of two substances having first order kinetics. For
instance ES001 can be used to model a nitrogen reaction with ammonia and
nitrate analogous to BOD and D.O. deficit, respectively.
- 6 -
-------
Basic Theory:
The distribution of a substance in a natural body of water is gov-
erned by the law of conservation of mass; that is, the rate of change
with respect to time of a mass within a volume must be equal to the net
influx of mass. For a one-dimensional estuary this may be expressed in
a general way by:
3c E 3 (K 3c. 2. 9c + v c
•rrr ° — r— (A T—) - -f -r— ±ZS
dt A dx 9x A dx
where
c = concentration of some mass (e.g., BOD or DO deficit)
[M/L3]*
t = time [T]
E = dispersion coefficient [L^/T] which incorporates dif-
fusion effects caused by turbulence, velocity gradients,
density currents as well as the mixing brought about by
the tidal velocity. ES001 does not model the concentra-
tion of substances within the tidal cycle but rather
looks at the distribution of the mass as it would behave
over a longer time span while incorporating the effects
of the tides in this coefficient.
A = cross-sectional area of the estuary [L^]. O'Connor has
developed equation 1 for variable areas, A(x),6. These
solutions are not suited for this general approach and
therefore, cross-sectional areas are approximated by con-
stant values in ES001.
x = distance along the estuary [L].
Q = net advective flow of the estuary [L-*/T]. As implied
in the definition of E, this does not include the tidal
velocity. In many estuaries this is best described as
the "fresh water flow."
±ZS = the various sources and sinks of the mass in question.
For DO a sink of oxygen may be a deposit of sludge on
the streambed which exerts a demand on the oxygen in
the estuary while a source might be the oxygen which
enters the water from the atmosphere.
*Letters in brackets designate units. M is mass, L is length and T is time,
- 7 -
-------
Program ES001 generates solutions for equation 1 when It has con-
stant coefficients (e.g., area is constant) and when it is at steady
state (3c/3t = 0). Equation (1) may then be written as follows for a
conservative substance like chlorides :
0 = E j-£ - V~+ Y ............... (2)
d 2 dx c
x
where
C = concentration of conservative substance [M/L3]
U = Q/A = net advective velocity [L/T]
Y = sources and sinks of conservative substance [M/(L3T)]
= 0 (no sources or sinks of conservative substance will
be considered in ES001)
For a non-conservative substance such as BOD, equation 1 can be
rewritten as:
2
0 = E -- - U 2 - RK-L + Yu ............ (3)
d z dx b
x
where
L = concentration of BOD [M/L3].
U = Q/A = net advective velocity [L/T] .
RK = removal coefficient of BOD [1/T], This represents the
rate at which BOD is removed from the system by a number
of means (usually assimilation by the organisms in the
estuary). As can be seen the assumption is made that
removal is of first order reaction kinetics. This holds
that the rate of reaction of a substance is directly pro-
portional to the concentration of that substance at any
time, t. This assumption is justified by the fact that
first order kinetics have been found to be an excellent
approximation to a large body of reactions in natural water
systems and that the equations are more amenable to solution.
Y, = additional sources and sinks of BOD [M/(L3T)].
b
For DO deficit, equation 1 can be rewritten as:
' ~
- 8 -
-------
where
D - DO deficit [M/L3] = C - DO
s
C = saturation value of dissolved oxygen in water at a specific
temperature and pressure [M/L3]
DO = concentration of dissolved oxygen [M/L3]
AK = reaeration coefficient [1/T]. The reaction rate represents
the rate at which oxygen enters the body of water from the
atmosphere. It is assumed to follow first order kinetics.
Y, = additional sources and sinks of DO deficit [M/(L3T)1
d
Equations (2), (3) and (4) can be integrated to yield:
5 - Ze(U/E)'X + J (5)
L = Be S'X + Ce V'X + YB (6)
D = GeSDX + YD (7)
where Z, J, B, C, G and H are constants to be evaluated by consideration
of the boundary conditions.
S and V Contain the parameters U, E, RK
SD and VD Contain the parameters U, E, AK
YB = The particular integral of the BOD
source and sink terms
YD = The particular integral of the
deficit source and sink terms and
which always includes the BOD solu-
tion equation
Equations (5), (6) and (7) were,developed on the assumption of con-
stant values of the parameters U, E, RK, and AK. These equations are
applicable to any section of an estuary where the assumption of constant
parameter values is valid. The estuary must be segmented at locations
where the value of a parameter changes, at locations of point waste
inputs and at locations of changes in the sources or sinks of material.
Figure 1 illustrates such segmentation. The appropriate equation
- 9 -
-------
(Equation (5), (6) or (7), depending on the substance being studied) i«
then applied to each of the sections formed by the segmentation proc«««.
The unknown coefficients are evaluated by considering the boundary
conditions at the junction points of adjacent segments. In eff«ct, th«
boundary conditions at each junction of sections couples the several
segments into a single system interrelated through the values of the
appropriate unknown coefficients. Considering the BOD solution as an
example, the sections formed by segmentation are coupled into a singla
system by evaluation of boundary conditions at the junction p»ints be-
tween sections. The values of the unknown coefficients B and C in any
section will reflect the phenomena occurring in all other s«cti*ns »f
the system.
The total number of boundary conditions that must be evaluated in
any system is equal to twice the number of segments in the systam.
There are three classes of boundary condition equations which can b«
specified. These are:
1. Infinity conditions:
At locations remote from a waste source, the anticipated
concentration attributable to the waste source should
approach or be equal to zero. This can be stated mathe-
matically by evaluating the appropriate equation, at posi-
tive or negative infinity.
2. Concentration equality:
At the junction of two segments, the concentration calcu-
lated by the equation in either segment must be the same.
Thus, the profile equation in the upstream segment must
equal the profile equation in the downstream segment at the
junction of the two segments.
3. Mass balance:
It is possible to develop a material balance around an
element at the junction of two sections. The quantity
of the substance entering the element must equal the
quantity of substance leaving the element.
Proper application of the three general classes of boundary condi-
tions indicated above will generate sufficient equations to evaluate
all unknown coefficients.
- 10 -
-------
CHANGE IN CROSS-
SECTIONAL AREA'
DAM
WASTE
INPUT
TRIBUTARY
fCHANGE IN
'DEPTH
H
5
DISTANCE
Fi.gure 1. Illustration of segmentation.
- 11 -
-------
At every junction of "n" segments (n ran it be greater than one), a
single mass balance equation, and (n - 1) concentration equations can
be specified. If one end of a section is not bounded by adjacent seg-
ments, the infinity boundary condition may be rpplied.
Evaluation of the boundary condition equations linking adjacent
segments results in a series of simultaneous equations equal to the num-
ber of unknown coefficients. It is convenient to arrange these unknown
equations into a matrix and employ standard techniques to evaluate each
of the unknown coefficients. In order to facilitate setting up th«
matrix of unknown coefficients, a series of gereral boundary condition
equations have been developed and tabulated in this report. The boundary
equations are presented in terms of the elements of the matrix and the
forcing functions. All specific boundary condition equations required
to mathematize the study area have been included in the tabulation.
Development of Water Quality Equations for an Fstuary:
The System
Figure 2. An Estunry.
For the purpone of analysis, the estuary is defined as that portion
of the river which is under the influence o^ tidal action in which the
dispersion factor is always significant. The advective term may or may
not be significant depending on the rate of fresh-water flow and the
cross-sectional area. As noted previously, the advective term is a net
- 12 -
-------
advection or fresh v/ater flow and does not -include the tidal velocity.
Any mixing due to the tidal velocity is incJnded in the bulk dispersion
rate, E, which also takes into account the effects of turbulence, velocity
gradients and density currents.
The particular system of ES001 is one dimensional in nature. Thi»
implies that the substance being modelled is assumed to be uniformly
distributed in the Y and Z direction* as illustrated in Figure 2.
Conservative Substance
A conservative substance is one which does not decay significamtly
when entering an estuary. As such, chlorides and some pesticides might
be considered in this category.
The equations for evaluation of conservative substances __•
can be developed by considering a differential element of an estuary
as shown in Figure 3:
ditive x
Figure 3. Differential E.:.. rut of an Estuary
(conservative s. \,r an;e) ,
A mass balance is obtained over the 1:ime interval At since the law
of conservation of mass applies. The bu?"..1Jup in the element is evaluated
from a consideration of the material enteiins and leaving through advec-
tive and dispersive transport.
-------
In this particular case, no source and sinks are considered as it
is assumed that the concentration of a conservative substance is only
affected by advective and dispersive forces.
Mass Balance:
entering element:
advection: Q- C-
dispersion: -E, A, 1
1 1 ~^—
dx
leaving element:
advection: Q7 (£- + ~~~~ AX)
dispersion: -E A9 -r— (£.. -1-
Z. £. Q X 1_
Equation:
f\ ri ?"
r1 A —- (r + -—
C2 A2 8x (?1 Dx
VOL = volume
assume
Q = Q = 0 , E = E = K ,
and divide by elemental volume, AAX, and simpjlfv
If = ' A
- 3.4 -
-------
at steady state QC/3t = 0)
A ¥x + ]
ox
Q/A = U
2
0 = F ~£ - U H (8)
as can be seen we have derived equation (2)„
Solution:
Equation (8) is a linear homogeneous second order differential
equation with constant coefficients whose general solution is
C - Ze(U/E)'X + J (9)
where
Z and J are unknown coefficients to be determined via boundary
conditions.
Biochemical Oxygen Demand (BOD)
BOD has historically been one of the major measures of the pollu-
tional strength of a water. It is usually defined as the amount of
oxygen required by organisms while stabilizing decomposable organic
matter under aerobic conditions. It is analyzed quantitatively by
placing a dose of the waste into a bottle whose contents simulate the
natural aerobic x;ater environment as closely as possible; that is,
water which is sufivated with oxygen, bar, the necessary nutrients,
has proper osmotic nnd temperature conditions, etc. The depletion in
oxygen is then rec./.ded over intervals of time and quantitative values
of BOD are calculated. Figure 4 shows the typical plot of these BOD
values for a water containing organic matter as well as some nitrogen.
- 15 -
-------
a
CD
H
03
o
0)
o
q
•H
Curve for TV Lai Demand
(Carbonaceous + Nitrogenous BOD)
Second Stage BOD
Curve for Carbonaceous Demand
Firs I Stage BOD
i i——i i i i i i
16 20 22 2U 26 26 30 12
Incubation Time In Davs
—r-
12
—r-
10
Figure A. BOD Curve for Typical Municipal Sewage.
As can be seen the reaction essentially takes place in two stages:
a carbonaceous stage (CBOD) during which about 70 to 80 percent of the
organic carbon is oxidized and a nitrogenous stage (NBOD) during which
biochemical oxidation of ammonia nitrogen and organic nitrogen is pre-
dominant. Since these reactions take place in a bottle and because
the actual reaction as it occurs in the stream is not fully understood,
there is some speculation as to the applicability of such test results
to actual water systems. At the present time, it is believed that highly
treated effluents will exert both first and second stage BOD simultaneously.
Untreated effluents will immediately exert a CBOD while the NBOD will
exert itself after a lag period.^ Another factor which is recognized
is that at low levels of DO (approximately <1.5 mg/1) nitrification prac-
tically ceases. The nitrogenous component can often be significant9
and the misunderstanding concerning its behavior is such that the user
must take care in modelling this parameter. The ES001 program is capable
of modelling NBOD as well as the conventional CBOD merely by inputting
the appropriate loadings and reaction rates for each component. The DO
deficit responses from each may be added to determine the total deficit
attributable to them. All theoretical development will be done in terms
of "BOD" which can signify either component.
- 16 -
-------
ES001 is also capable of analyzing either ultimate BOD or 5-day BOD.
A value of, F,
F =
ultimate BOD
5-day BOD
can be inputted for individual sections of the estuary when 5-day values
are being used. (When the ultimate values are used F is inputted as 1).
This option can be particularly useful in tin: verification of some models
as 5-day BOD's are much easier to obtain than ultimates. Of course, use
of this option depends on the types of waste involved and whether the
assumption of a uniform F for a section is realistic. ES001 will output
5-day BOD and the actual deficit in the strean when the option is exercised,
The equations for evaluation of BOD in an estuary can be developed by
considering a differential element of an estuar^ as shown in Figure 5.
I-- -si ti.ve ;•:
(2) // F~ Q*"
Figure 5. Differential Element of an Kstuary (BOD).
First, a simple case involving no waste loads will be examined to
illustrate the derivation of equation 3. Then a second more complex case
involving uniform loads will be derived. It is this second case which is
contained in ES001.
1) Case 1 - BOD with no waste loads
A mass balance is obtained, over the time interval At, since the law
of conservation of mass applies. The buildup in the element is evaluated
from a consideration of the material entering and leaving through advec-
tive and dispersive transport. The contribution of sources and sinks of
material within the element are also included in the inventory.
- 17 -
-------
Mass Balance: no waste load
entering element:
advection: Q L
3T
dispersion: - E1 A. 1
1 * 33T
leaving element:
advection: Q0 (L, + ~ AX)
/ J. OA
dispersion: - EZ ^ ^ (^ -f § AX)
sink:
BOD removal: VOL • L • RK
Equation:
HT _i = ^ ^ ^ ^
VOL = 0. L - E. A. 1 - n. (j + - AX)
' xrt-- .: I 9x
K2 A2 9X (L1 + 3X AX' - VOL-L"RK
assume
Q = Q1 = 02, E -
and divide by elemental volume, AAX, and siiiplity
3t A 3X 3X 3X
at steady state (9l,/3t = 0):
0 . _fi dL + E A _ Tj.RK
A dX ..,.2
- 18 -
-------
Q/A - u
ift
dx2
n uu
U dx"
- L-RK
(10)
a slightly more sophis-
Solutlon:
Equation (10) is a linear homogeneous second order differential
tion with constant coefficients whose general solution is
equa-
Sx
Ce
Vx
where
ditions
C0efflcients ^0 be determined via boundary
(11)
con-
S.V -
t
(Ha)
2) Case 2 - Now the effect of a uniform disch
of waste on the BOD will be examined.
WU(#/mile/day)
1 Ll +
Figure 6. BOD with Uniform Load.
- 19 -
-------
Mass Balance:
entering element:
Q1L1 - E1A1
leaving element:
Q2 (L, + AX) - E2 A2 j (LI + AX)
source: uniform load = WU-AX
sink: VOL-L-RK
Equation:
3L, „ ,T . 9L
3f vl "1 "1
VOL ~ = Q, L, - E, A ""1
_L r\ ••
3X
+ E0 A0 |rr (L. + |£ AX) - VOL-L-RK + WU-AX
2 2 dX 1 3X
assume
E - E = E , A = A = A , Q = Q = Q and
Ju — A, £• A. £
divide hy el.err.ontnl volume, AAX, set equation at steady state and simplify
0 dl, d L WIT
0 ~ V" T? + E —T ~ L'1-^K + •—-
dX ,2 A
dx
0/A = H
BL a .-, K l^ll .,. u ~ + L-RK (12)
A dx" dX
Solution:
Equation (12) is a linear nonhomogencr.us second order differential
equation with constant coefficients whoso general solution is the same
as in Case 1 /'equation (11))
L - Be SX + Ce VX
K
- 20 -
-------
The particular solution is
L = WU/(A-RK) (12a)
P
The total solution is therefore
T.TTT
(13)
H." IUN.
Dissolved Oxygen Deficit
Oxygen is soluble in water to a limited extent, being directly related
to the partial pressure of oxygen in the atmosphere. The saturation value
of oxygen is primarily a function of salinity and temperature. The dif-
ference between the saturation concentration and the actual concentrations
of oxygen is defined as the dissolved oxygen deficit.
Dissolved oxygen in natural waters can be employed to satisfy the
oxygen demand of pollutants, as measured by the BOD. In the development
of the equations for BOD, the rate of BOD removal was defined by the
parameter "RK". In a natural body of water, the rate of BOD removal does
not necessarily equal the rate at which oxygen is utilized to satisfy the
BOD. If organic matter, having a BOD, is removed by sedimentation, strip-
ping, or other physical means, oxygen is not consumed in this removal of
BOD. Since these phenomena can occur, it is necessary to employ a param-
eter in the dissolved oxygen material balance which indicates the rate
at which oxygen is utilized to satisfy BOD. This parameter is defined as
"DK". The parameter "DK" can be equal to, or less than the parameter "RK".
The rate at which oxygen is utilized is equal to the product of the
parameter "DK" and the ultimate BOD present. As the oxygen concentration
is reduced below the saturation level, the deficit is made up by atmos-
pheric oxygen going into solution. This phenomenon is termed atmospheric
reaeration. It occurs at a rate equal to the product of the dissolved
oxygen deficit and the reaeration coefficient, "AK".
The differential equation defining dissolved oxygen deficit may be
developed by constructing a material balance around a differential element
and including the various sources and sinks of dissolved oxygen within
the volume.
- 21 -
-------
Positive x
•
Flow - Q
Figure 7. Differential Element of an Estuary (DO deficit).
As in the development of the BOD equations, a simple case involving
no loadings will be derived first. It will be followed by the derivation
of the more sophisticated case which is used in program ES001. This
second case includes benthal demands and the effects of photosynthesis
and respiration of algae.
1) Case 1 - DO deficit with no waste load
The DO deficit is constructed in a similar fashion to the BOD. A bal-
ance is made between the amounts entering, leaving, and any sources or
sinks within the element.
Mass Balance:
entering element:
Q D - E A 1
11 11 7T~
leaving element:
_ ,~ . 3D
"^ AVN -p A _
TTTT AXJ - b_ A« -rrr
2 2
If
- 22 -
-------
source: deoxygenation due to BOD = DK-L -VOL
X
sink: reaeration = AK-D-VOL
Equation:
VOL
3 3D
+ E. A0 ^— (D. + -TTT AX) + DK-L -VOL - AK-D-VOL
f. 2. OA 1 OA X
E = \ = E2' A = Al = A2' Q = Ql = Q2 and
divide by elemental volume, AAX, set equation at steady state and simplify
2
DK-L = - E ^-4 + U ^- + AK-D ........... (14)
x , 2 dx
dx
The quantity "Lx" is the ultimate BOD concentration at a point X in the estuary
and is obtained from the BOD solution.
Solution:
Equation 14 is a linear nonhomogeneous second order differential
equation with constant coefficients whose general solution is:
D - GeSDX + He™ ................ (15)
g
where G and H are unknown coefficients to be determined via boundary
conditions.
SD, VD - [1± [1+]] ........... (15.)
- 23 -
-------
The particular solution is
DK-F
AK-RK
(BOD SOLUTION)
(16)
The BOD SOLUTION is equation (11), therefore, the complete solution is
D
Ge
SDX
+ Ke
™
(BOD SOLUTION) ... (17)
2) Case 2 - DO deficit with ben thai demand and algal effects
The complete dissolved oxygen deficit scheme of an estuary as pre-
sented in this model includes several sources and sinks: uniform waste
load, algae photosynthesis and respiration, and benthal demands.
Algae can influence the dissolved oxygen deficit profile. A daily cycle
of algae respiration can be approximated by a constant
respiration rate, :;RS;:, and an average daily photosynthesis rate of oxygen
production, "P" . The units of "RS" and "P" are the mass of oxygen used
or produced per day per unit volume.
Benthal demands result from sludge deposits on the streambed. A bottom
oxygen demand, BD, can be formulated, the units of which are mass of demand
per area of bottom per day.
Figure 8 will be referred to when deriving the complete DO deficit
solution used by this model.
Figure 8. Differential Element of Estuary for
DO deficit Approach.
- 24 -
-------
3X
Mass Balance:
entering element:
Q-i D., ~ E-, A.
11 13
leaving element:
3D 3 , 3D
Q (D + 1 AX) - E A ITT (D + %£ AX)
2 1 3Y" 2 2 3X 1 3X
sources: deoxygenation due to BOD = DK-L -VOL
X
respiration of algae: RS-VOL
benthal demand: BD-AX-WI
sinks: reaeration = AK-D-VOL
photosynthesis of algae = P-VOL
Equation:
M _ 3D 3D
3t 11 11 TTT; 2 1 9X
r
+ E_ A. -5^ (D. + |^ AX) + DK-L -VOL + RS-VOL
2. 2. OA 1 dX x
4- BD-AX'WI - AK-D-VOL - P-VOL
E = E± = E2, A = ^ = A2, 0 = Q1 = Q2 and
divide by elemental volume, AAX, set equation at steady state and simplify.
DK.L + (RS_P) + BD^AX^WI m _ E 9^D + u 3D + ^
x AAX *....£ oX
OA
note: ^4| - 4? (where HT is water depth)
A« AA Hi
- 25 -
-------
Therefore, the total equation is
2
DK-L + (RS-P) + = - E + U + AK-D .... (18)
Where
Solution:
Equation 18 differs from equation 14 in its nonhomogeneous part so
the general solutions must be the same:
D = GeSDX+HeVDX
g
The particular solutions for each source and sink are:
DK-F ^Y w
deoxygenation: - [Be +CeVA] + F-DK'WU
AK-RK-A
RS—P
respiration and photosynthesis: (20)
AK
benthal demand: • (21)
Letting RA = RS-P and calling it the "net algal effect", the complete
solution is therefore
_ SD-X _,_ „ VD-X . BD , -RA
D . Ge +He +_„___ + .„
, F.DK.M.J , DK sx . _ vx. „ /00.
[Be + Ce ] F ......... (22)
Table 1 presents a summary of the applicable equations for constant
cross-sectional area estuaries for several loading conditions and phen-
omena which normally influence water quality in an estuary.
- 26 -
-------
TABLE 1
SUMMARY OF EQUATIONS FOR STEADY - STATE SECTION ESTUARY ANALYSIS
Application
Upstream or
downstream
load
Upstream or
downstream
load with
bottom oxy-
gen demand and
algae
Uniform BOD
load WU (#/day/mi)
in section
Upstream or
downstream
load with uni-
form BOD
load, algae &
bottom deposit
in the section
BOD
„ SX „ VX
L = Be + Ce
T c SX . n VX
L = Be + Ce
L - EcSX 1 CeVX 1 W
L Ec ' Ce ' A-RK
L EaSX 1 CeVX 1 W
Dissolved Oxygen Deficit
_ „ SDX . „ VDX rDK-F , rR SX VX,
D = Ge -f He + [ .-•_ •„] [Be + Ce J
SDX n VDX BD ( RA ( ,DK-F , rEc,SX ( C£VX,
_ SDX TT VDX rDK-F , rn SX „ VX, DK-WU'F
"~ \J& '" ilfi T L*w^ -nV* [Jj© + LS j P ' £„ '
„ SDX . TT VDX . BD , RA . DK-WU-F . DK-F rr, SX „ VX,
HT-AK AK AK-RK'A AK-RK •*
-------
Other Applications:
Program ES001 was developed to model the Biochemical Oxygen Demand
(BOD), Dissolved Oxygen (DO) Deficit system which is presently the most
common measure of estuarine pollution. For this reason this documentation
focuses on this system. However, by analogy ES001 may also be employed
to model the distribution of other important substances. This is true
because the BOD-DO deficit system is representative of a general class of
reactions.
The BOD-DO deficit reaction might be generally described as a
"coupled system of two non-conservative substances which react with
first order kinetics." A schematic of this reaction is shown in
Figure 9:
Input :
Waste loads,
runoff loads,
etc.
Reaction in
Estuary System
Output:
BOD
response
of estuary.
Input :
BOD
response
of estuary
Reaction in
Estuary System
Output :
DO deficit
response
of estuary
Figure 9. BOD-DO Deficit System.
This oversimplified representation shows that in this particular
system, the basic input of waste load forces an output of BOD response
from the estuary system. This output in turn forms an input to the
system which forces a DO deficit response. Other systems might be re-
presented in analogous fashion as shown in Figure 10 if the assumptions
of the kinetics is -justifiable.
Input:
Polyphos-
phate load
Reaction in
Estuary System
*
Output :
Polyphos-
phate re-*
sponse of
estuary.
Input :
Polyphos-
phate re-
sponse of
estuary.
Reaction in
Estuary System
Output :
Orthophosphate
response of
estuary.
Figure 10. Polyphosphate-Orthophosphate System.
- 28 -
-------
Program ES001 is capable of modelling this type of reaction as well
as some subsets of it. A list of the options available using ES001 are
shown in Table 2.
Segmentation of Estuary Systems:
Equations have been developed in the preceding sections for estu-
aries having constant values of the parameters. Specifically, it was
assumed in each of the developments that the parameters AK, RK, DK, E,
A, and Q were constant. The solutions thus derived are valid for any
length of estuary in which the values of the individual parameters do
not change. At the location of changes in the value of individual parameters,
it is necessary to segment the system and apply the appropriate general solution
to each section.
Additional segments are required for any of the following condi-
tions which would cause a change in a parameter.
1. The magnitude of the dispersion phenomena may be substantially
reduced upstream of the point of maximum salt intrusion.
2. Fresh water flow may enter a system thus altering the advective
velocity and the dilution factor.
3. The cross-sectional area of the body of water may change, thus
altering the fresh water velocity and system volume.
4. Segmentation is required at the confluence or two or more bodies
of water.
For the specific BOD-DO deficit system:
5. Oxidation of nitrogenous compounds may not occur at the point
of waste introduction, or may require substantial oxidation
of carbon sources and the presence of high dissolved oxygen
concentrations. The system is segmented at the point of the
beginning of nitrogenous oxidation.
6. The rate of BOD removal RK may include the effects of settling,
stripping, etc., near the location of a raw discharge. In
this instance, the rate of 'BOD removal may be greater than the
rate of oxygen utilization, DK. Following the completion of
settling or stripping, the rate of BOD removal may be equal to
the rate of oxygen utilization. Segmentation is necessary
when the rate of BOD removal RK changes value.
- 29 -
-------
TABLE 2
System
Reaction Discussed in
This Documentation
Other Applications
Non-ccnservative substances
(coupled system, first order
reactive)
BOD-DO deficit
(p. 15-26)
Carbonaceous BOD-DO deficit*
Nitrogenous BOD-DO deficit*
Polyphosphate-orthophosphate,
etc.
o
I
Non-conservative substances
(single system, first order
reactive)
Conservative substances
(single system)
BOD (p. 15-21)
Conservative substances
(p. 13-14)
Coliform bacteria, etc.
Chlorides
Total Dissolved Solids
Pesticides (if their decay
rate is slow enough), «tc.
-ES001 has included the option to calculate each individual components of BOD and then add their DO
deficit to give a total deficit.
-------
7. Segmentation is required at the beginning and end of uniform
BOD loads, bottom oxygen demands or algal demands.
8. Toxic effects may begin or end resulting in changes in the
rate of BOD or nitrogenous oxidation. Segmentation is re-
quired when changes in rate occur.
9. The depth and/or velocity in an estuary may change, thus alter-
ing the value of AK, the reaeration coefficient.
10. Systems are segmented at the beginning and end of sections
where oxygen is produced by algae or rooted aquatic plants.
In addition, estuary systems must be sectioned at the location of
all waste inputs which are to be considered. There are two methods avail-
able to evaluate the influence of multiple waste inputs. It is possible
to construct a model which includes all waste inputs thus requiring system
segmentation at each input location.
A second method of evaluating multiple waste loads employs the su-
perposition technique which is valid for the linear equations describing
BOD and dissolved oxygen profiles. In the superposition method, the
effect of individual waste loads can be evaluated separately. The cumu-
lative effect of loads can be determined by adding the individual BOD
and dissolved oxygen deficit profiles.
Evaluation of Boundary Conditions:
General equations have been presented describing conservative sub-
stance, BOD and dissolved oxygen deficit profiles in estuary systems.
The appropriate equations are applicable to each section of the segmented
system. Since the differential equations which were developed in the
preceding sections are second degree, their solutions contain two arbi-
trary constants. The unknown coefficients Z and J in the conservative
equation, B and C in the BOD equation and G and H in the deficit equation
are evaluated by means of boundary conditions. Evaluations
of these conditions also provide a method of linking adjacent system
sections which were segmented in accordance with the criteria discussed
in the previous section. Since the solution in each segment contains two
unknown coefficients, the total number of boundary conditions that must
be evaluated in any system is therefore equal to twice the number of seg-
ments in the system.
- 31 -
-------
There are three classes of boundary conditions which are appropri-
ate:
1. Infinity conditions - At locations remote from a waste source,
the anticipated conservative substance, BOD and dissolved oxygen
deficit attributable to the waste source should approach or be
equal to zero. This can be stated mathematically by evaluating
the BOD or dissolved oxygen deficit equation, at positive or
negative infinity.
2. Concentration equality - At the junction of two segments, the
concentration calculated by the equation in either segment must
be the same. Thus, the profile equation in the upstream seg-
ment must equal the profile equation in the downstream segment
at the location of the junction.
3. Mass Balance - It is possible to develop a material balance
around an element at the junction of two sections. The quan-
tity of substance entering the element must equal the quantity
of BOD leaving the element.
Proper application of the three general classes of boundary condi-
tions will generate sufficient equations to evaluate all unknown coef-
ficients. In the following pages, examples of evaluation of boundary
conditions will be presented for estuary systems.
X = 0
W
In order to illustrate the specific application of boundary condi-
tions, a simplified estuary will be considered with a single waste
source of BOD located at x = 0. The estuary will be considered infin-
itely long in both directions and having constant cross-sectional area,
flow and parameters AK, RK, DK and E. The estuary must be divided into
two segments at the point of waste input. Section one extends upstream
from the point of waste input and section two extends dox-mstream from
the point of waste discharge. The BOD equations applicable to each
section are:
BOD: Section 1 (Upstream)
S X V X
LX = Be + C^e . (23)
- 32 -
-------
BOD; Section 2 (Downstream)
V X
C2e
Note that each unknown coefficient and roots S and V are subscripted
by section. In like manner the area, flow, reaction and dispersion are
subscripted by section. The sectioning procedure results in four (4)
unknown coefficients, thus requiring four (4) boundary conditions. The
boundary conditions applicable to the example system are:
Boundary Condition Location and Description
//I @ X = - °° L = 0
#2 @ X = + °° L = 0
#3 @ X = 0 LX = L2
#4 @ X = 0 mass balance
Boundary Condition //I
At large distances upstream from the waste source, the BOD present
in the estuary due to the waste discharges becomes small. Therefore,
as X approaches negative infinity (- °°), the BOD in section one
approaches zero.
L = 0 at X = -oo
S X V X
L = Be + C e
S., = positive number
V- = negative number
-f- c (- oo)
1 1
X -*• ^^^^^^ f\
- 33 -
-------
large number
If L = 0 at X = - «> C must =• 0
Equation (23) for section one reduces to:
S X
Ll " V
Boundary Condition #2
Using the same reasoning:
L = 0 at X = +
V X
= e » large number
e
0
Equation (24) for section two reduces to:
v v
o/L
L2 = V
Boundary Condition #3
At X = 0, the BOD concentration calculated by the equation for the
BOD profile in section one must equal the concentration as calculated
using the BOD profile equation in section two.
1^ = L9 at X = 0
S X V X
V ' V
- 34 -
-------
@ X = 0;
S1X
= e
V2X
and
(25)
Boundary Condition #4
Considering a cross-sectional plane at X = 0, a flux or material
balance can be obtained.
(- X)
W = point source
Material entering the section:
q.L- - E A
+ W
Material leaving the section:
Q?L. - EA2 dL2
22 2 2 dX™
Balance:
dL,
d.X
W =
dX
or
Q.L_
•
dL
dL
Q.L. + E A. l - E A9 2 - W
11 1:L~ 22
0 . . (26)
dX
dX
Since
Q = Q and
at X = 0
- 35 -
-------
L - EA 2 " W (27)
dX Z dX
Evaluating
^1 and ^2 at X = 0:
dX dX
dx
IT S1X
dLl = BS.e1
dx~ -l 1
X = 0; 1 = B S
dX ~r~
dX
dX
@ x = 0; 2 - C2V?
Substituting into Equation (27)
E3A1B1S1 ~ E2A2C2V2 = W
Grouping terras :
C2 ("E2A2V2) + Bl (E1A1S1) = W •••••...... (28)
Summary of Boundary Equations
B - C = 0 (concentration equality) ....... (29)
J_
-------
assuming
E = E = E and
= A
There are two boundary condition equations and two unknown coeffi-
cients remaining, thus it is possible to solve for the unknown coeffi-
cients.
Combining equations:
BI (SXEA) -
= w
(30A)
w
EA
and
C2 = W
L
l EA
W
W
slx
(31)
EA
w
(32)
The techniques for evaluating boundary conditions illustrated in the
previous example can also be employed for evaluating the upknown coeffi-
cients in the dissolved oxygen deficit or conser-'ative substance.
Additional Boundary Conditions:
There are two (2) special case boundary conditions wMch will be
illustrated. These are for a system bounded by a dam, and an estuary
entering a body of water with a fixed concentration.
Mass Balance at a Dam
->•£,
I A__
Dam
- 37 -
-------
Material entering the element:
QD • LD
Material leaving the element:
dL,
Ql • Ll - E1A1
Balance: (in = out)
QD • LD = Q L
dX
(33)
dX
The boundary condition can then be fully evaluated by substituting
the appropriate equations for L and dL^ into Equation (33).
1 air
Completely Mixed Bay with Exchange
Concentration Equality:
B
Mass Balance:
In this instance, the mass balance is obtained around the entire
bay.
¥
Tidal (R)
Exchange
- 38 -
-------
Note: VOL = Volume
Entering Bay:
Leaving Bay:
dL,
Q L - E A. "1 + W
-L J. ll. "T"™"
O"
OJ-,., + R-L -VOL
'.DO is
VOL
Balance:
dT,
dx"
'1 4- W =
c Q_ + VOL (R •!- RK) L
B B T
Evaluating L]_ and -T— - as functions of distance and the unknown coef-
ficients results in the required boundary equation.
An additional concept and phenomena has been introduced in the con-
sideration of the completely mixed bay. The tidal exchange "R" provides
a method of evaluating the transport of material out. of the bay into the
ocean by the tidal volume.
Tidal exchange "R" is defined bv:
R =
tidal volume
mean vr-lu^e of the bay
nuciber_pf tidal cycles
day
In the simplified estuary problem used as a deta.5J.ed e:-:r.mpln of
boundary condition evaluation, it was possible to solve f >:r the unknown
coefficients B .°nd c, manually. As the estuary beroir/ i; "or?, complex, it
is impractical to evaluate the coefficients manually. "r; This latter
case, it is convenient to arrange the unknown coefficients end boundary
condition equations into a matrix a:?.d employ compT.itai:ic:!al techniques
for numerical p/*alintion of fche unknown coefficient?, ?IT> boundary
condition equation (25) and (30A) from our simplified ecr'ttrvle can be
arranged into irr.tri:: form:
ff
I
EAS,
-1
-EAV,
Matrix
Vector of
Unknown
Coefficients
Vector of
- 39 -
-------
Another way of presenting this matrix is as shown below. Both
ways are used in this report.
B.C. No. Bl C2 Forcing Function
3 1-1 =0
4 EASjL -EAV2 = W
The basic method of establishing a matrix of unknown coefficients
can be employed regardless of the complexity of the system.
There are several general considerations which are useful in evalu-
ating the unknown coefficients and establishing matrices.
1. Evaluation of boundary conditions always results in a square
matrix which means that unique solutions for the values of the
unknown coefficients are available.
2. At the junction of "n" segments, there are always generated
one mass balance equation and (n - 1) concentration equality
equations.
3. The boundary conditions for the BOD and dissolved oxygen deficit
solutions are identical with two exceptions. First, the roots
of the BOD solution, S and V, are replaced by the roots of the
deficit solution, SD and VD when DO deficit calculations are
made. Second, in the case of boundary conditions in which flow
enters or leaves a completely mixed bay, a reaction rate for the
bay will appear as a part of a matrix element. This will change
from the BOD removal rate (RK) when solving for BOD to the
reaeration rate (AK) when solving for the deficit.
4. The right hand matrix of forcing functions are different for
the BOD and dissolved oxygen deficit solution.
- 40 -
-------
SAMPLE PROBLEM
WU (uniform load)
,
(3)
\ 4 1 V 1 1 1
(2)
i
\
(point
(1)
C
tf
load)
X = Xo
X = 0
Figure 11
The purpose of this problem is to apply the theory of the previous
sections to the relatively simple estuarine system shown in Figure 11.
This estuary is divided into three segments. Section 1 stretches
from minus infinity to a point which will be designated as X = 0 at
which a point load enters the system. Section 2 extends a distance, Xo,
downstream from the zero point. A uniform waste load is input to
this reach along its entire length. Finally, Section 3 extends from Xo
to positive infinity.
The following equations can then be written for each section for
BOD:
Section 1:
S X V X
= B e + C^e
(34)
Section 2:
L2 ' V
S2X
V2X
WU
(35)
Section 3:
L3 = V
S X
V X
(36)
- 41 -
-------
Six boundary conditions can then be applied:
1) L = 0 @ x = - °°
2) L = 0 @ X = + °°
Concentration equality equations:
3) LX = L2 @ X = 0
4) L2 - L3 @ X = Xo
Mass balance equations:
5) mass balance @ X = 0
6) mass balance @ X = Xo
These six conditions may now be applied to equations 34, 35 and
36 and will provide a means of evaluating the unknown coefficients.
As shown on page 33, the infinity boundary conditions may be applied
with the result that
@ X = - °° C = 0
(3 X = + °° B = 0
Since C- = 0 the following equations result at X = 0
Ll = Bl
2 2
Boundary Condition 3 may then be applied:
Ll = L2
A2-RK2
- 42 -
-------
Therefore ,
Bl - B2 - C2 = A .............. (37)
Since 63 = 0 the following equations result at X = Xo
S2X° V2X° WU
V Xo
Applying Boundary Condition 4
L2 ' L3
S2Xo V Xo V Xo
V + V * + Ax = V 3
Therefore,
S Xo V Xo V Xo
V + V - V
As can be seen two coefficients have been eliminated (C^ and 63) and
two equations have been generated. An additional two equations are re-
quired to obtain four equations for the four remaining unknowns (B^, 62,
G£» 63). The mass balance conditions provide these equations.
Boundary condition 5 @ X = 0
Q L - Q L + E A dLl - E A dL2 = W ....... (39)
22 11 i;Ld5r 22dX~
Again recall that
S X
Lx = I^e 1 ...................... (40)
S X V X
(41)
- 43 -
-------
Substituting 40 and 41 into 39
Bl(ElAlSrQl) + B2(Q2~E2A2S2) + C2(Q2~E2A2V2} = W ..... (42)
Similarly, boundary condition number 6 can be evaluated
- V2 + E2A2 * + ^3L3 - E3A3 l " °
dx dx
The resulting expression is
S Xo V Xo
B [e Z (EAS-Q)] + C [e Z
V Xo
+ C3 [e (Q3-E3A3V3)] = 0 ................. (43)
Equations 37, 38, 42 and 43 now contain 4 unknowns and can be written
in matrix form as in Table 3. All matrix elements can be obtained from
the geometry of the system. The forcing function is known and the values
of B , B , C and C can be obtained. By substituting them into equation
34, 35 and 36, the concentration of L can be computed for any point in
the estuary.
- 44 -
-------
TABLE 3
-1
-1
(ElAlVQl)
(Q2~E2A2S2}
(Q2~E2A2V2)
-e
S Xo
V Xo
V Xo
w:
W
wu
-------
The preceding sections have presented techniques for developing
equations which define the concentration profile of BOQ dissolved oxygen
deficit, and conservative substances. In addition a discussion of segment-
ation and methods of evaluation of boundary conditions and unknown coeffi-
cients for rather simple cases have been presented. As was illustrated
in the sample problem of the previous section, the mathematical formula-
tion of ES001 reduces to the following:
[A] • (X) = (b) (44)
where
[A] - is the matrix of known factors
(b) - is the forcing function vector
(X) - is the vector of unknown coefficients
This relationship is formed by applying appropriate boundary condi-
tions to the solutions to the water quality equations of the previous
sections. It is solved by inverting [A] which allows the calculation
of (X)
(X) = (b) • [A]"1 (45)
In the following sections the method of analysis to formulate equa-
tion 44 will be explained. Attention will be devoted to the actual
formulation of [A] and (b) via the specific boundary conditions of
ES001.
46 -
-------
Method for Evaluation and Presentation of
Boundary Conditions for Complex Estuaries:
The matrix of unknown coefficients generated for complex water quality
models can become large. Presentation of the specific values of the ele-
ments of the matrix for BOD and deficit solutions would be obscure. In
view of the above, and the fact that complex systems must be evaluated on
a computer, the following procedure can be adopted for presentation of the
water quality models.
- The BOD and dissolved oxygen deficit equations can be
presented in general form.
- A table can be developed defining the specific functional
forms of the general equations for the BOD and dissolved
oxygen deficit solutions.
- The concentration boundary condition may be generalized
in four (4) classes which cover all of the concentration
equality equations encountered in the present system.
The four classes of concentration equality boundary
equations can be presented in a form which is amenable
to matrix solution.
- The first derivative of the generalized BOD and dissolved
oxygen deficit concentration profile equations can be
defined in tabular form.
- Ten (10) generalized solutions to the mass balance boundary
condition equations can be presented. These equations
cover all of the mass balance conditions required in the
present systems. The mass balance equations have been
presented in a form amenable to matrix solution of the
simultaneous equations developed.
- The boundary conditions may be classified and tabulated.
- The specific boundary conditions from the general classi-
fication listing can be presented in tabular form for
each junction. The appropriate designation of upstream
and downstream must also be indicated.
- A model can then be constructed by joining segments.
- 47 -
-------
1) Presentation of equations in general form.
The formulas describing water quality in any specific segment (I)
can be written in the general form:
L(I) = B(I)f (X,I) 4- C(I)? (X,I) + YB(I) .... (46)
D(I) = G(I)f (X,I) + H(I)7 (X,I) + YD(X,I) . . . (47)
3. 3.
where :
(I) = Section number
(X,I) = Indicates that the expressions are a function
of both section and location "X" in the section.
L(I) = BOD in Section I.
D(I) = Dissolved oxygen deficit in Section I.
B(D
and = Unknown coefficients for the BOD equation.
and = Unknown coefficients for the deficit equations.
f (X, I) = Functional forms for the BOD equation. The
_ and specific forms are defined in Table 4.
fr(X,I)
YB(I) = Particular integral for the BOD solution as
defined in Table 4.
YD(X,I) - Particular integral for the deficit solution
as defined in Table 4.
fr(X,I), 7r(X,I), fa(X,I), and ?a(X,I) are functions of distances X and
segment number "I". YB(I) is a function of the segment number "I" alone,
and YD(X,I) is a function of segment and location.
- 48 -
-------
There are three (3) possible special cases which do not have the form
indicated above. They are:
1. Negative infinity boundary conditions where the only unknown
coefficient is the B or G coefficient.
L(I) = B(I)f (X,I) + YB(I) (48)
= G(I)f (X,I) + YD(I) (48a)
a
2. Positive infinity boundary conditions where the only unknown
coefficient is the C or H coefficient.
C(I)fr(X,I) + YB(I) (49)
H(I)T (X,I) + YD(I) (49a)
a
3. Completely mixed volumes.
L(I) = B(I) (50)
(50a)
For completely mixed volumes the BOD and deficit does not vary with
distance and can therefore be defined by constants B Ii and G {'I; respec-
tively.
- 49 -
-------
2) Specific functional forms of the general equations for the BOD and
D.O. deficit solution.
Function
fr(X,I)
Value
fa(x,D
fa(x,D
YB(I)
YD(X,I)
ev(i)x
SD(I)X
e
eVD(I)X
WU(I)
ED(I)
HT(I).AK(I) AK(I)
AK(I)-RK(I)
TABLE 4
-50-
-------
3) Concentration Equality Equations
flow
AO
Figure 12: One Dimensional Estuary
To illustrate the formulation of the concentration boundary condi-
tions in computational form the following example has been included.
Figure 12 is a one dimensional estuary with a change in cross-
sectional area at Xo. According to the concentration equality boundary
condition for two double sections
Ll = L2
@ X = X
fr(XQfl)
B(2) fr(XQ,2)
C(2)
Group terms:
fr(X,l)
YB(1) =
YB(2)
fr(X,l) - B(2) fr(X,2)
- C(2) f (X,2) = YB(2) - YB(1)
Consulting Table 4 we can then formulate the following matrix
denoting it as Boundary Condition "J"
BC#
J
B(l)
S(l)Xo
e
C(l)
V(l)Xo
e
B(2)
S(2)Xo
C(2)
V(2)Xo
-e
- 51 -
-------
with the forcing function for BOD
BFORC(J) = YB(2) - YB(1)
for DO deficit the matrix would be the same but with fa and fa replac-
ing fr and fr and
DFORC(J)
YD(X,2) - YD(X,1)
Boundary condition "J" can also be expressed in the following way
substituting the functional values shown below
fr(x,i)
fr(x,i)
fr(X,2)
fr(X,2)
B(2)
C(2)
Now, using the above notation, the boundary conditions of ES001
can be generalized. Thev will be set up for BOD with additional notes to
illustrate the procedure for DO deficit.
CLASS 1 - Boundary Condition - (Concentration Equality
Between Two Double Sections5^;
K
flow
Elements of matrix for Sections I and K and Boundary Condition J
BCD
c(D
B(K)
C(K)
f (X,I)
r
? (X,I)
r
-f (X,K)
r
-?r(X,K)
* Note - a double section is defined as a section which is bounded on both
ends by another section.
- 52 -
-------
BFORC(J) = YB(K) - YB(I)
DFORC(J) = YD(X,K) - YD(X,I)
For deficit Matrix, replace fr and fr by fa and fa, respectively.
CLASS 2 - Boundary Condition (Concentration Equality, Negative
Infinity Boundary for Section I, with C(I) = 0)
flow
Same as Class 1 Boundary Condition, but with
Forcing functions are the same as for Class 1 Boundary Condition.
CLASS 3 - Boundary Condition (Concentration Equality, Positive
Infinity Boundary for Section K, with B(K) = 0)
flow
(CO)
K
I
^fnii*
Same as Class 1 Boundary Condition, with
B(K)
Forcing functions are the same as for Class 1 Boundary Condition.
CLASS 4 - Boundary Condition (Concentration Equality,
between a double section and a completely
mixed volume)
flow
- 53 -
-------
BCD
fr(x,i)
c(D
B(KK)
= fr(X,I)
-1
BFORC(J) = YB(I)
DFORC(J) = YD(X,I)
For deficit matrix replace, fr and fr with fa and fa.
4) Specific functional forms of the derivatives according to x of the
general equation of the BOD and DO deficit solution.
Function
f-" /v" T \
1 A. • J. )
r
f"r(x,i)
t^V
Fa(x,i)
YB(I)
YD(X,I)
Value
S(I)eS(I)X
V(I)eV(I)X
SD(I)eSD(I)X |
1
VD(I)eVD(I)X |
Zero
DK(I} • Ffl') ' — "
j-^ivy-1-,/ L \J- J rp^T^-F ^ V T^ 1 P^T^-f- ^"V T *k 1
AK(I)-RK(I) r ' » r ' »
TABLE 5
54 -
-------
5) Mass Balance Equation
Additional classes of boundary conditions are required for the mass
balances. Before evaluation of mass balance boundary conditions, it is
necessary to define the first derivatives of the concentration profile
equations;
dL(l)
dx
B(I)fr(X,I)
C(I)fr(X,I)
YB(I)
= G(I)f
H(I)f
YD(X,I)
Table 5 presents the definitions of the functions f , f , f , f ,
•* * rraa
YD, and YD used in the present models.
Generalized solutions are available for the mass balances when sections
I, K, and KK are considered and the boundary equation number is defined
by J. It must be noted that Section I is defined as upstream, or more
negative than Section K. (That is, flow from Section I enters Section K.)
The direction of Section KK must be indicated.
CLASS 5 - Boundary Condition Class Balance,
Two Adjacent Sections)
1
K |
i
1 i
u...
I
' I
J
4
'.J
w
NOTE: The portion enclosed by the dotted lines in the figure denote
the element about which the mass balance is obtained as on
page 35.
- 55 -
-------
B(K)
C(K)
(X,I) - Q(I)fr(X,I)
- Q(I)fr(X,I)
Q(K)fr(X,K) - E(K)-A(K)-f^(X,K)
Q(K)fr(X,K) - E(K)-A(K)-fr(X,K)
BFORC(J) = W(n) (load entering at junction) +
Q(I)YB(I) - Q(K)YB(K)
DFORC(J) = Q(I)YD(X,I) - E(I)•A(I)-YD(X,I) +
*
E(K)-A(K)-YD(X,K) - Q(K)YD(X,K)
CLA.>S 6 - Boundary Condition (Mass Balance with Negative
Infinity Boundary for Section I makes C(I) = 0)
+»
V
1
K t
1
I t
i
.1 (-00)
1
t 1
Same as for Class 5 Boundary Condition with
CLASS 7 - Boundary Condition (Mass Balance with Positive
Infinity Boundary for Section K makes B(K) = 0)
(+
K
1 L J
W
- 56 -
-------
Same as for Class 5 Boundary Condition with
B(K)
CLASS 8 - Boundary Condition (Ilass Balance with Flow
Entering Sections K and KK from Section I)
C(I), B(I), C(K), and B(K) are same as for Class 5 Boundary Condition.
B(KK)
C(KK)
Q(KK)f (X,KK) - E(KK)-A(KK)-f'(X,KK)
Q(KK)f~ (X,KK) - E(KK)-A(KK)-7 (X,KK)
BFORC(J) = W(N) + Q(I)YB(I) - Q(K)YB(K) - Q(KK)YB(KK)
DFORC(J) = DFORC(J) in Class 5 Boundary, plus:
E(KK)-A(KK)-YD(X,KK) - Q(KK)YD(X,KK)
CLASS 9 - Boundary Condition (Mass Balance with Flow
Entering Section K from Sections I and KK)
- 57 -
-------
B(I), C(I), B(K) and C(K) are the same as for Boundary Condition 5.
B(KK) E(KK)-A(KK).f (X,KK) - Q(KK)-f (X,KK)
*
C(KK) E(KK)-A(KK)-I (X,KK) - Q(KK)-T (X,KK)
BFORC(J) = W(N) + Q(I)YB(I) + Q(KK)YB(KK) - Q(K)YB(K)
DFORC(J) = DFORC(J) as in Class 5 boundary condition,
plus:
Q(KK)-YD(X,KK) - E(KK) • A(KK) -YD(X,KK)
CLASS 10 - Boundary Condition (Mass Balance at Section
with Advective Flux, only) - DAM
B(K)
C(K)
K ||
1
1
1
DAM .
Q(K)-fr(X,K) - E(K)-A(K)-f~(X,K)
Q(K)-fr(X,K) - E(K)-A(K).fr(X,K)
BFORC(J) = Q(DAM)-L(DAM) - Q(K)YB(K)
DFORC(J) = E(K)-A(K)-YD(X,K) - Q(K)-YD(X,K) +
DEF(DAM)-Q(DAM)
L(DAM) and DEF(DAM) equal the concentration of BOD and dissolved
oxygen deficit in the flow over a dam while
Q(DAM) is the flow over the dam.
- 58 -
-------
B(I), C(I), B(K) and C(K) are the same as for Boundary Condition 5.
B(KK)
C(KK)
E(KK)-A(KK)-f'(X,KK) - Q(KK)-f (X,KK)
E(KK)-A(KK)-f (X,KK) - Q(KK)-f (X,KK)
BFORC(J)
W(N)
Q(KK)YB(KK) - Q(K)YB(K)
DFORC(J) = DFORC(J) as in Class 5 boundary condition,
plus: *
Q(KK)-YD(X,KK) - E(KK)•A(KK)•YD(X,KK)
CLASS 10 - Boundary Condition (Mass Balance at Section
with Advective Flux, only) - DAM
K ||
1 ^
1
l.DfJ
B(K)
C(K)
Q(K)-fr(X,K) - E(K)-A(K)-r(X,K)
Q(K)-fr(X,K) - E(K)-A(K).fr(X,K)
BFORC(J) = Q(DAM)-L(DAM) - Q(K)YB(K)
DFORC(J) = E(K)-A(K)-YD(X,K) - Q(K)•YD(X,K) +
DEF(DAM)-Q(DAM)
L(DAM) and DEF(DAM) equal the concentration of BOD and dissolved
oxygen deficit in the flow over a dam while
Q(DAM) is the flow over the dam.
- 58 -
-------
CLASS 11 - Boundary Condition (mass Balance with Section I
Leaving Completely Mixed Volume, KK)
B (KK) = VOL (KK) • RK (KK)
BFORC(J) = W (bay load) - Q(I)YB(I)
DFORC(J) = E(I)-A(I)-YD(X,I) - Q(I)-YD(X,I) +
F(KK)•DK(KK)B(KK)•VOL(KK)
For DO deficit change RK(KK) of B(KK) to AK(KK). Note that a
load should not be thought of as entering at the junction in this
case. This is so for all the bays since the balance is taken about
the entire bay and all loads are assumed to flow directly into the bay.
CLASS 12 - Boundary Condition (Mass Balance with Sections I and K
Leaving Completely Mixed Volume, KK)
- 59 -
-------
For two sections leaving a completely mixed volume, add B(K) and
C(K) $rtiich are identical in form to B(I) and C(I) of the Class 11
Boundary Condition)
To the forcing functions, add:
BFORC(J) plus - Q(K)-YB(K)
**
DFORC(J) plus E(K)-A(K)-YD(X,K) - Q(K)YD(X,K)
CLASS 13 - Boundary Condition (Mass Balance, Section I Entering
a Completely Mixed Volume, KK, with a Tidal Exchange
and Flow from the Volume to an Infinite Medium)
B(KK)
BFORC(J)
DFORC(J)
- Q(I)fr(X,I)
- Q(I)fr(X,I)
(Q(KK) + R(KK)VOL(KK) + VOL(KK)•RK(KK))
W(KK)
Q(I)YD(X,I) -
F(KK)•DK(KK)B(KK)•VOL(KK)
For DO deficit change R(KK) of AM(KK,J) to AK(KK).
- 60 -
-------
CLASS 14 - Boundary Condition (Mass Balance, Sections I and K Entering
a Completely Mixed Volume, KK, with a Tidal Exchange, R,
and Flow from the Volume to the Infinite Medium)
R
Q(KK)
.!/
,
!\
For two sections I and K, entering a completely mixed volume,
add B(K) and C(K), which are identical in form to B(I) and C(I) to the
Class 13 Boundary Condition.
To the Forcing Functions, add:
BFORC(J) plus + Q(K)YB(K)
DFORC(J) plus Q(K)YD(X,K) - E(K)A(K)YD(X,K)
6) Tabulation of Boundary Condition
The boundary condition classifications applicable to each combination
of sections are described in Table 6.
Using Table 6 the specific boundary condition can be tabulated for
each junction in the system. A model can then be constructed by joining
segments.
- 61 -
-------
TABLE 6
BOUNDARY CLASSIFICATIONS
Boundary
Classification Types of Applications
Concentration equality between two double sections
(when both sections have two (2) unknown coefficients
and are bounded at both ends.)
Concentration equality between a double and a negative
infinity section (when one section has two (2) unknown
coefficients and the second section has only B or G
coefficient.)
Concentration equality between a double and a positive
infinity section (when one section has two (2) unknown
coefficients and the second section has only C or H
coefficient.)
Concentration equality between a double section and a
completely mixed volume.
Mass balance between two double sections, both sections
having two (2) unknown coefficients.
Mass balance between a double section and a negative
infinity section when one section has two (2) unknown
coefficients and the second section has only B or
G coefficients.)
Mass balance between, a double section and a positive
infinity section (when one section has two (unknown)
coefficients and the second section has only a C or H
coefficient.
Mass balance for three (3) double sections with flow from
one section dividing into the two sections.
Mass balance for three (3) double sections with flow from
two sections combining into the third section.
- 62 -
-------
TABLE 6 (Continued)
BOUNDARY CLASSIFICATIONS
Boundary
Classification Types of Applications
10 Mass balance at a section with, only advective flux
entering (dam).
11 Mass balance between a double section and a
completely mixed volume. Flow leaves the bay
through the section.
12 Mass balance between two (2) double sections and
a completely mixed volume. Flow leaves the bay
through the sections.
13 Mass balance between a completely mixed volume
and a double section when flow enters bay through
the section.
14 Mass balance between a completely mixed volume
and two (2) double sections when flow enters bay through
two sections.
- 63 -
-------
NOMENCLATURE
Symbol Definition
A Cross-sectional area
AK Reaeration coefficient
B BOD constant (B associated with S)
BD Bottom oxygen demand
BOD Biochemical Oxygen Demand
c Concentration of some mass
C BOD constant (C associated with V)
CBOD Carbonaceous Component of BOD
C Saturation concentration of dissolved oxygen in
water at a particular temperature and solids
concentration
D Concentration of DO deficit
D General solution to the DO deficit concentration
s equation
DK Deoxygenation coefficient
DO Dissolved Oxygen
E Dispersion coefficient
F Ultimate to 5-day BOD ratio
G DO deficit constant (G associated with S)
H DO deficit constant (H associated with V)
HT Water depth
J Conservative substance constant
L Concentration of BOD
L, Concentration of BOD in bay
L. Concentration of BOD in flow over a dam
d
L General solution to the BOD concentration differ-
^ ential equation
Units
[1/T]
[M/L3]
[M/L2/T]
[M/L3]
[M/L3]
[M/L3]
IM/L3]
[M/LJ]
[M/L3]
[M/L3]
[1/T]
IM/L3]
[I/VT]
[M/L3]
[M/L3]
[L]
[M/L3]
[M/L3]
[M/L3]
[M/L3]
[M/L3]
- 64 -
-------
Symbol Definition Units
o
L Particular solution to the BOD concentration [M/L ]
differential equation (same as YB)
o
L Ultimate BOD of a water [M/L ]
N1BOD Nitrogenous component of BOD [M/L3]
Q Net advective flow of an estuary [L /T]
Q Net advective flow out of a bay [L3/T]
Q Flow over a dam [L /T]
P Algal photosynthesis effect [M/L3/T]
R Tidal exchange coefficient [1/T]
RA Net algal photosynthesis and respiration effect [M/L /T]
RK Removal coefficient of BOD [1/T]
RS Algal respiration effect [M/L /T]
2
±ZS Sum of sources and sinks [M/L /T]
S Argument of exponential in BOD solution associated [1/L]
with B
2- [i + U + ^]1/2]
^ IT
SD Argument of exponential in DO deficit solution [1/L]
associated with G
?=• [1 + [1 +
U
t Time [T]
U Net advective velocity [L/T]
V Argument of exponential in BOD solution [1/L]
associated with C
zt,
- 65 -
-------
Symbol
VD
VOL
W
WI
WU
b
YB
Y
c
YC
Yd
YD
Definition
Argument of exponential in DO deficit solution
associated with H
4AK-E
u2
1/2
] 1
Volume
Point waste load
Width
Uniform waste load
Distance along an estuary in direction of flow
Sources and sinks of BOD
Particular integral of BOD concentration differ-
ential equation
Sources and sinks of conservative substances
Particular integral of conservative concentration
differential equation
Sources and sinks of DO deficit
Particular integral of DO deficit differential
equation
Conservative substance constant
Concentration of conservative substance
Units
[1/L]
[M/T]
[L]
[M/L/T]
[L]
[M/L3/T]
[Tt/L3]
[M/L3/T]
[M/L3]
[M/L3/T]
[M/L3]
[M/L3]
[M/L3]
- 66 -
-------
SENSITIVITY
- 67 -
-------
System sensitivity may be defined as the relative change in water
quality associated with changes in the parameters or the geometry of
the system. There are three (3) significant areas where system sensi-
tivity has a bearing on the use and application of mathematical models
for water quality.
In the first instance, the relative sensitivity of the system to
changes in the value of a parameter can be employed to determine the
accuracy with which this parameter should be measured. The accuracy
with which a parameter should be measured will in turn dictate the
amount of money and effort to be spent in conducting field tests.
As an example, a system may be relatively insensitive to fresh water
flow and dispersion. Efforts to develop refined estimates of these
parameters can be minimized. If, on the other hand, the system is
sensitive to the nature and magnitude of waste inputs and to the value
of the reaeration coefficient, expenditures to develop realistic and
well-defined estimates of these parameters and inputs are justified.
In summary, information on the value of all system parameters, inputs
and geometry are required. Expenditures for developing highly refined
estimates of the various items should be viewed in terms of their in-
fluence on the water quality profiles; that is, the systems' sensitivity.
A second area where knowledge of the system sensitivity is impor-
tant is in developing a realistic set of consistent parameter values.
It is necessary to obtain calculated water quality profiles which are
in agreement with observed data in order to validate the model. Knowl-
edge of the influence of changes in parameters on calculated profiles
will assist in establishing the direction and magnitude of changes in
the assigned value of parameters and inputs. This will speed up the
process .of determining a consistent set of parameters which will pro-
vide agreement between calculated and observed profiles at various
temperature and flow conditions.
The third area where knowledge of system sensitivity is significant
is in selection and evaluation of the methods for improving water quality.
As an example, if water quality is insensitive to changes in fresh water
flow, this would indicate that low flow augmentation may not be a real-
istic method of improving water quality. On the other hand, the system
may be sensitive to tidal exchange and changes in system volume. This
indicates that consideration should be given to methods of increasing
tidal exchange and increasing the system volume as a means of improving
water quality.
It should be noted that the sensitivity of the system to a change
in the value of a single parameter will vary depending on the specific
value of all parameters; several general indications of system sensitiv-
ity can be defined. These are:
- 68 -
-------
1. The linear nature of the mathematical models results in
changes in the calculated BOD and dissolved oxygen deficit
profiles in direct proportion to changes in waste inputs.
That is, a 10 percent reduction in a waste discharge load
will result in a 10 percent reduction in the calculated
BOD and dissolved oxygen deficit profiles associated with
the load.
2. Steep concentration gradients will tend to be destroyed by
relatively small increases in the value of the dispersion
coefficient. It should be noted that the influence of
changes in the dispersion coefficient "E" are reduced when
relatively flat profiles (low concentration gradients) are
present.
3. The net influence on water quality resulting from changes
in the value of specific parameters will vary depending
upon the location of the point in the system being con-
sidered .
To illustrate this point, Figures 13 to 17 are presented. The Fig-
ures contain plots of a plane formed by the variable fresh water velocity
"U", the dispersion coefficient "E", and the calculated dissolved oxygen
deficit, "D". The values of RK and AK are constant and each Figure illus-
trates the variation of the calculated dissolved oxygen deficit at a
specific location for various combinations of values of "U" and "E". In
Figures 13 to 15, for locations upstream of the point of the waste dis-
charge, one notes that the maximum deficit always occurs at low values
of "U" and "E". The shape and magnitude of the changes in curvature of
the planes describe the influence of changing the parameters "U" and "E".
Figures 16 to 17 are the comparable planes for locations downstream from
the point of waste introduction. It will be noted from these Figures
that for a constant value of the dispersion coefficient, "E", the fresh
water velocity "U" at which the maximum deficit occurs increases as the
location studied becomes more remote from the waste input location.
That is, at a station located 10 miles from the point of waste introduc-
tion, the dissolved oxygen deficit increases as the fresh water velocity
increases, with the maximum deficit occurring at the station when the
flow velocity is in the order of two (2) miles per day.
The specific planes presented in this report were developed from
consideration of an infinite estuary. The geometry of the estuary may
have significant effects on the shape of the planes and in general, on
the sensitivity and behavior of the system.
An example of the influence of geometry on system sensitivity and
behavior is illustrated in Figures 18 and 19. In these Figures, a
four (4) section system is considered with no fresh water flow. In
Figure 18, the point of waste introduction is located between the two
- 69 -
-------
middle sections at X = 0. Case #1 on this Figure illustrated the cal-
culated dissolved oxygen deficit curve when each of the four sections
has the same cross-sectional area. In Case #2 of this Figure, the
cross-sectional area of the extreme section has been increased by a
factor of four (4). It can be observed that for Case #2, the entire
deficit profile is significantly reduced in magnitude and the maximum
deficit occurs behind the point of waste introduction, and away from
the section with the large cross-sectional area. Figure,19 presents
profiles for comparable system geometries. In this instance, the
waste input location is bounded by an infinity section. It should be
noted that for Case #1, in both Figures 18 and ,19 when the estuary
geometry is comparable in all sections, the location of the waste input
does not influence the magnitude of the maximum deficit, which is 5 ppm
and occurs at the point of waste input. In Case #2 for both Figures
13 and J9, the maximum deficit occurs behind the waste input, away from
the section with the large cross-sectional area. Note that the magni-
tude of the maximum deficit is lower in Case #2, as the waste load lo-
cation approaches the section with the large cross-sectional area.
Several points of information on sensitivity can be made from an
analysis of the information presented in Figures 18 and 19
1. The maximum dissolved oxygen deficit occurs at the point
of waste discharge for symmetrical systems with zero fresh
water flow.
2. The maximum deficit need not occur at the point of waste
discharge for unsymmetrical systems with zero flow.
3. The calculated dissolved oxygen deficit from a constant
waste input decreases as the liquid volume of the system
increases.
Figure 20 presents a plot of the plane relating the deficit to the
deoxygenation "DK" and reaeration "AK" coefficient. Note the marked
reduction in deficit associated with increasing reaeration coefficient.
L'Hospital's rule is applicable to the range where AK = RK as the par-
ticular integral in the deficit solution contains the difference in
the two coefficients (DK/AK-RK). The increase in deficit due to an in-
crease in the deoxygenation coefficient is also indicated on the Figure.
The range of the deoxygenation coefficient DK on Figure 20 is from 0.05
to 0.5 per day.
The previous discussions have indicated in a limited manner, the
response of estuary systems to changes in parameters, inputs and geometry.
The area of system sensitivity is relatively complex and requires addi-
tional effort before adequate definitions are available to define the
anticipated behavior of the waterways which are being studied in the pre-
sent program.
- 70 -
-------
It is suggested that as part of the analysis and use of the mathe-
matical models for water quality a systemic effort be made to develop
an understanding of the behavior and sensitivity of the system. Speci-
fically, it is suggested that for each of the models the following
program be developed:
1. Obtain first, best estimates of all system parameters,
inputs and geometry.
2. Using a unit load in all sections, develop profiles
employing the best estimates obtained in No. 1 above.
3. Repeat the procedure in item No. 2 for ranges of all
the system parameters. The best estimates may initially
be halved and doubled. (As an example, E = 10 best
estimate, also evaluate E = 20 and E = 5.)
A major challenge remains in developing a convenient graphical
method of presenting and displaying this information on system behavior
and sensitivity. This report has employed individual profiles and three-
dimensional graphs. In all cases, room for improvement in techniques
for displaying system sensitivity remain. The program for evaluating
system sensitivity should include provisions to provide information in
all three of the significant areas of application. Namely:
1. Determination of the variables and parameters which
require refined field measurements because the system
is sensitive to these items. This will provide infor-
mation on what must be measured, and the required
accuracy.
2. Determination of changes in calculated profiles to
assist in developing a consistent set of parameters
defining all of the phenomena having a major influence
on water quality.
3. Determination of the feasibility of alternate methods
of improving water quality based upon the systems'
response to changes in parameters, inputs, and geometry.
- 71 -
-------
SENSITIVITY PLANES
O.| X=- 8.0 MILES
2
FIGURE
- 72 -
-------
X= -6.0 MILES
X= -4.0 MILES
FIGURE 14
- 73 -
-------
FIGURE' 15
- 74 -
-------
D X = 4.0 MILES
U
30
FIGURE 16
- 75 -
-------
X= 6.0 MILES
U
FIGURE 17
- 76 -
-------
ZERO FLOW ANALYSIS
E
Q.
Q.
o 6
f.
UJ
X
o
Q 3
UJ
o 2
^
a?
Q
0
J I
10 8
20246
MILES
8 10
FIGURE is
- 77 -
-------
ZERO FLOW ANALYSIS
7
LJ
CD
x 4
o
LL! o
3
82
O)
0
CASE
CASE
W
„ -t
FLOW = 0
CASE I
Z
CASE 2
i i i
10 8
202
MILES
8 10
FIGURE 19
-------
DEOXYGENATION
COEFFICIENT
\/0 REOXYGENATION
Vb. COEFFICIENT
a
a
o
§
z
UJ
X
O
Q
LJ
21
o
en
CO
Q
E 10 Mi2/DAY
U .8 Mi/DAY
X = 0 Mi
DK( AXISK I/DAY)
INCREASE
INCREASING
ZONE
'OF La HOSPI
TALS RULE
APPLICATION
0.05
0.5
AK-(I/DAY)
FIGURE 20-
- 79 -
-------
DESCRIPTION OF THE COMPUTER PROGRAM
- 80 -
-------
ID JOB, IDATE,
IS2.IS1.NN,
IND, ICON,
FAC1.FAC2,
FAC3, ITYPE
Or-^T .
n
r~
DMI (I ,'j)'="
DM(J,I)
/ J *
2— ( J <
\ J*
yes
J+l \
JBM? \
1 /
DSOL(I)=
DFOR(I)
/I*
— H i 5
yes
1+1 \
JBM? )
•^—/
r
( SUBROUTINE DRIGHT
r
'
'
SUBROUTINE MATS
'
IJK=2
0
I
^_^ (^ SUBROUTINE MATINV J
) "
I
I
( SUBROUTINE BCGH J
)
)
print ITY
IDJ
1 ^
PE, '
OB, II
I
T SUBROUTINE PRI J
I
f SUBROUTINE PRII J
I
( SUBROUTINE PRJ )
N=
0 *
/~~i?T~\
' *\ ii.
\u±
IMAXJ) M
t+1 7
yes
r SUBROUTINE BDOUT J
t
SUBROUTINE BCGH
)
print
read
IDJOB ,
IDATE
IS2,
IS1.NN,
IND,
ICON,
FAC1,
FAC2,
FAC3,
ITYPE
/ IDJOB , IDATE
152 IS1 NN
IND.NOTH,
FAC1.FAC2,
FAC3, ITYPE
/-\
ll
SUBROUTINE MAMS
SUBROUTINE BCGH
SUBROUTINE DRIGHT
FIGURE 21: FLOW CHART OF MAIN PROGRAM OF ES001
81 -
-------
The ES001 computer program was written for an IBM 370 machine in
FORTRAN IV. In this section, the program will be described as follows:
1) a literal description of the main program accompanied by
a flow diagram will be presented;
2) individual literal descriptions of each of the subroutines
and functions used to support the main program will follow in
alphabetical order;
MAIN PROGRAM
A flow chart of the main program is shown in Figure 21 and a listing
in the appendix. As can be seen there are several branches due to a
number of decision operations which are made in the course of ESOOl's
execution. Most of these are due to the fact that the program has the
capability of modelling several general systems, as well as that it has
the option to test the response of a particular system to several sets
of loads by saving the inverted matrix from the first run in a series.
Since this is so, a straight description of the entire MAIN would be
cumbersome and might be confusing. To alleviate this, the description
will be presented in two parts. First, the program with the matrix
computations is described, and the section of the flow chart dealing
with it will be presented. Second, a run with no matrices computed will
be explained and the section of the main flow chart associated with it
will be presented.
[NOTE: A listing of ES002 and a discussion of that programs
specifications is included in the appendix. ES002
is designed to run on an IBM 1130 System. -For
this reason, some of its code will differ from that
of ES001 which is designed to run on an IBM 370
system. Therefore, the following sections will not
necessarily be an accurate description of the ES002
algorithm. However, as a large part of the programs
are similar this discussion may prove useful in working
with ES002.
As noted previously, all input, output and internal
computation for the programs are identical.]
- 82 -
-------
.
(
(.
c_
START
DMI(I,J)=
JBM =10° DM(J,I)
read < > ^
/IDJOB, IDATE,
IS2.IS1.NN,
IND, ICON,
FAC1.FAC2,
FAC3, ITYPE
/
-<
s
B
D
I "l" \
I± 2507 JV^-i
I -t-I+1 /
• 1
DX1(I)=0.
*
Pr
IMi
JMA
Lnt ^
IDJOB,
IDATE,
IS2,
IS1.NN
IND,
ICON,
FAC1,
FAC2,
FAC3,
ITYPE
y™.
1
^X~IS2+IS1
<=2*IS2+IS1
1
SUBROU1TNL DAIAR )
pri
-it |
ITYPE,
IDJOB ,
IDATEJ
SUBROUTINE PRI J
k
SUBROUTINE PRJ J
( SUBROUTINE DATAC J
<:
<
J print
yes
/ J ''•J+l \
1 no ( J
V j A / /V
1
(^ SUBROUTINE MATINV ^}
1
( SUBROUTINE BCGH J
r *<
print ITYPE,
IDJOB,
IDATE
( SUBROUTINE PRI J
1
( SUBROUTINE PRII )
L
( SUBROUTINE PRJ }
N= 0 *
/ I- 1 \
yes
( SUBROUTINE BDOUT J
*
yes
ICON +1 ;
print read | 3
IDJOB, / IDJOB, IDATE
U IDATE, « IS2,IS1,NN,
BSOL(I)= IS1,NN, FAC1.FAC2,
RFORfl^ IND, FAC3, ITYPE
no / i +-i+tf\
\ T -*- 1 /
FAC1,
FAC2,
FAC3,
ITYPE
U\
!i2 • NFLAG f" \ »
*^~/^ ^—1 ( SUBROUTINE BRIGHT 1
^Ajes
t
^,N< Q^, ygs J IJK= 1 J[ SUBROUTINE MATS ")
^= ^^ \ ^ -^
FIGURE 22: FLOW CHART OF MAIN PROGRAM RUN WITH MATRIX COMPUTATIONS
- 83
-------
1) Run with matrix computations (flow chart on page 83 ).
The program starts by setting JBM, which is the maximum allowable
number of boundary conditions, to 100. Then, the Job identifica-
tion card is read (see page 113 for description). After a do loop
seta BBXI and DDXI to zero, the values that were read from the
identification card are printed (as described on page!24_). After
IMAX (total number of sections) and JMAX (total number of boundary
conditions) are calculated, subroutine DATAR is called in order to
read in and print out the remaining data cards (as they appear as
card images) as well as to check the data for errors and print a
number of diagnostics. If such errors exist a value of the parameter
NFLAG >0 is returned to MAIN.
Subroutines PRI and PRJ, which list the data in tabular form, are
called followed by DATAC which sets up temperature correction factors,
conversion factors and does miscellaneous computations. At this
point, if NFLAG is greater than zero the program essentially bombs
out by returning to MAIN 014 where it reads the next set of data or
ends if no data set is available.
If no errors were discovered (NFLAG = 0) the program proceeds to
another decision where it is determined whether the model is a new
one for which new matrices should be computed and inverted, or one
for which the matrices have been previously computed but to which
different loads are to be applied.
If the latter is true the program merely omits the steps in which
the matrices are set up and inverted. (The actual mechanism is
described in the next section.)
If the former is true and this is a new model (NN = 0), the program
sets UK = 1. This index connotes that the succeeding steps are
for BOD. Subroutine MATS, which sets up the matrix (in this case
for BOD) and subroutine BRIGHT which sets up the BOD forcing function
are called followed by a do loop which redefines BFOR (BOD forcing
function), which was output from BRIGHT, as BSOL. It also transposes
the rows and columns of the matrix in preparation for the matrix
inverse. MATINV then inverts the matrix and solves for the unknown
coefficients. For output purposes BCGH then splits the coefficient
matrix into B's and C's.
At this point, the BOD portion of the solution has been essentially
computed. A decision is now made to determine whether a DO solution
is desired. If it is, IJK is set to 2, MATS is called to set up the deficit
matrix and DRIGHT generates the deficit forcing function. After a do loop
transforms DFOR to DSOL and transposes the elements of the deficit matrix
(as was done in the BOD), the matrix is inverted and BCGH is applied
to the solution of coefficients. If a deficit solution had not been
- 84 -
-------
desired, the program would skip from MAIN 052 to this point, thus,
omitting the deficit solution.
The particular model labels for this run are then printed and subroutine
PRI, PRII, and PRJ print the input data-in converted (or compatible) units
as well as several computed parameters (e.g., velocity). After this
N is set to zero and a do loop is used to print out the concentration
profile for each section. At this point if ICON is less than 2
(that is, if we are not dealing with a coupled system with two inputs
like carbonaceous and nitrogenous BOD), the program returns to
MAIN 014 to begin another model run. If we are dealing with a
coupled system with two inputs, ICON is incremented and a computed
go to statement is employed to direct the flow (in this case ICON
would now equal 3) to read and print the set of job identification
cards for the second input and returns to MAIN 026 and computes
and prints the BOD and DO solution as done previously and eventually
returns to the decision at MAIN 070. There it would be decided that
ICON would be incremented to 4. Then, using the computed go to of
MAIN 072 the program flows to MAIN 067 where the sum of the profiles
due to the two inputs is printed using BDOUT. ICON is finally incre-
mented to 5 and the computed go to returns the flow to MAIN 014 where
another run is started.
- 85 -
-------
START
read
ID JOB., IDATE,
IS2,IS1,NN,
IND, ICON,
FAC1,'FAC2,
FAC3 , ITYPE
!<. 250?
print ,
SUBROUTINE DATAR
print
ITYPE,
IDJOB,
IDAT
SUBROUTINE PRI
SUBROUTINE PRJ
yes
SUBROUTINE BRIGHT
SUBROUTINE MAMS
)"~*C
SUBROUTINE BCGH
J
print
ITYPE,
IDJOB,
C
c
c
SUBROUTINE PRI
SUBROUTINE PRII
D
SUBROUTINE PRJ
C
/ I* 1 \
K I<_ TMAXTf-uii —
\ Tn+i 7
1 yes
SUBROUTINE BDOUT
)
*
ICON =
ICON +1
^ 1>2>5 JL^ 4
print reatj |
O-
1 IDJOB,
IDATE ,
IS2,
IS1.NN,
IND,
ICON, .
FAC1,
FAC2,
FAC3,
ITYPE
^ — i
/ IDJOB , IDATE
( igT isi NN
IND.NOTH,
FAC1.FAC2,
FAC3, ITYPE
vO
f SUBROUTINE MAMS J
T
C SUBROUTINE DRIGHT ^
SUBROUTINE BCGH
FIGURE 23: FLOW CHART OF MAIN PROGRAM RUN WITHOUT MATRIX COMPUTATIONS
- 86
-------
2) Run without matrix computation (flow chart on page 86). This
is the case where succeeding runs are being executed on the same
system with different loads.
As can be seen from the flow chart the basic difference between this
and the previous flow is that this run skips the matrix computations.
This would occur in a case where several runs with different loads
would be run in succession and recomputation of the matrix would be
unnecessary. After the decision at MAIN 041, the program checks if
ICON is greater than 1. If this is so, the run is terminated, since
only one matrix can be saved by this program at a time (if it is
dealing with a double input system it would only save the second matrix
and erase the first). If ICON denotes that it is a single system or a
coupled system of one input, the program proceeds to BRIGHT, MAMS AND
BCGH which essentially set up and print the BOD matrix and forcing
function. If the system is single (ICON = 1) the program proceeds to
MAIN 063 where the data is printed in compatible units, the profile is
outputted and the flow is returned to MAIN 014 where the next model is
read. If the system is coupled (e.g., BOD and DO), DRIGHT, MAMS and BCGH
are applied to generate and print the DO deficit matrix and forcing function.
The da-ta is then printed in compatible units, the profiles are output
and the flow is returned to MAIN 014 where the next model is begun.
- 87 -
-------
SUBROUTINES
Subroutine BCGH
MATINV outputs the solution for the unknown coefficients in the
form of a single column of numbers. The purpose of BCGH is to trans-
form this column into sets of coefficients corresponding to each section
for output in PRI.
As was noted in the theory section . the mathematical
formulation for each double section generates two unknown coefficients
while each single section generates one coefficient. Utilizing this
fact, and the fact that the double sections have lower identification
numbers according to the definitions of this program (see page 118)t
BCGH converts the single column of coefficients into pairs of subscripted
B's and C's (for BOD) or G's and H's (for DO deficit) for each double
section. After the double sections are completed, the subroutine con-
tinues to designate the coefficients according to the following standards:
KFO
2 (- -)
3 (+ »)
4 (bay)
Coefficient
BOD
B
C
B
DO
Deficit
G
H
G
Subroutine BDOUT
BDOUT prints the concentration profile of the substance(s) being
modelled for a particular section (page 127). Values are printed
within the section for intervals of length PDEL or 1 mile if no PDEL
is stipulated.
Use is made of BODO, a subroutine which actually computes the con-
centration for each point and stores values of BOD and DO for coupled
systems with two components.
Subroutine BODO
Utilizing functions OM and YB, this subroutine generates and com-
putes the final BOD and DO deficit equations (as on page 27 ).
- 88 -
-------
The equations are calculated as follows:
1. BOD:
BOD = BBB =» OM (KFO,B,C,FR,FBR,I,X) + YB (I,X)
where
OM - is a function described on page 101. In this
particular case Oil is the general solution to
the BOD differential equation (equation 11,
page 19 ).
YB - is a function described on page 102 which
equals the particular solu-
tion to the BOD differential equation (equation
12a, page 21_) .
2. DO Deficit:
a) For double sections and infinity sections.
DO deficit = ODD = OM (KFO,G,H,FA,FBA,I,X) + ZZ
+ FAG(I) * OM (KFO,B,C,FR,FBR,I,X)
where
OM (KFO,G,H,FA,FBA,I,X) = In this case is the general
solution to the DO deficit differential equation.
ZZ + FAC(I) * OM (KFO,B,C,FR,FBR,I,X) - This terra is
the particular solution to the DO deficit general
solution. ZZ and FAC(I) were computed in DATAC
(page $2 ).
b) For completely mixed bays.
DO deficit = ODD = OM (KFO,G,II,FA,FBA,I,X)
Subroutine BRIGHT
BRIGHT sets up the BOD forcing functions for eac'i of the boundary
conditions. Utilizing a computed GO TO, the subroutine sets up the
following forcing functions for the appropriate boundary condition
(as delineated in Table 6 on page 62 ):
- 89 -
-------
Boundary Condition Forcing Function
1, 2, 3 YB(K,Y) - YB(I,X)
4 - YB(I,X)
5, 6, 7 W(J) + Q(I) * YB(I,X) - Q(K) * YB(K,Y)
W(J) + A(I) * YB(I,X) - Q(K) * YB(K,Y)
- Q(KK) * YB(KK,Z)
9 W(J) + Q(I) * YB(I,X) + Q(KK) * YB(KK,Z)
- Q(K) * YB(K,Y)
10 QD(J) * BDD(J) - Q(K) * YB(K,Y)
11 W(J) - Q(I) * YB(I,X)
12 W(J) - Q(I) * YB(I,X) - Q(K) * YB(K,Y)
13 W(J) + Q(I) * YB(K,Y)
14 W(J) + Q(I) * YB(K,Y) - Q(K) * YB(K,Y)
where:
W(J) = The point load at the junction in question.
Q(I,K or KK) = The advective flow in the subscripted section.
YB(I,K or KK; = This function is called to give the particular
X, Y or Z) solution for the BOD equation as delineated in
Table 3. on page 21_- The algorithm of function
YB is explained on page 104,
QD(J) = The advective flow over a dam.
BDD(J) = The BOD in the flow over a dam.
I,K or KK = Subscripts used to differentiate the sections
which make up a junction.
Since the cards must alxrays be inputted in order
(either upstream or downstream) I usually connotes
the first section, K the second and KK the third
section forming a junction.
- 90 -
-------
X, Y or Z = The distance of the junction from its reference
junction for each segment making up the junction.
X corresponds to segment I, Y to segment K, Z to
segment KK. (They also correspond to XSA, XSB and
XSC as noted on page 116 of the input description).
Finally, the forcing functions are printed.
Subroutine DATAC
This subroutine performs the following operations:
A. Temperature Correction Factors
AK(I) = AK(I) * FAC1 ** (TEMP(I) - 20.)
DK(I) = DK(I) * FAC2 ** (TEMP(I) - 20.)
RK(I) = RK(I) * FAC3 ** (TEMP(I) - 20.)
where
FAC1, FAC2 and FAC3 are temperature correction factors
which were read in on the job identification card.
B. Conversion Factors to Compatible Units
1. The following conversion factors are defined:
a) CA =• 1.0/(5280.0**2) - square feet to square miles
b) CV = l./(5280.0**3) - cubic feet to cubic miles
c) CF = 3600.0*24.O/(5280.**3) - cfs to cubic mile/day
d) CL = (l.OE + 6)/(62.4*5280.**3)
-Ibs/day to (mi3/day) (mg/1) for point loads
-Ibs/mile/day to (mi2/day) (mg/1) for uniform loads
e) CM = 12./39.37 - feet to meters
2. These factors are then applied as follows:
a) AREA(I) = AREA(I) * CA
b) Q(I) = Q(I) * CF
c) WU(I) = WU(I) * CL
- 91 -
-------
d) HT(I) = HT(I) * CM
e) VOL(I) = VOL(I) * CV
f) W(I) = W(I) * CL
g) QD(I) = QD(I) * CF
C. Miscellaneous variables are calculated:
a. FAC(I) = part of DO deficit forcing function -
DK(I) * FF(I)/(AK(I) - RK(I))
b. ZZ(I) = ((BD(I)/HT(I) + PR(I)
+ DK(I) * FF(I)/(RK(I) * AREA(I))/AK(I)
D. Forms and prints the symbolic BOD and DO deficit matrices
using PLACE.
E. S, V, SD, VD are calculated and placed in common.
F. If present calculation is only for DO, sets BOD to 0.
Subroutine DATAR
This subroutine reads the data cards and prints them as described
in the input description (pages 113 to 117). In the process of doing
this, it also searches for input errors, prints an appropriate error
message and aborts the run. Listed below are the errors and their
associated messages:
Error Message
Section #<0* Program aborts - no message printed
KFO is incon- Program aborts - no message printed
sistent with
section type
(double or
single)*
AK(I) = RK(I) ERROR IN AK.OR.RK followed by the section
number in which this occurred and the value
of AK and RK.
*Uses NTPS.
- 92 -
-------
Error Message
HT = 0 Prints ERROR IN HT followed by the section mmb«r
Junction cards Prints junction info followed by FAIL. If
are incorrectly correct prints PASS.
typed
Finally, the abbreviations for the 14 boundary conditions are assigned
via a data statement (as explained on page 125).
Subroutine BRIGHT
BRIGHT sets up the dissolved oxygen deficit forcing functions for
each of the boundary conditions. Utilizing a computed GO TO, the sub-
routine sets up the following forcing functions for the appropriate
boundary condition (as delineated in Table 6 , page 62_) '•
Boundary Condition Function
1, 2, 3 YD(K,Y) - YD(I,X)
4 - YD(I,X)
5, 6, 7 Q(K) A YD(I,X) - YDP(I,X) + YDP(K,Y)
- Q(K) A YD(K,Y)
8 Q(I) A YD(I,X) - YDP(I,X) + YDP(K,Y)
- Q(K) A YD(K,Y) + YDP(KK,Z) - Q(KK) * YD(KK.Z)
9 Q(I) A YD(I,X) - YDP(I,X) + YDP(K.Y)
- 0(K) A YD(K,Y) + Q(KK) * YD(KKZ) - YDP(KK.Z)
10 YDP(K.Y) - Q(K) A YD(K.Y) + QD(J) A DD(J)
11 YDP(I.X) - Q(I) A YD(I,X) + FF(K) A DK(K) A
B(K) A VOL(K)
12 YDP(I,X) - Q(K) A YD(I,X) + FF(KK) A
DK(KK) A B(KK) A VOL(KK) + YDP(K.Y)
A YD(K,Y)
- 93 -
-------
Boundary Condition Function
13 Q(I) * YD(I,X) - YDP(I.X)
+ FF(K) * DK(K) * B(K) * VOL(K)
1* Q(D * YD(I,X) - YDP(I.X)
+ FF(KK) * DK(KK) * B(KK) * VOL(KK)
+ Q(K) * YD(K,Y) - YDP(K,Y)
where:
Q(I,K, or KK) = The advective flow in the subscripted section.
FF(I,K or KK) = The ratio of the ultimate to the 5-day BOD.
DK(K or KK) = Deoxygenation rate.
VOL (K or KK) = Volume of section.
OD(J) = Flow over a dam.
DD(J) = DO deficit over a dam.
B (K or KK) = Constant which was generated in the BOD solution.
YD(I,K, or KK; = Function which generates the particular solution
X,Y,Z) to the DO deficit equation (as shown in Table 1 ,
on page 27 )•
YDP(I,K or KK: = Function which generates the derivative of the
X,Y or Z) particular solution to the DO deficit equation
multiplied times the dispersion coefficient and
the area.
Subroutine MAMS
This is only called if consecutive runs are being made with different
loadings. In such a case, the forming and inversion of the matrices
would be suppressed since these inverted matrices would already be in
core. MAMS merely multiplies the inverted matrices by the new forcing
functions which would be generated by the new loads.
- 94 -
-------
Subroutine MATINV
This is a standard matrix inversion subroutine using an elimination
mtthod to invert the matrix and calculate the unknown coefficients
utilizing the forcing function. If IND (printing indicator) is set to
2, it also prints out the original matrix and the forcing function and
then the inverted matrix and the coefficients.
Subroutine MATS
This subroutine sets up the individual elements of the BOD or DO
deficit matrix. It does this by evaluating each of the boundary condi-
tions, using NOSEC to subscript the column location of the elements.
In the case of double sections, two subscripts are allocated and in the
case of single sections one subscript is designated. The following
elements are computed for the associated boundary condition:
Boundary Condition
TM(IM,J)
TM(L,J)
TM(KM,J)
TM(M,J)
TM(L.J)
TM(KM,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(KM,J)
TM(M,J)
TM(L,J)
TM(KM,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(M,J)
Matrix Elements
FB(I,X)
F(K,Y)
FB(K,Y)
F(I,X)
F(K,Y)
FB(K,Y)
FB(I,X)
FB(K,Y)
FB(I,X)
1.0
FP(I,X) - Q(I) * F(I,X)
FBP(I,X) - Q(I) * FB(I,X)
Q(K) * F(K,Y) - FP(K,Y)
Q(K) * FB(K,Y) - FBP(K,Y)
FP(I,X) -
Q(K) * F(K,Y) - FP(K,Y)
Q(K) * FB(K,Y) - FBP(K.Y)
FP(I,X) -
FBP(I.X) - Q(I) * FB(L.X)
Q(K) * FB(K,Y) - FBP(K.Y)
- 95 -
-------
Boundary Condition
8
Matrix Elements
10
11
11 - BOD
11 - DO deficit
12
12 - BOD
12 - DO deficit
13
13 - BOD
13 - DO deficit
14
TM(IM,J)
TM(L,J)
TM(KM,J)
TM(M,J)
TM(N,J)
TM(KKM.J)
TM(IM,J)
TM(L,J)
TM(KM,J)
TM(M,J)
TM(KKM,J)
TM(N,J)
TM(KM,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(M,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(KM,J)
TM(M,J)
TM(N,J)
TM(N,J)
TM(IM,J)
TM(L,J)
TM(M,J)
TM(M,J)
TM(IM,J)
TM(L,J)
TM(KM,J)
TM(M,J)
FP(I,X) - Q(I) *
FBP(I,X) - Q(I) * FB(I,X)
Q(K) * F(K,Y) - FP(K,Y)
Q(K) * FB(K,Y) - FBP(K,Y)
Q(KK) * FB(KK,Z) - FBP(KK,Z)
Q(KK) * F(KK,Z) - FP(KK,Z)
FP(I,X) - Q(I) * F(I,X)
FBP(I.X) - Q(I) * FB(I,X)
Q(K) * F(K,Y) - FP(K,Y)
Q(K) * FB(K,Y) - FBP(K,Y)
FP(KK,Z) - Q(KK) * F(KK,Z)
FBP(KK,Z) - Q(KK) * FB(KK,Z)
Q(K) * F(K,Y) - FP(K,Y)
Q(K) * FB(K,Y) - FBP(K.Y)
* F(I,X) - FP(I,X)
Q(I) * FB(I,X) - FBP(I.X)
VOL(K) * RK(K)
VOL(K) * AK(K)
) - FP(I,X)
Q(I) * FB(I,X) - FBP(I.X)
Q(K) * F(K,Y) - FP(K,Y)
Q(K) * FB(K,Y) - FBP(K.Y)
VOL(I
-------
Boundary Condition Matrix Elements
14 ~ BOD TM(N,J) Q(KK) + VOL(KK) * (TIDEl(KK) +
RK(KK))
14 - DO deficit TM(N,J) = Q(KK) + VOL(KK) * (TIDEl(KK) +
AK(KK))
where
Q(I,K or KK) = The flow in a section.
VOL(K or KK) = The volume of a bay.
TIDE1(K or KK) = The tidal exchange coefficient of the section.
RK(K or KK) = BOD removal rate.
AK(K or KK) = Reaeration rate.
F, FB, FP, FBP are functions which are accessed via EXTERNAL and ENTRY
statements from FUN and PFUN.
For BOD*
F = eSX
FB = e™
FP = SeSX EA
FBP = Ve^ EA
For DO*
F = eSDX
FB =
FP = SDeSDX EA
FBP = VDe EA
*See pages 50 & 54
- 97 -
-------
Subroutine PLACE
This subroutine is called by DATAC to form the column location of the
elements of the symbolic matrix which is printed as part of the output
(p.!26_J. The rows correspond to the particular boundary conditions while
the columns are in terms of the coefficients to which the particular element
applies (this is determined using NOSEC).
Subroutine PRI
This subroutine prints parameter values as described in the output
description (p. 126.). The data is printed in a more structured form
than when originally printed via DATAR and headings are provided. It
consists of three portions which are accessed by individual entry state-
ments :
1) PRI - Prints the section data in either original or compatible
units.
2) PRII - Prints section parameters which have been computed within
program (e.g., S, V, U, ZZ, etc.).
3) PRJ - Prints the junction data in either original or compatible
units.
Subroutine SVD
This subroutine calculates the coefficients (S, V, SD, VD) in the
exponentials of the BOD and DO deficit solution of the differential
equations.
a. For the BOD solution the equations are as follows:
(1) When the velocity U = 0.
S = RK/E V = RK/E
(2) When the velocity U ^ 0.
S = U/2E (1+1 1 +
V = U/2E (1 - I 1 +
>4 U'
- 98 -
-------
b. For the DO deficit solution the equations remain the same as
above, except the quantities SD and VD are calculated and AK
replaces RK.
- 99 -
-------
FUNCTIONS
Function FUN
FUN utilizes entry statements to route various functions to PFUN
aa follows:
Entry FUN
FR(I,X) F(X,V,I,X)
FA(I,X) F(SD,VD,I,X)
FBR(I.X) FB(S,V,I,X)
FBA(I.X) FB(SD,VD,I,X)
FPR(I,X) FP(S,V,I,X)
FPA(I.X) FP(SD,VD,I,X)
FBPR(I.X) FBP(S,V,I,X)
FBPA(I.X) FBP(SD,VD,I,S)
YBP(I.X) YBPT(S,V,I,X)
The values of FUN are found in PFUN.
Function NOSECT
This is used to designate a subscript for each section so that any
operation associated with that section will be consistently labelled.
If the section is a double section
NOSECT = 2 * IL
where
IL = The section number (ISA, ISB, ISC as defined on page 116
of the input description).
If the section is a single section
NOSECT = 2 * IS2 + (IL - IS2)
- 100 -
-------
where
IS2 = The number of double sections.
Function NTPS
This function (called by subroutine DATAR) checks whether the sec-
tions are numbered correctly and whether the sections have a proper
value of KFO.
It checks for the following convention, and if they are not met, it
returns to DATAR with an appropriate value (>4) of NTPS which in turn
causes the main program to terminate the present run:
1) Checks if the section number >0.
2) Section numbers for double sections < IS2 (the total number of
double sections).
3) KFO for double sections = 1.
4) Section number for single sections
-------
where
BG(I) = B of the BOD solution or G or DO deficit.
CH(I) = C of the BOD solution of II of DO deficit.
F(I,X) and FB(I,X) = externalized function which through FUN and
PFUN result in the appropriate exponential
terms for the BOD or DO deficit solution or
its derivative.
Function PFUN
This function subprogram calculates the following values which are
accessed by entry statements:
Entry PFUN
F(SSD,WD,IX) EXP(X * SSD(I))
FB(SSD,WD,I,X) EXP(X * VVD(I))
YBT(SSD,WD,I,X) WU(I)/(AREA(I) * RK(I))
FP(SSD,WD,I,X) SSD(I) * EXP(X * SSD(I)) * E(I) * AREA(I)
FBP(SSD,WD,I,X) WD(I) * EXP(X * VVD(I)) * E(I) * AREA(I)
YBPT(SSD,WD,I,X) 0
These functions are used throughout ES001 but particularly in MATS
in calculating the elements of the matrices.
Function YB
YB = WU(I)/(AREA(I) * RK(I))
where
WU(I) = Uniform load to a section
AREA(I) = Cross-sectional area of the sectional
RK(I) = BOD removal rate of the section.
- 102 -
-------
Function YD
YD calculates the particular solution to the DO deficit equations
(as shown in Table 1 on page 27 ).
YD = FAC(I) * OM(KFO,B,C,FR,FBR,I,X) + ZZ(I)
where
FAC(I)* = DK(I) * FF(I)/(AK(I) - RK(I))
ZZ(I)* = ((BD(I)/HT(I)) + PR(I) + (DK(I) * FF(I)/(RK(I)
* AREA(I))/AK(I)
OM(KFO,B,C,FR,FBR,I,X) = a function (page 1Q1) which sets up the
BOD solution via FUN and PFUN.
Function YDP
YDP sets up the first derivative of the particular solution to the
DO deficit equations (as shown in Table 5 on page 54 ).
YDP = FAC(I) * OM(KFO,B,C,FPR,FBPR,I,X)
where
FAC(I) = Same as above in write-up of function YD.
OM(KFO,B,C,FPR,FBPR,I,X) = a function (page,101) -which sets up
the derivative of the BOD solution
via FUN and PFUN.
*NOTE: FAC(I) and ZZ(I) were calculated in DATAC.
- 103 -
-------
USER'S MANUAL
- 104 -
-------
Initial Steps in Formulating ES001
The Initial step in developing an analysis of water quality (in this
case BOD and DO deficit in particular) is a critical review of the avail-
able information on the body of water under study. Specifically, the
following physical, geophysical, and hydrological information must be
obtained and evaluated.
1. The mean water depth of the various bodies of water
should be determined.
2. The cross-sectional area of the various water bodies
should be plotted as a function of location.
3. The drainage area should be evaluated as a function
of distance, which provides an insight into fresh
water flow quantities. As an example, if a small
drainage area is associated with a section of estuary,
it is unrealistic to anticipate large fresh water
flows from run-off and normal ground water.
4. A qualitative evaluation should be made of the hydro-
graph of the system. This will assist in determining
when field sampling programs might be carried out.
5. If the fresh water flow rate does not significantly
change with distance, then the advective velocity of
the system can be calculated. On the other hand, if
area changes and/or fresh water flow inputs are such
that the system advective velocity varies, plots of
velocity with distance should be developed.
6. Knowing the tidal velocity in estuaries a plot of
"AK", the reaeration coefficient vs. distance, using
equation (49) to calculate "AK", is constructed.
7. The confluence of major water bodies should be
noted on the distance plots.
8. The BOD removal and oxidation coefficients RK and
DK should be evaluated from appropriate information.
9. The dispersion coefficient should be assigned, or
evaluated.
10. Waste sources should be located and indications of
the magnitude of loads obtained.
- 105 -
-------
11. Chloride or other conservative substance profiles
should be constructed. BOD and dissolved oxygen
deficit profiles as a function of distance should
be developed from any existing data.
It is usually advisable as a first step to construct a model incor-
porating as many simplifying assumptions as practical.
Having assessed the reliability of the available information on the
physical hydrodynamic and input data, the various distance plots are
compared in an attempt to minimize the number of segments required.
Changes in area and flow may assumed to occur at the confluence of two
bodies of water, thus reducing the number of segments required. Adjacent
waste inputs can be combined at their center of gravity to reduce the
number of segments. All of these steps are contingent upon the avail-
ability of data and the degree of complexity of the model justified by
the available information (additional discussion of segmentation can be
found on page 29 ) •
Discussion of System Parameters:
It is possible to evaluate the various parameters, U, E, RK, DK, R,
and AK from physical, hydrodynamic and observed data.
Advective Velocity (U)
The flow velocity which transports material out of the estuary can
be evaluated :
U = ...................... (46)
Dispersion (E)
A plot of the log of the concentration of a conservative substance
(usually chlorides) vs. distance results in a straight line. The
slope of the line is equal to:
Slope = U/E
therefore;
E = U/slope ................... (47)
- 106 -
-------
BOD Removal and/or Deoxygenation (RK and/or DK)
A plot of the log of the BOD vs. distance results in a straight line,
The slope of the line is defined by:
TT L . w . F
Slope = £ [1 ± ( 1 + fL_R|_E } ........
The value of RK can be obtained from the solution of the above equa-
tion.
If the rate of BOD removal changes, the semi-log plot of BOD vs.
distance will yield two (2) straight lines having different slopes.
The slope of the line adjacent to the waste input will provide a
measure of RK. The slope of the line farthest from the waste input
is generally employed to calculate DK. The calculated value of DK
must be equal to or less than RK. This analysis for calculating
the values of RK and DK is valid only for BOD data which results
from a single waste input. If the BOD profile is influenced by
several waste inputs which cannot be grouped into a single source
of BOD, the above analysis is not valid. A reasonable first assump-
tion for the range of values of RK and DK is 0.20 - 0.50 to the base
"e", with the two parameters assumed to be equal.
Reaeration Coefficient
This parameter can be calculated as follows :
AK = D ' U ................. (49)
H '
D = diffusivity in Equation 49
The velocity is the average velocity over a tidal cycle in estuaries.
Temperature Effects
Temperature influences the parameters AK, RK, and DK. Generally
accepted temperature corrections for these parameters are:
' "cao'c) 1-«2 ' ............ (50>
- 107 -
-------
RK = RK 1.04(t " 2°) (51)
(T) (20°C)
DK = DK 1.04^ ~ 2°) (52)
(T) (20°C)
RK = Nitrogenous BOD RK = RK 1.08^ ~ 2°) (53)
N(20°C)
In addition, temperature will also significantly influence the
saturation value of dissolved oxygen. Tables are available
indicating temperature and dissolved solids effects on the
saturation value of oxygen in water.
Tidal Exchange Coefficient
The tidal exchange, "R", provides a method of evaluating the trans-
port of material out of a bay into an infinite medium like the ocean
via the tidal volume (R is in units of I/tidal cycle)
R _ Tidal Volume (54^
Mean Volume of the Bay
Detailed discussion of the above parameters appears in several other
sources. ^'->
- 108 -
-------
Program Run Preparation;
A. Prior to punching any data cards for the program run of the
estuary that you are modelling, the following steps should
be kept in mind.*
1. Draw the symbolic river stretch (as in Figures 13 and 16
of the example problems).
2. Decide where to section your model to consider waste loads
(point, uniform, benthal, etc.)> physical parameters and
hydraulic parameters.
3. Number your sections (See Note 1, p. 118).
4. Determine your section types and define each junction as
to the type of boundary equation which is to be applied
there.
5. Determine your reference junction(s) (Note 2, p. 118) and
on your diagram place distances for each junction.
6. List all required parameters for section and junctions in
tabular form.
B. Punch the following input cards. (See Data Description.)
1. Identification Card (1 card/run).
2. Section Cards (A).
a. One for each section and in consecutive ascending
order.
*Reference will be made in this section to the example problems which
begin on Page 130.
- 109 -
-------
3. Section Cards (B).
a. One for each, section and in consecutive ascending
order.
4. Junction Card-Concentration.
a. One for each double junction, 2 for each triple
junction, and none for a dam.
5. Junction Card-Mass Balance.
a. One for each junction.
6. Repeat preceding 5 steps if you are using a model with
2 input components (i.e., carbonaceous and nitrogenous
load).
C. Place cards in following order to run program (also see
Figure 24).
1. Job Card.
2. Necessary JCL Cards.
3. Job ID Card.
- 110 -
-------
4. Section Card A for 1st Section
II II T> II l| ||
Section Card A for 2nd Section
Section Card A for Nth Section
II IT 13 II II II
5. Concentration Junction Card for Junction 1
Mass Balance Junction Card for Junction 1
Concentration Junction Card for Junction 2
Mass Balance Junction Card for Junction 2
Concentration Junction Card for Junction N
Mass Balance Junction Card for Junction N
6. /*
(NOTE: If there are more models or the run contains 2 input
components, cards 3-5 from above are placed before
the end of file (/*) card.
- Ill -
-------
Figure 24
Deck set up for run using a load module on EPA computer (OSI)
n total number of sections
m = total number of junctions
/*
JOB IDENTIFICATION CARD (SECOND RUN, IF ANY)
MASS BALANCE CARD FOR JUNCTION 1
CONCENTRATION EOUALITY CARD FOR JUNCTION I
MASS BALANCE -CARD FOR JUNCTION 1 = m-1
MASS BALANCE CARD FOR JUNCTION 1
CONCENTRATION EQUALITY CARD FOR JUNCTION i = 1
SECTION CARD (B) FOR SECTION I
SECTION CARD(A) FOR SECTION i
SECTION CARD (B) FOR SECTION i
SECTION CARD (A) FOR SECTION ! = 1
JOB IDENTIFICATION CARD (FIRST RUN)
DD
//FTftSFOOl
SV50UT=A
// EXEC PGM'MODEL,
DISP=SHR,DSN = EPA.R2.ES001
//JOBLIB DD UNIT=333n,VOL=SER=BCSni2,
/APRIORITY
//XXXESOOl JOB (A100,X02,8,5),NAME,MSOLEVF,L=1
112 -
-------
DATA DESCRIPTION & INPUT REQUIREMENTS (ALSO SEE DECK SET-UP)
Note: Superscripts refer to the notes
following this section.
CARD 1 (JOB IDENTIFICATION AND DATE CARD)
Column(s)
1-4
5
6-9
10
11-13
14-16
17-19
(12)
20-22
23-25
26-30
31-35
36-40
41-45
46-80
Math
Symbol
Variable
IDJOB
IDATE
IS2
IS1
NN
IND
ICON
FAC1
FAC2
FAC3
ITYPE
Description
Job name or no. (Arbitrary)
blank
Date of run (month, year)
blank
No. of double sections
No. of single sections
Indicator (0,1) which describes if:
0 - Current model is a new one.
1 - Current model is a repeat of the previous
one but the waste load magnitude or its
location have been changed.
Indicator (0,1) which causes the program to
print (1) or not to print (0) the DO deficit
and BOD matrices.
Indicator (0,1,2) for coupled system, single system,
and coupled with 2 inputs, respectively.
Temperature correction factor for AK.
Temperature correction factor for DK.
Temperature correction factor for RK.
blank
Description of system being modelled for
identification purposes.
Units
Format
A4
IX
A4
IX
13
13
13
13
13
F5.D
F5.D
F5.D
5X
35A1
-------
Card 2 (SECTION CARD A)
Column (s)
1
2
3-5 ,. ,
60 \-J- /
— o
9-10
11-13
/ 1 \
14-18 (2)
19-25
i 26-30J;:'
H 31-35 *• ;
£ 36-40
i 41-48
49-56
57"64(4)
65-72 (4)
73-76
77-80
Math
Symbol
A
RK
AK
F
Q
BD
WU
R
Variable
S
I
NS2
I
NJU(I)
XJU(I)
AREA(I)
RK(I)
AK(I)
FF(I)
Q(D
BD(I)
WU(I)
TIDEl(I)
IDJ
Description
Symbol to identify section card)
Used to help with deck set-up
Designation of Section Card 1
No. double sections (optional)
Section no. (in consecutive & ascending order)
blank
Arbitrary name assigned to upstream junction
of section
Upstream distance of junction NJU(I) from
reference junction
Cross-sectional area of section
BOD removal coefficient of section
Reaeration coefficient
Ultimate to 5-day BOD ratio
*Flow
*Bottom oxygen demand (Benthal)
*Uniform waste input
*Tidal exchange coefficient
blank
Job identification (must agree with
IDJOB-card 1)
Units
Mi.
2
Day'?"
Day
CFS
Gm/M2-Day
Lbs/ Day-Mi.
Format
IX
IX
13
13
A3
F5.0
F7.0
F5.0
F5.0
F5.0
F8.0
F8.0
F8.0
F8.0
A4
*0nly required when available or applicable, otherwise blank or 0.
-------
CARD 3 (SECTION CARD B)
Column (s)
1
2
3~5(1)
ft— R
9-10(5)
11-13
/ o \
, 14-18 (2)
M 19-25
, 26-30
31-35
36-40
41-48
49-56
57-64 ,..,.
65-72C '
73-76
77-80
Math
Symbol Variable
S
2
NS1
I
KFO(I)
NJD(I)
XJD(I)
Blank
DK DK(I)
E E(I)
T TEMP (I)
R-P PR (I)
HT HT(I)
VOL VOL (I)
PDEL(I)
ID (I)
Description
Symbol to identify section card
Used to help with deck set-up
Description of section card 2
No. single sections (optional)
Section no. (must agree with card A)
Code no. for type of section
1 - Double section
2 - Negative infinity section
3 - Positive infinity section
4 - Completely mixed volume (BAY)
Arbitrary name assigned to downstream junction
of section
Downstream distance of junction NJD(I)
from reference junction
BOD oxidation coefficient
Dispersion coefficient of section
Temperature of section
*Algal photosynthesis and respiration
Water depth
*Volume of section (BAYS)
Interval desired for output along section
blank
Section ID used to identify section in
output (arbitrary nomencl . )
Units Format
IX
IX
13
13
12
A3
Mi. F5.0
1 7X
Day" F5 . 0
Mi. 2 /Day F5.0
°C F5.0
Mg/l-Day F8.0
Ft. F8.0
Ft. F8.0
Mi . F8 . 0
A4
-------
CARD 4 (JUNCTION CARD - CONCENTRATION)
' I \
Math
Column(s) Symbol
1 (7)
2-3U;
4-6
fo\
7-9 ^ '
("9")
10-12 2
13-18UJ
/•q1)
1 Q_?1 v^'
c?')
22-27 V'
n rn
28-3og°>
31-36 v;
n "M
37-42 *• ' w
43-48 QD
49-54 DD
55-60 ID
61-76
77-80
Variable
J
I
JND(I)
JBV(I)
ISA(I)
XSA(I)
ISB(I)
XSB(I)
ISC(I)
XSC(I)
W(I)
QD(D
DD(I)
BDD(I)
IDJ
Description Units
Symbol to identify junction card (optional)
Junction card number
Junction identification (same as NJD & NJU
from section card)
Boundary type of J'th junction, concentration
'First1 section forming junction
Distance of 'first' sec. forming junction from its Mi.
own reference junction
'Second' section forming junction
Distance of 'second' section forming junction from Mi.
its own reference junction
*'Third' section forming junction
*Distance of 'third' section forming junction Mi.
from its own reference junction
(if triple junction)
*Point waste input Lbs./Day
*Flow over dam CFS
*Deficit at face of dam inflow QD MG/LIT.
* BOD in Flow OD MG/LIT.
blank
Job identification (same as card 1)
Format
IX
12
A3
13
13
F6.1
13
F6.1
13
F6.1
F6.0
F6.0
F6.0
F6.0
A4
*0nly required when available or applicable, otherwise blank or 0.
-------
CARD 5 (JUNCTION CARD - MASS BALANCE)
Math
Column(s) Symbol Variable Description Units Format
SAME AS CARD 4 EXCEPT:
2-3 ,„, I Junction card number 12
7-9 * . JBV(I) Boundary type of J'th junction, mass balance 13
ISC (I) *' Third' section forming junction 13
31-36 ' XSC(I) ^Distance of 'third' section of junction Mi. F6.1
from its own reference junction
(if triple junction)
"Only required when available or applicable, otherv/ise blank or 0.
-------
Notes on Input Cards
1. (Col. 6-8)
(a) The consecutive order must be observed because any listing
(by section) which is made by the computer program is in the
form 1,2,3, ...n(n= total number of sections) and it
overrides your numbering system unless it is in the above
form. For instance:
Your Input
(If sections are num-
bered as follows)
1 2
2 or 4
Program Listing
(Resulting Section
Numbering)
1
2
(b) In numbering your sections for a particular model the double
sections must have a lower number than the single. If this
is not observed the program prints out NCL (Non-Classifiable)
and the program will not execute. For instance it must be:
single
5
3
2
1
4
single
double•
2. (Col. 14-18)
The reference junction is a junction that can be arbitrarily chosen
to be the origin of your model and from which other sections can be
located by a distance (miles). This distance can be positive or
negative from the reference junction (where downstream is positive
and upstream is negative).
(a) The reference junction can be chosen at any point on your model.
(b) There can be as many reference junctions as you have sections,
however, when you choose such a junction you must decide which
sections it is applicable to and any reference to these sections
must be in terms of this reference junction.
- 118 -
-------
(c) To simplify matters, however, choose as few reference junctions
as possible and place them at strategic points.
(d) When deciding on your reference junction distance it is better
to start with 0 mileage since this mileage is placed in an
arithmetic expression of the form ek* and if a combination of
k and x becomes very large the solution might blow up.
(e) Junction AAA in example problem 1 and junction CCC in example
problem 2 are reference junctions.
3. (Col. 26-36) and (Col. 31-35)
These two variables must not be equal since their difference appears
in a denominator in the program. A warning is printed in regard to
this and corrective action must be taken.
4. (Col. 65-72)
Only applicable to a completely mixed section where flow goes into
the bay and is defined as:
TIDAL VOLUME * # of tidal cycles
mean volume of bay. day
Otherwise, equal to zero.
5. (Col. 9-10)
The negative infinity section is usually denoted by the uppermost
upstream section and the positive infinity section as the farthest
downstream section.
6. (Col. 65-72)
If PDEL is left blank the program sets its value to 1 mile.
7. (Col. 2-3)
Junction card numbers are in the form of 1, 2, 3, . . . n, (a) where
n = total number of equations which is equal to: 2* number of double
sections + number of single sections; (b) the junction card numbering
must begin with 1 and proceed in consecutive ascending order, either
upstream or downstream, to n. The order must be, a concentration
junction card followed by a mass balance junction card or vice versa.
This procedure is repeated as you move from junction to junction,
with the number punched on each succeeding card increased by 1 until
the n'th card.
- 119 -
-------
The boundary type is a code number indicating the particular equation(s)
to be applied at the junctions of interest depending on the types of
sections contiguous to that junction. These boundary types and their
numbers are as follows:
Concentration Equations
1 - Concentration equation for 2 double sections.
2 - Concentration equation for 1 double section and 1 negative in-
finity section.
3 - Concentration equation for 1 double section and 1 positive
infinity section.
4 - Concentration equation for 1 double section and a completely
mixed bay.
Mass Balance Equations
5 - Mass Balance equation for 2 double sections.
6 - Mass Balance equation for 1 double section and 1 negative
infinity section.
7 - Mass Balance equation for 1 double section and 1 positive
infinity section.
8 - Mass Balance equation for 3 double sections where flow from
one divides into the other two.
9 - Mass Balance for 3 double sections where flow from 2 com-
bines into the third.
10 - Mass Balance a section of advective flux only (dam).
11 - Mass Balance at a bay where flow leaves through the double
section.
12 - Mass Balance at a bay where flow leaves through 2 double
sections.
13 - Mass Balance at a bay where flow enters the bay from a double
section.
14 - Mass Balance at a bay where flow enters the bay from 2 double
sections.
- 120 -
-------
9. (Col. 10-12) and (Col. 19-21)
In the variable description it states, first and second section form-
ing a junction. This statement requires one to know not only the
actual section numbers (which are read from your model diagram and
are punched on the card) that form the junction but which of the two
is the "first" (Col. 10-12) and which is the "second" (Col. 19-21).
This first and second designation depends on section types (double,
negative infinity, etc.) that meet at the junctions. Listed below
are the 14 permissible section combinations and their card column
positions for any one junction card. The list is in terms of
symbols which are defined below:
D = double section, - °° = neg. inf., + °° = pos. inf., 0 = bay
Boundary
Type
(JBV)
1
2
3
4
First
Section
(ISA)
CD
— oo
cm
0
Second Third
Section Section
(ISB) (ISC)
a
cm
+ 00
0
5 Cn CD
6
7
8
9
10*
11
12
13
— oo
C=D
cm
cm
CD
CD
cm
cm
cm
-I- oo —
CD \ — \
r— \ CT3
a
0
cm o
0
14 C^3
*Both section numbers are the same.
- 121 -
-------
(NOTE: At a particular junction in any of the above cases where you have
only 2 double sections, the upstream-downstream or downstream-upstream
convention must be followed for the entire system being modeled.
In other words, if you are numbering from upstream to downstream,
at a particular junction of 2 double sections, the upstream section
is the "first section" and the downstream section the "second sec-
tion.").
10. (Col. 28-30)
When 3 double sections form a junction as in boundary conditions
8 and 9, the following placement of the sections must be adhered
to on the Junction Cards.
(a) Flow from 1 Section into 2 Sections
For Junction JT:
2 Concentration
Equality Cards
ISA
3
3
ISB
1
2
ISC
1 Mass Balance Card
x*
NOTE: Section 3 must be in the ISA location
(Cols. 10-12).
*The x stands" for either of the two other sections,
- 122 -
-------
(b) Flow from 2 Sections into 1 Section
For Junction JT:
2 Concentration
Equality Cards
1 Mass Balance Card
ISA
1
2
X*
ISB
3
3
3
ISC
-
X*
NOTE: Section 3 must be in the ISB location
(Cols. 19-21).
11. (Col. 37-42)
Not allowed at a dam junction.
12. (Col. 17-19)
This value cannot be 1 if a coupled system with 2 input components
is being modelled. The program will terminate if this is done.
*The x stands for either of the two other sections,
- 123 -
-------
OUTPUT
1. INPUT DATA LISTING
Th« 80 column card images of the input data are read in and th«n
printed in a format that spreads out the variables for easier read-
ing.
a. The Job Identification and Date Card is listed first.
b. The section cards are listed next. As this is done, the word
ERROR followed by the variable is printed when AK = RK* and/or
HT = 0*. (See Input Restrictions, Notes B(l) and B(2)).
c. Immediately following the above, the sections are listed and
identified by type according to the f ollox^ing:
(1) DBLE - A double section.
(2) S PI - A single positive infinity section.
(3) S NI - A single negative infinity section.
(4) S CM - A single completely mixed section (BAY).
In addition two error conditions are identified by MISS and NCL.
Their meanings are as follows:
*(1) MISS - Type of section (KFO) is not punched on the section
card.
*(2) N CL - Non-Classifiable. Either a section is not classified
correctly according to (KFO) or a double or single
section is not correctly classified according to
Note 1 (p.118).
d. The junction cards are then listed and parentheses are added to
more clearly identify the sections that meet at that junction
and the mile point of that junction.
e. Immediately following the previous listing, a line is listed for
each junction card. The following information is printed; equation
number (which is identified in the input as junction card number),
boundary type (JBV), abbreviation of the boundary type, and in
parentheses the section numbers and type of sections (see (c) above)
which meet at that junction.
*This condition will cause the run being done to be terminated.
- 124 -
-------
The abbreviations of the 14 boundary types are as follows:
(1) C 2D - Concentration equation for 2 double sections
(2) C DN - Concentration equation for 1 double section
and 1 negative infinity section
(3) C DP - Concentration equation for 1 double section
and 1 positive infinity section
(4) C DC - Concentration equation for 1 double section
and a completely mixed bay
(5) M 2D - Mass Balance equation for 2 double sections
(6) M DN - Mass Balance equation for 1 double section
and 1 negative infinity section
(7) M DP - Mass Balance equation for 1 double section
and 1 positive infinity section
(8) M3D1 - Mass Balance equation for 3 double sections
where flow from one divides into the other two
(9) M3D2 - Mass Balance for 3 double sections where flow
from 2 combines into the third
(10) M DM - Mass Balance a section of advective flux only (dam)
(11) MCM1 - Mass Balance at a bay where flow leaves through
the double section
(12) MCM2 - Mass Balance at a bay where flow leaves through
2 double sections
(13) IIICM - Mass Balance at a bay where flow enters the bay
from a double section
(14) M2CM - Mass Balance at a bay where flow enters the bay
from 2 double sections
In addition at the end of each line the word PASS or FAIL* is printed,
This denotes the correctness (PASS) or incorrectness (FAIL) of the
punching of the section numbers on the junction cards according to
the criteria of Note 9 (p. 121). The only exception being the case
of only double sections appearing on a junction card. When this
occurs the upstream-downstream rule of Note 9 must be followed and/or
the special instructions for triple junctions (Note 10, p. '^-22 ).
*This condition will cause the run to be terminated.
- 125 -
-------
2. TABLE OF SECTION AND JUNCTION PARAMETERS
a. The input data in the original units is listed in a more struc-
tured manner and with headings. The headings include the vari-
able identification according to the Input Data Description.,
abbreviations of the variable names, and their units. The param-
eters pertaining to sections are listed first. This listing
includes the section number, the section identification (ID),
the type of section from l(b) above, KFO, and the upstream and
downstream junction identifications (NJU and NJD) with the
actual mile point for each. This is followed by the data for
each section (e.g., flow, area, etc.).
b. Immediately following are the junction card variables. Included
on this Table are the equation numbers, the boundary type and
the abbreviation for the boundary, the junction identification (JND),
and in parentheses, the section numbers and the mile point of that
specific junction, followed by the data.
3. MATRIX PRINT-OUTS
a. A symbolic matrix that is formed from the simultaneous linear
equations is listed first. This matrix represents the positions
of both BOD and D.O. deficit elements according to the inputting
of the data cards and the segmentation of the estuary being
modelled. The elements are in terms of the various junction
identifications (JND's). There is also a row heading to denote
the equation numbers and the boundary type (JBV). Each column
corresponds to a particular coefficient of the solution.
*b. The Input BOD matrix is then printed. The elements of this are
calculated from the input data and the BOD equations. The last
row represents the elements in the BOD forcing function matrix.
*c. The inverted BOD matrix is then printed. At this point the units
of the various elements are in DAY/MI3. The last row now contains
the solutions for the unknown BOD coefficients. (B and C for each
section.)
*d. Immediately following the above the D.O. Deficit matrix is listed.
The elements of this are calculated from the input data and the
Deficit equations. The last row represents the elements in the
Deficit force function matrix.
*e. Finally the inverted D.O. Deficit matrix is listed. As in the
above the units of the various elements are in DAY/MI3. The
last row now contains the solutions for the unknown Deficit
coefficients. (G and H for each section.)
*The printing of these matrices is optional and is controlled by IND of
input card 1 (Job Identification and Date Card).
- 126 -
-------
4. TABLE OF SECTION AND JUNCTION PARAMETERS IN COMPATIBLE UNITS
a. The next portion of the print-out lists the same output as in
(2) above (structured form of the input data), but certain
conversions have been applied to some of the variables to make
the units compatible (see Subroutine DATA2).
b. In addition certain section variables, which have been computed
in the program, are also listed. The additional variables and
their definitions are as follows:
(1) U - The flow velocity in MI/DAY.
(2) S and V - The coefficients in the exponentials of the
BOD solution.
(3) SD and VD - The coefficient in the exponential of the
D.O. deficit solution.
(4) ZZ and FAC - Parts of the D.O. deficit forcing function.
BD . PR . DK*WU*FF
HT*AK AK AK*RK*A
_ DK*FF
FAC " AK-RK
(5) The unknown section coefficients that have been solved for
by the matrix inversion techniques.
(a) B and C - Values of the BOD coefficients in mg/1;
(b) G and H - Values of the D.O. deficit coefficients in mg/1.
*5. BOD AND D.O. DEFICIT SOLUTIONS
The resultant BOD and D.O. deficit values in mg/1 are listed by
ascending order of section and distance within sections according
to the printing interval (PDEL) designated in the input cards. In
addition, the section identification (ID) is printed on the left-
hand side of the listing and the junction identification (JND) at
both ends of the section is printed.
NOTE: If any negative BOD or D.O. deficit values are outputted or these
values are not equal at section junctions, check the program
restrictions and/or the input description of the variables.
*For a model that contains 2 waste input components such as carbonaceous
and nitrogenous wastes, the data for each component will be printed along
with its matrices and solutions as separate entities. The sum of their
responses (as outlined in 5) will then follow.
- 127 -
-------
RESTRICTIONS
NOTE: THE NUMBERS IN PARENTHESIS REFER TO THE NOTES ON DATA INPUT.
A. GENERAL
1. At the present there is a limit to the size of the estuary that
can be modelled.
a. In terms of junctions it means the maximum number of equa-
tions is 100.
b. In terms of sections the following expression gives the
limit. (2 x number of double sections + number of single
sections) < 100.
2. As with any model, ES001 has limits. When they are exceeded,
ES001 will either reject the problem or compute erroneous
values (e.g., negative values of BOD). The following guide-
lines are delineated to assist the programmer in formulating
a realistic, executable model.
If certain system parameters are exceeded the solution may blow
up (approach infinite values or negative values). Fortunately,
such extreme values are rarely encountered in actual estuaries.
The following table gives suggested ranges of parameters which
are within the limits of ES001.
Parameter Acceptable Values
E 1-20 mi2/day
Reaction rates < I/day
U < 1 fps
Section length < 20 miles
Table 7 Suggested parameter limitations for ES001.
These are by no means rigid constraints, as each particular config-
uration will have its exact limits, but care should be exercised when
approach these extremities.
For instance, as the upper end of the estuary is approached, the
dispersion coefficient (E) will decrease due to the absence of salinity
gradients which are present in the downstream end. If E gets small
enough it will cause the computation to blow up. This is due to the
fact that the elements of the matrices are primarily of the form:
- 128 -
-------
fi- [1 ± tl
u
where
X = mile point
U = velocity
K = reaction coefficient (AK in DO deficit solution or RK in
BOD solution).
When E gets too small, the exponent of the above expression gets very
large and in turn the entire expression gets very large. This will
often cause the solution to blow up.
In such a «.ase it would be advisable to shorten the section lengths
and minimize t'ie value of X by proper positioning of the reference junc-
tion. When a blowup occurs, the printout of the original matrix can
often be useful in determining where the cause of the error is by exam-
ining each of the elements to see which are inordinately large.
B. SECTIONS
1. For any one section the reaeration (AK) coefficient cannot equal
the removal (RK) coefficient (3) . *
2. For any section that has a benthal waste source the depth (HT)
cannot be zero (6).
3. For any run, the sections must be numbered consecutively with
the double sections having lower numbers than the single (1).
C. JUNCTIONS
1. When punching the junction cards, the order must be upstream-down
or vice versa (7) .
2. The column location of the section numbers on any junction card
has a specific order (9).
3. The program cannot handle anything greater than a triple junc-
tion and there are specific rules for a triple junction (10) .
* Numbers in parenthesis refer to notes beginning on page 118,
- 129 -
-------
EXAMPLE RUNS
A. The following two example RUNS will aid one in formulating a model
to run in a manner that conforms with the computer program, and thus
will result in the fewest errors when the first run is made.
To aid anyone with only slight familiarity with water quality models
the following definitions are made:
a. Double Section (i.e., Sec. 1, Fig.25 ) - A section that is
bounded on both ends by another section.
b. Single Sections - End sections which are only bounded on one
end by another section. The two types are:
(1) Infinity Sections - End sections where the BOD and DO
deficits approach zero as the length of the section
approaches infinity. When actually setting up a model
this is usually limited to a few miles since this is
the boundary of the system under consideration.
*(a) Negative Infinity (i.e., Sec. 2, Fig. 25) - An infin-
ity section that is by convention the furthest upstream
section.
*(b) Positive Infinity (i.e., Sec. 3, Fig. 25) - An infin-
ity section that is by convention the furthest down-
stream section.
(2) Completely Mixed Volume (i.e., Sec. 7, 8, Fig. 28) - An
example of this is a bay.
c Junctions
(1) Actual Junction - A river mile point where two or more
sections meet or at which a dam is located.
(a) Single Junction - The point in a double section at
which a dam is placed. (Fig.28 .)
(b) Double Junction (i.e., JUNG BBB, Fig. 28) - A junction
where two sections meet.
(c) Triple Junction (i.e.. JUNG CCC. Fig.28 ) - A junction
where three sections meet.
*There can be more than one if the estuarine system you are modelling
has tributaries entering it. The upstream and downstream on your dia-
gram is shown by the direction of the flow (Q).
- 130 -
-------
(2) End Section Junction (i.e., JUNG NBC, Fig. 25 ) - The river
mile points that define the limits of your model.
Now the following steps should be followed to prepare the 2 sample
problems :
1. Draw the river stretch in a symbolic manner as in Figure 25 for
the simple configuration, and Figure 28 for the more complex
problem. (Figure 28 is derived from Figure 26-)
2. Section your model to consider
a. Location of point, uniform and benthal sources and from
Figure 28 the different types are:
(1) Point Inputs - 5000 Ibs./day at junction DDD and
1000 Ibs./day in bay (Sec. 7) but since point waste
loads can only be inputted via the junction cards,
the load into the bay must be placed on the junction
cards for junction GGG.
(2) Uniform Waste Load - 500 Ibs./day /mi. in Section 1.
b. Physical parameters (i.e., cross-sectional areas).
c. Hydraulic parameters (i.e., flow, dispersion).
d. Other program characteristics as described in Restrictions
(p . 028) .
3. Number your sections so that they are consecutive and the double
sections have lower section numbers than the single.
4. Determine your reference junction (see page ng) (junction
CCC on Fig. 28,) and place routine mileages on your diagram from
this junction (nil6 5 at junction DDD.)
5. Finally, give each junction an I.D. and determine the section
types which are required on the input cards.
- 131 -
-------
EXAMPLE RUN #1 - SIMPLE ESTUARY
A. INTRODUCTION
PBC
BBB
AAA
NBC
3
(00)
(-00)
t
Flow
point load
20
10 0
DISTANCE (MILES)
-10
Figure 25 . A Simple Estuary.
Model the above estuarine system for dissolved oxygen and BOD
using the ES001 program.
- 132 -
-------
JOB ID AND DATE CARD
SIMP 0672 1 2
INPUT SECTION CARD LISTING
0 1.O2O 1.04O 1.040 EXAMPLE PROBLEM 1
1 1 1 AAA 10.000 18800.0 0.300 0.26* 1.250 375.0 0.0 0«0
211 BBS 20.000 13800.0 0.300 9.400 25.000 0.0 5.250 0.0
121 NBC 0.0 18800.0 0.300 0.264 1-250 375.0 0.0 0.0
222 AAA 10.000 13300.0
0.300
9.400 25.000
0.0 5.250 0.0
131 B3B 20.000 13800.0 0.300 0.264 1.250 375.0 0.0 0.0
233 PBC 30.000 10,0 0.300 9.400 25.000 0.0 5.2JO 0.0
SECTION 1 IS
DBLE
SECJION 2 IS S Nl
SECTION 3 IS S PI
INPUT JUNCTION CARD LISTING
0.0
0.0
0.0
0.0
0.0
0.0
SIMP
ooua
SIMP
UXST
SIMP
DNST
1 AAA
2 AAA
3 BBB
4 BBB
EQUATION
EQUATION
EQUATION
EQUATION
2 2 10.00
6 2 10.00
3 1 20.00
7 1 20.00
1 OF TYPE 2 C OH
2 OF TYPE 6 M Df;
3 OF TYPE 3 C OP
4 OF TYPE 7 M OP
1 10.00
1 10«00
3 20.00
3 20.00
( 2.S NI 1
< 2,S HI )
1 l.OBLE )
! 1>OBLE I
0 0.0
0 0.0
0 0.0
0 0.0
( l.DBLE
t l.DBLC
( 3.S PI
I 3.S PI
190000.0 0.0
100000.0 0.0
0.0 0.0
0.0 0.0
) i Of
) ( Oi
1 (Ot
) ! 0.
0.0 0.0
0.0 0.0
0.0 0.0
0*0 0.0
) PASS
) PASS
I PASS
) PASS
SIMP
SIMP
SIMP
SIMP
-------
MODFL RUN SIMP
COMPUTED 0672
EXAMPLE PROBLEM 1
DATA GIVEN IN ORIGINAL
UNITS
SECTION
K.FO JUNCTION ID AND DISTANCE
AT ENDS
REM.RATE REA.RATE
,1/OAYI
-LT/SDAT
TEMP
1
2
3
1
2
3
1
2
3
4
DOU8 DBLE
UXST S NI
DNST S PI
SECTION
DOUB DBLE
UXST S NI
DNST S PI
EQUATION
2 C ON AAA
6 M DN AAA
3 C DP BBB
7 M DP BBB
• -.til i i
-------
ROW JBV 1234
1 2 AAA AAA AAA
2 6 AAA AAA AAA
3 3 BBS BBB - BBB
4 7 BBB BBB - BBB
SYMBOLIC DO DEFICIT AND BOD MATRIX
-------
MODEL RUN SIMP
COMPUTED 0672
EXAMPLE PROBLEM 1
DATA GIVEN IN CCXPATIBLE UNITS
SECTION
KFO JUNCTION ID AND DISTANCE
AT DlDS
RK AH 0< E FF TEMP PDEL
RCM.RATE HCA.riATC PCC.nATC DIEP.RATE UUT/EfAY MILC~If!CR
(1/CAVI (1/DAYI (I/DAY) (MI2/DAYI (CENT) (MI)
1
2
3
DOUB
UXST
DNST
DBLE
S NI
S PI
SECTION
1
2
3
DOUB
UXST
DNST
DBLE
S NI
S PI
SECTION
1
2
3
DOUB
UXST
DNST
DBLE
S NI
S PI
EQUATION
1
2
3
*
2 C
6 M
3 C
7 M
DN AAA
DN AAA
DP BBB
DP BBB
1 (BBB.
2 (AAA,
3 (PBC»
CROSS
SEC-AREA
(MI»*2 )
0.674E-03
0.674E-03
0.674E-03
U
0.32640 0
0.32640 0
0.32640 0
20.0) (AAA, 10.0)
10.0) (NBC, 0.0)
30. C) (BBB, 20.0)
0 BD
FLOW BENTHL-DMD
(MI»*3/DAY) (GMS/M2-DAY)
0.220E-03 0.0
0.220E-03 0.0
0.220E-03 0.0
S V SO
.21518 -0.18045 0.19431
.21518 -0.13045 0.19431
,21518 -0.18045 0.19431
JUNCTION ID, SECTIONS
MEETING AT JUNCTION.
AND DISTANCE (MI
I 2 10.0)1
( 2 10.0) (
I 1 20.0)1
I 1 20.0)1
1 10.0)1 0 0.0)
1 10.0)1 0 0.0)
3 20.0) 1 0 0.0)
3 20.0)1 0 0.0)
0.365 0.291 0.365
0.365 C.291 0.365
0.365 0.291 0.365
PR TICE1
PHOTO-RESP TIDE-COEF
(MG/L-DAY)
0.0 0.0
0«0 0.0
0.0 0.0
VD B C
-0.15958 -0.00000 26.38181
-0.15958 0.50473 0.0
-0.15958 0.0 26.38179
W 00
PT-LOAD FLOW-DAX
••3/DAY) IMG/L) t.VI»»3/DAV)
0.01089 0.0
0.01089 0-0
0.3 0.0
0.0 0.0
9.400
9.400
9.400
HT
DEPTH
(MTRS)
0.160E 01
0.160E 01
0.160E 01
G
0.00000
4.31480
0.0
DO
DEF-CAM
(MO/LI
0.0
0.0
0.0
0.0
1.250
1.250
1.250
VOL
VOLUME
0.0
0.0
0.0
H
148.55461
0.0
148.55455
BDD
BOD-DAM
IMG/LI
0.0
0.0
0.0
0.0
25.000 0.0
25.000 0.0
25.000 0.0
WU/AREA»RK
UMF-LOAD
(MG/L)
0.0
0.0
0.0
ZZ FAC
0.0 -6.21
0.0 -6.21
0.0 -6«21
-------
SIMP
EXAMPLE PPOBLEM I
COMPUTED 0672
SECT.
NO.
1
1
1
1
1
1
1
1
1
1
1
SECT.
NAME
DOUB
OOUB
DOUB
DOUB
DOU8
DOUB
DOUB
DOUB
DOUB
DOUB
DOUB
D1ST.
MI
10.0
11.0
12.0
13.0
U.O
15.0
16.0
17.0
18.0
19.0
20.0
BOD
MG/L
*.3M2
3. 62*4
3.0260
2. 526*
2.1092
1.7610
1.4702
1.2275
1.02*8
0.8556
0.711.3
D.O DEFICIT
MG/L
?.1775
3.1632
3,1097
2.9816
2.8179
2.6327
2.1.368
2.2380
2.0*20
1.8528
1.6730
JUNCTION
ID
AAA
BBB
-------
SIMP
EXAMPLE PROBLEM J
COMPUTED 0672
SECT.
NO.
2
2
2
2,
2
2
2
2
2
2
2
SECT.
NAME
UXST
UXST
UXST
UXST
UXST
UXST
UXST
UXST
UXST
UXST
UXST
OIST.
MI
0.0
1.0
2.0
3.0
4 . 0
5.0
6.0
7f\
• 0
8.0
9.0
10.0
BOD
MG/L
0.5048
0.6260
0.7763
0.9626
1. 1937
1 .4803
1.8357
2.2764
2.8230
3.5007
A •* A I 5
** . JM 1 i
D.O DEFICIT ' JUNCTION
MG/L ID
— — — •-*— •
1-1822 NBC
1.3555
1.5'.67
1.7551
1.9785
2.2131
2.4525
2.6866
2.9010
3.0745
3.1775 444
-------
SIMP
EXAMPLE PROBLEM i
COMPUTED 0672
SCCT.
NO.
3
3
3
3
3
3
3
3
3
3
3
SECT.
NAME
DNST
DNST
DNST
DNST
DNST
DNST
DNST
DNST
DNST
DNST
DNST
DIST.
MI
20.0
21.0
22.0
23.0
24.0
25.0
26.0
27.0
28.0
29.0
30.0
OOP
MG/L
0.7143
0.5964
0.4979
0.4J57
0.3471
0.?898
0.2419
0.2020
0. 1686
0. 1409
0.1175
D.O DEFICIT JUNCTION
MG/L ,D
1.6730 BBB
1.5043
1 .3476
1.2032
1.0712
0.9511
0.8425
0.7446
0.6569
0.5784
0.5085 PRr
-------
EXAMPLE RUN NUMBER 2
A. INTRODUCTION
Figure26 . ECHO CREEK.
Echo Creek is a fictitious polluted estuary in the Virgin Islands
(see Fig. 26. Presently, two industries, United Mango, U.S. Pistachio,
and one village, Gaffe, are located within the estuarine system. United
Mango discharges its waste to Punch Bay which flows into Echo Creek via
a one mile stretch called Punch Run. U.S. Pistachio discharges directly
into the Echo at a point 5 miles downstream from its confluence with
Punch Run. Gaffe, which is situated upstream from the confluence has
no sanitary facilities and its waste eventually reaches the estuary in
an untreated form. Gaffe's waste in contrast to the industries' contains
an appreciable amount of settleable solids which have built up on the
streambed as a sludge deposit. This deposit exerts a benthal demand on
the dissolved oxygen concentration along the entire length of the estuary
adjacent to the town. Finally, a dam is located at Gaffe's upstream
boundary and a bay which is joined to the ocean is located at the down-
stream end of the Echo. The flow coming over the dam is slightly
polluted.
- 140 -
-------
A model is to be developed to determine the amount of waste treat-
ment required by the municipality and the industries to improve water
quality in the system to an acceptable level. Use of program ES001
will be employed to formulate such a model for BOD and DO deficit.
B. SEGMENTATION
Considering the confluence of Echo Creek and Punch Run as mile
point zero, the following values for the cross-sectional area of the
estuary are shown in Figure 27:
1800 -
1600 -
1400 -
I
1200 -
1000 -
800 -
600 -
400 -
200 -
if i i i »
15 10 5 0 -5 -10
Figure 27. Plot of Cross-sectional Area Versus
Distance Along Echo Creek.
- 141 -
-------
With this in mind the following segments and junctions were determined;
SEA*
EEE DDD CCC BBB DAM
FFF
GGG
BAY*
Figure 28. Segmentation of ECHO CREEK.
*NOTE: SEA and BAY are the names of pseudojunctions. They are meaning-
less but must be included when modelling bays.
- 142 -
-------
JOB ID AND DATE CARD
ECHO 0672 6 2
INPUT SECTION CARD LISTING
0 0
0 1.020 1.040 1.040
EXAMPLE PROBLEM 2
6 3 1 DAM -10.000 1000.0
2 l 1 BBB -5.000 1000.0
621 BBB -5.000 1000.0
2 2 1 CCC 0.0 1000.0
6 3 1 CCC 0.0 1000.0
231 DDD 5.000 1000.0
6 * 1 DDD 5.000 1000.0
2 * 1 EEE 7.500 1000.0
6 5 1 EEE 7.500 1500.0
2 5 1 FFF 10.000 1500.0
6 & 1 GGG -1.000 100.0
261 CCC 0.0 100.0
6 7 1 BAY -10.000 100.0
274 GGG -1.000 100.0
6 8 1 FFF 10.000 100.0
284 SEA 15.000 100.0
SECTION 1 is DBLE
SECTION 2 IS DBLE
SECTION 3 IS DBLE
0.275 0.100 i.2so 100.0 2.000 0.100E 03 0.0 ECHO
0.250 5.000 28.000 0.0 25.000 0.0 0.500 QAFF
0.250 0.100 ,.250 100.0 0.0 0.0 0.0 ECHO
0.250 5.000 28i000 0.Q 25iOCQ o%Q i>ooQ E^p
0 • 2 50 0-1QQ i-ien ,»«n
I<25° 11C-° 0-0 0.0 o.O ECHO
0.250 10.000 2a.OOO e.O 25.000 0.0 LOOO ECJC
0.250 0.100 1.750 iin n n r.
i./f5U 110.0 0.0 0.0 0.0 ECHO
0-250 !0.000 28,300 o.O 25.000 0.0 Il000 £CDN
0.250 0.100 1.250 no.o o.O 0.0 o.O ECHO
0.250 15.000 23.000 0.0 ,0.000 0.0 0.0 USPS
0.250 0.200 !.250 10.0 0.0 0.0 0.0 ECHO
0.250 5.000 30.UOO o.O 10.000 0.0 0.250 PUNC
0.250 0.050 ,.250 100.0 „.„ „.„ ^ ^
0-250 5.000 30.000 0.0 40.000 0.400E 06 1.000 PBAY
0.250 0.100 ,.250 110.0 0.0 0.0 0.500 ECHQ
0-250 5.000 28.000 0.0 40.000 0.120E 09 5.000 EC8Y
-------
SECTION 4 IS DBLE
SECTION 5 IS DBLE
SECTION 6 IS DBLE
SECTION 7 IS S CM
SECTION 8 IS S CM
INPUT JUNCTION CARD LISTING
1 DAM 10
2 BBB 1
3 BBB 5
4 CCC 1
5 CCC 1
6 CCC 9
7 ODD I
8 DDD 5
9 EEE 1
10 EEE 5
11 FFF 4
12 FFF 13
13 GGG 4
14 GGG 11
1 -10.00
1 -5.00
1 -5.00
2 0.0
6 0.0
2 0.0
3 5.00
3 5.00
4 7.50
4 7.50
5 10.00
5 10.00
6 -1.00
6 -1.00
1 •1°'00 ° O'O 0.0 100.0 0.2 0.2
2 •5'°° ° °'° o.o o.o o.o o.o
2 •5-°° ° o-o o.o o.o o.o o.o
3 °'° ° o.o o.o o.o o.o o.o
3 °*° ° o-o o.o o.o o.o o.o
3 °*° 6 °'° 0.0 0.0 0.0 0.0
4 5i°° ° o.o 5000.0 o.o o.o o.o
4 5'°° o o.o 5000.0 o.o o.o o.o
5 7'50 ° o.o o.o o.o o.o o.o
5 7>5° ° o.o o.o o.o o.o o.o
8 10*°° o o.o o'.o o.o o.o o.o
8 10'°° o o.o o.o o.o o.o o.o
7 -5.00 0 0.0 1000.0 0.0 0.0 0.0
7 -5.00 0 0.0 1000.U 0.0 0.0 O.o
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
ECHO
-------
EQUATION 1 OF TYPE 1O
EQUATION 2 OF TYPE 1
EQUATION 3 OF TYPE 5
EQUATION 4 OF TYPE I
EQUAT ION
EQUATION
EQUATION
EQUATION
EQUATION
EQUATION
EQUATION
EQUATION
EQUATION
EQUATION
5 OF TYPE 1
6 OF TYPE 9
7 OF TYPE 1
8 OF TYPE 5
9 OF TYPE 1
10 OF TYPE 5
11 OF TYPE 4
12 OF TYPE 13
13 OF TYPE 4
14 OF TYPE 11
M DM
C ZD
M 20
C 20
C 2D
M3D2
C 20
M 2D
C 2D
M 20
C DC
Ml CM
C DC
MCM1
( l.ODLE )'
( ItDBLE >
( l.DBLE I
< 2.DBLE )
( 6.DBLE )
< 2.08LE )
1 3.DBLE )
( 3.DRLE )
I 4, POLE 1
1 4.DBLE )
( StDBLE )
( 5.08LE )
( 6.DBLE )
( 6.D8LE )
( 1.D8LE I
( 2.08LE J
( 2.DBLE )
( 3.DBLE )
1 3.D5LE I
1 3.D3LE )
( 4.0BLE 1
( 4.DBLE )
( 5.DBL£ I
( 5.CB1.E )
( 8»S CM )
( a , s c M i
( 7iS CM )
( 7iS CV )
I O. 1 PASS
t Ot ) PASS
1 0. I PASS
I 0> I PASS
( 0 1 ) PASS
( 6iDBLE ) PASS
( 0 • J PASS
( 0> ) PASS
) PASS
I 0» ) PASS
1 0> ! PASS
I 0» ) PASS
( Ot ) PASS
-------
MODEL RUN ECHO
COMPUTED 0672
EXAMPLE PROBLEM 2
DATA GIVEN in ORIGINAL UNITS
SECTION
KFO JUNCTION ID AND DISTANCE
AT ENDS
AK
1 GAFF DBLE
2 CCUF> DBLE
3 ECJC DBLF
* ECDN DBLE
5 USPS DBLE
6 PUNC DBLE
7 PBAY S CM
8 ECBY S CM
SECTION
1 GAFF DBLE
2 ECUP DBLE
3 ECJC D8LE
* ECDN DBLE
5 USPS DBLE
6 PUNC DBLE
7 PBAY S CM
0 ECBY S CM
EQUATION
1 10 M DM DAM
2 1 C 2D BBB
3 5 M 20 BPB
* 1 C 2D CCC
1 (B3B, -5.0) (DAM, -10.0)
1 (CCC, 0,0) (ODD. -5.0)
1 (ODD,
i IEEE,
1 (FFF,
1 (CCC,
* (GGG,
4 (SEA,
CROSS
SEC-AREA
(FT«*2)
0.100E 04
0.100E 04
0.100E 04
0.100E 04
0.150E 04
0.100E 03
0.100E 03
0.100E 03
5>0) (CCC. 0.0)
7.5) (ODD, 5.0)
10.0) (EEE, 7.5)
0.0) (GGG, -1.0)
-1.0) (BAY, -10.0)
15.0) (FFF, 10.0)
0 eir>
U BD
FLOW BENTHL-DMD
(CFS) (GMS/M2-DAY)
0.100E 03 0.200E 01
0.100E 03 o.O
O.llOE 03 o.O
0.110E 03 0.0
O.llOE 03 0.0
0.100E 02 0.0
0.100E 03 o.O
O.llOE 03 o.O
0-275 0.100 0.250
o.;so c.-po c.-^n
v • A v. ^ ^ t t. 2 U
0.250 0.100 0.2TO
0-250 C.100 0.250
0.250 0.100 0.250
0.250 0.200 0.250
O.?50 0.050 0.250
0.250 0.100 0.250
PR TIDE1
PHOTO-RESP TIDE-COEF
(MG/L-OAYI
0.0 n n
U . 0
0.0 A 1-1
v • u 0.0
M f j
U" U U.O
O.O n r\
0.0
O.O rt r\
* • v U . 0
On n rt
u 0.0
fl M P rt rt
u« t. o.O
°*0 0.500E 00
5.000
5.000
10.000
10.000
15.000
5-000
5.000
5.000
HT
DEPTH
I FT. )
0.250E 02
0.250E 02
0.250E 02
0.250E 02
0.300E 02
0.100E 02
0.400E 02
0.400E 02
1.250 28.000 0.500
1.250 28.000 l.OCO
1-250 28.000 1.000
1.250 28.000 1.000
1.250 28.000 0.0
1.250 30.000 0.250
1.250 30.000 l-.OOO
1.250 28.000 5.000
VOL wu
VOLUME UNF-LOAD
(FT«*3I (LBS/DAY-MI)
°'° O.lOOE 03
0.0 o.O
0.0 o.O
0.0 o.O
o.o o.o
o.o o.o
0.400E 06 0.0
0.120E 09 0.0
JUNCTION ID, SECTIONS
MEETING AT JUNCTION,
AND DISTANCE
I 1-10. 0)( l-io. OK 0 o.ni
( 1 -5.0) (
1 1 -5. OK
1 2 0.0)1
2 -5.0) ( 0 0.0)
2 -5. OK 0 0.0)
3 0.0)1 0 0.0)
" CD
PT-LOAD FLOW-DAM
(LC5/DAY! (CCS)
o.o 100.00000
O.O n r\
«.u 0.0
o.o c.o
o.o o.o
DO
DEF-DAM
(fG/L)
0.20000
0.0
0.0
0.0
BDD
BOD-DAM
(MG/L)
0.20000
0.0
0.0
0.0
-------
5 1
6 9
7 1
8 5
9 1
10 5
11 <•
12 13
13 4
1* 11
C 20 CCC <
M3D2
C 20
M 2D
C 2D
M 2D
C DC
M1CM
C DC
MCM1
CCC (
ODD (
ODD (
EEE (
EEE (
FFF (
FFF (
GGG (
GGG (
6
2
3
3
A
4
5
5
6
6
O.O) (
0.0) (
5.0) I
5.0)1
7.5M
7.5) (
10.0) (
10.0) 1
-1.0)1
-1.0)1
a
3
4
4
5
3
8
8
7
7
O.O) 1
0.0) (
5.0)1
5.0)1
7.5) (
7.5)1
10.0) (
10.0) (
-5.0)1
-5.0) (
O
6
0
0
0
0
0
0
0
0
O.O)
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
C.O)
0.0)
O.O
0.0
sooo.oocoo
5000.00000
0.0
0.0
o.o
0.0
lOOu.UuOuO
1000.00000
O.O
0.0
0.0
0.0
0.0
0.0
0.0
0 . )
o.o
0.0
O.O
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
O.O
0.0
.0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-------
„, , SYMBOLIC 00 DEFICIT ANO BOD MATRIX
ROW JBV 1 2 3 « 5 ft 7 e 9 1O 11 12 13 I*
1 10 0AM 0AM -----_______
2 1 BBS BBB BBB BBS ----------
3 5 BBB BBB BBB BBB ---__-__,._
41 - - CCC CCC CCC CCC -_--___-
5 1 - CCC CCC - CCC CCC
6 9 - CCC CCC CCC CCC CCC CCC
7 1 - ODD OOP DDt> DDD ______
8 5 - DDD DDD DDO ODD ------
9 1 ------ EEE EEE EEE EtE
10 5 ------ EEE ECE EEE EEE
11 4 -------- FFF FFF - - - FFF
12 13 -------- FFF FFF - - - FFF
13 4 ---__--__. CGG GGC GOC
14 11 _--____._. GGO GGG GGG
-------
MODEL RUN ECHO
COMPUTED 0672
EXAMPLP PROBLEM 2
DATA GIVEN IN COMPATIBLE UMTS
SECTION KFO JUNCTION ID AND DISTANCE
AT ENDS
RK AK IK E FF TEMP PDEL
REM.RATE REA.RATE DEO.RATE DISP.RATE ULT/5DAY MILE-INCR
ll/DAY) U/DAYI ll/DAYI (MI2/DAY) (CENT) (Mil
1 GAFF DBLE
2 ECUP DBLE
3 ECJC DBLE
4 ECON DBLE
5 USPS DBLE
6 PUNC DBLE
7 PBAY S CM
8 ECBY S CM
SECTION
1 SAFF OBLE
2 ECUP DBLE
3 ECJC DBLE
4 ECDN DBLE
5 USPS DBLE
ft PUNC DBLE
7 PBAY S CM
8 ECBY S CM
SECTION
1 GAFF DBLE
2 ECUP DBLE
3 ECJC DBLE
4 ECDN DBLE
1 (BBB«
1 (CCC.
1 (DDDt
1 IEEE.
1 IFFF.
1 (CCC.
4 (GGG.
4 (SEA.
CROSS
SEC-AREA
(MI**2)
0.359E-04
0.359E-04
0.359E-04
0.359E-04
0.53BE-04
0.359E-05
0.359E-05
0.359E-05
U
1.63636 0
J.63636 0
1.80000 0
l.gOOOO 0
-5.0) (DAM. -10.0)
0.0) (BBB.
5.0) ICCC,
7.5) (DOO.
10. 0) IEEE.
0.0) (GGG.
-1.0) (BAY. -
15.0) (FFF.
-5.0)
0.0)
5.0)
7.5)
-1.0)
•10.0)
10.0)
0 BD
FLOW BENTHL-DMD
(MI»*3/DAY) (GMS/M2-DAY)
0.587E-04
0.587E-04
0.6'i6E-04
O.A46E-04
0.646E-04
0.587E-05
0.587E-04
0.646E-04
S V
.48309 -0.15581
.47219 -0.14492
.29570 -0.11570
.29570 -0.11570
0.2UOE 01
0.0
C.O
0.0
0.0
0.0
0.0
0.0
SD
0.38771
0.38771
0.23077
0.23077
0.376 0.117 0.342
0.342 0.117 0.342
0.342 0.117 0.3*2
0.342 0.117 0.342
0.342 0.117 0.342
0.370 0.244 0.370
0.370 0.061 O.?70
0.342 0.117 0.342
PR TIDE1
PHOTO-RESP TIDE-COEF
IMG/L-DAY)
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.500E 00
VD B C
-0.06044 -0.62723 -0.08684
-0.06044 1.43862 0.20615
-0.05077 C. 65495 0.98982
-0.05077 -0.18599 7.56833
5.000
5.000
10.000
10.000
15.000
5.000
5-000
5.000
HT
DEPTH
(MTRS)
0.762E 01
0.762E 01
0.762E 01
0.762E 01
0.914E 01
0.305E 01
0.122E 02
0.122E 02
G
-2*20547
4*04094
:«98756
-1.24445
1.250 28
1.250 28
1.250 28
1.250 28
1.250 28
1.250 30
1.250 30
1.250 28
VOL
VOLUME
0.0
0.0
0.0
0.0
0.0
0.0
0.272E-05
0.815E-03
H
-2.76414
1.49100
3.54438
16.75230
.000 0.500
.000 1.000
.000 1.000
.000 1.000
.000 0.0
.000 0.250
.000 1.000
.000 5.000
WU/AREA*RK
UNF-LOAD
(MG/L)
0.109E-04
0.0
0.0
0.0
0.0
0.0
0.0
0.0
ZZ FAC
5.18 -1.65
0.0 -1.90
0.0 -1.90
0.0 -1.90
-------
S USPS DBUE 1.2OOOO O. 19624 — O. 11624 O-137O1 -O.O57O1 -0 . 15956 5-17503 — 1-34O43 12.52564 O.O — 1.90
6 PUNC OBLE 1.63636 O. 48111 -0.15334 0.43848 -0.11120 -7.33666 8.98142 -30.76530 39.19626 0.0 -3.66
7 PBAY 5 CM 16.36363 3.29519 -0.02246 3.276*5 -O.OC372 5.94027 0.0 2-19978 . 0.0 0.0 -1.50
8 ECBY S CM 18.00000 3.61891 -O.C1891 3.60649 -0.00650 0.483013 0.0 0.8BBS7 0.0 0.0 -1.90
EQUATION JUNCTION ID. SECTIONS u
1
t
1
4
5
6
7
8
9
10
11
12
13
14
10
1
5
1
1
9
1
5
1
5
4
13
4
11
MEETING AT JUNCTION
AND DISTANCE
M DM
C 2D
M 20
C 2D
C 20
M3D2
C 2D
M 2D
C 2D
M 2D
C DC
M1CM
C DC
MCMi
DAM (
BBS (
BBB (
CCC (
CCC (
CCC (
ODD (
ODD (
EEE (
EEE I
FRF (
FFF (
GGG (
GGG I
1
1
1
2
6
2
3
3
4
4
5
5
6
6
-10.0) (
-5.0)1
-5.0)1
0.01 (
0.0) (
0.0)1
5.0)1
5.0) 1
7.5) I
7.5) (
10.0) (
10.0) (
-1.0) (
-1.01 (
1
2
2
3
3
3
4
4
5
5
8
8
7
7
-10.0) (
-5.0) (
-5.0)1
0.01 (
0.0)1
0.0) I
5.0) (
5.0) (
7.5) (
7.5) (
10.0) (
10.0) (
-5.0)1
-5.0) I
f
0
0
0
0
0
6
0
0
0
0
0
0
0
0
0.0)
0.0)
0.0)
0.01
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
0.0)
PT-LOAD
(MI*«3/3AY I (MG/L
0.0
0.0
0.0
0.0
0.0
0.0
0.90054
0.00054
0.0
0.0
0.0
0*0
O.OOC11
o.ooon
uu DO 6DD
FLOW-CAM DEF-PAM BOD-DAM
1 (MI«*3/DAYI IMG/U IMG/U
0.00006 0.20000 0.20000
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
-------
ECHO
EXAMPLE PROBLEM
COMPUTED 0672
SECT.
NO.
1
1
1
1
1
1
1
1
\
i
SECT.
NAME
GAFF
GAFF
GAFF
GAFF
GAFF
GAFF
GAFF
GAFF
GAFF
GAFF
GAFF
OIST.
MI
— — — ...
-10.0
-9.5
-9.0
-8.5
-8.0
-7.5
-7.0
-6.5
-6.0
-5.5
•» *\ rt
"!> • 0
BOD
V", / 1
*~o f L
•™— «— —
0.3890
0.4185
0 • 4 ^» ^> tt
0.4696
0.4913
0.5103
0.5267
0.5-02
0.5507
0.5579
0.5612
0.0 DEFICIT JUNCTION
MG/L JUNV.IION
I D
— - - 1. -••.
0.7683 OAM
0.8603
0.9503
1.C376
1.12ZO
1.2028
1 .2795
1.3515
1.4181
1.4785
1.5318 BBn
-------
ECHO
EXAMPL-E PROBLEM 2
COMPUTED 0672
SECT.
NO.
2
2
2
2
2
2
SECT.
NAME
ECUP
ECUP
ECUP
ECUP
ECUP
ECUP
OIST.
MI
-5.0
-4.0
-3.0
-2.0
-1.0
0.0
eoo
MG/L
0.5612
0.5857
0.6673
0.8350
1.1355
1 • 644 Q
D.O DEFICIT JUNCTION
MG/L ID
1.5318 BBB
1.6424
1.7816
1.9562
2.1676
2.4052 CCC
-------
ECHO
EXAMPLE PROBLEM ?
COMPUTED 0672
SECT.
NO.
3
3
3
3
3
3
SECT.
NAME
ECJC
ECJC
ECJC
tCJC
ECJC
ECJC
OIST.
MI
0.0
1.0
2.0
3.0
4.0
5.0
BOD
MG/L
1.6446
1.7620
1.9685
2.2H98
2.7606
3.42?9
D.O DEFICIT JUNCTION
MG/L ,D
2.4052 ccc
2.5229
2.6133
2.6625
2.6479
2.5347 r,Dn
-------
ECHO
EXAMPLE PROBLEM 2
COMPUTED 0672
SECT.
NO.
4
4
4
4
SECT.
NAME
ECDN
ECDN
ECDN
ECDN
DIST.
Ml
5.0
6.0
7.0
7.5
BOP
MG/L
3.A279
2.6831
1.8932
1.<.691
0.0 DEFICIT
KG/L
2.53<>7
2.2823
1.8632
1.6295
JUNCTION
ID
DDD
EEE
-------
ECHO
EXAMPLE PROBLEM 2
COMPUTED 0672
SCCT.
NO.
5
5
5
5
SECT.
NAME
USPS
USPS
USPS
USPS
OIST.
MI
7.5
8.5
9.5
10.0
BOO
MG/L
1.^691
1.0809
0.6B60
O.i.831
D.O DEFICIT
MG/L
1.6295
1 .35'.9
1.0572
0.8889
JLfNCT I OP4
ID
EEE
ITCC
-------
ECHO
EXAMPLE PROBLEM 2
COMPUTED 067?
SECT.
NO.
6
6
6
6
6
SECT.
NAME
PUNC
PUNC
PUNC
PUNC
PUNC
OIST.
MI
-1.0
-0.8
-0.5
•» f) 1
"*U • 3
0.0
BOO
MG/L
5.9*,03
".9655
3.9315
2.8283
! .f,ii."
D.O DEFICIT
MG/L
2.1997
2.2709
2.3254
2.3683
JUNCTION
10
GGG
ccc
-------
ECHO
EXAMPLE PROBLEM 2
COMPUTED 0672
SECT.
NO.
7
7
7
7
7
7
7
7
7
7
SECT.
NAME
PBAY
PBAY
PBAY
PBAY
PBAY
PBAY
PBAY
PEAY
PBAY
PBAY
DIST.
MI
-10.0
-9.0
-8.0
-7.0
-6.0
-5.0
-*.o
-3.0
-2.0
-1.0
BCD
MG/L
5,91.03
5.9*03
5.91.03
5.94C3
5.91.03
5.91.03
5. 91.C3
5.91.C3
5.9403
5 . 9 (. 0 3
0.0 DEFICIT
MG/L
2.1998
2.1998
2.1998
2.1991
2.1998
2 . 1 99 3
2.1998
2.1998
2. 199 =
2.1998
JUNCTION
ID
BAY
GGG
-------
ECHO
EXAMPLE PROBLEM 2
COMPUTCD 0672
SECT.
NO.
8
8
SECT.
NAME
ECBY
ECBY
OIST.
MI
10.0
15.0
ROD
MG/L
0.4831
0.4831
0.0 DEFICIT
MG/L
0.8B89
0.8809
JUNCTION
ID
FFF
SEA
-------
References
- 159 -
-------
1. Hydroscience, Inc.: Mathematical Models for Water Quality for the
Hudson-Champlain and Metropolitan Coastal Water Pollution Control
Project, FWPCA. Hydroscience, Inc., Westwood, N.J., April 1968.
2. Thomann, R.V., Systems Analysis and Water Quality Management,
Environmental Science Service Division, N.Y., N.Y., 1971,
p.p. 123-190.
3. O'Connor, D.J., Oxygen Balance of an Estuary. Jour. San. Eng.
Div. ASCE, Vol 86, No. SA3, May 1960, p.p. 35-55.
4. O'Connor, D.J., St. John, J.P., and Di Toro, D.M. , Water Quality
Analysis of the Delaware River Estuary. Jour. San. Eng. Div.
ASCE, Vol 94, No. SA6, Dec. 1968, p.p. 1225-1252.
5. Tracer, Inc.: Estuarine Modelling; An Assessment, Prepared
for Water Quality Office, Environmental Protection Agency.
Tracer, Inc., Austin, Texas, February 1971, p.p. 102-168.
6. O'Connor, D.J. and Di Toro, D.M., "The Solution of the Continuity
Equation in Cylindrical Coordinates with Dispersion and Advection
for an Instantaneous Release." Proceedings-Symposium on Diffusion
in Oceans and Fresh Waters, Lamont Geological Observatory.
Palisades, New York, 1965.
7. O'Connor, D.J., The Temporal and Spatial Distribution of Dissolved
Oxygen in Streams. Water Resources Research, Vol 3, No. 1, 1967,
p.p. 65-79.
8. "Effects of Pollution Discharges on the Thames Estuary", Water
Pollution Research, Technical Paper No. 11, Dept. of Sci. and
Ind. Res., Her Majesty's Stationary Office, 1964, 609 p.p. +
xxvll.
9. Courchaine, R.J., Significance of. Nitrification in Stream
Analysis - Effects on the Oxygen Balance, Journal Water Poll.
Control Fed., Vol 40, No. 5, Part 1, May 1968, p.p. 835-847.
10. Daily J. E. and Harleman D.R.F. Numerical Model for the Prediction
of Transient Water, ^ua^ity_,.in_Estuary_ Networks published by Md.T
October 1972.
11- O'Connor, D. J. and Mancini, J, "Water Quality Analysis of the New York
Harbor Complex" Journal WPCF, November 1972.
\L
-------
Appendix
ES001 Source Deck
(IBM 370)
-------
FORTRAN ..IV- 6.
0001
0002
0003
000*.
0005
0006
nno ?
0008
0009
0010
0011
0"12
0013
0014
001 ^
0016
nnj 7
OfHfl
0019
0020
0021
Or>2?
0023
nn2^
0025
00?^>
0027
r, n ? fl
0029
003n
0031
0032
0033
. 0034
0035
0036
0037
00*8
0039
. 0040
0041
LEVEL 2O MAIN DATE » 72258 l*/O3/
COMMON/HY/[SA<100),XSA<100>,IOJ03,AK(100),FF<100>,QUOO),BDUOO>,
-lWIJ(100),IMAX,IS2,ISe(100),XS8<100),aK(100),NJO(100),XJD<100>,
2S(100),JMAX,IS1,ISC(100),DK<100),NJU(100),XJU(100),KFO<100),
3VUOO),SD(100),B(100),G<100),AREA(100),TEMP(100),E(100),HT(100),
4JNO(100),W(100),VO(100),C(100),H(100),VOL(100),TIDE1(100),PR(100)
-5PDEL(100),Jt.iV(lOO)»00(100),DD(100),BDD(100),
6in(100),IDJ,FACUOO),NSV(100),U(100),IDATE,ZZ(100),XSC(100),JBM
COMMON/NE W/F AC 1«FAC2«FAC3»ITYPF(3'5)
DIMENSION tJM(100,100),DM(100,100),BMI(100,100),DMI<100,100)
DIMENSION BFOR(100),BSOL(100),DFOR(100). ,DSOL(1QO),BBXK250)
DIMENSION 00X1 (250), Mil 100), KJ( 100)
-EXTERNAL FR , FOR , FPR , FBPK , YD , YB , F A , FBA , FPA, FBPA
jnM=ioo
1 READ(5,10,END=999)IDJOB,IDATE,IS2,IS1,NN,IND,ICON,FAC1.FAC2,FAC3,
1ITYPE
DO 26 1=1,250
BBX! ( I )=0.
26 DOX1 I I }=0.
GO TO 25
_ 24 READ(5tlO,END=999)IDJOB,IDATEfIS2,ISltNN,INO,NOTH,FACltFAC2fFAC3f
11 TYPE
10 FORMATlA4,lX,A4,lX,5I3,Fc>.0,2F6.0t3X,35Al)
25 PRINT 157,lDjr)B,IDATE,IS2,ISl,NN,IND,ICON,FACl,FAC2,FAC3,ITYPE
L57 FORMAT! ' 1JOB ID AND DATE CARD 'Y1X , 20 ( 1H-J / IX « AS , IX » A8 , 51 5 , 3F8 .3.
11X/35A1 )
I NAX= I S2 + I SI
JKAX=2*I S2+IS1
CALL OATAR(NFLAC,MI,MJ)
PRINT 40, ITYPE, IDJ03, IOATE
30 FORMAT! ' 1 • .42X.35A1/' MODEL RUN ' , A8 , 2Xi ' COMPUT ED ' > 2X , A8 , 2X
1'OATA GIVEN IN COMPATIBLE UNITS'/I
40 FOrtMATl ' 1 ' ,42X,3I5A1/1 MODEL RUN ',A8,' COMPUTED ', 2X , A8, 2X,
130HDATA GIVEN IN ORIGINAL UNITS/)
CALL PR I ( N I ,MJ ,0 )
CALL PRJ(I"I,MJ,0)
P. Al 1 DAT AT I tf.DN )
IF(NFLAG) 6,6,5
*> PRINT 2nTNF'Ar,
20 FORMAT ( 1HO, 15, 2X.31HFRRORS FOUND COULD NOT CONTINUE /)
GO TO 1
6 IF(NN) 4,4,14
4 I I K = I
CALL KATS(FR,FRR,FPR,FBPR,IJK,BM)
CALL BRIGHT IBFOR)
DO 2 1=1, JBM
BSOLII )=BFOR(l)
DO 2 J=1,JBK
BMI ( [ , J J -BM ( J t I )
2 CONTINUE
CALL MATINV(BMI,JMAX,8SOL,1,BET,JMAX,IND,1)
CALL BCGH(BSOL,IS2, JMAX.KFO, B,C )
«,! _ PAGE
MAIN 001
MAIN 002
MAIN 003
MAIN 004
.MAIN 005
MA I N 006
MAIN 007
MA IN 008
MAIN 009
MAIN 010
MAIN Oil
MAIN 012
MA IN 013
MAIN 014
MAIN 015
MAIN 016
MAIN 017
MAIN 018
MAIN 019
MAIN 020
MAIN 021
MAIN 0?2
MAIN 023
MAIN 024
MAIN 025
MAIN 026
MAIN 027
MAIN 028
MAIN 029
MAIN 030
MAIN 031
MAIN 032
MAIN 033
MAIN 034
MAIN 035
MAIN 036
MAIN 037
MAIN 038
MAIN 039
MAIN 040
MAIN 041
MAIN 042
MAIN 043
MAIN 044
MAIN 045
MAIN 046
MAIN 047
MAIN 048
MAIN 049
MAIN 050
MAIN 051
-------
- FOR TRAM _tV_ &. LEVEL 2O
MAIN
DATE
72258
14/O3/41
PAGE OOO2
0042.
0043
. 0044
0045
0046
0047
004R
0049
___. ... 0050
0051
0052.
0053
0054
0055
.. 0056
0057
_. . .... 0053
0059
0060
0061
0062
0063
0064
0065
0066
0067
_ . . ... 0068
0069
0070
0071
007?
I F( ICON- 1)21, 22, 21
21 IJK=2
T.AI | MATS! FA» FPAiFPA tFBPAt IJKiOM)
CALL DRIGHT (DFOR)
nmi - 1 , .IRM
DSOL ( I )=OFOR( I )
nn 3 .1 = 1 , .IHM •
OMI (I,J) = DH( J,I)
3 CONTINUE
CALL MATINV(DMI,JMAX,OSOL,1,DET,JMAX,IND,2)
11 CALL RCGH(DSOL,IS2,JMAX,KFO,G,H)
22 PRINT 30, ITYPE, IDJOB, IDATE
CALL PRI (KL.MJ. 1 )
CALL PRI I (MI ,MJ, 1)
CALL PRJ(MI,MJ,1)
27 N=0
. . DO 13 1 = 1 i IMftX
13 CALL BDGUTl I ,N. I CON , BBX 1 , ODX 1 )
IF(ICQN-2)1,23,23
23 ICON=ICON+1
GC TOI1, 1,24, 27,1), ICON
14 IF( ICON-1 V28.28, 1
28 CALL BRIGHT(BFOR)
CALL MAM S(BMI,BFOR,BSOL,JMAX,JM AX, 1,100)
CAI 1 RCGhUBSOL. I S?f JMAX.KFO.B.C)
IF( ICON-1)29,22,1
_.. 29 CALL DRIGHT(DFOR)
C ALLMAMS(DMI, DFOR, DSOL, JMAX.JMAX, 1,1 00)
...... GO TO 11
999 STOP
END
MAIN
MA IN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
MAIN
052 ..... -- _-
053
054
055
056
057
058
059
060
061
062
063
064
065
066
067
068
069'
070
071
072
073
074
075
076
077
078
079
080
081
082
-------
FORTRAN— J-V
nnoi
0002
0003 _
0004
nnns
0006
0007
•0008
0009 ..
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
CL EVFI
i
i
i
i
i
!
1 1
1
12
13
2
20
BCGH
DATE
72258
14/03/41
PAGE 0001
.SUBROUTINE BCGH(SOL,IS2,JMAX.KFO,BGfCE)
DIMENSION SOL(l),BG11),CE(1),KFO(1)
= JMAX - 2 * IS2
D01N=1,IS2
N2=N*2
CE(N)=SOL(N2)
M=N2-1
QG(N)=SOL(M>
CONTINUE
0021=1,IS1
K=I+2*IS2
N=I+IS2
_KF=KFOIN)
GO TO(2,12il3,12),KF
BG(N)=SOL(K)
CElN)=0.0
GO TO 2
CE(N)=SOL(K)
BGLNJ=0.0-
CONTINUE
RETURN^ _ .. .. .._ .
END
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
BCGH
000
001
002
003
004
005
006
007
008
009
010
Oil
012
013
014
015
016
017
018
019
020
021
-------
—FORTRAN- IVL-G-LEVEU- 2O
BDOUT
DATE = 72258
PAGE O001
-CLQfll_
0002
-QOQ3.-
0004
0005
0006
_0007
0008
.0009
0010
0011
0012
SUBROUTINE BDOUT I I , N , I CON, BBX1, DDX1 ) _. ._
REAL*8 IDJOB,IDATE,ID,IDJ
COMMON/HY/ISA.( 100),XSA(100),IDJOB,AK(100),FF(100),Q(100),BD(100),
1WU1100),IMAX,IS2,ISB(100),XSB(100),RK(100),NJO<100),XJD(100),
.. _25(100J fJMAX,ISl,ISC(100),OKI 100),NJUl100),XJUI100),KFO(100),
3V(ICO).SD(IOO),B(100),Z(100),AREA(100),TEMP(100),E(100),HT(100),
___4JNOUOO),U100),VD(100),CI 100),H(1001,VOL(100),TIOEH100),PR( 100 )
5PDEL(100),JBV(100),QO(100),OD(100),BDD(100),
_.. 610(100),IDJ,FAC(100),NSV(100),U«1001,I DATE,II(100),XSC(100),JBM
COMMON/NEW/FAC1,FAC2,FAC3,I TYPE(35 I
DIMENSION BBX1I I ) ,DDXl(I)
IFIICON-4)40,41,41
_40..PRINT 30, ITYPE, IDJOB, IDATE
GO TO 44
41 PRINT 42
42 FORMAT!'1',40X,'COUPLED SYSTEM OF PREVIOUS 2 COMPONENTS'
PRINT 43, IDJOB,IDATE
43 rORMATIIOX,A3,21XOHCOKPUTED 5XA8//1HO,2X,'SECT. • , 1 4X ,
. 1'SECT.',15X,5HO1ST.,15X 3HBOD,13X4HO.O ,'DEFICIT1,12X,'JUNCTION'/
24X,'NO.',15X,'NAME',16X,' MI ',14X,' MG/L ',14X,' MG/L ',17X,MO'
. 33X, «—' ',14X,' ' ,14X,' ',13X,' •, 13X,11(1H-),12X,
48 I 1H-) )
BDOUTOOO
6DOUT001
BDQUT002
BDOUT 00 3
BDGUT004
BOOUT005
, BOOUT006
8DOUT007
BODUT008
BDOUT 009
BDCIUT010
BDOUTOll
BDOUT012
BDOUT013
BDOUT014
BDOUT015
8D(~IUT016
RDOUTOL7
BDOUT018
/BDOUTO 19
0013
0014
onis
0016
0017
0018
_ 0019
0020
0022
0023
0024
0025
0026
0027
0028
._ 0029
0030
— _ .... 0031
0032
0034
0035
0036
0037
0038
0039
0040
0041
0042
44
4
5
6
7
20
. 8
_ 9
32
31
14
10
. . .
IF(PDELU)) 5,5,4
DELTA = PUEL1I)
XS=XJD(I)
XE-XJU(I)
JJS=NJU(I)
J JE=NJO( I )
CALL BODO( I ,XE,BBX,DDX)
BDOUT021
BDOUTO?.2
BDOUT023
BDOUT024
BOOUT025
BDniJT026
BDOUT027
BDOUT028
BOGUT029
IFI iccN-4)6,7,7
DDX1 (N)=DOX1 (NH-ODX
8BX1 ( N)=8RXltN} +B6X
PRINT 20, I , ID( I ) ,XE,BBX, ODX, JJS
GO TO 8
PRINT 20, 1,101 I ),XE,BBX1(N),DDX1(N),JJS
FGRHAT(1HO,I5,16XA8,F16.1,2F20.4,17XA3)
XE=XE+OELTA
IF(XE-XS)9, 10,10
CALL BOOOI I ,XE,BBX,DDX! _
N=N+1
IFI ICON-4132, 14, 14 _________
BBX1 1N)=BBX1 (Nl+BBX
DDX1(N) = DOXKN)+DDX . ....
PRINT 20, I ,ID( I ) ,XE,BBX,DDX
GO TO 8 _
PRINT 20, I.IDU ),XE,BBX1(N),DDX1{N)
GO. TQ 8 ________________ .
CALL BODOI I ,XS,8BX,DDX)
IFUCON-4116,17,17
snouTon
BDOUT032
BDOUT033
FJOmiT034
BDOUT035
BDOUT036
BDOUT037
BOOUTnifi
BDOUT039
BOOUT040
BOOUT041
BDOUT042
BOOUT043
BDOUT044
BOOUT045
BDOUT046
BDOUT047
BOOUT048
BDOUT049
BDOUT050
BDOUT051
-------
FORTRAN-IV-G LEVEL - 2O
RDOUT --
DATE
72258
PAGE OOO2
00*4
0045
00*6
00*7
00*8
30
16-CDX1(N)=DDX11N)+DOX
BBX1(N)=*BBX1(N)+BBX
_. BDOUT052-
BDOUT053
BDOUT05*
BOOUT055
BOOUT056
PRINT 20,I,ID(I),XS,BBX,DDX,JJE
GO TO 18
— 17 PRINT 20, I,ID(I) ,XS,BBXKN) ,DDXUN),JJE
30 FORMAT(1H1,*3X,35A1/10X,A8,21X8HCOMPUTED 5XA8//1HO,2X,«SECT.' , 1*X,8DOUT057
L'.SECT. V,.15X,5HDIST. ,15X 3HBOO, L3X4HD.Q , ' DEF 1C 11' , 12X, • JUNCT I ON1 / BDOUT05B
2AX,'NO.•,15X,'MAME1t!6X,' Ml «,lAXt' MG/L '.l^X,' MG/L ',17X,•ID*/BDOUT039
33X, ' '.UX,' -• , 14X,' 'tlSX,' ',13X,11I1H-),12X, 8DOUT060
48(1H-)) BDOUT061
.-00^9 18-RETURN - . . .. BDOUT062
0050 END BDOUT063
-------
FDR TR AM T W C
nnn\
0002
0003
•
0004
0005
0006
Ofif)7
0008
0009
0010
0011 .
0012
nni^
0014
_S 00 15
SUBROUTINE BOD01 I , X , BBB , ODD)
REAL*0 IOJOB, IDATEi in, IDJ
CCMiMON/HY/ISA(100),XSA(100)tIDJOB,AK(
lWUUOO),IMAX,IS2,ISimoO),XS8(100),RK
2S ( 100) | JMAX > IS1 1 ISCl 100) r DK ( 100) iNJUf
3V(100),SO(100),B(100),G(100),AREA(100
4JNQ( inn) ,w( loo) ,VO( ino) ,<~ 110") tH( 10n)
5PDELUOO) , JEW (100) , QD( 100 ) , DD( 100 ) ,BO
6IO( 100) , IDJ ,FAC( 100) ,NSV( 100) fU( 100) ,
EXTERNAL FR, FBR, FA , F3A
BBM=OMtKFO,S,C.FR,FBR, I,X)
YDIX = Z7.( I )+FAC( I )*BnM
RPP=BPM+YP( i , x )
DDM=OM(KFO,G,H,FA,FBA,I,X)
-. - IF(KFOU)-4)1,2,2
1 DDD=DDM+YDIX
GC TO 3
2 DCO=DOM
•\ CriMTINIIE
RETURN
END . .
-DATE-= 72258 1W03/41
Bnon ooo
100),FF(100),Q(100),BD(100)
(100),NJ01100),XJD(100),
100),XJIM100),KFOUOO),
),TEMP(100),E(100),HT(100),
,VOL(100),TIDE1(100),PR(J.OO
0(100),
IOAT.E,ZZ.tlOO),XSC(100),JBM
BODO 001
, BOOO 002
BODO 003
BODO 004
BOOR 005
1,6000 006
BODO 007
BODO 008
BODO 009
BODO 010
BODO Oil
. BODO 012
BOOO 013
BOOO 014
BODO 015
BOOO 016
BOOO 017
BODO 018
BODO 019
BODO 020
_ PAGE O001._
-------
00
-------
FOR TRAN-.-IV-G -LEVEL 2O
DATAC
DATE
72258
14/03/^1
PAGE OOO1
rmni
0002
0003
0004
0005
0006
0007
0008
0009
0010
00 1 1
0012
On I A
0014
0015 _..
0016
0017
0018
Qn\q
0020
0021
0022
._.._ .. 0023
0024
nr^s
0026
0027
0028
nn?q
0030
non
0032
0013
0034
0035
0036
nm7
0038
0039
0040
0041
0042
0043
00*4..
0045
SUBROUTINE'DATAC(ICON)
REAL*8 IDJOB, IDATE, I0t IDJ
COMMON /HY/ISAt 100) »XSA( 100) •IDJOcl|AK(100)iFF(100)iQ(100)tBD(100)t
1WU(100),IMAX,IS2,ISB(100),XSB(100>,RK(100),NJD<100),XJD(100),
5 c / t n ni iMAv.rci Tcrfinn* ntffirmi-MiiirinrM-Yiiifinnt-tfpnMnni-
3V (100) ,50(100) ,R( 100) ,G( 100) , AREA! 100) , TEMPI 100), E( 100), HT ( 100),
4JN01 100) iW( 100 1 iVD( 100) , C ( 100) ,H( 100) , VOL ( 100) i TIDEl ( 100) , PR ( 100) ,
5POEL(100),JBV(100),OD(100),DD(100),BOO(100),
6ID(100),IDJ,FAC(100),NSV(100)tU(100),IDATE,ZZ(lOO),XSC(100),J8M
CUMMON/NEW/FACltFAC2fFAC3, I TYPE (35)
— DIMENSION LINE(IOO)
DATA LITERL/3H -/
DO^M -- I, I MA*
AKm=AK =qn I i )*TF
5 CONTINUE
. OH32 1=1, If AX
FACI I )=DK(I)*FF( I)/(AK(I)-RK( I.))
IK I ) =0.0
31 AA = BO( I ) /HTI I 1
IF(RMT>)!4rl4,l?
12 Bn=OK(I)*WUlI)*FF(I)/(RK(I)*AREA (I))
r, n in is
14 BB=0.
15 ZZm = (AA + PRm+BB)/AK(I) -
32 CONTINUE
pn i i~\ , IMAX
U(I )=Q(I)/AREA(I )
CALL SVD(U(I),E(I)tAK(I),sn(I)tvn(ri)
CALL SVDIUm.EU) tRKI I ) , S ( I ) , V ( I ) )
1 CONTINUE
PRINT 20, (J,J=1,JMAX)
?n F"RMAT( • I • ,38X, ' SYMBOL 1C Pfl PFFIf-TT AND RDD MATRIX' /4X. 'RDW, 2X,
1 ' JBV .24I4/13X.24I4)
SOI IS2M=1<;2*2 . .
00 3 JJ=1,JMAX
DATATOno
DATAC001
DATAC002
DATAC003
DAT AC004
OATAC005
DATAC006
DATAC007-
DATAC008
DATAC009
DATAC010
DATAC01 1
DATAC012 - . --
DATAC013
DATAC014
DATAC015
DATAC016
OATAC017
DATAC018
DATAC019
DATAC020
DATAC021
DATAC022
DATAC023
DATAC024 . „ ,
DATAC025
OATAC026
DATAC027
DATAC028
DATAC029
OATAC030
DATAC031
OATAC032
DATAC033
DATAC034
OATAC035
DATAC036
DATAC037
DATAC03R
DATAC039
DATAC040
DATAC041
DATAC042
DATAC043
DATAC044
DATAC045
DATAC046
OATAC047
DATAC048
DATAC049
DATAC050
DATAC051
-------
FORTRAN.
0(14 A
0047
_ 0048
0049
0051
0053
0054
0055
0056
0057
OOSFi
0059
0060
0001
0067
0063
0064
0065
J=INH(JJ)
DO 2 L=1,JBM
? L INF f 1 ) =L I TFRL
CALL PLACE! J,JJ,IS2H, ISA, LINE)
CALl PLACE ( J,JJ, I S2M, ISB, LINE )
CALL PLACEtJ, JJ, IS2M, ISC, LINE)
PRINT 30, JJ,JB'V( JJ), (LINFI I ) , 1 = 1, JMAX)
30 FORMAT(1HO,2I5,2X,24{1XA3)/13X,24(1X,A3))
3 f.ONTINUF
8 DO 9 1=1. JMAX
OD( I)=0.
set 11=0.
VD( I 1=0.
G( I 1=0.
HI I 1=0.
FAC( I )=0.
9 ZZ( I 1=0.
11 RFTURN
END
PAGE OOO2
DATAC052
OATAC053
DATAC054
DATAC055
DATAC056
DATAC057
DATAC058
OATAC059
DATAC060
DATAC061
OATAC062
DATAC063
DATAC064
DATAC065
DATAC066
DATAC067
DATAC068
DATAC069
DATAC070
OATAC071
-------
FORTRAM -IV- G_L£VEU 2O
DATAR
DATE
72258
PAGE 0001
-OOOl.
0002
0003-
0004
0005
0006
._0007_
0008
0009
. 0010
..0011..
0012
0013
0014
.. 0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033. __. _..
0034
_ 0035
0036
0037
0038
40
71
11
14
700
73
701
. 74
75
702
6
SUBROUTINE DATAR(NFLAG.MI,MJ) . _
REAL*8 IOJOB,IDATEfID,IDJ
__ COMMON/HY/ISA(100),XSA(100),IDJOB,AM100),FF(100),Q(100),BD(100),
1WUI100),IMAX,IS2,ISB(100),XSB(100),RK(100),NJD(100),XJD(100),
2S!100),JMAX,ISl,ISCUOO),DM100),NJU(100),XJU(100),KFO(100),
3V(100),Sr>(100),B(100),G(100),AREA(100),TEMP(100),EJ100),HT(100),
_4JNO(100),W(100),VD(100),C(100),H(100),VOL<100),TID~E1<100),PR(_100)
5PDEL(100) ,JBV(100),QO<1001,00(100),BDO(100),
6IDI 100) , IOJ,FAC(100),NSV1100),U(100),IDATE,ZZ(100),XSC(100),JBM
COMMON/NEW/FAC1 ,FAC2,FAC3,ITYPE(35)
DIMENSION !OJB(3,14),ISV(3),IDHJBI14),IDKFO(6),MI( 1 I ,MJ(1)
DIMENSION LI(3)',NDA(3)
__. DATA IDHJB/'C 20','C DN'.'C OP'.'C OC'.'M 20',«M DN'.'M DP',
1'M3D1•,'M302','M DM' , 'MCMl • , 'MCM2' , 'M1CM1 , 'M2CMV
... DATA lOKFG/'DBLE', ' S NI','S PI'.'S CM ' , ' MI SS ' , ' N CL'/
DATA IPASS/'PASSV, IFAIL/' FAIL' /, IBLNK/' '/
DATA IOJB/2*l ,5,2,1,5,1,3,5,1,4,5,2*1 ,5,2,1,5,1,3,5,
26*1,1,1,5,1,4,5,2*1,4,1,4,5, 2*1,4/
. ._. NFLAG = 0
KNAL3=I(31NK
PRINT 158
FORMAT!' INPUT SECTION CARD L I ST I NG'/ IX,26( 1 H- ) )
__. DO 6 I 1 = 1, IMftX
RE4D10,NS2, I,NUSD,NJU(I),XJUU),AREAlI),RMI),AMt),FF_CONTINUE_
00 72 11=1,IMAX
IST=NTPStII) .
FOR SECTION',13)
I .OR.RK I ,2X I10.2E20.4)
I ,9X,UO,2E20.4)
0040
MHII) = IDKFO(IST)
OATAROOO
DATAR001
OATAR002
DATAR003
DATAR004
DATAR005
.DATAR006
DATAR007
DATAR008
DATAR009
DATAR010
DATAR011
DATAR012
DATAR013
OATAR014
DATAR015
DATAR016
DATAR017
DATAR018
DATAR019
DATAR020
DATAR021
DATARO?2
DATAR023
DATARO?4
DATAR025
DATAR026
10ATAR027
DATAR028
1DATAR029
DATAR030
DATAR031
DATAR032
DATAR.033
DATAR034
OATAR035
DATAR036
DATAR037
DATAR038
DATAR039
DATAR040
DATAR041
OATAR042
OATAR043
DATAR044
DATAR045
DATAR046
DATAR047
DATAR048
DATAR049
OATAR050
DATAR051
-------
-FORTRAN-
QATAR
DATE =. 72258 . .
.J>AGE OO02
—0041
0042
-X)043
0044
-0045
0046
0047
0048
0049
no™
0051
0052
0053
0054.
0055
0056
0057
0058
0059
0060
0061
0^62
0063
0064
0065
0066
0067
0069
. 00/0
0071
0072
0073
0074
0075
-0076
0077
- 0078
0079
0080
0081
0082
0083
0084
50
8
81
82
AT
84
1
5
3
86
87
88
703
I
4
Fl
Fl
Ci
D
J
I
D
G
I
G
I
G
J
I
N
I
L
N
I
I
I
I
L
G
N
N
C
I
I
P
F
3
C
R
E
—72-PRINT .700, I I,IDKFOt1ST) .
PRINT 159
-159 FORMAT(« INPUT JUNCTION CARD LI STING'/IX,27(1H-))
DO 8 JJ=1,JMAX
-10 FORMAT(2X,2I3,I2, A3 ,F5.0,F7.0,3F5.0,4F8.0,4X,A4) .
20 FORHAT(1HO,3I4,3X,A3,F10.3,1X,F9.1,3F10.3,F9.1,F10.3,1X,610.3,
1F10.3.5X,A4/)-
REftD 50,J,JNO( J) ,JfW(J),ISA( J) ,XSA( J) , ISB(J) , XSB( J } , ISC ( J ) ,XSCU)
-_ IWtJ),OD
-------
-•FORTRAN—IV—G-L€VEU 2O
DRIGHT
DATE = 72258
1WO3/A1
PAGE OQO1
nnni
0002
0003
0004
000*5
0006
nnr>7 .
0008
OO'iq
0010
0011
0012
nnn
0014
0015
0016
nni v
ooin
n n i q
0020
0021
00?7
0023
nn?4
0025
0026 .
0027
00?8
0029
nn^o
0031
oni?
0033
0034
0035
nnift
0037
0038
0039
f)n<.0
SUBROUTINE ORIGHT( FDR )
DIMENSION FOR(l)
a £ A I * a ininn TDATC rn tni
COMMON/HY/ISA(100),XSA(100>,IDJ08,AK<100),FF(100),Q< 100 ) , BDt 100 ) ,
1UIIMOO).TVAV T^7 T ^ R M O H 1 V 9 R ( 1 O CM RKMHO) M 1 D I1 1 O ("M V 1 H 1 1 O CM •
2S(100),JMAX,IS1,ISC( 100) , DM 100 ) ,NJU< 100 ) , X JU ( 100 ) ,KFO(100),
3V(100)tSDtlOQ)iB(100)»G(100)pAHEA(100)iTEMP(LOO)fE(100)iHTtlOO)f
4JNO(100),W(100),VD(100),C(100),H(100) ,VOL( 100),TIOEl(100)fPR(100)
5PDELI 100) tJDVt 100) tQDl 100) .DD( 100) tBDD( 100) ,
6ID(100),IDJ,FAC(100),NSV(100),U(100),IDATE,ZZ1100),XSC(100)tJBM
DO 102 I- 1 t JBM
102 FCRI I )=0.0
PQ QQU=V,JMAX
JB=JBV( J)
[ - ISA ( J )
X=XSA( J)
K =ISO(J)
Y=XSB(J)
KK =i«;r(.n
Z=XSC(J)
GO TO ( 1, 1. It 4,5,5,5,5,5, 10t 11. 12i 13f 14) , JB
1 FOR( J)=YD(K,Y) -YDI I,X)
r,r in qq
4 FOR( J)=-YC( I,X)
r.n -in qq
5 FOR! J )=Q( I )*YD( I,X)- YDPtI,X)+ YDPU.Y)
1- Q(K)*Yn(K,Y)
IF (JB-3)9q,8,9
8 FHRIJ) =FOR(.I)+ YDP(KK.Z) - 0(KK)*YO IKK.Z)
GO Ti]99
9 FHR ! .1 ) =FRR [ .) ) +0 ( KK ) * YD ( KK , 7 )- YDP(KK.Z)
GO TO 99
10 FnR(.ll= YDP(K.Y)-OIK)*YO(K,Y)+GD( J )*DD( J)
GO TO 99
11 FCR[J)-YDPU,X>-Om*Yr)(I,X)+FF(K)*DMK)*B(K)*VOL
-------
0001
0002
OOQ'*
0004
. 0005
0006
00-^7
0008
0009
0010
0011
0012
nn\T.
0014
-. -.- 0015 „ - - . - ._
0016
0017
0018
OOIQ
0020
_ 0021
0022
0023.
00^4
nn?s
0026
0027
0028
0029
0030
rum
0032
0033
0034
0035
0036
0037
0038
0039
0040
p
R
c
1W
?s
3V
4 1
5P
IS I
E
F
F
R
E
F
F
- R
E
F
F
R
E
F
F
R
E
F
F
. P
E
F
f
F
E
f
f
F
1
1
[
1
I
I
1
FUN
DATE = 72258
1WO3/41
PAGE 0001
-FUNCTION FUN(ItX).... ....... . . ._ ____________ ..... ___ . ...... ___________
REAL*8 IDJOB, IDATE, ID, IDJ
COKMQN/HY/ISAt 100),XSA(100),IDJOB,AK(100),FF(100),Q(100),BD(100),
1WU(100),IMAX,IS2,ISB(100),XSB(100),RK(100),NJDUOO),XJD<100),
2S(100),JMAX,IS1,ISC( 100) , DK I 100 ) , NJU ( 100 ) , X JU ( 100 ) ,KFO(100),
3V(lOO)f SD( 100), B( 100), G( 100), AREA ( 100) , TEMPI 100) , E(IOO) ,HT< 100),
.4JNO(100),t<(100),VO{100),CilOO),HtlOO),VOL(100),T.IDEl(100),PR( 100 J
5PnEl(100),JBV(100>, 00(100), 0011001,600(100),
61 D( 100) ,IDJ,FAC<100),NSV(100),U(100),IDATE,Z2(100),XSC(100),JBM
ENTRY FR (I,X)
FUN=F(S,V, I,X)
= FUN
RETURN - _____ ___________ __________ ..... ______ . __ _______ - -
ENTRY FA (I,X)
FUN=F(SD,VD,I,X)
FA = FUN
IRN
ENTRY FBR (I,X)
= FB(S,V,I,X) .... ..... ___________ _ ...... _______ . .. .. -
Ffi<* ' = FUN
RETURN
ENTRY FBA ______________________ ------------------------- - -
YBP = FUN
— RETURN ______ ________ . - ...... - -
END -
FUN
FUN
FUN
FUN
FUN
FUN
, FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN
FUN'
FuN
FUN
FUN
..000
001
002
003
004
005
006
007
008
009
010
Oil
012
013
014
015
016
017
01.8
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
035
°36
°39
040
FUN
FUN
FUN
042
043
044
045
-------
FOUTRAN. IV-G -LEVEL 2O
MAMS
DATE- =. 72258
PAGE OO01
—0001
0002
_0003
0004
- 0005
0006
—£007 _
0008
_ 0009- _
0010
SUBROUTINE .MAMS -( -A,B.,C,L1, LZ, L3.NRD LM)
DIMENSION A(NROIK,l),B
-------
20
MATIMV
DATE -= -72258
-14/03/41
PAGE- 0001 ...
0001-- -
0002
0003
0004
0005 — -
0006
0007
0008
0009
0010
0011
001?
001^
0014
0015
0016
0017
0018
nn^q
0020
on?i
nn?2
0023
002', __.
0025
0026
0027
nn?g
0029
. 0030
0031
0032 . -
0033
0034
0035
0036
no 17
0058
003°
0040
0041
0042
0043
0044
00^5
0046
SUBROUTINE MATINV ( A >N » B » M,DETER.Mf JMAX * IND» IKND 1
DIMENSION IPIVOT(100),A(100fl),B(100,l), INDEX! 100, 2 ) , PIVOT ( 100)
IF( IND) 10, 10, 1003
1003 GO T0( 1004, 1006) , IKMO
1004 PRINT 1005 .. - -
1005 FORMAT( ' 1 ' ,51X,' INITIAL BOD MATRIX1
GC TO 1008
1006 PRINT 1007
1007 FORMAT! • 1 • ,48X , • INI TI AL 00 DEFICIT MATRIX'
1008 PRINT 1000, (J,J=l, JMAX )
1000 FORMAT! /4X.9I13/)
DO 1001 1=1, JMAX
inni PRINT ?nont i , ( A( i ,.j ) , j=i, JMAX)
PRINT 1009
-1009 FORMAT!/' FORCING FUNCTIONS'/)
PRINT 2000, I, (81 J,l ) , J = l, JMAX)
2000 FORMAT! I 3, IX, 9E 1 3. 5/ ( 4X , 9E 13. 5) )
10 OETERM=1.0
is nn ?o .1 = 1 ,N
20 IPIVOT(J)=0
3onn550i-i,N
c
£ SEARCH FOR PIVOT ELEMENT . ...
C
40 AMAX-0.0
45 DO 105 J=1,N
50 IF ( IPIVOT(J)-l) 60, 105, 60
60 PC 100 K=1,N
70 IF ( IPIVOT(K)-l) 80, 100, 740
80 IF ( ABSI ANAXJ-ABS! A 1 J,K) ) ) 85, 100, 100
fi5 IRRW = .I
90 ICDLUM=K
95 AMAX=A(J,K1
100 CONTINUE
105 CONTINUE
110 IPIVOT(ICOLUM)=IPIVOT( ICOLUM1+1
r.
C INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
C . ._.....
130 IF (IROW-TCOLUM) 140, 260, 140
140 DETERM=-DETfRM - - -. .
150 DO 200 L=1.N
\f-,r <;UAP-A ( IRT.W,! )
170 A ( IROW,L)=A( ICOLUM.L)
2nnfl(irnLiiM,L) = <;vjip . . . . . .
205 IF(M) 260, 260, 210
2in nn 2^n i = I , M . .
220 S'^AP = B( IRCW.L)
2^0 B 'lRnut|_ )-p( irOLUM.L)
250 B , ICOLUM,L) = SWAP
260 I !DEXt 1 , 1 ) = IROW . . __ - - —
270 I ,DEX(I,2) = ICOLUM
MAT IMOQO
MATIN001
MATIN002
MATIN003
MATIN004
MATIN005
MATTN006
MATIN007
MATINOOB
MATIN009
MATIN010
MATT NO 11
MAT I NO I 2
MATIN013
MATT NO 14
MATIN015
MATIN016
MATIN017
MATFN01 fl
MATIN019
MATIN020
MATIN021
MATIN022
MATIN023
MATIN024 „ ._. ... ...
MATINO?5
MATIN026 __._ .
MATIN027
MATINO?8
MAT IN029
MATIN030
MATIN031
MATIN032
MATIN033
MATIN034
MATIN035
MATIN036
MATIN037
MATIN038
MATIN039
MATIN040
MATIN041
MATIN042
MATIN043
MATIN044
MATIN045
MATIN046
MATIN047
MATIN048
MATIN049
MATIN050
MATIN051
-------
- FORTRAN-I.V_G_LEVEL-. 20
MATINV
BATE = 72258
PAGE 0002
0048
0049
0050
0051
..... 0052 L .
0053
0054
0055 . _..
0056
0057
0058
0059
0060
0061
0062
- _ 0063
0064
0065
0066
0067
0068
0069
. 0070
0071
. .... 0072
0073
on 74
0075
0076
0077
0078
0079
noao
0081
0082
0083
0084
0085
0086
0087
ooee_
0089
320
c
C
_ C
330
340
350
355
360
_ .370
C
C
C
300
390
.. - 400
420
430
450
_-. . 455
460
_. ._ 500
550
f.
C
._C
600
610
620
630
640
650
660
. ._ 670
70U
705
710
740
750
1010
1011
1012
1014
1015
1016
1002
U>17
PJVOTd )=A(ICOLUM,ICOLUM) .. ..
DETERM=OETERM*PIVOT(I)
DIVIDE PIVOT ROW BY PIVOT ELEMENT
AUCOLUM, ICOLL|M) = 1.0
00.350 L=1,N .
A(ICOLUM,L!*A( ICOLUM.D/PIVOTI I)
IF(M) 380, 300, 360
DO 370 L=1,M
6(1COLUM,L)=B(ICOLUM,L!/PIVOT(I)
REDUCE NON-PIVOT ROWS
DO 550 L1=1,N
IF(Ll-ICOLUH) 400, 550, 400
T = A (LI,ICOLUM)
AlLltICOLUM)=0.0
.DO 450 L=liN
A(L1,L)=A(L1,L)-A(ICOLUM,L)*T
IF(M) 550, 550, 460
00 500 L=1,M
BU 1, L)=3(L1,L)-B(ICOLUM,L)*T
CCNTINUE
INTERCHANGE COLUMNS
DC 710 1=1, N
L=N+1-I
IF I INDEX ( L, 1 1-INDEXI L,2) )
630, 710, 630
JCOLUM=INDEX(L,2)
DO 70S K=1,N
SWAP=A(K, JROW)
A (K, JROW)=A(K, JCOLUM) _ _
A (K, JCGLUMJ=SWAP
CCNTINUE ., .„ __ _______ .......... _______________________
CONTINUE
CHNITINUF _ .
CONTINUE
IF( INO) 10?0, 1020, 1010
GO TOl 101 1,1014) , IKNO
PRINT 1012 _ ________ . ____ ________________
FORMATJ///50X, • INVERTED BOO MATRIX"
GO TO 1016
PRINT 1015
FORMAT(///47X, • INVERTED DO DEFICIT MATRIX1
PRINT 1000, ( J, J=1,JMAX)
00_1002 I=1,JMAX ___ ______________
PRINT 2COOfIi«A(I,J)tJ=lf JMAX)
PRINT 1017
FORMATl/1 SOLUTIONS FOR UNKNOWN COEFFICIENTS'/)
MATIN052
MATIN053
MATIN054
MATIN055
MATIN056
MATIN057
MATIN058
MATIN059
MATIN060
MATIN061
MATIN062
MATIN063
MATIN064
MATIN065
MATIN066
MATIN067
MATIN068
MATIN069
MATIN070
MATIN071
MATIN072
MATIN073
MATIN074
MAT IN075
MATIN076
MATIN077
MATIN078
MATIN079
MAT IM080
MATIN081'
MATIN032
MATIN083
MAT IN084
MATIN085
MATIN086
MATIN087
MATIM088
MATHJ089
MATIN090
MATIN091
MATIN092
MATIN093
MAT1N094
MATIN095
MATIN096
MATIN097
MATIN098
MATIN099
MATIN100
MATIN101
MATIN102
MATIN103
-------
-- FORTRAN— IV-G-LEVEt- 20
._ . MATINV
DATE ._
-------
-PORTRAN—IV—G-LEVEL 20
MATS
DATE = 72258
14/03/41
PAGE 0001
J3004-
0002
-0003
0004
0005
0006
Qnn/
0008
0009 . _.. __.
0010
001 1
0012
nn^
0014
0015
0016
0017
.0018
nniq
0020
.. 0021
0022
0023
0024
nr>7S
0026
0027
0028
0079
0030
nmi
0032
0033
0034
0035
0036
0037
0038
0039
0040
.0041
0042
0043
0044
0045
0046-
E
r
c
t
101 1
[
]
j
i
»
>
i
i
i
i
i
i
i
?
•^
u
5
6
--SUBROUTINE MATS ( F , FB , FP , FBP , UK, TM )
REAL*8 IDJOB,IDATE,10,IOJ
COMMO^/HY/ISAf100),XSA!100),IDJOO,AK(100),FF(100),0(100),BD(100),
1WU(100),IMAX,IS2,ISB(,100),XSB(100),RK(100),NJD(100),XJD(100),
2S(100),JMAX,IS1,ISC<100),DK(100),NJU(100),XJUI 100).KF01100),
3V(100),SD(100),B(100),G(100),AKEA(100),TEMP(100),E(100),HT(100),
-4JNOI100),W(100),VD(100),CUOO),H<100),VOL(100),TIDE11100),P.R(100)
5PDEL(100),JBV(100),00(100),DD(1001,000(100),
-610(100), IDJ,FAC(100),NSV(100),U(100),IDATE,ZZUOO),XSC(100),JBM
EXTERNAL F,F8,FP,FBP
DIMENSION TM(100,1)
DO 101 1=1,JBM
00 101 J = 1,JBM ..
TM(I,J)=0.0
00199 J=1,JMAX -
JB=JBV(J)
=ISA(J)
= XSMJ)
= NQSECT.(I..,IS2J - — _ _
IM=L-1
K =ISB(J) . .
Y=XS8(J)
M=NOSECT(K ,1S2)
KK=M-I-
KK =ISC(J) _.. _. ._ . .. _ ...
Z=XSC(J)
N=NCSECT(KK ,IS2)
KKM= N-l
GOTO!1,2,3,4,5,6,7,8,9,10,11,12,13,14),JB
IM,J)=F(I,X)
TfML, J)=FB(I,X) . ___ . ..
I (/(KM, J)=-F(K,Y)
TfMM, J)=-FB(K,Y)
GO TO 99
TM(L,J)=F(I,XI
TMIKM,J)=-F(K,Y)
IMIM, JI=-FB(K,Y) -
GC TO 99
TM(IM,J)=F(I,X) ...
TK(L,J)=FB1I,XI
TM
-------
.FORTRAN.-.IV-G-J-EVEL 2O
MATS
DATE
72258
PAGE OOO2
"047
0048
0049
0050
0051
0052
_ __ . 00 S3
0054
0055
0056
0057
00 50
0059
0060
0061.
0062
0063
0064
0065 .
0066
. 0067 . . _
0068
0069 . .
0070
_- 0071,
0072
0073
0074
C075 .. . _ ._.
0076
... 0077 _
0078
0079 _._ .
0080
-. . _ 0081
0032
0033
0084
_ _ . 0085 ..
0086
0087 .
0038
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
IMIKM,J1-Q(K)*FIK,Y)-FP(K,Y)
TM(M,J)=0(K)*FB(K,Y)-FBP(K,Y)
GO TO 99
7 TMI IM,J)=FP( ,X)-Q( I )*F( I,X)
TN(L»J)=FBP( , X )-Q ( I ) *FB ( 1 , X )
TMIM ,J)=Q(K *F3(K, Yl-FBP (K,Y)
r,n TO 99
8 TM(IM,J)=FP( f Xl-QI I )*F( I ,X)
TM(L,J)=F8P( ,X)-Q(I)*FBU,X)
TM(KM,J)=Q(K *F(K,Y)-FP(K,Y)
TM(M,J)=Q(K)*rB(K,Y)-FBP(K,Y)
28 TM( NiJ) = Q(KK) *FB(KK,Z)-FBP(KK, Z)
TM(KKM,J)=Q(KK)*F(KK,Z)- FP(KK, Z)
GO TO 99
9 TMIM,J)=FPU,X)-Q(I)>t= FP(KK,Z)-Q(KK )*F (KK,Z)
TfM N,J)= F8P (KK.Z)-Q(KK ) *FB (KK,Z)
. .. GC TO 99
10 1C (KM, J )=C!K)*F(K, Y)- FP{K,Y)
TMM, J )=Q(K)*FB(K,Y)- _ FBPIK.Y)
GO TO 99
11TMIM,J)=QU)*F(I,X)- FP(ItX)
TIML, J)=0t [ )*FB( I ,X)- FBPtl.X)
JQB=JB-10
GC TOI 15, 16) , UK
_ 15 TM(M, J ) =V01_IK)*RK(K)
GO TO 99
.16 TMIM, J )=VCL (K) *AK(K)
GC TO 99
...12 Tf ( IM, J) =G( I)*F ( I ,X)-FP( I iX)
TIJ(L,J)=0(I *FB(I,X) - FBP(I.X)
TM(KM,J)=C(K)*FIK,Y)- FP(K,Y)
TM(M,J)=Q(K)*F8(K,Y)- FBP(K.Y)
00 T0( 17, 18) , UK
17 TM( N,J)=VOL(KK)*RK(KK)
GC TO 99
18 Tf( N, J)=VOL(KK)*AK(KK)
_ GO TO 99
13 TPUMtJ)= FP( I,X)-OII )*F( IfX)
TM(L,J)= FBP1 I ,X)-Q( I )*FB( I,X)
JB3=JB-12
GO T0( 19,20) f UK
19 TMIM, J)=Q(K)+VOL(K)*(TIDE11KH-RK(K) )
GO TO 99 "
20 TMtM, J)=Q(K) + VOLU)*(TID61(K)+AK(K) )
GC TO 99
14 TMIIM.J) = FPII,X)-QU) *F(I,X)
TM(L,J)= FBPII.Xl-QU )*FB( I ,X)
TM(KK,J)= FP(K,Y)-0(K)*F (K,Y)
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS'
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
0*5?
053
054
055
056
057
058
059
060
061
062
063
064
065
066
067
068
069
070
071
072
073
074
075
076
077
078
079
080
081
082
083
084
085
086
087
088
089
090
091
092
093
094
095
096
097
098
099
100
101
102
103
-------
0099
0100
0101
0102
0103
0104
0105
0106
0107
IV-C--LEVEL
21
22
99
199
2O _... MATS - - DATE -=72258 1^
GO TO(21,22),IJK
TM( N,J) =-QIKK)+VOL(KK)*(TlDEl(KK)+RK(KK) )
GO TO 99
TM( N,J) = 0(KK)*VQL(KK)*(TIDE1(KK)+AK(KK) ) . . ..
CONTINUE
_cnMTIK)iie
RETURN
END - ...
/03V ^.1
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
MATS
PAGE OOO3
104
105
106 . . .
107
108
109
1 ]O -
111
1 1 ?
-------
NOSECT DATE- .=-72258- . J.WO3/41 PAGE 0001
-FORTRAN--IV-G-L-EVEL - 2O
—OOOl FUNCTION - NOSECT-ttL,-Ii2J NOSECOOO
0002 IFUL-IS2) 1,1,2 NOSEC001
—0003 1 NOSECT=2*IL . . . . . ._. . NOSEC002
0004 RETURN NOSEC003
.1.0005- -2 NCSECT = 2*IS2*(Il.-IS2) .. - - NOSEC004
0006 RETURN NOSEC005
-.-OOO-T END ^ NOSEC006
-------
FOR TRAM
0001
0002
0003
0004
0005
0006
nnr>7 -
0008
: O009._..
0010
0011 -_
0012
nni -\
0014
.._. 0015. ._
0016
0017 ...
0018
noi9
-I-V— G-tEVEL 20 NTPS - .....DATE = 72258 14/03/41
FUNCTION NTPS (I) NTPS
REAL*8 IOJOB, IOATE, 10, IDJ
1WU(100),IMAX,IS2,ISB(100),XSB(100),RK(100),NJD(100),XJD(100),
3V(100),SD(100)Q ) ,w( 100) tVD(lOO) ,T (100) ,H( 100) ,VOL( 100) , T1DE1 ( 100) ,PR ' 100)
5PDEL(100),JBV(100),00(100),DO(100),BDO(100),
-: 6IO( 100), IDJ,FAC(100),NSV(100),U(100),IDATE,ZZ(100),XSC(100),JBM
IF( I) 1,1,2
1 NTPS =5
RETURN
2 I F ( i-j<;?) 4,4r7
4 IFIKF01 1 1-1) 5,6,5
5 NTPS =6
RETURN
6 NTPS =1
RETURN
7 1 F( T-..I MAX > 8 . 8 , S
8 K=KFO( I )
1FJK.GT. 1.ANO.K.LE.4) GO TO 9
GO TO 5
9 NTPS =K
RETURN
f_Nn .__._. _ __.
NTPS
NTPS
NTPS
NTPS
NTPS
,NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
NTPS
000.
001
002
003
004
005
006
007
008
009
010
Oil
012
013
014
015
016.
017
018
019
020
021
022
023
024_
PAGE 0001
-------
FORTRAN- I V-C—tEVEL.- 20
OM
.DATE—=- 72258-
1.4/03/41
_PAGE_ 0001
Qnni
0002
0003
0004
0005
0006 1
onn?
0008 3
0009 .
0010 2
00 11
0012 4
0013 5
0014
nnis
FUNCTION QM(KFO t B^tCHf Pt F6t I t X)
DIMENSION KFO(l) tBG(l) tCHtl)
EXTERNAL F»FB
M=KFO( I )
GC TO ( 1 »2» 3.4 ) t M
OM=BG< I )*F( I ,X)+CH( I )*FB( I ,X)
Gn Tn s
OM= CH«I)*FB(I,X)
GG TO 5
Of = BG( I )*F( If X)
GC TO 5
OK= BG(I)
CONTINUE
RETURN
END
QM
OM
OM
OM
OM
OM
OM
OM
OM
OK
OM
OM
OM
OM
OM
000 .. ..- ...
001
002
003
004
005 •
006 . . .. _._ .
007
008
009
OLO
Oil
012
013
014
-------
0001
0002
0003
0004
.
0005
0006
^007
0008
.... 0009
0010
001 1
0012 '
0013
0014
0016
nni?
0018
0020
0021
0022
002^
0024
fin?*;
0026
0028
0029
0030
firm
0032
0034
0036
0038
0039
0040
0<">M
0042
»-*— 0— «.C VC1. f
-------
CO"!
0002
OOO3
0004
0005
0006
nnny
0008
0009—
0010
0011
0012
0013
0014
1
2
3.
4
PLACE -
OATE-«--72258 —
14/03/41
PAGE 0001
-SUBROUTINE- PLACEUI, JJ-, IS2M*-I-&A.»JJLNE1-
DIMENSION LINE<41),ISA!1)
IS2 = IS2M/2 -
KS=ISA
-------
0001
0002
0003
0004
0005
0006
nno?
0008
0009
0010
00 I L
0012
nnn
0014
0015
0016
0017
0018
0019
0020
no?]
0022
nn;>-»
0024
0025
nr>26
0027
0028
0029
0030
e vet £u KKI - -- UAit— -r^^ao. - - i«»/u3/
—SUBROUTINE PR I (MI , LJ ,NK)
REAL*8 IDJOB, IDATE, 10, IDJ
rnMMnw/HY/i^AMnrM y^AMnni T n i OR A k f i nn* F P i i nn ^ n/inn^ Rnfinni
lWU(100),IMAX,IS2,ISB(100),XS(mOO),RK<100),NJD(100),XJD(100),
-.- 2SUCO),JMAX,IS1, I SCI 100), OKI 100) ,NJU( 100 ) , X JU ( 1 00 ) ,KFO(100),
3V(100),SD(100),B(100) ,G ( 100 ), AREA ( 100) , TEMPI 100) , El 100),HT( 100),
4JNQ(100),W(10G),VO(100),C(100),H(100),VOL(100),TIDE1(100),PR(100)
5PDEL(100),JBV(100),Qn(100),DD(100),BOO(100),
6IO(100),ICJ,FAC(100),NSV(100),U(100),IDATE,ZZ(100),XSC(100),JBM
COMMON/NEW/FAC1,FAC2,FAC3,ITYPE(35)
DIMENSION L J ( 1 ) t MI ( 1 )
PRINT 10
OH 1 1= 1, IMAX
10 FORMAT! ///5X7HSECT I ON , 4X , 4H KFO, 2X, ' JUNC T ION ID AND DISTANCE'
112X,2HRK,8X2HAK,8X2HnK,9XlHE,8X2HFF,6X,4HTEMP3X, 1PDELI/28X,'AT1,
2' ENDS' ,20X, 'REM. RATE ',« REA.RATE ',' DEO. RATE ',' DISP.RATE ',
3' ULT/5DAY '.AX,' (•; I LE- I NCR ' /56X , ' ( 1 /DAY) ' , 3x , ' ( 1 /DAY ) ' , 3X ,
4l{l/DAY)t,2X,t(MI2/DAY>',HX,1(CENT)1,2X,MMI)1/)
1 PRINT 20 , I , ID 1 I ) , MI 11 ) ,KFQ( I ) , NJDt I ) ,XJD( I ) , NJU ( I ) ,
1XJU( I ),RK( I),AK( I) ,OK( I ),E( I ), FF(I ), TEMPI I ),PDEL( I)
__20 FGRMAT(1HO,I2,1XA4,1XA4,3X,I<,,1X,2UX,1H(,A3,1H,F6.1,1H),2X),
16F10.3,F7.3)
I FINK 131 ,31,32
31 PRINT 30
_30 FORMAT (///5X7HSECTION,6X5HCROSS,10X1HO,13X2HBD,12X2HPR,10X5HTIDE1
11 lX,2HHT,ax,3HVOL, 12X,'WU'/17X, 'SEC-ARE A', 6X, • FLOW, 7X,'BENTHL-DM
_ _ 2' ,5X'PHOTC-RESP' .^X'TIDE-COEF'^X'DEPTH' ,5X' VOLUME1 ,8X, 'UNF-LOAD'
GO TO 34
..._32 PRINT 33
33 FORMAT! ///5X7HSECT I ON, 6X5HCROSS, 10X 1 HQ, 1 3X2HBD, 1 2X2HPR ,
llCX5HTIf)El,llX,2HHT, 8X , 3HVOL , 8X , ' WU/ AREA*RK'/17X,' SEC- AREA ',6X,
2 'FLOW t 7X, '6ENTHL-DMD1 ,5X, ' PHOTO-RESP1 ,4X, ' T IDE-COEF
.L 3, 7X, 'DEPTH1, 5X,' VOLUME ',8X, 'UNF-10AD'
34 IF(NK)5,5,7
5 PRINT 6
6 rORMAT(17X,' ( FT**2 ) • , 7X , ' ( CFS ) • , 5X , • ( GMS/M2-DAY ) • ,
J.4X, ' (MG/L-DAY) 1,20X,1(FT.)',4Xt'(FT**3)l,6X,' (LBS/ DAY-MI ) •
GO TO 9
7 PRINT fl
8 FnRMAT<17X,MMI**2)« ,4X, • ( MI**3/DAY) • ,2X,
l« (GMS/M2-DAY) ' ,4X, ' (MG/L-DAY) • ,20X, • (MTRS) ' t4X, • (MI**3) ' ,8X,
2' (MG/LP
q nn 2 I -1 , IMAX
2 PRINT 40, I, ID( I ) ,MI I I ) ,AREA( I ) ,Q( I ) ,BDI I ) ,PR( I ) ,TIOE1( I) ,
1HTI I ) , VOL ( I ) , WUI I )
RETURN
40 FORMAT! 1HO 12, 1XA4, lXA4,2E12.3,5X»E12.3,2XtE12.3tlXtE12.3flXt
12E12.3.3X.E12.3)
tNTRY PRII (HItLJtNK)
PRINT 50
DO 3 1=1, IMAX
50 rORMAT(///3X7HSECTION,8X,lHU,9XlHS,9XlHVf8X2HSD,
"*i
PR I
PRI
PR I
PRI
PRI
PRI
,PRI
PRI
PRI
PRI
PRI
PRI
PRI
,PRI
PRI
PRI
PRI
PRI
PR I
PRI
PRI
PRI
PRI
PRI
,PRI
DPRI
PRI
PRI
PRI
PRI
PRI
' PRI
PRI
PRI
PRI
PRI
_PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
non
001
002
003
004
00.5 .
006
007
008
009
010
Oil
012
013
014
015
016
017
01.8
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
-------
0031
0032
0033
0034
on^s
0036
0037
0038
00 VJ
C040
0041
0042
0043
0044
0045
0046
. 0047
lflX?HVri,<'XlHB,c>Xlnc,<3xiHG,c)XlHH,8X2H£/c>X'*HFAC/)
3 PRINT 60,I,ID( I), MI ( [) ,U( I), SI I), VI I),SD(I),VD{ I) ,B(I),C(I) ,
1GII ),H(I),ZZ( I)iFAC( I)
RETURN
60 FORMAT! 1HO,I2,1XA4,1XA4,9F10.5,2F8.2)
ENTRY PRJ ' , 2X, • ( M I **3/DAY )
li 4X, • (MG/L) ' ,5X, • IMG/L) '
16 DO 4 J=l.JKAX
4 PRINIT80,J,JBV(J),LJ(J),JNO(J),ISA(J),XSA(J),ISB(J),
lXSB(J),ISC(J)fXSC(J),W(J),OD(J),DD(J),BDD(J)
80 FCRMAT(1HO,I2,I3,1XA4,1XA3,1X,3(1H( , I 3, F5. 1, 1H ) ) , 6XF12. 5, 3X,
13F1P.51
RETURN
END
PRI
PRI
PRI
PRI
PRI
PRI
PR!
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
1 PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
PRI
052
053
054
055
056
057
C58
059
060
061
062
063
064
065
066
067
06fl
069
070
071
072
073
074
075
076
PAGE 0002
-------
PORTRAN-IV-G-L-EVEL- 20
SVO
.DATE
72258
PAGE 0001
onoi
0002
noo?
000*.
0005
0006
nnn?
0008
0009
0010
___ 0011
SUBROUTINE SVOIU.E, A.S, V)
IFfU) 2,1,2
1 S=SQRT (A/E)
V=-S
RETURN
2 SC=SORT (l.+4;*A*E/U**2)
S = P*U. + SQ)
. _V = P*(1.-SC).
RETURN
...__£ND_
SVD
SVD
SVD
SVD
SVD
SVD
SVD
SVD
SVD
SVD
SVD
000
001
002
003
004
005
006
007
008
009
010
-------
-J=ORTRAN_Ut_G—L-EVEL-- 20
SVD
0002
0003
0004
0005—.-
0006
0007
1
2
..SUBROUTINE. SVOtU.E. A,S, V).
IF(U) 2,1,2
.1 S = SORT (A/E)
V=-S
RETURN
SC = SORT (l.+4.-*A*E/U**2)
P. = .5*U/E ...... . __________ ______
0008
0009
0010
0011
__________ V=P*(1.-SC).
RETURN
_____________ END.
.DATE
72258
14/03/41
PAGE 0001
SVD
SVD
SVD
SVD
SVD
SVD
.SVD
SVD
SVD
SVO
SVD
000
001
002
003
004
005
006
007
008
009
010
-------
FORTRAN_lV_G-tEVEL 20
YB
DATE = 72258
14/03/41
PAGE 0001
nnm
0002
nnoi
0004
0005
0006
non?
0008
0009
0010
FDNf.T TON YR( 1 ,X)
REALMS IDJOB, IDATE, ID, IDJ
COMMON/HY/ISA(100),XSA(100),IDJOB,AK(100),FF(100),Q(100).BD(100)f
1WU( 1 00), I MAX, I 52, I SB (100) ,XSB(100) , RK ( 100) , N JD( 100) ,XJD(100),
2S (ICO) tJMAX, I SI t I SC ( 100) iDKI 100) ,NJU( 100) tXJUl 100) fKFO( 100) ,
3V < 1 00), SO ( 100) ,B(100),G( 1 00), AREA ( 100 ), TEMP ( 100 )• E ( 100 ) ,HT ( 100 ) ,
4JNO{100),W(100),VD(100),C(100).H(100)tVOL(100),TIDEl{100),PR<100)
5PDEL(100),JBV(100),OD(100),DDI100),BOD(100),
6IC(100),IDJtFAC(100)tNSV(100)tU(100)iIDATE,Z2(100),XSC(100)fJBM
IFIRM I ) )2,2il
1 YB = WU( I)/(AREA(I )*RK( I) )
YBR=YB
r,n TO i
2 YB=0.
3 RETURN
END
YB
YB
YB
YB
YB
YB
.YB
YB
YB
YB
YB
YB
YB
YB
YB
YB
000
001
002
003
004
005
006
007
008
009
010
Oil
012
013
014
015
-------
-FORTRAN-IV-G-LEVEL- 20
YD
DATE = 72258
1W03M1
PAGE 0001
-00 Oi.
0002
..0003
0004
0005 _..
0006
__ 0.0.03
-.FUNCTION YOU.X) _ .. _
REAL*8 IDJOB.IDATE.ID.IDJ
CQMMON/HY/ISA(100),XSA(lOO),IDJOB,AK(lOO),FF(100),0(lOO»,BD(100)i
IWU(IOO), IMAX,1S2,ISB(100),XSO(100»,RK(100),NJD(100),XJDf100),
-2S(100),JMAX,ISl,I SCI 100),DKt100),NJU(100),XJU(100),KFO(100),
3V(100),SO(100I,B(100),G(100),AREA!100).TEMP(100),E(100),HT(100),
_4JNO(100),W1100),VD(100),C(100),HI 100),VOL(100),TIDE1(100),PR(100]
5PDEL(100),JBV(100),QO(100),DD(100),BOD(100),
610(100),ICJ.FACf100),NSV(100),U(100), IDATE.ZZI1001,XSC(100),JBM
EXTERNAL FR.FBR
...YC = FAC(I)*OM(KFO,B,C,FR,FBR,I,X) + ^ Z 1 I )
RETURN
END _ __ _ ..._..
YD
YD
YD
YD
YD
YD
YD
YD
YD
YD
YD
YD
YD
000
001
002
003
004
DO'S
006
007
008
009
010
Oil
012
-------
-FORIRAN-IV-&-LEVEL- 20
YDP
DATE = 72258
1W03M1
PAGE 0001
0002
_OOQ3_
0004
...0005..
0006
__OOJ17_.
_£UNCTION YDP(ItX) . ... YDP 000
REAL*8 IDJOB.IOATEfIDiIOJ YDP 001
.. CDMMQN/HY/1SA(100),XSAt100),IDJOB,AK(100),FF(100),0(100!,RD<100), YDP 002
1WU(100),IM.AX,IS2,ISB(100),XSB(100),RK(100),NJD(100),XJD(100>, YDP 003
.2S(100>,JMAX,ISl,lSCliOO),DKUOO),NJU(100),XJU(100) ,KFO(100) , YDP 004
3V(100),SD(10p),B(100),G(100),AREA(100),TEMP(100),E(100)fHT(100), YDP 005
_4JNOtlOO)th'(lbo)tVO(100)iC(100),H(100) i VOL ( 100 ) , T I DEI ( 100) , PR (.100 ), YDP . 006
5PDEL(100),JHV(100),GDI 100),DDl100).BDDIIOO), YDP 007
6IDUOO),IDJ,FAC(100),NSV(100),U(100),IDATE,ZZ1100),XSC<100>,JBM YDP 008
EXTERNAL FPRtFBP.R YDP 009
. YDP = FACt I)*OM(KFO,B,C,FPR,FBPR,I,X) YDP 010
RETURN YDP Oil
_END _ _ _ YDP 012
------- |