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 
)

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

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

Referenced by ComputeTurbulentViscosity().

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