w3snl5md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3snl5md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           O. Gramstad             |
  !/                  |           Q. Liu                  |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         07-Jun-2021 |
  !/                  +-----------------------------------+
  !/
  !/    24-Sep-2013 : Origination.                        ( version 3.14 )
  !/    24-Sep-2013 : GKE for the regular wavenumbergrid  ( O. Gramstad  )
  !/                  (interpolation required)
  !/    02-Dec-2013 : GKE for WW3 logarithmic freq. grid  ( O. Gramstad  )
  !/                  (single grid point)
  !/    27-Feb-2019 : GKE for 2D applications.            ( version 7.13 )
  !/                                                      ( Q. Liu )
  !/    07-06-2021  : Merge into WW3 Github               ( version 7.13 )
  !/                                                      ( Q. Liu       )
  !/
  !  1. Purpose :
  !     Interface module for GKE (resonant & quasi-resonant four-wave
  !     interactions)
  !
  !  2. Variables and types :
  !
  !  3. Subroutines and functions :
  !
  !      Name         Type  Scope    Description
  !     -------------------------------------------------------------------
  !      W3SNL5       Subr. Public   Interface to gkeModule
  !      INSNL5       Subr. Public   Initialization routine
  !      CALC_WBTv2   Subr. Private  Calc. dominant wave breaking prob.
  !      INPOUT       Subr. Private  Point output
  !     -------------------------------------------------------------------
  !
  !  4. Future work: Dnl
  !/
  !/ ------------------------------------------------------------------- /
  IMPLICIT NONE
  !/
  ! Subrs.
  PUBLIC      :: w3snl5, insnl5
  PRIVATE     :: calc_wbtv2, inpout
  ! Vars.
  PRIVATE     :: nsel, psea, pnms
  !/ ------------------------------------------------------------------- /
  ! Parameter list
  INTEGER                              :: NSEL
  INTEGER, ALLOCATABLE, SAVE           :: PSEA(:)
  CHARACTER(LEN=10), ALLOCATABLE, SAVE :: PNMS(:)
  !
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3snl5(A, CG, WN, FMEAN, T1ABS, U10, UDIR, JSEA, &
       S, D, KURT)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           O. Gramstad             |
    !/                  |           Q. Liu                  |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         24-Apr-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Sep-2013 : Origination.                        ( version 3.14 )
    !/    24-Sep-2013 : GKE for resonant & quasi-resonant four-wave
    !/                  interactions                        ( O. Gramstad  )
    !/    27-Feb-2019 : Full implementation of GKE          ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/    21-Apr-2019 : Phase mixing option                 ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/    24-Apr-2019 : Phase mixing option (b_T)           ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/    02-May-2019 : Organize screen output & disable binary output
    !/                                                      ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/
    !/
    !  1. Purpose :
    !
    !     Interface to CalcQRSNL subr. of the GKE module. Please refer to
    !     -------------
    !     gkeModule.f90 for further details.
    !     -------------
    !
    !     ◆ Different times used in this module
    !
    !     |----o---------o----o--|-|--o-----o------o-----o------o----> (t)
    !     ^    ^         ^       ^ ^ T1ABS (absol. current time step)¹
    !     |    |         |       |
    !     |    |         |       v t0 (relat. time, previous time step)
    !     |    |<------->|
    !     |    | PM_IVAL (phase mixing interval, relat. time)
    !     |    |
    !     |    v PM_PREV (phase mixing, appear quasi-periodically)
    !     |              (relat. time)
    !     |
    !     v TBEG (absol. begining time, defined by ww3_shel.inp)
    !
    !     ¹ Because of using the dynamic integration scheme, T1ABS
    !       is related to, but not the same as, TIME in w3wdatmd.ftn
    !  2. Method :
    !
    !  3. Parameters :
    !
    !  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 :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav, tpi
    USE w3gkemd,   ONLY: calcqrsnl, qr_depth
    USE w3gdatmd,  ONLY: nk, nth, nspec, sig, th,            &
         gtype, rlgtype, clgtype,            &
         qr5dpt, qi5nnz, qi5pmx
    USE w3wdatmd,  ONLY: qi5tbeg, qr5tim0, qr5cvk0, qc5int0, &
         qr5tmix
    USE w3odatmd,  ONLY: flout, nopts, tosnl5, tolast, &
         iaproc, napout, screen
    USE w3parall,  ONLY: init_get_isea
    USE w3timemd,  ONLY: dsec21
    USE w3servmd,  ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: a(nth, nk)         ! N(θ, k)
    REAL, INTENT(IN)        :: cg(nk)             ! Cg(k)
    REAL, INTENT(IN)        :: wn(nk)             ! WN(k)
    REAL, INTENT(IN)        :: fmean              ! 1/T_{0, -1}
    INTEGER, INTENT(IN)     :: t1abs(2)           ! Absol. t1
    REAL, INTENT(IN)        :: u10                ! Wind velocity
    REAL, INTENT(IN)        :: udir               ! φ (in rad)
    INTEGER, INTENT(IN)     :: jsea               ! Local sea point count
    REAL, INTENT(OUT)       :: s(nth,nk),      &  ! Snl
         d(nth,nk),      &  ! Dnl
         kurt               ! Kurtosis
 
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    REAL, PARAMETER         :: btlow = 10., bthgh = 500.
    REAL                    :: t0rel, t1rel, tdel1, tdel2
    REAL                    :: cvk1(nspec), snl(nspec), dnl(nspec)
    REAL                    :: cvk0(nspec)
    COMPLEX                 :: inpqr0(qi5nnz)
    INTEGER                 :: ik, ith, ispec, isea, jloc
    INTEGER, ALLOCATABLE    :: pdiff(:)
    LOGICAL, SAVE           :: fstout = .true.
    REAL                    :: factor(nk), a2(nk, nth), s2(nk, nth)
    REAL                    :: pm_prev, pm_ival, pm_delt
    REAL                    :: wbt, btinv
    INTEGER                 :: iunt
#ifdef W3_S
    INTEGER, SAVE           :: ient = 0
#endif
 
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SNL5')
#endif
    !
    !/ ------------------------------------------------------------------- /
    ! Read in wave info. @ the previous time step t0
    ! Array initialization is done in w3wdat/w3setw (called by w3initmd)
    t0rel  = qr5tim0(jsea)              ! t0 (nsea)
    cvk0   = qr5cvk0(:, jsea)           ! Cvk (ns, nsea) @ t0
    inpqr0 = qc5int0(:, jsea)           ! Inpqr (nnz, nsea) @ t0
    !
    ! Calc. Relative time for T1ABS (QI5TBEG as the reference)
    t1rel  = dsec21(qi5tbeg, t1abs)     ! in second
    !
    ! W3WAVEMD: IF ( IT.EQ.0 ) DTG = 1 → T1REL = -1 (the first step of
    ! integration; TIME = TCALC/TOFRST, DTG = 0 →  1, QI5TBEG = TIME - 1)
    IF(t1rel < 0.) t1rel = 0.
    !
    ! Three options for phase mixing
    IF (qi5pmx .EQ. 0) THEN
      ! 1) 0: no phase mixing
#ifdef W3_TS
      IF (iaproc .EQ. napout)                                    &
           WRITE(screen, '(A, 2(I10.8, I7.6), E12.3)')                &
           " ⊚ → [WW3 SNL₅] QI5TBEG, T1ABS, T1REL:",              &
           qi5tbeg, t1abs, t1rel
#endif
    ELSE
#ifdef W3_TS
      IF (iaproc .EQ. napout)                                    &
           WRITE(screen, '(A, 2(I10.8, I7.6), E12.3)', advance='no')  &
           " ⊚ → [WW3 SNL₅] QI5TBEG, T1ABS, T1REL, T1REL[P]:",    &
           qi5tbeg, t1abs, t1rel
#endif
      !
      ! Calc. Phase mixing interval
      IF (qi5pmx .GT. 0) THEN
        ! 2) N: mix phase by every N characteristic wave periods
        IF (abs(fmean) < 1e-7) THEN      ! FMEAN may be 0.
          pm_ival = real(qi5pmx) * 1.  ! then, assume FMEAN = 1.
        ELSE
          pm_ival = real(qi5pmx) * (1. / fmean)
        END IF
        !
      ELSE IF (qi5pmx .LT. 0) THEN
        ! 3) < 0: mix phase based on dominant wave breaking probability bT
        ! Calc bT
        wbt   = calc_wbtv2(a, cg, wn, qr5dpt, u10, udir) ! [0, 1.]
        ! Mix phase by every 1/bT periods
        ! Odin used bT < 1/15. (0.066) → BTLOW = 15 and PM_IVAL > 150 s
        btinv = max(btlow, min(1./max(1e-6, wbt), bthgh))
        IF (abs(fmean) < 1e-7) THEN      ! FMEAN may be 0.
          pm_ival = btinv * 1.  ! then, assume FMEAN = 1.
        ELSE
          pm_ival = btinv * (1. / fmean)
        END IF
      END IF
      !
      ! Previous phase mixing time (relat. to TBEG)
      ! QR5TMIX has already been initialized in w3wdatmd as zero.
      pm_prev = qr5tmix(jsea)
      ! Update t1 if necessary
      pm_delt = t1rel - pm_prev
      IF (pm_delt .GE. pm_ival) THEN
        qr5tmix(jsea) = t1rel        ! relat. to TBEG → PM_PREV
        t1rel         = 0.
      ELSE
        t1rel         = pm_delt
      END IF
#ifdef W3_TS
      IF (iaproc .EQ. napout) THEN
        WRITE(screen, '(F9.1)') t1rel
        IF (qi5pmx .LT. 0 ) WRITE(screen, '(A, F6.3)') '↔ bT: ', wbt
      ENDIF
#endif
    END IF
    !
    ! Calc. Cvk1 from A (C(\bm{k}) = g N(k, θ) / k)
    DO ik = 1, nk
      DO ith = 1, nth
        ispec = ith + (ik-1) * nth
        cvk1(ispec) = a(ith, ik) / wn(ik) * grav
      END DO
    END DO
    !
    ! CalcQRSNL(nk, nth, sig, th, t0, t1, Cvk0, Cvk1, Inpqr0, Snl, Dnl, Kurt)
    ! Depth is needed for reading in kernels at the first run
    qr_depth = qr5dpt
    CALL calcqrsnl(nk, nth, sig(1:nk), th,     &
         t0rel, t1rel, cvk0, cvk1,   &
         inpqr0, snl, dnl, kurt)
    !
    ! Tranform back from C(k) to N(k)
    ! TODO D(ITH, IK) (See NL2 for reference)
    d = 0.0
    DO ik = 1, nk
      DO ith = 1, nth
        ispec = ith + (ik-1) * nth
        s(ith, ik) = snl(ispec) * wn(ik) / grav
      END DO
    END DO
    !
    ! Store wave info. @ t1 → t0
    qr5tim0(jsea)    = t0rel
    qr5cvk0(:, jsea) = cvk0
    qc5int0(:, jsea) = inpqr0
    !
    ! Point output (Snl term)
    ! First ouput action (Find nearest grid points & generate binary files)
    IF (fstout) THEN
      CALL inpout
      fstout = .false.
      IF (iaproc .EQ. napout) THEN
        WRITE(screen, *)
        WRITE(screen, '(A)')                          &
             ' ⊚ → [WW3 SNL₅] Point ouptut initialization'
        WRITE(screen, '(A, I4)')                      &
             ' ⊚ → [WW3 SNL₅] # of valid points: ', nsel
        WRITE(screen, *)
      END IF
    END IF
    !
    ! Calc FACTOR used for Jacobian tranformation from N(k, θ) to E(f, θ)
    factor = tpi / cg * sig(1:nk)
    !
    ! Regular grid & curvilinear grid
    IF ( ((gtype .EQ. rlgtype) .OR. (gtype .EQ. clgtype)) &
         .AND. flout(2) .AND. nsel .GT. 0) THEN
      tdel1 = dsec21(t1abs, tosnl5)
      tdel2 = dsec21(t1abs, tolast(:, 2)) ! not really useful since
      ! TOSNL5 can never catch
      ! TOLAST
      ! Output time
      IF (abs(tdel1) < 1e-6 .OR. abs(tdel2) < 1e-6) THEN
        ! JSEA→ ISEA
        CALL init_get_isea(isea, jsea)
        ! Find the loc of ISEA at PSEA (nearest sea grid point)
        IF (ALLOCATED(pdiff)) DEALLOCATE(pdiff); ALLOCATE(pdiff(nsel))
        pdiff = abs(psea(1:nsel) - isea)
        IF (any(pdiff .EQ. 0)) THEN
          jloc = minloc(pdiff, 1)
#ifdef W3_TS
          IF (iaproc .EQ. napout)                               &
               WRITE(screen, '(3A, I10.8, I7.6)')                    &
               '✓ Point output for |', pnms(jloc), '| @', t1abs
#endif
 
          ! N(θ, k) →  F(f, θ) & S(θ, k) →  S(f, θ)
          DO ith = 1, nth
            a2(:, ith) = a(ith, :) * factor
            s2(:, ith) = s(ith, :) * factor
          END DO
          ! NaN Check
          IF (hasnan(nk, nth, a2) .OR. hasnan(nk, nth, s2)) THEN
            IF (iaproc .EQ. napout)                            &
                 WRITE(screen, *) '★★★ Warning: find NaN in E(f, θ) &
                 & or Snl(f, θ) !'
          END IF
          ! unit no.
          iunt = 500 + jloc
          ! Store data (binary)
          !                 OPEN(IUNT, FILE='NL5_'//trim(PNMS(JLOC))//'_src.bin',  &
          !                      form='unformatted', convert=file_endian, ACCESS='stream',             &
          !                      STATUS='old', POSITION='append', ACTION='write')
          !                 WRITE(IUNT) T1ABS
          !                 WRITE(IUNT) KURT
          !                 WRITE(IUNT) A2
          !                 WRITE(IUNT) S2
          !                 CLOSE(IUNT)
          ! Store data (ascii)
          OPEN(iunt, file='NL5_'//trim(pnms(jloc))//'_src.dat', &
               form='formatted', status='old',                  &
               position='append', action='write')
          WRITE(iunt, '(I10.8, I7.6)') t1abs
          WRITE(iunt, '(ES11.3)')      kurt
          WRITE(iunt, 113) a2
          WRITE(iunt, 113) s2
          CLOSE(iunt)
          !
        END IF
      END IF
    END IF
    ! Format
113 FORMAT ((10es11.3))
    !/
    !/ End of W3SNL5 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3snl5
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE insnl5
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           Q. Liu                  |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         27-Feb-2019 |
    !/                  +-----------------------------------+
    !/
    !/    27-Feb-2019 : Origination.                        ( version 7.13 )
    !/                                                      ( Q. Liu )
    !/
    !  1. Purpose :
    !
    !     Initialization for the GKE module (Prepare wavenumber grid & kernel
    !     coefficients)
    !
    !  2. Method :
    !     See subrs. PrepKGrid & PrepKernelIO of gkeModule.f90
    !
    !  3. Parameters :
    !
    !  4. Subroutines used :
    !     ----------------------------------------------------------------
    !      Name          Type  Module      Description
    !     ----------------------------------------------------------------
    !      STRACE        Subr. W3SERVMD    Subroutine tracing.
    !      PrepKernelIO  Subr. gkeModule   KGrid & Kernel Coeff.
    !
    !  5. Called by :
    !     ----------------------------------------------------------------
    !      Name          Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOGR        Subr. W3IOGRMD Model definition file management.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gkemd,  ONLY: qr_depth, qr_oml, qi_disc, qi_kev, qi_nnz,  &
         qi_interp, prepkernelio
    USE w3gdatmd, ONLY: nk, nth, sig, th,                           &
         qr5dpt, qr5oml, qi5dis, qi5kev, qi5nnz,     &
         qi5ipl, qi5pmx
    USE w3odatmd, ONLY: iaproc, napout, screen
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: ient = 0
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INSNL5')
#endif
    !
    ! Set important parameters for GKE module (QR[I]5DPT/OML/DIS/KEV are
    ! defined in ww3_grid.inp, and QI5NNZ is not known yet)
    qr_depth = qr5dpt
    qr_oml   = qr5oml
    qi_disc  = qi5dis
    qi_kev   = qi5kev
    qi_interp= qi5ipl
    !
    ! Prepare (kx, ky) grid & kernel coefficients
    CALL prepkernelio(nk, nth, sig(1:nk), th, 'WRITE')
    !
    ! Store qi_NNZ to QI5NNZ (which will be used to initialize the
    ! QC5INT0 array)
    qi5nnz   = qi_nnz
    !
    ! Q. Liu (TODO)
    IF (iaproc .EQ. napout) THEN
      WRITE(screen, '(A, F6.1)') " ⊚ → [WW3 SNL₅]: water depth   : ", qr_depth
      WRITE(screen, '(A, F7.2)') " ⊚ → [WW3 SNL₅]: ω λc cut off  : ", qr_oml
      WRITE(screen, '(A, I4)'  ) " ⊚ → [WW3 SNL₅]: Discretiza.   : ", qi_disc
      WRITE(screen, '(A, I4)'  ) " ⊚ → [WW3 SNL₅]: GKE version   : ", qi_kev
      WRITE(screen, '(A, I12)' ) " ⊚ → [WW3 SNL₅]: # of quartets : ", qi_nnz
      WRITE(screen, '(A, I4)'  ) " ⊚ → [WW3 SNL₅]: interpol.     : ", qi_interp
      WRITE(screen, '(A, I4)'  ) " ⊚ → [WW3 SNL₅]: phase mixing  : ", qi5pmx
    END IF
    !/
    !/ End of INSNL5 ----------------------------------------------------- /
    !/
  END SUBROUTINE insnl5
  !/ ------------------------------------------------------------------- /
  FUNCTION calc_wbtv2 (A, CG, WN, DPT, U10, UDIR)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           Q. Liu                  |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         24-Apr-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Aug-2018 : Origination. (w3iogomd.ftn)         ( version 6.06 )
    !/                  Used for output parameter b_T       ( Q. Liu       )
    !/
    !/    24-Apr-2019 : Simplified for NL5                  ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/
    !  1. Purpose :
    !
    !     Estimate the dominant wave breaking probability b_T based on
    !     the empirical parameterization proposed by Babanin et al. (2001).
    !     From their Fig. 12, we have
    !
    !         b_T = 85.1 * [(εp - 0.055) * (1 + H_s/d)]^2.33,
    !
    !     where ε is the significant steepness of the spectral peak, H_s is
    !     the significant wave height, d is the water depth.
    !
    !     For more details, please see
    !         Banner et al.  2000: JPO,      30,  3145 -  3160.
    !         Babanin et al. 2001: JGR, 106(C6), 11569 - 11676.
    !
    !     See subr. CALC_WBT in w3iogomd.ftn for more details.
    !
    !/ ------------------------------------------------------------------- /
    USE w3dispmd, ONLY: wavnu1
    USE w3gdatmd, ONLY: nk, nth, sig, esin, ecos, dth, dsii
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    IMPLICIT NONE
    !
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)     :: a(nth, nk)       ! N(θ, k)
    REAL, INTENT(IN)     :: cg(nk)           ! Cg(k)
    REAL, INTENT(IN)     :: wn(nk)           ! WN(k)
    REAL, INTENT(IN)     :: dpt              ! water depth
    REAL, INTENT(IN)     :: u10              ! wind velocity
    REAL, INTENT(IN)     :: udir             ! wind dirc. (φ in rad)
    REAL                 :: calc_wbtv2
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: ient = 0
#endif
    !
    REAL, PARAMETER      :: beta  = 1.2
    !
    INTEGER              :: ik, ith
    REAL                 :: sinu, cosu, tc, tforce
    REAL                 :: esig(nk) ! E(σ)
    REAL                 :: factor, et, hs, etp, hsp, sigp, kp, &
         cgp, wstp, twbt
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'CALC_WBTv2')
#endif
    !
    ! Wind info. is required to select wind sea partition from the wave
    ! spectrum.
    !
    ! Following Janssen et al. (1989) and Bidlot (2001), spectral components
    ! are considered to be subject to local wind forcing when
    !
    !          c / [U cos(θ - φ)] < β,
    !
    ! where c is the phase velocity c = σ/k, φ is the wind direction, U is
    ! the wind speed U10, (sometimes approximated by U10≅ 28 * ust), β is
    ! the constant forcing parameter with β∈ [1.0, 2.0]. By default, we use
    ! β = 1.2 (Bidlot 2001).
    !
    sinu  = sin(udir)                    ! sinφ
    cosu  = cos(udir)                    ! cosφ
    !
    esig  = 0.                           ! E(σ)
    et    = 0.                           ! ΣE(σ)δσ
    etp   = 0.                           ! ΣE(σ)δσ at peak only
    !
    DO ik = 1, nk
      tc     = sig(ik) / wn(ik)        ! phase velocity c=σ/k
      factor = sig(ik) / cg(ik)        ! σ / cg
      factor = factor * dth            ! σ / cg * δθ
      !
      DO ith = 1, nth
        tforce = tc - u10 * (cosu*ecos(ith)+sinu*esin(ith)) &
             * beta
 
        IF (tforce .LT. 0.) THEN ! wind sea component
          esig(ik) = esig(ik) + a(ith, ik) * factor
        ENDIF
      ENDDO ! ITH
      !
    ENDDO ! IK
    !
    ! ESIG is E(σ) of the wind sea after filtration of any background swell.
    ! Now we need to get Hs & σp for the wind sea spectrum.
    ! Unlike w3iogomd.ftn, the tail energy is not added here.
    et = sum(esig * dsii)
    hs = 4. * sqrt(max(0., et))
    !
    ! Get σp from E(σ)
    ! FPOPT = 0 in w3iogomd.ftn: fp defined by Young (1999, p. 239)
    sigp = sum(esig**4. * sig(1:nk) * dsii) /  &
         max(1e-10, sum(esig**4. * dsii))
    IF (abs(sigp) < 1e-7) sigp = sig(nk) ! σp = 0
    !
    ! kp from σp (linear dispersion)
    CALL wavnu1 (sigp, dpt, kp, cgp)
    !
    !                         { /1.3σp         }1/2
    ! peak wave height Hp = 4 { |      E(σ) dσ }
    !                         { /0.7σp         }
    !
    DO ik = 1, nk
      IF ( (sig(ik) >= 0.7 * sigp) .AND. &
           (sig(ik) <= 1.3 * sigp) ) THEN
        etp = etp + esig(ik) * dsii(ik)
      ENDIF
    ENDDO ! IK
    hsp  = 4. * sqrt(max(0., etp))
    !
    ! significant steepness of the peak region εp
    !
    wstp = 0.5 * kp * hsp
    !
    ! Dominant wave breaking b_T
    !
    twbt = 85.1 * (max(0.0, wstp - 0.055) * (1 + hs/dpt))**2.33
    twbt = min(1.0, twbt)
    !
    calc_wbtv2 = twbt
 
    RETURN
    !/
    !/ End of  CALC_WBTv2 ------------------------------------------------ /
    !/
  END FUNCTION calc_wbtv2
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE inpout
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           Q. Liu                  |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         25-Mar-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Mar-2019 : Origination.                        ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/    27-Apr-2019 : Add the ascii option                ( Q. Liu       )
    !/
    !  1. Purpose :
    !
    !     Initialization for point output (Snl) [see also W3IOPP of w3iopomd]
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !  4. Subroutines used :
    !
    !  5. Called by :
    !     ----------------------------------------------------------------
    !      Name          Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SNL5        Subr. W3SNL5MD S_{nl} GKE
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: tpi
    USE w3gdatmd,  ONLY: nk, nth, sig, th, qr5dpt,         &
         flagll, xgrd, ygrd, mapsta, mapfs
    USE w3odatmd,  ONLY: nopts, ptnme, ptloc, iptint,      &
         iaproc, napout, screen
    USE w3servmd,  ONLY: dist_sphere
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER           :: ixs(4), iys(4), ix, iy, ipt, is, &
         jloc, jx, jy, isea, smap(4), iunt
    REAL              :: plon, plat, xlon, ylat, dist(4)
    !/ ------------------------------------------------------------------- /
    !/
    ! Initialize arrays
    IF (ALLOCATED(psea)) DEALLOCATE(psea); ALLOCATE(psea(nopts))
    IF (ALLOCATED(pnms)) DEALLOCATE(pnms); ALLOCATE(pnms(nopts))
    !
    nsel    = 0
    psea(:) = 0
    pnms(:) = 'null'
    dist(:) = -999.
    !
    DO ipt = 1, nopts
      ! Get lon & lat of this output point
      plon   = ptloc(1, ipt)
      plat   = ptloc(2, ipt)
      ! Get four indices surrounding the output point
      ixs(:) = iptint(1, :, ipt)
      iys(:) = iptint(2, :, ipt)
      DO is = 1, 4
        ! Get lon & lat of four corner points
        ix   = ixs(is)
        iy   = iys(is)
        xlon = xgrd(iy, ix)
        ylat = ygrd(iy, ix)
        ! Grid point status
        IF (mapsta(iy, ix) .EQ. 0) cycle
        ! Calc dist.
        IF (flagll) THEN
          dist(is) = dist_sphere(plon, plat, xlon, ylat)
        ELSE
          dist(is) = sqrt((plon - xlon)**2. + (plat - ylat)**2.)
        END IF
      END DO
      ! A sea point filter: there must be at least one sea grid point around
      ! the selected output location. [maybe not necessary since IOPP already
      ! checked this criterion]
      !
      IF (all(dist < 0.)) cycle
      ! Find the nearest sea grid point
      jloc = minloc(dist, 1, dist >= 0.)
      jx   = ixs(jloc)
      jy   = iys(jloc)
      isea = mapfs(jy, jx)
      ! Basic check
#ifdef W3_TS
      IF (flagll) THEN
        IF (iaproc .EQ. napout)                                    &
             WRITE(screen, "(A, 2F10.3, A, 2F10.3, A)")                 &
             '✗ (PLON, PLAT): (', plon, plat, ') | (XGRD, YGRD): (',&
             xgrd(jy, jx), ygrd(jy, jx), ')'
      ELSE
        IF (iaproc .EQ. napout)                                    &
             WRITE(screen, "(A, 2E10.3, A, 2E10.3, A)")                 &
             '✗ (PLON, PLAT): (', plon, plat, ') | (XGRD, YGRD): (',&
             xgrd(jy, jx), ygrd(jy, jx), ')'
      END IF
#endif
      ! Store ISEA
      nsel = nsel + 1
      psea(nsel) = isea
      pnms(nsel) = ptnme(ipt)
      ! Store Unit (Open & Write Binary files)
      iunt = 500 + nsel
      ! Binary
      !         OPEN(IUNT, FILE='NL5_'//trim(PNMS(NSEL))//'_src.bin',       &
      !              form='unformatted', convert=file_endian, ACCESS='stream', STATUS='replace', &
      !              ACTION='write')
      !         WRITE(IUNT) PLON, PLAT
      !         WRITE(IUNT) XGRD(JY, JX), YGRD(JY, JX)
      !         WRITE(IUNT) QR5DPT
      !         WRITE(IUNT) NK, NTH
      !         WRITE(IUNT) SIG(1:NK)/TPI  ! f, θ
      !         WRITE(IUNT) TH
      !         CLOSE(IUNT)
      ! Ascii
      OPEN(iunt, file='NL5_'//trim(pnms(nsel))//'_src.dat',       &
           form='formatted', status='replace', action='write')
      WRITE(iunt, '(2ES11.3)') plon, plat
      WRITE(iunt, '(ES11.3)' ) qr5dpt
      WRITE(iunt, '(2I5)')     nk, nth
      WRITE(iunt, 113) sig(1:nk)/tpi  ! f, θ
      WRITE(iunt, 113) th
      CLOSE(iunt)
      !
    END DO
    ! Format
113 FORMAT ((10es11.3))
    !
  END SUBROUTINE inpout
  !/ ------------------------------------------------------------------- /
 
  FUNCTION hasnan(NK, NTH, ARR2D)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           Q. Liu                  |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         25-Apr-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Apr-2019 : Origination.                        ( version 7.13 )
    !/                                                      ( Q. Liu       )
    !/
    !  1. Purpose :
    !     Check if the 2D array `ARR2D` contains NaN (see also w3gsrumd.ftn)
    !/
    IMPLICIT NONE
    !
    INTEGER, INTENT(IN)  :: nk, nth        ! # OF FREQ. & DIRC.
    REAL, INTENT(IN)     :: arr2d(nk, nth)
    LOGICAL              :: hasnan
    !/
    hasnan = .true.
    !
    IF ( all(arr2d .GE. -huge(arr2d(1, 1))) .AND.  &
         all(arr2d .LE.  huge(arr2d(1, 1))) ) THEN
      hasnan = .false.
    END IF
    !
    RETURN
    !/
  END FUNCTION hasnan
  !/ ------------------------------------------------------------------- /
  !/
  !/ End of module W3SNL5MD -------------------------------------------- /
  !/
END MODULE w3snl5md
!/ ------------------------------------------------------------------- /