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;
 
    const bool is_covered_kp1 = c_cflag(i,j,k+1).isCovered();
    const bool is_covered_km1 = c_cflag(i,j,k-1).isCovered();
 
    if (saturated && qc_index >= 0) {
 
        const amrex::Real theta_l = GetThetal(i, j, k  , cell_data, moisture_indices);
 
        inv_theta = one / theta_l;
 
        if (is_covered_kp1) {
 
            if (is_covered_km1) {
                return zero;
            }
 
            const amrex::Real theta_l_km1 = GetThetal(i, j, k-1, cell_data, moisture_indices);
 
            if (c_cflag(i,j,k-2).isCovered()) {
                dtheta_dz = (theta_l - theta_l_km1) * dzInv;
            } else {
                const 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) {
 
            const amrex::Real theta_l_kp1 = GetThetal(i, j, k+1, cell_data, moisture_indices);
 
            if (c_cflag(i,j,k+2).isCovered()) {
                dtheta_dz = (theta_l_kp1 - theta_l) * dzInv;
            } else {
                const 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 {
 
            const amrex::Real theta_l_kp1 = GetThetal(i, j, k+1, cell_data, moisture_indices);
            const amrex::Real theta_l_km1 = GetThetal(i, j, k-1, cell_data, moisture_indices);
            dtheta_dz = (theta_l_kp1 - theta_l_km1) * myhalf * dzInv;
        }
 
    } else if (qv_index >= 0) {
 
        const amrex::Real theta_v = GetThetav(i, j, k, cell_data, moisture_indices);
 
        inv_theta = one / theta_v;
 
        if (is_covered_kp1) {
 
            if (is_covered_km1) {
                return zero;
            }
 
            const amrex::Real theta_v_km1 = GetThetav(i, j, k-1, cell_data, moisture_indices);
 
            if (c_cflag(i,j,k-2).isCovered()) {
                dtheta_dz = (theta_v - theta_v_km1) * dzInv;
            } else {
                const 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) {
 
            const amrex::Real theta_v_kp1 = GetThetav(i, j, k+1, cell_data, moisture_indices);
 
            if (c_cflag(i,j,k+2).isCovered()) {
                dtheta_dz = (theta_v_kp1 - theta_v) * dzInv;
            } else {
                const 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 {
 
            const amrex::Real theta_v_kp1 = GetThetav(i, j, k+1, cell_data, moisture_indices);
            const amrex::Real theta_v_km1 = GetThetav(i, j, k-1, cell_data, moisture_indices);
            dtheta_dz = (theta_v_kp1 - theta_v_km1) * myhalf * dzInv;
        }
 
    } else {
 
        const amrex::Real theta = cell_data(i, j, k, rhoTheta_comp) / rho;
 
        inv_theta = one / theta;
 
        if (is_covered_kp1) {
 
            if (is_covered_km1) {
                return zero;
            }
 
            const amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
 
            if (c_cflag(i,j,k-2).isCovered()) {
                dtheta_dz = (theta - theta_km1) * dzInv;
            } else {
                const 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) {
 
            const amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
 
            if (c_cflag(i,j,k+2).isCovered()) {
                dtheta_dz = (theta_kp1 - theta) * dzInv;
            } else {
                const 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 {
 
            const amrex::Real theta_kp1 = cell_data(i, j, k+1, rhoTheta_comp) / cell_data(i, j, k+1, rho_comp);
            const amrex::Real theta_km1 = cell_data(i, j, k-1, rhoTheta_comp) / cell_data(i, j, k-1, rho_comp);
            dtheta_dz = (theta_kp1 - theta_km1) * myhalf * 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
 * @param[out] shear_is_valid whether the required shear stencil is available
 * @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,
                      bool& shear_is_valid)
{
    // Volume-fraction weighted reconstruction of a velocity component at the
    // centre of cell (ii,jj,kk) from the two staggered faces that bracket it in
    // the (di,dj) direction. Returns false when neither face carries any volume,
    // in which case there is no cell-centred value to be had.
    auto cell_avg = [=] (const amrex::Array4<const amrex::Real>& vfrac,
                         const amrex::Array4<const amrex::Real>& vel,
                         int ii, int jj, int kk, int di, int dj,
                         amrex::Real& q) -> bool
    {
        const amrex::Real f0   = vfrac(ii,    jj,    kk);
        const amrex::Real f1   = vfrac(ii+di, jj+dj, kk);
        const amrex::Real fsum = f0 + f1;
        if (fsum <= zero) { return false; }
        q = ((f0 > zero ? f0 * vel(ii,    jj,    kk) : zero) +
             (f1 > zero ? f1 * vel(ii+di, jj+dj, kk) : zero)) / fsum;
        return true;
    };
 
    // Vertical derivative of one velocity component at cell (i,j,k), using the
    // best stencil the EB geometry leaves available: centred when both k-1 and
    // k+1 are usable, otherwise one-sided towards whichever neighbour is (second
    // order when k-2 / k+2 is usable as well). Returns false only when no stencil
    // at all can be formed, so a single unusable neighbour degrades the order of
    // accuracy rather than discarding the correction.
    auto vert_deriv = [=] (const amrex::Array4<const amrex::Real>& vfrac,
                           const amrex::Array4<const amrex::Real>& vel,
                           int di, int dj, amrex::Real& dq_dz) -> bool
    {
        amrex::Real c_km1 = zero, c_kp1 = zero;
        const bool have_km1 = !c_cflag(i,j,k-1).isCovered() &&
                              cell_avg(vfrac, vel, i, j, k-1, di, dj, c_km1);
        const bool have_kp1 = !c_cflag(i,j,k+1).isCovered() &&
                              cell_avg(vfrac, vel, i, j, k+1, di, dj, c_kp1);
 
        // The centred difference never reads the centre cell, so it stays valid
        // even when the centre cell's own faces carry no volume.
        if (have_km1 && have_kp1) {
            dq_dz = myhalf * (c_kp1 - c_km1) * dzInv;
            return true;
        }
 
        // Every remaining stencil is one-sided, and those do need the centre cell.
        amrex::Real c_c;
        if (!cell_avg(vfrac, vel, i, j, k, di, dj, c_c)) { return false; }
 
        if (have_km1) {
            amrex::Real c_km2;
            if (!c_cflag(i,j,k-2).isCovered() &&
                cell_avg(vfrac, vel, i, j, k-2, di, dj, c_km2)) {
                dq_dz = (three * c_c - four * c_km1 + c_km2) * myhalf * dzInv;
            } else {
                dq_dz = (c_c - c_km1) * dzInv;
            }
            return true;
        }
 
        if (have_kp1) {
            amrex::Real c_kp2;
            if (!c_cflag(i,j,k+2).isCovered() &&
                cell_avg(vfrac, vel, i, j, k+2, di, dj, c_kp2)) {
                dq_dz = (-c_kp2 + four * c_kp1 - three * c_c) * myhalf * dzInv;
            } else {
                dq_dz = (c_kp1 - c_c) * dzInv;
            }
            return true;
        }
 
        // Both vertical neighbours are unusable: no stencil exists at all.
        return false;
    };
 
    amrex::Real du_dz, dv_dz;
    if (!vert_deriv(u_vfrac, u, 1, 0, du_dz) ||
        !vert_deriv(v_vfrac, v, 0, 1, dv_dz)) {
        shear_is_valid = false;
        return zero;
    }
    shear_is_valid = true;
 
    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_