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.

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