ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeDiffusivityYSUNew.cpp File Reference
#include "ERF_SurfaceLayer.H"
#include "ERF_DirectionSelector.H"
#include "ERF_Diffusion.H"
#include "ERF_Constants.H"
#include "ERF_TurbStruct.H"
#include "ERF_PBLModels.H"
#include "ERF_TileNoZ.H"
#include "ERF_MoistUtils.H"
Include 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)
 

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

Referenced by ComputeTurbulentViscosity().

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