UNITED STATES ENVIRONMENTAL PROTECTION AGENCY
REGION II
26 FEDERAL PLAZA
NEW YORK. NEW YORK 10007
documentation for
HAR03
"A computer program for the modeling of water quality
parameters in steady state multi-dimensional
natural aquatic systems"
SECOND EDITION
prepared by
Steven Chapra
George A. Nossa
Systems Analysis Section
Data Systems Branch
October, 1974
-------
TABLE OF CONTENTS
ABSTRACT
INTRODUCTION
THEORY
One dimensional analysis of a single substance T—2
One dimensional analysis of coupled substances T-10
Multi-dimensional analysis . T-13
Application of the theory via a computer program T-15
THE COMPUTER PROGRAM
Flow chart of the main program C-3
Description of the main program and the subroutines C-4
Listing of the source deck C-12
USERS MANUAL
Description of input U-4
Example problem U-15
Job Control Language U-46
Restrictions U-47
APPENDIX
References
List of figures and tables
Nomenclature
ACKNOWLEDGEMENTS
-------
ABSTRACT
HAR03 is a computer program which can be used
to model the steady-state distribution of water
quality variables for multi-dimensional bodies
of water. The technique underlying the program
is based on the conservation of mass and up to
two variables reacting in a feed forward fashion
with first order kinetics may be modeled.
HAR03 has been designed for an IBM 370 computer,
requires approximately 184 IC of core to compile
and was written for a Fortran IV G or H Level
compiler.
HAR03 has been designed with the Biochemical
Oxygen Demand (BOD) - Dissolved Oxygen system
in mind.
With minor modification the program may be used
to model other variables which are analogous to
this system such as chlorides, polyphosphate -
orthophosphate, coliform bacteria, etc.
In order to save core for small systems, HAR03
has been compiled in three versions. The first
version handles a system of up to 50 segments
and is designated HAR50; similarly HAR100 and
HAR200 handle a maximum of 100 and 200 segments
correspondingly. The only difference between
these versions is in the size of the arrays
defined in the programs.
ii
-------
INTRODUCTION
-------
APPENDIX 1 - LISTING
C
C THIS PROGRAM CONVERTS SYSTEMS PARAMETERS OF HAR03 VERSION 1 (APRIL 197:
C TO CONFORM WITH INPUT REQUIREMENTS OF HAR03 VERSION 2.0 (SEPT. 1974)
C FOR AREA, DISPERSION, FLOW AND INTERFACE NUMBER
C
DEFINE FILE 5(400,42,U,KK)
DIMENSION AREA(6),E(6]iQ{6),IARAY(6)
KK = 1
2 REAC!12,1XAREA( J) ,E(J) ,Q( J) , IARAY1 J) , J = l,6)
1 FORMAT ( 3 ( 3F6 . I i 12).)
IF(AREA(1) )10tll»I0
11 KSW=:KSW + 1
IF(KSWr2)I2t59,99
10 KSW = 0
12 WRITE 15 • KK )(AREA(J)» E ( J)fQ(J)*IARAY(J)* J=1,6)
GO TO 2
9 9 KXK=KK-1
KK= 1
DO 15 1= 1 »KXK
READ(5•KK)(AREA(J)»E(J),Q(J),IARAY(J),J=1,6)
WRITE (2, 3) (AREA't J)«E( J) *Q{ J)f IARAY t J ) ,J = 1,6)
3 FORMAT(3(2F6 • 21F6•Ct I 3) )
15 CONTINUE
CALL EXIT
END
A- 1
-------
The formulation of a mathematical model of any system is
greatly determined by two factors: the nature of the system itself
and the purposes and perspective of the investigator. The modeler
must strike a balance between objective reality and the subjectivity
of his needs to attain a successful analysis. This problem is further
compounded, when dealing with the high complexity of the natural world.
One of the more prevalent misconceptions among neophytes in the
field of water quality modeling is that there is one analytical technique
which is superior in depicting the water quality in a natural body of
water. This may be partially due to the fact that the field straddles
several more or less hard sciences and engineering disciplines and as
such can be perceived from a variety of perspectives.
For instance, hydrodynamicists, who are essentially interested in
the movement of fluids, often tend to emphasize the obviously important
effect of water motion on the transport of matter in a system. Ecologists
and aquatic biologists on the other hand stress the equally important
reactions between the community of organisms which populate the system.
The danger in these or in any particular approach comes from the automatic
exclusion or underestimation of viewpoints outside the area of expertise of
the modeler.
One of the older approaches to water quality modeling which rather
effectively incorporates a number of perspectives in representing the
causal relationships of stream pollution is that of the sanitary engineering
profession. Due to their interest in designing waste treatment facilities,
sanitary engineers were rather early introduced to the problems of wastes
and their impact on the environment. A classic study in this profession
was that done by Streeter and Phelps^- on the Ohio River in 1925.
By making a variety of simplifying assumptions in the hydrodynamic
and biological areas, these investigators arrived at a very utilitarian
approach to water quality analysis which still stands as a viable technique
for answering many questions about the relationship between pollution and
the aquatic environment of a stream. In the hydrodynamic area, they
assumed that the waste load was delivered by a pipe into a channel which
could be described as having constant geometrical dimensions and constant
flow. As well it was assumed that the pollutant was instantaneously mixed
in the lateral and vertical directions and that the simple continuity
equation, Q=AV, applied. From the biological standpoint, it was decided
that a chemical parameter upon which most species depend for life, namely,
dissolved oxygen could be modeled as an indicator of the health of the
biota. To do this, they had to use a measure of the oxygen demand of the
waste, the biochemical oxygen demand (BOD), as the input to the system
and formulated relationships between dissolved oxygen and BOD in terras
I- 2
-------
of first order kinetics. The result is what is now called the
Streeter-Phelps equation which in its basic form is:
K
„ r d / -K x/u -K x/u \ *1 T . x/u /T
D = H; — ( e r - e a )]L + Dea' (1-1)
K - K o o
a r
where:
D=dissolved oxygen deficit= D0S -DO
DO =saturation concentration of dissolved oxygen
DO =actual concentration of dissolved oxygen
Lo=initial concentration of BOD at point of introduction of waste
Do=initial concentration of dissolved oxygen deficit at point of
introduction of waste
Kr=BOD removal rate
Xd=Deoxygenation rate
Ka=Reaeration rate
U=stream velocity
x=distance downstream from point of introduction of waste
The results of the Streeter-Phelps equation is called the "D.O.
sag" and is illustrated in figure 1-1.
Figure 1-1: D.O. sag generated by
Streeter Phelps equation
-------
This environmental model is ideal for the. evaluation of various
treatment schemes as its basic control "variable is the waste input.
This emphasis on relating man's waste inputs to the aquatic environment
with the express purpose of managing the inputs and thus the water
quality is what typifies the sanitary engineering approach. This can
be contrasted with an aquatic, biologist who might be more interested in
the interraction batween the organisms with, a mind to prediction and
description rather than control.
Since 1925 various contributions and technical developments have
advanced the state of the art of water quality modeling. Particularly,
the advent of the digital computer has expanded study beyond the limits
experienced by Streeter and Phelps. One of the major influences from the
modern sanitary engineering perspective is Robert Thomann, formerly
technical director of the Delaware Estuary Comprehensive Study (DECS)
and now a professor at Manhattan College. The program documented here,
HAR03, is basically a representation of an approach developed by
Thomatm for the Delaware^ which has proven itself to be particularly
affective as a tool for water quality management.
As with Streeter and Phelps, the body of water is assumed to be at
a steady state in time, and BOD and dissolved oxygen are the water
quality variables (however, the program is general enough so that other
variables such as temper attire, coliform bacteria, etc.. which are
analagous to BOD and dissolved oxygen, could be modeled). Then the
system is broken up into finite volumes which are completely mixed.
By stipulating a priori the flow and dispersion across the interfaces
between volumes, the hydrodynamics of the system are input to the model
and by using macbeaiatics the concentrations of the variables in each
volumn can be calculated. Thus the approach allows extension from streams
to estuaries and lakes.
This report is in several parts. First, the theoretical background
is developed. This is done in terms of an estuary since it is one of
the more complex bodies of vat'er to which this program can be applied,
flowever, the concepts expressed can be readily applied to lakes, streams
and othgr bodies of water. Then, a description of the computer program
is given followed by a user's guide,
As a final note, KAR03 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 must 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 aquatic system or problem. 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.
1-4
-------
THEORY
T—1
-------
One dimensional analysis of a sipgle substance...
¦I
Thomanri's approach ±f to steady-state estuarine. water quality
modeling essentially applies a numerical solution technique to a
convective-diffusion equation for mass transport including decay
and source terms. This equation, in time variable form for a one-
dimensional channel of uniform cross-section, may. be written as:
3c 3c 92c
— + U— = E — lcc + (sources — s?Lnks) (T-l)
3t 9x 9x2
where
x=distance
t=time (on the scale of minutes to hours)
c=concentration of a constituent
E=longitudinal dispersion coefficient
U=tidal velocity
k=first order decay rate of constituent
A steady state in an estuary is a somewhat difficult concept to
define. The term literally means that the constituent under inves-
tigation and the parameters which describe it do not vary in the
time scale of the problem. As an estuary, in this report, is defined
as a body of water subject to tidal oscillations, it can be argued
that such a system will never reach a steady state but will con-
tinuously be in transition. Several investigators 3,4,5 have
treated the subject and one approach is to take average values over
a tidal cycle. The following relationship results:
9c 3c —3zc — . . .
— + Uf— = 12 - kc + (sources - sinks) (T-2)
3tt 9x 9x2
where
the bar above a parameter or variable designates
that the value is averaged over a tidal cycle
U.£ = net advective velocity. In many cases this can be
interpreted as the velocity caused by fresh water flow
entering the estuary
tt = time (on the scale of the time of a tidal cycle -
several hours)
E = a tidally averaged dispersion coefficient which includes
the dispersive effect caused by tidal motion
T-2
-------
A "steady state" for equation T-2 therefore means that.the
variables and parameters in question do not vary from tidal cycle
to tidal cycle. This situation is often approximated in the
summer and early fall when low rainfall minimizes the transients
caused by fresh water runoff to an estuary and meteorological
conditions are relatively constant. These periods are often
critical when investigating such a classical gauge of pollution as
dissolved oxygen.
With this as background, Thomann suggested that the system be
divided into volumes or segments as in the one dimensional system
depicted in figure T-l.
k-1
k+1
1
2
3
k
n—1
n
Figure T-l: Estaaify of n segments
Equation T-2 can then be written in the following form for
Section k (from this point on the bar which was used to designate that
a parameter or variable was tidally average will be dropped for convenience)
Vol
k3t.
8ck
+ Q„ 1, = E, A, 1
9 c,
fk\^- " Wkrr- ~ V01^ wk
(T-3)
where
k= a subscript which denotes that the variable or parameter is
an average for the segment k.
Q^k= average net advective flow for the segment.
=Ufk \
Vol^= volume of the segment = A^l^
A^= cross-sectional area of the segment
1^= length of the segment
wk= waste input to the segment
T-3
-------
At steady state c^./ t would equal zero and equation T-3
could be written as:
0 " WkT " " V°Vk + Mk (I-4)
dxz dx
Each term containing a differential in equation T-4 can be
approximated mathematically, as:
d2c'k
EkAk1k~~7" = Ek-1 ,k (ck-l~ck) + Ek*k+1 (ck+l~ck> ^T~5^
dx^
dck
^fk1lc = Qk,k+l(ak,k+lck+ek,k+lCk+l) " Q
dx
k-l,k(ak-l,kck-l +ek-l,kck)* * *(T 6)
a
where
E' . = E. .A. . /I (T-6a)
ij i.J i»J i,j
i,j = subscript designating the interface between segments
i and j.
1. • = The average length of adjacent segments
5 J
= )li + lj)/2
- • and B- • = weighting coefficients used to correct the
1'-' 1'-' approximations of equation T-6 for adjacent
segments of unequal length.
1.
3 .
l.+l.
i J
(T-7)
[-4
-------
*i,.1 -
1»J
(T—8)
Xi+1j
As well, the above weighting coefficients
can be used to insure meaningful results. It
can be shown that for all positive solutions
o?>l-E'/Q (T-8a)
It must be stated here that Thomann's
approximation in equation T—5 is strictly
correct for segments of equal length (and in
two and three dimensions for orthogonally shaped
segments). Therefore, care should be exercised
to avoid the placement of segments of highly
different lengths adjacent to each other or the
construction of very irregularly shaped segments.
In the event that this is unavoidable, adjustment
of E may be necessary to reflect the true effect
represented by equation T-5.
As well, it has been shown that a numerical
dispersion effect is implicit in this particular
difference scheme which.can be estimated by
¦Wrical " U1 Co " 1/2> W-8b'
This can be particularly significant in streams
where velocities are high and distortion may result.
It is therefore imperative that the modeller be always
aware of the numerical dispersion in this model and
make appropriate correction for it.
T-5
-------
Equations T-5 and T-6 can be substituted into (T-4):
0 = E k-l,k (ck-l"ck) + E k,k+l (ck+l"ck)
+ Qk-l,k (°W,kck-l +
" Qk,k+1 (otk,k+lck 6k,k+l ck+i>
" KkVolkCk + Wk • ¦ • • V
By grouping terms
("Qk-lskak-l,k"E k-l,k) ck-l
+ (Qk,k+laksk+l " Qk-l.k^k-l.k +- E k-l,k + E'lc,k+1 + VolkKk)ck
+ (Qk,k+].6k3k+l " E k,k+l> ck+l = Wk * (T—10)
Letting
ak,k-l = ~Qk-l,kak-l,k ~ E k-l,k (T—11)
ak,k = Qk,k+lak,k+l_Qk-l,kek-l,k + E k-l,k + E k,k+l + VolkKk ¦ •(T-l2)
ak,k+l = Qk,k+l3k,k+l ~E k,k+l ¦ • (T-13)
T-6
-------
The general equation for the kth segment is therefore
ak,k-l ck-l + ak, k ck + ak,k+lck+l = wk (T-14)
For the furthest upstream section (figure T-l, section 1), which
is a boundary of the system, a similar procedure would yield:
all C1 + a12 c2 = W1 + (QOl a01 + E Ol) . . (T 15)
or
I
a-^^ c^ + a^2c2 ~ ^ 1 • • • • 16)
1
where
o = a subscript to denote the values of the variables and parameters
beyond the upstream boundary.
For the furthest downstream segment (in figure T-l, segment k):
Si-l.n^-l + an,ncn Wn + ^-(^n,n+l^n,n+l + E n,n+l^cn+l • • • • (T-17)
or
an-l,ncn-l + an,ncn = w n • (T-18)
where
n + 1 = a subscript to denote the values of the variables and parameters
beyond the downstream boundary.
T-7
-------
The complete set of equations for the system can be written
as
a^c^ + ai2c2 ^ 0 +
a9jc^ + a22c2 a23c3 + 0 +
° + a32°2 + a33C3 + a34C4 +
+ 0
+
+
0
0
w.
w„
w„
0 0
0 + . . „ + 0 + a ,,c + a c = Wn
n,n•1 n+1 nn n
(T-19)
T-8
-------
At this point th.ere are n equations with 11 unknowns. Solutions
may be obtained either by solving the equations simultaneously or by
writing T-19 in matrix form and inverting it. This latter technique
is employed in HAR03. Equation T—19 in matrix form is
ail a12
a21 a22
a23 0
a32 a33 a34
c_ !,
0 | ' c
2
8 !
0 i ! c„
W,
w
w.
¦\
ian,n-l
n i
W
n i
(T-20)
or
[A] (c) = (W) . .
By inverting [A] a solution vector can be obtained:
(T-21)
(c) = [A]-1 (W)
(T-22)
T-9
-------
Besides allowing solution of equation (T-21)^ equation T-22 offers
an additional feature. The elements of the matrix IA]-1, which is called
the single system response matrix, represent the unit response of a particular
section ta the -addition of a unit load to another section. For instance,
the element a.2? of is the unit concentration response of section 2 to
a unit input of waste into section 3. This is a particularly useful
feature of Thomann's technique in that once a system response matrix is
generated for a particular set of parameters, the effect of various
load changes reflecting treatment strategies can be directly determined
from this matrix. It, therefore, provides a quick technique for judging
the trade-offs of various treatment alternatives.
One dimensional analysis of coupled substances...
Wv
' P
E3
~-k.
t>-L
Figure _-2: 'Schematic for a coupled
system of cao reactants
Figure T-2 is a schematic for a coupled system of two constituents,
c and b, which react with first order kinetics. Wc and Wb are waste
inputs of the constituents and the k's are the reaction rates which link
them.
In classical sanitary engineering practice the BOD-dissolved oxygen
system corresponds to the schematic with c = BOD and b = dissolved
oxygen deficit. The rate is referred to as the BOX' removal rate*
is a rate of reaeration from the atmosphere and k^ is a deoxygenation
rate. In the case where there is no loss o£ BOD via sedimentation or
another non-oxygen demanding removal, k-^ would equal .
The equations describing the conservation of ntass for c would be
exactly like that developed in the previous section. The equation for
b in a segment k at steady state would be
0 = E' (b ~b ) + E ("b -h )
k-1, k k-1 k k,k+l k+1 k
+• 0, , ¦ (a, , , b + 8 ' b )
k-1„k k-l,k k-1 k-l,k k
0 (a b + & b )
k,k+l k,k+l k k,l{+2 k+1
Vol X fc. + Vol K. + W
k ir;2 k k k fc
(T-23)
T-10
-------
T-21:
where
where
= BOD in section k
b^ = dissolved oxygen deficit in section k
IC^2 = reaeration rate in segment k
= deoxygenation rate in segment k
= waste load or source of dissolved oxygen deficit to
section k
By grouping terms this results in a matrix equation similar to
[B] (b) = [Vol K3] (c) + (Wb) (T-24)
[Vol K3] = an nxn diagonal matrix
[B] = similar in form to•[A] with the
exception that the diagonal elements
contain the expression Vol father than
Vol Kkl
(c) = solution vector to the single system (in this
case BOD) „
The solution is therefore
(b) = [B]-1[Vol K3](c) + lB]-l(Wb) (T-25)
In this program, three sources of dissolved oxygen deficit will be
considered:
Wv = -PMR + BD*H + Y (T-26)
T-ll
-------
where
PMR = the dissolved oxygen deficit due to the photosynthesis
and respiration of algae.
BD*H = the dissolved oxygen deficit due to bottom deposits.
Y - direct oxygen deficit from waste sources
One additional feature of the analysis should be noted. By
substituting equation T-22 into equation T-25 the following rela-
tionship results
-------
For boundaries where flow enters the section with a concentration c :
b
3ii = I (Qikaik+E'ik)+ViKi+Qii6ii+E,ii (T"33)
Wi = Wi +
For boundaries where flow leaves the- section to an area with a
concentration c, :
D
aii- = * (Qikaik + E'ik> + Vi + Qiiaii + E'ii
W- = wi + (E ii - Qiiaii) c
(T-36)
Similar to the previous analysis for one dimension, equation T-30
together with the appropriate boundary conditions can be incorporated
into a matrix.
lll a12 a13
i2i a22 a23
In
l2n
cl!
nl n2
or
nn
•- c
>
i
(c
! n
2|
. wn '
E 21
i
! }
I \
i w
i n
\ i
\ J:
[A] (©) = (W) (T-3'7)
with a solution
(c)=[A]"l(W) (T-38)
This may be extended to a coupled system of two reactants (e.g. BOD-DO
deficit) in a fashion analagous to that done on page
As well, the criterion that
a > 1 - E /Q ........
to insure positive results also applies.
(T-39)
1-14
-------
Application of the, theory via a computer program...
At this point, the solution has been reduced to two equations:
T-22 and T-26b for BOD and DO deficit, respectively.
(c) = [A]-1 (Wc) (T-22)
(b) = [C]"1 (Wc) + [B]-1 (Wb) (T-26b)
For HAR03, equation T-26b will not be used due., tp storage requirements
for the additional matrix, [C].* In its place, a form of equation
T-25 will be used
(b) = [B]"l ((Vol * K3) (c) + (Wb)) (T-40)
Aside from input-output considerations, the major part of the
algorithm of HAR03 therefore reduces to the setting up of the matrices
[A] and [B] and their respective forcing functions.
Since the only difference between matrix IAj and matrix [B] is that
their diagonals contain different reaction terms, a matrix excluding
reaction terms can be formulated. Called the system matrix, [AB], it
can. be transformed into [A] of [B] by merely adding the appropriate
reaction term to the diagonal.
note: another version of the program will be available in the
future which uses T-26b and which generates the valuable
total system response matrix,
T-15
-------
With this as a background, a very simple representation of
major steps of HAR03 can be presented:
1
c
START
step 1
Read data which is applicable to the
system as a whole (e.g. hydrodynamic and
__ physical data)
' "step 2
Cons true t the^sys t em matrix, AB i r
step 3
Read data which is applicable to the
specific constituent which is being model-
ed (e.g. rates, loads, boundary
.. v.,„„^QIld^t i pn.s )
step 4
Construction of the forcing functions
step 5
Construction of the specific constituent
matrix, AC
step 6
Calculation of the concentration of the
constituent by inverting the matrix, AC and
multiplying it by the forcing function
Figure T-.4: Simple Schematic of HAR03 computer
program
T-16
-------
Of the six major parts, step 2 and step 4 require further
discussion:
Step 2 - The construction of the system matrix IAB]
When constructing [AB], each section is considered individually
with each of its interfaces treated one at a time. The following
sequence is undertaken for each interface:
1) calculation of E' as in equation T-6a
^ij i :
^ (l.+i, )
Eij Ai :
E, =-_-J (T-41)
2) a value of alpha is calculated to adjust the advectivf
approximation for segments of unequal length.
for flow across the 1^
interface and into the a
Ji 1± +1.: (T-42)
section J
(0< o)
for flow across the _ 1
interface and out a ij
J
r:r r 1±+^ o)
3) the criterion for positive solutious is applied to alpha
(equation T-39) and if it is not met, an adjustment is made.
This is accomplished by the following test:
if a>l-E'/Q do not alter
if a
-------
3) calculation of the elements of the system matrix
The elements are calculated as in equation T-31 (excluding
the reaction term) and T-32.
When each of these three steps is done for each interface, the
procedure is repeated for each segment until all have been incorporated
into the matrix.
Step 4 — The construction of the forcing function
At present, the forcing functions which are calculated are
peculiar to the BOD - DO deficit system. From the theory, equations
T-34 and T-36 give their most basic form:
for Q > o Wi' = W± + (E'ijL - Q±i & ..) cb (T-36)
for Q
-------
2) Other forcing functions
of the form,
for BOD or chlorides
W. = any load to the section (e.g. waste loads, runoff
loads, etc.)
for DO deficit
W. = -PMR + BD*H + loads
1
as explained in equation T-26.
T-19
-------
THE COMPUTER PROGRAM
C-l
-------
The computer program from which HAR03 evolved was developed by
Hydroscience, Inc. for the Massachusetts Water Resources Commission-'-.
It was originally written for an IBM 1130 computer in FORTRAN IV and
utilized a Gauss-Seidel technique to solve the set of simultaneous
equations generated by the Thomann technique. The program was never
adequately documented and as such was not appropriate for general
distribution.
A source deck of the program was made available to the Environmental
Protection Agency by the State of Massachusetts and after modification
two new versions were published by E.P.A.^ entitled HAR01 and HAR02.
Compatible with an IBM 370 and an IBM 1130, respectively, these versions
retained the Gauss-Seidel solution technique and had the advantage of
being thoroughly documented.
The major difference between HAR03 and the previous versions is
in its use of a matrix inversion solution technique. As well, parts of
the algorithm have been redesigned to more accurately reflect Thomann's
technique and to improve its clarity and ease of.use.
In its present form HAR03 is programmed for an IBM 370/158 computer.
The core requirements for each version have been summarized below:
Version Maximum No. Storage
of segments Requirements
HAR50 50 76K
HAR100 100 148K
HAR200 200 410K
The compile and link, time for any of these versions is .46 minutes,
which includes optimization by the Fortran IV H level compiler. Without
optimization, compilation and linkage takes .24 minutes. Each of these
versions are stored as load modules at OSI in Bethesda, Md. (EPA computer
facility).
An estimate of the actual execution time is given in table C-l. Each
time is for what is called an estuarine additive coupled system (described
on page U-3) which analyzes the chlorides, carbonaceous and nitrogenous
BOD and the D.O. deficit of a system. This type of run performs five
matrix inversions which in larger systems takes up the largest amount of
t ime.
Number of segments execution time (minutes)
8
.02
30
.10
45
.27
64
.60
C-2
-------
SUBROUTINE DATA
Input General System Parameters which do not depend on the
constituent being modeled
| SUBROUTINE REVI (optional)
| Used to revise the parameters E and/or Q without altering
j any other elements of the input
SUBROUTINE FLOW
Computes a flow balance around each section to insure
continuity of flow
SUBROUTINE SETUP
Primarily sets up the System Matrix, [AB]
/~
READ INDIC & IEXIT
INDIC is used to designate the type
of constitment which will be
described by the subsequent set of
specific constituent cards
IEXIT is used to terminate the run
if so desired
STOP
SUBROUTINE BCC
Primarily introduces the boundary conditions into the algorithm
SUBROUTINE RATE
Reads in rates, constructs the specific constituent system
matrix [AC] and the forcing function
SUBROUTINE WRITE (optional)
I Writes matrix [AC]
±
SUBROUTINE MATH
Inverts [AC] and multiplies it by the forcing function
I SUBROUTINE WRITE (optional)
j Writes matrix [AC]**-1
SUBROUTINE STORE
Primarily writes final results and retains output from
this particular iteration which will be used in subsequent
iterations
FIGURE C-l: Flow chart of the main program of HAR03
0-3
-------
DESCRIPTION OF THE MAIN PROGRAM AND EACH OF THE SUBROUTINES
The Main Program-HARQ3...
HAR03 has been programmed using a modular approach. Wherever
possible, an attempt has been made to give each subroutine a unique
function. The main program was designed to coordinate the subroutines
to effectively analyze the water quality of a system. This coordination
is depicted in figure C-l which represents a development of the simple
schematic given in figure T-3. In figure C—1 certain details have been
omitted and if the user requires more detail he should consult the code
itself.
Subroutine DATA...
The primary purpose of this subroutine is to input parameters
which do not depend on the particular substance being modeled. These
include system parameters (number of sections, indicators, etc.) and
certain physical parameters such as cross-sectional area, net advective
flow, etc., which are only input once when modeling a water body. These
are in contrast to parameters such as loads and reaction rates which would
vary with the constituent being modeled.
The particular tasks accomplished by DATA are:
1) read and write system parameters
2) initialize certain matrices to zero
3) read and write interface parameters AREA, E and Q and place
them into the storage matrix ARRAY.
4) read and write the characteristic lengths (LA)
5) read the general section parameters depth (H), temperature (T),
and volume (VOL)
6) Apply the following conversion factors:
Q(MGD) = Q (CFS) * .6463 MGD/cfs
VOL(MG) = VOL (106ft3) * 7.48 gal/ft3
C-4
-------
At this point, a further word about the function of ARRAY is
in order. ARRAY is an N X 18 sized matrix which is used to save space
when storing the interface parameters AREA, E and Q. Used in conjunction
with IARAY, ARRAY avoids the use of three N X N matrices when handling
these variables. It does this as follows:
IARAY(I,J) contains the section number
of the section which forms
the interface with section I.
ARRAY(I,J)
ARRAY(I,J + 1)
ARRAY(I,J + 2)
J
contains the. area of interface
(I,IARAY(I,J))
contains the E of interface
(I,IARAY(I,J))
contains the Q of interface
(I,IARAY(I,J))
is an integer variable which
is incremented by 3 from 1 to 16
For computational purposes, every interface must be accounted
for in terms of the section to which it applies. This would involve
duplicate input — for instance, AREA (I,J) is the same as AREA (J,I).
To simplify input, only one of the values must be input and subroutine
DATA (statements DATA 075 to 091) fills out the remaining positions in
ARRAY.
Subroutine REVI...
This subroutine is used to revise the parameters E and/or Q
without altering any other elements of the input. This can often
be an advantage when doing sensitivity analyses. The indicator ICON
is employed to call the subroutine and designates which parameter is
to be changed.
REVI also prints the revised values in a labelled tabular form.
C-5
-------
Subroutine FLOW
FLOW basically computes a balance of the flows into and out
of a section. If there is more or less than a zero balance around
the section, the resultant value is printed as an excess..flow.
The purpose of this exercise is to both keep track of
flows and to insure that no mistakes were made when inputting the flow
regime. This can be particularly useful when working in two or three
dimensions.
Subroutine SETUP
The primary purpose of this subroutine is to set up the system
matrix, AB, as was generally outlined on pages T-17 to T-18 . The tasks
performed by SETUP are:
1) calculation of the bulk dispersion coefficient, E1
According to the theory:
E'i i = Ei,j Ai,j
where
E-£sj = the dispersion coefficient across the interface
A,- a = the cross-sectional area of the interface
-1- > J
li,j = the average length of the adjoining
sections = (1-^ + lj)/2
The program calculates E' (EPRIM) as follows:
EPR1M(J)=ARRAY(I,JJ)*ARRAY(I,JJ+l)*417.166/(L(K,I)+L(I,K))
C-6
-------
where
ARRAY(I,JJ) = A. . (ft2)
i,3
ARRAY(I,JJ+1) = E. ,(mi2/day)
1 9 J
L(I,K) = 1±
L(K,I) = lk
417.1166 - conversion factor which is computed as follows:
MO MGD
^"gnl
417.1166.
2) A value of alpha is calculated to adjust the advective term
for sections of different length. This value is checked against
the positivity criterion (equation T-39) and if it is not met,
alpha is recalculated (as in equation T-44)
Alpha(J) = 1.0-(EPRIM(J)/(2*ABS(ARRAY(I,JJ+2))))
which corresponds to
= 1 - E'/2Q
3) The elements of the system matrix [AB] are calculated
4) The values for EPRIM and ALPHA are printed.
5) The section parameters: depth, volume and temperature, are
printed.
C-7
-------
Subroutine BCC
The primary purpose of this subroutine is to introduce boundary
conditions into the algorithm. BCC marks the point at which the type,
of constituent being modeled becomes relevant. At the present time the
program distinguishes between three types of constituents:
specific type
in HAR03 which is analogous to
Chlorides
BOD
D.O. deficit
any conservative substance
a single constituent reacting with
first order kinetics
the second constituent oE a coupled,
feed-forward system reacting with first
order kinetics.
BCC performs the following tasks:
1) Seads the number of boundary conditions
2) If there are none it writes the message
ZERO 30UwDARi' CONDITION
3) if there are boundary conditions it reads-and writes
them with the appropriate labels (CL, BOD or DEF)
4} it calculates the part of the forcing function due
to the boundary condition as described on page
for Q < 0
FRCBC(I)=8.34*BC(KLAST)*(EPRIM(J)-ALPHA(J)^ARSAY(I,JJ+2))
for Q > Q
FRCBC(I)=8.34*BC(KLAST)*(EPRIM(J)-(l.0-ALPHA(J)*ARRAY(I,JJ+2))
C-8
-------
where
FRCBC(I) = forcing function
BC = boundary concentration
EPRIM = E' (MGD)
ARRAY(I, JJ+2) = Q (MGD)
8.34 = conversion factor of
due to the boundary condition
(mg/1)
(MGD*mg/l) to #/day
Subroutine RATE
RATE finishes preparation for calculation by completing the
construction of the specific constituent system matrix [AC] and
the appropriate forcing function. It does this by:
1) Reading the reaction rates (@ T=20° C) which are
to be added to the diagonal of [AB] to transform
it to [AC]. In the BOD-DO deficit system these rates
are equal to
a) The removal rate, Kr or Kn, for CBOD and NBOD,
respectively
b) the reaeration rate, Ka, for deficit
c) zero for chlorides
At the same time it inputs the appropriate temperature correction
factors which are used to adjust the rates via the following general
formula
Kt = K20 *9 (T_20)
where some typical values are
rate 0
Kr 1.047
Ka 1.024
2) the rates are multiplied by the section volumes and added
to the diagonal of [AB] converting it to [AC]
3) for BOD or chlorides the forcing function is formed by
reading the waste loads (LOAD) and adding them to FRCBC
which was calculated in subroutine BCC
-------
4) for D.O. deficit the forcing function is completed by
first reading the deoxygenation rate (K^) the benthal
deinand (BD) and the net photosynthesis effect (PMR).
Appropriate temperature correction factors are read and
applied to Kd and BD. The loads are read and then the
forcing function is formed (as in equation)
XK1(I)=8.34*V0L(I)*(XK1(I)*KD(I)*FL-PMR(I)
+BD(I)*H(I))+Y(I)
where
XK1(I) on the left hand side of the equation is the
deficit forcing function
XK1(I) on the right hand side of the equation is the
BOD concentration in the body of water as
calculated during the BOD iteration.
FL is a ratio relating ultimate to 5-day BOD
Y is the same as FRCBC as calculated in
subroutine BCC
then the deficit waste load is added.
XK1 = XK1 + LOAD
Subroutine WRITE
This subroutine writes the basic matrices before and after
inversion if this is required.
Subroutine MAIN
This subroutine inverts the matrix IAC] and multiplies it by-
the forcing function.
C-10
-------
Subroutine STORE...
Subroutine STORE performs a number of miscellaneous tasks:
1) writes the final concentrations for each
constituent being modeled
2) calculates the saturation value of oxygen ^
CS(I) = (1.-.000009*CL(I))*(14.652-.41022«T(I)
+ .0079910 * T(I) ** 2.-.000077774 * T(I)**3.)
where T(I) = the temperature of a section
CL(I) = the chloride concentration of the section
3) performs what can be called miscellaneous bookkeeping
functions. For instance, it retains values of a
particular constituent which would be necessary when
calculating another constituent. (e.g. the BOD is
an input to the calculation of DO deficit and therefore
must be retained) . It therefore plays an important r.ole
in the coordination of a particular configuration of
reactants.
C-ll
-------
c* *******#**#**************** *********~****#****«*#*$$$*******~#~**$***************~****»~*~* HAR03001
£$ 4 $$4*$$ $$$$**«$ $«4$:£*
-------
END
SUBROUTINE OATA
r
SUBROUTINE OATA
INPUTS DATA WHICH
CONSTITUENT BEING
HYORODYNAMIC DATA
WOULD NOT VARY WITH THE PARTICULAR
MODELED, I.E. PRIMARILY PHYSICAL AND
C
'C
C
C
"C
C
REAL L A( 10 0 » 6)
DIMENSION AREA(6)»E(6),Q(6),TITLE(10)
COMMON /ALG A E/NX,N, INDIC, I COM , MMM , NMK.
COMMON/BASS/MX, I PRNT,PRCBC{100)
COMNCN/COHOE/LA,SCALE(4),ARRAY( 100, 18)I ARAY< 100,6),NLA(100,6)
COMMON/PIKE/JCON,H( 100 ) ,T(100) ,VOL(100),A8< 100»100)
31 FORMAT(8E10.0)
100 FORMAT! 10A4,313,3X.4F7.0)
101 FORMAT (///28X» 'ZERO SEGMENT NUMBER IN INTERFACE 13, '-', 13,
1' COMPUTATION DISCONTINUED')
FORMAT(6(14,
I4»?X,F7.C,2X))
HAR03054
DATA 000
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
SEPT
DATA
DATA
DATA
SEPT
DATA
DATA
SEPT
DATA
DATA
DATA
DATA
'///« INTERFACE. LENGTH
INTERFACE LENGTH
14)
2C0
201 FORMAT(6(I3,F10.0))
300 FORMAT( 'I
1 CHARAC. LENGTHS OF SEGMENTS (FT)
IFACE LENGTH INTERFACE LENGTH
2 E LENGTH INTERFACE LENGTH')
1111 FORMAT(3(3F6.0,13))
2000 FORMAT( ' 1V/40X, 10A4////10X,*** SYSTEM PARAMETERS**',//10X,
1'NUMBER OF SECTIONS = 9 »I 4,25X » ' IPRNT = 8 ,I 4/10X,*JCON = •
2002 FORMATI//10X,•~~SCALE FACTORS**'//10X,
1 9 SCALE( 1 ) = ',F10.3,4X, 'SCALE(2 ) = ',F10.3,4X,
2'SCAL E(3) = F10.3.4X, 'SCALE(4) = '»F10.3)
2200 FORMAT( ,2X, ' INTERFACE',5X,'AREA',7X,'E',8X,* Q ',2(6X,'I NTERFACESEPT
1 • ,5X,•AREA",9X,•E',7X,'Q' )/,4X,'ROW-SEG ' ,4X,• (FT**2)* , 2X, SEPT
SEPT
SEPT
SEPT
SEPT
SEPT
DATA
INTERDATA
INTERFACDATA
DATA
DATA
DATA
SEPT
DATA
DATA
DATA
2' (MI**2/D)•,2X,•(CFS) ', 6X,•ROW-SEG•,5X , •(FT**2)•,2X,
3'(M I**2/D ) *
4 2X,'(CFS)•
3000 FORMAT( '
1
* _ i
13, • )
,2X,•(CFS)',5X,'ROW-SEG',5X,'(FT**2) ',2X,' (MI**2/D) •
)
{*,I3,'-«,I3,')*,2X,F8.0,4X,F6.3,F8.1,2(5X,'(,,I3,
,2X,F8.0»4X,F6.3,F8.1) )
DATA
£
C
c
READ AND WRITE THE SYSTEM PARAMETERS AND SCALE FACTORS
SCALE
READ(MX,100)TITLE,N,IPRNT,JCON,
WRITE(NX,2000)TITLE,N,IPRNT,JCON
WRITE(NX,2002)SCALE
C
C
C
INITIALIZE MATRICES TO ZERO
~ DATA
~DATA
~ DATA
DATA
SEPT
SEPT
DATA
DATA
~ DATA
~ DATA
~ DATA
001
002
003
004
005
006
007
008
009
010
011
74
013
014
015
74
017
018
74
020
021
022
023
024
025
026
027
029
031
74
036
037
038
74
74
74
74
74
74
74
042
043
044
045
046
74
74
049
050
051
052
053
C-13
-------
DO 1001 1=1,N
DO 1021 J = 1,18
1021 ARRAY(I,J)=0.0
DO 1001 J = 1, 6
LA (IfJ)=0.0
IARAY( I , J )=0
1001
c**=»*
c
c
c
c
c****
READ AND WRITE INTERFACE PARAMETERS AREA,
THEM INTO ARRAY
E AND Q AND PLACE
72
74
71
1002
1088
1012
DO 1002 1=1,N
READtMX,1111) (AREA(J),E(J ) ,
DO 1002 J = 1» 6
IF{ AREAU))72,71,72
IF{IARAY(I,J))74,74,7 1
WRITE (NX,101)I,IARAYfI,J )
CALL EXIT
J J =(J-l)* 3+1
ARRAY
-------
1032 CONTINUE DATA 100
C********************************************************************** DATA 101
C *DATA 102
C READ AND WRITE INTERFACE PARAMETER- CHARACTERISTIC LENGTH *DATA 103
C *DAT A 104
£********************************************************************** OATA 105
WRITE(NX, 300) DATA 106
L INE=0 SEPT 74
DO 500 1=1,N DATA 107
READ(MX,201)(NLA(I * J),LA(I,J), J=l,6) SEPT 74
DO 77 J=l,6 SEPT 74
77 LA( I,J)=LA(I,J)*SCALE(4) SEPT 74
WRITE(NX,200)(I»NLA(I,J),LA(I,J),J=l,6) SEPT 74
L INE = L INE + 1 SEPT 74
IF(FL0AT(LINE/54)-FL0AT(LINE)/54«)500,1098,1098 SEPT 74
1098 WRITE(NX,300) SEPT 74
LI N E = 0 SEPT 74
500 CONTINUE DATA 118
C********************************************************************** DATA 119
C *DAT A 120
C READ SECTION PARAMETERS- DEPTH, VOLUME AND TEMPERATURE *DATA 121
C *DATA 122
C********************************************************************** DATA 123
READ(MX,31 ) (H< I ) , 1 = 1,N ) DATA 124
READ(MXf 31) (T(I),VOL(I),1 = 1fN) DATA 125
IF(T(2) ) 14,13»14 DATA 126
13 DO 14 I = 1»N DATA 127
T( I ) = T(1) DATA 128
14 CONTINUE DATA 129
C********************************************************************** DATA 130
C *DAT A 131
C APPLY CONVERSION FACTORS DATA 132
C *DAT A 133
q********************************************************************** DATA 134
DO 12 1=1,N DATA 135
DO 23 J=1,6 DATA 136
JJ=(J—1)*3+1 DATA 137
23 ARRAY(I,JJ+2)=ARRAY(I,JJ+2)*.6463 DATA 138
12 VOL(I )=VOL(I )*7.48 DATA 139
RETURN DATA 140
END DATA 142
SUBROUTINE REVI REVI 000
q********************************************************************** REVI 001
C********************************************************************** REV I 002
C REVI 003
C SUBROUTINE REVI REVI 004
C REVI 005
C THIS SUBROUTINE CAN BE USED TO CHANGE THE PARAMETERS E AND/OR REVI 006
C Q WITHOUT ALTERING THE REST OF THE INPUT REVI 007
C REVI 008
q********************************************************************** REVI 009
f.********************************************************************** REVI 010
REAL LA{100,6) SEPT 74
DIMENSION Q(6),E(6),REVIS(12) REVI Oil
C- 15
-------
COMMON/ALGAE/NX,N,INDIC,ICON,MMM,NMK REVI 012
COMMON/BASS/MX,IPRNT,FRCBC(ICO) SEPT 74
COMMON/COHOE/LA,SCALE(4),ARRAY( 100, 18') , I ARAY( 100,6) ,NLA (100 ,6 ) SEPT 74
2 FORMAT!1H1,40X•REVISED PARAMETER LIST'//) REVI 014
5 FORMATt12F5.0) REVI 015
10 F0RMAT(//50X,'NEW DISPERSIONS 1,/6X,5(•INTERFACE•,3X,•E•,7X),•INTERREVI 016
1FACE' , 3X, • E 1 / 15X, 5 ( MMREVI 017
7I**2/DAY)•,9X) ,MMI**2/DAY)•//) REVI 018
13 F0RMAT(3X,6(3X,I4,,-*,I3»F9«3)) REVI 019
14 FORMAT(//50X'NEW FLOWS'/7X»6(•INTERFACE',4X•Q',6X)/18X,5(•( CFS )* »1REVI 020"
65X),'(CFS)•//) REVI 021
WRITE (NX,2) REVI 022
GO TO(8 , 9 ) ,I CON REVI 023
C *REVI 025
C CHANGE E'S ~REVI 026
C *REVI 027
8 DO 6 1=1,N REVI 029
I F( I—(I/2)*2)21»21,20 REVI 030
20 READ(MX,5)REVIS REVI 031
IP0S=0 REVI 032
21 DO 6 J= 1,6 REVI 033
IP0S=IP0S+1 REVI 034
JJ=(J-l)*3+1 REVI 035
6 ARRAY(I,JJ+l)=REVIS(IP0S)*SCALE(2) REVI 036
WRITE (NX,10) REVI 037
DO 11 1=1, N REVI 038
DO 12 J = 1,6 REVI 039
JJ=(J-1)*3+1 REVI 040
12 E(J )=ARRAY(I,JJ + 1)/SCALE(2) REVI 041
11 WRITE (NX,13)(I, IARAY (I,J), E(J), J=I, 6) REVI 042
GO TO 99 REVI 043
c *REVI 045
C CHANGE Q'S *REVI 046
C *REVI 047
I 048
9 DO 4 1=1,N REVI 049
IF( I-( I/2)*2)23, 23,24 REVI 050
24 READ(MX,5 )REVIS REV I 051
IP0S=0 REVI 052
23 DO 4 J = 1,6 REVI 053
IP0S=IP0S+1 REVI 054
JJ=(J-l)*3+l REVI 055
4 ARRAY{ I ,JJ + 2 )=REVIS(IPOS ) *SCALE(3)*0.6463 REVI 056
WRITE (NX,14) REVI 057
DO 15 1=1, N REVI 058
DO 16 J = 1 ,6 REVI 059
JJ=(J-1)*3+1 REVI 060
16 Q(J ) = ARRAY(I,JJ + 2)/(SCALE(3)*0.6463) REVI 061
15 WRITE (NX,13) (I• IARAY( I ,J ) , Q(J), J=l, 6) REVI 062
99 RETURN REVI 063
END REVI 064
C-16
-------
SUBROUTINE FLOW
C***************************************************************#$#****
C
SUBROUTINE FLOW
THIS SUBROUTINE CALCULATES A FLOW BALANCE AROUND EACH SECTION
TO INSURE THAT THE FLOW REGIME HAS BEEN INPUT CORRECTLY
C
C
C
C
C
C'**********************************************************************
c**********************************************************************
REAL LA(ICO,6)
COMMON/ALGAE/NX,N,INDIC,ICON,MMM,NMK
COMMON/COHOE/LA,SCALE(4)* ARRAY(100*18)* IARAY(100,6) , NLA{100,6)
NMN = 0
214
SECT ION
HAS AN EXCESS FLOW
OF
F12.5
216
37
29
27
213
FORMAT(*
1' CPS')
302 FORMAT(• 1' )
DO 213 1=1,N
GQ=ARRAY( 1,3 ) +ARRAY(I» 6 ) + ARR AYI 1,9)+ARRAY(I,12)+ARRAY(1,15)
1+ARRAY(1,18)
QQ=GG/•6463
IF!GC-,001)216,37,37
IF(QQ+.001)37,37,213
IF(NMN)29,29,27
NMN= 1
WRITE(NX,302 )
WRITE(NX,214 ) I,GQ
CONTINUE
RETURN
END
SUBROUTINE SETUP
C**********************************************************************
c**********************************************************************
c
C SUBROUTINE SETUP
C
C THE PRIMARY TASK OF SETUP IS TO CONSTRUCT THE MATRIX, AB
C
C**********************************************************************
^~~^~^~~~^~~~~~~~~X'****************************************************
REAL LA(100,6)
COMMON/ALGAE/NX,N,INDIC,ICON,MMM,NMK
COMMON/COHOE/LA,SCALE(4) » ARRAY ( 100 » 18) .* I ARAY ( 100,6) , NLA ( 100,6)
COMMON/PIKE/JCON,H(100),T(100),VOL(100),AB(100,100)
COMMON/BREAM/EPRIM(100),ALPHA(100)
FORMAT(6(13, '—•,I3,F6.0,F6.3,2X) )
FORMATC °1•,49X,"VALUES OF ALPHA AND
1PHA •)t/t6(8X,'(MGD)«,8X),//)
42 FORMAT(4X»I7,3F10«2)
23
24
EPRI.M® , / / / , 6 ( ' INTRFC EPR I M
43 FORMATl'IS3X,"SECTION
115X , ' ( C ) U0**6GAL)
WRITE (NX,24)
DO 28 1=1,N
TEMPERATURE
(FT)•//)
VOLUME
DEPTH'/
FLOW 000
FLOW 001
FLOW 002
FLOW 003
FLOW 004
FLOW 005
~FLOW 006
~FLOW 007
~FLOW 008
FLOW 009
FLOW 010
SEPT 74
FLOW 012
SEPT 74
FLOW 014
FLOW 015
FLOW 016
FLOW 017
FLOW 018
FLOW 019
FLOW 020
FLOW 021
FLOW 022
FLOW 023
FLOW 024
FLOW 025
FLOW 026
FLOW 027
FLOW 028
FLOW 029
FLOW 030
SETUPOOO
SETUP001
SETUP002
SETUP003
SETUP004
SETUP005
S6TUP006
SETUP007
SETUP008
SETUP009
SEPT 74
SETUP011
SEPT 74
SETUPO13
SETUP014
SEPT 74
ALSETUPO16
SETUP017
SETUPO18
SETUPO19
SETUP020
SETUP021
SETUP022
C-17
-------
DO 99 J= I, N SETUP023
99 A3 ( 11J)=0 « 0 SETUP024
S ETUP025
C ~SETUP026
C CALCULATE EPRIM AND ALPHA FOR EACH INTERFACE *SETUP027
C *S ETUP028
SETUP029
DO 20 J = 1,6 SETUP030
EPRIM(J)=0.0 SETUP031
ALPHA(J ) =0.0 SETUP032
JJ=(J-1)*3+1 SETUP033
K=IARAY(I,J) SETUP034
I F(K)20,20,10 SETUP035
C SEPT 74
C GENERATE (I,K) LENGTH ELEMENT SEPT 74
C SEPT 74
C ****£« *»ilr******t78,84,78 SEPT 74
84 XLEN2=LA(K,KK) SEPT 74
GO TO 83 SEPT 74
78 CONTINUE SEPT 74
83 EPRIM 22,20,22 SETUP056
C- 18
-------
2 2 AB( I , K ) = (1 .-ALPHA(J))*ARRAY(I,JJ + 2)-EPRIM(J)
20 CONTINUE
C**********************************************************************
c
C WRITE VALUES OF EPR IM AND ALPHA
C
£*$$$44* 4 4*4* * 44* 4 44* 44 ****** * ****** ***** **** 4*44444*** *44$** *444*4
WRITE(NX,23)
II, IARAY(1,1 ),EPRIM(I)»ALPHA(1 ) » I,IARAY(I,2),EPR IM(2)*ALPHA(2),
21,IARAY(1,3),EPRIM(3),ALPHA < 3),I,IARAY(I,4),EPR IM(4),ALPHA(4),
21, IARAY( 1,5),EPRIM(5)» ALPHA(5 ) , I,IARAY (I,6),EPRIM(6),AL PHA(6)
L INE = LINE + 1
IF(FLOAT(LINE/54)-FLOAT(LINE)/54.)28, 1028, 1028
1028 WRITE(NX,23)
L INE = 0
28 CONTINUE
C4*444**4**************************************************************
c
C WRITE SECTION PARAMETERS
C
Q444*444*4*****4*44******4***********4*******4*************************
WRITE(NX,43)
L IN E = 0
DO 399 1=1,N
WRIT E(NX,42) I *T(I)«VOL( I ) *H(I) '
LINE=LINE+1
I F(FLOAT(LINE/54)-FLOAT(LINE)/54.)399,1029,1029
1029 WRITE(NX,43)
L IN E = 0
399 CONTINUE
DO 32 1=1,N
32 H(I)=3«281/H(I)
RETURN
END
SUBROUTINE BCC
£44444444444444444444444*4***4*****************************************
£444444444**44*44*44*$444444*444444444444444444444:4444**********4*44**4
C
C SUBROUTINE BCC
C
C THIS SUBROUTINE INCORPORATES BOUNDARY CONDITIONS INTO THE
C SYSTEM
C
£4$444444444**44*4**$4444444444444**44444*44*4444*44*4*****************
C44444*444*4****4******************************************************
REAL LA(100,6)
DIMENSION ICOL(IOO),BC(100)
COMMON/ALGAE/NX,N,INDIC , I CON,MMM,NMK
C0MM0N/BASS/MX9IPRNT,FRCBC(100)
COMMON/BREAM/EPRIMl100),ALPHA(100)
C0MM0N/C0H0E/LA,SCALE(4),ARRAY(100,18),IARAY(100,6),NLA(100,6)
7 FORMAT(12)
44 FORMAT( 81® ,4CX» 'ZERO BOUNDARY CONCENTRATIONS' )
209 FORMAT( ' 1SEGPENT®,5X, 'CL BOUNDARY CONDITI ON(MG/L)
SETUP057
SETUP058
~SETUP059
*SETUP060
*SETUP061
*SETUP062
*S ETUP063
S ETUP064
SETUP065
SETUP066
SETUP067
SEPT 74
SEPT 74
SEPT 74
SEPT 74
S ETUP068
SETUP069
*SETUP070
*SETUP071
~SETUP072
SETUP073
SETUP074
SEPT
74
SEPT
74
SEPT
74
SEPT
74
SEPT
74
SEPT
74
SEPT
74
SEPT
74
SETUP076
SETUP077
SETUP078
SETUP079
BCC
000
BCC
001
BCC
002
BCC
003
BCC
004
BCC
005
BCC
006
BCC
007
BCC
008
BCC
009
BCC
010
SEPT
74
BCC
012
BCC
013
BCC
014
BCC.
016
SEPT
74
BCC
017
BCC
018
BCC
019
C-19
-------
1 /<14, 19X,F10.2)) BCC 020
210 FORMATE * lSEGf*ENT' ,5X, 'BOO BOUNDARY CCND I T ION ( MG/L 3 • BCC 021
1 /( 14, 19X,F10.2M BCC 022
211 FORMATf'1SEGMENT",5X,*DfcF BOUNDARY CONDITION{MG/L) » BCC 023
1 /(14,19X.F10.2)) BCC 024
1150 FORMAT(6( F8,0,I4]) BCC 025
£**4#*4*4444444444444444 + 444 44444 4444444*44444*44444444«4***4*4«44*4*#4 BCC 026
C *BCC 027
C READ AND WRITE BOUNDARY CONDITIONS *BCC 028
C *BCC 029
4*4**4**4*444 44*44444444 4 44444 444444 44*4444444**4*4*:fc* BCC 030
8 READtMX,7JNUMBC BCC 031
1F{NUMBC)31,31,36 BCC 032
31 WRITEINX,44) BCC 033
DO 103 1=1,50 BCC 034
1C0H 1)^0 BCC 035
103 BC(I)=0• BCC 036
GO TO 213 BCC 037
36 READtMX,1150)(BCfJ) , I COL (J ) , J = 1. NUMBC) BCC 038
IF( INDICJ206,207,208 BCC 039
206 WfUTE{NXt209) ( ICOU J J ,BCU) , J = i,NUMBC) BCC 040
GO TO 213 BCC 041
207 WRITE(MX,210M I COL ( J ) ,BC( J), J = l, NUMBC) BCC 042
GO TO 213 BCC 043
208 WRITE(NX,21l)(ICOL{J)»BC(J)« J=1»NUM8C) BCC 044
q*****4444****44*44444*4*4*444*444444*4444444*44***4444444*44*4*44***#* BCC 045
C *BCC 046
C CALCULATE THE PART OF THE FORCING FUNCTION CUE TO BOUNDARY COND. *SCC 047
C *BCC 048
£*£***4*4*4:44:* 4*4*4 44 *444 4# 44* 4***4 44 £*$*4 ***4 444 4* *4 44****4*********** BCC 049
213 KLAST = 1 BCC 050
DO 60 I •= I» N BCC 051
FRC 8C(I)=0.0 BCC 052
IF(ICOL(KLASTI-1)60,25,60 BCC 053
25 DO 30 J = 1,6 BCC 054
JJ = U-l)*3+l BCC 055
K=IARAYI I,J) SEPT 74
IF(K-I) 30,35,30 SEPT 74
30 CONTINUE BCC 058
£$$$4 $*$$4 $4 4*444*4************************************5EPT 74
c SEPT 74
C GENERATE |[,K! LENGTH ELEMENT SEPT 74
C SEPT 74
Q******$**4**44*4**4*******************************«*44*****************5EPT 74
35 00 77 KK=1,6 SEPT 74
IF(NLA« I ,KK )-K>77,81,77 SEPT 74
61 XL EN L = L A( I, KK ) SEPT 74
GO TO 82 SEPT 74
77 CONTINUE SEPT 74
£************** **** ** **4* 4 4*4 4444444 ************************************SEPT 74
C SEPT 74
C GENERATE tK, I ) LENGTH ELEMENT SEPT 74
C SEPT 74
£*********4**4**4*4************4*4#*4*********+*************************SEPT 74
C- 2 0
-------
8 2 DO 78 KK=1»6
I F(MLA{K t KK J-I)78»04,73
84 XL EN2=LA[K»KK)
GG TO 83
78 CONTINUE
83 A1_PHA< J 1=0.5
EPRIMJ )=ARRAYC I ,JJ}*ARRAY(I * J J+l)*417.1166/(KLEN2+XLEN1)
IF!-0.5+EPR!W* J )/ABS i ARRAYi + MO, 45,45
40 ALPHA? J ) = l.0-6PR.IM { JJ/( 2, 0*A8S{ ARRAY (I, JJ + 2 ) ) >
45 IF(ARRAY{ I ?JO+21)50* 55s 55
¦>0 FRCBCf I )=a.34*8CtKLAS~ >*SEPRIM{ J J-ALPHAl J ) * ARRAY I I»JJ*2))
K.LAST=Kt_AST+l
GC TO f>C
5 5 FRCBCI I J = 0«34*BClKLAST!*{EPRlfKJl — {l»0-ALPflA(J> >* ARitAYI I, JJ+2 3 3
KL 9 ST=S\L A ST-9-1
60 COkTlbSJE
10C0 RETURN
END
SUBROUTINE RATE
c
c subroutine rate
c
C THE PRIMARY PURPOSES CF THIS SUBROUTINE ARE TO
C SPECIFIC CONSTITUENT SYSTEM MATRfXt AC AND TG
C APPROPRIATE FORCING FUNCTIONS
C
REAL KA { LOOUKO< 100t»KSTORI 100J »LOAD
0I MEMS I ON YaOOJ,QD!iOO) ,PM«f 100)
COMMON/ALCAE/NX »N»IND1C»ICOMs
CQMMON/QASS/WX, IPRMT,FRCBCUQO}
CQHM0N/PlK.E/JCaM,M100niaCC! ,VQH LG01 *A61100» 8.00)
COMMON/TROUT 100 3 t X!<1»1001 *6C100* LI,4C(100,100i
FORHAH 16F5.Q J
CONSTRUCT THE
CGMPLET£ THE
1
12
16
17
13
34
35
36
FORMAT(FLO.
FORMAT(*
Do 13 3
LOAD =s
LOAD =f
t FlOr
sf 10.
PGUMD5/DAY
POUNDS/OAY
FOR
FOR
SECTION
SECTION
13!
131
BOUNDARY LOAD ("IGD*^&/L > «///)
BOD
FQRCAT t * CEF
FORMATS"Is I
FORMAT[* 1 SECTION CHL OK IDE
FCRMATi2Xpi5i12XVF14*2)
FORMAT{ ¦I THE TEMPERATURE CCRRECTICM FACTOR FCR THE BOD REMOVAL
IRATE = »,F7.3,///,° SECTION BOD REM RATE 800 WASTE LOAD BOD
LBOUNDARY LGAC/12X* <*T£MP C0RR*'t5Xf ?(tfGD*MG/L) !?7X,e (MGD*MG7l) • )
37 FORMAT!2X,15,F12„3,7X»F12.3,5X,F12.3)
38 FORWATCl THE TEMPERATURE CORRECTION FACTOR FOR THE REAERATION RARATE
ITE = *» F7« 3»/' THE 600 CORRECTION FACTOR (FL! = rsF7.3,/« THE RATE
1TEMPERATURE COSRECTIOM FACTOR FOR THE BENTHfl L DEMAND = *,F703,/' RATE
1 THE TEMP CORR FACTOR FOR KD = %F7„3,///• SECTION DEOX RATE RATE
1 REAEft RATE BOUND COHD LDAD TOTAL LOAD PMR BO */»12RATE
I W., 2 ( C*TEMP CC1RR* lMGQ*MG/t i %430 f 'MG/L/O *TEHP COR«=M J RATE
39 FORMATS 2X, 16i5K¥F6» 3 » ?X,F6»3 ,6X,FL2.3e3X,F 12* 3,2K9 F9«.3, 3X» F9»3} «ftTE
SEPT
SEPT
SEPT
SEPT
SEPT
SEPT
SEPT
BCC
BCC
BCC
BCC
BCC
acc
&cc
BCC
BCC
BCC
8CC
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
74
74
74
74
74
74
74
061
062
06 3
064
06 5
os.b
067
066
069
070
071
COO
001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
0La
019
02C
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
C~21
-------
C ~RATE 037
C READ IN DECAY RATE FOR BOD SYSTEM* REAERATION RATE FOR DO DEFIC IT~RATE 038
C SYSTEM RATE 039
C IF THE K•S ARE THE SAME FOR EVERY SECTION THEN JUST INPUT K(1) ~RATE 040
C READ AND APPLY TEMPERATURE CORRECTION FACTORS ~RATE 041
C ~RATE 042
£*$*******«****$**«*******$**#*******#*******************$**************£aje 043
I F(MMM-5)8,386,386 RATE 044
8 READ(MX,1)FAC1,(KA< I ) , 1 = 1,15) RATE 045
I F(K A(2 ) !202,3»202 RATE 046
3 DO 64 I = 2 , N RATE 047
64 KA(I ) =KA{I) RATE 048
GO TO 334 RATE 049
202 LMM=N-15 RATE 050
IF(LMM)334,334,203 RATE 051
203 READ(MX»1) (KA( I ), 1=16,NJ RATE 052
334 I F(MMM-3 ) 2,388, 2 RATE 053
388 DO 389 1=1,N RATE 054
389 KSTOR(I ) =KA( I) RATE 055
F AC4-F AC 1 RATE 056
GO TO 2 RATE 057
386 DO 387 1=1,N RATE 058
387 KA( I)=KSTOR(I ) RATE 059
FAC1=FAC4 RATE 060
2 DO 6 1=1,N RATE 061
K A( I )=KA( I ) ~FAC 1 ~~(T(I)-20.) RATE 062
RATE 063
c ~RATE 064
C CALCULATE THE SPECIFIC CONSTITUENT SYSTEM MATRIX, AC ~RATE 065
C ~RATE 066
£$$$$ 4 $4 4 RATE 067
DO 801 J=1,N RATE 068
801 AC( I,J )=A8( I,J) RATE 069
6 AC(1,1 ) = A8(I,I)+KA( I)~VOL(I) RATE 070
IF(INDIC)10,10,54 RATE 071
£*$ 4^ 4t 4t4r4c« *4c4t4^44c#*4^*4'4'$ ~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~R ATE 072
c ~RATE 073
C CALCULATE THE BOD OR CHLORIDE FORCING FUNCTION ~RATE 074
C ~RATE 075
Q4c$4^4^4^4t*4^**4^*4*:<:**4'«**4c4c4^**4'*4:**4'**********4'*4i***4:4c4<*4:*****4:*4'**4c*RATE 076
10 DO 9 1=1,N RATE 077
W(=*********~<=*~ *~~~~~~~~*~*~*~~~~~** ~~~~~~~RATE 089
C-22
-------
READ IN DEOXYGENATION RATE
IF THE KD * S ARE THE SAME FOR
READ IN PHOTOSYNTHESIS MINUS RESPI RATI ON(PMR),
READ IN BOTTOM DEMAND( BD) , GM./SQ. M./DAY
VALUES MUST BE READ IN FOR EACH SECTION
READ AND APPLY TEMPERATURE CORRECTION FACTORS
54 READ(MX,1)FAC2,FL,(KD
-------
GO TO 40
RATE
144
30
DO 304 1=1,N
RATE
145
304
W( I ) = W(I)/8.34
RATE
146
WRITE(NX,36)FAC1
RATE
147
DO 853 1=1,N
RATE
148
853
W( I ) = W( I)-FRCBC(I)
RATE
149
WRITE(NX,37)(I,KA(I)9W(I),FRCBC( I ) , 1 = 1»N)
RATE
150
DO 854 1=1,N
RATE
151
854
W( I ) = FRCBC( I )+W( I)
RATE
152
GO TO 40
RATE
153
31
DO 305 1=1,N
RATE
154
305
XK1(I) = XK1( I )/8.34
RATE
155
WRITE(NX,38)FAC1?FL»FAC3»FAC2
RATE
156
WRITE(NX,39)(I,KD!I ) ,KA(I ) ,FRCBC( I ) ,XK1(I ) ,PMR{ I),BD(I),I = 1,N)
RATE
157
40
IF(INDIC)843,84 3,288
RATE
158
288
DO 842 1=1,N
RATE
159
842
W ( I ) =XK 1 ( I)
RATE
160
843
DO 841 1=1,N
RATE
161
841
B ( I , 1 ) = W ( I )
RATE
162
RETURN
RATE
163
END
RATE
164
SUBROUTINE WRITE
WRITEOOO
c ***$$
c
c
c
c
c
800
803
1000
1005
1007
1012
1015
2000
1003
1018
1004
1006
1008
1001
SUBROUTINE WRITE
THIS SUBROUTINE OUTPUTS MATRIX AC AND MATRIX AC**-1
COMMON/ALGAE/NX,N,INDIC»I CON,MMM,NMK
COMMON/T ROUT/ W(100).XK1{100)•B<100.1),AC(100»100)
SYSTEM MATRIX')
SYSTEM MATRIX')
FORMAT('l',50X,'INITIAL
800 MATRIX9)
DO DEFICIT MATRIX')
BOD MATRIX*)
DO DEFICIT MATRIX*)
FORMAT(///48X,•INVERTED
FORMAT ( MX,9113/)
FORMAT(*1*,51X,'INITIAL
FORMAT ( ' 1 %48X, • INITIAL
FORMAT(///50X,«INVERTED
FORMAT(///47X,»INVERTED
FORMAT(I3,lXo9E13o5/(4X»9E13.5)5
JMAX=N
MM= 1
IF(NMK)1003,1003,1010
I F( INDIC)1018,ICOA,1006
WR I TE(NX » 800)
GO TO 1008
WRITE(NX,1005)
GO TO 1008
WRITEtNX,1007)
WRITE(NX, 1000) (J » J = Is JMAX)
DO 1001 1 = 1,JMA X
WRITE(NX,2000)I,(AC(I« J),J=1 , JMAX)
NMK=NMK+1
WRITEOOl
WRITE002
WRITE003
WRITE004
WRITE005
WRITE006
WRITE007
WRITE008
WRITE009
WRITE010
WRITE011
WRITE012
WRITE013
WRITE014
WRITE015
WRITE016
WRITE017
WRITE018
WRITE019
WRITE020
WRITE021
WRITE022
WRITE023
WRITE024
WRITE025
WRITE026
WRITE027
WRITE028
WRITE029
WRITE030
WRITE031
WRITE032
C-24
-------
RETURN
1010 IF< INDIO1019,1011,1014
1019 WR I T E(NX , 803 )
GO TO 1016
1011 WRITE(NX,1012)
GO TO 1016
1014 WRI T E (NX»1015)
1016 WRITE(NX,1000) (J,J=1,JMAX)
DO 1002 1 = 1, JMAX
1002 WRITE(NX,2000)I,(AC(I,J),J=l,JMAX)
NMK = 0
RETURN
END
SUBROUTINE MATN
SUBROUTINE MATN
THIS SUBROUTINE INVERTS THE MATRIX AC AND MULTIPLIES
IT BY THE FORCING FUNCTION
DIMENSION IPIVO(100),INDEX(100,2),PIVOT(100)
COMMON/ALGAE/NX,N,INDIC,ICON,MMM,NMK
COMMON/TROUT/W{100)fXKI(100),B(100,1),AC{100,100)
DO 900 1=1,N
B( I, 1 ) = B 4 1,1)/lOCCCOOo
DO 900 J=1» N
900 AC( I, J )=AC1I»J)/10CC000.
MM=1
DO 20 J = 1» N
20 IPIVO (J)=0
DO 550 1=1,N
C
C
C
SEARCH FOR PIVOT ELEMENT
60
80
85
100
105
IF
DO
IF
IF
60, 105, 60
100, 100
C
C
c
AMAX=0.0
DO 105 J=l,N
(IPIVO (J )-1S
100 M=1,N
(IPIVO (M ) — 1J 80, 100, 740
(ABS(AMAX)-ABS(AC(J,Ml))85,
IROW=J
ICOLU =M
AMAX=AC(J,M)
CONTINUE
CONTINUE
IPIVO (ICOLU )=IPIVO (ICOLU >~!
INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
IF (I ROW-ICOLU ) 150, 260, 150
WRITE033
WRITE034
WRITE035
WRITE036
WRITE037
WRITE038
WRITE039
WRITE040
WRITE041
WRITE042
WRITE043
WRITE044
WRITE045
MATN 000
MATN 001
MATN 002
MATN 003
MATN 004
MATN 005
MATN 006
MATN 007
MATN 008
MATN 009
MATN 010
MATN Oil
MATN 012
MATN 013
MATN 014
MATN 015
MATN 016
MATN 017
MATN 018
SEPT 74
MATN 021
SEPT 74
MATN 023
MATN 024
MAIN 025
SEPT 74
SEPT 74
SEPT 74
MATN 029
SEPT 74
MATN 031
MATN 032
SEPT 74
SEPT 74
MATN 035
MATN 036
SEPT 74
MATN 038
MATN 039
MATN 040
SEPT 74
C-25
-------
150
200
210
250
260
C
c
c
c
c
c
DO 200 L=*1>N
SWAP=AC(IROWeL)
AC«IROMrL)sAClICOLU»L3
AC!ICOLU»L 1 =SWAP
IF(MM I 2 601 2 60 s 210
DO 250 L=1s MM
SWAP=8(iROW,Li
BIIROVIrL)seUCOLU »U
BdCOLU jD^SWAP
INDEX 1111) = IROW
INDEXt1»2)=IC0LU
PIVOT(I)=ACIICOLU,ICOLU)
DIVIDE PIVOT ROW BY PIVOT ELEMENT
AC IICOLU,IC01U)*1.0
DO 350 L=1*N
350 AC?ICOLli.L >=AC(ICOLUsL)/PI VOTE I I
T F t )380f 3 80p 360
360 00 370 L=I» MM
370 BUCQLU ,L>=B(ICOLU »L)/PIVOTm
REDUCE NON-PIVGT ROWS
00 550 L1= i 9 N
1 F) L 1 —ICOLU ) 40Df 550, 400
S=AC(LlsICOLU)
AC(L1»ICDLU>=0.0
DO 450 L=19N
ACILUL )«AC< LUM-ACC ICOLU*L}*S
1 F ( MM) 550 j 550, <*60
DQ 500
saitl } = BILi»L>-BS 1CDJ.U ,L>*S
CONTINUE
380
400
450
460
500
550
C
C
C
630
705
710
740
901
INTERCHANGE COLUMNS
00 T10 1*1,N
L=N+1-I
IF ( INDEX(L, U-lNCEXiL* 25 >
JRQW-INDEXtL. 1)
JCOLU - INDEX(L t 2)
00 705 M= 1, N
SWAP=ACIM,JROW)
AC
-------
SUBROUTINE STORE
£4444444444*444** 44444444444444444444***4**4*******4**4*44*4*4*********
£*4444 444 44**4**4*4**4*****44 44*44*4**4************4**$****************
c
C SUBROUTINE STORE
C
C****44444444444444444444444444444444*44444444444444444444*444444444444
04444444444444444444444444*444444444444444444444444444444444444444444**
DIMENSION CL t100)9CSI100),R i100)
COMMON/ALGAE/NX»N» INDIC s, ICON, MMM,NMK
COMMON/PIKE/JCON,H(100),T{100),VOL(100>,AB(100t100)
COMWON/TROUT/WI 100)«XK1(100),B(100,1),AC!100,100)
FORMAT(aI SECTION DEFICIT(MG/L)1
FORMAT(®1 SECTION BOD IMG/L)8)
FORMAT( '1 SECTION CHLORIOES I MG/L)
FORMAT(3X,I4»3X,F14.3)
FORMAT(10X,13,5(4X,F9« 2))
FORMAT ( 11 SECTION CHLORIDES
ISATURAT ION DO',//)
817 FORMAT( ' I ® , # RESULTS
42
44
45
48
807
308
)
)
TEMPERATURE
FOR COMBINED
STOREOOO
ST0RE001
S T0RE002
ST0RE003
STORE004
STORE005
ST0RE006
ST0RE007
ST0RE008
ST0RE009
STOREOIO
ST0RE011
STOREQ12
ST0RE013
STOREO14
ST0RE015
STOREO16
ST0RE017
STOREO18
(CARBONAC EOUSTORE019
DEFICIT
TEMPERATURE
IC)
DEFICIT
(MG/L)
4S AND NITROGENOUS ) RUN'//)
828 FORMATt' SECTION CHLORIDES
1SATURATI ON DO 6,/,19X,'(MG/L)
1 (MG/L) IMG/DM
DO 805 1=1,N
805 XK1{I) = B( 1,1)
C ********************************** 4s ******* 4 4 * ^ ****** 4 4 4 4* 4 4 4 4 4 4 * 4 4 4 4 4 4
C
C WRITES FINAL CONCENTRATIONS
C
£444*4**444**********4**44*********4**44444***4**4****4*444*444*44*4***
11 IF(INDIC)850,840e41
850 WRITE(NX » 45 )
GO TO 643
840 WRITE!NX» 44 )
GO TO 643
41 WRITE(NX,42)
643 WRITE I NX,48)(I,XK1< I),I=1,N)
9999 GO T0(1800,1801,1802,1801,1803),MHM
1800 GO T0U199, 1199,813, 813),JCON
C* **************** **************:<<************ **************** **********
C
C CALCULATES THE SATURATION VALUE OF DO
C
c**********************************************************************
813 DO 815 1=1,N
CL( I ) = XK1(1)
815 CStI)=(l.-»000009*CL(I))*(14.652-,
IT II ) *»2.000077774*T(I)**3.)
MMM=MMM+1
RETURN
1199 DO 1139 I=1,N
CL(I)=0.
1139 CS(I ) = (l.-.000009*CL(I))*(14.6 52-,
>41022*T(I)+„0079910*
-41022*T ( I ) + «.0079910*
STOREO 20
ST0REO21
SEPT 74
SEPT 74
ST0RE023
S T0RE024
ST0RE025
ST0RE026
ST0RE027
ST0RE028
ST0RE029
ST0RE030
ST0REO31
ST0RE032
ST0RE033
ST0RE034
STORE035
ST0RE036
ST0RE037
ST0RE038
S T0RE039
ST0RE040
ST0RE041
ST0RE042
ST0RE043
ST0RE044
ST0RE045
ST0RE046
ST0RE047
ST0RE048
ST0RE049
ST0RE050
STORE051
ST0RE052
C-2 7
-------
IT tI )**2.-.000077774*T(I )*«3.)
ST0RE053
MMM=MMM*2
ST0RE054
RETURN
ST0RE055
1801
MMM=WMM+1
ST0RE056
RETURN
ST0RE057
1802
WRITE(NX.808)
ST0RE058
DO 806 1=1,N
STORE059
W(I) = CS(I )-XKl(I )
ST0RE060
WR I TE ( NX» 807 ) I ,CLU ) ,TC I )
eXKU I ) ,CS( I ) , Wl I )
ST0RE061
806
CONTINUE
ST0RE062
GO T0(93,1804,93,1804),JCON
ST0RE063
93
RETURN
ST0RE064
1804
DO 829 1=1,N
ST0RE065
829
R(I ) =XK1{ I)
ST0RE066
MMM=V. MM + 1
ST0RE067
RETURN
S T0RE068
1803
WRITE(NX»817)
ST0RE069
WRITE(NX,828)
S T0RE070
DO 881 1=1,N
ST0RE071
XK1(I ) = XK1(I ) + R (I)
ST0RE072
W
-------
USER'S MANUAL
U—1
-------
The manner in which the data deck for HAR03 is structured
is shown in figure U-I. As can be seen, the cards can be divided into
two basic groups:
1) General System Cards (GS cards) which are the same for
any constituent or substance being modieled and need only
be input once. They primarily contain the physical and
hydrodynamic parameters of the system.
2} Specific Constituent Cards (SC cards) which as the name
implies, are specific to the particular constituent being
modeled. These include boundary concentrations and
biological or chemical parameters such as rates and loads.
After reading the single set of GS cards, a variable number of
sets- of SC cards are read depending on the particular combination of
reactions of the constituents being modeled. Two indicators are used
to designate which combination is to be modeled:
1) INDIC: As can be seen from figure U-l, this indicator
is read in prior to each set of SC cards. It can
presently take on three values, -1, 0 and 1 to
designate chlorides, BOD and DO deficit respectively
2) JCOM: This indicator is read on the first GS card:
the Title Card. It designates a particular
combination or configuration of reactions and
is utilized by subroutine STORE to direct the
flow of the program.
Since HAR03 has been designed with the BOD-DO deficit system in
mind, it is capable of handling reaction configurations which are to
sortie extent peculiar to that system. However, it should be noted that
there are many substances which are analagous to BOD and dissolved oxygen
and can therefore be modeled using HAR03. With this in mind, the following
are the particular combinations of reactants which can be presently
modeled using the program:
1) Conservative substances: (Chlorides; total dissolved
solids; some pesticides and herbicides; etc) any sub-
stance which does not decay or can be approximated by
zero decay.
2) Single reactive substance: (BOD; coliform bacteria,
ammonia; etc) any substance which decays according
to first order kinetics
3) Coupled reactive substances: (BOD-DO deficit;
polyphosphate - orthophosphate; ammonia-nitrate). A
feed forward system of two substances reacting with
first order kinetics.
U-2
-------
4) Additive coupled substances: (.CBOD, NBOD-DO deficit)
this is peculiar to the. BOD-DO deficit system where
the BOD can be broken down into carbonaceous (CBOD)
and nitrogenous (NOD) components. Both components
deplete dissolved oxygen but at different rates so
ea<:h is computed separately and their deficits added
to determine the total deficit.
5) Sstuarine coupled reactive substances (Chlorides,
BOD-deficit) RAR03 has the additional option for
those estuaries where chlorides have been modeled.
A chloride run can procede the BOD-deficit so that the
resulting chloride concentrations can be used to
determine the saturation value of dissolved oxygen.
6) Estuarine additive coupled system (Chlorides, NBOD,
CBOD-deficit). The same as 5 but for the additive
system.
Figure U-2 summarizes the various deck setups for each of the
above combinations with the appropriate values of 1NDIC and JCON.
A description of input and output is then follox^ed by an example
problem including sample input and output to further clarify the procedure.
Finally, some restrictions on the use of the program are
-------
Figure U-l; Structure of data deck for HAR03
U-^
-------
Combinations of. reactants
ChlorIdes
(ConservatIve
subs tances)
BOD
(S i ngl e
reac ti ve
subs tance)
BOO
DEFICIT
(coupled
reac tl ve
subs tances)
CBOD
NBOD
DEFICIT
(addi tlve
coupled
subs tances)
ChlorIdes
BOD
DEFICIT
(es tuarIne
coupl ed)
ChlorIdes
CBOD
NBOD
DEFICIT
(es tuarIne
add I tl ve)
J CON
CARD SETUP
JCL
GS CARDS
3Z
^INDIC-IEXIT
CARD
(I NO IC=-l>
SC CARDS
(ChlorIdes)
/INDIC
-IEXIT
CARD
(IEXI1>1)
/
/*
JCL
JCL
GS CARDS
GS CARDS
INDIC-I EX I 7
CARD
( INDIC = 0 )
INDIC-IEXIT
CARD
- (INDIC = 0 )
/
/
SC CARDS
SC CARD f
(BOD)
(BOD)
flNDIC-lEXIT
( INDIC-IEXIT
CARD
CARD
(1ND1C=1)
<1ND1C = 1)
¦ir
/
SC CARDS
SC CARD 5
(deficit)
(def1c11)
INDIC-IEXIT
CARD
( I EX I T= 1)
INDIC-IEXIT
CARD
(IEXIT=1)
/*
JCL
JCL
JCL
GS CARDS
GS CARDS
GS CARDS
INDIC-IEXIT
CARD
(INDIC=-lJ
/ INDIC-IEXIT
CARD
( INDIC —1)
SC CARD S
(ChlorIdes)
SC CARDS
(Chlor I des)
INDIC-I EX I 1
CARD
(INDIC=0)
INDIC-IEXIT
CARD
(INDIC = 0 3
INDIC-IEXIT
CARD
(INDIC=0)
4r
/
/
/
SC CARDS
SC
CARD S
SC CARD S
(CBOD)
(BOD)
(CB.OD)
^INDIC-IEXIT
INDIC
-1 EXIT
( INDIC-IEXIT
CARD
CARD
CARD
( 1ND1C = 1)
(1ND1C=1)
(1ND1C=1)
/
/
SC CARD 5
SC
CARDS
SC CARDS
(def i cIt)
(deflei t)
(def1c i t)
flNBIC-lEXIT
( INDIC-IEXIT
CARD
CARD
(1ND1C=0)
(1ND1C=0)
/
/
SC CARD S
SC
CARDS
(NBOD)
(NBOD)
(INDIC-
IEXIT
( INDIC-IEXIT
CARD
CARD
(1ND1C=1)
(1ND1C=1)
/
JC CARDS
(defleit)
1
r
/INDIC-
IEXIT
CARD
(1 EX 1T=1)
r
/
/
»
SC CARD 5
(deficit)
INDIC-IEXIT
CARD
(I EX IT = 1)
INDIC-IEXIT
CARD
( I E X I T= 1J
/P
/*
Figure U-2
Deck setups for various combinations
of reactants using program HAR03.
u-S
-------
GENERAL SYSTEM CARDS
[TITLE CARD]
Columns Variable
(1-40) TITLE
(41-43) N
(44-46) IPRNT
(47-49)
(53-59)
(60-66)
(67-73)
(74-80)
J CON
SCALE(1)
SCALE(2)
SCALE(3)
SCALE(4)
Description
Ufiits
To be used to label the
run
The number of sections in
the model
An indicator,
IPRNT=1 for printout of
matrices
IPRNT=0 to suppress the
printout of
matrices
An indicator used to desig-
nate the type of system being
modeled
JC0N=1 BOD
Deficit
JC0N=2 CBOD
NBOD
Deficit
JC0N=3 Chlorides
BOD
Deficit
JC0N=4 Chlorides
CBOD
NBOD
Deficit
Area scale factor ft^/ou*
Dispersion scale factor mi^/day/ou
Flow scale factor cfs/ou
Length scale factor ft/ou
Format
A40
13
13
13
F7'
F7
F7
F7
*note: ou stands for the "original units" in which the parameter
to be converted by the scale factor is expressed. The purpose
of the scale factors is to allow the user to input parameters
in units which are different from those specified on the
following pages. For instance, according.to this program
length should be input as feet. However, it may be more
convenient to enter it in miles and set SCALE(4) to 5280.
The program would then internally convert the length from
miles to feet
U-6
-------
[INTERFACE PARAMETER CARDS (2 per section; total number = 2*N)]
Columns
Variable
Description
Units
Format
(1-6)
AREA
Cross sectional area of
interface with section
ft2
F6
(7-12)
E
Dispersion coefficient of
first interface with section
mi2/day
F6
(13-18)
Q
Flow * across the first
interface with section
cf s
F6
(19-21)
IARAY
Section forming the first
interface with section
13
(22-27)
AREA
Cross sectional area of
second interface
ft2
F6
(28-33)
E
Dispersion coefficient of
mi 2/day
second interface
F6
(34-39)
Q
Flow across the second
interface
cf s
F6
(40-42)
IARAY
Section forming the second
interface
13
(43-48)
AREA
Cross sectional area of third
interface
ft2
F6
(49-54)
E
Dispersion coefficient of
third interface
mi2/day
F6
(55-60)
Q
Flow across the third inter-
face
cf s
F6
(61-63)
IARAY
Section forming the third
interface
13
(note: The second interface parameter card is identical to
above, but for the fourth, fifth and sixth interfaces.
If the section in question has an interface which
forms a boundary across which mass is transported
(that is, a boundary condition) interface parameters
must be input for it. To do this, input the appropriate
AREA, E and Q for the interface and input the section's
number as the IARAY.
It is only necessary to input the parameters for a
particular interface once. In other words after
inputting the parameters of the interface of section 1
with section 2, it is not necessary to input the
parameters of the interface of section 2 with section 1.)
* flow out of- a section is positive; flow into a section is
negative.
U-7
-------
[CHARACTERISTIC LENGTHS](1. card per section; 1 length for each interface;
therfore up to 6 lengths per card)
Columns
Variable
Descriotion
Units
Format
1-3
JI
First section which forms
an
interface with section I
13
4-13
LA(I,JI)
Length of section I with respect
to section JI*
feet
F10
14-16
JJ
Second section which forms
an
interface with section I.
13
17-26
LA(I.,JJ)
Length of section I with
respect to section JJ
feet
F10
27-29
JK
Third section which forms
an interface with section
I
13
30-39
LA(I,JK)
Length of section I with
respect to section JK
feet
F10
40-42
JL
Fourth section which forms
an interface with section
I
13
43-52
LA(I,J[L)
Length of section I with
respect to section JL
feet
F10
53-55
JM
Fifth section which forms
an interface with section
I
13
56-65
LA(I,JM)
Length of section I with
respect to section JM
feet
F10
66-68
JN
Sixth section which forms
an interface with section
I
13
69-78
LA(I,JN)
Length of section I with
respect to section JN
feet
F10
* note: the following figure.illustrates the interpretation of
characteristic lengths as used in this program.
figure U-l: Depiction of characteristic lengths
U-8
-------
where
@ - LA(I,JX) and LA(I,JI)
(2) - LA(I,JJ), LA(I,JM) LA(I,JL) and LA(I,JN)
As can be seen, any other than orthogonally shaped segments
would present a problem as to the measurement of LA. It is therefore
reiterated that oddly shaped segmentation schemes should be avoided
except when absolutely necessary to depict the natural geometry of
the system.
[DEPTH CARDS]
Columns Variable
I-10
II-20
H(l)
H(2)
Description
Depth of 1st section
Depth of 2nd section
Units
ft
ft
Format
F10
F10
H(N)
Depth of Nth section
ft
F10
[TEMPERATURE-VOLUME CARDS]
I-10
II-20
21-30
31-40
T(l) Water temperature of 1st section °C
V0L(1) Volume of 1st section 106ft3
T(2) Temperature of 2nd section °C
V0L(2) Volume of 2nd section 10^ft3
F10
F10
F10
F10
T(N) Temperature of Nth section
VOL(N) Volume of Nth section
°C
106ft3
F10
F10
U-9
-------
[ICON CARD]
Columns Variable Description Units Format
(1-2) ICON An indicator. If desired 12
at this point changes in
the E's or Q's which were
read in previously, may be
made using ICON
If IC0N=0: No changes are
called for and the program
would continue by reading
the INDIC card as described
on page U-ll
IF IC0N=1: Would allow a
new set of dispersion co-
efficients to be read as on
"New E cards" (see below)
If IC0N=2: Would allow a
new set of flows to be read
as described on page U-i» (new
Q Cards")
[NEW E CARDS (OPTIONAL-MUST BE PRECEDED BY ICON CARD=1)]
Columns Variable Description Units Format
(1-5) Revised E's For 1st Interface with first section Mi^/day F5
(6-10) " " 2nd " " " " " "
(11-15) " " 3rd
(16-20) " " 4th
(21-25) " " 5th
(26-30) " " 6th
(31-35) " " 1st " " second "
(36-40) " " 2nd " " " " "
(41-45) " " 3rd " " " "
(46-50) " " 4th
(51-55) " " 5th
(56-60) " " 6th
IN SIMILAR FASHION, REVISIONS ARE INPUTTED FOR THE THIRD TO Nth SECTIONS
note: These cards must be followed by another "ICON CARD" (as explained
cxWcve). If revisions are to be made in flow, ICON must equal 2,
and the new values can be read as on page U H. If no more revisions
are necessary, ICON must equal 0 and the program will proceed to read
the "INDIC CARD" as explained on page U U
U-10
-------
[NEW Q CARDS (OPTIONAL-MUST BE PRECEDED BY ICON CARD=2)]*
Columns
(1-5)
(6-10)
(11-15)
(16-20)
(21-25)
(26-30)
(31-35)
(36-40)
(41-45)
(46-50)
(51-55)
(56-60)
Variable
Revised Q's
Description
Units
For 1st Interface with First Section cfs
Format
F5
2nd
3rd
4th
5th
6th
1st
2nd
3rd
4th
5 th
6 th
Second
IN SIMILAR FASHION, REVISIONS ARE INPUTTED FOR THE THIRD TO.Nth SECTIONS
These cards must be followed by another
"ICON CARD" (as explained on pageU-vC- ).
If revisions are to be made in dispersion,
ICON must equal 1 and the new values can
be read as described on page U-iC. If no
more revisions are necessary, ICON must
equal 0 and the program will proceed to
read the "INDIC CARD" as explained below.,
[INDIC CARD]
Columns Variable
(1-2)
INDIC
Description
An indicator used to denote whether
subsequent DATA will be for Chloride,
BOD or D.O. deficit
Units
Format
12
INDIC =1 for Chloride
INDIC = 0 for BOD
INDIC = 1 for D.O. deficit
(3-4) IEXIT An indicator which can either continue 12
or terminate the run
If IEXIT>0 The run will Terminate
IEXIT<0 The run will Continue
*¦ on both the "New E Cards" and the "New Q Cards" the scale, factors input on the
title card (page U-6) apply.
U-ll
-------
SPECIFIC CONSTITUENT CARDS
[NUMBER OF BOUNDARY CONDITIONS]
Columns Variable
(1-2) NUMBC
[BOUNDARY CONDITIONS]
(1-8) BC(1)
(9-12)
ICOL(l)
BC (2)
ICOL(2)
Description Units Format
Number of boundary interfaces across
which mass is transported 12
Concentration of water quality mg/1 F8
parameter at first boundary
Section adjoining first boundary 14
Concentration of water quality
parameter at second boundary mg/1 F8
Section adjoining second boundary 14
(61-68)
(69-72)
BC (6)
ICOL(6)
Concentration of water quality
parameter at 6th boundary
Section adjoining 6th boundary
mg/1
F8
14
Proceed with cards until a total of NUMBC boundary conditions have been read
[RATE CARD #1]***
(1-5)
(6-10)
FAC1
K(l)*
K(2)
Temperature correction factor for first 1/day F5
reaction rate
First reaction rate for first section ** 1/day F5
for Chlorides =0.0
for CBOD =CB0D removal rate(Kr)
for NBOD =NB0D removal rate(Kn)
for D.O. deficit =»Reaeration rate
First reaction rate for second section 1/day F5
(76-80) K(15) First reaction rate for Nth section 1/day
Proceed until N reaction rates have been input
F5
* note: If all the rates are the same, only one value need be input as K(l)
and the remaining values input as zero. The program has an internal
test to see if K(2) is zero. If this is so, the program automatically
sets all values of K equal to K(l). Only one card should be input
when exercising this option and besides the value of K(l) the tempera-
ture correction factor should still be input.
** note: rates at 20° C
*** see page U-13
-------
[RATE CARD #2] (Optional-for D.O. deficit only)
Columns Variable Description
(1-5)
(6-10)
(11-15)
(16-20)
FAC2
FL
kd(D*
KD (2)
Temperature correction factor for
second reaction rate (deoxygenation
rate)
Ratio of ultimate to 5-day BOD
Deoxygetiation rate for section 1**
Deoxygenation rate for section 2
Units
1/ day
1/day
Format
F5
F5
F5
F5
(76-80) KD(14)
Deoxygenation rate for section 1A
(proceed until N rates have been input)
[PHOTOSYNTHESIS RATE] (Optional for D.0. deficit only)***
1/ day
(1-5)
(6-10)
PMK(l)
FIE (2)
Het photosynthesis rate for section 1 mg/l/day
Uet photosynthesis rate for section 2 mg/l/day
F5
F5
F5
(76-80) PMR(16) Uet photosynthesis rate for section lfc mg/l/day
(proceed until N rates have been input)
[BEflTHAL DEOXYGENATION RATE] (Optional for D.O. deficit only)***
(1-5) FAC3
(6-10)
(11-15)
BD(1)
BD(2)
Temperature correction factor for
b entb al denar.d
Benthal demand for section, 1**
Benthal demand for section 2
F5
gnv/M^ / day E 5
gm/M2/day F5
(76-80) BD(15)
Benthal demand for section 15
(proceed until H rates have been input)
gm/M^/day F5
* As described on the previous page, if all rates are equal only KD(1) need be
input and all other KD's can be left blank.
** rates at T=20° C
*** note: In an additive coupled run which calculates the deficit due to CB0D, the
deficit due to NEOD and then adds them it is unnecessary to input the
reaeration rates, the photosynthesis rate and the benthal demand with the
1IB0D part of the run. All such cards should be omitted rather than inputting
them a second time. This is illustrated in the example problem which follows.
U-13
-------
Columns Variable Description Units Format
(1-10) LOAD LOAD to a section ///day F10
(11-13) ISEC Section to which LOAD is applied 13
(note - Individual LOAD CARDS are needed for each LOAD input.
The last LOAD should be followed by a blank card which
indicates that the previous card was the last LOAD to be
input. If there are no loads merely input a blank card).
U-14
-------
EXAMPLE PROBLEM
The following example problem is primarily intended to
illustrate the input and output of HAR03. It should be noted
that several simplifications have been made which would never
be used in an actual modeling study and as such this example
should not be used as a strict guide to applications of the
program.
figure U-3: Overview of a fictitious tidal bay
The tidal bay depicted in figure U-3 has a dam situated at
its head end and opens onto a large estuary at its mouth. As shown
in figure TJ-4, the bay can be generally divided into a deep section
on the north and shallow flats on the south. An investigation of
the bays hydrodynamics has indicated that the fresh water flow entering
the system over the dam makes its way to the sea along the northern
coast and that relatively small amounts of Xvrater are transferred
between the flats and the deeper section.
U-15
-------
figure U-4: Depth readings (feet) for the tidal bay
At present, a large municipality discharges raw sewage
directly to the southwest corner of the bay resulting in critical
depression of dissolved oxygen levels as shown in figure U-5. One
question which might be addressed using HAR03 would be whether dis-
charging the waste into the deeper area to the north would signi-
ficantly improve the water quality?
Utilizing the available data, the segmentation scheme shown
in figure U-6 was developed. Data for the segments is given in
tables U-l and U-2.
The results of the analysis is shoxjn in figure U-7. As can
be seen, discharge into the northern part of the bay significantly
improves water quality by taking advantage of the additional
dilution, transport and mixing which exist there.
A printout of the input cards is given in figure U-8 for
the run in which the load is discharged in the northern part of
the bay. The corresponding output is also attached.
U-16
-------
Figure U 6: Segmentation scheme for the bay
U—17
-------
Interfaces
1-1
1-2
2-3
3-4
4-5
4-6
6-7
5-7
5-8
8-8
Area (ft^)
Dispersion (mile^/day)
Flow (cfs)
Boundary condition (rag/1)
Chlorides
CBOD
NBOD
DO deficit
150.
2.
2.
1.
39600.
0.5
150.
95040.
.5
150.
100320.
1.0
150.
105600. 36960.
1.0 0.5
243.
21200. 31680. 99000. 1584
0.5 0.5 1.5 1.5
243. 243.
1000
.5
.5
Table U-l: Interface data for run with load
to section 4
Sections
1
2
3
4
5
6
7
8
Temperature (°c)
21.0
22.0
22.0
22.0
22.0
•4.0
24.-0
22.0
Depth (feet)
12.0
13.0
15.0
20.0
22.0
4.0
5.0
30.0
Volume (10 feet^)
83.64
271.8
418.18
557.57
613.32
111.51
139.30
836.35
CB0D removal rate (1/day)
.3
.3
.3
.3
.3
.3
.3
.3
NBOD removal rate (1/day)
.3
.3
.3
.3
.3
.2
.3
.3
reaeration rate (1/day)
.24
.23
.22
.17
.15
.16
.18
.12
CB0D deoxygenation rate (1/day)
.1
.1
.1
.1
.1
.1
.1
.1
NBOD deoxygenation rate (1/day)
.1
.1
.1
. 1
.1
.1
.1
.1
CBOD load {illday)
100000.
NBOD load (#/day)
100000.
Benthal load (gm/M^/day)
.05
Algal load (mg/l/day)
.25
Waste load flow (cfs)
93.
Table U-2: Section data with load to section 4
-------
Figure U-7a: Results of model run with, load to
segment 6.
Figure U-7b: Results of model run with load to
segment 4.
-------
. 1L
/—-¦
'1.08 i, .1
in
100000.
"I.C'B .1
tzz:..
4...
¦.:.l
".35
1.02'. ,24 .23 .22 .17 .15 .16 .18 .12
.05
.25
".5
XT
1 .5
! i i ! i
-ii:
-! !
d
I
K>
o
-XPQooo. .
...I.. 047.. . 35
L. -5
iifcrr
rr?
::r;ippo:r
¦i—5
.h-JLt
4 1.
613.32 '
83.64 2?-
13. 15.
8 1.0
6 1.0
III.51 24.
271.8 22.
20. 22.
139.39 22.
418.18 22.
4. 5.
836.35
557.57
30.
2 1.
J-
/ 11.
9 1.0
5.1. 0..
4 1.0
3 1.0
2 1.0
7 1.
I-...
•' {.
, ! . : ' . .
¦m. ._401.5 243. 8
/J 1.2 IT" ___ 7 ~ *'
/"' 83U68'T5"""
/l05.601.0 243. 536.96 .5
100.321.0 150.
95.04 .5
150.
1-- 1^150....... 139,6.. 0
EXAMPLE RUN (LOAD TO SEC
P.-5 .150. 2_...
:iI0N 4) 8
i I i i !
i
1000.
1.
5280.
Mi
¦.in !?
^:;!i ; < i
i ;M ^ •'1 ' ; F
i • : ¦ •' :. •" : : s —¦
! j! =:;
j
: i
i ii"
Hi;
! !; ;
j !:
|; ;' i:; |Figure U-8: Input
f * j_f deck for exair
run with loac
section 4
Ml :ii;
i :
oar
: :
.. ;:MMm
!UHiUJ
M H \t
.< : .
F
HHiip
-------
EXAMPLE.RUN (LOAD TO SECTION <»>.
**SYSTEH PARAMETERS**
NUMBER OF SECTIONS =
JCON = 4
IPRNT
IPNCH
**SCALE FACTORS**
SCALE!1)
1000.000 SCALE(2)
1.000
SCALE(3)
1.000 SCALEI4)
5280.000
-------
INTERFACE AREA E. 0 INTERFACE AREA E ft INTERFACE £lREA E 8 _
ROvi-SEG (FT»»2) IHI»*2/D) (CFS) ROW-SEC (FT»»2) (Ht**2/D) (CFS) ROW-SEG (FT«<>2) (MI«»2/D) (CFS)
I I- 1) 0. ._ 0.0 -150.0 ; ( 1- 2) .... 39600. 0.500 150.0 ( 1- 0) 0. 0.0 0.0
J 1- 0) 0. 0.0 0.0 I 1- 0) 0. 0.0 0.0 ( 1- 0) 0. 0.0 0.0
J_2r_ 3• 95040. 0,500 150.0 (_2-1) 39600. 0.500 ri50.-0 L_2~_QJ 0, 0• Q .0...0
< 2- 0) 0. 0.0 0.0 I 2- 0) 0. 0.0 0.0 ( 2- 0) 0. 0.0 0.0
« 3- 4) 100320. 1.000 150.0 ( 3- 2) . .. 95040. _ ... 0.500 -150.0 ... _. ( 3- 0) 0. 0.0 . ... 0.0
( 3- 0> 0. 0.0 0.0 ( 3- 0) 0. 0.0 0.0 I 3- 0) 0. 0.0 ' 0.0
(4- 5) 105,600. 1.000 243.0 ... ( 4- 6) 36960. 0.500 ,0.0 ( 4- 3) 100320. 1.000 ...-150.0
( 4- 0) 0. 0.0 0.0 (4-0) 0. 0.0 0.0 ( 4- 0) 0. 0.0 0.0
( 5- a) 99000. Li500 243.0 L_5z._7.> 31680. 0.500 0.0 1 5- 4) 105600. 1.000 -243.0
( 5- 0) 0. 0.0 0.0 ( 5- 0) 0. 0.0 0.0 ( 5- oi 0. 0.0 "o'.c"
( 6- 7) 21200. 0.500 0.0 .._ ( 6- 4) 36960. 0.500 0.0 ( 6- 0) 0. 0.0 0.0
« 6- 0) 0. 0.0 0.0 ( 6- 0) 0. 0.0 0.0 ( 6- 0) 0. 0.0 0.0
( 7- 5) 31680. 0.500 0.0 ( 7- 6) 21200. 0.500 0.0 (7-0) _ 0. 0.0 0.0
< 7- 0) 0. 0.0 0.0 (7-0) 0. 0.0 0.0 ( 7-0) ~ 0. 0.0 0.0
« 8- 8) 158400. 1.500 243.0 ( 8- 5) 99000. 1.500 -243.0 (8-0) 0. 0.0 0.0. _
• 8- 0) 0. 0.0 0.0 ( 8- 0) 0. 0.0 0.0 (8-0) ~ 0. 0.0 0.0
-------
CHARAC. LENGTHS OF SEGMENTS 1 FT>
_ interface" length
1- I 5280.
2— 1 5280.
3- 2 5280.
4-
3
5280.
4- 5
5280.
4- 6
5280.
4- 0
0.
4- 0
0.
4- 0
0.
5-
4
5280.
5- 8
5280.
5- 7
5280.
5- 0
0.
5- 0
0.
5- 0
0.
6-
4
5280.
6- 7
5280.
6- 0
0.
6- 0
0.
6- 0
0.
6- 0
0.
7-
5
5280.
7- 6
5280.
7- 0
0.
7- 0
0.
7- 0
0.
7- 0
0.
8-
5
5280.
8- 8
5280.
8- 0
0.
8- 0
0.
8- 0
0.
8- 0
0.
INTERFACE LENGTH INTERFACE LENGTH INTERFACE LENGTH INTERFACE LENGTH INTERFACE LENGTH
1- 2 5280. ..1- 0 0. I- 0 0. _ 1- 0 0. __ 1- 0 _ 0.
2- 3 5280. 2- 0 0. 2-0 0. 2-0 0. 2- 0 0.
3- 4 5280. 3- 0 0. 3- 0 0. 3- 0 0. 3- 0 0.
-------
SECTION * HAS AN EXCESS FLOW OF CPS
-------
VAtUES OF ALPHA ANO EPRtH
THTftFC EPR1K ALPHA INTRFC EPRIH ALPHA INTRFC GPRt H ALPHA INTRFC EPSUM ALPHA tNTRFC EPRIM ALPHA INTRFC EPRIH ALPHA
I MOD I IMGOI CHGDt (MGO) _ I.HGDI IHGOI
I- 1 0. I.COO 1- 2 782. 0,500 1- D 0. 0.0 1- 0 0. 0.0 1- C 0. C.C I- 0 D. 0.0
i- 3 1877. 0.500 2- 1 782. 0.500 .2- 0 0. 0.0 2- 0 0. 0.0 2- 0 0. 0.0 _ 2- 0 _ 0. 0.0
3- * 3963. 0.500 3- 2 1877. 0.500 3- 0 0. 0.0 3- 0 0. 0.0 3- 0 0. 0.0 3-0 0. 0.0
i- 5 4171. 0.500 6 730. 0.500 3 3963. 0.500
-------
SECTION TEMPERATURE VOLUME QEPTH.
(C) <10»»6GAL) (FT)
1
21.00
625.63
12.00
'J:
2
22.00
2033.06
13.00
3
22.00
3127.99
15.00
4
22.00
4170.62
20.00
¦
5
22.00
4587.63
22.00
6
24.00
834.09
4.00
7
2*1.00
1042.64
5.00
8
22.00
6255.89
30.00
-------
SEGMENT CL BOUNDARY COND IT ION(MG/L)
1 0.0
8 1000.00
-------
SECTION
-------
SECTION CHLORIDESIMG/LI
1 755.94 3
2 855.837
. 3 901.211
4 923.532
5 957.424
6 932. 199
7 947.311
8 __ 983.405_ _
-------
SEGMENT BOD BOUNDARY CONDITtON1HG/L)
i 2.00
8 0.50
-------
BOO LOAD > 100000. POUNOS/OAY FOR SECTION
-------
THE TEHPERATURE CORRECTION FACTOR FOR THE BOD REMOVAL RATE - 1.047
SECTION BOD REM RATE BOO WASTE LOAD BOO BOUNDARY LOAD
~TEMP CORR^ (MGD»MG/l> (MGD»MG/L>
1 0.366 0.0 193.890
2 0.384 0.0 0.0
3 0:384 0.0 0.0
4 0.384 11990.406 _ 0.0 _
5 0.384 0.0 0.0
6 0.421 0.0 0.0
7 0.421 0.0 0.0
8 0.384 CM) 4653.2 89
-------
SECTION BOD I HG/L)
1 0.899
2 1.035
3 1.538
* 2.257 _
5 1.112
6 1.335
7 0.846
8 0.638
J
-------
SEGMENT P6F BOUNOARY CONO8T10NIHG/L)
1 O.SO
8 _ ' 0.50
-------
THEJEMPERATURE .CORRECTION FACTOR FOR THE REAERATION RATE * 1.024
The OOD CORRECTION FACTOR (FL) = 1.000
THE TEMPERATURE CORRECTION FACTOR FOR THE BENTHAL OEMAND = 1.080
THE TEMP CORR FACTOR FOR KD » 1.047
SECTION OEOX RATE REAER RATE BOUND COND LOAO TOTAL LOAD PMR BD
® T e M P CORR» »TEHP CORR» (MGb«MG/L> (MGO*MG/L) MG/L/D »TEMP C(5rR«
I 0.366 ..0.246 _ _ 48.472 _ 254.649 0.0 0.0
2 0.384 0.241 0.0 807.254 0.0 0.0
3.... Oo 384 0.231 0.0 . 1845.740 0.0 0.0
4 0.384 0.178 0.0 3611.526 0.0 0.0
5 0.384 0.157 0.0 1957. 236 CU.0 0.0
6 0.421 0.176 0.0 514.864 0.0 0.068
7 0.421 0.198 0.0 110.375 _ 0.250 0.0
B 0.384 0.126 4653.289 6184.867 0.0 0.0
-------
SECTION OEFICIUMG/L)
1 U643
2 1 * 8 SB
3 2.013
4 1.9B7
•> 1.4b9
_6 £-000
7 1.492
8 0.930
-------
SECTION CHLORIDES TEMPERATURE DEFIC IT SATUPmTION p.g.
1 755.94 21.00 1.64 8.78 7.14
2 855.84 22.00 1.86 8.60 6.74
3 901.21 22.00 2.01 8.60 " 6.58
4 923^53 22.0,0 1.99 8.59 6.61
5 957.42 22.00 1.47 8.59 7.12
6 932.20 24.00 _ 2.00 8.26 6.26
7 947.31 24.00 1.49 8.26 6.77
8 983.40 22.00 0.93 8.59 7.66
-------
-------
BOO LOAD
100000
POUNDS/DAY FOR SECTION
-------
THE TEMPERATURE CORRECTION FACTOR FOR THE BOD REMOVAL RATE » l.OBO
SECTION BOO REM RATE 800 HASTE LOAO BOO 80UN0ARV LOAD
*TEMP CORR» (HGD»MG/l_) IMGO*MG/L)
1 0.10B 0.0 193.890
2 Q. U.7 0 .,0 0 ..0
3 0.117 0.0 0.0
4 0.117 11990.406 0.0
5 0.117 0.0 0.0
. 6 ... 0.136 0.0 „ , 0.0
7 0.136 0.0 0.0
8 0.117 Q_iO (M)
-------
SECTION BOO { MC/L )
1 2.369
2 2.635
3 3.098
A __ 3.615
5 1.832
6 2.731
T 1.930
B 0.681
-------
SEGMENT DEF BOUNDARY CONDITIONIHG/L)
1 0.0
-------
THE_TEHPERATURE CORRECTION FACTOR FOR THE REASRATION RATE » 1.024
THE BOO CORRECTION FACTOR (FLl = 1.000
THE TEMPERATURE CORRECT ION FACTOR FOR THE 8ENTHAL DEMAND = 1.080
THE TEMP CORR FACTDR FOR KD = 1.060
SECTION DEOX RATE REAER RATE BOUND COND LOAD TOTAL LOAD PMR BO.
~TEMP CORR* »TEMP"C0RR* IKGO«MG/L> (MGO*KG/L> MG/L/O »TEHP CORR»
1 0.108 0.246 0.0 160.041 . 0.0 0.0
2 0.117 0.241 0.0 624.963 O.D O.O
t 0.1J.J 0.231 .. _ D.D . 1130.«.
-------
SECTION DEFJLC1T ( MG/L )
1 0.965
2 _ 1.077
3 1.075
4 0.983
5 0.639
6 1.085
7 0.902
8 0.268
-------
RESULTS FOR CQKft I^EO
-------
JOB CONTROL LANGUAGE
The following JCL statements apply to users of HAR03 at OSI
computer facilities:
//INTHAR03 JOB (ACNT,INT,2,5),PROGNAME
//JOBLIB DD UNIT=3330,VOL=SER=BCS025,DISP=SHR,DSN=CHAPRA.xxxxx
//STEP1 EXEC PGM=MODEL,REGION=nnnK
//GO.FT07F001 DD SYSOUT=B
//GO.FT06F001 DD SYSOUT=A
//GO.FT05F001 DD *
Place data cards here
/*
Note:
xxxxx - HA100
for
100
segment
model
HA200
for
200
segment model
HA50
for
50
segment
model
nnn - 76 K
for
50
segment
model
148 K
for
100
segment
model
410 K
for
200
segment
model
INT - OSI user initials
ACNT - OSI account
U-46
-------
Restrictions
1) HAR03 is presently limited to a system of 200 segments
2) A section may have only one interface act as a boundary
3) Each section may have up to six interfaces
U-47
-------
APPENDIX
-------
References
Introduction
1. Streeter,,H.W. and Phelps, E.B., A Study of the Pollution and
Natural Purification of the Ohio River, III. Factors
Concerned in the Phenomena of Oxidation and Reaeration.
U.S. Public Health Service, Public Health Bulletin
Number 146, February, 1925, 75 pp. Reprinted by
U.S. Department of Health, Education and Welfare, PHS, 1958.
2. Thomann, R.V. Mathematical Model for Dissolved Oxygen
Journal of the Sanitary Engineering Division, ASCE, Vol.89
No. SA5, Oct. 1963, pp. 1-30.
Theory
1. Thomann, R.V. same as 2 above.
2. Thomann, R.V. Systems Analysis and Water Quality Management,
"Environmental Research and Applications, Inc., New York,
1971.
3. Arons, A.B. and Stommel, H„ , "A Mixing Length Theory of Tidal Flushing",
Transactions, American Geophysical Union, Vol. 32, No. 3,
June 1951, pp. 419-421.
4. Stommel, K., "Computation of Pollution iri a Vertically Mixed Estuary",
Sewage and Industrial Wastes, Vol. 25, No. 9, September, 1953,
pp. 1065-1071.
5. O'Connor, D.J., "Estuarine Distribution of Non-Conservative Substances",
Journal of the Sanitary Engineering Division, ASCE, Vol. 91,
No. SA1, February, 1965, pp. 23-42.
6. Thomann, R.V. System Analysis and Water Quality Management, p. 180
-------
The Computer Program
1) Development of Water Quality Model of Boston Harbor Prepared
for Massachusetts Water Resources Commission by Hydroscience, Inc.,
Westwood, N.J., July, 1971, 173 pp.
2) HAR01: A Steady State Estuarine Water Quality Computer Model
prepared by Data Systems Branch TJ.S.E.P.A., Region II, N.Y., N.Y.
3) Thomann, R.V. Systems Analysis and Water Quality Management p.96 & 103
4) Committee on Sanitary Engineering Research, Solubility of Atmospheric
Oxygen in Water, J. Sanit. Eng. Div. , Am Soc. Civil Engrs., SA 4,
41 (1960).
-------
LIST OF FIGURES
Introduction
figure page
1-1 D.O. sag generated by the Streeter - Phelps equation 1-3
Theory . . . . .
figure page
T-l One dimensional Estuary of N segments T-3
T-2 Schematic for a coupled system of two reactants T-10
T-3 Schematic of generalized two-dimensional system T-13
T-4 Simple schematic of HAR03 computer program T-16
Computer Program ....
figure page
C-l Flow chart of the main program of HAR03 C-3
Users Manual ....
figure page
U-l Structure of data deck for tt.AJR.03 U-4
U-2 Deck setups for various combinations of reactants using
program HAR.03 U-5
U-2a Depiction of characteristic lengths U-8
U-3 Overview of fictitious tidal bay U-15
U-4 Depth readings for the tidal bay U-16
U-5 Dissolved oxygen concentrations U-17
U-6 Segmentation scheme for the bay U-17
U-7a Results of model run with load to segment 6 U-19
U-7b Results of model run with load to segment 4 U-19
U-8 Input deck for example run with load to section 4 U-20
LIST OF TABLES
U-l Interface data for run with load to section 4 U-18
U-2 Section data with load to section 4 U-18
-------
NOMENCLATURE
SYMBOL
DEFINITION
UNITS *
Documentation
Computer
Program
AREA (J) or
ARRAY (I,JJ),
where
JJ=1,4,7,10,13,16
J=l, 6
1=1,N
Cross-sectional area
[A]
AC(I,J)
1=1, N
J=1,N
Single system matrix. A matrix which
represents a single substance which
decays with first order kinetics (see
L3/T
[A]-l
AC (I,J)
1=1,N
Single system response matrix. The r(M/L3)/
inverse of [A], this matrix represents (M/T)]
J=1,N
the response of the system to the input
of a load of a single substance decaying
with first order kinetics (see page
T-10)
[AB] or AB
AB(I, J)
1=1,N
J=1,N
The system matrix. Since matrices [A]
and [B] differ only in a reaction term
in their diagonals, a matrix without
such a term has been formulated and
is called the system matrix (see page
T-15)
L3/T
AC
AC(I,J)
1=1,N
J=1,N
The specific constituent matrix. When
the appropriate reaction term is added
to the diagonal of the system matrix
the matrix for a specific constituent
is generated. In the theory they were
called [A] or [B].
As well, the term AC is used in the
computer program to represent the
matrix after it has been inverted.
3 .
LJ/T
[(M/L3)/
(M/T)]
[B]
AC(I,J)
1=1,N
J=1,6
A matrix which represents the second
substance of coupled set of constituents
which react with first order kinetics
(see page T-ll)
M/L3
[BF1
AC(I,J)
1=1, N
J=1,6
The inverse of [B], this matrix represents
the response of the system to the input
of a load of the second substance of a
coupled set of constituents.
[(M/L3/
(M/T)]
BD
BD(.I)
1=1,N
Benthal oxygen demand
M/L2/T
BOD
Biochemical oxygen demand
M/L3
concentration of a constituent
M/L3
* units; T=time, L=length, M=mass
-------
SYMBOL
Documentation Computer
Program
DEFINITION
UNITS
[C] 1
The total system response matrix relates
[(M/L3/
input of the first substance in a coupled
set of reactants to the response of the
(M/T)]
second constituent.
CBOD
Component of BOD due primarily to oxygen
M/L3
utilized by suprophytic organisms as they
utilize the carbonaceous matter in a waste
D
Dissolved oxygen deficit = the saturation
value minus the actual concentration of
M/L
dissolved oxygen = D0S - DO
D0
The initial concentration of dissolved oxygen
deficit at the point introduction of x^aste
M/L3
into a stream.
DO
Concentration of dissolved oxygen
M/L3
D0 =
CS
Saturation concentration of dissolved
oxygen
M/L3
E(J)
ARRAY(I,
J J KL)
where
1=1 $
E
J=1 ,6
Longitudinal dispersion coefficient
l2/t
JJ=1,4,7
,10,
13,16
E1
EPRIM(J)
Bulk dispersion coefficient
J=1 6
=EA/1
l3/t
H
H
depth
L
k or K
First order reaction coefficient, (also
1/T
used to designate a section)
Ka
KA .
Reaeration rate (first order)
i/t
Kd
KD
Deoxygenation rate (first order)
1/T
Kr
KA
BOD removal rate (first order)
1/T
-------
SYMBOL DEFINITION
Documentation Computer
Program
UNITS
section.length
L0 Initial concentration of BOD at point of
introduction of a waste into a stream. M/L~
NBOD Nitrogenous component of BOD primarily M/L"^
caused be autotrophic bacteria utilizing
nitrogen in waste.
PMR PMR Net algal oxygen rate = photosynthesis M/L^/T
effect - respiration effect of algae on
oxyeen
Q Q(J) net advective flow L^/T
ARRAY(I,JJ+2)
where
J=l»6
1=1, N
JJ=1,4, 7,10,
13,16
t time
T(I) temperature
advective velocity L/T
-------
SYMBOL
Documentation
Computer
Program
DEFINITION
vol VOL Volume
W W Waste load of a constituent
Wb Deficit waste load
W„ BOD waste load
x Distance
Y Y Dissolved oxygen deficit waste load
a ALPHA weighting coefficient used in the
approximation scheme of the model.
3 1-ALPHA weighting coefficient = 1-a
-------
Acknowledgements
The material presented in this report is based on extensions and
revisions of the first version and that original work is gratefully
acknowledged. Steve Chapra has moved from E.P.A. and is presently an
Environmental Engineer for the Great Lakes Environmental Research
Laboratories at Ann Arbor, Michigan.
Thanks must be extended to Richard Winfield of the Manhattan College
Department of Environmental Engineering and Science and to Robert Braster,
Salvatore Nolfo and Kevin Bricke of E.P.A. for their thoughtful reviews
of the material contained in this report.
The typing of the report was very patiently and professionally
handled by Ms. Maryann LaBarbera.
G. A. N
------- |