#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3SBT9MD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA | !/ | M. Orzech NRL | !/ | W. E. Rogers NRL | !/ | FORTRAN 90 | !/ | Last update : 21-Nov-2013 | !/ +-----------------------------------+ !/ !/ 28-Jul-2011 : Origination. ( version 4.01 ) !/ 21-Nov-2013 : Preparing distribution version. ( version 4.11 ) !/ !/ Copyright 2009 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! 1. Purpose : ! ! Contains routines for computing dissipation by viscous fluid mud using ! Ng (2000) ! ! 2. Variables and types : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3SBT9 Subr. Public Fluid mud dissipation (Ng 2000) ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! CSINH Subr. ?? Complex sinh function ! CCOSH Subr. ?? Complex cosh function ! Z_WNUMB Subr. ?? Compute wave number from freq & depth ! ---------------------------------------------------------------- ! ! 5. Remarks : ! Historical information: ! This started as some equations (the "B" parameter equations ! in subroutine "Ng" below) in a standalone Fortran ! code written by Jim Kaihatu, December 2004. These were adapted by ! Erick Rogers for a simple model based on governing equation ! similar to SWAN, and installed in a full version of SWAN in ! March 2005 with an informal report in May 2005. Kaihatu provided ! a "patch" for the B equations May 2006. Mud code in SWAN v40.41A was ! finalized June 2006, and v40.51 August 2007. The code was applied ! to Cassino Beach ~Sep 2006. This work was presented at a conference ! in Brazil Nov 2006, and later published in Rogers and Holland ! (CSR 2009). The code was adapted for WW3 by Mark Orzech in Nov 2012 ! (he had installed the D&L routines as BT8 in July 2011). ! ! Reference: Ng, C.O.,2000. Water waves over a muddy bed: ! a two-layer Stokes’ boundary layer model. ! Coastal Engineering 40(3),221–242. ! ! 6. Switches : ! ! !/S Enable subroutine tracing. ! ! 7. Source code : !/ !/ ------------------------------------------------------------------- / !/ ! PUBLIC !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE W3SBT9(AC,H_WDEPTH,S,D,IX,IY) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA | !/ | M. Orzech NRL | !/ | W. E. Rogers NRL | !/ | FORTRAN 90 | !/ | Last update : 21-Nov-2013 | !/ +-----------------------------------+ !/ 28-Jul-2011 : Origination. ( version 4.01 ) !/ 21-Nov-2013 : Preparing distribution version. ( version 4.11 ) !/ ! 1. Purpose : ! ! Compute dissipation by viscous fluid mud using Ng (2000) ! (adapted from Erick Rogers code by Mark Orzech, NRL). ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! H_WDEPTH Real I Mean water depth. ! S R.A. O Source term (1-D version). ! D R.A. O Diagonal term of derivative (1-D version). ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! CALC_ND ! NG ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SRCE Subr. W3SRCEMD Source term integration. ! W3EXPO Subr. N/A Point output post-processor. ! GXEXPO Subr. N/A GrADS point output post-processor. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! Cg_mud calculation could be improved by using dsigma/dk instead ! of n*C. The latter is a "naive" method and its accuracy has ! not been confirmed. ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD, ONLY: NK,SIG,NSPEC,MAPWN USE W3IDATMD, ONLY: MUDT, MUDV, MUDD, INFLAGS1 USE CONSTANTS, ONLY: PI,GRAV,DWAT,NU_WATER USE W3SERVMD, ONLY: EXTCDE USE W3ODATMD, ONLY: NDSE !/S USE W3SERVMD, ONLY: STRACE !/ IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL, INTENT(IN) :: H_WDEPTH ! WATER DEPTH, DENOTED "H" IN NG (M) REAL, INTENT(IN) :: AC(NSPEC) ! ACTION DENSITY INTEGER, INTENT(IN) :: IX, IY REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ !/S INTEGER, SAVE :: IENT = 0 ! LOCAL VARIABLES REAL :: DMW(NK) REAL :: ROOTDG REAL :: SND REAL :: SND2 REAL :: WGD REAL :: CWAVE REAL :: KD_ROCK REAL :: CG_MUD REAL :: K_MUD REAL :: NWAVE_MUD REAL :: ND_MUD REAL :: SMUDWD(NK) ! DISSIPATION DUE TO MUD REAL :: CG_ROCK REAL :: K_ROCK REAL :: NWAVE_ROCK REAL :: ND_ROCK REAL :: KINVISM ! := THE KINEMATIC VISCOSITY OF THE MUD REAL :: KINVISW ! := KINEMATIC VISCOSITY OF WATER REAL :: RHOW ! := DENSITY OF WATER REAL :: RHOM ! := DENSITY OF MUD REAL :: DM ! := DEPTH OF MUD LAYER REAL :: ZETA ! := THIS IS ZETA AS USED IN NG PG. 238. IT IS THE ! RATIO OF STOKES' BOUNDARY LAYER THICKNESSES, ! OR DELTA_M/DELTA_W REAL :: GAMMA ! := THIS IS THE GAMMA USED IN NG PG. 238. THIS IS ! DENSITY(WATER)/DENSITY(MUD) REAL :: SBLTW ! := A FUNCTION OF VISCOSITY AND FREQ REAL :: SBLTM ! := A FUNCTION OF VISCOSITY AND FREQ REAL :: DTILDE ! := NORMALIZED MUD DEPTH = MUD DEPTH / DELTA_M, ! DELTA IS THE SBLT= SQRT(2*VISC/SIGMA) REAL :: ZTMP REAL :: KDCUTOFF REAL :: KD INTEGER :: IS LOGICAL :: INAN !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3SBT9') ! ! 0. Initializations ------------------------------------------------ * ! ! Ng (2000), Waves over soft muds. ! Based on code for SWAN created by Erick Rogers. ! Adapted for WW3 by Mark Orzech, Nov 2012. ! Initialize properties from mud fields if available IF (INFLAGS1(-2))THEN RHOM = MUDD(IX,IY) ELSE WRITE(NDSE,*)'RHOM NOT SET' CALL EXTCDE ( 1 ) ENDIF IF (INFLAGS1(-1)) THEN DM = MUDT(IX,IY) ELSE WRITE(NDSE,*)'DM NOT SET' CALL EXTCDE ( 2 ) ENDIF IF (INFLAGS1(0)) THEN KINVISM = MUDV(IX,IY) ELSE WRITE(NDSE,*)'KINVISM NOT SET' CALL EXTCDE ( 3 ) ENDIF ROOTDG = SQRT(H_WDEPTH/GRAV) WGD = ROOTDG*GRAV DO IS = 1, NK ! SND is dimensionless frequency SND = SIG(IS) * ROOTDG IF (SND .GE. 2.5) THEN ! ******* DEEP WATER ******* K_ROCK = SIG(IS) * SIG(IS) / GRAV CG_ROCK = 0.5 * GRAV / SIG(IS) NWAVE_ROCK = 0.5 ND_ROCK = 0. ELSE IF (SND.LT.1.E-6) THEN ! *** VERY SHALLOW WATER *** K_ROCK = SND/H_WDEPTH CG_ROCK = WGD NWAVE_ROCK = 1. ND_ROCK = 0. ELSE SND2 = SND*SND CWAVE = SQRT(GRAV*H_WDEPTH/(SND2+1./(1.+0.666*SND2 & +0.445*SND2**2 -0.105*SND2**3+0.272*SND2**4))) K_ROCK = SIG(IS)/CWAVE CALL CALC_ND(K_ROCK,H_WDEPTH,SND2,ND_ROCK) NWAVE_ROCK = 0.5*(1.0+2.0*K_ROCK*H_WDEPTH/SINH(2.0*K_ROCK*H_WDEPTH)) CG_ROCK= NWAVE_ROCK*CWAVE SND2=0 CWAVE=0 ENDIF KDCUTOFF = 10.0 ! hardwired (same as w3sbt8md) ! now that kh is known, we can use a definition of "deep" that is ! consistent with the definition used in sbot K_MUD=0.0 DMW(IS)=0.0 KD_ROCK = K_ROCK * H_WDEPTH ! KD_ROCK is used to determine whether we make the mud calculation IF((KD_ROCK.LT.KDCUTOFF).AND.(DM.GT.1.0E-5))THEN KINVISW=NU_WATER RHOW=DWAT ZETA=SQRT(KINVISM/KINVISW) GAMMA=RHOW/RHOM SBLTW=SQRT(2.0*KINVISW/SIG(IS)) SBLTM=SQRT(2.0*KINVISM/SIG(IS)) DTILDE=DM/SBLTM CALL NG(SIG(IS),H_WDEPTH,DTILDE,ZETA,SBLTM,GAMMA,K_ROCK,K_MUD, & DMW(IS)) ELSE ! IF ( KD_ROCK .LT. KDCUTOFF ) THEN K_MUD=K_ROCK END IF ! IF ( KD_ROCK .LT. KDCUTOFF ) THEN ! calculate cg_mud, nwave_mud here CWAVE=SIG(IS)/K_MUD ZTMP=2.0*K_MUD*H_WDEPTH IF(ZTMP.LT.70)THEN ZTMP=SINH(ZTMP) ELSE ZTMP=1.0E+30 ENDIF NWAVE_MUD=0.5*(1.0+2.0*K_MUD*H_WDEPTH/ZTMP) CG_MUD=NWAVE_MUD*CWAVE SND2 = SND*SND CALL CALC_ND(K_MUD,H_WDEPTH,SND2,ND_MUD) SND2=0 CWAVE=0 ! If we wanted to include the effects of mud on the real part of the ! wavnumber (as we do in SWAN), this is where we would do it. ! Set output variables k_out, cg_out, nwave_out, nd_out, dmw. !kinematics IF(MUD)THEN ! !kinematics K_OUT(IS) =K_MUD !kinematics CG_OUT(IS) =CG_MUD !kinematics NWAVE_OUT(IS)=NWAVE_MUD !kinematics ND_OUT(IS) =ND_MUD !kinematics ELSE ! USE ROCKY WAVENUMBER,ETC. !kinematics K_OUT(IS) =K_ROCK !kinematics CG_OUT(IS) =CG_ROCK !kinematics NWAVE_OUT(IS)=NWAVE_ROCK !kinematics ND_OUT(IS) =ND_ROCK !kinematics DMW(IS)=0.0 !kinematics ENDIF KD = K_MUD * H_WDEPTH IF ( KD .LT. KDCUTOFF ) THEN ! note that "IS" here is for the 1d spectrum SMUDWD(IS)=2.0*DMW(IS)*CG_MUD END IF ! NaN check: INAN = .NOT. ( DMW(IS) .GE. -HUGE(DMW(IS)) .AND. DMW(IS) & .LE. HUGE(DMW(IS)) ) IF (INAN) THEN WRITE(*,'(/1A/)') 'W3SBT9 ERROR -- DMW(IS) IS NAN' WRITE(*,*)'W3SBT9: RHOM, DM, KINVISM = ',RHOM, DM, KINVISM WRITE(*,*)'W3SBT9: IS,NK = ',IS,NK WRITE(*,*)'W3SBT9: H_WDEPTH,KD,KDCUTOFF = ',H_WDEPTH,KD, KDCUTOFF WRITE(*,*)'W3SBT9: K_MUD,CG_MUD,NWAVE_MUD = ',K_MUD,CG_MUD,NWAVE_MUD CALL EXTCDE (1) END IF END DO ! DO IS = 1, NK ! *** store the results in the DIAGONAL arrays D and S *** DO IS = 1,NSPEC ! note that "IS" here is for the directional spectrum (2d) D(IS) = -SMUDWD(MAPWN(IS)) END DO S = D * AC RETURN END SUBROUTINE W3SBT9 !/ ------------------------------------------------------------------- / SUBROUTINE NG(SIGMA,H_WDEPTH,DTILDE,ZETA,SBLTM,GAMMA,WK,WKDR,DISS) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | E. Rogers and M. Orzech | !/ | FORTRAN 90 | !/ | Last update : 21-Nov-2013 | !/ +-----------------------------------+ ! !/ 28-Jul-2011 : Origination. ( version 4.01 ) !/ 21-Nov-2013 : Preparing distribution version. ( version 4.11 ) !/ ! 1. Purpose : ! ! Compute dissipation by viscous fluid mud using Ng (2000) ! (adapted from Erick Rogers code by Mark Orzech, NRL). ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! SIGMA Real I radian frequency (rad) ! H_WDEPTH Real I water depth ! DTILDE Real I normalized mud depth ! ZETA Real I zeta as used in Ng ! SBLTM Real I mud Stokes boundary layer thickness ! GAMMA Real I gamma as used in Ng ! WK Real I wavenumber w/out mud ! WKDR Real O wavenumber w/mud ! DISS Real O dissipation rate ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! None. ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SBT9 Subr. W3SBT9MD Main routine (all freqs) ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! Calculations for the "B coefficients" came from a code by Jim Kaihatu ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! None. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ ------------------------------------------------------------------- / !/ ! IMPLICIT NONE ! INPUT VARIABLES : REAL, INTENT(IN) :: SIGMA ! radian frequency (rad) REAL, INTENT(IN) :: H_WDEPTH! water depth, denoted "h" in Ng (m) REAL, INTENT(IN) :: DTILDE ! normalized mud depth = mud depth / sbltm, ! delta is the sblt= sqrt(2*visc/sigma) REAL, INTENT(IN) :: ZETA ! this is zeta as used in Ng pg. 238. it is ! the ratio of stokes' boundary layer ! thicknesses, or sbltm/delta_w REAL, INTENT(IN) :: GAMMA ! this is the gamma used in Ng pg. 238. ! this is density(water)/density(mud) REAL, INTENT(IN) :: SBLTM ! sbltm is what you get if you calculate ! sblt using the viscosity of the mud, ! sbltm=sqrt(2*visc_m/sigma) ! .....also delta_m REAL, INTENT(IN) :: WK ! unmuddy wavenumber ! OUTPUT VARIABLES : REAL, INTENT(OUT) :: WKDR ! muddy wavenumber REAL, INTENT(OUT) :: DISS ! dissipation rate ! LOCAL VARIABLES : REAL :: B1 ! an Ng coefficient REAL :: B2 ! an Ng coefficient REAL :: B3 ! an Ng coefficient REAL :: BR ! an Ng coefficient REAL :: BI ! an Ng coefficient REAL :: BRP ! an Ng coefficient REAL :: BIP ! an Ng coefficient REAL :: DM ! MUD DEPTH, ADDED JUNE 2 2006 DM=DTILDE*SBLTM ! DTILDE=DM/SBLTM ! NOW CALCULATE Ng's B coefficients : see Ng pg 238 B1=GAMMA*(-2.0*GAMMA**2+2.0*GAMMA-1.-ZETA**2)*SINH(DTILDE)* & COSH(DTILDE)-GAMMA**2*ZETA*((COSH(DTILDE))**2+ & (SINH(DTILDE))**2)-(GAMMA-1.)**2*ZETA*((COSH(DTILDE))**2 & *(COS(DTILDE))**2+(SINH(DTILDE))**2*(SIN(DTILDE))**2)-2.0 & *GAMMA*(1.-GAMMA)*(ZETA*COSH(DTILDE)+GAMMA*SINH(DTILDE)) & *COS(DTILDE) B2=GAMMA*(-2.0*GAMMA**2+2.0*GAMMA-1.+ZETA**2)*SIN(DTILDE)* & COS(DTILDE) -2.0*GAMMA*(1.-GAMMA)*(ZETA*SINH(DTILDE)+GAMMA & *COSH(DTILDE))*SIN(DTILDE) B3=(ZETA*COSH(DTILDE)+GAMMA*SINH(DTILDE))**2*(COS(DTILDE))**2 & +(ZETA*SINH(DTILDE)+GAMMA*COSH(DTILDE))**2*(SIN(DTILDE))**2 BR=WK*SBLTM*(B1-B2)/(2.0*B3)+GAMMA*WK*DM BI=WK*SBLTM*(B1+B2)/(2.0*B3) BRP=B1/B3 ! "B_R PRIME" BIP=B2/B3 ! "B_I PRIME" ! now calculate dissipation rate and wavenumber DISS=-SBLTM*(BRP+BIP)*WK**2/(SINH(2.0*WK*H_WDEPTH)+2.0*WK*H_WDEPTH) WKDR=WK-BR*WK/(SINH(WK*H_WDEPTH)*COSH(WK*H_WDEPTH)+WK*H_WDEPTH) RETURN END SUBROUTINE NG !/ ------------------------------------------------------------------- / SUBROUTINE CALC_ND(KWAVE,H_WDEPTH,SND2,ND) !/ ------------------------------------------------------------------- / IMPLICIT NONE REAL, INTENT(IN) :: KWAVE REAL, INTENT(IN) :: H_WDEPTH REAL, INTENT(IN) :: SND2 REAL, INTENT(OUT) :: ND REAL :: FAC1 ! LOCAL REAL :: FAC2 ! LOCAL REAL :: FAC3 ! LOCAL REAL :: KND ! LOCAL KND = KWAVE*H_WDEPTH FAC1 = 2.*KND/SINH(2.*KND) FAC2 = SND2/KND FAC3 = 2.*FAC2/(1.+FAC2*FAC2) ND= FAC1*(0.5/H_WDEPTH - KWAVE/FAC3) END SUBROUTINE CALC_ND !/ ------------------------------------------------------------------- / !/ END MODULE W3SBT9MD