SurfaceLayer Class Reference

#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:

Public Types
enum class	FluxCalcType {   MOENG = 0 , DONELAN , 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 (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 MeshType &a_mesh_type, const TerrainType &a_terrain_type, const TurbChoice &a_turb_choice, 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	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)
void	fill_tsurf_with_coupled_sst (const int &lev)
void	fill_qsurf_with_qsat (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)
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)
template<typename FluxIter >
void	compute_fluxes (const int &lev, const int &max_iters, MultiFab &cons_in, 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 Attributes
amrex::Vector< amrex::Geometry >	m_geom
bool	m_rotate = false
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	cnk_visc {false}
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
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< 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.
DONELAN	Donelan 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
    DONELAN,   ///< Donelan 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 ( 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 MeshType &  a_mesh_type, const TerrainType &  a_terrain_type, const TurbChoice &  a_turb_choice, 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]	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]	a_mesh_type	mesh type
[in]	a_terrain_type	terrain representation
[in]	a_turb_choice	turbulence-model options
[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_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_ma(geom, (z_phys_nd[0] != nullptr), a_pp_prefix, a_mesh_type, a_terrain_type, eb_vec)
  {
      // We have a moisture model if Qv_prim is a valid pointer
      use_moisture = (Qv_prim[0].get());
 
      // Get roughness
      amrex::ParmParse pp("erf");
      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 == "donelan") {
              flux_type = FluxCalcType::DONELAN;
          } else 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!");
          }
      }
 
      // 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!");
      }
 
      // 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 (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!");
      }
 
      // 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);
          pp.queryAdd("most.charnock_viscosity", cnk_visc);
          if (cnk_a > 0) {
              amrex::Print() << "If there is water, Charnock relation with C_a="
                             << cnk_a << (cnk_visc ? " and viscosity" : "")
                             << " will be used" << std::endl;
          } else {
              amrex::Print() << "If there is water, Charnock relation with variable "
                  "Charnock parameter (COARE3.0)"
                             << (cnk_visc ? " and viscosity" : "") << " 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)"
                         << (cnk_visc ? " and viscosity" : "") << " 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);
      } 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
      // BC 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) {
          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 &	cons_in,
		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,2); // potential temperature
    const auto *const qvm_ptr  = m_ma.get_average(lev,3); // water vapor mixing ratio
    const auto *const tvm_ptr  = m_ma.get_average(lev,4); // virtual potential temperature
    const auto *const umm_ptr  = m_ma.get_average(lev,5); // horizontal 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 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 gtbx = mfi.tilebox(IntVect(0),ng);
 
        if (!l_use_eb && gtbx.smallEnd(2) != klo) { continue; }
 
        if (!l_use_eb) { gtbx.makeSlab(2,klo); }
 
        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 zref_arr = zref_ptr->array(mfi);
        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 ) noexcept
            {
                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, 0, 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);
                }
            });
        // 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

{
    est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
                     vars[lev][Vars::cons],m_lmask_lev[lev][0],
                     moisture_indices);
}

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

{
    Real eps = std::numeric_limits<amrex::Real>::epsilon();
    bool has_moisture = use_moisture;
    const int klo = m_geom[lev].Domain().smallEnd(2);
    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);
 
        // 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);
                }
                olen_arr(i,j,0)   = ( u_star_arr(i,j,0) * u_star_arr(i,j,0) * Thv ) /
                                    ( KAPPA * CONST_GRAV * t_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 klo = m_geom[lev].Domain().smallEnd(2);
    const auto& dxInv = m_geom[lev].InvCellSizeArray();
    for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
    {
        // 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
        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 = (m_rotate) ? Tau_lev[TauType::tau12]->array(mfi) : Array4<Real>{};
        auto t21_arr = (m_rotate) ? Tau_lev[TauType::tau21]->array(mfi) : Array4<Real>{};
 
        auto hfx1_arr = (m_rotate) ? xheat_flux->array(mfi) : Array4<Real>{};
        auto hfx2_arr = (m_rotate) ? yheat_flux->array(mfi) : Array4<Real>{};
        auto qfx1_arr = (m_rotate && xqv_flux) ? xqv_flux->array(mfi) : Array4<Real>{};
        auto qfx2_arr = (m_rotate && yqv_flux) ? 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 t_mean     = m_ma.get_average(lev,2);
        const auto *const q_mean     = m_ma.get_average(lev,3);
        const auto *const u_mag_mean = m_ma.get_average(lev,5);
 
        const auto um_arr  = u_mean->array(mfi);
        const auto vm_arr  = v_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);
 
        // 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
        //============================================================================
        Box bx = mfi.tilebox();
        if (bx.smallEnd(2) != klo) { continue; }
        bx.makeSlab(2,klo);
        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,
                                                 cons_arr, velx_arr, vely_arr,
                                                 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,0) = 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, Tflux, dxInv, zphys_arr,
                                   hfx1_arr, hfx2_arr, hfx3_arr);
            } else {
                hfx3_arr(i,j,klo) = 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,
                                                     cons_arr, velx_arr, vely_arr,
                                                     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, Qflux, dxInv, zphys_arr,
                                       qfx1_arr, qfx2_arr, qfx3_arr);
                } else {
                    qfx3_arr(i,j,k) = Qflux;
                }
            });
        } // custom
 
        if (!rotate) {
            // Rho*u flux
            //============================================================================
            Box bxx = surroundingNodes(bx,0);
            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,
                                                 cons_arr, velx_arr, vely_arr,
                                                 umm_arr, um_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,
                                                       cons_arr, velx_arr, vely_arr,
                                                       umm_arr, um_arr, u_star_arr);
                }
 
                t13_arr(i,j,k) = stressx;
                if (t31_arr) { t31_arr(i,j,k) = stressx; }
            });
 
            // Rho*v flux
            //============================================================================
            Box bxy = surroundingNodes(bx,1);
            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,
                                                 cons_arr, velx_arr, vely_arr,
                                                 umm_arr, vm_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,
                                                       cons_arr, velx_arr, vely_arr,
                                                       umm_arr, vm_arr, u_star_arr);
                }
 
                t23_arr(i,j,k) = stressy;
                if (t32_arr) { t32_arr(i,j,k) = stressy; }
            });
        } else {
            // 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,
                                                       cons_arr, velx_arr, vely_arr,
                                                       umm_arr, um_arr, u_star_arr);
                rotate_stress_tensor(i, j, k, 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 ||
            flux_type == FluxCalcType::DONELAN) {
            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);
                }
                olen_arr(i,j,0)   = ( u_star_arr(i,j,0) * u_star_arr(i,j,0) * Thv ) /
                                    ( KAPPA * CONST_GRAV * t_star_arr(i,j,0) );
            });
        }
 
    } // mfiter
 
    surface_diagnostic_source[lev]->FillBoundary(m_geom[lev].periodicity());
}

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

{
    // 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 t_mean     = m_ma.get_average(lev,2);
        // const auto *const q_mean     = m_ma.get_average(lev,3);
        const auto *const u_mag_mean = m_ma.get_average(lev,5);
 
        const auto um_arr  = u_mean->array(mfi);
        const auto vm_arr  = v_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);
 
        // 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,
                                                    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,
                                                        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,
                                                        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,
                                                        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,
                                                        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,
                                                        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,
                                                        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); }

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; }

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
 
    // Populate q_surf with qsat over water
    auto dz = m_geom[lev].CellSize(2);
    const int klo = m_geom[lev].Domain().smallEnd(2);
    for (MFIter mfi(*q_surf[lev]); mfi.isValid(); ++mfi)
    {
        Box gtbx = mfi.growntilebox();
 
        if (gtbx.smallEnd(2) != klo) { 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> {};
 
        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) {
                auto deltaZ = (z_arr) ? Compute_Zrel_AtCellCenter(i,j,k,z_arr) :
                                        myhalf*dz;
                auto Rho  = cons_arr(i,j,k,Rho_comp);
                auto RTh  = cons_arr(i,j,k,RhoTheta_comp);
                auto Qv   = cons_arr(i,j,k,RhoQ1_comp) / Rho;
                auto P_cc = getPgivenRTh(RTh, Qv);
                P_cc += Rho*CONST_GRAV*deltaZ;
                P_cc *= Real(0.01);
                erf_qsatw(t_surf_arr(i,j,k), P_cc, q_surf_arr(i,j,k));
            }
        });
    }
    q_surf[lev]->FillBoundary(m_geom[lev].periodicity());
}

Here is the call graph for this function:

fill_tsurf_with_coupled_sst()

void SurfaceLayer::fill_tsurf_with_coupled_sst

(

const int &

lev

)

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);
 
    // 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) { continue; }
 
        // NOTE: the coupled lane does not carry lateral ghost cells, so clamp
        //       into the valid box exactly as get_lsm_tsurf does. FillBoundary
        //       in update_fluxes picks up the interior 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 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>{};
 
        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) { return; }
 
            int li = amrex::min(amrex::max(i, i_lo), i_hi);
            int lj = amrex::min(amrex::max(j, j_lo), j_hi);
 
            if (has_valid && valid_arr(li,lj,k) == 0) { return; }
 
            t_surf_arr(i,j,k) = coupled_sst_arr(li,lj,k);
        });
    }
}

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
	)

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);
 
    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);
        auto lmask_arr  = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
                                                  Array4<int> {};
 
        ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 0;
            if (!is_land) {
                t_surf_arr(i,j,k) = sfc_sst;
            }
        });
    }
 
    t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
}

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;
                }
            });
        }
    }
    t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
}

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);
                t_surf_arr(i,j,k) = lsm_arr(li,lj,k);
            }
        });
    }
}

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(); }

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(); }

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(); }

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(); }

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(); }

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(); }

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(); }

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(); }

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(); }

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]; }

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) {
        moeng_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::DONELAN) {
        donelan_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::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

{
    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
    const int klo = m_geom[lev].Domain().smallEnd(2);
    Box bx_lo = u_star[lev]->boxArray().minimalBox();
    FArrayBox u_star_lo(bx_lo, 1); u_star_lo.setVal<RunOn::Device>(0);
    FArrayBox z_surf_lo(bx_lo, 1); z_surf_lo.setVal<RunOn::Device>(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) );
        });
    }
    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("erf");
 
      // 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) {
              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);
      }
 
      // 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 sst(K); prescribed-flux mode also uses H, LE, USTAR.
      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");
          }
 
          // 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; }
 
      // 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
      }
 
      // Attributes for MFs and FABs
      //--------------------------------------------------------
      // Create a 2D ba for planar terrain, 3D for EB terrain
      amrex::BoxArray ba = mf.boxArray();
      amrex::BoxArray ba_flux;
      amrex::IntVect ng{1,1,0};
 
      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(2,0); }
          ba_flux = amrex::BoxArray(std::move(bl2d));
      }
 
      const amrex::DistributionMapping& dm = mf.DistributionMap();
      const int ncomp = 1;
 
      // 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);
 
      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?
 
      bool use_sst = (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]);
      bool use_tsk = (!m_tsk_lev[lev].empty() && m_tsk_lev[lev][0]);
      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_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_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); }

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

{
    // 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 (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]) {
        fill_tsurf_with_sst_and_tsk(lev, elapsed_time_since_start_low);
    }
    if (m_use_sfc_sst) {
        // Set tsurf to time varying SST from sfc file
        fill_tsurf_with_sfc_sst(lev, elapsed_time);
    }
 
    // 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.
    fill_tsurf_with_coupled_sst(lev);
 
    // 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) { get_lsm_tsurf(lev); }
 
    // Fill interior ghost cells
    t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
 
    // 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);
                    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, cnk_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, 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,
                                                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,
                                                        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, cnk_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, 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,
                                                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,
                                                        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, cnk_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);
                    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);
                    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);
                    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) {
        const bool use_ref_theta = (theta_ref > 0);
        const Real l_inv_theta0  = (use_ref_theta) ? one / 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);
 
            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);
                }
            });
        }
    }
 
    u_star[lev]->FillBoundary(m_geom[lev].periodicity());
    t_star[lev]->FillBoundary(m_geom[lev].periodicity());
    q_star[lev]->FillBoundary(m_geom[lev].periodicity());
      olen[lev]->FillBoundary(m_geom[lev].periodicity());
}

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) {
        amrex::Error("YSU/MRF PBLH calc not implemented yet");
    }
}

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

cnk_visc

bool SurfaceLayer::cnk_visc {false}

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_final_low_time

double SurfaceLayer::m_final_low_time

private

m_geom

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

private

m_has_lsm_fluxes

bool SurfaceLayer::m_has_lsm_fluxes = false

private

Referenced by make_SurfaceLayer_at_level().

m_has_lsm_tsurf

bool SurfaceLayer::m_has_lsm_tsurf = false

private

Referenced by make_SurfaceLayer_at_level().

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

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