8 2 JUL 1974
Approved for Public Release;
Distribution Unlimited.
ENVPREDRSCHFAC
Technical Note No. 1-74
A VERTICALLY INTEGRATED
HYDRODYNAMICAL-NUMERICAL MODEL
(W, Hansen type)
MODEL DESCRIPTION AND OPERATING/RUNNING INSTRUCTIONS
By
OCEANOGRAPHY DEPARTMENT
ENVIRONMENTAL PREDICTION RESEARCH FACILITY
Part 1 of a series of four reports
prepared for:
ENVIRONMENTAL PROTECTION AGENCY
PACIFIC NORTHWEST ENVIRONMENTAL
RESEARCH LABORATORY
CORVALLIS, OREGON
JANUARY 1971
ENVIRONMENTAL PREDICTION RESEARCH FACILITY
NAVAL POSTGRADUATE SCHOOL
MONTEREY, CALIFORNIA 93940
-------
IINPI flCSTFTTp
SECURITY CLASSIFICATION OF This PAGE (When Data Entered)
REPORT DOCUMENTATION PAGE
READ INSTRUCTIONS
BEFORE COMPLETING FORM
t. REPORT NUMBER E N V P R E DRS C H FAC 2. GOVT ACCESSION NO.
Technical Note No. 1-74
1. RECIPIENT'S CATALOG NUMBER
4. TITLE (and Subtitle)
A Vertically Integrated Hydrodynamical-
Numerical Model (W. Hansen Type) Model
Description and Operating/Running
Instructions.
5. TYPE OF REPORT ft PERIOD COVERED
(. PERFORMING ORG. REPORT NUMBER
7 AUTHORf#;
Oceanography Department
Environmental Prediction Research Facilil
6. CONTRACT OR GRANT NUMBERfa)
y
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Environmental Prediction Research Facilit
Naval Postgraduate School
Monterey, California 93940
10. PROGRAM ELEMENT. PROJECT TASK
AREA ft WORK UNIT NUMBERS
y
11. CONTROLLING OFFICE NAME AND ADDRESS
Naval Air Systems Command
Department of the Navy
Washington. D.C. 20361
12. REPORT DATE
January 1974
11. NUMBER OF PAGES
69
14. MONITORING AGENCY NAME ft ADDRESSf/f different from Controlling Office)
IS. SECURITY CLASS, (ol thle report)
UNCLASSIFIED
Ita. DECLASSIFICATION/DOWNGRADING
SCHEDULE
16. DISTRIBUTION STATEMENT (ol thle Report)
Approved for public release; Distribution unlimited.
17. DISTRIBUTION STATEMENT (ol the mbetracl an farad In Block 30, If different from Report)
IB. SUPPLEMENTARY NOTES
19. KEY WOROS (Conffnua on reverie aid* II neceeeary and Identity by block number)
Hydrodynamical-Numerical Model Storm surge prediction
Transport and diffusion of pollutants Drift computations
Analysis and prediction of tides
20. ABSTRACT (Continue on ravaraa tide It neceeeary and Identity by block number)
This report presents an operational description of the
single-layer vertically integrated Hydrodynamical-Numerical
(HN) model. The model is an adaptation of one developed by
Professor Walter Hansen, University of Hamburg, and tested
over the course of two decades. The description includes
the setting up of the model for a given area and the selection
and preparation of input data. ^ brief ^ackgrouad «.f the
DD , j?nM73 1473 EDITION OF 1 NOV 68 IS OBSOLETE UNCLASSIFIED
S/N 0102-014-6601 — —— — „
.• SECURITY CLASSIFICATION OF THIS PAGe fWia« Dafa BnlaradJ
-------
UNCLASSIFIED
^LIJUITX CLASSIFICATION OF THIi FAseCHftan Dam BntarenQ
20. (continued)
hydrodynamical formulas and their finite difference forms are
also presented. Certain new features which have been added
to the program are described.
For coastal applications, it is now possible to run the
model with three open boundaries. The treatment of the
boundaries is described in this paper.
Another new feature makes possible the inclusion of the
permanent (thermohalire component) current as well as tidal
and wind currents. The new model also allows the specifica-
tion and subsequent modification of the mixed layer depth in
case there is thermal and haline stratification present in
the area. The wind current will be effective only in the
upper mixed layer.
The transport and diffusion equation is included in the
new program. Thus, the program can be used for computing the'
distribution of pollutants. In case of oil pollution or other
flotsam on the surface, the oil patch or flotsam is given an
additional push from the wind.
The program (see Appendix C) is written in FORTRAN and
thus will run on almost any computer. If the field size is
less than 1,000 grid points (e.g., 25x40) the model can be run
on a computer with a 32K memory. If the time step is not too
small, a 24-hr prediction requires about 30 minutes of computer
time on a 32K machine.
UNCLASSIFIED
SECURITY CLASSIFICATION OF THIS PMjefWhmi Data HnUr«<0
-------
FOREWORD
The sing1e-1ayer, vertically integrated Hydrodynamical-
Numerical (HN) model has been used at Environmental Prediction
Research Facility for analysis and prediction of tides, storm
surges, real-time currents and for the computation of drift
for search and rescue missions. Current components are
computed in short time intervals and, thus, the model is
ideally suited for analysis and prediction of the dispersal
and diffusion of pollutants. After testing various finite
difference forms for inclusion of diffusion into the model,
the one presented in this program description was found most
suited for models of small grid size.
The model was applied to three different size areas of
the New YorkBight. The results and verification of the
model are given in other reports in this series.
The project was funded by Environmental Protection
Agency, Pacific Northwest Environmental Research Laboratory,
Corvallis, Oregon. The project officer was Mr. R. J. Callaway,
Chief, Physical Oceanography Branch, National Coastal Pollu-
tion Research Program of the above mentioned EPA laboratory.
The project was carried out by the Oceanography
Department of Environmental Prediction Research Facility
under the supervision of Dr. Taivo Laevastu.
iii
-------
CONTENTS
FOREWORD iii
1. THE BASIC HYDRODYNAMICAL EQUATIONS AND THEIR
FINITE DIFFERENCE FORMS 1
2. THE BASIC HYDRODYNAMICAL-NUMERICAL (HN) MODEL .... 5
2.1 The Grid Net 5
2.2 Stability Criteria 7
2.3 Horizontal Viscosity 7
2.4 Bottom Friction 8
2.5 Input of Depths and Definition of Boundaries . . 9
2.6 Input of Tides
2.7 Input of Permanent (thermohaline) Current ... 10
2.8 Input of Wind and Derivation of "Layer and
a Half" Model 11
2.9 Treatment of Tidal Flats and Inundations .... 12
3. TREATMENT OF OPEN BOUNDARIES 13
4. COMPUTATION OF DIFFUSION AND DISPERSAL OF
POLLUTANTS 15
5. REVIEW OF THE PROGRAM AND OUTPUT FORMATS 23
REPORTS IN THIS SERIES 25
APPENDIX A - INPUT PARAMETERS, CARD FORMATS, AND
SET-UP OF TAPES 28
APPENDIX B - ABBREVIATIONS AND PARAMETERS USED IN
THE PROGRAM 36
APPENDIX C - PROGRAM NYBIGHT LISTING 41
v
-------
THE BASIC HYDRODYNAMICAL EQUATIONS AND
THEIR FINITE DIFFERENCE FORMS
The use of the numerical hydrodynamical method for the
computation of tides and currents was originally proposed by
Hansen in 1938.^ With the development and availability of
high-speed electronic computers the use of this method became
economically feasible and practically possible.
The derivation of vertically integrated hydrodynamical
equations is described by Hansen (1956),2 Slindermann ( 1966),3
Jensen, Weywadt and Jensen (1966),^ and Brettschneider (1967)
The following basic hydrodynamical equations are used in the
si ngle-1ayer model:
/ (x)
+ fv - W2u + £u /u2+v2 + g-jp- = X + — (1)
Hansen, W., 1938: Amplitudenverhaltnis und
Phasenunterschied der harmonischen Konstanten in der Nordsee.
Ann, d. Hydr. u Marit. Met. 66 (9):429-443.
2
Hansen, W., 1956: Theorie zur Errechnung des Wasserstandes
und der StrSmungen in Randmeeren nebst Anwendunqen. Tellus,
8:287-300.
3
Slindermann, J., 1966: Ein Vergleich zwischen der
ana.lytischen under der numerischen Berechung winderzengter
StrSmungen und Wasserstande in einem Modellmeer mit
Anwendungen auf die Nordsee. Mitteil. Inst. Meeresk., Hamburg
4:73 pp. tables & figures.
4
Jensen, H. D., S. Weywadt and A. Jensen, 1966: Forecasti ng
of storm surges in the North Sea. Part I. NATO Subcomm.
Oceanogr. Res. Techn. Rpt. 28 (Mimeo).
5
Brettschneider, G., 1967: Anwendung des hydrodynamisch-
numerischen Verfahrens zur Ermittlung der M„ - Mitschwinqung-
sgezeit der Nordsee. Mitteil. Inst. MeeresT., Hamburg, 7:65 pp
and figures.
- I -
-------
a_v
at
- fu - vV2v + £v /u2+v2 + g-^- = Y +
n y h
(2)
If + h (Hu) ~ hr (Hv)
9y
(3)
Equations (1) and (2) are vertically integrated equations
of motion, and equation (3) is an equation of continuity,
and (wind stress) are usually expressed as:
t(x) = A*Wy /w2 + W2
xx y
(4)
x(y) = A W/ d2 + W2
y x y
The bottom stress (friction) (t^)) in equations (1) and
(2) is:
^ u /u2 + v2; ^ v /x2 + v2 (5)
Analytical solutions to equations (1) through (3) are
of little value since exact solutions are possible only for
basins with a regular shape and constant depth and wind
distributions. However, Hansen (1956) and his later colla-
borators (e.g., SLidermann (1966) and Brettschnei der ( 1967))
have developed an explicit method for achieving time-dependent
solutions to these formulas using a finite difference approach.
The finite difference approximations for equations (1) through
(3) are:
St+T(n,m) = ct_T(n,m) - (n ,m)Ut(n,m)-HJ(n,m-l)
Ut(n,m-1) + Hj(n-1,m)Vt(n-1 ,m) - H^(n,m)Vt(n ,m) ^ ^
- 2 -
-------
Ut+2x(n,m) = j 1 - ^2Tr/Hj+2r(n,m) yOt(n,m)2+V*t(n,m)2|
Dt(n,m) + 2t fV^n.m) - p- | ;t+T(n ,m+l) - ct+T(n,m)| +
2TXt+2T(ri,m) (7)
Vt+2x(n ,m) = 1 -
2Tr/Hj+2x(n,m)
n ,m) 2+U*t (n ,m)2 |
^t (ri ,m) - 2t fU*t(n,m) - J ct+x( n ,m) - ct+x'(n+l ,m) j +
2xYt+2T(n,m) (8)
The following symbols are used in the preceding formulas
x,y space coordinates
t time
u,v components of velocity
h initial depth (when £=0)
C surface elevation
H total depth (H=h+c)
Hu I depths at u and v points respectively
Hvi
X,Y components of external forces
t(x) ,x(y) components of wind stress
g acceleration of gravity
f Coriolis parameter
r friction coefficient (bottom stress)
v coefficient of horizontal eddy viscosity
2
V Laplacian
X coefficient of friction (drag coefficient)
W ,W wind speeds
x y K
bottom stress
n ,m coordinates of the grid point
t half time step | ^ f|n^e difference equations
l grid length
- 3 -
-------
The averaged velocity and water elevation (sea level)
components are:
Ot(n,m) = a Ut( n ,m)+l^ j Ul { n-1 ,m) + Ut( n+1 ,m)+Ut( n ,m+l ) + Ul
(n,m-l)( . (n,m) and c ^(n »m) are analogous to Dt(n,m)
above. (9)
(The factor a can be interpreted as a "horizontal
viscosity parameter." Its normal value is 0.99).
Un(n,m) = | | Ut(n,m-1)+Ut(n+1Sm-1)+Ut(n,m) + Ut(n + 1,m) | (10)
i.
V (n,m) is analogous, to the U (n,m) above.
The time step is 2x . The total depth (H^.Hy) is
computed as:
hJ + 2t (n,m) = hu(n,m) + j
St+T(n,m) + £t+T(n,m+l)
(11)
The effects of wind (external force) are computed with
the following formulas:
*' = >«' /(wV ¦ (wV ,
8P
p 3 x
(12)
x w" / (wV + (w1)2
y X7 v y'
H
i 9P
J o
p 3 y
(13)
where PQ is atmospheric pressure
- 4 -
-------
2. THE BASIC HYDRODYNAMICAL-NUMERICAL
(HN) MODEL IN FORTRAN II
2.1 THE GRID NET
The grid net (staggered grid) is shown in Figure 1.
It consists of three different sets of grid points: (1) the
water elevation points (z) at the intersection of the grid;
(2) the point for u-velocity component to the right and (3)
the V-velocity component below the corresponding z point.
Each of these three points has the same coordinate designation
(n,m).
The geographical orientation of the grid is usually
selected so that the x axis (m coordinate) is parallel to
the entrance of the bay (or to the tidal input boundary). The
relationship between the geographic coordinates and computa-
tion coordinates is shown in Figure, 2. This orientation must
be considered in the input of wind directions, which are
introduced in computational grid coordinates. The direction
of the current is printed (and plotted) in geographical
direction.
The correction of the current direction to geographic
direction is made within the program. The geographic directior
is computed with the following formula:
0 = 9o°-0 + 0 . (14)
*geogr. vcomp. Frot '
where (8 is the desired geographical direction; 0
geogr. 3 3 r 'comp.
is the direction as computed in the program, and 0rot is the
rotation of the grid from its normal position, and this angle
is introduced as an input card. The selection of grid size
is described in the next section in connection with stability
criteria.
The coastline must pass through u and v points and not
through z points; thus, it becomes a "step line" along irregular
coasts.
- 5 -
-------
X Z-point, • U-point, o V-point
Figure 1. Scheme of the Grid Net
N,0°
- U,180°-
W,270°
/
- V, 270°
+V,90°
/
/
Angle of
rotation.
E,90°
•+u,o°
S,180°
Figure 2. Relations Between the Geographic Coordinates
and Computation Coordinates
- 6 -
-------
2.5 INPUT OF DEPTHS AND DEFINITION OF BOUNDARIES
Three sets of depth values must be read into the computer.
The first set is symbolic depth (sea and land table) at z
points: the z points over land are designated 0; at the
boundary on the coast, -1 (immediately outside the boundary
line); over the water 1; and at the input boundaries -2
(and -3). The second and third sets of depth values are read
at u and v points respectively. The depth values over the
land are 0; at the boundary points, -1; and over the sea,
the actual depth in centimeters.
All the boundary values should be defined if the program
is to work correctly. Thus, the boundaries along the coast
are defined as closed boundaries (i.e., no flow into or from
the land). This is done by setting either u or v components
at 0 on the coast (dependent on the direction of the coast)
by defining the depths at u (respectively v) points as -1 for
boundary identification purposes.
2.6 INPUT OF TIDES
The tidal heights are introduced at each time step at
the second line from the open boundary. To avoid computational
problems at coastal singularities, open boundaries should
preferably be along the X axis (M) on the top of the grid
(Y or N=2), or on the right hand side along (N) grid, whereby
the land is located to the left of this input line.
Normally, four tidal components are used for the input
of sea level elevation on the open boundary:
Z = A^cos (a.jt -
-------
In this equation, the speeds of the constituents (a)
must be given in radians per second; the epoch of the con-
stituent (k) is in radians, and the amplitude of the
constituent (A) is given in cm.
It should be noted that, by introducing only astronomical
tidal amplitudes at the open boundary, the possible effects
of atmospheric tidal components (caused by prevailing winds)
and sea-level deformation by permanent currents are neglected.
Furthermore, the tidal harmonic constants are derived from
tide recordings near the coast and their extrapolation over
deep water is uncertain. Experiments have shown, however,
that the model corrects the unsatisfactory input assumptions
quite rapidly and that reliable results can be extracted a
few grid lengths inside the input boundary. In some experi-
ments, the input tidal amplitudes have been made an inverse
function of depth. This somewhat improves the results near
the input boundary.
If the tidal constituents differ considerably at opposite
ends of the input boundary (or along a coast), a possibility
is provided in the program for the input of both sets of tidal
harmonics and their interpolation.
2.7 INPUT OF PERMANENT (THERMOHALINE) CURRENT
There are several possibilities for input of permanent
currents. The more correct (but time consuming) method is
to introduce u and v components of permanent currents at the
input boundary before tides are introduced. Computations are
then made until stability is reached in the area (a few hours
to 20 hours in real-time), at which time the u and v current
field is extracted, stored and added later in each time step,
while computing with tides and winds. This method is not
shown in the program given in the Appendix. The second and
fully satisfactory input of permanent currents is accomplished
by prescribing a slope to an open boundary, or a longitudinal
or transversal slope over the whole computational area.
- 10 -
-------
2.8 INPUT OF WIND AND DERIVATION OF "LAYER AND A HALF" MODEL
The wind speed is given on the input card in m secT^.
The wind direction must be given in grid coordinates and in
the direction toward which the wind is blowing. Thus, a
west wind is coming from -x toward +x and should be given an
angle of 0°; a north wind has an angle of 270°; east wind,
180°, and south wind, 90°. If the coordinate system is turned
(tilted) in relation to the geographic system, the angle of
tilt must be added.
In the program given in Appendix C, the atmospheric
pressure effects are neglected. This can be done if computa-
tions.are made over a small area. In larger area computations,
the effects of atmospheric pressure changes should be included
in the program. For this purpose, it can be assumed that the
sea level acts as an inverted barometer.
If the area on which the computation is made is deep or
has thermal and haline stratification, it can be assumed that
the wind current reaches only to the top of the thermocline.
An empirical approach to this effect is provided in the program.
An initial mixed layer depth must be specified and this initial
depth is changed in inverse proportion to sea-level changes.
The wind current is computed to the depth of the computed
mixed layer. When the depth is shallower than the mixed layer,
the wind current is computed to the bottom. This empirical
approach is not fully satisfactory in all instances as the
change of MLD set must be computed from divergence, which is
derived from transport computations. This is properly done
in multi-layer HN models (see Laevastu, 1 974a)^. Furthermore,
it should be borne in mind that the resulting computed current
is vertically integrated over total depth.
Laevastu, T., 1974a: A multilayer hydrodynamical
numerical model (W. Hansen type] Model description and
operating/running instructions. Report to EPA, Part 2.
-------
2.9 TREATMENT OF TIDAL FLATS AND INUNDATIONS
The symbolic depth values at z points can be used for
a variety of indicative purposes. In the computation for
areas with tidal flats, the z points can be used for testing
whether the flats are dry (e.g., using indicator -5). At the
same time, indicators of the absence of water must be used in
depth fields for u and v points. This is done in actual depth
fields (HGU and HGV in the program) by assigning a fixed
negative number (e.g., -3) to points that are temporarily dry.
A test on the flow from one to another grid point must be
performed in the case of extensive tidal flats.
If in storm-surge computations the extent of inundation
of bars and low coastal areas is required, negative depth
(i.e., height) values must be coded over land at u and v as
well as at z points. In this case the -1 boundary values also
get actual (but negative) height values.
- 12 -
-------
3. TREATMENT OF OPEN BOUNDARIES
The HN model has been used with success for the confuta-
tion of currents in straits and semi-closed seas with two or
more entrances and along an open coast. In these models the
principal entrance is the entrance from which the major tides
move into the area. This principal entrance should preferably
be parallel to the x coordinate at the top of the grid (usually
along N=2 gri d line).
Tidal heights are also introduced at the secondary open
boundary (which in the case of straits is usually the last
line of grid points along the x axis). Care should be taken
that the tidal epochs along the second input line have a
proper phase lag from the principal tidal entrance. The proper
phase lag can usually be obtained through a tidal harmonic
analysis from a location near the second entrance. If the
phase lag is correctly computed, no discrepancies will occur
between the z values at the second open boundary and the
preceding grid point z values.
The main requirement in models with multiple open bounda-
ries, which are applicable to open coasts, is that "non-input"
boundaries (i.e., boundaries which are not prescribed in each
time step) be kept open for in- and out-flow as dictated by
currents and sea-level changes inside the area. The first
successful method for treatment of the "along the coast" model
was termed the "lazyman method." In this method, the tides
were introduced at a boundary perpendicular to the coast so
that the coast was located at the left side in the computa-
tional area. The tidal height was made inversely proportional
to the depth along the input boundary. There was no special
treatment of the seaward open boundary, except that the
computations of u, v and z were extended toward this boundary
- 13 -
-------
as far as the finite difference approaches permitted (thus
the "lazyman method") and the missing values of U at the right
hand boundary and of V at the lower boundary were assumed to
be the same as at the previous column (m-1) and row (n-1)
respectively. The derivation of the open boundary values of
z, u, and v by linear extrapolation (from the gradients between
second and third grid points inside the area to the boundary)
has not given fully satisfactory results.
- 14 -
-------
4. COMPUTATION OF DIFFUSION AND DISPERSAL OF POLLUTANTS
As the HN models provide a current vector at each grid
point at each time step, one is tempted to solve the inter-
acting advection and diffusion problem using a time-stepping
finite difference method. It should be borne in mind that
diffusion (mixing) is irreversible, whereas advection is
reversible. This model treats the dispersion of those sub-
stances which can be readily mixed in shallow water or in
a surface mixed layer. If these problems can be properly
solved by the HN model, no difficulties are foreseen in
treating specific problems such as dispersion of oil, brackish
water, etc. The following were considered in designing the
dispersion portion of the model:
a. Vertical diffusion should be neglected (can be
included in multi-layer models).
b. A basic Lagrangian approach should be used because
of its application to computer solutions for grid
arrays in finite-difference time-stepping.
c. Eulerian modifications should be used for instanta-
neous releases where the size of the blob is small
in comparison to the mesh length.
d. Fickian diffusion with constant diffusity is used
since advection and diffusion equations are solved
separately in short, interlocking time-steps
(similar to the "two-step method" of Schonfeld and
Groen, 1961).8
e. The Austausch coefficient for eddy diffusion per se
must be related to grid length.
Schonfeld, J. C., and P. Groen, 1961: Mixing and
exchange processes. Radioactive waste disposal into the sea.
IAEA Safety Series 5:1-00-1 32;
- 15 -
-------
Diffusion in water bodies has been presented in many
formulas, a few of which are presented here (neglecting
verti cal di ffusi on) .
The general diffusion formula is:
32S , a2S 1 8S _ n
2 2 ~ ~
8x 9y A at
The basic dispersion formula:
ft ¦ V" f - Ix (Ux+PSx> - ly 'Uy+pSy'
The Fickian equation:
H = Y" I " "2S " l*(SuU> - |y(SvV)
where:
Y - addition (release); n - decay; S - concentration;
pS - concentration velocity component; t- time; K
*»y 2 ?
(coefficient) = 3avr (a - depth; vr = u + v ; $ = 0.003);
Su>Sy - concentration gradients in u and v direction; and
A - diffusion coefficient (Austausch coefficient).
The Lagrangian approach of diffusion in finite difference
q
form was adopted by Wolff, Hansen and Joseph (1970):
St+T M st ^ . 4tA\ + At /st st st
n,m n,m I J 02 \ n-1 ,m n+1 ,m n,m-l
+ - 4S*"
n,m + l n,m
where
t is time; t - time step; Z - grid size; S - concentration;
A - diffusion coefficient; n and m are grid coordinates.
Wolff, P. M., W. Hansen and J. Joseph, 1970: Investi-
gation and prediction of dispersion of pollutants in the
sea with Hydrodynamical-Numerical (HN) model. Mari ne
Pollution and Sealife (Fishing News Ltd., London): 1T6-150.
- 16 -
-------
The above finite equation deviates from the usual finite
difference diffusion formula only in the addition of the last
term (-4S), which makes the solution similar to the solution
of the "Laplacian" (KV S). Secondly, the advection is com-
puted linearly in finite difference form:
ct + T = ct
n ,m n ,m ~ T
.t + T
n ,m
(Sn,m " Sn,m±l) ~ T
Pt + T
n ,m
(st - St ^
\ n,m n±l.m I
where: Sn,m-1 or n,m+l (respectively n-1, n+1) are used,
depending on the direction (sign) of U and V.
The Lagrangian approach used by Wolff, Hansen and Joseph
(1970), though reproducing the diffusion process well, does
not conserve absolutely the amount of the dispersing substance
The following modified formula was, however, found to be
conservative, provided the proper A (Austausch coefficient)
is chosen and corresponds to the chosen time step and grid
size:
rt
.t + T _
n ,m
= S
n ,m
+ S
4tA ct
„2 n,m
+ S
tA
7
n ,m-1
n+1,m n ,m+l
- 4S
fst , +
y n-1 ,m
t \ + lA /V
n,mJ ^2 ^ n-1,m-
+ Sn-1 ,m+l + S(n+1 ,m-l) + S(n+l,m+l)j
It should be mentioned that the transport equation used
in EPRF programs is very similar to the "upwind" difference
scheme used in air pollution problems (Pandol fo e_t a_l_, 1971)!
10
Pandolfo, T. P., M. A. Atwa.ter and G. E. Anderson, 1971
Prediction by numerical models of transport and diffusion in
an urban boundary layer" Cent. Env. and Man, Inc., Hartford.
Rpt. 4082: 139 pp.
- 17 -
-------
This scheme requires that:
The Austausch coefficient (A) is a function of grid
size and time-step. In an experiment designed to investigate
the conservation of diffusing substances, one of the main
criteria was found to be a relationship between grid size and
time-step. This relation of A to grid size with a 1,000 sec
time-step is shown on Figure 3. The correct value of A is
found from this graph and from the relation: A = 1 A
(graph). With a small time-step, the A approaches that found
empirically by Okubo and Ozmidov (1 970),1 1 and Kullenberg
12
(1972). An idea of the proper A to be used in different grid
1 3
sizes can be obtained also from the Joseph and Sendner (1958)
formulation. There is still some slight uncertainty about the
dependence of the horizontal Austausch coefficient (K^) on
the length scale:
-3 k 2
Kh = k] x 10 0
Kullenberg (1972) gives the values for the coefficients
k 1 = 1 .3 and k2 = "1. 31 ; whereas, Okubo (1971)14 gives the
corresponding values as 1.03 and 1.15. Both lines are shown
in Figure 3.
Okubo, A., and R. V. Ozmidov, 1970: Empirical
dependence of the coefficient of horizontal turbulent dif-
fusion in the ocean on the scale of the phenomenon in question.
Izv. Atmosph. and Ocean. Phys. 6 (5):534-536.
1 2
Kullenberg, G., 1972: Apparent horizontal diffusion
in stratified vertical shear flow. Tel 1 us, 24(1 ): 1 7-28.
1 3
Joseph, J., and H. Sendner, 1958: Uber die horizontale
Diffusion im Meere. Dtsch. Hydrogr. Zeitschr. 11(2): 49-77.
1 4
Okubo, A., 1971: Oceanic diffusion diagrams. Deep
Sea Res. 18, (8): 789-802.
- 18 -
-------
N
- 0RUB0 & OZMI
DOV (1970)
*
vsn
\
X \j
LENBERG (1972]
QKU8Q (19
«r
\V
X
At=10t
\s
lOsec.-^
..
1 o3
102 10 1
GRID SIZE , km.
0 . 1
Figure 3. Dependence of Austausch Coefficient (A) on
grid size (At=1000 sec)
- 19 -
-------
If the dimensions of a diffusion blob are small in
relation to grid size (in the case of an instantaneous
release), it is preferable to use a Eulerian approach in the
initial time stage. In this treatment, the center of the
blob is advected by the currents and its position is recorded
in every time step. The concentrations of neighboring grid
points as well as the change of central concentration are
computed with, for example, the Joseph-Sendner (1958) formula
or other available formulations. This approach is repeated
until a number of surrounding grid points have some concen-
trations above a predetermined minimum level, whereafter the
Lagrangian treatment is continued.
It should be noted that if the grid size is large (e.g.,
1 km or larger), the results of diffusion computations are
not fully satisfactory as yet (see Laevastu e^t aj_, 1974c and
1974d).15&16
Decay (uptake by biota, remineralization, etc.) was not
included in the tests reported above; however, its inclusion
would not be difficult. Furthermore, "thermal pollution" can
be treated by including the heat exchange terms in the model,
(see Laevastu, and Hubert, 1965). ^ 7
1 5
Laevastu, T., in collaboration with M. Clancy and A.
Stroud, 1974: Currents, tides, and dispersal of pollutants
in Lower Bay and Approaches to New York with fine and medi~u"m
grid size hydrodynamical-numerical models"! Report to EPA,
Part 3.
^Laevastu, T. and R. Callaway in collaboration with A.
Stroud and M. Clancy, 1974: Computation of tides, currents
and dispersal of pollutants in the New York Bight from Block
Island to Atlantic City with large grid size, single and two-
layer hydrodynamical-numerical models" Report to EPA, Part 4
^Laevastu, T., and W. E. Hubert, 1967: Analysis and
prediction of the depth of the thermocline and near-surface
thermal structure. Navy Weather Research Facility Rpt.
36-0667-128:90 pp.
- 20 -
-------
1 8
Rammig (1971) used a somewhat different approach for
the computation of diffusion. His basic formula is;
ii + u ii + V ii - A 3*s + R a2s
» + u»x+»„ Ai7+Bi7
assuming that A = B, his finite difference form (for a square
grid) for the u component is:
SS!I = sS.. " u^,m(Sn,m+1 - ZSJ., + Sn\m+lV 2A*V
sS,m-i - »S,.+ sJ,m+i)+ 2,«'t2 (sS-i.. -«S.. +-sS+i..)
The V component is analogous to the above. This formula is
also contained with blocked comments in column one in.the
program in Appendix C.
The stability criterion of this diffusion formula is:
*2
At < JK
Rammig (1971) derived a new Austausch coefficient at each
point, which is a function of depth and current speed:
A = 0.324 H |u|
The values for A derived with the above formula might be
appropriate in shallow, turbulent rivers, but tend to be
somewhat too low for coastal waters. However, further tests
on this formula with actual data are required.
1 8
Rammig, H. G., 1971: Hydrodynamisch-Numerische
Untersuchungen liber Stromungsabhangi ge Hori zontal ausbrei tungen
von Stoffen in Model 1 fl lissen mi t Anwendungen auf die Elbe.
Mitteil. Inst. Meeresk., Hamburg 14:64 pp and figures.
- 21 -
-------
In case of computing the movement of drifting oil or
other flotsam at the surface, it is necessary to give addi-
tional movement to this flotsam by wind (leeway), i.e., the
oil at the surface is pushed with wind. The present model
computes the leeway as a percentage of the wind. This
percentage must be introduced with a control card.
In the program (see Appendix C) the initial position of
oil or drifting objects is introduced with a control card.
The drift due to tidal currents and the additional leeway
caused by wind are computed in each time step. The position
is printed out each time output is requested (e.g., hourly).
Any number of initial positions can be included in the program.
- 22 -
-------
5. REVIEW OF THE PROGRAM AND OUTPUT FORMATS
The computer program i ri FORTRAN is reproduced in
Appendix C. The control statments are contained in the first
subroutine, called Control Program. This is followed by
Subroutine JO 1 which is for reading the input data (which
are printed-out for checki ng. purposes) and ini ti al i zati on. of
the fields. If no u, v and z values are read from previous
computations, these fields and the intermediates (UPB.VPB)
are initialized 0 over land, but over the water a low initial
value (0.2) is, read-in to facilitate the initial computations.
The interpolation of tidal input is done next, if two sets of
tidal constituents are used.
The main computations are done in J02, which is repeated
once every time step. It consists essentially of input of
tides, smoothing and averaging of u, v and z fields, their
qomputation, and setting of boundary values.
If wind fields are specified, the subroutine of wind current
computation is called in J02. If diffusion computations are
desired, the diffusion and transport field (S(N,M)) is
computed before the wind current component computation.
However, additional drift of oil floating on the surface, as
well as drifting object computation, is done after computing
the wind current components.
The program contains an option for the computation of
rest currents (e.g., after a full tidal cycle) and two sets
of diffusion computation formulas. The formula from Rammig
(1971) is blocked with comment statements in the appended
program, as is' the computation of wind currents with mixed
layer depth..
The horizontal fields of water level and current speed
and direction are printed out every one or two hours, or in
any other desired time interval in Subroutine J03. Further,
the data for sel ected'speci al. points are written on tape 2
- 23; -
-------
every time field output is requested. The data at special
points are printed out at the very end of the computation
(J05).
A plotting program (not reproduced here) has been added
to the HN model to allow the plotting of results on the
1 9
Calcomp plotter (Bauer, 1969). A series of statements are
added to the basic program to capture the arrays of data
that are required to plot the water level and velocity
vectors. The plotting program is adequately documented by
Bauer (1969) and is not described further in this document.
It should be pointed out that the model must run six
to ten hours of real time (sometimes more, depending on the
size of the area and grid length) before equilibrium is
established and correct outputs can be taken.
The computations toward the coast (i.e., coast located
to the right or in the lower part of the grid field) have
not been successful yet. This is mainly due to the necessity
of complex treatment of singularities (boundary problems) at
the coast, which is not fully included in the present program.
19
Bauer, R. A., 1969: Modification to program Hansen.
Comp. Appl. , Inc., San Diego, MS Rpt.
- 24 -
-------
REPORTS IN THIS SERIES
A vertically integrated hydrodynamical-numerical
model (W. Hansen type): Model description and
operating/running instructions. (This report)
A multi-layer hydrodynamical-numerical model
(W. Hansen type): Model description and operating/
running instructions.
Computation of tides, currents and dispersal of
pollutants in Lower Bay and Approaches to New York
with fine and medium grid size hydrodynamical-
numerical models.
Computation of tides, currents, and dispersal of
pollutants in the New York Bight from Block Island
to Atlantic City with large grid size, single and
two-layer hydrodynamical-numerical models.
- 25 -
-------
APPENDIXES
APPENDIX A - INPUT PARAMETERS, CARD FORMATS,
AND SET-UP OF TAPES
APPENDIX B - ABBREVIATIONS AND PARAMETERS
USED IN THE PROGRAM
APPENDIX C - PROGRAM NYBIGHT LISTING
- 27 -
-------
APPENDIX A
INPUT PARAMETERS, FORMATS AND SET-UP OF TAPES
INPUT PARAMETERS AND DATA CARDS
Card 1 Format 3A8
TIT
TIT Title of the computation
Card 2
JA, NE, ME, KO, IUE, KKE, NURU, KI, LE, LU, LI, LOP, LUN, LA,
LIV, LEET
JA Indicator, whether computation starts from initially
0 fields (JA = 0: Z = 0, U = 0, V = 0) or from the
values of Z, U and V from tape).
NOTE: If JA = 1, T on card 4 could be specified as
the time of the computation of the input values or
any other adjusted "real" time.
NE,ME Delimeters of entire field (size of the computa-
ti onal gri d).
KO Number of points on the open boundary where tidal
constituents are introduced.
IUE Twice the number of wind fields.
KKE Number of wind field characteristics (A(I)) ' s
(=(IUE/2)*3)
NURU Number of selected special points
KI Indicator, whether one or two sets of harmonic
constants are used. In the latter case, the
constants are interpolated linearly to each input
point. KI = 0 - one set of constants; KI = 1 - two
sets of constants.
- 28 -
-------
LE Indicator, whether rest current computation over
i.
one tidal cycle (or any specified time interval)
is desired. LE = 0 - no rest current computation;
LE = 1 - rest current computation (specify beginning
and end times on card 4).
LU Indicator and number of input points of permanent
current. LU = 0 - no permanent current input;
LU = 12 - permanent current input at 12 points
(specify the N and M positions of input (NL, ML,
card 21)) .
LI Indicator for computation of diffusion and transport.
LI = 0 - no diffusion computation; LI = 1 - diffu-
sion computation.
LOP Indicator for plotting program; LOP = 0 - no plotting;
LOP = 1 - plotting required (requires additional
tape).
LUN Indicator for depth input. LUN = 0 - depths at U
and V point different and read from the cards;
LUN = 1 - depth for U and V points are identical
(changes have to be made at the boundaries).
LA Number of points on second open boundary (NG and
MG - i.e. , the N and M of the points must be
speci fi ed on card 17).
LIV Indicator for drift and oil drift calculations.
LIV = 0 - no drift calculations; LIV = 1 - object
drift calculation; LIV = 2 - oil drift calculation.
LEET Indicator, whether data to be put on tape and read
from it.
Card 3 Format 9F8.3
G, ALPHA, RBETA, CI , TIL
p
G Acceleration of gravity . (e.g., 978.0 cm sec )
ALPHA Smoothing parameter (e.g., 0.998).
RBETA Coefficient of geostrophic wind (e.g., 0.75)
- 29 -
-------
CI Wind indicator; 0 - no wind; 1 or 2 - wind (Specify
beginning and end times on card 4). The number
indicates the number of wind field changes.
TIL Angle between true N and 90° of computational grid
( + V).
Card 4 Format 9F8.0
DT, TE, TW, T1 , 12, SI, T, TRB, TRE
DT 1/2 time step (sec) (e.g., 30).
TE Length of computation (sec) (e.g., 180000)
TW Time when wind starts (sec) (e.g., 7200)
T1 Time interval between printouts (e.g., 3600 sec)
T2 Field output counter (0 if outputs desired from
the start of the computation; otherwise any other
delayed starting time, e.g., 72000 sec).
SI Time when wind stops (sec) (e.g., 180000)
T Time (initialized 0. If previous computations
made and z, u, v, read from tape, give the time
previous computation ended).
TRB Starting time for rest current computation. Used
also for starting time for drift computations.
TRE Ending time for rest current computation. Used
also for ending time for drift computations.
Card 5 Format 9F8.0
T4
T4 Time when data on Z, U and V to be put on tape.
Card 6 Format
DL, F, SIGMA, R, R0L, C
DL 1/2 space step (half
F Coriolis parameter (
6E12.4
grid size in cm, e.g., 185000)
.g . , 8.55 x 10"5)
30 -
-------
SIGMA
R
ROL
C
Card 7
AUS
AUS
Card 8
NK, MK
NK
MK
Coordinates of the points at the first.open boundary (input
boundary).
Card 9 (one or more cards Format 24.13
NZ, MZ
NZ N coordinate
MZ M coordinate
Coordinates of special selected points- where special outputs
are desired.
Card 10 Format 2413
NU, MU
NU N coordinate
MU M coordinate
Coordinates of the upper right and lower left corners of
the wind field.
Angular velocity, of M,.fide (1.4088 x 10~4^r (flpt^
used in enclosed program).
Fri cVi oW cofefcfix5! e^rft* t)03-)rl^ ^ 0 29fnB^ ^^
Density of the air (e.g., 1.1627 x 10"^) (not
used in enclosed3 pVo'gVam).
Drag coefficient (e.g., 3.2 x 10"®) (I)A
Format 4E10.2
Austausch coefficient (e.g., 1.0E5).
Format 2413
N coordinate
M coordinate
- 31 -
-------
Card 11
V(I)
Format 10A5
Names of the selected special points.
Card 12 Format 9F8.2
A (I)
A(I) Interval (in sec) between wind speed and direction
change. If only one wind speed is used, then
duration of the wind (SI - TW) (For the program
in Appendix C the time intervals between the wind
changes must be equal).
A(2) Wind speed (m sec"^)
A(3) Wind direction (in degrees in computation coordintes).
Card 13 Format 14F5.1
H
H Amplitudes of tidal constituents
Card 14
B
B Amplitudes of tidal constituents at second boundary.
Used only if two sets of tidal constituents are
used and interpolated, or tidal input is made at a
second boundary.
Card 15 Format 7F10.9
AL
AL Speeds of the tidal constituents (a) in radians
per second.
Card 16 Format 7F10.9
CAPA
CAPA Phase angles (k) of the tidal constituents.
- 32 -
-------
Card 17
CAPB
CAPB Phase angles (k) of the tidal constituents. Used
only if two sets of tidal constituents are used
and interpolated, or tidal input is made at a
second boundary.
Card ,18 (one or more cards) Format 2413
NG, MG
NG N coordinate
MG M coordinate
Coordinates of the points on the second open boundary.
Used only if two or more open boundaries are used in the area
Card 19 (several cards) Format 12F6.0
HTZ
Symbolic depths at the z points
Card 20 (several cards) Format 12F6.0
HTU • :
Depths at u points (cm)
Card 21 (several cards) Format 12F6.0
HTV
Depths at v points (cm)
Card 22 Format 2413
NL, ML
N and M coordinates of the input for permanent current.
Used only if permanent current input is included.
Card 23 (one or more) Format 14-F5;-. 1;
PERU
PERU U component of permanent current. Used only if
permanent current input is included.
- 33 -
-------
Card 24 (one or more) Format 14F5.1
PERV
PERV V component of permanent current
Card 25
POSN, POSM, DRI
POSN N coordinate of the
object.
POSM M coordinate of the
object.
DRI Percentage of wind
(1eeway)
Format 9F8.2
initial position of drifting
initial position of drifting
peed for additional drift
Card 26, 27 Plotting parameters. See special publication.
Used only if plotting program is included.
SET-UP OF TAPES
The input parameters are read into the computer after
the program, usually from cards (or from tape). Input tape
50 on FORTRAN 63 CO-OP monitor indicates card reader. For
running the program on the CDC 6500, no tapes are required
if no plotting is desired. The program and data input cards
are read with the card reader. A disk is assigned to tape 2
(intermediate tape for writing special point outputs).
Outputs go directly to the printer. In case outputs are
required for plotting, one tape is required at the end of
the program.
For running of the program on the CDC 1604 or other
smaller computer (with FORTRAN 60 or FORTRAN 63), the
following tape units are required:
Tape 1 is the FORTRAN system tape;
Tape 2 is the intermediate outputs of special points;
- 34 -
-------
Tape 3 can be for input of the program and data if no
card reader is used and may be used also for outputs for
piotting purposes ;
Tape 4 can be used as the printer if no printer is
available;
Tape 5 can be used for input and output of intermediate
data on u, v, and z.
- 35 -
-------
APPENDIX B
ABBREVIATIONS AND PARAMETERS USED IN THE PROGRAM
Note: Symbols marked with an * are further
explained in Appendix A
A(I)* Characteristics of the wind field
AH(I,J) Amplitudes of tidal constituents at individual
input points
AINC(I) Intermediate parameter in interpolation of
tidal constituents
AL(I)* Speed of tidal constituents
ALPHA* Smoothing parameter
AN G(N , M) Current direction
ARG(I,J) (at - k) Argument of tidal constituent
AUS* Austausch coefficient
A1 to A5 Intermediate parameters (see statement 60-65)
B(I)* Amplitudes of second set of tidal'constituents
BETA 1 - ALPHA
B1 to B5 Intermediate parameters (see statements 412-415)
C* Drag coefficient (3.2 x 10"^)
CAINC ) Intermediate parameters in interpolation of
CAINCU J tidal constituents
C A PA(I)* k of first set of tidal constituents
CAPAP(I,J) k of tidal constituents at individual input
points
CAPB(I)* k of second set of tidal constituents
CNCU Intermediate parameter in interpolation of
tidal constituents
CW A counter
CI* Wind indicator
C3 Intermediate parameter
DL* Half grid size
DRI* Leeway (drift speed in % of wind speed)
- 36 -
-------
DT* Half time step
ETC* |
/ Plotting parameters
ETH* j
H(I)* Amplitudes of first set of tidal constituents
HGU(N ,M) \
/ Actual depths at u and v points
HGV(N ,M))
HT Counter (in plotting)
HTU(N,M)* |
} Water depths at u and v points
HTV(N,M)*)
HTZ(N.M)* Symbolic depths at z points
HMLD Mixed layer depth
I A counter, delimeter
IU Number of wind fields
IUE* Twice the number of wind fields
J A counter
JA* Initiation indicator
K A counter
KI* Tidal input indicator
KK.E* Number of wind field characteristics
KO* Number of points at the open boundary
L A counter
LA* Number of points on second input boundary
LE* Indicator for rest current computations
LEET* Indicator for writing intermediate output
on the tape
LI* Indicator for computation of diffusion and
transport
LIV* Indicator for drift and oil drift calculations
LOP* Indicator for plotting program
LU* Indicator for permanent current input
LUN* Indicator for simplified depth input
M* Grid index (x)
MD* Coordinate for drift computation
- 37 -
-------
ME* Delimeter of grid size (x)
MEH ME-1
M6* M coordinate of input points at second open
b o u n d a ry
MK* M coordinate of input points at first open
boundary
ML* M coordinates of input points of permanent
current
MU* Wind field delimeter
MUM A counter
MZ* M coordinate for output at special points
Ml , M2 Indi ces
N* Gri d i ndex (y)
ND* Coordinate for drift computation
NE* Delimeter of grid size (y)
NEH NE-1
NG* N coordinates of input points at second open
boundary
NK* N coordinates of input points at first open
boundary
NL* N coordinates of input points of permanent
current
NON A counter
NU* Wind field delimeter
NURU* Number of special output points
NZ* N coordinate of output at special points
N1.N2 Indices
OTC(I) Intermediate for meteorological overtides
input at the coast
PERU(I)* )
> u and v components of permanent current
PERV(I)* )
POSM* )
/ M and N coordinates of the drifting object
POSN* )
- 38 -
-------
PREC(I ,J)
PROD
Intermediate for interpolation of tidal input
Percentage of wind speed for additional drift
of oil
Friction coefficient
Current speed (resultant)
Coefficient of geostrophic wind (0.75)
Density of the air ( 1 .1627xl0~3)
Intermediate (square root)
Concentration of diffusing substance
Time when wind stops
Angular velocity of M2 tide
Intermediates in diffusion computation:
U(N,M)
UDIS
Plotting parameters
Time counter
Intermediate for input of meteorological.
overtides at the coast
Length of computation
A counter
Tilt of the grid
Plotting parameters
Title
Starting time for .rest current computation (sec)
Ending time for rest current computation (sec)
Time when wind starts (sec)
Time interval between printouts
Field output counter
Time when intermediate u, v and z values to be
put on tape
U component of current
intermediate for drift computation
- 39 -
-------
UML
Intermediates for oil drift computation
Intermediate field for smoothing and averaging
U component of rest current
V component of current
Intermediates for smoothing
Intermediate for drift computation
Plotting parameter
Intermediates for averaging
Intermediates for oil drift computation
Intermediate field for smoothing and averaging
Intermediate for averaging
V component of rest current
Intermediates for averaging
Names of special output points
Intermediates for output at special points
U and v components of wind current
Decay of pollutant
Sea level
- 40 -
-------
APPENDIX C
PROGRAM NYBIGHT LISTING
PROGRAMNYBIGHT(INPUT,OUTPUT,TAPE50SINPUT,TaPE2«TaPe3»TaPE5)
DIMENSION HjZ(42.4742.47).HTV<42,47>,Z{42,47),U<«2,<7),V<42,
i47)«US(42«47)»VS(42*47),RAD(42'47)#ANG(42*47),XK(42»47),YK(42*47)i
2 NUl6),MU(6).NK(B'5).MK(85)« AH<5,85),CApAp<5,85>,PREC(5,85), ArG
3(5.85).PERU(65).PERV(85),NL(B5)»ML(85),N2(lp),M2(lo>»VKiO),UPB{42
4,47).VPB<42.47),st42,47>,A(85 J , V2 < 10), v3 < 10 >, V4 UO ), HGu < 42.4
57),HGV(42,47)#H(5).B(5)«AL(5>«CAPA(5)«CAPB<5)jAINC<5>.CAINC<5 J,
6 NG(40).MG(40),TlT(3).AUS(5),TAR(85),0TCt85)
7,USR(42*47),VSR(42,47)
COMMON ''TZ , MTU, HJV, Z # U, V, US, VS, RAD,.ANG, XK, YK,. NU# MU,HK,MK ,-A H, CAP A
1P,PREC,ARG, PERU,PERV>L,.ML,NZ,MZ., Vi, UPB,VPB,S,A , V2IV3,V4-,HGU.HG
2V/H, B, aL,CAPA » CAPS, A INC,CA INC,NG, MS, TIT,AUS,TAR, OTC»USR.VSR,
3T1.T2»T4.TE,TW,A1.A2,A3.A4,DT,BETA,ALPHA»TIR,DL.G,C»C3,^E,ME,NEH.M
4EH,JA, IUE, KKE, NURU, SI. T, RBETA,Cl, SIGMA, ROL, A5#RU, RVi KO* LI .'LE, LU.Kl
5,LOP,TRB,TRE,TIC, LUN» i>.TIL.LIV.LEET
6, POSN#POSM,DRI, STV, ETV, TINCV, STH, ETC, TINCH,CNCU,CM NCU
7,AA1,BB2,CC3,TSC,DL2
8,AA2,AA3,AA4,BB3,BB4'#BB5,9B6,0B7,BB8,CC4,CC5,CC6,Cc7.CC8
C CONTROLU PROGRAM
4 REWIND 2
5 CALL JOBl
1007 CONTINUE
6 CALL J0B3
1009 continue
1010 IF(LOP) 7.7,1011
1011 CALL JOB7
IF(T»TE)7.240» 24 0
240 IF(LOP)1016,1016,241
241 1F(LOP-B)1013.1016.1016
7 T2=T2+T1
8 TsT+Al
9 JO02 Control Program
10 IFCT-TEJ11.1012,1012 _ JOB6 and JOB7 are
11 lF(ABsF-Al/2,>6.6,12 for a special plot
12 GO TO 8 program, as is Tape 3
1012 IF(LOP)13(13,l0l3
1013 L0P=5
13 CALL JOB3
1014 IF(LOP)1016,1016,1015
1015 CALL JOB7
1016 END FILE 2
14 REWIND 2
15 CALL J0C5
1017 IF(LOP)*6,16,245
245 END FILE 3
- 41 -
-------
REWIND 3
1016 CONTINUE
16 STOP
END
SUBROUTINE JOB1
DIMENSION HT2(42,47),HtU(42,47),HTV<42,47),2(42,4?)iU<42,47),V<42,
147),US<42,47)»VS(42,47),RAD(42,47),ANG(42,47),XKt42,47),YK{42,47),
2 NU<6),MU(6),NK(85),MK(85),AH(5,85),CAPAP(5,85)»PREC(5,85).ARG
3(5,85) ,PERU(85),PERV(85>,NL<85>,ML<85),NZ(10),MZ(lO> »V1<10),UPB<42
4,47),VPB(42,47).S<42,47>,A(85) , V2( 10 ), V3(10 ), V4UO ), HGU( 42,4
57),HGV(^2.47),H(5),B(5)*AL(5),CAPA(5)»CAPB(5)»AINC(5)ICAINC(5).
6 NG(40),MG(40),Tlt(3),AUS(5),TAR(85),0TC(85)
7,USR<42,47),VSR(42,47)
COMMON HTZ.HTU,hTV,Z,U.V,US,V9iRAD,ANG,XK,YK, NUiMU,NK,MK, Ah, CAP A
IP.PREC.ARG,PERU.PERV.nl*ML,NZ,HZ,VI.UPB.VPB.S,A ,V2,V3,V4,HGU#HG
2v,h,b,aI!,capa.capbi ainc,cainc.ng,mg,tit,aus,tar,otc#usr,vsr,
3T1,T2#T4,TE,TW,A1jA2, A3,A4,DT,BETA,ALPHA,F",R,DL,GiC#C3,NE,ME, NEH,M
4EH, JA,.10E,KKE,NURU,S! .f,RBETA.Cl,SlGMA,R0L,A5,RU,«V»K0,LI ,LE,LU,K!
5,lop.trb,tre,tic, lun,l'a.til,liv,leet
6,P0SN.P0SM.DRl,STV,ETV.TINCV,5TH,ETC.TlNCH,CNCU.CAINCU
7,AA1.BB2,CC3,TSC,DL2
8,AA2,AA3.AA4,BB3,BB4#BB5,BB6(BB7,BB8,CC4,CC5,CC6iCc7#CC8
: reading or values and initiation
ETH » ETC
29 CoNTlNU*
30 FORMAT(2413)
31 F0RMAT(9F8.3)
32 FORMAT(9F8.0)
33 F0RMAT(6E12.4)
34 FORMAT(l5A5)
35 FQRMAT(9F8,2)
36 FORMAT(14F5.1)
37 FqrmAT(7F10.9)
38 FORMAT(12F6.0)
39 FORMAT(3A8)
55 FORMAT(ElO.2)
56 READ (50,39)(TIT(I).Uli3)
40 READ (50.30)JA.NE.ME,KO,IUE.KKE.NURU,KI,
1,LA,LUN,LIV,LEET
41 READ <50,31)G,ALPHA.RBETA.Cl,TIL-
42 READ (50,32)DT,TE.TW,T1,T2,SI,T.TRB.TRE
1006 IF(LEET) 43,43,1008
1008 READ <50,32)T4
43 READ <50,33)DL,F,SIGMA,R,ROL,C
57 READ (50,55)AUS(1)
44 READ (50.30)(NK(I),MK
-------
47
1030
1036
48
1032
1033
1034
1037
49
58
59
50
51
69
70
1020
1021
52
68
1022
53
READ (50,34)(VI(I),Isl,NURU>
NQNel
Cvi = 1,
READ (50*35)(A(i).I«NON«KKE)
CWoCwlT
NONckkE+l
kke=kke*non»i
IF(Ci-CW) 49,48.48
READ (50,36){H(l).Iol.5)
IF(KI>9y «50159
READ (50,36)
READ (50*37)(AL(I)«I^1«5)
READ (50,37>
IT(KI) 1020,1020,70
READ (50,37)(CAP0(I),1=1,5)
IF(IA)5?,52,102l . .
READ (50,30)(NG(I)(HG.M»i,ME>»NI«liNE)
MEH=ME-1
IFtLUN) 53,53,54
READ (50,38)((HtU•1al«LU)
READ (50,36) (PERV( 1 ) ,Ul.LU)
IF(JA) 1040,1040,93
READ (9?31)((Z(N,M),MBi,ME),NBl.NE)
READ (5,31)((U(N,M),M»l,ME),Nil.NE)
READ (5,31)((V(N«M),Hb1,ME:),Nb1*NE)
READ (5,32)t
1F(LIV)1005,1005,1041
IF(L1V-?) 1042,1005,1042
READ (50,35)POSN,POSM,DRI
lF{|_OP)&0i60«200
CAUL J0B6 .
DO 9967 N s 1,41
HTZCN,2i?»2.
HTZ«0.
HTV(30,17)a.l
DO 2031 N = 1»NE
DO 2031 M = 1,ME
IFtHTUtfcUM) ,QT,0 ) HTUtM«M)
IF(HTV(N,M) ,QT,0). HTV(NiM)
IF(HTU(N«M),IE«0 ,) GO tO 2033
IF(HTU(N.H)-200.>2032,2032,2033
Reading of
input data
}i
Calling of initial plotting
program; setting an input
boundary
B HTu(N.M) *
» HTV(NiM) *
185,201
185.201
Converting of
depths from
- fathoms to cm
and setting of
1 minimum depth.
-------
2032
2033
2034
2031
60
61
62
63
64
65
66
67
i 6.32
99
HTU(N,M> • 200,
IF(HTV(N«M),LE,0.) GO TO 2031
IF-200.> 20 34,2034,2C3l
HTV(NiM) e 200,
CONTINUE, .
B"ETA=(l7-ALPHA)/4.
A1s2.*DT
A2«F*a1
A3sR*Al
A4*DT/DL
A5sG*A4
C3iC*Al»10 000 •
NEH=NE-l
KKE a (IUE/2) * 3
TIC = 0 .
J
— Setting of some intermediate
parameters.
77
75
76
78
TSC
AA1
AA2
AA3
AA4
BB2
SB3
BB4
8B5
8B6
BB7
BBfl
CC3
CC4
CCS
CC6
CC7
CC6
o;
oT
o';
o'.'
o:
o:
o:
o;
o:
o;
o;
o;
o:
oi
o:
o:
o"
— _ - 0» —
DO 100 N«l,Nl
DO 100 M«1,ME
iF(HTZ(N,M))7BI78.80
Z(N.H)«0.
— Zeroing of transport arrays
VPBIN# H)«0.
79 GO TO 81
60 Z) 82 # 82»84
82 86,86.8-8
66 V(N'M ) *0 «
VPB(N»M)*0 .
87 GO TO 95
88 V(N» M)=0 ,2
Initialization of
arrays
- 44 -
-------
VPBCN,M>«0.2
95 XK(N,M>?0.
96 YK(N,M>=0,
S(N»M)=0 ,
US(N#H) b o.
VS(N.M) b 0.
USR(N,M)«0.
VSR(N,M)aO.
100 CONTINUE
1060 1F(LUN)101,1Q1«1061
1061 00 1064 N=1,NE — Special setting of
1062 DO 1064 Mai,ME depths if only one
1063 HtU(N,M)sHtv(N,M) set was read in.
1064 CONTINUE
101 IFCKj >1- 0.120,102 —-i
102 DO 107 J81.5
103 AINC(J)e(H(i>-B(I))/(K0-l)
1101 IF(CAPA J) = CAPAU) • CAINCU
1108 IF(CAPAP(I.J)-360.>1111.1111.110 9
1109 CAPAP (,I:, J) =CAPAP( I.J).360,
1110 G0 TO 113
1111 IF (CAPAP ( I# J) >1112#1112#113
1112 CAPAP(I.J)=CAPAP(I,J)*360,
113 CAINCU = CAINCU ~ CAlNC(I)
114 CONTINUE
115 QO TO 125
120 DO 124 'Jsi.KO;
121 DO 124 Is1,5
122 AH( I. J.) (;!')¦'
123 CAPAP(I,J)bCAPA(I)
124 CONTINUE
125 CONTINUE
1130 PRINT l20#(TIT(H,Iol.3)
106 DO 114 1 s 1.5
1065 CNCU=0.
1066 CAJnCusO,
109 DO 114 J s i#KO
110 i , b h< I) - CNCu.
111 CNCU = CNCU ~ AINC(I)
Interpolation of
tidal input para-
meters.
(If required)
- 45 -
-------
— Printing of input
1131 PRINT 181,DT,DL
1132 PRINT 182.F.AUSU)
1133 PRINT 1*3,ALPHA,Cl
1134 PRINT 184.TE.T1,TW.SI
1114 PRINT 1141.T2.TRB.tRE
1115 PRINT 1142.TIL
1116 PRINT 1143.KI.LE.LU.LI.L0P.LUN.LA.LIV
1117 IF{LIV)1120.1120.1118
1118 IF(L!V-2>lil9,1120.1119
1119 PRINT 1144, P0SN.P0SM.DRI
1120 CONTINUE
1135 PRINT 185. (NK( D.MKd). I»i.K0J *
1136 PRINT l86,H(l),H(2).H(3).He>).B(l),B(2).B(3)iB(4) for checking.
1137 PRINT 187,AL(1).AL(2).AL(3).Al(4) |
1138 PRINT 188,CAPA(1).CAPA(2>JCAPA(3).CAPA(4).CAPB(l),CAPB(2j,CApB{3).
lCAPB(4)
1140 MUM s xF IXF(Cl*3,)
1139 PRINT 189.(A(I),I«1.MUM)
1121 IF(L0P)1124.1124,1122
1122 PRINT 1145.STV.ETV.TINCV
1123 PRINT 1145,STH,ETH,T INCH
1124 CONTINUE
1125 IF(LU)1128,1128,1126
1126 PRINT 1*46, (NL(I ) .ML(M« I«1.LU)
1127 PRINT 1147. (PERU(I>.PERV(I>*!*1#LU>
1128 CONTINUE
180 F0RMAT(?H1,35x,3A8///)
181 FORMAT(5X,19HHALF TIME STEP DT c,F5.Q,5H
IE DL n ;E11,4.4H CM//>
182 FORMAT(5X.21HCORIOLIS PARAMETER .»,El2.4.5X,16HAUSTaUSCH COEF.b.
1E10,2//)
183 FORMAT(5X.27HSMOOTHING PARAMETER ALPHA = ,F6.3.20Xil'HWIINDICATO
1R Cl =»F8,3//)
184 F0RMAT(5X.13HC0MP. LENGTH .F8.0,5X.15HOUTPUT INTE*Vt .FB.0.5X.15HW
UN STARTS AT »F8.0.5X.11HWIND STOPS . F8, 0, 6HSEC0-MS* U)
185 FORMAT(5X»49H(N.M) COORDINATES OF INPUT POINTS AT THE BOUNDARY ,//
l5x.l2(lH(.I2.1H,.I2.1H).2x>//5X,12(1H(.I2.1H..12.Ih>#?X>//5X.5<1h<
2.I2.1H.*12.1H).2X)//)
186 F0RMAT(5X.29HAMPLITUDES OF TIDE CONSITUTES,5Xi8F10,1//)
187 FORMAT(5X.20HSEEDS OF. TIDE CONST,.5X.8F10.7.//)
188 F0RMAT(5X,2lHKAPpAS OF TIDE CONST,,5X.8Fl0,7,//>
189 FORMAT(5X.25HWIND.INTERV/.SPEED,DIR. .9F8,2i//5Xti2F8,2.7/)
1141 FORMAT(5X.25HWIND BEG..DRIFT BEG,.ENDS.3X.3F8.2.//)
1142 FORMAT(rX.l2HTlLT OF GR1D.5X.F4,Q.//)
1143 F0RMAT{5x.3HKIb.I3.3HLE*.I3.3WLUS.I3.3HLI",I3.4HL3P«#I3.4HLUN*.13.
l3HLAs,I3.4HLIVS.13,//)
1162 F0RMAT(///8X.15I8)
1163 F0RMAT(1X.I8.15F8.0)
1144 F0RMAT(5x'.25HSTARTING POS. OF Dr I FT , N>. F6 ,2. 2X, 2H*« # F 6 ,2 . 18HPERCEN
1TAGE OF WIND.F6.2.//)
1145 FORMAT((10X.15HPLOTTING PARAM,» 3F10,0/))
1146 FORMAT(5X.24HC00RD,OF PERM.CUPR.INPUT./5X.12(lH(.12'1H.iI2.1H).2X)
SEC.29X,2"HHALF G»ID S12
- 46 -
-------
l/5X>12(lH(.12,lH,f!2.1H),2X)/)
1147 FORMAT(EXi26HU AND V COMP.OF PERM,CURR,,/5X,8(1M(#F4.1#IH,,F4.1»lH
[lH(»F4.1»lW,iF4fl,iH)#2X)/5x#8(lH(iF4,lilH|,F4,l#lH)»2X
1)»2X)/5Xi8
2)/)
140 PRINT 160,
PRINT 1*1.
PR|NT 162,
PRINT 1^3i
PRlNtll62,
PRINTU63,
PRINf 164.
PRINT 161,
PRINT 162,
PRINT 163,
PRjNTll62:,
PRINT1163.
PRINT 165,
PRINT 161,
PRINT 1*2,
PRINT 163,
PRINT11&2.
PRINT1163,
DO 900 Nsl
DO 900 M«1
IF(HTU(N,M
141
142
143
144
145
146
147
148
149
150
151
N,N«1,16)
N,(HTZ
N,Nel7'.32>
N»(HtZ(N.M),M»i7.32),N«l,NE)
N,N«33,47)
N,.Nsl,NE)
N,Nsl,16)
N,(HTU(N,M),Mb1,16),Mb1,nE)
N, N=17,32)
N,(HTUfN,M)»Msl7,32),N=l,NE)
N,Ns33,47)
N,(HTU(N,M),Ms33,47>,N=i,NE)
N,N«l,l6)
N,al50000.
900 CONTINUE
160 FORMAT(5H1,16X,36HSYMB0LIC WATER
161 FORMATCI8,16F8.0)
162 F0RMAK///8X,1618)
163 FORMAf(l8,l6F8,0)
164 FORMAT(1H1#16X,32HWATER DEPTH
165 FqRMAT(1H1.16X,32HWATER DEPTH
170 DO 175 NbI.NE
171 DO 175 Mil,ME
173 HGU
174 HGV(N.M^HTV(N,M)
175 CONTINUE
1175 CONTINUE
RETURN
END
SUBROUTINE J0B2
DIMENSION HTZ
147),US<42,47>
2 NU'6),MU<
Setting of
maximum depth
DEPTH AT THE Z P01N'TSV/BX91618/1
THE U POINTS (CM)//8X,16I8/)
THE V POINTS (CM)//8X,16I8/)
Initial setting of variable depth
(42,47),HTU(42,47),HTV(42.47),Z(42,47),U(42,47),V(42,
,V?(42.47),RADM2»47),ANG(42i47),XK(42#47),VK(42.47).
< 6)»NK<83)1MK(85)* AM(5» 85),CAPAP(5|85)iPREC(5# 85), ARG
3(5,85),PERU(85),PERV(85)»NL(85),ML(85)#NZ(10)1MZ(l0 ^»VI(10),UPB'
4,47),VPB(42,47),S(42,47),A(85) ,V2(10),V3(io),V4<10),HGU<4
57),HGV(fl2'.47),H(5).B(5),AU(5).CAPA(5),CAPB(5),AINC(5),CAlNC(5).
„ i.A.iAi' kiAiiA; AMr/Ci VAO/oet HTPjQct
UPB(42
42,4
6 NG f 40 >.« HG c 40 > # 11 f (3>,AUS(5)iTAR(85),0TC(85)
7,USR(42i47),VSR(42,47)
- 47 -
-------
+ T
COMMON *TZ.HTU,HTV,Z,"U,V,'JS*VS*RAD,ANG,XK,YK( NU»MU#NK,MK,AH,CAPA
1P»PREC,ARG»PERU»PERV,NL*ML*NZ,MZ,V1,UPB,VPB*S,A *V2,V3,V4»HGUtHR
2V»H#B»AL*CAPA»CAPB,AINC'*CAINC,NG,MGiTIT*AUS|TARjQTc*USR»VSR,
3T1*T2»T4,TE»TW,A1»A2,A3jA4*DT,BETA*ALPHA*F,R,DL*G#C#C3*NE*ME,NEH,M
4EH, JA,IUE,KKE,NURU,S!,t,RBETA,Ci;SlGMA,R0L,A5,RU,RV#K0,Ll,LE.LU.K'!
5* LOP* TRB* TRE* TIC* LUM,LA,TlLiLIV,LEET
6*P0SN,P0SMj[1RI,STVjETV,TINCV,STh*ETC#TINCHiCNCU*CAINCU
7.AAl,BB2,CC3iTSC*DL2
8#AA2,AA3,AA4,BB3,BB4,BB5#BB6,BB7,8BB#CC4,CC5*CC6,Cc7,CC8
COMPUTATIONS
309 CONT'INUt
TD = 104400.
TG«TD-6100.
tgt=tg-td
Do 384 N«2,40
RN»N
Ds(RN-27)/39,
PY*FLOAT(NJ
ZtN*2)sH{l)*C0S(AL(l)*(TD*D*TQT)-CApA(l))*H{2>»C0S(At(2)*(TD+0»TGT
1)-CAPA(2) )+H(3)*C0S( AKSJttTD^D^TGTJ^CAPAtS) )*H(4)#C0S( AL<4)*(TD*D
2*TGT)-CAPA(4>>+H(5)*COS(AU(5)*< tD*D*TGT)-CAPA<5))
Z(N,2 ) s 0.86*Z(N» 2)
— Special boundary setting
Tidal input
384 CONTINUE
~O 386 Ni2,40
PZ*FLOATlN>
Z(N*2)'«E(N,2)-4./PZ
366 CONTINUE
DO 1215 M«2»44
PZ«FLOAT(M.)
Z(2*M)*£(2.M)-0.1*PZ
ZCl.M)a2{2,M)-0.08*P7
1215 CONTINUE —
390 DO 394 -NPl.NE
391 DO 394 Mil.ME
392 UPB(N»M)*0.
393 VPB(N*M)b0.
394 CONTINUE
500 DO 523 N.l.NE
501 DO 523 M»1.ME
IFtHTZtN.M)>523,523*503
IF(1-N)504.507,504
IF(HTZ(N-1*M>5 50 5,50 7,50 5
VAUP=Z(^-i.M>
506 GO TO 508
507 VAUP*Z(N.M>
IF(NE-N)509,512,509
IF(HTZ(N+l• M) )510.512,51'0
510 VALO»Z(N*l*M)
511 GO TO 513
312 VALO=Z(N*M>
513 JF(1-M)614,*16,514
502
503
504
505
508
509
— Computation of Z
- 48 -
-------
— Computation of Z
J
Computation of
real depths.
IF>1514,516,1514
VALE«Z(n.M-1)
GO TO 517
VALEsZ(K!.M>
IF(ME-M)518,521.518
IF(HTZ519.521.519
VARIsZ(N,M*l>
00 TO 522
VARIsZ(N.M>
VPB(N.H">bALPHA*Z(N»M)«BETA«(VAUP*VALO*VALE*VARI >
CONTINUE
DO 535 :"!s2,NE
DO 535 M«3.MEH
IF(HTZ(N,M)) 535,535.2533
lF(Hfz(NjM).3.>533,535,533
Z"CN#M>«VPB(N,M)-A4**U-HGU{N,H»1>*U(-V,M«i )*HQV
1*V(N-1.H>»HGV(N,M>*V(N.M))
535 CONTINUE
DO 1543 N s l.NEH
DO 1543 M = l.MEH,
IF(HTU(N.M> > 1541.1541.1540
HGU n HTU(N.M) ~+Z>/2, |_
IF(HTV(N,M>) 1543,1543,1542
HGV(N.M) s HTV(N.M>*(Z(N,JH)»Z(N*1.M) )/2,
CONTINUE
DO 1568 Nsl.NE
DO 1568 M«1,ME
IF(HTU(N.M)>1565,1565.1563
IF(HGu(N,M).10,>1564,1564,1565
HGU(N,M)1*3.
IF(HTV(N.H)>1568,1568,1566
fF-10,>1567,1567,1568
HGV(N,M)8-3.
CONTINUE
DO 585 Nc2.NEH
DO 585 Mi2.MEH
IF(HTU>585,585,543
IF(1-N>544,547,544
IF(HTU(N-1,M)>545.547,545
VUP=U(N*1,M>
GO TO 54B
VUP=U(N.M)
IF(NE-N>549,552.549
IF(HTU(N«1.M> >550 .552,550
VLO=U(N*l,M)
GO TO 553
VlO=U(N,M>
IF(1«M>554,557,554
IF{.HTU(N,M-i> >555,557,555
VLE = U(N-^f,l>
514
1514
515
516
517
518
519
m
522
523
530
531
532
2533
533
1537
1538
1539
1540
1541
1542
1943
1560
1561
1562
1563
1564
1565
1566
1567
1568
540
541
542
543
544
545
546
547
548
549
550
551
592
553
554
555
Checking of
"drying" of tidal
flats.
Computation of U
- 49 -
-------
556 GO TO 558
557 VLE = U(N»M)
558 IF(ME-M>559.562.559
559 IF(HTU(N.M*1>>560,562.560
560 VR!cU(N;'M + 1)
561 QO TO 563
362 VRI=U(N»M>
563 UPS(N,M)»ALPHA*u(N,M)~BETA*(VUP*VLO«.VLE + VRI)
564 IF(HtV(ti-l.M> >565.569.565 =
565 VUL=V(N-1,M>
566 I'F(H7V(N-1«M+1) >567»570»567
567 VURfV(N»l,M+l)
568 QO TO 571
569 VUL=V(N"1.M+1)
570 VUR=VUL
571 IF(HTV(N,M)>572.576.572
572 VLl=VCN,M>
573 IF(HtV(N.M+l>>574.577.574
574 VLR = V(Ni M*1>
575 GO TO 578
576 VLL=VCN.M*1>
577 VLRsVLL
578 VPB(N.M>8(VUL+VLL+VUR*VLR>/4,
585 CONTINUE
586 DO 595 N«2.NEH =
587 DO 595 M*2,MEH
588 IF(HTU(N.H)) 595.595.1570
1570 IF(HGu(N.M)-3,) 595,1571.589
1571 VPB(N.M>b0.
1572 GO TO 595
589 IF(UPB(N.M)*UPB(N.M> >595.590i59i
590 IF(VPB(N.M)#VPB(N.M> >595,595,591
591 ROOT»SQRTF{UPB*UPB*VPB*VPB >
592. Gr2=A3*600T
593 VPB(N.M)* (1, - G9Z/HGU(N,M>)#UPB ~ A2*V
1CZ(N,M+l> - Z(N,H>> ~ XK(N.H)
595 CONTINUE
1590 CONTINUE
3596 DO 3600 N»1,NE -
IF«HTU < N.1))3600,3600. 3597
3597 IF(HTU(N» 2)>3600.3600, 3598
3598 VPB(N.l>sVPB(N,2>*HTU(N«l>/HTU(NI,2>
3600 CONTINUE _
900 Do 918 N«1,NEH ~
901 DO 918 M 8 2,ME
902 IF(HTV(^,M)>918,918,903
903 IF(HTU(N,M-1)>904.908.904
904 VUL=U(N.M-1>
905 IF(HTU(N*1,M«1)>906,909»906
906 VLL=U
907 GO TO 910
Computation of
"Computation of
U.
PB(N»M> • A5*
Special boundary
settings.
- 50 -
-------
906
909
910
911
912
913
914
913
916
917
918
920
921
922
923
924
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
960
961
962
1990
1991
1992
963
964
965
966
967
VUL=U(N*1«M-l>
VUL=VUL
IF(HTU(N.H))911.915,911
VUR=U(N.H>
IF(HTU
GO TO 917
VUR=U(N*1,M>
vLr=vur
UPBe(VUL*VLL*VUR*VtR>/4,
CONTINUE _
GO TO 920
DO 924 Nal.NE
DO 924 M«1.HE
IF CHTUCN.M)>924.924.923
U(N#M) f VPB(N.M)
CONTINUE
Computation of u
Transfer of new u
into proper array,
— Computation of V
DO 954 N*1.NE
DO 954 Hsl.ME
IF >954,954.933
IF(1-N)934.937,934
IF(HfV{N-l.M)>935.937.935
VUP=V(N-1.H>
GO TO 938
VUP=V(N'H>
IF(NE«N*939.942.939
IF>940»942,940
VLO=V(N*l,M>
GO TO 943
vLosv(n;h>
1F<1-H>944,947,944
Ip>945.947.945
VLE=V(N.M.d
GO TO 948
VLE=V>950.952.950
vri=v(n;m*i>
GO TO 993
VRIsV(N.M)
VPB(N.Mi?»ALPHA*v(N.H>*BETA*(VUP*VlO*VI-E+VRI >
CONTINUE |
DO 969 N ¦ 2.NEH
DO 909 M»2.HEH
IF(HTV(N#H> > 969,969,1990
IF(HGv 969,1991,963
V(N # H)a 0,
GO TO 969
IF(VPB(N.H)*VPB(N.H>>969i964#965
IF(UPB(N.H>*UPB(N.H> >969.969,965
R00TcSQRTF(VPB*VPB*UP8)
1) - A2*UPB(N»H> • A5*(Z(N.M>
- 51 -
-------
~ 113 ——-
FLOW THROUGH
- special inflow setting.
the sections
1 - Z(N+1»M)) «• YK(N.M)
969 CONTINUE
V(26,40> ¦ V(26.40 )
: computations of net
DL.2 S 27 * DL
IF(T-TRB) 400,5250.5250
5250 IF(T-TRE) 5251,5251.400 ,
5251 AAlsAAl+(V(25,40)*HTV(25,40)*Al#Dl.2)/1000006l
BB2=BB2+(U(23#39)*HTU(23»39>*Al*DL2+U(24t39)*HTU(24#39)*Al*DL2
' ~U<25,39)*HTU(25.39)*A1#DL2)/100 0 00 0I
1
CC3 = CC3*(V(22*39)*HTV(22i39)*Al#DL2*V(23«39) *HTv<23i39>*Ai#DL2*
1V(24,39-)#HTV(24,39)*A1*DL2)/10 0 000 0 .
TRANSPORT THROUGH SECTIONS
TRANS(UORV.hTUORHTV,N,Mi KEY.NUMBER)
s
s
1960
i
400
401
402
403
TRANS(U#HTU'.36',7,2,5)
TRX*A1*DL2/1000000.0
TRANS(U,HTU.37,15.2#3)
TRX*Al*DL2/lflO0b0Ot0
TRANS(U»HTU'.32',27f 2.5)
TRX*A1*DL2/1000000 .C!
TRANS(U,HTU.23.39.2.3)
TRX*A1*DL2/1000000.0
TRANS(U/HTU,23,37*2.3)
TRX*A1*DL2/1000000.0
T^ANS(U.HtU.20.35>2«6)
TRX*A1*DL2/1000000#6
T ;ANS(U#HfU.17.27#2.10)
TRX*A1«DL2/1000060.0
TRANS(U#HfU»12.18*2.18)
TRX*A1*DL2/1000 000,6
TRANS(U,Hfu.6.9,2,28)
TRX*AI*DL2/I00p000,p
TRANS(V.HtV.23.37*1.3)
TRX^A1*DL2/1000000,0
TRANS(V.HtV.20,35.1.5)
TRX*A1*DL2/1000000,0
TRANS(V.HTV,17,27»1.13)
TRX«A1*DL2/1000000.0
TBANS(V»HtV,.l2.18»l»23)
TRX*A1*DL2/1000000.0
TfiANS(V.HTV.6.9.1,34)
TRX*A1*DL2/1000000.0
CONTINUE
COMPUTATION OF REST CURRENTS
IF410.410,401
IF(T~TRB)410.402.402
IF(T-TRE)403.403,410
TIC«f ip'+Al
DO 4090 N»1# NE
DO 4090 H°1#ME
IF(HTZ(N.M))4090.4 090.4O9i
TRX
AA2
TRX
AA3
TRX
AA4
f RX
BBS
TRX
BB4
TRX
BBS
TRX
BB6
TRX
BB7
TRX
BBS
TRX
CC4
TRX
CCS
TRX
CC6
TRX
CC7
TRX
CC8
x
s
X
s
3
S
8
3
S
s
r
s
Computations of currents
through special sections.
Accumulation for rest
current computation.
- 52 -
-------
4 091 USR (N, Ml sUSR (N,M)*(U2401,2401,2403
2401 IF(S(N,N.i))2402»2402.2403
2402 IF2404,2404.2403
2403 US(N#H)5US(NiM)*ui[N»M)i>2«DT
GO T0_ 2 405
2404 US(N,M>*0, .
2405 !F(S(N.H)}2406,2406.2408
240 6 IF(S*Q.
409 CONTINUB
B4b2•*Dt
DO 2610 N»1,NE
DO 2610 M«l,ME
IF(ABsFCUS(N,M)J/B4.1,) 2604,2601.2601
2601 IF(US(N»M)) 2602.2603.2603
2602 US(N,M>*»B4
GO TO 2604
2603 US(N,M)sB4
2604 IF(ABSFJ/84-1.> 2610•26q5,2605
2605 IF(VS(N.M)) 26o6.2607.26o7
2606 VS(N,M)?*B4
GO TO 2610
2607 VS(N,M)*B4
2610 CONTINUE-
410 IF(LI>984,984,411
411 CONTINUE
IFCT-IOBOO .) 412,4500 ,412
4500 CONTINUE
J
Setting of current
"accumulations" for
computations of
advection of pollutant
S(22,35)
S(22,27>
S(5,6)
S(4,28>
S(7>18)
S(6,39)
S(12,9)
S(12,32)
S(17,18>
S(17,37>
S(22,9)
S(27,18)
S(32,7)
S(34,27>
S(38,6)
S(38»l5>
10000.
10000.
10000,
10000.
10000.
10000.
10000.
10000.
10000,
10000,
10000.
10000,
10000.
10000.
10000.
10000,
Instantaneous
pollutants.
release of
53 -
-------
412 BlaDL'Dt.
AUS(1)*90000 .
413 B2 s (4;o#DT*Aus(l))/01
415 B5 = (AUS(l) * DT)/B1
416 Y1«0.
S<26»40> • S(26.40> ~ 3t continuous release of pollutants,
; COMPUTATION OF DISPERSION AND DIFFUSION
417 DO 3000 N*2»NEH
DO 3000 Ms2,MEH
IF(HTZ(N.M> >3000,3000.201
201 IF(HTU(N4H)>202.202.204
202 SSW=0,
GO TO 225
204 IF(U(N.M)>205.205.211
205 IF(US(N.M)~04)2209,2206. 2209
2206 IF(S(N.M*1>>2207.2207.2208
2207 SSH=S(N.M>/(DT*ABS(ll(N.M) >
US(N.,M>±0.
GO TO 225
2208 SSH=S(N'M*1>/(DT*ABS(U(N.M))
US(N.M > *6 «
GO TO 225
2209 SSH=-S(N.M*1>>/DL2
GO TO 225
211 IF(US(NiM)-B4)2213,2210.2210
2210 IF(S(N.M.l)>2211.2211.2212
2211 SSH=S(N.M>/(DT*A8S(U(N,M>>
US(N.M)*0 .
GO TO 225
2212 SSH=S(N.M*1>/(DT*ABS(IJ(N.M) >
US(N,M)*0,
GO TO 225
2213 SSH=(S(M,H)-S(N.M-1>>/DL2
225 IF(HTV(N.M>>2214,2214,2215
2214 IF(HtUCN.M)>3000.3000.2230
2230 SSV=0.
GO TO 4428
2215 IF(V(N.H))2216,2216,2221
2216 IF(VS(N.M)*84)2220,2217.2220
2218 SSV=S(N.M>/(DT#aBS(V(N,M> >
VS/(DT*ABS(V(N.H> >
VS(N» M > *0 .
GO TO 4428
2220 SSV=(S(N.H>-S(N-1.M>>/DL2
GO TO 4428
2221 IF(VS(N.M)-B4>2225,2222,2225
2222 IF(S(N + 1.M> >2223.2223,2224
2223 SSV«S(Nf'M>/(DT*ABS)
VS(N.M>±0.
GO TO 4428
-— Computation of advection
parameters.
- 54 -
-------
Computation of
and advection.
diffusion
2224 SSV=S(N*1.M>/(DT*ABS(V(N*M>)
VSCN,M)=0,
GO TO 4426
2225 SSVs(S(N,M)-S(N*l.M)I/DL2 ~~
4428 S(N»H)sSrN,M)-B2*S(N.M)*BS*{S(N»ljM>*S(N*l«M)*S(N«Mal)*.S(N,M+l)
1*4,*S(N»M)) |
2» >*SSH)-(DT*A0SF(V 2429,2429-300
2429 S(N,M> * 0,
3000 CONTINUE
300 CONTINUE __
C 415 B5«2,0*B4
C 416 DO 422 N«2,NEH
C 417 DO 422 H«2,MEH
C 410 IF(HGu«S(N.M)-Al*U{N.M)MSOJiMM)/B4*StN»M.l)/B4>*B2*+S(N+1,H)/2,)
C 421 S(N*M)sft(N*-M>-Al*V(N,,M)*(S(N+l'M>/B4*S(N»l*M>/B4)*B2MS'<'N»l«M')/2,*
C 2StN*'M)*S(N»ifM)/2»)*B2*tS(N«H«l)/2i"S(N»H)*S(NjM*l)^2t)
C 422 CONTINUE
904 IF(Cl)4'0i472» 470
470 IF(T-TWV472.472.471
471 CALL J0B4
472 CONTINUE
RETURN
END
FUNCTION TRANS(FV,FH,NtM,KEY»NUM>
DIMENSION FV(42.47)»FH(42i47)
C FV « U OR V COMPONENT OF VELOCITY
C FH s HEIGHT ON U OR V GRID
C N « INITIAL VALUE OF N.,,,MIN
C M « INITIAL VALUE OF M....MIN
C KEY b 1,,,section IS along constant N
C KEY = 2...SECTION IS ALONG CONSTANT M
C NUM = NUMBER OF POINTS ON THE SECTION
C
TOT « 0
Ml = M-l
N1 = N-J
GO T0(100 »200 ) KEY
100 CONTINUE
THIS IS A SECTION ALONG CONSTANT N '.,.V VELOCITY USED
00 150 k B 1,NUM'
TOT « TOT ~ FV(N,Ml*K)*FH(N.Ml + j<)
150 CONTINUE
tRANS = tOT
RETURN
— Function to compute
transport through section.
- 55 -
-------
200 CONTINUE
SECTION ALONG CONSTANT
1,NUM
~ FV(N1*K»M)*FH(N1+K»M>
• • I *
U VELOCltY USED
: THIS IS A
DO 250K a
TOT = TOT ~
25o continue
TRANS s TOT
RETURN
END
SUBROUTINE J0B3
DIMENSION HTZ<42,47>,HTU< 42.47).HTV<42,47>.Z<42,47),U<42,47),V(«2,
147>,US(42.47),VS(42.47),RAD(42#47),ANG(42#47),XK(42i47),yK(42,47).
2 NU'<6).MUC6).NK<85).MK»AH<5.85).CAPAP(5,85)iPReC(5,85>.ARG
3(5,85),PERU(85),pERV(fi51»NL(85).ML(85),NZ(10)iMZ(l0)iVl,V4<10),HGUt42,4
57),HQV(42,47),H(5>,8{5)»AL(5),CapA(5)#CAPB(5)*A1n3(5)#CAjnC(5),
6 NG(40),HG(40)«TlT(3)»AUS(5)tTAR(85),0TC(85)
7,USR(42,47),VSR(42,47>
COMMON yTZ.HTU,HTV.Z,U,V,US«VS»RAD»ANG»XK»YK» NU#mUiNK»MK.AH»CAPA
1P.PREC,ARG,PERU.PERV.NL•Ml,NZ.MZ#VI.UPfl.VPB.S.A ,V2,V3,V4,HGU,HG
2V,H.B.AL.CAPA»-CAPB.AlNc«CAlNCiNG* MG.TIT, aUS.TaR«0Tc#USR#VSR,
3Tl.T2.T2.TE,TW,Al.A2,A3.A4,DT.BETA.ALPHA.Ff R,DL'.G.C»C3f NE.ME.NEH.M
4EH, JA, IUE,KKE.NURlJ,Sr,T,RBETA,Ci,SlGMA,ROL, A5,RU#*0,II,LE,LU,KI
5,L0P.TRB,TRE,TIC, LUM.LA.TIL.LIV.LEET
6,P0SN.P0SM.DRI,STV,ETV,TINCV,STH,ETC,TINCH,CNCU,CAIMCU
7,AAl,BB2,CC3,TSC,DL?
8,AA2,AA3,AA4,BB3,F'B4,BB5,BB6,9B7,BB8#CC4fCC5,CC6lCc7,CC8
I WRITING OF OUTPUT FIELDS
31 F0RMAT(9F8.3>
32 FORMAT(9F8.0)
601 CONTINUE
610 PRINT602.T, (N.NbI',16)
611 PRlNT604,(N,(ZCN,M),Mal#16)#Nil.NE)
612 PRINt603.
613 PRINT604, ,M»33.47)»Ncl,NE)
1605 FORHAT(//8X.15 J 8 )
1604 FORMATUX, 18.15F8.0./)
614 DO 639 N*1.NE
617 DO 639 MbI.ME"
618 IF(V(N/n))620,619,620
619 IFCU(N.M))620 #622.620
620 RAD(N.M)bSQRTF(A8SF)**2*aBSF(U(N.M))**2)
621 GO TO 625
622 RAD(N,M)>0.
623 ANG(N,M)*999.
624 GO TO 639
625 ANG(N,M>«ABSF(ASINF))
626 ANG(N>M*«ANG*(36n./6,283lB53073)
627 IF(V(N.M)>629,628.628
628 IF(U(N|M)>631,636,636
629 IF(U(N,M))633.633.630
630 GO TO 635
- 56 -
-------
631
632
633
634
635
636
637
638
2638
639
650
651
652
653
,M>
1608
1609
660
661
662
670
1660
1661
1662
1664
1665
1666
1667
1650
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
640
641
PRJfiT
PRlNf
PRINT
print;
ANG(N*hJ«180,-ANG> ~ TIL
IF(ANQ(N« M))638»639,639
ANG
GO TO 637
CONTINUE
607,T»(N,Nb1,14>
696*(N,(RAD(N»M>.ANG(N
609,(N*Nol5»28)
6^8.(N,(RAD(NiM).ANG(N
609*(N,Nb29.42)
6Q8*(N*(RAD(N,M).ANG(N,
PRINT1609,(N,Na43,47)
PRINT16O0.(N,(RAD(N,M)»ANG(N»
FiDRHAf{(lX, !8,5(F4t0,lH« »F4, 0
Format(//8x. 5ie>
IF(LI) 1660.1660,661
PRINT 6*0,T.(N,Ntl,16)
PR|NT 6041 (N«(S(N»M)'*M8lii6),
PRINT 605*(N*N«17,32)
PRINT 604,(N,CS(N,M>,Mb17,32)
PRlNfl605.(N*Nf33.47)
PPINT1604.(N.(S(N,M).Ms33.47)
FpRHAT(lHl.i4X,22HCONCENTRATI
IF{LEEf> 1650.1650,166j
IF(T»T4) 1650*1662.1650
WRlTE(5;3l)( (Z(N.H).^{.ME),N
WRITE (Pi3l)((U(N.M).Mol.ME)»
WRITE (5,31).((V(N*M)«MalfME)«
WRITE (*,32)T
ENDFItE 5
IF(LOP-IO) 1630.640*1630
IF(LOP*5)640~1631,640
DO 1636 Nfl.NE
DO 1636 M«l,ME
tF1634,1635,1634
V(N,M)aVSR(N«M)/.(0,01*;!C)
IF(USR(N,M))1636,1637,1636
U,)
M=1,NE)
Printing of output
fields.
Nbi;
NE)
*Nb4Jne)
ON AFTER Tb¦* F9 # 0 •¦/¦/8X * 18•' 1518//)
¦iane;
NPl.NBr
N»l> N|)
Writing of outputs on
tape.
— Rest current output.
- 57 -
-------
Writing of data
points on tape.
from selected
642 M=MZ(I)
643 V2(I)=Z(N«M)
644 V3(I)=RAD(N.M)
643 V4(I)=ArjG(N.M)
646 CONTINUE
647 WRITE(2>T«(V2(I)iI«1#NURU)#(V3(j).I=1«NURU).(V4(I),I»liNURU)
602 F0RMAT
AND
608 FORMAT((I8,14(F4.0.1H.F4
609 FORMAT(fiX«17*1319/)
1670 FORMAT(///.5X.31HP0SIT ION
1/)
1651 IF(L!V)654•654,1652
1652 IF(LIV-2)1653,654.1653
1653 PRINT 1670.POSN.POSM
654 CONTINUE
3610 IF(T-TRB) 3620,3611,3611
3611 IF(T-TRE) 3612,3612,3620
3612 AAl=AAl/3600,
AA2 = A A 2/3 6 0 0 ,
AA3 » AA3/3600,
AA4 > AT4/3600.
BB2=BB2/3600,
BBS a BB3/3600 ,
BB4 s BB4/3600,
BBS e Ba5/3600,
BB6 b BB6/3600.
BB7 s BB7/3600,
BB8 s BB8/3600.
CC3=CC3/3600,
CC4 s CC4/3600,
Cf5/3600,
CC6/360 0 ,
CC7/3600 ,
CC8/3600,
5i 01» AA1,BB2,CC3
5101.AA2.AA3.AA4
5101.BB3.BB4.BB5
5101»BB6,BB7.8B8
5101.CC4,CC5,CC6
5101.CC7,CCS
0)/))
_0F DRIFTING OBJECT, Nb, f6 ,2,2Hm»,F6.2»//
Output of positions of drifting
objects.
5101
CC5 s
CC6 b
CC7 s
CCS =
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
FORMAT (J-H1.//10X.26HTRANSPORT. SECT I ON
15x.F1Q >0.//28x,2HC,,5x.Fl0,0.//>
Output of
sections.
transport through
J
A, M3" ,F1qi0«//28X,2HB,
AA1 s
AA2 s
AA3 =
AA4 b
BB2 s
0.
0.
0T
o;
07
- 58 -
-------
BBS
BB4
BBS
BB6
BB7
BBS
CC3
CC4
CC5
CC7
cca
o:
o:
o.
o;
o;
o;
o;
oT
o;
o:
o:
o:
Zeroing of transport array,
3620 CONTINUE
RETURN
END
SUBROUTINE jOB4
DIMENSION H7Z(42*47> •HTU(42«47)«HTV(42*47)f Z(42*47)fU(42«47> »V('42«
147) ,US(i'2.47),vS(42.47)#RAD(42,47),ANG(42,47)lXK(42#47),YK(42,47)<
2 NU(6)»HU(6)»NK(85)»MK(85)»AH(5»85>#CAPAP-<5«85)iPREC(5,85).ARG
3(5l85),PERU(fl5).pFRV(85)#NL(83).ML(85),NZ<10).MZ(10>«Vl(10)#UPB(42
4»47).VPB(42.-47),S(42.4 7).A(65) •V2<10>.V3(10)•¥4(10),HGU<42.4
57),HGv(42.47),H(5),B(S),AU(5)»CAPA(5)«CAPB(5).AlN"(3)iCAiNC(5),
6 NG«4O)(MQ(4O>it!t(3)-AUS(5).TAR(85)#OTC(05)
71USR(42 #47),VSR(42,47)
COMMON HTZ,HTU,HTV,Z,U,V,US»VS«RADiANG.XK.YK. NU#MJ#MK#MK,AH,CAPA
1P»PREC#ARG»PERU»PERV»N|.iML*NZ«MZ»V1#UPB#VPB»SiA ,V2,V3,V4,HGU#HG
2V,H,B.AL.CAPA,CAPB,AiNC.CAlNC.NG,HQ,TIT,AUS,TAR,0TC#USR,VSR,
3T1,T2,T:,TE,TWja1.A2,A3»A4,DT.BETA#ALPHA,F,R,DU#G»C»C3#NEiME,NgH.M
4EH,JA,IUE,KKE,NURU,SI,t,RBETA* Cl,SlGMAiROL,A5,RU,RV# KO, LI i LE, LU, KI
5, LOP, TRp, TRE, TIC, LUN.i^A.TjL.LlV.LEET
6,posn#posm»dri,stv,etv,tincv»sth.etc#tinch,cncu«cajncu
7,AA1,BBB* CC3« TSC« DL2
8,AA2,AA3,AA4,BB3,BB4fBB5,BB6,BB7,BB8,CC4,CC5,CC6#CC7#CC8
C WIND CURRENT COMPUTATION AND COMPUTATION OF
C DRIFT OF BOATS AND OIL
704 CONTINUE
705 if(Sl«T)745»745#706
706 KKE=(lU&/2>*3
707 lUoIUE/2.
708 K=XFIXF ((T-TW*A(1))/A(1))*1
709 L a K +1
715 Do 743 Ul, IU
716 NIbNU(2*I»l)
717 N2bNU(2*M
718 Ml»MU(2*J«l)
719 M2pMU(2*I)
720 DO 743 N « N1.N2
721 DO 743 H«M1,M2
722 IF(HTU(N.M>>733.733.723
723 IF(MGU(N.M)) 724,733,724
724 IP(HQU(N.N)*3.) 733.733.725
725
726
HMLD = 6500 . + 8«* Z(N?M)
IF(HMLD-MGU(N,M)) 730.727.727
Computation of wind current
component and a drift
component due to wind.
- 59 -
-------
727 HMLD = ?Gu(N,M)
730 XKiA(K>*COSF(A(L>>*100,
733 IF(Htv(N,M))743,743.734
734 IF(HGVCNfM)) 735,743.735
735 IF(HGv(N»M)4-3, ) 743.743.736
736 HHLD»2500.~8,*Z(N.M)
IF(HMLD-HGv(N»M) ) 740.737, 737
737 HMLD = &GV(NiM)
740 YK*C3#A*0,0174533)/HMLD
741 IF(LIV>743.743.742
742 VPB(N.M>bA(K)*SINF(A(L)>*100,
743 CONTINUE
744 Go TO 7^0
745 DO 749 N»1,NE
746 DO 749 MbI.ME
747 XK(N* M)e0,
748 YK(N.H)s0.
749 CONTINUE
750 CONTINUE
751 IF < LIV > 780,780.752
752 IF(T-TRB)780,753,753
753 IF(T-TRE>754,754,780
754 IFCLIV-^ > 770.755, 755
755 PRODs0.015
756 DO 766 N«2.NEH
757 DO 766 Ms2,MEH
758 IF(HTZCN.M)) 766.766.759
759 UMR=UPB(N,M)*PR0D*A4
760 UML=UPB(N,M-1>*PR0D*A4
761 VHBsVPB(N,M>*PR0D*A4
VM0?VPB(N"1.M>*PR0D*A4.
Computation of drift of
— objects and oil at the
surface (due to wind)
) )«-UML*(S(N#M»i)»S ).VHB*(S(N
S(N.M>sS(N»M>+UMR*.S)-VMO*(S(N-l«M)»S(N«M) >
764 IF(S(N,f)>765,766,766
765 S(N,M)s0 •
766 CONTINUE
770 lF(LlV-2> 771.780.780
771 ND*XFIXF(POSN)
772 MD*XF I Xf. (POSM >
773 UDIS = (A1 ~ (U(Nn,MD> ~ UPB(ND.MD) * DRI> > / 2, * DL
774 POSM=POSM*UDIS
775 VDIS = (A1 * (V(ND,MD>
776 POSNsPOSN+VDIS I
780 CONTINUE |
RETURN
END
SUBROUTINE J0B5
DIMENSION HtZ(42.47>.HtU(42.47>.HTV(42»47>.Z(42.47),U(42,47>.V(42.
147),US(<2.47),VS(42,47),RAD(42.47).ANG(42.47),XK(42#4 7),YK(42,47),
2 NU(6>,MU(6>.NK(85).MK(85>#AH(5.85).CAPAP{5i85>,PREC(5.85)fARG
3(5,85).PERU(85),PERV(85>.NL(85>,ML(85),NZ(10>.HZ(10 >#Vl(10 >,UPB(42
VPB(N)D.MD) * OR I >) / 2, * DL
60 -
-------
4#47)#VPB(42*47),S(42»47)»A(85) * V2( 10 ), V3 (10 >, V4 {10 ), HQU( 42, 4
57),HGv<42,47),Hc5>,B(5),AL(5),CAPA(5)#CAPB(5),AlNC'<5),CAINC(5).
6 NG(40)«MG(40)«tlT'(3)« AUS(5 > * TAR(85)»OTC185)
7,USR(42,47)fVSR{42,47)
COMMON HTZ,HTU,HTV,ZtU,V,US,VS,RAD#ANG,XK.YK, NUiMU,NK,MK,ah*CAPA
lP#PREC/-RG,PERU.PERV'.NL-iMLfN2.HZ,Vl,UPB,VPB',S.A , V2, V3, V4. HGU» HG
2V«H#B#AL»CaPA»CAP8#A!NC»CAINC.Ng»MQ*TIT»AUS,TARfOTc#USRiVSR,
3Tl#T2#TC#TE»TW,Ai,A2,A3*A4,DT»BETA*ALPHA»F,R,DL#GiC»C3,NE*ME,NEHtM
4EH,JA#I0E, KKE* NURU, SI, t, RBETA, Cj[,SlGMA#R0L»A5»RU,RV»K0#Ll»LE#LU«K!
5,L0P,TRB,TRE,TIC, LUN#LA*TlL,LIV.LEET
6,posn*posm.dri,stv,etv,tincv,sth,etc#tinch,cncu.caimcu
7,AAl.BB2,CC3,TSC,nL2
8j AA2,AA3,AA4,BB3, 8B4,BB5,RB6,9B7,BB8,CC4,CC5,CC6,Cc7,CCS
OUTPUT OF SPECIAL POINTS
802 CONTINUE
803 FORMAT(lHl»5X,ll4HWATER ELEVATION/)
804 f ORMAT(®X» 5W{SEC),4X,10(2X* 13,lH,, 13)/15X#5(2X, I3ilHf, J3).//).
805 FORMAT(F11.6#4x.10F9.2/15X,5F9,2)
806 F0RMAT(5X,14HEND OF PROGRAM) |
807 FORMATCJ5X,10/l5x*5(F4.0*lH,,F4,0)')
808 PRINT 803,(Vl(I),Isl.NURU) L Outnut at anecial noints
809 PRINT eo-4«*HZ', l«1,NvjRU>
813 PRINT 805,T. (V2(Ui IbI.NURU) [
814 PRINT 80 7,(V3(I),v4(I),Iel,NURU) |
815 IF (TE"T>816,816,810
816 PRINt 806
RETURN
END.
SUBROUTINE J066
DIMENSION HTZ(42.47).HtU(42,47>.HTV{42#47)#Z(42^47)#UU2.47>.V(42.
vincnajunt # in f # in i *
1^7),US(42,47)lvS(42,47),RAD{42»47),ANG(42#47)tX|<(<2^7)#yK{42,47),
2 NU(6)«MU(6>.NK(85).MK(85>»AH(5,85)«CAPAP(5»85>»PREC(5#85J.ARG
3(5»85),RERU(85)fpERV(85)»NL(85).ML(85)»NZ(lp)«MZ(l0)iVlCl0)iUPB(42
4,47>.VPB<42*47) .S(42* 47')f A (85) .V2(10)«V3(10>•V4<10)*HGU(42.4
57),HGv (42«47 ) # H'(5),B(5), AL (5), CAPA (5)«CAPB( 5) i A InC{ 9) iCA INC (5),
6 N6(40)«M6{40)«fIt(3)»AUS(5)*TAR(ft9)tQTC-(S5)
7 rUSR C 42# 47),VSR(42,47)
COMMON HTZ»HTU,HTV,Z*U.V,US«VS*RAD«ANG*XKiYK, NU«M^'NK^K»AH«CAPA
1P*PREC*£RG.PERU'.PERV.NL«ML*NZ#MZ,V1*UPB,VPB,$»A ,y2iV3,V4iMGU»HO
2V»H»B#AL»CAPA«CAPP»AINC»CAINC»NG»M0«TIT» AUS«TAR^.OTciUSR« VSR«
3U.T2,T4,TE,TW,A1.A2*A3*A4#DTiBETA»ALPHAiF,R,Dl.GiC«CJ#NE»ME#NEHfM
4EH# JArlUE* KKE* NURU, Si,t« RBETA« C^»SIGMA#ROL# A5*RU« &V»K0i L'l«LE#LU» Kl
5(L0P,TRB,TREiTIC. LUN»LA#TiL,LIV,LEET ,
6«P0SN*P0SM»DRI,STV«ETV,TINCV*STH#ETC# TlNCH;CNCU«CA JNCU
7,AA1,BB2,CC3.TSC,DL2 _
8,AA2,AA3,AA4,BB3,BB4,BB5#BB6,BB7,BB8,CC4,CC5,CC6,CC7#CC8-
ENTRY TO READ VECTOR ~ HEIGHT CARDS
ETH a ETC __
200 READ(50»215)STV*Etv»TlNCV
201 PRINT 216*STV*ETV*T1NCV
202 READ(50,215)STH.ETC,TiNCH
61 -
-------
Outputs for plotting.
203 PRINT 216,STH,ETC.TINCH
215 FORMAT(3F10•0 )
216 FORMAT(10X«3F10,0)
C PAUSE 3 IS FOR MOUNTING OF PLOTTING TAPE
C 217 PAUSE 3
218 WRITE (3)HTZ*V1,NURU,ME,NE
RETURN
END
SUBROUTINE JOB7
DlMENSlbN HTZ(42.47).HtU(42,47).HTV(42»47),Z<42,47)#U(42,47).V(42.
147),US(42".47)#VS(-42.47),RAD(42»47)#ANG(42,47)iXK<42i47),yK<42i47).
2 NUi6).MU(6)»NK(85).MK(85)»AH(5,85)#CAPAP(5.85)»PREC(5,85)iARG
StS^Sj^ERUfeSj^PERVtSS)! NL(83>,ML<85},NZ(ld).MZ{10)# V1(10),UP8(42
4#47)#VPB(42*47),S(42.47),A(85) »V2(10).V3,HQv<42M7>,H<5),B(5)#AL(5>.CaPA(5>»CAPB<5)iAInC(3)|CAINC(5>,
6 NG(40),MG(40)»TlT(3)»AUS(5)iTAR(85)|0TC<85)
7,USR(42»47),VSR(42,47)
COMMON HTZ»WTU»HTV«Z#U»V#US«VS#RAD#ANG*XK»YKj NU»Mll#NKiMK, AH. cap a
1P.PREC»ARG.PERU.PERV'.NL"«ML»NZ.MZ. Vl.UPB# VPB.S. A . V2,V3,V4.HGU.HG
2V.H#B»AL.CAPA,CAPB.AiNC.CAINC.NG#MQ.TIT,AUS»TAR'iOTciUSR#vSR,
3tl,T2.t*#TE»TW,Al.A2,A3.A4,DT.BETA#ALPHA»F,R,DL»GiC»C3,NE.ME,NEH.M
4EH,JA,JUE,KKE,NURU,SIJt»RBETA,Ci,SlGMA,R0L,A5,RU,3V#K0,Ll,LE,LU,K!
5,L0P,TRB.TRE.TIC, LUM.LA.TIL#LIV,LEET
6,P0SN.P0SM.DRI,STV,ETV#TlNCV,STH,ETC#TlNCH,CNCUiCAlNCU
7,AA1,BB2,CC3,TSC,DL2
8,AA2#AA3,AA4,BB3,BB4,RB5fBB6,BB7,BBS,CC4,CC5.CC6,CC7,CCB
ENtRY TO WRITE DATA ON TaP-E FOR VECTOR AND WATER HEIGHT PLOT
230 VEC=HT»0,
231 IF«Tl/2,>232.232.234
232 STVsSTV*TlNCV
233 VECS1.
234 IF(ABsF(T«STH).Tl/2.5 235,235,237
235 STH=STH+TINCH
236 HT«1.
237 IF(HT+VEC)238,239.238
238 WRITE (3)T,VEC,WT,RAD,ANG,Z
239 continue
REtURN
END
function asinf(*>
ASINF = ASIN(X)
RETURN
END
FUNCTION SINF(X)
SINF = SIN(X)
RETURN
END
FUNCTION COSF(X)
COSF = COS(X)
RETURN
END
FUNCTION SORTF(X)
SQRTF * SQRT(X)
RETURN
END
- 62 -
Outputs for plotting.
-------
FUNCTMJN ABSF(X>
ABSF = ABS(X)
REtURN
END
FUNCTION XFIXF(X)
XFlXF a 1F1XIX)
RETURN
END
------- |