ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Types | Public Member Functions | Public Attributes | Private Attributes | List of all members
SurfaceLayer Class Reference

#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:
Collaboration graph
[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 
978  {
979  MOENG = 0, ///< Moeng functional form
980  DONELAN, ///< Donelan functional form
981  CUSTOM, ///< Custom constant flux functional form
982  BULK_COEFF, ///< Bulk transfer coefficient functional form
983  ROTATE, ///< Terrain rotation flux functional form
984  RICO
985  };
surface_diagnostics::SurfaceDiagnosticSource::RICO
@ RICO

◆ MoistCalcType

enum SurfaceLayer::MoistCalcType
strong
Enumerator
ADIABATIC 
MOISTURE_FLUX 

Qv-flux specified.

SURFACE_MOISTURE 

Surface Qv specified.

993  {
994  ADIABATIC = 0,
995  MOISTURE_FLUX, ///< Qv-flux specified
996  SURFACE_MOISTURE ///< Surface Qv specified
997  };

◆ PBLHeightCalcType

enum SurfaceLayer::PBLHeightCalcType
strong
Enumerator
None 
MYNN25 
YSU 
MRF 
YSUNew 
1007 { None, MYNN25, YSU, MRF, YSUNew };
erf_wall_thermodynamics::TransferModel::None
@ None

◆ RoughCalcType

enum SurfaceLayer::RoughCalcType
strong
Enumerator
CONSTANT 

Constant z0.

CHARNOCK 
MODIFIED_CHARNOCK 
DONELAN 
WAVE_COUPLED 
999  {
1000  CONSTANT = 0, ///< Constant z0
1001  CHARNOCK,
1002  MODIFIED_CHARNOCK,
1003  DONELAN,
1004  WAVE_COUPLED
1005  };

◆ ThetaCalcType

enum SurfaceLayer::ThetaCalcType
strong
Enumerator
ADIABATIC 
HEAT_FLUX 

Heat-flux specified.

SURFACE_TEMPERATURE 

Surface temperature specified.

987  {
988  ADIABATIC = 0,
989  HEAT_FLUX, ///< Heat-flux specified
990  SURFACE_TEMPERATURE ///< Surface temperature specified
991  };

Constructor & Destructor Documentation

◆ SurfaceLayer()

SurfaceLayer::SurfaceLayer ( const amrex::Vector< amrex::Geometry > &  geom,
bool &  use_rot_surface_flux,
std::string  a_pp_prefix,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  Qv_prim,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  z_phys_nd,
const MeshType &  a_mesh_type,
const TerrainType &  a_terrain_type,
const TurbChoice &  a_turb_choice,
double  start_low_time,
double  final_low_time,
double  low_time_interval = 0.0,
const amrex::Vector< const eb_ * > &  eb_vec = {} 
)
inlineexplicit

Construct the surface-layer interface.

Parameters
[in]geomgeometry for all AMR levels
[in,out]use_rot_surface_fluxwhether to use rotated surface fluxes
[in]a_pp_prefixParmParse prefix used by MOST averages
[in]Qv_primprimitive water-vapor fields by level
[in]z_phys_ndnodal physical-height fields by level
[in]a_mesh_typemesh type
[in]a_terrain_typeterrain representation
[in]a_turb_choiceturbulence-model options
[in]start_low_timefirst available low-boundary-data time
[in]final_low_timefinal available low-boundary-data time
[in]low_time_intervallow-boundary-data time interval
[in]eb_vecoptional embedded-boundary geometry data
64  {})
65  : m_geom(geom),
66  m_rotate(use_rot_surface_flux),
67  m_start_low_time(start_low_time),
68  m_final_low_time(final_low_time),
69  m_low_time_interval(low_time_interval),
70  m_eb_vec(eb_vec),
71  m_terrain_type(a_terrain_type),
72  m_ma(geom, (z_phys_nd[0] != nullptr), a_pp_prefix, a_mesh_type, a_terrain_type, eb_vec)
73  {
74  // We have a moisture model if Qv_prim is a valid pointer
75  use_moisture = (Qv_prim[0].get());
76 
77  // Get roughness
78  amrex::ParmParse pp("erf");
79  pp.query("most.z0", z0_const);
80 
81  // Specify how to compute the flux
82  if (use_rot_surface_flux) {
83  flux_type = FluxCalcType::ROTATE;
84  } else {
85  std::string flux_string_in;
86  std::string flux_string{"moeng"};
87  auto read_flux = pp.query("surface_layer.flux_type", flux_string_in);
88  if (read_flux) {
89  flux_string = amrex::toLower(flux_string_in);
90  }
91  if (flux_string == "donelan") {
92  flux_type = FluxCalcType::DONELAN;
93  } else if (flux_string == "moeng") {
94  flux_type = FluxCalcType::MOENG;
95  } else if (flux_string == "rico") {
96  flux_type = FluxCalcType::RICO;
97  } else if (flux_string == "bulk_coeff") {
98  flux_type = FluxCalcType::BULK_COEFF;
99  } else if (flux_string == "custom") {
100  flux_type = FluxCalcType::CUSTOM;
101  } else {
102  amrex::Abort("Undefined MOST flux type!");
103  }
104  }
105 
106  // Include w* to handle free convection (Beljaars 1995, QJRMS)
107  pp.query("most.include_wstar", m_include_wstar);
108 
109  std::string pblh_string_in;
110  std::string pblh_string{"none"};
111  auto read_pblh = pp.query("most.pblh_calc", pblh_string_in);
112  if (read_pblh) {
113  pblh_string = amrex::toLower(pblh_string_in);
114  }
115  if (pblh_string == "none") {
116  pblh_type = PBLHeightCalcType::None;
117  } else if (pblh_string == "mynn25") {
118  pblh_type = PBLHeightCalcType::MYNN25;
119  } else if (pblh_string == "mynnedmf") {
120  pblh_type = PBLHeightCalcType::MYNN25;
121  } else if (pblh_string == "ysu") {
122  pblh_type = PBLHeightCalcType::YSU;
123  } else if (pblh_string == "mrf") {
124  pblh_type = PBLHeightCalcType::MRF;
125  } else {
126  amrex::Abort("Undefined PBLH calc type!");
127  }
128 
129  // Get surface temperature
130  auto erf_st = pp.query("most.surf_temp", surf_temp);
131  if (erf_st) { default_land_surf_temp = surf_temp; }
132 
133  // Get surface moisture
134  bool erf_sq = false;
135  if (use_moisture) { erf_sq = pp.query("most.surf_moist", surf_moist); }
136  if (erf_sq) { default_land_surf_moist = surf_moist; }
137 
138  // Custom type user must specify the fluxes
139  if (flux_type == FluxCalcType::CUSTOM ||
140  flux_type == FluxCalcType::RICO) {
141  theta_type = ThetaCalcType::HEAT_FLUX;
142  moist_type = MoistCalcType::MOISTURE_FLUX;
143  pp.get("most.ustar", custom_ustar);
144  pp.get("most.tstar", custom_tstar);
145  pp.get("most.qstar", custom_qstar);
146  pp.query("most.rhosurf", custom_rhosurf);
147  if (custom_qstar != 0) {
148  AMREX_ASSERT_WITH_MESSAGE(use_moisture,
149  "Specified custom MOST qv flux without moisture model!");
150  }
151  amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
152  << custom_ustar << " " << custom_tstar << " "
153  << custom_qstar << std::endl;
154 
155  // Bulk transfer coefficient (must specify coeffs and surface values)
156  } else if (flux_type == FluxCalcType::BULK_COEFF) {
157  pp.get("most.Cd", m_Cd);
158  pp.get("most.Ch", m_Ch);
159  pp.get("most.Cq", m_Cq);
160  pp.get("most.surf_temp", default_land_surf_temp);
161  pp.get("most.surf_moist", default_land_surf_moist);
162  amrex::Print() << "Using specified Cd, Ch, Cq for MOST = "
163  << m_Cd << " " << m_Ch << " "
164  << m_Cq << std::endl;
165 
166  // Specify surface temperature/moisture or surface flux
167  } else {
168  if (erf_st) {
169  theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
170  pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
171 
172  // Modify rate to be in units of K / s rather than K / hr
173  surf_heating_rate /= amrex::Real(3600.0); // [K/s]
174 
175  if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
176  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
177  }
178  } else {
179  pp.query("most.surf_temp_flux", surf_temp_flux);
180 
181  if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
182  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
183  }
184  if (std::abs(surf_temp_flux) >
185  std::numeric_limits<amrex::Real>::epsilon()) {
186  theta_type = ThetaCalcType::HEAT_FLUX;
187  } else {
188  theta_type = ThetaCalcType::ADIABATIC;
189  }
190  }
191 
192  if (erf_sq) {
193  moist_type = MoistCalcType::SURFACE_MOISTURE;
194  } else {
195  pp.query("most.surf_moist_flux", surf_moist_flux);
196  if (std::abs(surf_moist_flux) >
197  std::numeric_limits<amrex::Real>::epsilon()) {
198  moist_type = MoistCalcType::MOISTURE_FLUX;
199  } else {
200  moist_type = MoistCalcType::ADIABATIC;
201  }
202  }
203  }
204 
205  if (flux_type == FluxCalcType::RICO)
206  {
207  pp.query("most.rico.theta_z0", rico_theta_z0);
208  pp.query("most.rico.qsat_z0", rico_qsat_z0);
209  }
210 
211  // Make sure the inputs file doesn't try to use most.roughness_type
212  std::string bogus_input;
213  if (pp.query("most.roughness_type", bogus_input) > 0) {
214  amrex::Abort("most.roughness_type is deprecated; use "
215  "most.roughness_type_land and/or most.roughness_type_sea");
216  }
217 
218  // Specify how to compute the surface flux over land (if there is any)
219  std::string rough_land_string_in;
220  std::string rough_land_string{"constant"};
221  auto read_rough_land =
222  pp.query("most.roughness_type_land", rough_land_string_in);
223  if (read_rough_land) {
224  rough_land_string = amrex::toLower(rough_land_string_in);
225  }
226  if (rough_land_string == "constant") {
227  rough_type_land = RoughCalcType::CONSTANT;
228  } else {
229  amrex::Abort("Undefined MOST roughness type for land!");
230  }
231 
232  // Specify how to compute the surface flux over sea (if there is any)
233  std::string rough_sea_string_in;
234  std::string rough_sea_string{"charnock"};
235  auto read_rough_sea = pp.query("most.roughness_type_sea", rough_sea_string_in);
236  if (read_rough_sea) {
237  rough_sea_string = amrex::toLower(rough_sea_string_in);
238  }
239  if (rough_sea_string == "charnock") {
240  rough_type_sea = RoughCalcType::CHARNOCK;
241  pp.query("most.charnock_constant", cnk_a);
242  pp.query("most.charnock_viscosity", cnk_visc);
243  if (cnk_a > 0) {
244  amrex::Print() << "If there is water, Charnock relation with C_a="
245  << cnk_a << (cnk_visc ? " and viscosity" : "")
246  << " will be used" << std::endl;
247  } else {
248  amrex::Print() << "If there is water, Charnock relation with variable "
249  "Charnock parameter (COARE3.0)"
250  << (cnk_visc ? " and viscosity" : "") << " will be used"
251  << std::endl;
252  }
253  } else if (rough_sea_string == "coare3.0") {
254  rough_type_sea = RoughCalcType::CHARNOCK;
255  amrex::Print() << "If there is water, Charnock relation with variable "
256  "Charnock parameter (COARE3.0)"
257  << (cnk_visc ? " and viscosity" : "") << " will be used"
258  << std::endl;
259  cnk_a = -1;
260  } else if (rough_sea_string == "donelan") {
261  rough_type_sea = RoughCalcType::DONELAN;
262  } else if (rough_sea_string == "modified_charnock") {
263  rough_type_sea = RoughCalcType::MODIFIED_CHARNOCK;
264  pp.query("most.modified_charnock_depth", depth);
265  } else if (rough_sea_string == "wave_coupled") {
266  rough_type_sea = RoughCalcType::WAVE_COUPLED;
267  } else if (rough_sea_string == "constant") {
268  rough_type_sea = RoughCalcType::CONSTANT;
269  } else {
270  amrex::Abort("Undefined MOST roughness type for sea!");
271  }
272 
273  // use skin temperature instead of sea-surface temperature
274  // (wrfinput data may have lower resolution SST data)
275  pp.query("most.ignore_sst", m_ignore_sst);
276 
277  // If we're using the RANS k model, then we need to update the dirichlet
278  // BC based on the instantaneous u* and θ*; the turbulence modeling
279  // choices can vary per level but for now, assume that if specified then
280  // all levels are using the same RANS model.
281  m_update_k_rans = (a_turb_choice.rans_type == RANSType::kEqn &&
282  a_turb_choice.dirichlet_k == true);
283  if (m_update_k_rans) {
284  inv_Cmu2 = one / (a_turb_choice.Cmu0 * a_turb_choice.Cmu0);
285  theta_ref = a_turb_choice.theta_ref;
286  }
287 
288  } // constructor
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
pp
ParmParse pp("prob")
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_ASSERT_WITH_MESSAGE
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")
SurfaceLayer::theta_type
ThetaCalcType theta_type
Definition: ERF_SurfaceLayer.H:1010
SurfaceLayer::m_include_wstar
bool m_include_wstar
Definition: ERF_SurfaceLayer.H:1024
SurfaceLayer::m_rotate
bool m_rotate
Definition: ERF_SurfaceLayer.H:1019
SurfaceLayer::pblh_type
PBLHeightCalcType pblh_type
Definition: ERF_SurfaceLayer.H:1014
SurfaceLayer::m_final_low_time
double m_final_low_time
Definition: ERF_SurfaceLayer.H:1021
SurfaceLayer::use_moisture
bool use_moisture
Definition: ERF_SurfaceLayer.H:1047
SurfaceLayer::m_Cq
amrex::Real m_Cq
Definition: ERF_SurfaceLayer.H:1055
SurfaceLayer::m_eb_vec
amrex::Vector< const eb_ * > m_eb_vec
Definition: ERF_SurfaceLayer.H:1059
SurfaceLayer::rough_type_land
RoughCalcType rough_type_land
Definition: ERF_SurfaceLayer.H:1012
SurfaceLayer::z0_const
amrex::Real z0_const
Definition: ERF_SurfaceLayer.H:1025
SurfaceLayer::cnk_a
amrex::Real cnk_a
Definition: ERF_SurfaceLayer.H:1038
SurfaceLayer::m_Ch
amrex::Real m_Ch
Definition: ERF_SurfaceLayer.H:1054
SurfaceLayer::surf_temp
amrex::Real surf_temp
Definition: ERF_SurfaceLayer.H:1027
SurfaceLayer::m_start_low_time
double m_start_low_time
Definition: ERF_SurfaceLayer.H:1020
SurfaceLayer::rico_qsat_z0
amrex::Real rico_qsat_z0
Definition: ERF_SurfaceLayer.H:1045
SurfaceLayer::m_update_k_rans
bool m_update_k_rans
Definition: ERF_SurfaceLayer.H:1086
SurfaceLayer::surf_moist_flux
amrex::Real surf_moist_flux
Definition: ERF_SurfaceLayer.H:1032
SurfaceLayer::rough_type_sea
RoughCalcType rough_type_sea
Definition: ERF_SurfaceLayer.H:1013
SurfaceLayer::surf_moist
amrex::Real surf_moist
Definition: ERF_SurfaceLayer.H:1031
SurfaceLayer::m_ignore_sst
bool m_ignore_sst
Definition: ERF_SurfaceLayer.H:1057
SurfaceLayer::m_low_time_interval
double m_low_time_interval
Definition: ERF_SurfaceLayer.H:1022
SurfaceLayer::custom_qstar
amrex::Real custom_qstar
Definition: ERF_SurfaceLayer.H:1035
SurfaceLayer::custom_rhosurf
amrex::Real custom_rhosurf
Definition: ERF_SurfaceLayer.H:1036
SurfaceLayer::FluxCalcType::MOENG
@ MOENG
Moeng functional form.
SurfaceLayer::FluxCalcType::BULK_COEFF
@ BULK_COEFF
Bulk transfer coefficient functional form.
SurfaceLayer::FluxCalcType::CUSTOM
@ CUSTOM
Custom constant flux functional form.
SurfaceLayer::FluxCalcType::ROTATE
@ ROTATE
Terrain rotation flux functional form.
SurfaceLayer::FluxCalcType::DONELAN
@ DONELAN
Donelan functional form.
SurfaceLayer::MoistCalcType::SURFACE_MOISTURE
@ SURFACE_MOISTURE
Surface Qv specified.
SurfaceLayer::MoistCalcType::MOISTURE_FLUX
@ MOISTURE_FLUX
Qv-flux specified.
SurfaceLayer::depth
amrex::Real depth
Definition: ERF_SurfaceLayer.H:1040
SurfaceLayer::default_land_surf_moist
amrex::Real default_land_surf_moist
Definition: ERF_SurfaceLayer.H:1030
SurfaceLayer::rico_theta_z0
amrex::Real rico_theta_z0
Definition: ERF_SurfaceLayer.H:1044
SurfaceLayer::surf_temp_flux
amrex::Real surf_temp_flux
Definition: ERF_SurfaceLayer.H:1029
SurfaceLayer::m_geom
amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_SurfaceLayer.H:1018
SurfaceLayer::theta_ref
amrex::Real theta_ref
Definition: ERF_SurfaceLayer.H:1088
SurfaceLayer::custom_tstar
amrex::Real custom_tstar
Definition: ERF_SurfaceLayer.H:1034
SurfaceLayer::cnk_visc
bool cnk_visc
Definition: ERF_SurfaceLayer.H:1039
SurfaceLayer::surf_heating_rate
amrex::Real surf_heating_rate
Definition: ERF_SurfaceLayer.H:1028
SurfaceLayer::flux_type
FluxCalcType flux_type
Definition: ERF_SurfaceLayer.H:1009
SurfaceLayer::moist_type
MoistCalcType moist_type
Definition: ERF_SurfaceLayer.H:1011
SurfaceLayer::inv_Cmu2
amrex::Real inv_Cmu2
Definition: ERF_SurfaceLayer.H:1087
SurfaceLayer::custom_ustar
amrex::Real custom_ustar
Definition: ERF_SurfaceLayer.H:1033
SurfaceLayer::m_Cd
amrex::Real m_Cd
Definition: ERF_SurfaceLayer.H:1053
SurfaceLayer::default_land_surf_temp
amrex::Real default_land_surf_temp
Definition: ERF_SurfaceLayer.H:1026
SurfaceLayer::ThetaCalcType::SURFACE_TEMPERATURE
@ SURFACE_TEMPERATURE
Surface temperature specified.
SurfaceLayer::ThetaCalcType::HEAT_FLUX
@ HEAT_FLUX
Heat-flux specified.
SurfaceLayer::m_terrain_type
TerrainType m_terrain_type
Definition: ERF_SurfaceLayer.H:1060
SurfaceLayer::m_ma
MOSTAverage m_ma
Definition: ERF_SurfaceLayer.H:1061
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
TurbChoice::rans_type
RANSType rans_type
Selected RANS closure.
Definition: ERF_TurbStruct.H:599
TurbChoice::theta_ref
amrex::Real theta_ref
Reference potential temperature for stable stratification.
Definition: ERF_TurbStruct.H:588
TurbChoice::dirichlet_k
bool dirichlet_k
Whether TKE uses Dirichlet boundary treatment.
Definition: ERF_TurbStruct.H:601
TurbChoice::Cmu0
amrex::Real Cmu0
One-equation RANS Cmu0 coefficient.
Definition: ERF_TurbStruct.H:577

Member Function Documentation

◆ compute_fluxes() [1/2]

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

Compute MOST fluxes with a selected flux-iteration functor.

Parameters
[in]levlevel index
[in]max_itersmaximum MOST iteration count
[in,out]cons_inconserved state used by the flux computation
[in]most_fluxflux-iteration functor
[in]is_landwhether the land-surface branch is active

◆ compute_fluxes() [2/2]

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

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

Parameters
[in]levCurrent level
[in]max_itersMaximum iterations to use
[in]cons_inConserved state whose grids define the surface iteration
[in]most_fluxFlux-iteration functor used to compute ustar, tstar, qstar, and related fields
[in]is_landSelects whether land or sea cells are updated
286 {
287  // Pointers to the computed averages
288  const auto *const tm_ptr = m_ma.get_average(lev,2); // potential temperature
289  const auto *const qvm_ptr = m_ma.get_average(lev,3); // water vapor mixing ratio
290  const auto *const tvm_ptr = m_ma.get_average(lev,4); // virtual potential temperature
291  const auto *const umm_ptr = m_ma.get_average(lev,5); // horizontal velocity magnitude
292  const auto *const zref_ptr = m_ma.get_zref(lev); // reference height
293  const bool l_use_eb = (m_terrain_type == TerrainType::EB);
294 
295  const int klo = m_geom[lev].Domain().smallEnd(2);
296  IntVect ng = u_star[lev]->nGrowVect(); ng[2] = 0;
297 
298  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi)
299  {
300  Box gtbx = mfi.tilebox(IntVect(0),ng);
301 
302  if (!l_use_eb && gtbx.smallEnd(2) != klo) { continue; }
303 
304  if (!l_use_eb) { gtbx.makeSlab(2,klo); }
305 
306  auto u_star_arr = u_star[lev]->array(mfi);
307  auto t_star_arr = t_star[lev]->array(mfi);
308  auto q_star_arr = q_star[lev]->array(mfi);
309  auto t_surf_arr = t_surf[lev]->array(mfi);
310  auto q_surf_arr = q_surf[lev]->array(mfi);
311  auto olen_arr = olen[lev]->array(mfi);
312 
313  const auto tm_arr = tm_ptr->array(mfi);
314  const auto tvm_arr = tvm_ptr->array(mfi);
315  const auto qvm_arr = qvm_ptr->array(mfi);
316  const auto umm_arr = umm_ptr->array(mfi);
317  const auto zref_arr = zref_ptr->array(mfi);
318  const auto z0_arr = z_0[lev].array(mfi);
319 
320  // PBL height if we need to calculate wstar for the Beljaars correction
321  // TODO: can/should we apply this in LES mode?
322  const auto w_star_arr = (m_include_wstar) ? w_star[lev].get()->array(mfi) : Array4<Real> {};
323  const auto pblh_arr = (m_include_wstar) ? pblh[lev].get()->array(mfi) : Array4<Real> {};
324 
325  // Wave properties if they exist
326  const auto Hwave_arr = (m_Hwave_lev[lev]) ? m_Hwave_lev[lev]->array(mfi) : Array4<Real> {};
327  const auto Lwave_arr = (m_Lwave_lev[lev]) ? m_Lwave_lev[lev]->array(mfi) : Array4<Real> {};
328  const auto eta_arr = (m_eddyDiffs_lev[lev]) ? m_eddyDiffs_lev[lev]->array(mfi) : Array4<Real> {};
329 
330  // Land mask array if it exists
331  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
332  Array4<int> {};
333 
334  // Get EB flags if needed
335  const auto flag_arr = (l_use_eb) ? m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab()[mfi].const_array() : Array4<const EBCellFlag>{};
336 
337  if (!l_use_eb) {
338  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
339  {
340  if (( is_land && lmask_arr(i,j,0) == 1) ||
341  (!is_land && lmask_arr(i,j,0) == 0))
342  {
343  // NOTE: All 2D MFs so k index is always 0 from ba2d definition
344  most_flux.iterate_flux(i, j, 0, max_iters,
345  zref_arr, // set in most average
346  z0_arr, // updated if(!is_land)
347  umm_arr, tm_arr, tvm_arr, qvm_arr,
348  u_star_arr, // updated
349  w_star_arr, // updated if(m_include_wstar)
350  t_star_arr, q_star_arr, // updated
351  t_surf_arr, q_surf_arr, olen_arr, // updated
352  pblh_arr, // updated if(m_include_wstar)
353  Hwave_arr, Lwave_arr, eta_arr);
354  }
355  });
356  // EB
357  } else {
358  if (std::is_same<FluxIter, adiabatic_eb>::value ||
359  std::is_same<FluxIter, surface_temp_eb>::value ||
360  std::is_same<FluxIter, surface_flux_eb>::value) {
361  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
362  {
363  if (( is_land && lmask_arr(i,j,0) == 1) ||
364  (!is_land && lmask_arr(i,j,0) == 0))
365  {
366  if (flag_arr(i,j,k).isSingleValued()) {
367  most_flux.iterate_flux(i, j, k, max_iters,
368  zref_arr, // set in most average
369  z0_arr, // updated if(!is_land)
370  umm_arr, tm_arr, tvm_arr, qvm_arr,
371  u_star_arr, // updated
372  w_star_arr, // updated if(m_include_wstar)
373  t_star_arr, q_star_arr, // updated
374  t_surf_arr, q_surf_arr, olen_arr, // updated
375  pblh_arr, // updated if(m_include_wstar)
376  Hwave_arr, Lwave_arr, eta_arr);
377  }
378  }
379  });
380  } else {
381  amrex::Abort("FluxIter type not supported for EB");
382  }
383  }
384  }
385 }
value
Real value
Definition: ERF_HurricaneDiagnostics.cpp:30
get
pp get("wavelength", wavelength)
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
MOSTAverage::get_zref
amrex::MultiFab * get_zref(const int &lev) const
Definition: ERF_MOSTAverage.H:242
MOSTAverage::get_average
const amrex::MultiFab * get_average(const int &lev, const int &comp) const
Definition: ERF_MOSTAverage.H:235
SurfaceLayer::m_lmask_lev
amrex::Vector< amrex::Vector< amrex::iMultiFab * > > m_lmask_lev
Definition: ERF_SurfaceLayer.H:1077
SurfaceLayer::t_surf
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_surf
Definition: ERF_SurfaceLayer.H:1068
SurfaceLayer::q_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_star
Definition: ERF_SurfaceLayer.H:1065
SurfaceLayer::m_Lwave_lev
amrex::Vector< amrex::MultiFab * > m_Lwave_lev
Definition: ERF_SurfaceLayer.H:1083
SurfaceLayer::z_0
amrex::Vector< amrex::MultiFab > z_0
Definition: ERF_SurfaceLayer.H:1041
SurfaceLayer::w_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > > w_star
Definition: ERF_SurfaceLayer.H:1063
SurfaceLayer::u_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > > u_star
Definition: ERF_SurfaceLayer.H:1062
SurfaceLayer::q_surf
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_surf
Definition: ERF_SurfaceLayer.H:1069
SurfaceLayer::m_Hwave_lev
amrex::Vector< amrex::MultiFab * > m_Hwave_lev
Definition: ERF_SurfaceLayer.H:1082
SurfaceLayer::t_star
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_star
Definition: ERF_SurfaceLayer.H:1064
SurfaceLayer::m_eddyDiffs_lev
amrex::Vector< amrex::MultiFab * > m_eddyDiffs_lev
Definition: ERF_SurfaceLayer.H:1084
SurfaceLayer::olen
amrex::Vector< std::unique_ptr< amrex::MultiFab > > olen
Definition: ERF_SurfaceLayer.H:1066
SurfaceLayer::pblh
amrex::Vector< std::unique_ptr< amrex::MultiFab > > pblh
Definition: ERF_SurfaceLayer.H:1067
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:49
Here is the call graph for this function:

◆ compute_pblh() [1/2]

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

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

Parameters
[in]levlevel index
[in,out]varsstate variables used by the PBL-height calculation
[in]z_phys_cccell-centered physical-height field
[in]estPBL-height estimator functor
[in]moisture_indiceindices for moisture components

◆ compute_pblh() [2/2]

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

Compute PBL height with the supplied estimator.

Parameters
[in]levCurrent level
[in]varsLevel-indexed state MultiFabs passed to the estimator
[in]z_phys_ccCell-centered physical height used by the estimator
[in]estPBL height estimator functor
[in]moisture_indicesMoisture component indices used by the estimator
1354 {
1355  est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
1356  vars[lev][Vars::cons],m_lmask_lev[lev][0],
1357  moisture_indices);
1358 }
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:176
Here is the call graph for this function:

◆ compute_sfc_params_from_lsm_fluxes()

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

Derive MOST surface parameters from LSM fluxes.

Parameters
[in]levlevel index
[in,out]cons_inconserved state used by the surface-parameter computation

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

Parameters
[in]levCurrent level
[in]cons_inConserved state used to derive density, theta, and moisture at the surface
1066 {
1067  Real eps = std::numeric_limits<amrex::Real>::epsilon();
1068  bool has_moisture = use_moisture;
1069  const int klo = m_geom[lev].Domain().smallEnd(2);
1070  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
1071 
1072  Box vbx = mfi.validbox();
1073  if (vbx.smallEnd(2) != klo) { continue; }
1074  vbx.makeSlab(2,0);
1075 
1076  // Get CC state
1077  const Array4<const Real> cons_arr = cons_in.const_array(mfi);
1078 
1079  // Get SL params
1080  const auto u_star_arr = u_star[lev]->array(mfi);
1081  const auto t_star_arr = t_star[lev]->array(mfi);
1082  const auto q_star_arr = q_star[lev]->array(mfi);
1083  const auto olen_arr = olen[lev]->array(mfi);
1084 
1085  // Get LSM fluxes
1086  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
1087  Array4<int> {};
1088  auto lsm_t_flux_arr = Array4<Real> {};
1089  auto lsm_q_flux_arr = Array4<Real> {};
1090  auto lsm_tau13_arr = Array4<Real> {};
1091  auto lsm_tau23_arr = Array4<Real> {};
1092  // compute_sfc_params_from_lsm_fluxes consumes signed kinematic stress
1093  // components; their vector magnitude determines u_star^2.
1094  for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
1095  if (toLower(m_lsm_flux_name[n]) == "t_flux") { lsm_t_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
1096  if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
1097  if (toLower(m_lsm_flux_name[n]) == "tau13") { lsm_tau13_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
1098  if (toLower(m_lsm_flux_name[n]) == "tau23") { lsm_tau23_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
1099  }
1100 
1101  ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
1102  {
1103  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
1104  // Skip cells the LSM did not have a valid flux (lsm_undefined).
1105  if (is_land && lsm_t_flux_arr && lsm_t_flux_arr(i,j,0) < lsm_undefined) {
1106  Real rho = cons_arr(i,j,klo,Rho_comp);
1107  Real Thd = cons_arr(i,j,klo,RhoTheta_comp) / rho;
1108  Real qv = (has_moisture) ? cons_arr(i,j,klo,RhoQ1_comp) / rho : zero;
1109  Real Thv = Thd * (one + epsv*qv);
1110  Real tau = std::sqrt( lsm_tau13_arr(i,j,0)*lsm_tau13_arr(i,j,0)
1111  + lsm_tau23_arr(i,j,0)*lsm_tau23_arr(i,j,0) );
1112  u_star_arr(i,j,0) = amrex::max(std::sqrt(tau),eps);
1113  if (lsm_t_flux_arr(i,j,0)>=zero) {
1114  t_star_arr(i,j,0) = amrex::min(-lsm_t_flux_arr(i,j,0) / u_star_arr(i,j,0),-eps);
1115  } else {
1116  t_star_arr(i,j,0) = amrex::max(-lsm_t_flux_arr(i,j,0) / u_star_arr(i,j,0),eps);
1117  }
1118  if (lsm_q_flux_arr(i,j,0)>=zero) {
1119  q_star_arr(i,j,0) = amrex::min(-lsm_q_flux_arr(i,j,0) / u_star_arr(i,j,0),-eps);
1120  } else {
1121  q_star_arr(i,j,0) = amrex::max(-lsm_q_flux_arr(i,j,0) / u_star_arr(i,j,0),eps);
1122  }
1123  olen_arr(i,j,0) = ( u_star_arr(i,j,0) * u_star_arr(i,j,0) * Thv ) /
1124  ( KAPPA * CONST_GRAV * t_star_arr(i,j,0) );
1125  }
1126  });
1127  } // mfi
1128 }
epsv
constexpr amrex::Real epsv
Definition: ERF_Constants.H:53
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:63
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
lsm_undefined
constexpr amrex::Real lsm_undefined
Definition: ERF_Constants.H:35
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
rho
rho
Definition: ERF_InitCustomPert_Bubble.H:107
SurfaceLayer::m_lsm_flux_lev
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_flux_lev
Definition: ERF_SurfaceLayer.H:1079
SurfaceLayer::m_lsm_flux_name
amrex::Vector< std::string > m_lsm_flux_name
Definition: ERF_SurfaceLayer.H:1081
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:30
Here is the call graph for this function:

◆ compute_SurfaceLayer_bcs() [1/2]

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

Compute planar-terrain surface-layer flux boundary conditions.

Parameters
[in]levlevel index
[in]mfsstate and velocity fields used by the BC computation
[in,out]Tau_levstress fields to fill
[in,out]xheat_fluxx-face heat flux field
[in,out]yheat_fluxy-face heat flux field
[in,out]zheat_fluxz-face heat flux field
[in,out]xqv_fluxx-face moisture flux field
[in,out]yqv_fluxy-face moisture flux field
[in,out]zqv_fluxz-face moisture flux field
[in]z_physphysical-height field
[in]flux_compflux-computation functor

◆ compute_SurfaceLayer_bcs() [2/2]

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

Function to calculate MOST fluxes for populating ghost cells.

Parameters
[in]levCurrent level
[in]mfsState MultiFabs used to compute the boundary fluxes
[in,out]Tau_levDiffusive stress MultiFabs populated with surface stresses
[in,out]xheat_fluxx-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]yheat_fluxy-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]zheat_fluxz-face heat-flux MultiFab populated with vertical surface heat flux
[in,out]xqv_fluxx-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]yqv_fluxy-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]zqv_fluxz-face moisture-flux MultiFab populated when moisture is enabled
[in]z_physNodal physical height used to rotate terrain-following fluxes
[in]flux_compFlux-calculation functor used to compute scalar and momentum fluxes
518 {
519  bool rotate = m_rotate;
520  const int klo = m_geom[lev].Domain().smallEnd(2);
521  const auto& dxInv = m_geom[lev].InvCellSizeArray();
522  for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
523  {
524  // Get field arrays
525  const auto cons_arr = mfs[Vars::cons]->array(mfi);
526  const auto velx_arr = mfs[Vars::xvel]->array(mfi);
527  const auto vely_arr = mfs[Vars::yvel]->array(mfi);
528  const auto velz_arr = mfs[Vars::zvel]->array(mfi);
529 
530  // Diffusive stress vars
531  auto t13_arr = Tau_lev[TauType::tau13]->array(mfi);
532  auto t31_arr = (Tau_lev[TauType::tau31]) ? Tau_lev[TauType::tau31]->array(mfi) : Array4<Real>{};
533 
534  auto t23_arr = Tau_lev[TauType::tau23]->array(mfi);
535  auto t32_arr = (Tau_lev[TauType::tau32]) ? Tau_lev[TauType::tau32]->array(mfi) : Array4<Real>{};
536 
537  auto hfx3_arr = zheat_flux->array(mfi);
538  auto qfx3_arr = (zqv_flux) ? zqv_flux->array(mfi) : Array4<Real>{};
539 
540  auto olen_arr = olen[lev]->array(mfi);
541 
542  // Rotated stress vars
543  auto t11_arr = (m_rotate) ? Tau_lev[TauType::tau11]->array(mfi) : Array4<Real>{};
544  auto t22_arr = (m_rotate) ? Tau_lev[TauType::tau22]->array(mfi) : Array4<Real>{};
545  auto t33_arr = (m_rotate) ? Tau_lev[TauType::tau33]->array(mfi) : Array4<Real>{};
546  auto t12_arr = (m_rotate) ? Tau_lev[TauType::tau12]->array(mfi) : Array4<Real>{};
547  auto t21_arr = (m_rotate) ? Tau_lev[TauType::tau21]->array(mfi) : Array4<Real>{};
548 
549  auto hfx1_arr = (m_rotate) ? xheat_flux->array(mfi) : Array4<Real>{};
550  auto hfx2_arr = (m_rotate) ? yheat_flux->array(mfi) : Array4<Real>{};
551  auto qfx1_arr = (m_rotate && xqv_flux) ? xqv_flux->array(mfi) : Array4<Real>{};
552  auto qfx2_arr = (m_rotate && yqv_flux) ? yqv_flux->array(mfi) : Array4<Real>{};
553 
554  // Terrain
555  const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};
556 
557  // Get average arrays
558  const auto *const u_mean = m_ma.get_average(lev,0);
559  const auto *const v_mean = m_ma.get_average(lev,1);
560  const auto *const t_mean = m_ma.get_average(lev,2);
561  const auto *const q_mean = m_ma.get_average(lev,3);
562  const auto *const u_mag_mean = m_ma.get_average(lev,5);
563 
564  const auto um_arr = u_mean->array(mfi);
565  const auto vm_arr = v_mean->array(mfi);
566  const auto tm_arr = t_mean->array(mfi);
567  const auto qm_arr = q_mean->array(mfi);
568  const auto umm_arr = u_mag_mean->array(mfi);
569 
570  // Get derived arrays
571  const auto u_star_arr = u_star[lev]->array(mfi);
572  const auto t_star_arr = t_star[lev]->array(mfi);
573  const auto q_star_arr = q_star[lev]->array(mfi);
574  const auto t_surf_arr = t_surf[lev]->array(mfi);
575  const auto q_surf_arr = q_surf[lev]->array(mfi);
576  auto surface_source_arr = surface_diagnostic_source[lev]->array(mfi);
577 
578  // Get LSM fluxes
579  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
580  Array4<int> {};
581  auto lsm_t_flux_arr = Array4<Real> {};
582  auto soil_t_flux_arr = Array4<Real> {};
583  auto lsm_q_flux_arr = Array4<Real> {};
584  auto lsm_tau13_arr = Array4<Real> {};
585  auto lsm_tau23_arr = Array4<Real> {};
586  // LSM tau fields are cell-centered kinematic stresses [m2 s-2].
587  // Tau_lev tau13/tau23 are face-centered conservative stresses [N m-2].
588  for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
589  if (toLower(m_lsm_flux_name[n]) == "t_flux") { lsm_t_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
590  if (toLower(m_lsm_flux_name[n]) == "soil_t_flux") { soil_t_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
591  if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
592  if (toLower(m_lsm_flux_name[n]) == "tau13") { lsm_tau13_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
593  if (toLower(m_lsm_flux_name[n]) == "tau23") { lsm_tau23_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
594  }
595 
596  const bool has_lsm_t_flux = static_cast<bool>(lsm_t_flux_arr);
597  const bool is_custom = (flux_type == FluxCalcType::CUSTOM);
598  const bool is_rico = (flux_type == FluxCalcType::RICO);
599 
600 
601  // Rho*Theta flux
602  //============================================================================
603  Box bx = mfi.tilebox();
604  if (bx.smallEnd(2) != klo) { continue; }
605  bx.makeSlab(2,klo);
606  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
607  {
608  // Valid theta flux from LSM and over land. The LSM writes the
609  // lsm_undefined sentinel for cells it did not process (sea-ice /
610  // open water); fall back to MOST there instead of applying garbage.
611  Real Tflux;
612  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
613  const bool lsm_flux_is_valid = (lsm_t_flux_arr) ? (lsm_t_flux_arr(i,j,0) < lsm_undefined) :
614  false;
615  const bool has_land_and_flux = (is_land == 1 && lsm_flux_is_valid);
616  if (lsm_t_flux_arr && has_land_and_flux) {
617  // LSM flux MultiFabs store kinematic fluxes for MOST parameter
618  // updates. The applied hfx array stores the conservative RHS flux.
619  Tflux = cons_arr(i,j,k,Rho_comp) * lsm_t_flux_arr(i,j,0);
620  } else if (is_land == 2) { // no temperature flux within buildings
621  Tflux = zero;
622  } else {
623  Tflux = flux_comp.compute_t_flux(i, j, k,
624  cons_arr, velx_arr, vely_arr,
625  umm_arr, tm_arr, u_star_arr,
626  t_star_arr, t_surf_arr);
627  // NOTE: do NOT write the MOST-fallback flux back into lsm_t_flux_arr.
628  // Doing so flips a sentinel (water/unprocessed) cell to "valid LSM"
629  // on the next step, so a MOST-derived value is re-read as an LSM flux
630  // Only Noah-MP should populate the LSM cache.
631  }
632 
633  if (soil_t_flux_arr && is_land == 1) {
634  soil_t_flux_arr(i,j,k) = Tflux / cons_arr(i,j,k,Rho_comp);
635  }
636 
637  surface_source_arr(i,j,0) = surface_diagnostics::to_plot_value(
638  surface_diagnostics::classify_scalar_source(
639  is_custom, is_rico, is_land, has_lsm_t_flux, lsm_flux_is_valid));
640 
641  // Do scalar flux rotations?
642  if (rotate) {
643  rotate_scalar_flux(i, j, k, Tflux, dxInv, zphys_arr,
644  hfx1_arr, hfx2_arr, hfx3_arr);
645  } else {
646  hfx3_arr(i,j,klo) = Tflux;
647  }
648  });
649 
650  // Rho*Qv flux
651  //============================================================================
652  if (use_moisture) {
653  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
654  {
655  // Valid qv flux from LSM and over land (sentinel -> fall back to MOST)
656  Real Qflux;
657  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
658  const bool lsm_flux_is_valid = (lsm_q_flux_arr) ? (lsm_q_flux_arr(i,j,0) < lsm_undefined) :
659  false;
660  const bool has_land_and_flux = (is_land == 1 && lsm_flux_is_valid);
661  if (lsm_q_flux_arr && has_land_and_flux) {
662  // LSM flux MultiFabs store kinematic fluxes for MOST parameter
663  // updates. The applied qfx array stores the conservative RHS flux.
664  Qflux = cons_arr(i,j,k,Rho_comp) * lsm_q_flux_arr(i,j,0);
665  } else if (is_land == 2) { // no moisture flux within buildings
666  Qflux = zero;
667  } else {
668  Qflux = flux_comp.compute_q_flux(i, j, k,
669  cons_arr, velx_arr, vely_arr,
670  umm_arr, qm_arr, u_star_arr,
671  q_star_arr, q_surf_arr);
672  // NOTE: no writeback into lsm_q_flux_arr -- see the matching
673  // t_flux note above.
674  }
675 
676  // Do scalar flux rotations?
677  if (rotate) {
678  rotate_scalar_flux(i, j, k, Qflux, dxInv, zphys_arr,
679  qfx1_arr, qfx2_arr, qfx3_arr);
680  } else {
681  qfx3_arr(i,j,k) = Qflux;
682  }
683  });
684  } // custom
685 
686  if (!rotate) {
687  // Rho*u flux
688  //============================================================================
689  Box bxx = surroundingNodes(bx,0);
690  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
691  {
692  // Valid tau13 from LSM and over land. A side that is land but
693  // whose LSM flux is the sentinel (sea-ice / open water) is treated
694  // as non-LSM so that side uses the MOST stress instead.
695  Real stressx;
696  int is_land_hi = (lmask_arr) ? lmask_arr(i ,j,0) : 1;
697  int is_land_lo = (lmask_arr) ? lmask_arr(i-1,j,0) : 1;
698  const bool lsm_hi_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
699  static_cast<bool>(lsm_tau13_arr), is_land_hi == 1,
700  lsm_tau13_arr ? lsm_tau13_arr(i ,j,0) : zero, lsm_undefined);
701  const bool lsm_lo_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
702  static_cast<bool>(lsm_tau13_arr), is_land_lo == 1,
703  lsm_tau13_arr ? lsm_tau13_arr(i-1,j,0) : zero, lsm_undefined);
704  const bool has_land_and_flux_hi = (is_land_hi == 1 && lsm_hi_flux_is_valid);
705  const bool has_land_and_flux_lo = (is_land_lo == 1 && lsm_lo_flux_is_valid);
706  if (lsm_tau13_arr && (has_land_and_flux_hi || has_land_and_flux_lo)) {
707  const Real rho_hi = cons_arr(i ,j,k,Rho_comp);
708  const Real rho_lo = cons_arr(i-1,j,k,Rho_comp);
709  const Real most_stress = (!has_land_and_flux_hi || !has_land_and_flux_lo) ?
710  flux_comp.compute_u_flux(i, j, k,
711  cons_arr, velx_arr, vely_arr,
712  umm_arr, um_arr, u_star_arr) : zero;
713  const auto result = surface_layer_stress::combine_lsm_and_most_stress(
714  rho_lo, rho_hi, lsm_tau13_arr(i-1,j,0), lsm_tau13_arr(i,j,0),
715  has_land_and_flux_lo, has_land_and_flux_hi, most_stress);
716  stressx = result.face_stress;
717  // NOTE: do NOT write the MOST-fallback stress back into the
718  // cell-centered lsm_tau13_arr. This face-indexed ParallelFor
719  // touches cells (i) and (i-1), so each cell is written by two
720  // adjacent face threads in the same launch -> nondeterministic
721  // write-write race on GPU (ERF #3446). It also spuriously flips
722  // a sentinel (water/unprocessed) cell to "valid LSM" for the
723  // next step. The face stress is fully determined here; the LSM
724  // cache is (re)filled only by Noah-MP. Matches baseline 3ab899d3.
725  } else if (is_land_hi == 2 || is_land_lo == 2) { // no stress within buildings
726  stressx = zero;
727  } else {
728  stressx = flux_comp.compute_u_flux(i, j, k,
729  cons_arr, velx_arr, vely_arr,
730  umm_arr, um_arr, u_star_arr);
731  }
732 
733  t13_arr(i,j,k) = stressx;
734  if (t31_arr) { t31_arr(i,j,k) = stressx; }
735  });
736 
737  // Rho*v flux
738  //============================================================================
739  Box bxy = surroundingNodes(bx,1);
740  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
741  {
742  // Valid tau23 from LSM and over land (sentinel side -> MOST stress)
743  Real stressy;
744  int is_land_hi = (lmask_arr) ? lmask_arr(i,j ,0) : 1;
745  int is_land_lo = (lmask_arr) ? lmask_arr(i,j-1,0) : 1;
746  const bool lsm_hi_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
747  static_cast<bool>(lsm_tau23_arr), is_land_hi == 1,
748  lsm_tau23_arr ? lsm_tau23_arr(i,j ,0) : zero, lsm_undefined);
749  const bool lsm_lo_flux_is_valid = surface_layer_stress::lsm_flux_is_valid(
750  static_cast<bool>(lsm_tau23_arr), is_land_lo == 1,
751  lsm_tau23_arr ? lsm_tau23_arr(i,j-1,0) : zero, lsm_undefined);
752  const bool has_land_and_flux_hi = (is_land_hi == 1 && lsm_hi_flux_is_valid);
753  const bool has_land_and_flux_lo = (is_land_lo == 1 && lsm_lo_flux_is_valid);
754  if (lsm_tau23_arr && (has_land_and_flux_hi || has_land_and_flux_lo)) {
755  const Real rho_hi = cons_arr(i,j ,k,Rho_comp);
756  const Real rho_lo = cons_arr(i,j-1,k,Rho_comp);
757  const Real most_stress = (!has_land_and_flux_hi || !has_land_and_flux_lo) ?
758  flux_comp.compute_v_flux(i, j, k,
759  cons_arr, velx_arr, vely_arr,
760  umm_arr, vm_arr, u_star_arr) : zero;
761  const auto result = surface_layer_stress::combine_lsm_and_most_stress(
762  rho_lo, rho_hi, lsm_tau23_arr(i,j-1,0), lsm_tau23_arr(i,j,0),
763  has_land_and_flux_lo, has_land_and_flux_hi, most_stress);
764  stressy = result.face_stress;
765  // NOTE: no writeback into cell-centered lsm_tau23_arr -- see the
766  // matching tau13 note above (ERF #3446 write-write race + stale
767  // sentinel-becomes-valid). Face stress is complete here.
768  } else if (is_land_hi == 2 || is_land_lo == 2) { // no stress within buildings
769  stressy = zero;
770  } else {
771  stressy = flux_comp.compute_v_flux(i, j, k,
772  cons_arr, velx_arr, vely_arr,
773  umm_arr, vm_arr, u_star_arr);
774  }
775 
776  t23_arr(i,j,k) = stressy;
777  if (t32_arr) { t32_arr(i,j,k) = stressy; }
778  });
779  } else {
780  // All fluxes with rotation
781  //============================================================================
782  Box bxxy = convert(bx, IntVect(1,1,0));
783  ParallelFor(bxxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
784  {
785  Real stresst = flux_comp.compute_u_flux(i, j, k,
786  cons_arr, velx_arr, vely_arr,
787  umm_arr, um_arr, u_star_arr);
788  rotate_stress_tensor(i, j, k, stresst, dxInv, zphys_arr,
789  velx_arr, vely_arr, velz_arr,
790  t11_arr, t22_arr, t33_arr,
791  t12_arr, t21_arr,
792  t13_arr, t31_arr,
793  t23_arr, t32_arr);
794  });
795  }
796 
797  // For models that do not do iterations to yield u*/T*/q*,
798  // fill these values from the fluxes that were computed.
799 
800  // NOTE: For LSM, this has been handled in "compute_sfc_params_from_lsm_fluxes"
801  // NOTE: Fluxes here are for conserved quantities, we divide by rho
802  if (flux_type == FluxCalcType::BULK_COEFF ||
803  flux_type == FluxCalcType::DONELAN) {
804  constexpr Real eps = std::numeric_limits<Real>::epsilon();
805  bool l_use_moisture = use_moisture;
806  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/)
807  {
808  Real rho = cons_arr(i,j,klo,Rho_comp);
809  Real Thd = cons_arr(i,j,klo,RhoTheta_comp) / rho;
810  Real qv = (l_use_moisture) ? cons_arr(i,j,klo,RhoQ1_comp) / rho : zero;
811  Real Thv = Thd * (one + epsv*qv);
812 
813  Real tau = std::sqrt( t13_arr(i,j,klo)/rho * t13_arr(i,j,klo)/rho
814  + t23_arr(i,j,klo)/rho * t23_arr(i,j,klo)/rho );
815  u_star_arr(i,j,0) = amrex::max(std::sqrt(tau),eps);
816 
817  if (hfx3_arr(i,j,klo)>=zero) {
818  t_star_arr(i,j,0) = amrex::min(-hfx3_arr(i,j,klo) / (rho * u_star_arr(i,j,0)),-eps);
819  } else {
820  t_star_arr(i,j,0) = amrex::max(-hfx3_arr(i,j,klo) / (rho * u_star_arr(i,j,0)),eps);
821  }
822  if (!l_use_moisture) {
823  q_star_arr(i,j,0) = zero;
824  } else if (qfx3_arr(i,j,klo)>=zero) {
825  q_star_arr(i,j,0) = amrex::min(-qfx3_arr(i,j,klo) / (rho * u_star_arr(i,j,0)),-eps);
826  } else {
827  q_star_arr(i,j,0) = amrex::max(-qfx3_arr(i,j,klo) / ( rho * u_star_arr(i,j,0)),eps);
828  }
829  olen_arr(i,j,0) = ( u_star_arr(i,j,0) * u_star_arr(i,j,0) * Thv ) /
830  ( KAPPA * CONST_GRAV * t_star_arr(i,j,0) );
831  });
832  }
833 
834  } // mfiter
835 
836  surface_diagnostic_source[lev]->FillBoundary(m_geom[lev].periodicity());
837 }
tau12
@ tau12
Definition: ERF_DataStruct.H:38
tau23
@ tau23
Definition: ERF_DataStruct.H:38
tau33
@ tau33
Definition: ERF_DataStruct.H:38
tau22
@ tau22
Definition: ERF_DataStruct.H:38
tau11
@ tau11
Definition: ERF_DataStruct.H:38
tau32
@ tau32
Definition: ERF_DataStruct.H:38
tau31
@ tau31
Definition: ERF_DataStruct.H:38
tau21
@ tau21
Definition: ERF_DataStruct.H:38
tau13
@ tau13
Definition: ERF_DataStruct.H:38
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
rotate_scalar_flux
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void rotate_scalar_flux(const int &i, const int &j, const int &klo, const amrex::Real &flux, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &zphys_arr, const amrex::Array4< amrex::Real > &phi1_arr, const amrex::Array4< amrex::Real > &phi2_arr, const amrex::Array4< amrex::Real > &phi3_arr)
Definition: ERF_TerrainMetrics.H:611
rotate_stress_tensor
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void rotate_stress_tensor(const int &i, const int &j, const int &klo, const amrex::Real &flux, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &zphys_arr, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &w_arr, const amrex::Array4< amrex::Real > &tau11_arr, const amrex::Array4< amrex::Real > &tau22_arr, const amrex::Array4< amrex::Real > &tau33_arr, const amrex::Array4< amrex::Real > &tau12_arr, const amrex::Array4< amrex::Real > &tau21_arr, const amrex::Array4< amrex::Real > &tau13_arr, const amrex::Array4< amrex::Real > &tau31_arr, const amrex::Array4< amrex::Real > &tau23_arr, const amrex::Array4< amrex::Real > &tau32_arr)
Definition: ERF_TerrainMetrics.H:632
SurfaceLayer::surface_diagnostic_source
amrex::Vector< std::unique_ptr< amrex::MultiFab > > surface_diagnostic_source
Definition: ERF_SurfaceLayer.H:1073
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:177
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:179
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:178
surface_diagnostics::to_plot_value
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real to_plot_value(SurfaceDiagnosticSource source) noexcept
Definition: ERF_SurfaceDiagnosticSource.H:45
surface_diagnostics::classify_scalar_source
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SurfaceDiagnosticSource classify_scalar_source(bool is_custom, bool is_rico, bool is_land, bool has_lsm_flux, bool lsm_flux_is_valid) noexcept
Definition: ERF_SurfaceDiagnosticSource.H:61
surface_layer_stress::lsm_flux_is_valid
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool lsm_flux_is_valid(bool has_flux, bool is_land, amrex::Real flux, amrex::Real undefined)
Definition: ERF_SurfaceLayerStress.H:35
surface_layer_stress::combine_lsm_and_most_stress
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE FaceStressResult combine_lsm_and_most_stress(amrex::Real rho_low, amrex::Real rho_high, amrex::Real kinematic_low, amrex::Real kinematic_high, bool low_valid, bool high_valid, amrex::Real most_face_stress)
Definition: ERF_SurfaceLayerStress.H:64
Here is the call graph for this function:

◆ compute_SurfaceLayer_bcs_EB() [1/2]

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

Compute embedded-boundary surface-layer flux boundary conditions.

Parameters
[in]levlevel index
[in]mfsstate and velocity fields used by the BC computation
[in,out]Tau_levEB stress fields to fill
[in,out]xheat_fluxx-face heat flux field
[in,out]yheat_fluxy-face heat flux field
[in,out]zheat_fluxz-face heat flux field
[in,out]xqv_fluxx-face moisture flux field
[in,out]yqv_fluxy-face moisture flux field
[in,out]zqv_fluxz-face moisture flux field
[in]flux_compflux-computation functor

◆ compute_SurfaceLayer_bcs_EB() [2/2]

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

Function to calculate MOST fluxes for EB.

Parameters
[in]levCurrent level
[in]mfsState MultiFabs used to compute the EB boundary fluxes
[in,out]Tau_EBEB diffusive stress MultiFabs populated with surface stresses
[in,out]xheat_fluxx-face EB heat-flux MultiFab, currently unused
[in,out]yheat_fluxy-face EB heat-flux MultiFab, currently unused
[in,out]Hfx3_EBEB heat-flux MultiFab populated with scalar surface flux
[in,out]xqv_fluxx-face EB moisture-flux MultiFab, currently unused
[in,out]yqv_fluxy-face EB moisture-flux MultiFab, currently unused
[in,out]zqv_fluxz-face EB moisture-flux MultiFab, currently unused
[in]flux_compEB flux-calculation functor used to compute scalar and momentum fluxes
865 {
866  // Get EB flags for all centerings
867  const auto& cc_factory = m_eb_vec[lev]->get_const_factory();
868  const auto& cc_flags = cc_factory->getMultiEBCellFlagFab();
869  const auto& cc_vfrac = cc_factory->getVolFrac();
870 
871  const auto& u_factory = m_eb_vec[lev]->get_u_const_factory();
872  const auto& u_flags = u_factory->getMultiEBCellFlagFab();
873  const auto& u_vfrac = u_factory->getVolFrac();
874 
875  const auto& v_factory = m_eb_vec[lev]->get_v_const_factory();
876  const auto& v_flags = v_factory->getMultiEBCellFlagFab();
877  const auto& v_vfrac = v_factory->getVolFrac();
878 
879  const auto& w_factory = m_eb_vec[lev]->get_w_const_factory();
880  const auto& w_flags = w_factory->getMultiEBCellFlagFab();
881  const auto& w_vfrac = w_factory->getVolFrac();
882 
883  // EB does not currently have a cell-centered scalar-source classification.
884  // Keep the provenance mask missing rather than inventing face-aware
885  // semantics for the staggered stress path.
886  surface_diagnostic_source[lev]->setVal(
887  surface_diagnostics::to_plot_value(surface_diagnostics::SurfaceDiagnosticSource::Missing));
888 
889  for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
890  {
891  // Get flags for this box (all centerings)
892  const auto& cc_flag = cc_flags[mfi];
893  const auto& u_flag = u_flags[mfi];
894  const auto& v_flag = v_flags[mfi];
895  const auto& w_flag = w_flags[mfi];
896 
897  // Skip boxes that have no cut cells at any centering
898  if (cc_flag.getType() != FabType::singlevalued &&
899  u_flag.getType() != FabType::singlevalued &&
900  v_flag.getType() != FabType::singlevalued &&
901  w_flag.getType() != FabType::singlevalued
902  ) continue;
903 
904  // Get EB flag and volfrac arrays
905  auto const cc_flag_arr = cc_flag.const_array();
906  auto const u_flag_arr = u_flag.const_array();
907  auto const v_flag_arr = v_flag.const_array();
908  auto const w_flag_arr = w_flag.const_array();
909 
910  auto const cc_vfrac_arr = cc_vfrac.const_array(mfi);
911  auto const u_vfrac_arr = u_vfrac.const_array(mfi);
912  auto const v_vfrac_arr = v_vfrac.const_array(mfi);
913  auto const w_vfrac_arr = w_vfrac.const_array(mfi);
914 
915  // Get boundary normals only if cut cells exist
916  auto const bnorm_arr = (cc_flag.getType() == FabType::singlevalued) ?
917  cc_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
918  auto const u_bnorm_arr = (u_flag.getType() == FabType::singlevalued) ?
919  u_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
920  auto const v_bnorm_arr = (v_flag.getType() == FabType::singlevalued) ?
921  v_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
922  auto const w_bnorm_arr = (w_flag.getType() == FabType::singlevalued) ?
923  w_factory->getBndryNormal().const_array(mfi) : Array4<const Real>{};
924 
925  // Get field arrays
926  const auto cons_arr = mfs[Vars::cons]->array(mfi);
927  const auto velx_arr = mfs[Vars::xvel]->array(mfi);
928  const auto vely_arr = mfs[Vars::yvel]->array(mfi);
929  const auto velz_arr = mfs[Vars::zvel]->array(mfi);
930 
931  // Diffusive stress vars - t13 and t23 components for all grid types
932  auto u_t13_arr = Tau_EB[EBTauType::tau_eb13][EBGridType::xface]->array(mfi);
933  auto v_t13_arr = Tau_EB[EBTauType::tau_eb13][EBGridType::yface]->array(mfi);
934  auto w_t13_arr = Tau_EB[EBTauType::tau_eb13][EBGridType::zface]->array(mfi);
935 
936  auto u_t23_arr = Tau_EB[EBTauType::tau_eb23][EBGridType::xface]->array(mfi);
937  auto v_t23_arr = Tau_EB[EBTauType::tau_eb23][EBGridType::yface]->array(mfi);
938  auto w_t23_arr = Tau_EB[EBTauType::tau_eb23][EBGridType::zface]->array(mfi);
939 
940  auto hfx3_arr = Hfx3_EB->array(mfi);
941 
942  // Get average arrays
943  const auto *const u_mean = m_ma.get_average(lev,0);
944  const auto *const v_mean = m_ma.get_average(lev,1);
945  const auto *const t_mean = m_ma.get_average(lev,2);
946  // const auto *const q_mean = m_ma.get_average(lev,3);
947  const auto *const u_mag_mean = m_ma.get_average(lev,5);
948 
949  const auto um_arr = u_mean->array(mfi);
950  const auto vm_arr = v_mean->array(mfi);
951  const auto tm_arr = t_mean->array(mfi);
952  // const auto qm_arr = q_mean->array(mfi);
953  const auto umm_arr = u_mag_mean->array(mfi);
954 
955  // Get derived arrays
956  const auto u_star_arr = u_star[lev]->array(mfi);
957  const auto t_star_arr = t_star[lev]->array(mfi);
958  const auto t_surf_arr = t_surf[lev]->array(mfi);
959 
960  // Rho*Theta flux
961  //============================================================================
962  Box bx = mfi.tilebox();
963  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
964  {
965  if (cc_flag_arr(i,j,k).isSingleValued()) {
966  Real Tflux = flux_comp.compute_t_flux(i, j, k,
967  cons_arr, velx_arr, vely_arr, velz_arr,
968  umm_arr, tm_arr, u_star_arr,
969  t_star_arr, t_surf_arr,
970  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
971  bnorm_arr);
972  hfx3_arr(i,j,k) = Tflux;
973  }
974  });
975 
976  // Rho*u flux
977  //============================================================================
978  Box bxx = surroundingNodes(bx,0);
979  Box bxy = surroundingNodes(bx,1);
980  Box bxz = surroundingNodes(bx,2);
981  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
982  {
983  if (u_flag_arr(i,j,k).isSingleValued()) {
984  Real stressx = flux_comp.compute_u_flux(i, j, k,
985  cons_arr, velx_arr, vely_arr, velz_arr,
986  umm_arr, um_arr, u_star_arr,
987  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
988  cc_vfrac_arr, cc_flag_arr,
989  u_bnorm_arr, 0);
990  u_t13_arr(i,j,k) = stressx;
991  }
992  });
993  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
994  {
995  if (v_flag_arr(i,j,k).isSingleValued()) {
996  Real stressx = flux_comp.compute_u_flux(i, j, k,
997  cons_arr, velx_arr, vely_arr, velz_arr,
998  umm_arr, um_arr, u_star_arr,
999  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
1000  cc_vfrac_arr, cc_flag_arr,
1001  v_bnorm_arr, 1);
1002  v_t13_arr(i,j,k) = stressx;
1003  }
1004  });
1005  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
1006  {
1007  if (w_flag_arr(i,j,k).isSingleValued()) {
1008  Real stressx = flux_comp.compute_u_flux(i, j, k,
1009  cons_arr, velx_arr, vely_arr, velz_arr,
1010  umm_arr, um_arr, u_star_arr,
1011  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
1012  cc_vfrac_arr, cc_flag_arr,
1013  w_bnorm_arr, 2);
1014  w_t13_arr(i,j,k) = stressx;
1015  }
1016  });
1017 
1018  // Rho*v flux
1019  //============================================================================
1020  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
1021  {
1022  if (u_flag_arr(i,j,k).isSingleValued()) {
1023  Real stressy = flux_comp.compute_v_flux(i, j, k,
1024  cons_arr, velx_arr, vely_arr, velz_arr,
1025  umm_arr, vm_arr, u_star_arr,
1026  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
1027  cc_vfrac_arr, cc_flag_arr, u_bnorm_arr, 0);
1028  u_t23_arr(i,j,k) = stressy;
1029  }
1030  });
1031  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
1032  {
1033  if (v_flag_arr(i,j,k).isSingleValued()) {
1034  Real stressy = flux_comp.compute_v_flux(i, j, k,
1035  cons_arr, velx_arr, vely_arr, velz_arr,
1036  umm_arr, vm_arr, u_star_arr,
1037  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
1038  cc_vfrac_arr, cc_flag_arr, v_bnorm_arr, 1);
1039  v_t23_arr(i,j,k) = stressy;
1040  }
1041  });
1042  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
1043  {
1044  if (w_flag_arr(i,j,k).isSingleValued()) {
1045  Real stressy = flux_comp.compute_v_flux(i, j, k,
1046  cons_arr, velx_arr, vely_arr, velz_arr,
1047  umm_arr, vm_arr, u_star_arr,
1048  u_vfrac_arr, v_vfrac_arr, w_vfrac_arr,
1049  cc_vfrac_arr, cc_flag_arr, w_bnorm_arr, 2);
1050  w_t23_arr(i,j,k) = stressy;
1051  }
1052  });
1053  } // mfiter
1054 
1055 }
tau_eb23
@ tau_eb23
Definition: ERF_EBStruct.H:22
tau_eb13
@ tau_eb13
Definition: ERF_EBStruct.H:22
yface
@ yface
Definition: ERF_EBStruct.H:29
zface
@ zface
Definition: ERF_EBStruct.H:29
xface
@ xface
Definition: ERF_EBStruct.H:29
surface_diagnostics::SurfaceDiagnosticSource::Missing
@ Missing
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 
)

Fill surface moisture from saturation specific humidity.

Parameters
[in]levlevel index
[in]cons_inconserved state used to evaluate surface pressure
[in]z_phys_ndnodal physical-height field

Fill sea-surface moisture with saturation specific humidity.

Parameters
[in]levCurrent level
[in]cons_inConserved state used to derive pressure at the surface
[in]z_phys_ndNodal physical height used to compute terrain-relative surface height
1234 {
1235  // NOTE: We have already tested a moisture model exists
1236 
1237  // Populate q_surf with qsat over water
1238  auto dz = m_geom[lev].CellSize(2);
1239  const int klo = m_geom[lev].Domain().smallEnd(2);
1240  for (MFIter mfi(*q_surf[lev]); mfi.isValid(); ++mfi)
1241  {
1242  Box gtbx = mfi.growntilebox();
1243 
1244  if (gtbx.smallEnd(2) != klo) { continue; }
1245 
1246  auto t_surf_arr = t_surf[lev]->array(mfi);
1247  auto q_surf_arr = q_surf[lev]->array(mfi);
1248  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
1249  Array4<int> {};
1250  const auto cons_arr = cons_in.const_array(mfi);
1251  const auto z_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
1252  Array4<const Real> {};
1253 
1254  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1255  {
1256  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
1257  if (!is_land) {
1258  auto deltaZ = (z_arr) ? Compute_Zrel_AtCellCenter(i,j,k,z_arr) :
1259  myhalf*dz;
1260  auto Rho = cons_arr(i,j,k,Rho_comp);
1261  auto RTh = cons_arr(i,j,k,RhoTheta_comp);
1262  auto Qv = cons_arr(i,j,k,RhoQ1_comp) / Rho;
1263  auto P_cc = getPgivenRTh(RTh, Qv);
1264  P_cc += Rho*CONST_GRAV*deltaZ;
1265  P_cc *= Real(0.01);
1266  erf_qsatw(t_surf_arr(i,j,k), P_cc, q_surf_arr(i,j,k));
1267  }
1268  });
1269  }
1270  q_surf[lev]->FillBoundary(m_geom[lev].periodicity());
1271 }
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
erf_qsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:228
Compute_Zrel_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Zrel_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:389
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
Here is the call graph for this function:

◆ fill_tsurf_with_sst_and_tsk()

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

Fill surface temperature from available SST and TSK data.

Parameters
[in]levlevel index
[in]timeinterpolation time

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

Parameters
[in]levCurrent level
[in]elapsed_time_since_start_lowTime since the start of the lower-boundary data
1139 {
1140  int n_times_in_sst = static_cast<int>(m_sst_lev[lev].size());
1141 
1142  double dT = m_low_time_interval;
1143 
1144  int n_time_lo, n_time_hi;
1145  Real alpha;
1146 
1147  if (n_times_in_sst > 1) {
1148  n_time_lo = static_cast<int>( elapsed_time_since_start_low / dT);
1149  alpha = static_cast<Real>((elapsed_time_since_start_low - n_time_lo * dT) / dT);
1150 
1151  AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one);
1152 
1153  n_time_hi = n_time_lo + 1;
1154 
1155  // Do not over run the last sst file
1156  if (m_start_low_time + elapsed_time_since_start_low >= m_final_low_time) {
1157  n_time_lo = static_cast<int>(m_sst_lev[lev].size())-1;
1158  n_time_hi = n_time_lo;
1159  alpha = zero;
1160  }
1161 
1162  AMREX_ALWAYS_ASSERT( (n_time_lo >= 0) && (n_time_hi < m_sst_lev[lev].size()));
1163  } else {
1164  n_time_lo = 0;
1165  n_time_hi = 0;
1166  alpha = one;
1167  }
1168  AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one);
1169 
1170  Real oma = one - alpha;
1171 
1172  // Define a default land surface temperature if we don't read in tsk
1173  Real lst = default_land_surf_temp;
1174 
1175  bool use_tsk = (m_tsk_lev[lev][0]);
1176  bool ignore_sst = m_ignore_sst;
1177 
1178  const int klo = m_geom[lev].Domain().smallEnd(2);
1179 
1180  // Populate t_surf
1181  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
1182  {
1183  Box gtbx = mfi.growntilebox();
1184 
1185  if (gtbx.smallEnd(2) != klo) { continue; }
1186 
1187  auto t_surf_arr = t_surf[lev]->array(mfi);
1188 
1189  const auto sst_lo_arr = m_sst_lev[lev][n_time_lo]->const_array(mfi);
1190  const auto sst_hi_arr = m_sst_lev[lev][n_time_hi]->const_array(mfi);
1191 
1192  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
1193  Array4<int> {};
1194 
1195  if (use_tsk) {
1196  const auto tsk_arr = m_tsk_lev[lev][n_time_lo]->const_array(mfi);
1197  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1198  {
1199  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
1200  if (!is_land && !ignore_sst) {
1201  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
1202  + alpha * sst_hi_arr(i,j,k);
1203  } else {
1204  t_surf_arr(i,j,k) = tsk_arr(i,j,k);
1205  }
1206  });
1207  } else {
1208  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1209  {
1210  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
1211  if (!is_land) {
1212  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
1213  + alpha * sst_hi_arr(i,j,k);
1214  } else {
1215  t_surf_arr(i,j,k) = lst;
1216  }
1217  });
1218  }
1219  }
1220  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
1221 }
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
SurfaceLayer::m_sst_lev
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_sst_lev
Definition: ERF_SurfaceLayer.H:1075
SurfaceLayer::m_tsk_lev
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_tsk_lev
Definition: ERF_SurfaceLayer.H:1076
mp_wsm6::alpha
real(kind=kind_phys), parameter, public alpha
Definition: ERF_module_mp_wsm6.F90:44
Here is the call graph for this function:

◆ get_lmask()

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

Return the land-mask field.

Parameters
[in]levlevel index
928 { return m_lmask_lev[lev][0]; }

◆ get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf ( const int &  lev)

Fill surface temperature from the LSM surface-temperature field.

Parameters
[in]levlevel index

Fill surface temperature from land-surface-model data.

Parameters
[in]levCurrent level
1280 {
1281  const int klo = m_geom[lev].Domain().smallEnd(2);
1282  const bool has_sea_tsurf = (m_has_ocean_lsm_tsurf &&
1283  amrex::toLower(m_lsm_data_name[m_lsm_tsurf_indx]) == "t_surf");
1284  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
1285  {
1286  Box gtbx = mfi.growntilebox();
1287 
1288  if (gtbx.smallEnd(2) != klo) { continue; }
1289 
1290  // NOTE: LSM does not carry lateral ghost cells.
1291  // This copies the valid box into the ghost cells.
1292  // Fillboundary is called after this to pick up the
1293  // interior ghost and periodic directions.
1294  Box vbx = mfi.validbox();
1295  int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
1296  int j_lo = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1);
1297 
1298  auto t_surf_arr = t_surf[lev]->array(mfi);
1299  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
1300  Array4<int> {};
1301  const auto lsm_arr = m_lsm_data_lev[lev][m_lsm_tsurf_indx]->const_array(mfi);
1302 
1303  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1304  {
1305  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
1306  if ((!has_sea_tsurf && is_land) ||
1307  (has_sea_tsurf && !is_land)) {
1308  int li = amrex::min(amrex::max(i, i_lo), i_hi);
1309  int lj = amrex::min(amrex::max(j, j_lo), j_hi);
1310  t_surf_arr(i,j,k) = lsm_arr(li,lj,k);
1311  }
1312  });
1313  }
1314 }
SurfaceLayer::m_lsm_data_name
amrex::Vector< std::string > m_lsm_data_name
Definition: ERF_SurfaceLayer.H:1080
SurfaceLayer::m_lsm_tsurf_indx
int m_lsm_tsurf_indx
Definition: ERF_SurfaceLayer.H:1051
SurfaceLayer::m_lsm_data_lev
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_data_lev
Definition: ERF_SurfaceLayer.H:1078
SurfaceLayer::m_has_ocean_lsm_tsurf
bool m_has_ocean_lsm_tsurf
Definition: ERF_SurfaceLayer.H:1050
Here is the call graph for this function:

◆ get_mac_avg()

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

Return a MOST-average field.

Parameters
[in]levlevel index
[in]compcomponent index
863  {
864  return m_ma.get_average(lev, comp);
865  }
Here is the call graph for this function:

◆ get_olen()

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

Return the Obukhov length field.

Parameters
[in]levlevel index
847 { return olen[lev].get(); }

◆ get_pblh()

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

Return the planetary-boundary-layer-height field.

Parameters
[in]levlevel index
854 { return pblh[lev].get(); }

◆ get_q_star()

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

Return the moisture scale field.

Parameters
[in]levlevel index
840 { return q_star[lev].get(); }

◆ get_q_surf()

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

Return the surface-moisture field.

Parameters
[in]levlevel index
887 { return q_surf[lev].get(); }

◆ get_surface_diagnostic_source()

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

Return the surface-diagnostic provenance field.

Parameters
[in]levlevel index
902 { return surface_diagnostic_source[lev].get(); }

◆ get_t_star()

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

Return the temperature scale field.

Parameters
[in]levlevel index
833 { return t_star[lev].get(); }

◆ get_t_surf()

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

Return the surface-temperature field.

Parameters
[in]levlevel index
872 { return t_surf[lev].get(); }

◆ get_u_star()

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

Return the friction-velocity field.

Parameters
[in]levlevel index
819 { return u_star[lev].get(); }

◆ get_w_star()

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

Return the convective velocity scale field.

Parameters
[in]levlevel index
826 { return w_star[lev].get(); }

◆ get_z0()

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

Return the roughness-height field.

Parameters
[in]levlevel index
916 { return &z_0[lev]; }

◆ get_zref()

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

Return the minimum reference height for one level.

Parameters
[in]levlevel index
909 { return (m_ma.get_zref(lev))->min(0); }
Here is the call graph for this function:

◆ have_variable_sea_roughness()

bool SurfaceLayer::have_variable_sea_roughness ( )
inline

Return whether variable sea roughness is active.

921 { return m_var_z0; }
SurfaceLayer::m_var_z0
bool m_var_z0
Definition: ERF_SurfaceLayer.H:1042

◆ impose_SurfaceLayer_bcs()

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

Impose surface-layer boundary conditions for planar terrain.

Parameters
[in]levlevel index
[in]mfsstate and velocity fields used by the BC computation
[in,out]Tau_levstress fields to fill
[in,out]xheat_fluxx-face heat flux field
[in,out]yheat_fluxy-face heat flux field
[in,out]zheat_fluxz-face heat flux field
[in,out]xqv_fluxx-face moisture flux field
[in,out]yqv_fluxy-face moisture flux field
[in,out]zqv_fluxz-face moisture flux field
[in]z_physphysical-height field

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

Parameters
[in]levCurrent level
[in]mfsState MultiFabs used to compute the boundary fluxes
[in,out]Tau_levDiffusive stress MultiFabs populated with surface stresses
[in,out]xheat_fluxx-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]yheat_fluxy-face heat-flux MultiFab, used when rotated fluxes are enabled
[in,out]zheat_fluxz-face heat-flux MultiFab populated with vertical surface heat flux
[in,out]xqv_fluxx-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]yqv_fluxy-face moisture-flux MultiFab, used when rotated fluxes and moisture are enabled
[in,out]zqv_fluxz-face moisture-flux MultiFab populated when moisture is enabled
[in]z_physNodal physical height used to rotate terrain-following fluxes
413 {
414  if (flux_type == FluxCalcType::MOENG) {
415  moeng_flux flux_comp;
416  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
417  xheat_flux, yheat_flux, zheat_flux,
418  xqv_flux, yqv_flux, zqv_flux,
419  z_phys, flux_comp);
420  } else if (flux_type == FluxCalcType::DONELAN) {
421  donelan_flux flux_comp;
422  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
423  xheat_flux, yheat_flux, zheat_flux,
424  xqv_flux, yqv_flux, zqv_flux,
425  z_phys, flux_comp);
426  } else if (flux_type == FluxCalcType::ROTATE) {
427  rotate_flux flux_comp;
428  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
429  xheat_flux, yheat_flux, zheat_flux,
430  xqv_flux, yqv_flux, zqv_flux,
431  z_phys, flux_comp);
432  } else if (flux_type == FluxCalcType::RICO) {
433  rico_flux flux_comp(rico_theta_z0, rico_qsat_z0);
434  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
435  xheat_flux, yheat_flux, zheat_flux,
436  xqv_flux, yqv_flux, zqv_flux,
437  z_phys, flux_comp);
438  } else if (flux_type == FluxCalcType::BULK_COEFF) {
439  bulk_coeff_flux flux_comp(m_Cd, m_Ch, m_Cq);
440  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
441  xheat_flux, yheat_flux, zheat_flux,
442  xqv_flux, yqv_flux, zqv_flux,
443  z_phys, flux_comp);
444  } else if (flux_type == FluxCalcType::CUSTOM) {
445  custom_flux flux_comp(specified_rho_surf);
446  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
447  xheat_flux, yheat_flux, zheat_flux,
448  xqv_flux, yqv_flux, zqv_flux,
449  z_phys, flux_comp);
450  } else {
451  amrex::Abort("Unknown surface layer flux calculation type");
452  }
453 }
SurfaceLayer::specified_rho_surf
bool specified_rho_surf
Definition: ERF_SurfaceLayer.H:1037
SurfaceLayer::compute_SurfaceLayer_bcs
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)
bulk_coeff_flux
Definition: ERF_MOSTStress.H:2582
custom_flux
Definition: ERF_MOSTStress.H:2428
donelan_flux
Definition: ERF_MOSTStress.H:2264
moeng_flux
Definition: ERF_MOSTStress.H:2066
rico_flux
Definition: ERF_MOSTStress.H:2744
rotate_flux
Definition: ERF_MOSTStress.H:2910

◆ impose_SurfaceLayer_bcs_EB()

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

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

Parameters
[in]levlevel index
[in]mfsstate and velocity fields used by the BC computation
[in,out]Tau_levEB stress fields to fill
[in,out]xheat_fluxx-face heat flux field
[in,out]yheat_fluxy-face heat flux field
[in,out]zheat_fluxz-face heat flux field
[in,out]xqv_fluxx-face moisture flux field
[in,out]yqv_fluxy-face moisture flux field
[in,out]zqv_fluxz-face moisture flux field

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

Parameters
[in]levCurrent level
[in]mfsState MultiFabs used to compute the EB boundary fluxes
[in,out]Tau_EBEB diffusive stress MultiFabs populated with surface stresses
[in,out]xheat_fluxx-face EB heat-flux MultiFab, currently unused
[in,out]yheat_fluxy-face EB heat-flux MultiFab, currently unused
[in,out]Hfx3_EBEB heat-flux MultiFab populated with scalar surface flux
[in,out]xqv_fluxx-face EB moisture-flux MultiFab, currently unused
[in,out]yqv_fluxy-face EB moisture-flux MultiFab, currently unused
[in,out]zqv_fluxz-face EB moisture-flux MultiFab, currently unused
478 {
479  if (flux_type == FluxCalcType::MOENG) {
480  moeng_flux_eb flux_comp;
481  compute_SurfaceLayer_bcs_EB(lev, mfs, Tau_EB,
482  xheat_flux, yheat_flux, Hfx3_EB,
483  xqv_flux, yqv_flux, zqv_flux,
484  flux_comp);
485  } else {
486  amrex::Abort("Not implemented surface layer flux calculation type for EB");
487  }
488 }
SurfaceLayer::compute_SurfaceLayer_bcs_EB
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)
moeng_flux_eb
EB implementation of the Moeng surface-flux formulation.
Definition: ERF_EBMOSTStress.H:316

◆ init_tke_from_ustar()

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

Initialize TKE from the current surface friction velocity.

Parameters
[in]levlevel index
[in,out]consconserved state whose TKE component is initialized
[in]z_phys_ndnodal physical-height field
[in]tkefacscale factor applied to the initialized TKE
[in]zscalevertical decay scale

Initialize TKE from surface-layer friction velocity.

Parameters
[in]levCurrent level
[in,out]consConserved state whose RhoKE component is initialized
[in]z_phys_ndNodal physical height used to compute height above ground
[in]tkefacFactor multiplying ustar squared for the surface TKE value
[in]zscaleScale factor used to taper TKE with height
1375 {
1376  Print() << "Initializing TKE from surface layer ustar on level " << lev << std::endl;
1377 
1378  // Handle vertical decomposition by selectively copying into
1379  // a FArrayBox section on each rank. Then doing a reduce real sum
1380  // and broadcasting to each rank. No mask since all CC data
1381  const int klo = m_geom[lev].Domain().smallEnd(2);
1382  Box bx_lo = u_star[lev]->boxArray().minimalBox();
1383  FArrayBox u_star_lo(bx_lo, 1); u_star_lo.setVal<RunOn::Device>(0);
1384  FArrayBox z_surf_lo(bx_lo, 1); z_surf_lo.setVal<RunOn::Device>(0);
1385  Real* ustar_ptr = u_star_lo.dataPtr();
1386  Real* zsurf_ptr = z_surf_lo.dataPtr();
1387  for (MFIter mfi(cons); mfi.isValid(); ++mfi)
1388  {
1389  Box vbx = mfi.validbox();
1390  if (vbx.smallEnd(2) != klo) { continue; }
1391  vbx.makeSlab(2,0);
1392 
1393  auto const& u_star_arr = u_star[lev]->const_array(mfi);
1394  auto u_star_all = u_star_lo.array();
1395 
1396  auto const& z_phys_arr = z_phys_nd->const_array(mfi);
1397  auto z_surf_all = z_surf_lo.array();
1398 
1399  ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
1400  {
1401  u_star_all(i,j,0) = u_star_arr(i,j,0);
1402  z_surf_all(i,j,0) = fourth * ( z_phys_arr(i ,j ,klo) + z_phys_arr(i+1,j ,klo)
1403  + z_phys_arr(i ,j+1,klo) + z_phys_arr(i+1,j+1,klo) );
1404  });
1405  }
1406  ParallelDescriptor::ReduceRealSum(ustar_ptr, static_cast<int>(bx_lo.numPts()));
1407  ParallelDescriptor::ReduceRealSum(zsurf_ptr, static_cast<int>(bx_lo.numPts()));
1408 
1409  // Now work on all boxes (ustar has been filled above)
1410  constexpr Real small = Real(0.01);
1411  for (MFIter mfi(cons); mfi.isValid(); ++mfi)
1412  {
1413  Box vbx = mfi.validbox();
1414 
1415  auto const& u_star_arr = u_star_lo.const_array();
1416  auto const& z_surf_arr = z_surf_lo.const_array();
1417  auto const& z_phys_arr = z_phys_nd->const_array(mfi);
1418 
1419  auto const& cons_arr = cons.array(mfi);
1420 
1421  ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1422  {
1423  Real rho = cons_arr(i, j, k, Rho_comp);
1424  Real ust = u_star_arr(i, j, 0);
1425  Real tke0 = tkefac * ust * ust; // surface value
1426  Real zagl = Compute_Z_AtCellCenter(i, j, k, z_phys_arr) - z_surf_arr(i,j,0);
1427 
1428  // linearly tapering profile -- following WRF, approximate top of
1429  // PBL as ustar * zscale
1430  cons_arr(i, j, k, RhoKE_comp) = rho * tke0 * std::max(
1431  (ust * zscale - zagl) / (std::max(ust, small) * zscale),
1432  small);
1433  });
1434  }
1435 }
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
Compute_Z_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Z_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:361
Here is the call graph for this function:

◆ lmask_min_reduce()

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

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

Parameters
[in]lmaskland-mask field
[in]nghostnumber of ghost cells included in the reduction
938  {
939  int lmask_min = amrex::ReduceMin(lmask, nghost, [=] AMREX_GPU_HOST_DEVICE(
940  amrex::Box const& bx, amrex::Array4<int const> const& lm_arr) -> int
941  {
942  int locmin = std::numeric_limits<int>::max();
943  const auto lo = lbound(bx);
944  const auto hi = ubound(bx);
945  for (int j = lo.y; j <= hi.y; ++j) {
946  for (int i = lo.x; i <= hi.x; ++i) {
947  locmin = std::min(locmin, lm_arr(i, j, 0));
948  }
949  }
950  return locmin;
951  });
952 
953  return lmask_min;
954  }

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

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

Parameters
[in]levlevel index
[in]nlevsnumber of AMR levels
[in]mfvconserved and velocity MultiFabs for this level
[in]Theta_primprimitive potential-temperature field
[in]Qv_primprimitive water-vapor field
[in]Qr_primprimitive rain-water field
[in]z_phys_ndnodal physical-height field
[in]Hwavewave-height field
[in]Lwavewavelength field
[in]eddyDiffseddy-diffusivity field
[in]lsm_dataland-surface-model data fields
[in]lsm_data_namenames for lsm_data entries
[in]lsm_fluxland-surface-model flux fields
[in]lsm_flux_namenames for lsm_flux entries
[in]sst_levsea-surface-temperature data by time
[in]tsk_levskin-temperature data by time
[in]lmask_levland-mask data by time
328  {
329  // Update MOST Average
330  m_ma.make_MOSTAverage_at_level(lev, mfv,
331  Theta_prim, Qv_prim, Qr_prim,
332  z_phys_nd);
333 
334  // Get CC vars
335  amrex::MultiFab& mf = *(mfv[0]);
336 
337  amrex::ParmParse pp("erf");
338 
339  // Do we have a time-varying surface roughness that needs to be saved?
340  if (lev == 0) {
341  const int nghost = 0; // ghost cells not included
342  int lmask_min = lmask_min_reduce(*lmask_lev[0].get(), nghost);
343  amrex::ParallelDescriptor::ReduceIntMin(lmask_min);
344 
345  m_var_z0 = (lmask_min < 1) & (rough_type_sea != RoughCalcType::CONSTANT);
346  if (m_var_z0) {
347  std::string rough_sea_string{"charnock"};
348  pp.query("most.roughness_type_sea", rough_sea_string);
349  amrex::Print() << "Variable sea roughness (type " << rough_sea_string
350  << ")" << std::endl;
351  }
352  }
353 
354  if (m_eddyDiffs_lev.size() < lev+1) {
355  m_Hwave_lev.resize(nlevs);
356  m_Lwave_lev.resize(nlevs);
357  m_eddyDiffs_lev.resize(nlevs);
358 
359  m_lsm_data_lev.resize(nlevs);
360  m_lsm_flux_lev.resize(nlevs);
361 
362  m_sst_lev.resize(nlevs);
363  m_tsk_lev.resize(nlevs);
364  m_lmask_lev.resize(nlevs);
365 
366  // Size the MOST params for all levels
367  z_0.resize(nlevs);
368  u_star.resize(nlevs);
369  w_star.resize(nlevs);
370  t_star.resize(nlevs);
371  q_star.resize(nlevs);
372  t_surf.resize(nlevs);
373  q_surf.resize(nlevs);
374  surface_diagnostic_source.resize(nlevs);
375  olen.resize(nlevs);
376  pblh.resize(nlevs);
377  }
378 
379  // Get pointers to SST,TSK and LANDMASK data
380  int nt_tot_sst = sst_lev.size();
381  m_sst_lev[lev].resize(nt_tot_sst);
382  for (int nt(0); nt < nt_tot_sst; ++nt) {
383  m_sst_lev[lev][nt] = sst_lev[nt].get();
384  }
385  int nt_tot_tsk = static_cast<int>(tsk_lev.size());
386  m_tsk_lev[lev].resize(nt_tot_tsk);
387  for (int nt(0); nt < nt_tot_tsk; ++nt) {
388  m_tsk_lev[lev][nt] = tsk_lev[nt].get();
389  }
390  int nt_tot_lmask = static_cast<int>(lmask_lev.size());
391  m_lmask_lev[lev].resize(nt_tot_lmask);
392  for (int nt(0); nt < nt_tot_lmask; ++nt) {
393  m_lmask_lev[lev][nt] = lmask_lev[nt].get();
394  }
395 
396  // Get pointers to wave data
397  m_Hwave_lev[lev] = Hwave;
398  m_Lwave_lev[lev] = Lwave;
399  m_eddyDiffs_lev[lev] = eddyDiffs;
400 
401  // Get pointers to LSM data and Fluxes
402  int ndata = static_cast<int>(lsm_data.size());
403  int nflux = static_cast<int>(lsm_flux.size());
404  m_lsm_data_name.resize(ndata);
405  m_lsm_data_lev[lev].resize(ndata);
406  m_lsm_flux_name.resize(nflux);
407  m_lsm_flux_lev[lev].resize(nflux);
408  for (int n(0); n < ndata; ++n) {
409  m_lsm_data_name[n] = lsm_data_name[n];
410  m_lsm_data_lev[lev][n] = lsm_data[n];
411  const std::string lc_name = amrex::toLower(lsm_data_name[n]);
412  if (lc_name == "theta" || lc_name == "t_surf") {
413  m_has_lsm_tsurf = true;
414  m_lsm_tsurf_indx = n;
415  m_has_ocean_lsm_tsurf = (lc_name == "t_surf");
416  }
417  }
418  int n_valid_lsm_flux = 0;
419  bool has_soil_t_flux = false;
420  for (int n(0); n < nflux; ++n) {
421  m_lsm_flux_name[n] = lsm_flux_name[n];
422  m_lsm_flux_lev[lev][n] = lsm_flux[n];
423  if (m_lsm_flux_lev[lev][n]) { ++n_valid_lsm_flux; }
424  if (amrex::toLower(m_lsm_flux_name[n]) == "soil_t_flux") {
425  has_soil_t_flux = true;
426  }
427  }
428  AMREX_ALWAYS_ASSERT((n_valid_lsm_flux==0 || n_valid_lsm_flux>=4 ||
429  (n_valid_lsm_flux==1 && has_soil_t_flux)));
430  if (n_valid_lsm_flux>=4) { m_has_lsm_fluxes = true; }
431 
432  // Check if there is a user-specified roughness file to be read
433  std::string fname;
434  bool read_z0 = false;
435  if ( (flux_type == FluxCalcType::MOENG) ||
436  (flux_type == FluxCalcType::ROTATE)) {
437  int count = pp.countval("most.roughness_file_name");
438  if (count > 1) {
439  AMREX_ALWAYS_ASSERT(count >= lev+1);
440  pp.query("most.roughness_file_name", fname, lev);
441  read_z0 = true;
442  } else if (count == 1) {
443  if (lev == 0) {
444  pp.query("most.roughness_file_name", fname);
445  } else {
446  // we will interpolate from the coarsest level
447  fname = "";
448  }
449  read_z0 = true;
450  }
451  // else use z0_const
452  }
453 
454  // Attributes for MFs and FABs
455  //--------------------------------------------------------
456  // Create a 2D ba for planar terrain, 3D for EB terrain
457  amrex::BoxArray ba = mf.boxArray();
458  amrex::BoxArray ba_flux;
459  amrex::IntVect ng{1,1,0};
460 
461  if (m_terrain_type == TerrainType::EB) {
462  // Use full 3D BoxArray for EB terrain
463  ba_flux = ba;
464  ng = amrex::IntVect{1,1,1}; // Include z ghost cells
465  } else {
466  // Collapse to 2D for planar terrain
467  amrex::BoxList bl2d = ba.boxList();
468  for (auto& b : bl2d) { b.setRange(2,0); }
469  ba_flux = amrex::BoxArray(std::move(bl2d));
470  }
471 
472  const amrex::DistributionMapping& dm = mf.DistributionMap();
473  const int ncomp = 1;
474 
475  // Z0 heights FAB
476  //--------------------------------------------------------
477  z_0[lev].define(ba_flux, dm, ncomp, ng);
478  z_0[lev].setVal(z0_const);
479  if (read_z0) {
480  read_custom_roughness(lev, fname);
481  }
482 
483  // 2D MFs for U*, T*, T_surf
484  //--------------------------------------------------------
485  u_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
486  u_star[lev]->setVal(bogus_large_value);
487 
488  w_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
489  w_star[lev]->setVal(bogus_large_value);
490 
491  t_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
492  t_star[lev]->setVal(zero); // default to neutral
493 
494  q_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
495  q_star[lev]->setVal(zero); // default to dry
496 
497  olen[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
498  olen[lev]->setVal(bogus_large_value);
499 
500  pblh[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
501  pblh[lev]->setVal(bogus_large_value);
502 
503  t_surf[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
504  t_surf[lev]->setVal(default_land_surf_temp);
505 
506  q_surf[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
507  q_surf[lev]->setVal(default_land_surf_moist);
508 
509  surface_diagnostic_source[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
510  surface_diagnostic_source[lev]->setVal(
511  surface_diagnostics::to_plot_value(surface_diagnostics::SurfaceDiagnosticSource::Missing));
512 
513  // TODO: Do we want an enum struct for indexing?
514 
515  bool use_sst = (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]);
516  bool use_tsk = (!m_tsk_lev[lev].empty() && m_tsk_lev[lev][0]);
517  if (use_sst || use_tsk || m_has_lsm_tsurf) {
518  // Valid SST, TSK or LSM data; t_surf set before computing fluxes (avoids
519  // extended lambda capture) Note that land temp will be set from m_tsk_lev
520  // while sea temp will be set from m_sst_lev
521  theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
522 
523  // Pathways in fill_tsurf_with_sst_and_tsk
524  amrex::Print() << "Using MOST with specified surface temperature ";
525  if (m_has_ocean_lsm_tsurf) {
526  amrex::Print() << "(OceanSurf: t_surf)" << std::endl;
527  } else if (m_has_lsm_tsurf && !use_sst && !use_tsk) {
528  amrex::Print() << "(SLM: theta)" << std::endl;
529  } else {
530  // NOTE: SST from the LOW file populates TSK in update_sst_tsk.
531  // So if we have TSK, it contains everything and has been
532  // sanity checked for valid SST values.
533  if (use_tsk) { m_ignore_sst = true; }
534  if (use_tsk) {
535  amrex::Print() << "(land: TSK, ";
536  } else {
537  amrex::Print() << "(land: T0, ";
538  }
539  if (use_tsk && !use_sst) {
540  amrex::Print() << "sea: TSK)" << std::endl;
541  } else {
542  amrex::Print() << "sea: SST)" << std::endl;
543  AMREX_ALWAYS_ASSERT(m_sst_lev[lev][0]);
544  }
545  }
546  }
547  }
bogus_large_value
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
MOSTAverage::make_MOSTAverage_at_level
void make_MOSTAverage_at_level(const int &lev, const amrex::Vector< amrex::MultiFab * > &vars_old, 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)
Definition: ERF_MOSTAverage.cpp:104
SurfaceLayer::lmask_min_reduce
int lmask_min_reduce(amrex::iMultiFab &lmask, const int &nghost)
Definition: ERF_SurfaceLayer.H:936
SurfaceLayer::m_has_lsm_tsurf
bool m_has_lsm_tsurf
Definition: ERF_SurfaceLayer.H:1049
SurfaceLayer::m_has_lsm_fluxes
bool m_has_lsm_fluxes
Definition: ERF_SurfaceLayer.H:1048
SurfaceLayer::read_custom_roughness
void read_custom_roughness(const int &lev, const std::string &fname)
Definition: ERF_SurfaceLayer.cpp:1445
Here is the call graph for this function:

◆ read_custom_roughness()

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

Read custom roughness data for one level.

Parameters
[in]levlevel index
[in]fnameroughness-data file name

Read or interpolate custom roughness length data.

Parameters
[in]levCurrent level
[in]fnameRoughness file name; an empty name interpolates from level 0
1447 {
1448  // Read the file if we have it
1449  if (!fname.empty()) {
1450  // Only the ioproc reads the file
1451  Gpu::HostVector<Real> m_x,m_y,m_z0;
1452  if (ParallelDescriptor::IOProcessor()) {
1453  Print()<<"Reading MOST roughness file at level " << lev << " : " << fname << std::endl;
1454  std::ifstream file(fname);
1455  Real value1,value2,value3;
1456  while(file>>value1>>value2>>value3){
1457  m_x.push_back(value1);
1458  m_y.push_back(value2);
1459  m_z0.push_back(value3);
1460  }
1461  file.close();
1462 
1463  AMREX_ALWAYS_ASSERT(m_x.size() == m_y.size());
1464  AMREX_ALWAYS_ASSERT(m_x.size() == m_z0.size());
1465  }
1466 
1467  // Broadcast the whole domain to every rank
1468  int ioproc = ParallelDescriptor::IOProcessorNumber();
1469  int nnode = static_cast<int>(m_x.size());
1470  ParallelDescriptor::Bcast(&nnode, 1, ioproc);
1471 
1472  if (!ParallelDescriptor::IOProcessor()) {
1473  m_x.resize(nnode);
1474  m_y.resize(nnode);
1475  m_z0.resize(nnode);
1476  }
1477  ParallelDescriptor::Bcast(m_x.data() , nnode, ioproc);
1478  ParallelDescriptor::Bcast(m_y.data() , nnode, ioproc);
1479  ParallelDescriptor::Bcast(m_z0.data(), nnode, ioproc);
1480 
1481  // Copy data to the GPU
1482  Gpu::DeviceVector<Real> d_x(nnode),d_y(nnode),d_z0(nnode);
1483  Gpu::copy(Gpu::hostToDevice, m_x.begin(), m_x.end(), d_x.begin());
1484  Gpu::copy(Gpu::hostToDevice, m_y.begin(), m_y.end(), d_y.begin());
1485  Gpu::copy(Gpu::hostToDevice, m_z0.begin(), m_z0.end(), d_z0.begin());
1486  Real* xp = d_x.data();
1487  Real* yp = d_y.data();
1488  Real* z0p = d_z0.data();
1489 
1490  // Each rank populates it's z_0[lev] MultiFab
1491  const int klo = m_geom[lev].Domain().smallEnd(2);
1492  for (MFIter mfi(z_0[lev]); mfi.isValid(); ++mfi)
1493  {
1494  Box gtbx = mfi.growntilebox();
1495 
1496  if (gtbx.smallEnd(2) != klo) { continue; }
1497 
1498  // Populate z_phys data
1499  Real tol = Real(1.0e-4);
1500  auto dx = m_geom[lev].CellSizeArray();
1501  auto ProbLoArr = m_geom[lev].ProbLoArray();
1502  int ilo = m_geom[lev].Domain().smallEnd(0);
1503  int jlo = m_geom[lev].Domain().smallEnd(1);
1504  int ihi = m_geom[lev].Domain().bigEnd(0);
1505  int jhi = m_geom[lev].Domain().bigEnd(1);
1506 
1507  Array4<Real> const& z0_arr = z_0[lev].array(mfi);
1508  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
1509  {
1510  // Clip indices for ghost-cells
1511  int ii = amrex::min(amrex::max(i,ilo),ihi);
1512  int jj = amrex::min(amrex::max(j,jlo),jhi);
1513 
1514  // Location of nodes
1515  Real x = ProbLoArr[0] + ii * dx[0];
1516  Real y = ProbLoArr[1] + jj * dx[1];
1517  int inode = ii + jj * (ihi-ilo+2); // stride is Nx+1
1518  if (std::sqrt(amrex::Math::powi<2>(x-xp[inode])+amrex::Math::powi<2>(y-yp[inode])) < tol) {
1519  z0_arr(i,j,klo) = z0p[inode];
1520  } else {
1521  // Unexpected list order, do brute force search
1522  Real z0loc = zero;
1523  bool found = false;
1524  for (int n=0; n<nnode; ++n) {
1525  Real delta=std::sqrt(amrex::Math::powi<2>(x-xp[n])+amrex::Math::powi<2>(y-yp[n]));
1526  if (delta < tol) {
1527  found = true;
1528  z0loc = z0p[n];
1529  break;
1530  }
1531  }
1532  AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
1533  amrex::ignore_unused(found);
1534  z0_arr(i,j,klo) = z0loc;
1535  }
1536  });
1537  } // mfi
1538  } else {
1539  AMREX_ALWAYS_ASSERT(lev > 0);
1540 
1541  Print()<<"Interpolating MOST roughness at level " << lev << std::endl;
1542 
1543  // Create a BC mapper that uses FOEXTRAP at domain bndry
1544  Vector<int> bc_lo(3,ERFBCType::foextrap);
1545  Vector<int> bc_hi(3,ERFBCType::foextrap);
1546  Vector<BCRec> bcr; bcr.push_back(BCRec(bc_lo.data(),bc_hi.data()));
1547 
1548  // Create ref ratio
1549  IntVect ratio;
1550  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1551  ratio[idim] = m_geom[lev].Domain().length(idim) / m_geom[0].Domain().length(idim);
1552  }
1553 
1554  // Create interp object and interpolate from the coarsest grid
1555  MFInterpolater* interp = &mf_cell_cons_interp;
1556  interp->interp(z_0[0] , 0,
1557  z_0[lev], 0,
1558  1, z_0[lev].nGrowVect(),
1559  m_geom[0], m_geom[lev],
1560  m_geom[lev].Domain(),ratio,
1561  bcr, 0);
1562  }
1563 }
m_y
@ m_y
Definition: ERF_DataStruct.H:28
m_x
@ m_x
Definition: ERF_DataStruct.H:27
Coord::y
@ y
Coord::x
@ x
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
ERFBCType::foextrap
@ foextrap
Definition: ERF_IndexDefines.H:248

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

Set the surface-moisture field to a constant value.

Parameters
[in]levlevel index
[in]qsurfsurface moisture
895 { q_surf[lev]->setVal(qsurf); }

◆ set_t_surf()

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

Set the surface-temperature field to a constant value.

Parameters
[in]levlevel index
[in]tsurfsurface temperature
880 { 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 
)

Update surface fluxes and related surface-layer state.

Parameters
[in]levlevel index
[in]elapsed_timecurrent elapsed simulation time
[in]elapsed_time_since_start_lowelapsed time relative to low-data start
[in,out]cons_inconserved state used by the flux update
[in]z_phys_ndnodal physical-height field
[in]walldistwall-distance field
[in]max_itersmaximum MOST iteration count

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

Parameters
[in]levCurrent level
[in]elapsed_timeCurrent simulation time
[in]elapsed_time_since_start_lowTime since the start of the lower-boundary data
[in,out]cons_inConserved state, updated when RANS TKE is initialized from surface-layer data
[in]z_phys_ndNodal physical height used by terrain-aware surface calculations
[in]walldistWall distance used when updating boundary TKE
[in]max_itersMaximum iterations to use in the MOST flux solve
25 {
26  // Update with SST/TSK data if we have a valid pointer
27  if (!m_has_ocean_lsm_tsurf &&
28  !m_sst_lev[lev].empty() && m_sst_lev[lev][0]) {
29  fill_tsurf_with_sst_and_tsk(lev, elapsed_time_since_start_low);
30  }
31 
32  // Apply heating rate if needed
33  if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE &&
34  !m_has_ocean_lsm_tsurf) {
35  update_surf_temp(elapsed_time_since_start_low);
36  }
37 
38  // Update qsurf with qsat over sea
39  if (use_moisture) {
40  fill_qsurf_with_qsat(lev, cons_in, z_phys_nd);
41  }
42 
43  // Update land surface temp if we have a valid pointer
44  if (m_has_lsm_tsurf) { get_lsm_tsurf(lev); }
45 
46  // Fill interior ghost cells
47  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
48 
49  // Compute plane averages for all vars (regardless of flux type)
50  m_ma.compute_averages(lev);
51 
52  // NOTE: Do iterations to seed variables on the first step (LSM called post step)
53  // as well as compute values where invalid LSM fluxes may reside
54  //*******************************************************************************
55  // ***************************************************************
56  // Iterate the fluxes if moeng type
57  // First iterate over land -- the only model for surface roughness
58  // over land is RoughCalcType::CONSTANT
59  // ***************************************************************
60  if (flux_type == FluxCalcType::MOENG ||
61  flux_type == FluxCalcType::ROTATE) {
62  bool is_land = true;
63  // Do we have a constant flux for moisture over land?
64  bool cons_qflux = ( (moist_type == MoistCalcType::MOISTURE_FLUX) ||
65  (moist_type == MoistCalcType::ADIABATIC) );
66  if (m_terrain_type != TerrainType::EB) {
67  if (theta_type == ThetaCalcType::HEAT_FLUX) {
68  if (rough_type_land == RoughCalcType::CONSTANT) {
69  surface_flux most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
70  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
71  } else {
72  amrex::Abort("Unknown value for rough_type_land");
73  }
74  } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
75  if (rough_type_land == RoughCalcType::CONSTANT) {
76  surface_temp most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
77  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
78  } else {
79  amrex::Abort("Unknown value for rough_type_land");
80  }
81  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
82  (moist_type == MoistCalcType::ADIABATIC)) {
83  if (rough_type_land == RoughCalcType::CONSTANT) {
84  adiabatic most_flux(surf_temp_flux, surf_moist_flux);
85  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
86  } else {
87  amrex::Abort("Unknown value for rough_type_land");
88  }
89  } else {
90  amrex::Abort("Unknown value for theta_type");
91  }
92  // EB
93  } else {
94  if (theta_type == ThetaCalcType::HEAT_FLUX) {
95  if (rough_type_land == RoughCalcType::CONSTANT) {
96  surface_flux_eb most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
97  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
98  } else {
99  amrex::Abort("Unknown value for rough_type_land");
100  }
101  } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
102  if (rough_type_land == RoughCalcType::CONSTANT) {
103  surface_temp_eb most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
104  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
105  } else {
106  amrex::Abort("Unknown value for rough_type_land");
107  }
108  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
109  (moist_type == MoistCalcType::ADIABATIC)) {
110  if (rough_type_land == RoughCalcType::CONSTANT) {
111  adiabatic_eb most_flux(surf_temp_flux, surf_moist_flux);
112  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
113  } else {
114  amrex::Abort("Unknown value for rough_type_land");
115  }
116  } else {
117  amrex::Abort("Unknown value for theta_type");
118  }
119  } // EB
120  } // MOENG -- LAND
121 
122  // Update u*/T*/q*/L over land (iterations or from LSM fluxes)
123  if (m_has_lsm_fluxes && elapsed_time > zero) {
124  compute_sfc_params_from_lsm_fluxes(lev, cons_in);
125  }
126 
127  // ***************************************************************
128  // Iterate the fluxes if moeng type
129  // Next iterate over sea -- the models for surface roughness
130  // over sea are CHARNOCK, DONELAN, MODIFIED_CHARNOCK or WAVE_COUPLED
131  // NOTE: Sea surface fluxes are not supported for EB terrain
132  // ***************************************************************
133  if ((flux_type == FluxCalcType::MOENG ||
134  flux_type == FluxCalcType::ROTATE) &&
135  m_terrain_type != TerrainType::EB) {
136  bool is_land = false;
137  // NOTE: Do not allow default to adiabatic over sea (we have Qvs at surface)
138  // Do we have a constant flux for moisture over sea?
139  bool cons_qflux = (moist_type == MoistCalcType::MOISTURE_FLUX);
140  if (theta_type == ThetaCalcType::HEAT_FLUX) {
141  if (rough_type_sea == RoughCalcType::CHARNOCK) {
142  surface_flux_charnock most_flux(surf_temp_flux, surf_moist_flux,
143  cnk_a, cnk_visc, cons_qflux);
144  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
145  } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
146  surface_flux_mod_charnock most_flux(surf_temp_flux, surf_moist_flux,
147  depth, cons_qflux);
148  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
149  } else if (rough_type_sea == RoughCalcType::DONELAN) {
150  surface_flux_donelan most_flux(surf_temp_flux, surf_moist_flux,
151  cons_qflux);
152  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
153  } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
154  surface_flux_wave_coupled most_flux(surf_temp_flux, surf_moist_flux,
155  cons_qflux);
156  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
157  } else {
158  amrex::Abort("Unknown value for rough_type_sea");
159  }
160 
161  } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
162  if (rough_type_sea == RoughCalcType::CHARNOCK) {
163  surface_temp_charnock most_flux(surf_temp_flux, surf_moist_flux,
164  cnk_a, cnk_visc, cons_qflux);
165  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
166  } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
167  surface_temp_mod_charnock most_flux(surf_temp_flux, surf_moist_flux,
168  depth, cons_qflux);
169  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
170  } else if (rough_type_sea == RoughCalcType::DONELAN) {
171  surface_temp_donelan most_flux(surf_temp_flux, surf_moist_flux,
172  cons_qflux);
173  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
174  } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
175  surface_temp_wave_coupled most_flux(surf_temp_flux, surf_moist_flux,
176  cons_qflux);
177  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
178  } else {
179  amrex::Abort("Unknown value for rough_type_sea");
180  }
181 
182  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
183  (moist_type == MoistCalcType::ADIABATIC)) {
184  if (rough_type_sea == RoughCalcType::CHARNOCK) {
185  adiabatic_charnock most_flux(surf_temp_flux, surf_moist_flux,
186  cnk_a, cnk_visc);
187  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
188  } else if (rough_type_sea == RoughCalcType::MODIFIED_CHARNOCK) {
189  adiabatic_mod_charnock most_flux(surf_temp_flux, surf_moist_flux,
190  depth);
191  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
192  } else if (rough_type_sea == RoughCalcType::DONELAN) {
193  adiabatic_donelan most_flux(surf_temp_flux, surf_moist_flux);
194  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
195  } else if (rough_type_sea == RoughCalcType::WAVE_COUPLED) {
196  adiabatic_wave_coupled most_flux(surf_temp_flux, surf_moist_flux);
197  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
198  } else {
199  amrex::Abort("Unknown value for rough_type_sea");
200  }
201  } else {
202  amrex::Abort("Unknown value for theta_type");
203  }
204  } // MOENG -- SEA
205 
206  if (flux_type == FluxCalcType::CUSTOM || flux_type == FluxCalcType::RICO) {
207  if (custom_rhosurf > 0) {
208  specified_rho_surf = true;
209  u_star[lev]->setVal(std::sqrt(custom_rhosurf) * custom_ustar);
210  t_star[lev]->setVal(custom_rhosurf * custom_tstar);
211  q_star[lev]->setVal(custom_rhosurf * custom_qstar);
212  } else {
213  u_star[lev]->setVal(custom_ustar);
214  t_star[lev]->setVal(custom_tstar);
215  q_star[lev]->setVal(custom_qstar);
216  }
217  }
218 
219  if (m_update_k_rans) {
220  const bool use_ref_theta = (theta_ref > 0);
221  const Real l_inv_theta0 = (use_ref_theta) ? one / theta_ref : one;
222  const Real l_inv_Cmu2 = inv_Cmu2;
223  const int klo = m_geom[lev].Domain().smallEnd(2);
224  IntVect ng = u_star[lev]->nGrowVect(); ng[2] = 0;
225 
226  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi)
227  {
228  Box gpbx = mfi.tilebox(IntVect(0),ng);
229 
230  if (gpbx.smallEnd(2) != klo) { continue; }
231 
232  gpbx.makeSlab(2,klo);
233 
234  auto cons_arr = cons_in.array(mfi);
235  const auto& u_star_arr = u_star[lev]->const_array(mfi);
236  const auto& t_star_arr = t_star[lev]->const_array(mfi);
237  const auto& dist_arr = walldist->const_array(mfi);
238 
239  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
240  {
241  Real rho = cons_arr(i,j,k,Rho_comp);
242  if (t_star_arr(i,j,0) < -1e-8) {
243  // Only destabilizing buoyancy flux affects the boundary k
244  // tstar < 0 ==> B > 0
245  Real B = -CONST_GRAV * l_inv_theta0 * u_star_arr(i,j,0) * t_star_arr(i,j,0);
246  if (!use_ref_theta) {
247  B *= cons_arr(i,j,k,Rho_comp) /
248  cons_arr(i,j,k,RhoTheta_comp);
249  }
250 
251  // Axell & Liungman 2001, Eqn. 16
252  cons_arr(i,j,k,RhoKE_comp) = rho * l_inv_Cmu2 *
253  std::pow(
254  u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0)
255  + KAPPA * B * dist_arr(i,j,k),
256  two/three);
257  } else {
258  cons_arr(i,j,k,RhoKE_comp) = rho * l_inv_Cmu2 * u_star_arr(i,j,0) * u_star_arr(i,j,0);
259  }
260  });
261  }
262  }
263 
264  u_star[lev]->FillBoundary(m_geom[lev].periodicity());
265  t_star[lev]->FillBoundary(m_geom[lev].periodicity());
266  q_star[lev]->FillBoundary(m_geom[lev].periodicity());
267  olen[lev]->FillBoundary(m_geom[lev].periodicity());
268 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
MOSTAverage::compute_averages
void compute_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:925
SurfaceLayer::get_lsm_tsurf
void get_lsm_tsurf(const int &lev)
Definition: ERF_SurfaceLayer.cpp:1279
SurfaceLayer::fill_qsurf_with_qsat
void fill_qsurf_with_qsat(const int &lev, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd)
Definition: ERF_SurfaceLayer.cpp:1231
SurfaceLayer::compute_fluxes
void compute_fluxes(const int &lev, const int &max_iters, amrex::MultiFab &cons_in, const FluxIter &most_flux, bool is_land)
SurfaceLayer::update_surf_temp
void update_surf_temp(const double &time)
Definition: ERF_SurfaceLayer.H:781
SurfaceLayer::fill_tsurf_with_sst_and_tsk
void fill_tsurf_with_sst_and_tsk(const int &lev, const double &time)
Definition: ERF_SurfaceLayer.cpp:1137
SurfaceLayer::compute_sfc_params_from_lsm_fluxes
void compute_sfc_params_from_lsm_fluxes(const int &lev, amrex::MultiFab &cons_in)
Definition: ERF_SurfaceLayer.cpp:1064
adiabatic_charnock
Definition: ERF_MOSTStress.H:226
adiabatic_donelan
Definition: ERF_MOSTStress.H:437
adiabatic_eb
EB surface-layer model for adiabatic constant-roughness fluxes.
Definition: ERF_EBMOSTStress.H:14
adiabatic_mod_charnock
Definition: ERF_MOSTStress.H:336
adiabatic_wave_coupled
Definition: ERF_MOSTStress.H:522
adiabatic
Definition: ERF_MOSTStress.H:162
surface_flux_charnock
Definition: ERF_MOSTStress.H:738
surface_flux_donelan
Definition: ERF_MOSTStress.H:1003
surface_flux_eb
EB surface-layer model with prescribed surface fluxes and constant roughness.
Definition: ERF_EBMOSTStress.H:207
surface_flux_mod_charnock
Definition: ERF_MOSTStress.H:875
surface_flux_wave_coupled
Definition: ERF_MOSTStress.H:1127
surface_flux
Definition: ERF_MOSTStress.H:617
surface_temp_charnock
Definition: ERF_MOSTStress.H:1414
surface_temp_donelan
Definition: ERF_MOSTStress.H:1785
surface_temp_eb
EB surface-layer model with prescribed surface temperature and constant roughness.
Definition: ERF_EBMOSTStress.H:67
surface_temp_mod_charnock
Definition: ERF_MOSTStress.H:1604
surface_temp_wave_coupled
Definition: ERF_MOSTStress.H:1919
surface_temp
Definition: ERF_MOSTStress.H:1264
Here is the call graph for this function:

◆ update_mac_ptrs()

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

Update MOST-average field pointers.

Parameters
[in]levlevel index
[in]vars_oldold-time state variables
[in]Theta_primprimitive potential-temperature fields by level
[in]Qv_primprimitive water-vapor fields by level
[in]Qr_primprimitive rain-water fields by level
810  {
811  m_ma.update_field_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
812  }
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:293
Here is the call graph for this function:

◆ update_pblh()

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

Wrapper around compute_pblh.

Parameters
[in]levlevel index
[in,out]varsstate variables used by the PBL-height calculation
[in]z_phys_cccell-centered physical-height field
[in]moisture_indicesindices for moisture components

Update PBL height using the configured estimator.

Parameters
[in]levCurrent level
[in]varsLevel-indexed state MultiFabs passed to the PBL height estimator
[in]z_phys_ccCell-centered physical height used by the PBL height estimator
[in]moisture_indicesMoisture component indices used by the PBL height estimator
1329 {
1330  if (pblh_type == PBLHeightCalcType::MYNN25) {
1331  MYNNPBLH estimator;
1332  compute_pblh(lev, vars, z_phys_cc, estimator, moisture_indices);
1333  } else if (pblh_type == PBLHeightCalcType::YSU || pblh_type == PBLHeightCalcType::MRF) {
1334  amrex::Error("YSU/MRF PBLH calc not implemented yet");
1335  }
1336 }
SurfaceLayer::compute_pblh
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)
MYNNPBLH
Definition: ERF_PBLHeight.H:8

◆ update_sst_ptr()

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

Update one stored sea-surface-temperature pointer.

Parameters
[in]levlevel index
[in]itimetime-slice index
[in]sst_ptrsea-surface-temperature field pointer
963  {
964  m_sst_lev[lev][itime] = sst_ptr;
965  }

◆ update_surf_temp()

void SurfaceLayer::update_surf_temp ( const double &  time)
inline

Update prescribed surface temperature from the configured heating rate.

Parameters
[in]timeelapsed simulation time
782  {
783  if (m_has_ocean_lsm_tsurf) {
784  return;
785  }
786  if (surf_heating_rate != 0) {
787  int nlevs = static_cast<int>(m_geom.size());
788  for (int lev = 0; lev < nlevs; lev++) {
789  t_surf[lev]->setVal(surf_temp + surf_heating_rate * static_cast<amrex::Real>(time));
790  amrex::Print() << "Surface temp at t=" << time << ": "
791  << surf_temp + surf_heating_rate * time << std::endl;
792  }
793  }
794  }

◆ update_tsk_ptr()

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

Update one stored skin-temperature pointer.

Parameters
[in]levlevel index
[in]itimetime-slice index
[in]tsk_ptrskin-temperature field pointer
974  {
975  m_tsk_lev[lev][itime] = tsk_ptr;
976  }

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: