EPA-600/2-77-179C
August 1977
Environmental Protection Technology Series
PREDICTION OF MINERAL
IRRIGATION RETURN FLOW
Volume III. Simulation Model of
Conjunctive Use and Water
Quality for a River System or Basin
Robert S. Kerr Environmental Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Ada, Oklahoma 74820
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ENVIRONMENTAL PROTECTION TECH-
NOLOGY series. This series describes research performed to develop and dem-
onstrate instrumentation, equipment, and methodology to repair or prevent en-
vironmental degradation from point and non-point sources of pollution. This work
provides the new or improved technology required for the control and treatment
of pollution sources to meet environmental quality standards.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600-2-77-179C
August 1977
PREDICTION OF MINERAL QUALITY
OF IRRIGATION RETURN FLOW
VOLUME III
SIMULATION MODEL OF CONJUNCTIVE USE
AND WATER QUALITY FOR A
RIVER SYSTEM OR BASIN
USER'S MANUAL
by
Bureau of Reclamation
Engineering and Research Center
Denver, Colorado 80225
EPA-IAG-D4-0371
Project Officer
Arthur G. Hornsby
Source Management Branch
Robert S. Kerr Environmental Research Laboratory
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OKLAHOMA 74820
-------
DISCLAIMER
This report has been reviewed by the Robert S. Kerr
Environmental Research Laboratory, U.S. Environmental Protection
Agency, and approved for publication. Approval does not signify
that the contents necessarily reflect the views and policies of
the U.S. Environmental Protection Agency, nor does the mention of
trade names or commercial products constitute endorsement or
recommendation for use.
11
-------
FOREWORD
The Environmental Protection Agency was established to coordinate
administration of the major Federal programs designed to protect the
quality of our environment.
An important part of the Agency's effort involves the search for
information about environmental problems, management techniques and
new technologies through which optimum use of the Nation's land and
water resources can be assured and the threat pollution poses to the
welfare of the American people can be minimized.
EPA1s Office of Research and Development conducts this search
through a nationwide network of research facilities.
As one of these facilities, the Robert S- Kerr Environmental
Research Laboratory is responsible for the management of programs to:
(a) investigate the nature, transport, fate and management of pollutants
in groundwater; (b) develop and demonstrate methods for treating waste-
waters with soil and other natural systems; (c) develop and demonstrate
pollution control technologies for irrigation return flows; (d) develop
and demonstrate pollution control technologies for animal production
wastes; (e) develop and demonstrate technologies to prevent, control
or abate pollution from the petroleum refining and petrochemical
industries; and (f) develop and demonstrate technologies to manage
pollution resulting from combinations of industrial wastewaters or
industrial/municipal wastewaters.
This report contributes to the knowledge essential if the EPA is
to meet the requirements of environmental laws that it establish and
enforce pollution control standards which are reasonable, cost effective
and provide adequate protection for the American public.
William C. Galegar
Director
Robert S. Kerr Environmental
Research Laboratory
iii
-------
PREFACE
This report is one of a set which documents the development
and verification of a digital computer modeling effort to predict
the mineral quality changes in return flows occurring as a result
of irrigating agricultural lands. The set consists of five separate
volumes under one general title as follows:
"Prediction of Mineral Quality of Irrigation Return Flow"
Volume I. Summary Report and Verification
Volume II. Vernal Field Study
Volume III. Simulation Model of Conjunctive Use and Water
Quality for a River Basin System
Volume IV. Data Analysis Utility Programs
Volume V. Detailed Return Flow Salinity and Nutrient
Simulation Model
This set of reports represents the culmination of an effort
started in May 1969 by an interagency agreement between the U.S.
Bureau of Reclamation and the Federal Water Pollution Control
Administration on a joint research proposal on the "Prediction
of Mineral Quality of Return Flow Water from Irrigated Land."
This research project has had three different project identifica-
tion numbers during the project period. These numbers (13030 EII,
EPA-IAG-048-(D), and EPA-IAG-D4-0371) are given to avoid confusion
on the part of individuals who have previously tried to acquire
project reports for the earlier project numbers.
IV
-------
ABSTRACT
This volume of the report documents the development of a digital
computer coded simulation.model to predict the effect of irrigation
of agricultural lands on the resulting irrigation return flow
quality. The model is capable of simulating conjunctive uses of
water, however, validation for this purpose was not performed. The
model developed in this volume is much less rigorous than that
presented in Volume V, however, it can be used to provide an assess-
ment of water quality trends due to irrigation at much less cost
than the detailed model. A user's manual is included in the report.
This report was submitted in fulfillment of EPA-IAG-D4-0371 by
the Bureau of Reclamation, Engineering and Research Center, under
the sponsorship of the Environmental Protection Agency.
v
-------
CONTENTS
Introduction l
The Conjunctive Use Simulation Model 4
Main or Executive Program Titled ACUMEN 12
Subroutine PERUSE 33
Subroutine WRIT 69
Subroutine RESOPR 70
Function RESCAP 91
Function RESEL 94
Subroutine PWROPR 95
Subroutine DTABLE 110
Function NFIND 114
Subroutine RFCOMP 115
Subroutine CHEMGR 128
Subroutine CONVER 133
Subroutine VERNAL 137
Function ROOTWO 139
Function ROOT 141
Function FOFX 147
Subroutine SOILCO 148
Subroutine SETUP 163
Subroutine FRSTIM 171
Subroutine GYPEX 176
Subroutine EXCANA 180
Subroutine EXCAMG 186
Subroutine CALCAR 188
Subroutine SCRIBE 193
Function GAMMA 194
Mathematical Appendix of Soil Column Simulation 195
Computer Program Listing for Simulation Model 220
Executive Program ACUMEN 220
Subroutine PERUSE 228
Subroutine WRIT 237
Subroutine RESOPR 242
Function RESCAP 248
Function RESEL 249
Subroutine PWROPR 250
Subroutine DTABLE 254
Function NFIND 256
Subroutine RFCOMP 257
Subroutine CHEMGR 259
Subroutine CONVER 261
Subroutine VERNAL 262
vii
-------
CONTENTS - CONTINUED
Function ROOTWO 264
Function ROOT
Function FOFX
Subroutine SOILCO 271
Subroutine SETUP
Subroutine FRSTIM 274
Subroutine GYPEX 275
Subroutine EXCANA 278
Subroutine EXCAMG 28°
Subroutine CALCAR 281
Subroutine SCRIBE 283
Function GAMMA
Vlll
-------
FIGURES
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
NODE FEATURES
EXAMPLE RIVER BASIN
SYMBOLIC NODE STRUCTURE OF RIVER BASIN
EXAMPLE - NODE WITH SEVEN OPERATION SEQUENCES
GENERAL FLOW OF MODEL
ITERATIVE SCHEME - GENERAL
ITERATIVE SCHEME - DEMANDS
ITERATIVE SCHEME - SUBSURFACE DEMAND
ITERATIVE SCHEME - SURFACE RESERVOIR
ITERATIVE SCHEME - INFLOWS
ITERATIVE SCHEME - END OF LOOP LOGIC
ITERATIVE SCHEME - END OF LOOP LOGIC (CONT.)
FLOW CHART FOR SUBROUTINE PERUSE
EXAMPLE OPERATION SEQUENCE - NODE 2
EXAMPLE - RIVER BASIN WITH SIX NODES
STACKING OF DATA DECK FOR JREAD = 1-3
STACKING OF DATA DECK FOR JREAD = 4-5
AN EXAMPLE RULE CURVE
EXAMPLE CAPACITY - ELEVATION CURVE FOR SURFACE
RESERVOIR
FLOW CHART FOR SUBROUTINE RESOPR
FLOW CHART FOR SUBROUTINE RESOPR (CONT.)
FLOW CHART FOR SUBROUTINE RESOPR (CONT.)
SAMPLE PREDICTED CHARACTERISTIC CURVES
AN EXAMPLE POWER DEMAND CURVE
FLOW CHART FOR SUBROUTINE PWROPR
FLOW CHART FOR SUBROUTINE PWROPR (CONT.)
FLOW CHART FOR SUBROUTINE RFCOMP
FLOW CHART FOR SUBROUTINE RFCOMP (CONT.)
FLOW CHART FOR SUBROUTINE CHEMGR
FLOW CHART FOR SUBROUTINE CHEMGR (CONT.)
FLOW CHART FOR SUBROUTINE VERNAL
FLOW CHART FOR FUNCTION ROOT
LOGIC OF PISTON EFFECT OF SOLUTION MOVING
VERTICALLY THROUGH A SOIL COLUMN (THREE SEGMENTS)
SYMBOLIC NAMES USED IN SIMULATION OF SOIL COLUMN
FLOW CHART FOR SUBROUTINE SOILCO
FLOW CHART FOR SUBROUTINE SOILCO (CONT.)
FLOW CHART FOR SUBROUTINE SOILCO (CONT.)
FLOW CHART FOR SUBROUTINE SOILCO (CONT.)
FLOW CHART FOR SUBROUTINE SETUP
FLOW CHART FOR SUBROUTINE FRSTIM
FLOW CHART FOR SUBROUTINE GYPEX
FLOW CHART FOR SUBROUTINE EXCANA
FLOW CHART FOR SUBROUTINE CALCAR
7
8
9
10
11
16
17
18
19
20
21
22
34
59
61
63
66
73
78
80
81
82
96
99
100
101
123
124
130
131
138
144
154
155
156
157
158
159
165
172
177
182
189
IX
-------
No.
TABLES
Page
1 Variable Names or Data Arrays for Main Program 23
2 Model Input Data 35
3 Variable Names or Data Arrays for Subroutine RESOPR 83
4 Polynomials Used for Steinaker Reservoir 92
5 Variable Names or Data Arrays for Subroutine PWROPR 102
6 Variable Names or Data Arrays for Subroutine DTABLE 112
7 Variable Names or Data Arrays for Subroutine RFCOMP 125
8 Variable Names or Arrays for Subroutine CHEMGR 132
9 Variable Names or Arrays for Subroutine CONVER 135
10 Variable Names or Arrays for Function ROOTWO 140
11 Variable Names or Data Arrays for Function ROOT 145
12 Variable Names or Arrays for Subroutine SOILCO 160
13 Variable Names or Data Arrays for Subroutine SETUP 155
14 Variable Names or Data Arrays for Subroutine FRSTIM 173
15 Variable Names or Data Arrays for Subroutine GYPEX 173
16 Variable Names or Data Arrays for Subroutine EXCANA 133
17 Variable Names or Data Arrays for Subroutine EXCAMG 137
18 Variable Names or Data Arrays for Subroutine CALCAR 190
-------
INTRODUCTION
The majority of present day analyses or studies of river basin water
resources are made for the purpose of increasing certain yields of
an existing physical system or to determine the critical attributes
of an undesirable salinity profile.
The size and complexity of nearly all river basin studies prohibit
the use of manual efforts in obtaining solutions to the multitude of
variations of the physical system operation. Consequently, such
river basin studies are mathematically modeled as electronic computer
applications.
A number of mathematical models are available for use in studying
the physical systems of a river basin. Linear and nonlinear pro-
graming techniques are the basis for some of these models and yet
others are simulation models with varying degrees of complexity and
flexibility.
The availability of so many models makes it quite apparent that a
generalized model that would encompass all the objectives in studying
a river basin is an impossible task.
In the early stages of the conceptual design of this mathematical
model it was realized that limitations due to the varying objectives,
-------
model flexibility, and electronic computer capabilities were going
to be major factors in setting the design framework for the model.
A primary concern in the design of this model was the provision of
maximum flexibility. An examination of other models indicated that
most analysts and managers were interested in obtaining physically
feasible solutions. A physically feasible solution is defined as
the operation of a river basin physical system with a given set of
criteria and/or constraints such that the desired objective could be
accomplished if the system is simulated with this defined format for
"X" number of years in the future. The flexibility then becomes
important because several objectives might be included in the study
of a river basin, and certainly no two river basins are identical.
This examination also revealed a great interest in not only studying
the hydrology but also the changes in salinity within a given river
basin. These changes are often referred to as "salt loading."
Subsequent to the collection of model requirements and the perspective
gained in the aforementioned examination, it was concluded that an
overall model should be designed as two separate submodels;
namely, a data analysis submodel and a simulation submodel.
-------
The data analysis submodel as well as the simulation submodel are
complete and operational at this time.
In this volume the simulation submodel will be discussed in detail
and another volume will contain the description of the data analysis
submodel.
-------
THE CONJUNCTIVE USE SIMULATION MODEL
The conjunctive use capabilities of the model are those required to
simulate simultaneously the use of water resources within the river
basin from both surface resources and subsurface or ground-water
resources.
These capabilities include the resource magnitude as well as its
chemical quality. The chemical quality of the water resource is
further characterized in terms of eight inorganic constituents as
well as the total salts.
The overall simulation model has as a basis a nodal concept or
structure which facilitates the mathematical representation of a
river basin and a simple compact manner of performing calculations,
many of which are iterative in design.
A river basin can be studied as one node or as many as 20 nodes.
The model is designed for a maximum of 20 nodes; however, this
maximum is determined by the limitations of the computer system
being used. The number of nodes used in a particular river basin
study will be a decision the analyst must make on the basis of data
available, the number of response points desired, and the physical
features within the river basin. The node then as a common denominator
can be used to represent the simplest river basin study to one that
is quite complex.
-------
The node can include the simulation of one or all of the following
features:
a. Ten tributary inflows
b. Ten demands of water resources, both surface and subsurface
c. The operation of a surface reservoir
d. The operation of a power facility
e. The operation of a subsurface reservoir (aquifer)
f. The percolation of surface waters vertically through a
soil column
g. The operation of a pumping facility
h. The determination of return flows, both magnitude and
quality, when given consumptive use and conveyance losses.
A pictorial representation of the node features is shown on Figure 1 .
An example of a river basin with three nodes is shown on Figure 2.
The symbolic box representation of the three nodes and the symbolic
representation of the operation sequence of one node are shown on
Figures 3 and 4 , respectively.
The electronic computer application of the conjunctive use simulation
model consists of 24 subroutines or functions plus the executive or
main program. The FORTRAN listings included as part of this volume
for the main program as well as the subroutines and functions are
filled with comments at pertinent points. The extensive use of
-------
these comments is meant to aid in describing the flow of the model and
to provide information within the listings that would be helpful in
making program modifications or conversions to other computer
systems as either become necessary.
The general flow of the model as controlled by the main or executive
program (ACUMEN) is shown on Figure 5. Many different formats have
been used in the presentation of a user's manual. The format used
in this volume will contain one-page flow charts with the aid of
extensive narrative which should provide ease of reading and yet be
as comprehensive as possible.
A list of variable names and arrays and usages thereof will be
inserted at the end of each program, subroutine, or function descrip-
tion. A majority of variable names and data arrays used in the model
will be contained within several labeled common blocks.
-------
NODE FEATURES
SURFACE
RESERVOIR
SURFACE EVAPORATION
J BANK STORAGE AND/OR SEEPAGE
^ RECREATION AND/OR FLOOD CONTROL
I SEDIMENT INFLOW-
CONSUMPTIVE USE
IRRIGATION
POWER
PLANT
CONSUMPTIVE USE
MUNICIPAL
REQUIREMENT
SOIL
COLUMN o
10 SEGMENTS u!
WITH PISTON
EFFECT AND 1
LINEAR MIXING?
INDUSTRIAL
^
REQUIREMENT
FIGU RE I
-------
EXAMPLE RIVER BASIN
-------
SYMBOLIC NODE STRUCTURE OF RIVER BASIN
NODE 3
FIGURE 3
9
-------
EXAMPLE
NODE WITH SEVEN OPERATION SEQUENCES
+ 1
= Inflow
= Demands
FIGURE 4
10
-------
GENERAL FLOW OF MODEL
READ PERTINENT CONTROL
PARAMETERS AND
COEFFICIENTS FOR
MATHEMATICAL EQUATIONS
ASSIGN SEQUENCE OF
OPERATION AND NODE
STRUCTURE WITHIN
RIVER BASIN
READ DATA FOR ONE TIME
FRAME AND PERFORM
CONVERSIONS OF INPUT
UNITS TO COMPUTATION
UNITS
UPDATE PARAMETERS
AND PERTINENT ARRAYS
FOR NEXT TIME
FRAME
USE ITERATIVE
SCHEME TO BALANCE
EACH AND ALL NODES
IN RIVER BASIN
OUTPUT RESULTS FOR
THIS TIME FRAME
LAST TIME
FRAME IN
STUDY ?
-------
MAIN OR EXECUTIVE PROGRAM TITLED ACUMEN
The function of this program is to monitor the flow of the model and
to call various other subroutines or functions as a particular river
basin study dictates. The time frame used throughout the model is
1 month. Time frames other than 1 month can be used in the model
without modification to the concept. When historic data are avail-
able on a daily, weekly, or biweekly time frame, the rule curves
and power demand curves must be changed for each time frame with the
use of transaction cards (see Table 2) as part of the input data set.
The computational mode is essentially in terms of rate; i.e., cubic
feet per second. The conversion of volume (acre-feet) to rate
(cubic feet per second) requires the use of 1 acre-foot =
43,560 cubic feet with the proper transform for the time interval.
The period of study can be unlimited in terms of model capability
because data for one time frame only is stored in computer memory
during any one complete computational cycle.
The model can be used for nearly any river basin study without
modification to the main program or the functions and all the sub-
routines with the exception of the two subroutines which contain the
mathematical formulation of river basin operational criteria and/or
constraints and the subroutine which controls the report generation
or output requirements.
12
-------
The model will accept input data on punched cards or data on magnetic
tape. The input unit is assigned a variable name in the main program
and can be changed to fit the physical form of the input data. The
model requires 42,000 words of rapid access memory.
An earlier version of the simulation model contained a subroutine
titled "CONTRL" which has been discarded in favor of a set of control
cards read in as input and the related logic placed in the main
program. The control cards and logic allow the manipulation of
the data set (as did the subroutine CONTRL) without the usual
redistribution of the data set to accommodate an alternate study
whereby certain physical features of a river basin are omitted or
added as required.
The complete comprehension and appreciation of this capability can
only be acquired through the application of the model to more than
one river basin study.
General flow charts, as shown on Figures 6 to 12, indicate how
demands, inflows, subsurface to surface transfers, surface to sub-
surface transfers, subsurface to subsurface transfers (internodal),
surface reservoirs, and subsurface reservoirs (aquifers) are treated
in the iterative scheme.
13
-------
The iterative scheme begins with the node farthest upstream in the
river basin and completes the hydrologic balance along faith the chemistry
mixes in that node and continues downstream to the next node and
completes the hydrologic balance and chemistry mixes using both
nodes and in this manner proceeds to the most downstream node,
and therefore the whole river basin is balanced.
Whenever a node cannot be balanced hydrologically or a shortage
exceeds the preset shortage index, an appropriate error message will
be printed and computation will be terminated.
Each time a surface reservoir in an upstream node makes a release to
meet a deficit in a downstream demand, the iterative scheme begins
again at that point and updates the flow and chemistry of the river,
which updates the chemistry of all downstream demands as they are
affected by the upstream reservoir release. The reiteration through
the basin from point of incremental surface supply to the point of
deficit is made to update all chemistry with respect to constraining
salinity levels for demands and for selective withdrawals from both
surface and subsurface storage facilities.
The simulation of waters percolating vertically through a soil column
is the last operation completed in each node balance. The requirement
of the soil column simulation is determined by inspection of the
14
-------
volume of water stored in TTANK which acts symbolically as a ponding
pool on the surface for waters to be percolated.
The use of surface to subsurtace and subsurface to surface transfers
is an excellent mechanism for a "black box" approach to return flow
lag times and excess surface flows that cannot be expressed explicitly
because of inadequate data.
The report generating routines are called from the main program as
well as the routines containing the operational criteria. As
mentioned earlier in this volume, these routines must be designed
anew for each specific river basin to be studied.
The variable names or data arrays used in the main program are shown
in Table 1.
15
-------
ITERATIVE SCHEME-GENERAL
IS THIS A
NO / SURFACE
RESERVOIR ?
INITIALIZE
LOOP
LIMITS
SET LOOP
LIMITS FOR
NODES IN BASIN
SET LOOP LIMITS
FOR OPERATION
SEQUENCES IN
NODE
IS THIS A
DEMAND ?
YES
THIS AN
INFLOW
FIGURE 6
16
-------
ITERATIVE SCHEME-DEMANDS
SUBSURFACE
YES
CALL SUBROUTINE
DTABLE AND RANK,
IF ANY, UPSTREAM
RESERVOIRS THAT
MAY BE USED TO
SATISFY RIVER
DEFICIT
ARE THERE
UPSTREAM
RESERVOIRS TO
MEET DEFICIT
IS THIS A
SURFACE OR
SUBSURFACE
DEMAND
SURFACE
CAN FLOW IN
RIVER SATISFY
DEMAND
IS DEMAND
/FOR IRRIGATION
OR MUNICIPAL
,AND INDUSTRI
YES
RESET LOOP
LIMITS WITH NODE
AND SEQUENCE
OF THIS RESERVOIR
ROUTE RELEASE
DOWNSTREAM
THRU SYSTEM
SET CHEMISTRY
OF DEMAND =
CHEMISTRY OF
RIVER
SUBTRACT
DEMAND FROM
RIVER FLOW-IF
NEGATIVE SET
RIVER FLOW
AND CHEMISTRY
TO ZERO SET
DEMAND= RIVER
FLOW
CALL
SUBROUTINE
RFCOMP FOR
RETURN FLOW &
SHORTAGE COMP.
FIGURE 7
17
-------
ITERATIVE SCHEME-SUBSURFACE DEMAND
CAN AQUIFER SUPPLY
THIS DEMAND
CAN AN INTER-
NODAL TRANSFER
BE MADE FROM
AQUIFER IN
NODE ADJACENT
AND UPSTREAM
OF THIS NODE
TO MEET DEFICIT
IN DEMAND
PRINT ERROR
MESSAGE
MAKE TRANSFER
AND UPDATE THE
CAPACITY OF THE
TWO AQUIFERS.
UPDATE CHEMISTRY
OF AQUIFER IN
DEFICIT
HALT COMPUTATION
SET CHEMISTRY
OF DEMAND =
CHEMISTRY OF
AQUIFER AT
START OF TIME
FRAME
SUBTRACT DEMAND
FROM CAPICITY
OF AQUIFER
END
LOOP
FIGURE 8
18
-------
ITERATIVE SCHEME-SURFACE RESERVOHR
CALL SUBROUTINE RES OPR THE
FLOW CHART FOR THIS SUB-
ROUTINE IS INCLUDED IN A
LATER PORTION OF THIS VOLUME
END
LOOP
FltSURE 9
19
-------
TERATIVE SCHEME-INFLOWS
SUBSURFACE
IS THIS A
SURFACE OR
SUBSURFACE
INFLOW
SURFACE
MIX CHEMISTRY OF INFLOW
WITH CHEMISTRY OF
AQUIFER-UPDATE VOLUME
OF AQUIFER
MIX CHEMISTRY OF INFLOW
WITH CHEMISTRY OF
RIVER UPDATE VOLUME OF
FLOW IN RIVER
END
LOOP
FIGURE 10
20
-------
ITERATIVE SCHEME-END OF LOOP LOGIC
END
LOOP
LOOP INDEX
IS INCREMENTED
BY ONE
YES
IS SOIL COLUMN
TO BE SIMULATED
IS LOOP INDEX GREATER
THAN OR EQUAL TO
NUMBER OF SEQUENCES
IN NODE
YES
SIMULATE SOIL
COLUMN
INCREMENT NODE
COUNTER BY ONE
IS NODE LOOP
INDEX GREATER
THAN OR EQUAL
TO NUMBER OF
NODES
YES
SYSTEM IS IN
BALANCE FOR
THIS TIME FRAME
OUTPUT RESULTS
OF THIS TIME
FRAME
FIGURE II
21
-------
ITERATIVE SCHEME-END OF LOOP LOGIC (CONT.)
NO
INCREMENT TIME
FRAME BY ONE
IS TIME FRAME INDEX
I.E. MONTH AND YEAR
EQUAL TO ENDING
TIME FRAME I.E. MONTH
AND YEAR
YES
HALT
COMPUTATION
FIGURE 12
22
-------
Table 1
VARIABLE NAMES OR DATA ARRAYS FOR MAIN PROGRAM
Variable
name
Description
Remarks
or Units
BEGWR
CHEMBT
CHEMDM
CHEMGR
CHEMRV
CHEMTT
COMPA
COMP
CONST
CONVER
DELTAF
DEMAND
DTABLE
GRTANK
Data array for beginning of time
frame conditions for aquifer
Entry to subroutine CHEMGR
Data array for chemistry of
demands
Entry to subroutine CHEMGR
Entry to subroutine CHEMGR
Entry to subroutine CHEMGR
Computation aid
Computation aid
Computation aid
Entry to subroutine CONVER
Reservoir release made to meet a
downstream deficit
Data array for demands
Entry to subroutine DTABLE
Data array for volume and
chemistry of aquifer
See data
array BEGWR
See subroutine
CHEMGR
See data
array CHEMDM
See subroutine
CHEMGR
See subroutine
CHEMGR
See subroutine
CHEMGR
See subroutine
CONVER
Acre-feet
See data array
DEMAND
See subroutine
DTABLE
See data array
GRTANK
C* = Variable names shared in common blocks
23
-------
Table 1 - Continued
Variable
name
C*
Description
Remarks
or Units
1C
IFIR
INDY
IR
ISEQ
IS
IW
IYRBEG
IYREND
IYR
JREAD
KA
KFAK
KG
KNODE
KN
Index of number of upstream
reservoirs
Index used to run jobs back to
back redundant in this version
Index used in calling argument
for subroutine RFCOMP
Read unit designation
Data array for operation
sequences in node
Absolute value for index of
sequence operation
Write unit designation
Beginning year of study
Ending year of study
Year of study during computa-
tional cycle
Read index for subroutine PERUSE
Loop index
Data array to store ranking of
reservoirs
Computational aid
Node index
Loop index
See subroutine
RFCOMP
2 digits
i.e., 69,70
2 digits
2 digits
See subroutine
PERUSE
See subroutine
DTABLE
24
-------
Table 1 - Continued
Variable
name
C*
Description
Remarks
or Units
KPRI
K
KS
KUSE
KX
KZ
LINE
MINA
MIN
MOBEG
MODE
MOEND
MO
NCTRL
ND
NFIND
NPASS
Redundant in this version
Loop index
Computational aid
Data array for consumptive use
(equivalence to CONUSE)
Computational aid
Computational aid
Line counter - redundant
Beginning index for operational
sequence loop in iterative
scheme
Beginning index for node sequence
loop in iterative scheme
Beginning month of study
Computational aid
Ending month of study
Month or time frame of study
during computation
Data array for the control of
flow from node to node
Computational aid
Entry to function NFIND
Computational aid which counts
passes thru each reservoir
during one time frame
See data array
CONUSE
2 digits
i.e., JAN=01
2 digits
i.e., JUN=06
2 digits
i.e., DEC=12
See function
NFIND
25
-------
Table 1 - Continued
Variable
name
C*
Description
Remarks
or Units
NSEG
NSOIL
NS
NUMSEQ
OBSERV
OBSVAL
PERUSE
PREDIC
PRELM
QFLOW
QINFLO
QINRIV
RDATA
RESOPR
Number of segments in any one soil
column
Sequence of any soil column in
river basin (one per node)
Number of sequences of operation
in any one node used as loop
index in iterative scheme
Number of sequences of operation
in any one node
Total of observed flows and
chemistry as checks in any one
node
Same as OBSERV
Entry to read subroutine PERUSE
Predicted flows and chemistry
leaving each node
Computational aid to determine
possible reservoir release
An array used to transfer data as
a part of argument
An array used to store inflows to See data array
system QINFLO
An array used to store river flow
and chemistry at any point in
river basin
An array used to store perti-
nent reservoir data
Entry to reservoir operation
subroutine RESOPR
See subroutine
RESOPR
See subroutine
RESOPR
26
-------
Table 1 - Continued
Variable
name
ov?
Description
Remarks
or Units
RFCOMP
SOILCO
TEMPR
TEMP
TEST
TR
TTANK
VERNAL
WRIT
Entry to return flow subroutine
RFCOMP
Entry to monitor subroutine for
soil column simulation
An array used as a computational
aid
A computational array
A computational aid
Computational aid
An array used to store flows and
chemistry for water to be perco-
lated thru soil column
Subroutine containing the opera-
tion criteria for the VERNAL
RIVER BASIN
Subroutine used to output
results or generate report
See subroutine
RFCOMP
See subroutine
SOILCO
See data array
TTANK
This subroutine
will be re-
written for
each river
basin
This subroutine
will be re-
written for
each river
basin
27
-------
Table I - Continued
DATA ARRAYS IN MAIN PROGRAM
DATA ARRAY BEGWR (20,10)
BEGWR (J,l) = Volume of water stored in subsurface reservoir
(aquifer) at beginning of time frame
BEGWR (J,2) = Concentration of calcium in aquifer in ppm
BEGWR (J,3) = Concentration of magnesium
BEGWR (J,4) = Concentration of sodium
BEGWR (J,5) = Concentration of chloride
BEGWR (J,6) = Concentration of sulfate
BEGWR (J,7) = Concentration of bicarbonate
BEGWR (J,8) = Concentration of carbonate
BEGWR (J,9) = Concentration of nitrate
BEGWR (J,10) = Concentration of total salts (TDS)
J = Node sequence
DATA ARRAY CHEMDM (20,10,10)
CHEMDM (J,K,1) = Blank
CHEMDM (J,K,2) = Concentration of calcium in water used to
meet demand in ppm
CHEMDM (J,K,3) = Concentration of magnesium
CHEMDM (J,K,4) = Concentration of sodium
CHEMDM (J,K,5) = Concentration of chloride
CHEMDM (J,K,6) = Concentration of sulfate
CHEMDM (J,K,7) = Concentration of bicarbonate
28
-------
Table 1 - Continued
DATA ARRAY CHEMDM (20,10,10) - Continued
CHEMDM (J,K,8) = Concentration of carbonate
CHEMDM (J,K,9) = Concentration of nitrate
CHEMDM (J,K,10) = Concentration of total salts (TDS)
J = Node sequence K = Sequence of operation in node
DATA ARRAY DEMAND 20,10,16
DEMAND (J,K,1-12) = Monthly demand in acre-feet
DEMAND (J,K,13)
DEMAND (J,K,14)
DEMAND (J,K,15)
DEMAND (J,K,16)
J = Node sequence
= Source of supply
Alpha = SUR index = 1.0 for surface supply
Alpha = GWR index =2.0 for subsurface supply
= Use or transfer designation
Irrigation = IRR =1.0
Municipal and Industrial = M&I = 2.0
Transfer to aquifer = GWR = 3.0
Transfer to soil column (TTANK) = GWV =4.0
Transfer to river = SUR = 5.0
Transfer out of system = OUT =6.0
All units are in acre-feet
= Consumptive use flag
0 = No consumptive use
1 = Consumptive use
= Actual diversion made to satisfy demand in
acre-feet
K = Sequence of operation in node
DATA ARRAY GRTANK (20,10)
GRTANK (J,l)
GRTANK (J,2)
= Volume of water in subsurface storage (aquifer)
in acre-feet
= Concentration of calcium in aquifer in ppm
29
-------
Table 1 - Continued
DATA ARRAY GRTANK (20,10) - Continued
GRTANK (J,3) = Concentration of magnesium
GRTANK (J,4) = Concentration of sodium
GRTANK (J,5) = Concentration of chloride
GRTANK (J,6) = Concentration of sulfate
GRTANK (J,7) = Concentration of bicarbonate
GRTANK (J,8) = Concentration of carbonate
GRTANK (J,9) = Concentration of nitrate
GRTANK (J,10) = Concentration of total salts (TDS)
J = Node index
DATA ARRAY CONUSE (20,10,30) EQUIVALANCE KUSE (20,10,30)
CONUSE (J,K,1-12) = Monthly consumptive use in acre-feet
CONUSE (J,K,13) = Conveyance and farm losses in percent of
demand (see subroutine RFCOMP)
CONUSE (J,K,14) = Shortage index percent
CONUSE (J,K,15) = Shortage in acre-feet
CONUSE (J,K,19-30) = Return flow allocations and destinations
CONUSE (J,K,19,23,27)= Use or destination indicator
SUR = Return to river
GWR = Return to aquifer
GWV = Return to TTANK for percolation thru
soil column
CONUSE (J,K,20,24,28)= Node number of destination
CONUSE (J,K,21,25,29)= Operation sequence of destination
30
-------
Table 1 - Continued
DATA ARRAY CONUSE (20,10,30) EQUIVALANCE KUSE (20,10,30) - Continued
CONUSE (J,K,22,26,30)= Percent of total return flow allotted to
each destination or use
J = Node index
K = Sequence of operation in node
DATA ARRAY QINFLO
QINFLO (J,K,1)
QINFLO (J,K,2)
QINFLO (J,K,3)
QINFLO (J,K,4)
QINFLO (J,K,5)
QINFLO (J,K,6)
QINFLO (J,K,7)
QINFLO (J,K,8)
QINFLO (J,K,9)
QINFLO (J,K,10)
QINFLO (J,K,11)
QINFLO (J,K,12)
J = Node index
(20,10,12)
= Volume of inflow in acre-feet
= Concentration of calcium in ppm
= Concentration of magnesium
= Concentration of sodium
= Concentration of chloride
= Concentration of sulfate
= Concentration of bicarbonate
= Concentration of carbonate
= Concentration of nitrate
= Concentration of total salts (TDS)
= Destination of inflow
1.0 = To river
2.0 = To aquifer
3.0 = To ponding pool above soil column
(TTANK)
4.0 = To check point (observed flow and
chemistry at check point to compare
with predicted flow and chemistry
at check point or response point).
= Units of chemical concentrations
1 = MEQ/liter
2 = Parts per million (ppm)
K = sequence of operation in node
31
-------
Table I - Continued
DATA ARRAY TTANK (20,10)
TTANK (J,l) = Volume of water in ponding pool above surface
to percolate thru soil column in acre-feet
TTANK (J,2) = Concentration of calcium in ppm
TTANK (J,3) = Concentration of magnesium
TTANK (J,4) = Concentration of sodium
TTANK (J,5) = Concentration of chloride
TTANK (J,6) = Concentration of sulfate
TTANK (J,7) = Concentration of bicarbonate
TTANK (J,8) = Concentration of carbonate
TTANK (J,9) = Concentration of nitrate
TTANK (J,10) = Concentration of total salts (TDS)
J = Node index
DATA ARRAY BTANK (20,10)
BTANK (J,l) = Volume of water in tank immediately below soil
column and above aquifer
BTANK (J,2) = Concentration of calcium in ppm
BTANK (J,3) = Concentration of magnesium
BTANK (J,4) = Concentration of sodium
BTANK (J,5) = Concentration of chloride
BTANK (J,6) = Concentration of sulfate
BTANK (J,7) = Concentration of bicarbonate
BTANK (J,8) = Concentration of carbonate
BTANK (J,9) = Concentration of nitrate
BTANK (J,10) = Concentration of total salts (TDS)
J = Node index
32
-------
SUBROUTINE PERUSE
All read instructions for the model are contained within this sub-
routine. Also, most of the data arrays are filled or partially filled
in this subroutine. The read formats were designed in a manner that
will allow the input of several different types of data with the use
of the same format.
A general flow chart of this subroutine is shown on Figure 13.
The index used in the PERUSE subroutine is assigned the variable name
JREAD which is an integer number varying between 1 and 5. All one-
time reads made for any one particular simulation are made when JREAD
varies between 1 and 4.
Each of the card formats used to read the input data are shown in
Table 2. Many river basin studies will not have all of the data or
physical features included in the table of model input data. The
lack of data or physical features for any given river basin study
presents no problem as far as the read subroutine is concerned or
for that matter the flow of the model.
The control cards read when JREAD=1 are probably the most difficult
to explain, and yet once they are fully understood, they are simple
to apply and represent the versatility of the model in general. The
33
-------
FLOW CHART SUBROUTINE PERUSE
JREAD =
i
ENTER WITH ARGUMENT (JREAD)
JREAD =
1
JREAD = 3
JREAD =
JREAD =5
READ TITLE
CARD
i
READ NODE
CONTROL AND
NODE SEQUENCE
CARDS-END WITH
SIGNAL CARD
-"CONEND" IN
FIRST SIX
COLUMNS
RETURN
READ INITIAL
VALUES FOR
SURFACE
RESERVOIRS,
AQUIFERS,
TITLES AND
INITIAL SOIL
ANALYSIS
DATA CARDS
1
CONVERT INPUT
UNITS TO
COMPUTATION
UNITS
END WITH SIGNAL
CARD-"COLEND"
IN FIRST SIX
COLUMNS
I
READ
COEFFICIENTS FOR
MATHEMATICAL
EQUATIONS
DESCRIBING
RESERVOIR AREA-
CAPACITY, TAIL-
WATER, SEDIMENT,
TURBINE HORSE-
POWER, DISCHARGE
AND CONSTRAINTS,
FIRM POWER, AND
RULE CURVES,
EVAPORATION AND
BANK STORAGE
COEFFICIENTS-END
WITH SIGNAL CARD
"EQUEND" IN
FIRST SIX COL-
UMNS
'
READ INFLOW,
DEMAND AND
CONSUMPTIVE USE
DATA FOR FIRST
TIME FRAME WITH
TITLES FOR IN-
FLOWS AND
DEMANDS
1
CONVERTER INPUT
UNITS TO
COMPUTATION
UNITS
1
END WITH SIGNAL
CARD-"ENDATA"
IN FIRST SIX
COLUMNS
READ INFLOW,
DEMAND AND
TRANSACTION
CARDS FOR ALL
TIME FRAMES
OTHER THAN
FIRST TIME
FRAME
,
CONVERT INPUT
UNITS TO
COMPUTATION
UNITS
'
END WITH SIGNAL
CARD-"ENDATA"
IN FIRST SIX
COLUMNS
,
FIGURE 13
-------
Table 2
MODEL INPUT DATA
CARD
JREAD FORMAT COLUMNS
1 8A8 1-64
116 65-80
1 2R3 1-6
13 7-9
13 10-12
13 13-15
1 2R3 1-6
13 7-9
13 10-12
2013 13-72
VARIABLE
NAME
DESCRIPTION OR ARRAY
Title card, name of LTITLE
river basin study
Redundant in this KPRI
version of model
Read sentinel
Node number NCTRL(J,1)
Node number of NCTRL(J,2)
node immediately
downstream of
NCTRL(20,1)
Sequence number in NCTRL(j,3)
operational sequence
of downstream node
Read sentinel
Node number
Number of sequen- NUMSEQ(J)
tial operations in
node
Sequential ISEQ(J,K)
operation
Without sign =
inflows
With minus sign =
demands
Zero = surface
reservoir
IDENTIFIER
None
None
CONNOD
None
None
None
CONSEQ
None
None
None
J = Node sequence (1-20)
K = Operation sequence in node (1-20)
35
-------
Table 2 - Continued
JREAD FORMAT
1 2R3
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
13
CARD
COLUMNS DESCRIPTION
1-6 Read sentinel
7-9 Node number
13-15 Upstream node No.
16-18
19-21
22-24 " " "
25-27
28-30 " " "
31-33
4-36 " " "
37-39
40-42 " "
43-45 " "
46-48 " " "
49-51 " "
52-54
55-57
58-60 " "
61-63 " " "
64-66 " " "
67-69 " " "
VARIABLE
NAME
OR ARRAY IDENTIFIER
CONWAT
NDWTR(J,1)
NDWTR(J,2)
NDWTR(J,3)
NDWTR(J,4)
NDWTR(J,5)
NDWTR(J,6)
NDWTR(J,7)
NDWTR(J,8)
NDWTR(J,9)
NDWTR(J,10)
NDWTR(J,11)
NDWTR(J,12)
NDWTR(J,13)
NDWTR(J,14)
NDWTR(J,15)
NDWTR(J,16)
NDWTR(J,17)
NDWTR(J,18)
NDWTR(J,19)
36
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
13 70-72 Upstream node No. NDWTR(J,20)
J = Node index
2R3 1-6 Read sentinel CONEND
13 7-9 Beginning month of MOBEG 2 digits
study
13 10-12 Beginning year of IYRBEG 2 digits
study
13 13-15 Ending month of MOEND 2 digits
study
13 16-18 Ending year of IYREND 2 digits
study
37
-------
Table 2 - Continued
JREAD FORMAT
2 2R3
13
8A8
CARD
COLUMNS
1-6
10-12
16-79
DESCRIPTION
Read sentinel
Node number
Name of aquifer
VARIABLE
NAME
OR ARRAY IDENTIFIER
AQUFER
KAQUMN(J,l-8)
J = Node index
2 2R3
11
13
F6.0
F6.0
F6.0
F6.0
F6.0
F6.0
F6.0
F6.0
1-6
9
10-12
16-21
22-27
28-33
34-39
40-45
46-51
52-57
58-63
Read sentinel
Units indicator
1= MEQ/liter
2= PPM
Node number
Volume of aquifer
in AF * 10m
Calcium concen-
tration
Magnesium concen-
tration
Sodium concentra-
tion
Chloride concen-
tration
Sulfate concen-
tration
Bicarbonate con-
centration
Carbonate concen-
GRTANK
GRTANK(J,1)
GRTANK(J,2)
GRTANK (J, 3)
GRTANK(J,4)
GRTANK(J,5)
GRTANK(J,6)
GRTANK(J,7)
GRTANK (J, 8)
tration
38
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
F6.0 64-69 Nitrate concen- GRTANK(J,9)
tration
F6.0 70-75 Total salts (IDS) GRTANK(J,10)
12 76-77 Month
II 78 "m" a power of 10
12 79-80 Year
J = Node index
J =
2R3 1-6
13 10-12
8A8 16-79
Node index
Read sentinel
Node number
Name of reservoir KSURNM(J,1-8)
SURRES
2 2R3
12
13
F6.0
F6.0
F6.0
Ffi.n
1-6
7-8
10-12
16-21
22-27
28-33
34-39
Read sentinel
Operation sequence
number in node
Node number
Volume of reser-
voir in AF * 10m
Maximum reservoir
capacity AF * 10m
Minimum reservoir
release.cfs * 10m
Maximum release w/o
RESDAT
RDATA(J,2)
RDATA(J,1)
RDATA(J,11)
RDATA(J,12)
RDATA(J,13)
spill cfs * 10m
39
-------
Table 2 - Continued
VARIABLE
CARD NAME
JREAD FORMAT COLUMNS DESCRIPTION OR ARRAY IDENTIFIER
F6.0 40-45 Maximum release RDATA(J,14)
w/o spill cfs *
10m
F6.0 46-51 Sediment deposited RDATA(J,20)
in reservoir AF
* 10m
II 78 "m" a power of 10
J = Node index
2 2R3 1-6 Read sentinel RESCHM
II 9 Units indicator
1 = MEQ/Liter
2 = PPM
13 10-12 Node number
F6.0 22-27 Calcium concen- RDATA(J,101)
tration
F6.0 28-33 Magnesium concen- RDATA(J,102)
tration
F6.0 34-39 Sodium concentra- RDATA(J,103)
tion
F6.0 40-45 Chloride concen- RDATA(J,104)
tration
F6.0 46-51 Sulfate concen- RDATA(J,105)
tration
F6.0 52-57 Bicarbonate con- RDATA(J,106)
centration
F6.0 58-63 Carbonate concen- RDATA(J,107)
tration
40
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
F6.0 64-69 Nitrate concen- RDATA(J,108)
tration
F6.0 70-75 Total salts (TDS) RDATA(J,109)
J = Node index
J =
2R3 1-6 Read sentinel
12 7-8 Number of segments
in soil column
13 10-12 Node number
Node index
NSEG(J)
F5.0
tration in MEQ/
liter
33-37 Sulfate concen-
tration in MEQ/
liter
SOIL(J,K,5)
COLBEG
2 2R3 1-6
13 7-9
13 10-12
F5.0 13-17
F5.0 18-22
F5.0 23-27
F5.0 28-32
Read sentinel COLSEG
Node number
Segment number
Calcium concen- SOIL(J,K,1)
tration in MEQ/
liter
Magnesium concen- SOIL(J,K,2)
tration in MEQ/
liter
Sodium concentra- SOIL(J,K,3)
tion in MEQ/liter
Chloride concen- SOIL(J,K,4)
41
-------
Table 2 - Continued
VARIABLE
CARD NAME
JREAD FORMAT COLUMNS DESCRIPTION OR ARRAY IDENTIFIER
F500 38-42 Bicarbonate con- SOIL(J,K,6)
centration in MEQ/
liter
F5.0 43-47 Carbonate concen- SOIL(J,K,7)
tration in MEQ/
liter
F5.0 48-52 Nitrate concen- SOIL(J,K,8)
tration in MEQ/
liter
F3.0 53-55 Soil/water ratio SOIL(J,K,9)
of extract
F5.0 56-60 Volume of water in SOIL(J,K,10)
segment in milli-
liters
F5.0 61-65 Cation exchange SOIL(J,K,11)
capacity in MEQ/
100 gr
F5.0 66-70 Concentration of SOIL(J,K,12)
Gypsum Dihydrate in
units of MEQ/100 gr
Fl.O 71 Lime indicator SOIL(J,K,13)
0 = No lime present
1 = Lime present
F4.0 72-75 Bulk density of SOIL(J,K,14)
soil in segment in
units of grams/
cubic centimeter
F5.0 76-80 Length of soil SOIL(J,K,15)
segment in centi-
meters
J = Node index K = Operation sequence in node
42
-------
Table 2 - Continued
VARIABLE
CARD NAM£
JREAD FORMAT COLUMNS DESCRIPTION OR ARRAY IDENTIFIER
2R3 1-6 Read sentinel for COLEND
end JREAD=2
3 2R3 1-6 Read sentinel EQURES
14 7-10 Node number
13 11-13 Index "L" in array
sequence
F12.0 14-25 Base elevation in RDATA(J,K,L)
feet above msl
13 26-28 Index "L" in array
sequence
F12.0 29-40 Coefficient AQ in RDATA(J,K,L)
polynomial
13 41-43 Index "L" in array
sequence
F12.0 44-55 Coefficient AL in RDATA(J,K,L)
polynomial
13 56-58 Index "L" in array
sequence
F12.0 59-70 Coefficient A£ in RDATA(J,K,L)
polynomial
2
Polynomial - Capacity = A0 + Aj_H + A2H - where H = height above
base elevation
J = Node index K = Operation sequence in node
2R3 1-6 Read sentinel EQUSED
14 7-10 Node number
43
-------
Table 2 - Continued
VARIABLE
CARD NAME
JREAD FORMAT COLUMNS DESCRIPTION OR ARRAY IDENTIFIER
13 11-13 Index = 15
F12.0 14-25 Coefficient "A" in RDATA(J,K,15)
first equation
13 26-28 Index = 16
F12.0 29-40 Coefficient "B" in RDATA(J,K,16)
first equation
13 41-43 Index = 17
F12.0 44-55 Coefficient "A" in RDATA(J,K,17)
second equation
13 56-58 Index = 18
F12.0 59-70 Coefficient "B" in RDATA(J,K,18)
second equation
F10.0 71-80 Spline point for RDATA(J,K,19)
equations in cfs
-n
Equation - Sediment = AQ
J = Node index K = Operation sequence in node
2R3 1-6 Read sentinel EQUTWR
14 7-10 Node number
13 11-13 Index = 1
F12.0 14-25 Base elevation for TDATA(J,1)
tailwater curve in
feet above msl
13 26-28 Index = 2
F12.0 29-40 Coefficient AQ in TDATA(J,2)
polynomial
44
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
13 41-43 Index = 3
F12.0 44-55 Coefficient
polynomial
13 56-58 Index = 4
in TDATA(J,3)
F12.0 59-70 Coefficient A2 in TDATA(J,4)
polynomial
Polynomial - Tailwater elevation = Base elevation +
where H = height above base elevation
J = Node index
3 2R3
14
13
F12.0
13
F12.0
1-6
7-10
11-13
14-25
26-28
29-40
Read sentinel EQUEFF
Node number
Index = 1
Coefficient AQ in ENERGY (J-,1)
polynomial
Index = 2
Coefficient AL in ENERGY(J,2)
polynomial
13
41-43 Index = 3
F12.0 44-55 Coefficient A2 in ENERGY(J,3)
polynomial
Polynomial - Horsepower = A0+A]_H+A2H2 where H = Net head
J = Node index
45
-------
Table 2 - Continued
JREAD
3
FORMAT
2R3
14
13
F12.0
13
F12.0
13
F12.0
Polynomial -
3
2R3
14
13
F12.0
13
F12.0
CARD
COLUMNS
1-6
7-10
11-13
14-25
26-28
29-40
41-43
Discharge
where H =
1-6
7-10
11-13
14-25
26-28
29-40
VARIABLE
NAME
DESCRIPTION OR ARRAY
Read sentinel
Node number
Index = 5
Coefficient AQ in ENERGY(J,5)
polynomial
Index = 6
Coefficient A^ in ENERGY(J,6)
polynomial
Index = 7
Coefficient A in ENERGY(J,7)
polynomial
(cfs) at full gate = A0+A1H+A2H2
Net Head. J = Node index
Read sentinel
Node number
Index = 1
Number of genera- CAPBAB(J,1)
tor units in power-
plant
Index = 2
Generator rating CAPBAB(J,2)
IDENTIFIER
EQUEFF
EQUCAP
(KV-A)
13
41-43 Index = 3
46
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
F12.0 44-55
13 56-58
F12.0 59-70
Power factor
(dimensionless)
Index = 4
Percent of over-
load
CAPBAB(J,3)
CAPBAB(J,4)
J = Node index
3
2R3
14
13
1-6
7-10
11-13
Read sentinel
Node number
Index = 6
EQUCAP
F12.0 14-25 Generator effi-
ciency in %
CAPBAB(J,6)
13
26-28 Index = 7
F12.0 29-40 Maximum head in ft CAPBAB(J,7)
13
41-43 Index =
F12.0 44-55 Minimum head in ft CAPBAB(J,8)
J = Node index
Three cards with identical format are used for the firm power
requirements on a monthly time frame for a one-year period.
2R3 1-6 Read sentinel
14 7-10 Node number
13 11-13 Index = 1, 5, 9
EQUPOW
PDATA(J, Index)
47
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
F12.0
13
F12.0
13
FL2.0
13
F12.0
14-25 Firm power require-
ment in KWH for
months JAN, MAY,
SEPT
26-28 Index = 2, 6, 10
29-40 Firm power require- PDATA(J,Index)
ment in KWH for
months FEE, JUNE,
OCT
41-43 Index = 3, 7, 11
44-55 Firm power require- PDATA(J,Index)
ment in KWH for
months MARCH, JULY,
NOV
56-58 Index = 4, 8, 12
59-70 Firm power require-
ment in KWH for
months APRIL, AUG,
DEC
J - Node index
Three cards with identical format are used for the reservoir rule
curves. These curves can either be in terms of reservoir capacity
in acre-feet or reservoir water surface elevation in feet above mean
sea level on a monthly time frame for a one-year period.
2R3 1-6 Read sentinel
14 7-10 Node number
13 11-13 Index = 1, 5, 9
EQURUL
RULE(J,Index)
48
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
F12.0 14-25
13
F12.0
13
F12.0
13
F12.0
26-28
29-40
41-43
44-55
56-58
59-70
*F10.0 71-80
*0n last card only.
J = Node index
Reservoir rule
curve for months
JAN, MAY, SEPT
Index =2, 6, 10
Reservoir rule
curve for months
FEE, JUNE, OCT
Index = 3, 7, 11
Reservoir rule
curve for months
MAR, JULY, NOV
Index = 4, 8, 12
Reservoir rule
curve for months
APRIL, AUG, DEC
Indicator
0 = Capacities
1 = Elevations
RULE(J,Index)
RULE(J,Index)
RULE(J,13)
Three cards with identical format are used for the reservoir surface
evaporation losses in feet per acre of area on a monthly time frame
for a one-year period.
2R3 1-6 Read sentinel
14 7-10 Node number
13 11-13 Index = 1, 5, 9
EQUEVP
49
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
F12.0 14-25 Reservoir evapo-
ration in feet per
acre area for
months JAN, MAY,
SEPT
13 26-28 Index = 2, 6, 10
F12.0 29-40 Reservoir evapo-
ration in feet per
acre area for
months FEE, JUNE,
OCT
13 41-43 Index = 3, 7, 11
F12.0 44-55 Reservoir evapo-
ration in feet per
acre area for
months MAR, JULY,
NOV
13
F12.0
56-58 Index = 4, 8, 12
59-70 Reservoir evapo-
ration in feet per
acre area for
months APRIL, AUG,
DEC
CONEVP(J,Index)
CONEVP(J,Index)
CONEVP(J,Index)
CONEVP(J,Index)
J = Node index
Three cards with identical format are used for the reservoir changes
in bank storage in decimal units of change in capacity per monthly
time frame for a one-year period.
2R3 1-6 Read sentinel
14 7-10 Node number
EQUBNK
50
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
VARIABLE
NAME
OR ARRAY
IDENTIFIER
13 11-13 Index = J
F12.0 14-25 Reservoir bank BANKCO(J)
storage as a
decimal
13 26-28 Index = J
F12.0 29-40 Reservoir bank BANKCO(J)
storage as a
decimal
13 41-43 Index = J
F12.0 44-55 Reservoir bank BANKCO(J)
storage as a
decimal
13 56-58 Index = J
F12.0 59-70 Reservoir bank BANKCO(J)
storage as a
decimal
J = Node index
2R3 1-6 Read sentinel for
end of JREAD = 3
EQUEND
* 4
2R3
12
1-6
7-8
Read sentinel
Operation sequence
SURQIN
number in node
13
10-12 Node number
*This card is read for first time frame only. QINFLO card must
follow this card.
51
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
8A8 16-79 Alphanumeric title KQINNM(J,K,L)
of inflow
J = Node index K = Operation sequence in node
L = 1-1
4 or 2R3
5
12
II
13
R3
1-6 Read sentinel
7-f
10-12
13-15
F6.0 16-21
F6.0 22-27
F6.0 28-33
QINFLO
Operation sequence
number in node
Chemistry units
indicator
1 = MEQ/Liter
2 = PPM
Node number
Alpha destination
indicator
SUR = Surface
source = 1.0
GWR = Subsurface
source = 2.0
GWV = Redundant
in this version
CHK = Inflow
check point =
4.0
Volume of inflow
in AF * 10m
Calcium concen-
tration
Magnesium concen-
tration
QINFLO(J,K,12)
QINFLO(J,K,11)
QINFLO(J,K,1)
QINFLO(J,K,2)
QINFLO(J,K,3)
52
-------
Table 2 - Continued
VARIABLE
CARD NAME
JREAD FORMAT COLUMNS DESCRIPTION OR ARRAY IDENTIFIER
F6.0 34-39 Sodium concen- QINFLO(J,K,4)
tration
F6.0 40-45 Chloride concen- QINFLO(J,K,5)
tration
F6.0 46-51 Sulfate concen- QINFLO(J,K,6)
tration
F6.0 52-57 Bicarbonate con- QINFLO(J,K,7)
centration
F6.0 58-63 Carbonate con- QINFLO(J,K,8)
centration
F6.0 64-69 Nitrate concen- QINFLO(J,K,9)
tration
F6.0 70-75 Total salts in ppm QINFLO(J,K,10)
(TDS)
12 76-77 Month 2 digits
II 78 "m" a power of 10
12 79-80 Year 2 digits
J = Node index K = Operation sequence in node
* 4 2R3 1-6 Read sentinel SURDEM
12 7-8 Operation sequence
number in node
13 10-12 Node number
8A8 16-79 Alphanumeric title KDEMNM(J,K,L)
of demand
J = Node index K = Operation sequence in node L = 1-8
*This card is read for first time only. Demand card (i.e., DEMSUR
or DEMGWR) must follow this card.
53
-------
Table 2 - Continued
JREAD
4 or
5
FORMAT
R3
R3
CARD
COLUMNS
1-3
4-6
VARIABLE
NAME
DESCRIPTION OR ARRAY
Read sentinel
Read sentinel for DEMAND
IDENTIFIER
DEM
SUR =
Sur-
source of demand (J,K,13)
face demand
= 1.0
GWR = Sub-
surface
demand =
2.0
12 7-8 Operation sequence
number in node
II 9 Input indicator
where 1 indicates
consumptive card
to follow and 0
indicates no con-
sumptive use card
to follow
13 10-12 Node number
R3 13-15 Alpha indicator DEMAND
for destination (J,K,14)
of demand
IRR = Irrigation
= 1.0
MAI = Municipal or
industrial = 2.0
GWR = Aquifer =
3,0
GWV = Ponding pool
above soil column
(TTANK) =4.0
SUR = River =5.0
OUT = Demand leaves
river basin = 6.0
F6.0 16-21 Volume of demand DEMAND
in AF * lQm (J,K,L)
54
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
12 76-77 Month 2 digits
II 78 "m" a power of 10
J = Node index K = Operation sequence in node L = tnonth(l-12)
4 or
5
R6
1-6
Read
sentinel
CONUSE
Block
#1
13
12
F9.0
F3.0
F3.0
R3
13
7-9 Node number
10-11 Operation sequence
number in node
12-20 Volume of con-
sumptive use in
acre-feet
21-23 Conveyance and
farm losses in °L
24-26 Shortage index in
percent
27-29 Destination of
return flow.
Block #1 -
SUR = Return to
river
GWR = Return to
aquifer
GWV = Return to
ponding pool
above soil column
(TTANK)
30-32 Node number of
destination of
return
CONUSE
(J,K,L)
CONUSE
(J,K,13)
CONUSE
(J,K,14)
KUSE
(J,K,19)
KUSE
(J,K,20)
55
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY IDENTIFIER
12 33-34 Operation sequence
number in node of
return
Block
#1 13 35-37 Percent of return
J allocated to this
„!, destination
.>
R3 38-40 Second destina-
tion of return
flow (see Block
#1)
13 41-43 Node number of
destination of
Block return
#2
•>
*>
12 44-45 Operation sequence
number in node of
return
13 46-48 Percent of return
allocated to this
„ destination
^ R3 49-51 Third destination
return flow (See
Block #1)
13 52-54 Node number of
destination of
return
Block
#3 12 55-56 Operation sequence
number in node of
return
13 57-59 Percent of return
allocated to this
KUSE
(J,K,21)
KUSE
(J,K,22)
KUSE
(J,K,23)
KUSE
(J,K,24)
KUSE
(J,K,25)
KUSE
(J,K,26)
KUSE
(J,K,27)
KUSE
(J,K,28)
KUSE
(J,K,29)
KUSE
(J,K,30)
destination
56
-------
Table 2 - Continued
CARD
JREAD FORMAT COLUMNS
DESCRIPTION
VARIABLE
NAME
OR ARRAY
IDENTIFIER
12 76-77 Month 2 digits
12 79-80 Year 2 digits
J = Node sequence K = Operation sequence in node
2R3
1-6
Read sentinel
ENDATA
2R3 1-6 Read sentinel
13 10-12 Node number
F6.0 16-21 Index in array
F6.0 22-27 Variable or param-
eter
TRNACT
2R3
1-6
Read sentinel
ENDATA
57
-------
discussion of the main program includes pictorial representations of
a three node river basin and a pictorial representation of an opera-
tion sequence in a node. The flow from node to node is represented
with the use of the CONNOD cards; the data contained on these cards
are stored in the array NCTRL(20,3). An example of the use of the
NCTRL array is as follows for the six node river basin shown on
Figure 14.
Node number
NCTRL(1,1)=101
NCTRL(2,1)=102
NCTRL(3,1)=103
NCTRL(4,1)=104
NCTRL(5,1)=105
NCTRL(6,1)=106
Downstream
node number
NCTRL(1,2)=104
NCTRL(2,2)=103
NCTRL(3,2)=104
NCTRL(4,2)=106
NCTRL(5,2)=106
NCTRL(6,2)= 0
Operation
sequence
number
NCTRL(1,3)=3
NCTRL(2,3)=1
NCTRL(3,3)=1
NCTRL(4,3)=1
NCTRL(5,3)=5
NCTRL(6,3)=1
The next set of control cards CONSEQ, as shown in Table 2, controls
the sequence of operations in any particular node, and the data con-
tained on these cards are stored in the arrays NUMSEQ(20) and
ISEQ(20,20). An example node is shown on Figure 14, and the arrays
NUMSEQ and ISEQ would be filled as follows:
58
-------
EXAMPLE OPERATION SEQUENCE-NODE 2
I- Surface inflow from Node I
2- Tributary inflow
-I- Irrigation demand
3- Return flow to aquifer
-2- Transfer of flow aquifer to river
1- Inflow to river from aquifer
5- Tributary inflow
6- Tributary inflow
FIGURE 14
59
-------
NUMSEQ(l) = 9
ISEQ(l.l) = 1
ISEQ(L,2) = 2
ISEQ(1,3) = 0
ISEQ(L,4) = -1
ISEQ(1,5) = 3
ISEQ(1,6) = -2
ISEQ(1,7) = 4
ISEQ(1,8) = 5
ISEQ(1,9) = 6
The control cards CONWAT are used to designate which upstream nodes
can be used to satisfy a deficit of surface flow in a particular node.
It is obvious that the upstream nodes must have a surface reservoir
to supply such a deficit; however, for the purpose of completing these
data cards it is not important to know if the upstream nodes have the
reservoir feature.
The data contained on the CONWAT cards are stored in the array
NDWTR(20,20). The six-node river basin shown on Figure 15 will again
be used to indicate how the array NDWTR would be filled beginning with
Node #101 as node sequence 1. This array would appear as follows:
60
-------
EXAMPLE
RIVER BASIN WITH SIX NODES
NODE
NODE #102
NODE #103
NODE
NODE #105
NODE #106
N j-Operation sequence
number in node
(inflow) of outflow
from upstream node.
FIGURE 15
61
-------
NDWTR(1,1) = 0 (Node 101)
NDWTR(2,1) = 0 (Node 102)
NDWTR(3,1) = 102 (Node 103)
NDWTR(4,1) = 103"
NDWTR(4,2) = 102
-- (Node 104)
NDWTR(4,3) = 101J
NDWTR(5,1) = 0_ (Node 105)
NDWTR(6,1) = 105
NDWTR(6,2) = 104
NDWTR(6,3) = 103
NDWTR(6,4) = 102
NDWTR(6,5) = 101
— (Node 106)
The order in which the CONNOD cards are stacked is the node sequence
as used by the model. The node sequence begins with the node
farthest upstream in the river basin and continues with the next node
downstream until the ending node of the river basin is reached. The
end of this node represents the outflow from the river basin. The
six node example shown on Figure 15 has four branch nodes to the
mainflow nodes numbered 104 and 106. The stacking order for this
example would be nodes 101, 102, 103, 104, 105, and 106.
The title card and the CONNOD cards must be read prior to reading
the CONSEQ and CONWAT cards. The stacking order of the CONSEQ and
CONWAT cards is of no importance once the CONNOD cards are read.
The card stacking order for the read index JREAD = 1-3 is shown on
Figure 16. The stacking order shown for JREAD=2 is constrained by
the following sequences.
62
-------
STACKING OF DATA DECK FOR JREAD = 1-3
OTHER EQUATION CARDS
F I GURE 16
63
-------
The AQUFER card, which contains the title assigned to the aquifer
(subsurface reservoir) in a node, must precede the GRTANK card,
which contains the initial volume and chemistry of the particular
aquifer.
The STORES card, which contains the title assigned to a surface
reservoir in a node, must precede the RESDAT card, which contains
the initial conditions of the reservoir. The RESCHM card must follow
the RESDAT card and contains the initial chemistry of the reservoir.
The COLBEG card must precede the COLSEG cards. The number of segments
in a soil column in a node is contained on the COLBEG card. A
COLSEG card must be read for each and every segment of the soil
column and contain the initial conditions for each soil segment. The
COLSEG cards must be stacked such that the segment nearest the ground
surface is read first and the order continued until the lowermost
segment in the soil column is read last.
The stacking of the above described cards according to node sequence
is not required; however, each node cannot have more than one aquifer,
one surface reservoir, and one soil column. In the same note, a
node may contain as little as none or as many as all three of the
above described features.
64
-------
The equation cards read when the read index JREAD=3 are the coef-
ficients to mathematical equations used to represent the relationship
of any two variables. Mathematical equations are used to represent
the capacity and area of a reservoir as a function of reservoir water
surface elevation, tailwater as a function of reservoir discharge,
etc. These mathematical formulations will be discussed in greater
detail in later portions of this volume as they pertain to specific
subroutines or functions. Here again the stacking of the equation
cards with reference to node sequence is not required.
The next group of cards in the stacking order are those read when the
read index JREAD=4. The stacking diagram for this read index is
shown on Figure 17 . Many of types of data read when read index JREAD=4
will also be read when the read index JREAD=5.
When JREAD=4 the names or titles of all inflows and demands in the
river basin must precede the actual inflow or demand data cards.
This requirement pertains only to the data for the first time frame.
The SURQIN card contains the title or name of the inflow source.
The data, including volume and chemistry of inflow source, are
contained on the QINFLO card.
The SURDEM card contains the title or name of the demand requirement
from a surface or subsurface source. The DEMSUR card contains the
65
-------
STACKING OF DATA DECK FOR JREAD =4-5
ENDATA
TRNACT (JREAD = 5)
'ENDATA
COHUSE (OPTIONAL)(dREAD = ^ OR 5)
DEMGWR (JREAD = H OR 5)
'DEMSUR (JREAD =
CONUSE (OPTIONAL) (JREAD = H OR 5)
DEMSUR (JREAD = ^ OR 5)
SURDEM (JREAD =
QINFLO (JREAD = 1 OR 5)
SURQIN (JREAD =
FIGURE 17
66
-------
volume of demand from a surface source and an indicator if a consump-
tive use or transfer card is to follow; i.e., a CONUSE card. The
DEMGWR card has the same information as that of the DEMSUR card
except it relates to a demand from a subsurface source. When the
QINFLO, DEMSUR, DEMGWR, and CONUSE cards are read for any other time
frame, except the first time frame, they are read when the read
index JREAD=5 and are not preceded by the name cards SURQIN and
SURDEM. The maintenance of node sequence in stacking these cards
is not required.
All cards read during the first time frame are read when the read
index JREAD = 1-4. The transaction cards T.RNACT can be read during
any time frame other than the first time frame when read index
JREAD=5. The transaction cards are used to change one variable or
parameter per card for such physical features as surface reservoirs
and powerplants.
The variables used in this subroutine are either self-explanatory or
are defined in other portions of this volume.
The various arrays used to store titles or names of surface or sub-
surface facilities, sources of inflow, and demand requirements are
simply eight words of alphanumeric information representing these
various titles or names; i.e., KSURNM(20,8)= Titles of surface reser-
voirs and KAQUNM(20,8) = Titles of subsurface reservoirs of one per
67
-------
each node for 20 nodes. The arrays used to store the eight-word
titles of inflow sources and demand requirements are three dimensional
to accommodate ten inflows and ten demands per node; i.e.,
KQINNM(20,10,8) is used to store inflow titles and KDEMNM(20,10,8)
is used to store demand titles.
68
-------
SUBROUTINE WRIT
This subroutine was written for the report requirements for the
Vernal Project, Utah validation study.
A subroutine such as this will have to be written for each specific
river basin studied and in many cases several different reports will
be required for the same river basin. In this respect, some applica-
tions have been made where pertinent arrays were dumped on tape or
other peripheral equipment at the end of each time frame, and an
editing routine was developed for the desired report generation,
The description of the report generating routine for the Vernal Project
Study would not be of a benefit for the general use of the model
and, therefore, will not be included in this volume. However, a
listing of the WRIT subroutine for the Vernal Project is included
in the program listing for the model which is shown in a subsequent
portion of this volume.
69
-------
SUBROUTINE RESOPR
This subroutine is used to simulate the operation of a surface
reservoir in a node. No more than one reservoir can be included in
any one node.
The executive or main program calls this subroutine as it proceeds
through the iterative scheme. The main program recognizes a reser-
voir exists in a node as it inspects the ISEQ array which contains
the operation sequence of the node. A reservoir is designated with
a zero in the ISEQ array.
The first call to the RESOPR subroutine in any time frame represents
one of the three conditions in which the reservoir operation is
simulated. During this first call, or condition one, the reservoir
chemistry of the last time frame is stored in the array BEGRES. All
releases made from the reservoir during the current time frame are
made using the chemistry of the last time frame. The assumption is
that of a complete chemistry mix through the reservoir will occur in
a one-month time frame.
The sediment inflow from the last time frame is used to update the
total sediment deposited in the reservoir and also to update the
minimum capacity of the reservoir. It is obvious that this tech-
nique does not attempt to distribute the sediment deposited in any
one month time frame on the basis of the fluctuating water surface
levels.
70
-------
Under condition one or the first call, it is assumed that all inflows
to the system are transporting sediment in suspension; therefore, the
total sediment inflow to be used in the next time frame updates are
computed during the first call. The sediment inflow computation is
made using the following expression:
B
Qs = AQW
where Q = sediment inflow in tons per day;
Qw = water inflow in cubic feet per second.
The coefficients A and B are predetermined from correlation studies
made concerning the relationship of the sediment inflow rate and the
water inflow rate. The sediment inflow rate is converted to volume
as acre-feet per month using 75 pounds per cubic foot as the density
of sediment (including voids). The computed volume of sediment is
stored in array BEGRES for use in the next time frame.
The criteria for the operation of the reservoir during the first call
or condition one is to store or release water from the reservoir so
that the water surface level or reservoir capacity is as near as
possible to the target levels or capacities as designated by the
rule curve for the particular reservoir which is stored in the array
RULE. The rule curve represents an ideal use of the reservoir holding
capabilities for a one-year operation. Evacuations for flood control
or maintaining certain reservoir levels for power requirements can also
71
-------
be accomplished with the use of the rule curve. An example rule
curve is shown on Figure 18.
Constraints such as maximum or minimum release oftentimes prevent
the attainment of the target levels or capacities during the
first call to the reservoir for a particular time frame.
72
-------
AN EXAMPLE RULE CURVE
LU
LU
_jh-
LU
0
ac
LU
LU
a:
/
/
/
/
/
/
/
/
/
i
f
1
/
1
1
/
1
/
1
JAN
/
/
/
x
FEB
7
MAR
/
APR
/
/
>
MAY
/*
/ \
JUN
\
\
. \
\
\
JUL
,
\
\
\
i
AUG
\
\
\
\
\
SEP
\
N.
\
\
OCT
X
X
X
N
NOV
\
\
i
\
\
\
\
\
\
\
\
\
\
\
DEC
FIGURE 18
73
-------
All other calls to the reservoir for a particular time frame are for
the purpose of releasing water from the reservoir to meet a down-
stream deficit, which represents condition two, or for the purpose
of passing waters released from an upstream reservoir through this
reservoir to meet a downstream deficit, which represents condition
three.
When the reservoir operation is simulated for condition two, releases
are made to meet a downstream deficit. The releases made during this
call are added to the releases made in the previous call. The
constraints of maximum release or minimum capacity may prevent the
reservoir from satisfying the total downstream requirement, and in
such cases the maximum release possible is made and the reservoir is
flagged to bypass further condition two calls to the reservoir during
the particular time frame.
Condition three is the simulation of the reservoir operation where
waters released from an upstream reservoir are passed through a
particular reservoir which is downstream of the reservoir making the
release and upstream of the point of deficit. A call to the reservoir
under this condition will in most cases leave the reservoir in the
same state as it was in the previous call. However, when the reser-
voir release is updated, it may exceed the maximum release constraint
and in such an event the excess is considered as a spill for the
-------
purpose of continuity. It is possible that even though a spill is
indicated, the reservoir level may not be high enough to utilize the
spillway. An inspection of the reservoir level should be examined
and the indicated spill may in effect become a reservoir bypass.
The reservoir chemistry at the end of the previous time frame is used
as the chemistry of reservoir releases made during a particular time
frame. The chemistry of the reservoir is continually updated for
each and every inflow to the reservoir during a particular time frame.
Reservoir Evaporation
The evaporation losses from the reservoir water surface are in most
cases the mean pan evaporation rates for the immediate locale of the
reservoir. The mean rates are on a monthly basis and the variance
from year to year is not considered significant. The evaporation
rates as used by the model are in units of feet per acre of reservoir
water surface area.
The monthly evaporation losses in units of acre-feet are computed by
multiplying the mean reservoir water surface area by the evaporation
rate for the particular month. In this respect, the model can easily
be modified to include a more sophisticated simulation of evaporation
losses as more information or better techniques become available. The
mean water surface area is simply the mean of the area at the beginning
of a time frame and that of the ith iteration.
75
-------
The evaporation rates on a monthly basis are part of the model input
data set and are stored in the array CONEVP for a 1-year period.
Reservoir Bank Storage
In some cases surface reservoirs have significant water absorption
into the soils along the reservoir banks and in such cases this
phenomenon is called "bank storage." As the reservoir level increases
in elevation, the waters entering the bank are considered as addi-
tional reservoir storage. When the reservoir levels decrease in
elevation, waters stored in the banks are released through capillary
and/or other soil-water actions.
The simulation of this phenomenon is accomplished with the use of
bank storage coefficients which are based upon experimental data and
are stored on a monthly time frame in the array BANKCO. Considera-
tion of "bank storage" is usually reserved for large reservoirs. The
bank storage coefficient is multiplied by the change in storage (sign
included) to determine the effects of "bank storage." A slight
modification to this concept could be used to simulate reservoir
seepage losses.
Capacity-Elevation Curves
The model simulates the physical water storage capacity of a surface
reservoir with a series of second degree polynomials which are
stored as base elevation and three coefficients in the RDATA
76
-------
array. The first derivative of these polynomials is used to compute
the reservoir surface area. The polynomials relate reservoir
capacity as a function of depth above the base elevation. The
coefficients are determined from a stand alone computer program
titled "CARCAP." However, any curve fitting mechanism would suffice
and serve 'the same purpose.
An example of the capacity-elevation curve with an enlargement of
the splines (e) is shown on Figure 19. Twenty sets of second
degree polynomials can be used to represent the capacity-elevation
curve for the operating range of the surface reservoir. Discon-
tinuities do exist between sets of polynomials; however, none are
significant with respect to the general adequacy of the overall
simulation.
General Information
The model is designed to take advantage of the high calculation
speeds of the computer. Consequently, many iterations may be made
through each node during any one time frame before whole river basin
is balanced. All pertinent data arrays are updated both in volumes
of water and chemistry of water at the end of each iteration. It
follows then that several passes may be made through each reservoir
in the river basin and each pass or iteration will update the RDATA
array.
77
-------
EXAMPLE CAPACITY-ELEVAT I ON
CURVE FOR SURFACE RESERVOIR
o
<
Q_
<
O
ELEVATION
BASE ELEVATION
FIRST EQUATION
BASE ELEVATION
SECOND EQUATION
FIGURE 19
V_BASE ELEVATION
THIRD EQUATION
78
-------
A call will be made to the subroutine PWROPR for each pass through a
reservoir after the reservoir is in balance. In fact, the only
route to subroutine PWROPR which simulates the operation of a power-
plant is from subroutine RESOPR.
A general flow chart for the flow of the RESOPR subroutine is shown
on Figures 20 through 22. The list of variable names and arrays is
shown in Table 3.
79
-------
FLOW CHART FOR
SUBROUTINE RESOPR
ENTER WITH ARGUMENTS QFLOW = VOLUME
AND CHEMISTRY ARRAY OF WATERS
ENTERING OR LEAVING RESERVOIR,
KN = NODE INDEX, MO = MONTH AND
IW = WRITE UNIT DESIGNATION
NO
IS THIS A FIRST CALL OR
PASS FOR THIS RESERVOIR
IN THIS TIME FRAME
YES
UPDATE SEDIMENT
DEPOSIT AND
MINIMUM CAPACITY
DUE TO SEDIMENT
INFLOW OF LAST
TIME FRAME.
COMPUTE SEDIMENT
INFLOW FOR THIS
TIME FRAME. STORE
CHEMISTRY OF
RESERVOIR FROM LAST
TIME FRAME IN ARRAY
BEGRES
SET INITIAL VALUES
AND CONSTRAINTS FOR
THIS CALL OR PASS
FIGURE 20
80
-------
FLOW CHART FOR
SUBROUTINE RESOPR (CONT.)
NO
~C
IS QFLOW (I) EQUAL TO ZERO
YES
RETURN
LESS THAN.
IS QFLOW (I) GREATER THAN
OR LESS THAN ZERO
GREATER THAN
THIS REPRESENTS AN
INFLOW TO THE
RESttYOIR MIX
CHEMISTRY OF INFLOW
WITH'RESERVOIR
CHEMISTRY
YES
IS THIS THE FIRST
PASS OR CALL IN
TIME FRAME
NO
UPDATE RELEASE CHECK
MAXIMUM RELEASE
CONSTRAINT
IF EXCEEDED,
UPDATE SPILL
RETURN
FIGURE 21
81
-------
FLOW CHART FOR
SUBROUTINE RESQPR (CONT)
COMPUTE SEED CAPACITY
TO BEGIN ITERATIVE SCHEME
CAPACITY
IS TARGET LEVEL (RULE CURVE) IN
TERMS OF CAPACITY OR ELEVATION
ELEVATION
RESET SEED WITH
COMPUTED CAPACITY
CALL SUBROUTINE
RESCAP TO GET
CAPACITY
SET RELEASE AND STORAGE VALUES
COMPARE RELEASE WITH MAXIMUM AND
MINIMUM RELEASE CONSTRAINTS.
±
COMPUTE EVAPORATION AND BANK STORAGE
COMPUTE NEW CAPACITY IN SCHEME
COMPARE NEW CAPACITY WITH TARGET CAPACITY, AND
MINIMUM AND MAXIMUM CAPACITY CONSTRAINTS
COMPARE NEW CAPACITY WITH SEED
CAPACITY USING AN EPS I LOU TERM
NO
\
UPDATE VALUES
IN RDATA ARRAY
n«o traiLU* itKM
RETURN
BEEN EXCEEDED
YES
FIGURE 22
82
-------
Table 3
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE RESOPR
Variable
name C*
Description
Remarks
or units
AREA
A
BANKCO
BEGRES
B
BSTOR
CAPBAB
CAPMAX
CAPMIN
CAP
CAPX
CAPY
CHEMRR
Water surface area of reservoir
First coefficient of sediment
inflow equation
Coefficient for use in "bank
storage" computation
Storage array for reservoir
chemistry and sediment inflow at
start of time frame
Second coefficient of sediment
inflow equation
Bank storage value (with sign)
Test used to determine if a
power facility is part of
reservoir complex
Maximum capacity of reservoir
Minimum capacity of reservoir
Computational aid for reservoir
capacity
Seed capacity in balancing
reservoir water budget
Check capacity used in iterative
scheme to balance reservoir
water budget
Entry to subroutine CHEMGR
acres
See data array
BANKCO
See data array
BEGRES
acre-feet
See data array
CAPBAB in
PWROPR sub-
routine
acre-feet
acre-feet
acre-feet
acre-feet
acre-feet
See subroutine
CHEMGR
83
-------
Table 3 - Continued
Variable
name
C*
Description
Remarks
or units
CONEVP
CONST
DAYS
DELTAF
DIFF
ELEV
EPSLON
EVAP
IDAYS
IFIR
IMON
IW
J
KA
Evaporation rate for reservoir
water surface area in feet per
acre
Conversion factor to convert feet
cubed per second to acre-feet
per day = 1.983471
Days in a particular month
Additional release from a
reservoir to meet a downstream
deficit
Absolute difference between seed
(CAPX) and check capacity (CAPY)
Reservoir water surface
elevation
Term used to check convergence
of reservoir balance =
10 acre-feet
Reservoir water surface
evaporation
Data array with days of each
month in a year - FEE = 28 days
Computational aid
Month of calculation (particular
time frame)
Write unit designation
Loop index
Incremented loop index for
storage in RDATA array
See data array
CONEVP
acre-feet
acre-feet
Feet above mean
sea level
acre-feet
2 digits
84
-------
Table 3 - Continued
Variable
name
C*
Description
Remarks
or units
KNODE
KN
K
MARK
MODE
MO
NCTRL
NITER
PREA
PRECAP
PREL
PWROPR
QCFS
QFLOW
QSED
Node number
Part of argument list where
KN = Node index
Loop index
Computational aid
Computational aid
Month designation of particular
time frame
Node control array
Iteration counter for scheme to
balance reservoir water budget
Reservoir water surface area of
previous time frame
Reservoir capacity of previous
time frame
Reservoir water surface eleva-
tion of previous time frame
Enter to subroutine PWROPR
Inflow in cubic feet per second
Computational aid array used in
various arguments
Sediment inflow per time frame
See data array
NCTRL
acres
acre-feet
Feet above mean
sea level
See subroutine
PWROPR
Converted from
acre-feet per
month
acre-feet
85
-------
Table 3 - Continued
Variable
name
C*
Description
Remarks
or units
RDATA
RELAF
RELCFS
RES CAP
RESEL
RULE
SPILL
STORAF
TARGET
TEMCAP
Array used to store surface
reservoir parameters
Reservoir release
Reservoir release
Computational aid (reservoir
capacity)
Computational aid (Reservoir
elevation)
Reservoir rule curve (Ideal
operation levels on monthly
basis)
Reservoir spill
Reservoir storage per call or
pass
Target level (rule curve) as
capacity
Computational aid (reservoir
capacity)
See data array
RDATA
acre-feet
Cubic feet per
second
acre-feet
Feet above mean
sea level
See data array
RULE
acre-feet
acre-feet
acre-feet
acre-feet
86
-------
Table 3 - Continued
DATA ARRAYS IN SUBROUTINE RESOPR
DATA ARRAY BANKCO (20)
BANKCO (J)
J = Node index
= Bank storage coefficient for reservoir
in node "J" expressed as a decimal
DATA ARRAY BEGRES (20,11)
BEGRES (J,l)
BEGRES (J,2)
= Redundant (usually reserved for capacity of
reservoir)
BEGRES
BEGRES
BEGRES
BEGRES
BEGRES
BEGRES
BEGRES
BEGRES
BEGRES
-------
Table 3 - Continued
DATA ARRAY RULE
RULE (J,3)
RULE (J,4)
RULE (J,5)
RULE (J,6)
RULE (J,7)
RULE (J,8)
RULE (J,9)
RULE (J,10)
RULE (J,ll)
RULE (J,12)
RULE (J,13)
J = Node index
(20,13) - Continued
= Target level for March
= Target level for April
= Target level for May
= Target level for June
= Target level for July
= Target level for August
= Target level for September
= Target level for October
= Target level for November
= Target level for December
= Indicator of target level
0 = level in terms of capacity
1 = level in terms of elevation
DATA ARRAY RDATA (20,110)
RDATA (J,l) = Node number
RDATA (J,2)
RDATA (J,3)
RDATA (J,4)
= Sequence number
= Reservoir capacity at end of call in acre-feet
RDATA (J,5)
RDATA (J,6)
= Reservoir water surface elevation at end of
call in feet above mean sea level
= Reservoir water surface area at end of call
in acres
= Reservoir evaporation loss at end of call in
acre-feet
88
-------
Table 3 - Continued
DATA ARRAY RDATA (20,110) - Continued
RDATA (J,7)
RDATA (J,8)
RDATA (J,9)
RDATA (J,10)
RDATA (J,ll)
RDATA (J,12)
RDATA (J,13)
RDATA (J,14)
RDATA (J,15)
RDATA (J,16)
RDATA (J,17)
RDATA (J,18)
RDATA (J,19)
RDATA (J,20)
= Reservoir change in bank storage in acre-feet
(sign included) at end of call
= Reservoir spill (spillway release) in acre-
feet at end of call
= Reservoir release through outlets and turbines
in acre-feet at end of call
= Total reservoir release - RDATA (J,8) +
RDATA (J,9) in acre-feet at end of call
= Maximum reservoir capacity in acre-feet
= Minimum reservoir capacity in acre-feet
= Maximum reservoir release (w/o spill) in
cubic feet per second
= Minimum reservoir release in cubic feet per
second
= "A" coefficient for equation - sediment = f
(flow)
= "B^" coefficient for equation - sediment = f
(flow)
= "Ao" coefficient for equation - sediment = f
(flow)
= "62" coefficient for equation - sediment = f
(flow)
= Spline point for equations Qsj_ = A]_Qw^ and
Qs2 = A2Qw£B2 in cubic feet per second - units
on independent variable flow (Qw)
= Total sediment in reservoir at start of time
frame in acre-feet
RDATA (J,21-100) = Twenty base elevations plus coefficients for
20 second degree polynomials where -
capacity = f (elevation increment)
89
-------
Table 3 - Continued
DATA ARRAY RDATA (20,110) - Continued
RDATA (J,101)
RDATA (J,102)
RDATA (J,103)
RDATA (J,104)
RDATA (J,105)
RDATA (J,106)
RDATA (J,107)
RDATA (J,108)
RDATA (J,109)
RDATA (J,110)
J = Node index
= Concentration of calcium in reservoir at
end of call in ppm
= Concentration of magnesium
= Concentration of sodium
= Concentration of chloride
= Concentration of sulfate
= Concentration of bicarbonate
= Concentration of carbonate
= Concentration of nitrate
= Total salts (TDS) in ppm
= Total salts (TDS) in tons per acre-feet
90
-------
FUNCTION RESCAP
This function is used to determine the surface reservoir capacity
when the reservoir water surface elevation is given. The reservoir
surface elevation is passed as part of the argument along with the
node index and write unit designation.
As mentioned earlier in this volume the reservoir capacity is a
function of water surface elevation. This function is represented
as a second degree polynomial and several polynomials may be used
to describe this relationship from the minimum reservoir capacity
to the maximum capacity or the operating pool of the reservoir.
The data required for each of these polynomials are stored in the
RDATA array with four values entered for each polynomial. The
first of these values is the base elevation pertaining to the
particular polynomial, the second, third, and fourth values are the
first, second, and third coefficients where the polynomial is
represented in increasing powers of the independent variable.
An example of the RDATA array for the Steinaker Reservoir, which is
a surface reservoir in the Vernal Project study, is shown in
Table 4 .
91
-------
Table 4
Polynomials Used for Steinaker Reservoir
Base
elevation
5390
RDATA(KN,21)
5400
RDATA(KN,25)
5410
RDATA(KN,29)
5420
KDATA(KN,33)
5430
RDATA(KN,37)
5450
RDATA(KN,41)
5460
RDATA(KN,45)
First
0.0
RDATA(KN,22)
5.0
KDATA(KN,26)
40.0
RDATA(IN,30)
300.0
RDATA(KN,34)
1093.9286
RDATA(KN,38)
4565.0
KDATA(KN,42)
7188.75
RDATA(KN,46)
Coefficients
Second
0.0
KDATA(KN,23)
1.0
RDATA(KN,27)
6.0
RDATA(KN,31)
46.0
RDATA(KN,35)
113.6071
KDATA(KN,39)
233.0
RDATA(KN,43)
288.875
RDATA(KN,47)
Third
0.05
RDATA(KN,24)
0.25
RDATA(KN,28)
2.00
RDATA(KN,32)
3.35
RDATA(KN,36)
3.00
EDATA(KN,40)
2.90
RDATA(KN,44)
3.925
RDATA(KN,48)
5480
RDATA(KN,49)
5500
RDATA(KN,53)
14546.4286
KDATA(KN,50)
25215.3175
RDATA(KN,54)
444.3571
RDATA(KN,51)
626.7523
KDATA(KN,55)
4.45
RDATA(KN,52)
5.3717
RDATA(KN,56)
KN = Node index
92
-------
An example can be made of the function with the use of the data
shown in Table 4 . Given a reservoir water surface
elevation = 5423.7 feet above mean sea level find the reservoir
capacity that corresponds to this elevation.
A search is first made to find the first base elevation in the
KDATA array that is greater than the given elevation. In this
example that base elevation would be 5430.0. The location of this
base elevation is decreased by four, i.e., KDATA(KN,37) to
RDATA(KN,33) where the base elevation of 5420 is stored. The next
three higher locations contain the coefficients of the required
polynomials, i.e., RDATA(KN,34-36).
The reservoir capacity would be computed as follows:
Let: h = 5423.7 - 5420.0 = 3.7 feet
then: capacity = 300.0 + 46.0 (3.7) + 3.35 (3.7)(3.7) = 516 acre-ft
The above example is offered in lieu of a flow chart. The variable
names are self-explanatory and the data arrays have been previously
explained.
93
-------
FUNCTION RESEL
This function is the converse of Function RESCAP where the reser-
voir capacity is given, and it is necessary to find the reservoir
water surface elevation for the given capacity.
The argument list includes the capacity, reservoir water surface
area, node index, and write unit designation.
The reservoir elevation is determined by the solution of a quadratic
equation once the proper set of coefficients has been found for the
particular capacity. The first coefficient of the second degree
polynomial is actually the capacity at the particular base elevation.
The reservoir water surface area for the given capacity is deter-
mined by solving for the derivative dc/dh.
94
-------
SUBROUTINE PWROPR
This subroutine is used to simulate the operation of a hydro-
electric powerplant with any number of generating units. This
subroutine is called by the reservoir operation subroutine RESOPR.
In this respect all powerplants must have a forebay acting as a
reservoir or a reservoir proper. The constraint of one powerplant
per node is consequently the same as that of one surface reservoir
per node.
The argument list consists of the release made from the reservoir
or forebay in cubic feet per second, the node index, the it-u month,
the number of days in the ifc, month, the average reservoir or
forebay water surface elevation, the reservoir spill or bypass, and
the write unit designation.
The powerplant simulation assumes that all generating units have
identical capabilities and are all operated the same number of
hours at overload capacity and at full gate.
Simple relationships are derived from the manufacturer's performance
curves. These relationships are either first or second degree
polynomials relating horsepower and turbine discharge as functions
of net head. Figure 23 is a graphical representation of these
relationships. Consideration has been given to an optimum
95
-------
440
420
too
380
.360
£340
320
300
280
260
12 16 20 24 28 32 36 40 44 48 52 56 60
100 HORSEPOWER AT DESIGN HEAD OF 400 FT.
Generator
Rating-—
Max. head WO ft.-
Gate
/ Design Head
faneplate Rated
,'\ Head 365 ft.
i \
Mfr's Rated
Head 350 ft:
-Best Efficiency-
^-Full
/
Best Efficiency-
Mi n. Head 260 ft.
i i i i
20 30 40 50 60 70 800 1000 1200 1400 1600 70 80 90 100
1000 HORSEPOWER DISCHAR6E C.F.S. % EFFICIENCY
SAMPLE PREDICTED CHARACTERISTIC CURVES
FIGURE 23
96
-------
combination of operating at best gate and full gate as the net
head varies. However, when simulating with a 1-month time frame
and the assumptions necessary for other pertinent variables made
this consideration an invalid refinement.
Efficiency is computed as the ratio of actual horsepower to the
theoretical horsepower when the net head is less than the head
of maximum discharge. When the net head is equal to or greater
than the head at maximum discharge a constant efficiency is used.
The tailwater elevation is computed by entering a polynomial with
total discharge in cubic feet per second as the independent
variable. Coefficients for this polynomial are stored in the
TDATA array.
The net head is simply the average reservoir water surface eleva-
tion minus the tailwater elevation. Other minor head losses
can easily be inserted (if known) into the simulation with a slight
modification to the subroutine.
A firm power demand curve can be included at each powerplant site.
This curve is represented as kilowatt-hours per month. The curve
is assumed not to change from year to year; however, if changes
are desired they can be inserted with the use of transaction cards.
The curve is stored in the data array PDATA. A sample of a firm
97
-------
power demand curve is shown on Figure 24 . This type of curves can
be used to provide a desired distribution of a system power demand.
In this version of the powerplant simulation no provision is made
to increase reservoir or forebay release if power generated is
less than firm demand (iterative procedure) or to halt computation
if firm demand is not met. All power generated greater than firm
demand is considered as secondary power.
An iterative scheme is used in determining power generation as
the tailwater elevation and therefore the net head changes per
iteration. The epsilon term in the iterative scheme is the flow
through the turbines and is presently set at 5 cubic feet per second.
The number of hours of operation is determined as the ratio of
total water released from the reservoir to the total water released
through the turbines multiplied by the number of hours in the
particular month. Whenever this ratio is greater than 1.0 it is
set to 1.0 and the difference between the water released through the
turbines and the total water released from the reservoir is assumed
to be released through the outlet works.
All power generated is computed in units of kilowatt-hours.
A general flow chart for this subroutine is shown on Figures 25 and 26
and the list of variable names and arrays is shown in Table 5 .
98
-------
AN EXAMPLE POWER DEMAND CURVE
a
z
UJ
a
/
s
JAN
s\
\
FEB
^
\
\
\
\
\
MAR
\
\
^
\
\
APR
^
\
\
\
\
V
MAY
\ ,
\^r
JUN
/
s
JUL
/
/
f
/
AUG
/
r /
/
/
/
SEP
/
/
/•
OCT
/
/
/
/
NOV
/
/
DEC
MONTH
FIGURE 24
99
-------
FLOW CHART FOR
SUBROUTINE PWROPR
ENTER WITH ARGUMENTS ACFS = TOTAL WATER
RELEASED FROM RESERVOIR, JAP = NODE INDEX,
IMON = ITH MONTH, DAYS = DAYS IN ITH MONTH,
ELEVI = AVERAGE RESERVOIR WATER SURFACE
ELEVATION, SPILL = SPILL FROM RESERVOIR AND
IW = WRITE UNIT DESIGNATION
YES
IS THIS THE VERY FIRST TIME IN RUN THAT
ANY POWER PLANT IS TO BE SIMULATED
NO
INITIALIZE
CALL INDEX
FOR ALL
POWER PLANTS
YES
IS THIS THE VERY FIRST TIME IN RUN
THAT A PARTICULAR POWER PLANT IS TO
BE SIMULATED
NO
COMPUTE GENERATOR CAPABILITY
AND MAXIMUM HORSEPOWER OF
EACH TURBINE COMPUTE
MAXIMUM HEAD AND MAXIMUM
DISCHARGE AT MAXIMUM
HORSEPOWER. COMPUTE CONSTANT
EFFICIENCY AT HEADS
GREATER THAN MAXIMUM HEAD.
THESE ARE ALL ONE TIME
COMPUTATIONS
FOR EACH RUN FOR
EACH POWERPLANT
FIGURE 25
100
-------
FLOW CHART FOR
SUBROUTINE PKROPR (COHT.)
YES
IS HEAD LESS THAN
MINIMUM HEAD
NO
SET TAILWATER
ELEVATION EQUAL TO
STREAMBED ELEVATION
AND SET ALL COMPUTATIONAL
PARAMETERS EQUAL TO ZERO
RETURN
COMPUTE NUMBER OF
HOURS EACH UNIT
MUST OPERATE IN
MONTH. COMPUTE
POWER IN KWH AND
COMPARE WITH FIRM
DEMAND. COMPUTE
SECONDARY POWER
STORE ALL COMUTATIONAL
PARAMETERS INTO
ARRAY CAPBAB
RETURN
SET SEED FOR
TAILWATER ELEVATION
COMPUTATION
COMPUTE TAILWATER
ELEVATION AND
HEAD. ON BASIS OF
HEAD BEING
GREATER OR LESS
THAN HEAD AT
MAXIMUM HORSEPOWER
OR GREATER THAN
MAXIMUM HEAD
COMPUTE HORSEPOWER
AND DISCHARGE
MULTIPLY DISCHARGE
BY NUMBER OF UNITS AND
AND COMPARE WITH
SEED USED FOR THE
TAILWATER COMPUTATION
FIGURE 26
101
COMPARISON GREATER
THAN THE EPSLON
YES
RESET SEED
-------
Table 5
VARIABLE NAMES OR DATA ARRAYS FOR SUBROUTINE PWROPR
Variable
name
C*
Description
Remarks
or units
ACFS
CAPBAB
CA
CONST
DAF
DAYS
DENS
ELEVI
ENERGY
EPSLON
FOFX
GENCAP
GENEFF
GHEAD
Total release from reservoir or
forebay in cubic feet per second
Data array for pertinent powerplant
parameters
Temporary storage for coefficients
of polynomials
Computation constant to convert
cubic feet per second to acre-feet
per day
Computation constant = cubic
feet in 1 acre-foot
Days in a particular month
Computation constant = density
of water = 62.4 pounds per cubic
foot
Average water surface elevation
of reservoir or forebay
Coefficients for polynomial used
to compute discharge through
turbine at full gate as a function
of net head
Difference in flow through
turbines as iterative scheme
converges
Entry to subroutine to evaluate
polynomials
Generator capacity in kilowatts
Generator efficiency in percent
Gross head for turbine operation
See data array
CAPBAB
Feet above mean
sea level
See data array
ENERGY
5 cubic feet
per second
See subroutine
FOFX
Feet
102
-------
Table 5 - Continued
Variable
name
C*
Description
Remarks
or units
HDMAX
HDMIN
HOURS
HPF
HPMAX
HPTOKW
IFIR
IMON
IW
JAP
KA
KNODE
KP
K
KSET
KWH
KWHSEC
Maximum head for turbine operation
Minimum head for turbine operation
Hours of operation for each turbine
Computation constant to convert foot-
pounds per second to horsepower
Maximum horsepower for turbine
operation
Computation constant to convert
horsepower to kilowatts
Computation index designating entry
sequence to subroutine
Month designation of particular
time frame
Write unit designation
Node index
Loop index
Node number
Node index in loop
Loop index
Computation index of head with
respect to head at maximum horse-
power or at maximum head
Power generation for all units
combined in kilowatt-hours
Secondary or dump power generated
for all units combined in
kilowatt-hours
Feet
Feet
103
-------
Table 5 - Continued
Variable
name
C*
Description
Remarks
or units
MODE
MOLE
NCTRL
PDATA
RATIO
RCFS
RELAF
Computational index used for each
iteration into subroutine for
each node
Alpha designation of source of
call in obtaining root of equations
Control array - See previous
definition of this array
Data array containing firm power
requirements on a monthly basis
Ratio of total discharge released
from reservoir to total discharge
through all turbines combined
Second guess in iterative scheme
for power generated per unit in
cubic feet per second
Discharge through all turbines
combined in acre-feet
See subroutine
ROOTWO
See data array
PDATA
RKWH
Firm power demand for particular
month in kilowatt-hours
ROOTWO
SPILL
TCFS
TDATA
Entry into subroutine to
determine root of polynomial
Reservoir release not passed
through turbines or outlet works
Discharge through turbine at
maximum horsepower in cubic feet
per second
Coefficients for polynomial
describing tailwater elevation
as a function of discharge
through turbines and outlet
works in cubic feet per second
See subroutine
ROOTWO
Acre-feet
Full gate
See data array
TDATA
104
-------
Table 5 - Continued
Variable
name
C*
Description
Remarks
or units
TEFF
TEST
THEAD
THEORY
THP
TKW
TOTCFS
TOTKW
TSPILL
TWEL
UNITS
XT
Turbine efficiency in percent
Computational aid = average reservoir
water surface elevation minus stream-
bed elevation
Head at maximum horsepower
Computational aid representing
turbine horsepower assuming
100 percent effiency
Computed horsepower of turbine from
solution of polynomial - where
horsepower is a function of net head
Computation aid representing kilo-
watts generated per unit
Total discharge through all turbines
combined in powerplant in cubic feet
per second
Total generation of all units combined
in powerplant in kilowatts
Reservoir spill or bypass in units
of cubic feet per second. Used only
when spillway or bypass unit dis-
charge affects tailwater
Tailwater elevation
Number of units in powerplant,
i.e., turbines and generators
Independent variable representing
total discharge in cubic feet per
second in polynomial where tail-
water elevation is a function
of XT
Feet
Feet
Feet above mean
sea level
105
-------
Table 5 - Continued
DATA ARRAYS IN SUBROUTINE PWROPR
DATA ARRAY ENERGY (20,8)
This data array contains the coefficients for the polynomial where
turbine discharge is a function of net head in the first four locations,
i.e., ENERGY (J,l-4). The polynomial is as follows:
2 3
Discharge = AI + AH + Ah + Ah , where h = net head in feet. Then:
ENERGY (J,l) = A1
ENERGY (J,2) = A2
ENERGY (J,3) = A3
ENERGY (J,4) = A^
The second portion of the data array, i.e., ENERGY (J,5-8) contains the
coefficients for the polynomial where turbine horsepower is a function
of net head. The polynomial is as follows:
2 3
Horsepower = AX + A2h + Ash + A^h , where h = net head in feet. Then
ENERGY (J,5) = Ax
ENERGY (J,6) = A2
ENERGY (J,7) = A3
ENERGY (J,8) = A^
J = Node index
106
-------
Table 5 - Continued
DATA ARRAY TDATA (20,4)
This data array contains the coefficients for the polynomial where tail-
water elevation" is a function of total discharge through powerplant and/or
outlet works in cubic feet per second. The polynomial is as follows:
Tailwater elevation = A: + A2Q + A3Q2 + A^Q3, where Q = discharge in
cubic feet per second. Then:
TDATA (J,l) = A! (streambed elevation)
TDATA (J,2) = A2
TDATA (J,3) - A3
TDATA (J,4) = A^
J = Node index
* Tailwater elevation is feet above mean sea level
DATA ARRAY PDATA (20,12)
PDATA (J,l) = Firm power demand in kilowatt-hours for January
PDATA (J,2) = Firm power demand for February
PDATA (J,3) = Firm power demand for March
PDATA (J,4) = Firm power demand for April
PDATA (J,5) = Firm power demand for May
PDATA (J,6) = Firm power demand for June
PDATA (J,7) = Firm power demand for July
PDATA (J,8) = Firm power demand for August
PDATA (J,9) = Firm power demand for September
107
-------
Table 5 - Continued
DATA ARRAY PDATA (20,12) Continued
PDATA (J,10) = Firm power demand for October
PDATA (J,ll) = Firm power demand for November
PDATA (J,12) = Firm power demand for December
J = Node index
DATA ARRAY CAPBAB (20,20)
CAPBAB (J,l) = Number of units (turbines) in powerplant
CAPBAB (J,2) = Generator rating in kva
CAPBAB (J,3) = Power factor (dimensionless)
CAPBAB (J,4) = Percent overload allowed each generator
CAPBAB (J,5) = Generator capability in kilowatts
CAPBAB (J,6) = Generator efficiency in percent
CAPBAB (J,7) = Maximum net head in feet
CAPBAB (J,8) = Minimum net head in feet
CAPBAB (J,9) = Maximum horsepower of turbine at overload
CAPBAB (J,10) = Head at maximum horsepower (full gate) in feet
CAPBAB (J,ll) = Discharge through turbine at maximum horsepower
in cubic feet per second
CAPBAB (J,12) = Turbine efficiency at maximum horsepower
in percent
CAPBAB (J,13) = Hours in month each turbine operates at overload
CAPBAB (J,14) = Total power from powerplant for month in
kilowatt-hours
108
-------
Table 5 - Continued
DATA ARRAY CAPBAB (20,20) Continued
CAPBAB (J,15) = Total discharge through all units in acre-feet
CAPBAB (J,16) = Tailwater elevation in feet above mean sea level
CAPBAB (J,17) = Gross head (same as net head) in feet
CAPBAB (J,18) = Firm power generated in kilowatt-hours
CAPBAB (J,19) = Secondary or dump power generated in kilowatt-hours
CAPBAB (J,20) = Actual turbine efficiency in percent
J = Node index
109
-------
SUBROUTINE DTABLE
This subroutine is used to rank all the reservoirs within the river
basin for use in the determination of which reservoir should be
called first, second, third, etc., on the basis of the least loss in
elevation with the greatest incremental change in capacity. Ranking
in this manner is for the purpose of maintaining the maximum net
head for power generation.
When the ranking is completed the reservoirs are called in the
resulting sequence to satisfy a deficit in meeting a downstream
requirement. It is obvious that only reservoirs that are accessible
and upstream of a deficit will be ranked in the calling sequence.
The sequence in considering each node in the river basin, i.e.,
beginning with node farthest upstream and the operational node of
the simulation, will insure that none of the deficit calls to a
particular reservoir represents a first call to the reservoir.
Recalling the use of second degree polynomials to relate reservoir
capacity (dependent variable) to the reservoir water surface eleva-
tion and delta increments of elevation (independent variable), it
can easily be shown that the first derivative (^c/dh) is the
reservoir water surface area. Consequently, the reservoir with the
largest water surface area will have the highest rank.
110
-------
This subroutine in effect is a means of selecting withdrawals from
the various surface reservoirs or a decision table for the operation
of the reservoirs within the river basin.
Other factors may be considered in designing a decision table such
as that described for this subroutine; however, the logic used in
this subroutine has proven to be satisfactory for most river basin
simulations. In any event the task of modifying the logic would be
very simple if some other criteria were desired in setting a sequence
of ranking reservoirs for the purpose of selective withdrawals. The
variable names and arrays are shown in Table 6.
The argument list includes the node index, an array for the storage
of node numbers of nodes containing reservoirs upstream and accessible
to meet a deficit in a downstream node in a ranked sequence, an
index of the number of upstream reservoirs that can be used to meet
the deficit, and the number of operational sequences in the node
where the deficit exists.
Ill
-------
Table 6
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE DTABLE
Variable
name
C*
Description
Remarks
or units
AREA
1C
ISEQ
KFA
KN
K
KS
LTEMP
MAX
MIN
NCTRL
NDWTR
NFIND
Water surface area of reservoir
in acres
Index used as an indicator of
number of upstream reservoirs
Operational sequence array for
each node
Storage array of node numbers of
nodes upstream and accessible to
meet deficit
Node index of upstream node
Loop index
Index of number of operation
sequences in node of deficit
Computational aid
Computational aid used in rank-
ing technique (maximum)
Computational aid used in rank-
ing technique (minimum)
Control array for the simulation
sequence of a river basin
Control array for the simulation
sequence of a river basin
Entry to subroutine NFIND
See previous
definition of
this array (page 62)
See previous
definition of
this array (page 60)
See previous
definition of
this array (page 62)
See subroutine
NFIND
112
-------
Table 6 - Continued
Variable
name
C*
Description
Remarks
or units
RDATA C Array used to store reservoir
data
TEMP Computational aid
See previous
definition of
this array (Table 3)
ARRAY AREA (20) = Temporary storage array for reservoir areas
upstream and accessible to meet downstream
deficit
ARRAY KFA (20) = Storage array of node numbers of ranked reser-
voirs in descending order of reservoir areas
The dimension (20) represents the maximum number of nodes that can
be used in simulation model.
113
-------
FUNCTION NFIND
This function provides a simple method of determining the opera-
tional sequence of a node in the river basin. The argument list
contains only the node number. A search is made of the node control
array NCTRL to find the operational sequence of the node.
The need to find the operational sequence of the node occurs several
times in the simulation model; therefore, this simple function is
well justified.
114
-------
SUBROUTINE RFCOMP
This subroutine is among the most essential of all techniques used
in the simulation. It is primarily designed to simulate the magni-
tude and chemistry of return flow from either an irrigation or
municipal and industrial demand (diversion). The subroutine also
provides the mechanism to simulate exchanges of flow and chemistry
that must exist between surface and subsurface facilities to
maintain a hydraulic balance in each node.
An earlier version of this subroutine was more complex requiring
greater quantities of data and containing several more empirical
relationships in an attempt to define in a rigorous manner the
deposition and timing (lagging) of the excess or return flow. Many
applications were made using this early version, and it was concluded
that uncertainties concerning the deposition and timing of returnable
flows were often the result of lack of adequate data and, in some part,
lack of technology. Consequently, much of the sophistication of the
earlier version was discarded in favor of the simpler approach as
contained in this subroutine where the exchange mechanism was
adopted to "black box" many of the uncertainties and not destroy the
validity of the simulation.
The argument list includes an array where magnitude and chemistry
are stored for either incoming or outgoing waters, the node index,
115
-------
the operational sequence in the node, an index signifying the type
computation to be performed by the subroutine, the i., month, a
constant, and the write unit designation.
The application of the simulation model is highly dependent upon
an accurate determination of the types of uses of the water resources
and the deposition of the excess or return flow. A use as defined
in this simulation model is that water that is actually consumed to
support plant or animal life (in effect leaves the system or river
basin). An industrial requirement may or may not be considered a
use by this definition.
The supply of an irrigation or municipal and industrial need begins
with a point of diversion on a river, stream, or main service canal.
This diversion is transported by means of smaller canals either lined
or unlined to a farm area, municipality, or industrial site. As
these waters are being transported, losses occur by evaporation and
seepage into the soils of the unlined conveyances. Further losses
occur after the waters reach the farm area by seepage losses of
unlined laterals and by irrigating crops (plants) in excess of their
needs. Additional losses and excesses are also experienced after
waters reach a municipality. Most of the municipal losses are
collected in a sewer system and treated before returning to the
river basin. Simulation of the treatment facilities of these
116
-------
mu
micipal losses has not been included in this version of the model;
however, such a simulation is easily incorporated into the model if flow
through rates and types of treatment are known.
The preceding brief background information has been included in this
narrative to aid the user in understanding some of the philosophy
used in the design of this subroutine.
Consumptive Use
The consumptive use as treated in this subroutine is actually that
amount of water that leaves the river basin either in support of
plant life (crops) or in support of a municipality or industrial
function.
The consumptive use of a municipal and/or industrial complex must
be part of the input data in units of acre-feet per month. This
consumptive use is usually obtained from historical data concerning
present and past requirements and should include modifications
obtained from estimates of future population and industrial growth.
The consumptive use of irrigated area must also be part of the input
data in units of acre-feet per month. An irrigated area is usually
composed of several type crops with different areas for each, and
the sum of these areas can represent the total area of one irrigation
requirement.
117
-------
Three methods are most commonly used to compute the consumptive use
of an irrigated area and are stand-alone models or electronic com-
puter programs. Two of these common methods use the Blaney-Criddle
formula or a formula developed by R. L. Lowry, Jr., and A. F. Johnson.
The third method which was developed at the Engineering and Research
Center of the Bureau of Reclamation in Denver, Colorado, is a stand-
alone computer application titled EVPCOM (see reference 1). The
consumptive use which is a theoretical value is considered as being
supplied from soil moisture storage which is maintained by either
irrigation applications or precipitation. The model then treats the
changes in soil moisture storage as a lumped parameter on the premise
that such changes are adequately reflected whenever the consumptive
use requirement is met and on a monthly time frame the net effect of
soil moisture storage changes are insignificant or zero.
The Penman equation is used to determine the potential transpiration
rate in units of inches per day, and the crop curves of Jensen and
Haise (see reference 2) are used to modify potential transpiration
rate to reflect the amount of crop growth from the date of planting to
the date of harvest. These modified daily rates are summed for all crops
irrigated in the irrigated area for each day of each month. The
summed rates for each day of the month are again summed for the
month and converted to units of acre-feet per month. When precipi-
tation is significant, it can be included in the above computation
118
-------
or it can be treated as a tributary inflow to the node in which the
irrigated area is located.
The consumptive use for an irrigated area usually is the combination
of both the nonbeneficial (miscellaneous vegetation) and the beneficial
(crops) needs. The computation of the nonbeneficial use in model
developed by the Bureau of Reclamation is the same as that for the
beneficial use with the exception of the crop curve used which is
dependent upon the type of vegetation.
Demand (Diversion Requirement)
The demands are computed in the model as functions of the consumptive
use and the returnable losses. The returnable losses are obtained
from historical data and are expressed as percentages of the diver-
sion (demand). Conversely, if the consumptive use is known, the
diversion or demand can be determined by dividing the consumptive
use by the factor (1.0 - losses expressed as a decimal).
Return Flow
The return flow is simply the returnable losses and in the simulation
it is computed as the difference between the demand and the consumptive
use.
The chemistry of the return flow is determined using the premise
that the consumptive use (transpiration) leaves the river basin with
119
-------
zero chemistry and the concentrations of the waters diverted are
increased on this basis. As an example:
Let: V, = Volume of water diverted.
C, = Concentration of water diverted.
V9 = Volume of consumptive use.
C9 = Concentration of consumptive use = 0.
V_ = Volume of return flow.
Co = Concentration of return flow.
Then: V0 = V, - V
3 1 2
(C1V1
=
The array CONUSE (EQUIVALENCE KUSE) contains information designating
the destination (maximum of three) and the ratio of the total return
allocated to each destination.
Shortage Index
The shortage index is a function of the demand (diversion), and
whenever a demand cannot be supplied within the limits of the
shortage index an error message will be printed and the simulation
will be terminated.
120
-------
Summary
This subroutine has been designed to allow for maximum flexibility
in modification for variable demand growth without resorting to a
large increase in the data set. Also, if the demands for each month
do not change from year to year, the demand and consumptive use cards
are necessary for the first year only.
The use of the three destinations for return flow has proven to be a
valuable tool in obtaining a hydraulic balance between nodes, and
as a tool it can be used to improve the chemistry comparisons at
the end of each node.
It has been found that a check point or response point located as
far downstream in a node as possible is essential in obtaining valid
results with the use of the exchange mechanism between surface and
subsurface facilities. In fact, once a consistent criteria for the
hydraulic balance has been established, the chemistry or salinity
comparisons (observed versus predicted) will converge within reasonable
limits.
The flow charts for this subroutine are shown on Figures 27 and 28.
The description of variable names or data arrays is given in Table 7.
121
-------
SUBROUTINE RFCOMP
REFERENCES
1. Huntley, C. W., "Computation of Consumptive Use," EVPCOM,
Bureau of Reclamation Internal Document, June 1974.
2. Jensen, M. E., Wright, J. L., and Pratt, B. J., "Estimating
Soil Moisture Depletion from Climate, Crops and Soil Data,"
Sprinkler Irrigation, 3rd Edition, 1969, Sprinkler Irrigation
Association, Washington, D.C., Chapter 5, "Plant and Irrigation
Water Requirements."
122
-------
FLOW CHART FOR
SUBROUTINE RFCOMP
ENTER WITH ARGUMENTS
QFLOW = ARftAY CONTAINING
MAGNITUDE AND CHEMISTRY
OF INFLOWS OR OUTFLOWS,
NK = NODE INDEX, IP -
OPERATIONAL SEQUENCE IN
NODE, INDY = COMPUTATION
INDEX, MO - ITH MONTH,
CONST = MAGNITUDE OF
DEMAND OR TRANSFER,
IW = WRITE UNIT
INDY = I
DEMAND CAN BE
MET WITHOUT
SHORTAGE. SET
SHORTAGE = 0
INDY = 2
INDY = 33,5
DEMAND CANNOT BE
MET WITHOUT
SHORTAGE. COMPUTE
SHORTAGE ON BASIS
OF CONSUMPTIVE USE
1
SET SHORTAGE
NO
SHORTAGE GREATER
THAN INDEX
YES
COMPUTE MAGNITUDE AND
CHEMISTRY OF RETURN
FLOW. ALLOCATE RETURN
AND SET QINFLO ARRAY
ON BASIS OF ALLOCATION I
NODE AND OPERATIONAL
SEQUENCE IN NODE FOR
QINFLO
.E.
PRINT ERROR
MESSAGE
HALT ALL
SIMULATION
CALL EXIT
RETURN
FIGURE 27
123
-------
FLOW CHART FOR
SUBROUTINE RFCOMP (CONT.)
EXCHANGE MECHANISM SET UP
FLOW TRANSFERS NO RETURN
FLOW COMPUTATIONS
TRANSFER CHEMISTRY AND MAGNITUDE
OF FLOW TRANSFER INTO QINFLO
ARRAY AS TRANSFERS TO RIVER,
AQUIFER OR PONDING POOL
FOR PERCOLATION
RETURN
FIGURE 28
124
-------
Table 7
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE RFCOMP
Variable
name
Description
Remarks
or units
CONST
CONUSE
DEMAND
INDY
IP
IW
JK
KA
KDEMNM
KDES
KNODE
K
KSEQ
Constant used to store magnitude
of demand
Consumptive use array storing
various pertinent parameters
Array used to store magnitude of
demands on monthly basis in
acre-feet
Computation index
Operational sequence in node
Write unit designation
Node sequence in setting QINFLO
array
Loop index
Alpha description of demand or
diversion
Destination index of return and
transfer flows (3 character
alpha)
Node number of return flow
allocation and flow transfer
Loop index
Operational sequence in node of
return flow or flow transfer
Acre-feet
See previous
definition of
data array
CONUSE
See previous
definition of
data array
DEMAND
125
-------
Table 7 - Continued
Variable
name
C*
Description
Remarks
or units
KUSE
MO
NFIND
NK
PER
QFLOW
QINFLO
RATIO
RFLOW
SHORT
TEMCU
VOLA
VOLB
Integer equivalence of consump-
tive use array
1^ month in simulation
Entry to subroutine NFIND
Node index
Distribution ratio for return
flow allocation
Array containing magnitude and
chemistry of inflow or outflow
Array containing parameters of
tributary inflow
Decimal equivalent of conveyance
and farm losses or conveyance
and municipal or industrial
losses
Return or excess flow from a
diversion either irrigation or
municipal and/or industrial
demand
Actual shortage in supplying
ideal demand (diversion) in per-
cent of consumptive use
requirement or ideal demand
Computational aid in computing
consumptive use shortage
Computational aid for volume in
computing chemistry of return
flow
Same purpose as VOLA
Equivalence of
array CONUSE
See subroutine
NFIND
See data array
QFLOW
See previously
defined data
array QINFLO
See data array
RFLOW
Acre-feet
Acre-feet
126
-------
Table 7 - Continued
DATA ARRAY QFLOW (10)
QFLOW (1) = Magnitude of demand or transfer in acre-feet
QFLOW (2) = Concentration of calcium in parts per million
QFLOW (3) = Concentration of magnesium
QFLOW (4) = Concentration of sodium
QFLOW (5) = Concentration of chloride
QFLOW (6) = Concentration of sulfate
QFLOW (7) = Concentration of bicarbonate
QFLOW (8) = Concentration of carbonate
QFLOW (9) = Concentration of nitrate
QFLOW (10) = Concentration of total salts (TDS)
DATA ARRAY RFLOW (10)
RFLOW (1) = Magnitude of return flow in acre-feet
RFLOW (2-10) = Same as QFLOW (2-10)
127
-------
SUBROUTINE CHEMGR
This subroutine simulates the mixing of waters of different chemical
composition as the waters are intercepted at a common point. The
subroutine has seven entry points to facilitate the mixing calls
from several sources.
The argument list consists of an array containing the magnitude and
chemistry of both the incoming and outgoing waters, a constant, the
node index, and the operational sequence of the node.
All of the mixing equations contained in this subroutine are linear
as the following example indicates.
Let: Vj_ = Volume of water to be mixed from the first source.
C, . = Concentration in ppm of the jt^ constituent from
the first water source.
Vo = Volume of water to be mixed from the second source.
C-- = Concentration in ppm of the jt^ constituent from the
second water source.
Vo = Combined volume of first and second sources.
GOJ: = Combined concentration in parts per million of the
j^ constituent of the first and second sources.
Cv (V1CH + V2C2l)
JJ v3
It is a simple task to modify any of the mixing equations to include
nonlinear relationships if information and data justify such a
representation.
128
-------
A flow chart of this subroutine is shown on Figures 29 and 30 .
The variable names or arrays are shown in Table 8 .
Subroutine CHEMGR is a multiple entry device with six entry points,
These entry points are designated on Figures 29 and 30 with an
appropriate description as to function of each entry point.
129
-------
FLOW CHART FOfl
SUBROUTINE CHEMGR
ENTER WITH ARRAY
QFLOW = MAGNITUDE AND
CHEMISTRY OF RIVER,
INCOMING OR OUTGOING WATER,
CONST = CONSTANT, KN =
NODE INDEX, KS = OPERATIONAL
SEQUENCE IN NODE
ENTRY
CHEMGR
ENTRY
CHEMTT
ENTRY
CHEMRV
MIX INCOMING WATER
( QFLOW ) WITH WATER IN
AQUIFER. UPDATE
CHEMISTRY OF AQUIFER.
SET NEW VOLUME OF
WATER IN AQUIFER
MIX INCOMING WATER
( QFLOW ) WITH WATER IN
PONDING POOL ABOVE
SOIL COLUMN. UPDATE
CHEMISTRY OF PONDING
POOL. SET NEW VOLUME
OF WATER IN POOL
MIX INCOMING WATER
(QFLOW ) WITH WATER IN
RIVER. UPDATE
CHEMISTRY OF RIVER.
SET NEW VOLUME
OF WATER IN RIVER
RETURN
FIGURE 29
130
-------
FLOW CHART FOR
SUBROUTINE CHEMGR (CONT)
ENTRY
CHEMRS
EN
CHEI
TRY
MAQ
MIX INCOMING WATER
(QFLOW ) WITH WATER
IN SURFACE RESERVOIR.
UPDATE CHEMISTRY OF
RESERVOIR. SET VOLUME
OF RESERVOIR
MIX WATER IN
AQUIFER WITH
INCOMING WATER
( QFLOW ) PLACE
UPDATED CHEMISTRY
BACK INTO ( QFLOW)
SET NEW VOLUME
ENTRY
CHEMBT
MIX INCOMING WATER
( QFLOW) I.E. WATER
IN SYMBOLIC TANK
ABOVE AQUIFER
(USED IN SOIL
COLUMN SIMULATION)
WITH WATER IN
AQUIFER. UPDATE
CHEMISTRY OF
AQUIFER. SET
NEW VOLUME
RETURN
FIGURE 30
131
-------
Table 8
VARIABLE NAMES OR ARRAYS
FOR SUBROUTINE CHEMGR
Variable
name
C»t-
'*
Description
Remarks
or units
CONST
GRTANK
KA
KN
K
KS
QFLOW
QINFLO
RDATA
TTANK
VOLA
VOLB
VOLUME
Constant that represents magni-
tude of flow in argument list
An array that contains volume and
chemistry of aquifer
Loop index
Node index
Loop index
Operational sequence in node
An array in argument list with
magnitude and chemistry of flow
Array containing inflow
parameters
Array containing reservoir
parameters
Array containing volume and
chemistry of ponding pool above
soil column
Computational aid for volume
Computational aid for volume
Volume of combined flows,
reservoirs, etc., after mix
See previously
defined array
GRTANK
See previously
defined array
QFLOW
See previously
defined array
QINFLO
See previously
defined array
RDATA
See previously
defined array
TTANK
Acre-feet
Acre-feet
Acre-feet
132
-------
SUBROUTINE CONVER
This subroutine is used to convert units of concentration from
milliequivalents per liter to parts per million. Parts per million
are the units used throughout the simulation with the exception of
that portion simulating waters percolating vertically through a
soil column.
The argument list contains an array of the constituents to be con-
verted and an index signifying if the incoming units of concentration
are in parts per million or in milliequivalents per liter. When
this index is equal to one, the units are in milliequivalents per
liter.
The subroutine checks the balance of anions and cations with the use
of the concentrations of the eight constituents in units of equiva-
lents when such concentrations are given as part of the input data
set. The difference between the sum of cations and the sum of anions
(if any) is added to the concentrations of the calcium or the sulfate
constituent dependent upon which of the sums is larger. It is obvious
that not all of the constituents in a sample can be represented as
eight constituents; therefore, this check is not an explicit repre-
sentation of a mass balance. The use of this premise does not
extend a significant error and does provide continuity for subsequent
requirements for a mass balance. When the concentrations of
133
-------
constituents are not a part of the input data set, the balance of
ions is bypassed as well as the computation of total salts (TDS).
The total salts are computed with the use of the concentrations of
the eight constituents when such data is part of the data set. This
computation is performed in units of parts per million, and one-half
of the concentration of the bicarbonate constituent is used in the
computation.
The variable names or arrays are shown in Table 9.
134
-------
Table 9
VARIABLE NAMES OR ARRAYS
FOR SUBROUTINE CONVER
Variable
name C*
Description
Remarks
or units
DIFF
EQVWT
Difference between sum of
cations and anions
Equivalent weight array for
constituents
ppm
See data array
EQVWT
SUMA
SUMC
TDS
TEMP
Sum of anions
Sum of cations
Total salts
Incoming array of
constituents
ppm
ppm
ppm
See data array
TEMP
DATA ARRAY TEMP (15)
TEMP (2) = Concentration of calcium
TEMP (3) = Concentration of magnesium
TEMP (4) = Concentration of sodium
TEMP (5) = Concentration of chloride
TEMP (6) = Concentration of sulfate
TEMP (7) = Concentration of bicarbonate
TEMP (8) = Concentration of carbonate
TEMP (9) = Concentration of total salts
TEMP (10) = Total salts
DATA ARRAY EQVWT (8)
EQVWT (1) = Equivalent weight of calcium = 20.04
135
-------
Table 9 - Continued
DATA ARRAY EQVWT (8) - Continued
EQVWT (2) = Equivalent weight of magnesium = 12.16
EQVWT (3) = Equivalent weight of sodium =23.0
EQVWT (4) = Equivalent weight of chloride = 35.46
EQVWT (5) = Equivalent weight of sulfate = 48.03
EQVWT (6) = Equivalent weight of bicarbonate =30.0
EQVWT (7) = Equivalent weight of carbonate = 61.01
EQVWT (8) = Equivalent weight of nitrate = 62.008
136
-------
SUBROUTINE VERNAL
This subroutine is used for the simulation of the operational criteria
for the Vernal Project Study. A subroutine such as this will have to
be designed for each river basin simulated using this model.
The operational criteria differ considerably from one river basin
to another, and any attempt to generalize this portion of the model
would be futile.
It is believed that a near maximum use of all the model features
would require several subroutines, to include all the needs of a
complex operational criteria. In such a case it is recommended that
a monitor subroutine with one or two calls from the main or executive
program be utilized to maintain the flexibility of the model.
This subroutine as written contains the very simple operational
criteria for the three node Vernal Project. The flow chart of the
subroutine is shown on Figure 31.
137
-------
FLOW CHART FOR
SUBROUTINE VERNAL
NO
ENTER WITH ARRAY QFLOW
= MAGNITUDE AND CHEMISTRY
OF PREDICTED RIVER FLOW
AT A RESPONSE POINT, ARRAY
OBSERV = OBSERVED FLOW AT
THE SAME RESPONSE POINT,
KN = NODE INDEX, IS =
OPERATIONAL SEQUENCE
IN NODE, MO = MONTH INDEX.
SET INDICES "MAX", "MIN" AND
"KA" FOR EACH NODE. WHERE
MIN = OPERATION SEQUENCE OF
FIRST OBSERVED OUTFLOW AT RESPONSE
POINT IN NODE, MAX = OPERATION
SEQUENCE OF LAST OBSERVED
OUTFLOW AT RESPONSE POINT IN
NODE AND KA = OPERATION SEQUENCE
OF TRANSFER FROM RIVER TO AQUIFER.
OPERATION SEQUENCE OF TRANSFER
FROM AQUIFER TO RIVER = KA+I.
IS MAGNITUDE OF PREDICTED
FLOW AT RESPONSE POINT
GREATER THAN THE MAGNITUDE
OF OBSERVED FLOW
TRANSFER DIFFERENCE
AQUIFER TO RIVER
TRANSFER DIFFERENCE
RIVER TO AQUIFER
RETURN
FIGURE 31
138
-------
FUNCTION ROOTWO
This function was found very useful in the simulation because of the
several requirements throughout the model for the solution of a root for
quadratic equations. The roots of these quadratic equations by defini-
tion must be real and unequal because all the coefficients of the quad-
ratic are real. A criteria that no imaginary roots can exist is signaled
by a -write statement included in the function. The solution of the
quadratic equation for a single root is as follows:
-b +
X =
where: ax2 + bx + c = o and a ^ o
The argument list for this subroutine consists of the node index, opera-
tional sequence in the node, an array of the coefficients in decreasing
powers of "x", a six alpha character indentifier designating the origin
of the call and the write unit designation. The variable names and
arrays are shown in Table 10.
139
-------
Table 10
VARIABLE NAMES OR ARRAYS
FOR FUNCTION ROOTWO
Variable
name
C*
Description
Remarks
or units
CA
IW
KNODE
KSEG
MOLE
RAD
An array containing the coeffi-
cients of the quadratic equation
Write unit designation
Node index
Operational sequence in node
Six character alpha designation
of calling routine
Radical term in the solution of
a quadratic equation
See data array
GA
DATA. ARRAY CA(lO)
CA(l) = "a" coefficient
CA(2) = 'V coefficient
CA(3) = "c" coefficient
Quadratic equation = ax + bx + c
140
-------
FUNCTION ROOT
th
This subroutine is used to solve for a real root of an n degree poly-
nomial given that such a root lies between two limits and assuming that
only one real root exists within these limits.
The technique used is a combination of the Newton or Newton-Raphson
method where the first derivative is approximated (Regula falsi) and an
algorithm for the solution of a polynomial (Function FOFX). The sub-
routine was designed primarily for the solution of the more complex
functions; consequently, the necessity for approximating the first
derivative. In the model this subroutine is used to find roots of
those polynomials higher than second degree.
The argument list contains a constant, the upper limit of the root, the
lower limit of the root, a five alpha character designation for the
validity of the root, an array containing the coefficients of the poly-
nomial and a test constant for convergence.
A brief description of the procedure is:
given a polynomial, i.e.
„/ \ n n-j.
f(x) = ax -f a X + ... + a X +
J- 2 n
n-1
a
n+1
where n = degree of polynomial, n > 2
find a root X , where: X . < X < X
r7 mm — r — max
141
-------
The Nevton or Newton-Raphson method is formulated as follows:
f(X )
- - x - n
n f.(xn)
Let the approximation (Regula falsi) of f'(X) "be:
n
f(X ) - f(X n)
fi(X ) ~ n-1
v ;
X - X .
n n-1
Upon substitution then:
Xn+l " Xn " f(Xn) -
Xn -
or:
(X - X . )f(X )
v n n-1' v n'
n
Let the first seed (guess) be:
Xn-l ~ Xmin
X = X _ + 0.1 X .
n n-1 min
Begin iterative procedure with maximum of 2o iterations or until
[Xn-KL - Xn] ^ € '
This function calls the function FOFX for the evaluation of the polynomial
using the algorithm described in function FQFX,
When the root (X ) is within the limits of X . and X the alpha
r' man max ^
designator will be set as TRUE, when the root exceeds the limits the
142
-------
alpha designator will be set as FALSE, and when the number of iterations
exceeds 2o the alpha indicator will be set as NONE.
The method used to evaluate the polynomial requires the coefficients
of the polynomial to be passed in increasing powers of the independent
variable, i.e.:
f(x) = a + a X+aX2+...+aXn
012 n
where n = degree of polynomial.
The indexing format of the coefficient array in the computer applica-
tion prohibits the use of zero; therefore the function will actually
be represented as follows:
f(x) = a -*-aX + aX2-i-...+a,nXn
N i 2 3 n-H
If - f(x) = K = a +a X + a X2 + ... + a ...Xn
X 2 3 Xi»-L
then: f(x) = 0 = a -K + aX+aX2^-... +a ^Xn
^ ' i 23 n+1
The flow chart for this function is shown on Figure 32. The variable
names and arrays are shown in Table 11.
143
-------
FLOW CHART FOR
FUNCTION ROOT
ENTER WITH XMAX = UPPER
LIMIT OF ROQT, XMIN = LOWER
LIMIT OF ROOT, CONST = VALUE
OF POLYNOMIAL OTHER THAN ZERO,
MOLE = ALPH DESIGNATION OF
EXISTENCE OF ROOT, CA = ARRAY
OF COEFFICIENTS AND TESTA =
CONSTANT USED TO TEST
FOR CONVERGENCE.
SET SEED = XA
SET GUESS = XB = SEED + 0.I(XMAX)
CORRECT FIRST COEFFICIENT FOR
VALUE OF POLYNOMIAL OTHER THAN ZERO
EVALUATE POLYNOMIAL
IN TERMS OF "XA" AND "XB"
YES
IS ABSOLUTE DIFFERENCE
BETWEEN SEED AND GUESS
LESS THAN CONVERGENCE
CONSTANT (REGULA FALSI)
NO
/
1
f
RESET SEED
RESET GUESS
FIGURE 32
144
-------
Table 11
VARIABLE NAMES OR DATA ARRAYS
FOR FUNCTION ROOT
Variable
name
C*
Description
Remarks
or units
CA
CONST
DXDP
DX
EPSLON
FOFX
FXA
FXB
IW
K
MOEE
NK
ROOT
Array containing coefficients of
polynomial
Used when f(x) = K = CONST
Variable name for f(Xn) -
Variable name for X - X ..
n n-1
Variable name for [X ,- - X ]
n+1 n
Entry to function FOFX
Variable name for f(X , )
Variable name for f (X )
Write unit designation
Loop index
Five- character alpha designator
of root
Iteration counter
Variable name for X
See data array
CA
See Subroutine
FOFX
TESTA
XA
XB
XC
Convergence term for iterative
scheme
Variable name for X
n-1
Variable name for X
n
Variable name for X
n+1
145
-------
-Table_ll - Continued
VARIABLE NAMES OR DATA ARRAYS
FOR FUNCTION ROOT
Variable
name
C*
Description
Remarks
or units
XMAX
XMOT
Variable name for upper limit
of root
Variable name for lower limit
of root
DATA ARRAY CA (lO)
CA (l) = Coefficient for Xn
CA (2) = Coefficient for X11"1
CA(lO) = Coefficient for X°
The number of coefficients will always be one greater than the degree
of the polynomial; however, the coefficient for X must always be
stored in CA(l) and the index of the array CA will increase as the
power of "X" decreases in the polynomial.
146
-------
FUNCTION FOFX
This function was designed to evaluate any polynomial of degree n
from n = 1 to n = 9; however, in the model it is set to evaluate
polynomials from n = 1 to n = 4.
The argument list consists of the array containing the coefficients,
the value of "X" in the polynomial to "be evaluated and an indicator
used to clear the coefficient array after the evaluation of the
polynomial.
The function uses the following algorithm in the evaluation of the
polynomial:
given - 'f(x) = a + a X + a X2 + ...+aX
012 n
•which can be rewritten as follows:
f(x) = a +X a + X [a + X (a +...)]
Oil 2 3 J
Thus, (a X+a . ) can be generated and then multiplied by X and added
n n^J-
to a o» Tfr6 process is repeated until the final addition of a is
made.
It is obvious that the above algorithm will be valid when the indexing
of the coefficients begins with one rather than zero.
147
-------
SUBROUTINE SOILCO
This subroutine monitors several other subroutines that are used to
simulate the chemical exchange (inorganic) between the soil and
water as water percolates vertically through a soil column. This
subroutine is called from the main or executive program at the com-
pletion of the node balance at the end of each time frame.
The argument list consists of the node index, the write indicator,
the read indicator, and the number of segments in the soil column.
An early version of simulating the exchanges that take place in the
chemistry of waters percolating through a soil column was based upon
a model developed by Dr. Gordon R. Dutt, a staff member of the
University of Arizona. Subsequently, a more elaborate model was
developed at the University of Arizona by Drs. Gordon R. Dutt,
Marvin J. Shaffer, and William J. Moore for the predictions of the
ion exchanges between the applied waters and the soil complex. This
particular model is described in the University of Arizona Technical
Bulletin 196, dated October 1972 and titled "Computer Simulation
Model of Dynamic Bio-Physicochemical Processes in Soils."
A study was made of the Technical Bulletin 196 and the early version
of the soil column simulation was modified using the more recent
techniques as contained in the aforementioned bulletin. In this
148
-------
simulation model the interest is only in the prediction of changes in
salinity. Therefore, the soil column simulation predicts only the
changes in concentration of calcium, magnesium, sodium, sulfate,
carbonate, bicarbonate, and the concentrations of undissociated
calcium and magnesium sulfates. The soil column simulation does not
involve chemical reactions with the nitrate and chloride constituents.
Neglecting the chemical reactions of these two constituents does not
impair the ionic balance.
Most laboratory analyses of the soil extract report concentrations
of the constituents in solution in units of milliequivalents per
liter.
An equivalent weight is defined as the formula weight of a constituent
divided by its valence; i.e., SO/ has a formula weight of 96.06 grams
and has a valence of two; hence, the equivalent weight of SO^
(sulfate) is 96.06/2 = 48.03 grams.
The computational mode of the soil column simulation requires the
concentrations of the constituents be expressed in units of mols
(moles) per liter. A mol is defined as the formula weight of a
constituent per liter of solution. It follows from the above
definitions that valence is equal to equivalents per mol. The
dimensional characteristics of the conversion from milliequivalents
per liter to mols per liter are as follows:
149
-------
_3
Let: milliequivalents = meq = equivalents (10 )
and: (equivalents/liter) (mols/equivalent) = mols/liter
_3
then: (meq/liter/valence) (10 ) = mols/liter
In the solid phase (soil) concentrations are reported in milli-
equivalents per 100 grams of soil.
The concentrations of the solid phase (soil) in the computational
mode are in units of mols per gram. The dimensional characteristics
of this conversion are as follows:
o
Let: milliequivalents/100 grams = meq/grams(10 )
then: meq/grams(102) = equivalents (iQ-5)
grams
and: (equivalents/grams) (10 ) :—- ;—— =
equivalents/mol gram
then: (meq/100 grams) / (valence(10 )) = mols/gram
A very important parameter in the simulation of applied waters
percolating through a soil column is the moisture content of the soil.
It is obvious that this parameter changes for each small time increment
and is dependent upon the variations of applied waters for irriga-
tion, municipal or industrial requirements, precipitation, and the
soil composition. A study of the variations of applied water and
the homogeneity of the soil in the area to be simulated must be
150
-------
made to determine a reasonable value for the moisture content on a
monthly basis. The difficulty in the determination of moisture
content (field conditions) should not be under emphasized, or for
that matter the selection of a representative soil composition for
a large area.
In some river basins the soil column simulation is not required or
has little impact upon the prediction of changes in salinity
(chemistry), and in such circumstances .the soil column simulation
can be eliminated by simply not including any soil data in the
overall data set.
The soil extract in a laboratory analysis is commonly made by adding
one gram of water to one gram of oven-dried soil to produce a one to
one soil-water ratio. The definition of a soil extract represents
the solution separated from a soil suspension, and also by definition
moisture content is the ratio of the weight of water to the weight
of soil. The percentage moisture approximates the (weight of water/
weight of soil) (100).
Field capacity is defined as the moisture content in the soil 2 or
3 days after a thorough wetting of the soil profile by precipitation
or irrigation.
151
-------
It might be well to mention that most results of the analysis
of the chemical composition of waters are reported in units
of parts per million which is a weight to weight unit. To
convert parts per million to milliequivalents per liter, simply
divide parts per million by the equivalent weight of the constituent;
also, it can be easily shown that equivalents per million and
milliequivalents per liter are numerically equal.
The volume of water per soil segment which is part of the input data
is used to compute the moisture content. The bulk density, which
is defined as the ratio of the weight of oven-dried soil to the field
volume of the soil, and the length of the soil segment in centimeters
are also used to compute moisture content. The computation is as
follows:
Moisture content = (Weight of water (grams) /segment)
/ (Bulk density • length of segment
unit area of segment)
Dimensionally:
Moisture content = (grams If^O/segment)/^ grams of soil/cm^)
(cm2) (cni/segment) - Srams
grams soil
The movement of the percolating waters through the soil column is
on the basis of a "piston effect"; i.e., the solution contained in
152
-------
one segment is moved to the next lower segment and then equilibrated
with the soil characteristics of lower segment. A diagram of this
logic is shown on Figure 33.
A diagram with the symbolic names used in the soil column simulation
is shown on Figure 34. The flow chart is shown on Figures 35 through
38, and the table of variable names and arrays is shown in Table 12.
153
-------
LOGIC OF PISTON EFFECT OF SOLUTION MOVING VERTICALLY
THROUGH A SOIL COLUMN (THREE SEGMENTS)
INITIAL
SOIL ANALYSIS
fl APF
J L vi
APPLIED
WATER
^A/~^
(EQUILIBRATE)
BTANK
(MOVE)
(EQUILIBRATE)
(EQUILIBRATE)
BTANK
(MOVE)
(EQUILIBRATE)
BTANK
FIGURE 33
AQUIFER
(MOVE)
154
-------
SYMBOLIC NAMES USED
IN SIMULATION OF SOIL COLUMN
TTANK OR PONDING POOL
SOIL COLUMN
BTANK
GRTANK OR
AQUIFER
FIGURE 34
155
-------
YES
INITIALIZE
COUNTER
FOR WATER
APPLICATION
FLOW CHART FOR
SUBROUTINE SOILCO
ENTER WITH THE ARGUMENTS KNOOE
= NODE INDEX, IW = WRITE UNIT
INDEX, IR = READ UNIT INDEX
AND NSEGS = NUMBER OF SEGMENTS
IN THE SOIL COLUMN
IS THIS THE FIRST CALL
TO THIS SUBROUTINE
NO
BEGIN ITERATIVE
SCHEME TO
EQUILIBRATE ALL
SEGMENTS OF SOIL
COLUMN TO CONVERGE
ON COHC. OF CALCIUM
CALL SUBROUTINE
SETUP TO TRANSFER
DATA FROM ARRAY
"SOIL" TO SET OF
COMPUTATIONAL NAMES
COMPUTE MOISTURE
CONTENT IN PERCENT
IF FIRST CALL TO SOIL
COLUMN USE EXTRACT
RATIO
INITIALIZE THE
EQUILIBRATION
COUNTER
-------
FLOW CHART FOR
SUBROUTINE SOILCO (cont.)
YES
IS THIS THE FIRST CALL
FOR THE SIMULATION OF
A PARTICULAR SOIL COLUMN
NO
CALL SUBROUTINE
FIRSTIM TO COMPUTE
REQUIRED PARAMETERS
FROM THE ANALYSIS OF
THE SOIL EXTRACT
CALL SUBROUTINE
EXCANA TO COMPUTE
CHANGE IN CALCIUM
DUE TO EXCHANGE
WITH SODIUM
CALL SUBROUTINE
EXCAMG TO COMPUTE
CHANGE IN CALCIUM
DUE TO EXCHANGE
WITH MAGNESIUM
CALL SUBROUTINE
GYPEX TO COMPUTE
CHANGE IN CALCIUM
DUE TO STATE OF
GYPSUM
-------
FLOH CHART FOR
SUBROUTINE SOILCO (cont.)
HO
CALL SUBROUTINE
CALCAR TO COMPUTE
CHANGE IN CALCIUM
AS IT REACTS WITH
THE CARBONATES
DO ALL CHANGES
IN CALCIUM CONCENTRATION
(SOLUTION) CONVERGE
TO A SET LIMIT
YES
ALL SEGMENTS HAVE
BEEN EQUILIBRATED
STORE COMPUTATIONAL
SEQUENCE BACK INTO
ARRAY SOIL. ENTRY
RELOAD
NO
MOVE SOLUTION
INTO NEXT LOWER
SEGMENT-PLACE
APPLIED WATER
INTO UPPERMOST
SEGMENT AND
RESET LOOP INDEX
IN ITERATIVE
SCHEME
IS THIS AN INITIAL
EQUILIBRATION OF ALL
SOIL SEGMENTS
YES
FIGURE 37
158
-------
FLOW CHART FOR
SUBROUTINE SOILCO (cont.)
MOVE SOLUTION OF
ALL SEGMENTS INTO
NEXT LOWER SEGMENT.
ENTRY SMOVER. MIX
LOWERMOST SOLUTION
WITH SOLUTION IN
BTANK. REDUCE INDEX
OF ITERATIVE SCHEME
BY ONE
HAS APPLIED WATER PASSED
THROUGH ALL SEGMENTS
EQUILIBRATION IS
COMPLETE FOR APPLIED
WATER-MIX WATER
IN BTANK WITH
WATER IN AQUIFER
UPDATE VOLUME OF
AQUIFER
RETURN
FIGURE 38
159
-------
Table 12
VARIABLE NAMES OR ARRAYS
FOR SUBROUTINE SOILCO
Variable
name
C*
Description
Remarks
or units
BDENS
CALCAR
CA
CATEX
DEPTH
EPSLON
EXCAMG
EXCANA
EXTRAC
FRSTIM
GYPEX
Bulk density of soil =1.33
grams/cm^ (function of assumed
particle density of 2.66 grams/
cm ). Can vary with segment.
Entry to subroutine
Computational aid array to store
concentrations of calcium in
solution as convergent variable
Cation exchange capacity in
units of milliequivalents per
100 grams. Can vary with segment.
Depth of soil segment in
centimeters
Absolute difference of concen-
tration of calcium in solution as
exchanges and reactions converge
to some set limit
Entry to subroutine
Entry to subroutine
Ratio of water to soil as used in
soil analysis (laboratory)
Entry to subroutine
Entry to subroutine
See subroutine
CALCAR
See subroutine
EXCAMG
See subroutine
EXCANA
See subroutine
FRSTIM
See subroutine
GYPEX
160
-------
Table 12 - Continued
Variable
name
C*
Description
Remarks
or units
IFIR
I
ISW
ITER
IW
J
KNODE
K
KV
KWATER
L
MAX
MIN
NCT
NSEGS
PRAT10
RELOAD
Call index to this subroutine
Loop index
Switch used to determine state
of iterative scheme as simula-
tion of waters percolating
through soil column continues
Iteration counter for iterative
scheme used balancing each soil
segment
Write unit designation
Loop index
Node index
Loop index for iterative scheme
Loop limit
Counter for applied waters per-
colated through soil column
Loop index
Loop limit (maximum)
Loop limit (minimum)
Counter used for each
equilibration
Number of segments in soil
column
Moisture content in percent
(gravimetric)
Entry to subroutine
See subroutine
SETUP
161
-------
Table 12
Variable
name
C*
Description
Remarks
or units
SEGVOL
SETUP
SMOVER
WRATIO
Volume of water in segment in
milliliters
Entry to subroutine
Entry to subroutine
Moisture content-ratio of weight
of water to weight of soil
See subroutine
SETUP
See subroutine
SETUP
162
-------
SUBROUTINE SETUP
This subroutine was designed to simplify the transfer of data from
the permanent storage array "SOIL" to a transient set of variable
names used in the computational mode. The subroutine has three
entries including the main entry. The labels assigned to each entry
are indicative of the type function performed.
The argument list consists of the node index, the soil segment
sequence, the number of segments in the soil column, and a switch
index.
The main entry SETUP is used to convert units of input data of the
initial soil analysis to the computational units. This conversion
is only made once for each soil segment of each soil column for
each input; i.e., when switch index is set to zero. Subsequent to the
conversion the data from the "SOIL" array is loaded into the compu-
tational variable name sequence, and from this point on all data in
the "SOIL" array will be in the units required for computation.
A call to the entry SETUP is made from the monitor subroutine
SOILCO as is the case for the other two entries.
Entry SMOVER simulates the movement of the solution down through
the soil segment (piston effect). The first step is the conversion
163
-------
of the concentrations of the constituents in the solution of the
lowermost segments to parts per million and mixed (linear) with the
solution of the symbolically labeled BTANK. The second step is the
transfer of the solution from each of the soil segments to the next
lower soil segment. During the first transfer the applied water is
transferred to the uppermost segment.
Entry RELOAD is called after all segments have been equilibrated
following each move. Data is transferred back into the array "SOIL"
from the computational variable name sequence.
The flow chart for this subroutine is shown on Figure 39. The
variable names or data arrays are shown in Table 13.
164
-------
FLOW CHART FOR
SUBROUTINE SETUP
ENTER WITH ARGUMENT LIST
KN = NODE INDEX, KV = SOIL
SEGMENT INDEX, NSEGS =
NUMBER OF SOIL SEGMENTS
ISW = SWITCH INDEX
I
ENTRY
SETUP
ENTRY
SMOVER
ENTRY
RELOAD
NO
IS SWITCH
SET TO ZERO
YES
CONVERT INPUT
UNITS TO UNITS
FOR COMPUTATION
TRANSFER DATA
FROM SOIL ARRAY
TO COMPUTATIONAL
NAME SEQUENCE
MIX SOLUTION
OF LOWERMOST
SOIL SEGMENT
WITH SOLUTION
IN BTANK.
MOVE SOLUTION
OF EACH SEGMENT
TO NEXT LOWER
SEGMENT
MOVE APPLIED
WATER INTO
UPPERMOST SEGMENT
IF THIS IS FIRST
MOVE
TRANSFER DATA
FROM COMPUTATION
SEQUENCE TO
ARRAY SOIL
RETURN
FIGURE 39
165
-------
Table 13
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE SETUP
Variable
name
BDENS
BTANK
C*
C
C
Description
Bulk density of soil
Data array for symbolic tank
Remarks
or units
See previously
CATEX
EQVWT
EXCA
EXMG
EXNA
EXTRAC
GRTANK
GYPSUM
INTER
ISW
KG
KN
immediately below lowermost
segment of soil column
Cation exchange capacity
Array containing equivalent
weights of constituents
Exchangeable calcium in solid
phase in units of niols per gram
Exchangeable magnesium in solid
phase in units of mols per gram
Exchangeable sodium in solid
phase in units of mols per gram
Ratio of water to soil
Array containing data for
aquifer
Gypsum dihydrate (CaSO, • 2H20)
in units of mols per gram
Intercept for K1 where log K' =
-1.68 • log M - 4.46
Switch index used to denote mode
of simulation
Soil segment index
Node index
defined array
BTANK
Previously
defined
Previously
defined
Previously
defined
See mathematical
appendix
166
-------
Table 13 - Continued
Variable
name
C*
Description
Remarks
or units
KP
KR
K
KV
L
NSEGS
SEGVOL
SLCA
SLCL
SLC03
SLHC03
SLIME
SLNA
SLN03
SLS04
SOIL
Loop index
Loop index
Loop index
Soil segment index
Loop index
Number of soil segments in a
soil column
Volume of water in soil segment
in milliliters
Concentration of calcium in
solution in mols per liter
Concentration of chloride in
solution in mols per liter
Concentration of carbonate in
solution in mols per liter
Concentration of bicarbonate in
solution in mols per liter
Indicator designating presence
of lime in soil
Concentration of sodium in
solution in mols per liter
Concentration of nitrate in
solution in mols per liter
Concentration of sulfate in
solution in mols per liter
Array containing data related to
the simulation of a soil column
See data array
SOIL
167
-------
Table 13 - Continued
Variable
name
C*
Description
Remarks
or units
TTANK
UCAS04
UMGS04
VALNCE
Array containing data for ap-
plied water in ponding pool
Concentration of undissociated
calcium sulfate in mols per
liter
Concentration of undissociated
magnesium sulfate in mols per
liter
Array containing valence of
constituents
See previously
defined data
array TTANK
See data array
VALNCE
ARRAY SOIL (20,10,21)
SOIL (J,K,1) = Concentration of
tion in mols per
SOIL (J,K,2) = Concentration of
solution in mols
SOIL (J,K,3) = Concentration of
tion in mols per
SOIL (J,K,4) = Concentration of
solution in mols
SOIL (J,K,5) = Concentration of
tion in mols per
SOIL (J,K,6) = Concentration of
solution in mols
SOIL (J,K,7) = Concentration of
tion in mols per
SOIL (J,K,8) = Concentration of
tion in mols per
calcium in solu-
liter
magne s ium in
per liter
sodium in solu-
liter
chloride in
per liter
sulfate in solu-
liter
bicarbonate in
per liter
carbonate in solu-
liter
nitrate in solu-
liter
168
-------
Table 13 - Continued
ARRAY SOIL (20,10,21) - continued
SOIL (J,K,9) = Extract ratio from initial soil
analysis
Weight of
water/weight
of soil
SOIL (J,K,10)
SOIL (J,K,11)
Volume of water in soil segment in
units of milliliters
Cation exchange capacity in solid
phase (soil) in units of mols per
gram
SOIL (J,K,12) = Gypsum dihydrate in units of mols
per gram
SOIL (J,K,13) = Lime (CaC03) indicator in solid
phase
SOIL (J,K,14) =
SOIL (J,K,15) =
SOIL (J,K,16) =
SOIL (J,K,17) =
SOIL (J,K,18) =
SOIL (J,K,19) =
SOIL (J,K,20) =
SOIL (J,K,21) =
Bulk density of soil in units of
grams per cubic centimeter
Length of soil segment in centi-
meters
Exchangeable calcium in solid phase
in units of mols per gram
Exchangeable magnesium in solid
phase in units of mols per gram
Exchangeable sodium in solid phase
in units of mols per gram
Undissociated calcium sulfate in
solution in units of mols per liter
Undissociated magnesium sulfate in
solution in units of mols per liter
Intercept for K' (see mathematical
appendix)
1.0=present
0.0=not
present
J = Node index K = soil segment index
169
-------
Table 13 - Continued
ARRAY VALNCE (8)
VALNCE (1) = Valence of calcium = 2
VALNCE (2) = Valence of magnesium = 2
VALNCE (3) = Valence of sodium = 1
VALNCE (4) = Valence of chloride = 1
VALNCE (5) = Valence of sulfate = 2
VALNCE (6) = Valence of bicarbonate = 1
VALNCE (7) = Valence of carbonate = 2
VALNCE (8) = Valence of nitrate = 1
170
-------
SUBROUTINE FRSTIM
This subroutine is called only once per soil segment per soil
column. The overall simulation of the soil column will accommodate
ten segments per soil column.
The argument list consists of the node index, the segment index, and
the write unit designation.
The subroutine is used to compute the concentrations of the various
constituents such as calcium, magnesium, sodium, and the dissociated
calcium sulfate and magnesium sulfate all in the solution phase.
The exchangeable cations of the soil or solid phase are also computed
in this subroutine.
The subroutine is actually the equilibration of each soil segment
using the data from the initial soil analysis (laboratory). The
derivation of the polynomials as used in this subroutine is included
in the mathematical appendix.
The logic used in the subroutine is straightforward; a flow chart is
shown in Figure 40. The variable names or data arrays are shown in
Table 14.
171
-------
FLOW CHART FOR
SUBROUTINE FRSTIM
ENTER WITH KNODE =
NODE INDEX, KSEG =
SEGMENT INDEX AND
IW = WRITE UNIT
DESIGNATION
COMPUTE ACTIVITY COEFFICIENT
FOR DIVALENT ION
COMPUTE CONCENTRATION OF
UNCOMBINED SULFATE
COMPUTE CONCENTRATION OF
UNDISSOCIATED CALCIUM
AND MAGNESIUM SULFATES
COMPUTE ACTIVITY COEFFICENT
FOR MONOVALENT ION
COMPUTE CONCENTRATIONS OF
EXCHANGEABLE CATIONS
(CALCIUM, SODIUM, AND MAGNESIUM)
RETURN
FIGURE 40
172
-------
Table 14
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE FRSTIM
Variable
name
C*
Description
Remarks
or units
BDENS
CA
CATEX
CONST
EXCA
EXMG
EXNA
FS04
GAMMA1
GAMMA2
GAMMA
GSQR
IW
KNODE
K
Bulk density of soil
An array used as a computational
aid in passing coefficients of
polynomials
Cation exchange capacity
Computational aid
Exchangeable calcium in solid
phase
Exchangeable magnesium in solid
phase
Exchangeable sodium in solid
phase
Computation aid that represents
the uncombined sulfate
Activity coefficient for a
monovalent ion
Activity coefficient for a
divalent ion
Entry to Function GAMMA
(Activity coefficient for a
divalent ion) squared
Write unit designation
Node index
Loop index
See function
GAMMA
173
-------
Table 14 - Continued
Variable
name
C*
Description
Remarks
or units
KSEG
MOLE
ROOT
SLCA
SLMG
SLNA
SLS04
TECA
TECMG
TESTA
UCAS04
UMGS04
XMAX
XMIN
Segment index
Six alpha characters used in func-
tion ROOT to designate calling
subroutine
Entry to function ROOT
Concentration of calcium in
solution in mols per liter
Concentration of magnesium in
solution in mols per liter
Concentration of sodium in solu-
tion in mols per liter
Concentration of sulfate in solu-
tion in mols per liter
Thermodynamic equilibrium con-
stant for calcium sulfate. (See
Reference 1.)
Thermodynamic equilibrium con-
stant for magnesium sulfate.(See
Reference 1.)
Convergence term for use in
function ROOT
Concentration of undissociated
calcium sulfate in solution
Concentration of undissociated
magnesium sulfate in solution
Maximum limit of a root
Minimum limit of a root (real) in
a definite range
See previously
defined
subroutine ROOT
174
-------
SUBROUTINE FRSTIM
REFERENCES
1. Garrels, Robert, M. and Christ, Charles, L., "Solutions,
Minerals, and Equilibria," Harper and Row, 1965.
175
-------
SUBROUTINE GYPEX
This subroutine is used to simulate the state of gypsum dihydrate
during the equilibrium process or the exchanges between the solution
and solid phase (soil). Also simulated (computed) are the changes
in concentration of the ion pairs calcium and magnesium sulfates
(undissociated).
The argument list consists of the node index, segment index, the
write unit designation, and the moisture content in percent.
The logic and derivation of the equations used in this subroutine
are included in the mathematical appendix.
The concentration of calcium in solution at the completion of the
computations contained in this subroutine is stored and used as part
of the convergence criteria in the soil column simulation of the
equilibrium process.
The basic functions of this subroutine are indicated by the flow
chart shown on Figure 41. The variable names or data arrays are
shown in Table 15.
176
-------
FLOW CHART FOR
SUBROUTINE GYPEX
ENTER WITH ARGUMENT LIST
KNODE = NODE INDEX, KSEG =
SEGMENT INDEX, IW = WRITE
UNIT DESIGNATION, PRATIO =
MOISTURE CONTENT IN PERCENT
COMPUTE ACTIVITY COEFFICIENT
FOR DIVALENTION AND COEFFICIENTS
FOR POLYNOMIAL AND SOLVE
FOR ROOT OF QUADRATIC
UNDER
IS SOLUTION SUPERSATURATED
SUPER
1
\^ OR UNDERSATURATED ^
>
UPDATE GYPSUM
IN SOLUTION
PHASE
|
— — •"- •*
UPDATE GYPSUM
IN SOLID PHASE
UPDATE CALCIUM AND
SULFATE IN SOLUTION
.
NO
\
IS UNDISSOCIATED
CALCIUM SULFATE AT
MAXIMUM CONCENTRATION
YES
COMPUTE CHANGE IN
CONCENTRATION OF
UNDISSOCIATED
CALCIUM SULFATE
UPDATE CALCIUM
AND SULFATE IN
SOLUTION
COMPUTE CHANGE IN
CONCENTRATION OF
UNDISSOCIATED
MAGNESIUM SULFATE
UPDATE MAGNESIUM
AND SULFATE IN
SOLUTION
±
FIGURE 41
RETURN
177
-------
Table 15
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE GYPEX
Variable
name
C*
Description
Remarks
or units
CA
GAMMA2
GAMMA
GSQR
GYPSOL
GYPSUM
IW
KNODE
KSEG
MOLE
PRAT 10
ROOTWO
SLCA
SLMG
Computational array to aid in
passing coefficients of poly-
nomial
Activity coefficient - divalent
ion
Entry to Function GAMMA
Activity coefficient of divalent
ion squared
Concentration of gypsum in
solution
Concentration of gypsum in solid
phase
Write unit designation
Node index
Segment index
Six alpha character designation
of calling subroutine for
function ROOTWO
Moisture content in percent
(gravimetric)
Entry to function ROOTWO
Concentration of calcium in
solution
Concentration of magnesium in
solution
See Function
GAMMA
See previously
defined function
ROOTWO
178
-------
Table 15- Continued
Variable
name
Description
Remarks
or units
SLS04
UCAS04
UMGS04
W
X
Y
Concentration of sulfate in
solution
Concentration of undissociated
calcium sulfate in solution
Concentration of undissociated
magnesium sulfate in solution
Computational aid
Root of polynomial
Gypsum in solution
179
-------
SUBROUTINE EXCANA
This subroutine simulates the exchange between calcium and sodium
in both the solution and solid phases. The mathematical formulation
of this exchange is included in the mathematical appendix and has as
a basis the Gapon exchange equation.
The argument list consists of the node index, the segment index, the
write unit designation, and the moisture content in percent.
Experience has indicated that the exchange coefficient used in this
portion of the simulation is highly sensitive. The exchange coef-
ficient used in this subroutine is set at 0.707 (see Reference 1);
and on the basis of several applications, evidence is available that
strongly suggests a much higher exchange coefficient is required.
This evidence also suggests that the exchange coefficient is highly
dependent upon the soil type and varies within a wide range from soil
type to soil type.
An indication that the exchange coefficient is incorrect can be
detected by an examination of the net change in the calcium and sodium
concentrations as the equilibrium process passes through this portion
of the simulation. It would be helpful if experimental data were
available for comparison of the predicted and experimental results
to better establish a range of exchange coefficients for each general
type soil. However, the insertion of different values of the
180
-------
exchange coefficient can easily be accomplished within the model
structure without any modification to the computer program.
The concentration of calcium in solution is also the convergence
variable in this portion of the overall soil column simulation.
A flow chart of this subroutine is shown on Figure 42. The variable
names or data arrays are shown in Table 16.
181
-------
FLOW CHART FOR
SUBROUTINE EXCANA
ENTER WITH KNODE =
NODE INDEX, KSEG =
SEGMENT INDEX, IW =
WRITE UNIT DESIGNATION
AND PRATIO = MOISTURE
CONTENT IN PERCENT
(GRAVIMETRIC).
COMPUTE ACTIVITY COEFFICIENT
FOR MONOVALENT ION
COMPUTE COEFFICIENTS FOR
FOURTH DEGREE POLYNOMIAL
SET LIMITS FOR POSITIVE ROOT
AND SOLVE FOR ROOT
IS CONVERGENCE MET
X
s
NO
SET LIMITS FOR
NEGATIVE ROOT
AND SOLVE FOR ROOT
UPDATE CONCENTRATIONS
OF CALCIUM AND
SODIUM BOTH IN SOLUTION
AND SOLID PHASES
NO
RETURN
1
\
PRINT ERROR
MESSAGE AND
HALT COMPUTATION
IS CONVERGENCE MET
YES
FIGURE 42
182
-------
Table 16
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE EXCANA
Variable
name
Description
Remarks
or units
BETA
BETA2
CA
CONST
DEX
DEX2
EXCA
EXNA
GAMMA1
GSQR
ISW
IW
KNODE
K
Conversion factor used to convert
moles per gram to units mols per
liter
BETA as described above squared
Computational array to aid in
passing coefficients of poly-
nomials
A constant used in the argument
list for function ROOT (not
required in this subroutine)
Exchange coefficient
Exchange coefficient DEX squared
Exchangeable calcium in solid
phase
Exchangeable sodium in solid
phase
Activity coefficient for a
monovalent ion
(Activity coefficient GAMMA1)
squared
Computational aid used as a
switching device in changing
limits of desired root
Write unit designation
Node index
Loop index
183
-------
Table 16 - Continued
Variable
name
C*
Description
Remarks
or units
KSEG
MOLE
SLCA
SLNA
TESTA
T
XMAX
XMIN
X
Segment index
Six alpha characters used in
function ROOT to designate
calling subroutine
Concentration of calcium in
solution in mols per liter
Concentration of sodium in
solution in mols per liter
Convergence term for use in
function ROOT
Array used as a computational aid
Maximum limit of a root
Minimum limit of a root (real)
in a definite range
Root of polynomial or change in
calcium in mols per gram
184
-------
SUBROUTINE EXCANA
REFERENCES
1. Dutt, Gordon R., Shaffer, Marvin J., and Moore, William J.,
"Computer Simulation Model of Dynamic Bio-Physicochemical
Processes in Soils," Technical Bulletin 196, Agricultural
Experiment Station, The University of Arizona, October 1972.
185
-------
SUBROUTINE EXCAMG
This subroutine is used to simulate the exchange between calcium
and magnesium in both the solution and solid phases. The mathe-
matical formulation of this exchange is included in the mathematical
appendix.
The argument list consists of the node index, the segment index, the
write unit designation, and the moisture content in percent (gravimetric),
The exchange coefficient which has a numerical value of 0.67 (see
Reference 1) is generally accepted and has not presented any problems
in the application of this portion of the simulation.
The variable names or arrays are shown in Table 17.
The concentration of calcium in solution at the completion of the
computations contained in this subroutine is stored and used as part
of the convergence criteria in the soil column simulation of the
equilibrium process.
186
-------
Table 17
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE EXCAMG
Variable
name
C*
Description
Remarks
or units
BETA
CA
DEX
EXCA
EXMG
KNODE
KSEG
MOLE
PRAT 10
ROOTWO
SLCA
SLMG
X
Conversion factor used to convert
mols per gram to units of mols
per liter
Computational array to aid in
passing coefficients of
polynomials
Exchange coefficient
Exchangeable calcium in solid
phase
Exchangeable magnesium in solid
phase
Node index
Segment index
Six alpha characters used in
function ROOTWO to designate
calling subroutine
Moisture content
Entry to function ROOTWO
Concentration of calcium in
solution in mols per liter
Concentration of magnesium in
solution in mols per liter
Root of polynomial or change in
calcium and magnesium in mols
per gram
See previously
defined ROOTWO
187
-------
SUBROUTINE CALCAR
This subroutine is used to simulate the reaction of calcium with
the carbonates. The mathematical formulation of these reactions is
included in the mathematical appendix.
The argument list consists of the node index, the segment index, the
write unit designation, and the moisture content in percent
(gravimetric).
The initial input of data for the soil column simulation contains an
indicator signifying the presence or lack of presence of lime (CaCO )
in the soil or solid phase. Then intercept for the relationship of
K* is either computed or set as a constant dependent upon the lime
condition. The significance of K" is explained in the mathematical
appendix. The intercept is only computed once per segment per soil
column per run.
The concentration of calcium in solution at the completion of the
computations contained in this subroutine is stored and used as part
of the convergence criteria in the soil column simulation of the
equilibrium process.
The flow chart for this subroutine is shown on Figure 43, and the
variable names or arrays are shown in Table 18.
188
-------
FLOW CHART FOR
SUBROUTINE CALCAR
ENTER WITH ARGUMENT
LIST KNODE = NODE INDEX,
KSEG = SEGMENT INDEX,
IW = WRITE UNIT, PRATIO =
MOISTURE CONTENT IN
PERCENT
NO
IS LIME PRESENT
YES
SET INTERCEPT
TO A CONSTANT
NO ^""HAS INTERCEPT BEEN
COMPUTED
YES
COMPUTE INTERCEPT
COMPUTE CHANGE
IN CALCIUM OR
CALCUIM BICARBONATE
CONCENTRATION
UPDATE CONCENTRATIONS
OF CALCIUM AND CALCIUM
BICARBONATE IN SOLUTION
RETURN
FIGURE 43
189
-------
Table 18
VARIABLE NAMES OR DATA ARRAYS
FOR SUBROUTINE CALCAR
Variable
name
C*
Description
Remarks
or units
CA
CONST
GAMMA1
GAMMA2
GAMMA
INTER
ISW
IW
KNODE
KPRIME
K
KSEG
Computational array to aid in
passing coefficients of poly-
nomials
A constant used in the argument
list for function ROOT (not
required in this subroutine)
Activity coefficient for a
monovalent ion
Activity coefficient for a
divalent ion
Entry to function GAMMA
Intercept used in equation
for K1
Computational aid used as a
switching device in changing
limits of desired root
Write unit designation
Node index
Variable name for K1
Loop index
Segment index
See function
GAMMA
190
-------
Table 18 - Continued
Variable
name
C*
Description
Remarks
or units
MOLE
PRAT10
ROOT
SLCA
SLHC03
SLIME
TEMP
TESTA
XMAX
XMIN
Six alpha characters used in
function ROOT to designate
calling subroutine
Moisture content in percent
(gravimetric)
Entry to function ROOT
Concentration of calcium in
solution in mols per liter
Concentration of sodium bicar-
bonate in solution in mols per
liter
Indicator to designate presence
of lime in soil
Computational aid
Convergence term for use in
function ROOT
Maximum limit of root
Minimum limit of root (real) in a
definite range
Change in concentration of calcium
and bicarbonate in solution in
units of mols per liter
191
-------
SUBROUTINE CALCAR
REFERENCES
1. Dutt, Gordon, R., "Prediction of the Concentration of Solutes in
Soil Solutions for Soil Systems Containing Gypsum and Exchangeable
Ca and Mg," Soil Science Society of America, Proceedings,
Volume 26, Number 4, July-August 1962, pages 341-343.
192
-------
SUBROUTINE SCRIBE
This routine is used to output the values of pertinent parameters
of the soil-water equilibration.
The argument list consists of on the write unit designation.
The units of output are those in the computational mode. This
routine is an excellent diagnostic tool in the overall soil column
simulation.
The variable names have all been previously defined.
193
-------
FUNCTION GAMMA
This subroutine is used to compute the activity coefficient of a mono
or divalent ion. The mathematical formulation for this subroutine
is included in the mathematical appendix.
The argument "KVAL" represents the valence of ion.
194
-------
Ifethematical Appendix of Soil Column Simulation
Activity coefficient - Gamma (7)
Debye-Huckel theory states that:
r f«\ -K?&~
Log (r) « =
10 1 + p -JOT
where:
\ and p are temperature constants and at 27° c or k9° F,
\ = 0.509, P • 1.00
u SB ionic strength
Z a valence of ion species in solution
Substitution of constants for \ and p, and changing from Log.^ to Loge,
where Loge (10) = 2.30259, then:
** w - -~*^
taking logarithms of both sides -
- 1.17201 Z2 •v/u"
1+^
7 = I
i=n
where: ionic strength = u «= 1/2
j
where: C. = concentration of i , ion species in solution in units
of moles per liter
n ss number of ion species in solution
195
-------
Computation of concentration of undissociated calcium
sulfate (CaSO? in units of moles per liter in terms of GS
Begin derivation of an equation with one -unknown. Let unknown equal
concentration of the uncombined sulfate -
where: Cgo = concentration of sulfate in soil extract in units
T of moles per liter (total)
C,,- a concentration of uncombined sulfate (sulfate
k \
ion) in units of moles per liter
cr. an = concentration of sulfate that combines with
L-abOj,
calcium to form CaSOjL in solution in units
of moles per liter
concentration of sulfate that combines with
magnesium to form MgSo£ in solution in units
of moles per liter
In the use of equation (l) it becomes necessary to find expressions
for C »* C in temS of C
Solving for Cnr%an in terms of C :
k
Let - C_ = C + c
Ca
where: GC = concentration of calcium in soil extract in units
of moles per liter
CCa = concerrtra*iori o:f uncombined calcium (calcium ion)
in units of moles per liter
196
-------
CCaSO = concen'fcra'fcion of calcium that combines with sulfate
to form CaSO^ in solution in units of moles per liter
Given the thermodynamic equilibrium constant (K_ _A ) for CaSO,° as
a a CaSOj, H-
Ca * SO*
K°aSOr aCaS04 (2)
where:
a_ = activity of calcium ion
aSO ~ ac*ivi'fcy of sulfate ion
an so = ac^v^y °^ calcium sulfate (ion pair)
and
aCa = °Ca * 7Ca
ASSme
Substitution into equation (2) then
CCa ' 7Ca ' CSO,,
Solving for €„ using equation (3)
K" P
CaSOjk * CaSOu
n —i —
Ca "
n .0- .
Ca SOi
- .
Oi bO^
197
-------
then: KCaSO^ * CCaS0lf
or:
* 7S0lf '
7Ca
and: C
' 7Ca *
7Ca * SOu *
experimentally KCagQ = ^.9 • 10
The expression then for C g_ in terms of CSQ
is
rca
Computation of concentration of undissociated magnesi\jm stxlfate
in units of moles per liter in terms of C^
Solving for C in terms of C -
vhere: C = concentration of magnesixan in units of moles
per liter (total)
C = concentration of uncombined magnesium
(magnesium ion) in units of moles per liter
198
-------
concentration of magnesium that combines
with sulfate to form MgSQ^ in solution in
units of moles per liter
Given the thermodynamic equilibrium constant (K., — ) for MgSOj, as -
' ^ (5)
where: a^, » activity of magnesium ion
- activity of sulfate ion
a.. qn = activity of magnesium sulfate (ion pair)
and:
Assume 7^^ = LOO
substitution into equation (5) then -
SfeSO^
Solving for CM using equation (6)
199
-------
then:
and* P
Experimentally K_ = 5.9 . 10
= . .
The expression then for C,, „_ in terms of C__ is
CMg ' 7Mg ' 7SO, * CSO,
r -
^ _
4 5'9 ' 10
+ 7Mg •
Equations (k) and (?) are expressions of the undissociated calcium
sulfate and magnesium sulfate in terms of C . It now becomes
necessary to solve for C , which is the concentration of the
uncombined sulfate (sulfate ion) in units of moles per liter
substitution of equations (k) and (7) into equation (l)
200
-------
• 7Ca •
'SOijg, ~SO^ • M ^ 1Q-3 + y^ f 7^ ^ c^ ' 5.9 . ^-5
•y =7 =7 — ?„
7Ca 'SO^ 'ife 72
where: 72 = activity coefficient for a divalent ion.
Then:
/Ca.p-^'X __C- '7|
io"5 + yl ^~C~n 5.9
or:
A —
k.9 • IO" + 7l . C^ 5.9
(5.9 -10- + 7
expanding:
201
and
0 - C0rt (^.9 • 10"5 + 722 • C ) (5.9 • Hf5 + i\ •
- ' V5'9
-------
0 = 4.9 • 10-3 . 5.9 • ID"3 . C + 5.9 • ID'3 y\
-------
x3 + r2, I" (5.9 • io"3 + 4.9 . Hf3) + 72 (C + c „ c ) 1 ^
H ' 'o » Pa M«. ^Ort / A
+ T 4.9 • 10"3.5.9 • IO"3 + yz (5.9 • IO"3 Cn + 4.9
L
or
o . r + r I.OB .
The real root of the above polynomial (equation (8)) must be
0 5 X < CSQ
The solution of equation (8) will result in the concentration of
the uncombined sulfate ion SO^ in units of moles per liter. With
the value of Cqn known equations (4) and (7) can be used to determine
the concentrations of the undissociated calcium sulfate (CaSp£) and
magnesium sulfate
Exchangeable cation concentrations in the solid phase
The total exchangeable cation concentrations in the solid phase are
represented by the cation exchange capacity. This capacity is given
as part of the initial soil analysis in units of milliequivalents per
203
+ I 2.89 • IO"5 + -f (5.9 • IO"3 C_ + 4.9 • IO"3 CM - (8)
L 2 C^ %T
1.08 . IO"2 C. ) 1 X - 2.89 . IO"5 C
-------
100 grams and is converted to equivalents per gram for computational
purposes. Equation (9) shows the sum of the exchangeable cations
in the solid phase.
jj _ N + jj 4. u (a,)
where: N = Total cation exchange capacity in units or
equivalents per gram
w = Concentration of sodium ion in the soil or
Na
solid phase in units of equivalents per gram
Kp = Concentration of calcium ion in the soil or solid
phase in units of equivalents per gram
In equation (9) there are three unknowns. With the use of widely
used expressions the concentrations of magnesium and sodium will be shown
in terms of the concentration of calcium. An equation that is widely
used for the calcium-magnesium exchange is shown below as equation (10):
a,, H
fa PQ
\JO» __. Od.
where: a = the activity of the calcium ion in solution
= the activity of the magesium ion in solution
N_ and N,. have been defined above
Ca Mg
KI = 0.67 (experimentally derived)
204
-------
In a previous portion of this appendix it was shown that
aCa = CCa ' 72
._ .
Mg 2
where j^ represents activity coefficient for divalent ions.
Upon substitution into equation (10)
aCa _ NCa C'72
or:
Equation (ll) is an expression for the concentration of the
magnesium ion in the solid phase in terms of the concentration of
the calcium ion. With the use of the Gapon relationship of sodium
and calcium an equation can be derived to represent the concentration
of the sodium ion in the solid phase in terms of the concentration of
the calcium ion.
The Gapon relationship is as follows:
^a _ Na /12\
'/ - K2 if^ ^ '
where: a = activity of the calcium ion in solution
= CCa * I
a. = activity of sodium ion in solution
Ka
— r .7
- °Na * '1
205
-------
where 7 = activity coefficient for the monovalent ion
i
Upon substitution into equation (l2)
aNa
aCa Ca
or:
N . C.T . 7,
ca Na
K =2.0 (experimentally derived)
then:
_ NCa • CNa ' 7z
NNa -—(15)
Equation (15) is an expression for the concentration of the magnesium
ion in the solid phase in terms of the calcium ion in the solid phase.
Upon substitution of equations (ll) and (33) into equation (9) as
follows:
N_ . C.T .7 0.67 . w . CM
f^Q 1\TS3 1 I*Q M*T
TO _ ua j.Ma x M (^a jxig
2 o Jp - y — Ca C^
*° V * 7 Ca
Ca
or:
°Na * 7i 0.67 . n
+ 1.0 + -~ - JSL
2.0
NT
„ III • -J-.V • ~ n
2 O »/(^ -y t»
^. U \J\sn . /
' CMg
Ca
206
-------
All the values on the right side of equation (l^) can be obtained
from the data of the soil extract analysis; therefore, the concentra-
tion of the calcium ion in the solid phase can be determined. After
the determination of the concentration of the calcium ion equations (ll)
and (13) can be used to determine the concentrations of the magnesium and
sodium ions, respectively, as they exist in the solid state.
Solubility of CaCCU (lime) in water
In the University of Arizona Technical Bulletin 196, dated October 19?2
and titled "Computer Simulation Model of Dynamic Bio-Physicochemical
Processes in Soils," the dissociation (solubility) of CaCO, is shown
as:
CaCO, -»- Ca + CO,"
and the thermodynamic solubility product Ksp is given as:
Ksp = a_ . arr.
The concentration of CO, is a function of the COg partial pressure
and HC07 is usually the predominant form in which C02 occurs in a
soil-water system. Thus it is more convenient to consider the
following reaction:
HgCO, + CaCO, t Ca + 2R COj
where the thermodynamic equilibrium constant K^ CQ is expressed
as follows:
207
-------
aCa ' a HC
Oj
Ca *
2 *
HCOj *
72
i
where 7 = activity coefficient of unchanged species = 1.0.
_
Let . = C . C = — (15)
r2.r^ Ca HC05 r2.7^
Due to variations of C02 and the solubility product Ksp a relation-
ship vas developed from experimental data for the value Kf. This
relationship is shown in Technical Bulletin 196, page 13, equation (37)
as:
Log K« = -1.68 • Log M - ^M (l6)
where M = Moisture content of soil in percent.
With equations (15) and (l6) the changes in CaCO, solubility with
changing moisture contents in the soils can "be computed.
Let "X" = change in concentration of the divalent ion in solution
in units of moles per liter
Then: CL = C_ + X
Ca Ca
CL™ - C + 2X
HCOj HCOj
Substitution into equation (15) -
208
-------
' CCa '
or
-7.050206
X
= CCa ' + CCa ' C X
Collecting terms:
3 X X K
1 + N/U~
and: 0 = ^.0 X3 + ^.Oi CTT^4C_ ) X2-^ ^.0 CL CTT^ 4C2 )x-
Determination of change of gypsum from solution^to solid or from
solid to solution
The chemical formula for gyps-urn is CaSO^ . 2HgO and the reversible
reaction is as follows:
209
-------
C + sok* + 2H2°
Gyps-am can be present in the solid phase or can be precipitated from
a supersaturated solution. Gypsum is defined as a moderately soluble
salt whose solubility product "Ksp" is equal to 2.4 . 10 . The
solubility product is defined as follows:
KSP " aCa '
where: an = activity of calcium ion
l/a
s= activity of sulfate ion
From previous definitions
aCa - CCa ' \
Substitution into equation (l?)
Ksp . CCa . 72 . Cs^ . 72 (18)
To reiterate:
C = concentration of calcium ion in solution
Ca
in units of moles per liter
Cc_ <= concentration of sulfate ion in solution
bO^
in units of moles per liter
7 = activity coefficient of divalent ion
7
2
N^T
210
-------
Let "X" « change in concentration in calcium and sulfate ions.
If: X > 0 - solution is undersaturated
X = 0 - solution is saturated
X < 0 - solution is supersaturated
Upon substitution equation (l8) becomes:
Ksp . (Cca + X) (Cs + X) . Jl
or:
o = x2 + (cca + c_ ) x -
and:
^ Ca SO^ \
X -
2.0
Determination of change in concentration of ion pairs CaSOj^ and
Equation (3) defines the thermodynamic equilibrium constant - Kg
as follows:
2
CCa * 7Ca * CS01, * 7SOi, °Ca * CSOi, * 72
* (19)
and: K_ SQ = ^.9 ' 10" (experimentally derived)
Let: X = change of concentrations of calcium, sulfate, and calcium
sulfate in solution in units of moles per liter. An
increase in concentration in calcium and sulfate represents a
decrease in the concentration of calcium sulfate
211
-------
then equation (19) can be expressed as follows with the inclusion of "X
(C + X) ff- + TC) • y2
2 Ca
^.9 • 10 = rr—
- X)
or:
expanding:
0 = r2 Y2 + I y2 (r + c } + ^ Q * 10"^ I x + r2 C
70 X + | 7^ ^CCa + CgpJ + .y J-0 A + 72 ^a
Let:
" (CCa
CCa '
2a
It is easily shown that the maxinrum concentration of CaSOj, is a function
of the two experimental constants Ksp and 1C, OA as follows:
= CCa
, „ °Ca • C
and:
then: Ks 2 k in~^ 5
(max) = ——*— = -Z—* • • v = ^.9 • 10~ moles per liter
CaSOj^ ^.9 * 10
212
-------
Equation (6) defines the thenaodynamic equilibrium constant, K
as follows:
and: ^SOu = *^ * 10 (experimentally derived)
Again let: X = change of concentration of magnesium, sulfate, and
magnesium sulfate in solution in units of moles per
liter. An increase in concentration in magnesium and
sulfate represents a decrease in the concentration
of magnesium sulfate.
o
Using the same algebra as was used for the change in CaSO^, then:
0— -v Y^ a. -v f n 4. n ^ 4. "5 o • i c\* Y j. "/ p _ *5 o * i n"
- /„ A + /„ ^U»y • -LU
Let:
10"5
C ss
and
X =
Determination of calcium-sodium exchange
Equation (15) defines the relationship of calcium and sodium ions
as follows (Gapon):
213
-------
CHa
Squaring both sides of the above equation
JCa wMa '3. (21)
•v*^ /i . -v
\ CCa 72
Reiterating for terms in equation (2l)
!LT = concentration of sodium ion in solid phase
TSfa
(soil) in units of moles per gram
N_ = concentration of calcium ion in solid phase
oa
(soil) in units of moles per gram
C.7 = concentration of sodium ion in solution in units
Na
of moles per liter
C m concentration of calcium ion in solution in
units of moles per liter
7 and 7 = activity coefficients of monovalent and
divalent ions, respectively
214
-------
Let: X * change in concentration of calcium ion in the solid
phase in units of moles per gram
then if X > 0:
Concentration of calcium ion in solid phase decreases
Concentration of calcium ion in solution increases
Concentration of sodium ion in solid phase increases
Concentration of sodium ion in solution decreases
and in terms of "X" then:
NCa - NCa
CCa - CCa
CHa -
where:
p a conversion factor to convert moles per gram to
moles per liter and must have units of grams
per liter.
215
-------
Substitution into equation (22);
-------
The general equation for the above polynomial is:
where:
A . -
B * NCa
C .
a NCa CNa (Sa
E = 72 MJ; C_ - K2 N^ C5T
i Na Ca 3 Ca Na
Determination of Calciimi-Magnesitim Exchange
Equation (10) defines the relationship of calcium and magnesium as follows:
where:
-
_ _ Ca
a_ = the activity of the calcium ion in solution
Ca
the activity of the magnesium ion in solution
N = the concentration of the calcium ion in solid
Ca
(soil) phase in units of moles per gram
N = the concentration of the magnesium ion in solid
phase in units of moles per gram
217
-------
and: a_ = C_ • 7
Ca Ca 2
where:
Cn = the concentration of the calcium ion in
Ca
solution in units of moles per liter
C a the concentration of the magnesium ion in
solution in units of moles per liter
7 = activity coefficient of divalent ion
or:
7 N
2 _Ca /22)
m
Let: X = change in concentration of ions in solid phase
and: CCa=CCa+PX
c = c^ - ^x
N- = N_ - X
Ca Ca
f X
Then if: X > 0
Concentration of calcium ion in solution increases
Concentration of calcium ion in solid decreases
Concentration of magnesium ion in solution decreases
Concentration of magnesium ion in solid increases
218
-------
Substitution into equation (22)
"* 4» AY "V W V
^Ca * ^ 72 _ NCa " X
sis _K
-1" CCa»* -
Let: A = p (l - K )
or:
B =
C = C_ NM - K CM N_
Ca Mg i Mg Ca
„ - B + ^ B2 -
X B
or:
\C_ T pX/ \H»« "^ Xj = K \Q%, " &X) (N_ •• X)
v>a Mg i r^g Ca
expanding:
C-, K.. + (C_ 4- PH.. ) X + PA = K C., N^ - ((
Ca Mg Ca Mg i Mg Ca
Collecting terms:
P v-*- ** ^ / X^ "•" I p vli-.. "•" IV U— J "r K Cm. "r C^ I A - Ak. w,, *!„
\ . i K \ «„ n« / M_ Ca i Mg Ca
219
-------
COMPUTER PROGRAM LISTING FOR SIMULATION MODELING
PROGRAM ACUMEN
roooo THIS COMPUTER APPLICATION HAS BEEN DESIGNED AS A SYSTEMANALYSIS
roo»o AND CONJUNCTIVE USE MODEL EMULATING A NODAL CONCEPT IN DIVIDING
feooo A WATER BASIN INTO PERTINFNT FUNCTIONAL OR PESPONMVE CRITICAL
C«o«o AREAS ALLOWING AS MUCH FLEXIBILITY IN MODEL DESIGN AS WAS
roooo CONSIDERED PRACTICAL
COMMON / GFNPL / LTlTLE(B), NCTRLI20, 3), MQ*EG, IYRBEG, MOEND, I
AYRENC, NUMNnD, ISECM20. 20), KpRI, NUMSEQ(20M NDWTR(20, 20)
COMMON / RFSE" / RDATA(20, 110), RULF<20, 13', Cf>NFVP(20, 12), BA
ANKCO(?0), KSURNM(20, fi), OINFLOI20, 10. 12), DEMAMD120, 10, 16), *
PDEMNM(20, 1", 8), KQINNM(20, 10, 6) , KAQIjNM (20,P)
COMMON / PROPER / SOlL<20, 10, 21), RTANK(20, 10), TTANK120, 10),
1 GRTANK(20, 10) , MSEC (20) » CONUSE (20, 10, 30)
? , KUSE (20. 10, 30)
COMMON / POWER / CAPBAB(?o» 2o), TDATA(2o, A), EMERGY(2o, 8), PDA
ATA120, 12)
COMMON / SALT / QTOT120, 20), WQAL<20. 20)
COMMON / PAPAM / REGRES ( 20, 11 ) , PEGWR ( '0, 10 ) ,
1 PREPIC ( 20, 10 ) , OBSVAL (20,10), CHEMDM ( 20, 10, 10)
DIMENSION OlNRlV (10), NPASS (20), cRSERV (10), *FA (20) ,
i OFLOW no ) , TEMP no)
EQUIVALENCE ( CONUSE, KUSF )
DATA ( IFIR = 0 )
IF ( IFIR . NE. 0 ) Go To 10
ceooo THIS IS THE FIRST TIME PROGRAM HAS BEEN FNTfR^n FOR A PARTICULAR
r«o<>« JOB - SET READ AND WRITE UNITS AND READ SYSTEM STRUCTURE AND
PERTINENT U'lTI*LI 7 AT ION PARAMETERS
IR = 60
IW = 61
LINE = 0
DO 5 tf = 1, 4
JREAD = K
CALL PERUSE ( JPEAD , Iw , IR )
CONTINUE
IFIR = i
MO = MOBFG
IYR = IYRBEP
INITIAILIZE PASS INDEX OF EACH NODE FOR USE IN ITERATIVE SCHEME
10 DO 15 K = 1, 20
NPASS ( K ) = 0
15 CONTINUE
INITIALIZE STARTING INDICES FOR NODE (MIN) AND SEQUENCE ( MINA )
>«* LOOPS
MIN = 1
WRITE ( IW , 16 )
16 FORMAT ( 1H1, 8X )
BEGIN ITERATIVE SCHEME FOP EACH NODE WHERE NUMNOD = NUMBER OF NODES
20 CO 600 KN = MIN , NUMNoH
CHECK TO DETERMINE IF THIS IS THL FIRST PASS THROIIGH THIS NODE
220
-------
PROGRAM ACUMEN
C o o « o
25
50
C o o « o
60
Of>
70
80
10
.^
) = 0.0
) = 0.0
) = GRTANK
LET NS =
THIS IS THE FIRST PASS THROUGH THIS NODE »» INITIALIZE RIVER
FLOW AND CHFMISTRY ARRAY ( QINRIV ) AND UPDAT^ AQl'IFER ASSUME ONE
MONTH LAG IN CHFMISTRY MIX «« LET ARRAY PFGWR = BFGINNING"
CONDITIONS nF AQUIFER AT START OF TIKE FRAME
CO 2$ K = 1,
QlNRIy (K) =
BTANK ( KN, K
TTANK ( KN, K
BEGWR (KN, r ) = GRTANK ( KN, K )
CONTINUE
THIS PORTIOW OF THE SUBROUTINE IS USEH Pf)R ALL PASSFS
NUMBER OF SEQUENTIAL OPFRATIONS IN THJR NODE
NS = NUMSEQ ( KM)
LOOP THROUGH EACH SEQUENCE IN NODE
DO 500 KS = MINA , NS
IS = TAB? ( ISEO (KN, KS ) )
!«=• ( ISEQ (IfN, KS)) 50. 350, 400
CONST = DEMAND ( KN, Ifi, MO )
TEST = DEMAND < KN, is, 1*1
IF ( TEST .r-T. 2.0. AND. TEST . LT. 6.0 ) Go TO 2*0
THIS OPERATIONAL SEQUENCE REPRESENTS A DEMAND IN THE NODE
kHOSE SOURCp MAY BE FROM THE RIVER ( SURFACE PLOW ) OR FROM THE
AQUIFER ( SUBSURFACE FLOW ) TO SUPPLY EITHER AN IRRIGATION
CR A MUNICIPAL AND INDUSTRIAL REQUIREMENT
IF { pEMANp .GT. 1.0 ) GO TO i^O
CEMAND IS TH BE SUPPLIED FROM RIVER < SURFACE FLOu)
IF ! TEST . EG. 6.0 ) GO TO 230
TEMPR = QINPIV (1) - CONST
IF ( ABS (TFMPR).LT. 1.0 ) TEMPP = 0.0
IF ( TEMPR . LT. 0.0 ) GO TO 9o
DEMAND CAN RE FuLLY MET BY FLOW IN THE RlvER SET INDEX (INDY) = 1
INDY = 1
CALL RETURN FLOW AND SHORTAGE SUBROUTINE ( RFCOMP > TO COMPUTE
AND ALLOCATE RETURN FLOW MAGNITUDE AND CHEMISTRY AND COMPARE
SHORTAGE ( IF ANY ) WITH SHORTAGE INDEX *><> IF THE TEST FAILS THE
SHORTAGE CRITERIA RFCOMP WILL TERMlNftTE COMPUTATION
CALL RpCOMP ( QINPIV, KN, IS, INTY, MO, CONST, IW )
SET ACTUAL "IVEPSION AND CHEMISTRY OF WATERS DIVERTED
DEMAND < KN, IS, 16
QINRIV (1) = TEwpR
10
IS
) = CONST
K ) = QINRIV (K)
90
DO 80 K = 2,
CHEMDM ( KN,
CONTINUE
GO TO 500
DFMAND CAN »'CT BE FULLY SUPPLIED BY FLOW TN THE RTVER SEARCH
UPSTREAM RESERVOIRS FOR POSSIBLE RELEASES - MAKE RELEASES FROM
RESERVOIRS ON THE BASIS OF" AN IDEAL RANKING PROCESS-
CALL DTABLE f KN, KFA, 1C, KS )
CHECK TO DETERMINE IF RESERVOIRS EXIST ABOVE DEMAND POINT
221
-------
PROGRAM ACUMEN
THERE ARE NP UPSTREAM RESERVOIRS AVAILABLE FOR ADDITIONAL RELEASES
SET DIVERSION EQUAL TO FLOW IN RIVER AND CALL ROUTINE RFCOMP
r»«» AND ALLOCATION. IF SHORTAGE is NOT TOLERABLE ROUTTNE RFCOMP WILL
TERMINATE COMPUTATION
100 INDY = 2
CONST = QINPIV (1)
TEMPR = 0.0
GO TO 60
SEARCH FOP AVAILABLE UPSTREAM RFSERVCIRS FOR POSSTBLF ADDITIONAL
RELEASES TO MEET THIS DEFICIT
110 DO 140 K = 1, 1C
KFAK = KFA
KNODE = NFI»'D f KFAK )
CHECK TO DETERMINE IF THIS RESERVOIR IS CONSTRAINED IN MAKING
C<> LET PRELM = TOTAL RELEASE MADE DURING PREVIOUS PASS OR CALL
PRELM = RDATA ( KNODE «10 )
TEMPR = RELEASE REQUIRE" T0 MEET DEFICIT
QFLOW (1) = TEMPR
CALL RESOPR f QFLOW , KM, MO, IW )
LET DFLTAF = ADDITIONAL RELEASE (IF ANY )
DELTAF = RDATA ( KNODE i in ) - pPELM
IF < DEl-TAF • GT. 0.0 ) GO TO 120
NPASS ( KNOPE ) = 100
GO TO 140
ADDITIONAL RELEASE CAN RE MADE #» SFT LOOP INDICES FOR RE-ENTRY
X5<>« jNTO ITERATIVE SCHEME
1?0 TEMPR = TEMPR + DELTAF
NPASS ( KNO"E ) = 1
IF ( TEMPR . LT. 0.0 ) NPASS ( KNODE ) = 100
UPDATE VOLUME A"D CHEMISTRY OF RIVER FLOW BELOw RFSERVOlR
DO 1Z5 KA = 2, 10
QINRIV ( KA ) = QFLOW (KA)
125 CONTINUE
QINRIV (1) = RDATA ( KNODE » 10 )
ro«»o RESET LOOP INDICES
NS = NUMSEQ (KNODE )
DO 130 KA = 1, MS
KX = KA
IF ( ISEQ ( KNODE i |(A ) . EQ. 0 ) GO TO 135
130 CONTINUE
GO TO 1040
NOTE «<>«« A RESERVOIR CANNOT BE THE LAST SEQUFNCF QF OPERATIONS
jN A NODE **»«
IF ( KX . EO. NS ) GO TO 1050
MIN = KNODE
MINA = KX •»• 1
222
-------
PROGRAM ACuMEN
HO CONTINUE
» ALL RESERVOIRS HAVE BEEN SEARCHED AND NO ADDITIONAL RELFASES CAN
BE MADE TO MEFT DEFICIT "
GO TO 100
C«««« DEMAND IS TO BE SUPPLIED FROM AQUIFER ( SUBSU&FACF FLOW )
170 TEMPR = GRT»NK ( KN, 1 ) - CONST
IF ( TEST .FQ. 6.0 ) GO TO 340
IF ( ABS < TEMPR j .LT. i.o ) TEPPR = o.o
IF ( TEMPR . LT. 0.0 ) CO TO 220
DEMAND CAN BE FULLY SUPPLIED FROM AQUIFER ( SUBSURFACE FLOW )
SET INDEX (INBY> = 1
INDY = 1
CALL RETURN FLOW AND SHORTAGE SUBROUTINE ( RFCOMp ) TO COMPUTE
(•»«»« AND ALLOCATE RETURN FLOW MAGNITUDE ANf CHEMISTRY AND COMPARE
C««»» SHORTAGE ( IF ANY > WITH SHORTAGE INDEX ** IF THE TEST FAILS THE
f SHORTAGE CRITERIA RFCOMP WILL TERMINATE COMPUTATION
180 DO 185 K. = ?, 10
TEMP (K) = PERWR { KN» K »
185 CONTINUE
CALL RFCOMP t TEMP , KN, is, INDY, MO, CONST, iw i
GRTANK ( KN, 1 ) = TEMPR
C«o»« SFT ACTUAL DIVERSION AND CHEMISTRY OF WATFpS DlVEoTEn
190 DEMAND < KN , jS, 16 ) = CONST
CO 200 * = ?, 10
CHEMDM ( KN, is, K ) = BEGwR ( KN, K )
200 CONTINUE
GO TO 500
ro««« DEMAND CAN NOT gE FULLY SUPPLIED BY AQUIFER- rOMPMTE SHORTAGE
C«««o RETURN FLOW AND ALLOCATION ( IF SHORTAGE NOT TfUEPABi E RFCOMP
C»e»» WILL TERMINATE COMPUTATION )
?20 CONST = GRTANK ( KN, 1 )
INDY = 2
TEMPR = 0.0
GO TO 190
(•«««« THIS DEMAND REPRESENTS A SURFACE FLOW THAT IS TRANSFERRED OUT
C«««« OF SYSTEM
?30 TEMPR = QINPIV (1) - CONST
IF ( TEMPR . LT. 0.0 ) Go TO 1010
GO TO 70
<•««<»« NHEN DEMAND ( KN, IS. 14! THE DESTINATION OR USE INDEX IS GREATER
<•«««« THAN 2.0 SIMPLE TRANSFERS OF FLOW ARE TREATED AS DEMANDS SO THAT
C««»o ROUTINE RFCOMP WILL SET Up THE PROPER QINFLO ARRAY. THESE TRANSFERS
r«««« CAN HAVE RECOMPUTED VOLUMES (INPUT DATA > OR WILL BE COMPUTED IN
r««»s SIMULATION
?50 IF ( DEMAND ( KN, IS, 13 > -GT. 1.0 ) GO TO 290
t«« TRANSFERS WILL BE MADE FROM THE RIVFR (SURFACE FLOW) TO AQUIFER
r
-------
pRQGRAM ACUMEN
IF ( ABS (TFMPR) .LT. 1.0 ) TEMPR = 0.0
IF ( TEMPR . LT. 0.0 ) GO TO 1060
GO TO 60
TRANSFERS WjLL BE MADE FROM THE AQUIFER (SUBSURFACE FLOW > TO THE
(•»«<>» RIVER «« CONDITIONALLY IF THE AQUIFER IN ANY NODE CANNOT MEET THE
r»»»» TRANSFER REOuIRFMENT. AN INTERNODAL TRANSFER TO THE RIVER WILL
Co««o BE ALLOWED ONLY FROM THE NO"E IMMEDIATELY UPSTREAM OF THE NODE
C««»« IN DEFICIT
?90 TEMPR = GRTANK ( KN, 1 ) - CONST
IF ( ABS ( TEMPR ) . LT. 1.0 ) TEMPR = 0.0
IF ( KN. GT. 1 ) GO TO 295
IF ( TEMPR . LT. 0.0 ) GO TO 1070
GO TO l8o
SET QINFLO ARRAY »o LET NODE = KySE ( KN, TS' 24) AND SEQUENCE =
KUSE ( KN, TS , 25 ) . ALSO SET UPSTREAM AQUIFER NODE = KUSE
C«« ( KN, IS, 2" ) FOR INTERNODAL TRANSFER ( IF RFQUIRED )
?95 LN = KUSE ( KN, IS, 28 )
IF ( LN . EO. 0 ) GO TO 1070
KG = NFIND ( LN )
LN = KUSE ( KN, IS, 24 )
IF ( LN . EQ. o > GO TO 10^0
KP = NFIND < LN )
KZ = KUSE ( KN, IS , 25 >
IF ( TEMPR . LT. 0.0 ) GO TO 310
C»««« NO INTERNODAL TRANSFER REQUIRED
DO 300 K = 1, 10
QlNFLO ( KP, KZ, K ) = 0.0
CONTINUE
GO TO 180
AN INTERNODAL TRANSFER IS REQUIRED CHECK CONDITIONS OF AQUIFER IN
UPSTREAM NOPE
310 COMP = ABS ( TEMPR )
COMPA = GRTANK ( KG, 1 ) - COMP
IF ( COMPA . LT. 0.0 ) GO TO 1070
UPSTREAM AQUIFER CAN MFET nEFICIT SET QINFLO AND PRTANK ARRAYS
GRTANK ( KG, 1 ) = COMPA
DO 320 K = ?, 10
QINFLO < KP, KZ, K ) = SRTANK < KG , K >
32n CONTINUE
QINFLO ( KP, KZ, 1 ) = COMP
TEMPR = 0.0
GO TO 180
THIS DEMAND REPRESENTS A SUBSURFACE FLOW THAT IS TRANSFERRED
OUT OF SYSTFM
MO IF ( TEMPR . LT. 0.0 ) GO TO 1020
GO TO 190
THIS IS A RFSERVOIR UPn*TF.
CALL RESOPR ( QINRIV , KN , MO , IW )
MODE = 1
224
-------
PROGRAM ACUMEN
THIS PORTION OF THE PROGRAM INCLUDES THE VOLUME AND CHEMISTRY
C«»»« UPDATES FOR ALL INFLOWS TO THE NODE
40n DO 410 K 5 1, 10
TEMP (K! = "INFLO { KN, IS, K )
4lO CONTINUE
TEST = QINFLO
IF ! TEST . NE. 1.0 ) 60 TO 430
Coo«« THIS IS AN INFLOW TO THE RIVER
CALL CHEMRV ( QINRIV, CONST, KN, IS )
QlNRlVtl) = QINRIV (1) * CONST
GO TO 500
430 IF ( TEST . NF. 2.0 ) GO TO 440
foooo THIS IS AN INFLOW TO THE AQUIFER
CALL CHEMGR < TEMP, CONST, KN, IS )
GRTANK ( KN, D = GRTANK (KN, j ) * CONST
GO TO 500
44o IF ( TEST . NE. 3,0 ) GO TO 450
Coooo THIS IS AN INFLOW TO THE PONDING POOL
CALL CHEMTT ( TEMP, CONST, KN, IS )
TTANK ( KN, 1 ) = TTANK ( KN, I ) * CONST
GO TO 500
450 IF ( TEST .NE. 4.0 ) GO TO 1000
Ceeoo AT THIS POI^T IN THE MODFL CRITERIA FOR THF OPERATION OF THE
Coo»« SYSTEM IS RFQuIREO. THE OPERATIONAL CRITERIA SHOiiLD BE FORMULATED
AND INSERTED AS A SUBROUTINE OR A SERIES OF SUBROUTINES AND/OR
FUNCTIONS WHICH ARE MONITORED BY THE SUBROUTINE CALLED AT THIS
POINT. THE TECHNIQUE AS OUTLINED DOES NOT DESTROY THE
C«x»«o GENERALIZED CONCEPT OF THE MODEL.
C«o«o THIS AN OBSFRVEP CHECK POINT
C»o«e CALL vERNAL SUBROUTINE TO SIMULATE TRANSFERS nF FLOW
BETWEEN SURFACE AND SUBSURFACE FACILITIES TO OBTAIN HYDRAULIC
BALANCE IN MODE
CALL vERNAL ( KN, IS, MO, OBsERV, QINRIV )
CONTINUE
CHECK TO DET£RMlNE IF A SOiL COLUMN EXISTS IN THIS NODE
IF ( TTANK ( KN, 1 ) . LE. 0.0 ) GO TO 540
roo«« CHECK TO DETERMINE IF SOU- COLUMN DATA EXISTS FOR THTS NODE
NSOIL = NSEG (KN )
IF ( NSolL . FQ. 0 ) GO TO 1030
roo*>o SOIL cOLUMN DATA EX'STS PfRCOLATE WATERS IN TTANK THROUGH
SOIL COLUMN INTO BTANK
CALL SOILCO ( KN, Iw, IR, NSOIL )
UPDATE AQUIFER WITH LINEAR MIX WITH WATER IN BTANK
gTANK (KN, 1) = TTANK ( KN, 1 )
DO 5ZO K = 1 , 10
TEM» (K! = BTANK ( KN, K )
520 CONTINUE
CALL CONVER ( TFMP , 2 )
225
-------
PROGRAM ACUMEN
540
•550
•560
157(1
SET UP TRANSFERS Op OUTFLOW FpOM THIS NODE TO NEXT LOWE"
NODE - WHERF ND= NODE NUHBFR OF NEXT LOWER NODE AND N$ =
ENTERING SEQUENCE OF NExT LOWER NODE
ND = NCTRL ( KM, 2 )
IF < KN . EQ. NUMNQD )
ND = NFIND (ND
NS = NCTRL ( KN, 3 )
EO 550 K = 1, 10
QINFLO IND, NS, K)
CONTINUE
CO 570 K = 1 , 10
PREDIC ( KN, K ) = QINRIV (K)
OfiSVAL ( KN, K ) = OBSEPV (K)
CONTINUE
CONTINUE
CALL WRIT ( MO, IYR, Iw )
MO = MO «•
IF ( MO «
GO TO 560
= OBSERV (K)
1? ) GO TO 750
750
770
. LT. TYREND
GT. MOEND )
1 GO T.O 770
CALL EXIT
IW, IR )
inoo
1005
1010
1015
in2o
man
ln35
1040
1045
1055
1H60
1065
1070
( MO,
1
LF.
KO = 1
IYR = IYR * 1
IF ( IYR
IF ( MO
jREAD = 5
CALL PERUSE ( JREAD ,
GO TO 10
PROGRAM CONTROLLED ERROR MESSAGES ARE GENEPATrD HFRE
WRITE t Iw, 1005 ) TEST
FORMAT ( 1H1, 7X, 20HPFSTINATIQN INDEX = ,F3.0,12H IS IN ERRoR )
GO TO 1500
WRITE ( IW, 1015)
FORMAT(1H1,7X,45HDEMAND LEAVING SYSTEM FROM SURFACE IN DEFICIT)
GO TO 1500
WRITE ( IW, 1025 )
FORMAT(lHi,7x,48HDEMAND LEAVING SYSTEM FROM SUBSUPFACF IN DEFICIT)
GO TO 1500
WRITE ( TW , 1035 )
FORMAT ( 1H1, 7X, 25HNO SOIL COLUMN DATA FOUND )
GO TO 1500
WRITE ( IW , 1045 )
FORMAT
-------
PROGRAM ACUMEN
1500 WRITE ( IW, 1510 ) KN » Is
HlO FORMAT < 8X1 1EHFOR NODE = » I3« SXi 10HSFQ.NO. = , 13 )
CALL EXIT
END
227
-------
SUBROUTINE PERUSE UREAD , IW» IR )
COMMON / GENRL / LTlTLE<*>« NCTRH20, 3) , MODES' IYRREG,
AYREND, NUMNHD, !SEQ(20, 20! » KPRI, NUMSEQC?0), NDwTR(?0» 20)
COMMON / RFSER / RDATA(?0, 110), RULF(20, 13), CPNEVPI20, 1?), BA
ANKCo<20)» K*URNM<20, *), GINFLO<20, 10, 1?), HEMAMDI70, 10, 16), K
FMNM(20, in, ») , KQINNM120, 10, 8) , KACUNM (20, P)
COMMON / P"WER / CAP5AB, ENERGY(?0, 8), PDA
ATA<20, 12)
COMMON /
1 GRTANK(20,
?
10)
/ ?OIL<20, 10, Zl),
10), TTANKI20. 10),
, NSEG «?0> i cONUSE (20, 10» ^0
, KUSE (20« 10, 30)
EQUIVALENCE ( CONUSE, KUSE )
DIMENSION LTEMP122), CTEMPU5)
(•»««» SET READ ANP WRITE LOGICAL UNITS
t« uRITE INPUT LISTING HEADING AND VALUES OF JREAD
WRITE ( IW, 10 )
(1H1, 24X, 33H «<>« LISTING CF INPUT DATA
TW , 30 ) JREAD
<1H , RMJREAD = , 13 )
10 FORMAT
30 FORMAT
f.
/ / )
WHEN JREAD s 1 READ AND STORE SYSTEM STRUCTURE - CONCLUDE THIS
C«nn>*»READ WITH CARP CONTAINING THE WORD CONEND IN THE FIRST SIX COLUMNS
C*«»«*WHEN JREAD = 2 PEAD AND STORE TITLES AND INITIAL CONDITIONS OF
C««««»SURFACE RESPRVOIRS AND AQUIFERS ( SUBSURFACE RFSEPVOrRS).
C««*«»ALSO WHEN JREAP = 2 RE«H AND STORE INITIAL SOIL CONTITIONS AS
C«O«»«WELL AS SOIL COLUMN PHYSICAL FEATURES. CONCLUDE THIS READ WITH
C«««««CARD CONTAINING THE WORD COLEND IN TH£ FIRST SlX COLUMNS
C
C««»»»wHEN JREAD = 3 READ AND STORE MATHEMATICAL REPRESFNTATIONS OF
f«««»»PERTINENT SI'RFACE FACILITIES PHYSICAL FEATURES . CONCLUDE THIS
C«n>«««READ WITH CARD CONTAINING EQUEND IN THE FIRST SIX COLUMNS
r
C«o»0«WHEN JREAD = 4 READ INFLOWS AND DEMANDS WITH TITLES FOR THE FIRST Tl
C««««oFRAME ONLY NOTE »# TITLF MUST PRECEDE EACH INFLOW AND DEMAND DATA
r«««««CONCLUDE THIS READ WITH CARD CONTAINING THE WORD FNDATA IN THE
r*««#«FiRST six COLUMNS
c
C««#«»wHEN JREAD = 5 READ INFLOW AND DEMAND DATA CARDS 1USEQUF.NT TO THE
fi«««««F!RST TIME FRAMF AS WELL AS ANY TRANSACTION CARDS" CONCLUDE THTS
C*««READ WITH CARD CONTAINING THE WORD FNDATA IN THE FIRST SIX COLUMNS
60 TO ( 40, 220, 450, 600, 650 ) , JREAD
READ SYSTEM CONTROL CARDS
40 READ ( IR, ftO > LTlTLE, KPRl
" " <«A*, 116)
IW, 80 ) LTITlE i KpRI
<1H , 648, 116)
60 FORMAT
WRITE I
SO FORMAT
JK= 0
JJ s 0
100 READ ( IR, 120 )
IONE . ITWO , LTEMP
228
-------
SUBROUTINE PEPUSE (jREAD , IW, IR )
WRITE ( IW, UO ) IONE , ITWO , LTEMP
FORMAT QH , 2R3,22I3)
STORE SYSTEM CONTROL DATA
IF ( IONE . NE. 3RCON ) GO TO 1000
IF ( ITWO . NE. 3RNOD ) GO TO HO
JK = JK * 1
NUMNOD = JK
DO 150 K = !, 3
NCTRL (JK, r ) = LTEMP (K)
ISO CONTINUE
GO TO 100
160 IF ( ITWO . NE. 3RSEQ ) GO TO 190
JK = NFIND ( LTEMP (1! )
STORE NUMBER OF SEQUENCE POINTS IN A NODE
NUMSEQ UK) = LTEMP (2)
SS = LTEMP (2)
CHECK TO DETERMINE IF NUMBER OF SEQUENCE POINTS EXCEED 20
IF ( NS . GT. 20 ) GO TO 1020
DO l?o * = 1, NS
KX = K + 2
ISEQ (JK, K ) = LTEMP (KX)
170 CONTINUE
GO TO 100
190 IF ( ITWO . NE. 3RWAT ) GO TO 210
JK = NFIND ( LTEMP (1) )
>o°o STORE SYSTE^ FLOW ARRAY
DO 200 K = 1, 20
KX = K + 2
NDWTR ( JK, K ) = LTEMP (KX)
?00 CONTINUE
GO TO 100
710 IF ( ITWQ . NE. 3REND ) GO TO 1040
STORE DATES FOR BEGINNING AND END OF STUDY P£RlOD
MOeEG = LTEMP (1)
IYRBEG = LTEMP (2)
MOEND = LTEMp (3)
IYREND = LTFMP (4)
GO TO 2000
READ AND ST^RE INITIAL VALUES AND MATHEMATICAL REPRESENTATIONS
OF SYSTEM PHyslCAL FEATURES AS W?LL AS TITLES OF *AME IF
f.oooo APPLICABLE - SOIL COLUMN INITIAL ANT PHYSICAL DAT* ARE iNCLuDEr
IN THIS PORTION OF THE SUBROUTINE
?20 READ (IR, ?30> IONE, ITWO, (LTEMP(I), I = 1, 11)
??0 FORMAT (2R3.I2, IX, 13, R3, 8A8)
WRITE (IW, 240) IONE, iTWOi (LTEMP(I), I = 1, 11)
?40 FORMAT (1H , 2R3,I2, IX, 13, R3, 8*8)
IF ( IONE . NE. 3RAQU ! GO TO 260
STORE TITLE OF AQUIFER ( SUBSURFACE RESERVOIR )
JK = NFIND < LTFMP (2) )
229
-------
SUBROUTINE PERUSE (JREAD , IW, IP )
KX = K + 3
KAQUNM ( JK, K ) = LTEMP ( KX )
750 CONTINUE
GO TO 280
?60 IF ( IONE . NE. 3RSUR ) GO TO 370
STORE TITLE OF SURFACE "FSERVOIR
JK = NFIND ( LTFMP (2) )
DO 270 * = 1, 8
KSURNM ( JK, K ) = LTEM" (KX)
?70 CONTINUE
??.0 READ (IR, ?90) IONE, ITWO, (LTEMP(I), I = 1, 4>i (CTEMP(J), J = 1
A, 10) , (LTEMP (K) , K = 5, 7)
?90 FORMAT (2*3,12, II, 13, R3» 10F6.0, 12, II, 12)
WRITE (IW, 300) IONE, ITWO, (LTEMP(I), I = 1, 4), (CTE*PU), J =
Al, 10), (LTFMP(K), K =5, 7>
FORMAT Uh , 2P3,I2, II, 13, R3, 10F6.0, 12, II, 12)
IF ( IONE . NE. 3RRES ) Go TO 340
STORE DATA FOR SURFACE RESERVOIR IN RDATA ARRAY
JK r NFIND ( LTFMP <3>1
LET NExp REPRESENT POwER OF 10 USED IN REPORTING FLOw MAGNITUDES
NF.Xp = LTEMP (6)
FAC = 10.0 «« NFXP
RDATA (JK, 1 ) = LTEMP (3)
RDATA (JK, ? ) = LTEMP (1)
RDATA ( JK, 3) = CTFMP (1) « FAC
DO 310 K = ?, 5
KA = K + 9
RDATA ( JK, KA ) = CTEMP (K) * FAC
MO CONTINUE
RDATA ( JK, 20 ) = CTEMP (6) « FAC
CARD CONTAINING RESERVOIR CHEMISTRY DATA MUST ALWAYS FOLLOW
PHYSICAL DATA CARD - UNITS OF CONCENTRATION CAN BF IN MEQ/ LITER
OR PARTS PE" MILLION LTEMP (2) CONTAINS THIS UNIT? INDEX
READ (IR, ?90) IONE, ITWO, (LTEMP(I), I = 1, 4), (CTEMP(J), J = 1
A, 10), (LTEMP(K), K = 5, 7)
IF ( ITWO .NE. 3RCHM ) GO TO 1060
WRITE (IW, 300) IONE, ITWO, (LTEMP(I), I = 1, 4), (CTFMP(J), J =
Al, 10) , (LTFMP(K) , K = 5, 7)
JK = NFIND (LTEMP (3 ))
CHECK TO DETERMINE IF ONLY TOTAL SALTS IN PPM ARE REPORTED
IF ( CTEMP (10) . GT. 0.0 ) GO TO 320
CONCENTRATIONS OF CONSTITUENTS ARE RFPORTFD CONVERT TO
AND COMPUTE TOTAL SALTS ON BASIS OF CONCENTRATIONS
CALL CONVER ( CTEMP , LTEMP (2»
320 DO 330 K = 2, 10
KX = K + 99
RDATA (JK, KX ) = CTEMP (K)
330 CONTINUE
GO TO 220
230
-------
SUBROUTINE PERUSE (JREAD , IW, IR )
350
360
C«C«
380
390
4on
450
4.6 0
STORE INITIAL CAPACITY ANp CHEMISTRY FOR AQUIFER
JK = NFIND ( LT^MPOn
CHECK TO DETERMINE IF ONLY TOTAL SALTS IN PPM ARE REPORTED
IFICTEMP (10) . GT. 0.0 ) 60 TO 350
CONCENTRATIONS OF CONSTITUENTS ARE REPORTED CONVERT TO PPM
AND COMPUTE TOTAL SALTS ON BASIS OF CONCENTRATIONS
CALL CONVER ( CTEMP, LTEMP (2) !
DO 360 K = 2, 10
GRTANK ( JK, r ) = CTfM» (K)
CONTINUE
NEXP = LTEMP (ft)
FAC = 10-0 *« NFXP
GRTANK ( JK, 1 ) = CTEM^ (D « FAC
GO TO 220
RFAD SOIL DATA FOR SOIL COLUMN - EACH SEGMENT iN ?OlL COLUMN MUST
BE IN PPOPED SEQUENCE. SFGMENT NUMBER ONE IS SEGMENT NEAREST
GROUND SURFACE.
IF ( IONE . NE. 3RCOL ) GO TO 1080
IF ( ITW0 .FQ. 3REND ) GO TO *000
IF ( ITWO . NF.. 3RREG ) PO TO 1100
LSEG = 0
KNOtlE = LTE^P(2)
JK = NFIND (KNOflE )
NSEG ( JK) = LTFMP (1)
LNODE = KNOPE
MAX = NSEG (JK)
DO 400 J = 1, MAX
READ(IR,380) IONE, ITWO. (LTEMP (KA), KA=1,2), (CTFMP (KB) ,KB=1 , 15)
FORMAT(2R3,2I3,SF5.0,F3.0,3F5.0, Fl.O, F^t.n,F5.0)
WRITEdW, 390) TONE, ITWO, (LTEMP (KA), KA=!,?), (CTFMP (KB) )
FORMAT(1H ,?R3, 213, 8F5.0, F3.0, 3F5.0, Fl.O, FA.O, F5.0 )
IF ( IONE .NE. 3RCOL ) GO TO 1100
60 TO 1100
IF ( ITWO .NE. 3RSEG
KNOHE = LTEMP (1)
IF ( KNODE . NE. LNODE
LSEG = LSEG +
LSEG .NE
) GO TO 1100
LTEMP (2)) GO TO 1100
( JK,
I5
K
) = CTFMP (K)
1
IF
KS s LSEG
DO
SOIL
CONTINUE
GO TO 220
READ AND STORE MATHEMATICAL REPRESENTATIONS OF PHYSICAL FEATURES
OF SURFACE FACILITIES
READ ( IR, 460 ) IONE . ITW0 , KNOD , (LTEMP (K) , CTEMP (K) ,
K= 1, 4 ) , CTEMP <5>
( 2R?, 14, 4 ( 13, F12.0 ) , F10.0 )
IW, 46? > IONE, ITWO, KNOD , ( LT^MP 'K> , CTEMP(K) ,
K = 1,4 ), CTEMP (5)
FORMAT
WRITE l
231
-------
SUBROUTINE PERUSE (JRF.AD , IW. IP )
IF ( TONE .NE. 3REQU ) GO TO 1120
IF ( ITWO . EG. 3REND ) GO TO 2000
JK = NFIND ( KNOD )
DO 560 K = 1, 4
KX = LTEKP
GO TO 560
490 IF ( ITWO . NF. 3REFF ) GO TO 500
roooo STORE EQUATIONS OF HORSFPOwER AND DISCHAPGF ( AT FULL GATE)
ro«oo AS FUNCTION? OF NET HEAD
ENERGY ! JK, KX) = CTEMP (K)
GO TO 560
500 IF ( ITWO . NE. 3RCAP ) GO TO 510
Coo»o STORE CONSTRAINTS CONCFRNING TURBINF CAp^p i[_ I TI ES
CAPBAB (JK, KX) = CTEMP = CTEMP (K)
IF ( KX. EQ. 12 ) RULE < JK, 13 ) = CTEMP («5)
GO TO 560
530 Ip < ITWO . Nf . 3RPVP ) GO TO ?40
Coooo STORE RESERVOIR EVAPORATION COEFFICIENTS IN UNITS OF FEET PER
C««
-------
SUBROUTINE PERUSE UREAD , 1W, IR )
f « « » o
600
GO TO 450
READ AND STORF INFLOW, DEMAND, AND CONSUMPTIVE
TIME FRAME ONLY WHERE TITLE CARDS MUST PPFTEDE
. f,T. 4 ) GO TO 650
• 230 > IQNE , ITWO * < LTEMP ,
. 240 ) IONE , ITWO, ( LTFMP (K) , K
. EO. 3REND ) GO TO 2000
• NF. 3RSUR ) GO TO 1180
IF ( JREAD
READ ( IR
WRITE ( IW
IF ( IONE
( IONE
= MFIND
= LTEMP
( I TWO
USF CARDS FOR
DATA CARDS
K = 1, 11 )
= 1. 11 )
FIRST
IF
JK
IS
IF
)
(1)
NE
STORE TITLE OF
IS, K) = L^EMP «X)
610
e«o
62n
6.30
o««
65n
Of"
, 8
INFLOW SOURCE
IS, K ) = LTFMP
IF ( IONE . EQ. 3REND ) GO TO 2000
JK = NFlND ( LTFMP (3) )
IS = LTEMP (1)
MO = LTEMP (5!
RAISE MAGNITUDE OF DEMAND BY POWER OF 10 ( LTEMP
NEXP = LTFMP (6)
FAC = 10.0 «« NFXP
CTEMP 11) = CTEMP (1) * FAC
IF ( IONE . NE. 3RDEM ) GO TO 700
STORE DEMAND DATA WITH OR WITHOUT REQUIRED CONSUMPTIVE USE DATA
DEMAND ' JK' IS, MO ) = CT£MP <1>
<=TOPE SOURCE INDICATOR TN LOCATION 13 OF DEMAND APRAY
SYMBOLICALLY THERE ARE ONLY TWO SOURCES SUPDLYING DEMANDS
NAMELY , SURFACE AND SlJBSl)RFACE SOURCFs
IF ( ITWO . EQ. 3RSUR ) DEMAND ( JK, Is, 13 ) = 1.0
IF ! ITWQ . EQ. 3RGWR ) DEMAND < JK, IS, 13 ) = 2.0
IN LOCATION 14 OF T^E DFMAND ARRAY IS STORED THE USE IE
OR MUNICIPAL AND INDUSTRIAL OR THE DESTINATION IN A
FLOWS FROM FACILITY TO FACILITY
IF ( LTEMP (4) . EQ. 3RIRP > DEMAND I JK, IS, 14 )
IF ( LTEMP <4> • EQ. 3RMAI ) DEMAND ( JKi IS, 14 )
(CTEMf(j), j = 1
(CTEMPSJ), J =
(6) !
IRRIGATION
TRANSFER OF
= 1.0
= 2-0
233
-------
SUBROUTINE PERUSF UREAD , iw, IP s
IF ( LTEMP U) . EQ. 3RGWV )
IF ( LTEMP (4) , EQ. 3RSUR >
IF ( LTEMP (4) . EQ. SRfiUT )
LTEMP (2) CONTAINS INDICATOR
DEMAND t JK, is, 14 > = 4.0
DEMAND ( JK, IS, 14 ) = 5•0
DEMAND ( JK, IS* 14 ) = 6.0
DESIGNATING CONSUMPTIVE USE pATA
<-oo»o CR TRANSFER DATA WILL
IF ( LTEMP (2) .NE. 1
READ ( IR, f60 ) ICHK
OR WILL NOT FOLLOW PEMAND CARD
) GO TO 600
(LTEMP(K) ,K=1,2) , (CTEMP(KA) ,KA=l,3) ,
670
1 (LTEMP (KB) , KB = 3, 16 )
FORMAT ( Re, 13, 12, F9.0, 2F3.0, 3
1 IX, 12 )
WRITE ( IW, 670) ICHK,(LTFMP(K),K=1
1 (LTEMP (KB) , KB = 3, 16
FORMAT(1H ,P6,I3,I2, F9.0, 2F3.0, 3
1 IX, 12 )
CONSUMPTIVE USE OR TRANSFER DATA WILL
( R3, 13. I?, 13) , IftX, 12
2) . (CTEMPdcA) ,KA = 1,3) ,
)
( R3, 13, 12, 13) , 16X, 12
680
700
f « o oo
(*#<}<>«
f « a
710
BE READ AND
IF ( ICHK .
JK = NFIND
IS = LTEMP
MO = LTEMP
CONUSE (
CONUSE (
CONUSE (
EO 680 K
KA = K *
KUSE ( Jlf
JK
JK
JK
=
16
i ,
NE. 6RCONUSE ) GO TC
( LTEMP (1) )
(2)
(
,
,
,
3,
TS
15)
IS
IS
IS
*
1
1
MO
13
14
) =
) =
) =
CTFMP
CTEMP
CTEMp
(1)
(2)
(3)
14
*
if A
)
= LTEMP (K)
1160
CONTINUE
GO TO 600
IF ( IONE . NE. 3RQIN ) GO TO 740
STORE INFLOW DATA WITH OR WITHOUT BREAKDOWN Oc CHEMICAL
CONSTITUENTS - LTEMP (2) DESIGNATES UNITS OF INPUT I.E.
LlTFR. AND 2 = PPM
LOCATION 11
THE PARTICU
IF ( LTEMP
( LTEMP
( LTEMP
( LTEMP
710 K =
QINFLC
= 1
= 2
= 3
= it
.0
.0
.0
.0
1 = MEQ/
ON OF
(i)
IS'
LE.
1 )
0.0 )
= CTEMP
GO TO
(1)
600
CONSTTTUENTS
234
-------
SUBROUTINE PERUSE IJREAD
IR
730
f !> O <> 9
750
760
770
780
790
KOO
BIO
f O O it »
inoo
imo
1120
1P30
EO 730 K = ?, 10
QINFLO (JK, IS, K ) = CTEMp (K)
CONTINUE
QINFLO ( JK, IS, 12 ) = LTEMP (2)
GO TC 600
IF { IONE . NE. 3RTRN ) GO TO 1200
STORE DATA !N TRANSACTION CARDS oo* ONLY ONE VARIABLE OR "ApAMETFR
CHANGE CAN BE ACCOMMODATED PER CARD
KA = CTEMf (1)
IF { ITWO . NE. 3RRFS ) GO TO 750
RDATA ( JK, KA ) = CTEMP (2)
GO TO 600
IF ( ITWO . NE. 3RTDA I GO TO 760
TDATA ( JK, KA) = CTEMP (?)
GO TC 600
IF ( ITWO . NE. 3RPDA ) GO TO 770
PDATA ( JK, KA) = CTEMP (2)
GO TO 600
1050
jn6n
1070
IF ! ITWO .
ENERGY ( JK,
GO TO 600
IF ( ITWO .
CAPPAB < JK,
GO TO 600
IF ( ITWO .
BANKCO (JK)
GO TO 600
IF ( ITWO .
CONEVP ( JK,
GO TO 600
IF ( ITWO .
RULE t JK, *
GP TO 600
ALL PROGRAM
WRITE ( IW,
FORMAT ( iHj
CALL EXIT
WRITE ( IW ,
FORMAT ( 1H1
CALL EXIT
«RITE ! IW ,
FORMAT < IHI
CALL EXIT
WRITE ( IW ,
FORMAT ( 1H1
CALL FXIT
WRITE I IW ,
FORMAT ( jHi
CALL EXIT
NE. 3RENE ) GO TO 780
KA ) = CTEMP (2)
NE. 3RCAP ) GO TO 790
KA) = CTFMP (2)
NE. 3RBNK ) GO TO
= CTEMP (2)
NE. 3REVP ) GO TO 810
KA) = CTEMP (2)
NE. 3RRUL ) GO TO 1?00
A) = CTEMP (2)
CONTROLLED ERROR MESSAGES ARE GFNERATFD HERE
1010 )
, 7X, 26HTHIS IS NOT A CONTROL CARn 1
1030 )
, 7X, 36HNUMBER OF SEQUFNCE POINTS F.XCEEDS 20 )
1050 !
, 7X, 36HNO END SIGNAL FQUND FoR CONTROL 1ATA )
, 8X,
10
-------
SUBROUTINE PERUSE (jREAf) , IW, IR )
1110 FORMAT ( 1H1, 7*, ?qHSO!L COLUMN DATA OUT OF ORDEP )
CALL FXIT
11?0 WRITE ( IW, 1130 )
mo FORMAT t IHI, 7x, PPHEQUATION CARDS OUT OF ORDER )
CALL FXIT
ll'iO WRITE ( IW, 1150 !
lj-0 FORMAT ( jHi, 7x, 34HITWO INDEX NOT FOUND IN VALIP LIST )
CALL EXIT
1160 WRjTE ( IW , 1170)
1170 FORMAT ! 1H', 6X, 30HNO CONSUMPTIVE USE CARD FOUNp )
CALL EXIT
11PC WRITE ( IW , 1190 )
1190 FORMAT ( 1H1, 7X, 39HINFLOW OR DFMANE CARDS ARE OUT OF ORDER
CALL FXIT
)?00 WRITE t IW% 1?10 )
l?!n FORMAT ( 1H1, 7X, 34HTRANSACTI ON CARDS ARE OUT OF ORDER )
CALL EXIT
1??0 WRITE (IW , 1?30 !
1?30 FORMAT ( 1H1, 7", 34HDATA F°R JREAD = 2 IS OUT OF ORDER >
CALL EXIT
?nr,o RETURN
END
236
-------
SUBROUTINE WRIT ( MO, IYR , IW I
COMMON / GFNRL / LTITLE(P), NCTPH20, 3), *0*EG, IYRREG, MOENP, I
6YREND, NUMNOD, TSEQ(20. 20), KPRI, NUMSECI20), NDWTR(20, 20)
COMMON / RFSF.R / RDATA(?0, 110), RULEI20, 13), cnNEVP(20, 1?), BA
ANKC0120), KSURNMI20, 8), QINFLOI20, 10, 12), lEMANrHZO, 10, 16), K
PEEMNMI20, IP, 6), KQINNM(?0, 10, 8) , KAQUNM (?0,P)
COMMON / POWER / CAPEAB(20, 20), TDATAI20. 4> , ENERGYI20, 8), PDA
ATA(20, 12)
COMMON / PROPER / soiL(2o, 10, 211, BT*NK(?O, io>, TTAMMPO, 101,
1 GRTANM20, 10) , NSrG (?0) , CONUSE (20, 10, 30)
? , KUSE (20, 10, 30)
COMMON / PAPAM / BEGRES ( 20, 11 ) , BEGWp ( ?0, 10 ) ,
1 PREPlC ( 20, 10 ) , OBSVAL (20.10), CHFMDM ( ?0, 10, 10)
CIMFNSION MON (12) , LTEMP (ZO ) , TEMP (10)
EQUIVALENCE ( CONUSF, KU?E )
CATA (( MON (K) , < = 1, 12) = 3RJAN, 3RFEB, 3RMAP, 3RAPR, 3RMAy.
1 3RJCN, ^RJUL, 3RAUG, 3RSEp, 3PQCT, 3RNQV, 3RDEC 1
DATA ( iPAGr = 0 )
PONTH = MON (MO )
KYR = IYR + 1900
to 700 KN = i, NUMNOD
IPAOE = IPA^-E + 1
WRITE ( IW, 50) UTITLE (K) , K = 1, 8 ) , IPAGE
50 FORMAT! 1H1, 7X, 8A8, 3X, 19HNUMEEP CF NOHE5 =3 . 6X,
i PHPARE NO. , 13 / j
KNODE = NCTPL (KN,D
WRITE ( IW, 60 ) KNODE. MONTH, KYR
60 FORMAT ( 8X, 14HNODE NllMBFR = , 13, 5X, 9HMONTH OF . R3, 10X,
1 5HYEAP , 14 // ?8X, 6HVOLUME, 4X, 2HCA , ^X, 2HMG, 4X,
? 2HNA, AX, 2HCL, 4X, 3H504 , 2X, 4HHC03 , 3X, 3HC03 , 3X,
3 3HN03 , 3X, 11HTOTAL SALTS / 56x, 9HACRE FEFT,
4 9 !3X, 3HPPM ) , 3X, 7HTONS/AF / )
WRITE I IW, 70)
70 FORMAT ( 8X, 42HOPERATIONAL SEQUENCE OF SURFACE FACILITIES / )
95 GO TO ! 100, 11", 130 ) , KN
JOO fA = 4
MB = 5
PC = 9
MD = 10
GO TO 150
110 MA = 2
MB = 3
MC = 7
MB = 8
GO TO 150
130 MA= 2
MB = 3
MC = 5
KD = 6
150 NS = NUMSECl ( KN !
237
-------
SUBROUTINE WRIT ( MO, IYP. . IW )
DO 490 KS = 1, NS
IS = TABS ( ISEO ( KN, KS ) )
WRITE OUTPUT 0F LINES CONTROLLED BY INDEX(MA) (INFLOWS- DEMANDS)
IF ( KS. GT. HA ) GO TO 190
160 IF ( ISEQ ( KN, KS ).GT. 0 ) GO TO 450
IF ( ISEQ ( KN, KS >tLT. 0 ) 60 TO 300
GO TO 490
190 IF ( KS . NT. MF» ) GO TO 220
WRITE ( IW , 210 )
?10 FORMAT (/8X, 28HOBSERVEP OUTFLOWS FPCM NODE / )
??0 IF ( KS. LE. MC ) GO TO 450
IF ( KS • NE. MO 1 GO TO 250
WRITE ( IW , 230 )
?30 FORMAT (/8X, 4QHSUBSURFACE OPERATIONS AND FLOW TRANSFERS / )
LTEMP ( 1 ) = 6HAQUIFER
LTEMP t 2 ) = 8H CONDITT
LTEMP ( 3 ) = *HONS OF t.
LTEMP ( 4 ) s 8HAST TIME
LTEMP ( 5 ) s 8H FRAME
DO 240 K = 1, 10
KA = K «• 5
LTEMP (KA) = REGWR ( KN, K )
?40 CONTINUE
TONS = BEGWP ( KN, 10 ) » 1.36E-3
WRITE) IW, 4io ) (LTFMP (K) , K = 1, 15 ) , TONS
?50 IF ( KS . LF. NS ) GO TO 160
THIS PORTION OF THE SUBROUTINE SETS UP ALL DEMANDS FROM SYSTEM
N CORRECT OUTPUT MODE
300 DO 320 K = 1, 5
LTEMP (K) = KDEMNM ( KN, IS» K )
^20 CONTINUE
LTEMP (6) = DEMAND ( KN , IS, 16 )
330 DO 34Q K s ?, 10
KA = K * 5
LTEMP (KA) = CHEMDM ( KN, IS, K )
340 CONTINUE
TONS = CHEH"M (KN, IS, 10 ) » 1.36 E-3
THIS IS THE GENERAL OUTPUT STATEMENT FOR MOST OF THE OUTPUT
400 IF ( MARK .RT. 8 ) GO TO 408
IF ( LTEnp (6) .GT. 0 ) GO TO 408
TONS= 0.0
DO 405 * s f, 15
LTEMP (K) s 0
405 CONTINUE
408 WRITE ( IW, 410 > ( LTEMP (K) * K = 1, 15 ) , TON?
410 FORMAT ( 11X, 5ft8 » 5X, 19, 8 (2X, 14) , 2X, 15, ?X, F6.3 )
IF ( ISEQ ( KN, KS ) . GT. 0 ) GO TO 490
IF ( DEMAND ( KN, IS, i4 ) • GT. 2.Q ) GO TO 490
r»«»» THIS PORTIO« OF THE SUBROUTINE COMpUTFS AND OUTPUTS THE SHORTAGES
238
-------
SUBROUTINE WRIT ( MO, IYR
IW
42n
14X, 19 )
460
470
KN, IS , K )
KTEM = DEMAND ( KN , IS, MO )
LTEMP (17) = KTEM - LTEMP (6)
WRITE ( IW, 420 ) LTEMP (17 )
FORMAT ( 11*, 31HSHORTAGE FROM THE IDEAL DEMAND
GO TO 49n
THIS PORTION OF THE ROUTINE SETS UP ALL INFLOWS Tn THE
(•<>«<>« SYSTEM IN THE CORRECT OUTPUT MODE
45(1 EO 460 K = 1, 5
LTEMP (K) r KGINNM (
CONTINUE
LTEMP (6) = QINFLO ( KN, IS , 1 )
DO 470 K = ?, 10
KA = K •»• 5
LTEMP (KA) = QINFLO ( KN, IS, K )
CONTINUE
TONS = QINFLO < KN, is, 101 « 1.36E-3
GO TO 400
CONTINUE
BEGIN THE COMPARISON OF OBSERVED ANP PREDICTED ENP OF NODF
FLOWS AND CHEMISTRY
500
•520
530
540
WRITE
FORMAT
LTEMP
LTEMP
LTEMP
LTEMP
LTEMP
CO 530
KA = K
( IW , 520 )
( // BX, 16HCOMPARISON INnEX
( 1) = SHTQTAL OB
( 2) = 8HSERVED 0
( 3) = 8HUTFLOWS
( 4) = 8HFROM NOD
( 5) = SHE
K = 1, 10
5
) = O^SVAL ( KN , K )
( KN
LTEMP (KA
CONTINUE
TONS = OBSV^L
WRITL
LTEMP
LTEMP
LTEMP < 3)
LTfMP ( 4)
LTEMP ( 5)
EO 540 K =
KA = K * 5
LTEMP (KA ) = PREDK ( KN, K )
CONTINUE
10 ) »
: iw . 4in ) (LTEMP
(1) = 8HPREDICTF
( 2) = PHP OUTFLO
= BHW FROM T
= 8HHIS NODE
= RH
1, 10
1.36E-3
(K) , K =
1,1?) , TONS
TONS = PRED'C ( KN
WRITE
LTEMP
LTEMP
LTEMP
IW
1)
2)
3)
4)
5)
10 ! <
410 ) (LTEMP
RHSIMPLE D
8HTFFERENC
8HF (OPSER
8HVED - PR
RHFDICTED)
1.36E-3
(K) , K = 1
15 )
TON'S
239
-------
SUBROUTINE WRIT ( MO, IYR , IW )
TEMP - PREDIC i KN, K >
KA = r+ 5
LTEMP (KA ) = TFMP (K)
CONTINUE
TONS = TEMP (10) » 1.36F-3
WRITE ( IW , 4io ) (LTEMP (1C) , K = i, i5 ) , TONS
WRITE ( IW, 5f>0)
560 FORMAT ( / PX, 25HCHEMiCAL CHANGES IN NODE / )
LTEMP (1) = 8HOBSERVED
LTEMP (2) = 8H CHANGE
DO 570 K = ^, 5
LTEMP (K) = 6H
570 CONTINUE
DO 580 K = i, in
TEMP (K! = PBSVAL ( KN, K ) - QINFLO (KN, 1, K )
KA = K + 5
LTEMP (KA) - TE^P (K)
5BO CONTINUE
TONS = TEMp (10) « 1.36E-3
LTEMP (6) = 0
WRITE ( IW , 410 ) (LT^MP (K) , K = 1 , 1 e, ) , TONS
LTEMP (1) = BHPOEPICTE
LTEMp (2) = 8HD CHANGE
DO 590 K = ] , m
TEMP (K) = PRFDIC (KN, K ) - QINFLO (KN, 1, K )
KA = K + 5
LTEMP (KA) = TEMP (K)
5
-------
SUBROUTINE WRIT ( MO, IYR . IW )
CONTINUE
LTEMP (6) = 0
TONS = TEMP (10 ) « 1.3fcE-3
WRITE ( IW , 410 ) ILTEMP (K) , K = 1« 15 ) « TONS
700 CONTINUE
?(inn RETURN
END
241
-------
SUBROUTINE RESOPR ( QFLOW , KN, MO, IW )
COMMON / GENRL / LTlTLF(P), NCTPLI20, 3), MOBEG, IYRBEG, MOEND, I
AYRENE. NUMNOD, TSEQ120, 20)i KPRI» NUHSEO(?0), NDWTR120, 20)
COMMON / RFSER / RDATA
DATA ( (IDAYS(K),K = 1,12) = 31,28,31'^0,3l,30,31,31.30,11,30.31)
CATA ( IMON = 0) , ( CONST = 1.963471) , ( IFlR = 0 ) ,
1 ( EPSLON = 1.0E1 )
THIS SUBROUTINE SIMULATES THE OPERATION OF A "ESERVOIR WITH OR
WITHOUT A POWER PLANT - ASSUMES CHEMISTRY MIX ( LINEAR) FOR EVERY
INFLOW IE WHEN QFLOW (1) IS GREATER THAN zERO . HOWEVER, THE
UPDATE is SIMULATED ONCF PER MONTHLY TIME FRAME. THE SEDIMENT
INpLOW AND INCREMENTAL DEPOSIT IS COMPUTED DURING THE INITIAL
ITERATION ONLY- THIS ASSUMES THAT SUBSEQUENT INFLOWS ARE FROM
UPSTREAM RESERVOIR RELEASES AND CONSIDERED SEDIMENT FREE . THE
MINIMUM RESFRvOIR CAPACITY IS INCREASED AT THE BEGINNING OF FACH
TIME FRAME ON THE BASIS OF THE SEDIMENT PEPOSlTED IN THE PREVIOUS
••>« TIME FRAME. IN THE ARGUMENT LIST - QFLOW IS AN ARRAY CONTAINING
*o VOLUME AND CHEMISTRY OF RESERVOIR INFLOW OR OUTFLOW, KN = NODE
»« INDEX- M0 = MONTH , Iw = LOGICAL WRITE UNIT
KNODE = NCT<>L ( KNi 1)
IF ( IFlR . NF. 0 ) GO TO 20
CO 10 J = 1, 20
CO 10 K = 1, 11
BEGRES ( J, K) = 0.0
10 CONTINUE
IFIR = 1
CHECK TO DETERMINE IF THIS IS THE FIRST CALL FOR ANY RESERvOTR
IN THIS TIMF FRAMF
IF ( IMON . EQ. MO ) GO TO 35
INITIALIZE <"ALL COUNTERS FOR ALL RESERVOIRS
DO 30 K = 1, 20
MODE (K) = 0
so CONTINUE
DAYS = IDAY? (MO)
IMON = MO
CHECK TO DETERMINE IF THIS IS THE INTITIAL CALL FOR THIS
REsERyOlR IF SO UPDATE RESERVOIR CHEMISTRY AND
35 IF ( MODE (KN) . NE- 0 ) GO TO 100
STQRE RESERVOIR CHEMISTRY oF PREVIOUS TIME FRAME
EO 40 K s 2, 10
KA s K * 99
BEGRES
-------
SUBROUTINE PESOPR ! QFLOW , KN, HO, IW )
BEGRES ( KN, ] ) = RDATA (KN, 3 )
UPDATE MINIMUM PESERVHIP CAPACITY ON BASIS OF SEDIMENT DEPOSITED
C*o DURING PREVIOUS TIME FRAMF
RDATA ( KN, 1?) = RDATA ( KN, 12 ) * BEGRFS ( KN, 11 )
UPDATE TOTAL SEDIMENT DEPOSITED
RDATA ( KN, 20) = RDATA ( KN, 20 ) * BEGRFS < KN, 11 )
DURING INITIAL ITERATION QFLOW (1) MUST EQUAL TO PR GREATCP THAN
ZERO ( MUST NOT REPRESENT A RELEASE )
IF ( QFLOW (1) . LT. 0.0 ) GO TO 1000
COMPUTE SEDIMENT INFLOW AS A FUNCTION OF WATFR INFLOW WHE"F
QSED (TONS/nAY) = A « QCF5 «« B
QSED = o.o
A = RDATA (KN, 15)
CHECK TO DETERMINE IF SFDIMENT is TO RE CONSTDERET FOP THIS
RESERVOIR
IF ( A . EQ. 0.0 ) GO TO 60
IF ( QFLOW (1! . EQ. 0.0 ) &0 TC 60
B = RDATA ( KN, 16)
QCFS = QFLOW (1) / CONST/ DAYS
TEST TO DETERMINE IF THERF ARE TWO FUNCTIONS FOP SEDIMENT AND
WHICH OF THF TWO IS REQUIRED FOR THF SEDIMENT COMPUTATION
IF ( QCFS . LT. RDATA (KN, 19 ) ! GO TO 50
A = RDATA (KN, 17)
B = RDATA (KN, IB )
QSED = A * OOFS *« B
CONVERT FROM TONS PER D&Y TO ACRE FEET PFR MONTH
QSED = QSFD « DAYS « 6.1218E-4
BEGRES < KN, 11 ) = QSED
CHECK TO DETERMINE IF RFSERVOIR AREA AND WATER SURFACE ELEVATION
ARE INITIALIZED
IF ( RDATA ( KN, 4 ) ,
CAp = RDATA (KN, 3 )
IF ( CAP . LE' 0.0 )
RDATA (KN, 4) = RFSEL
RDATA (KN, M = AREA
INITIALIZE PESFRVOlR EVAPORATION, CHANGE IN BANK STORAGE, SPlLL,
RELEASE THROUGH OUTLETS AND TUR?INES AND TOTAL RELEASE
BO 90 K = 6, 10
RDATA (KN, * ) = 0.0
CONTINUE
STORE RESERVOIR CAPACITY, WATER SURF/sCE FLEVflTION AND AREA AT
END OF LAST TIME FRAME AS INITIAL VALUES FOR THIS TIME FRAME
PRECAP (KN) = RHATA (KN, 3)
pREL (KN! = RDATA (KN, A)
PARFA (KN) = RDATA (KN. 51
SET TARGET CAPACITY FOR THIS RESERVOIR FOR THIS TIME FRAME
THEsE TARGET CAPACITIES ARE READ IN AS RULE CURVF*.
TARGET = RULE ! KN, MO )
IF TARGET UFVEL IS IN TERMS 0^ WATER SURFACE FLEViTI0N FIND
r«o<>«
50
(•««««
60
. C-T. 0.0 ) GO TO PO
GO TO 1020
( CAP, AREA, KN, IW )
9n
90
mo
243
-------
SUBROUTINE RESOPR ( QFi_OW , KM, MO* IW
f o
C 00
THIS IS THE FIRST CALL FOR THIS RESERVOIR IN THIS TIME FRAME
AND REPRESENTS CONDITION (1) - SET STORAGE ( STORAF) AND RELEASE
( REALF ) AS INITIAL VALUES IN ACRE FEET
STQRAF = QFLQW (1)
RELAF = RELWIN « CONST « DAYS
COMPUTE A SEED FOR RESERVOIR CAPACITY TO START ITERATIVE SCHEME
TO BALANCE "ESFRvOIR BUDGET - LET NITER = NUM°ER OF ITERATIONS
REQUIRED TO OBTAIN BALANCE
CAPX = TEMCAP «• STORAF - RELAF
NITER = 0
NITER = NITFR f 1
COMPARE SEEP CAPACITY W?TH TARGET CAPACITY
IF 1 CAPX . LF. TARGET 1 C-0 TO 220
SEED CAPACITY IS GREATER THAN TARGET CAPACITY - INCREASE RELEASF
DELTAF = CAPX - TARGET
RELAF = RELAF + DELTAF
SET CAPACITY FQUAL TO TARGET CAPACITY
CAPX = TARGET
244
-------
SUBROUTINE RESOPR ( QFLOW , KN, MO, IW )
RELEASE
RELCFS = RELAF / CONST / DAYS
IF ( RELCFS . LE. RELMAX ) GO TO 240
RELEASE IS r,RFATER THAN MAXIMUM SET RELEASE EQUAL TO MAXIMUM
RELCFS = RELMAX
RELAF = PELTS * CONST » DAYS
MARK = 1
COMPARE CAPACITY WITH MAXIMUM CAPACITY
IF ( CAPX . LE- CAPMAX ) GO TO 260
foo»o CAPACITY IS GREATER THAN MAXIMUM CAPACITY INCREASF RELEASE IF
Cooo« POSSIBLE ANP SET CAPACITY EQUAL TO MAXIMUM CAPACITY
DELTAF = CA^x - CAPMAx
CAPx = CAPM&X
IF ( MARK . ED. 1 ) GO TO 280
RELAF = RFLAF t- DELTAF
GO TO 220
<-ooo« COMPARE CAPACITY WITH MINIMUM CAPACITY
?60 IF ( CAPx -r-E. CAPMIN ! GO TO 300
(-0000 CApACITY IS LESS THAN MINIMUM CAPACITY DECREASE RELEASE
DELTAF = CAPMIN - CAPX
RELAF = RELAF - DELTAF
RELCFS = PELAF / CONST / PAYS
SET CAPACITY TO MINIMUM CAPACITY
CAPX = CAPMIN
COMPARE REDUCED REl EASE WITH MINIMUM RELEASF
IF ( RELCFS .GE. RELMIN ) GO TO 300
IF MORE THAN TWO ITERATIONS HAVE BEFN MADE AND THF RELEASE IS
LESS THAN MINIMUM RELEASE HALT COMPUTATIONS AS A PROGRAM
CONTROLLED FRROP
IF { NITER . LT. 3 ) GO TO 300
GO TO 1040
CAPACITY ANP RELEASE ARF AT MAXIMUM - SET SPlLL
?80 SPILL = DELTAF
>oo0 COMPUTE RESFRvOIR EVAPORATION i^vAP) AND CHANGE IN BANK STORAGE
( BSTOR ) ( IF APPLICABLE >
300 ELEV = RESEL t CAPX, AREA , K.N, iw )
EVAP = ( ARFA + PAREA (KN) ) » 5-OE-l « CONEVP (*N)
BSTOR = 0.0
IF ( BANKCO (KN) • LE. 0-0 I GO TO 320
BSTOR = ( CAPX - PRECAP (KN) ) ° BANKCO (KN)
REFINE SEED FOR CAPACITY
320 CAPy = TEMCAP + STORAF - RELAF - EvAP - ESTOR
CHECK FOR CONVERGENCE LIMIT OR BALANCED RESERVOIR BUDGET
DIFF = ABS ( CAPX - CAPY )
IF ( DIFF .LE. FPSLON ) GO TO 350
RESERVOIR NOT IN BALANCE RESET SEED START NEW ITERATION-
CHECK NUMBER OF ITERATIONS AGAINST LIMIT
IF ( NITER . GT. 10 ) GO TO 1060
CApX = CAPY
245
-------
SUBROUTINE RESOPR ( QFLOW , KN, pO, IW )
f" O
f « c <* •:>
f o o o
&AQ
f O
-------
SUBROUTINE RESO°R < QF|.OW , KN, ^0, !W )
RDATA (KN, 10) = PDATA (KN, 9 ) * RDATA (IfN, A )
GO TO <»00
RELEASE EXCEEDS MAXIHUM RELEASE SET RELEASE FQUAL TO MAXIMUM
RELEASE AND ENTER ITERATIVE SCHEME
660 RELAF = RELMAX » CONST « DAYS
STORAF = QFLOW (1) - RELAF
SET SEED CAPACITY FOR ITERATIVE SCHEME
CAPX = TEMCAP + EVAP + SPILL * BSTOR
GO TO 200
f«««o THIS IS A RFLEASE REQUIRED FROM THIS PFSERVOlR TO MEFT A
r«««o DOWNSTREAM PEQUIREMENT . CONDITION (2)
70ri RELAF = ARS ( QFuOW (1) )
ro««« CHECK CONSTRAINTS TO DETERMINE IF THIS RfLEASr IS POSSIBLF
IF ( SPILL . GT. o-O > GO TO gQOO
IF ( TEMCAP . LE. CAPMIN ) GO TO 2000
r««<>« RELEASE is POSSIBLE SFT STORAGE AND RELEASE IN ACPE FEET
AND ENTER ITERATIVE SCHEMF
STOPAF = 0.0
GO TO 660
PROGRAM CONTROLLED ERRO" MESSAGES ARE GENERATED HFRE
inOo WRITE ( IW , 1010 ) QFLOW (1)
iniO FORMAT ( iHi, 7X, jTHlNlTIAL INFLOW = , FiQ-O , ?X,
] 21HNfGATTVE INFLOW ERROR )
GO TO 1500
1020 WRITE ! IW, 1030 ) CAP
in3o FORMAT ( IHi , 7X, HHCAPACITY = , FlO.O , 3X,
1 23HNEGATIVE CAPACITY ERROR )
GO TO 1500
1040 WRITE ( IW , 1050) PELcFS
I05o FORMAT ( iHj, 7X, 10HRELFASE = , FlS.l , 3X, 17HANP IS LFSS THAN
1 15H*INIMUM RELEASE )
GO TO 1500
1060 WRITE ( IW, 1070 1 NITER
FORMAT t iHi, 7x, ^AHNUMBER OF ITERATIONS EXCEED LIMIT AND ARE =
1 13 )
GO TO 1500
1500 WRITE ( IW, 1510 ) KNODE , MO, MODE (KN)
151(1 FORMAT < 8X, 12HIN NODE NO. , 13, 0HMONTH = , T3, 5X,
1 12HCALL NO. = , 13, 3X, 18HCOMPUTATION HALTED )
CALL EXIT
?nrio RETURN
END
247
-------
FUNCTION RESCAP ( ELEV, KN. IW )
(-0000 THIS FUNCTION IS USED TO COMPUTE RESERVOIR CAPACITIES AS A SECOND
DEGREE POLYNOMIAL WHERE CAPACITY is A FUNCTTON OF ELEVATION
COMMON X GENRL / LTlTLF(fl), NCTPL(?0, 3), MO*FG, IYRBEG, MOEND, I
IYREND, NUMNOD, isEotzo. 201, KRRI, NUMSEQ(20), NOWTRIPO, ?oj
COMMON / RFSER / RDATA^O, no>, RtiLr(20, 131, c^MFVPiao* i?)% RA
4NKCOt?0)< K*URNM(20* »' » OINFLO(20, 10* I?)-. DEMAMD(20» I0i 16'' *
?DEMNM(20, 10, 8), KQINNM(?0, 10, 8) , KAODNM (20,«)
DIMENSION C* (10)
EO 10 K = 1, 5
CA (K.) = 0.0
10 CONTINUE
CO 30 K = 2*, 97. 1
KX = K
TEST = R.DATA ! ifN, K )
IF ( TEST .FQ. 0.0 ) GO TO 60
IF ( ELEV . LT. TEST ) GO TO 60
30 CONTINUE
NODF = NCTRL (KN, 1 )
WRITE ( IW, 50 ) NODE . FLEV
50 FORMAT ( 1H1, 7X, 41HTHE EQUATION ARRAY FOR RESERVOIR IN NODF
1 13, ^OH WAS EXCEEHEn FOR ELEVAT^N = ,F10.2 )
CALL EXIT
60 KW = KX - 4
EO 80 K = 1, 3
KY = KW * K
CA (K) = RDATA ( K.N, KY)
SO CONTINUE
x = ELEV
X = ELEV - PDATA ( KN, KW )
RESCAP = FOFx ( CA , X« 1 >
RETURN
END
248
-------
FUNCTION RESEL ( CAP , AREA , «N, IW )
C«o«o THIS FUNCTION IS USFD TO COMPUTE ELEVATION AND ARFA FOR
r«o«« A RESERVOIR ON THE BASIS OF A SECOND DEGREE POLYNOMIAL WHERE
f»o»o CAPACITY IS A FUNCTION OF ELEVATION AND AREA Is THE FIRST
C»o«o DERIVATIVE OF THE FUNCTION
COMMON / GPNRL / LTITLFU), NCTRL(20, 3), VOBEG, IYRBEG, MOEND, I
AYREND, NUMNOD, isEo<2o, 201, KPRI, NUMSEQ^O), NDWTRIZO, 201
COMMON / RFSFR / RDATA120, 110), RULM20, 13), C(20, 1?), RA
ANKCO(?0>, KSURNM120, fl), OINFLOI20, 10, 12), DEMAMDI20, 10, 16U K
PDEMNM(20, 10, 8), KQINNM(20, 10, 8) , KAGUNM (20,»)
ElMENSION CF dO)
DO 50 K = 2^, 98,
WRITE ( Iw, 60 ) NODE « CAP
60 FORMAT < IHI, ?x, ^IHTHE EQUATION ARRAY FOR RFSERVOIR IN NODE
1 13, 30H WAS FXCEFDEU FOR CAPACITY OF , F10.0 )
CALL EXIT
70 KW = KX - 1
CO 90 K = It 3
KY = KW - K
CF (K) = RDATA ( KN, KY )
90 CONTINUE
CF (3) = CF (3) - CAP
KOLE = 6REL<"AP
XR = ROOTWO ( KN, 1, CF, MOLE , IW )
AREA = t 2.r' « CF (1) » XR ) + CF (?)
KR = KX - 5
RESEL = RPATA ( KN, KR ) + XR
RETURN
END
249
-------
SUBROUTINE PWROPR (ACFS, jAP, IMON, DAYS, ElEVl, SPILL* IW )
cocoo THIS SUBROUTINE SIMULATES THE OPERATION PF A POWEP PLANT WITH ANY
(-0000 NUMBER OF UNITS ( GENERATORS ) WHERE EACH AND EVEPY UNIT MUST
(-00*0 HAVE IDENTICAL <"Ap ABILITIES AND CONSTRAINTS . THE ROUTINE
Coooo [ESIGNED ASSUMES UNITS fiRp ALWAYS OPERATING AT pULL GATE.
Cocoa THE ROUTINE CAN EASILY BE MODIFIED FCR OPERATION AT BEST RATE OR
Coooo A COMBINATION OF BOTH GATE SETTINGS
Cocoo ALTHOUGH THIS SUBROUTINE SIMULATES ONE UNIT IN A POWER PLANT ANP
roooo THE NUMBER OF HOURS OF OPERATION ARE IDENTICAL , o» IT is NOT
<-oooo NECESSARILY ASSUMED THAT THEY BE THE SAME HOURS U' TIME
COMMON / GFNRL / LTITLE<8), NCTRL<20, 3), MQBEG, IYRBEG, MOEND, I
AYREND, NUMNOD, TSFQ(20, 20), KpRI, NUMSEQI20), NDWTR120, 20)
COMMON / RpSER / RDATA(20, 110), RULr(20, 13), CONFVP<20, 1?>, BA
ANKC0120), KSURNMI20, P.), OINFLOI20, 10, 12), ^FMAMDtPO, 10, 16>» *
BEEMNM120, in, 8), K3INNMI20, 10, 8) , KAQUNM (20, 8)
COMMON / pnWFR / CApBAB(20< 20), TDATA(20, 4), ENERGYI20, 8), PDA
ATA(20, 12)
DIMENSION CA (10) , MODE (20)
DATA (HPTOKw = T.i(57E - 1), (DENS = 6.24E1), (HPF = 5-5E2), (DAF
A= 4.356E4) , (IFIR = 0 )
CHECK TO DETERMINE IF THIS IS THE FIRST CALL TO THIS SUBROUTINE
FOR THIS JOB
IF ( IFIR .NE. 0 ) GO TO 20
INITIALIZE CALL INDEX FOR EACH POWER PLANT
DO 10 K = 1, 20
MODE (K) = 0,
IQ CONTINUE
IFIR = 1
20 KP = JAP
Coo«o CHECK TO DETERMINE IF THIS IS THE FIRST CALL FOR PARTICULAR POWER
Coooo PLANT
IF ( MODE (KP) . NE. 0 ) GO TO 70
Cooo« COMPUTE GENERATOR CAPABILITY ( GCNCAP = GENERATOR RATING (KV-A) o
POWER FACTOP o OVERLOAD (PERCENT) / 100.0 > IN »
Coooo COEFFICIENTS ARE ARRANGED IN DECREASING POWERS OF THE
Coooo INDEPENDENT VARIABLE »« POLYNOMIAL IS EITHER SFCOND DEGREE OR
Coooo LINFAR
CA (1) = ENFRGY ( KP,7 )
CA (2) = ENERGY < KP,S >
CA (3) = ENFRGY ( KP,5 )
Coooo CHECK TO DETERMINE IF POLYNOMIAL IS LINEAR
IF ( CA (1) . EQ. 0.0 ) GO TO 30
coooo POLYNOMIAL IS NONLINEAR
250
-------
SUBROUTINE PWROPR (ACFS, jAP,
PAYS, FLEV1, SPILL» I* )
MOLE = 6RPWROPR
THEAD = ROOfwO ( KP, 1 , CA, MOLE, ! v, )
GO TO 40
C«o POLYNOMIAL TS LINEAR
30 THEAD = ( HPMAX - CAC3) ) / CA(2)
(•*«»« COMPUTE MAXIMUM DISCHARGE ( TCFS ) FOR ONE UNIT AT MAXIMUM
<•««»» HORSEPOWER USING POLYNOMIAL WHERE DISCHARGE IS A FUNCTION OF
C THE COMPUTATION OF HORSEPOWER
40 DO 50 K = 1, 10
CA (K) = 0.0
50 CONTINUE
DO 60 K = 1, 4
CA (K) = ENFRGY ( KP, K)
60 CONTINUE
TCFS = FOFX ( CA , THEAD , 0 )
COMPUTE CONSTANT EFFICIENCY FOR HEAHS GREATER THAN HEAD AT
MAXIMUM DISCHARGE
TEFF = HPMAX » HPF / ( TCFS * THEAD « DENS )
STORE THESE ONE TIME COMPUTED VALUES IN ARRAY CAPPAB
CAPBAP ( KP» 9 ) = HPMAX
CAPBAB ( KP, 10 ) = THEAD
CAPBAP ( KP, 11 ) = TCFS
CApBAB ( KP, 12 ) = TEFF
BEGIN SIMULATION
70 TEST = ELEV1 - TDATA ( KP.
HDMIN = CAP*AB ( KP, ft )
HDMAX = CApBAB ( Kp, 7 )
HPMAX = CAPBAP ( KP* 9 )
THEAD = CAP"AB (KP, 10 )
TCFS = CAPBAB ( KP, 11 )
TEFF = CAPBAB (KP, 12 )
CONST = DAYS « 2.4E1 « ^.6E3 / "AF
CLEAR VARIABLES CONTAINED IN ARRAY CAPBAP THAT ARF TO BE COMPUTED
C«««o DURING THIS CALL
DO 80 K = 13, 20
CApBAB
-------
SUBROUTINE PWPQPR ao LESS THAN OR GRFATpR THAN HFAD AT MAXIMUM HORSE P|">WER
C«o«o GROSS HEAD = RESERVOIR WATER SURFACE ELEVATION - TAIL WATER
KSET = 0
GHEAD = ELEV1 - TWEL
IF ( GHEAD . GT. THEAD ) GO TO 120
GROSS HEAD TS LESS THAN HEAD AT MAXIMUM HORSEPOWER - USING
POLYNOMIAL COMPUTE DISCHARGE AT THIS HEAD
DO 110 K = 1 , <>° COMPARE GROSS HFAD TO MINIMUM HEAD
:iAO JF ( GHEAD . LT. HDMjN ) G" TO 1000
IF ( KSET . NE. 0 ) GO TO 1&0
ro HORSEPOWER
TEFF = THP / THFORY
GO TO 170
160 THP = HPMAX
r««x>o CONVERT HORSEPOWER TO KILOWATTS
170 TKW = THP » HPTOKW
C.««»» COMPUTE TOTAL POWER IN KILOWATTS ANfl TOTAL DISCHARGE IN CFS FOR
POWER PLANT WHERE EACH yNlT IN PLANT IS ASSUMED ITFNTICAL IN EVERY
RESPECT
TOTKw = TKW « UNITS
252
-------
SUBROUTINE PWROPR (ACF?, JAP, IHON, DAYS, ELEvl, SPILL, IW )
C«
COMPUTE HOURS OP OPERATION PER UNIT AS THE RATio np TOTCFS TO
WATER RELEASED FROM RESERVOIR ( ACFS IN ARGUMENT LIST )
RATIO = ACFS / TOTCFS
IF ( RATIO . GT. 1.0 ) RATIO = 1.0
HOURS = RATIO « 2.4E1 « DAYS
COMPUTE POWER IN KILOWATT-HOURS FOR WHOLE POWER PLANT
KWH = TOTKW * HOURS
COMPUTE TOTAL FLOW THROUGH POWER PLANT IN ACRE FEFT
RELAF = TOTCFS » CONST « DAYS
COMPARE POw^R GENERATED WITH POwER REQUIREMENT { FIRM DEMAND)
FOR THIS TIME FRAME
RKWH = PDATA i KP, IMONJ
IF ( KWH . LT. RKWH ) GO TO 1040
POWER GENERATED IS GREATER THAN FIRM REQUIREMENT «« COMPUTE
SECONDARY POWER GENERATED
KWHSEC = RKWH - KWH
KwHSEC = KwH - PKWH
SET VARlABLFs IN CAPBAB ARRAY
CAPBAP ( KP, 13 ) = HOURS
CAPBAB t KP, 14 ) = KWH
CAPBAB KP, 15 ) = RELAF
CAPBAB KP, 16 ) = TWEL
CAPBAB KP, 17 ) = GHEAD
CAPBAB KP, 16 ) = RKWH
CAPBAB KP, 19 ) = KWHSEC
CApBAB Kp, 20 ) = TEFF
GO TO 2000
WRITE ( IW
FORMAT (
KNODE = NCTPL
WRITE ( IW
1030 FORMAT ( 8
CALL EXIT
WRITE ( IW , 1050 )
FORMAT t IHI, ?x, 40HP0wER GENERATED is LESS THAN FIRM DEMAND
GO TO 1020
RETURN
END
inoo
inlo
1020
1010 )
, 7X, 36HGROSS HEAD IS LESS THAN MINIMUM HEAD )
KP * 1 )
1030 ) KNODE ,
11HFOR NODE =
IMON
, I3» 5X, 1&HFOP MONTH NO. =
13 I
253
-------
SUBROUTINE DTABLE ( JK, KFA, 1C, KS )
COMMON / GENRL / LTITLE(S), NCTPH20, 3), MOBEG. IYR6EG. MOEND, 1
AYREND, NUMNOD, ISEO(20. 20), KPRI, NUMSEQI20), NDWTR(?0, 20)
COMMON / RESER / RDATA(?0, 110), RULE120, 13), CONEVP(20, 12), BA
ANKCO« DEFICIT
1C = 0
GO TO 65
A RESERVOIR EXISTS IN THIS NODE AND SHOULD BE INCLUDED TN RANKING
60 1C = 1
IN THIS MODFL RESERVOIR CAPACITY = AO + AI o H + A2 « H «» 2 ,
THE FIRST DERIVATIVE DC/DH= Alf2»A2oH = RFSERVQlR ARFA . LET
C»o«« DC/DH = AREA ANP RANK WITH RESPECT TO AREA. A»EA = RDATA ( JK, 5 )
(-0000 JN THE RDATA ARRAY
AREA (1C ) = RDATA ( JK, 5 )
KFA ( 1C ) = NCTRt ( JK, 1 )
roooo fcUH THE USF OF TH£ ARRAY NDWTR WHICH CONTAINS THF NODE SYSTEM
(-0000 STRUCTURE SFARCH FOR UPSTREAM RESERVOIRS THAT CAN SERVE THIS NODF
65 DO 70 K = 1, 20
IF ( NDWTR ( JK, K) . LF. 0 ) GC To 80
NODE = NDWTR i JK, K>
KN = NFIND ( NODE )
CHECK TO DETERMINE IF A RESERVOIR EXISTS IN NODE
IF ( RDATA ( KN, 3 ) . EQ. 0.0 ) GO TO 70
1C = 1C •»• 1
AREA (1C ) = PDATA (KN, 5 )
KFA (1C ) = NDWTR ( JK, K)
70 CONTINUE
RANK UPSTREAM RESERVOIRS WITH RESPECT TO APFA
80 MAX = 1C
MlN = 0
90 MlN = MlN + 1
IF ( MlN . PE. MAX ) GO TO 2000
DO 100 K = MJN , MAX
IF ( AREA (MIN ) . GE. AREA (K) ) GO TO 100
254
-------
SUBROUTINE UTABLE ( JK, KFA» 1C, KS )
I TEMP s UFA ( MIN )
AREA (MlN ) = AREA (K)
UFA (MIN ) = Kp6
-------
FUNCTION NHND (NODE)
COMMON / GFNRL / LTITLEU)* NCTPL
-------
SUBROUTINE RFCOMP ( QFLOW, NK, IP, INDY, MO, CONST, IW !
COMMON / PROPER / SOlL(2f), 10, 21), BTANMPO. 10), TTANKI20, 10),
1 GRTANKI20, 10) , NSEG (20) , CONUSE (20, 10, 30)
2 , KUSE (20. 10, 30)
COMMON / RFSER / RDATA120, 110), RULF(20, 13), CONEVP<20, 12), BA
10, 12), DEMAND<20» 10, 16), KsY
»o ARE THOSE CONSIDERED IDEAL OR FULL SUPPLY COMD'JTEP AT THE
Coooo POINT OF DIVERSION ANP INCLUDE CONVEYANCE LOSSES AND FARM LOSSES
C«o«« CONSUMPTIVE USE VALUES REPRESENT A LuMPED PF.NEFICIAL AND NON
BENEFICIAL REQUIREMENT. QFLOW = VOLUME ANH CHFMISTRY OF
SOURCE, NK AND IP ARE THE NODE AND SEQUENCE INDICES RESPECTIVELY
INDY = USE INDEX , MO = MONTH ANC IW = WRITE UNIT.
IF ( INDY .PT. 5 ) GO TO 1010
GO TO { 8Q» 200, 300» 300, 300) , INDY
DEMAND HAS PEEN MET SET SHORTAGE I PERCENT ) EQUAL TO ZERO
*0 CONUSE ( NK, IP, 15 ) = 0.0
teo»« BEGIN TO COMPUTE RETURN FLOW ITS ALLOCATION ANn CHEMISTRY
RFLOw (1) = CONST - CONUSE
PER
} =
) QINFLOI JK, KSEQ, 11 >
1.0
?.o
IF ( KDES . EQ. 3RGWV ) QlNFLO< JK,
11 )
257
-------
SUBROUTINE R.FCOMP ( QFLOW, NK, IP, INDY, MO, CONST, IW )
IRRIGATION OR A MUNICIPAL AND INDUSTRIAL USE
IF ( DEMAND ( NK, IP, 14) . EQ. 1.0 ) KUSE ( NK, IP, 17 ) = 3RIRR
IF ( DEMAND < NK, IP, 14) . EQ. 2.0 ) KUSE ( NK, IP, l~f > = 3RMAI
GO TO 120
DEMAND HAS WOT PEEN MET COMPUTE SHORTAGE IN PERCENT
?00 RATIO = CONUSE < NK, IP, 13 ) « l.OF-2
TEMCU = CONST *> ( i.o - RATIO )
SHORT = ( i.o - ( TEMCU / CONUSF INK, IP, MO») <» i.op2
C«oo« COMPARE SHORTAGE WITH SHORTAGE INDEX
IF ( SHORT . GT. CONUSE < NK, IP, 14 )) GO TO 1000
C«o»« SHORTAGE IS LESS THAN SHORTAGE INDEX COMPUTE RETUPN FLOW
RFLOW (1) = CONST - TEMCU
CONUSE ( NK, IP, 15) = SHORT
GO TO 90
WHEN THE INPEX ( INDY ) IS GREATER THAN 2 THE ROUTINE TpFATS
THE DEMANDS AS TRANSFERS OF FLOW WjTHOUT RETURN FLOW OR
CONSUMPTIVE USE
300 KDES = KUSE (NK, IP, 19)
KNODE = KUSF CNK, IP, 20!
KSEQ = KUSE ( NK, jP, 21 )
JP = NFIND ( KNODF)
IF ( CONST . GT. o-O > 6° T° 315
EO 310 K = 1, 10
QlNFLO< Jp» KSEQ, K 1 = 0.0
310 CONTINUE
GO TO 33o
315 GlNFLOl JP, KSEO, 1 ) = CONST
DO 320 K = 2, 10
QINFLO ( JP, KSFQ, K ) = QFLOW (K)
32o CONTINUE
330 IF ( KDES . EQ. 3RSUR ) QlNFLOC JK, KSEQ, 11 ) = 1.0
IF t KDES . £Q. 3RGWR j QlNFLOl JK, KSEQ, 11 > = 7.0
IF ( KDES . EQ. 3RGWV ) QINFL01 JK, KSFQ, 11 ) = .t.O
GO TO 2000
inOQ WRITE ( IW, 1005) (KDEMNM ( NK, IP, K) , K= 1, 8), NK, IP
1005 FORMAT f 1H1, 7X, 32HSHORTAGE INDEX IS VlOLATFD FOR , RA8 /
1 13HFOR NODE NO. , 13, 5X, 12HSEGMENT NO. . 13 )
CALL EXIT
1010 WRITE ( iw, 1015 ) INDY , NK, IP
inl5 FORMAT < 1H1, 7HINDY = , 13, 5X, 24HWHICH IS GREATER THAN 5 ,
1 13HpOR NODE NO. , 13, 5X« IZHSEGM^NT NO. , 13 )
CALL EXIT
?nno RETURN
END
258
-------
SUBROUTINE CHEMGR ( QFLOW , CONST , KM, KS !
COMMON / RFSER / RDATA(20. 110), RULES20, 13>, CONEVPt20, 12). 8A
ANKCO<20), KSURNM<20, 8), QlNFLO{20, 10, 12), DEMAND<20, 10% 16), K
BDEMNM(20, 10, 8), KQINNM(20, 10, 8) , KAQUNM (20,*)
COMMON / PROPER / SOIL(20, 10. 21), PTANK(20, 10), TTANK(?0, 10),
1 GRTANK<20, 10) , NSEG (20> , CCNUSE (?0, 10, 30)
? , KUSE (20, 10, 30)
EQUIVALENCE ( CnNUSE, KUSF >
DIMENSION QFLOW (10)
o THIS SUBRnuTINE IS USED TO COMPUTE A LINEAR Mix OF CHEMISTRY
C»o»« FOR FLOW ADnED TO THE RIVFR ( SURFACE) , TO THE AOUIFER (SUBSURFACF)
, To THE PODDING POOL ABQVE SOIL COLUMN, ANTl FLQW ADDED To THE
RESERVOIR.
QFLOW = FLOW IN RIVER ARRAY , KN AND KS ART INDICFS FOR NODF AND
SEQUENCE RESPECTIVELY, CONST = DEMAND OR PlyFRslON
r«o»o THIS ENTRY MIXE5 RIVER FLOW WITH AQUIFER
VOLA = CONST
VOLB = GRTANK (KN, 1 )
VOLUME = VOLA + VOLB
DO 5U K = 2, 10
GRTANK (KN, K ) = I QFLOW (K) « VOLA * GRTANK (KN, K ) « VOLB )
1 / VOLUME
50 CONTINUE
GO TO 2000
roc«o THIS ENTRY MIXES RIVER FLOW WITH PONCING POOL
ENTRY CHEMTT
VOLA = CONST
VOLB = TTANK (KN, 1 )
VOLUME = VOLA + VOLB
DO 60 K = 2, 10
TTANK (KN, K ) = ( QFLOW (K) « VOLA + TTANK (KN, K ) « VOLB )
1 / VOLUME
*0 CONTINUE
GO TO 2000
C»«« THIS ENTRY MTXES FLOW TO RIVER WITH RIVER
ENTRY CHEMRV
VOLA s QINFLOt KN, KS, 1 >
VOLB = QFLOW (1)
VOLUME = VOLA + VOLB
2?Lo! I*? =',QlSFLO(KN,KS,K) » VOLA * QFLOW (K) « VOLB ) / VOLUME
70 CONTINUE
r«M,.. 5S5S°EN?R? MjXES RESERVOIR CHEMISTRY WjTH RIVER CHEMISTRY FOR
C«o«« DEFICIT RELFASE
ENTRY CHEMRS
VOLA = QFLOW (1)
VOLB = CONST
VOLUMF = VOLA * VOLP
DO go K = ?» 10
259
-------
SUBROUTINE CHEMGP ( QFLOW , CONST , KN, KS )
QFLCK (K) = ( QFLOW (K) « yOLA * RDATA ( KM, KA) « vOLB ) / VOL|JMF
'IXES PIVER FLOW ( AS INFLOW TO RESERVOIR) WITH
C«««G RESERVOIR CHEMISTRY
ENTRY CHEHRR
VOLA = QFLOW (1)
VOLR = RDATA ( KN, 3 )
VOLUME = VOLA + VOLR
DO 110 K = ?, 9
KA = Y. * 99
RDATA ( KN, KA ) =(QFLOW (K) « VOLA + RDATA (
-------
SUBROUTINE CONVfR ( TEMP i KR )
THIS SUBROUTINE IS USED TO CONVERT INPUT UNITS TO PPM AND
ALSO TO CHECK BALANCE BETWEEN CATIONS AND ANIONS
DIMENSION TFMP (15) , EQVWT (8)
DATA ((EQVWMio , K=l,8) = 20.04,12.1 6*23.,35.46,4?.03,30.»61. 01,
1 62.008 )
IF ( KR . Nc. 1 ) GO Tn 80
UNITS ARE IN MILL I FQuI VALENTS PER LITER CHECK SUM OF CATION*
AND ANIONS
SUM CATIONS
10 SUMC = 0.0
DO 20 K = 2« 4
SUMC = SUMC + TFMP (K)
20 CONTINUE
SUM ANIONS
SUMA = 0.0
DO 30 K = 5»9
SUMA = SUMA + T^MP (K)
30 CONTINUE
DIFF = SUMC - SUMA
IF ( DIFF . GT. 0.0 1 GO TO 40
TEMP (2) = TEMp (?) + ABS ( DIFF )
GO TO 45
40 TEMP (6) = T£MP(6! + DIFF
CONVERT FRO" MILL I EQUIVALENTS PER LITER TO PARTS PER MILLION
45 DO 50 K = 2, 9
KA = K - 1
TEMP (K) = T£*p (K) « EQVWT ( KA )
50 CONTINUE
COMPUTE TOOL DISSOLVED SOLIDS FROM ANALYSIS
TDS = 0.0
DO 70 K = 2* 9
IF ( K • EO. 7 ) GO TO 60
TDS = TDS + TEMP
-------
SUBROUTINE VERNAL ( KN, IS* MO. OBSERV * QFLOW )
f-oooo THIS SUBROUTINE SIMULATES THE ASSUMED CRITERIA FOP EACH NODE
(-0000 IN OBTAINING A HYDROLOGIC BALANCE WITH SIMPLE TRANSFERS FROM
foooo SURFACE TO SUBSURFACE OR THF CONVERSE. THE ARGUMENT LIST IS
(-0000 AS FOLLOWS - KN AND IS ARE THE NODE AND SEQUENCE INDICES
Coooo RESPECTIVELY. 0«SF.VR IS THE ARRAY CONTAINING THE OBSERVED OUTFLOWS
coooo QFLOW IS THE ARPAY CONTAINING THE RlvFR VOLUMF ANn CHEMISTRY
COMMON / RFSER / RDATA(20, 110), RULF(?0, 13), CONEV»(20, 12), BA
ANKCO<20>, KSURNM(20, 8), QlNFLO(20, 10* 121* OEMAND<20* 10* 16)* K
PDEMNM(20, 10* 8), KQINNMI20, 10, 8) * KAQUNM (20,P)
DIMENSION OPSERV (10) . QFLOW (1C)
GO TO ( 50, 70, 90 ) , KN
c««oo THESE ARE THE TPANSFER PARAMETERS FOR NODE 101
5n MIN = 3
MAX = 7
KA = 4
GO TO 110
roooo THESE ARE THE TPANSFER PARAMETERS FOR NODF 10?
70 MIN = 2
MAX = 6
KA = 3
GO TO 110
C«ooo THESE ARE THE TRANSFER PARAMETERS FOR NODE 10?
90 MIN = 2
MAX = 4
KA = 3
IF ( IS .GT. MIN )
110
130
150
GO TO 150
DO 130 K = 1, 10
CBSERV (K) = QINFLO ( KN, IS, K)
CONTINUE
GO TO 2000
VOLA = OBSEPV (1)
VOLB = QINFLO (KN, IS, 1)
170
VOLUME = VOLA *
DO 170 K = ?, 10
OBSERV
-------
SUBROUTINE VERNAL ( KN, IS, MO, OfiSERV , QFLOW )
DEMAND (KN, KB, MO ) = ABs (DIFF)
?000 RETURN
END
263
-------
FUNCTION ROOTWO ( KNODF i KSEG , CA , MOLF , Iw )
fOOOO
KNOPE = NODF
50
ROOT OF (X) USING Qt)ADRATlC WHERE
X = ( - B + SQRT ( P»B _ 4 « A « C ) / 2 A
INDEX AND KSEG r SEGMENT INDEX CA = ARRAy OF COEFFICIENTS
IN DECREASING POWERS OF (x) MOLE = ALPHA DFFINITION OF (x)
DIMENSION CA (10)
RAD = CA (2) » fA (2) - ( 4.0 * CA (1) « CA (3) )
IF ( RAD . LE. 0.0 ) GO TO 1000
ROOTWO = ( - CA (2) f SORT ( RAD )) /(2.0 » CA (1))
V.RITE ( iw , so > MOLE , CAUI , CAm * CAOJ , PAD , POOTWO
FORMAT ( 1H1, flX, R6, 2X,
» F15.8, 2X,
9HCA (1)
9»CA (2)
9HCA (3)
9HRAHCAL
9HROOTWO
F15.8,
F15.8,
F15.8
F15.8
9X.
GO TO 2000
inOQ WRITE
?0A0
( IW , 1005
FORMAT(1H1» 8X,
1 5H ZFRO
? 4HFOP
RETURN
END
KSEG i MCLE
KNoDE
25HPROBABLY NOT A REAL ROOT , 5X,
, 5X, 7HNODE = , 13 , 5X, 6HSEG =
R6, 2X, 11HCOMPUTATION )
11HSET ROOT TO
, 13, 5X,
264
-------
FUNCTION Ron? (CONST, XMAX, XMIN 1 MOLE , CA, TF.STA )
C»«M>« THIS FUNCTION WILL COMPUTF THE REAL ROOT IF ONE FXISTS
f.«««« IN THE RANGF XMIN - XMAX INCLUSIVE. IF A ROOT EXISTS IN THIS RANGE
NOLF WILL BF SET = TRUE, IF NO ROOT EXISTS , HOLF WILL BF SET =
FALSE. F (x) = CONST AND COEFFICIENTS CA (1) » CA (N + l ! ARE
TRANSFERRED THROUGH ARGUMENT LIST IN INCREASING POWERS OF THE
INDEPENDENT VARlABLEi IN THE FUNCTION FOFX NK = « , SO IF A
A HIGHER THAN FOURTH DEGRFE POLYNOMIAL IS TO SOLvFD RESET NK .
IF NO ROOT TS FOUND AFTFR TWENTY ITERATIONS NO ROOT WILL RE
CONSIDERED *NP HOLE WILL BE SET = NONE.
DIMENSION CAUO)
XA s i.E-1 » XMftX
0.0 ) XA = l.E-1 «
f
C
C
C
C
C
C
IF ( XMIN . LT
XB = XA + XA
IW * 61
WRITE < IW , 5 ) XMAX , XMIN
5 FORMAT t 1H1, 8X, 7HX»
-------
FUNCTION ROOT (CONST, XMAX, XMIN t MOLE t CA, TESTA )
END
266
-------
FUNCTION Fo^XtCF^, (CAM )
C«««<> THIS FUNCTION EVALUATES ANY POLYNOMIAL FROM 1ST nEGRFE TO
<•«««« FOURTH DEGRFE INCLUSIVE 1 SEE SETTING OF NIC ) IF A HIGHER DFRREf
C«««<> POLYNOMIAL tS TO BE EVALUATED RESET NK WHERE NK - 1 = DEGREE
DIMENSION (T(10)
FOFX = 0.0
NK = *>
?0 FOFX = FOFX * CF(NK)
FOFX = FOFX « X
NIC s NK - 1
IFINK.GT.1) GO TO 30
FOFX = FOFX * CF(l)
IF t KAM . NE. o ) GO TO ?ooo
BO <»0 < - 1» 10
CF (K) = 0.0
4o CONTINUE
70^0 RETURN
END
267
-------
SUBROUTINE SOILCO < KNODE , IW , IR , NSFGS)
f O O O O
r«ooo THIS SUBROUTINE IS USED A? A MONITOR SUBROUTINE FOR ALL SUB-
(-ofioo ROUTINES USFD IN THE SIMULATION OF PERCOLATING APPLIED WATERS
r«o VERTICALLY THROUGH A SOIL COLUMN USING PISTON DISPLACEMENT OF
t»o«o SOLUTION FROM SEGMENT TO SEGMENT. THE FIRST STEP IS TO CHECK
r»o«« FOR FIRST CALL TO THIS SUBROUTINE — IF FIRST CALL CLFAR ARRAY
foooo ITER WHERE TALLY IS KEPT FOR NUMBER CF TIMES SEGMFNT HAS *EEN
tcooo EQUILIBRATE? FOR EACH APPLICATION OF WATERS TO BF PERCOLATED
Coo<><> AND ALSO ARRAY KwATER WHERE TALLY IS KEPT 0F EACH APPLICATION
f.»o«o OF NEW WATER.
COMMON /SOLi'TF / SLCA, SLMG, SLNA, SLCL , SLS04, ?LHC03, SLC03,
] SLN03 ,UCAS04 » UMGS04
COMMON /SOLTD/ CATEX, GYPSUM, SLIME
-------
SUBROUTINE SOILCO ( KNODE , IW , IR , NSEGS)
PRATIO = WRATIO o 1.E2
IF ! KWATER (KNODE ) . GT. 0 ) GO TO 50
C««»IL AM.UYSIS
CALL SUBROUTINE TO EQUILIBRATE INITIAL SOIL ANALYSIS FC'P EACH SEGMENT
CALL FRSTIM 1 KNODE , K, I w >
C CALL SCRIBE (IW)
CALL SUBROUTINE TO CALCULATE SOLUTION - PRECIPITATION Or GYPSUM <>'•>
CAS04 . 2H20
50 CALL GYPEX ( KNODE , K, IW , PRATIO )
f CALL SCRIBE UW)
C CALL SCRIBE (IW)
CA (I) = SLCA
CALL SUBROUTINE To CALCULATE EXCHANGE C? CALCIUM AfjH So^lUM
CALL F.XCANA ( KNODE i K, IW , PRAKIO )
C CALL SCRIBE (IW)
CA (2) = SLCA
CALL SUBROUTINE Tn CALCULATE EXCHANGE isf CAICIUH A^D MAGNESIUM
CALL EXCAMG ( KNODE , K, IW , PRAT1C )
C CALL SCRIBE (IW)
CA (3) = SLCA
CALL SUBROUTINE TO CALCULATE REACTION Of7 CALCIUi-i WIYn ('MU)O,,AIK -
C,a OHIL -
C SYSTEM USINR CONCENTRATION OF CAICIUM IN SOLUTION »5 1ME
CONVERGENT VARIABLE
C WRITE ( IW, 60 ) NCT , ( TA (M) , M - 1, s )
c 60 FORMAT ( IHI, ex, IOHCOUNTFR = , iir , / f ( ?
C WRITE ( IW . 60 ) NCT, (CA (M) , M = 1, 4 )
r 80 FORMAT ( // 8X« 10HCOUNTER = , 12 / A ( 8X1 fl*.*) /
r 1 !7HENDING CONDITIONS / )
C CALL SCRIBE (IW)
?00 CONTINUE
THE SOIL COLUMN HAS BEEN EQUILIBRATED FOR THIS PARTICULAR
ITERATION oo TEST FIRST TO DETERMINE IF IrilS ITERATION wAs FOR
THE EQUILIBRATION OF THE ORIGINAL SOIL ANALYSIS
IF ( KWATER ( KNODE ) . EO. 0 ) GO TO 300
269
-------
SUBROUTINE SOILCO ( KNODE , IW , IR , NS?GS)
Cooflo CF SURFACE WATERS
GO TO 310
r«e*« MOVE APPLIED WATER INTO FIRST SEGMENT ANT EQiJILI PR ATE
300 MN = 0
MAX = NSEGS
r«o«o UPDATE COUNTER OF NUMBER OF APPLICATIONS OF SURFACE WATERS
KWATER ( KNOBF ) = KWAT^R (KNODE ) + 1
ISW = 3
310 ("IN = MIN + I
IF ( MIN . GT. NSEGS ) GO TO 350
KV = MIN
CALL SMOVER ( KNODE , KV , NSEGS , ISW )
GO TO 35
r«o«« THIS IS THE LAST EQUILIBRATION FOR THIS SOIL COLUMN FOR THIS TIME
C«o»w FRAME »» MIX SOLUTION iN LOWEREMCST SEGMENT WITH PTANK o°«« SET
roo»-r, VOLUME OF BTANK
350 ISw = 4
CALL SMOVER ( KNODE , KV , NSEGS » ISW )
GO TO 2000
?000 RETURN
END
270
-------
SUBROUTINE
i KN , KV, NSEGS,
(•.««*>« THIS SUBROUTINE HAS MULTIPLE ENTRIES AND CONTROLS THE TRANSFER
C«« OF DATA TO AND FROM THF COMPUTATIONAL VARIABLE NAME SEQUENCE AS
C<>»<><> THIS DATA RELATES TO THAT STORED IN THE ARRAY TITLFD SOIL. THE
C.«» SUBROUTINE ALSO SIMULATES THE VERTICAL DISPLACEMENT OF SOLUTION
FROM ONE SEGMENT TO THE NFxT LOWER SEGMENT (PlSTC* EFFECT). KN AND
KV IN THE ARGUMENT LIST RFpRESENT THE NODE INDEX AND SEGMENT INDEX
RESPECTIVELV.
COMMON/PROPER / SOIL (20,10.21), ETANK , TTANK (20,10),
1 GRTANK (?0,10)
COMMON /SOLI'TE / SLCA, SLMG, SLMA, SLCL , SLS"4, CL^C03' ?LCQ3«
1 SLN03 ,UCAS04 , UMGS04
COMMON /SOLID/ CATEX, GY^snM, SLIME ,PDFN'S, EXCA ,EXNA, EXMG ,
i FXTRAC , DEPTH, SEGVOL , INTER
DATA ((EQVWT(K), K = l,8) = 20.4 ,12 , 16,23. ,35 .46,48.03,30.,61.01,
1 f>2.0n8 )
DIMENSION EQVWT (fl) , VALNCE (8)
RFAL INT^R
DATA( (VALNCF (M, K =1,8) = 2=, 2-, 1 • , 1 • , 2 • ,1 • , 2• , 1 • )
THE MAIN ENTRY TO THIS SUBROUTINE CONVERTS I>N!TS OF INITIAL SOIL
ANALYIS INTO COMPUTAT IONAL UNITS
CONVERSION OF MILlIFQUIVALENTS PER LITER (MFQ/L) TC COMPUTATIONAL UNITS
OF MOLS PEP LITFR (SOLUTION)- MOLS/LITER =(PE3/L) /(VALENCE * 1.E3)
IF ( ISW . ME. 0 ) Go TO 70
DO 60 K = 1,6
SOIL (KN,KV,K) = SOIL (KN,KV,K) / (VALNCE (K) <> 1.E3!
CONTINUE
MILL EQUIVALENTS PER too GRAMS TO COMPUTATION UMITS or MOLS/GRAM
WHEN jSW = 0 CONVERSION IS NECESSARY FOH INJTIAL SOIL ANALYSIS
PER GRAM (SOLID)- KOLS/r-RAM = (MEO/ ] OOGRAMS) /(VALENCE <> 1.E5)
SOIL (KN,KV,11) = SOIL (KN,KV,11)/ (VALNCE<4)<> l.FS)
SOIL (KN,KV,12) = SOIL (KN,KV,12)/ (VALNCE15)0 l.FS)
WHEN jSW NOT F.GUAL ZERO COMPUTATIONAL U^lTS RETAINED NO
CONVERSION REQUIRED
IN COMPUTATIONAL VARIABLF NAME SEQUENCE
50
7(1
PLACE DATA OF
SLCA
SLMG
SLNA
SLCL
SLS04
SLHC03
SLC03
SLN03
EXTRAC
SEGvOL
CATEX
GYPSUM
SLIME
BDENS
DEPTH
=
s
—
s
s
s
s
:
s
2
s
s
z
-
-
SO'L
SOIL
SOTL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOTL
SOIL
SOIL
SOIL
SOIL ARRAY
(KN,
(KN.
(KN.
(KN.
(KN,
(KN.
(if N.
(KN.
(KN.
(KN.
(KN.
(KN.
(KN.
(KN.
,KV,
iKV,
,KV,
,KV,
,KV,
,KV,
iKV,
,KV,
,KV,
i K v
.KV,
,KV,
,KV,
,KV,
,KV,
1)
2)
3)
4)
5)
6)
7)
6)
9)
10)
11)
12)
13)
14)
IS)
271
-------
SUBROUTINE SETUP ( KN « KV, NSEGS, ISW )
EXMG = SOIL (KN,KV,l7>
ExNA = SOTL (KN,KV,18)
UCAS04 = SOIL CKN»KV,19»
UMGS04 = SOIL (KN,KV,20)
INTER = SOIL (KN,KV,2l)
GO TO 2000
C<> THIS ENTRY MOVES OF SOLUTION OF BOTTOMMOST SEGMENT INTO BTANK ARRAY
tooor, USING A LlNFAR MIX AFTER CONVERTING UNITS OF COMPUTATIONAL MODE TO
Co<5«c PARTS PER MILLIONfPPM) WHERE PPM = MOLES PER LITFR » EQUIVALENT
coooo WEIGHT OF SPECIFS « VAtENCF OF SPECIES * 1.E3 . ALSO MOVES SOLUTION
C.«
(SOIL (
-i
TO NEXT
.K)» EQVWT(K)
LOWER SEGMENT
VALNCF(K)
= SOIL (KN,KR,19>
= SOIL UN,KR»20!
TO ?000
INTO UPPERMCIST SEGMFNT
/ (EOVWT(K) « VALNCE(K) « 1.E3)
SOIL (KN,KP,i9!
SOIL (KN,KP,20!
GO TO 95
120 IF ( KV. GT. 1 > GO
PLACE APPLIED WATFRS
DO 130 K = 1, 8
L = K + 1
SOlLUNi 1'K) =
130 CONTINUE
GO TO 2000
AT THE COMPLETION OF EQILIRRATIOK OF FACH SEGMENT
<-«««« ARRAY SOIL By TRANSFER OF DATA FROM COMPUTATIONAL
C««<>« SEQUENCE. ALL UNITS ARE IN COMPUTATIONAL UNITS.
ENTRY RELOAD
1) = SLCA
?} - SLMG
3) = SLNA
4) = SLCL
THIS ENTRY RELOAnS
VARIABLE NAME
SOIL (KN,KV,
SOIL (KN.KV,
SOIL (KN,KV,
SOIL (KN,KV,
272
-------
SUBROUTINE SETUP ( KN , Kv, NSEGS, ISW )
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
SOIL
(KN
(KN
(KN
(KN
(KN
(KN
(KN
(KN
(KN
(KN
(KN
»KV<
»KV,
iKVi
,KV,
»KV<
iKV,
,KV,
tKV,
,KV,
,KV.
,KV.
- 6!
i 7)
i 8)
i!2>
,13)
i!6)
,17)
>!«>
,19!
i20)
.21)
=
-
-
-
—
-
=
—
-
-
-
SLHC03
SLC03
SLN03
GYPSUM
SLIME
EXCA
FXMG
EXNA
UCAS04
UMGSo^
INTER
RETURN
END
273
-------
SUBROUTINE FRSTTM ( KNODE , KSEG i Iw )
r<>oco IS USED ONLY ONCE FOR EACH SOIL SEGMFNT AND USES r>ATA FROM INITIAL
roooo SOIL ANALSIS.
COMMON /SOLUTF / SLCA, SLMG, SLNA, SLCL , SLS04, ?LHC03, SLC03,
1 SLN03 ,UCAS04 , UMGS04
COMMON /SOLID/ CATEXi GYPSUM, SLIME ,BDF.NS» EXCA ,EXNA, EXMG *
I EXTRAC , DEPTH, ^SEGVOL , INTFR
DIMENSION CA (10)
DATA (TECA= ROOT OF A THIRD DEGREE POLNOMIAL . LET THF^ODyNA^1 jc EQUILIBRIUM
C.«*°a CONSTANTS FOR UNDISSOCI ^TFD CALCIUM SULFATE AND MAGNESIUM SULFATF
(•oaoo = A.9E-3 AND 5.9E-3 RESPECTIVELY
CHFCK FIRST TO DETERMINE IF SULFATE IS IN SOLUTION - IF NO SULFATE IS
r*o«« SOLUTION BYPASS COMPUTATIONS WHERE SULFATE IS A VARIABLE
IF(SLSOA.LE.O.O) GO TO 60
CALCULATE COEFFICIENTS FOR THIRD DEGREE POLYNOMIAL
CO 50 K = 1,10
CAU) = 0.0
50 CONTINUE
GSQR = GAMMA2 <> GAMMA2
CA!4> = GSQP * PSQR
CA(3) = GSQR *MTECA + TECMG) * GSQR c GSQP « (SLCA + SLMG -SLS04)
CA(2) = TECAo TECMG + GSQR o(TECMG« SLCA + TECA« SLMG - SLS04 «•
1 (TECA + TECMG) )
CA(1> = - SLSOA ° TECA * TECMG
roo«o CALCULATE UNCOMBINED S04 AS ROOT GREATER. THAN ?ERH AND LESS THAN
Coo«o SLS04
CONST = 0.0
CALCULATE A POSITIVE POOT CRF.AL> IN RANGE XMIN - XMAX
XMIN = o.o
X^AX = SLSO<*
rocoo SET CONVERGFNCE CRITERION (TESTA) SO THAT F(X) IS LESS THAN l.E-6
TESTA = l.E-6
MOLE = 6RFRSTIM
Fs04 = ROOT (CONsT,XMAX, XMlN, MCLE , CA , TEsTA )
IF(MOLE.NE.5RTRUE ) GO TO 900
CALCULATE CONCENTRATION OF UNDISSOCIATED CAS04 (MOLS/LITER)
UCAS04 = SLCA « FS04 « GSQR / (TECA -t- GSQR « FsO't)
CALCULATE UNCOMBINED CA (MOLES/ LITER)
SLCA = SLCA - UCAS04
CALCULATE CONCENTRATION OF UNPISSOCIATED MGS04 (MOLS/ LITER)
UMGS04 ~ SLMG « FsO^ « GSQR / (TECMG * GSQR <> FS04)
CALCULATE UNCOMBINED MG (MOLS / LITER)
274
-------
SUBROUTINE FRSTIM ( KNODE . KSEG , Ih )
CHANGE CONCENTRATIONS O^ DIVALENT IONS FROM MOLS PER LITER TO
EQUIVALENTS PER LITER
SLMG = SLMG « 2.
SLCA r SLCA « 2.
GO TO 70
60 FS04 = 0.0
CALCULATE EXCHANGEABLE CALCIUM (EXCA) IN EQUIVALENTS PFP GRAM
(•»««« LET EXCHANGE CONSTANT CA-MG = o.e? AND FOR CA-NA = 2.0
CALCULATE ACTIVITY COEFFICIENT FOR MONOVALENT ION (GAMMai)
70 GAMMAI = GA^MA (1)
EXCA = SLNA « PAMMA1/12-0 « SQRTtSLCA « GAMMA2)) + 1.0 * 6.7E-1
1 « SLMG /SLCA
EXCA = CATE* / EXCA
CALCULATE EXCHANGEABLE SODIUM IEXNA; IN EQUIVALENTS PEP GRAM
EXNA = EXCA « SLNA » GAMMAI / (2 = 0 <» SQRT (SLCA « GAMMA2))
CALCULATE EXCHANGEABLE MAGNESIUM (EXMGJ IN EQUIVALENTS PER GRAM
EXMG = CATEX - FXCA - EXNA
£»««<> CHANGE UNITS FROM EQUIVALENTS PER GRAM AND EQUIVALENTS PER
C« ERROR MESSAGES ARE "RTNTET HERE
900 WRITE (IW,910) ^NODE , KSEG
9lo FORMAT (1H1,1AHERROP IN NODE ,13 ,2X ,1?H5EGMENT NO. , 13,
1 15HBECAUSE OF os»0 )
<56Q WRITE (IW,965) FS04 , SLS04
965 FORMAT (8X i 20HROOT IS OUT OF RANGE i 7HROOT = »F10.5 i
} 8X,15HRANGE = 0.0 TO » F10.5 )
CALL EXIT
?OCO RETURN
END
275
-------
SUBROUTINE PYPEX (KNODE, KSEG, IV,, PRATIO )
r<>o««»o Cr C-YDSUM ( CA50402H20 ) AND THE REQUIREMENTS OF CALCIUM
f.xtoo AM: MAGNESIUM FOR THE UNDT.SSOCIATED ION FAIRS CASH'* AND
/SOLUTE / SLCA, <^LMG, SLNA, SLTL , SLS04, ?LHr03, SLCOj,
I SLN03 ,UCAS04 , UM^SOA
COMMON /SOLTD/ CATEX, GYPSUM, SLIME ,BDENS, EXCA ,EXNA, EXMG ,
i FXTRAC , DEPTH, SEGVOL , JNTFP
CALCULATE ACTIVITY COFFFICIENCT OF DIVALENT ION (GAMMA?!
GAM^A? = GAHMA (2)
GSQR = GAMMA2 « GAMMA2
CALCULATE CHANGES IN CONCFNTR AT I ONS OF CALCIUM ANH SULFATF (X). GIVEN THp
f'.>o^o SOLUBILITY PRODUCT (KSP) OF GYPSUM = 2-4E-5. STLVF FOR (X) USING
CO QUADRATIC rORM OF KSP = ((SLCA + X11SLSC4 +X!)c t-AMMA? <»»2 . WHEN
r is ZERO o°
roooo NEGATIVE SOLUTION IS SATURATED OP SUPf RSAT'JR AT FD . COMDUTF COEFFIClF^
c<>o«o FOR QUADRATIC EQUATION
AO ct ( i > = i.o
CA(2> = SLCA + SLSOA
CA!3) = (SLTA » SLSCA) - (2.AE-5 / GSQR)
KOLF = 6RGYPSUM
CnrvFPT GYPSUM (SnLID) FROM MOLS/GRAM TO HOLES PEP LITEP - PRATIO IS
r<>o»-:> MOISTURE COMTFNT IN PFR^ENCT - MCLS/LITER = (1.0 /PRATlO)<> 1.E5
GYPSOL = GYPSUM « (1.0 /RPATIO) « l.E"5
X = ROOTWO (KNODE,KSEG,CA, MOLE, IW)
IF(X.LT.O.O) GO TO 60
CCNTITICN-SOLUTION IS UNDERSATUPED
Y = GYPSOL - X
IF (Y.LE.0.0) GO TO 50
fOM:ITION - ENOUGH GVSUM IN SOLID STATE TO SATURATF SOLUTION - UPDATE
reo«« CALCIUM AND SULFATE CONCFNTRAT IONS •
SLCA = SLCA + X
SLSOA = SLS04 t- x
GYPSOL = Y
GO TO 70
CONDITION-NOT ENOUGH GYPSUM IN SOLID STATE TO SATURATE SOLUTION DISSOLVE
C«««
-------
SUBROUTINE GYPEX (KNODE, KsEG, IW, PRAflO )
CONDITION - UCAS04 MOT AT MAXIMUM CONCENTRATION
to = 4.987E-3 - UCAS04
IF (W.LT.GYPSOL) 60 TO 80
CONDITION-NOT ENOUGH RESIDUAL GYPSUM TO BRING UC«S04 Tn MAXIMUM CONCENTRAT
UCAS04 = UCAS04 + GYPSOL
GYPSOL = O.n
GO TO 100
CONDITION - ENOUGH RESIDUAL GYPSUM TO MEET MAXIMUM CONCENTRATION OF UCAS04
80 GYPSOL = GYPSOL - W
UCAS04 = 4.°87E-3
CONVFRT-RESIDUAL GYPSUM FOR MOLS/LITER BACK TO MiUS/GRAw
90 GYPSUM = GYPSOL « PRATlO / 1.E5
GO TO 200
CALCULATE CHANGES IN CALCIUM AND SULFATE CONCENTRATIONS (X) TO SATISFY
UNDISSOCIATED CA504 (UCAS04) - WHERE EQUATION OF FQUILIBRIUM IS
(UCASo^ -X ) K = GAMMA?»*?«SLCA + X HSLSo** + X)S 6ND K = THERMODYNAM
EQUILIBRIUM CONSTANT = 4.
UPDATE CONCENTRATIONS OF CAS04, CA AND S04
UCASOA = UCAS04 - X
SLCA = SLCA + X
SLS04 = SLS04 + X
CALCULATE CHANGES IN MAGNESIUM ANH SULFATE CONCENTRATIONS (X! TO SATISFY
(•««•»« UNDISSOCIATFD MGSO-* (UMGS04)- WHERE EQUATION OF EQUILIBRIUM IS
C»o«ff (UMGS04 -X ) K = GAMMA2»°2((SLMG +X JISLSO^ +x)) AND K = THERMODYNAM
EQUILIBRIUM CONSTANT = «5.9E-3 - COEFFICIENTS FOR SECOND DEGREE
EQUATION ARF -
?00 CA(l) = GSQR
CA(2) = GSQ" » (SLMG * SLSQ^) * 5.9E-3
CA(3) = (GSOR » SLMG * SLS04! -(5.9f-3 '•> UMGS04)
MOLE = 6RUMGS04
X = ROOTWO fKNODF, KSEG, CA» MOLE* Ito)
UPDATE CONCENTRATIONS OF Mf,S04, MG AND S04
UMGSOA = UMGS04 - X
SLMG = SLMG * X
SLS04 = SLSIH * x
RETURN
END
277
-------
SUBROUTINE ExCANA ( KNODE , KSEG » Iw PRATIO )
roa«o THIS SUBROUTINE CALCULATES THE SODIU^ (MONOVALFNT) EXCHANGE WITH
ro««o CALCIUM(DIVALENT) . THE G«PON EQUATION ExPRFSSING THF
RELATIONSHIP OF CALCIUM AND SODIUM IN TERMS OF COCFNTRATION5
^AS USED. LFT (X) REPRESENT THE CHANGF IN CONCENTP AT ION OF THE
DIVALENT lQt» IN THE SOLTD pHASE IN UNITS OF MOLS -L"R GRAM.
CALCULATE COEFFICIENTS rOR THE FOURTH DEGREF POLYNOMIAL-- BETA
IS THE CONVERSION FACTO" TO CONVERT MOLS PFR GRAM TO MOLS
f«o«o PER LITER AND HAS UNITS OF GRAMS OF SOLID PER LITFR 0*7 SOLUTION
<-o««» EXCHANGE CCFFFICIENT K (DFX> = 7.07F-1
COMMON /SOLUTE / SLCA» SLMG» SLNA, SLCL • SLSO^, SLHCM* $LC03<
1 SLN03 ,UCAS04 , UMGS04
COMMON /SOLID/ CATEX, GYPSUM* SLIME ,PDENS, EXCA ,EXNA, EXMG ,
1 FXTRAC , DFPTH, SEGVOL , INTFR
DIMENSION CA (10)
DIMENSION TU)
DEX = 7-07E-1
ro« IF THERE IS NO SODIUM IN SOLUTION OR SOLID FH.'.SF. EXIT ROUTINE
IF ( SLNA .FQ. o.O. AND.EXNA . EC. 0.0 ) GO TO 2000
CALCULATE MQNOVALENT ACTIVITY COEFFICIENT GAMMA1
GAMMA1 = GA^MA (1)
DO 50 K = 1, 10
CA (K) = o-O
50 CONTINUE
BETA = ( 1.0 / PRATIO) « ltf-5
ARRANGE ORDER OF COEFFICIENTS IN INCREASING POWERS OF (X!
GSQR = GAMMAl « GAMMA1
BETA2 = BETA « BETA
DEX2 = DEX « DEX
CA(1) = ( GSQP # EXNA » EXNA « SLCA ) - ( DFX?
1
« SLNA « SLNA )
CA(2) = GSQ" « EXNA « ( 4.0 » SLCA
» FXCA o EXCA
1
EXC«
+ BETA <> EXNA )
EXCA >
2.0 « DEX?«
CA(3)
I
CA(4)
» DEX2« BETA «
« SLNA« SLNA
GSQR ••
» i SLNA * 2.0 » BETA
= 4.0 * GSQR « ( SLCA * BETA « FXNA ) -
EXCA « ( BFTA « EXCA * 2.0 * SLNA )
= 4.0 « BETA <> ( 2.0 « DEx2 « ExCA * BETA
1 SLNA )
CA(5) = - 4.0«BFTA2 « DFX2
ISW = 1
CALCULATE A POSITIVE ROOT (REAL) IN RANGE XMIN - XMAx
XMIN = 0.0
XKAX = EXCA
IF ( EXCA .GT. FXNA
60 MOLE = 6PEXCANA
C«o«*> SET CONVERGENCE CRITERION (TESTA) SO THAT F(x) IS LESS THAN l.E-R
TESTA = l.E-8
X = ROOT ( CONST, XMAX , XMIN , MOLE
IF ( MOLF . EQ . 5RTRUE 1 GO TO 66
GO TO ( 75, 68 ) , ISW
Coo«o TEST ROOT TO INSURE NO UPDATED CONCENTRATIONS ARE LESS THAN ZERO
) XMAX = EXNA
CA , T^STA )
278
-------
SUBROUTINE EXCANA ( KNODE , KSEG i IW, PRATIO )
T(2) = EXCA - X
TO) = SLNA - t 2.0 « BET* « X )
T(4) = SLCA •»• ! BETA « X )
DO 70 K = 1,4
IF ( T (K) . LE.0.0 ) GO TO 75
70 CONTINUE
GO TO 80
75 IF ( ISW . EQ. 2 ) GO TO 1000
CALCULATE A NEGATIVF »OOT (REAL) IN RANGE XMIN - XMAx
XMIN = - XMAX
XMAx= 0-0
ISW = 2
GO TO 60
c«« UPDATE CALCIUM AND SODIUM CONCENTRATIONS IN SOLUTION AND SOLin
C« PHASES
80 EXNA = T(l)
EXCA = T(?)
SLNA = T(3)
SLCA = T(4>
GO TO 2000
inno WRITE ( iw
(
(
1(105 FORMAT
1
WRITE
IClO FORMAT ( // 8X,
1
CALL EXIT
?nr.O RETURN
END
1005 ) KNODE , KSEG
, 4?HNO REAL ROOT EXISTS IN KANr,£ FOR CA-NA FOR
1H1
7HNODE = , 13 , 5X, 10HSEGMENT = ,
IW, 1010 ) X, XMIN, XMAX, MClE
4HX = . F20.8 , 3X, 7HX"IN = ,
7HXMAX = , F20.6, ?X, 7HMQLE =
, 3X,
> )
279
-------
SUBROUTINE E.XCAMG t KNODE, KSEG, Iw, PRAiIO )
COMMON /SOLUTE / SLCA, SLMo, SLNA, SLCL 1 SLS04, 5LHCC3. SLC03,
I SLN03 ,UCAS04 , UMGSO*,
COMMON /SOLTD/ CATEX, GYPSUM, SLIME ,PDENS, EXCA ,EXNA, EXMG ,
i EXTRAC , DEPTH, SEGVOL , INTER
C««»« THIS SUBROUTINE CALCULATES THE MAGNESIUM (DIVALENT) EXCHANGE WITH
(•oooo CALCIUM (DIVALENT) . THE COMMON EQUATION EXPRESSING TH£ RELATiONSHtp
<-<>o<>ft OF CALCIUM AND MAGNESIUM IN TERMS Of CONCENTRATIONS WAS USED.
<-<>«<><> LET (x) REPRESENT THE CHANGE IN CONCENTRATION OF THE DIVALENT IPN
r<:-o<><> IN THE SoLIr PHASE IN UNITS QF MOLS PER GRAM. CALCULATE THE
(-0000 COEFFICIENTS OF SECOND DEGREE POLYNOMIAL = *ETA = CONVERSION
r« FACTOR TO CONVERT MOLS PER GRAM TO H.cLs PER LITER AND HAS UNITS OF
r»«<><>« K (DEx) = 6.7E-1
DIMENSION CA (10)
DEX = 6.7E-1
BETA = ( l.n / PRATIO ) <> 1.E5
CA (1) = BET* « ( ).0 - DEx)
CA(2)=BETA *> ( T1EX <> EXCA t- EXMG ) + HEX » SL^G + SLCA
CA(3) = SLCA « FXMG - ( DEX « SLMG « EXCA )
MOLE = 6RCALMAG
X = ROOTwO t KNODE , KSEG , CA , MOLE , Iw )
c«« UPDATE CONCENTRATIONS or CALCIUM ANP MAGNESIUM IN SOLUTION AND
SOLID PHASES
SLCA = SLCA + ( BETA « X )
SLMG = SLMG - ( BETA <> x '
EXCA = EXCA - X
EXMG = EX"G + X
RETURN
END
280
-------
SUBROUTINE CALCAR ( KNODE, K$EG, IW, PRATln )
C<>« WITH THE CARBONATES. IF CACOj ( LIMt ) is PRESENT IN THE SOLID
C<>»<>« PHASE THE INDICATOR SLIME WILL RE SET TO OK'F , IF LlMF IS NOT
C»« PRESENT SLIME WILL BE SET TO ZERO . LOG10 (KPRIME )= -1.68 »
C»o LOG10 (PRATTO) - LOGIC (INTER) , WHFRF pPATIc = MOISTURE CONTENT
IN PERCENT, KPRIMF = SLCA <> SLHCOS tt SLHCOS « GAMMAI « GAMMAI
« GAMMA2 ANT1 INTER = INTERCEPT, THE EXPRESSION FOP KPPIME WAS
DERIVED USING A SIMPLE LINEAR REGRESSION ( BY DlJTT ET AL )
C<>« WITH A LOGARITHMIC TRANSFORM.
COMMON /SOLUTE / SLCA, SLMG, SLN'A, SLCL , SLS04, SLHCC3, SLC03,
1 SLN03 ,UCASQ4 , UMGS04
COMMON /SOLTD/ CATEX, GYPSUM, SLIME .P.DF.NS, EXCA ,EXNA, EXMG ,
1 FXTRAC , DEPTH, SEGVOL , INTFR
DIMENSION CA (in)
REAL KPRIME , INTER
CHFCK IF LIME IS PRESENT IN SOLID PHASE
IF ( SLIME .EQ. 0.0 ) RO TO 50
C»oos LIME IS PRESENT CALCULATE KPRIME
IF (INTER • NE • 0.0 ) GO TO 43
40 GAMMA] = GAMMA d)
GAMMA? = GAMMA (2)
KPRIME = SLCA « SLHC03 «SLHCo3 » GAMMAl » GAMWA1 * GAMMA 2
COMPUTE INTERCEPT ONCF ONLY pOR E^CH SOIL SEGMENT
INTER = KPRTME <> 1.66)
43 KPRIME = INTER / IPRATIO »» 1.6« )
DO 45 K = 1,10
CA (K) = 0.0
45 CONTINUE
GO TO 55
coo»« LIME IS NOT PRESENT IN SOLID PHASE ASSUME TNTER = 4.52E2
5o INTER = 4.5?E2
KPRIME = INTER / ! PRATTO «» 1.68)
GO TO 100
CALCULATE SOLUBILITY OF LIME ! CAC03) 9 USE EXPRESSION K = SLcA
*SLHC03« SLHC03 »» 2 » GAMMAl o<>2 / SLH2C03 ( THEN LET KPRIME
= K
-------
SUBROUTINE CALCfiR ( KNODE, KSEG, !w,
5RTRUE
TESTA = l.E-6
X = ROOT ( CONST
GC TO ( 67, 6fl
67 IF ( MOLE .FQ.
CALCULATE A NEGATTyF
XMAX = 0.0
XMIN = - XM-AX
ISW = 2
GO TO 60
68 IF ( MOLE .NE.
C»<> UPDATE CONCENTRATIONS
C.oooo IN SOLUTION
69 SLCA = SLCA + X
SLHC03 = SLHC03 + X
70 IF (SLIME . GT. 0.0 )
100 TEMP = SLCA » SLHC03
IF ( TEMP .r-E. KrRU't
GO TO 2000
incc V,RITE ( iw,
1005 FORMATI1H1,
1
CALL EXIT
?niQ RETURN
END
XKAx
) , ISW
SRTRUE ) r-C TO f.9
ROOT (REAL) IN RANGE
XMIN , MOLE i CA , TESTA i
XMIN - XMAX
) GO TO 1000
OF CALCIUM (SLCA)
GO TO 2000
> 5LUC03 <> GAHMA1
! GO 70 43
AND BICARBONATE (SLHC03)
G A M M A 1
1005 ) KNODE ,
-------
SUBROUTINE SCRIBE (IW)
COMMON /SOLUTE / SLCA, sLMGi SLNA, sLCL 9
1 SLN03 ,UCAS04 , UMGS04
COMMON /SOLID/ CATEX, GYPSUM, SLIME ,RDENS
i EXTRAC , DFPTH, SEGVOL , INTER
REAL INTER
WRITE ( IW
SLS04, t;LMCC3, SLC03,
. EXN'A, EX KG ,
10 )
1
SLCA, SLM&, SLNA, SLCL i SL504 <,
SLCCn, EXTRAC, SEGVOL, CATFX, GYBSliM, ?:
EXMG, UtASOt
SLHC03
SLND3.
10 FORMAT
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
RETURN
END
6X,
ex,
ex,
ex,
ex,
ex,
ex,
ex,
px,
ex,
ex,
ex,
ex,
ex,
ex,
fiX,
ex,
ex,
ex,
ex,
8X,
DFPTH , E
INTER
9HSLCA
9 H S L M <••
9HSLNA
9HSLCL
9HSLSH4
9HSLHC03
9HSLN03
9HSLC03
9HEXTRAC
9HSFGVOL
9HCATEX
9HGYPSUM
9HSL IMF
9HBDEMS
9HDEPTH
9HEXCA
9HEXNA
9HEXMG
9HUCAS04
9HUMGSOA
9HINTER
.XCA,
— r
S
r
r
r
=
~
=
r
r
»
"Si
—
s
=
=
SI
r
-
—
r
EXN
F20
F?0
F2U
F20
F20
F20
F20
F20
F20
F20
F20
-F20
F?0
F?0
F20
F20
F20
F?0
F2o
F20
F20
A,
.f:/
, £,/
.e/
.e/
.e/
.e/
.e/
.e/
.e/
.a/
.e/
.e/
.a/
.e/
.e/
.a/
.e/
»6/
.S/
.e/
.e/
283
-------
FUNCTION GAMMA ( KVAL )
C ooooe
roooo THIS FUNCTION CALCULATES THE ACTIVITY COEFFICIENT (GAMMA)
foooo FOR MONOvALFNT ION ( K.VAL = 1) AND DIVALENT ION ( KVAL = 2 )
C eoo-jo
COMMON /SOLHTF / SLCA* SLMG, SLNAi SLCL . SLS04, UHCC3, SLC03*
1 SLN03 ,UCAS04 , UMGSOi*
CALCULATE THE IONIC STRENGTH OF ALL SPECIES IN SOLUTION (U) WHERE
C«o«o u= 1/2 (SUMMATION 0^ CONCENTRATION OF IQN ( MOLS/LlTFR) TIMES
C«««o VALENCE OF ION SQUARED - FOR ALL SPECIES)
U = 2.0 « (SLCA + SLMG + SLS04 f SLCC3) + 5-E-l » (SLNA + SLCL i
1 SLN03 * SLHC03 )
CALCULATE ACTIVITY COEFFICIENT «»«o«*»«»»««>oo«« ION USING DEBYE-
C««<>« HUCKEL FORMULA AS A FUNCTION OF lONjC STRENGTH flNn VALENCE
COMP = 0-0
COMp = SORT (U)
GC TO ( 5o, 100) , KVAL
50 GAMMA = EXP ( -1.17201 <* CO^P / ( 1.0 + COMP )!
GO TO 200
100 GAMMA r EXP ( -4.66^0^ <> COMP / ( 1.0 + COMP ))
?00 RETURN
END
284
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
REPORT NO.
EPA-600/2-77-179C
3. RECIPIENT'S ACCESSIOf*NO.
4, TITLE AND SUBTITLE
PREDICTION OF MINERAL QUALITY OF IRRIGATION RETURN
FLOW, VOLUME III, Simulation Model of Conjunctive Use
and Water Quality for a River System or Basin
5. REPORT DATE
August 1977 issuing date
6. PERFORMING ORGANIZATION CODE
. AUTHOR(S)
8. PERFORMING ORGANIZATION REPORT NO
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Bureau of Reclamation
Engineering and Research Center
Denver, Colorado 80225
10. PROGRAM ELEMENT NO.
1HB617
11. CONTRACT/GRANT NO.
EPA-IAG-D4-0371
12. SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Lab.-Ada, OK
Office of Research and Development
U.S. Environmental Protection Agency
Ada, Oklahoma 74820
13. TYPE OF REPORT AND PERIOD COVERED
Final
14. SPONSORING AGENCY CODE
EPA/600/15
15. SUPPLEMENTARY NOTES
VOLUMES I, II, IV, V
(EPA-600/2-77-179a, 179b, 179d, 179e)
16. ABSTRACT
The development and evaluation of modeling capability to simulate and predict the
effects of irrigation on the quality of return flows are documented in the five
volumes of this report. The report contains two different modeling packages which
represent different levels of detail and sophistication. Volumes I, II, and IV
pertain to the model package given in Volume III. Volume V contains the more
sophisticated model. User's manuals are included in Volumes III and V.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS
COS AT I Field/Group
Mathematical Model, digital simulation,
scheduling, Irrigated land, Evapotrans-
piration, Agriculture, Agronomy, water
pollution, water loss
Irrigation Return Flow
02 C/D
13. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (ThisReport)
Unclassified
?1. NO. OF PAGES
295
20. SECURITY CLASS (Thispage)
Unclassified
22. PRICE
EPA Form 2220-1 (9-73)
285
U.S. GOVERNMENT PRINTING OFFICE: 1977-757-056/6551! Region No. 5-1 1
------- |