ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeDiffusivityMRF.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 "ERF_PBLScaleAwareBlending.H"
Include dependency graph for ERF_ComputeDiffusivityMRF.cpp:

Functions

void ComputeDiffusivityMRF (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)
 

Function Documentation

◆ ComputeDiffusivityMRF()

void ComputeDiffusivityMRF ( 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 
)

Compute vertical diffusivity using the Medium-Range Forecast (MRF) boundary layer scheme.

Parameters
[in]xvelx-velocity field.
[in]yvely-velocity field.
[in]cons_inInput conserved variables.
[out]eddyViscosityEddy viscosity and diffusivity coefficients.
[in]geomGrid geometry.
[in]turbChoiceTurbulence model options.
[in]SurfLayerSurface layer data.
[in]use_terrain_fitted_coordsUse terrain-fitted coordinates.
[in]use_moistureEnable moisture components.
[in]levelCurrent AMR level.
[in]bc_ptrBoundary condition records.
[in]z_phys_ndNodal physical heights.
[in]z_phys_ccCell-centered physical heights.
[in]moisture_indicesIndices for moisture variables.
47 {
48  /*
49  ============================================================================
50  Medium-Range Forecast (MRF) Boundary Layer Parameterization Scheme
51  ============================================================================
52 
53  Implementation of the MRF (Medium Range Forecast) boundary layer scheme
54  based on Hong and Pan (1996), with enhanced moisture handling and
55  cloud-aware stability corrections developed for ERF.
56 
57  PRIMARY REFERENCES:
58  -------------------
59  - Hong, S. Y., and H.-L. Pan, 1996: Nonlocal Boundary Layer Vertical
60  Diffusion in a Medium-Range Forecast Model. Monthly Weather Review, 124,
61  2322-2339. https://doi.org/10.1175/1520-0493(1996)124<2322:NBLVDI>2.0.CO;2
62 
63  - WRF Reference Implementation: module_bl_mrf.F
64  https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F
65 
66  - Hong, S. Y., Y. Noh, and J. Dudhia, 2006: A new vertical diffusion
67  package with an explicit treatment of entrainment processes. Monthly
68  Weather Review, 134, 2318-2341. https://doi.org/10.1175/MWR3250.1
69  (YSU scheme, used for free atmosphere mixing)
70 
71  CORE ALGORITHM (Hong & Pan 1996):
72  ----------------------------------
73  1. Predictor-Corrector Bulk Richardson Number (Rib) PBL Height Diagnosis
74  - Rib = (g*z*(θ_v(z) - θ_v(0))) / (U²*θ_v(0))
75  - h = min(z where Rib > Rib_critical)
76 
77  2. Nonlocal Countergradient Flux Corrections
78  - HGAMT: thermal excess from sensible heat flux
79  - HGAMQ: moisture excess from latent heat flux
80  - VPERT: virtual potential temperature perturbation
81 
82  3. Stability-Dependent Mixing Lengths in PBL
83  - K_m = ρ * w_* * κ * z * (1 - z/h)²
84  - K_t = K_m / Pr_t, K_q = K_m / Pr_q
85 
86  4. Free Atmosphere Mixing via Richardson Number
87  - Uses YSU stability functions (Hong et al. 2006)
88  - Ri-dependent diffusivity above PBL
89 
90  ENHANCEMENTS IN ERF IMPLEMENTATION:
91  -----------------------------------
92 
93  A. Virtual Potential Temperature (θ_v) Treatment
94  - Proper handling of moisture effects on buoyancy
95  - Important for accurate PBL height and convective strength
96  - Not explicitly detailed in Hong & Pan (1996) but standard practice
97 
98  B. Cloud-Aware Stability Function Adjustments (OPTIONAL)
99  - Physically motivated extension beyond Hong & Pan (1996)
100  - Detects clouds via qc + qi > threshold (0.1 g/kg)
101  - In stable layers with clouds:
102  * Reduces stability damping (φ reduced by 10-20%)
103  * Represents suppressed turbulence due to stable stratification
104  - In unstable layers with clouds:
105  * Slightly increases instability enhancement
106  * Represents latent heat release effects
107  - Disabled by default for backward compatibility
108  - Reference: Conceptually similar to WRF's IMVDIF cloud handling
109 
110  C. Explicit Moisture Sign Convention Handling
111  - q_star uses WRF convention: negative for upward flux (evaporation)
112  - HGAMQ calculation: -const_b * u_* * q_* / w_*
113  - Produces positive HGAMQ for unstable (evaporating) conditions
114  - Corrected with MAX(HGAMQ, 0) to prevent negative values
115  - Land/water discrimination: HGAMQ zeroed over water surfaces
116 
117  D. Improved VPERT (Virtual Potential Temperature Perturbation)
118  Formulation
119  - VPERT = max(HGAMT + 0.61*θ*HGAMQ, 0)
120  - Enhancement vs WRF: Does not limit VPERT to GAMCRT
121  - WRF limits VPERT after adding moisture term (VPERT = min(VPERT, GAMCRT))
122  - This is physically incorrect: moisture contribution can exceed GAMCRT
123  - ERF only limits HGAMT (sensible heat), allowing larger VPERT when
124  both sensible and latent heat create strong surface heating
125  - This produces more realistic PBL heights in moist convection
126 
127  E. Prandtl Number with Stability Correction
128  - Pr_t = φ_t/φ_m + const_b * κ * sf
129  - Bounded: 0.5 ≤ Pr_t ≤ 4.0
130  - Consistency with WRF: const_b = 7.8, sf = 0.1
131  - Reference: Hong et al. (2006), Equation A17
132 
133  PARAMETER DEFAULTS (Matching WRF):
134  -----------------------------------
135  - const_b = 7.8 (heat flux weight in countergradient)
136  - sf = 0.1 (surface layer / PBL height ratio in stability functions)
137  - GAMCRT = 3.0 K (maximum heat countergradient)
138  - GAMCRQ = 2.0e-3 kg/kg (maximum moisture countergradient)
139  - Rib_crit = 0.5 (critical bulk Richardson number for PBL height)
140  - Pr_min = 0.5, Pr_max = 4.0 (Prandtl number bounds)
141  - Kmin = 0.1 m²/s, Kmax = 300 m²/s (diffusivity bounds, Hong & Pan 1996)
142 
143  STABILITY FUNCTIONS:
144  --------------------
145  Unstable (L < 0, HOL < 0):
146  phi_m = (1 - 8 * sf * h/L)^(-1/3)
147  phi_t = (1 - 16 * sf * h/L)^(-1/2)
148 
149  Stable (L > 0, HOL > 0):
150  phi_m = phi_t = 1 + 5 * sf * h/L
151 
152  MIXING ABOVE PBL (Free Atmosphere):
153  -----------------------------------
154  Uses YSU scheme (Hong et al. 2006, Appendix A) for Richardson number
155  dependent mixing. This avoids MRF oscillations in stable conditions.
156 
157  Gradient Richardson number: Ri_g = (g/theta_v) * (dtheta_v/dz) / ((du/dz)^2 + (dv/dz)^2)
158 
159  For Ri_g > 0 (stable):
160  f_m = 1 / ((1 + 5*Ri_g)^2)
161  f_t = 1 / ((1 + 5*Ri_g)^2)
162  Pr_t = 1 + 2.1*Ri_g, bounded to [0.25, 4.0]
163 
164  For Ri_g < 0 (unstable):
165  f_m = 1 - 8*Ri_g / (1 + 1.746*sqrt(-Ri_g))
166  f_t = 1 - 8*Ri_g / (1 + 1.286*sqrt(-Ri_g))
167 
168  NUMERICAL CONSIDERATIONS:
169  -------------------------
170  - GPU-optimized with parallel_for loops
171  - Safe bounds on shear and vertical derivatives
172  - Gradient Richardson number limited: -100 <= Ri_g <= 100
173  - Avoids division by zero with minimum shear threshold (1e-10)
174  - theta_v_klo and theta_v guarded against zero/negative values
175 
176  TESTING & VALIDATION:
177  ---------------------
178  - Tested against WRF_ideal unstable/stable ABL cases
179  - Backward compatible when enhancements are disabled
180  - Cloud detection threshold: qc + qi > 1e-4 kg/kg (0.1 g/kg)
181  */
182 
183  // Domain extent in z-dir
184  int klo = geom.Domain().smallEnd(2);
185  int khi = geom.Domain().bigEnd(2);
186 
187 #ifdef _OPENMP
188 #pragma omp parallel if (Gpu::notInLaunchRegion())
189 #endif
190  for (MFIter mfi(eddyViscosity, TileNoZ()); mfi.isValid(); ++mfi) {
191 
192  // Box operated on must span fill domain in z-dir
193  const Box& gbx = mfi.growntilebox(IntVect(1,1,0));
194  AMREX_ALWAYS_ASSERT( gbx.smallEnd(2) == klo &&
195  gbx.bigEnd(2) == khi );
196 
197  // create flattened boxes to store PBL height and related quantities
198  const GeometryData gdata = geom.data();
199  const Box xybx = PerpendicularBox<ZDir>(gbx, IntVect{0, 0, 0});
200 
201  // Pass 1 (predictor): PBLH with base surface temperature, no VPERT
202  // Pass 2 (wstar/VPERT): compute wstar, HGAMT, HGAMQ, VPERT from predictor height
203  // Pass 3 (corrector): PBLH with VPERT-enhanced surface temperature (WRF-consistent)
204  // WRF reference (module_bl_mrf.F lines 813-964):
205  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L813-L964
206  FArrayBox pbl_height_predictor(xybx, 1, The_Async_Arena()); // Pass 1: base t_layer_v
207  FArrayBox pbl_height_corrector(xybx, 1, The_Async_Arena()); // Pass 3: VPERT-enhanced
208  IArrayBox pbl_index(xybx, 1, The_Async_Arena());
209  IArrayBox pbl_index_zero_ri(xybx, 1, The_Async_Arena()); // Index for zero-Ri diagnostic pass
210  FArrayBox hgamt_fab(xybx, 1, The_Async_Arena()); // Store HGAMT/h (normalized countergradient)
211  FArrayBox hgamq_fab(xybx, 1, The_Async_Arena()); // Store HGAMQ/h (normalized countergradient)
212  FArrayBox wstar_fab(xybx, 1, The_Async_Arena()); // Convective velocity scale
213  FArrayBox vpert_fab(xybx, 1, The_Async_Arena()); // Virtual temperature perturbation VPERT
214  const auto& pblh_pred_arr = pbl_height_predictor.array(); // predictor (base t_layer_v)
215  const auto& pblh_corr_arr = pbl_height_corrector.array(); // corrector (VPERT-enhanced)
216  const auto& pbli_arr = pbl_index.array();
217  const auto& pbli_zero_arr = pbl_index_zero_ri.array(); // Zero-Ri diagnostic PBL index
218  const auto& hgamt_arr = hgamt_fab.array();
219  const auto& hgamq_arr = hgamq_fab.array();
220  const auto& wstar_arr = wstar_fab.array();
221  const auto& vpert_arr = vpert_fab.array();
222 
223  // Get some data in arrays
224  const auto& cell_data = cons_in.const_array(mfi);
225  const auto& uvel = xvel.const_array(mfi);
226  const auto& vvel = yvel.const_array(mfi);
227 
228  const Real Ribcr = turbChoice.pbl_mrf_Ribcr;
229  const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
230  const auto& t_star_arr = SurfLayer->get_t_star(level)->const_array(mfi);
231  const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
232  const auto& t10av_arr = SurfLayer->get_mac_avg(level, 2)->const_array(mfi);
233  const auto& q10av_arr = SurfLayer->get_mac_avg(level, 3)->const_array(mfi);
234  const auto& lmask_arr = (SurfLayer->get_lmask(level)) ?
235  SurfLayer->get_lmask(level)->const_array(mfi) :
236  Array4<int>{};
237  // Only retrieve z_phys_nd array if terrain-fitted coordinates are in use
238  const Array4<Real const> z_nd_arr = use_terrain_fitted_coords ? z_phys_nd->array(mfi)
239  : Array4<Real const>{};
240  const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
241 
242  //
243  // PASS 1 (PREDICTOR): Compute PBL height using base surface virtual temperature.
244  // No thermal excess (VPERT = 0). Produces pblh_pred_arr and pbli_arr.
245  // WRF reference (module_bl_mrf.F lines 813-842):
246  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L813-L842
247  //
248  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
249  {
250  const Real t_layer = t10av_arr(i, j, 0);
251  const Real moisture_fraction = use_moisture ? q10av_arr(i, j, 0) : Real(0);
252  const Real t_layer_v = t_layer * (one + epsv * moisture_fraction);
253 
254  Real zval, Rib;
255  int kpbl = klo;
256 
257  {
258  zval = (use_terrain_fitted_coords)
259  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
260  : (kpbl + myhalf) * gdata.CellSize(2);
261  const Real theta_v = GetThetav(i, j, kpbl, cell_data, moisture_indices);
262  // FIX: guard theta_v_klo against zero to prevent NaN in Rib
263  const Real theta_v_klo = amrex::max(GetThetav(i, j, klo, cell_data, moisture_indices), Real(1.0));
264  const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
265  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
266  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
267  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
268  const Real ws2 = amrex::max(ws2_raw, Real(1.0));
269  Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
270  }
271 
272  Real zval0 = zval, Rib0 = Rib;
273  bool above_critical = false;
274  while (!above_critical && ((kpbl + 1) <= khi)) {
275  zval0 = zval;
276  Rib0 = Rib;
277  kpbl += 1;
278 
279  zval = (use_terrain_fitted_coords)
280  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
281  : (kpbl + myhalf) * gdata.CellSize(2);
282  const Real theta_v = GetThetav(i, j, kpbl, cell_data, moisture_indices);
283  // FIX: guard theta_v_klo against zero to prevent NaN in Rib
284  const Real theta_v_klo = amrex::max(GetThetav(i, j, klo, cell_data, moisture_indices), Real(1.0));
285  const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
286  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
287  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
288  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
289  const Real ws2 = amrex::max(ws2_raw, Real(1.0));
290  Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
291  above_critical = (Rib >= Ribcr);
292  }
293 
294  const Real pblh_emp = (use_terrain_fitted_coords)
295  ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr)
296  : myhalf * gdata.CellSize(2);
297  const Real z_max = (use_terrain_fitted_coords)
298  ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr)
299  : (khi + myhalf) * gdata.CellSize(2);
300  const Real pblh_max = Real(0.9) * z_max;
301  const Real pblh_min = amrex::max(pblh_emp, Real(10.0));
302 
303  if (above_critical) {
304  // FIX: guard against division by zero when Rib == Rib0
305  Real pblh_interp;
306  const Real rib_diff = Rib - Rib0;
307  if (std::abs(rib_diff) > Real(1.0e-10)) {
308  pblh_interp = zval0 + (zval - zval0) / rib_diff * (Ribcr - Rib0);
309  } else {
310  // If Rib not changing, use zval0 as PBL height
311  pblh_interp = zval0;
312  }
313  pblh_pred_arr(i, j, 0) = amrex::max(amrex::min(pblh_interp, pblh_max), pblh_min);
314  pbli_arr(i, j, 0) = kpbl;
315  } else {
316  pblh_pred_arr(i, j, 0) = pblh_min;
317  pbli_arr(i, j, 0) = klo + 1;
318  }
319  });
320 
321  // FIX: guard against null dereference of q_star when use_moisture is false (dry case)
322  // get_q_star can return null pointer when the field is not allocated
323  amrex::MultiFab* q_star_mf = SurfLayer->get_q_star(level);
324  Array4<Real const> q_star_arr = (q_star_mf != nullptr) ? q_star_mf->const_array(mfi)
325  : Array4<Real const>{};
326 
327  const Real const_b = turbChoice.pbl_mrf_const_b;
328  const Real sf = turbChoice.pbl_mrf_sf;
329  constexpr Real prmin = Real(0.5);
330  constexpr Real prmax = Real(4.0);
331  constexpr Real GAMCRT = Real(3.0);
332  constexpr Real GAMCRQ = Real(2.e-3);
333  const bool enable_mrf_countergradient = turbChoice.enable_mrf_countergradient;
334  const bool enable_mrf_unbounded_vpert = turbChoice.enable_mrf_unbounded_vpert;
335 
336  //
337  // PASS 2 (WSTAR / VPERT): Compute wstar, HGAMT, HGAMQ, VPERT using the
338  // predictor PBL height (pblh_pred_arr). VPERT will be fed into the corrector
339  // Rib search in Pass 3 to raise the effective surface temperature.
340  // WRF reference (module_bl_mrf.F lines 857-880):
341  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L857-L880
342  //
343  ParallelFor(xybx, [=,zero_d=zero] AMREX_GPU_DEVICE(int i, int j, int) noexcept
344  {
345  const Real t_layer = t10av_arr(i, j, 0);
346  Real obuk_val = l_obuk_arr(i, j, 0);
347  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
348  obuk_val = (obuk_val >= zero_d) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
349  }
350 
351  // HOL computed from predictor height (WRF uses predictor height for wstar/VPERT)
352  const Real HOL = sf * pblh_pred_arr(i, j, 0) / obuk_val;
353  const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
354  const Real one_quarter = Real(0.25);
355  const Real phiM = (obuk_val > 0)
356  ? (1 + 5 * HOL_bounded)
357  : std::pow(amrex::max(1 - 16 * HOL_bounded, Real(0.01)), -one_quarter);
358  const Real phiM_safe = amrex::max(phiM, Real(0.01));
359 
360  // wstar = u* / phi_m
361  // Absolute bounds [0.01, 5.0] m/s: prevent division-by-zero (floor) and
362  // free-convection blow-up (ceiling), independent of u* for u*->0 safety.
363  Real wstar = u_star_arr(i, j, 0) / phiM_safe;
364  wstar = amrex::max(wstar, Real(0.01));
365  wstar = amrex::min(wstar, Real(5.0));
366 
367  bool SFCFLG = (obuk_val <= zero_d);
368  const Real HGAMT = (SFCFLG && enable_mrf_countergradient)
369  ? amrex::min(-const_b * u_star_arr(i, j, 0) * t_star_arr(i, j, 0) / wstar, GAMCRT)
370  : zero_d;
371 
372  Real HGAMQ = Real(0);
373  if (SFCFLG && use_moisture && enable_mrf_countergradient) {
374  const Real q_star = q_star_arr(i, j, 0);
375  const Real HGAMQ_calc = -const_b * u_star_arr(i, j, 0) * q_star / wstar;
376  HGAMQ = amrex::max(amrex::min(HGAMQ_calc, GAMCRQ), Real(0));
377 
378  if (lmask_arr) {
379  bool is_land = (lmask_arr(i,j,0) == 1);
380  if (!is_land) HGAMQ = zero_d;
381  }
382 
383  if (moisture_indices.qv >= 0) {
384  Real qv_klo = cell_data(i, j, klo, moisture_indices.qv) / cell_data(i, j, klo, Rho_comp);
385  Real T_klo = getTgivenRandRTh(cell_data(i, j, klo, Rho_comp),
386  cell_data(i, j, klo, RhoTheta_comp), qv_klo);
387  Real p_klo = getPgivenRTh(cell_data(i, j, klo, RhoTheta_comp), qv_klo) * Real(0.01);
388  Real qsat_klo = zero_d;
389  erf_qsatw(T_klo, p_klo, qsat_klo);
390  Real rh_klo = (qsat_klo > Real(1.0e-10)) ? (qv_klo / qsat_klo) : Real(0);
391  if (rh_klo > Real(0.95)) {
392  Real rh_scaling = amrex::max(zero_d, (Real(1) - rh_klo) / Real(0.05));
393  HGAMQ *= rh_scaling;
394  }
395  }
396  }
397 
398  // Compute VPERT = HGAMT + 0.61*theta*HGAMQ (virtual temperature perturbation)
399  // This will be added to the surface temperature in the corrector Rib search.
400  // WRF Reference: module_bl_mrf.F lines 879-880
401  if (pbli_arr(i, j, 0) <= klo + 1 || !enable_mrf_countergradient) {
402  vpert_arr(i, j, 0) = zero_d;
403  } else {
404  const Real VPERT_raw = HGAMT + epsv * t_layer * HGAMQ;
405  const Real VPERT_capped = enable_mrf_unbounded_vpert
406  ? VPERT_raw
407  : amrex::min(VPERT_raw, GAMCRT);
408  vpert_arr(i, j, 0) = amrex::max(VPERT_capped, zero_d);
409  }
410  });
411 
412  //
413  // PASS 3 (CORRECTOR): Recompute PBL height with VPERT-enhanced surface temperature.
414  // theta_s = theta_va + VPERT (Hong & Pan 1996, Eq. 4)
415  // This is the WRF-consistent corrector pass. pblh_corr_arr is used for the
416  // K-profile shape function K = rho*wstar*kappa*z*(1-z/h)^2.
417  // WRF reference (module_bl_mrf.F lines 932-964):
418  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L932-L964
419  //
420  ParallelFor(xybx, [=,one_d=one]
421  AMREX_GPU_DEVICE(int i, int j, int) noexcept
422  {
423  const Real t_layer = t10av_arr(i, j, 0);
424  const Real moisture_fraction = use_moisture ? q10av_arr(i, j, 0) : Real(0);
425  // VPERT-enhanced surface virtual temperature (WRF THERMAL variable)
426  const Real t_layer_v = t_layer * (one + epsv * moisture_fraction)
427  + vpert_arr(i, j, 0);
428 
429  Real obuk_val = l_obuk_arr(i, j, 0);
430  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
431  obuk_val = (obuk_val >= zero) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
432  }
433 
434  int kpbl = klo;
435  Real zval, Rib; // FIX: removed uninitialized zval0, Rib0
436  {
437  zval = (use_terrain_fitted_coords)
438  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
439  : (kpbl + myhalf) * gdata.CellSize(2);
440  const Real theta_v = GetThetav(i, j, kpbl, cell_data, moisture_indices);
441  // FIX: guard theta_v_klo against zero to prevent NaN in Rib
442  const Real theta_v_klo = amrex::max(GetThetav(i, j, klo, cell_data, moisture_indices), one_d);
443  const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
444  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
445  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
446  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
447  const Real ws2 = amrex::max(ws2_raw, one_d);
448  Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
449  }
450  Real zval0 = zval, Rib0 = Rib; // FIX: initialize here, mirrors Pass 1
451 
452  bool above_critical = false;
453  while (!above_critical && ((kpbl + 1) <= khi)) {
454  zval0 = zval;
455  Rib0 = Rib;
456  kpbl += 1;
457 
458  zval = (use_terrain_fitted_coords)
459  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
460  : (kpbl + myhalf) * gdata.CellSize(2);
461  const Real theta_v = GetThetav(i, j, kpbl, cell_data, moisture_indices);
462  // FIX: guard theta_v_klo against zero to prevent NaN in Rib
463  const Real theta_v_klo = amrex::max(GetThetav(i, j, klo, cell_data, moisture_indices), one_d);
464  const Real ws2_raw = fourth * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
465  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
466  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
467  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
468  const Real ws2 = amrex::max(ws2_raw, one_d);
469  Rib = CONST_GRAV * zval * (theta_v - t_layer_v) / (ws2 * theta_v_klo);
470  above_critical = (Rib >= Ribcr);
471  }
472 
473  const Real pblh_emp = (use_terrain_fitted_coords)
474  ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr)
475  : myhalf * gdata.CellSize(2);
476  const Real z_max = (use_terrain_fitted_coords)
477  ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr)
478  : (khi + myhalf) * gdata.CellSize(2);
479  const Real pblh_max = Real(0.9) * z_max;
480  const Real pblh_min = amrex::max(pblh_emp, Real(10.0));
481 
482  if (above_critical) {
483  // FIX: guard against division by zero when Rib == Rib0
484  Real pblh_interp;
485  const Real rib_diff = Rib - Rib0;
486  if (std::abs(rib_diff) > Real(1.0e-10)) {
487  pblh_interp = zval0 + (zval - zval0) / rib_diff * (Ribcr - Rib0);
488  } else {
489  // If Rib not changing, use zval0 as PBL height
490  pblh_interp = zval0;
491  }
492  pblh_corr_arr(i, j, 0) = amrex::max(amrex::min(pblh_interp, pblh_max), pblh_min);
493  pbli_arr(i, j, 0) = kpbl;
494  } else {
495  pblh_corr_arr(i, j, 0) = pblh_min;
496  pbli_arr(i, j, 0) = klo + 1;
497  }
498  });
499 
500  //
501  // PASS 4 (WSTAR RECOMPUTE): Recompute wstar, HGAMT, HGAMQ using the corrected
502  // PBL height (pblh_corr_arr) to ensure internal consistency between the K-profile
503  // amplitude and the countergradient fluxes. HOL = sf*h/L uses pblh_corr_arr.
504  // WRF reference (module_bl_mrf.F lines 857-880):
505  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L857-L880
506  //
507  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
508  {
509  //const Real t_layer = t10av_arr(i, j, 0);
510  Real obuk_val = l_obuk_arr(i, j, 0);
511  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
512  obuk_val = (obuk_val >= Real(0)) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
513  }
514 
515  // HOL now uses corrected PBLH for full internal consistency
516  const Real HOL = sf * pblh_corr_arr(i, j, 0) / obuk_val;
517  const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
518  const Real one_quarter = Real(0.25);
519  const Real phiM = (obuk_val > 0)
520  ? (1 + 5 * HOL_bounded)
521  : std::pow(amrex::max(1 - 16 * HOL_bounded, Real(0.01)), -one_quarter);
522  const Real phiM_safe = amrex::max(phiM, Real(0.01));
523 
524  // Absolute bounds [0.01, 5.0] m/s
525  Real wstar = u_star_arr(i, j, 0) / phiM_safe;
526  wstar = amrex::max(wstar, Real(0.01));
527  wstar = amrex::min(wstar, Real(5.0));
528  wstar_arr(i, j, 0) = wstar;
529 
530  bool SFCFLG = (obuk_val <= Real(0));
531  const Real HGAMT = (SFCFLG && enable_mrf_countergradient)
532  ? amrex::min(-const_b * u_star_arr(i, j, 0) * t_star_arr(i, j, 0) / wstar, GAMCRT)
533  : Real(0);
534 
535  Real HGAMQ = Real(0);
536  if (SFCFLG && use_moisture && enable_mrf_countergradient) {
537  const Real q_star = q_star_arr(i, j, 0);
538  const Real HGAMQ_calc = -const_b * u_star_arr(i, j, 0) * q_star / wstar;
539  HGAMQ = amrex::max(amrex::min(HGAMQ_calc, GAMCRQ), Real(0));
540 
541  if (lmask_arr) {
542  bool is_land = (lmask_arr(i,j,0) == 1);
543  if (!is_land) HGAMQ = Real(0);
544  }
545 
546  if (moisture_indices.qv >= 0) {
547  Real qv_klo = cell_data(i, j, klo, moisture_indices.qv) / cell_data(i, j, klo, Rho_comp);
548  Real T_klo = getTgivenRandRTh(cell_data(i, j, klo, Rho_comp),
549  cell_data(i, j, klo, RhoTheta_comp), qv_klo);
550  Real p_klo = getPgivenRTh(cell_data(i, j, klo, RhoTheta_comp), qv_klo) * Real(0.01);
551  Real qsat_klo = Real(0);
552  erf_qsatw(T_klo, p_klo, qsat_klo);
553  Real rh_klo = (qsat_klo > Real(1.0e-10)) ? (qv_klo / qsat_klo) : Real(0);
554  if (rh_klo > Real(0.95)) {
555  Real rh_scaling = amrex::max(Real(0), (Real(1) - rh_klo) / Real(0.05));
556  HGAMQ *= rh_scaling;
557  }
558  }
559  }
560 
561  if (pbli_arr(i, j, 0) <= klo + 1) {
562  hgamt_arr(i, j, 0) = Real(0);
563  hgamq_arr(i, j, 0) = Real(0);
564  } else {
565  const Real pblh = pblh_corr_arr(i, j, 0);
566  // FIX: guard against division by zero when pblh is extremely small
567  if (pblh > Real(1.0e-10)) {
568  hgamt_arr(i, j, 0) = (enable_mrf_countergradient) ? HGAMT / pblh : Real(0);
569  hgamq_arr(i, j, 0) = (enable_mrf_countergradient && use_moisture) ? HGAMQ / pblh : Real(0);
570  } else {
571  hgamt_arr(i, j, 0) = Real(0);
572  hgamq_arr(i, j, 0) = Real(0);
573  }
574  }
575  });
576 
577  //
578  // PASS 5 (ZERO-RI): Diagnostic PBL height with Ribcr=0, VPERT-enhanced surface temp.
579  // Used optionally (pbl_mrf_use_zero_ri_extent) to extend the nonlocal mixing region.
580  // WRF reference (module_bl_mrf.F lines 932-964):
581  // https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F#L932-L964
582  //
583  constexpr Real Ribcr_zero = Real(0);
584  ParallelFor(xybx, [=,one_d=one]
585  AMREX_GPU_DEVICE(int i, int j, int) noexcept
586  {
587  const Real t_layer = t10av_arr(i, j, 0);
588  const Real moisture_fraction = use_moisture ? q10av_arr(i, j, 0) : zero;
589  const Real t_layer_v = t_layer * (one + epsv * moisture_fraction);
590  const Real t_layer_v_enhanced = t_layer_v + vpert_arr(i, j, 0);
591 
592  int kpbl_zero = klo;
593  Real zval_zero, Rib_zero;
594  {
595  zval_zero = (use_terrain_fitted_coords)
596  ? Compute_Zrel_AtCellCenter(i, j, kpbl_zero, z_nd_arr)
597  : (kpbl_zero + myhalf) * gdata.CellSize(2);
598  const Real theta_v = GetThetav(i, j, kpbl_zero, cell_data, moisture_indices);
599  // FIX: guard theta_v_klo against zero to prevent NaN in Rib
600  const Real theta_v_klo = amrex::max(GetThetav(i, j, klo, cell_data, moisture_indices), one_d);
601  const Real ws2_raw = fourth * ( (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) *
602  (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) +
603  (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) *
604  (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) );
605  const Real ws2 = amrex::max(ws2_raw, one_d);
606  Rib_zero = CONST_GRAV * zval_zero * (theta_v - t_layer_v_enhanced) / (ws2 * theta_v_klo);
607  }
608 
609  bool above_critical_zero = false;
610  while (!above_critical_zero && ((kpbl_zero + 1) <= khi)) {
611  kpbl_zero += 1;
612  zval_zero = (use_terrain_fitted_coords)
613  ? Compute_Zrel_AtCellCenter(i, j, kpbl_zero, z_nd_arr)
614  : (kpbl_zero + myhalf) * gdata.CellSize(2);
615  const Real theta_v = GetThetav(i, j, kpbl_zero, cell_data, moisture_indices);
616  // FIX: guard theta_v_klo against zero to prevent NaN in Rib
617  const Real theta_v_klo = amrex::max(GetThetav(i, j, klo, cell_data, moisture_indices), one_d);
618  const Real ws2_raw = fourth * ( (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) *
619  (uvel(i, j, kpbl_zero) + uvel(i + 1, j, kpbl_zero)) +
620  (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) *
621  (vvel(i, j, kpbl_zero) + vvel(i, j + 1, kpbl_zero)) );
622  const Real ws2 = amrex::max(ws2_raw, one_d);
623  Rib_zero = CONST_GRAV * zval_zero * (theta_v - t_layer_v_enhanced) / (ws2 * theta_v_klo);
624  above_critical_zero = (Rib_zero >= Ribcr_zero);
625  }
626 
627  if (above_critical_zero) {
628  pbli_zero_arr(i, j, 0) = kpbl_zero;
629  } else {
630  pbli_zero_arr(i, j, 0) = klo + 1;
631  }
632  });
633 
634  // -- Compute diffusion coefficients --
635 
636  const Array4<Real>& K_turb = eddyViscosity.array(mfi);
637 
638  bool c_ext_dir_on_zlo = ((bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir));
639  bool c_ext_dir_on_zhi = ((bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir));
640  bool u_ext_dir_on_zlo = ((bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir));
641  bool u_ext_dir_on_zhi = ((bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir));
642  bool v_ext_dir_on_zlo = ((bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir));
643  bool v_ext_dir_on_zhi = ((bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir));
644 
645  const auto& dxInv = geom.InvCellSizeArray();
646  const Real dz_inv = geom.InvCellSize(2);
647  const int izmin = geom.Domain().smallEnd(2);
648  const int izmax = geom.Domain().bigEnd(2);
649 
650  // Blending parameters captured as scalars for GPU lambda.
651  const amrex::Real l_blend_length = turbChoice.pbl_blend_length;
652  const amrex::Real l_blend_cs = turbChoice.pbl_blend_cs;
653  const amrex::Real l_blend_cmax = turbChoice.pbl_blend_c_max;
654  //const bool l_use_smag_ceil= turbChoice.pbl_blend_use_smag;
655  // dx: use horizontal spacing at this level (assume dx = dy for regular grids).
656  const amrex::Real l_dx = geom.CellSize(0);
657 
658  ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
659  {
660  Real obuk_val = l_obuk_arr(i, j, 0);
661  if (std::abs(obuk_val) < amrex::Real(1.0e-10)) {
662  obuk_val = (obuk_val >= Real(0)) ? amrex::Real(1.0e-10) : amrex::Real(-1.0e-10);
663  }
664 
665  const Real zval = (use_terrain_fitted_coords)
666  ? Compute_Zrel_AtCellCenter(i, j, k, z_nd_arr)
667  : (k + myhalf) * gdata.CellSize(2);
668  const Real rho = cell_data(i, j, k, Rho_comp);
669  // FIX: skip ghost cells with uninitialized (zero) density — prevents
670  // division-by-zero inside GetThetav at lateral ghost cells of gbx
671  if (rho <= Real(0)) {
672  K_turb(i, j, k, EddyDiff::Mom_v) = Real(0);
673  K_turb(i, j, k, EddyDiff::Theta_v) = Real(0);
674  K_turb(i, j, k, EddyDiff::Q_v) = Real(0);
675  K_turb(i, j, k, EddyDiff::HGAMT_v) = Real(0);
676  K_turb(i, j, k, EddyDiff::HGAMQ_v) = Real(0);
677  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = Real(0);
678  return;
679  }
680  const Real met_h_zeta = (use_terrain_fitted_coords)
681  ? Compute_h_zeta_AtCellCenter(i, j, k, dxInv, z_nd_arr) : Real(1);
682  const Real dz_terrain = met_h_zeta / dz_inv;
683 
684  constexpr Real qc_threshold = Real(1.0e-4); // Cloud water/ice threshold (kg/kg)
685  Real qc_mix = Real(0);
686  Real qi_mix = Real(0);
687  if (use_moisture) {
688  if (moisture_indices.qc >= 0) {
689  qc_mix = cell_data(i, j, k, moisture_indices.qc) / rho;
690  }
691  if (moisture_indices.qi >= 0) {
692  qi_mix = cell_data(i, j, k, moisture_indices.qi) / rho;
693  }
694  }
695  const Real total_qcloud = qc_mix + qi_mix;
696  const bool has_cloud = turbChoice.enable_mrf_cloud_adjustment && (total_qcloud > qc_threshold);
697 
698  const int pbli_extent = turbChoice.pbl_mrf_use_zero_ri_extent ? pbli_zero_arr(i, j, 0) : pbli_arr(i, j, 0);
699 
700  if (k < pbli_extent) {
701  bool SFCFLG = (obuk_val <= Real(0));
702 
703  const Real HOL = sf * pblh_corr_arr(i, j, 0) / obuk_val;
704  const Real HOL_bounded = amrex::max(amrex::min(HOL, Real(100.0)), Real(-100.0));
705 
706  const Real one_quarter = Real(0.25);
707  const Real phiM = (obuk_val > 0)
708  ? (1 + 5 * HOL_bounded)
709  : std::pow(amrex::max(1 - 16 * HOL_bounded, Real(0.01)), -one_quarter);
710  const Real phit = (obuk_val > 0)
711  ? (1 + 5 * HOL_bounded)
712  : std::pow(amrex::max(1 - 16 * HOL_bounded, Real(0.01)), -Real(0.5));
713 
714  Real phit_cloud = phit;
715  Real phiM_cloud = phiM;
716  if (has_cloud && obuk_val > Real(0)) {
717  Real reduction_factor = Real(1) - Real(0.15) * amrex::min(total_qcloud / qc_threshold, Real(1));
718  phiM_cloud = Real(1) + Real(5.0) * reduction_factor * sf * pblh_corr_arr(i, j, 0) / obuk_val;
719  phit_cloud = Real(1) + Real(5.0) * reduction_factor * sf * pblh_corr_arr(i, j, 0) / obuk_val;
720  } else if (has_cloud && obuk_val <= Real(0)) {
721  Real cloud_boost = Real(1.0) + Real(0.05) * amrex::min(total_qcloud / qc_threshold, Real(1));
722  phiM_cloud = std::pow(amrex::max(Real(1) - Real(16.0) * HOL_bounded / cloud_boost, Real(0.01)), -one_quarter);
723  phit_cloud = std::pow(amrex::max(Real(1) - Real(16.0) * HOL_bounded / cloud_boost, Real(0.01)), -Real(0.5));
724  }
725 
726  const Real phiM_eff = phiM_cloud;
727  const Real phit_eff = phit_cloud;
728 
729  Real Prt_base = phit_eff / phiM_eff;
730  const Real Prt = amrex::min(amrex::max(Prt_base + const_b * KAPPA * sf, prmin), prmax);
731 
732  const Real wstar = wstar_arr(i, j, 0);
733 
734  if (SFCFLG) {
735  // K-profile: K = rho * wstar * kappa * zrel * (1 - zrel/pblh_rel)^2
736  // WRF Reference: module_bl_mrf.F L976-978
737  const Real z_sfc = (use_terrain_fitted_coords)
738  ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr)
739  : Real(0);
740  const Real zrel = zval - z_sfc;
741  const Real pblh = pblh_corr_arr(i, j, 0);
742  const Real pblh_rel = pblh - z_sfc;
743  const Real zfac = amrex::max(Real(1) - zrel / pblh_rel, Real(1.0e-8));
744 
745  K_turb(i, j, k, EddyDiff::Mom_v) = rho * wstar * KAPPA * zrel * zfac * zfac;
746  K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prt;
747 
748  if (turbChoice.mrf_moistvars) {
749  Real Prq_base = phit_eff / phiM_eff;
750  const Real Prq = amrex::min(amrex::max(Prq_base + const_b * KAPPA * sf, prmin), prmax);
751  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prq;
752  } else {
753  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
754  }
755  } else {
756  const Real lambda = Real(150.0);
757  const Real lscale = (KAPPA * zval * lambda) / (KAPPA * zval + lambda);
758  Real dthetadz, dudz, dvdz;
759  ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
760  c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo,
761  u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz,
762  dudz, dvdz, moisture_indices);
763 
764  const Real dudz_safe = (k < izmax) ? dudz : Real(0);
765  const Real dvdz_safe = (k < izmax) ? dvdz : Real(0);
766  const Real wind_shear = dudz_safe * dudz_safe + dvdz_safe * dvdz_safe;
767  const Real wind_shear_safe = std::max(wind_shear, Real(1.0e-8));
768 
769  // FIX: guard theta_v against zero/negative in grad_Ri denominator
770  const Real theta_v = amrex::max(GetThetav(i, j, k, cell_data, moisture_indices), Real(1.0));
771  const Real dtheta_v_dz = dthetadz;
772 
773  Real grad_Ri = CONST_GRAV / theta_v * dtheta_v_dz / wind_shear_safe;
774  grad_Ri = std::max(std::min(grad_Ri, Real(100.0)), -Real(100.0));
775 
776  const Real grad_Ri_safe = amrex::max(grad_Ri, -Real(100.0));
777  Real Pr_rich = Real(1) + Real(2.1) * grad_Ri;
778  const Real fm = (grad_Ri_safe > 0)
779  ? Real(1) / ((Real(1) + Real(5.0) * grad_Ri_safe) * (Real(1) + Real(5.0) * grad_Ri_safe))
780  : 1 - 8 * grad_Ri_safe / (1 + Real(1.746) * std::sqrt(amrex::max(-grad_Ri_safe, Real(0))));
781  const Real ft = (grad_Ri_safe > 0)
782  ? Real(1) / ((Real(1) + Real(5.0) * grad_Ri_safe) * (Real(1) + Real(5.0) * grad_Ri_safe))
783  : 1 - 8 * grad_Ri_safe / (1 + Real(1.286) * std::sqrt(amrex::max(-grad_Ri_safe, Real(0))));
784  const Real rl2wsp = rho * lscale * lscale * std::sqrt(wind_shear);
785 
786  Pr_rich = std::max(amrex::Real(0.25), std::min(Pr_rich, Real(4.0)));
787 
788  // In the stable regime fm and ft are the same function, so the
789  // Prandtl number has to be applied explicitly or momentum and heat
790  // mix identically. WRF (module_bl_mrf.F) takes heat as primary and
791  // scales momentum up by it:
792  // XKZH = DK/(1+5*RI)**2 ; XKZM = XKZH*PRNUM
793  // In the unstable regime fm and ft already differ (1.746 vs 1.286)
794  // and WRF applies no Prandtl factor there.
795  const Real Pr_mom = (grad_Ri_safe > 0) ? Pr_rich : Real(1);
796 
797  K_turb(i, j, k, EddyDiff::Mom_v) = rl2wsp * fm * Pr_mom;
798  K_turb(i, j, k, EddyDiff::Theta_v) = rl2wsp * ft;
799  if (use_moisture && turbChoice.mrf_moistvars) {
800  K_turb(i, j, k, EddyDiff::Q_v) = rl2wsp * ft;
801  } else {
802  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
803  }
804  }
805  } else if (k >= pbli_extent) {
806  const Real lambda = Real(150.0);
807  const Real lscale = (KAPPA * zval * lambda) / (KAPPA * zval + lambda);
808  Real dthetadz, dudz, dvdz;
809  ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
810  c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo,
811  u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz,
812  dudz, dvdz, moisture_indices);
813 
814  const Real dudz_safe = (k < izmax) ? dudz : Real(0);
815  const Real dvdz_safe = (k < izmax) ? dvdz : Real(0);
816  const Real wind_shear = dudz_safe * dudz_safe + dvdz_safe * dvdz_safe;
817  const Real wind_shear_safe = std::max(wind_shear, Real(1.0e-8));
818 
819  // FIX: guard theta_v against zero/negative in grad_Ri denominator
820  const Real theta_v = amrex::max(GetThetav(i, j, k, cell_data, moisture_indices), Real(1.0));
821  const Real dtheta_v_dz = dthetadz;
822 
823  Real grad_Ri = CONST_GRAV / theta_v * dtheta_v_dz / wind_shear_safe;
824  grad_Ri = std::max(std::min(grad_Ri, Real(100.0)), -Real(100.0));
825 
826  const Real grad_Ri_safe = amrex::max(grad_Ri, -Real(100.0));
827  Real Pr = Real(1) + Real(2.1) * grad_Ri;
828  const Real fm = (grad_Ri_safe > 0)
829  ? Real(1) / ((Real(1) + Real(5.0) * grad_Ri_safe) * (Real(1) + Real(5.0) * grad_Ri_safe))
830  : 1 - 8 * grad_Ri_safe / (1 + Real(1.746) * std::sqrt(amrex::max(-grad_Ri_safe, Real(0))));
831  const Real ft = (grad_Ri_safe > 0)
832  ? Real(1) / ((Real(1) + Real(5.0) * grad_Ri_safe) * (Real(1) + Real(5.0) * grad_Ri_safe))
833  : 1 - 8 * grad_Ri_safe / (1 + Real(1.286) * std::sqrt(amrex::max(-grad_Ri_safe, Real(0))));
834  const Real rl2wsp = rho * lscale * lscale * std::sqrt(wind_shear);
835 
836  Pr = std::max(amrex::Real(0.25), std::min(Pr, Real(4.0)));
837 
838  // In the stable regime fm and ft are the same function, so the Prandtl
839  // number has to be applied explicitly or momentum and heat mix
840  // identically. WRF (module_bl_mrf.F) takes heat as primary and scales
841  // momentum up by it:
842  // XKZH = DK/(1+5*RI)**2 ; XKZM = XKZH*PRNUM
843  // In the unstable regime fm and ft already differ (1.746 vs 1.286) and
844  // WRF applies no Prandtl factor there.
845  const Real Pr_mom = (grad_Ri_safe > 0) ? Pr : Real(1);
846 
847  K_turb(i, j, k, EddyDiff::Mom_v) = rl2wsp * fm * Pr_mom;
848  K_turb(i, j, k, EddyDiff::Theta_v) = rl2wsp * ft;
849  if (use_moisture && turbChoice.mrf_moistvars) {
850  K_turb(i, j, k, EddyDiff::Q_v) = rl2wsp * ft;
851  } else {
852  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
853  }
854  }
855 
856  // Scale-aware blending for grey-zone resolution (Boutle et al. 2014).
857  // ERF_PBLScaleAwareBlending.H. Gated by pbl_blend_length > 0.
858  // Applied to Theta_v and Q_v. Mom_v is not modified.
859  //
860  // K_turb holds rho*K [kg/m/s], but the ceiling applied inside
861  // pbl_kh_blend_and_cap is a bare diffusivity [m^2/s]. Hand the helper the
862  // kinematic value and restore the density weighting on the way out; passing
863  // rho*K straight through would make the effective cap scale like 1/rho and
864  // hence increase with height (issue #3580).
865  if (l_blend_length > 0.0) {
866  const Real rho_inv = Real(1.0) / rho;
867 
868  // For MRF, use power-law ceiling (SmnSmn not available)
869  K_turb(i, j, k, EddyDiff::Theta_v) = rho * pbl_kh_blend_and_cap(
870  K_turb(i, j, k, EddyDiff::Theta_v) * rho_inv,
871  l_dx, l_blend_length, l_blend_cs, l_blend_cmax,
872  amrex::Real(-1.0), false);
873 
874  K_turb(i, j, k, EddyDiff::Q_v) = rho * pbl_kh_blend_and_cap(
875  K_turb(i, j, k, EddyDiff::Q_v) * rho_inv,
876  l_dx, l_blend_length, l_blend_cs, l_blend_cmax,
877  amrex::Real(-1.0), false);
878  }
879 
880  // Limit diffusion coefficients to physical bounds
881  // Hong & Pan (1996): Kmin=0.1, Kmax=300 m^2/s (module_bl_mrf.F lines 1014-1025)
882  // Hong et al. (2006): Kmin=ckz*dz*rho, Kmax=1000 m^2/s (high-res option)
883  Real rhoKmin, rhoKmax;
884  if (turbChoice.pbl_mrf_highres_bounds) {
885  constexpr Real ckz = Real(0.001);
886  constexpr Real Kmax = Real(1000.0);
887  rhoKmin = ckz * dz_terrain * rho;
888  rhoKmax = rho * Kmax;
889  } else {
890  constexpr Real Kmin = Real(0.1);
891  constexpr Real Kmax = Real(300.0);
892  rhoKmin = rho * Kmin;
893  rhoKmax = rho * Kmax;
894  }
895 
896  K_turb(i, j, k, EddyDiff::Mom_v) = std::max(
897  std::min(K_turb(i, j, k, EddyDiff::Mom_v), rhoKmax), rhoKmin);
898  K_turb(i, j, k, EddyDiff::Theta_v) = std::max(
899  std::min(K_turb(i, j, k, EddyDiff::Theta_v), rhoKmax), rhoKmin);
900  K_turb(i, j, k, EddyDiff::Q_v) = std::max(
901  std::min(K_turb(i, j, k, EddyDiff::Q_v), rhoKmax), rhoKmin);
902  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = pblh_corr_arr(i, j, 0);
903 
904  if (k < pbli_extent) {
905  K_turb(i, j, k, EddyDiff::HGAMT_v) = hgamt_arr(i, j, 0);
906  K_turb(i, j, k, EddyDiff::HGAMQ_v) = hgamq_arr(i, j, 0);
907  } else {
908  K_turb(i, j, k, EddyDiff::HGAMT_v) = Real(0);
909  K_turb(i, j, k, EddyDiff::HGAMQ_v) = Real(0);
910  }
911  });
912 
913  // FOEXTRAP top and bottom ghost cells
914  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
915  {
916  K_turb(i, j, klo-1, EddyDiff::Mom_v ) = K_turb(i, j, klo, EddyDiff::Mom_v );
917  K_turb(i, j, klo-1, EddyDiff::Theta_v) = K_turb(i, j, klo, EddyDiff::Theta_v);
918  K_turb(i, j, klo-1, EddyDiff::Q_v ) = K_turb(i, j, klo, EddyDiff::Q_v );
919  K_turb(i, j, klo-1, EddyDiff::HGAMT_v) = K_turb(i, j, klo, EddyDiff::HGAMT_v);
920  K_turb(i, j, klo-1, EddyDiff::HGAMQ_v) = K_turb(i, j, klo, EddyDiff::HGAMQ_v);
921  K_turb(i, j, klo-1, EddyDiff::Turb_lengthscale) = K_turb(i, j, klo, EddyDiff::Turb_lengthscale);
922  K_turb(i, j, khi+1, EddyDiff::Mom_v ) = K_turb(i, j, khi, EddyDiff::Mom_v );
923  K_turb(i, j, khi+1, EddyDiff::Theta_v) = K_turb(i, j, khi, EddyDiff::Theta_v);
924  K_turb(i, j, khi+1, EddyDiff::Q_v ) = K_turb(i, j, khi, EddyDiff::Q_v );
925  K_turb(i, j, khi+1, EddyDiff::HGAMT_v) = K_turb(i, j, khi, EddyDiff::HGAMT_v);
926  K_turb(i, j, khi+1, EddyDiff::HGAMQ_v) = K_turb(i, j, khi, EddyDiff::HGAMQ_v);
927  K_turb(i, j, khi+1, EddyDiff::Turb_lengthscale) = K_turb(i, j, khi, EddyDiff::Turb_lengthscale);
928  });
929  }// mfi
930 }
epsv
constexpr amrex::Real epsv
Definition: ERF_Constants.H:53
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:63
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
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
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
pbl_kh_blend_and_cap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pbl_kh_blend_and_cap(amrex::Real K_h, amrex::Real dx, amrex::Real L_blend, amrex::Real C_s, amrex::Real c_max, amrex::Real SmnSmn, bool use_smag) noexcept
Apply scale-aware blending and ceiling to a single K_h value.
Definition: ERF_PBLScaleAwareBlending.H:132
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::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
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
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_blend_length
amrex::Real pbl_blend_length
Boutle blending length L [m]. 0 = off.
Definition: ERF_TurbStruct.H:722
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_mrf_cloud_adjustment
bool enable_mrf_cloud_adjustment
Whether MRF cloud-aware stability adjustments are enabled.
Definition: ERF_TurbStruct.H:728
TurbChoice::pbl_mrf_const_b
amrex::Real pbl_mrf_const_b
MRF constant used to compute PBL height.
Definition: ERF_TurbStruct.H:717
TurbChoice::mrf_moistvars
bool mrf_moistvars
Whether MRF applies turbulence to moisture variables.
Definition: ERF_TurbStruct.H:726
TurbChoice::pbl_blend_cs
amrex::Real pbl_blend_cs
Smagorinsky coeff for K_h ceiling.
Definition: ERF_TurbStruct.H:723
TurbChoice::pbl_blend_c_max
amrex::Real pbl_blend_c_max
Power-law ceiling coeff [m^(2/3)/s].
Definition: ERF_TurbStruct.H:724
TurbChoice::enable_mrf_countergradient
bool enable_mrf_countergradient
Whether MRF countergradient corrections are enabled.
Definition: ERF_TurbStruct.H:727
TurbChoice::pbl_mrf_highres_bounds
bool pbl_mrf_highres_bounds
Whether MRF applies high-resolution grid-dependent diffusivity bounds.
Definition: ERF_TurbStruct.H:729
TurbChoice::pbl_mrf_Ribcr
amrex::Real pbl_mrf_Ribcr
Critical bulk Richardson number for the MRF PBL scheme.
Definition: ERF_TurbStruct.H:716
TurbChoice::pbl_mrf_sf
amrex::Real pbl_mrf_sf
MRF surface flux value used to compute PBL height.
Definition: ERF_TurbStruct.H:718

Referenced by ComputeTurbulentViscosity().

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