ERF_PrognosticCloudFraction.H

Go to the documentation of this file.

#ifndef ERF_PROGNOSTIC_CLOUD_FRACTION_H_
#define ERF_PROGNOSTIC_CLOUD_FRACTION_H_
 
#include <AMReX_GpuControl.H>
#include <AMReX_Math.H>
#include <AMReX_FArrayBox.H>
#include <cmath>
 
/**
 * @file ERF_PrognosticCloudFraction.H
 * @brief Prognostic cloud fraction diagnosis from relative humidity and cloud water.
 *
 * Implements RH/qc-based cloud fraction diagnosis for use in TwoStream radiation.
 *
 * Cloud fraction diagnosis formula:
 *   cf_rh(k) = max(0, min(1, (rh(k) - rh_min) / (rh_max - rh_min)))  [linear ramp]
 *   cf_qc(k) = min(1, qc / qc_scale(k))  [qc contribution]
 *   cf_diag(k) = min(1, cf_rh(k) + cf_qc(k))  [combined with saturation]
 *
 * References:
 * -----------
 */
 
/**
 * @brief Compute relative humidity from water vapor mixing ratio.
 *
 * Approximate RH from qv and temperature using saturation mixing ratio.
 * The saturation vapour pressure over liquid water is the Magnus formula
 * with the WMO (2008) coefficients, applied to the temperature in Celsius.
 *
 * @param[in] qv Water vapor mixing ratio [kg/kg]
 * @param[in] T Temperature [K]
 * @param[in] P Pressure [Pa]
 * @return Relative humidity [0, 1]
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real compute_relative_humidity (amrex::Real qv, amrex::Real T, amrex::Real P)
{
    // Guard against invalid inputs
    if (T <= 0.0 || P <= 0.0 || qv < 0.0 || !amrex::Math::isfinite(qv) ||
        !amrex::Math::isfinite(T) || !amrex::Math::isfinite(P)) {
        return 0.0;
    }
 
    // Saturation vapour pressure over liquid water from the Magnus formula,
    // with the WMO (2008) coefficients for the temperature in Celsius:
    //   e_sat = e0 * exp(a * tc / (tc + b)),  tc = T - T0
    // and the saturation mixing ratio qsat = epsilon * e_sat / max(e_sat, P - e_sat).
    const amrex::Real e0 = 611.2;  // Pa (saturation vapor pressure at 0°C)
    const amrex::Real T0 = 273.15; // K
    const amrex::Real a = 17.62;
    const amrex::Real b = 243.12;  // °C
 
    // tc + b <= 0 only below about 30 K, where e_sat tends to 0. Otherwise
    // arg < a, so e_sat stays below e0 * exp(a) and cannot overflow.
    const amrex::Real tc = T - T0;
    const amrex::Real arg = (tc + b > amrex::Real(0.0)) ? a * tc / (tc + b) : amrex::Real(-100.0);
    const amrex::Real e_sat = e0 * std::exp(arg);
 
    // Saturation mixing ratio. The denominator is floored at e_sat, as in
    // erf_qsatw, so qsat stays positive and bounded by epsilon when e_sat
    // approaches or exceeds P (hot or very low-pressure air).
    const amrex::Real epsilon = 0.622;  // ratio of molecular weights
    amrex::Real qsat = epsilon * e_sat / amrex::max(e_sat, P - e_sat);
    if (qsat <= 0.0) {
        qsat = 1.0e-6;  // fallback when e_sat underflows (unphysically cold air)
    }
 
    // RH = qv / qsat, clamped to [0, 1]
    amrex::Real rh = qv / qsat;
    if (rh < 0.0) rh = 0.0;
    if (rh > 1.0) rh = 1.0;
 
    return rh;
}
 
/**
 * @brief Diagnose cloud fraction from relative humidity and cloud water.
 *
 * Computes diagnosed cloud fraction from RH and qc:
 *   - RH contribution: linear ramp from 0 at rh_min to 1 at rh_max
 *   - qc contribution: scaled by qc_scale coefficient
 *   - Combined: saturated blend (cf ≤ 1)
 *
 * This is a GPU-safe inline function for device-side kernels.
 *
 * @param[in] rh Relative humidity [0, 1]
 * @param[in] qc Cloud liquid water mixing ratio [kg/kg]
 * @param[in] rh_min Minimum RH threshold [0, 1]
 * @param[in] rh_max Maximum RH threshold [0, 1]
 * @param[in] qc_scale Scaling coefficient for qc contribution
 *
 * @return Diagnosed cloud fraction [0, 1]
 *
 * @note All inputs are assumed to be validated and finite.
 * @note Output is always in [0, 1], finite, and safe for further use.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real diagnose_cloud_fraction_from_rh_qc (
    amrex::Real rh, amrex::Real qc,
    amrex::Real rh_min, amrex::Real rh_max,
    amrex::Real qc_scale)
{
    // Guard against invalid inputs
    if (!amrex::Math::isfinite(rh) || !amrex::Math::isfinite(qc) ||
        !amrex::Math::isfinite(rh_min) || !amrex::Math::isfinite(rh_max) ||
        !amrex::Math::isfinite(qc_scale)) {
        return 0.0;  // Safe fallback on invalid input
    }
 
    // Clamp RH and qc to valid ranges
    if (rh < 0.0) rh = 0.0;
    if (rh > 1.0) rh = 1.0;
    if (qc < 0.0) qc = 0.0;
 
    amrex::Real cf = 0.0;
 
    // RH contribution: linear ramp from 0 at rh_min to 1 at rh_max
    if (rh_max > rh_min) {
        // Linear interpolation
        amrex::Real cf_rh = (rh - rh_min) / (rh_max - rh_min);
        if (cf_rh < 0.0) cf_rh = 0.0;
        if (cf_rh > 1.0) cf_rh = 1.0;
        cf += cf_rh;
    } else if (rh >= rh_min) {
        // rh_max == rh_min: step function
        cf += 1.0;
    }
 
    // qc contribution: qc_scale is the cloud water [kg/kg] at which the
    // liquid-water term alone saturates, so cf_qc = qc / qc_scale, capped at 1.
    if (qc_scale > 0.0) {
        amrex::Real cf_qc = qc / qc_scale;
        if (cf_qc > 1.0) cf_qc = 1.0;
        cf += cf_qc;
    }
 
    // Saturate combined cloud fraction to [0, 1]
    if (cf < 0.0) cf = 0.0;
    if (cf > 1.0) cf = 1.0;
 
    // Final sanity check
    if (!amrex::Math::isfinite(cf)) {
        cf = 0.0;
    }
 
    return cf;
}
 
#endif // ERF_PROGNOSTIC_CLOUD_FRACTION_H_