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