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

Referenced by ComputeTurbulentViscosity().

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