ERF_TwoStreamColumn.H File Reference

Per-column two-stream radiation kernels. More...

#include <algorithm>
#include <cmath>
#include <type_traits>
#include <AMReX_Array4.H>
#include <AMReX_Box.H>
#include <AMReX_Geometry.H>
#include <AMReX_GpuControl.H>
#include <AMReX_Math.H>
#include <AMReX_REAL.H>
#include <ERF_Constants.H>
#include <ERF_EOS.H>
#include <ERF_IndexDefines.H>
#include <ERF_RadStruct.H>
#include <ERF_TwoStreamSW.H>
#include <ERF_TwoStreamLW.H>
#include <ERF_PrognosticCloudFraction.H>
#include <ERF_AerosolOpticalDepth.H>
#include <ERF_OrbCosZenith.H>
#include <ERF_SurfaceTemperature.H>
#include <ERF_TerrainMetrics.H>

Include dependency graph for ERF_TwoStreamColumn.H:

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

Go to the source code of this file.

Classes
struct	TwoStreamParams
	Trivially copyable subset of RadChoice used by the column kernels. More...

Namespaces
	TwoStreamScratch
	Per-column scratch the sweep keeps in a caller-provided Array4.

Functions
TwoStreamParams	make_two_stream_params (const RadChoice &rc, const amrex::Real rdOcp)
	Copy the column-kernel parameters out of a RadChoice. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	get_temperature_from_rhotheta (amrex::Real rho_theta, amrex::Real rho, amrex::Real qv=0.0)
	GPU-safe helper to compute absolute temperature from (rho, rho*theta). More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool	valid_surface_temperature (amrex::Real value)
	Check whether a heterogeneous absolute surface temperature is valid. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	resolve_surface_albedo_sw (int i, int j, const amrex::Array4< const amrex::Real > *hetero_alb_sw, const TwoStreamParams &rad_choice, bool has_hetero_alb)
	Resolve per-column shortwave surface albedo from hetero field or fallback. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	resolve_surface_emissivity_lw (int i, int j, const amrex::Array4< const amrex::Real > *hetero_emiss_lw, const TwoStreamParams &rad_choice, bool has_hetero_emiss)
	Resolve per-column longwave surface emissivity from hetero field or fallback. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	resolve_surface_temp_k (int i, int j, const amrex::Array4< const amrex::Real > *lsm_t_sfc, bool has_lsm_t_sfc, const amrex::Array4< const amrex::Real > *seb_t_sfc, bool has_seb_t_sfc, const amrex::Array4< const amrex::Real > *surface_layer_theta, bool has_surface_layer, const TwoStreamParams &rad_choice, bool &from_surface_layer)
	Resolve the per-column surface-temperature boundary condition. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool	is_cloud_level (amrex::Real z_center, const TwoStreamParams &rad_choice)
	GPU-safe helper to determine whether a layer falls within the cloud band [cloud_base_height_m, cloud_top_height_m]. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	get_qv_from_state (int i, int j, int k, const amrex::Array4< const amrex::Real > &state_arr)
	Water-vapor mixing ratio qv = RhoQv / Rho at (i,j,k). More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	get_inverse_exner (int i, int j, int k, const amrex::Array4< const amrex::Real > &state_arr)
	Inverse Exner function 1/pi at (i,j,k) from the state, with pi = (p / p_0)^(R_d/c_p) evaluated by getExnergivenRTh(). Used to convert a temperature tendency dT/dt into the potential-temperature tendency dtheta/dt = (dT/dt) / pi that the RhoTheta source term expects. Returns 1 for unphysical input. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	get_qc_from_state (int i, int j, int k, const amrex::Array4< const amrex::Real > &state_arr)
	Cloud-water mixing ratio qc = RhoQc / Rho at (i,j,k). More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	diagnose_tau_dynamic (amrex::Real tau_base, amrex::Real qv, amrex::Real qc, amrex::Real coeff_qv, amrex::Real coeff_qc)
	Moisture-dependent (dynamic) per-layer optical depth. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	tau_layer_value (amrex::Real z_center, amrex::Real tau_base, const TwoStreamParams &rad_choice, bool apply_cloud)
	GPU-safe helper to compute the per-layer optical depth at level k, given the base (clear-sky) optical depth and cloud-layer parameters. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	select_scattering_props (amrex::Real z_center, const TwoStreamParams &rad_choice, bool apply_cloud, amrex::Real &omega, amrex::Real &g)
	GPU-safe helper to select the single-scattering albedo and asymmetry factor to use for level k's diffuse SW calculation. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	diagnose_cloud_fraction_prognostic (int i, int j, int k, const amrex::Array4< const amrex::Real > &state_arr, const TwoStreamParams &rad_choice)
	GPU-safe helper to diagnose prognostic cloud fraction from per-level relative humidity and cloud liquid water. More...
amrex::Box	two_stream_scratch_box (const amrex::Box &bx)
	Box of the scratch FArrayBox a sweep over box bx needs. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	diagnose_layer_tau (int i, int j, int k, amrex::Real dz_layer, amrex::Real z_center, const amrex::Array4< const amrex::Real > &state_arr, amrex::Real tau_base, bool is_sw, bool cloudy, const TwoStreamParams &rad_choice)
	GPU-safe helper to assemble the total per-layer optical depth at level k for either the SW or the LW band. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	diagnose_layer_optics (int i, int j, int k, amrex::Real dz_layer, amrex::Real z_center, const amrex::Array4< const amrex::Real > &state_arr, amrex::Real tau_base, bool is_sw, bool cloudy, const TwoStreamParams &rad_choice, amrex::Real &tau, amrex::Real &omega, amrex::Real &g)
	Optical depth, single-scattering albedo and asymmetry factor of layer k for one band, for either optical-depth model. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	vertical_two_stream_sweep (int i, int j, const amrex::Box &bx, amrex::Real dz_uniform, const amrex::Array4< const amrex::Real > &state_arr, const TwoStreamParams &rad_choice, bool cloudy, const amrex::Array4< amrex::Real > &qheating_arr, amrex::Real &max_heating_rate, amrex::Real &sw_surface_flux, amrex::Real &sw_up_toa, amrex::Real &lw_net_surface, amrex::Real &lw_up_toa, amrex::Real &sw_down_toa, const amrex::Array4< const amrex::Real > &z_phys_nd, const amrex::Array4< amrex::Real > &scratch, bool has_hetero_alb_sw=false, const amrex::Array4< const amrex::Real > *hetero_alb_sw=nullptr, bool has_hetero_emiss_lw=false, const amrex::Array4< const amrex::Real > *hetero_emiss_lw=nullptr, bool has_lsm_t_sfc=false, const amrex::Array4< const amrex::Real > *lsm_t_sfc=nullptr, bool has_seb_t_sfc=false, const amrex::Array4< const amrex::Real > *seb_t_sfc=nullptr, bool has_surface_layer=false, const amrex::Array4< const amrex::Real > *surface_layer_theta=nullptr, bool has_latlon=false, const amrex::Array4< const amrex::Real > *lat_arr=nullptr, const amrex::Array4< const amrex::Real > *lon_arr=nullptr, const amrex::Array4< amrex::Real > *rad_flux_out=nullptr)
	GPU-safe per-column vertical integration kernel for two-stream radiation, computing either the clear-sky or cloudy-column fluxes and per-level heating rates, depending on the cloudy flag. More...

Variables
constexpr int	TwoStreamScratch::DZ = 0
constexpr int	TwoStreamScratch::Z = 1
constexpr int	TwoStreamScratch::TAU_SW = 2
constexpr int	TwoStreamScratch::F_DIR = 3
constexpr int	TwoStreamScratch::A = 4
constexpr int	TwoStreamScratch::S = 5
constexpr int	TwoStreamScratch::LW_UP = 6
constexpr int	TwoStreamScratch::LW_DN = 7
constexpr int	TwoStreamScratch::TAU_LW = 8
constexpr int	TwoStreamScratch::NCOMP = 9

Detailed Description

Per-column two-stream radiation kernels.

Everything needed to evaluate one (i,j) column of the two-stream model: temperature from the equation of state, per-column surface properties, per-layer optical depth assembly, and vertical_two_stream_sweep(), which performs the SW and LW sweeps and writes per-level heating rates.

Vertical orientation follows ERF: k = kmin is the surface layer and k = kmax the top layer. SW sweeps downward from kmax to kmin; LW sweeps downward (TOA -> surface) and then upward (surface -> TOA) on layer interfaces. Layer temperature is obtained from rho*theta through the Exner function.

The functions are header-only and GPU-safe so that the column physics can be exercised directly by unit tests (Tests/Unit/Radiation) as well as by the level driver in ERF_TwoStreamRadiation.cpp (TwoStreamRadiation::advance).

Function Documentation

diagnose_cloud_fraction_prognostic()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real diagnose_cloud_fraction_prognostic	(	int	i,
		int	j,
		int	k,
		const amrex::Array4< const amrex::Real > &	state_arr,
		const TwoStreamParams &	rad_choice
	)

GPU-safe helper to diagnose prognostic cloud fraction from per-level relative humidity and cloud liquid water.

Computes cloud fraction from RH and qc using: cf_rh(k) = linear ramp from 0 at rh_min to 1 at rh_max cf_qc(k) = qc_scale * qc(k) cf(k) = min(1, cf_rh + cf_qc) [saturated blend]

Temperature comes from the equation of state (Exner function) and pressure from getPgivenRTh(), both evaluated with the local qv. Returns 0 when cloud_fraction_prog_enable is false.

Parameters

[in]	i,j,k	Grid indices
[in]	state_arr	State array proxy (contains Rho, RhoTheta, qv, qc)
[in]	rad_choice	Radiation parameters (prognostic cloud fraction settings)

Returns

Diagnosed cloud fraction [0, 1] if enabled; 0 if disabled

{
    if (!rad_choice.cloud_fraction_prog_enable) {
        return 0.0;
    }
 
    const amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
    const amrex::Real qc = get_qc_from_state(i, j, k, state_arr);
 
    amrex::Real rho       = state_arr(i, j, k, Rho_comp);
    amrex::Real rho_theta = state_arr(i, j, k, RhoTheta_comp);
    // Unreachable from TwoStreamRadiation::advance, which refuses a state with
    // a non-finite or non-positive density or rho*theta before the sweep;
    // kept so a direct kernel caller (unit tests) cannot divide by zero.
    if (rho <= 0.0 || !amrex::Math::isfinite(rho)) rho = 1.0;
    if (rho_theta <= 0.0 || !amrex::Math::isfinite(rho_theta)) rho_theta = 288.15;
 
    // Temperature (via the Exner function) and pressure from the equation of state
    amrex::Real T = get_temperature_from_rhotheta(rho_theta, rho, qv);
    amrex::Real P = getPgivenRTh(rho_theta, qv);
    if (P <= 0.0 || !amrex::Math::isfinite(P)) P = p_0;
 
    const amrex::Real rh = compute_relative_humidity(qv, T, P);
 
    return diagnose_cloud_fraction_from_rh_qc(
        rh, qc,
        rad_choice.cloud_fraction_rh_min,
        rad_choice.cloud_fraction_rh_max,
        rad_choice.cloud_fraction_qc_scale);
}

Referenced by diagnose_layer_optics(), and diagnose_layer_tau().

Here is the call graph for this function:

Here is the caller graph for this function:

diagnose_layer_optics()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void diagnose_layer_optics	(	int	i,
		int	j,
		int	k,
		amrex::Real	dz_layer,
		amrex::Real	z_center,
		const amrex::Array4< const amrex::Real > &	state_arr,
		amrex::Real	tau_base,
		bool	is_sw,
		bool	cloudy,
		const TwoStreamParams &	rad_choice,
		amrex::Real &	tau,
		amrex::Real &	omega,
		amrex::Real &	g
	)

Optical depth, single-scattering albedo and asymmetry factor of layer k for one band, for either optical-depth model.

Per-layer model (default): the optical depth comes from diagnose_layer_tau() and the scattering properties from the clear-sky or cloud-band inputs through select_scattering_props(), exactly as before.

Mass model (tau_model = mass), shortwave: each constituent contributes an extinction optical depth from the layer mass path, dry absorption rho dz sw_kabs_dry (omega = 0) Rayleigh rho dz sw_kscat_dry (omega = 1, g = 0) water vapor rho dz sw_kabs_vapor qv (omega = 0) cloud water rho dz sw_kext_cloud qc (sw_cloud_omega, sw_cloud_g) plus the prescribed cloud-band enhancement (with the cloud-band scattering inputs) and the absorbing moisture-coefficient and aerosol additions. The layer properties are the extinction-weighted mixtures tau = sum tau_i, omega = sum omega_i tau_i / tau, g = sum g_i omega_i tau_i / sum omega_i tau_i. Longwave uses diagnose_layer_tau(), whose clear-sky base is the mass path rho dz (lw_kabs_dry + lw_kabs_vapor qv + lw_kabs_cloud qc) in this model.

Parameters

[in]	i,j,k	Grid indices (kmin is the surface layer).
[in]	dz_layer	Thickness of this layer [m].
[in]	z_center	Height of the layer center above the surface [m].
[in]	state_arr	State array (read-only).
[in]	tau_base	Per-layer-model clear-sky optical depth for this band.
[in]	is_sw	true for the shortwave band, false for longwave.
[in]	cloudy	true for the cloudy-column evaluation.
[in]	rad_choice	Column-kernel parameters.
[out]	tau	Optical depth of the layer [unitless].
[out]	omega	Single-scattering albedo of the layer (0 for LW).
[out]	g	Asymmetry factor of the layer (0 for LW).

{
    omega = 0.0;
    g = 0.0;
    if (!is_sw || rad_choice.tau_model != TauModel::Mass) {
        tau = diagnose_layer_tau(i, j, k, dz_layer, z_center, state_arr, tau_base, is_sw, cloudy, rad_choice);
        if (is_sw) {
            select_scattering_props(z_center, rad_choice, cloudy, omega, g);
        }
        return;
    }
 
    // Shortwave mass model: constituent extinction optical depths.
    amrex::Real rho = state_arr(i, j, k, Rho_comp);
    if (!(rho > 0.0) || !amrex::Math::isfinite(rho)) rho = 0.0;
    const amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
    const amrex::Real qc = get_qc_from_state(i, j, k, state_arr);
    const amrex::Real path = rho * dz_layer;
 
    const amrex::Real tau_dry_abs = path * rad_choice.sw_kabs_dry;
    const amrex::Real tau_rayleigh = path * rad_choice.sw_kscat_dry;
    const amrex::Real tau_vapor = path * rad_choice.sw_kabs_vapor * qv;
    const amrex::Real tau_cloud = path * rad_choice.sw_kext_cloud * qc;
 
    amrex::Real ext = tau_dry_abs + tau_rayleigh + tau_vapor + tau_cloud;
    amrex::Real sca = tau_rayleigh + rad_choice.sw_cloud_omega * tau_cloud;
    amrex::Real gsca = rad_choice.sw_cloud_g * rad_choice.sw_cloud_omega * tau_cloud;
 
    // Prescribed cloud band (optionally scaled by the prognostic cloud
    // fraction) with the cloud-band scattering inputs.
    amrex::Real band = tau_layer_value(z_center, 0.0, rad_choice, cloudy);
    if (rad_choice.cloud_fraction_prog_enable && cloudy &&
        rad_choice.tau_profile_type == TauProfileType::CloudLayer &&
        is_cloud_level(z_center, rad_choice)) {
        band = diagnose_cloud_fraction_prognostic(i, j, k, state_arr, rad_choice) * rad_choice.cloud_tau_per_layer;
    }
    if (band > 0.0) {
        ext += band;
        sca += rad_choice.cloud_single_scattering_albedo * band;
        gsca += rad_choice.cloud_asymmetry_factor * rad_choice.cloud_single_scattering_albedo * band;
    }
 
    // Absorbing additions (moisture coefficients, aerosol): everything that
    // diagnose_layer_tau() adds on top of a zero base, minus the band.
    const amrex::Real additions = diagnose_layer_tau(i, j, k, dz_layer, z_center, state_arr,
                                                     0.0, /*is_sw=*/true, cloudy, rad_choice) - band;
    if (additions > 0.0) ext += additions;
 
    if (!amrex::Math::isfinite(ext) || ext <= 0.0) {
        tau = 0.0;
        return;
    }
    tau = ext;
    omega = sca / ext;
    if (omega > 1.0) omega = 1.0;
    if (omega < 0.0) omega = 0.0;
    g = (sca > 0.0) ? gsca / sca : 0.0;
    if (g > 1.0) g = 1.0;
    if (g < -1.0) g = -1.0;
}

Referenced by vertical_two_stream_sweep().

Here is the call graph for this function:

Here is the caller graph for this function:

diagnose_layer_tau()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real diagnose_layer_tau	(	int	i,
		int	j,
		int	k,
		amrex::Real	dz_layer,
		amrex::Real	z_center,
		const amrex::Array4< const amrex::Real > &	state_arr,
		amrex::Real	tau_base,
		bool	is_sw,
		bool	cloudy,
		const TwoStreamParams &	rad_choice
	)

GPU-safe helper to assemble the total per-layer optical depth at level k for either the SW or the LW band.

The contributions are applied in the same order for both bands:

1. clear-sky base value, plus the cloud-layer enhancement when the column is evaluated as "cloudy" and the level lies inside the cloud band;
2. optional moisture-dependent (dynamic) term from qv and qc;
3. optional prognostic cloud fraction, which replaces the cloud-band enhancement by cf(k) * cloud_tau_per_layer;
4. optional prescribed aerosol term (constant, exponential or table).

Parameters

[in]	i,j,k	Grid indices (k increases upward; kmin is the surface layer).
[in]	dz_layer	Thickness of this layer [m].
[in]	z_center	Height of the layer center above the surface [m] (cloud-band detection and the aerosol profile).
[in]	state_arr	State array (read-only).
[in]	tau_base	Clear-sky optical depth per layer for this band (for LW, replaced by the mass-path value when lw_mass_absorption_enable or tau_model = mass is set).
[in]	is_sw	true for the shortwave band, false for longwave.
[in]	cloudy	true for the cloudy-column evaluation.
[in]	rad_choice	Radiation parameters.

Returns

Optical depth of layer k [unitless].

{
    // Longwave option: gray optical depth from the layer mass path,
    //   tau = rho dz (k_dry + k_vapor qv + k_cloud qc),
    // which makes the column optical depth independent of the vertical
    // resolution and gives water vapor and cloud water a real greenhouse
    // effect. Replaces the fixed tau_lw_per_layer as the clear-sky base.
    amrex::Real base = tau_base;
    if (!is_sw && (rad_choice.lw_mass_absorption_enable || rad_choice.tau_model == TauModel::Mass)) {
        amrex::Real rho = state_arr(i, j, k, Rho_comp);
        if (!(rho > 0.0) || !amrex::Math::isfinite(rho)) rho = 0.0;
        const amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
        const amrex::Real qc = get_qc_from_state(i, j, k, state_arr);
        base = rho * dz_layer * (rad_choice.lw_kabs_dry
                                 + rad_choice.lw_kabs_vapor * qv
                                 + rad_choice.lw_kabs_cloud * qc);
        if (!amrex::Math::isfinite(base) || base < 0.0) base = 0.0;
    }
 
    // Base (clear-sky) value, plus cloud-band enhancement for cloudy columns.
    amrex::Real tau = tau_layer_value(z_center, base, rad_choice, cloudy);
 
    // Dynamic (moisture-dependent) contribution.
    const bool dynamic_enabled = is_sw ? rad_choice.tau_sw_dynamic_enable
                                       : rad_choice.tau_lw_dynamic_enable;
    if (dynamic_enabled) {
        const amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
        const amrex::Real qc = get_qc_from_state(i, j, k, state_arr);
        tau = is_sw ? diagnose_tau_dynamic(tau, qv, qc, rad_choice.tau_sw_coeff_qv, rad_choice.tau_sw_coeff_qc)
                    : diagnose_tau_dynamic(tau, qv, qc, rad_choice.tau_lw_coeff_qv, rad_choice.tau_lw_coeff_qc);
    }
 
    // Prognostic cloud fraction: scale the cloud-band enhancement by cf(k).
    if (rad_choice.cloud_fraction_prog_enable && cloudy &&
        rad_choice.tau_profile_type == TauProfileType::CloudLayer &&
        is_cloud_level(z_center, rad_choice)) {
        // tau already holds the full cloud-band enhancement (this is a cloudy
        // column inside the band) plus the dynamic term; scale only the
        // enhancement by cf(k) and leave the dynamic term in place.
        amrex::Real cf_prog = diagnose_cloud_fraction_prognostic(i, j, k, state_arr, rad_choice);
        tau += (cf_prog - 1.0) * rad_choice.cloud_tau_per_layer;
    }
 
    // Prescribed bulk aerosol contribution (added on top of everything above).
    if (rad_choice.aerosol_enable) {
        amrex::Real tau_aerosol = 0.0;
        if (rad_choice.aerosol_profile_type == AerosolProfileType::Constant) {
            tau_aerosol = diagnose_tau_aerosol_constant(rad_choice.aerosol_tau_per_layer);
        } else if (rad_choice.aerosol_profile_type == AerosolProfileType::Exponential) {
            tau_aerosol = diagnose_tau_aerosol_exponential(z_center, dz_layer,
                                                           rad_choice.aerosol_tau_surface,
                                                           rad_choice.aerosol_scale_height_m);
        } else if (rad_choice.aerosol_profile_type == AerosolProfileType::Table) {
            tau_aerosol = diagnose_tau_aerosol_table(k);
        }
        tau += tau_aerosol;
    }
 
    return tau;
}

Referenced by diagnose_layer_optics(), and vertical_two_stream_sweep().

Here is the call graph for this function:

Here is the caller graph for this function:

diagnose_tau_dynamic()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real diagnose_tau_dynamic	(	amrex::Real	tau_base,
		amrex::Real	qv,
		amrex::Real	qc,
		amrex::Real	coeff_qv,
		amrex::Real	coeff_qc
	)

Moisture-dependent (dynamic) per-layer optical depth.

tau = tau_base + coeff_qv * qv + coeff_qc * qc

clamped to [0, 100]. With both coefficients zero this returns tau_base unchanged, so the static configuration is reproduced exactly. Used for both the SW and the LW band with the respective coefficient pair.

Parameters

[in]	tau_base	Optical depth before the moisture contribution [unitless].
[in]	qv	Water-vapor mixing ratio [kg/kg].
[in]	qc	Cloud-water mixing ratio [kg/kg].
[in]	coeff_qv	Optical depth per unit qv [unitless per kg/kg].
[in]	coeff_qc	Optical depth per unit qc [unitless per kg/kg].

Returns

Optical depth for this layer [unitless].

{
    if (!(qv > 0.0) || !amrex::Math::isfinite(qv)) qv = 0.0;
    if (!(qc > 0.0) || !amrex::Math::isfinite(qc)) qc = 0.0;
    amrex::Real tau = tau_base + coeff_qv * qv + coeff_qc * qc;
    if (!amrex::Math::isfinite(tau)) return tau_base;
    if (tau < 0.0) tau = 0.0;
    if (tau > 100.0) tau = 100.0;
    return tau;
}

Referenced by diagnose_layer_tau().

Here is the caller graph for this function:

get_inverse_exner()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_inverse_exner	(	int	i,
		int	j,
		int	k,
		const amrex::Array4< const amrex::Real > &	state_arr
	)

Inverse Exner function 1/pi at (i,j,k) from the state, with pi = (p / p_0)^(R_d/c_p) evaluated by getExnergivenRTh(). Used to convert a temperature tendency dT/dt into the potential-temperature tendency dtheta/dt = (dT/dt) / pi that the RhoTheta source term expects. Returns 1 for unphysical input.

{
    amrex::Real rho_theta = state_arr(i, j, k, RhoTheta_comp);
    if (!(rho_theta > 0.0) || !amrex::Math::isfinite(rho_theta)) return 1.0;
    const amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
    const amrex::Real pi = getExnergivenRTh(rho_theta, RdoCp, qv);
    return (amrex::Math::isfinite(pi) && pi > 0.0) ? 1.0 / pi : 1.0;
}

Referenced by vertical_two_stream_sweep().

Here is the call graph for this function:

Here is the caller graph for this function:

get_qc_from_state()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_qc_from_state	(	int	i,
		int	j,
		int	k,
		const amrex::Array4< const amrex::Real > &	state_arr
	)

Cloud-water mixing ratio qc = RhoQc / Rho at (i,j,k).

Returns 0 when the state carries no cloud-water component, or when the stored values are non-finite or negative.

{
    if (state_arr.nComp() <= RhoQ2_comp) return 0.0;
    amrex::Real rho    = state_arr(i, j, k, Rho_comp);
    amrex::Real rho_qc = state_arr(i, j, k, RhoQ2_comp);
    if (!(rho > 0.0) || !amrex::Math::isfinite(rho) || !amrex::Math::isfinite(rho_qc)) return 0.0;
    amrex::Real qc = rho_qc / rho;
    return (amrex::Math::isfinite(qc) && qc > 0.0) ? qc : 0.0;
}

Referenced by diagnose_cloud_fraction_prognostic(), diagnose_layer_optics(), and diagnose_layer_tau().

Here is the caller graph for this function:

get_qv_from_state()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_qv_from_state	(	int	i,
		int	j,
		int	k,
		const amrex::Array4< const amrex::Real > &	state_arr
	)

Water-vapor mixing ratio qv = RhoQv / Rho at (i,j,k).

Returns 0 when the state carries no moisture components, or when the stored values are non-finite or negative.

{
    if (state_arr.nComp() <= RhoQ1_comp) return 0.0;
    amrex::Real rho    = state_arr(i, j, k, Rho_comp);
    amrex::Real rho_qv = state_arr(i, j, k, RhoQ1_comp);
    if (!(rho > 0.0) || !amrex::Math::isfinite(rho) || !amrex::Math::isfinite(rho_qv)) return 0.0;
    amrex::Real qv = rho_qv / rho;
    return (amrex::Math::isfinite(qv) && qv > 0.0) ? qv : 0.0;
}

Referenced by diagnose_cloud_fraction_prognostic(), diagnose_layer_optics(), diagnose_layer_tau(), get_inverse_exner(), and vertical_two_stream_sweep().

Here is the caller graph for this function:

get_temperature_from_rhotheta()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_temperature_from_rhotheta	(	amrex::Real	rho_theta,
		amrex::Real	rho,
		amrex::Real	qv = 0.0
	)

GPU-safe helper to compute absolute temperature from (rho, rho*theta).

ERF stores dry density and dry potential temperature. The pressure follows from the equation of state, p = p_0 * (R_d * rho * theta_m / p_0)^gamma, theta_m = theta * (1 + R_v/R_d * qv), and the absolute temperature is recovered through the Exner function, T = theta * (p / p_0)^(R_d / c_p), which getTgivenRandRTh() evaluates as p / (R_d * rho * (1 + R_v/R_d * qv)).

Parameters

[in]	rho_theta	RhoTheta component [K·kg/m^3]
[in]	rho	Density [kg/m^3]
[in]	qv	Water-vapor mixing ratio [kg/kg] (0 for dry air)

Returns

Temperature [K], clamped to [100, 400]; 288.15 for unphysical input.

{
    if (!(rho > 0.0) || !amrex::Math::isfinite(rho) ||
        !(rho_theta > 0.0) || !amrex::Math::isfinite(rho_theta)) {
        return 288.15;  // Defensive: fallback to standard T
    }
    if (!(qv >= 0.0) || !amrex::Math::isfinite(qv)) {
        qv = 0.0;
    }
 
    amrex::Real T = getTgivenRandRTh(rho, rho_theta, qv);
    if (!amrex::Math::isfinite(T) || T <= 0.0) {
        return 288.15;  // Defensive: fallback
    }
 
    // Defensive clipping to a sensible terrestrial range
    T = std::max(T, amrex::Real(100.0));
    T = std::min(T, amrex::Real(400.0));
    return T;
}

Referenced by diagnose_cloud_fraction_prognostic(), and vertical_two_stream_sweep().

Here is the call graph for this function:

Here is the caller graph for this function:

is_cloud_level()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_cloud_level	(	amrex::Real	z_center,
		const TwoStreamParams &	rad_choice
	)

GPU-safe helper to determine whether a layer falls within the cloud band [cloud_base_height_m, cloud_top_height_m].

Parameters

[in]	z_center	Height of the layer center above the surface [m].
[in]	rad_choice	Radiation parameters.

Returns

true if this layer is inside the configured cloud band.

{
    return (z_center >= rad_choice.cloud_base_height_m &&
            z_center <= rad_choice.cloud_top_height_m);
}

Referenced by diagnose_layer_optics(), diagnose_layer_tau(), select_scattering_props(), and tau_layer_value().

Here is the caller graph for this function:

make_two_stream_params()

TwoStreamParams make_two_stream_params ( const RadChoice &  rc, const amrex::Real  rdOcp  )

inline

Copy the column-kernel parameters out of a RadChoice.

{
    TwoStreamParams p;
    p.sw_enabled = rc.sw_enabled;
    p.lw_enabled = rc.lw_enabled;
    p.tau_per_layer = rc.tau_per_layer;
    p.tau_lw_per_layer = rc.tau_lw_per_layer;
    p.tau_profile_type = rc.tau_profile_type;
    p.cloud_base_height_m = rc.cloud_base_height_m;
    p.cloud_top_height_m = rc.cloud_top_height_m;
    p.cloud_tau_per_layer = rc.cloud_tau_per_layer;
    p.single_scattering_albedo = rc.single_scattering_albedo;
    p.asymmetry_factor = rc.asymmetry_factor;
    p.cloud_single_scattering_albedo = rc.cloud_single_scattering_albedo;
    p.cloud_asymmetry_factor = rc.cloud_asymmetry_factor;
    p.surface_albedo_sw = rc.surface_albedo_sw;
    p.surface_albedo_sw_diffuse = rc.surface_albedo_sw_diffuse;
    p.surface_emissivity_lw = rc.surface_emissivity_lw;
    p.t_sfc_default = rc.rad_t_sfc;
    p.rdOcp = rdOcp;
    p.tau_sw_dynamic_enable = rc.tau_sw_dynamic_enable;
    p.tau_lw_dynamic_enable = rc.tau_lw_dynamic_enable;
    p.tau_sw_coeff_qv = rc.tau_sw_coeff_qv;
    p.tau_sw_coeff_qc = rc.tau_sw_coeff_qc;
    p.tau_lw_coeff_qv = rc.tau_lw_coeff_qv;
    p.tau_lw_coeff_qc = rc.tau_lw_coeff_qc;
    p.cloud_fraction_prog_enable = rc.cloud_fraction_prog_enable;
    p.cloud_fraction_rh_min = rc.cloud_fraction_rh_min;
    p.cloud_fraction_rh_max = rc.cloud_fraction_rh_max;
    p.cloud_fraction_qc_scale = rc.cloud_fraction_qc_scale;
    p.aerosol_enable = rc.aerosol_enable;
    p.aerosol_profile_type = rc.aerosol_profile_type;
    p.aerosol_tau_per_layer = rc.aerosol_tau_per_layer;
    p.aerosol_scale_height_m = rc.aerosol_scale_height_m;
    p.aerosol_tau_surface = rc.aerosol_tau_surface;
    p.tau_model = rc.tau_model;
    p.sw_kabs_dry = rc.sw_kabs_dry;
    p.sw_kscat_dry = rc.sw_kscat_dry;
    p.sw_kabs_vapor = rc.sw_kabs_vapor;
    p.sw_kext_cloud = rc.sw_kext_cloud;
    p.sw_cloud_omega = rc.sw_cloud_omega;
    p.sw_cloud_g = rc.sw_cloud_g;
    p.lw_mass_absorption_enable = rc.lw_mass_absorption_enable;
    p.lw_kabs_dry = rc.lw_kabs_dry;
    p.lw_kabs_vapor = rc.lw_kabs_vapor;
    p.lw_kabs_cloud = rc.lw_kabs_cloud;
    // Sun: a fixed cosine when given; the date-dependent quantities (calday,
    // declination and the irradiance of the date) are filled by the driver.
    p.solar_dynamic = !(rc.fixed_solar_zenith_angle > 0.0);
    p.cos_zenith_fixed = rc.fixed_solar_zenith_angle;
    p.S0 = (rc.fixed_total_solar_irradiance >= 0.0) ? rc.fixed_total_solar_irradiance : 1360.9;
    p.lat_cons_rad = rc.rad_cons_lat * PI / 180.0;
    p.lon_cons_rad = rc.rad_cons_lon * PI / 180.0;
    return p;
}

Referenced by TwoStreamRadiation::advance().

Here is the caller graph for this function:

resolve_surface_albedo_sw()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real resolve_surface_albedo_sw	(	int	i,
		int	j,
		const amrex::Array4< const amrex::Real > *	hetero_alb_sw,
		const TwoStreamParams &	rad_choice,
		bool	has_hetero_alb
	)

Resolve per-column shortwave surface albedo from hetero field or fallback.

Precedence:

1. If hetero_alb_sw array available and value at (i,j) is finite ∈ [0,1], use it
2. Otherwise, use rad_choice.surface_albedo_sw (already clamped by init_params)
3. Hard default: 0.3

Parameters

[in]	i,j	Column index
[in]	hetero_alb_sw	Heterogeneous SW albedo field (may be nullptr)
[in]	rad_choice	Radiation parameters with fallback surface_albedo_sw
[in]	has_hetero_alb	true if hetero_alb_sw is available

Returns

SW albedo in [0, 1]

{
    amrex::Real alb = rad_choice.surface_albedo_sw;  // Default fallback (already clamped)
 
    if (has_hetero_alb && hetero_alb_sw != nullptr && hetero_alb_sw->contains(i, j, 0)) {
        amrex::Real hetero_val = (*hetero_alb_sw)(i, j, 0, 0);
        if (amrex::Math::isfinite(hetero_val) && hetero_val >= 0.0 && hetero_val <= 1.0) {
            alb = hetero_val;
        }
    }
 
    return alb;
}

Referenced by vertical_two_stream_sweep().

Here is the caller graph for this function:

resolve_surface_emissivity_lw()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real resolve_surface_emissivity_lw	(	int	i,
		int	j,
		const amrex::Array4< const amrex::Real > *	hetero_emiss_lw,
		const TwoStreamParams &	rad_choice,
		bool	has_hetero_emiss
	)

Resolve per-column longwave surface emissivity from hetero field or fallback.

Precedence:

1. If hetero_emiss_lw array available and value at (i,j) is finite ∈ [0,1], use it
2. Otherwise, use rad_choice.surface_emissivity_lw (already clamped by init_params)
3. Hard default: 0.99

Parameters

[in]	i,j	Column index
[in]	hetero_emiss_lw	Heterogeneous LW emissivity field (may be nullptr)
[in]	rad_choice	Radiation parameters with fallback surface_emissivity_lw
[in]	has_hetero_emiss	true if hetero_emiss_lw is available

Returns

LW emissivity in [0, 1]

{
    amrex::Real emiss = rad_choice.surface_emissivity_lw;  // Default fallback (already clamped)
 
    if (has_hetero_emiss && hetero_emiss_lw != nullptr && hetero_emiss_lw->contains(i, j, 0)) {
        amrex::Real hetero_val = (*hetero_emiss_lw)(i, j, 0, 0);  // Assume single component
        if (amrex::Math::isfinite(hetero_val) && hetero_val >= 0.0 && hetero_val <= 1.0) {
            emiss = hetero_val;
        }
    }
 
    return emiss;
}

Referenced by vertical_two_stream_sweep().

Here is the caller graph for this function:

resolve_surface_temp_k()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real resolve_surface_temp_k	(	int	i,
		int	j,
		const amrex::Array4< const amrex::Real > *	lsm_t_sfc,
		bool	has_lsm_t_sfc,
		const amrex::Array4< const amrex::Real > *	seb_t_sfc,
		bool	has_seb_t_sfc,
		const amrex::Array4< const amrex::Real > *	surface_layer_theta,
		bool	has_surface_layer,
		const TwoStreamParams &	rad_choice,
		bool &	from_surface_layer
	)

Resolve the per-column surface-temperature boundary condition.

The candidates remain separate until this per-cell resolver runs. This is important for an LSM field that exists globally but contains an undefined value in an individual column.

Precedence for each column:

1. valid LSM absolute temperature;
2. valid prognostic SEB absolute temperature;
3. valid SurfaceLayer potential temperature (marked for Exner conversion);
4. the scalar absolute-temperature fallback, erf.rad_t_sfc.

Only the SurfaceLayer candidate is converted by the caller. The LSM and prognostic SEB candidates already satisfy the absolute-temperature contract.

Parameters

[in]	i,j	Column index
[in]	lsm_t_sfc	LSM absolute-temperature field (may be nullptr)
[in]	has_lsm_t_sfc	true if the LSM field is available
[in]	seb_t_sfc	prognostic SEB absolute-temperature field (may be nullptr)
[in]	has_seb_t_sfc	true if the prognostic SEB field is available
[in]	surface_layer_theta	SurfaceLayer potential-temperature field (may be nullptr)
[in]	has_surface_layer	true if the SurfaceLayer field is available
[in]	rad_choice	Radiation parameters with the fallback t_sfc_default
[out]	from_surface_layer	true when the returned value is SurfaceLayer theta

Returns

Surface temperature candidate [K]

{
    amrex::Real t_surf = rad_choice.t_sfc_default;  // erf.rad_t_sfc (already validated)
    from_surface_layer = false;
 
    if (has_lsm_t_sfc && lsm_t_sfc != nullptr && lsm_t_sfc->contains(i, j, 0)) {
        const amrex::Real lsm_value = (*lsm_t_sfc)(i, j, 0, 0);
        if (valid_surface_temperature(lsm_value)) {
            return lsm_value;
        }
    }
 
    if (has_seb_t_sfc && seb_t_sfc != nullptr && seb_t_sfc->contains(i, j, 0)) {
        const amrex::Real seb_value = (*seb_t_sfc)(i, j, 0, 0);
        if (valid_surface_temperature(seb_value)) {
            return seb_value;
        }
    }
 
    if (has_surface_layer && surface_layer_theta != nullptr &&
        surface_layer_theta->contains(i, j, 0)) {
        const amrex::Real theta = (*surface_layer_theta)(i, j, 0, 0);
        if (valid_surface_temperature(theta)) {
            from_surface_layer = true;
            return theta;
        }
    }
 
    return t_surf;
}

Referenced by vertical_two_stream_sweep().

Here is the call graph for this function:

Here is the caller graph for this function:

select_scattering_props()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void select_scattering_props	(	amrex::Real	z_center,
		const TwoStreamParams &	rad_choice,
		bool	apply_cloud,
		amrex::Real &	omega,
		amrex::Real &	g
	)

GPU-safe helper to select the single-scattering albedo and asymmetry factor to use for level k's diffuse SW calculation.

When this column evaluation applies the cloud-layer enhancement (apply_cloud == true, tau_profile_type == CloudLayer, and level k falls within the cloud band), the cloud scattering properties (cloud_single_scattering_albedo, cloud_asymmetry_factor) are used. Otherwise, the clear-sky scattering properties (single_scattering_albedo, asymmetry_factor) are used. Both default to 0.0, so by default this function always yields omega == 0.0 and the layers neither scatter the direct beam nor reflect diffuse light; the diffuse field then consists of the surface-reflected beam only.

Parameters

[in]	z_center	Height of the layer center above the surface [m].
[in]	rad_choice	Radiation parameters.
[in]	apply_cloud	Same flag passed to tau_layer_value(); true for the cloudy-column evaluation, false for the clear-sky column evaluation.
[out]	omega	Selected single-scattering albedo for this level.
[out]	g	Selected asymmetry factor for this level.

{
    bool use_cloud_props = apply_cloud &&
        rad_choice.tau_profile_type == TauProfileType::CloudLayer &&
        is_cloud_level(z_center, rad_choice);
 
    if (use_cloud_props) {
        omega = rad_choice.cloud_single_scattering_albedo;
        g = rad_choice.cloud_asymmetry_factor;
    } else {
        omega = rad_choice.single_scattering_albedo;
        g = rad_choice.asymmetry_factor;
    }
}

Referenced by diagnose_layer_optics().

Here is the call graph for this function:

Here is the caller graph for this function:

tau_layer_value()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real tau_layer_value	(	amrex::Real	z_center,
		amrex::Real	tau_base,
		const TwoStreamParams &	rad_choice,
		bool	apply_cloud
	)

GPU-safe helper to compute the per-layer optical depth at level k, given the base (clear-sky) optical depth and cloud-layer parameters.

When rad_choice.tau_profile_type == Constant, returns tau_base unchanged (byte-identical). When == CloudLayer, adds rad_choice.cloud_tau_per_layer whenever the level height falls within [cloud_base_height_m, cloud_top_height_m].

Parameters

[in]	z_center	Height of the layer center above the surface [m].
[in]	tau_base	Clear-sky optical depth per layer.
[in]	rad_choice	Radiation parameters.
[in]	apply_cloud	If false, always returns tau_base (used for the clear-sky column computation even when cloud_fraction > 0).

Returns

Optical depth for this layer [unitless].

{
    if (!apply_cloud || rad_choice.tau_profile_type != TauProfileType::CloudLayer) {
        return tau_base;
    }
    if (is_cloud_level(z_center, rad_choice)) {
        return tau_base + rad_choice.cloud_tau_per_layer;
    }
    return tau_base;
}

Referenced by diagnose_layer_optics(), and diagnose_layer_tau().

Here is the call graph for this function:

Here is the caller graph for this function:

two_stream_scratch_box()

amrex::Box two_stream_scratch_box ( const amrex::Box &  bx )

inline

Box of the scratch FArrayBox a sweep over box bx needs.

{
    amrex::Box sbx(bx);
    sbx.growHi(2, 1);
    return sbx;
}

Referenced by TwoStreamRadiation::advance().

Here is the caller graph for this function:

valid_surface_temperature()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool valid_surface_temperature

(

amrex::Real

value

)

Check whether a heterogeneous absolute surface temperature is valid.

{
    return amrex::Math::isfinite(value) && value > 0.0 && value < lsm_undefined;
}

Referenced by resolve_surface_temp_k().

Here is the caller graph for this function:

vertical_two_stream_sweep()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void vertical_two_stream_sweep	(	int	i,
		int	j,
		const amrex::Box &	bx,
		amrex::Real	dz_uniform,
		const amrex::Array4< const amrex::Real > &	state_arr,
		const TwoStreamParams &	rad_choice,
		bool	cloudy,
		const amrex::Array4< amrex::Real > &	qheating_arr,
		amrex::Real &	max_heating_rate,
		amrex::Real &	sw_surface_flux,
		amrex::Real &	sw_up_toa,
		amrex::Real &	lw_net_surface,
		amrex::Real &	lw_up_toa,
		amrex::Real &	sw_down_toa,
		const amrex::Array4< const amrex::Real > &	z_phys_nd,
		const amrex::Array4< amrex::Real > &	scratch,
		bool	has_hetero_alb_sw = false,
		const amrex::Array4< const amrex::Real > *	hetero_alb_sw = nullptr,
		bool	has_hetero_emiss_lw = false,
		const amrex::Array4< const amrex::Real > *	hetero_emiss_lw = nullptr,
		bool	has_lsm_t_sfc = false,
		const amrex::Array4< const amrex::Real > *	lsm_t_sfc = nullptr,
		bool	has_seb_t_sfc = false,
		const amrex::Array4< const amrex::Real > *	seb_t_sfc = nullptr,
		bool	has_surface_layer = false,
		const amrex::Array4< const amrex::Real > *	surface_layer_theta = nullptr,
		bool	has_latlon = false,
		const amrex::Array4< const amrex::Real > *	lat_arr = nullptr,
		const amrex::Array4< const amrex::Real > *	lon_arr = nullptr,
		const amrex::Array4< amrex::Real > *	rad_flux_out = nullptr
	)

GPU-safe per-column vertical integration kernel for two-stream radiation, computing either the clear-sky or cloudy-column fluxes and per-level heating rates, depending on the cloudy flag.

Vertical orientation. ERF's vertical index increases upward: k = kmin is the layer adjacent to the surface and k = kmax is the layer adjacent to the top of the domain (TOA for this model). Layer k spans the interfaces m = k - kmin (bottom) and m = k - kmin + 1 (top), so interface m = 0 is the surface and m = nlev is the TOA.

Per (i,j) column:

1. SW: Beer-Lambert direct beam on the interfaces, then the two-stream diffuse field (upward and downward streams) by the adding method with the surface albedo as the lower boundary condition; the per-level SW heating rate from the net-flux divergence is written to qheating_arr(i,j,k,0).
2. LW: sweep downward from the TOA (F_down = 0) to the surface, then upward from the surface (F_up = eps * sigma * T_s^4 + (1 - eps) * F_down) to the TOA, storing both interface profiles. The per-level LW heating rate from the net-flux divergence is written to qheating_arr(i,j,k,1).
3. Scalar diagnostics (max heating rate, surface fluxes) are returned for the reduction in the caller.

Layer temperature is obtained from (rho, rho*theta, qv) through the equation of state, i.e. including the Exner function, so LW emission uses absolute temperature rather than potential temperature.

Integrates per-column heterogeneous surface properties (albedo, emissivity, surface temperature) from optional fields with robust fallback to RadChoice scalar parameters.

Parameters

[in]	i,j	Column indices
[in]	bx	Computational box (cell-centered, full vertical extent)
[in]	dz_uniform	Uniform vertical cell size [m], used when z_phys_cc is unavailable. Passed in rather than taken from a Geometry because Geometry::CellSize() is host-only and this runs on device.
[in]	state_arr	Array proxy to state data (read-only)
[in]	rad_choice	Column-kernel parameters (TwoStreamParams, built from RadChoice)
[in]	cloudy	If true and tau_profile_type == CloudLayer, apply the cloud-layer optical depth enhancement (and cloud scattering properties).
[out]	qheating_arr	Component 0 receives the SW and component 1 the LW radiative tendency of potential temperature, dtheta/dt = (dT/dt) / pi [K/s], at every level k in [kmin, kmax]. This is the convention the RhoTheta source term (ERF_MakeSources.cpp) and the RRTMGP path use.
[out]	max_heating_rate	Maximum |Q_sw|+|Q_lw| of the stored theta tendencies in this column
[out]	sw_surface_flux	SW absorbed by the surface, (1 - alb_dir) F_dir + (1 - alb_dif) F_dif [W/m^2]
[out]	sw_up_toa	Upwelling (reflected) SW at the top of the atmosphere [W/m^2]
[out]	lw_net_surface	Net LW (up - down) at the surface [W/m^2]
[out]	lw_up_toa	Upwelling LW at the top of the atmosphere (outgoing longwave) [W/m^2]
[out]	sw_down_toa	Incident SW at the top of the atmosphere, S0 cos(zenith) or zero at night [W/m^2]
[in]	z_phys_nd	Optional nodal physical heights (layer interfaces on a nonuniform grid)
[in]	has_hetero_alb_sw	true if hetero_alb_sw is available
[in]	hetero_alb_sw	Optional per-column SW surface albedo field
[in]	has_hetero_emiss_lw	true if hetero_emiss_lw is available
[in]	hetero_emiss_lw	Optional per-column LW surface emissivity field
[in]	has_lsm_t_sfc	true if the LSM absolute-temperature field is available
[in]	lsm_t_sfc	Optional LSM absolute surface-temperature field [K]
[in]	has_seb_t_sfc	true if the prognostic SEB absolute-temperature field is available
[in]	seb_t_sfc	Optional prognostic SEB absolute surface-temperature field [K]
[in]	has_surface_layer	true if the SurfaceLayer potential-temperature field is available
[in]	surface_layer_theta	Optional SurfaceLayer potential-temperature field [K]
[in]	has_latlon	true if per-column latitude and longitude fields are available
[in]	lat_arr	Optional per-column latitude [degrees] (dynamic sun)
[in]	lon_arr	Optional per-column longitude [degrees, east positive] (dynamic sun)
[out]	rad_flux_out	Optional 4-component flux output in RRTMGP's level layout: at index k the lower interface of layer k, and at kmax + 1 the top-of-atmosphere interface (the array must extend one cell above the column); component 0 is the upward SW, 1 the downward SW (direct plus diffuse), 2 the upward LW and 3 the downward LW [W/m^2]

{
    // Grid bounds: kmin is the surface layer, kmax the top layer.
    const int kmin = bx.smallEnd(2);
    const int kmax = bx.bigEnd(2);
    const int nlev = kmax - kmin + 1;
 
    // The scratch box must hold nlev + 1 interface entries for this column,
    // and so must the flux output when it is wanted.
    AMREX_ASSERT(scratch.contains(i, j, kmin) && scratch.contains(i, j, kmax + 1));
    AMREX_ASSERT(rad_flux_out == nullptr ||
                 (rad_flux_out->contains(i, j, kmin) && rad_flux_out->contains(i, j, kmax + 1)));
 
    // Per-column scratch accessors (see TwoStreamScratch). Layer fields are
    // indexed by layer m in [0, nlev), interface fields by m in [0, nlev].
    auto dz_level        = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::DZ); };
    auto z_level         = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::Z); };
    auto tau_sw_level    = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::TAU_SW); };
    auto F_dir           = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::F_DIR); };
    auto A               = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::A); };
    auto S               = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::S); };
    auto F_lw_up_iface   = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::LW_UP); };
    auto F_lw_down_iface = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::LW_DN); };
    auto tau_lw_level    = [&] (int m) -> amrex::Real& { return scratch(i, j, kmin + m, TwoStreamScratch::TAU_LW); };
 
    // Physical constants
    const amrex::Real sigma  = stefan_boltzmann;  // Stefan-Boltzmann [W/(m^2·K^4)]
    const amrex::Real cp_air = Cp_d;              // Dry-air specific heat at constant pressure [J/(kg·K)]
 
    // ------------------------------------------------------------------
    // Layer thickness dz_level(m) and layer-centre height z_level(m) above
    // the surface. On a stretched or terrain-following grid both come from
    // the interface heights: the height of face k over cell (i,j) is the
    // mean of its four nodes in z_phys_nd, so the thickness is the distance
    // between the layer's two interfaces (not the centre-to-centre spacing,
    // which differs by the stretch ratio) and the top layer gets its own
    // thickness. Without z_phys_nd the spacing is uniform. The heating
    // divergence dF/(rho cp dz) and the mass path rho dz both use this
    // thickness; z_level locates the cloud band and the aerosol profile.
    // ------------------------------------------------------------------
    if (z_phys_nd) {
        auto z_face = [&] (int k) -> amrex::Real {
            return 0.25 * (z_phys_nd(i, j, k) + z_phys_nd(i + 1, j, k) +
                           z_phys_nd(i, j + 1, k) + z_phys_nd(i + 1, j + 1, k));
        };
        const amrex::Real z_surface = z_face(kmin);
        for (int m = 0; m < nlev; ++m) {
            const amrex::Real z_lo = z_face(kmin + m);
            const amrex::Real z_hi = z_face(kmin + m + 1);
            const amrex::Real dz_computed = z_hi - z_lo;
            dz_level(m) = (dz_computed > 0.0 && amrex::Math::isfinite(dz_computed)) ? dz_computed : dz_uniform;
            z_level(m)  = 0.5 * (z_lo + z_hi) - z_surface;
        }
    } else {
        for (int m = 0; m < nlev; ++m) {
            dz_level(m) = dz_uniform;
            z_level(m)  = (static_cast<amrex::Real>(m) + 0.5) * dz_uniform;
        }
    }
 
    // Cosine of the solar zenith angle: the fixed value of
    // erf.fixed_solar_zenith_angle, or the sun's position over this column at
    // the calendar time of the call, with the same formula RRTMGP uses
    // (orbital_cos_zenith without the interval average). The latitude and
    // longitude are the column's own where the grid carries them, else the
    // erf.rad_cons_lat/lon constants.
    amrex::Real cos_zenith;
    if (rad_choice.solar_dynamic) {
        const amrex::Real lat_rad = (has_latlon && lat_arr != nullptr && lat_arr->contains(i, j, 0))
                                        ? (*lat_arr)(i, j, 0) * PI / 180.0 : rad_choice.lat_cons_rad;
        const amrex::Real lon_rad = (has_latlon && lon_arr != nullptr && lon_arr->contains(i, j, 0))
                                        ? (*lon_arr)(i, j, 0) * PI / 180.0 : rad_choice.lon_cons_rad;
        cos_zenith = static_cast<amrex::Real>(
            orbital_cos_zenith_instant(rad_choice.calday, lat_rad, lon_rad, rad_choice.declin));
    } else {
        cos_zenith = rad_choice.cos_zenith_fixed;
    }
 
    // Top-of-atmosphere irradiance of this call (fixed, or 1360.9 W/m^2
    // times the Earth-Sun distance factor of the date; set by the driver).
    const amrex::Real S0 = rad_choice.S0;
    const amrex::Real tau_sw_base = rad_choice.tau_per_layer;
    const amrex::Real tau_lw_base = rad_choice.tau_lw_per_layer;
 
    amrex::Real local_max_heating = 0.0;
 
    // Zero-initialize this column's heating rate output (covers the
    // sw_enabled=false and lw_enabled=false cases), and the interface fluxes
    // when they are wanted (a disabled band or night leaves zeros).
    for (int k = kmin; k <= kmax; ++k) {
        qheating_arr(i, j, k, 0) = 0.0;
        qheating_arr(i, j, k, 1) = 0.0;
    }
    if (rad_flux_out != nullptr) {
        for (int k = kmin; k <= kmax + 1; ++k) {
            for (int comp = 0; comp < 4; ++comp) { (*rad_flux_out)(i, j, k, comp) = 0.0; }
        }
    }
    sw_down_toa = (rad_choice.sw_enabled && cos_zenith > 0.0) ? S0 * cos_zenith : 0.0;
 
    // ========================================================================
    // SHORTWAVE. Interfaces m = 0 (surface) .. nlev (TOA); layer m lies
    // between interfaces m and m+1.
    //  1. Direct beam F_dir(m) by Beer-Lambert from the TOA downward.
    //  2. Adding method, upward pass: A(m) is the albedo of everything below
    //     interface m for diffuse light, S(m) the upward diffuse flux at m
    //     produced by the direct beam illuminating everything below m. The
    //     surface starts the recursion with A(0) = diffuse albedo and
    //     S(0) = direct albedo * F_dir(0).
    //  3. Downward pass from the TOA (no incident diffuse flux): diffuse
    //     downward d(m), diffuse upward u(m) = A(m) d(m) + S(m), and the net
    //     flux F_dir + d - u whose divergence gives the heating rate.
    // ========================================================================
    amrex::Real tau_sw_cum = 0.0;          // Column SW optical depth accumulated from the top
    amrex::Real sw_surface_absorbed = 0.0; // absorbed direct + diffuse at the surface
    amrex::Real sw_up_at_toa = 0.0;        // u(nlev)
    // Direct-beam albedo (per column, LSM field or fallback) and diffuse albedo
    // (surface_albedo_sw_diffuse, or the direct value when negative).
    const amrex::Real alb_dir = resolve_surface_albedo_sw(i, j, hetero_alb_sw, rad_choice, has_hetero_alb_sw);
    const amrex::Real alb_dif = (rad_choice.surface_albedo_sw_diffuse >= 0.0)
                                    ? rad_choice.surface_albedo_sw_diffuse : alb_dir;
 
    if (rad_choice.sw_enabled && cos_zenith > 0.0) {
 
        // 1. Per-layer optical depth and the direct beam on the interfaces.
        F_dir(nlev) = S0 * cos_zenith;
        for (int k = kmax; k >= kmin; --k) {
            const int m = k - kmin;
            amrex::Real omega_m = 0.0, g_m = 0.0;
            diagnose_layer_optics(i, j, k, dz_level(m), z_level(m), state_arr, tau_sw_base,
                                  /*is_sw=*/true, cloudy, rad_choice, tau_sw_level(m), omega_m, g_m);
            tau_sw_cum += tau_sw_level(m);
            F_dir(m) = compute_sw_direct_flux(tau_sw_cum, S0, cos_zenith);
        }
 
        // Two-stream properties of layer m (evaluated in both passes).
        auto layer_props = [&](int m) {
            amrex::Real tau_m = 0.0, omega = 0.0, g = 0.0;
            diagnose_layer_optics(i, j, kmin + m, dz_level(m), z_level(m), state_arr, tau_sw_base,
                                  /*is_sw=*/true, cloudy, rad_choice, tau_m, omega, g);
            return compute_sw_layer_two_stream(tau_sw_level(m), omega, g, cos_zenith);
        };
 
        // 2. Upward pass of the adding method.
        A(0) = alb_dif;
        S(0) = alb_dir * F_dir(0);
        for (int m = 0; m < nlev; ++m) {
            const TwoStreamLayerSW L = layer_props(m);
            amrex::Real denom = 1.0 - L.R_dif * A(m);
            if (denom < 1.0e-12) denom = 1.0e-12;
            A(m + 1) = L.R_dif + L.T_dif * L.T_dif * A(m) / denom;
            S(m + 1) = L.R_dir * F_dir(m + 1)
                     + L.T_dif * (S(m) + A(m) * L.T_dir * F_dir(m + 1)) / denom;
        }
 
        // 3. Downward pass: fluxes on the interfaces and layer heating rates.
        amrex::Real d_above = 0.0;          // No diffuse flux incident at the TOA
        amrex::Real u_above = S(nlev);
        sw_up_at_toa = u_above;
        amrex::Real F_net_above = F_dir(nlev) + d_above - u_above;
        if (rad_flux_out != nullptr) {
            (*rad_flux_out)(i, j, kmax + 1, 0) = u_above;      // SW up at the top of the atmosphere
            (*rad_flux_out)(i, j, kmax + 1, 1) = F_dir(nlev);  // SW down at the top: the incident beam
        }
        for (int k = kmax; k >= kmin; --k) {
            const int m = k - kmin;
            const TwoStreamLayerSW L = layer_props(m);
            amrex::Real denom = 1.0 - L.R_dif * A(m);
            if (denom < 1.0e-12) denom = 1.0e-12;
 
            const amrex::Real d_m = (L.T_dif * d_above + L.T_dir * F_dir(m + 1) + L.R_dif * S(m)) / denom;
            const amrex::Real u_m = A(m) * d_m + S(m);
            const amrex::Real F_net_m = F_dir(m) + d_m - u_m;
            if (rad_flux_out != nullptr) {
                (*rad_flux_out)(i, j, k, 0) = u_m;              // SW up at the lower interface of layer k
                (*rad_flux_out)(i, j, k, 1) = F_dir(m) + d_m;   // SW down, direct plus diffuse
            }
 
            amrex::Real rho = state_arr(i, j, k, Rho_comp);
            if (rho <= 0.0 || !amrex::Math::isfinite(rho)) rho = 1.0;
 
            // dT/dt from the flux divergence, stored as dtheta/dt = (dT/dt)/pi
            // (the RhoTheta source term and RRTMGP use the same convention).
            const amrex::Real Q_sw = compute_sw_heating_rate(F_net_above, F_net_m,
                                                             dz_level(m), rho, cp_air)
                                     * get_inverse_exner(i, j, k, state_arr);
            qheating_arr(i, j, k, 0) = Q_sw;
            local_max_heating = std::max(local_max_heating, std::abs(Q_sw));
 
            d_above = d_m;
            F_net_above = F_net_m;
            if (m == 0) {
                sw_surface_absorbed = (1.0 - alb_dir) * F_dir(0) + (1.0 - alb_dif) * d_m;
            }
        }
    }
 
    // ========================================================================
    // LONGWAVE: interface flux profiles, index m = 0 at the surface and
    // m = nlev at the TOA.
    // ========================================================================
 
    if (rad_choice.lw_enabled) {
        // Downward sweep: TOA -> surface. Also caches the per-layer LW optical
        // depth for the upward sweep.
        F_lw_down_iface(nlev) = 0.0;  // No incoming LW from space
        for (int k = kmax; k >= kmin; --k) {
            const int m = k - kmin;
 
            amrex::Real rho       = state_arr(i, j, k, Rho_comp);
            amrex::Real rho_theta = state_arr(i, j, k, RhoTheta_comp);
            if (rho <= 0.0 || !amrex::Math::isfinite(rho)) rho = 1.0;
            if (rho_theta <= 0.0 || !amrex::Math::isfinite(rho_theta)) rho_theta = 288.15;
            amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
            amrex::Real T_layer = get_temperature_from_rhotheta(rho_theta, rho, qv);
 
            tau_lw_level(m) = diagnose_layer_tau(i, j, k, dz_level(m), z_level(m),
                                                 state_arr, tau_lw_base, /*is_sw=*/false, cloudy,
                                                 rad_choice);
 
            F_lw_down_iface(m) = compute_lw_flux_down(F_lw_down_iface(m + 1), T_layer, sigma, tau_lw_level(m));
        }
 
        // Upward sweep: surface -> TOA. The surface emits eps * sigma * T_s^4
        // and reflects the fraction (1 - eps) of the downwelling flux that
        // reaches it (gray surface, Kirchhoff's law).
        {
            bool t_sfc_from_surface_layer = false;
            amrex::Real t_surface = resolve_surface_temp_k(
                i, j, lsm_t_sfc, has_lsm_t_sfc, seb_t_sfc, has_seb_t_sfc,
                surface_layer_theta, has_surface_layer, rad_choice,
                t_sfc_from_surface_layer);
            // The surface layer carries a potential temperature (MOST works in
            // theta); the emission needs the temperature, so convert with the
            // physical surface pressure diagnosed from the lowest atmospheric
            // cell.
            if (t_sfc_from_surface_layer) {
                const amrex::Real rho = state_arr(i, j, kmin, Rho_comp);
                const amrex::Real rho_theta = state_arr(i, j, kmin, RhoTheta_comp);
                const amrex::Real qv = get_qv_from_state(i, j, kmin, state_arr);
                const amrex::Real delta_z = z_level(0);
                const amrex::Real pressure = erf_surface_temperature::pressure_at_surface(
                    rho, rho_theta, qv, delta_z);
                if (!erf_surface_temperature::theta_to_temperature(
                        t_surface, pressure, rad_choice.rdOcp, t_surface)) {
                    t_surface = rad_choice.t_sfc_default;
                }
            }
            amrex::Real emiss_lw  = resolve_surface_emissivity_lw(i, j, hetero_emiss_lw, rad_choice, has_hetero_emiss_lw);
            F_lw_up_iface(0) = emiss_lw * compute_thermal_intensity(t_surface, sigma)
                             + (1.0 - emiss_lw) * F_lw_down_iface(0);
        }
        for (int k = kmin; k <= kmax; ++k) {
            const int m = k - kmin;
 
            amrex::Real rho       = state_arr(i, j, k, Rho_comp);
            amrex::Real rho_theta = state_arr(i, j, k, RhoTheta_comp);
            if (rho <= 0.0 || !amrex::Math::isfinite(rho)) rho = 1.0;
            if (rho_theta <= 0.0 || !amrex::Math::isfinite(rho_theta)) rho_theta = 288.15;
            amrex::Real qv = get_qv_from_state(i, j, k, state_arr);
            amrex::Real T_layer = get_temperature_from_rhotheta(rho_theta, rho, qv);
 
            F_lw_up_iface(m + 1) = compute_lw_flux_up(F_lw_up_iface(m), T_layer, sigma, tau_lw_level(m));
        }
 
        // Per-level LW heating rate from the net-flux divergence across
        // each layer (bottom interface m, top interface m+1).
        for (int k = kmin; k <= kmax; ++k) {
            const int m = k - kmin;
            amrex::Real rho = state_arr(i, j, k, Rho_comp);
            if (rho <= 0.0 || !amrex::Math::isfinite(rho)) rho = 1.0;
            if (rad_flux_out != nullptr) {
                (*rad_flux_out)(i, j, k, 2) = F_lw_up_iface(m);    // LW up at the lower interface of layer k
                (*rad_flux_out)(i, j, k, 3) = F_lw_down_iface(m);  // LW down
                if (k == kmax) {
                    (*rad_flux_out)(i, j, kmax + 1, 2) = F_lw_up_iface(nlev);    // outgoing LW at the top
                    (*rad_flux_out)(i, j, kmax + 1, 3) = F_lw_down_iface(nlev);  // zero: nothing comes in from space
                }
            }
 
            amrex::Real F_net_top = F_lw_up_iface(m + 1) - F_lw_down_iface(m + 1);
            amrex::Real F_net_bot = F_lw_up_iface(m)     - F_lw_down_iface(m);
 
            // dT/dt from the net-flux divergence, stored as dtheta/dt = (dT/dt)/pi.
            amrex::Real Q_lw = compute_lw_heating_rate(F_net_top, F_net_bot, dz_level(m), rho, cp_air)
                               * get_inverse_exner(i, j, k, state_arr);
            qheating_arr(i, j, k, 1) = Q_lw;
 
            amrex::Real Q_sw_here = qheating_arr(i, j, k, 0);
            local_max_heating = std::max(local_max_heating, std::abs(Q_sw_here) + std::abs(Q_lw));
        }
    }
 
    // ========================================================================
    // SURFACE AND DIAGNOSTICS
    // ========================================================================
    if (rad_choice.sw_enabled) {
        // Absorbed direct plus diffuse (see the SW block above)
        sw_surface_flux = sw_surface_absorbed;
        sw_up_toa = sw_up_at_toa;
    } else {
        sw_surface_flux = 0.0;
        sw_up_toa = 0.0;
    }
 
    if (rad_choice.lw_enabled) {
        amrex::Real F_lw_up_sfc   = F_lw_up_iface(0);
        amrex::Real F_lw_down_sfc = F_lw_down_iface(0);
        lw_up_toa = F_lw_up_iface(nlev);
        lw_net_surface = F_lw_up_sfc - F_lw_down_sfc;
    } else {
        lw_net_surface = 0.0;
        lw_up_toa = 0.0;
    }
 
    max_heating_rate = local_max_heating;
}

Referenced by TwoStreamRadiation::advance().

Here is the call graph for this function:

Here is the caller graph for this function: