tests_zgrid_metrics Function

  logical function tests_zgrid_metrics()
    use modglobal, only : kb, ke, kh, zf, zh, dzf, dzh, dzfi, dzhi, dzfiq, &
                          dzfi5, dzhiq, dzf2, dzh2i, dzfc, dzfci, dzhci, &
                          delta, dxf, dy, ipoiss, POISS_FFT2D, POISS_FFT2D_2DECOMP, BCzp
    implicit none

    real :: tol, stretch
    integer :: k, i, nfail
    logical :: stretched

    nfail = 0
    tol = 1000.*epsilon(zh(kb))*max(1., zh(ke+1))

    call expect(abs(zh(kb)) <= tol, 'zh(kb) is 0')

    do k = kb, ke
      call expect(dzf(k) > 0., 'dzf positive')
      call expect(abs(dzf(k) - (zh(k+1) - zh(k))) <= tol, 'dzf = zh(k+1)-zh(k)')
      ! the midpoint identity: 0.5*(w(k)+w(k+1)) style averages depend on it
      call expect(abs(zf(k) - 0.5*(zh(k) + zh(k+1))) <= tol, 'zf is the midpoint of its faces')
      if (k > kb) call expect(zf(k) > zf(k-1), 'zf strictly increasing')
    end do

    call expect(abs(dzh(kb) - dzf(kb)) <= tol, 'dzh(kb) ghost convention')
    do k = kb + 1, ke
      call expect(abs(dzh(k) - (zf(k) - zf(k-1))) <= tol, 'dzh = zf(k)-zf(k-1)')
      call expect(abs(dzh(k) - 0.5*(dzf(k-1) + dzf(k))) <= tol, 'dzh is the mean of adjacent dzf')
    end do
    call expect(abs(dzf(kb-1) - dzf(kb)) <= tol, 'dzf bottom ghost convention')
    call expect(abs(dzf(ke+1) - dzf(ke)) <= tol, 'dzf top ghost convention')

    do k = kb, ke + kh
      call expect(abs(dzfi(k)*dzf(k) - 1.) <= tol, 'dzfi = 1/dzf')
      call expect(abs(dzhi(k)*dzh(k) - 1.) <= tol, 'dzhi = 1/dzh')
      call expect(abs(dzfiq(k) - 0.25*dzfi(k)) <= tol, 'dzfiq = dzfi/4')
      call expect(abs(dzhiq(k) - 0.25*dzhi(k)) <= tol, 'dzhiq = dzhi/4')
      call expect(abs(dzfi5(k) - 0.5*dzfi(k)) <= tol, 'dzfi5 = dzfi/2')
      call expect(abs(dzf2(k) - dzf(k)*dzf(k)) <= tol*dzf(k), 'dzf2 = dzf**2')
      call expect(abs(dzh2i(k) - dzhi(k)*dzhi(k)) <= tol*dzhi(k)*dzhi(k), 'dzh2i = dzhi**2')
    end do

    ! kappa-scheme copies exist only when a kappa advection scheme is selected
    if (allocated(dzfc)) then
      do k = kb, ke
        call expect(abs(dzfc(k) - dzf(k)) <= tol, 'dzfc matches dzf on the interior')
        call expect(abs(dzfci(k)*dzfc(k) - 1.) <= tol, 'dzfci = 1/dzfc')
        call expect(abs(dzhci(k) - dzhi(k)) <= tol*dzhi(k), 'dzhci matches dzhi on the interior')
      end do
    end if

    do k = kb, ke
      i = lbound(delta, 1)
      call expect(abs(delta(i,k) - (dxf(i)*dy*dzf(k))**(1./3.)) <= tol, &
                  'delta = (dxf*dy*dzf)**(1/3)')
    end do

    ! a stretched grid must only be combined with a Poisson solve that treats the
    ! vertical direction by Gaussian elimination; the FFT branches hard-code
    ! dzfi(1) in their vertical eigenvalues (modpois: "Assumes equidistant")
    stretch = maxval(dzf(kb:ke))/minval(dzf(kb:ke))
    stretched = stretch > 1. + 100.*epsilon(stretch)
    if (stretched) then
      call expect((ipoiss == POISS_FFT2D .and. BCzp == 1) .or. ipoiss == POISS_FFT2D_2DECOMP, &
                  'stretched grid requires ipoiss=0 with BCzp=1 or ipoiss=3 (vertical Gaussian elimination)')
    end if

    if (myid == 0) then
      write(*,'(A,ES10.2,A,L1)') 'zgrid metrics: max/min dzf = ', stretch, ', stretched = ', stretched
      if (nfail == 0) then
        write(*,*) 'PASS: all vertical grid metric checks passed'
      else
        write(*,'(A,I0,A)') 'FAIL: ', nfail, ' vertical grid metric checks failed'
      end if
    end if

    tests_zgrid_metrics = (nfail == 0)

  contains

    subroutine expect(condition, label)
      logical, intent(in) :: condition
      character(*), intent(in) :: label
      if (.not. condition) then
        nfail = nfail + 1
        if (myid == 0) write(*,'(A,A,A,I0)') 'FAIL zgrid metric: ', label, ' at k=', k
      end if
    end subroutine expect

  end function tests_zgrid_metrics