ERF_RichardsonNumber.H

Go to the documentation of this file.

/** \file ERF_RichardsonNumber.H
 *
 * 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.
 */
 
#ifndef ERF_RICHARDSON_NUMBER_H_
#define ERF_RICHARDSON_NUMBER_H_
 
#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"
 
/**
 * Check if grid cell is saturated based on cloud liquid water content.
 *
 * @param[in] i cell i-index
 * @param[in] j cell j-index
 * @param[in] k cell k-index
 * @param[in] cell_data  conserved state (ρ, ρθ, ρq_v, ρq_c, etc.)
 * @param[in] qc_index conserved cloud-water component index
 * @return 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.
 */
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)
{
    if (qc_index >= 0) {
        amrex::Real rho = cell_data(i, j, k, Rho_comp);
        amrex::Real qc = cell_data(i, j, k, qc_index) / rho;
        return (qc > amrex::Real(1.0e-8));  // 0.01 g/kg threshold
    }
    return false;
}
 
/**
 * Compute moist Brunt-Väisälä frequency squared (N²_moist).
 *
 * @param[in] i cell i-index
 * @param[in] j cell j-index
 * @param[in] k cell k-index
 * @param[in] dzInv  inverse vertical grid spacing (1/Δz) [1/m]
 * @param[in] const_grav  gravitational acceleration g [m/s²]
 * @param[in] cell_data  conserved state
 * @param[in] moisture_indices  moisture species indices
 * @return 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²]
 */
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)
{
    const int rho_comp       = Rho_comp;
    const int rhoTheta_comp  = RhoTheta_comp;
    const int qv_index = moisture_indices.qv;
    const int qc_index = moisture_indices.qc;
 
    amrex::Real rho = cell_data(i, j, k, rho_comp);
 
    bool saturated = IsSaturated(i, j, k, cell_data, qc_index);
 
    amrex::Real inv_theta;
    amrex::Real dtheta_dz;
 
    if (saturated && qc_index >= 0) {
        amrex::Real theta_l    = GetThetal(i, j, k  , cell_data, moisture_indices);
        amrex::Real theta_l_kp1= GetThetal(i, j, k+1, cell_data, moisture_indices);
        amrex::Real theta_l_km1= GetThetal(i, j, k-1, cell_data, moisture_indices);
 
        inv_theta = one / theta_l;
        dtheta_dz = myhalf * (theta_l_kp1 - theta_l_km1) * dzInv;
 
    } else if (qv_index >= 0) {
        amrex::Real theta_v    = GetThetav(i, j, k  , cell_data, moisture_indices);
        amrex::Real theta_v_kp1= GetThetav(i, j, k+1, cell_data, moisture_indices);
        amrex::Real theta_v_km1= GetThetav(i, j, k-1, cell_data, moisture_indices);
 
        inv_theta = one / theta_v;
        dtheta_dz = myhalf * (theta_v_kp1 - theta_v_km1) * dzInv;
 
    } else {
        amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
        amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
        amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
 
        inv_theta = one / theta;
        dtheta_dz = myhalf * (theta_kp1 - theta_km1) * dzInv;
    }
 
    return const_grav * inv_theta * dtheta_dz;
}
 
/**
 * Compute moist Brunt-Vaisala frequency squared with embedded-boundary one-sided differences.
 *
 * @param[in] i cell i-index
 * @param[in] j cell j-index
 * @param[in] k cell k-index
 * @param[in] c_cflag cell-centered EB flags
 * @param[in] dzInv inverse vertical grid spacing
 * @param[in] const_grav gravitational acceleration
 * @param[in] cell_data conserved state
 * @param[in] moisture_indices moisture species indices
 * @return moist Brunt-Vaisala frequency squared
 */
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)
{
    const int rho_comp       = Rho_comp;
    const int rhoTheta_comp  = RhoTheta_comp;
    const int qv_index = moisture_indices.qv;
    const int qc_index = moisture_indices.qc;
 
    amrex::Real rho = cell_data(i, j, k, rho_comp);
 
    bool saturated = IsSaturated(i, j, k, cell_data, qc_index);
 
    amrex::Real inv_theta;
    amrex::Real dtheta_dz;
 
    bool is_covered_kp1 = c_cflag(i,j,k+1).isCovered();
    bool is_covered_km1 = c_cflag(i,j,k-1).isCovered();
 
    if (saturated && qc_index >= 0) {
        amrex::Real theta_l    = GetThetal(i, j, k  , cell_data, moisture_indices);
        amrex::Real theta_l_kp1= GetThetal(i, j, k+1, cell_data, moisture_indices);
        amrex::Real theta_l_km1= GetThetal(i, j, k-1, cell_data, moisture_indices);
 
        inv_theta = one / theta_l;
 
        if (is_covered_kp1) {
            amrex::Real theta_l_km2= GetThetal(i, j, k-2, cell_data, moisture_indices);
            dtheta_dz = (three*theta_l - four*theta_l_km1 + theta_l_km2) * myhalf * dzInv;
        } else if (is_covered_km1) {
            amrex::Real theta_l_kp2= GetThetal(i, j, k+2, cell_data, moisture_indices);
            dtheta_dz = (-theta_l_kp2 + four*theta_l_kp1 - three*theta_l) * myhalf * dzInv;
        } else {
            dtheta_dz = myhalf * (theta_l_kp1 - theta_l_km1) * dzInv;
        }
 
    } else if (qv_index >= 0) {
        amrex::Real theta_v    = GetThetav(i, j, k  , cell_data, moisture_indices);
        amrex::Real theta_v_kp1= GetThetav(i, j, k+1, cell_data, moisture_indices);
        amrex::Real theta_v_km1= GetThetav(i, j, k-1, cell_data, moisture_indices);
 
        inv_theta = one / theta_v;
 
        if (is_covered_kp1) {
            amrex::Real theta_v_km2= GetThetav(i, j, k-2, cell_data, moisture_indices);
            dtheta_dz = (three*theta_v - four*theta_v_km1 + theta_v_km2) * myhalf * dzInv;
        } else if (is_covered_km1) {
            amrex::Real theta_v_kp2= GetThetav(i, j, k+2, cell_data, moisture_indices);
            dtheta_dz = (-theta_v_kp2 + four*theta_v_kp1 - three*theta_v) * myhalf * dzInv;
        } else {
            dtheta_dz = myhalf * (theta_v_kp1 - theta_v_km1) * dzInv;
        }
 
    } else {
        amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
        amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
        amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
 
        inv_theta = one / theta;
 
        if (is_covered_kp1) {
            amrex::Real theta_km2 = cell_data(i, j, k-2, rhoTheta_comp) / cell_data(i, j, k-2, rho_comp);
            dtheta_dz = (three*theta - four*theta_km1 + theta_km2) * myhalf * dzInv;
        } else if (is_covered_km1) {
            amrex::Real theta_kp2 = cell_data(i, j, k+2, rhoTheta_comp) / cell_data(i, j, k+2, rho_comp);
            dtheta_dz = (-theta_kp2 + four*theta_kp1 - three*theta) * myhalf * dzInv;
        } else {
            dtheta_dz = myhalf * (theta_kp1 - theta_km1) * dzInv;
        }
    }
 
    return const_grav * inv_theta * dtheta_dz;
}
 
/**
 * Compute vertical shear squared (S²_vert).
 *
 * @param[in] i cell i-index
 * @param[in] j cell j-index
 * @param[in] k cell k-index
 * @param[in] dzInv inverse vertical grid spacing (1/Δz) [1/m]
 * @param[in] u face-centered x-velocity (u) array with ngrow≥1 in z
 * @param[in] v face-centered y-velocity (v) array with ngrow≥1 in z
 * @return 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²]
 */
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::Real u_c_km1 = myhalf * (u(i, j, k-1) + u(i+1, j, k-1));
    amrex::Real u_c_kp1 = myhalf * (u(i, j, k+1) + u(i+1, j, k+1));
 
    amrex::Real v_c_km1 = myhalf * (v(i, j, k-1) + v(i, j+1, k-1));
    amrex::Real v_c_kp1 = myhalf * (v(i, j, k+1) + v(i, j+1, k+1));
 
    amrex::Real du_dz = myhalf * (u_c_kp1 - u_c_km1) * dzInv;
    amrex::Real dv_dz = myhalf * (v_c_kp1 - v_c_km1) * dzInv;
 
    amrex::Real S2_vert = du_dz * du_dz + dv_dz * dv_dz;
 
    return S2_vert;
}
 
/**
 * Compute vertical shear squared with embedded-boundary one-sided differences.
 *
 * @param[in] i cell i-index
 * @param[in] j cell j-index
 * @param[in] k cell k-index
 * @param[in] c_cflag cell-centered EB flags
 * @param[in] u_vfrac x-face volume fractions
 * @param[in] v_vfrac y-face volume fractions
 * @param[in] dzInv inverse vertical grid spacing
 * @param[in] u face-centered x-velocity
 * @param[in] v face-centered y-velocity
 * @return vertical shear squared
 */
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::Real du_dz;
    amrex::Real dv_dz;
 
    amrex::Real u_c = ( u_vfrac(i,j,k) * u(i, j, k) + u_vfrac(i+1,j,k) * u(i+1, j, k) )
                    / ( u_vfrac(i,j,k) + u_vfrac(i+1,j,k) );
    amrex::Real v_c = ( v_vfrac(i,j,k) * v(i, j, k) + v_vfrac(i,j+1,k) * v(i, j+1, k) )
                    / ( v_vfrac(i,j,k) + v_vfrac(i,j+1,k) );
 
    if (c_cflag(i,j,k+1).isCovered()) {
        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) )
                            / ( u_vfrac(i,j,k-1) + u_vfrac(i+1,j,k-1) );
        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) )
                            / ( u_vfrac(i,j,k-2) + u_vfrac(i+1,j,k-2) );
        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) )
                            / ( v_vfrac(i,j,k-1) + v_vfrac(i,j+1,k-1) );
        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) )
                            / ( v_vfrac(i,j,k-2) + v_vfrac(i,j+1,k-2) );
        du_dz = (three*u_c - four*u_c_km1 + u_c_km2) * myhalf * dzInv;
        dv_dz = (three*v_c - four*v_c_km1 + v_c_km2) * myhalf * dzInv;
 
    } else if (c_cflag(i,j,k-1).isCovered()) {
        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) )
                            / ( u_vfrac(i,j,k+1) + u_vfrac(i+1,j,k+1) );
        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) )
                            / ( u_vfrac(i,j,k+2) + u_vfrac(i+1,j,k+2) );
        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) )
                            / ( v_vfrac(i,j,k+1) + v_vfrac(i,j+1,k+1) );
        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) )
                            / ( v_vfrac(i,j,k+2) + v_vfrac(i,j+1,k+2) );
        du_dz = (-u_c_kp2 + four*u_c_kp1 - three*u_c) * myhalf * dzInv;
        dv_dz = (-v_c_kp2 + four*v_c_kp1 - three*v_c) * myhalf * dzInv;
    } else {
        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) )
                            / ( u_vfrac(i,j,k-1) + u_vfrac(i+1,j,k-1) );
        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) )
                            / ( u_vfrac(i,j,k+1) + u_vfrac(i+1,j,k+1) );
        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) )
                            / ( v_vfrac(i,j,k-1) + v_vfrac(i,j+1,k-1) );
        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) )
                            / ( v_vfrac(i,j,k+1) + v_vfrac(i,j+1,k+1) );
        du_dz = myhalf * (u_c_kp1 - u_c_km1) * dzInv;
        dv_dz = myhalf * (v_c_kp1 - v_c_km1) * dzInv;
    }
 
    amrex::Real S2_vert = du_dz * du_dz + dv_dz * dv_dz;
 
    return S2_vert;
}
 
/**
 * Compute moist Richardson number.
 *
 * @param[in] N2_moist  moist Brunt-Väisälä frequency squared [1/s²]
 * @param[in] S2_vert   vertical wind shear squared [1/s²]
 * @return Ri_moist = N²_moist / S²_vert [dimensionless]
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
ComputeRichardson(amrex::Real N2_moist, amrex::Real S2_vert)
{
    amrex::Real S2_safe = amrex::max(S2_vert, amrex::Real(1.0e-10));
    return N2_moist / S2_safe;
}
 
/**
 * Mason (1989) stability function for turbulent mixing suppression.
 *
 * @param[in] Ri  gradient Richardson number [dimensionless]
 * @param[in] Ri_crit  critical Richardson number (default 0.25)
 * @return f(Ri) ∈ [0, 1] [dimensionless]
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
StabilityFunction(amrex::Real Ri, amrex::Real Ri_crit)
{
    amrex::Real stability_factor = amrex::Real(0.);
    if (Ri >= Ri_crit) {
        stability_factor = amrex::Real(0.);
    } else if (Ri > amrex::Real(0.) && Ri < Ri_crit) {
        stability_factor = std::sqrt(amrex::max(amrex::Real(0.), amrex::Real(1.) - Ri / Ri_crit));
    } else {
        stability_factor = amrex::Real(1.);
    }
    return stability_factor;
}
 
#endif // ERF_RICHARDSON_NUMBER_H_