w3sbt8md.F90

Go to the documentation of this file.

 
 
#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
#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
 
  !/ ------------------------------------------------------------------- /
  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