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 STRESS C CALC. BOTTOM SHEAR STRESS DUE TO WIND WAVES AND CURRENTS C USES GRANT/MADSEN/GLENN WAVE MODEL C REVISION DATE: July 2, 1999 c This Subroutine uses the CDCALC subroutine rovided by Rich Signell c and introduced in this code by Quamrul Ahsan, QA on 7/2/99 c The routine Calculates effective CD due to current and waves using c the lowest layer velocity from ECOM and wave information c from WAVESMB or WAVEDON Model (PERIOD = mean wave period of energy c density spectrum, UBOT = rms value of the maxima of the orbital c velocity near the bottom (m/s)) c This routine is a simplified form of Grant & Madsen (1979) c which assumes colinear waves and currents. This is described c in "Effect of wave-current interaction on steady wind-driven c circulation in narrow, shallow embayments", Signell et al, c JGR, 1990, v 95, n c6, 9671--9678. C************************************************************** CNG 04272010 Changed to make more wavedon like, and for W&D. c Common block REFLE needed for Grant & Madsen routines CNG04272010 Common /REFLE/ ZR, Z0C, USTARC, USTACW, USTARW Common /REFLE/ ZR, Z0C, USTARC, USTACW, USTARW, PHICOL,COSPHICOL ! NG04272010 Include wave-current angle c USTARW = the friction velocity due to waves only c USTACW = the friction velocity due to combined wave and currents Real OMEGA, KBR INCLUDE 'comdeck' REAL FSMWSTRS(IM,JM) !NG06092008 CNG12212010 Not used Data IFIRST /0/ PI = 3.141593 TPI= 6.283185 HPI= 1.570796 chli DO J=1,JM DO I=1,IM FSMWSTRS(I,J)=FSM(I,J) ENDDO ENDDO DO J=1,JMM1 DO I=1,IMM1 IF(DUMWAD(I,J)+DUMWAD(I+1,J)+DVMWAD(I,J)+DVMWAD(I,J+1).EQ.0) . FSMWSTRS(I,J)=0.0 ! NG12212010 FOR WAD ENDDO ENDDO DO N=1,NUMQBC IC=IQC(N) JC=JQC(N) FSMWSTRS(IC,JC)=0.0 ENDDO CNG06082012 DO N=1,NUMEBC CNG06082012 FSMWSTRS(IETA(N),JETA(N))=0.0 ! NG12212010 ALSO OPEN BOUNDARIES CNG06082012 ENDDO chli C C CALC. TOTAL BOTOM SHEAR STRESS C DO 30 J=2,JMM1 DO 30 I=2,IMM1 IF (FSMWSTRS(I,J).EQ.0.0) GOTO 30 C C NO SHEAR STRESS IF D < 10 cm C CNG IF (DT(I,J).LT.0.10) THEN ! CNG04272010 W&D CNG11012013 IF (DT(I,J).LE.AMAX1(0.10,WETMIN)) THEN ! CNG04272010 W&D !!!!!!!!!!!crm IF (DT(I,J)-WETMIN.LE.0.10) THEN ! CNG11012013 W&D !!!!!!!!!!!crm TAU(I,J,KB)=0.0 !!!!!!!!!!!crm GOTO 30 !!!!!!!!!!!crm ENDIF IF(FSMW(I,J).EQ.0.0.OR.WVHT(I,J).LT.0.30.OR. . DT(I,J).LT.1.0) THEN !!!!!!!RMarsooli crm IF(WVPD(I,J).LE.0.1.OR.WVHT(I,J).LE.0.01) THEN CNG04272010 FOR SED TRAN AS PER PRAVI'S NOTE IN BOTTAU TAU(I,J,KB)=10000.*CBC(I,J)*QBAR(I,J)*QBAR(I,J) TAU(I,J,KB)=10000.*CBC(I,J)/VARBF(I,J)*QBAR(I,J)*QBAR(I,J) GOTO 30 ENDIF C C CALC. WAVE FREQUENCY (1/s) C !!!!!!!!!!!!! SMALLVALL added by RMarsooli OMEGA=TPI/(WVPD(I,J)+SMALLVAL) ! CNG04272010 MOVE IT HERE CNG04272010C CNG04272010C CALC. WAVE NUMBER (1/m) CNG04272010 WILL USE WNUMYOU DISPERSION SOLVER FROM NOW ON, CONSISTENT WITH WAVEDON WAVEK=WNUMYOU(OMEGA/TPI,DT(I,J),GRAV) CNG04272010C ASSUME OMEGAR=OMEGA CNG04272010C CALC. SHALLOW WATER WAVE SPEED (m/s) CNG04272010C CNG04272010 C0=SQRT(GRAV*DT(I,J)) CNG04272010C CNG04272010C ESTIMATE WAVE NO. CNG04272010C CNG04272010 WAVEK0=TPI/(C0*WVPD(I,J)) CNG04272010C CNG04272010C ITERATE FOR WAVE NO. CNG04272010C CNG04272010 TRAT=TPI*TPI/(GRAV*WVPD(I,J)*WVPD(I,J)) CNG04272010C CNG04272010C USE NEWTON'S METHOD CNG04272010C CNG04272010 DO 40 N=1,20 CNG04272010 TANHKH=TANH(WAVEK0*DT(I,J)) CNG04272010 FKH=WAVEK0*TANHKH-TRAT CNG04272010 FPKH=TANHKH+WAVEK0*DT(I,J)*(1.-TANHKH*TANHKH) CNG04272010 WAVEK=WAVEK0-FKH/FPKH CNG04272010C CNG04272010C CHECK FOR CONVERGENCE CNG04272010C CNG04272010 DELTA=ABS(WAVEK/WAVEK0-1.) CNG04272010 IF (DELTA.LE.0.0001) GOTO 45 CNG04272010 WAVEK0=WAVEK CNG04272010 40 CONTINUE CNG04272010 45 CONTINUE C C CALC. BOTTOM ORBITAL AMPLITUDE (m) C ABM=WVHT(I,J)/(2.*SINH(WAVEK*DT(I,J))) + 0.00001 C C CALC. BOTTOM ORBITAL VELOCITY (m/s) C UB2=ABM*OMEGA + 0.00001 C C HEIGHT OF VELOCITY ABOVE BOTTOM C ZR=DZ(KBM1)*DT(I,J)/2. IF(TOR.EQ.'BAROTROPIC') ZR=0.5*DT(I,J) ! FOR BAROTROPIC THINK OF IT AS MID-POINT? ZR=AMIN1(ZR,0.4) !RMarsooli CRM KBR=30.*Z0B KBR=30.*Z0B(I,J) C C GET CURRENTS FROM MODEL AT HEIGHT ZR ABOVE BED & CONVERT TO CGS UNITS C crm URSPD = SQRT(UBAR(I,J,KBM1)*UBAR(I,J,KBM1)+ crm + VBAR(I,J,KBM1)*VBAR(I,J,KBM1)) * 100. URSPD0 = QBAR(I,J) + SMALLVAL !RMarsooli, 2016 URSPD = URSPD0*100.0 !RMarsooli, 2016 CNG04272010 WE SHOULD ADD, FOR NON-COLINEAR WAVES: crm PHICOL=((-(WVDR(I,J)-90.)-ANG(I,J))*PI/180.- crm . ATAN2(VBAR(I,J,KBM1),UBAR(I,J,KBM1))) C RMarsooli, 2016 COSPHICOL=0.0 c IF(WAVEDYN.EQ.'MELLOR') THEN c VEL_CX_TEMP = 0.5*(UA(I,J)+UA(I+1,J)) c VEL_CY_TEMP = 0.5*(VA(I,J)+VA(I,J+1)) c VEL_CX = VEL_CX_TEMP*COSANG(I,J)-VEL_CY_TEMP*SINANG(I,J) !w.r.t west-east axis c VEL_CY = VEL_CX_TEMP*SINANG(I,J)+VEL_CY_TEMP*COSANG(I,J) !w.r.t south-north axis c VEL_MAG=SQRT(VEL_CX**2+VEL_CY**2)+SMALLVAL c IF(VEL_CX.eq.0.0) THEN c ANGLE_C=0.5*3.141593*SIGN(1.0,VEL_CY) c ELSE c ANGLE_C=ATAN2(VEL_CY,VEL_CX) c ENDIF c ANGLE_W=thtav(I,J) c PHICOL=ANGLE_W-ANGLE_C c COSPHICOL=COS(PHICOL) c ENDIF C C CALC. SHEAR VELOCITY DUE TO WAVES & CURRENTS C NEAR BOTTOM CURRENT VELOCITY C SET MINIMUM CURRENT VELOCITY TO 1 mm/s C CRM URSPD=AMAX1(URSPD,0.1) C C CALL GRANT & MADSEN ROUTINE (CGS UNITS) C TO GET EFFECTIVE DRAG COEFFICIENT DUE TO WAVES AND CURRENTS C CALL DRGCO(UB2*100,KBR*100.,OMEGA,URSPD,0.,CDE,Z0) CNG CALL CDNEW(UB2*100,KBR*100.,OMEGA,URSPD,CDE,Z0) C C DO NOT LET DRAG COEFFICIENT BE LESS THAN CBC C CNG06092008 CBC(I,J) = AMAX1(CDE,BFRIC) C IF (I.EQ.23.AND.J.EQ.26) C . write (*,*) 'TIME:',TIME,'CDE:',CDE,'CBC',CBC(I,J) ! NG LITTLE CHECK FOR DRAG ! IF(I.EQ.93.and.J.EQ.15.AND.MOD(INTX,10).EQ.0.) THEN ! write (1003,'(5(A,F12.9),4(A,F7.3))') ! . 'TIME:',TIME,' CDE:',CDE,' CBC',CBC(I,J), ! . ' UCUR:',URSPD/100.,' UWAV:',UB2, ! . ' DT:',DT(I,J),' WVPD:',WVPD(I,J), ! . ' WVHT:',WVHT(I,J), 'WVN:',WAVEK ! NG LITTLE CHECK FOR DRAG ! ENDIF CBC(I,J)=AMIN1(AMAX1(CDE,CBC(I,J)),0.5) ! NG05172010 UPDATE, BUT DO NOT ALLOW CBC TO EXCEED 0.5 !UNITY. TAU(I,J,KB)=AMAX1(USTARC,USTACW) TAU(I,J,KB)=TAU(I,J,KB)*TAU(I,J,KB) 30 CONTINUE C C FORCING TAU VALUE = 0 AT CELLS CONNECTED BY RIVERS (IC,JC) C By Quamrul QA 10/28/98 C DO 120 N=1,NUMQBC IC=IQC(N) JC=JQC(N) TAU(IC,JC,KB) = 0.0 120 CONTINUE RETURN END C***************************************************************** Subroutine DRGCO(UB,KBR,OMEGA,UR,DIA,CD,Z0) CNG I don't think we should save FW ! SAVE FW C ****************************************** C SUBROUTINE DRGCO DRAG COEFFICIENT C W. D. GRANT W.H.O.I. C GENERATED FEB 14, 1984 C LATEST REVISION APR 1, 1984 C Determines the drag coefficient for mean C flow. C SUBROUTINES CALLED: C WFRIC C MOVKBR C ****************************************** c If you want movable bed effects, specify KBR=100. and call c this routine with a realistic sand diameter (DIA). c If you DON'T want movable bed effects, call this routine c with DIA = 0. Real KBRAB, KBR CNG04272010 Common /REFLE/ ZR,z0c,ustarc,ustacw,ustarw Common /REFLE/ ZR,z0c,ustarc,ustacw,ustarw,PHICOL,COSPHICOL ! CNG04272010 Included wave-current angle If (KBR.EQ.100.0) Then Call MOVKBR(DIA,UB,OMEGA,KBRAB) KBR = KBRAB * (UB/OMEGA) Else CNG04272010 KBRAB is the Nikuradse bed roughness Ks over the CNG04272010 semi-excursion distance, UBOT(WAVE)/OMEGA CNG04272010 where UBOT(WAVE)=pi*Ho/(Tsinh(kh)). CNG04272010 In scientific literature, it is usually noted as (ks/A). KBRAB = KBR / (UB/OMEGA) End If Call WFRIC(KBRAB,FW) C WRITE(6,*) 'KBR ',KBR,' KBRAB ',KBRAB,' UB ',UB,' OMEGA ', C * OMEGA, C * ' FW ',FW USTARW = SQRT(0.5*FW) * UB C ------------------------------------------ C INITIAL DRAG COEF ESTIMATE DELTAW = 2.0 * 0.4 * USTARW / OMEGA Z0 = KBR / 30.0 Z0C = Z0 * ((DELTAW/Z0)**0.8) CPL PRAVI & QA 09/16/99 CD = (0.4/ALOG(ZR/Z0C)) ** 2.0 CD = (0.4/ALOG(ZR*100./Z0C)) ** 2.0 !ZR converted to cm C WRITE(6,*) 'INITIAL ESTIMATE OF CD ',CD C WRITE(6,*) 'Z0 ',Z0 C WRITE(6,1000) C1000 FORMAT(//,' I CD USTARC UAUB ALPHA USTACW DELTAW ', C * ' B Z0C',/) c write(1002,*)'SUBRO DRGCO = ',Z0 CNG Return CNG Entry CDNEW(UB,KBR,OMEGA,UR,CD,Z0) C ------------------------------------------ C UPDATE CD C WRITE(4001,*) 'INSIDE CDNEW---- ZR ',ZR C WRITE(4001,*) 'FW ',FW C WRITE(4001,*) 'Z0 ',Z0 C WRITE(4001,*) 'UR ',UR C WRITE(4001,*) 'UB ',UB,' KBR ',KBR,' OMEGA ',OMEGA C ----- INSERT SIMPLE VERSION OF DRAG COEFF. SUBROUTINE ------- c write(1002,*)'ENTRY CDNEW = ',Z0 Do 100 I = 1, 10 CD0 = CD USTARC = SQRT(CD) * UR C UAUB= (CD/(0.64*FW))*(UR/UB)**2.0 C ALPHA= (1.0+UAUB)**2.0 C USTACW= ((0.5*FW*ALPHA)**0.5)*UB USTAW2 = 0.5 * FW * UB ** 2 USTAW = SQRT(USTAW2) USTACW = SQRT(USTAW2+USTARC**2) CNG04272010 ADDED ANGLE FOR NON-COLINEAR WAVES CNG04272010 TOO BIG USTACW = SQRT(USTAW2+USTARC**2+(2.*USTAW*USTARC*COS(PHICOL))) USTACW = SQRT(USTAW2+USTARC**2+(2.*USTAW*USTARC*COSPHICOL)) !RMarsooli !write (*,*) USTAW,USTARC,USTACW,PHICOL DELTAW = 2.0 * 0.4 * USTACW / OMEGA B = 1.0 - (USTARC/USTACW) CNG04272010 NOW B<0 IS POSSIBLE CNG04272010 If (B.LT.0) Then CNG04272010 Write (6,*) 'B LESS THAN 0' CNG04272010 Stop CNG04272010 End If Z0C = Z0 * ((DELTAW/Z0)**B) CPL PRAVI & QA 09/16/99 CD = (0.4/ALOG(ZR/Z0C)) ** 2.0 CD = (0.4/ALOG(ZR*100./Z0C)) ** 2.0 !ZR converted to cm DIFF = ABS(CD-CD0) If (DIFF.LE.1.E-07) Go To 110 C WRITE(6,'(' ITERATION =',I2,' CD =',F10.6)') I, CD C WRITE(6,1010) I, CD, USTARC, UAUB, ALPHA, USTACW, DELTAW, C * B, Z0C 100 Continue C 110 WRITE(6,*) 'CD ',CD 110 Continue Return 5000 Format (I2,8(1X,F7.4)) End Subroutine WFRIC(KBRAB,FW) C ****************************************** C SUBROUTINE WFRIC WAVE FRICTION C W. D. GRANT W.H.O.I. C GENERATED FEB 14, 1984 C LATEST REVISION AUG 21, 1986 C Determines wave friction factor, FW, from C relative roughness. C ****************************************** Real KBRAB C ------------------------------------------ C REGION 3: KBRAB .GT. 1.0 If (KBRAB.GT.1.0) Then ! NG If (A/ks<1) FW = 0.23 C ------------------------------------------ C REGION 2: KBRAB .LT. 0.08 ! NG If (A/ks>12.5) Else If (KBRAB.LT.0.08) Then FW = 0.13 * KBRAB ** 0.40 C ------------------------------------------ C REGION 1: 0.08 .LE. KBRAB .LE. 1.0 Else ! NG For intermediate A/ks FW = 0.23 * KBRAB ** 0.62 End If Return End Subroutine MOVKBR(DIA,UB,OMEGA,KBRAB) C ****************************************** C SUBROUTINE MOVKBR MOVEABLE BED ROUGHNESS C W. D. GRANT W.H.O.I. C GENERATED APR 1, 1984 C LATEST REVISION APR 1, 1984 C Calculates the relative roughness over C moveable beds approximately after Grant & C Madsen 1982. C All units are CGS. C This subroutine assumes quartz sand. C ****************************************** Real KBRAB, KRIPAB, KBREDAB C ------------------------------------------ C CONSTANTS S = 2.65 G = 980.6 VISC = 0.011 C ------------------------------------------ C COMMON FACTORS SUBWT = (S-1.0) * G * DIA SSTAR = (DIA/(4.0*VISC)) * SQRT(SUBWT) PHICR = 0.054 AB = UB / OMEGA KBRAB = DIA / AB RE = UB * AB / VISC TRANRO = 0.56 / SQRT(RE) If (KBRAB.LT.TRANRO) Then FW = 0.09 * (1.0/(RE**0.2)) Else Call WFRIC(KBRAB,FW) End If TAUWB = 1.0 * 0.5 * FW * (UB**2.0) PHIPRI = TAUWB / SUBWT RATIO = PHIPRI / PHICR BREAK = 1.8 * (SSTAR**0.6) If (RATIO.LE.BREAK) Then STEEP = 0.16 * (1.0/(RATIO**0.04)) HTAB = 0.22 * (1.0/(RATIO**0.16)) Else STEEP = 0.28 * ((SSTAR**0.6)*(1.0/RATIO)) HTAB = 0.48 * ((SSTAR**0.8)*(1.0/(RATIO**1.5))) End If KRIPAB = 30.0 * HTAB * STEEP KBREDAB = 27.2 * KBRAB * ((SQRT(RATIO)-0.7)**2.0) KBRAB = KRIPAB + KBREDAB Write (6,5000) KRIPAB, KBREDAB, KBRAB Return 5000 Format (/,5X,' KRIPAB =',F10.6,/,5X,' KBREDAB =',F10.6,/,5X, * ' KBRAB =',F10.6) End