ERF_TwoStreamSW.H File Reference

Shortwave (solar) radiation kernels: Beer-Lambert direct beam and the two-stream layer solution for the diffuse field. More...

#include <AMReX_GpuControl.H>
#include <AMReX_Math.H>
#include <AMReX_REAL.H>
#include <cmath>

Include dependency graph for ERF_TwoStreamSW.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	TwoStreamLayerSW
	Two-stream reflectance and transmittance of one homogeneous layer. More...

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	compute_sw_direct_flux (amrex::Real tau_cumulative, amrex::Real S0, amrex::Real cos_zenith)
	Compute Beer-Lambert direct-beam flux at a given optical depth. More...
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)
	Compute the shortwave heating rate of a layer from the net (downward positive) flux at its top and bottom interfaces. More...
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)
	Compute the two-stream layer solution for shortwave radiation. More...

Detailed Description

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.

Function Documentation

compute_sw_direct_flux()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_sw_direct_flux	(	amrex::Real	tau_cumulative,
		amrex::Real	S0,
		amrex::Real	cos_zenith
	)

Compute Beer-Lambert direct-beam flux at a given optical depth.

F_dir = S0 * cos_zenith * exp(-tau_cumulative / cos_zenith)

Parameters

[in]	tau_cumulative	Optical depth from the top of the atmosphere [unitless].
[in]	S0	Solar constant at the top of the atmosphere [W/m^2].
[in]	cos_zenith	Cosine of the solar zenith angle [unitless].

Returns

Downwelling direct-beam flux [W/m^2], or 0 if cos_zenith <= 0 (night).

{
  if (cos_zenith <= 0.0) {
    return 0.0;
  }
  return S0 * cos_zenith * std::exp(-tau_cumulative / cos_zenith);
}

Referenced by vertical_two_stream_sweep().

Here is the caller graph for this function:

compute_sw_heating_rate()

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
	)

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.

Parameters

[in]	flux_top	Net downward flux at the top of the layer [W/m^2].
[in]	flux_bot	Net downward flux at the bottom of the layer [W/m^2].
[in]	dz	Vertical thickness of the layer [m]. Must be positive. For terrain-aware grids, pass dz = z_cc(k+1) - z_cc(k).
[in]	rho	Density [kg/m^3]. Must be positive.
[in]	cp	Specific heat at constant pressure [J/(kg·K)]. Must be positive.

Returns

Heating rate [K/s]. Positive values indicate warming.

{
  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;
}

Referenced by vertical_two_stream_sweep().

Here is the caller graph for this function:

compute_sw_layer_two_stream()

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
	)

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).

Parameters

[in]	tau	Optical depth of the layer [unitless].
[in]	omega	Single-scattering albedo of the layer, in [0, 1].
[in]	g	Asymmetry factor of the layer, in [-1, 1].
[in]	cos_zenith	Cosine of the solar zenith angle, in (0, 1].

Returns

Layer reflectances and transmittances (see TwoStreamLayerSW).

{
  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;
}

Referenced by vertical_two_stream_sweep().

Here is the caller graph for this function: