•
I (
-c-•
•
•
Delft UniversityDepartment of Civil Engineeringof Technology Group of Fluid Mechanics•
•
•
•
1mplementation (I) of the numericalshallow water wave hindeast model H1SWA
•
N. Booij T.H.C. Herbers L.H. Holthuijsen Report No. 6-84•
•
•
•
•
PROJECT REPORT Delft University of Technology Departaent of Civil Engineering Group of Fluid MechanicsProject title GEOMOR wave model (HISWA)
•
Project description
Development of a two-dimensional model to hindcast spectral wave parameters in an estuary with tidal flats on the basis of bottomtopography, current and wind data
•
•
CustomerRijkswaterstaat
Deltadienst, Afdeling Kustonderzoek van Alkemadelaan 400
2597 AT THE HAGUE, the Netherlands
•
represented by J. v. MarleProject leader dr. ir. L.H. Holthuijsen
•
work carried out by T.H.C. Herbers dr. ir. N. Booijdr. ir. L.H. Holthuijsen
•
ConclusionThe propagation of waves has been imp le-mented (except diffraction) and tested. The input program is under design. The output program will be developed in two versions.
•
•
Status of report confidential, progress report•
City/date: Delft, June 14, 1984•
•
TABLE OF CONTENTSpage
•
1. Introduction2. Status of programs 2
•
3. Tests for wave propagation 44. Dissipation of wave energy due to currents 10
• References 13
Appendix I - User's manual of HISWA
•
Appendix 11 - Subroutines of COMPUAppendix 111 - System documentation of HISWA
•
•
•
•
•
•
•
- 1-•
1. Introduction
•
This report describes the progresswater wave hindcast model HISWA since the last progress report~n the development of the shallow(Herbers et al.,March 21, 1984).•
The status of the development of the three programs of HISWA ~s g~ven~n section 2. Tests of the second of these programs (COMPU) are given in section 3 (wave propagation). The original design of this program is supplemented with a method to account for wave energy dissipation ~n opposing currents. It is described in section 4.•
•
•
•
•
•
•
•
•
2-•
2. Status of programsThe HISWA model is planned to consist of three programs :
•
a. PREP, which is an input model, planned to be a modified version of CREDIZ of Rijkswaterstaatb. COMPU, which is a computation program to determine wave propagation and wave generation and dissipation, under development
•
c. OUTP, which is an output program, planned to be a modified version of the UITV program of CREDIZ but possibly a separately developed program (see below).•
The modification of PREP of CREDIZ to suit HISWA has been started by changing the text of the user's manual of CREDIZ (Appendix I). The modi-fications in the source text of PREP will follow in a straightforward manner from this revised manual.•
•
Most of the subroutines of COMPU have been coded and tested (see Appen-dix 11). The subroutines that have been tested relate primarily to the retrieval of input data and the propagation of the waves (excluding dif-fraction). The subroutines controlling diffraction, wave generation and dissipation will be implemented when the propagation tests are concluded. Modifications 1n the system documentation produce a new version of this documentation which is provided in Appendix 111.
•
•
The modifications of UITV have started in the sense that the nature of these modifications has been studied. The conclusion is that the re-quired modifications are extensive and possibly outside the .scope of developing HISWA. Considering the corresponding time-consuming revision of UITV we suggest a two-track approach which, on the one hand will pro-vide print and plot output of HISWA during the test stage, and, on the
other hand, will hopefully result in the modified version of UITV inte-grated in the final HISWA model. Two output programs will therefore be developed
•
•
a. OUTI, a preliminary print and plot .subroutine to be used during the tests of HISWA
•
3-•
b. OUT2, a modified version of UITV, planned to be integrated 1n the final version of HISWA.•
If OUT2 is not finalizedthe output program of HISWA.by DecemberThis would be weIl within31, 1984, OUT 1 will be releasedthe terms of as agreement of the development of HISWA.•
•
•
•
•
•
•
•
•
•
4-•
3. Tests for wave propagationThe program COMPU has been used to compute several situations of wave propagation
current
(case la and Ib) (case lc)
(case 2a, band c)
•
shoalingrefraction
•
For some of these situations the numerical results are compared withanalytical solutions. The situations are•
a. the main wave direction perpendicular to the bottom contours and a narrow directional distribution (cos16 8)b. the main wave direction perpendicular to the bot tom contours and a wide directional distribution (cos2 8)
•
c. the main direction at an incidence angle of 22.50 relative to the normal of the bottom contours.•
a. a following current increasing with travel distance b. an opposing current increasing with travel distance c. a slanting current increasing with travel distance.•
•
The slope of the plane beach is 1:120. Both the following and opposing current increase from 0 mIs at the "upwind" boundary of the grid to
I mIs at the "downwind" boundary. The slanting current increases from
o
mIs at the right hand boundary to 4 mIs at the left hand boundary. The difference in direction with the waves is 22,50•Allother relevant numerical information for the tests is given in table I and illustrations of these conditions are given in fig. I.
•
•
The results of the computations are given in figs. 2 through 4 for the cross-section y
=
1200 m parallel to the x-axis. The agreement with the analytical solutions in cases la, band c is good and the results seem to be realistic in the cases 2a, band c.•
5-•
% General
computational grid size 960 m x 2400 m
mesh size 40 m x ]00 m
depth grid size 1080 mx 2700 m
mesh size ]20 mx 300 m
minimum direction of propagation % -450 (-600 for case ]b)
maximum direction of propagation +450 (+600 for case ]b)
directional increment ]00 (13~0 for case ]b) % Boundar~ conditions (long side of grid)
Plane beach
---~ ]a ]b ]c
significant wave height (m) 0.40 0.39 0.38 average wave frequency (Hz) 0.5 0.5 0.5
main wave direction 00 00 22.50
directional width a = ,%X% ]6 2 ]6
depth (m) 3 3 3
Current
---case 2a 2b 2c
significant wave height (m) 0.40 0.40 0.40 - 0.35%% average wave frequency (Hz) 0.5 0.5 0.5
main wave direction 00 00 22.50
directional width a,= %%% ]6 ]6 ]6
depth (m) ].5 1.5 ].5
•
•
•
•
•
•
Table 1. Wave, bottom and current conditions for the propagation tests.% directions relative to x-axis, which is aligned with the shortest
axis of the grid
%% wave action
=
constant=
0.0032 m2/s•
ae
directional distribution is appliedtt% a cos
•
•
•
•
6-•
•
1. Plane beaeh ( area 960 m x 2400 m )( beach -0 m
•
d=3 m la 0° lb 0° Ic 22.5° a=
16 a=
2 a = 16•
2. Currents ( area 960 m x 2400 m )r
T
r
t
J
.
!
T
t f'T
t
t
t,
J,T
f =1.5 m'"
......
..,....
,
•
2a () 0° 2b (~ 0° 2e 22.5° a=
16 a=
16 a=
16•
Fig. 1 Review of wave bottom and eurrent conditions in the propa-gation tests of COMPU .
•
•
•
•
•
- 7 -H (m) s•
-- - x case la { shoaling lb { shoaling lc (refraction ) cf' cosïse distribution ) IIcos2 e distribution )
f.
r-model o case"
5
f I/
v! ~,.."'
""
.)lr .-:'"-
"
-:.r: - ~ - ~ - ..._ ~
-d-ll:a.
-"'-
_X-:;;..:':-:f5"" -. • - -. ~. =:.--=----=-
.
~.
_-="-=><---c:::tII' -.. .. --X- _.><:.... _ ~ x :.._"'"~""" -"(JIJ~ ~~_c ~_- _~ CJt
..__
_
( lines : linear wavetheory )
-
--
-
.
case•
.0L---~---~---L---~-- ~ ~~ ~~ 500 X{m) Cle(rad) 0 0 0 0 0 0.
5
0 0 0 0 0 0 0 0 0 • ~ x....
.
..
0 lC l{ Je .-)t A.
0
500 X{m) :e(rad) 4-
-
.
.
--~
-
-
.
-
-....
_
_.
.
'
.
-
.-
..
-
..
...
..
--
..
-
--.
......: -.
.
... _
...
....,-...
_
~
.
...
,...
..
.3 ( case 1c )•
.... ...'...
.2 ,.
.
\•
,",
• 1 '\ \, •·
o~---L---~---L---~---~---L-
~~
~
500 X{m) ( case 1a , 1b )Fig. 2 Significant wave height Hs ' directional spreading Cle and main
wave direction
ë
for plane beach (shoaling and refraction, no current).•
•
•
•
•
•
•
•
•
•
•
•
H (m) s
• case 2a ( following current ) X case 2b ( opposing current )
o case 2c ( slanting current )
I I X , I
,
x
, , 1.0 ,x ., "'x '" ...l(.. x-_,__x1(.---
-
--
.
-
---
---.
---.
---
.
- --
.
- -- ----0 - - -0 --_0__ -0- --0 -- -0 - - -0 - --=~-:.~::..&-:..:.-::&-:.:::.a-=- --.- - -..,- L": ~ -: - 7-- _0 0__ -0 0 _--0 __ -0 _•-0- __ 0 0 . --
--.
...
-__x x- __x- __, ~- _-x - -"x - _x---)<. - - X
-
.
---
.
---
.
---
.
-
-
-
.
-
-
-
.
-
--
.
--
-
.
00 Io
500 X (m) .5.
---
.
-
--
.
-
-
...
· - - -. - - -. - - -. - - -::::.:; -= : ~ : -= ~ ~:.=-~:~
:
-
_
-:
:
==
~: ~
~:
:
:=:
:
::
~
:
:
=~
=
:
~
= ~
:
=
=
=
.
: ~:
:
~
:
~
0- - _0- -- 0 0 Co: .,.._....'-~==.::.~ ...:~--= ....-:..~=
=-Jl::..:::"x - --)<- - - x.- - -x- - - ...- _~-- - x- -- x-_ -x - _-)(-- _,1(__ -x., X _ -x.-__><__ ->0. - -x- -_-<- - - •o
500 1000Fig. 3 Significant wave height Hsand directional spreading
oe
for following • opposing•
4 mis•
•
•
wave direction•
•
o
2 mis•
•
•
•
•
C>
o
mis•
9-_..,.-
---
__.
---0.36 0.35---
0.35---
--
---
--
-__..-
__..-
---
-
--0.37TH
s 0.36 0.35---
---
--
----
--X(m)•
- 10-•
4. Dissipation of wave energy due to currents•
It is assumed that in a situation with currents opposite to the direction of wave propagation the degree of wave breaking increases at an increasing rate. In the present model a source term is built to account for this effect.
In the subroutine ITWN a critical wave frequency w is computed, which is c
the maximum frequency for which a solution exists for the wave number k (eq. 1, fig. 5.).
•
-+ -+ {gk tanh (kd)}%
(1) w - k • U - = 0 -+with U current vector
•
-+k wave number vector d water depth
•
•
gk tanh (kd) }i
•
k•
Fig. 5 Critical wave frequency Wc•
It is assumed that the frequencies for which w > w loose energy. The ceffect of this loss on the total energy E and the average frequency w
o 0
is model led with two relaxation models
•
dEodt = (2)
dbr
•
11-•
dw 0=
-(w - w ) (3) dt T2 0 00 dbr•
00 in which E=
f E(w) dw 0 0 00 w A(w) dw I 00 w=
f f A(w) dw 0 0 0 E=
~CE(W) dw 00 0•
w w Woo=
f Cw A(w) dw Ifc A(w) dw 0 0 A=
E/cr , cr=
{gk tanh (k~}1/2•
The time constants Tl and T2 have to be determined empirically. Initially
-1
it is assumed that Tl = T2 = W
o
•
The values of (E - E ) and of (w -w) are related directly to wand w
o 00 0 00 c 0
if a standard form of the spectrum is assumed. We choose here the spec-trum, as presented by Kruseman (I976, fig. 6).
E(w)
•
E max E(w)=
0 E max E max ( wlw - b ) ~ -5 ( wlw ) m for w < bw = m for bw < w < m=
wm for w>
w == ···m•
bw w m m wcFig. 6 Kruseman spectrum ( 1976 )
•
•
If the mean wave action frequency
w w
and w by f ~w E(w) dw I f CE(w)
00 0 0
stages in GONO, Sanders, 1976) then w
o
The (analytical) expressions for (E
-o -4 (E - E )
=
0.3335 (w Iw) E for o 00 c 0 0 00 w is approximated by f w E(w) dw/E 000dwand b
=
0.7 (as for most of the growth=1.17w. m Eoo) and (w - w ) o 00 are w c > 0.8547 wo
•
•
•
- 12-•
= {-5.704(w /w )2 + 6.825(w /w ) - 1.042} E for 0.5983 w<
w ~ coc 0 0 0 c 0.8547 w o E for w ~ 0.5983 w 0 c 0 1 - 0.4444 (w /w )-3 and (w - w ) w-
{ co} w for w>
0.8547 w•
o 00 o. - 0.6250 (w /w )-4 0 C 0 C 0 0.6667 (w /w )3 - 0.5983=
w _ { C 0 o (w /w )2 C 0 (w /w - 0.5983)2 C 0 + 0.07138 } w for o•
0.62 w<
w ~ 0.8547 w o C 0 w - w o C for wC ~0.62 w0•
If E = 0 o quency is dw h O. t en -dt . as dbr transporteà wheng~ven the value 0 50 that a constant wave
fre-all wave energy has been dissipated.
•
•
•
•
•
•
•
•
- ]3-•
References•
]. Herbers, T.H.C., Booij, N. and Holthuijsen, L.H., ]984, The design of a numerical shallow water wave hindeast model, Delft University of Technology, Report no. 3-84.
•
2. Kruseman, P., ]976, Twee praktische methoden voor het maken van ver-wachtingen voor golfkomponenten met perioden tussen ]0 en 25 seconden nabij Hoek van Holland, De Bilt, KNMI, WR 76-].
3. Sanders,J.W., ]976, A growth-stage sealing model for the wind-driven sea, DeutscheHydrographischeZeitschrift, Band 29, heft 4, pp. 136-16].
•
•
•
•
•
•
•
•
•
•
Appendix r User's manual of HrSWA•
•
•
•
•
•
•
•
•
•
•
HIS~A user manual page 1
•
•
•
•
by N. BooijI.H.C. HerbersL.H. Holthuijsen
•
Contents
•
1. Physical Background [ not yet available]2. Facilities of the program [not yet available]
3. Run instructions [ to be provided by O.l.V.]
•
4. Command descriptions [ partly available]5. Error messages [not yet available]
6. A sample problem [not yet available]
•
•
•
•
•
HIS~A user manual page 2
•
The followi~g commands are available to users of HISWA:
•
GeneralPROJECT TEST
•
RUNSTOPcommands
Title of the problem to be computed
requests the output of intermediate
testing purposes. starts a computation. end of user's input.
results for
•
Commands for model description
SET sets values of certain general parameters
GRID position and size of computational grid
BOTIOM position and size of bottom grid
READ read bottom and/or current velocity values
INC defines incident wave field
BOUND boundary conditions at lateral sides of
computational grid
WIND wind influence on wave field
BREAK breaking of waves
FRIC bottom friction
Dlf diffraction
the
•
Output commands
Net yet defined; existing version will always generate the
same type of output, irrespective of user commands.
•
•
r---,
I I I PROJECT 'TITLEl' I I I I I I 'TITLE2' ( I I I 'TITLE3' I I I L ~•
A description of the run is given in the'TIILE2' and 'TITLE3'. Each of these is
long. The lines will be reproduced in
program.
3 lines 'TITLEl', max. 56 characters
the output by the
•
r---,
I
RUN I
I
L .I
•
Starts a computation. After this one can modify certaindata, for instance the direction or frequency of the
incident wave, and start a second computation by giving RUN
again.
•
HISWA user manual page 3
•
r---,
I I
I STOP- I
I I
L .J
•
Marks the end of the users input~•
r---,
I I
I SET (LEVEL] [GRAV] I
I I
L- --.J
•
Assigns values to the (constant) water level in the region
(in m). and to the gravitational acceleration (in m/s2).
The depths used in the wave propagation computation are
equal to the difference between tha water level detined by
this command. and the bottom level read by the command READ. Init: LEVEL=O •• GRAV=9.81
r---,
•
GRID [XLEN] [YLEN ] (SECTOR] [MX] [MY] [MD] &I FIXED [XUC] [YUC] [ALUC] I
I 1
<
>
I ROl [XUR] [YUR] [XCR] [YCR ] I
( I
•
L -_~
•
Defines the computational grid. Ihe X-axis of the grid is
the computational direction, the Y-axis is normal to this
direction. The X-axis should be roughly in the mean wave
propagation direction. The orientation of the grid can
either be FIXED, or ROIATING. In the latter case the X-axis
is made coincident with the average direction of the
incoming waves; this is allowed only it •••
•
[MY]
[MD]
Meaning of the parameters:
[XLEN] length of the grid in X-direction (in m).
[YLEN] length of the grid in Y-direction (in m).
[SECTOR] directional interval of propagation directions for
which the wave energy density will be determined (in
degrees)• Ihis sector must be smaller than 180
degrees. It is assumed that the X-direction is in
the middle of the sector.
number of steps in X-direction. In view of
numerical stability the program will check whether
[XLEN ]/[MX] is smalIer than [YLENV[ MY] divided by
tan(O.5*[SECTOR]). If not, an error message is
printed.
number of steps in Y-direction.
number of subdivisions of the directional interval,
so [SECIOR]/[MD] is the spectral discretization.
•
[MX]
•
•
•
•
HIS~A user manual page 4
[XUC] [YUC] [ I\LUC ]
•
•
[XUR) [YUR] [XCR] [YCR]Fixed grid: position of computational grid in terms
m) •
the origin of the of user coordinates (in
"
Fixed grid: orientation of with respect to the.X-axis system (angle in counterclockwise) •
Rotating grid: position of the user coordinates (in m).
the computational grid of the user coordinate degrees, measured rotation point in
"
Rotating grid: position of the rotation point in computational grid coordinates (in m).
"
r---
·
---,
I I
I BOTTOM [XUB] [YUB] {ALUB] (MX] [MY] & I
I I I I I [DX] [DY] [DEPX ] I I I L ~
•
•
Defines the position and size of current field is present the current same points as the bottom.the bottom grid~ If a has to be given at the Parameters:
[XUB] position of the origin of the bottom grid In user
coordinates (in m).
•
[YUB] [ALUB]•
(MX][HY] [OX] [DY1
[DEPX]•
"
orientation of the bottom grid with respect to the X-axis of the user coordinate system (angle in degrees, measured counterclockwise).
number of meshes in X-direction of the bottom grid. number of meshes in Y-direction of the bottom grid. mesh size in the bottom grid (in m).
mesh size in the bottom grid (in m); default: equal
to [DX].
depth outside the bottom grid (in computational region contains points bottom grid, in these points the depth equal to [DEPX]. Default: [DEPX]=1000.
m); if the outside the will be taken
r---,
I I I I BOTTOM I f I I I II READ
<
I -) SEP I )&
II I I I I I I I CURRENT
<
)
I I I I I COMB I II
(I
I
I I I II [FAC] [IDFM] (IDLA] 'FORMAT' (PRINT) I
I ---- I
L _;t
•
•
•
HISUA user manual page 5
•
Uith this command the bottom configuration and/or thecurrent field is read by the program.
•
The bottom levels or the current velocity components are
read from the file according to lay-out index [lOLA], and
multiplied by [FAC]. Oefaults: [FAC]=l. and [rOLA]=l.
[IOLA]=l means that the numbers appear line by line, with
each line starting on a new input line. (IOLA]=2 means the
same order, however, a new line must not necessarily appear
on a new input line. (IOLA]=3 means that the numbers áppear
column by column, with a new column starting on a new input
line, =4 means the same order, however, a new column must
not always start on a new input line.
•
•
The format in which the numbers are written must also be
given; Free Format ([IOFM]=O) is default. The format number
(IOFH] is interpreted as follows: 5: format (20F5), 6:
(l2F6), 8: (IOF8), 2: FORMAT must be provided by the user,
such as: '(lOX,12F5)'.
Regarding the input on the file the rules tor ordinary
Fortran input apply; i.e. continuation marks are not
allowed. Simply continue on next line if a line is not long
enough to hold all the numbers.
•
In reading with free format, if a value z(i,j)is left unchanged. The initial values of the
With free format empty fields, repetition
closure of a line by a slash, can be used.
is absent, it bottom are O.
factors, and
•
Remark: HISUA uses the level of the bottom with respect to acommon datum. Iherefore a positive depth will usually meana negative bottom level. If necessary the sign of the value
can be corrected by assigning a negative value to the number
(FAC
i-•
One canfollowing command:obtain a contour plot of the bottom levelsOUTPUT ISO TITLE='BOTTOM CONTOURS, INTERVAL 5 H'
by the
•
r---I I I I INC I I I I I I L _
---
,
•
I PARAH I<
I I I I I(HSIG] [PER] [OIR] [OSPR] lI
I I I I I I I I
---~
FILE [NF] SPECTR I (<
I<
[ E ] [FREQ] )•
This command definesin not compulsory; if the commandthe incident wave field.is absent itIhis commandis assumed that there are no incident waves.Ihe incident waves can be given either in PARAMETRIC or in
SPECIRAL form.
•
HISWA user manual page 6
•
In the PARAMETRIC case the user prescribes: [HSIG], significat wave height, [PER], the neam wave period, [DIR], the mean direction of the waves, and [DSPR], which is a parameter for the directional spread, defined as follows. The program assumes that the incident waves are distributed over the directions according to the function{Cos(Iheta»**[DSPR]. The default is: [DSPR]=2.
•
•
In the SPECIRAL case the spectrum is prescribed either by giving a set of values for energy density and mean period, in which case the same spectrum is assumed for the whole boundary x=O, or by reading the values from a file, in which case the spectrum must be given for each point of this boundary.
.
'
The numbers [E] and [FREQ] must be given for each spectral direction, starting with the lowest value of Iheta, so [MD]+l occurrences of this pair must be present. Ihe command GRID must be given prior to INC in this case. [E] is the spectral energy density tor that direction, and [FREQ] is the mean frequency for the same direction.
•
With the option FILE [NF] the spectral values are read from a file; the file reference number.[NF] must be given. For each point on the boundary x=O, starting with the point y=O the whole spectrum must be present on the file.
•
r---·---~
I I I I LEFI I I INC I I I I ( I I BOUNDARY<
)
<
I I I RIGHI I f FILE [NF] I I I I I I L - --4•
Ihis command determines the boundary condition at the lateral side of the computational region. If this command is absent it is assumed that no waves enter the region from outside. The command has no effect on the waves leaving the computational region; these are always absorbed by the boundary.•
Ihe RIGHT boundary is the boundary y=O, the LEFT boundary is y=[MY J*[DY]. left boundary+---•
computational direction -)•
+---right boundaryINC: the wave spectrum at the boundary incoming waves) is the same as that of
(at least the the wave field
•
HISUA user manual page 1
•
entering through x=O, which is defined by the command INC. Applicable only in the case of INC PARAMETRIC~•
FILE [NF]: the spectra are read from a file. For each point of the boundary the whole spectrum must be present on the file, 50 also values for wave directions pointing out of ·the
region must be available on the file, although these do not have any effect on the computation•
•
.
r---,
I I
I UIND [VEL] [DIR] [A] [B] [Cl [0] I
I I
L ~
•
A souree term due to wind influence is added to the energy balance eguation. Only the values for [VEL], the wind velocity (at an altitude of 10 m, unit: mIs) and [DIR] the direction of the wind with respect to the user coordinate system (in degrees). are required. For the other parameters reasonable values are assumed by the program. These values are empirical. Ref: •••
r---
---•
I1 BREAKING [GAMS] [GAMD] [CB] [ALFA] (,
J
I ( FREQ [A] [B] )
I
---_._----~
•
Upon this command a source term For each coefficient a default program; usually it suffices to BRE FREQ.due to breaking is added. value is assumed by the give the command: BRE or:
•
If the keyword FREQ is not present. the the breaking process does not influence waves. If FREQ is given, the change of to the change of wave energy, by means Ref: •••program assumes that the freguency of the freguency is related of a power-formula.
•
For meaning of parameters see CREDI2 manual.r---,
I I
I FRICTION [fll] [FS ] (FREQ [A] [B]) I
t I
•
L-- ~
Upon this command a source term due to bottom friction is added. For each coefficient a default value is assumed by the program; usually it suffices to give the commande FRIC or: FRrC FREQ. Ref: •••
•
For meaning of parameters see CREDIZ manual.•
HISUA user manual page 8
•
r--- -
1 I I DIFFRACTION I L _---~
•
Upon this command diffraction termsIt is stressed that diffraction
roughly.
are added to the model.
is approximated rather
Command not yet formulated.
•
•
•
•
•
•
•
•
•
•
•
Appendix 11 - Subroutines of COMPU.
•
•
•
•
•
•
•
•
•
•
•
Subroutines of COMPU
x = available and/or performed
•
•
design/documentation I
name souree text tests remarks
OPENF X X X VERSIE X X INPOOL X X X REQDA X X WRCOM X X X WRCOMX X X X WRPOOL X X X ADPOOL X X X WRDUMP X X X HEAD X STARTB X X X WAVPA X X X BOCUR X X X INPDC X X X
ITWN X X X modified subroutine
VWPRO X X X NUMSC X X X PREDT X X X DISPA X FRABRE X TERMD X X X TRSY X X X TRST X X X DIFT X SWIND X SBOT X SSURF X SDBR X new subroutine SUMDE X WRIRE X STRACE X X X MSGERR X X X COPYCH X X X
•
•
•
•
•
•
•
•
•
•
•
Appendix 111 - System ,documentation of HISWA•
•
•
•
•
•
•
•
•
•
•
•
l'AtiLEOF CON1'I::NTSpage
1. Introduction 1
1.1 General characteristics ot t ne model 1
•
1.2 Computer programs 21.3 Documentation 3
2. structure of the program COM1'U !>
3. Uescription of subroutines 9
•
3.1 Sub routine U1't;Nl'· 93.2 Subroutine VEHSIE 9
3.3 subrout.ine INPOOL 9
3.4 Subroutine HEUDA 9
3.5 Subrout.ines UHCUM and IoIHCOMX 9
3.6 Subroutine WH1'OOL 9
•
3.7 sub routine AD1'UUL 93.8 Subroutine YHOUM1' 1U 3.9 Subroutine HI::AO 10 3.10 Subroutine SIARTti 10 3.11 Subroutine WAVPA 11 3.12 Subroutine BOCUH 12
•
3.13 Subroutine lN1'DC 13 3.14 Subroutine IIWN 1'1 3.15 Subroutine VWPRO 17 3.16 Subrout.ine NUMSC 18 3.17 Subroutine PHEDl' 2U 3.18 Subroutine DISPA 21•
3.19 Subroutine f'RABRI:: 23 .3.20 Subroutine IëHHD 2'1 3.21 Subrout.ine THSY 25 3.22 Subroutine IRSI 26 3.23 Subroutineuir
r
27 3.24 Subroutine SWINU 28•
3.25 Subroutine SBOT 3U 3.26 Subroutine SSUHt' 3U 3.27 Subroutine S08H 31 3.28 Subroutine SUMDt: 33 3.29 Subroutine WHIHE 3'1 3.30 Subroutine STHACt: 3'1 3.31 Subroutine MSliEHH 35•
3.32 Subrout.ine C01'YCH 3!> 4. Storage of data 3bij.1 Dynamic data pool 3b
4.2 Common bloc~s 31
ij.3 Files 40
•
Heferences 41•
•
1
•
1. INrHODUCTION•
l·! ~~~l
~haracteristi~Q1 lhg
IDQ~In this document the des~gn of a numerical shallow water waves hindcast model named HISWA is presented. Ihis model is expected to provide realistic estimates of the wave conditions in the Oosterschelde. lt is a directionally decoupled parametrie model containing bot tom refraction, wave growth, dissipation aue to wave breaking (surf zone and deep water) and bottom friction as weIl as a simple representation of diffraction effects. further the effects of currents on refraction, wind generation and bottom dissipation is included •
for the mathematical formulation of tnis mOdel reterence is made to Holthuijsen and Booij (1983).
Two balance equations in the parameters AO (frequency integrated wave action) and WO (mean wave action frequency), containing gradients in three dimensions x, y and & (wave direction) , are solved :
.
'•
t_(CXO.AO)+~(CïO.AO)+!_(CeO.AO+Cdif{CYO.~AO alXen
ae ex•
-CXO.~AO»=SO-&Q.g_WO (1) ~'i dowo
dl' ~_(CXO.WO.AO)+~(C'iO.WO.AO)+~(CSO.WO.AO alX ~ï éVS•
+Cdif(C'iO.t_(WO.AO)-CXO.!_{WO.AO»)=~Q a)X éVYdo
( 2) in eq. 1. and 2 :•
CXO, cro, cao are the components in X, 'i resp. edirection ot wave action transport velocity
do
is the relative frequency•
50 is the souree term including wind generation, bottom friction and wave breaKing (surf zone and deep water)
Cdif is a diffraction coefficient
(expressions for tnese terms are given in chapter 3)
• A numerical grid is defined in three dimensions x, y and
e
(fig. 1). The direction of wave propagatione
is defined as the angle bet ween the wave number vector and the positive x-axis, measured counter-clocKwise.•
•
2•
y kIl ---I • I..
I • I ••
I • I • I • I • I • ••
I • I • I • l • I __ 1 ___ - -______1-x
o
t n*DX ILX•
fig. 1 t.ne computational region•
•
•
•
•
•
•
rhe computation progresses in the positive x-direction and propagation of wave energy is limited to a sector defined by ea and eb around the direction of the x-axis. The computations are carried out line by line with an explicit predictor-corrector scheme. Ihe number of corrector steps is tree but two steps are sufficient to obtain a stable scheme. Lines are defined parallel to the y and e axis.
Besiàe the computational grid described above two other grids are used in the model HISWA :
- a problem grid in which the user defines his problem (in x-y plane)
- a bottom grid containing the bottom topography and current field (in x-y plane)
l.l
ComQ~l~~ QLogrêID~rne model HISWA consists of three computer programs : PREP input preparation and control part
COMPU computational part DUIP output of results
Here the design of the computational part CDMPU will be considered. rhe programs PH~P and DUIP are planned to be adjusted versions of the programs PREP and UIIV of the refraction/diffraction model CREDIZ of Rijkswaterstaat.
INPUT :
In the program COMPU 1nstructions, definitions of grids, coefficients etc. (formulated in PREP) and arrays containing bottom topography and current fieIÓ are read from a file. PtWCt::OURE:
•
3
• mean action frequency WO are determineo at the boundary x=O.
~urther depths, currents, wave numbers and wave propagation
velocity components are computeo at this boundary. for every
new line (x=n*ox) AO and WO are obtained through the
application of the numerical scheme described in section
1.1.
• OUIPUI :
for every line the following results are written to a file
that will be reaO by the program UUTP :
•
wave actLon, -frequency, rela tive frequency, group
velocity, wave number and components of wave action
transport velocity (in every grid point in the y-8
plane)
leakage of energy through the
e
boundary (for everyvalue of y)
•
dissipation of energy oue to bottom friction, surt- anddeep water breaking (for every value ot y)
the fraction of breaking waves (tor every value of y)
•
1•3A description ofQQ.~!!!.~nll.liQD. the computer program CUMPU is given inchapters 2 through q.
In chapter 2 the structure of the program COMPU 1s
explained. Ihe relations between the subroutines in COMPU
are shown in block diagrams. An example of a block diagram
is given below.
•
•
A 1 1---t l I 1--- D B I 1---c
fig. 2 block diagram
•
•
In the program or subrutine A, subroutines Band D are
called. further subroutine C is calleo in subroutine B.
Ihe structure ot a program or subroutine is presented in
Nassi-SchneiOermann diagrams. for convenience the
conventional construction in the left part of fig. 3 is
replaced by the one on the right.
•
•
•
•
•
•
•
•
•
•
•
•
•
I if ••••••••• I I T f I1---:---1
I •••••• I ••••••• I I t I I I I I I I t if ••• thenI
:---1
I I ••••• I(---:---1
I else II
:---1
1
I •••••
I
fig. 3 representation of a conditional statement
Descriptions of the various subroutines are given in chapter
3. !he sequence in which the subroutines are discussed
corresponds with the place in the program COMPU at which the
subroutine is called for the first time. Parameter lists of
the subroutines are described in which input- and output
parameters are denoted by (1) resp. (0).
In chapter ij the storage of variables and arrays in common
bloc~s and files is described. A flexible handling of
computer storage.necessary for the considerable number of
arrays in COMPU is obtained through the application of a
Dynamic Data Pool.
In this report reference is made to tne system documentation
ot CREDIZ tor a detailed description. in so tar as
subroutines and other facilities of the system CREDIZ are
•
5
• 2. SIRUCrUHE Of THE PROGRAM COMPU
The computational part COMPU of the moael HIS~A is a FORTRAN
program eonsisting of a main program and varieus
subroutines. Input is read from a file named INSTH and
output is written te a file namea HEKHES (seet10n 4.3).
• fig. 1 shows a aiagram ot the main program.
---
-
--
-
---:
•
---I
~RLL UPENf epen all neeessary files 1read aimension of pool and testparameter I
trom file INSTR I
---I
CALL lNPOUL initialize the pool and fill 1
it with empty arrays I
---I
CALL WHCUM reaa eommon bloeks and pool I
arrays from file INSTA
•
---~---•
CALL WRCOMX write eommon bleeks and pool
I arrays to file REKHES
1---I CALL ADPOUL enlarge the dimensions of pool arrays
J---I CRLL HEAU write a title
1---I CALL STARTS eompute wave parameters on line
-1 x=O and write the results to file REKRES
.1---
Itor every line do I
I
:---1
1 1 CALL NUMSC eompute wave parameters on this I
I l line and write the results to file HEKHES I
•
•
.
.
----_
..
.
---_
.
_---_.
.
fig. 1 diagram of the main program
•
•
A bioCK diagram showing the relations between the various
subroutines of COMPU is given in fig. 2. Separate diagrams
for the subroutines ~RCüM, ~RCOMX, WAVPA and IEHMUE are
ineluded in fig_ 3 through 6. The subroutines STRACE, COPYCH and MSGEHH are eallea in various parts of the program.
•
•
•
COMPU I 1---OPt:Nf I t I r---VERSIE I I I---INPOOL I I Ir
---HEQUAr
I---URCUM I ( I 1---I I---UHCOMX I 1 I (---l I---P.DPOOL I , I l---H1:;QUA I I---HEAO I I---SIARIB Ir
I , I l---WAVPAt
t I 1r
[
,
---I ( I I---WHIRE I I---NUMSCr
(---PH1:;OI t ( I---UISPA I ,;r
~---fRAfjHE t l---IEHMD I' ~ I 1---Ir
---SUMOl::r
t
r---WAVPAr
,
l t---I 1---WtUHI::•
•
•
•
•
•
•
•
bopen all necessary files
a message is given at the moment the model HlSWA is generated
initialization of the dynamic data pool
expansion Of the dynamic data pool
a major part of the common bloCKS is written to and read from a file
ditto except for the common block UITVUA
reduction or expansion of a pool array
expansion of the dynamic data pool
a title is written above the output
computation of wave parameters at the boundary x=O
computation of wave parameters at a line in the computational grid
results of computations are written to a file
computation of wave parameters at a new line
wave act ion and frequency are computed throug linear extrapolation from the former lines
computation of parameters in source terms
the fraction of breaking waves is computed
the terms of the two balance equations are computed in a grid point
evaluation of balance equations yields the wave act ion and -frequency in a grid point
computation of wave parameters at a line in the computational grid
•
fig. 2 relations between subroutines in CUMPUresults of computations are'written to a file
•
7
•
WHCOM
t
I---WH~OOL a pool array is read from or written
I'
to a fileI
r
f l---AUPOUL reduction or expansion of a pool array
l I
r
I---REQDA expansion of the dynamic data poolr
l---WHUUMP an array is printed
•
fig. 3 relations between subroutines in WHCOM
•
WRCOMXII---WRPOOL a pool array is read from or wr itten
r
t
to a filef ~
l I---AOPOOL reduction or expansion of a pool array
t I
( I---HEQDA expansion of the dynamic data pool
f
l---WHUUMP an array is printed
•
fig. q relations oetween subroutines in WHCOMX
•
WAVPA f l---BOCUH , 1 fr
t l---INPDC l,
t---ITWN r: I I---VWPKOdepths and currents are determined
at a line
•
depth and currenta bot tom grid pointare determined inwave number and relative frequency are
computed in a grid point
•
computedpropagationin a gridvelocitypointcomponents arefig. ~ relations between subroutines in WAVPA
•
rEHMD t 1---If<SY I (---IRSr I 1---DIl"Ir
(---SWINU f (---SBOr t 1---SSUHf l (---!iOtiHdivergence of transport in y-direction
•
divergence of transport in e-directiondiffraction terms
wind generation souree terms
bottom dissipation souree terms
surf breaKing sou ree terms
•
deep water breaKing souree termsfig. 6 relations between subroutines in TERMOE
•
•
•
•
9 3. OESCRIPTION Of SU~ROUTINESIn this chapter the subroutines of the program COMPU are
described. A number of these subroutines is copied trom the
model CREOIZ with a few adjustments in the source text. Only
a short description is given of the function of these
subroutines and reference is made to the documentation of
CtiJ:;OIZ.
1.1
~y'QLQutin~~NfIn this subroutine all necessary files are opened in order
to reserve input/output buffers. This action is ta~en in
connection to repeated calls of the standard routine REQDA. ibis subroutine is copied from CRt;DIZ.
1.~
~y'brQ~tin~~~~A message is printed at the moment (time and date) the model
HISWA is generated.
.. This subroutine is copied trom CREDIZ.
1·1
~llQu!:i!l~~QQb
Tne dynamic data pool is initialized by this subroutine.The
dimension ot the pool is determined from the common variabie
N~OOL (NPOOL
*
1U24). A number of empty arrays is initiated(SO in he program COMPU) •
lNPOOL is copied trom CRt;DIZwith minor adjustments.
•
•
•
•
•
•
•
1.~
SUQLQylin~ Rt;QQaIhe standard routine Rt;QUA,copiedtrom CREUIZ, is used for
expansion of the dynamic data pool.
1.2
~~QLQ~1in~ URCOtl2nQ WRCOMXA major part of the eommon blocks is written to and read
from a file by the subroutines WRCOM and WRCOMX. This is
necessary for the communication between the programs PREP,
COMPU and OUIP. the difference between WHCOM and WRCOMX is
the fact that WRCOMX doesn't read or write the eommon bloek
UrIVUA, containing instructions and information for the
program OUTP.
WHCDM and WRCOMX are copied trom CtiEOIZ with minor
adjustments.
1-2
~~QLQu1in~~BPQQb
rhe subroutine WtiPOOL reads or writes a pool array
(unformatted) trom resp. to a file.
WHPOOL is copied from CHEOIZ with minor adjustments_
1.1
~~bqH!11n~ ADPQQbIhis routine is calleO
program for shrinking
Oynamic data pool.
Subroutine ADPOOL is adjustments_ by WRCOM,WRCOMX or expansion of and in the an array in main the
•
10
• 1·!! ~~!ll:Q~!.ingWHD
!ll:1f
Ihe contents of an array is printed by the subroutine
WtWUMP.
WHUJMP is copied from CHEOlZ with minor adjustments.
•
•
•
•
•
•
•
•
•
•
·function:A title is printed above the output of the program COMPU.
Ihe callof this subroutine is :
CALL HEAU(strl,str2,str3) Parameters :
strl,strl,str3 are character strings containing a title
function :
In the subroutine START~ the wave conditions at the boundary x=O are determined.
Method :
Ihe boundary condition is represented here by the
directional integrated wave action A(y) and -frequency W(y).
Ihe distribution of A(y) over the directions is obtained by
a cos**m directional distribution. The parameter m is read
trom file INSTto{.
m
AO(l{,tt) =cvcos (e-v).A(Y) for ,~-vl<90 deg.
(1)
=0 for le-vl290 deg.
wo(l{,e) =W(0 (2)
v is the wind direction defined as the angle between the
wind velocity vector and the positive x-axis, measured
counter-clockwise.
rhe normalization coefticient c, consisting
functions is determined from approximations
Abramowitz and Stegun (1965).
Oepths and currents are determined
Wave number, relative frequency.
components of propagation velocities
grid point in the l{-e plane.
Next these parameters are written to the file REKHES.
Other options for the representation of boundary wave
conditions can be implemented. .
of
gamma-g1ven in
at the boundary x=o.
group velocity and
•
11
•
Structure :*SIARIB*
:---_.
•
I determine wind direction relative to theI computational grid1---I compute the direct10nal distribution of waves
1---I 1f option 1 is chosen then
I
:---I for every y do
1
:---I for every S do
,
:---I I compute wave action AO and
I I wave frequency WO
•
•
1---:---:---:---1
•
I else if option 2 is chosen then I
:---~---I
1 I •••••••••••••••••••• I
1---:---(
I CALL WAVPA determine wave parameters I
f at boundary x=O I
1---1
1 CALL WRIR~ write wave parameters to file REKRES I
._---_.
.
.
•
l'hecallof this subroutine is :CALL STARTS
•
function :In this subroutine wave numbers KO, relative frequency 00,
group velocity
e~o
and propagation velocity componentsexo,
ero
and CSO are computed at a line IX in the computationalgrid.
•
Method :In order to evaluate these parameters first the depths 0 and
current velocity components UX,
ur
at line IX aredete rm Lned ,
•
•
•
•
12
• Structure :
:---~---:
•
I if predictor is passed or line is boundary I
I x=o then I
I
:---1
( if predictor is passed then I
I
:---1
I I move arrays containing depths and currents 1
I 1 l at line IX to arrays w1th oid values I
I
1---1
I I CALL BOCUH cornpute depths and currents I
I I at line IX 1
I---~---I
J for every y doI
I :---( 1 I tor everye
doI
I
:---1
I I I CALL lIWN calculate wave number and 1
I f
t
relative frequency 1I
1
t---I
I I I CALL VWPHO calculate group velocity and I I t I components of propagation velocity f
•
•
•
._---_._---_.
. .
.
.
Ihe callof this subroutine is : CALL WAVPA(IX)
•
parameter :IX (I) line in computational grid at which wave parameters are computed (X=(IX-1)DX)
J.!1
~ubrQgt1~ §QkQft•
function :In this subroutine depths D(Y) and current velocity components UX(Y), UY(Y) are determined at a line in the computational grid.
•
Method :
ror every grid point, bottom grid coordinates are computed and depth and current velocity components are determined through bilinear interpolation.
•
•
•
13 • Structure : *BOCUH*---:
tor every y do:---f
•
t determine bottom grid coordinates of point I
1---1
I CALL lNPUC determine depth and current f
t in point I
(---1
I
if current is on then1
(
:---1
, t .determine current relati ve to
I t computational grid
•
--_._-_._---_.
.
.
.
•
Ine c311 of this subroutine isCALL BOCUR (.IX)
parameter :
IX (I) line in computational grid at which depths
and currents are computed (X=(IX-l) .DX)
•
function :
Oepth and current velocity components are computed in a
pOint given in bot tom grid coordinates (IB,JB).
•
Method :A bilinear interpolation is carried out
points in the bottom grid. If point
outside the bot tom grid then a constant
is assume d ,
with the surrounding
(lB,JB) is Iocated
depth and no current
•
•
•
•
•
• Structure :•
'
.
•
*INPDC* .:---_.
1 if point is located in bottom grid then I I I 1 1 I I I 1 I
:---·1 compute c eptn
1---I if deptn is positive and current is on then
I
:---I I compute current components
r---:---, else I :---1 t current is o ,
,
_---_._---
.
.
--else:---r depth is constant value outside bottom grid
l and current is O.
:---:---•
Ine callofCALL lNPDCthis subroutine is :function :
In IIWN tne wave number KO and relative frequency do in a
grid point (IX,Iï,rI) is determined.
•
Method :Ihe current component in tne direction of wave propagation
is determined.
•
u
= Ux.cos& + Uy.sin& ( 1)lf U ) 0 then the wave number KO is computed through a
Newton-Hapns~n iteration process, applied to eq. 2.
•
1/2
r ~
WO-KO(UX.cos&+UY.sine)-(g.KO.tanh(KO.D»=
0 ( 2)Ihis procedure requires an estimate of the wave number KO as
a start value. Here the value of KO on line IX-1 is used as
an estimate for KO. At the boundary x=O the following
approximation of KO is applied.
•
•
•
g.KO 1 = (3) 2 110 2 1/2•
1~
• Tne relative freguency 00 is calculated witn eg. 4.
1/2
0'0= (g.KO.tanh (KO.D»
• 11 U ( 0 f1rst the freguency Wc is determined wnicn is the highest
freguency capable to transport wave energy against the current U.
For this purpose eg. ~ is solved through a Newton-Haphson iteration
pr oces s, 1
o
•
~ - U + 01.( ---- + --- ) =0 (S) l.Kl sinh(2.Kl.U) witn 1/2•
dl
=
(g.Kl.tanh (Kl.U» Wc is given by Wc=
Kl.U +0'1
•
If WO<
Wc (eg. 2 na s acomputed with the procedure
2 has no s olu t tonj then Kl
KO resp.
oo ,
(b) (7)
solution) then KO and 00 are
described above. If WO
>
Wc (eg.and
Ol
are used as estimates of•
Acomponentssou ree termwith freguencySdbr is>
Wcintroduced(section to3.21).dissipate wavefor points outside the bot tom grid KO is computed with eg. 8.
2
1<0=1010/g (8)
•
1t a negative depth is encountered then KU and ÓO are given thevalues -1.0 resp.
o
.
o.
•
•
•
•
•
16
• Structure :
:---:
•
( if point is located outside bottom grid then
1
1
:---1
1 t compute wave number and relative frequency 1
1---:---1
•
I else if depth is negative then I
1
:---1
1 1 wave number is -l ••relative frequency is
o.
I1
---:---1
1 else 1
I
:---1
I 'give an estimate tor the wave number Ke 1
1
,---1
1 I compute current component in direction of 1
I I wave propagation U
r
1
1---1
1
'i
t U<
0 th en 11
,
:---1
I I 1 compute estimate tor wave number Kl 1
1
I
t---I
I I 1 compute tunction G I
I t
1---1
I 1 1 for i= 1 to ~O while G
>
accuracy dO 1I I I. :
---1
I
I
1
1
compute deri va ti ve of G (1
,
I
1---1
1 I I I compute wave number Kl I
1
1
1
1---1
I I l 1compute tunct Lon G I
1 I 1---: ---(
I I t compute critical frequency I
1 t ---: ---.---(
I I else I
f
I
:---1
I I I cri tical trequency is -1. 1
1
1---:---1
( I it critical trequency
<
0 or>
wave frequency II
,
:---1
I I t compu te function f' I
1
I
&---1
I 1 t tor i
=
1 to 50 while F>
accuracy do 11
I
I
;---1
1
r
I I compute ce rrvatrve of f I1
I
I
1---1
1 t I I compute wave nunoe r Ke I
I
I
I
I---~---l
1 Ir
I compute function F II
I
·
l---: ---1
1 I I else 11
1
t:---~---I
1 I f· I wave number is Kl 1•
•
•
•
•
•
._-_._-_._-_.---_.
.
.
.
.
.
•
•
11
• Ihe call o~ this subroutine is : CALL II~N(IX,lï.lT)
•
•
•
Parameters : IX (I) I ï (I) II (I) ) )coordinates of point in computational grid
Function :
In subroutine VWPRO the components of wave action C~O are determined in computational grid.
group velocity CGO transport velocity
exo.
a point (lX.IY.II)
and the
cr o
and in the Method :Ihe relations for the parameters mentioned above used in this model are :
1 0
.,
CGO = 00 (----
+---
) (1)2.KO sinh (2.KO.0)
exo
= CGO.cos& +ux
(2)cro
= CGO .sin& +ur
(3)•
0'0 il)0 ëJ)U «vUXeux
•
•
•
•
•
C~O= - ---(-sin~--+cos9--)-cose(-sine---+cos&---) sinh(2.KO.D)ax
iYax
in
Q)ur
mur
-sin9(-sin&---+cos&---) mX éi)Y (U)Ihe terms containing current velocity components UX and
ur
in eg. 2 through 4 are omitted if no current is present. The termceo
is evaluated intermediate lines IX and IX+l. Uerivatives of depth and current are determined through a central difference scheme (after the predictor step).IC negative depths are encountered then all velocity components are given the value 0 ••
•
18 • $tructure :•
•
•
•
•
•
•
•
•
---_.
.
.
it depth is negative then
:---I give velocity components and derivatives
l the value
o.
---:---else
:---l compute group ve:---locity
t---~---I if predictor step is passeo thenr
:---I
r
compute depth derivatives and currentI t deriva tives
1---·-.
.
I compute components of wave action transport f
I velocity I
l (C&O only if line ~ boundary x=O) I
--_._---_.
.
'
.
Ihe callof this subroutine is : CALL VUPHO{IX,IY,II) l'arameters
.
.
IX (1.) ) lY (I) ) lT (I)coordinates of point in computational grid E'unction:
In the subroutine NUMSC wave parameters are computed at a new line IX+1 in the computational grid.
Method :
Ihe following numerical scheme is applied :
~stimates tor the wave action AO and -frequency WO at line lX+l are obtained through a linear extrapolation from the lines IX-l and IX (predictor step). Uith these estimates the other wave parameters at line IX.l can be determined. Linear interpolation between the lines IX and IX+l yields the wave parameters at line IX+1/2, necessary for the corrector step. Ihe corrector step (which can be repeated several times) consists of an explicit differential scheme,applied to the two balance equations described in section 1.1.
Ihe amount of energy lost through dissipation (fU) and lea~age of energy through the boundaries &a and &b (fL) is kept.
•
19
• Hesults of the computations are written to the file REKRES.
Structure :
*NUMSC*
•
._---_.
·
.
•
I CALL PR~DT predictor estimates for wave action t
1
and -frequency on line IX+1I
1---1
1 move contents of arrays with new values of wave 1
1 number.relative and critical frequency and propaga-I
I
tion velocity components to arrays with Oid values1
1---1
I CALL WAVPA compute wave parameters at line IX+l I
1---1
t determine depths and currents 1
1
intermediate lines IX and IX+11
1---1
1 for every corrector step do 1
I
:---1
.
1
I
determine wave parameters intermediate linesI
1 1 . IX and IX+l 1
I
1---1
I I CALL OISPA compute parameters in dissipation II I terms I
I
1---1
1 t for every
r
do I1
1
:---(
1
I
1
if last corrector step is in ptogress tnenI
I
I
t:---1
1 I t ,initialize Iea~age and dissipation
I I t 1 in point x,y I
1
I
[---:---1
I
I
1
if depth is positive then1
(
1
I
---1
I I l for every ij do I
1
I
1 :---1
f t I I CALL T~RMD compute terms of the I
1 1 t t two balance equations 1
I '1 t
---t
( 'I I CALL SUMUE determine wave action I
I I t I and frequenoy I
I 1 I ---:---(
I I t 1f last corrector step in progress thenl
I
(l
:---1
I til compute leakage and dissipation in I
I t I t point x,y I
I
1---:---:---:---1
I t CALL WAVPA compute wave parameters on line IX+1!
I~--:---r
I CALL WHIHE write results of line lX+1 to REKH~S I
•
•
•
•
•
•
·
---_
.
_---_.
·
.
•
•
•
20
• Ihe callof this subroutine is :
CAL!. NUMS((IX)
Parameter :
IX (I) wave parameters are determined at line IX+1 in the computational grid (X=IX.DX)
•
f'unction :
Estimates for the wave action AO and -frequency UO at line
IX+l (predictor step) are determined in this subroutine.
•
Method :Ihe predictor is a simple extrapolation procedure. AO and WO at line lX+1 are determined as follows :
•
AOIX+l = 2.AOIX - AOIX-l (1.)IX+l IX IX-1
WO = 2.1010 - WO (2)
lf a negative deptn is encountered then AO and WO are given the value 0 ••
•
•
•
•
•
•
•
21 • Structure : :::PREOI*._---_.
.
.
•
I if line is boundary x=o then
I
:---1
I ,wave action and -frequency on new line are
I • given the values on the oid line I
1---:---1
I else I
I
:---1
I I for every y do
I
I
:---1
I I I if depth is positive then
I
,
,
:---1
I 1 I for every
e
doI I t : --- t
I I f I move wave action and -frequency
1
I
I
,
to arrays with oId values and1 t I I compute new values
I
,
1---:---:---I l I e1se
I I I
:---, I I for every
e
dof I I :---I I I r move wave action and -frequency
I I I I to arrays with oid values and
I
t
I - 1 give new values tne valueo.
.
'
•
•
:---:---:---:---:---•
The callof this subroutine is :CALL PREOT (LX) Parameter :
IX (I) wave parameters are determined at line lX+1 in the computational grid (X=IX.OX)
•
•
function :
In this subroutine parameters at line 1X+1/2 are determined, necessary for the evaluation of the dissipation terms in the two balance equations.
•
Method :
rhe following pararueters are determined : orbitaI velocity at the bottom Ubot current velocity at the bottom Ucur wave energy density ~t
local maximum wave height Hm
the fraction of breaking waves Qb
Yor these parameters the foIlowing relations are used :
•
•
22•
2 2 9 .KO .60.AO 1/2 Ubot :: (ue.r ---
)
& l. WO •cosn(K 0 .0) (1)•
2 2 1/2Ucur :: (UX + Uï )
Et
=
De.[O'O.AO&
(2)
(3)
•
_-1HlII :: 0.88.KO .tanh(t.KO.O/0.88) , KO
=
1 --·LKONe e
(4)•
l-Qb Et=
-8.---2 (evaluated in ~HA~R~) (5) lnQb Hm Structure :•
*UISPA*._---.
.if bottom dissipation is on then
:---•
t for every y do
I
:---I I cornpute orbital velocity at the bottolll
1---:---•
t if current is on then t t I I I • tor every y do t 1 If compute current veloc1ty at the
1 bot tom
._---.
:---_._-_._-_.---
. .
.
if surf breaking is on then
._---.
•
t for every y do
I
:---t I cornpute wave energy density
I
l---t I compute local maximum wave height
(
I
l---I I I CALL fRAtlHE cornpute traction of breaking
I I e waves
•
:-
-
-:---:---The callof this subroutine is :
CALL DlSPA