#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      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
!/ ------------------------------------------------------------------- /
      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
!/S      USE W3SERVMD, ONLY: STRACE
!/
      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
!/
!/S      INTEGER, SAVE           :: IENT = 0
     
      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.))

!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3SBT8')
!
! 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

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

      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

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

      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