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

Referenced by ComputeTurbulentViscosity().

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