The design of a numerical shallow water wave hindcast model

(1)

(2)

•

Delft University of Technology Department of Civil Engineering Group of Fluid Mechanics

Report No. 3-84

•

The design of a numerical

shallow water wave hindeast model

• T

.

H.C.

Herbers N. Booij L.H. Holthuijsen

•

(3)

•

PROJECT REPORT Delft University of Technology Department of Civil Engineering Group of Fluid Mechanica

•

Project title GEOMOR wave model

•

Project description Development of a two-dimensional

model to hindcast spectral wave parameters in an estuary with tidal flats on the basis of bot-tom topography, current- and

wind data

•

Customer Rijkswaterstaat

Deltadienst

Afdeling Kustonderzoek Van Alkemadelaan 400

2597 AT The Hague, The Netherlands

•

repreaented by

J.

van Marle

Project leader dr. ir. L.H. Holthuijsen

•

work carried out by _{T.H.C. Herbers}

dr. ir. N. Booij

dr. ir. L.H. Holthuijsen

•

Conclusion

HISWA is a directionally decoupled parametric wave hindcast model containing bottom- and current re-fraction, difre-fraction, wave growth and dissipation. The design for this hindcast model is presented with emphasis on the system docu-mentation of the computational program of HISWA.

•

Status of report Confidential, progress report

•

City/date: Delft, March 21, 1984

(4)

•

TABLE OF CONTENIS

page

1. Introduetion 1

1.1 General eharacteristies of the model 1

•

1.2 Computer programs 2

1.3 Doeumentation 3

2. Strueture of the program COMPU 5

3. Deseription of subroutines 8

•

3.1 Subroutine OPENF 8

3.2 Subroutine VERSIE 8

3.3 subroutine INPOOL 8

3.4 Subroutine REQDA 8

3.5 Subroutines ~RCOH and ~RCOHX 8

3.6 Subroutine ~RPOOL 8

•

3.7 Subroutine ADPOOL 8 3.8 Subroutine ~RDUMP 9 3.9 Subroutine HEAO 9 3.10 Subroutine STARTB 9 3.11 Subroutine IJAVPA 10 3.12 Subroutine BOCUR 11

•

3.13 Subroutine INPDC 12 3.14 Subroutine ITWN 13 3.15 Subroutine VIJPRO 14 3.16 Subroutine NUMS 16 3.17 Subroutine PREDT 18 3.18 Subroutine DISPA 20

•

3.19 Subroutine FRABRE 20 3.20 Subroutine TERMOE 21 3.21 Subroutine TRSY 22 3.22 Subroutine TRST 23 3.23 Subroutine DIFT 24 3.24 Subroutine SWIND 25

•

3.25 Subroutine SEOT 27 3.26 Subroutine SSURF 27 3.27 Subroutine SUMDE 28 3.28 Subroutine WRIRE 30 3.29 Subroutine STRACE 30 3.30 Subroutine MSGERR 30

•

3.31 Subroutine COPYCH 30 4. Storage of data 32

4.1 Dynamic data pool 32

4.2 Common b.Loc k s 33

4.3 Files 36

•

Referenees 37

•

(5)

•

1

•

1. INTRODUCTION

• 1·1 Ge~~l

fhg!2fl~ri~lic~

Qf lhg mQ~gl

In this report the design of a numerical shallow water waves

hindcast model named HIS~A is presented. This model is

expected to provide realistic estimates of the wave

conditions in the Oosterschelde. It is a directionally

decoupled paraffietric model containing bot tom refraction,

wave growth, dissipation due to wave breaking and bottom

friction as weIl as a simple representation of diffraction

effects. Further the effects of ~urrents on refraction, wind

generation and bottom dissipation is included.

For the mathematical formulation of this model reference is

made to Holthuijsen and Booij (1983).

Two balance eguations in the parameters AO (frequency

integrated wave action) and ~O (mean wave action freguency),

containing gradients in three dimensions x, y and

&

(wave

direction) , are solved :

•

!_(CXO.AO)~~_(eYO.AO)+~_(eeO.Ao~edif(eYO.~_AO

~X

~Y

~e

mX

•

-eXO.~_AO»=SO-aQ.g_~OWY ~O ~O dI (1) !_(eXO.WO.AO) +~{CYO.~O.AO)+~(eeO.WO.AO ~X

~y

we

•

+Cdif(CYO.~(WO.AO)-eXO.~_(WO.AO»)=~Q~SOWX

wY

~O

(2) in eg. 1 and 2 :

• exo,

directionCYO,of wave

eso

areactionthe componentstransport velocityin X, Y resp.

e

ao

is the relative frequency

SO is the source term including wind generation,bottom

friction and surf breaking

•

Cdif is a diffraction coefficient

(expressions for these terms are given in chapter 3)

•

A numerical grid is defined in three dimensions x, y and

&

(fig. 1). Ihe direction of wave propagation

e

is defined as

the angle between the wave number vector and the positive

x-axis, measured counter-clockwise.

•

(6)

•

2

•

y

LXI

~---I--I • 1 1 • 1 I • 1 1 • I 1 • I I • 1 1 • 1 I • I I • I • I • I • I • I I • I __I I__

x

Ol n=::OX lLX

•

fig. 1 t~e computational region

•

The 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. The number of corrector steps is free but two steps are sufficient to cbtain a stabIe scheme. Lines are defined parallel to the y and

e

axis.

•

Beside the computationöl 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)

•

• 1·~

~!!ill.Y.!.~1: Ql:Q9Lill!l~

The model HlSWA consists of three computer programs : PREP input preparation and contrcl part

COMPU computational part DUTP output of results

Here the design of the computational part COMPU will be considered. The programs PREP and OUTP are planned to be adjusted versions of the programs PREP and UITV of the refraction/diffraction model CREDIZ of Rijkswaterstaat.

•

INPUT :

In the program COMPU instructions, definitions of grids, coefficients etc. (formulated in PREP) and arrays containing bottom topography and current field are read from a file. PROCEDURE :

First the values of the frequency integrated action AO and

.

'

,

(7)

•

3

•

mean action frequency WO are determined at the boundary x=O.

Further depths, currents, wave numbers and wave propagation

velocity components are computed at this boundary. For every

new line (x=n*dx)' Aa a nd WO are obtained through the

application of the numerical scheme described in section

1.1.

OUTPUT :

For every line the following results are written to a file

that will be read by the program OUIP :

•

wave action,-frequency, relative frequency,group

velocity, wave number and ccmpcnents of wave action

transport velocity (in every grid point in the y-e

plane)

leakage of energy through the e boundary (for every

value of y)

•

dissipation of energy due to bottom friction and surf

breaking (for every value of y)

the fraeticn of breaking waves (for every value of y)

•

A1.3description

Do~Ym~nlali2n

of the computer program COMPU is given in

ehapters 2 through 4.

In ehapter 2 the structure of the program COMPU is

explained. The relations between the subroutines in COMPU

are shown in block diagrams. An example of a block diagram

is given below.

•

Il. 1 1---I I I 1--- 0 B I 1--- C

fig. 2 bloek diagram

•

In the program or subrutine A, subroutines Band D are

called. Further subroutine C is called in subroutine B.

The structure of a program or subrcutine is presented in

Nassi-Schneidermann diagrams. For convenienee the

conventional construction in the left part of fig. 3 is

replaced by the one on the right.

•

(8)

•

ij

•

I if ••••••••• I T 'F I

1---:---

-

--1

I •••••• I ••••••• I

1

r

I

til I ( I I if ••• then I

I

:---

-

--

-

1

( I· •••• I

1---:---1

1

else

I

1 :---1

I I. •••• I

•

fig. 3 representation of a conditicnal statement

Descriptions of the various subroutines are given in chapter 3. The 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. CC).

In chapter ij the storage of variables and arrays in common

blocks and files is described. A flexible handling of

• computer storage,necessary for the considerabie number of arrays in COMPU is obtained through the application of a Dynamic Data Pool.

In this report referrence is made to the system documentation of CREDIZ for a detailed description, in so far as subroutines and other facilities of the system CREDIZ

• are implemented in the present model.

•

(9)

•

5

• 2. STRUCTURE OF THE PROGRAM COMPU

The eomputational part COMPU of the medel HIS~A is a FORTRAN program eonsisting of a main program and various subroutines. Input is read from a file named INSTR and output is written to a file named REKRES (seetion 4.3). • Fig. 1 shows a diagram of the main pregram •

.

_---_.

.

•

1 CALL OPENF open all neeessary files 1

1---I

1 read dimension of pool and testparameter 1

1 from file INSIR 1

1---I

1 CALL INPOOL initialize the pool and fill 1

1 it with empty arrays 1

1---I

1 CALL WRCOM read eommon bloeks and"pool 1

1 arrays from file INSTR 1

1---(

1 CALL WRCOMX write eommon bloeks and pool 1

1 arrays to file REKRES 1

1---(

1 CALL AOPOOl enlarge the dimensions of pool arrays 1

1---I

I CALL HEAO write a title 1

1---I

1 CALL STARTB eompute wave parameters on line 1 1 x=O and write the results to file REKRES 1

1---I

1 for every line do 1

I

;---1

I

1

CALL NUMS eompute wave parameters on this

1

1 I

line and write the results to file REKRES

1

•

• ._---_._---_.

.

fig. 1 diagram cf the main program

•

_{A bloek} _{diagram showing the} _{relaticns between} _{the various}

subroutines of COMPU is given in fig. 2. Separate diagrams for the subroutines WRCOM,WRCOMX,WAVPA and TERMOE are ineluded in fig. 3 through 6. The subroutines STRACE.COPYCH • and MSGERR are ealled in various parts of the program.

•

(10)

•

COMPU

1

I---OPENF I , I

,

,---VERSIE 1 (---INPOOL 1 , ( (---REQDA

1

(---\JRCOM I I I 1---1 I---\JRCOMX I I 1 1---I l---P.DPOOL ( (---REQDA ---HEAO ---STARTB

,

( l---\JAVPA I I [ 1---( I---WRIRE ---NUMS I I---PREOT I I I---DISPA ( ( I I---FRABRE I I---TERMDE I I I 1---I (---SUMDE

1

I I---\JAVPA I I , 1---I ,---WRIRE

•

• '

.

•

6

open all neeessary files

a messa ge is given at the moment the model HISWA is generated

initialization of the dynamie data pool expansion of the dynamie data pool

a major part of the eommon bloeks is written to and read from a file

ditto exee~t for the eommon bloek UITVDA

reduetion or expansion of a pool array expansicn cf the dynamie data pool a title,is written above the output eomputation of wave parameters at the boundary x=O

eomputation of wave parameters at a line in the eomputational grid

results of eomputations are written to a fil eomputation of wave parameters at a new line wave action and frequency are eomputed t~rou ~inear extrapolation trom the former lines eomputation ot parameters in souree term~ the fraction of breaking waves is computed the terms of the two balance equations are eomputed in a grid point

evaluation of balanee equations yields t~e wave action and -frequency in a grid point eomputation of wave parameters at a line in the computational grid

results of computations are written to a til • fig. 2 relations between subroutines in COMPU

(11)

•

7

•

WRCOM

I

I---WRPOOL a pool array is read from or writter

I to a file

I

I---ADPOOL reduction or expansion of a pool array I

I---REQDA expansion of the dynamic data pool

•

fig. 3 relations bet ween subroutines in ~RCOM

•

WRCOMXI

I---WRPOOL a pool array is read from or written

I to a file

I

I---ADPOOL reduction or expansion of a pool array I

I---REQDA expansion of the dynamic data pool

•

fig. ij relations between subroutines in WRCOMX

•

WAVPA I l---BOCUR I I I I I I---INPDC I I I---IIWN I I I---VWPRO

depths and currents are determined at a line

depth and current are determined in a bottom grid point

•

wave number and relative frequency are

computed in a grid point

propagation velocity components are computed in a grid point

• fig. 5 relations bet ween subroutines in WAVPA

•

TERMDE I I---TRSY I I---TRST I (---DIFT I I---SWIND ( I---SBOT I I---SSURF

divergence of transport in y-direction

•

divergence of transport in e-direction

diffraction terms

wind generation source terms bottorn dissipation souree terms surf breaking source terms

fig. 6 relations between subroutines in TERM DE

•

(12)

•

8

• 3. DESCRIPTION OF SUBROUTINES

In this ehapter the subroutines of the program COMPU are deseribed. A number of these subroutines is eopied from the model CREDIZ with a few adjustments in the seuree text. Only a short deseription is given of the funetion of these • subroutines and referenee is made te the doeumentation of

CREDIZ.

1.1 ~~Q~ln§

QPENF

In this subroutine all necessary files are opened in order to reserve input/output buffers. This aetion is taken in • eonneetion to repeated ealls of the standard routine REQDA.

This subroutine is eopied from CREDIZ.

~.l

Su!ll.Qy!ing'yERSI~

A message is printed at the moment (time and date) the model HIS~A is generated.

• This subroutine is eopied from CREDI2.

1.1

~~!.ing lNPOOl

The dynamie data pool is initialized by this subroutine.The dimension of the pool is determined from the eommon variabIe NPOOL (NPOOL

*

1024). A number of empty arrays.is initiated • (50 in he program COMPU) •

INPOOL is eopied from CREDIZ with miner adjustments.

1.!!_

~~broutin~ BmDA

Ihe standard routine REQDA,copied from CREDIZ, is used for expansion of the dynamie data pool.

• 1.~

Su.Q.[outi!l~~RCOt1gIl!!HE~QMX

A major part of the common bloeks is written to and read from a file by the subroutines ~RCOM and ~RCOMX. This is necessary for the eommunication between the programs PREP, COMPU and OUIP. The differenee between ~RCOM and ~RCOMX is the fact that ~RCOMX doesn't read or write the eommon bloek UITVDA, eontaining instructions and information for the program OUTP.

~RCOM and ~RCOMX are eopied from CREDIZ with minor adjustments.

•

1.2

SubrQY!ing ~RP001

The subroutine WRPOOL reads or writes a pool array (unformatted) from resp. to a file.

~RPOOL is eopied from CREDIZ with miner adjustments.

• 1.1 ~QIQYling

This routine is8DPQQ~ealled by WRCOM,WRCOMX and in the program for shrinking or expansion of an array in dynamie data poel.

Subroutine ADPOOL is adjustments.

main the eopied from CREDIZ with a few

•

(13)

•

9

• J•ft

~!2LQ1!11n~

!i

RD

u

1:1f

The contents of an array is printed by the subroutine IoJRDUMP.

~RDUMP is copied from CREDIZ with minor adjustments.

•

Function:

A title is printed above the output of the program COMPU. The callof this subroutine is :

CALL HEAD(str1,str2,str3)

•

Parameters :

str1,str2,str3 are character strings containing a title

• Function :

In the subroutine STARTB the wave conditions at the boundary x=O are determined.

•

Method :

.The boundary condition is represented here by the directional integrated wave action A{y) and -frequency ~(y). The distribution of A{y) over the directions is obtained by a cos**rndirectional distribution. Ihe parameter rn is read from file INSTR.

•

m

AO (Y

,e)

=c vcos (9-v).A(Y) =0

for 19-vl<90 deg.

tor le-vl~90 deg.

( 1)

•

~O (Y

,e)

=~(Y) (2)

•

v is the wind direction defined as the angle between the wind velocity vector and the positive x-axis, measured counter-clockwise.

The norrnalization coefficient c, consisting functions is determined trom approximations Abramowitz and Stegun (1965).

Depths and currents are determined .

Wave number, relative frequency, components of propagation velocities grid point in the Y-9 plane.

Next these parameters are written to the file REKRES.

Other options for the representation of boundary wave conditions can be implemented.

of gamma-given in

•

at the boundary x=O. group velocity and are computed in every

•

(14)

•

10

•

Struct ure e *STARTS:::

._---_.

.

•

1 determine ~ind direction relative to the1 com~utational grid I

---1

compute the directional distribution of waves I

---1

if option 1 is chosen then I

:---1

( for every y do I

I

:---1

1

tor every e do

I

:---1

I

compute ~ave act ion AO and

1

f 1 1 wave frequency

~o

1

•

• ---:---:---:---1

else if option 2 is chosen then

:---1

••••••••••••••••••••

---:---1

•

CALL ~AVPA determine ~ave parametersat boundary x=O 11

1---1

1 CALL ~RIRE write wave parameters to file REKRES I

._---_.

.

•

The callof this subroutine is : CALL STARTB

•

Function :

In this subroutine wave numbers KO, relative frequency 00, group velocity eGO and propagation velocity components CXO, CYO and

ceo

are computed at a line IX in the computational grid.

•

Method :

In order to evaluate these parameters first the depths D and current velocity components UX, UY at line IX are determined.

The subroutine IT~N requires an estimate of the wave number KO as input parameter. 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.

•

(15)

11

•

g.KO 1 = (1) 2 WO 2 1/2 (tanh(WO .D/g» Structure :

•

*WAVPA*

._---:

•

if predictor is passed or line is boundary 1 x=O then I

:---1

1

if predictor is passed then

1 I

:---1

, 1 move arrays containing depths and currents I

1

at line IX to arrays ~ith old values

1 1---1

I CALL BOCUR compute depths and currents I'

, at line IX I

---1

for every y do

1

• :---1

•

I for every

e

do I

1 :---1

I 'give an estimate for the ~ave number 1 I

1---1

I I CAlL ITYN caIcuIate wave number and I

, I relative frequency I

i

1---1

I I CALL VWPRO caIcuIate group velocity and I I I components of prcpagation velocity I

• ._-_._-_._---_.

.

The callof this subroutine is :

CALL WAVPA (IX)

•

parameter :

IX (I) line in computational grid at'which ~ave parameters are computed .(X=(IX-1)DX) 1.12 ~ub[Qyting BOCU!!

•

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 :

For every grid point, bottom grid cocrdinates are computed and depth and current velocity components are determined through bilinear interpolation.

•

(16)

12 • Structure : :,":BOCUR*

._---_.

.

t or every y do

:---1

•

I determine bottom grid coordinates of point I

1---1

I CALL INPDC determine depth and current I

I (if current is on) in point I

1---1

I

if current is on then

1 I

:---1

1 I

determine current relative to

I

1

computational grid

I

• ._---_._---_.

.

. .

.

•

The callof this subroutine is : CALL BOCUR (IX)

•

parameter ;

IX (I) line in computational grid at which depths and currents are computed (X:(IX-1).DX)

Function :

Depth and current velocity components are computed in a point given in bottom grid coordinates (IB,JB).

•

Method :

A bilinear interpolation is carried out points in the bottom grid. If point outside the bottorn grid then a consta~t is assumed ,

with the surrounding (IB,JB) is located depth and no current

•

(17)

-13 • Structure : :{:INPDC*

:---:

•

I if point is located in bottom grid then I I ---(

I compute depth I

I

---1

I if depth is positive and cur~ent is on then 1

I

:---~---I

I I compute current components I

I

---:---1

I else if depth is negative and current is on then I

I

:---1

I I current is

o.

I

1--- ---:---1

I else

J

:---1

I 1 depth is constant value outside bottom grid I

I

(---1

I I if cur rent is on then I

I

:---~---I

I I I current is

o.

I

:---:---:---I

•

The callof this subroutine is : CALL INPDC

•

Function :

In ITWN the wave number KO and relative freguency 00 in a grid point (IX,lY,II) is determined.

Method :

• KO is computed through a Newton-Raphson iteration process, applied to eg. 1.

1/2

F

=.

WO - K 0 (U X •cos

e

+ UY •5in

e) -

(g•K0 •tanh (K0 •D))

=

0 (1)

•

(the second term in eg. 1 falls off if no current is present) For points outside the bottom grid KO is computed with eg. 2.

2

KO=WO Ig (2)

•

00 is calculated with eg. 3.

1/2

ero::

{g.K 0 • ta nh (K 0 • 0) ) ( 3)

•

(18)

14

• If a negative depth is encountered then KO and 00 are given the values -1.0 resp.

o.o.

Structure

•

*ITWN*

._---_.

.

•

I if point is 'located outside bottom grid then

I

:---1

I

1

compute wave number and relative frequency

1 1---:---1

I else if depth is negative then 1

I

:---1

1 I

wave number is -l.,relative frequency is

o.

I

1---:---1

else I

:---1

••

I compute function f I

1---1

1

for i= 1 to 50 while f ) accuracy do

I

1 :-~---~---I

I I compute derivative of f I

I

1---1

1 I

compute wave number

1 ,

1---1

1 I compute function f 1

1---:---1

I

compute relative frequency

1

• ---:---:

•

CALL ITWN (IY,IT) Parameters :

.'

IY (I)

IT (I)

) coordinates of point in computational grid

•

Function :

In subroutine VWPRO the group velocity CGO components of wave action transport velocity CXO,

eeo

are determined in a point (IX,IY,IT) computational grid.

and the CYO and in the

• Method :

The relations for the parameters mentioned above used in this model are

•

(19)

•

15 1 D CGO

=

0'0(

----

+

_---

) (1) 2.K0 sinh(2.KO.D) CXO ::: _{CGO.cose + UX} ₍₂₎

CYO ::: _{CGO. sine + UY} ₍₃₎

•

~O iD ~D ~UX ~UX

CSO= - ---(-sine--+cosS--)-ccsS(-sine---+cosS---)

sinh (2.KO.D) ex (Dy ~X àlY

•

~UY (VUY

-sin9(-sine---+cosS---) (4)

ex

àl'I

The terms containing current velocity components UX and UY in eg. 2 through 4 are omitted if no current is present.

• The term CSO is evaluated intermediate lines IX and IX+1. Derivatives of depth and current are determined through a central àifference scheme (after the predictor step).

If negative depths are encountered then all velocity components are given the value 0 ••

•

Structure : *VWPRO~:

:---:

•

if depth is negative then

:---1

1

1 give velocity components and derivatives I

1 the value

o ,

I

---:---else

:---•

I compute group velocity

1---if preàictor step is passeà then

:---1 compute àepth àerivatives and current I derivatives (if current is on)

•

1---:---I compute components of wave action transport

I velocity

I (CSO only if line 1- boundary x=ü)

:---:---•

The callof this subroutine is : CALL VIJPRO(IX,1'1.IT)

Parameters :

•

(20)

-16

•

IX (I) IY (I) IT (I) ) )

coordinates of point in computational grid

Function :

In the subroutine NUMS wave parameters are computed at a new line IX+1 in the computational grid.

• Method :

The following numerical scheme is applied :

Estimates for the wave action AO and -frequency WO at line IX+1 are obtained through a linear extrapolation from the lines IX-1 and IX (predictor step). With these estimates the other wave parameters at line IX+1 can be determined. Linear

• interpolation between the lines IX and IX+1 yields the wave parameters at line IX+1/2, necessary tor the corrector step. The 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.

The amount of energy lost through dissipation (FO) and

• leakage of energy through the boundaries Sa and Sb (FL) is

kept ,

Results of the computations are written to the file REKRES.

•

.:;.

•

(21)

-17 .• Structure : *NUMS*

---:

•

CALL PREDT predictor estimates tor wave act ion

1

and -frequency on line IX+1 1

---1

move contents of arrays with new values of wave 1

number,relative frequency and prcpagation velocity ( components to arrays with old values 1

---1

CALL ~AVPA compute wave parameters at line IX+1 1

---1

determine depths and currents (if current is on) ( intermediate lines IX and IX+1

1 ---1

tor every corrector step do I

• :---1

•

1 àetermine wave parameters intermediate lines 1

I IX and IX+1 1

1---1

1

CALL DISPA compute parameters in dissipation

1

J

1

terms

1

1 1---1

(

1

for every y do

1 I

1 :---1

1 I

1

if last corrector step is in progress then

I

l

:---1

1 I

I

initialize leakage and dissipation

I

1 I

in poin t x,y

1

1 1---:---1

1

1 I

if depth is positive then

1

1 I

1 :---1

I

1

for every

e

do

I

1

1 :---1

I

1 I

1

CALL TERMDE compute terms of the

I

1 I

1

two balance equations

I

1

1 1---1

1

1 I

CALL SUMDE deterrninewave action

I

I I I I I and frequency I

1 [

I

1---:---1

1 I

I

1

if last corrector step in progress thent

I

:---1

1

1 I

I

compute leakage and dissipation in

1

1 I

1

1 I

point x,y

1

I

1---:---:---:---1

1 I

CALL WAVPA compute wave parameters on line IX+11

1---:---1

I CALL ~RIRE write results of line IX+1 to REKRES 1

•

• ._---_.

.

The callof this subroutine is : CALL NUH S (IX)

•

(22)

18

•

Parameter :

IX (I) wave parameters are determined at line IX+1

in the computational grid (X=IX.DX)

1.11

~ubrQyti~ PR~Q!

• Function :

Estimates for the wave action AO and -frequency

wo

at line IX+1 (predictor step) are determined in this subroutine. Method :

The predictor is a simple extrapolation procedure. AO and ~O

.' at line IX+1 are determined as follows :

•

IX+1 IX IX-1

AO = 2.A 0

-

AO (1)

IX+1 IX IX-1

WO = 2.WO

-

\JO (2)

If a negative depth is encountered then AO and \JO are given the value 0 •• • Structure : *PREDT*

._---_.

.

•

if line is boundary

x=o

then 1

:---t

I wave action and -frequency on new line are 1

I

given the values on the old line

1 1---:---1

I else 1

1 ---1

I for every y do I

1 :---1

1 1 if depth is positive then I

I

1 :---1

I I f0r e verY

e

dol

1

1 :---1

I I I move wave act ion and -frequency I I I (to arrays with old values and I

I I I compute new values I

1---:---:---1

1 else I

1 :---1

1 I for every

e

do 1

I

:---1

I

1

move wave action and -frequency I I I to arrays with old values and I I 1 give new values the value O.

.

'

•

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

.

. .

.

•

(23)

..

19

• The callof this subroutine is : CALL PREOT (IX)

Parameter : IX (I)

•

wave parameters are determined at line IX+l in the computational grid (X=IX.OX)

J.18 ~ubrQy!iDg OISP8 Function :

In this subroutine parameters at line IX+1/2 are determined, necessary for the evaluation of the dissipation terms in the

• two balance equations. Method :

The following parameters are determined : orbital velocity at the bottom Ubot current velocity at the bottom Ucur

• wave energy density Et

loçal maximum wave height Hm

the fracticn of breaKing waves Qb

For these parameters the following relations are used :

Ubot 2 2

= (

O&.L

-~_:~~_:~~:~~-& 2 WO •cosh(K 0.0) 1/2 ) (1)

•

Ucur

=

(UX + UY )2 2 1/2 Et

=

De.Ioo.~o

e

(2) (3) -1 1

•

Hm

=

O.BB.KO .tanh(t.KO .0/0 .88)

,

Ka = -- J:.KO (4)

Ne e

1-Qb Et

=

-8.---

(evaluated in FRABRE) (5) 2

.

'

lnQb Hm

•

(24)

,.

20

•

Structure : *DISPA*

---_

.

•

if bottom dissipation is on then I

:---1

1

for every y do

1 I

:---~---I

1 I

compute orbital velocity at the bottom

1 1---:---1

I if current is on then 1

1 :---1

1 for every y do 1

I

:---1

I 1 compute current velocity at the 1

I I' bot tom 1

• ---:---:---:---1

if surf breaking is on then

:---1

•

for every y do I

:---~---I

I compute wave energy density I

1---1

I compute local maximum wave height . 1

1---1

1 CALL FRABRE compute fraction of breakingl

I waves I

• ---:---:---:

• CALL DISPA

1·12 ~!.llrrQill~ rnlillE~

Function :

In this subroutine the fraction of breaking waves in a point

• x.y in the computational grid (Qb) is computed. Method :

A Newton-Raphson iteration process is applied to eq. 1.

•

F - 1-Qb+8.---.lnQb : 0Et

2

(1)

Hm

•

An initial estimate of Qb is obtained from approximation of this implicit relation between 8.Et/Hm**2. The parameters Et and Hm in eq. determined in subroutine DISPA.

a crude Qb and 1 are

•

(25)

-21 • Structure : *fRABRE*

._---_.

.

•

Et 1 if

8.--- )

1.0 then 1 2 1 Hm I

:---1

I fraction of breaking waves Qb is 1 I

:---:---1

•

'I elseI if

8.---

Et

<

.15 then

1

2

1 Hm

1 :---1

I

1

fraction of breaking waves Qb is 0

I

1---:---1

1 eIse I

1 :---(

I I determine an estimate for the fraction of I 1 breaking waves Qb I I

---f

1 compute function f I

I

---1

1 for I

=

1 to 50 while f ) accuraey do I

I

:---1

I 1 eompute derivative of funetion f 1

I

1---1

I I compute fraction of breaking waves Qb I

1

1---1

I

compute function·F

1

•

• ._--

--_._---_.

.

• CALL FRABRE (IY) Parameter :

IY'(I) y-coordinate of point in whieh the fraetion of breaking wa~es Qb is computed

•

Function :

In this subroutine the terms of the t~o evaluated in the point IX+1/2,IY,II. and dWO/dT are split up in eomponents bottom- and surf dissipation.

balans equations are The souree terms 50 of wind generation,

•

(26)

...

22 • Structure : *TERMDE*

---_.

• 0

•

1

CALL TRSY compute transportation terms in

I

1

y-direction

1 1---1

1

CALL TRST compute transportation terms in

1 I

e-direction

1 1---1

I if diffraction is on then 1

1 :---1

I

1

CALL DIFT compute diffraction terms

I---:---~---1 else

I

0---

.

-

---

.

1

give diffraction terms the value O.

J---:---~---1

if wind generation is on then

I

:---1

I

CALL S\HND compute wind generation terms

• 1---:---•

I else

I

:---~---1

I

give wind generation terrnsthe value

o.

1 • __ ---

.

•

1 if bottom dissipation is on then

1 :---1

1 I

CALL SBOT compute bottom dissipation terms

1 1---:---1

1

else

1 I

:---1

1

give bottom dissipation terms the value O.

1 1---:---1

I

if surf breaking is on then

1 I

:---

-

--

-

---1

1 1 CALL SSURF compute surf breaking terms I

1---:---1

1

else.

I

:---1

I I give surf breaking terms the value O.

• :---:---:

•

The callof this subroutine is : CALL TERMDE(IY,IT)

•

Parameters :lY (I)

IT (I)

•

(27)

-23

• Function :

The transportation terms of the two talance equations in y-direct ion :

~ (CYO.AO) and ~ (CYO.WO.AO)

~Y wY

•

are determined in this subroutine. Method :

A conservative central difference scheme is applied :

•

IX+l/2,IY+l,IT IX+l/2,IY-l,IT

f - f

(l ) 2.dY

Energy entering the computational region through the • boundaries y=o and Y=LY is not taken into account. At these

boundaries somewhat different schernes are used.

Structure :

•

*TRSY*

._---_.

.

•

if point is located on boundary and wave energy I

is entering the computational region then I

:---1

1

give transportation terms in y-direction

1 I

1

the value

o , 1

1---:---1

I

else

1

I :---[

I

1

compute transportation terms in y-direction

I

--_._---

. .

.

•

CALL TRSY (IY, IT)

•

ParametersIY (I) :

IT (I)

•

Function :

The transportation terms of the two balance equations in e-direction :

L

cceo ,

AO) and ~_

rceo • wo

.AO)

~e

•

(28)

-24

• are determined in this subroutine. Method :

A conservative central difference scheme is applied :

•

fIX+1/2,IY,IT+1 - fIX+1/2,IY,IT-1

wf: ---w& 2.d&

(1)

Energy entering the computational region through the boundaries &=&a and &:&b is not taken into account. At these • boundaries somewhat different schemes are used. The leakag~

through these boundaries ICeOI.AO.OO is kept.

Structure :

•

*TRST*

._---.

if point is located on boundary and wave energy is entering the computational region then

:---•

lI give transportation terms in a-directionthe value

o.

1---:---else

:---•

I compute transportation terms in 9-direction

1---I if point is located on boundary then

I

:---I I compute leakage of wave energy

:---:---:---•

The callof this subroutine is : CALL TRST (IY,IT)

Parameters :

IY (I)

IT (I)

~•~1

~y!2I.Qu

ti~

!2l[I

Function :

The diffraction terms in the two balance equations :

•

!_(Cdif(CYO.~_AO-CXO.~_AO» éi)& wX wY ~_(Cdif (CYO.~_(WO.AO)-CXO.~_(WO.AO») w& WX éi)Y and

•

(29)

25

• are determined in subroutine DIFT. Method :

Derivatives in x,y and ~ direction are approximated by central difference schemes :

Ö)f IX+1,IY,IT f IX,IY,IT -f

• = ---

(1) Ö)X dX

•

fIX+1/2,IY+1,IT -f IX+1/2,IY-1,IT

= ---

(2) Ö)Y 2.dY

IX+1/2,IY,IT+1 IX+1/2·,IY,Il-1

Ö)f f +f

• -

---

( 3)

2.d9

At the boundaries in the

y-e

plane somewhat different schemes are applied.

• The callof this subroutine is : CALL DIFT (IY,IT)

•

Parameters : IY (I)

IT (I)

coordinates of point in computational grid

Function :

In this subroutine the wind generation components

•

SO and

wind

!i_WO

dt wind

of the source terms in the two balance equations are determined.

•

Method :

The following relations are used for the terms mentioned above :

•

UlO 3 (1) SO

=

Babcd wind 9

•

(30)

'.

26

•

2 9 9 b2 g_WO

=

a2 (----.SO /B), dt 2 3 wind

wind UlO UlO

(2)

• with

2 9

E is the d imensionless wave e nerçy density EO.---- ,EO=O"O.Aa

ij

UlO

•

B is the directional distribution of waves,EO{Y,O)=B(O) .Et (Y) UlO is the wind velocity at lam elevaticn relative

to the current velocity a,b,c and d are coefficients derived from literature

a2,b2 are coefficients to be determined empirically' Eg. 1 and 2 hold only for growing waves (E<aB).In the case E)aB then the wind generation terms are assumed to be

o.

•

Structure : *SWINO*

• :-~---:

• l

if current is on then

I

1 :---1

lidetermine wind speed relative to current I

1---:---1

I compute dimensionless wave energy density E I

l---l

I if E < aB then I

I

:---1

1 1 compute wind generation terms 1

1---:---1

1 else I

I

:---1

1

give wind generation terms the value

o.

• :---:---:

• CALL SWIND(IY.IT) Parameters :

IY (I)

) coordinates of point in computational grid IT (I)

•

Function :

The bottom dissipation terms

•

(31)

27

•

50

bot

and .9._\.10

dt bot

are computed in this subroutine. Method :

• The following relations are applied :

50

=

-101 .aO.AO bot bot (1) dWO 2 =1010 .a3(g -2 3 b3 .1010.101

.co

.AO) bot (2)

•

dt bot with

•

g.KO.CGO 101

=

---.(Cfw.Ubot+Cfc.Ucur) bot 2 2n.WO .cosh(KO.D) (3)

The second term in eq. 3 is omitted if no current is

present. In this formulation the effects of currents on

bottom dissipation are included in the same way dissipation

• due to wave orbital veloeities is determined. The same

procedure with somewhat different relations is applied in

the model CREDIZ. The terms

•

\.Ibot

and Q_WO

dt bot

are determined in point IX+l/2,IY,II in the computational

grid. The wave action AO in the linear term SObot is

included implicitly in the two balance equations (section

3.27) •

• The coefficients Cfw,Cfc,a3 and b3 have to be determined, by

the user, empirically.

CALL SBOT (IY,IT)

•

Parameters :

IY (I)

) coordinates of pOint in computational grid

IT (I)

•

Function :

The terms representing dissipation cf wave energy due to

surf break ing

•

(32)

•

28

•

SO and .9._\olO

dt surf

surf

are determined in this subroutine.

Method :

• Relations for these terms applied in this model are : SO

=

-\ol .rrO.AO

surf surf

(1)

• evo

=

\.l02.a4.(g-2•\.lO.\ol3 •era.AO) b4

dt surf surf ( 2) with 2 surf re Hm

=

O(l.-.Qb.\.lO.---2 Et (3)

•

The terms

•

iJ.surf and Q_\olO dt surf

are determined in point IX+1/2,IY,Il in the computational grid. The wave action AO in the linear term SOsurf is included implicitly in the two balance equations (section

3.27) •

• The coefficient 1 is of order 1 while the coefficients a4 and b4 should be determined empirically.

The callof this subroutine is : CALL SSURF(IY,IT)

• Parameters : IY (I)

•

Function :

In this subroutine the wave action Aa and -frequency \.l0 in the grid point IX+1,IY,IT are determined.

Method :

• The two balance equations 1 and 2 are solved.

•

(33)

•

29

•

50

~ wind

--(CXO.AO)

=

-(transportation+diffraction terms) +

---~X ~O

•

-lW +W ---.---)AO1 dWO (1)

bot surf WO dt

50

~ WO. wind

--(CXO.WO.AO)

=

-(transportation+diffraction terms) +

---• ~X

ao

- (W + W )WO .AO (2)

bot surf For brevity the transportation- and diffraction terms as

• weIl as the source term dWO/dt have net been written in full in eg. 1 and 2.

The numerical scheme applied to these eguations is:

IX+l,IY,IT IX+l,IY,IT IX,IY,IT IX,IY,IT

• CXO f - CXO f dX

•

H IX+l,IY,IT IX,IY,IT

=

G - -Cf + f ) 2 (3) with :

f represents the terms AO (eg.l) or WO.AO (eg.2)

G contains the non-linear terms of eg. 1 and 2 H contains the linear terms of eg. 1 and 2

• Further the dissipation of wave energy (Wbot+Wsurf)AO.OO.de is determined in point IX+l/2,IY,IT in the computational grid. Structure :

• SUMDE

._---_.

.

•

I determine wave action and -freguency I

1---1

I if last corrector step is in progress then (

I

:---1

I I compute dissipation of wave energy

._-_._---_.

. .

.

The callof this subroutine is : CALL SUMDE (IY,IT)

•

(34)

•

30

I

.

1 I I Parameters : IY (I)

•

Function :

In the subroutine ~RIRE arrays containing wave parameters at a line in the computational grid are written to the file REKRES.

•

Structure : *lolRIRE*

:---:

•

1 if line 1 boundary x=O then

I

:---1

I I write leakage of energy fL,dissipation of 1 1 1 energy FD,fraction of breaking waves QB and thel

1

wave action transport velocity component C90

I

J J to the file REKRES 1

1---:---1

( write depth D,current velocity components UX,UY (ifl

1 current is off tnen f111 arrays with 0) ,wave act ionI I AO,-frequency ~O,relative frequency ~O,wave number I I KO,group velocity eGO and wave action transport I ( velocity components eXO,CYO to the file REKRES I

• :---:

•

The callof this subroutine is : CALL ~RIRE (IX)

Parameter :

IX (I) results of line IX in the computational grid

• are written to the file REKRES

•

This subroutine, called at the start of every subroutine and the main program, provides a message of the entry of this subroutine resp. program.

STRACE is copied from the model CREDI2.

The subroutine MSGERR, called when an error is encountered

• during tne executä on of the program eCHPU, provides an error message.

MSGERR is copied from eREDI2.

•

(35)

31

• In subroutine COPYCH character strings are copied to real variables and back. It is copied from CREDIZ.

•

(36)

•

• 4. STORAGE OF DATA

32

As in CREDIZ a dynamic data pool is used to obtain an

I

.

efficient and flexible sto r a qe , cf arrays. \.Jith the subroutine ADPOOL the dimension of an array can be extended or reduced. The structure of the pool is the same as in CREDIZ.

An element of an array A is found by :

A(I)=POOL(IA+I), IA is the adres of array A A two-dimensional array is stored row by row :

• eg. : array B(1:n,l :m) storage in pool :

8 (1,1), ••• ,B (n,1). ,8 (1,2) , ••• ,B (n ,2) ,B (l,m), •••,B (n ,m)

thus B(k,I)=POOL(IB+(I-l)n+k), IE is the adres of array B For more detailed information on the structure of the dynamic data pool, reference is made to the system

• documentation of CREDIZ.

The following arrays of the program COMPU are included in the pool :

~--~---r---~---I numberl name ( adres ~--~---r---~---I description

• t---t---t---t---

· ---

----•

1 2 3 4 5 6 1 8 9 10

•

11 12 13 14

•

15 16

•

I I 11 1

1

18 I I 19 I 20 I 21 I 22 I 23

•

I DEP I VX I VY I \.JAO I \.JFO I \.JKO I RFO I CGO I CXO I I CYO I I \.JA t \.JF I \.JK I RF I I CG I CX I I CY I ( CT I I \.JAN I WFN I \.JKN I RFN I CGN I lDEP I IVX

1

IVY I l\.JAO I I\.JFO (,l\.JKO I IRFO I ICGO I ICXO I I ICYO I I I\.JA I I\.JF I IWK I IRF I I lCG I ICX I I ICY I I ICT I I IWAN ( I\.JFN I H1KN I IRFN ( ICGN I depths

( x-component current velocity I y-component current velocity

I wave action (old line) I wave frequency (old line) I wave number (old line)

I relative frequency '(old line) 1 group velocity (old line) 1 x-component wave action

I transport velocity (old line) I y-component wave action

I transport velocity (old line)

I wave action (between oid and new line) I wave frequency (between old and new line) I wave number (between old and new line) I relative frequency (between oid

land new line)

I group velocity (between old and new line)

I x-component wave action transport

I velocity (between old and new line) I y-component wave action transport I velocity (between old and new line) I a-component wave action transport

I velocity (between old and new line)

I wave action (new line) I wave frequency (new line) I wave number (new line)

I relative frequency (new line)

(37)

I CXN I ICXN x-component wave action

I I transport velocity (new line)

I CYN I ICYN y-component wave action

I I transport velocity (new line)

I DO 1'100 derivatives in bottom geometry l DC I IOC derivatives in current field I HM I IHM local maximum wave height I QB I IQB fraction of breaking waves

I ET I IET directionally integrated wave energy

t

FD I IFD dissipation of wave energy I FL I IFL leakage of wave energy

I UBOT I IUBOT orbital velecity near the bottom I UCUR I IUCUR current velocity near the bottom I DEO I IDEO depths (old line)

I UXO I IUXO x-component current velocity

I I (old line)

1

UYO

I

IUYO y-component current velocity

L I (old line)

IDEM I IDEM depths (between old and new line) 1 UXM r IUXM x-component current velocity

I I (between old and new line) I UYM I IUYM y-component current velocity

I I (between ci

e

and new line) I DEN f IOEN depths (new line)

1

UXN

I

IUXN x-component current velocity

I I (nev Iine)

I UYN I IUYN y-component current velocity

I I I (new line) _____ ~ ~ .L .__ . _

•

I 24 I l 2S I I 26 I 27

1

28 I 29 I 30 I 31 1 32 I 33 , 34 I 3S I 36

•

37

•

38 39 40

•

4142 "3 33

•

table 1 arrays in dynamic data pool

4.2 Common blocks A number of common blocks are defined in whIch-principal~ata for the model HISWA is included.Each of these blocks contains a certain categcry of information :

•

(38)

•

34

• r---T---,

r name I description I

~---t---

· ---.-1

•

I TITEL lOEPROS I REKROS I I TRANSf I 1 NUMS l TERMOE I TESTOA I I FYSPAR I UITVOA I I REFNRS I LEESOA I POOL I

•

I titLe

llocation and dirnensions of bottom grid llocation and dimensions

I of computational grid

I transformation coefficients between I different grids

I information on the numerical scheme I terms of the balance equations

I information for output control I (especially the tracing of errors) I physical parameters

I information tor the ~rogram OUTP I (not used by the program COMPU) I data set reference numbers

I information for reading data I references to arrays in the

r dynamic data pool

L . .~ .1

tabl~ 1 description of common blocks

• Most of the common blocks are copied from the model CREDIZ. In the following only the changes with regard to the common blocks in CREDIZ are discussed. for more information reference is made to the systerndocumentation of CREDIZ. TITEL

•

DEPROS as in CREDIZthe elements AKX,CCGX and WKCX are omitted REKROS

•

TRANSf

.

'

NUMS

•

the following elements are added :

MTR number of grid points in a-direction TETAA

) boundaries of computational grid in &-direction TETAB

OT grid size in a-direct ion as in CREDIZ

contains the following elements : NCOR number of corrector steps IPRE indicator of predictor step ICOR indicator of corrector step

IOaW option for the representation of the boundary condition at x=O

ICUR switch for the introduction of current IOIf switch for the introduction of diffraction

IBOT switch for the introduction of bottom dissipation ISURF switch for the introduction of surf breaking WOlP wind direction relative to problem grid

(39)

•

35

.

'

UlO wind velocity at 10 m. elevation

U10C wind velocity at 10 m. elevation relative to current ADIR coefficient for directional distribution of waves AT (300) wave action at boundary x=O

WT(300) wave frequency at boundary x=O BDIR(30} directional distribution of waves CDIF diffraction coefficient

PWIND(10} parameters of wind generation (a,b,c,d,a2,b2) PBOT(5} parameters of bottom dissipation (Cfw,Cfc,a3,b3) PSURF(5} parameters of surf breaking ( i, ,a4,b4)

TERMDE contains the following elements :

• FYA

•

) transportation terms in y-direction FYF

FTA

) transportation terms in a-direction

•

FTFDIFA

) diffraction terms DIFF

~INDA

) wind generation terms

•

WINDFWBOT

} bottom dissipation terms BOTF

WSURF

) surf breaking terms SURFF

• TDIS dissipation of wave energy TLEAK leakage of energy

TESIDA as in CREDIZ

FYSPAR the following elements are omitted :

• IM,DEP1,DEP2,DEP3,UX1,UX2,UX3,UY1,UY2,UY3, AK1,AK2,AK3,CCG1,CCG2,CCG3,WKC1,UKC2,WKC3, ISTA1,ISTA2,ISTA3,SIGMA1,SIGMA2,SIGMA3, SINH1,SINH2,SINH3,AMPL1,AMPL2,AMPL3 the following elements are added :

DEP depth in a point in the bottom grid

• UXC

} current velocity components relative to computational gri UYC

UNU wave number

•

UITVDAREFNRS as in CREDIZthe elements HULPF1 and HULPF2 are omitted(for the time being) LEESDA as in CREDIZ

POOL this block will be adjusted to the new construction

•

(40)

•

36

•

of the dynamic data pool·(section 4.1)

In the model HISWA a number of files are used that serve as communication tools between the programs PREP, COMPU and

• OUTP. Two of these files are used by the program COMPU : a) INSTR

This file contains instructions formulated in PREP and data necessary for the computations carried out in the program COMPU.

contents :

• NPOOL number of pool arrays ITEST test parameter

common bLocks :

TITEL, DEPROS, REKROS, IRANSF,

NUMS, FYSPAR, REFNRS, UIIVDA, TESTDA arrays :

•

DEP,VX,VY,OUIREQ,OUTOA

b) REKRES

In file REKRES the results of the ccmputations carried out in the program COMPU are stored. This file will be read by the program OUIP.

contents :

•

common b Lo cks

TITEL, OEPROS, REKROS, TRANSF, NUM S, FYSPAR , REFNRS

arrays

..

DEP, VX, VY

DEN, UXN, UYN, WAN, WFN, RFN,

line x=O RFN, WKN, CGN, CXN, CYN

fL, FO, QB, CT ) line x= (i-1/2) .dx ) )

DEN, UXN, UYN, WAN, WFN, RFN,

)i=I,t-'XR-) line x=i.dx )

RFN, WKN, CGN, CXN, CYN )

•

Each array in the files described abcve is preceded by the dimension of the array. Ihe same conventions as described

• under 4.1 are applied with regard to he storage of two-dimensional arrays. Files are read and written unformatted.

•

(41)

•

37

• REFERENCES

Abramowitz,M. and Stegun,I.A.,1968,

Handbook of mathematica1 functions,Dover publications

• CREDIZ01 system documentation,1984 (in Dutch)

Holthuijsen,L.H. and Booij,N.,1983,Selection and formulation of a numerical shallow water wave hindcast model,Delft University of Technology,Report no 1783