SurfaceModel Class Reference

Interface between ERF and land-surface and urban models. More...

#include <ERF_SurfaceModel.H>

Collaboration diagram for SurfaceModel:

Classes
struct	Field
	Describes a registered land/urban field mapping. More...
struct	RadiationField

Public Member Functions
	SurfaceModel (int nlevs, const amrex::Vector< amrex::BoxArray > &ba, const amrex::Vector< amrex::Geometry > &geom, const amrex::Vector< amrex::DistributionMapping > &dm, const SolverChoice &, amrex::Vector< amrex::Vector< std::unique_ptr< amrex::iMultiFab >>> &)
	Constructs a surface-model interface for the supplied AMR levels. More...
void	initialize_for_level (int lev, const amrex::BoxArray &ba, const amrex::Geometry &geom, const amrex::DistributionMapping &dm, const amrex::Vector< std::unique_ptr< amrex::iMultiFab >> &lmask_lev, const amrex::Vector< amrex::BCRec > &domain_bcs_type, const amrex::Vector< amrex::IntVect > &refRatio)
	Initializes storage and boundary data for an AMR level. More...
void	set_model_data (const int lev, const amrex::Vector< amrex::MultiFab * > model_data, const amrex::Vector< std::string > &data_names, SurfaceModelType type)
	Registers model variables used in weighted averaging. More...
void	set_model_fluxes (const int lev, const amrex::Vector< amrex::MultiFab * > model_fluxes, const amrex::Vector< std::string > &flux_names, SurfaceModelType type)
	Registers model fluxes as additional weighted-average inputs. More...
void	set_model_fields (SurfaceModelType type, const amrex::Vector< int > &field_indices, const bool use_fluxes=true)
	Selects the model fields that participate in weighted averaging. More...
void	apply_weight_average (int lev, const amrex::MultiFab *lsm_data, amrex::MultiFab *lsm_weighted, const amrex::MultiFab *urban_data, amrex::MultiFab *urban_weighted)
	Applies land and urban weights into separate output fields. More...
void	calculate_weight_average (int lev, amrex::MultiFab *const urban_frac, bool update_derived=true)
	Computes land and urban weights and applies them to registered fields. More...
bool	are_fluxes ()
	Returns whether the interface exports fluxes rather than MOST variables. More...
void	request_surface_outputs (const bool persistent=true)
	Enables and materializes the surface outputs consumed by ERF. More...
void	request_surface_layer_outputs ()
	Enables blended surface outputs needed by the SurfaceLayer. More...
void	mark_fields_valid ()
	Marks the registered surface fields as filled by an actual model advance. More...
bool	fields_are_valid () const
	Returns whether the surface models have advanced at least once. More...
amrex::MultiFab *	get_ustar (const int &lev)
	Returns the weighted horizontal momentum output for a level. More...
amrex::MultiFab *	get_tstar (const int &lev)
	Returns the weighted thermal output for a level. More...
amrex::MultiFab *	get_qstar (const int &lev)
	Returns the weighted moisture output for a level. More...
amrex::MultiFab *	get_tsurf (const int &lev)
	Returns the weighted surface-temperature output for a level. More...
SurfaceFluxView	get_surface_flux_view (const int lev, const amrex::MFIter &mfi) const
	Returns provider or blended surface-flux views for one tile. More...
const amrex::MultiFab *	get_weighted_model_data (const int lev, SurfaceModelType type, const int field_idx) const
	Returns weighted data for a selected model field. More...
amrex::MultiFab *	get_wavg_factors (const int &lev)
	Returns the land and urban weighting factors for a level. More...
void	activate_all_field_maps (const bool persistent=true)
	Activates all registered mapped fields for an output consumer. More...
void	write_output (const int &finest_lev, const amrex::Real &time, const std::string &plot_prefix, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &ref_ratio)
	Writes registered surface-model fields to plotfile output. More...
void	register_field_map (std::string name, const std::pair< int, int > &lsm_urb_map, bool fill_boundary=false)
	Registers a field map using land and urban field indices. More...
void	register_field_map (std::string name, amrex::Vector< amrex::MultiFab * > &lsm_lev_mf, amrex::Vector< amrex::MultiFab * > &urb_lev_mf, bool fill_boundary=false)
	Registers a field map using per-level land and urban pointers. More...
void	set_field_map_pointers (const std::string &name, int lev, amrex::MultiFab *lsm_mf, amrex::MultiFab *urban_mf)
	Updates a pointer-based field mapping for one level. More...
void	register_radiation_input (const std::string &name, const std::pair< int, int > &lsm_urb_map)
	Registers a canonical radiation input mapping. More...
void	register_radiation_inputs (const std::unordered_map< std::string, std::pair< int, int >> &input_map)
	Registers canonical radiation input mappings. More...
void	register_radiation_output (const std::string &name, const std::pair< int, int > &lsm_urb_map)
	Registers a canonical radiation output mapping. More...
void	register_radiation_outputs (const std::unordered_map< std::string, std::pair< int, int >> &output_map)
	Registers canonical radiation output mappings. More...
amrex::MultiFab *	get_field (const std::string &name, int lev=0)
	Retrieves a registered field by its common name. More...
const amrex::Vector< const amrex::MultiFab * >	get_radiation_fields (int lev=0)
	Returns fields registered under the RRTMGP radiation names. More...
const amrex::Vector< amrex::MultiFab * >	get_radiation_output_fields (int lev=0)
	Returns canonical radiation output destinations for one level. More...
amrex::MultiFab *	get_radiation_output_field (int lev, const std::string &name)
	Returns the destination of one canonical radiation output by name. More...
void	distribute_radiation_outputs (int lev)
	Distributes updated canonical radiation outputs to other providers. More...
void	distribute_radiation_output (int lev, int output_index)
	Distributes one updated canonical radiation output to other providers. More...
void	WriteCheckpoint (const std::string &checkpointname)
	Writes surface-model state to a checkpoint file. More...
void	ReadCheckpoint (const std::string &checkpointname)
	Restores surface-model state from a checkpoint file. More...
void	validate_radiation_output_layout (int lev, const amrex::MultiFab *mf) const
	Validates the layout contract for a canonical radiation output. More...
void	calculate_simple_average (int lev, amrex::MultiFab *const urban_frac)
	Computes simple land and urban weights from the urban fraction. More...
void	build_radiation_fieldlist (int lev)
	Builds the ordered radiation-field list for one AMR level. More...
void	weight_average_fields (int lev, amrex::MultiFab *const)
	Applies the current land and urban weights to selected model fields. More...
void	activate_field_map (const std::string &name, const bool persistent=true)
	Activates a mapped field and allocates its storage on all levels. More...
void	deactivate_transient_field_maps ()
	Disables and releases mapped fields used only by the last output. More...
bool	is_field_mapped (int lev, SurfaceModelType type, int field_idx, const amrex::MultiFab *mf) const
	Checks whether a model field participates in a registered mapping. More...
void	weight_model_field (int lev, const amrex::MultiFab *source, amrex::MultiFab *weighted, SurfaceModelType type)
	Writes a weighted copy of a model field. More...

Static Public Member Functions
static void	GotoNextLine (std::istream &is)
	Advances an input stream to the next line. More...

Private Member Functions
void	update_provider_mode (const int lev)
	Updates the provider mode after model data are registered. More...
void	ensure_weight_factors (const int lev)
	Creates weighting factors when an external caller requests them. More...
void	ensure_surface_outputs (const int lev)
	Allocates surface-output fields for one level. More...
void	release_transient_surface_outputs ()
	Releases surface outputs that have no persistent consumer. More...

Private Attributes
int	m_nlevs
bool	m_use_urban = false
bool	m_use_land = false
bool	m_export_fluxes = true
bool	m_export_fluxes_set = false
amrex::GpuArray< bool, 2 >	m_output_fields_registered {{false, false}}
bool	m_fields_are_valid = false
bool	m_weights_updated = false
amrex::Vector< SurfaceProviderMode >	m_provider_mode
amrex::Vector< amrex::BoxArray >	m_ba
amrex::Vector< amrex::BoxArray >	m_ba2d
amrex::Vector< amrex::Geometry >	m_geom
amrex::Vector< amrex::DistributionMapping >	m_dmap
amrex::Vector< amrex::Geometry >	m_geom2d
amrex::Vector< amrex::iMultiFab * >	m_lmask
amrex::Vector< amrex::MultiFab * >	urban_frac_lev
amrex::Vector< amrex::Vector< amrex::MultiFab * > >	lsm_data_lev
amrex::Vector< amrex::Vector< amrex::MultiFab * > >	urban_data_lev
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > >	weighted_lsm_data_lev
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > >	weighted_urban_data_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	u_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	t_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	q_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	t_surf
amrex::Vector< int >	m_surface_outputs_requested
bool	m_surface_outputs_enabled = false
bool	m_surface_layer_outputs_enabled = false
amrex::Vector< amrex::MultiFab * >	m_last_urban_frac
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	wavg
amrex::Vector< int >	lsm_fields
amrex::Vector< int >	urban_fields
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > >	fields
std::unordered_map< std::string, Field >	fieldmap
amrex::Vector< std::string >	m_lsm_names
amrex::Vector< std::string >	m_urban_names
amrex::Vector< amrex::Vector< const amrex::MultiFab * > >	rad_fields
amrex::Vector< amrex::Vector< amrex::MultiFab * > >	rad_output_fields
std::unordered_map< std::string, RadiationField >	radiation_input_map
std::unordered_map< std::string, RadiationField >	radiation_output_map

Static Private Attributes
static const std::vector< std::string >	radnames = {"tskin", "emiss", "albedo_vis", "albedo_nir", "albedo_vis_diff", "albedo_nir_diff"}
static const std::vector< std::string >	rad_output_names
static const std::vector< std::string >	field_names

Detailed Description

Interface between ERF and land-surface and urban models.

SurfaceModel collects data and fluxes from the LSM and urban models, performs the land-urban weight averaging and exposes the resulting fields for use by ERF boundary-condition and physics components.

Models may register data and mappings independently, without requiring a fixed ordering or one-to-one correspondence between their fields.

Model-specific field names can also be registered under a common name for consumers such as RRTMGP and MOST.

This class supports its own checkpoint state, ensuring the mappings are restored correctly on restart. This class also provides a surface plotfile for weight-averaged model data for the fields registered.

Constructor & Destructor Documentation

SurfaceModel()

SurfaceModel::SurfaceModel ( int  nlevs, const amrex::Vector< amrex::BoxArray > &  ba, const amrex::Vector< amrex::Geometry > &  geom, const amrex::Vector< amrex::DistributionMapping > &  dm, const SolverChoice &  , amrex::Vector< amrex::Vector< std::unique_ptr< amrex::iMultiFab >>> &    )

inlineexplicit

Constructs a surface-model interface for the supplied AMR levels.

Parameters

nlevs	Number of AMR levels initially available.
ba	Cell-centered box arrays for each level.
geom	Geometry for each level.
dm	Distribution mappings for each level.
solverChoice	ERF solver configuration.
lmask_lev	Land-mask data for each level.

        : m_nlevs(nlevs), m_ba(ba), m_geom(geom), m_dmap(dm) {
 
        u_star.resize(m_nlevs);
        t_star.resize(m_nlevs);
        q_star.resize(m_nlevs);
        t_surf.resize(m_nlevs);
 
        wavg.resize(m_nlevs);
        m_provider_mode.resize(m_nlevs, SurfaceProviderMode::None);
        m_output_fields_registered = {false, false};
 
        m_ba2d.resize(m_nlevs);
        m_geom2d.resize(m_nlevs);
        m_lmask.resize(m_nlevs);
 
 
        urban_frac_lev.resize(m_nlevs);
        lsm_data_lev.resize(m_nlevs);
        urban_data_lev.resize(m_nlevs);
        weighted_lsm_data_lev.resize(m_nlevs);
        weighted_urban_data_lev.resize(m_nlevs);
 
        rad_fields.resize(m_nlevs);
        rad_output_fields.resize(m_nlevs);
        m_surface_outputs_requested.resize(m_nlevs, false);
        m_last_urban_frac.resize(m_nlevs, nullptr);
 
        // set initial fields for each model to none
        lsm_fields = amrex::Vector<int>(5, -1);
        urban_fields = amrex::Vector<int>(5, -1);
    }

Member Function Documentation

activate_all_field_maps()

void SurfaceModel::activate_all_field_maps

(

const bool

persistent = true

)

Activates all registered mapped fields for an output consumer.

{
    for (const auto& entry : fieldmap) {
        activate_field_map(entry.first, persistent);
    }
}

activate_field_map()

void SurfaceModel::activate_field_map	(	const std::string &	name,
		const bool	persistent = true
	)

Activates a mapped field and allocates its storage on all levels.

Parameters

name	Common mapped field name.

{
    auto field_it = fieldmap.find(name);
    if (field_it == fieldmap.end()) { return; }
 
    Field& field = field_it->second;
    field.active = true;
    field.persistent_consumer = field.persistent_consumer || persistent;
    if (field.mf_ind == -1) {
        field.mf_ind = static_cast<int>(fields.size());
        fields.emplace_back(m_nlevs);
    }
 
    for (int lev = 0; lev < m_nlevs; ++lev) {
        if (fields[field.mf_ind][lev] == nullptr) {
            fields[field.mf_ind][lev] = std::make_unique<amrex::MultiFab>(
                m_ba2d[lev], m_dmap[lev], 1, IntVect(1,1,0));
            fields[field.mf_ind][lev]->setVal(0.0);
        }
    }
}

Referenced by get_field(), and initialize_for_level().

Here is the caller graph for this function:

apply_weight_average()

void SurfaceModel::apply_weight_average	(	int	lev,
		const amrex::MultiFab *	lsm_data,
		amrex::MultiFab *	lsm_weighted,
		const amrex::MultiFab *	urban_data,
		amrex::MultiFab *	urban_weighted
	)

Applies land and urban weights into separate output fields.

A null land or urban source is skipped, but both sources cannot be null. The source fields are not modified and the fields may use different underlying grids.

Parameters

lev	AMR level whose weights are used.
lsm_data	Land-model source data to weight.
lsm_weighted	Destination for weighted land data.
urban_data	Urban-model source data to weight.
urban_weighted	Destination for weighted urban data.

{
    bool valid_land = (lsm_data != nullptr);
    bool valid_urban = (urban_data != nullptr);
 
    AMREX_ASSERT_WITH_MESSAGE(valid_land || valid_urban, "Need at least one pointer to apply weights");
    AMREX_ASSERT_WITH_MESSAGE(!valid_land || lsm_weighted != nullptr,
                              "Need a destination for weighted land data");
    AMREX_ASSERT_WITH_MESSAGE(!valid_urban || urban_weighted != nullptr,
                              "Need a destination for weighted urban data");
    AMREX_ASSERT_WITH_MESSAGE(lsm_data == nullptr || lsm_data != lsm_weighted,
                              "Weighted land data must use separate storage");
    AMREX_ASSERT_WITH_MESSAGE(urban_data == nullptr || urban_data != urban_weighted,
                              "Weighted urban data must use separate storage");
 
    // make sure the weights have been calculated before applying them
    AMREX_ALWAYS_ASSERT(m_weights_updated);
 
    // Apply weights separately to support different underlying grids between lsm_data and urban_data
    if (valid_land) {
        weight_model_field(lev, lsm_data, lsm_weighted, SurfaceModelType::LAND);
    }
 
    if (valid_urban) {
        weight_model_field(lev, urban_data, urban_weighted, SurfaceModelType::URBAN);
    }
}

Here is the call graph for this function:

are_fluxes()

bool SurfaceModel::are_fluxes ( )

inline

Returns whether the interface exports fluxes rather than MOST variables.

Returns

True when weighted outputs are fluxes.

{ return m_export_fluxes; }

build_radiation_fieldlist()

void SurfaceModel::build_radiation_fieldlist ( int  lev )

inline

Builds the ordered radiation-field list for one AMR level.

Parameters

lev	AMR level to update.

                                             {
        amrex::Print() << "  --- building radiation field list" << std::endl;
 
        // loop over fields and check against names of radiation variables
        for (int li = 0; li < static_cast<int>(radnames.size()); ++li) {
            amrex::Print() << "    - Checking for radiation variable '" << radnames[li] << "'" << std::endl;
 
            amrex::MultiFab* mf = nullptr;
            auto rad_it = radiation_input_map.find(radnames[li]);
            if (rad_it != radiation_input_map.end()) {
                const int land_idx = rad_it->second.map.first;
                const int urban_idx = rad_it->second.map.second;
                const SurfaceProviderMode mode = m_provider_mode[lev];
                const bool use_land = mode == SurfaceProviderMode::LandOnly ||
                                      mode == SurfaceProviderMode::Both;
                const bool use_urban = mode == SurfaceProviderMode::UrbanOnly ||
                                       mode == SurfaceProviderMode::Both;
                const bool valid_land = use_land && land_idx >= 0 &&
                    lev < static_cast<int>(lsm_data_lev.size()) &&
                    land_idx < static_cast<int>(lsm_data_lev[lev].size()) &&
                    lsm_data_lev[lev][land_idx] != nullptr;
                const bool valid_urban = use_urban && urban_idx >= 0 &&
                    lev < static_cast<int>(urban_data_lev.size()) &&
                    urban_idx < static_cast<int>(urban_data_lev[lev].size()) &&
                    urban_data_lev[lev][urban_idx] != nullptr;
                if (valid_land && valid_urban) {
                    if (!rad_it->second.weighted[lev]) {
                        rad_it->second.weighted[lev] = std::make_unique<amrex::MultiFab>(
                            m_ba2d[lev], m_dmap[lev], 1, amrex::IntVect(1,1,0));
                    }
                    auto& weighted = *rad_it->second.weighted[lev];
                    const auto& land = *lsm_data_lev[lev][land_idx];
                    const auto& urban = *urban_data_lev[lev][urban_idx];
                    for (amrex::MFIter mfi(weighted); mfi.isValid(); ++mfi) {
                        const amrex::Box bx = mfi.tilebox();
                        const auto weights = wavg[lev]->const_array(mfi);
                        const auto land_arr = land.const_array(mfi);
                        const auto urban_arr = urban.const_array(mfi);
                        // LSM surface fields may be stored at k <= 0; use the top valid LSM
                        // plane, which matches k=0 when the LSM keeps those fields synchronized
                        // at the surface (such as SLM). Urban surface fields are stored at k=0.
                        const int land_k = land.box(mfi.index()).bigEnd(2);
                        const int urban_k = urban.box(mfi.index()).smallEnd(2);
                        auto output = weighted.array(mfi);
                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                            output(i,j,k) = land_arr(i,j,land_k) *
                                weights(i,j,0,SurfaceModelType::LAND) +
                                urban_arr(i,j,urban_k) *
                                weights(i,j,0,SurfaceModelType::URBAN);
                        });
                    }
                    weighted.FillBoundary(m_geom[lev].periodicity());
                    mf = rad_it->second.weighted[lev].get();
                } else if (valid_land) {
                    mf = lsm_data_lev[lev][land_idx];
                } else if (valid_urban) {
                    mf = urban_data_lev[lev][urban_idx];
                }
            }
            if (mf) {
                rad_fields[lev].push_back(mf);
            } else {
                amrex::Print() << "     - WARNING: did not find radiation variable, setting to null to use default RRTMGP value!" << std::endl;
                rad_fields[lev].push_back(nullptr);
            }
        }
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            static_cast<int>(rad_fields[lev].size()) == static_cast<int>(radnames.size()),
            "Did not find all radiation fields");
    }

Referenced by get_radiation_fields().

Here is the call graph for this function:

Here is the caller graph for this function:

calculate_simple_average()

void SurfaceModel::calculate_simple_average	(	int	lev,
		amrex::MultiFab *const	urban_frac
	)

Computes simple land and urban weights from the urban fraction.

Parameters

lev	AMR level to update.
urban_frac	Urban fraction used to compute the weights.

{
    const SurfaceProviderMode mode = m_provider_mode[lev];
 
    // If no urban fraction is available, use the only enabled model as the
    // complete surface.  A missing fraction is ambiguous when both models
    // are enabled, so fail explicitly rather than silently choosing a model.
    if (urban_frac == nullptr) {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            mode != SurfaceProviderMode::Both,
            "Urban fraction is required when both land and urban models are enabled");
    }
 
    if (mode != SurfaceProviderMode::Both) {
        m_weights_updated = true;
        return;
    }
 
    for (MFIter mfi(*wavg[lev], TileNoZ()); mfi.isValid(); ++mfi)
    {
        Box tbx = mfi.tilebox();
        // weights are defined only on k = 0
        tbx = tbx.makeSlab(2, 0);
 
        auto urban_frac_arr = (urban_frac != nullptr) ? urban_frac->const_array(mfi) : Array4<const Real>{};
        auto weights_arr = wavg[lev]->array(mfi);
 
        ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
        {
            // Weights are proportional to urban fraction coverage in the current cell.
            weights_arr(i, j, 0, SurfaceModelType::URBAN) = urban_frac_arr(i, j, 0);
            weights_arr(i, j, 0, SurfaceModelType::LAND) = 1.0 - urban_frac_arr(i, j, 0);
        });
    }
 
    m_weights_updated = true;
}

Here is the call graph for this function:

calculate_weight_average()

void SurfaceModel::calculate_weight_average	(	int	lev,
		amrex::MultiFab *const	urban_frac,
		bool	update_derived = true
	)

Computes land and urban weights and applies them to registered fields.

Parameters

lev	AMR level to update.
urban_frac	Urban fraction used to compute the weights.

{
    m_last_urban_frac[lev] = urban_frac;
    const SurfaceProviderMode mode = m_provider_mode[lev];
    // An urban model can be enabled on only a subset of AMR levels.  Those
    // inactive levels have no surface-provider data to update.
    if (mode == SurfaceProviderMode::None) {
        return;
    }
 
    // Calculate weighted averages between Land and Urban.  Output consumers
    // can reuse the weights computed during the timestep update.
    if (update_derived || !m_weights_updated) {
        calculate_simple_average(lev, urban_frac);
    }
 
    const bool use_urban = mode == SurfaceProviderMode::UrbanOnly ||
                           mode == SurfaceProviderMode::Both;
    const bool use_land = mode == SurfaceProviderMode::LandOnly ||
                          mode == SurfaceProviderMode::Both;
    const bool use_weights = mode == SurfaceProviderMode::Both;
    const bool use_fluxes = m_export_fluxes;
    const bool output_fields_registered =
        (use_land && m_output_fields_registered[SurfaceModelType::LAND]) ||
        (use_urban && m_output_fields_registered[SurfaceModelType::URBAN]);
 
    amrex::MultiFab* const outputs[] = {u_star[lev].get(), t_star[lev].get(), q_star[lev].get(), t_surf[lev].get()};
 
    const int nfields = (m_export_fluxes) ? 5 : 4;
 
    if (m_surface_outputs_requested[lev] && output_fields_registered && !use_weights) {
        const auto& provider_fields = use_land ? lsm_fields : urban_fields;
        const auto& provider_data = use_land ? lsm_data_lev[lev] : urban_data_lev[lev];
        auto source_for = [&] (const int field) -> const amrex::MultiFab* {
            if (field >= static_cast<int>(provider_fields.size()) ||
                provider_fields[field] < 0 ||
                provider_fields[field] >= static_cast<int>(provider_data.size())) {
                return nullptr;
            }
            return provider_data[provider_fields[field]];
        };
 
        const amrex::MultiFab* source0 = source_for(0);
        const amrex::MultiFab* source1 = source_for(1);
        const amrex::MultiFab* source2 = source_for(2);
        const amrex::MultiFab* source3 = source_for(3);
        const amrex::MultiFab* source4 = source_for(4);
 
        for (MFIter mfi(*u_star[lev], TileNoZ()); mfi.isValid(); ++mfi)
        {
            const Box tbx = mfi.tilebox();
            const auto source0_arr = source0 ? source0->const_array(mfi) : Array4<const Real>{};
            const auto source1_arr = source1 ? source1->const_array(mfi) : Array4<const Real>{};
            const auto source2_arr = source2 ? source2->const_array(mfi) : Array4<const Real>{};
            const auto source3_arr = source3 ? source3->const_array(mfi) : Array4<const Real>{};
            const auto source4_arr = source4 ? source4->const_array(mfi) : Array4<const Real>{};
            const int source0_k = source0 ? (use_land ? source0->box(mfi.index()).bigEnd(2)
                                                      : source0->box(mfi.index()).smallEnd(2)) : 0;
            const int source1_k = source1 ? (use_land ? source1->box(mfi.index()).bigEnd(2)
                                                      : source1->box(mfi.index()).smallEnd(2)) : 0;
            const int source2_k = source2 ? (use_land ? source2->box(mfi.index()).bigEnd(2)
                                                      : source2->box(mfi.index()).smallEnd(2)) : 0;
            const int source3_k = source3 ? (use_land ? source3->box(mfi.index()).bigEnd(2)
                                                      : source3->box(mfi.index()).smallEnd(2)) : 0;
            const int source4_k = source4 ? (use_land ? source4->box(mfi.index()).bigEnd(2)
                                                      : source4->box(mfi.index()).smallEnd(2)) : 0;
 
            auto u_arr = u_star[lev]->array(mfi);
            auto t_arr = t_star[lev]->array(mfi);
            auto q_arr = q_star[lev]->array(mfi);
            auto tsurf_arr = t_surf[lev]->array(mfi);
            ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                if (use_fluxes) {
                    u_arr(i, j, k, 0) = source0_arr ? source0_arr(i, j, source0_k) : 0.0;
                    u_arr(i, j, k, 1) = source1_arr ? source1_arr(i, j, source1_k) : 0.0;
                    t_arr(i, j, k) = source2_arr ? source2_arr(i, j, source2_k) : 0.0;
                    q_arr(i, j, k) = source3_arr ? source3_arr(i, j, source3_k) : 0.0;
                    tsurf_arr(i, j, k) = source4_arr ? source4_arr(i, j, source4_k) : 0.0;
                } else {
                    u_arr(i, j, k, 0) = source0_arr ? source0_arr(i, j, source0_k) : 0.0;
                    u_arr(i, j, k, 1) = 0.0;
                    t_arr(i, j, k) = source1_arr ? source1_arr(i, j, source1_k) : 0.0;
                    q_arr(i, j, k) = source2_arr ? source2_arr(i, j, source2_k) : 0.0;
                    tsurf_arr(i, j, k) = source3_arr ? source3_arr(i, j, source3_k) : 0.0;
                }
            });
        }
 
        u_star[lev]->FillBoundary(m_geom[lev].periodicity());
        t_star[lev]->FillBoundary(m_geom[lev].periodicity());
        q_star[lev]->FillBoundary(m_geom[lev].periodicity());
        t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
    }
 
    // Output weighted surface fluxes into ustar, tstar, qstar and surface temperature into tsurf
    // TODO: make sure grids of urban and LSM inputs match
    if (m_surface_outputs_requested[lev] && output_fields_registered && use_weights) {
    if (!use_fluxes) {
        u_star[lev]->setVal(0.0, 1, 1, 0);
    }
    for (int field = 0; field < nfields; ++field)
    {
        int output_field = (use_fluxes && field > 0) ? field - 1 : field;
        int comp = (use_fluxes && field < 2) ? field : 0;
 
        // whether we have a valid LSM multifab
        bool valid_land = (use_land &&
                           lsm_fields[field] != -1 &&
                           lsm_data_lev[lev][lsm_fields[field]]);
 
        // whether we have a valid urban multifab
        bool valid_urban = (use_urban &&
                            urban_fields[field] != -1 &&
                            urban_data_lev[lev][urban_fields[field]]);
 
        for (MFIter mfi(*u_star[lev], TileNoZ()); mfi.isValid(); ++mfi)
        {
            Box tbx = mfi.tilebox();
 
            const auto weights_arr = use_weights
                ? wavg[lev]->const_array(mfi) : Array4<const Real>{};
 
            // Outputs for surface boundary condition
            auto output_arr = outputs[output_field]->array(mfi);
            const amrex::MultiFab* lsm_mf = valid_land ? lsm_data_lev[lev][lsm_fields[field]] : nullptr;
            const amrex::MultiFab* urban_mf = valid_urban ? urban_data_lev[lev][urban_fields[field]] : nullptr;
            const int lsm_khi = valid_land ? lsm_mf->box(mfi.index()).bigEnd(2) : 0;
            const int urban_klo = valid_urban ? urban_mf->box(mfi.index()).smallEnd(2) : 0;
            auto lsm_data_arr = valid_land ? lsm_mf->const_array(mfi) : Array4<const Real>{};
            auto urban_data_arr = valid_urban ? urban_mf->const_array(mfi) : Array4<const Real>{};
 
            ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                if (valid_land && use_land && !use_urban) {
                    output_arr(i, j, k, comp) = lsm_data_arr(i, j, lsm_khi);
                } else if (valid_urban && use_urban && !use_land) {
                    output_arr(i, j, k, comp) = urban_data_arr(i, j, urban_klo);
                } else {
                    Real land = (lsm_data_arr) ? lsm_data_arr(i, j, lsm_khi) * weights_arr(i, j, 0, SurfaceModelType::LAND) : 0.0;
                    Real urb  = (urban_data_arr) ? urban_data_arr(i, j, urban_klo) * weights_arr(i, j, 0, SurfaceModelType::URBAN) : 0.0;
                    output_arr(i, j, k, comp) = land + urb;
                }
 
            });
        }
 
        outputs[output_field]->FillBoundary(comp, 1, m_geom[lev].periodicity());
    }
    }
 
    // Write weighted copies only when both providers contribute.  With one
    // provider, get_weighted_model_data() returns the provider-owned source.
    if (update_derived && use_weights && use_land) {
        for (int field = nfields; field < static_cast<int>(lsm_fields.size()); ++field) {
            if (lsm_fields[field] == -1) continue;
            if (is_field_mapped(lev, SurfaceModelType::LAND, lsm_fields[field],
                                lsm_data_lev[lev][lsm_fields[field]])) continue;
            const int field_idx = lsm_fields[field];
            if (static_cast<int>(weighted_lsm_data_lev[lev].size()) <= field_idx) {
                weighted_lsm_data_lev[lev].resize(lsm_data_lev[lev].size());
            }
            if (weighted_lsm_data_lev[lev][field_idx] == nullptr) {
                const auto* source = lsm_data_lev[lev][field_idx];
                weighted_lsm_data_lev[lev][field_idx] = std::make_unique<MultiFab>(
                    source->boxArray(), source->DistributionMap(), source->nComp(),
                    source->nGrowVect());
            }
            weight_model_field(lev, lsm_data_lev[lev][field_idx],
                               weighted_lsm_data_lev[lev][field_idx].get(),
                               SurfaceModelType::LAND);
        }
    }
 
    if (update_derived && use_weights && use_urban) {
        for (int field = nfields; field < static_cast<int>(urban_fields.size()); ++field) {
            if (urban_fields[field] == -1) continue;
            if (is_field_mapped(lev, SurfaceModelType::URBAN, urban_fields[field],
                                urban_data_lev[lev][urban_fields[field]])) continue;
            const int field_idx = urban_fields[field];
            if (static_cast<int>(weighted_urban_data_lev[lev].size()) <= field_idx) {
                weighted_urban_data_lev[lev].resize(urban_data_lev[lev].size());
            }
            if (weighted_urban_data_lev[lev][field_idx] == nullptr) {
                const auto* source = urban_data_lev[lev][field_idx];
                weighted_urban_data_lev[lev][field_idx] = std::make_unique<MultiFab>(
                    source->boxArray(), source->DistributionMap(), source->nComp(),
                    source->nGrowVect());
            }
            weight_model_field(lev, urban_data_lev[lev][field_idx],
                               weighted_urban_data_lev[lev][field_idx].get(),
                               SurfaceModelType::URBAN);
        }
    }
 
    weight_average_fields(lev, urban_frac);
 
    if (!update_derived) { return; }
 
    for (auto& entry : radiation_input_map) {
        const int land_idx = entry.second.map.first;
        const int urban_idx = entry.second.map.second;
        const bool valid_land = use_land && land_idx >= 0 &&
            land_idx < static_cast<int>(lsm_data_lev[lev].size()) &&
            lsm_data_lev[lev][land_idx] != nullptr;
        const bool valid_urban = use_urban && urban_idx >= 0 &&
            urban_idx < static_cast<int>(urban_data_lev[lev].size()) &&
            urban_data_lev[lev][urban_idx] != nullptr;
        if (!(valid_land && valid_urban)) { continue; }
 
        if (entry.second.weighted[lev] == nullptr) {
            entry.second.weighted[lev] = std::make_unique<MultiFab>(
                m_ba2d[lev], m_dmap[lev], 1, IntVect(1,1,0));
        }
        MultiFab& output = *entry.second.weighted[lev];
        const MultiFab& land = *lsm_data_lev[lev][land_idx];
        const MultiFab& urban = *urban_data_lev[lev][urban_idx];
        for (MFIter mfi(output, TileNoZ()); mfi.isValid(); ++mfi) {
            const Box bx = mfi.tilebox();
            const auto weights = wavg[lev]->const_array(mfi);
            const auto land_arr = land.const_array(mfi);
            const auto urban_arr = urban.const_array(mfi);
            // LSM surface fields may be stored at k <= 0; use the top valid LSM
            // plane, which matches k=0 when the LSM keeps those fields synchronized
            // at the surface (such as SLM). Urban surface fields are stored at k=0.
            const int land_k = land.box(mfi.index()).bigEnd(2);
            const int urban_k = urban.box(mfi.index()).smallEnd(2);
            auto out = output.array(mfi);
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                out(i,j,k) = land_arr(i,j,land_k) * weights(i,j,0,SurfaceModelType::LAND) +
                             urban_arr(i,j,urban_k) * weights(i,j,0,SurfaceModelType::URBAN);
            });
        }
        output.FillBoundary(m_geom[lev].periodicity());
    }
}

Here is the call graph for this function:

deactivate_transient_field_maps()

void SurfaceModel::deactivate_transient_field_maps

(

)

Disables and releases mapped fields used only by the last output.

{
    for (auto& entry : fieldmap) {
        Field& field = entry.second;
        if (field.persistent_consumer) { continue; }
        field.active = false;
        if (field.mf_ind == -1) { continue; }
        for (auto& mf : fields[field.mf_ind]) {
            mf.reset();
        }
    }
}

distribute_radiation_output()

void SurfaceModel::distribute_radiation_output	(	int	lev,
		int	output_index
	)

Distributes one updated canonical radiation output to other providers.

Parameters

lev	AMR level to update.
output_index	Index in the canonical radiation output list.

{
    AMREX_ALWAYS_ASSERT(output_index >= 0 &&
                        output_index < static_cast<int>(rad_output_names.size()));
    auto it = radiation_output_map.find(rad_output_names[output_index]);
    if (it == radiation_output_map.end()) { return; }
 
    const int land_idx = it->second.map.first;
    const int urban_idx = it->second.map.second;
    const SurfaceProviderMode mode = m_provider_mode[lev];
    const bool use_land = mode == SurfaceProviderMode::LandOnly ||
                          mode == SurfaceProviderMode::Both;
    const bool use_urban = mode == SurfaceProviderMode::UrbanOnly ||
                           mode == SurfaceProviderMode::Both;
    MultiFab* lsm = (use_land && land_idx >= 0 && land_idx < static_cast<int>(lsm_data_lev[lev].size()))
        ? lsm_data_lev[lev][land_idx] : nullptr;
    MultiFab* urban = (use_urban && urban_idx >= 0 && urban_idx < static_cast<int>(urban_data_lev[lev].size()))
        ? urban_data_lev[lev][urban_idx] : nullptr;
 
    // The LSM destination is always primary when it is available.  Urban
    // receives the updated surface plane only when both destinations exist.
    if (lsm && urban && lsm != urban) {
        for (MFIter mfi(*urban, TileNoZ()); mfi.isValid(); ++mfi) {
            const Box& source_box = (*lsm)[mfi].box();
            const Box& target_box = (*urban)[mfi].box();
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                source_box.smallEnd(2) <= 0 && source_box.bigEnd(2) >= 0 &&
                target_box.smallEnd(2) <= 0 && target_box.bigEnd(2) >= 0,
                "Radiation output destinations must contain the k=0 surface plane");
            //
            // These are the *grown* boxes, so the two destinations need not carry the same
            // ghost vector.  Loop over the intersection of the two surface planes rather than
            // over the target alone: indexing the source at the target's indices would read
            // out of bounds wherever the urban halo reaches past the LSM one.
            //
            const Box copy_slab = makeSlab(source_box, 2, 0) & makeSlab(target_box, 2, 0);
            if (copy_slab.isEmpty()) { continue; }
            const auto source = (*lsm)[mfi].array();
            auto target = (*urban)[mfi].array();
            amrex::ParallelFor(copy_slab, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                target(i, j, k) = source(i, j, k);
            });
        }
    }
}

Here is the call graph for this function:

distribute_radiation_outputs()

void SurfaceModel::distribute_radiation_outputs

(

int

lev

)

Distributes updated canonical radiation outputs to other providers.

Parameters

lev	AMR level to update.

{
    for (int i = 0; i < static_cast<int>(rad_output_names.size()); ++i) {
        distribute_radiation_output(lev, i);
    }
}

ensure_surface_outputs()

void SurfaceModel::ensure_surface_outputs ( const int  lev )

inlineprivate

Allocates surface-output fields for one level.

Parameters

lev	AMR level whose outputs should be materialized.

    {
        if (u_star[lev] != nullptr) { return; }
 
        const amrex::IntVect ng(1,1,0);
        u_star[lev] = std::make_unique<amrex::MultiFab>(m_ba2d[lev], m_dmap[lev], 2, ng);
        t_star[lev] = std::make_unique<amrex::MultiFab>(m_ba2d[lev], m_dmap[lev], 1, ng);
        q_star[lev] = std::make_unique<amrex::MultiFab>(m_ba2d[lev], m_dmap[lev], 1, ng);
        t_surf[lev] = std::make_unique<amrex::MultiFab>(m_ba2d[lev], m_dmap[lev], 1, ng);
        u_star[lev]->setVal(0.0);
        t_star[lev]->setVal(0.0);
        q_star[lev]->setVal(0.0);
        t_surf[lev]->setVal(300.0);
    }

Referenced by initialize_for_level(), request_surface_layer_outputs(), request_surface_outputs(), and update_provider_mode().

Here is the caller graph for this function:

ensure_weight_factors()

void SurfaceModel::ensure_weight_factors ( const int  lev )

inlineprivate

Creates weighting factors when an external caller requests them.

Parameters

lev	AMR level whose weighting factors should be materialized.

    {
        if (wavg[lev] != nullptr) { return; }
 
        const amrex::IntVect ng(1,1,0);
        wavg[lev] = std::make_unique<amrex::MultiFab>(m_ba2d[lev], m_dmap[lev], 2, ng);
        const SurfaceProviderMode mode = m_provider_mode[lev];
        wavg[lev]->setVal(mode != SurfaceProviderMode::UrbanOnly ? 1.0 : 0.0,
                          SurfaceModelType::LAND, 1, 0);
        wavg[lev]->setVal(mode == SurfaceProviderMode::UrbanOnly ? 1.0 : 0.0,
                          SurfaceModelType::URBAN, 1, 0);
    }

Referenced by get_wavg_factors(), and update_provider_mode().

Here is the caller graph for this function:

fields_are_valid()

bool SurfaceModel::fields_are_valid ( ) const

inline

Returns whether the surface models have advanced at least once.

Returns

False until the first advance, while the registered fields still hold their initial values rather than anything the models computed.

{ return m_fields_are_valid; }

get_field()

amrex::MultiFab* SurfaceModel::get_field ( const std::string &  name, int  lev = 0  )

inline

Retrieves a registered field by its common name.

Parameters

name	Registered field name.
lev	AMR level to query.

Returns

Pointer to the registered field, or nullptr if it is unknown.

                                                                  {
        if (auto field = fieldmap.find(name); field != fieldmap.end()) {
            activate_field_map(name);
            return fields[field->second.mf_ind][lev].get();
        } else {
            return nullptr;
        }
    }

Here is the call graph for this function:

get_qstar()

amrex::MultiFab* SurfaceModel::get_qstar ( const int &  lev )

inline

Returns the weighted moisture output for a level.

Parameters

lev	AMR level to query.

Returns

Pointer to the weighted q-star or vapor-flux field.

                                              {
        return q_star[lev].get();
    }

get_radiation_fields()

const amrex::Vector<const amrex::MultiFab*> SurfaceModel::get_radiation_fields ( int  lev = 0 )

inline

Returns fields registered under the RRTMGP radiation names.

The resulting map is cached after the first call for reuse on subsequent timesteps. When an expected field is missing, it is marked nullptr for the RRTMGP implementation to apply its fallback values.

Parameters

lev	AMR level to query.

Returns

Radiation fields in the order expected by RRTMGP.

                                                                               {
        // check if we need to build the radiation list (for first time step,)
        if (rad_fields[lev].size() == 0) {
            build_radiation_fieldlist(lev);
        }
        return rad_fields[lev];
    }

Here is the call graph for this function:

get_radiation_output_field()

MultiFab * SurfaceModel::get_radiation_output_field	(	int	lev,
		const std::string &	name
	)

Returns the destination of one canonical radiation output by name.

Parameters

lev	AMR level to query.
name	Canonical radiation output name (aborts on an unknown name).

Returns

The registered destination, or nullptr when none is registered.

{
    const auto it = std::find(rad_output_names.begin(), rad_output_names.end(), name);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != rad_output_names.end(),
                                     "Unknown canonical radiation output: " + name);
    return get_radiation_output_fields(lev)[static_cast<int>(it - rad_output_names.begin())];
}

Here is the call graph for this function:

get_radiation_output_fields()

const Vector< MultiFab * > SurfaceModel::get_radiation_output_fields

(

int

lev = 0

)

Returns canonical radiation output destinations for one level.

Parameters

lev	AMR level to query.

Returns

One destination per canonical radiation output.

{
    if (!rad_output_fields[lev].empty()) { return rad_output_fields[lev]; }
    rad_output_fields[lev].resize(rad_output_names.size(), nullptr);
    const SurfaceProviderMode mode = m_provider_mode[lev];
    const bool use_land = mode == SurfaceProviderMode::LandOnly ||
                          mode == SurfaceProviderMode::Both;
    const bool use_urban = mode == SurfaceProviderMode::UrbanOnly ||
                           mode == SurfaceProviderMode::Both;
    for (int i = 0; i < static_cast<int>(rad_output_names.size()); ++i) {
        auto it = radiation_output_map.find(rad_output_names[i]);
        if (it == radiation_output_map.end()) { continue; }
        const int land_idx = it->second.map.first;
        const int urban_idx = it->second.map.second;
        const bool valid_land = use_land && land_idx >= 0 &&
            land_idx < static_cast<int>(lsm_data_lev[lev].size()) &&
            lsm_data_lev[lev][land_idx];
        const bool valid_urban = use_urban && urban_idx >= 0 &&
            urban_idx < static_cast<int>(urban_data_lev[lev].size()) &&
            urban_data_lev[lev][urban_idx];
        if (valid_land) {
            validate_radiation_output_layout(lev, lsm_data_lev[lev][land_idx]);
            rad_output_fields[lev][i] = lsm_data_lev[lev][land_idx];
        } else if (valid_urban) {
            rad_output_fields[lev][i] = urban_data_lev[lev][urban_idx];
        }
        if (valid_urban) {
            validate_radiation_output_layout(lev, urban_data_lev[lev][urban_idx]);
        }
    }
    return rad_output_fields[lev];
}

get_surface_flux_view()

SurfaceFluxView SurfaceModel::get_surface_flux_view ( const int  lev, const amrex::MFIter &  mfi  ) const

inline

Returns provider or blended surface-flux views for one tile.

Parameters

lev	AMR level to query.
mfi	Tile whose flux arrays are requested.

Returns

Read-only flux views selected for the level's provider mode.

    {
        AMREX_ALWAYS_ASSERT(m_export_fluxes);
 
        const SurfaceProviderMode mode = m_provider_mode[lev];
        if (mode == SurfaceProviderMode::Both) {
            AMREX_ALWAYS_ASSERT(u_star[lev] != nullptr);
            AMREX_ALWAYS_ASSERT(t_star[lev] != nullptr);
            AMREX_ALWAYS_ASSERT(q_star[lev] != nullptr);
            return {u_star[lev]->const_array(mfi, 0),
                    u_star[lev]->const_array(mfi, 1),
                    t_star[lev]->const_array(mfi),
                    q_star[lev]->const_array(mfi)};
        }
 
        const auto& provider_fields = (mode == SurfaceProviderMode::LandOnly)
            ? lsm_fields : urban_fields;
        const auto& provider_data = (mode == SurfaceProviderMode::LandOnly)
            ? lsm_data_lev[lev] : urban_data_lev[lev];
        auto get_provider_field = [&] (const int field) -> amrex::Array4<const amrex::Real> {
            if (field >= static_cast<int>(provider_fields.size())) {
                return {};
            }
            const int index = provider_fields[field];
            if (index < 0 || index >= static_cast<int>(provider_data.size()) ||
                provider_data[index] == nullptr) {
                return {};
            }
            return provider_data[index]->const_array(mfi);
        };
 
        return {get_provider_field(0), get_provider_field(1),
                get_provider_field(2), get_provider_field(3)};
    }

Here is the call graph for this function:

get_tstar()

amrex::MultiFab* SurfaceModel::get_tstar ( const int &  lev )

inline

Returns the weighted thermal output for a level.

Parameters

lev	AMR level to query.

Returns

Pointer to the weighted t-star or heat-flux field.

                                              {
        return t_star[lev].get();
    }

get_tsurf()

amrex::MultiFab* SurfaceModel::get_tsurf ( const int &  lev )

inline

Returns the weighted surface-temperature output for a level.

Parameters

lev	AMR level to query.

Returns

Pointer to the weighted surface-temperature field.

                                              {
        return t_surf[lev].get();
    }

get_ustar()

amrex::MultiFab* SurfaceModel::get_ustar ( const int &  lev )

inline

Returns the weighted horizontal momentum output for a level.

Parameters

lev	AMR level to query.

Returns

Pointer to the weighted u-star or momentum-flux field.

                                              {
        return u_star[lev].get();
    }

get_wavg_factors()

amrex::MultiFab* SurfaceModel::get_wavg_factors ( const int &  lev )

inline

Returns the land and urban weighting factors for a level.

Parameters

lev	AMR level to query.

Returns

Pointer to the weighting-factor field.

                                                     {
        ensure_weight_factors(lev);
        return wavg[lev].get();
    }

Here is the call graph for this function:

get_weighted_model_data()

const amrex::MultiFab* SurfaceModel::get_weighted_model_data ( const int  lev, SurfaceModelType  type, const int  field_idx  ) const

inline

Returns weighted data for a selected model field.

Parameters

lev	AMR level to query.
type	Provider type owning the field.
field_idx	Provider field index.

Returns

Weighted field or the source field for a single provider, or nullptr if it has not been selected.

                                                                               {
        const SurfaceProviderMode mode = m_provider_mode[lev];
        const bool single_provider =
            (mode == SurfaceProviderMode::LandOnly && type == SurfaceModelType::LAND) ||
            (mode == SurfaceProviderMode::UrbanOnly && type == SurfaceModelType::URBAN);
        if (single_provider) {
            const auto& selected = (type == SurfaceModelType::LAND)
                ? lsm_fields : urban_fields;
            const int first_output_field = m_export_fluxes ? 5 : 4;
            bool is_weighted_field = false;
            for (int field = first_output_field;
                 field < static_cast<int>(selected.size()); ++field) {
                if (selected[field] == field_idx) {
                    is_weighted_field = true;
                    break;
                }
            }
            if (!is_weighted_field) {
                return nullptr;
            }
 
            const auto& source = (type == SurfaceModelType::LAND)
                ? lsm_data_lev[lev] : urban_data_lev[lev];
            if (field_idx < 0 || field_idx >= static_cast<int>(source.size())) {
                return nullptr;
            }
            return source[field_idx];
        }
 
        const auto& weighted = (type == SurfaceModelType::LAND)
            ? weighted_lsm_data_lev[lev] : weighted_urban_data_lev[lev];
        if (field_idx < 0 || field_idx >= static_cast<int>(weighted.size())) {
            return nullptr;
        }
        return weighted[field_idx].get();
    }

GotoNextLine()

void SurfaceModel::GotoNextLine ( std::istream &  is )

static

Advances an input stream to the next line.

taken from ERF_Checkpoint, reuse?

Parameters

is	Input stream to advance.

{
    constexpr std::streamsize bl_ignore_max { 100000 };
    is.ignore(bl_ignore_max, '\n');
}

initialize_for_level()

void SurfaceModel::initialize_for_level ( int  lev, const amrex::BoxArray &  ba, const amrex::Geometry &  geom, const amrex::DistributionMapping &  dm, const amrex::Vector< std::unique_ptr< amrex::iMultiFab >> &  lmask_lev, const amrex::Vector< amrex::BCRec > &  domain_bcs_type, const amrex::Vector< amrex::IntVect > &  refRatio  )

inline

Initializes storage and boundary data for an AMR level.

Existing coarser-level data are interpolated when a finer level is initialized.

Parameters

lev	AMR level to initialize.
ba	Cell-centered box array for the level.
geom	Geometry for the level.
dm	Distribution mapping for the level.
lmask_lev	Land-mask data used by the level.
domain_bcs_type	Boundary conditions for interpolation.
refRatio	Refinement ratios between adjacent levels.

    {
        if (lev >= m_nlevs) {
            m_nlevs = lev+1;
 
            u_star.resize(m_nlevs);
            t_star.resize(m_nlevs);
            q_star.resize(m_nlevs);
            t_surf.resize(m_nlevs);
 
            wavg.resize(m_nlevs);
            m_provider_mode.resize(m_nlevs, SurfaceProviderMode::None);
 
            m_dmap.resize(m_nlevs);
            m_ba.resize(m_nlevs);
            m_ba2d.resize(m_nlevs);
            m_geom.resize(m_nlevs);
            m_geom2d.resize(m_nlevs);
            m_lmask.resize(m_nlevs);
 
            urban_frac_lev.resize(m_nlevs);
            lsm_data_lev.resize(m_nlevs);
            urban_data_lev.resize(m_nlevs);
            weighted_lsm_data_lev.resize(m_nlevs);
            weighted_urban_data_lev.resize(m_nlevs);
 
            rad_fields.resize(m_nlevs);
            rad_output_fields.resize(m_nlevs);
            m_surface_outputs_requested.resize(m_nlevs, m_surface_outputs_enabled);
            m_last_urban_frac.resize(m_nlevs, nullptr);
        }
 
        m_ba[lev] = ba;
        m_geom[lev] = geom;
        m_dmap[lev] = dm;
 
        // Create a 2D ba, dm, & ghost cells
        amrex::BoxList bl2d = ba.boxList();
        for (auto& b : bl2d) {
            b.setRange(2, 0);
        }
        m_ba2d[lev] = amrex::BoxArray(std::move(bl2d));
 
        amrex::RealBox dom2d = geom.ProbDomain();
        dom2d.setHi(2, geom.CellSize(2));
        m_geom2d[lev].define(makeSlab(m_geom[lev].Domain(), 2, 0), dom2d, geom.Coord(), geom.isPeriodic());
 
        m_lmask[lev] = lmask_lev[0].get();
 
        weighted_lsm_data_lev[lev].clear();
        weighted_urban_data_lev[lev].clear();
 
        u_star[lev].reset();
        t_star[lev].reset();
        q_star[lev].reset();
        t_surf[lev].reset();
        wavg[lev].reset();
        if (m_surface_outputs_requested[lev]) {
            ensure_surface_outputs(lev);
        }
        m_provider_mode[lev] = SurfaceProviderMode::None;
 
        // Recreate registered fields on the current level layout. This is also
        // needed when an existing AMR level is remade with a new grid layout.
        for (int field = 0; field < static_cast<int>(fields.size()); ++field) {
            if (static_cast<int>(fields[field].size()) < m_nlevs) {
                fields[field].resize(m_nlevs);
            }
            fields[field][lev].reset();
        }
        for (const auto& entry : fieldmap) {
            if (entry.second.active) {
                activate_field_map(entry.first);
            }
        }
        for (auto& entry : radiation_input_map) {
            if (static_cast<int>(entry.second.weighted.size()) < m_nlevs) {
                entry.second.weighted.resize(m_nlevs);
            }
            entry.second.weighted[lev].reset();
        }
 
        // Cached radiation fields point into the registered storage above.
        // Rebuild the cache when this level's fields have been replaced.
        rad_fields[lev].clear();
        rad_output_fields[lev].clear();
 
        // Interpolate the new finer level using coarse data
        if (lev > 0) {
            amrex::Real time_for_fp = zero; // This is not actually used
            amrex::Vector<amrex::Real> ftime    = {time_for_fp, time_for_fp};
            amrex::Vector<amrex::Real> ctime    = {time_for_fp, time_for_fp};
            amrex::IntVect ng_od(0, 0, 0); // ng outside domain
            amrex::Interpolater* mapper = &amrex::cell_cons_interp;
            amrex::Vector<amrex::MultiFab*> fmf;
            amrex::Vector<amrex::MultiFab*> cmf;
            if (u_star[lev] != nullptr && u_star[lev-1] != nullptr) {
            amrex::InterpFromCoarseLevel(*u_star[lev], u_star[lev-1]->nGrowVect(), ng_od,
                                         *u_star[lev-1], 0, 0, 2,
                                         m_geom[lev-1], m_geom[lev],
                                         refRatio[lev-1], &amrex::cell_cons_interp,
                                         domain_bcs_type, BCVars::cons_bc);
            amrex::InterpFromCoarseLevel(*t_star[lev], t_star[lev-1]->nGrowVect(), ng_od,
                                         *t_star[lev-1], 0, 0, 1,
                                         m_geom[lev-1], m_geom[lev],
                                         refRatio[lev-1], &amrex::cell_cons_interp,
                                         domain_bcs_type, BCVars::cons_bc);
            amrex::InterpFromCoarseLevel(*q_star[lev], q_star[lev-1]->nGrowVect(), ng_od,
                                         *q_star[lev-1], 0, 0, 1,
                                         m_geom[lev-1], m_geom[lev],
                                         refRatio[lev-1], &amrex::cell_cons_interp,
                                         domain_bcs_type, BCVars::cons_bc);
            amrex::InterpFromCoarseLevel(*t_surf[lev], t_surf[lev-1]->nGrowVect(), ng_od,
                                         *t_surf[lev-1], 0, 0, 1,
                                         m_geom[lev-1], m_geom[lev],
                                         refRatio[lev-1], &amrex::cell_cons_interp,
                                         domain_bcs_type, BCVars::cons_bc);
 
            fmf = {u_star[lev  ].get(), u_star[lev  ].get()};
            cmf = {u_star[lev-1].get(), u_star[lev-1].get()};
            amrex::FillPatchTwoLevels(*u_star[lev].get(), u_star[lev]->nGrowVect(), ng_od,
                                      time_for_fp, cmf, ctime, fmf, ftime,
                                      0, 0, 2, m_geom[lev-1], m_geom[lev],
                                      refRatio[lev-1], mapper, domain_bcs_type,
                                      BCVars::cons_bc);
            fmf = {t_star[lev  ].get(), t_star[lev  ].get()};
            cmf = {t_star[lev-1].get(), t_star[lev-1].get()};
            amrex::FillPatchTwoLevels(*t_star[lev].get(), t_star[lev]->nGrowVect(), ng_od,
                                      time_for_fp, cmf, ctime, fmf, ftime,
                                      0, 0, 1, m_geom[lev-1], m_geom[lev],
                                      refRatio[lev-1], mapper, domain_bcs_type,
                                      BCVars::cons_bc);
            fmf = {q_star[lev  ].get(), q_star[lev  ].get()};
            cmf = {q_star[lev-1].get(), q_star[lev-1].get()};
            amrex::FillPatchTwoLevels(*q_star[lev].get(), q_star[lev]->nGrowVect(), ng_od,
                                      time_for_fp, cmf, ctime, fmf, ftime,
                                      0, 0, 1, m_geom[lev-1], m_geom[lev],
                                      refRatio[lev-1], mapper, domain_bcs_type,
                                      BCVars::cons_bc);
            fmf = {t_surf[lev  ].get(), t_surf[lev  ].get()};
            cmf = {t_surf[lev-1].get(), t_surf[lev-1].get()};
            amrex::FillPatchTwoLevels(*t_surf[lev].get(), t_surf[lev]->nGrowVect(), ng_od,
                                      time_for_fp, cmf, ctime, fmf, ftime,
                                      0, 0, 1, m_geom[lev-1], m_geom[lev],
                                      refRatio[lev-1], mapper, domain_bcs_type,
                                      BCVars::cons_bc);
            }
 
            for (int field = 0; field < static_cast<int>(fields.size()); ++field) {
                if (fields[field][lev] == nullptr || fields[field][lev-1] == nullptr) { continue; }
                amrex::IntVect ngv = fields[field][lev-1]->nGrowVect(); ngv[2] = 0;
                amrex::InterpFromCoarseLevel(*fields[field][lev], ngv, ng_od,
                                             *fields[field][lev-1], 0, 0, 1,
                                             m_geom[lev-1], m_geom[lev],
                                             refRatio[lev-1], &amrex::cell_cons_interp,
                                             domain_bcs_type, BCVars::cons_bc);
                fmf = {fields[field][lev  ].get(), fields[field][lev  ].get()};
                cmf = {fields[field][lev-1].get(), fields[field][lev-1].get()};
                amrex::FillPatchTwoLevels(*fields[field][lev].get(), ngv, ng_od,
                                          time_for_fp, cmf, ctime, fmf, ftime,
                                          0, 0, 1, m_geom[lev-1], m_geom[lev],
                                          refRatio[lev-1], mapper, domain_bcs_type,
                                          BCVars::cons_bc);
            }
        }
    }

Here is the call graph for this function:

is_field_mapped()

bool SurfaceModel::is_field_mapped	(	int	lev,
		SurfaceModelType	type,
		int	field_idx,
		const amrex::MultiFab *	mf
	)		const

Checks whether a model field participates in a registered mapping.

Parameters

lev	AMR level to query.
type	Land or urban model owning the field.
field_idx	Index of the model field.
mf	Field pointer to test.

Returns

True when the field is part of a registered mapping.

{
    for (const auto& entry : fieldmap) {
        const Field& field = entry.second;
 
        const int mapped_idx = (type == SurfaceModelType::LAND)
            ? field.map.first : field.map.second;
        if (mapped_idx == field_idx) {
            return true;
        }
 
        // Pointer-based mappings have no model-field index.  Match the
        // registered pointer for the current level instead.
        if (field.map.first == -1 && field.map.second == -1 && mf != nullptr) {
            const auto& ptrs = (type == SurfaceModelType::LAND)
                ? field.lsm_ptr : field.urb_ptr;
            if (lev < static_cast<int>(ptrs.size()) && ptrs[lev] == mf) {
                return true;
            }
        }
    }
 
    for (const auto& entry : radiation_input_map) {
        const auto& map = entry.second.map;
        const int mapped_idx = (type == SurfaceModelType::LAND)
            ? map.first : map.second;
        if (mapped_idx == field_idx) {
            return true;
        }
    }
 
    return false;
}

Here is the call graph for this function:

mark_fields_valid()

void SurfaceModel::mark_fields_valid ( )

inline

Marks the registered surface fields as filled by an actual model advance.

{ m_fields_are_valid = true; }

ReadCheckpoint()

void SurfaceModel::ReadCheckpoint

(

const std::string &

checkpointname

)

Restores surface-model state from a checkpoint file.

Parameters

checkpointname

Checkpoint file name.

{
    request_surface_outputs(false);
    activate_all_field_maps(false);
    auto check_start = amrex::second();
 
    amrex::Print() << " Reading SurfaceModel checkpoint " << std::endl;
 
    auto checkpoint_error = [&] (const std::string& message) {
        amrex::Abort("SurfaceModel checkpoint '" + checkpointname + "': " + message);
    };
 
    // Header
    const std::string File(checkpointname + "/SurfaceModel_Header");
    if (!amrex::FileExists(File)) {
        checkpoint_error("missing SurfaceModel_Header");
    }
 
    Vector<char> fileCharPtr;
    ParallelDescriptor::ReadAndBcastFile(File, fileCharPtr);
    std::string fileCharPtrString(fileCharPtr.dataPtr(), fileCharPtr.size());
    std::istringstream is(fileCharPtrString, std::istringstream::in);
 
    std::string line;
 
    // read in title line
    if (!std::getline(is, line) || line != "Checkpoint file for SurfaceModel") {
        checkpoint_error("invalid or truncated header title");
    }
 
    // read in number of levels
    int chk_nlevs = -1;
    if (!(is >> chk_nlevs) || chk_nlevs < 1) {
        checkpoint_error("invalid number of levels");
    }
 
    // read in flags
    bool chk_use_urban = false;
    bool chk_use_land = false;
    bool chk_export_fluxes = false;
    bool chk_weights_updated = false;
    bool chk_fields_are_valid = false;
    if (!(is >> chk_use_urban >> chk_use_land >> chk_export_fluxes >> chk_weights_updated
             >> chk_fields_are_valid)) {
        checkpoint_error("invalid flags in header");
    }
    GotoNextLine(is);
 
    if (chk_nlevs != m_nlevs) {
        checkpoint_error("level count mismatch: checkpoint has " + std::to_string(chk_nlevs) +
                         ", current SurfaceModel has " + std::to_string(m_nlevs));
    }
    if (chk_use_urban != m_use_urban || chk_use_land != m_use_land ||
        chk_export_fluxes != m_export_fluxes) {
        checkpoint_error("model configuration flags do not match the current SurfaceModel");
    }
 
    // Read box arrays into temporaries so current runtime geometry is not
    // overwritten before it has been validated.
    amrex::Vector<amrex::BoxArray> chk_ba(chk_nlevs);
    amrex::Vector<amrex::BoxArray> chk_ba2d(chk_nlevs);
    for (int lev = 0; lev < chk_nlevs; ++lev) {
        chk_ba[lev].readFrom(is);
        if (!is) checkpoint_error("invalid 3D BoxArray for level " + std::to_string(lev));
        GotoNextLine(is);
    }
    GotoNextLine(is);
    for (int lev = 0; lev < chk_nlevs; ++lev) {
        chk_ba2d[lev].readFrom(is);
        if (!is) checkpoint_error("invalid 2D BoxArray for level " + std::to_string(lev));
        GotoNextLine(is);
    }
    GotoNextLine(is);
 
    // The checkpointed fields are read back on the decomposition they were written with and
    // then redistributed onto the live one, so the two BoxArrays need not match box for box.
    // They do have to cover the same index space: ERF turns regridding of level 0 on by itself
    // when a restart uses more ranks than level 0 has boxes (see ERF::restart), which changes
    // the boxes but not the region they tile.
    const std::string coverage_hint =
        "; the SurfaceModel checkpoint can be read back on a different decomposition, but not "
        "on a different domain -- check that the restart uses the same grid extents and "
        "refinement as the run that wrote the checkpoint";
    amrex::Vector<amrex::DistributionMapping> chk_dmap(chk_nlevs);
    for (int lev = 0; lev < chk_nlevs; ++lev) {
        if (chk_ba[lev].minimalBox() != m_ba[lev].minimalBox()) {
            checkpoint_error("3D domain coverage mismatch at level " + std::to_string(lev) +
                             coverage_hint);
        }
        if (chk_ba2d[lev].minimalBox() != m_ba2d[lev].minimalBox()) {
            checkpoint_error("2D domain coverage mismatch at level " + std::to_string(lev) +
                             coverage_hint);
        }
        chk_dmap[lev].define(chk_ba2d[lev]);
    }
 
    // Read in LSM fields
    if (!std::getline(is, line)) {
        checkpoint_error("missing LSM field list");
    }
    amrex::Vector<int> chk_lsm_fields;
    std::istringstream lsm_stream(line);
    int field_idx = -1;
    while (lsm_stream >> field_idx) chk_lsm_fields.push_back(field_idx);
 
    // Read in Urban fields
    if (!std::getline(is, line)) {
        checkpoint_error("missing urban field list");
    }
    amrex::Vector<int> chk_urban_fields;
    std::istringstream urban_stream(line);
    while (urban_stream >> field_idx) chk_urban_fields.push_back(field_idx);
    GotoNextLine(is);
 
    if (chk_lsm_fields != lsm_fields || chk_urban_fields != urban_fields) {
        checkpoint_error("model field lists do not match the current SurfaceModel");
    }
 
    // Read number of mapped fields
    int nfields = 0;
    if (!(is >> nfields) || nfields < 0) {
        checkpoint_error("invalid mapped-field count");
    }
    GotoNextLine(is);
 
    // Read any mapped fields
    if (nfields != static_cast<int>(fieldmap.size())) {
        checkpoint_error("mapped-field count does not match the current SurfaceModel");
    }
 
    // The field mapping should already be created before the restart, but verify consistency with the file
    amrex::Vector<std::string> checkpoint_field_names;
    for (int i = 0; i < nfields; i++) {
        if (!std::getline(is, line)) {
            checkpoint_error("missing mapped-field entry " + std::to_string(i));
        }
        std::istringstream lis(line);
 
        std::string field_name;
        int chk_mf_ind = -1;
        int lsm_ind = -1;
        int urb_ind = -1;
        int fill_bound_int = -1;
        if (!(lis >> field_name >> chk_mf_ind >> lsm_ind >> urb_ind >> fill_bound_int) ||
            chk_mf_ind < 0 || lsm_ind < -1 || urb_ind < -1 ||
            (fill_bound_int != 0 && fill_bound_int != 1)) {
            checkpoint_error("invalid mapped-field entry for '" + field_name + "'");
        }
        const bool fill_bound = (fill_bound_int != 0);
 
        auto field_it = fieldmap.find(field_name);
        if (field_it == fieldmap.end()) {
            checkpoint_error("checkpoint mapping '" + field_name + "' is not registered");
        }
        for (const auto& seen_name : checkpoint_field_names) {
            if (seen_name == field_name) {
                checkpoint_error("duplicate mapped-field entry for '" + field_name + "'");
            }
        }
        checkpoint_field_names.push_back(field_name);
        const auto &field = field_it->second;
        if (field.mf_ind < 0 || field.mf_ind >= static_cast<int>(fields.size())) {
            checkpoint_error("current mapped-field storage is invalid for '" + field_name + "'");
        }
        if (static_cast<int>(fields[field.mf_ind].size()) != m_nlevs) {
            checkpoint_error("current mapped-field level storage is invalid for '" + field_name + "'");
        }
        if (field.map.first != lsm_ind || field.map.second != urb_ind ||
            field.fill_bound != fill_bound) {
            checkpoint_error("mapped-field metadata mismatch for '" + field_name + "'");
        }
        if (lsm_ind == -1 && urb_ind == -1) {
            if (static_cast<int>(field.lsm_ptr.size()) != m_nlevs ||
                static_cast<int>(field.urb_ptr.size()) != m_nlevs) {
                checkpoint_error("pointer mapping has the wrong number of levels for '" + field_name + "'");
            }
            for (int lev = 0; lev < m_nlevs; ++lev) {
                if (fields[field.mf_ind][lev] == nullptr) {
                    checkpoint_error("mapped-field storage is missing at level " +
                                     std::to_string(lev) + " for '" + field_name + "'");
                }
            }
        } else {
            for (int lev = 0; lev < m_nlevs; ++lev) {
                const bool invalid_land_index =
                    !lsm_data_lev[lev].empty() &&
                    lsm_ind >= static_cast<int>(lsm_data_lev[lev].size());
                const bool invalid_urban_index =
                    !urban_data_lev[lev].empty() &&
                    urb_ind >= static_cast<int>(urban_data_lev[lev].size());
                if (invalid_land_index || invalid_urban_index) {
                    checkpoint_error("mapped-field index is out of range for '" + field_name + "'");
                }
                if (fields[field.mf_ind][lev] == nullptr) {
                    checkpoint_error("mapped-field storage is missing at level " +
                                     std::to_string(lev) + " for '" + field_name + "'");
                }
            }
        }
        amrex::ignore_unused(chk_mf_ind);
    }
 
    const std::string prefix = "SurfaceModel_";
 
    IntVect ng = IntVect(1,1,0);
    auto validate_multifab_header = [&] (const std::string& name, int lev, int ncomp) {
        const std::string header_name = name + "_H";
        if (!amrex::FileExists(header_name)) {
            checkpoint_error("missing MultiFab header '" + header_name + "'");
        }
 
        Vector<char> header_data;
        ParallelDescriptor::ReadAndBcastFile(header_name, header_data);
        std::istringstream header_stream(
            std::string(header_data.dataPtr(), header_data.size()), std::istringstream::in);
        VisMF::Header header;
        if (!(header_stream >> header)) {
            checkpoint_error("invalid MultiFab header '" + header_name + "'");
        }
        if (header.m_ncomp != ncomp || header.m_ngrow != ng ||
            !(header.m_ba == chk_ba2d[lev])) {
            checkpoint_error("MultiFab layout mismatch for level " + std::to_string(lev) +
                             " in '" + name + "'");
        }
        if (static_cast<amrex::Long>(header.m_fod.size()) != chk_ba2d[lev].size()) {
            checkpoint_error("MultiFab FAB count mismatch for level " + std::to_string(lev) +
                             " in '" + name + "'");
        }
        const std::string data_dir = VisMF::DirName(name);
        for (const auto& fab_file : header.m_fod) {
            if (!amrex::FileExists(data_dir + fab_file.m_name)) {
                checkpoint_error("missing MultiFab data file '" + data_dir + fab_file.m_name + "'");
            }
        }
    };
 
    // Validate all required data headers before changing any runtime fields.
    for (int lev = 0; lev < m_nlevs; ++lev) {
        validate_multifab_header(MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "ustar"), lev, 2);
        validate_multifab_header(MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "wavg"), lev, 2);
        validate_multifab_header(MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "tstar"), lev, 1);
        validate_multifab_header(MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "qstar"), lev, 1);
        validate_multifab_header(MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "tsurf"), lev, 1);
        for (const auto &field : fieldmap) {
            validate_multifab_header(
                MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "_f_" + field.first), lev, 1);
        }
    }
 
    //
    // Move a field that was read on the checkpointed layout onto the live one.  When the two
    // layouts are identical this is the same local copy it always was, so an unchanged restart
    // stays bitwise identical; otherwise the data is redistributed.
    //
    // The redistribution is done in two passes because the checkpoint carries ghost cells and
    // a single ParallelCopy over grown source boxes would leave the outcome up to the order in
    // which overlapping sources happen to be applied.  The first pass seeds everything,
    // including the ghost cells outside a non-periodic physical boundary, which no amount of
    // valid-region copying can reach.  The second pass then lays the valid data over the top,
    // so wherever a cell is covered by a valid source cell that value wins.
    //
    auto redistribute = [&] (MultiFab& dst, const MultiFab& src, int ncomp, int lev)
    {
        if (dst.boxArray() == src.boxArray() && dst.DistributionMap() == src.DistributionMap()) {
            MultiFab::Copy(dst, src, 0, 0, ncomp, ng);
            return;
        }
        dst.ParallelCopy(src, 0, 0, ncomp, ng, ng);
        dst.ParallelCopy(src, 0, 0, ncomp, IntVect(0), ng, m_geom2d[lev].periodicity());
    };
 
    for (int lev = 0; lev < m_nlevs; lev++) {
        {
            MultiFab mf(chk_ba2d[lev],chk_dmap[lev],2,ng);
            VisMF::Read(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "ustar"));
            redistribute(*(u_star[lev]),mf,2,lev);
 
            VisMF::Read(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "wavg"));
            ensure_weight_factors(lev);
            redistribute(*wavg[lev],mf,2,lev);
        }
 
        MultiFab mf(chk_ba2d[lev],chk_dmap[lev],1,ng);
 
        VisMF::Read(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "tstar"));
        redistribute(*t_star[lev],mf,1,lev);
 
        VisMF::Read(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "qstar"));
        redistribute(*q_star[lev],mf,1,lev);
 
        VisMF::Read(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "tsurf"));
        redistribute(*t_surf[lev],mf,1,lev);
 
        for (auto &field : fieldmap) {
            VisMF::Read(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "_f_" + field.first));
            redistribute(*(fields[field.second.mf_ind][lev]), mf, 1, lev);
        }
    }
 
    m_weights_updated = chk_weights_updated;
    // Restore whether the surface models had already integrated a step when the checkpoint
    // was written.  Without this a restart would spend its first step treating the
    // checkpointed u*/t*/q* as unfilled and fall back to the MOST values instead.
    m_fields_are_valid = chk_fields_are_valid;
 
    deactivate_transient_field_maps();
    release_transient_surface_outputs();
 
    auto check_end = amrex::second() - check_start;
    ParallelDescriptor::ReduceRealMax(check_end,ParallelDescriptor::IOProcessorNumber());
    amrex::Print() << "    SurfaceModel Checkpoint load time = " << check_end << " seconds." << '\n';
}

register_field_map() [1/2]

void SurfaceModel::register_field_map	(	std::string	name,
		amrex::Vector< amrex::MultiFab * > &	lsm_lev_mf,
		amrex::Vector< amrex::MultiFab * > &	urb_lev_mf,
		bool	fill_boundary = false
	)

Registers a field map using per-level land and urban pointers.

Parameters

name	Common name used to retrieve the mapped field.
lsm_lev_mf	Land-model fields, indexed by AMR level.
urb_lev_mf	Urban-model fields, indexed by AMR level.
fill_boundary	Whether to fill mapped-field boundaries.

{
    // Register a field using explicit data pointers, rather than indices into lsm_data and urban_data
    AMREX_ALWAYS_ASSERT(lsm_lev_mf.size() > 0);
    amrex::Print() << " adding mapping between LSM<->Urban mfs <"<<lsm_lev_mf[0] << "," << urb_lev_mf[0] << "> to common surface name " << name << std::endl;
 
    AMREX_ALWAYS_ASSERT(lsm_lev_mf.size() == urb_lev_mf.size());
    if (fieldmap.find(name) == fieldmap.end()) {
 
        Field field;
        field.map = std::pair<int,int>(-1,-1);
        field.mf_ind = -1;
        field.fill_bound = fill_boundary;
        field.lsm_ptr = lsm_lev_mf;
        field.urb_ptr = urb_lev_mf;
 
        // Storage is allocated when a consumer activates this mapping.
        amrex::Vector<std::unique_ptr<amrex::MultiFab>> mf_lev(m_nlevs);
        fields.push_back(std::move(mf_lev));
        field.mf_ind = fields.size() - 1;
 
        amrex::Print() << "    -- created at ind = " << field.mf_ind << std::endl;
 
        fieldmap.insert({name, field});
    } else {
        return;
    }
}

Here is the call graph for this function:

register_field_map() [2/2]

void SurfaceModel::register_field_map	(	std::string	name,
		const std::pair< int, int > &	lsm_urb_map,
		bool	fill_boundary = false
	)

Registers a field map using land and urban field indices.

Parameters

name	Common name used to retrieve the mapped field.
lsm_urb_map	Pair of land and urban field indices.
fill_boundary	Whether to fill mapped-field boundaries.

{
    amrex::Print() << " adding mapping between LSM<->Urban fields <"<<lsm_urb_map.first << "," << lsm_urb_map.second << "> to common surface name " << name << std::endl;
 
    AMREX_ALWAYS_ASSERT(!(lsm_urb_map.first == -1 && lsm_urb_map.second == -1));
    if (fieldmap.find(name) == fieldmap.end()) {
 
        Field field;
        field.map = lsm_urb_map;
        field.mf_ind = -1;
        field.fill_bound = fill_boundary;
 
        // Storage is allocated when a consumer activates this mapping.
        amrex::Vector<std::unique_ptr<amrex::MultiFab>> mf_lev(m_nlevs);
        fields.push_back(std::move(mf_lev));
        field.mf_ind = fields.size() - 1;
 
        amrex::Print() << "    -- created at ind = " << field.mf_ind << std::endl;
 
        fieldmap.insert({name, field});
    } else {
        return;
    }
}

Here is the call graph for this function:

register_radiation_input()

void SurfaceModel::register_radiation_input	(	const std::string &	name,
		const std::pair< int, int > &	lsm_urb_map
	)

Registers a canonical radiation input mapping.

Parameters

name	Canonical radiation input name.
lsm_urb_map	Pair of land and urban model field indices.

{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(map.first == -1 && map.second == -1),
                                     "A radiation input must have a provider mapping");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        std::find(radnames.begin(), radnames.end(), name) != radnames.end(),
        "Unknown canonical radiation input: " + name);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(map.first >= -1 && map.second >= -1,
                                     "Radiation input indices must be -1 or nonnegative");
    if (map.first >= 0) {
        AMREX_ALWAYS_ASSERT(map.first < static_cast<int>(m_lsm_names.size()));
    }
    if (map.second >= 0) {
        AMREX_ALWAYS_ASSERT(map.second < static_cast<int>(m_urban_names.size()));
    }
    RadiationField field;
    field.map = map;
    field.weighted.resize(m_nlevs);
    radiation_input_map[name] = std::move(field);
    for (auto& fields_at_level : rad_fields) { fields_at_level.clear(); }
}

Here is the call graph for this function:

register_radiation_inputs()

void SurfaceModel::register_radiation_inputs

(

const std::unordered_map< std::string, std::pair< int, int >> &

input_map

)

Registers canonical radiation input mappings.

Parameters

input_map

Canonical names mapped to land and urban indices.

{
    for (const auto& entry : input_map) {
        register_radiation_input(entry.first, entry.second);
    }
}

register_radiation_output()

void SurfaceModel::register_radiation_output	(	const std::string &	name,
		const std::pair< int, int > &	lsm_urb_map
	)

Registers a canonical radiation output mapping.

Parameters

name	Canonical radiation output name.
lsm_urb_map	Pair of land and urban model field indices.

{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(map.first == -1 && map.second == -1),
                                     "A radiation output must have a provider mapping");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        std::find(rad_output_names.begin(), rad_output_names.end(), name) != rad_output_names.end(),
        "Unknown canonical radiation output: " + name);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(map.first >= -1 && map.second >= -1,
                                     "Radiation output indices must be -1 or nonnegative");
    if (map.first >= 0) {
        AMREX_ALWAYS_ASSERT(map.first < static_cast<int>(m_lsm_names.size()));
    }
    if (map.second >= 0) {
        AMREX_ALWAYS_ASSERT(map.second < static_cast<int>(m_urban_names.size()));
    }
    for (int lev = 0; lev < m_nlevs; ++lev) {
        if (map.first >= 0 && map.first < static_cast<int>(lsm_data_lev[lev].size())) {
            validate_radiation_output_layout(lev, lsm_data_lev[lev][map.first]);
        }
        if (map.second >= 0 && map.second < static_cast<int>(urban_data_lev[lev].size())) {
            validate_radiation_output_layout(lev, urban_data_lev[lev][map.second]);
        }
    }
    RadiationField field;
    field.map = map;
    radiation_output_map[name] = std::move(field);
    for (auto& fields_at_level : rad_output_fields) { fields_at_level.clear(); }
}

Here is the call graph for this function:

register_radiation_outputs()

void SurfaceModel::register_radiation_outputs

(

const std::unordered_map< std::string, std::pair< int, int >> &

output_map

)

Registers canonical radiation output mappings.

Parameters

output_map

Canonical names mapped to land and urban indices.

{
    for (const auto& entry : output_map) {
        register_radiation_output(entry.first, entry.second);
    }
}

release_transient_surface_outputs()

void SurfaceModel::release_transient_surface_outputs ( )

inlineprivate

Releases surface outputs that have no persistent consumer.

    {
        if (m_surface_outputs_enabled) { return; }
        for (int lev = 0; lev < m_nlevs; ++lev) {
            if (m_surface_layer_outputs_enabled &&
                m_provider_mode[lev] == SurfaceProviderMode::Both) {
                continue;
            }
            m_surface_outputs_requested[lev] = false;
            u_star[lev].reset();
            t_star[lev].reset();
            q_star[lev].reset();
            t_surf[lev].reset();
        }
    }

request_surface_layer_outputs()

void SurfaceModel::request_surface_layer_outputs ( )

inline

Enables blended surface outputs needed by the SurfaceLayer.

    {
        m_surface_layer_outputs_enabled = true;
        for (int lev = 0; lev < m_nlevs; ++lev) {
            if (m_provider_mode[lev] == SurfaceProviderMode::Both) {
                m_surface_outputs_requested[lev] = true;
                ensure_surface_outputs(lev);
            }
        }
    }

Here is the call graph for this function:

request_surface_outputs()

void SurfaceModel::request_surface_outputs ( const bool  persistent = true )

inline

Enables and materializes the surface outputs consumed by ERF.

    {
        m_surface_outputs_enabled = m_surface_outputs_enabled || persistent;
        for (int lev = 0; lev < m_nlevs; ++lev) {
            m_surface_outputs_requested[lev] = true;
            ensure_surface_outputs(lev);
        }
    }

Here is the call graph for this function:

set_field_map_pointers()

void SurfaceModel::set_field_map_pointers	(	const std::string &	name,
		int	lev,
		amrex::MultiFab *	lsm_mf,
		amrex::MultiFab *	urban_mf
	)

Updates a pointer-based field mapping for one level.

Parameters

name	Common name of the mapped field.
lev	AMR level to update.
lsm_mf	Land-model field pointer, or nullptr.
urban_mf	Urban-model field pointer, or nullptr.

{
    auto field_it = fieldmap.find(name);
    AMREX_ALWAYS_ASSERT(field_it != fieldmap.end());
 
    Field& field = field_it->second;
    AMREX_ALWAYS_ASSERT(field.map.first == -1 && field.map.second == -1);
    if (static_cast<int>(field.lsm_ptr.size()) < m_nlevs) {
        field.lsm_ptr.resize(m_nlevs, nullptr);
        field.urb_ptr.resize(m_nlevs, nullptr);
    }
    AMREX_ALWAYS_ASSERT(lev < m_nlevs);
    field.lsm_ptr[lev] = lsm_mf;
    field.urb_ptr[lev] = urban_mf;
}

Here is the call graph for this function:

set_model_data()

void SurfaceModel::set_model_data ( const int  lev, const amrex::Vector< amrex::MultiFab * >  model_data, const amrex::Vector< std::string > &  data_names, SurfaceModelType  type  )

inline

Registers model variables used in weighted averaging.

Parameters

lev	AMR level associated with the data.
model_data	Pointers to the model variables.
data_names	Names corresponding to model_data.
type	Whether the data belong to the land or urban model.

    {
        if (type == SurfaceModelType::LAND) {
            lsm_data_lev[lev].resize(model_data.size());
            weighted_lsm_data_lev[lev].clear();
            weighted_lsm_data_lev[lev].resize(model_data.size());
            if (lev == 0) {
                // only set names on first level, since they are the same for all levels
                m_lsm_names.resize(model_data.size());
                AMREX_ALWAYS_ASSERT(data_names.size() == model_data.size());
                for (int i = 0; i < static_cast<int>(data_names.size()); ++i) {
                    m_lsm_names[i] = data_names[i];
                }
            }
 
            for (int i = 0; i < static_cast<int>(lsm_data_lev[lev].size()); ++i) {
                AMREX_ALWAYS_ASSERT(model_data[i]);
                lsm_data_lev[lev][i] = model_data[i];
            }
            m_use_land = true;
            update_provider_mode(lev);
            for (auto& fields_at_level : rad_fields) { fields_at_level.clear(); }
            for (auto& fields_at_level : rad_output_fields) { fields_at_level.clear(); }
            for (auto& entry : radiation_input_map) {
                if (lev < static_cast<int>(entry.second.weighted.size())) {
                    entry.second.weighted[lev].reset();
                }
            }
        } else if (type == SurfaceModelType::URBAN) {
            urban_data_lev[lev].resize(model_data.size());
            weighted_urban_data_lev[lev].clear();
            weighted_urban_data_lev[lev].resize(model_data.size());
            if (lev == 0) {
                // only set names on first level, since they are the same for all levels
                m_urban_names.resize(model_data.size());
                AMREX_ALWAYS_ASSERT(data_names.size() == model_data.size());
                for (int i = 0; i < static_cast<int>(data_names.size()); ++i) {
                    m_urban_names[i] = data_names[i];
                }
            }
            for (int i = 0; i < static_cast<int>(urban_data_lev[lev].size()); ++i) {
                AMREX_ALWAYS_ASSERT(model_data[i]);
                urban_data_lev[lev][i] = model_data[i];
            }
            m_use_urban = true;
            update_provider_mode(lev);
            for (auto& fields_at_level : rad_fields) { fields_at_level.clear(); }
            for (auto& fields_at_level : rad_output_fields) { fields_at_level.clear(); }
            for (auto& entry : radiation_input_map) {
                if (lev < static_cast<int>(entry.second.weighted.size())) {
                    entry.second.weighted[lev].reset();
                }
            }
        }
    }

Here is the call graph for this function:

set_model_fields()

void SurfaceModel::set_model_fields ( SurfaceModelType  type, const amrex::Vector< int > &  field_indices, const bool  use_fluxes = true  )

inline

Selects the model fields that participate in weighted averaging.

When using fluxes, the first five elements map to u-stress, v-stress, heat flux, vapor flux, and surface temperature. Otherwise, the first four elements map to ustar, tstar, qstar, and surface temperature. A field index of -1 skips that field.

Parameters

type	Whether the fields belong to the land or urban model.
field_indices	Indices of the fields to average.
use_fluxes	Whether the selected fields are fluxes.

    {
        const int nfields = static_cast<int>(field_indices.size());
        if (!use_fluxes) {
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nfields >= 4, "Must have at least four fields to average (ustar, tstar, qstar, tsurf)");
        } else {
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nfields >= 5, "Must have at least five fields to average (u-stress, v-stress, hfx, qfx, tsurf)");
        }
 
        // Latch the export mode off the first model that registers, then make sure every
        // subsequent model agrees: either all of them supply fluxes, or all of them supply
        // u*, t*, q* directly.  Taking the mode from the caller rather than assuming the
        // default is what makes the non-flux path reachable at all.
        if (!m_export_fluxes_set) {
            m_export_fluxes = use_fluxes;
            m_export_fluxes_set = true;
        }
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(use_fluxes == m_export_fluxes,
            "All registered surface models must export either fluxes or u*/t*/q*, not a mix");
 
        const int nregistered = (type == SurfaceModelType::LAND)
            ? static_cast<int>(m_lsm_names.size())
            : static_cast<int>(m_urban_names.size());
        for (const int field : field_indices) {
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                field == -1 || (field >= 0 && field < nregistered),
                "Model field index must be -1 or within the registered provider data range");
        }
 
        if (type == SurfaceModelType::LAND) {
            lsm_fields = amrex::Vector<int>(field_indices);
            m_output_fields_registered[SurfaceModelType::LAND] = true;
            for (auto& fields_at_level : weighted_lsm_data_lev) {
                fields_at_level.clear();
            }
        } else if (type == SurfaceModelType::URBAN) {
            urban_fields = amrex::Vector<int>(field_indices);
            m_output_fields_registered[SurfaceModelType::URBAN] = true;
            for (auto& fields_at_level : weighted_urban_data_lev) {
                fields_at_level.clear();
            }
        }
    }

Here is the call graph for this function:

set_model_fluxes()

void SurfaceModel::set_model_fluxes ( const int  lev, const amrex::Vector< amrex::MultiFab * >  model_fluxes, const amrex::Vector< std::string > &  flux_names, SurfaceModelType  type  )

inline

Registers model fluxes as additional weighted-average inputs.

Parameters

lev	AMR level associated with the fluxes.
model_fluxes	Pointers to the model flux fields.
flux_names	Names corresponding to model_fluxes.
type	Whether the fluxes belong to the land or urban model.

    {
        amrex::Vector<amrex::MultiFab*>* model_data = nullptr;
        amrex::Vector<std::string>* model_names = nullptr;
 
        if (type == SurfaceModelType::LAND) {
            model_data = &lsm_data_lev[lev];
            model_names = &m_lsm_names;
        } else if (type == SurfaceModelType::URBAN) {
            model_data = &urban_data_lev[lev];
            model_names = &m_urban_names;
        }
 
        AMREX_ALWAYS_ASSERT(model_data != nullptr);
        AMREX_ALWAYS_ASSERT(model_names != nullptr);
        AMREX_ALWAYS_ASSERT(flux_names.size() == model_fluxes.size());
 
        const int data_size = static_cast<int>(model_data->size());
        model_data->resize(data_size + static_cast<int>(model_fluxes.size()));
        if (lev == 0) {
            model_names->resize(data_size + static_cast<int>(model_fluxes.size()));
        }
 
        for (int i = 0; i < static_cast<int>(model_fluxes.size()); ++i) {
            AMREX_ALWAYS_ASSERT(model_fluxes[i]);
            (*model_data)[data_size + i] = model_fluxes[i];
            if (lev == 0) {
                (*model_names)[data_size + i] = flux_names[i];
            }
        }
    }

Here is the call graph for this function:

update_provider_mode()

void SurfaceModel::update_provider_mode ( const int  lev )

inlineprivate

Updates the provider mode after model data are registered.

Parameters

lev	AMR level whose provider mode should be updated.

    {
        const bool has_land = !lsm_data_lev[lev].empty();
        const bool has_urban = !urban_data_lev[lev].empty();
        if (has_land && has_urban) {
            m_provider_mode[lev] = SurfaceProviderMode::Both;
            ensure_weight_factors(lev);
            if (m_surface_layer_outputs_enabled) {
                m_surface_outputs_requested[lev] = true;
                ensure_surface_outputs(lev);
            }
        } else if (has_land) {
            m_provider_mode[lev] = SurfaceProviderMode::LandOnly;
            wavg[lev].reset();
        } else if (has_urban) {
            m_provider_mode[lev] = SurfaceProviderMode::UrbanOnly;
            wavg[lev].reset();
        } else {
            m_provider_mode[lev] = SurfaceProviderMode::None;
        }
    }

Referenced by set_model_data().

Here is the call graph for this function:

Here is the caller graph for this function:

validate_radiation_output_layout()

void SurfaceModel::validate_radiation_output_layout	(	int	lev,
		const amrex::MultiFab *	mf
	)		const

Validates the layout contract for a canonical radiation output.

Parameters

lev	AMR level associated with the destination.
mf	Radiation output destination to validate.

{
    if (mf == nullptr) { return; }
    BoxList horizontal_boxes = mf->boxArray().boxList();
    // include the ghost cells here for SLM: surface values are exchanged at
    // k = 0 which is a ghost cell for SLM.
    const int k_grow = mf->nGrowVect()[2];
    for (Box& box : horizontal_boxes) {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            box.smallEnd(2) <= 0 && box.bigEnd(2) + k_grow >= 0,
            "Radiation output destination must contain the k=0 surface plane "
            "in its valid or ghost region");
        box.setRange(2, 0);
    }
    const BoxArray horizontal_ba(std::move(horizontal_boxes));
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        horizontal_ba == m_ba2d[lev] && mf->DistributionMap() == m_dmap[lev],
        "Radiation output destination must match the level horizontal layout");
}

Here is the call graph for this function:

weight_average_fields()

void SurfaceModel::weight_average_fields	(	int	lev,
		amrex::MultiFab * const
	)

Applies the current land and urban weights to selected model fields.

Parameters

lev	AMR level to update.
urban_frac	Urban fraction used to compute the weights.

{
    const SurfaceProviderMode mode = m_provider_mode[lev];
    const bool use_land = mode == SurfaceProviderMode::LandOnly ||
                          mode == SurfaceProviderMode::Both;
    const bool use_urban = mode == SurfaceProviderMode::UrbanOnly ||
                           mode == SurfaceProviderMode::Both;
 
    for (auto &field : fieldmap)
    {
        if (!field.second.active) { continue; }
        int mf_idx = field.second.mf_ind;
        AMREX_ASSERT(mf_idx != -1);
 
        int lsm_idx = field.second.map.first;
        int urb_idx = field.second.map.second;
 
        // whether we have a valid LSM multifab
        bool valid_land = (use_land &&
                           lsm_idx != -1 &&
                           lsm_data_lev[lev][lsm_idx]);
 
        // whether we have a valid urban multifab
        bool valid_urban = (use_urban &&
                            urb_idx != -1 &&
                            urban_data_lev[lev][urb_idx]);
 
        bool use_mf = false;
        if (lsm_idx == -1 && urb_idx == -1) {
            // use explicit MF ptrs rather than indices
            use_mf = true;
            valid_land = (use_land && field.second.lsm_ptr[lev]);
            valid_urban = (use_urban && field.second.urb_ptr[lev]);
        }
 
        for (MFIter mfi(*fields[mf_idx][lev], TileNoZ()); mfi.isValid(); ++mfi)
        {
            Box tbx = mfi.tilebox();
            const auto weights_arr = (mode == SurfaceProviderMode::Both)
                ? wavg[lev]->const_array(mfi) : Array4<const Real>{};
 
            // Calculate weight average into output
            auto output_arr = fields[mf_idx][lev]->array(mfi);
            const amrex::MultiFab *lsm_mf = (use_mf) ? field.second.lsm_ptr[lev] : (lsm_idx != -1 ? lsm_data_lev[lev][lsm_idx] : nullptr);
            const amrex::MultiFab *urb_mf = (use_mf) ? field.second.urb_ptr[lev] : (urb_idx != -1 ? urban_data_lev[lev][urb_idx] : nullptr);
 
            auto lsm_data_arr = (valid_land) ? lsm_mf->const_array(mfi) : Array4<const Real>{};
            auto urban_data_arr = (valid_urban) ? urb_mf->const_array(mfi) : Array4<const Real>{};
 
            if (valid_land && !valid_urban) {
                // use solely land value
                ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    output_arr(i, j, k) = lsm_data_arr(i, j, k);
                });
            } else if (!valid_land && valid_urban) {
                // use solely urban value
                ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    output_arr(i, j, k) = urban_data_arr(i, j, k);
                });
            } else {
                ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    Real land = (lsm_data_arr) ? lsm_data_arr(i, j, k) * weights_arr(i, j, 0, SurfaceModelType::LAND) : 0.0;
                    Real urb  = (urban_data_arr) ? urban_data_arr(i, j, k) * weights_arr(i, j, 0, SurfaceModelType::URBAN) : 0.0;
 
                    output_arr(i, j, k) = land + urb;
                });
            }
        }
        if (field.second.fill_bound) {
            fields[mf_idx][lev]->FillBoundary(m_geom[lev].periodicity());
        }
        //outputs[output_field]->FillBoundary(comp, 1, m_geom[lev].periodicity());
    }
}

Here is the call graph for this function:

weight_model_field()

void SurfaceModel::weight_model_field	(	int	lev,
		const amrex::MultiFab *	source,
		amrex::MultiFab *	weighted,
		SurfaceModelType	type
	)

Writes a weighted copy of a model field.

Parameters

lev	AMR level whose weights are used.
source	Model-owned source field.
weighted	SurfaceModel-owned destination field.
type	Provider type owning the source field.

{
    AMREX_ASSERT(source != nullptr);
    AMREX_ASSERT(weighted != nullptr);
 
    for (MFIter mfi(*source, TileNoZ()); mfi.isValid(); ++mfi)
    {
        Box tbx = mfi.tilebox();
        const auto source_arr = source->const_array(mfi);
        auto weighted_arr = weighted->array(mfi);
        const bool use_weights = m_provider_mode[lev] == SurfaceProviderMode::Both;
        const auto weights_arr = use_weights
            ? wavg[lev]->const_array(mfi) : Array4<const Real>{};
 
        ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            weighted_arr(i, j, k) = use_weights
                ? source_arr(i, j, k) * weights_arr(i, j, 0, type)
                : source_arr(i, j, k);
        });
    }
 
    weighted->FillBoundary(m_geom[lev].periodicity());
}

Here is the call graph for this function:

write_output()

void SurfaceModel::write_output	(	const int &	finest_lev,
		const amrex::Real &	time,
		const std::string &	plot_prefix,
		const amrex::Vector< int > &	level_steps,
		const amrex::Vector< amrex::IntVect > &	ref_ratio
	)

Writes registered surface-model fields to plotfile output.

Parameters

finest_lev	Finest AMR level to write.
time	Simulation time associated with the output.
plot_prefix	Output plotfile prefix.
level_steps	Number of steps represented at each level.
ref_ratio	Refinement ratios between adjacent levels.

{
    request_surface_outputs(false);
    activate_all_field_maps(false);
    std::string plotfilename = amrex::Concatenate(plot_prefix + "2D_", level_steps[0], 5);
 
    const int nfields = (m_export_fluxes) ? 5 : 4;
    // Surface outputs, land mask, urban fraction, and mapped fields.
    const int noutput = nfields + 2 + fieldmap.size();
    IntVect ng(0, 0, 0);
 
    amrex::Vector<std::string> varnames(nfields);
 
    amrex::Vector<amrex::MultiFab> fab(finest_lev+1);
    for (int lev = 0; lev <= finest_lev; lev++) {
        calculate_weight_average(lev, m_last_urban_frac[lev], false);
        amrex::MultiFab* const outputs[] = {u_star[lev].get(), t_star[lev].get(), q_star[lev].get(), t_surf[lev].get()};
        fab[lev].define(outputs[0]->boxArray(), m_dmap[lev], noutput, ng);
 
        //fab[lev].setVal(0.0);
 
        // Output weighted surface fluxes into ustar, tstar, qstar and surface temperature into tsurf
        // TODO: make sure grids of urban and LSM inputs match
        for (int field=0; field < nfields; field++)
        {
            int output_field = (m_export_fluxes && field > 0) ? field - 1 : field;
            int comp = (m_export_fluxes && field < 2) ? field : 0;
 
            MultiFab::Copy(fab[lev], *(outputs[output_field]), comp, field, 1, ng);
 
            varnames[field] = field_names[field];
        }
 
        int nout = nfields;
 
        MultiFab lmask_tmp = amrex::ToMultiFab(*m_lmask[lev]); // iMultiFab -> MultiFab
        MultiFab::Copy(fab[lev], lmask_tmp, 0, nout, 1, ng);
        nout++;
        MultiFab::Copy(fab[lev], *get_wavg_factors(lev), SurfaceModelType::URBAN, nout, 1, ng);
        nout++;
 
        if (lev == 0) {
            varnames.push_back("lmask");
            varnames.push_back("urb_frac");
        }
 
        // Add any mapped fields
        for (auto &field : fieldmap) {
            MultiFab::Copy(fab[lev], *(fields[field.second.mf_ind][lev]), 0, nout, 1, ng);
            if (lev == 0) {
                varnames.push_back(field.first);
            }
            nout++;
        }
 
        AMREX_ALWAYS_ASSERT(varnames.size() == noutput);
    }
 
    amrex::WriteMultiLevelPlotfile(plotfilename, finest_lev+1, GetVecOfConstPtrs(fab), varnames, m_geom2d, time, level_steps, ref_ratio);
    deactivate_transient_field_maps();
    release_transient_surface_outputs();
}

Here is the call graph for this function:

WriteCheckpoint()

void SurfaceModel::WriteCheckpoint

(

const std::string &

checkpointname

)

Writes surface-model state to a checkpoint file.

Parameters

checkpointname

Checkpoint file name.

{
    request_surface_outputs(false);
    activate_all_field_maps(false);
    auto check_start = amrex::second();
 
    // write header
    if (ParallelDescriptor::IOProcessor()) {
 
        amrex::Print() << " Writing SurfaceModel checkpoint " << std::endl;
 
        std::string HeaderFileName(checkpointname + "/SurfaceModel_Header");
        VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
        std::ofstream HeaderFile;
        HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
        HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out   |
                                                std::ofstream::trunc |
                                                std::ofstream::binary);
        if(! HeaderFile.good()) {
            FileOpenFailed(HeaderFileName);
        }
 
        HeaderFile.precision(17);
 
        // write out title line
        HeaderFile << "Checkpoint file for SurfaceModel\n";
 
        // write out number of levels
        HeaderFile << m_nlevs << "\n";
 
        // write out flags
        HeaderFile << m_use_urban << "\n";
        HeaderFile << m_use_land << "\n";
        HeaderFile << m_export_fluxes << "\n";
        HeaderFile << m_weights_updated << "\n";
        HeaderFile << m_fields_are_valid << "\n";
        HeaderFile << "\n";
 
        // Write box arrays
        for (int lev = 0; lev < m_nlevs; lev++) {
            m_ba[lev].writeOn(HeaderFile);
            HeaderFile << '\n';
 
        }
        HeaderFile << '\n';
        for (int lev = 0; lev < m_nlevs; lev++) {
            m_ba2d[lev].writeOn(HeaderFile);
            HeaderFile << '\n';
        }
        HeaderFile << '\n';
 
        // Write out fields
        // LSM fields
        for (int i = 0; i < static_cast<int>(lsm_fields.size()); ++i) {
            HeaderFile << lsm_fields[i] << " ";
        }
        HeaderFile << '\n';
 
        // Urban fields
        for (int i = 0; i < static_cast<int>(urban_fields.size()); ++i) {
            HeaderFile << urban_fields[i] << " ";
        }
        HeaderFile << '\n';
 
        // Field mapping
        HeaderFile << '\n';
        HeaderFile << fieldmap.size() << "\n";
        for (auto &field : fieldmap) {
            // pointer fields are reconstructed at load
            HeaderFile << field.first << " " << field.second.mf_ind << " " << field.second.map.first << " " << field.second.map.second << " " << field.second.fill_bound << '\n';
        }
        HeaderFile << '\n';
    }
 
    amrex::ParallelDescriptor::Barrier();
 
    const std::string prefix = "SurfaceModel_";
 
    // Radiation input fields are derived from provider data and are intentionally
    // rebuilt after restart rather than written to the checkpoint.
    for (int lev = 0; lev < m_nlevs; lev++) {
        IntVect ng(1,1,0);
        calculate_weight_average(lev, m_last_urban_frac[lev], false);
 
        {
            MultiFab mf(m_ba2d[lev],m_dmap[lev],2,ng);
            MultiFab::Copy(mf,*u_star[lev],0,0,2,ng);
            VisMF::Write(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "ustar"));
 
            MultiFab::Copy(mf,*get_wavg_factors(lev),0,0,2,ng);
            VisMF::Write(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "wavg"));
        }
 
        MultiFab mf(m_ba2d[lev],m_dmap[lev],1,ng);
 
        MultiFab::Copy(mf,*t_star[lev],0,0,1,ng);
        VisMF::Write(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "tstar"));
 
        MultiFab::Copy(mf,*q_star[lev],0,0,1,ng);
        VisMF::Write(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "qstar"));
 
        MultiFab::Copy(mf,*t_surf[lev],0,0,1,ng);
        VisMF::Write(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "tsurf"));
 
        for (auto &field : fieldmap) {
            MultiFab::Copy(mf, *(fields[field.second.mf_ind][lev]), 0, 0, 1, ng);
            VisMF::Write(mf, MultiFabFileFullPrefix(lev, checkpointname, "Level_", prefix + "_f_" + field.first));
        }
    }
 
    auto check_end = amrex::second() - check_start;
    deactivate_transient_field_maps();
    release_transient_surface_outputs();
    ParallelDescriptor::ReduceRealMax(check_end,ParallelDescriptor::IOProcessorNumber());
    amrex::Print() << "    SurfaceModel Checkpoint write time = " << check_end << " seconds." << '\n';
}

Member Data Documentation

field_names

const std::vector<std::string> SurfaceModel::field_names

inlinestaticprivate

Initial value:

= {
        "uflux",
        "vflux",
        "hfx",
        "qfx",
        "tskin",
    }

fieldmap

std::unordered_map<std::string, Field> SurfaceModel::fieldmap

private

Referenced by get_field(), and initialize_for_level().

fields

amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab> > > SurfaceModel::fields

private

Referenced by get_field(), and initialize_for_level().

lsm_data_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceModel::lsm_data_lev

private

Referenced by build_radiation_fieldlist(), get_surface_flux_view(), get_weighted_model_data(), initialize_for_level(), set_model_data(), set_model_fluxes(), SurfaceModel(), and update_provider_mode().

lsm_fields

amrex::Vector<int> SurfaceModel::lsm_fields

private

Referenced by get_surface_flux_view(), get_weighted_model_data(), set_model_fields(), and SurfaceModel().

m_ba

amrex::Vector<amrex::BoxArray> SurfaceModel::m_ba

private

Referenced by initialize_for_level().

m_ba2d

amrex::Vector<amrex::BoxArray> SurfaceModel::m_ba2d

private

Referenced by build_radiation_fieldlist(), ensure_surface_outputs(), ensure_weight_factors(), initialize_for_level(), and SurfaceModel().

m_dmap

amrex::Vector<amrex::DistributionMapping> SurfaceModel::m_dmap

private

Referenced by build_radiation_fieldlist(), ensure_surface_outputs(), ensure_weight_factors(), and initialize_for_level().

m_export_fluxes

bool SurfaceModel::m_export_fluxes = true

private

Referenced by are_fluxes(), get_surface_flux_view(), get_weighted_model_data(), and set_model_fields().

m_export_fluxes_set

bool SurfaceModel::m_export_fluxes_set = false

private

Referenced by set_model_fields().

m_fields_are_valid

bool SurfaceModel::m_fields_are_valid = false

private

Referenced by fields_are_valid(), and mark_fields_valid().

m_geom

amrex::Vector<amrex::Geometry> SurfaceModel::m_geom

private

Referenced by build_radiation_fieldlist(), and initialize_for_level().

m_geom2d

amrex::Vector<amrex::Geometry> SurfaceModel::m_geom2d

private

Referenced by initialize_for_level(), and SurfaceModel().

m_last_urban_frac

amrex::Vector<amrex::MultiFab*> SurfaceModel::m_last_urban_frac

private

Referenced by initialize_for_level(), and SurfaceModel().

m_lmask

amrex::Vector<amrex::iMultiFab*> SurfaceModel::m_lmask

private

Referenced by initialize_for_level(), and SurfaceModel().

m_lsm_names

amrex::Vector<std::string> SurfaceModel::m_lsm_names

private

Referenced by set_model_data(), set_model_fields(), and set_model_fluxes().

m_nlevs

int SurfaceModel::m_nlevs

private

Referenced by initialize_for_level(), release_transient_surface_outputs(), request_surface_layer_outputs(), request_surface_outputs(), and SurfaceModel().

m_output_fields_registered

amrex::GpuArray<bool, 2> SurfaceModel::m_output_fields_registered {{false, false}}

private

Referenced by set_model_fields(), and SurfaceModel().

m_provider_mode

amrex::Vector<SurfaceProviderMode> SurfaceModel::m_provider_mode

private

Referenced by build_radiation_fieldlist(), ensure_weight_factors(), get_surface_flux_view(), get_weighted_model_data(), initialize_for_level(), release_transient_surface_outputs(), request_surface_layer_outputs(), SurfaceModel(), and update_provider_mode().

m_surface_layer_outputs_enabled

bool SurfaceModel::m_surface_layer_outputs_enabled = false

private

Referenced by release_transient_surface_outputs(), request_surface_layer_outputs(), and update_provider_mode().

m_surface_outputs_enabled

bool SurfaceModel::m_surface_outputs_enabled = false

private

Referenced by initialize_for_level(), release_transient_surface_outputs(), and request_surface_outputs().

m_surface_outputs_requested

amrex::Vector<int> SurfaceModel::m_surface_outputs_requested

private

Referenced by initialize_for_level(), release_transient_surface_outputs(), request_surface_layer_outputs(), request_surface_outputs(), SurfaceModel(), and update_provider_mode().

m_urban_names

amrex::Vector<std::string> SurfaceModel::m_urban_names

private

Referenced by set_model_data(), set_model_fields(), and set_model_fluxes().

m_use_land

bool SurfaceModel::m_use_land = false

private

Referenced by set_model_data().

m_use_urban

bool SurfaceModel::m_use_urban = false

private

Referenced by set_model_data().

m_weights_updated

bool SurfaceModel::m_weights_updated = false

private

q_star

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceModel::q_star

private

Referenced by ensure_surface_outputs(), get_qstar(), get_surface_flux_view(), initialize_for_level(), release_transient_surface_outputs(), and SurfaceModel().

rad_fields

amrex::Vector<amrex::Vector<const amrex::MultiFab*> > SurfaceModel::rad_fields

private

Referenced by build_radiation_fieldlist(), get_radiation_fields(), initialize_for_level(), set_model_data(), and SurfaceModel().

rad_output_fields

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceModel::rad_output_fields

private

Referenced by initialize_for_level(), set_model_data(), and SurfaceModel().

rad_output_names

const std::vector<std::string> SurfaceModel::rad_output_names

inlinestaticprivate

Initial value:

= {
        "cos_zenith_angle", "sw_flux_dn", "sw_flux_dn_dir_vis",
        "sw_flux_dn_dir_nir", "sw_flux_dn_dif_vis", "sw_flux_dn_dif_nir",
        "lw_flux_dn"}

radiation_input_map

std::unordered_map<std::string, RadiationField> SurfaceModel::radiation_input_map

private

Referenced by build_radiation_fieldlist(), initialize_for_level(), and set_model_data().

radiation_output_map

std::unordered_map<std::string, RadiationField> SurfaceModel::radiation_output_map

private

radnames

const std::vector<std::string> SurfaceModel::radnames = {"tskin", "emiss", "albedo_vis", "albedo_nir", "albedo_vis_diff", "albedo_nir_diff"}

inlinestaticprivate

Referenced by build_radiation_fieldlist().

t_star

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceModel::t_star

private

Referenced by ensure_surface_outputs(), get_surface_flux_view(), get_tstar(), initialize_for_level(), release_transient_surface_outputs(), and SurfaceModel().

t_surf

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceModel::t_surf

private

Referenced by ensure_surface_outputs(), get_tsurf(), initialize_for_level(), release_transient_surface_outputs(), and SurfaceModel().

u_star

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceModel::u_star

private

Referenced by ensure_surface_outputs(), get_surface_flux_view(), get_ustar(), initialize_for_level(), release_transient_surface_outputs(), and SurfaceModel().

urban_data_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceModel::urban_data_lev

private

Referenced by build_radiation_fieldlist(), get_surface_flux_view(), get_weighted_model_data(), initialize_for_level(), set_model_data(), set_model_fluxes(), SurfaceModel(), and update_provider_mode().

urban_fields

amrex::Vector<int> SurfaceModel::urban_fields

private

Referenced by get_surface_flux_view(), get_weighted_model_data(), set_model_fields(), and SurfaceModel().

urban_frac_lev

amrex::Vector<amrex::MultiFab*> SurfaceModel::urban_frac_lev

private

Referenced by initialize_for_level(), and SurfaceModel().

wavg

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceModel::wavg

private

Referenced by build_radiation_fieldlist(), ensure_weight_factors(), get_wavg_factors(), initialize_for_level(), SurfaceModel(), and update_provider_mode().

weighted_lsm_data_lev

amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab> > > SurfaceModel::weighted_lsm_data_lev

private

Referenced by get_weighted_model_data(), initialize_for_level(), set_model_data(), set_model_fields(), and SurfaceModel().

weighted_urban_data_lev

amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab> > > SurfaceModel::weighted_urban_data_lev

private

Referenced by get_weighted_model_data(), initialize_for_level(), set_model_data(), set_model_fields(), and SurfaceModel().

---

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

• Source/BoundaryConditions/ERF_SurfaceModel.H
• Source/BoundaryConditions/ERF_SurfaceModel.cpp