meke_lateral Subroutine

   pure subroutine meke_lateral(nx, ny, sdt, kh_bg, k4, khmeke_fac, &
                                kh_flux_enabled, dy_cu, dx_cv, idxCu, idyCv, &
                                iareaT, i_mass, mass, kh_diff, &
                                uflux, vflux, del2, meke)
      !! Harmonic-mass Laplacian diffusion of MEKE (+ optional biharmonic).
      !! Flux-form, conservative on a closed domain (interior faces only;
      !! array-edge faces carry zero flux).
      !!   Kh_u = max(0,kh_bg) + khmeke_fac*0.5*(kh_i+kh_{i+1}), CFL-capped 0.25
      !!   uflux = Kh_u*(dy_cu*idxCu)*[2 m_i m_{i+1}/(m_i+m_{i+1}+eps)]*(E_i-E_{i+1})
      !!   E += sdt*iareaT*i_mass*((uflux_{i-1}-uflux_i)+(vflux_{j-1}-vflux_j))
      !! Biharmonic: del2 = iareaT*(d uflux' + d vflux') with the bare-gradient
      !! flux uflux' = (dy_cu*idxCu)*(E_{i+1}-E_i); then a harmonic-mass flux
      !! of del2 with CFL cap 0.3 and E += that divergence (additive).
      integer, intent(in) :: nx, ny
      real(wp), intent(in) :: sdt, kh_bg, k4, khmeke_fac
      logical, intent(in) :: kh_flux_enabled
      real(wp), intent(in) :: dy_cu(nx + 1, ny)
      real(wp), intent(in) :: dx_cv(nx, ny + 1)
      real(wp), intent(in) :: idxCu(nx + 1, ny)
      real(wp), intent(in) :: idyCv(nx, ny + 1)
      real(wp), intent(in) :: iareaT(nx, ny)
      real(wp), intent(in) :: i_mass(nx, ny)
      real(wp), intent(in) :: mass(nx, ny)
      real(wp), intent(in) :: kh_diff(nx, ny)
      real(wp), intent(inout) :: uflux(nx + 1, ny)
      real(wp), intent(inout) :: vflux(nx, ny + 1)
      real(wp), intent(inout) :: del2(nx, ny)
      real(wp), intent(inout) :: meke(nx, ny)
      integer :: i, j
      real(wp) :: kh_u, kh_v, hm, geo, inv_max, k4_u, k4_v, mke

      ! ---------- Biharmonic (computed first; tendency added after diffusion). ----------
      if (k4 >= 0.0_wp) then
         ! bare-gradient flux into uflux/vflux workspaces (units m^2/s^2).
         do concurrent(j=1:ny, i=1:nx + 1)
            uflux(i, j) = 0.0_wp
         end do
         do concurrent(j=1:ny, i=2:nx)
            uflux(i, j) = (dy_cu(i, j)*idxCu(i, j))*(meke(i, j) - meke(i - 1, j))
         end do
         do concurrent(j=1:ny + 1, i=1:nx)
            vflux(i, j) = 0.0_wp
         end do
         do concurrent(j=2:ny, i=1:nx)
            vflux(i, j) = (dx_cv(i, j)*idyCv(i, j))*(meke(i, j) - meke(i, j - 1))
         end do
         do concurrent(j=1:ny, i=1:nx)
            del2(i, j) = iareaT(i, j)*((uflux(i + 1, j) - uflux(i, j)) + &
                                       (vflux(i, j + 1) - vflux(i, j)))
         end do
         ! harmonic-mass flux of del2 with K4 (CFL cap 0.3).
         do concurrent(j=1:ny, i=1:nx + 1)
            uflux(i, j) = 0.0_wp
         end do
         do concurrent(j=1:ny, i=2:nx) local(k4_u, geo, hm, inv_max)
            geo = dy_cu(i, j)*idxCu(i, j)
            inv_max = 64.0_wp*sdt*(geo*max(iareaT(i - 1, j), iareaT(i, j)))**2
            k4_u = k4
            if (k4_u*inv_max > 0.3_wp) k4_u = 0.3_wp/inv_max
            hm = 2.0_wp*mass(i - 1, j)*mass(i, j)/((mass(i - 1, j) + mass(i, j)) + MASS_NEGLECT)
            uflux(i, j) = (k4_u*geo)*hm*(del2(i, j) - del2(i - 1, j))
         end do
         do concurrent(j=1:ny + 1, i=1:nx)
            vflux(i, j) = 0.0_wp
         end do
         do concurrent(j=2:ny, i=1:nx) local(k4_v, geo, hm, inv_max)
            geo = dx_cv(i, j)*idyCv(i, j)
            inv_max = 64.0_wp*sdt*(geo*max(iareaT(i, j - 1), iareaT(i, j)))**2
            k4_v = k4
            if (k4_v*inv_max > 0.3_wp) k4_v = 0.3_wp/inv_max
            hm = 2.0_wp*mass(i, j - 1)*mass(i, j)/((mass(i, j - 1) + mass(i, j)) + MASS_NEGLECT)
            vflux(i, j) = (k4_v*geo)*hm*(del2(i, j) - del2(i, j - 1))
         end do
         do concurrent(j=1:ny, i=1:nx) local(mke)
            mke = sdt*(iareaT(i, j)*i_mass(i, j))* &
                  ((uflux(i, j) - uflux(i + 1, j)) + (vflux(i, j) - vflux(i, j + 1)))
            del2(i, j) = mke   ! stash the biharmonic tendency in del2
         end do
      end if

      ! ---------- Laplacian (harmonic-mass) diffusion. ----------
      if (kh_flux_enabled) then
         do concurrent(j=1:ny, i=1:nx + 1)
            uflux(i, j) = 0.0_wp
         end do
         do concurrent(j=1:ny, i=2:nx) local(kh_u, geo, hm, inv_max)
            geo = dy_cu(i, j)*idxCu(i, j)
            kh_u = max(0.0_wp, kh_bg) + khmeke_fac*0.5_wp*(kh_diff(i - 1, j) + kh_diff(i, j))
            inv_max = 2.0_wp*sdt*(geo*max(iareaT(i - 1, j), iareaT(i, j)))
            if (kh_u*inv_max > 0.25_wp) kh_u = 0.25_wp/inv_max
            hm = 2.0_wp*mass(i - 1, j)*mass(i, j)/((mass(i - 1, j) + mass(i, j)) + MASS_NEGLECT)
            uflux(i, j) = (kh_u*geo)*hm*(meke(i - 1, j) - meke(i, j))
         end do
         do concurrent(j=1:ny + 1, i=1:nx)
            vflux(i, j) = 0.0_wp
         end do
         do concurrent(j=2:ny, i=1:nx) local(kh_v, geo, hm, inv_max)
            geo = dx_cv(i, j)*idyCv(i, j)
            kh_v = max(0.0_wp, kh_bg) + khmeke_fac*0.5_wp*(kh_diff(i, j - 1) + kh_diff(i, j))
            inv_max = 2.0_wp*sdt*(geo*max(iareaT(i, j - 1), iareaT(i, j)))
            if (kh_v*inv_max > 0.25_wp) kh_v = 0.25_wp/inv_max
            hm = 2.0_wp*mass(i, j - 1)*mass(i, j)/((mass(i, j - 1) + mass(i, j)) + MASS_NEGLECT)
            vflux(i, j) = (kh_v*geo)*hm*(meke(i, j - 1) - meke(i, j))
         end do
         do concurrent(j=1:ny, i=1:nx) local(mke)
            mke = sdt*(iareaT(i, j)*i_mass(i, j))* &
                  ((uflux(i, j) - uflux(i + 1, j)) + (vflux(i, j) - vflux(i, j + 1)))
            meke(i, j) = meke(i, j) + mke
         end do
      end if

      ! add the biharmonic tendency (computed above, stashed in del2).
      if (k4 >= 0.0_wp) then
         do concurrent(j=1:ny, i=1:nx)
            meke(i, j) = meke(i, j) + del2(i, j)
         end do
      end if
   end subroutine meke_lateral