WW3/w3gkemd_8F90_source.html

#include "w3macros.h"

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

module w3gkemd

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

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |          Odin Gramstad            |

  !/                  |          Qingxiang Liu            |

  !/                  |                                   |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         03-Jun-2021 |

  !/                  +-----------------------------------+

  !/

  !/    26-May-2014 : Origination.                        ( version 3.14 )

  !/                  Intially Dr. O. Gramstad implemented his GKE method

  !/                  in WW3 v3.14 which worked well for the single-grid

  !/                  point duration-limited test.

  !/

  !/    09-Nov-2018 : Fully implemented in WW3            ( version 7.13 )

  !/                                                      ( Q. Liu       )

  !/    16-Apr-2019 : Add save attribute explicitly       ( version 7.13 )

  !/                                                      ( Q. Liu       )

  !/    18-Apr-2019 : Add the bilinear interp. option     ( version 7.13 )

  !/                                                      ( Q. Liu       )

  !/    08-Jul-2019 : Use kind=8 for qi_nnz               ( Q. Liu       )

  !/    01-Apr-2020 : Boundary conditions                 ( Q. Liu       )

  !/    03-Jun-2021 : Merge into the WW3 Github           ( version 7.12 )

  !/                                                      ( Q. Liu       )

  !/

  !  1. Purpose:

  !     Calculate the (resonant & quasi/near-resonant) four-wave nonlinear

  !     interaction term S_{nl} according to the generalized kinetic

  !     equation (GKE) developed in Gramstad and Stiassnie (2013).

  !

  !     References:

  !     Gramstad and Stiassnie (2013), JFM, 818, 280-303 (hereafter GS13)

  !     Gramstad and Babanin (2016),    OD,  66, 509-526 (hereafter GB16)

  !     Liu et al. (2021),             JFM, 910, A50     (hereafter LGB21)

  !     &

  !     Annenkov and Shrira (2006),    JFM, 561, 181-207 (*)

  !     Annenkov and Shrira (2015),    JPO,  45, 807-812

  !     Annenkov and Shrira (2018),    JFM, 844, 766-795 (*)

  !     &

  !     Shrira and Annenkov (2013),    Book Ch., 239-281

  !     Annenkov and Shrira (2016),    Book Ch., 159-178

  !

  !     (*) Note that equations therein contain typos.

  !

  !  2. Subroutines and functions :

  !

  !     [Part 1]: Kernel Function

  !

  !     Calculate the kernel function T_{0, 1, 2, 3} for the Zakharov

  !     Equation.

  !

  !     References:

  !     Krasitskii (1994), JFM, 272, 1 - 20 (hereafter K94)

  !     Janssen    (2009), JFM, 637, 1 - 44 (hereafter J09)

  !     Mei et al. (2005), ch. 14, 882 - 884

  !

  !     Based on my own observation, Odin has closely followed the

  !     equations presented in the appendix (A.1/2) of J09.

  !

  !     ----------------------------------------------------------------

  !      Name                Type  Scope    Description

  !     ----------------------------------------------------------------

  !      QFunc               Func. Private  q = ω^2 / g

  !      VpFunc              Func. Private  V^{(+)}_{1, 2, 3}

  !      VmFunc              Func. Private  V^{(-)}_{1, 2, 3}

  !      UFunc               Func. Private  U_{1, 2, 3, 4}

  !      TFunc               Func. Public   T_{1, 2, 3, 4}

  !     ----------------------------------------------------------------

  !

  !     [Part 2]: Find Quartets (total number & configurations)

  !

  !     References:

  !     Annenkov and Shrira (2015),       JPO, 45, 807-812

  !     Hasselmann and Hasselmann (1981), Exact-NL/DIA report

  !     Hasselmann and Hasselmann (1985), JPO, 15, 1369-1377

  !

  !     ----------------------------------------------------------------

  !      Name                Type  Scope    Description

  !     ----------------------------------------------------------------

  !      FindQuartetNumber   Subr. Private  Total No. of Quartets

  !      FindQuartetConfig   Subr. Private  Config.   of Quartets

  !

  !      [Part 3]: Sparse matrix (storage, operation)

  !

  !      References:

  !      Saad (1994) SPARSKIT: a basic tool kit for sparse matrix

  !                            compuation (version 2)

  !

  !     ----------------------------------------------------------------

  !      Name                Type  Scope    Description

  !     ----------------------------------------------------------------

  !      CooCsrInd           Subr. Private  COO to CSR format

  !      ASymSmatTimVec      Subr. Private  (A±A^T)∙X, where (A±A^T) is an

  !                                         (anti)symmetric sparse matrix

  !

  !      [Part 4]: GKE Integral (main subrs.)

  !

  !      References:

  !      Gramstad and Stiassnie (2013),    JFM, 818, 280-303 (hereafter GS13)

  !      Gramstad and Babanin (2016),       OD,  66, 509-526 (hereafter GB16)

  !      Liu et al. (2021),                JFM, 910, A50     (hereafter LGB21)

  !      Janssen (2003),                   JPO,  33, 863-884 (hereafter J03)

  !      Janssen (2009),                   JFM, 637,   1- 44 (hereafter J09)

  !      Annenkov and Shrira (2013),       JFM, 726, 517-546

  !

  !      Hasselmann and Hasselmann (1981), Exact-NL/DIA report

  !      Hasselmann and Hasselmann (1985), JPO, 15, 1369-1377

  !      van Vledder (2006),               CE,  53, 223-242

  !      Tolman (2013),                    OM,  70,  11- 24

  !

  !     ----------------------------------------------------------------

  !      Name                Type  Scope    Description

  !     ----------------------------------------------------------------

  !      PrepKGrid           Subr. Private  (σ, θ) to (kx, ky)

  !      PrepKernelIO        Subr. Private  Read/Write Quartet Cfg file

  !      BiInterpWT          Subr. Private  Calc. interp. weights

  !      CalcQRSNL           Subr. Public   GKE Transfer Integral

  !

  !  3. Future work (TODO)

  !   * The current version only works for a constant-depth application (

  !     either deep or finite deep). Extension of this module to be

  !     applicable to varying-depth cases may be pursued in the future.

  !

  !   * Dnl   -- diagonal term

  !   * βnpqr -- nonlinear stokes correction

  !/

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

  !

  ! Public parameters

  !

  ! * `qi_` denotes the variable is integer number

  ! * `qr_` ...                     real    ...

  ! * `qs_` ...                     string

  ! * `ql_` ...                     logical

  ! * `qc_` ...                     complex

  !/

  implicit none

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

  public  :: prepkernelio, calcqrsnl

  !

  public  :: qr_depth, qr_oml, qi_disc, qi_kev, qi_nnz, qi_interp


  private :: qfunc, vpfunc, vmfunc, ufunc, tfunc,           &

       findquartetnumber, findquartetconfig,          &

       coocsrind, asymsmattimvec,                     &

       prepkgrid, biinterpwt

  !

  private :: qs_ver, qi_lrb, qr_eps, qr_grav,               &

       qr_pi, qr_tpi, qr_dmax, qc_iu, qs_cfg,         &

       qr_kx, qr_ky, qr_dk, qr_om, qi_nrsm,           &

       qi_nn, qi_pp, qi_qq, qi_rr,                    &

       qr_k4x, qr_k4y, qr_om4, qr_dom,                &

       qr_tkern, qr_tkurt,                            &

       qi_iccos, qi_ircsr, qr_sumqr, qr_sumnp,        &

       qi_bind, qr_bwgh, qr_wn1

  !

  private :: qi_bound, qr_fpow, qr_bdry

  !

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

  real                       :: qr_depth                    ! Real water depth d (m)


  real                       :: qr_oml                      ! λ cut off factor

  ! λ ≤ 0 →  quartets far

  ! from resonance will not

  ! be excluded.


  integer                    :: qi_disc                     ! Discretization of GKE

  ! 0: continuous (like Exact-NL, WRT)

  ! 1: discrete   (see GS13)


  integer                    :: qi_kev                      ! Version of KE

  ! 0: GKE

  ! 1: KE from J03


  integer                    :: qi_interp                   ! Interp. option

  ! 0: Nearest bin

  ! 1: Bilinear Interp


  integer, parameter         :: qi_bound= 1                 ! Boundary conditions

  ! 0: no bound

  ! 1: tail extension


  real, parameter            :: qr_fpow= -5.                ! E(f) tail power law

  !

  character(len=50), parameter                              &

       :: qs_ver  = 'gkev0'           ! version number/str

  integer, parameter         :: qi_lrb  = 4                 ! 4 bytes

  real, parameter            :: qr_eps  = epsilon(100.0)    ! Smallest positive

  ! value supported by the

  ! compiler (e.g., gfortran

  ! → 1.19E-7)


  real, parameter            :: qr_grav = 9.806             ! Gravational acc (m/s^2)

  real, parameter            :: qr_pi   = 3.141592653589793 ! π

  real, parameter            :: qr_tpi  = 2 * qr_pi         ! π * 2

  real, parameter            :: qr_dmax = 3000.0            ! Maximum allowed water

  complex, parameter         :: qc_iu   = (0.0, 1.0)        ! complex unit `i`

  !

  character(len=100)         :: qs_cfg                      ! File name for quartet/kernel

  !

  real, allocatable, save    :: qr_kx(:), qr_ky(:),   &     ! kx, ky (2D grid → 1D vector)

       qr_dk(:), qr_om(:)          ! Δ\vec{k}, ω,

  !

  integer(kind=8)            :: qi_nnz                      ! # of quartets

  integer                    :: qi_nrsm                     ! # of rows of SMat

  integer, allocatable, save :: qi_nn(:), qi_pp(:),   &     ! Index for Quartets

       qi_qq(:), qi_rr(:)

  real, allocatable, save    :: qr_k4x(:), qr_k4y(:), &     ! kx, ky, ω for 4th wave

       qr_om4(:)

  real, allocatable, save    :: qr_dom(:),            &     ! Δω

       qr_tkern(:),          &     ! Kernel `T`

       qr_tkurt(:)                 ! Kurtosis `T`

  integer, allocatable, save :: qi_iccos(:),          &     ! col index of CooCsr

       qi_ircsr(:)                 ! row begining index of

  ! Csr sparse matrix

  real, allocatable, save    :: qr_sumqr(:),          &     ! Σ over Q, R

       qr_sumnp(:, :)              ! Σ over P

  !

  integer, allocatable, save :: qi_bind(:, :)               ! Bilinear interp. (index and

  real, allocatable, save    :: qr_bwgh(:, :)               ! weight)

  real, allocatable, save    :: qr_bdry(:)                  ! Boundary weight

  real, allocatable, save    :: qr_wn1(:)                   ! wavenumber k(nk)

  !/

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

contains

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

  !/ [Part 1]

  !/

  function qfunc(kx, ky)

    !/

    !/    19-Dec-2011 : Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/

    !  1. Purpose: Define q = ω^2 / g (i.e., q in K94 & J09)

    !

    !  2. Method:

    !     For wind-generated ocean surface waves, the dispersion relation

    !     reads

    !         ω^2 = g k tanh(kd),

    !     where g is the gravtional acceleration, ω is the radian frequency,

    !     d is the water depth. Hence,

    !

    !             / k = √(k_x**2. + k_y**2.) for deep-water

    !         q = |

    !             \ k tanh(kd)               for finite-deep water (e.g.,

    !                                            d < 2500.)

    !

    !/

    implicit none

    !

    real, intent(in) :: kx, ky   ! x, y components of wavenumber

    ! vector (kx, ky)

    real             :: qfunc    ! Returned function

    !/

    qfunc = sqrt(kx*kx + ky*ky)  ! deep-water case (q = k)

    !

    ! Odin used qr_dmax = 2500.

    !

    if (qr_depth > qr_eps .and. qr_depth < qr_dmax) then

      qfunc = qfunc * tanh(qfunc * qr_depth)  ! finite-deep

    end if

    !

    return

    !/

  end function qfunc

  !/

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

  !/

  function vpfunc(k0x, k0y, k1x, k1y, k2x, k2y)

    !/

    !/    19-Dec-2011: Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/

    !/

    !  1. Purpose:

    !     Calculate the second-order coefficient V^{(+)}_{1, 2, 3} of J09,

    !     which corresponds to U^{(3)}_{0, 1, 2} of K94.

    !

    !     ◆ V^{(+)}_{1, 2, 3} differs from U^{(3)}_{0, 1, 2} by a factor

    !       of 1/2π --- this is because the wave spectrum F(k) used in K94

    !       and J09 differ by a fator of (1/2π)^2.

    !

    !/

    implicit none

    !

    real, intent(in) :: k0x, k0y, k1x, k1y, k2x, k2y ! 3 waves

    real             :: vpfunc                       ! V^{(+)}_{1, 2, 3}

    !

    real             :: q0, q1, q2                   ! q for 3 waves

    !/

    ! Call Q function here

    q0 = qfunc(k0x, k0y)

    q1 = qfunc(k1x, k1y)

    q2 = qfunc(k2x, k2y)

    !

    ! Odin has ignored √g here because it will be absorbed/vanish when we

    ! calculate the kernel function T. I, however, included √g here for

    ! clarity.

    ! V^{(+)}_{1, 2, 3}

    !

    vpfunc = sqrt(1.0/32.0) * (                                       &

         (k0x*k1x + k0y*k1y + q0*q1) * sqrt(sqrt(qr_grav*q2 / (q0*q1)))&

         + (k0x*k2x + k0y*k2y + q0*q2) * sqrt(sqrt(qr_grav*q1 / (q0*q2)))&

         + (k1x*k2x + k1y*k2y + q1*q2) * sqrt(sqrt(qr_grav*q0 / (q1*q2))))

    !

    return

    !/

  end function vpfunc

  !/

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

  !/

  function vmfunc(k0x, k0y, k1x, k1y, k2x, k2y)

    !/

    !/    19-Dec-2011 : Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/

    !  1. Purpose:

    !     Calculate the second-order coefficient V^{(-)}_{1, 2, 3} of J09,

    !     which corresponds to U^{(1)}_{0, 1, 2} of K94.

    !

    !     ◆ V^{(-)}_{1, 2, 3} differs from U^{(1)}_{0, 1, 2} by a factor

    !       of 1/2π

    !

    !/

    implicit none

    !

    real, intent(in) :: k0x, k0y, k1x, k1y, k2x, k2y ! 3 waves

    real             :: vmfunc                       ! V^{(-)}_{1, 2, 3}

    !

    real             :: q0, q1, q2                   ! q for 3 waves

    !/

    ! Call Q function here

    q0 = qfunc(k0x, k0y)

    q1 = qfunc(k1x, k1y)

    q2 = qfunc(k2x, k2y)

    !

    ! V^{(-)}_{1, 2, 3}

    !

    vmfunc = sqrt(1.0/32.0) * (                                       &

         (k0x*k1x + k0y*k1y - q0*q1) * sqrt(sqrt(qr_grav*q2 / (q0*q1)))&

         + (k0x*k2x + k0y*k2y - q0*q2) * sqrt(sqrt(qr_grav*q1 / (q0*q2)))&

         + (k1x*k2x + k1y*k2y + q1*q2) * sqrt(sqrt(qr_grav*q0 / (q1*q2))))

    !

    return

    !/

  end function vmfunc

  !/

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

  !/

  function ufunc(k0x, k0y, k1x, k1y, k2x, k2y, k3x, k3y)

    !/

    !/    19-Dec-2011 : Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/

    !  1. Purpose:

    !     Calculate the intermediate quantity (i.e., U_{1, 2, 3, 4} in J09,

    !     V_{0, 1, 2, 3} in K94) for the third-order coefficient (i.e.,

    !     W^{(2)}_{1, 2, 3, 4} in J09, V^{(2)}_{0, 1, 2, 3} in K94).

    !

    !     ◆ U_{1, 2, 3, 4} differs from V_{0, 1, 2, 3} by a factor of

    !       (1/2π)^2.

    !

    !/

    implicit none

    !

    real, intent(in) :: k0x, k0y, k1x, k1y,       &

         k2x, k2y, k3x, k3y        ! 4 waves

    real             :: ufunc                     ! U_{1, 2, 3, 4}

    !

    real             :: q0, q1, q2, q3            ! q for 4 waves

    !/

    ! Call Q function here

    q0 = qfunc(k0x, k0y)

    q1 = qfunc(k1x, k1y)

    q2 = qfunc(k2x, k2y)

    q3 = qfunc(k3x, k3y)

    !

    ! U_{1, 2, 3, 4}

    !

    ufunc = (1.0/16.0) * sqrt(sqrt(q2*q3 / (q0*q1))) *              &

         (2.0*((k0x*k0x + k0y*k0y) * q1 + (k1x*k1x + k1y*k1y) * q0)-&

         q0*q1*( qfunc(k0x+k2x, k0y+k2y) + qfunc(k1x+k2x, k1y+k2y)+&

         qfunc(k0x+k3x, k0y+k3y) + qfunc(k1x+k3x, k1y+k3y) ))

    !

    return

    !/

  end function ufunc

  !/

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

  !/

  function tfunc(k0x, k0y, k1x, k1y, k2x, k2y, k3x, k3y)

    !/

    !/    19-Dec-2011 : Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/

    !  1. Purpose:

    !     Calculate the Kernel function for the four-wave interaction, i.e.,

    !     (T_{1, 2, 3, 4}, \widetilde{V}^{(2)}_{0, 1, 2, 3} in K94).

    !     ◆ T from J09 and K94 differ by a factor of (1/2π)^2.

    !

    !     Odin's comment:

    !     Kernel function for all combination that are not Stokes correction.

    !     I.e. n0 != n2 and n0 != n3

    !/

    implicit none

    !

    real, intent(in) :: k0x, k0y, k1x, k1y, &

         k2x, k2y, k3x, k3y        ! 4 waves

    real             :: tfunc                     ! T_{1, 2, 3, 4}

    !

    ! Virtual-state interaction: two free waves generate a virtual state

    ! consisting of bound waves, which then decays into a different set of

    ! free waves (see J09)

    !

    real             :: om0, om1, om2, om3, &     ! ω for 4 waves

         om02,               &     ! ω_{0-2}

         om13,               &     ! ω_{1-3}

         om12,               &     ! ω_{1-2}

         om03,               &     ! ω_{0-3}

         om0p1,              &     ! ω_{0+1}

         om2p3                     ! ω_{2+3}

    !

    real             :: l14, l23, l56,      &

         w                         ! W^{(2)}_{1, 2, 3, 4} in J09

    ! or

    ! V^{(2)}_{0, 1, 2, 3} in K94

    !/

    ! Initilization

    om0p1 = 0.

    om2p3 = 0.

    w     = 0.

    l14   = 0.

    l23   = 0.

    l56   = 0.

    !

    ! Get ω from q: q = ω^2 / g →  ω = √(qg)

    ! Odin has ignored √g here because it will be absorbed/vanish when we

    ! calculate the kernel function T (V / ω). I, however, included √g here for

    ! clarity.

    !

    ! ω for four free waves

    om0 = sqrt(qr_grav * qfunc(k0x, k0y))

    om1 = sqrt(qr_grav * qfunc(k1x, k1y))

    om2 = sqrt(qr_grav * qfunc(k2x, k2y))

    om3 = sqrt(qr_grav * qfunc(k3x, k3y))

    !

    ! ω for other combined waves

    !

    om02 = sqrt(qr_grav * qfunc(k0x-k2x, k0y-k2y))

    om13 = sqrt(qr_grav * qfunc(k1x-k3x, k1y-k3y))

    om12 = sqrt(qr_grav * qfunc(k1x-k2x, k1y-k2y))

    om03 = sqrt(qr_grav * qfunc(k0x-k3x, k0y-k3y))

    !

    if (abs(k0x+k1x) > qr_eps .or. abs(k0y+k1y) > qr_eps) then

      !           k₀ + k₁ = k₂ + k₃ = 0., ω_{0+1} = 0.,

      !           V^{(-)}_{0+1, 0, 1} ~ 1/ω_{0+1} = NaN, L56 = NaN

      om0p1 = sqrt(qr_grav * qfunc(k0x+k1x, k0y+k1y))

      om2p3 = sqrt(qr_grav * qfunc(k2x+k3x, k2y+k3y))

    end if

    !

    ! W^{(2)}_{1, 2, 3, 4} [Call U function here] for direct interaction

    !

    w = ufunc(-k0x, -k0y, -k1x, -k1y,  k2x,  k2y,  k3x,  k3y) +     &

         ufunc( k2x,  k2y,  k3x,  k3y, -k0x, -k0y, -k1x, -k1y) -     &

         ufunc( k2x,  k2y, -k1x, -k1y, -k0x, -k0y,  k3x,  k3y) -     &

         ufunc(-k0x, -k0y,  k2x,  k2y, -k1x, -k1y,  k3x,  k3y) -     &

         ufunc(-k0x, -k0y,  k3x,  k3y,  k2x,  k2y, -k1x, -k1y) -     &

         ufunc( k3x,  k3y, -k1x, -k1y,  k2x,  k2y, -k0x, -k0y)

    !

    ! First & Fourth lines for virtual-state interaction in J09

    !

    l14 = vmfunc(k0x, k0y, k2x, k2y, k0x-k2x, k0y-k2y) *            &

         vmfunc(k3x, k3y, k1x, k1y, k3x-k1x, k3y-k1y) *            &

         (1.0/(om2 + om02 - om0) + 1.0/(om1 + om13 - om3)) +       &

         vmfunc(k1x, k1y, k3x, k3y, k1x-k3x, k1y-k3y) *            &

         vmfunc(k2x, k2y, k0x, k0y, k2x-k0x, k2y-k0y) *            &

         (1.0/(om3 + om13 - om1) + 1.0/(om0 + om02 - om2))

    !

    ! Second & Third lines for virtual-state interaction in J09

    !

    l23 = vmfunc(k1x, k1y, k2x, k2y, k1x-k2x, k1y-k2y) *            &

         vmfunc(k3x, k3y, k0x, k0y, k3x-k0x, k3y-k0y) *            &

         (1.0/(om2 + om12 - om1) + 1.0/(om0 + om03 - om3)) +       &

         vmfunc(k0x, k0y, k3x, k3y, k0x-k3x, k0y-k3y) *            &

         vmfunc(k2x, k2y, k1x, k1y, k2x-k1x, k2y-k1y) *            &

         (1.0/(om3 + om03 - om0) + 1.0/(om1 + om12 - om2))

    !

    ! Fifth & Sixth lines for virtual-state interaction in J09

    !

    if (abs(k0x+k1x) > qr_eps .or. abs(k0y+k1y) > qr_eps) then

      !           k₁ + k₂ = k₃ + k₄ = 0., ω_{1+2} = 0.,

      !           V^{(-)}_{1+2, 1, 2} ~ 1/ω_{1+2} = NaN, L56 = NaN

      l56 = vmfunc(k0x+k1x, k0y+k1y, k0x, k0y, k1x, k1y) *        &

           vmfunc(k2x+k3x, k2y+k3y, k2x, k2y, k3x, k3y) *        &

           (1.0/(om0p1 - om0 - om1) + 1.0/(om2p3 - om2 - om3)) + &

           vpfunc(-k0x-k1x, -k0y-k1y, k0x, k0y, k1x, k1y) *      &

           vpfunc(-k2x-k3x, -k2y-k3y, k2x, k2y, k3x, k3y) *      &

           (1.0/(om0p1 + om0 + om1) + 1.0/(om2p3 + om2 + om3))

    end if

    !

    ! T_{1, 2, 3, 4}

    !

    tfunc =  w - l14 - l23 - l56

    !

    return

    !/

  end function tfunc

  !/

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

  !/ [Part 2]

  !/

  subroutine findquartetnumber(ns, kx, ky, om, oml, nnz)

    !/

    !/    19-Dec-2011 : Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/    02-Apr-2020 : Boundary conditions (< kmin, > kmax)( Q. Liu      )

    !/

    !  1. Purpose:

    !     Find the total number of quartets (resonant and quasi/near-resonant

    !     four waves) satisfying the criteria below:

    !

    !     1) \vec{k₁} + \vec{k₂} = \vec{k₃} + \vec{k₄}

    !

    !     2) Δω = |ω₁ + ω₂ - ω₃ - ω₄| <= λc \min(ω₁, ω₂, ω₃, ω₄)

    !        - that is, quartets far from the resonance is excluded for

    !          saving the computational cost.

    !

    !     3) For a given 2D frequency-direction grid (k_i, θ_j, i = 1, ...,

    !        NK, j = 1, ..., NTH) consisting of NS points (NS = NK * NTH),

    !        we will first reshape the 2D spectral grid (k_i, θ_j)

    !        into a 1D wavenumber vector (k_{xl}, k_{yl}, l = 1, ..., NS).

    !        Afterwards, we should have

    !        3.a) l₂ >= l₁

    !        3.b) l₃ ≠  l₁, l₃ ≠ l₂ (otherwise the third and fourth wave

    !                                components will be the same as the

    !                                first and second)

    !        3.c) l₄ >= l₃

    !        3.d) l₄ ≠  l₁, l₄ ≠ l₂

    !        3.e) k₄ >= k_{min}, k₄ <= k_{max}

    !

    !        Note that `l` here only denotes the index of a specific wave

    !        component inside the 1D wavenumber vector array. For k₄, its

    !        index l₄ is not exact, and is just approximated by the index

    !        of its closest wave component.

    !

    !     4) If we store the located quartets in a 2D large sparse matrix,

    !        which can be organized as

    !               |K K K |

    !               |∩ ∩ ∩ |

    !               |3 q l₃| 1 2 . . . N 1 2 . . . N . . . . . . 1 2 . . . N

    !               |4 r l₄| 1 1 1 1 1 1 2 2 2 2 2 2 . . . . . . N N N N N N

    !               |∪ ∪ ∪ |

    !        -------

    !        K {1 2}

    !        K {n p}

    !        K {l₁l₂}

    !        -------

    !           1 1

    !           2 1                                  (2,1,N,3) ✗⁴    ✗²(2,1,3,N)

    !           . 1                           col > row  → ▲

    !           . 1

    !           . 1

    !           N 1

    !           1 2                                  (1,2,N,3) ✗³   [★ (1,2,3,N)]

    !           2 2

    !           . 2

    !           . 2                        ⊚ ← row = l₁ + (l₂ - 1) * NS

    !           . 2                        ↑

    !           N 2                        col = l₃ + (l₄ - 1) * NS

    !           . .

    !           . .

    !           . .

    !           . .

    !           . .

    !           . .  (N,3,2,1) ✓⁴        ✓³(N,3,1,2)

    !           1 N              ▼ ← col < row

    !           2 N

    !           . N  (3,N,2,1) ✓²        ☆ (3,N,1,2)

    !           . N

    !           . N

    !           N N

    !

    !        where `N` shown above denotes `NS`, not `NK`, therefore the shape

    !        of this large sparse matrix is (NS*NS, NS*NS).

    !

    !        Only quartets with col > row (highlighted by ▲ , i.e.,

    !                           ---------

    !        elements in the upper trianglar matrix) are selected because,

    !        for example, ★ (1, 2, 3, NS) & ☆ (3, NS, 1, 2) essentially

    !        refer to the same quartet.

    !

    !     To sum up, criteria 1) and 2) are kinetic, whereas 3) and 4) are

    !     enforced to avoid duplicating quartets since

    !     ★  (k₁, k₂, k₃, k₄)

    !

    !     & [symmetric]

    !     ✗² (k₂, k₁, k₃, k₄) ← filterd by 3.a)

    !     ✗³ (k₁, k₂, k₄, k₃) ← filterd by 3.c)

    !     ✗⁴ (k₂, k₁, k₄, k₃) ← filterd by 3.a) and 3.c)

    !

    !     & [antisymmetric]

    !     ☆  (k₃, k₄, k₁, k₂) ← filterd by 4)

    !     ✓² (k₃, k₄, k₂, k₁)

    !     ✓³ (k₄, k₃, k₁, k₂)

    !     ✓⁴ (k₄, k₃, k₂, k₁)

    !

    !     are essentially the same.

    !

    !     ◆ criteria 3.b) and 3.d) exclude two quartets:

    !        / k₁ = k₃, k₂ = k₄ (k₁, k₂, k₁, k₂)

    !        \ k₁ = k₄, k₂ = k₃ (k₁, k₂, k₂, k₁)

    !        →  singular points for the nonlinear transfer integral as

    !           T_{1, 2, 1, 2} or T_{1, 2, 2, 1} ~ 1 / 0 = NaN

    !

    !       van Vledder (2006, p. 231) argued that the first quadruplet had

    !       negligible contribution to the total transfer rate. Similarly,

    !       for the symmetric reason, the contribution from the second

    !       quadruplet is also very limited.

    !

    !     ◆ We should keep in mind that the Snl term for wave component

    !       3 in ★ (i.e., k₃) and in ☆ (i.e., k₁) are the same.

    !

    !     ◆ Although the other 7 quartets are not counted here, their

    !       contributions to the nonlinear transfer rates should not be

    !       ignored as the interval of the 6D integration starts from

    !       -∞ and ends at +∞ !

    !

    !     More details can be found in Appendix of LGB21.

    !

    !     See also references:

    !     Hasselmann and Hasselmann (1981)  Exact-NL/DIA report

    !     Hasselmann and Hasselmann (1985), JPO, 15, 1369 - 1377.

    !     Annenkov and Shrira (2015),       JPO, 45, 807-812

    !

    !/

    implicit none

    !

    integer, intent(in)  :: ns           ! length of 1D wavenumber

    ! vector, ns = nk * nth

    real, intent(in)     :: kx(ns),   &

         ky(ns),   &  ! (kx, ky) components

         om(ns)       ! ω or σ

    real, intent(in)     :: oml          ! cut-off value λc for the

    ! quasi-resonant criterion 2)

    !

    integer(kind=8), intent(out)      &

         :: nnz          ! total number of quartets

    ! i.e., nonzero values

    ! in the large-sparse matrix

    ! illustrated above

    !

    ! Local parameters

    real                 :: k(ns)        ! scalar/mag k

    integer              :: i1, i2, i3, i4, row, col

    real                 :: k4x, k4y, k4, om4, kmin, kmax, dom

    !/

    ! Scalar wavenumber (i.e., magnitude)

    k    = sqrt(kx*kx + ky*ky)

    kmin = minval(k)

    kmax = maxval(k)

    !

    ! Boundary conditions: include k4 beyond kmin & kmax

    if (qi_interp .eq. 1 .and. qi_bound .eq. 1) then

      kmin = kmin / 9.  ! 1/3 fmax

      kmax = kmax * 9.  ! 3   fmax

    end if

    !

    ! Start to find the quartets: \vec{k_j}, j = 1, 2, 3 are chosen at the

    ! grid points, and \vec_{k_4} is found by

    !     \vec{k_4} = \vec{k_1} + \vec{k_2} - \vec{k_3}

    !

    nnz = 0

    !

    do i1 = 1, ns

      !           criterion 3.a) ← starting from i1

      do i2 = i1, ns

        do i3 = 1, ns

          !                   criterion 3.b)

          if (i3 .ne. i1 .and. i3 .ne. i2) then

            !                       criterion 1)

            k4x = kx(i1) + kx(i2) - kx(i3)

            k4y = ky(i1) + ky(i2) - ky(i3)

            k4  = sqrt(k4x*k4x + k4y*k4y)

            !

            ! wavenumber k4 falls outside the grid (criterion 3.e)

            if (k4 >= kmin .and. k4 <= kmax) then

              !                           ω = √(qg) & Δω

              om4 = sqrt(qr_grav * qfunc(k4x, k4y))

              dom = abs(om(i1) + om(i2) - om(i3) - om4) / &

                   min(om(i1), om(i2), om(i3), om4)

              !                           criterion 2)

              if (oml <= qr_eps .or. dom <= oml) then

                i4 = minloc((kx - k4x)*(kx - k4x) +     &

                     (ky - k4y)*(ky - k4y), 1)

                !                               criterion 3.d)

                if (i4 .ne. i1 .and. i4 .ne. i2) then

                  !                                   criterion 3.c)

                  if (i4 >= i3) then

                    row = i1 + ns * (i2-1)

                    col = i3 + ns * (i4-1)

                    !                                       criterion 4)

                    if (col > row) then

                      nnz = nnz + 1

                    end if

                  end if

                end if

              end if

            end if

          end if

        end do

      end do

#ifdef W3_TS

      write(*, *) '→ nnz = ', nnz

#endif

    end do

    !/

  end subroutine findquartetnumber

  !/

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

  !/

  subroutine findquartetconfig(ns, kx, ky, om, oml, nnz,     &

       NN, PP, QQ, RR,               &

       k4x, k4y, om4)

    !/

    !/    19-Dec-2011 : Origination.                        ( O. Gramstad )

    !/

    !/    09-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/    02-Apr-2020 : Boundary conditions (< kmin, > kmax)( Q. Liu      )

    !/

    !  1. Purpose:

    !     Find all the quartets that we are interested in. Initially I thought

    !     we may merge this subroutine and the subroutine above (i.e.,

    !     FindQuartetNumber) in such a way that we first initialize a large

    !     array like Quartet(HNum), where HNum is a huge integer (something

    !     like 0.5*NS**4). But it quickly turned out this was a very naive

    !     idea because for the wavenumber grid (k, θ) used by 3G spectral

    !     wave models, in general NS~O(10^2-3), then NS^4~O(10^8-12). Thus,

    !     HNum becomes really very very very huge, and then we may have

    !     the integer/memory overflow problem.

    !

    !     Based on the above-mentioned, we must split the whole process:

    !     1) find the total number of quartets with FindQuartetNumber, `nnz`

    !     2) allocate arrays with the known `nnz`, and store the wavenumber

    !        and ω for k₄

    !

    !     For more details, see the header of the subr. FindQuartetNumber.

    !

    !/

    implicit none

    !

    integer, intent(in)  :: ns           ! length of 1D wavenumber

    ! vector, ns = nk * nth

    real, intent(in)     :: kx(ns),   &

         ky(ns),   &  ! (kx, ky) components

         om(ns)       ! ω or σ

    real, intent(in)     :: oml          ! cut-off value λc for the

    ! quasi-resonant criterion 2)

    integer(kind=8), intent(in)       &

         :: nnz          ! total number of quartets

    ! returned from the subr.

    ! FindQuartetNumber

    !

    integer, intent(out) :: nn(nnz),  &  ! index of k₁

         pp(nnz),  &  !          k₂

         qq(nnz),  &  !          k₃

         rr(nnz)      !          k₄ in the 1D

    ! wavenumber vector [1 - NS]

    real, intent(out)    :: k4x(nnz), &

         k4y(nnz), &  ! x, y comp. of k₄

         om4(nnz)     ! ω₄

    !

    ! Local parameters

    real                 :: k(ns)        ! scalar/mag k

    integer              :: i1, i2, i3, i4, row, col, s

    real                 :: k4xt, k4yt, k4t, om4t, kmin, kmax, dom

    !/

    ! Scalar wavenumber (i.e., magnitude)

    k    = sqrt(kx*kx + ky*ky)

    kmin = minval(k)

    kmax = maxval(k)

    !

    ! Boundary conditions: include k4 beyond kmin & kmax

    if (qi_interp .eq. 1 .and. qi_bound .eq. 1) then

      kmin = kmin / 9.  ! 1/3 fmax

      kmax = kmax * 9.  ! 3   fmax

    end if

    !

    ! Start to find the quartets: \vec{k_j}, j = 1, 2, 3 are chosen at the

    ! grid points, and \vec_{k_4} is found by

    !     \vec{k_4} = \vec{k_1} + \vec{k_2} - \vec{k_3}

    !

    ! s: count of quartets. This time the total number of quartets `nnz` is

    ! already known from `FindQuartetNumber`.

    !       nnz = 0

    s = 0

    !

    do i1 = 1, ns

      !           criterion 3.a) ← starting from i1

      do i2 = i1, ns

        do i3 = 1, ns

          !                   criterion 3.b)

          if (i3 .ne. i1 .and. i3 .ne. i2) then

            !                       criterion 1)

            k4xt = kx(i1) + kx(i2) - kx(i3)

            k4yt = ky(i1) + ky(i2) - ky(i3)

            k4t  = sqrt(k4xt*k4xt + k4yt*k4yt)

            !

            ! wavenumber k4 falls outside the grid (criterion 3.e)

            if (k4t >= kmin .and. k4t <= kmax) then

              !                           ω = √qg & Δω

              om4t = sqrt(qr_grav * qfunc(k4xt, k4yt))

              dom  = abs(om(i1) + om(i2) - om(i3) - om4t)/&

                   min(om(i1), om(i2), om(i3), om4t)

              !                           criterion 2)

              if (oml <= qr_eps .or. dom <= oml) then

                i4 = minloc((kx - k4xt)*(kx - k4xt) +   &

                     (ky - k4yt)*(ky - k4yt), 1)

                !                               criterion 3.d)

                if (i4 .ne. i1 .and. i4 .ne. i2) then

                  !                                   criterion 3.c)

                  if (i4 >= i3) then

                    row = i1 + ns * (i2-1)

                    col = i3 + ns * (i4-1)

                    !                                       criterion 4)

                    if (col > row) then

                      !                                           nnz    = nnz + 1

                      s      = s + 1 ! Find 1 quartet

                      !

                      nn(s)  = i1    ! Store index

                      pp(s)  = i2

                      qq(s)  = i3

                      rr(s)  = i4

                      !

                      k4x(s) = k4xt  ! k₄, ω₄

                      k4y(s) = k4yt

                      om4(s) = om4t

                      !

                    end if

                  end if

                end if

              end if

            end if

          end if

        end do

      end do

    end do

    !

    ! Check consistency of s and nnz

    if (s .ne.  nnz) then

      write(*, 1001) 'FindQuartetConfig'

      call exit(1)

    end if

    !

    ! Formats

1001 FORMAT(/' *** GKE ERROR IN gkeModule : '/  &

         '     Subr. ', a, ': The number of Quartet Configs. does not match NNZ!'/)

    !/

  end subroutine findquartetconfig

  !/

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

  !/ [Part 3]

  !/

  subroutine coocsrind (nrow, nnz, ir, jc, ind_translate, iao)

    !/

    !/    12-Sep-2012 : Origination.                        ( version 3.14 )

    !/                  Based on coocsr of SPARKIT          ( O. Gramstad  )

    !/

    !/    16-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/

    !  1. Purpose:

    !     It becomes clear from subr. FindQuartetNumber & FindQuartetConfig

    !     that we are faced with a problem of large sparse matrice when we

    !     manipulate the huge set of quartets. By sparse matrix we mean

    !     only a `relatively small number` of its matrix elements are nonzero.

    !

    !     For saving time or memory space, a sparse matrix is usually stored

    !     in some compressed formats in the computer memory. Two among those

    !     formats, COO & CSR are relevant here in our application:

    !     1) The coordinate format (COO) --- the simplest storage scheme

    !        For a given sparse matrix `A` (N, N) with NNZ nonzero elements,

    !        the COO format consists of 3 arrays:

    !        * a (nnz): real nonzero values of A in `any order`

    !        * ir(nnz): row indices of these nonzero values

    !        * jc(nnz): column indices

    !

    !     2) The Compressed Sparse Row format (CSR)

    !        The CSR format is the basic format used in SPARSKIT, consisting

    !        of three arrays as well

    !        * a (nnz): real nonzero values of A stored row by row from row

    !                   1 to row N

    !        * jc(nnz): column indices in `any order`

    !        * ia(N+1): the index of the first nonzero element at this

    !                   corresponding row in the array a and jc, that is

    !                   ia(i) provides the position in a & jc where the i-th

    !                   row starts.

    !

    !     This subroutine converts the sparse matrix (nrow, nrow) in the COO

    !     format, as represented by (ir, jc) to the CSR format, as represented

    !     by (ind_translate, iao).

    !

    !     N.B.:

    !     This subr. neither needs the real value array in the COO format,

    !     nor returns the real value array in the CSR format. Alternatively,

    !     it returns the tranformed index (ind_translate) from COO to CSR.

    !     With such indices, we have

    !     *) a_csr  = a_coo(ind_translate)

    !     *) jc_csr = jc_coo(ind_translate)

    !

    !     References:

    !     Youcef Saad, 1994, SPARSKIT: a basic tool kit for sparse matrix

    !         compuation (version 2, `coocsr` therein)

    !     See also Numerical Recipe in Fortran (ch. 2.7, p. 71)

    !/

    implicit none

    !

    integer, intent(in)  :: nrow               ! # of rows of sparse matrix

    integer(kind=8), intent(in)     &

         :: nnz                ! # of nonzero elements

    integer, intent(in)  :: ir(nnz)            ! COO row

    integer, intent(in)  :: jc(nnz)            ! COO col

    integer, intent(out) :: ind_translate(nnz) ! indices from COO to CSR

    integer, intent(out) :: iao(nrow+1)        ! CSR iao

    !

    ! Local parameters

    integer              :: i, j, k, k0, iad

    !/

    ! Determine the number of non-zeros in each row (iao(i), i = 1, ..., nrow,

    ! will be the # of nonzero elements at the i-th row), whereas

    ! iao(nrow+1) = 0

    iao(1:nrow+1) = 0

    do k = 1, nnz

      iao(ir(k)) = iao(ir(k)) + 1 ! row by row

    end do

    !

    ! Find the positions that correspond to the first value in each row.

    ! Now iao(i) is the position where the i-th row starts, and

    ! iao(nrow+1) = 1 + nnz

    k = 1

    do j = 1, nrow+1

      k0     = iao(j)  ! num_i, # of nonzero in this row

      iao(j) = k       ! starting pos

      k      = k + k0  ! k = Σnum_i, where i <= j

    end do

    !

    ! Go through the structure once more. Fill in ind_translate

    do k = 1, nnz

      i = ir(k)        ! coo row

      j = jc(k)        ! coo col

      !

      ! When i-th row is encountered by the first time, iad = iao(i) denotes

      ! the starting position for this row. Afterwards, iao(i) is added by 1

      ! when i-th row arrives every time. In the end, iao(i) records the

      ! starting position for the (i+1)-th row. However, the last element of

      ! iao remains unchanged, i.e., iao(nrow+1) = iao(nrow) = 1 + nnz

      iad                = iao(i)

      ind_translate(iad) = k

      iao(i)             = iad + 1

    end do

    !

    ! Shift back IAO.

    do j = nrow, 1, -1

      iao(j+1) = iao(j)

    end do

    iao(1) = 1

    !

    return

    !/

  end subroutine coocsrind

  !/

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

  !/

  subroutine asymsmattimvec (n, a, ja, ia, x, y, Symb)

    !/

    !/    07-Sep-2012 : Origination.                        ( version 3.14 )

    !/                  Based on amux & atmux of SPARKIT    ( O. Gramstad  )

    !/

    !/    16-Nov-2018 : Prepare WW3 distribution            ( Q. Liu      )

    !/    19-Feb-2018 : Add `Symb` keyword                  ( Q. Liu      )

    !/

    ! 1. Purpose:

    !    --------> Symb = -1 (antisymmetric)

    !    Calculate the dot product of an antisymmetric CSR sparse matrix

    !    and a vector X.

    !

    !    An antisymmetric (skew-symmetric) matrix is a square matrix `B`

    !    whose transpose equals to its negative, i.e.,

    !        B^T = -B

    !

    !    ◆ Do not be confused by the name of this subr. The coming-in CSR

    !      sparse matrix `a` is not symmetric or antisymmetric. In our case,

    !      `a` is a upper triangular sparse matrix, and we are acturally

    !      calculating the dot product of `a - a^T` and `x`, where

    !      'a - a^T' is an antisymmetric matrix due to the symmetry of

    !      four-wave nonlinear interactions (dN₁/dt = -dN₃/dt).

    !

    !    This operation is in essence the dot product of two common dense

    !    matrix/vector, such as

    !        M(n, 1)  = A(n, n) * X(n, 1)

    !        or

    !        M_{i, 1} = Σa(i, j)  * x(j, 1)

    !

    !    For the transposed array A^T,

    !        N_{i, 1} = Σat(i, j) * x(j, 1)

    !                 = Σ a(j, i) * x(j, 1)

    !    Alternatively, we can exchange the index of i, j for easy

    !    understanding:

    !        N_{j, 1} = Σat(j, i) * x(i, 1)

    !                 = Σ a(i, j) * x(i, 1)

    !

    !    Finally, Y = M - N = A * X - A^T * X

    !

    !    --------> Symb = 1 (Symmetric)

    !    Same as above but for Y = M + N = A * X + A^T * X

    !/

    implicit none

    !

    integer, intent(in)   :: n           ! # of rows/cols

    !

    real, intent(in)      :: a(:)        ! CSR a (nnz)

    integer, intent(in)   :: ja(:)       ! CSR ja(nnz)

    integer, intent(in)   :: ia(n+1)     ! CSR ia(n+1)

    !

    real, intent(in)      :: x(n)        ! vector of the same length

    real, intent(out)     :: y(n)        ! return product y = B * x

    real, intent(in)      :: symb        ! -1 for minus, 1 for plus

    !

    ! Local parameters

    integer               :: i, k

    real                  :: t

    !/

    ! Initilization

    y(1:n) = 0.0

    !

    do i = 1, n

      t = 0.0

      do k = ia(i), ia(i+1)-1

        !

        ! M_{i, 1} =  Σa(i, j) * x(j, 1)

        t        = t + a(k) * x(ja(k))

        !

        !±N_{j, 1} = ±Σa(i, j)  * x(i, 1)

        y(ja(k)) = y(ja(k)) + symb * a(k)*x(i)

      end do

      y(i) = y(i) + t

    end do

    ! The final Y = M ± N = A * x ± A^T * x

    !

    return

    !/

  end subroutine asymsmattimvec

  !/

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

  !/ Part 4

  !/

  subroutine prepkgrid(nk, nth, dpt, sig, th)

    !/

    !/    04-Dec-2018 : Origination.                        ( Q. Liu      )

    !/    04-Dec-2018 : Based on `z_cmpcg` & `z_wnumb` of serv_xnl4v5.f90

    !/                                                      ( Q. Liu      )

    !/    01-Apr-2019 : Add the option using WAVNU1         ( Q. Liu      )

    !/

    ! 1. Purpose:

    !    Compute wave number k for a given discrete frequency grid and water

    !    depth based on the dispersion relation for the linear wave theory

    !        ω^2 = gk tanh(kd)

    !

    !    ◆ In WW3, the radian frequency grid ω is invariant.

    !

    !    ◆ It is desired that the GKE module should be independent from WW3

    !      as much as possible. So I decided not to directly obtain

    !      `NK, NTH, SIG, WN, CG, DSII` from WW3

    !

    ! 2. Method

    !    ✓ dispopt = 0

    !    Finite depth linear dispersion relation, using a Pade approximation

    !    (Hunt, 1988) [see WRT serv_xnl4v5.f90]

    !

    !    ✓ dispopt = 1 for WAVNU1

    !/

    USE w3dispmd, ONLY: wavnu1

    !

    implicit none

    !

    integer, intent(in)    :: nk          ! # of frequencies

    integer, intent(in)    :: nth         ! # of directions

    real, intent(in)       :: dpt         ! water depth (m)

    real, intent(in)       :: sig(nk)     ! radian frequency σ

    real, intent(in)       :: th(nth)     ! θ (rad) [equally spaced,

    ! but may start from non-zero

    ! value]

    !

    integer, parameter     :: dispopt = 1 ! dispersion relation

    !

    integer                :: ik, ith, jkth, ns

    real                   :: x, xx, y, omega

    real                   :: k, cg, dsii, angr, dth

    real                   :: esin(nth), ecos(nth)

    !/

    ! Initialization

    ns     = nk * nth

    ! Allocation of qr_kx/ky/dk/om/wn1 was done in PrepKernelIO)

    qr_kx  = 0. ! ns

    qr_ky  = 0.

    qr_dk  = 0.

    qr_om  = 0.

    qr_wn1 = 0. ! nk

    !

    ! Calc Δθ, cosθ, sinθ [θ is equally spaced]

    dth = qr_tpi / real(nth)

    !

    do ith = 1, nth

      angr      = th(ith)

      esin(ith) = sin(angr)

      ecos(ith) = cos(angr)

      !

      if (abs(esin(ith)) .lt. 1.e-5) then

        esin(ith) = 0.

        if (ecos(ith) .gt. 0.5) then

          ecos(ith) = 1.      ! θ = 0.

        else

          ecos(ith) = -1.     ! θ = π

        end if

      end if

      !

      if (abs(ecos(ith)) .lt. 1.e-5) then

        ecos(ith) = 0.

        if (esin(ith) .gt. 0.5) then

          esin(ith) = 1.     ! θ = π/2

        else

          esin(ith) = -1.    ! θ = π * 3/2

        end if

      end if

    end do

    !

    do ik = 1, nk

      if (dispopt .eq. 0) then

        ! Calc k & Cg (`z_cmpcg` & `z_wnumb` of serv_xnl4v5.f90)

        omega   = sig(ik)**2.0/qr_grav

        y       = omega*dpt

        xx      = y*(y+1.0/(1.0+y*(0.66667+y*(0.35550+y*(0.16084+y*(0.06320+y* &

             (0.02174+y*(0.00654+y*(0.00171+y*(0.00039+y*0.00011))))))))))

        x       = sqrt(xx)

        k       = x/dpt

        !

        if(dpt*k > 30.0) then

          cg  = qr_grav/(2.0*sig(ik))

        else

          cg  = sig(ik)/k*(0.5+dpt*k/sinh(2.0*dpt*k))

        end if

        !

      else if (dispopt .eq. 1) then

        ! Calc k & cg (WAVNU1 from WW3)

        call wavnu1(sig(ik), dpt, k, cg)

      end if

      qr_wn1(ik) = k ! Store k in qr_wn1 ('ll used for interp.)

#ifdef W3_TS

      write(*, *) 'σ, k, cg: ', sig(ik), k, cg

#endif

      ! Calc Δσ

      if (ik .eq. 1) then

        dsii = 0.5 * (sig(2) - sig(1))       ! first bin

      else if (ik .eq. nk) then

        dsii = 0.5 * (sig(nk) - sig(nk-1))   ! last bin

      else

        dsii = 0.5 * (sig(ik+1) - sig(ik-1)) ! interm. bin

      end if

      ! Calc Kx, Ky

      do ith = 1, nth

        jkth        = ith + (ik - 1) * nth

        qr_kx(jkth) = k * ecos(ith)

        qr_ky(jkth) = k * esin(ith)

        ! Calc Δ\vec{k} = k Δk Δθ = k Δσ/cg Δθ

        qr_dk(jkth) = k * dsii / cg * dth

        qr_om(jkth) = sig(ik)

      end do

    end do

#ifdef W3_TS

    write(*, *) 'qr_kx: ', qr_kx

    write(*, *) 'qr_ky: ', qr_ky

    write(*, *) 'qr_dk: ', qr_dk

    write(*, *) 'qr_om: ', qr_om

#endif

    !

    return

    !/

  end subroutine prepkgrid

  !/

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

  !/

  subroutine prepkernelio(nk, nth, sig, th, act)

    !/

    !/    04-Dec-2018 : Origination                         ( Q. Liu      )

    !/    04-Dec-2018 : Extracted from Odin's subr. `calcQRSNL`

    !/                                                      ( Q. Liu      )

    !/

    ! 1. Purpose:

    !    Read & Write the pre-computed kernel coefficients `T` for a given

    !    discrete wavenumber grid and water depth.

    !

    !    For a typical 2D finite-depth wave model application, the wavenumber

    !    grid varies according to water depth. Consequently, the quartet

    !    configuration and interactive kernel coefficients will change as

    !    well.

    !

    !    Therefore, it seems extremely difficult to handle a 2D varied-depth

    !    application as the total number of quartets (qi_nnz) and thus the

    !    array size of `Inpqr0` vary [see CalcQRSNL]. Initializing a 2D

    !    array `Inpqr0` with a fixed size of (qi_nnz, nsea) becomes impossible.

    !

    !    So currently we are limiting ourself to deep-water or constant

    !    finite-deep cases.

    !

    !/

    USE constants, ONLY: file_endian


    implicit none

    !

    integer, intent(in)          :: nk             ! # of frequencies

    integer, intent(in)          :: nth            ! # of directions

    real, intent(in)             :: sig(nk)        ! radian frequency (rad)

    real, intent(in)             :: th(nth)        ! θ (rad) [equally spaced,

    ! but may start from non-zero

    ! value]

    character(len=*), intent(in) :: act            ! 'read' or 'write'

    !

    ! Local parameters

    integer                      :: ns, iq, i1, i3, icol

    integer(kind=8)              :: rpos           ! reading position

    integer, allocatable         :: irow_coo(:), & ! row of coo mat

         icootcsr(:)    ! index for coo → csr

    !/

    ! Initilization

    ns      = nk * nth

    qi_nrsm = ns * ns

    ! → Be very careful that the size of `qi_irCsr` is not qi_nnz !

    if (allocated(qi_ircsr)) deallocate(qi_ircsr); allocate(qi_ircsr(qi_nrsm+1))

    if (allocated(qr_sumqr)) deallocate(qr_sumqr); allocate(qr_sumqr(qi_nrsm))

    if (allocated(qr_sumnp)) deallocate(qr_sumnp); allocate(qr_sumnp(ns, ns))

    ! qr_dk/om

    if (allocated(qr_dk))    deallocate(qr_dk);    allocate(qr_dk(ns))

    if (allocated(qr_om))    deallocate(qr_om);    allocate(qr_om(ns))

    !

    ! Determine water depth for the whole module, which will be used by

    ! `T` & `Q` func.

    qr_depth  = max(0., min(qr_depth, qr_dmax))

    qi_disc   = max(0,  min(qi_disc, 1))

    qi_kev    = max(0,  min(qi_kev, 1))

    qi_interp = max(0,  min(qi_interp, 1))

    if (qi_disc .eq. 1) qi_interp = 0

    !

    ! Determine the name for the binary file which stores the quartet

    ! configuration and the corresponding kernel coefficient ['gkev?_d????.cfg]

    ! constant-depth or deep water

    write(qs_cfg, "(A, '_d', I4.4, '.cfg')") trim(qs_ver), int(qr_depth)

    !

    if (trim(act) == 'WRITE') then

      ! Calc KGrid → [qr_kx/ky/dk/om/wn]

      if (allocated(qr_kx))     deallocate(qr_kx);    allocate(qr_kx(ns))

      if (allocated(qr_ky))     deallocate(qr_ky);    allocate(qr_ky(ns))

      if (allocated(qr_wn1))    deallocate(qr_wn1);   allocate(qr_wn1(nk))

      call prepkgrid(nk, nth, qr_depth, sig, th)

      ! Find total # of quartets → [qi_nnz]

      call findquartetnumber(ns, qr_kx, qr_ky, qr_om, qr_oml, qi_nnz)

      ! Find Quartet Config. → [qi_NN/PP/QQ/RR & qr_k4x/k4y/om4]

      if (allocated(qi_nn))     deallocate(qi_nn);    allocate(qi_nn(qi_nnz))

      if (allocated(qi_pp))     deallocate(qi_pp);    allocate(qi_pp(qi_nnz))

      if (allocated(qi_qq))     deallocate(qi_qq);    allocate(qi_qq(qi_nnz))

      if (allocated(qi_rr))     deallocate(qi_rr);    allocate(qi_rr(qi_nnz))

      !

      if (allocated(qr_k4x))    deallocate(qr_k4x);   allocate(qr_k4x(qi_nnz))

      if (allocated(qr_k4y))    deallocate(qr_k4y);   allocate(qr_k4y(qi_nnz))

      if (allocated(qr_om4))    deallocate(qr_om4);   allocate(qr_om4(qi_nnz))

      !

      call findquartetconfig(ns, qr_kx, qr_ky, qr_om, qr_oml, qi_nnz, &

           qi_nn, qi_pp, qi_qq, qi_rr,              &

           qr_k4x, qr_k4y, qr_om4)

      !

      ! Calc Kernel `T`

      if (allocated(qr_tkern))  deallocate(qr_tkern); allocate(qr_tkern(qi_nnz))

      if (allocated(qr_tkurt))  deallocate(qr_tkurt); allocate(qr_tkurt(qi_nnz))

      if (allocated(qr_dom))    deallocate(qr_dom);   allocate(qr_dom(qi_nnz))

      !

      do iq = 1, qi_nnz

        qr_tkern(iq) = tfunc(qr_kx(qi_nn(iq)), qr_ky(qi_nn(iq)),&

             qr_kx(qi_pp(iq)), qr_ky(qi_pp(iq)),&

             qr_kx(qi_qq(iq)), qr_ky(qi_qq(iq)),&

             qr_k4x(iq)      , qr_k4y(iq)      )

      end do

      ! Calc Kernel coeff. for Kurtosis

      qr_tkurt = qr_tkern * sqrt(qr_om(qi_nn) * qr_om(qi_pp) * qr_om(qi_qq) * qr_om4)

      ! Calc Δω (Remove very small Δω; Δω=0 →  resonant quartets)

      qr_dom = qr_om(qi_nn) + qr_om(qi_pp) - qr_om(qi_qq) - qr_om4

      ! TODO: should we use double precision for qr_dom

      ! Note for GNU compiler, qr_eps~1.2E-7 (single prec.) & ~2.2E-16 (double).

      ! The values above are also true for the intel compiler.

      ! sin(Δωt) / Δω is very different for Δω = 0 and Δw~1E-7 when t is large.

      where(abs(qr_dom) < qr_eps) qr_dom = 0.0

      !

      ! Calc interp. weight if necessary

      if (qi_interp .eq. 1) then

        if (allocated(qi_bind)) deallocate(qi_bind); allocate(qi_bind(4, qi_nnz))

        if (allocated(qr_bwgh)) deallocate(qr_bwgh); allocate(qr_bwgh(4, qi_nnz))

        if (qi_bound .eq. 1 ) then

          if (allocated(qr_bdry)) deallocate(qr_bdry); allocate(qr_bdry(qi_nnz))

        end if

        call biinterpwt(nk, nth, qr_wn1, th)

      end if

      !

      deallocate(qr_kx, qr_ky)

      deallocate(qr_k4x, qr_k4y, qr_om4)

      if (qi_interp .eq. 1) deallocate(qr_wn1)

      !

      ! Sparse matrix index conversion [icCos shared by two formats: COO & CSR]

      if (allocated(qi_iccos))  deallocate(qi_iccos); allocate(qi_iccos(qi_nnz))

      if (allocated(irow_coo))  deallocate(irow_coo); allocate(irow_coo(qi_nnz))

      if (allocated(icootcsr))  deallocate(icootcsr); allocate(icootcsr(qi_nnz))

      !

      irow_coo = qi_nn + (qi_pp - 1) * ns

      qi_iccos = qi_qq + (qi_rr - 1) * ns

      !

      ! FindQuartetConfig stores the quartet row by row in a discontinuous order,

      ! so we need keep icooTcsr & qi_irCsr

      call coocsrind(qi_nrsm, qi_nnz, irow_coo, qi_iccos, icootcsr, qi_ircsr)

      !

      ! Reorder index & arrays [coo → crs]

      qi_nn    = qi_nn(icootcsr)     ! used for calc. action prod.

      qi_pp    = qi_pp(icootcsr)

      qi_qq    = qi_qq(icootcsr)

      qi_rr    = qi_rr(icootcsr)

      qr_tkern = qr_tkern(icootcsr)

      qr_tkurt = qr_tkurt(icootcsr)

      qr_dom   = qr_dom(icootcsr)    ! Δω

      !

      if (qi_interp .eq. 1) then     ! bilinear interp. weight

        qi_bind = qi_bind(:, icootcsr)

        qr_bwgh = qr_bwgh(:, icootcsr)

        if (qi_bound .eq. 1) qr_bdry = qr_bdry(icootcsr)

      end if

      !

      qi_iccos = qi_iccos(icootcsr)

      deallocate(irow_coo, icootcsr)

      !

      ! Construct the sum vectors [used for 6D integration]

      ! Σ over Q, R [qr_sumQR]

      qr_sumqr = 2.0

      do i3 = 1, ns

        !              i3 == i4

        icol = i3 + (i3 - 1) * ns

        qr_sumqr(icol) = 1.0

      end do

      ! Σ over P [qr_sumNP]

      qr_sumnp = 1.0

      do i1 = 1, ns

        !               i1 == i2

        qr_sumnp(i1, i1) = 0.5

      end do

      !

      ! WRITE KGrid & Kernel into qs_cfg

      write(*, *) '[W] Writing |', trim(qs_cfg), '| ...'

      open(51, file=trim(qs_cfg), form='unformatted', convert=file_endian, &

           access='stream', status='replace')

      !

      ! It is not necessary to store `ns` since `ns = nk * nth`

      write(51) nk, nth, sig, th                  ! (f, θ) grid

      write(51) qr_depth, qr_oml, qi_disc,        &

           qi_kev, qi_interp                 ! parameters

      write(51) qr_om, qr_dk

      write(51) qi_nnz

      write(51) qi_nn, qi_pp, qi_qq, qi_rr

      write(51) qr_tkern, qr_tkurt, qr_dom

      write(51) qi_iccos, qi_ircsr

      write(51) qr_sumqr, qr_sumnp

      !

      if (qi_interp .eq. 1) write(51) qi_bind, qr_bwgh

      if ( (qi_interp .eq. 1) .and. (qi_bound .eq. 1) ) &

           write(51) qr_bdry

      close(51)

      ! Screen Test

#ifdef W3_TS

      write(*, *) "[W] qr_depth: ", qr_depth

      write(*, *) "[W] qr_oml  : ", qr_oml

      write(*, *) "[W] qi_disc : ", qi_disc

      write(*, *) "[W] qi_kev  : ", qi_kev

      write(*, *) "[W] qr_om   : ", qr_om

      write(*, *) "[W] qr_dk   : ", qr_dk

      write(*, *) "[W] The total number of quartets is ", qi_nnz

      write(*, *) '[W] qi_NN   : ', qi_nn

      write(*, *) '[W] qi_PP   : ', qi_pp

      write(*, *) '[W] qi_QQ   : ', qi_qq

      write(*, *) '[W] qi_RR   : ', qi_rr

      write(*, *) '[W] qr_TKern: ', qr_tkern

      write(*, *) '[W] qr_TKurt: ', qr_tkurt

      write(*, *) '[W] qr_dom  : ', qr_dom

      write(*, *) '[W] qi_icCos: ', qi_iccos

      write(*, *) '[W] qi_irCsr: ', qi_ircsr

      write(*, *) '[W] Σ_QR    : ', qr_sumqr(1: qi_nrsm: ns+1)

      write(*, *) '[W] Σ_P     : ', (qr_sumnp(iq, iq), iq = 1, ns)

#endif

      !

    else if (trim(act) == 'READ') then

      write(*, *) '⊚ → [R] Reading |', trim(qs_cfg), '| ...'

      open(51, file=trim(qs_cfg), form='unformatted', convert=file_endian, &

           access='stream', status='old')

      ! nk, nth, sig, th can be skipped by using pos

      rpos = 1_8 + qi_lrb * (2_8 + nk + nth)

      read(51, pos=rpos) qr_depth, qr_oml, qi_disc,   &

           qi_kev, qi_interp

      !

      ! read ω & Δ\vec{k}

      read(51) qr_om, qr_dk

      ! read total # of quartets

      read(51) qi_nnz

      write(*, *) "⊚ → [R] The total number of quartets is ", qi_nnz

      write(*, *)

      ! allocate arrays

      if (allocated(qi_nn))     deallocate(qi_nn);    allocate(qi_nn(qi_nnz))

      if (allocated(qi_pp))     deallocate(qi_pp);    allocate(qi_pp(qi_nnz))

      if (allocated(qi_qq))     deallocate(qi_qq);    allocate(qi_qq(qi_nnz))

      if (allocated(qi_rr))     deallocate(qi_rr);    allocate(qi_rr(qi_nnz))

      !

      if (allocated(qr_tkern))  deallocate(qr_tkern); allocate(qr_tkern(qi_nnz))

      if (allocated(qr_tkurt))  deallocate(qr_tkurt); allocate(qr_tkurt(qi_nnz))

      if (allocated(qr_dom))    deallocate(qr_dom);   allocate(qr_dom(qi_nnz))

      !

      if (allocated(qi_iccos))  deallocate(qi_iccos); allocate(qi_iccos(qi_nnz))

      !

      read(51) qi_nn, qi_pp, qi_qq, qi_rr

      read(51) qr_tkern, qr_tkurt, qr_dom

      read(51) qi_iccos, qi_ircsr

      read(51) qr_sumqr, qr_sumnp

      !

      if (qi_interp .eq. 1) then

        if (allocated(qi_bind)) deallocate(qi_bind); allocate(qi_bind(4, qi_nnz))

        if (allocated(qr_bwgh)) deallocate(qr_bwgh); allocate(qr_bwgh(4, qi_nnz))

        read(51) qi_bind, qr_bwgh

        !

        if (qi_bound .eq. 1) then

          if (allocated(qr_bdry)) deallocate(qr_bdry); allocate(qr_bdry(qi_nnz))

          read(51) qr_bdry

        end if

        !

      end if

      !

      close(51)

      ! Screen Test

#ifdef W3_TS

      write(*, *) "[R] qr_depth: ", qr_depth

      write(*, *) "[R] qr_oml  : ", qr_oml

      write(*, *) "[R] qi_disc : ", qi_disc

      write(*, *) "[R] qi_kev  : ", qi_kev

      write(*, *) "[R] qr_om   : ", qr_om

      write(*, *) "[R] qr_dk   : ", qr_dk

      write(*, *) "[R] The total number of quartets is ", qi_nnz

      write(*, *) '[R] qi_NN   : ', qi_nn

      write(*, *) '[R] qi_PP   : ', qi_pp

      write(*, *) '[R] qi_QQ   : ', qi_qq

      write(*, *) '[R] qi_RR   : ', qi_rr

      write(*, *) '[R] qr_TKern: ', qr_tkern

      write(*, *) '[R] qr_TKurt: ', qr_tkurt

      write(*, *) '[R] qr_dom  : ', qr_dom

      write(*, *) '[R] qi_icCos: ', qi_iccos

      write(*, *) '[R] qi_irCsr: ', qi_ircsr

      write(*, *) '[R] Σ_QR    : ', qr_sumqr(1: qi_nrsm: ns+1)

      write(*, *) '[R] Σ_P     : ', (qr_sumnp(iq, iq), iq = 1, ns)

#endif

    end if

    !/

  end subroutine prepkernelio

  !/

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

  !/

  subroutine biinterpwt(nk, nth, wn, th)

    !/

    !/    19-Apr-2019 : Origination                         ( Q. Liu      )

    !/    19-Apr-2019 : Extracted from a few subrs. of mod_xnl4v5.f90

    !/                                                      ( Q. Liu      )

    !/    01-Apr-2020 : Boundary conditions                 ( Q. Liu      )

    !/

    ! 1. Purpose:

    !    Calculate weights for the bilinear interpolation.

    !

    ! 2. Method:

    !    See also Fig. 9 of van Vledder (2006, CE) and mod_xnl4v5.f90 (WRT).

    !    [q_t13v4,  q_weight, q_makegrid

    !/

    implicit none

    !

    integer, intent(in) :: nk

    integer, intent(in) :: nth

    real, intent(in)    :: wn(nk)  ! k

    real, intent(in)    :: th(nth) ! θ

    !

    integer             :: iq, jku, jk4, jk4p, jth4t, jth4, jth4p

    real                :: dth, aref, k4t, angr, &

         r_jk, r_jth, delk, w_k4, w_th4

    real                :: kmin, kmax, k4r

    real                :: qr_kpow

    !

    ! Initialization

    qi_bind = 0

    qr_bwgh = 0.

    if (qi_bound .eq. 1) qr_bdry = 1.

    !

    ! Get power law for F(k) from qi_fpow for E(f)

    !     E(f) df = F(k) dk →  F(k) ~ f^n * cg = k^{(n-1)/2}

    !     N(k) = F(k) / ω ~ k^{n/2-1}

    !     C(k) = N(k) / k ~ k^{n/2-1 -1}

    qr_kpow = qr_fpow / 2. - 2.


    ! Kmin & Kmax

    kmin = minval(wn)

    kmax = maxval(wn)

    !

    ! In general, th(nth) in [0, 2π). Note however, it is not the case when

    ! the first directional bin defined in ww3_grid.inp (RTH0) is not zero.

    dth     = qr_tpi / real(nth)

    aref    = th(1)

    !

    ! qr_k4x(nnz), qr_k4y(nnz), wn(nk) are already available for use

    do iq = 1, qi_nnz

      k4r = sqrt(qr_k4x(iq)**2. + qr_k4y(iq)**2.)  ! k₄

      angr= atan2(qr_k4y(iq), qr_k4x(iq))          ! θ₄ [-π, π]

      ! Boundary

      if (qi_bound .eq. 1) then

        k4t = max(kmin, min(k4r, kmax))

      else

        k4t = k4r ! already bounded in [kmin, kmax]

      end if

      !

      ! Layout of surrouding four (f, θ) grid points

      !

      !          (θ)↑

      !             ↑ ₄             ₃

      !       jth4+1 ▪ ----------- ▪

      !              |             |

      !              |      r_jth) |

      !     w-       |---✗ (r_jk,  |

      !     t|       |   |         |

      !     h|       |₁  |         |₂

      !     4-  jth4 ▪ ----------- ▪  → → (k)

      !             jk4            jk4+1

      !              |---|

      !               wk4

      !

      ! i) θ index (counted counterclockwisely)

      r_jth = (angr - aref) / dth + 1.

      jth4t = floor(r_jth)              ! 'll be revised later

      w_th4 = r_jth - real(jth4t)       ! dirc. weight

      !

      jth4  = mod(jth4t-1+nth, nth) + 1 ! wrap around 2π

      jth4p = mod(jth4t+nth, nth) + 1

      !

      ! ii) k index (counted in an ascending order). Note, as required in

      ! FindQuartetConfig, k4T >= kmin & k4T <= kmax are already satisfied.

      ! Thus, the resulted jkU will be in [1, nk].

      ! Two special cases:

      !           /  1,  k4T = kmin

      !     jkU = |

      !           \ NK,  k4T = kmax or k4T in (wn(nk-1), kmax)

      !

      jku = 1

      do while (k4t > wn(jku))

        jku  = jku + 1

        if (jku > nk) exit ! impossible in our case

      end do

      !

      if (jku .eq. 1) then           ! k4T = kmin

        r_jk = 1.

      else                           ! k4T in (kmin, kmax]

        delk = wn(jku) - wn(jku-1) ! Δk

        r_jk = real(jku - 1.) + (k4t - wn(jku-1)) / delk

      end if

      ! Parse r_jk

      jk4   = floor(r_jk)            ! in [1, nk]

      w_k4  = r_jk  - real(jk4)

      jk4p  = min(jk4+1, nk)         ! k4T = kmax ← min func.

      !

      ! Store indices (in 1D vector; jkth = ith + (ik-1) * nth)

      qi_bind(1, iq) = jth4  + (jk4  - 1) * nth

      qi_bind(2, iq) = jth4  + (jk4p - 1) * nth

      qi_bind(3, iq) = jth4p + (jk4p - 1) * nth

      qi_bind(4, iq) = jth4p + (jk4  - 1) * nth

      !

      ! Store weights

      qr_bwgh(1, iq) = (1. - w_k4) * (1. - w_th4)

      qr_bwgh(2, iq) = w_k4 * (1. - w_th4)

      qr_bwgh(3, iq) = w_k4 * w_th4

      qr_bwgh(4, iq) = (1. - w_k4) * w_th4

      !

      ! Note that the qi_bind & qr_bwgh do not make full sense when

      ! k4 < kmin (k indices are not correct at all) or k4 > kmax (k index = NK)

      ! because we have capped k4T in between kmin and kmax.

      ! But no need to worry about this because

      !     1) C(k) = 0. when k < kmin

      !     2) C(k) = C(NK) * power decay when k > kmax

      !

      if (qi_bound .eq. 1) then

        if (k4r < kmin) then

          qr_bdry(iq) = 0.

        else if (k4r > kmax) then

          qr_bdry(iq) = (k4r/kmax)**qr_kpow

        end if

      end if

      !

    end do

    !/

  end subroutine biinterpwt

  !/

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

  !/

  subroutine calcqrsnl(nk, nth, sig, th,         &

       t0, t1, Cvk0, Cvk1,       &

       Inpqr0, Snl, Dnl, Kurt)

    !/

    !/    09-Dec-2018 : Origination                         ( Q. Liu      )

    !/    09-Dec-2018 : Extracted from Odin's subr. `calcQRSNL`

    !/                                                      ( Q. Liu      )

    !/    10-Jun-2019 : Include Janssen's KE properly       ( Q. Liu      )

    !/    07-Jun-2021 : Switch off the cal. of kurtosis (!|KT|)

    !/                                                      ( Q. Liu      )

    !/

    ! 1. Purpose:

    !    Calculate the nonlinear transfer rates for a given frequency

    !    grid and given action density spectrum C(\vec{k}).

    !

    !    According to J09 and GS13, C(\vec{k}) is given by

    !             /

    !        m0 = | F(\vec{k}) d \vec{k}

    !             /

    !

    !        F(\vec{k}) = ω C(\vec{k}) / g,

    !

    !    whereas the wave action density spectrum used in WW3 is given by

    !        F(\vec{k}) d \vec{k} = F(k, θ)   dk dθ

    !                             = N(k, θ) ω dk dθ

    !

    !    Thus, we have

    !        C(\vec{k}) = N * g / k.

    !

    ! 2. Method

    !    See GS13 & GB16 for all the details.

    !

    !    ◆ t0, t1 here are time `relative` to the begining time of the

    !      simulaiton, rather than the `absolute` time

    !

    !    ◆ Cvk0, Cvk1, Snl, Dnl shoud be organized/stored in the same way as

    !      qr_kx, qr_ky, qr_dk, qr_om

    !/

    implicit none

    !

    integer, intent(in)        :: nk             ! # of frequencies

    integer, intent(in)        :: nth            ! # of directions

    real, intent(in)           :: sig(nk)        ! radian frequency (rad)

    real, intent(in)           :: th(nth)        ! θ (rad) [equally spaced,

    ! but may start from non-zero

    !

    real, intent(inout)        :: t0             ! previous time step

    real, intent(inout)        :: cvk0(nk*nth)   ! Action density @ t0

    complex, intent(inout)     :: inpqr0(:)      ! I(t) @ t0

    !

    real, intent(in)           :: t1             ! current  time step

    real, intent(in)           :: cvk1(nk*nth)   ! Action density @ t1

    ! ... C(\vec{k})

    !

    real, intent(out)          :: snl(nk*nth)    ! Snl = dC/dt

    real, intent(out)          :: dnl(nk*nth)    ! Dnl

    real, intent(out)          :: kurt           ! Kurtosis

    !

    ! Local parameters

    !

    real                       :: delt           ! Δt

    logical, save              :: flread = .true.

    integer                    :: num_i, ns

    real                       :: dvk0(nk*nth),& ! Odin's discrete Cvk @ t0

         dvk1(nk*nth)   ! ...     @ t1

    real, allocatable, save    :: cvk0_r(:),   & ! C₄      @ t0

         cvk1_r(:)      !         @ t1

    real, allocatable, save    :: fnpqr0(:),   & ! C prod. @ t0

         fnpqr1(:),   & ! C prod. @ t1

         mnpqr (:)      ! δC_1/δt * δk_1

    complex, allocatable, save :: inpqr1(:),   & ! I(t) @ t1

         etau(:),     & ! exp(iΔωt)

         edelt(:)       ! exp(iΔωΔt)

    !

    real, allocatable, save    :: mnp1d(:), mnp2d(:, :)

    real                       :: secm2          ! Second-order moment²

    !/

    !

    ! Initilization

    ns      = nk * nth

    qi_nrsm = ns * ns

    !

    ! Only constant depth is allowed now. Accordingly, we only need a single

    ! binary config file which provides the wavenumber grid and kernel

    ! coefficients.

    !

    ! Read quartets & kernel coefficients in

    if (flread) then

      ! Only read data once

      call prepkernelio(nk, nth, sig, th, 'READ')

      flread = .false.

      !           write(*, *) "⊚ → [R] FLag for Reading Kernels becomes |", FlRead, "|"

      ! Allocate arrays

      ! ✓ A variable with the SAVE attribute retains its value and definition,

      ! association, and `allocation` status on exit from a procedure

      !

      if (allocated(fnpqr0)) deallocate(fnpqr0); allocate(fnpqr0(qi_nnz))

      if (allocated(fnpqr1)) deallocate(fnpqr1); allocate(fnpqr1(qi_nnz))

      if (allocated(etau  )) deallocate(etau  ); allocate(etau(qi_nnz))

      if (allocated(edelt )) deallocate(edelt ); allocate(edelt(qi_nnz))

      if (allocated(mnpqr )) deallocate(mnpqr ); allocate(mnpqr(qi_nnz))

      if (allocated(inpqr1)) deallocate(inpqr1); allocate(inpqr1(qi_nnz))

      !

      if (allocated(mnp1d))  deallocate(mnp1d);  allocate(mnp1d(qi_nrsm))

      if (allocated(mnp2d))  deallocate(mnp2d);  allocate(mnp2d(ns, ns))

      !

      if (qi_disc .eq. 0) then

        if (allocated(cvk0_r)) deallocate(cvk0_r); allocate(cvk0_r(qi_nnz))

        if (allocated(cvk1_r)) deallocate(cvk1_r); allocate(cvk1_r(qi_nnz))

      end if

      !

    end if

    !

    ! Screen output (check whether the kernel data are stored in memory)

#ifdef W3_TS

    write(*, *) "◆ qr_depth      :", qr_depth

    write(*, *) "◆ qr_oml        :", qr_oml

    write(*, *) "◆ qi_disc       :", qi_disc

    write(*, *) "◆ qi_kev        :", qi_kev

    write(*, *) "◆ qi_nnz        :", qi_nnz

    write(*, *) "◆ qi_NN(:10)    :", qi_nn(:10)

    write(*, *) "◆ qr_TKern(:10) :", qr_tkern(:10)

    write(*, *) "◆ qr_TKurt(:10) :", qr_tkurt(:10)

    write(*, *) "◆ qr_sumQR(:10) :", qr_sumqr(:10)

#endif

    !

    num_i = size(inpqr0)

    if (num_i .ne. qi_nnz) then

      write(*, 1001) 'CalcQRSNL'

      call exit(1)

    end if

    !

    ! Start to calc. Snl term

    if (qi_disc == 0) then

      ! Define ΔC = dC/dt * Δk Δt, we have ΔC₁ = ΔC₂ = -ΔC₃ = -ΔC₄ (Δt can be

      ! removed by taking the unit time)

      !

      ! Cvk0/1_R (bilinear interp. or nearest bin)

      if (qi_interp .eq. 0) then

        cvk0_r = cvk0(qi_rr)

        cvk1_r = cvk1(qi_rr)

        !

      else if (qi_interp .eq. 1) then

        cvk0_r = qr_bwgh(1, :) * cvk0(qi_bind(1, :)) + &

             qr_bwgh(2, :) * cvk0(qi_bind(2, :)) + &

             qr_bwgh(3, :) * cvk0(qi_bind(3, :)) + &

             qr_bwgh(4, :) * cvk0(qi_bind(4, :))

        !

        cvk1_r = qr_bwgh(1, :) * cvk1(qi_bind(1, :)) + &

             qr_bwgh(2, :) * cvk1(qi_bind(2, :)) + &

             qr_bwgh(3, :) * cvk1(qi_bind(3, :)) + &

             qr_bwgh(4, :) * cvk1(qi_bind(4, :))

        !

        if (qi_bound .eq. 1) then

          cvk0_r = cvk0_r * qr_bdry

          cvk1_r = cvk1_r * qr_bdry

        end if

        !

      end if

      !

      ! F = [C₃ C₄ (C₁ + C₂) - C₁ C₂ (C₃ + C₄)] dk₂ dk₃ dk₄ ∙ dk₁

      ! dk₄ vanishes with the δ function

      fnpqr0 = (cvk0(qi_qq) * cvk0_r      * (  &

           cvk0(qi_nn) + cvk0(qi_pp) ) -  &

           cvk0(qi_nn) * cvk0(qi_pp) * (  &

           cvk0(qi_qq) + cvk0_r      )) * &

           qr_dk(qi_nn) * qr_dk(qi_pp) * qr_dk(qi_qq)

      !

      fnpqr1 = (cvk1(qi_qq) * cvk1_r      * (  &

           cvk1(qi_nn) + cvk1(qi_pp) ) -  &

           cvk1(qi_nn) * cvk1(qi_pp) * (  &

           cvk1(qi_qq) + cvk1_r      )) * &

           qr_dk(qi_nn) * qr_dk(qi_pp) * qr_dk(qi_qq)

      !

    else if (qi_disc == 1) then

      ! Used in GS13 & GB16

      ! F = [C₃dk₃ C₄dk₄ (C₁dk₁ + C₂dk₂) - C₁dk₁ C₂dk₂ (C₃dk₃ + C₄dk₄)]

      ! It seems the bilinear interpolation for this discretization approach

      ! is not very meaningful.

      dvk0   = cvk0 * qr_dk

      fnpqr0 = dvk0(qi_qq) * dvk0(qi_rr) * ( &

           dvk0(qi_nn) + dvk0(qi_pp) ) - &

           dvk0(qi_nn) * dvk0(qi_pp) * ( &

           dvk0(qi_qq) + dvk0(qi_rr) )

      !

      dvk1   = cvk1 * qr_dk

      fnpqr1 = dvk1(qi_qq) * dvk1(qi_rr) * ( &

           dvk1(qi_nn) + dvk1(qi_pp) ) - &

           dvk1(qi_nn) * dvk1(qi_pp) * ( &

           dvk1(qi_qq) + dvk1(qi_rr) )

      !

    end if

    !|KT|! Calc m2 for Kurtosis estimation ((2.6) of Annekov & Shrira (2013))

    !|KT|        SecM2 = sum(Cvk1 * qr_om * qr_dk) ** 2.

    !

    !       write(*, *) '.... Input args: t0, t1 :', t0, t1

    if (abs(t1) < qr_eps) then

      ! t1 = 0.0 [essentially I₁ = 0 → I₀ = 0]

      t0     = 0.0

      cvk0   = cvk1

      inpqr0 = (0.0, 0.0)  ! \int_{0}^{0} dt  = 0

      snl    = 0.0

      dnl    = 0.0

      kurt   = 0.0

    else

      ! t1 ≠ 0.0

      delt   = t1 - t0

      if (delt < 0.0) then

        write(*, 1002) 'CalcQRSNL'

        call exit(2)

      end if

      etau   = exp(qc_iu * cmplx(qr_dom * t1))       ! exp(iΔωt)

      edelt  = exp(qc_iu * cmplx(qr_dom * delt))     ! exp(iΔωΔt)

      !

      ! ◆ Calc. I₁: note here I₁ = I(t₁) dk₁ dk₂ dk₃ for both qi_disc = 0/1

      if (qi_kev .eq. 0) then

        ! GKE from GS13, GB16

        inpqr1 = inpqr0 + cmplx(0.5 * delt) *  &

             conjg(etau)                *  &       ! exp(-iΔωt)

             (cmplx(fnpqr0) * edelt + cmplx(fnpqr1))

      else if (qi_kev .eq. 1) then

        ! KE from J03 (Fnpqr1 is taken outside the time integral; Fnpqr0 is not

        ! used in this case; and the real part of Inpqr1 is sin(Δωt)/Δω, and

        ! the imaginary part is [1 - cos(Δωt)] / Δω

        ! Approximation used before

        !               Inpqr1 = Inpqr0 + cmplx(0.5 * DelT) *  &

        !                        conjg(ETau) * (EDelT + 1)

        !

        where (abs(qr_dom) < qr_eps)

          ! Δω = 0., sin(Δωt)/Δω ~ t, [1 - cos(Δωt)] / Δω ~ 0

          inpqr1 = cmplx(t1, 0.)

        elsewhere

          ! Δω ≠ 0., cacl. sin(Δωt)/Δω & [1 - cos(Δωt)] / Δω directly

          ! TODO: the sign of cos is not clear yet.

          inpqr1 = cmplx(sin(qr_dom * t1) / qr_dom,     &

               (1 - cos(qr_dom * t1)) / qr_dom)

        end where

      end if

      ! ◆ Snl [Tranfer Integal]

      if (qi_kev .eq. 0) then

        ! GKE from GS13, GB16

        mnpqr  = 4.0 * (qr_tkern ** 2.) * real(etau * inpqr1)

      else if (qi_kev .eq. 1) then

        ! KE from J03

        !               Mnpqr  = 4.0 * (qr_TKern ** 2.) * Fnpqr1 * real(ETau * Inpqr1)

        mnpqr  = 4.0 * (qr_tkern ** 2.) * fnpqr1 * real(inpqr1)

      end if

      ! Calc. Σ over Q, R [Mnpqr is a upper triangular sparse matrix]

      ! dN₁/dt = - dN₃/dt →  anti-symmetric array operation

      ! Mnp1D  = (Mnpqr - Mnpqr^{T}) × S_{qr}

      call asymsmattimvec(qi_nrsm, mnpqr, qi_iccos, qi_ircsr, qr_sumqr, mnp1d, -1.0)

      ! Calc. Σ over P [Mnp2D is a upper triangular matrix]

      ! dN₁/dt = dN₂/dt →  symmetric array operation

      ! Snl    = {Σ (Mnp + Mnp^{T}) ⊙ S_{p}} / d\vec{k₁}

      mnp2d  = reshape(mnp1d, (/ns, ns/))

      snl    = sum((mnp2d + transpose(mnp2d)) * qr_sumnp, 2) / qr_dk

      ! ◆ Conservation Check

#ifdef W3_TS

      write(*, '(A, E15.3)') '   ← {WW3 GKE } ΣSnl(k) * dk: ', sum(snl * qr_dk)

#endif

      !

      ! ◆ Dnl [Diagonal term]  <TODO>

      !   i) it is easy to calculate Dnl for Janssen's KE (but we may

      !      have to abandon the sparse array approach)

      !  ii) it is challenging to get Dnl for GKE.

      dnl  = 0.0

      kurt = 0.0

      !

      !|KT|! ◆ Kurtosis

      !|KT|            if (qi_kev .eq. 0) then

      !|KT|! GKE from GS13, GB16

      !|KT|                Mnpqr  = -3.0 / SecM2 * qr_TKurt * aimag(ETau * Inpqr1)

      !|KT|            else if (qi_kev .eq. 1) then

      !|KT|! KE from J03 (here the imaginary part becomes [1 - cos(Δωt)] / Δω

      !|KT|!               Mnpqr  = -3.0 / SecM2 * qr_TKurt * Fnpqr1 * aimag(ETau * Inpqr1)

      !|KT|                Mnpqr  = -3.0 / SecM2 * qr_TKurt * Fnpqr1 * aimag(Inpqr1)

      !|KT|            end if

      !|KT|! Calc. Σ over Q, R [Mnpqr is a upper triangular sparse matrix]

      !|KT|! symmetric array operation Mnp1D  = (Mnpqr - Mnpqr^{T}) × S_{qr}

      !|KT|            call ASymSmatTimVec(qi_nrsm, Mnpqr, qi_icCos, qi_irCsr, qr_sumQR, Mnp1D, 1.0)

      !|KT|            Mnp2D  = reshape(Mnp1D, (/ns, ns/))

      !|KT|            Kurt   = sum((Mnp2D + transpose(Mnp2D)) * qr_sumNP)

      !

      ! I₁ → I₀ for next computation (time step)

      t0     = t1

      cvk0   = cvk1

      inpqr0 = inpqr1

    end if

    !

    !       write(*, *) '.... Output args: t0, t1 :', t0, t1

    !

    ! Formats

1001 FORMAT(/' *** GKE ERROR IN gkeModule : '/  &

         '     Subr. ', a, ': the stored total number of quartets &

         & and the size of Inpqr0 do not match !'/)

1002 FORMAT(/' *** GKE ERROR IN gkeModule : '/  &

         '     Subr. ', a, ': t0 ≤ t1 is not satisfied !'/)

    !/

  end subroutine calcqrsnl

  !/

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

end module w3gkemd

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