#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:

[legend]

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_qsurf_with_qsat (const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)

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)

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)

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)

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)

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	surf_heating_rate {0}

amrex::Real	surf_temp_flux {0}

amrex::Real	default_land_surf_moist {zero}

amrex::Real	surf_moist

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

bool	m_has_ocean_lsm_tsurf = false

int	m_lsm_tsurf_indx = -1

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

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

RayleighDampingType::None

@ None

◆ 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

                                                                  {})
   : 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.query("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.query("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.query("most.include_wstar", m_include_wstar);
  
       std::string pblh_string_in;
       std::string pblh_string{"none"};
       auto read_pblh = pp.query("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!");
       }
  
       // Get surface temperature
       auto erf_st = pp.query("most.surf_temp", surf_temp);
       if (erf_st) { default_land_surf_temp  = surf_temp; }
  
       // Get surface moisture
       bool erf_sq = false;
       if (use_moisture) { erf_sq = pp.query("most.surf_moist", surf_moist); }
       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.query("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.query("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.query("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.query("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.query("most.rico.theta_z0", rico_theta_z0);
           pp.query("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.query("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.query("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.query("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.query("most.charnock_constant", cnk_a);
           pp.query("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.query("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.query("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_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]	most_flux	structure to iteratively compute ustar and tstar

 {
     // 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_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
	)

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

 {
     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_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,out]	mfs	MultiFabs to populate
[in]	eddyDiffs	Diffusion coefficients from turbulence model
[in]	flux_comp	structure to compute 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>{};
  
         // 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 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]) == "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);
                 if (lsm_t_flux_arr) {
                     lsm_t_flux_arr(i,j,0) = 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);
                     if (lsm_q_flux_arr) {
                         lsm_q_flux_arr(i,j,0) = Qflux / cons_arr(i,j,k,Rho_comp);
                     }
                 }
  
                 // 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;
                     if (!has_land_and_flux_lo) { lsm_tau13_arr(i-1,j,0) = result.low_kinematic_cache; }
                     if (!has_land_and_flux_hi) { lsm_tau13_arr(i  ,j,0) = result.high_kinematic_cache; }
                 } 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);
                     if (lsm_tau13_arr) {
                         // LSM tau13 is cell-centered kinematic stress [m2 s-2].
                         // Tau_lev tau13 is a face-centered conservative stress [N m-2].
                         // MOST compute_u_flux returns the conservative face stress.
                         lsm_tau13_arr(i  ,j,0) = surface_layer_stress::conservative_to_kinematic_stress(
                             stressx, cons_arr(i  ,j,k,Rho_comp));
                         lsm_tau13_arr(i-1,j,0) = surface_layer_stress::conservative_to_kinematic_stress(
                             stressx, cons_arr(i-1,j,k,Rho_comp));
                     }
                 }
  
                 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;
                     if (!has_land_and_flux_lo) { lsm_tau23_arr(i,j-1,0) = result.low_kinematic_cache; }
                     if (!has_land_and_flux_hi) { lsm_tau23_arr(i,j  ,0) = result.high_kinematic_cache; }
                 } 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);
                     if (lsm_tau23_arr) {
                         // LSM tau23 is cell-centered kinematic stress [m2 s-2].
                         // Tau_lev tau23 is a face-centered conservative stress [N m-2].
                         // MOST compute_v_flux returns the conservative face stress.
                         lsm_tau23_arr(i,j  ,0) = surface_layer_stress::conservative_to_kinematic_stress(
                             stressy, cons_arr(i,j  ,k,Rho_comp));
                         lsm_tau23_arr(i,j-1,0) = surface_layer_stress::conservative_to_kinematic_stress(
                             stressy, cons_arr(i,j-1,k,Rho_comp));
                     }
                 }
  
                 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);
                 });
         }
     } // 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_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,out]	mfs	MultiFabs to populate
[in]	eddyDiffs	Diffusion coefficients from turbulence model
[in]	flux_comp	structure to compute fluxes

 {
     // Get EB flags for all centerings
     const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
     const auto& u_flags = m_eb_vec[lev]->get_u_const_factory()->getMultiEBCellFlagFab();
     const auto& v_flags = m_eb_vec[lev]->get_v_const_factory()->getMultiEBCellFlagFab();
     const auto& w_flags = m_eb_vec[lev]->get_w_const_factory()->getMultiEBCellFlagFab();
  
     const auto& cc_vfrac = m_eb_vec[lev]->get_const_factory()->getVolFrac();
     const auto& u_vfrac = m_eb_vec[lev]->get_u_const_factory()->getVolFrac();
     const auto& v_vfrac = m_eb_vec[lev]->get_v_const_factory()->getVolFrac();
     const auto& w_vfrac = m_eb_vec[lev]->get_w_const_factory()->getVolFrac();
     const auto& cc_bnorm = m_eb_vec[lev]->get_const_factory()->getBndryNormal();
     const auto& u_bnorm = m_eb_vec[lev]->get_u_const_factory()->getBndryNorm();
     const auto& v_bnorm = m_eb_vec[lev]->get_v_const_factory()->getBndryNorm();
     const auto& w_bnorm = m_eb_vec[lev]->get_w_const_factory()->getBndryNorm();
  
     // 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);
         auto const bnorm_arr = cc_bnorm.const_array(mfi);
         auto const u_bnorm_arr = u_bnorm.const_array(mfi);
         auto const v_bnorm_arr = v_bnorm.const_array(mfi);
         auto const w_bnorm_arr = w_bnorm.const_array(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 - 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 (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,
                                                         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:

◆ 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
	)

 {
     // 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_sst_and_tsk()

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

 {
     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

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

◆ get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf ( const int & lev )

 {
     const int klo = m_geom[lev].Domain().smallEnd(2);
     const bool has_sea_tsurf = (m_has_ocean_lsm_tsurf &&
                                 amrex::toLower(m_lsm_data_name[m_lsm_tsurf_indx]) == "t_surf");
     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 ((!has_sea_tsurf && is_land) ||
                 (has_sea_tsurf && !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 m_ma.get_average(lev, comp);
   }

Here is the call graph for this function:

◆ get_olen()

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

inline

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

◆ get_pblh()

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

inline

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

◆ get_q_star()

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

inline

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

◆ get_q_surf()

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

inline

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

◆ get_surface_diagnostic_source()

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

inline

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

◆ get_t_star()

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

inline

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

◆ get_t_surf()

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

inline

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

◆ get_u_star()

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

inline

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

◆ get_w_star()

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

inline

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

◆ get_z0()

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

inline

652 { return &z_0[lev]; }

◆ get_zref()

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

inline

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

654 { return m_var_z0; }

SurfaceLayer::m_var_z0

bool m_var_z0

Definition: ERF_SurfaceLayer.H:750

◆ 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
	)

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

Parameters

[in]	lev	Current level
[in,out]	mfs	MultiFabs to populate
[in]	eddyDiffs	Diffusion coefficients from turbulence model

 {
     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) {
         custom_flux flux_comp(specified_rho_surf);
         compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
                                  xheat_flux, yheat_flux, zheat_flux,
                                  xqv_flux, yqv_flux, zqv_flux,
                                  z_phys, flux_comp);
     } 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
	)

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

Parameters

[in]	lev	Current level
[in,out]	mfs	MultiFabs to populate
[in]	eddyDiffs	Diffusion coefficients from turbulence model

 {
     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)`
	)

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

◆ lmask_min_reduce()

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

inline

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

◆ 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

   {
       // 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.query("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);
  
           // 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;
  
       // 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;
               m_has_ocean_lsm_tsurf = (lc_name == "t_surf");
           }
       }
       int n_valid_lsm_flux = 0;
       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; }
       }
       AMREX_ALWAYS_ASSERT((n_valid_lsm_flux==0 || n_valid_lsm_flux>=4));
       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.query("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) {
           // Valid SST, TSK or LSM 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_ocean_lsm_tsurf) {
               amrex::Print() << "(OceanSurf: t_surf)" << std::endl;
           } 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)" << std::endl;
               } else {
                   amrex::Print() << "sea: SST)" << std::endl;
                   AMREX_ALWAYS_ASSERT(m_sst_lev[lev][0]);
               }
           }
       }
   }

Here is the call graph for this function:

◆ read_custom_roughness()

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

 {
     // 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_q_surf()

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

inline

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

◆ set_t_surf()

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

inline

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

◆ 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`
	)

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

Parameters

[in]	lev	Current level
[in]	max_iters	maximum iterations to use

 {
     // Update with SST/TSK data if we have a valid pointer
     if (!m_has_ocean_lsm_tsurf &&
         !m_sst_lev[lev].empty() && m_sst_lev[lev][0]) {
         fill_tsurf_with_sst_and_tsk(lev, elapsed_time_since_start_low);
     }
  
     // Apply heating rate if needed
     if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE &&
         !m_has_ocean_lsm_tsurf) {
         update_surf_temp(elapsed_time_since_start_low);
     }
  
     // 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 (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

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

 {
     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_sst_ptr()

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

inline

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

◆ update_surf_temp()

void SurfaceLayer::update_surf_temp ( const double & time )

inline

   {
       if (m_has_ocean_lsm_tsurf) {
           return;
       }
       if (surf_heating_rate != 0) {
           int nlevs = static_cast<int>(m_geom.size());
           for (int lev = 0; lev < nlevs; lev++) {
               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

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

Referenced by update_surf_temp().

◆ 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_has_ocean_lsm_tsurf

bool SurfaceLayer::m_has_ocean_lsm_tsurf = false

private

Referenced by make_SurfaceLayer_at_level(), and update_surf_temp().

◆ 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

◆ 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_zref(), make_SurfaceLayer_at_level(), 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_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}

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

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

◆ specified_rho_surf

bool SurfaceLayer::specified_rho_surf {false}

private

◆ surf_heating_rate

amrex::Real SurfaceLayer::surf_heating_rate {0}

private

Referenced by update_surf_temp().

◆ surf_moist

amrex::Real SurfaceLayer::surf_moist

private

◆ surf_moist_flux

amrex::Real SurfaceLayer::surf_moist_flux {0}

private

◆ surf_temp

amrex::Real SurfaceLayer::surf_temp

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

Public Types

Public Member Functions

Public Attributes

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ FluxCalcType

◆ MoistCalcType

◆ PBLHeightCalcType

◆ RoughCalcType

◆ ThetaCalcType

Constructor & Destructor Documentation

◆ SurfaceLayer()

Member Function Documentation

◆ compute_fluxes() [1/2]

◆ compute_fluxes() [2/2]

◆ compute_pblh() [1/2]

◆ compute_pblh() [2/2]

◆ compute_sfc_params_from_lsm_fluxes()

◆ compute_SurfaceLayer_bcs() [1/2]

◆ compute_SurfaceLayer_bcs() [2/2]

◆ compute_SurfaceLayer_bcs_EB() [1/2]

◆ compute_SurfaceLayer_bcs_EB() [2/2]

◆ fill_qsurf_with_qsat()

◆ fill_tsurf_with_sst_and_tsk()

◆ get_lmask()

◆ get_lsm_tsurf()

◆ get_mac_avg()

◆ get_olen()

◆ get_pblh()

◆ get_q_star()

◆ get_q_surf()

◆ get_surface_diagnostic_source()

◆ get_t_star()

◆ get_t_surf()

◆ get_u_star()

◆ get_w_star()

◆ get_z0()

◆ get_zref()

◆ have_variable_sea_roughness()

◆ impose_SurfaceLayer_bcs()

◆ impose_SurfaceLayer_bcs_EB()

◆ init_tke_from_ustar()

◆ lmask_min_reduce()

◆ make_SurfaceLayer_at_level()

◆ read_custom_roughness()

◆ set_q_surf()

◆ set_t_surf()

◆ update_fluxes()

◆ update_mac_ptrs()

◆ update_pblh()

◆ update_sst_ptr()

◆ update_surf_temp()

◆ update_tsk_ptr()

Member Data Documentation

◆ cnk_a

◆ cnk_visc

◆ custom_qstar

◆ custom_rhosurf

◆ custom_tstar

◆ custom_ustar

◆ default_land_surf_moist

◆ default_land_surf_temp

◆ depth

◆ flux_type

◆ inv_Cmu2

◆ m_Cd

◆ m_Ch

◆ m_Cq

◆ m_eb_vec

◆ m_eddyDiffs_lev

◆ m_final_low_time

◆ m_geom

◆ m_has_lsm_fluxes

◆ m_has_lsm_tsurf

◆ m_has_ocean_lsm_tsurf

◆ m_Hwave_lev

◆ m_ignore_sst

◆ m_include_wstar

◆ m_lmask_lev