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)
 
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.0 → full mixing (correct!)
  • In stable regions: ∂θ_v/∂z > 0 → Ri > 0 → f(Ri) < 1.0 → reduced mixing
  • In neutral regions: Ri ≈ 0 → f(Ri) ≈ 1.0 → 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]i,j,kcell indices
[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²]

131 {
132  const int rho_comp = Rho_comp;
133  const int rhoTheta_comp = RhoTheta_comp;
134  const int qv_index = moisture_indices.qv;
135  const int qc_index = moisture_indices.qc;
136 
137  amrex::Real rho = cell_data(i, j, k, rho_comp);
138 
139  bool saturated = IsSaturated(i, j, k, cell_data, qc_index);
140 
141  amrex::Real inv_theta;
142  amrex::Real dtheta_dz;
143 
144  if (saturated && qc_index >= 0) {
145  amrex::Real theta_l = GetThetal(i, j, k , cell_data, moisture_indices);
146  amrex::Real theta_l_kp1= GetThetal(i, j, k+1, cell_data, moisture_indices);
147  amrex::Real theta_l_km1= GetThetal(i, j, k-1, cell_data, moisture_indices);
148 
149  inv_theta = 1.0 / theta_l;
150  dtheta_dz = 0.5 * (theta_l_kp1 - theta_l_km1) * dzInv;
151 
152  } else if (qv_index >= 0) {
153  amrex::Real theta_v = GetThetav(i, j, k , cell_data, moisture_indices);
154  amrex::Real theta_v_kp1= GetThetav(i, j, k+1, cell_data, moisture_indices);
155  amrex::Real theta_v_km1= GetThetav(i, j, k-1, cell_data, moisture_indices);
156 
157  inv_theta = 1.0 / theta_v;
158  dtheta_dz = 0.5 * (theta_v_kp1 - theta_v_km1) * dzInv;
159 
160  } else {
161  amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
162  amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
163  amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
164 
165  inv_theta = 1.0 / theta;
166  dtheta_dz = 0.5 * (theta_kp1 - theta_km1) * dzInv;
167  }
168 
169  return const_grav * inv_theta * dtheta_dz;
170 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
rho
rho
Definition: ERF_InitCustomPert_Bubble.H:106
GetThetal
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetal(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:101
GetThetav
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:72
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:80
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
LsmVar_MM5::theta
@ theta
Definition: ERF_MM5.H:20
MoistureComponentIndices::qv
int qv
Definition: ERF_DataStruct.H:107
MoistureComponentIndices::qc
int qc
Definition: ERF_DataStruct.H:108

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 
)
181 {
182  const int rho_comp = Rho_comp;
183  const int rhoTheta_comp = RhoTheta_comp;
184  const int qv_index = moisture_indices.qv;
185  const int qc_index = moisture_indices.qc;
186 
187  amrex::Real rho = cell_data(i, j, k, rho_comp);
188 
189  bool saturated = IsSaturated(i, j, k, cell_data, qc_index);
190 
191  amrex::Real inv_theta;
192  amrex::Real dtheta_dz;
193 
194  bool is_covered_kp1 = c_cflag(i,j,k+1).isCovered();
195  bool is_covered_km1 = c_cflag(i,j,k-1).isCovered();
196 
197  if (saturated && qc_index >= 0) {
198  amrex::Real theta_l = GetThetal(i, j, k , cell_data, moisture_indices);
199  amrex::Real theta_l_kp1= GetThetal(i, j, k+1, cell_data, moisture_indices);
200  amrex::Real theta_l_km1= GetThetal(i, j, k-1, cell_data, moisture_indices);
201 
202  inv_theta = 1.0 / theta_l;
203 
204  if (is_covered_kp1) {
205  amrex::Real theta_l_km2= GetThetal(i, j, k-2, cell_data, moisture_indices);
206  dtheta_dz = (3.*theta_l -4.*theta_l_km1 + theta_l_km2) * 0.5 * dzInv;
207  } else if (is_covered_km1) {
208  amrex::Real theta_l_kp2= GetThetal(i, j, k+2, cell_data, moisture_indices);
209  dtheta_dz = (-theta_l_kp2 + 4.*theta_l_kp1 - 3.*theta_l) * 0.5 * dzInv;
210  } else {
211  dtheta_dz = 0.5 * (theta_l_kp1 - theta_l_km1) * dzInv;
212  }
213 
214  } else if (qv_index >= 0) {
215  amrex::Real theta_v = GetThetav(i, j, k , cell_data, moisture_indices);
216  amrex::Real theta_v_kp1= GetThetav(i, j, k+1, cell_data, moisture_indices);
217  amrex::Real theta_v_km1= GetThetav(i, j, k-1, cell_data, moisture_indices);
218 
219  inv_theta = 1.0 / theta_v;
220 
221  if (is_covered_kp1) {
222  amrex::Real theta_v_km2= GetThetav(i, j, k-2, cell_data, moisture_indices);
223  dtheta_dz = (3.*theta_v -4.*theta_v_km1 + theta_v_km2) * 0.5 * dzInv;
224  } else if (is_covered_km1) {
225  amrex::Real theta_v_kp2= GetThetav(i, j, k+2, cell_data, moisture_indices);
226  dtheta_dz = (-theta_v_kp2 + 4.*theta_v_kp1 - 3.*theta_v) * 0.5 * dzInv;
227  } else {
228  dtheta_dz = 0.5 * (theta_v_kp1 - theta_v_km1) * dzInv;
229  }
230 
231  } else {
232  amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
233  amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
234  amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
235 
236  inv_theta = 1.0 / theta;
237 
238  if (is_covered_kp1) {
239  amrex::Real theta_km2 = cell_data(i, j, k-2, rhoTheta_comp) / cell_data(i, j, k-2, rho_comp);
240  dtheta_dz = (3.*theta - 4.*theta_km1 + theta_km2) * 0.5 * dzInv;
241  } else if (is_covered_km1) {
242  amrex::Real theta_kp2 = cell_data(i, j, k+2, rhoTheta_comp) / cell_data(i, j, k+2, rho_comp);
243  dtheta_dz = (-theta_kp2 + 4.*theta_kp1 - 3.*theta) * 0.5 * dzInv;
244  } else {
245  dtheta_dz = 0.5 * (theta_kp1 - theta_km1) * dzInv;
246  }
247  }
248 
249  return const_grav * inv_theta * dtheta_dz;
250 }

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]
371 {
372  amrex::Real S2_safe = amrex::max(S2_vert, 1.0e-10);
373  return N2_moist / S2_safe;
374 }

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]i,j,kcell indices
[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} = 0.5 * (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²]

285 {
286  amrex::Real u_c_km1 = 0.5 * (u(i, j, k-1) + u(i+1, j, k-1));
287  amrex::Real u_c_kp1 = 0.5 * (u(i, j, k+1) + u(i+1, j, k+1));
288 
289  amrex::Real v_c_km1 = 0.5 * (v(i, j, k-1) + v(i, j+1, k-1));
290  amrex::Real v_c_kp1 = 0.5 * (v(i, j, k+1) + v(i, j+1, k+1));
291 
292  amrex::Real du_dz = 0.5 * (u_c_kp1 - u_c_km1) * dzInv;
293  amrex::Real dv_dz = 0.5 * (v_c_kp1 - v_c_km1) * dzInv;
294 
295  amrex::Real S2_vert = du_dz * du_dz + dv_dz * dv_dz;
296 
297  return S2_vert;
298 }

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 
)
310 {
311  amrex::Real du_dz;
312  amrex::Real dv_dz;
313 
314  amrex::Real u_c = ( u_vfrac(i,j,k) * u(i, j, k) + u_vfrac(i+1,j,k) * u(i+1, j, k) )
315  / ( u_vfrac(i,j,k) + u_vfrac(i+1,j,k) );
316  amrex::Real v_c = ( v_vfrac(i,j,k) * v(i, j, k) + v_vfrac(i,j+1,k) * v(i, j+1, k) )
317  / ( v_vfrac(i,j,k) + v_vfrac(i,j+1,k) );
318 
319  if (c_cflag(i,j,k+1).isCovered()) {
320  amrex::Real u_c_km1 = ( u_vfrac(i,j,k-1) * u(i, j, k-1) + u_vfrac(i+1,j,k-1) * u(i+1, j, k-1) )
321  / ( u_vfrac(i,j,k-1) + u_vfrac(i+1,j,k-1) );
322  amrex::Real u_c_km2 = ( u_vfrac(i,j,k-2) * u(i, j, k-2) + u_vfrac(i+1,j,k-2) * u(i+1, j, k-2) )
323  / ( u_vfrac(i,j,k-2) + u_vfrac(i+1,j,k-2) );
324  amrex::Real v_c_km1 = ( v_vfrac(i,j,k-1) * v(i, j, k-1) + v_vfrac(i,j+1,k-1) * v(i, j+1, k-1) )
325  / ( v_vfrac(i,j,k-1) + v_vfrac(i,j+1,k-1) );
326  amrex::Real v_c_km2 = ( v_vfrac(i,j,k-2) * v(i, j, k-2) + v_vfrac(i,j+1,k-2) * v(i, j+1, k-2) )
327  / ( v_vfrac(i,j,k-2) + v_vfrac(i,j+1,k-2) );
328  du_dz = (3.*u_c - 4.*u_c_km1 + u_c_km2) * 0.5 * dzInv;
329  dv_dz = (3.*v_c - 4.*v_c_km1 + v_c_km2) * 0.5 * dzInv;
330 
331  } else if (c_cflag(i,j,k-1).isCovered()) {
332  amrex::Real u_c_kp1 = ( u_vfrac(i,j,k+1) * u(i, j, k+1) + u_vfrac(i+1,j,k+1) * u(i+1, j, k+1) )
333  / ( u_vfrac(i,j,k+1) + u_vfrac(i+1,j,k+1) );
334  amrex::Real u_c_kp2 = ( u_vfrac(i,j,k+2) * u(i, j, k+2) + u_vfrac(i+1,j,k+2) * u(i+1, j, k+2) )
335  / ( u_vfrac(i,j,k+2) + u_vfrac(i+1,j,k+2) );
336  amrex::Real v_c_kp1 = ( v_vfrac(i,j,k+1) * v(i, j, k+1) + v_vfrac(i,j+1,k+1) * v(i, j+1, k+1) )
337  / ( v_vfrac(i,j,k+1) + v_vfrac(i,j+1,k+1) );
338  amrex::Real v_c_kp2 = ( v_vfrac(i,j,k+2) * v(i, j, k+2) + v_vfrac(i,j+1,k+2) * v(i, j+1, k+2) )
339  / ( v_vfrac(i,j,k+2) + v_vfrac(i,j+1,k+2) );
340  du_dz = (-u_c_kp2 + 4.*u_c_kp1 - 3.*u_c) * 0.5 * dzInv;
341  dv_dz = (-v_c_kp2 + 4.*v_c_kp1 - 3.*v_c) * 0.5 * dzInv;
342  } else {
343  amrex::Real u_c_km1 = ( u_vfrac(i,j,k-1) * u(i, j, k-1) + u_vfrac(i+1,j,k-1) * u(i+1, j, k-1) )
344  / ( u_vfrac(i,j,k-1) + u_vfrac(i+1,j,k-1) );
345  amrex::Real u_c_kp1 = ( u_vfrac(i,j,k+1) * u(i, j, k+1) + u_vfrac(i+1,j,k+1) * u(i+1, j, k+1) )
346  / ( u_vfrac(i,j,k+1) + u_vfrac(i+1,j,k+1) );
347  amrex::Real v_c_km1 = ( v_vfrac(i,j,k-1) * v(i, j, k-1) + v_vfrac(i,j+1,k-1) * v(i, j+1, k-1) )
348  / ( v_vfrac(i,j,k-1) + v_vfrac(i,j+1,k-1) );
349  amrex::Real v_c_kp1 = ( v_vfrac(i,j,k+1) * v(i, j, k+1) + v_vfrac(i,j+1,k+1) * v(i, j+1, k+1) )
350  / ( v_vfrac(i,j,k+1) + v_vfrac(i,j+1,k+1) );
351  du_dz = 0.5 * (u_c_kp1 - u_c_km1) * dzInv;
352  dv_dz = 0.5 * (v_c_kp1 - v_c_km1) * dzInv;
353  }
354 
355  amrex::Real S2_vert = du_dz * du_dz + dv_dz * dv_dz;
356 
357  return S2_vert;
358 }

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]i,j,kcell indices
[in]cell_dataconserved state (ρ, ρθ, ρq_v, ρq_c, etc.)
[in]moisture_indicesindices for moisture species
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.

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

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]
387 {
388  if (Ri >= Ri_crit) {
389  return 0.0;
390  } else if (Ri > 0.0) {
391  return std::sqrt(1.0 - Ri / Ri_crit);
392  } else {
393  return 1.0;
394  }
395 }

Referenced by ComputeTurbulentViscosityLES(), and ComputeTurbulentViscosityLES_EB().

Here is the caller graph for this function: