SurfaceLayer Class Reference

#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:

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={})
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	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)
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)
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	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_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)
void	update_pblh (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars, amrex::MultiFab *z_phys_cc, const MoistureComponentIndices &moisture_indices)
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_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::Real	get_zref (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)
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)}
amrex::Real	default_land_surf_temp {amrex::Real(300.)}
amrex::Real	surf_temp {amrex::Real(-1.)}
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_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
amrex::Real	m_rdOcp = RdoCp
MOSTAverage	m_ma
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< 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

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

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 = {}  )

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

                                                                 {})
  : 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_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)
  {
      // 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);
      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)");
      }
 
      // Get surface temperature.  surf_temp and surf_moist are declared with negative
      // sentinels (see below) so that "did the user set this" is a property of the
      // value rather than of the queryAdd return value, which only reports whether the
      // key existed before the call and so stops being meaningful once anything has
      // parsed the key.  Both most.surf_temp and most.surf_moist are also parsed by
      // ERF_InputSoundingData.H, so the two sites would poison each other otherwise.
      pp.queryAdd("most.surf_temp", surf_temp);
      const bool erf_st = (surf_temp > amrex::Real(0));
      if (erf_st) { default_land_surf_temp  = 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);
          pp.get("most.surf_temp", default_land_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];
 
    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 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);
 
        // 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<Real> {};
        auto lsm_q_flux_arr = Array4<Real> {};
        auto lsm_tau13_arr  = Array4<Real> {};
        auto lsm_tau23_arr  = Array4<Real> {};
        // compute_sfc_params_from_lsm_fluxes consumes signed kinematic stress
        // components; their vector magnitude determines u_star^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]->array(mfi); }
            if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
            if (toLower(m_lsm_flux_name[n]) == "tau13")  { lsm_tau13_arr  = m_lsm_flux_lev[lev][n]->array(mfi); }
            if (toLower(m_lsm_flux_name[n]) == "tau23")  { lsm_tau23_arr  = m_lsm_flux_lev[lev][n]->array(mfi); }
        }
 
        ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
        {
            int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
            // Skip cells the LSM did not have a valid flux (lsm_undefined).
            if (is_land && lsm_t_flux_arr && lsm_t_flux_arr(i,j,0) < lsm_undefined) {
                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(i,j,0)*lsm_tau13_arr(i,j,0)
                                    + lsm_tau23_arr(i,j,0)*lsm_tau23_arr(i,j,0) );
                u_star_arr(i,j,0) = amrex::max(std::sqrt(tau),eps);
                if (lsm_t_flux_arr(i,j,0)>=zero) {
                    t_star_arr(i,j,0) = amrex::min(-lsm_t_flux_arr(i,j,0) / u_star_arr(i,j,0),-eps);
                } else {
                    t_star_arr(i,j,0) = amrex::max(-lsm_t_flux_arr(i,j,0) / u_star_arr(i,j,0),eps);
                }
                if (lsm_q_flux_arr(i,j,0)>=zero) {
                    q_star_arr(i,j,0) = amrex::min(-lsm_q_flux_arr(i,j,0) / u_star_arr(i,j,0),-eps);
                } else {
                    q_star_arr(i,j,0) = amrex::max(-lsm_q_flux_arr(i,j,0) / 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
	)

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
	)

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<Real> {};
        auto soil_t_flux_arr = Array4<Real> {};
        auto lsm_q_flux_arr = Array4<Real> {};
        auto lsm_tau13_arr  = Array4<Real> {};
        auto lsm_tau23_arr  = Array4<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]->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]->array(mfi); }
            if (toLower(m_lsm_flux_name[n]) == "tau13")       { lsm_tau13_arr   = m_lsm_flux_lev[lev][n]->array(mfi); }
            if (toLower(m_lsm_flux_name[n]) == "tau23")       { lsm_tau23_arr   = m_lsm_flux_lev[lev][n]->array(mfi); }
        }
 
        const bool has_lsm_t_flux = static_cast<bool>(lsm_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();
 
        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 = (lsm_t_flux_arr) ? (lsm_t_flux_arr(i,j,0) < lsm_undefined) :
                                                              false;
            const bool has_land_and_flux = (is_land == 1 && lsm_flux_is_valid);
            if (lsm_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) * lsm_t_flux_arr(i,j,0);
            } else if (is_land == 2) { // no temperature flux within buildings
                Tflux = zero;
            } 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 (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 = (lsm_q_flux_arr) ? (lsm_q_flux_arr(i,j,0) < lsm_undefined) :
                                                                  false;
                const bool has_land_and_flux = (is_land == 1 && lsm_flux_is_valid);
                if (lsm_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) * lsm_q_flux_arr(i,j,0);
                } else if (is_land == 2) { // no moisture flux within buildings
                    Qflux = zero;
                } 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.
                }
 
                // 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>(lsm_tau13_arr), is_land_hi == 1,
                    lsm_tau13_arr ? lsm_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>(lsm_tau13_arr), is_land_lo == 1,
                    lsm_tau13_arr ? lsm_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 (lsm_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, lsm_tau13_arr(i-1,j,0), lsm_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>(lsm_tau23_arr), is_land_hi == 1,
                    lsm_tau23_arr ? lsm_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>(lsm_tau23_arr), is_land_lo == 1,
                    lsm_tau23_arr ? lsm_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 (lsm_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, lsm_tau23_arr(i,j-1,0), lsm_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.

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

{ 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);
    for (int ib = 0; ib < ba.size(); ++ib) {
        if (ba[ib].smallEnd(2) != k_ground) { cols.needed = true; }
    }
    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.");
    }
 
    // The runs of cells in z, each as one box; those 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.
    const BoxArray ba_joined = join_boxes_stacked_in_z(ba);
    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_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

{
    // NOTE: We have already tested a moisture model exists
    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 Real rdOcp = m_rdOcp;
    amrex::Gpu::DeviceScalar<int> d_conversion_failed(0);
    int* conversion_failed = d_conversion_failed.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);
                const Real rho_qv = cons_arr(i,j,k,RhoQ1_comp);
                if (!amrex::Math::isfinite(rho) || rho <= Real(0.0) ||
                    !amrex::Math::isfinite(rho_theta) || !amrex::Math::isfinite(rho_qv)) {
                    if (authoritative) {
                        amrex::Gpu::Atomic::Max(conversion_failed, 1);
                    }
                    return;
                }
                const Real qv = moist ? rho_qv / rho : Real(0.0);
                if (!amrex::Math::isfinite(qv)) {
                    if (authoritative) {
                        amrex::Gpu::Atomic::Max(conversion_failed, 1);
                    }
                    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 if (authoritative) {
                    amrex::Gpu::Atomic::Max(conversion_failed, 1);
                }
            }
        });
    }
    amrex::Gpu::streamSynchronize();
    int conversion_failed_host = d_conversion_failed.dataValue();
    amrex::ParallelDescriptor::ReduceIntMax(conversion_failed_host);
    if (conversion_failed_host != 0) {
        amrex::Abort("SurfaceLayer fill_qsurf_with_qsat: failed to convert the authoritative "
                     "surface potential temperature to absolute temperature.");
    }
    fill_planar_boundary(lev, *q_surf[lev]);
}

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<int> d_conversion_failed(0);
    int* conversion_failed = d_conversion_failed.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)) {
                amrex::Gpu::Atomic::Max(conversion_failed, 1);
                return;
            }
            Real qv = Real(0.0);
            if (moist) {
                if (!have_rho_qv) {
                    amrex::Gpu::Atomic::Max(conversion_failed, 1);
                    return;
                }
                const Real rho_qv = cons_arr(li,lj,klo,RhoQ1_comp);
                if (!amrex::Math::isfinite(rho_qv)) {
                    amrex::Gpu::Atomic::Max(conversion_failed, 1);
                    return;
                }
                qv = rho_qv / rho;
                if (!amrex::Math::isfinite(qv)) {
                    amrex::Gpu::Atomic::Max(conversion_failed, 1);
                    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 {
                amrex::Gpu::Atomic::Max(conversion_failed, 1);
            }
        });
    }
    amrex::Gpu::streamSynchronize();
    int conversion_failed_host = d_conversion_failed.dataValue();
    amrex::ParallelDescriptor::ReduceIntMax(conversion_failed_host);
    if (conversion_failed_host != 0) {
        amrex::Abort("SurfaceLayer fill_tsurf_with_coupled_sst: failed to convert the authoritative "
                     "coupled sea-surface temperature to potential temperature.");
    }
}

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<int> d_conversion_failed(0);
    int* conversion_failed = d_conversion_failed.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)) {
                    amrex::Gpu::Atomic::Max(conversion_failed, 1);
                    return;
                }
                Real qv = Real(0.0);
                if (moist) {
                    if (!have_rho_qv) {
                        amrex::Gpu::Atomic::Max(conversion_failed, 1);
                        return;
                    }
                    const Real rho_qv = cons_arr(li,lj,klo,RhoQ1_comp);
                    if (!amrex::Math::isfinite(rho_qv)) {
                        amrex::Gpu::Atomic::Max(conversion_failed, 1);
                        return;
                    }
                    qv = rho_qv / rho;
                    if (!amrex::Math::isfinite(qv)) {
                        amrex::Gpu::Atomic::Max(conversion_failed, 1);
                        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 {
                    amrex::Gpu::Atomic::Max(conversion_failed, 1);
                }
            }
        });
    }
 
    amrex::Gpu::streamSynchronize();
    int conversion_failed_host = d_conversion_failed.dataValue();
    amrex::ParallelDescriptor::ReduceIntMax(conversion_failed_host);
    if (conversion_failed_host != 0) {
        amrex::Abort("SurfaceLayer fill_tsurf_with_sfc_sst: failed to convert the authoritative "
                     "sea-surface temperature to potential temperature.");
    }
 
    fill_planar_boundary(lev, *t_surf[lev]);
}

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) {
            const auto tsk_arr = m_tsk_lev[lev][n_time_lo]->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 && !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_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

{ return m_lmask_lev[lev][0]; }

get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf

(

const int &

lev

)

Fill surface temperature from the LSM surface-temperature field.

Parameters

[in]

lev

level index

Fill surface temperature from land-surface-model data.

Parameters

[in]

lev

Current level

{
    const int klo = m_geom[lev].Domain().smallEnd(2);
    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.
        Box vbx  = mfi.validbox();
        int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
        int j_lo = vbx.smallEnd(1); 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 lsm_arr = m_lsm_data_lev[lev][m_lsm_tsurf_indx]->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) {
                int li = amrex::min(amrex::max(i, i_lo), i_hi);
                int lj = amrex::min(amrex::max(j, j_lo), j_hi);
                const Real lsm_value = lsm_arr(li,lj,k);
                if (amrex::Math::isfinite(lsm_value) && lsm_value > Real(0.0) &&
                    lsm_value < lsm_undefined) {
                    t_surf_arr(i,j,k) = lsm_value;
                }
            }
        });
    }
}

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

{ 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

{ return m_ma.get_plane_average(lev); }

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.

{ return m_ma.get_navg(); }

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

{ 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

{ 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

{ 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

{ 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

{ 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

{ 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

{ 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

{ 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

{ 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

{ return &z_0[lev]; }

Referenced by ERF::FillPlot2DVars().

Here is the caller graph for this function:

get_zref()

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

inline

Return the minimum reference height for one level.

Parameters

[in]

lev

level index

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

{ return m_var_z0; }

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
	)

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 = 0.1; // 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);
    } 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);
    } 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);
    } 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);
    } else if (flux_type == FluxCalcType::CUSTOM) {
        const bool fluxes_include_rho = specified_rho_surf || m_use_sfc_fluxes;
        custom_flux flux_comp(fluxes_include_rho);
        compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                 xheat_flux, yheat_flux, zheat_flux,
                                 xqv_flux, yqv_flux, zqv_flux,
                                 z_phys, flux_comp);
    } 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:

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

{ return m_ma.time_avg_is_initialized(lev); }

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.

{ return m_ma.do_time_averaging(); }

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_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;
          }
      }
 
      // 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]);
          if (lc_name == "theta" || lc_name == "t_surf") {
              m_has_lsm_tsurf  = true;
              m_lsm_tsurf_indx = n;
          }
      }
      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.");
 
      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) {
          // 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);
            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());
        }
 
        // 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);
 
        // 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:

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

{ m_use_coupled_sst = active; }

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

{ m_ma.set_time_avg_initialized(lev); }

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

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

Here is the call graph for this function:

set_pblh()

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

{
    AMREX_ASSERT(pblh[lev]);
    amrex::MultiFab::Copy(*pblh[lev], pblh_in, 0, 0, 1, 0);
    fill_planar_boundary(lev, *pblh[lev]);
}

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

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

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

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

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
    if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
        update_surf_temp(elapsed_time_since_start_low);
    }
 
    // 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 qsurf with qsat over sea
    if (use_moisture) {
        fill_qsurf_with_qsat(lev, cons_in, z_phys_nd);
    }
 
    // Update land surface temp if we have a valid pointer
    if (m_has_lsm_tsurf && zlo) get_lsm_tsurf(lev);
 
    // 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);
 
    // 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;
 
            // overwrite the custom_ustar/tstar/qstar values with the new values and
            // use the existing pathway to set u*,t*,q* with or without a custom_rhosurf
            // note - when m_use_sfc_fluxes=true, custom_flux has specified_rho_surf=true,
            // so there is no rho factor here
            custom_ustar = std::sqrt(sfc_ustar); // convert tau from file to u*
            custom_tstar = sfc_tflux / Cp_d;
            custom_qstar = sfc_qflux / L_v;
        }
 
        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);
        }
    }
 
    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 (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

                                                                              {
      m_sst_lev[lev][itime] = sst_ptr;
  }

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

                                                                              {
      m_tsk_lev[lev][itime] = tsk_ptr;
  }

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

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_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(), 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().

m_rdOcp

amrex::Real SurfaceLayer::m_rdOcp = RdoCp

private

m_rotate

bool SurfaceLayer::m_rotate = false

private

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

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_sfc_fluxes

bool SurfaceLayer::m_use_sfc_fluxes = false

private

Referenced by make_SurfaceLayer_at_level().

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

moist_type

MoistCalcType SurfaceLayer::moist_type {MoistCalcType::ADIABATIC}

Referenced by make_SurfaceLayer_at_level().

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(), and update_surf_temp().

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

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

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

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

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