!> @file
!> @brief Contains routines for computing dissipation by viscous fluid mud using
!>  Ng (2000).
!>
!> @author M. Orzech
!> @author W. E. Rogers
!> @date   21-Nov-2013
!>

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!>
!> @brief Contains routines for computing dissipation by viscous fluid
!>  mud using Ng (2000).
!>
!> @author M. Orzech
!> @author W. E. Rogers
!> @date   21-Nov-2013
!>
!> @copyright Copyright 2009-2022 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.
!>
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
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Compute dissipation by viscous fluid mud using Ng (2000)
  !>  (adapted from Erick Rogers code by Mark Orzech, NRL).
  !>
  !> @param[in]  AC        Action density.
  !> @param[in]  H_WDEPTH  Mean water depth.
  !> @param[out] S         Source term (1-D version).
  !> @param[out] D         Diagonal term of derivative (1-D version).
  !> @param[in]  IX
  !> @param[in]  IY
  !>
  !> @author M. Orzech
  !> @author W. E. Rogers
  !> @date   21-Nov-2013
  !>
  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

  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Compute dissipation by viscous fluid mud using Ng (2000).
  !>
  !> @details Adapted from Erick Rogers code by Mark Orzech, NRL.
  !>
  !> @param[in]  SIGMA     Radian frequency (rad).
  !> @param[in]  H_WDEPTH  Water depth, denoted "h" in Ng (m).
  !> @param[in]  DTILDE    Normalized mud depth.
  !> @param[in]  ZETA      The ratio of stokes' boundary layer.
  !> @param[in]  SBLTM     Sbltm is what you get if you calculate sblt using
  !>                       the viscosity of the mud
  !> @param[in]  GAMMA     The gamma used in Ng page 238, density(water)/density(mud).
  !> @param[in]  WK        Unmuddy wavenumber.
  !> @param[out] WKDR      Muddy wavenumber.
  !> @param[out] DISS      Dissipation rate.
  !>
  !> @author E. Rogers
  !> @author M. Orzech
  !> @date   21-Nov-2013
  !>
  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

  !/ ------------------------------------------------------------------- /
  !>
  !> @brief NA
  !>
  !> @param[in]  KWAVE
  !> @param[in]  H_WDEPTH
  !> @param[in]  SND2
  !> @param[out] ND
  !>
  !> @author NA
  !> @date   NA
  !>
  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
