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 BED_STRESS C CALCULATES WAVE-ENHANCED BOTTOM FRICTION COEFFICIENT C Reza Marsooli, 2015 INCLUDE 'comdeck' c SAVE REAL*8 COEF_KS,DELTAWC,DELTAWCmax,ORBVEL REAL*8 VEL_CX_TEMP,VEL_CY_TEMP,VEL_CX,VEL_CY REAL*8 UNIT_CX,UNIT_CY REAL*8 DIR_W,UNIT_WX,UNIT_WY REAL*8 ANGLE_C,ANGLE_W,ANGLE_WC,COSANGLE_WC REAL*8 UFACE,VFACE,VEL_CMAG REAL*8 COEF_AB,RATIO_ABKS,FRICTION_FW,USTAR_W REAL*8 TAU_C REAL*8 FRICTION_FCINI,FRICTION_FC0,FRICTION_FC REAL*8 USTAR_C,USTAR_WC,COEF_BETA,COEF_KBC,DIFF_FRICTION DO J=2,JM-1 DO I=2,IM-1 C BED PHYSICAL ROUGHNESS COEF_KS=Z0B(I,J) C REFERENCE HEIGHT, WAVE-CURRENT BOTTOM BOUNDARY LAYER THICKNESS DELTAWCmax=0.2 !m cc DELTAWC=AMIN1( DELTAWCmax , AMAX1(2.*COEF_KS , 0.01*DT(I,J)) ) DELTAWC=AMIN1(DELTAWCmax , AMAX1(2.*COEF_KS , 0.05*DT(I,J))/30.) C CALCULATE ONLY FOR NON-BLOCKED CELLS AND WET CELLS IF(FSM(I,J).GT.0.0.AND.DHMWAD(I,J).GT.0.0) THEN C WAVE ORBITAL VELOCITY OMEGA=2.0*3.141593/AMAX1(WVPD(I,J),0.05) ORBVEL=OMEGA*WVHT(I,J)/(2.*SINH(WN(I,J)*DT(I,J))) + 0.000001 C CALCULATE MAGNITUDE OF NEAR-BOTTOM CURRENT IF(TOR.EQ.'BAROTROPIC') THEN UFACE=0.5*(UA(I,J)+UA(I+1,J)) VFACE=0.5*(VA(I,J)+VA(I,J+1)) ELSE UFACE=0.5*(U(I,J,KBM1)+U(I+1,J,KBM1)) VFACE=0.5*(V(I,J,KBM1)+V(I,J+1,KBM1)) ENDIF VEL_CMAG=SMALLVAL+SQRT(UFACE**2+VFACE**2) IF (ORBVEL.LT.0.01.OR.VEL_CMAG.LT.0.01) GOTO 10 C FIND THE ANGLE (IN RADIANTS) BETWEEN WAVES AND CURRENTS VEL_CX_TEMP = 0.5*(UA(I,J)+UA(I+1,J)) VEL_CY_TEMP = 0.5*(VA(I,J)+VA(I,J+1)) VEL_CX = VEL_CX_TEMP*COSANG(I,J)-VEL_CY_TEMP*SINANG(I,J) !w.r.t west-east axis VEL_CY = VEL_CX_TEMP*SINANG(I,J)+VEL_CY_TEMP*COSANG(I,J) !w.r.t south-north axis IF(VEL_CX.eq.0.0) THEN ANGLE_C=0.5*3.141593*SIGN(1.0,VEL_CY) ELSE ANGLE_C=ATAN2(VEL_CY,VEL_CX) ENDIF ANGLE_W=thtav(I,J) ANGLE_WC=ANGLE_W-ANGLE_C COSANGLE_WC=COS(ANGLE_WC) c CONVERT THE LENGTH SCALE FROM M TO CM COEF_KS=COEF_KS*100. DELTAWC=DELTAWC*100. ORBVEL=ORBVEL*100. VEL_CMAG=VEL_CMAG*100. C-----------Calculate Wave-Enhanced Bottom Friction Coef. & Shear Stress--- C EXCURSION AMPLITUDE COEF_AB=ABS(ORBVEL/(OMEGA+SMALLVAL)) C RATIO BETWEEN EXCURSION AMPLITUDE AND BED ROUGHNESS RATIO_ABKS=COEF_AB/(COEF_KS+SMALLVAL) C RATIO BETWEEN BED ROUGHNESS AND EXCURSION AMPLITUDE RATIO_KSAB=COEF_KS/(COEF_AB+SMALLVAL) C WAVE FRICTION FACTOR - METHOD OF JONSSON 1966 CC IF(RATIO_ABKS.LE.1.57) THEN CC FRICTION_FW=0.3 CC ELSE CC FRICTION_FW=EXP(-6.+5.2/RATIO_ABKS**0.19) CC ENDIF C WAVE FRICTION FACTOR - METHOD OF SIGNELL ET AL 1990 IF(RATIO_KSAB.LT.0.08) THEN FRICTION_FW=0.13*RATIO_KSAB**0.4 ELSEIF(RATIO_KSAB.LT.1.0) THEN FRICTION_FW=0.23*RATIO_KSAB**0.62 ELSE FRICTION_FW=0.23 ENDIF C CALCULATE WAVE-ALONE SHEAR VELOCITY USTAR_W=SQRT(0.5*FRICTION_FW)*ORBVEL ! CM/S C C CALCULATE INITIAL WAVE-CURRENT FRICTION FACTOR (NO WAVE EFFECTS) FRICTION_FCINI=2.*(0.41/ALOG(30.*DELTAWC/COEF_KS))**2 C LIMIT FRICTION FACTOR FRICTION_FCINI=AMIN1(1.0,FRICTION_FCINI) C INITIALIZE THE ITERATION PROCEDURE FRICTION_FC0=FRICTION_FCINI FRICTION_FC =FRICTION_FCINI ITERCOUNT=0 DIFF_FRICTION=999.99 C START THE ITERATION PROCESS TO CALCULATE COMBINED WAVE-CURRENT FRICTION DO WHILE ((ITERCOUNT.LT.20).AND.(DIFF_FRICTION.GT.0.0002)) ITERCOUNT=ITERCOUNT+1 C CALCULATE CURRENT-ALONE SHEAR VELOCITY USTAR_C=SQRT(0.5*FRICTION_FC)*VEL_CMAG ! CM/S C CALCULATE COMBINED WAVE-CURRENT SHEAR VELOCITY USTAR_WC=SQRT( USTAR_C**2+USTAR_W**2 . +2.*USTAR_C*USTAR_W*COSANGLE_WC ) ! CM/S C CALCULATE APPARNET BOTTOM ROUGHNESS COEF_BETA=1.-USTAR_C/(USTAR_WC+SMALLVAL) COEF_KBC=COEF_KS*(24.*(USTAR_WC/ORBVEL) . *RATIO_ABKS)**COEF_BETA C UPDATE THE REFERENCE THICKNESS cc DELTAWC=1.5*0.41*USTAR_WC/OMEGA cc DELTAWC=AMIN1(DELTAWCmax*100. , AMAX1(2.*COEF_KBC , DELTAWC)) C UPDATE WAVE-CURRENT FRICTION FACTOR (WITH WAVE EFFECTS INCLUDED IN COEF_KBC) FRICTION_FC=2.*(0.41/ALOG(30.*DELTAWC . /(COEF_KBC+SMALLVAL)))**2 DIFF_FRICTION=ABS(FRICTION_FC-FRICTION_FC0) FRICTION_FC0=FRICTION_FC ENDDO CBC(I,J)=AMIN1(AMAX1(FRICTION_FC,CBC(I,J)),1.0) IF(ITERCOUNT.GT.19) WRITE(4155,'(F15.8,2I5,3F15.8)') . THOUR,I,J,COEF_KBC,DELTAWC,FRICTION_FC if(i.eq.53.and.j.eq.57) write(133,'(9f15.8)') . THOUR,COEF_KBC,DELTAWC, . ORBVEL,FRICTION_FW,FRICTION_FC,USTAR_W,USTAR_C,USTAR_WC ENDIF 10 CONTINUE ENDDO ENDDO RETURN END !SUBROUTINE BED_STRESS