#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:

[legend]

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

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

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

enum class	PBLHeightCalcType { None , MYNN25 , YSU }

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 TerrainType &a_terrain_type, amrex::Real start_bdy_time=0.0, amrex::Real bdy_time_interval=0.0)

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< amrex::MultiFab * > lsm_flux, 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 amrex::Real &time, int max_iters=500)

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

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)

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)

void	fill_tsurf_with_sst_and_tsk (const int &lev, const amrex::Real &time)

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 int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp)

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 int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp)

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

void	update_surf_temp (const amrex::Real &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)

amrex::Real	get_zref ()

amrex::FArrayBox *	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)

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 PBLHeightEstimator >
void	compute_pblh (const int &lev, Vector< Vector< MultiFab >> &vars, MultiFab *z_phys_cc, const PBLHeightEstimator &est, int RhoQv_comp, int RhoQc_comp, int RhoQr_comp)

Public Attributes
FluxCalcType	flux_type {FluxCalcType::MOENG}

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

amrex::Real	m_start_bdy_time

amrex::Real	m_bdy_time_interval

bool	m_include_wstar = false

amrex::Real	z0_const {0.1}

amrex::Real	default_land_surf_temp = 300.

amrex::Real	surf_temp

amrex::Real	surf_heating_rate {0}

amrex::Real	surf_temp_flux {0}

amrex::Real	custom_ustar {0}

amrex::Real	custom_tstar {0}

amrex::Real	custom_qstar {0}

amrex::Real	cnk_a {0.0185}

bool	cnk_visc {false}

amrex::Real	depth {30.0}

amrex::Vector< amrex::FArrayBox >	z_0

bool	m_var_z0 {false}

bool	use_moisture

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< 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< amrex::MultiFab * >	m_Hwave_lev

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

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

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–489. https://doi.org/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/10.1142/2975.

Member Enumeration Documentation

◆ FluxCalcType

enum SurfaceLayer::FluxCalcType

strong

Enumerator
MOENG	Moeng functional form.
DONELAN	Donelan functional form.
CUSTOM	Custom constant flux functional form.
ROTATE	Terrain rotation flux functional form.

                              {
         MOENG = 0,      ///< Moeng functional form
         DONELAN,        ///< Donelan functional form
         CUSTOM,         ///< Custom constant flux functional form
         ROTATE          ///< Terrain rotation flux functional form
     };

◆ PBLHeightCalcType

enum SurfaceLayer::PBLHeightCalcType

strong

Enumerator
None
MYNN25
YSU

                                   {
         None, MYNN25, YSU
     };

◆ 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 TerrainType &	a_terrain_type,
		amrex::Real	start_bdy_time = `0.0`,
		amrex::Real	bdy_time_interval = `0.0`
	)

inlineexplicit

     : m_geom(geom),
       m_rotate(use_rot_surface_flux),
       m_start_bdy_time(start_bdy_time),
       m_bdy_time_interval(bdy_time_interval),
       m_ma(geom,(z_phys_nd[0]!=nullptr),a_pp_prefix,a_terrain_type)
     {
         // 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 == "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 {
             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; }
  
         // Custom type user must specify the fluxes
         if (flux_type == FluxCalcType::CUSTOM) {
             theta_type = ThetaCalcType::HEAT_FLUX;
             pp.query("most.ustar", custom_ustar);
             pp.query("most.tstar", custom_tstar);
             pp.query("most.qstar", custom_qstar);
             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;
  
         // Specify surface temperature or surface flux
         } else {
             if (erf_st) {
                 theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
                 pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
                 surf_heating_rate = surf_heating_rate / 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;
                 }
             }
         }
  
         // 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!");
         }
     } // constructor

Here is the call graph for this function:

Member Function Documentation

◆ compute_fluxes()

template<typename FluxIter >

void SurfaceLayer::compute_fluxes	(	const int &	lev,
		const int &	max_iters,
		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
  
     for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         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 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 z0_arr  = z_0[lev].array();
  
         // 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> {};
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             if (( is_land && lmask_arr(i,j,k) == 1) ||
                 (!is_land && lmask_arr(i,j,k) == 0))
             {
                 most_flux.iterate_flux(i, j, k, max_iters,
                                        z0_arr, umm_arr, tm_arr, tvm_arr, qvm_arr,
                                        u_star_arr, w_star_arr,  // to be updated
                                        t_star_arr, q_star_arr,  // to be updated
                                        t_surf_arr, olen_arr,    // to be updated
                                        pblh_arr, Hwave_arr, Lwave_arr, eta_arr);
             }
         });
     }
 }

◆ 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 int	RhoQv_comp,
		const int	RhoQc_comp,
		const int	RhoQr_comp
	)

◆ 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,
		int	RhoQv_comp,
		int	RhoQc_comp,
		int	RhoQr_comp
	)

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

◆ 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);
  
         // Get LSM fluxes
         auto lmask_arr    = (m_lmask_lev[lev][0])    ? m_lmask_lev[lev][0]->array(mfi) :
                                                        Array4<int> {};
         auto lsm_flux_arr = (m_lsm_flux_lev[lev][0]) ? m_lsm_flux_lev[lev][0]->array(mfi) :
                                                        Array4<Real> {};
  
         // 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)
         {
             Real 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 (rotate) {
                 rotate_scalar_flux(i, j, k, Tflux, dxInv, zphys_arr,
                                    hfx1_arr, hfx2_arr, hfx3_arr);
             } else {
                 hfx3_arr(i,j,klo) = Tflux;
                 int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
                 if (is_land && lsm_flux_arr) {
                     lsm_flux_arr(i,j,k) = Tflux;
                 }
             }
         });
  
         // Rho*Qv flux
         //============================================================================
         // TODO: Generalize MOST q flux
         if ((flux_type == FluxCalcType::CUSTOM) && use_moisture) {
             ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 Real 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, t_surf_arr);
  
                 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
  
         // Rho*u flux
         //============================================================================
         Box bxx = surroundingNodes(bx,0);
         ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             Real stressx = flux_comp.compute_u_flux(i, j, k,
                                                     cons_arr, velx_arr, vely_arr,
                                                     umm_arr, um_arr, u_star_arr);
  
             if (rotate) {
                 rotate_stress_tensor(i, j, k, stressx, 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);
             } else {
                 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)
         {
             Real stressy = flux_comp.compute_v_flux(i, j, k,
                                                     cons_arr, velx_arr, vely_arr,
                                                     umm_arr, vm_arr, u_star_arr);
  
             // NOTE: One stress rotation for ALL the stress components
             if (!rotate) {
                 t23_arr(i,j,k) = stressy;
                 if (t32_arr) { t32_arr(i,j,k) = stressy; }
             }
         });
     } // mfiter
 }

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 amrex::Real &	time
	)

 {
     int n_times_in_sst = m_sst_lev[lev].size();
  
     int n_time_lo, n_time_hi;
     Real alpha;
  
     if (n_times_in_sst > 1) {
         // Time interpolation
         Real dT = m_bdy_time_interval;
         Real time_since_start = time - m_start_bdy_time;
         int n_time = static_cast<int>( time_since_start /  dT);
         n_time_lo = n_time;
         n_time_hi = n_time+1;
         alpha = (time_since_start - n_time * dT) / dT;
         AMREX_ALWAYS_ASSERT( (n_time >= 0) && (n_time < (m_sst_lev[lev].size()-1)));
     } else {
         n_time_lo = 0;
         n_time_hi = 0;
         alpha     = 1.0;
     }
     AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
  
     Real oma   = 1.0 - 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]);
  
     // Populate t_surf
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         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) {
                     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());
 }

◆ get_lmask()

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

inline

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

◆ get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf ( const int & lev )

 {
     for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
  
         // TODO: 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. Is there
         //       a better approach?
         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][0]->const_array(mfi);
  
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
             if (is_land) {
                 int li = amrex::min(amrex::max(i, i_lo), i_hi);
                 int lj = amrex::min(amrex::max(j, j_lo), j_hi);
                 t_surf_arr(i,j,k) = lsm_arr(li,lj,k);
             }
         });
     }
 }

◆ get_mac_avg()

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

inline

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

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

◆ get_pblh()

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

inline

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

◆ get_q_star()

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

inline

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

◆ get_t_star()

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

inline

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

◆ get_t_surf()

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

inline

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

◆ get_u_star()

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

inline

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

◆ get_w_star()

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

inline

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

◆ get_z0()

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

inline

474 { return &z_0[lev]; }

◆ get_zref()

amrex::Real SurfaceLayer::get_zref ( )

inline

471 { return m_ma.get_zref(); }

MOSTAverage::get_zref

amrex::Real get_zref() const

Definition: ERF_MOSTAverage.H:104

Here is the call graph for this function:

◆ have_variable_sea_roughness()

bool SurfaceLayer::have_variable_sea_roughness ( )

inline

476 { return m_var_z0; }

SurfaceLayer::m_var_z0

bool m_var_z0

Definition: ERF_SurfaceLayer.H:552

◆ 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 {
         custom_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);
     }
 }

◆ 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< amrex::MultiFab * >	lsm_flux,
		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 ( (lev == 0) || (lev > nlevs-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);
             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 = 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 = 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 nvar = lsm_data.size();
         m_lsm_data_lev[lev].resize(nvar);
         m_lsm_flux_lev[lev].resize(nvar);
         for (int n(0); n<nvar; ++n) {
             m_lsm_data_lev[lev][n] = lsm_data[n];
             m_lsm_flux_lev[lev][n] = lsm_flux[n];
         }
  
         // 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) ) {
             read_z0 = pp.query("most.roughness_file_name", fname);
         }
  
         // Attributes for MFs and FABs
         //--------------------------------------------------------
         // Create a 2D ba, dm, & ghost cells
         amrex::BoxArray ba  = mf.boxArray();
         amrex::BoxList bl2d = ba.boxList();
         for (auto& b : bl2d) {
             b.setRange(2,0);
         }
         amrex::BoxArray ba2d(std::move(bl2d));
         const amrex::DistributionMapping& dm = mf.DistributionMap();
         const int ncomp   = 1;
         amrex::IntVect ng = mf.nGrowVect(); ng[2]=0;
  
         // Z0 heights FAB
         //--------------------------------------------------------
         amrex::Arena* Arena_Used = amrex::The_Arena();
 #ifdef AMREX_USE_GPU
         Arena_Used = amrex::The_Pinned_Arena();
 #endif
         amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng);
         bx.setSmall(2,0);
         bx.setBig(2,0);
         z_0[lev].resize(bx,1,Arena_Used);
         z_0[lev].setVal<amrex::RunOn::Device>(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>(ba2d,dm,ncomp,ng);
         u_star[lev]->setVal(1.E34);
  
         w_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
         w_star[lev]->setVal(1.E34);
  
         t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
         t_star[lev]->setVal(0.0); // default to neutral
  
         q_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
         q_star[lev]->setVal(0.0); // default to dry
  
         olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
         olen[lev]->setVal(1.E34);
  
         pblh[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
         pblh[lev]->setVal(1.E34);
  
         t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
         t_surf[lev]->setVal(default_land_surf_temp);
  
         // TODO: Do we want lsm_data to have theta at 0 index always?
         //       Do we want an external enum struct for indexing?
         if (m_sst_lev[lev][0] || m_tsk_lev[lev][0] || m_lsm_data_lev[lev][0]) {
             // 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;
         }
     }

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 are on the coarsest level
     if (lev==0) {
         // Only the ioproc reads the file and broadcasts
         if (ParallelDescriptor::IOProcessor()) {
             Print()<<"Reading MOST roughness file: "<< fname << std::endl;
             std::ifstream file(fname);
             Gpu::HostVector<Real> m_x,m_y,m_z0;
             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();
  
             // Copy data to the GPU
             int nnode = m_x.size();
             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();
  
             // Populate z_phys data
             Real tol = 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 klo = 0;
             int ihi = m_geom[lev].Domain().bigEnd(0);
             int jhi = m_geom[lev].Domain().bigEnd(1);
  
             // Grown box with no z range
             Box xybx = z_0[lev].box();
             xybx.setRange(2,0);
  
             Array4<Real> const& z0_arr = z_0[lev].array();
             ParallelFor(xybx, [=] 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(std::pow(x-xp[inode],2)+std::pow(y-yp[inode],2)) < tol) {
                     z0_arr(i,j,klo) = z0p[inode];
                 } else {
                     // Unexpected list order, do brute force search
                     Real z0loc = 0.0;
                     bool found = false;
                     for (int n=0; n<nnode; ++n) {
                         Real delta=std::sqrt(std::pow(x-xp[n],2)+std::pow(y-yp[n],2));
                         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;
                 }
             });
         } // Is ioproc
  
         int ioproc = ParallelDescriptor::IOProcessorNumber();
         ParallelDescriptor::Barrier();
         ParallelDescriptor::Bcast(z_0[lev].dataPtr(),z_0[lev].box().numPts(),ioproc);
     } else {
         // 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
         Interpolater* interp = &cell_cons_interp;
         interp->interp(z_0[0]  , 0,
                        z_0[lev], 0,
                        1, z_0[lev].box(),
                        ratio, m_geom[0], m_geom[lev],
                        bcr, 0, 0, RunOn::Gpu);
     }
 }

Referenced by make_SurfaceLayer_at_level().

Here is the caller graph for this function:

◆ update_fluxes()

void SurfaceLayer::update_fluxes	(	const int &	lev,
		const amrex::Real &	time,
		int	max_iters = `500`
	)

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

Parameters

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

 {
     // Update SST data if we have a valid pointer
     if (m_sst_lev[lev][0]) fill_tsurf_with_sst_and_tsk(lev, time);
  
     // TODO: we want 0 index to always be theta?
     // Update land surface temp if we have a valid pointer
     if (m_lsm_data_lev[lev][0]) 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);
  
     // ***************************************************************
     // 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;
         if (theta_type == ThetaCalcType::HEAT_FLUX) {
             if (rough_type_land == RoughCalcType::CONSTANT) {
                 surface_flux most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else {
                 amrex::Abort("Unknown value for rough_type_land");
             }
         } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
             update_surf_temp(time);
             if (rough_type_land == RoughCalcType::CONSTANT) {
                 surface_temp most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else {
                 amrex::Abort("Unknown value for rough_type_land");
             }
         } else if (theta_type == ThetaCalcType::ADIABATIC) {
             if (rough_type_land == RoughCalcType::CONSTANT) {
                 adiabatic most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else {
                 amrex::Abort("Unknown value for rough_type_land");
             }
         } else {
             amrex::Abort("Unknown value for theta_type");
         }
     } // MOENG -- LAND
  
     // ***************************************************************
     // 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
     // ***************************************************************
     if (flux_type == FluxCalcType::MOENG ||
         flux_type == FluxCalcType::ROTATE) {
         bool is_land = false;
         if (theta_type == ThetaCalcType::HEAT_FLUX) {
             if (rough_type_sea == RoughCalcType::CHARNOCK) {
                 surface_flux_charnock most_flux(m_ma.get_zref(), surf_temp_flux, cnk_a, cnk_visc);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
                 surface_flux_mod_charnock most_flux(m_ma.get_zref(), surf_temp_flux, depth);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::DONELAN) {
                 surface_flux_donelan most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
                 surface_flux_wave_coupled most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else {
                 amrex::Abort("Unknown value for rough_type_sea");
             }
  
         } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
             update_surf_temp(time);
             if (rough_type_sea == RoughCalcType::CHARNOCK) {
                 surface_temp_charnock most_flux(m_ma.get_zref(), surf_temp_flux, cnk_a, cnk_visc);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
                 surface_temp_mod_charnock most_flux(m_ma.get_zref(), surf_temp_flux, depth);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::DONELAN) {
                 surface_temp_donelan most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
                 surface_temp_wave_coupled most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else {
                 amrex::Abort("Unknown value for rough_type_sea");
             }
  
         } else if (theta_type == ThetaCalcType::ADIABATIC) {
             if (rough_type_sea == RoughCalcType::CHARNOCK) {
                 adiabatic_charnock most_flux(m_ma.get_zref(), surf_temp_flux, cnk_a, cnk_visc);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
                 adiabatic_mod_charnock most_flux(m_ma.get_zref(), surf_temp_flux, depth);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::DONELAN) {
                 adiabatic_donelan most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, most_flux, is_land);
             } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
                 adiabatic_wave_coupled most_flux(m_ma.get_zref(), surf_temp_flux);
                 compute_fluxes(lev, max_iters, 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) {
         u_star[lev]->setVal(custom_ustar);
         t_star[lev]->setVal(custom_tstar);
         q_star[lev]->setVal(custom_qstar);
     }
 }

◆ 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

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

MOSTAverage::update_field_ptrs

void update_field_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)

Definition: ERF_MOSTAverage.cpp:246

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 int	RhoQv_comp,
		const int	RhoQc_comp,
		const int	RhoQr_comp
	)

 {
     if (pblh_type == PBLHeightCalcType::MYNN25) {
         MYNNPBLH estimator;
         compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQc_comp, RhoQr_comp);
     } else if (pblh_type == PBLHeightCalcType::YSU) {
         amrex::Error("YSU PBLH calc not implemented yet");
     }
 }

◆ update_surf_temp()

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

inline

     {
         if (surf_heating_rate != 0) {
             int nlevs = m_geom.size();
             for (int lev = 0; lev < nlevs; lev++)
             {
                t_surf[lev]->setVal(surf_temp + surf_heating_rate*time);
                amrex::Print() << "Surface temp at t=" << time
                    << ": "
                    << surf_temp + surf_heating_rate*time
                    << std::endl;
             }
         }
     }

Member Data Documentation

◆ cnk_a

amrex::Real SurfaceLayer::cnk_a {0.0185}

private

Referenced by SurfaceLayer().

◆ cnk_visc

bool SurfaceLayer::cnk_visc {false}

private

Referenced by SurfaceLayer().

◆ custom_qstar

amrex::Real SurfaceLayer::custom_qstar {0}

private

Referenced by SurfaceLayer().

◆ custom_tstar

amrex::Real SurfaceLayer::custom_tstar {0}

private

Referenced by SurfaceLayer().

◆ custom_ustar

amrex::Real SurfaceLayer::custom_ustar {0}

private

Referenced by SurfaceLayer().

◆ default_land_surf_temp

amrex::Real SurfaceLayer::default_land_surf_temp = 300.

private

Referenced by make_SurfaceLayer_at_level(), and SurfaceLayer().

◆ depth

amrex::Real SurfaceLayer::depth {30.0}

private

Referenced by SurfaceLayer().

◆ flux_type

FluxCalcType SurfaceLayer::flux_type {FluxCalcType::MOENG}

Referenced by make_SurfaceLayer_at_level(), and SurfaceLayer().

◆ m_bdy_time_interval

amrex::Real SurfaceLayer::m_bdy_time_interval

private

◆ m_eddyDiffs_lev

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

private

Referenced by make_SurfaceLayer_at_level().

◆ m_geom

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

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_include_wstar

bool SurfaceLayer::m_include_wstar = false

private

Referenced by SurfaceLayer().

◆ m_lmask_lev

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

private

Referenced by get_lmask(), and make_SurfaceLayer_at_level().

◆ m_lsm_data_lev

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

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

◆ m_start_bdy_time

amrex::Real SurfaceLayer::m_start_bdy_time

private

◆ m_tsk_lev

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

private

Referenced by make_SurfaceLayer_at_level().

◆ m_var_z0

bool SurfaceLayer::m_var_z0 {false}

private

Referenced by have_variable_sea_roughness(), and make_SurfaceLayer_at_level().

◆ olen

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

private

Referenced by get_olen(), and make_SurfaceLayer_at_level().

◆ pblh

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

private

Referenced by get_pblh(), and make_SurfaceLayer_at_level().

◆ pblh_type

PBLHeightCalcType SurfaceLayer::pblh_type {PBLHeightCalcType::None}

Referenced by SurfaceLayer().

◆ q_star

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

private

Referenced by get_q_star(), and make_SurfaceLayer_at_level().

◆ rough_type_land

RoughCalcType SurfaceLayer::rough_type_land {RoughCalcType::CONSTANT}

Referenced by SurfaceLayer().

◆ rough_type_sea

RoughCalcType SurfaceLayer::rough_type_sea {RoughCalcType::CHARNOCK}

Referenced by make_SurfaceLayer_at_level(), and SurfaceLayer().

◆ surf_heating_rate

amrex::Real SurfaceLayer::surf_heating_rate {0}

private

Referenced by SurfaceLayer(), and update_surf_temp().

◆ surf_temp

amrex::Real SurfaceLayer::surf_temp

private

Referenced by SurfaceLayer(), and update_surf_temp().

◆ surf_temp_flux

amrex::Real SurfaceLayer::surf_temp_flux {0}

private

Referenced by SurfaceLayer().

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

◆ theta_type

ThetaCalcType SurfaceLayer::theta_type {ThetaCalcType::ADIABATIC}

Referenced by make_SurfaceLayer_at_level(), and SurfaceLayer().

◆ u_star

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

private

Referenced by get_u_star(), and make_SurfaceLayer_at_level().

◆ use_moisture

bool SurfaceLayer::use_moisture

private

Referenced by SurfaceLayer().

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

private

Referenced by make_SurfaceLayer_at_level(), and SurfaceLayer().

◆ z_0

amrex::Vector<amrex::FArrayBox> 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

◆ PBLHeightCalcType

◆ RoughCalcType

◆ ThetaCalcType

Constructor & Destructor Documentation

◆ SurfaceLayer()

Member Function Documentation

◆ compute_fluxes()

◆ compute_pblh() [1/2]

◆ compute_pblh() [2/2]

◆ compute_SurfaceLayer_bcs() [1/2]

◆ compute_SurfaceLayer_bcs() [2/2]

◆ fill_tsurf_with_sst_and_tsk()

◆ get_lmask()

◆ get_lsm_tsurf()

◆ get_mac_avg()

◆ get_olen()

◆ get_pblh()

◆ get_q_star()

◆ get_t_star()

◆ get_t_surf()

◆ get_u_star()

◆ get_w_star()

◆ get_z0()

◆ get_zref()

◆ have_variable_sea_roughness()

◆ impose_SurfaceLayer_bcs()

◆ lmask_min_reduce()

◆ make_SurfaceLayer_at_level()

◆ read_custom_roughness()

◆ update_fluxes()

◆ update_mac_ptrs()

◆ update_pblh()

◆ update_surf_temp()

Member Data Documentation

◆ cnk_a

◆ cnk_visc

◆ custom_qstar

◆ custom_tstar

◆ custom_ustar

◆ default_land_surf_temp

◆ depth

◆ flux_type

◆ m_bdy_time_interval

◆ m_eddyDiffs_lev

◆ m_geom

◆ m_Hwave_lev

◆ m_include_wstar

◆ m_lmask_lev

◆ m_lsm_data_lev

◆ m_lsm_flux_lev

◆ m_Lwave_lev

◆ m_ma

◆ m_rotate

◆ m_sst_lev

◆ m_start_bdy_time

◆ m_tsk_lev

◆ m_var_z0

◆ olen

◆ pblh

◆ pblh_type

◆ q_star

◆ rough_type_land

◆ rough_type_sea

◆ surf_heating_rate

◆ surf_temp

◆ surf_temp_flux

◆ t_star

◆ t_surf

◆ theta_type

◆ u_star

◆ use_moisture

◆ w_star

◆ z0_const