!> @file
!> @brief Contains routines for computing dissipation by viscous fluid mud using
!>  Dalrymple and Liu (1978) "Thin Model".
!>
!> @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
!>  Dalrymple and Liu (1978) "Thin Model".
!>
!> @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 W3SBT8MD
  !/
  !/                  +-----------------------------------+
  !/                  | 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
  !     Dalrymple and Liu (1978) "Thin Model".
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SBT8    Subr. Public   Fluid mud dissipation (Dalrymple & Liu, 1978)
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      CSINH     Subr.   ??     Complex sinh function
  !      CCOSH     Subr.   ??     Complex cosh function
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !     Historical information:
  !        This started as a matlab script provided to Erick Rogers by Tony
  !        Dalrympyle Sep 2006. Erick Rogers converted to Fortran and put
  !        it into the SWAN model May 2007. Mark Orzech adapted the code for
  !        WW3 and added it to NRL code repository July-Dec 2011.
  !        Erick Rogers brought it over to the NCEP repository May 2013
  !        and has been updating and maintaining it there.
  !
  !     Reference: Dalrymple, R.A., Liu,P.L.-F.,1978:
  !                Waves over soft muds :a 2-layer fluid model.
  !                Journal of Physical Oceanography, 8, 1121–1131.
  !
  !  6. Switches :
  !
  !     !/S  Enable subroutine tracing.
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  PUBLIC
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Compute dissipation by viscous fluid mud using Dalrymple and Liu (1978).
  !>
  !> @details "Thin Model" (adapted from Erick Rogers code by Mark Orzech, NRL).
  !>
  !> @param[in]  AC        Action density spectrum (1-D).
  !> @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 W3SBT8(AC,H_WDEPTH,S,D,IX,IY)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA      |
    !/                  |           M. Orzech     NRL       |
    !/                  |           W. E. Rogers  NRL       |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         21-Nov-2013 |
    !/                  +-----------------------------------+
    !/
    !/    20-Dec-2004 : Origination.                        ( version 3.06
    !/    23-Jun-2006 : Formatted for submitting code for   ( version 3.09 )
    !/                  inclusion in WAVEWATCH III.
    !/
    !  1. Purpose :
    !
    !     Compute dissipation by viscous fluid mud using Dalrymple and Liu (1978)
    !     "Thin Model" (adapted from Erick Rogers code by Mark Orzech, NRL).
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       AC        R.A.  I   Action density spectrum (1-D)
    !       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.
    !     ----------------------------------------------------------------
    !
    !  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 W3ODATMD, ONLY: NDSE
    USE W3SERVMD, ONLY: EXTCDE
    USE W3DISPMD, ONLY: WAVNU1
#ifdef W3_S
    USE W3SERVMD, ONLY: STRACE
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)    :: H_WDEPTH ! water depth
    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

    COMPLEX    :: K
    COMPLEX    :: SHH
    COMPLEX    :: CHH
    COMPLEX    :: SHD
    COMPLEX    :: CHD
    COMPLEX    :: LAM1
    COMPLEX    :: LAM2
    COMPLEX    :: CHLAM2
    COMPLEX    :: SHLAM2
    COMPLEX    :: A1
    COMPLEX    :: A2
    COMPLEX    :: A3
    COMPLEX    :: B1
    COMPLEX    :: B2
    COMPLEX    :: B3
    COMPLEX    :: B4
    COMPLEX    :: C0
    COMPLEX    :: TESTALF1
    COMPLEX    :: TESTALF2
    COMPLEX    :: ALF1
    COMPLEX    :: ALF2
    COMPLEX    :: PSI1
    COMPLEX    :: PSI2
    COMPLEX    :: M1
    COMPLEX    :: M0
    COMPLEX    :: C(4,3)
    COMPLEX    :: C41A
    COMPLEX    :: C42A
    COMPLEX    :: C43A
    COMPLEX    :: B(4)
    COMPLEX    :: CD
    COMPLEX    :: HH
    COMPLEX    :: DD
    COMPLEX    :: GG
    COMPLEX    :: FM1
    COMPLEX    :: KM1
    COMPLEX    :: FP
    COMPLEX    :: I
    COMPLEX    :: F
    COMPLEX    :: KMUD

    REAL    :: BET0
    REAL    :: KINVISW
    REAL    :: RHOW
    REAL    :: EXPH
    REAL    :: A
    REAL    :: K_UNMUD
    REAL    :: SIGMA       ! radian frequency (rad)
    REAL    :: SMUDWD(NK)  ! dissipation due to mud
    REAL    :: KMIMAG(NK)  ! imag part of kmud
    REAL    :: KD
    REAL    :: KDCUTOFF
    REAL    :: CWAVE
    REAL    :: ZTMP
    REAL    :: NWAVE_MUD
    REAL    :: CG_MUD
    REAL    :: KCHECK
    REAL    :: KTHRESHOLD
    REAL    :: RHOM
    REAL    :: KINVISM
    REAL    :: THICKM
    REAL    :: CG_UNMUD

    INTEGER :: ICOUNT
    INTEGER :: IK
    INTEGER :: IS

    PARAMETER (I=(0.,1.))

    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL STRACE (IENT, 'W3SBT8')
#endif
    !
    ! 0.  Initializations ------------------------------------------------ *
    !
    !     Dalrymple and Liu, Waves over soft muds:  1978.
    !     Thin layer solution.
    !     Matlab code provided by Tony Dalrymple
    !     Converted to Fortran by Erick Rogers

    ! 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
      THICKM = MUDT(IX,IY)
    ELSE
      WRITE(NDSE,*)'THICKM NOT SET'
      CALL EXTCDE ( 2 )
    ENDIF
    IF (INFLAGS1(0)) THEN
      KINVISM = MUDV(IX,IY)
    ELSE
      WRITE(NDSE,*)'KINVISM NOT SET'
      CALL EXTCDE ( 3 )
    ENDIF

    RHOW=DWAT        ! Density of seawater
    KINVISW=NU_WATER
    KDCUTOFF = 10.0
    KTHRESHOLD=1.0E-9

    A=1.0

    ! initialize matrix diagonal contributions
    D = 0.0
    S = 0.0

    IF ( THICKM>0.0 .AND. RHOM>0.0 .AND. KINVISM>0.0 ) THEN

      SMUDWD = 0.0

      ! *** loop over frequencies
      DO IK = 1,NK

        SIGMA = SIG(IK)

        !     un-muddy wave number, to start things off
        CALL WAVNU1(SIGMA,H_WDEPTH,K_UNMUD,CG_UNMUD)
        K=K_UNMUD

        !     start iterative loop

        DO ICOUNT=1,20  ! *** May need more ***

          CALL CSINH(K*H_WDEPTH,SHH)
          CALL CCOSH(K*H_WDEPTH,CHH)
          CALL CSINH(K*THICKM,SHD)
          CALL CCOSH(K*THICKM,CHD)

          !   define lambdas
          LAM1=SQRT(K*K-I*SIGMA/KINVISW)
          LAM2=SQRT(K*K-I*SIGMA/KINVISM)

          !   define hyperbolics on lamda2, lamda1
          CALL CCOSH(LAM2*THICKM,CHLAM2)
          CALL CSINH(LAM2*THICKM,SHLAM2)

          !   define exp decay
          EXPH=EXP(-LAM1*H_WDEPTH)

          !   define a1, a2, a3
          A1=-LAM2*SHD/K+SHLAM2
          A2=-CHD+CHLAM2
          A3=SHD-LAM2*SHLAM2/K

          !   define b1, b2, b3, b4
          B1=LAM1*SHH/K-CHH
          B2=LAM2*A2*SHH/K+A1*CHH
          B3=-A3*SHH+A2*CHH
          B4=LAM1*SHH/K+CHH

          ! define c0
          C0=B4*EXPH-(LAM1*LAM1+K*K)/(2*K*K)

          ! define beta0
          BET0=-EXPH/C0

          ! define alfa1, alfa2
          TESTALF1=-RHOM*KINVISM*(LAM2*LAM2+K*K)*(-LAM2/K)*CHD/K  &
               -2*RHOM*KINVISM*LAM2*CHLAM2-(RHOM-RHOW)*GRAV*(I*(A1)/SIGMA)
          TESTALF2=-RHOM*KINVISM*(LAM2*LAM2+K*K)*(-1)*SHD/K  &
               -2*RHOM*KINVISM*LAM2*SHLAM2-(RHOM-RHOW)*GRAV*(I*(A2)/SIGMA)
          ALF1=-TESTALF1
          ALF2=-TESTALF2

          ! define psi1, psi2
          PSI1=2*K*(-LAM2/K)*SHD+(LAM2*LAM2+K*K)*SHLAM2/K
          PSI2=2*K*(-1)*CHD+(LAM2*LAM2+K*K)*CHLAM2/K

          ! define M1, MO
          M1=I*RHOW*SIGMA/K-2*RHOW*KINVISW*K
          M0=B1+(LAM1*LAM1+K*K)*EXPH/(K*K)

          ! matrix coefficients (eq. 22)
          C(1,1)=LAM1*(BET0*M0+1)*SHH/K+(BET0*M0-1)*CHH+M0/C0+EXPH
          C(1,2)=(LAM1*BET0*B2+LAM2*A2)*SHH/K+(BET0*B2+A1)*CHH+B2/C0
          C(1,3)=(LAM1*BET0*B3/K-A3)*SHH+(BET0*B3+A2)*CHH+B3/C0

          ! matrix coefficients (eq. 23)
          C(2,1)=LAM1*(BET0*M0+1)*M1*CHH/K+(BET0*M0-1)*M1*SHH  &
               -2*RHOW*KINVISW*LAM1*M0/C0+2*RHOW*KINVISW*LAM1*EXPH
          C(2,2)=(LAM1*BET0*B2+LAM2*A2)*M1*CHH/K+(BET0*B2+A1)*M1*SHH &
               -2*RHOW*KINVISW*LAM1*B2/C0
          C(2,3)=(LAM1*BET0*B3/K-A3)*M1*CHH+(BET0*B3+A2)*M1*SHH  &
               -2*RHOW*KINVISW*LAM1*B3/C0

          ! matrix coefficients (eq. 21)
          C(3,1)=2*K*RHOW*KINVISW*(BET0*M0-1)+RHOW*KINVISW  &
               *(LAM1*LAM1+K*K)*(1-BET0*M0)/K
          C(3,2)=2*K*RHOW*KINVISW*(BET0*B2+A1)-RHOW*KINVISW  &
               *(LAM1*LAM1+K*K)*BET0*B2/K-I*RHOM*KINVISM*PSI1
          C(3,3)=2*K*RHOW*KINVISW*(BET0*B3+A2)-RHOW*KINVISW  &
               *(LAM1*LAM1+K*K)*BET0*B3/K-I*RHOM*KINVISM*PSI2

          ! matrix coefficients (eq.19)
          C41A=LAM1*M1*(BET0*M0+1)/K+2*RHOW*KINVISW*LAM1  &
               +2*RHOW*KINVISW*LAM1*BET0*M0
          C42A=M1*(LAM1*BET0*B2+LAM2*A2)/K  &
               +2*RHOW*KINVISW*LAM1*BET0*B2+ALF1
          C43A=M1*(LAM1*BET0*B3/K-A3)+2*RHOW*KINVISW*LAM1*BET0*B3+ALF2

          ! method 1
          C(4,1)=C41A*C(3,1)/C41A-C(3,1)
          C(4,2)=C42A*C(3,1)/C41A-C(3,2)
          C(4,3)=C43A*C(3,1)/C41A-C(3,3)

          ! force terms......righthand side
          B(1)=-I*SIGMA*A
          B(2)=RHOW*GRAV*A
          B(3)=0
          B(4)=0
          !  coefficients
          CD=-(C(3,3)-C(3,2)*C(4,3)/C(4,2))/C(3,1)
          HH=B(2)/(C(2,1)*CD -C(2,2)*(C(4,3)/C(4,2))+C(2,3))
          DD=CD*HH
          GG=-C(4,3)*HH/C(4,2)

          !  find k
          F=C(1,1)*DD+C(1,2)*GG+C(1,3)*HH-B(1)

          IF(ICOUNT.EQ.1)THEN
            FM1=F
            KM1=K
            K=K*(.995)+.001*I
            KCHECK=100.0
          ELSE
            KCHECK=ABS(IMAG(K)-IMAG(KM1))
            IF((F.EQ.FM1).OR.(K.EQ.KM1).OR.(KCHECK<KTHRESHOLD))THEN
              ! notes: I have noticed that if iterations are not stopped early enough, NaNs result
              KMUD=K
              EXIT
            END IF
            FP=(F-FM1)/(K-KM1)
            KM1=K
            FM1=F
            K=K-0.8*F/FP
          ENDIF
        END DO

        KMUD=K

        ! KD calc: not that important: just used to determine whether we make the mud calculation
        KD = REAL(KMUD) * H_WDEPTH

        IF ( KD .LT. KDCUTOFF ) THEN

          ! Notes: It would be better to have CG_MUD stored for each freq.
          CWAVE=SIGMA/REAL(KMUD)
          ZTMP=2.0*REAL(KMUD)*H_WDEPTH
          IF(ZTMP.LT.70)THEN
            ZTMP=SINH(ZTMP)
          ELSE
            ZTMP=1.0E+30
          ENDIF
          NWAVE_MUD=0.5*(1.0+2.0*REAL(KMUD)*H_WDEPTH/ZTMP)
          CG_MUD=NWAVE_MUD*CWAVE
          !        --- compute fluid mud-induced wave dissipation
          SMUDWD(IK)=2.0*IMAG(KMUD)*CG_MUD

          !        --- store imaginary part of KMUD
          KMIMAG(IK)=IMAG(KMUD)

        END IF !     IF ( KD .LT. KDCUTOFF ) THEN

      END DO ! IK

      !     *** store the results in the DIAGONAL array D ***
      DO IS = 1,NSPEC
        D(IS) = -SMUDWD(MAPWN(IS))
      END DO

    END IF !   IF ( THICKM>0.0 & RHOM>0.0 & KINVISM>0.0 ) THEN

    S = D * AC

    RETURN
    !/
    !/ End of W3SBT8 ----------------------------------------------------- /
    !/
  END SUBROUTINE W3SBT8

  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Complex hyperbolic sin (sinh).
  !>
  !> @param[in]  C
  !> @param[out] CS
  !>
  !> @author NA
  !> @date   NA
  !>
  SUBROUTINE CSINH(C,CS)
    COMPLEX, INTENT(IN) ::  C
    COMPLEX, INTENT(OUT) :: CS
    X = REAL(C)
    Y = AIMAG(C)
    CS = CMPLX(SINH(X) * COS(Y), SIN(Y) * COSH(X))
    RETURN
  END SUBROUTINE CSINH

  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Complex hyperbolic cos (cosh).
  !>
  !> @param[in]  C
  !> @param[out] CC
  !>
  !> @author NA
  !> @date   NA
  !>
  SUBROUTINE CCOSH(C,CC)
    COMPLEX, INTENT(IN) ::  C
    COMPLEX, INTENT(OUT) :: CC
    X = REAL(C)
    Y = AIMAG(C)
    CC = CMPLX(COSH(X) * COS(Y), SIN(Y) * SINH(X))
    RETURN
  END SUBROUTINE CCOSH

  !/ ------------------------------------------------------------------- /
  !/
END MODULE W3SBT8MD
