ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeDiffusivityYSUNew.cpp File Reference
#include <algorithm>
#include <AMReX_Math.H>
#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 <cmath>
Include dependency graph for ERF_ComputeDiffusivityYSUNew.cpp:

Functions

void ComputeDiffusivityYSUNew (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 *qheating_rates, const MultiFab *terrain_blank)
 

Function Documentation

◆ ComputeDiffusivityYSUNew()

void ComputeDiffusivityYSUNew ( 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 *  qheating_rates,
const MultiFab *  terrain_blank 
)

Compute vertical eddy viscosity coefficients using the Yonsei University (YSU) boundary layer scheme.

Parameters
[in]xvelX-velocity field.
[in]yvelY-velocity field.
[in]cons_inInput conservative variables.
[out]eddyViscosityMultiFab to store computed eddy viscosity and countergradient terms.
[in]geomGeometry used for grid spacings and domain extent.
[in]turbChoiceTurbulence model configuration and parameters.
[in]SurfLayerPointer to surface layer data.
[in]use_terrain_fitted_coordsUse terrain-fitted coordinates if true.
[in]use_moistureInclude moisture in the diffusivity calculation.
[in]levelCurrent AMR level.
[in]bc_ptrBoundary condition records.
[in]vert_onlyReserved flag for vertical-only computation.
[in]z_phys_ndNodal physical height field.
[in]z_phys_ccCell-centered physical height field.
[in]moisture_indicesIndices for moisture variables in the state vector.
[in]qheating_ratesOptional heating rates for cloud-top mixing.
60 {
61  const StratType pbl_strat = turbChoice.pbl_strat_type;
62  /*
63  ============================================================================
64  Yonsei University (YSU) Boundary Layer Parameterization Scheme
65  ============================================================================
66 
67  Implementation of the YSU boundary layer scheme based on:
68  - Hong, S.-Y., Y. Noh, and J. Dudhia, 2006: A new vertical diffusion
69  package with an explicit treatment of entrainment processes.
70  Monthly Weather Review, 134, 2318-2341. [HND06]
71  https://doi.org/10.1175/MWR3250.1
72  - Hong, S.-Y., 2010: A new stable boundary-layer mixing scheme and
73  its impact on the simulated East Asian summer monsoon.
74  Quarterly Journal of the Royal Meteorological Society, 136, 1481-1496. [H10]
75  https://doi.org/10.1002/qj.665
76  - WRF Reference Implementation: module_bl_ysu.F
77  https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F
78 
79  CORE ALGORITHM (HND06):
80  -----------------------
81  1. Three-Pass Bulk Richardson Number (Rib) PBL Height Diagnosis
82  Pass 1: Ribcr=f(surface type), base surface θ_v (predictor)
83  Pass 2: Same Ribcr, enhanced surface θ_v + VPERT (corrector)
84  Pass 3: Ribcr=0 diagnostic, gives mixed-layer depth h_ze
85 
86  2. Nonlocal Countergradient Flux Corrections
87  HGAMT = -const_b * u_* * θ_* / wscale (capped at GAMCRT=3K)
88  HGAMQ = -const_b * u_* * q_* / wscale (capped at GAMCRQ=2e-3)
89  VPERT = max(HGAMT + 0.61*θ*HGAMQ, 0)
90 
91  3. K-Profile Below PBL (HND06 Eq. 2):
92  K_m = ρ * wscale * κ * z * (1 - z/h)^pfac
93  K_t = K_m / Pr_t
94  wscale = (u_*³ + phifac * κ * wstar_conv³ * (1-zfac))^(1/3)
95 
96  4. Explicit Entrainment at PBL Top (WRF bl_ysu.F90 lines 831-925):
97  K_entr = ρ * |we| * dz at k = kpbl (applied only when we < 0)
98  wm3 = wstar³ + 5*u_*³
99  bfxpbl = -0.15 * θv(klo)/g * wm3 / pblh
100  dthvx = max(θv(kpbl+1) - θv(kpbl), 1e-2)
101  we = max(bfxpbl/dthvx, -sqrt(wm3^(2/3)))
102 
103  5. Free Atmosphere / Stable PBL (H10, YSU Appendix A):
104  Ri_g-dependent mixing with grid-adaptive length scale
105  λ = min(max(0.1*dz, 30m), 300m)
106  l = λ * κz / (λ + κz)
107  Optional QNSE stability functions (Sukoriansky et al. 2005):
108  - If enable_qnse_stable_functions=true, replace linear stable stability
109  functions with bounded rational QNSE functions:
110  phi_m = (1 + qnse_am * zeta) / (1 + qnse_bm * zeta)
111  phi_h = (1 + qnse_ah * zeta) / (1 + qnse_bh * zeta)
112  where zeta = sf * h/L, avoiding unbounded growth at large stability
113 
114  PARAMETER DEFAULTS (Matching WRF-YSU):
115  ----------------------------------------
116  const_b = 7.8 (CFAC: countergradient coefficient)
117  sf = 0.1 (SFCFRAC: surface layer depth fraction)
118  phifac = 8.0 (PHIFAC: convective wscale weight)
119  pfac = 2.0 (profile shape exponent)
120  GAMCRT = 3.0 K (max heat countergradient)
121  GAMCRQ = 2e-3 (max moisture countergradient)
122  Ribcr land= 0.25 (critical bulk Richardson number over land)
123  QNSE default coefficients: qnse_am/qnse_bm/qnse_ah/qnse_bh = 2.5/0.2/2.5/0.2
124  */
125 
126  // Domain extent in z-dir
127  int klo = geom.Domain().smallEnd(2);
128  int khi = geom.Domain().bigEnd(2);
129  const int izmin = klo;
130  const int izmax = khi;
131 
132  const Real dz = geom.CellSize(2);
133  const Real dz_inv = geom.InvCellSize(2);
134  const auto& dxInv = geom.InvCellSizeArray();
135 
136  // Collect YSUNew-computed PBLH for writing back to SurfaceLayer.
137  // NOTE: OpenMP is not used here; if re-enabled, this MultiFab must remain
138  // outside any omp parallel structured block to avoid undefined behaviour.
139  MultiFab pblh_mf(eddyViscosity.boxArray(), eddyViscosity.DistributionMap(), 1, 0);
140  pblh_mf.setVal(0.0);
141 
142  // PBLH smoothing reads one column per pass outside the cells it writes, so the
143  // work arrays have to carry that many columns of halo and the PBLH passes have
144  // to fill them: otherwise the stencil reads off the end of the tile and the
145  // answer depends on the decomposition. With smoothing off this is the single
146  // column of halo the diffusivity kernels already use, so nothing changes.
147  const int ng_pblh = (turbChoice.enable_pblh_smoothing)
148  ? std::max(1, turbChoice.pblh_smoothing_passes) : 1;
149  if (ng_pblh > 1) {
150  const int ng_avail = std::min({cons_in.nGrowVect()[0], cons_in.nGrowVect()[1],
151  xvel.nGrowVect()[0], xvel.nGrowVect()[1],
152  yvel.nGrowVect()[0], yvel.nGrowVect()[1],
153  SurfLayer->get_u_star(level)->nGrowVect()[0],
154  SurfLayer->get_u_star(level)->nGrowVect()[1],
155  SurfLayer->get_olen(level)->nGrowVect()[0],
156  SurfLayer->get_olen(level)->nGrowVect()[1],
157  (terrain_blank) ? terrain_blank->nGrowVect()[0] : 1000,
158  (terrain_blank) ? terrain_blank->nGrowVect()[1] : 1000});
159  if (ng_pblh > ng_avail) {
160  amrex::Abort("erf.pblh_smoothing_passes = " + std::to_string(turbChoice.pblh_smoothing_passes)
161  + " needs " + std::to_string(ng_pblh) + " halo columns, but the state and "
162  "surface-layer arrays carry only " + std::to_string(ng_avail)
163  + "; reduce erf.pblh_smoothing_passes to at most " + std::to_string(ng_avail));
164  }
165  }
166 
167 #ifdef _OPENMP
168 #pragma omp parallel if (Gpu::notInLaunchRegion())
169 #endif
170 
171  for (MFIter mfi(eddyViscosity, TileNoZ()); mfi.isValid(); ++mfi) {
172 
173  // Box operated on must span fill domain in z-dir
174  const Box& gbx = mfi.growntilebox(IntVect(1,1,0));
175  AMREX_ALWAYS_ASSERT( gbx.smallEnd(2) == klo &&
176  gbx.bigEnd(2) == khi );
177 
178  //
179  // PREDICTOR: Compute the height of the PBL without thermal excess
180  // From Hong et al. 2006, Eqns. 1 & 2:
181  //
182  // h = Rib_cf * theta_va * | U(h) |^2 / (g * (theta_v(h) - theta_s))
183  //
184  // where the surface virtual potential temperature is
185  //
186  // theta_s = theta_va + theta_T
187  //
188  // and here the thermal excess theta_T = zero
189  //
190  // WRF reference (module_bl_ysu.F lines 100-130):
191  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F#L100-L130
192 
193  // create flattened boxes to store PBL height and related quantities
194  const Box xybx = PerpendicularBox<ZDir>(gbx, IntVect{0, 0, 0});
195 
196  // The PBLH passes run on the region the smoothing stencil will consume
197  // (see ng_pblh above). growntilebox() is no use here: it does not grow at
198  // an interior tile edge, so a tile in the middle of a box gets no halo at
199  // all -- which is exactly how the stencil came to read off the end of the
200  // array. Grow the tile box explicitly instead, which costs a little
201  // duplicated work in the overlaps. With smoothing off these are gbx and
202  // xybx and nothing is duplicated.
203  const Box gbx_work = (turbChoice.enable_pblh_smoothing)
204  ? amrex::grow(mfi.tilebox(), IntVect(ng_pblh,ng_pblh,0)) : gbx;
205  const Box xybx_work = PerpendicularBox<ZDir>(gbx_work, IntVect{0, 0, 0});
206  const Box xybx_tile = PerpendicularBox<ZDir>(mfi.tilebox(), IntVect{0, 0, 0});
207  FArrayBox pbl_height_corrector(xybx_work, 1, The_Async_Arena());
208  IArrayBox pbl_index(xybx_work, 1, The_Async_Arena());
209  IArrayBox pbl_index_zero_ri(xybx_work, 1, The_Async_Arena()); // Index for zero-Ri diagnostic pass
210  FArrayBox hgamt_fab(xybx_work, 1, The_Async_Arena()); // Store HGAMT/h (normalized countergradient)
211  FArrayBox hgamq_fab(xybx_work, 1, The_Async_Arena()); // Store HGAMQ/h (normalized countergradient)
212  FArrayBox hgamu_fab(xybx_work, 1, The_Async_Arena()); // Store HGAMU (YSU u countergradient)
213  FArrayBox hgamv_fab(xybx_work, 1, The_Async_Arena()); // Store HGAMV (YSU v countergradient)
214  FArrayBox wstar_fab(xybx_work, 1, The_Async_Arena()); // Convective velocity scale computed with pblh_corr
215  FArrayBox vpert_fab(xybx_work, 1, The_Async_Arena()); // Virtual temperature perturbation VPERT for Pass 3
216  FArrayBox entr_fab(xybx_work, 1, The_Async_Arena()); // Per-column entrainment diffusivity
217  IArrayBox cloud_top_fab(xybx_work, 1, The_Async_Arena()); // Cloud-top cell index for top-down mixing
218  FArrayBox wstar3_down_fab(xybx_work, 1, The_Async_Arena()); // Top-down convective velocity scale cubed
219  FArrayBox sflux_fab(xybx_work, 1, The_Async_Arena()); // Virtual kinematic heat flux (sensible + latent)
220  FArrayBox wstar3_fab(xybx_work, 1, The_Async_Arena()); // Convective velocity scale cubed per column
221  FArrayBox zol1_fab(xybx_work, 1, The_Async_Arena()); // Monin-Obukhov stability parameter at first level
222  FArrayBox sfcflg_fab(xybx_work, 1, The_Async_Arena()); // Surface stability flag: 1=unstable/neutral, 0=stable
223  // REQUIRED: zero-initialize all FArrayBoxes before GPU kernels read them.
224  // The_Async_Arena() does not zero memory; reading uninitialized data
225  // corrupts rib_enhan_arr (via vpert_arr) and SFCFLG (via sflux_arr).
226  vpert_fab.setVal<RunOn::Device>(zero);
227  sflux_fab.setVal<RunOn::Device>(zero);
228  wstar3_fab.setVal<RunOn::Device>(zero);
229  wstar3_down_fab.setVal<RunOn::Device>(zero);
230  sfcflg_fab.setVal<RunOn::Device>(zero);
231  hgamt_fab.setVal<RunOn::Device>(zero);
232  hgamq_fab.setVal<RunOn::Device>(zero);
233  hgamu_fab.setVal<RunOn::Device>(zero);
234  hgamv_fab.setVal<RunOn::Device>(zero);
235  wstar_fab.setVal<RunOn::Device>(zero);
236  entr_fab.setVal<RunOn::Device>(zero);
237  pbl_height_corrector.setVal<RunOn::Device>(zero);
238  const auto& pblh_corr_arr = pbl_height_corrector.array();
239  const auto& pbli_arr = pbl_index.array();
240  const auto& pbli_zero_arr = pbl_index_zero_ri.array(); // Zero-Ri diagnostic PBL index
241  const auto& hgamt_arr = hgamt_fab.array();
242  const auto& hgamq_arr = hgamq_fab.array();
243  const auto& hgamu_arr = hgamu_fab.array();
244  const auto& hgamv_arr = hgamv_fab.array();
245  const auto& wstar_arr = wstar_fab.array(); // Stored convective velocity for use in K-profile
246  const auto& vpert_arr = vpert_fab.array(); // Stored VPERT for Pass 3 diagnostic loop
247  const auto& entr_arr = entr_fab.array(); // Stored entrainment diffusivity
248  const auto& cloud_top_arr = cloud_top_fab.array(); // Cloud-top cell index for top-down mixing
249  const auto& wstar3_down_arr = wstar3_down_fab.array(); // Top-down wstar cubed
250  const auto& sflux_arr = sflux_fab.array(); // Virtual kinematic heat flux
251  const auto& wstar3_arr = wstar3_fab.array(); // Convective velocity scale cubed
252  const auto& zol1_arr = zol1_fab.array(); // Monin-Obukhov stability parameter
253  const auto& sfcflg_arr = sfcflg_fab.array(); // Surface stability flag (1=unstable/neutral, 0=stable)
254 
255  // Get some data in arrays
256  const auto& cell_data = cons_in.const_array(mfi);
257  const auto& uvel = xvel.const_array(mfi);
258  const auto& vvel = yvel.const_array(mfi);
259 
260  //const Real Ribcr = turbChoice.pbl_mrf_Ribcr; // Now surface-type-dependent in YSU
261  //const Real f0 = turbChoice.pbl_mrf_coriolis_freq;
262  const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
263  const auto& t_star_arr = SurfLayer->get_t_star(level)->const_array(mfi);
264  const auto& q_star_arr = SurfLayer->get_q_star(level)->const_array(mfi);
265  const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
266  const auto& t10av_arr = SurfLayer->get_mac_avg(level, 3)->const_array(mfi);
267  const auto& q10av_arr = SurfLayer->get_mac_avg(level, 4)->const_array(mfi);
268  const auto& ws10av_arr = SurfLayer->get_mac_avg(level, 6)->const_array(mfi); // 10m wind speed for Rossby number
269  const auto& z0_arr = SurfLayer->get_z0(level)->const_array(mfi); // Roughness length for Rossby number
270  //const auto& t_surf_arr = SurfLayer->get_t_surf(level)->const_array(mfi);
271  // Get land/water mask for proper handling of moisture countergradient
272  const auto& lmask_arr = (SurfLayer->get_lmask(level)) ?
273  SurfLayer->get_lmask(level)->const_array(mfi) :
274  Array4<int>{};
275 
276  // --------------------------------------------------------------------
277  // Immersed-boundary awareness (turbChoice.pbl_ib_aware), as in the MRF
278  // scheme: per column the first fluid cell above the immersed solid
279  // (ksurf) and its height (zib); the surface height z_sfc of every
280  // kernel below becomes zib, indices at the surface become ksurf, the
281  // diffusivities vanish inside the solid, and the surface scales of a
282  // column with solid cells are a neutral log law at its top. With the
283  // option off, or without immersed cells, everything reduces to the
284  // original arithmetic.
285  // --------------------------------------------------------------------
286  const bool l_ib = turbChoice.pbl_ib_aware && (terrain_blank != nullptr);
287  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(l_ib && use_terrain_fitted_coords),
288  "erf.pbl_ib_aware is not supported with terrain-fitted coordinates");
289  const Array4<Real const> blank_arr = l_ib ? terrain_blank->const_array(mfi) : Array4<Real const>{};
290  IArrayBox ksurf_fab(xybx_work, 1, The_Async_Arena());
291  FArrayBox zib_fab(xybx_work, 1, The_Async_Arena());
292  FArrayBox pblh_floor_fab(xybx_work, 1, The_Async_Arena()); // floor of the corrected height per column, reapplied after the smoothing
293  FArrayBox us_eff_fab(xybx_work, 1, The_Async_Arena()), ts_eff_fab(xybx_work, 1, The_Async_Arena());
294  FArrayBox qs_eff_fab(xybx_work, 1, The_Async_Arena()), ol_eff_fab(xybx_work, 1, The_Async_Arena());
295  FArrayBox t10_eff_fab(xybx_work, 1, The_Async_Arena()), q10_eff_fab(xybx_work, 1, The_Async_Arena());
296  FArrayBox ws10_eff_fab(xybx_work, 1, The_Async_Arena()), z0_eff_fab(xybx_work, 1, The_Async_Arena());
297  const auto& ksurf_arr = ksurf_fab.array();
298  const auto& zib_arr = zib_fab.array();
299  const auto& pblh_floor_arr = pblh_floor_fab.array();
300  const auto& us_eff_arr = us_eff_fab.array(); const auto& ts_eff_arr = ts_eff_fab.array();
301  const auto& qs_eff_arr = qs_eff_fab.array(); const auto& ol_eff_arr = ol_eff_fab.array();
302  const auto& t10_eff_arr = t10_eff_fab.array(); const auto& q10_eff_arr = q10_eff_fab.array();
303  const auto& ws10_eff_arr = ws10_eff_fab.array(); const auto& z0_eff_arr = z0_eff_fab.array();
304  {
305  const Real dz_ib = geom.CellSize(2);
306  const Real z0_ib = turbChoice.pbl_ib_z0;
307  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
308  {
309  int ks = klo;
310  if (l_ib) {
311  for (int kk = klo; kk <= khi; ++kk) { if (blank_arr(i, j, kk) >= Real(0.5)) { ks = kk + 1; } }
312  if (ks > khi) { ks = khi; }
313  }
314  ksurf_arr(i, j, 0) = ks;
315  zib_arr(i, j, 0) = (ks - klo) * dz_ib;
316  if (ks > klo) {
317  const Real u_top = myhalf * (uvel(i, j, ks) + uvel(i + 1, j, ks));
318  const Real v_top = myhalf * (vvel(i, j, ks) + vvel(i, j + 1, ks));
319  const Real ws_top = std::sqrt(u_top * u_top + v_top * v_top);
320  us_eff_arr(i, j, 0) = amrex::max(KAPPA * ws_top / std::log(myhalf * dz_ib / z0_ib), Real(1.0e-3));
321  ts_eff_arr(i, j, 0) = Real(0);
322  qs_eff_arr(i, j, 0) = Real(0);
323  ol_eff_arr(i, j, 0) = Real(1.0e10);
324  t10_eff_arr(i, j, 0) = cell_data(i, j, ks, RhoTheta_comp) / cell_data(i, j, ks, Rho_comp);
325  q10_eff_arr(i, j, 0) = (use_moisture && moisture_indices.qv >= 0)
326  ? cell_data(i, j, ks, moisture_indices.qv) / cell_data(i, j, ks, Rho_comp) : Real(0);
327  ws10_eff_arr(i, j, 0) = ws_top;
328  z0_eff_arr(i, j, 0) = z0_ib;
329  } else {
330  us_eff_arr(i, j, 0) = u_star_arr(i, j, 0);
331  ts_eff_arr(i, j, 0) = t_star_arr(i, j, 0);
332  qs_eff_arr(i, j, 0) = q_star_arr(i, j, 0);
333  ol_eff_arr(i, j, 0) = l_obuk_arr(i, j, 0);
334  t10_eff_arr(i, j, 0) = t10av_arr(i, j, 0);
335  q10_eff_arr(i, j, 0) = q10av_arr(i, j, 0);
336  ws10_eff_arr(i, j, 0) = ws10av_arr(i, j, 0);
337  z0_eff_arr(i, j, 0) = z0_arr(i, j, 0);
338  }
339  });
340  }
341  const Array4<Real const> z_nd_arr = z_phys_nd->array(mfi);
342  const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
343  // Get qheating_rates if provided (for LW radiation coupling to top-down mixing)
344  const Array4<Real const> qheat_arr = (qheating_rates != nullptr)
345  ? qheating_rates->const_array(mfi)
346  : Array4<Real const>{};
347  const bool has_qheating_rates = (qheating_rates != nullptr);
348 
349 
350  // ========================================================================
351  // PRE-COMPUTE Rib ARRAYS FOR GPU PERFORMANCE
352  // ========================================================================
353  // Pre-compute both base (without VPERT) and enhanced (with VPERT) Richardson numbers
354  // in a single vectorized kernel to reduce divergent control flow in subsequent PBLH loops.
355  // This improves GPU performance by eliminating serial while-loops in favor of simple
356  // array lookups during the three PBLH pass loops.
357  //
358  const bool enable_ysu_liquid_theta = turbChoice.enable_ysu_liquid_theta;
359  FArrayBox rib_base_fab(gbx_work, 1, The_Async_Arena()); // Rib with base t_layer_v
360  FArrayBox rib_enhan_fab(gbx_work, 1, The_Async_Arena()); // Rib with enhanced t_layer_v+VPERT
361  const auto& rib_base_arr = rib_base_fab.array();
362  const auto& rib_enhan_arr = rib_enhan_fab.array();
363 
364  BL_PROFILE_VAR("YSUNew_Rib_Precompute", prof_rib_precomp);
365  ParallelFor(gbx_work, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
366  {
367  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
368  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
369  const amrex::Real t_layer = t10_eff_arr(i, j, 0);
370  const amrex::Real q_frac = use_moisture ? q10_eff_arr(i, j, 0) : zero;
371  const amrex::Real t_layer_v = t_layer * (one + epsv * q_frac);
372  // On first call, vpert_arr is initialized to zero
373  const amrex::Real t_enh = t_layer_v + vpert_arr(i,j,0);
374  const amrex::Real z_sfc = (use_terrain_fitted_coords)
375  ? Compute_Zrel_AtCellCenter(i,j,ksrf,z_nd_arr) : zib;
376  const amrex::Real zval = (use_terrain_fitted_coords)
377  ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
378  : (k + myhalf) * dz;
379  const amrex::Real zrel = amrex::max(zval - z_sfc, amrex::Real(1.0e-4));
380  const amrex::Real theta_v = (use_moisture && enable_ysu_liquid_theta)
381  ? GetThetavl(i,j,k,cell_data,moisture_indices)
382  : GetThetav(i,j,k,cell_data,moisture_indices);
383  const amrex::Real theta_v_klo = (use_moisture && enable_ysu_liquid_theta)
384  ? GetThetavl(i,j,ksrf,cell_data,moisture_indices)
385  : GetThetav(i,j,ksrf,cell_data,moisture_indices);
386  const amrex::Real ws2_raw = fourth * ((uvel(i,j,k)+uvel(i+1,j,k))*(uvel(i,j,k)+uvel(i+1,j,k))
387  + (vvel(i,j,k)+vvel(i,j+1,k))*(vvel(i,j,k)+vvel(i,j+1,k)));
388  // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
389  // to better represent shear associated with surface-layer turbulence at low wind speeds.
390  // Reference: Vogelezang, D.H.P., and A.A.M. Holtslag, 1996: Evaluation and model impacts of
391  // alternative boundary-layer height formulations. Boundary-Layer Meteorology, 81, 245–269.
392  const amrex::Real ws2 = (turbChoice.enable_vh96_shear_correction)
393  ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
394  : amrex::max(ws2_raw, amrex::Real(1.0));
395  rib_base_arr(i,j,k) = CONST_GRAV * zrel * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
396  rib_enhan_arr(i,j,k) = CONST_GRAV * zrel * (theta_v - t_enh) / (ws2 * theta_v_klo);
397  });
398  BL_PROFILE_VAR_STOP(prof_rib_precomp);
399 
400  // ========================================================================
401  // PRE-COMPUTE surface flux (sflux), stability flag (sfcflg), and convective
402  // velocity scale (wstar3) per column, following WRF-YSU formulation
403  // ========================================================================
404  // WRF Reference: module_bl_ysu.F lines 610-680
405  // These quantities are needed for:
406  // - Stability-dependent functions throughout the PBLH passes
407  // - Countergradient flux calculations (HGAMT, HGAMQ)
408  // - K-profile computations
409  // Computing them once per column improves GPU performance and ensures consistency.
410  //
411  BL_PROFILE_VAR("YSUNew_SurfFlux_Precompute", prof_sflux_precomp);
412  ParallelFor(xybx_work, [=, zero_d=zero] AMREX_GPU_DEVICE(int i, int j, int) noexcept
413  {
414  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
415  const Real rho_sfc = cell_data(i, j, ksrf, Rho_comp);
416  const Real t_layer = t10_eff_arr(i, j, 0); // Dry potential temperature at z1
417 
418  // WRF: rhox = psfc/(rd*tx(1)*tvcon), approximate with cell-bottom density
419  const Real rhox = rho_sfc;
420 
421  // Virtual temperature factor: ep1 = Rv/Rd - 1 = 0.608
422  constexpr Real ep1 = amrex::Real(0.608);
423  constexpr Real cp_air = amrex::Real(1004.0);
424 
425  // WRF bl_ysu.F90 line 613: sflux = hfx/(rho*cp) + qfx/rho*ep1*theta
426  // Approximate hfx and qfx from MOST quantities:
427  // hfx [W/m^2] = -rho * cp * u_star * t_star
428  // qfx [kg/m^2/s] = -rho * u_star * q_star
429  const Real ustar = us_eff_arr(i, j, 0);
430  const Real tstar = ts_eff_arr(i, j, 0);
431  const Real qstar = use_moisture ? qs_eff_arr(i, j, 0) : zero;
432  const Real hfx = -rhox * cp_air * ustar * tstar; // W/m^2
433  const Real qfx = -rhox * ustar * qstar; // kg/m^2/s
434  const Real sflux = hfx / (rhox * cp_air) + qfx / rhox * ep1 * t_layer;
435  sflux_arr(i, j, 0) = sflux;
436 
437  // WRF sfcflg: TRUE when unstable/neutral (sflux > 0), FALSE when stable
438  // WRF bl_ysu.F90 line 615: if(br(i).gt.0.0) sfcflg(i) = .false.
439  // In ERF, br > 0 corresponds to sflux < 0 (stable), so:
440  const bool sfcflg = (sflux > zero);
441  sfcflg_arr(i, j, 0) = sfcflg ? one : zero; // Store as scalar field for unified use
442 
443  // WRF: govrth = g/thx(1) — uses DRY theta at first level
444  const Real govrth = CONST_GRAV / t_layer;
445 
446  // WRF bl_ysu.F90 lines 664-668: wstar3 from sflux
447  // First-guess pblh needed for wstar3 — use nominal 500m for initial estimate
448  // wstar3 is refined in subsequent loop after corrector pblh is available
449  const Real pblh_guess = amrex::Real(500.0);
450  const Real bfx0 = amrex::max(sflux, zero_d);
451  const Real wstar3 = govrth * bfx0 * pblh_guess;
452  wstar3_arr(i, j, 0) = wstar3;
453 
454  // WRF bl_ysu.F90 line 676: wscale = (ust3 + phifac*karman*wstar3*0.5)^(1/3)
455  // The 0.5 factor is the column-average of (1-zfac) over the PBL depth
456  constexpr Real phifac = amrex::Real(8.0);
457  const Real ust3 = ustar * ustar * ustar;
458  Real wscale = std::cbrt(ust3 + phifac * KAPPA * wstar3 * amrex::Real(0.5));
459  // WRF bounds: min(wscale, ust*16), max(wscale, ust/5)
460  wscale = amrex::min(wscale, ustar * amrex::Real(16.0));
461  wscale = amrex::max(wscale, ustar / amrex::Real(5.0));
462  wstar_arr(i, j, 0) = wscale; // Store column-average wscale for HGAMT/HGAMQ
463 
464  // zol1: Monin-Obukhov stability parameter at first level
465  // WRF bl_ysu.F90 lines 651-662: zol1 = max(br*fm*fm/fh, rimin)
466  // Approximate: zol1 = z1 / L_obuk
467  Real obuk_val = ol_eff_arr(i, j, 0);
468  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
469  obuk_val = (obuk_val >= zero) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
470  }
471  const Real zl1 = (use_terrain_fitted_coords)
472  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
473  : (ksrf + myhalf) * dz;
474  Real zol1 = zl1 / obuk_val;
475  zol1 = amrex::max(zol1, amrex::Real(-100.0));
476  if (sfcflg) {
477  zol1 = amrex::min(zol1, amrex::Real(-1.0e-8));
478  } else {
479  zol1 = amrex::max(zol1, amrex::Real(1.0e-8));
480  }
481  zol1_arr(i, j, 0) = zol1;
482  });
483  BL_PROFILE_VAR_STOP(prof_sflux_precomp);
484 
485  // ========================================================================
486  // PASS 1 — PREDICTOR: Compute initial PBL index with surface-type-dependent
487  // critical Richardson number (Ribcr).
488  // ========================================================================
489  // WRF Reference: module_bl_ysu.F (Hong, Noh & Dudhia 2006, lines 1-100)
490  // Ribcr = 0.25 over land; Rossby-number-dependent over water (line 150-170).
491  // Rib = (g*z/θv0) * (θv(z) - θv_surf) / ws² (lines 180-200)
492  // References: Hong et al., Mon. Wea. Rev., 134, 2318-2341 (2006)
493  //
494 
495  BL_PROFILE_VAR("YSUNew_PBLH_Passes", prof_pblh);
496  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
497  {
498  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
499  // Determine surface-type-dependent critical Richardson number
500  // Over land: Ribcr = 0.25; over water: Ribcr depends on Rossby number
501  Real Ribcr;
502  {
503  bool over_land = (!lmask_arr) || (lmask_arr(i, j, 0) == 1);
504  if (over_land) {
505  Ribcr = turbChoice.pbl_ysu_land_Ribcr;
506  } else {
507  const Real z0 = z0_eff_arr(i, j, 0);
508  const Real ws_layer = ws10_eff_arr(i, j, 0);
509  const Real Rossby = ws_layer / (turbChoice.pbl_ysu_coriolis_freq * amrex::max(z0, Real(1.0e-4)));
510  Ribcr = amrex::min(Real(0.16) * std::pow(Real(1.0e-7) * Rossby, -Real(0.18)), Real(0.3));
511  }
512  }
513 
514  // Scan rib_base_arr(i,j,ksrf..khi) to find first crossing of Ribcr
515  // GAP 9: Seed from surface bulk Richardson number (WRF bl_ysu.F90 lines 620-621)
516  // WRF seeds the loop with br (surface bulk Ri) so ksrf is checked first.
517  int kpbl = ksrf;
518  Real Rib = rib_base_arr(i,j,ksrf); // Seed from surface level
519  bool above_critical = (Rib >= Ribcr); // Check if already above critical
520 
521  // Scan from ksrf+1 if not already above critical
522  for (int kk = ksrf+1; !above_critical && kk <= khi; ++kk) {
523  if (rib_base_arr(i,j,kk) >= Ribcr) {
524  kpbl = kk;
525  above_critical = true;
526  break;
527  }
528  kpbl = kk; // keep updating so kpbl = khi if never exceeded
529  }
530  pbli_arr(i, j, 0) = kpbl;
531  });
532  BL_PROFILE_VAR_STOP(prof_pblh);
533 
534  //
535  // FIRST ITERATION: Compute HGAMT/HGAMQ/VPERT using Pass 1 PBL height
536  // ========================================================================
537  // This produces initial estimates of countergradient corrections and VPERT,
538  // which will then be used to recompute rib_enhan_arr for a more accurate
539  // corrector pass.
540  //
541 
542  // (Reuse existing corrector loop code below, first pass...)
543 
544  //
545  // CORRECTOR: Apply YSU nonlocal countergradient correction (HGAMT/HGAMQ)
546  // following HND06 with countergradient modifications.
547  //
548  // This step raises the effective surface virtual potential temperature
549  // to account for convective contribution to buoyancy, leading to a
550  // higher diagnosed PBL height in unstable conditions.
551  //
552  // WRF reference:
553  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F
554  //
555  // WRF constants (module_bl_ysu.F lines 100-120):
556  // RLAM = 150. (mixing length scale in free atm., m)
557  // PRMIN = 0.5 (minimum Prandtl number)
558  // PRMAX = 4.0 (maximum Prandtl number)
559  // BRCR = 0.25 (bulk Richardson number critical value over land)
560  // CFAC = 7.8 (= const_b, surface layer factor in countergradient term)
561  // PHIFAC = 8.0 (phifac: convective wscale weight)
562  // SFCFRAC= 0.1 (surface layer depth fraction)
563  // GAMCRT = 3.0 (maximum heat countergradient, K)
564  // GAMCRQ = 2.e-3 (maximum moisture countergradient, kg/kg)
565  //
566  // WRF formulas (module_bl_ysu.F lines 220-260):
567  // GAMFAC = CFAC / (rho * wscale)
568  // HGAMT = min(GAMFAC * HFX / CPM, GAMCRT)
569  // HGAMQ = min(GAMFAC * QFX, GAMCRQ) [if over water, HGAMQ = 0]
570  // VPERT = HGAMT + EP1 * theta * HGAMQ
571  // THERMAL += max(VPERT, 0) [raises effective surface virtual pot. temp.]
572  //
573  // YSU K-profile below PBL (module_bl_ysu.F lines 190-220):
574  // K_mom = rho * wscale * kappa * z * (1 - z/h)^pfac
575  // K_theta = K_mom / Prt
576  // K_q = K_mom / Prq (Prq ~ Prt for moisture)
577  //
578  //
579  // PASS 2 — CORRECTOR: Apply surface-type-dependent Ribcr with VPERT enhancement
580  // following HND06 (Hong, Noh & Dudhia 2006), with countergradient corrections.
581  //
582  // This pass raises the effective surface virtual potential temperature
583  // to account for convective contribution to buoyancy, resulting in
584  // higher diagnosed PBL height in unstable conditions.
585  //
586  // WRF reference (module_bl_ysu.F lines 150-180):
587  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F
588  //
589  // WRF constants (module_bl_ysu.F lines 50-75):
590  // RLAM = 150. (mixing length scale in free atm., m)
591  // PRMIN = 0.5 (minimum Prandtl number)
592  // PRMAX = 4.0 (maximum Prandtl number)
593  // BRCR = 0.25 (over land, bulk Richardson number critical value)
594  // CFAC = 7.8 (= const_b, surface layer factor in countergradient term)
595  // PFAC = 2.0 (power law exponent for PBL height factor)
596  // SFCFRAC= 0.1 (surface layer depth fraction)
597  // GAMCRT = 3.0 (maximum heat countergradient, K)
598  // GAMCRQ = 2.e-3 (maximum moisture countergradient, kg/kg)
599  //
600  // WRF formulas (module_bl_ysu.F lines 200-250):
601  // GAMFAC = CFAC / (rho * wstar)
602  // HGAMT = min(GAMFAC * HFX / CPM, GAMCRT)
603  // HGAMQ = min(GAMFAC * QFX, GAMCRQ) [if over water, HGAMQ = 0]
604  // VPERT = HGAMT + EP1 * theta * HGAMQ
605  // THERMAL += max(VPERT, 0) [raises effective surface virtual pot. temp.]
606  //
607  // WRF diffusion coefficients below PBL (module_bl_ysu.F lines 300-350):
608  // K_mom = rho * wstar * kappa * z * (1 - z/h)^2
609  // K_theta = K_mom / Prt
610  // K_q = K_mom / Prq (Prq similar to Prt for moisture)
611  //
612  const Real const_b = turbChoice.pbl_mrf_const_b;
613  const Real sf = turbChoice.pbl_mrf_sf;
614  //constexpr Real prmin = Real(0.5);
615  //constexpr Real prmax = Real(4.0);
616  constexpr Real GAMCRT = Real(3.0); // WRF GAMCRT: max heat countergradient (K)
617  constexpr Real GAMCRQ = Real(2.e-3); // WRF GAMCRQ: max moisture countergradient (kg/kg)
618  const bool enable_ysu_countergradient = turbChoice.enable_ysu_countergradient;
619  const bool enable_ysu_unbounded_vpert = turbChoice.enable_mrf_unbounded_vpert;
620  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
621  {
622  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
623  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
624  // Determine surface-type-dependent critical Richardson number (same as Pass 1)
625  // Over land: Ribcr = 0.25; over water: Ribcr depends on Rossby number
626  Real Ribcr;
627  {
628  bool over_land = (!lmask_arr) || (lmask_arr(i, j, 0) == 1);
629  if (over_land) {
630  Ribcr = turbChoice.pbl_ysu_land_Ribcr;
631  } else {
632  const Real z0 = z0_eff_arr(i, j, 0);
633  const Real ws_layer = ws10_eff_arr(i, j, 0);
634  const Real Rossby = ws_layer / (turbChoice.pbl_ysu_coriolis_freq * amrex::max(z0, Real(1.0e-4)));
635  Ribcr = amrex::min(Real(0.16) * std::pow(Real(1.0e-7) * Rossby, -Real(0.18)), Real(0.3));
636  }
637  }
638 
639  /*
640  const Real t_layer = t10_eff_arr(i, j, 0);
641  const Real moisture_fraction = use_moisture ? q10_eff_arr(i, j, 0) : zero;
642  const Real t_layer_v = t_layer * (one + epsv * moisture_fraction);
643  */
644 
645  const amrex::Real z_sfc_col = (use_terrain_fitted_coords)
646  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
647  : zib;
648 
649  int kpbl = ksrf;
650  Real zval0 = zero, Rib0 = zero;
651  // GAP 9: Seed from surface bulk Richardson number (WRF bl_ysu.F90 lines 620-621)
652  // WRF seeds the loop with br (surface bulk Ri) so ksrf is checked first.
653  Real Rib = rib_enhan_arr(i,j,ksrf); // Seed from surface level
654  zval0 = (use_terrain_fitted_coords)
655  ? Compute_Zrel_AtCellCenter(i,j,ksrf,z_nd_arr)
656  : (ksrf + myhalf) * dz;
657  Rib0 = Rib;
658  bool above_critical = (Rib >= Ribcr); // Check if already above critical
659 
660  // Scan rib_enhan_arr to find first crossing of Ribcr with linear interpolation
661  for (int kk = ksrf+1; !above_critical && kk <= khi; ++kk) {
662  if (rib_enhan_arr(i,j,kk) >= Ribcr) { kpbl = kk; above_critical = true; break; }
663  zval0 = (use_terrain_fitted_coords)
664  ? Compute_Zrel_AtCellCenter(i,j,kk,z_nd_arr)
665  : (kk + myhalf) * dz;
666  Rib0 = rib_enhan_arr(i,j,kk);
667  kpbl = kk;
668  }
669 
670  // Compute bounds for pblh clamping
671  const Real z_sfc = (use_terrain_fitted_coords)
672  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
673  : zib;
674  const Real dz_terrain = (use_terrain_fitted_coords)
675  ? (Compute_Zrel_AtCellCenter(i, j, ksrf + 1, z_nd_arr) - z_sfc)
676  : dz;
677  const Real z_max = (use_terrain_fitted_coords)
678  ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr)
679  : (khi + myhalf) * dz;
680  const Real pblh_max = Real(0.9) * z_max;
681 
682  amrex::Real pblh_min;
683  if (turbChoice.enable_ysu_terrain_pblh_floor && use_terrain_fitted_coords) {
684  const amrex::Real met_h_klo = Compute_h_zeta_AtCellCenter(i, j, ksrf, dxInv, z_nd_arr);
685  const amrex::Real dz0 = met_h_klo / dz_inv;
686  pblh_min = amrex::max(z_sfc_col + Real(0.5)*dz0, Real(10.0));
687  } else {
688  pblh_min = amrex::max(z_sfc + Real(0.5) * dz_terrain, Real(10.0));
689  }
690 
691  // Linear interpolation from previous level if crossing detected
692  if (kpbl < khi && rib_enhan_arr(i,j,kpbl) >= Ribcr) {
693  const Real zval = (use_terrain_fitted_coords)
694  ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr)
695  : (kpbl + myhalf) * dz;
696  Rib = rib_enhan_arr(i,j,kpbl);
697  Real pblh_interp = zval0 + (zval - zval0) / (Rib - Rib0) * (Ribcr - Rib0);
698  pblh_corr_arr(i, j, 0) = amrex::max(amrex::min(pblh_interp, pblh_max), pblh_min);
699  } else {
700  pblh_corr_arr(i, j, 0) = pblh_min;
701  }
702  pbli_arr(i, j, 0) = kpbl;
703 
704  // Initialize countergradient fluxes for now
705  hgamt_arr(i, j, 0) = zero;
706  hgamq_arr(i, j, 0) = zero;
707  wstar_arr(i, j, 0) = zero;
708 
709  });
710  // Debug output (disabled):
711  /*
712  amrex::Print() << "PBL height computed for MRF scheme at level "
713  << pblh_corr_arr(2, 2, 0)
714  << std::endl;
715  amrex::Print() << "PBL Temp:" << t_surf_arr(2, 2, 0) << " "
716  << t10av_arr(2, 2, 0) << std::endl;
717  */
718 
719  //
720  // Recompute convective velocity scale (wstar) and countergradient fluxes
721  // (HGAMT/HGAMQ) using the corrector PBL height (pblh_corr_arr).
722  //
723  // Background: The corrector PBL height is now available. The Monin-Obukhov stability
724  // functions and countergradient fluxes depend strongly on the PBL height estimate,
725  // so these quantities must be recomputed after the corrector pass to ensure
726  // consistency between the PBL height diagnosis and the K-profile mixing lengths.
727  //
728  // Note: Previous versions computed both a predictor PBL height (without countergradient
729  // corrections) and a corrector PBL height (with countergradient effects). The inconsistency
730  // between these two estimates led to unrealistic mixing intensity since the stability
731  // parameter HOL = sf*h/L differs. The predictor pass has been removed; only the
732  // corrector PBL height (pblh_corr_arr) is now computed and used throughout.
733  //
734  // WRF implements WSCALE (convective velocity) and countergradient corrections
735  // only once before both the corrector loop and K-profile computations, ensuring
736  // consistency. ERF implements three passes: this loop performs the needed computation
737  // after the corrector pass is complete.
738  //
739  // YSU (HND06): WSCALE = u* / φ_m(h/L)
740  // Countergradient: HGAMT = min(CFAC * u* * θ*, GAMCRT), where CFAC=7.8, GAMCRT=3K
741  // WRF Reference: module_bl_ysu.F lines 220-250
742  const bool enable_ysu_sat_limiter = turbChoice.enable_ysu_sat_limiter;
743  const bool enable_ysu_rad_tend_limiter = turbChoice.enable_ysu_rad_tend_limiter;
744  const amrex::Real ysu_rad_tend_limiter_magnitude = turbChoice.ysu_rad_tend_limiter_magnitude;
745 
746  // ========================================================================
747  // Cloud-top detection for top-down mixing (H10 Section 3b)
748  // ========================================================================
749  const bool enable_ysu_topdown = turbChoice.enable_ysu_topdown;
750  const amrex::Real ysu_qcloud_threshold = turbChoice.ysu_qcloud_threshold;
751 
752  if (enable_ysu_topdown && use_moisture) {
753  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
754  {
755  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
756  cloud_top_arr(i, j, 0) = -1; // Initialize: no cloud found
757 
758  // Search from PBL top downward for cloud-top cell
759  int kpbl = pbli_arr(i, j, 0);
760  for (int kk = kpbl - 1; kk >= ksrf; --kk) {
761  amrex::Real qc_kk = zero, qi_kk = zero;
762  if (moisture_indices.qc >= 0)
763  qc_kk = cell_data(i, j, kk, moisture_indices.qc) / cell_data(i, j, kk, Rho_comp);
764  if (moisture_indices.qi >= 0)
765  qi_kk = cell_data(i, j, kk, moisture_indices.qi) / cell_data(i, j, kk, Rho_comp);
766 
767  if (qc_kk + qi_kk > ysu_qcloud_threshold) {
768  cloud_top_arr(i, j, 0) = kk;
769  break;
770  }
771  }
772  });
773  } else {
774  // Disable cloud-top detection
775  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
776  {
777  cloud_top_arr(i, j, 0) = -1;
778  });
779  }
780 
781  ParallelFor(xybx_work, [=, zero_d=zero, one_d=one, CONST_GRAV_d=CONST_GRAV] AMREX_GPU_DEVICE(int i, int j, int) noexcept
782  {
783  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
784  const Real t_layer = t10_eff_arr(i, j, 0);
785  Real obuk_val = ol_eff_arr(i, j, 0);
786  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
787  obuk_val = (obuk_val >= zero) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
788  }
789 
790  // Compute Monin-Obukhov stability parameter using corrector PBL height
791  // HOL = sf * h / L, where L is Monin-Obukhov length scale
792  // WRF Reference: module_bl_ysu.F lines 200-210
793  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F#L200-L210
794  const Real HOL = sf * pblh_corr_arr(i, j, 0) / obuk_val;
795  const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
796  const Real one_quarter = Real(1.0) / Real(4.0);
797  const Real phiM = (obuk_val > 0)
798  ? (1 + 5 * HOL_bounded)
799  : std::pow(
800  amrex::max(1 - 16 * HOL_bounded, Real(0.01)),
801  -one_quarter);
802  const Real phiM_safe = amrex::max(phiM, Real(0.01));
803 
804  // Convective velocity scale (wscale = u*/phi_m), now computed with pblh_corr_arr
805  // wscale is the characteristic turbulent velocity in the boundary layer.
806  // Bounds (u*/5 to 16*u*) prevent unrealistic values in very weak or strong convection.
807  // WRF Reference: module_bl_ysu.F L210-215
808  Real wscale = us_eff_arr(i, j, 0) / phiM_safe;
809  wscale = amrex::max(wscale, us_eff_arr(i, j, 0) / Real(5.0)); // Mechanical turbulence floor
810  wscale = amrex::min(wscale, Real(16.0) * us_eff_arr(i, j, 0)); // Free convection ceiling
811 
812  // Compute top-down convective velocity scale wstar3_down (H10 Eq. 12)
813  // This represents turbulence driven by radiative cooling at cloud top
814  amrex::Real wstar3_down = zero;
815  if (enable_ysu_topdown && cloud_top_arr(i, j, 0) >= ksrf) {
816  int k_cloud_top = cloud_top_arr(i, j, 0);
817 
818  // Longwave cooling flux at cloud top (W/m^2 equivalent in K*m/s)
819  // When RRTMGP is active, read from qheating_rates if available.
820  // Otherwise use zero (feature disabled when radiation is off).
821  // Compute LRAD by integrating the LW heating rate (component index 1)
822  // over the PBL column from surface (ksrf) to cloud top (k_cloud_top)
823  amrex::Real LRAD = zero;
824  if (has_qheating_rates) {
825  // Integrate LW heating rate from surface to cloud top
826  // qheating_rates contains heating rates (K/s); negative values represent cooling
827  for (int kk = ksrf; kk <= k_cloud_top; ++kk) {
828  amrex::Real ldz = (use_terrain_fitted_coords)
829  ? (z_nd_arr(i, j, kk+1) - z_nd_arr(i, j, kk))
830  : dz;
831  // Accumulate cooling (negative heating) integrated over height
832  // LW component is at index 1
833  LRAD += -qheat_arr(i, j, kk, 1) * ldz;
834  }
835  }
836 
837  // YSUNew radiative tendency limiter/smoothing
838  // Apply optional finite guards, bounds checking, and smoothing
839  amrex::Real LRAD_raw = LRAD; // Store raw value for diagnostics
840  amrex::Real LRAD_limited = LRAD;
841 
842  if (enable_ysu_rad_tend_limiter && has_qheating_rates) {
843  // Guard against NaN/Inf in the raw heating rate
844  if (!amrex::Math::isfinite(LRAD_raw)) {
845  LRAD_limited = zero; // Safe fallback: no radiative forcing
846  } else {
847  // Apply the magnitude limiter/bounds
848  // Clamp LRAD to [-limiter_magnitude, +limiter_magnitude]
849  const amrex::Real lim_mag = ysu_rad_tend_limiter_magnitude;
850  LRAD_limited = amrex::min(LRAD_raw, lim_mag);
851  LRAD_limited = amrex::max(LRAD_limited, -lim_mag);
852  }
853  // Use limited value for subsequent wstar computation
854  LRAD = LRAD_limited;
855  }
856 
857  // Top-down convective velocity (H10 Eq. 12):
858  // wstar_down^3 = g/theta * LRAD/(rho*cp) * pblh
859  const amrex::Real t_local = cell_data(i, j, k_cloud_top, RhoTheta_comp)
860  / cell_data(i, j, k_cloud_top, Rho_comp);
861  const amrex::Real rho_ct = cell_data(i, j, k_cloud_top, Rho_comp);
862  wstar3_down = amrex::max(CONST_GRAV_d / t_local * LRAD
863  / (rho_ct * amrex::Real(1004.0))
864  * pblh_corr_arr(i, j, 0), zero_d);
865  }
866  wstar3_down_arr(i, j, 0) = wstar3_down;
867 
868  // Recompute wstar3 and wscale with corrected pblh_corr_arr
869  const Real bfx0_corr = amrex::max(sflux_arr(i, j, 0), zero_d);
870  const Real t_dry = t10_eff_arr(i, j, 0);
871  const Real wstar3_corr = (CONST_GRAV / t_dry) * bfx0_corr * pblh_corr_arr(i, j, 0);
872  wstar3_arr(i, j, 0) = wstar3_corr;
873 
874  // Recompute wscale with corrected pblh:
875  const Real ust3 = us_eff_arr(i, j, 0) * us_eff_arr(i, j, 0) * us_eff_arr(i, j, 0);
876  Real wscale_corr = std::cbrt(ust3 + amrex::Real(8.0) * KAPPA * wstar3_corr * amrex::Real(0.5));
877  wscale_corr = amrex::min(wscale_corr, us_eff_arr(i, j, 0) * amrex::Real(16.0));
878  wscale_corr = amrex::max(wscale_corr, us_eff_arr(i, j, 0) / amrex::Real(5.0));
879  wstar_arr(i, j, 0) = wscale_corr;
880 
881  // Compute HGAMT with corrected wscale and WRF convention
882  // WRF bl_ysu.F90 line 687: gamfac = bfac/rhox/wscale
883  // WRF bl_ysu.F90 line 688: hgamt(i) = min(gamfac*hfx(i)/cp, gamcrt)
884  // Use sflux-based sfcflg for consistency with the precompute kernel (Fix 3)
885  bool SFCFLG = (sfcflg_arr(i, j, 0) > zero); // TRUE when unstable/neutral (sflux > 0)
886  const Real rho_sfc = cell_data(i, j, ksrf, Rho_comp);
887  const Real hfx_col = -rho_sfc * amrex::Real(1004.0) * us_eff_arr(i, j, 0) * ts_eff_arr(i, j, 0);
888  const Real gamfac = const_b / (rho_sfc * wscale_corr);
889  const Real HGAMT = (SFCFLG && enable_ysu_countergradient)
890  ? amrex::min(gamfac * hfx_col / amrex::Real(1004.0), GAMCRT)
891  : zero;
892 
893  // Compute HGAMQ with corrected wscale and WRF convention
894  // WRF bl_ysu.F90 line 689: hgamq(i) = min(gamfac*qfx(i), gamcrq)
895  const Real qfx_col = -rho_sfc * us_eff_arr(i, j, 0) * qs_eff_arr(i, j, 0);
896  const Real HGAMQ_raw = (SFCFLG && use_moisture && enable_ysu_countergradient)
897  ? gamfac * qfx_col
898  : zero;
899  Real HGAMQ = amrex::min(amrex::max(HGAMQ_raw, zero_d), GAMCRQ);
900 
901  // Land/water surface discrimination
902  if (lmask_arr && SFCFLG && use_moisture && enable_ysu_countergradient) {
903  bool is_land = (lmask_arr(i,j,0) == 1);
904  if (!is_land) HGAMQ = zero;
905  }
906 
907  // Saturation-Aware HGAMQ limiter
908  if (enable_ysu_sat_limiter && moisture_indices.qv >= 0 && SFCFLG && use_moisture && enable_ysu_countergradient) {
909  Real qv_klo = cell_data(i, j, ksrf, moisture_indices.qv) / cell_data(i, j, ksrf, Rho_comp);
910  Real T_klo = getTgivenRandRTh(cell_data(i, j, ksrf, Rho_comp),
911  cell_data(i, j, ksrf, RhoTheta_comp),
912  qv_klo);
913  Real p_klo = getPgivenRTh(cell_data(i, j, ksrf, RhoTheta_comp), qv_klo) * Real(0.01);
914  Real qsat_klo = zero;
915  erf_qsatw(T_klo, p_klo, qsat_klo);
916  Real rh_klo = (qsat_klo > Real(1.0e-10)) ? (qv_klo / qsat_klo) : zero;
917  if (rh_klo > Real(0.95)) {
918  Real rh_scaling = amrex::max(zero_d, (one - rh_klo) / Real(0.05));
919  HGAMQ *= rh_scaling;
920  }
921  }
922  // Store HGAMT/h and HGAMQ/h for implicit solver (normalized by corrected PBL height)
923  // Also compute VPERT = max(HGAMT + 0.61*θ*HGAMQ, 0) for Pass 3 use
924  // WRF Reference: module_bl_ysu.F lines 230-240 (THERMAL update with VPERT)
925  if (pbli_arr(i, j, 0) <= ksrf + 1) {
926  hgamt_arr(i, j, 0) = zero;
927  hgamq_arr(i, j, 0) = zero;
928  hgamu_arr(i, j, 0) = zero;
929  hgamv_arr(i, j, 0) = zero;
930  vpert_arr(i, j, 0) = zero;
931  } else {
932  const Real pblh = pblh_corr_arr(i, j, 0);
933  hgamt_arr(i, j, 0) = (enable_ysu_countergradient) ? HGAMT / pblh : zero;
934  hgamq_arr(i, j, 0) = (enable_ysu_countergradient && use_moisture) ? HGAMQ / pblh : zero;
935 
936  // Compute YSU momentum countergradient terms (hgamu, hgamv)
937  // WRF bl_ysu.F90 lines 694-696:
938  // brint = -15.9*ust^2/wspd * wstar3/wscale^4
939  // hgamu = brint * ux(1), hgamv = brint * vx(1)
940  //
941  // NOTE: WRF stores the raw (m/s) values here and divides by hpbl at
942  // the point of use (module_bl_ysu.F: -hgamu(i)/hpbl(i)). ERF
943  // instead normalizes at the source, exactly as HGAMT/HGAMQ are
944  // normalized above, so that all four HGAM* components share one
945  // contract: a countergradient *gradient* (units of phi per m)
946  // that the implicit solvers multiply by rho*alpha to form a flux.
947  // See EddyDiff::HGAMU_v in ERF_IndexDefines.H.
948  hgamu_arr(i, j, 0) = zero;
949  hgamv_arr(i, j, 0) = zero;
950  if (SFCFLG && enable_ysu_countergradient) {
951  const Real wspd_sfc = ws10_eff_arr(i, j, 0);
952  const Real ustar = us_eff_arr(i, j, 0);
953  wscale = wstar_arr(i, j, 0); // column-average wscale
954  const Real wstar3 = wstar3_arr(i, j, 0); // convective velocity scale cubed
955 
956  const Real wscale4 = amrex::max(wscale * wscale * wscale * wscale,
957  amrex::Real(1.0e-10));
958  const Real brint = amrex::Real(-15.9) * ustar * ustar
959  / amrex::max(wspd_sfc, amrex::Real(1.0e-4))
960  * wstar3 / wscale4;
961 
962  // ux(i,1) and vx(i,1) in WRF = cell-center velocity at lowest level (ksrf in ERF)
963  const Real u_klo = myhalf * (uvel(i, j, ksrf) + uvel(i+1, j, ksrf));
964  const Real v_klo = myhalf * (vvel(i, j, ksrf) + vvel(i, j+1, ksrf));
965  hgamu_arr(i, j, 0) = brint * u_klo / pblh; // 1/s
966  hgamv_arr(i, j, 0) = brint * v_klo / pblh; // 1/s
967  }
968 
969  // VPERT for Pass 3 diagnostic: unnormalized (not divided by pblh)
970  // This represents the virtual temperature perturbation at surface
971  // When enable_ysu_unbounded_vpert=true, use unbounded VPERT (ERF enhanced)
972  // When enable_ysu_unbounded_vpert=false, cap at GAMCRT (WRF-compatible)
973  // WRF Reference: module_bl_ysu.F lines 230-240
974  // GAP 8: Include WRF height limiter (WRF bl_ysu.F90 line 689-690)
975  if (enable_ysu_countergradient) {
976  // Compute height limiter: min(zl1/(sfcfrac*pblh), 1.0)
977  // This reduces VPERT contribution when the first cell is coarse relative to thin PBLs
978  const Real zl1_col = (use_terrain_fitted_coords)
979  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
980  : (ksrf + myhalf) * dz;
981  constexpr Real sfcfrac_h = amrex::Real(0.1); // WRF SFCFRAC
982  const Real height_lim = amrex::min(zl1_col / (sfcfrac_h * pblh), one_d);
983 
984  const Real VPERT_raw = HGAMT + amrex::Real(0.61) * t_layer * HGAMQ;
985  const Real VPERT_capped = enable_ysu_unbounded_vpert
986  ? VPERT_raw
987  : amrex::min(VPERT_raw, GAMCRT);
988  vpert_arr(i, j, 0) = amrex::max(VPERT_capped, zero_d) * height_lim;
989  } else {
990  vpert_arr(i, j, 0) = zero;
991  }
992  }
993 
994  });
995 
996  // ========================================================================
997  // RECOMPUTE rib_enhan_arr WITH UPDATED VPERT
998  // ========================================================================
999  // BUG FIX: The enhanced Rib array was computed with vpert_arr=0 during the
1000  // initial precompute phase (line 226). Now that VPERT has been computed and
1001  // stored in vpert_arr (lines 709-761), we must RECOMPUTE rib_enhan_arr to
1002  // incorporate the now-nonzero VPERT. This ensures Pass 3 (zero-Ri diagnostic)
1003  // uses the correctly enhanced Richardson numbers, matching WRF methodology
1004  // where thermal excess is computed between predictor and corrector passes.
1005  //
1006  // Reference: Hong et al. (2006) HND06, WRF module_bl_ysu.F lines 150-170
1007  //
1008  BL_PROFILE_VAR("YSUNew_Rib_Recompute", prof_rib_recomp);
1009  ParallelFor(gbx_work, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1010  {
1011  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1012  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
1013  const amrex::Real t_layer = t10_eff_arr(i, j, 0);
1014  const amrex::Real q_frac = use_moisture ? q10_eff_arr(i, j, 0) : zero;
1015  const amrex::Real t_layer_v = t_layer * (one + epsv * q_frac);
1016  // Now vpert_arr contains the computed thermal perturbation from the corrector pass
1017  const amrex::Real t_enh = t_layer_v + vpert_arr(i,j,0);
1018  const amrex::Real z_sfc = (use_terrain_fitted_coords)
1019  ? Compute_Zrel_AtCellCenter(i,j,ksrf,z_nd_arr) : zib;
1020  const amrex::Real zval = (use_terrain_fitted_coords)
1021  ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
1022  : (k + myhalf) * dz;
1023  const amrex::Real zrel = amrex::max(zval - z_sfc, amrex::Real(1.0e-4));
1024  const amrex::Real theta_v = (use_moisture && enable_ysu_liquid_theta)
1025  ? GetThetavl(i,j,k,cell_data,moisture_indices)
1026  : GetThetav(i,j,k,cell_data,moisture_indices);
1027  const amrex::Real theta_v_klo = (use_moisture && enable_ysu_liquid_theta)
1028  ? GetThetavl(i,j,ksrf,cell_data,moisture_indices)
1029  : GetThetav(i,j,ksrf,cell_data,moisture_indices);
1030  const amrex::Real ws2_raw = fourth * ((uvel(i,j,k)+uvel(i+1,j,k))*(uvel(i,j,k)+uvel(i+1,j,k))
1031  + (vvel(i,j,k)+vvel(i,j+1,k))*(vvel(i,j,k)+vvel(i,j+1,k)));
1032  // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
1033  const amrex::Real ws2 = (turbChoice.enable_vh96_shear_correction)
1034  ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
1035  : amrex::max(ws2_raw, amrex::Real(1.0));
1036  // Recompute ONLY rib_enhan_arr with updated vpert_arr; leave rib_base_arr unchanged
1037  rib_enhan_arr(i,j,k) = CONST_GRAV * zrel * (theta_v - t_enh) / (ws2 * theta_v_klo);
1038  });
1039  BL_PROFILE_VAR_STOP(prof_rib_recomp);
1040 
1041  // ========================================================================
1042  // RE-RUN CORRECTOR PASS 2 WITH UPDATED rib_enhan_arr
1043  // ========================================================================
1044  // Now that rib_enhan_arr has been recomputed with the nonzero VPERT values,
1045  // we need to re-scan for the PBL height using the corrector Ribcr criterion.
1046  // This corrects the earlier Pass 2 which used zero-VPERT Rib values.
1047  //
1048  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
1049  {
1050  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1051  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
1052  // Determine surface-type-dependent critical Richardson number (same as before)
1053  Real Ribcr;
1054  {
1055  bool over_land = (!lmask_arr) || (lmask_arr(i, j, 0) == 1);
1056  if (over_land) {
1057  Ribcr = turbChoice.pbl_ysu_land_Ribcr;
1058  } else {
1059  const Real z0 = z0_eff_arr(i, j, 0);
1060  const Real ws_layer = ws10_eff_arr(i, j, 0);
1061  const Real Rossby = ws_layer / (turbChoice.pbl_ysu_coriolis_freq * amrex::max(z0, Real(1.0e-4)));
1062  Ribcr = amrex::min(Real(0.16) * std::pow(Real(1.0e-7) * Rossby, -Real(0.18)), Real(0.3));
1063  }
1064  }
1065 
1066  // Rescan rib_enhan_arr with updated values
1067  int kpbl = ksrf;
1068  Real zval0 = zero, Rib0 = zero;
1069  Real Rib = rib_enhan_arr(i,j,ksrf);
1070  zval0 = (use_terrain_fitted_coords)
1071  ? Compute_Zrel_AtCellCenter(i,j,ksrf,z_nd_arr)
1072  : (ksrf + myhalf) * dz;
1073  Rib0 = Rib;
1074  bool above_critical = (Rib >= Ribcr);
1075 
1076  for (int kk = ksrf+1; !above_critical && kk <= khi; ++kk) {
1077  if (rib_enhan_arr(i,j,kk) >= Ribcr) { kpbl = kk; above_critical = true; break; }
1078  zval0 = (use_terrain_fitted_coords)
1079  ? Compute_Zrel_AtCellCenter(i,j,kk,z_nd_arr)
1080  : (kk + myhalf) * dz;
1081  Rib0 = rib_enhan_arr(i,j,kk);
1082  kpbl = kk;
1083  }
1084 
1085  // Compute bounds
1086  const Real z_sfc = (use_terrain_fitted_coords)
1087  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
1088  : zib;
1089  const Real dz_terrain = (use_terrain_fitted_coords)
1090  ? (Compute_Zrel_AtCellCenter(i, j, ksrf + 1, z_nd_arr) - z_sfc)
1091  : dz;
1092  const Real z_max = (use_terrain_fitted_coords)
1093  ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr)
1094  : (khi + myhalf) * dz;
1095  const Real pblh_max = Real(0.9) * z_max;
1096 
1097  amrex::Real pblh_min;
1098  if (turbChoice.enable_ysu_terrain_pblh_floor && use_terrain_fitted_coords) {
1099  const amrex::Real met_h_klo = Compute_h_zeta_AtCellCenter(i, j, ksrf, dxInv, z_nd_arr);
1100  const amrex::Real dz0 = met_h_klo / dz_inv;
1101  const amrex::Real z_sfc_col = (use_terrain_fitted_coords)
1102  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
1103  : zib;
1104  pblh_min = amrex::max(z_sfc_col + Real(0.5)*dz0, Real(10.0));
1105  } else {
1106  pblh_min = amrex::max(z_sfc + Real(0.5) * dz_terrain, Real(10.0));
1107  }
1108  pblh_floor_arr(i, j, 0) = pblh_min;
1109 
1110  // Update PBL height with the enhanced Rib values
1111  if (kpbl < khi && rib_enhan_arr(i,j,kpbl) >= Ribcr) {
1112  const Real zval = (use_terrain_fitted_coords)
1113  ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr)
1114  : (kpbl + myhalf) * dz;
1115  Rib = rib_enhan_arr(i,j,kpbl);
1116  Real pblh_interp = zval0 + (zval - zval0) / (Rib - Rib0) * (Ribcr - Rib0);
1117  pblh_corr_arr(i, j, 0) = amrex::max(amrex::min(pblh_interp, pblh_max), pblh_min);
1118  } else {
1119  pblh_corr_arr(i, j, 0) = pblh_min;
1120  }
1121  pbli_arr(i, j, 0) = kpbl;
1122  });
1123 
1124  // Apply PBLH spatial smoothing if enabled (Seibert et al. 2000 methodology)
1125  // Seibert et al. (2000): Review and intercomparison of operational methods
1126  // for the determination of the mixing height. Atmospheric Environment, 34, 1001-1027.
1127  // Spatial smoothing removes unphysical grid-to-grid noise from discrete Rib-crossing detection
1128  if (turbChoice.enable_pblh_smoothing) {
1129  // Smooth the tile's own columns, reading the halo the passes above filled.
1130  // The result is the same however the domain is split.
1131  ApplyPBLHSmoothing(pbl_height_corrector, xybx_tile,
1132  turbChoice.pblh_smoothing_weight,
1133  turbChoice.pblh_smoothing_passes,
1134  geom.Domain(), geom.periodicity());
1135  // The smoothed height stays at or above the column's own floor (half
1136  // a cell above its surface, at least 10 m): beside lower neighbours a
1137  // column over a building (erf.pbl_ib_aware) would otherwise be handed
1138  // a height below its first cell, and the K profile a negative depth.
1139  ParallelFor(xybx_work, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept {
1140  pblh_corr_arr(i, j, 0) = amrex::max(pblh_corr_arr(i, j, 0), pblh_floor_arr(i, j, 0));
1141  });
1142  }
1143 
1144  // Copy corrected PBL height into pblh_mf for SurfaceLayer storage.
1145  // pbl_height_corrector only covers this tile (grown by one), so loop over
1146  // the tile, not the valid box: with tiling in x/y the valid box reaches
1147  // past it. Under TileNoZ the tile spans the full column, and pblh_mf has
1148  // no ghost cells, so the tiles together still fill every cell.
1149  {
1150  auto pblh_out = pblh_mf.array(mfi);
1151  const Box& tbx = mfi.tilebox();
1152  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
1153  pblh_out(i, j, k) = pblh_corr_arr(i, j, 0);
1154  });
1155  }
1156  // ========================================================================
1157  // Extension scan using liquid potential temperature (WRF bl_ysu.F90 lines 733-769)
1158  // ========================================================================
1159  // Scan upward from Pass 2 PBLH using theta_li to extend kpbl for
1160  // cloud-topped boundary layers. Triggered when enable_ysu_topdown and
1161  // use_moisture are both true.
1162  //
1163  if (enable_ysu_topdown && use_moisture) {
1164  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
1165  {
1166  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1167  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
1168  const int kpblold = pbli_arr(i,j,0); // Starting PBLH index from Pass 2
1169  bool definebrup = false;
1170 
1171  // Surface liquid-theta virtual potential temperature reference
1172  const Real thermalli = GetThetavl(i, j, ksrf, cell_data, moisture_indices);
1173 
1174  // Upward scan using liquid-theta stability criterion
1175  for (int kk = kpblold; kk < khi; ++kk) {
1176  // Velocity shear for Richardson number calculation
1177  const Real ws2_raw = fourth * ((uvel(i,j,kk)+uvel(i+1,j,kk))*(uvel(i,j,kk)+uvel(i+1,j,kk))
1178  + (vvel(i,j,kk)+vvel(i,j+1,kk))*(vvel(i,j,kk)+vvel(i,j+1,kk)));
1179  // Vogelezang & Holtslag (1996): Add shear correction term to denominator instead of ad-hoc floor
1180  const Real ws2 = (turbChoice.enable_vh96_shear_correction)
1181  ? (ws2_raw + turbChoice.vh96_shear_const_b * us_eff_arr(i,j,0) * us_eff_arr(i,j,0))
1182  : amrex::max(ws2_raw, amrex::Real(1.0));
1183 
1184  // Elevation at level kk
1185  const Real z_sfc = (use_terrain_fitted_coords)
1186  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
1187  : zib;
1188  const Real zval_kk = (use_terrain_fitted_coords)
1189  ? Compute_Zrel_AtCellCenter(i, j, kk, z_nd_arr)
1190  : (kk + myhalf) * dz;
1191  const Real zrel_kk = amrex::max(zval_kk - z_sfc, amrex::Real(1.0e-4));
1192 
1193  // Liquid-theta virtual potential temperature at current level and surface
1194  const Real thlix_kk = GetThetavl(i, j, kk, cell_data, moisture_indices);
1195  const Real thlix_klo = GetThetavl(i, j, ksrf, cell_data, moisture_indices);
1196 
1197  // Bulk Richardson number using liquid-theta potential temperature
1198  const Real bruptmp = CONST_GRAV * zrel_kk * (thlix_kk - thermalli) / (ws2 * thlix_klo);
1199 
1200  // Stability check: threshold is zero for extension scan
1201  const bool stable = (bruptmp >= zero);
1202 
1203  if (definebrup) {
1204  pbli_arr(i,j,0) = kk;
1205  definebrup = false;
1206  }
1207 
1208  if (!stable) {
1209  // Continue scanning while Richardson number indicates mixing
1210  definebrup = true;
1211  }
1212  }
1213 
1214  // Bound result to domain top
1215  pbli_arr(i,j,0) = amrex::min(pbli_arr(i,j,0), izmax);
1216  });
1217  }
1218 
1219  //
1220  // PASS 3 — ZERO-RI DIAGNOSTIC: Compute PBL index with Ribcr=0.0 criterion
1221  //
1222  // Background: WRF employs three distinct PBL height estimates (module_bl_ysu.F lines 150-200):
1223  // Pass 1: Ribcr = surface-type-dependent (predictor, base surface temperature)
1224  // Pass 2: Ribcr = surface-type-dependent (corrector, enhanced surface temperature with VPERT)
1225  // Pass 3: Ribcr = 0.0 diagnostic (uses enhanced surface temperature, neutral criterion)
1226  //
1227  // The third pass with Ribcr=0 produces the PBL height where Richardson number
1228  // first becomes neutral (Rib >= 0). This is the "zero-Richardson" diagnostic
1229  // commonly used in observation-based PBL depth estimates (module_bl_ysu.F lines 190-200).
1230  //
1231  // Use of three passes: The corrector PBL height (Pass 2) is used for
1232  // computing the mixing intensity K = ρ*wstar*κ*z*(1-z/h)², while
1233  // the diagnostic height (Pass 3, Ribcr=0) determines the vertical extent
1234  // of the nonlocal mixing region. This separation provides:
1235  // - Realistic mixed-layer extent in convective conditions
1236  // - Mixing formula consistent with corrector diagnostics
1237  // - Consistency with WRF methodology
1238  //
1239  // Reference: Hong et al. (2006), module_bl_ysu.F
1240  //
1241  constexpr Real Ribcr_zero = zero; // Zero critical Richardson number for diagnostic pass
1242  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
1243  {
1244  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1245  // GAP 9: Seed from surface bulk Richardson number (WRF bl_ysu.F90 lines 620-621)
1246  // WRF seeds the loop with br (surface bulk Ri) so ksrf is checked first.
1247  int kpbl_zero = ksrf;
1248  Real Rib = rib_enhan_arr(i,j,ksrf); // Seed from surface level
1249  bool above_critical = (Rib >= Ribcr_zero); // Check if already above critical (Ribcr=0)
1250 
1251  // Scan rib_enhan_arr(i,j,ksrf..khi) for Ribcr=0 crossing (no interpolation needed)
1252  for (int kk = ksrf+1; !above_critical && kk <= khi; ++kk) {
1253  if (rib_enhan_arr(i,j,kk) >= Ribcr_zero) { kpbl_zero = kk; above_critical = true; break; }
1254  kpbl_zero = kk;
1255  }
1256  pbli_zero_arr(i, j, 0) = kpbl_zero;
1257  });
1258  BL_PROFILE_VAR_STOP(prof_pblh);
1259 
1260  // Compute entrainment diffusivity at PBL top cell (WRF bl_ysu.F90 lines 831-925)
1261  // Uses WRF's buoyancy-flux-based entrainment velocity (we), which accounts for
1262  // the actual inversion strength (dthvx) at the PBL top.
1263  {
1264  const bool enable_ysu_entrainment = turbChoice.enable_ysu_entrainment;
1265  const bool enable_ysu_cloud_pblh = turbChoice.enable_ysu_cloud_pblh;
1266 
1267  ParallelFor(xybx, [=, zero_d=zero] AMREX_GPU_DEVICE(int i, int j, int) noexcept
1268  {
1269  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1270  entr_arr(i,j,0) = zero;
1271  if (!enable_ysu_entrainment) return;
1272 
1273  const int kpbl = pbli_arr(i,j,0);
1274  const Real rho_kpbl = cell_data(i, j, kpbl, Rho_comp);
1275  const Real pblh = pblh_corr_arr(i,j,0);
1276  const Real wscale = wstar_arr(i,j,0);
1277  const Real ustar = us_eff_arr(i, j, 0);
1278  const Real ustar3 = ustar * ustar * ustar;
1279 
1280  // WRF bl_ysu.F90 lines 831-839: buoyancy-flux-based entrainment velocity
1281  // wm3 = wstar3 + 5*ust3 (combined convective + mechanical velocity scale)
1282  const Real wstar3_col = wstar3_arr(i,j,0);
1283  const Real wm3 = wstar3_col + amrex::Real(5.0) * ustar3;
1284  const Real wm2 = std::pow(wm3, amrex::Real(2.0) / amrex::Real(3.0));
1285 
1286  // bfxpbl = -0.15 * thvx(ksrf)/g * wm3 / hpbl
1287  // thvx at ksrf: use liquid-theta equivalent if enabled (WRF thlix path)
1288  const Real thvx_klo = (use_moisture && enable_ysu_liquid_theta)
1289  ? GetThetavl(i, j, ksrf, cell_data, moisture_indices)
1290  : GetThetav(i, j, ksrf, cell_data, moisture_indices);
1291  const Real bfxpbl = amrex::Real(-0.15) * thvx_klo / CONST_GRAV * wm3 / pblh;
1292 
1293  // dthvx = max(thvx(kpbl+1) - thvx(kpbl), tmin) [inversion strength at PBL top]
1294  // Guard against kpbl+1 exceeding domain top
1295  const int kpbl_p1 = amrex::min(kpbl + 1, izmax);
1296  const Real thvx_kpbl = (use_moisture && enable_ysu_liquid_theta)
1297  ? GetThetavl(i, j, kpbl, cell_data, moisture_indices)
1298  : GetThetav(i, j, kpbl, cell_data, moisture_indices);
1299  const Real thvx_kpbl_p1 = (use_moisture && enable_ysu_liquid_theta)
1300  ? GetThetavl(i, j, kpbl_p1, cell_data, moisture_indices)
1301  : GetThetav(i, j, kpbl_p1, cell_data, moisture_indices);
1302  const Real dthvx = amrex::max(thvx_kpbl_p1 - thvx_kpbl, amrex::Real(1.0e-2));
1303 
1304  // we = max(bfxpbl / dthvx, -sqrt(wm2))
1305  // we is negative for true (downward) entrainment
1306  const Real we = amrex::max(bfxpbl / dthvx, -std::sqrt(wm2));
1307 
1308  // dz at kpbl level (terrain-adjusted)
1309  const Real met_h = (use_terrain_fitted_coords)
1310  ? Compute_h_zeta_AtCellCenter(i, j, kpbl, dxInv, z_nd_arr) : one;
1311  const Real dz_kpbl = met_h / dz_inv;
1312 
1313  // Apply entrainment only when we < 0 (downward entrainment occurring)
1314  const Real K_cap = Real(5.0) * rho_kpbl * wscale * KAPPA * pblh * Real(0.01);
1315 
1316  // Cloudy entrainment correction (WRF bl_ysu.F90 lines 840-875)
1317  // Applied when qc+qi exceeds threshold at cell kpbl-1 (layer below PBL top)
1318  // Uses cloud liquid water content and buoyancy jump for entrainment efficiency.
1319  // Radiation coupling is not required; wstar3_2 and hgamt2 remain zero.
1320  Real we_final = we;
1321  Real K_entr_final = (we < zero) ? rho_kpbl * (-we) * dz_kpbl : zero;
1322 
1323  if (enable_ysu_topdown && use_moisture) {
1324  const int k_below = amrex::max(kpbl - 1, ksrf);
1325  const Real rho_k = cell_data(i, j, k_below, Rho_comp);
1326  Real qc_below = (moisture_indices.qc >= 0)
1327  ? cell_data(i, j, k_below, moisture_indices.qc) / rho_k : zero;
1328  Real qi_below = (moisture_indices.qi >= 0)
1329  ? cell_data(i, j, k_below, moisture_indices.qi) / rho_k : zero;
1330 
1331  constexpr Real cloud_thresh_entr = amrex::Real(0.01e-3); // 0.01 g/kg
1332 
1333  if ((qc_below + qi_below) > cloud_thresh_entr && kpbl >= ksrf + 2) {
1334  // Cloud correction path: liquid water at cell k_below triggers adjustment
1335 
1336  // Entrainment efficiency computed from liquid water content and
1337  // buoyancy jump between layers. Limited to 0.4 by WRF methodology.
1338  // Uses theta-li jump from k_below to kpbl+1 for inversion strength.
1339 
1340  // Liquid-theta virtual potential temperature at relevant levels
1341  const int k_p2 = amrex::min(kpbl + 1, izmax);
1342  const Real thlix_kbelow = GetThetavl(i, j, k_below, cell_data, moisture_indices);
1343  const Real thlix_kp2 = GetThetavl(i, j, k_p2, cell_data, moisture_indices);
1344 
1345  // Buoyancy jump using liquid-theta across two layers
1346  const Real dthvx_li = amrex::max(thlix_kp2 - thlix_kbelow, amrex::Real(0.1));
1347 
1348  // Entrainment efficiency from cloud liquid water
1349  // ent_eff = min(0.2 * 8 * (xlv/cp) * qc_below / dthvx_li, 0.4)
1350  constexpr Real xlv_over_cp = amrex::Real(2.5e6) / amrex::Real(1004.0);
1351  const Real ent_eff = amrex::min(amrex::Real(0.2) * amrex::Real(8.0)
1352  * xlv_over_cp * qc_below / dthvx_li,
1353  amrex::Real(0.4));
1354 
1355  // Entrainment velocity using cloudy inversion strength
1356  // Replaces clear-sky dthvx with dthvx_li in we calculation
1357  const Real we_cloud = amrex::max(bfxpbl / dthvx_li, -std::sqrt(wm2));
1358 
1359  // Cloud-top entrainment velocity contribution
1360  // bfxpbl_cloud = -ent_eff * max(sflux, 0) (surface buoyancy flux only, no radsum)
1361  const Real bfx0_cloud = amrex::max(sflux_arr(i, j, 0), zero_d);
1362  const Real bfxpbl_cloud = -ent_eff * bfx0_cloud;
1363  const Real we_cloud_top = amrex::max(bfxpbl_cloud / dthvx_li,
1364  -std::sqrt(std::pow(wm3, amrex::Real(2.0)/amrex::Real(3.0))));
1365 
1366  // Summed entrainment velocity
1367  we_final = we_cloud + we_cloud_top;
1368 
1369  // Recompute entrainment coefficient with adjusted entrainment velocity
1370  K_entr_final = (we_final < zero) ? rho_kpbl * (-we_final) * dz_kpbl : zero;
1371 
1372  // Note: WRF zeroes xkzm at k_below when cloud is detected.
1373  // In ERF, entrainment is applied only at kpbl level, making this implicit.
1374  }
1375  }
1376 
1377  entr_arr(i,j,0) = amrex::min(K_entr_final, K_cap);
1378 
1379  // Cloud PBLH extension: check for cloud at PBL top cell and extend if present
1380  if (enable_ysu_cloud_pblh && use_moisture) {
1381  const int kpbl_current = pbli_arr(i, j, 0);
1382 
1383  // Check for cloud condensate at PBL top cell
1384  amrex::Real qc_kpbl = zero, qi_kpbl = zero;
1385  if (moisture_indices.qc >= 0)
1386  qc_kpbl = cell_data(i, j, kpbl_current, moisture_indices.qc) / cell_data(i, j, kpbl_current, Rho_comp);
1387  if (moisture_indices.qi >= 0)
1388  qi_kpbl = cell_data(i, j, kpbl_current, moisture_indices.qi) / cell_data(i, j, kpbl_current, Rho_comp);
1389 
1390  // Extend PBL top by one cell if cloud present exceeds threshold
1391  if ((qc_kpbl + qi_kpbl) > ysu_qcloud_threshold) {
1392  pbli_arr(i, j, 0) = amrex::min(kpbl_current + 1, izmax);
1393  }
1394  }
1395  });
1396  }
1397 
1398  // Compute K_down for top-down cloud-driven mixing (H10 Eq. 11)
1399  FArrayBox K_down_fab(gbx, 1, The_Async_Arena());
1400  K_down_fab.setVal<RunOn::Device>(zero);
1401  const auto& K_down_arr = K_down_fab.array();
1402 
1403  if (enable_ysu_topdown) {
1404  ParallelFor(gbx, [=, zero_d=zero] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1405  {
1406  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1407  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
1408  K_down_arr(i,j,k) = zero;
1409  if (k >= pbli_arr(i,j,0) || k < ksrf) return; // only inside PBL, never inside an immersed solid
1410 
1411  const amrex::Real rho = cell_data(i,j,k,Rho_comp);
1412  const amrex::Real pblh = pblh_corr_arr(i,j,0);
1413  const amrex::Real zval = (use_terrain_fitted_coords)
1414  ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
1415  : (k + myhalf) * dz;
1416  const amrex::Real z_sfc = (use_terrain_fitted_coords)
1417  ? Compute_Zrel_AtCellCenter(i,j,ksrf,z_nd_arr) : zib;
1418  const amrex::Real zrel = zval - z_sfc;
1419  const amrex::Real pblh_rel = amrex::max(pblh - z_sfc, amrex::Real(1.0e-4));
1420 
1421  // Top-down K-profile (H10 Eq. 11): K_down = rho * wstar_down * kappa * (h-z) * (z/h)^2
1422  // Use stored top-down convective velocity scale (computed from LRAD in the kernel above)
1423  // wstar3_down_arr(i,j,0) = 0 when LRAD=0 (radiation not coupled), so K_down will be zero
1424  const amrex::Real wstar3_down_col = wstar3_down_arr(i,j,0);
1425  const amrex::Real wstar_down_eff = (wstar3_down_col > zero)
1426  ? std::cbrt(wstar3_down_col)
1427  : zero;
1428  const amrex::Real zfac_up = amrex::max(zrel / pblh_rel, zero_d);
1429  K_down_arr(i,j,k) = rho * wstar_down_eff * KAPPA
1430  * amrex::max(pblh_rel - zrel, zero_d)
1431  * zfac_up * zfac_up;
1432  });
1433  }
1434 
1435  // -- Compute diffusion coefficients --
1436 
1437  const Array4<Real>& K_turb = eddyViscosity.array(mfi);
1438 
1439  // Dirichlet flags to switch derivative stencil
1440  bool c_ext_dir_on_zlo = ((bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir));
1441  bool c_ext_dir_on_zhi = ((bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir));
1442  bool u_ext_dir_on_zlo = ((bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir));
1443  bool u_ext_dir_on_zhi = ((bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir));
1444  bool v_ext_dir_on_zlo = ((bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir));
1445  bool v_ext_dir_on_zhi = ((bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir));
1446 
1447 
1448  BL_PROFILE_VAR("YSUNew_Kprofile", prof_kprof);
1449  ParallelFor(gbx, [=, wstar3_arr_cap=wstar3_arr, zol1_arr_cap=zol1_arr, sfcflg_arr_cap=sfcflg_arr,
1450  zero_d=zero, one_d=one, two_d=two] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1451  {
1452  const int ksrf = ksurf_arr(i, j, 0); // local surface: first fluid cell of the column
1453  const Real zib = zib_arr(i, j, 0); // its height above the domain bottom
1454  // Guard: skip lateral ghost cells that may have uninitialized density.
1455  const Real rho_guard = cell_data(i, j, k, Rho_comp);
1456  if (rho_guard <= Real(0)) {
1457  K_turb(i, j, k, EddyDiff::Mom_v) = Real(0);
1458  K_turb(i, j, k, EddyDiff::Theta_v) = Real(0);
1459  K_turb(i, j, k, EddyDiff::Q_v) = Real(0);
1460  K_turb(i, j, k, EddyDiff::HGAMT_v) = Real(0);
1461  K_turb(i, j, k, EddyDiff::HGAMQ_v) = Real(0);
1462  K_turb(i, j, k, EddyDiff::HGAMU_v) = Real(0);
1463  K_turb(i, j, k, EddyDiff::HGAMV_v) = Real(0);
1464  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = Real(0);
1465  return;
1466  }
1467 
1468  Real obuk_val = ol_eff_arr(i, j, 0);
1469  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
1470  obuk_val = (obuk_val >= zero) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
1471  }
1472 
1473  const Real zval = (use_terrain_fitted_coords)
1474  ? Compute_Zrel_AtCellCenter(i, j, k, z_nd_arr)
1475  : (k + myhalf) * dz;
1476  const Real rho = cell_data(i, j, k, Rho_comp);
1477  if (k < ksrf) { // inside an immersed solid: no diffusivity
1478  K_turb(i, j, k, EddyDiff::Mom_v) = zero_d;
1479  K_turb(i, j, k, EddyDiff::Theta_v) = zero_d;
1480  K_turb(i, j, k, EddyDiff::Q_v) = zero_d;
1481  K_turb(i, j, k, EddyDiff::HGAMT_v) = zero_d;
1482  K_turb(i, j, k, EddyDiff::HGAMQ_v) = zero_d;
1483  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = zero_d;
1484  return;
1485  }
1486  const Real met_h_zeta = (use_terrain_fitted_coords)
1487  ? Compute_h_zeta_AtCellCenter(i, j, k, dxInv, z_nd_arr) : one;
1488  const Real dz_terrain = met_h_zeta / dz_inv;
1489 
1490  // Cloud-aware diffusion: compute cloud fraction from qc and qi
1491  // Cloud water threshold for detection (threshold ~ 0.01 g/kg, matching WRF)
1492  //
1493  // CLOUD-AWARE ENHANCEMENTS (ERF Extension):
1494  // This is an optional feature NOT in HND06. It adjusts the
1495  // stability functions in cloudy layers to improve representation of:
1496  // 1. Reduced turbulence near cloud tops (stable layers with clouds)
1497  // 2. Enhanced mixing from latent heat release (unstable layers with clouds)
1498  // 3. Fog/stratus evolution in marine boundary layers
1499  //
1500  // Physical justification:
1501  // - Clouds modify vertical buoyancy structure through radiative effects
1502  // - Latent heat release enhances convective mixing
1503  // - Cloud-top entrainment is distinct from clear-air turbulence
1504  //
1505  // Reference concept: Similar to WRF's IMVDIF cloud-aware parameterization
1506  // (Bretherton & Park 2009, see WRF module_bl_mynn.F).
1507  //
1508  // Disable via enable_ysu_countergradient=false for strict HND06.
1509  constexpr Real qc_threshold = Real(1.0e-5); // 0.01 g/kg in mixing ratio (aligned with WRF)
1510  Real qc_mix = zero;
1511  Real qi_mix = zero;
1512  if (use_moisture) {
1513  // qc (mixing ratio) = qc (density) / rho
1514  if (moisture_indices.qc >= 0) {
1515  qc_mix = cell_data(i, j, k, moisture_indices.qc) / rho;
1516  }
1517  if (moisture_indices.qi >= 0) {
1518  qi_mix = cell_data(i, j, k, moisture_indices.qi) / rho;
1519  }
1520  }
1521  const Real total_qcloud = qc_mix + qi_mix;
1522  const bool has_cloud = (total_qcloud > qc_threshold);
1523 
1524  // Select PBL extent index based on configuration:
1525  // - Default (pbl_mrf_use_zero_ri_extent=false): use pbli_arr (Ri=0.5 corrector)
1526  // This matches WRF behavior for better comparability
1527  // - Alternative (pbl_mrf_use_zero_ri_extent=true): use pbli_zero_arr (Ri=0)
1528  // This provides a ~30-60% taller mixing region for convective ABL cases
1529  const int pbli_extent = turbChoice.pbl_mrf_use_zero_ri_extent ? pbli_zero_arr(i, j, 0) : pbli_arr(i, j, 0);
1530 
1531  if (k < pbli_extent) {
1532  // Within PBL: use nonlocal mixing with diagnostic stability functions
1533  // WRF reference (module_bl_ysu.F lines 190-220):
1534  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F#L190-L220
1535  //
1536  // PBL extent selection:
1537  // - Default (pbl_ysu_use_zero_ri_extent=false): uses pbli_arr (Ri=0.5 corrector)
1538  // Matches WRF behavior for K-profile mixing
1539  // - Alternative (pbl_ysu_use_zero_ri_extent=true): uses pbli_zero_arr (Ri=0)
1540  // Provides physically appropriate mixed-layer depth with extended mixing region
1541  // The mixing intensity is governed by pblh_corr_arr used in the formula K = ρ*wscale*κ*z*(1-z/h)^pfac.
1542  // This two-level approach ensures realistic PBL height behavior across
1543  // different stability regimes while maintaining stable mixing coefficients.
1544  //
1545  // Key physics: YSU nonlocal scheme represents updrafts/downdrafts by
1546  // countergradient fluxes, enabling faster PBL growth than local schemes.
1547  //
1548  // Implement SFCFLG stable-side gating (WRF L808, L872-884)
1549  // When stable (BR > 0, i.e., sflux < 0), skip nonlocal PBL mixing
1550  // and use free-atmosphere Richardson scheme throughout PBL.
1551  // Use sflux-based sfcflg for unified definition (Fix 3)
1552  bool SFCFLG = (sfcflg_arr_cap(i, j, 0) > zero); // TRUE when unstable/neutral (sflux > 0)
1553 
1554  // Stability function phiM for momentum - BUSINGER-DYER form (WRF)
1555  // Unstable (L < 0): phiM = (1 - 16*sf*h/L)^(-1/4) [APHI16=16, exponent -1/4]
1556  // Stable (L > 0): phiM = 1 + 5*sf*h/L
1557  // Bound HOL to [-100, 100] to prevent numerical issues in extreme stability
1558  const Real HOL = sf * pblh_corr_arr(i, j, 0) / obuk_val;
1559  const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
1560 
1561  const Real one_quarter = Real(1.0) / Real(4.0);
1562  const Real phiM = (obuk_val > 0)
1563  ? (1 + 5 * HOL_bounded)
1564  : std::pow(
1565  amrex::max(1 - 16 * HOL_bounded, Real(0.01)),
1566  -one_quarter);
1567 
1568  // Stability function phit for heat/temperature - BUSINGER-DYER form
1569  // Unstable (L < 0): phit = (1 - 16*sf*h/L)^(-1/2)
1570  // Stable (L > 0): phit = 1 + 5*sf*h/L (same as phiM for stability)
1571  // Reference: HND06, Table 1; WRF L857-861
1572  const Real phit = (obuk_val > 0)
1573  ? (1 + 5 * HOL_bounded)
1574  : std::pow(
1575  amrex::max(1 - 16 * HOL_bounded, Real(0.01)),
1576  -one / two);
1577 
1578  // Cloud-aware adjustment: In cloudy regions, adjust stability damping
1579  // because clouds modify buoyancy oscillations through latent heating
1580  // and radiation. This improves representation of:
1581  // - Stratus-topped boundary layers (reduced mixing with clouds)
1582  // - Cumulus-capped boundary layers (enhanced mixing with convection)
1583  // This is a physically motivated extension not in original HND06.
1584  Real phit_cloud = phit;
1585  Real phiM_cloud = phiM;
1586  if (has_cloud && obuk_val > zero) {
1587  // Stable layers with clouds: reduce stability enhancement
1588  // Cloud presence reduces oscillations, damping is less effective
1589  // Reduction factor: up to 15-20% where cloud water exceeds threshold
1590  Real reduction_factor = one - Real(0.15) * amrex::min(total_qcloud / qc_threshold, one_d);
1591  // Reduce the phiM enhancement in stable conditions
1592  phiM_cloud = one + Real(5.0) * reduction_factor * sf * pblh_corr_arr(i, j, 0) / obuk_val;
1593  phit_cloud = one + Real(5.0) * reduction_factor * sf * pblh_corr_arr(i, j, 0) / obuk_val;
1594  } else if (has_cloud && obuk_val <= zero) {
1595  // Unstable layers with clouds: slightly enhance instability
1596  // Clouds warm the boundary layer through latent heat release
1597  // Boost factor: up to 5% where cloud water exceeds threshold
1598  Real cloud_boost = Real(1.0) + Real(0.05) * amrex::min(total_qcloud / qc_threshold, one_d);
1599  // Ensure numerically stable exponentiation: base must be in [0.01, 1]
1600  // Businger-Dyer form: phi_m = (1 - 16*HOL)^(-1/4) for unstable conditions
1601  phiM_cloud = std::pow(
1602  amrex::max(one_d - Real(16.0) * HOL_bounded / cloud_boost, Real(0.01)),
1603  -one_quarter);
1604  phit_cloud = std::pow(
1605  amrex::max(one_d - Real(16.0) * HOL_bounded / cloud_boost, Real(0.01)),
1606  -one_d / two_d);
1607  }
1608 
1609  // Use cloud-aware stability functions
1610  const Real phiM_eff = phiM_cloud;
1611  const Real phit_eff = phit_cloud;
1612 
1613  // SECTION A: Full Prandtl Number Calculation (WRF bl_ysu.F90 lines 948-971)
1614  // Three-component formula matching WRF exactly:
1615  // (a) prfac = bfac*karman*sfcfrac = 6.8*0.4*0.1 = 0.272
1616  // (b) prfac2 = free-convection correction that reduces Pr for strong convection
1617  // (c) prnumfac = height-dependent exponential blending
1618  //
1619  // Reference: https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F#L948-L971
1620 
1621  // Step 1: Compute WRF constants (WRF lines 948-950)
1622  constexpr Real bfac = amrex::Real(6.8); // WRF BFAC (not 7.8)
1623  constexpr Real sfcfrac = amrex::Real(0.1); // SFCFRAC
1624  const Real conpr = bfac * KAPPA * sfcfrac; // = 0.272 (not 0.312)
1625 
1626  // Step 2: Compute height-dependent zq_kp1 for prnumfac calculation
1627  const Real zq_kp1_prandtl = zval + myhalf * dz_terrain;
1628 
1629  // Step 3: Compute prfac (surface layer Prandtl correction)
1630  // WRF bl_ysu.F90 line 951: prfac = conpr if SFCFLG, else 0
1631  const Real prfac = SFCFLG ? conpr : zero;
1632 
1633  // Step 4: Compute prfac2 (free-convection Prandtl correction)
1634  // WRF bl_ysu.F90 lines 962-963:
1635  // prfac2 = 15.9*(wstar3+wstar3_2)/ust3 / (1 + 4*karman*(wstar3+wstar3_2)/ust3)
1636  const Real wstar3_col = wstar3_arr_cap(i, j, 0);
1637  const Real ust3 = us_eff_arr(i, j, 0) * us_eff_arr(i, j, 0) * us_eff_arr(i, j, 0);
1638  // Top-down wstar3 term (wstar3_2 in WRF bl_ysu.F90 line 949)
1639  // Non-zero only when top-down convection is active (requires radiation coupling for LRAD)
1640  const Real wstar3_2 = wstar3_down_arr(i, j, 0);
1641 
1642  Real prfac2 = zero;
1643  if (SFCFLG && ust3 > amrex::Real(1.0e-10)) {
1644  const Real wstar_tot3 = wstar3_col + wstar3_2;
1645  prfac2 = amrex::Real(15.9) * wstar_tot3 / ust3
1646  / (one + amrex::Real(4.0) * KAPPA * wstar_tot3 / ust3);
1647  }
1648 
1649  // Step 5: Compute prnumfac (height-dependent exponential blend)
1650  // WRF bl_ysu.F90 lines 964-965:
1651  // prnumfac = -3*(max(zq(k+1) - sfcfrac*hpbl, 0))^2 / hpbl^2
1652  const Real pblh = pblh_corr_arr(i, j, 0);
1653  const Real sfclayer = sfcfrac * pblh;
1654  const Real zdiff = amrex::max(zq_kp1_prandtl - sfclayer, zero_d);
1655  const Real prnumfac = amrex::Real(-3.0) * zdiff * zdiff / (pblh * pblh);
1656 
1657  // Step 6: Compute prnum0 base Prandtl (WRF lines 966-970)
1658  // prnum0 = phiH/phiM + prfac
1659  constexpr Real prmin_wrf = amrex::Real(0.25); // WRF minimum (different from ERF)
1660  constexpr Real prmax_wrf = amrex::Real(4.0);
1661  Real prnum0 = (phit_eff / phiM_eff) + prfac;
1662  prnum0 = amrex::min(amrex::max(prnum0, prmin_wrf), prmax_wrf);
1663 
1664  // Step 7: Compute heat Prandtl with free-convection correction and height blend
1665  // WRF uses prnum0_heat = prnum0 / (1 + prfac2*karman*sfcfrac) * exp(prnumfac)
1666  Real prnum0_heat = prnum0 / (one + prfac2 * KAPPA * sfcfrac);
1667  prnum0_heat = amrex::min(amrex::max(prnum0_heat, prmin_wrf), prmax_wrf);
1668  const Real Prt = one + (prnum0_heat - one) * std::exp(prnumfac);
1669 
1670  // Step 8: Compute moisture Prandtl (uses prnum0 WITHOUT the free-convection correction)
1671  // WRF uses prnum_q = prnum0 * exp(prnumfac)
1672  const Real prnum_q = one + (prnum0 - one) * std::exp(prnumfac);
1673 
1674  // Use pre-computed wstar from the dedicated recomputation loop (lines 465-545).
1675  // wstar_arr was computed with pblh_corr_arr to ensure consistency
1676  // between countergradient diagnostics and K-profile calculations.
1677  //const Real wstar = wstar_arr(i, j, 0);
1678 
1679  // SFCFLG gating: WRF skips nonlocal mixing in stable PBL (SFCFLG=.FALSE., BR>0, obuk_val>0)
1680  // In stable conditions, use free-atmosphere Richardson mixing instead.
1681  // WRF Reference: module_bl_ysu.F lines 200, 220-260
1682  if (SFCFLG) {
1683  // PHASE 13 CHANGES: Use layer-top interface height (zq_kp1) and
1684  // level-dependent wscalek for proper K-profile formulation.
1685  // WRF Reference: module_bl_ysu.F lines 943-961
1686 
1687  // Step 1: Compute zq_kp1 (top interface of cell k)
1688  // zq(k+1) in WRF = cell center + dz/2
1689  const Real zq_kp1 = zval + myhalf * dz_terrain;
1690 
1691  // Step 2: Get surface height and first-level height (zl1)
1692  /*
1693  const Real z_sfc = (use_terrain_fitted_coords)
1694  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
1695  : zib;
1696  */
1697  const Real zl1 = (use_terrain_fitted_coords)
1698  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
1699  : (ksrf + myhalf) * dz;
1700 
1701  // Step 3: Compute zfac using zq_kp1 and zl1 (WRF bl_ysu.F90 line 943)
1702  // zfac = min(max((1-(zq(k+1)-zl1)/(hpbl-zl1)), zfmin), 1.)
1703  constexpr Real zfacmin = amrex::Real(1.0e-8);
1704  const Real pblh_rel = amrex::max(pblh_corr_arr(i, j, 0) - zl1, amrex::Real(1.0e-4));
1705  const Real zfac = amrex::min(
1706  amrex::max(one_d - (zq_kp1 - zl1) / pblh_rel, zfacmin), one_d);
1707 
1708  // Step 4: Compute level-dependent wscalek (for K-profile, NOT for HGAMT/HGAMQ)
1709  // WRF stores wstar3 per column; we now have it from Phase 12
1710  const Real ust3_wscale = us_eff_arr(i, j, 0) * us_eff_arr(i, j, 0) * us_eff_arr(i, j, 0);
1711 
1712  Real wscalek_val;
1713  // Unstable/neutral: wscalek = (ust3 + phifac*karman*wstar3*(1-zfac))^(1/3)
1714  wscalek_val = std::cbrt(ust3_wscale + amrex::Real(8.0) * KAPPA * wstar3_col * (one - zfac));
1715 
1716  // Step 5: Compute K_m using zq_kp1 and wscalek_val (WRF line 961)
1717  // WRF: xkzm(i,k) = wscalek(k)*karman*zq(k+1)*zfac(k)**pfac
1718  // The mixing length is the height above the column's own
1719  // surface (zib with erf.pbl_ib_aware, zero otherwise).
1720  constexpr Real ckz_pbl = Real(0.001);
1721  const Real K_base = ckz_pbl * dz_terrain * rho;
1722  constexpr Real pfac = amrex::Real(2.0);
1723  K_turb(i, j, k, EddyDiff::Mom_v) = K_base + rho * wscalek_val * KAPPA * (zq_kp1 - zib) * std::pow(zfac, pfac);
1724 
1725  // Apply Prandtl number to get heat and moisture diffusivity
1726  K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prt;
1727  if (turbChoice.ysu_moistvars) {
1728  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Mom_v) / prnum_q;
1729  } else {
1730  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
1731  }
1732 
1733  // Add top-down contribution (H10 Eq. 11)
1734  if (enable_ysu_topdown && k < pbli_arr(i, j, 0)) {
1735  K_turb(i, j, k, EddyDiff::Mom_v) += K_down_arr(i, j, k);
1736  K_turb(i, j, k, EddyDiff::Theta_v) += K_down_arr(i, j, k);
1737  }
1738  } else {
1739  // SECTION B: Stable PBL (sfcflg=false, k < pbli_extent)
1740  // WRF uses wscalek = ust/phi_m(zq(k+1)/L), NOT Richardson mixing inside PBL
1741  // Richardson mixing is ONLY for k >= pbli_extent (free atmosphere)
1742  // WRF Reference: bl_ysu.F90 lines 951-957
1743 
1744  // Step 1: Compute zq_kp1 (top interface of cell k)
1745  const Real zq_kp1_stable = zval + myhalf * dz_terrain;
1746 
1747  // Step 2: Get first-level height (zl1)
1748  const Real zl1_stable = (use_terrain_fitted_coords)
1749  ? Compute_Zrel_AtCellCenter(i, j, ksrf, z_nd_arr)
1750  : (ksrf + myhalf) * dz;
1751 
1752  // Step 3: Compute zfac for stable PBL
1753  constexpr Real zfacmin_stable = amrex::Real(1.0e-8);
1754  const Real pblh_rel_stable = amrex::max(pblh_corr_arr(i, j, 0) - zl1_stable, amrex::Real(1.0e-4));
1755  const Real zfac_stable = amrex::min(
1756  amrex::max(one_d - (zq_kp1_stable - zl1_stable) / pblh_rel_stable, zfacmin_stable), one_d);
1757 
1758  // Step 4: Compute stable wscalek using phi_m
1759  // WRF bl_ysu.F90 lines 951-957: wscalek = ust / phi_m(zq(k+1)/L)
1760  const Real zol1_stable = zol1_arr_cap(i, j, 0); // stored from Phase 12
1761  const Real zol_ratio = (zq_kp1_stable - zib) / (zl1_stable - zib); // zq(k+1) / zl1, above the column's surface
1762  const Real phim_stable_arg = zol1_stable * zol_ratio; // (z/L) for level k+1
1763  // Enable QNSE stable functions if requested, otherwise use default linear form
1764  const Real enable_qnse_d = (turbChoice.enable_qnse_stable_functions) ? Real(1.0) : Real(0.0);
1765  const Real qnse_am_d = turbChoice.qnse_am;
1766  const Real qnse_bm_d = turbChoice.qnse_bm;
1767  const Real phim_stable = (enable_qnse_d > Real(0.5))
1768  ? ((one + qnse_am_d * phim_stable_arg) / (one + qnse_bm_d * phim_stable_arg))
1769  : (one + amrex::Real(5.0) * phim_stable_arg); // stable: phi_m = 1 + 5*(z/L)
1770  const Real wscalek_stable = amrex::max(
1771  us_eff_arr(i, j, 0) / amrex::max(phim_stable, amrex::Real(0.01)),
1772  amrex::Real(0.001));
1773 
1774  // Step 5: Compute K_m for stable PBL using wscalek
1775  constexpr Real ckz_pbl_stable = Real(0.001);
1776  const Real K_base_stable = ckz_pbl_stable * dz_terrain * rho;
1777  constexpr Real pfac_stable = amrex::Real(2.0);
1778  K_turb(i, j, k, EddyDiff::Mom_v) = K_base_stable + rho * wscalek_stable * KAPPA * (zq_kp1_stable - zib) * std::pow(zfac_stable, pfac_stable);
1779 
1780  // Step 6: Apply Prandtl number for stable PBL
1781  // For stable, prfac=0 (from step 3 of SECTION A)
1782  // prnum = phiH/phiM without height blending adjustment
1783  const Real prnum_stable = one + (prnum0 - one) * std::exp(prnumfac);
1784  K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / prnum_stable;
1785 
1786  if (turbChoice.ysu_moistvars) {
1787  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Mom_v) / prnum_q;
1788  } else {
1789  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
1790  }
1791  }
1792  } else if (k >= pbli_extent) {
1793  // Free atmosphere above PBL: use local Richardson number-dependent mixing
1794  // with H10 (Hong 2010) grid-adaptive lengthscale
1795  // WRF Reference: module_bl_ysu.F; H10 Section 3a
1796  // H10 = Hong, S.-Y., 2010: A new stable boundary-layer mixing scheme. QJRMS, 136, 1481-1496.
1797  const Real lambda_min = Real(30.0);
1798  const Real lambda_max = Real(300.0);
1799  const Real lambdadz = amrex::min(amrex::max(Real(0.1) * dz_terrain, lambda_min), lambda_max);
1800  const Real lscale = (lambdadz * KAPPA * zval) / (lambdadz + KAPPA * zval);
1801  Real dthetadz, dudz, dvdz;
1802  ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
1803  c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo,
1804  u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz,
1805  dudz, dvdz, moisture_indices, pbl_strat);
1806 
1807  // This branch is the free atmosphere above the PBL, so it only
1808  // reaches the first fluid cell of the column when the PBL index
1809  // collapses to it. The resolved gradients are unusable there
1810  // (ERF #4037), so fall back on the MOST profile. The effective
1811  // surface scales and the local surface index are used so that an
1812  // immersed surface is handled the same way as the domain bottom.
1813  if (k == ksrf) {
1814  const Real zrel_sl = amrex::max(zval - zib, Real(1.0e-4));
1815  const Real theta = cell_data(i, j, k, RhoTheta_comp) / rho;
1816  const Real qv = (moisture_indices.qv >= 0) ?
1817  cell_data(i, j, k, moisture_indices.qv) / rho : zero;
1818  PBLSurfaceLayerGradient sl;
1819  sl.u_star = us_eff_arr(i, j, 0);
1820  sl.tstar_v = ComputeVirtualTStarPBL(ts_eff_arr(i, j, 0),
1821  use_moisture ? qs_eff_arr(i, j, 0) : zero,
1822  theta, qv, use_moisture);
1823  sl.zval = zrel_sl;
1824  sl.zeta = zrel_sl / obuk_val;
1825  ApplySurfaceLayerGradientsPBL(sl, dthetadz, dudz, dvdz);
1826  }
1827 
1828  // Apply boundary safeguards to avoid numerical instability in calm conditions above PBL
1829  const Real dudz_safe = (k < izmax) ? dudz : zero;
1830  const Real dvdz_safe = (k < izmax) ? dvdz : zero;
1831 
1832  // Wind shear magnitude with safety threshold (WRF approach, optimized for GPU precision)
1833  // Using separate safety threshold avoids numerical issues with very small wind shear
1834  const Real wind_shear = dudz_safe * dudz_safe + dvdz_safe * dvdz_safe;
1835  const Real wind_shear_safe = std::max(wind_shear, Real(1.0e-8));
1836 
1837  // Use virtual potential temperature (θ_v) for Richardson number stability calculation.
1838  // This correctly accounts for moisture effects on buoyancy.
1839  // WRF Reference: module_bl_ysu.F uses THVX (virtual potential temperature).
1840  // For moist air, θ_v = θ * (1 + 0.61*q_v - q_l - q_i) ≈ θ * (1 + 0.61*q_v)
1841  const Real theta_v = GetThetav(i, j, k, cell_data, moisture_indices);
1842  const Real dtheta_v_dz = dthetadz;
1843 
1844  // Gradient Richardson number: Ri_g = (g/θ_v) * (dθ_v/dz) / (shear²)
1845  // Reference: WRF module_bl_ysu.F line ~450-456, Hong et al. 2006, Eqn. A18
1846  //
1847  // For STABILITY ANALYSIS:
1848  // Ri_g > 0.5: typically considered strongly stable (turbulence suppressed)
1849  // 0 < Ri_g < 0.5: weakly stable
1850  // Ri_g < 0: unstable (turbulence enhanced)
1851  // Ri_g < -100: very strong instability (apply lower bound for safety)
1852  //
1853  // Bound values from below (-100.0) and above (100.0)
1854  // to prevent extreme floating-point scales from causing numerical instability
1855  Real grad_Ri = CONST_GRAV / theta_v * dtheta_v_dz / wind_shear_safe;
1856  grad_Ri = std::max(std::min(grad_Ri, Real(100.0)), -Real(100.0));
1857 
1858  // GAP 10: In-cloud moist Richardson number modification
1859  // WRF bl_ysu.F90 lines 992-1000 (imvdif == 1, which is hardcoded in WRF)
1860  // Only apply when both the current AND next cell are in-cloud (qc+qi > 0.01e-3 kg/kg)
1861  if (use_moisture && moisture_indices.qc >= 0 && k < izmax) {
1862  const Real rho_k = cell_data(i,j,k, Rho_comp);
1863  const Real rho_kp1 = cell_data(i,j,k+1,Rho_comp);
1864  Real qc_k = (moisture_indices.qc >= 0) ? cell_data(i,j,k, moisture_indices.qc)/rho_k : zero;
1865  Real qc_kp1 = (moisture_indices.qc >= 0) ? cell_data(i,j,k+1,moisture_indices.qc)/rho_kp1 : zero;
1866  Real qi_k = (moisture_indices.qi >= 0) ? cell_data(i,j,k, moisture_indices.qi)/rho_k : zero;
1867  Real qi_kp1 = (moisture_indices.qi >= 0) ? cell_data(i,j,k+1,moisture_indices.qi)/rho_kp1 : zero;
1868 
1869  constexpr Real cloud_thresh = amrex::Real(0.01e-3); // 0.01 g/kg in kg/kg
1870  if ((qc_k + qi_k) > cloud_thresh && (qc_kp1 + qi_kp1) > cloud_thresh) {
1871  // WRF bl_ysu.F90 lines 997-1000:
1872  // qmean = 0.5*(qvx(k)+qvx(k+1))
1873  // tmean = 0.5*(tx(k)+tx(k+1)) [temperature, not theta]
1874  // alpha = xlv*qmean/rd/tmean
1875  // chi = xlv*xlv*qmean / (cp*rv*tmean*tmean)
1876  // ri = (1+alpha)*(ri - g^2/(ss*tmean*cp) * (chi-alpha)/(1+chi))
1877  const Real qv_k_mri = (moisture_indices.qv >= 0)
1878  ? cell_data(i,j,k, moisture_indices.qv) / rho_k : zero;
1879  const Real qv_kp1_mri = (moisture_indices.qv >= 0)
1880  ? cell_data(i,j,k+1,moisture_indices.qv) / rho_kp1 : zero;
1881  // Use proper temperature (not theta) matching WRF's tx(i,k) — Fix 2
1882  const Real T_k = getTgivenRandRTh(rho_k, cell_data(i,j,k, RhoTheta_comp), qv_k_mri);
1883  const Real T_kp1 = getTgivenRandRTh(rho_kp1, cell_data(i,j,k+1,RhoTheta_comp), qv_kp1_mri);
1884  const Real tmean = myhalf * (T_k + T_kp1);
1885  const Real qmean = myhalf * (qv_k_mri + qv_kp1_mri);
1886  constexpr Real xlv = amrex::Real(2.5e6); // latent heat of vaporization (J/kg)
1887  constexpr Real rd = amrex::Real(287.0); // gas constant dry air (J/kg/K)
1888  constexpr Real rv = amrex::Real(461.5); // gas constant water vapor (J/kg/K)
1889  constexpr Real cp = amrex::Real(1004.0); // specific heat dry air (J/kg/K)
1890  const Real alpha = xlv * qmean / (rd * tmean);
1891  const Real chi = xlv * xlv * qmean / (cp * rv * tmean * tmean);
1892  // Moist Ri correction:
1893  grad_Ri = (one + alpha) * (grad_Ri
1894  - CONST_GRAV * CONST_GRAV / wind_shear_safe / tmean / cp
1895  * (chi - alpha) / (one + chi));
1896  grad_Ri = amrex::max(amrex::min(grad_Ri, amrex::Real(100.0)), amrex::Real(-100.0));
1897  }
1898  }
1899 
1900  // YSU stability functions (Hong et al. 2006, MWR, Appendix A)
1901  // Reference: https://doi.org/10.1175/MWR3250.1
1902  // See equations A19-A20
1903  //
1904  // Prandtl number in stable regime:
1905  // Pr = 1 + 2.1 * Ri_g (Eqn. A19)
1906  // Stability functions for momentum and heat:
1907  // For Ri_g > 0 (stable):
1908  // f_m = 1 / ((1 + 5*Ri_g)²) (Eqn. A20a)
1909  // f_t = 1 / ((1 + 5*Ri_g)²) (Eqn. A20b, same as f_m for this scheme)
1910  // For Ri_g < 0 (unstable):
1911  // f_m = 1 - 8*Ri_g / (1 + 1.746*sqrt(-Ri_g)) (Eqn. A20d)
1912  // f_t = 1 - 8*Ri_g / (1 + 1.286*sqrt(-Ri_g)) (Eqn. A20c)
1913  // Protect against numerical errors causing negative grad_Ri
1914  const Real grad_Ri_safe = amrex::max(grad_Ri, -Real(100.0)); // Bound negative values
1915  const Real fm = (grad_Ri_safe > 0)
1916  ? one / ((one + Real(5.0) * grad_Ri_safe) * (one + Real(5.0) * grad_Ri_safe))
1917  : 1 - 8 * grad_Ri_safe / (1 + Real(1.746) * std::sqrt(amrex::max(-grad_Ri_safe, zero_d))); // Eqn. A20b/d
1918  const Real ft = (grad_Ri_safe > 0)
1919  ? one / ((one + Real(5.0) * grad_Ri_safe) * (one + Real(5.0) * grad_Ri_safe))
1920  : 1 - 8 * grad_Ri_safe / (1 + Real(1.286) * std::sqrt(amrex::max(-grad_Ri_safe, zero_d))); // Eqn. A20a/c
1921  const Real rl2wsp = rho * lscale * lscale * std::sqrt(wind_shear);
1922 
1923  K_turb(i, j, k, EddyDiff::Mom_v) = rl2wsp * fm;
1924  // H10 formulation: for stable conditions, K_t = K_m / Pr
1925  if (grad_Ri_safe > 0) {
1926  Real Pr = amrex::min(amrex::max(one + Real(2.1) * grad_Ri, Real(0.25)), Real(4.0));
1927  K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Pr;
1928  } else {
1929  K_turb(i, j, k, EddyDiff::Theta_v) = rl2wsp * ft;
1930  }
1931  // WRF YSU: moisture diffusivity matches heat in free atmosphere
1932  // Physics: Above PBL, Richardson number mixing applies same Ri_g
1933  // dependent mixing to both heat and moisture. The heat scaling (ft or K_m/Pr)
1934  // is used for moisture, which implicitly assumes equal turbulent
1935  // Prandtl and Schmidt numbers in the stable/unstable free atmosphere.
1936  // Q diffusivity matches heat in free atmosphere (same Ri_g-dependent mixing)
1937  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
1938  }
1939  // Note: The conditions (k < pbli_extent) and (k >= pbli_extent) are exhaustive,
1940  // so the else clause would be unreachable. All cells reach K-bounds clipping below.
1941 
1942  // Add entrainment diffusivity at PBL top cell (HND06 Eq. 6)
1943  // Applied only at k == pbli_arr(i,j,0) (the first cell above the PBL)
1944  if (k == pbli_arr(i,j,0)) {
1945  K_turb(i,j,k,EddyDiff::Mom_v) += entr_arr(i,j,0);
1946  K_turb(i,j,k,EddyDiff::Theta_v) += entr_arr(i,j,0);
1947  if (turbChoice.ysu_moistvars) {
1948  K_turb(i,j,k,EddyDiff::Q_v) += entr_arr(i,j,0);
1949  }
1950  }
1951 
1952  // Limit diffusion coefficients to physical bounds
1953  // These bounds ensure numerical stability and prevent unrealistic diffusivity
1954  // that could violate conservation principles or cause solver issues.
1955  //
1956  // HND06 Conservative Limits (module_bl_ysu.F lines 300-320):
1957  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_ysu.F#L300-L320
1958  // Kmin = 0.1 m²/s, Kmax = 300 m²/s
1959  //
1960  // These limits were calibrated for global forecast models and found to be
1961  // robust across a wide range of atmospheric conditions. They prevent:
1962  // - Excessive diffusion in calm conditions (Kmax bound)
1963  // - Numerical instability from too-small diffusion (Kmin bound)
1964  // - CFL violations from over-diffusive mixing
1965  //
1966  // Alternative (Higher Limits) H10 formulation (QJRMS Appendix A):
1967  // Kmin = ckz * dz * rho (where ckz = 0.001), Kmax = 1000 m²/s
1968  // These higher limits allow greater mixing in free atmosphere and are
1969  // more appropriate for high-resolution simulations.
1970  Real rhoKmin, rhoKmax;
1971  if (turbChoice.pbl_ysunew_highres_bounds) {
1972  // Hong et al. 2006, MWR, Appendix A: Higher limits for free atmosphere
1973  // These are recommended for high-resolution simulations (Δz < 100 m)
1974  // where the free atmosphere grid can resolve small-scale mixing.
1975  constexpr Real ckz = Real(0.001);
1976  constexpr Real Kmax = Real(1000.0);
1977  rhoKmin = ckz * dz_terrain * rho;
1978  rhoKmax = rho * Kmax;
1979  } else {
1980  // HND06 Conservative limits (default, used in global forecasts)
1981  // These provide good results for typical meteorological scales (Δz > 100 m)
1982  // and are the standard in legacy WRF configurations.
1983  constexpr Real Kmin = Real(0.1);
1984  constexpr Real Kmax = Real(300.0);
1985  rhoKmin = rho * Kmin;
1986  rhoKmax = rho * Kmax;
1987  }
1988 
1989  #ifdef ERF_USE_WINDFARM
1990  // Wind farm coupling: if TKE from Fitch/EWP model is available and positive at this cell,
1991  // augment wscale to reflect enhanced turbulence from turbine wakes.
1992  // This requires mf_vars_windfarm to be passed in (added to function signature below).
1993  // For now, add a placeholder comment:
1994  // TODO: Accept optional MultiFab* windfarm_tke argument. When non-null and
1995  // windfarm_tke(i,j,k) > 0, augment K_turb by rho * sqrt(windfarm_tke(i,j,k)) * lscale.
1996  #endif
1997 
1998  // SECTION C: Asymmetric Background Diffusivity Floor (WRF bl_ysu.F90 lines 141-142, 972-974)
1999  // Applied INSIDE the PBL only (k < pbli_extent)
2000  // WRF parameters:
2001  // xkzminm = 0.1 m²/s (background for momentum, added to xkzm)
2002  // xkzminh = 0.01 m²/s (background for heat and moisture, added to xkzh and xkzq)
2003  if (k < pbli_extent) {
2004  constexpr Real xkzminm = amrex::Real(0.1); // background for momentum
2005  constexpr Real xkzminh = amrex::Real(0.01); // background for heat/moisture
2006  K_turb(i, j, k, EddyDiff::Mom_v) += rho * xkzminm;
2007  K_turb(i, j, k, EddyDiff::Theta_v) += rho * xkzminh;
2008  K_turb(i, j, k, EddyDiff::Q_v) += rho * xkzminh;
2009  }
2010  // NOTE: The free-atmosphere branch (k >= pbli_extent) keeps the existing
2011  // ckz*dz*rho floor from the WRF high-res bounds parameterization.
2012 
2013  K_turb(i, j, k, EddyDiff::Mom_v) = std::max(
2014  std::min(K_turb(i, j, k, EddyDiff::Mom_v), rhoKmax), rhoKmin);
2015  K_turb(i, j, k, EddyDiff::Theta_v) = std::max(
2016  std::min(K_turb(i, j, k, EddyDiff::Theta_v), rhoKmax), rhoKmin);
2017  K_turb(i, j, k, EddyDiff::Q_v) = std::max(
2018  std::min(K_turb(i, j, k, EddyDiff::Q_v), rhoKmax), rhoKmin);
2019  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = pblh_corr_arr(i, j, 0);
2020 
2021  // Store countergradient correction terms (HGAMT/h and HGAMQ/h)
2022  // Use the selected PBL extent index (pbli_arr or pbli_zero_arr based on pbl_mrf_use_zero_ri_extent)
2023  // IMPORTANT — units convention for HGAMT_v / HGAMQ_v / HGAMU_v / HGAMV_v:
2024  // HGAMT_v and HGAMQ_v are stored normalized by pblh (i.e. HGAMT/h,
2025  // HGAMQ/h), matching the MRF scheme's convention. HGAMU_v/HGAMV_v
2026  // store the signed brint*u/v momentum countergradient term directly
2027  // (no pblh normalization, since brint already carries the 1/pblh^0
2028  // scaling from the WRF formulation — see hgamu_arr/hgamv_arr above).
2029  //
2030  // None of these four fields require an extra 1/dz (dz_inv) factor
2031  // when consumed by the implicit solver. In ERF_ImplicitDiff_T.cpp
2032  // (ImplicitDiffForStateLU_T and ImplicitDiffForMomLU_T), the RHS
2033  // countergradient contribution is:
2034  // RHS += Fact * rhoAlpha * gam / met_h_zeta
2035  // where `Fact = implicit_fac * dt * dz_inv` ALREADY contains the
2036  // single grid-spacing factor required. Multiplying by dz_inv a
2037  // second time double-applies the grid spacing and
2038  // produces an incorrect (grid-resolution-dependent) correction.
2039  // This was verified against Exec/CanonicalTests/ABL/MRF_YSUNew_Enhancements/
2040  // canonical/veryunstable_cgcheck_abl — the extra dz_inv breaks the
2041  // countergradient balance check in that case. Do NOT re-add dz_inv
2042  // to the gam_hi/gam_lo terms in ERF_ImplicitDiff_T.cpp.
2043  if (k < pbli_extent) {
2044  // Inside PBL: store the normalized countergradient terms
2045  K_turb(i, j, k, EddyDiff::HGAMT_v) = hgamt_arr(i, j, 0);
2046  K_turb(i, j, k, EddyDiff::HGAMQ_v) = hgamq_arr(i, j, 0);
2047  K_turb(i, j, k, EddyDiff::HGAMU_v) = hgamu_arr(i, j, 0);
2048  K_turb(i, j, k, EddyDiff::HGAMV_v) = hgamv_arr(i, j, 0);
2049  } else {
2050  // Outside PBL: zero countergradient
2051  K_turb(i, j, k, EddyDiff::HGAMT_v) = zero;
2052  K_turb(i, j, k, EddyDiff::HGAMQ_v) = zero;
2053  K_turb(i, j, k, EddyDiff::HGAMU_v) = zero;
2054  K_turb(i, j, k, EddyDiff::HGAMV_v) = zero;
2055  }
2056  });
2057  BL_PROFILE_VAR_STOP(prof_kprof);
2058  // Debug print disabled for production runs.
2059  // amrex::Print()<<" Turbulent Viscosity at cell "<<K_turb(2, 2, 2, EddyDiff::Mom_v)<<" "<<pblh_corr_arr(2, 2, 0)<<std::endl;
2060  // FOEXTRAP top and bottom ghost cells
2061  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
2062  {
2063  K_turb(i, j, klo-1, EddyDiff::Mom_v ) = K_turb(i, j, klo, EddyDiff::Mom_v );
2064  K_turb(i, j, klo-1, EddyDiff::Theta_v) = K_turb(i, j, klo, EddyDiff::Theta_v);
2065  K_turb(i, j, klo-1, EddyDiff::Q_v ) = K_turb(i, j, klo, EddyDiff::Q_v );
2066  K_turb(i, j, klo-1, EddyDiff::HGAMT_v) = K_turb(i, j, klo, EddyDiff::HGAMT_v);
2067  K_turb(i, j, klo-1, EddyDiff::HGAMQ_v) = K_turb(i, j, klo, EddyDiff::HGAMQ_v);
2068  K_turb(i, j, klo-1, EddyDiff::HGAMU_v) = K_turb(i, j, klo, EddyDiff::HGAMU_v);
2069  K_turb(i, j, klo-1, EddyDiff::HGAMV_v) = K_turb(i, j, klo, EddyDiff::HGAMV_v);
2070  K_turb(i, j, klo-1, EddyDiff::Turb_lengthscale) = K_turb(i, j, klo, EddyDiff::Turb_lengthscale);
2071  K_turb(i, j, khi+1, EddyDiff::Mom_v ) = K_turb(i, j, khi, EddyDiff::Mom_v );
2072  K_turb(i, j, khi+1, EddyDiff::Theta_v) = K_turb(i, j, khi, EddyDiff::Theta_v);
2073  K_turb(i, j, khi+1, EddyDiff::Q_v ) = K_turb(i, j, khi, EddyDiff::Q_v );
2074  K_turb(i, j, khi+1, EddyDiff::HGAMT_v) = K_turb(i, j, khi, EddyDiff::HGAMT_v);
2075  K_turb(i, j, khi+1, EddyDiff::HGAMQ_v) = K_turb(i, j, khi, EddyDiff::HGAMQ_v);
2076  K_turb(i, j, khi+1, EddyDiff::HGAMU_v) = K_turb(i, j, khi, EddyDiff::HGAMU_v);
2077  K_turb(i, j, khi+1, EddyDiff::HGAMV_v) = K_turb(i, j, khi, EddyDiff::HGAMV_v);
2078  K_turb(i, j, khi+1, EddyDiff::Turb_lengthscale) = K_turb(i, j, khi, EddyDiff::Turb_lengthscale);
2079  });
2080  }// mfi
2081  // Write YSUNew-computed PBLH back into SurfaceLayer so Beljaars correction
2082  // and diagnostics can use it, and update_pblh no longer aborts for YSUNew type.
2083  // REGRID NOTE: On regrid, SurfaceLayer is reallocated and PBLH returns to sentinel.
2084  // The driver must call ComputeDiffusivityYSUNew (or an equivalent bootstrap pass)
2085  // before consuming PBLH in update_fluxes() or Beljaars correction after any regrid.
2086  SurfLayer->set_pblh(level, pblh_mf);
2087 }
epsv
constexpr amrex::Real epsv
Definition: ERF_Constants.H:52
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:67
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:68
getTgivenRandRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:46
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:42
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:43
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
klo
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
use_moisture
const bool use_moisture
Definition: ERF_InitCustomPert_ABL.H:71
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
erf_qsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:264
GetThetav
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(const int &i, const int &j, const int &k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:74
GetThetavl
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetavl(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:94
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
ComputeVerticalDerivativesPBL
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)
Definition: ERF_PBLModels.H:304
ApplyPBLHSmoothing
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)
Apply spatial smoothing to PBLH field using 5-point stencil.
Definition: ERF_PBLModels.H:717
ApplySurfaceLayerGradientsPBL
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void ApplySurfaceLayerGradientsPBL(const PBLSurfaceLayerGradient &sl, amrex::Real &dthetadz, amrex::Real &dudz, amrex::Real &dvdz)
Definition: ERF_PBLModels.H:423
ComputeVirtualTStarPBL
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)
Definition: ERF_PBLModels.H:386
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Compute_Zrel_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Zrel_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:758
Compute_h_zeta_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:197
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:109
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:92
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:108
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:300
EddyDiff::Theta_v
@ Theta_v
Definition: ERF_IndexDefines.H:253
EddyDiff::Turb_lengthscale
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:257
EddyDiff::HGAMU_v
@ HGAMU_v
Definition: ERF_IndexDefines.H:265
EddyDiff::Q_v
@ Q_v
Definition: ERF_IndexDefines.H:256
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:252
EddyDiff::HGAMQ_v
@ HGAMQ_v
Definition: ERF_IndexDefines.H:264
EddyDiff::HGAMT_v
@ HGAMT_v
Definition: ERF_IndexDefines.H:263
EddyDiff::HGAMV_v
@ HGAMV_v
Definition: ERF_IndexDefines.H:266
LsmVar_SLM::ustar
@ ustar
Definition: ERF_SLM.H:28
LsmVar_SLM::ks
@ ks
Definition: ERF_SLM.H:74
LsmVar_SLM::qstar
@ qstar
Definition: ERF_SLM.H:30
LsmVar_SLM::tstar
@ tstar
Definition: ERF_SLM.H:29
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:26
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:31
SLM_DST::alpha
@ alpha
Definition: ERF_SLM.H:151
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:218
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:219
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
module_model_constants::cp
real(c_double), parameter cp
Definition: ERF_module_model_constants.F90:22
module_model_constants::xlv
real(c_double), parameter xlv
Definition: ERF_module_model_constants.F90:57
MoistureComponentIndices::qi
int qi
cloud ice
Definition: ERF_DataStruct.H:251
MoistureComponentIndices::qv
int qv
water vapor
Definition: ERF_DataStruct.H:249
MoistureComponentIndices::qc
int qc
cloud liquid water
Definition: ERF_DataStruct.H:250
PBLDerivativeDzInv_T
Functor for inverse vertical spacings for terrain-following grids using cell-center heights.
Definition: ERF_PBLModels.H:663
PBLSurfaceLayerGradient
Definition: ERF_PBLModels.H:366
PBLSurfaceLayerGradient::tstar_v
amrex::Real tstar_v
virtual temperature scale
Definition: ERF_PBLModels.H:368
PBLSurfaceLayerGradient::zeta
amrex::Real zeta
stability coordinate z/L
Definition: ERF_PBLModels.H:370
PBLSurfaceLayerGradient::u_star
amrex::Real u_star
friction velocity
Definition: ERF_PBLModels.H:367
PBLSurfaceLayerGradient::zval
amrex::Real zval
height of the cell center above the surface
Definition: ERF_PBLModels.H:369
TurbChoice::enable_mrf_unbounded_vpert
bool enable_mrf_unbounded_vpert
Whether MRF leaves VPERT unlimited by GAMCRT.
Definition: ERF_TurbStruct.H:996
TurbChoice::pbl_mrf_use_zero_ri_extent
bool pbl_mrf_use_zero_ri_extent
Whether MRF uses the Ri=0 K-profile extent.
Definition: ERF_TurbStruct.H:997
TurbChoice::pbl_strat_type
StratType pbl_strat_type
Definition: ERF_TurbStruct.H:819
TurbChoice::enable_vh96_shear_correction
bool enable_vh96_shear_correction
Whether Vogelezang & Holtslag (1996) shear-correction term is enabled.
Definition: ERF_TurbStruct.H:974
TurbChoice::enable_ysu_entrainment
bool enable_ysu_entrainment
Whether YSU entrainment-layer parameterization is enabled.
Definition: ERF_TurbStruct.H:962
TurbChoice::enable_ysu_sat_limiter
bool enable_ysu_sat_limiter
Whether YSU applies a saturation limiter to moisture countergradient terms.
Definition: ERF_TurbStruct.H:958
TurbChoice::enable_ysu_cloud_pblh
bool enable_ysu_cloud_pblh
Whether YSU cloud-based PBL-height detection is enabled.
Definition: ERF_TurbStruct.H:964
TurbChoice::pbl_mrf_const_b
amrex::Real pbl_mrf_const_b
MRF constant used to compute PBL height.
Definition: ERF_TurbStruct.H:982
TurbChoice::pblh_smoothing_weight
amrex::Real pblh_smoothing_weight
Center-cell weight in PBLH smoothing stencil (must be in [0,1]).
Definition: ERF_TurbStruct.H:978
TurbChoice::ysu_qcloud_threshold
amrex::Real ysu_qcloud_threshold
Cloud liquid water threshold for YSUNew [kg/kg].
Definition: ERF_TurbStruct.H:966
TurbChoice::qnse_bm
amrex::Real qnse_bm
Definition: ERF_TurbStruct.H:1006
TurbChoice::enable_qnse_stable_functions
bool enable_qnse_stable_functions
Definition: ERF_TurbStruct.H:1003
TurbChoice::enable_ysu_countergradient
bool enable_ysu_countergradient
Whether YSU countergradient corrections are enabled.
Definition: ERF_TurbStruct.H:954
TurbChoice::enable_pblh_smoothing
bool enable_pblh_smoothing
Whether spatial smoothing of diagnosed PBLH is enabled.
Definition: ERF_TurbStruct.H:976
TurbChoice::pblh_smoothing_passes
int pblh_smoothing_passes
Number of PBLH smoothing iterations to apply.
Definition: ERF_TurbStruct.H:977
TurbChoice::pbl_ib_z0
amrex::Real pbl_ib_z0
Definition: ERF_TurbStruct.H:993
TurbChoice::enable_ysu_topdown
bool enable_ysu_topdown
Whether YSU top-down mixing is enabled.
Definition: ERF_TurbStruct.H:960
TurbChoice::pbl_ib_aware
bool pbl_ib_aware
Definition: ERF_TurbStruct.H:992
TurbChoice::ysu_rad_tend_limiter_magnitude
amrex::Real ysu_rad_tend_limiter_magnitude
Definition: ERF_TurbStruct.H:972
TurbChoice::pbl_ysu_land_Ribcr
amrex::Real pbl_ysu_land_Ribcr
Critical bulk Richardson number over land for stable YSU conditions.
Definition: ERF_TurbStruct.H:948
TurbChoice::vh96_shear_const_b
amrex::Real vh96_shear_const_b
Vogelezang & Holtslag (1996) shear-correction constant b.
Definition: ERF_TurbStruct.H:975
TurbChoice::qnse_am
amrex::Real qnse_am
Definition: ERF_TurbStruct.H:1005
TurbChoice::pbl_mrf_sf
amrex::Real pbl_mrf_sf
MRF surface flux value used to compute PBL height.
Definition: ERF_TurbStruct.H:983
TurbChoice::pbl_ysunew_highres_bounds
bool pbl_ysunew_highres_bounds
Whether YSUNew applies high-resolution grid-dependent diffusivity bounds.
Definition: ERF_TurbStruct.H:969
TurbChoice::pbl_ysu_coriolis_freq
amrex::Real pbl_ysu_coriolis_freq
Coriolis frequency used by YSU-family PBL schemes.
Definition: ERF_TurbStruct.H:939
TurbChoice::ysu_moistvars
bool ysu_moistvars
Whether YSU applies turbulence to moisture variables.
Definition: ERF_TurbStruct.H:968
TurbChoice::enable_ysu_terrain_pblh_floor
bool enable_ysu_terrain_pblh_floor
Whether YSU applies a terrain-following PBL-height floor.
Definition: ERF_TurbStruct.H:956
TurbChoice::enable_ysu_liquid_theta
bool enable_ysu_liquid_theta
Whether YSU uses liquid-water virtual potential temperature for stability.
Definition: ERF_TurbStruct.H:952
TurbChoice::enable_ysu_rad_tend_limiter
bool enable_ysu_rad_tend_limiter
Definition: ERF_TurbStruct.H:971

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:
Here is the caller graph for this function: