w3sdb1md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3sdb1md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           J. H. Alves             |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         29-May-2009 |
  !/                  +-----------------------------------+
  !/
  !/    25-Apr-2007 : Origination of module.              ( version 3.11 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/
  !/    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 :
  !
  !     Dummy slot for bottom friction source term.
  !
  !  2. Variables and types :
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SDB1    Subr. Public   Battjes and Janssen depth-induced
  !                               breaking.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !     See subroutine documentation.
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !     See subroutine documentation.
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  PUBLIC
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3sdb1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |                        FORTRAN 90 |
    !/                  |           J. H. Alves             |
    !/                  |           H. L. Tolman            |
    !/                  !           A. Roland               |
    !/                  | Last update :         08-Jun-2018 |
    !/                  +-----------------------------------+
    !/
    !/    25-Apr-2007 : Origination of module.              ( version 3.11 )
    !/    08-Jun-2018 : Add DEBUGDB1.                       ( version 6.04 )
    !/    03-Apr-2019 : Rewrite in terms of energy density (A. Roland,version 6.07)
    !/    03-Apr-2019 : Add Thornton & Guza, 1983          (A. Roland,version 6.07)
    !/
    !  1. Purpose :
    !
    !     Compute depth-induced breaking using Battjes and Janssen bore
    !     model approach
    !
    !  2. Method : Battjes & Janssen (1978),
    !
    !           Sbr   = Dtot/Etot*WA = D * WA
    !           Dtot  = 0.25*alpha*Qb*fm*Hmax²
    !           fm    = sigma/2Pi
    !           BB    = Hrms²/Hmax² = 8Etot/Hmax²
    !           D     = Dtot/Etot = BJALFA * sigma / pi * Qb/BB = 2 * BJALFA * fm * Qb/BB
    !
    !           AR: only valid for Hrms .le. Hm, Qb .le. 1, otherwise, in the degenrative regime it is
    !           due to Qb > 1 that all wave are broken and Hrms .le. Hmax
    !           MLIM can be used to enforce this conditions, source term will smoothly converge to this limit.
    !
    !     Where CDB   = SDBC1 = BJALFA (defaults to BJALFA = 1)
    !                   modified via ww3_grid namelist parameter BJALFA
    !           HM    = GAMMA * DEP
    !           GAMMA = SDBC2 defaults to 0.73 (mean Battjes/Janssen value)
    !                   modified via ww3_grid namelist parameter BJGAM
    !
    !     And QB is estimated by iterations using the nonlinear expression
    !
    !           1 - QB = HRMS**2
    !           ------   -------
    !            ln QB    HM**2
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.  I   Action density spectrum (1-D)
    !       EMEAN   Real  I   Mean wave energy.
    !       FMEAN   Real  I   Mean wave frequency.
    !       WNMEAN  Real  I   Mean wave number.
    !       DEPTH   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 :
    !
    !       STRACE   Subroutine tracing (!/S switch).
    !
    !  5. Called by :
    !
    !       W3SRCE   Source term integration.
    !       W3EXPO   Point output post-processor.
    !       GXEXPO   GrADS point output post-processor.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Note that the Miche criterion con influence wave growth.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/Tn  Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, sdbc1, sdbc2, fdonly, fssource, dden
    USE w3odatmd, ONLY: ndst
    USE w3gdatmd, ONLY: sig
    USE w3odatmd, only : iaproc
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T0
    USE w3arrymd, ONLY: prt2ds
#endif
#ifdef W3_T1
    USE w3arrymd, ONLY: outmat
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: IX ! Local grid number
    REAL, INTENT(IN)        :: A(NSPEC)
    REAL, INTENT(INOUT)     :: EMEAN, FMEAN, WNMEAN, DEPTH
    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)
    REAL, INTENT(IN)        :: CG(NK)
    LOGICAL, INTENT(OUT)    :: LBREAK
    INTEGER                 :: ITH, IK, IWB
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IS
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    real*8                    :: hm, bb, arg, q0, qb, b, cbj, hrms, eb(nk)
    real*8                    :: aux, cbj2, ratio, s0, s1, thr, br1, br2, fak
    REAL                      :: ETOT, FMEAN2
#ifdef W3_T0
    REAL                    :: DOUT(NK,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SDB1')
#endif
    !
    ! 0.  Initialzations ------------------------------------------------- /
    !     Never touch this 4 lines below ... otherwise my exceptionhandling will not work.
    s = 0.
    d = 0.
 
    thr = dble(1.e-15)
    IF (sum(a) .LT. thr) RETURN
 
    iwb = 1
    !
#ifdef W3_T
    WRITE (ndst,9000) sdbc1, sdbc2, fdonly
#endif
    !
    ! 1.  Integral quantities. AR: make sure mean quantities are computed, need to move upward
    !
    etot = 0.
    fmean2 = 0.
    DO ik=1, nk
      eb(ik) = 0.
      DO ith=1, nth
        eb(ik) = eb(ik) + a(ith+(ik-1)*nth)
      END DO
    END DO
    DO ik=1, nk
      eb(ik) = eb(ik) * dden(ik) / cg(ik)
      etot  = etot  + eb(ik)
    END DO
    DO ik=1, nk
      fmean2 = fmean2 + eb(ik) * sig(ik)
    END DO
    fmean2 = fmean2 / etot * tpiinv
    !
    ! 2do compute wlmean
    !
    ! 1.a. Maximum wave height
    ! 1.a.1. Simple limit
    !
    IF ( fdonly ) THEN
      hm     = dble(sdbc2) * dble(depth)
    ELSE
      !
      ! 1.a.2. Miche style criterion
      !
      hm     = dble(sdbc2) / dble(wnmean) * tanh( dble(wnmean) * max(depth,0.) )
    END IF
    !
    !AR: Add Dingemans ...
    ! 1.b. Hrms and ratio Hrms / Hmax
    !
    hrms = dsqrt(8.d0 * dble(emean))
    IF ( hm .GT. thr) THEN
      bb     = hrms * hrms / ( hm * hm )
      b      = dsqrt(bb)
    ELSE
      bb     = 0.d0
      b      = 0.d0
    END IF
    !
    ! 2. Fraction of breaking waves -------------------------------------- /
    ! 2.a. First guess breaking fraction
    !
    IF ( b .LE. 0.5d0 ) THEN
      q0     = 0.d0
    ELSE IF ( b .LE. 1.d0 ) THEN
      q0     = ( 2.d0 * b - 1.d0 ) ** 2
    END IF
    !
    ! 2.b. Iterate to obtain actual breaking fraction
    !
    IF ( b .LE. 0.2d0 ) THEN
      qb     = 0.d0
    ELSE IF ( b .LT. 1.d0 ) THEN
      arg    = exp(( q0 - 1.d0 ) / bb )
      qb     = q0 - bb * ( q0 - arg ) / ( bb - arg )
      DO is=1, 3
        qb     = exp((qb-1.)/bb)
      END DO
    ELSE
      qb = 1.0 - thr
    END IF
    !
    ! 3. Estimate the breaking coefficient ------------------------------- /
    !
    cbj  = 0
    IF (iwb == 1) THEN
      IF ( ( bb .GT. thr) .AND. ( abs( bb - qb ) .GT. thr) ) THEN
        IF ( bb .LT. 1.0) THEN
          cbj = 2 * dble(sdbc1) * qb * dble(fmean) / bb
        ELSE
          cbj = 2 * dble(sdbc1) * dble(fmean) * bb ! AR: degenerative regime, all waves must be .le. Hmax, we just smoothly let the excessive energy vanish by * BB.
        END IF
      ELSE
        cbj = 0.d0
      ENDIF
      d = - cbj
      s = d * a
    ELSE IF (iwb == 2) THEN
      IF (etot .GT. thr) THEN
        hrms = sqrt(8*emean)
        fak  = (1+4./sqrt(pi)*(b*bb+1.5*b)*exp(-bb)-erf(b))
        cbj  = -sdbc1*sqrt(pi)/16.*fmean*hrms**3/depth/etot
      ELSE
        cbj  = 0.
      ENDIF
      d = - cbj
      s = d * a
    ENDIF
 
    IF (cbj .GT. 0.) THEN
      lbreak = .true.
    ELSE
      lbreak = .false.
    ENDIF
 
#ifdef W3_DEBUGRUN
    IF (ix == debug_node) THEN
      WRITE(*,'(A200)') 'IX, DEPTH, CBJ, BB, QB, SDBC1, SDBC2, FMEAN, FMEAN2, HS'
      WRITE(*,'(I10,20F20.10)') ix, depth, cbj, bb, qb, sdbc1, sdbc2, fmean, fmean2, 4*sqrt(etot)
    ENDIF
#endif
    !
    ! ... Test output of arrays
    !
#ifdef W3_T0
    DO ik=1, nk
      DO ith=1, nth
        dout(ik,ith) = d(ith+(ik-1)*nth)
      END DO
    END DO
    CALL prt2ds (ndst, nk, nk, nth, dout, sig, '  ', 1.,    &
         0.0, 0.001, 'Diag Sdb', ' ', 'NONAME')
#endif
    !
#ifdef W3_T1
    CALL outmat (ndst, d, nth, nth, nk, 'diag Sdb')
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3SDB1 : PARAMETERS :',2f7.3,l4)
#endif
    !/
    !/ End of W3SDB1 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sdb1
  !/
  !/
  !/ End of module W3SDB1MD -------------------------------------------- /
  !/
END MODULE w3sdb1md