Nothing Special   »   [go: up one dir, main page]

Rec Erc 79 07

Download as pdf or txt
Download as pdf or txt
You are on page 1of 92

REC-ERC-79-7

Engineering and Research Center


Bureau of Reclamation

.
July 1979
MS-230 (l-76)
Bureau of Reclamation
TECHNICAL RI EPORT STANDARD TITLE PAG
1. REPORT NO. 3. RECIPIENT’S CATALOG NO.
REC-ERC-79-7
4. TfTLE AND SUBTITLE 5. REPORT DATE

Gate Stroking July 1979


6. PERFORMING ORGANIZATION CODE

7. AUTHOR(S) 8. PERFORMING ORGANIZATION


REPORT NO.
H. T. Falvey
REC-ERC-79-7
P. C. Luning
9. PERFORMING ORGANIZATION NAME AND ADDRESS 10. WORK UNIT NO.

Engineering and Research Center


Bureau of Reclamation 11. CONTRACT OR GRANT NO.

Denver, Colorado 80225


13. TYPE OF REPORT AND PERIOD
COVERED
12. SPONSORING AGENCY NAME AND ADDRESS

Same

14. SPONSORING AGENCY CODE

15. SUPPLEMENTARY NOTES

16. ABSTRACT

Gate stroking consists of a series of continuous or discontinuous gate motions which produce a
desired water surface profile in a canal. The first obvious application of gate stroking within the
Bureau was on the Granite Reef Aqueduct, Central Arizona Project. The mathematical
development of the technique is outlined and methods for treating the unique conditions found
on Bureau aqueducts are presented. Sufficient computer documentation is provided to permit
application of the program to other aqueducts. The technique can be applied either to the entire
aqueduct or to the component reaches. If applied to the aqueduct, the output consists of both
gate and pump schedules. When applied to reaches only, the set of gate schedules is produced.

7. KEY WORDS AND DOCUMENT ANALYSIS

3. DESCRIPTORS--/ *gate control/ *open channels/ *computer models/ pumps/ water surface
xofiles/ aqueducts/ scheduling/ mathematical analysis

1. lDENTfF IERS-- / Central Arizona Project/ Granite Reef Aqueduct, Arizona/ Salt-Gila
\queducts/ *gate stroking/ method of characteristics
:. COSATI Field/Group 13 COWRR: 1313
8. DlSTRIBUTlON STATEMENT 19. SECURITY CLASS 21. NO. OF PAGE
ITHIS REPORT)
\voilobfe from the Notional Technical Information Service, Operations
UNCLASSIFIED
86
division. Sfxingfield, Vlrginio 22151.
20. SECURITY CLASS 22. PRICE
REC-ERC-79-7

GATE STROKING

by
H. T. Falvey
P. C. Luning

duly 1979

Hydraulics Branch
Division of Research and
Division of Data Processing
Engineering and Research Center
Denver, Colorado Sl METRIC

UNITED STATES DEPARTMENT OF THE INTERIOR * BUREAU OF RECLAMATION


ACKNOWLEDGEMENTS

Funds for this study were made available from the Central
Arizona Project through coordination efforts by the E&R Center
Automation Team. The careful review of this report and the
continuous support by engineers in the Operations Division,
Central Arizona Project is gratefully acknowledged.

The information contained in this report


regardingcommercialproductsorfirms may not
be used for advertising or promotional purposes
and is not to be construed as an endorsement of
any product or firm by the Bureau of
Reclamation.
CONTENTS

Page

introduction ....................... 1
Summary and conclusions ................... 1
Analytical development ................... 2

Basic equations ...................... 2


Numerical method .................... 5

Boundary conditions .................... 5

Reach versus aqueduct computation ..............


Changes in prism cross section ................ i
Turnouts ........................ 7
Gates ......................... 7
Siphons ........................ 8
Specifications constraints .................. 9

Method of utilization .................... 10


Validity of the program ............. : ..... 11
Program description ..................... 11

Background ....................... 11
Problem formulation .................... 11
Initial conditions ..................... 12
Boundary conditions .................... 12

Depth schedules .................... 12


Discharge schedules ................... 12
Turnout schedules .................... 12
Pumps ........................ 12
Gate positions ..................... 13
Inverted siphons .................... 13
Boundaries where’no control structure is present ......... 13

Physical descriptors .................... 13


Modular construction ................... 13
output ........................ 13
Development and structure of the program ............ 13
Program details ...................... 14

Bibliography ....................... 15

APPENDIX

Gate stroking computer program listing . . . . . . . . . . . . . . 17

i
CONTENTS-Continued

FIGURES
.

Figure Page
1 Characteristic lines .................. 2
Computational regions ................
i Region 1 ..................... :
4 Region 2 ...................... 3
Region 3 ..................... 3
: Interpolation scheme ................. 4
7 Region 4 ..................... 4
Region 5 ..................... 4
ii Dead-band operation with a pump ............
10 Turnouts and change of cross-sectional area ......... :
11 Computational method across pool boundaries ........ 7
12 Gate definition sketch ................ 8
13 Coefficient of free and submerged Tainter-gate discharge
forr/a= 1.5 ................... 8
14 Discharge coefficient as a function of submergence ...... 8
15 Dead-band operation of gates .............. 8
16 Typical siphon layout ................ 9
17 l!oss coefficients .................. 10

ii
INTRODUCTION In addition, the aqueduct is designed without any
wasteways or rere9ulating reservoirs. When put
The change from one discharge to another in into operation, the canal is to be operated at
open-channel flow always creates a disturbance nearly the design capacity. Thus, there is very
in the water-surface elevation. The magnitude of little margin for errors and the control of tran-
the disturbance is related to the manner in which sients is a significant factor. These considera-
the change is accomplished. In some cases, the tions’coupled with a desire to miiiimize onpeak
disturbances persist for long periods with large pumping costs led to the development of thegate
amplitudes. In 1969, Wylie [l]” developed a stroking technique within the Bureau.
method to control the disturbances within certain
limits. The objective of the method was to Wylie’s original concept was followed very
produce a predetermined variation in the water- closely in the development which follows;
surface elevation at one location in a canal pool however, provision had to be made to include
by properly varying the boundary conditions at structures like turnouts, siphons, and free
each end of the pool. flowing tunnels in the Bureau’s computer
program. In addition, a technique had to be
Since pool boundaries are usually gates, the developed in which the initial conditions were
technique of water-surface control is called not always at steady state.
“gate stroking.” The term comes from a similar
procedure used in closed conduits known as The purpose of this report is to outline the
“valve stroking.” Although the term “stroking” mathematical development; present the methods
as used in this context cannot be found in the used in handling the unique conditions found on
dictionary, its definition can be implied from the Bureau aqueducts; illustrate how the gate
several meanings commonly listed. Thus, strok- stroking is utilized; and provide sufficient
ing is “any of a series of continuous or documentation for the computer program, ap-
discontinuous efforts to do, produce, or ac- pendix A, so that, it can be applied to other
complish something, especially a successful aqueducts.
result.” Based upon its usage and the general
meaning of stroking, gate stroking can be defined
as, “‘a continuous or series of discontinuousgate SUMMARY AND CONCLUSIONS
motions which produce a predetermined water-
surface variation in a canal. ‘. Gate stroking is a series of continuous or
discontinuous gate motions which produce a
Since the original paper by Wylie, the number of desired water-surface profile in a canal. Ttie first
gate stroking applications has been minimal. obvious application of gate stroking within the
O’Loughlin [2] and Gientke [3] are two of the few Bureau was on the Granite Reef Aqueduct,
who have referred to application of gate stroking Central Arizona Project.
for specific installations. Perhaps one of the
reasons-there are not more applications is that The mathematical development of the technique
gate stroking requires a scheduled type of is outlined and methods for treating the unique
operation with centralized control. The number of conditions found on Bureau aqueducts are
systems meeting this requirement is rather presented. Sufficient computer documentation is
limited. provided to permit .application of the program to
other aqueducts.
The first obvious application of gate stroking
within the Bureau was on the Granite Reef The technique can be applied either to the entire
Aqueduct, Central Arizona Project. This aqueduct aqueduct or to the component reaches. If applied
consists of several canal reaches separated by to the aqueduct, the output consists of both gate
pumping stations. The aqueduct is intended to be and pump schedules. When applied to reaches
controlled by a computer-assisted remote- only the set of gate schedules is produced.
control system. It anticipated that all delivery
schedules can be reasonably estimated several The most nebulous relationship in the technique
days in advance of the actual need for the water. is the gate discharge equation. Additional
research to better define the gate discharge
l Numbers in brackets refer to literature cited in the coefficients is required.
bibliography.
ANALYTICAL DEVELOPMENT
Basic Equations

Any computational scheme which calculates the


unsteady water-surface profilescan be used with
the concept of gate stroking. The method chosen
by Wylie is the method of characteristics. Two
computational schemes are used with the
method of characteristics. These are (a) the grid
of characteristics and (b) the method of specified
time increments. The grid of characteristics
method was chosen for use in the computer
program. This method is more accurate because
interpolation at the interior of the grid is avoided. x-
Since the computations for each segment are
performed independently of the other segments, Figure 1 .-Characteristic lines.
the grid does not have to intersect specified
points on the boundary. Thus, the usual problem
of joining computations at a boundary, when Using the notation of Wylie, the intersection
using the grid of characteristics, does notarise in point P presents a location at which the solution
the gate stroking application. of the variables x, y, Q and t is theoretically
possible. If conditions are known simultaneously
The equations of motion and continuity in an at either of the R and at either of the S points,
open channel can be expressed asfour particular then a numerical integration of the four char-
total differential equations. These are: acteristic equations will result in a solution
of the variables at P.

J-G
A dt
+g(S-So)=0 (1) Several computational schemes are used in the
computer program. Each of these is dependent
I ct upon which combination of points R and S is
known. The schemes can be related to specific
-z-t
dx
dr
Q
A 1 (2) regions in the x-t plane, figure 2.

gTdv
--t
-d--A dt
dx Q
dt=s;i--

where g =.acceleration of gravity


= top width of water prism
A7 = cross-sectional area of water prism
= flow depth
: = time
= discharge at section
so = friction slope
so = bottom slope
X = horizontal distance
The first two equations are valid along the
positive characteristic 0, figure 1. The last two
are valid along the negative characteristic. Figure Z.-Computational regions.
Each of the computational regions requires FLOW
specific
borders.
information from at least one of its
Each computational scheme then
I INITIAL CONDITIONS -

conveys information to the other borders. In


general, borders occur along any time-distance
line which borders a region, such as A-B,
B-C, C-A, figure 3. Certain borders, however,
have special significance. These occur at specific
locations where depth variations are known or
desired, and are known as boundaries. The canal
segment between boundaries is called a “pool.”
The pool boundaries may or may not coincide
with hydraulic barriers, such as gates. The canal
segment between boundaries at which dis- Figure X-Region 1.
charges are specified, is known as a “reach.”
Several pools may be contained within a reach. FLOW
C
BOUNDARY CONDITIONS
The sequence of the computations for each of the Either depth or discharge
regions which comprise a pool is described in the (not both) specified for
lime between t=O at A
following paragraphs.

Region 1, figure 3.-The initial conditions, that is,


depths and velocities along the entire length of
the pool, must be defined at time equal zero. The
variables at all points within the region can be
determined by successively extending the grid
forward in time. This region is also known 8s the
domain of dependence since conditions within
the region 8re uniquely determined by the
conditions specified on one boundary. When the Figure 4.-Region 2.
computations are completed, the value of the
four variables will have been determined along D
the border A-B-C.
C
Region 2, figure 4.-In this region, conditionsare
known along the border A-B. At the point P,
the value of the Variable x is also known.
s
However, to solve for the conditions at P using
the positive characteristic, one additional vari- I
able must be specified. If the downstream : R’
boundary is 8 pump or 8 delivery point, the F
8
discharge will be specified. As will be seen later, A
if the boundary is 8 gate, the discharge will also
be specified to maintain continuity. Therefore,
along the boundary A-C, .the discharge must be I I
Upstream Downstream
specified, whereas, time and depth will be
determined. Interior points P, are computed from
known values of R, and St as was done in Figure I.-Region 3.
region 1. Upon completion of the computations,
the values of all variables are known along the
border B-C-A. necessary to specify all four of the variables. To
be compatible with the solution from region 2,
Region 3, figure 5.-This region is also 8 domain the values of the four variablesat point A muat be
of dependence since the point P is determined indentical with those 81 point C in region 2. The
from conditions given at the downstream other boundary values along the line A-D are
boundary. Thus, along the boundary A-D, it is determined from the prescribed discharge and

3
water level schedules. When the computations computational procedure has been provided,
are finished, the conditions are known along the figure 8.
borders A-B and C-D. To determine the
conditions along B-C an interpolation is neces- Region 5, figure B.-In this region, the discharge
sary. The interpolation scheme used involves a is specified along boundaries A-D and B-C. The
linear interpolation between points on the computations proceed in a fashion similar to
positive characteristics which cross the boundary those in region 3. However, at the downstream
(A to B, fig. 6). This is different than the scheme boundary B-C, distance and discharge as a
used by Wylie [l] who interpolated between the function of time are known. This computational
points of the diamond shape formed by the grid procedure is required only when the upstream
which crossed the boundary (C to 0, fig. 6). boundary discharge is scheduled.
Neither method seems to have any particular
advantage over the other.

Computational

Figure 7.-Region 4.
L Poromctarr such OS time, depth, and
discharge at thi! point determined by
l$ey mergPolotlon between v0lusr

Figure &-Interpolation scheme.

The computations are extended in time, as


required, (region 3B), by the same computational
scheme used in region 3, figure 2. Region 38 has
no theoretical significance. It is merely a device
for extending the computations forward to an
arbitrary time t without requiring additional
computer memory when t is large.

Region 4, figure 7.-According to Wylie [l), this


region is known as the domain of influence since
it is influenced by the values of the conditions at
point A. The computations proceed from known
values on the borders A-B and C-A, when the
computations have been completed, conditions
along the boundary C-B will have been
determined.

If the discharge at the upstream end of the most


upstream pool is specified, then the depth in that
pool cannot be specified. In this case, a fifth Figure E.-Region 5.
Numerical Method In these equations

The four characteristic equations were inte-


grated using the trapezoidal rule. The trapezoidal
rule is given by

YdX=& (Y. +2Y, +2Y2 + BOUNDARY CONDITIONS


/
2
X0 Reach Versus Aqueduct Computation

. . . 2Y”-, + Y/l) (5) The Granite Reef Aqueduct consists of several


reaches separated by pumping stations. Each
whereh=X= X, -X0 reach consists of several pools which are
separated by control gates. Two methods of
applying the gate stroking technique are
Ifn= l.then possible. For one, the water level variations over
an operating cycle within each reach are
Xl specified. In addition, the operation at each
pumping station is specified. The other technique
YdX=$ (Y. + Y,) (6) involves specifying water level changes at each
s pool in the entire aqueduct. Pumps are treated
X0
like gates having dead bands, figure 9. This type
of operation could properly be called “pump
Applying this rule to the characteristicequations stroking.”
gives
The difference between reach and aqueduct
2-l computations is .mainly one of operating phil-
+ J- (Yp- Y,, + (VP - V,) + osophy. If a pump and water level schedule are
2 ( CR cp >
both specified, then the water level in at least one
pool in a reach cannot be specified. That is, the
f (SR + sp - 24-J (‘p - ?R)= 0 (7‘1

(Xp- Xf?)= vp+vR


2 +CR; ‘I’ >(tp- rR) (8)

-2-l
2
‘+$ (Y/y Ys,+(V,- Vs)+
( cs P >

4 (S, + s, - 2S(g up - 5) = O (9)


L
I

VP+ v, -- cp+ c. (‘p- Cs) (10) TIME -

(XP- X,) =
t 2 2 Figure 9.-Dead-band operation with a pump.
level in one pool must be allowed to seek its own The exact solution of this case requires two
level or “float.” On the other hand, if only the additional equations for the characteristic which
water level variations and delivery schedules are crosses the change in section. These are (1) an
specified, then the necessary pump schedules to equation which defines thedistance between the
achieve this can be uniquely determined. For the known point (A or B) and the change in section,
first case, the pump schedules are of primary and (2) an energy equation across the change in
importance and water level fluctuations in at section. This requires that every step in the
least one pool must be tolerated. In the second computational matrix of pools with changes in
case, water level fluctuations are controlled section be checked to determine whether a
everywhere and the pumping schedule is the crossing has occurred. If one has occurred, then a
flexible component. set of six simultaneous equations must be solved
using an iterative scheme. The iterations are
For a floating water level, the computations pass necessary since the area, top width, and
from a predictive to an analysis stage. The hydraulic radius at the unknown point C are
discharges at each end of the pool are specified functions of the unknown depth at that point.
and the water levels at the boundaries are
calculated, figure 8. This computer program fails If the cross-sectional changes are not too abrupt,
when there is a very small flow into the pool. As the differences in the slopes of the characteristic
the water level drops, the wave celerity at the change can be ignored. This approximation
approaches zero. Since the reciprocal of the requires the solution of five simultaneous
celerity is used in the method of characteristics, equations. Studies h8Ve not been performed to
the range of permissible values for the variable determine what constitutes a large change in
may be exceeded in the computer with small flow section.
depths.
As an alternate approach, the pool can be divided
Changes in Prism Cross Section into subpools. The computations in each subpool
proceed as if each were a pool. The computations
In the derivation of the characteristics method, it begin at the downstream subpool in which the
was assumed that the section was prismoidal in depth must be assumed along the region 3
each pool. This assumption does not hold true if a boundary, (fig. 108). As a consequence, both
change in cross-sectional shape occurs in the the depth and discharge are determined just
pool. For large changes in shape, the char- downstream of the change in cross section. If
acteristic line which crosses the change can continuity of flow across the change is main-
deviate greatly from a straight line (point A to tained, then application of the energy equation
point C, fig. 10A). wiri yield the depths on the upstream side of the

x- 0
*.col SaCI

A. Actual Characteristic Grid 8. Simplified Computational Scheme

Figure IO.-Turnouts and change of cross-sectionel area.


transition. These depths can be compared with the Turnouts
depths resulting from the upstream subpool
computations of region 2. If the differences The concepts developed for changes in cross
between the two are within acceptable limits, the section also apply to turnouts. In this case, the
depth variation assumed at the downstream pool continuity of flow also includes the proper
boundary is acceptable. If the two are outside of accounting for the turnout discharge. Theory and
acceptable limits, the downstream pool boundary experiment indicate that the turnout flow does
depths must be respecified. Generally, rapid not affect the energy head [41. Therefore, energy
changes in water depths do not occur; therefore, is assumed constant across the turnout.
iteration has not proven to be necessary.
Gates
The largest errors with this procedure develop at
the last computation in region 2 of the upstream The general form of the equation for flow under a
subpool, (fig. 1OB). The effect of errors along the free flowing radial (Tainter) gate can be obtained
timeline O-A can be minimized by a relatively from Bernoulli’s equation [5]. It is expressed as:
simple procedure. Using continuity of flow at the
intermediate boundary, the flow depths are
determined in the upstream subpool region 2. 0 = C,bB 4 2g (A.y + y/Zg) (12)
Then, using the energy equation across the
intermediate boundary with known values of where C, = Contraction coefficient
flow and depth from region 4 (and in some cases, b = gate opening
region 3) a slightly different upstream flow depth = gate width
is determined at the intermediate boundary in 6 = difference between upstream and
region 2. This depth based on energy con- downstream water depths
siderations is substituted for that based on v, = upstream mean velocity
continuity considerations. A linear interpolation 9 = acceleration of gravity
is required at the transition, (fig. 11). since the
grid points from the downstream subpool regions If the downstream depth exceeds the gate open-
3 and 4 do not match with those from the ing, the flow is said to be submerged. For this
upstream subpool region 2. case, the equation must be modified to include
losses that occur on the downstream side of the
With the present program it has been assumed gate. It has been argued that the parameters Cc
that the flow passes from one cross section to.the and Vf:/2g are governed by several linear terms
other with no loss of energy. However, the energy which are characterized by the flow pattern [6].
equation could include a transition loss if desired. Thus, it is possible to express the discharge as

0 = C,bi3 (13)
i- &JY,
where Cd = Discharge coefficient
Yl = Upstream depth
The discharge coefficient is a function of the
submergence y/b, trunnion height a/b, and
upstream depth y/b ratios, (figs. 12 and 13).

The Corps of Engineers (7) uses the following


form of the equation for submerged flow

Q = CSB Y’v ’ 2g (A Y + y/29) (14)

UPSTREAM POOL ’ DOWNSTNEAY POOL


where C, = Submerged discharge coefficient.
Porametars such os time. depth, and daschorqe at this
point
at A ond
determined
8. The
by lineor
constant
intqrpolotion
energy qquotion
bctwqcn
is opphqd
vqbms With this equation, the discharge coefficient
at turnouts and chon es in cross-sectional orro. Tha varies linearly on a log-log scale from a value of
qote cquotion is opp Pmd ot control structurqs.
0.04 at y/b = 20 to 0.35 at y/b = 2.5. The
Figure 11 .-Computational method across pool boundaries. effect of the radius of curvature of the sector
II -

-1
1.0 -

,-

I
////////////////////// .I -

Figure 1LGate definition sketch.


.I -

0.6
l-

Figure 14.-Discharge coefficient as a function of


cd submergence.

a dead-band type of operation is preferred. With


this, a gate movement is initiated when the gate
stroking predicts excursions which exceed the
dead band. If this occurs, the time at which the
0 0.1 42 0.3 0.4 QI ob 0.7 Q6 0.9 I.0 1.1 gate motion should have started is calculated.
Yl This time is a function of the gate speed, (fig. 15).
6 Future gate structure installations may include
multi- and variable-speed gate motors. With
Figure 13.-Coefficient of free and submerged radial gate these features the gates could follow the
discharge for r/a = 1.5 openings specified by the gate stroking tech-
nique more closely.
which forms the gate and the trunnion height has
not been systematically investigated for this form Siphons
of the equation.
Normally, siphons are located at the upstream
The equation used in the computer program is end of a pool with a control gate located on the
based upon the equation from the Corps of upstream end of the siphon, (fig. 16).
Engineers, using the gate opening instead of the
downstream depth.

Q = CBb j/2g(Ay + v;2/2g) (15)

where C = discharge coefficient


= Csy2/b
Since the spread in the data is so large (fig. 14)
detailed model studies of the specific structures
are recommended. With detailed studies, the
simpler form of the discharge relationship given
by equation,(l3) should be used.
TINE-
In current practice, it is impractical for a gate to
follow the timewise change in openings that are
specified by the gate stroking technique. Instead. Figure 15..Dead-band operation of gates.

8
Figure 1 B.-Typical siphon layout.

The passage of a pressure wave through the The subscripts 1 and 2 refer to locations in the
siphon actually requires a finite length of time. canal immediately upstream and downstream of
However, if the wave travel time is small relative the siphon, respectively. The loss coefficient
to the time steps used in the open channel includes all singular and frictional effects. It can
portion, the siphon can be approximated with a be approximated by
lumped parameter. That is, it is only necessary to
treat the siphon as an isolated loss which occurs
with no time delays. The lumped parameter
approximation is valid when

I-<AT where A = cross-sectional area of siphon


a
f = Darcy-Weisbach friction factor
L = .siphon length
where L = siphon length D = equivalent hydraulic diameter of
a = wave velocity siphon = four times hydraulic radius
AT = computational time increment in C, =contraction coefficient of entrant
open channel portion of conduit. flow
C, = expansion coefficient of exit flow
Normally, siphons are ehorter than 2 km and
wave speeds are of the order of 1000 m/s. These
values indicate that the lumped parameter The contraction and expansion coefficients can
approximation is good for computational time be accurately estimated from the ratio of the
increments exceeding 20 seconds in which AT siphon area to the cross-sectional area of the
exceeds the L/a ratio by an order of magnitude. respective water prisms using figure 17.
Ignoring inertial effects, the equation of motion Specifications ‘Constraints
for flow in the siphon is given by
v;’ - v: t y, t z2 + KQ2 (16)
With gate stroking, two boundary conditions are
unknown (not specified) in each pool. These two
2g +y1+z’ - 2g unknowns are reduced to one, at boundaries that
where V = velocity require a continuity of discharge. Based upon
y = flow these two premises, it is possible to construct a
z = distance from datum series of general rules or constraints which
g = local acceleration of gravity define the computational procedures to be
,
K = loss coefficient followed with gate stroking in canal reaches.
.6

.4

-0 .2 .4 .6 .s 1.0
At/As

A. Contraction Coefficient B. Expansion Coefficient

Figure 1 ‘I.-Loss coefficients.

These are as follows: METHOD OF UTILIZATION

1. In a reach there must exist at least one The Gate Stroking Model is intended to be
pool in which two conditions are specified interfaced with two other types of programs. The
(de@t and discharge). This pool is de- first of these is a scheduling program which
fined as the pivot pool. determines the desirable water-surface eleva-
tions and pumping schedules for the daily opera-
2. The computations start at a pivot pool tion of an aqueduct. Two types of scheduling
and progress in both directions (up- programs have been developed for the Central
stream and downstream) away from the Arizona Project. These are the Constant Volume
pivot pool. If the pivot pool is at the end of Model and the Water Power Optimization Model.
a reach, the computations progress in one The Constant Volume Model assumes that some
direction only. The present program does target storage value can be maintained in each
not handle the condition of a pivot pool aqueduct pool during both onpeak and offpeak
located other than at the end of a reach. periods. The Water Power Optimization Model
on the other hand varies the storage value in
3. More than one pivot pool can exist in a each pool to optimize the onpeak and offpeak
reach; however, for the entire reach, ex- pumping schedules.
actly N+l conditions must be specified,
The second type of program with which the Gate
where N = number of pools.
Stroking Model was intended to interface is an
analysis program. The analysis program actually
4. A region 5 computation (fig. 8) must be
simulates the canal flow. For that reason. it is
performed when the specified condition
called the Aqueduct Simulation Model.
is at the far end of the pool from the .
direction in which the computations,are Since the Aqueduct Simulation Model contains
peceeding. detailed flow data concerning the aqueduct at
any point in time, it is used to obtain initial from two to five iterations for convergence,
conditions for the Gate Stroking Model. Then a depending upon the length of time modeled and
schedule of operations is developed for a the number and magnitude of gate movements.
subsequent time interval based on either the Unfortunately, for a canal of any length and for
Constant Volume or the Water Power Optimiza- time periods of the order of a day, computer
tion Model. This schedule of operations consists resource use can be significant. Investigations
of changes in deliveries, possible pump opera- are being continued in an effort to improve the
tions, and specified water level fluctuations convergence of the technique.
within the pools. Using this information, the Gate
Stroking Model determines the gate and, in some
cases, the pump schedules for the time interval PROGRAM DESCRIPTION
under consideration. These schedules include
both the gate openings and timing for the Background
operation of the gates and pumps. Basically, the
changes in flow rates through the pumps and The program (listed in appendix) determines the
gates should be continuous functions. However, gate stroking schedules for a series of pools in
in current practice, all of these changes actually a canal reach. If pools are separated by pumps,
take place in discrete amounts at discrete times. the pump schedules may either be determined or
Therefore, it is necessarytodetermine the effects specified. The program is of a general nature and
of the discrete operations through application may be used to generate gate and, if desired,
of the Aqueduct Simulation Model to the time pump schedulesforvirtuallyanycanal. Modifica-
interval being considered. At the end of the inter- tions were made to the program for the purpose
val the entire process is repeated. The scheduled of expediting a series of operations studies on
time interval is one day in the case of the Granite Granite Reef and Salt-Gila Aqueducts, and
Reef Aqueduct. separate versions of the program were used for
those studies. The program description below
pertains to the general version of the program.
VALIDITY OF THE PROGRAM
The purpose of the program is to determine gate
As a check on the validity of the gate stroking schedules which will produce desired boundary
program, the program was restructured to act as conditions (that is, depths and discharges) in a
an analysis-type program. Several comparisons canal. The differential equations describing the
were made between the Aqueduct Simulation flow apply to free surface flow in a prismatic
Model and the restructured Gate Stroking Model. canal section. Following Wylie [l], the method of
It was found that the restructured Gate Stroking characteristics is used to solve the equations
Model, using the grid of characteristics, con- within the interior of a computational segment.
verges to the steady-state conditions determined Abbott [9] gives an excellent description of the
by the standard step method. Conversely, the significance of characteristic functions as they
method of specified time increments, as used by pertain to the equations of state of water
Shand [B], allows the water surface to drift rather conveyance systems.
than reach steady state. Depth differences from
several tenths to over 0.3 m were found in the Problem Formulation
various comparisons.
The present gate stroking technique has been
The technique for restructuring the Gate formulated such that:
Stroking Model was to make an initial estimate of 1. The computations proceed from one pool to
depths downstream of gates, calculate the state another in an upstream direction.
variables in the pool upstream of the gate using 2. Initial conditions along the entire length
that estimate, calculate the state variables in the of the aqueduct must be specified.
pool downstream of the gate using the upstream 3. At the downstream end of the most down-
values just calculated, recalculate the upstream stream pool:
pool using the new downstream values, etc. This a. The discharge must be specified for all
process was continued until the water surface time.
along the length of the canal stayed within some b. The depth must be specified for all time
convergence limits from one iteration to the next. greater than one wave travel time in the
The technique is stable and generally requires pool.

11
4, At intermediate pools, the depth must be particular pool, that is, the most upstream pool in
specified for time greater than one wave the reach or a pool with no control structure at its
travel time in the pool. downstream end, the schedule will be ignored by
5. To obtain a gate schedule for all pools it is the program. However, schedules must be
necessary to project the most downstream provided for each and every pool. Depths are to be
pool depth and discharge schedules at least provided on the hour, and transitions from one
one wave travel time for the reach into the depth to another are made linearly by the
future. If pump stroking is involved, the program. The specific format for input is
most downstream pool schedules must be described in subroutine READIT, which may be
projected at least one wave travel for the altered to suit the needs of the user.
entire length of the aqueduct.
6. The discharge at all intermediate pools is a Discharge schedules.- must be provided for the
dependent variable; that is, the discharge downstream boundaries of all pools. The dis-
cannot be specified. charge schedule at the upstream boundaryof the
most upstream pool must also be furnished to
Initial Conditions provide consistency with the input/output for-
mats of the Constant Volume and Aqueduct Sim-
The state of the system must be described ulation Models. The program wilt ignore all but
completely at some specific time, usually at a the most upstream and downstream schedules
zero time reference. The system may be in a and will maintain discharge continuity through-
steady-state condition (all time derivatives equal out the length of the reach. Provision is made in
to zero) or in an unsteady, or transient, condition. the program to stroke the pumps at the upstream
The program will accept initial conditions in any end of the reach or to stroke the gate which con-
of three formats, individually for each pool. The nects the reach to a reservoir. In these cases, the
formats differ only in the number of points at appropriate upstream discharge schedules will
which the initial conditions are specified. The be ignored. Discharges are to be specified on the
state variables, x(position along the length of the hour, and transitions from one discharge to
pool), y (water depth), and 0 (water discharge another are made in one step at the time of
quantity), may be specified at zero time at up to40 change.
equally spaced locations along the length of the
pool. Turnout schedules.-must be provided for all
turnouts. Turnouts are considered to be at the
If N = 1, the state variables must be specified at downstream ends of pools. They may occur
either end of the pool. A backwater calcula- separately or in conjunction with a control
tion is then performed by the program to deter- structure such as a check gate or siphon. The pool
mine the state variables at 11 equally spaced number, time of change of discharge (on the
points in the pool. hour), and the change in discharge (positive
discharge is outflow from the canal) are input for
If 2 < N 5 10, linear interpolations of the state up to 48 discharge changes for each turnout. The
variables using the values at x,and f,, only initial turnout discharge is assumed zero, and the
are calculated at 11 equally spaced pomts in specified values are added to the current
the pool. discharge at the specified times. No energy
losses are considered at turnouts.
If 11 < N 5 40, the values of the state var-
Pumps.- at the downstream end of the reach are
iables provided by the user are accepted as the
treated as scheduled discharge devices. At the
initial condition.
upstream end of the reach, pumps are treated as
scheduled discharge devices if pump stroking is
No checks are made at boundaries for con-
not employed.
sistency across boundaries.
With pump stroking, start and stop times for
Boundary Conditions pumps within the reach are calculated using a
dead band for discharge. That is, pumps will be
Depth schedules.- must be provided for the started or shut down when the discharge change
downstream boundaries of all pools. The pro- calculated by the program equals or exceeds the
gram will accept schedules for a 45hour period. incremental discharge change specified for the
If a depth schedule is not appropriate for a pump in the physical descriptors for the pumps.

13
4fte positions. -are determined using the gate Modular Construction
equation which is
The program has been constructed in a modular
Q=C,bB {w fashion. It is anticipated that expanded cap
abilities will be easily accomplished by this
where C,, =
discharge coefficient approach. For example, circular sections were
b =
gate opening added by defining a circular pool type and adding
I3 =
gate width the appropriate code in the area, top width, and
Wf mean velocity in the rectangular
= wetted perimeter function routines.
section upstream of the gate
= depth in the rectangular section output
y,
upstream of the gate
Y2 = depth in the rectangular section The output file is formated to be compatible for
downstream of the gate use as input to the Aqueduct Simulation Model.

Gate movements are performed only when the Development and Structure of the Program
calculated gate position has changed from the
previous position by an amount equal to or The gate stroking program was developed within
exceeding the dead-band value for that gate. Gate rather severe time constraints to perform real-
structures are treated as having one gate. The time operation studies for the Granite Reef and
width of the rectangular gate section is specified Salt-Gila Aqueducts. Previously published ma-
in the variable BOTGATE. Independent operation terial dealt almost exclusively with techniques
of multiple gates within a structure is not a for the solution of the flow equations within a few
feature of the program. Thistype of operation can computational segments. The equations of state
be simulated, however, by modifying the dis- and their solution for a system of pools present a
charge coefficient. The modified coefficient for challenge which is an order of magnitude greater
single gate operation within a structure having than the problem of solving the open channel
multiple gates is given by equations in a single segment. Some comments
cd B,
are appropriate concerning the general structure
c; = and the evolution of the present program.
B
where Bg = bottom width of a single gate The program is written in FORTRAN and is
implemented on a CDC CUBER 74 computer
Inverted siphons. -with or without control gates system. Comment cards”within the program
at the upstream end of the siphons are treated as explain the computational procedures. The
zero length friction effects. The invert drop from assumption is made that anyone using the
the gate structure or siphon inlet to the siphon program will have a good understanding of free
outlet and the friction loss factor are specified for surface flow; an understanding of the various
each siphon. mathematical techniques encountered in the
solution of the equations of state of the system;
Boundaries where no tontrol structure is and a reasonable, but not expert, understanding
present.+that is, changes in cross section, or of digital computing techniques. Within this
turnouts) are treated as points at which nb framework, the program may be utilized to assist
energy losses occur. Depths upstream of these in the investigation of flow characteristics of
boundaries are calculated using the energy virtually any canal configuration.
equation.
The program ,does not contain coding to
Physical Descriptors accommodate all possible structures found in all
canals. It also does not handle all possible canal
The constant descriptors of the reach, that is, cross sections. However, the modular structure
stationing at pool ends, cross-section para- of the program enables the user to add new
meters, gate parameters, roughness, and siphon structure types and pool types in a straight-
parameters, are set in DATA statements in the forward fashion. Likewise, from the multitude of
main program. Comment cards in the program possible operational philosophies for a canal, one
describe the various parameters. was chosen as appropriate for the Granite Reef

13
Aqueduct, and that one was implemented in the across the turnout so the effect of large diver-
program. The philosophy chosen is that pump sions is minimized.
schedules are predetermined at both ends of the
reach and the most upstream pool in the reach As was discussed earlier, various equations have
will be uncontrolled with respect to depth been used to describe the state of a canal system
variations. The same computational techniques at gate structures. Mathematically, flows through
can be used to “float” the downstream pool by the gates can be very sensitive to small
calculating in a downstream direction or to pivot differences in depths on the upstream and
on an interior pool and calculate both upstream downstream sides of the gates. Accordingly,
and downstream of that pool. The modular convergence criteria for the various iterative
construction of the program allows the user to schemes in the program should be examined
implement these schemes in a relatively straight- using the gate parameters particular to the
forward fashion. Provision .was also made for system under investigation. An adequate com-
stroking the most upstream control structures. promise between program run times and long-
The input and output for the program are term stability of the solution must be reached. For
contained in the routines READIT and WRlTElT self-consistency, identical methods of computing
and may be modified easily to suit the user. the discharge through gates should be used in
both the gate stroking and the Aqueduct
Various techniques were investigated for the Simulation Model programs. For this reason,
treatment of the system at control structuresand equations more closely approximating the
turnouts. The generalized program evolved out of equations in Shand’s [8] program were used in
experience gained with these various tech- the Granite Reef and Salt-Gila versions of the
niques. Some observations relative to the gate stroking program. Due to the finite number
evolution may be beneficial. of gate movements and different computational
techniques, differences will occur between the
Since turnouts generally divert only a small predicted and “actual” water surfaces and
portion of the total discharge, the program used discharges, even with identical physical param-
for the Granite Reef study merely adjusted the eters and gate equations. These differences
discharge within a computational segment when should, however, be small. Large differences if
a characteristic crossed the turnout location. This they occur, can be attributed to interpolation
treatment appears adequate when the turnout errors in one or both of the programs.
diversion is small and when the turnout is not
near the control structures at the ends of the Program Details
pool. The advantage of this technique is that depth
continuity is maintained at a turnout location for The equations upon which the coding is based
all time. lf the turnout discharge were a large are described in the preceding sections. The
percentage of the total discharge, a simultaneous friction slope is computed using English cus-
depth correction could be made at the turnout tomary units. The conversion to SI metric units
using the energy equation. Difficulties in making can be accomplished by deleting the constant
the proper discharge correction were en- 1.49 in Subprogram S and by changing the value
countered when the turnout fell between a of the gravitational constant GGRAV to metric
boundary and the grid point upstream or . units.
downstream of the boundary. Since the method The purposes of the subroutines and function
of characteristics is stable for rather large time subprograms are as follows:
increments and, consequently, the x spacing of PROGRAM GSM - Main calling program. Canal
grid points can be rather large, the chance for parameters are input through DATA state-
a turnout between a boundary and the next grid ments located in this program
point is fairly high. Thus, for the Salt-Gila study SUBROUTINE READIT - Routine to read in dis-
and the general program, turnouts are treated charges, depths, turnout schedules, and
as boundaries. The disadvantage of this approach initial conditions
is that depths at the turnout for one wave-travel SUBROUTINE STROKER - Routine to solve gate
time in the upstream pool are not calculated stroking problem for free surface flow
identically in the upstream and downstream SUBROUTINE SOLVER - Routine to solve for X, t,
pools. Apparently, serious discrepancies do not 0 and y at point P using values at point R on
occur for normal operation schemes. The energy positive characteristic and values at points on
equation is used to calculate a corrected depth negative characteristic
FUNCTION C - Function to calculate wave celer- Division, American Society of Civil Engineers,
ity Vol. 100, No. HY3, p. 405-424, March 1974.
FUNCTION R - Function to calculate hydraulic [4] Chow, V. T., Open-Channel Hydraulics,
radius McGraw-Hill, p. 328, 1959.
FUNCTION A - function to calculate cross-sec- [5] Ippen, A. T., “Channel Transitions and Con-
tional area in trapezoidal, horseshoe, or trols,” Engineering Hydraulics, edited by
circular channels H. Rouse, p. 536, 1950.
FUNCTION TW - Function to calculate top width [6] Metzler, D. E., “A Model Study of Tainter-Gate
in trapezoidal, horseshoe, or circular channels Operation,” M.S. Thesis, Iowa State Univer-
FUNCTION S - Function to calculate energy slope sity, Ames, 1948.
using Manning’s equation [7] U.S. Army Engineer Waterway Experiment
FUNCTION P - Function to calculate wetted peri- Station, “Hydraulic Design Criteria,” Chart
meter in horseshoe, trapezoidal, or circular 320-8, 1959.
channels [8] Shand, M. J., “Final Report - Automatic
SUBROUTINE BOUNDARY - Calculates t, 0, and Downstream Control Systems for Irrigation
y at an upstream boundary and t and y at a Canals,” Hydraulic Engineering Laboratory,
downstream boundary College of Engineering, University of Cali-
FUNCTION HDN - Function to determine down- fornia, Berkeley, Report HEL-8-4, August 1971.
stream boundary depth [9] Abbott, M. B., An Introduction to the Method
FUNCTION QDN - Function to determine dis- of Characteristics, American Elsevier, New
charge boundary condition at downstream York, 1966.
end of pools
FUNCTION QUPST - Function to specify dis-
charge at upstream end of most upstream pool
SUBROUTINE TRNOUT - Routine to specify turn-
out discharge
SUBROUTINE GATEMO - Routine to calculate
gate openings and pump schedules
SUBROUTINE GATEY - Routine to calculate
depth in gate section from channel depth
FUNCTION CD - Routine to calculate discharge
coefficient as a function of gate opening to
upstream depth ratio
SUBROUTINE FOLOWER - Routine to follow
stroking motions to determine start times and
final positions of gates and pump schedules
SUBROUTINE WRITEIT - Routine to merge turn-
out, gate, and pump schedules into input file
for unsteady model

BIBLIOGRAPHY
.!
[l] Wylie, E. B., “Control of Transient Free-
Surface Flow,” Journal of the Hydraulics
Division, American Society of Civil Engineers,
Vol. 95, No. HYl, p. 347-361, January 1969.
[2] O’Loughlin, E. M., “Application of Unsteady
Flow Analysis to Operation Decisions in Long
Aqueducts,” International Commission on
Irrigation and Drainage, Eighth Congress,
Bulgaria, R. 16, p. 28.2.207-28.2.223 1972.
[3] Gientke, F. J., “Transient Control in Lower
Sacramento River,” Journal of the Hydraulics
APPENDIX

Gate Stroking Computer Program Listing


,
.1
PROGRAM GSM 74/74 OPT= 1 FIN 4.6+426 79/61/99. 14.41.14 PAGE 1

1 PROGRAM GSM(1NPU1=64,0UTPU1=64,TAPE1~64,TAPE2=54,1APE3~64,
TAPE4=64,TAPE5=64.TAPE6r64,TAPE6~64,TAPE9=64)
C GENERALIZED GATE-STROKING PROGRAM

5 THE GATE-STROKING TECHNIOUE IS TAKEN FROM BEN WYLIE5


PAPER CONTROL OF TRANSIENT FREE-SURFACE FLOW IN
JAN., 1969, PROC OF ASCE, JOURNAL OF THE HYDRAULICS
DIVISION

10 THIS PROGRAM IS IN ENGLISH UNITS


THE. GRAVITATIONAL CONSTANTr’GGRAV, AND THE FACTOR
1.49 IN MANNINGS N IN SUBPROGRAM S ARE THE ONLY
CHANGES REGUIRED FOR METRIC CONVERSION
SOME INTERMEDIATE OUTPUTS MAY BE AFFECTED BY A
15 METRIC CONVERSION, AS FORMAT F10.3 IS GENERALLY
USED, SO IT WOULD PAY TO CHECK AND SEE THAT NO
SIGNIFICANCE WAS LOST IN ROUNDING

THE PROGRAM IS SET UP TO INPUT PUMP SCHEDULES AND


20 DEPTH SCHEDULES AT 0, 0700 AND 2300.HOURS.
THE INPUT IS IN ROUTINE READIT, AN0 MAY HAVE TO BE
CHANGED. THE OUTPUT IS IN ROUTINE WRITEIT, AND MAY
HAVE TO BE CHANGED TO SUIT THE USERS REOUIREMENTS.

25 TAPE1 IS INTERMEDtATE T.0 AND Y FILE AT POOL BOUNDARIES


TAPE2 15 INTERMEDIATE T AND GATE OPENING FILE
TAPE3 15 INPUT FILE CONTAINING PUMP, DEPTH AND TURNOUT SCHEDULES
TAPE4 IS SCRATCH FILE USED TO OVERLAY COMMON
TAPES IS INTERMEDIATE START TIME AND GATE POSITION FILE
30 TAPE6 IS OUTPUT FILE CONTAINING PROGRAM MESSAGES TO USER
TAPE6 15 INPUT FILE CONTAINING INITIAL CONDITIONS
IN UNSTEADY MODEL FORMAT
TAPE9 15 OUTPUT FILE IN UNSTEADY MODEL FORMAT FOR INPUT
TO UNSTEADY MODEL
35
COMBON /ALL/
PHYSICAL PARAMETERS
. NMANN.SO,GGRAV,PI,6OTTOM(50),SIDSLDP(5O),RADIUS(5O).6DTGATE(51),
MANN(SO),6OTGRAD(5O),SYFONDY(5O),SYFNLOS(5O),5PEED(5l),PUMPW
40 : (51~.0E06AN0(51),NGATES,GATEC0(11,51),051NVT(50)~U51NVT(51),
WORKING VARIABLES
. T~N.XEND,TINC,NXINC,NPOOLS,IPOOL,TYP~L,STABEG(5O),
STAEND(SO),Tl,Yl,STUDY,
CANAL DESCRIPTION
45 . UPSTA(~~)~DWNSTA(~O),NOPOOLS,POOLTYP(~O),STRTYPE(~~)
. TMAXI.NUMTRNS,USRES,DSRES

REAL ‘NYANN ,MANN


XNTEGER TYPPOOL,PDOLTYP,STRTYPE,STUDY
50
SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SWROUTINES
TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE 6OUND
COMMON /SCRATCH/
N0X(5D),XSTART(40,50),YSTART(40,50),0START(40,50),
55 : HSCHED(46,5O),XINIT(4O),YINIT(4O),OINIT(4O),H(46),
DUMUVl(7169)
’ DIMENSION SCRATCH(lS7BT)
PROGRAM GSY 74/74 OPT=1 FTN 4.0+420 79/01/09. 14.41.14 PABE 2

EDUIVALENCE (SCRAfCH(l).NDX(l))

60 : STUDV TELLS THE PROGRAM IF THE MOST UPSTREAM


C STRUCTURE, EITHER A PUMP OR A GATE, IS To BE
C STROUED OR NOT
C STUDY=1 IS NO STROKING (D IS SCHEDULED)
C STUDY=2 IS STROKING
65 DATA (STUDV=2)
C

6 POOLS ARE NUM6ERED FROM DOWNSTREAM TO UPSTREAM IN THE REACH


C ARRAYS ARE DIMENSIONED FOR 50 POOLS
70 C ARRAYS ARE DIMENSIONED FOR 40 GRID POINTS PER POOL
C
C SET MAX TIME FOR THE REACH
C MAX TIME MUST BE AT LEAST STUDY TIME (USUALLY 24 HOURS)
C PLUS WAVE TRAVEL TIME FOR THE REACH
76 C TIME IS IN SECONDS
DATA (TMAXI=19600.)
C STATIONING AT UPSTREAM ENDS OF POOLS
C ARRAY IS UPSTA(POOL)
DATA (uPSTA=270494.4.234654,4,196105.6,162146.6,134270.4
80 .45+0.0)
C STAiIONING AT DOWNSTREAM ENDS OF POOLS
C ARRAY IS OWNSTA(PWL)
0;:; ~0tj~STA.304972.6,2704D4.4.234654.~.196105.6,162140.6

85 c N&ER DF’PDDLS IN THE REACH


DATA (NOPOOLS=l)
C TYPES OF POOLS
POOLTVP=O THIS POOL NOT USED
: POOLTVP+l TRAPEZOIDAL CHANNEL
90 POOLTVP=2 HORSESHOE TUNNEL
: PDDLTVP13 CIRCULAR TUNNEL
/y C ARRAY IS POOLTVP(POOL)
DATA (PWLTVP=S*1,45*0I
C MANNINGS N
95 C ARRAV IS MANN( POOL) ‘.
DATA (MANN=Li*0.016,45*0.0)
C 6DTTDM SLOPE - SO
C ARRAY IS BOTGRAD(PWLI
DATA (60TGRA0~5+0.00006.45*0.0)
100 C STRUCTURE TVPES AT UPSTREAM ENDS OF POOLS
ARRAV IS STRTVPE(PWL)
: STRTVPE(S1) IS 11051 DDUNSTREAM STRUCTURE
STATVPE=O POOL NOT USED
: SlRlVPE’l PUMP
101 SlRTVPE=2 NORMAL CHECK GATE
: SIRTYPE= SYPHDN WITH CHECK GATE
C STRlVPEF4 SVPHDN WITH No CONTROL STRUCTURP
C STRTVPEFS No CONTROL STRUCTURE
DATA (STRTVPE=3*2,3,2,45*0,2I
110 C TRAPEZOIDAL CHANNEL PROPERTIES - BOTTDM WIDTH AND SIDESLOPE
C SIoESLoPE=O.O GIVES A RECTANGULAR SECTION
C
DATA (6DTTOM=S*24.0,45+D.9)
DATA (SIOSLOP=S*1.5,46*0.0I
.4
PROGRAM G5M 74/n OPT=1 FTN 4.6+428 70/0l/09. (4.41.14 PAG5 3

115 C CHECK GATE PROPERTIES


C POOL NUMBERING 15 O/5 TO U/5
C GATE IS AT U/S END Of POOL
; PARAMETERS FOR GATE Af,O/S END OF REACH IN ARRAy(51)

120 C SET UPSTREAM AN0 DGWNSTREAM CONSTANT LEVEL


C RESERVOIR OEPTHS IF GATE5 AT U/S OR O/S END5 Of REACH
DATA (USRES=lB.O),(DSRES-14.0)
C
C CHECK GATE BOTTOM WIOTH
125 DATA (601GATE=5+36.0,45*0.0,36.0)
C GATE MOTOR SPEED
DATA (SPEED=5*0.75,45+0.0,0.75)
C GATE MOTION DEADBAND
DATA (DED6ANO=5*0.5,45*0.0,0.5)
130 C INVERT DROP PROM GATE SECTIDTt TO CHANNEL ON D/S SIDE OF GATE
C Xf SYPHON. LUMP INVERT DROPS ON D/S SIDE IN SYFONDY BELOW
If GATE AT D/S END OF REACH, DEPTH IN DSRES IS IN
GATE SECTION - No o/s INVERT DROP DR sypHffl PARAMETERS ARE
C INPUT FOR THAT GATE
135 DATA (DSINVT=3*O.T,O.O.O.l,45*O.O)
C INVERT DROP CROM CHANNEL TO GATE SECTION ON U/S SIDE OF GATE
C If GATE AT U/S END OF REACH, DEPTH IN USRES IS IN
C GATE SECTION - NO U/S INVERT DROP IS INPUT FOR THAT GATE
C U/S INVERT DROP 15 ALLOYED FOR GATE AT D/S END
140 C Of REACH, IF ANY
DATA (USINVT=50+0.0,0.0)
C GATE COEFFICIENTS
C ARRAY X5 GAlECO(l-ll.PODL)
C SEE ROUTINE GATEMO FOR THEIR USE IN GATE EOUATIDN
145 DATA (GATECOP
. .71..606,.67,.655,.645,.644,.659,.693.
. .756..662,1.019
. ,.71..666,.67,.655,.645,.644,.659,.693,
. .756..862.1.019
150 . ..71,.666,.67,.655,.645,.644,.659,.693,
. .756, .662.1.019
. ..71..666,.67,.655,.645,.644,.659,.693,
. .756..662,1.019
. ..71..666,.67,.655,.645,.644,.659..693,
155 . .756..662.1.019
. I495~0.0~~
. ..71..666,.67,.655,.645,.644,.659,.693,
.756..662,1.019)
C PUMP PARAMETERS
160 C PLMPOG(POOL) IS OtSCHARGE INCREMENT FOR PUMPS AT U/5 END Of PGDL
C USED IN ROUTINE FOLOWER TO DETERMINE PUMP SCHEDULES
C PUMPDfJ(61) IS AT D/S END OF REACH
DATA (PUMPDG~4+0.0,125.,45+0.0,125.)
C
165 C TUNNEL PROPERTIES
C RADIUS IS HEIGHT Of TUNNEL, INVERT TO CREST, FOR HOR5ESHGE TUNNEL
C RADIUS IS RADIUS Of TUNNEL FOR CIRCULAR TUNNEL
c ARRAY IS RADIUS(POOL)
DATA (RADIUSp50+0.0)
170 C
C SYPHON PARAMETER5
PROGRAM GSM 74/74 OPT= 1 FTN 4.6+420 7o/01/09, 14.41614 PAGE 4

C ARRAYS ARE SYFONDY(POOL) AND SYFNLDS(PDOL)


C INVERT DROP IN SYPHDN
DATA (SYFDNDY=3*0.0,3.24,0.0,4S*O.G)
,175 C SYPHDN LOSSES
DATA (SYFNLOS*3*0.0,3.21E-7,0.0,45*D.D)
C TURNOUTS
C TURNOUTS CAN BE Al D/S END OF ANY POOL
C NUMBER OF TURNOUTS IN THE REACH
100 DATA (NUMTRNS=O)
C
DATA fGGRAV=32.2),(PI=3.14lS927)
C
REWIND 6
10s REWIND 1
REWIND 2
REWIND 3
REWIND 4
REWIND 5
190 REWIND B
REWIND 9
C
C READ IN FLOWS, DEPTHS AND TURNOUT SCHEDULES
C AND INITIAL CONDITIONS
195 CALL READIT
C
C SET TIME INCREMENT FOR CALCULATIONS AT BOUNDARIES (IN SECONDS)
C 5 MIN X30 HOURS ,STAYS IN ARRAY LIMITS
C
J
200 TINC’300.
3
C SET NUMBER DC POOLS IN THIS REACH
NPDOLS=NOPODLS
C SET BEGINNING AND ENDING STATIONING FOR POOLS
Do 10 I=l.NPDOLS
205 STABEG(I).UPSTA(I)
STAEND(I)*DWNSTA(I)
CONTINUE-
C’:ET MAXIMUM CALCULATION TIME FOR THE REACH -
C -DEL TIME, NOT COMPUTE TIME
210 TMAX=fhlAXI
C SOLVE FOR EACH POOL
C CALCULATIONS PROCEED IN UPSTREAM DIRECTION
DO 100 I=l,NPOOLS
IPooL=I
215 C SET POOL LENGTH AND TYPE
XEND=STAEND(IPOOL)-STABEG(IPODL)
TYPPODL=PDDLTYP( IPDOL)
C SET INITIAL CONDITIONS FOR THIS POOL
C HSCHED IS D/S DEPTH SCHEDULE
220 C HSCHED WAS READ IN BY READIT
C SCHEDULE IS FOR 48 HOURS
w 20 J8T.40m
H(J)=HSCHED(J,IPDOL)
CONTINUE
225 C’iET NUMBER OF INCREMENTS ALONG LENGTH OF POOL
C CALCULATE BACXWATER CURVE IF ONLY ONE POINT IS GIVEN -
c THE NON-ZERO VALUE FOR y MAY BE AT EITHER END OF POOL
. C AND THE BACKWATER CURVE WILL BE CALCULATED TO THE OTHER END
PROGRAM GSM 74/74 OPT= 1 FTN 4.6+426 79/01/09. 14.41.14 rr AGB 5

NXINC=NOX(IPOOL)
230 C CALCULATE BACKWATER CURVE, IF DESIRED, FOR 11 POINTS IN POOL
IF(NXIN~.EQ.2.AND.YSTART(l,IPOOL).EQ.0.0) CALL EAKWA,rR(l)
IF(NXINC.EQ.2.AND.YSTARl(2.IPOOL).EQ.O.O) CALL BAKWAtR(2)
C LINEAR WATER SURFACE IF LESS THAN 11 POINTS
C IF BAKWATR WAS CALLED, NXINC WAS SET TO 11
235 IF(NXINC.GE.11) GO TO 40
WRITE(6.9000)
OELX=XSTART(NXINC,IPOOL)-XSTARl(l,IPOOL)
OELY=YSTART(NXINC,IPOOL)-YSTART(l,IPOOL)
DELQ=OSTART(NXINC,IPOOL)-OSTART(l,IPOOL)
240 C LINEAR INTERPOLATION FROM ZERO TO XENO IN 10 EQUAL STEPS
00 30 L=l.ll
FL=L-1
XSTAl?l(L.IPOOL).XSTARi(l,IPOOL)+FL*OELX/lO.O
YSTART(L.IPOOL)=YSlART(1.IPOOL)+FL*DELY/lO.O
245 tJ~~~~~;~,IPOOL)=QS’ART(l,IPOOL)+FL*DELQ/lO.O
30
NXINC=ll
NDX(IPOOL)=il
CONTINUE
250 C4:“ECK THAT ARRAYS ARE NOT OVERRUN
WRITE(6.9010) NXINC
IF(NXINC.GT.40) STOP
C MEASURE X FROM ZERO AT UPSTREAM END OF POOL
DO 50 L=l.NXlNC
255 XINIf(L)=XSTART(L,IPDOL)-STABEG(IPOOL)
YINIl(L).YSTART(L,IPOOL)
QINIl(L)=Q5TART(L,IPOOL)
50 CONTINUE
WRITE(6.9040) XINIT(l),XlNIT(NXINC)
250 XINIT(l)=O.O
XINIT(NXINC)=XEND
C SAVE COMMON ON TAPE4
WRITE(4,9030) (SCRATCH(ISCR),ISCR=l,B45ti
REWIND 4
265 C SOLVE THIS POOL
CALL STROKER
c RETRIEVE corwo~ FOR &XT POOL
;:~t;;,f030) (SCRATCH(ISCR), ISCR=I ,645O)

270 WRIlE(6.9020) IPOOL


100 CONTINUE
REWIND 1
C CALCULATE GATE OPENINGS AND PUMP DISCHARGES FOR 24 HOUR PERlOD
CALL GATEMO
275 C CALCULATE GATE MOVEMENTS FROM GATE OPENINGS AND CALCULATE PUMP
CALL FOLOWER
C OUTPUT CALCULATED DATA FOR UNSTEADY MODEL INPUT
CALL kRrTE IT
WRrTEfU;9050)
260 REWIND 6
REWIND 9
STOP
C FORMATS
9000 FORMAT(lX.43HINTERPOLATING WATER SURFACE FROM END POINTS)
265 9010 FORMAT(1X,6HNXINC=,IlO,l9li MAX ALLOWED IS 40 )
PROGRAM GSM 74/74 OPT= 1 FtN 4.6+429 79/01/09. 14.41.14 PAGE 6

9020 FORMAT(lX,9HENO POOL,IS/)


9030 FORMAT(5020)
9040 FORMAT(lX.6HX1NIT=.2FlO.2)
9050 FORMAT(lX,19HNORMAL TERMINATION )
290 EN0

-1
SUBROUTINE BAKdATR 74/74 OPT= 1 FTN 4.6+420 79/01/09. 14.41.14 PAGE 1

SUBROUTINE BAKUAll?(ICHOOSE)
C ROUTINE TO CALCULATE INITIAL BACKWATER CURVES
C GIVEN DOWNSTREAM Y AND 0 (ICHOOSE=l)
C OR UPSTREAM Y AND 0 (ICHOOSE=2)
C SECOND ORDER RUNGA-KUTTA METHOD ON
C DY/DX=(SO-S)/(l-Q**2*l/A**3*G)
DIMENSION X(40)
COMMON /ALL/
C PHYSICAL PARAMETERS
10 NMANN,SO,GGRAV,PI,9OTTOM(5O),StDSLOP(SO),RADIUS(SO),9OTGATE(Sl
: -MANN(SO),BOTGRAO(50),SYFONDY(SO),SYFNLOS(5O),SPEED(5l),PUMPW
(51),OEOBAND(5l),NGATES.GATECO(ll,5l),DSINVT(5O),USINVT(5l),
C WORKING VARIABLES
. TMAX.XENO,fINC,NXINC,NPOOLS,IPOOL.tYPKK)L,STABEG(5O),~
IS STAENO(SO),Tt ,Yl ,STUDV,
C CANAL DESCRIPTION
. UPSfA(50),DWNSTA(SO),NOPOOLS,POOLTYP(5O),STRTYPE(5l),
TMAXt.NUMTRNS,USRES,OSRES
c *
20 REAL NMANN ,MANN
INTEGER TVPPOOL,POOLTYP,STRTYPE,STUDY
C
C SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES
C TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
25 COMMON /SCRATCH/
N0X(50).XSTART(40,50),YS1ART(40,50),0START(40,50),
: HSCHEO(4B,50),XINIT(4O),YINIT(4O~,OtNIT(4O~,H(4B~,
DUMMYl(7l89)
’ REAL K
30 WRITE(6.70)
C OIVIDE POOL INTO (N-l)*10 SEGMENTS
c STORE ONLY EVERY TENTH VALUE OF XN ~140 VN
XEN0=XSTART(2,IPOOL)-XSTART(l,IPOOLt
N=ll
35 S0=BOTGRAO(IPOOL)
TYPPOOL=POOLTYP( IPOOit
NOX( IPOOL)-11
._-._..-
urtur.11 ..
IF(ICHOOSE.E0.i) GO TO 5
40 C CALCULATE IN DOWNSTREAM DIRECTION
OO=OSTART(l,IPOOL)
YO=YSTART(1,tpOOL)
FNFN- 1
DX=XEND/(FN+10.)
45 X(1)=0.0
L=O
VN=YO
00 2 1=2,N
L=L+l
50 00 1 J=l.lO
FJ=J
XN=X( Lt+FJ+OX
Al=A(yN)
TWlpTW(YN)
55 tF(TWl.LE.0.) W TO 60

{~;:;%i!:‘~ TO 40
SUBROUTfNE BAKYATR 74/74 OPT= 1 FTN 4.6+426 79/01/0@. 14.41.14 PAGE 2

Vl =QO/Al
IF(A6S(Vl*“‘-Vc2).Lf.0.00l) GO TO 50
60 Rl=R(YN)
Sl=S(Vl.R P)
Fl~(SO-Sl)/(l.-QO**2*tW1/(Al+*3+GGRAV))
Y2=YN+PX*Fl
A2=AIY2)
65 TW2=TW(Y2)
IF(TWP.LE.0.) GO TO 60
VC2=GGRAV+A2/TW2
IF(VCP.LE.0.) GO TO 40
V2*QO/A2
70 IF(ABS(V2*+2-VC2).LT.O.O01) GO TO 56
R2=R( Y2)
S2=Sl.V2,R2)
F2=(SO-S2)/(l.-QO+*2+TW2/(A2**3+GGRAV))
K=O.S*(Fl+Fl)
7s YNsYN+bX*K
1 CONTINUE
XI L+l )+XN
XjiARi( L+l , IPOOL)=XN+STABEG( IPOOL)
YSTART(L+l ,XPOOL)=YN
60 QSTART(L+l.IPOOL)=09
2 CONTINUE
XSTARl(l.IPOOL).STA6EG(~POOL)
XSTART(N.IPOOL)=STAENlMIPOOL)
RETURN
65 C CALCULATE IN UPSTREAM DXRECTION
6 CONTINUE
QO~QSTART(2,IPOOL)
YO=YSTART(2,IPOOL)
FNtN- 1
90 OX*-XENO/( FN+lO.)
X(N) =XENO
L=N+l
tN=YO
00 20 ts2.N
95 L=L-1
00 10 J=l,lO
F&J
XN=Xl L)+FJ+OX
Al=A( YN)
100 TWlmTW(YN)
IF(TWl.LE.0.) GO TO 60
VC2.GGRAV+Al/TYl
IF(VC2.LE.0.) GO TO 40
Vl;QO/Al
105 1F(A6S(V1**2-Vc2).LT.0.901) GO TO 56
Rl=Rl YN)
Sl=S(vl .A1 )
F1=(SO-Sl)/(l.-0O*+2+TWl/(Al**3*GGRAV))
Y2riN+OX+ii
110 A2=Al Y2)
TW2=.lW(Y2)
IF(TWP.LE.9.) GO TO 60
VC2*GGRAV+Al/TY2
IF(VC2.LE.0.) GO TO 40
-1
SUBROUTINE BAKWATF 74/74 OPT.1 -0” - .a - FTN 4.5+420 79/01/09. 14.41.14 PAGE 3

115 V2*00/A2
IF(ABS(V2**2-VC2).LT.0.001) GO TO 50
R2=R(Y2)
S2=S(V2.R2)
F2=(SO-S2J/(l.-Q0**2*TW2/(A2*+3*GGRAVJJ
120 K=O.S*(Fl+FZJ
YN=YN+OX*K
10 CONTINUE
X( L-l J=XN
XSTART( L-l , IPOOL)=XN+STABEG( IPOOL)
125 YSTARl(L-1 ,IPOOLJ=YN
QSTARf(L-l,IPOOL)=QO
20 CONTINUE
XSTARf(l,IPOOL)~STABEG(~POOL)
XSTAffT(N,IPOOLJ~STAEND(IPOOLJ
130 YSTART(N,~POOL)=YO
QSTARl(N,IPOOL)=QO
RETURN
40 CONTINUE
WRITE(6.45) YN,XN
135 45 FORMAT(lX.22HZERO OR NEG V++2 AT Y99F10.3,
3~ X=,FlO.B,12H IN BAKWATER J
STOP
50 CONTINUE
WRITE(6,55) YN,XN
140 55 FORM4T(lX.l4HCR1TICAL DEPTH,F10.3,3H AT,F10.3,
. 12H IN BAKWATER J
STOP
60 CONT t NUE
WRITE(6,65) YN,XN
145 65 FORMAT(lX,32HtERO OR NEG TM !N BAKWATER AT Y=,F10.3,
3H Xr.Fl0.3)
70 ’ FbRMAT(lX.4lHCALCULATING BACKWATER CURVE FOR THtS POOL
STOP
-_.-
LNP
SUBROUTINE READ! 1 74/74 OPT= 1 FtN 4.5+428 79/01/00. 14.41.14 PAOE t

1 SUBROUTINE READIT
C ROUTINE TO READ IN FLOWS AND DEPTHS AND TURNOUT SCHEDULES

5 COMMON /ALL/
C PHYSICAL PARAMETERS
NMANN.SO.GGRAV.Pt,BOTTOM(5O),SIDSLOP(5O),RADtUS(5O),BOTOATE(5l),
: MANN(50),BDTGRAD(5O),SYFONDY(SO),SYFNLOS(5O),SPEED(5l),PUMP~
. (51~.DEDBAND(5l).NGATES.GATECO(11,~1),DStNVT(5O),UStNVT(51),
10 C WDRKING VARIABLES
. TMAX.XEND.T~NC.NXINC.NPODLS.~PODL,TYP~L,STABEG(~O),
. STAEND(SO).Tl,Yl,STUDY.
C CANAL DESCRIFTIDN -
. JPS~A(~~~.DWNSTA(~~).NOPOOLS.POOLTYP(~O),STRTYPE(~~),
15 TMAXI;NUMTRNS,USRES;DSRES
c -
REAL NMANN ,MANN
INTEGER TVPPDDL,PDOLTYP.STRTYPE,STUDY

20 :
C SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES
C TO KEEP MEMORY WITHIN SDME SORT DF REASONABLE BOUND
COMMON /SCRATCH/
NDXc50),XSTART(40,50),YSTART(40,50),0START(40,50),
25 : HSCHED(48.50),
XtNtT(40),YtNtT(4O),GtNtT(4O),H(4B),
4 : DD(48),W(48),07(48.50),
. DUMNVl(4673)

30 : TAPE 3 HAS UPSTREAM AND DDWNSTREAM 0 SCHEDULES, DEPTH


C SCHEDULES AND TURNOUT SCHEDULES AT
C 110. 0700 AND 2200
C SCHEDULES MUST BE SUPPLIED AT EVERY POOL BDUNDARY,
C EVEN WHEN THEV CANNOT BE USED, E.G., AT POOL BOUNDARtES
35 C WHERE ND CONTROL STRUCTURE IS PRESENT
C
REWIND 3
C READ OVER 4 LINE HEADING
AEAD(3.9000)
40 MAXFNOPDOLS
CREADATHOURSO,7,22
C READ U/S PUMP SCHEDULE
C INDEX iS HOUR+1
READ(3.9030) OU(l),OU(8),OU(23)
45 BACKSPACE 3
C STORE IN D/S TO U/S ORDER
C TAPE 3 IS IN U/S TO D/S ORDER
tNDX=MAX
Do 30 I=l,MAX
50 READ(3.9040) HSCHED(1,tNDX),HSCHED(B,tNDX),HSCHED(23,tNDfl)
INDX= INDX-.l
30 CDNTINUE
C READ D/S PUMP SCHEDULE
C READ AT HOURS 0. 0700 AND 2200
55 C INDEX IS HDUR+l-
READ(3.9050) OD(l),OD(‘B),OD(23)
C FILL IN REMAINING HOURS OF DAY
SUBROUTINE READIT 74/‘14 OPT= 1 FTN 4.6+426 79/01/09. 14.41414 PAGE a

DO 40 1x2.7
QUlI)=QUlI-1)
60 QDtI)=QDlI-1)
DO 40 J=l,MAX
MSCHED~I.J)=HSCHED(I-I,~)
40 CONTINUE
DO 50 1=9,22
66 QU(I)=QU(I-1)
QDll)=QD(I-1)
DO 50 J=l,MAX
HSCHED~I,IJ)=HSCHED(I-l,sJ)
50 CONTINUE
70 QU(24)=QUl23)
QD(24)*00(23)
DO 60 J=l.MAX
HSCHEOC24,J)=HSCHED(23,~)
CONTINUE
75 C%NTINUE SCHEDULES FOR SECOND 24 HOURS
DO 70 111.24
QUlI+24)=QUl24)
QD(1+24)=90(24)
DO 70 J=l.MAX
60 HSCHED(I+24,J)=HSCHED(24.J)
CONTINUE
C’:T(HDuR.pOOL) IS TURNOUT Q
C TURNO;T6;S1A: ;;S END OF POOL

65 DO 60 J:l:46
QT(J.I)+O.O
CONTINUE
C?RE THERE ANY TURNOUTS IN THIS REACH
C NUYTRNS IS SET IN A DATA STATEMENT IN MAIN PRGQRAM
90 IF(NUMTRNS.EQ.0) GD TO 140
C READ TURNOUT SCHEDULES
C MAXTURN IS LAST INDEX IN THIS REACH
C READ OVER 3 LINE HEADING FOR SCHEDULE
READ1 3.9060)
9s MAXTURN*NUMTRNS
DD 130 INDX.1 ,MAXTURN
C READ NUMBER OF CHANGES AT THIS TURNOUT
READ(3,SOlO) LPGDL,NCHANGE
C READ TIME OF CHANGE AND DELTA Q
100 C CHANGES MAY ONLY BE MADE DN THE HOUR FROM 0000 TO 2!300 WOUAS
DG 120 I=l,NCHANGE
REA013,9070) NTIYE,DELO
C TIME INDEX IS HOUR+1
tTIYE=NTIME+l
106 QT(ITINE,LPOOL)=QT(ITXME,LPGOL)+DELQ
MINFITIME+l
IF(MIN.GT.24) GG TO 100
C SET 0 AT TURNOUT FOR THE REMAtNDER OF THE DAY
DO 90 JaMIN,
110 QT(d.LPDOL)=QT(d-1,LPOOL)
90 CONTINUE
100 CONTINUE
C CONTINUE SCHRDULE FOR SECOND 24 HOURS
DO 110 ‘d=l,24
SUBROUTINE READIT ?4/74 OPT= t FtN 4.6+428 79/@!jOS. 14.41.14 PAQE 3

11s OT(Jt24:1POOL)~OT(24,lPooL)
110 CONTINUE
120 CONTINUE
130 CONTINUE
140 CONTINUE
120 TAPE 8 HAS INPUT FOR UNSTEADY MODEL
READ INITIAL CONDITIONS FROM THIS FILE
COPY TAPE E TO TAPE 9 FOR USE AS INPUT TO
UNSTEADY MODEL WITH GATE SCHEDULES INSE :RTED
AFTER THIS PROGRAM HAS CALCULATED THEM
125 AND WITH TURNOUT SCHEDULES INSERTED
REWIND 8
REWIND 9
READ 3 TITLE CARDS AN0 STUDY START TIME CAR0
Do 150 1X1.4
13c REAO(E.SOeft)
KRITE(9.9080)
YRITE(6rSO80)
lb0 CONTINUE
C READ INITIAL CONDITIWS
13s KPOOL=NOPOOLS
C READING POOLS FROM UPSTREAM TO DOWNSTREAM
C STORING POOLS FROM D/S TO U/S
C X VALUES RUN FROM U/S TO D/S IN BOTH CASES
C MUST HAVE CARDS FOR U/S AND D/S ENDS OF ALL POOLS
140 C IF ONLY U/S AND D/S POINTS ARE GIVEN, THEN
IF Y U/S l 0.0, CALCULATE BACKWATER CURVE FROM b/S END
: IF Y O/S = 0.0, CALCULATE BACKWATER CURVE FROM U/S END
1F Y U/S AND V D/S .NE. 0.0, LINEAR WATER SURFACE AND
: LINEAR 0 CHANGE BETWEEN ENOS OF POOL
145 C
170 XNDX= 1
READ(B.9110) VSTART(XNDX,KPOOL),,OSTARt(tNOX,KPOOL1
,XSTART(INDX,KPOOL)
’ YR1TE(9,9110) YSTART(INDX.KPOOL),OSTART(INDX,KPOOL)
150 .XSTART(INOX,KPOOL)
C SEA&” FOR END OF POOL
XENO=DWNSTA(KPOOL)
100 INDX=INDX+l
REA0lE,SllO) VSTART(INDX,KPOOL),QSTART(INDX,KPOOL)
155 .XSTART(INDX,KPOOL)
‘URXTEl9,9110) YSTART(INDX.KPOOL),OSTART(tNDX,KPOOL)
,XSTART(INDX,KPOOL)
IF(XSTART(INDX,KPOOL).EQ.XENO) GO TO 190
IF(XSTART(tNDX,KPOOL).LT.XEND) GO TO 180
160 WRITE(6.9120) XSTART(INDX,KPOOL),XEND
STOP
190 CONT I NUE
C STORE THE NUMBER OF 1NITXAL CONDITION POfNTS
NDX(KPOOL)-INDX
165 C HAVE WE READ DOWN TO POOL I, THE MOST b/S POOL tN tHE REACH
!F(KPOOL.EO.l) RETURN
KPOOL=KPOOL-1
C READ NEXT D/S POOL
W TO 170
170
:
SUBROUTXNE READ11 74/?4 OPT. 1 FtN 4.6+426 79/01/09. i4.41.14 PAQE 4

C FORMATS
9000 FORMAT(lX///lX)
9010 FORMAf(I4,IS)
175 9030 FORMAT(l7X,Fl1.2,2F16.2 1
9040 FORM4T(lX/l7X,F11.2,2Fl6.2/lX)
9050 FORMAT( 17X,F11.2,2Fl6.2 )
9060 FDRYAf(lX//lX)
9070 FORMAT(f9,3X,F10.2)
160 9060 FORMAT(40H
9090’ F:i:AT(IlO) 1

9100 FORMATtlOH ,110,2OH


40H 1
165 9110’ FORMAT(Fl0.3,lOH IF10.3,
30H ,F10.2)
9120’ FORMAf(lX.13HFDUNO StAtION,Fl0.3,
. 21H BEFORE FINDIN XEND=,FlO.B,lOH IN READIt )
END
SUBROUTINE STROKER 74/74 OPT= 1 FtN 4.9+420 79/01/09, i4.4l.14 PA09 1

1 SUBROUTINE STROKER
C ROUTINE TO SOLVE GATE STROKING PROBLEM FOR OPEN
C CHANNEL FLOW AS PER WYLIE
C THIS PROGRAM IS SET UP FOR Q AN0 Y
5 C
COMMON /sOLV/ XP,YP,~P,QP.XR,YR,TR,QR,XS,YS,TS,QS
C
COMMON /ALL/
C PHYSICAL PARAMETERS
10 . NMANN,SO.OGRAV.PI.BOtfOM(5O),SIOSLOP(5O),RAOIUS(5O),BOtaAtE(I1),
. MANN(50),BOTGRAO(5O),SYFONOY(5O),SYFNLOS(SO),SPEEO(51),PUMP~
(5l).OEOBANO(51),NGAfES,OITECO(11,51),OSINVT(50),USINVt(S1).
C WORKING VARIABLES
. fRAX.XENO,TINC,NXINC.NPW~S,IPOOL,TYPPOOL,StA9LQ(SO). _ . .
15 STAENO(50),Tl,~l,StU~y,
C CANiL DESCRIPTION
. UPSf*(50),OWNSTA(SO),NOPOOLS,POOLTYP(S0),STRTYPE(9l),
. TYAXI,NWTRNSrUSRES.OSRES
C
20 REAL NMANN ,MANN
INTEGER TYPPOOL,POOLTyP,STRTYPE,STUOY
C

: SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES


29 C TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
CGMNON /SCRATCH/ XCP1(40),YCP1(40),1CP1(40),QCP1(40)
. .XCM1(40).YCM1(40),TCM1(40),QCM1(40)
. ,XCP2(40),YCP2(40).fCP2(4O),QCP2(40)
. .XON2(4O),YON2(4O),TON2(4O),QON2(40)
30 . ,XCF3(4O),YCP3(40),TCP3(4O),QCP3(40)
. ,XCM3(4O),YCM3(4O).TCM3(4O),QCM3(40)
. ,XON3(4O),YON3(4O),TON3(4O),QON3(40)
. ,XUF3(4O),YUP3(4O),TUP3(4O),QUP3(40)
. ,XCP3B(40),YCP3B(40),TCP3B(40),QCP38(40;
3.? . ,XON3B(4O),YON3B(4O),TON3B(4O),QON3B(4O)
. .XUP~~(~O),YUP~B(~O),~~P~B(~O),QUP~B(~O)
. .XCP4(Bl),YCP4(Bl).TCP4(Bl),QCP4(Bl)
. .XCM4(Bl),YCM4(Bl),TCM4(9l),QCN4(8l)
. .XUP4(B1),YUP4(Bl),TUP4(Bl),QUP4(Bl)
40 . .XCP5(Bl),YCPS(Bl),TCP5(8l),QCP5(Bl)
.OU~~NYl(5394),XlNIT(40).Y1NIT(40),QINIT(40),OUMMY2(2544)
: ,INOX(B5),TINOX(BS),KOWN,TOWN(5OO),QOWN(5OO),YOWN(SOO)
. .KUF,TUP(500),OUP(5OO),YUP(5OO)
.KMAX.SAVETUP(!lOO).SAVEYUP(5OO),SAVEQUP(500)
45 ’ DIMENSION ISEONCE(5,2),X(BS),Y(B5),T(BS),Q(BB)
C ISEONCE IS SEQUENCE Tb CALCULATE SOLUTION REGIONS

: MAP OF SOLUTION REGIONS FOR ALL BUT MOST U/S POOL


C REGION 38 IS REPEATED AS NECESSARY TO GET TO TYAX
50 C Q AND Y ARE GIVEN AT D/S BGUNOARY AND Q AND Y ARE
C CALCULATED AT U/S BGUNOARY
l
E: . +
C+ * +
55 . +
g: l *
C+ l +
. .

SUBROUTINE STROKER 74/74 FTN 4.6+428 79/0~/09. 14.41.14 PAOE 2

l l
:r (I
60 fit*+ l
C+ l
+ l
:r l
+
65 :r *
C+ +
C+ l l
*
EL .
70 C+ .
. l
I: .
l
75 :r l
l
:: .+
4 l
:r .
*
80 :: *
C+ l
.
:: l *
C+ l
85 c*+ l
c ++**++~,**++++*++*++**********~************

: MAP OF SOLUTION REOIONS FOR MOST U/S POOL


C IF 0 INTO POOL IS SCHEDULED (STUDY=i)
90 C
ii 0 1s GIVEN AT BOTH BOUNDARIES AND Y 1S CALCULATED
C FOR BOTH BOUNDARIES

L .
95 l
.: : +
*
:: l
l
100 :: +
l
;: .
C* .
.
105 E: I l
l l
:r l l
C+ l l
l 2 l
110 ;r l l .
. l l
I: l 1 . l
C**
c .*+**~*+++***.+*.*++~**~~***~~************~
SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.8+428 79/01/09. 14.41.14 PAGE 3

115

c ISEQNCEI~ EXITS THE cokwT~TIoN ~oof~


DATA (!SEQNCE=1,2,3,4,~,1.1,5.818)
C LAST=2 FOR MOST UPSTREAM POOL
120 C IF Q IS SCHEDULED FOR U/S END OF THAT POOL
C FLOAT DEPTH FOR MGST UPSTREAM POOL
C IF Q IS SCHEDULED
LAST=1
C IPOOL IS THE POOL BEING CALCULATED
125 C NPOOLS IS THE NUMBER OF THE MOST U/S POOL
IF(IPOOL.EQ.NPOOLS.AND.STUDY.EQ.l) LAST=2
C TIMINC IS THE TIME INCREMENT USE0 FOR CALCULATIONS
C IT MAY BE INCREASED IF 40 POINTS IN REGION 3 YONT GET
C THE SOLUTION TO X.0
130 C OR DECREASED IF POOL IS VERY SHORT
TIMINC=TINC
C SET THE NUMBER OF CWPUTATION POINTS FOR EACH REGION
Nl=NXlNC
N2=Nl
135 N4rNl+N2-1
NS=N4
C KUP IS UPSTREAM INDEX FOR STORING T,.Q AND Y
KUPsN4
IF(LAST.EQ.2) KUP=l
140 C KDii W;OrSTREAM INDEX FOR STORING T, 0 AND Y

DO 9iCl IPASS*l,S
ISEQ=lSEQNCE(IPASS,LAST)
GO 70(100,200,300,400,500,900,900,900),1SEQ
146 C X, Y. T AND Q ARE WORKING ARRAYS
C . ..CPJ AND . ..CMJ ARRAYS ARE X, Y. T AND 0 VALUES
C ALONG + AND - CHARACTERISTICS IN REGION J
c ’
C REGION 1 - INITIAL CONDITIONS
150 100 CONTINUE
C XINIT. ETC.. ARE INITIAL CONDITIONS AT TfME l 0
00 110 I=1 ,Nl
x(I)=XINIT(:)
Y(I)=VINIT(I)
155 T(I)=O. ’
Q(I)=QINIT(:)
110 CONTINUE
X(N1 )=XEND
XCPl(O=X(l)
160 vcPl(l)=Y(l)
TCPl(1 )=T( 1)
QCPl ( l ) =Q( l )
XCMl(l)=X(Nl)
YCMlfl)=Y(Nl)
165 TCMl(l)=T(Nt)
QCMl(l)=Q(Nl)
MaNI -1
DO 160 1=2,Nl
DO 150 J=l ,M
170 C IP, IR AND IS ARE P. R AND S POINT INDICES IN X, Y, t
C AND Q WHERE P IS POINT TO BE SOLVED FOR,
. .

l . .

SUBROUTINE STROKER 74/?4 OPT= 1 FTN 4.6+429 79/01/09. 14.41.14 PAGE 4

C R IS KNOWN POINT ON + CHARACTERISTIC AN0


C S IS KNOWN POINT ON - CHARACTERISTIC
C P, R. S NOTATION FOLLOWS WYLIE
175 IP=J
IR=J
IS=J+1
C SET R AN0 S POINT VALUES
XR=XI IR)
190 YR=Y( IR)
TR=T( IR)

x ::;
YGY( IS)
185 TS+T( IS)
OS=O( IS)
C SOLVE FOR POINT P
CALL SOLVER
C STORE NEW VALUES
190 X( IP!=XP
Y(IP)=YP
T(IP)=TP
O(IP)=OF
C STORE EN0 POINTS IN REGION 1 BOUNDING CHARACTERISTfCS
195 IF(J.EO.1) GO TO 130
120 IF(J.EG.M, GO TO 140
GO TO 150
C STORE IN + CHARACTERISTIC
130 CONTINUE
200 XCPl(I)=XP
YCPl( I)=YP
TCPl(I)=TP
GCPI f I ) =OP
C LAST POINT IS ALSO ON - CHARACTERISTIC
205 C CHECK FOR THIS SITUATION
GO TO 120
C STORE IN - CHARACTERISTIC
140 CONTINUE
XCMl(I)*XP
210 YCMl(I)=YP
TCMl(I)=TP
OCMl~I)=OP
150 CONTINUE
M=M- 1
215 160 CONTINUE
WRITE~6,9000)ISEG
GO TO 900
c REGION CONTINU;
2 - SPECIFIED D/S
200
220 XDN2l l)=XCMl(l)
YDN2(l)=YCMl(l)
TDN2(t)=TCRll(l)
ODN2~l,=OCMl(l)
DO 230 I=P,Nl
225 x(l)=xcMl(I)
Y(l)=YCMl( I)
T(t)=TCMl(I)
0(1)=OC~l(I)
SUBROUTINE STROKER ?4/74 OPT= 1 FTN 4.0+428 T9/01/0@. 14.41.14 PAW 5

M=I
230 Do 210 J=2.M
IR=J-1
IS*d
IP8J
C SET R-AN6 S POINT VALUES
235 XR=X( IR)
YR=+i IRj
lR=t( IR)
OR-01 IR)
C CHECK FOR BOUNDARY
240 IF(J.EO.Mj GO TO 220
XS=Xf IS)
YS=Y( IS)
fS=l( IS)
OS=O( IS1
245 C SOLVE‘FOli POiNl P
CALL SOLVER
C STORE NEW VALUES
X(IP)=XP
Y(IP)=YP
250 T(IP)=TP

%z%
XP=XEND
C BNDRY( 1) IS- D/S BOUNDARY - 0 SPECIFIED, C+ CHARACTERISTIC
255 CALL eNDRY
c STORE IIOUNDARY tiA,iUEs
X(IP,=XP
Y(IP).YP
7( IP)=lP
260 O( IPJ=OP
C STORE IN D/S ARRAY
XDN2(l)=XP
YDN2(1)=YP
TDN2(I)=TP
265 ODN2 I I ) =OP
230 CONTINUE
C STORE IN REGION 2 BOUNDING CHARACTERISTIC
DO 240 I=1 ,N2
XCP2(1)-X(I)
270 YCP211)=Y(I)
TCPPf I)=T( I)
OCP2(1)=0(1)
240 CONTINUE
c STORE IN TDwN.oDwN,YDWN FOR GATE MOTION CALCULATIONS LATER
I 275 DO 250 I-l .NP
TDwN(I)=lDN2(1)
OOwN~I)=OON2(I)
YDwN(I)=YDNP(I)
250 CONTINUE
280 C UPDATE D/S INDEX FOR SAVING T, 0 AND Y
UW=N2
C STORE ENDING TIME AND DEPTW FOR INTERPOLATION
C IN POOL 2 AND FOLLOWING
Tl=TDN2(N2)
285 Yl=YDNS(N2)
l .

SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.6+428 79/01/09. 14.41.14 PAGE 6

URITEf6.9OOO)ISEQ
GO TO 600
C REGION 3 - H AND 0 SPECIFIED D/S
300 CONTINUE
290 SET 1 PREVIOUS FOR STORING BOUNDARY VALUES ONLY
EVERY TINC TIME ADVANCE
TUPPREV=O.
TDNPREV=O.
COMPARE WAVE TRAVEL TIME WITH TIMINC TO SEE IF
295 POOL IS TOO SHORT TO USE TIMINC
YCHECK=YDNZ(NZ)
CELERTY=C(YCHECK)
FN=Nt-1
DXCHECK=XEND/FN
300 DX3=CELERTY*TIMINC
IF(DX3/DXCHECK.GT.Z.) TIMINC=DXCHECK/CELERTY+lr
IF(DX3/DXCHECK.Gl.Z.) WRITE(6.9040) TIMINC
IF TIMINC TOO SHORT, TREAT AS ZERO LENGTH POOL
IF(TIMINC.LT.O.S*TINC) 00 TO 810
305 NUMBER OF POINTS IN 38 REGION
N30=0
NCROSS CHECKS THAT REGION 3 CROSSES X-0
NCROSS=O
MATCH WITH REGION 2 END POINT ON D/S BDUNDARY
cd
310 XDN I( t )=XEND
U TDN I( 1 )*TDNZ(N2)
YDN II 1 )=YDNZ(N2)
QDN II 1 )=ODNZ(NZ)
Do 104 182.40
315 N3=
XDN I( I)=XEND
TDN l(I)=TDNJ(I-l)+TIMINC
IF( ‘DN3(I).GE.TMAX) GO TO 302
YDN I(I)=HDN(TDN~(I))
320 ODN II I)=GDN(TDN3(1))
GD ‘0 304
302 TDN I( I)=TMAX
YDN I( I )=HDN(TMAX)
I( I )=GDN(TMAX)
325 2 ‘0 306
304 CON INUE
N3B 1
306 CON INUE
YRITE(6~9010) N3
1
,330 IF(TDNJ(l).GE.TDN3(N3)) WRITE(6,9020)TDN3(1),TDN3(N3)
I : TDN3(1).GE.TDN3(NB)) STOP
311)=XON3fN31
ycP3fij=YDr43iNjj
i TCP3(l)=TDN3(N3)
_ 1 33s GCP3(l)=ODN3(N3)
XCM3(1)*XDN3(1)
YCM3(1)=YDN3(1)
TCWB(l)=TDN3(1)
OCM3(1)=0DN3(1)
340 DO 308 I=1 ,N3
X(I)=XDN3( I)
Y(I)=YDNB( I)
SIJ6ROUfINE STROKER 74/74 OPT= 1 FTN 4.6+426 7Q/Ot/OQ. 14.41.t4 PAGE 7

T(I)=TDNI( I)
Q(I)=QDN3(1)
345 306 CONTINUE
C K=UPStREAM X=0 IFitERSECt INDEX
K=l
M=N3- 1
N=M
350 00 322 I=1 ,Y
C NALLX CHECKS THAT ALL XP ARE LESS THAN ZERO
NALLX=O
DO 320 d=l ,N
I P=J
355 IS-J
IR=,t+l
C SET- R AN0 S POINT VALUES
XR=X( IR)
YR=Y( IR)
360 TR=f( IR)
QR=Q( IR)
XS*X( IS)
YS=Y( IS)
TS=T( IS)
365 QS*Q( IS)
C SOLVE FOR POINT P
CALL SOLVER
C SET NALLX IF ANY XP HAS NOT CROSSED X-0
IF(XF.GE.O.0) NALLX=l
370 C CHECK FOR CROSSING X-0
1F(XF.LT.O.O.AND.XS.GE.O.O) GO TO 316
310 CONTINUE
C STORE NEW VALUES
X(IP)=XP
375 Y(IP)=YP
T(IPl=TP
Q(IPb=QP
C CHECK FOR END POINTS
IF(J.EQ.l) GO TO 314
360 312 1FtJ.EQ.N) GO TO 318
GO TO 320
C STORE IN C-
314 CONT I RUE
XCM3T I+l).XP
365 YCM3Tl+l)=YP
TCY3T I+1 ).TP
QCM3( I+1 )=QP
C LAST POINT MAY ALSO BE ON + CHARACTERISTIC
C CHECK FOR THIS SITUATION
390 GO TO 312
C STORE IN C+
316 CONTINUE
XCP3l x+l)*XP
YCPIT 1+1 )=YP
395 TCP3t I+l)=TP
QCP3(?+l)=QP
GO TO- 320
C INTERPOLATE UPSTREAM VALUES
C Al X=0

/--- -.

.
su0~OUtlNE c*ofiKER 74/74 OPT= 1 FTN 4.6+428 79/01/09. 14.41.14 PAOE 8

400 318 CON1 INUE


C SET NCROSS TO SHOW THAT X-0 WAS CROSSED
NCROSS*l
DELX=XP-XS
DX=XP
4OS RAlIO=DX/DELX
XUP3(K)=O.O
YUP3(K)=YP-RATIO+(YP-YS)
TUP3lK)=TP-RATIO+(TP-TS)
p;kV;K)=PP-RAT’“(OP-QS)
410
G; TO 310
320 CONTINUE
C HAVE ALL POINTS CROSSED X=0 - IF SO STOP CALCULATtNQ
IF(NALLX.EQ.0) GO TO 324
415 N-N-1
322 CONTINUE
324 CONTINUE
C IF DID NOT CROSS X=0, STOP
tF(NCROSS.EQ.0) GO TO 346
420 C FIND C+ X=0 INTERCEPT .
DO 326 1=2,N3
IF(XCP3(I).GT.O.O) GO TO 326
DELX=XCP3(1)-XCP3(1-1)
DX=XCP3(1)
425 RATIO=DX/DELX
XUPBlK)=O.O
YUP3~K)=YCP3(~)-RAtIO*(YCP3(I)-YCP3(t-l))
TUP3~K).lCP3(I)-RAT1O*(TCP3(I)-TCP3(l-l))
430 ~l$‘3~K~;~CP3(t)-RATIO*(QCP3(I)-QCP3(1-1))

326 CONTINUE
326 CONTINUE
C SGRT UPSTREAM BOUNDARY ON T
DO 330 1.c.n
43s T(I)=TUPB(I)
TINOX(I)=l.
330 INDXf I)=0
ILOW=l
TLOW=T(l I
440 DO 336 J=l ,K
DO 332 181 ,K
~F(TINDX(I).LT.O.O) GO TO 332
IF(TfI).GE.TLOW) GO TO 332
I LOW= I
445 TLOW=T(I)
332 CONTlNUc
1NDXfJ)~ILOW
TINDX(ILOW)=-1.
DO 334 L-1 ;n
450 IF(TINDX(L).LT.O.O) GO TO 334
ILOW=L
TLOW=T( L)
GO TO 336
334 CONTINUE
45s 336 CONTINUE
DO 336 I=1 ,K
sumoutxNE STROKER 74/74 OPT=1 FTN 4.6+428 79/01/09. 14.41.t4 PAN 0

X(x)rXUP3( 1)
Y(I)=YUp3(11
T(I)=TUPI(I)
460 O(I)=OUP3(1)
338 CON1 INUE
DO 340 1.l.K
J=INDX(I)
XUPI( I )=X(J)
465 YUPO( I )=Y(J)
TUPJ( I )=T(d)
OUP3l I bO(J)
CONTINUE
C3”s:ORE IN TDWN,ODWN,YDWN FOR GATE MOTION CALCULATIONS LATER
470 DO 342 I=l,N3
C UPDATE D/S INDEX TO SAVE 1, 0 AND Y
KDWN=KQWN+t
J=ltDWN
TDWN(J)=TDN3(1)
47s ODWN(J)=OON3(1)
YDWN(J)=YDN3(1)
IF(I.EO.l) GO TO 342
C SAVE ONLY AFTER TINC TIME ADVANCE
IF(TDN3(I).GE.TONPREV+TINC) GO TO 341
490 KDWN=KDWN-1
GO TO 342
C UPDATE TDNPREV (T PREVIOUS)
341 TDNPREV=TON3(f)‘O.D1
342 CONTINUE
405 UR~TE(6r9000)ISEO
C STORE IN TUP. OUP AND YUP FOR FUNCTIONS OON AND HDN FOR NEXT POOL
DO 344 i=i,K
C UPDATE U/S .INDEX TO SAVE T, 0 AND Y
KUP=KUP+l
490 N=KUP
TUP(N)*TUPB(I)
OUP(N)=OUpB(I)
YUP(N).YUP3(1)
IF(I.EO.1) GO TO 344
495 C SAVE ONLY AFTER TINC TIME ADVANCE
IF(TUP~(I).GE.TUPPREV+TINC) GO TO 343
KUP-KUP-1
GO TO 344
C3y;DATE TUPPREV (T PREVIOUS)
500 TUPpREV=fUP3(1)-0.01
CONTINUE
C’::ECK TO SEE IF 38 REGION MUST BE USED TO GET TO TMAX
IF(N3B.EG.0) GO 10 809
C INITIALiZE FOR 38 REGION
505 XDN3B( 1 )=XEND
TDN3fl( l).lDN3(N3) .
YDNBO( 1 )=YDN3(N3)
ODN3B(l)=ODN3(N3)
C START.CALCULATIONS IN 38 REGION
510 GO TO 350
340 CONTINUE
WRITE(6.9030) 1POOL
C INCREASE TIME INCREMENT TO GET ACROSS POOL 1N 40 POINTS

_~-- .----- * .- ---..

.
SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.8+428 79/01/09. 14.41.14

TIM1NC=TIMINC+TINC
515 IF(N3.Lt.49) STOP
WRIfE(6.9040)tIMINC
GO TO 300
C 8 ALCULATE 36 REGION
30 CONTINUE
520 NAL ,X=0
DO 154 112.40
N30 I
XDN IB(I)=XEND
TON 1El(I)=lDN38(1-l)+TIMINC
525 lF( DN3B(I).GE.TMAX) GO TO 352
YDN 161 I)=HDN(TDNBB(I))
IB!I)=QDN(TDN~~(I))
:iF 0 354
352 TO 135( I).TMAX
530 YDN ll3(1)=MDN(TMAX)
QDN let I )=QDN(TMAX)
GO ‘0 356
3s4 CON ‘INUE
356 ‘INUE
535 WRI ‘E(6,9050)N38
IF( ‘DN38(1).GE.TDN3B(N38)) WR:TE(~,9060)TDN3B(~),TDN3B(N38)
IF( ‘DN3B(l).GE.TDN3B(N3ll)) STOP
XCP lB( l)=XDN3B(N38)
YCP 1st l)=YDN3B(N3B)
540 TCP IB( l)=TDN3B(N3B)
lB( 1 )=QDN3B(N3B)
zip I58 I=1 ,N3B
X(1 =XDN3E(I)
Y(I =YDNBE( I)
545 T(1 =:rDN3B( I)
OfI =QDNBB(I)
358 CON ‘INUE
C K=UPSTRE M X=0 INTERSECT 1NDEX

550 z::, l-l


N=N IO-1
I68 I-1 ,M
EL .X=0
L-N Is+1
555 Do 166 J-1 ,N
L-L -1

::: :-1
560 C RET Rx& i 5 POINT VALUES
XR=
YR= :;::;
‘(*RI
;::
X5= 11 :t;
‘( 15;
rs: *(Is)
I( IS)
c 50L”EQf: t POINT P
CAL . SOLVER
SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.6+426 79/01/09. 14.41.14 PAGE 11

C SET NALLX IF ANY XP HAS NOT CROSSED X=0


XF(XP.GE.O.0) NALLX=l
C
C CHECK FOR CROSSING X=0
575 XF(XP.LT.O.O.AND.XR.GE.O.0) GO TO 364
360 CONTINUE
C STORE NEW VALUES
X( IP)=XP
Y(lP)=VP
560 T(IP)=TP
V(IP)=VP
C CHECK FOR END POINTS
IF(L.EV.N36) GO TO 362
GO TO 366
565 C STORE IN C+
362 . CONTINUE
XCP3B(I+l)=XP
VCPIB( x+1 )=YP
TCPBB(I+l)*TP
590 a&Pg’;;;“oP
C INTERPOLATE UPSTREAM VALUES
c AT X=0
364 i0Ntl~uE
595 DELX=XP-XR *
DX=XP
RATIO=OX/DELX
XUP3B(K)=O.O
wP3B(K)=VP-RATIO+(VP-YR)
500 TUP3B(K)=TP-RATlO*(TP-TR)
y;(K)=VP-RATlO*(VP-OR)

GO TO 360
366 CONTINUE
.w C HAVE ALL POINTS CROSSED X-0 - XF SO, STOP CALCULATIONS
XlL(NALLX,EV.9) GO TO 379
IF(l:.GE:N3.) N=N-1
iF(~iGE.Y3j GO. TO 366
x(l-)ixcP6(.1+1)
610 v(l).YcP3(1+1
T( 1 )=fCP3( +l
v(1)=vcF3( f *lb 4
366 CONT rNUE
370 CONTINUE
815 C SORT UPSTREAM BOUNDARV ON T
K=K-1
DO 372 I=l,K
T(I)=TUP36(:)
TINDX(I)=1.
620 372 INDXI f )=O
ILOW=l
TLOW=l(l)
DO 376 J=l ,K
DO 374 I=1 ,K
625 tF(TINDX(I).LT.9.9) GO TO 374
:F(T(I).GE.TLOW) GO TO 374
.LOW=X
SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.6+426 79/01/09. 14.41.14 PAGE 12

TLOW=T(I)
374 CONTINUE
630 INDXIJl=ILOW
TINDX(ILOW)=-1.
DO 376 L=l,K
IF(tINDX(L).LT.O.O) GO TO 376
ILOW=L
636 fLOW=T(L)
GO TO 376
376 CON1 INUE
376 CONTINUE
DO 380 I=l,K
640 x(I)=xuPse(I)
Y(I)=YUP36(1)
T(I)=lUP3B(I)
OlI)=OUP36(1)
360 CONTINUE
645 DO 382 111 ,K
J=INDX(I)
XUP3B(I)=X(J)
YUP3E(I)=Y(J)
TUP3B( I )=1(J)
650 OlJP3R(I)=O(J)
362 CONTINUE
c ‘STORE IN TDWN.ODWN,YDWN FOR GATE MOTION CALCULATIONS LATER
DO 304 I=l,N36
C UPDATE O/S INDEX FOR SAVING T, Q AND Y
656 KDWN=KDwN+l
J=KDWN
TDWNI J)=TDNIB( I)
ODWN(J)=ODN3B(I)
YDWN(J)=YDNBB(I)
660 IF(I.EO.l) GO TO 364’
C SAVE ONLY AFTER TINC TIME ADVANCE
IF(TDN3B(I).GE.TDNPREV+TINC) GO TO 363
KDWN=KDWN-1
GO TO 364
665 363 TDNPREV=lDN36(1)-0.01
364 CONTINUE
WRITE(6.9OOO)ISEO
C STORE IN TUP, OUP AND YUP FOR FUNCTIONS ODN AND YDN FOR NEXT POOL
DO 386 I.l,K
370 C UPDATE U/S INDEX FOR SAVING T, 0 AN0 Y
KUPsKUP+l
N=KUP
TUP(N)=TUPJB(I)
OUP(N)=OUP36(1)
675 YUP(N)=YUPBB(I)
IF(I.EO.l) GO TO 366.
C SAVE ONLY AFTER TINC TIME ADVANCE
IF(TUP3B(I).GE.TUPPREV+TINC) GO TO 365
KUP=KUP-1
660 GO TO 386
365 TUPPREV=TUP3B(I)-0.01
366 CONTINUE
C SAVE TIME FOR MAX-FOR NEXT POOL
SAVTIME=TUP(KUP)
SUBROUTINE STROKER 74/74 OPT= 1 FtN 4.6+426 iS/6l/O!#. 14.41.14 PAM 13

695 C STORE IN C+3


00 389 1-t .N3
XCPB( I)=XCP3B( I)
YCP3lI)=YCP3B(I)
TCP3(I)=lCP36(X)
690 QcP3(t)=DCP36(t)
366 CON1 I NUE
IF(tCP36(l).GE.TMAX) GO TO 600
C CLEAR FOR ANOTHER PASS THRU 36 BAN0
00 390 1st ,40
695 XUP3B(I)=O.
YUPOB(I)=O.
TUP39( 1) =O .
OUP30(1)=0.
XON3B(I)*O.
700 YON36( t)=O.
TDNOB(I)=O.
OON38( I)=O.
390 CONTINUE
XDN39( 1 )=XEND
705 YDN3E(l)*YCP38(1)
TDN3E(l)=TCP36(1)
w$y;;;DcP36(1)

C REGION 4 - NO BOUNDARY CONDITlDNS


710 400 CONT I WE
C LOAD C+ FROM REGION5 1 AND 2
P M=N2
P DO 410 I=l,Y
J=N2-t+l
715 XCP4(1)=XCP2(J)
YCP411)=YCP2(J)
TCP41 I )=TCP2(J)
OCP4(l)=DCP2(J)
410 CONTINUE
720 N=N4
J=Nl
DO 420 I=M,N
XCP4( I )=XCPl(J)
YCP4(Il=YCPl(J)
725 TCP4f I)=lCPl(J)
ocP4lI)+ocPl(J)
J=J- I
420 CONTINUE
C INTERPOLATE C-AT SAME X AS C+
730 XCM4ll )=XCP4(1)
Ycwfl)=YCP4(t)
TCM4( t J=TCP4(1)
~~a;‘;yP;‘t 1

735 DO 430 J:2:N3


tF(XCM3(d).GT.XCP4(1)) 00 IU 430
DX=XCP4(1)-XCW3(J)
DELX=XCM3(d-I)-XCM3(d)
AATIO=DX/OLLX
740 XCrns(I)=XCP4(!)
YCM~(I).YCM~(IJ)+RATIO*(YCY~(J-I)-YCMB(d))
..-.-..---__-- ___. -_..-..~ ~.._.. --

SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.6+420 79/01/09. 14.41.14 PAGE 14

745 430
440 CONTINUE
CK- INDEX FOR UPSTREAM BOUNDARY
K=l
DO 450 I=l,N4
750 X(I)=XCP4(I)
Y(I)=YCp4( I)
T(I)=lCPS( I)
D(I)=DCp4(1)
450 CONTINUE
755 c SET x=0, T=O UPSTREAM
XUP4(1)=0.0
YUP411 )=YCP4(N4)
TuP4I 1 )=TCP4(N4)
OUP4(t)=DCp4(N4)
760 DO 400 1*2,N4
X(l)=XCM4(I)
Y(l)=YCM4(1)
T(l)=lCM4(1)
O(l)=OCM4(I)
765 DO 470 J=2,N4
IP=J
IR=J-1
IS=J
C SET R AN0 S POINT VALUES
770 XR=X( IR)

r:=;l ::j
OR=O( IA)

775 E’CI is;


& 15)
OS=O( IS)
C SOLVE FOR POINT P
CALL SOLVER
760 C STORE NEW VALUES
x(IP)=xP
Y(IP)=YP
T(IP)=TP
o(IPJ=op
765 C CHECK FOR X=0 CROSSING
IF(XP.Ll.O.O.ANO.XR.GE.O.O) GO TO 464
GO TO 470
c INTERPOLATE AT X=0 ON C+
460 CON1 INUE
790 K=K+t
0X80. o-XP
DELX=XR-XP
RAtIO=DX/OELX
XUP4lK)=O.O
795 YUP4(K).YP+RAlIO+(YR-YP)
TUP4(K)=lP+RAtIO*(lR-lp)
QUP4~K)=OP+RAfIO*(OR-Op)
470 CONTINUE
SUBROLJtlNE STROKER 74/74 OPT= 1 FTN 4.6+428 79fU111B. 14.41.14 PAOE 15

480 CONTINUE
000 WRITE~6.9000)ISEO
C STORE IN TUP, OUP AND YUP FOR USE IN FUNCT1ONS ODN AN0 HDN IN NEXT
DO 490 I=1 ,K
TUP(l)+TUP4(1)
ouP(I)=ouP4(I)
005 YUP(I)=YUP4(1)
490 CONTINUE
GO TO 800
C REGION 5 - 0 SPECIFIED ON BOTH BOUNDARIES
500 CONTINUE
010 C SET TPREV FOR STORING BOUNDARY VALUES ONLY EVERY
C TINC TIME ADVANCE
TPREV=-(TINC+O.Ol)
C LOAD C+ FROM REGIONS 1 AND 2
DO 510 I*l,Nl
015 XCP51I)=XCPl(I)
YCPS(I)=YCPl(I)
TCPS(I)=TCPl(I)
ocP5~I)=ocPl(I)
510 CONTINUE
020 M=Nl
DO 520 I=1 ,N2
XCPStM)=XCPZ(I)
YCPSlM)=YCP2(1)
TCPS(H)=TCP2( I)
025 fyl%~““ocP2”’

520 C;NT INUE


00 530 I=1 ,NS
X(X)=XCPS( I)
030 Y(I)=YCP5(1)
T(I)=TCPS(I)
O(I)=OCP5(1)
530 CON1 I NUE
C STORE EN0 POINTS
935 540 CONTINUE
YUP(KUP)=Y(l)
TUP(KUP)=T(l)
OUP(KUP)=O(l)
YOWNtKDWN)=Y(NS)
040 lDWNlKDWN)=T(NS)
ODWN~KOWNl~O(N5)
C STOP AFTER. TMAX
IF(TtNS).GE.TMAX) GO TO 590
C STORE ONLY EVERY TIME T HAS ADVANCED BY TINC
a45 C TO KEEP NUMBER OF POlNlS IN YUP, YDWN, ETC.
C APPROXIMATELY THE SAME &S IN THE OTHER POOLS
DELTAt= T(l)-TPREV
IF(OELTAT.LT.TINC) GO TO 550
TPREV=T(l)
850 KUP=KUP+l
KDWN=KDWN+l
550 CONTINUE
DO 506 I=1 ,N5
IP=I
055 IR-I-l
. .

-.- ~.- .~~_ ----_ .--

SUBROUTINE STROKER 74/74 OPT= 1 FtN 4.6+426 79/01/09. 14.41.14 PAGE 16

IS=I+l
IF(I.EO.l) GO TO 560
IF(I.EO.NS) GO TO 570
C SET R AN0 S POINT VALUES
660 XR=X( IR)
YR=Y( IR)
TR=T( IR)
QR=Gf IR)
XS=X(IS)
065 YS=Yl IS)
fS=T( SS)
OS-0 rsi
C SOLVE FOR POINT P
CALL SOLVER
070 C STORE NEW VALUES
X(IP =XP
Y(IP =YP
f( IP =TP
O(IP )=OP
675 GO TO 560
C UPSTREAM BOUNDARY
660 X%X( IS)
YS=Y( IS)
TS=Tl IS)
060
~sp’KS’
CA;L’BNDRY(O)
C
ii STORE NEW VALUES
065 X(IP)=XP
Y(IP)=YP
f( IP)=TP
O( IP)=OP ,
GO TO 560
a90 C6pSfREAM BOUNDARY

KI :t;
TR:,, IR)
QR=Gt IR)
095 XP=XENO
CALL BNORY(1)

: STORE
NEW VALUES
X( IP)=XP
900 Y(IPI=YP
T(IP)=lP
O(IPb=OP
560 CONTINUE
GO TO 540
905 590 CONTINUE
WRITE(6.9000)ISEO
600 CONTINUE
GO TO 900
C TIME INCREMENT FOR CALCULATIONS TOO SHORT
910 C FOR REALISTIC RUN TIMES, SO TREAT AS ZERO
C LENGTH POOL
610 CONTINUE
SUBROUTINE STROKER 74/74 OPT= 1 FTN 4.B+42B 79/01/09. 14.41.14 PAQE 17

WRIfE(6.9120).
C CALCULATE 0/S VALUES AT TINC XNTERVAL5 FROM END OF REQION 2
915 C TO TIME TMAX
TP=TDN2(N2)-TINC
00 820 I=N2,500
TP=TP+TINC
IF(TP.Gf.TMAX) TP=TMAX
920 TOWNlI)=TP
DOWNlI)=ODN(TP)
YDWN(I)=HDN(TP)
KSAVE=I
IF(TP.GE.TNAX) GO TO 830
925 820 CONTINUE
C SAVE SAME VALUES FOR U/S END OF POOL
830 CONT INUE
KMAX=KSAVE
DO 040 I=1 ,KSAVE
930 TUP( I )=TDWN(I)
DUP(I)=ODWN(I)
YUP( 1 ,=YDWN(I)
840 CONT I NUE
SAVTlME=TMAX
935 KUP=KMAX
KDWN=KMAX
900 CONTINUE
C SAVE TMAX FOR NEXT POOL
TMAX=SAVTIYE
940 C SAVE KUP FOR USE IN FUNCTION5 ODN AND HDN IN NEXT POOL
KMAX-KUP
C SAVE T. O AND Y FOR MATCHING IN NEXT POOL
MAXPl=KMAX
DO 1000 I=l,MAXPT
945 SAVETUP(I)=TUP(I)
SAVEOUp(I)=OUp(I)
SAVEYUP(I)=YUP(I)
1000 CONTlNUE
C
950 C
C WRITE OUT 1. 0 AND Y AT BOTH ENDS OF POOL
C FOR USE IN GATE OPENING CALCULATIONS
C STRTYPE IS TYPE OF CONTROL STRUCTURE
C DEPTHS ARE IN THE CHANNEL, NOT IN THE GATE SECTION
955 WRITE(l,9070) IPOOL,KDWN
wRITE(i,soeo) (TowN(I),I=~.K~wN)
wRIrEci.9oeo) (oowN(I),I=~,K~wN)
WRITE(l.9080) (YOWN(I),I=l.KOWN)
WRITE(6,9090)KDWN
960 WRITE(6.9lOO)KUP
WRIfE(l.9~10) IPOOL,KUP, STRTYPE(IpDOL)
WRITE(l,9OBo) (TUP(I),I=l,KUP)
WRITE(l.9OBo) (OUP(I),I=l.KUP)
~~;‘~~1,9080’ (VUP(I),I-1,KUP)
965
C FORMATS
9000 FORMAT(llH END REGION,15 )
9010 FORMAT(5H N3= ,110)
9020 FORMAT(lX,4SHTRANSIENT TRAVEL TLME EXCEEDS TIME AVAILABLE /
SUBROUTINE STROKER 74174 OPT= 1 FTN 4.8+420 79/01/09. 14.41.14 PAGE 18

970 2BHTO TMAX. COMPUTATION HALTED.


/lX.6Ht2MAX=,Fl0.3,iOX,5tilMAX=,FlO.3//)
9030’ FORMAT(lX,32HfIME SPAN TOO SHORT TO CALCULATE,
43H AT UPSTREAM BOUNDARY IN REGION 3 FOR POOL ,!J/)
9040 - FDRMAT(lX,OOHNEW TIME IkCREMENT 8,FlO.l)
975 9050 FDRMAT(5H N38= ,110)
9080 FORMAT(lX.17HSfOP 1N 38. tYIN=,F10.3.BH lMAX=,FlO.B//)
9070 FORMAT(lX,311O,lOX)
9000 FORMAl(lX,lOF1D.3)
9090 FORMAT(lX,5HKOWN=.I5//)
990 9100 FORMAT(lX,SHKUP- rIS//)
9110 FORMAT(lX.4110)
9120 FORMAf(lX.2BHTREAtING AS ZERO LENGTH POOL )
END
SUBROUTINE SOLVER 74/74 OPT= 1 FTN 4.6+420 79/01/09. 14.41.14 PAGC 1

1 SUBROUTINE SOLVER
C ROUTINE TO SOLVE FOR X,t,G AND Y Al POINT p
C USING VALUES AT POINT R ON + CHARACTERISTiC
C AND VALUES Al POINT S ON - CHARACTERISTIC
s C METHOD AND NOTATION FOLLOW WYLIE
C
COMMON /SOLV/ XP,YP,fP,GP.XR,yR,TR.QR,XS,YS,TS,QS
C
COMMON /ALL/
10 C PHYSICAL PARAMETERS
. NMANN,SO.GGRAV,PI,BOTlOM(5O),SIOSLOP(5O),RADIUS(5O),BOTGATE(5l
MANN(5O),9OTGRAD(5O),SYFONOY(5O),SYFNLOS(5O),SPEEO(5l),PUMP~
: (51),DEDBANO(51),NGATES,GATECO(ll,5l),OS~NVT(5O),USXNVT(5l),
C WORKING VARIABLES
15 . TMAX,XEND,~INC,NXINC,NPOOLS,IPOOL,TYPPWL.STABEG(SO),
STAEND(SO),Tl,Y1,STUDY,
C CANAL DESCRIPTION
. UPS~A(~O),DWNSTA(SO),NOPOOLS,POOLT~P(~O),STRTYPE(~~),
TMAXI.NUMTRNS,USRES,DSRES
20 c l

REAL NMANN ,MANN


INTEGER TYPPOOL,POOLTYP,STRTYPEIS~U~Y
C

25 CR=Cf YR)
CS=Cf VS)

‘JI
0 30 2~=:1:~;““’

AS=Af VS)
VS=QS/AS
SS=sfVS.RS)
35 VP=(VR+VS)/2.
VPPREV=vp
YP=(VR+YS)/2.
YPPREV=yp
Tpm((v~+c~)*tR-(VS-CS)rtS-XA+XS)/(vR+CR-vS+CS)
40 xpp(vR+CR)*(TP-TR)+XR
SO=BOlGRAD(XpOOL)
AP=A( VP)

C lCOUN:P;:~~:PLOOPING IF SOLUTION DOES NOT CONVERGE


45 ICOUNT=*
10 CONTINUE
RP=R( VP)
APpA( VP)
VP*Gp/AP
50 SP=S(Vp,RP)
CP=C( VP)
TP=(2.*(XS-XR)+TR*(VP+CP+VR+CR)-TS*(VPICP+VS-CS))/
(VR+CR-(VS-CS)+l.*CP)
CS=GGRAV/2.+(1./CR+l./CP)
55 C3=VR+C4*YR-GGRAV/2.*(SR+SP-2.+SO)*(TP-TR)
C2=GGRAV/2.*(l./CS+l./CP)
Cl=VS-C2+YS-GGRAV/2..(SS+SP-?.+SO)*(TP-TS)
SUBROUTINE SOLVER 74/‘14 OPT= 1 FTN 4.6+420 t9/01/09. 14.41.14 PAGE 2

YPPREV=YP
YP=(C3-cl)/(c2+c4)
60 VPPREV=VP
VP=C3-C4.Y P
AP=At VP)
QP=VP*AP
ICOUNl=ICOUNt+l
65 IF(ICOUNT.GT.50) WRIlE(6,9000)
IF(ICOUNT.GT.SO) STOP
IF(At?S(YPPREV-YP).GT.O.OOl) GO TO 10
IF(AOS(VPPREV-VP).GT.O.OOl) GO TO 10
XP=XR+((VP+VR)/2.+(CR+CP)/2.)*(TP-TR)
70 C CHECK THAT FLOW IS NOT SUPERCRITICAL
CELRITY=C(YP)
AREA=A(YP)
IF(AOS(CELRITY+AREA/OP).LE.1.) WRITE(6,9010) XP,YP,TP,QP
IF(ABS(CELRITY*AREA/oP).LE.l.) STOP
75 RETURN
C FORMATS
9000 FORMAT(lX.14HSTOP IN SOLVER)
9010 FORMAT(lX.22HSUPERCRITICAL FLOW AT
.. /5X,3HX= ,FlO.l/SX,3HY* ,FlO.l/5X,
80 . MT= .FlO.l/5X,3HO= ,FlO.l/)
FUNCTION C 74174 OPT= 1 FtN 4.(1+428 PAQE 1

FUNCTION C(Y)
C CELERITV
C
COMMON /ALL/
C PHYSICAL PARAMETERS
. NMANN.SO.GGRAV,PI,BOfTOM(W),SIOSLGP(SO),RAO1US(SO),BOTGATE(Sl),
. MANN(5O),BOTGRAD(50),SYFONOV(5O),SVFNLOS(5O),SPEEO(5l),PUMP~
(51).OEDBANO(S1),NGATES,GATECO(ll,5l),OSINVT(5O),USINVT(5l),
C WORliING VARIABLES
10 . TMAX.XENO,TINC.NXINC,NPWLS,IPOOL,TVPPWL,STABEG(5O),
STAENO(50) ,Tl ,Vl ,STUOV.
C CANiL OESCRIPTION
. UPSTA(5O),OWNSTA(5O),NOPWLS,POOLTVP(50),STRTVPE(Sl),
TMAXI.NUMTRNS,USRES,OSRES
lb c *
REAL NMANN ,MANN
INTEGER TYPPOOL,POOLTiP,STRTYPE,STUDV

:
20 C A=AREA
C TY-TOP WIDTH
C.SQRT(GGRAV*A(V)/TW(Y))
RETURN

FUNCTION R 74/74 OPT= 1 FtN 4.9+429 T@/01/09. 14.41.14

FUNCTION R(Y)
C HYDRAULIC RADIUS
C A-AREA
C PWETTED PERIMETER
yl;~““’
END
FUNCTION A 74/74 OPT= 1 FTN 4.0+428 79/01/09. 14.41.14 PAGE 1

FUNCTION A(Y)
C CROSS SECTIONAL AREA
C POOL IS TRAPEZOIDAL CHANNEL - tYPPOOL=l
C OR HORSESHOE TUNNEL - TYPPOOL=Z
C OR CIRCULAR TUNNEL - TYPPOOL=O
C
C BOTTOM WIDTH AND SIDESLOPE OF TRAPEZOIOAL CHANNEL
C ARE SPECIFIED IN BOTtOM(POOL) AND SIOSLOP(POOL)
C
10 COMMON /ALL/
C PHYSICAL PARAMETERS
. NMANN,SO,GGRAV,pI,BOTTOM(5O),SIOSLOp(5O),RAOIUS(5O),BOTGATE(5l),
. MANN(~O),BOTGRAO(~O),SYFONOY(~O),SYFNLOS(~O),SPEEO(~I),PUMPO~
(~~~.OEO~ANO(~~),NGATES,GATECO(~~~S~),DSINVT(~O),USINVT(~~),
IS C WORKING VARIABLES
. TMAX.XENO,TINC,NX1NC,NPOOLS,IPOOLrTYpPWL,STABEG(50),
STAEND(5O),Tl.Y1.STUOY,
C CANAL DESCRIPTION
. . UPSTA(50).0WNSTA(5O),NOPOOLS,POOLTYP(5O),STRTYPE(Sl),
20 TMAXI.NUMTRNS,USRES,OSRES
c ’
REAL NMANN .MANN
INTEGER TYPPOOL,POOLTYP,STRTYPE,STUDY

25 :
IF(Y.LE.O.0) YRITE(B,9001) IPOOLJ
IF(Y.LE.O.0) STOP
IF(TYPPOOL.EO.1) GO TO 50
IF(TYPPOOL.EfJ.3) GO TO 60
30 C IN HORSESHOE TUNNEL SECTION
C RIO IS INVERT TO CREST DISTANCE
RAO~RAOIUS(IPOOL)
C CHECK FOR LOWEST ARC
Yl=RAD/4.*(3.-SGRT(7.))
35 IF(Y.GT.Yl) GO TO 10
A=(Y-RAO)*SORT(RAO**Z-(Y-RAO)**Z)
+RAO**Z*ASIN((Y-RAO)/RAO)+PI*RAO**2/2.
RETURN
C CHECK FOR AT SPRINGLINE
40 10 IF(ABS(Y-RA0/2.).GT.0.00001) GO TO 20
YZ-RAD/Z.
As(Yl-RAO)*SORT(RAO+*2-(Yl-RAO)**Z)
. +RAO*+2+ASIN((Yl-RAO)/RAO)+PI*RAO**2/2.
. -(Yl-YZ)*SQRT(RAO**Z-(Yl-Y2)**2)
45 -RAO**Z*ASIN((Yl-YZ)/RAO)-RAO*(Y-Yl)
’ RETURN
C
C CHECK FOR BELOW SPRINGLINE
20 IF(Y.GT.RAO/Z.) GO TO 30
50 YZ=RAO/Z.
A*(Yl-RAO)*SORT(RAO+*2-(Yl-RAO)**Z)
+RAO*+Z*ASIN((Yl-RAO)/RAO)+pI*RAO**2/2.
: +(Y-Y2)+SORT(RAO+*2-(Y-Y2)**2)
. +RAO+*Z*ASIN((Y-YZ)/RAO)
55 . -(Yl-YZ)*SQRT(RAO**Z-(Yl-Y2)**2)
. R;~~~;*2’ASIN((Yl-Y2)/RAO)-RAO*(Y-Yl)
FUNCTION A 74/74 OPT= 1 FTN 4.6+426 79/01/09. 14.4t.14 PAGE 2

C CHECK FOR FILLED TUNNEL


30 IF(Y.GT.RAD) GO TO 40
60 Yl*RA0/2.
A=(Yl-RAD)*SORT(RAD+*2-(Yl-RAD)+,2)
. +RADf*2+ASIN((Yl-RAD)/RAD)+PI*RAD**2/2.
. -(Yl-Y2)*SORT(RAD*+2-(Yl-Y2)**2)
. -RAD*+2+ASIN((Yl-Y2)/RAD)-RAD*(Y2-Yt)
65 +(Y-Y2)+SORT((RAO/2.)**2-(Y-Y2)**2)
: R;;;R;/2.)**2*ASIN((Y-Y2)/(RAD/2.))

40 WRITE(6,9000) Y,IPOOL
STOP
70 C TRAPEZOIDAL CHANNEL SECTION
50 CONTINUE
~;~~~~TOH(IPOOL)+S1DSLOP(!POOL)*Y)*Y

C CIRCULAR TUNNEL SECTION


75 60 CONTINUE
RAD=RADIUS(IPOOL)
IF(A6S(Y-RAD).GT.O.OOOOl) GO TO 70
C HALF-FULL
A=PI*RAD*+2/2.
60 RETURN
70 IF(Y.GT.RAD) GO TO 60
CJl C LOilER HALF OF SECTION
ul XFRAD-Y
THETA=ACOS(X/RAD)
65 AstHETA*RAD*+2-X+RAD*SIN(THETA)
RETURN
C UPPER HALF OF SECTION
60 IF(Y.GT.Z.+RAD) GO TO 90
X=Y-RAD
90 THETA=ACOS(X/RAD)
A=(PI-THETA)*RAD*+2+X*RAD*SIN(THETA)
RETURN
90 WRITE(6,9000) Y,IPOOL
STOP
96 C FORMAT
9000 FORMAT(lX.36HDEPTH GREATER THAN TUNNEL HEIGHT At
3H Y=,F10.2.9H IN POOL ,151
9001’ FORMAT(lX,l3HDEPTH IN POOL,I5,2W -,FlO.l/
, lX,POHSTOP IN AREA ROUTINE )
100 END
FUNCllON TW 74/74 OPT= 1 FTN 4.6+420 79/01/m. 14.41.14 PAM 1

1 FUNCTION TIM(Y)
c TOP UIOTH
C
COMMON /ALL/
5 C PHYSICAL PARAMETERS
. NMANN,SO,GGRAV,PI,BOTlOM(50),SIOSLOp(5~~,RADIuS(5O),~OTGATE(5~~,
MANN(~~),BOTGRAO(~~),SYFONDY(~O),SYFNLOS(~O~,SPEED(~~~,~U~~~
: (5l),OEOBAND(51),NGATES,GAfECO(ll,5l),DSINVl(5O),US1NVT(5~~,
C WORKING VARIABLES
10 . fMAX,XENO,TINC,NXINC,NPOOLS,IPOOL,TYpp#)L,STABEG(5O),
STAEN0(50),Tl,Yl,STUDY,
C CANAL DESCRIPTION
. UPSTA(~~),OWNSTA(~~),NOPOOLS,POOLTYP(~O~,STRTY~E(~~~,
. TMAXI.NUMTRNS,USRES,DSRES
15 C
REAL NMANN ,MANN
INTEGER TYPP0OL,POOLTYP,STRTYPE,STUDY

20 POOL IS TRAPEZOIDAL CHANNEL - TYPpOOL=l


OR HORSESHOE TUNNEL - TYPPOOL=2
OR CIRCULAR TUNNEL - TYPPOOL.=3

BOTTOM WIDTH AND SIDESLOPE OF TRAPEZOIDAL CHANNEL


25 ARE SPECIFIED Ii BOTTOM AND SIDSLOP
IF(TVPPOOL.EG.1) GO TO 50
_... ~--
IFIlVPPOOL.EO.31 GO TO 60
IN HORSESHOE TUNNEL SECTION
RID IS XNVERT TO CREST DISTANCE
30 RADtRAOIUS(IPOOL1
CHECK FOR LOWEST ARC
Vl=RAO/4.+13.-SORT(7.11
IF(Y.GT.Ylj GO TO i0 --
TWs2.*SGRT(RAD*+2-(RAD-Y)**2)
35 RETURN
C ,$lEcK FOR AT SPRINGLINE
IF(ABS(Y-RAD/2.).GT.O.GOOOl) GO TO 20
TW=RAO
RETURN
“0 C CHECK FOR BELOW SPRINGLINE
20 Y2+RAO/2.
IF(Y.GT.Y2) GO TO 30
THETA=ASIN((Y2-Y)/RAo)
TW=RAO*(2.*COS(THETA)-1.)
45 RETURN
C CHECK FOR FILLED TUNNEL
30 IF(Y.GT.RAD) GO TO 40
THETA.ASIN((Y-Y21/(RAD/2.1)
TYIRAO*COS(THETA~
50 RETURN
40 WRITE(6,900G) Y,IpOOL
STOP
C TRAPELOIOAL CHANNEL
50 CONTINUE
55 TW~BOllOM( :p00L)+2.+SIDSL0p( 1p00L)+Y
RETURN
C IN CIRCULAR TUNNEL
FUNCTION TW 74/74 OPT= 1 FtN 4.6+428 79/01/09. 14.41.14 PAGE 2

60 CONTINUE
RAO=RADIUS(IPOOL)
60 IF(ABS(Y-RAD).Gt.O.OOOOl) GO TO 70
c HALF-FULL
tWr2. l RAD
RETURN
IF(Y.GT.RAD) GO TO 80
65 C’IhER HALF OF SECTION
X=RAD-Y
lHETA=ACOS(X/RAD)
fW=2.*RAD*SIN(THEfA)
RETURN
70 C UPPER HALF OF SECTION
80 IF(Y.Gf.2.*RAD) GO TO SO
X=Y-RID
THETA=ACOS(X/RAD)
TW=2.+RAD*SIN(THETA)
7s RETURN
90 WRITE(6.9000) Y,IPOOL
STOP
C FORMAT
9000 FORMAT(lX.3OHDEPTH GREATER THAN TUNNEL HEIGHT At
00 3H Y=.F10.2.9H IN POOL ,151
END
FUNCTION S 74/14 OPT= 1 FTN 4.6+429 7S/Ol/OS. 14.41.14 PAGE 1

FUNCTION S(V,R)
C ENERGY GRADE
C
COMMON /ALL/
C PHYSICAL PARAMETERS
. NMANN.SO,GGRAV,PI,~O~TOM(~O),SIDSLOP(~O),RADIUS(SO),~O~GATE(~~),
. MANN(SO),BOTGRAD(5O),SYFONDY(5O),SYFNLOS(5O),SPEED(5l),PUMP~
(S~,.OEDSAND(S~),NGATES,GATECO(~~,S~),OSINVT(SO),USINVT(~~),
C YORlilNG VAAtABLES
10 . TMAX.XEND,TINC,NX~NC,NPOOLS,IPOOL,TYPPWL,STA~EG(SO),
STAEND(50),Tl,Yl,STUDY,
C CANAL DESCRIPTION
. UPSlA(5O~~DWNSTA(SO),NOPOOLS,POOLTYP(SO),StRTYPE(Sl),
TMAXI.NUMTRNS,USRES,DSRES
15 c -
REAL NMANN,MANN
INTEGER TYPPOOL,POOLtYP,STRTYPE,STUDY

:
20 C V IS VELOCITY
C R IS HYDRAULIC RADIUS
C GET MANNINGS N FOR THIS POOL
NMANN=MANN(IPOOL)
S-(V*NMANN/l.49)**2/R**(4./3.)
25 RETURN
END

n
x)
FUNCTION P 74/74 OPT= 1 FTN 4.6+423 79/oc/o9. 14.41.14 PAM 1

FUNCTION P(Y)
C WETTED PERIMETER
C
COMMON /ALL/
C PHYSICAL PARAMETERS
. NMANN.SO.GGRAV,pI,3OTTOM(5O),SIOSLOp(50),RAOIUS(5O),6OTGATE(5l),
MANt~(50).3OlGRAD(5O),SYFONDY(5O),SYFNLOS(5O),SPEED(5l),PUMP~
: (51~.DEDBAND(51),NGATES,GATECD(ll,5l),OSINVT(3O),USINVT(5l),
C WORKING VARIABLES ’
10 . TMAX.XEND,lINC,NXINC,NPOOLS,IPOOL,TYPPOOL,STA3EG(3O),
STAEND(SO),Tl,Yl,STUDY,
C CAN;L DESCRIPTION
. UP5~A(5O),OWNSTA(5O),NOpOOLS,POOLTYP(50),5TRTYPE(5~),
TMAXI,NUMTRNS,USRES,DSRES
15 c *
REAL NMANN ,MANN
INTEGER TYPPOOL,POOLTYP,STRTYPE,STUDY

: POOL IS TRAPEZOIDAL CHANNEL - TYPPOOL=l


20 C OR HORSESHOE TUNNEL - TYPPOOL=2
C OR CIRCULAR TUNNEL - TYPPOOL=B

: BOTTOM WIDTH AND SIDESLDPE OF TRAPEZOIDAL CHANNEL


C ARE SPECIFIED IN 3OTTDM(pOOL) AND SIDSLOp(pOOL)
23 IF(TYPpOOL.EO.1) GO TO $0
IF(TYPPOOL.EO.2) GO TO 60
C IN TUNNiL SECTION
C RAD IS INVERT TO CREST DISTANCE
RAO=RADIUS(1pODL)
30 C CHECK FOR LOWEST ARC
Yl=RAD/4.*(3.-SQRT(7.))
IF(Y.GT.Yl) GO TO 10
P=2.*RAD+ACOS((RAD-Y)/RAO)
RETURN
35 10 IF(ABS(Y-RAD/2.).GT.0’.00001) GO TO 20
Y2=RAD/2.
PI~.‘RAD*(ACOS((RAD-Y~)/RAD)+PI/~.-ACOS((Y~-Y~)/RA~))
RETURN
C CHECK FOR BELOW SPRINGLINE
40 20 Y28RAD/2.
IF(Y.GT.Y2) GO TO 30
Pm2.+RAD*(ACOS((RAD-Yl)/RAO)+ACOS((Y2-Y)/RAD)
. R;:~~~W2-Yl)/tW)

45 C CHECK FOR FILLED TUNNEL


30 IF(Y.GT.RAD) GO TO 40
P~~.~RAD+(ACOS((RAD-Y~)/RAD)+PI/~.-ACO~((Y~-Y~)/RAD))
. +RAD*ASIy((Y-Yl)/(RAD/2.))
RETURN
50 40 YRITE(6,9000)Y,1pOOL
STOP
C TRAPEZOIDAL CHANNEL
50 CONTINUE
P=3OffOM(IPOOL)+Y*2.*SORTfSIDSLOP(IPOOL)**2+l.)
55 RETURN
C IN CIRCULAR TUNNEL
60 CONTINUE
FUNCTION P 74/74 OPT= 1 FtN 4.6+426 79/01/09. 14.41.14 PAGE 2

RA bRADIUS( IPOOLJ
A6S(Y-RAOJ.Gt.O.OOOOlJ GO TO 70
60 C HALF-:: .L
P= ‘l*RAD
RE ‘URN
70 Y.GT.RADJ GO TO 60
c LOSIER*t ILF OF SECTION
65 IAD-Y
:;, LTA=ACOS(X/RAD)
P= I. l RAD*THETA
‘URN
c UPPEIE LLF OF SECTION -
70 60 tF ‘Y.GT.Z.+RAOJ GO TO 90
~-RID
t;; :TA=ACOS(X/RADJ
Pm- 1. *RAWI PI-THETA)
RETURN .
76 90 YRITE(6,9000J Y,IPOOL
STOP
C FORMAT
9000 FORMAT(lX.39HDEPTH GREATER THAN TUNNEL HEIGHT AT
3H Y=.FlO.2,SH IN POOL ,lSJ
66 ’ END
SUBROUTINE BNORV 74/74 OPT= 1 FtN 4.6+429 79/01/09d 14.41.14 PAQE 1

SUBROUTINE 9NORV(IUPROWN)
C UPSTREAM OR DOWNSTREAM BOUNOARV
C IUPRDWN ZERO FOR UPSTREAM BOUNDARY
C 1UPRDWN ONE FOR DOWNSTREAM BOUNDARY
C 0 IS SPECIFIED IN FUNCTIONS GUPST AN0 QDN
C NOTATION FOLLOWS WYLIE
C
COMMON /SOLV/ XP,VP,TP.OP,XR,YR,fR,GR,XS,VS,TS,QS
C.
10 COMMON /ALL/
C PHYSICAL PARAMETERS
. NMANN.SO,GGRAV,PI,9OTtOM(5O),SIDSLOP(5O),RADIUS(5O),9OTGATE(5l),
. MANN(50).9OTGRAD(5O),SVFONDV(50),SVFNLOS(5O),SPEED(5~),PUMPW
. _ I51).DED9AND(51),NGATES,GATEC0(11,51),DS~NVT(50),US1NVT(51),
15 C WORKiN ViRIABLES -
. tMAx,XEND,tINC,NXINC,NPOOlS,IPOOl,tyPP~l,STA9EG(59~,
StAEND(SO),Tl,Yl,STUDY,
C CAN;1 DESCRIPTION
. UPStA(5O),DWNSTA(50),NOPbOLS,POOltyP(5O),STRTyPE(5~),
20 . TMAXI.NUMTRNS,USRES,DSRES

REAL NMANN ,MANN


INTEGER TVPPOOl,POOLTYP,STRTYPE,STUDY
C
25 C
IF(IUPRDWN.EQ.0) GO TO 20
C DOWNSTREAM BOUNDARY
C INTERSECT C+ WITH BOUNDARY
CR=C( VR)
30 RR=R( VR)
AR=A(YR)
VRsQR/AR
SR=SlVR,RR)
V=VR
35 b;;;(;P-XR)/(VR+CR)
I
Y*VR
X=XEND
SO=9OtGRAO(1POOl)
40 AP=A(V)

%:b”
DiREV=CJ
10 CONT I NUE
45 VPREV’Y
RP=R( V)
SP=S(V.RP)
CP=C(Y)
Y=YR+(~.*(VR-V)/GGRAV-(SR+SP-2.*SO)*(T-TR))
50 . /(l ./CR+1 ./CP)
RP=R(Y)
SP=S(V.RP)
CP=C(V)
AP=i(vj*’
55 V=VR-GGRAV/2.*(SR+SP-2.*SO)*(T-TR)
-GGRAV/l.*(l./CR+l./CP)*(Y-YR)
’ TPREV=t
SUBROUTINE BNDRY 74/14 OPT=1 FTN 4.6+42C 79/01/09b 14.41.14 PAGE L

T=2.*(XP-XR)/(V+VR+Ca+cP)+TR
OPREV=O
60
:‘:K~T’
I;(A%(OPREV-O).GT.O.OOl) Go TO 10
IF(ABS(TpREV-t).Gt.O.OOlI GO TO 10
IF(AES(YPREV-Yj.GT.O.OOlj Go TO 16
85 VP=V
TP=T
YP=Y
QP-0 ~~
c CORRECT FOR CHANGE IN CHANNEL CROSS SECTION AT ~0 GATE o/s BOUNDARY
70 IF(IPOOL.EO.l) RETURN
C YBEFORE IS D/S Y CALCULATED FROM INITIAL CONDITIONS
C AND CONTINUITY OF 0
C YAFTER IS D/S Y CALCULATED FROM ENERGY BALANCE
c YBEFoRE AND YAFTER SHOULD BE NEARLY EOUAL FOR
75 C A VALID SOLUTION ACROSS THE CHANGE IN SECTION
IF(SlRTYPE(IPOOL-l).EQ.4) WRITE(6.9000) YP
IF(STRlYPE(IPOOL-l).E0.4) YP=HDN(T)
IF(slRTYPE(rPOOL-l).EQ.4) YRITE(6,9010) yP
IF(SlRlYPE(1POOL-l).EO.5) YRITE(6.9000) YP
80 IF(SlRTYPE(IPOOL-l).EQ.S) YP=HDN(T)
IF(sTRTYPE(IPOOL-0.EO.B) WRITE(B.9010) VP
RETURN
C UPSTREAM BOUNDARY
C INTERSECT C- WITH BWNDARY
85 20 CONTINUE
CS=Cl VSI
RS=Rf YS)
AS=A,cYS)
vS=OS/AS
90 ss=s(vS,RS)
V=VS
w&(;P-xs)/(vs-cs)
i

IOTGRAD ( I POOL)
I(V)
IpST(T)
‘AP
100 :v=o
30 ‘INUE
:v=y
I(Y)
;(V.RP)
105 CP=C(Y)
Y=YS-(2.*(VS-V)/GGRAV-(SS+sp-2..s9)*(T-Ts))
/(I ./cs+1 ./CP)
‘=R(V)
‘=S(V,RP)
110 W(Y)
‘=A(Y)
IVs-GGRAV/2.+(SS+SP-2.*sO)*(T-Ts)
~~RAV/~.+(~./CS+I./CP)~(~-YS)
‘REV=T
SUBROUTINE BNDRY 74/74 OPT= 1 FTN 4.6+429 79/0t/09. 14.41.14 PAGE 3

115 1=2.*(xP-xs)/(v+vs-c+cP)+Ts
QPREV-0
O=QUPST(T)
VrQ/AP
IF(ABS(QPREV-O).Gl.O.OOl) GO TO 30
120 IF(ABS(lPREV-T).Gt.O.OOl) GO TO 30
XF(ABS(YPREV-Y).GT.O.OOl) GO TO 30
VP=V
TP=T
YP=Y
125 OP.0
RETURN
C FORMATS
9000 FORMAT(lX.9HY BEFORE=,F10.3)
9010 FORMAl(lX.OHY AFTER=,FlO.B/)
130 END
FUNCTION HDN 14/74 OPT= 1 FfN 4.0+426 7P~01/09. 14.41.14 PAN i

1 FUNCTION HDN(T)
C DOiiNSfREAM BOUNDARY DEPTH
C T IS TIME IN SEC
C
5 COMMON /ALL/
C PHYSICAL PARAMETERS
NMANN.SO,GGRAV,PI,6OTTOM(SO),SIDSLDP(5O),RAD~US(5O),6OTQATE(5l),
: MANN(50).8DTGRAO(50),SVFONOV(5O),SVFNLOS(5O),5PEED(5l),PUMP~
(S1~,DED6AND(51).NGATES,GATECO(ll,5l),DS~NVT(SO),USlNVT(5l),
10 C WORKING VARIABLES
. TMAX,XENO,TINC,NXINC,NPDOLS,IPOOL,TVP~L,STA6EG(SO),
STAEND(50).Tl,Vl,STUDV.
C CANAL DESCRIPTION
. UPSTA(5O),DWNSTA(5O),NOPOOLS,POOLTVP(50),5TRTYPE(5l),
15 TMAXI,NUMTRNS.USRES,DSRES
c ’
-REAL NMANN ,MANN
INTEGER TVPPDOL,PDDLTVP,STRTVPE,5TUDV

20 :
C SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINE5
C TO KEEP MEMORY WITHXN SOME SORT OF REASONASLE BOUND
COMMON /SCRATCH/
, Dum:.1v 1 ( 8570)
25 . ,H(46)
. .DUKMV2(5666)
. , KMAX
. .SAVETUP(500),5AVEVUP(56O),SAVEDUP(5DD)

30 : MAKE2NERGV BALANCE DEPTH CHANGE AT ND GATE CHANGE IN CROSS-SECTIDN


IF(IPOOL.ED.l) GO TO ID
)F(SlRTVPE(IPOOL-l).EG.4) GO TO 40
IF(SlRTVPE(IPOOL-l).EG.S) GO TO 40
10 CONTINUE
35 TIMEHR=T/3600.
ITIME=TIMEHR
ITIME=ITIME+l t
IF(ITIME.LT.l) ITIME=l
IF(lTIhlE.GT.48) ITIME=
40 IF(IlIME.GT.l)-GO TO 20
C Tl.vl IS LAST COMPUTATION POINT ON D/S BOUNDARV
C IN REGION 2
C MATCH WITH Tl,Vl IN 1ST HR
C THIS ASSUMES Tl.LT. 1 HR
45 C CHECK WAVE TRAVEL TIME IF IN DOUBT
IF(Tl.GE.3600.) WRITE(6.9000) Tl
DELV=H(21-Vl
DELT=3600.-Tl
DT=T-11
50 DV=DELV*DT/DELT
HDNFVl+DV
RETURN
C LINEARLY JOIN H(ITIYE) AND H(ITIYE+l)
20 IF(ITIME.EO.46) GO TO 30
55 DELV=ti(ITIME+l)-H(ITIME)
DELT=3600.
I=ITIME-1
FUNCTION HDN 74/74 OPT= 1 FTN 4.6+426 76/01/09. 14.41.14 PAGE 2

TBEG= I
TBEG=l8EG*3600.
80 DT=T-TEEG
DY=DELY+DT/DELT
HDN=HfITIME)+DY
RETURN
C HOLD DEPTH CONSTANT AFTER 47tH HOUR
66 30 HDN=H( ITIME)
RETURN
C MODIFIED NEWTONS METHOD FOR ENERGY BALANCE
C4~P+VUP**2/2*G~YDOWN+VDOWN+*2/2*G
CONTINUE
70 C SET DOWNSTREAM VALUES
C FIND DEPTH tN DOWNSTREAM SECTION
DO 50 I=2,KMAX
J=I
IF(T.LE.SAVETUP(I)) GG TO 60
76 60 CONTINUE
c T GREATER THAN TMAX - USE FINAL DEPTH
YD=SAVEYUP(KMAX)
GO TO 70
C INTERPOLATE FOR Y
60 60 CONTINUE
DEl.Y=SAVEYUP(d)-SAVEYUP(J-l)
DELT.SAVETUP(sJ)-SAVETUP(J-l)
DT=T-SAVETUP(J-1)
DY+OELY*DT/DELT
66 YD=SAVEYUP(J-l)+DY
C GET AREA AND 0 IN OOWNSTREAM SECTION
70 __._ --
CONTINUE
C GDN RETURNS 0 UPSTREAM OF TURNOUT
OO=ODNIT1
_-- 7

SO C GET TURNOUT‘Os
CALL tRNOUT(DELG,T).
ODSFOO-DEL0
OUS=OO
C SET IPOOL FOR AREA ROUTINE
96 IPoOL=lPOOL-1
TYPPOOL=POOLTYP( IPOOL)
C IF SYPHON, CORRECT DEPTH FOR SYPHON LOSSES AND INVERT
IF(STRTYPE(IPOOL-l).NE.4) GO TO 76
C SYPHON EFFECTS ONLY - NO CHANGE IN CROSS-SECTION OF CHANNEL
100 C CHANGE IN CROSS-SECTION, IF ANY, IS TAKEN INTO ACCOUNT LATER
Y=YD
ABELOW=A(YD)
YSYPHON=SYFONDY(IPWL)
SYPHON=SYFNLOS(IPOOL)
106 ClrOD~**2/(2.+00RAV+ABELOY+*2)+YD-YSYP~+SY~ON*ODS**2
RAO=RADIUS(IPOOL)
71 CONTINUE
IF(TYPPOOL.EG.2.AND.Y.GE.RAD) YFRAD-0.1
IF(TYFPOOL.EO.3.AND.Y.GE.2.*RAD) Y=2.+RAD-0.1
110 IF(Y.LT.O.0) Y-O.1
F=Y+ODS+*2/(2.+GGRAV+A(Y)*+2)-Cl
FPRIME=l.-GDS++P/(GGRAV*A(Y)**3)+TW(Y)
YPREV=Y
Y=Y-F/FPRIHE
FUNCTION HON 74/74 OPT=1 FtN 4.6+428 79/01/09. 14.41.14 PAGE 3

115 IF(ASS(Y-YPREV).Gl.O.OOl) GO TO 71
YO=Y
C CALCULATE EFFECT DUE TO CHANGE XN CROSS-SECTION
75 CONTINUE
AD=A( YD)
120 C RESET IPODL FOR U/S POOL
IPDOL=IPDOL+l
lYPPOOL=PODLfYP(1P00L)
YU=YD
C TCOUNT STOPS PROGRAM IF NO CONVERGING
125 ICOUNf=l
RAD=RADIUS(IPOOL)
80 CONTINUE
YPR iv=vu
F=- ‘~~YD+(l./(2.*GGRAV)+(ODS++2/AD**2-OUS++2/A(YU).~2))
130 FPR =-l.+QUS+*2*TW(YU)/(GGRAV+A(YU)++3)
YU+ ‘U-F/FPRIME
IF( ‘U.LT~O.1) YU=O.l
tF( ‘VPPOOL.NE.2) GO TO 90
lF( ‘U.GE.RAD) YU=RAD-0.1
135 90 CON ‘INUE
IF( ‘VPPOOL.NE.3) GO TO 100
IF( UM;.2.*RAD) YLbl.*RAD-0.1
100 CON
ICO INT=ICOUNt+l
140 CDUNf.GT.20) WRITE(6,9010)
::I CDUNT.GT.20) STOP
IF{ rBS(YU-YPREV).GT.O.OOl) GO TO 80
HDN IVU
RETURN
145 C FORMATS
9000 FORMAT(lX.3HTl=,F10.3,7H tN HDN )
9010 FDRMAT(lX,llHSTOP XN HDN )
FUNCTION GDN 74/74 OPT= 1 FtN 4.6+426 rs/ol/oo. 14.41.14 PAQE 1

FUNCTION GON(t)
C DDiINSTREAM FLOW AT EACH POOL
C SPECIFY G AT POOL I IN DD
C Q CONTINUITY IS MAINTAINED FOR ALL OTHER POOLS
C T IS TIME IN SEC
C
COMMON /ALL/
C PHYSICAL PARAMETERS
NMANN.SO,GGRAV,PI,9DTTDM(50),SIDSLOP(SO),RADIUS(SO),BOTGATE(bl),
10 : MANN(50).6OTGRAD(SO),SYFONDY(SO),SYFNLDS(SO),SPEED(5l),PUMP~
. (51~.DED6AND(51).NGATES.GATECD(11r5l),DSINVT(50).USINVT(5l).
C WORKING VARIABLES
. TMAX.XENDtTINC,NXINC,NPDOLS,IPOOLrTYPPWL,STA5EG(~9),
. STAEND(5O),Tl,Yl,STUDY,
15 C CANAL DESCRIPTION
. UPSTA(5O),DWNSTA(50),NDPODLS,PDOLlYP(5O),STRTYPE(5l),
. TMAXI.NUMTRNS,USRES,DSRES
C
REAL NMANN ,MANN
20 INTEGER TYPPOOL,POOLTYP,STRTYPE,STUDY

:
C SCRATCH AREA TO BE OVERLAID iN THE VARIOUS SUBROUTINES
C TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
25 COMMON /SCRATCH/
DUMMVl(6616),OD(46),DUMMY2(5620),KMAX,
: SAVETUP(5OO),SAVEYUP(5OO),SAVEOUP(5O9)
C

30 : USE
SAVETUP. SAVEGUP’FOR SECOND ANU FOLLOWING POOLS
IF(IPOOL.GE.2) GO TO 10
C 0 STEPS Al EACH HOUR
lIMEHR=T/3600.
ITIHE=TI.M5HR
35 ITIME=ITIME+l
IF(ITIME.LT.l) lTIME=l
IF(ITIME.GT.46) ITIME-
DDN=OO( ITIME)
C CHECK FOR TURNOUT
40 CALL fRNOUT(DELG,T)
ODN=OON+DE LO
RETURN
C 0 FOR SECOND AND FOLLDWINQ POOLS IN EACH REACH
C LINEAR INTERPOLATION OF UPSTREAM 0 FROM PREVIOUS POOL
45 10 CONTINUE
LIYFKMAX-1
DO 20 I=1 ,LIk
K=I
IF(T.GE.SAVETUP(I).AND.T.LT.SAVETUP(I+l)) GO TO 30
50 20 CDNT INUE
IF(T.EG.SAVETUP(KMAX)) GO To 30
C LET T RUN BEYOND TMAX FOR MDST UPSTREAM POOL FOR LAST POINT
IF((T.GT.SAVETUP(K)).AND.(IPDOL.EG.NPOOLS)) GG TO 30
YRITE(6,9000) T,IPWL
55 STOP
30 CONTINUE
DT-T-SAVETUP(K)
FUNCTION QDN 74/74 OPTS 1 FTN 4.6+428 79/01/09. 14.44.14. PAQE 2

DELt=SAVElUP(K+l)-SAVETUP(K)
RAT10=01/DEL1
60 QDN=SAVEOUP(K)+RAtIO+(SIVEQUP(K+l I-SAVEQUP ‘(K))
C CHECK FOR TURNOUT
CALL tRNOUT(DELQ,T)
QDNsQON+DE LO
RETURN
65 C FORMAT
BOO0 FORMAT(lX,16HQ(T) NOT FOUND AT ,f '1 ,0.3*
13H SEC FOR POOL,IS/)
END
FUNCTION OUPST 74/74 OPT= 1 FTN 4.6+429 79/01/09, 14.41.14 PAGE 1

1 FUNCTION OUPSt(t)
C UPSTREAM FLOW
C SPECIFY 0 -IN QU FOR MOST UPSTREAM POOL
; t IS TIME IN SEC
5
COMMON /ALL/
C PHYSICAL PARAMETERS
NMANN,SO.GGRAV,PI,BOttOM(5O),SIDSLOP(~O),RAOIUS(5O),BOtGAtE(~~),
: MANN(SO),~OTGRAD(SO),SYFONDY(SO),SYFNLOS(5O),SPEED(5l),PUMPW
(Sl!.DEDBAND(5l),NGAtES.OAfECO(ll,5l),DSINVt(5O),USINVt(5l),
C WORKING VARIABLES
. tMAX.XEND,tINC,NXINC,NPOOLS,IPOOL.TYPPWL,StABEG(5G),
StAEND(5O),tl,Yl,StUDY,
C CANAL DESCRIPTION
. UPSTA(5O),DWNStA(SO),NOPOOLS,POOLTYP(5O),StRtYPE(5~),
tMAXI,NUMtRNS,USRES,OSRES
c *
REAL NMANN ,MANN
INTEGER TYPPOOL,POOL~YP,S~R~YPE,S~UDY
C

E SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES


C TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
COMMON /SCRATCH/
25 DWM~1(8666),QU(4E),DUMMV2(7073)
C
C 0 STEPS AT EACH HOUR
tIMEHR=t/3600.
ITlME=TIMEHR
30 ItIME=ItIME+l
IF(ITIME.Lt.1) ltIME=1
IF(ITIME.Gt.4B) ItIME
@JiSl=OU(ItIME)
RETURN
35 END
SUBROUTINE TRNOUT 74/74 OPT=1 FTN 4.6+429 79/01/09. 14.41.14 PAN 1

1 SUEROUTINE TRNOUT(DElO,f)
C ROUTINE TO GET TURNOUT 0 FOR USE IN FUNCTION OON

Ii COMMON /All/
C PHYSICAL PARAMETERS
NMANN.SO,GGRAV,PX.BOflW(SO),SIOSlOP(50),RAD1US(5O),BOtQATE(bl),
: MANN(5O),BOtGRAD(5O),SYFONDY(5O),SYFNlOS(5O),SPEED(5~),PU~P~
(Sl,,OEDBANO(Sl),NGATES,GATECO(ll,5l),OSINVT(5O),USINVT(Sl),
10 C WORliING VARIABLES
. TMAX.XEND,TINC,NX1NC,NPOOlS,IPODl,TYPPWl,STABEG(50),
STAENO(50).Tl,Yl,STUOY,
C CAN;1 DESCRIPTION
. UPSTA(5O),DWNSTA(5O),NOPOOlS,POOlTYP(5O),STRTYPE(Sl),
15 TMAXI.NUMTRNS,USRES,DSRES
c *
REAL NMANN ,MANN
INTEGER TYPPOOl,POOlTYP,STRTYPE,STUDY
C
20
; SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES
C TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
COhlMON /SCRATCH/
0UMUY1(9714),OT(49,50),0UMMY2(4673)
25 C 1N:;IAlIZE TURNOUT 0
DElO=O.O

: NUMTRNS IS NUMBER OF TURNOUTS 1N THIS REACH


IF(NUMTRNS.EO.0) RETURN
30 TIYEHR=T/3600.
ITIME=TIMEHR
ITIME=ITIML+l
IF(tTIME.lT.l) ITIME=
IF(ITIME.GT.49) ITIME=
35 DElO=QT(ITIME,IPDOl)
RETURN
END
SUBROUTINE GATEMO 74/74 OPT= 1 FtN 4.6+428 79/01/09. 14.41.14 PAGE 1

1 SUBROUTINE GAtEino
C ROUTINE TO CALCULATE GATE OPENINGS
C AND PUMP SCHEDULES, IF ANY

5
COMMON /ALL/
c PHYS ICAL PARAMETERS
NMANN,SO.GGRAV,PI.BOttOM(5O).SIDSLOP(5O),RADIUS(5O),BOtGAtE(5l
MANN~5O).BOtGRAD(5O),SYFONDY(5O),SYFNLOS(5O),SPEED(5l),PUMP~
10 (51~.DEDBAND(51),NGAtES,GAtECO(ll,5l),DSINVt(5O),USINVt(5l),
c WORli,ING VARIABLES
tMAX,XEND,tINC,NXINC.NPOOLS,IPOOL,tYpp#)L,StABEG(5O),
StAEND(5O),t1,Yl,StUIlY,
C CAN; ,L- DESCRIPTION
--_.
15 UPSlA(5O,~DUNStA(5O),NOPOOLS,POOLtYp(5O),stRtYPE(Sl),
tMAXI.NUMtRNS,USRES,DSRES

REAL NMANN ,MANN


INTEGER tYPPOOL,POOLtYP,StRtYPE,StUOY
20
;
C SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES
C TO KEEP MEMORY WITHIN SWE SORT OF REASONABLE BOUND
COMMON /SCRATCH/
25 . ouMMYl(l2ooo)
. .10(50i),00(~00),YD(500)
. ,tu(5oo).ou(5oo),Yu(5oo)
,WMYY21787)

30
dRItE(6.9000j
C READ U/S 1.Q.H OF FIRST STRUCTURE
c IF THIS IS A GATE, THE DEPTH o/s OF THE GATE MUST
C HAVE BEEN SPECIFIED IN DSRES
35 C IF THIS IS A PUMP, SKIP CALCULATIONS AS 0 WAS SCHEDULED
READ(1.9010) LMAX
REAO(1.9020) (fU(I),I=l,LMAX)
READ(1.9020) (GU(I),I=l,LMAX)
READ(l,S 20) (YU(I),I=l,LMAX)
40 C WST O/S STRU I TURE IS IN ARRAYS(51)
IF(StRtYPE(5l).EG.l) GO TO 10
C CALCULATE GATE OPENINGS FOR GATE OPENING INTO
C CONSTANT LEVEL RESERVOIR At D/S END OF REACH
YDN=DSRES
45 IPWL=O
&JPoGL=l
KMAX=O
_-...
fSlRTYP=2
WRItE(6.9030) IPOOL,KMAX,IStRtYP
50 WRXtEf6.9030) JPOOL.LMAX
C CALCUiAik iiAiE 6PkNINGS it TIMES TttINC
tMAX=tU( LYAX)
C LIMIT TO 24 HOURS
IF(tMAX.Gt.86400.) tMAX=66400.
55 t=o.o
C CALCULATE-AT 5 MINUTE INTERVALS
tINC*300.
SUBROUTINE GATEMO 74/74 OPT= 1 FfN 4.6+426 79/01/09. 14.41.14 PAGE 2

CONilNUt
C41NTERPOLATE UPSTREAM
60 M=LMAX-1
DO 7 I*l,M
K=I
tF(t.,GE.TU(I).ANO,T.LE.TU(I+l)) GO TO 8
7 CONTINUE
65 6 CONTINUE
Ot=t-TO(K)
DELf=fU(K+l)-tU(K)
RATt0=DT/OELT
QUP=OU(K)+RATtO*(VU(K+l)-GU(K))
70 C 0 THRU GATE IS 0 U/S OF GATE LESS TURNOUT 0 AT GATE
IPOOL=l -
CALL TRNOUT(0ELV.T)
IPOoL=O
ODN=OUp-DE LO
75 YUPFYU(K)+RATIO+(YU(K+l)-YU(K))

: x:0 SIDE
Ct$N~E,IN
OF GATE
INVERT ELEVATION COMING INTO GATE

EG:-UStNVT(S1)
90 SYPHON=O.O
YSYPHON=O. 0
1POOL=IPoOL+l
TYPPOOL=POOLTYP(IP00L)
IUPROWN= 1
85 CALL GATEY(VUP,YUP,EC,EG,SYPHON,YSYPnON,YGATE,GDN,IUPRDYN)
.I YUPFYGATE
3
IPOOL=IPOOL-1
C CHECK THAl U/S DEPTH IS GREATER THAN D/S DEPTH
IF(YUP.LT.YDN) GO TO 500
90
: CALCULATE GATE OPENING
G=YUp
tF(ABS(YUP-YDN).LT.O.Ol) GO TO 11
9 GPREV=G
95 GOY.G/YUP
IF(GOY.LT.O.0) GOYFO.0
IF(GOY.GT.l.0) GOY=l.O
C SET IPOOL FOR COEFFICIENT ROUTINE
IPOOL=
100 COEFF=CD(GOY)
C RESET IPOOL
IPOOL=O
AON=BOTGATE(Sl)*YDN
C SET IPOOL FOR AREA ROUTINE
105 IPOoL=IPOOL+l
TYPPOOL=POOLTYP(IP00L)
AUP=YUP*BOTGATE(Sl)
C RESET IPOOL FOR D/S POOL
IPWL=IPWL-1
110 G=QON/(BOTGATE(ll )+COEFf !*SGRl ‘(2, 5l GGRAV* (YUP-YDN)
+QDN**2/AUP**2))
IF(G.GT.YUP) G=YUP
G=O:S*(G+GPREV)
IF(ABS(G-GPREV).GT.O.OOl 100 TO 9
SUBROUTINE GATEMO 74/14 OPT= 1 FTN 4.6+429 79/01/09, 14.41.14 PAGE 3

115 11 CONTINUE
WRIfE(2,9040) f.G,YUP,YDN
T=T+T INC
1F(f.LE.TMAX) GO TO 4
C READ IN DOWNSTREAM AND UPSTREAM T,Q,H OF GATE
120 C INDEX IPOOL IS ON DOWKSTREAM SIDE OF GATE

: GATE TYPE IS DATA ON DOWNSTREAM SIDE OF GATE

: ISTRTYP=O POOL NOT USED


125 C =l PUMPS
92 NORMAL GATE
E 93 SYPHDN DOWNSTREAM OF GATE
94 SYPHON WITH NO GATE
: =5 NO GATE
130 C
10 CONTINUE
READ(1.9030) IPOOL,KMAX,ISTRfYP
1F(EOF(l)) 400,20
20 CONTINUE
135 WRITE(6.9030) IPODL,KMAX,1STRTYP
READ(l.9020) (TD(I),I=l,KMAX)
REA0(1,9020) (~D(I),I=T,KKAx)
REA0(1,9020) (YO(I),I=l,KMAX)
READf 1.9030) JPOOL,LMAX
140 IF(EDF(1)) 35.30
30 WRITE(6.9030) JPDOL,LMAX
READ11.9020) (TU(I),I=l,LMAX)
i4 READ11.9020) (QU(I).I=l,LMAX)
READf1.9020) (YU(I),I=l,LMAX)
145 GO TO 36
C F1NISHED IF LAST STRUCTURE IS NOT STROKED
35 IF(SlUDY.EO.l) GO TO 400
36 CONTINUE
C SKIP ND GATE BOUNDARIES
150 IF(ISTRTYP.EQ.4.0R.ISTRTYP,E0.5) GO TO 10
C CALCULATE GATE OPENINGS OR PUMP DISCHARGES AT TIMES T+TINC
lMAX=lU( LMAX)
C LIMIT TO 24 HOURS
IF(TMAX.GT.86400.) TMAXiBS400.
155 T=O.O
C CALCULATE Al 5 YlNUTE INTERVALS
TINC=300.
40 CONTINUE
C INTERPOLATE DOWNSTREAM
160 MaKhlAX-1
Do 50 1=l,Y
K-1
IF(T.GE.TO(I).AND.T.LT.TD(I+l)) GO TO 60
CONTINUE
165 it8 CONTINUE
DT=T-fD(K)
DELT=fD(K+l)-TD(K)
RATIO=DT/DELT
QDN=OD(K)+RATIO*(OD(K+l)-QD(K))
170 YDN*YD(K)+RATIO+(YD(K+l)-YD(K))
c INTERPOLATE UPSTREAM
SUBROUTINE GATEMo 74174 OPl=l FtN 4.6+426 70/0~/00. 14.41.14 PAGE 4

C SET U/.5 OEPTH JUST IN CASE THIS IS A GATE At MOST U/S END Of REACH
YUP=USRES
IF(tPOOL.EO.NPOOLS) Go TO 85
175 MrLMAX-1

i?t’” *=leM
I;(T.GE.TU(t).AND.T.LE.fU(t+l)) Go TO 66
70 CONTINUE
180 80 CONTINUE
DT=T-TU(K)
DELT=TU(K+l)-TU(K)
RATID=DT/DELT
DUP=GU(K)+RATIO+(GU(K+l)-OU(K))
165 YUP=YU(K)+RATIO*(YU(K+l)-YU(K))
65 CONTINUE
IS THIS A PUMP
tF(tSTRTYP.ED.1) GO TO 600
IF(ISfRTYP.NE.2) GO TO 200
190 CALCULATE Y IN NORMAL GATE STRUCTURE
SET CHANNEL X-SECTION TYPE FOR AREA ROUTINE
TYPPOOL=PODLTYP( IPOOL)
IPOOL HAS BEEN READ IN ABOVE
ZERO OUT SYPHON PARAMETERS
195 SYPHON=O.O
YSYPHON=O. 0
D/S SIDE OF GATE
INVERT DROP TO TRAPEZOIOAL CHANNEL
EC IS CHANNEL INVERT ELEVATION
200 EG IS GATE INVERT ELEVATION

EC=O. 0
EG=DSINVl( IPOOL)
IUPRDYN r0 FOR D/S SIDE OF GATE, =l FOR U/S SIDE
206 IUPRDWN=O
CALL GATEY(ODN,YDN,EC,EG,SYPMON,YSYPHON,YGATE,OON,IUPRDUN)
YDN=YGATE
u/sSIDE OF GATE
IF(IPDOL.EO.NPOOLS) GO TO 310
210 SET CHANGE IN INVERT ELEVATION COMING INTO GATE
EC=O.O .
EG=-UStNVl(IPOOL)
IPooL=tPooL+1
TYPPDOL=POOLTYP(lPOOL)
215 IUPRDWNsl
CALL GATEY(OUP,YUP,EC,EG,SYPHON,YSYPHON,YGATE,GDN,IUPRDWN)
. _. -. _
VUV.VGAlP
IPooL=IPooL-1
TYPPODL=POOLTYP(IPOoL)
220 GO TO 310
200 IF(tSTRTYP.NE.3) Go TO 10
C SYPHON
C
C SET CHANNEL X-SECTION TYPE FOR AREA ROUTINE
225 TYPPOOL=POOLTYP(IPOOL)
C D/S SIDE OF GATE
C SET SYPHON HEAD LOSS AND SYPHON DROP
SYPWON=SYFNLOS(tPOOL)
SUBROUTINE GATEMO 74/74 OPT= 1 FTN 4.6+420 ro/ot/oo, (4.41.14 PAQE S

YsYPHON=SYFONDY(IPOOL)
230 C SET GATE AND CHANNEL ELEVATIONS
C INVERT DROP FROM GATE THRU SYPHON LUMPED IN SYFDNDY
EG=O. 0
EC=O. 0
IUPRDWN=O
239 CALL GATEY(ODN,YDN,EC,EG,SYPHON, ,YSYPHDN.YGAl , IUPRDYN)
YDN=YGATE
c u/s SIDE OF GATE
lF(IPODL.EQ.NPDOLS) GO TO 310
C SET CHANGE IN INVERT ELEVATION GOING INTO GATE
240 SYPHON=O ..O
YSYPHON*O. 0
EC-O. 0
EG=-USINVT(IPOOL)
IPOOL=IPODL+l
245 TYPPOOL.POOLTYP(IPOOL)
IUPRDWN* 1
CALL GATEY(QUP.YUP,EC,EG,SYPHON,YSYPHON,YGATE,OD~,tUPRDUN)
YUPsYGATE
IPDOL=IPOOL-1
250 TYPPOOL=POOLTYP(IP00L)
310 CONTINUE
C CHECK THAT U/S DEPTH IS GREATER THAN D/S DEPTH
J IF(YUP.LT.YDN) GO TO $00
n
255 ; CALCULATE GATE OPENING
C DEPTHS AND OS ARE INSIDE GATE SECTION
G*YUF
IF(AElS(YUP-YDN).LT.O.Ol) GO TO 330
320 GPREv=G
290 GOY=G/YUP
lF(GOY.LT.0.0) GDYFD.0
IF(GOY.GT.l.O) GDY=l.O
COEFF=CD(GOY)
TYPPOOL=PDOLTYP(IPDOL)
299 _..
ADN~YDN*BOTGATE(IPOOL)
C GET AREA IN U6iTREAii‘SIDEmilF GATE SECTION
AUP=YUP*9OTGATE(IPDOL)
G=QDN/(9OTGATE(XPDOL)*COEFF+SORT(2.*OaRAV~(YUP-YDN)
+QDN+*2/AUP**2))
270 IF(G.GT.YUP) G=YUP
G=O.S*(G+GPREV)
IF(AGS(G-GPREV).GT.O.OOl) GO TO 320
330 CONTINUE
WRITE(2.9040) T.G,YUP,YDN
275 T-T+TINC
IF(T.LE.TYAX) GO TO 40
GO-TO 10
400 REWIND 1
REWIND 2
290 RETURN
C REVERSE FLOW THRU GATE - STOP
500 CONT I RUE
WRlTE(6,9050) fWOL,JPDOL,XSTRTYP,T,YUP~YDN

2es
SUBROUTINE GATEMO 74/74 OPT= 1 FTN 4.6+426 79/01/09. 14.41.14 PAQE 6

600 CON1 INUE


WRIfE(2.9040) T,DDN,YUP,YDN
T=T+TINC
IF(T.LE.TMAX) GO TO 40
290 GO TO 10
c FORMATS
9000 FIhlAT(lX.9HIN GATEMO/lX.
30H POOL KMAX STR TYPE/)
9010’ FORMAl(lX.lOX,IlO)
295 9020 FORMAl(lX,1OF10.3)
9030 FORMAf(lX.4110)
9040 FORMAfilX,4FlD;3)
9050 FORMATflX,43liU/S DEPTH LESS THAN D/S DEPTH BETWEEN POOLS,
. 15,4H AND,fS/lX,9HGATE TYPE,XS/lX,SHAT T*,F10.3
300 . /lX.lOHU/S DEPTH=,F10.3/1X,lOHD/S DEPTH-,FlO.B)
END
SUBROUTINE GATEY 74/74 OPT= 1 FTN 4.6+420 7o/01/0e. 14.41.14 PAGE 1

1 SUBROUTINE GATEY(OCHANNL,YCHANNL,ECHANNL,EGATE,SYPHON,YSYPHON,
YGATE.QGATE,IUPRDWN)
C ROUTINE TO CALCULATE DEPTH IN GATE STRUCTURE FROM CHANNEL DEPTH
C
5 COMMON /ALL/
C PHYSICAL PARAMETERS
. NMA~~.S~.GGRAV,~I,BOTTOM(~O),SIDSLOP(~O),RADIUS(~O),BOTGATE(~I),
MANN(~D).~OTGRAD(~O),SYFONDY(~O),SYFNLOS(SO),SPEED(~~),PUMP~
: (~~,.OED~AND(§~).NGATES,GATECO(~~,S~),DSINVT(SO),USIN~T(~~),
10 C WORKING VARIABLES
. TMAX.XENO,yINC,NXINC,NPOOLS,IPOOL,TYPPWL,STABEG(5D),
. STAEND(SO),Il.Yl,STUDY,
C CANAL DESCRIPTION
. UPSlA(50),DWNSTA(50),NOPOOLS,PDOLTYP(5O),STRTYPE(~l),
15 TMAXI.NUMTRNS,USRES,DSRES
c -
REAL NMANN,MANN
INTEGER TYPPOOL,POOLTYP,STRTYPE,STUOY
C
20 e
C IUPRDWNIO FOR D/S SIDE OF GATE, -1 FOR U/S SIDE
IF(IUPRDWN.EQ.1) GO TO 20
C
i D/S SIOE OF GATE
25 C
Y=YCHANNL
ACHANNL=A( Y)
Cl=~CHANNL**2/(2.*GGRAV*ACHANNL+*2)+YCHANNL+ECHANNL-EGAtE
-YSYPHDN+SYPHON*QCHANNL**2
30 C GATE 5OTTOM WIDTH
B=BOTGATE(IPODL)
C2=QCHANNL**2/(2.*GGRAV*B**2)
CON1 INUE
C%E RECTANGULAR SECTION At GATE
35 IF(Y.LE.0.) Y=O.l
F=Y+C2/Y+*2-Cl
FPRIlE=I .-2.+c2/y**3
YPREV=Y
Y=Y-F/FPRIME
40 IF(ABS(Y-YPREV).GT.O.OOl) Go TO 10
YGATE=Y
RETURN

: U/S SIDE OF GATE


45
r2;Poo, AND TYPPDDL WERE SET IN GATEMD FOR THE U/S POOL
CONTINUE
YFYCHANNL
ACHANNL=A( Y)
50 C RESET TO D/S POOL FORGATE PROPERTIES
IPOOL=IPDDL-1
IF(IPODL.EQ.0) IpOOL=Sl
TYPPOOL=PDOLlYP(IpooL)
Clt~CHANNL*+2/(2.*GGRAV*ACHANNL++1)+YCHANNL+ECHANNL-EGATE
55 C GATE BOTTOM WIDTH
B=BOTGATE( IPOOL)
C2rQGATE+*2/(2.+GGRAV+Bs+2)
SUBROUTINE GATEY 74174 OPT* 1 FTN 4.6+428 79/01/09. 14.41.14 PAGE 2

C RESET TO U/S POOL BEFORE RETURN TO GATEMO


IF(IPOOL.EO.51) IPODL=O
60 IPDOL=IPOOL+l
TYPPOOL=POOLtyP( XPOOL)
CONTINUE
C%E RECTANGULAR SECTION AT GATE
IF(Y.LE.0.) Y=O.l
65 F=Y+C2/Y**2-Cl
FPRIMEsl.-S.+C3/Y**3
YPREV=Y
Y-Y-F/FPRIYE
~F~AS~(Y-YPREV).GT.~.~~~) GO To 30
70 YGATE=Y
RETURN
END
FUNCTION CD 74/74 OPT= 1 FtN 4.5+428 70/01/00, 14.41.14 PAW 1

..
1 FUNCTION CD(GOY)
ROUTINE TO CALCULATE GATE COEFFICIENT AS FUNCTION
: OF GATE OPENING/YUPStREAM
C
5 COMMON /ALL/
C PHYSICAL PARAMETERS
. NMANN.SO,GGRAV,P~.~O~TOW(~O),SIOSLOP(SO),RADIUS(~O),~OTGATE(S~),
. MANN(~O).~OTGRAD(~O),SYFOKDY(~O),SYFNLOS(~O),SPEED(~~),PUMP~
(~~~.OEDBAND(~~).NGATESIOA~ECO(~~~S~),DSINVT(~O),U~INVT(~~),
10 C WORiING VARIABLES
. ~THAX.XENO,TINC,NXINC,NPOOLS.IPOOL,TYPWOL,STA5EG(5O),
STAEN0(50),Tl,Yl,STUDY.
C CANiL DESCRIPTION
. UPSTA(5O),DWNSTA(5O),NOPOOlS,POOLTYP(50),STRTYPE(5l),
15 . TMAXI,NUMTRNS,USRES,DSRES
C
REAL NMANN #MANN
INTEGER TYPPOOL.POOLTYP,STRTYPE,STUDY
C
20
E FIT PARABOLA TO 3 POINTS NEAR GOY
K=GOY*lO.+l.
IF(K.LT.l) K=l
IF(K.GT.9) K=9
25 OFFSEt*K
X=GOY*lO.-OFFSET
I=IPOOL
Al~0.5*(GATECO(K,I)-2.*GATECO(K+l,I)+GATECO(K+2,I))
~1+0.5+(GATECO(K+2,1)-GATECO(K.1))
30 Cl =GATECO( K+l , I)
CD=Al*X+*2+5l*X+Cl
RETURN
END
SUBROUTINE FOLOWER 74/74 OPT= 1 FTN 4.5+428 79/01/09* 14.41.14 PAGE 1

1 SUBROUT IN5 FOLOWER


ROUTINE TO FOLLOW GATE STROKING GATE MOTIONS TO DETERMINE
c” START TIMES AND FINAL POSITIONS OF GATES

5 : OR TO TURN PUMPS ON AND OFF

C SPEED IS GATE SPEED IN Ft/MIN


DEDBAND IS DIFFERENCE ALLOWED IN GATE OPENING BETWEEN
E ACTUAL POSITION AND STROKING SOLUTION BEFORE MOVING GATE
C
IO C PUMPDO IS INCREMENTAL DISCHARGE FOR PUMPS
COMMDN /ALL/
C PHYSICAL PARAMETERS
. NMANN.SO.GGRAV.PI,BOTTOM(5O),SIDSLOP(50),RADIUS(5O),5OTGATE(Sl),
. MANN(50),5OTGRAO(50),SYFONDY(5O),SYFNLOS(5O),SPEED(5l),PUMP~
15 (5l~~DEO5AND(51).NGATES,GATECO(ll,5l),OSXNVT(S0),USfNVt(~l),
C WORhNG VARIABLES
, TMAX.XENO,TINC,NXINC,NPOOLS,IPOOL,TYP~L,STA~LG(~O),
. STAEND(5O).Tl.Yl.STUDY.
C CANAL DESCRIFTION
20 . UPSTA(5O),DWNSTA(5O),NOPOOLS,POOLTYP(5O),STRTYP~(5l),
. TMAXt,NUMTRNS,USRES,DSRES
C
REAL NYANN .MANN
INTEGER TYPPOOL.POOLTYP,STRTYPE,STUOY
25
c”
SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES
: TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
COMMON /SCRATCH/
30 . 0wMY1(12000)
. .T(500).G(5OO),TMIN(5OO),GPOSN(5OO)
. .DUMMY2(1757)
C
EOUIVALENCE 0 AND G SO PUMP STROKING VARIABLE NAMES MAKE
35 : . SOME SENSE
XF THIS 1S.A PUMP, GATEMO WROTE 0; IF IT IS A GATE, GATEMD
E WROTE GATE OPENING, G
C
DIMENSION O(5OO),OLAST(500)
40 EDUIVALENCE (O(l),G(l)),(OLAST(1).GPOSN(1)),(PUMPG,GATEOPN)

E
C
REWIND 2
45 REWIND 5
IPooL=Sl
C READ IN T AND GATE OPENING
NOATE=
10 1-l
50 IEOF=O
REA0~2.9010) T(I).GtI)
IF(EOF(2)) 70.20
20 I=I+l
READf2.9010) T(l),G(I)
55 IF(EOF(2)) 70,30
30 IF(T(I).NE.O.O) GO TO 20
BACKSPACE 2
SUBROUTINE FOLOWER 74/74 OPT= 1 FTN 4.6+429 79/01/09. 14.41.14 PAOE 2

40 MAXal-
c CALCULATE START 1 IMES
60 C INITIAL OPENING
c TIME IN MINUTES NOW
TMINl 1 )=O.D
GPOSN(l)=G(l)
GATEOPN=GPOSN(l)
65 I=2
C FIND POOL D/S OF THIS GATE
42 IF(SlRlYPE(IPOOL).ED.2) GD TO 45
IF(SlRTYPE(IPOOL).EO.3) GO TO 45
C SKIP STROKING MOST D/S PUMP
70 IF(IPOOL.EQ.51) GO TO 44
IF(slRTYPE(IPOOL).EO.1) GO TO 80
44 CONTINUE
IF(IPOOL.EO.51) IPOOL=D
IPOOL~IPooL+1
75 GO TO 42
45 CONTlNUE
DO 50 J=l.MAX
~F(A~S(GATEOPN-Q(IJ)).LT.DEO~AND(IPDOL)) GO TO 59
c MOVE GATE
60 DELY=ABS(GATEOPN-G(d))
OELT=OELY/SPEED(IPOOL)
TMINfI)=T(tJ)/60.-DELT
GPOSN(I)=G(J)
GATEOPN=G( d)
65 x=1+1
50 CONTINUE
NOMOVE=I-1
WRITE(6.9020) NGATE.NOMDVE
C SAVE FOR ROUTINE WRITElT TO OUTPUT IN REVERSE ORDER
90 WRITE(5.9020) NGATErNOMOVE
NGATES=NGATE
DO 60 I=l,NDMOVE
WRITE(5.9010) TYIN(I),GPOSN(I),SPEED(IPOOL)
60 CONTINUE
95 WRITE(6.9000) DEDBAND(IPOOL)
NGATE=NGATE+l
IF( IPOOL.EO.51) IPODL=O
IPOOL=IPODL+l
IF(IEOF.EQ.0) GO TO 10
100 REWIND 5
RETURN
t0 IEOF= 1
GO TO 40
c PUMP STROKING
105 80 CONTINUE
C SET INITIAL DISCHARGE TO N*PUMPDO
N=O
65 CONTINUE
FN=N
110 DTEST=PUMPQ-FN+PUMPDQ(IPOOL)
IF(DTEST.LT.O.6) GO TO 90
N=N+l
GO TO 65
90 N=N-1
SUBROUTINE FOLOWER 74/14 OPT= 1 FTN 4.6+428 79/01/09. 14.41.14 PAW 3

115 FN=N
PUMPQ=FN*PUMPOG(IPOOL)
GLAST(l)=PUMPG
DO 110 J=2,MAX
IF(ARS(PUt.lPO-G(J)).Lt.PUMPOG(IPOOL)) GO TO 110
120 C TURN PUMP ON OR OFF
IF(PUMPO.LT.G(J)) GO TO 100
C TURN PUMP OFF
TMINII)=f(J)/60.
95 CON1 INUE
125 GLAST(l)=PUMPG-PUMPbG(IPOOL)
PUMPO=OLASt( I)
IF(ABS(PUMPG-G(d)).GE.PUMPOG(IPOOL)) GO TO 95
I=I+l
GO TO 110
130 C TURN PUMP ON
100 CONTINUE
TMlkl 1 )=Y J)/80.
105 CONTINUE
OLAST~I)=PUMPG+PUMPOG(IPGOL)
135 PUMPO=OLAST(I)
IF(AOS(PUMPG-O(J)).GE.PUMPOG(IPOOL)) GO TO 105
1=X+1
110 CONTINUE
NGMOVE=I-1
140 WRITE(6.9030) NGATE,NOMOVE
WRITE(5.9030) NGATE,NOMOVE
NGATES=NGATE
DO 120 I=l,NOMGVE
WRITEf5.9040) TYIN(I),OLAST(I)
145 120 CONTINUE
WRITE(6.9050) PUMPGG(IPOOL)
NGATE=NGATE+l
IPOOL=IPooL+1
lF(lEOF.EG.0) 00 TO 10
150 REWIND 5
RETURN
C FORMAT5
9000 FORYAT(lX.24HGATE MOVEMENT OEADBAND n ,F10.3,3H FT)
9010 FORWAT(lX,3F10.3)
155 9020 FORMAf(lX,4HGATE.I5,IlO,2X,9HMOVEMENTS)
9030 FORMAT(lX.4HPUMP,l5,110,2X,9HCHANGES )
9040 FORMAf(lX,2F10.3,1OX)
9050 FORMAT(lX.14HPUMP DELTA 0 =,F10.3)
EN0
SUBROUTINE WRITE11 74/74 OPT= 1 FTN 4.6+428 79/01/09. 14.41.14 PAGE 1

1 SUBROUTINE WRITE11
C ROUTINE TO MERGE TURNOUT SCHEDULES, GATE SCHEDULES
C AND PUMP SCHEDULES INTO INPUT FILE ‘FOR UNSTEAOV MODEL
C FILE ON TAPE 6 IS POSITIONED AT END OF INITIAL CONDITIONS
5 C OUTPUT IS FROM UPSTREAM TO DOWNSTREAM
C
COMMON /ALL/
C PHVSICAL PARAMETERS
. NMANN,SO,GGRAV,PI,BDTTOM(SO),SIDSLOP(50),RADIUS(5O),9OTGATE(Sl),
10 MANN(5O).~DTGRAD(SO),SVFONDV(SO),SVFNLOS(5O),SPEED(Sl),PuMP~
: (5l,.DED~AND(Sl,.NGATES,GATECO(ll,5l),DSINVT(5O),USINVT(5l),
C WORKING VARIABLES
TMAX.XEND,TINC,NXINC,NPODLS,IPOOLltYPPM)L,STABEG(5O),
: STAEND(SO),Tl,Vl.STUDV,
15 C CANAL DESCRIPTION
. UPSTA(50),DWNSTA(5O),NDPODLS,PDDLTVP(SO),STRTVPE(Sl),
TMAXI,NUMTRNS,USRES,DSRES
c -
REAL NMANN,MANN
20 INTEGER TVPPOOL,POOLTVP,STRTVPE,STUDV
C

E SCRATCH AREA TO BE OVERLAID IN THE VARIOUS SUBROUTINES


C TO KEEP MEMORY WITHIN SOME SORT OF REASONABLE BOUND
25 CWMON /SCRATCH/
DUMNlV1(8618),0D(48),0U(48),0t(48,50),DUMMV2(886),
: TMOVE(SOO),GO(SOO),SPED(500),DUMMV3(2287)
C
DIMENSION O(500)
30 C EOUIVALENCE GATE OPENING, GO, FOR GATES WITH 0 FOR
C PUMPS SO VARIABLE NAMES WILL MAKE SOME SENSE.
C Go AND 0 ARE IN THE SAME FIELD ON TAPE 5.
C
EOUIVALENCE (GO(l),O(l))
35 C
C START AT UPSTREAM END
NGATE=NGATES
JPDoL= 1
K=NOPOOLS
40 C WAS U/S STRUCTURE A GATE OR A PUMP
IF(STRTVPE(K).EO.l) Go TO 50
c GATE
C FIN6 GATE SCHEDULE ON TAPE 5
- _ 5
REWIND
45 10 READ~5.9110) IGfiTE,IMOVE
Do 20 I=l.IMOVE
READ(5.9120) TMOVE(I),GO(I),SPED(I)
20 CONTINUE
IF(IGATE.NE.NGATE) GO TO 10
50 c SET FOR NEXT D/S STRUCTURE
NGATE=NGATE-1
WRITE(9.9130) JPOOL,STRTVPE(K),IMOVE
DO 30 I=l,IMOVE
WRITE(9.9140) T~VE(I),GO(I),SPED(I)
55 30 CONT 1NUE
GO TO 120
c PUMP
SUBROUTINE WRltEIf 74/74 OPT= 1 FTN 4.6+426 79/or/oo. 14.41.14 PAGE 2

60 CONTINUE
C WAS PUMP SCHEDULED OR STROKE0
60 fF(STUOY.EG.2) GO TO 60
C SCHEDULED
C EXTRACT Q CHANGES FROM QU
XMOVE=l
TMOVEfl)=O.O
65 O(l)=OU(l)
dti 60 J+2,46
IF(Ou(J).EO.OU(J-1)) GO TO 60
IMOVE= IMOVE+l
TMOVE(IMOVE)=(J-1)*60
70
60 itiNT It&
WRITE(9.9130) JPOOL,STRTYPE(K),IMOV6
00 70 I-1, IMOVE
WRITE(9.9150) TMOVE(I),O(X)
75 70 CONTINUE
GO TO 120
C STROKED
CONTINUE
C”!:IND PUMP SCHEOULE ON TAPE 3
60 REWIND 5
90 REAOI5.9119) IGATE,IMOVE
00 100 I-1,IMOVE
REAOt5.9120) TMOVE(I),O(I)
100 CONTINUE
65 IF(XGATE.NE.NGATE) GO TO 90
C SET FOR NEXT O/S STRUCTURE
NGATE=NGATE-1
WRITE(9.9130) JPOOL,STRTYPE(K),IMOVE
00 110 I=lrIMOVE
90 WRITE(9,9150) THOVE(I),O(I)
110 CONTINUE
c IS THERE A TURNOUT AT D/s END OF POOL
120 CONTINUE
DO 130 1=1,46
9s IF(Qf(I,K).NE.O.O) GO TO 140
130 CONTINUE
GO TO 200
C,E;TRACT TURNOUT Q CHANGES.FROH QT
CONTINUE
100 IMOVE= 1
TMOVE(l)=O.O
O(l)=OT(l,K)
DO 150 h2.46
IF(Qltd,K).EO.QT(d-1.K)) GO TO 150
10s IMOVE= IMOvE+l
TMOVE(IMOVE)=(J-I)*60
~‘J;:;=O”J,W
150
ITRI NTYP=99
110 WRITE(9.9130) JPOOL;XTRNTYP,IMOVE
DO 160 181 ,IMOVE
WRITE(9.9+10) TMOVE(S),O(l)
160 CONTINUE
C
SUBROUTINE WRITEIT 74/74 OPT= 1 FTN 4.6+426 79/01/0s. 14.41*14 PAGE 3

115 C INTERIOR BOUNDARIES


C
200 CONTINUE
REWIND 5
K=K-1
120 JPOOL=JPOOL+l
IF(K.EO.0) GO TO 900
C IS THERE A STRUCTURE AT U/S END OF THIS POOL
IF(STRTYPE(K).EQ.4) GO TO 300
IF(SfRTVPE(K).EO.S) GO TO 300
125 IF(STRTVPE(K).NE.l) GO TO 240
c PUMPS
C FIN0 PUMP SCHEDULE ON TAPE 5
210 REAO(5,9110) IGATE.IMOVE
DO 220 I=1 ,IYOVE
130 REAOf5.9120) TMOVE(I).O(I)
220 CONTINUE
IF(IGATE.NE.NGATE) GO TO 210
C SET FOR NEXT D/S STRUCTURE
NGATE=NGATE-1
135 WRITE(9.9130) JPOOC,STRTvPE(K),IMOvE
DO 230 I=1 ,IMOVE
WRITE(9.9l50) TMOVE(I),O(I)
230 CONTINUE
co to 300
140 C GATE -- -
240 CONTINUE
C2;N0 GATE SCHEDULE ON TAPE 5
REAOl5.9110) IGATE,IMOVE
DO 260 I=1 ,IMOVE
14s REAO(5.9120) TMOVE(I),GO(I),SPEDtI)
260 CONTINUE
IF(IGATE.NE.kriATE) GO TO 250
C SET FOR NEXT O/S STRUCTURE
NGATE=NGATE-1
150 WRITE(9,9130) JPOOL;STRTVPE(K),IMOVE
DO 270 I=l.IMOVE
WRITE(9.9140) TMOVE(I).GO(I).SPED(1)
270 CONTINUE
C3;t THERE A TURNOUT AT D/S END OF THIS POOL
155 CONTINUE
DO 310 X=1,46
IF(OT(I.K).NE.O.O) GO TO 320
310 CONTINUE
GO TO 200
160 C3EyRACT CHANGES IN TURNOUT Q FROM QT
CONTINUE
IMOVE=
TMOVEfl)=O.O
Ofl)=OT(l.K)
160 DO 330 Ja2.46
IFtOT(d.K).EQ.Ot(J-1.M)) GO TO 330
IMOvE=IMOvE+l
TMOVE( IMOVE)=(J-1)*60
O(IMOvE)=GT(J,K)-OT(J-1.K)
170 330 CONTINUE
ITRNTVP=99
SUBROUTINE WRITEIT 74/74 OPT=1 FTN 4.6+420 79/01/09. 14.41.44 PAGE 4

WRITE(9.9130) JPOOL,ITRNTYP,IMOVE
00 340 I=1 ,IMOVE
WRITE(9.9150) TMGVE(I),G(I)
175 340 CONTINUE
GO TO 200
C
C MOST O/S STRUCTURE
C
100 900 CONTINUE
C Q WAS SCHEOULEO IN GO
C IS THIS A PUMP
IF(STRTYPE(5l).NE.l) GO TO 950
c PUMP
195 C EXTRACT CHANGES IN 0 FROM QO
IMOVE=l
TMOVE(l)=O.O
0(1)=00(1)
00 910 J+2,49
190 IF(OO(J).EO.OO(J-1)) GO TO 910
IhlOVE=IMOVE+l
TMOVE(IMOVE)=(J-l)+BO
O( ImovE,=Go(d)
910 CONTINUE
195 WRITE(9.9130) JPOOL,STRTYPE(Sl),IMOVE
00 920 I=1 *IMOVE
WRITE(9,9150) TMOVE(I),G(I)
920 CONTINUE
co GO TO 1000
cn 200 C GATE
950 CONTINUE
REAO(5.9110) IGATE,IMOVE
WRITE(9.9130) JPOOL,STRTYPE(51),IMOVE
00 960 I=l,IMOVE
205 REAO(5.9120) TMOVE(I),GO(I),SPEO(I)
WRITE(9.9140) TMOVE(I),GO(I),SPEO(I)
960 CONTINUE
1000 CON7 tNUE
REWIND 5
210 REWIND 9
RETURN
9110 FORHAT(5X, 15,110)
9120 FORMAT(lX,BF10.3)
9130 FORMAT(3110)
215 9140 FORMAT(3FlG.3)
9150 FORMAT(2FlO.3,lOX)
END

You might also like