WW3/w3sbt9md_8F90_source.html

#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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    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

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif


    !  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


    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SBT9')

#endif

    !

    ! 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