#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:

[legend]

Classes
struct	PBLHColumns

Public Types
enum class	FluxCalcType { MOENG = 0 , CUSTOM , BULK_COEFF , ROTATE , RICO }

enum class	ThetaCalcType { ADIABATIC = 0 , HEAT_FLUX , SURFACE_TEMPERATURE }

enum class	MoistCalcType { ADIABATIC = 0 , MOISTURE_FLUX , SURFACE_MOISTURE }

enum class	RoughCalcType { CONSTANT = 0 , CHARNOCK , MODIFIED_CHARNOCK , DONELAN , WAVE_COUPLED }

enum class	PBLHeightCalcType { None , MYNN25 , YSU , MRF , YSUNew }

Public Member Functions
	SurfaceLayer (amrex::Orientation face, const amrex::Vector< amrex::Geometry > &geom, bool &use_rot_surface_flux, std::string a_pp_prefix, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qv_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &z_phys_nd, const amrex::Vector< amrex::Vector< amrex::Real >> &zlevels_stag, const MeshType &a_mesh_type, const TerrainType &a_terrain_type, const TurbChoice &a_turb_choice, amrex::Real a_rdOcp, double start_low_time, double final_low_time, double low_time_interval=0.0, const amrex::Vector< const eb_ * > &eb_vec={}, SurfaceModel *surf_model=nullptr, bool a_has_immersed_buildings=false)

void	make_SurfaceLayer_at_level (const int &lev, int nlevs, const amrex::Vector< amrex::MultiFab * > &mfv, std::unique_ptr< amrex::MultiFab > &Theta_prim, std::unique_ptr< amrex::MultiFab > &Qv_prim, std::unique_ptr< amrex::MultiFab > &Qr_prim, std::unique_ptr< amrex::MultiFab > &z_phys_nd, amrex::MultiFab Hwave, amrex::MultiFab Lwave, amrex::MultiFab eddyDiffs, amrex::Vector< amrex::MultiFab > lsm_data, amrex::Vector< std::string > lsm_data_name, amrex::Vector< amrex::MultiFab * > lsm_flux, amrex::Vector< std::string > lsm_flux_name, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &sst_lev, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &tsk_lev, amrex::Vector< std::unique_ptr< amrex::iMultiFab >> &lmask_lev)

void	set_immersed_surface_params (const int &lev)

void	update_fluxes (const int &lev, const double &elapsed_time, const double &elapsed_time_since_start_low, amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &walldist, int max_iters=100)

template<typename FluxIter >
void	compute_fluxes (const int &lev, const int &max_iters, amrex::MultiFab &cons_in, const FluxIter &most_flux, bool is_land)

void	init_tke_from_ustar (const int &lev, amrex::MultiFab &cons, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const amrex::Real tkefac=one, const amrex::Real zscale=amrex::Real(700.0))

void	fill_planar_boundary (const int &lev, amrex::MultiFab &mf)

amrex::Real	surface_sum (const int &lev, const amrex::MultiFab &mf, int comp=0) const

void	impose_SurfaceLayer_bcs (const int &lev, amrex::Vector< const amrex::MultiFab * > mfs, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Tau_lev, amrex::MultiFab xheat_flux, amrex::MultiFab yheat_flux, amrex::MultiFab zheat_flux, amrex::MultiFab xqv_flux, amrex::MultiFab yqv_flux, amrex::MultiFab zqv_flux, const amrex::MultiFab z_phys, amrex::MultiFab mf_surface_fluxes=nullptr)

void	impose_SurfaceLayer_bcs_EB (const int &lev, amrex::Vector< const amrex::MultiFab * > mfs, amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab >>> &Tau_lev, amrex::MultiFab xheat_flux, amrex::MultiFab yheat_flux, amrex::MultiFab zheat_flux, amrex::MultiFab xqv_flux, amrex::MultiFab yqv_flux, amrex::MultiFab zqv_flux)

template<typename FluxCalc >
void	compute_SurfaceLayer_bcs (const int &lev, amrex::Vector< const amrex::MultiFab * > mfs, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Tau_lev, amrex::MultiFab xheat_flux, amrex::MultiFab yheat_flux, amrex::MultiFab zheat_flux, amrex::MultiFab xqv_flux, amrex::MultiFab yqv_flux, amrex::MultiFab zqv_flux, const amrex::MultiFab z_phys, const FluxCalc &flux_comp, amrex::MultiFab mf_surface_fluxes=nullptr)

template<typename FluxCalc >
void	compute_SurfaceLayer_bcs_EB (const int &lev, amrex::Vector< const amrex::MultiFab * > mfs, amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab >>> &Tau_lev, amrex::MultiFab xheat_flux, amrex::MultiFab yheat_flux, amrex::MultiFab zheat_flux, amrex::MultiFab xqv_flux, amrex::MultiFab yqv_flux, amrex::MultiFab zqv_flux, const FluxCalc &flux_comp)

void	compute_sfc_params_from_lsm_fluxes (const int &lev, amrex::MultiFab &cons_in)

void	set_sfc_params_from_file_fluxes (const int &lev, const amrex::MultiFab &cons_in, const amrex::Real &tflux_in, const amrex::Real &qflux_in, const amrex::Real &sfc_tau)

void	fill_tsurf_with_sst_and_tsk (const int &lev, const double &time)

void	fill_tsurf_with_sfc_sst (const int &lev, const double &time, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	fill_tsurf_with_coupled_sst (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	fill_tsurf_with_abs_surf_temp (const int &lev, const double &time, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	fill_tsurf_with_skin_temperature (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	fill_qsurf_with_skin_moisture (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	fill_qsurf_with_qsat (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	set_pblh (const int &lev, const amrex::MultiFab &pblh_in)

void	update_sfc_time_index (const amrex::Real &time)

amrex::Real	interpolate_sfc_column (const amrex::Real &time, int col) const

void	get_lsm_tsurf (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

void	update_pblh (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars, amrex::MultiFab *z_phys_cc, const MoistureComponentIndices &moisture_indices)

void	fill_pblh_from_coarser (const int &lev)

void	inject_pblh_from_coarser (const int &lev, amrex::MultiFab &fine) const

template<typename PBLHeightEstimator >
void	compute_pblh (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars, amrex::MultiFab *z_phys_cc, const PBLHeightEstimator &est, const MoistureComponentIndices &moisture_indice)

void	read_custom_roughness (const int &lev, const std::string &fname)

void	update_surf_temp (const double &time)

void	update_mac_ptrs (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_old, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Theta_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qv_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qr_prim)

amrex::MultiFab *	get_u_star (const int &lev)

amrex::MultiFab *	get_w_star (const int &lev)

bool	computes_w_star () const

bool	computes_pblh () const

amrex::MultiFab *	get_t_star (const int &lev)

amrex::MultiFab *	get_q_star (const int &lev)

amrex::MultiFab *	get_olen (const int &lev)

amrex::MultiFab *	get_pblh (const int &lev)

const amrex::MultiFab *	get_mac_avg (const int &lev, int comp)

bool	mac_avg_is_time_averaged () const

int	get_num_mac_avg () const

bool	mac_avg_is_initialized (const int &lev) const

void	set_mac_avg_initialized (const int &lev)

amrex::MultiFab *	get_mac_avg_ptr (const int &lev, int comp)

amrex::Vector< amrex::Real >	get_mac_plane_avg (const int &lev) const

bool	set_mac_plane_avg (const int &lev, const amrex::Vector< amrex::Real > &pavg)

amrex::MultiFab *	get_t_surf (const int &lev)

void	set_skin_temperature (const int &lev, const amrex::MultiFab *t_skin)

std::string	skin_temperature_conflict () const

void	set_skin_moisture (const int &lev, const amrex::MultiFab *moisture, bool use_resistances=false, amrex::Real hs=amrex::Real(0))

std::string	skin_moisture_conflict () const

bool	set_default_roughness (amrex::Real z0)

amrex::Real	roughness () const

bool	use_sfc_fluxes () const

bool	rotates_surface_fluxes () const

void	set_t_surf (const int &lev, const amrex::Real tsurf)

amrex::MultiFab *	get_q_surf (const int &lev)

void	set_q_surf (const int &lev, const amrex::Real qsurf)

amrex::MultiFab *	get_surface_diagnostic_source (const int &lev)

amrex::MultiFab *	get_zref (const int &lev)

amrex::Real	get_zref_min (const int &lev)

amrex::MultiFab *	get_z0 (const int &lev)

bool	have_variable_sea_roughness ()

amrex::iMultiFab *	get_lmask (const int &lev)

int	lmask_min_reduce (amrex::iMultiFab &lmask, const int &nghost)

void	update_sst_ptr (const int lev, const int itime, amrex::MultiFab *sst_ptr)

void	update_tsk_ptr (const int lev, const int itime, amrex::MultiFab *tsk_ptr)

void	update_coupled_sst_ptr (const int lev, amrex::MultiFab sst_ptr, amrex::iMultiFab valid_ptr)

void	set_coupled_sst_active (const bool active)

void	set_surface_layer_faces (const amrex::GpuArray< int, AMREX_SPACEDIM *2 > &active_faces)

void	fill_lateral_surface_parameter_ghosts (const int &lev, amrex::MultiFab *selected_field=nullptr)

template<typename FluxIter >
void	compute_fluxes (const int &lev, const int &max_iters, MultiFab &, const FluxIter &most_flux, bool is_land)

template<typename FluxCalc >
void	compute_SurfaceLayer_bcs (const int &lev, Vector< const MultiFab * > mfs, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab xheat_flux, MultiFab yheat_flux, MultiFab zheat_flux, MultiFab xqv_flux, MultiFab yqv_flux, MultiFab zqv_flux, const MultiFab z_phys, const FluxCalc &flux_comp, MultiFab mf_surface_fluxes)

template<typename FluxCalc >
void	compute_SurfaceLayer_bcs_EB (const int &lev, Vector< const MultiFab * > mfs, Vector< Vector< std::unique_ptr< MultiFab >>> &Tau_EB, [[maybe_unused]] MultiFab xheat_flux, [[maybe_unused]] MultiFab yheat_flux, MultiFab Hfx3_EB, [[maybe_unused]] MultiFab xqv_flux, [[maybe_unused]] MultiFab yqv_flux, [[maybe_unused]] MultiFab zqv_flux, const FluxCalc &flux_comp)

template<typename PBLHeightEstimator >
void	compute_pblh (const int &lev, Vector< Vector< MultiFab >> &vars, MultiFab *z_phys_cc, const PBLHeightEstimator &est, const MoistureComponentIndices &moisture_indices)

Static Public Member Functions
static amrex::Vector< amrex::Vector< amrex::Real > >	read_cols (const std::string &fname, const int skip_nlines=1)

Public Attributes
FluxCalcType	flux_type {FluxCalcType::MOENG}

ThetaCalcType	theta_type {ThetaCalcType::ADIABATIC}

MoistCalcType	moist_type {MoistCalcType::ADIABATIC}

RoughCalcType	rough_type_land {RoughCalcType::CONSTANT}

RoughCalcType	rough_type_sea {RoughCalcType::CHARNOCK}

PBLHeightCalcType	pblh_type {PBLHeightCalcType::None}

Private Member Functions
void	define_pblh_columns (const int &lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)

Private Attributes
amrex::Orientation	m_face

std::string	m_pp_prefix

amrex::Vector< amrex::Geometry >	m_geom

bool	m_rotate = false

amrex::GpuArray< int, AMREX_SPACEDIM *2 >	m_surface_layer_faces {}

double	m_start_low_time

double	m_final_low_time

double	m_low_time_interval

bool	m_include_wstar = false

amrex::Real	z0_const {amrex::Real(0.1)}

bool	m_z0_given {false}

amrex::Real	default_land_surf_temp {amrex::Real(300.)}

amrex::Real	surf_temp {amrex::Real(-1.)}

bool	surf_temp_is_absolute {false}

amrex::Real	surf_heating_rate {0}

amrex::Real	surf_temp_flux {0}

amrex::Real	default_land_surf_moist {zero}

amrex::Real	surf_moist {amrex::Real(-1.)}

amrex::Real	surf_moist_flux {0}

amrex::Real	custom_ustar {0}

amrex::Real	custom_tstar {0}

amrex::Real	custom_qstar {0}

amrex::Real	custom_rhosurf {0}

bool	specified_rho_surf {false}

amrex::Real	cnk_a {amrex::Real(0.0185)}

bool	smooth_flow_visc {true}

amrex::Real	depth {amrex::Real(30.0)}

amrex::Vector< amrex::MultiFab >	z_0

bool	m_var_z0 {false}

amrex::Real	rico_theta_z0 {amrex::Real(298.0)}

amrex::Real	rico_qsat_z0 {amrex::Real(0.001)}

bool	use_moisture

bool	m_has_lsm_fluxes = false

bool	m_has_lsm_tsurf = false

int	m_lsm_tsurf_indx = -1

bool	m_lsm_tsurf_is_theta = false

bool	m_use_lsm_skin_temperature = false

bool	m_use_sfc_fluxes = false

bool	m_use_sfc_sst = false

int	sfc_time_ind = 0

amrex::Vector< amrex::Vector< amrex::Real > >	sfc

amrex::Real	sfc_qflux = zero

amrex::Real	sfc_tflux = zero

amrex::Real	sfc_ustar = zero

bool	m_use_coupled_sst = false

amrex::Real	m_Cd = zero

amrex::Real	m_Ch = zero

amrex::Real	m_Cq = zero

bool	m_ignore_sst = false

amrex::Vector< const eb_ * >	m_eb_vec

TerrainType	m_terrain_type

bool	m_has_immersed_buildings = false

amrex::Real	m_rdOcp = RdoCp

MOSTAverage	m_ma

SurfaceModel *	m_surf_model = nullptr

bool	use_surface_model = false

bool	surf_model_fluxes = false

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

amrex::Vector< std::unique_ptr< amrex::MultiFab > >	w_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 > >	olen

amrex::Vector< std::unique_ptr< amrex::MultiFab > >	pblh

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

amrex::Vector< std::unique_ptr< amrex::MultiFab > >	q_surf

amrex::Vector< PlanarBoundary >	m_planar_bndry

amrex::Vector< PBLHColumns >	m_pblh_columns

amrex::Vector< std::unique_ptr< amrex::MultiFab > >	surface_diagnostic_source

amrex::Vector< amrex::Vector< amrex::MultiFab * > >	m_sst_lev

amrex::Vector< amrex::Vector< amrex::MultiFab * > >	m_tsk_lev

amrex::Vector< amrex::Vector< amrex::iMultiFab * > >	m_lmask_lev

amrex::Vector< amrex::MultiFab * >	m_coupled_sst_lev

amrex::Vector< amrex::iMultiFab * >	m_coupled_sst_valid_lev

amrex::Vector< const amrex::MultiFab * >	m_skin_tsurf_lev

amrex::Vector< const amrex::MultiFab * >	m_skin_moisture_lev

amrex::Vector< int >	m_skin_moisture_resistances_lev

amrex::Vector< amrex::Real >	m_skin_moisture_hs_lev

amrex::Vector< amrex::Vector< amrex::MultiFab * > >	m_lsm_data_lev

amrex::Vector< amrex::Vector< amrex::MultiFab * > >	m_lsm_flux_lev

amrex::Vector< std::string >	m_lsm_data_name

amrex::Vector< std::string >	m_lsm_flux_name

amrex::Vector< amrex::MultiFab * >	m_Hwave_lev

amrex::Vector< amrex::MultiFab * >	m_Lwave_lev

amrex::Vector< amrex::MultiFab * >	m_eddyDiffs_lev

bool	m_update_k_rans = false

amrex::Real	inv_Cmu2 = zero

amrex::Real	theta_ref = zero

Detailed Description

Abstraction layer for different surface layer schemes (e.g. MOST, Cd)

van der Laan, P., Kelly, M. C., & Sørensen, N. N. (2017). A new k-epsilon model consistent with Monin-Obukhov similarity theory. Wind Energy, 20(3), 479–amrex::Real(489.) https://doi.org/amrex::Real(10.1002)/we.2017

Consistent with Dyer (1974) formulation from page 57, Chapter 2, Modeling the vertical ABL structure in Modelling of Atmospheric Flow Fields, Demetri P Lalas and Corrado F Ratto, January 1996, https://doi.org/amrex::Real(10.1142)/amrex::Real(2975.)

Member Enumeration Documentation

◆ FluxCalcType

enum SurfaceLayer::FluxCalcType

strong

Enumerator
MOENG	Moeng functional form.
CUSTOM	Custom constant flux functional form.
BULK_COEFF	Bulk transfer coefficient functional form.
ROTATE	Terrain rotation flux functional form.
RICO

                            {
     MOENG = 0,  ///< Moeng functional form
     CUSTOM,     ///< Custom constant flux functional form
     BULK_COEFF, ///< Bulk transfer coefficient functional form
     ROTATE,     ///< Terrain rotation flux functional form
     RICO
   };

◆ MoistCalcType

enum SurfaceLayer::MoistCalcType

strong

Enumerator
ADIABATIC
MOISTURE_FLUX	Qv-flux specified.
SURFACE_MOISTURE	Surface Qv specified.

                             {
     ADIABATIC = 0,
     MOISTURE_FLUX,      ///< Qv-flux specified
     SURFACE_MOISTURE    ///< Surface Qv specified
   };

◆ PBLHeightCalcType

enum SurfaceLayer::PBLHeightCalcType

strong

Enumerator
None
MYNN25
YSU
MRF
YSUNew

1783 { None, MYNN25, YSU, MRF, YSUNew };

LSM_Type::None

@ None

Definition: ERF_LSM_Types.H:6

◆ RoughCalcType

enum SurfaceLayer::RoughCalcType

strong

Enumerator
CONSTANT	Constant z0.
CHARNOCK
MODIFIED_CHARNOCK
DONELAN
WAVE_COUPLED

                             {
     CONSTANT = 0, ///< Constant z0
     CHARNOCK,
     MODIFIED_CHARNOCK,
     DONELAN,
     WAVE_COUPLED
   };

◆ ThetaCalcType

enum SurfaceLayer::ThetaCalcType

strong

Enumerator
ADIABATIC
HEAT_FLUX	Heat-flux specified.
SURFACE_TEMPERATURE	Surface temperature specified.

                             {
     ADIABATIC = 0,
     HEAT_FLUX,          ///< Heat-flux specified
     SURFACE_TEMPERATURE ///< Surface temperature specified
   };

Constructor & Destructor Documentation

◆ SurfaceLayer()

SurfaceLayer::SurfaceLayer	(	amrex::Orientation	face,
		const amrex::Vector< amrex::Geometry > &	geom,
		bool &	use_rot_surface_flux,
		std::string	a_pp_prefix,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	Qv_prim,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	z_phys_nd,
		const amrex::Vector< amrex::Vector< amrex::Real >> &	zlevels_stag,
		const MeshType &	a_mesh_type,
		const TerrainType &	a_terrain_type,
		const TurbChoice &	a_turb_choice,
		amrex::Real	a_rdOcp,
		double	start_low_time,
		double	final_low_time,
		double	low_time_interval = `0.0`,
		const amrex::Vector< const eb_ * > &	eb_vec = `{}`,
		SurfaceModel *	surf_model = `nullptr`,
		bool	a_has_immersed_buildings = `false`
	)

inlineexplicit

Construct the surface-layer interface.

Parameters

[in]	face	orientation of face (for wall geometries)
[in]	geom	geometry for all AMR levels
[in,out]	use_rot_surface_flux	whether to use rotated surface fluxes
[in]	a_pp_prefix	ParmParse prefix used by MOST averages
[in]	Qv_prim	primitive water-vapor fields by level
[in]	z_phys_nd	nodal physical-height fields by level
[in]	zlevels_stag	nominal staggered z levels by level
[in]	a_mesh_type	mesh type
[in]	a_terrain_type	terrain representation
[in]	a_turb_choice	turbulence-model options
[in]	a_rdOcp	configured Rd/cp exponent for T-theta conversions
[in]	start_low_time	first available low-boundary-data time
[in]	final_low_time	final available low-boundary-data time
[in]	low_time_interval	low-boundary-data time interval
[in]	eb_vec	optional embedded-boundary geometry data
[in]	surf_model	optional weighted surface-model interface
[in]	a_has_immersed_buildings	whether building immersed forcing is configured

                                                                  {},
                          SurfaceModel* surf_model = nullptr,
                          bool a_has_immersed_buildings = false)
   : m_face(face),
     m_pp_prefix(a_pp_prefix),
     m_geom(geom),
     m_rotate(use_rot_surface_flux),
     m_start_low_time(start_low_time),
     m_final_low_time(final_low_time),
     m_low_time_interval(low_time_interval),
     m_eb_vec(eb_vec),
     m_terrain_type(a_terrain_type),
     m_has_immersed_buildings(a_has_immersed_buildings),
     m_rdOcp(a_rdOcp),
     m_ma(face, geom, (z_phys_nd[0] != nullptr), a_pp_prefix, a_mesh_type, a_terrain_type,
          zlevels_stag, eb_vec),
     m_surf_model(surf_model)
   {
       // We have a moisture model if Qv_prim is a valid pointer
       use_moisture = (Qv_prim[0].get());
  
       // Keep standalone SurfaceLayer instances compatible with the historical
       // single-face behavior. ERF replaces this with the complete active-face
       // set after constructing all SurfaceLayer objects.
       m_surface_layer_faces[static_cast<int>(face)] = 1;
  
       // Get roughness
       amrex::ParmParse pp(a_pp_prefix);
       m_z0_given = pp.contains("most.z0");
       pp.queryAdd("most.z0", z0_const);
  
       // Specify how to compute the flux
       if (use_rot_surface_flux) {
           flux_type = FluxCalcType::ROTATE;
       } else {
           std::string flux_string_in;
           std::string flux_string{"moeng"};
           auto read_flux = pp.queryAdd("surface_layer.flux_type", flux_string_in);
           if (read_flux) {
               flux_string = amrex::toLower(flux_string_in);
           }
           if (flux_string == "moeng") {
               flux_type = FluxCalcType::MOENG;
           } else if (flux_string == "rico") {
               flux_type = FluxCalcType::RICO;
           } else if (flux_string == "bulk_coeff") {
               flux_type = FluxCalcType::BULK_COEFF;
           } else if (flux_string == "custom") {
               flux_type = FluxCalcType::CUSTOM;
           } else {
               amrex::Abort("Undefined MOST flux type!");
           }
       }
  
       if (flux_type == FluxCalcType::BULK_COEFF ||
           flux_type == FluxCalcType::CUSTOM ||
           flux_type == FluxCalcType::RICO) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow(),
               "BULK_COEFF, CUSTOM, and RICO surface-layer fluxes are supported only on the z-low face.");
       }
  
       // Include w* to handle free convection (Beljaars 1995, QJRMS)
       pp.queryAdd("most.include_wstar", m_include_wstar);
  
       std::string pblh_string_in;
       std::string pblh_string{"none"};
       auto read_pblh = pp.queryAdd("most.pblh_calc", pblh_string_in);
       if (read_pblh) {
           pblh_string = amrex::toLower(pblh_string_in);
       }
       if (pblh_string == "none") {
           pblh_type = PBLHeightCalcType::None;
       } else if (pblh_string == "mynn25") {
           pblh_type = PBLHeightCalcType::MYNN25;
       } else if (pblh_string == "mynnedmf") {
           pblh_type = PBLHeightCalcType::MYNN25;
       } else if (pblh_string == "ysu") {
           pblh_type = PBLHeightCalcType::YSU;
       } else if (pblh_string == "mrf") {
           pblh_type = PBLHeightCalcType::MRF;
       } else {
           amrex::Abort("Undefined PBLH calc type!");
       }
  
       if (pblh_type != PBLHeightCalcType::None || m_include_wstar) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow(),
               "MOST PBL-height calculation and wstar correction are supported only on the z-low face.");
       }
  
       // The w* correction is computed from the PBL height, so it needs a scheme that
       // actually diagnoses one.  With pblh_calc = "none" the pblh MultiFab keeps the
       // bogus_large_value it was initialized with and calc_wstar turns that into a
       // convective velocity scale of ~1e50, which destroys the surface fluxes.
       if (m_include_wstar && (pblh_type == PBLHeightCalcType::None)) {
           amrex::Abort("erf.most.include_wstar requires a PBL height: set "
                        "erf.most.pblh_calc (MYNN25 is the only scheme implemented)");
       }
  
       if (m_surf_model) {
           use_surface_model = true;
           amrex::Print() << " Using fluxes from surface model!" << std::endl;
           surf_model_fluxes = m_surf_model->are_fluxes();
       }
  
       // Get surface temperature, as a potential temperature (most.pot_surf_temp) or an
       // absolute temperature (most.abs_surf_temp); most.surf_temp is the deprecated
       // spelling of most.pot_surf_temp.  These are read with query, not queryAdd, so
       // that whether a key was given is the query return value: queryAdd would add a
       // key that is absent, and these keys are also read by ERF_InputSoundingData.H and
       // by the SurfaceLayer of every face, so the first reader would make the key look
       // given to the rest.  A key that is given must be positive; it is never ignored.
       {
           amrex::Real pot_st = amrex::Real(-1.);
           amrex::Real abs_st = amrex::Real(-1.);
           amrex::Real old_st = amrex::Real(-1.);
           const bool has_pot = pp.query("most.pot_surf_temp", pot_st);
           const bool has_abs = pp.query("most.abs_surf_temp", abs_st);
           const bool has_old = pp.query("most.surf_temp",     old_st);
           if (int(has_pot) + int(has_abs) + int(has_old) > 1) {
               amrex::Abort("Specify only one of most.pot_surf_temp, most.abs_surf_temp "
                            "and the deprecated most.surf_temp");
           }
           if ((has_pot && !(pot_st > amrex::Real(0))) ||
               (has_abs && !(abs_st > amrex::Real(0))) ||
               (has_old && !(old_st > amrex::Real(0)))) {
               amrex::Abort("most.pot_surf_temp, most.abs_surf_temp and most.surf_temp "
                            "must be positive temperatures [K]");
           }
           if (has_old) {
               amrex::Print() << "WARNING: most.surf_temp is deprecated and is read as a surface "
                              << "potential temperature; use most.pot_surf_temp (potential) or "
                              << "most.abs_surf_temp (absolute) instead" << std::endl;
               pot_st = old_st;
           }
           if (has_abs) {
               surf_temp = abs_st;
               surf_temp_is_absolute = true;
           } else {
               surf_temp = pot_st;
           }
       }
       const bool erf_st = (surf_temp > amrex::Real(0));
       // With an absolute temperature this is only the placeholder until update_fluxes
       // converts it with the surface pressure.
       if (erf_st) { default_land_surf_temp  = surf_temp; }
  
       if (surf_temp_is_absolute) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow() && m_terrain_type != TerrainType::EB,
               "most.abs_surf_temp is supported only on the z-low face without EB terrain");
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               flux_type != FluxCalcType::CUSTOM && flux_type != FluxCalcType::RICO &&
               flux_type != FluxCalcType::BULK_COEFF,
               "most.abs_surf_temp is not supported with custom, rico or bulk_coeff fluxes; "
               "use most.pot_surf_temp");
       }
  
       // Get surface moisture
       bool erf_sq = false;
       if (use_moisture) {
           pp.queryAdd("most.surf_moist", surf_moist);
           erf_sq = (surf_moist >= amrex::Real(0));
       }
       if (erf_sq) { default_land_surf_moist = surf_moist; }
  
       // Custom type user must specify the fluxes
       if (flux_type == FluxCalcType::CUSTOM ||
           flux_type == FluxCalcType::RICO) {
           theta_type = ThetaCalcType::HEAT_FLUX;
           moist_type = MoistCalcType::MOISTURE_FLUX;
           pp.get("most.ustar", custom_ustar);
           pp.get("most.tstar", custom_tstar);
           pp.get("most.qstar", custom_qstar);
           pp.queryAdd("most.rhosurf", custom_rhosurf);
           if (custom_qstar != 0) {
               AMREX_ASSERT_WITH_MESSAGE(use_moisture,
               "Specified custom MOST qv flux without moisture model!");
           }
           amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
                          << custom_ustar << " " << custom_tstar << " "
                          << custom_qstar << std::endl;
  
       // Bulk transfer coefficient (must specify coeffs and surface values)
       } else if (flux_type == FluxCalcType::BULK_COEFF) {
           pp.get("most.Cd", m_Cd);
           pp.get("most.Ch", m_Ch);
           pp.get("most.Cq", m_Cq);
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(erf_st,
               "bulk_coeff fluxes need most.pot_surf_temp");
           default_land_surf_temp = surf_temp;
           pp.get("most.surf_moist", default_land_surf_moist);
           amrex::Print() << "Using specified Cd, Ch, Cq for MOST = "
                          << m_Cd << " " << m_Ch << " "
                          << m_Cq << std::endl;
  
       // Specify surface temperature/moisture or surface flux
       } else {
           if (erf_st) {
               theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
               pp.queryAdd("most.surf_heating_rate", surf_heating_rate); // [K/h]
  
               // Modify rate to be in units of K / s rather than K / hr
               surf_heating_rate /= amrex::Real(3600.0);        // [K/s]
  
               if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
                   amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
               }
           } else {
               pp.queryAdd("most.surf_temp_flux", surf_temp_flux);
  
               if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
                   amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
               }
               if (std::abs(surf_temp_flux) >
                   std::numeric_limits<amrex::Real>::epsilon()) {
                   theta_type = ThetaCalcType::HEAT_FLUX;
               } else {
                   theta_type = ThetaCalcType::ADIABATIC;
               }
           }
  
           if (erf_sq) {
               moist_type = MoistCalcType::SURFACE_MOISTURE;
           } else {
               pp.queryAdd("most.surf_moist_flux", surf_moist_flux);
               if (std::abs(surf_moist_flux) >
                   std::numeric_limits<amrex::Real>::epsilon()) {
                   moist_type = MoistCalcType::MOISTURE_FLUX;
               } else {
                   moist_type = MoistCalcType::ADIABATIC;
               }
           }
       }
  
       if (theta_type == ThetaCalcType::HEAT_FLUX ||
           (theta_type == ThetaCalcType::ADIABATIC &&
            moist_type == MoistCalcType::ADIABATIC)) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow(),
               "HEAT_FLUX and ADIABATIC surface-layer fluxes are supported only on the z-low face.");
       }
  
       if (flux_type == FluxCalcType::RICO)
       {
           pp.queryAdd("most.rico.theta_z0", rico_theta_z0);
           pp.queryAdd("most.rico.qsat_z0", rico_qsat_z0);
       }
  
       // Make sure the inputs file doesn't try to use most.roughness_type
       std::string bogus_input;
       if (pp.queryAdd("most.roughness_type", bogus_input) > 0) {
           amrex::Abort("most.roughness_type is deprecated; use "
                        "most.roughness_type_land and/or most.roughness_type_sea");
       }
  
       // Specify how to compute the surface flux over land (if there is any)
       std::string rough_land_string_in;
       std::string rough_land_string{"constant"};
       auto read_rough_land =
           pp.queryAdd("most.roughness_type_land", rough_land_string_in);
       if (read_rough_land) {
           rough_land_string = amrex::toLower(rough_land_string_in);
       }
       if (rough_land_string == "constant") {
           rough_type_land = RoughCalcType::CONSTANT;
       } else {
           amrex::Abort("Undefined MOST roughness type for land!");
       }
  
       // Allow for smooth-flow limit?
       pp.queryAdd("most.smooth_flow_viscosity", smooth_flow_visc);
       amrex::Print() << "The smooth-flow limit will be included for variable roughness models over sea: " << smooth_flow_visc << "\n";
  
       // Specify how to compute the surface flux over sea (if there is any)
       std::string rough_sea_string_in;
       std::string rough_sea_string{"charnock"};
       auto read_rough_sea = pp.queryAdd("most.roughness_type_sea", rough_sea_string_in);
       if (read_rough_sea) {
           rough_sea_string = amrex::toLower(rough_sea_string_in);
       }
       if (rough_sea_string == "charnock") {
           rough_type_sea = RoughCalcType::CHARNOCK;
           pp.queryAdd("most.charnock_constant", cnk_a);
           if (cnk_a > 0) {
               amrex::Print() << "If there is water, Charnock relation with C_a="
                              << cnk_a << " will be used" << std::endl;
           } else {
               amrex::Print() << "If there is water, Charnock relation with variable "
                   "Charnock parameter (COARE3.0) will be used" << std::endl;
           }
       } else if (rough_sea_string == "coare3.0") {
           rough_type_sea = RoughCalcType::CHARNOCK;
           amrex::Print() << "If there is water, Charnock relation with variable "
               "Charnock parameter (COARE3.0) will be used" << std::endl;
           cnk_a = -1;
       } else if (rough_sea_string == "donelan") {
           rough_type_sea = RoughCalcType::DONELAN;
       } else if (rough_sea_string == "modified_charnock") {
           rough_type_sea = RoughCalcType::MODIFIED_CHARNOCK;
           pp.queryAdd("most.modified_charnock_depth", depth);
           if (depth < amrex::Real(10.0) || depth > amrex::Real(100.0) ) {
               amrex::Print() << "Specified depth of " << depth
                              << " is outside the valid range of [10 100], resetting to bounds now."
                              << std::endl;
           }
           // Limiter based upon the fit range in Jiménez & Dudhia
           depth = amrex::min(amrex::max(depth,amrex::Real(10.0)),amrex::Real(100.0));
       } else if (rough_sea_string == "wave_coupled") {
           rough_type_sea = RoughCalcType::WAVE_COUPLED;
       } else if (rough_sea_string == "constant") {
           rough_type_sea = RoughCalcType::CONSTANT;
       } else {
           amrex::Abort("Undefined MOST roughness type for sea!");
       }
  
       // use skin temperature instead of sea-surface temperature
       // (wrfinput data may have lower resolution SST data)
       pp.queryAdd("most.ignore_sst", m_ignore_sst);
  
       // If we're using the RANS k model, then we need to update the dirichlet
       // wall value of k (written into the first cell above the wall and held
       // there through the step) based on the instantaneous u* and θ*; the turbulence modeling
       // choices can vary per level but for now, assume that if specified then
       // all levels are using the same RANS model.
       m_update_k_rans = (a_turb_choice.rans_type == RANSType::kEqn &&
                          a_turb_choice.dirichlet_k == true);
       if (m_update_k_rans) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow(),
               "RANS surface-layer k updates are supported only on the z-low face.");
       }
       if (m_update_k_rans) {
           inv_Cmu2 = one / (a_turb_choice.Cmu0 * a_turb_choice.Cmu0);
           theta_ref = a_turb_choice.theta_ref;
       }
  
   } // constructor

Member Function Documentation

◆ compute_fluxes() [1/2]

template<typename FluxIter >

void SurfaceLayer::compute_fluxes	(	const int &	lev,
		const int &	max_iters,
		amrex::MultiFab &	cons_in,
		const FluxIter &	most_flux,
		bool	is_land
	)

Compute MOST fluxes with a selected flux-iteration functor.

Parameters

[in]	lev	level index
[in]	max_iters	maximum MOST iteration count
[in,out]	cons_in	conserved state used by the flux computation
[in]	most_flux	flux-iteration functor
[in]	is_land	whether the land-surface branch is active

◆ compute_fluxes() [2/2]

template<typename FluxIter >

void SurfaceLayer::compute_fluxes	(	const int &	lev,
		const int &	max_iters,
		MultiFab &	,
		const FluxIter &	most_flux,
		bool	is_land
	)

Function to compute the fluxes (u^star and t^star) for Monin Obukhov similarity theory

Parameters

[in]	lev	Current level
[in]	max_iters	Maximum iterations to use
[in]	cons_in	Conserved state whose grids define the surface iteration
[in]	most_flux	Flux-iteration functor used to compute ustar, tstar, qstar, and related fields
[in]	is_land	Selects whether land or sea cells are updated

 {
     // Pointers to the computed averages
     const auto *const tm_ptr      = m_ma.get_average(lev,3); // potential temperature
     const auto *const qvm_ptr     = m_ma.get_average(lev,4); // water vapor mixing ratio
     const auto *const tvm_ptr     = m_ma.get_average(lev,5); // virtual potential temperature
     const auto *const umm_ptr     = m_ma.get_average(lev,6); // horizontal velocity magnitude
     const auto *const uw_mag_mean = m_ma.get_average(lev,7); // x/z velocity magnitude
     const auto *const vw_mag_mean = m_ma.get_average(lev,8); // y/z velocity magnitude
     const auto *const zref_ptr    = m_ma.get_zref(lev);      // reference height
     const bool l_use_eb = (m_terrain_type == TerrainType::EB);
  
     const int dir = m_face.coordDir();
     // The source fields and Taus use the same BoxArray ordering and
     // DistributionMapping. Full-state fields follow cons_in, while the
     // collapsed SurfaceLayer/MOSTAverage fields follow u_star.
     int sm_index = 0;
     if (m_face.isLow()) {
         sm_index = m_geom[lev].Domain().smallEnd(dir);
     } else {
         sm_index = m_geom[lev].Domain().bigEnd(dir);
     }
  
     // Use the 2-D surface mask as the iterator so ranks participate only when
     // their grids coincide with the selected face.
     for (MFIter mfi(*m_lmask_lev[lev][0]); mfi.isValid(); ++mfi)
     {
         // Face ownership is a property of the full grid, not an individual
         // tile. Keep the tile box for the kernel below, but test the valid
         // box so all tiles of a face-owning grid are processed.
         Box vbx = mfi.validbox();
         const Box tbx = mfi.tilebox();
         Box gtbx = mfi.growntilebox();
  
         // Since lmask is used in the MFIter, its z extent is collapsed.  The
         // lateral kernels need the full z column; z-face ownership is instead
         // determined from the original 3-D surface-copy mapping.
         if (dir != 2 || l_use_eb) {
             gtbx.setSmall(2, m_geom[lev].Domain().smallEnd(2));
             gtbx.setBig(2, m_geom[lev].Domain().bigEnd(2));
         }
  
         const bool owns_surface = l_use_eb || dir != 2 ||
             m_planar_bndry[lev].is_surface_copy(mfi.index());
         if (!owns_surface) {
             continue;
         }
  
         if (dir != 2) {
             if (m_face.isLow()) {
                 if (vbx.smallEnd(dir) != sm_index ||
                     tbx.smallEnd(dir) != sm_index) {
                     continue;
                 }
             } else {
                 if (vbx.bigEnd(dir) != sm_index ||
                     tbx.bigEnd(dir) != sm_index) {
                     continue;
                 }
             }
         }
  
         if (m_face.isLow()) {
             gtbx.setSmall(dir, sm_index);
             gtbx.setBig(dir, sm_index);
         } else {
             gtbx.setSmall(dir, sm_index);
             gtbx.setBig(dir, sm_index);
         }
  
         // X/Y faces still need the full valid z column
         if (!l_use_eb && dir == 2) {
             gtbx.makeSlab(2, sm_index);
         } else {
             gtbx.setBig(2, m_geom[lev].Domain().bigEnd(2));
         }
  
         // The mask iterator can have more lateral ghost cells than the
         // surface-layer fields, and an EB FAB may be decomposed in z. Keep
         // the face selection above, but never launch outside the target FAB.
         gtbx &= u_star[lev]->fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         auto u_star_arr = u_star[lev]->array(mfi);
         auto t_star_arr = t_star[lev]->array(mfi);
         auto q_star_arr = q_star[lev]->array(mfi);
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto q_surf_arr = q_surf[lev]->array(mfi);
         auto olen_arr   = olen[lev]->array(mfi);
  
         const auto tm_arr   = tm_ptr->array(mfi);
         const auto tvm_arr  = tvm_ptr->array(mfi);
         const auto qvm_arr  = qvm_ptr->array(mfi);
         const auto umm_arr  = umm_ptr->array(mfi);
         const auto vwmm_arr = (dir == 0) ? vw_mag_mean->array(mfi) : Array4<Real>{};
         const auto uwmm_arr = (dir == 1) ? uw_mag_mean->array(mfi) : Array4<Real>{};
         const auto zref_arr = zref_ptr->array(mfi);
  
         // umm depending on face direction (YZ, XZ, XY)
         const auto dir_umm_arr = ((dir == 0) ? vwmm_arr : ((dir == 1) ? uwmm_arr : umm_arr));
         const auto z0_arr   = z_0[lev].array(mfi);
  
         // PBL height if we need to calculate wstar for the Beljaars correction
         // TODO: can/should we apply this in LES mode?
         const auto w_star_arr = (m_include_wstar) ? w_star[lev].get()->array(mfi) : Array4<Real> {};
         const auto pblh_arr   = (m_include_wstar) ? pblh[lev].get()->array(mfi) : Array4<Real> {};
  
         // Wave properties if they exist
         const auto Hwave_arr = (m_Hwave_lev[lev]) ? m_Hwave_lev[lev]->array(mfi) : Array4<Real> {};
         const auto Lwave_arr = (m_Lwave_lev[lev]) ? m_Lwave_lev[lev]->array(mfi) : Array4<Real> {};
         const auto eta_arr   = (m_eddyDiffs_lev[lev]) ? m_eddyDiffs_lev[lev]->array(mfi) : Array4<Real> {};
  
         // Land mask array if it exists
         auto lmask_arr    = (m_lmask_lev[lev][0])    ? m_lmask_lev[lev][0]->array(mfi) :
                                                        Array4<int> {};
  
         // Get EB flags if needed
         const auto flag_arr = (l_use_eb) ? m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab()[mfi].const_array() : Array4<const EBCellFlag>{};
  
         if (!l_use_eb) {
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 // always check for land mask at k=0 even if on other faces
                 if (( is_land && lmask_arr(i,j,0) == 1) ||
                     (!is_land && lmask_arr(i,j,0) == 0))
                 {
                     // NOTE: All 2D MFs so k index is always 0 from ba2d definition
                     most_flux.iterate_flux(i, j, k, max_iters,
                                         zref_arr,                            // set in most average
                                         z0_arr,                              // updated if(!is_land)
                                         dir_umm_arr, tm_arr, tvm_arr, qvm_arr,
                                         u_star_arr,                          // updated
                                         w_star_arr,                          // updated if(m_include_wstar)
                                         t_star_arr, q_star_arr,              // updated
                                         t_surf_arr, q_surf_arr, olen_arr,    // updated
                                         pblh_arr,                            // updated if(m_include_wstar)
                                         Hwave_arr, Lwave_arr, eta_arr);
                 }
             });
         // EB
         } else {
             if (std::is_same<FluxIter, adiabatic_eb>::value ||
                 std::is_same<FluxIter, surface_temp_eb>::value ||
                 std::is_same<FluxIter, surface_flux_eb>::value) {
                 ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
                     if (( is_land && lmask_arr(i,j,0) == 1) ||
                         (!is_land && lmask_arr(i,j,0) == 0))
                     {
                         if (flag_arr(i,j,k).isSingleValued()) {
                             most_flux.iterate_flux(i, j, k, max_iters,
                                                 zref_arr,                            // set in most average
                                                 z0_arr,                              // updated if(!is_land)
                                                 umm_arr, tm_arr, tvm_arr, qvm_arr,
                                                 u_star_arr,                          // updated
                                                 w_star_arr,                          // updated if(m_include_wstar)
                                                 t_star_arr, q_star_arr,              // updated
                                                 t_surf_arr, q_surf_arr, olen_arr,    // updated
                                                 pblh_arr,                            // updated if(m_include_wstar)
                                                 Hwave_arr, Lwave_arr, eta_arr);
                         }
                     }
                 });
             } else {
                 amrex::Abort("FluxIter type not supported for EB");
             }
         }
     }
 }

Here is the call graph for this function:

◆ compute_pblh() [1/2]

template<typename PBLHeightEstimator >

void SurfaceLayer::compute_pblh	(	const int &	lev,
		amrex::Vector< amrex::Vector< amrex::MultiFab >> &	vars,
		amrex::MultiFab *	z_phys_cc,
		const PBLHeightEstimator &	est,
		const MoistureComponentIndices &	moisture_indice
	)

Compute planetary-boundary-layer height with the selected estimator.

Parameters

[in]	lev	level index
[in,out]	vars	state variables used by the PBL-height calculation
[in]	z_phys_cc	cell-centered physical-height field
[in]	est	PBL-height estimator functor
[in]	moisture_indice	indices for moisture components

◆ compute_pblh() [2/2]

template<typename PBLHeightEstimator >

void SurfaceLayer::compute_pblh	(	const int &	lev,
		Vector< Vector< MultiFab >> &	vars,
		MultiFab *	z_phys_cc,
		const PBLHeightEstimator &	est,
		const MoistureComponentIndices &	moisture_indices
	)

Compute PBL height with the supplied estimator.

Parameters

[in]	lev	Current level
[in]	vars	Level-indexed state MultiFabs passed to the estimator
[in]	z_phys_cc	Cell-centered physical height used by the estimator
[in]	est	PBL height estimator functor
[in]	moisture_indices	Moisture component indices used by the estimator

 {
     const MultiFab& cons = vars[lev][Vars::cons];
     const iMultiFab* lmask = m_lmask_lev[lev][0];
  
     // The estimator scans each box from its lowest cell to its highest and writes the planar
     // pblh of that box, so every box it is given must start at the ground.  Grids that hold
     // such boxes only -- full height or not -- go straight to it.  Any other grids (boxes
     // stacked in z, or boxes aloft) go through columns: the runs of cells that start at the
     // ground, each as one box (see define_pblh_columns).
     if (static_cast<int>(m_pblh_columns.size()) <= lev) { m_pblh_columns.resize(lev+1); }
     if (m_pblh_columns[lev].ba != cons.boxArray() ||
         m_pblh_columns[lev].dm != cons.DistributionMap()) {
         define_pblh_columns(lev, cons.boxArray(), cons.DistributionMap());
     }
     const PBLHColumns& cols = m_pblh_columns[lev];
  
     // A refined level whose grids do not run from the ground to the top of the domain in every
     // column cannot diagnose the boundary layer there: a level that ends below the top would
     // return a height capped by its grids, and a level that does not reach the ground has no
     // column to scan (it was left at zero).  Such a level takes the PBL height of the next
     // coarser level at each of its columns.  Levels are updated coarse to fine, so that height
     // is the one set at the start of the enclosing coarse step: with subcycling, the later fine
     // substeps read a coarse height up to one coarse step old.  A refined level that spans the
     // full height diagnoses its own.
     if (lev > 0 && !cols.full_height && m_face.coordDir() == 2 && m_face.isLow() &&
         m_terrain_type != TerrainType::EB) {
         fill_pblh_from_coarser(lev);
         return;
     }
  
     if (!cols.needed) {
         est.compute_pblh(m_geom[lev], z_phys_cc, pblh[lev].get(), cons, lmask, moisture_indices);
         return;
     }
  
     // Zero is the estimator's own value for a height it did not find.  It stays on the planar
     // boxes over which no box of this level reaches the ground: all of them on a level that
     // lies entirely aloft.
     pblh[lev]->setVal(zero);
     if (cols.ba_col.empty()) { return; }
  
     const Periodicity period = m_geom[lev].periodicity();
  
     // The estimator reads the density, the potential temperature, the TKE and the moisture
     // species that enter theta_v.  Those live in [0, RhoKE_comp] and in the moist window, so
     // the columns carry the state up to the highest of them and no further: the species above
     // it (the number concentrations of a two-moment scheme, the non-water species) are the
     // bulk of a moist state and are never read here.  The span is contiguous rather than the
     // two pieces it is made of, so that every component of cons_col is filled by the copies
     // below -- a component the estimator reads must never be one this routine left unset --
     // and the components keep their place, since the estimator indexes the state by
     // component number.
     int q_hi = -1;
     for (const int q : {moisture_indices.qv, moisture_indices.qc, moisture_indices.qi,
                         moisture_indices.qr, moisture_indices.qs, moisture_indices.qg}) {
         AMREX_ALWAYS_ASSERT((q < 0) || ((q > RhoKE_comp) && (q < cons.nComp())));
         q_hi = std::max(q_hi, q);
     }
     const int ncomp_col = std::max(RhoKE_comp+1, q_hi+1);
  
     // The halo the columns need: in x and y the ghost cells of pblh, which the estimator
     // writes and so reads the state over, and in z the one cell above the top of each column
     // that the scan's k+1 reads end in.  The state must hold that halo for the copies below
     // to have anything to take it from.
     const IntVect ng_pblh = pblh[lev]->nGrowVect();
     const IntVect ng_col = elemwiseMax(ng_pblh, IntVect(AMREX_D_DECL(0,0,1)));
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(cons.nGrowVect().allGE(ng_col),
         "erf.most.pblh_calc = MYNN25 on grids that do not all start at the ground needs the "
         "state to carry at least the ghost cells of the surface-layer fields, and one in z.");
  
     // Every valid cell of a column is a valid cell of this level, and every ghost cell of a
     // column is a valid or a ghost cell of the box that holds the cell next to it, so the two
     // passes below leave no cell of the columns unset.  Ghost cells go first (they hold the
     // physical boundary values and, next to a coarser level, the values interpolated from
     // it), then the valid cells, so that every cell this level owns comes from the box that
     // owns it and not from a neighbour's ghost cell.
     MultiFab cons_col(cols.ba_col, cols.dm_col, ncomp_col, ng_col);
     for (const IntVect& ng_src : {ng_col, IntVect(0)}) {
         cons_col.ParallelCopy(cons, 0, 0, ncomp_col, ng_src, ng_col, period);
     }
  
     std::unique_ptr<MultiFab> zcc_col;
     if (z_phys_cc) {
         const IntVect ng_z = z_phys_cc->nGrowVect();
         zcc_col = std::make_unique<MultiFab>(cols.ba_col, cols.dm_col, 1, ng_z);
         zcc_col->ParallelCopy(*z_phys_cc, 0, 0, 1, ng_z, ng_z, period);
         zcc_col->ParallelCopy(*z_phys_cc, 0, 0, 1, IntVect(0), ng_z, period);
     }
  
     std::unique_ptr<iMultiFab> lmask_col;
     if (lmask) {
         lmask_col = std::make_unique<iMultiFab>(cols.ba_col2d, cols.dm_col, 1, ng_pblh);
         lmask_col->setVal(1);
         lmask_col->ParallelCopy(*lmask, 0, 0, 1, elemwiseMin(lmask->nGrowVect(), ng_pblh), ng_pblh, period);
         lmask_col->ParallelCopy(*lmask, 0, 0, 1, IntVect(0), ng_pblh, period);
     }
  
     MultiFab pblh_col(cols.ba_col2d, cols.dm_col, 1, ng_pblh);
     est.compute_pblh(m_geom[lev], zcc_col.get(), &pblh_col, cons_col, lmask_col.get(), moisture_indices);
  
     // Onto every planar box.  The ghost cells of the columns go first, for the ghost cells of
     // pblh outside the domain.  They also reach valid cells of pblh over which no box of this
     // level starts at the ground (the estimator fills the ghost cells of a column next to such
     // a gap from ghost data), so those are set back to zero before the valid cells of the
     // columns are copied.  Every planar copy of a cell then holds the same value, which makes
     // the FillBoundary that ends this well defined despite the duplicate boxes.
     pblh[lev]->ParallelCopy(pblh_col, 0, 0, 1, ng_pblh, ng_pblh, period);
     pblh[lev]->setVal(zero, 0, 1, 0);
     pblh[lev]->ParallelCopy(pblh_col, 0, 0, 1, IntVect(0), IntVect(0), period);
     pblh[lev]->FillBoundary(period);
 }

Here is the call graph for this function:

◆ compute_sfc_params_from_lsm_fluxes()

void SurfaceLayer::compute_sfc_params_from_lsm_fluxes	(	const int &	lev,
		amrex::MultiFab &	cons_in
	)

Derive MOST surface parameters from LSM fluxes.

Parameters

[in]	lev	level index
[in,out]	cons_in	conserved state used by the surface-parameter computation

Compute surface-layer parameters from land-surface-model fluxes.

Parameters

[in]	lev	Current level
[in]	cons_in	Conserved state used to derive density, theta, and moisture at the surface

 {
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         static_cast<int>(m_face) == Orientation::zlo(),
         "LSM surface-layer parameters are supported only on the z-low face.");
  
     Real eps = std::numeric_limits<amrex::Real>::epsilon();
     bool has_moisture = use_moisture;
     const bool use_surface_model_fluxes = use_surface_model && surf_model_fluxes;
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const auto *const umm_ptr  = m_ma.get_average(lev,6); // horizontal velocity magnitude
     const auto *const zref_ptr = m_ma.get_zref(lev);     // reference height
     for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
  
         Box vbx = mfi.validbox();
         if (vbx.smallEnd(2) != klo) { continue; }
         vbx.makeSlab(2,0);
  
         // The ghost columns are set too, so that they hold the same kind of value as the
         // valid columns: compute_fluxes writes MOST values over the grown box, and
         // fill_planar_boundary below refills only the interior and periodic ghosts, so a
         // coarse/fine or domain-edge ghost would otherwise keep a MOST value next to
         // LSM values.  The LSM fluxes do not carry lateral ghost cells, so a ghost column
         // takes the fluxes of the nearest valid column, as in get_lsm_tsurf, with its own
         // atmospheric state.
         const NearestColumn col(vbx);
         const int idx = mfi.index();
         Box gbx = amrex::grow(vbx, IntVect(AMREX_D_DECL(u_star[lev]->nGrow(0),
                                                         u_star[lev]->nGrow(1), 0)));
         Box cons_slab = cons_in.fabbox(idx);
         cons_slab.setRange(2, 0);
         gbx &= u_star[lev]->fabbox(idx);
         gbx &= z_0[lev].fabbox(idx);
         gbx &= umm_ptr->fabbox(idx);
         gbx &= zref_ptr->fabbox(idx);
         gbx &= cons_slab;
         if (m_lmask_lev[lev][0]) { gbx &= m_lmask_lev[lev][0]->fabbox(idx); }
  
         // Get CC state
         const Array4<const Real> cons_arr = cons_in.const_array(mfi);
  
         // Get SL params
         const auto u_star_arr = u_star[lev]->array(mfi);
         const auto t_star_arr = t_star[lev]->array(mfi);
         const auto q_star_arr = q_star[lev]->array(mfi);
         const auto olen_arr   = olen[lev]->array(mfi);
  
         const auto umm_arr  = umm_ptr->array(mfi);
         const auto zref_arr = zref_ptr->array(mfi);
         const auto z0_arr   = z_0[lev].array(mfi);
  
         // Get LSM fluxes
         auto lmask_arr      = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                       Array4<int> {};
         auto lsm_t_flux_arr = Array4<const Real> {};
         auto lsm_q_flux_arr = Array4<const Real> {};
         auto lsm_tau13_arr  = Array4<const Real> {};
         auto lsm_tau23_arr  = Array4<const Real> {};
         // compute_sfc_params_from_lsm_fluxes consumes signed kinematic stress
         // components; their vector magnitude determines u_star^2.
         if (use_surface_model_fluxes) {
             const SurfaceFluxView surf_flux = m_surf_model->get_surface_flux_view(lev, mfi);
             lsm_t_flux_arr = surf_flux.t_flux;
             lsm_q_flux_arr = surf_flux.q_flux;
             lsm_tau13_arr  = surf_flux.tau13;
             lsm_tau23_arr  = surf_flux.tau23;
         } else {
             for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
                 if (toLower(m_lsm_flux_name[n]) == "t_flux") { lsm_t_flux_arr = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                 if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                 if (toLower(m_lsm_flux_name[n]) == "tau13")  { lsm_tau13_arr  = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                 if (toLower(m_lsm_flux_name[n]) == "tau23")  { lsm_tau23_arr  = m_lsm_flux_lev[lev][n]->const_array(mfi); }
             }
         }
  
         ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
         {
             int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
             const int li = col.i(i);
             const int lj = col.j(j);
             // Only ordinary land (mask 1) uses LSM fluxes. Immersed buildings and
             // terrain (mask 2) retain their existing surface-layer parameters.
             // Skip cells where the LSM does not have a valid value (finite and below
             // lsm_undefined) for every one of the four fluxes; these keep the MOST values.
             if (is_land == 1 && lsm_t_flux_arr &&
                 is_valid_lsm_value(lsm_t_flux_arr(li,lj,0)) &&
                 is_valid_lsm_value(lsm_q_flux_arr(li,lj,0)) &&
                 is_valid_lsm_value(lsm_tau13_arr (li,lj,0)) &&
                 is_valid_lsm_value(lsm_tau23_arr (li,lj,0))) {
                 const Real t_flux = lsm_t_flux_arr(li,lj,0);
                 const Real q_flux = lsm_q_flux_arr(li,lj,0);
                 Real rho = cons_arr(i,j,klo,Rho_comp);
                 Real Thd = cons_arr(i,j,klo,RhoTheta_comp) / rho;
                 Real qv  = (has_moisture) ? cons_arr(i,j,klo,RhoQ1_comp) / rho : zero;
                 Real Thv = Thd * (one + epsv*qv);
                 Real tau = std::sqrt( lsm_tau13_arr(li,lj,0)*lsm_tau13_arr(li,lj,0)
                                     + lsm_tau23_arr(li,lj,0)*lsm_tau23_arr(li,lj,0) );
                 u_star_arr(i,j,0) = amrex::max(std::sqrt(tau),eps);
                 if (t_flux>=zero) {
                     t_star_arr(i,j,0) = amrex::min(-t_flux / u_star_arr(i,j,0),-eps);
                 } else {
                     t_star_arr(i,j,0) = amrex::max(-t_flux / u_star_arr(i,j,0),eps);
                 }
                 if (q_flux>=zero) {
                     q_star_arr(i,j,0) = amrex::min(-q_flux / u_star_arr(i,j,0),-eps);
                 } else {
                     q_star_arr(i,j,0) = amrex::max(-q_flux / u_star_arr(i,j,0),eps);
                 }
                 Real tstv = t_star_arr(i,j,0)*(one + epsv*qv) + epsv*Thd*q_star_arr(i,j,0);
                 tstv = (tstv >= zero) ? amrex::max(tstv, eps) : amrex::min(tstv, -eps);
                 olen_arr(i,j,0) = ( u_star_arr(i,j,0) * u_star_arr(i,j,0) * Thv ) /
                                   ( KAPPA * CONST_GRAV * tstv );
                 z0_arr(i,j,0)   = Compute_roughness(zref_arr(i,j,0),   olen_arr(i,j,0),
                                                     umm_arr(i,j,0), u_star_arr(i,j,0));
             }
         });
     } // mfi
 }

Here is the call graph for this function:

◆ compute_SurfaceLayer_bcs() [1/2]

template<typename FluxCalc >

void SurfaceLayer::compute_SurfaceLayer_bcs	(	const int &	lev,
		amrex::Vector< const amrex::MultiFab * >	mfs,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	Tau_lev,
		amrex::MultiFab *	xheat_flux,
		amrex::MultiFab *	yheat_flux,
		amrex::MultiFab *	zheat_flux,
		amrex::MultiFab *	xqv_flux,
		amrex::MultiFab *	yqv_flux,
		amrex::MultiFab *	zqv_flux,
		const amrex::MultiFab *	z_phys,
		const FluxCalc &	flux_comp,
		amrex::MultiFab *	mf_surface_fluxes = `nullptr`
	)

Compute planar-terrain surface-layer flux boundary conditions.

Parameters

[in]	lev	level index
[in]	mfs	state and velocity fields used by the BC computation
[in,out]	Tau_lev	stress fields to fill
[in,out]	xheat_flux	x-face heat flux field
[in,out]	yheat_flux	y-face heat flux field
[in,out]	zheat_flux	z-face heat flux field
[in,out]	xqv_flux	x-face moisture flux field
[in,out]	yqv_flux	y-face moisture flux field
[in,out]	zqv_flux	z-face moisture flux field
[in]	z_phys	physical-height field
[in]	flux_comp	flux-computation functor

◆ compute_SurfaceLayer_bcs() [2/2]

template<typename FluxCalc >

void SurfaceLayer::compute_SurfaceLayer_bcs	(	const int &	lev,
		Vector< const MultiFab * >	mfs,
		Vector< std::unique_ptr< MultiFab >> &	Tau_lev,
		MultiFab *	xheat_flux,
		MultiFab *	yheat_flux,
		MultiFab *	zheat_flux,
		MultiFab *	xqv_flux,
		MultiFab *	yqv_flux,
		MultiFab *	zqv_flux,
		const MultiFab *	z_phys,
		const FluxCalc &	flux_comp,
		MultiFab *	mf_surface_fluxes
	)

Function to calculate MOST fluxes for populating ghost cells.

Parameters

[in]	lev	Current level
[in]	mfs	State MultiFabs used to compute the boundary fluxes
[in,out]	Tau_lev	Diffusive stress MultiFabs populated with surface stresses
[in,out]	xheat_flux	x-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]	yheat_flux	y-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]	zheat_flux	z-face heat-flux MultiFab populated with vertical surface heat flux
[in,out]	xqv_flux	x-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]	yqv_flux	y-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]	zqv_flux	z-face moisture-flux MultiFab populated when moisture is enabled
[in]	z_phys	Nodal physical height used to rotate terrain-following fluxes
[in]	flux_comp	Flux-calculation functor used to compute scalar and momentum fluxes

 {
     bool rotate = m_rotate;
  
     const int dir = m_face.coordDir();
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         !mfs.empty() && mfs[0] != nullptr,
         "Surface-layer stress computation requires a conserved-state MultiFab.");
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         Tau_lev[TauType::tau13] != nullptr && Tau_lev[TauType::tau23] != nullptr,
         "tau13 and tau23 are required by surface-layer stress computation.");
     int sm_index = 0;
     if (m_face.isLow()) {
         sm_index = m_geom[lev].Domain().smallEnd(dir);
     } else {
         sm_index = m_geom[lev].Domain().bigEnd(dir);
     }
  
     // These requirements are invariant over the loop below. Check them once
     // before acquiring per-tile array views.
     if (dir == 0) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
             Tau_lev[TauType::tau31] != nullptr,
             "tau31 is required when imposing an x-face surface layer.");
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
             Tau_lev[TauType::tau21] != nullptr,
             "tau21 is required when imposing an x-face surface layer.");
     } else if (dir == 1) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
             Tau_lev[TauType::tau32] != nullptr,
             "tau32 is required when imposing a y-face surface layer.");
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
             Tau_lev[TauType::tau12] != nullptr,
             "tau12 is required when imposing a y-face surface layer.");
     }
  
     const int klo = sm_index;
     const auto& dxInv = m_geom[lev].InvCellSizeArray();
     const Box domain = m_geom[lev].Domain();
     const int xlo_node = domain.smallEnd(0);
     const int xhi_node = domain.bigEnd(0) + 1;
     const int ylo_node = domain.smallEnd(1);
     const int yhi_node = domain.bigEnd(1) + 1;
     const int zlo_node = domain.smallEnd(2);
     const int zhi_node = domain.bigEnd(2) + 1;
     const bool xlo_surface = m_surface_layer_faces[Orientation::xlo()];
     const bool xhi_surface = m_surface_layer_faces[Orientation::xhi()];
     const bool ylo_surface = m_surface_layer_faces[Orientation::ylo()];
     const bool yhi_surface = m_surface_layer_faces[Orientation::yhi()];
     const bool zlo_surface = m_surface_layer_faces[Orientation::zlo()];
     const bool zhi_surface = m_surface_layer_faces[Orientation::zhi()];
  
     for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
     {
         // Skip boxes that do not own the selected face before acquiring any
         // arrays from collapsed or staggered layouts.  Interior boxes may
         // share the same collapsed coordinates but do not contain surface
         // data for this boundary.
         Box bx = mfi.tilebox();
         const Box valid_bx = mfi.validbox();
         if (m_face.isLow()) {
             if (valid_bx.smallEnd(dir) != klo ||
                 bx.smallEnd(dir) != klo) {
                 continue;
             }
             bx.setBig(dir, bx.smallEnd(dir));
         } else {
             if (valid_bx.bigEnd(dir) != klo ||
                 bx.bigEnd(dir) != klo) {
                 continue;
             }
             bx.setSmall(dir, bx.bigEnd(dir));
         }
  
         // Get field arrays
         const auto cons_arr  = mfs[Vars::cons]->array(mfi);
         const auto velx_arr  = mfs[Vars::xvel]->array(mfi);
         const auto vely_arr  = mfs[Vars::yvel]->array(mfi);
         const auto velz_arr  = mfs[Vars::zvel]->array(mfi);
  
         // Output stresses:
         //            T     Q      U     V
         // X-faces: hfx1   qfx1   t21   t31
         // Y-faces: hfx2   qfx2   t12   t32
         // Z-faces: hfx3   qfx3   t13   t23
  
         // Output stresses nodal locations:
         //            T     Q
         // X-faces: hfx1   qfx1(1,0,0)   t21(V)(1,1,0)   t31(W)(1,0,1)
         // Y-faces: hfx2   qfx2(0,1,0)   t12(U)(1,1,0)   t32(W)(0,1,1)
         // Z-faces: hfx3   qfx3(0,0,1)   t13(U)(1,0,1)   t23(V)(0,1,1)
  
         // Diffusive stress vars
         auto t13_arr =  Tau_lev[TauType::tau13]->array(mfi);
         auto t31_arr = Tau_lev[TauType::tau31]
                      ? Tau_lev[TauType::tau31]->array(mfi) : Array4<Real>{};
  
         auto t23_arr =  Tau_lev[TauType::tau23]->array(mfi);
         auto t32_arr = Tau_lev[TauType::tau32]
                      ? Tau_lev[TauType::tau32]->array(mfi) : Array4<Real>{};
  
  
         auto hfx3_arr = zheat_flux->array(mfi);
         auto qfx3_arr = (zqv_flux)  ? zqv_flux->array(mfi)   : Array4<Real>{};
  
         auto olen_arr   = olen[lev]->array(mfi);
  
         // Rotated stress vars
         auto t11_arr = (m_rotate) ? Tau_lev[TauType::tau11]->array(mfi) : Array4<Real>{};
         auto t22_arr = (m_rotate) ? Tau_lev[TauType::tau22]->array(mfi) : Array4<Real>{};
         auto t33_arr = (m_rotate) ? Tau_lev[TauType::tau33]->array(mfi) : Array4<Real>{};
         auto t12_arr = Tau_lev[TauType::tau12]
                      ? Tau_lev[TauType::tau12]->array(mfi) : Array4<Real>{};
         auto t21_arr = Tau_lev[TauType::tau21]
                      ? Tau_lev[TauType::tau21]->array(mfi) : Array4<Real>{};
  
         auto hfx1_arr = (m_rotate || dir == 0) ? xheat_flux->array(mfi) : Array4<Real>{};
         auto hfx2_arr = (m_rotate || dir == 1) ? yheat_flux->array(mfi) : Array4<Real>{};
         auto qfx1_arr = (xqv_flux && (m_rotate || dir == 0)) ? xqv_flux->array(mfi) : Array4<Real>{};
         auto qfx2_arr = (yqv_flux && (m_rotate || dir == 1)) ? yqv_flux->array(mfi) : Array4<Real>{};
  
         // Terrain
         const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};
  
         // Get average arrays
         const auto *const u_mean      = m_ma.get_average(lev,0);
         const auto *const v_mean      = m_ma.get_average(lev,1);
         const auto *const w_mean      = m_ma.get_average(lev,2);
         const auto *const t_mean      = m_ma.get_average(lev,3);
         const auto *const q_mean      = m_ma.get_average(lev,4);
         const auto *const u_mag_mean  = m_ma.get_average(lev,6);
         const auto *const uw_mag_mean = m_ma.get_average(lev,7);
         const auto *const vw_mag_mean = m_ma.get_average(lev,8);
         const auto *const zref_ptr   = m_ma.get_zref(lev);
  
         const auto um_arr  = u_mean->array(mfi);
         const auto vm_arr  = v_mean->array(mfi);
         const auto wm_arr  = w_mean->array(mfi);
         const auto tm_arr  = t_mean->array(mfi);
         const auto qm_arr  = q_mean->array(mfi);
         const auto umm_arr = u_mag_mean->array(mfi);
         const auto vwmm_arr = (dir == 0) ? vw_mag_mean->array(mfi) : Array4<Real>{};
         const auto uwmm_arr = (dir == 1) ? uw_mag_mean->array(mfi) : Array4<Real>{};
  
         // umm depending on face direction (YZ, XZ, XY)
         const auto dir_umm_arr = ((dir == 0) ? vwmm_arr : ((dir == 1) ? uwmm_arr : umm_arr));
  
         const auto zref_arr = zref_ptr->array(mfi);
         const auto z0_arr   = z_0[lev].array(mfi);
  
         // Get derived arrays
         const auto u_star_arr = u_star[lev]->array(mfi);
         const auto t_star_arr = t_star[lev]->array(mfi);
         const auto q_star_arr = q_star[lev]->array(mfi);
         const auto t_surf_arr = t_surf[lev]->array(mfi);
         const auto q_surf_arr = q_surf[lev]->array(mfi);
         auto surface_source_arr = surface_diagnostic_source[lev]->array(mfi);
  
         // Get LSM fluxes
         auto lmask_arr      = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                       Array4<int> {};
         auto lsm_t_flux_arr  = Array4<const Real> {};
         auto soil_t_flux_arr = Array4<Real> {};
         auto lsm_q_flux_arr = Array4<const Real> {};
         auto lsm_tau13_arr  = Array4<const Real> {};
         auto lsm_tau23_arr  = Array4<const Real> {};
         // LSM tau fields are cell-centered kinematic stresses [m2 s-2].
         // Tau_lev tau13/tau23 are face-centered conservative stresses [N m-2].
         for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
             if (toLower(m_lsm_flux_name[n]) == "t_flux")      { lsm_t_flux_arr  = m_lsm_flux_lev[lev][n]->const_array(mfi); }
             if (toLower(m_lsm_flux_name[n]) == "soil_t_flux") { soil_t_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
             if (toLower(m_lsm_flux_name[n]) == "q_flux")      { lsm_q_flux_arr  = m_lsm_flux_lev[lev][n]->const_array(mfi); }
             if (toLower(m_lsm_flux_name[n]) == "tau13")       { lsm_tau13_arr   = m_lsm_flux_lev[lev][n]->const_array(mfi); }
             if (toLower(m_lsm_flux_name[n]) == "tau23")       { lsm_tau23_arr   = m_lsm_flux_lev[lev][n]->const_array(mfi); }
         }
  
         // Get provider or weight-averaged LSM+Urban fluxes through SurfaceModel.
         const SurfaceFluxView surf_flux = (use_surface_model && surf_model_fluxes)
             ? m_surf_model->get_surface_flux_view(lev, mfi) : SurfaceFluxView{};
         const auto surf_tflux_arr = surf_flux.t_flux;
         const auto surf_qflux_arr = surf_flux.q_flux;
         const auto surf_uflux_arr = surf_flux.tau13;
         const auto surf_vflux_arr = surf_flux.tau23;
         const bool use_surface_model_fluxes = use_surface_model && surf_model_fluxes;
         const auto t_flux_arr = use_surface_model_fluxes ? surf_tflux_arr : lsm_t_flux_arr;
         const auto q_flux_arr = use_surface_model_fluxes ? surf_qflux_arr : lsm_q_flux_arr;
         const auto surface_tau13_arr = use_surface_model_fluxes ? surf_uflux_arr : lsm_tau13_arr;
         const auto surface_tau23_arr = use_surface_model_fluxes ? surf_vflux_arr : lsm_tau23_arr;
         const bool has_lsm_t_flux = static_cast<bool>(t_flux_arr);
         const bool is_custom = (flux_type == FluxCalcType::CUSTOM);
         const bool is_rico   = (flux_type == FluxCalcType::RICO);
  
  
         // Rho*Theta flux
         //============================================================================
         const bool is_low_face = m_face.isLow();
  
         amrex::Array4<amrex::Real> surface_flux_arr;
         if (mf_surface_fluxes != nullptr) {
             surface_flux_arr = mf_surface_fluxes->array(mfi);
         }
  
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             // Valid theta flux from LSM and over land. The LSM writes the
             // lsm_undefined sentinel for cells it did not process (sea-ice /
             // open water); fall back to MOST there instead of applying garbage.
             Real Tflux;
             int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
             const bool lsm_flux_is_valid = (t_flux_arr) ? (t_flux_arr(i,j,0) < lsm_undefined) :
                                                               false;
             const bool has_land_and_flux = (is_land == 1 && lsm_flux_is_valid);
             if (t_flux_arr && has_land_and_flux) {
                 // LSM flux MultiFabs store kinematic fluxes for MOST parameter
                 // updates. The applied hfx array stores the conservative RHS flux.
                 Tflux = cons_arr(i,j,k,Rho_comp) * t_flux_arr(i,j,0);
                 if(surface_flux_arr) {
                     surface_flux_arr(i,j,0,0) = Tflux;
                 }
             } else if (is_land == 2) { // no temperature flux within buildings
                 Tflux = zero;
                 if(surface_flux_arr) {
                     surface_flux_arr(i,j,0,0) = Tflux;
                 }
             } else {
                 Tflux = flux_comp.compute_t_flux(i, j, k, dir,
                                                  cons_arr, velx_arr, vely_arr, velz_arr,
                                                  dir_umm_arr, tm_arr, u_star_arr,
                                                  t_star_arr, t_surf_arr);
                 // NOTE: do NOT write the MOST-fallback flux back into lsm_t_flux_arr.
                 // Doing so flips a sentinel (water/unprocessed) cell to "valid LSM"
                 // on the next step, so a MOST-derived value is re-read as an LSM flux
                 // Only Noah-MP should populate the LSM cache.
                 if (surface_flux_arr) {
                     surface_flux_arr(i,j,0,0) = Tflux;
                 }
             }
  
             if (soil_t_flux_arr && is_land == 1) {
                 soil_t_flux_arr(i,j,k) = Tflux / cons_arr(i,j,k,Rho_comp);
             }
  
             surface_source_arr(i,j,k) = surface_diagnostics::to_plot_value(
                 surface_diagnostics::classify_scalar_source(
                     is_custom, is_rico, is_land, has_lsm_t_flux, lsm_flux_is_valid));
  
             // Do scalar flux rotations?
             if (rotate) {
                 rotate_scalar_flux(i, j, k, dir, Tflux, dxInv, zphys_arr,
                                    hfx1_arr, hfx2_arr, hfx3_arr);
             } else {
                 // swap sign for upward faces
                 if (!is_low_face && dir != 2) {
                     Tflux = -Tflux;
                 }
  
                 // write out to corresponding face
                 if (dir == 0) {
                     if (!is_low_face) {
                         hfx1_arr(i+1,j,k) = Tflux;
                     } else {
                         hfx1_arr(i,j,k) = Tflux;
                     }
                 } else if (dir == 1) {
                     if (!is_low_face) {
                         hfx2_arr(i,j+1,k) = Tflux;
                     } else {
                         hfx2_arr(i,j,k) = Tflux;
                     }
                 } else {
                     if (!is_low_face) {
                         hfx3_arr(i,j,k+1) = Tflux;
                     } else {
                         hfx3_arr(i,j,k) = Tflux;
                     }
                 }
             }
         });
  
         // Rho*Qv flux
         //============================================================================
         if (use_moisture) {
             ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 // Valid qv flux from LSM and over land (sentinel -> fall back to MOST)
                 Real Qflux;
                 int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
                 const bool lsm_flux_is_valid = (q_flux_arr) ? (q_flux_arr(i,j,0) < lsm_undefined) :
                                                                   false;
                 const bool has_land_and_flux = (is_land == 1 && lsm_flux_is_valid);
                 if (q_flux_arr && has_land_and_flux) {
                     // LSM flux MultiFabs store kinematic fluxes for MOST parameter
                     // updates. The applied qfx array stores the conservative RHS flux.
                     Qflux = cons_arr(i,j,k,Rho_comp) * q_flux_arr(i,j,0);
                     if(surface_flux_arr) {
                         surface_flux_arr(i,j,0,1) = Qflux;
                     }
                 } else if (is_land == 2) { // no moisture flux within buildings
                     Qflux = zero;
                     if(surface_flux_arr) {
                         surface_flux_arr(i,j,0,1) = Qflux;
                     }
                 } else {
                     Qflux = flux_comp.compute_q_flux(i, j, k, dir,
                                                      cons_arr, velx_arr, vely_arr, velz_arr,
                                                      dir_umm_arr, qm_arr, u_star_arr,
                                                      q_star_arr, q_surf_arr);
                     // NOTE: no writeback into lsm_q_flux_arr -- see the matching
                     // t_flux note above.
                     if (surface_flux_arr) {
                         surface_flux_arr(i,j,0,1) = Qflux;
                     }
                 }
  
                 // Do scalar flux rotations?
                 if (rotate) {
                     rotate_scalar_flux(i, j, k, dir, Qflux, dxInv, zphys_arr,
                                        qfx1_arr, qfx2_arr, qfx3_arr);
                 } else {
                     // swap sign for upward faces
                     if (!is_low_face && dir != 2) {
                         Qflux = -Qflux;
                     }
  
                     // write out to corresponding face
                     if (dir == 0) {
                         if (!is_low_face) {
                             qfx1_arr(i+1,j,k) = Qflux;
                         } else {
                             qfx1_arr(i,j,k) = Qflux;
                         }
                     } else if (dir == 1) {
                         if (!is_low_face) {
                             qfx2_arr(i,j+1,k) = Qflux;
                         } else {
                             qfx2_arr(i,j,k) = Qflux;
                         }
                     } else {
                         if (!is_low_face) {
                             qfx3_arr(i,j,k+1) = Qflux;
                         } else {
                             qfx3_arr(i,j,k) = Qflux;
                         }
                     }
                 }
             });
         } // custom
  
         if (!rotate) {
             // Rho*u flux
             //============================================================================
             const IntVect stressx_nodal = (dir == 2) ? IntVect(1,0,1) : IntVect(1,1,0);
             Box bxx = convert(bx, stressx_nodal);
             const int stressx_face_index = is_low_face
                                          ? m_geom[lev].Domain().smallEnd(dir)
                                          : m_geom[lev].Domain().bigEnd(dir) + 1;
             bxx.setRange(dir, stressx_face_index);
             ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 // Valid tau13 from LSM and over land. A side that is land but
                 // whose LSM flux is the sentinel (sea-ice / open water) is treated
                 // as non-LSM so that side uses the MOST stress instead.
                 Real stressx;
                 int is_land_hi = (lmask_arr) ? lmask_arr(i  ,j,0) : 1;
                 int is_land_lo = (lmask_arr) ? lmask_arr(i-1,j,0) : 1;
                 const bool lsm_hi_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
                     static_cast<bool>(surface_tau13_arr), is_land_hi == 1,
                     surface_tau13_arr ? surface_tau13_arr(i  ,j,0) : zero, lsm_undefined);
                 const bool lsm_lo_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
                     static_cast<bool>(surface_tau13_arr), is_land_lo == 1,
                     surface_tau13_arr ? surface_tau13_arr(i-1,j,0) : zero, lsm_undefined);
                 const bool has_land_and_flux_hi = (is_land_hi == 1 && lsm_hi_flux_is_valid);
                 const bool has_land_and_flux_lo = (is_land_lo == 1 && lsm_lo_flux_is_valid);
                 if (surface_tau13_arr && (has_land_and_flux_hi || has_land_and_flux_lo)) {
                     const Real rho_hi = cons_arr(i  ,j,k,Rho_comp);
                     const Real rho_lo = cons_arr(i-1,j,k,Rho_comp);
                     const Real most_stress = (!has_land_and_flux_hi || !has_land_and_flux_lo) ?
                         flux_comp.compute_u_flux(i, j, k, dir,
                                                  cons_arr, velx_arr, vely_arr, velz_arr,
                                                  dir_umm_arr, um_arr, vm_arr, wm_arr, u_star_arr) : zero;
                     const auto result = surface_layer_stress::combine_lsm_and_most_stress(
                         rho_lo, rho_hi, surface_tau13_arr(i-1,j,0), surface_tau13_arr(i,j,0),
                         has_land_and_flux_lo, has_land_and_flux_hi, most_stress);
                     stressx = result.face_stress;
                     // NOTE: do NOT write the MOST-fallback stress back into the
                     // cell-centered lsm_tau13_arr. This face-indexed ParallelFor
                     // touches cells (i) and (i-1), so each cell is written by two
                     // adjacent face threads in the same launch -> nondeterministic
                     // write-write race on GPU (ERF #3446). It also spuriously flips
                     // a sentinel (water/unprocessed) cell to "valid LSM" for the
                     // next step. The face stress is fully determined here; the LSM
                     // cache is (re)filled only by Noah-MP. Matches baseline 3ab899d3.
                 } else if (is_land_hi == 2 || is_land_lo == 2) { // no stress within buildings
                     stressx = zero;
                 } else {
                     stressx = flux_comp.compute_u_flux(i, j, k, dir,
                                                        cons_arr, velx_arr, vely_arr, velz_arr,
                                                        dir_umm_arr, um_arr, vm_arr, wm_arr, u_star_arr);
                 }
  
                 // write out to corresponding face
                 if (dir == 0) {
                     t21_arr(i,j,k) = stressx;
                     if (t12_arr &&
                         (!ylo_surface || j != ylo_node) &&
                         (!yhi_surface || j != yhi_node)) {
                         t12_arr(i,j,k) = stressx;
                     }
                 } else if (dir == 1) {
                     t12_arr(i,j,k) = stressx;
                     if (t21_arr &&
                         (!xlo_surface || i != xlo_node) &&
                         (!xhi_surface || i != xhi_node)) {
                         t21_arr(i,j,k) = stressx;
                     }
                 } else {
                     t13_arr(i,j,k) = stressx;
                     if (t31_arr &&
                         (!xlo_surface || i != xlo_node) &&
                         (!xhi_surface || i != xhi_node)) {
                         t31_arr(i,j,k) = stressx;
                     }
                 }
             });
  
             // Rho*v flux
             //============================================================================
             const IntVect stressy_nodal = (dir == 0) ? IntVect(1,0,1) : IntVect(0,1,1);
             Box bxy = convert(bx, stressy_nodal);
             const int stressy_face_index = is_low_face
                                          ? m_geom[lev].Domain().smallEnd(dir)
                                          : m_geom[lev].Domain().bigEnd(dir) + 1;
             bxy.setRange(dir, stressy_face_index);
             ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 // Valid tau23 from LSM and over land (sentinel side -> MOST stress)
                 Real stressy;
                 int is_land_hi = (lmask_arr) ? lmask_arr(i,j  ,0) : 1;
                 int is_land_lo = (lmask_arr) ? lmask_arr(i,j-1,0) : 1;
                 const bool lsm_hi_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
                     static_cast<bool>(surface_tau23_arr), is_land_hi == 1,
                     surface_tau23_arr ? surface_tau23_arr(i,j  ,0) : zero, lsm_undefined);
                 const bool lsm_lo_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
                     static_cast<bool>(surface_tau23_arr), is_land_lo == 1,
                     surface_tau23_arr ? surface_tau23_arr(i,j-1,0) : zero, lsm_undefined);
                 const bool has_land_and_flux_hi = (is_land_hi == 1 && lsm_hi_flux_is_valid);
                 const bool has_land_and_flux_lo = (is_land_lo == 1 && lsm_lo_flux_is_valid);
                 if (surface_tau23_arr && (has_land_and_flux_hi || has_land_and_flux_lo)) {
                     const Real rho_hi = cons_arr(i,j  ,k,Rho_comp);
                     const Real rho_lo = cons_arr(i,j-1,k,Rho_comp);
                     const Real most_stress = (!has_land_and_flux_hi || !has_land_and_flux_lo) ?
                         flux_comp.compute_v_flux(i, j, k, dir,
                                                  cons_arr, velx_arr, vely_arr, velz_arr,
                                                  dir_umm_arr, um_arr, vm_arr, wm_arr, u_star_arr) : zero;
                     const auto result = surface_layer_stress::combine_lsm_and_most_stress(
                         rho_lo, rho_hi, surface_tau23_arr(i,j-1,0), surface_tau23_arr(i,j,0),
                         has_land_and_flux_lo, has_land_and_flux_hi, most_stress);
                     stressy = result.face_stress;
                     // NOTE: no writeback into cell-centered lsm_tau23_arr -- see the
                     // matching tau13 note above (ERF #3446 write-write race + stale
                     // sentinel-becomes-valid). Face stress is complete here.
                 } else if (is_land_hi == 2 || is_land_lo == 2) { // no stress within buildings
                     stressy = zero;
                 } else {
                     stressy = flux_comp.compute_v_flux(i, j, k, dir,
                                                        cons_arr, velx_arr, vely_arr, velz_arr,
                                                        dir_umm_arr, um_arr, vm_arr, wm_arr, u_star_arr);
                 }
  
                 // write out to corresponding face
                 if (dir == 0) {
                     t31_arr(i,j,k) = stressy;
                     if (t13_arr &&
                         (!zlo_surface || k != zlo_node) &&
                         (!zhi_surface || k != zhi_node)) {
                         t13_arr(i,j,k) = stressy;
                     }
                 } else if (dir == 1) {
                     t32_arr(i,j,k) = stressy;
                     if (t23_arr &&
                         (!zlo_surface || k != zlo_node) &&
                         (!zhi_surface || k != zhi_node)) {
                         t23_arr(i,j,k) = stressy;
                     }
                 } else {
                     t23_arr(i,j,k) = stressy;
                     if (t32_arr &&
                         (!ylo_surface || j != ylo_node) &&
                         (!yhi_surface || j != yhi_node)) {
                         t32_arr(i,j,k) = stressy;
                     }
                 }
             });
         } else {
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dir == 2 && m_face.isLow(), "Stress rotation only supported for zlo face");
             // All fluxes with rotation
             //============================================================================
             Box bxxy = convert(bx, IntVect(1,1,0));
             ParallelFor(bxxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 Real stresst = flux_comp.compute_u_flux(i, j, k, dir,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, um_arr, vm_arr, wm_arr, u_star_arr);
                 rotate_stress_tensor(i, j, k, dir, stresst, dxInv, zphys_arr,
                                      velx_arr, vely_arr, velz_arr,
                                      t11_arr, t22_arr, t33_arr,
                                      t12_arr, t21_arr,
                                      t13_arr, t31_arr,
                                      t23_arr, t32_arr);
                 });
         }
  
         // For models that do not do iterations to yield u*/T*/q*,
         // fill these values from the fluxes that were computed.
  
         // NOTE: For LSM, this has been handled in "compute_sfc_params_from_lsm_fluxes"
         // NOTE: Fluxes here are for conserved quantities, we divide by rho
         if (flux_type == FluxCalcType::BULK_COEFF) {
             constexpr Real eps  = std::numeric_limits<Real>::epsilon();
             bool l_use_moisture = use_moisture;
             ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/)
             {
                 Real rho = cons_arr(i,j,klo,Rho_comp);
                 Real Thd = cons_arr(i,j,klo,RhoTheta_comp) / rho;
                 Real qv  = (l_use_moisture) ? cons_arr(i,j,klo,RhoQ1_comp) / rho : zero;
                 Real Thv = Thd * (one + epsv*qv);
  
                 Real tau = std::sqrt( t13_arr(i,j,klo)/rho * t13_arr(i,j,klo)/rho
                                     + t23_arr(i,j,klo)/rho * t23_arr(i,j,klo)/rho );
                 u_star_arr(i,j,0) = amrex::max(std::sqrt(tau),eps);
  
                 if (hfx3_arr(i,j,klo)>=zero) {
                     t_star_arr(i,j,0) = amrex::min(-hfx3_arr(i,j,klo) / (rho * u_star_arr(i,j,0)),-eps);
                 } else {
                     t_star_arr(i,j,0) = amrex::max(-hfx3_arr(i,j,klo) / (rho * u_star_arr(i,j,0)),eps);
                 }
                 if (!l_use_moisture) {
                     q_star_arr(i,j,0) = zero;
                 } else if (qfx3_arr(i,j,klo)>=zero) {
                     q_star_arr(i,j,0) = amrex::min(-qfx3_arr(i,j,klo) / (rho * u_star_arr(i,j,0)),-eps);
                 } else {
                     q_star_arr(i,j,0) = amrex::max(-qfx3_arr(i,j,klo) / ( rho * u_star_arr(i,j,0)),eps);
                 }
                 Real tstv = t_star_arr(i,j,0)*(one + epsv*qv) + epsv*Thd*q_star_arr(i,j,0);
                 tstv = (tstv >= zero) ? amrex::max(tstv, eps) : amrex::min(tstv, -eps);
                 olen_arr(i,j,0) = ( u_star_arr(i,j,0)  * u_star_arr(i,j,0) * Thv ) /
                                   ( KAPPA * CONST_GRAV * tstv );
                 z0_arr(i,j,0)  = Compute_roughness(zref_arr(i,j,0),   olen_arr(i,j,0),
                                                     umm_arr(i,j,0), u_star_arr(i,j,0));
             });
         }
  
     } // mfiter
  
     fill_planar_boundary(lev, *surface_diagnostic_source[lev]);
 }

Here is the call graph for this function:

◆ compute_SurfaceLayer_bcs_EB() [1/2]

template<typename FluxCalc >

void SurfaceLayer::compute_SurfaceLayer_bcs_EB	(	const int &	lev,
		amrex::Vector< const amrex::MultiFab * >	mfs,
		amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab >>> &	Tau_lev,
		amrex::MultiFab *	xheat_flux,
		amrex::MultiFab *	yheat_flux,
		amrex::MultiFab *	zheat_flux,
		amrex::MultiFab *	xqv_flux,
		amrex::MultiFab *	yqv_flux,
		amrex::MultiFab *	zqv_flux,
		const FluxCalc &	flux_comp
	)

Compute embedded-boundary surface-layer flux boundary conditions.

Parameters

[in]	lev	level index
[in]	mfs	state and velocity fields used by the BC computation
[in,out]	Tau_lev	EB stress fields to fill
[in,out]	xheat_flux	x-face heat flux field
[in,out]	yheat_flux	y-face heat flux field
[in,out]	zheat_flux	z-face heat flux field
[in,out]	xqv_flux	x-face moisture flux field
[in,out]	yqv_flux	y-face moisture flux field
[in,out]	zqv_flux	z-face moisture flux field
[in]	flux_comp	flux-computation functor

◆ compute_SurfaceLayer_bcs_EB() [2/2]

template<typename FluxCalc >

void SurfaceLayer::compute_SurfaceLayer_bcs_EB	(	const int &	lev,
		Vector< const MultiFab * >	mfs,
		Vector< Vector< std::unique_ptr< MultiFab >>> &	Tau_EB,
		[[maybe_unused] ] MultiFab *	xheat_flux,
		[[maybe_unused] ] MultiFab *	yheat_flux,
		MultiFab *	Hfx3_EB,
		[[maybe_unused] ] MultiFab *	xqv_flux,
		[[maybe_unused] ] MultiFab *	yqv_flux,
		[[maybe_unused] ] MultiFab *	zqv_flux,
		const FluxCalc &	flux_comp
	)

Function to calculate MOST fluxes for EB.

Parameters

[in]	lev	Current level
[in]	mfs	State MultiFabs used to compute the EB boundary fluxes
[in,out]	Tau_EB	EB diffusive stress MultiFabs populated with surface stresses
[in,out]	xheat_flux	x-face EB heat-flux MultiFab, currently unused
[in,out]	yheat_flux	y-face EB heat-flux MultiFab, currently unused
[in,out]	Hfx3_EB	EB heat-flux MultiFab populated with scalar surface flux
[in,out]	xqv_flux	x-face EB moisture-flux MultiFab, currently unused
[in,out]	yqv_flux	y-face EB moisture-flux MultiFab, currently unused
[in,out]	zqv_flux	z-face EB moisture-flux MultiFab, currently unused
[in]	flux_comp	EB flux-calculation functor used to compute scalar and momentum fluxes

 {
     const int dir = m_face.coordDir();
     // Get EB flags for all centerings
     const auto& cc_factory = m_eb_vec[lev]->get_const_factory();
     const auto& cc_flags = cc_factory->getMultiEBCellFlagFab();
     const auto& cc_vfrac = cc_factory->getVolFrac();
  
     const auto& u_factory = m_eb_vec[lev]->get_u_const_factory();
     const auto& u_flags = u_factory->getMultiEBCellFlagFab();
     const auto& u_vfrac = u_factory->getVolFrac();
  
     const auto& v_factory = m_eb_vec[lev]->get_v_const_factory();
     const auto& v_flags = v_factory->getMultiEBCellFlagFab();
     const auto& v_vfrac = v_factory->getVolFrac();
  
     const auto& w_factory = m_eb_vec[lev]->get_w_const_factory();
     const auto& w_flags = w_factory->getMultiEBCellFlagFab();
     const auto& w_vfrac = w_factory->getVolFrac();
  
     // EB does not currently have a cell-centered scalar-source classification.
     // Keep the provenance mask missing rather than inventing face-aware
     // semantics for the staggered stress path.
     surface_diagnostic_source[lev]->setVal(
         surface_diagnostics::to_plot_value(surface_diagnostics::SurfaceDiagnosticSource::Missing));
  
     for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
     {
         // Get flags for this box (all centerings)
         const auto& cc_flag = cc_flags[mfi];
         const auto& u_flag = u_flags[mfi];
         const auto& v_flag = v_flags[mfi];
         const auto& w_flag = w_flags[mfi];
  
         // Skip boxes that have no cut cells at any centering
         if (cc_flag.getType() != FabType::singlevalued &&
             u_flag.getType() != FabType::singlevalued &&
             v_flag.getType() != FabType::singlevalued &&
             w_flag.getType() != FabType::singlevalued
         ) continue;
  
         // Get EB flag and volfrac arrays
         auto const cc_flag_arr = cc_flag.const_array();
         auto const u_flag_arr = u_flag.const_array();
         auto const v_flag_arr = v_flag.const_array();
         auto const w_flag_arr = w_flag.const_array();
  
         auto const cc_vfrac_arr = cc_vfrac.const_array(mfi);
         auto const u_vfrac_arr = u_vfrac.const_array(mfi);
         auto const v_vfrac_arr = v_vfrac.const_array(mfi);
         auto const w_vfrac_arr = w_vfrac.const_array(mfi);
  
         // Get boundary normals only if cut cells exist
         auto const bnorm_arr = (cc_flag.getType() == FabType::singlevalued) ?
             cc_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
         auto const u_bnorm_arr = (u_flag.getType() == FabType::singlevalued) ?
             u_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
         auto const v_bnorm_arr = (v_flag.getType() == FabType::singlevalued) ?
             v_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
         auto const w_bnorm_arr = (w_flag.getType() == FabType::singlevalued) ?
             w_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
  
         // Get field arrays
         const auto cons_arr  = mfs[Vars::cons]->array(mfi);
         const auto velx_arr  = mfs[Vars::xvel]->array(mfi);
         const auto vely_arr  = mfs[Vars::yvel]->array(mfi);
         const auto velz_arr  = mfs[Vars::zvel]->array(mfi);
  
         // Diffusive stress vars - t13 and t23 components for all grid types
         auto u_t13_arr = Tau_EB[EBTauType::tau_eb13][EBGridType::xface]->array(mfi);
         auto v_t13_arr = Tau_EB[EBTauType::tau_eb13][EBGridType::yface]->array(mfi);
         auto w_t13_arr = Tau_EB[EBTauType::tau_eb13][EBGridType::zface]->array(mfi);
  
         auto u_t23_arr = Tau_EB[EBTauType::tau_eb23][EBGridType::xface]->array(mfi);
         auto v_t23_arr = Tau_EB[EBTauType::tau_eb23][EBGridType::yface]->array(mfi);
         auto w_t23_arr = Tau_EB[EBTauType::tau_eb23][EBGridType::zface]->array(mfi);
  
         auto hfx3_arr = Hfx3_EB->array(mfi);
  
         // Get average arrays
         const auto *const u_mean     = m_ma.get_average(lev,0);
         const auto *const v_mean     = m_ma.get_average(lev,1);
         // const auto *const w_mean     = m_ma.get_average(lev,2);
  
         const auto *const t_mean     = m_ma.get_average(lev,3);
         // const auto *const q_mean     = m_ma.get_average(lev,4);
         const auto *const u_mag_mean = m_ma.get_average(lev,6);
         const auto *const uw_mag_mean = m_ma.get_average(lev,7);
         const auto *const vw_mag_mean = m_ma.get_average(lev,8);
  
         const auto um_arr  = u_mean->array(mfi);
         const auto vm_arr  = v_mean->array(mfi);
         // const auto wm_arr  = w_mean->array(mfi);
         const auto tm_arr  = t_mean->array(mfi);
         // const auto qm_arr  = q_mean->array(mfi);
         const auto umm_arr = u_mag_mean->array(mfi);
         const auto vwmm_arr = (dir == 0) ? vw_mag_mean->array(mfi) : Array4<Real>{};
         const auto uwmm_arr = (dir == 1) ? uw_mag_mean->array(mfi) : Array4<Real>{};
  
         // umm depending on face direction (YZ, XZ, XY)
         const auto dir_umm_arr = ((dir == 0) ? vwmm_arr : ((dir == 1) ? uwmm_arr : umm_arr));
  
         // Get derived arrays
         const auto u_star_arr = u_star[lev]->array(mfi);
         const auto t_star_arr = t_star[lev]->array(mfi);
         const auto t_surf_arr = t_surf[lev]->array(mfi);
  
         // Rho*Theta flux
         //============================================================================
         Box bx = mfi.tilebox();
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (cc_flag_arr(i,j,k).isSingleValued()) {
                 Real Tflux = flux_comp.compute_t_flux(i, j, k,
                                                     cons_arr, velx_arr, vely_arr, velz_arr,
                                                     dir_umm_arr, tm_arr, u_star_arr,
                                                     t_star_arr, t_surf_arr,
                                                     u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                     bnorm_arr);
                 hfx3_arr(i,j,k) = Tflux;
             }
         });
  
         // Rho*u flux
         //============================================================================
         Box bxx = surroundingNodes(bx,0);
         Box bxy = surroundingNodes(bx,1);
         Box bxz = surroundingNodes(bx,2);
         ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (u_flag_arr(i,j,k).isSingleValued()) {
                 Real stressx = flux_comp.compute_u_flux(i, j, k,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, um_arr, u_star_arr,
                                                         u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                         cc_vfrac_arr, cc_flag_arr,
                                                         u_bnorm_arr, 0);
                 u_t13_arr(i,j,k) = stressx;
             }
         });
         ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (v_flag_arr(i,j,k).isSingleValued()) {
                 Real stressx = flux_comp.compute_u_flux(i, j, k,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, um_arr, u_star_arr,
                                                         u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                         cc_vfrac_arr, cc_flag_arr,
                                                         v_bnorm_arr, 1);
                 v_t13_arr(i,j,k) = stressx;
             }
         });
         ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (w_flag_arr(i,j,k).isSingleValued()) {
                 Real stressx = flux_comp.compute_u_flux(i, j, k,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, um_arr, u_star_arr,
                                                         u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                         cc_vfrac_arr, cc_flag_arr,
                                                         w_bnorm_arr, 2);
                 w_t13_arr(i,j,k) = stressx;
             }
         });
  
         // Rho*v flux
         //============================================================================
         ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (u_flag_arr(i,j,k).isSingleValued()) {
                 Real stressy = flux_comp.compute_v_flux(i, j, k,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, vm_arr, u_star_arr,
                                                         u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                         cc_vfrac_arr, cc_flag_arr, u_bnorm_arr, 0);
                 u_t23_arr(i,j,k) = stressy;
             }
         });
         ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (v_flag_arr(i,j,k).isSingleValued()) {
                 Real stressy = flux_comp.compute_v_flux(i, j, k,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, vm_arr, u_star_arr,
                                                         u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                         cc_vfrac_arr, cc_flag_arr, v_bnorm_arr, 1);
                 v_t23_arr(i,j,k) = stressy;
             }
         });
         ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             if (w_flag_arr(i,j,k).isSingleValued()) {
                 Real stressy = flux_comp.compute_v_flux(i, j, k,
                                                         cons_arr, velx_arr, vely_arr, velz_arr,
                                                         dir_umm_arr, vm_arr, u_star_arr,
                                                         u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
                                                         cc_vfrac_arr, cc_flag_arr, w_bnorm_arr, 2);
                 w_t23_arr(i,j,k) = stressy;
             }
         });
     } // mfiter
  
 }

Here is the call graph for this function:

◆ computes_pblh()

bool SurfaceLayer::computes_pblh ( ) const

inline

Do we actually compute the PBL height? If not then the field returned by get_pblh holds only the value it was initialized with, and must not be reported as a diagnostic.

1349 { return (pblh_type != PBLHeightCalcType::None); }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ computes_w_star()

bool SurfaceLayer::computes_w_star ( ) const

inline

Do we actually compute w*? If not then the field returned by get_w_star holds only the value it was initialized with, and must not be reported as a diagnostic.

1343 { return m_include_wstar; }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ define_pblh_columns()

void SurfaceLayer::define_pblh_columns	(	const int &	lev,
		const amrex::BoxArray &	ba,
		const amrex::DistributionMapping &	dm
	)

private

Build the columns on which compute_pblh runs the PBL-height estimator when the grids of a level do not all start at the ground (see PBLHColumns). Called when the grids change.

Parameters

[in]	lev	Current level
[in]	ba	BoxArray of the state at this level
[in]	dm	DistributionMapping of the state at this level

 {
     PBLHColumns& cols = m_pblh_columns[lev];
     cols = PBLHColumns{};
     cols.ba = ba;
     cols.dm = dm;
  
     const int k_ground = m_geom[lev].Domain().smallEnd(2);
     const int k_top    = m_geom[lev].Domain().bigEnd(2);
     for (int ib = 0; ib < ba.size(); ++ib) {
         if (ba[ib].smallEnd(2) != k_ground) { cols.needed = true; }
     }
  
     // The runs of cells in z, each as one box
     const BoxArray ba_joined = join_boxes_stacked_in_z(ba);
     for (int ib = 0; ib < ba_joined.size(); ++ib) {
         if (ba_joined[ib].smallEnd(2) != k_ground || ba_joined[ib].bigEnd(2) != k_top) { cols.full_height = false; }
     }
     if (!cols.needed) { return; }
  
     // With EB the surface-layer fields live on the 3D grids and the estimator writes their
     // k = 0 plane, which a box that does not start at the ground does not hold
     if (m_terrain_type == TerrainType::EB) {
         Abort("erf.most.pblh_calc = MYNN25 with EB needs every grid at level " + std::to_string(lev) +
               " to start at the bottom of the domain: the PBL height is written into the lowest "
               "plane of each grid.  Choose grids that are not decomposed in z "
               "(amr.max_grid_size_z) and refined regions that reach the ground.");
     }
  
     // Those of the runs of cells in z that start at the ground are the columns.  A column goes
     // to the rank that owns its lowest corner cell, which keeps most of the copies to and from
     // the columns on the rank.
     BoxList bl_col(IndexType::TheCellType());
     BoxList bl_col2d(IndexType::TheCellType());
     Vector<int> pmap;
     for (int ib = 0; ib < ba_joined.size(); ++ib) {
         const Box& b = ba_joined[ib];
         if (b.smallEnd(2) != k_ground) { continue; }
         const auto& owners = ba.intersections(Box(b.smallEnd(), b.smallEnd()));
         AMREX_ALWAYS_ASSERT(!owners.empty());
         Box b2d(b);
         b2d.setRange(2, k_ground);
         bl_col.push_back(b);
         bl_col2d.push_back(b2d);
         pmap.push_back(dm[owners[0].first]);
     }
     if (bl_col.isEmpty()) { return; }
  
     cols.ba_col   = BoxArray(std::move(bl_col));
     cols.ba_col2d = BoxArray(std::move(bl_col2d));
     cols.dm_col   = DistributionMapping(std::move(pmap));
  
     // These hold no data.  AMReX drops the communication metadata of a BoxArray and
     // DistributionMapping pair with the last FabArray built on it, and the MultiFabs that
     // compute_pblh builds on the columns are temporaries.
     cols.hold_col.define(cols.ba_col, cols.dm_col, 1, 0, MFInfo().SetAlloc(false));
     cols.hold_col2d.define(cols.ba_col2d, cols.dm_col, 1, 0, MFInfo().SetAlloc(false));
 }

Here is the call graph for this function:

◆ fill_lateral_surface_parameter_ghosts()

void SurfaceLayer::fill_lateral_surface_parameter_ghosts	(	const int &	lev,
		amrex::MultiFab *	selected_field = `nullptr`
	)

Fills surface layer data for lateral faces. A normal FillBoundary on lateral faces would bring in data from other ranks or interior grids that do not own data on the lateral face.

Parameters

[in]	lev	level index
[in,out]	selected_field	if non-null, fill only this field; otherwise fill all lateral surface fields

 {
     const int dir = m_face.coordDir();
     AMREX_ALWAYS_ASSERT(dir < 2);
  
     const iMultiFab& surface_mask = *m_lmask_lev[lev][0];
     const int face_index = m_face.isLow()
         ? m_geom[lev].Domain().smallEnd(dir)
         : m_geom[lev].Domain().bigEnd(dir);
     const auto& mask_ba = surface_mask.boxArray();
  
     IntVect period = m_geom[lev].periodicity().intVect();
     period[dir] = 0; // The selected wall is not a periodic source plane.
     const Periodicity tangential_periodicity(period);
  
     Vector<MultiFab*> fields;
     if (selected_field) {
         fields.push_back(selected_field);
     } else {
         const Vector<MultiFab*> all_fields{
             u_star[lev].get(), t_star[lev].get(), q_star[lev].get(), olen[lev].get(),
             t_surf[lev].get(), q_surf[lev].get(), pblh[lev].get(),
             surface_diagnostic_source[lev].get()};
         for (MultiFab* field : all_fields) {
             if (field) { fields.push_back(field); }
         }
         if (m_include_wstar && w_star[lev]) { fields.push_back(w_star[lev].get()); }
     }
  
     for (MultiFab* field : fields) {
         // Each parameter field owns its own index type and BoxArray. Build
         // the compact layout from that field rather than reusing u_star's
         // layout for all parameters.
         BoxList face_boxes;
         Vector<int> face_pmap;
         Vector<int> source_indices;
         const auto& source_ba = field->boxArray();
         const auto& source_pmap = field->DistributionMap().ProcessorMap();
         for (int ibox = 0; ibox < mask_ba.size(); ++ibox) {
             const Box& mask_box = mask_ba[ibox];
             const bool owns_face = m_face.isLow()
                 ? mask_box.smallEnd(dir) == face_index
                 : mask_box.bigEnd(dir) == face_index;
             if (!owns_face) { continue; }
  
             face_boxes.push_back(source_ba[ibox]);
             face_pmap.push_back(source_pmap[ibox]);
             source_indices.push_back(ibox);
         }
         if (face_boxes.isEmpty()) { continue; }
  
         BoxArray face_ba(std::move(face_boxes));
         DistributionMapping face_dm(std::move(face_pmap));
         MultiFab face_values(face_ba, face_dm, 1, field->nGrowVect());
  
         // Selectively gather only face-owned FABs. A direct FillBoundary on
         // field would also treat interior-grid FABs as valid sources because
         // all lateral FABs collapse onto the same wall plane.
         for (MFIter fmfi(face_values, false); fmfi.isValid(); ++fmfi) {
             const int compact_index = fmfi.index();
             const int source_index = source_indices[compact_index];
             const auto src = field->const_array(source_index);
             const auto dst = face_values.array(fmfi);
             // Preserve the source FAB's existing ghosts.  A nonperiodic
             // physical ghost has no FillBoundary source, so leaving the
             // temporary FAB uninitialized would replace a valid local value
             // with garbage at mixed-face corners.
             Box source_fab = field->boxArray()[source_index];
             source_fab.grow(field->nGrowVect());
             const Box copy_box = fmfi.fabbox() & source_fab;
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                 !copy_box.isEmpty(),
                 "Temporary surface-layer source copy has no valid overlap.");
             ParallelFor(copy_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 dst(i,j,k) = src(i,j,k);
             });
         }
         Gpu::streamSynchronize();
  
         face_values.FillBoundary(tangential_periodicity);
  
         for (MFIter fmfi(face_values, false); fmfi.isValid(); ++fmfi) {
             const int compact_index = fmfi.index();
             const int source_index = source_indices[compact_index];
             const auto src = face_values.const_array(fmfi);
             const auto dst = field->array(source_index);
             Box source_fab = field->boxArray()[source_index];
             source_fab.grow(field->nGrowVect());
             const Box copy_box = fmfi.fabbox() & source_fab;
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                 !copy_box.isEmpty(),
                 "Mapped surface-layer ghost copy has no valid overlap with its source FAB.");
             ParallelFor(copy_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 dst(i,j,k) = src(i,j,k);
             });
         }
     }
     Gpu::streamSynchronize();
 }

Here is the call graph for this function:

◆ fill_pblh_from_coarser()

void SurfaceLayer::fill_pblh_from_coarser ( const int & lev )

Set the PBL height of a refined level whose grids do not run from the ground to the top of the domain in every column from the next coarser level, column by column (see compute_pblh). Public because it launches a device lambda.

Parameters

[in] lev Current level (> 0)

Set the PBL height of a refined level from the next coarser level: each column of this level takes the height of the coarse column that holds it. Used for a level whose grids do not run from the ground to the top of the domain in every column (see compute_pblh).

Parameters

[in] lev Current level (> 0)

 {
     AMREX_ALWAYS_ASSERT(lev > 0 && pblh[lev]);
     inject_pblh_from_coarser(lev, *pblh[lev]);
     fill_planar_boundary(lev, *pblh[lev]);
 }

Here is the call graph for this function:

◆ fill_planar_boundary()

void SurfaceLayer::fill_planar_boundary	(	const int &	lev,
		amrex::MultiFab &	mf
	)

Fill the ghost cells of a planar surface-layer MultiFab, and the valid region of its uncomputed copies when the 3D BoxArray is split in z (see PlanarBoundary).

Parameters

[in]	lev	level index
[in,out]	mf	planar MultiFab to fill

Fill the ghost cells of a planar surface-layer MultiFab.

The planar MultiFabs hold one box per 3D box, so a 3D BoxArray split in z gives duplicate planar boxes of which only the surface copy is computed; a FillBoundary could then fill a ghost cell from an uncomputed copy (see PlanarBoundary). With the split, the valid region of the uncomputed copies is filled as well. With EB terrain the fields are 3D and FillBoundary is well defined.

Parameters

[in]	lev	Current level
[in,out]	mf	Planar MultiFab to fill

 {
     // PlanarBoundary handles both z-low and z-high z-collapsed arrays split
     // in z. EB fields use ordinary FillBoundary because their 3-D layout has
     // no duplicate collapsed surface boxes. Lateral planar fields need a
     // selective exchange so an interior grid cannot provide wall data.
     if (m_terrain_type == TerrainType::EB) {
         mf.FillBoundary(m_geom[lev].periodicity());
     } else if (m_face.coordDir() != 2) {
         fill_lateral_surface_parameter_ghosts(lev, &mf);
     } else {
         m_planar_bndry[lev].fill(mf, m_geom[lev].periodicity());
     }
 }

◆ fill_qsurf_with_qsat()

void SurfaceLayer::fill_qsurf_with_qsat	(	const int &	lev,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Fill surface moisture from saturation specific humidity.

Parameters

[in]	lev	level index
[in]	cons_in	conserved state used to evaluate surface pressure
[in]	z_phys_nd	nodal physical-height field

Fill sea-surface moisture with saturation specific humidity.

Parameters

[in]	lev	Current level
[in]	cons_in	Conserved state used to derive pressure at the surface
[in]	z_phys_nd	Nodal physical height used to compute terrain-relative surface height

 {
     const int dir = m_face.coordDir();
     int sm_index;
     if (m_face.isLow()) {
         sm_index = m_geom[lev].Domain().smallEnd(dir);
     } else {
         sm_index = m_geom[lev].Domain().bigEnd(dir);
     }
  
     // Populate q_surf with qsat over water. The selected face is the only
     // authoritative slab; state/terrain ghosts are useful when valid, but a
     // bad halo must not turn into a collective conversion failure.
     const Real dz = m_geom[lev].CellSize(2);
     const int ng_z = amrex::min(cons_in.nGrowVect()[2],
                                 amrex::min(t_surf[lev]->nGrowVect()[2],
                                            q_surf[lev]->nGrowVect()[2]));
     const bool moist = use_moisture;
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     const Real rdOcp = m_rdOcp;
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
     // Use the 2-D surface mask as the iterator so ranks participate only when
     // their grids coincide with the selected face.
     for (MFIter mfi(*m_lmask_lev[lev][0]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
         Box tbx = mfi.validbox();
         const Box tilebx = mfi.tilebox();
  
         // Since lmask is used in the MFIter, its Z dimension is 0. These are
         // temporary geometry boxes, so use the physical domain's Z range for
         // indexing the 3-D surface/state arrays.
         tbx.setSmall(2, m_geom[lev].Domain().smallEnd(2));
         tbx.setBig(2, m_geom[lev].Domain().bigEnd(2));
         gtbx.setSmall(2, m_geom[lev].Domain().smallEnd(2));
         gtbx.setBig(2, m_geom[lev].Domain().bigEnd(2));
  
         if (dir == 2) {
             gtbx.makeSlab(2, sm_index);
         } else {
             gtbx.grow(2, ng_z);
             gtbx.setSmall(dir, sm_index);
             gtbx.setBig(dir, sm_index);
         }
  
         if (dir == 2) {
             if (m_terrain_type != TerrainType::EB &&
                 !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
                 continue;
             }
         } else {
             if (tbx[m_face] != sm_index || tilebx[m_face] != sm_index) {
                 continue;
             }
         }
  
         // The mask iterator can have more ghosts than the surface fields.
         // Limit every array accessed by the kernel to its corresponding FAB.
         gtbx &= t_surf[lev]->fabbox(mfi.index());
         gtbx &= q_surf[lev]->fabbox(mfi.index());
         gtbx &= cons_in.fabbox(mfi.index());
         if (z_phys_nd) {
             // z_phys_nd is nodal; the cell box whose corner nodes all lie in its FAB
             // is the nodal box converted to cells (one fewer cell on the high side).
             gtbx &= amrex::convert(z_phys_nd->fabbox(mfi.index()), IntVect::TheCellVector());
         }
         if (gtbx.isEmpty()) { continue; }
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto q_surf_arr = q_surf[lev]->array(mfi);
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
         const Box source_box = cons_in.boxArray()[mfi.index()];
         const int src_i_lo = source_box.smallEnd(0);
         const int src_i_hi = source_box.bigEnd(0);
         const int src_j_lo = source_box.smallEnd(1);
         const int src_j_hi = source_box.bigEnd(1);
         const int src_k_lo = source_box.smallEnd(2);
         const int src_k_hi = source_box.bigEnd(2);
         const bool z_face = (dir == 2);
         const bool low_face = m_face.isLow();
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
             if (!is_land) {
                 const bool authoritative = i >= src_i_lo && i <= src_i_hi &&
                                            j >= src_j_lo && j <= src_j_hi &&
                                            k >= src_k_lo && k <= src_k_hi;
                 const Real rho = cons_arr(i,j,k,Rho_comp);
                 const Real rho_theta = cons_arr(i,j,k,RhoTheta_comp);
                 if (!amrex::Math::isfinite(rho) || rho <= Real(0.0) ||
                     !amrex::Math::isfinite(rho_theta)) {
                     erf_surface_temperature::record_conversion_failure(fail, authoritative, erf_surface_temperature::ConversionError::bad_density,
                         i, j, k, rho, rho_theta);
                     return;
                 }
                 Real qv = Real(0.0);
                 if (moist) {
                     if (!have_rho_qv) {
                         erf_surface_temperature::record_conversion_failure(fail, authoritative, erf_surface_temperature::ConversionError::no_moisture_comp,
                             i, j, k);
                         return;
                     }
                     const Real rho_qv = cons_arr(i,j,k,RhoQ1_comp);
                     if (!amrex::Math::isfinite(rho_qv)) {
                         erf_surface_temperature::record_conversion_failure(fail, authoritative, erf_surface_temperature::ConversionError::bad_rho_qv,
                             i, j, k, rho, rho_theta, rho_qv);
                         return;
                     }
                     qv = rho_qv / rho;
                     if (!amrex::Math::isfinite(qv)) {
                         erf_surface_temperature::record_conversion_failure(fail, authoritative, erf_surface_temperature::ConversionError::bad_qv,
                             i, j, k, rho, rho_theta, rho_qv);
                         return;
                     }
                 }
                 Real delta_z = Real(0.0);
                 if (z_face) {
                     if (z_arr) {
                         const Real z_cc = Compute_Z_AtCellCenter(i,j,k,z_arr);
                         const Real z_face_local = low_face
                             ? Compute_Z_AtWFace(i,j,k,z_arr)
                             : Compute_Z_AtWFace(i,j,k+1,z_arr);
                         delta_z = z_cc - z_face_local;
                     } else {
                         delta_z = low_face ? myhalf*dz : -myhalf*dz;
                     }
                 }
                 const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                     rho, rho_theta, qv, delta_z);
                 Real t_surface = t_surf_arr(i,j,k);
                 if (erf_surface_temperature::theta_to_temperature(
                         t_surface, pressure, rdOcp, t_surface)) {
                     erf_qsatw(t_surface, pressure * Real(0.01), q_surf_arr(i,j,k));
                 } else {
                     erf_surface_temperature::record_conversion_failure(fail, authoritative, erf_surface_temperature::ConversionError::bad_theta_to_temperature,
                         i, j, k, rho, rho_theta,
                         moist ? cons_arr(i,j,k,RhoQ1_comp) : Real(0.0),
                         t_surface, pressure, delta_z);
                 }
             }
         });
     }
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure("SurfaceLayer fill_qsurf_with_qsat",
                                                          lev, d_fail.dataValue());
     fill_planar_boundary(lev, *q_surf[lev]);
 }

Here is the call graph for this function:

◆ fill_qsurf_with_skin_moisture()

void SurfaceLayer::fill_qsurf_with_skin_moisture	(	const int &	lev,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Overwrite the land surface mixing ratio from the external skin temperature and the moisture field of set_skin_moisture, with q_air the reference-height mixing ratio the flux uses. Without resistances: beta * q_sat(T_skin, p_s) + (1 - beta) * q_air (erf_surface_moisture::surface_mixing_ratio) with beta the soil-water factor, so the moisture flux is beta times the potential one. With resistances: erf_surface_moisture::two_source_surface_mixing_ratio of the aerodynamic resistance (from the last u* and Obukhov length; neutral before the first flux), the canopy resistance with its vapour-deficit factor, the soil resistance and the soil pore air's relative humidity. Needs set_skin_temperature and the averages of this step.

Parameters

[in]	lev	level index
[in]	cons_in	conserved state (surface pressure and q_air from the lowest cell)
[in]	z_phys_nd	nodal heights, or nullptr on a flat mesh

Overwrite the land surface mixing ratio with that of a surface of moisture availability beta at the external skin temperature (see the declaration). Like the skin temperature, written on the land cells of the surface slab, the halo taking the nearest valid column.

Parameters

[in]	lev	Current level
[in]	cons_in	Conserved state (surface pressure and q_air from the lowest cell)
[in]	z_phys_nd	Nodal heights, or nullptr on a flat mesh

 {
     using namespace erf_surface_moisture;
     const MultiFab& skin = *m_skin_tsurf_lev[lev];
     const MultiFab& moisture = *m_skin_moisture_lev[lev];
     // Indexed with the MFIter of q_surf, so the layouts must agree.
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         skin.boxArray() == q_surf[lev]->boxArray() &&
         skin.DistributionMap() == q_surf[lev]->DistributionMap() &&
         moisture.boxArray() == q_surf[lev]->boxArray() &&
         moisture.DistributionMap() == q_surf[lev]->DistributionMap() &&
         moisture.nComp() >= NumComponents,
         "Skin temperature or moisture layout does not match the surface-layer layout.");
  
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const Real dz = m_geom[lev].CellSize(2);
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     const bool use_resistances = m_skin_moisture_resistances_lev[lev] != 0;
     const Real hs = m_skin_moisture_hs_lev[lev];
     const auto *const qvm_ptr  = m_ma.get_average(lev, 4); // reference-height mixing ratio
     const auto *const umm_ptr  = m_ma.get_average(lev, 6); // horizontal velocity magnitude
     const auto *const zref_ptr = m_ma.get_zref(lev);       // reference height
     const similarity_funs sfuns{};
     constexpr Real wsmin = most_min_wind_speed; // the flux kernels' WSMIN
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
  
     for (MFIter mfi(*q_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
         gtbx &= q_surf[lev]->fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         const Box vbx = mfi.validbox();
         const int i_lo = vbx.smallEnd(0); const int i_hi = vbx.bigEnd(0);
         const int j_lo = vbx.smallEnd(1); const int j_hi = vbx.bigEnd(1);
  
         auto q_surf_arr = q_surf[lev]->array(mfi);
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
         const auto skin_arr = skin.const_array(mfi);
         const auto moist_arr = moisture.const_array(mfi);
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
         const auto qvm_arr  = qvm_ptr->const_array(mfi);
         const auto umm_arr  = umm_ptr->const_array(mfi);
         const auto zref_arr = zref_ptr->const_array(mfi);
         const auto z0_arr   = z_0[lev].const_array(mfi);
         const auto ustar_arr = u_star[lev]->const_array(mfi);
         const auto olen_arr = olen[lev]->const_array(mfi);
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const int li = amrex::min(amrex::max(i, i_lo), i_hi);
             const int lj = amrex::min(amrex::max(j, j_lo), j_hi);
             const int is_land = (lmask_arr) ? lmask_arr(li,lj,0) : 1;
             if (!is_land) { return; }
             if (!have_rho_qv) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::no_moisture_comp, li, lj, klo);
                 return;
             }
  
             // Surface pressure from the lowest cell, as for the skin temperature.
             const Real rho = cons_arr(li,lj,klo,Rho_comp);
             const Real qv_cell = cons_arr(li,lj,klo,RhoQ1_comp) / rho;
             const Real delta_z = z_arr
                 ? Compute_Z_AtCellCenter(li,lj,klo,z_arr) -
                   Compute_Z_AtWFace(li,lj,klo,z_arr)
                 : myhalf*dz;
             const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                 rho, cons_arr(li,lj,klo,RhoTheta_comp), qv_cell, delta_z);
             const Real t_skin = skin_arr(li,lj,k);
             const Real q_air = qvm_arr(li,lj,k);
             if (!amrex::Math::isfinite(qv_cell)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_qv, li, lj, klo,
                     rho, cons_arr(li,lj,klo,RhoTheta_comp), cons_arr(li,lj,klo,RhoQ1_comp));
                 return;
             }
             if (!amrex::Math::isfinite(q_air)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_q_air, li, lj, klo,
                     rho, cons_arr(li,lj,klo,RhoTheta_comp), cons_arr(li,lj,klo,RhoQ1_comp),
                     t_skin, pressure, delta_z, q_air);
                 return;
             }
             if (!erf_surface_temperature::valid_pressure(pressure)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_pressure, li, lj, klo,
                     rho, cons_arr(li,lj,klo,RhoTheta_comp), cons_arr(li,lj,klo,RhoQ1_comp),
                     t_skin, pressure, delta_z);
                 return;
             }
             if (!erf_surface_temperature::valid_temperature(t_skin)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_skin_temperature, li, lj, klo,
                     rho, cons_arr(li,lj,klo,RhoTheta_comp), cons_arr(li,lj,klo,RhoQ1_comp),
                     t_skin, pressure, delta_z);
                 return;
             }
             Real q_sat = Real(0.0);
             erf_qsatw(t_skin, pressure * Real(0.01), q_sat);
  
             Real q_surf_value = Real(0.0);
             if (use_resistances) {
                 // The aerodynamic resistance the moisture flux will see (the surface_temp
                 // kernel's): the last u* and Obukhov length, or neutral before the first flux.
                 // That neutral guess takes the plane-averaged wind, not the w*-enhanced wind
                 // the kernel uses with include_wstar, so on the first step of such a run in
                 // light wind it is not the kernel's r_a; it is only a first guess, and from
                 // then on u* carries the enhancement.
                 const Real r_a = surface_layer_aerodynamic_resistance(
                     sfuns, zref_arr(li,lj,k), z0_arr(li,lj,k), ustar_arr(li,lj,k),
                     olen_arr(li,lj,k), umm_arr(li,lj,k), wsmin, KAPPA, bogus_large_value);
                 const Real r_c = canopy_resistance(
                     moist_arr(li,lj,k,CanopyResistanceWithoutVPD), hs, q_sat, q_air);
                 q_surf_value = two_source_surface_mixing_ratio(
                     r_a, r_c, moist_arr(li,lj,k,SoilResistance),
                     moist_arr(li,lj,k,VegetationFraction), q_sat,
                     moist_arr(li,lj,k,GroundRelativeHumidity), q_air);
             } else {
                 q_surf_value = surface_mixing_ratio(moist_arr(li,lj,k,SoilFactor), q_sat, q_air);
             }
             if (!amrex::Math::isfinite(q_surf_value)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_q_surf, li, lj, klo,
                     rho, cons_arr(li,lj,klo,RhoTheta_comp), cons_arr(li,lj,klo,RhoQ1_comp),
                     t_skin, pressure, delta_z, q_air);
                 return;
             }
             q_surf_arr(i,j,k) = q_surf_value;
         });
     }
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure(
         "SurfaceLayer fill_qsurf_with_skin_moisture", lev, d_fail.dataValue());
 }

Here is the call graph for this function:

◆ fill_tsurf_with_abs_surf_temp()

void SurfaceLayer::fill_tsurf_with_abs_surf_temp	(	const int &	lev,
		const double &	time,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Set t_surf from the prescribed absolute surface temperature (most.abs_surf_temp, plus most.surf_heating_rate times the time), converted to the potential temperature the surface layer works in with the pressure at the surface diagnosed from the lowest cell, as for coupled SST. Like update_surf_temp, which it replaces in that case, it writes every cell of the surface.

Parameters

[in]	lev	level index
[in]	time	elapsed simulation time
[in]	cons_in	conserved state used to evaluate surface pressure
[in]	z_phys_nd	nodal physical-height field

Set t_surf from the prescribed absolute surface temperature, converted to potential temperature at the surface pressure.

Parameters

[in]	lev	Current level
[in]	time	Elapsed time, for the heating rate
[in]	cons_in	Conserved state (surface pressure from the lowest cell)
[in]	z_phys_nd	Nodal heights, or nullptr on a flat mesh

 {
     const Real t_abs = surf_temp + surf_heating_rate * static_cast<Real>(time);
     if (surf_heating_rate != 0) {
         amrex::Print() << "Surface temp (absolute) at t=" << time << ": " << t_abs << std::endl;
     }
  
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const Real dz = m_geom[lev].CellSize(2);
     const bool moist = use_moisture;
     const Real rdOcp = m_rdOcp;
     // As in fill_tsurf_with_coupled_sst: a moist run whose state lacks the vapour
     // component aborts below rather than reading past the components.
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
  
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
         gtbx &= t_surf[lev]->fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         // The halo takes the pressure of the nearest valid column; fill_planar_boundary
         // in update_fluxes then fills the interior and periodic ghosts.
         const Box vbx = mfi.validbox();
         const int i_lo = vbx.smallEnd(0); const int i_hi = vbx.bigEnd(0);
         const int j_lo = vbx.smallEnd(1); const int j_hi = vbx.bigEnd(1);
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const int li = amrex::min(amrex::max(i, i_lo), i_hi);
             const int lj = amrex::min(amrex::max(j, j_lo), j_hi);
  
             // The state is checked before it is divided by, as in fill_tsurf_with_coupled_sst,
             // so that a bad state reaches the diagnostic rather than a floating-point trap.
             const Real rho = cons_arr(li,lj,klo,Rho_comp);
             const Real rho_theta = cons_arr(li,lj,klo,RhoTheta_comp);
             if (!amrex::Math::isfinite(rho) || rho <= Real(0.0) ||
                 !amrex::Math::isfinite(rho_theta)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_density,
                     li, lj, klo, rho, rho_theta);
                 return;
             }
             Real rho_qv = Real(0.0);
             Real qv = Real(0.0);
             if (moist) {
                 if (!have_rho_qv) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::no_moisture_comp,
                         li, lj, klo);
                     return;
                 }
                 rho_qv = cons_arr(li,lj,klo,RhoQ1_comp);
                 if (!amrex::Math::isfinite(rho_qv)) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_rho_qv,
                         li, lj, klo, rho, rho_theta, rho_qv);
                     return;
                 }
                 qv = rho_qv / rho;
                 if (!amrex::Math::isfinite(qv)) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_qv,
                         li, lj, klo, rho, rho_theta, rho_qv);
                     return;
                 }
             }
             const Real delta_z = z_arr
                 ? Compute_Z_AtCellCenter(li,lj,klo,z_arr) -
                   Compute_Z_AtWFace(li,lj,klo,z_arr)
                 : myhalf*dz;
             const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                 rho, rho_theta, qv, delta_z);
             Real theta = t_surf_arr(i,j,k);
             if (erf_surface_temperature::temperature_to_theta(
                     t_abs, pressure, rdOcp, theta)) {
                 t_surf_arr(i,j,k) = theta;
             } else {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_temperature_to_theta,
                     li, lj, klo, rho, rho_theta, rho_qv, t_abs, pressure, delta_z);
             }
         });
     }
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure(
         "SurfaceLayer fill_tsurf_with_abs_surf_temp", lev, d_fail.dataValue());
 }

Here is the call graph for this function:

◆ fill_tsurf_with_coupled_sst()

void SurfaceLayer::fill_tsurf_with_coupled_sst	(	const int &	lev,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Overwrite surface temperature with coupled ocean SST where the coupler covers the cell.

Runs after fill_tsurf_with_sst_and_tsk so that the lower-boundary data is the base layer: land cells, and water cells with no ocean donor, keep the value written there.

Parameters

[in] lev level index

Overwrite surface temperature with coupled ocean SST where covered.

Parameters

[in] lev Current level

 {
     // No coupler has handed us anything yet. Whatever fill_tsurf_with_sst_and_tsk
     // wrote stands, which is the correct answer for one-way and uncoupled runs.
     if (m_coupled_sst_lev.empty() || !m_coupled_sst_lev[lev]) { return; }
  
     // The loop below iterates t_surf and indexes the coupled arrays with the
     // same MFIter, so the layouts must agree. They do for planar terrain, where
     // t_surf is grids[lev] flattened with setRange(2,0) -- the same construction
     // GetOceanToAtmosSurfaceLayout reports. Under EB terrain t_surf keeps the
     // full 3D BoxArray and they would not, so fail loudly rather than read the
     // wrong fab.
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         m_coupled_sst_lev[lev]->boxArray()       == t_surf[lev]->boxArray() &&
         m_coupled_sst_lev[lev]->DistributionMap() == t_surf[lev]->DistributionMap(),
         "Coupled SST layout does not match the surface-layer layout.");
  
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const Real dz = m_geom[lev].CellSize(2);
     const bool moist = use_moisture;
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     const Real rdOcp = m_rdOcp;
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
  
     // Absent coverage information we must assume nothing is covered: silently
     // treating the whole field as valid is how an uncovered cell ends up holding
     // the remap's zero fill.
     const bool has_valid = (m_coupled_sst_valid_lev[lev] != nullptr);
  
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
  
         // NOTE: the coupled lane does not carry lateral ghost cells, so clamp
         //       into the physical source box exactly as the fallback path does.
         //       FillBoundary in update_fluxes picks up the interior and
         //       periodic directions.
         const Box source_box = cons_in.boxArray()[mfi.index()];
         const Box donor_box = m_coupled_sst_lev[lev]->boxArray()[mfi.index()];
         const int source_i_lo = source_box.smallEnd(0);
         const int source_i_hi = source_box.bigEnd(0);
         const int source_j_lo = source_box.smallEnd(1);
         const int source_j_hi = source_box.bigEnd(1);
         const int donor_i_lo = donor_box.smallEnd(0);
         const int donor_i_hi = donor_box.bigEnd(0);
         const int donor_j_lo = donor_box.smallEnd(1);
         const int donor_j_hi = donor_box.bigEnd(1);
         const int donor_k = donor_box.smallEnd(2);
         gtbx &= t_surf[lev]->fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
         const auto coupled_sst_arr = m_coupled_sst_lev[lev]->const_array(mfi);
         auto const& valid_arr = has_valid ? m_coupled_sst_valid_lev[lev]->const_array(mfi)
                                           : Array4<const int>{};
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const int li = amrex::min(amrex::max(i, source_i_lo), source_i_hi);
             const int lj = amrex::min(amrex::max(j, source_j_lo), source_j_hi);
             const int si = amrex::min(amrex::max(li, donor_i_lo), donor_i_hi);
             const int sj = amrex::min(amrex::max(lj, donor_j_lo), donor_j_hi);
             int is_land = (lmask_arr) ? lmask_arr(li,lj,0) : 1;
             if (is_land) { return; }
  
             if (!has_valid || valid_arr(si,sj,donor_k) == 0) { return; }
  
             const Real rho = cons_arr(li,lj,klo,Rho_comp);
             const Real rho_theta = cons_arr(li,lj,klo,RhoTheta_comp);
             if (!amrex::Math::isfinite(rho) || rho <= Real(0.0) ||
                 !amrex::Math::isfinite(rho_theta)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_density,
                     li, lj, klo, rho, rho_theta);
                 return;
             }
             Real qv = Real(0.0);
             if (moist) {
                 if (!have_rho_qv) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::no_moisture_comp,
                         li, lj, klo);
                     return;
                 }
                 const Real rho_qv = cons_arr(li,lj,klo,RhoQ1_comp);
                 if (!amrex::Math::isfinite(rho_qv)) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_rho_qv,
                         li, lj, klo, rho, rho_theta, rho_qv);
                     return;
                 }
                 qv = rho_qv / rho;
                 if (!amrex::Math::isfinite(qv)) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_qv,
                         li, lj, klo, rho, rho_theta, rho_qv);
                     return;
                 }
             }
             const Real delta_z = z_arr
                 ? Compute_Z_AtCellCenter(li,lj,klo,z_arr) -
                   Compute_Z_AtWFace(li,lj,klo,z_arr)
                 : myhalf*dz;
             const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                 rho, rho_theta, qv, delta_z);
             Real theta = t_surf_arr(i,j,k);
             if (erf_surface_temperature::temperature_to_theta(
                     coupled_sst_arr(si,sj,donor_k), pressure, rdOcp, theta)) {
                 t_surf_arr(i,j,k) = theta;
             } else {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_temperature_to_theta,
                     li, lj, klo, rho, rho_theta,
                     moist ? cons_arr(li,lj,klo,RhoQ1_comp) : Real(0.0),
                     coupled_sst_arr(si,sj,donor_k), pressure, delta_z);
             }
         });
     }
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure("SurfaceLayer fill_tsurf_with_coupled_sst",
                                                          lev, d_fail.dataValue());
 }

Here is the call graph for this function:

◆ fill_tsurf_with_sfc_sst()

void SurfaceLayer::fill_tsurf_with_sfc_sst	(	const int &	lev,
		const double &	time,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Fill surface temperature interpolated from time varying SST file

Parameters

[in]	lev	level index
[in]	time	interpolation time

 {
     update_sfc_time_index(elapsed_time);
     const Real sfc_sst = interpolate_sfc_column(elapsed_time, 1);
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const Real dz = m_geom[lev].CellSize(2);
     const bool moist = use_moisture;
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     const Real rdOcp = m_rdOcp;
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
  
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         // A z-split planar field has one copy for every stacked 3-D box. Only
         // the copy whose source box touches the physical z-low face owns the
         // text-SST conversion and may read the 3-D state or terrain.
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
  
         const Box source_box = cons_in.boxArray()[mfi.index()];
         gtbx &= t_surf[lev]->fabbox(mfi.index());
         gtbx &= cons_in.fabbox(mfi.index());
         if (z_phys_nd) {
             gtbx &= amrex::convert(z_phys_nd->fabbox(mfi.index()), IntVect::TheCellVector());
         }
         if (m_lmask_lev[lev][0]) {
             gtbx &= m_lmask_lev[lev][0]->fabbox(mfi.index());
         }
         if (gtbx.isEmpty()) { continue; }
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
         const int i_lo = source_box.smallEnd(0);
         const int i_hi = source_box.bigEnd(0);
         const int j_lo = source_box.smallEnd(1);
         const int j_hi = source_box.bigEnd(1);
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const int li = amrex::min(amrex::max(i, i_lo), i_hi);
             const int lj = amrex::min(amrex::max(j, j_lo), j_hi);
             int is_land = (lmask_arr) ? lmask_arr(li,lj,0) : 0;
             if (!is_land) {
                 const Real rho = cons_arr(li,lj,klo,Rho_comp);
                 const Real rho_theta = cons_arr(li,lj,klo,RhoTheta_comp);
                 if (!amrex::Math::isfinite(rho) || rho <= Real(0.0) ||
                     !amrex::Math::isfinite(rho_theta)) {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_density,
                         li, lj, klo, rho, rho_theta);
                     return;
                 }
                 Real qv = Real(0.0);
                 if (moist) {
                     if (!have_rho_qv) {
                         erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::no_moisture_comp,
                             li, lj, klo);
                         return;
                     }
                     const Real rho_qv = cons_arr(li,lj,klo,RhoQ1_comp);
                     if (!amrex::Math::isfinite(rho_qv)) {
                         erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_rho_qv,
                             li, lj, klo, rho, rho_theta, rho_qv);
                         return;
                     }
                     qv = rho_qv / rho;
                     if (!amrex::Math::isfinite(qv)) {
                         erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_qv,
                             li, lj, klo, rho, rho_theta, rho_qv);
                         return;
                     }
                 }
                 const Real delta_z = z_arr
                     ? Compute_Z_AtCellCenter(li,lj,klo,z_arr) -
                       Compute_Z_AtWFace(li,lj,klo,z_arr)
                     : myhalf*dz;
                 const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                     rho, rho_theta, qv, delta_z);
                 Real theta = t_surf_arr(i,j,k);
                 if (erf_surface_temperature::temperature_to_theta(sfc_sst, pressure, rdOcp, theta)) {
                     t_surf_arr(i,j,k) = theta;
                 } else {
                     erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_temperature_to_theta,
                         li, lj, klo, rho, rho_theta,
                         moist ? cons_arr(li,lj,klo,RhoQ1_comp) : Real(0.0),
                         sfc_sst, pressure, delta_z);
                 }
             }
         });
     }
  
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure("SurfaceLayer fill_tsurf_with_sfc_sst",
                                                          lev, d_fail.dataValue());
  
     fill_planar_boundary(lev, *t_surf[lev]);
 }

Here is the call graph for this function:

◆ fill_tsurf_with_skin_temperature()

void SurfaceLayer::fill_tsurf_with_skin_temperature	(	const int &	lev,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Overwrite the land surface temperature with an external skin temperature.

The skin is an absolute temperature [K]; it is converted to the potential temperature the surface layer works in with the pressure at the surface diagnosed from the lowest cell, as for coupled SST. Runs last among the surface-temperature writers, so on land it replaces what they wrote.

Parameters

[in]	lev	level index
[in]	cons_in	conserved state used to evaluate surface pressure
[in]	z_phys_nd	nodal physical-height field

Overwrite the land surface temperature with an external absolute skin temperature, converted to potential temperature at the surface pressure.

Parameters

[in]	lev	Current level
[in]	cons_in	Conserved state (surface pressure from the lowest cell)
[in]	z_phys_nd	Nodal heights, or nullptr on a flat mesh

 {
     const MultiFab& skin = *m_skin_tsurf_lev[lev];
     // Indexed with the MFIter of t_surf, so the layouts must agree: both are the
     // level's grids flattened to the surface for planar terrain.
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         skin.boxArray() == t_surf[lev]->boxArray() &&
         skin.DistributionMap() == t_surf[lev]->DistributionMap(),
         "Skin temperature layout does not match the surface-layer layout.");
  
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const Real dz = m_geom[lev].CellSize(2);
     const bool moist = use_moisture;
     const Real rdOcp = m_rdOcp;
     // As in fill_tsurf_with_coupled_sst: a moist run whose state lacks the vapour
     // component aborts below rather than reading past the components.
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
  
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
         gtbx &= t_surf[lev]->fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         // The skin is written on valid cells only, so the halo takes the nearest
         // valid column; fill_planar_boundary in update_fluxes then fills the
         // interior and periodic ghosts.
         const Box vbx = mfi.validbox();
         const int i_lo = vbx.smallEnd(0); const int i_hi = vbx.bigEnd(0);
         const int j_lo = vbx.smallEnd(1); const int j_hi = vbx.bigEnd(1);
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
         const auto skin_arr = skin.const_array(mfi);
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const int li = amrex::min(amrex::max(i, i_lo), i_hi);
             const int lj = amrex::min(amrex::max(j, j_lo), j_hi);
             const int is_land = (lmask_arr) ? lmask_arr(li,lj,0) : 1;
             if (!is_land) { return; }
  
             const Real rho = cons_arr(li,lj,klo,Rho_comp);
             const Real rho_theta = cons_arr(li,lj,klo,RhoTheta_comp);
             if (moist && !have_rho_qv) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::no_moisture_comp, li, lj, klo);
                 return;
             }
             const Real qv = moist ? cons_arr(li,lj,klo,RhoQ1_comp) / rho : Real(0.0);
             const Real delta_z = z_arr
                 ? Compute_Z_AtCellCenter(li,lj,klo,z_arr) -
                   Compute_Z_AtWFace(li,lj,klo,z_arr)
                 : myhalf*dz;
             const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                 rho, rho_theta, qv, delta_z);
             Real theta = t_surf_arr(i,j,k);
             if (!amrex::Math::isfinite(qv)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_qv, li, lj, klo, rho, rho_theta,
                     moist ? cons_arr(li,lj,klo,RhoQ1_comp) : Real(0.0));
             } else if (erf_surface_temperature::temperature_to_theta(
                            skin_arr(li,lj,k), pressure, rdOcp, theta)) {
                 t_surf_arr(i,j,k) = theta;
             } else {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_temperature_to_theta,
                     li, lj, klo, rho, rho_theta,
                     moist ? cons_arr(li,lj,klo,RhoQ1_comp) : Real(0.0),
                     skin_arr(li,lj,k), pressure, delta_z);
             }
         });
     }
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure(
         "SurfaceLayer fill_tsurf_with_skin_temperature", lev, d_fail.dataValue());
 }

Here is the call graph for this function:

◆ fill_tsurf_with_sst_and_tsk()

void SurfaceLayer::fill_tsurf_with_sst_and_tsk	(	const int &	lev,
		const double &	elapsed_time_since_start_low
	)

Fill surface temperature from available SST and TSK data.

Parameters

[in]	lev	level index
[in]	time	interpolation time

Fill surface temperature from SST/TSK lower-boundary data.

Parameters

[in]	lev	Current level
[in]	elapsed_time_since_start_low	Time since the start of the lower-boundary data

 {
     int n_times_in_sst = static_cast<int>(m_sst_lev[lev].size());
  
     double dT = m_low_time_interval;
  
     int n_time_lo, n_time_hi;
     Real alpha;
  
     if (n_times_in_sst > 1) {
         n_time_lo = static_cast<int>( elapsed_time_since_start_low /  dT);
         alpha = static_cast<Real>((elapsed_time_since_start_low - n_time_lo * dT) / dT);
  
         AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one);
  
         n_time_hi = n_time_lo + 1;
  
         // Do not over run the last sst file
         if (m_start_low_time + elapsed_time_since_start_low >= m_final_low_time) {
             n_time_lo = static_cast<int>(m_sst_lev[lev].size())-1;
             n_time_hi = n_time_lo;
             alpha     = zero;
         }
  
         AMREX_ALWAYS_ASSERT( (n_time_lo >= 0) && (n_time_hi < m_sst_lev[lev].size()));
     } else {
         n_time_lo = 0;
         n_time_hi = 0;
         alpha     = one;
     }
     AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one);
  
     Real oma   = one - alpha;
  
     // Define a default land surface temperature if we don't read in tsk
     Real lst = default_land_surf_temp;
  
     bool use_tsk = (m_tsk_lev[lev][0]);
     bool ignore_sst = m_ignore_sst;
  
     const int klo = m_geom[lev].Domain().smallEnd(2);
  
     // Populate t_surf
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         if (gtbx.smallEnd(2) != klo) { continue; }
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
  
         const auto sst_lo_arr = m_sst_lev[lev][n_time_lo]->const_array(mfi);
         const auto sst_hi_arr = m_sst_lev[lev][n_time_hi]->const_array(mfi);
  
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
  
         if (use_tsk) {
             // TSK carries the SST from the low file over water (update_sst_tsk), so it is
             // interpolated in time between the same two slices as the SST.  Written as
             // lo + alpha * (hi - lo) so that a cell whose TSK is the same in both slices,
             // as it is over land, keeps exactly that value.
             const auto tsk_lo_arr = m_tsk_lev[lev][n_time_lo]->const_array(mfi);
             const auto tsk_hi_arr = (m_tsk_lev[lev][n_time_hi]) ? m_tsk_lev[lev][n_time_hi]->const_array(mfi)
                                                                 : tsk_lo_arr;
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
                 if (!is_land && !ignore_sst) {
                     t_surf_arr(i,j,k) = oma   * sst_lo_arr(i,j,k)
                                       + alpha * sst_hi_arr(i,j,k);
                 } else {
                     t_surf_arr(i,j,k) = tsk_lo_arr(i,j,k)
                                       + alpha * (tsk_hi_arr(i,j,k) - tsk_lo_arr(i,j,k));
                 }
             });
         } else {
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
                 if (!is_land) {
                     t_surf_arr(i,j,k) = oma   * sst_lo_arr(i,j,k)
                                       + alpha * sst_hi_arr(i,j,k);
                 } else {
                     t_surf_arr(i,j,k) = lst;
                 }
             });
         }
     }
     fill_planar_boundary(lev, *t_surf[lev]);
 }

Here is the call graph for this function:

◆ get_lmask()

amrex::iMultiFab* SurfaceLayer::get_lmask ( const int & lev )

inline

Return the land-mask field.

Parameters

[in] lev level index

1639 { return m_lmask_lev[lev][0]; }

◆ get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf	(	const int &	lev,
		const amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd
	)

Fill surface temperature from the LSM surface-temperature field, converted to potential temperature at the surface pressure when the field is an absolute temperature.

Parameters

[in]	lev	level index
[in]	cons_in	conserved state (surface pressure from the lowest cell)
[in]	z_phys_nd	nodal heights, or nullptr on a flat mesh

Fill surface temperature from land-surface-model data.

By default the LSM field is copied into t_surf unchanged. With most.use_lsm_skin_temperature, t_surf is a potential temperature: an LSM skin temperature that is an absolute temperature (m_lsm_tsurf_is_theta false) is converted at the surface pressure, as fill_tsurf_with_skin_temperature does. When the LSM also supplies the surface fluxes, t_surf follows the skin only where all four fluxes are valid: those are the cells whose fluxes and u*, t*, q*, L come from the LSM, so the MOST fluxes of the remaining cells keep the surface temperature they had.

Parameters

[in]	lev	Current level
[in]	cons_in	Conserved state (surface pressure from the lowest cell)
[in]	z_phys_nd	Nodal heights, or nullptr on a flat mesh

 {
     const int klo = m_geom[lev].Domain().smallEnd(2);
  
     if (!m_use_lsm_skin_temperature) {
         // Default: the LSM field is copied into t_surf as it is.
         for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
         {
             Box gtbx = mfi.growntilebox();
  
             if (gtbx.smallEnd(2) != klo) { continue; }
  
             // NOTE: LSM does not carry lateral ghost cells.
             //       This copies the valid box into the ghost cells.
             //       Fillboundary is called after this to pick up the
             //       interior ghost and periodic directions.
             const NearestColumn col(mfi.validbox());
  
             auto t_surf_arr = t_surf[lev]->array(mfi);
             auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                       Array4<int> {};
             const auto& lsm_tsurf = *m_lsm_data_lev[lev][m_lsm_tsurf_indx];
             // get the top-most index of the LSM to use as the surface temperature
             const auto& lsm_box = lsm_tsurf.box(mfi.index());
             const int lsm_khi = lsm_box.bigEnd(2);
             const auto lsm_arr = lsm_tsurf.const_array(mfi);
  
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
                 if (is_land) {
                     const Real lsm_value = lsm_arr(col.i(i),col.j(j),lsm_khi);
                     if (amrex::Math::isfinite(lsm_value) && lsm_value > Real(0.0) &&
                         lsm_value < lsm_undefined) {
                         t_surf_arr(i,j,k) = lsm_value;
                     }
                 }
             });
         }
         return;
     }
  
     const Real dz = m_geom[lev].CellSize(2);
     const bool moist = use_moisture;
     const bool have_rho_qv = cons_in.nComp() > RhoQ1_comp;
     const Real rdOcp = m_rdOcp;
     const bool convert = !m_lsm_tsurf_is_theta;
     const bool gate_on_fluxes = m_has_lsm_fluxes;
     const bool use_surface_model_fluxes = use_surface_model && surf_model_fluxes;
     amrex::Gpu::DeviceScalar<erf_surface_temperature::ConversionFailure>
         d_fail(erf_surface_temperature::ConversionFailure{});
     auto* fail = d_fail.dataPtr();
  
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         // Only the copy whose source box touches the z-low face may read the
         // state there; fill_planar_boundary in update_fluxes fills the others.
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
         gtbx &= t_surf[lev]->fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         // NOTE: LSM does not carry lateral ghost cells.
         //       This copies the valid box into the ghost cells.
         //       Fillboundary is called after this to pick up the
         //       interior ghost and periodic directions.
         const NearestColumn col(mfi.validbox());
  
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                   Array4<int> {};
         const auto& lsm_tsurf = *m_lsm_data_lev[lev][m_lsm_tsurf_indx];
         // get the top-most index of the LSM to use as the surface temperature
         // this is -1 for SLM, but could be different for other models?
         const auto& lsm_box = lsm_tsurf.box(mfi.index());
         const int lsm_khi = lsm_box.bigEnd(2);
         const auto lsm_arr = lsm_tsurf.const_array(mfi);
         const auto cons_arr = cons_in.const_array(mfi);
         const auto z_arr    = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
                                             Array4<const Real> {};
  
         // The same four fluxes, and the same validity test, as
         // compute_sfc_params_from_lsm_fluxes
         auto lsm_t_flux_arr = Array4<const Real> {};
         auto lsm_q_flux_arr = Array4<const Real> {};
         auto lsm_tau13_arr  = Array4<const Real> {};
         auto lsm_tau23_arr  = Array4<const Real> {};
         if (gate_on_fluxes) {
             if (use_surface_model_fluxes) {
                 const SurfaceFluxView surf_flux = m_surf_model->get_surface_flux_view(lev, mfi);
                 lsm_t_flux_arr = surf_flux.t_flux;
                 lsm_q_flux_arr = surf_flux.q_flux;
                 lsm_tau13_arr  = surf_flux.tau13;
                 lsm_tau23_arr  = surf_flux.tau23;
             } else {
                 for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
                     if (toLower(m_lsm_flux_name[n]) == "t_flux") { lsm_t_flux_arr = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                     if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                     if (toLower(m_lsm_flux_name[n]) == "tau13")  { lsm_tau13_arr  = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                     if (toLower(m_lsm_flux_name[n]) == "tau23")  { lsm_tau23_arr  = m_lsm_flux_lev[lev][n]->const_array(mfi); }
                 }
             }
         }
         const bool have_fluxes = lsm_t_flux_arr && lsm_q_flux_arr && lsm_tau13_arr && lsm_tau23_arr;
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const int li = col.i(i);
             const int lj = col.j(j);
             int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
             if (!is_land) { return; }
  
             const Real lsm_value = lsm_arr(li,lj,lsm_khi);
             if (!amrex::Math::isfinite(lsm_value) || lsm_value <= Real(0.0) ||
                 lsm_value >= lsm_undefined) {
                 return;
             }
  
             if (gate_on_fluxes) {
                 if (!have_fluxes ||
                     !is_valid_lsm_value(lsm_t_flux_arr(li,lj,0)) ||
                     !is_valid_lsm_value(lsm_q_flux_arr(li,lj,0)) ||
                     !is_valid_lsm_value(lsm_tau13_arr (li,lj,0)) ||
                     !is_valid_lsm_value(lsm_tau23_arr (li,lj,0))) {
                     return;
                 }
             }
  
             if (!convert) {
                 t_surf_arr(i,j,k) = lsm_value;
                 return;
             }
  
             const Real rho = cons_arr(li,lj,klo,Rho_comp);
             const Real rho_theta = cons_arr(li,lj,klo,RhoTheta_comp);
             if (moist && !have_rho_qv) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::no_moisture_comp, li, lj, klo);
                 return;
             }
             const Real qv = moist ? cons_arr(li,lj,klo,RhoQ1_comp) / rho : Real(0.0);
             const Real delta_z = z_arr
                 ? Compute_Z_AtCellCenter(li,lj,klo,z_arr) -
                   Compute_Z_AtWFace(li,lj,klo,z_arr)
                 : myhalf*dz;
             const Real pressure = erf_surface_temperature::pressure_at_boundary_from_cell(
                 rho, rho_theta, qv, delta_z);
             Real theta = t_surf_arr(i,j,k);
             if (!amrex::Math::isfinite(qv)) {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_qv, li, lj, klo, rho, rho_theta,
                     moist ? cons_arr(li,lj,klo,RhoQ1_comp) : Real(0.0));
             } else if (erf_surface_temperature::temperature_to_theta(
                            lsm_value, pressure, rdOcp, theta)) {
                 t_surf_arr(i,j,k) = theta;
             } else {
                 erf_surface_temperature::record_conversion_failure(fail, true, erf_surface_temperature::ConversionError::bad_temperature_to_theta,
                     li, lj, klo, rho, rho_theta,
                     moist ? cons_arr(li,lj,klo,RhoQ1_comp) : Real(0.0),
                     lsm_value, pressure, delta_z);
             }
         });
     }
     amrex::Gpu::streamSynchronize();
     erf_surface_temperature::abort_on_conversion_failure("SurfaceLayer get_lsm_tsurf",
                                                          lev, d_fail.dataValue());
 }

Here is the call graph for this function:

◆ get_mac_avg()

const amrex::MultiFab* SurfaceLayer::get_mac_avg	(	const int &	lev,
		int	comp
	)

inline

Return a MOST-average field.

Parameters

[in]	lev	level index
[in]	comp	component index

   {
       return m_ma.get_average(lev, comp);
   }

Here is the call graph for this function:

◆ get_mac_avg_ptr()

amrex::MultiFab* SurfaceLayer::get_mac_avg_ptr	(	const int &	lev,
		int	comp
	)

inline

Return a MOST-average field for modification (restart).

Parameters

[in]	lev	level index
[in]	comp	component index

1420 { return m_ma.get_average(lev, comp); }

Here is the call graph for this function:

◆ get_mac_plane_avg()

amrex::Vector<amrex::Real> SurfaceLayer::get_mac_plane_avg ( const int & lev ) const

inline

Return the filtered plane averages, which hold the filter state for the plane and EB averaging policies.

Parameters

[in] lev level index

1428 { return m_ma.get_plane_average(lev); }

MOSTAverage::get_plane_average

amrex::Vector< amrex::Real > get_plane_average(const int &lev) const

Definition: ERF_MOSTAverage.H:301

Here is the call graph for this function:

◆ get_num_mac_avg()

int SurfaceLayer::get_num_mac_avg ( ) const

inline

Return the number of MOST-average components.

1398 { return m_ma.get_navg(); }

MOSTAverage::get_navg

int get_navg() const

Definition: ERF_MOSTAverage.H:270

Here is the call graph for this function:

◆ get_olen()

amrex::MultiFab* SurfaceLayer::get_olen ( const int & lev )

inline

Return the Obukhov length field.

Parameters

[in] lev level index

1370 { return olen[lev].get(); }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_pblh()

amrex::MultiFab* SurfaceLayer::get_pblh ( const int & lev )

inline

Return the planetary-boundary-layer-height field.

Parameters

[in] lev level index

1377 { return pblh[lev].get(); }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_q_star()

amrex::MultiFab* SurfaceLayer::get_q_star ( const int & lev )

inline

Return the moisture scale field.

Parameters

[in] lev level index

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

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_q_surf()

amrex::MultiFab* SurfaceLayer::get_q_surf ( const int & lev )

inline

Return the surface-moisture field.

Parameters

[in] lev level index

1591 { return q_surf[lev].get(); }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_surface_diagnostic_source()

amrex::MultiFab* SurfaceLayer::get_surface_diagnostic_source ( const int & lev )

inline

Return the surface-diagnostic provenance field.

Parameters

[in] lev level index

1606 { return surface_diagnostic_source[lev].get(); }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_t_star()

amrex::MultiFab* SurfaceLayer::get_t_star ( const int & lev )

inline

Return the temperature scale field.

Parameters

[in] lev level index

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

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_t_surf()

amrex::MultiFab* SurfaceLayer::get_t_surf ( const int & lev )

inline

Return the surface-temperature field.

Parameters

[in] lev level index

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

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_u_star()

amrex::MultiFab* SurfaceLayer::get_u_star ( const int & lev )

inline

Return the friction-velocity field.

Parameters

[in] lev level index

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

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_w_star()

amrex::MultiFab* SurfaceLayer::get_w_star ( const int & lev )

inline

Return the convective velocity scale field.

Parameters

[in] lev level index

1337 { return w_star[lev].get(); }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_z0()

amrex::MultiFab* SurfaceLayer::get_z0 ( const int & lev )

inline

Return the roughness-height field.

Parameters

[in] lev level index

1627 { return &z_0[lev]; }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

◆ get_zref()

amrex::MultiFab* SurfaceLayer::get_zref ( const int & lev )

inline

Return the 2D MF of zref heights.

Parameters

[in] lev level index

1613 { return m_ma.get_zref(lev); }

Here is the call graph for this function:

◆ get_zref_min()

amrex::Real SurfaceLayer::get_zref_min ( const int & lev )

inline

Return the minimum reference height for one level.

Parameters

[in] lev level index

1620 { return (m_ma.get_zref(lev))->min(0); }

Here is the call graph for this function:

◆ have_variable_sea_roughness()

bool SurfaceLayer::have_variable_sea_roughness ( )

inline

Return whether variable sea roughness is active.

1632 { return m_var_z0; }

SurfaceLayer::m_var_z0

bool m_var_z0

Definition: ERF_SurfaceLayer.H:1832

◆ impose_SurfaceLayer_bcs()

void SurfaceLayer::impose_SurfaceLayer_bcs	(	const int &	lev,
		amrex::Vector< const amrex::MultiFab * >	mfs,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	Tau_lev,
		amrex::MultiFab *	xheat_flux,
		amrex::MultiFab *	yheat_flux,
		amrex::MultiFab *	zheat_flux,
		amrex::MultiFab *	xqv_flux,
		amrex::MultiFab *	yqv_flux,
		amrex::MultiFab *	zqv_flux,
		const amrex::MultiFab *	z_phys,
		amrex::MultiFab *	mf_surface_fluxes = `nullptr`
	)

Impose surface-layer boundary conditions for planar terrain.

Parameters

[in]	lev	level index
[in]	mfs	state and velocity fields used by the BC computation
[in,out]	Tau_lev	stress fields to fill
[in,out]	xheat_flux	x-face heat flux field
[in,out]	yheat_flux	y-face heat flux field
[in,out]	zheat_flux	z-face heat flux field
[in,out]	xqv_flux	x-face moisture flux field
[in,out]	yqv_flux	y-face moisture flux field
[in,out]	zqv_flux	z-face moisture flux field
[in]	z_phys	physical-height field

Wrapper to impose Monin Obukhov similarity theory fluxes by populating ghost cells.

Parameters

[in]	lev	Current level
[in]	mfs	State MultiFabs used to compute the boundary fluxes
[in,out]	Tau_lev	Diffusive stress MultiFabs populated with surface stresses
[in,out]	xheat_flux	x-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]	yheat_flux	y-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]	zheat_flux	z-face heat-flux MultiFab populated with vertical surface heat flux
[in,out]	xqv_flux	x-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]	yqv_flux	y-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]	zqv_flux	z-face moisture-flux MultiFab populated when moisture is enabled
[in]	z_phys	Nodal physical height used to rotate terrain-following fluxes

 {
     if (flux_type == FluxCalcType::MOENG) {
         amrex::Real wsmin = most_min_wind_speed; // TODO: change for different faces
         const Box& domain = m_geom[lev].Domain();
         moeng_flux flux_comp(wsmin, m_face.isLow(),
                              domain.smallEnd(2), domain.bigEnd(2));
         compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                  xheat_flux, yheat_flux, zheat_flux,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  z_phys, flux_comp, mf_surface_fluxes);
     } else if (flux_type == FluxCalcType::ROTATE) {
         rotate_flux flux_comp;
         compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                  xheat_flux, yheat_flux, zheat_flux,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  z_phys, flux_comp, mf_surface_fluxes);
     } else if (flux_type == FluxCalcType::RICO) {
         rico_flux flux_comp(rico_theta_z0, rico_qsat_z0);
         compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                  xheat_flux, yheat_flux, zheat_flux,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  z_phys, flux_comp, mf_surface_fluxes);
     } else if (flux_type == FluxCalcType::BULK_COEFF) {
         bulk_coeff_flux flux_comp(m_Cd, m_Ch, m_Cq);
         compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                  xheat_flux, yheat_flux, zheat_flux,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  z_phys, flux_comp, mf_surface_fluxes);
     } else if (flux_type == FluxCalcType::CUSTOM) {
         custom_flux flux_comp(specified_rho_surf);
         compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                  xheat_flux, yheat_flux, zheat_flux,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  z_phys, flux_comp, mf_surface_fluxes);
     } else {
         amrex::Abort("Unknown surface layer flux calculation type");
     }
 }

◆ impose_SurfaceLayer_bcs_EB()

void SurfaceLayer::impose_SurfaceLayer_bcs_EB	(	const int &	lev,
		amrex::Vector< const amrex::MultiFab * >	mfs,
		amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab >>> &	Tau_lev,
		amrex::MultiFab *	xheat_flux,
		amrex::MultiFab *	yheat_flux,
		amrex::MultiFab *	zheat_flux,
		amrex::MultiFab *	xqv_flux,
		amrex::MultiFab *	yqv_flux,
		amrex::MultiFab *	zqv_flux
	)

Impose surface-layer boundary conditions for embedded-boundary terrain.

Parameters

[in]	lev	level index
[in]	mfs	state and velocity fields used by the BC computation
[in,out]	Tau_lev	EB stress fields to fill
[in,out]	xheat_flux	x-face heat flux field
[in,out]	yheat_flux	y-face heat flux field
[in,out]	zheat_flux	z-face heat flux field
[in,out]	xqv_flux	x-face moisture flux field
[in,out]	yqv_flux	y-face moisture flux field
[in,out]	zqv_flux	z-face moisture flux field

Wrapper to impose Monin Obukhov similarity theory fluxes by populating ghost cells.

Parameters

[in]	lev	Current level
[in]	mfs	State MultiFabs used to compute the EB boundary fluxes
[in,out]	Tau_EB	EB diffusive stress MultiFabs populated with surface stresses
[in,out]	xheat_flux	x-face EB heat-flux MultiFab, currently unused
[in,out]	yheat_flux	y-face EB heat-flux MultiFab, currently unused
[in,out]	Hfx3_EB	EB heat-flux MultiFab populated with scalar surface flux
[in,out]	xqv_flux	x-face EB moisture-flux MultiFab, currently unused
[in,out]	yqv_flux	y-face EB moisture-flux MultiFab, currently unused
[in,out]	zqv_flux	z-face EB moisture-flux MultiFab, currently unused

 {
     if (flux_type == FluxCalcType::MOENG) {
         moeng_flux_eb flux_comp;
         compute_SurfaceLayer_bcs_EB(lev, mfs, Tau_EB,
                                  xheat_flux, yheat_flux, Hfx3_EB,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  flux_comp);
     } else {
         amrex::Abort("Not implemented surface layer flux calculation type for EB");
     }
 }

◆ init_tke_from_ustar()

void SurfaceLayer::init_tke_from_ustar	(	const int &	lev,
		amrex::MultiFab &	cons,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const amrex::Real	tkefac = `one`,
		const amrex::Real	zscale = `amrex::Real(700.0)`
	)

Initialize TKE from the current surface friction velocity.

Parameters

[in]	lev	level index
[in,out]	cons	conserved state whose TKE component is initialized
[in]	z_phys_nd	nodal physical-height field
[in]	tkefac	scale factor applied to the initialized TKE
[in]	zscale	vertical decay scale

Initialize TKE from surface-layer friction velocity.

Parameters

[in]	lev	Current level
[in,out]	cons	Conserved state whose RhoKE component is initialized
[in]	z_phys_nd	Nodal physical height used to compute height above ground
[in]	tkefac	Factor multiplying ustar squared for the surface TKE value
[in]	zscale	Scale factor used to taper TKE with height

 {
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         static_cast<int>(m_face) == Orientation::zlo(),
         "TKE initialization from surface-layer ustar is supported only on the z-low face.");
  
     Print() << "Initializing TKE from surface layer ustar on level " << lev << std::endl;
  
     // Handle vertical decomposition by selectively copying into
     // a FArrayBox section on each rank. Then doing a reduce real sum
     // and broadcasting to each rank. No mask since all CC data
     //
     // Pinned so the reduction below can read them on the host; the device
     // still reaches pinned memory, so the loops here are unaffected.
     const int klo = m_geom[lev].Domain().smallEnd(2);
     Box bx_lo = u_star[lev]->boxArray().minimalBox();
     FArrayBox u_star_lo(bx_lo, 1, The_Pinned_Arena()); u_star_lo.setVal<RunOn::Host>(0);
     FArrayBox z_surf_lo(bx_lo, 1, The_Pinned_Arena()); z_surf_lo.setVal<RunOn::Host>(0);
     Real* ustar_ptr = u_star_lo.dataPtr();
     Real* zsurf_ptr = z_surf_lo.dataPtr();
     for (MFIter mfi(cons); mfi.isValid(); ++mfi)
     {
         Box vbx = mfi.validbox();
         if (vbx.smallEnd(2) != klo) { continue; }
         vbx.makeSlab(2,0);
  
         auto const& u_star_arr = u_star[lev]->const_array(mfi);
         auto        u_star_all = u_star_lo.array();
  
         auto const& z_phys_arr = z_phys_nd->const_array(mfi);
         auto        z_surf_all = z_surf_lo.array();
  
         ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
         {
             u_star_all(i,j,0) = u_star_arr(i,j,0);
             z_surf_all(i,j,0) = fourth * ( z_phys_arr(i  ,j  ,klo) + z_phys_arr(i+1,j  ,klo)
                                        + z_phys_arr(i  ,j+1,klo) + z_phys_arr(i+1,j+1,klo) );
         });
     }
     Gpu::streamSynchronize(); // the fills above are async, the reduction is not
     ParallelDescriptor::ReduceRealSum(ustar_ptr, static_cast<int>(bx_lo.numPts()));
     ParallelDescriptor::ReduceRealSum(zsurf_ptr, static_cast<int>(bx_lo.numPts()));
  
     // Now work on all boxes (ustar has been filled above)
     constexpr Real small = Real(0.01);
     for (MFIter mfi(cons); mfi.isValid(); ++mfi)
     {
         Box vbx  = mfi.validbox();
  
         auto const& u_star_arr = u_star_lo.const_array();
         auto const& z_surf_arr = z_surf_lo.const_array();
         auto const& z_phys_arr = z_phys_nd->const_array(mfi);
  
         auto const& cons_arr   = cons.array(mfi);
  
         ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             Real rho  = cons_arr(i, j, k, Rho_comp);
             Real ust  = u_star_arr(i, j, 0);
             Real tke0 = tkefac * ust * ust; // surface value
             Real zagl = Compute_Z_AtCellCenter(i, j, k, z_phys_arr) - z_surf_arr(i,j,0);
  
             // linearly tapering profile --  following WRF, approximate top of
             // PBL as ustar * zscale
             cons_arr(i, j, k, RhoKE_comp) = rho * tke0 * std::max(
                 (ust * zscale - zagl) / (std::max(ust, small) * zscale),
                 small);
         });
     }
 }

Here is the call graph for this function:

◆ inject_pblh_from_coarser()

void SurfaceLayer::inject_pblh_from_coarser	(	const int &	lev,
		amrex::MultiFab &	fine
	)		const

The PBL height of the next coarser level, injected onto the grids of this level over all of their cells, ghost cells included. Public because it launches a device lambda.

Parameters

[in]	lev	Current level (> 0)
[out]	fine	MultiFab on the grids of pblh[lev]

Write into fine, a MultiFab on the grids of pblh[lev], the PBL height of the next coarser level over all of its cells, ghost cells included (see fill_pblh_from_coarser).

Parameters

[in]	lev	Current level (> 0)
[out]	fine	PBL height of level lev-1, injected onto the grids of level lev

 {
     AMREX_ALWAYS_ASSERT(lev > 0 && pblh[lev-1]);
     const Geometry& geom_f = m_geom[lev];
     const Geometry& geom_c = m_geom[lev-1];
     const IntVect rr(AMREX_D_DECL(static_cast<int>(std::lround(geom_c.CellSize(0) / geom_f.CellSize(0))),
                                   static_cast<int>(std::lround(geom_c.CellSize(1) / geom_f.CellSize(1))),
                                   1));
  
     const MultiFab& crse = *pblh[lev-1];
  
     // Each fine column takes the height of the coarse column that holds it (injection, no
     // interpolation): the fine height is then exactly the coarse diagnostic, and independent of
     // the decomposition, at the price of being constant over each block of rr x rr fine columns.
     // The coarse level's own height varies at that resolution too.  Interpolating would reach
     // further into the coarse halo than these planar fields hold: a bilinear stencil for a fine
     // ghost cell at a non-periodic domain boundary wants a second coarse cell outside the
     // domain, and pblh carries one ghost cell.  A smooth cap would need a wider halo there, or
     // a one-sided stencil.
     //
     // The coarse heights over the boxes of this level, with enough ghost cells to cover theirs.
     // Ghost cells go first, then the valid cells, so every cell the coarse level owns comes
     // from the box that owns it.  The estimator writes pblh over its grown box (it clips only
     // in z), so the coarse ghost cells the first copy reads hold a diagnosed height, outside a
     // non-periodic domain too.
     const IntVect ng_f = fine.nGrowVect();
     const IntVect ng_c(AMREX_D_DECL((ng_f[0] + rr[0] - 1) / rr[0], (ng_f[1] + rr[1] - 1) / rr[1], ng_f[2]));
     MultiFab crse_on_fine(amrex::coarsen(fine.boxArray(), rr), fine.DistributionMap(), 1, ng_c);
     crse_on_fine.setVal(zero);
     for (const IntVect& ng_src : {crse.nGrowVect(), IntVect(0)}) {
         crse_on_fine.ParallelCopy(crse, 0, 0, 1, elemwiseMin(ng_src, ng_c), ng_c, geom_c.periodicity());
     }
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
     for (MFIter mfi(fine, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
         const Box bx = mfi.growntilebox();
         const Array4<Real> f = fine.array(mfi);
         const Array4<Real const> c = crse_on_fine.const_array(mfi);
         const int rx = rr[0];
         const int ry = rr[1];
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             f(i,j,k) = c(amrex::coarsen(i,rx), amrex::coarsen(j,ry), k);
         });
     }
 }

Here is the call graph for this function:

◆ interpolate_sfc_column()

Real SurfaceLayer::interpolate_sfc_column	(	const amrex::Real &	time,
		int	col
	)		const

Interpolates the SFC/SST data at the given column and time

Parameters

[in]	time	elapsed time
[in]	col	column index of file data

 {
     if (sfc.empty() || sfc[0].empty()) { return zero; }
     if (sfc[0].size() == 1) { return sfc[col][0]; }
  
     const Real t0 = sfc[0][sfc_time_ind];
     const Real t1 = sfc[0][sfc_time_ind+1];
     const Real x0 = sfc[col][sfc_time_ind];
     const Real x1 = sfc[col][sfc_time_ind+1];
  
     if (elapsed_time < t0) {
         return x0;
     }
  
     if (t0 == t1 || elapsed_time > t1) {
         return x1;
     }
  
     const Real dt = (elapsed_time - t0) / (t1 - t0);
     return x0 + (x1 - x0) * dt;
 }

◆ lmask_min_reduce()

int SurfaceLayer::lmask_min_reduce	(	amrex::iMultiFab &	lmask,
		const int &	nghost
	)

inline

Compute the minimum land-mask value over valid and optional ghost cells.

Parameters

[in]	lmask	land-mask field
[in]	nghost	number of ghost cells included in the reduction

   {
       int lmask_min = amrex::ReduceMin(lmask, nghost, [=] AMREX_GPU_HOST_DEVICE(
                       amrex::Box const& bx, amrex::Array4<int const> const& lm_arr) -> int
                       {
                           int locmin = std::numeric_limits<int>::max();
                           const auto lo = lbound(bx);
                           const auto hi = ubound(bx);
                           for (int j = lo.y; j <= hi.y; ++j) {
                               for (int i = lo.x; i <= hi.x; ++i) {
                                   locmin = std::min(locmin, lm_arr(i, j, 0));
                               }
                           }
                           return locmin;
                       });
  
       return lmask_min;
   }

Referenced by make_SurfaceLayer_at_level().

Here is the caller graph for this function:

◆ mac_avg_is_initialized()

bool SurfaceLayer::mac_avg_is_initialized ( const int & lev ) const

inline

Return whether the time filter at this level holds meaningful history.

Parameters

[in] lev level index

1405 { return m_ma.time_avg_is_initialized(lev); }

MOSTAverage::time_avg_is_initialized

bool time_avg_is_initialized(const int &lev) const

Definition: ERF_MOSTAverage.H:278

Here is the call graph for this function:

◆ mac_avg_is_time_averaged()

bool SurfaceLayer::mac_avg_is_time_averaged ( ) const

inline

Return whether the MOST averages are filtered in time.

1393 { return m_ma.do_time_averaging(); }

MOSTAverage::do_time_averaging

bool do_time_averaging() const

Definition: ERF_MOSTAverage.H:265

Here is the call graph for this function:

◆ make_SurfaceLayer_at_level()

void SurfaceLayer::make_SurfaceLayer_at_level	(	const int &	lev,
		int	nlevs,
		const amrex::Vector< amrex::MultiFab * > &	mfv,
		std::unique_ptr< amrex::MultiFab > &	Theta_prim,
		std::unique_ptr< amrex::MultiFab > &	Qv_prim,
		std::unique_ptr< amrex::MultiFab > &	Qr_prim,
		std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		amrex::MultiFab *	Hwave,
		amrex::MultiFab *	Lwave,
		amrex::MultiFab *	eddyDiffs,
		amrex::Vector< amrex::MultiFab * >	lsm_data,
		amrex::Vector< std::string >	lsm_data_name,
		amrex::Vector< amrex::MultiFab * >	lsm_flux,
		amrex::Vector< std::string >	lsm_flux_name,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	sst_lev,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	tsk_lev,
		amrex::Vector< std::unique_ptr< amrex::iMultiFab >> &	lmask_lev
	)

inline

Allocate and initialize surface-layer data for one AMR level.

Parameters

[in]	lev	level index
[in]	nlevs	number of AMR levels
[in]	mfv	conserved and velocity MultiFabs for this level
[in]	Theta_prim	primitive potential-temperature field
[in]	Qv_prim	primitive water-vapor field
[in]	Qr_prim	primitive rain-water field
[in]	z_phys_nd	nodal physical-height field
[in]	Hwave	wave-height field
[in]	Lwave	wavelength field
[in]	eddyDiffs	eddy-diffusivity field
[in]	lsm_data	land-surface-model data fields
[in]	lsm_data_name	names for lsm_data entries
[in]	lsm_flux	land-surface-model flux fields
[in]	lsm_flux_name	names for lsm_flux entries
[in]	sst_lev	sea-surface-temperature data by time
[in]	tsk_lev	skin-temperature data by time
[in]	lmask_lev	land-mask data by time

   {
       // Update MOST Average
       m_ma.make_MOSTAverage_at_level(lev, mfv,
                                      Theta_prim, Qv_prim, Qr_prim,
                                      z_phys_nd);
  
       // Get CC vars
       amrex::MultiFab& mf = *(mfv[0]);
  
       amrex::ParmParse pp(m_pp_prefix);
  
       // Do we have a time-varying surface roughness that needs to be saved?
       if (lev == 0) {
           const int nghost = 0; // ghost cells not included
           int lmask_min = lmask_min_reduce(*lmask_lev[0].get(), nghost);
           amrex::ParallelDescriptor::ReduceIntMin(lmask_min);
  
           m_var_z0 = (lmask_min < 1) & (rough_type_sea != RoughCalcType::CONSTANT);
           if (m_var_z0) {
               AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                   m_face.coordDir() == 2 && m_face.isLow(),
                   "Variable sea roughness surface-layer fluxes are supported only on the z-low face.");
  
               std::string rough_sea_string{"charnock"};
               pp.queryAdd("most.roughness_type_sea", rough_sea_string);
               amrex::Print() << "Variable sea roughness (type " << rough_sea_string
                              << ")" << std::endl;
           }
       }
  
       if (m_eddyDiffs_lev.size() < lev+1) {
           m_Hwave_lev.resize(nlevs);
           m_Lwave_lev.resize(nlevs);
           m_eddyDiffs_lev.resize(nlevs);
  
           m_lsm_data_lev.resize(nlevs);
           m_lsm_flux_lev.resize(nlevs);
  
           m_sst_lev.resize(nlevs);
           m_tsk_lev.resize(nlevs);
           m_lmask_lev.resize(nlevs);
  
           m_coupled_sst_lev.resize(nlevs, nullptr);
           m_coupled_sst_valid_lev.resize(nlevs, nullptr);
  
           // Size the MOST params for all levels
           z_0.resize(nlevs);
           u_star.resize(nlevs);
           w_star.resize(nlevs);
           t_star.resize(nlevs);
           q_star.resize(nlevs);
           t_surf.resize(nlevs);
           q_surf.resize(nlevs);
           surface_diagnostic_source.resize(nlevs);
           olen.resize(nlevs);
           pblh.resize(nlevs);
           m_planar_bndry.resize(nlevs);
       }
  
       // Get pointers to SST,TSK and LANDMASK data
       int nt_tot_sst = sst_lev.size();
       m_sst_lev[lev].resize(nt_tot_sst);
       for (int nt(0); nt < nt_tot_sst; ++nt) {
           m_sst_lev[lev][nt] = sst_lev[nt].get();
       }
       int nt_tot_tsk = static_cast<int>(tsk_lev.size());
       m_tsk_lev[lev].resize(nt_tot_tsk);
       for (int nt(0); nt < nt_tot_tsk; ++nt) {
           m_tsk_lev[lev][nt] = tsk_lev[nt].get();
       }
       int nt_tot_lmask = static_cast<int>(lmask_lev.size());
       m_lmask_lev[lev].resize(nt_tot_lmask);
       for (int nt(0); nt < nt_tot_lmask; ++nt) {
           m_lmask_lev[lev][nt] = lmask_lev[nt].get();
       }
  
       // Get pointers to wave data
       m_Hwave_lev[lev] = Hwave;
       m_Lwave_lev[lev] = Lwave;
       m_eddyDiffs_lev[lev] = eddyDiffs;
  
       // Text-file driven surface forcing modes. The file always contains at
       // least time(day) and absolute sst(K); the value is normalized to the
       // canonical SurfaceLayer theta field at the z-low surface.
       pp.queryAdd("most.use_sfc_fluxes", m_use_sfc_fluxes);
       pp.queryAdd("most.use_sfc_sst", m_use_sfc_sst);
       if (m_use_sfc_fluxes && m_use_sfc_sst) {
           amrex::Abort("Only one of most.use_sfc_fluxes and most.use_sfc_sst may be enabled");
       }
       if (m_use_sfc_fluxes || m_use_sfc_sst) {
           if (m_use_sfc_fluxes) {
               AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!use_surface_model, "Cannot use surface flux forcing with LSM/Urban model!");
               // The file's H and LE already carry the density, and the stress
               // takes the local density, so a surface density would count it
               // twice in H and LE (see set_sfc_params_from_file_fluxes).
               // The value is tested, not the key: flux_type = custom has already
               // queryAdd-ed most.rhosurf with its default of 0 (use the first cell).
               amrex::Real rhosurf_in = zero;
               pp.query("most.rhosurf", rhosurf_in);
               AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rhosurf_in <= zero,
                   "most.rhosurf cannot be combined with most.use_sfc_fluxes");
           }
           if (m_terrain_type == TerrainType::EB) {
               amrex::Abort("Text-file surface forcing is not supported with EB terrain");
           }
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow(),
               "most.use_sfc_fluxes and most.use_sfc_sst are supported only on the z-low face.");
  
           // load sensible and latent heat fluxes from sfc to prescribe
           std::string sfc_file = "";
           pp.queryAdd("most.sfc_file", sfc_file);
           if (sfc_file.empty()) {
               amrex::Abort("most.sfc_file must be set when using text-file surface forcing");
           }
  
           // sfc contains: time(day) sst(K) H(W/m2) LE(W/m2) TAU(m2/s2)
           sfc = read_cols(sfc_file, 1);
  
           const int min_cols = m_use_sfc_fluxes ? 5 : 2;
           if (static_cast<int>(sfc.size()) < min_cols) {
               amrex::Abort("Surface forcing file does not contain the required number of columns");
           }
  
           // shift time column in days to be relative to current elapsed time
           const amrex::Real start_day = sfc[0][0];
           for (int i = 0; i < static_cast<int>(sfc[0].size()); ++i) {
               sfc[0][i] = 86400.0 * (sfc[0][i] - start_day);
           }
  
           if (m_use_sfc_sst) {
               if (theta_type == ThetaCalcType::HEAT_FLUX || surf_heating_rate != zero) {
                   amrex::Abort("most.use_sfc_sst cannot be combined with prescribed heat flux or surf_heating_rate");
               }
               if (moist_type == MoistCalcType::MOISTURE_FLUX) {
                   amrex::Abort("most.use_sfc_sst cannot be combined with prescribed moisture flux");
               }
               theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
               amrex::Print() << "Using MOST with prescribed SST from most.sfc_file '" << sfc_file << "' over sea" << std::endl;
           }
  
           if (m_use_sfc_fluxes) {
               flux_type = FluxCalcType::CUSTOM;
               amrex::Print() << "Using MOST with prescribed time-varying surface fluxes from '" << sfc_file << "'" << std::endl;
           }
       }
  
       // The constructor's check of most.abs_surf_temp sees only most.flux_type; the
       // modes known here replace the prescribed temperature too.  The surface-flux file
       // makes the fluxes custom, and the SST/TSK and sfc-file SST fills use the
       // prescribed value as their land temperature, which is a potential temperature there.
       if (surf_temp_is_absolute) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_use_sfc_fluxes,
               "most.abs_surf_temp cannot be combined with most.use_sfc_fluxes");
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               (m_sst_lev[lev].empty() || !m_sst_lev[lev][0]) && !m_use_sfc_sst,
               "most.abs_surf_temp cannot be combined with SST/TSK or sfc-file SST input");
       }
  
       // By default an LSM field named "tsurf", "theta" or "t_surf" is copied into
       // t_surf unchanged. With this option the LSM skin temperature ("tsurf" from
       // SLM, "t_sfc" from Noah-MP) is taken as an absolute temperature and
       // converted to potential temperature at the surface pressure.
       pp.queryAdd("most.use_lsm_skin_temperature", m_use_lsm_skin_temperature);
  
       // Get pointers to LSM data and Fluxes
       int ndata = static_cast<int>(lsm_data.size());
       int nflux = static_cast<int>(lsm_flux.size());
       m_lsm_data_name.resize(ndata);
       m_lsm_data_lev[lev].resize(ndata);
       m_lsm_flux_name.resize(nflux);
       m_lsm_flux_lev[lev].resize(nflux);
       for (int n(0); n < ndata; ++n) {
           m_lsm_data_name[n]     = lsm_data_name[n];
           m_lsm_data_lev[lev][n] = lsm_data[n];
           const std::string lc_name = amrex::toLower(lsm_data_name[n]);
           // With most.use_lsm_skin_temperature, "tsurf" (SLM) and "t_sfc" (Noah-MP)
           // are absolute skin temperatures and are converted to potential
           // temperature in get_lsm_tsurf; "theta" and "t_surf" are taken to be
           // potential temperature already.
           if (lc_name == "tsurf" || (m_use_lsm_skin_temperature && lc_name == "t_sfc") ||
               lc_name == "theta" || lc_name == "t_surf") {
               m_has_lsm_tsurf  = true;
               m_lsm_tsurf_indx = n;
               m_lsm_tsurf_is_theta = (lc_name == "theta" || lc_name == "t_surf");
           }
       }
       if (m_use_sfc_sst && m_use_coupled_sst) {
           amrex::Abort("most.use_sfc_sst cannot be combined with an ocean LSM t_surf input");
       }
       int n_valid_lsm_flux = 0;
       bool has_soil_t_flux = false;
       for (int n(0); n < nflux; ++n) {
           m_lsm_flux_name[n]     = lsm_flux_name[n];
           m_lsm_flux_lev[lev][n] = lsm_flux[n];
           if (m_lsm_flux_lev[lev][n]) { ++n_valid_lsm_flux; }
           if (amrex::toLower(m_lsm_flux_name[n]) == "soil_t_flux") {
               has_soil_t_flux = true;
           }
       }
       AMREX_ALWAYS_ASSERT((n_valid_lsm_flux==0 || n_valid_lsm_flux>=4 ||
                            (n_valid_lsm_flux==1 && has_soil_t_flux)));
       if (n_valid_lsm_flux>=4) { m_has_lsm_fluxes = true; }
  
       const bool use_sst = (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]);
       const bool use_tsk = (!m_tsk_lev[lev].empty() && m_tsk_lev[lev][0]);
  
       // Check if there is a user-specified roughness file to be read
       std::string fname;
       bool read_z0 = false;
       if ( (flux_type == FluxCalcType::MOENG) ||
            (flux_type == FluxCalcType::ROTATE)) {
           int count = pp.countval("most.roughness_file_name");
           if (count > 1) {
               AMREX_ALWAYS_ASSERT(count >= lev+1);
               pp.query("most.roughness_file_name", fname, lev);
               read_z0 = true;
           } else if (count == 1) {
               if (lev == 0) {
                   pp.queryAdd("most.roughness_file_name", fname);
               } else {
                   // we will interpolate from the coarsest level
                   fname = "";
               }
               read_z0 = true;
           }
           // else use z0_const
       }
       if (read_z0) {
           AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
               m_face.coordDir() == 2 && m_face.isLow(),
               "Custom MOST roughness is supported only on the z-low face.");
       }
  
       // LSM flux arrays are planar and compute_sfc_params_from_lsm_fluxes only
       // writes the k=0 slab, while EB MOST consumes surface parameters at
       // arbitrary cut-cell k; no valid planar-to-cut-cell mapping exists yet.
       if (!erf_surface_layer::planar_sources_supported_for_terrain(
               m_terrain_type, use_sst, use_tsk, m_use_coupled_sst,
               m_has_lsm_tsurf, m_has_lsm_fluxes, read_z0)) {
           amrex::Abort(
               "EB SurfaceLayer does not support planar SST/TSK, coupled SST, LSM surface "
               "temperature or fluxes, or custom/file-driven roughness; no mapping exists "
               "from those planar inputs to arbitrary EB cut cells.");
       }
  
       // Attributes for MFs and FABs
       //--------------------------------------------------------
       // Create a 2D ba for planar terrain, 3D for EB terrain
       const int dir = m_face.coordDir();
       int sm_index;
       if (m_face.isLow()) {
           sm_index = m_geom[lev].Domain().smallEnd(dir);
       } else {
           sm_index = m_geom[lev].Domain().bigEnd(dir);
       }
  
       amrex::BoxArray ba = mf.boxArray();
       amrex::BoxArray ba_flux;
       amrex::IntVect ng{1,1,0};
  
       // The lateral-wall implementation uses the two-dimensional land-mask
       // layout to identify grids on the requested face.  Lateral walls are
       // intentionally restricted to grids spanning the complete z domain;
       // z-decomposed and partial-height grids are not supported.
       if (dir != 2) {
           const int dom_lo_z = m_geom[lev].Domain().smallEnd(2);
           const int dom_hi_z = m_geom[lev].Domain().bigEnd(2);
           for (int ibox = 0; ibox < ba.size(); ++ibox) {
               AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                   ba[ibox].smallEnd(2) == dom_lo_z && ba[ibox].bigEnd(2) == dom_hi_z,
                   "Surface layer boundaries on x/y faces require Cartesian grids that "
                   "span the full level z domain; partial-height refined grids and grids "
                   "decomposed in z are not supported. Set erf.max_grid_size_z accordingly.");
           }
       }
  
       if (m_terrain_type == TerrainType::EB) {
           // Use full 3D BoxArray for EB terrain
           ba_flux = ba;
           ng = amrex::IntVect{1,1,1};  // Include z ghost cells
       } else {
           // Collapse to 2D for planar terrain
           amrex::BoxList bl2d = ba.boxList();
           for (auto& b : bl2d) {
             b.setRange(dir,sm_index);
           }
           ba_flux = amrex::BoxArray(std::move(bl2d));
           // Lateral faces need state-width ghosts in the tangential z
           // direction.  A z face only needs the original one-cell x/y halo;
           // kernels using the wider state/mask boxes clip to the FAB bounds.
           if (dir == 2) {
               ng = amrex::IntVect{1,1,0};
           } else {
               ng = mf.nGrowVect();
               ng[dir] = 0;
           }
       }
  
       const amrex::DistributionMapping& dm = mf.DistributionMap();
       const int ncomp = 1;
  
       // Surface copies of the planar boxes (see fill_planar_boundary)
       // PlanarBoundary handles duplicate boxes created by a z-split. The x/y
       // layouts use selective lateral exchange instead.
       if (m_terrain_type != TerrainType::EB && dir == 2) {
           const int ksurface = m_face.isLow()
               ? m_geom[lev].Domain().smallEnd(2)
               : m_geom[lev].Domain().bigEnd(2);
           m_planar_bndry[lev].define(ba, ba_flux, dm, ksurface, m_face.isLow());
       }
  
       // Z0 heights FAB
       //--------------------------------------------------------
       z_0[lev].define(ba_flux, dm, ncomp, ng);
       z_0[lev].setVal(z0_const);
       if (read_z0) {
           read_custom_roughness(lev, fname);
       }
  
       // 2D MFs for U*, T*, T_surf
       //--------------------------------------------------------
       u_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       u_star[lev]->setVal(bogus_large_value);
  
       w_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       w_star[lev]->setVal(bogus_large_value);
  
       t_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       t_star[lev]->setVal(zero); // default to neutral
  
       q_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       q_star[lev]->setVal(zero); // default to dry
  
       olen[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       olen[lev]->setVal(bogus_large_value);
  
       pblh[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       pblh[lev]->setVal(bogus_large_value);
  
       t_surf[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       t_surf[lev]->setVal(default_land_surf_temp);
  
       q_surf[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       q_surf[lev]->setVal(default_land_surf_moist);
  
       const amrex::iMultiFab& surface_mask = *m_lmask_lev[lev][0];
       AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
           surface_mask.boxArray().size() == mf.boxArray().size(),
           "Surface-layer mask and state must have the same number of boxes.");
       AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
           surface_mask.DistributionMap() == mf.DistributionMap(),
           "Surface-layer mask and state must have identical ownership.");
       AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
           surface_mask.boxArray().size() == u_star[lev]->boxArray().size(),
           "Surface-layer mask and parameters must have the same number of boxes.");
       AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
           surface_mask.DistributionMap() == u_star[lev]->DistributionMap(),
           "Surface-layer mask and parameters must have identical ownership.");
  
       set_immersed_surface_params(lev);
  
       surface_diagnostic_source[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
       surface_diagnostic_source[lev]->setVal(
           surface_diagnostics::to_plot_value(surface_diagnostics::SurfaceDiagnosticSource::Missing));
  
       // TODO: Do we want an enum struct for indexing?
  
       if (use_sst || use_tsk || m_has_lsm_tsurf || m_use_coupled_sst || use_surface_model) {
           // Valid SST, TSK, LSM or coupled-ocean data; t_surf set before computing
           // fluxes (avoids extended lambda capture) Note that land temp will be set
           // from m_tsk_lev while sea temp will be set from m_sst_lev
           theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
  
           // Pathways in fill_tsurf_with_sst_and_tsk
           amrex::Print() << "Using MOST with specified surface temperature ";
           if (m_has_lsm_tsurf && !use_sst && !use_tsk) {
               amrex::Print() << "(LSM: " << m_lsm_data_name[m_lsm_tsurf_indx] << ")";
           } else if (!use_sst && !use_tsk) {
               amrex::Print() << "(land: T0, sea: none)";
           } else {
               // NOTE: SST from the LOW file populates TSK in update_sst_tsk.
               //       So if we have TSK, it contains everything and has been
               //       sanity checked for valid SST values.
               if (use_tsk) { m_ignore_sst = true; }
               if (use_tsk) {
                   amrex::Print() << "(land: TSK, ";
               } else {
                   amrex::Print() << "(land: T0, ";
               }
               if (use_tsk && !use_sst) {
                   amrex::Print() << "sea: TSK)";
               } else {
                   amrex::Print() << "sea: SST)";
                   AMREX_ALWAYS_ASSERT(m_sst_lev[lev][0]);
               }
           }
           // The coupler is layered on top of whatever the above selected: it
           // overwrites only the water cells it actually covers, so the pathway
           // named above remains the value for land and for uncovered water.
           if (m_use_coupled_sst) {
               amrex::Print() << " + coupled ocean SST where covered";
           }
           amrex::Print() << std::endl;
       }
   }

Here is the call graph for this function:

◆ read_cols()

amrex::Vector< amrex::Vector< amrex::Real > > SurfaceLayer::read_cols	(	const std::string &	fname,
		const int	skip_nlines = `1`
	)

static

Reads columns of data from a text file, returning each column in a vector.

Parameters

[in]	fname	path to text file
[in]	skip_nlines	number of lines to skip before reading data (e.g, header lines)

Returns: Vector containing each column in the file as a vector

 {
     std::ifstream ifs(fname);
     if (!ifs.is_open())
     {
         amrex::Error("Error opening input file " + fname);
     }
  
     amrex::Vector<amrex::Vector<amrex::Real>> col_data;
     std::string line;
     int nlines = 0;
     int ncols = -1;
  
     const auto print_err = [](const std::string &err_fname, int lineno, int cols, int expected_cols) {
         amrex::Error("Error reading file '" + err_fname + "': expected line " +
                      std::to_string(lineno) + " to have " + std::to_string(expected_cols) +
                      " columns, but got " + std::to_string(cols));
     };
  
     while (std::getline(ifs, line))
     {
         nlines++;
         if (nlines <= skip_nlines) continue;
         if (line.empty()) continue;
  
         std::istringstream iss(line);
  
         amrex::Real tmp;
         // Get the number of columns in the file
         if (ncols == -1) {
             int j = 0;
             while (iss >> tmp) {
                 col_data.push_back(amrex::Vector<amrex::Real>());
                 j+= 1;
             }
  
             ncols = j;
             iss = std::istringstream(line);
         }
  
         int j = 0;
         while (iss >> tmp) {
             // verify each line has the same number of columns
             if (j >= ncols) {
                 print_err(fname, nlines, j+1, ncols);
             }
             col_data[j].push_back(tmp);
             j+= 1;
         }
  
         // throw error if there are fewer columns in the line than expected
         if (j != ncols) {
             print_err(fname, nlines, j, ncols);
         }
     }
  
     ifs.close();
  
     return col_data;
 }

Referenced by make_SurfaceLayer_at_level().

Here is the caller graph for this function:

◆ read_custom_roughness()

void SurfaceLayer::read_custom_roughness	(	const int &	lev,
		const std::string &	fname
	)

Read custom roughness data for one level.

Parameters

[in]	lev	level index
[in]	fname	roughness-data file name

Read or interpolate custom roughness length data.

Parameters

[in]	lev	Current level
[in]	fname	Roughness file name; an empty name interpolates from level 0

 {
     // Read the file if we have it
     if (!fname.empty()) {
         // Only the ioproc reads the file
         Gpu::HostVector<Real> m_x,m_y,m_z0;
         if (ParallelDescriptor::IOProcessor()) {
             Print()<<"Reading MOST roughness file at level " << lev << " : " << fname << std::endl;
             std::ifstream file(fname);
             if (!file.is_open()) {
                 Abort("Could not open the MOST roughness file \"" + fname + "\" given as "
                       "erf.most.roughness_file_name -- please check the file name in the inputs file");
             }
             Real value1,value2,value3;
             while(file>>value1>>value2>>value3){
                 m_x.push_back(value1);
                 m_y.push_back(value2);
                 m_z0.push_back(value3);
             }
             file.close();
  
             AMREX_ALWAYS_ASSERT(m_x.size() == m_y.size());
             AMREX_ALWAYS_ASSERT(m_x.size() == m_z0.size());
  
             if (m_x.empty()) {
                 Abort("The MOST roughness file \"" + fname + "\" holds no (x, y, z0) triples");
             }
         }
  
         // Broadcast the whole domain to every rank
         int ioproc = ParallelDescriptor::IOProcessorNumber();
         int nnode = static_cast<int>(m_x.size());
         ParallelDescriptor::Bcast(&nnode, 1, ioproc);
  
         if (!ParallelDescriptor::IOProcessor()) {
             m_x.resize(nnode);
             m_y.resize(nnode);
             m_z0.resize(nnode);
         }
         ParallelDescriptor::Bcast(m_x.data() , nnode, ioproc);
         ParallelDescriptor::Bcast(m_y.data() , nnode, ioproc);
         ParallelDescriptor::Bcast(m_z0.data(), nnode, ioproc);
  
         // The loop below indexes the file's values as nodes of this level's grid, so a file
         // holding fewer values than the grid has nodes would be read past its end
         {
             const Box& dom = m_geom[lev].Domain();
             const int nnode_needed = dom.bigEnd(0) + dom.bigEnd(1) * (dom.length(0)+1) + 1;
             if (nnode < nnode_needed) {
                 Abort("The MOST roughness file \"" + fname + "\" holds " + std::to_string(nnode) +
                       " values but the grid at level " + std::to_string(lev) + " needs at least " +
                       std::to_string(nnode_needed) + " -- is this file written for this grid?");
             }
         }
  
         // Copy data to the GPU
         Gpu::DeviceVector<Real> d_x(nnode),d_y(nnode),d_z0(nnode);
         Gpu::copy(Gpu::hostToDevice, m_x.begin(), m_x.end(), d_x.begin());
         Gpu::copy(Gpu::hostToDevice, m_y.begin(), m_y.end(), d_y.begin());
         Gpu::copy(Gpu::hostToDevice, m_z0.begin(), m_z0.end(), d_z0.begin());
         Real* xp  = d_x.data();
         Real* yp  = d_y.data();
         Real* z0p = d_z0.data();
  
         // Each rank populates it's z_0[lev] MultiFab
         const int klo = m_geom[lev].Domain().smallEnd(2);
         for (MFIter mfi(z_0[lev]); mfi.isValid(); ++mfi)
         {
             Box gtbx = mfi.growntilebox();
  
             if (gtbx.smallEnd(2) != klo) { continue; }
  
             // Populate z_phys data
             Real tol = Real(1.0e-4);
             auto dx = m_geom[lev].CellSizeArray();
             auto ProbLoArr = m_geom[lev].ProbLoArray();
             int ilo = m_geom[lev].Domain().smallEnd(0);
             int jlo = m_geom[lev].Domain().smallEnd(1);
             int ihi = m_geom[lev].Domain().bigEnd(0);
             int jhi = m_geom[lev].Domain().bigEnd(1);
  
             Array4<Real> const& z0_arr = z_0[lev].array(mfi);
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
             {
                 // Clip indices for ghost-cells
                 int ii = amrex::min(amrex::max(i,ilo),ihi);
                 int jj = amrex::min(amrex::max(j,jlo),jhi);
  
                 // Location of nodes
                 Real x = ProbLoArr[0]  + ii  * dx[0];
                 Real y = ProbLoArr[1]  + jj  * dx[1];
                 int inode = ii + jj * (ihi-ilo+2); // stride is Nx+1
                 if (std::sqrt(amrex::Math::powi<2>(x-xp[inode])+amrex::Math::powi<2>(y-yp[inode])) < tol) {
                     z0_arr(i,j,klo) = z0p[inode];
                 } else {
                     // Unexpected list order, do brute force search
                     Real z0loc = zero;
                     bool found = false;
                     for (int n=0; n<nnode; ++n) {
                         Real delta=std::sqrt(amrex::Math::powi<2>(x-xp[n])+amrex::Math::powi<2>(y-yp[n]));
                         if (delta < tol) {
                             found = true;
                             z0loc = z0p[n];
                             break;
                         }
                     }
                     AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
                     amrex::ignore_unused(found);
                     z0_arr(i,j,klo) = z0loc;
                 }
             });
         } // mfi
     } else {
         AMREX_ALWAYS_ASSERT(lev > 0);
  
         Print()<<"Interpolating MOST roughness at level " << lev << std::endl;
  
         // Create a BC mapper that uses FOEXTRAP at domain bndry
         Vector<int> bc_lo(3,ERFBCType::foextrap);
         Vector<int> bc_hi(3,ERFBCType::foextrap);
         Vector<BCRec> bcr; bcr.push_back(BCRec(bc_lo.data(),bc_hi.data()));
  
         // Create ref ratio
         IntVect ratio;
         for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
             ratio[idim] = m_geom[lev].Domain().length(idim) / m_geom[0].Domain().length(idim);
         }
  
         // Create interp object and interpolate from the coarsest grid
         MFInterpolater* interp = &mf_cell_cons_interp;
         interp->interp(z_0[0]  , 0,
                        z_0[lev], 0,
                        1, z_0[lev].nGrowVect(),
                        m_geom[0], m_geom[lev],
                        m_geom[lev].Domain(),ratio,
                        bcr, 0);
     }
 }

Referenced by make_SurfaceLayer_at_level().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ rotates_surface_fluxes()

bool SurfaceLayer::rotates_surface_fluxes ( ) const

inline

Return whether the surface fluxes are rotated with the terrain slope (erf.use_rotate_surface_flux), splitting them over the x, y and z faces.

1576 { return m_rotate; }

◆ roughness()

amrex::Real SurfaceLayer::roughness ( ) const

inline

The constant roughness length [m] the levels take.

1565 { return z0_const; }

◆ set_coupled_sst_active()

void SurfaceLayer::set_coupled_sst_active ( const bool active )

inline

Declare that an ocean coupler will supply SST for this run.

Must be called before make_SurfaceLayer_at_level, which uses it to select ThetaCalcType::SURFACE_TEMPERATURE.

Parameters

[in] active whether coupled SST is configured

1718 { m_use_coupled_sst = active; }

◆ set_default_roughness()

bool SurfaceLayer::set_default_roughness ( amrex::Real z0 )

inline

The land roughness length [m] to use when erf.most.z0 is not given (the two-stream balance's, from Noah-MP's tables). Call before the levels are made, which set z0 from it.

Taking it also replaces the erf.most.z0 entry that the constructor's queryAdd put in the ParmParse table. job_info (ParmParse::dumpTable) then records the roughness the run used.

Returns: whether the roughness was taken: an erf.most.z0 given in the inputs stands.

   {
       if (m_z0_given) { return false; }
       z0_const = z0;
       amrex::ParmParse pp(m_pp_prefix);
       pp.remove("most.z0");
       pp.add("most.z0", z0_const);
       return true;
   }

Here is the call graph for this function:

◆ set_immersed_surface_params()

void SurfaceLayer::set_immersed_surface_params ( const int & lev )

Set benign surface scales in immersed columns (land mask 2), including ghosts. Called after allocation, checkpoint restoration, and surface-scale updates. Preserve PBL height, roughness, surface temperature, and surface moisture.

Parameters

[in] lev level index

 {
     // The land mask describes ground columns, not the full height of a lateral wall.
     if (m_face.coordDir() != 2 ||
         (m_terrain_type != TerrainType::ImmersedForcing && !m_has_immersed_buildings)) {
         return;
     }
     const auto& surface_mask = *m_lmask_lev[lev][0];
     for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi) {
         Box bx = mfi.fabbox();
         // Limit the halo to the mask's available ghost cells.
         Box mask_columns = surface_mask.fabbox(mfi.index());
         mask_columns.setRange(2, bx.smallEnd(2), bx.length(2));
         bx &= mask_columns;
  
         const auto mask_arr = surface_mask.const_array(mfi);
         const auto u_star_arr = u_star[lev]->array(mfi);
         const auto w_star_arr = w_star[lev]->array(mfi);
         const auto t_star_arr = t_star[lev]->array(mfi);
         const auto q_star_arr = q_star[lev]->array(mfi);
         const auto olen_arr = olen[lev]->array(mfi);
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
             if (mask_arr(i,j,0) == 2) {
                 // MOST and LSM updates skip immersed columns. Finite, inactive
                 // scales keep their initialization sentinels out of stress/TKE
                 // calculations; immersed forcing handles roof and wall physics.
                 u_star_arr(i,j,k) = zero;
                 w_star_arr(i,j,k) = zero;
                 t_star_arr(i,j,k) = zero;
                 q_star_arr(i,j,k) = zero;
                 olen_arr(i,j,k) = immersed_neutral_obukhov_length;
             }
         });
     }
 }

Referenced by make_SurfaceLayer_at_level().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ set_mac_avg_initialized()

void SurfaceLayer::set_mac_avg_initialized ( const int & lev )

inline

Declare the time filter at this level to hold meaningful history (restart).

Parameters

[in] lev level index

1412 { m_ma.set_time_avg_initialized(lev); }

MOSTAverage::set_time_avg_initialized

void set_time_avg_initialized(const int &lev)

Definition: ERF_MOSTAverage.H:289

Here is the call graph for this function:

◆ set_mac_plane_avg()

bool SurfaceLayer::set_mac_plane_avg	(	const int &	lev,
		const amrex::Vector< amrex::Real > &	pavg
	)

inline

Restore the filtered plane averages from a checkpoint; returns false if the checkpoint does not hold what this run expects.

Parameters

[in]	lev	level index
[in]	pavg	filtered plane averages, one per average component

1437 { return m_ma.set_plane_average(lev, pavg); }

MOSTAverage::set_plane_average

bool set_plane_average(const int &lev, const amrex::Vector< amrex::Real > &pavg)

Definition: ERF_MOSTAverage.H:315

Here is the call graph for this function:

◆ set_pblh()

void SurfaceLayer::set_pblh	(	const int &	lev,
		const amrex::MultiFab &	pblh_in
	)

Store the PBL height a PBL scheme (MRF, YSU, YSUNew) has diagnosed.

pblh_in has no ghost cells, so the ghost cells of pblh are filled here: those that another grid of this level or a periodic image holds by fill_planar_boundary; those at a coarse/fine boundary from the next coarser level (the same injection as fill_pblh_from_coarser), so that they agree with the coarse state the other surface fields there are computed from; those outside a non-periodic domain from the nearest column of the domain (after fill_planar_boundary, so that a corner beside another grid takes that grid's height); and the rest – where the coarser level has no height yet – from the nearest column of their own grid.

Parameters

[in]	lev	Current level
[in]	pblh_in	PBL height, on the grids of the eddy diffusivities (3D, no ghost cells)

 {
     // NOTE: pblh[lev] is a 2D MF (except with EB)
     //       pblh_in   is a 3D MF built off the eddydiffs BA
     AMREX_ALWAYS_ASSERT(pblh[lev]);
     MultiFab& dst = *pblh[lev];
  
     // The boxes of the two are paired by index below
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         pblh_in.size() == dst.size() && pblh_in.DistributionMap() == dst.DistributionMap(),
         "SurfaceLayer::set_pblh: the PBL height must be on the grids of the surface layer");
  
     std::unique_ptr<MultiFab> crse;
     if (lev > 0 && pblh[lev-1]) {
         crse = std::make_unique<MultiFab>(dst.boxArray(), dst.DistributionMap(), 1, dst.nGrowVect());
         inject_pblh_from_coarser(lev, *crse);
     }
     // Only the columns inside the domain (or a periodic image of it) take the coarse height
     const Geometry& geom = m_geom[lev];
     const NearestColumn in_domain(geom.Domain());
     const bool per_x = geom.isPeriodic(0);
     const bool per_y = geom.isPeriodic(1);
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
     for (MFIter mfi(dst, TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
         const Box gtbx = mfi.growntilebox();
         const Box& vbx = mfi.validbox();
         const Box& sbx = pblh_in.box(mfi.index());
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
             sbx.smallEnd(0) == vbx.smallEnd(0) && sbx.bigEnd(0) == vbx.bigEnd(0) &&
             sbx.smallEnd(1) == vbx.smallEnd(1) && sbx.bigEnd(1) == vbx.bigEnd(1),
             "SurfaceLayer::set_pblh: the PBL height must be on the grids of the surface layer");
  
         // Clamp to the source box, which is the one we read
         const NearestColumn col(sbx);
         const int klo = sbx.smallEnd(2);
         const int khi = sbx.bigEnd(2);
  
         const Array4<      Real> pblh_dst = dst.array(mfi);
         const Array4<const Real> pblh_src = pblh_in.const_array(mfi);
         const Array4<const Real> pblh_crse = (crse) ? crse->const_array(mfi) : Array4<const Real>{};
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             const int li = col.i(i);
             const int lj = col.j(j);
             const int lk = amrex::min(amrex::max(k, klo), khi);
             Real h = pblh_src(li,lj,lk);
             const bool inside = (per_x || in_domain.i(i) == i) && (per_y || in_domain.j(j) == j);
             if (pblh_crse && inside && (li != i || lj != j) && pblh_is_diagnosed(pblh_crse(i,j,k))) {
                 h = pblh_crse(i,j,k);
             }
             pblh_dst(i,j,k) = h;
         });
     }
     fill_planar_boundary(lev, dst);
  
     // A ghost cell outside a non-periodic domain but beside another grid of this level (a
     // corner) was left above with the nearest column of its own grid, which fill_planar_boundary
     // does not overwrite; give it the nearest column of the domain, which is inside this fab and
     // now holds the height of the grid that owns it
     if (!per_x || !per_y) {
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(dst, TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
             const Box gtbx = mfi.growntilebox();
             const Array4<Real> pblh_dst = dst.array(mfi);
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 const int ci = (per_x) ? i : in_domain.i(i);
                 const int cj = (per_y) ? j : in_domain.j(j);
                 if (ci != i || cj != j) {
                     pblh_dst(i,j,k) = pblh_dst(ci,cj,k);
                 }
             });
         }
     }
 }

Here is the call graph for this function:

◆ set_q_surf()

void SurfaceLayer::set_q_surf	(	const int &	lev,
		const amrex::Real	qsurf
	)

inline

Set the surface-moisture field to a constant value.

Parameters

[in]	lev	level index
[in]	qsurf	surface moisture

1599 { q_surf[lev]->setVal(qsurf); }

◆ set_sfc_params_from_file_fluxes()

void SurfaceLayer::set_sfc_params_from_file_fluxes	(	const int &	lev,
		const amrex::MultiFab &	cons_in,
		const amrex::Real &	tflux_in,
		const amrex::Real &	qflux_in,
		const amrex::Real &	sfc_tau
	)

Set kinematic u*, t*, and q* from the fluxes interpolated from most.sfc_file.

Parameters

[in]	lev	level index
[in]	cons_in	conserved state supplying the first-cell density
[in]	tflux_in	sensible heat flux H [W/m2]
[in]	qflux_in	latent heat flux LE [W/m2]
[in]	sfc_tau	kinematic surface stress TAU [m2/s2]

Set u*, t*, and q* from the fluxes interpolated from most.sfc_file.

The file's H and LE [W/m2] are dynamic fluxes, but its TAU [m2/s2] is kinematic. All three are stored kinematic, as most.ustar/tstar/qstar are without most.rhosurf, and custom_flux multiplies each by the local density: H/Cp_d and LE/L_v are divided here by the first-cell density that custom_flux multiplies back, so they are applied unchanged, and the applied stress is rho*TAU.

Parameters

[in]	lev	Current level
[in]	cons_in	Conserved state supplying the first-cell density
[in]	tflux_in	Sensible heat flux H [W/m2]
[in]	qflux_in	Latent heat flux LE [W/m2]
[in]	sfc_tau	Kinematic surface stress TAU [m2/s2]

 {
     u_star[lev]->setVal(std::sqrt(sfc_tau));
  
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const Real rho_tflux = tflux_in / Cp_d;
     const Real rho_qflux = qflux_in / L_v;
  
     for (MFIter mfi(*t_star[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         // Only the copy of a z-split planar field on the z-low face reads the 3-D
         // state; fill_planar_boundary at the end of update_fluxes fills the others.
         if (gtbx.smallEnd(2) != klo ||
             !m_planar_bndry[lev].is_surface_copy(mfi.index())) {
             continue;
         }
         gtbx &= cons_in.fabbox(mfi.index());
         if (gtbx.isEmpty()) { continue; }
  
         const auto cons_arr   = cons_in.const_array(mfi);
         const auto t_star_arr = t_star[lev]->array(mfi);
         const auto q_star_arr = q_star[lev]->array(mfi);
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             const Real rho = cons_arr(i,j,klo,Rho_comp);
             t_star_arr(i,j,k) = rho_tflux / rho;
             q_star_arr(i,j,k) = rho_qflux / rho;
         });
     }
 }

Here is the call graph for this function:

◆ set_skin_moisture()

void SurfaceLayer::set_skin_moisture	(	const int &	lev,
		const amrex::MultiFab *	moisture,
		bool	use_resistances = `false`,
		amrex::Real	hs = `amrex::Real(0)`
	)

inline

Take the land surface mixing ratio of level lev from the skin set by set_skin_temperature (fill_qsurf_with_skin_moisture), or stop doing so (nullptr). The owner keeps the field alive and resets the pointer whenever it reallocates it.

Parameters

[in]	lev	level index
[in]	moisture	per-column field in the components of erf_surface_moisture::Component, on the surface-layer layout
[in]	use_resistances	form q_surf from the soil resistance and pore humidity and, where the field's VegetationFraction is positive, the canopy's (two_source_surface_mixing_ratio), not from the soil-water factor alone: the two-stream balance passes true whenever it has a soil type
[in]	hs	the canopy's vapour-deficit parameter [1/(kg/kg)]; 0 without vegetation

   {
       if (lev >= static_cast<int>(m_skin_moisture_lev.size())) {
           m_skin_moisture_lev.resize(lev + 1, nullptr);
           m_skin_moisture_resistances_lev.resize(lev + 1, 0);
           m_skin_moisture_hs_lev.resize(lev + 1, amrex::Real(0));
       }
       m_skin_moisture_lev[lev] = moisture;
       m_skin_moisture_resistances_lev[lev] = use_resistances ? 1 : 0;
       m_skin_moisture_hs_lev[lev] = hs;
   }

◆ set_skin_temperature()

void SurfaceLayer::set_skin_temperature	(	const int &	lev,
		const amrex::MultiFab *	t_skin
	)

inline

Take the land surface temperature of level lev from an external absolute skin temperature [K] on the surface-layer layout, or stop doing so (nullptr). update_fluxes reads it; the owner keeps it alive and resets the pointer whenever it reallocates the field.

Parameters

[in]	lev	level index
[in]	t_skin	absolute skin temperature, or nullptr

   {
       if (lev >= static_cast<int>(m_skin_tsurf_lev.size())) {
           m_skin_tsurf_lev.resize(lev + 1, nullptr);
       }
       m_skin_tsurf_lev[lev] = t_skin;
   }

◆ set_surface_layer_faces()

void SurfaceLayer::set_surface_layer_faces ( const amrex::GpuArray< int, AMREX_SPACEDIM *2 > & active_faces )

inline

Set the domain faces that use the surface-layer boundary condition.

The face set is used to suppress transpose stress writes at shared edges; the face-normal component remains owned by its corresponding face path.

Parameters

[in] active_faces whether a specific x/y/z lo/hi face is enabled

   {
       m_surface_layer_faces = active_faces;
   }

◆ set_t_surf()

void SurfaceLayer::set_t_surf	(	const int &	lev,
		const amrex::Real	tsurf
	)

inline

Set the surface-temperature field to a constant value.

Parameters

[in]	lev	level index
[in]	tsurf	surface temperature

1584 { t_surf[lev]->setVal(tsurf); }

LsmVar_SLM::tsurf

@ tsurf

Definition: ERF_SLM.H:27

◆ skin_moisture_conflict()

std::string SurfaceLayer::skin_moisture_conflict ( ) const

inline

Why this surface layer cannot take its surface mixing ratio from the skin (set_skin_moisture), or an empty string when it can: everything skin_temperature_conflict asks, plus a moisture model and a moisture flux computed from a surface mixing ratio (erf.most.surf_moist given).

   {
       const std::string temperature = skin_temperature_conflict();
       if (!temperature.empty()) { return temperature; }
       if (!use_moisture) {
           return "the run has no moisture model";
       }
       if (moist_type != MoistCalcType::SURFACE_MOISTURE) {
           return "the surface layer does not compute its moisture flux from a surface "
                  "mixing ratio (set erf.most.surf_moist, not a flux)";
       }
       return "";
   }

Here is the call graph for this function:

◆ skin_temperature_conflict()

std::string SurfaceLayer::skin_temperature_conflict ( ) const

inline

Why this surface layer cannot take its surface temperature from an external skin temperature, or an empty string when it can: the skin sets t_surf, so the surface layer must compute its heat flux from t_surf (a surface temperature given, no heating rate) on a non-EB lower surface without rotated surface fluxes, and no land model may own t_surf.

   {
       if (m_terrain_type == TerrainType::EB) {
           return "the surface layer is on EB terrain";
       }
       // Rotated, the flux is split over hfx1/hfx2/hfx3 and the balance, which removes
       // the vertical-face flux hfx3 only, would lose cos(slope) of what the air gains.
       if (m_rotate) {
           return "erf.use_rotate_surface_flux is set";
       }
       if (theta_type != ThetaCalcType::SURFACE_TEMPERATURE) {
           return "the surface layer does not compute its heat flux from a surface "
                  "temperature (set erf.most.pot_surf_temp or erf.most.abs_surf_temp, not a "
                  "flux or custom/bulk fluxes)";
       }
       if (surf_heating_rate != amrex::Real(0)) {
           return "erf.most.surf_heating_rate is set";
       }
       if (m_has_lsm_tsurf || use_surface_model) {
           return "a land-surface or surface model supplies the surface temperature";
       }
       return "";
   }

Referenced by skin_moisture_conflict().

Here is the caller graph for this function:

◆ surface_sum()

Real SurfaceLayer::surface_sum	(	const int &	lev,
		const amrex::MultiFab &	mf,
		int	comp = `0`
	)		const

Sum a planar field over the valid cells of the surface, counting each surface cell once.

On a level whose grids are split in the surface-normal direction the planar BoxArray holds one duplicate box per stacked 3D box (see PlanarBoundary), so a plain sum over the planar MultiFab counts every surface cell once per stacked box. Only the computed surface copies are read here, so the duplicates need not have been filled by fill_planar_boundary, and the gather needs no communication. On EB terrain the fields are 3D without duplicates, and the lowest plane of the domain is summed. Only for a surface layer on a z face.

Parameters

[in]	lev	level index
[in]	mf	cell-centered planar MultiFab on this surface's planar BoxArray
[in]	comp	component to sum

 {
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_face.coordDir() == 2,
         "SurfaceLayer::surface_sum is only defined for a surface layer on a z face");
     if (m_terrain_type == TerrainType::EB) {
         // EB fields live on the full 3D BoxArray without duplicates: sum the lowest plane
         return sumToLine(mf, comp, 1, m_geom[lev].Domain(), 2)[0];
     }
     // A face-centered planar field (the MOST velocity averages) shares the planar BoxArray
     // but not its index type, and neighbouring surface boxes would then share a face, which
     // this sum would count twice
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mf.boxArray().ixType().cellCentered(),
         "SurfaceLayer::surface_sum is only defined for a cell-centered planar MultiFab");
     const PlanarBoundary& pb = m_planar_bndry[lev];
     MultiFab surface_copy(pb.surface_boxes(), pb.surface_dm(), 1, 0);
     pb.gather_surface(mf, surface_copy, comp, 0, 1);
     return surface_copy.sum(0);
 }

Here is the call graph for this function:

◆ update_coupled_sst_ptr()

void SurfaceLayer::update_coupled_sst_ptr	(	const int	lev,
		amrex::MultiFab *	sst_ptr,
		amrex::iMultiFab *	valid_ptr
	)

inline

Hand over the coupled sea-surface temperature and its per-cell coverage flag.

Both pointers stay owned by the caller and must remain valid, and on the same layout, until replaced. Passing a null sst_ptr retracts the coupled lane, so the lower-boundary SST/TSK value stands everywhere again.

Parameters

[in]	lev	level index
[in]	sst_ptr	coupled sea-surface temperature [K]
[in]	valid_ptr	per-cell coverage flag; nonzero where the coupler supplied a value. Null means no cell is covered.

                                                           {
       m_coupled_sst_lev[lev]   = sst_ptr;
       m_coupled_sst_valid_lev[lev] = valid_ptr;
   }

◆ update_fluxes()

void SurfaceLayer::update_fluxes	(	const int &	lev,
		const double &	elapsed_time,
		const double &	elapsed_time_since_start_low,
		amrex::MultiFab &	cons_in,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	walldist,
		int	max_iters = `100`
	)

Update surface fluxes and related surface-layer state.

Parameters

[in]	lev	level index
[in]	elapsed_time	current elapsed simulation time
[in]	elapsed_time_since_start_low	elapsed time relative to low-data start
[in,out]	cons_in	conserved state used by the flux update
[in]	z_phys_nd	nodal physical-height field
[in]	walldist	wall-distance field
[in]	max_iters	maximum MOST iteration count

Wrapper to update ustar and tstar for Monin Obukhov similarity theory.

Parameters

[in]	lev	Current level
[in]	elapsed_time	Current simulation time
[in]	elapsed_time_since_start_low	Time since the start of the lower-boundary data
[in,out]	cons_in	Conserved state, updated when RANS TKE is initialized from surface-layer data
[in]	z_phys_nd	Nodal physical height used by terrain-aware surface calculations
[in]	walldist	Wall distance used when updating boundary TKE
[in]	max_iters	Maximum iterations to use in the MOST flux solve

 {
     bool zlo = (int) m_face == Orientation::zlo();
     // Update with SST/TSK data if we have a valid pointer.
     //
     // This runs even when an ocean coupler is active: it is the only writer of
     // t_surf over land, and it is the fallback for the water cells the coupler
     // does not cover. Coupled SST is applied below and only where the coupler
     // actually supplied a value, so the lower-boundary data is the base layer
     // rather than an alternative to it.
     if (zlo && !m_sst_lev[lev].empty() && m_sst_lev[lev][0]) {
         fill_tsurf_with_sst_and_tsk(lev, elapsed_time_since_start_low);
     }
     if (zlo && m_use_sfc_sst) {
         // Set tsurf to time varying SST from sfc file
         fill_tsurf_with_sfc_sst(lev, elapsed_time, cons_in, z_phys_nd);
     }
  
     // Apply heating rate if needed.  The rate is applied from the start of the
     // simulation, so it is given the elapsed simulation time, not the time since the
     // start of the lower-boundary data: that includes start_time, which with
     // erf.start_datetime and no wrflowinp file is the epoch time of the start.
     if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
         update_surf_temp(elapsed_time);
         if (surf_temp_is_absolute) {
             fill_tsurf_with_abs_surf_temp(lev, elapsed_time, cons_in, z_phys_nd);
         }
     }
  
     // Overwrite the covered water cells with coupled ocean SST. This must come
     // after update_surf_temp, which is a whole-domain setVal, and before
     // fill_qsurf_with_qsat, which derives sea-surface humidity from t_surf.
     if (zlo) {
         fill_tsurf_with_coupled_sst(lev, cons_in, z_phys_nd);
     }
  
     // Update land surface temp if we have a valid pointer
     if (zlo && m_has_lsm_tsurf) {
         get_lsm_tsurf(lev, cons_in, z_phys_nd);
     }
  
     // An external skin temperature (the two-stream surface energy balance's) owns
     // the land surface temperature; last, so it replaces what was written above.
     if (zlo && lev < static_cast<int>(m_skin_tsurf_lev.size()) && m_skin_tsurf_lev[lev]) {
         fill_tsurf_with_skin_temperature(lev, cons_in, z_phys_nd);
     }
  
     // Update qsurf with qsat over sea
     if (use_moisture) {
         fill_qsurf_with_qsat(lev, cons_in, z_phys_nd);
     }
  
     // Fill interior ghost cells
     fill_planar_boundary(lev, *t_surf[lev]);
  
     // Compute plane averages for all vars (regardless of flux type)
     m_ma.compute_averages(lev);
  
     // The skin's surface mixing ratio over land (the two-stream balance's), relative to the
     // reference-height mixing ratio just averaged, which the moisture flux then uses.
     if (zlo && use_moisture &&
         lev < static_cast<int>(m_skin_tsurf_lev.size()) && m_skin_tsurf_lev[lev] &&
         lev < static_cast<int>(m_skin_moisture_lev.size()) && m_skin_moisture_lev[lev]) {
         fill_qsurf_with_skin_moisture(lev, cons_in, z_phys_nd);
         fill_planar_boundary(lev, *q_surf[lev]);
     }
  
     // NOTE: Do iterations to seed variables on the first step (LSM called post step)
     //       as well as compute values where invalid LSM fluxes may reside
     //*******************************************************************************
     // ***************************************************************
     // Iterate the fluxes if moeng type
     // First iterate over land -- the only model for surface roughness
     // over land is RoughCalcType::CONSTANT
     // ***************************************************************
     if (flux_type == FluxCalcType::MOENG ||
         flux_type == FluxCalcType::ROTATE) {
         bool is_land = true;
         // Do we have a constant flux for moisture over land?
         bool cons_qflux = ( (moist_type == MoistCalcType::MOISTURE_FLUX) ||
                             (moist_type == MoistCalcType::ADIABATIC) );
         if (m_terrain_type != TerrainType::EB) {
             if (theta_type == ThetaCalcType::HEAT_FLUX) {
                 if (rough_type_land == RoughCalcType::CONSTANT) {
                     surface_flux most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_land");
                 }
             } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
                 if (rough_type_land == RoughCalcType::CONSTANT) {
                     surface_temp most_flux(surf_temp_flux, surf_moist_flux, cons_qflux,
                                            m_face.coordDir(), m_face.isLow());
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_land");
                 }
             } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
                        (moist_type == MoistCalcType::ADIABATIC)) {
                 if (rough_type_land == RoughCalcType::CONSTANT) {
                     adiabatic most_flux(surf_temp_flux, surf_moist_flux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_land");
                 }
             } else {
                 amrex::Abort("Unknown value for theta_type");
             }
         // EB
         } else {
             if (theta_type == ThetaCalcType::HEAT_FLUX) {
                 if (rough_type_land == RoughCalcType::CONSTANT) {
                     surface_flux_eb most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_land");
                 }
             } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
                 if (rough_type_land == RoughCalcType::CONSTANT) {
                     surface_temp_eb most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_land");
                 }
             } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
                        (moist_type == MoistCalcType::ADIABATIC)) {
                 if (rough_type_land == RoughCalcType::CONSTANT) {
                     adiabatic_eb most_flux(surf_temp_flux, surf_moist_flux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_land");
                 }
             } else {
                 amrex::Abort("Unknown value for theta_type");
             }
         } // EB
     } // MOENG -- LAND
  
     // Update u*/T*/q*/L over land (iterations or from LSM fluxes)
     if (m_has_lsm_fluxes && elapsed_time > zero) {
         compute_sfc_params_from_lsm_fluxes(lev, cons_in);
     }
  
     // ***************************************************************
     // Iterate the fluxes if moeng type
     // Next iterate over sea -- the models for surface roughness
     // over sea are CHARNOCK, DONELAN, MODIFIED_CHARNOCK or WAVE_COUPLED
     // NOTE: Sea surface fluxes are not supported for EB terrain
     // ***************************************************************
     if ((flux_type == FluxCalcType::MOENG ||
          flux_type == FluxCalcType::ROTATE) &&
         m_terrain_type != TerrainType::EB) {
         bool is_land = false;
         // NOTE: Do not allow default to adiabatic over sea (we have Qvs at surface)
         // Do we have a constant flux for moisture over sea?
         bool cons_qflux = (moist_type == MoistCalcType::MOISTURE_FLUX);
             if (theta_type == ThetaCalcType::HEAT_FLUX) {
                 if (rough_type_sea == RoughCalcType::CHARNOCK) {
                     surface_flux_charnock most_flux(surf_temp_flux, surf_moist_flux,
                                                     cnk_a, smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
                     surface_flux_mod_charnock most_flux(surf_temp_flux, surf_moist_flux,
                                                         depth, smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::DONELAN) {
                     surface_flux_donelan most_flux(surf_temp_flux, surf_moist_flux,
                                                    smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
                     surface_flux_wave_coupled most_flux(surf_temp_flux, surf_moist_flux,
                                                         smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_sea");
                 }
  
             } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
                 if (rough_type_sea == RoughCalcType::CHARNOCK) {
                     surface_temp_charnock most_flux(surf_temp_flux, surf_moist_flux,
                                                     cnk_a, smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
                     surface_temp_mod_charnock most_flux(surf_temp_flux, surf_moist_flux,
                                                         depth, smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::DONELAN) {
                     surface_temp_donelan most_flux(surf_temp_flux, surf_moist_flux,
                                                    smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
                     surface_temp_wave_coupled most_flux(surf_temp_flux, surf_moist_flux,
                                                         smooth_flow_visc, cons_qflux);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_sea");
                 }
  
             } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
                     (moist_type == MoistCalcType::ADIABATIC)) {
                 if (rough_type_sea == RoughCalcType::CHARNOCK) {
                     adiabatic_charnock most_flux(surf_temp_flux, surf_moist_flux,
                                                  cnk_a, smooth_flow_visc);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
                     adiabatic_mod_charnock most_flux(surf_temp_flux, surf_moist_flux,
                                                      depth, smooth_flow_visc);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::DONELAN) {
                     adiabatic_donelan most_flux(surf_temp_flux, surf_moist_flux, smooth_flow_visc);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
                     adiabatic_wave_coupled most_flux(surf_temp_flux, surf_moist_flux, smooth_flow_visc);
                     compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
                 } else {
                     amrex::Abort("Unknown value for rough_type_sea");
                 }
             } else {
                 amrex::Abort("Unknown value for theta_type");
             }
     } // MOENG -- SEA
  
     if (flux_type == FluxCalcType::CUSTOM || flux_type == FluxCalcType::RICO) {
         if (m_use_sfc_fluxes) {
             // update custom surface fluxes interpolated from file
             update_sfc_time_index(elapsed_time);
             sfc_tflux = interpolate_sfc_column(elapsed_time, 2);
             sfc_qflux = interpolate_sfc_column(elapsed_time, 3);
             sfc_ustar = interpolate_sfc_column(elapsed_time, 4);
  
             amrex::Print() << " ABLMOST: Interpolating SHF and LHF at time "
                         << elapsed_time
                         << ": SHF = " << sfc_tflux
                         << " (W/m^2) LHF = " << sfc_qflux
                         << " (W/m^2) TAU = " << sfc_ustar
                         << " (m^2/s^2)" << std::endl;
  
             set_sfc_params_from_file_fluxes(lev, cons_in, sfc_tflux, sfc_qflux, sfc_ustar);
         } else if (custom_rhosurf > 0) {
             specified_rho_surf = true;
             u_star[lev]->setVal(std::sqrt(custom_rhosurf) * custom_ustar);
             t_star[lev]->setVal(custom_rhosurf * custom_tstar);
             q_star[lev]->setVal(custom_rhosurf * custom_qstar);
         } else {
             u_star[lev]->setVal(custom_ustar);
             t_star[lev]->setVal(custom_tstar);
             q_star[lev]->setVal(custom_qstar);
         }
     }
  
     // Custom fluxes can write whole fields; reset immersed columns before
     // the RANS boundary TKE calculation or any subsequent stress calculation.
     set_immersed_surface_params(lev);
  
     if (m_update_k_rans) {
         // Clamped divisor: the select is if-converted, so 1/theta_ref runs even
         //   when theta_ref = 0 and would trip fpe_trap_zero (see ERF_SetupDiff.H)
         const bool use_ref_theta = (theta_ref > 0);
         const Real inv_theta_ref = one / amrex::max(theta_ref, std::numeric_limits<Real>::min());
         const Real l_inv_theta0  = (use_ref_theta) ? inv_theta_ref : one;
         const Real l_inv_Cmu2    = inv_Cmu2;
         const int klo = m_geom[lev].Domain().smallEnd(2);
         IntVect ng = u_star[lev]->nGrowVect(); ng[2] = 0;
  
         for (MFIter mfi(cons_in); mfi.isValid(); ++mfi)
         {
             Box gpbx = mfi.tilebox(IntVect(0),ng);
  
             if (gpbx.smallEnd(2) != klo) { continue; }
  
             gpbx.makeSlab(2,klo);
             gpbx &= cons_in.fabbox(mfi.index());
             gpbx &= walldist->fabbox(mfi.index());
             if (gpbx.isEmpty()) { continue; }
  
             auto cons_arr = cons_in.array(mfi);
             const auto& u_star_arr = u_star[lev]->const_array(mfi);
             const auto& t_star_arr = t_star[lev]->const_array(mfi);
             const auto& dist_arr   = walldist->const_array(mfi);
  
             ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 Real rho = cons_arr(i,j,k,Rho_comp);
                 if (t_star_arr(i,j,0) < -1e-8) {
                     // Only destabilizing buoyancy flux affects the boundary k
                     // tstar < 0 ==> B > 0
                     Real B = -CONST_GRAV * l_inv_theta0 * u_star_arr(i,j,0) * t_star_arr(i,j,0);
                     if (!use_ref_theta) {
                         B *= cons_arr(i,j,k,Rho_comp) /
                              cons_arr(i,j,k,RhoTheta_comp);
                     }
  
                     // Axell & Liungman 2001, Eqn. 16
                     cons_arr(i,j,k,RhoKE_comp) = rho * l_inv_Cmu2 *
                         std::pow(
                             u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0)
                             + KAPPA * B * dist_arr(i,j,k),
                         two/three);
                 } else {
                     cons_arr(i,j,k,RhoKE_comp) = rho * l_inv_Cmu2 * u_star_arr(i,j,0) * u_star_arr(i,j,0);
                 }
             });
         }
     }
  
     // If using Land and/or Urban models, then overwrite u_star and t_star with values calculated from the surface models
     //
     // Gate on whether the surface models have actually advanced, not on the lower-boundary-data
     // clock: elapsed_time_since_start_low is offset by (start_time - start_low_time), so it is
     // already positive at step 0 when the wrflowinp series starts before erf.start_time (which
     // would let the still-uninitialized u*/t*/q* overwrite the MOST values), and it never becomes
     // positive at all when the series starts after it.
     if (use_surface_model && !surf_model_fluxes && m_surf_model->fields_are_valid()) {
         for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
         {
             Box gtbx = mfi.growntilebox();
             // The surface-model fields need not carry lateral ghost cells, so a ghost
             // column takes the nearest valid column's values, as in
             // compute_sfc_params_from_lsm_fluxes; fill_planar_boundary below then fills
             // the interior and periodic ghosts.
             const NearestColumn col(mfi.validbox());
  
             auto u_star_arr = u_star[lev]->array(mfi);
             auto t_star_arr = t_star[lev]->array(mfi);
             auto q_star_arr = q_star[lev]->array(mfi);
             auto olen_arr   = olen[lev]->array(mfi);
  
             // Land mask array if it exists
             auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                      Array4<int> {};
  
             auto lsm_tstar_arr = Array4<Real> {};
             auto lsm_qstar_arr = Array4<Real> {};
             auto lsm_ustar_arr = Array4<Real> {};
             auto lsm_olen_arr  = Array4<Real> {};
             if (use_surface_model) {
                 auto mf = m_surf_model->get_field("tstar", lev);
                 lsm_tstar_arr = (mf) ? mf->array(mfi) : Array4<Real> {};
                 mf = m_surf_model->get_field("qstar", lev);
                 lsm_qstar_arr = (mf) ? mf->array(mfi) : Array4<Real> {};
                 mf = m_surf_model->get_field("ustar", lev);
                 lsm_ustar_arr = (mf) ? mf->array(mfi) : Array4<Real> {};
                 mf = m_surf_model->get_field("olen", lev);
                 lsm_olen_arr = (mf) ? mf->array(mfi) : Array4<Real> {};
             }
  
             ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 // Overwrite ustar, tstar, qstar, and olen computed above with LSM values
                 if (lmask_arr(i,j,0) == 1) {
                     const int li = col.i(i);
                     const int lj = col.j(j);
                     if (lsm_tstar_arr) t_star_arr(i,j,k) = lsm_tstar_arr(li,lj,0);
                     if (lsm_qstar_arr) q_star_arr(i,j,k) = lsm_qstar_arr(li,lj,0);
                     if (lsm_ustar_arr) u_star_arr(i,j,k) = lsm_ustar_arr(li,lj,0);
                     if (lsm_olen_arr)  olen_arr(i,j,k)   = lsm_olen_arr(li,lj,0);
                 }
             });
         }
  
     }
  
     if (m_terrain_type == TerrainType::EB || m_face.coordDir() == 2) {
         fill_planar_boundary(lev, *u_star[lev]);
         if (m_include_wstar) { fill_planar_boundary(lev, *w_star[lev]); }
         fill_planar_boundary(lev, *t_star[lev]);
         fill_planar_boundary(lev, *q_star[lev]);
         fill_planar_boundary(lev, *olen[lev]);
     } else {
         // The ordinary FillBoundary is unsafe for lateral walls: the
         // collapsed layout contains FABs from interior grids whose surface
         // parameters were never computed. Exchange only after all flux
         // iterations have written the selected face FABs, and copy the
         // communicated values back only to those FABs.
         fill_lateral_surface_parameter_ghosts(lev);
     }
 }

Here is the call graph for this function:

◆ update_mac_ptrs()

void SurfaceLayer::update_mac_ptrs	(	const int &	lev,
		amrex::Vector< amrex::Vector< amrex::MultiFab >> &	vars_old,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	Theta_prim,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	Qv_prim,
		amrex::Vector< std::unique_ptr< amrex::MultiFab >> &	Qr_prim
	)

inline

Update MOST-average field pointers.

Parameters

[in]	lev	level index
[in]	vars_old	old-time state variables
[in]	Theta_prim	primitive potential-temperature fields by level
[in]	Qv_prim	primitive water-vapor fields by level
[in]	Qr_prim	primitive rain-water fields by level

   {
       m_ma.update_field_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
   }

Here is the call graph for this function:

◆ update_pblh()

void SurfaceLayer::update_pblh	(	const int &	lev,
		amrex::Vector< amrex::Vector< amrex::MultiFab >> &	vars,
		amrex::MultiFab *	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices
	)

Wrapper around compute_pblh.

Parameters

[in]	lev	level index
[in,out]	vars	state variables used by the PBL-height calculation
[in]	z_phys_cc	cell-centered physical-height field
[in]	moisture_indices	indices for moisture components

Update PBL height using the configured estimator.

Parameters

[in]	lev	Current level
[in]	vars	Level-indexed state MultiFabs passed to the PBL height estimator
[in]	z_phys_cc	Cell-centered physical height used by the PBL height estimator
[in]	moisture_indices	Moisture component indices used by the PBL height estimator

 {
     if (pblh_type == PBLHeightCalcType::MYNN25) {
         MYNNPBLH estimator;
         compute_pblh(lev, vars, z_phys_cc, estimator, moisture_indices);
     } else if (pblh_type == PBLHeightCalcType::YSU || pblh_type == PBLHeightCalcType::MRF || pblh_type == PBLHeightCalcType::YSUNew)
     {
         //amrex::Error("YSU/MRF PBLH calc not implemented yet");
         // PBLH for MRF/YSU is computed inside ComputeDiffusivityMRF
         // and written back via set_pblh(). Nothing to do here.
     }
 }

◆ update_sfc_time_index()

void SurfaceLayer::update_sfc_time_index ( const amrex::Real & time )

Updates current time index for interpolating data from SFC/SST file.

Parameters

[in] time elapsed time

 {
     if (sfc.empty() || sfc[0].size() < 2) { return; }
  
     Real t1 = sfc[0][sfc_time_ind+1];
     while (elapsed_time >= t1)
     {
         int prev_index = sfc_time_ind;
         sfc_time_ind = std::min(sfc_time_ind + 1, int(sfc[0].size() - 2));
         t1 = sfc[0][sfc_time_ind+1];
         if (prev_index == sfc_time_ind) {
             break;
         }
     }
 }

◆ update_sst_ptr()

void SurfaceLayer::update_sst_ptr	(	const int	lev,
		const int	itime,
		amrex::MultiFab *	sst_ptr
	)

inline

Update one stored sea-surface-temperature pointer.

Parameters

[in]	lev	level index
[in]	itime	time-slice index
[in]	sst_ptr	sea-surface-temperature field pointer

                                                                                {
       AMREX_ALWAYS_ASSERT(itime >= 0 && itime < static_cast<int>(m_sst_lev[lev].size()));
       m_sst_lev[lev][itime] = sst_ptr;
   }

Here is the call graph for this function:

◆ update_surf_temp()

void SurfaceLayer::update_surf_temp ( const double & time )

inline

Update prescribed surface temperature from the configured heating rate.

Parameters

[in] time elapsed simulation time

   {
       // NOTE: this is a whole-domain setVal, so it overwrites the SST/TSK fill
       //       done earlier in update_fluxes. Coupled SST is applied after this
       //       call and therefore still wins on the water cells it covers.
       // An absolute surface temperature is written, converted, by fill_tsurf_with_abs_surf_temp
       if (surf_temp_is_absolute) { return; }
       if (surf_heating_rate != 0) {
           // Use the actual size of t_surf, not m_geom.size(), which is always
           // max_level+1 and so runs past the levels that exist.  t_surf is sized for
           // all levels up front but filled one level at a time, so we also have to
           // skip the entries that have not been allocated yet.
           int nlevs = static_cast<int>(t_surf.size());
           for (int lev = 0; lev < nlevs; lev++) {
               if (!t_surf[lev]) { continue; }
               t_surf[lev]->setVal(surf_temp + surf_heating_rate * static_cast<amrex::Real>(time));
               amrex::Print() << "Surface temp at t=" << time << ": "
                              << surf_temp + surf_heating_rate * time << std::endl;
           }
       }
   }

◆ update_tsk_ptr()

void SurfaceLayer::update_tsk_ptr	(	const int	lev,
		const int	itime,
		amrex::MultiFab *	tsk_ptr
	)

inline

Update one stored skin-temperature pointer.

Parameters

[in]	lev	level index
[in]	itime	time-slice index
[in]	tsk_ptr	skin-temperature field pointer

                                                                                {
       AMREX_ALWAYS_ASSERT(itime >= 0 && itime < static_cast<int>(m_tsk_lev[lev].size()));
       m_tsk_lev[lev][itime] = tsk_ptr;
   }

Here is the call graph for this function:

◆ use_sfc_fluxes()

bool SurfaceLayer::use_sfc_fluxes ( ) const

inline

Return whether prescribed surface fluxes are active.

1570 { return m_use_sfc_fluxes; }

Member Data Documentation

◆ cnk_a

amrex::Real SurfaceLayer::cnk_a {amrex::Real(0.0185)}

private

◆ custom_qstar

amrex::Real SurfaceLayer::custom_qstar {0}

private

◆ custom_rhosurf

amrex::Real SurfaceLayer::custom_rhosurf {0}

private

◆ custom_tstar

amrex::Real SurfaceLayer::custom_tstar {0}

private

◆ custom_ustar

amrex::Real SurfaceLayer::custom_ustar {0}

private

◆ default_land_surf_moist

amrex::Real SurfaceLayer::default_land_surf_moist {zero}

private

Referenced by make_SurfaceLayer_at_level().

◆ default_land_surf_temp

amrex::Real SurfaceLayer::default_land_surf_temp {amrex::Real(300.)}

private

Referenced by make_SurfaceLayer_at_level().

◆ depth

amrex::Real SurfaceLayer::depth {amrex::Real(30.0)}

private

◆ flux_type

FluxCalcType SurfaceLayer::flux_type {FluxCalcType::MOENG}

Referenced by make_SurfaceLayer_at_level().

◆ inv_Cmu2

amrex::Real SurfaceLayer::inv_Cmu2 = zero

private

◆ m_Cd

amrex::Real SurfaceLayer::m_Cd = zero

private

◆ m_Ch

amrex::Real SurfaceLayer::m_Ch = zero

private

◆ m_coupled_sst_lev

amrex::Vector<amrex::MultiFab*> SurfaceLayer::m_coupled_sst_lev

private

Referenced by make_SurfaceLayer_at_level(), and update_coupled_sst_ptr().

◆ m_coupled_sst_valid_lev

amrex::Vector<amrex::iMultiFab*> SurfaceLayer::m_coupled_sst_valid_lev

private

Referenced by make_SurfaceLayer_at_level(), and update_coupled_sst_ptr().

◆ m_Cq

amrex::Real SurfaceLayer::m_Cq = zero

private

◆ m_eb_vec

amrex::Vector<const eb_*> SurfaceLayer::m_eb_vec

private

◆ m_eddyDiffs_lev

amrex::Vector<amrex::MultiFab*> SurfaceLayer::m_eddyDiffs_lev

private

Referenced by make_SurfaceLayer_at_level().

◆ m_face

amrex::Orientation SurfaceLayer::m_face

private

Referenced by make_SurfaceLayer_at_level().

◆ m_final_low_time

double SurfaceLayer::m_final_low_time

private

◆ m_geom

amrex::Vector<amrex::Geometry> SurfaceLayer::m_geom

private

Referenced by make_SurfaceLayer_at_level().

◆ m_has_immersed_buildings

bool SurfaceLayer::m_has_immersed_buildings = false

private

◆ m_has_lsm_fluxes

bool SurfaceLayer::m_has_lsm_fluxes = false

private

Referenced by make_SurfaceLayer_at_level().

◆ m_has_lsm_tsurf

bool SurfaceLayer::m_has_lsm_tsurf = false

private

Referenced by make_SurfaceLayer_at_level(), and skin_temperature_conflict().

◆ m_Hwave_lev

amrex::Vector<amrex::MultiFab*> SurfaceLayer::m_Hwave_lev

private

Referenced by make_SurfaceLayer_at_level().

◆ m_ignore_sst

bool SurfaceLayer::m_ignore_sst = false

private

Referenced by make_SurfaceLayer_at_level().

◆ m_include_wstar

bool SurfaceLayer::m_include_wstar = false

private

Referenced by computes_w_star().

◆ m_lmask_lev

amrex::Vector<amrex::Vector<amrex::iMultiFab*> > SurfaceLayer::m_lmask_lev

private

Referenced by get_lmask(), and make_SurfaceLayer_at_level().

◆ m_low_time_interval

double SurfaceLayer::m_low_time_interval

private

◆ m_lsm_data_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceLayer::m_lsm_data_lev

private

Referenced by make_SurfaceLayer_at_level().

◆ m_lsm_data_name

amrex::Vector<std::string> SurfaceLayer::m_lsm_data_name

private

Referenced by make_SurfaceLayer_at_level().

◆ m_lsm_flux_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceLayer::m_lsm_flux_lev

private

Referenced by make_SurfaceLayer_at_level().

◆ m_lsm_flux_name

amrex::Vector<std::string> SurfaceLayer::m_lsm_flux_name

private

Referenced by make_SurfaceLayer_at_level().

◆ m_lsm_tsurf_indx

int SurfaceLayer::m_lsm_tsurf_indx = -1

private

Referenced by make_SurfaceLayer_at_level().

◆ m_lsm_tsurf_is_theta

bool SurfaceLayer::m_lsm_tsurf_is_theta = false

private

Referenced by make_SurfaceLayer_at_level().

◆ m_Lwave_lev

amrex::Vector<amrex::MultiFab*> SurfaceLayer::m_Lwave_lev

private

Referenced by make_SurfaceLayer_at_level().

◆ m_ma

MOSTAverage SurfaceLayer::m_ma

private

Referenced by get_mac_avg(), get_mac_avg_ptr(), get_mac_plane_avg(), get_num_mac_avg(), get_zref(), get_zref_min(), mac_avg_is_initialized(), mac_avg_is_time_averaged(), make_SurfaceLayer_at_level(), set_mac_avg_initialized(), set_mac_plane_avg(), and update_mac_ptrs().

◆ m_pblh_columns

amrex::Vector<PBLHColumns> SurfaceLayer::m_pblh_columns

private

◆ m_planar_bndry

amrex::Vector<PlanarBoundary> SurfaceLayer::m_planar_bndry

private

Referenced by make_SurfaceLayer_at_level().

◆ m_pp_prefix

std::string SurfaceLayer::m_pp_prefix

private

Referenced by make_SurfaceLayer_at_level(), and set_default_roughness().

◆ m_rdOcp

amrex::Real SurfaceLayer::m_rdOcp = RdoCp

private

◆ m_rotate

bool SurfaceLayer::m_rotate = false

private

Referenced by rotates_surface_fluxes(), and skin_temperature_conflict().

◆ m_skin_moisture_hs_lev

amrex::Vector<amrex::Real> SurfaceLayer::m_skin_moisture_hs_lev

private

Referenced by set_skin_moisture().

◆ m_skin_moisture_lev

amrex::Vector<const amrex::MultiFab*> SurfaceLayer::m_skin_moisture_lev

private

Referenced by set_skin_moisture().

◆ m_skin_moisture_resistances_lev

amrex::Vector<int> SurfaceLayer::m_skin_moisture_resistances_lev

private

Referenced by set_skin_moisture().

◆ m_skin_tsurf_lev

amrex::Vector<const amrex::MultiFab*> SurfaceLayer::m_skin_tsurf_lev

private

Referenced by set_skin_temperature().

◆ m_sst_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceLayer::m_sst_lev

private

Referenced by make_SurfaceLayer_at_level(), and update_sst_ptr().

◆ m_start_low_time

double SurfaceLayer::m_start_low_time

private

◆ m_surf_model

SurfaceModel* SurfaceLayer::m_surf_model = nullptr

private

◆ m_surface_layer_faces

amrex::GpuArray<int, AMREX_SPACEDIM*2> SurfaceLayer::m_surface_layer_faces {}

private

Referenced by set_surface_layer_faces().

◆ m_terrain_type

TerrainType SurfaceLayer::m_terrain_type

private

Referenced by make_SurfaceLayer_at_level(), and skin_temperature_conflict().

◆ m_tsk_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceLayer::m_tsk_lev

private

Referenced by make_SurfaceLayer_at_level(), and update_tsk_ptr().

◆ m_update_k_rans

bool SurfaceLayer::m_update_k_rans = false

private

◆ m_use_coupled_sst

bool SurfaceLayer::m_use_coupled_sst = false

private

Referenced by make_SurfaceLayer_at_level(), and set_coupled_sst_active().

◆ m_use_lsm_skin_temperature

bool SurfaceLayer::m_use_lsm_skin_temperature = false

private

Referenced by make_SurfaceLayer_at_level().

◆ m_use_sfc_fluxes

bool SurfaceLayer::m_use_sfc_fluxes = false

private

Referenced by make_SurfaceLayer_at_level(), and use_sfc_fluxes().

◆ m_use_sfc_sst

bool SurfaceLayer::m_use_sfc_sst = false

private

Referenced by make_SurfaceLayer_at_level().

◆ m_var_z0

bool SurfaceLayer::m_var_z0 {false}

private

Referenced by have_variable_sea_roughness(), and make_SurfaceLayer_at_level().

◆ m_z0_given

bool SurfaceLayer::m_z0_given {false}

private

Referenced by set_default_roughness().

◆ moist_type

MoistCalcType SurfaceLayer::moist_type {MoistCalcType::ADIABATIC}

Referenced by make_SurfaceLayer_at_level(), and skin_moisture_conflict().

◆ olen

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::olen

private

Referenced by get_olen(), and make_SurfaceLayer_at_level().

◆ pblh

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::pblh

private

Referenced by get_pblh(), and make_SurfaceLayer_at_level().

◆ pblh_type

PBLHeightCalcType SurfaceLayer::pblh_type {PBLHeightCalcType::None}

Referenced by computes_pblh().

◆ q_star

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

private

Referenced by get_q_star(), and make_SurfaceLayer_at_level().

◆ q_surf

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::q_surf

private

Referenced by get_q_surf(), make_SurfaceLayer_at_level(), and set_q_surf().

◆ rico_qsat_z0

amrex::Real SurfaceLayer::rico_qsat_z0 {amrex::Real(0.001)}

private

◆ rico_theta_z0

amrex::Real SurfaceLayer::rico_theta_z0 {amrex::Real(298.0)}

private

◆ rough_type_land

RoughCalcType SurfaceLayer::rough_type_land {RoughCalcType::CONSTANT}

◆ rough_type_sea

RoughCalcType SurfaceLayer::rough_type_sea {RoughCalcType::CHARNOCK}

Referenced by make_SurfaceLayer_at_level().

◆ sfc

amrex::Vector<amrex::Vector<amrex::Real> > SurfaceLayer::sfc

private

Referenced by make_SurfaceLayer_at_level().

◆ sfc_qflux

amrex::Real SurfaceLayer::sfc_qflux = zero

private

◆ sfc_tflux

amrex::Real SurfaceLayer::sfc_tflux = zero

private

◆ sfc_time_ind

int SurfaceLayer::sfc_time_ind = 0

private

◆ sfc_ustar

amrex::Real SurfaceLayer::sfc_ustar = zero

private

◆ smooth_flow_visc

bool SurfaceLayer::smooth_flow_visc {true}

private

◆ specified_rho_surf

bool SurfaceLayer::specified_rho_surf {false}

private

◆ surf_heating_rate

amrex::Real SurfaceLayer::surf_heating_rate {0}

private

Referenced by make_SurfaceLayer_at_level(), skin_temperature_conflict(), and update_surf_temp().

◆ surf_model_fluxes

bool SurfaceLayer::surf_model_fluxes = false

private

◆ surf_moist

amrex::Real SurfaceLayer::surf_moist {amrex::Real(-1.)}

private

◆ surf_moist_flux

amrex::Real SurfaceLayer::surf_moist_flux {0}

private

◆ surf_temp

amrex::Real SurfaceLayer::surf_temp {amrex::Real(-1.)}

private

Referenced by update_surf_temp().

◆ surf_temp_flux

amrex::Real SurfaceLayer::surf_temp_flux {0}

private

◆ surf_temp_is_absolute

bool SurfaceLayer::surf_temp_is_absolute {false}

private

Referenced by make_SurfaceLayer_at_level(), and update_surf_temp().

◆ surface_diagnostic_source

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::surface_diagnostic_source

private

Referenced by get_surface_diagnostic_source(), and make_SurfaceLayer_at_level().

◆ t_star

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

private

Referenced by get_t_star(), and make_SurfaceLayer_at_level().

◆ t_surf

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

private

Referenced by get_t_surf(), make_SurfaceLayer_at_level(), set_t_surf(), and update_surf_temp().

◆ theta_ref

amrex::Real SurfaceLayer::theta_ref = zero

private

◆ theta_type

ThetaCalcType SurfaceLayer::theta_type {ThetaCalcType::ADIABATIC}

Referenced by make_SurfaceLayer_at_level(), and skin_temperature_conflict().

◆ u_star

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

private

Referenced by get_u_star(), and make_SurfaceLayer_at_level().

◆ use_moisture

bool SurfaceLayer::use_moisture

private

Referenced by skin_moisture_conflict().

◆ use_surface_model

bool SurfaceLayer::use_surface_model = false

private

Referenced by make_SurfaceLayer_at_level(), and skin_temperature_conflict().

◆ w_star

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::w_star

private

Referenced by get_w_star(), and make_SurfaceLayer_at_level().

◆ z0_const

amrex::Real SurfaceLayer::z0_const {amrex::Real(0.1)}

private

Referenced by make_SurfaceLayer_at_level(), roughness(), and set_default_roughness().

◆ z_0

amrex::Vector<amrex::MultiFab> SurfaceLayer::z_0

private

Referenced by get_z0(), and make_SurfaceLayer_at_level().

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

Source/BoundaryConditions/ERF_SurfaceLayer.H
Source/BoundaryConditions/ERF_SurfaceLayer.cpp

Classes

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ FluxCalcType

◆ MoistCalcType

◆ PBLHeightCalcType

◆ RoughCalcType

◆ ThetaCalcType

Constructor & Destructor Documentation

◆ SurfaceLayer()

Member Function Documentation

◆ compute_fluxes() [1/2]

◆ compute_fluxes() [2/2]

◆ compute_pblh() [1/2]

◆ compute_pblh() [2/2]

◆ compute_sfc_params_from_lsm_fluxes()

◆ compute_SurfaceLayer_bcs() [1/2]

◆ compute_SurfaceLayer_bcs() [2/2]

◆ compute_SurfaceLayer_bcs_EB() [1/2]

◆ compute_SurfaceLayer_bcs_EB() [2/2]

◆ computes_pblh()

◆ computes_w_star()

◆ define_pblh_columns()

◆ fill_lateral_surface_parameter_ghosts()

◆ fill_pblh_from_coarser()

◆ fill_planar_boundary()

◆ fill_qsurf_with_qsat()

◆ fill_qsurf_with_skin_moisture()

◆ fill_tsurf_with_abs_surf_temp()

◆ fill_tsurf_with_coupled_sst()

◆ fill_tsurf_with_sfc_sst()

◆ fill_tsurf_with_skin_temperature()

◆ fill_tsurf_with_sst_and_tsk()

◆ get_lmask()

◆ get_lsm_tsurf()

◆ get_mac_avg()

◆ get_mac_avg_ptr()

◆ get_mac_plane_avg()

◆ get_num_mac_avg()

◆ get_olen()

◆ get_pblh()

◆ get_q_star()

◆ get_q_surf()

◆ get_surface_diagnostic_source()

◆ get_t_star()

◆ get_t_surf()

◆ get_u_star()

◆ get_w_star()

◆ get_z0()

◆ get_zref()

◆ get_zref_min()

◆ have_variable_sea_roughness()

◆ impose_SurfaceLayer_bcs()

◆ impose_SurfaceLayer_bcs_EB()

◆ init_tke_from_ustar()

◆ inject_pblh_from_coarser()

◆ interpolate_sfc_column()

◆ lmask_min_reduce()

◆ mac_avg_is_initialized()

◆ mac_avg_is_time_averaged()

◆ make_SurfaceLayer_at_level()

◆ read_cols()

◆ read_custom_roughness()

◆ rotates_surface_fluxes()

◆ roughness()

◆ set_coupled_sst_active()

◆ set_default_roughness()

◆ set_immersed_surface_params()

◆ set_mac_avg_initialized()

◆ set_mac_plane_avg()

◆ set_pblh()

◆ set_q_surf()

◆ set_sfc_params_from_file_fluxes()

◆ set_skin_moisture()