Approved for Public Release;
Distribution Unlimited.
ENVPREDRSCHFAC
Technical Note No. 2-
A MULTILAYER
HYDRODYNAMICAL-NUMERICAL MODEL
(W. Hansen type)
MODEL DESCRIPTION AND OPERATING/RUNNING INSTRUCTIONS
By
TAIVO LAEVASTU
ENVIRONMENTAL PREDICTION RESEARCH FACILITY
Part 2 of a series of four reports
prepared for:
ENVIRONMENTAL PROTECTION AGENCY
PACIFIC NORTHWEST ENVIRONMENTAL
RESEARCH LABORATORY
CORVALLIS, OREGON
JANUARY 1974
ENVIRONMENTAL PREDICTION RESEARCH FACILITY
NAVAL POSTGRADUATE SCHOOL
MONTEREY, CALIFORNIA 93940

-------
UNCI asst ft Fn
SECURITY CLASSIFICATION OF THIS PAGE r»Wian Data Enfrtd)
REPORT DOCUMENTATION PAGE
READ INSTRUCTIONS
BEFORE COMPLETING FORM
\. REPORT NUMBER E N V P R E D RS C H F A C 2. GOVT ACCESSION NO.
Technical Note No. 2-74
3. RECIPIENT'S CATALOG NUMBER
4. TITLE (end Subtitle)
A Multilayer Hydrodynamical-Numerical
Model (W. Hansen): Model Description and
Operating/Running Instructions
3. TYPE OF REPORT ft PERIOO COVERED
•• PERFORMING ORG. REPORT NUMBER
7. author^;
Taivo Laevastu
«. CONTRACT OR GRANT NUMBERO)
9. PERFORMING ORGANIZATION NAME AND AOORESS
Environmental Prediction Research Facility
Naval Postgraduate School
MontereyCalifornia 93940
10. PROGRAM ELEMENT. PROJECT TASK.
AREA ft WORK UNIT NUMBERS
11. CONTROLLING OFFICE NAME ANO ADDRESS
Naval Air Systems Command
Department of the Navy
Washinaton. D.C 20361
12. REPORT DATE
January 1974
IS. NUMBER OF PAGES
60
M. MONITORING AGENCY NAME ft ADDRESSf/f dUftrtnt from Controlling Office)
15. SECURITY CLASS, (of thl» report)
UNCLASSIFIED
11a. DECLASSIFICATION/DOWNGRADING
SCHEDULE
16. DISTRIBUTION ST ATEMENT (of thi* Report)
Approved for Public Release; Distribution Unlimited.
17. DISTRIBUTION STATEMENT (ot th» mbetrmct entered In Block 20. II dJ//aran* from Ropori)
IS. SUPPLEMENTARY NOTES
19. KEY WORDS (Conttnuo on rsvarss elde If neceeemjj and Identity bjr block number.)
Multilayer Hydrodynamical-Numerical Model
Dispersion and diffusion of pollutants
20. ABSTRACT (Contlnuu on revere* elde II neceeemry and Identify by block ntmbmt)
The multilayer Hydrodynami cal-N-wner-i cal (HN) model (W. Hansen
type) is, in many respects, similar to Hansen's single-layer HN
model. The hydrodynamical formulas present vertically integrated
motion through given layers. The lower layer is partly driven by
friction, which is proportional to the mean density relations of
the layers,and partly by horizontal pressure gradients caused by
the inclination of the surface._as well as of the Interface.
DD ,	1473 EDITION OF 1 NOV 65 IS OBSOLETE	UNCLASSIFIED
S/N 0102-014-6601 |	,•		Um,LHJJiriLU	.
' SECURITY CLASSIFICATION OF THIS PAGE (Whan Dmt* Knf?»d.

-------
	UNCLASSIFIED	
¦.lUJHITY CLASSIFICATION OF THIS PAGEftVhan Dale Entered)
20. (Continued)
The finite difference forms are similar to those used in single-
layer models, except for the additional terms dictated by the
second layer. The model uses a continuity equation in its direct,
physically simple and correct form.
In respect to the inputs, the model is also similar to the
single-layer HN model, with the exception of additional inputs
of second layer thickness and internal tides in the second layer.
Besides the basic formulas of the model, this report describes
the numerical model in FORTRAN. The description of the input
parameters, abbreviations used, and the program itself are given
in the Appendixes.
The verification of the model is described in another report.
The model described in FORTRAN has not been fully optimized
forroutine production runs.
UNCLASSIFIED	
SECURITY CLASSIFICATION OF THIS PAGEfHTiwi Date Entered)

-------
FOREWORD
This documentation of a two-layer hydrodynamical-
Numerical (HN) model is the second in a series of reports
from Environmental Prediction Research Facility (EPRF), to
Environmental Protection Agency (EPA).
This two-layer HN model is in many respects similar to
the single-layer HN model of type Walter Hansen and reference
is often made to the latter in this report. The model was
applied to the New York Bight and the results and verification
are given in another report. It could be pointed out that
similar models with three and six layers have been applied at
EPRF for other areas of the Northern Hemisphere oceans.
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 Pollution
Research Program of the above mentioned EPA laboratory.
The project was carried out by the Oceanography Department
of the Environmental Prediction Research Facility under the
supervision of Dr. Taivo Laevastu.
iii

-------
CONTENTS
FOREWORD	111
1.	THE BASIC MULTILAYER HYDRODYNAMICAL EQUATIONS AND
THEIR FINITE DIFFERENCE FORMS 	 		1
2.	THE ESSENTIALS OF THE MULTILAYER HN MODEL		6
2.1	The Stability Criteria and Horizontal
Viscosity Coefficient 		6
2.2	The Grid Net and the Inputs 		7
2.3	Treatment of Open Boundaries		7
2.4	Other General Information, Including Notes
on Verification				8
3.	REVIEW OF THE FORTRAN PROGRAM	10
3.1	The Control Program	10
3.2	Subroutine J02	10
3.3	Subroutine J03	11
3.4	Subroutine Vector 		12
3.5	Subroutine J04			12
3.6	Subroutine J05	13
3.7	The Transport Function	16
3.8	Subroutine SOI	16
3.9	Subroutine S08	16
3.10	Subroutine J0B6, J0B7 and the Function
Subroutines	16
REFERENCES	17
APPENDIX A - INPUT PARAMETERS AND FORMATS 		20
APPENDIX B - LIST OF SYMBOLS		23
APPENDIX C - PROGRAM NYAREA LISTING 		29
v

-------
I. BASIC MULTILAYER HYDRODYNAMICAL EQUATIONS AND
THEIR FINITE DIFFERENCE FORMS
The layer-by-layer vertically integrated hydrodynamical
equations were proposed by Professor Walter Hansen (personal
communication). They are analogous to the well-proven single-
layer HN model of the Walter Hansen type.
C, - c2 + Hul(Ulx) ~ Hvl (Vly) = 0
i2 * Hu2 (U2x) +¦ Hv2 (V2 ) = 0
01	+ r /U.lZ < u, . fvi - g ? = k'*'
"ul	lx
02	+ *" Xz * U2 " fM " 9 $ ?1* ' 9 (1" p7» ?2x *
»¦» « r *^11 + VI ^ \n 4.	r = V^y)
'vl
$1 + r /ui^ + vi VI + fUl - g 5ly = k'
V2 + r M* * V2* y2 + m . g £ . g (1. ^ „
where:
?1 - surface elevation;
?2 - deviation of MLD (mixed layer depth) from mean
Ul, VI - u,v components in upper layer
U2, V2 - u.v components in lower layer
r - friction coefficient
- I -

-------
f - Coriolis parameter
g - acceleration of gravity
H - depth
p 1 , p2 - dens i ti es
«(*) ^ k(y) - external forces
There are two inter-dependent continuity equations, one
for each layer. These compute the change of sea level, the
change of depth, and thickness of the layers.
The equations of motion for each layer are vertically
integrated through this given layer. The lower layer is
driven by internal friction and by pressure gradients.
The finite difference forms used for computation are
essentially the same as in W. Hansen's single-layer model,
except for the additional terms dictated by the presence of
a second layer.
_ t + T	~t-T	_ T_ flit	,,t	_ nt
^(n,m,l) ~ ^(n,m,l) I u(n,m,l) (n,m,l) ~ u(n,m-l,l)
t	t	t	t
U(n,m-l,l) + Hv(n-l,m,l) V(n-l,m,l) " Hv(n,m,l)
^(n,m,l) ^ " T ^u(n,m,2) ^(n,m,2) ~ ^u(n,m-l,2)
U(n,m-1,2) + Hv(n-1 ,m,2) V(n-l,m,2) " Hv(n,m,2)
V / n \ }
( n,m,2)
_t + T	_ ~ t - T	_ T rut	lit	_ ijt
?(n,m,2) ?(n,m,2) " I u(n,m,2) (n,m,2) " u(n,m-l,2)
U(n,m-1,2) + Hv(n-l,m,2) + V(n-l,m,2) " Hv(n,m,2)
V(n,m ,2) ^
- 2 -

-------
"(n.m.l) + 2Tf Vi)] /»fn>|n>1)2 + U(JtB(1)2)
o_r	Tfl • r rt^T	rt+T	-i
*TT (n ,m ,1) " £ ic( n ,m ,1) " 5(n+l ,m t.l )J
+9TVt+2T
(n,m)
{l-C2Tr/H*J^2)] /»Jn,.>2)Z - ^>m,2)Z)
'2Tf U(n,m,2) " ^ {p7 Cc(rMi,l) " C2j " 5(^t ,m,2)]}

-------
The computation of 0 and U is done as in the single-layer
model
°(n,m) = a U(n,m) + {U( n-1 ,m) + U(n + 1 ,m) + U(n,m+1)
+ U(n ,m-l)}
^( n,m) = I { U(n,m-1) + U(n+l,m-1) + U(n,m) + u(n+l,m)}
"" *
The computation of V and V are analogous to the correspond-
ing U computations above. The actual depth (and thickness of
the layer) is computed with the following equations:
ut + 2t	, ¦ h	, 1 r rt + T	. _t + T
u ( n,m , 1) u ( n ,m , 1) 2 ?(n,m,l) ^(n,m+l,l)
ut+2t	_ .	, J. j rt + T	t + T	-i
v(n,m,l) ~ v(n,m,l) 2 ?(n,m,l) ^(n+l,m,l)
The following symbols were used in the finite difference
formulas above:
space coordi nates
ti me
components of velocity
indicators of U and V point (location) in
the grid
initial depth (when c = 0)
surface elevation (for the second layer, it
indicates the deviation of the depth of the
layer from its prescribed mean value)
total depth (H = h + 5)
depths at u and v points respectively
x,y
t
U,V
u ,V
h
C
- 4 -

-------
X,Y	components of external forces
g	acceleration of gravity
f	Coriolis parameter
r	friction coefficient (bottom stress for lower
1ayer)
a	coefficient of horizontal eddy viscosity
(also acts as a smoothing coefficient)
n,m	coordinates of the arid point, 1 and 2
indicate the first (surface) or second
(deep) layer
t	half time step
Z	half grid length
p.| ,p2 densities of first and second layer
- 5 -

-------
2. ESSENTIALS OF THE MULTILAYER HN MODEL
This chapter reviews briefly the inputs to the model,
the differences between the single-layer and multilayer model,
and the treatment of open boundaries. As most of the essential
features of this multilayer model are similar to the single-
layer model, the description of the latter should be consulted
(see Part 1 of this series).
2.1 THE STABILITY CRITERIA AND HORIZONTAL VISCOSITY
COEFFICIENT
The same stability criteria as for the single-layer model
applies for the multilayer model:
0. 5£
0.5 At < / 2gH
- 3 max
The Hmax is the maximum depth or thickness of any of the
layers to be expected during the computation. Thus, if
original data indicates that the maximum depth (thickness)
in the lower layer is 520 m, one should add, for example,
80 m for a contingent increase of this depth due to possible
upwards movement of the upper boundary of the layer during
computati ons.
The horizontal viscosity coefficient (also as a smoothing
coefficient) has the same value for the upper layer as in
single layer models (0.98 to 0.99). However, it has been
found necessary to decrease this coefficient for the lower
layers (0.90 to 0.92) in order to smooth out some disturbances
in the interface, especially at the initial state of the com-
putation. Using a higher value for the coefficient results
in longer, more costly computations, as output should be
taken when these initial disturbances have dissipated. The
use of a lower value for the coefficient saves computer time.
The final results are not affected to any great extent.
- 6 -

-------
A lower value of the "bottom friction" coefficient can
be used for the upper layer, as it acts as an internal friction
coefficient. There is no data available as to what this
coefficient should be. Preliminary numerical experiments have
indicated a plausible value of 0.0015.
2.2	THE GRID NET AND THE INPUTS
The same grid net as in single-layer models is used. The
sea and land table and the depths at u and v points are also
introduced as in the single-1ayer model.
¦ _ The depth of the mixed layer (the thickness of the upper
layer) is given with a data statement as an initial flat field.
However, if its topography is known, it can be read-in from
data cards.
The program computes the deviation of the sea surface
elevation from a flat field, and the deviation of the depth of
the upper boundary of the lower layer either from a flat field
or from a prescribed field. The program allows the disappear-
ance of the lower layer (e.g., interception with the bottom)
and/or of the upper layer (e.g., upwelling). The geographical
reorientation of the grid follows the same roles as in the
single 1ayer model.
It has been found that the model behaves best when the
open boundary is chosen such that the iterations are forward-
stepped from the boundary.
2.3	TREATMENT OF OPEN BOUNDARIES
The input of tides at a primary open boundary is done in
every time step for the upper layer the same way as in a single
layer model, using tidal harmonic constants.
Z = Aj cos (o^ t-K-j) + Agcos (agt-Kg) ...
- 7 -

-------
The internal tides must also be introduced into the
second layer. Usually the magnitudes of the internal tidal
waves and their phase lag is unknown. On the basis of avail-
able time series recordings of surface-to-bottom temperature
profiles, it can be generalized that in most areas over the
continental shelf, the tidal internal wave is 4 to 8 times
higher than the surface tidal wave, and lags 1 to 3 hours
behind the surface wave. This preliminary knowledge has been
utilized for the input of internal tides at the open boundary.
The other open boundaries (other than with tidal input)
must be left entirely open for inflow and outflow. Many ways
have been tested to treat these boundaries. It could be
mentioned that the bounds of various Do-loops should be properly
adjusted in treatment of open boundaries in the model, depending
on the orientation of the grid and input boundary. The missing
values of u at the right hand boundary and v at the lower
boundary can be assumed to be the same as at the previous
column (m-1) and row (n-1 ) respectively. The method of
"lubricated wall" for open boundaries is described in paragraph
3.6.
The permanent thermohaline currents are introduced into
the model by prescribing properly precomputed inclinations
(pressure gradient) at the boundaries. The manner of its
application varies with the application area of the model.
This method of input allows the input of the permanent current
independently in two layers; for example, for simulation of
estuarine circulation (i.e., outflow at the surface, inflow in
lower layer). The outflow of rivers must be prescribed at the
gridpoints in or off the river at each time step.
2.4 OTHER GENERAL INFORMATION, INCLUDING NOTES ON
VERIFICATION
The multilayer HN model can be used for the same tasks
(and more) as the single-layer model. The dispersion and
diffusion is handled in this model in the same way as in the
single-layer model, except that it can be treated separately
in two layers.
- 8 -

-------
The two-layer model was originally set up by Professor
Walter Hansen during his visit to Monterey a few years ago.
The original version was programmed by Mr. R. Mann, utilizing
to a large extent, a single-layer program originally written
by Professor J. Slindermann. Since then considerable changes,
improvements, and adaptions have been made by the author.
The program has been tested and compared in a number of
areas. However, the results have not been published previously,
except for some results with a two-layer model for the Strait
of Gibraltar and some with a three-layer model for the Labrador
Sea by Hamilton, Williams and Laevastu (1973).
The two-layer model can also be partly verified with
mareograph data. Contrary to single-layer model verification,
however, total and thorough verification of the two-layer model
requires data on currents and mixed layer depth fluctuations.
A description of the results of the application of a
two-layer model to the New York Bight and its partial verifi-
cation are described in Part 4 of this series.
- 9 -

-------
3. REVIEW OF THE FORTRAN PROGRAM
The program of the two-layer model, presented in Appendix
C,is written in FORTRAN and will run, with very few changes,
on any machine which has a FORTRAN compiler. The necessary
changes are pointed out in the below description of this model.
The main peculiarity of this model (which might vary from
machine to machine) is the reading and writing of extended
core on CDC 6000 and 7000 series machines.
3.1	THE CONTROL PROGRAM
The set-up of the program is simple and convenient. For
example, the dimension, common and equivalence statements are
essentially the same in all subroutines, though all the
variables are not necessarily used in all subroutines. Some
equivalences are used to reduce the core requirements. Some
of the variables, especially in the second layer, are put into
a second common block. An extended core map in the form of
comment cards is given in the control program. The rest of
the control program calls-in various subroutines and counts
for the real time. At the end, it rewinds the tapes used for
various outputs which will be described later. It should be
noted that tapes 4 and 5 are used for a particular plotting
program which is not documented here.
3.2	SUBROUTINE 002
This subroutine reads-in data from control cards as well
as the basic data of depth and sea/land table. The control
cards and their formats are described in Appendix A to this
report. A small Do-loop (DO 702) is for conversion of depth
from fathoms to centimeters in case the depths are read from
a navigation chart in fathoms.
The next series of print statements is for printing the
input data for checking purposes, especially for checking the
digitization of the depth fields. It should be noted that
- 10 -

-------
errors are often made in digitization as well as in the set-
ting of boundaries along the coast. These boundary settings
must be carefully checked, and any necessary corrections made
prior to running the programs. Following the printing, a
Do-loop (DO 900) is used for setting the maximum depth. It
often happens that in the area of interest there is a section
of deep water, the maximum depth of which is used to compute
the required time step.
Do-loops DO 26 and DO 25 are for setting the thickness
of the upper layer. In this particular model, a flat mixed"
layer depth over the whole area is used in the initial state.
However, if the topography of the mixed layer depth for a
given area is known or is constructed on the basis of avail-
able knowledge, it can be read-in with control cards just as
depths are read-in. Finally, Subroutine J02 contains the
setting of some computational parameters, which includes the
zeroing of three fields (AA1, BB2 and CC3). These are used
later to compute the transport through predetermined sections.
3.3 SUBROUTINE J03
Subroutine J03 is used for the initialization of the
fields. In the beginning, it contains a possibility to
read-in Z, U, and V values from a previous computation which
might have been interrupted and the values written on the
tape. It serves the same purpose as the so-called "checkpoint"
feature in CDC machines. The next section is for the initiali-
zation of actual depth which is mainly necessary in case the
Z, U and V values were read from a previous computation. This
subroutine also contains several calls for reading and writing
fields into extended core. This must be changed if the program
is run on a machine which does not have the extended core. At
the end of the subroutine there are Do-loops for zeroing a
number of fields and for initial setting of counters. It
should be noted that if the Z, U, and V values are continued
-ll-

-------
from a previous computation, the counters TIC and TIT must
be set as they were at the end of previous computation. This
depends very much on the set up of the program and will be
described in later paragraphs.
3.4	SUBROUTINE VECTOR
This subroutine is called by Subroutine J04. This serves
to compute the current direction and speed from computed U and
V components. This is done only during the output cycle which
can be hourly, or any desired time interval specified on the
control card.
It should be noted that in this particular program system,
subroutine calls are in FORTRAN II (e.g., ABSF, SQRTF and
others). The so-called "F functions" do not exist in many
newer FORTRAN compilers. There are two possibilities to over-
come this. The first, which offers the greatest savings in
machine time, is to recut the cards leaving out the F at the
end of the compiler subroutines. The second possibility is
to make special function subroutines at the end of the program,
as is the case with this program. Furthermore, it should be
noted that statement #433 contains (at the end) a number +180.
This is the tilt angle of the grid in relation to the north-
south orientation. This problem is discussed in the descrip-
tion of the single-layer model and the value should be changed
according to the tilt taken in a particular program.
3.5	SUBROUTINE J04
This subroutine is for output of the computed fields.
How often it is called in real time is determined by the
parameter T1 which is read-in from an input card in Subroutine
J02. It should be noted that when a new area with new dimen-
sions is printed out, the dimension statements as well as the
Do-loops must be changed correspondingly. First, the water
level variations (astronomic tides and meteorological tides,
or the water level changes caused by winds) are printed. The
- 12 -

-------
second printed field is the deviation of the mixed layer depth
from the mean value which was originally given to this para-
meter. The subsequent two fields give the direction and speed
of the current in the upper and lower layers respectively.
The program contains an option to compute the so-called
"rest current" after a full tidal cycle (ca., 24 hours, 50
minutes). The following Do-loop and few additional statements
are for the preparation of rest currents for the output (using
the same printing statements and formats as in printing the
current direction and speed).
The following little section with Do-loops 9 and 1,000
is for the output of data at selected points on the tape which
will be printed later with Subroutine S08. The positions and
names of the special points from which output is desired is
put in by a control card in J02. Thereafter follows a section
for output of transport through selected sections (in cubic
meters per hour) and zeroing the sections after the output.
Finally, Subroutine J04 includes the printing-out of concen-
tration fields (e.g., pollution) if dispersion and diffusion
computations are made.
3.6 SUBROUTINE J05
This is the main computation subroutine and is called in
at every time step. First it contains the mean Z and Z2L
fields (smoothing of the"height" fields). It should be noted
that, in statement 22, different smoothers are used for upper
and lower layers. A little heavier smoothing in lower layers
has been found necessary, whereas the smoothing coefficient
for the upper layers is normally 0.98 or less, in the lower
layers 0.92 or 0.90 is used. Do-loop 70 is the main Do-loop
for the computation of Z and Z2L values at grid points over
the sea. After the computation of Z and Z2L fields, the outer
open boundaries are set in this particular model, as can be
seen from this boundary setting in the program. The open
- 13 -

-------
boundaries presented in the program in Appendix C are essen-
tially closed for in- and out-flow, but they move up and down
in the fashion of the so-called "lubricated wall." Do-loop
278 is for the input of tides at the primary open boundary.
Preceding this Do-loop, and within the Do-loop itself, a
specific adjustment in this model is made to compute spring
tides of the area. A time lag of 5,400 seconds is also intro-
duced which was required in a double input of tides in another
run of this model but is not documented here. Thereafter
follows a specific adjustment of boundaries to allow some
permanent flow through the areas as known from empirical
measurements. This is a very special input and varies from
one model to another, the description of which will be found
in describing the results of any output from the model run.
As the new elevation values (here, water level and mixed
layer depth changes) are computed, new layer thickness and
depth values must be computed for use in various current
component computations. This is done in Do-loop 85 which is
essentially the same as found in Subroutine J03. Thereafter,
Do-loop 760 contains a check of whether some of the points
have gone dry or the lower or upper layer have vanished due
to sinking of rising of the mixed layer depth.
Before the computation of U and V components, Subroutine
SOI is called if the winds are used as input in the model.
The following section contains the computation of the
U component. It starts with the computation of 0 and V* and
the corresponding values for the lower layers. Before comput-
ing the U and U2L in Do-loops 412 and 1050, the boundaries for
0 and V* are set. It should be noted that U and U2L are
needed from the previous time step for the computation of V.
Thereafter, the boundaries of U and U2L are set, U* and U2L*
computed and the new U and U2L values are transferred from the
intermediate fields (ZM and ZM2) into the proper fields.
- 14 -

-------
Before the computation of V, and V2L are computed, bound-
aries are set and the V components are computed in Do-loops
428 and 1100.
It can be noted that a special additional inflow is set
in every time step in this program at point 26/40, which is
the location of The Narrows on the large scale grid. Such
additional flows can and should be set if outflow from rivers
is known or can be estimated.
After setting the V boundaries, a following section
contains the computation of net flow through predetermined
sections. In this program the net flow is output hourly
during a full tidal cycle. Two approaches which are essen-
tial ly identlcal are used in this section for the computation
of transport. In one approach (sections AA1 and BB2).the
transports are computed by listing the corresponding depth
and current values at desired grid points. The other part uses
a special transport function for computation. This function
is described later in this document. Before coming to. the
computation of dispersion and .diffusion, the rest currents, are
computed and accumulated in special fields during a full tidal
cycle.'
The dispersion and diffusion of pollutants is usually
computed throughout one tidal cycle. Both the instantaneous
source and continuous source can be used and examples of
setting the sources are shown in the program. There are
several different possible approaches to compute the disper-
sion and diffusion. If the grid size is small, the diffusion
can be computed as indicated in this program. On the other
hand, if the grid size is large, the approach used in this
program is unsatisfactory. In this case the central position
of the instantaneous point source is followed and the shape
of the resulting fields is computed in a special program which
uses the output of the current fields from this program and
computes the shape of the pollution blob. At the very end of
- 15 -

-------
the subroutine, another method of computing diffusion (a
simplified Rammig method) is indicated. It gives essentially
the same results as the subroutine above.
3.7	THE TRANSPORT FUNCTION
This function is used to compute transport through the
section as mentioned earlier. It contains several comments
which makes the function easily understandable.
3.8	SUBROUTINE SOI
Subroutine SOI is used for computing the surface wind
current component. It contains the provision to compute the
surface current from surface to bottom (in case of shallow
water) or through the upper layer, or the mixed layer depth
can be used. In the latter case, the last variable in state-
ments 5 and 7 must be changed to HMLD. This is a very simple
but quite satisfactory wind current computation. There are
much more complex and, theoretically, more accurate subroutines
now available in EPRF which use, among other methods, a
prescribed pressure field to obtain the wind.
3.9	SUBROUTINE S08
This subroutine is called at the very end of the program
for output of Z, U and V values in both layers at desired
special points.
3.10	SUBROUTINE J0B6, J0B7 AND THE "F FUNCTIONS"
These latter two subroutines are used for output on a
special plotting tape and could be omitted in case this plot-
ting program is not available. In this case, corresponding
calls for the subroutine should be deleted from the control
program. The "F function" subroutines are indicated at the
very end of the program but could be deleted as described
earlier if the few statements containing these in the program
are repunched.
- 16 -

-------
REFERENCES
Hamilton, G.D., F.R. Williams and T. Laevastu, 1973: Computa-
tion of Tides and Currents With Multilayer Hydrodynamical
Numerical (MHN) Models With Several Open Boundaries.
Proc. of the 1973 Summer Computer Simulation Conference,
Montreal , Canada, 634-639"!
Hansen, W., 1966: The Reproduction of the Motion in the Sea
by Means of Hydrodynamical-Numerical Methods. Mittei1.
Inst. Meeresk., Hamburg, 5:57 pp. and figures.
Hansen, W., Personal communication (Including unpublished
version of the two-layer model development).
Reports in this series:
Part 1: A vertically integrated hydrodynamical-numerical
model (W. Hansen type): Model description and
operati ng/runni ng instructi ons.
Part 2: A multi-layer hydrodynamical-numerical model
(W. Hansen type): Model description and operating/
running instructions. (This report.)
Part 3: 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.
Part 4: 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.
- 17 -

-------
APPENDIXES
APPENDIX A
APPENDIX B
APPENDIX C
-	INPUT PARAMETERS AND FORMATS
-	LIST OF SYMBOLS
-	PROGRAM NYAREA LISTING
- 19 -

-------
APPENDIX A
INPUT PARAMETERS AND FORMATS
THE DATA INPUT CARDS
Card 1 Format 2413
JA, NE, ME, IZE, IUE, KKE , NURU, NG, NS, LI
JA	Indicator, J A = 1 , read initial values; JA=0,
z=0, u=0, v=0
NE,ME Delimeters of entire field (size of the
computational grid)
IZE	Number of points on the open boundary where
tidal heights 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
NG	Number of points at the second open input
boundary
NS	(not used)
LI	Parameter for computation of dispersion and
diffusion of pollutants: LI=0 - no computation;
LI=1 - computation to be performed
Card 2 Format 9F8.3
AFGN, G,	ALPHA, RBETA, CI
AFGN	Problem number
2
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.65)
CI	Wind indicator: 0 - no wind; 1 or 2 - wind
- 20 -

-------
Card 3 Format 9F8.0
DT, TE,	TW, T1, T2, SI, T, (T3)
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., 7200 sec)
SI	Time when wind stops (sec) (e.g., 180000)
T	Time (initialized 0, if previous computations
were made and s, u and v read from tape, give
the time previous computation ended)
(T3)	If plot desired, time when fields are required
for plotting
Card 4 Format 6E12.4.
DL, F, SIGMA, R, R0L, C
1/2 space step (half grid size in cm - e.g.,
185000)
Coriolis parameter (e.g., 8.55x10"®)
Angular velocity of Mg tide (1.4088x10"^)
Friction coefficient (0.003)
Density of the air (e.
-------
First the coordinates of the points at the first open boundary
(input boundary) and therafter the coordinates of special
selected points where special outputs are desired.
Card 6	Format 2413
NA, MA	N and M coordinates at second input boundary
Card 7	Format 2413
NU, MU
NU	N coordinate
MU	M coordinate
Coordinates of the upper right and lower left corners of the
wind field
Card 8
V( I)
Card 9
A (I)
A (1)
A (2)
A (3)
Card 10
HTZ
Card 11
HTU
Card 12
HTV
Cards 13 and 14
Input cards for plotting parameters
Format 10A5
Names of the selected special points
Format 9F8.2
Time when wind starts (sec)
Wi nd speed (m sec"^)
Wind direction (in degrees in computation
coordi nates)
Format 12F6.0 (usually more than 10 cards)
Symbolic depths at the z points
Format 12F6.0 (usually more than 10 cards
Depths at u points
Format 12F6.0 (usually more than 10 cards)
Depths at v points
- 22 -
-------
APPENDIX B
LIST OF SYMBOLS
A	Characteristics of wind field (A1 - time when wind
starts (sec); A2 - wind speed (m sec"'); A3 - wind
direction (degrees, counted clockwise from 0° grid
coordinate)
AAT	A array for computation of transport through given
section (BB2, CCS, etc. serve the same purpose)
AFGN	Program number
ALPHA Smoothing parameter (0.96 to 0.998)
ALP	Smoothing, parameter for second layer (Z2L) (0.88
to 0.90)
ANG>	Compass direction of the current (upper and lower
ANG21	layers)
ARG	Argument of tidal component, (Cos (at - <))
(primary input boundary)
AUS	Austausch coefficient
AT	2 x DT
Ab	F x AT
A3	R x A1
A4	DT/DL
A5	G x A4
B	Intermediate in tidal input
BB2	See AAT
BETA	(1 - ALPHA)/4
BET	(T - ALP)/4
B4	2 * DL
- 23 -
-------
CON	Argument of tidal component, (Cos (at - <))
(secondary input boundary)
CI	Wind indicator; 0 = no wind; 1 = wind
DIR	Compass direction of rest current
DL	Half grid size (cm)
DL2	Same as B4, used in diffusion computation
DT	Half time step
F	Coriolis parameter
G	Acceleration of gravity
GRZ	A3 x JlW2 + ZST2
HGU	Thickness of upper layer at u points
HGU2	Depth of lower layer (thickness of lower layer)
at u poi nts
HGV	Thickness of upper layer at v points
HGV2	Depth of lower layer, at v points
HMLD	Mixed layer depth (not used in this program)
HTL	Initial thickness of upper layer (cm)
HTU	Initial depth at u points (0 = land; -1 = boundary)
HTU2	Water depth (initial) in lower layer, at u points
HTV	Initial depth at v points (0 = land; -1 = boundary)
HTV2	Water depth (initial) in lower layer at v points
HTZ	Sea-land table (0 = land; -1 = boundary; 1 = sea)
I	Index
IUE	Twice the number of wind fields
IZE	Number of boundary points at the primary input
boundary
- 24- -
-------
J	Index
JA	Input indicator
LEN	A counter
M	Field index in x direction
MA	M coordinate of tidal input at the secondary boundary
ME	Total width gf the grid (x direction).
MEH	ME - 1
MU	M coordinate of the wind field delimeter
MZ	M coordinate of primary open boundary (and
M coordinate of special point)
N	Field index in y direction
NA	N coordinate of tidal input at the secondary boundary
NE	Total length of the grid (y direction).
NEH	NE - 1
NG	Number of points at secondary input boundary
NURU	Number of special points
NU	N coordinate of the wind field delimeter
NZ	N coordinate of primary open boundary (and N
coordinate at special point)
PZ	Intermediate parameter in adjustment of boundaries
RAD	Resultant speed for upper layer
RAD2	Resultant speed for lower layer
RB	Intermediate parameter in tidal input
RBETA	Coefficient (fcf geostrophic wind
REST	Resultant speed for rest current (after a full tidal
cycle)
- Z5 -
-------
RH01	Density for top layer
RH02	Density for lower layer
ROL	Density of the air (not used)
S	Concentration field (e.g., pollutant)
SI	End of wind input (sec)
SIGMA	Angular speed of tide (not used)
SSH	Intermediate for computing advection of
SSV	pollutants
T	Time (counter) (sec)
T1	Output interval (sec)
T2	Initial time for first output (sec)
TIC	Index (used in summation of rest current)
TIT	A counter
TRB	Start of rest current output or pollution input
TRE	End of rest current output or pollution input
TW	Time when wind starts (sec)
U	U component of velocity for upper layer
US	Accumulator of u component for computation of
advection of pollutant
USR	Accumulator for computation of rest currents
111	Amplitudes of tidal components at primary
U2	boundary
U2L	U component of velocity for lower layer
U3	Amplitudes of tidal components at secondary
U4	boundary
V	V component of speed in upper layer
VS	Accumulator of v component for computation of
advection of pollutant
- 26 -
-------
VSR
VI
V2
V2L
V	3
V4
V5
V	6
V7
WERTL
WERTL2
WERTO
WERTOL
WERT0L2
WERTOR
WERT0R2
WERTR
WERTR2
WERTU
WERTUL
WERTUL2
WERTUR
WERTUR2
WERTU2
WURTZEL
X
XK

Accumulator for computation of rest current
Names of special points
N and M components of special points
V component of speed in lower layer
Speed and direction in upper layer at
special points
Mixed layer depth at special points (thickness
of upper layer)
Speed and direction in lower layer at special
poi nts
Intermediates for smoothing and averaging
Intermediates for smoothing and averaging
/ZM2 + ZST2
Intermediate tidal parameter at primary open boundary
x component of wind current
- 27 -
-------
Y	Intermediate tidal parameter at secondary open
boundary
YK	y component of wind current
Z	Water level (elevation) (cm)
ZM	Intermediate for smoothing (upper layer)
ZM2	Intermediate for smoothing (lower layer)
ZST	Intermediate for averaging (upper layer)
ZST2	Intermediate for averaging (lower layer)
Z1	Amplitudes of tidal components at primary
Z2	boundary
Z2L	Deviation of mixed layer depth from mean
Z3	Amplitudes of tidal components at secondary
Z4	boundary
- 28 -
-------
APPENDIX C
PROGRAM NYAREA LISTING
PROGRAM NYAREA^INPUT #0UTPUT,TAPfc2,TAPb4,TAPfcialN'UI *TAPb5,TAPfc6)
TWO LAYfcR MODtL.H H,
DIMENSION Z(«2»«7>»U<42,4 7), V ( «2. 4/ ) # *M< 42,47), ZS r (¦»<, 4' > • Hl»U < 42.4
l7),HGV<42,47>,HTZ(4ii,47>,HTU(4*,47>,HTV(42,47>,XM4, YK(42,47>
2,A(12).NZ(90),MZ(V0)»NU(10),MUilU)» ZlUO) i Z2( B(J >» "1('U j« V*( *0 ) »UK
320)#U2(2U)»Xt8),ARG(8),HAD{42,'»7}»ANljil42»4/),V«J(2li)»v4tiio),
4HTL<42,47)#NA(60),MA(60),Y(8),23(8U)^Z4(tt0),U4(8Q),JM'OU),CON(8)
5»V5{20)*V6(2U)»V7(2U)»US»M2»47'J#VS142»47)jS(42,47)
6,ANG2(42,47)#HAD2(4*,47>
7,USR(42,47),VSR{42,47)
tOUlVALENCE(ZM,RAU)»(ZST#ANG», (XK*£2Ui (HGU.MAU2)# ( T*,U*L>,
1(HGV,ANU2), (ZMZ.HTV), (ZS>T2,HTUJ
COMMON P,U,V#ZM,ZM2j42#4/)#ZsiT1«J(42,47).
lHGU2(42»47),HGV2(42,47),HTU2(^#4 7j,MTV2(42,47j,NU(?U),iD(SU)
U0MM0N/LAYER2/V2L#HGU2#HGV2,HTy2iHTV2;N0#MDlAAi,Bi2»tU^
ECS MAP	_
c
ZM
1
c
ZST
1V75
c
XK
3949
c
YK
5923
c
HGU
7897
c
HGV
98*1
c
ZM2
11845
c
ZST2
13819
c
RAD
15793
c
ang
17767
c
RAD2
236»«
c
ANG2
25663
c
Z2L
31585
c
U2L
33559
c
HTV
35533
c
HTU
37J>U7
Extended core map for a
particular model.
REWIND 2
CALL. J02
CALL JO!?
Call jou6
IF(T ) 6 # £, 1
1 CALI J04
29
-------
6 CALL jOt'7
8 CONTINUE
T2=T2+T1
2	T = T +A1
CALL JO^
"IF (T-TE)J,5,5
3	IF(ABsF^T-T2)-A1/2,)1,1,4
4	GO TO 2
5	LEN = b
L0P = 5
CALL J04
CALL JOV7
END FILE: 2
RtWlND 2
fcND FILL 6
REWIND i
CALL SO«
END F I LI 4
REWIND 4
END F I LI" f
REWIND !>
STOP
END
SUBROUTINE JU2
DIMENSION Z<4if,47),U(42,47)lV(42»4/),ZM(42,47),ZSI(^,4/),HliU(42,4
17) ,HGV( 42, 47 ),HTZ(42,47) ,HTU( 42,4/),h| V(42,4/),XK(4^4/j# YK( 42, 4 7 )
2,A(l2)»NZ(90)»MZ(90>|NU(10)|MUil0)#Zi(b0)#Z2(Hu)|VlUu)#V2(20),Ul(
320),u2(20)»X(8)#AHG(8),hAD(42|47),AnU(42#47)»V (20 ) , V&(2U ) , V7 (2U ) , US< 42, 47l, VS I 42, 47) , Si( 42» 4/j
6»ANG2(42»47)#HAD2(42»47)
7,USR(42,4/),VSR(42,47)
E0UIVALENCE(ZM>RAD)#(ZSr#ANG)ilXK»Z2L)i(HGU»HAU2>»(T^,y2L)»
1(HGV,ANG2), (ZM'2,HTV)#(ZST2, HTJ)
COMMON Z,U.V,ZM,ZST.HGU#HGV,HT,BB6# BB7,BB8,CC4,CCS>,CU6, «C/,CCd
DIMENSION Z2L(42»47)#U2L(4^#47i»V<»L(4k!i4/)#ZMi;i42»40»23T2(42»47),
lHGU2 ( 42, 4 7 ),HGV2( 42, 4 7 ),HTU2f4i!,4/'), flfV2(42,4/),NU(3u) »iD(SU>
COMMON/LAYER2/V2L#HGU2»HGV2#HTU2»H I V^iND,MD, AAi,Bd^,oC,J
DATA HTL1/25UU,/
1U0 FORMAT(6F10.2)
101	FORMAT(10A5)
102	FORMAT(*413)
103	F0RMAT(9F«.3)
104	FORMAT(9Fb,0)
10b F0RMAT(6El2.4)
106 F 0 R M A T ( ?• 2 F- 6 . 0 )
— Control program. JOB6 and JOB7
are for a particular plotting
program. Tapes 4 and 5 are for
the same plotting program.
- 30 -
-------
107 F0RMAT(12F6.1>
112	F0RMAT(!>X,19HHALF TIME STEP UT=,F5,o»SH SbC,29X,2UH1ALF SPACE STEP
1 DL=,E11,4.4H CM///)	"
113	F0RMAT(5X,27HSM00TH1NG PA"AMbTbR ALPHA" , K6 , 3 # 29*t i .i-"1 <* ICT I ON CObFF
llCENT R*,bll,4///)	.......
114	K0RMAT(SX,22HC0RI0L1S PARAMETEH F=#Ell14#/H 1/StJ,XOX,49MANGULAR
lVELOCITY OF THfc M2 TIDE SI GMAn, Ell# 4#/H l/SEC////>) ~
115	FORMAT(PX|13HH1ND IS CALM,/////)
116	F 0RMAT(9X,21HTHE WIND BbGlNS AFTbRiF/,0,6H SbU ,, 19*» 1 J^A JK ¦ ObNSl T
1Y=,E11.4,VH ti/CM**d///)
117	F0RMAT(9X#23HRfeDUCTi0N FACTOR WE T As , ^ 3, 26X, 2SHIHA CueF^ I (J I bNT L
1AMBDA = #fc11,4»11H (CONSTANT)/////)
120	F0RMAT(t>X#54M< N, M) COORUlNATbS OK THb I POINTS ON T^t U-'EN bOUNUA«
1Y//J>X,12(1H(, 12,1H,, I2,lH),2X)//5x;i2(iH(, 12,1H, , 12,Ihj,2X)//)
2//5X »12 (1H ( #l2#lH,»l2»lH),2X)/>5X»l2UH(#l2»lH#»U,lM)»iX)
3//5x»9(lH(,12,1H,#I2,1H),2X)
121	F0RMAT(lHl,l4xi36HSYM80HC WATtR DbPfH AT THb L P JI S///BX, 14 I 8/>
122	F0RMAT(///8X#4I8//)
23	F0HMAT(lHl,14X,33HMATbR DEPTH AT THE U2 POINTS CC1)///BX*1418/)
123	F0RMAT(1H1.14X,33HWATER DEPTH AJ THE U1 POINTS (Cl>///BX, 14.18/)
24	FORMAT(lHl»14X,33HWATfcR DEPTH AT	THb V2 POINTS CC1)<'/4*j1418/)
124	F0RMAT(lHl,14X,33HWATER DEPTH Aj	THb VI POINTS  (NZ(I>, MZ(1 *39), I = l, NUKU)	control cards
READ (3,102) CNACI),MA{i)#
READ	( 3' i.02 ) (N U ( I) ,MU( 1), I»l, 1 Jb)
READ	(3,101)C VI(I),lal«NURU)
READ (3'100> (A( I ),1¦!,KKE)
READ	(3,106) ((HTZ(N|M)#H«i»liE)#N«l#Nfc)
READ	(3,106) ((HTU(N,M),Mal,MEi,Nsi,Nbj
READ	(3,106) ((HTV(N,M),Msl,MEi,N*i,Nb)
UO 9967 Nal,41
HTZ(N,i>«0.
HTZ(N,2)8-2.
Special input boundary setting.
9967 CONTINUE
700 j
I
I
I
i F(HTU(«",M))/01, 701*703
Conversion of depths In fathoms to cm.
703	HTU(N.M>sHTU(N;M)«1B»,2U1
701	IF(HTV(N»M))702»702»7U4
704	HtV(N,M^bHtv(N,M)*1B9,2U1
702	CONTINUE		_J
- 31 -
-------
PRINT 11*,DT»UL		
PRINT 1L3,ALPHA,R
PRINT 114,F,SIGMA
IF(C1)4,3#4
3	PRINT 115
GO TO 5
4	PRINT 116,TW,HOI,
PRINT 117, RBfc T A » C
PRINT 133, (A( I ), Isl.KKE)
5	PRINT 120,,1=1,UE>
8 PRINT 121f(N»Nslf14)
PRINT 130,(N,(HTZ(N,M),M = 1,14),N = l,Nt)
PRINT 222,,2B)	|
(HTU(N,M),M = li>,28;,N = l,Nt)
N = 29,42 )	|
(HTU(N,M),Ms2V,4 2),N=l,Nb)
N = 43 , 4 / )	' |
(HTU(N,M),M=4j,^7>,N=l,Nt)
N=1,14)
(HTV(N#M),M=l,14),N=l,Nfc)
Ncl5,2b>	|
(HTV(N,M),M=15,28),Nsl,Nfc)
Ns29,4
-------
10
12
14
2U1
lb
17
19
21
202
25
6*2
656
660
661
HTU2(N#y) = -1,0
GO TO 15
HT(J2(N,M) s U.U
GO TO lb
HTU2(N,M) = HTU(N,M) • hTL(N'»M>
IF(HTU2(N,M>)201,201,15
HTU2 ( N»f1) = 0 .
IF CHTV
HTV2(N,r> = U.O
GO TO 2»
HTV2(N«M) S H TV(N» M) - HTL(N # M J
IF (HTV2 >2112,202,25
HTV2(N»H)=0.
CONTINUE
FORMAT(///bX,14Ib/)
F ORMAT(ib,14Pa;0)
FORMATtI9,5I»/>
FORMAT(?5,'5FB,0)
HETAs(l.-ALPHA)/4,
LOP = l
T4 = T1
a1b2,*DT
A2sF*Al
A3cR*Al
A4sDT/DL
A5sG*A4'
C^sC*Al*lU00U .
AA1=Q.
682=0,
CC3=0,
ne!h=ne-i
MEHsME-1
NEHH=NE"2
MEHHsME-2
RETURN
END
If)
2# A(12)
320)
4HTL(42i
5,V5(20:.	. , 	
6,ANG2(42»47)»RAD2(42|47)
7,USR(42,47),VSR{42,47>
fcQUIVALEfiCE, (HUU,RAU2)» ( riWU^U ,
1(HGV,AN!32), (ZM2,HTV),(ZbT2,HTU)
COMMON Z, U,Vr 7m 7CT-ur.ii.Uf;u.ur/
1V1,
21
31
fcQUI VALENCE , (HUU,RAU2)» (TIWU2U ,
L(HGV,AN!32), (ZM2,HTV),(ZbT2,HTU)
COMMON Z,U,V,ZM,ZST,HGU#HGV,HT*,HTU,HTV,XK,YK,A,N*,'UiNJ,MU#Zl,*2i
LV1#V2»U1#U2. X,ARGlVi#V4«HTL#N4»MA,Y,2i,Z4,U3,U4,i: JN
?,V5,V6,V7,US» VS#S,USR,Vt>R
S»F,G,Sl!fMAfALPHA#R»ROL#HBETA,C#DL»t)T, T # T1, T2 ,Nfc, M = , 11, K*E, WE, IUE,
- 33 -
-------
4TW,BETA»Al,A2/A3,A4,Ab#Cl,C3»Sl#NUHU,NURV»JA»NfcH, ifciiNfcHH#MfcHH,HOC
5,NG,T3,5TV#ETV,	T INCV,STHf ETH, TlNCHi THAN# T4, AF"GN,Ticii-fcNiNb
6,TRB	» TRt		
7,AA2»AA3,AA4,BB3,BB^|BB5,BB6,Hb7»BtiB,CC4»CC5,CC6,JC/iCC0
DIMENSION Z2L(42,47),U2L(42>47) , VifL ( » 4 7 ) ,ZM2<42>~4/ J i 2 i ( 4 2 » 4 / ) i
lHG|j2(42#47),H[iV2(42»4 7)»HTu2C4i ,M = l,Mb j ,N = jL,.Nb j
READ (3/100) <(V(N,M),M = l,ME),N = l,Nfcj
GO TO 11
7 CONTINUE		
IF(JA)11,410,11
41U UO 12 N* 1# NE
DO 12 M=1,ME
IF(HTZ(N,M))2U1,200.2U1
2U0 Z(N,M)=0.
Z2L(N,M) = 0,U
GO TO 2L: 2
*Ul Z(N » M ) = 0 , ^
Z2L(N»M) = 0,2
2 U 2 IF(HTU(N,M))2U3,203,204
203 U(N#H) =0 ,
U2L(N,M) = 0,U
GO TO 2U!j
2 U 4 U(N,M)=0 .2
U2L(N,M{ = 0,2
20b IF(HTV(N,M ))2U6,2U6,207
206 V(N/M)=U.
V2L(N,M) = 0,U
GO TO If
*07 v(N,H)=0 ,2
V2L(N» M' = 0,2
12 CONTINUE
CAUL WRlTtC(HTU,3/5U7ilV74J
CALL RE: DfcC(ZST2,13BlV,l974)
CALL WRITbC(HTV,3i?533ilV74)
CALL READbC(^M2,1184t>fiy74)
DO 208 N=l#Nfc
UO 208 M=l,Mb
IF(HTZ(N,M))2OV,2U9,210
209	ZM(N,M)=0,
ZM2(N» M ) = 0,U
ZST < N» H ^ =U ,
ZST2(N» M ) = 0.0
GO TO 20B
210	ZM(N» M) = U,2
ZM2(N» M ? = 0,2
ZSTM?=U.2
Z S T 2(N,M) = u.2
— Reading of Z, U and V
values from previous
computation.
Initialization of fields (O or
a small value over the sea).
- 34 -
-------
(Z2LCM»M) ~ 4*U)/2.
I
)/2.
2U8 CONTINUE
11 CONTINUE
80 0 CAL.L WRITfcC5,5,4
5	HGU(N#M)sO,
GO TO 5U1
4 HGU(N,M)8HTL(N,M)*03,502
502	HGU2(N» M) = HTU2(N,M) +
GO TO 13
503	HGU2CN.M) = U.0
13	ir(HTV("'#M) )b,8,6
8 HGV=0.
GO TO 5"3
6	HGV
-------
BB5
Bh>6
BB7
UB8
CC4
CCjj
CC6
CC7
CC8
0 .
O:
0 .
0;
0.
0;
0.
o7
Setting some counters.
CALL WRlTtC(ANG2,25663,1974)
CALL WRITbC=ANG(N,M)M.S60./6l2831B53[)/3j
Subroutine for computation
of current direction and
speed from U and V
components.
U(V(N.M)>143,142,142
142	IF(U(N,M))50U»433,433
143	If (U(N.f') >501»&01,144
144	GO TO 502
5UQ ANG(N»M)sl80,-ANG(NiM)
GO TO 4J3
5 01 ANG(N,M)=18Q.«-ANG(N,M>
GO TO 433
502 ANG(N,M)=^60,-ANG(N,M)
433 ANG(N,M)s<90|-ANG(M,M))*180.
5 U 3 IF(ANG(N,M))14*,432,432
145	ANG(N,M)=360.+ANG(N,M>
GO TO 503
432 L'ONT I NUt
RETURN
END
SUBROUTINb J04
UI MENS I ON Z<^^,4/ ) ,U< 42, 47) , V( 42*4/ ) #ZM( 42,4/) #ZS r ( 4 / ) ,HiiU( 42. 4
- 36 -
-------
17),HGV<42#47),HTZ(42,47),HTU(42,47J,HTV#Ui(
J20),U2(20),X(«)»ARG(8),KAD(4i»>57)»ANG(42#4 7),VJ(20)» V^UO)»
4HTL(42*:7),NA(60) «MA(60)» Y(8)|23(SU)*^4(tiU)iU^(BQ),^(ttl))iCUN(B)
5,V5<20>#V6<2U>7v7(2U),Ui»<«2#47),VSl4*,47j,S(42,47j
6,ANG2<42,47),RAD2<42,47>
7,USR(42i47),VSR(42,47)
EQUIVALENCE(ZM#RAU). (ZSr, ANG), IXK,^2L), (HGU.HAU2)# ( TiS,u#
KHGV.ANU2), (ZM2,HTV), (ZbT2#HfU>	' "
COMMON Z,U.V»ZM#ZST#HGU«'HGV»HTi.HTU#Ml V»XK,YK,Af N*,H.NJ,MU»Z1»*2»
1V1#V2.U1»U2*X#ARG,V3,V4#HTL#NA»MA#Y,43#Z4,U3#U4,CJN
2#V5,V6»V7|US#VS#S»USR»Vt>R
3,F.G,SlKMA,ALPHA,R,K0L»WBETA,C»DL#UT#T#Ti»T2#Nfc,M=,Ie#K*E» UE# IUE«
4TW.BETA#Al,A2,A3,A4,A5,C1,C3#SJ,NJNU»NURV#JA,NbH,iE"#NfcHH»MbHH,HQC
5#NG,T3,STV,ETV;TINCV,STH,ETH, T INCH i THAN, T4 , Af GN, T I C» Ut=N '$ N!>
6,TRB»	TR|	" "			
7,AA2,AA3,AA4,bB3,BB4,BB5,BB6,Bd7,dH8,CC4,CC5,CC6,w'C/#CCa
DI MENS I N Z2L(42#4 7),U2L(42#47i,V2L(42;47),ZM2l42;4/)»"2ST2(42,4/).
1HGU2(42#47),HGV2(42.#47) ,HTU2C42,4i),Hrv2(42,47),NU{^U)|rtD(5U)
C0MM0N/UAYER2/V2L»HGU2«HGV2#HTU2;HiV2iND#MU,AAi,Bd2»tC3
CALL WR»TbC,2b)
PRINT 500 #(N# (Z(N,M),M8l5,2H),Nsl,l\|E)
PRINT 5-2'» (N«Ns29,42>
PRINT 500,(N,(Z(N,M),Ms29,42),Nsl,NE)
PRINT 501#(N,N«43#47)
PRINT 500#(N,(Z(N,M),Ms29,42),Nsl,ME)
PRINT 5U3#(N#(Z(N,M),M*43#47),Nsl,NE)
PRINT 500,(N#(Z(N#M)#M#29#42),Nal,NEj
-Printing of Z values
494	F0RMAT(lHl,l4X,64HRtSULTANT CUWRbN rrTPEfcU (CM/Skw) ANU DlRfcCUUN
1(DEG,) AFJEH T=F9tU,/5X#I7#i319/)'
495	FORMAT(//5>X# 17,419/)
496	FORMAT((I5,14(F4I0,1H,,M,0)/))
497	FORMAT((I5,5(K4.0,1H# #F4,0)/))
498	Format(//5x#17#1319/)
499	F0RMAT(1H1.14X,25HWATER LEVfcL 1 CM> AKbR T«,~ 9, 0/3X# W# 131 8)
500	Format(i5,i4FBto/)
501	F0RMAT(//5X«i/#4i8/)
502	F0RMATC/5X# 17,1318/)
503	FORMAT(I5#5F8,0/)
PRINT 599,T»(N,Nsl#14)		
PRINT 500#(N,
- 37 -
-------
Printing of Z2L
values.
PRINT	5»0 2 , 
PRINT	50Q,(N»(Z2L(N»tf)#M=lb#28J,N=l»Nb)
PRINT	502,iNsi»lMb)
PRINT	5(J 11 (N,Ns43,4/>
PRINT	5'3>(N»(Z2L(N»M),M=43»47J »N=l»Nfc)	
554 F0RMAT(lRl,14X,64HDEEP LAY". CURRENT * SFEFO (CM/SEu) ANU D I H fc C T I UN
1(DEG, ) AFTER T = ,F9 , 0,/5X,j7#10 I 9/)
559 F 0RMAT(1H1,14X,2JHMLD DEPTH CM" AFfEK T = , F V . 2/5 X , I 7 . A 3 .1 0 )
1634 CALL VECTOR(HTZ,U,V,RaD, AN(j,NE»ME) "
PRINT 494,T.(N,Nsl,14>
Nt)
PRINT 4'6|(N,(RAD(N,M),ANG(N',MJ,M = 1,14),N = 1;NE>
PRINT 49B,(N,N=15,2«)	j
PRINT 4*6, (N, (RAD(N.tf) ,ANG(N,MJ ,M = lt>.2«) #N:{,
PRINT 49B,(N,Ns29,42)	|
PRINT'496,(N,('RAD(N,M),ANG(NfM)#M = 29#4 2),Nsl,Nt)
PRINT 4V5,(N,Ns43,4/)	|
PRINT 49/,(N,(RAU(N,M)fANG(NlMJ,M=43#47),Nsl#Nfc)
IF(LOP--0>1605,1U,163b
1635 CALL VECTOR(HTZ,U2L#V2L#RAD2»ANG2,Nb,Mb)
PRINT 5541T »(N,Nsl,14)
PRINT 496, (N#(RAD2(N>M),AMG2(\)/M)iM = 1,14)iN = 1,NE)
PRINT 498,163U,1U,1&30 ~		" 		1
1630	IF/(0,U1*TIC>
U(N*M)=USR(N»M)/(0,U1*T1C>
1633 CONTINUE
L0P=1Q
GO TO 1634
1U UO 9 Is'.NURU
N = NZ ( I ~>S9 )
M=MZ(I~39)
V2(I)=ZlN,M)
V3(I)=RAU(N < M )
V4 { I ) = Ar|G ( N » M )
9 CONTINUE
Printing of directions
and speed of upper
and lower layers.
Preparation of rest currents
for output.
Output of data from selected points
on tape to be printed later with
WRITE (2) T, ( VZ< I ) » I sliNURU) # < V0( J ) , lsl,NURU) » < V4( I ) , lU.NURU)
DOIOOOI=1,NURU
N = NZ(I +39 )
M = MZ( I -t-39)
Vb( I )=Z'U(N.N)
V6(I)=RAD2(N,M)
V 7(I) =ANG2(N,M)
10U0 CONTINUE
subroutine S08
- 38 -
-------
4611
4612
WRITE (?)T, ( V»( I )i I«lf NURU)« (V0( J) • I ai#fJuHU)» ( V7( I ) #|81#NURU)
TRB=5460U, " —'—1
TRE=144UQU, . . ,
lF
5lUl#Bti6,BB7#dBB
5101,CC4,CC!>,C;C<)
9*oi,ccy;ccB
Output of transport through
selected sections (m^/hour)
S101 FOPMAT(lHl,//10Xt.26HTHANS'POHTf SECTION A, M3s , Fiu u,//28X, 2HB,,
15X,F10.U»//28X,2HC, #5X, HO. 0 ,"/?)"
4620
AA1 = 0
•
AA2 s
0.
A A3 =
0i
A A4 =
0 .
BB2 = 0
•
BBS =
0.
BB4 =
0.
BBS =
0.
BB6 :
0.
BB7 =
0.
BSe =
o;
CC3 = 0
• i
CC4 =
0.
CC5 =
0.
CC6 =
0.
CC7 =
0.
CC8 =

	 Setting of sections to 0
CONTINUE		
Print 614,t» (n,nsi,.i4)	|
PRINT 5U0, (N, (S(N#M>,Msl,.14>,N = l.,Nb)
PRINT 502,{N,Nal5,ZB) " ~ |
39 -
-------
PRINT	500,,Ms43,4 7 ),Nel,Nb|_
Printing of concentration
fields (eg. pollution).
614 FORMAT<1H1,14X,37HCUNCENTRaTI3IM UF LUWtR LAYER AFfE* T = # F V . U , / / bX ,
118,1518 /) '
452 RETURN
fcND
SUBROUTINE J 1)5
DIMENSION Z(42,47),U(42»47),V(42i4/)iZM(42,4/),ZSr(^,4/),HbU(4t?,4
17),HGV(42,47 ),HTZ<4ii,47),HTU(43,47j,rtTV(42,47),XKi42i'»7j#YK(42,47)
2,A(liJ).NZ(90)»MZ(90>#NU(10J#Mlj(lU)»Zi(ao>.22(Bu)»ViUu').Vi!(^0)»Ul(
32 0),U2(20),X(8),AnG
b, Vb(20)» V6(2U),V7(2U),Ut>(42,4 7V, VS(4^,47) ,S(42,4/)
6# ANG2(42»47) #RAD2<42#47)
7jUSR(42»4 7),VSR(42,47)
fcOUI VALENCE ( ZM,RAD) # (ZS1 » ANG) , ( XK» C2L), (HGU,RAD2 ) » ( Ti\,U«!L) »
1 ( HGV, ANG2 ) , ( ZM2, HTV) , ( Zt>T2, H1UJ
COMMON Z,U,V.ZM,ZST#HGU»HGV»HU,HTU,H!V#XK,YK,A,N*,',U#NJ,MU«Zl»Z2#
1V1,V2,U1»U2.X,ARG,v^,V4»HtL,NAVma,Y,Z5,Z4.U3,U4,CJN
2,	V5,V6»V7,US,VS.S.USR,VbR
3,F.G.SIGMA,ALPHA,R,HOL,HdfcTA,C»	DL.UT,T,Tl,T2,Nb,M=, It,K*E, IZE,IUE,
4TW,BETA»Al#A2»A3,A4,A5#Cl,C3iSI,NJHU#NUHV«JA»NbH, lE^tHH, MbHH, POC
5, NG, T3, STV,ETV, TlNCV,STH,ETH, jIncH; THA.N, T4, AFGN, T ».C»l-cN,"Nb
6# TRB,TRh
7, AA2, AAJ, AA4,BB3,bB4,BB5,BB6, HU7,BU8,CC4,CC5,CC6, -'C{ #CCc5
U I MENS I UN Z2LU2#47) ,U2L<42#47) ,V2L.(42i 4 / >#ZM*(42#4')#23t*(42,4/>•
lHGu2(42>47),HGv2(42,47),Hru2(42,47),RTV2t42,47),NU(3u)i^U(!?U)
C0MM0N/UAYER2/V2L.HGU2»HGV2#HTU2#HrVii#ND»MU,AAl,Bj2«oCj
Data RH.l,RH02/l,U2b,l.U29/
CJ05 COMP
CALL WRITtC(HAD,157V3,lV74)
CALL WRl.TbCCANG, 17767,1*74)
CALL HRITEC(HAD2,23089,1974)
CALL WR,l,TEC(ANG2,25063, 1974)
Call REAUfcc(ZM;i,i9/4)
CALL READtC(ZST,197.5,19/4)
CALL REAUbCCHGU,789/,19/4)
Call RE'«UfcC(HGv,9a7iii9/4>
CALL WRlTbC(HTV,355«S3, 1V 7 4 )
CALL READbC(ZM2,11845,1974)
ALP = U , 90	~~~
BET=(1. " ALP ) /4 ,
DO 3110 Nsl.NE
DO 3110 M = 1# Mb
ZST(N,M)=o.
3110 ZM2(N,M)=U.
DO 4 70 Ns2,NfcH
DO 470 M=2»MbH
IfCHTZ«N,M))710,470,75
710 IF(HTZ(N#M)+2.)470#75,4/0
- 40 -
-------
Computation of Z and Z2L
("smoothing").
75 IF(1-N)28 j 6 j 20
28	IF(HTZ>57,6,57
57	WERT0 = Z(N«1,M> '
WERT02 » Z2L
HERT02 f Z2L(N,M)
7	JF(N£-N>1B,9#1B
18 IF(HTZ(N*1,M))58,9,58
58	WERTU=Z
WERTU2 ' Z2L(N,M)
10 IF(1-M)29,12,29
29	IF(HTZ(,i#M-l))59#12»59
59	WERTL=Z
GO TO 13
12 WfcRTL=Z(N,M)
WERTL2 = Z2L(N#M)
U IF(ME-M)l/,15il7
V 1F(HTZ(^»M+1))60,15,6U
60	WERTR=Z(N#M*1)
WERTR2 = Z2L(N,M+1)
GO TO 22
15 WERTR=Z(N,M) .
WERTR2 = Z2L a ALP •Z2t
470 CONTINUE		
UO 70 N = 2»NE 	Computation of Z and Z2L.
DO 70 M=3,ME
604 IF(HTZ(N,M))^U,70,149
149 ZCN»M)aZM(N»M)"A4*(HGU(N.M)*U(N,M)-HtiUCN#M-l)*U(NtMwl)*r,GV(N-l,M)*
1V(N-1#M)-HGV(N,M)*V o ZM2 (N | M) - A4» < HGU21 N, rl) *U2L (N, M ) - HGJ2t)*U2l(N,
1M-1) ~ HGV2(N-1,M)*V2U(N-1#M) - HGV25»M)
Z(
Z2L<2»M>sZ2L<«S#M)
Z2L(i,M)sZ2L'(2,M)
Z2LC2,M)aZ2l(4,M)
570 CONTINUE "	
TDel04400,*T
TG»TD-5400,
TGT=TG-TD
UO 78 N~2, 40
RBcN
a=/38t
Setting of boundaries for Z and Z21.
Input of tides at one open boundary.
- 41 -
-------
Z, y *CU:> (
1, 000145416* (	TGT ) - 4 , / 65 6 ) ~ / . / ( , 0 6 U 13 7 §5 « • ( T LUb *
2	)-4 . 2124 )+6 . 6*C0S ( , UOOO'29o8* ( I D*ti* TUT ) -1, 034 j ) +6 , j *CUS (
3	, 000067585* ( TD*B* TGT >-2,4256)
UN,2) - U,7 U *Z(N,2)
78 CONTINUE
DO 380 f' = 2,4U
PZ=FLOAi(N)
Z(N,2)=Z(N,2)-4,/PZ
3b0 CONTINUE
DO 390 s2,44
PZsFLOAT(M)
Z(2,M)=Z(2,M)-0.1*Pl
Z(l,M)=4(2,M)-0.0b*HZ
390 CONTINUE
TGT
Specific adjustaments to boundaries.
CALL WRiTfcC(ZM2,11845,iV74)
CALL WRITfcC(ZST2,13ttl9,lV74)
CALL RECUfcC(HTU,375U7,1V74)
CALL REAUbC(HTV#35533'1974)
20U3 DO 8b N = 1»NEH
DO 85 M = 1» MEH
IF (HTU(f!»M) )76,76,81
HGU(N.M)=HTL(N»M)*(Z(N,M)*^(n,h*i))/£,
IF  >/7,
HI
76
y 3
n
8 7
82
84
85
761
762
7bJ
764
765
7 66
767
76 0
161
162
103
3112
HGU2(N»r>=HTU2(N,M)*(Z2L(N,M;*Z.2L85,85,84
HGV2(N,fi)=HTV2(N,M)*(Z2L(N,M)*^2L(N+l,M) )/2
CONTINUE
DO 760 rj = 1, N b
DO 760 s 1, Mb
IF(HGU(N,M)-3U,)761,761,762
HGU(N,M)=0.
IF(HGU2(NiM)-30,>76.3, 763, 764
HGU2(N, )=0.
IF(HGV(N,M)-4U,)765,765»766
HG V(N» M i sU,
IF(HGV2(N,M)-3(J. )76/»76/»760
HGV2(N,H)=0.
CONTINUE
lF1U3» 162,1~62
CALL SOI
CALL WRiTtC(HTU#3/5U7,lV74)
CALL WRITEC(HTV,35533,1974)
CALL READEC(ZST2,13819,1974)
CALL READEC(ZM2,11845,1974)
DO 3112 N=l,Nb
DO 3112 M = 1» Mb
ZST2(N,M) = 0 .
UO 512 :=2,NtH
Computations of
actual depths and
layer thickness.
Calling of wind subroutine (SOI)
- 42 -
-------
Computation of U
(smoothing).
and U2L
DO £>12 M«2,MEH
IF(HUU(N,M) )!>12,512i413
413 IF(1-N)12«,1U6,12«
IF(HQU(N#H)>S12,512i413
.413 IF <1-N)128,106,12d
128	iF(HGU(N-l,M) >1*7,106,1&7
157 WERTOsU(N-l.M)'
WtRT02= UZl
WERT02 " U2L(N,M)
107	lF(Nfc-N>118,109,118
118 IF(HGU(N*1,H>)i58,109,158
lt>8 WERTU = U(N*1,M)
WERTU2 e U2L(N*1,M >
GO TO 110
109	WERTU=U >159,112,1I>9
159	WERTL=U(N,M-1)*
WERTL2 r U2L(N|M.1>
GO TO 113
112	WERTL=U(N#M>
WERTL2 P U2L,160
160	WERTR=U '
W^RTR2 * U2U(N»M#1>
GO TO 2i
US WERTR = U(N|M>
WERTR2 s U2L61,32,61
61	WERTOLsV(N-l#M)
WERT0L2 s V2L(N-1,M>
.IF(HGV(N-1,M*1)>62,35,62
62	WERT0RsV(N-l,M*lj
WERT0R2 = V2LCN-1#M*1)
GO TO 36
32 WERTOL=V(N-liM*l)
WERT0L2 » V2L«N-1|M*1)
35	WERT0RsWkRT0U
WERT0R2 » WERT0L2
36	IF CHGV C N,M)>63,38,63
63	WERTUL=V(N,M)
WERTUL2 ¦ V2L(N, M)
IF(HGV(N,M*1)>64,41,64
64	WERTURSV(N,M*1>
WERTUR2 B V2U
GO TO 42
38 WERTUL=V(N,M*1>
J
* *
Computation of V and V2L .
- 4-3 -
-------
41
42
51 2
WERTUL2 s V2L(N,M*1>
WERTURsWERTUL
WERTUR2 = WEKTUL2
ZST(N»M) = (WEHT0L*WEWTUL4'WERTU(?*w[:RTUHil/4|
ZST2(N»M> = (WERT 0L2 ~ W E R t U L 2 ~ WtRfuH2 '~ WtKTUK*!) t 4.n
CONTlNUr	' -	I ~
	 Setting of boundaries.
4/3
513
51/
518
94
14
DO 473 Ms2,Mt
ZM2(NE,n)=ZM2(NEM#M)
ZM (NE»M)=ZM (NEH,M)
ZST(NE»f)=ZST(N E HiM)
ZST2(NE,M)=ZST2(NEH,M>
CONTINUE
CALL WRITEC(Z2L,31585,1974)
CALL REAUECtXK,3949,1974)
UO 412 Ns2,Nfc
DO 412 M=2,MbH
I K { HGU 412,412,513
lF(ZM(N,M)*ZM(N,M))412,!>17,t)10
IF(ZST(N,M)*ZST (N,M) >412,412,518
WURZEL=SQRTF*ZM(N»M>*2ST(N#M)*^ST)
URZ=A3*"UWZfcL
IF(HGUCN.M)-UHZ>94,V4,14
GRZ=HGU(N,M)
Computation of U and U2L.
ZM(N.M) = (l.-liKZ/HGU(N#M>)*ZM{\#M)+A2*^!>T(N,M)-A5*lZl!Yil1*lW(N,M))
1*XK(N,M»
412 CONTINUE
CALL WRITfcC(XK,3949,1974)
CALL PEADEC(Z2L»415«5,1974)
UO 1050 N=2,Nt
DO 1050 M = 2» MEH
IF(HTU2(N»M)) i05U,1050»10U5
IF(ZM2(N,M)#ZM2(N,M)) 1U50,1U1U,1015
IF(ZST2(N,M>*ZST2) 1050,lU50,ioi?
WURZEL = SQRTF **2 >
GRZ = A«3*WURZtU
IF(HGU2(N,M)'GRZ) 1U3U » 1030 »1025
GRZ=HGU2(N,M)"
ZM2WVI2) - A 5 M i ,
10 U 5
1010
1015
1030
1025
A2*ZS'12(	- A5*>
1050 CONTINUt
UO 4 /4 f =2 , NJt
ZN2(N#Mt)eZM2(N.MEH>
ZM(N,ME>=ZM(N,MEH>
CONTINUE
DO 528 N=1,NEH
DO 528 M=2,MtH
IF(HGV(N#M)>528,528,429
IF(HGU( #M-lj >§5,44,65
WERTOL=ll(N,M*l j
WERT0L2 = U2L(N,M-1>
IF(HGU(N*1,M-1> >66,46, 66
66 WERTUL=U(N*1,M-1)
474
429
65
Setting of boundaries.
- 44 -
-------
44
46
47
67
68
*2
>3
528
* *
Computation of U and U2L .
163
170
171
172
l7^
174
430
WERTUL2 = U2L(N+1,M-1)
GO TO 47
WERT0L = U{N«.1,M-i)
WERT0L2 s U2L(N+1#M-1J
WERTUL=WtRTOL
WERTUL2 s WEKT0L2
IF(HGU >67,49,6/
WERTOR=U(N,M)
WERT0R2 a U2L>68,52,68
WERTURSU430,430,163
IF(HGUC ,M)>170,170,171
U(N,H)s0.
C?0 TO 172
U(N>M)sZM(N»M>
|F(HGU2(N,M)>173,173,174
U2L(N»M>so ,
GO TO 4^0
U2L=ZM2(N,M>
continue
Transfer of
into proper
new U and
fields
U2L
432
229
257
206
20/
218
2*0
209
CALL WRlTfcC(HTU,37507,1974)
CALL READEC(ZST2,13B19,1974>
UO 431 NsJJ, NEH
DO 431 Ms2,MEH
IF >431,431,432
IF CI-M>220,206,228
IF (HGV (' -1,M> >257,206,257
WERTO=V
WERTQ2 F V2L(N-1,M>
GO TO 207
WERTO=V
WERT02 = V2L(N,M>
IFIN6-N)218,2U9.,218
IF (HGV (r!*l,M> >258, 209,258
WERTU=V(N*1,M>
WERTU2 = V2L(N*1,M>
GO TO 210
WERTU=V(N,M>
Computation of
(smoothing).
V and V2L
- 45 -
-------
WERTU2 = V2L(N,M)
1F(1-M)229,212, 229
IF(HGV(N,M-i))259,2l2,259
WfcRTL=V(N,M-l)
WERTL2 = V2L(N,M-D
^0 TO 213
WERTL=V(N,M)
WERTL2 s V2L(N,M)
IF(ME-M)217,215,217
iFtHGvtN.M*!))260,2l5(2bO
WERTR=V(N,M*1)
WtRTR2 = V2L(N,M*i)
GO TO 2?
WERTR=V*bErA«(WfcR10 + WtRTU*WtRTL*WtRT-
Z S T 2 ( N # E > = ZST2(N,MbHJ_
CONTlNUf
CALL WRlTtC(U2L.3J5&9,1V74)
CALL READbCfYK,5923,1974)
DO 428 .', = 2,NbH
DO 428. M=2,Mb
210
229
259
212
213
21/
260
215
25
475
Setting of boundaries.
529
617
618
IF(HGV(N»M)) 428,428 , 529
IF(ZM(N#M)»ZM(N,-M)) 42b,017,618
3 U 0
415
Computation of V.
A2*ZST)^
lF(ZST(N,M)*ZST(N,M))42tt,4^8,618
WURZELsSQRTF(ZM(N#M>*ZM(N,M)»ZST(N#MMZST>
liRZ = A3*":URZEL
IF(HGV(N,M)-GRZ)300#300»415
GRZ=HGVlN|M)
V(N»M) = <1 ,-GRZ/HGV(N,M)MZM(N,M)
lYK(N.M)			[
428 CONTINUE
V	( 26, 40 ) = V ( 26, 4 0 ) + 1,3 	setti Qf
CALL WRITbC(YK,5923,1974)	B v
CALL RE:DfcC(U2L»3«J559, 1974 ) _
DO 1100 N=2»NEH
DO 1100 M=2»Mb
IF(HGV2(N,M)) 1100,1100*1055
IF(ZM2(N»M)*ZM2(N,M)) 1100,10611,10 65
IF(ZST2(N,M)*ZST2(N,M)) 110 0,1100.1063
WURZkL = SQRTF (ZM2(N,M)**2 ~ Z 5T2 (IN > M > **2 )
GRZ = A3«WURZbL	I- Computation of V2L.
1F'(HGV2(N,M) - GRZ) 108U,10&U,10/5
GRZ=HGV2(N,M)	I
V	2l(N » M ) = >*ZM2(N»M) - A2«ZbT2(M,1) • A5MRH01/
1RH02)*(Z(N,M) - Z ( M ~ 1» M ) ) - A 5 * ( ; . - KHOl/RHU2 > * U2L I \l« 1) " Z2L
-------
S(17»37) s 1UOU0,
S(22,9) a 1UUOO,
S(27«18) b 1UUOO,
S(32,7) s 1UUOO,
S(34,27) c 1UQD0i
S(38,6) s 1UOOO,
s(ie#i5) s mooo,	
S(26» 40 ) B S(26» 4U) "~ 3,
4412	01«DL*DU
UU2=2.*^L
AUS=5.0b4
4413	B2 = <4,U«DT*AUS >/Bl
414 04s2.*DL
4415 B5 = (A'JS * UT)/bl
416 Y1*U.
S(26,40) a S (26,40) ~ 3,
; COMPUTATION Of ~DISPfcRSIU.N anU UlKFJSiUN
4417 UO 3000 N°2,NtH
DO 3000 Ms2#MtH
IF(HTZCN,M)>3O00,3OU0,2Ul
201	IF(HTU(!;» H)) 2 U 2,2 0 2., 2.0 4
202	SSH=0.
GO TO 225
204	IF(U(N,flj >205.205,211
205	IF2207,22117,2208
22U7	SSH = S CN #N)/{DT•ABS(U( N# M ) ) J
US(N # M)*0.
Go TO 225
2208	SSH = S(N»H*1>//DL2
GO TO 225
211 IF(US2213,2210#2210
2210	IF(S(N.p-1))22il,2211,2212
2211	SSH = S (N » M ) / (D.T*ABS ( U.( N#;M ).) )
US(M#M>*0,
GO TO 225
2212	SSH = S(N«n-l)/(DT*ABb(U(N,M)).
US(N»M)s0 t
GO TO 2e5
2213	SSH=(S(!:#M)-b(N,M-l ) )/-DL2
225 IF{HTV(N#MJ>2214,2214,2215
2214	IF (HTU(N,M) )30 00, 30UO ,22.30
2230 SSV=0.
GO TO 442B
2215	fF(V(N.*)>2216,2216,2221
2216	IF(VS(N,M>*84)2220,2217,2220
2217	JF(S(N-1,M))2218,2218,2219
2218	SSVsS(N,M)/(DT*ABS(V(N,M)))
VS(N#M)r0,
GO TO 442B
continuous source.
Computation of dispersion
and diffusion, of pollutants.
- 49 -
-------
*219 SSV=S(N-1,M)/(UT*ABS(V(N,M)))
VS(N,M)=U,
GO TO 4*28
2220	SSV=(S(N»M)-S(N-1,M>>/DL2
GO TO 4428
2221	IF ( VS( N,M)-B4)222b,2222, 222b
2222	IF(S(N*',M)>2223,2223,2224
2224 SSV = S(N,M)/(UT*ABS(V(N,M)) )
VS(N,M)"0,
GO TO 4428
2224	SSV = S(N*1»M)/(UT*ABS>(V(N»M) ) )
VS(N»M>=0,
GO TO 4 42®
2225	SSV=(S(N»M)-S(N*1»M) )/DL2
4428	S(N,M)=i!(IM,«)-b2*S(N,M)*Bb*(Ji(N-l,ri)+S(N*l,M)>!»(NiM*i)*»(N,M«.l)
1-4,*S(N.M))
2-(DT*ABsF(U2L(N»M))«SSH)-(L)T*AbSMV2l.(N#M))*SSV)*Yl
(S(N-l,M-l	>*b(N*l»M*l> )
4429	IF (S(N» M) ) 24 29,2429,4Qu0
2429 S ( N« M ) = U.
.3000 CONTINUE		
C 415 B5 = 2.0*84
C 416 DO 422 N=2,NbH
C 417 DO 422 Ms2,MbH
C 418 IF(HGU(N,M))42U,420,419
C 4i9 S(N.Mj=5(N,M)-Al*U(N,M>*(S(N',H*l)/a4-S(N#M-l)/b4)*B^*(b(N,M-l)/2l-
C 420 IF(HGV(N,M))422,422,421
C lS(N#M)*S(N,M*lj/2l)-»B2*(S(IM-l,M)/2t-iiCN#M)*S(N*l,1)/^,)
C 421 S(N.M)= (N,M)-Al*V(N,M)*(S(N*l»M)/S4-bi(N-l,M)/B4)*B^i(b(N-l,M)/2,-
C 2S(N«M)+S(N+l«M)/2, ) *B2* C S ( N, M -1)/2 , - b ( N» M ) *S ( N , M* 1) / it,')
C 422 CONTIMUE
904 CONTINU!
11 RETURN
bND
F_ UNCTION TRANS(FV,FH#N,M,K£V,NUM)
D I MENS I :-.N F V ( 4 2, 4 7 ) , F H ( 4 2, 4 7 )
c
c
c
c
c
c
c
c
f-V = u OR v COMPONENT 01- VbLUClTY
FH r HEIGHT ON u OR V GHIU"
N = INITIAL VALUE OF N , « . , M I N
M = INITIAL VALUE Oh M,,.,M I N
KEY = 1...SECTION IS ALONG CONSTANT N
KEY = r...SECTION IS ALONG CONSTANT M
NUM = NUMBER"OF POINTS ON THE SECTION
C
TOT = 0
Ml = M-i
N1 = N-l
Function for computation of
transport through sections.
C
GO TO(100,20U) KEY
C
100 CONTINUE
- 50 -
-------
c THIS IS A SECTION ALONG CONSTANT N . ,,V VfcLOC IT V* ^Scu.vn* <]..?<* 16
DO 150 K s 1,NUM	"		 -
TOT s TOT ~ FV{N;Ml>Kr*>H(»#Mi*K)
150 CONTINUE
TRANS = TOT
RETURN
C
200 CONTINUE
C "THIS IS A SECTION ALONG CONSTANT M..lt U VbLOCjTy U^tJ
UO 250K a X,NUH	" " ~
TOT « TOT ~ >'V(Nl*K.M)#t-H(Nl+<.H>
250 CONTINUE
TRANS = TOT
RETURN
end
SUBROUTINE SU1
U I MENS I UN Z(42,47 ),U(42,47),v<«2,4V)UH<42,4/)(Zsi <«*,4/),HUU<42,4
1?>»HQV(42#47 j,HTZ'(4'a,47).HTU(45f 4/J,nrV(4Z,.47),XK«4-, NU < 10 )# MU < 10 >, Zi ( t»0 ). L2 ( bu ) ,*1 i*0 j,V2( 2q>, Ul<
320 ).U2< 20 ),XC»> # A«G *B>,MADc42/«»-7 J #ANUC42#4 7 J, ViCZli) »v4(^0)#
4HTL(42#:7),NA(60)#MA(60)»Y(B),Z3(au)i^4(aO j,U4(Bo J, Jjj«OU)iCUN«8)
5,V&	< 20),V6(2U)^ V7 < 2U >,Ub(42,47 ),VSI 42, 4 7j,s(42,4 7 )
6,ANG2(4*,47)»RAU2<42,47)
7,USR(4Z.4/),VSsR(42,47)
kQUlVALENC£(ZM,RAD),(^Sf»ANG), 1 XK Wdl), ( HUU, HAU2) # ( Ti\-,y2i..)»
1(HGV.ANR2).(ZM2,HTV),(ZST2.H1U)
COMMON Z,U,V.ZM,ZST#HGU«HGV.HT2,Hry,rilI V, XK,YK,A,N<,*£« N J, MU>.Zl, ^2#
1V1# V2,U1»U2>X,ARG# V«S#V41 HTL»nA»MA» Z4#U3|U4,C.JN
2,V5.	V6,V7,US,VS,S,USR, Vi>R
3,F,G,SICMA,ALPHA,R,ROL,HBETA,C,DL»UT,T,T1,T2,Nfc,M	= ,,K*E» Ut,IUE,
4TW,BETA,Al,A2,A3,A4,Ab,Cl,C3.Sl ,NUKU,NURV. JA,NtH,1t'i7NbHH,MtHH,POC
S,NG.T3.i' TV,ETV,TINCV,STH,EfH, T INCH. THAN, T4, AfliN, J 1C« L t(V» N!>
6 %TRB #TRfc
7|AA2,AA3#AA4»yB3lBB4#bB5»HB6#Bb7#dbafCC4#CC5lCC6twC{#CC? .
01 MENS I ON Z2L(42,47)#U2L(42# 47),V2L(42#4/),ZM2142»4/)•L»T2(42i4/)i
lHGU2(42J47)#HGV2(42,4;),HTu2l42,4/),Htv2t42l4/)#N«J(3U)»'iDl>U)
COMMON/LAYER2/V2L«HUU2,HGV2#HT02,Hf V'^;n13,MU# AAi,B^2»i:^3
C SOI WIND
CALL WRITEC(Z2L,3l5bS,iV74T
CALL WR;'TEC0 HMLDs2000,4.8,*Z(N,M)
51	IF
-------
61 HMLD=HGV{N,M)	1
7	YK(N,M)=C3*A(2)*#2*SlN(A(3)*U,Ul/43)/rlGV(N,M)
8	CONTINUE
GO TO 1?
1U DO 11 N»1,NE
UO 11 M*1,ME
XK(N,M)»0,
11	YK(N » M ) = Q i		
12	call writecrARG(8),KAU<42,57)fANli(42#4/)lvJ(i;U).v4{(42#47J,VSl42#4/),b(42|4/l
6»ANG2(42»47>»HAU2(42#47)
7,USR(42,47),VSR(42,47)
EQUI VALENCE(ZM,RAD) ,(ZS1» ANG),IXK # Z2L ) »(HGU,KAU2)»( Ti\,u*L> •
l{HGVfANGZ),(ZMZ,HTV)#(ZST2.HrU'
COMMON Z,U» V,ZM,ZST,HGU,HGV,HTi,HTU,HTV,XK, YK, A, N^ , •'U, NJ, mu, Zl ,/¦>!.
1V1,V2,U1,U2-, X,ARG,V3, V4#HTL.NA#'ma, Y,Z3,Z4,U3,U4,CJN
2, V5, V6, V7,US, VS,S,USR, VSR
3»F",G, SIGMA, ALPHA, H,ROL»HBETA,C»DL,UT« r#Tl.T21Nt,M=,lc>K^t,Ufc,lUb,
4TW,BETA» Al, A2» A3, A4# Ab,Cl,C3»Sl iNUKU.NURV# JA,NtH, ifc^NbriH, MfcHH, HOC
5,	NG, T3, STV.ETV, TINCV, STH,ETH, T INCH," THAN, T4, AhGN, T »C«~l£n." Nb
6,TRB	» TRt
7#AA2,AA3,AA4,BB3,BB^;lU(5U)
C0MM0N/(-AYER2/V2Li HUU2. HGV2, HTU2» H I V2," ND, MD, A A l, H*2» oC3
C SOB OUTPUT FOR SPECIAL POINTS
131	FORMAT(lHi,5X»li5HWATER LtVEL lCM),ShU CURRENT SP = b^ (t V^bC ), D I RE
1CTI0N (DEG).MLD (M),DfeEP CURRENT SPEbj (CM/SEC) A VQ "uIHcC TI UN (UEG
2)///7X#5HTIMfc.4XilO(2X,A7)V
132	FORMAT((6X,5H(SEC),4X«1U(2X«I3#1H,,I3)//))
133	FORMATUll.0i*X,10FV,2)
140 FORMAT((15X,1U(F4,0,1H,,F4,0>//15X»1UK9.2//15X,10(F5»U , 1H, A 4,0)/)
1)
134	H ORM A T (5 X > 14 HE ND OF PROGRAM) 	^ Output of values from special
NURV = 0	points.
PRINT 131,,I=1,NUNU)
PRINT 132#(NZ(I*39)»MZ(1 ~ 3 9)»J=1,NURU)
1 READ (2 ) T , (V2(Ij,1=1,NURU), ( V<3 ( I ) »I a 1,NUKU)» 
READ (,c ) T , ( V 5 ( I ) , IS1,NURU ),(V£( I ),I=1,NUHU>, < V 7(I),J = i,NURU)
PRINT 133,T,(V2(I),1*1,NURU)
PRINT 1*0, (V3( I ) , V 4 ( I ), V 5(I),V6( I ) # V / t I ), I51,NURU)
- 52 -
-------
IF(Tt-T >2,2,X
2 PRINT 134
RETURN
end
SUBROUTINE JOB6
DIMENSION Z<4ii747),U<42#47), V< 42,4/), ZM( 42,47 ),ZS I ( , 4 / ), HliU < 4*, 4
17)#HGv(42,47) ,HTZ (42,47 ),MjU( 42,47 ),HTV( 42,47 ):,XK( 42; 47 j, rK(4i!, 47)
2»A(12)«NZ(90)'t1Z(9o)*NU(10)f MJ(lii)«ZX(tlO)*Z2(Bu)*Vii,u2(20),X(B),ArG{8)#hAD(42,57),An9(42,47),V^(20)#V4(^0),
4HTL(42«"7),NA(60),MA(60)#Y(B),23(8u)»^4(dO),U4(8o),J^tOU),CUN(8)
5.V5<20),V6(2U),V7<2U),US<42,	47')',yS<4i>,47),S<42,47i "
6#ANG2(42,47),RAD2<42,47) "	"
7,USR<42#47),VSR(42,47)
EQUIVALENCE(ZM,RAD),(ZS T, ANG),IXK,Z2L),(HGU,HAU2)#(TK,U*L>:'
l(HGv,ANG2),(ZM2,HTV),(ZbT2,HTj)
COMMON Z,U,V,Zh,ZST#HBU»HQV#MTi»MTU,HTV,XK,YK# A,N4,'l(t#NU,MU»Zl,Z2,
lVl,V2,UX»U2,X,ARG#V»#V4,HTL,NAVMA,r,23,Z4,U3,U4,CJN "
2»V5,V6,V7#US,VS,S,USR,VSR
3,F ,G.Sl!fMA#ALPHA,R,ROl,,RBETA,C#DL#UT, \, T1, T2, Nfc, Mc, 11, k*E , UE ,\ UE,
4TW,BETA»A1,A2»A3,A4,A&,C1,C3»SI,NUKU,NURV»JA,NfcH,1E"iNfc*H,MbHH,POC
S#NG,T3,^TV,EtV,TlNCV,STH,EfH,TiNCH;TKAN»T4, AKliN, T IC'LeNiNb
6,TRB«TRE
7#AA2,AA3,AA4,HB3,BB4,BB5,BB6,3b7,aB8,UC4,CC5,CU6,oC/,CUd
DIMENSION Z2L(42,47),U2U(42,47), V2L(<»ii 47),ZM^(42M' > •'2STi* ( 42, 4 / ) #
lHGu2(42,47),HGv2(42,47),Hru2(4i,47)#RTV2(42,47)#nJ(?U),^D(50)
COMMON/LAYER2/V2L,HGU2,HGV2,HTU2#HtV^;ND,MD, AAi,B52»<-^«J
ENTRY TO READ VECTOR ~ HEIGHf UAHUS
ETH = ETC
ifOO READ	(3, 215) iTV11 T V# T INCV	\ A plotting sub-
201	PRINT 216,STV»ETV,T1NCV	x routine.
202	READ	(3,*l!>)!>TH,tTC»Tl "¦JUH
203	PRINT 216,STH,ETC,TINCH
21b FORMAT(CF10•0)
216 FORMAT(10X,3F10,0)
PAUSE 3 IS FOR MOUNTING OF Pl_0iT1NU TAPE
218 WRITE	<*) HTZ,Vl,NUKU,ME,Mfc
RETURN
END
SUBROUTINE JOB7
U 11
17)
2,
320 )
4HTL ( 42 »«»' ; , « i « o * #-n " » " u » » ¦ \ ~ - v *-*•'- • z w • w - - •-».
5,V5(20)»V6(20)»v7(ZU),US(42,47),VSt4Z,47),S(42,47)
6,ANG2(42,47),RAD2(4ii#47)
7,USR(42,47),VSR(42,47)
EQUI VALENCE(ZM,RAU),C 2S1,ANG),UK,*2L),(HGU,HAU2)#(Ti\,y^L>#
ItHGV.ANrZ),(ZM2,HTV)»(ZST2,HiJ)
COMMON Z,U, V#ZM,ZST»HGU«HGV,HTj£,HTO,HTV#XK, YK, A,NNJ,MU,Zl,^2»
1V1#V2,U1»U2,X#ARG,V3,V4,HTL,,NA»mA#Y#23,Z4,u3,U4,CJN
2,V5,V6#V7,US,VS,S,USR,VSR
3,F,G#SIHMA,	ALPHA ,H,K0L#WBETA,C»DL»yT,1,Tl»T2,Nb,M = ,[t,K|
-------
4TW,BETA#Al,A2.A3,A4iA!5,Cl,C3»SI , NUHU, NUR V , J A , NfcH, tfci,NfcHH,MfcHH,HOC
5,NG, T3,;.TV,ETV, T INC	t'JNJUH; THAN, T4,A^n# TI C # LtN," Nb
6 » TRB # TRE
7,AA2,AA3,AA4fBB3,BB4,BB:>#HB£.,4ri7tfciab,CC4,CC!>,Ci;6,JC'#CCeJ
UI MENS I HN Z 2 L ( 4 2 i 4 7 ) ,U2L< *2« 7 J , V ( 4 2 » 4 7 > , ZM2\ 42 m M » L ST 2 ( 42 » 4 / ) .
lHGU2(42»47),HGv2(42,4/),Hr!j^(^,4/J>'1[V2(zl2,4/),NLJ{3u),^U(5U)
COMMON/tA YER2/V2L,HliU2> HGV2. HTU?, H IV2»NU*MU,AAi,bi2'^3
C	ENTRY TO WRITE DATA ON lAPfc Y OK VfcCTUK A N U WATbR ^E1^' PL^T
230	VEC=HT=0. '
231	IF(ABSF(T-STV)-Tl/2,)232»232,234
232	STV=STV*TINCV
233	yEC = l,		^ A plotting subroutine.
234	lF(ABSF(T-STH)-Tl/2, > 233 , 2ib . 2-17
235	STH = STH*T INCH
236	HTsl.
2 3 7 IF'(HT + VEC)23tt#239,238
238	WRITE	(4)T,VEC,HT.RAD|ANU#*
WRITE (S)T,VbC(HT,RAD2,ANG2#ZZL
239	CONTINUl:
RETURN
END
FUNCTION ASINF(X)
AS 1 Nf = AS I N C X )
return
end
FUNCTION COSF(X)
cosf = toscx)
RETURN
END
FUNCTION SQRTF(X)
SQRTF = SURT(X)
RETURN
END
Function sinmx)
sinf = t> i N ( x)
RETURN
END
FUNCTION ABSF(X)
ABSF = ABS(X >
RETURN
- 54 -
-------