ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_RichardsonNumber.H File Reference
#include <AMReX_GpuQualifiers.H>
#include <AMReX_REAL.H>
#include <cmath>
#include "ERF_Constants.H"
#include "ERF_IndexDefines.H"
#include "ERF_DataStruct.H"
#include "ERF_MoistUtils.H"
Include dependency graph for ERF_RichardsonNumber.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

AMREX_GPU_DEVICE AMREX_FORCE_INLINE bool IsSaturated (int i, int j, int k, const amrex::Array4< const amrex::Real > &cell_data, int qc_index)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeN2 (int i, int j, int k, amrex::Real dzInv, amrex::Real const_grav, const amrex::Array4< const amrex::Real > &cell_data, const MoistureComponentIndices &moisture_indices)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeN2_EB (int i, int j, int k, const amrex::Array4< const amrex::EBCellFlag > &c_cflag, amrex::Real dzInv, amrex::Real const_grav, const amrex::Array4< const amrex::Real > &cell_data, const MoistureComponentIndices &moisture_indices)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeVerticalShear2 (int i, int j, int k, amrex::Real dzInv, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeVerticalShear2_EB (int i, int j, int k, const amrex::Array4< const amrex::EBCellFlag > &c_cflag, const amrex::Array4< const amrex::Real > &u_vfrac, const amrex::Array4< const amrex::Real > &v_vfrac, amrex::Real dzInv, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, bool &shear_is_valid)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeRichardson (amrex::Real N2_moist, amrex::Real S2_vert)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real StabilityFunction (amrex::Real Ri, amrex::Real Ri_crit)
 

Detailed Description

Functions for computing moist Richardson number and Mason (1989) stability function for Smagorinsky turbulence closure with conditional instability correction.

PHYSICAL BACKGROUND:

In cloud-resolving simulations, standard Smagorinsky models overpredict turbulent mixing in conditionally unstable regions (cloud updrafts) where buoyancy-driven instability occurs. This module implements the moist Richardson number correction following Mason (1989) and Stevens et al. (2005) for marine stratocumulus.

KEY EQUATIONS:

1) Moist Brunt-Väisälä frequency: N²_moist = (g/θ_v) * ∂θ_v/∂z (unsaturated) N²_moist = (g/θ_l) * ∂θ_l/∂z (saturated, in clouds)

2) Moist Richardson number: Ri_moist = N²_moist / S²_vert where S²_vert = (∂u/∂z)² + (∂v/∂z)²

3) Mason (1989) stability function: f(Ri) = 0 if Ri ≥ Ri_crit (complete suppression) = sqrt(1 - Ri/Ri_crit) if 0 < Ri < Ri_crit (partial suppression) = 1 if Ri ≤ 0 (no suppression, unstable/neutral)

4) Corrected vertical eddy viscosity: μ_t,v = ρ(C_s·Δ)²|S| · f(Ri_moist)

PHYSICAL INTERPRETATION:

  • In cloud updrafts: ∂θ_l/∂z < 0 → Ri < 0 → f(Ri) = 1 → full mixing (correct!)
  • In stable regions: ∂θ_v/∂z > 0 → Ri > 0 → f(Ri) < 1 → reduced mixing
  • In neutral regions: Ri ≈ 0 → f(Ri) ≈ 1 → standard Smagorinsky

CUDA-CONSERVATIVE DESIGN:

  • All functions are explicit device functions (not lambdas)
  • Use simple scalar parameters (no struct captures in device context)
  • Geometry passed via local scalars (dzInv), not complex objects
  • Reuse existing ERF GetThetav/GetThetal from ERF_MoistUtils.H

REFERENCES:

  • Mason, P. J. (1989): Large-eddy simulation of the convective atmospheric boundary layer. J. Atmos. Sci., 46, 1492-1516.
  • Stevens et al. (2005): Evaluation of large-eddy simulations via observations of nocturnal marine stratocumulus. Mon. Wea. Rev., 133, 1443-1462.
  • Smagorinsky, J. (1963): General circulation experiments with the primitive equations. Mon. Wea. Rev., 91, 99-164.

Function Documentation

◆ ComputeN2()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeN2 ( int  i,
int  j,
int  k,
amrex::Real  dzInv,
amrex::Real  const_grav,
const amrex::Array4< const amrex::Real > &  cell_data,
const MoistureComponentIndices &  moisture_indices 
)

Compute moist Brunt-Väisälä frequency squared (N²_moist).

Parameters
[in]icell i-index
[in]jcell j-index
[in]kcell k-index
[in]dzInvinverse vertical grid spacing (1/Δz) [1/m]
[in]const_gravgravitational acceleration g [m/s²]
[in]cell_dataconserved state
[in]moisture_indicesmoisture species indices
Returns
N²_moist [1/s²]

PHYSICS:

Uses centered finite differences: ∂/∂z ≈ (f_{k+1} - f_{k-1}) / (2Δz)

Unsaturated: N² = (g/θ_v) · ∂θ_v/∂z where θ_v = θ(1 + 0.61q_v - q_c) accounts for vapor buoyancy and condensate loading

Saturated: N² = (g/θ_l) · ∂θ_l/∂z where θ_l = θ - (L_v/c_p·Π)·q_c is the liquid water potential temperature

CONDITIONAL INSTABILITY:

In cloud updrafts, ∂θ_l/∂z can be negative due to release of latent heat during ascent. This gives N² < 0, which is physically correct for convective instability and should enhance (not suppress) mixing.

UNITS:

Input: dzInv [1/m], const_grav [m/s²], θ [K] Output: N² [1/s²]

135 {
136  const int rho_comp = Rho_comp;
137  const int rhoTheta_comp = RhoTheta_comp;
138  const int qv_index = moisture_indices.qv;
139  const int qc_index = moisture_indices.qc;
140 
141  amrex::Real rho = cell_data(i, j, k, rho_comp);
142 
143  bool saturated = IsSaturated(i, j, k, cell_data, qc_index);
144 
145  amrex::Real inv_theta;
146  amrex::Real dtheta_dz;
147 
148  if (saturated && qc_index >= 0) {
149  amrex::Real theta_l = GetThetal(i, j, k , cell_data, moisture_indices);
150  amrex::Real theta_l_kp1= GetThetal(i, j, k+1, cell_data, moisture_indices);
151  amrex::Real theta_l_km1= GetThetal(i, j, k-1, cell_data, moisture_indices);
152 
153  inv_theta = one / theta_l;
154  dtheta_dz = myhalf * (theta_l_kp1 - theta_l_km1) * dzInv;
155 
156  } else if (qv_index >= 0) {
157  amrex::Real theta_v = GetThetav(i, j, k , cell_data, moisture_indices);
158  amrex::Real theta_v_kp1= GetThetav(i, j, k+1, cell_data, moisture_indices);
159  amrex::Real theta_v_km1= GetThetav(i, j, k-1, cell_data, moisture_indices);
160 
161  inv_theta = one / theta_v;
162  dtheta_dz = myhalf * (theta_v_kp1 - theta_v_km1) * dzInv;
163 
164  } else {
165  amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
166  amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
167  amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
168 
169  inv_theta = one / theta;
170  dtheta_dz = myhalf * (theta_kp1 - theta_km1) * dzInv;
171  }
172 
173  return const_grav * inv_theta * dtheta_dz;
174 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:42
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:43
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:74
GetThetal
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetal(const int &i, const int &j, const int &k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:130
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
IsSaturated
AMREX_GPU_DEVICE AMREX_FORCE_INLINE bool IsSaturated(int i, int j, int k, const amrex::Array4< const amrex::Real > &cell_data, int qc_index)
Definition: ERF_RichardsonNumber.H:82
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:26
MoistureComponentIndices::qv
int qv
water vapor
Definition: ERF_DataStruct.H:249
MoistureComponentIndices::qc
int qc
cloud liquid water
Definition: ERF_DataStruct.H:250

Referenced by ComputeTurbulentViscosityLES(), and ComputeTurbulentViscosityLES_EB().

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

◆ ComputeN2_EB()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeN2_EB ( int  i,
int  j,
int  k,
const amrex::Array4< const amrex::EBCellFlag > &  c_cflag,
amrex::Real  dzInv,
amrex::Real  const_grav,
const amrex::Array4< const amrex::Real > &  cell_data,
const MoistureComponentIndices &  moisture_indices 
)

Compute moist Brunt-Vaisala frequency squared with embedded-boundary one-sided differences.

Parameters
[in]icell i-index
[in]jcell j-index
[in]kcell k-index
[in]c_cflagcell-centered EB flags
[in]dzInvinverse vertical grid spacing
[in]const_gravgravitational acceleration
[in]cell_dataconserved state
[in]moisture_indicesmoisture species indices
Returns
moist Brunt-Vaisala frequency squared
198 {
199  const int rho_comp = Rho_comp;
200  const int rhoTheta_comp = RhoTheta_comp;
201  const int qv_index = moisture_indices.qv;
202  const int qc_index = moisture_indices.qc;
203 
204  amrex::Real rho = cell_data(i, j, k, rho_comp);
205 
206  bool saturated = IsSaturated(i, j, k, cell_data, qc_index);
207 
208  amrex::Real inv_theta;
209  amrex::Real dtheta_dz;
210 
211  const bool is_covered_kp1 = c_cflag(i,j,k+1).isCovered();
212  const bool is_covered_km1 = c_cflag(i,j,k-1).isCovered();
213 
214  if (saturated && qc_index >= 0) {
215 
216  const amrex::Real theta_l = GetThetal(i, j, k , cell_data, moisture_indices);
217 
218  inv_theta = one / theta_l;
219 
220  if (is_covered_kp1) {
221 
222  if (is_covered_km1) {
223  return zero;
224  }
225 
226  const amrex::Real theta_l_km1 = GetThetal(i, j, k-1, cell_data, moisture_indices);
227 
228  if (c_cflag(i,j,k-2).isCovered()) {
229  dtheta_dz = (theta_l - theta_l_km1) * dzInv;
230  } else {
231  const amrex::Real theta_l_km2 = GetThetal(i, j, k-2, cell_data, moisture_indices);
232  dtheta_dz = (three*theta_l - four*theta_l_km1 + theta_l_km2) * myhalf * dzInv;
233  }
234 
235  } else if (is_covered_km1) {
236 
237  const amrex::Real theta_l_kp1 = GetThetal(i, j, k+1, cell_data, moisture_indices);
238 
239  if (c_cflag(i,j,k+2).isCovered()) {
240  dtheta_dz = (theta_l_kp1 - theta_l) * dzInv;
241  } else {
242  const amrex::Real theta_l_kp2 = GetThetal(i, j, k+2, cell_data, moisture_indices);
243  dtheta_dz = (-theta_l_kp2 + four*theta_l_kp1 - three*theta_l) * myhalf * dzInv;
244  }
245 
246  } else {
247 
248  const amrex::Real theta_l_kp1 = GetThetal(i, j, k+1, cell_data, moisture_indices);
249  const amrex::Real theta_l_km1 = GetThetal(i, j, k-1, cell_data, moisture_indices);
250  dtheta_dz = (theta_l_kp1 - theta_l_km1) * myhalf * dzInv;
251  }
252 
253  } else if (qv_index >= 0) {
254 
255  const amrex::Real theta_v = GetThetav(i, j, k, cell_data, moisture_indices);
256 
257  inv_theta = one / theta_v;
258 
259  if (is_covered_kp1) {
260 
261  if (is_covered_km1) {
262  return zero;
263  }
264 
265  const amrex::Real theta_v_km1 = GetThetav(i, j, k-1, cell_data, moisture_indices);
266 
267  if (c_cflag(i,j,k-2).isCovered()) {
268  dtheta_dz = (theta_v - theta_v_km1) * dzInv;
269  } else {
270  const amrex::Real theta_v_km2 = GetThetav(i, j, k-2, cell_data, moisture_indices);
271  dtheta_dz = (three*theta_v - four*theta_v_km1 + theta_v_km2) * myhalf * dzInv;
272  }
273 
274  } else if (is_covered_km1) {
275 
276  const amrex::Real theta_v_kp1 = GetThetav(i, j, k+1, cell_data, moisture_indices);
277 
278  if (c_cflag(i,j,k+2).isCovered()) {
279  dtheta_dz = (theta_v_kp1 - theta_v) * dzInv;
280  } else {
281  const amrex::Real theta_v_kp2 = GetThetav(i, j, k+2, cell_data, moisture_indices);
282  dtheta_dz = (-theta_v_kp2 + four*theta_v_kp1 - three*theta_v) * myhalf * dzInv;
283  }
284 
285  } else {
286 
287  const amrex::Real theta_v_kp1 = GetThetav(i, j, k+1, cell_data, moisture_indices);
288  const amrex::Real theta_v_km1 = GetThetav(i, j, k-1, cell_data, moisture_indices);
289  dtheta_dz = (theta_v_kp1 - theta_v_km1) * myhalf * dzInv;
290  }
291 
292  } else {
293 
294  const amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
295 
296  inv_theta = one / theta;
297 
298  if (is_covered_kp1) {
299 
300  if (is_covered_km1) {
301  return zero;
302  }
303 
304  const amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
305 
306  if (c_cflag(i,j,k-2).isCovered()) {
307  dtheta_dz = (theta - theta_km1) * dzInv;
308  } else {
309  const amrex::Real theta_km2 = cell_data(i, j, k-2, rhoTheta_comp) / cell_data(i, j, k-2, rho_comp);
310  dtheta_dz = (three*theta - four*theta_km1 + theta_km2) * myhalf * dzInv;
311  }
312 
313  } else if (is_covered_km1) {
314 
315  const amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
316 
317  if (c_cflag(i,j,k+2).isCovered()) {
318  dtheta_dz = (theta_kp1 - theta) * dzInv;
319  } else {
320  const amrex::Real theta_kp2 = cell_data(i, j, k+2, rhoTheta_comp) / cell_data(i, j, k+2, rho_comp);
321  dtheta_dz = (-theta_kp2 + four*theta_kp1 - three*theta) * myhalf * dzInv;
322  }
323 
324  } else {
325 
326  const amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
327  const amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
328  dtheta_dz = (theta_kp1 - theta_km1) * myhalf * dzInv;
329  }
330  }
331 
332  return const_grav * inv_theta * dtheta_dz;
333 }
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:39
four
constexpr amrex::Real four
Definition: ERF_NumericalConstants.H:40
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36

Referenced by ComputeTurbulentViscosityLES_EB().

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

◆ ComputeRichardson()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeRichardson ( amrex::Real  N2_moist,
amrex::Real  S2_vert 
)

Compute moist Richardson number.

Parameters
[in]N2_moistmoist Brunt-Väisälä frequency squared [1/s²]
[in]S2_vertvertical wind shear squared [1/s²]
Returns
Ri_moist = N²_moist / S²_vert [dimensionless]
507 {
508  amrex::Real S2_safe = amrex::max(S2_vert, amrex::Real(1.0e-10));
509  return N2_moist / S2_safe;
510 }

Referenced by ComputeTurbulentViscosityLES(), and ComputeTurbulentViscosityLES_EB().

Here is the caller graph for this function:

◆ ComputeVerticalShear2()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeVerticalShear2 ( int  i,
int  j,
int  k,
amrex::Real  dzInv,
const amrex::Array4< const amrex::Real > &  u,
const amrex::Array4< const amrex::Real > &  v 
)

Compute vertical shear squared (S²_vert).

Parameters
[in]icell i-index
[in]jcell j-index
[in]kcell k-index
[in]dzInvinverse vertical grid spacing (1/Δz) [1/m]
[in]uface-centered x-velocity (u) array with ngrow≥1 in z
[in]vface-centered y-velocity (v) array with ngrow≥1 in z
Returns
S²_vert = (∂u/∂z)² + (∂v/∂z)² [1/s²]

DISCRETIZATION:

ERF uses Arakawa C-grid: u defined at (i+1/2, j, k) (x-faces) v defined at (i, j+1/2, k) (y-faces) w defined at (i, j, k+1/2) (z-faces) cell centers at (i, j, k)

To get shear at cell center (i,j,k), interpolate velocities to center then difference: ∂u/∂z|_{i,j,k} ≈ (u_{i,j,k+1} - u_{i,j,k-1}) / (2Δz) where u_{i,j,k} = myhalf * (u_{i-1/2,j,k} + u_{i+1/2,j,k})

UNITS:

Input: dzInv [1/m], u, v [m/s] Output: S²_vert [1/s²]

370 {
371  amrex::Real u_c_km1 = myhalf * (u(i, j, k-1) + u(i+1, j, k-1));
372  amrex::Real u_c_kp1 = myhalf * (u(i, j, k+1) + u(i+1, j, k+1));
373 
374  amrex::Real v_c_km1 = myhalf * (v(i, j, k-1) + v(i, j+1, k-1));
375  amrex::Real v_c_kp1 = myhalf * (v(i, j, k+1) + v(i, j+1, k+1));
376 
377  amrex::Real du_dz = myhalf * (u_c_kp1 - u_c_km1) * dzInv;
378  amrex::Real dv_dz = myhalf * (v_c_kp1 - v_c_km1) * dzInv;
379 
380  amrex::Real S2_vert = du_dz * du_dz + dv_dz * dv_dz;
381 
382  return S2_vert;
383 }

Referenced by ComputeTurbulentViscosityLES(), and ComputeTurbulentViscosityLES_EB().

Here is the caller graph for this function:

◆ ComputeVerticalShear2_EB()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeVerticalShear2_EB ( int  i,
int  j,
int  k,
const amrex::Array4< const amrex::EBCellFlag > &  c_cflag,
const amrex::Array4< const amrex::Real > &  u_vfrac,
const amrex::Array4< const amrex::Real > &  v_vfrac,
amrex::Real  dzInv,
const amrex::Array4< const amrex::Real > &  u,
const amrex::Array4< const amrex::Real > &  v,
bool &  shear_is_valid 
)

Compute vertical shear squared with embedded-boundary one-sided differences.

Parameters
[in]icell i-index
[in]jcell j-index
[in]kcell k-index
[in]c_cflagcell-centered EB flags
[in]u_vfracx-face volume fractions
[in]v_vfracy-face volume fractions
[in]dzInvinverse vertical grid spacing
[in]uface-centered x-velocity
[in]vface-centered y-velocity
[out]shear_is_validwhether the required shear stencil is available
Returns
vertical shear squared
411 {
412  // Volume-fraction weighted reconstruction of a velocity component at the
413  // centre of cell (ii,jj,kk) from the two staggered faces that bracket it in
414  // the (di,dj) direction. Returns false when neither face carries any volume,
415  // in which case there is no cell-centred value to be had.
416  auto cell_avg = [=] (const amrex::Array4<const amrex::Real>& vfrac,
417  const amrex::Array4<const amrex::Real>& vel,
418  int ii, int jj, int kk, int di, int dj,
419  amrex::Real& q) -> bool
420  {
421  const amrex::Real f0 = vfrac(ii, jj, kk);
422  const amrex::Real f1 = vfrac(ii+di, jj+dj, kk);
423  const amrex::Real fsum = f0 + f1;
424  if (fsum <= zero) { return false; }
425  q = ((f0 > zero ? f0 * vel(ii, jj, kk) : zero) +
426  (f1 > zero ? f1 * vel(ii+di, jj+dj, kk) : zero)) / fsum;
427  return true;
428  };
429 
430  // Vertical derivative of one velocity component at cell (i,j,k), using the
431  // best stencil the EB geometry leaves available: centred when both k-1 and
432  // k+1 are usable, otherwise one-sided towards whichever neighbour is (second
433  // order when k-2 / k+2 is usable as well). Returns false only when no stencil
434  // at all can be formed, so a single unusable neighbour degrades the order of
435  // accuracy rather than discarding the correction.
436  auto vert_deriv = [=] (const amrex::Array4<const amrex::Real>& vfrac,
437  const amrex::Array4<const amrex::Real>& vel,
438  int di, int dj, amrex::Real& dq_dz) -> bool
439  {
440  amrex::Real c_km1 = zero, c_kp1 = zero;
441  const bool have_km1 = !c_cflag(i,j,k-1).isCovered() &&
442  cell_avg(vfrac, vel, i, j, k-1, di, dj, c_km1);
443  const bool have_kp1 = !c_cflag(i,j,k+1).isCovered() &&
444  cell_avg(vfrac, vel, i, j, k+1, di, dj, c_kp1);
445 
446  // The centred difference never reads the centre cell, so it stays valid
447  // even when the centre cell's own faces carry no volume.
448  if (have_km1 && have_kp1) {
449  dq_dz = myhalf * (c_kp1 - c_km1) * dzInv;
450  return true;
451  }
452 
453  // Every remaining stencil is one-sided, and those do need the centre cell.
454  amrex::Real c_c;
455  if (!cell_avg(vfrac, vel, i, j, k, di, dj, c_c)) { return false; }
456 
457  if (have_km1) {
458  amrex::Real c_km2;
459  if (!c_cflag(i,j,k-2).isCovered() &&
460  cell_avg(vfrac, vel, i, j, k-2, di, dj, c_km2)) {
461  dq_dz = (three * c_c - four * c_km1 + c_km2) * myhalf * dzInv;
462  } else {
463  dq_dz = (c_c - c_km1) * dzInv;
464  }
465  return true;
466  }
467 
468  if (have_kp1) {
469  amrex::Real c_kp2;
470  if (!c_cflag(i,j,k+2).isCovered() &&
471  cell_avg(vfrac, vel, i, j, k+2, di, dj, c_kp2)) {
472  dq_dz = (-c_kp2 + four * c_kp1 - three * c_c) * myhalf * dzInv;
473  } else {
474  dq_dz = (c_kp1 - c_c) * dzInv;
475  }
476  return true;
477  }
478 
479  // Both vertical neighbours are unusable: no stencil exists at all.
480  return false;
481  };
482 
483  amrex::Real du_dz, dv_dz;
484  if (!vert_deriv(u_vfrac, u, 1, 0, du_dz) ||
485  !vert_deriv(v_vfrac, v, 0, 1, dv_dz)) {
486  shear_is_valid = false;
487  return zero;
488  }
489  shear_is_valid = true;
490 
491  amrex::Real S2_vert = du_dz * du_dz + dv_dz * dv_dz;
492 
493  return S2_vert;
494 }
WSM6Ind::q
@ q
Definition: ERF_WSM6.H:273
module_mp_morr_two_moment::di
real(c_double), private di
Definition: ERF_module_mp_morr_two_moment.F90:203

Referenced by ComputeTurbulentViscosityLES_EB().

Here is the caller graph for this function:

◆ IsSaturated()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE bool IsSaturated ( int  i,
int  j,
int  k,
const amrex::Array4< const amrex::Real > &  cell_data,
int  qc_index 
)

Check if grid cell is saturated based on cloud liquid water content.

Parameters
[in]icell i-index
[in]jcell j-index
[in]kcell k-index
[in]cell_dataconserved state (ρ, ρθ, ρq_v, ρq_c, etc.)
[in]qc_indexconserved cloud-water component index
Returns
true if saturated (q_c > 1e-8 kg/kg), false otherwise

PHYSICAL JUSTIFICATION: Threshold of 1e-8 kg/kg (0.01 g/kg) is typical for distinguishing cloudy from clear regions in LES. Below this threshold, condensate is negligible and air is effectively unsaturated.

85 {
86  if (qc_index >= 0) {
87  amrex::Real rho = cell_data(i, j, k, Rho_comp);
88  amrex::Real qc = cell_data(i, j, k, qc_index) / rho;
89  return (qc > amrex::Real(1.0e-8)); // 0.01 g/kg threshold
90  }
91  return false;
92 }
MicVar_SatAdj::qc
@ qc
Definition: ERF_SatAdj.H:42

Referenced by ComputeN2(), and ComputeN2_EB().

Here is the caller graph for this function:

◆ StabilityFunction()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real StabilityFunction ( amrex::Real  Ri,
amrex::Real  Ri_crit 
)

Mason (1989) stability function for turbulent mixing suppression.

Parameters
[in]Rigradient Richardson number [dimensionless]
[in]Ri_critcritical Richardson number (default 0.25)
Returns
f(Ri) ∈ [0, 1] [dimensionless]
523 {
524  amrex::Real stability_factor = amrex::Real(0.);
525  if (Ri >= Ri_crit) {
526  stability_factor = amrex::Real(0.);
527  } else if (Ri > amrex::Real(0.) && Ri < Ri_crit) {
528  stability_factor = std::sqrt(amrex::max(amrex::Real(0.), amrex::Real(1.) - Ri / Ri_crit));
529  } else {
530  stability_factor = amrex::Real(1.);
531  }
532  return stability_factor;
533 }

Referenced by ComputeTurbulentViscosityLES(), and ComputeTurbulentViscosityLES_EB().

Here is the caller graph for this function: