'Jp.iteo c 2"-;;s -•r'rcn 2 Mew Jtirsev "''lew York
Er~".';rc"r"»r ;a; PrfitSCtion Z3 ~-3'1--r3! -'•.•;: 3 ?'., ir'o ; •; j
New /ork NV 10007 Virgin isiancs
FEDBAKS33 - A COMPUTER PROGRAM FOR THE MODELLING OF FIRST ORDER
CONSECUTIVE REACTIONS WITH FEEDBACK UNDER A STEADY
STATE MULTIDIMENSIONAL NATURAL AQUATIC SYSTEM
PROGRAM DOCUMENTATION AND USERS GUIDE
George A. Nossa
Chief, Environmental Systems Section
Data Systems Branch
Management Division
U.S. Environmental Protection Agency
New York, New York 10007
November 1978
-------
-1-
TABLB OF CONTENTS
Abstract ------_______-_______________--- 2
Introduction ------____----------__-------- 2
Theory --------_____----___-_-_-_------- 2
Acknowledgements ---------------------------- 6
References --------_-_---------___-------- 6
Description of the Main Program and Accompanying subroutines ------ 7
Program Flowchart --------------------------- 8
Source Program Listing -------------------------16
Input Description --------------------------- 3Q
Restrictions ------------------------------49
Application to the Delaware Estuary ------------------50
Description of Area of Study -----------------51
Nitrogenous Point Source Loadings --------------53
Program Output ----'--------------------54
Graphs of Computed vs Observed Values ------------ 72
Program Execution at the COMNET facility ----------------83
Acknowledgements ---------------------------6,84
References to the Computer Program -------------------85
APPENDIX A
Listing of Input data for Delaware Estuary Verification - - - 86
-------
-2-.
FEDBAK03 - A COMPUTER PROGRAM FOR THE. .MODELLING OF FIRST ORDER CONSECUTIVE REACTIONS WITH
FEEDBACK UNDER A'STEADY STATE' MULTIDIMENSIONAL NATBRAL AOTJATIC SYSTEM *
George A. Nossa
Environmental Engineer
Data Systems Branch
U.S. Environmental Protection Agency
New York, N.Y.
ABSTRACT
The computer model described 1s used to compute the
steady-state distribution of water quality variables
undergoing consecutive reactions with feedback and
following first order kinetics. The program has
been developed 1n a general form but 1s specifically
applicable to the reactions observed by nitrogenous
species and the associated dissolved oxygen uptake
1n the natural environment.
The basis for this model is the theory of conserva-
tion of mass. The approach used to solve the equa-
tions is.a finite difference scheme developed by
Thomann 14»5/, which has been shown to be a very
effective tool 1n the field of water quality manage-
ment.
INTRODUCTION
A computer model has been developed to serve as a
useful tool 1n the prediction of water quality para-
meters which react under first order kinetics and as
a system of consecutive" reactions .where any parameter
can react 1n .a feedback fashion. The problem setting
assumes an aquatic environment 1n which steady state
conditions can be applied
Lets consider a system of five reactants:
Y
f
KAI
"1
*
"13-
f
14
SCHEMATIC
FIGURE 1
OF FIVE REACTANT SYSTEM
In the above system« the consecutive feedfoward and
feedback reaction rates are presented. All the
possible reaction loops for the first reactant have
also been Included. This particular system will be
used in the theory development.
THEORY
The general estuarine advection/dispersion equation
may be written as:
SN = E62N - U6N - KN + WN(X)
fit "fix? «x
where:
N - concentration of constituent
t = time, 1n tidal cycles
E = tidally averaged dispersion coefficient
which includes the dispersive effects of
tidal motion
U = net advectlve velocity
K = first order decay coefficient of constitu-
ent N
W(x) - direct discharges of N
Assuming steady state conditions: oN/dt
equation (1) becomes:
0 = E d2N - U dN - KN + W(x)
0 and
(2)
Direct solutions for the above equation have been
computed by O'ConnorO .2), a computer program has
been documented by EPA, Region II (3) which uses this
technique to solve for water quality parameters.
A second solution approach developed by
solves the above differential equation directly by
replacing the derivatives with finite-difference
approximations. This approach is used by computer
program HAR03; documented by EPA, Region III6) to
analyze systems of consecutive reactions. The
Thomann solution technique will also be the approach
followed in program FEDBAK03.
If we take the first reactant N-j on figune 1 and we
incorporate all the feedback loops, equation (2)
becomes :
0 = E d2N] . U dNi
--
where Kij are the appropriate first order reaction
rate constants. The incorporation of the reaction
term KH allows for the first order decay of this
first component out of the system.
For the other remaining reactants equation (2) would
be written as:
Q = _ z _ K22N2+K12N2+K32N3+K42N4+
dx
dx
K52N5+W2
<53N5+W3
(5)
* Reproduced from: Environmental Modelling and Simulation, USEPA Office of
Research and Development, Washington, D.C., July 1976 (EPA600/9-76-016)
-------
. U dN4 .
dx~
K54N5+W4
.(6)
. U
<45VW5
•(7)
The Thomann solution approach divides the system
into completely mixed segments, as illustrated for a
one dimensional estuary in figure 2.
FIGURE 2
SEGMENTED ONE DIMENSIONAL ESTUARY
where segments 1 and 'n1 each form an interface with
the boundaries. Equation (3) for the "1" segment on
figure (2) can be written as:
3x2"
dx
(8)
The derivatives on the equation above ean be replaced
by finite-difference approximations giving:
dx
where: E'1fj
(9B)
(9C)-
i i = subscripts used to denote the
interface -between adjacent segments
1 and- j
I 4 a average length of segments 1 and
1>J j
a., j and-*j j are weight factors to
correct the concentrations from
equation 9B
(10)
(ID
In order-for the ^final-concentrations to be positive,
it is required that:
02)
-3-
Substltuting equations (9A) and (9B) Into equation
(8) yields:
,1N2>1+V1K31 f1N3,
(13)
Grouping terms in the above equation yields:
Letting:
»i.i-r-Qi-i,iai-i.rEi-i,i
ai ,i=(5i ,i+l°i .i+T^i-l ,iB1-l
(14)
The general equation for the 1th segment becomes:
,TN4+V1K51 ,1H5. . (18)
The use of this finite difference approximation
scheme has a numerical dispersion, which can be
approximated as (5)
num
.......................... (ISA)
Where £.,^15 the -numerical dispersion. This is par-
ticularly important to stream applications where
advective'-velocities may be high and this effect may
lead to distorted results. -
Extension to multi -dimensional analysis:
If we consider a grid of orthogonally straped
sections such as the one illustrated in figure 3:
k4
FIGURE 3-
HYPOTHETICAtfTWO'-DIMENSIONAL SYSTEM
-------
Foll owing the convention that flows entering a sec-
tion is negative and flow out of a section is posi-
tive, a mass balance due to the transport and disper-
sion of material from section 1 to all surrounding
sections k(s) 1s:(4»5,8)
=°=
dt
(1=1,2. ..n)
t1N5§i
(19)
This equation 1s the equivalent to equation (13) for
the one dimensional case. The generalization of the
advection term is possible since:
and
(20)
(21)
Using equation (19), if the terms containing the de-
pendent variable N-|tf are grouped on the left hand
side and the direct loads of this component and the
terms for the formation of NT due to other components
are placed on the right hand side, one obtains :fe;)
where
(22)
(22A)
For sections where flow enters a section from the
boundary with a concentration C|j
... (23)
and the forcing function at the boundary is added to
the direct loads at that section by:
(24)
For sections forming a boundary, where flow leaves
this section to an area-with a concentration CD:
••' (25)
and
The set of equations for component _Ni for n number of
spacial sections in the system described in figure 1
would be:
allNl ,l+a!2N1 ,
>4+. . .+alnN1 >n=
N5,l
a21Nl ,l+a22Nl ,2+a23Nl ,3+
,n
,2
,2+a33Nl ,
,n=wl ,3+
an!Nl ,l
+a+- • -
,2an3
,n
(29)
(30)
In matrix notation equations (27) thru (30) can be
written as:
[A1](N1)-(W1MVK21](N2)+[VK31](N3)H-[VK41]
(N4KVK51](N5) ........................ (31)
where:
[A]] 1s a square matrix of n order, containing the
a's as defined on equations (22A) and (22B),
note that the main diagonal has the reaction
rate constant K-p
(N-|),(N2),...(N5) are nxl vectors of the reactant
over all sections
(W-]) 1s an nxl vector of the waste loads for react-
ant NT for all sections
[VK2l],[VK31],[VK4l] and [VK^leach of these is an
nxn diagonal matrix of the section volume and
the first order reaction coefficient at that
segment.
A similar analysis as above for the second reactant
N2 on figure 1 yields:
[A2](N2)=(W2)+[VK12](N1)+[VK32](N3)-i-[VK42]
(N4KVK52](N5) ........................ (32)
where:
[A2] is an nxn matrix similar to [AiJ, but the
main diagonal contains the reaction rate
constant K22
(W2) is an nxl vector of direct waste loads for
component N2 over the n sections
For the other reactants N3, N4,
are -generated:
similar equations
d:
[AS] (NaHwaHv*! 3] (NT )+[VK23] (N2)+[VK43]
(N4KVK53](N5) ........................ (33)
[A4](N4)=(W4)+[VK14](N1)H-[VK24](N2)+[VK34]
(N3KVK54](N5) ........................ (34)
(N3)+[VK45](N4) ........................ (35)
The above matrix equations (31) thru (35) can be
written as a matrix of Matrices: (5, 8)
N5,2
-------
-CVK21J-[VK31]-[VK4lJ-[VK5l!
IA3 3-[VK43]-[VK35]
-[VK24J-[VK34] [A4 J-[VK54]
K15] -CVKzsHVKssKVMs] [A5 J
(N2
<»3
(N4
(N5
(N)
(W).
(W2)
(W3)
(W4)
(37)
re [A~] is the 5n x 5n matrix above and (R) and
are 5n x 1 vectors. The solution of the five
ctants over all the spadal sections are given by
(N) = [AT1 (W) ., ......................... (38)
II cation of the theory by the computer program
program described follows a modular approach 1n
:h the user specifies to the main line program the
ions desired, and subroutines are called accord-
y to perform specific tasks.
steps to be accomplished can be summarized as:
nput the physical characteristics of the system;
lamely, the geometry, temperature, hydro! ogic
:haracteristics, reaction schemes and correspond-
ng reaction rates.
alculate E' and °>'s for all the sections as
escribed on equations (9C) and (10). In order to
andle the constraint stated on equation (12), the
rogram tests the expression:
if such is the case
is recalculated as:
C40)
h places °^j well within the tolerable range.
et up the system matrix [Af] by computing its
lements as given on equations (22A) -and {22B);-
t should be noted that the difference between the
A-j] matrices in equations (31) thru -{35 ) is the
ddition of a separate main diagonal term KtKjj'f.
n order to conserve space, this matrix is set up "
ithout this term and during the creation of the
atrix of matrices [A], the appropriate VfKjj^
erm is added.
et up the matrix of matrices [AJ, this 1s done by
ombining an offline disk file containing all the
iKj terms for the system and -the [A^] matrix.-
nput of direct discharges into the system and
oundary concentrations, and from these compute
he system s ource -vector
-------
-6-
If we assume these reactions to follow first order
kinetics, the system can readily be solved using
program FEDBAK03. The computation of dissolved oxy-
gen deficit can be accomplished two ways. A deficit
"species" can be defined (noted Ng above), the decay
of which is the reaeration rate, Ka. The reaction
schemes producing deficit are then defined, and the
corresponding reaction rate would be the product of
the stochiometric coefficient by the reaction rate
of the reaction using up oxygen. A second method
to compute deficit concentrations as done for com-
ponent NI in equations (22) thru (31), one obtains
[B](D)1-3.43[VK23](N2)+1.14[VK34](H3) (43)
where [B] is a matrix similar to the [A-j] matrices
of equations (31)'thru (35), except that the mafn
diagonal term has the reaeration rate Ka instead of
*11i (N2) and (N?) are nxl vectors of the steady-
state concentration of these reactants. (D)^ is an
nxl vector of the deficit concentrations over all
segments due to the oxydation of ammonia and nitrite.
The solution to the deficit concentration over all
space 1s given by:
(D)i=3.43[VK23](N2)[B]-1+1.14[VK34](N3)[B]-l.(44)
This method is used to compute deficit in program
FEDBAK, by using the optional subroutine.
The application to nitrification and dissolved oxygen
deficit assumed first order kinetics for the bacter-
ial reactions. This should be confirmed by laboratory
studies, or the nature of the system should be care-
fully considered. This computer model has been found
to be very useful as a predictive tool and in provi-
ding insights to the behavior of nitrogen species in
the aquatic environment. On new applications this
will ultimately depend on the applicatively of the
underlying assumptions to the system of interest.
ACKNOWLEDGEMENTS
The author is grateful to Steve Chapra of the Great
Lakes Environmental Research Laboratories and Richard
W1nf1eld of the Manhattan College Department of Envi-
ronmental Engineering and Science for their review
and comments on this paper.
The data for test applications of this model were made
available by Dr. Richard Tortoriello of the Delaware
River Basin Commission. His input to this project is
gratefully acknowledged.
REFERENCES
(1) O'Connor, D.J., "Oxygen Balance of an Estuary".
Jour. San. Enq. Div. ASCE. Vol 86, May 1960, pp 35-55
(2) O'Connor, D.J., "The Temporal and Spacial Distri-
bution of Dissolved Oxygen in Streams". Water Resour-
ces Research,-Vol 3, No. 1, 1967, pp 65-79.
(3) Chapra, S.C. and Gordimer S., ES001 A Steady Sta-
te. One Dimensional, Estuarine Hater Quality Model.
USEPA,Region II, New York, N.Y. September 1973.
(4) Thomann, R.V., "Mathematical Model for Dissolved
Oxygen". Jour. San. Eng.-Div., ASCE, Vol 89, No. SA5,
October 1963, pp 1-30.
(5) Thomann, R.V., System Analysis and water Quality
Management, Environmental Science Division, New York,
1971.
(6) Chapra, S.C. and Nossa, G.A. HAR03 A Computer
Program for the Modelling of Water Quality Parameters
in*Steady State Multidimensional Natural Aquatic Sys-
tem., USEPA, Region II, New York,'N.Y. October 1974
(7) Stratton, F.E. and Me Carty, P.L., "Prediction of
Nitrification Effects on the Dissolved Oxygen Balance
of Streams" Eny. Sci. and Tech.. Vol. 1, No. 5, May
1967, pp 405^4Tcn
(8) O'Connor, D.J., Thomann, R.V., and D1 Toro, D.M.,
Dynamic Water Quality Forecasting and Management.
USEPA, Office of Research and Development, Wash.,D.C.
Report No. EPA-660/3-73-009, August 1973
-------
-7-
DESCRIPTION OF THE MAIN PROGRAM AND
ACCOMPANYING SUBROUTINES
-------
-8-
Mt-v/RE fc.
R.£At>
CALL S«j
TO «»»
L
^x^
:D g«T{^.
,X1
C*t-<- iu«. DOOCf
To e.t»M«srtC BlF'tCtT
T*rAL "OCfi^lf /*>
BACU .3£cri£>A*
-------
-9-
The Main Program
FEDBAK03 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 various subroutines used in
the model. This inter-relationship among the various subroutines is illus-
trated in figure 6.
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 segments, indicators, etc.) and certain physical par-
ameters such as cross-sectional area, net advective flow, etc., which are
only input once when modeling a water body. These are in contrast to par-
ameters such as loads and reaction rates which would vary with the constitu-
ent 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 chloride concentrations (CHLOR), volume (VOL), and temperature
(TEMPER) at each segment
6) Apply the following conversion factors:
Q(MGD) = (CFS) * .6463 MGD/cfs
VOL(MG) = (106ft3) * 7.48 gal/ft3
-------
-10-
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) contains the area of interface
(I,IARAY(I,J))
ARRAY(I,J + 1) contains the E of interface
(I,IARAY(I,J))
ARRAY(I,J + 2) contains the Q of interface
(I,IARAY(I,J))
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 fills out
the remaining positions in ARRAY.
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.
-------
-11-
5 ubrout in e S ETUP . . .
The primary purpose of this subroutine is to set up the individual
flownet matrices (A) , (described in equations 31 through 35 in the theory
section) and output these on a direct access device. The tasks performed
by SETUP are:
1) Calculation of the bulk dispersion coefficient, E1
According to the theory:
E. . A. .
„•
E =
Where
E. . = the dispersion coefficient across the interface
!f J
A. . = the cross-sectional area of the interface
1 1 D
1. . = the average length of the adjoining
1'11 sections = (1. + l.)/2
The program calculates E ' (EPRIM) as follows :
EPRIM( J) =AKRAY (I , JJ) * ARRAY (I , JJ+1) *417 . 166/ (LA(K, I) +LA (I ,K) )
Where
ARRAY(I,JJ) = A. . (ft2)
1 / 3
ARRAY (I, JJ+1) = E. . (mi2/day)
1 1 J
LA(I,K) = 1.
LA(K,I) = lk
To conserve computer memory these lengths are stored as a N X 6
sized array and used in conjunction with array NLA to identify
the N,I interface. This again avoids the use of an N X N matrix
to handle this variable.
417.1166 - conversion factor which is computed
as follows :
Original Units _ Conversion Factor _ New Units
ft2*mile2 (5280 ft) 2 x 7.481 gal y 2 x MG MGD
day*ft*2 (mile) ^ ft3 10 gal
-------
-12-
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 12) and if it is not met, alpha is
recalculated.
Alpha(J) = 1.0-(EPRIM(J)/(2*ABS(ARRAY(I,JJ+2))))
Which corresponds to
o=l- E'/2Q
3) The elements of the system matrix (A) are calculated and saved on
mass storage
4) The values for EPRIM and ALPHA are printed
5) The section parameters: volume, temperature and chloride con-
centration are printed.
-------
-13-
SUBROUTINE SOURCE
This subroutine will read the reaction schemes and corresponding
reaction rate constants for each system segment and construct the (VK)
system matrices. It also reads all direct discharges and boundary con-
ditions, to form the system waste vector (W). These steps are accomplished
in the following sequence:
1. The reaction coefficients are read and are corrected to the tem-
perature in the section by the formula:
Ki - Ki 9(t-20)
For the nitrification reactions, values of 9 are discussed in
reference (2). The corrected reaction rates are then printed.
2. The (VK) matrix is formed by multiplying by the corresponding
section volumes and saved in mass storage.
3. The direct discharges of the various reactants into the system are
read and used to construct the system source vector. The contribu-
tion to the source vector due to the boundary conditions are calcu-
lated from the formulas given in equations (24) and (26).
4. The system source vector is then printed out.
-------
-14-
Subroutine MATRIX
This subroutine will:
(1) Combine the (VK) matrices on mass storage with the (A) matrices to
form the system matrix of matrices (A) given in equation (37)
(2) Optionally print the (A) matrix
(3) Read mile points corresponding to each of the segments
Subroutine MATN
This subroutine inverts the (A) matrix and multiplies it by the source
vector (W) giving the solution vector (N). This solution' vector represents
the component concentrations at the various segments. The printed output
includes the mile points corresponding to the beginning of each segment.
Subroutine SENS
This subroutine will perform system sensitivity analysis on the reaction
rate constants and/or direct discharges into the system. The logical se-
quence of this subroutine is:
For Changes in direct discharges;
1.) Read and print revised loadings for section
2.) Modify original source vector to reflect the new loadings
3.) Multiply new waste vector by the inverse of system matrix
4.) Print new solution vector
5.) Optionally, the deficit analysis subroutine KBOD may be called
For Changes in reaction rates;
For changes of direct discharges and reaction rates, the first two steps
above take place, the procedure below is then followed:
1.) The matrix of matrices (A) is re-assembled using the (A)and (VK)
matrices on mass storage files.
2.) The revised reaction constant(s) at all segments in the system are
read and temperature corrected. _
3.) The above values are substituted in the (A) matrix. (Note: Sub-
sequent calls to this subroutine make the changes to the original (A)
matrix and not from the one modified in this step.). For cummulative
effects this new (A) matrix should be copied to the mass storage file
containing the old (A) matrix._
4.) Optionally print the revised (A) matrix
This subroutine would then return to the main line program which in turns
calls the MATN and optionally the D.O. deficit analysis routine. Then if
desired it returns back to this routine for multiple runs of sensitivity
analysis.
-------
-15-
Subroutine DODEF
This subroutine will compute Dissolved Oxygen Deficit due to selected
reaction schemes. For each section of the system this routine will de-
termine the deficit due to each reactant and the resulting sum total
deficit in the system. The individual steps of this routine are:
(1) Read the reaction coefficients at each section, and temperature cor-
rect these by the expression:^
KA(I)=KA(I)*1.02488(TEMPER(I)-20.)
(2) The temperature and reaeration rates are printed
(3) The saturation level of dissolved oxygen at each system segment is
computed from the expression. ' '
CS (!) = (!. -.000009*CHLOR(I))*(14.652-.41022*TEMPER(I)+-
1.0079910*TEMPER(I)**2-.000077774*TEMPER(I)**3)
Where:
TEMPER(I) - is the temperature of the section
CHLOR(I) - is the chloride concentration (mg/1) of the section
(4) A system deficit matrix is formed by adding the VKa term to the
(A) matrix on the mass storage file (N.B. the VK term is omitted
for flexibility on the Flownet matrix (A) on mass storage)
(5) The above deficit matrix is inverted and optionally printed. The
structure of the program logic requires a separate array to f&rm the
system deficit matrix and compute its inverse. In the repeated use
of this subroutine during system sensitivity analysis, this inverse
of the deficit matrix is saved for subsequent use.
(6) Background contributions of deficit such as benthal demand or direct
sources are read into the section, the deficit source vector is
formed and multiplied by the inverse of the deficit matrix giving
the system deficit response due to direct loads. These results are
then printed.
(7) Next, the number of reactions and the reaction paths which produce
deficit as well as -the-"stochicmefcric coefficients are read.
(8) The (VKij) matrix for the deficit reactions is formed by accessing
the large (A) on mass storage.
(9) The deficit contribution to each reaction scheme is computed as
shown in equation (44) for the oxydation of ammonia.
(10) The resulting deficit and the dissolved oxygen concentration is
printed and accumulated for the calculation of total deficit.
-------
-16-
SOURCE PROGRAM LISTING
NOTE: The source code and the sample test case can be obtained in a
computer readable format for a nominal fee ($20.00 at present) by
writing to the office originating this document. This fee is waived
for federal agencies or local governmental planning organizations.
-------
-17-
******************************************************************
THIS PROGRAM SOLVES A SYSTEM OF FIRST ORDER CONSECUTIVE REACTION
WITH FEEDBACK IN A MULTI-DIMENSIONAL AQUATIC ENVIRONMENT
*** FEDBAK03 ***
THIS VERSION 2.0 CAN SOLVE UP TO A 60 SEGMENT SYSTEM, WHERE THE NO.
OF COMPONENTS * NUMBER OF SEGMENTS IS EQUAL OR LESS THAN 180
******************************************************************
INTEGER OUT
DIMENSION TITLEC20)
COMMON/BLK1/NR,N,NP,INPUT,OUT
COMMON/OPT I ON/IPRINT,ISENS,KWDEF,KBOD,NRUNS
DEFINE FILE 1(60,60,U,IREC1)
DEFINE FILE 3(180,180,U,IREC3)
*************************************************************** *
*
ASSIGN DEVICE NUMBERS TO CARD READER, LINE PRINTER ... *
*
******************************************************************
INPUT=8
OUT=6
IDODSW=0
READ*INPUT,200)TITLE
00 FORMAT (20A4)
N= NO. OF SEGMENTS (SECTIONS)
NR = NO. OF COMPONENTS (REACTANTS)
*
*
*
*
READ(INPUT,9999)N,NR
9 FORMAT(2I2)
WRITE (OUT,201) TITLE,N,NR
01 FORMAT(1H1,20X,20A4//28X,'NUMBER OF SECTIONS = »,13,10X,'NUMBER OF
1 REACTION COMPONENTS = f,I2,//)
NP=N*NR
READ PROGRAM OPTIONS
READUNPUT,1101)IPRINT,ISENS,K8QD,KWDEF
01 FORKAT(4I2)
******************************************************************
*
READ SYSTEM GEOMETRY — DISPERSIONS *
*
CALL DATA
******************************************************************
*
CHECK FLOW REGIMES *
*
CALL FLOW
******************************************************************
MAIN001
MAIN002
MAIN003
MAING04
MAIN005
MAIN006
MAIN007
MAIN008
MAIN009
MAIN010
MAINOll
MAIN012
MAIN013
MAIN014
MAIN015
MAIN016
MAIN017
MAIN018
MAIN019
MAIN020
MAIN021
MAIN022
MAIN023
MAIN024
MAIN025
MAIN026
MAIN027
MAIN028
MAIN029
MAIN030
MAIN031
MAIN032
MAIN033
MAIN034
MAIN035
MAIN036
MAIN037
MAIN038
MAIN039
MAIN04Q
MAIN041
MAIN042
MAIN043
MAIN044
MAIN045
MAIN046
MAIN047
MAIN048
MAIN049
MAIN050
MAIN051
MAIN052
MAIN053
MAIN054
MAIN055
MAIN056
MAIN057
MAIN058
MAIN059
-------
-18-
SET UP SYSTEM MATRIX
#
*
*
CALL SETUP
******************************************************************
*
READ AND PRINT LOADS + REACTION RATES — COMPUTE SOURCE VECTOR *
*
******************************************************************
CALL SOURCE
******************************************************************
*
NOW FORM LARGE MATRIX OF KV TERMS AND (A) MATRICES *
*
******************************************************************
CALL MATRIX
COMPUTE INVERSE OF LARGE MATRIX * SYSTEM CONCENTRATION VECTOR
*********a
CALL MATN
*********5
TEST FOR DISSOLVED OXYGEN DEFICIT ANALYSIS
IFCKBOD.GT.OJCALL DODEF
40 IF( ISENSUOO, 100,801
READ NUMBER OF SENSITIVITY RUNS DESIRED
********************:
101 READ( INPUT,30DNRUNS
101 FORMATU2)
102 CALL SENS
104 CALL MATN
510 IF(KBOD.GT.O)CALL DODEF
IF(NRUNS.GT.O)GO TO 302
LOO CALL EXIT
END
SUBROUTINE DATA
*
*
*
*
*
INPUTS DATA WHICH
CONSTITUENT BEING
HYDRODYNAMIC DATA
SUBROUTINE DATA
WOULD NOT VARY WITH THE
MODELED, I.E. PRIMARILY
PARTICULAR
PHYSICAL AND
*
*
*
*
*
*
*
>*********************************************************************
I*******:
INTEGER OUT
REAL LA(60,6)
DIMENSION AREA(6),E(6),Q(6)
COMMON/8LK1/NR,N,NP, INPUT,OUT
COMMON/BLK3/VOL(60),TEMPER(60),CHLOR(60)
MAIN060
MAIN061
MAIN062
MAIN063
MAIN064
MAIN065
MAIN066
MAIN067
MAIN068
MAIN069
MAIN070
MAIN071
MAIN072
MAIN073
MAIN074
MAIN075
MAIN076
MAIN077
MAIN078
MAIN079
MAIN080
MAIN081
MAIN082
MAIN083
MAIN084
MAIN085
MAIN086
MAIN087
MAIN088
MAIN089
MAIN090
MAIN091
MAIN092
MAIN093
MAIN094
MAIN095
MAIN096
MAIN097
MAIN098
MAIN099
MAINIOO
MAIN101
DATA001
DATA002
DATA003
DATA004
DATA005
DATA006
DATA007
DATA008
DATA009
OATA010
DATA011
DATA012
DATA013
DATA014
DATA015
DATA016
DATA017
-------
-19-
INTERFACE », 13, •-
13.
COMNON/COHOE/LA,SCALE(4),ARRAY(60,18)
COMMON/OPT I ON/1 PR INT,ISENS,KWDEF,KBOD
COMMON/UT1/IARAY(60,6)
COMMON/BLK4/NLA(60,6)
31 FORMAT(SFIO.O)
CO FORMAT14F7.0)
01 FORMAT (///28X, 'ZERO SEGMENT NUMBER IN
1« COMPUTATION DISCONTINUED')
00 FORMAT(6(I4,' -•,14,2X,F7.0,2X))
01 FORMAT(6(I3,F10.0))
00 FORMAT! »1
1 CHARAC. LENGTHS OF SEGMENTS (FT)
IFACE LENGTH INTERFACE LENGTH
2E LENGTH INTERFACE LENGTH')
11 FORMAT(3(3F6.0,I3))
02 FORMAT(//10X,'**SCALE FACTORS**'//10X,
l'SCALE(l) = ',F10.3,4X,'SCALE(2) = »,F10.3,4X,
2'SCALEt3) = •,F10.3,4X,«SCALE(4) = «,F10.3 )
00 FORMAT!»1«,2X,'INTERFACE',5X,'AREA',7X,'E',8X,«Q •,2(6X,'INTERFACE
l^SX,'AREA*,9Xt>E>*7XttQ*)/,4X,'ROW-SEG',4X,•(FT**2)', 2X,
2MMI**2/D)',2X,MCFS)', 6X, • ROW-SEG', 5X, • (FT**2 )',2X,
3 M MI**2/D>»,2X,MCFS)',5X,» ROW-SEG', 5X, ' ( FT**2 ) ' ,2X, ' (MI**2/D)' ,
4 2X,MCFS)» )
100 FORMAT* • (•,I3,f-',I3,M',2X,F8.0,4X,F6.3,F8.1,2(5X,M',I3,
1 '-'.IS,' )',2X^8.0,4X^6.3, F8.D)
•///' INTERFACE LENGTH
INTERFACE LENGTH
INTER
INTERFAC
READ AND WRITE THE SYSTEM PARAMETERS AND SCALE FACTORS
READ
-------
-20-
02 ARRAY(I,JJ+2)=Q(J)*SCALE(3)
******** $*:). $ $ ** * ******:**:$*;**#**:£**:$: * $**$#**********************
PROPAGATE AREAS,DISPERSIONS AND FLOWS FOR INTERFACES AT
-------
-21-
I**********$$********************************************************
IF(KBOD.GT.O)READ(INPUT,31)(CHLOR(K),K=1,N)
READtINPUT,31) tTEMPER!I),VOL(I),I=1,N)
IFITEMPER(2))14,13,14
.3 DO 14 1 = 1,N
TEMPERU)=TEMPER{1)
L4 CONTINUE
APPLY CONVERSION FACTORS TO FLOW AND VOLUME
DO 12 1=1,N
DC 23 J=U6
JJ=(J-1)*3.+ 1
>3 ARRAYU,JJ+2) = ARRAY(I,JJ+2)*.6463
L2 VOL{I)=VOL(I)*7.48
RETURN
END
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(60,6)
INTEGER OUT
COMMON/8LK1/NR,N,NP,INPUT,OUT
COMMON/COHOE/LA,SCALE(4),ARRAY(60,18)
NMN=0
SECTION ',13,' HAS AN EXCESS
FLOW
14 FORMATi1
!• CFS»)
02 FORMAT(*1«)
DO 213 1=1,N
QQ=ARRAYU,3)+ARRAY(I,6) + ARRAY(I,9)+ARRAYU,12)
1+ARRAY(I,18)
GQ=CQ/.6463
IF4QQ-.001)216,37,37
16 IF
-------
-22-
************^
REAL LA(60,6) SETUP011
INTEGER OUT SETUP012
DIMENSION EPRIM(6),ALPHA(6) SETUP013
COMMON/BLK1/NR,N,NP,INPUT,OUT,AB(60,60) SETUP014
COMMON/BLK3/VOL(60),T{60),CHLOR(60) SETUP015
CCMMON/COHOE/LA,SCALE(4),ARRAY(60,18) SETUP016
COMMON/OPTION/I PR INT,I SENS,KWDEF,KBOD SETUP017
COMMON/UT1/IARAY(60,6) SETUP018
COMMON/BLK4/NLA(60,6) SETUP019
23 FORMAT(6(I3,'-«,I3,F6.0,F6.3,2X)) SETUP020
24 FORMAT('i',49X, 'VALUES OF ALPHA AND EPR IM» ,///,6( • INTRFC EPRIM ALSETUP021
1PHA M,/,6(8X,'(MGD)',8X> } SETUP022
42 FORMAT(4X,I4,9X,F12.2,16X,F5.2,23X,F10.0) SETUP023
43 FORMAT('1',3X,'SECTION »,8X,'VOLUME (MGAL)',8X, SETUP024
1'TEMPERATURE (CENT.)',8X,'CHLORIDE CONC (MG/L)' ) SETUP025
44 FCRMAT('l*,3X,iSECTION »,8X,'VOLUME (MGAL)«,8X, SETUP026
1'TEMPERATURE (CENT.)') SETUP027
00 FORMAT(///,10X,'SEGMENT (•,13,*-•,13,') IS A BOUNDARY - BUT NO LENSETUP028
IGTH WAS SPECIFIED') SETUP029
01 FORMAT(//,10X,'PROGRAM TERMINATING DUE TO ',13,' BOUNDARY CONDITIOSETUP030
INS WITH MISSING LENGTHS'1 SETUP031
WRITE (OUT,24) SETUP032
IEXIT=0 SETUP033
DO 28 1=1,N SETUP034
DO 99 J=1,N SETUP035
99 AB(I,J)=0.0 SETUP036
*SETUP038
CALCULATE EPRIM AND ALPHA FOR EACH INTERFACE *SETUP039
*SETUP040
DO 20 J=l,6 SETUP042
EPRIM(J)=0.0 SETUP043
ALPHA(J)=0.0 SETUP044
JJ=(J-1)*3+1 SETUP045
K=IARAY(I,J) SETUP046
IF(K)20,20,10 SETUP047
*SETUP049
GENERATE (I,K) LENGTH ELEMENT *SETUP050
*SETUP051
10 XLEN1=0.0 SETUP053
DO 77 KK=1,6 SETUP054
IF(NLA(I,KK)-K)77,81,77 SETUP055
81 XLEN1=LA(I,KK> SETUP056
GO TO 82 SETUP057
77 CONTINUE SETUP058
!****«****************************************************************SETUP059
*SETUP060
GENERATE (K,I) LENGTH ELEMENT *SETUP061
*SETUP062
82 XLEN2=0.0 SETUP064
DC 78 KK=1,6 SETUP065
IF(NLA(K,KK)-I)78,84,78 SETUP066
84 XLEN2=LA(K,KK) SETUP067
-------
WRITE (VK) COEFFICIENTS
STORAGE FILE 3
-25-
FOR CORRESPONDING REACTION ON MASS
L01
95
100
DO 100 L=1,NS
IP=NS*(J-1)+L
READ(3'IP)(B(II),II=1,NP)
JP=NS*(I-l)-»-L
B(JP)=V(L)*K(L)
IF(I-J)101. 95,101
B(JP)=-B(JP)
WRITE(3'IPMBUC),IC=1,NP)
CONTINUE
GO TO 200
READ DIRECT COMPONENT DISCHARGES INTO THE SYSTEM
89 DO 61 JK=1,NP
61 W
-------
-26-
GO TO 3
SOURC105
GENERATE LA(K,K) LENGTH ELEMENT
*SOURC107
*SOURC108
*SOURC109
35
BOUNDARY CONDI
XLEN1=0.0
DO 77 KK=1,6
IF(NLAU,KK)-L)77,88,77
CONTINUE
WRITEtOUT,102)
CALL EXIT
FCRMAT{//,10X,'PRQGRAM TERMINATING DUE TO A
INS WITH MISSING LENGTHS')
88 XLEN1=LA
-------
-27-
DO 60 L=1,NS MATRX019
LL = NS *(J-1)+L MATRX020
»********************* jjc*^************************************** ****** *MATRX 021
*MATRX022
FORMATION OF THE LARGE MATRIX WHICH CONSISTS OF KV AND A MATRICES.*MATRX023
*MATRX024
60
93
20
L15
70
80
89
12
39
JOO
21
20
B(LL)=B(LL)+A(L)
DO 93 JC=l,NP
AC(JJ,JC)=B(JC)
CONTINUE
CONTINUE
IFUPRINT)90,90,70
WRITEIR,(AC(IR,JC),JC=1,NP)
FORMAT(IX,'ROW NO.',13,IX,10F12.2/<12X,10F12. 2))
CONTINUE
MATRX026
MATRX027
MATRX028
MATRX029
MATRX030
MATRX031
MATRX032
MATRX033
MATRX034
MATRX035
MATRX036
MATRX037
MATRX038
MATRX039
*MATRX041
PTMIL IS AN ARRAY CONSISTING OF THE MILEAGE POINTS OF THE SECTION*MATRX042
*MATRX043
90
>00
READ (INPUT, 200HPTMIL(I),I = 1,NS)
FORMAT (8F10.0)
RETURN
END
SUBROUTINE MATN
SUBROUTINE MATN
THIS SUBROUTINE INVERTS THE MATRIX AC AND MULTIPLIES
IT BY THE FORCING FUNCTION WS(I) GIVING SOLUTION VECTOR BU>
INTEGER OUT
COMMON/UT1/IPIVO(180),PIVOT(180)
CCMMON/UT2/INOEX<180,2)
COMMON/BLK1/N1,N2,N,INPUT,OUT,AC(180,180)
COMMON/OPTION/1 PR INT
COMMON/BLK4/PTMIL(60),B(180>
COMMON/BLK5/WS(180)
DO 900 1=1,N
DO 900 J=1,N
AC(ItJ)-AC(I»J)/1000000.
DO 20 J=1,N
IPIVO(J)=0
DO 21 1=1,2
INDEXfJ,I)=0
CONTINUE
DO 550 1=1,N
SEARCH FOR PIVOT ELEMENT
MATRX045
MATRX046
MATRX047
MATRX048
MATN001
MATN002
MATN003
MATN004
MATN005
MATN006
MATN007
MATN008
MATN009
MATN010
MATN011
MATN012
MATN013
MATN014
MATN015
MATN016
MATN017
MATN018
MATN019
MATN020
MATN021
MATN022
MATN023
MATN024
MATN025
MATN026
MATN027
MATN028
MATN029
-------
-28-
AMAX=0.0
DO 105 J=1,N
IF (IPIVO (J)-l)
SO
BO
85
DO
05
50
00
60
CO
50
150
>30
105
740
901
60, 105, 60
DO 100 M=1,N
IF (IPIVG (M)-l) 80, 100, 740
IF (ABS(AMAX)-ABS(AC(J,M)))85,
IROW=J
ICOLU =M
AMAX=ACU,M)
CONTINUE
CONTINUE
IPIVO (ICOLU )=IPIVO (ICOLU
100, 100
INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
IF UROW-ICOLU > 150, 260, 150
DO 200 L=1,N
SWAP=AC(IROW,L)
AC(IROW,L)=AC(ICOLU,L)
AC(ICOLU,L)=SWAP
INDEX(I,1)=IRQW
INDEX(I,2)=ICOLU
PIVOT(I)=AC( ICOLU, ICOLU)
DIVIDE PIVOT ROW BY PIVOT ELEMENT
AC(ICOLU,ICOLU)=1.0
DO 350 L=1,N
50 AC ( ICOLU, L) = AC( ICOLU, L)/PIVOT( I)
REDUCE NON-PIVOT ROWS
DO 550 L1=1,N
IF(L1-ICOLU ) 400, 550, 400
S=AC
-------
-30-
30 READUNPUT,31)LOAD,ISEG,IREACT
31 FORMAT(F10.0,212)
IF
-------
-31-
2 'THETA',5X,13I8,/,4(13X,13I8,/M
DC 115 1=1,NS
READ(1»IJROW
DC 29 J=1,N
JJ=NS*(J-1)+ I
READ(3«JJ)ROW2
DO 60 L=1,NS
LL = NS *{J-1)+L
FORMATION OF THE LARGE MATRIX WHICH CONSISTS OF KV AND ROWMATRICES*
*
60 ROW2(LL)=RCW2(LL)+ROWU)
DO 93 JC=1,NP
29 CONTINUE
15 CONTINUE
DO 25 I=1,NRATES
READdNPUT,2>II,KK,THETA
2 FORMAT(2I2,1X,F5.3)
READ REVISED DECAY RATE KIII.KK) FOR ALL SEGMENTS
IFdI.GT.O.AND.KK.GT.OJGO TO 21
WRITE(OUT,5)
5 FOfcMAT(« PROGRAM TERMINATED DURING SYSTEM SENSITIVITY ANALYSIS
1DUE TO ERRORS IN REACTION RATES DATA INPUT' )
CALL EXIT
21 READdNPUT,3)(K(L),L=l,NS)
3 FORMAT(SFIO.O)
TEMPERATURE CORRECT NEW DECAY CONSTANTS
IF(THETA.LE.O)GO TO 82
DO 83 L=ltNS
83 K(L)=K-20)
82 WRITE(OUT,4)II,KK,THETA,(K
-------
-32-
13
16
>5
ro
30
12
39
30
99
NOTE — PROGRAM ASSUMES CHANGES ARE MADE FROM *** ORIGINAL ***
SYSTEM MATRIX, IF CUMMULATIVE EFFECTS ARE DESIRED
SAVE THIS NEW SYSTEM MATRIX ON OFFLINE DISK FILE 1
IROW=IROW+1
ICOL=ICOL-H
CONTINUE
CONTINUE
NRUNS=NRUNS-1
IF(IPRINT)90,90,70
WRITEtOUT.80)
FORMAT(1H1,///,43X,'REVISED SYSTEM MATRIX BEFORE SOLUTION*,///)
WRITE(OUT,89)(I,1=1,NP)
FORMAT(/,12X,10112,/)
DO 39 IR=1,NP
WRITE(OUT,12)IR,(AUR,JC),JC=1,NP)
FORMAT(IX,•ROW NO.•,13,IX,10F12.2/<12X,10F12.2)>
CONTINUE
RETURN
CALL EXIT
END
SUBROUTINE DODEF
* SENS155
* SENS156
* SENS157
* SENS158
SENS159
SENS160
SENS161
SENS162
SENS163
SENS164
SENS165
SENS166
SENS167
SENS168
SENS169
SENS170
SENS171
SENS172
SENS173
SENS174
SENS175
SENS176
OODEF001
*DODEF002
THIS SUBROUTINE WILL COMPUTE DISSOLVED OXYGEN DEFICIT
BY SELECTING COMPONENTS CAUSING DEFICIT WITHIN THE SYSTEM
*DODEF004
*DODEF005
*DODEF006
THE OUTPUT WILL BE D.O. DEFICIT DUE TO EACH COMPONENT AND *DODEF007
TOTAL DEFICIT DUE TO ALL COMPONENTS IN EACH SECTION OF THE SYSTEM*DODEF008
10
INTEGER OUT
REAL KA
DIMENSION STOCH{10»,ICOEF1(10),ICOEF2(101,BMAT(60,60),CS(60)
DIMENSION VK(60),RES(60),PROD(60),DODSUM(60),SROW(60),WDODEF(60)
COMMON/UT1/IPIVO(180),PIVOT(180)
COMMON/UT2/INDEX(180,2),KA(60)
COMMON/BLK1/NR,N,NP,INPUT,OUT
COMMON/BLK3/V(60),TEMPERC60),CHLOR160)
COMMON/BLK4/PTMIL(60),XK1(180)
COMMON/OPTION/IPRINT,ISENS,KWDEF
EQUIVALENCE
-------
-33-
*DODEF038
IF KA ARE SAME IN ALL SECTIONS, INPUT ONLY VALUE OF FIRST SECTION*OODEF039
*DODEF040
******** ************#*#£*****$**$**$*$*$*****************************DODEF041
IF(KA(2))2t2,3 DODEF042
2 DO 4 1=2,N DODEF043
KA(I)=KA(1) DODEF044
4 CONTINUE DODEF045
*********************************************************************QQDEF046
*DODEF047
TEMPERATURE CORRECT REAERATION RATES . *DODEF048
*********************************************************************OODEF050
DO 6 I=l,N DODEF051
KA(I)=KA(I)*1.024**(TEMPER(I)-20.) DODEF052
6 CONTINUE DODEF053
***********************************#*********************************OQDEF054
*DODEF055
PRINT REAERATION RATES *DODEF056
*DODEF049
*DODEF057
*********************************************************************DQQEFO58
3 WRITE(OUT,8) DOOEF059
8 FORMATC1',40X,'DISSOLVED OXYGEN DEFICIT ANALYSIS',//,5X, DODEF060
l'SEGMENT',10X,'TEMPERATURE (CENT.»»,8X,'REAERATION RATE (I/DAY)' JDODEF061
WRITE(OUT,9J{I,TEMPER(I),KA{I),I=1,N) DODEF062
9 FORMAT(5X,I4,14X,F6.2,21X,F12.6) DODEF063
*********************************************************************DQDEF064
*DODEF065
COMPUTE SATURATION LEVEL OF DISSOLVED OXYGEN IN EA. SECTION *DODEF066
*DQDEF067
DO 7 1 = 1,N DODEF069
CS(I)=(l.-.000009*CHLORtI))*<14.652--41022*TEMPERU)+ DODEF070
1.0079910*TEMPER(IJ**2.-.000077774*TEMPER(I)**3.) DODEF071
7 CONTINUE DODEF072
COMPUTE DEFICIT MATRIX *DOOEF073
DO 101 1=1,N DODEF074
READ(1§ IMBMATdf J)v J-ltN) DODEF075
BMATCI,I)=BMAT(I,I)+V(I)*KA(I) DODEF076
01 CONTINUE DODEF077
*DODEF079
PRINT DISS. OXYGEN DEF. MATRIX IF DESIRED *DODEF080
*DODEF081
*********************************************************************DOOEP082
IF(IPRINT)90,90,70 DODEF083
70 WRITECOUT,280) DODEF084
80 FORMATdHl,///,SOX,'DEFICIT MATRIX BEFORE INVERSION',///) DODEF085
WRITE(OUT,89)tI,I=l,N) DODEF086
89 FORMAT(/,12X,10112,/) DODEF087
DO 39 IR=1,N DODEF088
WRITE(OUT,162)IR,
-------
-34-
10 DC 900 1 = 1,N
DC 900 J=l,N
10 8MATU, J) = 8MATU,J)/1000000.
DC 20 J=1,N
IPIVO(J)=0
DO 21 1=1,2
>1 INDEX(J,I)=0
!0 CONTINUE
DO 550 1=1,N
SEARCH FOR PIVOT ELEMENT
AMAX=0.0
00 105 J=1,N
IF (IPIVO (J)-l)
60, 105, 60
iO DO 100 M=1,N
IF (IPIVO (M)-l) 80, 100, 740
30 IF (ABS(AMAX)-ABS(BMAT(J.M)))85, 100, 100
35 IROW=J
ICOLU =M
AHAX=BMATCJ,M)
30 CONTINUE
35 CONTINUE
IPIVO UCOLU )=IPIVO (ICOLU )+l
INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
IF UROW-ICOLU ) 150, 260, 150
50 DO 200 L=1,N
SWAP=BMAT(IROW,L)
BPAT(IROW,L)=BMAT(ICOLU,L)
00 BMAT(ICOLU,L)=SWAP
60 INDEX(I,1)=IRQW
INDEX(I,2)=ICOLU
PIVOT(I)=BMAT(ICOLU,ICOLU)
DIVIDE PIVOT ROW BY PIVOT ELEMENT
BMATlICOLU,ICOLUJ=1.0
DO 350 L=1,N
50 BMAT{ICOLU,L)=BMAT(ICOLU,L)/PIVOT(I)
REDUCE NON-PIVOT ROWS
DC 550 L1=1,N
IFCL1-ICOLU ) 400, 550, 400
00 S=BMAT(L1,ICOLU)
BMAT(L1,ICOLU)=0.0
DO 450 L=1,N
50 BMAT(L1,L)=BMAT(L1,L)-BMAT(ICOLU,L)*S
50 CONTINUE
INTERCHANGE COLUMNS
DO 710 1=1,N
L=N+1-I
IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630
i30 JROW=INDEX(L,1)
JCOLU =INDEX(L,2)
DQDEF097
DODEF098
DODEF099
DODEF100
DODEF101
DODEF102
DODEF103
DODEF104
DODEF105
*DODEF106
*DODEF107
*DODEF108
DODEF109
DODEF110
DODEF111
DODEF112
DODEF113
DODEF114
DODEF115
DODEF116
DODEF117
DODEF118
DODEF119
DODEF120
*DODEF121
*DODEF122
*DODEF123
DODEF124
DODEF125
DODEF126
DODEF127
OODEF128
DODEF129
DODEF130
DODEF131
*DODEF132
*DODEF133
*DODEF134
DODEF135
OODEF136
DOOEF137
*DODEF138
*DODEF139
*OODEF140
DODEF141
DODEF142
DODEF143
DODEF144
DODEF145
DODEF146
DODEF147
*DODEF148
*OQDEF149
*DODEF150
DQDEF151
DODEF152
DODEF153
DODEF154
DODEF155
-------
-35-
00 705 M=1,N
SWAP=BMAT(M,JROW)
BNAT(M,JROW)=BMAT(M,JCOLU)
B!«AT(M,JCOLU)=SWAP
705 CONTINUE
710 CONTINUE
740 CONTINUE
DO 901 1=1,N
DO 901 J=1,N
901 BMATU,J) = BMATU, JJ/1000000.
IOODSW=999
PRINT INVERSE OF DEFICIT MARTIX IF DESIRED
IF(IPRINT)139,139,370
370 WRITE(OUT,281)
281 FORMAT(1H1,//,40X,«INVERSE OF DEFICIT MATRIX1,///}
WRITE(OUT,89X1,1=1,N)
DO 339 IR=1,N
WRITE
-------
-36-
ICOL=(ICOEF1UI)-1)*N+1
LAST=IROW+N-1
ISTRT=1
00 104 J=IROW,LAST
READ < 3' J)(SROW(K),K=lfNP)
VK(ISTRT) = ABS(SROW(ICOD)
ISTRT=ISTRT-H
ICQL=ICOL+1
.04 CONTINUE
MULTIPLY (VK)* *STOCH (II) *VK ( I )
LOC=LOC+1
L55 CONTINUE
MULTIPLY CVK)*STOCH*(SOL) BY INVERSE OF DEFICIT MATRIX
DO 106 1=1, N
PRODI I)=0.0
DO 106 J=1,N
PROD ( I )=PRCD ( I KBMAT ( I , J )*RES i J )
106 CONTINUE
PRINT RESULTING DEFICIT FOR REACTION
WRITE (OUT, 12 ) ICOEF1 (II), ICOEF2( 1 1 ) ,STOCH( 1 1 J
12 FORMAT!//, 20X, 'DEFICIT DUE TO REACTION* , 12, • TO ',12,
1 10X,fSTOCHIOMETRIC COEFFICIENT OF • ,F6.2, //,5X,
1 'SECTION*, 5X,»MILE PTf,5X, 'DEFICIT CONCENTRATION (MG/L)',5X,
1 fD.O. SATURATION {HG/L )', 5X, 'DISSOLVED OXYGEN (MG/L)r )
DO 108 1=1, N
DO=CS(IJ-PROD(I)
WRITE(OUT,13)I,PTMIL! I ) ,PROD( I ),CS( I ),00
108 CONTINUE
13 FORMAT(7X,I3,8X,F5.1,13X,F7.3,25X,F7.3,19X,F7.3)
SUM REACTION DEFICIT FOR TOTAL DEFICIT COMPUTATION
DC 107 1=1, N
DODSUM(I)=DODSUM(I)+PRCDII )
107 CONTINUE
103 CONTINUE
PRINT TOTAL DISSOLVED OXYGEN DEFICIT
DODEF229
DODEF230
DODEF231
DODEF232
DODEF233
'DODEF234
*DODEF235
*OODEF236
*DODEF237
^DODEF238
DODEF239
DODEF240
DODEF241
DODEF242
DODEF243
'DODEF244
*DODEF245
*DODEF246
*DODEF247
*DQDEF248
DODEF249
DODEF250
DODEF251
DODEF252
DODEF253
DODEF254
DODEF255
DODEF256
DODEF257
DODEF258
'DODEF259
*DODEF260
*DODEF261
*DODEF262
CDODEF263
DODEF264
DODEF265
DODEF266
DODEF267
'DODEF268
*DQDEF269
*DODEF270
*DOOEF271
WRITE(OUT,14)(ICOEF1(I),ICOEF2(I), I=1,NRDEF)
DODEF273
-------
1. TITLE CARD
-39-
COLUMNS
1-80
VARIABLE
Title
FORMAT
20A4
SYSTEM SIZE CARD
COLUMNS
1- 2
3- 4
PROGRAM OPTIONS
variable.
COLUMNS
1- 2
3- 4
5- 6
7- 8
SYSTEM GEOMETRY
COLUMNS
1- 7
8-14
15-21
22-28
VARIABLE
N
NR
- If option
VARIABLE
IPRINT
ISENS
KBOD
KWDEF
CONTROL CARD
VARIABLE
SCALE (1)
SCALE (2)
SCALE O)
SCALE (4)
FORMAT
12
12
desired, enter
FORMAT
12
12
12
12
FORMAT
F7.0
F7.0
F7.0
F7.0
DESCRIPTION
Users description of run
DESCRIPTION
Number of system segments
Number of components
(reactants)
a positive integer in corresponding
DESCRIPTION
Print system matrix and its
inverse
Perform sensitivity analysis
of system
Perform dissolved oxygen defi-
cit analysis
Control variable to read
direct (background) sources of
deficit in the deficit anal-
ysis subroutine
DESCRIPTION
Area scale factor ft /ou*
Dispersion scale factor
mi /day/ou
Flow scale factor CFS/ou
Length scale factor ft/ou
*note: ou stands for the "original units" in which the parameter to be con-
verted 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
-------
-40-
5.' INTERFACE PARAMETER CARDS (_2 per segment; total number = 2*N)
Format Description
Cross sectional area of
F6.0 interface with, section
Dispersion coefficient of
F6.0 first interface with section
Flow- * across the first
F6.0 interface with section
Section forming the first
13 interface with section
Cross sectional area of
F6.0 second interface
Dispersion coefficient of
F6.0 second interface
Flow across the second
F6.0 interface
Section forming the second
13 interface
Cross sectional area of third
F6.0 interface
Dispersion coefficient of
F6.0 third interface
Flow across the third inter-
F6.0 face
Section forming the third
13 interface
Note: The second interface parameter card is identical to
above, but for the fourth, fifth and sixth interfaces.
If the segment 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 segment 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 segment 1
with segment 2, it is not necessary to input the
parameters of the interface of segment 2 with segment 1.
INTERFACE
Columns
1- 6
7-12
13-18
19-21
22-27
28-33
34-39
40-42
43-48
49-54
55-60
61-63
PARAMETER (
Variable
AREA
E
Q
IARAY
AREA
E
Q
IARAY
AREA
E
Q
IARAY
Units
ft
mi /day
cfs
ft
mi /day
cfs
ft
mi /day
cfs
* flow out of a segment is positive;
negative.
flow into a segment is
-------
-41-
6.
CHARACTERISTIC LENGTHS (1 card per section; 1 length for each interface;
therefore up to 6 lengths per card)
Columns Variable
1- 3
4-13
14-16
17-26
27-29
30-39
40-42
43-52
53-55
56-65
66-68
69-78
NLA(I,JI)
LA(I,JI)
NLA(I,JJ)
LA(I,JJ)
NLA(I,JK)
LA(I,JK)
NLA(I,JL)
LA(I,JL)
NLA(I,JM)
LA(I,JM)
NLA(I,JN)
LA(I,JN)
13
F10.0
13
F10.0
13
F10.0
13
F10.0
13
F10.0
13
F10.0
Format Description
First segment which forms an
interface with segment I
Length of segment I with
respect to section JI*
Second segment which forms an
interface with segment I.
Length of segment I with
respect to segment JJ
Third segment which forms
an interface with segment I
Length of segment I with
respect to segment JK
Fourth segment which forms
an interface -with segment I
Length of segment I with
respect to segment JL
Fifth segment which forms
an interface with segment I
Length of segment I with
respect to segment JM
Sixth segment which forms
an interface with segment I
Length of segment I with
respect to segment JN
Units
feet
feet
feet
feet
feet
feet
* note: the following figure illustrates the interpretation of
characteristic lengths as used in this program.
JN
JJ
JK
JI
JL
JM
Figure 7: Depiction of characteristic lengths
-------
-42-
Where
@ - LA(I,JK) and LA(I,JI)
(D - 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.
Note also that for sections which form boundaries across which
mass is transported, the characteristic length of the segment with respect
to the boundary must be inputted. This length is read in as LA(I,I), where
I is the boundary section. (see equations 23 and 25).
7. SECTION CHLORIDE CONCENTRATIONS. Optional. This data is only necessary if
KBOD>1 (see card type #3)
Columns Variable Format
1-10
11-20
CHLOR(l)
CHLOR(2)
F10.0
F10.0
Description
Chloride concentration (mg/1) at the
1st section
Chloride concentration at the second
section
71-80 CHLOR(8) F10.0 Chloride concentration (mg/1) at the
8th section
The card above is repeated as necessary until all chloride concentrations
for the system have been entered.
8. TEMPERATURE/VOLUME CARDS
Columns Variable Format Description
1-10 TEMPER(l) F10.0 Water temperature (°C) of 1st section
11-20 VOL(l) F10.0 Volume (10 ft3) of 1st section
21-30 TEMPER(2) F10.0 Water temperature (°C) of 2nd section
31-40 VOL(2) F10.0 Volume (106 ft3) of 2nd section
71-80
VOL(4)
F10.0
Volume of 4th section
This card is repeated until the temperature and volume for all sections have
been entered. If the temperature is constant in the system, it is only
necessary to enter this parameter for the 1st section, as the program will
then test for the 2nd section temperature and since it is blank it will
propagate the value of the 1st section to all the sections in the system.
-------
9.
-43-
REACTION RATE COEFFICIENTS
REACTION SCHEME CARD.
Consider the reactions:
-f D
These reactants are represented numerically in the program as follows:
The reaction rates will have indices referring to these reactant numbers:
K
K14
The program will read these reaction indices as follows:
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-2 I
3-4 J
6-10 THETA
FORMAT
12
12
F 5.3
1st reaction index of reaction rate
constant Kij
2nd reaction index of reaction rate
constant Kij
temperature correction coefficient; the
program assumes the rate constants are
given at 20°C and can be converted to
any temperature by the equation:
Kij(T°C) = Kij(20°C) * THETA**(TEMP-20).
If this field is left blank, the program
will not correct the reaction rates for
temperature.
10. REACTION RATE CONSTANTS (Units: I/day)
COLUMNS
1-10
11-20
21-30
VARIABLE
NAME
K(l)
K(2)
K(3)
FORMAT
F10.0
F10.0
F10.0
DESCRIPTION
reaction rate constant for 1st segment
reaction rate constant for 2nd segment
reaction rate constant for 3rd segment
71-80
K(8)
F10.0
reaction rate constant for 8th segment
-------
-44-
The card image above is repeated until the reaction rate constant Kij;
where i, j are the indices specified on card #9, have been read. The
reaction rates for a new reaction are then entered by preparing a new #9
card with the appropriate indices for the new reaction, then reaction
rate cards follow with the reaction rate constant for all segments.
This procedure is repeated until the reaction rate coefficients have been
read for all reaction schemes. A blank card is inserted after the last
set of reaction coefficients, which signals an end to this type of input.
11. DIRECT DISCHARGES
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-10 LOAD F10.0 Waste load into segment (Ib/day)
11-12 ISEC 12 segment number where waste load is
being discharged
13-14 IREAC 12 component (number) making up waste
load
These cards are repeated for all point sources into the system, the input
is terminated by a blank card.
12. BOUNDARY CONDITION
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-10 BC F10.0 Boundary concentration (mg/1)
11-12 ISYS 12 reaction component number
13-14 ISECT 12 section number
These cards are repeated for each reaction component at all boundary
segments. The last card will again be a blank card to terminate this type
of input.
13. MILEAGE POINT OF SEGMENTATION
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-10 PTMIL(4) F10.0 Mile point designation for interface at
end of first segment
11-20 PTMIL(2) F10. mile point at end of second segment
71-80 PTMIL(8) F10. mile point at end of eighth segment
-------
-45-
These cards are again repeated if necessary to input remaining mile points.
The program will read the same number of mile points as the number of seg-
ments defined for the system on card image #2.
This mile point is used only to locate the final concentrations in the system;
it will not be used for any numerical computation such as the length of a
segment.
DISSOLVED OXYGEN DEFICIT ANALYSIS (KBODX) ON OPTION CARD)
14. SEGMENT REAERATION RATE
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-5 KA(1) F5.3 Reaeration rate (I/day) for first segment
6-10 KA(2) F5.3 Reaeration rate for second segment
76-80 KA(16) ' F5.3 Reaeration rate for sixteenth segment
This card image is repeated if necessary to read reaeration rate constants
for remaining segments.
15. BACKGROUND DEFICIT DUE TO BACKGROUND SOURCES (KWDEF> 0 ON OPTION CARD)
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-10 WDODEF(l) F10.0 Direct discharge of deficit (Ib/day) into
the 1st segment
11-20 WDODEF(2) F10.0 Direct discharge of deficit (Ib/day) into
the 2nd segment
71-80 WDODEF(8) F10.0 Direct discharge of deficit (Ib/day) into
the 8th segment
The above card is repeated until these direct discharges have been defined
for the whole system.
-------
-46-
1-6. SELECTING REACTION SCHEMES PRODUCING DEFICIT
COLUMNS
1- 5
11-12
13-14
15-20
VARIABLE
NAME
NRDEF
ICOEFl(l)
ICOEF2(1)
STOCK(1)
FORMAT
15
12
12
F6.0
DESCRIPTION
Number of reactions producing deficit in
the system
1st index of reaction rate constant
IJ
2nd index of reaction rate constant
Stochiometric coefficient used in converting
component concentration to equivalent
oxygen utilized
The subscripting of the "variable names" in the above cards refers to the
first reaction producing deficit, the same procedure is utilized for all
other deficit producing reactions:
21-22
23-24
25-30
ICOEF1(2)
ICOEF2(2)
STOCK (2)
12
12
F6.0
1st index of reaction rate constant K
2nd index of reaction rate constant K.
Stochiometric coefficient for second
reaction
IJ
'IJ
71-72 ICOEFK7) 12
73-74 ICOEF2(7) 12
75-80 STOCH(7) F6.0
1st index of reaction rate constant K
IJ
2nd index of reaction rate constant K
Stochiometric coefficient for seventh reaction
NOTE: The reaeration rates constants and the reaction schemes producing
deficit are read only once in the program. During consecutive runs; such
as in system sensitivity analysis, there is an internal switch in the pro-
gram to bypass this input and also the computation of the system deficit
matrix inverse after the initial use of the DODEF subroutine.
-------
-47-
5YSTEM SENSITIVITY ANALYSIS
(ISENSX) ON OPTION CARD)
17. MAIN PROGRAM CONTROL VARIABLE
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-2 NRUNS 12 Number of sensitivity analysis runs to be
performed
18. SUBROUTINE CONTROL VARIABLE
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-2 NRATES 12 Number of reaction rate schemes to be
revised
3-4 NLOADS 12 Number of loads in system source vector
to be changed
19. CHANGES TO WASTE LOADS (NLOADS > 0)
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1-10 LOAD F10.0 New waste load (Ib/day) into segment
11-12 ISEG 12 Segment number of modified waste load
13-14 IREACT 12 Component number of waste load
This card image is repeated until the number of new waste loads specified
on card #18' as NLOADS have been read.
CHANGES TO REACTION RATE COEFFICIENTS (NRATES > 0)
20. REACTION RATE SCHEME TO BE MODIFIED
VARIABLE
COLUMNS NAME FORMAT DESCRIPTION
1- 2 II 12 1st index of reaction rate Kn, „,
J.X f KK
3- 4 KK 12 2nd index of reaction rate K
X J. f JS-K.
6-10 THETA F5.3 Reaction rate temperature correction
coefficient
-------
21. NEW REACTION RATE CONSTANTS
VARIABLE
COLUMNS NAME
FORMAT
-48-
DESCRIPTION
1-10
11-20
K(2)
F10.0 New reaction rate constant K
for 1st segment.
F10.0 New reaction rate constant K
segment
11,KK (see above)
11,KK for 2nd
71-80
K(8)
F10.0 New reaction rate constant K
segment.
11,KK for 8th
This card image is repeated until the reaction rate constant specified in
card #20 as K has been read in for all segments,
iX f KK
The card image #20 and the sequence of card image #20 is repeated until the
number reaction schemes specified on card #19 (NRATES) has been satisfied.
A new sensitivity run (if NRUNS > I) would begin by reading a new card image
#18 (see schematic diagram for various inputs).
-------
-49-
KESTRICTIONS
1.) FEDBAK03 as presently written is limited to a system of sixty sections
2.) Each section can have a maximum of six interfaces.
3.) Only one interface of a particular section may act as a boundary.
4.) The maximum number of reactants is a variable, which when multiplied
by the number of sections, it cannot exceed 180.
Restrictions (1) and (4) can be easily expanded by minor modifications
to the computer program. The application of this model to the nitro-
gen cycle by assuming first order kinetics of the bacterial reactions
is a major assumption which should be verified by laboratory studies.
-------
-50-
APPLICATION TO THE DELAWARE ESTUARY
-------
-51-
The area modelled in the Delaware Estuary is the 86 mile stretch of tidal
river between Trenton, New Jersey and Liston Point, Delaware. Figure 8.
shows this area and the corresponding system segmentation which was
followed.
This area was extensively investigated during the 1960's by the Delaware
Estuary Comprehensive Study. The purpose of this study was to determine
the water quality objectives for the estuary, the cost and benefits of
these objectives, and the treatment levels necessary to achieve these
objectives. The survey data used in this test case was gathered by Hydro-
science Inc. for the Delaware River Basin Commission. The intent of this
section is to demonstrate the user of this model its application by re-verify-
ing the results published in the forementioned report. Interested readers are
suggested to review the Hydroscience report for an in depth discussion on the
rational for the reaction coefficients, as well as the differences between
computed and observed values.
-------
-52-
DELAWARE ESTUARY
COMPREHENSIVE STUDY
SECTIONS FOR
MATHEMATICAL MODEL
U. S. ENVIRONMENTAL PROTECTION AGENCY
Region II New York. New York
FIGURE
-------
TABLE 2
Source
Del River at Trenton
Trenton
Philadelphia N. E.
Camden
Philadelphia S. E.
Texaco
Philadelphia S. W.
Mobil
Repauno
C.D.C.A.
Chester
Sinclair
Sun Oil
Wilmington
Carneys Pt.
Chambers
POINT SOURCE
DRBC
River Mile
134.0
132.0
104.2
97.9
96.6
93.9
90.7
87.5
85.5
85.0
80.8
80.6
78.8
72.2
71.2
69.8
NITROGEN LOADINGS
Seg
No
1
1
9
13
14
15
16
17
17
17
18
18
19
21
21
22
Organic N
(#N/day)
9,000
740.
6,620.
2,280.
4,210.
v-
4,320.
.-•
1 — :-
1,000.
1,000.
8,250.
TO THE DELAWARE ESTUARY
Ammonia-N Nitrite N
(#N/day) (#N/day)
1,000.
7,020. 10.
8,930. 10.
1,870. 330.
9,230.
1,000.
13,780.
3,600.
18,000.
1,230. 40.
1,270.
2,000.
7,900.
• 5,100. 250.
Nitrate N Total N
(#N/day) (#/day)
16,000 26,000
40. 7,810.
110. 15,670.
700. 5,180.
13,440.
1,000.
18,100.
3,600. i
Ul
U)
18,000. '
50. 2,320.
2,270.
2,000.
7,099.
800. 14,400.
3.0,500. ) 30,500.
' >
In order to reflect nitrogen inputs from natural runoff and smaller tributaries, it was assumed that 1,0 mg/1
ammonia nitrogen were contributed along the length of the estuary with an associated flow of 20 CFS/miles,
This is equivalent to a uniform loading of 107 Ibs/day mile" as shown on the system waste vector of the computer
program.
-------
NUN'BCK OF SECTIONS = 30 NUMBER OF REACTION COMPONENTS =
**SC~ALE FACTORS**~~
rr^nnrcrm sraxrnnr^ rrc'c'c rrsrirrTr^ tTvvv STTSTIIA) =—rcuo.oco
-------
SECTION
SECTION
SECTION
SECTION.
SECTION
. _. SECTION
SECTION
SECTION-
SECTION
SECTION
SECTION
SECTION
'SECTION
SECTION.,
SECTION
SECTION.
SECTION
SECTION
SECTION
SECTION_
S'ECTION
SECTION.
SECTION
SECTION.
SECTION
S.ECTIJQN_
SECTION
SECTION
S'ECTION
1 HAS AN EXCESS FLOW
_2__H A S _AN_£XCESS FLOW
3 HAS AN EXCESS FLOW
Oh 0.303000^ 03 CFS
0.126000E 03 CFS
OF 0.108000E 03 CFS
4_HA$. AN. EXCESS. FLOW Of _Q.999983£ 01_CFS_
5 HAS AN EXCESS FLOW OF 6.330000E 02 CFS
6_HAS. AH. EXCESS^ FLOW OF..D.9000QOE.D2 CF$._
7 HAS AN EXCESS FLOW OF -.130000E 03 CFS
9 .H6_S_JN EXCESS FLOW OF QT?sono3FL_PZ_rFS_
9 HAS AN EXCESS FLOW OF 0.319997E 02 CFS
..10_HAS__AN_EXCESS__E.LQ.H D.F__Qt200000E_.D2_CFS_
11 HAS AN EXCESS FLOW OF 0.230000E 03 CFS
12 HAS AN EXCESS FLOW OF 0.459999E 02 CFS
OF O.A99992E 01 CFS
J)F 0.100002E 02 CFS
OF 0.649996E 02 CFS
_QF__0.7090POE_03 CFS_
OF 0.369000E 03 CFS
_QF__jO. 1130QOE 03 CFS_
OF 0.749999E 02 CFS
J)F 0.320001F 02 CFS
13 HAS AN EXCESS FLOW
_14_HAS_AN_EXCESS FLOW
15 HAS AN EXCESS FLOW
_16 HAS AN exCESS_FigW_
17 HAS AN EXCESS FLOW
18 HASSAN EXCESS FLOW
19 HAS AN EXCESS FLOW
20 HAS AN EXCESS FLOW
21 HAS AN EXCESS FLOW
_2 2 H A S_ A N E X C E S S_£L QW_
OF 0.110001E 02 CFS
OF 0.3190COE 03 CFS
23 HAS AN EXCESS FLOW . OF 0.199981E 01 CFS
. 2 -V _H AS_ Aji . E X C E S S__F L OW OF 0.20 0019E_0 l_C£S_
25 HAS AN EXCESS FLOW OF 0.299972E 01 CFS
_2_6_HAS AN EXCESS FLOW OF 0.600020E 01 CFS
27 HAS AN EXCESS FLOW
..2_8_HAS_AJ\LEXCES5_ELOiL
29 HAS AN EXCESS FLOW
30 HAS AN EXCESS FLOW
OF 0.999983E 01 CFS
OF 0.300010E 01 CFS
OF 0.699999E 02 CFS
OF 0.250000E 02 CFS
01
-j
10
9
e
-------
2D
00
40
0
e
30
10
10
oo
130
?)
o
| 0
130 120
110
100
90
80
70
60
50
e
o
n
O-JUy 30, 664
O-AUG K). 1964
• -AUC H. I9M
OWPC DATA
ALL LWS
FLDW= JOOO CFS
TEMP -=26 -ZT C
19
I
Q
O
a
e
s r^
j_
110- CO 90 80
DISTARCE-MILES.
70
60
50
i
^j
NJ
FIGURE 9: Observed vs Computed (Dashed Line) Nutrient Profiles, August 1964
-------
20
o>
6 1.0
'f
00
5.0
J30
30
to
i
1.0'
130
• 8
El 0
i-
Q
•
120
no
100
90
80
70
60
50
I
8
T
o
9
•
•
O
5
3
f
O
o
o
9
o
a
•
o
120
110
100
90 80
DISTANCE-MLES
70
60
50
I
~j
ui
FIGURE 10: Observed vs Computed (Dashed Line) Nutrient Profiles, August 1964
-------
-74-
Th e following output pages are a second method of computing deficit profiles
for the Delaware Estuary, where the computation of dissolved oxygen deficit
is accomplished by defining a deficit "species" as part of the system.
In this case it is defined as the NS component. Since ammonia and nitrite
(N2 and N3 correspondingly) are the species producing deficit, a reaction
scheme for these two species is defined. The reaction rates for each of
these is the product of the stochiometric coefficient for the oxydation
reaction times the reaction rate constant of the species that consumes
oxygen.
It should be noted that this method computes total deficit and not the con-
tribution of each component. For a large system, this technique significantly
increases the amount of computer time and main memory required, since it in-
creases the order of the various matrices by the number of segments in the
system.
The pages that follow begin with the "System Reaction Rate Constants", since
the preceding pages would be identical to the example shown previously.
-------
-83-
Execution at the WCC facility
The two versions of the program have been stored as core-image
modules on a partitioned dataset at the EPA Washington Computer
Center Facility (COMNET). Users of this facility may execute the
programs by using the JCL stream listed below.
//EPAIII JOB ('AAAA1,MIII,,'FEDBAK03',B,,20),TIME=(5,30)
//JOBLIB DD UNIT=3330-1,VOL=SER=USER66,DISP=OLD/
// DSNAME=CN.EPAGAN.R2TW.WATER.MODELS
//SI EXEC PGM=PPPPPPPP,REGION=300K
//GO.FT01F001 DD UNIT=SYSDA,DSNAME=+FILE1,DISP=(NEW DELETE),
// SPACE=(720,(720,5)),DCB=(RECFM=FS,BLKSIZE=720,LRECL=720)
//GO.FT03F001 DD UNIT=SYSDA,DSNAME=+FILE3,OISP*(NEW,DELETE),
// SPACE=(1280,(400,5)),DCB=(RECFM=FS,BLKSIZE=1280,LRECL=1280)
//GO.FT06F001 DD SYSOUT=A
//GO.FT07F001 DD SYSOUT=B
//GO.FT09F001 DD SYSOUT=A,DCB=BLKSIZE=120
//GO.FT08F001 DD *
** DATA CARDS GO HERE **
/*
WHERE III = USER'S INITIALS
AAAA = ACCOUNT NUMBER
PPPPPPPP=FEDBAK3A FOR 120 SEGMENT MODEL
PPPPPPPP=FEDBAK3B FOR 180 SEGMENT MODEL
-------
-84-
ACKNOWLEDGEMENTS
Some material presented in this report is based on extensions of the HAR03(4)
computer model, originally documented by Steve Chapra.
The typing was very professionally handled by Ms. Gloria Kasporwitz.
-------
-85-
REFERENCFS FOR THE COMPUTER PROGRAM
(1) Committee on Sanitary Engineering Research, Solubility of Atmos-
pheric Oxygen in Water, J. Sanit. Eng. Div., Am. Soc. Civil Engrs.,
SA 4, 41 (1960)
(2) Nitrification of the Delaware Estuary, prepared for the Delaware
River Basin Commission by Hydroscience, Inc., Westwood, N.J., June
1969.
(3) Thomann, R.V., Systems Analysis and Water Quality Management,
Environmental Science Division, New York, 1971.
(4) Chapra, S.C. and Nossa G.A. HAF.g3_ A.Computer Program for the Modelling
of Water Quality Parameters in .a. Steady State Multi dimensional
Aquatic System., USEPA Region II, New York, N.Y., October 1974.
-------
-86-
APPENDIX A
LISTING OF INPUT DATA FOR VERIFICATION OF DELAWARE ESTUARY