ERF_ComputeDiffusivityMRF.cpp File Reference

#include <algorithm>
#include "ERF_SurfaceLayer.H"
#include "ERF_DirectionSelector.H"
#include "ERF_Diffusion.H"
#include "ERF_Constants.H"
#include "ERF_TurbStruct.H"
#include "ERF_PBLModels.H"
#include "ERF_TileNoZ.H"
#include "ERF_MoistUtils.H"
#include "ERF_PBLScaleAwareBlending.H"

Include dependency graph for ERF_ComputeDiffusivityMRF.cpp:

Functions
void	ComputeDiffusivityMRF (const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const std::unique_ptr< MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices, const MultiFab *terrain_blank)

Function Documentation

ComputeDiffusivityMRF()

void ComputeDiffusivityMRF	(	const MultiFab &	xvel,
		const MultiFab &	yvel,
		const MultiFab &	cons_in,
		MultiFab &	eddyViscosity,
		const Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const std::unique_ptr< MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices,
		const MultiFab *	terrain_blank
	)

Compute vertical diffusivity using the Medium-Range Forecast (MRF) boundary layer scheme.

Parameters

[in]	xvel	x-velocity field.
[in]	yvel	y-velocity field.
[in]	cons_in	Input conserved variables.
[out]	eddyViscosity	Eddy viscosity and diffusivity coefficients.
[in]	geom	Grid geometry.
[in]	turbChoice	Turbulence model options.
[in]	SurfLayer	Surface layer data.
[in]	use_terrain_fitted_coords	Use terrain-fitted coordinates.
[in]	use_moisture	Enable moisture components.
[in]	level	Current AMR level.
[in]	bc_ptr	Boundary condition records.
[in]	z_phys_nd	Nodal physical heights.
[in]	z_phys_cc	Cell-centered physical heights.
[in]	moisture_indices	Indices for moisture variables.

{
    const StratType pbl_strat = turbChoice.pbl_strat_type;
    /*
    ============================================================================
    Medium-Range Forecast (MRF) Boundary Layer Parameterization Scheme
    ============================================================================
 
    Implementation of the MRF (Medium Range Forecast) boundary layer scheme
    based on Hong and Pan (1996), with enhanced moisture handling and
    cloud-aware stability corrections developed for ERF.
 
    PRIMARY REFERENCES:
    -------------------
    - 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
 
    - WRF Reference Implementation: module_bl_mrf.F
      https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F
 
    - Hong, S. Y., Y. Noh, and J. Dudhia, 2006: A new vertical diffusion
      package with an explicit treatment of entrainment processes. Monthly
      Weather Review, 134, 2318-2341. https://doi.org/10.1175/MWR3250.1
      (YSU scheme, used for free atmosphere mixing)
 
    CORE ALGORITHM (Hong & Pan 1996):
    ----------------------------------
    1. Predictor-Corrector Bulk Richardson Number (Rib) PBL Height Diagnosis
       - Rib = (g*z*(θ_v(z) - θ_v(0))) / (U²*θ_v(0))
       - h = min(z where Rib > Rib_critical)
 
    2. Nonlocal Countergradient Flux Corrections
       - HGAMT: thermal excess from sensible heat flux
       - HGAMQ: moisture excess from latent heat flux
       - VPERT: virtual potential temperature perturbation
 
    3. Stability-Dependent K-Profile in the PBL (used in EVERY stability regime)
       - K_m = ρ * w_s * κ * z * (1 - (z - z_1)/(h - z_1))², z_1 = first level height
       - K_t = K_m / Pr_t, and K_q = K_t (WRF diffuses moisture with the heat
         coefficient XKZH; module_bl_mrf.F L1069, L1180)
       - Stability enters only through w_s = u_* / φ_m (WRF's WSCALE); as in WRF there
         is no regime test below the PBL top (module_bl_mrf.F L968-986).  The
         Ri-based mixing ERF used to apply inside stable/neutral PBLs is gone with
         it; see KNOWN LIMITATION below for what that costs in stable layers.
 
    4. Free Atmosphere Mixing via Richardson Number
       - Uses YSU stability functions (Hong et al. 2006)
       - Ri-dependent diffusivity above the PBL top only (module_bl_mrf.F L988 on)
 
    ENHANCEMENTS IN ERF IMPLEMENTATION:
    -----------------------------------
 
    A. Virtual Potential Temperature (θ_v) Treatment
       - Proper handling of moisture effects on buoyancy
       - Important for accurate PBL height and convective strength
       - Not explicitly detailed in Hong & Pan (1996) but standard practice
 
    B. Cloud-Aware Stability Function Adjustments (OPTIONAL)
       - Physically motivated extension beyond Hong & Pan (1996)
       - Detects clouds via qc + qi > threshold (0.1 g/kg)
       - In stable layers with clouds:
         * Reduces stability damping (φ reduced by 10-20%)
         * Represents suppressed turbulence due to stable stratification
       - In unstable layers with clouds:
         * Slightly increases instability enhancement
         * Represents latent heat release effects
       - Disabled by default for backward compatibility
       - Reference: Conceptually similar to WRF's IMVDIF cloud handling
 
    C. Explicit Moisture Sign Convention Handling
       - q_star uses WRF convention: negative for upward flux (evaporation)
       - HGAMQ calculation: -const_b * u_* * q_* / w_*
       - Produces positive HGAMQ for unstable (evaporating) conditions
       - Corrected with MAX(HGAMQ, 0) to prevent negative values
       - Land/water discrimination: HGAMQ zeroed over water surfaces
 
    D. Improved VPERT (Virtual Potential Temperature Perturbation)
       Formulation
       - VPERT = max(HGAMT + 0.61*θ*HGAMQ, 0)
       - Enhancement vs WRF: Does not limit VPERT to GAMCRT
       - WRF limits VPERT after adding moisture term (VPERT = min(VPERT, GAMCRT))
       - This is physically incorrect: moisture contribution can exceed GAMCRT
       - ERF only limits HGAMT (sensible heat), allowing larger VPERT when
         both sensible and latent heat create strong surface heating
       - This produces more realistic PBL heights in moist convection
 
    E. Prandtl Number with Stability Correction
       - Pr_t = φ_t/φ_m + const_b * κ * sf
       - Bounded: 0.5 ≤ Pr_t ≤ 4.0
       - Consistency with WRF: const_b = 7.8, sf = 0.1
       - Reference: Hong et al. (2006), Equation A17
 
    PARAMETER DEFAULTS (Matching WRF):
    -----------------------------------
    - const_b = 7.8 (heat flux weight in countergradient)
    - sf = 0.1 (surface layer / PBL height ratio in stability functions)
    - GAMCRT = 3.0 K (maximum heat countergradient)
    - GAMCRQ = 2.0e-3 kg/kg (maximum moisture countergradient)
    - Rib_crit = 0.5 (critical bulk Richardson number for PBL height)
    - Pr_min = 0.5, Pr_max = 4.0 (Prandtl number bounds)
    - Kmin = 0.1 m²/s, Kmax = 300 m²/s (diffusivity bounds, Hong & Pan 1996)
 
    STABILITY FUNCTIONS:
    --------------------
    Unstable (L < 0, HOL < 0):
      phi_m = (1 - 8 * sf * h/L)^(-1/3)
      phi_t = (1 - 16 * sf * h/L)^(-1/2)
 
    Stable (L > 0, HOL > 0):
      Default: phi_m = phi_t = 1 + 5 * sf * h/L
      QNSE (if enable_qnse_stable_functions=true):
        phi_m = (1 + qnse_am * zeta) / (1 + qnse_bm * zeta)
        phi_t = (1 + qnse_ah * zeta) / (1 + qnse_bh * zeta)
        where zeta = sf * h/L, and default coefficients are from
        Sukoriansky, Galperin & Perov (2005): qnse_am=2.5, qnse_bm=0.2,
        qnse_ah=2.5, qnse_bh=0.2. QNSE functions remain bounded as zeta→∞,
        avoiding over-mixing suppression in very stable conditions.
 
    MIXING ABOVE PBL (Free Atmosphere):
    -----------------------------------
    Uses YSU scheme (Hong et al. 2006, Appendix A) for Richardson number
    dependent mixing. This avoids MRF oscillations in stable conditions.
    It is applied above the PBL top only, as in WRF; below the top the
    K-profile of item 3 above is used in every stability regime.
 
    Gradient Richardson number: Ri_g = (g/theta_v) * (dtheta_v/dz) / ((du/dz)^2 + (dv/dz)^2)
 
    For Ri_g > 0 (stable):
      f_m = 1 / ((1 + 5*Ri_g)^2)
      f_t = 1 / ((1 + 5*Ri_g)^2)
      Pr_t = 1 + 2.1*Ri_g, bounded to [0.25, 4.0]
 
    For Ri_g < 0 (unstable):
      f_m = 1 - 8*Ri_g / (1 + 1.746*sqrt(-Ri_g))
      f_t = 1 - 8*Ri_g / (1 + 1.286*sqrt(-Ri_g))
 
    NUMERICAL CONSIDERATIONS:
    -------------------------
    - GPU-optimized with parallel_for loops
    - Safe bounds on shear and vertical derivatives
    - Gradient Richardson number limited: -100 <= Ri_g <= 100
    - Avoids division by zero with minimum shear threshold (1e-10)
    - theta_v_klo and theta_v guarded against zero/negative values
 
    KNOWN LIMITATION (stable layers):
    ---------------------------------
    With the K-profile applied in every regime, this scheme inherits WRF MRF's
    over-mixing of stable boundary layers: on GABLS1 its peak K_m is roughly
    2.1x YSUNew.  That is the behaviour Hong et al. (2006) set out to fix in YSU,
    so stable cases that need the weaker mixing should use the YSUNew scheme
    (erf.pbl_type = YSUNew).
 
    TESTING & VALIDATION:
    ---------------------
    - Tested against WRF_ideal unstable/stable ABL cases
    - Backward compatible when enhancements are disabled
    - Cloud detection threshold: qc + qi > 1e-4 kg/kg (0.1 g/kg)
    */
 
    // Domain extent in z-dir
    int klo = geom.Domain().smallEnd(2);
    int khi = geom.Domain().bigEnd(2);
 
 
 
// Collect MRF-computed PBLH for writing back to SurfaceLayer
MultiFab pblh_mf(eddyViscosity.boxArray(), eddyViscosity.DistributionMap(), 1, 0);
pblh_mf.setVal(0.0);
 
// PBLH smoothing reads one column per pass outside the cells it writes, so the
// planar work arrays have to carry that many columns of halo and the PBLH passes
// have to fill them: otherwise the stencil reads off the end of the tile and the
// answer depends on the decomposition. With smoothing off this is the single
// column of halo the diffusivity kernels already use, so nothing changes.
const int ng_pblh = (turbChoice.enable_pblh_smoothing)
                  ? std::max(1, turbChoice.pblh_smoothing_passes) : 1;
if (ng_pblh > 1) {
    const int ng_avail = std::min({cons_in.nGrowVect()[0], cons_in.nGrowVect()[1],
                                   xvel.nGrowVect()[0],    xvel.nGrowVect()[1],
                                   yvel.nGrowVect()[0],    yvel.nGrowVect()[1],
                                   SurfLayer->get_u_star(level)->nGrowVect()[0],
                                   SurfLayer->get_u_star(level)->nGrowVect()[1],
                                   SurfLayer->get_olen(level)->nGrowVect()[0],
                                   SurfLayer->get_olen(level)->nGrowVect()[1],
                                   (terrain_blank) ? terrain_blank->nGrowVect()[0] : 1000,
                                   (terrain_blank) ? terrain_blank->nGrowVect()[1] : 1000});
    if (ng_pblh > ng_avail) {
        amrex::Abort("erf.pblh_smoothing_passes = " + std::to_string(turbChoice.pblh_smoothing_passes)
                   + " needs " + std::to_string(ng_pblh) + " halo columns, but the state and "
                     "surface-layer arrays carry only " + std::to_string(ng_avail)
                   + "; reduce erf.pblh_smoothing_passes to at most " + std::to_string(ng_avail));
    }
}
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
 
    for (MFIter mfi(eddyViscosity, TileNoZ()); mfi.isValid(); ++mfi) {
 
        // Box operated on must span fill domain in z-dir
        const Box& gbx = mfi.growntilebox(IntVect(1,1,0));
        AMREX_ALWAYS_ASSERT( gbx.smallEnd(2) == klo &&
                             gbx.bigEnd(2)   == khi );
 
        // create flattened boxes to store PBL height and related quantities
        const GeometryData gdata = geom.data();
        const Box xybx = PerpendicularBox<ZDir>(gbx, IntVect{0, 0, 0});
 
        // The PBLH passes run on the region the smoothing stencil will consume
        // (see ng_pblh above). growntilebox() is no use here: it does not grow at
        // an interior tile edge, so a tile in the middle of a box gets no halo at
        // all -- which is exactly how the stencil came to read off the end of the
        // array. Grow the tile box explicitly instead, which costs a little
        // duplicated work in the overlaps. With smoothing off this is xybx and
        // nothing is duplicated.
        const Box gbx_work  = (turbChoice.enable_pblh_smoothing)
                            ? amrex::grow(mfi.tilebox(), IntVect(ng_pblh,ng_pblh,0)) : gbx;
        const Box xybx_work = PerpendicularBox<ZDir>(gbx_work, IntVect{0, 0, 0});
        const Box xybx_tile = PerpendicularBox<ZDir>(mfi.tilebox(), IntVect{0, 0, 0});
 
        // Pass 1 (predictor): PBLH with base surface temperature, no VPERT
        // Pass 2 (wstar/VPERT): compute wstar, HGAMT, HGAMQ, VPERT from predictor height
        // Pass 3 (corrector): PBLH with VPERT-enhanced surface temperature (WRF-consistent)
        //   WRF reference (module_bl_mrf.F lines 813-964):
        //   https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L813-L964
        FArrayBox pbl_height_predictor(xybx_work, 1, The_Async_Arena());  // Pass 1: base t_layer_v
        FArrayBox pbl_height_corrector(xybx_work, 1, The_Async_Arena());  // Pass 3: VPERT-enhanced
        IArrayBox pbl_index(xybx_work, 1, The_Async_Arena());
        IArrayBox pbl_index_zero_ri(xybx_work, 1, The_Async_Arena());  // Index for zero-Ri diagnostic pass
        FArrayBox hgamt_fab(xybx_work, 1, The_Async_Arena());  // Store HGAMT/h (normalized countergradient)
        FArrayBox hgamq_fab(xybx_work, 1, The_Async_Arena());  // Store HGAMQ/h (normalized countergradient)
        FArrayBox wstar_fab(xybx_work, 1, The_Async_Arena());  // Convective velocity scale
        FArrayBox vpert_fab(xybx_work, 1, The_Async_Arena());  // Virtual temperature perturbation VPERT
        FArrayBox pblh_floor_fab(xybx_work, 1, The_Async_Arena());  // Floor of the corrected height per column, reapplied after the smoothing
        const auto& pblh_pred_arr   = pbl_height_predictor.array();  // predictor (base t_layer_v)
        const auto& pblh_corr_arr   = pbl_height_corrector.array();  // corrector (VPERT-enhanced)
        const auto& pblh_floor_arr  = pblh_floor_fab.array();
        const auto& pbli_arr        = pbl_index.array();
        const auto& pbli_zero_arr   = pbl_index_zero_ri.array();  // Zero-Ri diagnostic PBL index
        const auto& hgamt_arr       = hgamt_fab.array();
        const auto& hgamq_arr       = hgamq_fab.array();
        const auto& wstar_arr       = wstar_fab.array();
        const auto& vpert_arr       = vpert_fab.array();
 
        // Get some data in arrays
        const auto& cell_data = cons_in.const_array(mfi);
        const auto& uvel = xvel.const_array(mfi);
        const auto& vvel = yvel.const_array(mfi);
 
        const Real Ribcr = turbChoice.pbl_mrf_Ribcr;
        const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
        const auto& t_star_arr = SurfLayer->get_t_star(level)->const_array(mfi);
        const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
        const auto& t10av_arr  = SurfLayer->get_mac_avg(level, 3)->const_array(mfi);
        const auto& q10av_arr  = SurfLayer->get_mac_avg(level, 4)->const_array(mfi);
        const auto& lmask_arr = (SurfLayer->get_lmask(level)) ?
                                SurfLayer->get_lmask(level)->const_array(mfi) :
                                Array4<int>{};
        // Only retrieve z_phys_nd array if terrain-fitted coordinates are in use
        const Array4<Real const> z_nd_arr = use_terrain_fitted_coords ? z_phys_nd->array(mfi)
                                                                      : Array4<Real const>{};
        const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
 
        //
        // PASS 1 (PREDICTOR): Compute PBL height using base surface virtual temperature.
        // No thermal excess (VPERT = 0). Produces pblh_pred_arr and pbli_arr.
        // WRF reference (module_bl_mrf.F lines 813-842):
        // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L813-L842
        //
        // FIX: guard against null dereference of q_star when use_moisture is false (dry case)
        // get_q_star can return null pointer when the field is not allocated
        amrex::MultiFab* q_star_mf = SurfLayer->get_q_star(level);
        Array4<Real const> q_star_arr = (q_star_mf != nullptr) ? q_star_mf->const_array(mfi)
                                                               : Array4<Real const>{};
 
        // --------------------------------------------------------------------
        // Immersed-boundary awareness (turbChoice.pbl_ib_aware): per column,
        // the first fluid cell above the immersed solid (ksurf) and its height
        // above the domain bottom (zib). Every height below is measured from
        // it, the diffusivities vanish inside the solid, and the surface
        // scales of a column with solid cells are a neutral log law at its
        // top instead of the ground surface layer's values (which come from
        // cells inside the solid there). With the option off, or without
        // immersed cells, ksurf = klo, zib = 0 and the scales are copies, so
        // the arithmetic below is unchanged.
        // --------------------------------------------------------------------
        const bool l_ib = turbChoice.pbl_ib_aware && (terrain_blank != nullptr);
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(l_ib && use_terrain_fitted_coords),
            "erf.pbl_ib_aware is not supported with terrain-fitted coordinates");
        const Array4<Real const> blank_arr = l_ib ? terrain_blank->const_array(mfi) : Array4<Real const>{};
        IArrayBox ksurf_fab(xybx_work, 1, The_Async_Arena());
        FArrayBox zib_fab(xybx_work, 1, The_Async_Arena());
        FArrayBox us_eff_fab(xybx_work, 1, The_Async_Arena()), ts_eff_fab(xybx_work, 1, The_Async_Arena());
        FArrayBox qs_eff_fab(xybx_work, 1, The_Async_Arena()), ol_eff_fab(xybx_work, 1, The_Async_Arena());
        FArrayBox t10_eff_fab(xybx_work, 1, The_Async_Arena()), q10_eff_fab(xybx_work, 1, The_Async_Arena());
        const auto& ksurf_arr = ksurf_fab.array();
        const auto& zib_arr   = zib_fab.array();
        const auto& us_eff_arr = us_eff_fab.array();  const auto& ts_eff_arr = ts_eff_fab.array();
        const auto& qs_eff_arr = qs_eff_fab.array();  const auto& ol_eff_arr = ol_eff_fab.array();
        const auto& t10_eff_arr = t10_eff_fab.array(); const auto& q10_eff_arr = q10_eff_fab.array();
        {
            const Real dz_ib = geom.CellSize(2);
            const Real z0_ib = turbChoice.pbl_ib_z0;
            const Array4<Real const>& uv_ib = xvel.const_array(mfi);
            const Array4<Real const>& vv_ib = yvel.const_array(mfi);
            const bool have_qstar = (q_star_mf != nullptr);
            ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
            {
                int ks = klo;
                if (l_ib) {
                    for (int kk = klo; kk <= khi; ++kk) { if (blank_arr(i, j, kk) >= Real(0.5)) { ks = kk + 1; } }
                    if (ks > khi) { ks = khi; }
                }
                ksurf_arr(i, j, 0) = ks;
                zib_arr(i, j, 0)   = (ks - klo) * dz_ib;
                if (ks > klo) {
                    // Column with immersed solid: neutral log law at its top.
                    const Real u_top = myhalf * (uv_ib(i, j, ks) + uv_ib(i + 1, j, ks));
                    const Real v_top = myhalf * (vv_ib(i, j, ks) + vv_ib(i, j + 1, ks));
                    const Real ws_top = std::sqrt(u_top * u_top + v_top * v_top);
                    us_eff_arr(i, j, 0)  = amrex::max(KAPPA * ws_top / std::log(myhalf * dz_ib / z0_ib), Real(1.0e-3));
                    ts_eff_arr(i, j, 0)  = Real(0);
                    qs_eff_arr(i, j, 0)  = Real(0);
                    ol_eff_arr(i, j, 0)  = Real(1.0e10);
                    t10_eff_arr(i, j, 0) = cell_data(i, j, ks, RhoTheta_comp) / cell_data(i, j, ks, Rho_comp);
                    q10_eff_arr(i, j, 0) = (use_moisture && moisture_indices.qv >= 0)
                                         ? cell_data(i, j, ks, moisture_indices.qv) / cell_data(i, j, ks, Rho_comp) : Real(0);
                } else {
                    us_eff_arr(i, j, 0)  = u_star_arr(i, j, 0);
                    ts_eff_arr(i, j, 0)  = t_star_arr(i, j, 0);
                    qs_eff_arr(i, j, 0)  = have_qstar ? q_star_arr(i, j, 0) : Real(0);
                    ol_eff_arr(i, j, 0)  = l_obuk_arr(i, j, 0);
                    t10_eff_arr(i, j, 0) = t10av_arr(i, j, 0);
                    q10_eff_arr(i, j, 0) = use_moisture ? q10av_arr(i, j, 0) : Real(0);
                }
            });
        }
 
        ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            const int  ksrf = ksurf_arr(i, j, 0);   // local surface: first fluid cell of the column
            const Real zib  = zib_arr(i, j, 0);     // its height, subtracted from every height below
            amrex::ignore_unused(zib);
            const Real t_layer = t10_eff_arr(i, j, 0);
            const Real moisture_fraction = use_moisture ? q10_eff_arr(i, j, 0) : Real(0);
            const Real t_layer_v = t_layer * (one + epsv * moisture_fraction);
 
            Real zval, Rib;
            int kpbl = ksrf;
 
            {
                zval = (use_terrain_fitted_coords)
                     ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
                     : ((kpbl + myhalf) * gdata.CellSize(2) - zib);
                const Real theta_v     = GetThetav(i, j, kpbl, cell_data, moisture_indices);
                const Real theta_v_klo = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), Real(1.0));
                const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
                                                (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
                // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
                // to better represent shear associated with surface-layer turbulence at low wind speeds.
                // Reference: Vogelezang, D.H.P., and A.A.M. Holtslag, 1996: Evaluation and model impacts of
                // alternative boundary-layer height formulations. Boundary-Layer Meteorology, 81, 245–269.
                const Real ws2 = (turbChoice.enable_vh96_shear_correction)
                                ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
                                : amrex::max(ws2_raw, Real(1.0));
                Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
            }
 
            Real zval0 = zval, Rib0 = Rib;
            bool above_critical = false;
            while (!above_critical && ((kpbl + 1) <= khi)) {
                zval0 = zval;
                Rib0  = Rib;
                kpbl += 1;
 
                zval = (use_terrain_fitted_coords)
                     ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
                     : ((kpbl + myhalf) * gdata.CellSize(2) - zib);
                const Real theta_v     = GetThetav(i, j, kpbl, cell_data, moisture_indices);
                const Real theta_v_klo = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), Real(1.0));
                const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
                                                (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
                // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
                // to better represent shear associated with surface-layer turbulence at low wind speeds.
                const Real ws2 = (turbChoice.enable_vh96_shear_correction)
                                ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
                                : amrex::max(ws2_raw, Real(1.0));
                Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
                above_critical = (Rib >= Ribcr);
            }
 
            const Real pblh_emp = (use_terrain_fitted_coords)
                                ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
                                : myhalf * gdata.CellSize(2);
            const Real z_max = (use_terrain_fitted_coords)
                             ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr)
                             : ((khi + myhalf) * gdata.CellSize(2) - zib);
            const Real pblh_max = Real(0.9) * z_max;
            const Real pblh_min = amrex::max(pblh_emp, Real(10.0));
 
            if (above_critical) {
                Real pblh_interp;
                const Real rib_diff = Rib - Rib0;
                if (std::abs(rib_diff) > Real(1.0e-10)) {
                    pblh_interp = zval0 + (zval - zval0) / rib_diff * (Ribcr - Rib0);
                } else {
                    // If Rib not changing, use zval0 as PBL height
                    pblh_interp = zval0;
                }
                pblh_pred_arr(i, j, 0) = amrex::max(amrex::min(pblh_interp, pblh_max), pblh_min);
                pbli_arr(i, j, 0) = kpbl;
            } else {
                pblh_pred_arr(i, j, 0) = pblh_min;
                pbli_arr(i, j, 0) = ksrf + 1;
            }
        });
 
        const Real const_b = turbChoice.pbl_mrf_const_b;
        const Real sf = turbChoice.pbl_mrf_sf;
        constexpr Real prmin = Real(0.5);
        constexpr Real prmax = Real(4.0);
        constexpr Real GAMCRT = Real(3.0);
        constexpr Real GAMCRQ = Real(2.e-3);
        const bool enable_mrf_countergradient = turbChoice.enable_mrf_countergradient;
        const bool enable_mrf_unbounded_vpert = turbChoice.enable_mrf_unbounded_vpert;
 
        //
        // PASS 2 (WSTAR / VPERT): Compute wstar, HGAMT, HGAMQ, VPERT using the
        // predictor PBL height (pblh_pred_arr). VPERT will be fed into the corrector
        // Rib search in Pass 3 to raise the effective surface temperature.
        // WRF reference (module_bl_mrf.F lines 857-880):
        // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L857-L880
        //
        ParallelFor(xybx_work, [=,zero_d=zero] AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            const int  ksrf = ksurf_arr(i, j, 0);   // local surface: first fluid cell of the column
            const Real zib  = zib_arr(i, j, 0);     // its height, subtracted from every height below
            amrex::ignore_unused(zib);
            const Real t_layer  = t10_eff_arr(i, j, 0);
            Real obuk_val = ol_eff_arr(i, j, 0);
            if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
                obuk_val = (obuk_val >= zero_d) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
            }
 
            // HOL computed from predictor height (WRF uses predictor height for wstar/VPERT)
            const Real HOL = sf * pblh_pred_arr(i, j, 0) / obuk_val;
            const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
            // |HOL|, so the base 1 + 16*HOL_abs of the unstable arm below equals its
            // 1 - 16*HOL there and is at least 1 for every HOL; a guard written inside the
            // arm gets folded away by the optimiser and the pow hoisted (see sqrt_neg_Ri).
            const Real HOL_abs = std::abs(HOL_bounded);
            const Real one_quarter = Real(0.25);
            const Real phiM = (obuk_val > 0)
                            ? (1 + 5 * HOL_bounded)
                            : std::pow(1 + 16 * HOL_abs, -one_quarter);
            const Real phiM_safe = amrex::max(phiM, Real(0.01));
 
            // wstar = u* / phi_m
            // Absolute bounds [0.01, 5.0] m/s: prevent division-by-zero (floor) and
            // free-convection blow-up (ceiling), independent of u* for u*->0 safety.
            // BOOTSTRAP NOTE: pblh_pred_arr must be finite (not sentinel) before this
            // kernel runs. On cold start, pblh_pred_arr is computed in Pass 1 above;
            // on regrid, SurfaceLayer PBLH may be sentinel — callers must re-run
            // ComputeDiffusivityMRF before consuming wstar or Beljaars correction after regrid.
            // Deardorff (1972) convective velocity scale — well-behaved as u_star -> 0.
            // wstar_conv = (g/theta_v0 * zi * wthv0)^(1/3)
            const Real theta_v0 = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), Real(1.0));
            const Real wthv0 = amrex::max(-us_eff_arr(i,j,0) * ts_eff_arr(i,j,0), Real(0.0)); // surface buoyancy flux, unstable only
            const Real zi = amrex::max(pblh_pred_arr(i, j, 0), Real(10.0));  // predictor: Pass 3 has not run yet
            const Real wstar_conv = std::cbrt(CONST_GRAV / theta_v0 * zi * wthv0);
 
            Real wstar_shear = us_eff_arr(i, j, 0) / phiM_safe;
 
            // Blend: standard convective/shear combination (WRF YSU-style),
            // weighted by cube sum so either term can dominate smoothly
            Real wstar = std::cbrt(wstar_shear*wstar_shear*wstar_shear +
                                wstar_conv*wstar_conv*wstar_conv);
 
            wstar = amrex::max(wstar, Real(0.01));
            wstar = amrex::min(wstar, Real(5.0));
 
            bool SFCFLG = (obuk_val <= zero_d);
            const Real HGAMT = (SFCFLG && enable_mrf_countergradient)
                             ? amrex::min(-const_b * us_eff_arr(i, j, 0) * ts_eff_arr(i, j, 0) / wstar, GAMCRT)
                             : zero_d;
 
            Real HGAMQ = Real(0);
            if (SFCFLG && use_moisture && enable_mrf_countergradient) {
                const Real q_star = qs_eff_arr(i, j, 0);
                const Real HGAMQ_calc = -const_b * us_eff_arr(i, j, 0) * q_star / wstar;
                HGAMQ = amrex::max(amrex::min(HGAMQ_calc, GAMCRQ), Real(0));
 
                if (lmask_arr) {
                    bool is_land = (lmask_arr(i,j,0) == 1);
                    if (!is_land) HGAMQ = zero_d;
                }
 
                if (moisture_indices.qv >= 0) {
                    Real qv_klo = cell_data(i, j, ksrf, moisture_indices.qv) / cell_data(i, j, ksrf, Rho_comp);
                    Real T_klo = getTgivenRandRTh(cell_data(i, j, ksrf, Rho_comp),
                                                  cell_data(i, j, ksrf, RhoTheta_comp), qv_klo);
                    Real p_klo = getPgivenRTh(cell_data(i, j, ksrf, RhoTheta_comp), qv_klo) * Real(0.01);
                    Real qsat_klo = zero_d;
                    erf_qsatw(T_klo, p_klo, qsat_klo);
                    Real rh_klo = (qsat_klo > Real(1.0e-10)) ? (qv_klo / qsat_klo) : Real(0);
                    if (rh_klo > Real(0.95)) {
                        Real rh_scaling = amrex::max(zero_d, (Real(1) - rh_klo) / Real(0.05));
                        HGAMQ *= rh_scaling;
                    }
                }
            }
 
            // Compute VPERT = HGAMT + 0.61*theta*HGAMQ (virtual temperature perturbation)
            // This will be added to the surface temperature in the corrector Rib search.
            // WRF Reference: module_bl_mrf.F lines 879-880
            if (pbli_arr(i, j, 0) <= ksrf + 1 || !enable_mrf_countergradient) {
                vpert_arr(i, j, 0) = zero_d;
            } else {
                const Real VPERT_raw = HGAMT + epsv * t_layer * HGAMQ;
                const Real VPERT_capped = enable_mrf_unbounded_vpert
                                        ? VPERT_raw
                                        : amrex::min(VPERT_raw, GAMCRT);
                vpert_arr(i, j, 0) = amrex::max(VPERT_capped, zero_d);
            }
        });
 
        //
        // PASS 3 (CORRECTOR): Recompute PBL height with VPERT-enhanced surface temperature.
        // theta_s = theta_va + VPERT  (Hong & Pan 1996, Eq. 4)
        // This is the WRF-consistent corrector pass. pblh_corr_arr is used for the
        // K-profile shape function K = rho*wstar*kappa*z*(1-z/h)^2.
        // WRF reference (module_bl_mrf.F lines 932-964):
        // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L932-L964
        //
        ParallelFor(xybx_work, [=,one_d=one]
                    AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            const int  ksrf = ksurf_arr(i, j, 0);   // local surface: first fluid cell of the column
            const Real zib  = zib_arr(i, j, 0);     // its height, subtracted from every height below
            amrex::ignore_unused(zib);
            const Real t_layer  = t10_eff_arr(i, j, 0);
            const Real moisture_fraction = use_moisture ? q10_eff_arr(i, j, 0) : Real(0);
            // VPERT-enhanced surface virtual temperature (WRF THERMAL variable)
            const Real t_layer_v = t_layer * (one + epsv * moisture_fraction)
                                 + vpert_arr(i, j, 0);
 
            Real obuk_val = ol_eff_arr(i, j, 0);
            if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
                obuk_val = (obuk_val >= zero) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
            }
 
            int kpbl = ksrf;
            Real zval, Rib;
            {
                zval = (use_terrain_fitted_coords)
                     ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
                     : ((kpbl + myhalf) * gdata.CellSize(2) - zib);
                const Real theta_v     = GetThetav(i, j, kpbl, cell_data, moisture_indices);
                const Real theta_v_klo = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), one_d);
                const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
                                                (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
                // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
                const Real ws2 = (turbChoice.enable_vh96_shear_correction)
                                ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
                                : amrex::max(ws2_raw, one_d);
                Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
            }
            Real zval0 = zval, Rib0 = Rib;
 
            bool above_critical = false;
            while (!above_critical && ((kpbl + 1) <= khi)) {
                zval0 = zval;
                Rib0  = Rib;
                kpbl += 1;
 
                zval = (use_terrain_fitted_coords)
                     ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
                     : ((kpbl + myhalf) * gdata.CellSize(2) - zib);
                const Real theta_v     = GetThetav(i, j, kpbl, cell_data, moisture_indices);
                const Real theta_v_klo = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), one_d);
                const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
                                                (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
                                                (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
                // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
                const Real ws2 = (turbChoice.enable_vh96_shear_correction)
                                ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
                                : amrex::max(ws2_raw, one_d);
                Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
                above_critical = (Rib >= Ribcr);
            }
 
            const Real pblh_emp = (use_terrain_fitted_coords)
                                ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
                                : myhalf * gdata.CellSize(2);
            const Real z_max = (use_terrain_fitted_coords)
                             ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr)
                             : ((khi + myhalf) * gdata.CellSize(2) - zib);
            const Real pblh_max = Real(0.9) * z_max;
            const Real pblh_min = amrex::max(pblh_emp, Real(10.0));
 
            if (above_critical) {
                Real pblh_interp;
                const Real rib_diff = Rib - Rib0;
                if (std::abs(rib_diff) > Real(1.0e-10)) {
                    pblh_interp = zval0 + (zval - zval0) / rib_diff * (Ribcr - Rib0);
                } else {
                    // If Rib not changing, use zval0 as PBL height
                    pblh_interp = zval0;
                }
                pblh_corr_arr(i, j, 0) = amrex::max(amrex::min(pblh_interp, pblh_max), pblh_min);
                pbli_arr(i, j, 0) = kpbl;
            } else {
                pblh_corr_arr(i, j, 0) = pblh_min;
                pbli_arr(i, j, 0) = ksrf + 1;
            }
            pblh_floor_arr(i, j, 0) = pblh_min;
        });
 
 
        // Apply PBLH spatial smoothing if enabled (Seibert et al. 2000 methodology)
        // Seibert et al. (2000): Review and intercomparison of operational methods
        // for the determination of the mixing height. Atmospheric Environment, 34, 1001-1027.
        // Spatial smoothing removes unphysical grid-to-grid noise from discrete Rib-crossing detection
        if (turbChoice.enable_pblh_smoothing) {
            // Smooth the tile's own columns, reading the halo the passes above
            // filled. The result is the same however the domain is split.
            // The stencil mixes neighbouring columns, so it runs on the
            // absolute height (the stored convention): the column's surface
            // height zib goes on over the work box first and comes off
            // after, and the result is kept at or above the column's floor.
            // zib is zero without erf.pbl_ib_aware.
            ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept {
                pblh_corr_arr(i, j, 0) += zib_arr(i, j, 0);
            });
            ApplyPBLHSmoothing(pbl_height_corrector, xybx_tile,
                             turbChoice.pblh_smoothing_weight,
                             turbChoice.pblh_smoothing_passes,
                             geom.Domain(), geom.periodicity());
            ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept {
                pblh_corr_arr(i, j, 0) = amrex::max(pblh_corr_arr(i, j, 0) - zib_arr(i, j, 0), pblh_floor_arr(i, j, 0));
            });
        }
 
        // Copy corrected PBL height into pblh_mf for SurfaceLayer storage.
        // pbl_height_corrector only covers this tile (grown by one), so loop over
        // the tile, not the valid box: with tiling in x/y the valid box reaches
        // past it. Under TileNoZ the tile spans the full column, and pblh_mf has
        // no ghost cells, so the tiles together still fill every cell.
        {
            auto pblh_out = pblh_mf.array(mfi);
            const Box& tbx = mfi.tilebox();
            // The scheme measures its height from the top of an immersed column
            // (zib); the stored field is the absolute height above the domain
            // bottom, the convention of YSUNew and of every consumer of
            // SurfaceLayer::get_pblh(). zib is zero without erf.pbl_ib_aware.
            ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                pblh_out(i, j, k) = pblh_corr_arr(i, j, 0) + zib_arr(i, j, 0);
            });
        }
        //
        // PASS 4 (WSTAR RECOMPUTE): Recompute wstar, HGAMT, HGAMQ using the corrected
        // PBL height (pblh_corr_arr) to ensure internal consistency between the K-profile
        // amplitude and the countergradient fluxes. HOL = sf*h/L uses pblh_corr_arr.
        // WRF reference (module_bl_mrf.F lines 857-880):
        // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L857-L880
        //
        ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            const int  ksrf = ksurf_arr(i, j, 0);   // local surface: first fluid cell of the column
            const Real zib  = zib_arr(i, j, 0);     // its height, subtracted from every height below
            amrex::ignore_unused(zib);
            //const Real t_layer  = t10_eff_arr(i, j, 0);
            Real obuk_val = ol_eff_arr(i, j, 0);
            if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
                obuk_val = (obuk_val >= Real(0)) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
            }
 
            // HOL now uses corrected PBLH for full internal consistency
            const Real HOL = sf * pblh_corr_arr(i, j, 0) / obuk_val;
            const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
            // |HOL|, so the base 1 + 16*HOL_abs of the unstable arm below equals its
            // 1 - 16*HOL there and is at least 1 for every HOL; a guard written inside the
            // arm gets folded away by the optimiser and the pow hoisted (see sqrt_neg_Ri).
            const Real HOL_abs = std::abs(HOL_bounded);
            const Real one_quarter = Real(0.25);
            // Enable QNSE stable functions if requested, otherwise use default linear form
            const Real enable_qnse_d = (turbChoice.enable_qnse_stable_functions) ? Real(1.0) : Real(0.0);
            const Real qnse_am_d = turbChoice.qnse_am;
            const Real qnse_bm_d = turbChoice.qnse_bm;
            const Real phiM = (obuk_val > 0)
                            ? (enable_qnse_d > Real(0.5)
                               ? (1 + qnse_am_d * HOL_bounded) / (1 + qnse_bm_d * HOL_bounded)
                               : (1 + 5 * HOL_bounded))
                            : std::pow(1 + 16 * HOL_abs, -one_quarter);
            const Real phiM_safe = amrex::max(phiM, Real(0.01));
 
            // Absolute bounds [0.01, 5.0] m/s
            // Deardorff (1972) convective velocity scale — well-behaved as u_star -> 0.
            // wstar_conv = (g/theta_v0 * zi * wthv0)^(1/3)
            const Real theta_v0 = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), Real(1.0));
            const Real wthv0 = amrex::max(-us_eff_arr(i,j,0) * ts_eff_arr(i,j,0), Real(0.0)); // surface buoyancy flux, unstable only
            // Corrected height: this pass exists to make the K-profile amplitude and the
            // countergradient fluxes consistent with the height Pass 3 produced, and HOL
            // above already uses it.  (Pass 2 has only the predictor height to work with.)
            const Real zi = amrex::max(pblh_corr_arr(i, j, 0), Real(10.0));
            const Real wstar_conv = std::cbrt(CONST_GRAV / theta_v0 * zi * wthv0);
 
            Real wstar_shear = us_eff_arr(i, j, 0) / phiM_safe;
 
            // Blend: standard convective/shear combination (WRF YSU-style),
            // weighted by cube sum so either term can dominate smoothly
            Real wstar = std::cbrt(wstar_shear*wstar_shear*wstar_shear +
                                wstar_conv*wstar_conv*wstar_conv);
            //Real wstar = us_eff_arr(i, j, 0) / phiM_safe;
            wstar = amrex::max(wstar, Real(0.01));
            wstar = amrex::min(wstar, Real(5.0));
            wstar_arr(i, j, 0) = wstar;
            bool SFCFLG = (obuk_val <= Real(0));
            // WRF caps the heat countergradient at GAMCRT and then floors it at zero:
            //   HGAMT(I)=MIN(GAMFAC*HFX(I)/CPM(I),GAMCRT)   (module_bl_mrf.F L874)
            //   HGAMT(I)=MAX(HGAMT(I),0.0)                  (module_bl_mrf.F L880)
            // The floor is not redundant here.  The unstable branch is selected on
            // sign(L), but the numerator is the sensible flux -u_* theta_* alone: over
            // water with a small downward sensible flux and a large latent flux, L is
            // negative while -u_* theta_* is not.  Without the floor the term would then
            // add to the down-gradient flux rather than oppose it.  VPERT in Pass 2 is
            // deliberately built from the unfloored value, as in WRF (L877 precedes L880).
            const Real HGAMT = (SFCFLG && enable_mrf_countergradient)
                             ? amrex::max(amrex::min(-const_b * us_eff_arr(i, j, 0)
                                                     * ts_eff_arr(i, j, 0) / wstar, GAMCRT), Real(0))
                             : Real(0);
 
            Real HGAMQ = Real(0);
            if (SFCFLG && use_moisture && enable_mrf_countergradient) {
                const Real q_star = qs_eff_arr(i, j, 0);
                const Real HGAMQ_calc = -const_b * us_eff_arr(i, j, 0) * q_star / wstar;
                HGAMQ = amrex::max(amrex::min(HGAMQ_calc, GAMCRQ), Real(0));
 
                if (lmask_arr) {
                    bool is_land = (lmask_arr(i,j,0) == 1);
                    if (!is_land) HGAMQ = Real(0);
                }
 
                if (moisture_indices.qv >= 0) {
                    Real qv_klo = cell_data(i, j, ksrf, moisture_indices.qv) / cell_data(i, j, ksrf, Rho_comp);
                    Real T_klo = getTgivenRandRTh(cell_data(i, j, ksrf, Rho_comp),
                                                  cell_data(i, j, ksrf, RhoTheta_comp), qv_klo);
                    Real p_klo = getPgivenRTh(cell_data(i, j, ksrf, RhoTheta_comp), qv_klo) * Real(0.01);
                    Real qsat_klo = Real(0);
                    erf_qsatw(T_klo, p_klo, qsat_klo);
                    Real rh_klo = (qsat_klo > Real(1.0e-10)) ? (qv_klo / qsat_klo) : Real(0);
                    if (rh_klo > Real(0.95)) {
                        Real rh_scaling = amrex::max(Real(0), (Real(1) - rh_klo) / Real(0.05));
                        HGAMQ *= rh_scaling;
                    }
                }
            }
 
            if (pbli_arr(i, j, 0) <= ksrf + 1) {
                hgamt_arr(i, j, 0) = Real(0);
                hgamq_arr(i, j, 0) = Real(0);
            } else {
                const Real pblh = pblh_corr_arr(i, j, 0);
                if (pblh > Real(1.0e-10)) {
                    hgamt_arr(i, j, 0) = (enable_mrf_countergradient) ? HGAMT / pblh : Real(0);
                    hgamq_arr(i, j, 0) = (enable_mrf_countergradient && use_moisture) ? HGAMQ / pblh : Real(0);
                } else {
                    hgamt_arr(i, j, 0) = Real(0);
                    hgamq_arr(i, j, 0) = Real(0);
                }
            }
        });
        //
        // PASS 5 (ZERO-RI): Diagnostic PBL height with Ribcr=0, VPERT-enhanced surface temp.
        // Used optionally (pbl_mrf_use_zero_ri_extent) to extend the nonlocal mixing region.
        // WRF reference (module_bl_mrf.F lines 932-964):
        // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L932-L964
        //
        constexpr Real Ribcr_zero = Real(0);
        ParallelFor(xybx, [=,one_d=one]
                    AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            const int  ksrf = ksurf_arr(i, j, 0);   // local surface: first fluid cell of the column
            const Real zib  = zib_arr(i, j, 0);     // its height, subtracted from every height below
            amrex::ignore_unused(zib);
            const Real t_layer  = t10_eff_arr(i, j, 0);
            const Real moisture_fraction = use_moisture ? q10_eff_arr(i, j, 0) : zero;
            const Real t_layer_v = t_layer * (one + epsv * moisture_fraction);
            const Real t_layer_v_enhanced = t_layer_v + vpert_arr(i, j, 0);
 
            int kpbl_zero = ksrf;
            Real zval_zero, Rib_zero;
            {
                zval_zero = (use_terrain_fitted_coords)
                          ? Compute_Zrel_AtCellCenter(i, j, kpbl_zero, z_nd_arr)
                          : ((kpbl_zero + myhalf) * gdata.CellSize(2) - zib);
                const Real theta_v     = GetThetav(i, j, kpbl_zero, cell_data, moisture_indices);
                const Real theta_v_klo = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), one_d);
                const Real ws2_raw = fourth * ( (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) *
                                                (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) +
                                                (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) *
                                                (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) );
                // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
                const Real ws2 = (turbChoice.enable_vh96_shear_correction)
                                ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
                                : amrex::max(ws2_raw, one_d);
                Rib_zero = CONST_GRAV * zval_zero * (theta_v - t_layer_v_enhanced) / (ws2 * theta_v_klo);
            }
 
            bool above_critical_zero = false;
            while (!above_critical_zero && ((kpbl_zero + 1) <= khi)) {
                kpbl_zero += 1;
                zval_zero = (use_terrain_fitted_coords)
                          ? Compute_Zrel_AtCellCenter(i, j, kpbl_zero, z_nd_arr)
                          : ((kpbl_zero + myhalf) * gdata.CellSize(2) - zib);
                const Real theta_v     = GetThetav(i, j, kpbl_zero, cell_data, moisture_indices);
                const Real theta_v_klo = amrex::max(GetThetav(i, j, ksrf, cell_data, moisture_indices), one_d);
                const Real ws2_raw = fourth * ( (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) *
                                                (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) +
                                                (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) *
                                                (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) );
                // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
                const Real ws2 = (turbChoice.enable_vh96_shear_correction)
                                ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
                                : amrex::max(ws2_raw, one_d);
                Rib_zero = CONST_GRAV * zval_zero * (theta_v - t_layer_v_enhanced) / (ws2 * theta_v_klo);
                above_critical_zero = (Rib_zero >= Ribcr_zero);
            }
 
            if (above_critical_zero) {
                pbli_zero_arr(i, j, 0) = kpbl_zero;
            } else {
                pbli_zero_arr(i, j, 0) = ksrf + 1;
            }
        });
 
        // -- Compute diffusion coefficients --
 
        const Array4<Real>& K_turb = eddyViscosity.array(mfi);
 
        bool c_ext_dir_on_zlo = ((bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir));
        bool c_ext_dir_on_zhi = ((bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir));
        bool u_ext_dir_on_zlo = ((bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir));
        bool u_ext_dir_on_zhi = ((bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir));
        bool v_ext_dir_on_zlo = ((bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir));
        bool v_ext_dir_on_zhi = ((bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir));
 
        const auto& dxInv = geom.InvCellSizeArray();
        const Real dz_inv = geom.InvCellSize(2);
        const int izmin   = geom.Domain().smallEnd(2);
        const int izmax   = geom.Domain().bigEnd(2);
 
        // Blending parameters captured as scalars for GPU lambda.
        const amrex::Real l_blend_length = turbChoice.pbl_blend_length;
        const amrex::Real l_blend_cs     = turbChoice.pbl_blend_cs;
        const amrex::Real l_blend_cmax   = turbChoice.pbl_blend_c_max;
        //const bool        l_use_smag_ceil= turbChoice.pbl_blend_use_smag;
        // dx: use horizontal spacing at this level (assume dx = dy for regular grids).
        const amrex::Real l_dx = geom.CellSize(0);
 
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            const int  ksrf = ksurf_arr(i, j, 0);   // local surface: first fluid cell of the column
            const Real zib  = zib_arr(i, j, 0);     // its height, subtracted from every height below
            amrex::ignore_unused(zib);
            Real obuk_val = ol_eff_arr(i, j, 0);
            if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
                obuk_val = (obuk_val >= Real(0)) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
            }
 
            const Real zval = (use_terrain_fitted_coords)
                            ? Compute_Zrel_AtCellCenter(i, j, k, z_nd_arr)
                            : ((k + myhalf) * gdata.CellSize(2) - zib);
            const Real rho = cell_data(i, j, k, Rho_comp);
            // Guard: skip lateral ghost cells that may have uninitialized density.
            // Prevents division-by-zero inside GetThetav at lateral ghost cells of gbx.
            // REGRID NOTE: On regrid, SurfaceLayer is reallocated and PBLH returns to sentinel.
            // The driver must call ComputeDiffusivityMRF (or an equivalent bootstrap pass)
            // before consuming PBLH in update_fluxes() or Beljaars correction after any regrid.
            if (rho <= Real(0) || k < ksrf) {   // no diffusivity inside an immersed solid
                K_turb(i, j, k, EddyDiff::Mom_v)   = Real(0);
                K_turb(i, j, k, EddyDiff::Theta_v) = Real(0);
                K_turb(i, j, k, EddyDiff::Q_v)     = Real(0);
                K_turb(i, j, k, EddyDiff::HGAMT_v) = Real(0);
                K_turb(i, j, k, EddyDiff::HGAMQ_v) = Real(0);
                K_turb(i, j, k, EddyDiff::Turb_lengthscale) = Real(0);
                return;
            }
            const Real met_h_zeta = (use_terrain_fitted_coords)
                                  ? Compute_h_zeta_AtCellCenter(i, j, k, dxInv, z_nd_arr) : Real(1);
            const Real dz_terrain = met_h_zeta / dz_inv;
 
            constexpr Real qc_threshold = Real(1.0e-4);  // Cloud water/ice threshold (kg/kg)
            Real qc_mix = Real(0);
            Real qi_mix = Real(0);
            if (use_moisture) {
                if (moisture_indices.qc >= 0) {
                    qc_mix = cell_data(i, j, k, moisture_indices.qc) / rho;
                }
                if (moisture_indices.qi >= 0) {
                    qi_mix = cell_data(i, j, k, moisture_indices.qi) / rho;
                }
            }
            const Real total_qcloud = qc_mix + qi_mix;
            const bool has_cloud = turbChoice.enable_mrf_cloud_adjustment && (total_qcloud > qc_threshold);
 
            const int pbli_extent = turbChoice.pbl_mrf_use_zero_ri_extent ? pbli_zero_arr(i, j, 0) : pbli_arr(i, j, 0);
 
            if (k < pbli_extent) {
                const Real HOL = sf * pblh_corr_arr(i, j, 0) / obuk_val;
                const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
                // |HOL|, so the base 1 + 16*HOL_abs of the unstable arm below equals its
                // 1 - 16*HOL there and is at least 1 for every HOL; a guard written inside the
                // arm gets folded away by the optimiser and the pow hoisted (see sqrt_neg_Ri).
                const Real HOL_abs = std::abs(HOL_bounded);
 
                const Real one_quarter = Real(0.25);
                // Enable QNSE stable functions if requested, otherwise use default linear form.
                const Real enable_qnse_d = (turbChoice.enable_qnse_stable_functions) ? Real(1.0) : Real(0.0);
                const Real qnse_am_d = turbChoice.qnse_am;
                const Real qnse_bm_d = turbChoice.qnse_bm;
                const Real qnse_ah_d = turbChoice.qnse_ah;
                const Real qnse_bh_d = turbChoice.qnse_bh;
                const Real phiM = (obuk_val > 0)
                                ? (enable_qnse_d > Real(0.5)
                                   ? (1 + qnse_am_d * HOL_bounded) / (1 + qnse_bm_d * HOL_bounded)
                                   : (1 + 5 * HOL_bounded))
                                : std::pow(1 + 16 * HOL_abs, -one_quarter);
                const Real phit = (obuk_val > 0)
                                ? (enable_qnse_d > Real(0.5)
                                   ? (1 + qnse_ah_d * HOL_bounded) / (1 + qnse_bh_d * HOL_bounded)
                                   : (1 + 5 * HOL_bounded))
                                : std::pow(1 + 16 * HOL_abs, -Real(0.5));
 
                Real phit_cloud = phit;
                Real phiM_cloud = phiM;
                if (has_cloud && obuk_val > Real(0)) {
                    Real reduction_factor = Real(1) - Real(0.15) * amrex::min(total_qcloud / qc_threshold, Real(1));
                    // Apply cloud reduction on top of QNSE if enabled
                    if (enable_qnse_d > Real(0.5)) {
                        phiM_cloud = ((1 + qnse_am_d * HOL_bounded * reduction_factor) / (1 + qnse_bm_d * HOL_bounded));
                        phit_cloud = ((1 + qnse_ah_d * HOL_bounded * reduction_factor) / (1 + qnse_bh_d * HOL_bounded));
                    } else {
                        phiM_cloud = Real(1) + Real(5.0) * reduction_factor * sf * pblh_corr_arr(i, j, 0) / obuk_val;
                        phit_cloud = Real(1) + Real(5.0) * reduction_factor * sf * pblh_corr_arr(i, j, 0) / obuk_val;
                    }
                } else if (has_cloud && obuk_val <= Real(0)) {
                    Real cloud_boost = Real(1.0) + Real(0.05) * amrex::min(total_qcloud / qc_threshold, Real(1));
                    phiM_cloud = std::pow(Real(1) + Real(16.0) * HOL_abs / cloud_boost, -one_quarter);
                    phit_cloud = std::pow(Real(1) + Real(16.0) * HOL_abs / cloud_boost, -Real(0.5));
                }
 
                const Real phiM_eff = phiM_cloud;
                const Real phit_eff = phit_cloud;
 
                Real Prt_base = phit_eff / phiM_eff;
                const Real Prt = amrex::min(amrex::max(Prt_base + const_b * KAPPA * sf, prmin), prmax);
 
                // wstar is WRF's WSCALE: u_*/phi_m in a stable or neutral column (the
                // convective part is floored at zero by the surface buoyancy flux), and
                // the cube-sum blend with w_* when the surface flux is upward.
                const Real wstar = wstar_arr(i, j, 0);
 
                // K-profile: K = rho * wstar * kappa * z * (1 - (z - zl1)/(h - zl1))^2
                //
                // The first-level offset ZL1 belongs to ZFAC alone: the linear factor
                // is the full height above the surface, ZQ(I,K).  Subtracting zl1 from
                // both drives K_m to exactly zero in the first fluid cell.
                // ERF_ComputeDiffusivityYSUNew.cpp:1699-1720 splits them the same way.
                //
                // zval already measures height from the top of an immersed column (zib
                // was subtracted where it was formed), so it is that height directly;
                // zl1 is zval evaluated at the column's first fluid cell, so zval - zl1
                // is exactly zero there in both the terrain and the non-terrain path.
                const Real zl1 = (use_terrain_fitted_coords)
                               ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
                               : ((ksrf + myhalf) * gdata.CellSize(2) - zib);
                const Real pblh = pblh_corr_arr(i, j, 0);
                // The Rib passes floor pblh at max(zl1, 10), so pblh - zl1 is exactly
                // zero for any column whose height lands on that floor with a first
                // layer thicker than 20 m -- and zval - zl1 is zero too at k == ksrf,
                // making the quotient 0/0.  Guard it, as YSUNew does (line 1701).
                const Real pblh_rel = amrex::max(pblh - zl1, Real(1.0e-4));
                const Real zfac = amrex::max(Real(1) - (zval - zl1) / pblh_rel, Real(1.0e-8));
 
                K_turb(i, j, k, EddyDiff::Mom_v)   = rho * wstar * KAPPA * zval * zfac * zfac;
                K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prt;
                K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
            } else if (k >= pbli_extent) {
                // Free atmosphere above the PBL top.  This is the only branch left that
                // reads resolved vertical gradients, but they are not evaluated anywhere
                // near the surface.
                const Real lambda = Real(150.0);
                const Real lscale = (KAPPA * zval * lambda) / (KAPPA * zval + lambda);
                Real dthetadz, dudz, dvdz;
                ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
                                              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, pbl_strat);
 
                const Real dudz_safe = (k < izmax) ? dudz : Real(0);
                const Real dvdz_safe = (k < izmax) ? dvdz : Real(0);
                const Real wind_shear = dudz_safe * dudz_safe + dvdz_safe * dvdz_safe;
                const Real wind_shear_safe = std::max(wind_shear, Real(1.0e-8));
                const Real theta_v     = amrex::max(GetThetav(i, j, k,   cell_data, moisture_indices), Real(1.0));
                const Real dtheta_v_dz = dthetadz;
 
                Real grad_Ri = CONST_GRAV / theta_v * dtheta_v_dz / wind_shear_safe;
                grad_Ri = std::max(std::min(grad_Ri, Real(100.0)), -Real(100.0));
 
                const Real grad_Ri_safe = amrex::max(grad_Ri, -Real(100.0));
                // sqrt(|Ri|), which equals the sqrt(-Ri) of the unstable arm below and is valid
                // for every Ri. Written as max(-Ri, 0) inside the arm, the optimiser used the
                // arm's condition to drop the max and then hoisted the bare sqrt above the
                // selection, where Ri can be positive; written as min(Ri, 0) it folded the
                // sqrt into the arms of the min. Either way an invalid operation whose
                // result is discarded but whose flag kills the run under
                // amrex.fpe_trap_invalid. There is nothing to fold in the absolute value.
                const Real sqrt_neg_Ri = std::sqrt(std::abs(grad_Ri_safe));
                Real Pr = Real(1) + Real(2.1) * grad_Ri;
                const Real fm = (grad_Ri_safe > 0)
                              ? Real(1) / ((Real(1) + Real(5.0) * grad_Ri_safe) * (Real(1) + Real(5.0) * grad_Ri_safe))
                              : 1 - 8 * grad_Ri_safe / (1 + Real(1.746) * sqrt_neg_Ri);
                const Real ft = (grad_Ri_safe > 0)
                              ? Real(1) / ((Real(1) + Real(5.0) * grad_Ri_safe) * (Real(1) + Real(5.0) * grad_Ri_safe))
                              : 1 - 8 * grad_Ri_safe / (1 + Real(1.286) * sqrt_neg_Ri);
                const Real rl2wsp = rho * lscale * lscale * std::sqrt(wind_shear);
 
                Pr = std::max(amrex::Real(0.25), std::min(Pr, Real(4.0)));
 
                // In the stable regime, fm and ft are the same function, so the Prandtl
                // number has to be applied explicitly or momentum and heat mix
                // identically.
                //
                // In the unstable regime, fm and ft already differ (1.746 vs 1.286) and
                // WRF applies no Prandtl factor there.
                const Real Pr_mom = (grad_Ri_safe > 0) ? Pr : Real(1);
 
                K_turb(i, j, k, EddyDiff::Mom_v)   = rl2wsp * fm * Pr_mom;
                K_turb(i, j, k, EddyDiff::Theta_v) = rl2wsp * ft;
                K_turb(i, j, k, EddyDiff::Q_v)     = K_turb(i, j, k, EddyDiff::Theta_v);
            }
 
            // Scale-aware blending for grey-zone resolution (Boutle et al. 2014).
            // ERF_PBLScaleAwareBlending.H. Gated by pbl_blend_length > 0.
            // Applied to Theta_v and Q_v. Mom_v is not modified.
            //
            // K_turb holds rho*K [kg/m/s], but the ceiling applied inside
            // pbl_kh_blend_and_cap is a bare diffusivity [m^2/s].  Hand the helper the
            // kinematic value and restore the density weighting on the way out; passing
            // rho*K straight through would make the effective cap scale like 1/rho and
            // hence increase with height (issue #3580).
            if (l_blend_length > 0.0) {
                const Real rho_inv = Real(1.0) / rho;
 
                // For MRF, use power-law ceiling (SmnSmn not available)
                K_turb(i, j, k, EddyDiff::Theta_v) = rho * pbl_kh_blend_and_cap(
                    K_turb(i, j, k, EddyDiff::Theta_v) * rho_inv,
                    l_dx, l_blend_length, l_blend_cs, l_blend_cmax,
                    amrex::Real(-1.0), false);
 
                K_turb(i, j, k, EddyDiff::Q_v) = rho * pbl_kh_blend_and_cap(
                    K_turb(i, j, k, EddyDiff::Q_v) * rho_inv,
                    l_dx, l_blend_length, l_blend_cs, l_blend_cmax,
                    amrex::Real(-1.0), false);
            }
 
            // Limit diffusion coefficients to physical bounds
            // Hong & Pan (1996): Kmin=0.1, Kmax=300 m^2/s (module_bl_mrf.F lines 1014-1025)
            // Hong et al. (2006): Kmin=ckz*dz*rho, Kmax=1000 m^2/s (high-res option)
            Real rhoKmin, rhoKmax;
            if (turbChoice.pbl_mrf_highres_bounds) {
                constexpr Real ckz  = Real(0.001);
                constexpr Real Kmax = Real(1000.0);
                rhoKmin = ckz * dz_terrain * rho;
                rhoKmax = rho * Kmax;
            } else {
                constexpr Real Kmin = Real(0.1);
                constexpr Real Kmax = Real(300.0);
                rhoKmin = rho * Kmin;
                rhoKmax = rho * Kmax;
            }
 
            K_turb(i, j, k, EddyDiff::Mom_v) = std::max(
                std::min(K_turb(i, j, k, EddyDiff::Mom_v), rhoKmax), rhoKmin);
            K_turb(i, j, k, EddyDiff::Theta_v) = std::max(
                std::min(K_turb(i, j, k, EddyDiff::Theta_v), rhoKmax), rhoKmin);
            K_turb(i, j, k, EddyDiff::Q_v) = std::max(
                std::min(K_turb(i, j, k, EddyDiff::Q_v), rhoKmax), rhoKmin);
            K_turb(i, j, k, EddyDiff::Turb_lengthscale) = pblh_corr_arr(i, j, 0) + zib;   // absolute height, as stored
 
            // IMPORTANT — units convention for HGAMT_v / HGAMQ_v:
            // These are stored as (countergradient term / pblh), i.e. already
            // divided by the PBL height, NOT as a raw physical gradient needing
            // an extra 1/dz scaling downstream.
            //
            // The implicit solver (ERF_ImplicitDiff_T.cpp, ImplicitDiffForStateLU_T)
            // multiplies this stored value by `Fact = implicit_fac * dt * dz_inv`
            // and `rhoAlpha / met_h_zeta` ONLY — it must NOT apply an additional
            // dz_inv factor. `Fact` already carries the single grid-spacing
            // factor needed to convert the countergradient correction into a
            // flux-divergence contribution consistent with the K*d(phi)/dz terms
            // it is being added alongside. Adding a second dz_inv
            // double-counts the grid spacing and is
            // dimensionally incorrect — do NOT "fix" this by inserting dz_inv
            // into the RHS_a += Fact * gam_hi/lo terms in ERF_ImplicitDiff_T.cpp.
            if (k < pbli_extent) {
                K_turb(i, j, k, EddyDiff::HGAMT_v) = hgamt_arr(i, j, 0);
                K_turb(i, j, k, EddyDiff::HGAMQ_v) = hgamq_arr(i, j, 0);
            } else {
                K_turb(i, j, k, EddyDiff::HGAMT_v) = Real(0);
                K_turb(i, j, k, EddyDiff::HGAMQ_v) = Real(0);
            }
        });
 
        // FOEXTRAP top and bottom ghost cells
        ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
        {
            K_turb(i, j, klo-1, EddyDiff::Mom_v  ) = K_turb(i, j, klo, EddyDiff::Mom_v  );
            K_turb(i, j, klo-1, EddyDiff::Theta_v) = K_turb(i, j, klo, EddyDiff::Theta_v);
            K_turb(i, j, klo-1, EddyDiff::Q_v    ) = K_turb(i, j, klo, EddyDiff::Q_v    );
            K_turb(i, j, klo-1, EddyDiff::HGAMT_v) = K_turb(i, j, klo, EddyDiff::HGAMT_v);
            K_turb(i, j, klo-1, EddyDiff::HGAMQ_v) = K_turb(i, j, klo, EddyDiff::HGAMQ_v);
            K_turb(i, j, klo-1, EddyDiff::Turb_lengthscale) = K_turb(i, j, klo, EddyDiff::Turb_lengthscale);
            K_turb(i, j, khi+1, EddyDiff::Mom_v  ) = K_turb(i, j, khi, EddyDiff::Mom_v  );
            K_turb(i, j, khi+1, EddyDiff::Theta_v) = K_turb(i, j, khi, EddyDiff::Theta_v);
            K_turb(i, j, khi+1, EddyDiff::Q_v    ) = K_turb(i, j, khi, EddyDiff::Q_v    );
            K_turb(i, j, khi+1, EddyDiff::HGAMT_v) = K_turb(i, j, khi, EddyDiff::HGAMT_v);
            K_turb(i, j, khi+1, EddyDiff::HGAMQ_v) = K_turb(i, j, khi, EddyDiff::HGAMQ_v);
            K_turb(i, j, khi+1, EddyDiff::Turb_lengthscale) = K_turb(i, j, khi, EddyDiff::Turb_lengthscale);
        });
    }// mfi
    // Write MRF-computed PBLH back into SurfaceLayer so Beljaars correction
    // and diagnostics can use it, and update_pblh no longer aborts for MRF type.
    SurfLayer->set_pblh(level, pblh_mf);
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function: