C ********************************************************************** C * * C * SOFTWARE LICENSING * C * * C * This program is free software; you can redistribute * C * it and/or modify it under the terms of the GNU * C * General Public License as published by the Free * C * Software Foundation, either Version 2 of the * C * license, or (at your option) any later version. * C * * C * This program is distributed in the hope that it * C * will be useful, but without any warranty; without * C * even the implied warranty of merchantability or * C * fitness for a particular purpose. See the GNU * C * General Public License for more details. * C * * C * A copy of the GNU General Public License is * C * available at http://www.gnu.org/copyleft/gpl.html * C * or by writing to the Free Software Foundation, Inc.,* C * 59 Temple Place - Suite 330, Boston, MA 02111, USA. * C * * C ********************************************************************** SUBROUTINE WAVEDON C******************************PROGRAM WAVE*************************** C C PURPOSE: TO PREDICT SURFACE GRAVITY WAVES FROM THE WIND. C C VARIABLES : C THERE ARE SIX MAJOR ARRAYS USED TO DESCRIBE THE WAVE FIELD. C THEY ARE DIMENSIONED (IM,JM). THE SIX ARRAYS ARE: C C XMOM - THE X-COMPONENT OF MOMENTUM (units of m*s) C YMOM - THE Y-COMPONENT OF MOMENTUM C CO - THE MAGNITUDE OF THE MEAN PHASE SPEED (m/s) C SIG - THE WAVE HEIGHT STANDARD DEVIATION (m) C CS - THE COSINE OF THE MEAN WAVE DIRECTION C SN - THE SINE OF THE MEAN WAVE DIRECTION C C THE THREE OTHER ARRAYS ARE: C C WSX2D and WSY2D - THE X AND Y COMPONENTS OF THE WIND (m/s) c positive wsx2d is a westerly wind and c positive wsy2d is a southerly wind. C C HISTORY : C WRITTEN BY J.R. BENNETT, 1982, GLERL, ANN ARBOR,MI., BY C MODIFYING A PROGRAM WRITTEN BY M.A. DONELAN, NATIONAL WATER C RESEARCH INSTITUTE,BURLINGTON,ONTARIO. C C ********************************************************************* C MODIFIED TO INCORPORATE IN HYDROQUAL'S ECOM MODEL IN CURVILINEAR GRID C By QUAMRUL AHSAN, QA, October, 1998 C C HydroQual, Inc C 1 Lethbridge Plaza C Mahwah, NJ 07430 C e-mail qahsan@hydroqual.com C ********************************************************************* C INCLUDE 'comdeck' C C PHYSICAL CONSTANTS C SAVE REAL FSMWAD(IM,JM) !NG12202010 FOR WAD FOR WAVES DATA VK/0.41/ ! NG12172009 Von Karman to 0.41 from 0.4 chli DATA RHOAIR,RHOH2O,GAMMA,CBF/1.2233,1000.,0.028,0.0001/ c DATA RHOAIR,RHOH2O,GAMMA,CBF/1.2233,1000.,0.028,0.00005/ csfan change gamms to 0.028*2 CNG09192006 DATA RHOAIR,RHOH2O,GAMMA,CBF/1.2233,1000.,0.056,0.00005/ CNG09192006 Changed gamma back to 0.028, as ~3% (Lin, et.al. 2002 Continental Self Research pg 2677) CNG09192006 and changed formulation below for consistency CNG09192006 Also, a factor of 1/2 was brought in the bottom friction equation for consistency with linear CNG09192006 wave theory. Thus, the opening parameter CBF is also adjusted. CNG05152008 DATA RHOAIR,RHOH2O,GAMMA,CBF/1.2233,1000.,0.028,0.0002/ DATA RHOAIR,RHOH2O,GAMMA,CBF/1.2233,1000.,0.028,0.0002/ ! 2.8% is the standard !DATA RHOAIR,RHOH2O,GAMMA,CBF/1.2233,1000.,0.05,0.0002/ ! 5% of the total momentum to water goes to waves... upper Donelan limit C C COMPUTE CONSTANTS C c write(*,*) 'start wavedon, nr= ', nr CNG09192006 Changed gamma back to 0.028, as ~3% (Lin, et.al. 2002 Continental Self Research pg 2677) CNG09192006 and changed formulation below for consistency CNG09192006 FAC2=0.5*GAMMA*RHOAIR/(RHOH2O*G) CNG04232008 HOWEVER THERE IS A 0.5 IN THE ORIGINAL CODE, WITH A 0.028 FAC2=GAMMA*RHOAIR/(RHOH2O*G) CNG09192006 At sea level: g=9.780327[1+0.0053024sin^2(lat)-0.0000058sin^2(2lat)] G=GRAV ! NG04272010 Feed from ecom3d for consistency FAC2=GAMMA*RHOAIR/(RHOH2O*G) ! 2.8% of the wave momentum is in PI= 3.141593 TPI=6.283185 GTPI=G/TPI C-------------------------------------------------------------------------- C SET INITIAL CONDITIONS THIS SHOULD GO TO first.f C IF(NR.EQ.0)THEN DO J=1,JM DO I=1,IM XMOM(I,J)=0.0 YMOM(I,J)=0.0 CO(I,J)=0.0 CS(I,J)=1.0 SN(I,J)=0.0 SIG(I,J)=0.0 WN(I,J)=1. ENDDO ENDDO chli DO J=1,JM DO I=1,IM FSMW(I,J)=FSM(I,J) ENDDO ENDDO DO N=1,NUMQBC IC=IQC(N) JC=JQC(N) FSMW(IC,JC)=0.0 ENDDO c csjf let FSMW on the open boundary = 0 to void calculation on the c open boundary. Wave parameters on the open boundary are forced c DO N=1,NUMEBC FSMW(IETA(N),JETA(N))=0.0 ENDDO DO N=1,NUMECR ! SPECIAL BOUNDARY CORNER FSMW(ICHARD(N),JCHARD(N))=0.0 ENDDO c WRITE(IUNIT,*) 'print FSMW' CNG call print(fsm,fsm,im,jm,var,1.0,10,'SRC') ENDIF chli C C THIS PROGRAM USES TWO TIME STEPS. THE LARGE TIME STEP (DTI) C IS SPECIFIED BY THE USER AND CONTROLS THE FREQUENCY AT WHICH DATA C CAN BE PRINTED OR STORED FOR LATER ANALYSIS BY WOUTP, THE USER SUP- C PLIED OUTPUT ROUTINE. THE SMALL TIME STEP (DTT) IS CALCULATED BY C THE PROGRAM FROM THE FORMULA FOR NUMERICAL STABILITY OF THE UPWIND C FINITE DIFFERENCE SCHEME. THE NUMBER OF SMALL TIME STEPS PER LARGE C TIME STEP IS NDTT; IT IS STORED IN COMMON BLOCK /TPARM/. C INSIDE EACH SMALL TIME STEP THERE ARE TWO LOOPS OVER THE SPATIAL C COORDINATES I AND J. IN THE FIRST LOOP THE WAVE ADVECTION TERMS ARE C CALCULATED. IN THE SECOND LOOP THE WIND INPUT OF MOMENTUM IS ADDED C AND THE DIAGNOSTIC VARIABLES ARE CALCULATED. C C CALCULATE SMALL TIME STEP BASED ON CONSTANT WIND FROM MIDDLE OF TIME STEP C IF(CREADWBC) THEN ! Use read-in boundary condition WMAX=0. CNG09192006simplified NMAX=0 CNG09192006simplified DMAX=99999. CNG05312012 MOVE THIS BEFORE THE TIME STEP ASSIGNMENT TO ACCOMODATE CNG CALCULATION OF "EQUIVALENT-WIND" OF OBC SWELL FOR TIME STEP CNG12212010 MOVE NUMEBC ASSIGNMENT OUT OF THE 60 LOOP. FSMW ENSURES THESE ARE NOT CHANGED. DO N=1,NUMEBC CNG09192006 WHBFLX: WAVE HEIGHT STANDARD DEVIATION (Ho/4) IN METERS SIG(IETA(N),JETA(N))=WHBFLX(N) CNG09192006 WDBFLX: MOMENTUM FLUX DIRECTION, W.R.T. KSI1-DIRECTION AS CNG09192006 DONELAN AND SCHWAB'S CONVENTION. Wave is moving towards the CNG09192006 direction not from, counter clockwise is positive. CNG01312008 CS(IETA(N),JETA(N))=COSD(WDBFLX(N)*PI/180.) CNG01312008 SN(IETA(N),JETA(N))=SIND(WDBFLX(N)*PI/180.) CS(IETA(N),JETA(N))=COS(WDBFLX(N)*PI/180.) !NG01312008 standard FORTRAN does not have COSD SN(IETA(N),JETA(N))=SIN(WDBFLX(N)*PI/180.) !NG01312008 standard FORTRAN does not have SIND CNG09192006 Note that open boundary period WPBFLX or its phase speed, CO, CNG09192006 does not enter the calculations presently, so it is not required. WVPD(IETA(N),JETA(N))=WPBFLX(N) !NG06012009 Store it anyway for display of open BC WHNG=4*WHBFLX(N) WMAX=AMAX1(WMAX,0.0065*WHNG*WHNG*WHNG-0.23*WHNG*WHNG+4.*WHNG+2.)! CNG12212010 BASED ON AN UPPER BOUND OF BEAUFORT WMAX=AMAX1(WMAX,GTPI*WPBFLX(N))! CNG12212010 BASED ON LIMIT OF PHASE SPEED ENDDO CNG05312012 MOVE THIS BEFORE THE TIME STEP ASSIGNMENT TO ACCOMODATE ELSE ! assign from internal solution (open boundary condition). WMAX=0. DO N=1,NUMEBC CNG10292013 FOR NO-GRADIENT BOUNDARY CONDITION, TAKE BOUNDARY FROM PRIOR INNER SOLUTION: WHBFLX(N)=SIG(ICON(N),JCON(N)) WPBFLX(N)=WVPD(ICON(N),JCON(N)) WDBFLX(N)=90. - WVDR(ICON(N),JCON(N)) + - ANG(ICON(N),JCON(N))*180./PI ! KSI-1 DIR CNG10292013 AND THEN DO THE SAME AS ABOVE SIG(IETA(N),JETA(N))=WHBFLX(N) CS(IETA(N),JETA(N))=COS(WDBFLX(N)*PI/180.) !NG01312008 standard FORTRAN does not have COSD SN(IETA(N),JETA(N))=SIN(WDBFLX(N)*PI/180.) !NG01312008 standard FORTRAN does not have SIND WVPD(IETA(N),JETA(N))=WPBFLX(N) !NG06012009 Store it anyway for display of open BC WHNG=4*WHBFLX(N) WMAX=AMAX1(WMAX,0.0065*WHNG*WHNG*WHNG-0.23*WHNG*WHNG+4.*WHNG+2.)! CNG12212010 BASED ON AN UPPER BOUND OF BEAUFORT WMAX=AMAX1(WMAX,GTPI*WPBFLX(N))! CNG12212010 BASED ON LIMIT OF PHASE SPEED ENDDO ENDIF DS=99999. CNG if (nr.eq.0) then c call print(wsx2d,fsm,im,jm,var,1.0,10,'SRC') c call print(wsy2d,fsm,im,jm,var,1.0,10,'SRC') CNG endif c write(1007) wsx2d,wsy2d DO J=2,JMM1 DO I=2,IMM1 IF(FSMW(I,J).NE.0.) THEN IF(DT(I,J).LE.0) THEN WRITE(*,*)'WAVE ERROR',I,J,'DT<=0',DT(I,J) WRITE(IUPRT,*)'WAVE ERROR',I,J,'DT<=0',DT(I,J) CALL PRINTS(DRHOX,DRHOY,TRNU,TRNV) STOP ENDIF DS =AMIN1(DS,SQRT(H1(I,J)**2 + H2(I,J)**2)) WMAX=AMAX1(WMAX,SQRT(WSX2D(I,J)**2 + WSY2D(I,J)**2) + ,ABS(CO(I,J))) c write(*,*) 'wmax,win,co(i,j),wsx2d(i,j),wsy2d(i,j)', c & wmax,win,co(i,j),wsx2d(i,j),wsy2d(i,j),i,j CNG09192006simplified DS = AMIN1(DMAX,DS) CNG09192006simplified DMAX = DS ENDIF ENDDO ENDDO CNG Calculate number of timesteps: sqrt(2)*maxwindspeed*timeinterval/smallesthypotenuse + 2 NDTT=IFIX((1.414*WMAX*DTI*NWAVE)/DS)+2 CNG09192006simplified NDTT=MAX0(NMAX,NDTT) DTT=DTI*FLOAT(NWAVE)/NDTT CNG WRITE (*,*) NDTT,DTT DTTGAM=FAC2*DTT C C BEGIN LOOP OVER SMALL TIME STEPS C c write(1005,*) 'time=',time c write(1005,'(8f10.5)') (WHBFLX(N),n=1,numebc) c write(1005,'(8f10.5)') (WDBFLX(N),n=1,numebc) c write(*,*) 'win,ds',win,ds c write(*,*) 'wmax,dtt,ndtt,nwave,dti = ',wmax,dtt,ndtt,nwave,dti c cBased on sfan's (SIT) input, force sig and cs, sn on open boundaries c(hli,03/19/04) c DO I=1,IMM1 DO J=1,JMM1 FSMWAD(I,J)=DUMWAD(I+1,J)+DUMWAD(I,J)+DVMWAD(I,J+1)+DVMWAD(I,J) ! NG12202010 DEFINE WAVE WAD MASK HERE 0=DRY ENDDO ENDDO DO 60 NDT=1,NDTT C C FIRST LOOP OVER I,J: CALCULATE THE WAVE ADVECTION TERMS C c call smooth(sig,fsmw,im,jm,1) c call smooth(cs,fsmw,im,jm,1) c call smooth(sn,fsmw,im,jm,1) DO 70 I=2,IMM1 DO 70 J=2,JMM1 IF(FSMW(I,J)*FSMWAD(I,J).EQ.0) GOTO 70 ! NG12202010 CALCULATION ONLY FOR INTERNAL, NON-DRY, NON-BOUNDARY ELEMENTS DTTDS=DTT*0.25/(H1(I,J)*H2(I,J)) C C X-COMPONENT UPWIND SCHEME C IF(CS(I,J).GE.0.0) CNG03022007 For Thin Dams . XFLXD=H2(I,J)*((SIG(I,J)*CS(I,J))**2) CNG03022007 For Thin Dams . -H2(I-1,J)*((SIG(I-1,J)*CS(I-1,J))**2) . XFLXD= H2(I,J)*((SIG(I,J)*CS(I,J))**2) . +DUM(I,J)*DUMWAD(I,J)* . (-H2(I-1,J)*((SIG(I-1,J)*CS(I-1,J))**2)) IF(CS(I,J).LT.0.0) CNG03022007 For Thin Dams . XFLXD=H2(I+1,J)*((SIG(I+1,J)*CS(I+1,J))**2) CNG03022007 For Thin Dams . -H2(I,J)*((SIG(I,J)*CS(I,J))**2) . XFLXD=DUM(I+1,J)*DUMWAD(I+1,J)* . (H2(I+1,J)*((SIG(I+1,J)*CS(I+1,J))**2)) . -H2(I,J)*((SIG(I,J)*CS(I,J))**2) IF(SN(I,J).LT.0.0) CNG03022007 For Thin Dams . YFLXD=H1(I,J+1)*(CS(I,J+1)*SN(I,J+1)*SIG(I,J+1)**2) CNG03022007 For Thin Dams . -H1(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) . YFLXD=DVM(I,J+1)*DVMWAD(I,J+1)* . (H1(I,J+1)*(CS(I,J+1)*SN(I,J+1)*SIG(I,J+1)**2)) . -H1(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) IF(SN(I,J).GE.0.0) CNG03022007 For Thin Dams . YFLXD=H1(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) CNG03022007 For Thin Dams . -H1(I,J-1)*(CS(I,J-1)*SN(I,J-1)*SIG(I,J-1)**2) . YFLXD= H1(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) . +DVM(I,J)*DVMWAD(I,J)* . (-H1(I,J-1)*(CS(I,J-1)*SN(I,J-1)*SIG(I,J-1)**2)) SIP=SIG(I+1,J) SIM=SIG(I-1,J) H2F=H2 (I+1,J) H2B=H2 (I-1,J) CNG03022007 For Thin Dams IF(FSMW(I+1,J).EQ.0.) THEN CNG03022007 For Thin Dams SIP=2.*SIG(I,J)-SIG(I-1,J) CNG03022007 For Thin Dams H2F=2.*H2 (I,J)-H2 (I-1,J) IF(FSMW(I+1,J)*DUM(I+1,J)*DUMWAD(I+1,J).EQ.0.) THEN CNG IF(FSM(I+1,J)*DUM(I+1,J)*DUMWAD(I+1,J).EQ.0.) THEN SIP=2.*SIG(I,J)-DUM(I,J)*DUMWAD(I,J)*SIG(I-1,J) H2F=2.*H2 (I,J)-DUM(I,J)*DUMWAD(I,J)*H2 (I-1,J) ENDIF CNG03022007 For Thin Dams IF(FSMW(I-1,J).EQ.0.) THEN CNG03022007 For Thin Dams SIM=2.*SIG(I,J)-SIG(I+1,J) CNG03022007 For Thin Dams H2B=2.*H2 (I,J)-H2 (I+1,J) IF(FSMW(I-1,J)*DUM(I,J)*DUMWAD(I,J).EQ.0.) THEN CNG IF(FSM(I-1,J)*DUM(I,J)*DUMWAD(I,J).EQ.0.) THEN SIM=2.*SIG(I,J)-DUM(I+1,J)*DUMWAD(I+1,J)*SIG(I+1,J) H2B=2.*H2 (I,J)-DUM(I+1,J)*DUMWAD(I+1,J)*H2 (I+1,J) ENDIF XMOM(I,J)=XMOM(I,J)-DTTDS*(XFLXD+YFLXD+0.25*(H2F*SIP**2 . -H2B*SIM**2)) C C Y-COMPONENT UPWIND SCHEME C IF(CS(I,J).LT.0.0) CNG03022007 For Thin Dams . XFLXD=H2(I+1,J)*(CS(I+1,J)*SN(I+1,J)*SIG(I+1,J)**2) CNG03022007 For Thin Dams . -H2(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) . XFLXD=DUM(I+1,J)*DUMWAD(I+1,J)* . (H2(I+1,J)*(CS(I+1,J)*SN(I+1,J)*SIG(I+1,J)**2)) . -H2(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) IF(CS(I,J).GE.0.0) CNG03022007 For Thin Dams . XFLXD=H2(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) CNG03022007 For Thin Dams . -H2(I-1,J)*(CS(I-1,J)*SN(I-1,J)*SIG(I-1,J)**2) . XFLXD= H2(I,J)*(CS(I,J)*SN(I,J)*SIG(I,J)**2) . +DUM(I,J)*DUMWAD(I,J)* . (-H2(I-1,J)*(CS(I-1,J)*SN(I-1,J)*SIG(I-1,J)**2)) IF(SN(I,J).GE.0.0) CNG03022007 For Thin Dams . YFLXD=H1(I,J)*((SN(I,J)*SIG(I,J))**2) CNG03022007 For Thin Dams . -H1(I,J-1)*((SN(I,J-1)*SIG(I,J-1))**2) . YFLXD= H1(I,J)*((SN(I,J)*SIG(I,J))**2) . +DVM(I,J)*DVMWAD(I,J)* . (-H1(I,J-1)*((SN(I,J-1)*SIG(I,J-1))**2)) IF(SN(I,J).LT.0.0) CNG03022007 For Thin Dams . YFLXD=H1(I,J+1)*((SN(I,J+1)*SIG(I,J+1))**2) CNG03022007 For Thin Dams . -H1(I,J)*((SN(I,J)*SIG(I,J))**2) . YFLXD=DVM(I,J+1)*DVMWAD(I,J+1)* . (H1(I,J+1)*((SN(I,J+1)*SIG(I,J+1))**2)) . -H1(I,J)*((SN(I,J)*SIG(I,J))**2) SJP=SIG(I,J+1) SJM=SIG(I,J-1) H1F=H1 (I,J+1) H1B=H1 (I,J-1) CNG03022007 For Thin Dams IF(FSMW(I,J+1).EQ.0.) THEN IF(FSMW(I,J+1)*DVM(I,J+1)*DVMWAD(I,J+1).EQ.0.) THEN CNG IF(FSM(I,J+1)*DVM(I,J+1)*DVMWAD(I,J+1).EQ.0.) THEN SJP=2.*SIG(I,J)-DVM(I,J)*DVMWAD(I,J)*SIG(I,J-1) H1F=2.*H1 (I,J)-DVM(I,J)*DVMWAD(I,J)*H1 (I,J-1) ENDIF CNG03022007 For Thin Dams IF(FSMW(I,J-1).EQ.0.) THEN IF(FSMW(I,J-1)*DVM(I,J)*DVMWAD(I,J).EQ.0.) THEN CNG IF(FSM(I,J-1)*DVM(I,J)*DVMWAD(I,J).EQ.0.) THEN SJM=2.*SIG(I,J)-DVM(I,J+1)*DVMWAD(I,J+1)*SIG(I,J+1) H1B=2.*H1 (I,J)-DVM(I,J+1)*DVMWAD(I,J+1)*H1 (I,J+1) ENDIF YMOM(I,J)=YMOM(I,J)-DTTDS*(XFLXD+YFLXD+0.25*(H1F*SJP**2 . -H1B*SJM**2)) C 70 CONTINUE c call smooth(xmon,fsmw,im,jm,1) c call smooth(ymon,fsmw,im,jm,1) C C SECOND LOOP OVER I,J: CALCULATE THE WIND INPUT OF MOMENTUM C AND COMPUTE THE DIAGNOSTIC VARIABLES (CS,SN,S,C) C DO 90 J=2,JMM1 DO 90 I=2,IMM1 IF(FSMW(I,J)*FSMWAD(I,J).EQ.0) GOTO 90 ! NG12202010 CALCULATION ONLY FOR INTERNAL, NON-DRY, NON-BOUNDARY ELEMENTS C C CALCULATE THE VERTICAL MOMENTUM FLUX IN THE 2 DIRECTIONS, C THE WIND AND THE WAVE FIELD CNG09192006 NOTE WIND IS ROTATED TO KSI-1 AND KSI-2 DIRECTIONS FROM E-W AND N-S CNG09192006 TO FORM THE REQUIRED DOT PRODUCT OF WAVES.WIND, AS IN THE ORIGINAL CNG09192006 SCHWAB-BENNETT CODE. C c wsx2d and wsy2d are E-W and N-S components. positive u is a westerly c wind and positive v is a southerly wind. CNG04232008 UPDATED FORMULATION TO BRING IT CLOSER TO THE SCHWAB ORIGINAL GLERL-51 CNG04232008 ALSO TRIED THE GLERL-51 FORMULATION (NOT CODE), BUT IT CRASHES ON HIGH WINDS UWIND=(WSX2D(I,J)*COS(ANG(I,J))+ + WSY2D(I,J)*SIN(ANG(I,J))) VWIND=(-WSX2D(I,J)*SIN(ANG(I,J))+ + WSY2D(I,J)*COS(ANG(I,J))) WSPDSQ=UWIND*UWIND+VWIND*VWIND WNDSPD=SQRT(WSPDSQ)+1.E-8 ! NG12172009 increased from 1.E-10 CSWIND=UWIND/WNDSPD SNWIND=VWIND/WNDSPD CTHE1=CSWIND*CS(I,J)+SNWIND*SN(I,J) B=CO(I,J)*.83/WNDSPD SG=AMAX1(0.005,SIG(I,J)*ABS(CTHE1)) ! NG The magnitude of covariance enters the drag coeff... ! so that if wind is normal to waves, minimal drag. Z0WAVE=AMIN1(SG,49.999)/5. ! NG12172009 Original surface roughness Z0 in GLERL CNG12172009 TRIED, BUT GAMMA SEEMS TO NEED ADJUSTMENT TO 6-7%... Z0WAVE=1.38E-4*4.*SIG(I,J)*(WNDSPD/CO(I,J))**2.66 +1.E-5 ! NG12172009 Z0 from Donelan 1990, a function of inverse wave age DRAG1=(VK/ALOG(10./Z0WAVE))**2 CNG03232007 NG's TAKE OF SCHWAB'S PUBLISHED METHOD (GLERL-51) BASED ON FORMULA PAGE 3. CRASHES. CNG03232007 AA=SQRT(1+B*B-2*B*CTHE1) CNG03232007 WI1=DRAG1*AA/WNDSPD CNG03232007 WI2=-B*WI1 AA=1.-B*CTHE1 WI1=DRAG1*AA*ABS(AA) AA=CTHE1-B ! AS PER ORIGINAL WI2=DRAG1*AA*ABS(AA) CNG09192006 VFLXA=WI1*CSWIND+WI2*CS(I,J) CNG09192006 VFLYA=WI1*SNWIND+WI2*SN(I,J) CNG09192006 INCLUDE * WSPSSQ HERE, NOT IN DTFAC CNG09192006 ADD SKIN FRICTION CONTRIBUTION DRAGS = 0.7E-03 ! SKIN DRAG COEFFICIENT VFLXA = DRAGS*WNDSPD*UWIND + + (WI1*CSWIND+WI2*CS(I,J)) * WSPDSQ VFLYA = DRAGS*WNDSPD*VWIND + + (WI1*SNWIND+WI2*SN(I,J)) * WSPDSQ C C COMPUTE BOTTOM FRICTION LOSS FOR SHALLOW WATER USING STANDARD DEVIATION C OF WAVEHEIGHT FOR AMPLITUDE AND WAVENUMBER OF PEAK ENERGY FREQUENCY C IF(DT(I,J)*WN(I,J).LT.10.) THEN CNG09192006 UBOT=SIG(I,J)*WN(I,J)*CO(I,J)/SINH(WN(I,J)*DT(I,J)) CNG09192006 A 1/2 FACTOR WAS BROUGHT IN THE EQUATION FOR CONSISTENCY WITH LINEAR CNG09192006 WAVE THEORY. THUS, THE OPENING PARAMETER CBF IS ALSO ADJUSTED. CNG09192006 10.0 is just an arbitrary limit for when 1/sinh->0 UBOT=SIG(I,J)*WN(I,J)*CO(I,J)/2./(SINH(WN(I,J)*DT(I,J))+1.E-10) ELSE UBOT=0. ENDIF CNG09192006 DTFAC=DTTGAM*WSPDSQ DTFAC=DTTGAM XMOM(I,J)=XMOM(I,J)+DTFAC*VFLXA-DTT*CBF*UBOT*UBOT*CS(I,J) YMOM(I,J)=YMOM(I,J)+DTFAC*VFLYA-DTT*CBF*UBOT*UBOT*SN(I,J) 90 CONTINUE C c using laplacian smoother every time step c call smooth(xmon,fsmw,im,jm,1) c call smooth(ymon,fsmw,im,jm,1) DO 91 J=2,JMM1 DO 91 I=2,IMM1 IF(FSMW(I,J)*FSMWAD(I,J).EQ.0) GOTO 91 ! NG12202010 CALCULATION ONLY FOR INTERNAL, NON-DRY, NON-BOUNDARY ELEMENTS C C COMPUTE THE NEW COSINES AND SINES; AND COMPUTE THE PHASE C SPEED (C) AND WAVE HEIGHT STANDARD DEVIATION (S). CNG09192006 This uses Hasselmann et al. (1973) JONSWAP spectrum: CNG09192006 Joint North Sea Wave Observation Project CM=SQRT(XMOM(I,J)**2+YMOM(I,J)**2)+1.E-10 CS(I,J)=XMOM(I,J)/CM SN(I,J)=YMOM(I,J)/CM CNG12172009 BIG ERROR: UWIND and VWIND were not defined as (I,J) matrices? Added... UWIND=(WSX2D(I,J)*COS(ANG(I,J))+ + WSY2D(I,J)*SIN(ANG(I,J))) VWIND=(-WSX2D(I,J)*SIN(ANG(I,J))+ + WSY2D(I,J)*COS(ANG(I,J))) CNG09212006 UU is the wind in the direction of the wave momentum flux CNG09212006 (simply rotated). Then, if wind opposing the wave, do 150. UU=UWIND*CS(I,J)+VWIND*SN(I,J) IF (UU.LE.0.0) GO TO 150 FM=0.01788735*(((UU*UU)/(CM*CM*CM+1.E-10))**0.14285714) IF((UU*FM).GT.0.760545) GOTO 160 150 FM=(CM*14343.09)**(-.3333333333) 160 CONTINUE CNG09212006 THROW THIS AFTER IF STATEMENT:SIG(I,J)=SQRT(CM * GTPI/FM) CNG09192006 WN(I,J)=WNUM(FM,DT(I,J)) CNG12022008 WN(I,J)=WNUM(FM,DT(I,J),G) ! NG09192006 INCLUDE G AS INPUT WN(I,J)=WNUMYOU(FM,DT(I,J),G) ! NG12022008 USED BETTER SOLUTION IF(DT(I,J)*WN(I,J).LT.3.14) THEN CNG09192006 LINEAR WAVE THEORY. STRAIGHT FROM THE DISPERSION EQUATION. CNG04272010 w^2=gktanh(kh) and c=w/k. CNG09192006 3.14 is just the limit when tanh->1 CO(I,J)=SQRT(G*TANH(WN(I,J)*DT(I,J))/(WN(I,J)+1.E-10)) ELSE CO(I,J)=GTPI/(FM+1.E-10) ! NG09192006 classic limit w^2=gk, g/(2pif)=g/w=w/k ENDIF CNG09192006 SIG(I,J)=SQRT(CM * GTPI/FM) CNG12202010 SIG(I,J)=SQRT(CM * CO(I,J)) ! THROWN HERE, AFTER CO IS CALCULATED SIG(I,J)=AMAX1(0.,SQRT(CM * CO(I,J))) ! DO NOT ALLOW NEGATIVE HEIGHTS C 91 CONTINUE c call smooth(sig,fsmw,im,jm,1) c call smooth(cs,fsmw,im,jm,1) c call smooth(sn,fsmw,im,jm,1) C C END LOOP OVER SMALL TIME STEPS C c cBased on sfan's (SIT) input, force sig and cs, sn on open boundaries c(hli,03/19/04) c c DO N=1,NUMEBC c SIG(IETA(N),JETA(N))=WHBFLX(N) c CS(IETA(N),JETA(N))=COSD(WDBFLX(N)) c SN(IETA(N),JETA(N))=SIND(WDBFLX(N)) c FM=1/(WPBFLX(N)+1.E-10) c WN(IETA(N),JETA(N))=WNUM(FM,DT(IETA(N),JETA(N))) c IF(DT(IETA(N),JETA(N))*WN(IETA(N),JETA(N)).LT.3.14) THEN c CO(IETA(N),JETA(N))=SQRT(G*TANH(WN(IETA(N),JETA(N)) c & *DT(IETA(N),JETA(N)))/WN(IETA(N),JETA(N))) c ELSE c CO(IETA(N),JETA(N))=GTPI/FM c ENDIF c XMOM(IETA(N),JETA(N))=0.0 c YMOM(IETA(N),JETA(N))=0.0 c ENDDO C END 60 CONTINUE C C CALCULATE WAVE PARAMETERS (Significant Wave Height, Wave Period C and Wave Direction) C DO 21 J=1,JM DO 21 I=1,IM CNG03092007PreserveOpenBC IF(FSMW(I,J).EQ.0.) GO TO 21 CNG IF(FSM(I,J)*FSMWAD(I,J).EQ.0) THEN !NG12212010 IF(FSM(I,J).EQ.0) THEN !NG12212010 WVHT(I,J)=-99999.0 ! mask out exterior points and currently dry points (but not OBC, unless dry!) GOTO 21 ENDIF ! 2 below worked but took out WAD and Open BC points in output at least... ! WVHT(I,J)=-999.0 ! IF(FSMW(I,J)*FSMWAD(I,J).EQ.0) GOTO 21 WVHT(I,J)=4.*SIG(I,J) 21 CONTINUE CPL DO J=1,JM DO I=1,IM CNG IF(FSMW(I,J)*FSMWAD(I,J).NE.0) THEN ! NG12212010 if interior non-boundary, or wet (now), then IF(FSMW(I,J).NE.0) THEN ! NG12212010 if interior non-boundary, or wet (now), then CNG10202008 HBREAK=0.78*DT(I,J) ! Outdated McCowan, 1894 expression HBREAK=0.8261*DT(I,J) ! NG10202008 Use Longuet-Higgins and Fenton, 1974 !WVHT(I,J)=AMIN1(WVHT(I,J),HBREAK) WVHT(I,J)=AMIN1(ABS(WVHT(I,J)),HBREAK) ! NG12212010 FOR PREVIOUSLY DRY MASKED-OUT POINTS THAT ARE NOW WET USE ABS SIG(I,J)=WVHT(I,J)/4. ! NG09212006 update SIG too! Be mindful of the 21 loop! END IF END DO END DO C C SAVE WAVE PERIODS C DO 22 J=1,JM DO 22 I=1,IM IF(FSMW(I,J).EQ.0) GOTO 22 ! NG12202010 CALCULATION ONLY FOR INTERNAL, NON-BOUNDARY ELEMENTS CNG IF(FSMW(I,J).EQ.0) GOTO 22 ! NG12202010 CALCULATION ONLY FOR INTERNAL, NON-DRY, NON-BOUNDARY ELEMENTS CNG04272010 THIS WAS DEEP WATER APPROXIMATION: WVPD(I,J)=TPI*CO(I,J)/GRAV ! NG09192006 CNG04272010 INSTEAD USE: c=w/k => w=ck => 2pi/T=ck => T=2pi/(ck) CNG12202010 WVPD(I,J)=TPI/(CO(I,J)*WN(I,J)+1.E-10) WVPD(I,J)=AMAX1(0.,TPI/(CO(I,J)*WN(I,J)+1.E-10)) ! DO NOT ALLOW NEGATIVE PERIODS IF (FSMWAD(I,J).EQ.0) WVPD(I,J)=0. ! ZERO OUT PERIOD OF DRY POINTS FOR OUTPUT 22 CONTINUE C C SAVE WAVE DIRECTIONS C DO 23 J=1,JM DO 23 I=1,IM CNG03092007forOpenBoundaryPlotting IF(FSMW(I,J).EQ.0.) GO TO 23 CNG03092007forOpenBoundaryPlotting UU=XMOM(I,J) CNG03092007forOpenBoundaryPlotting VV=YMOM(I,J) CNG IF(FSM(I,J)*FSMWAD(I,J).EQ.0) GOTO 23 ! NG12202010 CALCULATION ONLY FOR ALL WET ELEMENTS, INCLUDING OBC IF(FSM(I,J).EQ.0) GOTO 23 ! NG12202010 CALCULATION FOR ALL WET AND DRY ELEMENTS, INCLUDING OBC UU=CS(I,J) VV=SN(I,J) WVDR(I,J) = 0. IF(UU .EQ. 0. .AND. VV .EQ. 0.) GOTO 23 C Donelan and Schwab's convention of wave direction w.r.t. zeta1 axis c wave is moving towards the direction not from, counter clock is positive WVDR(I,J)=ATAN2(VV,UU)*180./PI ! w.r.t. zeta1-direction CNG IF(WVDR(I,J) .LT. 0.) WVDR(I,J) = WVDR(I,J) + 360. CNG IF(WVDR(I,J) .GT. 360.) WVDR(I,J) = WVDR(I,J) - 360. c Direction with respect to east (x-direction) WVDR(I,J) = WVDR(I,J) + ANG(I,J)*180./PI C IN stress.f and in SMB Model WAVE DIR IS IN w.r.t North Clock Wise WVDR(I,J)= 90.-WVDR(I,J) WVDR(I,J) = AMOD(WVDR(I,J),360.) IF(WVDR(I,J).LE.0.) WVDR(I,J) = WVDR(I,J) + 360. c final wave direction is w.r.t North clock wise and towards (oceanographic) 23 CONTINUE C CNG03072008 NR=NR+1 NR=1 ! stop counter C HARDWIRE for Checking Purposes************ c WRITE(1001)TIME c WRITE(1001)WVHT c WRITE(1001)WVPD c WRITE(1001)WVDR c write(1001) time,wvht(720,79),wvht(574,89),wvht(320,111) c & ,wvht(144,26),wvpd(720,79),wvpd(574,89),wvpd(320,111) c & ,wvpd(144,26),wvdr(720,79),wvdr(574,89) c & ,wvdr(320,111),wvdr(144,26) C 540 FORMAT(1X, 'DEPTH RELATIVE TO MEAN WATER LEVEL') c write(*,*) 'end wavedon, nr= ', nr c write(1008) wsx2d,wsy2d RETURN END FUNCTION WNUM(F,DD,G) ! NG09192006 Include G as input CNG09192006 FUNCTION WNUM(F,DD) CNG I HAVE CHECKED THIS VERSUS CETN-1-17 6/85 Coastal Engineering CNG Technical Note; Direct Methods for Calculating Wavelength. CNG Provides 0.1% accuracy. C C APPROXIMATE SOLUTION OF WAVE DISPERSION EQUATION C C F = FREQUENCY (HZ) C DD = DEPTH (M) C WNUM = WAVENUMBER (1/M) C C REFERENCE: HUNT, J.N. 1979. DIRECT SOLUTION OF WAVE DISPERSION EQUATION. C JWPCOD, ASCE, 105(WW4):457-459. CNG ACTUALLY, AS REFINED BY CHEN AND THOMPSON, 1985. Iterative and CNG Pade solution for the water wave dispersion relation CERC-85-4. C DATA TPI/6.283185/ CNG09192006 DATA G/9.806/ Y=DD*(TPI*F)**2/G X=Y*(Y+1./(1.+Y*(0.6522+Y*(0.4622+Y*Y*(0.0864+Y*0.0675))))) WNUM=SQRT(X)/(DD+1.E-10) RETURN END FUNCTION WNUMYOU(F,DD,G) ! NG12022009 CNG CNG APPROXIMATE SOLUTION OF WAVE DISPERSION EQUATION CNG CNG F = FREQUENCY (HZ) CNG DD = DEPTH (M) CNG WNUMYOU = WAVENUMBER (1/M) CNG CNG Reference: Zai-Jin You, 2008 CNG "A close approximation of wave dispersion relation for CNG direct calculation of wavelength in any coastal water depth." CNG Equations 3 and 2 (or 3 and 4). CNG Applied Ocean Research CNG Provides 0.01% (0.1%) maximum relative error. CNG CNG Note that the YOU (2) equation can be analytically rewritten BUT CNG the resulting simpler expression becomes numerically indefinite CNG because it involves ~inf/~inf expressions like cosh*sinh. CNG DATA TPI/6.283185/ HK0=DD*TPI*F*TPI*F/G X0=SQRT(HK0)*(1.+HK0/6.+HK0*HK0/30.) ! YOU (3) Y=X0*( (HK0+X0/COSH(X0)*X0/COSH(X0)) / + (X0*TANH(X0)+X0/COSH(X0)*X0/COSH(X0)) ) ! YOU (2) NEWTON-RAPHSON CNG Y=HK0/TANH(X0) ! YOU (4), SOMEWHAT FASTER ~ HUNT IN ACCURACY AND SPEED WNUMYOU=Y/(DD+1.E-10) RETURN END