tests_sgs_statistics Function

  logical function tests_sgs_statistics()
    use modglobal,    only : ib, ie, ih, jb, je, jh, kb, ke, kh, runmode, &
                             dx, dzf, dzfi, dzh, dzhi, zf, zh
    use modfields,    only : initfields, um, vm, wm, thlm, u0, v0, w0, thl0, up, wp, thlp
    use modsubgrid,   only : initsubgrid, diffc, diffu, diffw, ekh, ekm
    use modboundary,  only : initboundary
    use modibm,       only : initibm, createmasks
    use modstatsdump, only : compute_sgs_fluxes
    use initfac,      only : readfacetfiles
    use decomp_2d,    only : zstart

    implicit none

    real, parameter :: acoef = 0.37   ! u   = acoef*z^2
    real, parameter :: bcoef = 0.21   ! w   = bcoef*x + dcoef*z
    real, parameter :: ccoef = 0.53   ! thl = ccoef*z^2
    real, parameter :: dcoef = 0.43
    real, parameter :: nuconst = 0.17 ! constant ekm for the manufactured solution
    real, parameter :: kaconst = 0.29 ! constant ekh for the manufactured solution

    real, allocatable :: usgs(:,:,:), vsgs(:,:,:), wsgs(:,:,:)
    real, allocatable :: thlsgs(:,:,:), qtsgs(:,:,:)
    real, allocatable :: sv1sgs(:,:,:), sv2sgs(:,:,:), sv3sgs(:,:,:), sv4sgs(:,:,:)
    real, allocatable :: zc(:), zw(:), got(:,:,:), exp(:,:,:)
    logical :: all_passed
    real    :: trunc_tol, scale
    integer :: i, j, k, gi

    if (myid == 0) then
      write(*, '(A)') '================================================'
      write(*, '(A, I8)') 'runmode = ', runmode
      write(*, '(A)') 'tests_sgs_statistics: SGS STATISTICS TEST'
      write(*, '(A)') '------------------------------------------------'
      write(*, '(A)') 'Using stretched-grid case 300 to validate modstatsdump SGS fluxes'
    end if

    call initfields
    call initboundary
    call initsubgrid
    call readfacetfiles
    call initibm
    call createmasks

    allocate(usgs  (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))
    allocate(vsgs  (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))
    allocate(wsgs  (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))
    allocate(thlsgs(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))
    allocate(qtsgs (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))
    allocate(sv1sgs(ib:ie,jb:je,kb:ke+kh))
    allocate(sv2sgs(ib:ie,jb:je,kb:ke+kh))
    allocate(sv3sgs(ib:ie,jb:je,kb:ke+kh))
    allocate(sv4sgs(ib:ie,jb:je,kb:ke+kh))
    allocate(got(ib:ie,jb:je,kb:ke+kh), exp(ib:ie,jb:je,kb:ke+kh))

    ! Cell-centre heights, extended one level below the surface by mirroring
    ! (zh(kb) = 0), which is the level the kb stencils actually reach into.
    allocate(zc(kb-kh:ke+kh), zw(kb-kh:ke+kh))
    do k = kb, ke+kh
      zc(k) = zf(k)
      zw(k) = zh(k)
    end do
    zc(kb-kh) = zf(kb) - dzh(kb)
    zw(kb-kh) = zh(kb) - dzf(kb)

    all_passed = .true.

    !----------------------------------------------------------------------
    ! (a) Manufactured solution on the stretched grid
    !----------------------------------------------------------------------
    ekm = nuconst
    ekh = kaconst

    do k = kb-kh, ke+kh
      do j = jb-jh, je+jh
        do i = ib-ih, ie+ih
          gi = zstart(1) + i - 1
          um(i,j,k)   = acoef * zc(k)**2
          vm(i,j,k)   = 0.
          wm(i,j,k)   = bcoef * (real(gi) - 0.5) * dx + dcoef * zw(k)
          thlm(i,j,k) = ccoef * zc(k)**2
        end do
      end do
    end do

    call compute_sgs_fluxes(usgs,vsgs,wsgs,thlsgs,qtsgs,sv1sgs,sv2sgs,sv3sgs,sv4sgs)

    ! (a1) exact discrete identity: the second-order stencils are exact for a
    !      quadratic in z and a linear in x, so no truncation error enters here.
    do k = kb, ke+kh
      got(:,:,k) = usgs(ib:ie,jb:je,k)
      exp(:,:,k) = -nuconst * (acoef*(zc(k) + zc(k-1)) + bcoef)
    end do
    if (.not. compare_real_3d('usgs discrete', got, exp, 1.e-10)) all_passed = .false.

    do k = kb, ke+kh
      got(:,:,k) = thlsgs(ib:ie,jb:je,k)
      exp(:,:,k) = -kaconst * ccoef * (zc(k) + zc(k-1))
    end do
    if (.not. compare_real_3d('thlsgs discrete', got, exp, 1.e-10)) all_passed = .false.

    ! wsgs is the normal stress tau_33 = -2*K_m*dw/dz at the CELL CENTRE, so with
    ! dw/dz = dcoef it is exactly -2*nu*dcoef everywhere. If wsgs were the
    ! vertical part of the diffw tendency instead (units m/s^2, not m^2/s^2) this
    ! would be identically zero for a linear w and the check would fail.
    do k = kb, ke
      got(:,:,k) = wsgs(ib:ie,jb:je,k)
      exp(:,:,k) = -2. * nuconst * dcoef
    end do
    got(:,:,ke+kh) = 0.; exp(:,:,ke+kh) = 0.
    if (.not. compare_real_3d('wsgs discrete', got, exp, 1.e-10)) all_passed = .false.

    ! (a2) second-order consistency with the analytic flux at the w-levels.
    !      The stretched grid gives zf(k)+zf(k-1) - 2*zh(k) = (dzf(k)-dzf(k-1))/2,
    !      so anything worse than that bound (a dzf/dzh swap, a missing factor,
    !      a sign error) is caught here.
    trunc_tol = 0.
    do k = kb+1, ke+kh
      trunc_tol = max(trunc_tol, abs(dzf(k) - dzf(k-1)))
    end do
    trunc_tol = 0.5 * trunc_tol

    do k = kb, ke+kh
      got(:,:,k) = usgs(ib:ie,jb:je,k)
      exp(:,:,k) = -nuconst * (2.*acoef*zh(k) + bcoef)
    end do
    if (.not. compare_real_3d('usgs vs analytic', got, exp, &
                              1.05*nuconst*acoef*trunc_tol + 1.e-10)) all_passed = .false.

    do k = kb, ke+kh
      got(:,:,k) = thlsgs(ib:ie,jb:je,k)
      exp(:,:,k) = -kaconst * 2.*ccoef * zh(k)
    end do
    if (.not. compare_real_3d('thlsgs vs analytic', got, exp, &
                              1.05*kaconst*ccoef*trunc_tol + 1.e-10)) all_passed = .false.

    !----------------------------------------------------------------------
    ! (c) Sign regression: monotone positive shear must give a negative
    !     (down-gradient) SGS momentum flux. One-line guard on the convention.
    !----------------------------------------------------------------------
    if (maxval(usgs(ib:ie,jb:je,kb+1:ke)) >= 0.) then
      if (myid == 0) write(*,'(A,1X,A,1X,ES12.4)') 'FAIL', 'usgs sign', &
           maxval(usgs(ib:ie,jb:je,kb+1:ke))
      all_passed = .false.
    end if
    if (maxval(thlsgs(ib:ie,jb:je,kb+1:ke)) >= 0.) then
      if (myid == 0) write(*,'(A,1X,A,1X,ES12.4)') 'FAIL', 'thlsgs sign', &
           maxval(thlsgs(ib:ie,jb:je,kb+1:ke))
      all_passed = .false.
    end if

    !----------------------------------------------------------------------
    ! (b) Consistency with the solver's diffusion terms.
    !     u, v and thl vary in z only and w varies in x only, so the x- and
    !     y-terms of diffu/diffc vanish identically and what is left is exactly
    !     minus the vertical divergence of the SGS fluxes.
    !----------------------------------------------------------------------
    do k = kb-kh, ke+kh
      do j = jb-jh, je+jh
        do i = ib-ih, ie+ih
          gi = zstart(1) + i - 1
          ekm(i,j,k)  = 0.5 + 0.4*sin(1.7*real(gi) + 0.9*real(j) + 0.31*real(k))
          ekh(i,j,k)  = 0.6 + 0.3*cos(0.7*real(gi) - 1.1*real(j) + 0.53*real(k))
          um(i,j,k)   = sin(0.41*real(k)) + 0.3*cos(0.17*real(k))
          vm(i,j,k)   = cos(0.23*real(k)) - 0.2*sin(0.61*real(k))
          wm(i,j,k)   = 0.7*sin(0.29*(real(gi) - 0.5)*dx)
          thlm(i,j,k) = 288. + sin(0.37*real(k)) + 0.4*cos(0.13*real(k))
        end do
      end do
    end do
    u0   = um
    v0   = vm
    w0   = wm
    thl0 = thlm

    call compute_sgs_fluxes(usgs,vsgs,wsgs,thlsgs,qtsgs,sv1sgs,sv2sgs,sv3sgs,sv4sgs)

    up   = 0.
    thlp = 0.
    call diffu(up)
    call diffc(ih,jh,kh,thl0,thlp)

    do k = kb, ke
      got(:,:,k) = up(ib:ie,jb:je,k)
      exp(:,:,k) = -(usgs(ib:ie,jb:je,k+1) - usgs(ib:ie,jb:je,k)) * dzfi(k)
    end do
    got(:,:,ke+kh) = 0.; exp(:,:,ke+kh) = 0.
    scale = max(maxval(abs(exp)), 1.)
    if (.not. compare_real_3d('diffu vs div(usgs)', got, exp, 1.e-10*scale)) all_passed = .false.

    do k = kb, ke
      got(:,:,k) = thlp(ib:ie,jb:je,k)
      exp(:,:,k) = -(thlsgs(ib:ie,jb:je,k+1) - thlsgs(ib:ie,jb:je,k)) * dzfi(k)
    end do
    got(:,:,ke+kh) = 0.; exp(:,:,ke+kh) = 0.
    scale = max(maxval(abs(exp)), 1.)
    if (.not. compare_real_3d('diffc vs div(thlsgs)', got, exp, 1.e-10*scale)) all_passed = .false.

    ! Same idea for the normal stress: with u = v = 0 and w a function of z only,
    ! the horizontal terms of diffw vanish and what is left is exactly minus the
    ! vertical divergence of wsgs. This is the check that pins wsgs to a flux
    ! (m^2/s^2) rather than a tendency (m/s^2).
    do k = kb-kh, ke+kh
      do j = jb-jh, je+jh
        do i = ib-ih, ie+ih
          um(i,j,k) = 0.
          vm(i,j,k) = 0.
          wm(i,j,k) = sin(0.33*real(k)) + 0.5*cos(0.19*real(k))
        end do
      end do
    end do
    u0 = um
    v0 = vm
    w0 = wm

    call compute_sgs_fluxes(usgs,vsgs,wsgs,thlsgs,qtsgs,sv1sgs,sv2sgs,sv3sgs,sv4sgs)

    wp = 0.
    call diffw(wp)

    got = 0.; exp = 0.
    do k = kb+1, ke
      got(:,:,k) = wp(ib:ie,jb:je,k)
      exp(:,:,k) = -(wsgs(ib:ie,jb:je,k) - wsgs(ib:ie,jb:je,k-1)) * dzhi(k)
    end do
    scale = max(maxval(abs(exp)), 1.)
    if (.not. compare_real_3d('diffw vs div(wsgs)', got, exp, 1.e-10*scale)) all_passed = .false.

    deallocate(usgs, vsgs, wsgs, thlsgs, qtsgs, sv1sgs, sv2sgs, sv3sgs, sv4sgs, got, exp, zc, zw)

    if (all_passed .and. myid == 0) then
      write(*, '(A)') '------------------------------------------------'
      write(*, '(A)') 'ALL TESTS PASSED: tests_sgs_statistics'
      write(*, '(A)') '================================================'
    else if ((.not. all_passed) .and. myid == 0) then
      write(*, '(A)') '------------------------------------------------'
      write(*, '(A)') 'TESTS FAILED: tests_sgs_statistics'
      write(*, '(A)') '================================================'
    end if

    tests_sgs_statistics = all_passed

  end function tests_sgs_statistics