TwoStreamRadiation Class Reference

The two-stream radiation model (erf.radiation_model = TwoStream). More...

#include <ERF_TwoStreamRadiation.H>

Collaboration diagram for TwoStreamRadiation:

Classes
struct	FluxDiag
struct	OrbitalCache

Public Member Functions
	TwoStreamRadiation ()=default
void	resize (int nlevs_max)
	Size the per-level containers; no allocation happens here. More...
bool	active () const
	True once a level has been defined with the model selected. More...
const amrex::MultiFab *	prognostic_surface_temperature_state (int lev) const noexcept
	Read-only view of the per-level prognostic surface-temperature state. More...
void	define_level (int lev, const RadChoice &rad_choice, amrex::Real rdOcp, const amrex::BoxArray &ba2d, const amrex::DistributionMapping &dm, const amrex::BoxArray &ba, const amrex::Box &domain, bool supply_land_forcing=false)
	Allocate the 2D surface fields of one level on the given horizontal BoxArray and fill them with the scalar defaults. A no-op unless rad_choice.rad_type == TwoStream. More...
void	advance (int lev, int nstep, amrex::Real time, amrex::Real dt_step, const std::string &call_site, const amrex::MultiFab &cons_old, const amrex::MultiFab *z_phys_nd, const amrex::Geometry &geom, LandSurface &lsm, const amrex::Vector< const amrex::MultiFab * > &radiation_inputs, bool noahmp_active, amrex::MultiFab *qheating, amrex::MultiFab *rad_fluxes, const amrex::MultiFab *t_surf, const amrex::MultiFab *sfc_sens_flux, const amrex::MultiFab *sfc_laten_flux, const amrex::MultiFab *lat_m, const amrex::MultiFab *lon_m, double epoch_time, bool have_datetime)
	One radiation call at one level. More...
amrex::MultiFab *	seb_t_sfc (int lev)
	The prognostic surface state of one level, or nullptr when the SEB is off or the level has not been built. More...
amrex::MultiFab *	seb_q_sfc (int lev)
const amrex::MultiFab *	seb_hfx (int lev) const
	The sensible and latent heat fluxes [W/m^2] the surface energy balance of one level last used, or nullptr when the SEB is off or the level has not been built. Read-only: for output. More...
const amrex::MultiFab *	seb_lh (int lev) const
const amrex::MultiFab *	seb_surface_moisture (int lev)
	What level lev's surface hands the surface layer for its surface mixing ratio (erf.radiation.seb_surface_layer_uses_moisture), per column, in the components of erf_surface_moisture::Component: the soil-water factor (seb_moisture_availability of the current q_s); with erf.radiation.seb_soil_type, the bare-soil resistance and the relative humidity of its pore air at the skin temperature; and with erf.radiation.seb_vegetation_type as well, the vegetated fraction and the canopy resistance without its vapour-deficit factor (from the net shortwave the balance holds divided by (1 - surface albedo) – the sweep's with seb_use_radiation_fluxes, seb_sw_flux_default otherwise – and the skin temperature). Absent pieces are a fraction of 0 and resistances of 1e6 s/m. nullptr when the option is off or the level has not been built. Recomputed on every call, on the valid cells; the halo is not filled. More...
bool	seb_prognostic_active () const
	True when the force-restore state is being evolved and so needs transferring. More...
bool	supplies_land_forcing (int lev) const noexcept
	True when define_level allocated the land-model forcing fields on lev. More...
const amrex::MultiFab *	land_forcing_sw_dn (int lev) const noexcept
	The surface forcing of the last sweep on lev (2D: i,j), or nullptr when the level does not supply it: total downwelling shortwave [W/m^2], downwelling longwave [W/m^2] and the cosine of the solar zenith angle floored at zero [1]. lsm_undefined until a sweep has run on the level, so a copy made before one is caught by the land model's missing-input check instead of reading as zero. More...
const amrex::MultiFab *	land_forcing_lw_dn (int lev) const noexcept
const amrex::MultiFab *	land_forcing_cos_zenith (int lev) const noexcept
void	write_land_forcing (int lev, amrex::MultiFab *sw_flux_dn, amrex::MultiFab *lw_flux_dn, amrex::MultiFab *cos_zenith_angle) const
	Copy the surface forcing of the last sweep on lev into a land-surface model's radiation inputs (SurfaceModel's "sw_flux_dn", "lw_flux_dn" and "cos_zenith_angle" outputs). A null destination is skipped. Each destination must lie on this level's horizontal grids; see copy_surface_plane. Before the level's first sweep this copies the lsm_undefined sentinel, not zero. More...
void	supply_canopy_forcing (int lev, int m_top)
	Have every later sweep on lev write what the building faces need besides the interface fluxes (erf.ibseb.radiation = two_stream; ERF_TwoStreamCanopyForcing.H): the direct beam on a horizontal surface at the interfaces 0 .. m_top of each column, counted from the surface, and the cosine of the zenith. Allocates one field on the level's horizontal grids given those interfaces in z and one 2D field. Call after define_level(); a later define_level() of the level (a regrid) drops the request. More...
bool	supplies_canopy_forcing (int lev) const noexcept
	True when supply_canopy_forcing() was called for lev since its last define_level(). More...
TwoStreamCanopyView	canopy_forcing (int lev, int nstep) const noexcept
	The beam and the cosine of the zenith of lev as the sweep of step nstep wrote them (the view's fluxes are left for the caller, who owns rad_fluxes), or an empty view when the level does not supply them or has not swept at that step: the faces must not take a previous step's radiation, nor the zeros of a level whose radiation came from its parent. More...
void	write_checkpoint (int lev, const std::string &checkpointname) const
	Write the prognostic surface state of one level (SEB on). More...
void	read_checkpoint (int lev, const std::string &restart_chkfile, const amrex::Periodicity &period)
	Read it back when the checkpoint carries it; else keep defaults. More...

Private Attributes
const RadChoice *	m_rad = nullptr
amrex::Real	m_rdOcp = RdoCp
OrbitalCache	m_orbit
amrex::Vector< std::unique_ptr< RadiationDiagnostics > >	m_diag
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_alb_sw
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_emiss_lw
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_t_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_sw_flux_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_lw_flux_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_hfx_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_lh_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_grdflx_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_q_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_t_deep
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_q_deep
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_moisture
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_sw_dn_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_lw_dn_sfc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_cos_zenith
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_canopy
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_canopy_cosz
amrex::Vector< int >	m_canopy_m
amrex::Vector< int >	m_canopy_step
amrex::Vector< FluxDiag >	m_flux_diag

Detailed Description

The two-stream radiation model (erf.radiation_model = TwoStream).

Owns everything the model keeps between steps: the per-level surface property fields it falls back on when no land-surface model supplies them, the surface-energy-balance (SEB) fields and force-restore state, and the domain-mean flux diagnostics of the last column sweep. ERF holds one instance and calls it only after testing solverChoice.rad_type == RadiationType::TwoStream; every method here checks the same predicate and returns without acting otherwise.

Lifecycle:

• resize(nlevs) once at construction (before the inputs are read);
• define_level(lev, radChoice, rdOcp, ba2d, dm, ba, domain, supply_land_forcing) whenever a level is (re)built;
• advance(...) twice per step: "pre_dycore" runs the column sweep on the old state and writes the heating rates, "post_dycore" reuses that sweep's flux diagnostics and advances the SEB state by dt_step;
• write_land_forcing(...) after a sweep, when a land-surface model integrates on the radiation this model computes (see supply_land_forcing);
• supply_canopy_forcing(...) once, when the building faces of the immersed-boundary balance take their radiation from the sweep (erf.ibseb.radiation = two_stream), and canopy_forcing(...) after each sweep;
• write_checkpoint / read_checkpoint for the prognostic surface state.

The heating-rate MultiFab itself stays in ERF, since RRTMGP writes the same 2-component (SW, LW) array and the RhoTheta source term and the plotfile read it without caring which model filled it.

Constructor & Destructor Documentation

TwoStreamRadiation()

TwoStreamRadiation::TwoStreamRadiation ( )

default

Member Function Documentation

active()

bool TwoStreamRadiation::active ( ) const

inline

True once a level has been defined with the model selected.

                                       {
        return (m_rad != nullptr) && m_rad->enabled;
    }

Referenced by seb_hfx(), seb_lh(), seb_prognostic_active(), seb_q_sfc(), and seb_t_sfc().

Here is the caller graph for this function:

advance()

void TwoStreamRadiation::advance	(	int	lev,
		int	nstep,
		amrex::Real	time,
		amrex::Real	dt_step,
		const std::string &	call_site,
		const amrex::MultiFab &	cons_old,
		const amrex::MultiFab *	z_phys_nd,
		const amrex::Geometry &	geom,
		LandSurface &	lsm,
		const amrex::Vector< const amrex::MultiFab * > &	radiation_inputs,
		bool	noahmp_active,
		amrex::MultiFab *	qheating,
		amrex::MultiFab *	rad_fluxes,
		const amrex::MultiFab *	t_surf,
		const amrex::MultiFab *	sfc_sens_flux,
		const amrex::MultiFab *	sfc_laten_flux,
		const amrex::MultiFab *	lat_m,
		const amrex::MultiFab *	lon_m,
		double	epoch_time,
		bool	have_datetime
	)

One radiation call at one level.

Parameters

lev	level
nstep	step number, used as the diagnostics row index
time	time logged with the diagnostics row [s]
dt_step	step size, used by the SEB Euler update [s]
call_site	"pre_dycore" (column sweep + diagnostics) or "post_dycore" (cached diagnostics + SEB update)
cons_old	old-state conserved variables of this level
z_phys_nd	nodal heights (layer interfaces), or nullptr on a uniform grid
geom	geometry of this level
lsm	land-surface model (queried for per-column fields)
radiation_inputs	canonical surface-property inputs, in fixed order
noahmp_active	whether Noah-MP diagnostics may be queried
qheating	2-component (SW, LW) heating-rate MultiFab of this level, written at the pre-dycore call; may be nullptr
rad_fluxes	4-component (SW up, SW down, LW up, LW down) interface fluxes of this level in RRTMGP's level layout (index k the lower interface of layer k, the top of the atmosphere in the z-ghost cell above the top layer), written at the pre-dycore call; may be nullptr
t_surf	the surface layer's surface temperature, or nullptr; used as the longwave boundary condition where no land-surface field exists, as RRTMGP does
sfc_sens_flux	the conservative heat flux (rho w'theta', z faces) the zlo surface layer applies to this level, or nullptr when it applies none; the surface energy balance's H where no land-surface model supplies one (see ERF_SEBTurbulentFlux.H)
sfc_laten_flux	the moisture flux (rho w'qv', z faces) likewise, for LE; nullptr without a moisture model
lat_m	per-column latitude [degrees], or nullptr (then erf.rad_cons_lat); only the WRF/metgrid inits fill it
lon_m	per-column longitude [degrees], likewise
epoch_time	seconds since 1970-01-01 00:00 UTC of this call (start_time + time), the RRTMGP time base; double, as RRTMGP keeps it, since a float cannot resolve a time of day at ~1.7e9 s
have_datetime	true when start_datetime (or a WRF/metgrid start date) fixed the calendar; without it a dynamic sun cannot be placed and the call aborts

{
    BL_PROFILE("TwoStreamRadiation::advance()");
 
    // Only proceed if TwoStream radiation is enabled
    if (!active()) { return; }
    const RadChoice& rad_choice = *m_rad;
 
    // A level on which no box holds a whole column has nothing for the sweep or the surface
    // state to work on; ERF::advance_radiation interpolates its fields from the parent
    // instead. Returning here covers the post-dycore call too, which comes from ERF::Advance
    // and does not test for it. Per box, and identical to level_needs_interpolation() there:
    // if the two tests disagreed, one of them would try to sweep a level the other had
    // already interpolated.
    if (lev > 0) {
        const BoxArray& lev_ba = cons_old.boxArray();
        const Box& dom = geom.Domain();
        for (int ibox = 0; ibox < lev_ba.size(); ++ibox) {
            const Box& b = lev_ba[ibox];
            if (b.smallEnd(2) != dom.smallEnd(2) || b.bigEnd(2) != dom.bigEnd(2)) { return; }
        }
    }
 
    // ---- Contract checks. Each of these would otherwise surface as a wrong
    // heating rate or an out-of-bounds read several routines away.
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(call_site == "pre_dycore" || call_site == "post_dycore",
        "TwoStreamRadiation::advance: call_site must be pre_dycore or post_dycore");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lev < m_alb_sw.size() && m_alb_sw[lev] != nullptr,
        "TwoStreamRadiation::advance called on a level define_level has not built");
    // The sweep indexes the 2D surface fields with the MFIter of the 3D state,
    // so both must have the same number of boxes on the same ranks.
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        m_alb_sw[lev]->boxArray().size() == cons_old.boxArray().size() &&
        m_alb_sw[lev]->DistributionMap() == cons_old.DistributionMap(),
        "TwoStreamRadiation: the 2D surface fields and the 3D state are laid out differently");
    if (qheating != nullptr) {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(qheating->nComp() == 2 &&
                                         qheating->boxArray() == cons_old.boxArray() &&
                                         qheating->DistributionMap() == cons_old.DistributionMap(),
            "TwoStreamRadiation: qheating must be the 2-component (SW, LW) field on the state's grids");
    }
    if (rad_choice.seb_enable) {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_sw_flux_sfc[lev] && m_lw_flux_sfc[lev] && m_hfx_sfc[lev] &&
                                         m_lh_sfc[lev] && m_grdflx_sfc[lev] && m_q_sfc[lev] &&
                                         m_t_deep[lev] && m_q_deep[lev],
            "TwoStreamRadiation: seb_enable is set but the SEB fields were not allocated");
    }
    // The land-model forcing is read off the interface fluxes the sweep writes. Reached
    // only on a level that will sweep: the per-box interpolation return above comes first,
    // and a level whose radiation comes from its parent has no sweep and no rad_fluxes to
    // read here. Keep that order -- TwoStreamRadiationDriver.LandForcingIsUndefinedUntilASweep
    // advances such a level and expects it to return untouched.
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!supplies_land_forcing(lev) || rad_fluxes != nullptr ||
                                     call_site == "post_dycore",
        "TwoStreamRadiation: supplying a land-surface model's forcing needs the rad_fluxes array");
    if (call_site == "post_dycore" && rad_choice.seb_prognostic_enable) {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt_step > 0.0 && std::isfinite(dt_step),
            "TwoStreamRadiation: the force-restore update needs a positive, finite dt_step");
    }
 
    // The surface energy balance runs on every level that sweeps, each evolving its
    // own force-restore state.
    //
    // These fields are per level (see define_level), so two levels covering the same
    // column are not advancing one shared state twice -- they are forming two
    // independent estimates of one surface. Left alone those estimates drift apart,
    // which would give the longwave boundary condition two different values for the
    // same ground. ERF averages t_sfc and q_sfc down after the finer levels have
    // advanced, so a coarse cell carries the mean of its fine children and the two
    // agree where they overlap.
    //
    // A level that does not sweep has no fluxes of its own and returns earlier.
    //
    // NOTE for anyone resolving a merge here: the development side of this line still
    // carries `&& (lev == 0)`, the single-level restriction. Taking that side back --
    // which a "resolve by accepting theirs" does -- silently turns this whole branch
    // into a no-op: the surface energy balance returns to level 0 only while the tests,
    // docs and PR all still claim otherwise. Keep this form.
    const bool seb_active = rad_choice.seb_enable;
    // The force-restore state is only the longwave boundary when no external
    // surface-temperature provider owns that boundary on this level. The
    // canonical radiation input covers SurfaceModel providers such as SLM;
    // Noah-MP can also supply t_sfc directly when no canonical input is set.
    // Orthogonal to the level question above: this asks WHO owns the boundary,
    // that one asks WHETHER this level sweeps at all.
    const bool has_external_surface_temperature =
        (radiation_inputs.size() > 0 && radiation_inputs[0] != nullptr) ||
        (noahmp_active && lsm_has_field(lsm, lev, "t_sfc"));
    // The column kernel would substitute placeholders (rho = 1, rho*theta of
    // 288 K) for a non-finite or non-positive density or rho*theta and carry
    // on, hiding a corrupt state behind plausible heating rates. Refuse such
    // a state here instead.
    //
    // This runs every step, so it is one pass over the state and one
    // collective, not the three that contains_nan() plus two MultiFab::min()
    // calls would cost. A non-finite value is mapped to -infinity so that the
    // same minimum answers both questions: -infinity means non-finite, and any
    // other value <= 0 means non-positive. The reduction is collective, so
    // every rank takes the same branch.
    if (call_site == "pre_dycore") {
        ReduceOps<ReduceOpMin, ReduceOpMin> state_ops;
        ReduceData<amrex::Real, amrex::Real> state_data(state_ops);
        using StateTuple = typename decltype(state_data)::Type;
        for (MFIter mfi(cons_old, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
            const auto& arr = cons_old.const_array(mfi);
            state_ops.eval(mfi.tilebox(), state_data,
                [=] AMREX_GPU_DEVICE (int i, int j, int k) -> StateTuple {
                    constexpr amrex::Real neg_inf = -std::numeric_limits<amrex::Real>::infinity();
                    const amrex::Real rho = arr(i,j,k,Rho_comp);
                    const amrex::Real rth = arr(i,j,k,RhoTheta_comp);
                    return {amrex::Math::isfinite(rho) ? rho : neg_inf,
                            amrex::Math::isfinite(rth) ? rth : neg_inf};
                });
        }
        auto state_tuple = state_data.value(state_ops);
        amrex::Real mins[2] = {amrex::get<0>(state_tuple), amrex::get<1>(state_tuple)};
        ParallelDescriptor::ReduceRealMin(mins, 2);
        const amrex::Real rho_min = mins[0];
        const amrex::Real rth_min = mins[1];
        if (!std::isfinite(rho_min) || !std::isfinite(rth_min)) {
            amrex::Abort("TwoStreamRadiation: the state handed to the column sweep at level " +
                         std::to_string(lev) + ", step " + std::to_string(nstep) +
                         " has a non-finite density or rho*theta");
        }
        if (!(rho_min > 0.0) || !(rth_min > 0.0)) {
            amrex::Abort("TwoStreamRadiation: the state handed to the column sweep at level " +
                         std::to_string(lev) + ", step " + std::to_string(nstep) +
                         " has a non-positive density (min " + std::to_string(rho_min) +
                         ") or rho*theta (min " + std::to_string(rth_min) + ")");
        }
    }
 
    // The column sweep runs once per step, at the pre-dycore call. The
    // post-dycore call would sweep the same old state (vars_old is not
    // swapped until the next step), so it reuses the cached flux diagnostics
    // and only advances the surface-energy-balance state.
    const bool do_sweep = (call_site != "post_dycore");
 
    // One diagnostics writer per level for the life of the run, so the
    // header-written flag and the (step, call_site, time) duplicate guard carry
    // over between calls. Per level rather than one for the run: the levels of a
    // hierarchy report the same step, call_site and time, so a single shared
    // guard would drop every level but the first. All of them append to the one
    // CSV, whose `level` column tells the rows apart.
    if (!m_diag[lev]) {
        m_diag[lev] = std::make_unique<RadiationDiagnostics>(
            rad_choice.verbosity, rad_choice.diag_file, lev,
            rad_choice.diag_enable, rad_choice.diag_stdout_enable,
            rad_choice.diag_tagged_enable, rad_choice.diag_regtest_line_enable,
            rad_choice.diag_csv_enable, rad_choice.diag_callsite_mode,
            rad_choice.diag_dedup_tol);
    }
    RadiationDiagnostics& rad_diag = *m_diag[lev];
 
    // ========================================
    // GPU-Safe ParallelFor Implementation with Cloud Fraction
    // Blending, Diffuse (Scattering) SW Flux, and Per-Level Heating Rate
    // Output (qheating_rates MultiFab)
    // ========================================
 
    // Initialize global diagnostics
    amrex::Real SW_surface = 0.0;
    amrex::Real SW_TOA = 0.0;
    amrex::Real SW_up_TOA = 0.0;
    amrex::Real LW_net_surface = 0.0;
    amrex::Real LW_up_TOA = 0.0;
    amrex::Real heating_rate_max = 0.0;
    amrex::Real seb_residual_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
    amrex::Real seb_residual_max  = std::numeric_limits<amrex::Real>::quiet_NaN();
 
    //  Prognostic SEB surface temperature and moisture diagnostics
    amrex::Real t_s_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
    amrex::Real t_s_max  = std::numeric_limits<amrex::Real>::quiet_NaN();
    amrex::Real q_s_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
    amrex::Real q_s_max  = std::numeric_limits<amrex::Real>::quiet_NaN();
 
    // Get state at this level (conservative variables: density, RhoTheta, etc.)
    const MultiFab& state_cons = cons_old;
 
    // Only compute radiation if we have valid state data
    if (state_cons.nComp() > 0 ) {
 
    // Trivially copyable parameter set for the device lambdas below
    // (RadChoice itself holds std::string members and cannot be captured),
    // with the sun of this call from the inputs shared with RRTMGP. The
    // incident SW at the top (SW_TOA) is a domain mean formed by the sweep,
    // since with a calendar sun it varies across the columns.
    TwoStreamParams ts_params = make_two_stream_params(rad_choice, m_rdOcp);
    if (do_sweep) { set_solar_state(ts_params, rad_choice, m_orbit, epoch_time, have_datetime, lev, nstep); }
 
        // Host-side storage for reduction results (will be set by device-side reduction)
        amrex::Real max_heating_global = 0.0;
        amrex::Real sw_surface_sum = 0.0;
        amrex::Real sw_up_toa_sum = 0.0;
        amrex::Real lw_net_sum = 0.0;
        amrex::Real lw_up_toa_sum = 0.0;
        amrex::Real sw_toa_sum = 0.0;
        amrex::Long n_columns_total = 0;
 
        // SEB residual diagnostics
        amrex::Real seb_residual_sum = 0.0;
        amrex::Long n_seb_columns = 0;
 
        // cloud fraction used to blend clear-sky and cloudy-column results.
        // cloud_fraction == 0.0 (default) means only the clear-sky column is
        // ever evaluated, and the blend below reduces to F = F_clear exactly.
        amrex::Real cloud_fraction = rad_choice.cloud_fraction;
 
        // qheating is ERF's 2-component (SW, LW) heating-rate MultiFab of this
        // level. If it is not allocated yet, the sweep still runs for the
        // diagnostics and skips the per-level heating write.
        MultiFab* qheating_mf = qheating;
 
        // Surface radiative fluxes for the SEB. Precedence: an LSM field when
        // the LSM exposes one; otherwise, with seb_use_radiation_fluxes, the
        // fluxes the column sweep below computes at the surface (written per
        // column by the sweep and left in place for the post-dycore call);
        // otherwise the scalar defaults.
        const bool sw_flux_from_rad = seb_active &&
                                      rad_choice.seb_use_radiation_fluxes &&
                                      !(noahmp_active && lsm_has_field(lsm, lev, "sav"));
        const bool lw_flux_from_rad = seb_active &&
                                      rad_choice.seb_use_radiation_fluxes &&
                                      !(noahmp_active && lsm_has_field(lsm, lev, "fira"));
 
        if (seb_active) {
            fill_or_copy_seb_field(m_alb_sw[lev].get(), lsm, lev, lsm_broadband_albedo_name, rad_choice.surface_albedo_sw);
            fill_or_copy_seb_field(m_emiss_lw[lev].get(), lsm, lev, "sfc_emis", rad_choice.surface_emissivity_lw);
 
            // In prognostic mode, do not seed this state from an LSM field.
            // The column sweep applies external-provider precedence directly,
            // and the post-dycore update advances this state only when
            // TwoStream owns the surface-temperature boundary.
            if (!rad_choice.seb_prognostic_enable) {
                fill_or_copy_seb_field(m_t_sfc[lev].get(), lsm, lev, "t_sfc", rad_choice.rad_t_sfc);
            }
 
            // Net absorbed shortwave: Noah-MP splits it into the canopy (sav)
            // and ground (sag) parts. Net longwave: Noah-MP's fira is the
            // net flux to the atmosphere (positive up); the SEB wants the
            // absorbed flux, so the sign flips.
            if (!sw_flux_from_rad) {
                fill_or_copy_seb_field(m_sw_flux_sfc[lev].get(), lsm, lev, "sav",
                                       rad_choice.seb_sw_flux_default, 1.0, "sag");
            }
            if (!lw_flux_from_rad) {
                fill_or_copy_seb_field(m_lw_flux_sfc[lev].get(), lsm, lev, "fira",
                                       rad_choice.seb_lw_flux_default, -1.0);
            }
            // H and LE: a land-surface field of that name, else the flux the
            // surface layer applies to the air (so the ground loses what the air
            // receives), else the scalar defaults; see ERF_SEBTurbulentFlux.H.
            // No land model exposes "hfx" or "lh" today, so without one the
            // surface layer is the source. G is Noah-MP's grdflx when present.
            //
            // At the post-dycore call, the one the force-restore update uses, the
            // surface layer's flux is the one it applied during this step. At the
            // pre-dycore call it is the previous step's, which only the residual
            // diagnostic reads.
            const int surface_k = geom.Domain().smallEnd(2);
            fill_seb_turbulent_flux(*m_hfx_sfc[lev], SEBTurbulentFlux::Sensible,
                                    lsm_field(lsm, lev, "hfx"), sfc_sens_flux, surface_k,
                                    rad_choice.seb_turbulent_flux_source,
                                    rad_choice.seb_hfx_default);
            fill_seb_turbulent_flux(*m_lh_sfc[lev], SEBTurbulentFlux::Latent,
                                    lsm_field(lsm, lev, "lh"), sfc_laten_flux, surface_k,
                                    rad_choice.seb_turbulent_flux_source,
                                    rad_choice.seb_lh_default);
            fill_or_copy_seb_field(m_grdflx_sfc[lev].get(), lsm, lev, "grdflx", rad_choice.seb_grdflx_default);
 
            // Gate q_sfc fill on prognostic mode: same reasoning as t_sfc.
            if (!rad_choice.seb_prognostic_enable) {
                fill_or_copy_seb_field(m_q_sfc[lev].get(), lsm, lev, "noahmp_water_vapor_mixing_ratio_2m_vegetated", rad_choice.seb_q_sfc_default);
            }
            // No LSM exposes a deep-soil temperature or moisture in kg/kg by
            // name (Noah-MP's smstav / smstot are soil-moisture availability
            // and total column water), so the reservoir values are the
            // scalar defaults.
            m_t_deep[lev]->setVal(rad_choice.seb_t_deep_default);
            m_q_deep[lev]->setVal(rad_choice.seb_q_deep_default);
        }
 
        // The column sweep integrates the whole atmospheric column in one
        // kernel, so it needs every k of a box in a single pass. MFIter tiling
        // would hand it partial columns (the default CPU tile size splits z),
        // so this loop is deliberately untiled and works on valid boxes; the
        // horizontal ParallelFor below still provides the parallelism.
        const Box& rad_domain = geom.Domain();
        if (do_sweep) {
        for (MFIter mfi(state_cons, false); mfi.isValid(); ++mfi)
        {
            const Box& bx = mfi.validbox();
 
            // A box that does not span the domain vertically would give this
            // column solver a truncated atmosphere: no beam from above, no
            // cooling to space. ERF only decomposes in z when max_grid_size_z
            // is smaller than the domain, so refuse that configuration rather
            // than return heating rates that look plausible and are wrong.
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                bx.smallEnd(2) == rad_domain.smallEnd(2) &&
                bx.bigEnd(2)   == rad_domain.bigEnd(2),
                "TwoStream radiation requires grids that span the domain in z; "
                "set amr.max_grid_size_z to at least amr.n_cell in z");
            const auto& state_arr = state_cons.const_array(mfi);
            // Geometry::CellSize() is host-only; read it here and hand the
            // value to the device sweep.
            const amrex::Real dz_uniform_lev = geom.CellSize(2);
 
            // Nodal heights for the layer thicknesses on a stretched or
            // terrain-following grid (nullptr on a uniform grid).
            Array4<const amrex::Real> z_phys_nd_arr;
            if (z_phys_nd != nullptr) {
                z_phys_nd_arr = z_phys_nd->const_array(mfi);
            }
 
            // Per-column scratch of the sweep (interfaces need nlev + 1 entries),
            // shared by the clear and cloudy evaluations of the same column.
            FArrayBox scratch_fab(two_stream_scratch_box(bx), TwoStreamScratch::NCOMP);
            Array4<amrex::Real> scratch_arr = scratch_fab.array();
 
            //  Wire LSM surface property fields or use standalone fallback MultiFabs
            // Priority:
            // 1. If LSM is active (lsm.Get_DataIdx() returns >=0), use real LSM fields
            // 2. Otherwise, use standalone fallback MultiFabs (allocated and constant-filled from RadChoice scalars)
            // The resolve_surface_*() helpers implement the full precedence chain with finite guards
 
            // SW albedo. With Noah-MP, its broadband albedo (see lsm_broadband_albedo_name)
            // comes first: the canonical SurfaceModel input radiation_inputs[2] is the
            // visible direct-beam band of RRTMGP's four (Noah-MP's sfc_alb_dir_vis), not a
            // broadband value. Otherwise that canonical input (SLM), then the standalone
            // field. resolve_surface_albedo_sw takes the value per column where it is in
            // [0, 1] and otherwise the erf.radiation.surface_albedo_sw default, so the
            // lsm_undefined placeholder over water, before the first land step and at
            // night (Noah-MP has no albedo without sunlight) falls back to the default.
            bool has_hetero_alb_sw = false;
            Array4<const amrex::Real> hetero_alb_sw_arr;
            {
                const MultiFab* lsm_albedo = noahmp_active
                                           ? lsm_field(lsm, lev, lsm_broadband_albedo_name)
                                           : nullptr;
                if (lsm_albedo != nullptr) {
                    hetero_alb_sw_arr = lsm_albedo->const_array(mfi);
                    has_hetero_alb_sw = true;
                } else if (radiation_inputs.size() > 2 && radiation_inputs[2]) {
                    hetero_alb_sw_arr = radiation_inputs[2]->const_array(mfi);
                    has_hetero_alb_sw = true;
                } else if (m_alb_sw[lev]) {
                    hetero_alb_sw_arr = m_alb_sw[lev]->const_array(mfi);
                    has_hetero_alb_sw = true;
                }
            }
 
            // LW emissivity: Try LSM field "sfc_emis"
            bool has_hetero_emiss_lw = false;
            Array4<const amrex::Real> hetero_emiss_lw_arr;
            {
                std::string varname_emiss = "sfc_emis";
                int lsm_idx = noahmp_active ? lsm.Get_DataIdx(lev, varname_emiss) : -1;
                if (radiation_inputs.size() > 1 && radiation_inputs[1]) {
                    hetero_emiss_lw_arr = radiation_inputs[1]->const_array(mfi);
                    has_hetero_emiss_lw = true;
                } else if (lsm_idx >= 0) {
                    auto lsm_ptr = lsm.Get_Data_Ptr(lev, lsm_idx);
                    if (lsm_ptr) {
                        hetero_emiss_lw_arr = lsm_ptr->const_array(mfi);
                        has_hetero_emiss_lw = true;
                    }
                } else if (m_emiss_lw[lev]) {
                    hetero_emiss_lw_arr = m_emiss_lw[lev]->const_array(mfi);
                    has_hetero_emiss_lw = true;
                }
            }
 
            // Keep the surface-temperature sources separate until the column
            // kernel resolves them per cell. An LSM field can exist globally
            // while carrying an undefined sentinel in an individual column;
            // that column must still fall through to SurfaceLayer theta.
            // The kernel's precedence is valid LSM absolute temperature,
            // valid prognostic SEB absolute temperature, SurfaceLayer theta,
            // then the scalar absolute-temperature fallback.
            bool has_lsm_t_sfc = false;
            Array4<const amrex::Real> lsm_t_sfc_arr;
            bool has_seb_t_sfc = false;
            Array4<const amrex::Real> seb_t_sfc_arr;
            bool has_surface_layer = false;
            Array4<const amrex::Real> surface_layer_theta_arr;
            {
                std::string varname_t_sfc = "t_sfc";
                int lsm_idx = noahmp_active ? lsm.Get_DataIdx(lev, varname_t_sfc) : -1;
                if (radiation_inputs.size() > 0 && radiation_inputs[0]) {
                    lsm_t_sfc_arr = radiation_inputs[0]->const_array(mfi);
                    has_lsm_t_sfc = true;
                } else if (lsm_idx >= 0) {
                    auto lsm_ptr = lsm.Get_Data_Ptr(lev, lsm_idx);
                    if (lsm_ptr) {
                        lsm_t_sfc_arr = lsm_ptr->const_array(mfi);
                        has_lsm_t_sfc = true;
                    }
                }
                // Prognostic SEB is a level-wide alternative to an authoritative
                // external surface-temperature provider. When one owns this
                // level's boundary, m_t_sfc is not advanced and is not offered
                // as a per-cell fallback.
                if (!has_lsm_t_sfc && rad_choice.seb_prognostic_enable &&
                    seb_active && m_t_sfc[lev]) {
                    seb_t_sfc_arr = m_t_sfc[lev]->const_array(mfi);
                    has_seb_t_sfc = true;
                }
                if (t_surf != nullptr) {
                    surface_layer_theta_arr = t_surf->const_array(mfi);
                    has_surface_layer = true;
                }
            }
 
            // Per-column latitude and longitude for the calendar sun (only
            // the WRF and metgrid initialisations fill them; otherwise the
            // erf.rad_cons_lat/lon constants apply).
            const bool has_latlon = (lat_m != nullptr) && (lon_m != nullptr);
            Array4<const amrex::Real> lat_arr, lon_arr;
            if (has_latlon) {
                lat_arr = lat_m->const_array(mfi);
                lon_arr = lon_m->const_array(mfi);
            }
 
            // Interface fluxes for ERF's rad_fluxes (SW up, SW down, LW up,
            // LW down at each layer's lower interface, plus the top of the
            // atmosphere in the z-ghost cell above the column), the level
            // layout RRTMGP writes. The cloudy evaluation goes to a scratch
            // FArrayBox on the same grown box and is blended like the
            // heating rates.
            const bool write_fluxes = (rad_fluxes != nullptr);
            Array4<amrex::Real> rad_flux_clear_arr;
            FArrayBox rad_flux_cloudy_fab;
            Array4<amrex::Real> rad_flux_cloudy_arr;
            if (write_fluxes) {
                rad_flux_clear_arr = rad_fluxes->array(mfi);
                if (cloud_fraction > 0.0) {
                    rad_flux_cloudy_fab.resize(two_stream_scratch_box(bx), 4);
                    rad_flux_cloudy_arr = rad_flux_cloudy_fab.array();
                }
            }
 
 
            // Surface flux arrays the sweep fills for the SEB when asked to.
            Array4<amrex::Real> sw_sfc_out;
            Array4<amrex::Real> lw_sfc_out;
            if (sw_flux_from_rad) sw_sfc_out = m_sw_flux_sfc[lev]->array(mfi);
            if (lw_flux_from_rad) lw_sfc_out = m_lw_flux_sfc[lev]->array(mfi);
 
            // What a land-surface model integrates on, written by the same kernel: the
            // surface-interface values of the (blended) fluxes, and the sun the sweep used.
            const bool supply_land = supplies_land_forcing(lev);
            Array4<amrex::Real> sw_dn_out, lw_dn_out, coszen_out;
            if (supply_land) {
                sw_dn_out  = m_sw_dn_sfc[lev]->array(mfi);
                lw_dn_out  = m_lw_dn_sfc[lev]->array(mfi);
                coszen_out = m_cos_zenith[lev]->array(mfi);
            }
 
            // What the building faces take (erf.ibseb.radiation = two_stream) besides the
            // interface fluxes, written by the same kernel: the direct beam at the canopy's
            // interfaces from the scratch of each evaluation, weighted as the blend weights
            // them, and the cosine of the zenith the sweep used.
            const bool supply_canopy = supplies_canopy_forcing(lev);
            const int  m_canopy_top = supply_canopy ? m_canopy_m[lev] : 0;
            const bool sw_on = rad_choice.sw_enabled;
            Array4<amrex::Real> canopy_out, canopy_cosz_out;
            if (supply_canopy) {
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(write_fluxes,
                    "TwoStreamRadiation: the canopy forcing of the building faces needs the rad_fluxes array");
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_canopy_top <= bx.length(2),
                    "TwoStreamRadiation: the canopy reaches above the top of the column");
                canopy_out = m_canopy[lev]->array(mfi);
                canopy_cosz_out = m_canopy_cosz[lev]->array(mfi);
            }
 
            // Create a 2D box for (i,j) iteration over the horizontal extent
            // One GPU thread per (i,j) column; k-loop is sequential within each thread
            const auto& lo = bx.loVect();
            const auto& hi = bx.hiVect();
            Box xy_box(IntVect(lo[0], lo[1], 0), IntVect(hi[0], hi[1], 0));
 
            // Count columns in this box for later averaging
            amrex::Long n_cols = static_cast<amrex::Long>(bx.length(0)) *
                                 static_cast<amrex::Long>(bx.length(1));
            n_columns_total += n_cols;
 
            // Clear-sky-column heating rates are written directly
            // into the real qheating_rates MultiFab when available.
            // Fall back to a throwaway local FArrayBox otherwise (keeps the
            // kernel call GPU-safe even if qheating_rates isn't allocated).
            FArrayBox qheating_fallback_fab;
            Array4<amrex::Real> qheating_clear_arr;
            if (qheating_mf != nullptr) {
                qheating_clear_arr = qheating_mf->array(mfi);
            } else {
                qheating_fallback_fab.resize(bx, 2);
                qheating_clear_arr = qheating_fallback_fab.array();
            }
 
            // Cloudy-column heating rates always go into a scratch
            // FArrayBox; only used/blended when cloud_fraction > 0.
            FArrayBox qheating_cloudy_fab;
            Array4<amrex::Real> qheating_cloudy_arr;
            if (cloud_fraction > 0.0) {
                qheating_cloudy_fab.resize(bx, 2);
                qheating_cloudy_arr = qheating_cloudy_fab.array();
            }
 
            // GPU-safe reduction using ReduceOps (per-column results aggregated on device)
            amrex::Real max_heating_box = 0.0;
            amrex::Real sw_sum_box = 0.0;
            amrex::Real sw_up_sum_box = 0.0;
            amrex::Real lw_sum_box = 0.0;
            amrex::Real lw_up_sum_box = 0.0;
            amrex::Real sw_toa_sum_box = 0.0;
 
            // Device-side reduction: max heating, sums of the surface and
            // top-of-atmosphere fluxes
            ReduceOps<ReduceOpMax, ReduceOpSum, ReduceOpSum, ReduceOpSum, ReduceOpSum, ReduceOpSum> reduce_ops;
            ReduceData<amrex::Real, amrex::Real, amrex::Real, amrex::Real, amrex::Real, amrex::Real> reduce_data(reduce_ops);
 
            using ReduceTuple = typename decltype(reduce_data)::Type;
 
         // Launch parallel kernel over (i,j) columns
            reduce_ops.eval(xy_box, reduce_data,
                [=] AMREX_GPU_DEVICE (int i, int j, int /*k_unused*/) -> ReduceTuple
                {
                    // Clear-sky column (always evaluated; this is the sole
                    // contributor when cloud_fraction == 0.0, matching earlier behavior)
                    amrex::Real max_heating_clear = 0.0;
                    amrex::Real sw_flux_clear = 0.0;
                    amrex::Real sw_up_clear = 0.0;
                    amrex::Real lw_net_clear = 0.0;
                    amrex::Real lw_up_clear = 0.0;
                    amrex::Real sw_toa_clear = 0.0;
                    amrex::Real cos_zenith_col = 0.0;
                    vertical_two_stream_sweep(
                        i, j, bx, dz_uniform_lev, state_arr, ts_params, /*cloudy=*/false,
                        qheating_clear_arr,
                        max_heating_clear, sw_flux_clear, sw_up_clear, lw_net_clear, lw_up_clear,
                        sw_toa_clear,
                        z_phys_nd_arr, scratch_arr,
                        has_hetero_alb_sw, &hetero_alb_sw_arr,
                        has_hetero_emiss_lw, &hetero_emiss_lw_arr,
                        has_lsm_t_sfc, &lsm_t_sfc_arr,
                        has_seb_t_sfc, &seb_t_sfc_arr,
                        has_surface_layer, &surface_layer_theta_arr,
                        has_latlon, &lat_arr, &lon_arr,
                        write_fluxes ? &rad_flux_clear_arr : nullptr,
                        &cos_zenith_col);
 
                    // The scratch holds this evaluation's beam until the cloudy one overwrites it.
                    if (supply_canopy) {
                        two_stream_canopy_beam(i, j, bx.smallEnd(2), m_canopy_top, scratch_arr,
                                               TwoStreamScratch::F_DIR, sw_on, cos_zenith_col,
                                               cloud_fraction, false, canopy_out);
                    }
 
                    amrex::Real max_heating_col = max_heating_clear;
                    amrex::Real sw_flux_col = sw_flux_clear;
                    amrex::Real sw_up_col = sw_up_clear;
                    amrex::Real lw_net_col = lw_net_clear;
                    amrex::Real lw_up_col = lw_up_clear;
 
                    // Cloudy column only needs to be evaluated when there is a
                    // nonzero cloud fraction; this keeps the cloud_fraction==0
                    // path numerically and computationally identical.
                    if (cloud_fraction > 0.0) {
                         amrex::Real max_heating_cloudy = 0.0;
                         amrex::Real sw_flux_cloudy = 0.0;
                         amrex::Real sw_up_cloudy = 0.0;
                         amrex::Real lw_net_cloudy = 0.0;
                         amrex::Real lw_up_cloudy = 0.0;
                         amrex::Real sw_toa_cloudy = 0.0;
                         vertical_two_stream_sweep(
                            i, j, bx, dz_uniform_lev, state_arr, ts_params, /*cloudy=*/true,
                            qheating_cloudy_arr,
                            max_heating_cloudy, sw_flux_cloudy, sw_up_cloudy, lw_net_cloudy, lw_up_cloudy,
                            sw_toa_cloudy,
                            z_phys_nd_arr, scratch_arr,
                            has_hetero_alb_sw, &hetero_alb_sw_arr,
                            has_hetero_emiss_lw, &hetero_emiss_lw_arr,
                            has_lsm_t_sfc, &lsm_t_sfc_arr,
                            has_seb_t_sfc, &seb_t_sfc_arr,
                            has_surface_layer, &surface_layer_theta_arr,
                            has_latlon, &lat_arr, &lon_arr,
                            write_fluxes ? &rad_flux_cloudy_arr : nullptr);
 
                        // Blend clear-sky and cloudy-column results
                        sw_flux_col = (1.0 - cloud_fraction) * sw_flux_clear +
                                      cloud_fraction * sw_flux_cloudy;
                        sw_up_col = (1.0 - cloud_fraction) * sw_up_clear +
                                    cloud_fraction * sw_up_cloudy;
                        lw_net_col = (1.0 - cloud_fraction) * lw_net_clear +
                                     cloud_fraction * lw_net_cloudy;
                        lw_up_col = (1.0 - cloud_fraction) * lw_up_clear +
                                    cloud_fraction * lw_up_cloudy;
                        max_heating_col = std::max(max_heating_clear, max_heating_cloudy);
                        if (supply_canopy) {
                            two_stream_canopy_beam(i, j, bx.smallEnd(2), m_canopy_top, scratch_arr,
                                                   TwoStreamScratch::F_DIR, sw_on, cos_zenith_col,
                                                   cloud_fraction, true, canopy_out);
                        }
 
                        // Blend per-level heating rates in place
                        // into qheating_clear_arr (which is the real output
                        // MultiFab when qheating_mf != nullptr).
                        int kmin = bx.smallEnd(2);
                        int kmax = bx.bigEnd(2);
                        for (int k = kmin; k <= kmax; ++k) {
                            for (int comp = 0; comp < 2; ++comp) {
                                amrex::Real q_clear_val = qheating_clear_arr(i, j, k, comp);
                                amrex::Real q_cloudy_val = qheating_cloudy_arr(i, j, k, comp);
                                qheating_clear_arr(i, j, k, comp) =
                                    (1.0 - cloud_fraction) * q_clear_val +
                                    cloud_fraction * q_cloudy_val;
                            }
                        }
                        if (write_fluxes) {
                            // nlev + 1 levels: the top interface sits at kmax + 1
                            for (int k = kmin; k <= kmax + 1; ++k) {
                                for (int comp = 0; comp < 4; ++comp) {
                                    rad_flux_clear_arr(i, j, k, comp) =
                                        (1.0 - cloud_fraction) * rad_flux_clear_arr(i, j, k, comp) +
                                        cloud_fraction * rad_flux_cloudy_arr(i, j, k, comp);
                                }
                            }
                        }
                    }
 
                    // Surface fluxes for the SEB: absorbed shortwave, and the
                    // absorbed longwave, which is minus the net (up - down).
                    if (sw_flux_from_rad) sw_sfc_out(i, j, 0) = sw_flux_col;
                    if (lw_flux_from_rad) lw_sfc_out(i, j, 0) = -lw_net_col;
 
                    // The land model's forcing, after the clear/cloudy blend above.
                    if (supply_land) {
                        two_stream_land_forcing(i, j, bx.smallEnd(2), rad_flux_clear_arr,
                                                cos_zenith_col, sw_dn_out(i, j, 0),
                                                lw_dn_out(i, j, 0), coszen_out(i, j, 0));
                    }
 
                    // The sun the building faces must share (their beam is in canopy_out).
                    if (supply_canopy) {
                        canopy_cosz_out(i, j, 0) = (cos_zenith_col > amrex::Real(0.0)) ? cos_zenith_col : amrex::Real(0.0);
                    }
 
                    // The incident TOA flux is the same for both evaluations.
                    // Return tuple for reduction
                    return {max_heating_col, sw_flux_col, sw_up_col, lw_net_col, lw_up_col, sw_toa_clear};
                }
            );
 
            // Copy results from device to host
            amrex::Gpu::synchronize();
            auto reduce_tuple = reduce_data.value(reduce_ops);
            max_heating_box = amrex::get<0>(reduce_tuple);
            sw_sum_box = amrex::get<1>(reduce_tuple);
            sw_up_sum_box = amrex::get<2>(reduce_tuple);
            lw_sum_box = amrex::get<3>(reduce_tuple);
            lw_up_sum_box = amrex::get<4>(reduce_tuple);
            sw_toa_sum_box = amrex::get<5>(reduce_tuple);
 
            // Accumulate box results into global results
            max_heating_global = std::max(max_heating_global, max_heating_box);
            sw_surface_sum += sw_sum_box;
            sw_up_toa_sum += sw_up_sum_box;
            lw_net_sum += lw_sum_box;
            lw_up_toa_sum += lw_up_sum_box;
            sw_toa_sum += sw_toa_sum_box;
 
        }
        // Every accumulator above is rank-local. Reduce before forming means
        // and maxima, so the diagnostics describe the whole domain and do
        // not change with the decomposition. These are collective calls; the
        // conditions around them are input-driven and identical on all ranks.
        if (do_sweep) {
            amrex::Real sums[5] = {sw_surface_sum, sw_up_toa_sum, lw_net_sum, lw_up_toa_sum, sw_toa_sum};
            ParallelDescriptor::ReduceRealSum(sums, 5);
            sw_surface_sum = sums[0]; sw_up_toa_sum = sums[1];
            lw_net_sum = sums[2];     lw_up_toa_sum = sums[3];
            sw_toa_sum = sums[4];
            ParallelDescriptor::ReduceLongSum(n_columns_total);
            ParallelDescriptor::ReduceRealMax(max_heating_global);
            // The faces may take this sweep's canopy forcing in this step, and only in it.
            if (supplies_canopy_forcing(lev)) { m_canopy_step[lev] = nstep; }
        }
        } // do_sweep
 
        // The heating rates go straight into the RhoTheta source term, so a
        // non-finite value here becomes a non-finite state one step later.
        if (do_sweep && qheating_mf != nullptr &&
            (qheating_mf->contains_nan(0, 2, 0) || qheating_mf->contains_inf(0, 2, 0))) {
            amrex::Abort("TwoStreamRadiation: non-finite heating rate after the column sweep at level " +
                         std::to_string(lev) + ", step " + std::to_string(nstep));
        }
 
         // Warn if diagnostic is requested but SEB infrastructure isn't enabled
        if (rad_choice.seb_diagnostic_enable && !rad_choice.seb_enable && lev == 0) {
            static bool warned_seb_misconfig = false;
            if (!warned_seb_misconfig && ParallelDescriptor::IOProcessor()) {
                Print() << "WARNING: erf.radiation.seb_diagnostic_enable=true but "
                           "seb_enable=false; SEB residual diagnostics will report NaN. "
                           "Set erf.radiation.seb_enable=true to enable SEB field "
                           "population.\n";
                warned_seb_misconfig = true;
            }
        }
        // Compute SEB residual diagnostics if enabled
        if (rad_choice.seb_diagnostic_enable && seb_active) {
            seb_residual_max = 0.0;
            // Second loop over boxes to compute SEB residual from populated SEB
            // MultiFabs. Untiled: the work below is per surface column, and a
            // tiled iteration would visit each (i,j) once per z tile.
            for (MFIter mfi(state_cons, false); mfi.isValid(); ++mfi) {
                const Box& bx = mfi.validbox();
                const auto& lo = bx.loVect();
                const auto& hi = bx.hiVect();
                Box xy_box(IntVect(lo[0], lo[1], 0), IntVect(hi[0], hi[1], 0));
 
                // Get SEB field arrays
                Array4<const amrex::Real> sw_flux_arr = m_sw_flux_sfc[lev]->const_array(mfi);
                Array4<const amrex::Real> lw_flux_arr = m_lw_flux_sfc[lev]->const_array(mfi);
                Array4<const amrex::Real> hfx_arr = m_hfx_sfc[lev]->const_array(mfi);
                Array4<const amrex::Real> lh_arr = m_lh_sfc[lev]->const_array(mfi);
                Array4<const amrex::Real> grdflx_arr = m_grdflx_sfc[lev]->const_array(mfi);
 
                // Count columns and compute residuals
                amrex::Long n_cols_box = static_cast<amrex::Long>(bx.length(0)) *
                                        static_cast<amrex::Long>(bx.length(1));
                amrex::Real residual_sum_box = 0.0;
                amrex::Real residual_max_box = 0.0;
 
                // GPU-safe reduction for SEB residuals
                ReduceOps<ReduceOpSum, ReduceOpMax> seb_reduce_ops;
                ReduceData<amrex::Real, amrex::Real> seb_reduce_data(seb_reduce_ops);
 
                using SEBReduceTuple = typename decltype(seb_reduce_data)::Type;
 
                seb_reduce_ops.eval(xy_box, seb_reduce_data,
                    [=] AMREX_GPU_DEVICE (int i, int j, int /*k_unused*/) -> SEBReduceTuple {
                        amrex::Real sw_net = sw_flux_arr(i, j, 0);
                        amrex::Real lw_net = lw_flux_arr(i, j, 0);
                        amrex::Real hfx = hfx_arr(i, j, 0);
                        amrex::Real lh = lh_arr(i, j, 0);
                        amrex::Real grdflx = grdflx_arr(i, j, 0);
 
                        // Compute residual using helper function
                        amrex::Real residual = diagnose_seb_residual(sw_net, lw_net, hfx, lh, grdflx);
 
                        // Return sum and abs(max) of residual
                        return {residual, std::abs(residual)};
                    }
                );
 
                // Copy results from device to host
                amrex::Gpu::synchronize();
                auto seb_reduce_tuple = seb_reduce_data.value(seb_reduce_ops);
                residual_sum_box = amrex::get<0>(seb_reduce_tuple);
                residual_max_box = amrex::get<1>(seb_reduce_tuple);
 
                // Accumulate into global results
                seb_residual_sum += residual_sum_box;
                seb_residual_max = std::max(seb_residual_max, residual_max_box);
                n_seb_columns += n_cols_box;
            }
            ParallelDescriptor::ReduceRealSum(seb_residual_sum);
            ParallelDescriptor::ReduceRealMax(seb_residual_max);
            ParallelDescriptor::ReduceLongSum(n_seb_columns);
        }
 
        //  Prognostic SEB surface temperature and moisture evolution
        // Only advance the force-restore state when TwoStream owns the
        // surface-temperature boundary at this level.
        //
        // With the surface radiation taken from the sweep, also wait for this level's
        // first sweep. Until then m_sw_flux_sfc / m_lw_flux_sfc hold what define_level
        // left, and advancing on them would apply the scalar defaults for that step. That
        // is the first step of a level built by interp_atmos_from_coarse, whose pre-dycore
        // sweep ERF::advance_radiation skips (the level's state is not yet consistent),
        // while this post-dycore call still arrives.
        const bool seb_radiation_ready = m_flux_diag[lev].valid ||
                                         !(sw_flux_from_rad || lw_flux_from_rad);
        if (rad_choice.seb_prognostic_enable && seb_active &&
            call_site == "post_dycore" && seb_radiation_ready) {
            if (!has_external_surface_temperature) {
                // No external provider owns the boundary; advance TwoStream's state.
 
                // Initialize diagnostics for T_s and q_s
                amrex::Real t_s_sum = 0.0;
                amrex::Real t_s_max_val = -std::numeric_limits<amrex::Real>::max();
                amrex::Real q_s_sum = 0.0;
                amrex::Real q_s_max_val = -std::numeric_limits<amrex::Real>::max();
                amrex::Long n_prog_columns = 0;
 
                // Third loop over boxes for prognostic SEB update. Untiled for
                // the same reason as above, and here it also matters for
                // correctness: the force-restore update is applied in place, so
                // visiting a column twice would advance it twice in one step.
                for (MFIter mfi(state_cons, false); mfi.isValid(); ++mfi) {
                    const Box& bx = mfi.validbox();
                    const auto& lo = bx.loVect();
                    const auto& hi = bx.hiVect();
                    Box xy_box(IntVect(lo[0], lo[1], 0), IntVect(hi[0], hi[1], 0));
 
                    // Get SEB field arrays (read-only)
                    Array4<const amrex::Real> sw_flux_arr = m_sw_flux_sfc[lev]->const_array(mfi);
                    Array4<const amrex::Real> lw_flux_arr = m_lw_flux_sfc[lev]->const_array(mfi);
                    Array4<const amrex::Real> hfx_arr = m_hfx_sfc[lev]->const_array(mfi);
                    Array4<const amrex::Real> lh_arr = m_lh_sfc[lev]->const_array(mfi);
                    Array4<const amrex::Real> grdflx_arr = m_grdflx_sfc[lev]->const_array(mfi);
                    Array4<const amrex::Real> t_deep_arr = m_t_deep[lev]->const_array(mfi);
                    Array4<const amrex::Real> q_deep_arr = m_q_deep[lev]->const_array(mfi);
 
                    // Get SEB state arrays (read-write for prognostic update)
                    Array4<amrex::Real> t_s_arr = m_t_sfc[lev]->array(mfi);
                    Array4<amrex::Real> q_s_arr = m_q_sfc[lev]->array(mfi);
                    // Count columns and prepare for reductions
                    amrex::Long n_cols_box = static_cast<amrex::Long>(bx.length(0)) *
                                             static_cast<amrex::Long>(bx.length(1));
                    amrex::Real t_s_sum_box = 0.0;
                    amrex::Real t_s_max_box = -std::numeric_limits<amrex::Real>::max();
                    amrex::Real q_s_sum_box = 0.0;
                    amrex::Real q_s_max_box = -std::numeric_limits<amrex::Real>::max();
 
                    // GPU-safe update for prognostic T_s and q_s with reductions
                    ReduceOps<ReduceOpSum, ReduceOpMax, ReduceOpSum, ReduceOpMax> prog_reduce_ops;
                    ReduceData<amrex::Real, amrex::Real, amrex::Real, amrex::Real> prog_reduce_data(prog_reduce_ops);
 
                    using ProgReduceTuple = typename decltype(prog_reduce_data)::Type;
 
                    prog_reduce_ops.eval(xy_box, prog_reduce_data,
                            [=,  C_s=rad_choice.seb_surface_heat_capacity,
                            tau=rad_choice.seb_restore_timescale_s,
                            d_s=rad_choice.seb_moisture_layer_depth_m,
                            tau_q=rad_choice.seb_moisture_restore_timescale_s,
                            t_min=rad_choice.seb_prognostic_t_min_k,
                            t_max=rad_choice.seb_prognostic_t_max_k,
                            q_min=rad_choice.seb_prognostic_q_min,
                            q_max=rad_choice.seb_prognostic_q_max]
                            AMREX_GPU_DEVICE (int i, int j, int /*k_unused*/) -> ProgReduceTuple {
                            amrex::Real t_s_old = t_s_arr(i, j, 0);
                            amrex::Real q_s_old = q_s_arr(i, j, 0);
 
                            // Read forcing data
                            amrex::Real sw_net = sw_flux_arr(i, j, 0);
                            amrex::Real lw_net = lw_flux_arr(i, j, 0);
                            amrex::Real hfx = hfx_arr(i, j, 0);
                            amrex::Real lh = lh_arr(i, j, 0);
                            amrex::Real grdflx = grdflx_arr(i, j, 0);
                            amrex::Real t_deep_val = t_deep_arr(i, j, 0);
                            amrex::Real q_deep_val = q_deep_arr(i, j, 0);
 
                            // Compute SEB residual
                            amrex::Real seb_res = diagnose_seb_residual(sw_net, lw_net, hfx, lh, grdflx);
 
                            // Compute tendencies
                            amrex::Real dT_s_dt = prognostic_dTs_dt(seb_res, t_s_old, t_deep_val,
                                                                     C_s, tau);
                            amrex::Real dq_s_dt = prognostic_dqs_dt(lh, q_s_old, q_deep_val,
                                                                     d_s, tau_q);
 
                            // Explicit Euler update over the step size
                            amrex::Real t_s_new = t_s_old + dt_step * dT_s_dt;
                            amrex::Real q_s_new = q_s_old + dt_step * dq_s_dt;
 
                            // Clamp to valid ranges
                            t_s_new = amrex::max(t_min, amrex::min(t_max, t_s_new));
                            q_s_new = amrex::max(q_min, amrex::min(q_max, q_s_new));
 
                            // Write back updated values (this modifies the device array)
                            t_s_arr(i, j, 0) = t_s_new;
                            q_s_arr(i, j, 0) = q_s_new;
 
                            // Return for reduction: sum T_s, max T_s, sum q_s, max q_s
                            return {t_s_new, std::abs(t_s_new), q_s_new, std::abs(q_s_new)};
                        });
 
                    // Copy results from device to host
                    amrex::Gpu::synchronize();
                    auto prog_reduce_tuple = prog_reduce_data.value(prog_reduce_ops);
                    t_s_sum_box = amrex::get<0>(prog_reduce_tuple);
                    t_s_max_box = amrex::get<1>(prog_reduce_tuple);
                    q_s_sum_box = amrex::get<2>(prog_reduce_tuple);
                    q_s_max_box = amrex::get<3>(prog_reduce_tuple);
 
                    // Accumulate into global results
                    t_s_sum += t_s_sum_box;
                    t_s_max_val = std::max(t_s_max_val, t_s_max_box);
                    q_s_sum += q_s_sum_box;
                    q_s_max_val = std::max(q_s_max_val, q_s_max_box);
                    n_prog_columns += n_cols_box;
                }
 
                // The force-restore kernels return a zero tendency on any
                // non-finite input, so a NaN here means the state itself was
                // corrupted (a bad restart file, or an overwrite elsewhere).
                if (m_t_sfc[lev]->contains_nan(0, 1, 0) || m_q_sfc[lev]->contains_nan(0, 1, 0)) {
                    amrex::Abort("TwoStreamRadiation: non-finite surface temperature or moisture "
                                 "after the force-restore update at level " + std::to_string(lev) +
                                 ", step " + std::to_string(nstep));
                }
 
                {
                    amrex::Real sums[2] = {t_s_sum, q_s_sum};
                    ParallelDescriptor::ReduceRealSum(sums, 2);
                    t_s_sum = sums[0]; q_s_sum = sums[1];
                    ParallelDescriptor::ReduceRealMax(t_s_max_val);
                    ParallelDescriptor::ReduceRealMax(q_s_max_val);
                    ParallelDescriptor::ReduceLongSum(n_prog_columns);
                }
 
                // Compute mean values from sums
                if (n_prog_columns > 0) {
                    t_s_mean = t_s_sum / static_cast<amrex::Real>(n_prog_columns);
                    t_s_max = t_s_max_val;
                    q_s_mean = q_s_sum / static_cast<amrex::Real>(n_prog_columns);
                    q_s_max = q_s_max_val;
                } else {
                    t_s_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
                    t_s_max = std::numeric_limits<amrex::Real>::quiet_NaN();
                    q_s_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
                    q_s_max = std::numeric_limits<amrex::Real>::quiet_NaN();
                }
            } else {
                // An external provider owns the boundary; leave the unused
                // TwoStream state untouched and do not report it as active.
                t_s_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
                t_s_max = std::numeric_limits<amrex::Real>::quiet_NaN();
                q_s_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
                q_s_max = std::numeric_limits<amrex::Real>::quiet_NaN();
            }
        }
        // equivalent to a single-column value for spatially UNIFORM atmospheres
        // (as in the current SW_ClearSky_Analytical / LW_Isothermal RegTests).
        // Cloud layer and scattering tests are ALSO spatially uniform (identical
        // cloud/tau/scattering parameters applied to every column), so the
        // domain-averaged value still equals the true single-column flux there.
        // True horizontal heterogeneity (e.g., patchy clouds varying by column)
        // remains deferred to future work.
        if (do_sweep) {
            if (n_columns_total > 0) {
                const amrex::Real inv_n = 1.0 / static_cast<amrex::Real>(n_columns_total);
                SW_surface     = sw_surface_sum * inv_n;
                SW_TOA         = sw_toa_sum * inv_n;
                SW_up_TOA      = sw_up_toa_sum * inv_n;
                LW_net_surface = lw_net_sum * inv_n;
                LW_up_TOA      = lw_up_toa_sum * inv_n;
            }
            heating_rate_max = max_heating_global;
            m_flux_diag[lev] = FluxDiag{true, SW_surface, SW_TOA, SW_up_TOA,
                                        LW_net_surface, LW_up_TOA,
                                        heating_rate_max};
        } else {
            const FluxDiag& cached = m_flux_diag[lev];
            if (!cached.valid) {
                // No sweep has run on this level yet, so there is nothing to report. This
                // happens on the first step of a level built by interp_atmos_from_coarse:
                // ERF::advance_radiation interpolates that level's fields and returns
                // before the pre-dycore sweep, while this post-dycore call still arrives
                // from ERF::Advance. Emitting here would write a row of zeros that reads
                // exactly like a real zero-flux result. A nested patch already contributes
                // no rows for the same reason; this is the same rule for a level that
                // simply has not swept yet.
                return;
            }
            SW_surface       = cached.SW_surface;
            SW_TOA           = cached.SW_TOA;
            SW_up_TOA        = cached.SW_up_TOA;
            LW_net_surface   = cached.LW_net_surface;
            LW_up_TOA        = cached.LW_up_TOA;
            heating_rate_max = cached.heating_rate_max;
        }
 
        // Compute SEB residual mean from sum
        if (rad_choice.seb_diagnostic_enable && seb_active && n_seb_columns > 0) {
            seb_residual_mean = seb_residual_sum / static_cast<amrex::Real>(n_seb_columns);
        } else {
            // When feature is disabled, use NaN for backward compatibility
            seb_residual_mean = std::numeric_limits<amrex::Real>::quiet_NaN();
            seb_residual_max = std::numeric_limits<amrex::Real>::quiet_NaN();
        }
 
    }
 
    // Logging output
    if (rad_choice.verbosity >= 1 && ParallelDescriptor::IOProcessor()) {
        Print() << "Radiation diagnostics at step " << nstep << ":\n"
                << "  SW TOA = " << SW_TOA << " W/m^2\n"
                << "  SW surface = " << SW_surface << " W/m^2\n"
                << "  SW up (TOA) = " << SW_up_TOA << " W/m^2\n"
                << "  LW net (surface) = " << LW_net_surface << " W/m^2\n"
                << "  LW up (TOA) = " << LW_up_TOA << " W/m^2\n"
                << "  Max heating rate = " << heating_rate_max << " K/s\n";
        if (rad_choice.seb_diagnostic_enable && std::isfinite(seb_residual_mean)) {
            Print() << "  SEB residual (mean) = " << seb_residual_mean << " W/m^2\n"
                    << "  SEB residual (max) = " << seb_residual_max << " W/m^2\n";
        }
        if (rad_choice.seb_prognostic_enable && std::isfinite(t_s_mean)) {
            Print() << "  Surface temperature (mean) = " << t_s_mean << " K\n"
                    << "  Surface temperature (max) = " << t_s_max << " K\n"
                    << "  Surface moisture (mean) = " << q_s_mean << " m^3/m^3\n"
                    << "  Surface moisture (max) = " << q_s_max << " m^3/m^3\n";
        }
    }
 
    rad_diag.append(nstep, time, call_site, SW_surface, SW_TOA,
                    SW_up_TOA, LW_net_surface, LW_up_TOA, heating_rate_max,
                    seb_residual_mean, seb_residual_max,
                    t_s_mean, t_s_max, q_s_mean, q_s_max);
}

Here is the call graph for this function:

canopy_forcing()

TwoStreamCanopyView TwoStreamRadiation::canopy_forcing ( int  lev, int  nstep  ) const

inlinenoexcept

The beam and the cosine of the zenith of lev as the sweep of step nstep wrote them (the view's fluxes are left for the caller, who owns rad_fluxes), or an empty view when the level does not supply them or has not swept at that step: the faces must not take a previous step's radiation, nor the zeros of a level whose radiation came from its parent.

                                                                                         {
        TwoStreamCanopyView v;
        if (supplies_canopy_forcing(lev) && m_canopy_step[lev] == nstep) {
            v.beam = m_canopy[lev].get();
            v.cos_zenith = m_canopy_cosz[lev].get();
        }
        return v;
    }

Here is the call graph for this function:

define_level()

void TwoStreamRadiation::define_level	(	int	lev,
		const RadChoice &	rad_choice,
		amrex::Real	rdOcp,
		const amrex::BoxArray &	ba2d,
		const amrex::DistributionMapping &	dm,
		const amrex::BoxArray &	ba,
		const amrex::Box &	domain,
		bool	supply_land_forcing = false
	)

Allocate the 2D surface fields of one level on the given horizontal BoxArray and fill them with the scalar defaults. A no-op unless rad_choice.rad_type == TwoStream.

Also enforces the model's one grid requirement: every box of ba must span domain in z, because the column sweep applies the top-of-atmosphere and surface boundary conditions at the ends of a box. Aborting here means a bad refinement patch is rejected when the level is built (including on regrid), not at the first radiation call.

Parameters

lev	level
rad_choice	the model's inputs
ba2d	horizontal BoxArray the 2D surface fields live on
dm	distribution mapping shared with the level's state
ba	the level's 3D BoxArray, checked for the z span
domain	the level's domain Box
supply_land_forcing	allocate the downwelling SW, downwelling LW and cos(zenith) fields the sweep fills for a land-surface model (SolverChoice::rad_feeds_lsm()); off, the sweep does no extra work and the fields do not exist

{
    if (!rad_choice.enabled) { return; }
    m_rad = &rad_choice;
    m_rdOcp = rdOcp;
 
    // The sweep integrates a whole column in one pass inside a single box, applying the
    // top-of-atmosphere boundary condition above the top layer and the surface condition
    // below the bottom one. So the requirement is per box: every box must span the domain
    // in z. Levels that do not meet it come in several shapes -- a shallow nest, a level
    // tagged at different heights in different horizontal regions, grids decomposed in z --
    // and they all have the same consequence, that no box holds a column to sweep.
    //
    // Above level 0 that is not an error: ERF::advance_radiation interpolates the level's
    // heating rates and fluxes from its parent, the route RRTMGP takes for a nested patch,
    // and level_needs_interpolation() there applies this same per-box test so the two agree.
    //
    // On level 0 there is no parent to interpolate from, so it is fatal. Level 0 always
    // covers the domain, which makes a box that does not span z there exactly the
    // decomposed-in-z case -- and the remedies below are the ones that apply to it.
    if (lev == 0) {
        for (int ibox = 0; ibox < ba.size(); ++ibox) {
            const Box& b = ba[ibox];
            if (b.smallEnd(2) != domain.smallEnd(2) || b.bigEnd(2) != domain.bigEnd(2)) {
                amrex::Abort("erf.radiation_model = TwoStream cannot run on a level 0 whose "
                             "grids are decomposed in z: a box there holds only part of a "
                             "column, and level 0 has no parent to interpolate from. Set "
                             "amr.max_grid_size_z to at least the number of cells in z, or "
                             "leave amr.no_box_split_dir at its default of 2 so grids are "
                             "never split in z.");
            }
        }
    }
 
    // 2D surface fields on the horizontal BoxArray, one ghost cell in x and y
    const IntVect ng_sfc{1,1,0};
    m_alb_sw[lev]   = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
    m_emiss_lw[lev] = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
    m_t_sfc[lev]    = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
    m_alb_sw[lev]->setVal(rad_choice.surface_albedo_sw);
    m_emiss_lw[lev]->setVal(rad_choice.surface_emissivity_lw);
    m_t_sfc[lev]->setVal(rad_choice.rad_t_sfc);
 
    if (rad_choice.seb_enable) {
        m_sw_flux_sfc[lev] = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_lw_flux_sfc[lev] = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_hfx_sfc[lev]     = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_lh_sfc[lev]      = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_grdflx_sfc[lev]  = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_q_sfc[lev]       = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_t_deep[lev]      = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_q_deep[lev]      = std::make_unique<MultiFab>(ba2d, dm, 1, ng_sfc);
        m_sw_flux_sfc[lev]->setVal(rad_choice.seb_sw_flux_default);
        m_lw_flux_sfc[lev]->setVal(rad_choice.seb_lw_flux_default);
        m_hfx_sfc[lev]->setVal(rad_choice.seb_hfx_default);
        m_lh_sfc[lev]->setVal(rad_choice.seb_lh_default);
        m_grdflx_sfc[lev]->setVal(rad_choice.seb_grdflx_default);
        m_q_sfc[lev]->setVal(rad_choice.seb_q_sfc_default);
        m_t_deep[lev]->setVal(rad_choice.seb_t_deep_default);
        m_q_deep[lev]->setVal(rad_choice.seb_q_deep_default);
        if (rad_choice.seb_surface_layer_uses_moisture) {
            m_moisture[lev] = std::make_unique<MultiFab>(ba2d, dm,
                                                         erf_surface_moisture::NumComponents, ng_sfc);
            m_moisture[lev]->setVal(0.0);
        } else {
            m_moisture[lev].reset();
        }
    }
 
    // The land-model forcing (see the members). It holds the lsm_undefined sentinel until a
    // sweep on this level fills it. Zero would be a valid forcing (night, no sky), so a copy
    // made before any sweep -- a level whose advance() returned early, or a caller that
    // drifts out of step with the sweep -- would hand the land model a 0 K sky that its
    // first-land-step check (NOAHMP::Advance_With_State) accepts. The sentinel is what
    // that check looks for.
    if (supply_land_forcing) {
        m_sw_dn_sfc[lev]  = std::make_unique<MultiFab>(ba2d, dm, 1, 0);
        m_lw_dn_sfc[lev]  = std::make_unique<MultiFab>(ba2d, dm, 1, 0);
        m_cos_zenith[lev] = std::make_unique<MultiFab>(ba2d, dm, 1, 0);
        m_sw_dn_sfc[lev]->setVal(lsm_undefined);
        m_lw_dn_sfc[lev]->setVal(lsm_undefined);
        m_cos_zenith[lev]->setVal(lsm_undefined);
    } else {
        m_sw_dn_sfc[lev].reset();
        m_lw_dn_sfc[lev].reset();
        m_cos_zenith[lev].reset();
    }
    // A rebuilt level has new grids: the canopy forcing, if the faces asked for it,
    // must be asked for again on them (the faces do not regrid; ERF::init_ibseb()).
    m_canopy[lev].reset();
    m_canopy_cosz[lev].reset();
    m_canopy_m[lev] = -1;
    m_canopy_step[lev] = -1;
    m_flux_diag[lev] = FluxDiag{};
}

land_forcing_cos_zenith()

const amrex::MultiFab* TwoStreamRadiation::land_forcing_cos_zenith ( int  lev ) const

inlinenoexcept

                                                                                        {
        return supplies_land_forcing(lev) ? m_cos_zenith[lev].get() : nullptr;
    }

Here is the call graph for this function:

land_forcing_lw_dn()

const amrex::MultiFab* TwoStreamRadiation::land_forcing_lw_dn ( int  lev ) const

inlinenoexcept

                                                                                   {
        return supplies_land_forcing(lev) ? m_lw_dn_sfc[lev].get() : nullptr;
    }

Here is the call graph for this function:

land_forcing_sw_dn()

const amrex::MultiFab* TwoStreamRadiation::land_forcing_sw_dn ( int  lev ) const

inlinenoexcept

The surface forcing of the last sweep on lev (2D: i,j), or nullptr when the level does not supply it: total downwelling shortwave [W/m^2], downwelling longwave [W/m^2] and the cosine of the solar zenith angle floored at zero [1]. lsm_undefined until a sweep has run on the level, so a copy made before one is caught by the land model's missing-input check instead of reading as zero.

                                                                                   {
        return supplies_land_forcing(lev) ? m_sw_dn_sfc[lev].get() : nullptr;
    }

Here is the call graph for this function:

prognostic_surface_temperature_state()

const amrex::MultiFab* TwoStreamRadiation::prognostic_surface_temperature_state ( int  lev ) const

inlinenoexcept

Read-only view of the per-level prognostic surface-temperature state.

                                                                                                     {
        return (lev >= 0 && lev < static_cast<int>(m_t_sfc.size()))
            ? m_t_sfc[lev].get() : nullptr;
    }

read_checkpoint()

void TwoStreamRadiation::read_checkpoint	(	int	lev,
		const std::string &	restart_chkfile,
		const amrex::Periodicity &	period
	)

Read it back when the checkpoint carries it; else keep defaults.

{
    // Older checkpoints do not carry the state; then the defaults set by
    // define_level stand. Every level, matching write_checkpoint. A checkpoint
    // written before the SEB ran on fine levels simply has no file for them, which
    // the FileExists test below already handles.
    //
    // Each field is read on the grids in the file and redistributed onto the grids the
    // run uses, which differ when the restart re-made level 0 (issue 4225). The fields
    // are written with their ghost cells and an ordinary restart has always restored
    // those, so the halo is offered too.
    if (!active() || !m_rad->seb_enable) { return; }
    const std::string tsfc_name =
        MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "TwoStream_TSfc");
    if (m_t_sfc[lev] && amrex::FileExists(tsfc_name + "_H")) {
        erf_checkpoint::read_checkpoint_multifab(*m_t_sfc[lev], tsfc_name, 0, 0, 1,
                                                 m_t_sfc[lev]->nGrowVect(), true, period);
    }
    const std::string qsfc_name =
        MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "TwoStream_QSfc");
    if (m_q_sfc[lev] && amrex::FileExists(qsfc_name + "_H")) {
        erf_checkpoint::read_checkpoint_multifab(*m_q_sfc[lev], qsfc_name, 0, 0, 1,
                                                 m_q_sfc[lev]->nGrowVect(), true, period);
    }
}

Here is the call graph for this function:

resize()

void TwoStreamRadiation::resize

(

int

nlevs_max

)

Size the per-level containers; no allocation happens here.

{
    m_alb_sw.resize(nlevs_max);
    m_emiss_lw.resize(nlevs_max);
    m_t_sfc.resize(nlevs_max);
    m_sw_flux_sfc.resize(nlevs_max);
    m_lw_flux_sfc.resize(nlevs_max);
    m_hfx_sfc.resize(nlevs_max);
    m_lh_sfc.resize(nlevs_max);
    m_grdflx_sfc.resize(nlevs_max);
    m_sw_dn_sfc.resize(nlevs_max);
    m_lw_dn_sfc.resize(nlevs_max);
    m_cos_zenith.resize(nlevs_max);
    m_q_sfc.resize(nlevs_max);
    m_t_deep.resize(nlevs_max);
    m_q_deep.resize(nlevs_max);
    m_moisture.resize(nlevs_max);
    m_canopy.resize(nlevs_max);
    m_canopy_cosz.resize(nlevs_max);
    m_canopy_m.resize(nlevs_max, -1);
    m_canopy_step.resize(nlevs_max, -1);
    m_flux_diag.resize(nlevs_max);
    m_diag.resize(nlevs_max);
}

Referenced by ERF::ERF_shared().

Here is the caller graph for this function:

seb_hfx()

const amrex::MultiFab* TwoStreamRadiation::seb_hfx ( int  lev ) const

inline

The sensible and latent heat fluxes [W/m^2] the surface energy balance of one level last used, or nullptr when the SEB is off or the level has not been built. Read-only: for output.

                                                               {
        return (active() && m_rad->seb_enable && lev < m_hfx_sfc.size()) ? m_hfx_sfc[lev].get() : nullptr;
    }

Here is the call graph for this function:

seb_lh()

const amrex::MultiFab* TwoStreamRadiation::seb_lh ( int  lev ) const

inline

                                                              {
        return (active() && m_rad->seb_enable && lev < m_lh_sfc.size()) ? m_lh_sfc[lev].get() : nullptr;
    }

Here is the call graph for this function:

seb_prognostic_active()

bool TwoStreamRadiation::seb_prognostic_active ( ) const

inline

True when the force-restore state is being evolved and so needs transferring.

                                                      {
        return active() && m_rad->seb_enable && m_rad->seb_prognostic_enable;
    }

Here is the call graph for this function:

seb_q_sfc()

amrex::MultiFab* TwoStreamRadiation::seb_q_sfc ( int  lev )

inline

                                                     {
        return (active() && m_rad->seb_enable && lev < m_q_sfc.size()) ? m_q_sfc[lev].get() : nullptr;
    }

Here is the call graph for this function:

seb_surface_moisture()

const MultiFab * TwoStreamRadiation::seb_surface_moisture

(

int

lev

)

What level lev's surface hands the surface layer for its surface mixing ratio (erf.radiation.seb_surface_layer_uses_moisture), per column, in the components of erf_surface_moisture::Component: the soil-water factor (seb_moisture_availability of the current q_s); with erf.radiation.seb_soil_type, the bare-soil resistance and the relative humidity of its pore air at the skin temperature; and with erf.radiation.seb_vegetation_type as well, the vegetated fraction and the canopy resistance without its vapour-deficit factor (from the net shortwave the balance holds divided by (1 - surface albedo) – the sweep's with seb_use_radiation_fluxes, seb_sw_flux_default otherwise – and the skin temperature). Absent pieces are a fraction of 0 and resistances of 1e6 s/m. nullptr when the option is off or the level has not been built. Recomputed on every call, on the valid cells; the halo is not filled.

{
    if (!active() || !m_rad->seb_enable || !m_rad->seb_surface_layer_uses_moisture ||
        lev >= static_cast<int>(m_moisture.size()) || !m_moisture[lev] || !m_q_sfc[lev]) {
        return nullptr;
    }
    using namespace erf_surface_moisture;
    const Real q_wilt = m_rad->seb_soil_moisture_wilt;
    const Real q_fc = m_rad->seb_soil_moisture_fc;
 
    // With a soil type: the soil's parameters for the bare-soil resistance of the top
    // seb_moisture_layer_depth_m, vegetation or not. With a vegetation type as well (which
    // needs a soil type): Noah-MP's canopy parameters. RadChoice::init_params checked both
    // categories.
    const NoahMPVegetationParams* veg = (m_rad->seb_vegetation_type != 0)
                                      ? noahmp_vegetation_params(m_rad->seb_vegetation_type) : nullptr;
    const NoahMPSoilParams* soil = (m_rad->seb_soil_type != 0)
                                 ? noahmp_soil_params(m_rad->seb_soil_type) : nullptr;
    const bool has_soil = (soil != nullptr);
    const bool vegetated = (veg != nullptr && has_soil);
    const Real f_veg = vegetated ? m_rad->seb_vegetation_fraction : Real(0.0);
    const Real lai = m_rad->seb_leaf_area_index;
    const Real rs_min = vegetated ? veg->rs_min : Real(0.0);
    const Real rs_max = vegetated ? veg->rs_max : Real(0.0);
    const Real rgl = vegetated ? veg->rgl : Real(0.0);
    const Real t_opt = vegetated ? veg->t_opt : Real(0.0);
    const Real smc_max = has_soil ? soil->smc_max : Real(0.0);
    const Real bb = has_soil ? soil->bb : Real(0.0);
    const Real psi_sat = has_soil ? soil->psi_sat : Real(0.0);
    const Real d1 = m_rad->seb_moisture_layer_depth_m;
    const Real resistance_exponent = noahmp_soil_resistance_exponent;
 
    MultiFab& out = *m_moisture[lev];
    const MultiFab& q_s = *m_q_sfc[lev];
    const MultiFab& t_s = *m_t_sfc[lev];
    // The shortwave down at the surface, from the net the balance holds and the albedo.
    const MultiFab* sw_net = m_sw_flux_sfc[lev].get();
    const MultiFab* albedo = m_alb_sw[lev].get();
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(out, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        const auto out_arr = out.array(mfi);
        const auto q_arr = q_s.const_array(mfi);
        const auto t_arr = t_s.const_array(mfi);
        const auto sw_arr = (vegetated && sw_net) ? sw_net->const_array(mfi) : Array4<const Real>{};
        const auto alb_arr = (vegetated && albedo) ? albedo->const_array(mfi) : Array4<const Real>{};
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            const Real soil_factor = seb_moisture_availability(q_arr(i, j, k), q_wilt, q_fc);
            out_arr(i, j, k, SoilFactor) = soil_factor;
            out_arr(i, j, k, VegetationFraction) = f_veg;
            Real r_c = no_flux_resistance;
            Real r_soil = no_flux_resistance;
            if (vegetated) {
                const Real a = (alb_arr) ? alb_arr(i, j, k) : Real(0.0);
                const Real sw_down = (sw_arr && a < Real(1.0)) ? sw_arr(i, j, k) / (Real(1.0) - a) : Real(0.0);
                r_c = seb_canopy_resistance_without_vpd(rs_min, rs_max, rgl, t_opt, lai, sw_down,
                                                        t_arr(i, j, k), soil_factor);
            }
            Real rh_ground = Real(1.0);
            if (has_soil) {
                r_soil = seb_soil_evaporation_resistance(q_arr(i, j, k), smc_max, q_wilt, bb, d1,
                                                         resistance_exponent);
                rh_ground = seb_soil_surface_relative_humidity(q_arr(i, j, k), smc_max, psi_sat,
                                                               bb, t_arr(i, j, k));
            }
            out_arr(i, j, k, CanopyResistanceWithoutVPD) = r_c;
            out_arr(i, j, k, SoilResistance) = r_soil;
            out_arr(i, j, k, GroundRelativeHumidity) = rh_ground;
        });
    }
    return m_moisture[lev].get();
}

Here is the call graph for this function:

seb_t_sfc()

amrex::MultiFab* TwoStreamRadiation::seb_t_sfc ( int  lev )

inline

The prognostic surface state of one level, or nullptr when the SEB is off or the level has not been built.

Exposed so ERF can average t_sfc and q_sfc down after the finer levels have advanced. Every level evolves its own copy, and without that transfer the two estimates of one surface drift apart – which would give the longwave boundary condition two different values for the same ground. ERF owns refRatio and the level geometries; this class does not, and should not grow them for one call.

                                                     {
        return (active() && m_rad->seb_enable && lev < m_t_sfc.size()) ? m_t_sfc[lev].get() : nullptr;
    }

Here is the call graph for this function:

supplies_canopy_forcing()

bool TwoStreamRadiation::supplies_canopy_forcing ( int  lev ) const

inlinenoexcept

True when supply_canopy_forcing() was called for lev since its last define_level().

                                                                        {
        return lev >= 0 && lev < static_cast<int>(m_canopy.size()) && m_canopy[lev] != nullptr;
    }

Referenced by canopy_forcing().

Here is the caller graph for this function:

supplies_land_forcing()

bool TwoStreamRadiation::supplies_land_forcing ( int  lev ) const

inlinenoexcept

True when define_level allocated the land-model forcing fields on lev.

                                                                      {
        return lev >= 0 && lev < static_cast<int>(m_sw_dn_sfc.size()) && m_sw_dn_sfc[lev] != nullptr;
    }

Referenced by land_forcing_cos_zenith(), land_forcing_lw_dn(), and land_forcing_sw_dn().

Here is the caller graph for this function:

supply_canopy_forcing()

void TwoStreamRadiation::supply_canopy_forcing	(	int	lev,
		int	m_top
	)

Have every later sweep on lev write what the building faces need besides the interface fluxes (erf.ibseb.radiation = two_stream; ERF_TwoStreamCanopyForcing.H): the direct beam on a horizontal surface at the interfaces 0 .. m_top of each column, counted from the surface, and the cosine of the zenith. Allocates one field on the level's horizontal grids given those interfaces in z and one 2D field. Call after define_level(); a later define_level() of the level (a regrid) drops the request.

Parameters

lev	level, which must have been defined; its sweep asserts that the columns reach m_top
m_top	highest interface the faces read (IBFaceSet::top_sample_interface()), in [0, number of cells in z]

{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(active(),
        "TwoStreamRadiation::supply_canopy_forcing: erf.radiation_model is not TwoStream");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lev >= 0 && lev < static_cast<int>(m_alb_sw.size()) && m_alb_sw[lev],
        "TwoStreamRadiation::supply_canopy_forcing: level " + std::to_string(lev) + " has not been defined");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_top >= 0,
        "TwoStreamRadiation::supply_canopy_forcing: the top of the canopy must be at or above the surface");
    // The level's grids flattened (the 2D fields' boxes, on their distribution), given the
    // interfaces 0 .. m_top in z; the sweep asserts that the column reaches m_top.
    BoxList canopy_boxes = m_alb_sw[lev]->boxArray().boxList();
    for (Box& b : canopy_boxes) { b.setRange(2, 0, m_top + 1); }
    // The beam alone: the diffuse light, the longwave and the ground's fluxes are the
    // sweep's interface fluxes, which ERF keeps (rad_fluxes); only the beam and the
    // cosine of the zenith are lost after the sweep.
    m_canopy[lev] = std::make_unique<MultiFab>(BoxArray(std::move(canopy_boxes)), m_alb_sw[lev]->DistributionMap(), 1, 0);
    m_canopy_cosz[lev] = std::make_unique<MultiFab>(m_alb_sw[lev]->boxArray(), m_alb_sw[lev]->DistributionMap(), 1, 0);
    // Never read before a sweep has written them: canopy_forcing() returns null until then.
    m_canopy[lev]->setVal(0.0);
    m_canopy_cosz[lev]->setVal(0.0);
    m_canopy_m[lev] = m_top;
    m_canopy_step[lev] = -1;
}

Here is the call graph for this function:

write_checkpoint()

void TwoStreamRadiation::write_checkpoint	(	int	lev,
		const std::string &	checkpointname
	)		const

Write the prognostic surface state of one level (SEB on).

{
    // The force-restore state is the only part of this model a restart must
    // carry; without it T_s and q_s restart from the scalar defaults.
    //
    // Every level: each one evolves its own force-restore state now, so each has
    // state worth restoring. A coarse level's values where a finer level covers it
    // are the average of that finer level's, restored on both sides.
    if (!active() || !m_rad->seb_enable) { return; }
    if (m_t_sfc[lev]) {
        VisMF::Write(*m_t_sfc[lev],
                     MultiFabFileFullPrefix(lev, checkpointname, "Level_", "TwoStream_TSfc"));
    }
    if (m_q_sfc[lev]) {
        VisMF::Write(*m_q_sfc[lev],
                     MultiFabFileFullPrefix(lev, checkpointname, "Level_", "TwoStream_QSfc"));
    }
}

write_land_forcing()

void TwoStreamRadiation::write_land_forcing	(	int	lev,
		amrex::MultiFab *	sw_flux_dn,
		amrex::MultiFab *	lw_flux_dn,
		amrex::MultiFab *	cos_zenith_angle
	)		const

Copy the surface forcing of the last sweep on lev into a land-surface model's radiation inputs (SurfaceModel's "sw_flux_dn", "lw_flux_dn" and "cos_zenith_angle" outputs). A null destination is skipped. Each destination must lie on this level's horizontal grids; see copy_surface_plane. Before the level's first sweep this copies the lsm_undefined sentinel, not zero.

{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(supplies_land_forcing(lev),
        "TwoStreamRadiation::write_land_forcing: level " + std::to_string(lev) +
        " was not defined to supply a land-surface model's forcing");
    if (sw_flux_dn != nullptr)       { copy_surface_plane(*m_sw_dn_sfc[lev], *sw_flux_dn); }
    if (lw_flux_dn != nullptr)       { copy_surface_plane(*m_lw_dn_sfc[lev], *lw_flux_dn); }
    if (cos_zenith_angle != nullptr) { copy_surface_plane(*m_cos_zenith[lev], *cos_zenith_angle); }
}

Here is the call graph for this function:

Member Data Documentation

m_alb_sw

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_alb_sw

private

m_canopy

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_canopy

private

Referenced by canopy_forcing(), and supplies_canopy_forcing().

m_canopy_cosz

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_canopy_cosz

private

Referenced by canopy_forcing().

m_canopy_m

amrex::Vector<int> TwoStreamRadiation::m_canopy_m

private

m_canopy_step

amrex::Vector<int> TwoStreamRadiation::m_canopy_step

private

Referenced by canopy_forcing().

m_cos_zenith

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_cos_zenith

private

Referenced by land_forcing_cos_zenith().

m_diag

amrex::Vector<std::unique_ptr<RadiationDiagnostics> > TwoStreamRadiation::m_diag

private

m_emiss_lw

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_emiss_lw

private

m_flux_diag

amrex::Vector<FluxDiag> TwoStreamRadiation::m_flux_diag

private

m_grdflx_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_grdflx_sfc

private

m_hfx_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_hfx_sfc

private

Referenced by seb_hfx().

m_lh_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_lh_sfc

private

Referenced by seb_lh().

m_lw_dn_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_lw_dn_sfc

private

Referenced by land_forcing_lw_dn().

m_lw_flux_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_lw_flux_sfc

private

m_moisture

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_moisture

private

m_orbit

OrbitalCache TwoStreamRadiation::m_orbit

private

m_q_deep

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_q_deep

private

m_q_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_q_sfc

private

Referenced by seb_q_sfc().

m_rad

const RadChoice* TwoStreamRadiation::m_rad = nullptr

private

Referenced by active(), seb_hfx(), seb_lh(), seb_prognostic_active(), seb_q_sfc(), and seb_t_sfc().

m_rdOcp

amrex::Real TwoStreamRadiation::m_rdOcp = RdoCp

private

m_sw_dn_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_sw_dn_sfc

private

Referenced by land_forcing_sw_dn(), and supplies_land_forcing().

m_sw_flux_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_sw_flux_sfc

private

m_t_deep

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_t_deep

private

m_t_sfc

amrex::Vector<std::unique_ptr<amrex::MultiFab> > TwoStreamRadiation::m_t_sfc

private

Referenced by prognostic_surface_temperature_state(), and seb_t_sfc().

---

The documentation for this class was generated from the following files:

• Source/Radiation/TwoStream/ERF_TwoStreamRadiation.H
• Source/Radiation/TwoStream/ERF_TwoStreamRadiation.cpp