w3gkemd Module Reference

Functions/Subroutines
subroutine, public	prepkernelio (nk, nth, sig, th, act)
subroutine, public	calcqrsnl (nk, nth, sig, th, t0, t1, Cvk0, Cvk1, Inpqr0, Snl, Dnl, Kurt)

Variables
real, public	qr_depth
real, public	qr_oml
integer, public	qi_disc
integer, public	qi_kev
integer, public	qi_interp
integer(kind=8), public	qi_nnz

Function/Subroutine Documentation

calcqrsnl()

subroutine, public w3gkemd::calcqrsnl	(	integer, intent(in)	nk,
		integer, intent(in)	nth,
		real, dimension(nk), intent(in)	sig,
		real, dimension(nth), intent(in)	th,
		real, intent(inout)	t0,
		real, intent(in)	t1,
		real, dimension(nk*nth), intent(inout)	Cvk0,
		real, dimension(nk*nth), intent(in)	Cvk1,
		complex, dimension(:), intent(inout)	Inpqr0,
		real, dimension(nk*nth), intent(out)	Snl,
		real, dimension(nk*nth), intent(out)	Dnl,
		real, intent(out)	Kurt
	)

Definition at line 1625 of file w3gkemd.F90.

    !/
    !/    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 !'/)
    !/

References prepkernelio(), qi_disc, qi_interp, qi_kev, qi_nnz, qr_depth, and qr_oml.

Referenced by w3snl5md::w3snl5().

prepkernelio()

subroutine, public w3gkemd::prepkernelio	(	integer, intent(in)	nk,
		integer, intent(in)	nth,
		real, dimension(nk), intent(in)	sig,
		real, dimension(nth), intent(in)	th,
		character(len=*), intent(in)	act
	)

Definition at line 1202 of file w3gkemd.F90.

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

References file(), constants::file_endian, qi_disc, qi_interp, qi_kev, qi_nnz, qr_depth, and qr_oml.

Referenced by calcqrsnl(), and w3snl5md::insnl5().

Variable Documentation

qi_disc

integer, public w3gkemd::qi_disc

Definition at line 172 of file w3gkemd.F90.

  integer                    :: qi_disc                     ! Discretization of GKE

Referenced by calcqrsnl(), w3snl5md::insnl5(), and prepkernelio().

qi_interp

integer, public w3gkemd::qi_interp

Definition at line 180 of file w3gkemd.F90.

  integer                    :: qi_interp                   ! Interp. option

Referenced by calcqrsnl(), w3snl5md::insnl5(), and prepkernelio().

qi_kev

integer, public w3gkemd::qi_kev

Definition at line 176 of file w3gkemd.F90.

  integer                    :: qi_kev                      ! Version of KE

Referenced by calcqrsnl(), w3snl5md::insnl5(), and prepkernelio().

qi_nnz

integer(kind=8), public w3gkemd::qi_nnz

Definition at line 209 of file w3gkemd.F90.

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

Referenced by calcqrsnl(), w3snl5md::insnl5(), and prepkernelio().

qr_depth

real, public w3gkemd::qr_depth

Definition at line 165 of file w3gkemd.F90.

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

Referenced by calcqrsnl(), w3snl5md::insnl5(), prepkernelio(), and w3snl5md::w3snl5().

qr_oml

real, public w3gkemd::qr_oml

Definition at line 167 of file w3gkemd.F90.

  real                       :: qr_oml                      ! λ cut off factor

Referenced by calcqrsnl(), w3snl5md::insnl5(), and prepkernelio().