UNITED STATES ENVIRONMENTAL PROTECTION AGENCY
REGION II
26 FEDERAL PLAZA
NEW YORK. NEW YORK 10007
documentation for
HAR03 J
"A computer program for the modeling of water quality
parameters in steady state multi-dimensional
natural aquatic systems"
prepared by
Steven Chapra
Systems Analysis Sect
Data Systems Branch
April, 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-rdimensional analysis, 		T-13
Application pf 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
Restrictions					U-46
APPENDIX
References
List pt figures and tables
Nomenclature
ACKNOWLEDGEMENTS
i

-------
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 K of core to compile
and was written for a Fortran IV G 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.
At present, HAR03 can accommodate up to 100
completely mixed segments.
ii

-------
INTRODUCTION

-------
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 qf 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 terms
1-2

-------
of first order kinetics. The result is what is now called the
Streeter-Phelps equation which in its basic form is:
D = ^	 ( e"Kr x/u - e~Ka x/u ) ] L + D e~Ka x/u 	(1-1)
•K - K	o o
a r
where:
D=dissolved oxygen deficit= D0R -DO
D0g =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
Kd=Deoxygenation rate
Ka=Reaeration rate
U=stream velocity
x=distance downstream from point of introduction of waste
The result of the Streeter-Phelps equation is called the "D.O.
sag" and is illustrated in figure I—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 between 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
Thomann for the Delaware2 which has proven itself to be particularly
effective 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 temperature, 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 mathematics the concentrations of the variables in each
volutnn 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 water to which this program can be applied.
However, 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, HAR03 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 iis 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-l

-------
One dimensional analysis of a single substance...
Thomann's approach 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:
— + U— = E-— - kc + (sources - sinks) 	(T-l)
3t 3x 3x2
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 systejn 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:
3 c". 3c" —32 "c —	• t x
— + Ufi— = ^kc + (sources - sinks) 	(T—2)
3tt 8x 3x2
where
the bar above a parameter or variable designates
that the value is averaged over a tidal cycle
Uj = 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 pf 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


n
1
z
J

k—1
k

n-1
n+1
Figure T-l:- Estaary. of a 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)
3ct
Vol
8ci
3 c,
k3t +	~ EkAk1k- VVk* wk	(t-3)
where
k= a subscript which denotes that the variable or parameter is
an average for the segment k.
Qfk= 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 Sc^/St would equal zero and equation T-3
could be written as:
0 -	" ^fkH— " *kVolkck + \	w-4)
dxz	dx
Each term containing a differential in equation T-4 can be
approximated mathematically, as:
d2ck
= Ek-l,k  + Ek,k+1 	
i,j l.+l.
i J
T-4

-------

= 1 - a

(T-8)
1.
J
As well, the above weighting coefficients
can be used to insure meaningful results. It
can be shown that for all positive solutions
ct>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 shqwn that a numerical
dispersion effect is implicit in this particular
difference scheme which.can be estimated by
^numerical * U1 <«" V»	(«)
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):
k-l,k (ck-l"Qk) + E k,k+l (ck+rck>
Qk-l,k (ak-l,kck-l + 6k-l,kQk)
" Qk,k+1 (ak,k+lck + 8k,k+l ck+i>
" KkVolkCk + Wk	
By grouping terms
(~Qk-l,kak-l,k"E k-l,k) pk-l
+ (Qk,k+lak,k+l " Qk-l,k3k-l,k + E k-l,k + E'k,k+1 + VolkKk)ck
+ (Qk,k+lBk,k+l " E k,k+l) ck+l = Wk • • • • ¦	CT-10)
Letting
ak,k-l = ~Qk-l,kqk-l,k ~ E k-l,k * *		CT-11)
ak,k = Qk,k+lak,k+l"Qk-l,kek-l,k + E k-l,k + E k,k+l + VolkKk * •(T"12)
ak,k+1 = Qk,k+]A,k+l ~E k,k+l	(T-13)
T-t6

-------
The general equation for the kth segment is therefore
ak,k-l ck_i + 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 pl + a12 c2 = W1 + (QOl «01 + E Ol) cD .. • • • • ••• • • • (T-15)
or
all C1 + a12c2 = W'l				 • (T-16)
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):
an-l,nCn-l + an,ncn ~ Wn + ^-<^n,n+l^n,n+l + ^ n.n+l^Si+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
llcl + a12c2 + 0
21ci + a22c2 + a23c3
0 +
0 +
a32C2 +a33C3 +a34C4+
+ 0
+ 0
+ 0
w.
w„
w„
+ . . . + 0 + an>n+1cn+1 + anncn = Wn
T-8

-------
At this point there are ti equations with n 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
a32 a33 a34

-------
Besides allowing solution of equation (T-21)s equation T-22 offers
an additional feature. The elements of the matrix [A]--*-, which is called
the single system response matrix, represent the unit response of a particular
section to the addition of a unit load to another section. For instance,
the element &2-\ of [A]-"'" 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...
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 k^ is referred to as the BOD removal rate,
k2 is a rate of reaeration from the atmosphere and k^ is a deoxygenation
rate. In the case where there is no loss of BOD via sedimentation or
another non-oxygen demanding removal, k-^ would equal k3.
The equations describing the conservation of mass 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
V
Figure T-2: Schematic for a coupled
system of two reactants
-b ) + E
k,k+l^k+l ^k^
Q (a b + g	b )
k,k+1 k,k+l k k,k+1 k+1
(T-23)
T-10

-------
where
= BOD in section k
= dissolved oxygen deficit in section k
Kk2 = reaeration rate in segment k
= deoxygenation rate in segment k
= waste load or source of dissolved oxygen deficit to
section k
T-21:
where
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 rather than
Vol Kfc!
(c) = solution vector to the single system (in this
case BOD).
The solution is therefore
(b) = [B]_1[Vol K3](c) +[B]'l(Wb) 	 (T-25)
In this program, three sources of dissolved oxygen deficit will be
considered:
W, = -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
isubstituting. equation 1V22'. into equation T-25 the following rela-
tionship results
(b) = [B]-1 (Vol K^) [A]"1 (Wc) + [B]"l (Wb)	(T-26a)
or after expressing the terra (Vol K^) as a diagonal matrix and multiplying
the matrixes the following results
(b) = IC]-1 (Wc) + [B]"1 (Wj,) 			(T-26b)
As can be seen [C]~l is a particularly useful entity in that it directly
relates the deficit in the body of water to the waste load of BOD being
discharged. It is called the total system response matrix.
T-12

-------
Multi-dimensional analysis.
An extension of the previous analysis:to more than one-dimension is
relatively straightforward.

k

k
i
k

k

Figure T-3: schematic of generalized two dimensional system.
Figures T-3 illustrates "a two dimensional system. Using the convention
that flow entering a section is negative and flow out of a section is positive,
an equation can be written representing the conservation of mass in segment i:
V,
.-4^ ~ 0 = I -Q., (a., c.+B.iC, )+E* .. (ci-c . ) -V .K • C .+W .. .
i "t	t. xik ik l ik k' ik k i l i i l
k
. . . . (T-27)
It should be noted that the generalization of the advection term
can be stated in this way due to the fact that
~^ikaik ~ ^ki^ki		
and	• 			 . . 		
(T-27) may be rewritten as
(T-28)
(T-29)
iici + £ aikci=Wi "(1=1,2...n)
where
_ E
+ E ^) + ViK±
aii ~ k wik"ik T * ik
aik ^ik^ik~E ik
(T-30)
(T-31)
(T-32)
T-13

-------
For boundaries where flow enters the section with a concentration c :
b
aii " k (Qilc0'ik+E'ik)+VlKl"H,li6"+E'ii	(T"33)
Wi = Wi +  + Vi +	+ E'ii			(T-35)
«;= wi +  1 - E /Q .......
to insure positive results also applies.
(T-39)
T-1A

-------
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) = [CI"1 (Wc) + [B]"1 (Wb) . . . . . 		 (T-26b)
For HAR03, equation T-26b will not be used due.tep 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 [A] 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] or IB] 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, [C]-l.
T-15

-------
With this as a background, a.very simple representation of the
major steps of HAR03 can be presented:
START
step 1
Read data which is applicable to the
system as a whole (e.g. hydrodynamic and
	physical data)		
step 2
Construct the system matrix.
AB
¦ step 3
Read data which is applicable to the
specific constituent which is being model-
ed (e.g. rates, loads, boundary
BAk
step 4
Construction of the forcing functions
step 5
Construction of the specific constituent
matrix
step
Calculation of the concentration of the
constituent"by inverting the matrix, AC and
multiplying it by the forcing function
Figure T-^: 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 IAB], 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
®ii Ai4
E'i:j = ' 3 ^ 	(T-41)

2) a value of alpha is calculated to adjust the advective
approximation for segments of unequal length.
for flow across the	1^
interface and into the
a ji 1± +1. . . . 		(T-42)
section	1 J
(Q< o)
for flow across the	_	lj
"face fnd °Ut	" 1J	li + 1,	(T-43)
of the section	1 J
(Q>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 V * w± + (E'u ~ Q±i k±i) -<*b	(T~36>
for Q < o W.' = W •_ + (E'±i - Qii (X .) c b	(T-34)
As can be seen, the forcing function consists of two basic
parts:
1) Boundary condition forcing function
of the form, [E'^ - Qii ( a i± or8ii)]cb
the bulk dispersion across the interface
between the section and the boundary,
the flow across that interface
or = assumed initially to be 0.5 and then
subsequently tested as in equation
T-39 and if necessary recalculated
according to equation T-44.
cb = the concentration across the boundary.
Aside from this, this part of the forcing
function does not vary as to the constituent
being modeled.
where
E'ii -
Qi:
T-7I8

-------
2) Other forcing functions
of the form,
for BOD or chlorides
y. = any load to the section (e.g. waste loads, runoff
1	loads, etc.)
for DO deficit
W£ = -PMR + BD*H. + loads
as explained in equation T-26-
7-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 arid to improve its clarity and ease of use.
In its present form it is programmed for the IBM 370/155 where it
requires 184 K of storage and will run on both a G and H level compiler'.
On the G Level it now requires approximately .26 minutes to compile.
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 time.
Number of segments
8
30
45
64
execution time (minutes)
.02
.10
.27
.60
At present, HAR03 • can handle up to 100 segments.
The description of the program which follows consists of:
a)	a flow chart of the main program
b)	brief descriptions of the main program and each of the subroutines
c)	a listing of the source deck
C—2

-------
START
SUBROUTINE DATA
Input General System Parameters-*which do not depend on the
constituent being modeled
I	SUBROUTINE REVI (optional)
| Used to revise the parameters E and/or Q without altering
I 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, JAB]
READ INDIC & IEXIT
INDIC is used to designate the type
of constitaefit 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)
' Writes matrix [AC]
SUBROUTINE MATN
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
C-3

-------
DESCRIPTION OP THE MAIM PROGRAM AND EACH OF THE SUBROUTINES
The Main Pr pgrara-HAR03¦..
HAKO3 has been programmed using a modular approach. Wherever
possible, ar. attempt has been made to give e&cj:i 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-l certain details have been
omitted and if the user requires isote detail rie 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 feeing modeled. These
include system parameters (number of sections, indicators, etc.) and
certain physical parameters such as cross-sectional area, net advective
flaw, etc, , wliich are only input once when modeling a water body. -Triese
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~j initialise certain r&atrices to zero
3) read and writs interface parameters AKEA, E and Q and place
them into the storage oiatrix ARRAY.
i) read and write the characteristic, lengths (LA)
5)	read the geaeral section parameters depth (H), temperature (T>,
and volume (VOL)
6)	Apply the following conversion factors:
QCKED) = Q fCFS) * .6463 MGO/cfs
VOL(MG) = VOL (105ft3) * 7.48 gal/ft3
C-4

-------
At tliis 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)
ARRAY(I,J)
ARRAY (I, J + 1)
ARRAY(I,J + 2)
J
contains the section number
of the section which forms
the interface with section I.
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
' ¦¦ ¦ 			« i *
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, E'
According to the theory:
where
E'i i = Ei,j Ai,.j
,J 1- •
E^j = the dispersion coefficient across the interface
A,- i = the cross-sectional area of the interface
» J
li,j = the average length of the adjoining
sections = (1^ + lj)/2
The program calculates E1 (EPRIM) as follows:
EPRIM(J)=ARRAY(I,JJ)*ARRAY(I,JJ+1)*417.166/(L(K,I)+L(I,K))
C-6

-------
where
ARRAY(I.JJ) = A. . (ft2)
1 J J
ARRAY(I,JJ+1) = Ei}j(mi2/day)
L(I,K) = 1±
L(K,I) = lk
417.1166 - conversion.factor which is computed as follows:
Original Units	Conversion Factor	New Units
? 2	•>
ft -mile	(5280 ft)z x 7.481 gal x 2 X_MG	MGD
day'"f t"-':2	(mile.) ^	f. t3	10^~gal
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 riot 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 [ABJ 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 HARP3	which is analogous to	
Chlorides	any conservative substance
BOD	a single constituent reacting with
first order kinetics
D.O. deficit	the second constituent of a coupled,
feed-forward system reacting with first
order kinetics.
BCC performs the following tasks:
1)	Reads the number of boundary conditions
2)	If there are none it writes the message
ZERO BOUNDARY 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)*ARRAY(I,Jj+2))
for Q > 0
FRCBC(I)=8.34*BC(KLAST)*(EPRIM(J)-(1.0-ALPHA(J)*ARRAY(I,JJ+2))
C-8

-------
where
FRCBC(I) = forcing function due to the boundary condition
BC = boundary concentration (mg/1)
EPRIM = E* (HGD)
ARRAY(I,JJ+2) = Q (MGD)
8.34 = conversion factor of (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
Kj - K20 *8 (T"20)
where some typical values are 3
rate	9
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
C-9

-------
4) for D.O. deficit the forcing function is completed by
first reading the deoxygenation rate {K
-------
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) = (l.-.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 cole
in the coordination of a particular configuration of
reactants.
C-ll

-------
£***#*
£****«
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
COMPUTER PROGRAM HAR03
HAR03 IS A COMPUTER PROGRAM WHICH CAN BE USED TO MODEL THE
STEADY-STATE DISTRIBUTION OF WATER QUALITY VARIABLES FOR
MULTI-DIMENSIONAL BODIES OF
THE PROGRAM IS BASED ON THE
OF MASS AND UP TO TWO WATER
WATER. THE TECHNIQUE UNDERLYING
PRINCIPLE OF THE CONSERVATION
QUALITY VARIABLES REACTING
IN A FEED
MODELED
FORWARD FASHION WITH FIRST ORDER KINETICS MAY BE
HAR03 HAS BEEN DESIGNED FOR
APPROXIMATELY	K OF CORE
G LEVEL COMPILER
AN IBM 370/155 COMPUTER, REQUIRES
AND WAS WRITTEN FOR A FORTRAN IV
THIS VERSION
SEGMENTS
CAN ACCOMODATE UP TO 100 COMPLETELY MIXEO
COMMON/ALGAE/NX,N»iNDICtICON t MMM,NMK
COMMON/BASS/MX, I PRNT t.FRCBC ( 100)
1 FORMAT(12)
12 FORMAT(CE10.4)
20 FORMAT(312)
MX = 5
NX = 6
NMK=0
CALL DATA
25	READ(MX ,1)ICON
IF(I CON)57,57,26
26	CALL REVI
GO TO 25
57 CONTINUE
CALL FLOW
CALL SETUP
MMM= 1
93 READ(MX,20)INDIC» I EX I T
IF(IEXIT)8,8,1000
8 CALL BCC
CALL RATE
IF(IPRNT)10,10,1003
1003 CALL WHITE
10	CALL MATN
IF(IPRNT)11,11,1013
1013 CALL WRITE
11	CALL STORE
GO TO 93
1000 STOP
END
HAR03000
HAR0300I
HARO 3002
HAR03003
HAR0300*
HAR03005
HARO 3006
HAR03007
HAR03008
HARO 3009
HAR03010
HAR0301I
HARO 3012
HAR03013
HAR030 14-
HAR03015
HAR03016
HARO 3017
HAR03018
HAR03019
HAR03020
HAR0302I
HAR03022
HAR03023
HAR03024
HAR03025
HAR03026
HAR0302 7
HAR03028
HAR03029
HAR030 30
H A R 0 30 31
HAR03032
HAR03033
HAR03034
HAR03035
HAR03036
HAR03037
HAR03038
HAR03039
HAR03040
HAR03041
HAR0304-2
HAR03043
HAR03044
HAR03045
HAR03046
HAR03047
HARO 3048
HAR03049
HAR03050
HAR03051
HAR03052
HAR0305 3
HAR0305<»

-------
INPUTS DATA WHICH
CONSTITUENT BEING
HYDRODYNAMIC DATA
WOULD NOT VARY WITH THE
MODELED, I.E. PRIMARILY
PARTICULAR
PHYSICAL AND
SUBROUTINE DATA
£44444444444444444444444444444*4**4*4 4 *******4*4 (-'*4444****44444*4444444
£44444 4444 4444 44444 444444 4444**4*********4*44444**4*4**************444 4
C
C	SUBROUTINE OATA
0
C
C
C
c
£44 44 4444444 44 **444** 4444444444*********************************** *4444
£44444444444444*4444444444444444444444444*********4********************
REAL LA( 100 , iOO)
DIMENSION AREA(6),E(6),Q(6),TITLE<10)
COMMON/ALGAE/NX» N,1NDIC»I CON,MM*,NMK
COMMON/BASS/MX,IPRNT,FRCBC(100)
COMMON/COHUE/LA,SCALE(4),ARRAY{100,L8),IARAY(100,6>
COMMON/P IKE/JC0N,H( 100) ,T( 100) ,VQL (100),AB(100,100)
FORMAT! 8 F10 . 0 )
FORMAT!10A4,4I3,4E7.0)
31
100
101
FORMAT {///28X, 'ZEHO SEGMENT
1» COMPUTATION DISCONTINUED')
200	FORMAT!3X,6(12,,-',l2»3X»F7.0»6X))
201	F0RMATI6II3,FIO.O>)
300 FORMA T I	«1
1 CHARAC. LENGTHS OF SEGMENTS (FT)
1FACE LENGTH INTERFACE LENGTH
2E LENGTH INTERFACE LENGTH")
1110	FORMAT! 3 (3F8.0,12))
1111	F0RMAT13I3F6.0,12))
DATA
OATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
OATA
DAT A
DATA
DATA
///• INTERFACE LENGTH INTERDATA
INTERFACE LENGTH INTERFACDATA
DATA
DATA
DATA
NUMBER IN INTERFACE 12,
12,
1112 FORMAT(1
2000 FORMAT!1
1« NUMBER
1
RAW INPUT DATA FOR INTERFACE PARAMETERS'///)DATA
l'//40X,10A4////10X, '~~SYSTEM PARAMETERS*** ,//10X,
OF SECTIONS = ',14,25X » * IPRNT = ',14/10X,•JCON = ',14,
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
DATA
~	DATA
~DATA
~	DATA
DATA
DATA
DATA
DATA
DATA
~	OATA
C	INITIALIZE MATRICES TO ZERO	*DATA
C	*DAT A
£4 4 444 444444 44444 4444444444444444444444444444444444444444444444444 444*4 DATA
DO 1001 1=1,N	OATA
DO 1020 J=1,N	DATA
1020 LA(I,J)s0•0	DATA
239X,•IPNCH = •,14)
2001	F0RMAT(//10X,•~~BOD RATIOS^*'// lOX.'FLC = »tF6.3,
138X, ' F.LN = * ,F6.3)
2002	FORMAT(//10X, '**SCALE FACTORS***//10X,
1'SCALE( 1 ) = ',FL0.3,4X,'SCALE 12) = '.F10.3.4X,
2'SCALE(3) = ',F10.3.4X,'SCALE!4) = '.F10.3)
2200 FORMAT! • 1' , 3!4X, 'INTERFACE *,4X, 'AREA' ,7X,'E',6X,'Q,,4X)/3(
I5X,'ROW-SEG 
-------

DO 1021 J=1»18
DATA
058
1021
ARRAY I I,J ) =0. 0
DATA
059

00 1001 J= 1,6
DATA
060
1001
I ARAY( I ,J ) =0
DATA
06 L
C**********************************************************************
OATA
062
c

*DAT A
063
c
READ AND WRITE INTERFACE PARAMETERS AREA, E AND U AND PLACE
*0 AT A
064
c
THEM INTO ARRAY
~ DATA
065
c

*DAT A
066
$ *$ * *« *£ $ £$ * $ «*$ *** 4* *4 « $ $4 « **£ 4 $ $ 44 $ $ 4* $#$* $$$
DATA
067
1031
DO 1002 1=1,N
DATA
068

READ(MX,1111) (AREA(J),E(J),Q(J)~ IARAY(I,J),J=1,6)
DATA
069
1007
DO 1002 J= 1, 6
DATA
070

JJ=(J-1)*3+l
DATA
071

ARRAY( I,JJ)=AREA(J)*SCALE< 1)
DATA
072

ARRAY(l,JJ+l)=E(J)*SCALE(2)
DATA
073
1002
ARRAY!1,JJ+2)=Q(J>*SCALE<3)
DATA
074

DO 1008 1=1,N
DATA
075

DO 1008 JJ -1, 6
DATA
076

IF(IARAY( I,JJ) ) 100 8,1008,1088
DATA
077
1088
IF(IARAYt I,JJ )-I) 1012,1008, 1012
DATA
078
1012
JJJ=(JJ-1)*3+1
DATA
079

DO 1004 J= 1,6
DATA
080

JK=(J-1 1*3*1
DATA
081

NK=IARAY(I,JJ )
DATA
082

IF(IAR AY(NK, J)>1005,1006, 1005
DATA
083
1006
ARRAY( NK , JK)=ARRAY( I,JJJ)
DATA
084

ARRAY(NK , JK+1)=ARRAY,ARRAY(1,12),
DATA
097

1 I , IARAY(I,5) ,ARRAY(I,13).ARRAY( I ,14),ARRAY! 1,15),
DATA
098

1 I.IARAY(1,6).ARRAY(I,16).ARRAY(1,17),ARRAY(I,1e)
DATA
099
1032
CONTINUE
DATA
100
£*$#**#***4************************************************************
DATA
101
C

~ DATA
102
c
READ AND WRITE INTERFACE PARAMETER- CHARACTERISTIC LENGTH
~ DATA
103
c

~ DATA
104
£$4*6*$**4<********************************** **********~**********#*****
DATA
105

WRITE* NX, 300)
DATA
106

DO 500 1=1,N
DATA
107

READ (MX,2 0 1 ) J I» L A( I , J I ),JJ,LA ( I ,JJ),JK,LA(I ,JK) ,
DATA
106

1 JL,LA(I,JL),JM,LA(I,JM),JN,LA(|, JN)
DATA
109

LA(I,J I )=LA(1,J I)*SCALE(4)
DATA
110

LA(I,J J)=LA(I,JJ)*SCALE(4)
DATA
111

LAU,JK)=LA(I,JK)*SCALE(4)
DATA
112

LA(I,JL )=LA( I ,JL)*SCALE(4)
DATA
113

LA(I,JM ) =LA(I,JM)*SCALE(4)
DATA
114

LA(I,JN)=LA(I,JN)*SCALE(4)
DATA
115

-------
WRIT E(NXi200 > I,J I,LA(I,J I),I,JJ,LA(I,JJ),I,JK,LA(I,JK),	DATA	116
1 I, JL,LA{I,JL>,I,JH,LA(l,JM),I,JN,LA(I,JN)	DATA	117
500 CONTINUE	OATA	118
£********************************************************************** DATA	119
C	«DATA	120
C READ SECTION PARAMETERS- OEPTH, VOLUME AND TEMPERATURE	*DATA	121
C	«DATA	122
C********************************************************************** DATA	123
READ(MX,31 MH(l) ,1 = 1 ,N]	DATA	124
READ(MX,31) (TCI),VOL<1),1=1,N)	DATA	125
IF(T(2)>14,13,14	DATA	126
13	DO 14 1=1,N	DATA	127
T(I)=T(1)	DATA	128
14	CONTINUE	DATA	129
£********************************************************************** OATA	130
C	*OATA	131
C APPLY CONVERSION FACTORS	DATA	132
C	~DAT A	133
£*******************************************************,*************** DATA	134
DO 12 1=1,N	DATA	135
DO 23 J = 1,6	DATA	136
JJ=(J-1) *3+1	DATA	137
23 ARRAYU , JJ+2)=ARRAY( I, JJ+2)*.6463	DATA	138
12 VOL < I)=VOL( I)*7.48	DATA	139
RETURN	DATA	140
999 STOP	DATA	141
ENO	DATA	142

-------
SUBROUTINE REVI	REVI COO
£######*##*$#*##*######***********#*##*#**«$*#######*##*#######*#**#*** REVI OOl
£**#*#*#****«*«***#«*##**##*#######*********#******##*#***##*##'*####**# REVI 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
REV I 009
C#*##**#*##########*#**####*#######*###########*#*####*###****###*##### REV I 010
DIMENSION Q(6),E(6),REVIS(12)	REVI Oil
COMMON/ALGAE/NX,N,I NO IC,ICON,MMM,NMK	REVI 012
COMMON/COHOE/L A,SCALE(4).ARRAY(100»18)*IARAY(100.6)	REVI 013
2 FORMAT(1H1,40X'REVISED PARAMETER LIST'//)	REVI 014
5	FORMAT(12F5.0)	REVI 015
10	FORMAT<<#$,G>##*###*##*#*#######*#*##################################$ft£VI 048
9	DO 4 1=1,N	REVI 049
1	F( l-( 1/2)^2)23,23,24	REVI 050
24 READ(MX,5)REVIS REVI 051
IP0S=0	REVI 052
23 DO 4 J=1,6	REVI 053
1P0S=IP0S+1	REVI 054
JJ=(J-1)#3+1	REVI 055
4 ARRAY( I,JJ+2)=REVIS( IP0S)*SCALE(3)»0.646 3	REVI 056
WRITE (NX,14)	REVI 057

-------
00 15 1=1, N	REVI 058
DO 16 J= 1,6	REVI 059
JJ=(J-l)*3*l	REVI 060
16 Q(J)=ARRAY(l,JJ+2)/(SCALE(3)*0.6463)	REVI 061
15 WRITE (NX,13)( I, IARAY(I,J>, Q(J)f J = l, 6)	REVI 062
99 RETURN	REVI 063
END	REVI 064

-------
q*****
<:*~~«*
c
c
c
c
c
c
c*****
214
302
216
37
29
27
213
SUBROUTINE FLOW
*****************************************************************
**************************************************** *************
SUBROUTINE FLOW
THIS SUBROUTINE CALCULATES A FLOW BALANCE AROUND EACH SECTION
TO INSURE THAT THE FLOW REGIME HAS BEEN INPUT CORRECTLY
* ****************************************************************
*****************************************************************
REAL LA<100,100)
COMMON/ALGAE/NX,N,INDIC,ICON * MMM,NMK
COMMQN/COHOE/LA,SCALE(4),ARRAY(100#18),IARAYl100,6)
NMN=0
FORMAT( '	SECTION ',13,' HAS AN EXCESS FLOW OF »,F12.5,
L' CFS' )
FORMAT('1•)
DO 213 I = I, N
0Q=ARRAY(1,3)+ARRAY!I» 6)*ARRAY(I» 9)+ARRAYII,12) +ARRAY!I,15)
L +ARRAY(1,18)
QQ=QQ/.6463
IF1QQ-.0011216,37,37
IF(QQ+.001)37,37,213
IF(NMN)29,29,27
NMN= 1
WRITE(NX,302)
WRITE(NX,214)I,QQ
CONTINUE
RETURN
END
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
~	FLOW
~	FLOW
~	FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
FLOW
000
001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016.
017
018
019
020
021
022
023
024
025
026
027
028
029
030

-------
SUBROUTINE SETUP
£********>!! *************************************************************
£***«**********««***********#***********«**********************«*******
C
C	SUBROUTINE SETUP
C
C	THE PRIMARY TASK OF SETUP IS TO CONSTRUCT THE MATRIX, AB
C
£******:«:** *««*«** *********************************** ******** ***********
£*** ****** **** ***************** * ******4****** * ******** * ****************
REAL LAI 100 ,100)
COMMON/ALGAE/NX,N,INOIC,I CON »MMM,NMK
CQMMON/COHOE/LA,SCALE!4) ,ARRAY(100,18] ,IARAY(100,6 I
COMMON/PIKE/JCON,H!100),T(100)~VOL!100)~AB(100,100)
C0MM0N/8REAM/EPR[M(100),ALPHA(100)
23	FORMAT{ lX,6(I2t,-,,I2.F7.0,F6.3,3X>)
24
SETUPOOO
SETUPOOl
SETUP002
SETUP003
SETUP004
SETUP005
SETUP006
SETUP007
SETUP008
SETUP009
SETUP010
SETUP011
SETUP012
SETUP013
SETUP014
SETUP015
ALSETUP016
SETUP017
SETUP018
SETUP019
SETUP020
SETUP021
SETUP022
SETUP023
SE TUP024
FORMAT! ' 1 • ,49X, ¦ VALUES OF ALPHA AND EPRI M* ,///,6 ( • INTRFC EPRIM
1PHA • ) ,/,6!8X,*(MGOJ»,8X),//)
42	FORMAT 1 4X,I7,3F10.2)
43	FORMAT(* IS3X,'SECTION TEMPERATURE VOLUME	DEPTH*/
115X.MC1 <10#*6GAL) {FT}*//)
WRITE (NX,24)
DO 2 8 I = 1»N
DO 99 J=1»N
99 AB(ItJ >=0.0
C********************************************************************** SETUP02 5
C	*SETUP026
C	CALCULATE EPRIM AND ALPHA FOR EACH INTERFACE	*SETUP027
C	* SE TUP02 8
q****************************************************************** **** SETUP029
DO 20 J = 1» 6 SETUP030
EPRIM!J)=0.0 SETUP031
ALPHA(J)=0.0 SETUP032
JJ=I J—1)*3+1 SETUP033
K=IARAY( I, J ) SETUP034
IF t K ) 20» 20, 10 SE TLfP035
10 EPRIM(J ) =ARRAY(I,J J)*ARRAY(I,JJ+1)*417.1166/(LA IK,I)+LA18,19,19	SETUP050
18	A8( I, I ) = Ati ( I, I)+£PRIM(J)+(1.0-ALPHA(J) )*ARRAY( I, JJ + 2)	SETUP05 1
IF!I-KJ21.20,21 SETUP052
21	AB! l,K)=ALPHA(J)*ARRAY!I,J J + 2) - EPR IM!J)	SETUP053
GO TO 20 SETUP054
19	AB! I,I ) = AB( I,IJ+ALPHA!J)«ARRAY(I,JJ + 2)+ EPRIM(J)	SETUP055
IF(I-K) 22,20,22 SETUP056
22	AB
-------
20 CONTINUE	SETUP058
£t*********** ft******************** **************************** ********** SETUP059
c	*SETUP060
c	WRITE VALUES OF EPRIM AND ALPHA	*SETUP061
c	*SE TUP062
£*********«***********¦**** ************** ****** ******** ******************$£ TUP 06 3
WRITE(NX,23 )	SETUP064
II,IARAYII,1),EPRIMI1),ALPHA(1),1,IARAY,EPRtMi 4) ,ALPHA14), SETUP066
21,1ARAY(1,5),FPRIM{5), ALPHA i 5),1,1ARAY(I,6),EPRIM(6),ALPHA(6)	SETUP067
28 CONTINUE	SETUP068
£********************************************************************** SETUP069
c	~SETUP070
c	WRITE SECTION PARAMETERS	*SETUP07l
c	* SETUP072
£$***********#********************«*******~**************************«* SETUP073
WRITE(NX,43)	SETUPO 74
WRITE!NX,42)(!,T(I),VOL(II,H(I),	1=1,N)	SETUP075
DO 32 1=1,N	SETUPC76
32 HCI)=3.281/H(I»	SETUP077
RETURN	SETUP078
END	SETUP079

-------
subroutine 6CC
BCC
000

acc
001
£**#$**#$**********4#****«*«****#***#*•#*************«**#«**
BCC
002
c
BCC
003
C SUBROUTINE 8CC
BCC
004
C
BCC
005
C THIS SUBROUTINE INCORPORATES BOUNDARY CONDITIONS IWO THE
BCC
006
C SYSTEM
BCC
007
C
BCC
008

BCC
009

BCC
010
REAL LA( lOOt100 )
BCC
ou
DIMENSION ICOLt100),BC<100)
BCC
012
COMMON/ALGAE/NX,N,INDIC,I CON,MKM,NMK
BCC
013
COMMON/BASS/MX,IPRNT,FRCBC(100)
BCC
014
COMMON/COHOE/LAjSCALE(4}tARRAYt100,left1ARAYIlQQt6)
BCC
015
COMMON/BREAM/EPRIM(100),ALPHA(100)
BCC
016
7 FORMAT < 12)
BCC
01.7
44 FORMAT* * 1 ' , 40X, • ZERO BOUNDARY CONCENTRATIONS*}
BCC
018
209 FORMAT t'1 SEGMENT' i 5X,1CL BOUNDARY CONDITION(MG/L> •
BCC
019
L /(I4,I9X,FLG~2])
acc
C20
210 FOftMATi • ISEGHENT • , 5* , ' B^D EQlNO.VRY CDNOI T 1 DM MG/L 1 •
acc
C2I
1 t {14,19X,F10.2))
BCC
022
211 FORMAT* •1SEGMENT',5X, 'DEF BOUNDARY CONOITI ON[MG/L> '
BCC
023
I /U4,19X,F10.2M
BCC
024
1150 FORMAT 16 I FB.0,14))
BCC
025
£*****»******************«****************************************#44**
BCC
026
c
~ BCC
027
C READ AND WRITE BOUNDARY CONDITIONS
~ BCC
028
C
~ BCC
029
c# ~ ~ ~ + 4 * * * 4 * 4= * 4 £ <¦ *r * 4 4:t4** 444,444444c «****# 4443*44444ft
BCC
049
213 KLAS T=I
BCC
050
00 60 r=l,N
BCC
051
FRCBCI N=0.0
BCC
052
IF i ICOL(KLAST)—I )60 > 2 5,60
BCC
053
25 00 30 J=I,6
BCC
054
JJ=
-------
30
CONTINUE
ecc
058
35
ALPHA* J)=0 . 5
8CC
059

EPRIM< J >=ARRAY(1,JJ)~ARRAY(I~J J*1)*417.1166/(LA(M,I)+LA(I»M) )
BCC
060

IF(-0.5+EPRIM(J)/ABS(ARRAY(I,JJ+2)))40,45,45
BCC
061
40
ALPHA(J) = 1.0-EPRIM(J)/(2.0*ABSIARRAY(I,JJ + 2)) )
BCC
062
45
IF(ARRAY( I,JJ+2 >) 50,5 5,55
BCC
063
50
FRCBC(I>=8.34*BC(KLAST»*IEPRIM
BCC
067

KLAST=KLAST+1
BCC
068
60
CONTINUE
BCC
069
1000
RETURN
8CC
070

END
BCC
071

-------
CONSTRUCT THE
COMPLETE THE
SUBROUTINE RATE
£*£************************ ************** ******4* ******************** *9
£***~*****~********************»*****************»*********»#**********
c
c	SUBROUTINE RATE
C
C	THE PRIMARY PURPOSES OF THIS SUBROUTINE ARE TO
C	SPECIFIC CONSTITUENT SYSTEM MATRIX, AC AND TO
C	APPROPRIATE FORCING FUNCTIONS
C
q**********************************************************************
c**********************************************************************
REAL KA< 100) ,JO)< 100) ,KSTOR( 100) , LOAD
01 MENS ION Y{£&),BD( PMR( ssr)
COMMON/ALGAE/NX, N,INDIC,I CON,MMM,NMK
COMMON/BASS/MX,IPRNT,FRCBC(100)
COMMON/PIKE/JCON,H(100),T(100),VOL(100),AB(100,100)
COMMON/TROUT/W(100),XK1(100),B{100,1),AC(100t100)
FORMAT(16F5.0)
1
12 FORMAT!F10.0,13)
16	FORMAT(• BOO LOAD ='tFl0.0,'
17	FORMAT(• OEF LOAD =',F10.0,'
18	FORMATt•1•)
34	FORMAT( * 1 SECTION CHLORIDE
35	FORMATt2X,15,12X.F14.2)
36	FORMAT('1 THE TEMPERATURE
POUNOS/DAY
POUNDS/DAY
FOR
FOR
SECTION
SECTION
13)
13)
BOUNDARY LOAD (MGD*MG/L)•///)
CORRECTION FACTOR FOR THE BOD REMOVAL
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
IRATE = ' ,F7.3,///,• SECTION BOD REM RATE BOD WASTE LOAD BOD RATE
1B0UNDARY LOAD • / 1 2X,•*TEMP CORR*', 5X ,•(MGD*MG/L)'»7X,•(MGD*MG/L)•) RATE
37	FORMAT(2X,I5,F12.3,7X,F12.3»5X,F12.3)	RATE
38	FORMAT(•1 THE TEMPERATURE CORRECTION FACTOR FOR THE REAERAT I ON RARATE
1TE = •,F7•3,/' THE BOD CORRECTION FACTOR (FL) = ',F7.3,/' THE RATE
1TEMPERATURE CORRECTION FACTOR FOR THE BENTHAL DEMAND =
1 THE TEMP CORR FACTOR FOR KD = • ,F7.3,///»- SECTION
1 REAER RATE BOUND COND LOAD TOTAL LOAD	PMR
»,F7.3,/' RATE
DEOX RATE RATE
BD */,12RAT E
IX,2(•~TEMP CORR* •),2(
(MGD*MG/L)•,4X)t'MG/L/D *TEMP CORR««)
RATE
39 FORMATt 2X,16,5X,F6.3,7X,F6.3»6X,F12.3,3X,F12.3»2X,F9.3»3X,F9.3) RATE
Q# *RATE
C	READ AND APPLY TEMPERATURE CORRECTION FACTORS	*RATE
C	~RATE
q***********************************************************************kate
I F(MMM—518,386,386
8 READtMX,l)FACl,(KA( I)»1 = 1,15)
I F(KA(2))202,3,202
3 DO 64 1=2,N
64 KA( I)=KA(1)
GO TO 334
202	LMM=N—15
IF(LMM)334,334,203
203	READ!MX,1)	( KA(I),I = 16,N)
334 IF(MMM-3)2,388,2
388	DO 389 1=1,N
389	KSTOR ( I •) =KA ( I )
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
RATE
000
001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
052
053
054
055

-------
FAC4=FAC 1	RATE 056
GO TO 2	RATE 057
386	DO 387 1=1,N	RATE 058
387	KA( I )=KSTORlI)	RATE 059
FAC1-FAC4	RATE 060
2 00 6 1=1,N	RATE 061
KA{ I)=KA(I)*FAC 1 *~(T(I)-20.)	RATE 062
Q********************************************************************** RATE 063
C	*R ATE 064
C CALCULATE THE SPECIFIC CONSTITUENT SYSTEM MATRIX, AC	*RATE 065
C	«RATE 066
£********************************************************************** RATE 067
DO 801 J = 1» N	RATE 068
801 AC ( I, J ) - Afj ( I , J )	RATE 069
6 AC( I, I)=AB( I, I )+KA( U*VOL( I J	RATE 070
I F( IN0ICJ10,10,54	RATE 071
Q* ******************************** ***** ****** ****** ** * *** ********** *****RAT £ 072
C	*RATE 073
C CALCULATE THE BOD OR CHLORIDE FORCING FUNCTION	*RATE 074
C	*RATE 075
q***********************************************************************rATE 076
10 DO 9 1=1,N	RATE 077
W(I ) =FRCBC(I)	RATE 078
9 FRCBC(I)=FRCBC(I) /8•34	RATE 079
LLL=0	RATE 080
1100 READ(MX,12)L0AD,ISEC	RATE 081
IF( ISEC) 13,22,13	RATE 082
13	I F(LLL >28, 28, 14	RATE 083
28 WRITE(NX,18)	RATE 084
14	W( ISEC)=W
-------
56	BD< 1 ) = BD( I )*FAC3 **( T( I 1-20.')
GO TO 57
58 DO 53 1=1,N
PMR(I)=0.
53 BD( I )=0.
57	DO 33 1=1,N
KD-20.)
C *~+.*+ + + + *****~~*** + *** ~~***~ ~***»* ************
c
C	READ AND WRITE LOADS
C	CALCULATE THE DEFICIT FORCING FUNCTION
C
Q*4t*$4******************** *********************
XK1 { I 1 = 8. 34*VOL ( I )*< XKK I )*KD( I >*FL -PMR(
33 CONTINUE
LLL = 0
101 READ(MX,12IL0AD,ISEC
I F (ISEC)2l,22,21
21	IF(LLL ) 326,3 26, 325
326 WRI TE(NX,18)
325 XK1( ISEC)=XK1( ISECULOAD
LLL=LLL+1
WRITE(NX, 17)LOAD,ISEC
GO TO 101
22	IF(INDIC)29,30,31
29	DO 303 1=1,N
303	W( I)=Wl I ) /8.34
WRITE(NX,34)
WRIT£(NX,35)(I,W(I),1=1,N)
GO TO 40
30	DO 304 1 = 1,N
304	W(I)=W(I 1/8.34
WRITE(NX,36)FAC1
DO 853 1=1,N
853	W( I)=W(I)-FRC8C
-------
subroutine write
£*****###$*«9$*r«T#*$****«#*«* ««* «c#«*«**«*««*** £**$$***««*«*
£**$ $*$*4 #£** ft 4* «$#$*#***« ft******# *4#* * ft 4t« *4 ****** #**#£* ****** *********
C
C	SUBROUTINE WRITE
C
C THIS SUBROUTINE OUTPUTS MATRIX AC ANO MATRIX AC**-1
C
£******44444444*4444*4444***44**4444*****4*44*4*******4*4*4**4******44*
£******4444&4*44*444*444******44*444*44**44444'<'*4******4444*$****4<#***4
COMMON/ALGAE/NX,N,INDIC,ICON,MMM,NMK
COMMON/TROUT/W(100),XK1<100),8(100»1),AC C100~100)
800
803
1000
1005
1007
1012
1015
2000
FORMAT(•I',50X,
FORMAT(///48X,•
• INITIAL
INVERTED
SYSTEM
SYSTEM
MATRIX')
MATRIX')
FORMAT{
FORMAT(
FORMATC
/4X,9I13/)
«1',51X,•INITIAL
1',48X,'INITIAL
FORMAT(///50Xt•INVERTED
FORMAT(///47X,~INVERTED
FORMAT ( 13, lX,9F.13.5/14Xt96i3.5) )
JMAX=N
BOO MATRIX')
DO DEFICIT MATRIX*)
BOD MATRIX*)
DO DEFICIT MATRIX*)
MR -I TE000
WR ITEOOl
WRITE002
WRITE003
WRIT E004
WRITE005
WRITE006
WRIT E007
WRITE008
WR I TE009
WRITE010
WRITEOU
WRITE012
MR ITE013
WR1TE014
WRITE015
WRITE016
HRITE017
MRITE018
WRITE019
WRITE020

MM=1

WRITE021

IF(NMK)1003,1003,1010
WRITE022
1003
IF( INDIC>1018,
1004,1006
WRITE023
1018
WR I T E(NX,800)

WR ITE024

GO TO 1008

WRITE025
1004
WRITE(NX,1005)

WRITE026

GO TO 1008

WR1TE027
1006
WR I TE(NX*1007)

WR ITE028
1008
WR I TE(NX,10 00)
(J » J=11JMAX}
WRITE029

DO 1001 I¦= 1, JMAX
WRITE030
1001
WR I TE(NX,2000)
1, 1 AC(I»J),J = ltJMAX)
WRITE031

NMK = NMK«-1

WRITE032

RETURN

WRITE033
1010
IF( INDIC)1019,
1011,1014
WR I T £034
1019
WR ITE(NX,803)

WRITE035

GO TO 1016

WRITE036
1011
WR ITE(NX,1012)

WRITE037

GO TO 1016

WR ITE038
1014
WRIT6(NX,1015)

WRITE039
1016
WRITE< NX,1000)
(J,J=1,JMAX)
WRITE040

DO 1002 1 = 1»JMAX
WRITE041
1002
WRITE(NX, 2000)I,{AC(I,J),J=1,JMAX)
WR ITE042

NMK = 0

WRIT £043

RETURN

WR ITE044

END

WR ITE045

-------

SUBROUTINE MAIN

MATN
000
£********£*******************************«*****************************
MATN
001
C***** ****** *4 *************************************************** ****««
MATN
002
C


MATN
00*
c
SUBROUTINE MATN
MATN
004
c


MATN
005
c
THIS SUBROUTINE INVERTS THE MATRIX AC
AND MULTIPLIES
MATN
006
c
IT BY THE FORCING FUNCTION

MATN
007
c


MATN
COS
£****«**********************44**********************************4******
MATN
009
£*»****#***************************************************************
MATN
010

DIMENSION IPIVO(1001,INOEX(I00,2),PIVOT<100)
MATN
Oil

COMMON/ALGAE/NX,N»INOIC,ICON,MMM,NMK

MATN
012

COMMON/TROUT/W(100),XK1(100),B<100,1),
AC(100,100)
MATN
013

DO 900 I = 1» N

MATN
014

B( I , 1) =6(I11)/1000000.

MATN
015

00 900 J=1» N

MATN
016
900
AC{ 11 J ) = AC (I,J J/1000000.

MATN
017

MH=1

MATN
018
10
DETER =1.0

MATN
019
15
00 20 J = 1, N

MATN
020
20
IPIVO 1 J)=0

MATN
021
30
DO 550 1=1,N

MATN
022
C


MATN
023
C
SEARCH FOR PIVOT ELEMENT

MATN
024
C


MATN
025
AO
AMAX=0•0

MATN
026
45
00 105 J=i,N

MATN
027
50
IF (IPIVO (J)-1) 60, 105, 60

MATN
023
60
DO 100 M=1,N

MATN
029
70
IF (IP IVQ (M) — 1) 80, 100, 740

MATN
030
80
IF (ABS(AMAX)-ABS(AC(J,M)))85, 100, 100
MATN
031
85
1R0W=J

MATN
032
90
ICOLU =M

MATN
033
95
AMAX= AC( J » M )

MATN
03 4
100
CONTINUE

MATN
035
105
CONTINUE

MATN
036
110
IPIVO (ICOLU )=IPIVO (ICOLU )+l

MATN
037
C


MATN
038
c
INTERCHANGE ROWS TO PUT PIVOT ELEMENT
ON DIAGONAL
MATN
039
c


MATN
040
130
IF (I ROW-ICOLU ) 140, 260, 140

MATN
041
140
DETER =-OETER

MATN
042
150
00 200 L = I,N

MATN
043
160
SWAP=AC(I ROW t L)

MATN
044
170
AC(I ROW,L) = AC(ICOLUtL)

MATN
045
200
AC(ICOLU,L1 = SWAP

MATN
046
205
IF(MM)260, 260, 210

MATN
047
210
DO 250 L= I» MM

MATN
048
220
SWAP-B(I ROW,L)

MATN
049
230
B(.IROW,L)=B( ICOLU ,L)

MATN
050
250
B(ICOLU i L)= SWAP

MATN
051
260
INDEXt I, I ) = IROW

MATN
052
270
INOEXl I , 2 ) = ICOLU

MATN
053
310
PIVOTt [ 1 =AC 
-------


HATN
058
330
AC(ICOLU,ICOLU)=1.0
MATN
059
340
00 350 L = 1» N
MATN
060
350
AC IICOLU,LI= AC
MATN
061
355
IF(MM)380, 380, 360
MATN
062
360
DO 370 L = 1» MM
MATN
063
370
BUCOLU ,L)=8(ICOLU ,L)/PIV0TIU
MATN
064


MATN
065

REDUCE NON-PIVOT ROWS
MATN
066


MATN
067
380
DO 550 L 1=1,N
MATN
068
390
IF(L1-IC0LU ) 400, 550, 400
MATN
069
400
S = AC ( L 1»ICULU)
MATN
070
420
ACIL1, ICULU)-=0.0
MATN
071
430
DO 450 L=1,N
MATN
072
450
AC(L 1, L ) =AClLhL)—ACl ICOLU,L)*S
MATN
073
455
I FI MM J 550, 550, 460
MATN
074
460
DO 500 L = 1 ,MM
MATN
075
500
81 LI,L )=B(L1»L)-B(ICOLU ,L)*S
MATN
076
550
CONTINUE
MATN
077


MATN
078

INTERCHANGE COLUMNS
MATN
079


MATN
080
600
DO 710 I=1,N
HATN
081
610
L=N+1—1
MATN
082
620
IF IINDEX(L »15-1 NOEX(L,2)) 630, 710, 630
MATN
083
630
JROW=INDEX{L,1)
MATN
084
640
JCOLU =INDEX IL,2)
MATN
085
650
DO 705 M=1TN
MATN
086
660
SVJAP = AC( M, JROW)
MATN
087
670
AC(M,JROW)=AC(M,JCOLU)
MATN
088
700
AC(M,JCOLU) = SWAP
MATN
089
705
CONTINUE
MATN
090
710
CONTINUE
MATN
091
740
CONTINUE
MATN
092
750
CONTINUE
MATN
093

DO 901 1=1,N
MATN
094

DO 901 J=1,N
MATN
095
901
ACtI,J)=AC(I,J3/1000000.
MATN
096

RETURN
MATN
097

END
MATN
098

-------
SUBROUTINE STORE
C*******«***************************************** *********************
£****«*#******$********************************************************
c
C	SUBROUTINE STORE
C
C**********************************************************************
£****************************«**«#**********#*«*************~**********
DIMENSION CL(100),CS (100)>R(100)
COMMON/ALGAE/NX,N,I NOIC»IC ON, MMM, NMK
COMMON/PIKE/JC0N,H(100) ,T(100).VOL(100),AB(100,100)
COMMON/TROUT/W( 100),XK1<100),B(100,1),AC I100,100)
42 FORMAT(•1 SECTION	DEFICIT(MG/L)')
44	FORMAT(•1 SECTION	BOD(MG/L)')
45	FORMAT(01 SECTION	CHLORIOES(MG/L)•)
48 FORMAT!3X,I4,3X,F14.3)
807	FORMAT(10 X,I3,5(4X,F9.2))
808	FORMAT('1	SECTION CHLORIDES
TEMPERATURE
RESULTS FOR COMBINED
TEMPERATURE
(C)
STOREOOO
STOREOOI
ST0REC02
ST0RE003
ST0RE004
ST0REGC5
ST0RE006
STORE007
STORE008
ST0RE009
ST0RE010
ST0RE01I
STOREO12
ST0RE013
ST0RE014
ST0RE015
STOREO16
STOREO17
ST0RE018
(CARBONACEOUSTORE019
ST0RE020
DEFICIT
DEFICIT
(MG/L)
1SATURATION	DO',//)
817 FORMAT{• 1« , •
4S ANO NITROGENOUS ) RUN*//)
828 FORMAT('	SECTION CHLORIDES
ISATURA TION	DO•,/,19X,•(MG/L)
1 (MG/L)	(MG/L)*)
DO 805 1 = 1,N
805 XK1 ( I ) = (3( 1,1)
C**********************************************************************
c
C WRITES FINAL CONCENTRATIONS
C
C**********************************************************************
11 IF(INDIC)850,840,41
850 WR IT E ( NX , 45)'
GO TO 643
840 WRITE(NX,44)
GO TO 643
41 WRITE(NX,42)
643 WRITE(NX,48)(I,XK1(I)~I-1,N)
9999 GO TOI1800,1801,1802,1801,1803),MMM
1800	GO T0U199, 1199, 813,813), JCON
C**********************************************************************
C
C	CALCULATES THE SATURATION VALUE OF DO
C
Q**********************************************************************
813 DO 815 1=1,N
CL(I)=XK1(I)
815 CS( I ) = ( l.-.000009*CL(I))*(14.652-.41022*T(I)+.0079910*
IT(I )**2.-.000077774*T( I )**3. )
MMM=MMM+1
RETURN
1199 DO 1139 1 = 1,N
CL{I)=0.
1139 CS(I)=(l«-.000009*CL
-------

RETURN
ST0RE058
1802
WRITE(NX,808)
ST0RE059

DO 806 1 = 1,N
ST0RE060

Wl I )=CS(I) - XK1{I)
ST0RE06I

URITE(NX,607) I,CL(I) tT(I) tXKl(I) tCS( I) iWU)
ST0RE062
806
CONTINUE
STOR6063

GO TO<93,1804,93,1804),JCON
STORE064
93
RETURN
ST0RE065
1804
DO 829 1=1,N
ST0RE066
829
R( I ) =XKI( I )
ST0RE067

1
ST0RE068

RETURN
STORE069
1803
WRIT E(NX,817)
ST0RE070

WRITE(NX,828)
ST0RE071

00 881 1=1,N
ST0RE072

XK 1 ( I)=XK1(I )+R< I )
ST0RE073

W(I)=CS( I) — XK I ( I )
ST0RE074

writeinx, 807) j ,clu) ,tc J >» xKim ,csm »wm
STORE075
881
CONTINUE
ST0RE076

RETURN
ST0RE077

END
ST0RE078

-------
USER'S MANUAL
U-l

-------
The manner in which, the data deck for HAR03 is structured
is shown in figure U-J. 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 modeled 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) JCON: 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
some 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
each is computed separately and their deficits added
to determine the total deficit.
5)	Estuarine coypled reactive substances (Chlorides,
BOD-deficit) HAR03 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 INDIC and JCON.
A description of .input and output is then followed by an example
problem including sample input and output to further clarify the procedure.
Finally, some restrictions on the use of the program are
given.
U-3

-------
^ START ^

/ JCL

CARDS

	V —

/ TITLE

CARD.


/ INTERFACE

PARAMETER

CARDS


/ LENGTH

CARDS


/ DEPTH

CARDS


( TEMPERATURE

VOLUME

CARDS


/ ICON

( END OF FILE
21-> C/*)
CARD
STOP ^

0


( NUMBER OF
BOUNDARY
CND.
CARD
B.C.
CONCENTRAT.
CARDS
l(4efIcit)
s 0
(chlorides,
BOD)
s


RATE
CARDS

to



/
PROTOSYNTE.

RATE

CARDS
I
BENTHAL
RATE
CARDS
LOAD
CARDS
Figure U-l: Structure of data deck for HAR03
U-U

-------
Combinations of reactants
Chlorides	BOD
(Conservative (Single
substances)	reactive
substance)
BOD
DEFICIT
(coupled
reactive
subs tances)
CBOD
NBOD
DEFICIT
(addl tl ve
coupled
subs tances)
Chlorides
BOD
DEFICIT
(es tuar I ne
coupled)
ChlorIdes
CBOD
NBOD
DEFICIT
(es tuar I ne
addltlve)
J CON
CARD SETUP
JCL
GS CARDS
^INDIC-IEXIT
CARD
(INDIC°-l)
SC CARDS
(ChlorIdes)
/INDIC
-1 EX IT
CARD
(IEXIB>1)
¦<
r
/

/*
JCL
GS CARDS
flNDIC-IEXI T
CARO
(INDIC-0)
INDIC-IEXIT
CARD
(IEXIT-1)
JCL
JCL
JCL
GS CARDS
GS CARDS
INDIC-IEXIT
GARD
(INDIC-0)
GS CARDS
INDIC-IEXIT
CARD
(INDIC--1)
.SC CARD S
(ChlorIdes)
INDIC-IEXI T
CARD
(INDIC-0)
S INDIC-IEXIT
~ CARD
(INDIC-0)
+¦

4-




/
/

/

SC CARDS

SC CARD <

SC CARDS

SC CARD S
(BOD)

(BOD)

(CBOD)

(BOD)
+
+ ...


flNDIC-IEXIT

( INDIC-IEXIT

^INDIC-IEXIT

( INDIC-IEXIT
CARD

CARD

CARD

CARD
(INDIC-1)

(INDIC-1)

(INDIC-1)

(INDIC-1)


¦4r

/
/ ¦
/
/
SC CARDS

SC CARD 5

SC CARD $

SC CARDS
(deficit)

(deficit)

(deflclt)
_L

(deficit)
INBIC-IEXIT
CARD
(INDIC-0)
SC CARD J
(NBOD)
INDIC-IEXIT
CARD
(INDIC-1)
:s:
INDIC-IEXIT
CARD
(IEXIT-1)
SC CARDS
(deficit)

'
(^INDIC-IEXIT
CARD
(IEXIT-1)

r
/
/
•
INDIC-IEXIT
CARD
(IEXIT-1)
/*
JCL
/


GS CARDS

/
INDIC-IEXIT

CARD

(INDIC--1)

/


SC CARDS

(Chlor1des)

/
INDIC-IEXIT

CARD

(INDIC-0)

/


SC CARD $

(CBOD)

y
INDIC-IEXIT

CARD

(INDIC-1)

/


SC CARDS

(deficit)


INDIC-IEXIT

CARD

(INDIC-0)
. ±
/


SC CARDS

(NBOD)
( INDIC-IEXIT
CARD
(INDIC-l)
+
SC CARD f
(defIc11)
INDIC-IEXIT
CARD
(I EX I T-l)
/*
Figure U-2
Deck setups for various combinations
of reactants using program HAR03.
o-5

-------
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)
Pascription	Units	Format
To be used to label the
run	A40
The number of sections in
the model	13
An indicator,	13
IPRNT=1 for printout of
matrices
IPRNT=0 to suppress the
printout of
matrices
An indicator used to desig-	13
nate the type of system being
modeled
J C0N=1 BOD
Deficit
JCON"2 CBOD
NBOD
Deficit
JC0N=3 Chlorides
BOD
Deficit
JC0N=4 Chlorides
CBOD
NBOD
Deficit
Area scale factor	ft^/ou* F7
Dispersion scale factor mi^/day/ou F7
Flow scale factor	cfs/ou F7
Length scale factor	ft/ou	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-20)
IARAY
Section forming the first




interface with section

12
(21-26)
AREA
Cross sectional area of




second interface
ft2
F6
(27-32)
E
Dispersion coefficient of
mi2/day



second interface
F6
(33-38)
Q
Flow across the second




interface
cf s
F6
(39-40)
IARAY
Section forming the second




interface

12
(41-46)
AREA
Cross sectional area of third




interface
ft2
F6
(47-52)
E
Dispersion coefficient of




third interface
ml2/day
F6
(53-58)
Q
Flow across the third inter-




face
cf s
F6
(59-60)
IARAY
Section forming the third




interface

12
(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
Description
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,JL)
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.



j

JN
JK
JL
J
' 1
*
VI
JJ
! ^
MD I
JM
JI

figure U-l: Depiction of characteristic lengths
U-8

-------
where
Q - LA(I,JK) and LA(I,JI)
Q) - LA(I,JJ), LA(I.JH) IA(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 lt>6ft3
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
(1-2)
Variable
ICON
Description	Units
An indicator. If desired
at this point changes in
the, E's or Q's which were
read in previously, may be
made using ICON
Format
12
If ICON=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-U (new
Q Cards'')
[NEW E CARDS (OPTIONAL-MUST BE PRECEDED BY ICON CARD=1)]
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 E's
II
II
II
11
ir
it
ii
ir
it
ii
ii
Description	Units
For 1st Interface with first section Mi
2nd
-*2/day
Format
F5
"3rd
4	th
5	th
6	th
1st
2nd
3rd
4	th
5	th
6th
IT
11
If
II
II
second
I*
M
11
H
II
IN SIMILAR FASHION, REVISIONS ARE INPUTTED FOR THE THIRD TO Nth SECTIONS
note: These cards must be followed by another "ICON CARD" (as explained
o-Wove)* If revisions are to be made in flow, ICON must equal 2,
and the new values can be read as on page U-\\. 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-ll .
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
4	th
5	th
6	th
1st
2nd
3rd
4	th
5	th
6	th
?f
ff
Second
ff
II
II
IN SIMILAR FASHION, REVISIONS ARE INPUTTED FOR THE THIRD TO Nth SECTIONS
These cards must be followed by another
"ICON CARD" (as explained on page U-IO ).
If revisions are to be made in dispersion,
ICON must equal 1 and the new values can
be read as described on page U-IO. 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
(3-4)
IEXIT
Description
An indicator used to denote whether
subsequent DATA will be for Chloride,
BOD or D.O. deficit
INDIC = -1 for Chloride
INDIC = 0 for BOD
INDIC = 1 for D.O. deficit
An indicator which can either continue
or terminate the run
Units
Format
12
12
If IEXIT>0
IEXIT<0
The run will Terminate
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)
(9-12)
BC (1)
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)
FACl
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 NBOp	=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 t.his 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
U-12

-------
[RATE CARD #2] (Optional-for D.O. deficit only)
Columns Variable
(1-5)	FAC2
(6-10)
(11-15)
(16-20)
FL
KD(1)*
KD(2)
Description	Units Format
Temperature correction factor for	F5
second reaction rate (deoxygenation
rate)
Ratio of ultimate to 5-day BOD	F5
Deoxygenation rate for section 1**	1/day F5
Deoxygenation rate for section 2	1/day F5
(76-80) KD(14)
Deoxygenation rate for section 14
(proceed until N rates have been input)
[PHOTOSYNTHESIS RATE] (Optional for D.O. deficit only)***
1/day
(1-5)
(6-10)
PMR(l)
PMR(2)
Net photosynthesis rate for section 1 mg/l/day
Net photosynthesis rate for section 2 mg/l/day
F5
F5
F5
(76-80) PMR(16)	Net photosynthesis rate for section i£ mg/l/day
(proceed until N rates have been input)
[BENTHAL DEOXYGENATION RATE] (Optional for D.O. deficit only)***
(1-5)
(6-10)
(11-15)
FAC3
BD(1)
BD(2)
Temperature correction factor for
benthal demand
Benthal demand for section 1**
Benthal demand for section 2
F5
gm/M^/day F5
gm/M^/day F5
(76-80) SD(15)
Benthal demand for section 15
(proceed until N 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 KD1s 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 NBOD and then adds them it is unnecessary to input the
reaeration rates, the photosynthesis rate and the benthal demand with the
NBOD 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

-------
[LOAD CARDS]
Columns Variable	Description	Units	Format
(1-10) LOAD	LOAD to a section	it/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 arp 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 U-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 water are transferred
between the flats and the deeper section.
U-15

-------
figure U-A: 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 usingIHAR03 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 shown 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 qards 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-5: Dissolved oxygen concentrations (ing/1)
Figure U-6. Segmentation scheme for th.e bay
U—17

-------
1-1 1-2
Area (ft2)	39600
Dispersion (mile2/day)	0.5
Flow (cfs)	150. 150.
Boundary condition (mg/1)
Chlorides
CBOD	2.
NBOD	2.
DO deficit	1.
cj
i
M
00
Temperature (°c)
Depth (feet)
Volume (10^ feet^)
CBOD removal rate (1/day)
NBOD removal rate (1/day)
reaeration rate (1/day)
CBOD deoxygenaticr. rate (1/day)
NBOD deoxygenaticr.. rate (1/day)
CBOD load (#/day)
NBOD load (#/day)
Bent'nal lead (gm/M2/day)
Algal load (mg/l/day)
Waste load flow (cfs)
Interfaces
2-3
3-4
4-5
4-6
6-7
5-7
5-8
8-8
95040, 100320. 105600. 36960. 21200. 31680. 99000. 158400.
.5	1,0	1.0	0.5	0,5 0.5 1.5 1.5
150. 150, 243.	243. 243.
1000.
.5
Table
U-l:
Interface
data for
run with
load




to section
4





Sections





1
2
3
4
5
6
7
8
21.0
22.0
22.0
22.0
22.0
4.0
24.0
22.0
12.0
13.0
15.0
20.0
22.0
4.0
5.0
30.0
83.64
271.8
418.18
557.57
613.32
111.51
139.30
836.35
.3
.3
.3
.3
.3
.3
.3
.3
.3
.3
.3
.3
.3
.2
.3
.3
.24
.23
.22
.17
.15
.16
.18
.12
.1
. i
.1
.1
.1
= 1
, 1
.1

1
: i

:1
= 1
. 1
.1
100000.
100000.
93.
,05
25
Table U-2: Section data with load to section 4

-------
Figure U-7a: Results of model rurv with, load to
segment 6.
Figure U-7b: Results of model run with load to
segment 4.

-------

-------
EXAMPLE.RUN (LOAD TO SECTION *)
•~SYSTEM PARAMETERS**	
NUMBER _OF SECTIONS 5	«	 		 _ 	 IPRNT -	 0
JCON 3 4	IPNCH - 0
•~SCALE FACTORS**	.	
SCALE! 1) - 1000.000 SCALE(2) -	1.000 SCAIEO) *	1.000 SCALEU] « 5280.000

-------















i INTERFACE
AREA
E
Q
INTERFACE
AREA
c
Q
INTERFACE
AREA
F
a
ROM-:
SEC
(FT»»2)

0.
0.0
O
•
o
1
8- 0J
0.
0.0
¦o.-o
I
8- 0J
0.
o
•
o
0.0

-------
CHARAC. LENGTHS OF SEGMENTS (FT)
INTERFACE LENGTH INTERFACE	LENGTH INTERFACE	LENGTH	INTERFACE LENGTH INTERFACE LENGTH INTERFACE LENGTH
1-	1 5280. _ 1- 2	5280.	1-0	0. 	 1-0	0.	1- 0	0. 	 1- 0	0.
2-	I 5260. 2- 3	5280. *2-0	0.	2-0 0. 2- 0 0."	2- 0	0.
3-	2 5280.	3— »	5280.	3- 0	0.	3-0	0.	3- 0	0.	3- 0	
4-	3 5280. 4- 5	5280. «- 6	5280.	4- 0 0. 0 0.	0	0.
5-	4 _ 5280. 5- 8	5280.	 5- 7	5280.	5- 0	0.	5-_ 0	0.	 5- 0	 _ 0.
6-	* 5280. 	 6- 7.	5280. 6- 0	0.	6- 0 0." 6- 0 0.	6- 0	0.
	7-_A	 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. 6- 0 ~ 0. " 8-0	0.

-------
SECTION 4 HAS AN EXCESS FLO* OF <)2.999<)« CFS

-------
VALUES OF ALPHA AND EPRIH
INTRfC EPRIH ALPHA INTRFC EPRIH ALPHA INTRFC EPRIM ALPHA ~ INTRFC EPRIH ALPHA INTRFC EPRIH ALPHA~" INTRFC EPRIH ALPHA
	 	IM60I 		 (HCO)	(HGO)	_	(HCD) 		 _(HGO) . 	 __.tMGO) .
I-	I 0.	1.000 1-	2	782.	0.500	1- 0	0. 0.0 I- 0 0. 0.0 1-	0 0. 0.0	1- 0 0.	0.0
	1877.	0.500 __ 2-	1	702.	0.500	, 2- 0 0. 0.0 	 2-0 0. 0.0 	 2-	0	0. 0.0	 2- 0 	0.	0.0
i-	4 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
	5 . >171.	0.500	4-	6	730.	0.500	_ 4- 3	3963. 0.500 4- 0	0. 0.0 __ 4-	0	0. 0.0	4- 0	0.	0.0
5-	8 5866.	0.500 5- 7	626.	0.500	5- 4	4171. 0.500 " 5- 0	0. 0.0 " 5-	0 "o.'o.o" ~5- 0 " 0.	0.0
S-	7 419.	0.500 6-	4	730.	0.500	6- 0	0. 0.0	6- 0	0. 0.0	6-	0	0. 0.0	6- 0 0.	0.0
7-	5 626.	0.500 7-	6	419.	0.500	7- 0	0. 0.0 7- 0 0. 0.0 7-	0 0. 0.0	7- 0 6'.	0.0
9- 8	9385^0.500	5866._0.500 	8- 0 _ 0. 0.0 _ 8- 0	0. 0.0	8-_ 0 	0._0.0	8- 0	0.	0.0_

-------
SECTION TEMPERATURE VOLUME DEPTH
IC1 I10»»6GAL)
1	: 21.00 625.63	_ 12.00
2	22.00 ' 2033.06	13.00
J	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	_ 24.00 _ 1042.64	5.00
8	22.00 6255.89	30.00

-------

-------
SECTION CHLORIDE BOUNDARY LOAD 1HGD*MG/L>
	1			0.0 ¦
2	0.0
	0.0
0.0
	.0.0
0.0
	0.0
9306581.00

-------
SECTION CHLORIDES IKG/LI
t	755.9^3
	2	 _ 855.837
3	901.211
	*	923.532 	
5	957.424
	6	 932.199
7	947.311
	8	9fl3.405	

-------

-------
9 BOO LOAD - lOOOOO. POUNOS/OAY FOR SECTION

-------
THE TEMPERATURE CORRECTION FACTOR FOR THE BOD REHOVAL RATE - 1.047'
SECTION BOO REM RATE BOO WASTE LOAD BOO BOUNDARY LOAD
		~TEMP CORR«	(MGD*MG/L ) _ (MG0*MG/L> _
1	0.366	0.0	193.890
	2	o. 38*	o._o	o.jo	
3	0.384	0.0	0.0
	A	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,394		4653.289

-------
SECTION	800( MG/L I
1	0.899
	2	1.035
3	1.938
*	__ 2. 257
5	1.112
	6	 1.335
T 0.846
8		0.638

-------
SEGMENT DEF BOUNDARY CONDITIONING/!.!
1 0.50
_8	L	0.50

-------
THE TEMPERATURE CORRECTION FACTOR FOR THE REAERATION RATE = 1.02*
THE OOO CORRECTION FACTOR (FL) » I.000
THE TEMPERATURE CORRECTION FACTOR FOR THE BENTHAL DEMAND * 1.080
THE TEMP CORR FACTOR FOR KD » 1.0*7
SECTION DEOX RATE REAER RATE BOUND COND LOAD TOTAL LOAD PHR	80	
•TEMP CORR* *TEHP CORR* 
-------
SECTION	OEM C ]TIHG/L )
1 1.643
.. 2	1.B58

3	2-012
4	1.987



5	1.469
6	2.000



7	1.492
8	0.930



















-------
SECTION CHLORIDES TEMPERATURE DEFICIT SATURATION	DO
1
2
3
4
755.94
655.84
901.21
923.53
21.00
22.00
22.00
22.00
957.42
932.20
947.31
963.40
22.00
24.00
24.00
22.00
1.64
1.86
2.01
1.9 9_
1.47
2.00
1.49
0.93
8.78
8.60
8.60
8.59
8.59
8.26
8.26
8.59
7.14
6.74
6.58
JL5_6JL
7.12
6.26
6.77
7.66

-------

-------
BOO LOAD - 100000. POUNOS/OAY FOR SECTION 4

-------
THE TEMPERATURE CORRECTION FACTOR FOR THE BOD REHOVAL RATE » 1.060
SECTION BOD REM RATE BOD MASTE LOAD BOD BOUNDARY LOAD
	 	»TEMP CORR*	(MGD'MG/L) 	tMGO*MG/L I		
1	0.108	0.0	193.890
	2	0.117	0±jQ	Otfi	
1	0.117	0.0	0.0
		0.117	.	11990 • $06	0.0	
i	0.117	0.0	0.0
	6	0.136		 0.0		0.0 _ _
T 0.136 0.0 0.0
	S	0.117	SW!	£*a	

-------
SECTION	BODIHG/L»
1	2.369
	2	 2.635 _
3	3.096
	*	3.615	
5 1.832
	6	2.731
7	1.930
	8	0.681

-------
THE TEMPERATURE CORRECTtON FACTOR FOR THE REAERATION RATE =» 1.024
THE 800 CORRECTION FACTOR t Fl) = 1.000
THE TEMPERATURE CORRECTION FACTOR FOR THE 8ENTHAL DEMAND = 1.080
THE TEMP CORR FACTOR FOR KO = 1.080
SECTION DEOX RATE REAER RATE 80UND CONO LOAD TOTAL LOAO	PMR	B.O
•TEMP CORR* *TEMP CORR* (MCD»MC/L»	tMGO*MG/L) MG/L/O *TEMP CORR*
	1	0. 108	 0.246 	0.0 	 160.041 	0.0 	0.0
2 0.117 0.241 0.0 624.963 0.0 0.0
.,..3 	 0. 1J7 ._ 	.0.231.... 	 0.0 ........ 1130.466 	 0.0	 . _ 0.0
4	0.117	0.178 *"	0.0 " 1758.574	0.0 " 0.0
	5	0.117	0.157	0±0	S80.310	0*3	&..J3.
6	0.136 0.176 0.0 309.920 0.0 0.0
__7	0.136	0.198 			0.0	 273.810 __ 0.0	0.0
6	0.117	0.126	0.0	497.244 ' 0.0	0.0

-------
SEGHENT DEF BOUNDARY CONDITION(HC/L)
1	0.0
	&	;	cue	

-------
SECTION	DEFICIT ( H6/L 1
1	0.965
	2	1.077
3	1.075
_4	0.983
5	0.639
	6	t .085
7 0.902
_8	0.268

-------
RESULTS FOR COHBINED (CARBONACEOUS AND NITROGENOUS ) RUN
SECTION CHLORIDES TEMPERATURE DEFICIT SATURATION	00
	CMG/L)	 (C) __ _(MG/L) _ (MG/L)	< MG/L )
1	755.94	21.00	2.61	8.78	6.17
	2	855.84	22.00	2.94	8^60	i.ftfc	
3	901.21 22.00 3.09 8.60 5.51
	4	 923.53 	 22.00 		_2.9T 	 _ 8.59	5.62
S 957.42 22.00 2.11 8.S9 6.48
	fc	 932.20 	 24.00 	3.08 _ _ 8.26	5.18
7	• 947.31	24.00	2.39	8.26	5.87
	8	983.40	22.00	1.20	8.59	7.39

-------
Restrictions
1)	HAR03 is presently limited to a system of 100 segments
2)	A section may have only one interface act as a boundary
3)	Each section may have up to six interfaces
U-46

-------
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, H., "Computation of Pollution in 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)	HARQ1: A Steady State Estuarine Water Quality Computer Model
prepared by Data Systems Branch U.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 	
fig^e	page
T—1 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 HAR03	U-4
U-2 Deck setups for various combinations of reactants using
program HAR03	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 lpad to.section 4	U-18
U-2 Section data with load to section 4	U-18

-------
NOMENCLATURE
SYMBOL
Documentation
Computer
Program
DEFINITION
UNITS *
A
AREA (J) or
ARRAY (I,JJ),



where
Cross-sectional area
L 2

4,7,10,13,16



J=l,6



1=1, N


[A]
AC (I, J)
Single system matrix. A matrix which
L3/T

1=1,N
represents a single substance which

J=1,N
decays with first order pinetics (see
page T-9)

[A]-l
AC(I,J)
1=1,N
J=1,N
Single system response matrix. The KM/L^)/
inverse of [A], this matrix represents (M/T)]
the response of the system to the input
of a load of a single substance decaying


with first order pinetics (see page
T-10)

[AB] or AB
AB(I, J)
The system matrix. Since matrices fAl


1=1,N
and [B] differ only in a reaction term
l3/t

J=1,N
in their diagonals, a matrix without
such a term has been formulated and
is called the system matrix (see page
T-15)
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 invented.
3 .
L /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 pinetics
(see page T-ll)
M/L3
[B]"1
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.
t (M/L3/
(M/T)]
BD
BD(I)
1=1,N
Benthal oxygen demand
M/L2/T
BOD

Biochemical oxygen demand
M/LJ
c
* units; T=time,
concentration of a constituent
L=length, M=mass
M/L^

-------
	SYMBOL	
Documentation Computer
Program
DEFINITION
UNITS
[C] '"1

The total system response matrix relates
input of the first substance in a coupled
set of reactants to the response of the
second constituent.
[(M/L3/
(M/T)]
CBOD

Component of BOD due primarily to oxygen
utilized by suprophytic organisms as they
utilize the carbonaceous matter in a waste
M/L3
D

Dissolved oxygen deficit ¦ the saturation
value minus the actual concentration of
dissolved oxygen = D0S - DO
M/L3
Do

The initial concentration of dissolved oxygen
deficit at the point introduction of waste
into a stream.
M/L3
DO

Concentration of dissolved oxygen
M/LJ
DOs
cs
Saturation concentration of dissolved
oxygen
M/L3
E
E(J)
ARRAY(I,JJ "M)
where
1=1 $
J=1,6 Longitudinal dispersion coefficient
JJ=1,4,7,10,
13,16
L2/T
E'
EPRIM(J)
J=1 6
Bulk dispersion coefficient
=EA/1
l3/t
H
H
depth
L
k or K

First order reaction coefficient, (also
used to designate a section)
1/T
*a
KA
Reaeration rate (first order)
1/T
*d
KD
Deoxygenation rate (first order)
1/T
Kr
KA
BOD removal rate (first order)
1/T

-------
SYMBOL	 DEFINITION
Documentation Computer
Program
UNITS
1
section length
L
Lo
Initial concentration of BOD at point of
introduction of a waste into a stream.
M/L3
NBOD
Nitrogenous component of BOD primarily
caused be autotrophic bacteria utilizing
nitrogen in waste.
M/L3
PMR
PMR Net algal oxygen rate = photosynthesis
effect — respiration effect of algae on
oxygen
m/l3/t
Q
Q(J) net advective flow
ARRAY(I,JJ+2)
where
J=1,6
1=1 ,N
JJ=1,4,7,10,
13,16
l3/t
t
time
T
T
T(I) temperature
° C
U
advectlve velocity
L/T

-------
SYMBOL
Documentation
Computer
Program
DEFINITION
UNITS
vol
VOL
Volume
L
V
W
Waste load of a constituent
M/T
Wb

Deficit waste load
M/T
wc

BOD waste load
M/T
X

Distance
L
Y
Y
Dissolved oxygen deficit waste load
M/L3
a
ALPHA
weighting coefficient used in the
approximation scheme of the model.

6	1—ALPHA	weighting coefficient = 1-a

-------
Acknowledgements
Thanks must be extended to Richard Wlnfield of the Manhattan College
Department of Environmental Engineering and Science and to Robert Braster,
Salvatore Nolfo and Kevin Bricke of E.P.A. for th^ir thoughtful reviews
of the material contained in this report.
The typing of the report was very patiently and professionally
handled by Ms. Maryann LaBarbera.
S«C * C •

-------