ERF_PBLModels.H

Go to the documentation of this file.

#ifndef ERF_PBLMODELS_H_
#define ERF_PBLMODELS_H_
 
#include <functional>
#include <limits>
#include <memory>
 
#include <AMReX_MultiFab.H>
#include <AMReX_Periodicity.H>
#include "ERF_Constants.H"
#include "ERF_DataStruct.H"
#include "ERF_MOSTUtils.H"
 
/**
 * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
 * Mellor-Yamada-Janjić (MYJ) boundary layer scheme.
 * Reference: Janjić, Z. I., 1994: The Step-Mountain Eta Coordinate Model:
 * Further Developments of the Convection, Viscous Sublayer, and Turbulence
 * Closure Schemes. Monthly Weather Review, 122, 927–945.
 * https://doi.org/10.1175/1520-0493(1994)122<0927:TSMECM>2.0.CO;2
 *
 * @param[in] xvel Velocity in x-dir
 * @param[in] yvel Velocity in y-dir
 * @param[in] cons_in Cell center conserved quantities
 * @param[out] eddyViscosity Holds turbulent viscosity
 * @param[in] geom Problem geometry
 * @param[in] turbChoice Container with turbulence parameters
 * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
 * @param[in] use_moisture If we have microphysics enabled
 * @param[in] level Current level
 * @param[in] bc_ptr Pointer to array with boundary condition info
 * @param[in] vert_only Only compute vertical eddy diffusivities
 * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
 * @param[in] z_phys_cc Physical location of cell centers
 */
void
ComputeDiffusivityMYJ (double dt,
                       const amrex::MultiFab& xvel,
                       const amrex::MultiFab& yvel,
                       amrex::MultiFab& cons_in,
                       amrex::MultiFab& eddyViscosity,
                       const amrex::Geometry& geom,
                       const TurbChoice& turbChoice,
                       std::unique_ptr<SurfaceLayer>& SurfLayer,
                       bool use_terrain_fitted_coords,
                       bool use_moisture,
                       int level,
                       const amrex::BCRec* bc_ptr,
                       bool /*vert_only*/,
                       const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                       const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                       const MoistureComponentIndices& moisture_indices);
 
/**
 * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
 * Mellor-Yamada-Nakanishi-Niino Level amrex::Real(2.5) PBL scheme
 *
 * @param[in] xvel Velocity in x-dir
 * @param[in] yvel Velocity in y-dir
 * @param[in] cons_in Cell center conserved quantities
 * @param[out] eddyViscosity Holds turbulent viscosity
 * @param[in] geom Problem geometry
 * @param[in] turbChoice Container with turbulence parameters
 * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
 * @param[in] use_moisture If we have microphysics enabled
 * @param[in] level Current level
 * @param[in] bc_ptr Pointer to array with boundary condition info
 * @param[in] vert_only Only compute vertical eddy diffusivities
 * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
 * @param[in] z_phys_cc Physical location of cell centers
 */
void
ComputeDiffusivityMYNN25 (const amrex::MultiFab& xvel,
                          const amrex::MultiFab& yvel,
                          const amrex::MultiFab& cons_in,
                          amrex::MultiFab& eddyViscosity,
                          const amrex::Geometry& geom,
                          const TurbChoice& turbChoice,
                          std::unique_ptr<SurfaceLayer>& SurfLayer,
                          bool use_terrain_fitted_coords,
                          bool use_moisture,
                          int level,
                          const amrex::BCRec* bc_ptr,
                          bool /*vert_only*/,
                          const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                          const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                          const MoistureComponentIndices& moisture_indices);
 
/**
 * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
 * Mellor-Yamada-Nakanishi-Niino Level amrex::Real(2.5) PBL scheme
 *
 * @param[in] xvel Velocity in x-dir
 * @param[in] yvel Velocity in y-dir
 * @param[in] cons_in Cell center conserved quantities
 * @param[out] eddyViscosity Holds turbulent viscosity
 * @param[in] geom Problem geometry
 * @param[in] turbChoice Container with turbulence parameters
 * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
 * @param[in] use_moisture If we have microphysics enabled
 * @param[in] level Current level
 * @param[in] bc_ptr Pointer to array with boundary condition info
 * @param[in] vert_only Only compute vertical eddy diffusivities
 * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
 * @param[in] z_phys_cc Physical location of cell centers
 */
void
ComputeDiffusivityMYNNEDMF (const amrex::MultiFab& xvel,
                            const amrex::MultiFab& yvel,
                            const amrex::MultiFab& cons_in,
                            amrex::MultiFab& eddyViscosity,
                            const amrex::Geometry& geom,
                            const TurbChoice& turbChoice,
                            std::unique_ptr<SurfaceLayer>& SurfLayer,
                            bool use_terrain_fitted_coords,
                            bool use_moisture,
                            int level,
                            const amrex::BCRec* bc_ptr,
                            bool /*vert_only*/,
                            const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                            const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                            const MoistureComponentIndices& moisture_indices);
 
/**
 * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
 * Yonsei University PBL scheme
 *
 * @param[in] xvel Velocity in x-dir
 * @param[in] yvel Velocity in y-dir
 * @param[in] cons_in Cell center conserved quantities
 * @param[out] eddyViscosity Holds turbulent viscosity
 * @param[in] geom Problem geometry
 * @param[in] turbChoice Container with turbulence parameters
 * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
 * @param[in] use_moisture If we have microphysics enabled
 * @param[in] level Current level
 * @param[in] bc_ptr Pointer to array with boundary condition info
 * @param[in] vert_only Only compute vertical eddy diffusivities
 * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
 * @param[in] z_phys_cc Physical location of cell centers
 */
void
ComputeDiffusivityYSU (const amrex::MultiFab& xvel,
                       const amrex::MultiFab& yvel,
                       const amrex::MultiFab& cons_in,
                       amrex::MultiFab& eddyViscosity,
                       const amrex::Geometry& geom,
                       const TurbChoice& turbChoice,
                       std::unique_ptr<SurfaceLayer>& SurfLayer,
                       bool use_terrain_fitted_coords,
                       bool use_moisture,
                       int level,
                       const amrex::BCRec* bc_ptr,
                       bool /*vert_only*/,
                       const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                       const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                       const MoistureComponentIndices& moisture_indices);
/**
 * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
 * Medium Range Forecast (MRF) boundary layer scheme based on Hong and Pan (1996).
 * Reference: Hong, S. Y., and H.-L. Pan, 1996: Nonlocal Boundary Layer Vertical
 * Diffusion in a Medium-Range Forecast Model. Monthly Weather Review, 124,
 * 2322-2339. https://doi.org/10.1175/1520-0493(1996)124<2322:NBLVDI>2.0.CO;2
 *
 * @param[in] xvel Velocity in x-dir
 * @param[in] yvel Velocity in y-dir
 * @param[in] cons_in Cell center conserved quantities
 * @param[out] eddyViscosity Holds turbulent viscosity
 * @param[in] geom Problem geometry
 * @param[in] turbChoice Container with turbulence parameters
 * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
 * @param[in] use_moisture If we have microphysics enabled
 * @param[in] level Current level
 * @param[in] bc_ptr Pointer to array with boundary condition info
 * @param[in] vert_only Only compute vertical eddy diffusivities
 * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
 * @param[in] z_phys_cc Physical location of cell centers
 */
void
ComputeDiffusivityMRF (const amrex::MultiFab& xvel,
                       const amrex::MultiFab& yvel,
                       const amrex::MultiFab& cons_in,
                       amrex::MultiFab& eddyViscosity,
                       const amrex::Geometry& geom,
                       const TurbChoice& turbChoice,
                       std::unique_ptr<SurfaceLayer>& SurfLayer,
                       bool use_terrain_fitted_coords,
                       bool use_moisture,
                       int level,
                       const amrex::BCRec* bc_ptr,
                       bool /*vert_only*/,
                       const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                       const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                       const MoistureComponentIndices& moisture_indices,
                       const amrex::MultiFab* terrain_blank = nullptr);
 
/**
 * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
 * WRF Yonsei University (YSU) PBL scheme.
 *
 * Based on:
 *   Hong, Noh &amp; Dudhia (2006), MWR, https://doi.org/10.1175/MWR3250.1
 *   Hong (2010), QJRMS, https://doi.org/10.1002/qj.665
 *
 * Features: three-pass bulk-Richardson PBLH diagnosis, nonlocal
 * countergradient fluxes, explicit entrainment at PBL top, and
 * grid-adaptive Richardson-number mixing in free atmosphere.
 *
 * @param[in]  xvel   x-direction velocity
 * @param[in]  yvel   y-direction velocity
 * @param[in]  cons_in  cell-center conserved quantities
 * @param[out] eddyViscosity  turbulent eddy viscosity
 * @param[in]  geom   problem geometry
 * @param[in]  turbChoice  turbulence parameter container
 * @param[in]  SurfLayer  Monin-Obukhov surface layer object
 * @param[in]  use_terrain_fitted_coords  flag for terrain coordinates
 * @param[in]  use_moisture  flag for moisture active
 * @param[in]  level  AMR level
 * @param[in]  bc_ptr  boundary condition array
 * @param[in]  vert_only  compute vertical diffusivity only
 * @param[in]  z_phys_nd  physical node heights (terrain)
 * @param[in]  z_phys_cc  physical cell-center heights
 * @param[in]  moisture_indices  component index mapping for moisture variables
 * @param[in]  qheating_rates  radiation heating rates (SW, LW components); optional
 */
void
ComputeDiffusivityYSUNew (const amrex::MultiFab& xvel,
                          const amrex::MultiFab& yvel,
                          const amrex::MultiFab& cons_in,
                          amrex::MultiFab& eddyViscosity,
                          const amrex::Geometry& geom,
                          const TurbChoice& turbChoice,
                          std::unique_ptr<SurfaceLayer>& SurfLayer,
                          bool use_terrain_fitted_coords,
                          bool use_moisture,
                          int level,
                          const amrex::BCRec* bc_ptr,
                          bool /*vert_only*/,
                          const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                          const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                          const MoistureComponentIndices& moisture_indices,
                          const amrex::MultiFab* qheating_rates = nullptr,
                          const amrex::MultiFab* terrain_blank = nullptr);
 
/**
 * Function for computing vertical derivatives for use in PBL model
 *
 * @param[in] u velocity in x-dir
 * @param[in] v velocity in y-dir
 * @param[in] cell_data conserved cell center vars
 */
struct PBLDerivativeDzInv
{
    // Each value is the inverse physical distance between the two samples
    // used by its corresponding derivative stencil.
    amrex::Real interior;
    amrex::Real dirichlet_lo;
    amrex::Real dirichlet_hi;
};
 
/**
 * Compute vertical derivatives of potential temperature and velocity for PBL models.
 *
 * @param i x-index
 * @param j y-index
 * @param k z-index
 * @param[in] uvel x-velocity
 * @param[in] vvel y-velocity
 * @param[in] cell_data conserved cell-center variables
 * @param[in] izmin lower z-index limit
 * @param[in] izmax upper z-index limit
 * @param[in] dz_inv inverse physical distance spacings
 * @param[in] c_ext_dir_on_zlo Dirichlet condition on z-low face for cell data
 * @param[in] c_ext_dir_on_zhi Dirichlet condition on z-high face for cell data
 * @param[in] u_ext_dir_on_zlo Dirichlet condition on z-low face for u-velocity
 * @param[in] u_ext_dir_on_zhi Dirichlet condition on z-high face for u-velocity
 * @param[in] v_ext_dir_on_zlo Dirichlet condition on z-low face for v-velocity
 * @param[in] v_ext_dir_on_zhi Dirichlet condition on z-high face for v-velocity
 * @param[out] dthetadz computed vertical derivative of potential temperature
 * @param[out] dudz computed vertical derivative of x-velocity
 * @param[out] dvdz computed vertical derivative of y-velocity
 * @param[in] moisture_indices moisture variable component indices
 */
/**
 * Thermodynamic variable whose vertical gradient the PBL schemes use for
 * buoyancy, selected by erf.thermal_stratification.
 *
 * @param[in] strat stratification variable to return
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
GetPBLStratVar (int i, int j, int k,
                const StratType strat,
                const amrex::Array4<const amrex::Real>& cell_data,
                const MoistureComponentIndices& moisture_indices)
{
    if (strat == StratType::theta) {
        return cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
    } else if (strat == StratType::thetal) {
        return GetThetal(i,j,k,cell_data,moisture_indices);
    }
    return GetThetav(i,j,k,cell_data,moisture_indices);
}
 
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
ComputeVerticalDerivativesPBL (int i, int j, int k,
                               const amrex::Array4<const amrex::Real>& uvel,
                               const amrex::Array4<const amrex::Real>& vvel,
                               const amrex::Array4<const amrex::Real>& cell_data,
                               const int izmin,
                               const int izmax,
                               const PBLDerivativeDzInv& dz_inv,
                               const bool c_ext_dir_on_zlo,
                               const bool c_ext_dir_on_zhi,
                               const bool u_ext_dir_on_zlo,
                               const bool u_ext_dir_on_zhi,
                               const bool v_ext_dir_on_zlo,
                               const bool v_ext_dir_on_zhi,
                               amrex::Real& dthetadz,
                               amrex::Real& dudz,
                               amrex::Real& dvdz,
                               const MoistureComponentIndices& moisture_indices,
                               const StratType strat = StratType::thetav)
{
    // ext_dir ghost values hold the value at the physical boundary, which is
    // half a ghost-center spacing from the adjacent cell center. Use a
    // two-point, first-order boundary-to-center difference with inverse
    // distance 2/dz.
    if ( k==izmax && c_ext_dir_on_zhi ) {
        dthetadz = ( GetPBLStratVar(i,j,k+1,strat,cell_data,moisture_indices)
                    -GetPBLStratVar(i,j,k  ,strat,cell_data,moisture_indices) )*dz_inv.dirichlet_hi;
    } else if ( k==izmin && c_ext_dir_on_zlo ) {
        dthetadz = ( GetPBLStratVar(i,j,k  ,strat,cell_data,moisture_indices)
                    -GetPBLStratVar(i,j,k-1,strat,cell_data,moisture_indices) )*dz_inv.dirichlet_lo;
    } else {
        dthetadz = ( GetPBLStratVar(i,j,k+1,strat,cell_data,moisture_indices)
                    -GetPBLStratVar(i,j,k-1,strat,cell_data,moisture_indices) )*dz_inv.interior;
    }
 
    if ( k==izmax && u_ext_dir_on_zhi ) {
        dudz = myhalf*( uvel(i  ,j,k+1) - uvel(i  ,j,k)
                       +uvel(i+1,j,k+1) - uvel(i+1,j,k) )*dz_inv.dirichlet_hi;
    } else if ( k==izmin && u_ext_dir_on_zlo ) {
        dudz = myhalf*( uvel(i  ,j,k) - uvel(i  ,j,k-1)
                       +uvel(i+1,j,k) - uvel(i+1,j,k-1) )*dz_inv.dirichlet_lo;
    } else {
        dudz = myhalf*( uvel(i  ,j,k+1) - uvel(i  ,j,k-1)
                       +uvel(i+1,j,k+1) - uvel(i+1,j,k-1) )*dz_inv.interior;
    }
 
    if ( k==izmax && v_ext_dir_on_zhi ) {
        dvdz = myhalf*( vvel(i,j  ,k+1) - vvel(i,j  ,k)
                       +vvel(i,j+1,k+1) - vvel(i,j+1,k) )*dz_inv.dirichlet_hi;
    } else if ( k==izmin && v_ext_dir_on_zlo ) {
        dvdz = myhalf*( vvel(i,j  ,k) - vvel(i,j  ,k-1)
                       +vvel(i,j+1,k) - vvel(i,j+1,k-1) )*dz_inv.dirichlet_lo;
    } else {
        dvdz = myhalf*( vvel(i,j  ,k+1) - vvel(i,j  ,k-1)
                       +vvel(i,j+1,k+1) - vvel(i,j+1,k-1) )*dz_inv.interior;
    }
}
 
/**
 * Monin-Obukhov quantities needed to rebuild the vertical gradients in the first
 * cell above a surface-layer boundary.
 */
struct PBLSurfaceLayerGradient
{
    amrex::Real u_star  {zero};  ///< friction velocity
    amrex::Real tstar_v {zero};  ///< virtual temperature scale
    amrex::Real zval    {zero};  ///< height of the cell center above the surface
    amrex::Real zeta    {zero};  ///< stability coordinate z/L
};
 
/**
 * Virtual temperature scale, thetav_* = theta_*(1 + epsv qv) + epsv theta q_*.
 * The surface virtual heat flux is <w'thetav'> = -u_* thetav_*.
 *
 * @param[in] t_star temperature scale
 * @param[in] q_star moisture scale
 * @param[in] theta potential temperature
 * @param[in] qv water vapor mixing ratio
 * @param[in] use_moisture whether the moisture term is active
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
ComputeVirtualTStarPBL (amrex::Real t_star,
                        amrex::Real q_star,
                        amrex::Real theta,
                        amrex::Real qv,
                        bool use_moisture)
{
    amrex::Real tstar_v = t_star * (one + epsv*qv);
    if (use_moisture) { tstar_v += epsv * theta * q_star; }
    return tstar_v;
}
 
/**
 * Overwrite the vertical gradients in the first cell above the surface with the
 * MOST profile gradients,
 *   |dU/dz|    = u_*      phi_m(zeta) / (kappa z)
 *   dthetav/dz = thetav_* phi_h(zeta) / (kappa z)
 *
 * With zlo.type = surface_layer the state carries foextrap and the velocities
 * hoextrap at the lower boundary, so the centered stencil in
 * ComputeVerticalDerivativesPBL reads ghost values that hold no surface
 * information there: the potential-temperature gradient comes out half the
 * one-sided value and the velocity gradient has no dependence on the friction
 * velocity at all (ERF #4037). The MOST profile used here is the one the surface
 * stress and heat flux were themselves derived from.
 *
 * Every PBL scheme that consumes these gradients uses only the shear magnitude
 * dudz^2 + dvdz^2, so that magnitude is carried in dudz and dvdz is left at zero
 * rather than being split by wind direction.
 *
 * @param[in] sl Monin-Obukhov inputs for this column
 * @param[out] dthetadz vertical derivative of virtual potential temperature
 * @param[out] dudz vertical derivative of x-velocity (carries the shear magnitude)
 * @param[out] dvdz vertical derivative of y-velocity (set to zero)
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
ApplySurfaceLayerGradientsPBL (const PBLSurfaceLayerGradient& sl,
                               amrex::Real& dthetadz,
                               amrex::Real& dudz,
                               amrex::Real& dvdz)
{
    const similarity_funs sfuns{};
    const amrex::Real inv_kappa_z = one / (KAPPA * sl.zval);
    dudz     = sl.u_star  * sfuns.calc_phi_m(sl.zeta) * inv_kappa_z;
    dvdz     = zero;
    dthetadz = sl.tstar_v * sfuns.calc_phi_h(sl.zeta) * inv_kappa_z;
}
 
/**
 * Function for computing the QKE source terms (NN09, Eqn. 5).
 *
 * @param[in] u velocity in x-dir
 * @param[in] v velocity in y-dir
 * @param[in] cell_data conserved cell center vars
 * @param[in] cell_prim primitive cell center vars
 * @param[in] K_turb turbulent viscosity
 * @param[in] domain box of the whole domain
 * @param[in] dz_inv inverse spacings for interior and Dirichlet-boundary differences
 * @param[in] pbl_mynn_B1_l a parameter
 * @param[in] theta_mean average theta
 * @param[in] moisture_indices moisture variable component indices
 * @param[in] tau13 rho<u'w'> on z-faces; the surface-layer stress at the zlo face
 * @param[in] tau23 rho<v'w'> on z-faces; the surface-layer stress at the zlo face
 * @param[in] hfx_z rho<w'theta'> on z-faces; the surface-layer heat flux at the zlo face
 * @param[in] qfx1_z rho<w'qv'> on z-faces; the surface-layer moisture flux at the zlo face
 * @param[in] surface_layer_on_zlo a SurfaceLayer supplies the fluxes through the zlo face
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
ComputeQKESourceTerms (int i, int j, int k,
                       const amrex::Array4<const amrex::Real>& uvel,
                       const amrex::Array4<const amrex::Real>& vvel,
                       const amrex::Array4<const amrex::Real>& cell_data,
                       const amrex::Array4<const amrex::Real>& cell_prim,
                       const amrex::Array4<const amrex::Real>& K_turb,
                       const amrex::Box& domain,
                       const PBLDerivativeDzInv dz_inv,
                       amrex::Real pbl_mynn_B1_l,
                       const amrex::Real theta_mean,
                       const MoistureComponentIndices& moisture_indices,
                       bool c_ext_dir_on_zlo,
                       bool c_ext_dir_on_zhi,
                       bool u_ext_dir_on_zlo,
                       bool u_ext_dir_on_zhi,
                       bool v_ext_dir_on_zlo,
                       bool v_ext_dir_on_zhi,
                       const amrex::Array4<const amrex::Real>& tau13,
                       const amrex::Array4<const amrex::Real>& tau23,
                       const amrex::Array4<const amrex::Real>& hfx_z,
                       const amrex::Array4<const amrex::Real>& qfx1_z,
                       bool surface_layer_on_zlo)
{
    amrex::Real source_term = zero;
 
    int izmin = domain.smallEnd(2);
    int izmax = domain.bigEnd(2);
 
    const amrex::Real rho = cell_data(i,j,k,Rho_comp);
 
    // Notes:
    // - We transport TKE = myhalf*QKE rather than QKE, so the RHS terms do not
    //   have a factor of two.
    // - Transport terms due to turbulence and pressure are included when
    //   DiffusionSrcForState_* is called from ERF_slow_rhs_post.
    // - Eddy diffusivities are updated at the beginning of each time step only.
 
    const bool use_surface_fluxes = ( surface_layer_on_zlo && (k == izmin) &&
                                      tau13 && tau23 && hfx_z );
 
    if (use_surface_fluxes)
    {
        // The first cell center lies inside the unresolved surface layer. With
        // zlo.type = surface_layer the state carries foextrap and the velocities
        // hoextrap at the lower boundary, so the stencil used below reads ghost
        // values that hold no surface information: the resolved gradients there
        // are not a discretization of anything (ERF #4037). Build the production
        // terms from the fluxes the model actually carries instead:
        //
        //   tau13, tau23  = rho<u'w'>, rho<v'w'>      (z-faces, staggered in x,y)
        //   hfx_z, qfx1_z = rho<w'theta'>, rho<w'qv'> (z-faces)
        //
        // Two different heights are wanted here, so the two are kept apart:
        //   - u_*, thetav_* and L are surface-layer parameters and have no
        //     cell-centered counterpart, so they come from the zlo face alone;
        //   - the production terms belong at the cell center, so their flux
        //     factors average the two z-faces bounding this cell.
        // Only the profile shape, phi(zeta_1)/(kappa z_1), ties the two together.
 
        const amrex::Real theta = cell_data(i,j,k,RhoTheta_comp) / rho;
        const amrex::Real qv    = (moisture_indices.qv >= 0) ?
                                  cell_data(i,j,k,moisture_indices.qv) / rho : zero;
 
        // ---- Surface-layer parameters, on the zlo face ----
 
        // |tau|/rho = u_*^2
        const amrex::Real t13_s = myhalf * ( tau13(i,j,k) + tau13(i+1,j,k) );
        const amrex::Real t23_s = myhalf * ( tau23(i,j,k) + tau23(i,j+1,k) );
        const amrex::Real tau_s = std::sqrt(t13_s*t13_s + t23_s*t23_s);
        const amrex::Real ustar2 = tau_s / rho;
        const amrex::Real ustar  = std::sqrt(ustar2);
 
        // Surface virtual heat flux rho<w'theta_v'> (Stull Eqn. 4.4.5d)
        amrex::Real hfxv_s = hfx_z(i,j,k) * (one + epsv*qv);
        if (qfx1_z) { hfxv_s += epsv * theta * qfx1_z(i,j,k); }
 
        // Height of this cell center above the surface. dz_inv.dirichlet_lo is by
        // construction the inverse of the boundary-to-first-center distance, so it
        // carries the right height for uniform, stretched and terrain-fitted grids.
        const amrex::Real z1 = one / dz_inv.dirichlet_lo;
 
        // Obukhov length, from the surface fluxes. A vanishing stress (a building
        // face, say) leaves the shear production zero whatever zeta is, so the
        // neutral value is the safe one to carry into phi_m there.
        const amrex::Real thetav = theta * (one + epsv*qv);
        const amrex::Real eps    = std::numeric_limits<amrex::Real>::epsilon();
        amrex::Real zeta = zero;
        if (std::abs(hfxv_s) > eps && ustar > eps) {
            const amrex::Real l_obukhov = -( thetav * rho * ustar2 * ustar )
                                        / ( KAPPA * CONST_GRAV * hfxv_s );
            zeta = z1 / l_obukhov;
        }
 
        // phi_m, the companion of the psi_m that SurfaceLayer uses
        const similarity_funs sfuns{};
        const amrex::Real phi_m = sfuns.calc_phi_m(zeta);
 
        // ---- Production terms, at the cell center ----
 
        // Cell-centered stress: average the two z-faces bounding this cell
        const amrex::Real t13_c = fourth * ( tau13(i,j,k  ) + tau13(i+1,j,k  )
                                           + tau13(i,j,k+1) + tau13(i+1,j,k+1) );
        const amrex::Real t23_c = fourth * ( tau23(i,j  ,k) + tau23(i,j  ,k+1)
                                           + tau23(i,j+1,k) + tau23(i,j+1,k+1) );
 
        // Cell-centered virtual heat flux. Both faces are from this stage: the
        // QKE sources are added after RhoQ1 is diffused (see AddTurbKESources).
        amrex::Real hfxv_c = myhalf * ( hfx_z(i,j,k) + hfx_z(i,j,k+1) ) * (one + epsv*qv);
        if (qfx1_z) {
            hfxv_c += epsv * theta * myhalf * ( qfx1_z(i,j,k) + qfx1_z(i,j,k+1) );
        }
 
        // Shear production, -<u'w'> dU/dz - <v'w'> dV/dz. MOST sets the direction
        // of the mean shear anti-parallel to the surface stress and its magnitude
        // to u_* phi_m/(kappa z_1), so the dot product is the cell-centered stress
        // projected onto the surface stress. This reduces to |tau| u_* phi_m /
        // (kappa z_1) when the stress does not turn across the cell.
        if (tau_s > eps) {
            const amrex::Real tau_proj = (t13_c*t13_s + t23_c*t23_s) / tau_s;
            source_term += std::max(tau_proj,amrex::Real(0.)) * ustar * phi_m / (KAPPA * z1);
        }
 
        // Buoyancy production, which needs no gradient at all
        source_term += (CONST_GRAV/theta_mean) * hfxv_c;
    }
    else
    {
        // Compute some relevant derivatives
        amrex::Real dthetadz, dudz, dvdz;
        ComputeVerticalDerivativesPBL(i, j, k,
                                      uvel, vvel, cell_data, izmin, izmax, dz_inv,
                                      c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                      u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                      v_ext_dir_on_zlo, v_ext_dir_on_zhi,
                                      dthetadz, dudz, dvdz,
                                      moisture_indices);
 
        // Second-order turbulent fluxes, e.g.:
        // -<uw> = L q SM dU/dz  (NN09, Eqn. 18)
        //       = Kmv/rho dU/dz
 
        // Shear Production
        source_term += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);
 
        // Buoyancy Production
        source_term -= (CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;
    }
 
    // Dissipation (NN09, Eqn. 12)
    amrex::Real qke = two * cell_prim(i,j,k,PrimKE_comp);
    if (std::abs(qke) > zero) {
        source_term -= rho * std::pow(qke,amrex::Real(1.5)) /
                       (pbl_mynn_B1_l * K_turb(i,j,k,EddyDiff::Turb_lengthscale));
    }
 
    return source_term;
}
 
/**
 * @brief Functor for inverse vertical spacings with constant grid spacing.
 */
struct PBLDerivativeDzInv_N
{
    amrex::Real value;
 
    /**
     * Return inverse vertical spacings for constant grid spacing.
     * @return PBLDerivativeDzInv inverse spacing values
     */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    PBLDerivativeDzInv operator() (int, int, int) const noexcept
    {
        return {myhalf*value, two*value, two*value};
    }
};
 
/**
 * @brief Functor for inverse vertical spacings for stretched grids using a spacing array.
 */
struct PBLDerivativeDzInv_S
{
    const amrex::Real* dz;
    int klo;
    int khi;
 
    /**
     * Return inverse vertical spacings for stretched grids.
     * @param k z-index
     * @return PBLDerivativeDzInv inverse spacing values
     */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    PBLDerivativeDzInv operator() (int, int, int k) const noexcept
    {
        // Adjacent cell centers are separated by half the sum of their widths.
        const amrex::Real dz_lo = dz[(k > klo) ? k-1 : k];
        const amrex::Real dz_hi = dz[(k < khi) ? k+1 : k];
        return {two / (dz_lo + two*dz[k] + dz_hi),
                two / dz[k],
                two / dz[k]};
    }
};
 
/**
 * @brief Functor for inverse vertical spacings for terrain-following grids using cell-center heights.
 */
struct PBLDerivativeDzInv_T
{
    amrex::Array4<const amrex::Real> z_cc;
 
    /**
     * Return inverse vertical spacings for terrain-following grids.
     * @param i x-index
     * @param j y-index
     * @param k z-index
     * @return PBLDerivativeDzInv inverse spacing values
     */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    PBLDerivativeDzInv operator() (int i, int j, int k) const noexcept
    {
        return {one / (z_cc(i,j,k+1) - z_cc(i,j,k-1)),
                two / (z_cc(i,j,k  ) - z_cc(i,j,k-1)),
                two / (z_cc(i,j,k+1) - z_cc(i,j,k  ))};
    }
};
 
/**
 * @brief Apply spatial smoothing to PBLH field using 5-point stencil.
 *
 * Applies a simple 5-point spatial averaging to reduce grid-to-grid noise in
 * the diagnosed PBLH field. This is standard post-processing practice for
 * discrete Rib-crossing detection as described in Seibert et al. (2000):
 * "Review and intercomparison of operational methods for the determination of
 * the mixing height", Atmospheric Environment, 34, 1001-1027.
 *
 * Stencil: PBLH_smooth(i,j) = w * PBLH(i,j) + (1-w)/4 * [PBLH(i+-1,j) + PBLH(i,j+-1)]
 *
 * The stencil reads one column outside the region it writes, so the result on a
 * cell depends on columns that a single tile -- or a single box -- does not own.
 * The caller therefore names the region it wants smoothed (xybx_valid, normally
 * the planar tile box) inside a work array carrying `passes` columns of halo
 * around it, filled from the same inputs. Pass p writes the region grown by
 * passes-1-p and reads the one grown by passes-p, so each pass consumes one
 * column of halo and the last lands exactly on xybx_valid. Every read stays
 * inside the array, and the result does not depend on how the domain is divided
 * into boxes or tiles.
 *
 * At a non-periodic domain edge there is no column to read, so the stencil reuses
 * the cell's own value (reflective/Neumann). Across a periodic edge the halo holds
 * the wrapped columns and is read directly, which is why the periodicity has to be
 * passed in: clamping there would fold an edge into a domain that has none.
 *
 * @param[in,out] pblh_fab   FArrayBox holding the 2D PBLH field (component 0),
 *                           sized to contain xybx_valid grown by `passes` in x and y
 * @param[in] xybx_valid     Region to smooth (x,y indices only)
 * @param[in] weight         Center cell weight in stencil (must be in [0,1])
 * @param[in] passes         Number of smoothing iterations
 * @param[in] domain         Domain box, for the reflective edge condition
 * @param[in] periodicity    Domain periodicity, so periodic edges read the halo
 */
inline void
ApplyPBLHSmoothing (amrex::FArrayBox& pblh_fab,
                    const amrex::Box& xybx_valid,
                    const amrex::Real weight,
                    const int passes,
                    const amrex::Box& domain,
                    const amrex::Periodicity& periodicity)
{
    if (passes <= 0) { return; }
 
    // The caller must supply the halo the passes will consume; see the note above.
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        pblh_fab.box().contains(amrex::grow(xybx_valid, amrex::IntVect(passes,passes,0))),
        "ApplyPBLHSmoothing: the PBLH work array needs `passes` columns of halo "
        "around the region being smoothed");
 
    // Use a temporary to hold intermediate results
    amrex::FArrayBox pblh_temp(pblh_fab.box(), 1, amrex::The_Async_Arena());
    auto pblh = pblh_fab.array();
    auto pblh_tmp = pblh_temp.array();
 
    const auto& dom_lo = amrex::lbound(domain);
    const auto& dom_hi = amrex::ubound(domain);
 
    const bool per_x = periodicity.isPeriodic(0);
    const bool per_y = periodicity.isPeriodic(1);
 
    const amrex::Real wt_side = (one - weight) / four; // Weight for each of 4 neighbors
 
    // Apply smoothing passes
    for (int pass = 0; pass < passes; ++pass) {
        const int halo = passes - 1 - pass;
        const amrex::Box bx = amrex::grow(xybx_valid, amrex::IntVect(halo,halo,0));
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            // Reflective (Neumann) at a non-periodic domain edge: reuse the cell's
            // own value. Everywhere else -- including across a periodic edge, where
            // the halo holds the wrapped column -- read the neighbour.
            const int i_xlo = (per_x || i > dom_lo.x) ? (i - 1) : i;
            const int i_xhi = (per_x || i < dom_hi.x) ? (i + 1) : i;
            const int j_ylo = (per_y || j > dom_lo.y) ? (j - 1) : j;
            const int j_yhi = (per_y || j < dom_hi.y) ? (j + 1) : j;
 
            // 5-point stencil
            pblh_tmp(i, j, 0) = weight * pblh(i, j, 0)
                              + wt_side * (pblh(i_xlo, j, 0) + pblh(i_xhi, j, 0) +
                                          pblh(i, j_ylo, 0) + pblh(i, j_yhi, 0));
        });
 
        // Copy back only what this pass computed: outside bx the array still holds
        // the previous pass, which no later pass reads.
        pblh_fab.copy<amrex::RunOn::Device>(pblh_temp, bx, 0, bx, 0, 1);
    }
}
 
 
/**
 * Copy of a field read by the PBLH passes, widened to `ng` horizontal ghost cells.
 *
 * With erf.pblh_smoothing_passes = n the PBLH passes read their inputs n columns outside
 * the tile (see ApplyPBLHSmoothing), but some of those inputs carry fewer ghost cells: the
 * surface-layer fields have one on the bottom surface. The PBL schemes read such fields
 * from a copy made here instead.
 *
 * The valid cells and the field's own ghost cells are copied unchanged, so a run that
 * needs no more halo than the field has gives the same answer as reading the field
 * itself. `fill` then fills the new ghost cells that a box of the level covers, inside
 * the domain or across a periodic edge, the same way the field's owner fills its own
 * (for a surface-layer field, SurfaceLayer::fill_planar_boundary). Ghost cells no box
 * covers -- outside a non-periodic domain edge, or next to a coarser level -- take the
 * value of the nearest cell of the field's original halo. The smoothing never reads a
 * column outside a non-periodic edge; next to a coarser level this is a zeroth-order
 * extrapolation of the surface scales, while the state there is the interpolated one.
 *
 * @param[in] src   Field to widen
 * @param[in] ng    Ghost cells wanted in x and y; z keeps those of src
 * @param[in] fill  Callable taking the widened MultiFab& and filling its ghost cells
 *                  from the valid cells of the other boxes of the level
 *
 * `fill` is a std::function, not a template parameter: nvcc does not accept the device
 * lambda below inside a function template instantiated with a lambda local to the caller.
 */
inline std::unique_ptr<amrex::MultiFab>
WidenForPBLH (const amrex::MultiFab& src, const int ng,
              const std::function<void(amrex::MultiFab&)>& fill)
{
    const amrex::IntVect ng_src = src.nGrowVect();
    const amrex::IntVect ng_dst(ng, ng, ng_src[2]);
    const int ncomp = src.nComp();
 
    // A value no field holds, marking the cells nothing below fills
    constexpr amrex::Real unset = std::numeric_limits<amrex::Real>::max();
 
    auto dst = std::make_unique<amrex::MultiFab>(src.boxArray(), src.DistributionMap(),
                                                 ncomp, ng_dst);
    dst->setVal(unset);
    amrex::MultiFab::Copy(*dst, src, 0, 0, ncomp, ng_src);
    fill(*dst);
 
    for (amrex::MFIter mfi(*dst); mfi.isValid(); ++mfi) {
        // The cells holding src's own data; never written below, so they can be read
        const amrex::Box have = amrex::grow(mfi.validbox(), ng_src);
        const auto lo = amrex::lbound(have);
        const auto hi = amrex::ubound(have);
        auto a = dst->array(mfi);
        amrex::ParallelFor(mfi.fabbox(), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            if (have.contains(i, j, k)) { return; }
            const int ic = amrex::min(amrex::max(i, lo.x), hi.x);
            const int jc = amrex::min(amrex::max(j, lo.y), hi.y);
            for (int n = 0; n < ncomp; ++n) {
                if (a(i, j, k, n) == unset) { a(i, j, k, n) = a(ic, jc, k, n); }
            }
        });
    }
    return dst;
}
 
#endif