ERF_TwoStreamSW.H

Go to the documentation of this file.

#ifndef ERF_TWO_STREAM_SW_H_
#define ERF_TWO_STREAM_SW_H_
 
#include <AMReX_GpuControl.H>
#include <AMReX_Math.H>
#include <AMReX_REAL.H>
#include <cmath>
 
/**
 * @file ERF_TwoStreamSW.H
 * @brief Shortwave (solar) radiation kernels: Beer-Lambert direct beam and
 * the two-stream layer solution for the diffuse field.
 *
 * The column model (ERF_TwoStreamColumn.H) treats the solar radiation as a
 * direct beam plus a diffuse field with upward and downward streams:
 *
 * - The direct beam is attenuated by Beer-Lambert,
 *     F_dir(m) = S0 * mu0 * exp(-tau_cum(m) / mu0),
 *   where tau_cum is the optical depth from the top of the atmosphere to
 *   interface m and mu0 the cosine of the solar zenith angle.
 *
 * - Each layer converts part of the direct beam into diffuse radiation
 *   (scattering) and reflects/transmits the diffuse streams. For a
 *   homogeneous layer with optical depth tau, single-scattering albedo omega
 *   and asymmetry factor g, compute_sw_layer_two_stream() returns the
 *   two-stream reflectance and transmittance for diffuse incidence (R_dif,
 *   T_dif) and the diffuse reflectance and transmittance generated by direct
 *   incidence (R_dir, T_dir), together with the direct transmittance
 *   T_noscat = exp(-tau/mu0). The gamma coefficients are those of the
 *   practical improved flux method (Zdunkowski et al. 1980), as used by
 *   RRTMGP; the layer solution follows Meador and Weaver (1980), Eqs. 14-15.
 *
 * - The layers are combined with the surface by the adding method: the
 *   surface reflects the fraction alpha of the direct and diffuse flux that
 *   reaches it, and the reflected radiation is scattered and absorbed again
 *   on its way up. Heating rates follow from the divergence of the net flux
 *   F_dir + F_diff_down - F_diff_up.
 *
 * With omega == 0 in every layer the diffuse field reduces to the reflected
 * direct beam only, and the absorbed surface flux reduces exactly to
 * (1 - alpha) times the Beer-Lambert direct beam.
 *
 * References:
 * -----------
 * - Meador, W. E., and W. R. Weaver, 1980: Two-stream approximations to
 *   radiative transfer in planetary atmospheres: A unified description
 *   of existing methods and a new improvement. J. Atmos. Sci., 37,
 *   630-643.
 * - Zdunkowski, W. G., R. M. Welch, and G. Korb, 1980: An investigation of
 *   the structure of typical two-stream methods for the calculation of
 *   solar fluxes and heating rates in clouds. Beitr. Phys. Atmos., 53,
 *   147-166.
 * - Toon, O. B., C. P. McKay, T. P. Ackerman, and K. Santhanam, 1989:
 *   Rapid calculation of radiative heating rates and photodissociation
 *   rates in inhomogeneous multiple scattering atmospheres. J. Geophys.
 *   Res., 94, 16287-16301.
 */
 
/**
 * @brief Compute Beer-Lambert direct-beam flux at a given optical depth.
 *
 *   F_dir = S0 * cos_zenith * exp(-tau_cumulative / cos_zenith)
 *
 * @param[in] tau_cumulative Optical depth from the top of the atmosphere [unitless].
 * @param[in] S0 Solar constant at the top of the atmosphere [W/m^2].
 * @param[in] cos_zenith Cosine of the solar zenith angle [unitless].
 *
 * @return Downwelling direct-beam flux [W/m^2], or 0 if cos_zenith <= 0 (night).
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real
compute_sw_direct_flux(amrex::Real tau_cumulative, amrex::Real S0,
                       amrex::Real cos_zenith)
{
  if (cos_zenith <= 0.0) {
    return 0.0;
  }
  return S0 * cos_zenith * std::exp(-tau_cumulative / cos_zenith);
}
 
/**
 * @brief Compute the shortwave heating rate of a layer from the net
 * (downward positive) flux at its top and bottom interfaces.
 *
 *   Q = (F_net_top - F_net_bot) / (dz * rho * cp)   [K/s]
 *
 * Energy converging into the layer warms it. Returns 0 for unphysical
 * inputs (dz, rho, cp <= 0) or a non-finite result.
 *
 * @param[in] flux_top Net downward flux at the top of the layer [W/m^2].
 * @param[in] flux_bot Net downward flux at the bottom of the layer [W/m^2].
 * @param[in] dz Vertical thickness of the layer [m]. Must be positive.
 *              For terrain-aware grids, pass dz = z_cc(k+1) - z_cc(k).
 * @param[in] rho Density [kg/m^3]. Must be positive.
 * @param[in] cp Specific heat at constant pressure [J/(kg·K)]. Must be positive.
 *
 * @return Heating rate [K/s]. Positive values indicate warming.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real
compute_sw_heating_rate(amrex::Real flux_top, amrex::Real flux_bot,
                        amrex::Real dz, amrex::Real rho, amrex::Real cp)
{
  if (dz <= 0.0 || rho <= 0.0 || cp <= 0.0) {
    return 0.0;
  }
  // Flux divergence: (W/m^2) / (m) = W/m^3
  amrex::Real flux_divergence = (flux_top - flux_bot) / dz;
  // Heating rate: (W/m^3) / (kg/m^3 * J/(kg*K)) = K/s
  amrex::Real heating = flux_divergence / (rho * cp);
  if (!amrex::Math::isfinite(heating)) {
    return 0.0;
  }
  return heating;
}
 
/**
 * @brief Two-stream reflectance and transmittance of one homogeneous layer.
 *
 * All quantities are fractions of the incident flux:
 * - R_dif, T_dif: reflectance and transmittance for diffuse incidence.
 * - R_dir, T_dir: diffuse flux reflected upward / transmitted downward per
 *   unit direct-beam flux incident at the top of the layer (the scattered
 *   part only; the surviving direct beam is T_noscat).
 * - T_noscat = exp(-tau / mu0): direct-beam transmittance.
 *
 * For a non-absorbing layer (omega = 1) R_dif + T_dif = 1 and
 * R_dir + T_dir + T_noscat = 1. For a non-scattering layer (omega = 0)
 * R_dif = R_dir = T_dir = 0 and T_dif = exp(-2 tau) (diffusivity factor 2).
 */
struct TwoStreamLayerSW
{
  amrex::Real R_dif;
  amrex::Real T_dif;
  amrex::Real R_dir;
  amrex::Real T_dir;
  amrex::Real T_noscat;
};
 
/**
 * @brief Compute the two-stream layer solution for shortwave radiation.
 *
 * Gamma coefficients (practical improved flux method, Zdunkowski et al.
 * 1980, as in RRTMGP):
 *   gamma1 = (8 - omega * (5 + 3 g)) / 4
 *   gamma2 = 3 omega (1 - g) / 4
 *   gamma3 = (2 - 3 g mu0) / 4
 *   gamma4 = 1 - gamma3
 *   k      = sqrt(gamma1^2 - gamma2^2)
 *
 * Diffuse incidence (Meador and Weaver 1980):
 *   R_dif = gamma2 (1 - e^{-2 k tau}) / D
 *   T_dif = 2 k e^{-k tau} / D,      D = k (1 + e^{-2 k tau}) + gamma1 (1 - e^{-2 k tau})
 *
 * Direct incidence (Meador and Weaver 1980, Eqs. 14-15), with
 * alpha1 = gamma1 gamma4 + gamma2 gamma3, alpha2 = gamma1 gamma3 + gamma2 gamma4:
 *   R_dir = omega / (D (1 - k^2 mu0^2)) *
 *           [ (1 - k mu0)(alpha2 + k gamma3) - (1 + k mu0)(alpha2 - k gamma3) e^{-2 k tau}
 *             - 2 (k gamma3 - alpha2 k mu0) e^{-k tau} T_noscat ]
 *   T_dir = - omega / (D (1 - k^2 mu0^2)) *
 *           [ (1 + k mu0)(alpha1 + k gamma4) T_noscat - (1 - k mu0)(alpha1 - k gamma4) e^{-2 k tau} T_noscat
 *             - 2 (k gamma4 + alpha1 k mu0) e^{-k tau} ]
 *
 * The removable singularity at k mu0 = 1 is avoided by nudging mu0 slightly,
 * and the results are clipped to the physical range (non-negative, and
 * R_dir + T_dir <= 1 - T_noscat).
 *
 * @param[in] tau Optical depth of the layer [unitless].
 * @param[in] omega Single-scattering albedo of the layer, in [0, 1].
 * @param[in] g Asymmetry factor of the layer, in [-1, 1].
 * @param[in] cos_zenith Cosine of the solar zenith angle, in (0, 1].
 * @return Layer reflectances and transmittances (see TwoStreamLayerSW).
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
TwoStreamLayerSW
compute_sw_layer_two_stream(amrex::Real tau, amrex::Real omega, amrex::Real g,
                            amrex::Real cos_zenith)
{
  TwoStreamLayerSW L{0.0, 1.0, 0.0, 0.0, 1.0};
 
  if (!(tau > 0.0) || !amrex::Math::isfinite(tau)) {
    return L;  // Empty layer: transparent to everything
  }
  if (!(cos_zenith > 0.0)) {
    // Night: no direct beam; the diffuse properties are still defined but
    // the caller does not use them.
    L.T_noscat = 0.0;
  }
 
  amrex::Real w0 = omega;
  if (!(w0 > 0.0) || !amrex::Math::isfinite(w0)) w0 = 0.0;
  if (w0 > 1.0) w0 = 1.0;
  amrex::Real asym = g;
  if (!amrex::Math::isfinite(asym)) asym = 0.0;
  if (asym > 1.0) asym = 1.0;
  if (asym < -1.0) asym = -1.0;
  amrex::Real mu0 = (cos_zenith > 0.0) ? cos_zenith : 1.0;
  if (mu0 > 1.0) mu0 = 1.0;
 
  const amrex::Real gamma1 = (8.0 - w0 * (5.0 + 3.0 * asym)) / 4.0;
  const amrex::Real gamma2 = 3.0 * w0 * (1.0 - asym) / 4.0;
  const amrex::Real gamma3 = (2.0 - 3.0 * asym * mu0) / 4.0;
  const amrex::Real gamma4 = 1.0 - gamma3;
 
  const amrex::Real alpha1 = gamma1 * gamma4 + gamma2 * gamma3;
  const amrex::Real alpha2 = gamma1 * gamma3 + gamma2 * gamma4;
 
  // k^2 = gamma1^2 - gamma2^2 >= 0; guard the conservative limit omega = 1.
  amrex::Real k_sq = (gamma1 - gamma2) * (gamma1 + gamma2);
  if (k_sq < 1.0e-12) k_sq = 1.0e-12;
  const amrex::Real k = std::sqrt(k_sq);
 
  const amrex::Real E  = std::exp(-k * tau);
  const amrex::Real E2 = E * E;
  const amrex::Real D  = k * (1.0 + E2) + gamma1 * (1.0 - E2);
  const amrex::Real RT = 1.0 / D;
 
  L.R_dif = RT * gamma2 * (1.0 - E2);
  L.T_dif = RT * 2.0 * k * E;
  if (L.R_dif < 0.0) L.R_dif = 0.0;
  if (L.T_dif < 0.0) L.T_dif = 0.0;
  if (L.R_dif + L.T_dif > 1.0) {
    const amrex::Real s = 1.0 / (L.R_dif + L.T_dif);
    L.R_dif *= s;
    L.T_dif *= s;
  }
 
  if (!(cos_zenith > 0.0)) {
    return L;
  }
 
  // Direct beam: nudge mu0 away from the removable singularity k mu0 = 1.
  amrex::Real k_mu = k * mu0;
  if (std::abs(1.0 - k_mu * k_mu) < 1.0e-4) {
    k_mu = (k_mu < 1.0) ? (1.0 - 1.0e-2) : (1.0 + 1.0e-2);
  }
  const amrex::Real T_noscat = std::exp(-tau / mu0);
  L.T_noscat = T_noscat;
 
  const amrex::Real RT2 = w0 * RT / (1.0 - k_mu * k_mu);
  const amrex::Real k_g3 = k * gamma3;
  const amrex::Real k_g4 = k * gamma4;
 
  amrex::Real R_dir = RT2 * ((1.0 - k_mu) * (alpha2 + k_g3)
                             - (1.0 + k_mu) * (alpha2 - k_g3) * E2
                             - 2.0 * (k_g3 - alpha2 * k_mu) * E * T_noscat);
  amrex::Real T_dir = -RT2 * ((1.0 + k_mu) * (alpha1 + k_g4) * T_noscat
                              - (1.0 - k_mu) * (alpha1 - k_g4) * E2 * T_noscat
                              - 2.0 * (k_g4 + alpha1 * k_mu) * E);
 
  if (!amrex::Math::isfinite(R_dir) || R_dir < 0.0) R_dir = 0.0;
  if (!amrex::Math::isfinite(T_dir) || T_dir < 0.0) T_dir = 0.0;
  const amrex::Real budget = 1.0 - T_noscat;  // scattered + absorbed fraction
  if (R_dir + T_dir > budget) {
    const amrex::Real s = (R_dir + T_dir > 0.0) ? budget / (R_dir + T_dir) : 0.0;
    R_dir *= s;
    T_dir *= s;
  }
  L.R_dir = R_dir;
  L.T_dir = T_dir;
  return L;
}
 
#endif // ERF_TWO_STREAM_SW_H_