<%BANNER%>
  • TABLE OF CONTENTS
HIDE
 Title Page
 Front Matter
 Table of Contents
 Introduction and overview of the...
 Wave model in greater detail
 Circulation model in greater...
 Variations of the program
 Variables used in the computer...






Documentation of the wave-induced circulation model
CITATION SEARCH THUMBNAILS PDF VIEWER PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00073939/00001
 Material Information
Title: Documentation of the wave-induced circulation model
Physical Description: ii, 39 leaves ; 28 cm.
Creator: Winer, Harley Stanford 1944- ( Author, Primary )
Publisher: Civil and Coastal Engineering, University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: August 1988
 Subjects
Subjects / Keywords: Coastal Engineering
Aerospace Engineering, Mechanics, and Engineering Science thesis Ph. D.
Dissertations, Academic -- Aerospace Engineering, Mechanics, and Engineering Science -- UF
Tidal currents -- Mathematical models
Wave equation
Waves -- Mathematical models
 Notes
Funding: This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
 Record Information
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 24655525
System ID: UF00073939:00001

Downloads

This item has the following downloads:

( PDF )


Table of Contents
    Title Page
        Title Page
    Front Matter
        Front Matter
    Table of Contents
        Table of Contents
    Introduction and overview of the program
        Page 1
        Common statements
            Page 1
            Page 2
            Page 3
        Brief description of the Subroutines
            Page 4
            Page 5
            Page 6
            Page 7
            Page 8
            Page 9
    Wave model in greater detail
        WAVRAD subroutine
        Page 10
        AMPCO Subroutine
            Page 11
        J=2,N DO LOOP
            Page 12
        J=1 and J=N Boundaries
            Page 13
        RAD Subroutine
            Page 14
            Page 15
    Circulation model in greater detail
        Page 16
        XCOEF subroutine
            Page 17
            Page 18
        YCOEF subroutine
            Page 19
            Page 20
        CIRC subroutine
            Page 21
            Page 22
    Variations of the program
        Page 23
        WINER.BARR directory
        Page 24
        Page 25
        Page 26
        WINER.GROIN directory
            Page 27
            Page 28
        WINER.CLOSED directory
            Page 29
        WINER.PLOT directory
            Page 30
            Page 31
    Variables used in the computer code
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
        Page
Full Text



UFL/COEL/MP-88/1


DOCUMENTATION OF THE
WAVE-INDUCED CIRCULATION MODEL











By

Harley S. Winer


Version 1.0


August, 1988












DOCUMENTATION OF THE
WAVE-INDUCED
CIRCULATION MODEL




by




HARLEY S. WINER


Version 1.0
August 1988
















TABLE OF CONTENTS




1 INTRODUCTION AND OVERVIEW OF THE PROGRAM ........... 1

1.1 Common Statements .................... ........... 1

1.2 Brief Description of the Subroutines ........ ........ .... 4

2 THE WAVE MODEL IN GREATER DETAIL ...... . ....... 10

2.1 The WAVRAD Subroutine .................... ........ 10

2.2 The AMPCO Subroutine .................... ........ 11

2.2.1 The J=2,N DO LOOP .................... ...... 12

2.2.2 The J=1 and J=N Boundaries ................ .... 13

2.3 The RAD Subroutine .................... .......... 14

3 THE CIRCULATION MODEL IN GREATER DETAIL ............. 16

3.1 The XCOEF Subroutine .................... ... ...... 17

3.2 The YCOEF Subroutine ... .......... ............. 19

3.3 The CIRC Subroutine ..... ............. ...... .... 20

3.4 The FLOOD subroutine ............................... 21

4 VARIATIONS OF THE PROGRAM .... ....... ............... 23

4.1 The (WINER.BARR] Directory .................. ....... 23

4.2 The [WINER.GROIN] Directory ......................... 27

4.3 The [WINER.CLOSED] Directory ......... ............. 29

4.4 The [WINER.PLOT] Directory .. .............. .. .. .. 30

APPENDICES

A VARIABLES USED IN THE COMPUTER CODE ......... ...... 32















CHAPTER 1
INTRODUCTION AND OVERVIEW OF THE PROGRAM



This report is to document the programs that I have developed for my dissertation.

There are actually several separate programs in separate directories. They are named PRO,

GROIN, and BREAK. Each of these programs is written for a Cartesian grid. The PRO

program has no internal barriers. GROIN is an adaptation of PRO so as to include the

internal boundary conditions that would exist at a shore perpendicular barrier. Similarly,

BREAK is an adaptation of PRO for the case of a shore parallel breakwater. In both cases

the barrier is represented as infinitesimally thin and impermeable. The two adaptations

GROIN and BREAK will be described after the mechanics of PRO are explained. Also an

adaptation for closed basins is included.

The main program PRO calls many subroutines which in turn call other subroutines.

All of the subroutines used are listed in the file PRO.OPT which is also used for the linking

command. The linking command is

$LINK PRO/OPT

Each of the subroutines listed in PRO.OPT will be described separately. For reference the

file PRO.OPT is printed below.
PRO,INPUT,TRIDA,XCOEF,XCOJ1,XCOJN,YCOEF,YCVAR,EXCOEF
UP,CHECK,PRINT,WAVRAD,WAVNUM,AMPCO,DISP,FLOOD,SHEAR2
CTRIDA,RAD,WDIS,CIRC,AMPOUT,WRITE

1.1 Common Statements

Each of the subroutines has as the second line

INCLUDE 'DIM.FOR'















CHAPTER 1
INTRODUCTION AND OVERVIEW OF THE PROGRAM



This report is to document the programs that I have developed for my dissertation.

There are actually several separate programs in separate directories. They are named PRO,

GROIN, and BREAK. Each of these programs is written for a Cartesian grid. The PRO

program has no internal barriers. GROIN is an adaptation of PRO so as to include the

internal boundary conditions that would exist at a shore perpendicular barrier. Similarly,

BREAK is an adaptation of PRO for the case of a shore parallel breakwater. In both cases

the barrier is represented as infinitesimally thin and impermeable. The two adaptations

GROIN and BREAK will be described after the mechanics of PRO are explained. Also an

adaptation for closed basins is included.

The main program PRO calls many subroutines which in turn call other subroutines.

All of the subroutines used are listed in the file PRO.OPT which is also used for the linking

command. The linking command is

$LINK PRO/OPT

Each of the subroutines listed in PRO.OPT will be described separately. For reference the

file PRO.OPT is printed below.
PRO,INPUT,TRIDA,XCOEF,XCOJ1,XCOJN,YCOEF,YCVAR,EXCOEF
UP,CHECK,PRINT,WAVRAD,WAVNUM,AMPCO,DISP,FLOOD,SHEAR2
CTRIDA,RAD,WDIS,CIRC,AMPOUT,WRITE

1.1 Common Statements

Each of the subroutines has as the second line

INCLUDE 'DIM.FOR'







2

where DIM.FOR is a file that has all of the dimension statements that are used as common

statements. The reason for having all of the common statements in a separate file is to

facilitate changing the size of the domain. This will be useful if the program is to be run

on a larger computer for a larger domain. At present the dimension of the domain is 100

by 100. This is because the VAX 750 has limited storage.

The DIM.FOR file is printed below.
COMPLEX AA,BB,CC,DD,VVV,AMP,AMPO
COMMON/A/U(100,100),UNEW(100,100)
COMMON/B/ETA(100,100),ETANEW(100,100),ETADT(100,100)
COMMON/C/V(100,100),VNEW(100,100)
COMMON/D/H(100,100),TD(100,100),TD2(100,100)
COMMON/E/A(200),B(200),C(200),D(200),VV(200)
COMMON/H/WAVEHGT(100,100)
COMMON/I/TBX(100,100),TBY(100,100)
COMMON/J/SXX(100,100),SXY(100,100),SYY(100,100)
COMMON/K/G,DT,DX,DY,FRIC,T,RO,IWET(100)
COMMON/L/AA(100),BB(100),CC(100),DD(100),VVV(100)
COMMON/LL/AMP(100,100),AMPO(100)
COMMON/M/RK(100,100),RKO(100),CG(100,100),SIG(100,100)
COMMON/M/KBREAK(100)
COMMON/O/EX(100,100),EY,EMAX
COMMON/P/THETA(100),THETAH

Appendix A lists each of these variables along with many other variables in alphabetical

order with a brief description.

The wave model part of the program solves for the complex amplitude at each of the

grid points. The complex amplitude array is named AMP. The complex amplitude on the

first grid row is named AMPO. The AMPO array is used so as to enable a smooth start-up.

AMP on the first row is given by a time start function times AMPO. This is explained on

page 59-60 of my dissertation. The arrays AA, BB, CC, and DD are the column vectors

for the tridiagonal matrix equation obtained from the parabolic wave equation. These are

complex. VVV is the solution vector for the tridiagonal matrix equation obtained from the

parabolic wave equation. AA, BB, CC, DD, and VVV which are complex are distinct from

the column arrays A, B, C, D, and VV which are the column vectors and solution vector

for the non-complex tridiagonal matrix equation obtained in the circulation model.








3

The circulation model solves for the mean water surface elevation and the x and y

components of the depth integrated mean current. U and V are the values of the x and y

components of the current and ETA is the mean surface elevation measured upwards from

the still water line. UNEW, VNEW and ETANEW are the new values of the current and

surface elevation which are determined by the CIRC subroutine. ETADT is the time rate

of change of ETA. This is used at the off shore boundary for one of the non-linear advective

acceleration terms. This will be explained later. A, B, C, and D are the column vectors for

the tridiagonal matrix equation obtained from the governing circulation equations. VV is

the solution vector. It should be noted that A, B, C, D, and VV have a dimension that is

twice that of the variables U, V, and ETA. This is because the tridiagonal matrix equation,

which is solving for a staggered grid system, solves for both U and ETA in the x direction

and for V and ETA in the y direction.

H is the depth measured downward from the still water line and TD is the total depth

equal to H plus ETA. TD2 is also the total depth and is equal to TD when ever TD is

positive. If TD is zero or less than zero, TD2 is set to an arbitrary value of 1 millimeter

of water. The TD2 array is used in the wave model so as to avoid the necessity of internal

boundary conditions.

WAVEHGT is the wave height which is actually twice the absolute value of the complex

amplitude AMP. The WAVEHGT array is used in the SHEAR2 subroutine to determine

the bottom shear stresses which are designated by the array names TBX and TBY. TBX is

named for Tau Bottom in the X direction and TBY is for Tau Bottom in the Y direction.

Two separate arrays are used since each is positioned at a different grid edge.

SXX, SXY, and SYY are the radiation stress components which are computed in the

RAD subroutine. RK is the wave number, CG is the group velocity, SIG is the relative

frequency, and RKO is an averaged wave number which is obtained by averaging RK across

a grid row. RK, CG, SIG, and RKO are calculated in the subroutine WAVNUM. THETA is

the wave angle at each grid row assuming a plane beach without currents as calculated by








4

Snell's law. THETAH is the wave angle calculated by Snell's law for the position half way

between grid centers ( i.e. at the grid edge) at the location of the breakwater. THETAH is

used only in the BREAK program.

G, DT, DX, DY, FRIC, T, and RO are the gravitational acceleration constant, delta

t, delta x, delta y, the friction factor, the wave period and p the water density respectively.

IWET is an array which indicates the position of the shoreline for each of the grid columns.

This array is used to allow the shoreline to move with raising or lowering of the water

level. This mechanism will be discussed further in the section (3.4 explaining the FLOOD

subroutine.

KBREAK is an index which is used in the WDIS subroutine to indicate whether or not

dissipation due to breaking should be applied. This is discussed further in the description

of the subroutine WDIS below.

EX, and EY are eddy viscosity coefficients which are calculated in the subroutine EX-

COEF. EMAX is the maximum value that these coefficients may be. EMAX is used so as

to place an upper limit on the eddy viscosity coefficients and thus avoid numerical stability

problems that may arise from too large a coefficient. This is important since the lateral

mixing terms are represented explicitly in the model and thus place a constraint on the size

of either the delta t increment or the lateral mixing coefficient.

1.2 Brief Description of the Subroutines

INPUT(M,N,WAVAMP,THETAO,TIME,ILAST) is called by the main program. It is

called only once at the beginning of the run. This subroutine prompts for the name of

the input file, reads the file to set the size of the computational domain M and N the grid

spacing DX and DY, the time step increment DT, the start-up time constant TIME, the

incident wave height HO at the offshore grid row, the wave period T, the wave angle at

the offshore grid row THETAO, the friction coefficient FRIC, the maximum for the lateral

mixing coefficient EMAX, and the maximum number of iterations ILAST. For the special

programs GROIN and BREAK the input data file also contains the position of the structure.








5

All of the above are separated by commas in the first line of the input file, so that the read

statement is in free format. A sample input file is IN.PL in the [WINER.MODELJ directory.

The subroutine then prompts for the name of the depth file, opens it and reads the M by N

array of the depth H. H is measured positive downward from the still water line. Thus land

above the still water line will be represented by a negative value. It is important that M be

large enough to allow for flooding of the shoreline as set-up is produced. After reading the

H depth array the subroutine determines the values of IWET(J) for all values of J=1,N.

This is simply done by determining for each grid column the first grid from the shoreline

to have a positive depth. Values for IMIN and IMAX the minimum and maximum of the

IWET array are also determined. The values of IWET(J), IMIN, and IMAX determine the

limits of the actual computation in the circulation model. The program updates the values

of IWET, IMIN, and IMAX with each iteration by calls to the FLOOD subroutine from

the main program. The subroutine initializes the values of U, V, and ETA to be zero. The

final task of the INPUT subroutine is to determine the values of AMPO. This is done as

follows. First a call is made to the subroutine DISP to determine the wavenumber on the

offshore grid row. Then the real part of AMPO is given by HO cos(Y RK1 cos(THETAO))

and the imaginary part by HO sin(Y RK1 cos(THETAO)), where RK1 is the wave number

along the first grid row assuming no currents, and Y is the distance along the y coordinate

of the grid center from the reference y=0 position, and HO and THETAO are the wave

height and wave angle at the offshore grid row.

TRIDA(IF,L) is a subroutine that is called by the CIRC subroutine to solve the tridi-

agonal matrix equation using a double sweep method. IF and L are the first and last indices

of the range of the column vectors A, B, C, and D which are passed in common statements.

IF is usually 1 and L is usually 2*IWET(J) or 2*N+1.

XCOEF(M,J) is called by the CIRC subroutine to determine the A, B, C, and D column

vectors for each grid column in the x direction part of the alternating direction implicit

method solution of the governing circulation equations. It is called sequentially for J=2,N-








6

1. It determines all of the coefficients from I=1 to I=IWET(J). Special adaptations of

XCOEF are made for any case where the J grid column is adjacent to either the domain

boundary or an internal boundary imposed by the presence of a barrier. The A, B, and C

arrays are the coefficients of the tridiagonal matrix equation, this is the logic for the naming

of the XCOEF subroutine.

XCOJ1(M,J) is a special case of XCOEF for J=1. This subroutine incorporates the

lateral boundary conditions.

XCOJN(M,J) is also a special case of XCOEF but for J=N. Again, this subroutine

incorporates the lateral boundary conditions.

YCOEF(I,N) is called by the CIRC subroutine to determine the A, B, C, and D column

vectors for each grid row in the y direction part of the alternating direction implicit method

solution of the governing equations. This subroutine is called for values of I less than IMIN.

It determines the coefficients for the I grid row from J=1 to J=N.

YCVAR(I,N) is an adaption of the YCOEF subroutine for the case of I equal to or

greater than IMIN. The adaptation is necessary so as to include the conditions that are

imposed by adjacent dry grids. (Dry is used to refer to grids where the depth is such

that it is above the water level. Similarly the term wet is used to refer to grids that have

water.) This subroutine is also called from the CIRC subroutine. This subroutine is full of

if statement to determine the relationship of the grid to the shoreline. If the neighboring

grid in the y direction is dry then the V velocity between the two is set to zero. If the I+1

grid is dry then the edge of the grid is considered the shore and three point derivatives are

taken using a zero value at the shoreline. (This is done for x derivatives of V and SXY).

EXCOEF(M,N) is the subroutine that determines the lateral mixing coefficients. It is

called from the main program.

UP(M,N) is called from the main program. It does the simple task of updating the values

of U, V, and ETA setting these values equal to those of UNEW, VNEW, and ETANEW.

This subroutine also computes the ETADT terms for the first two grid rows.








7

CHECK(KFLAG,M,N) is called from the main program prior to the call to UP. It

checks for convergence by comparing the new values (UNEW, VNEW, ETANEW) with the

old values (U, V, ETA). Convergence is arbitrarily set to be when the percentage difference

between the two is less than an arbitrarily small epsilon EPS for each grid and for each

unknown U, V, and ETA. The KFLAG term is a flag that the main uses to decide whether

to continue for another iteration or to exit the program.

PRINT(M,N) is a subroutine that prints the values of U, V, and ETA. This subroutine

was used primarily in the development stage.

WAVRAD(M,N,THETAO,WAVHT) is the wave model subroutine. It is called from

the main program. It in turn makes calls to the subroutines WAVNUM, DISP, AMPCO,

CTRIDA, and RAD. The purpose of this subroutine (and those that it calls) is to determine

the radiation stress terms. This subroutine is the same as WAVE shown in the flow chart

of the model on page 11 of the dissertation. WAVHT is the start-up function C(t) times

WAVAMP which is one half the incident wave height. WAVHT times AMPO(J) gives the

value of AMP(1,J) to start the wave model marching.

WAVNUM(M,N,THETAO) is a subroutine that calculates the wavenumber RK, the

averaged wave number across a grid row RKO, the group velocity CG, and the relative

frequency SIG. A Newton method is used to solve for RK using equation (3.76) on page

44 of the dissertation. RKO is obtained by a simple averaging while CG and SIG are

directly obtained from the appropriate formulas. The WAVNUM subroutine also solves for

the THETA array using Snell's law. To do this calls are made to DISP to determine the

wavenumber in the absence of any currents. The assumption is implicit in the model that at

the lateral boundaries there are no significant on-offshore currents and that the topography

has no y variations. Thus Snell's law for refraction applies at the lateral boundaries.

DISP(HH,RRKK) solves for the wavenumber in the absence of currents. HH is the

depth of the water column and RRKK is the wavenumber that is returned.

AMPCO(I,N) determines the coefficients of the tridiagonal matrix equation which re-








8

suits from the finite differencing of the parabolic wave equation. The governing equation

used is equation (4.26) on page 52 of the dissertation. It is finite difference using a Crank-

Nicholson scheme and is written as equation (4.27) on page 53 of the dissertation. AMPCO

is called sequentially for I=1,IMAX-1 from the WAVRAD subroutine.

CTRIDA(IF,L) solves the tridiagonal matrix equation using a double sweep algorithm.

CTRIDA is essentially the same as TRIDA except that it is for complex variables.

RAD(M,N) calculates the radiation stress components SXX, SXY, and SYY from the

gradients in the complex amplitude AMP using equation (2.57) on page 29 of the disserta-

tion.

WDIS(I,J,TTTDD,CCGG,AATEMP,WAVDIS) is called from the AMPCO subroutine.

It computes the value of the dissipation coefficient in the governing wave equation using

equation (3.97) on page 47 of the dissertation. TTTDD is the total depth and CCGG is

the group velocity. These values are passed in the argument rather than through common

statements because of the need to compute the dissipation at the grid edge in the BREAK

program. AATEMP is complex and is the value of the complex amplitude on the I row or

the value of the complex amplitude calculated explicitly on the 1+1 grid by use of equation

(4.29) on page 54 of the dissertation.

CIRC(M,N) is the circulation model. It is called directly from the main program. It

makes calls to XCOJ1, XCOEF, XCOJN, YCOEF, YCVAR, and TRIDA.

AMPOUT(M,N) writes out the absolute value of the complex amplitude to a data file

named AMP.LST and writes out the value of the instantaneous water surface elevation due

to the wave field to another data file named ETA.LST.

FLOOD(M,N) determines whether the boundary of the computation domain at the

shoreline should be increased or reduced due to changes in the water level. This is done

with a lot of IF statements to determine whether a grid has gone dry or if the elevation

in a neighboring grid is sufficient to wet the shoreward grid. The flooding procedure is

done only in the x direction along grid columns. The procedure is done in such a way as








9

to insure continuity of water mass in the domain. When it is determined that the water

surface elevation in one grid is sufficient to flood the next grid an arbitrarily small amount

of water is taken from the one grid and placed in the flooded grid. This is done so as to

avoid a zero water depth in the computational domain (the wave routine will react to a zero

water depth). If the depth in a grid becomes negative, then the deficient amount is taken

from the neighboring grid to produce a zero water column and the grid is removed from the

computation domain.

SHEAR2(M,N) calculates the TBX and TBY arrays using the formulation given by

equations (2.15) and (2.16) on page 20 of the dissertation.

WRITE(M,N) creates a file OUT.LST and writes out the ETA, U and V arrays. These

values are written in free format. First ETA is written for J=1,N for each grid row from

I=1 through I=IMAX. Then U is written the same and then V is written for J=1,N+1 for

each grid row. The values of the IWET array are also written so that the shoreline and

shore values can be found.















CHAPTER 2
THE WAVE MODEL IN GREATER DETAIL



2.1 The WAVRAD Subroutine

The wave model part of the program is the WAVRAD subroutine. It makes calls

to the following subroutines: WAVNUM, AMPCO, CTRIDA, and RAD. The WAVNUM

subroutine in turn makes calls to DISP and the AMPCO subroutine makes calls to WDIS.

Within the WAVRAD subroutine the first do loop is to set the KBREAK array equal

to zero. This array is used by the WDIS subroutine to determine whether dissipation

due to wave breaking should be included. If KBREAK(J) is equal to zero than the wave

dissipation is set to zero. If KBREAK(J) is one than dissipation is a positive value. The

value of KBREAK is determined by the WDIS subroutine. The value remains at zero until

the wave height (twice the absolute value of the complex amplitude) becomes equal to or

greater than .78 time the total water depth, and is than set to one. Once set to one, the

value of KBREAK(J) remains one until the wave height becomes less than .4 times the water

depth. The rational for this is found in section 3.3 starting on page 45 of the dissertation

or the reader can go to the original work of Dally cited in the dissertation.

Within the same do loop setting KBREAK equal to zero the value of AMP on the first

grid row is determined by the expression AMP(1,J)=WVAMP*AMPO(J), where WVAMP

is passed from the main program. Some clarification is needed concerning the use of several

variables such as WAVAMP and WVAMP and the array WAVEHGT. WAVAMP is equal

to one half of the input incident wave height. It is set by INPUT and remains unchanged

for the entire program. WVAMP is determined in the main and is the time function C(t)

of equation (4.45) on page 60 times WAVAMP. The time function C(t) has the name UU in

the main. WAVEHGT is the array of the wave height determined in the RAD subroutine








11

and is merely twice the absolute value of the complex amplitude AMP.

The WAVRAD subroutine then calls WAVNUM to calculate the RK, RKO, CG, THETA,

and SIG arrays. After the call to WAVNUM there is the following do loop
DO 20 I=1,IMAX-1
CALL AMPCO(I,N)
CALL CTRIDA(1,N) (2.1)
DO 20 J=1,N
20 AMP(I+1,J)=VVV(J)
This do loop is the meat of the wave model. It is the marching grid row by grid row

of the parabolic equation. The AMPCO subroutine determines the column vectors AA,

BB, CC, and DD of the tridiagonal matrix equation. The AA, BB, CC, and DD vectors

are determined very straight forward. Grouping terms, equation (4.27) on page 53 of the

dissertation is written as

AA(J) A'1 + BB(J) A* 1 + CC(J) A'- = DD(J) (2.2)

This equation then is the definition of the AA, BB, CC, and DD column vectors. CTRIDA

then solves the tridiagonal matrix equation producing the solution vector VVV. And then

the solution vector is assigned to the AMP values on the next or I+1 grid row until values

of the AMP array are determined for every grid in the IMAX by N domain.

The last line in WAVRAD before the return to the main program is a call to RAD

to calculate the radiation stress arrays according to equation (2.57) on page 29 of the

dissertation.

2.2 The AMPCO Subroutine

As explained above the AMPCO subroutine determines the AA, BB, CC, and DD

column vectors of the tridiagonal matrix equation resulting from the finite differencing of

the governing parabolic equation. The values of these vectors are determined for each grid

1 through N. Special consideration is given for the J values of 1 and N so as to incorporate

boundary conditions. Thus the subroutine has one set of computations for J=1, then a

do loop for J=2,N and then another computation for J=N. The do loop for J=2,N will be

examined first.










2.2.1 The J=2,N DO LOOP

The first eight lines of the do loop are to obtain spatially centered velocities. This is

necessary since the velocities U and V are at the grid edges while all other variables are

at grid centers. The eight values determined are VIJ, VIP1J, VIJP1, VIP1JP1, VIJM1,

VIP1JM1, UIJ, and UIP1J. The logic of these code names is the velocity and the position.

For example, VIP1JM1 is the V velocity at the I plus 1 and J minus 1 grid (VI0 ), while

UIP1J is the U velocity at the I plus 1, J (U+l) grid center. These values are the average

of the two values at the opposite grid edges. On the next two lines CGPU is [C, cos 0 + U]}

and CGP1U is [C, cos 0 + Uj+1.

The next several lines compute some intermediary terms which are used for short hand

purposes (not necessarily for computational efficiency). These terms are listed below

A2 = I
A2= (C, cos 0 + U)


B2 =
C, cos +UI s+1 + (C Cos +U

oC2 =o7'

C cos 0 +U). 4 DY
DY'


D2 a M+1 ^ -1
2 = C,cos + U 4DY

c i+1 i+1
PAA = s+l + D1
2 Dy2(C, cos0 + U),+I


PBB = + / i-1
2 DY (C, cos + U)}+'


PCC = +1 V
cC = ( )+ +
2 DY2(C, cos 0 + U)


PDD= ( +
2 DY(C, cos + U)
S(kC,(l- cos2 )'
CSQ 2(C,cosO+U)].










SP = kC,(1 -cos2 ) i+
CSQP1 2(CcosB + U) j

CS = ('i- kc)cosO'

CSP1 = (ki+ k+'1)cos8'i+


The PAA, PBB, PCC, and PDD terms are used for the double y derivative as expressed by

equation (4.28) on page 54 of the dissertation.

The next three lines are
AMPIJ=AMP(I,J)
TTDD=TD(I,J) (2.3)
CCG=CG(I,J)
These three lines are used to set values to be passed in the argument of the call to the

WDIS subroutine. These values are passed in the argument rather than through common

statements because for some applications such as a shore parallel breakwater the WDIS

subroutine is called for the location at the grid edge where the depth is not an array. Also

the second call to WDIS within the J=2,N do loop is for a temporary value of the wave

amplitude that is not retained as an array.

The first call to WDIS returns the wave dissipation value W1J. This is the wi term.

Next an intermediate value of A}+1 is obtained using the explicit formulation of equation

(4.29) on page 54 of the dissertation. This intermediate value of A.+' is designated as

ATEMP in the program. Using this intermediate ATEMP value a second call to WDIS is

made returning the value of W2J, the w*+ term.

Finally the AA(J), BB(J), CC(J), and DD(J) terms are determined using all of the

short hand terms listed above.

2.2.2 The J=1 and J=N Boundaries

The J=1 and J=N values of the AA, BB, CC, and DD arrays are determined in a

manner very similar to the determination in the above do loop except that certain boundary

conditions are incorporated. First a value of ANG is determined by solving equation (4.47)

on page 60 of the dissertation for m. Equation (4.47) is finite difference on the I grid row










with known values of AMP. Thus

i(A's AX)
m = (2.4)
DY (A;+A;)
2
or in code ANG=-AII*(AMP(I,2)-AMP(I,1))/(DY*(AMP(I,2)+AMP(I,1))/2). There is of

course a safety valve to avoid division of zero so that if Ai + Ai is zero ANG is set to zero.

The same procedure is used for J=N using A's and A'N_1 to determine a value for ANG.

It is not really necessary to calculate a value for m since it would be easier to use

m = kosin0o = constant throughout. There is no advantage to computing m at each

boundary for each grid row. However that is how it currently is and perhaps if other

boundary conditions such as longshore gradients in the water surface due to tidal circulation

are added in the future, the present method may be needed.

Once m is determined then the coefficient BVAL is determined

I imn-
I 2+ (2.5)

This term is used to eliminate At using equation (4.35) on page 55 of the dissertation. For

J=N A',+ is eliminated using equation (4.36) on page 55 of the dissertation.

Other boundary conditions in the J=1 and J=N computation are that there are no

gradients of the current, depth, wave number, and group velocity terms across the lateral

boundaries. This means that values of SIG, CG, RK, and V for grids outside the compu-

tational domain (but adjacent to the boundary) are eliminated in terms of the value of the

term at the grid inside the domain (and adjacent to the boundary).

2.3 The RAD Subroutine

In the RAD subroutine the first task is to compute the WAVEHGT array which will

be used in the SHEAR2 subroutine. Then the radiation stress terms are computed using

equations (2.57) and (2.59) on page 29 of the dissertation. There is a do loop J=2,N-1 for

the first I=1 grid row which uses a forward difference derivative in x because of the off shore

boundary. Then there is the main do loop I=2,IMAX-1 with an internal do loop J=2,N-1

using central difference derivatives. Finally there is a do loop for J=2,N-1 for each of the








15

IWET(J) grids using a three point derivative incorporating the zero wave height condition

at the shore line.

The following are some of the terms that are calculated for intermediary use

TWOKH = 2kh
aB aA
DADX = -- +ik cos 9A
OB* aA*
DACONDX = i cos0A'
ax ax
aB aA
DADY =
ay ay
aB* aA*
DACONDY --=
By ay
AC = BI = AA*

Using these values the SXX, SYY, and SXY terms are then computed using equation

(2.57) on page 29 of the dissertation. The code is straight-forward and should not be difficult

to decipher using the terms above.

The do loop which calculates the radiation stress terms on the grid adjacent to the shore

line uses a three point x derivative. The three points are the zero amplitude at the shore

edge of the IWET(J) grid and the two amplitude values at the IWET(J) and IWET(J)-1

grid centers. Using a Taylor series to expand the terms about the IWET(J) grid center it

can be shown that to second order the x derivative of aA is written as
(_3AIWET(J) AWET(J)-1
(-3A AJ (2.6)
3AX

with a similar differencing for the complex conjugate.

The final do loop of the RAD subroutine defines the radiation stress terms for the J=1

grid as being equal to those at the J=2 grid, and similarly the J=N terms as being equal to

those at the J=N-1 grid. This is an extension of the no gradients in the y direction boundary

condition. This is done since the central difference y derivative can not be used at the J=1

and J=N grids. Technically this is not a correct application of the boundary condition.

Perhaps more correct would be to use the = imA boundary condition. However I don't

think that this introduces any significant error.















CHAPTER 3
THE CIRCULATION MODEL IN GREATER DETAIL



The governing equations for the circulation model are the depth integrated time aver-

aged horizontal equations of momentum and continuity given by equations (2.1), (2.2), and

(2.3) on page 16 of the dissertation. When all of the assumptions concerning friction and

mixing are included, and the numerical ADI scheme described in section 4.1 is used these

equations are finite difference as follows. In the x direction the momentum equation is

UNEWi, Ui ETANEW^j ETANEW,_-,
AT AX

+Ui, Utj U 4 + 4 (V + Vi- + Vij+ + Vi-,+) UY ,2' Uij-
1 SXXi SXXi,-i, SXYYi_ I- SXYi- 1+ SXYi-j+_ SXYi-_ij-_
+-- AX + 2AY

+p-TBX,,UNEWij EY 1 + ii- = 0 (3.1)


and the continuity equation is

ETANEWi ETAj (DUNEW);+1, (DUNEW)i, ((DV),i+i DV) (
AT AX AY

where D is the depth at the grid edge obtained by averaging the values at the two appropriate

grid centers.

For the y direction the momentum equation is difference as

VNEWi. Vi,j ETANEWij ETANEWij-1
AT AY
I Vi+1, Vi-j V+lj+ Vi'-1
+4 (U,j + U-lj + uj+l + Ui-1j+1) 2AX J 2AY
1 (SXY.+i SXYi-_j + SXY+.i_,-1 SXY-.ij-1 SYYi, SYYi.i-1
+- .
'pD 2AX S Ay
V1 V + 1,- 2V' V + V = 0
+ TBYijVNEWi+ + EKX., + ,j + .i =0 (3.3)










and the continuity equation as

ETANEWi, ETANEW7j (DVNEW).ij+i (bVNEW);j (DV)i,i+i (DV)ii
AT + AY AY
(3.4)

where EXij is the averaged EX term. Note that the (Dv)'y i term is subtracted

from the right hand side of (3.2) and added to the right hand side of (3.4). Thus if the two

equations (3.2) and (3.4) are added, a fully implicit two dimensional continuity equation

results.

The basic subroutines in the circulation model are XCOEF and YCOEF which deter-

mine the coefficient vectors of the tridiagonal matrix.

3.1 The XCOEF Subroutine

If equations (3.1) and (3.2) are written as


A2i-2ETANEW*(I-1,J) + B2;-2UNEW(I,J) + C2i-2ETANEW'(!,J) = D2i-


A2.-1UNEW(I,J) + B2i-1ETANEW'(I,J) + C2j-iUNEW(I+1,J) = D2i-1


(3.5)


(3.6)


then


A(2*I-2) =

B(2*I-2) =

C(2*I-2) =

D(2*I-2) =








A(2*I-1) =

B(2*I-1) =


-g* RX

1 + AT TBX(I, J)/(pD)

g* RX
Ui+1,i Ui-li
U(I, J) AT vij.i -1

- -(V, X + Vij + Vij+ + Vi-lj+l) Ui+ i-1
AT (SXXijjj+1 X,-i SXYi,+i SXYij-j + SXY.-ij+1 SXYi-ij_-i
pD AX 2AY
+AT Ui+ 2Ui + Uii-
Ay2

- RX D,

1.0










C(2*I-1) = RX *Di+,,i

D(2*I-1) = ETAij RY *(Vi+1~i,-+ VjjD,) (3.7)

where RX = A and RY = AT

The above are the basis for all of the subroutines that have an X as the first letter in the

name. Variations of the above are found near the boundaries or next to internal boundaries

due to the presence of structures. Some examples of this follow.

In the XCOEF subroutine the offshore boundary condition of i7 = 0 is given in the code

as
A(1) = 0.0
B(1) = 1.0
C(1) = 0.0
D(1) = 0.0

When (3.7) is written for 1=2 there is the problem that Ui-,1 is not defined since it

is outside of the computational domain which goes offshore only as for as the center of the

first grid row. The Ui-Ij term, which is used in the advective acceleration terms, is avoided

by using the continuity equation to solve for a.

aU +1 [ U auD av*D (3.8)
z D5 at -z ay

The right hand side of the above equation is thus used for the x derivative of U in the

computation of the D(2) term.

The shore boundary condition.of no flow into the shore is represented as UIWET(J)+,J =

0. In code this is
A(2*IWET(J)) = 0.0
B(2*IWET(J)) = 1.0
C(2*IWET(J)) = 0.0
D(2*IWET(J)) = 0.0

Once the XCOEF subroutine is understood the other X-series subroutines are easy to

decipher. Each one is merely an adaptation to make sure that y derivatives and depth

averages are not taken across a barrier or outside of the computational domain. XCOJ1 is

written for the special case of J=1. XCOJN is for the case of J=N. In the [WINER.GROIN]








19

director there are also the special subroutines XCJGR and XCJGR1 which are written

for the case of a shore perpendicular barrier between the J=JGR and J=JGR+1 grid

columns. These two subroutines are for the two grid columns adjacent to the barrier.

In the [WINER.CLOSED] directory there is a subroutine XCOJ33 which is in the GPRO

program which simulates the lab experiment of Gourlay as reported in the dissertation.

This subroutine is written since the J=33 grid column is next to a boundary for the first 30

grids. Also in the IWINER.CLOSED] directory there is the XCO2 subroutine written for

J=34,N since computations are only for the 1=31 through I=IWET(J) grids.

3.2 The YCOEF Subroutine

If equations (3.3) and (3.4) are written as


A2i-1ETANEW(I,J-1) + B2i-1VNEW(I,J) + C2i-IETANEW(I,J) = D2,i-


A2,VNEW(I,J) + B2y ETANEW(I,J) + C2iVNEW(I,J+1) = D2i


(3.9)


(3.10)


then


A(2*J-1) =

B(2*J-1) =

C(2*J-1) =

D(2*J-1) =








A(2*J) =

B(2*J) =

C(2*J) =

D(2*J) =


-g* RY

1 + AT TBY(I, J)/(pD)

g*RY

V(I, J) AT Vi vi,+l 2ViY-1
2AY
S(Ui, + Ui+l, + Ui,-1 + Ui+1.-1) +l Vj
AT (SYYij SYY,i-I SXYSXY+-XYi-l, + SXYi+li-1 SXYi-1,y-1
pD AY 2AX
+AT EX Vi+1,i 2V1j" + Vi-li
AX2

-RY *Di

1.0

RY Dij+

ETANEWI* + RY (Vij+1Dij+1 VjDiaj) (3.11)








20

The above is the essence of the YCOEF subroutine. For J=l and for J=N+1 the D(2*J-

1) term is coded so as to incorporate the lateral boundary condition of no y gradients of the

circulation terms at the boundary. This means that the y derivative terms are not included

for D(1) and D(2*N+1).

As with the X-series of subroutines so with the Y-series of subroutines. Any subroutine

with a Y as the first letter is a variation of YCOEF. The YCVAR subroutine, which is

named for the variable shoreline position is the same as the YCOEF subroutine except

that it has many IF statements to determine the position of the shoreline (computational

domain boundary) in relation to the present grid. For each grid, if the boundary is adjacent

to the grid in the x direction then the x derivative terms for the radiation stress SXY, the

mixing and the advective acceleration are all written using a three point derivative that

incorporates the zero longshore velocity (V=0) and zero wave height (SXY=0) at the shore

grid edge.

In the [WINER.BARR] directory YCOJJJ and YCJJJ1 are written for the grid row

on the offshore and onshore sides respectively of a barrier that is on the shore edge of the

I=JJJ grid row.

In the (WINER.CLOSED] directory which is written for a closed basin the V=0 lateral

boundary condition is incorporated by having D(1)=0 and D(2*N+1)=0. YCOI1 is used

for the I=1 grid row and incorporates the U(1,J)=0 offshore boundary condition (which is

used for a closed lab basin). YCI31 is written for the GPRO program to incorporate the

boundary of the basin from J=34 through J=N.

3.3 The CIRC Subroutine

The CIRC subroutine varies in each of the directories. The variation is but to include

the calls to the necessary special subroutines when appropriate. CIRC sequentially calls

the X- series subroutines. After each call to a X-series subroutine a call is made to TRIDA

to solve the tridiagonal matrix equation and then the solution vector VV is assigned to the

appropriate variables. After the X-series of calls the CIRC subroutines sequentially calls








21

the Y-series of subroutines making a call to TRIDA to solve the tridiagonal matrix equation

and the solution vector VV is assigned to the appropriate variables.

It should be noted that the sequential calls to the X-series and then to the Y-series

subroutines is done only because of the linear nature of the VAX computer. The calls to the

XCOEF could be vectorized if a computer with parallel processors is used. Likewise the calls

to the Y-series could be vectorized. However the X-series of calls must be completed before

any calls to the Y-series is made since the Y-series uses ETANEW* which is determined in

the x sweep.

3.4 The FLOOD subroutine

The FLOOD subroutine is called directly from the main program after the call to UP.

This discussion of the FLOOD subroutine is included in the chapter on the circulation model

since the FLOOD subroutine controls the shoreline boundary of the circulation model.

The FLOOD subroutine is quite simple in concept but a little cumbersome in its code.

The subroutine goes grid column by grid column from J=1 through J=N to determine

whether the shore boundary should be moved. First it determines whether the last grid in

the grid column has gone dry. This is simply whether the total depth TD(IWET(J),J) is

less than zero. If it is IWET(J), the number of grids in the J grid column, is reduced by

one. Further the amount of the water deficiency in the IWET(J+1),J grid (which used to

be the IWET(J),J grid) is taken from the IWET(J),J grid and added to the IWET(J+1),J

grid so that the total depth in the IWET(J+1),J grid is zero. This is done so as to maintain

conservation of water mass in the computational domain. The above requires corrections

to both the ETA and TD values.

If the total depth in the IWET(J),J grid is not less than zero the subroutine checks to

see if the water level in the IWET(J),J grid is greater than the elevation of the neighboring

IWET(J+1),J grid. This is a simple IF statement

IF((ETA(IWET(J),J)+H(IWET(J),J)).GT.O.O)THEN

which determines this since ETA and H are measured positive in opposite directions. If








22

the above IF statement is true then IWET(J) is incremented by one and an arbitrary one

millimeter (.001 meters) of water is taken from the IWET(J-1),J grid and placed in the

newly wetted IWET(J),J grid. Again this is done for both the ETA and TD values.

At present the flooding mechanism is only in the x direction. This requires a lot of if

statements. To allow for flooding in the y direction would require many more IF statements

and it is doubtful whether this would be worth while since with a parabolic wave model

wave energy is essentially propagated in only one direction. If other wave models are to be

included in the model, it may be worth implementing a two directional flooding mechanism.
















CHAPTER 4
VARIATIONS OF THE PROGRAM



All of the above described subroutines and the executable program can be found in

the IWINER.MODEL] directory (or else have been put on tape in which case see Sid-

ney or Suburna for their location). The program in the [WINER.MODELJ directory is

for cases with open lateral boundaries and no internal boundaries due to the presence of

structures. Variations of the main program in [WINER.MODEL) have been written for

the following cases: for closed (reflective) lateral boundaries, for shore parallel breakwa-

ters and for shore perpendicular breakwaters. The programs for closed lateral bound-

aries are in the [WINER.CLOSED) directory. Those for a shore parallel breakwater are

in the [WINER.BARR] directory and those for a shore perpendicular breakwater are in

[WINER.GROIN]. This chapter describes the adaptations made in each of these directo-

ries. The last section of this chapter also describes some of the plotting programs found in

the [WINER.PLOT] directory.

4.1 The [WINER.BARRJ Directory

In this directory the program is modified to include the internal boundary conditions

found with a shore parallel breakwater that is represented as an infinitesimally thin imper-

meable barrier. The structure of the program is essentially the same as described in the

preceding chapters. Many of the subroutines are exactly the same. Those that are changed

are described below. The main program is named BARR.FOR and the executable file is

BARR.EXE. The option file is BARR.OPT and has a listing of all the subroutines used in

the BARR program. The linking command is

$LINK BARR/OPT








24

In the BARR.OPT file there is listed the new subroutines ACOH1, ACOH2, WVNMHA,

YCOJJJ, and YCJJJ1. WVRAD is a modification of WAVRAD, CIRBARR is a modifi-

cation of CIRC and RADBARR is a modification of RAD. With the exception of slight

changes to XCOEF and INPUT and the main program, the other files are unchanged.

The DIM.FOR file that is included in all the subroutines remains unchanged except for

the addition of several arrays that are passed through the following common statements.
COMMON/Q/RKHA(100),RKOHA,CGHA(100),SIGHA(100)
COMMON/R/REFLEC(100,100),AMPSAV(100),AMPH(100),REFL(100)

SIGHA(100), CGHA(100), and RKHA(100) are the relative frequency, the group velocity,

the wave number arrays at the position half way between the I=JJJ and I=JJJ+1 grid

centers, which is where the breakwater is located. RKOHA is the averaged wave number

at the shore side of the I=JJJ grid row. These values are computed in the WVNMHA

subroutine. REFLEC, AMPSAV, AMPH, and REFL are declared complex. AMPH and

AMPSAV are the complex amplitude at the grid edge position of the breakwater. These

are used in the WVRAD and RADBARR subroutines. REFLEC and REFL are not used

in the program. They remain in the common statement as a remnant of an unsuccessful

attempt to add reflection from the breakwater to the model.

The INPUT subroutine is modified so as to read and then pass the position of the

breakwater. The input data file has two additional values that are read by the INPUT

subroutine. These are JJJ and NBB which are passed in the argument of the INPUT

subroutine back to the main program. JJJ is the position of the breakwater in the x

direction. The breakwater is on the shore side of the I=JJJ grid row, between the I=JJJ and

I=JJJ+1 grid rows. NBB is the number of grids for one half the length of the breakwater.

Thus the breakwater length is equal to 2*NBB*DY. The rest of the INPUT subroutine is

unchanged from that described previously. Within the [WINER.BARR] directory there is

a sample run available. The input file IN.BARR and depth file PL.DEP is a test case of

normal incident waves with a breakwater 100 meters offshore (JJJ=10, DX=10m for a plane

beach with 20 wet grids in the x direction) and 100 meters in length (NBB=5, DY=10).








25

Once NBB is passed back to the main it is used to determine the fixed end positions

of the breakwater N1 and N2. Thus the breakwater is on the shore side of the I=JJJ grid

row and is from the J=N1 through the J=N2 grids. JJJ, N!, and N2 are then passed in

the arguments of the calls to various subroutines. From the main these are passed in the

WVRAD and the CIRBARR subroutines.

The WVRAD subroutine is an adaptation of the WAVRAD subroutine described above

with the modification to include the barrier. This subroutine makes calls to WAyNUM

which is unchanged and WVNMHA which determines the wave characteristics (k, C,a, k,

and 0) at the position of the edge of the grid row where the barrier is located. The WVN-

MHA subroutine does this by taking the depth at the grid edge to be the average of the

two adjacent grids. The WVNMHA subroutine has the same structure and logic as the

WAVNUM subroutine. Next the WVRAD subroutine initializes the values of KBREAK(J)

and AMP(1,J) the same as was done in the WAVRAD subroutine described previously.

Next is a DO loop for I=1,JJJ-1 which calls AMPCO then CTRIDA and then assigns the

solution vector to the AMP values on the next grid row. This gives the values of the AMP

array up to and including the I=JJJ grid row.

Then a call is made to ACOH1 which has the same structure as AMPCO except that

DX is set to one half DX. ACOH1 determines the coefficients of the tridiagonal matrix

equation in terms of the known AMP values and the wave characteristics (k, k, a, Cg) on the

I=JJJ grid row and the wave characteristics on the grid edge (the RKHA, CGHA, SIGHA,

and RKOHA values). Then a call is made to CTRIDA so that the solution vector is the

complex amplitude at the grid edge half way between the I=JJJ and I=JJJ+1 grid centers.

The solution vector VVV is then assigned to the AMPH array. Then for J=N1,N2 the

value of AMPH(J) is assigned to AMPSAV(J) and then AMPH(J) is set to zero. This is

done so that AMPSAV(J) can be used in the RADBARR subroutine to get a three point x

derivative on the offshore side of the breakwater and so that AMPH(J) can be used in the

ACOH2 subroutine to advance the parabolic model to the I=JJJ+1 grid center.








26

Next the WVRAD subroutine calls the ACOH2 subroutine. The ACOH2 subroutine

has the same structure as the AMPCO subroutine. The difference is that DX is again

one half DX and the coefficients are in terms of the known AMPH values and the wave

characteristics on the grid edge (the RKHA, CGHA, SIGHA, and RKOHA values) and the

wave characteristics (k, k, a, C.) on the I=JJJ+1 grid row. A call is then made to CTRIDA

and the solution vector is assigned to the AMP values on the I=JJJ+1 grid row. With the

values of AMP known on the I=JJJ+1 grid row the rest of the AMP array is determined

by a do loop for I=JJJ+1,IMAX-1 which calls AMPCO then CTRIDA and then assigns

the solution vector to the AMP values on the next grid row.

The final execution in the WVRAD subroutine is a call to RADBARR which is the

same as the RAD subroutine previously described except that there is the additional do

loops for J=N1,N2 for both I=JJJ and I=JJJ+1. These do loops have the same structure

as other computations in the RAD subroutines except that DADX and DACONDX are

written as three point derivatives. For the I=JJJ grid row the three point derivatives are

written using the AMPSAV array and on the I=JJJ+1 grid row the derivative are written

using a zero value at the lee side of the breakwater.

The CIRBARR subroutine differs from the CIRC subroutine previously described only

in that for the y sweep calls to YCOEF are made for J=2,JJJ-1 and then a call is made to

YCOJJJ and then YCJJJ1 which are adaptations of YCOEF for the I=JJJ and I=JJJ+1

grid rows respectively. These two subroutines are just like YCOEF except that for values

of J=N1,N2 changes are made so that x derivatives (V, -i, "r) are not taken across the

boundary. For the derivatives involving V the assumption is made that the barrier has a

free-slip condition on the V velocity which enables the use of the !' = 0 condition at the

barrier to eliminate the value on the other side of the barrier from the derivative. For the

derivatives of S., a reflective = 0 condition is used on the offshore side of the barrier

and on the lee side the zero wave condition (S=y = 0) at the barrier is used in a three point

derivative.








27

The only other change in the circulation model is a modification to the XCOEF sub-

routine. JJJ, N1, and N2 are passed in the argument of the call to XCOEF. The subroutine

is the same as previously described except for two if statements at the end to determine

if J is equal to or greater than N1 and if it is also equal to or less than N2. If it is, then

A(2*JJJ), C(2*JJJ), and D(2*JJJ) are set to zero and B(2*JJJ) is set to one. These are

the tridiagonal matrix coefficients for the equation U(JJJ+1,J)=0.

4.2 The JWINER.GROIN] Directory

This directory contains a modification of the program found in the (WINER.MODEL]

directory. The modification is for the case of a shore perpendicular breakwater that is rep-

resented as an infinitesimally thin impermeable barrier located between the J=JGR and

J=JGR+1 grid columns. Most of the subroutines are the same as in the [WINER.MODELJ

directory or with minor modifications. The main program is GROIN.FOR and the exe-

cutable file is also GROIN. A list of all the subroutines is found in GROIN.OPT and the

link command is

$LINK GROIN/OPT

The only new subroutines are AMPLL and AMPLLL which are modifications of AM-

PCO to include the internal boundary conditions at the shore perpendicular breakwater,

and XCJGR and XCJGR1 which are modifications of the XCOEF subroutine to incorporate

the boundary conditions adjacent to the breakwater.

The INPUT subroutine is modified to read the position and length of the groin (groin

and breakwater and jetty are used interchangeably in this section) and then pass this info

to the main through the argument of the call. Thus the argument of the INPUT subroutine

has two additional variables JGR and LLL. JGR is the J position of the breakwater. It is

located between the J=JGR and the J=JGR+1 grid columns. JGR is read directly from

the input data file. LLL is the grid row position of the end of the jetty. The input data file

has two more variables in this directory than in the [WINER.MODELI directory. The first

is JGR and the second is the length of the groin GRLEN. The length of the groin should be








28

in increments of DX. LLL is then determined after the IWET array is determined by the

definition

LLL=IWET(JGR)-INT(GRLEN/DX)+1

the LLL and JGR variables are then sent in the arguments of the CIRC and WAVRAD

calls so that the internal boundary conditions can be met.

The WAVRAD subroutine has the values of LLL and JGR passed in the argument of

the call statement. This subroutine has the same structure as the subroutine with the same

name in the [WINER.MODEL] directory. The WAVNUM, AMPCO, CTRIDA, and WDIS

subroutines are unchanged. The first do loop to call AMPCO in the WAVRAD subroutine

is for I=1,LLL-2. This moves the solution for the AMP array up to the LLL-1 grid row

which is the grid row just offshore of the end of the groin. Two adaptations of the AMPCO

subroutine are then used to solve for the rest of the domain. AMPLL which is specifically for

the I=LLL-1 call which determines the tridiagonal matrix coefficients for the I=LLL-1 grid

row which does not have an internal boundary condition, and the I=LLL grid row which

does have an internal boundary condition. The other adaptation of the AMPCO subroutine

is AMPLLL which is called sequentially for I=LLL,IMAX-1 to complete the solution for

the AMP array.

The internal boundary condition is handled in the following manner. When the parabolic

equation (equation 4.27 on page 53 of the dissertation) is written for J=JGR, there are terms

for the J=JGR+1 grid on the other side of the barrier. The barrier is treated as perfectly

reflective. Thus any quantity on the other side of the barrier is eliminated from the equa-

tion in terms of the quantity on the proper side of the barrier. Similarly for the J=JGR+I

grid, values at the J=JGR grid are eliminated in terms of the value at the J=JGR+1 grid.

Essentially the barrier at the groin is treated the same as the reflective boundary condition

as discussed in section 4.4 of the dissertation.

In the RAD subroutine the values of LLL and JGR are passed in the argument of the

call. In RAD a special do loop is written for the two grid columns on each side of the groin to








29
insure that the y derivatives are not taken across the barrier. Thus for I=LLL,IWET(JGR)

on the J=JGR grid column the y derivatives of AMP (and similarly for the complex conju-

gate) are written as

(AMP(I,J)-AMP(I,J-1))/(2*DY)

and for the J=JGR+1 grid column they are written as

(AMP(I,J+1)-AMP(I,J))/(2*DY)

These two expressions result from using the reflective condition to eliminate the value of

AMP on the other side of the barrier in terms of the value of AMP on the proper side of

the barrier.

The CIRC subroutine is altered in that LLL and JGR are passed in the argument of

the call, and that the two subroutines XCJGR and XCJGR1 are added to determine the

coefficients for the x sweep for the special cases of J=JGR and J=JGR+1 respectively.

The XCJGR and XCJGR1 subroutines are adaptations of XCOEF so as to avoid taking

derivatives across the impermeable barrier. The YCOEF and YCVAR subroutines are also

altered in that LLL and JGR are passed in the argument. If I is equal to or greater than

LLL the internal boundary condition V(I,JGR+1)=0 is imposed by setting the appropriate

coefficients to zero in the same manner as was done for the XCOEF subroutine in the

previous section.

4.3 The [WINER.CLOSED] Directory

In this directory there are two executable programs PRO and GPRO. This directory

has adaptations of the program for laboratory basin studies where there lateral boundary

conditions are zero velocity across the boundaries and also a similar condition on the offshore

boundary. The discussion here will be brief since the reader by now should understand how

adaptations are made to the program.

For the PRO program the changes are minimal, the AMPCO subroutine is changed to

have a reflective condition at J=1 and J=N. This is easily accomplished by setting BVAL








30

equal to one. The YCOEF and YCVAR subroutines are similarly changed to incorporate

the V=0 condition for J=1 and J=N+1. The CIRC subroutine also has a call to a new

subroutine YCOI1 which is for the y sweep for the special case of I=1 which must be

computed since the ETA values on the offshore grid row are not constrained to being zero.

The GPRO program is a special adaptation used to simulate the experiments of Gourlay

as described in section 5.6 of the dissertation. By looking at figure 5.30 on page 98 of the

dissertation it is seen that the computational domain is somewhat complicated. The CIRC

subroutine thus has special calls to XCOJ33 for the special case of J=33, to XCO2 which is

an adaptation of XCOEF for the cases of J=34,N since the computation begins at 1=31, and

to XCOJN2 which differs from XCOJN in that the computations begin at 1=31. Also there

is a call to YCOI31 for the special case of 1=31 which includes boundary conditions. For

the XC series of subroutines the indices of the coefficients is shifted since the computations

begin with U(1,J) rather than ETA(1,J) as in the other directories.

IN.PL is a data file and PL.DEP is a depth file for a case of a plane beach at angle

to the wave maker in a rectangular basin. GOUR.IN and GOUR.DEP are the input and

depth files for the Gourlay experiment.

4.4 The [WINER.PLOT] Directory

In this directory there are several files which use the NCAR and QTEK plotting rou-

tines. If the reader is familiar with these packages there should be no problem using the

plotting programs in this directory. The NCAR programs are only as fortran files. This is

because in the NCAR package the dimension statement must be the exact same as the di-

mension of the array that is sent in the argument of the call to the NCAR subroutine. This

is remedied by having a parameter statement as the first statement of the NCAR fortran

files. These files then must by edited to change the parameters M and N to the appropriate

size, compiled and then linked. The linking command is simply

SLINK filnam










$@[KIRBY]NCAR


command has already been executed. The user could put such a command in his LO-

GIN.COM file in the main directory so as to automatically execute it. This command

makes the NCAR library a part of the system library as far as linking command are con-

cerned.

The NCAR fortran files in the [WINER.PLOT] directory are ETA, SURFACE, VEL,

STREAM, and QVEC. Each of these programs reads an output file in the same format as it

was written by the WRITE subroutine. The ETA program reads only the ETA array and

makes a call to the NCAR subroutine EZCNTR to produce a contour plot of the surface

elevation. The ETA program can also be used to read the AMP.LST and ETA.LST files to

get contours of the wave height and of the instantaneous water surface due to the wave field.

The SURFACE file makes a call to SRFACE. This file was obtained from Rengi Phillip so

I cannot describe the logic. DX and DY should be changed as needed along with the M

and N in the parameter statement. Otherwise the program prompts for the data file (which

must be entered with a single quote on each side) and the eye position of the viewer, etc.

The VEL program reads the output file to get the ETA U and V arrays, averages the U and

V from two opposite grid edges to get grid centered velocity arrays (called UC and VC) and

calls the NCAR subroutine EZVEC to a plot of the velocity vectors. The QVEC program

is a minor change to the VEL program. The QVEC program also reads a depth file to get

a depth array and then creates the total depth array by adding the ETA array to the depth

array. Then the UC and VC arrays are multiplied by the total depth array so as to create

arrays of the x and y components of the total discharge. A call to EZVEC then gives a plot

of the vectors of the total discharge.
















APPENDIX A
VARIABLES USED IN THE COMPUTER CODE



A is a column vector of coefficients for-the tridiagonal matrix equation resulting from

the governing circulation equations. It is dimensioned (200).
SV \ti+1
A2 in the AMPCO subroutine is co, +u)l

AA is a column vector of coefficients for the tridiagonal matrix equation resulting from

the parabolic equation. It is dimensioned (100).

AB in the AMPOUT subroutine is the absolute value of the amplitude.

AII is an imaginary number, the square root of minus one

AMP is the complex amplitude. It is dimensioned (100,100).

AMPIJ is used in the AMPCO subroutine and is the value of Ai,j which is passed in

the argument of the WDIS call.

AMPO is the complex amplitude on the first grid row. It is dimensioned (100)

ANG is used in the AMPCO subroutines and is the m or longshore component of the

wavenumber.

ATEMP is a complex variable used in the AMPCO series subroutines. It is actually

the value of the complex amplitude as computed by an explicit finite differencing of the

parabolic wave equation. The value of ATEMP is then passed to the WDIS subroutine in

order to compute the amount of wave energy dissipation to be used. It is computed using

equation (4.29) on page 54 of the dissertation.

ATMP is used in the AMPCO subroutines and is imrn

B is a column vector of coefficients for the tridiagonal matrix equation resulting from

the governing circulation equations. It is dimensioned (200).








33

Cg co.i+s 1 #+u CGco .e+U
B2 in the AMPCO subroutine is ---'---> v--)-cau
( + +(C'
BB is a column vector of coefficients for the tridiagonal matrix equation resulting from

the parabolic equation. It is dimensioned (100).

BETA is a variable that is used both the TRIDA and the CTRIDA subroutines. It is

complex in CTRIDA and not complex in TRIDA. It is not passed by any common statements

and is used only in the subroutines for the double sweep solution. In CTRIDA it must be

dimensioned at least as large as N. In TRIDA it must be at least 2N+1 or 2M+1, which

ever is larger.

BVAL is used at the lateral boundaries in the AMPCO subroutines.It is the coefficient

of equation (4.35) or (4.36) on page 55 of the dissertation.

C is a column vector of coefficients for the tridiagonal matrix equation resulting from

the governing circulation equations. It is dimensioned (200).

C2 in the AMPCO subroutine is cose+u) DY'

CC is a column vector of coefficients for the tridiagonal matrix equation resulting from

the parabolic equation. It is dimensioned (100).

CG is the group velocity and is dimensioned (100,100). The values of the CG array are

determined in the WAVNUM subroutine using the relationship CG= n times the celerity,

where
1 1 2RK TD2
2 ( sinh(2RK TD2)(A1)
CGPU in the AMPCO routine is (C, + U),ij.

CGP1U in the AMPCO routine is (C, + U)i+1,,

CNN is a constant used in the EXCOEF subroutine it is the N of equation (2.24) on

page 23 of the dissertation. It has arbitrarily been set to .08.

CS in the AMPCO subroutine is (ki k}) cos 8'

CSP1 in the AMPCO subroutine is (kci+ k, +) cos 8i+1

CSQ in the AMPCO subroutine is (c,.+)'








34

CSQP1 in the AMPCO subroutine is ( Cg -co, )+1

CXP in the AMPOUT subroutine is ef kd

D is a column vector of the right hand side for the tridiagonal matrix equation resulting

from the governing circulation equations. It is dimensioned (200).

D2 in the AMPCO subroutine is (c, co +U) i )

DACONDX in the RAD subroutine is the x derivative of the complex conjugate of

Aei Idfz

DACONDY in the RAD subroutine is the y derivative of the complex conjugate of

Aei f ;dz

DADX in the RAD subroutine is the x derivative of AeiJf dz

DADY in the RAD subroutine is the y derivative of Aei f idz

DD is a column vector of right hand side for the tridiagonal matrix equation resulting

from the parabolic equation. It is dimensioned (100).

DT is the time step increment used in the circulation model.

DX is the grid spacing in the x direction.

DY is the grid spacing in the y direction.

EMAX is the maximum value that the eddy viscosity coefficients may be. EMAX is

used so as to place a cap on the eddy viscosity coefficients and thus avoid numerical stability

problems that may arise from too large a coefficient. This is important since the lateral

mixing terms are represented explicitly in the model and thus place a constraint on the size

of either the delta t increment or the lateral mixing coefficient.

EPS and EPSS are arbitrary e used in the DISP, CHECK, and WAVNUM subroutines.

ETA is the mean surface elevation measured upwards from the still water line. It is

dimensioned (100,100).

ETABAR in the AMPOUT subroutine is the value of the instantaneous water surface

due to the waves obtained from equation (3.72) on page 42 of the dissertation.

ETANEW is the (100,100) array of the newly calculated mean water surface elevation.








35
ETADT is the time rate of change of ETA. It is dimensioned (100,100), but only values

for the first two offshore grid rows are used. The values of ETADT are determined in the

UP subroutine prior to reassigning the values of ETA. The values of ETADT are used in

the XCOEF series of subroutines to determine the non- linear term. This is necessary

since U(1,J) is not calculated and thus a centered x derivative at the U(2,J) position is not

directly obtainable. Instead the continuity equation is solved for a.

(U 1 / OTD OTDV\
OU ETADT + U aOx ) (A.2)
5x TD ( ay

EX is an eddy viscosity coefficient that is dimensioned (100,100). Its values are cal-

culated in the subroutine EXCOEF according to a Longuet-Higgens formulation given by

equation (2.24) on page 23 of the dissertation.

EY is an eddy viscosity coefficient which is equal to EMAX.

FILNAM is used in the INPUT subroutine. It is declared as LOGICAL*1 FILNAM(30)

and is used to enter the name of the input data file and the depth data file.

FRIC is a constant friction coefficient that is an input variable.

G is the gravitational acceleration constant.

GAMMA is a variable that is used both the TRIDA and the CTRIDA subroutines.

It is complex in CTRIDA and not complex in TRIDA. It is not passed by any common

statements and is used only in the subroutines for the double sweep solution. In CTRIDA

it must be dimensioned at least as large as N. In TRIDA it must be at least 2N+1 or 2M+1,

which ever is larger.

H is the water depth measured positive down from the still water line. It is dimensioned

(100,100)

IMAX is the largest value in the IWET array. The computational domain of the wave

model is IMAX by N. Also the output is in the form of IMAX by N arrays.

IMIN is the smallest value in the IWET array.

ITER is a count of the number of interactions.

ITERST is a constant that is used in the main to determine the smoothly increasing








36

values of the wave height during the start-up of the model. It is the same as Cl in equation

(4.45) on page 60 of the dissertation except that it is expressed in units of iterations.

IWET, which is dimensioned (100), is the number of grids (for each grid column) from

the offshore boundary to the shoreline. The values of the IWET array are changed by the

FLOOD subroutine with the rise or fall of the mean water level.

KBREAK is an index which is used in the WDIS subroutine to indicate whether or

not dissipation due to breaking should be applied. If KBREAK(J) is zero then the WDIS

subroutine returns a zero dissipation coefficient. If KBREAK(J) is one, a positive dissipation

coefficient is returned. The WDIS subroutine controls the value of KBREAK(J). If the wave

height becomes greater than .78 times the depth KBREAK(J) is set to one and remains so

until the wave height drops below a value of .4 times the local depth. The rational for this

is in Daily's dissertation cited in my dissertation.

KFLAG is an indicator that is controlled by the CHECK subroutine to indicate con-

vergence. KFLAG is set to zero at the beginning of the CHECK subroutine and then any

non-converge will set KFLAG to 2 and return to the main.

M is the number of grids in the x direction. It should be noted that M must be

sufficiently large to allow for flooding of the shoreline. The actual computations will take

place over a number of grids somewhat less than M.

N is the number of grids in the y direction.

OMEGA is the wave frequency = 2*PI/T where PI is r=3.14159

PAA in the AMPCO subroutine is used as part of the differencing of the (CCGA)


term. It is 2 +Y I
2 DY2(C, coa +U)j.+
PBB in the AMPCO subroutine is used as part of the differencing of the (CCGA)


term. It is
2DY2(C, cos +U).
PCC in the AMPCO subroutine is used as part of the differencing of the (CCo )


term. It is N
2 DY2(C, cos+U)








37

PDD in the AMPCO subroutine is used as part of the differencing of the (CCa)V


term. It is i DYi-
2DY2(C cos +U)')
RK is the wavenumber array. It is dimensioned (100,100) and is determined in the

WAVNUM subroutine using a Newton method solution to equation (3.76) on page 44 of the

dissertation.

w = a + k cos OU (A.3)

where a is the relative frequency given by


a2 = gk tanh kh (A.4)


RKO is the average of the wave number across a grid row. It is dimensioned (100).

RKN is the wave number along the first grid row. It is used in the INPUT subroutine.

RO is the fluid density. It has arbitrarily been set to 1000 kilograms per cubic meter.

The value of the fluid density is not important since in the places it is used it is cancelled

out by use in both the denominator and the numerator.

RX is used in all of the coefficient subroutines that are called from CIRC. It is =DT/DX

RY is used in all of the coefficient subroutines that are called from CIRC. It is =DT/DY

SIG is the relative wave frequency and is dimensioned (100,100). It is determined in

the WAVNUM subroutine using a Newton method solution to the dispersion relationship

in the absence of currents.

SUM in the AMPOUT subroutine is f kdz

SXX is the S., radiation stress term. It is dimensioned (100,100) and is determined in

the RAD subroutine.

SXY is the Sz, radiation stress term. It is dimensioned (100,100) and is determined in

the RAD subroutine.

SYY is the Syy radiation stress term. It is dimensioned (100,100) and is determined in

the RAD subroutine.

T is the wave period which is an input variable.








38
TBX is the bottom shear stress in the x direction. It has a dimension of (100,100). It

is named for Tau Bottom in the X direction.

TBY is the bottom shear stress in the y direction. It has a dimension of (100,100). It

is named for Tau Bottom in the Y direction.

TD is the total depth of the water column. It is dimensioned (1000,100) and is equal to

H plus ETA. The values of the TD array are updated each iteration in the UP subroutine.

TD2 is also the total depth of the water column and is dimensioned (100,100). It is

equal to TD except for the grids where TD is equal to or less than zero (TD should never

be less than zero, but just in case by error it is, the .LE. if statement is used). If TD is

zero or negative TD2 is set to an arbitrarily small value of 1 millimeter. The TD2 array is

updated in the WAVNUM subroutine and is used in all of the wave model subroutines.

THETA is the wave angle at each grid row assuming a plane beach without currents as

calculated by Snell's law. It is dimensioned (100).

THETAH is the wave angle calculated by Snell's law for the position half way between

grid centers ( i.e. at the grid edge) at the location of the breakwater. THETAH is used

only in the BREAK program and in the GPRO program.

THETAO is the input value of the wave angle at the offshore grid row location.

TIME is the time in minutes for the start-up of the model. It is actually the C1 constant

of equation (4.45) on page 60 of the dissertation.

U is the (100,100) array of x component of the depth integrated mean current.

UIJ is a grid centered Ui, used in the AMPCO subroutine.

UIP1J is a grid centered Ui+,j used in the AMPCO subroutine.

UNEW is the (100,100) array of the newly calculated values of the x component of the

depth integrated mean current.

UU is used in the main program and is actually the C(t) of equation (4.45) on page 60

of the dissertation. It is used to gradually increase the wave height so as to avoid a shock

start- up.









39

V is the (100,100) array of y component of the depth integrated mean current.

VIJ is a grid centered Vi,j used in the AMPCO subroutine.

VIIM1 is a grid centered Vej- used in the AMPCO subroutine.

VIJP1 is a grid centered Vijj used in the AMPCO subroutine.

VIP1J is a grid centered V,i.j used in the AMPCO subroutine.

VIP1JM1 is a grid centered V+i,j-1 used in the AMPCO subroutine.

VIP1JP1 is a grid centered V+,,j+l used in the AMPCO subroutine.

VNEW is the (100,100) array of the newly calculated values of the y component of the

depth integrated mean current.

VV is the solution vector for the circulation model. It is determined by the TRIDA

subroutine and in the CIRC subroutine (which called TRIDA) the VV vector is assigned

to the circulation values of UNEW or VNEW and ETANEW. It should be noted that VV

has a dimension that is twice that of the variables U, V, and ETA. This is because the

tridiagonal matrix equation, which is solving for a staggered grid system, solves for both U

and ETA in the x direction and for V and ETA in the y direction.

WAVAMP is the amplitude of the incident wave on the offshore grid row. It is equal

to one half of WAVGHT which is the input wave height. WAVAMP is set in INPUT and

passed through the call statement to the main program.

WAVGHT is the input wave height.

WAVEHGT, dimensioned (100,100), is the wave height which is actually twice the ab-

solute value of the complex amplitude AMP. The WAVEHGT array is used in the SHEAR2

subroutine to determine the bottom shear stresses which are designated by the array names

TBX and TBY.

WVAMP is a variable that it is passed in the argument of the WAVRAD subroutine. It is

some value from 0.0-1.0 times WAVAMP. The value 0.0-1.0 is actually C(t) as determined by

equation (4.45) on page 60 of the dissertation. WVAMP times AMPO(J) is used to initialize

AMP on the first grid row to start the row marching of the parabolic wave scheme.




University of Florida Home Page
© 2004 - 2011 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Powered by SobekCM