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

Referenced by ComputeTurbulentViscosity().

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