Home

REF/DIF 1 - University of Delaware

1. i 250 300 500 Figure 21 Waves interacting with a rip current Shoreline at right Wave height contours 69 4 4 Obliquely Incident Waves on a Plane Beach 70 5 REF DIF 1 Program Listing This section provides the listing for the Fortran code for REF DIF 1 which has been produced using the noweb documentation standard noweb provides a programming environment which allows the programmer to specify the operation of a block of code in full typeset detail after which the actual Fortran or C code is spelled out This procedure places a high premium on the use of a highly structured and modularized programming technique 71 5 1 Refraction Diffraction Model REF DIF 1 Version 3 0 REF DIF 1 calculates the forward scattered wave field in regions with slowly varying depth and current including the effects of refraction and diffraction The program is based on the parabolic equation method Physical effects included in the present version include 1 Parabolic approximation a Minimax approximation given by Kirby 1986b 2 Wave nonlinearity choice of a Linear b Composite nonlinear approximate model of Kirby and Dalrymple 1986b c Stokes nonlinear model of Kirby and Dalrymple 1983a 3 Wave breaking a Model of Dally Dean and Dalrymple 1985 4 Absorbing structures and shorelines a Thin film model surrounded by a natural surfzone Kirby and Dalrymple 1986a 5 Energy dissipation any of
2. J Geophys Research 90 11917 11927 Dalrymple R A J T Kirby and P A Hwang 1984a Wave diffraction due to areas of energy dissipation J Waterway Port Coastal and Ocean Engineering 110 67 79 Dalrymple R A J T Kirby and D W Mann 1984b Wave propagation in the vicinity of islands Proc of the 16 Offshore Tech Conf No 4675 Houston May Dalrymple R A 1988 A model for the refraction of water waves J Waterway Port Coastal and Ocean Engineering 114 423 435 Dalrymple R A 1991 REFRACT A refraction program for water waves Version 2 0 Report CACR 91 09 Center for Applied Coastal Research Dept of Civil Engrng Univ of Delaware Newark Dean R G and R A Dalrymple 1984 Water Wave Mechanics for Engineers and Scientists Englewood Cliffs Prentice Hall Djordjevic V D and L G Redekopp 1978 On the development of packets of surface gravity waves mov ing over and uneven bottom Z Angew Math and Phys 29 950 962 164 Hales L Z and Herbich J B 1972 Tidal inlet current ocean wave interaction Proc 13th ICCE 669 688 Hedges T S 1976 An empirical modification to linear wave theory Proc Inst Civ Eng 61 575 579 Houston J R 1981 Combined refraction diffraction of short waves using the finite element method Applied Ocean Res 3 163 170 Jonsson I G and O Skovgaard 1979 A mild slope wave equation and i
3. a Turbulent bottom friction damping b Porous bottom damping c Laminar boundary layer damping 6 Lateral boundary conditions either of a Reflective condition b Open boundary condition Kirby 1986c 7 Input wave field either of a Model specification of monochromatic or directional wave field b Input of initial row of data from disk file 8 Output wave field a Standard output b Optional storage of last full calculated row of complex amplitudes 72 The documentation of present program is contained in Combined Refraction Diffraction Model REF DIF 1 Version 3 0 Documentation and User s Manual James T Kirby Robert A Dalrymple and Fengyan Shi Report No CACR 02 02 Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware March 2002 2002 James T Kirby Robert A Dalrymple and Fengyan Shi This program is free software you can redistribute it and or modify it under the terms of the GNU Gen eral Public License as published by the Free Software Foundation either version 2 of the License or at your option any later version This program is distributed in the hope that it will be useful but WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY of FITNESS FOR A PARTICULAR PURPOSE See the GNU general Public License for more details You should have recieved a copy of the GNU General Public License along with t
4. 1 iinput 1 1 ioutput 1 input values of iinput ioutput 1 standard i e not starting from previous run 2 if starting from previous run 1 standard not saving restart data 2 if saving restart data 147 read iinput ioutput write iun 3 115 115 format input value of isurface 1 isurface 0 no surface picture generated 1 isurface 1 surface picture generated read isurface if isurface eq 0 fname6 write iun 4 nml ingrid if ispace eq 1 write iun 4 nml inmd if iinput eq 1 then c C write wavesl portion of indat dat c write iun 3 input iwave 1 discrete 2 directional spread read iwave write iun 3 input nfreq 4 of frequencies read nfregs write iun 4 nml wavesla if iwave eq 2 then write enter central direction thetO0 read thetO endif do 113 ifreq 1 nfreqs c Ex line 10 iinput 1 c write iun 3 109 109 format input wave period and tide stage read freqs ifreq tide ifreq c Ex line 11 iwave 1 iinput 1 c if iwave eq 1 then write iun 3 110 110 format input of waves per frequency nwavs read nwavs ifreq c G line 12 iwave 1 iinput 1 c do 111 iwavs 1 nwavs ifreq write iun 3 input amplitude and direction read amp ifreg iwavs dir ifreg iwavs 111 continue 148 else iwave 2 iin
5. 141 do 113 ifreq 1 nfreqs c line 10 iinput 1 write 109 109 format input wave period and tide stage read freqs ifreq tide ifreq EX line 11 iwave 1 iinput 1 if iwave eq 1 then write 110 110 format input of waves per frequency nwavs read nwavs ifreq Ex line 12 iwave 1 iinput 1 do 111 iwavs 1 nwavs ifreq write input amplitude and direction read amp ifreg iwavs dir ifreg iwavs 111 continue else c Er iwave 2 iinput 1 c write 112 112 format input en density and on next line directional 1 spreading factor read edens ifreq read nwavs ifreq nseed 500 endif 113 continue if iwave eq 1 write 10 nml waveslb if iwave eq 2 write 10 nml waveslc endif if iinput eq 2 then c GE line 9 iinput 2 c write input wave period and tide stage read freqin tidein 142 write 10 nml waves2 endif close 10 stop end 143 7 2 datgenv26 f datgen C C cx C C c c c c c c c c c c c c C cx c c c C cx datgenv26 f This program generates input data files for several exampl applications of REF DIF 1 In particular the first four cases listed here correspond to the four test cases shown in the User s Manual James T Kirby kirby udel edu 302 831 2438 FAX 302 831 12
6. As a result the output from REF DIF 1 has been almost completely restructured in Version 2 5 Values of wave height wave angle water depth and in the near future radiation stress components are stored in separate files at the reference grid resolution The complex amplitude data needed to construct a surface image is stored at the computational resolution The program surface f interpolates this data onto a regularly spaced rectangular grid and stores the surface image These files are described in section 2 9 Finally the user may store an estimate of the magnitude of the bottom velocity in a file bottomu dat 3 1 7 Changes Appearing in Version 3 0 Changes to version 3 0 stem from an increase in the number of variables provided as output and from initial work to make the program compatible with the NOPP Nearshore Community Model modeling system Some of the philosophy of that system is discussed in Section 1 2 The provisions in version 2 5 aimed at providing compatibility with the LRSS system have been dropped This has eliminated the separate need for the files infile f and infile2 f and the single option infile fis now imbedded where needed in the main code All references to the LRSS HDF libraries have been eliminated along with the alternate Makefile 23 3 2 Overview of Operating Manual This section provides a description of the program structure section 3 3 followed by some notes on prob lems which are likely to be enc
7. a m j dconv iu j 1 n nd write i n 3 x m dconv iu do 201 j 1 n a 1 3 a continue return format m rj psibar x f10 2 format 501 10 4 warning Ursell number f10 4 encountered at 1 grid location i6 16 format 1 should be using Stokes Hedges model 1 water format 1 used referenc grid row ir 13 117 pis phase psibar f20 4 ntype 1 due to shallow x direction subdivisions end 118 5 8 CTRIDA Tridiagonal matrix solution by double sweep algorithm Present subroutine adopted from the subroutine described in Carnahan Luther and Wilkes Applied Numerical Methods Wiley 1969 The original subroutine has been modified to handle complex array coefficients and solution values Input and output are a b c coefficients of row in tridiagonal matrix d right hand side vector of matrix equation U solution vector ii l beginning and end indices of positions in the dimensioned range of the column vector Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby September 1984 refdif subroutine ctrida ii l a b c d v include param h complex a iy b iy c iy d iy v iy beta iy gamma iy c Compute intermediate vectors beta and gamma beta ii b ii gamma ii d ii beta
8. 0 closed l open ibe O program picks x spacing 1 is minimum constant or variable x spacing 0 for constant 103 104 105 106 108 read mdc do 103 iko 1 mr 1 md iko mdc continue else write 104 format input md i for i l to mr 1 read md i i 1 mr 1 endif write 106 format input if 1 turbulent if 2 porous if 3 write standard choice 1 0 O0 read iff 1 iff 2 iff 3 write input isp subgrid features standard read isp write 108 format input values of iinput ioutput laminar 0 1 iinput 1 standard i e not starting from previous run ye 2 if starting from previous run 1 ioutput 1 standard not saving restart data J 2 if saving restart data read iinput ioutput write 115 115 format input value of isurface 1 isurface 0 no surface picture generated 1 isurface 1 surface picture generated read isurface if isurface eq 0 fname6 write 10 nml ingrid if ispace eq 1 write 10 nml inmd if iinput eq 1 then c EX write wavesl portion of indat dat c write input iwave 1 discrete 2 directional spread read iwave write input nfreg of frequencies read nfreqs write 10 nml wavesla if iwave eq 2 then write enter central direction thet0 read thetO endif
9. 112 c crested This routine was heavily modified by Raul Medina University of c Cantabria It is further modified by Arun Chawla to take out a one sided c derivative that introduced an asymmetry bias do 15 j 1 n if a m j eq 0 0 then akx2 0 else akx2 aimag clog a m J endif if a m 1 3 eq 0 0 then akx1 0 else akx1 aimag clog a m 1 3 endif if abs akx2 akx1 gt pi then akx sign 2 pi abs akx1 tabs akx2 dx akx1 else akx akx2 akx1 dx endif if j eq 1 then if a m j 1 eq 0 0 then aky2 0 else aky2 aimag clog a m 3 1 endif if a m j eq 0 0 then akyl 0 else akyl aimag clog a m J endif if abs aky2 akyl gt pi then aky sign 2 pi abs aky1 abs aky2 dy akyl else aky aky2 aky1 dy endif endif if j eq n then if a m j eq 0 0 then aky2 0 else aky2 aimag clog a m J endif if a m j 1 eq 0 0 then 113 akyl 0 else akyl aimag clog a m j 1 endif if abs aky2 akyl gt pi then aky sign 2 pi abs aky1 abs aky2 dy akyl else aky aky2 aky1 dy endif endif if j gt 1 and j 1t n then if a m j 1 eq 0 0 then aky2 0 else aky2 aimag clog a m 3 1 endif if a m j 1 eq 0 0 then akyl 0 else akyl aimag clog a m j 1 endif if abs aky2 akyl gt pi then aky sign 2 pi abs aky1 abs aky2 2 dy aky1 else ak
10. Pass Sxy i 1 3 1 3 1 yt iav TeS ee 1 3 1 3 yfl ys aye 1 e 1 3 1 3 radiation stresses for 3D circulation model Y d 2 dyr 1 2 dyr 1 2 dyr 1 2 dyr 1 1 dyr 1 1 dyr 1 1 dyr 1 1 dyr dxr dxr dxr dxr dxr dxr dxr dxr dxr dxr dxr dxr do i 1 nx wave aux 2 Pass_WaveNum i j Depth_wave i Jj capg aux SINH aux Sm i 3 1 16 dr i 3 1 capg grav amp Pass_Height i j Pass_Height i Jj Sp 1 3 1 16 dr 1 3 grav Pass_Height i j Pass_Height i j amp capg 1 enddo enddo Inside the surfzone do j 1 ny wave do i 1 nx wave if Pass ibrk i j eq 1 then aux 2 Pass WaveNum i j Depth wave i Jj capg aux SINH aux Sm i j 1 16 dr i j 1 capg grav amp Pass Height i j Pass Height i Jj t1 dr i j 0 06 Pass Height i j Pass C i j Pass C i j endif enddo enddo do j 1 ny wave do i 1 nx wave Pass_Sxx_body i j cos Thetar i j cos Thetar i Sm i jJ Pass_Sxy_body i j sin Thetar i j cos Thetar i j Sm i Jj Pass_Syy_body i j sin Thetar i j sin Thetar i j Sm i Jj Pass Sxx surf i j 2Sp i j Pass Sxy surf i j 0 Pass Syy surf i j 2Sp i j enddo enddo C print sxx sxy and Syy ibrk if fnamelO ne then do i 1 nx wave write 10 100 Pass Sxx i j j 1 ny_wave
11. end do do i 1 nx wave write 11 100 Pass Sxy i j J 1 ny_wave end do do i 1 nx wave write 12 100 Pass Syy i j J 1 ny_wave end do 131 endif c write Fx Fy if fnamel3 ne then do i 1 nx wave write 13 100 Pass Wave Fx i j j l ny_wave end do do i 1 nx wave write 14 100 Pass Wave Fy i j j 1 ny_wave end do endif Cc write ibrk if fnamel9 ne then do i 1 nx wave write 19 100 Pass ibrk i j j 1 ny_wave end do endif c write Quv if fnamel5 ne then do i 1 nx wave write 15 100 Pass MassFluxU i j j l ny_wave end do do i 1 nx wave write 16 100 Pass MassFluxV i j j l ny_wave end do endif c write tb if fnamel7 ne then do i 1 nx wave write 17 100 Pass_ubott i j j 1 ny_wave end do endif Cc write local radiation stresses for 3d circ model if fname21 ne then do i 1 nx wave write 21 100 Pass Sxx body i j j l ny_wave end do do i 1 nx wave write 22 100 Pass Sxy body i j J 1 ny_wave end do do i 1 nx wave 132 101 100 write 23 1 end do do i Lr 00 write 24 100 end do do i Lr write 25 100 end do do i Lr write 27 100 end do endif format 501 13 format 501 10 4 return end nx_wave nx_wave Pass_Syy_body nx_wave Pass_Sxx_sur Pass_Sxy_sur Fh Pass Syy sur
12. heightb 3 sgrt 0 25 height j nd 2 0 5 height j 2 1 0 25 height jtnd 2 end do endif Print out abs a at grid reference points mm1 m 1 write 5 205 ir 1 mm1 write 5 202 x m dconv iu psibar Wave heights on height dat write 26 203 heightb j dconv iu j 1 n nd Pass wave height mks unit Fengyan 02 04 2002 do 3j 1 n nd jje 3 1 nd 1 Pass Height ir 1 jj heightb 3 enddo Pass wave number do 3j 1 n nd 33 3 1 ndt 1 Pass_WaveNum irt 1 jj k m 3 enddo Wave angles on angle dat do 18 j 1 n 116 c 18 thet j continue write 7 203 1 80 thet 3 pi thet 3 j 1 n nd Pass angle Fengyan 02 04 2002 do 3j 1 n nd JJS enddo 1 nd 1 Pass_Theta ir 1 jj Water depths on depth dat Pass wave cg Bottom velocities on 19 write 8 203 do 3j 1 n nd JJS Pass_C ir 1 33 sig m Pass_ibrk ir 1 33 ibr enddo 1 nd 1 Pass_Cg ir 1 33 cg m Jj if fname7 ne then do 19 j 1 n nd bottomu m continue write 17 203 endif 3 bottomu m 3j cabs a thet j j k m j bottomu dat d m 3 dconv iu j 1 n nd c and ibrk Fengyan 02 04 2002 m 3 bottomu m 3 dconv iu j 1 n nd Write out reference grid data on disk file iun 3 Roll back solution to first grid level 201 202 203 204 205 write iu un 3
13. separate works But when you distribute the same sections as part of a whole which is a work based on the Program the distribution of the whole must be on the terms of this License whose permissions for other licensees extend to the entire whole and thus to each and every part regardless of who wrote it Thus it is not the intent of this section to claim rights or contest your rights to work written entirely by you rather the intent is to exercise the right to control the distribution of derivative or collective works based on the Program In addition mere aggregation of another work not based on the Program with the Program or with a work based on the Program on a volume of a storage or distribution medium does not bring the other work under the scope of this License You may copy and distribute the Program or a work based on it under Section 2 in object code or executable form under the terms of Sections 1 and 2 above provided that you also do one of the following a Accompany it with the complete corresponding machine readable source code which must be distributed under the terms of Sections 1 and 2 above on a medium customarily used for software interchange or b Accompany it with a written offer valid for at least three years to give any third party for a charge no more than your cost of physically performing source distribution a complete machine readable copy of the corresponding source code to be distri
14. A carefully controlled set of waves measurements has been made in the laboratory for this case see Berkhoff et al 1982 Kirby 19862 2 The wave pattern represents the case of ray crossing in the refraction method and thus the computed results indicate present method s utility in situations where ray tracing breaks down 3 The example gives a thorough test of the accuracy of the large angle and composite nonlinearity for mulations The topography to be studied is shown in Figure 15 Details of the calculation of the topography may be found in Kirby 1986a For the case of an incident plane wave we have performed a run of the model using input data correspond ing to the experiment of Berkhoff et al 1982 The run was done using the full model with the Stokes Hedges composite nonlinearity 59 Incident wave T Isec 0 ci crees Q E 0 45 Moso 5 x m lo l 15 20 Figure 15 Bottom contours and computational domain for the experiment of Berkhoff et al 1982 Experi mental data on transects 1 8 4 2 1 The Input Data Files The input data file indat dat for the present case follows Sfnames same as previous example Send Singrid mr 100 nr 100 iu 1 ntype icur tbe dxr 0 2500000 dyr 0 2500000 dt 10 00000 ispace nd iff isp iinput ioutput 1 Send Swavesla iwave 1 nfreqs Send Swaveslb freq
15. T continue return end 5 10 RANDI Generate a floating point pseudo random number between 0 and 1 by the multiplicative congruential method see Knuth D E 1969 p 155 Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by Robert A Dalrymple January 1986 refdif function randl x ix ifix 32767 x irand mod 4573 ix 6923 32767 randl float irand 32767 return end 121 5 11 RDFACT compute the factorial of xi Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by Robert A Dalrymple January 1986 refdif function rdfact xi prod 1 if xi gt 1 then do 17 ii 22 int xi prod prod float 11 17 continue endif rdfact prod return end 5 12 BNUM Compute the combination n n n in Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by Robert A Dalrymple January 1986 refdif subroutine bnum in n bn xin in xt rdfact xin xb rdfact float n rdfact float in n bn xt xb return end 122 5 13 ACALC Calculate the normalization factor a for the directional spectrum such that f ds cos 0 2 5 qg 1 where in code 6 thmax Center for Applied Coastal Research Department of Civil and
16. dx do 403 3 1 n y 3 float j 1 dy dy 2 do 404 j 2 n 1 g sqrt y 3 6 096 y 3 do 404 s m if x i a 67 g d i m 0 4572 ope 10 67 g and x i le asta do 67 g x i 25 if x i gt 18 29 g d i e 0 1524 a do 405 i 1 m d i 1 d i 2 d i n d i n 1 continue endif if itype eq 8 then 18 29 g d i j c c surface piercing breakwater c 451 452 write input m n dx dy read m n dx dy breakwater tip radius xt 1 5 yt 4 do 451 j 1 n y 3 float 3 1 dy continue do 452 i 1 m x 1 float 1 1 dx do 452 j 1 n if x 1 1t xt then r sqrt x 1 xt 2 y 3 yt 2 dep 66 r 37 else dep 66 abs y 3 yt 37 endif if dep gt 36 dep 36 d i j dep continue endif if itype eq 9 then generate a pair of breakwaters with rounded heads on each 154 c side of a channel Gr xt beginning location for trunk of breakwater EX sl slope of the sides of the breakwater c write input m n dx dy xt sl do read m n dx dy xt sl do w float n 1 dy do 503 i 1 m x i float i 1 x if x i le xt then do 501 j 1 n y 3 float 3 1 dy if j 1t n 2 then resqrt x i xt 2 y j 2 if r eq 0 then cr 0 else cr abs x i xt r endif dep s1 r 05 Gx add bulbous head dep dep sl sqrt r cr else r sqrt x i xt 2 y 4 w 2 if r eq 0 then cr 0 else cr abs x i
17. dx sig itl Jj 1 2 ifilt damp 1 3 complx 1 0 alphn dx cp2 i j cmplx deltal dx v i 1 3 v 1 3 2 dy 1 b1 u2 bet i 3 u itl 4 v itl 4 u 1 3 v 1 3 dy sig i 1 3 1 1 deltal u2 u i41 j41 v i41 j41 u i 3 1 v 1 3 1 12 u i 1 3 v 1 1 3 2 dy sig 1 1 3 1 dx b1 bet i sig i 1 j v i 1 j sig i j v i j 2 dy sig i 1 j41 1 cmp1x 2 b1 dy dy k 141 3 k 1 3 1 b1 bet 1 3 dx 2 dy dy 2 dx 2 dy dy pv 1 1 3 1 pv i RA E o b1 sig 1 1 3 v 1 1 3 dy sig 1 1 3 1 ys 107 k 141 3 1k i j ifilt damp i Jj cp3 i 3 cmp1x deltal dx v i 1 3 v 1 3 2 dy b1 u2 bet 1 3 u 1 1 3 v 1 1 3 u 1 3 v 1 3 dy sig i 1 3 1 deltal u2 u 1 1 3 1 v 1 1 3 1 u 1 3 1 v 1 3 1 12 u 1 1 3 v 1 1 3 2 dy sig 1 1 3 1 1dx b1 bet 1 3 sig 1 1 3 v 1 1 3 sig i j v i j 1 2 dy sig 1 1 3 1 cmplx 2 b1 dy dy k 1 1 3 k i j de E 1 b1 bet 1 3 dx 2 dy dy deltap i j dx 2 dy dy pv i 1 3 1 pv itl j 1 sig itl j 1 4 cmp1x 0 1 b1 sig i 1 3 v i 1 3 1 dy sig 1 1 3 1 k 1 1 3 k 1 3 ifilt damp i 3 an po RA dv 1 3 sig i 1 3 sig 4 si M uu nud a0 k i 3 cg i j u i 3 12 cmplx 0 1 omeg b
18. dyr 11 continue 12 continue c Set number of x points and define x values 99 f ispace eq 0 then c ispace 0 program sets subdivisions do 13 j 21 n dref d 1 3 tide ifreg f dref 1t 0 001 dref 0 001 call wvnum dref u 1 3 fregs ifreq itime k 1 3 eps icdw 1 Jj 13 continue npts 0 sumk 0 do 14 j 1 n if d 1 3 9t 0 05 then sumk sumk k 1 3 npts npts 1 endif 14 continue kb 1 sumk float npts alw 2 pi kb 1 anw dxr alw np ifix 5 anw if np lt 1 np 1 md ir min ix 1 np if np gt ix 1 write 5 100 ir endif o ispace 1 user specified subdivision m md ir 1 dx dxr float md ir do 15 i 1 m x 1 xr 1r f loat 1 1 dx 15 continue E interpolate values on m row do 16 j 1 n nd d m j dr irtl j nd 1 Intp_eta_Wave ir 1 3 1 nd 1 u m j ur ir 1 nd 1 v m vr ir 1 j 1 nd 1 16 continue do 18 jj 2 nr do 17 j 1 nd 1 jjj nd jj 2 j 1 d m 333 dr ir 1 jj dr ir 1 jj 1 y jjj dyr 1 yr 33 dr 1r 1 33 1 yr 3 d us dyr u m 333 ur r1 jj ur r 1 33 1 y 333 dyr lt iyr jj ur irr1 jj 1 yr jj 1 ur ir l1 3 dyt v m 333 vr r 1 j3 vr ir 1 j3j 1 y 333 dyr 100 1 yr 33 vr ir 1 33 1 T1 continue 18 continue c Interpolate values in do 19 i 2 m 1 do 19 j 1 n yr jj 1 vr ir 1 j3j dyr x dire
19. so existing indat dat files will still work properly 3 1 3 Changes Appearing in Version 2 2 Version 2 2 includes a revised version of a dissipation filter which is used to damp out noise after the onset of breaking in the numerical computations This filter has been found to be much more suitable in applications to field situations than was the original filter In addition Version 2 2 also provides the capability to input the first row of complex amplitude A from a new external data file wave dat The procedures for specifying wave dat are described in section 2 4 of the manual 3 1 4 Changes Appearing in Version 2 3 Version 2 3 provides a provision to save the last subdivided row of complex valued amplitudes computed on the last model grid row in file owave dat This option is potentially useful if the user wants to perform a run in several segments owave dat could then be used to store the intermediate calculation required to initialize a subsequent model run The provision for using this option is described in section 2 5 Use of this option does not affect any internal model calculations 3 1 5 Changes Appearing in Version 2 4 Version 2 4 incorporates two revisions to the basic model scheme The first revision is an extension to the model equations to handle the Minimax approximation of Kirby 1986b as well as the Pad large angle approximation Booij 1981 Kirby 1986a This algorithm is not used yet in the released version of the
20. xt r endif dep s1 r 05 add bulbous head Cc dep dep sl sqrt r cr endif if dep gt do dep do if dep 1t 05 dep 05 d i j dep 501 continue else do 502 j 1 n if j 1t n 2 then yt 0 else yt w endif dep s1 abs y 3 yt 05 if dep gt do dep do d i j dep 502 continue end if 503 continue endif if itype eq 10 then c 155 generate a breakwater with rounded head with an orientation alpha degrees to the x axis xh yh locus of breakwater head sl slope of the sides of the breakwater do constant depth section 551 write input m n dx dy xh yh alpha sl do read m n dx dy xh yh alpha s1 do al alpha 3 1415927 180 write alpha al co cos al si sin al do 551 i 1 m x 1 float 1 1 dx do 551 j 1 n y 1 float 3 1 dy xp x i xh co y 1 yh si yp x i xh si y 1 yh co see if we are in front of the trunk if xp 1t 0 0 then r sqrt xp xp yp yp dep s1 r 10 else dep s1 abs yp 10 endif if dep gt do dep do d i j dep continue endif return end c error function hasting s method c function erfjk x dimension a 5 a 1 0 254830 a 2 0 284497 a 3 1 421414 a 4 1 453152 a 5 1 061405 t 1 1 0 327591 x rf3k 1 exp x 2 t a 1 t a 2 t a 3 t a 4 t a 5 return end 156 7 3 surface f surface c c c s
21. 0 icur 1 ibc 0 ibc 1 86 5 3 INWAVE Read in wave parameters Variable definitions iinput ioutput if iinput 1 iwave nfreqs freqs tide if iwave 1 nwaus amp dir if wave 2 edens nsp n nseed if iinput 2 freqs tide determine method of specifying the first row of computational values 1 input values from indat dat according to value of wave 2 input complex amplitude values from file wave dat determine whether last row of complex amplitudes are stored in separate file owave dat 1 extra data not stored 2 extra data stored in file owave dat input wave type 1 input discrete wave amplitudes and directions 2 read in dominant direction total average energy density and spreading factor number of frequency components to be used separate model run for each component wave frequency for each of n freqs runs tidal offset for each of n freqs runs number of discrete wave components at each of nfregs runs initial amplitude of each component direction of each discrete component in or degrees from the x direction not presently recommended total amplitude variance density over all directions at each frequency spreading factor in cos 0 directional distribution stored in nwavs Seed value for random number generator Integer value between 0 and 9999 Wave frequency for one run Tidal offset for one run All data is entered using the namelist convention C
22. 19716 Coded by James T Kirby October 1984 August 1997 refdif subroutine con ifreq ir include param h include common h c Constants eps 1 0e 05 g 9 80621 c Calculate constants do 1 i 1 m do 1 j 1 n call wvnum d i j u i j freqs ifreq itime k i j eps icdw i Jj Sig i j freqs ifreg itime k i j u i jJ akd k i j d i j q 1 3 1 akd sinh akd cosh akd 2 p i j q i j g tanh akd k i j dd 1 3 cosh 4 akd 8 2 tanh akd 2 8 sinh akd 4 bottomu 1 3 g9 k 1 3 2 freqs ifreq itime cosh akd 1 continue is shallow probably due to using artificial current fields We now qaaaqaa positive number if not and flag the change to the user if icur eq 1 then do 2 i 1 m do 2 j 1 n cgroup sqrt p i j q i j3 advect cgrouptu i j 103 We seem to have occasional problems in adverse currents when the water check to see if the total advection velocity is negative make it a small Cc c if advect le 0 then advect 0 1 utemp advect cgroup write 5 100 ir i j u i j utemp u i j utemp endif continue endif Calculate the dissipation term wl call diss Calculate the mean kb on each row 10 11 100 do 11 i 1 m npts 0 sumk 0 do 10 j 1 n if d i j gt 0 05 then sumk sumk k i j npts npts 1 endif continue if npts eq 0 then kb i k 1 1 else kb i sumk float npts en
23. Environmental Engineering University of Delaware Newark DE 19716 Coded by Robert A Dalrymple January 1986 refdif subroutine acalc thmax nsp a itn 2 nsp call bnum itn nsp bn a thmax bn 2 itn 1 sum 0 do 10 ik 1 nsp ki ik 1 call bnum itn ki bn 10 sum sum bn sin float nsp k1 thmax float nsp ki a a sum 2 itn 2 return end 123 5 14 WVNUM Calculate wavenumber k according to the form c ku gk tanh kd where d local water depth s 0 absolute frequency g gravitational acceleration constant u x component of ambient current eps e tolerance for iteration convergence i j indices in finite difference grid icdw switch 0 no convergence failures encountered 1 at least one convergence failure Solution is by Newton Raphson iteration using Eckart s approximation as a seed value Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby September 1984 refdif subroutine wvnum d u s k eps icdw i Jj include param h common ref2 dr ixr iyr ur ixr iyr vr ixr iyr iun 8 iinput t ioutput real k kn c Constants g 9 806 pi 3 1415927 k s s g sqrt tanh s s d g c Newton Raphson iteration do 1 ii 1 20 f s s 2 s k utk k u u g k tanh k d fp 2 s u 2 k u u g tanh k d g k d cosh k d 2 kn k f fp if abs kn k kn lt eps go to 2 k kn 1 continu
24. Legal Information This section provides information pertaining to copyright and warranty This information pertains to the free standing REF DIF 1 program and to the Nearshore Community Model framework as a whole 1 1 1 Limits to Liability REF DIF 1 is provide as is without representation as to its fitness for any purpose and without warranty of any kind either express or implied including without limitation and implied warranties of merchantability and fitness for a particular purpose The University of Delaware shall not be liable for any damages in cluding special indirect incidental or consequential damages with respect to any claim arising out of or in connection with the use of this software even if it has been or is hereafter advised of the possibility of such damages 1 1 2 The GNU Public License REF DIF 1 has been developed by James T Kirby Robert A Dalrymple and Fengyan Shi for the purpose of computing the combined refraction diffraction of monochromatic surface water waves in the coastal envi ronment The program is copyright C 2002 by James T Kirby Robert A Dalrymple and Fengyan Shi The particulars of the GNU Public License follow TERMS AND CONDITIONS FOR COPYING DISTRIBUTION AND MODIFICATION 0 This License applies to any program or other work which contains a notice placed by the copyright holder saying it may be distributed under the terms of this General Public License The Program be low ref
25. Wave breaking e td ee a e dan Ge ee he ace ebd ge 14 2 4 Wave Climate cecep emm eS ee ae eee ee be oS am Pee ee eS 14 2 4 1 Monochromatic waves lee 14 2 4 20 Discrete directional waves not presently recommended 15 2 4 3 Directional spectrum not presently recommended 15 2 57 Model Output mes eto aues Red Ide RAS Bhd uS NOE Te SO Bet IRE eng 16 23 11 Complex Amplitude 4 spun a XX ves ES ebd uk 16 2 5 2 Wave Heights and Angles e 16 2 5 3 Radiation Stresses and Forcing Terms 2200 16 2 5 4 Wave Induced Mass Flux 000000000000 17 2 5 5 Velocity Moments for Bottom Stress Calculation 18 2 6 Numerical Development an e vase y e Ast BOY RES SORE IR 18 2 6 1 Crank Nicolson Technique 0 0 0 0 0 0000000004 18 2 6 2 Initial and Lateral Boundary Conditions llle 19 2 0 3 SUBEEIS s oe on yA RRR A 19 2 6 4 Damping of Spurious Computational Modes oo o 19 3 USER S MANUAL 20 3 1 REF DIF 1 Revision History 000000000000 000008 20 3 1 1 Changes Appearing in Version 2 0 o a 20 3 1 2 Changes Appearing in Version 2 1 o o 21 3 1 3 Changes Appearing in Version 2 2 a 21 3 1 4 Changes Appearing in Version 2 3 ees 21 3 1 5 Changes Appearing in Version 2 4 es 21 3 1 6 Changes Appearing in Version 2 5 22e 21 3 1 7 Changes Appearin
26. agk Cg 4 U A i p V i i i UV 5 b A D APA3 Za 4 tp v5l 2 k CTY A A A b 2iwU 2icV 20V oO d oO y oO xy 22 id iok 2 detoneeena 8 aafe 8 e kwU a 1 6 where ky k p U Ay OE 7 p 2k p U 0 Ai a b1 8 A 1 2a1 2b4 9 k N 01 b 10 k and w is a dissipation factor discussed in the next section The coefficients ap a and b depend on the aperture width chosen to specify the minimax approximation see Kirby 1986 The combination Q0 1 Qa 0 5 b 0 11 recovers Radder s approximation while the choices ag L a 0 75 TERRE S 12 recover the approximation of Booij 1981 The values of ay a1 and b used for the Minimax approximation depends on the range of angles to be considered Kirby 1986b found that the values for a maximum angular range of 70 gave reasonable results over the range of angles typically used The corresponding coefficient values for this choice are ao 0 994733 a 0 890065 bi 0 451641 13 Equation 6 is the model equation used in REF DIF 1 At present the model still utilizes the Pad approxi mation form based on the coefficients in 12 testing is underway to extend the model to the minimax model with coefficients 13 In the previous two equations the dispersion relationship relating the angular frequen
27. and construct various plots This program uses the quiver routine from Matlab 4 2 James T Kirby Center for Applied Coastal Research University of Delaware Newark DE 19716 kirby udel edu 302 831 2438 FAX 302 831 1228 JP o9 o9 o9 o oe o o o AP AP AP oe oe 11 27 94 oe Read data files load height dat load angle dat load depth dat dx input enter dx dy input enter dy sz size height oe Compute x y vectors x dx 1 sz 1 dx y dy 1 sz 2 dy Constructing scales arrow plot for wave heights and directions Compute x y components of arrows oe oe DX height cos pi angle 180 DY height sin pi angle 180 A Now do contours of wave height over depth contours figure 1 clf hold off cf contour x y height clabel cf manual xlabel x ylabel y hold on contour x y depth axis equal oe Now overlay scaled arrows on contours of wave height figure 2 clf hold off cf contour x y height clabel cf manual xlabel x ylabel y hold on quiver x y DX DY axis equal oe Now plot the surface 160 isurf input do you want to plot the surface l yes if isurf 1 load surf dat figure 3 clf hold off cf contour x y surf clabel cf manual xlabel x ylabel y hold on contour x y depth axis equal end 161 8 FREQUENTLY ASKED QUESTIONS 1 What is the maximum cur
28. common dims x ixr y iyr dimension iun 3 dimension amp ncomp ncomp dir ncomp ncomp tide ncomp 1 freqs ncomp edens ncomp nwavs ncomp write iun 3 1 format parabolic model in rectangular 1 grid KR KKK se s x XA ug 2 input type of bottom desired 3 l surface piercing island 4 2 bbr submerged shoal 5 3 arthur rip current 6 4 test case planar bottom 7 8 9 ih 2 3 5 radder 1979 configuration 2 6 grazing incidence on linear caustic T whalin s channel 8 surface piercing breakwater 9 channel 10 breakwater read itype if itype eq 1 then c surface piercing island c 101 write iun 3 101 format surface piercing island 150 102 103 104 105 106 write iun 3 102 format input m n dx dy depth period read m n dx dy dep t write iun 3 m n dx dy dep t write iun 3 103 format input crest height x semiaxis y semiaxis read hb xa ya write iun 3 hb xa ya sig 2 3 1415927 t do 104 i 1 m x i float i 1 dx continue do 105 j 1 n y 3 float 3 1 0 5 dy continue xc xa 3 dx do 106 i 1 m do 106 j 1 n d i j dep 1 sqrt x i xc xa 2 y 3 ya 2 1 hb if d i j gt dep d i j dep continue end if if itype eq 2 then c c bbr submerged shoal LL 152 153 154 iu
29. difference equation nonlinear terms are approximated on a first pass by using the A values Once the A 1 terms are computed the equation is solved again for A 1 using now the just calculated values in the nonlinear terms This two pass iterative method insures that the nonlinearities in the model are treated accurately Kirby and Dalrymple 19832 The solution proceeds by moving one grid row in the x direction incrementing 7 by one and using the two pass implicit implicit technique determining the complex amplitude A 1 for all the values of j on this row Marching in the wave direction these calculations are repeated until all the A are known for all and 7 While it appears that the Crank Nicolson procedure could be time consuming given that there is a matrix inversion for each grid row the coefficient matrix size is only 3 by n and the matrix inversion procedure is in fact very fast The procedure is economical on storage as only the values for the rows 2 and i 1 are necessary at each calculation 2 6 2 Initial and Lateral Boundary Conditions The initial condition is vital for the parabolic model The furthest seaward grid row corresponding to 2 1 is taken as constant depth and the incident wave s is prescribed here This wave is then propagated over the bathymetry by the model The various initial conditions were discussed in the section Wave Climate As in the solution of any differential equation in a domain th
30. direction and the y axis pointing towards the left 163 9 REFERENCES Arthur R S 1950 Refraction of shallow water waves the combined effect of currents and underwater topography Trans AGU 31 549 552 Berkhoff J C W 1972 Computation of combined refraction diffraction Proc 13th Int Conf Coastal Engrg Vancouver Berkhoff J C W N Booij and A C Radder 1982 Verification of numerical wave propagation models for simple harmonic linear waves Coastal Engineering 6 255 279 Bettess P and O C Zienkiewicz 1977 Diffraction and refraction of surface waves using finite and infinite elements Int J for Numerical Methods in Engrg 1 1271 1290 Booij N 1981 Gravity Waves on Water with Non uniform Depth and Current Doctoral dissertation Tech nical University of Delft The Netherlands 131 pp Booij N 1983 A note on the accuracy of the mild slope equation Coastal Engineering 7 191 203 Carnahan B H A Luther and J O Wilkes 1969 Applied Numerical Methods Wiley Chawla A Ozkan Haller H T and Kirby J T 1998 Spectral model for wave transformation over irregular bathymetry Journal of Waterway Port Coastal and Ocean Engineering 124 189 198 Chu V C and C C Mei 1970 On slowly varying Stokes waves J Fluid Mech 41 873 887 Dally W R R G Dean and R A Dalrymple 1985 Wave height variations across beaches of arbitrary profile
31. do ire 1 num data dis ire freqs 1 ire enddo do ire 1 num total between Master dt N Interval CallWave ire AX UPDATE INTERVAL partial between int between nstart int between 1 if nstart 1t num_data then freqs 1 ire dis nstart 1 partial dis nstart 1 partial else freqs 1 ire dis num data endif enddo c amp do ire 1 num data dis ire amp 1 ire enddo do ire 1 num_total between Master_dt N_Interval_CallWave ire X UPDATE_INTERVAL partial between int between nstart int between 1 if nstart 1t num_data then amp 1 ire dis nstart 1 partial dis nstart 1 partial else amp 1 ire dis num_data endif 90 enddo c angle 100 101 102 103 104 105 106 107 108 do ir dis e l1 num data ire dir 1 ire enddo do dr bet e num total ween Master dt N Interval CallWave ire par nst if dir els dir end enddo retur forma forma forma forma forma forma 1 fre forma de 5 deb forma wav forma 1n i end UPDATE INTERVAL tial between int between art int between 1 nstart 1t num_data then 1 ire dis nstart 1 partial dis nstart 1 partial e 1 ire dis num data if n t 1514 t 501 f10 4 t 1 20x input section wave data values t iwave 1 discrete wave amps and directions t
32. i Jj 126 aux 2 wave_number Depth_wave i j capg aux SINH aux Sm 1 3 1 16 1 capg grav rho amp Pass_Height i j Pass_Height i j Sp i j 1 16 grav rho amp Pass Height i j Pass_Height 1 3 capg enddo enddo Inside the surfzone do j 1 ny wave do i 1 nx wave if Pass ibrk i j eq 1 then wave number 2 pi period Pass C i Jj aux 2 wave number Depth wave i J capg aux SINH aux Sm 1 3 1 16 1 capg grav rho amp Pass Height i j Pass Height i Jj E 0 06 Pass_Height i j Pass_C i j Pass_C i j endif enddo enddo do j 1 ny wave do i 1 nx wave Pass Sxx i j cos Thetar i j E cos Thetar i j Sm i j Sp i jJ Pass_Sxy i j Ssin Thetar i j cos Thetar i j Sm i 3 Pass_Syy i j sin Thetar i j sin Thetar i j Sm i j Sp i Jj if Pass ibrk i j eq 0 then Pass MassFlux i j grav Pass Height i j Pass Height i j Pass C i j BO else Pass MassFlux i j grav Pass Height i j Pass Height i j Pass C i j BO tPass Height i j Pass C i j 0 06 endif enddo enddo transition using spline function only for x direction e goto 911 look for lst breakpoint along each each j suppose wave from left side dx dxr do j 1 ny wave index 1 iflag 0 127 do i 1 nx wave if Pass ibrk i j 1 and iflag 0 then index i iflag 1 endif enddo index index 1 slope Depth Wave index 1 j Depth wave
33. i nx_wave Pass_Wave_Fx 1 Pass_Wave_Fy 1 enddo i 1 j 1 Pass Wave Fx i Pass Wave Fy i i nx wave j 1 Pass Wave Fx i Pass Wave Fy i i 1 j ny_wave Pass Wave Fx i Pass Wave Fy i i nx wave j ny wave Pass Wave Fx i Pass Wave Fy i and rho do j 1 ny wave do i 1 nx wave Pass Wave Fx i Pass Wave Fy i enddo enddo local i j Pass F Pass Sxy i r i n jJ Pass_ F Pass Sxy i J Pass Pass_Syy i i j 1 3 1 J Pass_ F Pass Sxy i J Pass j Pass_ j Pass_ Outside the surfzone do j 1 ny_wave Pass_ F Pass Sxy i Pass Pass Syy i Pass Pass Syy i Pass F Pass Sxy i Pass Pass Syy i Sxy i Pass F Pass Sxy i Pass Pass Syy i _Sxy i Pass Syy i Sxx i 1 j J 1 Pass_Sxy i Sxy i 1 j J 1 Pass_Syy 1 Sxx i pat 1 Sxy i J 1 i j Pass Sxy i rd Pass Syy i Sxx i tl j Pass oe jt1 Pass_Sxy i Sxy i 1 j Pass ET J3 1 Pass_Syy i Sxx i 3 1 13 Pass Sxy i 13 Pass Syy i J JEL Sxx i 1 jJ Pass_Sxx J Pass_Sxy i 3 1 Sxy 1 1 3 Pass_Sxy 1 J Pass_Syy 1 3 1 Sxx i 13 i j Pass Sxy i rd Pass Syy i Jal 13 3 1 rho rho Wave Fx Wave Ea 130 Pass Sxx i j Jl Pass_Sxy i j Jal Pass Sxx i JL Pass Sxy i jb Pass Sxx i Pass Sxy i Pass Sxx i
34. if an input error is detected 36 e ismooth A value of ismooth 1 causes the program to attempt a smoothing operation applied to the output waveheight wave angle radiation stresses and forcing term calculations This is intended to mimic the fact that these quantities when computed as statistical averages in a random sea state would not have the localized spatial variability that they have in a representative monochromatic sea The averaging is intended to provide a more realistic estimate of the model output in the context of a random sea state The chosen averaging window is however ad hoc and is likely to be optimal only in a limited range of conditions This option is not viewed as a substitute for running a real spectral model such as REF DIF S e dxr dyr Reference grid x spacing and y spacing which are assumed to be uniform over the entire reference grid e dt Depth tolerance value e ispace nd ispace is switch controlling subdivisions ispace 0 program attempts its own x subdivisions ispace 1 user specifies x subdivisions e nd nd is the number of y direction subdivisions needed in either case e ffl iff 2 iff 3 Dissipation switches iff 1 1 turn on turbulent boundary layer iff 2 1 turn on porous bottom damp ing iff 3 1 turn on laminar boundary layers No damping if all values are zero e isp Switch for user specified sub grid specifications isp 0 no subgrids to be read isp 1 subg
35. iwave 2 directional spreading model chosen t the model is to be run for i3 separate quency components t wave component i2 amplitude f8 4 rection f8 4 t frequency component i2 e period f8 4 sec tidal offset f8 4 t total variance density f8 4 spreading factor 2 seed number i5 91 5 4 MODEL This subroutine is the control level for the actual wave model Data read in during inref and inwave is conditioned and passed to the wave model This routine is executed once for each frequency component specified in inwave The wave model is split in three parts which are run sequentially for each reference grid row grid subroutine performs the interpolation of depth and current values con calculate the constants needed by the finite difference scheme fdcalc perform the finite difference calculations calculate_wave forcing calculate radiation stresses depth integrated forcing and local forcing Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby November 1984 Last modified March 2002 refdif subroutine model include param h include common h include pass h integer i j dimension dthi 31 thi 31 thet iy C DON t THINK THESE ARE NEEDED dimension sxy iy sxx iy syy iy c Constants g 9
36. j dr i 1 j dr i j 1 dr 1 3 1 4 if abs dcheck dr i j gt dt write 5 104 dr i j i 3 dt 8 continue if icur eq 1 then do 9 i 1 mr do 9 j 1 nr if dr i 3 1e 0 0 go to 9 fr ur 1 3 ur 1 3 vr i j vr 1 3 g dr 1 3 if fr gt l write s 105 1 3 r 9 continue endif c Establish coordinates for reference grid do 10 ir 1 mr xr ir float ir 1 dxr 10 continue do 11 jr 1 nr yr jr float jr 1 dyr 11 continue c Establish y coordinates for interpolated grid n nd nr 1 1 dy dyr float nd do 12 j 1 n y j 1loat j 1 dy 12 continue c Check friction values E iff 1 1 turbulent boundary layer damping everywhere C iff 2 1 porous bottom damping everywhere c iff 3 1 laminar boundary layer damping everywhere do 13 i 1 3 if iff i ne 0 and iff i ne 1 iff i 0 ES continue write 5 116 1f f 1 i 1 3 c Specify whether or not user specified subgrids are to be read in during c model operation C lisp 0 no subgrids specified 84 isp 1 subgrids to be read in later from unit 3 if isp eq 0 write 5 117 if isp eq 1 then write 5 118 open unit 3 file fname3 endif if isp eq 1 and ispace eq 0 write 5 113 if isp eq 0 then do 14 ir 1 mr do 14 jr 1 nr isd ir jr 0 14 continue else C CHECK UNIT NUMBER HERE SUBDAT N I EDS TO BE OPENED TOO do 15 ir 1 mr 1
37. problem section Several things must be kept in mind while constructing a namelist oriented data file namelist does not allow a part of an array to be read or written It is therefore imperative that the number of values being entered into any dimensioned variable in a namelist group be equal to the dimension of that variable in REF DIF 1 It is thus highly recommended 40 that any program being used to construct the indat dat file should use the param h file to specify parameters In any event the user should consult the datgenv26 f or indat createv26 f programs to get a feel for how the file is constructed The file may also be easily constructed by hand In constructing the indat dat in the provided programs we have followed the most restrictive convention that was found which is that the items in a namelist group must be read in in the order in which they are specified in the namelist statement Most compilers will allow an arbitary ordering of variables within each namelist group 41 3 8 Program Input Reference Grid and Subgrid Data This section describes two files which provide arrays of data on the various grids Reference Grid Data File This file usually named refdat dat but freely chosen as fname on input is accessed as logical device number iun 1 Its contents consist of the arrays of depth dr x velocity ur and y velocity vr at the reference grid points This file is accessed only once per model run and its entire co
38. sig i j 1 cmplx 2 b1 dy dy k 1 1 3 1k i j bl bet 1 3 dx 2 dy dy deltap i j dx J 2 108 1 2 dy dy 1 pv i j pv i j 1 sig i j 1 ifilt damp i j fl i j tanh k i j d i 3 5 2 i j k i j d i j3 sinh k i j d i 3 4 c Constants defining the parabolic model angular aperture 70 degree minimax coefficients a0 0 994733 al 0 890065 b1 0 451641 qaaaqadaaa c Pad e coefficients a0 1 0 al 0 75 b1 0 25 c Additional constants u2 1 0 kap 0 78 gam 0 4 omeg freqs ifreq itime pi 3 1415927 ci cmp1x 0 1 c cdamp is the constant for the noise suppression formulation It can be c changed cdamp 0 00025 alphn 0 deltal al b1 delta2 1 2 a1 2 b1 c Initialize breaking index if ir 1 if ir eq 1 then ifilt 0 do 100 j 1 n ibr 3 0 wb 1 3 cmp1x 0 0 100 continue endif c Solution for m grid blocks in reference block ir do 200 i 1 m 1 109 ipl i 1 it 1 ih 1 C r h s of matrix equation rhs 1 cmp1x 0 0 do 1 3 2 n 1 rhs 3 c1 1 3 a i a 3 a i j 1 c3 1 5 a i j 1 l dx w a 1 3 wb 1 3 a 2 1 dx cmplx 0 iun anl de j k i yt k 1 j dd i 1 1 float ibr e cabs a i 3 2 ati 3 2 1 dx cmplx 0 l f1 Quae c M anl 1 sig i AO me 13 k 1 3 k 1 3 cabs a 1 3 2 1 dd i ah keen Ue 3 d i Dez k i j cabs a i 3 1
39. sxys dat surface stress part of local radiation stresses Syy The values are always zero If REF DIF 1 is given a null string as the input for this file name no file is generated e fname26 syys dat surface stress part of local radiation stresses Syy If REF DIF 1 is given a null string as the input for this file name no file is generated Variable definitions 78 mr nr dar dyr ju dt ispace nd md mr 1 ntype icur ibc ismooth dr ur Ur reference grid dimensions max izr yr grid spacing for reference grid physical unit descriptor 12mks 2 english default value is 1 mks units depth tolerence value to check for anomalous depth values switch to control grid subdivision 0 program attempts its own subdivisions 1 user specifies subdivisions y direction subdivision space 0 or 1 must be 1t y nr 1 x direction subdivisions if space 1 must be e ix 1 nonlinearity control parameter 0 linear model 1 Stokes matched to Hedges in shallow water 2 Stokes throughout switch to tell program if current data is to be used and read on input 0 no input current data 1 input current data to be read program defaults to cur 0 boundary condition switch 0 use closed lateral boundaries 1 use open lateral conditions program defaults to bc 0 artificial smoothing switch 0 no additional artificial smoothing 1 a fairly wide averaging wind
40. tanh k i j d i ps Lyre 2 1 continue rhs n cmp1x 0 0 c Return here for iterations 2 if it eq 1 ii i if it eq 2 ii ipl c Establish boundary conditions if ibc eq 1 then ksthl real 2 a i 2 a i 1 a i 2 a i 1 dy cmplx 19 zT ksth2 real 2 a i n a i n 1 a i n a i n 1 dy lecmplx 0 1 Er cmplx 1 ksthl dy 2 c 1 cmp1x 1 ksth1 dy 2 c n emplx 1 ksth2 dy 2 c n cmplx 1 ksth2 dy 2 else c 1 cmpl1x 1 0 c 1 bc 1 c n cmplx 1 0 c n bc n er c Calculate dissipation in rows where breaking occurs do 3 j 1 n if ibr j eq 1 wb 2 j cmplx 1 0 0 15 cg ip1 j 1 1 gam d ip1 3 2 cabs a 11 3 2 d ipl 3 if ibr j eq 0 wb 2 3 cmpl1x 0 0 3 continue 110 c Coefficients for forward row do 4 3 2 n 1 ac 3 cp3 i j bc 3 cp1 1 3 dx 2 w 1 1 3 wb 2 3 cmp1x 0 an anl 1 sig i 1 7 1 k 141 3 k 1 1 3 dd 1 1 3 cabs a 11 3 2 dx 2 1 cmp1x 0 an 1 an1l sig 1 1 3 dx 2 1 1 i 1 3 1k 1 1 3 k 1 1 3 cabs a 11 3 2 dd 1 1 3 tanh k 1 1 3 1d 1 1 3 2 1 1 3 k 141 3 cabs a 11 3 tanh k 141 3 d itl 3 13 Ta cc 3j cp2 i j 4 continue c Update solution one step call ctrida 1 n ac bc cc rhs sol do 5 j 1 n a ipl j sol j sol j cmp1x 0 0 5 continue if it eq 2 or ih eq 2 go to 8 c Check
41. the onset of breaking the REF DIF 1 model is able to represent waves both outside and inside of a surf zone The wave breaking algorithm is always active in the model Large surface piercing islands and causeways which would have surf zones are handled by the thin film technique of Dalrymple Kirby and Mann 1984b and Kirby and Dalrymple 1986a This procedure permits the easy computation of wave heights around arbitrarily shaped islands by replacing islands with shoals of extremely shallow depth 1 cm The wave breaking routine reduces the wave heights over the shoal to less than one half centimeter which results in a wave which carries negligible energy and therefore no longer affects any physical processes Thus the REF DIF 1 model does not distinguish between islands and deeper water computationally However the model output clearly shows the influence of the islands as will be shown in section 3 Examples of wave breaking and the combined refraction diffraction model appear in Kirby and Dalrymple 1986a and Dalrymple et al 1984b 2 4 Wave Climate 2 4 1 Monochromatic waves While the REF DIF 1 model is typically used with monochromatic wave trains propagating in one given direction there is no intrinsic restriction to this case As an example for a given frequency the wave di rection is determined by the initial wave height distribution provided by the user on the offshore grid row corresponding to x 0 As this row is parall
42. with one group corresponding to each value of 1 in isd The program accesses the data file by x row ir and then by y column jr For the above example the subgrid array groups should thus be stored in the order of the coordinate pairs 4 3 4 4 5 3 5 4 The dimensions of each subgrid are given by m md ir 1 and ns nd 1 The borders of adjacent subgrids share the same common boundary points Each group of subgrid data should be written using a format similar to write unit d i j j 1 ns i 1 m then if icur 1 write unit u i j j 1 ns i 1 m write unit v i j j 1 ns i 1 m format 20 10 4 using the appropriate value of m for the grid block in question Data may be in MKS or English units depending on the value of iu in the input data file The integer array isd is read by inref and the individual subgrids are read by grid 3 9 Program Output This section discusses output of two forms output sent to a log file and array output stored in disk files 3 9 1 Output log file Log file output consists of three types title and header information which is printed in order to identify the run and the operating perameters runtime error messages either warnings or terminal error messages and information about calculations on each grid row The default name for the file is refdifl log Header Information At the start of each run a set of two title pages are printed Page 1 ide
43. work based on the Program and copy and distribute such modifications or work under the terms of Section 1 above provided that you also meet all of these conditions a You must cause the modified files to carry prominent notices stating that you changed the files and the date of any change b You must cause any work that you distribute or publish that in whole or in part contains or is derived from the Program or any part thereof to be licensed as a whole at no charge to all third parties under the terms of this License c If the modified program normally reads commands interactively when run you must cause it when started running for such interactive use in the most ordinary way to print or display an announcement including an appropriate copyright notice and a notice that there is no warranty or else saying that you provide a warranty and that users may redistribute the program under these conditions and telling the user how to view a copy of this License Exception if the Program itself is interactive but does not normally print such an announcement your work based on the Program is not required to print an announcement These requirements apply to the modified work as a whole If identifiable sections of that work are not derived from the Program and can be reasonably considered independent and separate works in themselves then this License and its terms do not apply to those sections when you distribute them as
44. 0 0000 wave component 1 amplitude 0 5000 direction 0 0000 Figure 8 Sample title page 2 Dimensional quantities are printed on the title page in the units used for input Title page gives an indication whether quantities are in MKS or English units Run time Error Messages REF DIF 1 performs some data checking and checking of calculations during a run This checking may result in warnings or terminal errors which are beyond calculation errors which would lead to standard FORTRAN error messages A list of possible errors and the resulting messages follow 1 Reference grid dimensions were specified as being too large on input mr gt ixr and or nr gt iyr Message dimensions for reference grid too large stopping Action Program stops Error occurs in inref 2 User specifies a y direction subdivision nd which will cause the number of y grid points n to exceed the maximum iy Message y direction subdivision too fine maximum number of y grid points will be exceeded Exe cution terminating Action Program stops 45 Error occurs in inref User specifies an x direction subdivision on one of the grid blocks ir which exceeds the maximum amount ix 1 As a result the dimension of the subdivided grid will be too large Message x direction subdivision too fine on grid block ir execution terminating Action program stops Error occurs in inref A depth value occurs in the reference grid which differs fro
45. 1 bet 1 3 u 1 1 3 u i j sig i j 14 cmp1x 0 1 omeg b1 3 u 1 1 3 u i j dx v 1 1 3 1 lv i 3 1 v itl j 1 v i j 1 4 dy sig i j k i 1 3 k 1 3 1 cmplx 2 b1 dy dy k 1 1 3 k 1 3 b1 bet 1 3 ldx 2 dy dy deltap 1 3 dx 2 dy dy pv 1 3 1 12 pv i j pv i j 1 sig i j cmplx 1 0 omeg delta2 1 3 u 1 1 3 u 1 3 2 sig 1 3 ci omeg a0 1 k i j tu i 3 dx siq i 3 2 ifilt damp i j cmplx 1 0 alphn dx c2 1 3 cmp1x deltal dx v i 1 3 v i j 2 dy b1 u2 bet 1 3 u itl 3 v 1 1 3 u i j v i j3 dy sig i j 1 1 deltal u2 u 1 1 3 1 v 1 1 3 1 u i j 1 v i j 1 1 2 u i j v i j 2 dy sig i j 1 4 b1 sig i j v i j 1 dy k 141 3 k 1 3 sig 1 3 1 dx b1 bet 1 4 lsig itl j v itl 3 sig i j v i j 2 dy sig i j 1l lemplx 2 b1 dy dy k itl 3 k 1 3 bl bet i j dx 2 dy dy deltap i j dx 2 dy dy pv i j 1 pv i 3 sig i j 1 1ifilt damp i j c3 i 4 cmplx deltal dx v i 1 3 v 1 3 2 dy b1 u2 bet 1 3 1 u 141 3 v 141 3 u i j v i j dy sig i j 1 a 1 u itl j 1 v 1 1 3 1 u 1 3 1 v 1 3 1 u i j v i j 12 dy sig i j 1 4 b1 sig 1 3 v 1 3 dye cred dy i 1k i j sig i j 1 dx b1 bet 1 3 sig 1 1 3 v 1 1 3 lsig i j v i j 2 dy
46. 1 m 100 n 100 c20 cos 20 3 1415927 180 s20 sin 20 3 1415927 180 dx 0 25 dy 0 25 do 151 i 1 m x i float i 1 dx do 152 3 1 n y 3 float 3 1 dy do 154 i 1 m do 154 j 1 n xp x 1 10 5 c20 y 3 s20 yp x 1 10 5 s20 y 3 3 E20 test yp 4 2 xp 3 2 if xp 1t 5 84 d i j 20 45 if xp ge 5 84 d 1 3 0 45 0 02 xp 5 84 if test gt 1 go to 153 10 10 d i1 3 d 1 3 0 5 sqrt 1 yp 5 2 xp 3 75 2 1 0 3 continue continue end if 151 if itype eq 3 then c arthur 1952 rip current c 201 202 203 iu 1 icur 1 dx 5 0 dy 5 0 slope 0 02 m 100 n 100 do 201 i 1 m x i float i 1 dx continue do 202 j 1 n y j float 2 3 n 1 dy 2 continue xm x m do 203 j 1 n do 203 i 1 m d i1 3 xm x i dx slope continue endif if itype eq 4 then c c planar bottom test case c 251 252 253 write iun 3 input m n dx dy depth period read m n dx dy dep period write iun 3 input bottom slope read xm sig 2 3 1415927 period do 251 i 1 m do 251 j 1 n d i j dep xm float i 1 dx continue do 252 i 1 m x i float i 1 dx do 253 j 1 n y j float j 1 dy endif if itype eq 5 then c c radder 1979 configuration 2 c write iun 3 input m n dx dy depth read m n dx dy dep iu 1 ra
47. 133 i j i j i j i j ny wave ny wave ny wave ny wave 5 16 SPLINE1 refdif subroutine splinel nx ny xin yin n y2 This subroutine will compute the cubic spline for an array include pass h integer i j real xin nx_max yin nx_max y2 nx_max u nx_max y2 1 0 u 1 0 do i 2 n 1 sig xin 1 xin i 1 xin i 1 xin 1 p sig y2 i 1 2 y2 1 sig 1 p u i 6 dd y Sen yin i 1 a x i xin i 1 xin i 1 xin i 1 sig u i 1 p end do an 0 un 0 y2 n un gn u n 1 qn y2 n 1 1 do kk n 1 1 1 y2 kk y2 kk y2 kk 1 u kk end do return end 134 5 17 SPLINT1 This subroutine will compute values from the cubic spline refdif TI SUBROUTINE SPLINT1 nx ny xin yin y2 n xout yout index jump j include pass h integer i j real xin nx max yin nx max y2 nx max xout nx max real yout nx max ny max integer index jump do i 1 index 1 yout i j yin i enddo do i index jump 1 nx yout 1 3 yin i J ump 1 enddo do i index index jump KLO 1 KHI N 1 IF KHI KLO GT 1 TH kk KHI KLO 2 IF xin kk GT xout i THEN KHI kk ELSE KLO kk ENDIF GOTO 1 ENDIF hh xin KHI xin KLO IF hh EQ 0 then print Bad xin input stop endif aa xin KHI xout i hh bb xout i xin KLO hh yout i j aa yin KLO bb yin KHI A aa 3 aa y
48. 2 KLO bb 3 bb y2 KHI hh 2 6 E N end do return end 135 6 ADDITIONAL MODEL COMPONENTS 6 1 REF DIF 1 has been used as a wave driver wave module in a comprehensive community model developed with the support of the National Ocean Partnership Program NOPP Basically a master program is used as an interface to link three modules i e wave module circulation module and sediment module The mas ter program handles module coupling control internal data transfer and interpolation extrapolation between modules as well as input of control parameters and result output Passing variables are included in pass h master f The passing variables in REF DIF 1 are Depth_Wave water depth on Wave_Grid from the sediment module Pass_Sxx Pass_Sxy and Pass_Syy radiation stresses Pass_Sxx_body Pass Sxy body Pass_Syy_body local radiation stresses body part Pass Sxx surf Pass Sxy surf Pass Syy surf local radiation stresses surface part Pass Wave Fx Pass Wave Fy wave forcing Pass MassFluxU Pass MassFlux V mass flux Pass Diss dissipation caused by wave breaking Pass WaveNum wave number Pass Theta wave angle Pass ubott bottom velocity Pass Height wave height Pass C phase velocity Pass Cg group velocity Pass Period wave period Pass ibrk wave breaking index 1 breaking O nonbreaking Intp U Wave
49. 28 Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 January 1991 revised July 1994 for REF DIF 1 version 2 5 Revised February 2002 for version 2 6 Last revision 2 28 02 include param h dimension iun 4 md ixr common ref dr ixr iyr ur ixr iyr vr ixr iyr mr nr dxr dyr litype common ind iu ntype icur ibc ispace nd iff isp iinput liwave nfreqgs freqs tide nwavers amp dir edens common dims x ixr y iyr dimension iff 3 dimension amp ncomp ncomp dir ncomp ncomp tide ncomp 1 freqs ncomp edens ncomp nwavs ncomp character 255 fnamel fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell fnamel2 fnamel3 fnamel4 fnamel5 fnamel6 fnamel 7 fnamel8 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 data fname2 refdat dat name3 subdat dat fname4 wave dat name5 refdifl log fname6 height dat name7 angle dat fname8 depth dat name9 fnamel0 sxx namell sxy dat fnamel2 syy dat namel3 fx dat fnamel4 fy dat Fh FH Fh Fh Fh Fh Fh 144 fnamel5 qx dat fnamel6 qy dat fnamel7 tbx dat fnamel8 tby dat fnamel9 ibrk dat fname20 owave dat fname21 fname22 fname23 fname24 fname25 fname26 namelist ingrid mr nr iu ntype icur
50. 80621 rho 1000 pi 3 1415927 eps 1 0e 05 c find updated wave information itime max 1 Time Master N Interval CallWave Master dt itime min itime num total print wave itime itime height 2 amp 1 itime c Execute model once for each frequency 92 ifreq is the controlling index value do 200 ifreq 1 nfreqs psibar 0 write 5 203 ifreg period 2 pi freqs ifreq itime ismooth 1 compute window width if ismooth eq 1 then al0 g period period 2 pi nrs int 20 al0 dyr 1 2 write 5 window half width nrs endif c Specify initial nonlinear parameters for each run if ntype eq 0 an 0 if ntype ne 0 an 1 if ntype ne 2 anl 0 if ntype eq 2 anl 1 c Calculate the mean kb on the first row for use in c specifying initial conditions npts 0 sumk 0 do 10 jr 1 nr d 1 jr dr 1 jr ttide ifreq Intp_eta_Wave 1 jr call wvnum d 1 jr ur 1 jr fregs ifreg itime k 1 jr eps licdw 1 1 if d 1 jr gt 0 05 then sumk sumk k 1 jr npts npts 1 endif pass wave number the first row Pass_WaveNum 1 3r k 1 31 pass c cg on the first row sig 1 jr freqs ifreg itime k 1 jr ur 1 jr akd k 1 jr dr 1 jr q 1 3r 1 akd sinh akd cosh akd 2 p 1 jr q 1 jr g tanh akd k 1 31 Pass_C 1 jr sig 1 jr k 1 jr 93 c c c Cc c Pass Cg 1 31 sqrt p 1 jr q 1 jr 10 continue kb 1 sumk float npt
51. Combined Refraction Diffraction Model REF DIF 1 Version 3 0 Documentation and User s Manual James T Kirby Robert A Dalrymple and Fengyan Shi Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Research Report NO CACR 02 02 September 2002 Supercedes Version 2 5 Report No 94 22 1994 Contents 1 INTRODUCTION 1 1 1 dTegalInformation ia tal ee Eb ee e td y A 1 1 1 1 Limits to Liability o 0 00 00 0000000000000 1 1 1 2 The GNU Public License ie eee RA AI A A 1 1 2 Notes on Using REF DIF 1 in the NOPP Nearshore Community Model System 1 3 Document and Source Code Generation using NOWEB 6 2 THEORETICAL BACKGROUND 6 2 1 Wave Models orden be pure A Se Oe Soke we pe oe s 8 2 1 1 Mildslopeequation 2 5 58 55 4 mm RS RE eee a xEE 8 212 Diftracthonimodels oap weno Ge URS Sura Roue cda pe 9 2 1 3 Nonlinear combined refraction diffraction models 9 2 1 4 Wave current interaction models lees 10 2 24 ASSUMPLODS sg gb eee ee ee ee Bel o eB a a ee eed ix 12 2 3 Energy Dissipation sc ar Yeon em PAE eS ee he ew phe i fee 13 2 3 General fori s Li scd gw bw c Oy IS we ad Be eS Sha at 13 2 3 2 Laminar surface and bottom boundary layers o o 13 2 3 3 Turbulent bottom boundary layer lee 13 2 3 4 Porous Sand i G aa c Sue uuum ieget ts de 13 23 3
52. JI piso Uy adz 40 h 1 f toco f memet 41 Outside the surfzone S ag and Sa ap Can be determined using sinusoidal wave theory b 2 a H 1 42 Sag TA 1 G 42 o 656 pol G1 4 ag 4 055709 G 1 43 Inside the surfzone S o and S ag can be written as 2 b 2 C A h Sap eagpgH D E Tuy d 44 Sas bag pg Bo 1 45 For the depth integrated time averaged circulation model the short wave forcing terms are calculated by OS ag Fo 46 98 46 2 5 4 Wave Induced Mass Flux Outside the surfzone the wave volume flux is given by gH ka B ui 47 Quo TUE 47 where ka is the wave number vector in direction za and C is the value of the phase velocity Inside the surfzone the wave volume flux is given by Svendsen 1984 H C Ah ka wa Bot ml 4 Q E nm 48 17 2 5 5 Velocity Moments for Bottom Stress Calculation 2 6 Numerical Development 2 6 1 Crank Nicolson Technique The parabolic model is conveniently solved in finite difference form In order to accomplish this the study area bathymetry must be input as a grid with the a y directions divided into rectangles of Ax and Ay sizes The complex amplitude A x y will then be sought at each grid and therefore we can keep track of A by denoting its location not by x y but by i j where x i 1 Aa and y j 1 Ay Now we have to determine the values of A i j which satisfy Eq 6 for all i between 1 and m and for all j betw
53. THE POSSIBILITY OF SUCH DAMAGES END OF TERMS AND CONDITIONS 1 2 Notes on Using REF DIF 1 in the NOPP Nearshore Community Model System REF DIF 1 has been commonly used as a wave driver in conjunction with a number of wave induced nearshore circulation models At present a comprehensive community model is under development with the support of the National Ocean Partnership Program NOPP As this code is developed small adjust ments will be needed in the REF DIF 1 program in order to accomodate the needs of the overall modelling system Changes which are transparent to the users of REF DIF 1 as a stand alone program will not trigger a revision of the program documentation Notes on using REF DIF 1 in the context of the comprehensive system will appear here The NOPP model system is now undergoing preliminary development and documentation and will be described separately 1 3 Document and Source Code Generation using NOWEB The program source and documentation for REF DIF 1 are maintained using NOWEB which is described at http www eecs harvard edu nr noweb 2 THEORETICAL BACKGROUND The propagation of water waves over irregular bottom bathymetry and around islands involves many pro cesses shoaling refraction energy dissipation and diffraction Until recently only very approximate mod els existed to predict the wave behavior due to these effects This manual describes the weakly nonlinear combined refraction and diffraction
54. a 1 jj thet j enddo endif c Store first row of wave heights on unit CH if fname6 ne then write 26 202 2 cabs a 1 j dconv iu J 96 aky sign 2 pi abs aky1 abs aky2 dy aky1 002 ECK UNIT 12 tynna endif C Pass wave height mks unit Fengyan 02 04 2002 do 3j 1 n nd 33 3 1 nd 1 Pass height 1 33 2 cabs a 1 3 Pass_ibrk 1 33 0 enddo c Store surface on file fname9 if fname9 ne then x 1 0 write 9 n write 9 y j j 1 n write 9 x 1 write 9 real a 1 3 j 71 n endif c Now execute model for the ifreq frequency over each of mr c grid blocks ir is the controlling index value do 100 ir 1 mr 1 c Establish interpolated grid block for segment lirl call grid ifreq ir c If ir 1 write initial values on iun 3 if ir eq 1 then write 5 201 x 1 dconv iu psibar endif c Calculate constants for each grid block call con ifreq ir c Perform finite difference calculations call facalc ifreg ir c Grid block ir done print output and go to next grid 100 continue if ioutput eq 2 then write 20 a m 3 3 1 n 97 endif c Termination for the surface dat file if fname9 ne then x 1 100 write 9 x 1 endif c Calculate radiation stresses using new formula roller and spline c transition fyshi 02 04 2002 Pass period 2 pi freqs 1 itime call
55. ame9 Enter output file name write enter output file name in single quotes read fileout Read number of y direction points from surface dat read 10 ny read 10 y j j 1 ny write number of y points ny write maximum y y ny Read surface data do 10 i 1 100000 read 10 xold i if xold i 1t 0 go to 20 read 10 surfold i j j 1 ny continue continue m i 1 write number of x points in file m write maximum x xold m dy y 2 y 1 dx dy write grid spacing x and y in new image dy nx int xold m dx 1 write number of x points in interpolated image nx do 25 3 1 ny jj ny j 1 surface 1 3 surfold 1 3 3 continue x 1 0 do 40 i 2 nx 1 x 1 float i 1 dx 158 30 35 40 45 50 51 do 35 ii 1 m 1 if xold ii le x i and xold iitl gt x i then fac x i xold ii xold ii 1 xold ii do 30 3 1 ny jj ny 3 1 surface i j 1 fac surfold ii jj fac surfold ii 1 Jjj continue endif continue continue do 45 3 1 ny surface nx J surfold m 3 continue open 11 file fileout do 50 i 1 nx write 11 51 surface 1 3 J3J 1 ny continue close 11 stop format 501 10 4 end 159 7 4 refdifplot m refdifplot refdifplot m Script file to read in wave height wave angle water depth and surface data from refdifl output
56. an have many times the resolution of the reference grid The principal purpose of the subgrid is to provide enough computational points to the numerical model to preserve accuracy The user specifies the number of subgrid divisions in the y direction with the parameter nd If nd 1 then the subgrid spacing in the y direction is the same as the reference grid If nd 2 then the model uses twice as many computational points in the y direction as there are in the reference grid In the propagation direction x the model will automatically determine the subgrid spacing if ispace has been set to unity Otherwise the user provides the subgrid spacing using the input mr which permits variable spacing in the x direction For subgrids the input flag isp must be set to one and an array isd must be specified 2 6 4 Damping of Spurious Computational Modes When the large angle parabolic approximation is used as a basis for the computation of wave fields around islands the presence of wave breaking and resulting sharp lateral variations in wave height leads to the generation of high wavenumber spectral components in the computed complex amplitude A in the lateral y direction Kirby 1986a has shown that these components have propagation velocities which can become large in an unbounded fashion as a result they can propagate across the grid filling the computational 19 domain with high wavenumber noise Previous versions of REF DIF 1 have attem
57. an nr data values need to appear in the file wave dat An insufficent number of data points will simply trigger an end of file type of error during the read process The user may have any number of reasons for wanting to compute A externally to the program Several possibilities are e Running tests with complex initial conditions as in waves through a breakwater gap etc e Testing a directional distribution model which is different from the model chosen here e Specifying a nearly planar wave field which has height and angle variations along the y direction To serve as an example consider the case where you wish to specify a wave field having varying height 2a y and direction 0 y along an offshore boundary having nonuniform depth h y Assuming that the local wavenumber k y has been determined from the dispersion relation 2 T gk tanh kh 51 T where T is the wave period then the complex amplitude A can be specified by Aly a y e 52 where the phase function 4 y is given by v y i k y sin 0 y dy 53 35 3 7 Program Input Model Control and Wave Data This section discusses the structure of the input data file indat dat which is read in through logical device number iun 5 This file contains all the information needed to control the operations of the program and all of the input wave data Reference grid values of depth dr and currents ur and vr are treated in the following section Data is read fro
58. and Intp V Wave current velocity interpolated from the circulation module Intp eta Wave surface elevation interpolated from the circulation module Details can be found in the master program documentation Note that using the master program does not affect input and output of REF DIF 1 136 6 2 param h param h is used to define array lengths for all dimension parameters 6 3 common h common h includes all the common blocks used in the REF DIF 1 6 4 pass h pass h is used to define all of the passing variables used in the nearshore community model See details in the master program documentation 137 7 PROGRAMS FOR GENERATING AND POST PROCESSING DATA FILES This section provides listings for the following programs e indat createv26 f used to create the indat dat file based solely on input from the user e datgenv26 f used to create Version 2 5 indat dat and refdat dat files for specific examples e surface f used to convert the data in file surface dat to a regularly spaced ascii formatted array repre senting an instantaneous picture of the water surface e refdifplot m sample Matlab program illustrating the reading and plotting of the output data Note that the versions of the programs supplied with the program distribution may be slightly updated relative to the codes listed here 138 7 1 indat createv26 f indat create c c c c c c c c c c c c c c c c
59. aries and the external region An exception to this rule occurs when two subgrids adjoin each other Then care should be taken that the data along the common boundary match These common boundary point data are duplicated in the input data since each subgrid data set includes the boundaries 32 IR I JR 1 y IRSA ar Ren IR JR x O DR IR I JR MD IR 4 for example for example Figure 5 Interpolation of depth data 33 Reference Grid Sub grid Feature User Specified Group of Subgrid Blocks Grid Block IR JR IR JR 1 IR I JR l Sub grid Feature Depth Contours IR JR IR I UR Reference Grid Points MD IR 5 Figure 6 User defined subgrids 34 Instructions and formats for creating the subgrid data file subdat dat are included in section 2 8 3 6 User Specifi cation of Complex Amplitude on First Grid Row This section discusses the option of inputing the values of complex amplitude A for the first grid row from an external data file This option is invoked by setting the data value iinput 2 in line 8 of indat dat See the following section In the event that the option is chosen the data should be stored in a file wave dat The format for writing data to the file should be similar to complex a iy write a j Jj 1 n The data will be read in by the model subroutine The data is read in after grid subdivision in the y direction and hence n rather th
60. artificial island run is presented in two forms First Table 1 provides values of wave heights at the measurement locations indicated in Figure 12 ir or Juan 3 1 28 Table 1 Calculated wave heights at measurement locations In addition to this we have constructed contour plots of instantaneous surface elevation and wave height these are shown in Figures 13 and 14 respectively Contour elevations for wave height are in increments of 5 ft 56 1800F 1600F 1400F 1200F I 2000 kJ A 1800 i d O 1600 O G30 Q 1400 A P E ges 1200 o 0 gt 1000 800 1g 600 400 200 0 Figure 14 Artificial island example contours of wave height 58 4 2 Wave Focussing by a Submerged Shoal In this example we study the propagation of an initially plane wave over a submerged shoal resting on a plane beach This example has been chosen for several reasons 1
61. as been told to generate its own subgrid spacings ispace 0 possible incompatibility in any or all subgrid blocks Action None by program Should restart unit with correct ispace isp values Error occurs in inref While calculating its own subdivision spacings the model may try to put more division in a reference grid block than is allowed by by dimension ix If this occurs the program uses the maximum number 46 of subdivisions allowed ix 1 but prints a message indicating that the reference grid spacing is too large with respect to the waves being calculated This problem may be circumvented by increasing the size of ix in parameter statements Message model tried to put more spaces than allowed in grid block ir Action Program performs fixup and continues Model resolution and accuracy may be poor and a finer reference grid or increased value of ix in the parameter statements should be used Error occurs in grid 8 While using the Stokes wave form of the model ntype 2 the model may encounter large values of the Ursell number indicating that the water is too shallow for that model to be appropriate The cutoff point recognized by the program is 4 h kh 0 5 Message Warning Ursell number u encountered at grid location i j should be using Stokes Hedges model ntype 1 due to shallow water Action The program should be re run with the composite nonlinear model Error occurs in fdcalc 9 The New
62. ased on a Minimax principle which further extends the range of validity of the model equations The wave current version of the resulting model is included here as an option and may be chosen by modifying the choice of coefficients in the FDCALC see equations 12 and 13 2 1 2 Diffraction models In contrast to the mild slope model which is valid for varying bathymetry researchers in the area of wave diffraction were developing models for constant bottom applications For example Mei and Tuck 1980 developed a simple parabolic equation for wave diffraction and applied it to the diffraction of waves by a slender island Their equation is 0A i OA oS 3 Ox 2k Oy where A is a complex amplitude related to the water surface displacement by n Ae kx ot 4 Yue and Mei 1980 using a multiple scales approach developed a nonlinear form of this equation which accurately predicts the propagation of a third order Stokes wave A striking result of their numerical exper iments was the development of Mach stem reflection due to the reflection of obliquely incident waves from a breakwater This phenomenon is uniquely a nonlinear effect and not predictable from a linear refraction theory The parabolic model described below combines the essential features of the two approaches described above The variable depth features of the mild slope equation along with extensions to include effects of wave current interaction are retained but the mod
63. at since the model does not calculate reflected waves the predicted wave heights at points and 2 will be identical Therefore computations may be started arbitrarily close to the leading edge of the island base in the absence of any current field distorted by the island s presence 4 1 1 Setting up the Model First the island topography was established Note that the region of the island above the water line is not treated explicitly in the computations We therefore represented the island in the input data as a right circular cone with a peak height of 153 33 ft and a base radius of 400 ft The model will truncate the island in order to create a thin film over the exposed portions see Figure 11 32 Thin Film Region Created by Model Y M Figure 11 Representation of the island geometry in the program Next the reference grid spacing was chosen Since the physical region to be modelled is small we picked a reference grid with fine enough resolution so that no subdivision of the reference grid will be required Using a wave period of 10 seconds and a depth of 60 ft we use the dispersion relationship 4n 7 gk tanh kh 55 to calculate kh 97 giving a wave length of L 389ft L is just slightly less than rp 400 ft The grid spacings dxr and dyr 20 ft were chosen for the reference grid giving approximately 20 points per wavelength away from the island We use 100 x 100 storage locations for the refere
64. at user will be supplying ubgrids isp 1 am has been told to generate its own subgrid 120x General format format format format format format format end spacings ispace 0 1 possible incompatibility in any or all subgrid blocks format physical unit switch iu il input in mks units format physical unit switch iu il input in english units format 1 7 switches for dissipation terms d ne cl Pu turbulent boundary layer qt doy porous bottom q orf es DE laminar boundary layer format 20x Refraction Diffraction Model for 120x Weakly Nonlinear Surface Water Waves 120x REF DIF 1 Version 2 6 March 2002 120x Center for Applied Coastal Research 120x Department of Civil Engineering 120x University of Delaware 120x Newark Delaware 19716 110x James T Kirby Robert A Dalrymple and Fengyan Shi 120x Copyright C 2002 James T Kirby 120x REF DIF 1 comes with ABSOLUTELY NO WARRANTY 120x and is copyrighted under the provisions of the GNU Public License isp 0 no user defined subgrids isp 1 user defined subgrids to be read y direction subdivision according to nd i3 no current values read from input files current values read from data files closed reflective lateral boundaries open lateral boundaries icur
65. ation write to the Free Software Foundation we sometimes make exceptions for this Our decision will be guided by the two goals of preserving the free status of all derivatives of our free software and of promoting the sharing and reuse of software generally NO WARRANTY BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE THERE IS NO WARRANTY FOR THE PROGRAM TO THE EXTENT PERMITTED BY APPLICABLE LAW EXCEPT WHEN OTH ERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND OR OTHER PARTIES PRO VIDE THE PROGRAM AS IS WITHOUT WARRANTY OF ANY KIND EITHER EXPRESSED OR IMPLIED INCLUDING BUT NOT LIMITED TO THE IMPLIED WARRANTIES OF MER CHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU SHOULD THE PRO GRAM PROVE DEFECTIVE YOU ASSUME THE COST OF ALL NECESSARY SERVICING REPAIR OR CORRECTION IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER OR ANY OTHER PARTY WHO MAY MODIFY AND OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE BE LIABLE TO YOU FOR DAM AGES INCLUDING ANY GENERAL SPECIAL INCIDENTAL OR CONSEQUENTIAL DAM AGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF
66. ber of spaces each reference grid cell is divided into rather than the number of extra points being inserted Several restrictions are placed on the choice of nd and md s First the maximum dimensions of the subdivided grid cell is given by the parameters ix and iy This implies that any md can be at most ix 1 and nd can be at most iy 1 nr 1 The maximum number of added spaces may be increased by increasing ix and iy in the parameter statements Further the y subdivision specified by nd is applied uniformly along each grid row for the full extent to the reference grid No provision is made for variable grid spacing in the y direction Grid spacing in the x direction is arbitrary so md s may differ arbitrarily for each grid block The user must specify the single value of nd in the input data Two choices may be made regarding ma s however 1 The user may let the program calculate md The program proceeds by calculating an average wavenum ber along the reference grid row and uses this to estimate the wavelength in the grid block The program then chooses a subdivision so that at least 5 grid points per wavelength will occur in the grid block 30 Reference Grid kx D Grid Row Currently Being Used For Calculations Figure 4 Sample grid subdivision 31 The program checks to see that this desired number of subdivisions does not exceed the maximum If 1t does the program reduces the number to the maximum and prints a warn
67. buted under the terms of Sections 1 and 2 above on a medium customarily used for software interchange or c Accompany it with the information you received as to the offer to distribute corresponding source code This alternative is allowed only for noncommercial distribution and only if you received the program in object code or executable form with such an offer in accord with Subsection b above The source code for a work means the preferred form of the work for making modifications to it For an executable work complete source code means all the source code for all modules it contains plus any associated interface definition files plus the seripts used to control compilation and installation of the executable However as a special exception the source code distributed need not include anything that is normally distributed in either source or binary form with the major components compiler kernel and so on of the operating system on which the executable runs unless that component itself accompanies the executable If distribution of executable or object code is made by offering access to copy from a designated place then offering equivalent access to copy the source code from the same place counts as distribution of the source code even though third parties are not compelled to copy the source along with the object code 4 You may not copy modify sublicense or distribute the Program except as expressly provided under
68. c c indat createv26 f This program generates an indat dat file by asking the operator a series of questions This file is intended to make life a little easier its function is just as easily carried out manually if you are used to the form of the indat dat file James T Kirby Center for Applied Coastal Research University of Delaware Newark DE 19716 302 831 2438 FAX 302 831 1228 kirby udel edu Last revision 02 28 02 program indatcreate include param h dimension md ixr dconv 2 iff 3 dimension freqs ncomp edens ncomp nwavs ncomp dimension amp ncomp ncomp dir ncomp ncomp tide ncomp character 255 fnamel fname2 fname3 fname4 fname5 fname6 fname7 fnames fname9 fnamel0 fnamell fnamel2 fnamel3 fnamel4 fnamel5 fnamel6 fnamel7 fnamel8 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 namel indat dat fname2 refdat dat name3 subdat dat fname4 wave dat name5 refdifl log fname6 height dat name7 angle dat fname8 depth dat name9 fnamel0 sxx namell sxy dat fnamel2 syy dat namel3 fx dat fnamel4 fy dat namel5 qx dat fnamel6 qy dat namel7 tbx dat fnamel8 tby dat namel9 ibrk dat name20 owave dat name21 fname22 fname23 fname24 name25 fname26 data Fh Fh Fh Fh Fh Ph Fh rh Fh Fh Fh Fh th namelist ingrid mr nr iu ntype icur ib
69. c ismooth dxr dyr ispace nd iff isp iinput ioutput inmd md 139 fnames fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell fnamel2 fnamel13 fnamel4 fnamel5 fnamel6 fnamel7 fnamel8 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 wavesla iwave nfreqs waveslb freqs tide nwavs amp dir waveslc thet0 freqs tide edens nwavs nseed waves2 fregin tidein open unit 10 file fnamel write enter name for dat file containing reference grid in 1 single quotes read fname2 write enter name for ouput data file read fname2 write 10 nml fnames Enter control data write read write enter grid dimensions mr mr nr enter grid spacings dxr read wri read wri wri read te EE te dxr dyr dt input iu l mks iu input dispersion relationship nr dyr and depth tolerance dt 2 english ntype O linear ntype write input lateral write read ibc wri aid read write read T if ispace eq 0 write read if iansl write be l user chooses Ky 25 input ispace ispace input nd nd 4 y divisions go to 105 iansl eq 0 then input constant md 140 boundary condition l composite 2 stokes
70. calculate wave forcing c Model complete for the ifreq frequency component go to the next frequency C component 200 continue E Runs completed for all frequencies Return to end of main program return 201 format x f10 2 psibar f20 4 202 format 501 10 4 203 format 1 20x model execution frequency 1 component i14 end 98 5 5 GRID Interpolate the depth and current grids for reference grid block ir Velocities u and v are now set to zero in thin film areas Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby October 1984 August 1997 refdif subroutine grid ifreq ir include param h include common h include pass h integer i j c Constants pi 3 1415927 eps 1 0e 05 c Perform y interpolation on reference grid c Interpolate first row do 10 j 1 n nd d 1 j 2dr ir j 1 nd 1 Intp eta Wave ir j 1 nd 1 u 1 3j zur ir j 1 nd 1 v 1 3 vr ir j 1 nd 1 10 continue do 12 jj 22 nr do 11 j 1 nd 1 jjj nd 33 2 341 d 1 jjj dr ir jj dr ir jj 1 y 333 dyr 1 yr 33 dr ir 33 1 dew ea m Serdar dyr u 1 433 ur ir jj ur ir 35 1 y 333 dyr 1 yr 33 ur ir 33 1 aes aS mae ure dyr v 1 333 vr ir jj vr ir 35 1 y 333 dyr lal ye 45 orn o 3 1 Grid
71. controiling parameters zem Ti mh D m i m TA i mi El E ba El IL ni mt Ba xecute model reduencyY coma Figure 1 REF DIF 1 refdifl program level 25 Figure 2 REF DIF 1 model subroutine level 26 4 model Called from WaveModule model controls execution of the computational part of the program For each frequency component specified in the input model performs the following series of operations a Initialize program by calculating the incident wave field on the first grid row b Then for each grid block in the reference grid e Call grid to perform the grid interpolation specified in the input data e Call con to calculate constants on the interpolated grid e Call fdcalc to perform the numerical integration of the parabolic equation over the interpo lated subgrid The model execution is then complete Control is returned to refdif1 5 grid Called by model grid performs the required interpolation over a single grid block of the reference grid as specified in the input data The interpolation is performed as described in section 2 3 grid checks to see whether a user specified subgrid feature should be read in and reads it in from logical unit number iun 2 if needed The interpolated depth grid is then corrected for tidal offset and checked for surface piercing features These features are modified using the thin film approach see Kirby and Dalrymple 1986a Control is return
72. ction d i j d m j d 1 3 x i dxr x m d 1 3 x 1 d u i j u m j u 1 3 x i dxr x m u 1 j x 1 u v i j v m j v 1 3 dxr x m v 1 3 x 1 v 19 continue c Add in user specified grid subdivisions do 30 jr 1 nr 1 f isd ir jr eq 1 js nd jr 1 nd jf jstnd read 3 101 d i if icur eq 1 then read 3 101 u i read 3 101 v i endif do 31 i 1 m do 31 j23s jf d i j d i u i j u i v i j v i Bil continue end if 30 continue c Add tidal offset to all rows and establish thin film rj j js jf 3 dconv iu 3 dconv iu 3 dconv iu then i 1 m rj j js jf i 1 m rj j js jf i 1 m Set c current speed to zero in thin film area do 20 i 1 m do 20 3 1 n d i j d i if d i Tq E 0 001 d i j Jes Pc DEED v i j endif SSSL 00 0 0 20 continue c Interpolation complete return to J tide ifreg then model 101 read from unit 3 return 100 format model tried to put more spaces than allowed in 1 grid block i3 101 format 501f10 4 end 102 5 6 CON Subroutine calculates constants for reference grid block ir Program now checks for the existance of blocking currents and reduces the flow velocity to remove the blocking if it occurs Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE
73. cy of the wave the depth and the wave number is changed to reflect the Doppler shift due to currents The new form of eq 2 is w kU gk tanh kh 14 where the absolute frequency c is related to the intrinsic frequency c by w 0 kU 15 where the assumption that the wave is primarily travelling in the x direction has been used 11 2 2 Assumptions The REF DIF 1 model in parabolic form has a number of assumptions inherent in it and it is necessary to discuss these directly These assumptions are 1 Mild bottom slope The mathematical derivation of the model equations assumes that the variations in the bottom occur over distances which are long in comparison to a wave length For the linear model Booij 1983 performed a comparison between an accurate numerical model and the mild slope model for waves shoaling on a beach He found that for bottom slopes up to 1 3 the mild slope model was accurate and for steeper slopes it still predicted the trends of wave height changes and reflection coefficients correctly 2 Weak nonlinearity Strictly the model is based on a Stokes perturbation expansion and is therefore restricted to applications where Stokes waves are valid A measure of the nonlinearity is the Ursell parameter which is given as U HL ih 16 When this parameter exceeds 40 then the Stokes solution is no longer valid In order to provide a model which is valid in much shallower water a heuristic dispersi
74. d dep 0 116 dm 0 1379 dep e0 dep dm rad ak0 2 3 1415927 0 288 rad 152 301 302 c 303 sig2 9 806 ak0 tanh ak0 dep sig sqrt sig2 do 301 i 1 m x 1 float i 1 dx continue do 302 j 1 n y 3 float 3 1 dx continue do 303 i 1 m do 303 j 1 n r sqrt x i x if ix rad dx 1 2 y 3 y n t1 2 VARN if r gt rad d i j dep if r le rad d i j dmte0 r r if r le rad d i j dmte0 r continue endif if itype eq 6 then grazing incidence on caustic kirby and dalrymple 1983 351 352 353 write iun 3 input m n dx dy depth period read m n dx dy dep per pi 3 1415927 sig 2 pi per d2 2 dep alph atan 0 02 thet 25 pi 180 b d2 dep tan alph tt tan thet do 351 j 1 n y j float j 1 dy continue do 352 i 1 m x i float i 1 dx continue do 353 i 1 m do 353 j 1 n if y 3 1t y n x 1 tt d 1 3 dep 1f y 3 9e y n x 1 tt d 1 3 dep cos thet tan alph 1 x 1 tt y 3 y n 1f y 3 gt y n x i tt b cos thet d 1 3 d2 continue endif if itype eq 7 then whalin s channel 1971 401 write iun 3 401 format whalins channel input wave period read period write iun 3 period 153 402 403 404 405 m 100 n 74 iu 1 dx 242424242 dy 33866666 4 pi 3 1415927 sig 2 pi period do 402 i 1 m x i float 1 1
75. dal jet is shown in Figure 19 This photograph taken from a paper by Hales and Herbich 1972 is provided for guidance in interpreting the contour plot of surface elevations provided below 4 3 1 Setting Up the Model We choose a grid spacing of dxr 5m and dyr 5m We choose mr 100 and nr 100 giving an offshore and longshore extent of 495m The most shoreward grid row is established 5m from shore giving a depth range of 10m to 0 1m Arthur s wave period is retained and we use an initial wave amplitude of 0 1m The input conditions are a single normally incident wave no user specified grid subdivision and one frequency component 64 A H 233 3y0v MOUs JINVISIO IDEE PI e E A tT a oer 1334 11430 O NL RE QA eo Figure 18 Pattern of orthogonals and wave crests for waves in presence of rip currents refraction approxi mation from Arthur 1950 65 o_O o o Teoro T H LE G v 1 L60 O0 0 N uNsuno SEX NY OINO DNIIVOVdON NIVEL ZAVM TT ora tidal jet from Hales and Herbich 1972 Figure 19 Wave pattern on an ebb 66 The input data file for the present case follows Sfnames Same as previous example Send Singrid mr 100 nr 100 iu 1 ntype icur ibc dxr 5 000000 dyr 5 000000 dt 10 00000 ispace nd iff isp iinput ioutput 1 Send Swavesla iwave 1 nfreqs Send Swaveslb freqs 8 000000 tid
76. dif continue return format at grid block i3 location i i4 lfound a negative blocking velocity u f8 2 lu f8 2 end 104 j72 i4 model and reduced it to 5 7 FDCALC Perform the Crank Nicolson finite difference calculations on grid block ir Method is the implicit implicit iteration used by Kirby and Dalrymple 1983 Parameters for use in determining the minimax approximation are defined here 60 degree minimax coefficients a0 0 998214 al 0 854229 b1 0 383283 70 degree minimax coefficients a0 0 994733 al 0 890065 bl 0 451641 80 degree minimax coefficients a0 0 985273 al 0 925464 b1 0 550974 Pad coefficients default a0 1 al 0 75 bl 0 25 Small angle coefficients Radder s approximation a0 1 al 5 b1 0 0 Coded by James T Kirby October 1984 January 1992 July 1992 There is an unexplained and odd behavior in the minimax model when waves around islands are com puted For this reason the program distributed here has the coefficients for the Pad model refdif subroutine fdcalc ifreq ir 105 include param h include common h include pass h integer i j real kap ksthl ksth2 complex cl c2 c3 cpl cp2 cp3 ci damp complex ac iy bc iy cc iy rhs iy sol iy dimension thet iy urs iy height iy heightb iy dimension sxy iy sxx iy syy iy sbxx iy sbxy i
77. e 124 t 2 pi sqrt g k tanh k d k u write 5 100 i j k u d f t icdw 1 return 2 k kn return 100 format wavenumber iter failed to converge on row i10 1 column i10 T k f15 8 u f15 8 quf d f15 8 f 2 f15 8 1 t f15 8 end 125 5 15 Calculate wave forcing This subroutine is used to calculate radiation stresses according to the formulas used in Shorecirc model The roller effect is considered and results are transited with cubic SPLINE function to get smooth forcing around the breaking line Notice that the transition is only used in one direction on shore Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by Fengyan Shi 02 04 2002 refdif subroutine calculate_wave_forcing include param h include common h include pass h integer i j real Sm Nx Max Ny Max Sp Nx Max Ny Max SXxfake Nx_Max Sxy fake Nx_Max Syyfake Nx Max Qwfake Nx Max xrr Nx Max y2 Nx Max Xfake Nx Max sxx Nx Max Ny Max SXY Nx Max Ny Max syy Nx Max Ny Max qwtrs Nx Max Ny Max lhetar Nx Max Ny Max RR BS m m pi 3 1415926 grav 9 8 B0 1 8 rho 1030 0 period Pass_period c calculate Sxx Sxy Syy bes Outside the surfzone do j 1 Ny wave do i 1 nx wave Thetar i j pi 180 Pass Theta i Jj enddo enddo do j 1 ny wave do i 1 nx wave wave number Pass WaveNum
78. e 0 0000000E 00 nwavs 1 amp 0 5000000 dir 0 0000000E 00 Send H OO Oo G2 PES p o x o ll pa The input data files may be constructed for this test case using the program datgenv25 f 4 3 2 Model Results Results for this case are limited to plots of surface contours and wave height contours These plots were constructed using the information stored in the data files depth dat height dat surface dat The plots are given in Figures 20 and 21 Note that the plots only cover the region 51 ir lt 100 26 jr lt 75 in order to show greater detail in the wave pattern over the rip current The plots show a shoaling plane wave which approaches the beach with no distortion until the wave begins to interact with the rip current The rip causes a focussing of waves and the formation of discontinuities in the wave crests as in the photograph in Figure 19 67 350 I 4 300 TI E250 uz E Y O gt a 200 150 i j 250 300 450 500 Figure 20 Waves interacting with a rip current Shoreline at right Surface displacement contours 68 350r 300r E 250 200r 1505
79. e that contradict the conditions of this License they do not excuse you from the conditions of this License If you cannot distribute so as to satisfy simultaneously your obligations under this License and any other pertinent obligations then as a consequence you may not distribute the Program at all For example if a patent license would not permit royalty free redistribution of the Program by all those who receive copies directly or indirectly through you then the only way you could satisfy both it and this License would be to refrain entirely from distribution of the Program If any portion of this section is held invalid or unenforceable under any particular circumstance the balance of the section is intended to apply and the section as a whole is intended to apply in other circumstances It is not the purpose of this section to induce you to infringe any patents or other property right claims or to contest validity of any such claims this section has the sole purpose of protecting the integrity of the free software distribution system which is implemented by public license practices Many people have made generous contributions to the wide range of software distributed through that system in reliance on consistent application of that system it is up to the author donor to decide if he or she is willing to distribute software through any other system and a licensee cannot impose that choice This section is intended to make thoroughly c
80. e directional distribution is discretized into 31 compo nents each with an amplitude characteristic of the waves in that particular directional band These discrete waves are then assigned random phases and summed as in Eq 21 Note that this directional spectrum is for a given frequency and not for a continuous distribution of frequencies as in a true directional spectrum At the present time the REF DIF 1 model can only calculate waves at a single frequency per calculation The model will compute numerous frequencies per computer run set nfregs greater than one however they are not superimposable as the wave wave interactions between different frequencies are not included Using the linear mode ntype 0 and superimposing the results will provide a linear approximation to a directional spectrum The problem of computing the shoreward evolution of a directional spectrum of refracting diffracting and breaking waves is addressed in the model REF DIF S Chawla et al 1998 This model is an enhancement to REF DIF 1 which allows for the simultaneous computation of many wave components as a simulation of a random sea 2 5 Model Output 2 5 1 Complex Amplitude The complex amplitude A is output as a complex variable 2 5 2 Wave Heights and Angles Wave heights are calculated using H 2 4 Spatial smoothing of wave height is carried out if ismooth 0 arctan Ka kz k where k represents the weighted average of k along the trans
81. e lateral boundary conditions are important There are several ways to treat the boundaries however none of the presently existing boundary conditions result in the total transmission of scattered waves Therefore for the REF DIF 1 model a totally reflecting condition is generally used for each side j 1 and n This requires that the specification of the model grid be done with care as the reflection of the incident wave from the lateral boundaries can propagate into the region of interest rapidly and result in erroneous results In general the width of the model should be such that no reflection occurs until far downwave of the re gion of interest As a precaution a graphical representation of the computed wave field should be examined to determine where the reflection from the boundaries is important By using the switch ibc partially transmit ting boundaries can be used Kirby 1986c In general this boundary condition will result in less reflection in the model domain however since some reflection will occur it is recommended that runs be carried out with the reflecting boundary conditions in order to assess the regions potentially affected by reflection from the model boundaries 2 6 3 Subgrids In order to reduce the amount of data input and yet provide the user the ability to prescribe the fine scale bathymetry in areas of interest REF DIF 1 utilizes a coarse scale user specified reference grid and a fine scale subgrid which c
82. e pattern on an ebb tidal jet from Hales and Herbich 1972 Waves interacting with a rip current Shoreline at right Surface displacement contours Waves interacting with a rip current Shoreline at right Wave height contours 1 INTRODUCTION REF DIF 1 is presently used by hundreds of researchers practicing engineers and planners worldwide The program is freely distributed through the web site http chinacat coastal udel edu kirby programs refdif refdif html and links to various activities using the program are provided The program is provided without warranty and under the copyright model of the Free Software Foundation as detailed below Work on the present upgrade of REF DIF 1 is supported by the National Ocean Partnership Program NOPP through the project Development and Verification of a Comprehensive Community Model for Physi cal Processes in the Nearshore Ocean described at http chinacat coastal udel edu kirby NOPP index html The main goal in the upgrade to Version 3 0 has been to provide compatibility between REF DIF 1 and the Nearshore Community Model system Similar compatibility is being provided for the spectral model REF DIF S Chawla et al 1998 and a time dependent refraction diffraction model developed by Kennedy and Kirby 2002 Each of these models will be documented independently and will be provided as free standing programs and as Nearshore Community Model components 1 1
83. ear dispersion terms in the governing equations 2 1 4 Wave current interaction models Booij 1981 using a Lagrangian approach developed a version of the mild slope equation including the influence of current This model is a weak current model in that the currents are assumed to be small and any products of currents are neglected as small Kirby 1984 presented the corrected form of this mild slope model A nonlinear correction and the ability to handle strong currents were added by Kirby and Dalrymple 1983b and results for waves interacting with a current jet were obtained Their equation is Cg U A VAy i k k l Cy U A 4 e S E a oO i k 97 9 VA o DIAI A 0 5 where p CC and k reference wave number taken as the average wave number along the y axis and U is the mean current velocity in the z coordinate direction and V is in the y direction The nonlinear term includes D which is cosh 4kh 4 8 2 tanh kh 8sinh kh Kirby 19862 rederived the above equation for a wide angle parabolic approximation which allows the study of waves with larger angles of wave incidence with respect to the x axis This more accurate equation was used as the basis for earlier versions of REF DIF 1 The equation has been extended to include the more accurate minimax approximation Kirby 1986b for the present version of REF DIF 1 The revised governing equation is given by Cy U Az 2A1V Ay i k
84. ed to model 6 con Called by model con calculates various constants for the reference grid created by grid Control Is returned to model 7 fdcalc Called from model fdcalc performs the integration of the governing parabolic equation over the grid defined in grid The coefficients of the finite difference form of the parabolic equation are developed according to the Crank Nicolson method A complete description of the equations and the treatment of nonlinearities may be found in Kirby 1986a and Kirby and Dalrymple 1986b The sequence of steps in fdcalc is as follows a An implicit step is performed to update complex amplitude A along an entire grid row b The model checks for the start or stop of breaking on the updated row c If the status of breaking changes the model recomputes the breaking wave dissipation coefficient d Then if nonlinearity is being used or breaking status at any point along the row has changed the model computes a new estimate of A on the updated row based on values obtained during the previous iteration This series of operations is performed for each row in the subdivided ir grid block until the end of the grid defined in grid is reached Control is then returned to model which passes on to the next ir 1 reference grid block 27 8 ctrida Utility routine which is called by fdcalc to perform the double sweep elimination to solve the implicit set of equations 9 diss Called by co
85. een 1 and n The procedure involves expressing all the derivatives in the x y directions in terms of the complex amplitude at various grid points For example the forward difference representation of dena Ay at location 2 7 49 If a forward difference is used for the x direction and a central difference representation centered at 2 is used for the second derivatives in the lateral direction for all the derivatives in Eq 6 then an explicit finite difference equation results for A 1 This equation can be solved directly for all the A j 1 2 3 n for a given 2 provided appropriate lateral boundary conditions are prescribed This explicit representation is not as accurate as an implicit scheme and therefore an implicit Crank Nicolson procedure is used for the amplitude calculations For a given i row the Crank Nicolson scheme can be written Qi 1k bAi41 C i L1 j 1 dAj 541 eAi fAi j 1 50 where the coefficients a b c d e f involve variable complex and nonlinear terms The amplitudes on the left hand side of this equation are unknown while the terms on the right hand side are known from either the previous calculation or from the initial boundary condition on 1 and n This equation is solved for all the A 1 j 2n 1 and i fixed at once by a tridiagonal matrix solution procedure Carnahan Luther and Wilkes 1969 adapted for complex valued coefficients Due to the nonlinearity of the finite
86. efraction Diffraction Model REF DIF 1 Version3 0 72 3 2 INRER i148 ci ad br ae se sugueimi TI 5 2 1 Read file names from namelist eee 81 5 35 INWAVE cade S04 ad e Me ooh ee ek ee a be p Royan E ee pu eS 87 5 4 MODEL see o a oe GL uuum iege eed he ts hs 92 I GRID eek wee ee ee be ee ee ee oe ee e 99 3 07 CONE 3 sce noe shod be Oy eS Se Shae ee RUP a PA ee da Bie dee oe ee Ed 103 IAN EDCALE aA cect ta le E EAS RAR Fes s RUE pO URN S 105 5 7 FDCALC statement functions en 107 5 8 SCTRIDA onu Soe oe a A och eee RUP ep PA ee ae BR dee oe te eo 119 DIDIER 120 LO RAND i we doc ad a ee ek E alada ebd 121 SJHLRDEACT cacc s oi a ee E ep e Bie dee pum Hes 122 2 12 BNUML ara ta dme Te EU VCI de Ro ee etos UO Te Bet UPC ER 122 SJ3 ACAEG elm IE eges UEM EUR PEE Faccio o ae Bae ad ae a 123 5 14 WN NUM eu ah oS a ee ae om Rd a RUE p e que ee Es 124 DAS Calculate Wave TOFCIDOS 2 4 4 vere RE RA But s EE PO SEP OOM EOS 126 5416 SPLINE zu os e GREGIS Sig A A Ru EE e RU 134 DE SBEINIEE Sae Aaa e fede EUR eR A A osos ae eee is HE za 135 6 ADDITIONAL MODEL COMPONENTS 136 Oil MASter actos E eve ele A a mein 136 6 2 parame 4 1s s ee AA o a ey a Oe 137 6 3 COMMON Mam ass eee ee a ee eA ee eel eee ee we Se e seh x 137 64 DASS A BO ae A pump BA OR Ua Sur AS 137 7 PROGRAMS FOR GENERATING AND POST PROCESSING DATA FILES 138 TA INndat Cregiev 207 safe a v ke proe a euet hel ne 139 7 2 dat
87. el is developed in parabolic form and in terms of a complex amplitude A 2 1 3 Nonlinear combined refraction diffraction models Kirby 1983 using a Lagrangian approach and Kirby and Dalrymple 19832 with a multiple scales tech nique developed the predecessor to the REF DIF 1 model which bridged the gap between the nonlinear diffraction models and the linear mild slope equation This model can be written in several forms depending on the application The hyperbolic form for time dependent applications and the elliptic form for steady state problems require the use of boundary conditions on all sides of the model domain This is a difficult requirement as the reflected wave at a boundary is not generally known a priori These models however have the advantage that there is no restriction on the wave direction A detailed comparison of results of the weakly nonlinear model of Kirby and Dalrymple 1983a to labo ratory data was shown by Kirby and Dalrymple 1984 The laboratory test conducted at the Delft Hydraulics Laboratory by Berkhoff Booij and Radder 1982 consisted of determining the wave amplitude over a shoal on a sloping bottom While results predicted by ray tracing techniques were shown by Berkhoff Booij and Radder to be very poor the agreement between the weakly nonlinear model and the laboratory data was ex cellent Comparisons between linear and nonlinear parabolic models clearly showed the importance of the nonlin
88. el to the y axis the wave is generally prescribed as A 0 y Age Q4 where Ag is the given wave amplitude and is the wave number in the y direction The is related to the wave number k by the relationship ksin 0 where 0 is the angle made by the wave to the x axis This case is obtained by using the data switches wave and nwavs set to one 2 4 2 Discrete directional waves not presently recommended For several waves with different directions at a given frequency the following relationship could be used for the initial wave condition nwavs A 0 y Y Ane 25 n l The REF DIF 1 model is equipped to calculate the wave field produced by this boundary condition for a large number of user supplied A s and 0 s up to 50 This mode is accessed by iwave 1 and nwavs set to the number of discrete waves to be used 2 4 3 Directional spectrum not presently recommended Often a cos 0 directional spread is used with a given frequency component This can be done with REF DIF 1 by specifying iwave equal to 2 and nwavs to the value of n The total energy at frequency c is 27 0 0 Elo En cos 0 26 In order to avoid the problem of waves propagating at large angles to the propagation direction 09 the directional distribution of energy is automatically truncated to include only those directions which contain more than 10 of the total energy To prescribe the initial conditions for the model th
89. ence grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel6 qy dat short wave flux in y direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated 77 e fnamel7 bottomu dat magnitude of bottom velocity at reference grid points If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel8 not used at present e fmamel9 ibrk dat wave breaking index If REF DIF 1 is given a null string as the input for this file name no file is generated e fname20 owave dat complex amplitude on last row for ioutput 2 If REF DIF 1 is given a null string as the input for this file name no file is generated e fname21 sxxb dat body stress part of local radiation stresses If REF DIF 1 is given a null string as the input for this file name no file is generated e fname22 sxyb dat body stress part of local radiation stresses Sz If REF DIF 1 is given a null string as the input for this file name no file is generated e fname23 syyb dat body stress part of local radiation stresses Syy If REF DIF 1 is given a null string as the input for this file name no file is generated e fname24 sxxs dat surface stress part of local radiation stresses 54 If REF DIF 1 is given a null string as the input for this file name no file is generated e fname25
90. enter for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby Oct 1984 Sept 1989 Jan 1991 July 1994 November 1994 February 2002 refdif subroutine inwave include param h include inc namelist wavesla iwave common h lude pass h nfreqs 87 waveslb update_interval num_data fregs tide nwavs amp dir waveslc thet0 freqs tide edens nwavs nseed waves2 fregin tidein poumon real dis numdata integer i j pi 3 1415927 c Values of iinput ioutput already entered in namelist statement in c inref if iinput ne 1 and iinput ne 2 then write 5 invalid value chosen for iinput check indat dat stop endif if ioutput ne 1 and ioutput ne 2 then write 5 invalid value chosen for ioutput check indat dat stop endif if iinput eq 1 then write 5 iinput 1 program specifies initial row of a c Enter iwave nfreqs for iinput 1 read 1 nml wavesla write 5 102 c Enter data for case of iinput 1 iwave 1 if iwave eq 1 then read 1 nml waveslb write 5 103 endif c Enter data for case of iinput 1 iwave 2 if iwave eq 2 then read 1 nml 2waveslc write 5 104 endif write 5 105 nfregs if iwave eq 2 then 88 thet0 thet0 pi 180 endif For each frequency enter the wave period and tidal offs
91. ers to any such program or work and a work based on the Program means either the Program or any derivative work under copyright law that is to say a work containing the Program or a por tion of it either verbatim or with modifications and or translated into another language Hereinafter translation is included without limitation in the term modification Each licensee is addressed as gt you Activities other than copying distribution and modification are not covered by this License they are outside its scope The act of running the Program is not restricted and the output from the Program is covered only if its contents constitute a work based on the Program independent of having been made by running the Program Whether that is true depends on what the Program does You may copy and distribute verbatim copies of the Program s source code as you receive it in any medium provided that you conspicuously and appropriately publish on each copy an appropriate copy right notice and disclaimer of warranty keep intact all the notices that refer to this License and to the absence of any warranty and give any other recipients of the Program a copy of this License along with the Program You may charge a fee for the physical act of transferring a copy and you may at your option offer warranty protection in exchange for a fee You may modify your copy or copies of the Program or any portion of it thus forming a
92. et do 3 ifreq 1 nfreqs do itime 1 num_data write 5 107 ifreq fregs ifreg itime tide ifreq Convert angles to radians If If If fregs ifreq itime 2 pi freqs ifreq itime tide ifreq tide ifreq dconv iu iwave 1 read the number of discrete components if iwave eq 1 then do 1 iwavs 1 nwavs ifreq write 5 106 iwavs amp ifreq itime dir ifreg itime dir ifreq itime dir ifreq itime pi 180 amp ifreq itime amp ifreq itime dconv iu continue endif enddo iwave 2 read the parameters for each frequency if iwave eq 2 then seed float nseed 9999 write 5 108 edens ifreq nwavs ifreq nseed dir ifreg 1 thet0 edens ifreq edens ifreq dconv iu 2 endif continue endif iinput 2 read in wave period and tidal offset if linput eq 2 then read 1 nml waves2 freqs 1 1 fregin tide 1 tidein 89 write 5 iinput 2 user specifies a in wave dat nfreqs 1 write 5 102 write 5 wave period freqs 1 1 sec write 5 tidal offset tide 1 freqs 1 1 2 pi freqs 1 1 tide 1 tide 1 dconv iu endif C re arrange freqs amp and angle num total int total time Master dt N Interval CallWave print total time Master dt N interval CallWave print total number of wave updated num total if num total ge numdata then print you should make a larger numdata in param h endif
93. file fname2 O EEES skip above after the first call wave module open unit 5 file fname5 c Open optional output files c Store wave height if fname6 ne open 26 file fname6 c Store wave angle if fname7 ne open 7 file fname7 c Store water depth if fname8 ne open 8 file fname8 c Store water surface if fname9 ne open 9 file fname9 c Store depth integrated radiation stresses if fnamelO ne then open 10 file fnamel0 open 11 file fnamell open 12 file fnamel2 endif c Store depth integrated forcing if fnamel3 ne then open 13 file fnamel3 open 14 file fnamel4 endif c Store total wave induced mass flux if fnamel5 ne then open 15 file fnamel5 open 16 file fnamel6 endif 81 c c Store bottom velocity moments if fnamel7 ne then open 17 file fnamel7 open 18 file fnamel8 endif Store breaking index if fnamel9 ne open 19 file fnamel9 Store complex amplitude on last row if fname20 ne open 20 file fname20 Store radiation stresses for 3d circulation model Sum es A X Ru fname21 ne fname22 ne if open file fname21 3 if fname23 ne Tft if ift 21 open 22 file fname22 open 23 file fname23 24 25 27 fname24 ne fname25 ne fname26 ne open file fname24 open file fname25 open file fname26 print headers
94. for start or stop of breaking in each row do 6 j 1 n urs j 2 cabs a ipl j 6 continue isavel 0 isave2 0 do 7 3 1 n iset 0 ireset 0 if urs 3 d ip1 3 9t kap and ibr 3 eq 0 iset 1 if iset eq 1 then ibr 3 1 isavel 1 end if if urs 3 d ip1 3 1t gam and ibr 3 eq 1 ireset 1 if ireset eq 1 then ibr 3 0 isave2 1 end if 111 continue ih 2 c Redo initial calculation if breaking status changes if isavel eq 1 or isave2 eq 1 go to 2 continue if it eq 2 go to 9 it 2 go to 2 continue c For Stokes model alone ntype eq 2 test to s whether Ursell C parameter is too large 11 if ntype eq 2 then do 11 j 1 n urs 3 cabs a ip1 3 d ip1 3 k ip1 3 d ip1 3 2 if urs 3 9t 0 5 write 5 204 urs 3 1 continue end if c Roll back breaking dissipation coefficient at each row 12 do 12 j 1 n wb 1 3 wb 2 3 continue c Calculate reference phase function for surface plotting psibar psibar kb ip1 kb 1 dx 2 c Store plotted surface if requested if fname9 ne then write 9 x ipl write 9 real a ipl j cexp cmplx 0 psibar j71 n endif c Start filter if breaking is occuring do 13 j 1 n if ibr j eq 1 ifilt 1 13 continue 200 continue c Calculate wave angles at reference grid rows Note angles are not well c defined in a directional multicomponent sea or where waves become short
95. g in Version 3 0 eA 23 3 2 Overview of Operating Manudl eA 24 3 3 Program Outline and Flow Chart o a 24 3 4 Special Installation Instructions 000000000220 28 3 5 Computational Grids and Grid Interpolation a 28 3 5 1 Grid Subdivision 2 scu e E A ea a e 30 3 5 2 User specifled S bgrids so e RAR ER E Rer A IN 32 3 6 User Specification of Complex Amplitude on First Grid Row 35 3 7 Program Input Model Control and Wave Data lee 36 3 8 Program Input Reference Grid and Subgrid Data o 42 3 9 Program p Output ia a a a a RE Bag 43 39 1 Output log file 2n ines A A A eed A deus 43 3 92 Stored Outputs 24s st Pe Sat pe nk Rei eR OR a Pe ee Dee 48 4 EXAMPLE CALCULATIONS 50 4 1 Waves Around an Artificial Island lee 51 4 Ll Setting up the Model ia oo ra e dde Be ee eek dk 52 4 1 2 The InputData Files e o 55 4 1 3 Model R sults 20 20 000400 a a GR ts o 56 4 2 Wave Focussing by a Submerged Shoal o llle 59 4 21 TheInputData Files o e o 61 42 2 Model Results s v xus PAE RR A Bee ORS ien 62 4 3 Waves Interacting with a Rip Current e 64 4 3 1 Setting Up the Model e 64 43 23 Model Results estat os ar e e 67 4 4 Obliquely Incident Waves on a Plane Beach o e 70 111 5 REF DIF 1 Program Listing 71 5 1 R
96. genv20 fr es 2g p Ep SESS PS xU xS SE bee ER ESSE 144 T3 Ssurf c f svi scu Rb eere AE E AAA A quee iude eus xe 157 TAE vefdufplotgn u sm sme eR Eri debe SES SE rs delata 160 8 FREQUENTLY ASKED QUESTIONS 162 9 REFERENCES 164 List of Figures 0 ON QN ta Bb un REP DIF 1 refdifl programlevel o e REF DIF 1 model subroutine level 2 ee Reference grid notation o ooa Sample orid subdivision me v vod A A E Interpolation of depth data e User defined subgrids en Sample title papel iussa o RA A RS A PE ok Sample title pageZ sos or c RAO Be e SO Pte a Artificial island geometry ee Locations for wave height measurements 0 00000 0004 Representation of the island geometry in the prograM lll ln Measurement points in relation to reference grid o o 000002 ee Artificial island example contours of instantaneous surface elevation Artificial island example contours of wave height len Bottom contours and computational domain for the experiment of Berkhoff et al 1982 Experimental data on transects 1 8 een Results for waves propagating over a submerged shoal surface elevation contours Results for waves propagating over a submerged shoal wave height contours Pattern of orthogonals and wave crests for waves in presence of rip currents refraction ap proximation from Arthur 1950 2 2 ee Wav
97. gle if fname7 ne close 7 Water depth if fname8 ne close 8 74 c Water surface if fname9 ne close 9 c Radiation stresses if fnamelO ne then c cl cl lose 10 lose 11 lose 12 endif c Depth integrated forcing if fnamel3 ne then El c Lose 13 lose 14 endif c Total wave induced mass flux if fnamel5 ne then close 15 close 16 endif c Bottom velocity moments if fnamel7 ne then close 17 close 18 endif c Breaking index if fnamel9 ne close 19 c Stored complex amplitude on last row if fname20 ne close 20 c local radiation stresses if fname21 ne close 21 if fname22 ne close 22 if fname23 ne close 23 if fname24 ne close 24 if fname25 ne close 25 if fname26 ne close 27 C All done return 75 end 76 5 2 INREF Subroutine reads in and checks dimensions and values for large scale reference grid Wave parameters for the particular run are read in later by subroutine inwave The following unit device numbers are assumed e fnamel e fname2 e fname3 e fname4 e fname5 e fname 6 e fname7 e fname e fname9 indat dat automatically assigned to the namelist input data file indat dat refdat dat depth and u v on the reference grid subdat dat user specified subgrids wave dat user specified complex amplitude on ro
98. h Annalen 47 317 374 U S Army Coastal Engineering Research Center 1973 Shore Protection Manual Vol I Wiegel R L 1962 Diffraction of waves by a semi infinite breakwater J Hydraulic Div 88 27 44 Yue D K P and C C Mei 1980 Forward diffraction of Stokes waves by a thin wedge J Fluid Mech 99 33 52 166
99. he rest of indat dat is as follows wavesla namelist group e iwave iwave is switch for wave field type iwave 1 discrete wave components iwave 2 directional spread ing model not presently recommended e nfreqs nfreqs is the number of frequency components to be run Maximum is value of ncomp in param h The remainder of the file depends on the choice of iwave For iwave 1 waves1b namelist group 39 e freqs ncomp Wave period for each frequency component ncomp values must be given e tide ncomp Tidal offset for each frequency component ncomp values must be given e nwavs ncomp Number of wave components for each frequency ncomp values must be given e amp ncomp ncomp Amplitude not height for each component wave e dir ncomp ncomp Direction in degrees relative to x axis for each wave component For iwave 2 waveslc namelist group e 1het0 The central direction for the model spectrum e freqs ncomp tide ncomp As above e edens ncomp Variance density m or ft for each frequency component e nwavs ncomp Directional spreading factor the factor n in cos 0 2 e nseed The seed value for the random number generator between 0 and 9999 For iinput 2 the remainder of the data file is waves2 namelist group e freqin tidein Wave period and tidal offset for the single frequency component Examples of indat dat data files are given in the example
100. his program if not write to the Free Software Foundation Inc 59 Temple Place Suite 330 Boston MA 02111 1307 USA or get it from http www gnu org copyleft gpl html Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 refdif subroutine WaveModule c Define array dimension bounds include param h c Internal common blocks include common h c Common block data passing to Master program include pass h integer i j c master start 0 or 1 for initialization 73 Q if Master_Start eq 1 then write wave module initialization else write call wave module endif Constants which provide for conversion between MKS and English units on input and output dconv 1 1 dconv 2 0 30488 dconv2 1 1 dconv2 2 14 594 Read control parameters and reference grid data call inref Read control parameters and initializing wave data if Master_Start ge 0 then call inwave close 1 endif Pass program control to subroutine model For each frequency component specified in inwave model xecutes the model throughout the entire grid and then reinitializes the model for the next frequency if Master Start le 0 then call model endif All done Close output data files if open and close statements are being used Log file close 5 Wave height if fname6 ne close 26 Wave an
101. hoses to specify the md ir values himself the switch space has to be set to unity In this case the user again has a choice The user can specify isp 0 then the program will perform interpolations The user can set isp 1 in which case he she has to specify a subgrid in the file subdat dat 4 What are the guidelines for the specification of the number of subdivisions md ir in the x direction It is recommended that the number of subdivisions md ir be determined by the program However if the user wants to specify md ir s the user should make sure that the subdivided grid is at least as fine as the program would have determined In order to achieve this the user can first run the program by letting it pick its own subdivisions Then the user can chose his her subdivision using the programs subdivisions as a guideline 5 What are the guidelines for the choice of the number of subdivisions nd in the y direction nd should be chosen such that the final subdivided reference grid cells have increments in the z and y directions that are fairly close in size This is best determined by running REF DIF 1 once to check how many x direction subdivisions are introduced in the most finely subdivided grid and then subdividing in y by a corresponding amount 162 6 Where is the origin of the grid The origin of the domain is always chosen to be the lower right hand corner of the domain with the x axis pointing in the propagation
102. ibc ismooth dxr dyr dt ispace nd iff isp iinput ioutput inmd md fnames fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell fnamel2 fnamel13 fnamel4 fnamel5 fnamel6 fnamel7 fnamel18 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 wavesla iwave nfreqs waveslb freqs tide nwavs amp dir waveslc thet0 freqs tide edens nwavs nseed waves2 fregin tidein setup the logical devices for input keyboard input iun 2 output file refdat dat iun 3 screen output use 0 for sun 3 for pc iun 4 output file indat dat iun 2 20 iun 3 0 iun 4 24 fnamel indat dat open iun 4 file fnamel initialize all entries for indat dat prior to generating the depth grid iu 0 ntype 0 icur 0 ibc 0 ispace 0 nd 1 if1 0 if2 0 if3 0 isp 0 iinput 0 iwave 0 nfreqs 0 nwaves 0 145 c c open file for reference grid data c open iun 2 file fname2 establish depth grid call depth iun c c calculate constants c dt 10 call con write reference grid data 100 do 1 i 1 mr write iun 2 100 dr i j j 1 nr continue if icur eq 1 then do 2 i 1 mr write iun 2 100 ur i j j71 nr continue do 3 i 1 mr write iun 2 100 vr 1 3 3 1 n1r continue endif close iun 2 format 501f10 4 generation of file indat dat write d
103. ii iipl ii 1 do 1 i iip1 1 beta 1 b 1 a 1 c i 1 beta i 1 gamma i d 1 a 1 gamma 1 1 beta i 1 continue c Compute solution vector v v 1 gamma 1 last 1 11 do 2 k 1 last i 1 k v 1 gamma 1 c 1 v 1 1 beta i 2 continue return end 119 5 9 DISS Subroutine calculates the dissipation at a single grid point based on values of the switch iw at that point Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby October 1984 refdif subroutine diss include param h include common h real nu cp kd c Statement function Sq i j sqrt nu 2 sig i j c Constants nu 1 3e 06 cp 4 5e 11 g 9 80621 pi 3 1415927 c Value of f here is value assuming ltau f 8 u 2 c Sf 4 f w Sf wS is the wave friction factor f 0 01 4 0 do 1 j 1 n do 1 i 1 m w i j cmplx 0 0 kd k i j d i j e If iff 1 1 use turbulent boundary layer damping if iff 1 eq 1 w 1 3 2 f cabs a 1 J sig i j k i 3 1 sinh 2 kd sinh kd 3 pi c If iff 2 1 add porous bottom damping if iff 2 eq 1 w i j w i j g k i j cp nu cosh kd 2 1 cmplx 1 0 120 G If iff 3 1 add boundary layer damping if iff 3 eq 1 w i j 2w i j 42 k i j sig i j sq i 3 1 1 cosh kd 2 cmp1x 1 1 sinh 2 kd
104. ination of the amount of wave energy penetrating an island chain or calculation of the sheltering and hence the disturbance of the littoral processes by an island situated near a shoreline They are not intended to replace diffraction theories currently in use for wave force calculations The weakly nonlinear combined refraction and diffraction model described here denoted REF DIF 1 is based on a Stokes expansion of the water wave problem and includes the third order correction to the wave phase speed The wave height is known to second order Liu and Tsay 1984 It should be noted that it is not a complete third order theory as all the third order terms are not retained Known ambient currents which effect the height and direction of wave propagation are input for the model and enable it to predict waves where currents may be strong The application of the theoretical model to practical situations involves the use of a parabolic approxima tion which restricts the model to cases where the wave propagation direction is within 60 of the assumed wave direction and the use of finite difference techniques for the wave amplitude which results in tridiagonal matrices which are computationally very fast to invert The REF DIF 1 model is described in detail in this manual which also documents the application of the model to actual examples and provides explicit descriptions of the input and output 2 1 Wave Models 2 1 1 Mild slope equatio
105. index 3 dx jJump int 0 25 1 slope 1 1 Depth wave index j dx outside the surfzone do i 1 index Sxxfake i Pass Sxx i j Sxyfake i Pass Sxy i j Syyfake 1 Pass_Syy 1 J Owfake i Pass_MassFlux i Jj xfake 1 1 1 dx enddo inside the surfzone do i index jump nx wave Sxxfake i jump l Pass Sxx i Jj Sxyfake 1 Jjump 1 Pass_Sxy 1 J Syyfake 1 Jjump 1 Pass_Syy i j Qwfake i Jump 1 Pass_MassFlux i Jj xfake i jumpt1 i 1 dx enddo Do the cubic spline to for roller transition do i 1 nx_wave xrr i i 1 dx enddo call splinel nx wave ny wave xfake Sxxfake nx wave jump y2 call splintl nx wave ny wave xfake Sxxfake y2 nx wave jump xrr amp sxx index jump jJ call splinel nx wave ny wave xfake Sxyfake nx wave jump y2 call splintl nx wave ny wave xfake Sxyfake y2 nx wave jump xrr amp sxy index jump Jj call splinel nx wave ny wave xfake Syyfake nx wave jump y2 call splintl nx wave ny wave xfake Syyfake y2 nx wave jump xrr amp syy index jump j call splinel nx wave ny wave xfake Qwfake nx wave jump y2 call splintl nx wave ny wave xfake Qwfake y2 nx wave jump xrr amp qwtrs index jump j enddo do j 1 ny wave do i 1 nx wave Pass Sxx i j sxx i j Pass_Sxy i j sxy i j 128 Pass_Syy i J syy i j Pass_MassFlux i j qwtrs i j Pass_MassF1luxU i j Pass_MassFlu
106. ing message indicating that the grid block can not be subdivided finely enough Computed results in this case must be regarded as being suspect The number of subdivisions used is indicated on the output The user chooses this option by setting the switch ispace 0 on input 2 The user may specify each value of md ir from ir 1 to mr 1 This is done by setting the switch ispace 1 on input and the values of md are then read in from the input data file Note that this choice is necessary if the user defined subgrids discussed below are to be used since the user will be inputting subgrids with pre specified spacings As it is presently written the program will only print a warning if it encounters subgrids with ispace 0 chosen extensive garbling of input data may result After the subdivided grid block is established the program uses this grid as the actual computation grid New values of depth d and ambient current u and v are calculated at the extra grid points by fitting a twisted surface to the reference grid using linear interpolation An example of the resulting bottom topography for a single grid cell is shown in Figure 5 3 5 2 User specified Subgrids In some applications an important topographic feature may be present at a subgrid scale within the reference grid Examples include artificial islands shoals borrow pits etc which are superimposed on an otherwise slowly varying topography which is represented by the sort of grid reso
107. ints If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel8 not used at present e fmamel9 ibrk dat wave breaking index If REF DIF 1 is given a null string as the input for this file name no file is generated e fname20 owave dat complex amplitude on last row for ioutput 2 If REF DIF 1 is given a null string as the input for this file name no file is generated e fname21 sxxb dat body stress part of local radiation stresses Sz If REF DIF 1 is given a null string as the input for this file name no file is generated e fname22 sxyb dat body stress part of local radiation stresses Syy If REF DIF 1 is given a null string as the input for this file name no file is generated e fname23 syyb dat body stress part of local radiation stresses Syy If REF DIF 1 is given a null string as the input for this file name no file is generated e fname24 sxxs dat surface stress part of local radiation stresses Sz If REF DIF 1 is given a null string as the input for this file name no file is generated e fname25 sxys dat surface stress part of local radiation stresses Syy The values are always zero If REF DIF 1 is given a null string as the input for this file name no file is generated e fname26 syys dat surface stress part of local radiation stresses Syy If REF DIF 1 is given a null string as the input for this file name no file is generated If iinput 1 t
108. is exact Numerous authors have applied the mild slope model to various examples primarily using finite element techniques See for example Jonsson and Skovgaard 1979 Bettess and Zienkiewicz 1977 and Houston 1981 For the linear mild slope equation Radder 1979 developed a parabolic model which had several ad vantages over the elliptic form presented by Berkhoff First the boundary condition at the downwave end of the model area is no longer necessary and secondly very efficient solution techniques are available for the finite difference form of the model Radder used a splitting matrix approach which involves separating the wave field into a forward propagating wave and a backward propagating wave and then neglecting the backward scattered wave which is justified in most applications as only the forward propagating wave is used for design Radder s approximation for derivatives transverse to the wave direction results in a restriction on his parabolic model the waves must propagate within 45 of the assumed wave direction Booij 1981 also developed a splitting of the elliptic equation but his procedure included more terms in the approximation to the lateral derivative and therefore his procedure enables the parabolic model to handle waves propagating within 60 of the assumed direction Booij s procedure is usually used in the REF DIF 1 model More recently Kirby 1986b has developed an extension to the Booij approximation b
109. kb 1 sin thi i thetO0 1 y 3 dir ifreg 1i 1 2 6 continue 7 continue endif endif c If iinput 2 read a from data file fname4 if iinput eq 2 then open 4 file fname4 read 4 a 1 3 j71 n close 4 c Calculate wave angle on the first row CCC NO THIS DOESN T WORK IT HAS A BACKWARDS DERIVATIV E CCC NEW Tony add this part that calculates angle at the first row do 15 j 1 n if a m 3 EQ0 0 0 then akx2 0 else akx2 aimag clog a m J endif if a m 1 j EQ 0 0 then akx1 0 else akxl aimag clog a m 1 3 95 endif if abs akx2 akx1 GT pi then akx sign 2 pi abs akx1 abs akx2 dx akx1 else akx akx2 akx1 dx endif if j NE n then if a m j 1 EQ 0 0 then aky2 0 lse ky2 aimag clog a m 3 1 ndif f a m j EQ 0 0 then kyl 0 lse kyl aimag clog a m j ndif lse f a m j EQ 0 0 then ky2 0 lse ky2 aimag clog a m j ndif f a m j 1 EQ 0 0 then kyl 0 oo 0 MY FO o 0 w oO FO WM O WM F E n 0 kyl aimag clog a m j 1 endif endif if abs aky2 aky1 GT pi then else aky aky2 aky1 dy endif thet 3 atan2 aky akx kb m thet 3 180 thet 3 pi 15 continue write 7 202 thet 3 j 1 n nd C END NEW C pass angle on first row Fengyan 02 04 2 do 3j 1 n nd jj 3 1 nd 1 Pass Thet
110. l array variables in namelist makes it necessary that the dimension of the array in the program generating the data be the same as the dimension in the program reading the data For this reason we have isolated the parameter statement in a file param h which is then used to dimen sion all of the programs This file may be edited in isolation after which all programs which are to be used pre and post processing as well as REF DIF 1 should be recompiled For UNIX users this updating is automated by the included Makefile Stored output data files Prior to version 2 5 output was directed to two files outat dat was used primarily to store the complex 22 amplitude data which could later be used to construct either a wave height field or an image of the instanta neous surface Data was stored at the reference grid spacing In instances where a large amount of internal subdividing was being done this procedure was inadequate for the construction of a picture of the surface since the surface undulations are not resolved at the reference grid spacing The remainder of output data in older versions was directed either to the screen or to a file rundat dat This included header information and error log as well as x location reference phase height and wave angle at each reference grid point It has been difficult to use this information conveniently since the appearance of a warning or error message in the output could disturb the file format
111. lear what is believed to be a consequence of the rest of this License 8 If the distribution and or use of the Program is restricted in certain countries either by patents or by copyrighted interfaces the original copyright holder who places the Program under this License may add an explicit geographical distribution limitation excluding those countries so that distribution is permitted only in or among countries not thus excluded In such case this License incorporates the limitation as if written in the body of this License The Free Software Foundation may publish revised and or new versions of the General Public License from time to time Such new versions will be similar in spirit to the present version but may differ in detail to address new problems or concerns Each version is given a distinguishing version number If the Program specifies a version number of this License which applies to it and any later version you have the option of following the terms and conditions either of that version or of any later version published by the Free Software Foundation If the Program does not specify a version number of this License you may choose any version ever published by the Free Software Foundation If you wish to incorporate parts of the Program into other free programs whose distribution conditions are different write to the author to ask for permission For software which is copyrighted by the Free Software Found
112. lution appropriate to tidal models An illustration of such a feature is shown in Figure 6 where a poorly resolved feature occupies portions of four reference grid cells For cases such as this the program includes the option for the user to input user defined subdivided grid cells in order to specify these features at the level of the computational grid The use of user defined subgrids implies that the user will be choosing the grid spacing of the subgrids Since grid spacings must be uniform across a grid block the user should choose the option ispace 1 and specify the values of md in the input file indat dat Then the program will look for a data array of the correct pre specified dimension when it reads the subgrid data from the file subdat dat unit number iun 2 Several aspects of the subgrid data should be noted Data for depth d and current velocities u and v need to be specified at each of the subgrid points These data are input in the same units and with the same datum for depth as the data for the reference grid Also note that the subgrid includes the points on its boundaries for example the single subgrid shown in Figure 6 has a dimension of 6 by 6 The data values on the outer borders of the subgrids should be setup to match with the linearly inter polated values in the region external to the user defined subgrid If the data do not vary linearly along the subgrid boundaries there may be some mismatch between the subgrid bound
113. m iun 5 in both the inref and inwave subroutines The data is arranged in several lists and put in the indat dat file using the namelist convention The various namelist groups are defined according to the following prototype namelist statement namelist ingrid mr nr iu ntype icur ibc ismooth dxr dyr dt ispace nd iff isp iinput output inmd md fnames fnamel fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell fnamel2 fnamel3 fnamel4 fnamel5 fnamel6 fnamel7 fnamel8 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 wavesla iwave nfreqs waveslb update_interval num_data fregs tide nwavs amp dir waveslc thet0 freqs tide edens nwavs nseed waves2 fregin tidein The definition of each input variable follows ingrid namelist group e mr nr Reference grid dimensions Maximum values are ixr iyr respectively iu is switch for physical units u 1 MKS u 2 English Program defaults to u 1 if an error is made on input e ntype ntype is switch for nonlinearity ntype 0 linear model ntype 1 composite model ntype 2 Stokes wave model e icur icur is switch for input current data icur 0 no currents input icur 1 currents input icur defaults to a value of zero if an input error is detected e ibc ibc is the boundary condition switch ibc 0 closed boundaries ibcz1 open boundaries ibc defaults to a value of zero
114. m the average of its neighbors by more than the tolerance value dt specified on input This is basically a data checking feature Printed values are in meters Message Depth dr m at reference grid location ir jr differs from the average of its neighbors by more than dt m Execution continuing Action None by program Data in file refdat dat should be corrected if wrong Error occurs in inref An ambient current value occurs which implies that the flow would be supercritical at the given loca tion This serves as both a check for anomalously large current values and an indicator of possible subsequent computational problems Message ambient current at reference grid location irjr is supercritical with froude number froude number execution continuing Action None by program Data in file refdat dat should be corrected 1f wrong Error occurs in inref If the user specifies that predetermined subgrids are to be read in while at the same time telling the program to performs its own subdivisions the computed dimensions of the subgrid may be different than those of the subgrid included in the input Runs requiring user specified subgrids should choose the ispace 1 option If an incompatible set of dimensions occurs the program will either garble the input array or run out of data Message Warning input specifies that user will be supplying specified subgrids isp 1 while pro gram h
115. model initialy developed by Kirby and Dalrymple 1983a which incor porates all of the effects mentioned above The practical analysis of the refraction of water waves has generally been carried out in the past using ray tracing techniques This technique does not include wave diffraction and therefore it is inaccurate whenever diffraction effects are important Often due to complexities in the bottom topography wave tracing diagrams have many intersecting wave rays which leads to difficulties in interpretation as the theory predicts infinite wave heights at these locations Recently finite difference refraction models have been developed which have the advantage of providing wave heights and directions on a model grid rather than on irregularly spaced rays see for example Dalrymple 1988 1991 The diffraction of water waves around simple structures such as an offshore breakwater has been obtained analytically for a constant water depth Sommerfeld 1896 Diagrams of the wave heights in the vicinity of such a structure have been presented by the Corps of Engineers 1973 also Wiegel 1962 For a cylindrical structure MacCamy and Fuchs 1954 presented the constant depth solution These solutions give not only the wave height transmitted past the structure but also the scattered or reflected wave radiating away from the structure Generalized versions of these diffraction problems using numerical techniques and the Green s function me
116. n The problem of water waves propagating over irregular bathymetry in arbitrary directions is a three dimensional problem and involves complicated nonlinear boundary condition Very few solutions to the three dimensional problem exist and those that do are only for flat bottoms In one horizontal dimension sophisticated mod els by Chu and Mei 1970 and Djordjevic and Redekopp 1978 predict the behavior of Stokes waves over slowly varying bathymetry In order to simplify the problem in three dimensions Berkhoff 1972 noted that the important properties of linear progressive water waves could be predicted by a weighted vertically integrated model The vertical integration reduces the problem to only the two horizontal dimensions x and y Berkhoff s equation is known as the mild slope equation It is written in terms of the surface displacement n x y The equation in terms of horizontal gradient operator is C Vn CCgVan on 0 1 Here C vy g k tanh kh the wave celerity Cy C 1 2kh sinh2kh 2 the group velocity where the local water depth is h x y and g is the acceleration of gravity The local wave number k x y is related to the angular frequency of the waves c and the water depth h by the linear dispersion relationship o gktanh kh 2 The model equation 1 is an approximation however it is quite good even for moderately large local bottom slopes see Booij 1983 In both deep and shallow water it
117. n diss calculates frictional dissipation coefficients based on values of the switches read in by inref 10 wvnum Called by model grid and con wvnum performs a Newton Raphson solution of the linear wave current dispersion relation to obtain values of the wavenumber k 11 rand1 Called by model This function is a simple random number generator used to initialize the random wave phases if the directional spreading model is being used 12 acalc Called by model acalc normalizes the directional spectrum energy density over a 90 sector 13 bnum Called by acalc bnum computes the Bernoulli number n k n k 14 fact Called by bnum fact computes the factorial n of an integer n 3 4 Special Installation Instructions Several features of the program REF DIF 1 may require some modification or customizing during program installation on various systems REF DIF 1 is written using the features of FORTRAN 77 No use is made of vectorized solution techniques so the program should be useable on a wide range of systems with little or no modification If the program is to be used on machines with no upward limit placed on the size of the compiled executable code only one variable has to be checked during the initial installation of the program This is the logical device number for the controlling input date file This number may vary from system to system It appears in the program as iun 5 logical device number for controlling in
118. n 1 0 was found to not be sufficiently random and a new version is supplied with Version 2 0 ibc Open Boundary Condition Parameter Version 2 0 contains an option to use open lateral boundary conditions which are designed to be rea sonably transparent both to entering and exiting waves if the topography near the boundary is reasonably uniform The theory for these conditions are contained in the reference Kirby 1986c It is recommended that initial runs for a particular site be done with the default reflective conditions in order to see the magnitude of the boundary effects and then to use the open conditions for final runs The open boundary condition is invoked by choosing a value of ibc 1 in the input data icur No Currents Parameter A new parameter icur is included in the input data and determines whether or not the program is to read in current values from the input data files This change provides the option of not having to include zero current values on input whenno currents are being considered 20 3 1 2 Changes Appearing in Version 2 1 Version 2 1 represents a minor modification The change consists of a revision of the input file format struc tures for the file indat dat The formats for this file have been replaced by free format read statements so that the user can enter data separated by comments without regard to the column structure Note that data entered in the previously defined formats will still be read properly
119. n Eng 112 78 93 Kirby J T and Dalrymple R A 1986b An approximate model for nonlinear dispersion in monochro matic wave propagation models Coast Eng 9 545 561 Kirby J T Tjan K and Shi F 2002 Damping of high frequency noise in the large angle parabolic equation method manuscript 165 Krommes J A 1992 The web system of structured software design and documentation for C C Fortran Ratfor and TEX draft report Liu P L F and R A Dalrymple 1984 The damping of gravity water waves due to percolation Coastal Engineering Liu P L F and T K Tsay 1984 On weak reflection of water waves J Fluid Mech 131 59 71 MacCamy R D and R A Fuchs 1954 Wave forces on piles a diffraction theory Tech Memo 69 Beach Erosion Board Medina R 1991 personal communication Mei C C and E O Tuck 1980 Forward scattering by thin bodies SIAM J Appl Math 39 178 191 Mitsuyasu H F Tasai T Suhara S Mizuno M Ohkusu T Honda and K Rikishii 1975 Observations of the directional spectrum of ocean waves using a cloverleaf buoy J Physical Oceanography 5 750 760 Phillips O M 1966 The Dynamics of the Upper Ocean Cambridge University Press Radder A C 1979 On the parabolic equation method for water wave propagation J Fluid Mech 95 159 176 Sommerfeld A 1896 Mathematische theorie der diffraction Mat
120. n grid rows ir and ir 1 The reference points are separated by spacings dxr and dyr which are uniform in the x and y directions The spacings dxr and dyr may have arbitrary independent values Values of dr ur and vr constitute the reference grid data Section 2 8 describes the required input data file for these quantities The computational grid for any particular application should be chosen with care Since REF DIF 1 tries to use at least 5 points per wavelength the length of the computational domain in the propagation direction is restricted with the given parameter statements to just over 200 wavelengths This range can be extended by increasing ixr and or ix in the parameter statements The width of the model domain should be chosen 29 such that interference from the boundaries does not affect the study area If the reflective boundary conditions are used the extent of boundary influence is usually obvious particularly when viewing 2 dimensional plots of the amplitude envelope A z y 3 5 1 Grid Subdivision The only major feature of REF DIF 1 which is not described in Chapter 1 is the ability to subdivide the given reference grid into a more finely subdivided computational grid This would usually be done in cases where the reference grid spacing is too large to be used directly for calculations In this case the user may specify how the reference grid is to be subdivided or the user may tell the program to attempt its own s
121. n of this chapter describes in full the model s application to a specific problem Following a description of the problem and an indication of the type of results desired the input data files for the program are displayed and explained These data files are then used to run the program REF DIF 1 with no job specific modifications to the program involved Program output is then presented in such a way as to adequately indicate the results although in application individual users may wish to alter the nature of the presentation of output data The output for the various examples has been presented using some plotting programs which are external to the main body of the supplied program REF DIF 1 These specialized programs have been included in order to provide some guidance in reading the data files generated by the main program However plotting routines are likely to vary from one computer system to another The extra programs are therefore likely to be extensively machine specific to the systems on which the computations were performed Section 3 1 presents calculations of waves around an artificial surface piercing island This example makes particular use of the breaking wave thin film and shallow water dispersion relation capabilities of the model Section 3 2 provides calculations for waves propagating over a submerged elliptic shoal resting on a plane beach This example has been studied experimentally and provides a means for checking the acc
122. nce grid indicating a model of approximately 5r by 5r in x and y We sited the island center at x 460 ft and y 10ft where x and y are measured from the computational grid corner The island and measurement points are shown in relation to the grid in Figure 12 53 jr 100 90 80 70 60 50 40 30 20 10 L O O o L o o o L 1 1 L L 0 20 40 60 80 100 Figure 12 Measurement points in relation to reference grid 54 4 1 2 The Input Data Files One run of the model was performed using the specified input conditions The input data file indat dat for the run follows The reference grid data was stored in file refdat dat Sfnames fname2 refdat dat fname3 subdat dat fname4 wave dat fname5 owave dat A fname6 surface dat fname7 bottomu dat fname8 angle dat fname9 d fnamel0 refdifl log fnamell height dat fnamel2 sxx dat T fnamel3 sxy dat fname14 syy dat fnamel5 depth dat end Singrid mr 100 nr 100 iu 2 ntype icur ibe dxr 20 00000 dyr 20 00000 dt 10 00000 ispace nd iff isp iinput ioutput 1 Send Swavesla iwave 1 nfreqs 1 Send Swaveslb freqs 10 00000 tide 0 0000000E 00 nwavs 1 amp 14 00000 dir 0 00000001 Send ll oO ll HO Oq oO o x o Eri o o 55 4 1 3 Model Results The output for the
123. ntents are read in by subroutine inref If icur 0 only the depth data dr need to be specified Data for this file should be written in the following format do 1 ir 1 mr write unit dr ir jr jr 1 nr 1 continue then if icur 1 do 2 ir 1 mr write unit ur ir jr jr 1 nr 2 continue do 3 ir 1 mr write unit vr ir jr jr 1 nr 3 continue format 501f10 4 The data may be in either MKS or English units set the units switch iu in the iun 5 data file indat dat accordingly Subgrid Data File This file named subdat dat is accessed as logical device number iun 2 If no user defined subgrids are to be read in this file may be omitted The file consists of two parts an integer array of 1 s and 0 s indicating which reference grid cells are to be defined by the user and then a sequence of groups of arrays of d u and v one group for each subgrid The integer array isd is dimensioned mr 1 by nr 1 with one point for each spacing in the reference grid The array contains a O if that cell is not to be user defined and a one if it is For example the mr nr 7 6 reference grid shown in Figure 6 has four cells which are to be read in as user defined subgrids The array isd should be written in the data file first using the format do 1 ir 1 mr 1 write unit isd ir jr jr 1 nr 1 42 T continue format 1514 Now a group of arrays of d u and v must be entered into the data file
124. ntifies the program and then prints out messages identifying the parameters which set up the model run as read in by subroutine inref A sample title page is shown in Figure 7 Title page 2 gives the input wave conditions which were read in by inwave A sample title page 2 is given in Figure 8 43 Refraction Diffraction Model for Weakly Nonlinear Surface Water Waves REF DIF 1 Version 2 5 Center for Applied Coastal Research Department of Civil Engineering University of Delaware Newark Delaware 19716 James T Kirby and Robert A Dalrymple November 1994 input section reference grid values reference grid dimensions mr 100 nr 100 reference grid spacings dxr 5 0000 dyr 5 0000 physical unit switch iu 1 input in mks units icur 1 current values read from data files ibc 0 closed reflective lateral boundaries ispace 0 chosen program will attempt its own reference grid subdivisions y direction subdivision according to nd 1 ntype 1 stokes model matched to hedges model switches for dissipation terms 0 turbulent boundary layer 0 porous bottom 0 laminar boundary layer isp 0 no user defined subgrids iinput 1 program specifies initial row of a Figure 7 Sample title page 1 44 input section wave data values iwave 1 discrete wave amps and directions the model is to be run for 1 separate frequency components frequency component 1 wave period 8 0000sec tidal offset
125. o included and shows that the model is quite accurate in predicting the wave field 62 20 Figure 17 Results for waves propagating over a submerged shoal wave height contours 63 4 3 Waves Interacting with a Rip Current An example of waves which are normally incident on a planar beach and interact with a steady rip current flowing offshore from the beach has been included here in order to show the effects of wave current interac tion in the model This example was first used by Arthur 1950 to illustrate the effects of currents and depth changes acting together on results of ray tracing schemes However it also provides an important example of the usefulness of the combined refraction diffraction scheme since it represents a case where ray tracing breaks down due to the crossing of wave rays The velocity field studied by Arthur is shown together with computed wave orthogonals in Figure 18 Denoting a coordinate x pointed offshore the opposite of the x we will be using the velocity distribution 1s given by U 0 022952 e 76 2 2 4 7 62 2 56 V 0 2188 2 2 16 2 e 1762 P evt y 76 22 sign y 57 where the velocities are in m sec In terms of z the bottom topography is given by nia 02x 58 Arthur ran his calculations for a wave period of T 8 seconds A photograph of the wave field created in a laboratory study of waves interacting with an ebb ti
126. o you want to create indat dat yes 1 read ians if ians eq 1 then open iun 4 file indat dat write fnames portion of namelist file 146 write iun 4 nml fnames write ingrid portion of namelist file 102 103 104 105 106 107 108 if iu eq 0 then write iun 3 input iu l mks 2 english read endif iu write iun 3 input dispersion relationship ntype O linear write iun 3 l composite 2 stokes write iun 3 input lateral boundary condition ibc O closed write iun 3 l open read write iun 3 102 format read write read input ispace 0 program picks x spacing l user choses ispace input nd y divisions 1 is minimum nd if ispace eq 0 go to 105 write iun 3 constant or variable x spacing 0 for constant read if iansl eq 0 then write read do 103 iko 1 mr 1 md iko mdc continue else write iun 3 104 format read endif iansl input constant md mdc input md i for i 1 to mr 1 md i i 1 mr 1 write iun 3 106 format input if 1 turbulent if 2 porous if 3 laminar write iun 3 standard choice 1 0 O read iff 1 iff 2 iff 3 write iun 3 107 format read input isp subgrid features standard 0 isp write iun 3 108 format
127. on output write 5 120 write 5 106 if Master_Start ge 0 then Read control data from unit 1 read 1 nml ingrid if ispace eq 1 read 1 nml inmd write 5 107 mr nr dxr dyr if iu eq 1 write 5 114 iu if iu eq 2 write 5 115 iu if icur eq 0 write 5 200 if icur eq 1 write 5 201 if ibc eq 0 write 5 202 if ibc eq 1 write 5 203 if ispace eq 0 write 5 108 if ispace eq 1 write 5 109 write 5 119 nd 82 Cc c if ntype eq 0 write 5 110 if ntype eq 1 write 5 111 if ntype eq 2 write 5 112 Check input from unit 1 if mr gt ixr or nr gt iyr then write 5 dimensions for reference grid too large stopping call exit 1 not good for intel compiler stop end if if iu ne 1 and iu ne 2 iu 1 dt dt dconv iu dxr dxr dconv iu dyr dyr dconv iu if dt eq 0 dt 2 if nd gt ifix float iy 1 float nr 1 then write 5 102 call exit 1 not good for intel compiler stop endif if ispace eq 1 then test 0 do 1 i 1 mr 1 if md i gt ix 1 then write 5 103 1 test 1 endif continue if test eq 1 then call exit 1 not good for intel compiler stop endif endif endif skip above after the first call wave module do 7 i 1 mr do 7 j 1 nr dr i j Depth Wave i j ur i j Intp_U_Wave i j vr i j Intp_V_Wave i j continue Check for large depth changes and large currents in reference grid data 83 do 8 i 2 mr 1 do 8 3 2 nr 1 dcheck dr i41
128. on output are meaningless if multiple direction compo nents are being used and for single component runs become meaningless if the waves become short crested or the crests become significantly curved Finally a file surface dat is generated if isurface is set to one This file provides the same type of infor mation that was put in outdat dat in older versions of the program with the exception that the present version stores data for every computational point in the domain 2 y j 1 n Then for each x position in the computational grid the program stores the following information 3 x 4 ae ij j l n 48 Since it is not usually known a priori how many steps will be taken in the x direction this file is ended by writing in a negative value of x This file is processed by the program surface fto give data on a regularly spaced grid An anotated listing of the REF DIF 1 program code is given in Appendix A Various preprocessing and postprocessing programs are also listed in Appendix B for normal usage and Appendix C for LRSS usage 49 4 EXAMPLE CALCULATIONS This chapter presents calculations performed using the combined refraction diffraction model REF DIF 1 The problems studied here were chosen as representative tests of various features of the model Further examples illustrating the use of the computational schemes upon which the program is based may be found in the technical report by Kirby 1983 Each sectio
129. on relationship developed by Hedges 1976 is provided as an option in the model This relationship between the frequency and the water depth is c gk tanh kh 1 A h 17 In shallow water this equation matches that of a solitary wave while in deep water it asymptotically approaches the linear wave result neglecting real amplitude dispersive effects For this reason a model with a dispersion relationship which is a smooth patch between the Hedges form valid in shallow water and the Stokes relationship valid in deep water is used This hybrid model is described in Kirby and Dalrymple 1986b There are as a result of the different dispersion relationships possible three options in REF DIF 1 1 a linear model 2 a Stokes to Hedges nonlinear model and 3 a Stokes model Of these options the second will cover a broader range of water depths and wave heights than the others 3 The wave direction is confined to a sector 4 70 to the principal assumed wave direction due to the use of the minimax wide angle parabolic approximation of Kirby 1986b 2 3 Energy Dissipation 2 3 1 General form Energy dissipation in the model occurs in a number of ways depending on the situation being modelled An energy loss term due to Booij 1981 and expanded by Dalrymple et al 1984a permits the model to treat bottom frictional losses due to rough porous or viscous bottoms surface films and wave breaking The linear form of
130. ools beyond the usual Fortran compiler We caution that namelist is not a standard Fortran 77 feature however we have not found any compilers yet which do not provide this feature Revision to indat dat file structure The most readily apparent change to the long term user of REF DIF 1 is the change to the use of namelist to structure the indat dat data file The structure of the file and the meaning of each input variable are described in section 2 7 and the file for each example is given in Chapter 3 The program datgen f supplied with older versions of the program has been updated and renamed to dat genv25 f The new version produces the namelist formatted indat dat file Since long term users are likely to have their own versions of old indat dat files we have also provided a new program indat convertv25 f which converts old indat dat files to new indat new files which should then be renamed A quick test of the integrity of this new data file convention could be made by using the old datgen to generate an indat dat converting it to indat new and then comparing it to the file indat dat generated by the new datgenv25 The two resulting files should be identical Use of param h file to dimension arrays Changing dimensions in REF DIF 1 in the past has involved a careful search through a number of sub routines to get all parameter statements revised properly This has not been a pleasant process In addition the specification of severa
131. ountered during the installation and use of the program on different computer systems section 3 4 Section 3 5 presents the two levels of grid information used by the program Section 3 6 describes the option of reading in the first row of complex amplitude values A from an additional external data file wave dat The input data file structure is then discussed The program reads data in two essentially separate groups The first group of data establishes the size of the model grid and gives the wave conditions to be studied this group is discussed in section 3 7 The second group of data gives the reference grid data values and defines any user specified subgrids this is discussed in section 3 8 The structure of the program output is discussed in section 3 9 A listing of the program is included in sections 3 and 6 Section 7 provides listings of several pre and postprocessing programs provided with the distribution of REF DIF 1 3 3 Program Outline and Flow Chart The model REF DIF 1 is organized in one main program WaveModule and fourteen subroutines The pro gram does not contain calls to any external package and should be free standing on any system The program should be a legal code for any compiler which accepts the full FORTRAN 77 standard language A possible exception would be the appearance of namelist statements which are not part of the FORTRAN 77 standard but which are provided by all compilers we have checked REF DIF 1 is struct
132. ow is applied to output results depths at reference grid points gt 0 submerged areas lt 0 elevation above surface datum x velocities at reference grid points only entered if cur 1 y velocities at reference grid points only entered if icur 1 Data is entered in namelist format from the data file indat dat Subroutine is called from WaveModule and returns control to calling location unless a fatal error is encountered during input data checking Center for Applied Coastal Research Department of Civil and Environmental Engineering University of Delaware Newark DE 19716 Coded by James T Kirby October 1984 Version 2 6 Revised February August 2002 refdif subroutine inref incl incl ude ude param h common h 79 include pass h integer i j namelist ingrid mr nr iu ntype icur ibc ismooth dxr dyr dt ispace nd iff isp iinput ioutput inmd md fnames fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell fnamel2 fnamel13 fnamel4 fnamel5 fnamel6 fnamel7 fnamel8 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 if Master_Start ge 0 then c Constants g 9 80621 c Specify name of namelist data file fnamel indat dat open unit 1 file fnamel 80 5 2 1 Read file names from namelist refdif read 1 nml fnames c Open standard input and output files E open unit 2
133. plitude data for constructing an image of the instantaneous water surface at the computational resolution If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel0 sxx dat Sz components at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel1 sxy dat Sz components at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel12 syy dat Syy components at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel3 fx dat depth integrated wave forcing in x direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel4 fy dat depth integrated wave forcing in y direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel15 qx dat short wave flux in x direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated 38 e fnamel6 qy dat short wave flux in y direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel7 bottomu dat magnitude of bottom velocity at reference grid po
134. program The second revision is the replacement of the add on noise filtering algorithm with an algorithm that functions as an imbedded part of the model equation and finite difference approximation as described by Kirby 1993 Additionally several inconsistencies appearing in input error checking have been corrected A more robust algorithm for computing wave angles due to Medina 1991 has been added 3 1 6 Changes Appearing in Version 2 5 A substantial number of changes appear in Version 2 5 These changes represent a combination of data format changes and enhanced post processing with the addition of several new computed output variables Since the use of Matlab is becoming more widespread we have included a m script which we are presently using to present computed results 21 Many of the changes in version 2 5 were made in order to include the model in the Littoral Remote Sens ing System LRSS developed for the U S Navy In cases where the changes in data formats were consistent with the vast majority of available Fortran 77 compilers as in the use of namelist formats the changes were made directly in REF DIF 1 In cases where the standard required the use of a non typical format as in the use of the machine independent binary HDF formats for large arrays then the data transfer to and from LRSS is handled using a pre and post processing layer The intent is that the normal user of REF DIF 1 should not need access to any t
135. pted to control the high frequency noise generated by the breaking process using smoothing filters which are applied to the computed complex amplitude in between each computational row In version 2 5 of REF DIF 1 this procedure was replaced by a new algorithm described by Kirby et al 2002 The damping is built into the computational algorithm and is turned on automatically if breaking has started anywhere in the computational domain 3 USER S MANUAL This section provides the operating manual for the program REF DIF 1 Section 4 provides sample docu mentation and calculations for example problems A separate Fortran program datgenv26 f is provided which generates the input data files for these as well as a number of additional examples In addition the programs indat createv26 f which assists the user in constructing the indat dat file is provided 3 1 REF DIF 1 Revision History 3 1 1 Changes Appearing in Version 2 0 Several changes have been made to the program released as Version 1 0 of REF DIF 1 These changes are outlined here for convenience Directional spectra application The routine used to specify a directional spectrum as the initial condition for the wave calculations has been extensively revised and tested The present model is based on a Mitsuyasu type spreading factor and apportions wave directions and energy density in randomly sized directional bins In addition the original random number generator supplied with Versio
136. put 1 113 write iun 3 112 format input en density and on next line directional 1 spreading factor read edens ifreq read nwavs ifreq nseed 500 endif continue if iwave eq 1 write iun 4 nml waveslb if iwave eq 2 write iun 4 nml waveslc endif if iinput eq 2 then line 9 iinput 2 write iun 3 input wave period and tide stage read freqin tidein write iun 4 nml waves2 endif close iun 4 endif stop end c subroutine con include param h common ref d ixr iyr u ixr iyr v ixr iyr m n dx dy itype common ind iu ntype icur ibc ispace nd iff isp iinput liwave nfreqgs freqs tide nwavers amp dir edens common dims x ixr y iyr dimension iff 3 dimension amp ncomp ncomp dir ncomp ncomp tide ncomp 1 freqs ncomp edens ncomp nwavs ncomp 149 do 1 i 1 m do 1 3 1 n if icur eq 1 and itype eq 3 then xp x m x i u i j 2 0 02295 exp xp 76 2 2 2 exp y 3 7 62 2 1 2 xp v 1 3 0 2188 2 xp 76 2 2 exp xp 76 2 2 2 lerf3k abs y 3 107 76 y 3 abs y 3 else u i j v i j endif continue 0 0 return end subroutine depth iun c include param h common ref d ixr iyr u ixr iyr v ixr iyr m n dx dy itype common ind iu ntype icur ibc ispace nd iff 3 isp iinput liwave nfreqs freqs tide nwavers amp dir edens
137. put data file initialized near the top of subroutine inref The supplied version of REF DIF 1 has this value set to iun 5 5 All output is directed to data files with pre chosen unit numbers Many systems require that access to disk files be initialized and terminated by the use of open and close statements in the program code Since the parameter list of the open statement varies from system to system the user should take care that the open statements are compatible with the system software being used The open statements appear near the top of refdif1 and inref The corresponding close statements appear near the end of refdif1 3 5 Computational Grids and Grid Interpolation The reference grid terminology is defined in Figure 3 The grid consists of a mesh of points with dimensions mr x nr in x and y The values of mr and nr must be less than or equal to the parameters ixr and yr whose 28 i NN JR I IR I e 3 4 5 Grid Block 3 Grid Row 3 Figure 3 Reference grid notation values are set in the parameter statement in param h Reference grid data of depth dr and ambient current components ur and vr are defined at the grid points The program assumes that the x y coordinate system is established with the origin at grid point ir jr 1 1 In this manual we make the distinction between the terms grid row ir which is the row of points jr 1 nr at location ir and grid block ir which is the physical space betwee
138. read 3 100 isd ir jr jr 1 nr 1 15 continue endif c Input done return to main program return 100 format 1514 101 format 501 f10 4 102 format y direction subdivision too fine 1 maximum number of y grid points will be exceeded execution terminating 103 format x direction subdivision too fine on grid block 1 2x i3 execution terminating 104 format depth 2x f7 2 m at reference grid location 12 2x i3 differs from the average of its neighbors by 1 more than 2x f7 2 m execution continuing 105 format ambient current at reference grid location 1 2 2x i3 is supercritical with froude number f7 4 1 execution continuing 106 format 0 20x input section reference grid values 107 format reference grid dimensions mr 13 nr 13 reference grid spacings dxr f8 4 r dyr f8 4 108 format ispace 0 chosen program will attempt its own l reference grid subdivisions 109 format ispace 1 chosen subdivision spacings will be 1 input as data 85 110 111 112 113 114 115 116 120 aL Le 118 119 200 201 202 203 a p format format format format warn 1 specified s 1 while progr ntype 0 linear model ntype 1 stokes model matched to hedges model ntype 2 stokes model ing input specifies th
139. rent this model can handle The program can handle currents that are of the same order as the wave speed However opposing currents which are strong enough to stop the waves cause a singularity in the wave height due to the absence of reflection effects and no provision for steepness limited breaking These extensions will be added in the future 2 In which situations should the user specify totally reflecting lateral boundaries bc 0 and when should partially transmitting lateral boundaries bc 1 be used It is common practice to use the closed boundary condition as the amount of reflection can be identified and the domain width can be chosen so that there are no reflections in the region of interest The final run may be made using partially transmitting boundaries but it should be noted that there will still be a certain amount of reflected waves present in the domain 3 What is the difference between user specified subdivisions and user specified subgrid User specified subdivisions and user specified subgrids are two completely independent options The user has to specify the desired number of subdivisions in the y direction nd However the user has an option between specifying the number of subdivisions in the x direction md ir or having the program specify those If the program specifies md ir the switch space has to be set to zero and the program performs interpolations of the depth and velocity grids If the user c
140. rids will be read e iinput iinput specifies whether the program or the user will generate the first row of complex amplitude A values If iinput 1 the program constructs A based on the input conditions specified as follows If iinput 2 the user must specify A in an external data file wave dat e ioutput ioutput specifies whether the last computed row of complex amplitudes A are to be stored in file owave dat A value of ioutput 1 skips this option ioutput 2 turns the option on inmd namelist group 37 e md ixr If ispace 1 the x direction subdivisions are inserted here one for each reference grid block The array md must be dimensioned according to the value of ixr in the main program regardless of how many values actually are being used fnames namelist group e fnamel indat dat automatically assigned to the namelist input data file indat dat e fname2 refdat dat depth and u v on the reference grid e fname3 subdat dat user specified subgrids e fname4 wave dat user specified complex amplitude on row 1 for iinput 2 e fname 5 refdifl log run log for refdifl program e fname 6 height dat wave heights at reference grid locations This file is always generated e fname7 angle dat wave directions 0 in degrees at reference grid points This file is always generated e fname depth dat tide corrected depths at reference grid locations This file is always generated e fname9 surface dat complex am
141. s Establish initial wave conditions for the ifreq frequency if iinput eq 1 then Compute wave from data given in indat dat if iwave eq 1 then iwave eq 1 discrete components specified do 3 j 1 n a 1 3 cmp1x 0 0 do 2 iwavs 1 nwavs ifreq thet j dir ifreq itime 180 pi a 1 j3 a 1 3 amp ifreg itime cexp cmp1x 0 kb 1 sin ldir ifreq itime y 3 continue 3 continue if fname7 ne write 7 202 thet j j 1 n nd pass angle Fengyan 02 04 2002 do j 1 n nd jj j 1 nd 1 Pass_ enddo else Theta 1 33 thet j iwave eq 2 directional spreading model not recommended sp float nwavs ifreq nsp nwavs ifreq thmax pi 4 call acalc thmax nsp al edens ifreq sqrt edens ifreq al nn 31 li nn 1 2 1 seed randl seed Compute randomly distributed Delta theta s sum0 0 do 12 i 1 nn seed randl seed 94 dthi 1 seed sum0 sum0 seed 12 continue xnorm 2 thmax sum0 do 101 i 1 nn dthi i dthi i xnorm 101 continue thi0 thmax do 4 i 1 nn thi0 thi0 dthi i thi 1 thi0 dthi i 2 dth dthi i amp ifreq i edens ifreq sqrt dth sqrt cos thi 1 dth 2 1 2 nsp cos thi i dth 2 2 nsp 4 continue do 5 i 1 nn ipl i 1 seed randl seed dir ifreq ip1 2 pi seed 100 5 continue do 7 3 1 n a 1 3 cmp1x 0 0 do 6 i 1 nn a 1 j a 1 j amp ifreqg i cexp cmplx 0
142. s 000000 tide 0 0000000E 00 nwavs 1 amp 3200000E dir 0000000E 00 Send HHO O H Oo oobr o x o ll pi pi N ll o The supplied program datgenv25 f may be used to generate the depth grid refdat dat 61 4 2 2 Model Results The output files depth dat height dat surface dat for this example have been used to construct plots of instantaneous surface elevation and wave height In this case wave heights have been non dimensionalized using the incident wave height The resulting plots are shown in Figures 16 and 17 The plots show the effect of wave focussing over the shoal area and show that the prediction of the wave field beyond the shoal does not involve the problem of caustic or singularity formation common to ray tracing algorithms used to model similar situations 20r 15r P4 E gt 10 f J 5r j f 0 ly 0 5 Figure 16 Results for waves propagating over a submerged shoal surface elevation contours Comparisons of the predictions of this model using both the Stokes wave nonlinearity and the composite model of Kirby and Dalrymple 1986b are given in Kirby and Dalrymple 1986b a comparison to laboratory data is als
143. the mild slope equation with dissipation is A idPA T gt A 18 2 809 79 18 where the dissipation factor w is given by a number of different forms depending on the nature of the energy dissipation The factor w is the energy dissipation divided by the energy and has the units of time 5 2 3 2 Laminar surface and bottom boundary layers At the water surface and at the bottom boundary layers occur due to the action of viscosity For a contam inated surface resulting from surface films natural or otherwise a significant amount of damping occurs which is dependent on the value of the fluid viscosity From Phillips 1966 the surface film damping is ek Gas i i tanh kh where v is the kinematic viscosity The term under the square root sign is related to the thickness of the boundary layer which is generally small At the bottom the boundary layer damping is ote N i 20 By setting the input switch iff 3 in the REF DIF 1 model surface and bottom damping are computed at all locations in the model 2 3 3 Turbulent bottom boundary layer In the field the likely wave conditions are such that the bottom boundary layer is turbulent In this case an alternative specification of the energy dissipation must be provided Utilizing a Darcy Weisbach friction factor f the dissipation term can be shown to be 20kf A 1 i 3rsinh 2kh sinh kh en See Dean and Dalrymple 1984 In order to implement this damping term the val
144. this License Any attempt otherwise to copy modify sublicense or distribute the Program is void and will automatically terminate your rights under this License However parties who have received copies or rights from you under this License will not have their licenses terminated so long as such parties remain in full compliance 5 You are not required to accept this License since you have not signed it However nothing else grants you permission to modify or distribute the Program or its derivative works These actions are prohibited by law if you do not accept this License Therefore by modifying or distributing the Program or any work based on the Program you indicate your acceptance of this License to do so and all its terms and conditions for copying distributing or modifying the Program or works based on it 6 Each time you redistribute the Program or any work based on the Program the recipient automatically receives a license from the original licensor to copy distribute or modify the Program subject to these terms and conditions You may not impose any further restrictions on the recipients exercise of the rights granted herein You are not responsible for enforcing compliance by third parties to this License 7 Tf as a consequence of a court judgment or allegation of patent infringement or for any other reason not limited to patent issues conditions are imposed on you whether by court order agreement or otherwis
145. thod have yielded very powerful procedures for wave force calculation for cases where the drag force is much smaller than the inertia force In order to incorporate diffractive effects the general practice has been to suspend refraction in areas where diffraction is dominant and only permit diffraction there using Sommerfeld s analytic solution for a flat bottom Far from the diffraction area refraction is resumed This ad hoc technique clearly is inaccurate but does permit the inclusion of diffraction in an approximate way Combined refraction diffraction models include both effects explicitly thus permitting the modelling of waves in regions where the bathymetry is irregular and where diffraction is important Regions where wave rays cross due to local focussing or where caustics are caused by other means are treated correctly by the models and no infinite wave heights are predicted The models developed in parabolic form do not predict the waves which are scattered upwave that is waves which are reflected directly back the way they came are not modelled and are neglected This means in general wave reflection phenomena are not reproduced correctly Combined refraction diffraction models are uniquely suited for the calculation of wave heights and wave direction in areas where one or both of these effects are present Examples include the determination of wave heights in a bay given the offshore wave heights periods and directions and determ
146. ton Raphson iteration for wavenumber k may not converge in the specified number of steps This may occur for waves on strong opposing currents Message WAVENUMBER FAILED TO CONVERGE ON ROW T COLUMN J K last iterated value of wavenumber D depth T period calculated from last iterated value of k U x direction velocity F value of objective function should be 0 for convergence Action Program continues with last iterated value of k Computed results are of questionable accuracy Error occurs in wvnum Log of Calculations For each frequency component the program starts a new page of output and indicates the number of the component in the input stack The model then prints the x position the value of the reference phase function and the number of x direction subdivisions used for each reference grid row 47 3 9 2 Stored Output The program stores a set of data files with a resolution of the reference grid points Each of these files is written using a common format statement and may be read using the format do i 1 mr read unit 100 variable i j j 1 nr end do 100 format 200 f10 4 The available files are height dat wave height depth dat water depth with tide correction included angle dat wave angle in degrees sxx dat sxy dat syy dat radiation stresses not available yet bottomu dat magnitude of bottom velocity if requested It should be noted that the wave directions given
147. ts application to tsunami calcula tions Mar Geodesy 2 41 58 Kirby J T 1983 Propagation of weakly nonlinear surface water waves in regions with varying depth and current ONR Tech Rept 14 Res Rept CE 83 37 Department of Civil Engineering University of Delaware Newark Kirby J T and R A Dalrymple 1983a A parabolic equation for the combined refraction diffraction of Stokes waves by mildly varying topography J Fluid Mech 136 543 566 Kirby J T and R A Dalrymple 1983b The propagation of weakly nonlinear waves in the presence of varying depth and currents Proc XX Congress 1 A H R Moscow Kirby J T 1984 A note on linear surface wave current interaction J Geophys Res 89 745 747 Kirby J T and R A Dalrymple 1984 Verification of a parabolic equation for propagation of weakly non linear waves Coastal Engineering 8 219 232 Kirby J T 1986a Higher order approximations in the parabolic equation method for water waves J Geophys Res 91 933 952 Kirby J T 1986b Rational approximations in the parabolic equation method for water waves Coastal Engineering 10 355 378 Kirby J T 1986c Open boundary condition in parabolic equation method J Waterway Port Coast and Ocean Eng 112 460 465 Kirby J T and Dalrymple R A 1986a Modeling waves in surfzones and around islands J Waterway Port Coast and Ocea
148. ubdivisions The maximum reference grid dimensions have been set fairly arbitrarily and may be changed by modi fying the parameter statements at the beginning of each routine in the program The grid is large enough to comply with typical grids for various tidal computational models which may be used to specify the ambient currents and get small enough so that data values do not occupy too much internal storage It is anticipated that the spacings dxr and dyr if based on such a model may be too large for an accurate integration of the parabolic model to be performed Consequently the model has been arranged so that the reference grid block between any two adjacent reference grid rows may be subdivided down to a finer mesh in order to provide sufficient resolution in the computational scheme This subdivision process is performed internally in the program with an exceptional feature to be described below and may be controlled by the user or the program Remember that the computational scheme proceeds by marching in the x direction and therefore the only reference grid information required at any particular step is the data on row ir where computation starts and on row r 1 where computation ends The fine subdivided mesh is set up on the intervening grid block An example of a particular subdivision of a grid block is shown in Figure 4 Here the choices of x direction subdivision md ir 5 and nd 2 are illustrated with md and nd being the num
149. ue of f 0 01 was assumed In REF DIF 1 if the input data switch iff 1 1 then turbulent damping is computed at all locations 2 3 4 Porous sand Most sea bottoms are porous and the waves induce a flow into the bed This results in wave damping due to the Darcy flow in the sand For beds characterized by a given coefficient of permeability Cp the damping can be shown to be gkCp 1 i w cosh kh The coefficient of permeability Cp has the units of m and is of order 4 5 x 107 m Liu and Dalrymple 22 1984 show for very permeable sands that the damping is inversely related to Cp and a different w term must be used However this case is not likely to occur in nature Porous bottom damping is computed in REF DIF 1 when iff 2 1 2 3 5 Wave breaking For wave breaking the model is more complicated Dally et al 1985 showed that the rate of loss of wave energy flux is dependent on the excess of energy flux over a stable value This model has been tested for laboratory data for a number of different bottom slopes and predicts the wave height in the surf zone extremely well Kirby and Dalrymple 19862 show that the dissipation due to wave breaking is given as ip KC 1 Qm 7 23 where K and y are empirical constants determined by Dally et al to be equal to 0 017 and 0 4 respectively Here the wave height H 2 A By using this dissipation model and a breaking index relation H gt 0 78h to determine
150. uracy of the model calculations It also provides an example of the type of results provided by a combined refraction diffraction model in a situation where ray tracing predicts a strong convergence of wave rays with resulting singularities in the prediction of wave height Section 3 3 provides example calculations for the case of waves shoaling on a plane beach and interacting with a rip current This example illustrates the wave current interaction feature of the model 50 4 1 Waves Around an Artifi cial Island The first example involves the calculation of the wave field around an artificial surface piercing island of the type used in offshore operations The island is circular with a base radius of 400 ft and a crest elevation of 80 ft above the flat seabed The island radius at the crest is 160 ft leading to a side slope of 1 3 The water depth around the island is taken to be 60 ft The island geometry is shown in Figure 9 Incident Wave Figure 9 Artificial island geometry 51 Figure 10 Locations for wave height measurements The wave conditions to be studied are given by e Wave height H 28 ft e Wave period T 10sec e No currents The model was run with a depth of 60 ft away from the island and a tidal offset of 0 0 ft The required set of wave predictions consisted of wave height at 12 locations as indicated in Figure 10 The spacing between the points are in units of the base radius r 400 ft Note th
151. ured in two levels a main level which reads in and checks input data and then starts the operation of the wave model level and the model level itself which performs the actual finite difference calculations The flow charts for the two levels are given in Figures 1 and 2 A short description of each routine in the model follows 1 WaveModule Main program now defined to be a subroutine controls the calls to inref and inwave to read in data and to model which performs the actual calculations No calculations are performed by this routine 2 inref Called by WaveModule inref reads in data controlling reference grid dimensions and the grid interpolation scheme from logical device number iun 5 and reads in the reference grid values of depth dr x direction velocity ur and y direction velocity vr from logical device number iun 1 Some data checking is performed If data is read in in English units inref converts it to MKS units using the dconv factor Output files are initialized A description of the data files may be found in sections 2 6 2 8 of this manual At the end of the subroutine control is returned to refdifl1 3 inwave Called by WaveModule inwave reads in data specifying the initial wave field along the first row of reference grid points Data is read from logical device number iun 5 Conversion to MKS units is performed for data read in in English units Control is returned to refdif1 24 iNREF Read and check
152. urface f c c This program converts the file usually surface dat containing e an instantaneous snapshot of the water surface at the computational er grid spacing to a regularly spaced ascii file suitable for directly GR reading into Matlab format c Cx James T Kirby Er Center for Applied Coastal Research Ex University of Delaware Er Newark DE 19716 e kirby udel edu 302 831 2438 FAX 302 831 1228 c EX Last revision 12 22 94 c c program surfacev26 include param h integer i j k m n nx ny iswap integer iret iout real surface iy iy real x iy y iy dx dy xold iy surfold iy iy character 255 fileout integer idimsizes 2 character 255 fnamel fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell fnamel2 fname13 fnamel4 fname15 fnamel16 fnamel7 fnamel18 fnamel19 fname20 fname21 fname22 fname23 fname24 fname25 fname26 namelist ingrid mr nr iu ntype icur ibc ismooth dxr db ispace nd iff isp iinput ioutput inmd md fnames fname2 fname3 fname4 fname5 fname6 fname7 fname8 fname9 fnamel0 fnamell1 fname12 fnamel13 fnamel4 fnamel5 fnamel6 fnamel 7 fnamel8 fnamel9 fname20 fname21 fname22 fname23 fname24 fname25 fname26 wavesla iwave nfreqs waveslb freqs tide nwavs amp dir waveslc thet0 freqs tide edens nwavs nseed waves2 fregin tidein open 8 file indat dat 157 c 10 20 25 read 8 nml fnames open 10 file fn
153. verse y direction kz ky represent wave 1 Wave angles are calculated by number components in Cartesian coordinates x y 2 5 3 Radiation Stresses and Forcing Terms For a depth integrated time averaged circulation model such as SHORECIRC the radiation stress Sag is defined as 1 Sag l pag puuouwg dz 9 pgh 27 ho where Uwa is the horizontal shortwave induced velocity h represents mean water surface elevation The generalized radiation stress tensor is given by Sag apSm ap Sp 28 where cos 9 sin cos Pant sin cos 0 sin 0 29 The scalars Sm and 5 are defined as Ou pu2 dz 30 ho 1 E Sp f pida pg 31 ho where n C h Outside the surfzone the results for sinusoidal waves are used and thus 30 and 31 may be written as 1 Sm AU PA G 32 Sy l HG 33 where 2kh __ 4 sinh 2kh Gm Inside the surfzone those results are augmented with the contribution from the roller using the results from Svendsen 1984 Thus Sm and Sp are determined as Sm p9H 35 pet gh H L and Sp 59H Bo 36 where A is the roller area and Bo is the wave shape parameter defined by T Bo 7 Da 37 For a time averaged 3D circulation model radiation stresses are defined by the surface stress 55 and the body stress 9 a3 88 following S25 agp agpw2 38 S o PU 39 The body stress and surface stress may be evaluated by 1 r S
154. w 1 for iinput 2 refdif1 log run log for refdifl program height dat wave heights at reference grid locations This file is always generated angle dat wave directions 0 in degrees at reference grid points This file is always generated depth dat tide corrected depths at reference grid locations This file is always generated surface dat complex amplitude data for constructing an image of the instantaneous water surface at the computational resolution If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel0 sxx dat Sz components at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel1 sxy dat Sz components at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fname12 syy dat Syy components at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel3 fx dat depth integrated wave forcing in x direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel4 fy dat depth integrated wave forcing in y direction at reference grid locations If REF DIF 1 is given a null string as the input for this file name no file is generated e fnamel15 qx dat short wave flux in x direction at refer
155. x i j cos Thetar i j Pass MassFluxV i j Pass MassFlux i j sin Thetar i j enddo enddo 911 continue cdiss cubott gam 0 3 do j 1 ny wave do i 1 nx wave Pass Diss i j 20 125 grav Pass Height i Jj x Pass Height i j 0 15 Pass_Cg 1 3 Depth_wave i j ie 1 gam Depth_wave i J Pass_Height i j 2 Pass_ibrk i j Pass ubott i J Pass_Height i j pi period sinh Depth_wave i j 2 pi period Pass C i j enddo enddo C depth integerated short wave forcing do j 2 ny wave 1 do i 2 nx wave 1 Pass Wave Fx i j Pass Sxx i 1 j Pass Sxx i 1 j 2 Pass Sxy i j 1 Pass Sxy i j 1 2 dyr Pass Wave Fy i j Pass Sxy i 1 j Pass Sxy i 1 j 2 Pass Syy i j 1 Pass_Syy i j 1 2 dyr enddo enddo do i 2 nx wave 1 j 1 Pass Wave Fx i j Pass Sxx i 1 j Pass Sxx i 1 j 2 Pass Sxy i j 1 Pass Sxy i j 1 dyr Pass Wave Fy i j Pass Sxy i 1 j Pass Sxy i 1 j 2 Pass Syy i j 1 Pass_Syy i j 1 dyr enddo do i 2 nx wave 1 j ny wave Pass Wave Fx i j Pass Sxx i 1 j Pass Sxx i 1 j 2 Pass Sxy i j Pass Sxy i j 1 1 dyr Pass Wave Fy i j Pass Sxy i 1 j Pass Sxy i 1 j 2 Pass Syy i j Pass Syy i j 1 1 dyr enddo do j 22 ny wave 1 129 dxr dxr dxr dxr dxr dxr i 1 Pass Wave Fx i Pass Wave Fy i enddo iJi j do j 2 ny_wave 1
156. y aky2 aky1 2 dy endif endif thet 3 atan2 aky akx kb m 15 continue do Estimation of Wave Height 17 j 1 n height 3 2 cabs a m 3 17 continue Spatial smoothing of wave height if nd gt 1 if jh do nd NE 1 then int nd 2 j 1 n nd 114 heightb 3 0 if j eq 1 then do jj 1 1 3h heightb 1 heightb 1 height jj 2 end do heightb 1 sqrt heightb 1 3h 1 endif if j eq n then do jj n jh n heightb n heightb n height 33 2 end do heightb n sqrt heightb n jht1 endif if j gt 1 and j 1t n then do jj j jh jtjh heightb j heightb j height 33 2 end do heightb 3 sqrt heightb 3 2 3h 1 endif end do endif if nd eq 1 then do j 1 n heightb 3 height 3 end do endif c If ismooth 1 perform a much more drastic smoothing over a wider c footprint if ismooth eq 1 then do j 1 n nd height 3 heightb 3 end do 115 Cc First employ a simple average over a wide base to get rid of most noise and most medium scale variability do j 1 nrs nd n nrs nd nd heightb 3 0 do is nrs nrs 1 heightb 3 heightb 3 height j is nd 2 end do heightb 3 sqrt l heightb 3 float 2 nrs 1 end do Then emply a weighted thr point average to get rid of small scale noise do 3j 1 n nd height 3 heightb j end do do 3j 1 nd n nd nd
157. y sbyy iy 106 5 7 1 FDCALC statement functions The following code provides the statement functions used in establishing the tridiagonal matrix structure used in fdcalc refdif cg 1 3 sqrt p 1 3 q i Jj pv 1 3 p 1 3 v 1 3 v 1 3 bet i j 4 k i 1 j k i j dx k i 1 j k i j 2 ARA A Gi ANA k i j p i j u i j 2 dx k 141 3 k 1 3 2 p 1 1 3 p i j u i 1 j 2 u i j 2 dv i j cg i 1 3 u i 1 j sig i 1 3 cg i j u i 3 sig i j deltal dx v i 1 3 1 sig 1 1 3 1 v i j 1 sig i j 1 v it 1 j 1 sig i 1 j 1 v i j 1 sig i j 1 2 dy damp 1 3 2 ci cdamp cg i 1 j u i 1 3 cg i j u i j dy dy k 1 1 3 2 k 1 3 2 deltap i j al bl kb i k i j cpl 1 3 cg 1 1 3 u 1 1 3 cmp1x 1 dx kb 1 1 a0 k 1 1 3 1 cmp1x 1 0 cg 1 3 u 1 3 dv 1 3 sig i l j sig i 3 4 14 2 omeg cmp1x 0 1 b1 bet 1 3 u 1 1 3 u 1 3 sig i 1 j 1 4 omeg b1 cmp1x 0 1 3 u 1 1 3 u 1 3 dx v 1 1 3 1 Lev i j 1 v 1 1 3 1 v 4 3 1 4 dy sig i 1 3 k 41 3 L k i 3 1 cmp1x 2 b1 dy dy k 141 3 k 1 3 b1 bet i 3 dx 2 L dy dy l deltap i j dx 2 dy dy pv 1 1 3 1 2 pv 1 1 3 lpv i 1 j 1 sig i 1 j cmplx 1 0 omeg delta2 3 u i 1 j 1 u 1 3 2 sig 1 1 3 ci omeg a0 1 k itl 3 u 1 1 3

REF/DIF 1 - University of Delaware

Contents

Download Pdf Manuals

Related Search

Related Contents