ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
SurfaceLayer Class Reference

#include <ERF_SurfaceLayer.H>

Collaboration diagram for SurfaceLayer:

Public Types

enum class  FluxCalcType {
  MOENG = 0 , DONELAN , CUSTOM , BULK_COEFF ,
  ROTATE , RICO
}
 
enum class  ThetaCalcType { ADIABATIC = 0 , HEAT_FLUX , SURFACE_TEMPERATURE }
 
enum class  MoistCalcType { ADIABATIC = 0 , MOISTURE_FLUX , SURFACE_MOISTURE }
 
enum class  RoughCalcType {
  CONSTANT = 0 , CHARNOCK , MODIFIED_CHARNOCK , DONELAN ,
  WAVE_COUPLED
}
 
enum class  PBLHeightCalcType { None , MYNN25 , YSU , MRF }
 

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, amrex::Real start_low_time, amrex::Real final_low_time, amrex::Real low_time_interval=zero)
 
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 amrex::Real &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< 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 eb_ &ebfact)
 
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< 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 eb_ &ebfact, const FluxCalc &flux_comp)
 
void fill_tsurf_with_sst_and_tsk (const int &lev, const amrex::Real &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 amrex::Real &time)
 
void update_mac_ptrs (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_old, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Theta_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qv_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qr_prim)
 
amrex::MultiFab * get_u_star (const int &lev)
 
amrex::MultiFab * get_w_star (const int &lev)
 
amrex::MultiFab * get_t_star (const int &lev)
 
amrex::MultiFab * get_q_star (const int &lev)
 
amrex::MultiFab * get_olen (const int &lev)
 
amrex::MultiFab * get_pblh (const int &lev)
 
const amrex::MultiFab * get_mac_avg (const int &lev, int comp)
 
amrex::MultiFab * get_t_surf (const int &lev)
 
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::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< 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, [[maybe_unused]] const eb_ &ebfact, 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
 
amrex::Real m_start_low_time
 
amrex::Real m_final_low_time
 
amrex::Real 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_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
 
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< 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

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 
646  {
647  MOENG = 0, ///< Moeng functional form
648  DONELAN, ///< Donelan functional form
649  CUSTOM, ///< Custom constant flux functional form
650  BULK_COEFF, ///< Bulk transfer coefficient functional form
651  ROTATE, ///< Terrain rotation flux functional form
652  RICO
653  };

◆ MoistCalcType

Enumerator
ADIABATIC 
MOISTURE_FLUX 

Qv-flux specified.

SURFACE_MOISTURE 

Surface Qv specified.

661  {
662  ADIABATIC = 0,
663  MOISTURE_FLUX, ///< Qv-flux specified
664  SURFACE_MOISTURE ///< Surface Qv specified
665  };

◆ PBLHeightCalcType

Enumerator
None 
MYNN25 
YSU 
MRF 
675 { None, MYNN25, YSU, MRF };

◆ RoughCalcType

Enumerator
CONSTANT 

Constant z0.

CHARNOCK 
MODIFIED_CHARNOCK 
DONELAN 
WAVE_COUPLED 
667  {
668  CONSTANT = 0, ///< Constant z0
669  CHARNOCK,
670  MODIFIED_CHARNOCK,
671  DONELAN,
672  WAVE_COUPLED
673  };

◆ ThetaCalcType

Enumerator
ADIABATIC 
HEAT_FLUX 

Heat-flux specified.

SURFACE_TEMPERATURE 

Surface temperature specified.

655  {
656  ADIABATIC = 0,
657  HEAT_FLUX, ///< Heat-flux specified
658  SURFACE_TEMPERATURE ///< Surface temperature specified
659  };

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,
amrex::Real  start_low_time,
amrex::Real  final_low_time,
amrex::Real  low_time_interval = zero 
)
inlineexplicit
47  : m_geom(geom),
48  m_rotate(use_rot_surface_flux),
49  m_start_low_time(start_low_time),
50  m_final_low_time(final_low_time),
51  m_low_time_interval(low_time_interval),
52  m_ma(geom, (z_phys_nd[0] != nullptr), a_pp_prefix, a_mesh_type, a_terrain_type)
53  {
54  // We have a moisture model if Qv_prim is a valid pointer
55  use_moisture = (Qv_prim[0].get());
56 
57  // Get roughness
58  amrex::ParmParse pp("erf");
59  pp.query("most.z0", z0_const);
60 
61  // Specify how to compute the flux
62  if (use_rot_surface_flux) {
64  } else {
65  std::string flux_string_in;
66  std::string flux_string{"moeng"};
67  auto read_flux = pp.query("surface_layer.flux_type", flux_string_in);
68  if (read_flux) {
69  flux_string = amrex::toLower(flux_string_in);
70  }
71  if (flux_string == "donelan") {
73  } else if (flux_string == "moeng") {
75  } else if (flux_string == "rico") {
77  } else if (flux_string == "bulk_coeff") {
79  } else if (flux_string == "custom") {
81  } else {
82  amrex::Abort("Undefined MOST flux type!");
83  }
84  }
85 
86  // Include w* to handle free convection (Beljaars 1995, QJRMS)
87  pp.query("most.include_wstar", m_include_wstar);
88 
89  std::string pblh_string_in;
90  std::string pblh_string{"none"};
91  auto read_pblh = pp.query("most.pblh_calc", pblh_string_in);
92  if (read_pblh) {
93  pblh_string = amrex::toLower(pblh_string_in);
94  }
95  if (pblh_string == "none") {
97  } else if (pblh_string == "mynn25") {
99  } else if (pblh_string == "mynnedmf") {
101  } else if (pblh_string == "ysu") {
103  } else if (pblh_string == "mrf") {
105  } else {
106  amrex::Abort("Undefined PBLH calc type!");
107  }
108 
109  // Get surface temperature
110  auto erf_st = pp.query("most.surf_temp", surf_temp);
111  if (erf_st) { default_land_surf_temp = surf_temp; }
112 
113  // Get surface moisture
114  bool erf_sq = false;
115  if (use_moisture) { erf_sq = pp.query("most.surf_moist", surf_moist); }
116  if (erf_sq) { default_land_surf_moist = surf_moist; }
117 
118  // Custom type user must specify the fluxes
123  pp.get("most.ustar", custom_ustar);
124  pp.get("most.tstar", custom_tstar);
125  pp.get("most.qstar", custom_qstar);
126  pp.query("most.rhosurf", custom_rhosurf);
127  if (custom_qstar != 0) {
129  "Specified custom MOST qv flux without moisture model!");
130  }
131  amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
132  << custom_ustar << " " << custom_tstar << " "
133  << custom_qstar << std::endl;
134 
135  // Bulk transfer coefficient (must specify coeffs and surface values)
136  } else if (flux_type == FluxCalcType::BULK_COEFF) {
137  pp.get("most.Cd", m_Cd);
138  pp.get("most.Ch", m_Ch);
139  pp.get("most.Cq", m_Cq);
140  pp.get("most.surf_temp", default_land_surf_temp);
141  pp.get("most.surf_moist", default_land_surf_moist);
142  amrex::Print() << "Using specified Cd, Ch, Cq for MOST = "
143  << m_Cd << " " << m_Ch << " "
144  << m_Cq << std::endl;
145 
146  // Specify surface temperature/moisture or surface flux
147  } else {
148  if (erf_st) {
150  pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
151  surf_heating_rate = surf_heating_rate / amrex::Real(3600.0); // [K/s]
152  if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
153  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
154  }
155  } else {
156  pp.query("most.surf_temp_flux", surf_temp_flux);
157 
158  if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
159  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
160  }
161  if (std::abs(surf_temp_flux) >
164  } else {
166  }
167  }
168 
169  if (erf_sq) {
171  } else {
172  pp.query("most.surf_moist_flux", surf_moist_flux);
173  if (std::abs(surf_moist_flux) >
176  } else {
178  }
179  }
180  }
181 
183  {
184  pp.query("most.rico.theta_z0", rico_theta_z0);
185  pp.query("most.rico.qsat_z0", rico_qsat_z0);
186  }
187 
188  // Make sure the inputs file doesn't try to use most.roughness_type
189  std::string bogus_input;
190  if (pp.query("most.roughness_type", bogus_input) > 0) {
191  amrex::Abort("most.roughness_type is deprecated; use "
192  "most.roughness_type_land and/or most.roughness_type_sea");
193  }
194 
195  // Specify how to compute the surface flux over land (if there is any)
196  std::string rough_land_string_in;
197  std::string rough_land_string{"constant"};
198  auto read_rough_land =
199  pp.query("most.roughness_type_land", rough_land_string_in);
200  if (read_rough_land) {
201  rough_land_string = amrex::toLower(rough_land_string_in);
202  }
203  if (rough_land_string == "constant") {
205  } else {
206  amrex::Abort("Undefined MOST roughness type for land!");
207  }
208 
209  // Specify how to compute the surface flux over sea (if there is any)
210  std::string rough_sea_string_in;
211  std::string rough_sea_string{"charnock"};
212  auto read_rough_sea = pp.query("most.roughness_type_sea", rough_sea_string_in);
213  if (read_rough_sea) {
214  rough_sea_string = amrex::toLower(rough_sea_string_in);
215  }
216  if (rough_sea_string == "charnock") {
218  pp.query("most.charnock_constant", cnk_a);
219  pp.query("most.charnock_viscosity", cnk_visc);
220  if (cnk_a > 0) {
221  amrex::Print() << "If there is water, Charnock relation with C_a="
222  << cnk_a << (cnk_visc ? " and viscosity" : "")
223  << " will be used" << std::endl;
224  } else {
225  amrex::Print() << "If there is water, Charnock relation with variable "
226  "Charnock parameter (COARE3.0)"
227  << (cnk_visc ? " and viscosity" : "") << " will be used"
228  << std::endl;
229  }
230  } else if (rough_sea_string == "coare3.0") {
232  amrex::Print() << "If there is water, Charnock relation with variable "
233  "Charnock parameter (COARE3.0)"
234  << (cnk_visc ? " and viscosity" : "") << " will be used"
235  << std::endl;
236  cnk_a = -1;
237  } else if (rough_sea_string == "donelan") {
239  } else if (rough_sea_string == "modified_charnock") {
241  pp.query("most.modified_charnock_depth", depth);
242  } else if (rough_sea_string == "wave_coupled") {
244  } else if (rough_sea_string == "constant") {
246  } else {
247  amrex::Abort("Undefined MOST roughness type for sea!");
248  }
249 
250  // use skin temperature instead of sea-surface temperature
251  // (wrfinput data may have lower resolution SST data)
252  pp.query("most.ignore_sst", m_ignore_sst);
253 
254  // If we're using the RANS k model, then we need to update the dirichlet
255  // BC based on the instantaneous u* and θ*; the turbulence modeling
256  // choices can vary per level but for now, assume that if specified then
257  // all levels are using the same RANS model.
258  m_update_k_rans = (a_turb_choice.rans_type == RANSType::kEqn &&
259  a_turb_choice.dirichlet_k == true);
260  if (m_update_k_rans) {
261  inv_Cmu2 = one / (a_turb_choice.Cmu0 * a_turb_choice.Cmu0);
262  theta_ref = a_turb_choice.theta_ref;
263  }
264 
265  } // constructor
constexpr amrex::Real one
Definition: ERF_Constants.H:7
ParmParse pp("prob")
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")
ThetaCalcType theta_type
Definition: ERF_SurfaceLayer.H:678
bool m_include_wstar
Definition: ERF_SurfaceLayer.H:692
bool m_rotate
Definition: ERF_SurfaceLayer.H:687
PBLHeightCalcType pblh_type
Definition: ERF_SurfaceLayer.H:682
bool use_moisture
Definition: ERF_SurfaceLayer.H:715
amrex::Real m_Cq
Definition: ERF_SurfaceLayer.H:721
RoughCalcType rough_type_land
Definition: ERF_SurfaceLayer.H:680
amrex::Real z0_const
Definition: ERF_SurfaceLayer.H:693
amrex::Real cnk_a
Definition: ERF_SurfaceLayer.H:706
amrex::Real m_Ch
Definition: ERF_SurfaceLayer.H:720
amrex::Real m_start_low_time
Definition: ERF_SurfaceLayer.H:688
amrex::Real surf_temp
Definition: ERF_SurfaceLayer.H:695
amrex::Real rico_qsat_z0
Definition: ERF_SurfaceLayer.H:713
bool m_update_k_rans
Definition: ERF_SurfaceLayer.H:746
amrex::Real surf_moist_flux
Definition: ERF_SurfaceLayer.H:700
RoughCalcType rough_type_sea
Definition: ERF_SurfaceLayer.H:681
amrex::Real surf_moist
Definition: ERF_SurfaceLayer.H:699
bool m_ignore_sst
Definition: ERF_SurfaceLayer.H:723
amrex::Real custom_qstar
Definition: ERF_SurfaceLayer.H:703
amrex::Real custom_rhosurf
Definition: ERF_SurfaceLayer.H:704
@ MOENG
Moeng functional form.
@ BULK_COEFF
Bulk transfer coefficient functional form.
@ CUSTOM
Custom constant flux functional form.
@ ROTATE
Terrain rotation flux functional form.
@ DONELAN
Donelan functional form.
@ SURFACE_MOISTURE
Surface Qv specified.
@ MOISTURE_FLUX
Qv-flux specified.
amrex::Real depth
Definition: ERF_SurfaceLayer.H:708
amrex::Real default_land_surf_moist
Definition: ERF_SurfaceLayer.H:698
amrex::Real rico_theta_z0
Definition: ERF_SurfaceLayer.H:712
amrex::Real surf_temp_flux
Definition: ERF_SurfaceLayer.H:697
amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_SurfaceLayer.H:686
amrex::Real theta_ref
Definition: ERF_SurfaceLayer.H:748
amrex::Real custom_tstar
Definition: ERF_SurfaceLayer.H:702
amrex::Real m_final_low_time
Definition: ERF_SurfaceLayer.H:689
bool cnk_visc
Definition: ERF_SurfaceLayer.H:707
amrex::Real surf_heating_rate
Definition: ERF_SurfaceLayer.H:696
FluxCalcType flux_type
Definition: ERF_SurfaceLayer.H:677
MoistCalcType moist_type
Definition: ERF_SurfaceLayer.H:679
amrex::Real inv_Cmu2
Definition: ERF_SurfaceLayer.H:747
amrex::Real custom_ustar
Definition: ERF_SurfaceLayer.H:701
amrex::Real m_Cd
Definition: ERF_SurfaceLayer.H:719
amrex::Real default_land_surf_temp
Definition: ERF_SurfaceLayer.H:694
@ SURFACE_TEMPERATURE
Surface temperature specified.
@ HEAT_FLUX
Heat-flux specified.
amrex::Real m_low_time_interval
Definition: ERF_SurfaceLayer.H:690
MOSTAverage m_ma
Definition: ERF_SurfaceLayer.H:725
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
RANSType rans_type
Definition: ERF_TurbStruct.H:413
bool dirichlet_k
Definition: ERF_TurbStruct.H:415
amrex::Real Cmu0
Definition: ERF_TurbStruct.H:391
amrex::Real theta_ref
Definition: ERF_TurbStruct.H:402
Here is the call graph for this function:

Member Function Documentation

◆ compute_fluxes() [1/2]

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

◆ compute_fluxes() [2/2]

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

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

Parameters
[in]levCurrent level
[in]max_itersmaximum iterations to use
[in]most_fluxstructure to iteratively compute ustar and tstar
233 {
234  // Pointers to the computed averages
235  const auto *const tm_ptr = m_ma.get_average(lev,2); // potential temperature
236  const auto *const qvm_ptr = m_ma.get_average(lev,3); // water vapor mixing ratio
237  const auto *const tvm_ptr = m_ma.get_average(lev,4); // virtual potential temperature
238  const auto *const umm_ptr = m_ma.get_average(lev,5); // horizontal velocity magnitude
239  const auto *const zref_ptr = m_ma.get_zref(lev); // reference height
240 
241  const int klo = m_geom[lev].Domain().smallEnd(2);
242  IntVect ng = u_star[lev]->nGrowVect(); ng[2] = 0;
243 
244  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi)
245  {
246  Box gtbx = mfi.tilebox(IntVect(0),ng);
247 
248  if (gtbx.smallEnd(2) != klo) { continue; }
249 
250  gtbx.makeSlab(2,klo);
251 
252  auto u_star_arr = u_star[lev]->array(mfi);
253  auto t_star_arr = t_star[lev]->array(mfi);
254  auto q_star_arr = q_star[lev]->array(mfi);
255  auto t_surf_arr = t_surf[lev]->array(mfi);
256  auto q_surf_arr = q_surf[lev]->array(mfi);
257  auto olen_arr = olen[lev]->array(mfi);
258 
259  const auto tm_arr = tm_ptr->array(mfi);
260  const auto tvm_arr = tvm_ptr->array(mfi);
261  const auto qvm_arr = qvm_ptr->array(mfi);
262  const auto umm_arr = umm_ptr->array(mfi);
263  const auto zref_arr = zref_ptr->array(mfi);
264  const auto z0_arr = z_0[lev].array(mfi);
265 
266  // PBL height if we need to calculate wstar for the Beljaars correction
267  // TODO: can/should we apply this in LES mode?
268  const auto w_star_arr = (m_include_wstar) ? w_star[lev].get()->array(mfi) : Array4<Real> {};
269  const auto pblh_arr = (m_include_wstar) ? pblh[lev].get()->array(mfi) : Array4<Real> {};
270 
271  // Wave properties if they exist
272  const auto Hwave_arr = (m_Hwave_lev[lev]) ? m_Hwave_lev[lev]->array(mfi) : Array4<Real> {};
273  const auto Lwave_arr = (m_Lwave_lev[lev]) ? m_Lwave_lev[lev]->array(mfi) : Array4<Real> {};
274  const auto eta_arr = (m_eddyDiffs_lev[lev]) ? m_eddyDiffs_lev[lev]->array(mfi) : Array4<Real> {};
275 
276  // Land mask array if it exists
277  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
278  Array4<int> {};
279 
280  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
281  {
282  if (( is_land && lmask_arr(i,j,0) == 1) ||
283  (!is_land && lmask_arr(i,j,0) == 0))
284  {
285  // NOTE: All 2D MFs so k index is always 0 from ba2d definition
286  most_flux.iterate_flux(i, j, 0, max_iters,
287  zref_arr, // set in most average
288  z0_arr, // updated if(!is_land)
289  umm_arr, tm_arr, tvm_arr, qvm_arr,
290  u_star_arr, // updated
291  w_star_arr, // updated if(m_include_wstar)
292  t_star_arr, q_star_arr, // updated
293  t_surf_arr, q_surf_arr, olen_arr, // updated
294  pblh_arr, // updated if(m_include_wstar)
295  Hwave_arr, Lwave_arr, eta_arr);
296  }
297  });
298  }
299 }
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+myhalf) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=one) { Real dT=T_pert *(std::cos(PI *L)+one)/two;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
pp get("wavelength", wavelength)
amrex::MultiFab * get_zref(const int &lev) const
Definition: ERF_MOSTAverage.H:108
const amrex::MultiFab * get_average(const int &lev, const int &comp) const
Definition: ERF_MOSTAverage.H:105
amrex::Vector< amrex::Vector< amrex::iMultiFab * > > m_lmask_lev
Definition: ERF_SurfaceLayer.H:737
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_surf
Definition: ERF_SurfaceLayer.H:732
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_star
Definition: ERF_SurfaceLayer.H:729
amrex::Vector< amrex::MultiFab * > m_Lwave_lev
Definition: ERF_SurfaceLayer.H:743
amrex::Vector< amrex::MultiFab > z_0
Definition: ERF_SurfaceLayer.H:709
amrex::Vector< std::unique_ptr< amrex::MultiFab > > w_star
Definition: ERF_SurfaceLayer.H:727
amrex::Vector< std::unique_ptr< amrex::MultiFab > > u_star
Definition: ERF_SurfaceLayer.H:726
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_surf
Definition: ERF_SurfaceLayer.H:733
amrex::Vector< amrex::MultiFab * > m_Hwave_lev
Definition: ERF_SurfaceLayer.H:742
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_star
Definition: ERF_SurfaceLayer.H:728
amrex::Vector< amrex::MultiFab * > m_eddyDiffs_lev
Definition: ERF_SurfaceLayer.H:744
amrex::Vector< std::unique_ptr< amrex::MultiFab > > olen
Definition: ERF_SurfaceLayer.H:730
amrex::Vector< std::unique_ptr< amrex::MultiFab > > pblh
Definition: ERF_SurfaceLayer.H:731
@ ng
Definition: ERF_Morrison.H:48
Here is the call graph for this function:

◆ compute_pblh() [1/2]

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

◆ compute_pblh() [2/2]

template<typename PBLHeightEstimator >
void SurfaceLayer::compute_pblh ( const int &  lev,
Vector< Vector< MultiFab >> &  vars,
MultiFab *  z_phys_cc,
const PBLHeightEstimator &  est,
const MoistureComponentIndices moisture_indices 
)
913 {
914  est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
915  vars[lev][Vars::cons],m_lmask_lev[lev][0],
916  moisture_indices);
917 }
@ cons
Definition: ERF_IndexDefines.H:158
Here is the call graph for this function:

◆ compute_SurfaceLayer_bcs() [1/2]

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

◆ compute_SurfaceLayer_bcs() [2/2]

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

Function to calculate MOST fluxes for populating ghost cells.

Parameters
[in]levCurrent level
[in,out]mfsMultiFabs to populate
[in]eddyDiffsDiffusion coefficients from turbulence model
[in]flux_compstructure to compute fluxes
413 {
414  bool rotate = m_rotate;
415  const int klo = m_geom[lev].Domain().smallEnd(2);
416  const auto& dxInv = m_geom[lev].InvCellSizeArray();
417  for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
418  {
419  // Get field arrays
420  const auto cons_arr = mfs[Vars::cons]->array(mfi);
421  const auto velx_arr = mfs[Vars::xvel]->array(mfi);
422  const auto vely_arr = mfs[Vars::yvel]->array(mfi);
423  const auto velz_arr = mfs[Vars::zvel]->array(mfi);
424 
425  // Diffusive stress vars
426  auto t13_arr = Tau_lev[TauType::tau13]->array(mfi);
427  auto t31_arr = (Tau_lev[TauType::tau31]) ? Tau_lev[TauType::tau31]->array(mfi) : Array4<Real>{};
428 
429  auto t23_arr = Tau_lev[TauType::tau23]->array(mfi);
430  auto t32_arr = (Tau_lev[TauType::tau32]) ? Tau_lev[TauType::tau32]->array(mfi) : Array4<Real>{};
431 
432  auto hfx3_arr = zheat_flux->array(mfi);
433  auto qfx3_arr = (zqv_flux) ? zqv_flux->array(mfi) : Array4<Real>{};
434 
435  // Rotated stress vars
436  auto t11_arr = (m_rotate) ? Tau_lev[TauType::tau11]->array(mfi) : Array4<Real>{};
437  auto t22_arr = (m_rotate) ? Tau_lev[TauType::tau22]->array(mfi) : Array4<Real>{};
438  auto t33_arr = (m_rotate) ? Tau_lev[TauType::tau33]->array(mfi) : Array4<Real>{};
439  auto t12_arr = (m_rotate) ? Tau_lev[TauType::tau12]->array(mfi) : Array4<Real>{};
440  auto t21_arr = (m_rotate) ? Tau_lev[TauType::tau21]->array(mfi) : Array4<Real>{};
441 
442  auto hfx1_arr = (m_rotate) ? xheat_flux->array(mfi) : Array4<Real>{};
443  auto hfx2_arr = (m_rotate) ? yheat_flux->array(mfi) : Array4<Real>{};
444  auto qfx1_arr = (m_rotate && xqv_flux) ? xqv_flux->array(mfi) : Array4<Real>{};
445  auto qfx2_arr = (m_rotate && yqv_flux) ? yqv_flux->array(mfi) : Array4<Real>{};
446 
447  // Terrain
448  const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};
449 
450  // Get average arrays
451  const auto *const u_mean = m_ma.get_average(lev,0);
452  const auto *const v_mean = m_ma.get_average(lev,1);
453  const auto *const t_mean = m_ma.get_average(lev,2);
454  const auto *const q_mean = m_ma.get_average(lev,3);
455  const auto *const u_mag_mean = m_ma.get_average(lev,5);
456 
457  const auto um_arr = u_mean->array(mfi);
458  const auto vm_arr = v_mean->array(mfi);
459  const auto tm_arr = t_mean->array(mfi);
460  const auto qm_arr = q_mean->array(mfi);
461  const auto umm_arr = u_mag_mean->array(mfi);
462 
463  // Get derived arrays
464  const auto u_star_arr = u_star[lev]->array(mfi);
465  const auto t_star_arr = t_star[lev]->array(mfi);
466  const auto q_star_arr = q_star[lev]->array(mfi);
467  const auto t_surf_arr = t_surf[lev]->array(mfi);
468  const auto q_surf_arr = q_surf[lev]->array(mfi);
469 
470  // Get LSM fluxes
471  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
472  Array4<int> {};
473  auto lsm_t_flux_arr = Array4<Real> {};
474  auto lsm_q_flux_arr = Array4<Real> {};
475  auto lsm_tau13_arr = Array4<Real> {};
476  auto lsm_tau23_arr = Array4<Real> {};
477  for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
478  if (toLower(m_lsm_flux_name[n]) == "t_flux") { lsm_t_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
479  if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
480  if (toLower(m_lsm_flux_name[n]) == "tau13") { lsm_tau13_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
481  if (toLower(m_lsm_flux_name[n]) == "tau23") { lsm_tau23_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
482  }
483 
484 
485  // Rho*Theta flux
486  //============================================================================
487  Box bx = mfi.tilebox();
488  if (bx.smallEnd(2) != klo) { continue; }
489  bx.makeSlab(2,klo);
490  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
491  {
492  // Valid theta flux from LSM and over land
493  Real Tflux;
494  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
495  if (lsm_t_flux_arr && is_land) {
496  Tflux = lsm_t_flux_arr(i,j,0);
497  } else {
498  Tflux = flux_comp.compute_t_flux(i, j, k,
499  cons_arr, velx_arr, vely_arr,
500  umm_arr, tm_arr, u_star_arr,
501  t_star_arr, t_surf_arr);
502 
503  }
504 
505  // Do scalar flux rotations?
506  if (rotate) {
507  rotate_scalar_flux(i, j, k, Tflux, dxInv, zphys_arr,
508  hfx1_arr, hfx2_arr, hfx3_arr);
509  } else {
510  hfx3_arr(i,j,klo) = Tflux;
511  }
512  });
513 
514  // Rho*Qv flux
515  //============================================================================
516  if (use_moisture) {
517  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
518  {
519  // Valid qv flux from LSM and over land
520  Real Qflux;
521  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
522  if (lsm_q_flux_arr && is_land) {
523  Qflux = lsm_q_flux_arr(i,j,0);
524  } else {
525  Qflux = flux_comp.compute_q_flux(i, j, k,
526  cons_arr, velx_arr, vely_arr,
527  umm_arr, qm_arr, u_star_arr,
528  q_star_arr, q_surf_arr);
529  }
530 
531  // Do scalar flux rotations?
532  if (rotate) {
533  rotate_scalar_flux(i, j, k, Qflux, dxInv, zphys_arr,
534  qfx1_arr, qfx2_arr, qfx3_arr);
535  } else {
536  qfx3_arr(i,j,k) = Qflux;
537  }
538  });
539  } // custom
540 
541  if (!rotate) {
542  // Rho*u flux
543  //============================================================================
544  Box bxx = surroundingNodes(bx,0);
545  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
546  {
547  // Valid tau13 from LSM and over land
548  Real stressx;
549  int is_land_hi = (lmask_arr) ? lmask_arr(i ,j,0) : 1;
550  int is_land_lo = (lmask_arr) ? lmask_arr(i-1,j,0) : 1;
551  if (lsm_tau13_arr && (is_land_hi || is_land_lo)) {
552  stressx = zero;
553  if (!is_land_hi || !is_land_lo) {
554  stressx += myhalf * flux_comp.compute_u_flux(i, j, k,
555  cons_arr, velx_arr, vely_arr,
556  umm_arr, um_arr, u_star_arr);
557  }
558  if (is_land_hi) {
559  stressx += myhalf * lsm_tau13_arr(i ,j,0);
560  }
561  if (is_land_lo) {
562  stressx += myhalf * lsm_tau13_arr(i-1,j,0);
563  }
564  } else {
565  stressx = flux_comp.compute_u_flux(i, j, k,
566  cons_arr, velx_arr, vely_arr,
567  umm_arr, um_arr, u_star_arr);
568  }
569 
570  t13_arr(i,j,k) = stressx;
571  if (t31_arr) { t31_arr(i,j,k) = stressx; }
572  });
573 
574  // Rho*v flux
575  //============================================================================
576  Box bxy = surroundingNodes(bx,1);
577  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
578  {
579  // Valid tau13 from LSM and over land
580  Real stressy;
581  int is_land_hi = (lmask_arr) ? lmask_arr(i,j ,0) : 1;
582  int is_land_lo = (lmask_arr) ? lmask_arr(i,j-1,0) : 1;
583  if (lsm_tau23_arr && (is_land_hi || is_land_lo)) {
584  stressy = zero;
585  if (!is_land_hi || !is_land_lo) {
586  stressy += myhalf * flux_comp.compute_v_flux(i, j, k,
587  cons_arr, velx_arr, vely_arr,
588  umm_arr, vm_arr, u_star_arr);
589  }
590  if (is_land_hi) {
591  stressy += myhalf * lsm_tau23_arr(i,j ,0);
592  }
593  if (is_land_lo) {
594  stressy += myhalf * lsm_tau23_arr(i,j-1,0);
595  }
596  } else {
597  stressy = flux_comp.compute_v_flux(i, j, k,
598  cons_arr, velx_arr, vely_arr,
599  umm_arr, vm_arr, u_star_arr);
600  }
601 
602  t23_arr(i,j,k) = stressy;
603  if (t32_arr) { t32_arr(i,j,k) = stressy; }
604  });
605  } else {
606  // All fluxes with rotation
607  //============================================================================
608  Box bxxy = convert(bx, IntVect(1,1,0));
609  ParallelFor(bxxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
610  {
611  Real stresst = flux_comp.compute_u_flux(i, j, k,
612  cons_arr, velx_arr, vely_arr,
613  umm_arr, um_arr, u_star_arr);
614  rotate_stress_tensor(i, j, k, stresst, dxInv, zphys_arr,
615  velx_arr, vely_arr, velz_arr,
616  t11_arr, t22_arr, t33_arr,
617  t12_arr, t21_arr,
618  t13_arr, t31_arr,
619  t23_arr, t32_arr);
620  });
621  }
622  } // mfiter
623 }
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
@ tau12
Definition: ERF_DataStruct.H:32
@ tau23
Definition: ERF_DataStruct.H:32
@ tau33
Definition: ERF_DataStruct.H:32
@ tau22
Definition: ERF_DataStruct.H:32
@ tau11
Definition: ERF_DataStruct.H:32
@ tau32
Definition: ERF_DataStruct.H:32
@ tau31
Definition: ERF_DataStruct.H:32
@ tau21
Definition: ERF_DataStruct.H:32
@ tau13
Definition: ERF_DataStruct.H:32
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
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
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
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_flux_lev
Definition: ERF_SurfaceLayer.H:739
amrex::Vector< std::string > m_lsm_flux_name
Definition: ERF_SurfaceLayer.H:741
@ xvel
Definition: ERF_IndexDefines.H:159
@ zvel
Definition: ERF_IndexDefines.H:161
@ yvel
Definition: ERF_IndexDefines.H:160
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< 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 eb_ ebfact,
const FluxCalc &  flux_comp 
)

◆ compute_SurfaceLayer_bcs_EB() [2/2]

template<typename FluxCalc >
void SurfaceLayer::compute_SurfaceLayer_bcs_EB ( const int &  lev,
Vector< const MultiFab * >  mfs,
Vector< 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,
[[maybe_unused] ] const eb_ ebfact,
const FluxCalc &  flux_comp 
)

Function to calculate MOST fluxes for populating ghost cells.

Parameters
[in]levCurrent level
[in,out]mfsMultiFabs to populate
[in]eddyDiffsDiffusion coefficients from turbulence model
[in]flux_compstructure to compute fluxes
646 {
647  const int klo = m_geom[lev].Domain().smallEnd(2);
648  // const auto& dxInv = m_geom[lev].InvCellSizeArray();
649  for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
650  {
651  // Get field arrays
652  const auto cons_arr = mfs[Vars::cons]->array(mfi);
653  const auto velx_arr = mfs[Vars::xvel]->array(mfi);
654  const auto vely_arr = mfs[Vars::yvel]->array(mfi);
655 
656  // Diffusive stress vars
657  auto t13_arr = Tau_EB[EBTauType::tau_eb13]->array(mfi);
658  auto t23_arr = Tau_EB[EBTauType::tau_eb23]->array(mfi);
659 
660  auto hfx3_arr = Hfx3_EB->array(mfi);
661 
662  // Get average arrays
663  const auto *const u_mean = m_ma.get_average(lev,0);
664  const auto *const v_mean = m_ma.get_average(lev,1);
665  const auto *const t_mean = m_ma.get_average(lev,2);
666  // const auto *const q_mean = m_ma.get_average(lev,3);
667  const auto *const u_mag_mean = m_ma.get_average(lev,5);
668  const auto *const k_indx = m_ma.get_k_indices(lev);
669 
670  const auto um_arr = u_mean->array(mfi);
671  const auto vm_arr = v_mean->array(mfi);
672  const auto tm_arr = t_mean->array(mfi);
673  // const auto qm_arr = q_mean->array(mfi);
674  const auto umm_arr = u_mag_mean->array(mfi);
675  const auto k_arr = k_indx->const_array(mfi);
676 
677  // Get derived arrays
678  const auto u_star_arr = u_star[lev]->array(mfi);
679  const auto t_star_arr = t_star[lev]->array(mfi);
680  const auto t_surf_arr = t_surf[lev]->array(mfi);
681 
682  // Rho*Theta flux
683  //============================================================================
684  Box bx = mfi.tilebox();
685 
686  if (bx.smallEnd(2) != klo) { continue; }
687  bx.makeSlab(2,klo);
688  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
689  {
690  int mk = k_arr(i,j,0);
691 
692  Real Tflux = flux_comp.compute_t_flux(i, j, mk,
693  cons_arr, velx_arr, vely_arr,
694  umm_arr, tm_arr, u_star_arr,
695  t_star_arr, t_surf_arr);
696  hfx3_arr(i,j,mk) = Tflux;
697  });
698 
699  // Rho*u flux
700  //============================================================================
701  Box bxx = surroundingNodes(bx,0);
702  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
703  {
704  int mk = k_arr(i,j,0);
705 
706  Real stressx = flux_comp.compute_u_flux(i, j, mk,
707  cons_arr, velx_arr, vely_arr,
708  umm_arr, um_arr, u_star_arr);
709  t13_arr(i,j,mk) = stressx;
710  });
711 
712  // Rho*v flux
713  //============================================================================
714  Box bxy = surroundingNodes(bx,1);
715  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
716  {
717  int mk = k_arr(i,j,0);
718 
719  Real stressy = flux_comp.compute_v_flux(i, j, mk,
720  cons_arr, velx_arr, vely_arr,
721  umm_arr, vm_arr, u_star_arr);
722  t23_arr(i,j,mk) = stressy;
723  });
724  } // mfiter
725 }
@ tau_eb23
Definition: ERF_EBStruct.H:16
@ tau_eb13
Definition: ERF_EBStruct.H:16
const amrex::iMultiFab * get_k_indices(const int &lev) const
Definition: ERF_MOSTAverage.H:111
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 
)
818 {
819  // NOTE: We have already tested a moisture model exists
820 
821  // Populate q_surf with qsat over water
822  auto dz = m_geom[lev].CellSize(2);
823  const int klo = m_geom[lev].Domain().smallEnd(2);
824  for (MFIter mfi(*q_surf[lev]); mfi.isValid(); ++mfi)
825  {
826  Box gtbx = mfi.growntilebox();
827 
828  if (gtbx.smallEnd(2) != klo) { continue; }
829 
830  auto t_surf_arr = t_surf[lev]->array(mfi);
831  auto q_surf_arr = q_surf[lev]->array(mfi);
832  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
833  Array4<int> {};
834  const auto cons_arr = cons_in.const_array(mfi);
835  const auto z_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
836  Array4<const Real> {};
837 
838  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
839  {
840  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
841  if (!is_land) {
842  auto deltaZ = (z_arr) ? Compute_Zrel_AtCellCenter(i,j,k,z_arr) :
843  myhalf*dz;
844  auto Rho = cons_arr(i,j,k,Rho_comp);
845  auto RTh = cons_arr(i,j,k,RhoTheta_comp);
846  auto Qv = cons_arr(i,j,k,RhoQ1_comp) / Rho;
847  auto P_cc = getPgivenRTh(RTh, Qv);
848  P_cc += Rho*CONST_GRAV*deltaZ;
849  P_cc *= Real(0.01);
850  erf_qsatw(t_surf_arr(i,j,k), P_cc, q_surf_arr(i,j,k));
851  }
852  });
853  }
854  q_surf[lev]->FillBoundary(m_geom[lev].periodicity());
855 }
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:31
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
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
const Real dz
Definition: ERF_InitCustomPert_Bubble.H:25
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:171
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
Here is the call graph for this function:

◆ fill_tsurf_with_sst_and_tsk()

void SurfaceLayer::fill_tsurf_with_sst_and_tsk ( const int &  lev,
const amrex::Real time 
)
730 {
731  int n_times_in_sst = m_sst_lev[lev].size();
732 
734 
735  int n_time_lo, n_time_hi;
736  Real alpha;
737 
738  if (n_times_in_sst > 1) {
739  n_time_lo = static_cast<int>( elapsed_time_since_start_low / dT);
740  alpha = (elapsed_time_since_start_low - n_time_lo * dT) / dT;
741 
743 
744  n_time_hi = n_time_lo + 1;
745 
746  // Do not over run the last sst file
747  if (m_start_low_time + elapsed_time_since_start_low >= m_final_low_time) {
748  n_time_lo = m_sst_lev[lev].size()-1;
749  n_time_hi = n_time_lo;
750  alpha = zero;
751  }
752 
753  AMREX_ALWAYS_ASSERT( (n_time_lo >= 0) && (n_time_hi < m_sst_lev[lev].size()));
754  } else {
755  n_time_lo = 0;
756  n_time_hi = 0;
757  alpha = one;
758  }
760 
761  Real oma = one - alpha;
762 
763  // Define a default land surface temperature if we don't read in tsk
765 
766  bool use_tsk = (m_tsk_lev[lev][0]);
767  bool ignore_sst = m_ignore_sst;
768 
769  const int klo = m_geom[lev].Domain().smallEnd(2);
770 
771  // Populate t_surf
772  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
773  {
774  Box gtbx = mfi.growntilebox();
775 
776  if (gtbx.smallEnd(2) != klo) { continue; }
777 
778  auto t_surf_arr = t_surf[lev]->array(mfi);
779 
780  const auto sst_lo_arr = m_sst_lev[lev][n_time_lo]->const_array(mfi);
781  const auto sst_hi_arr = m_sst_lev[lev][n_time_hi]->const_array(mfi);
782 
783  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
784  Array4<int> {};
785 
786  if (use_tsk) {
787  const auto tsk_arr = m_tsk_lev[lev][n_time_lo]->const_array(mfi);
788  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
789  {
790  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
791  if (!is_land && !ignore_sst) {
792  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
793  + alpha * sst_hi_arr(i,j,k);
794  } else {
795  t_surf_arr(i,j,k) = tsk_arr(i,j,k);
796  }
797  });
798  } else {
799  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
800  {
801  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
802  if (!is_land) {
803  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
804  + alpha * sst_hi_arr(i,j,k);
805  } else {
806  t_surf_arr(i,j,k) = lst;
807  }
808  });
809  }
810  }
811  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
812 }
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_sst_lev
Definition: ERF_SurfaceLayer.H:735
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_tsk_lev
Definition: ERF_SurfaceLayer.H:736
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
616 { return m_lmask_lev[lev][0]; }

◆ get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf ( const int &  lev)
859 {
860  const int klo = m_geom[lev].Domain().smallEnd(2);
861  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
862  {
863  Box gtbx = mfi.growntilebox();
864 
865  if (gtbx.smallEnd(2) != klo) { continue; }
866 
867  // NOTE: LSM does not carry lateral ghost cells.
868  // This copies the valid box into the ghost cells.
869  // Fillboundary is called after this to pick up the
870  // interior ghost and periodic directions.
871  Box vbx = mfi.validbox();
872  int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
873  int j_lo = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1);
874 
875  auto t_surf_arr = t_surf[lev]->array(mfi);
876  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
877  Array4<int> {};
878  const auto lsm_arr = m_lsm_data_lev[lev][m_lsm_tsurf_indx]->const_array(mfi);
879 
880  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
881  {
882  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
883  if (is_land) {
884  int li = amrex::min(amrex::max(i, i_lo), i_hi);
885  int lj = amrex::min(amrex::max(j, j_lo), j_hi);
886  t_surf_arr(i,j,k) = lsm_arr(li,lj,k);
887  }
888  });
889  }
890 }
int m_lsm_tsurf_indx
Definition: ERF_SurfaceLayer.H:717
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_data_lev
Definition: ERF_SurfaceLayer.H:738
Here is the call graph for this function:

◆ get_mac_avg()

const amrex::MultiFab* SurfaceLayer::get_mac_avg ( const int &  lev,
int  comp 
)
inline
600  {
601  return m_ma.get_average(lev, comp);
602  }
Here is the call graph for this function:

◆ get_olen()

amrex::MultiFab* SurfaceLayer::get_olen ( const int &  lev)
inline
595 { return olen[lev].get(); }

◆ get_pblh()

amrex::MultiFab* SurfaceLayer::get_pblh ( const int &  lev)
inline
597 { return pblh[lev].get(); }

◆ get_q_star()

amrex::MultiFab* SurfaceLayer::get_q_star ( const int &  lev)
inline
593 { return q_star[lev].get(); }

◆ get_q_surf()

amrex::MultiFab* SurfaceLayer::get_q_surf ( const int &  lev)
inline
607 { return q_surf[lev].get(); }

◆ get_t_star()

amrex::MultiFab* SurfaceLayer::get_t_star ( const int &  lev)
inline
591 { return t_star[lev].get(); }

◆ get_t_surf()

amrex::MultiFab* SurfaceLayer::get_t_surf ( const int &  lev)
inline
604 { return t_surf[lev].get(); }

◆ get_u_star()

amrex::MultiFab* SurfaceLayer::get_u_star ( const int &  lev)
inline
587 { return u_star[lev].get(); }

◆ get_w_star()

amrex::MultiFab* SurfaceLayer::get_w_star ( const int &  lev)
inline
589 { return w_star[lev].get(); }

◆ get_z0()

amrex::MultiFab* SurfaceLayer::get_z0 ( const int &  lev)
inline
612 { return &z_0[lev]; }

◆ get_zref()

amrex::Real SurfaceLayer::get_zref ( const int &  lev)
inline
610 { 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
614 { return m_var_z0; }
bool m_var_z0
Definition: ERF_SurfaceLayer.H:710

◆ impose_SurfaceLayer_bcs()

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

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

Parameters
[in]levCurrent level
[in,out]mfsMultiFabs to populate
[in]eddyDiffsDiffusion coefficients from turbulence model
320 {
322  moeng_flux flux_comp;
323  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
324  xheat_flux, yheat_flux, zheat_flux,
325  xqv_flux, yqv_flux, zqv_flux,
326  z_phys, flux_comp);
327  } else if (flux_type == FluxCalcType::DONELAN) {
328  donelan_flux flux_comp;
329  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
330  xheat_flux, yheat_flux, zheat_flux,
331  xqv_flux, yqv_flux, zqv_flux,
332  z_phys, flux_comp);
333  } else if (flux_type == FluxCalcType::ROTATE) {
334  rotate_flux flux_comp;
335  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
336  xheat_flux, yheat_flux, zheat_flux,
337  xqv_flux, yqv_flux, zqv_flux,
338  z_phys, flux_comp);
339  } else if (flux_type == FluxCalcType::RICO) {
341  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
342  xheat_flux, yheat_flux, zheat_flux,
343  xqv_flux, yqv_flux, zqv_flux,
344  z_phys, flux_comp);
345  } else if (flux_type == FluxCalcType::BULK_COEFF) {
346  bulk_coeff_flux flux_comp(m_Cd, m_Ch, m_Cq);
347  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
348  xheat_flux, yheat_flux, zheat_flux,
349  xqv_flux, yqv_flux, zqv_flux,
350  z_phys, flux_comp);
351  } else if (flux_type == FluxCalcType::CUSTOM) {
352  custom_flux flux_comp(specified_rho_surf);
353  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
354  xheat_flux, yheat_flux, zheat_flux,
355  xqv_flux, yqv_flux, zqv_flux,
356  z_phys, flux_comp);
357  } else {
358  amrex::Abort("Unknown surface layer flux calculation type");
359  }
360 }
bool specified_rho_surf
Definition: ERF_SurfaceLayer.H:705
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)
Definition: ERF_MOSTStress.H:2181
Definition: ERF_MOSTStress.H:2064
Definition: ERF_MOSTStress.H:1931
Definition: ERF_MOSTStress.H:1768
Definition: ERF_MOSTStress.H:2308
Definition: ERF_MOSTStress.H:2436

◆ impose_SurfaceLayer_bcs_EB()

void SurfaceLayer::impose_SurfaceLayer_bcs_EB ( 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 eb_ ebfact 
)

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

Parameters
[in]levCurrent level
[in,out]mfsMultiFabs to populate
[in]eddyDiffsDiffusion coefficients from turbulence model
380 {
382  moeng_flux flux_comp;
383  compute_SurfaceLayer_bcs_EB(lev, mfs, Tau_EB,
384  xheat_flux, yheat_flux, Hfx3_EB,
385  xqv_flux, yqv_flux, zqv_flux,
386  ebfact, flux_comp);
387  } else {
388  amrex::Abort("Not implemented surface layer flux calculation type for EB");
389  }
390 }
void compute_SurfaceLayer_bcs_EB(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 eb_ &ebfact, const FluxCalc &flux_comp)

◆ 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) 
)
925 {
926  Print() << "Initializing TKE from surface layer ustar on level " << lev << std::endl;
927 
928  // Handle vertical decomposition by selectively copying into
929  // a FArrayBox section on each rank. Then doing a reduce real sum
930  // and broadcasting to each rank. No mask since all CC data
931  const int klo = m_geom[lev].Domain().smallEnd(2);
932  Box bx_lo = u_star[lev]->boxArray().minimalBox();
933  FArrayBox u_star_lo(bx_lo, 1); u_star_lo.setVal<RunOn::Device>(0);
934  FArrayBox z_surf_lo(bx_lo, 1); z_surf_lo.setVal<RunOn::Device>(0);
935  Real* ustar_ptr = u_star_lo.dataPtr();
936  Real* zsurf_ptr = z_surf_lo.dataPtr();
937  for (MFIter mfi(cons); mfi.isValid(); ++mfi)
938  {
939  Box vbx = mfi.validbox();
940  if (vbx.smallEnd(2) != klo) { continue; }
941  vbx.makeSlab(2,0);
942 
943  auto const& u_star_arr = u_star[lev]->const_array(mfi);
944  auto u_star_all = u_star_lo.array();
945 
946  auto const& z_phys_arr = z_phys_nd->const_array(mfi);
947  auto z_surf_all = z_surf_lo.array();
948 
949  ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
950  {
951  u_star_all(i,j,0) = u_star_arr(i,j,0);
952  z_surf_all(i,j,0) = fourth * ( z_phys_arr(i ,j ,klo) + z_phys_arr(i+1,j ,klo)
953  + z_phys_arr(i ,j+1,klo) + z_phys_arr(i+1,j+1,klo) );
954  });
955  }
956  ParallelDescriptor::ReduceRealSum(ustar_ptr, bx_lo.numPts());
957  ParallelDescriptor::ReduceRealSum(zsurf_ptr, bx_lo.numPts());
958 
959  // Now work on all boxes (ustar has been filled above)
960  constexpr Real small = Real(0.01);
961  for (MFIter mfi(cons); mfi.isValid(); ++mfi)
962  {
963  Box vbx = mfi.validbox();
964 
965  auto const& u_star_arr = u_star_lo.const_array();
966  auto const& z_surf_arr = z_surf_lo.const_array();
967  auto const& z_phys_arr = z_phys_nd->const_array(mfi);
968 
969  auto const& cons_arr = cons.array(mfi);
970 
971  ParallelFor(vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
972  {
973  Real rho = cons_arr(i, j, k, Rho_comp);
974  Real ust = u_star_arr(i, j, 0);
975  Real tke0 = tkefac * ust * ust; // surface value
976  Real zagl = Compute_Z_AtCellCenter(i, j, k, z_phys_arr) - z_surf_arr(i,j,0);
977 
978  // linearly tapering profile -- following WRF, approximate top of
979  // PBL as ustar * zscale
980  cons_arr(i, j, k, RhoKE_comp) = rho * tke0 * std::max(
981  (ust * zscale - zagl) / (std::max(ust, small) * zscale),
982  small);
983  });
984  }
985 }
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
rho
Definition: ERF_InitCustomPert_Bubble.H:106
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
620  {
621  int lmask_min = amrex::ReduceMin(lmask, nghost, [=] AMREX_GPU_HOST_DEVICE(
622  amrex::Box const& bx, amrex::Array4<int const> const& lm_arr) -> int
623  {
624  int locmin = std::numeric_limits<int>::max();
625  const auto lo = lbound(bx);
626  const auto hi = ubound(bx);
627  for (int j = lo.y; j <= hi.y; ++j) {
628  for (int i = lo.x; i <= hi.x; ++i) {
629  locmin = std::min(locmin, lm_arr(i, j, 0));
630  }
631  }
632  return locmin;
633  });
634 
635  return lmask_min;
636  }

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
284  {
285  // Update MOST Average
287  Theta_prim, Qv_prim, Qr_prim,
288  z_phys_nd);
289 
290  // Get CC vars
291  amrex::MultiFab& mf = *(mfv[0]);
292 
293  amrex::ParmParse pp("erf");
294 
295  // Do we have a time-varying surface roughness that needs to be saved?
296  if (lev == 0) {
297  const int nghost = 0; // ghost cells not included
298  int lmask_min = lmask_min_reduce(*lmask_lev[0].get(), nghost);
299  amrex::ParallelDescriptor::ReduceIntMin(lmask_min);
300 
301  m_var_z0 = (lmask_min < 1) & (rough_type_sea != RoughCalcType::CONSTANT);
302  if (m_var_z0) {
303  std::string rough_sea_string{"charnock"};
304  pp.query("most.roughness_type_sea", rough_sea_string);
305  amrex::Print() << "Variable sea roughness (type " << rough_sea_string
306  << ")" << std::endl;
307  }
308  }
309 
310  if (m_eddyDiffs_lev.size() < lev+1) {
311  m_Hwave_lev.resize(nlevs);
312  m_Lwave_lev.resize(nlevs);
313  m_eddyDiffs_lev.resize(nlevs);
314 
315  m_lsm_data_lev.resize(nlevs);
316  m_lsm_flux_lev.resize(nlevs);
317 
318  m_sst_lev.resize(nlevs);
319  m_tsk_lev.resize(nlevs);
320  m_lmask_lev.resize(nlevs);
321 
322  // Size the MOST params for all levels
323  z_0.resize(nlevs);
324  u_star.resize(nlevs);
325  w_star.resize(nlevs);
326  t_star.resize(nlevs);
327  q_star.resize(nlevs);
328  t_surf.resize(nlevs);
329  q_surf.resize(nlevs);
330  olen.resize(nlevs);
331  pblh.resize(nlevs);
332  }
333 
334  // Get pointers to SST,TSK and LANDMASK data
335  int nt_tot_sst = sst_lev.size();
336  m_sst_lev[lev].resize(nt_tot_sst);
337  for (int nt(0); nt < nt_tot_sst; ++nt) {
338  m_sst_lev[lev][nt] = sst_lev[nt].get();
339  }
340  int nt_tot_tsk = static_cast<int>(tsk_lev.size());
341  m_tsk_lev[lev].resize(nt_tot_tsk);
342  for (int nt(0); nt < nt_tot_tsk; ++nt) {
343  m_tsk_lev[lev][nt] = tsk_lev[nt].get();
344  }
345  int nt_tot_lmask = static_cast<int>(lmask_lev.size());
346  m_lmask_lev[lev].resize(nt_tot_lmask);
347  for (int nt(0); nt < nt_tot_lmask; ++nt) {
348  m_lmask_lev[lev][nt] = lmask_lev[nt].get();
349  }
350 
351  // Get pointers to wave data
352  m_Hwave_lev[lev] = Hwave;
353  m_Lwave_lev[lev] = Lwave;
354  m_eddyDiffs_lev[lev] = eddyDiffs;
355 
356  // Get pointers to LSM data and Fluxes
357  int ndata = static_cast<int>(lsm_data.size());
358  int nflux = static_cast<int>(lsm_flux.size());
359  m_lsm_data_name.resize(ndata);
360  m_lsm_data_lev[lev].resize(ndata);
361  m_lsm_flux_name.resize(nflux);
362  m_lsm_flux_lev[lev].resize(nflux);
363  for (int n(0); n < ndata; ++n) {
364  m_lsm_data_name[n] = lsm_data_name[n];
365  m_lsm_data_lev[lev][n] = lsm_data[n];
366  if (amrex::toLower(lsm_data_name[n]) == "theta") {
367  m_has_lsm_tsurf = true;
368  m_lsm_tsurf_indx = n;
369  }
370  }
371  for (int n(0); n < nflux; ++n) {
372  m_lsm_flux_name[n] = lsm_flux_name[n];
373  m_lsm_flux_lev[lev][n] = lsm_flux[n];
374  }
375 
376  // Check if there is a user-specified roughness file to be read
377  std::string fname;
378  bool read_z0 = false;
379  if ( (flux_type == FluxCalcType::MOENG) ||
381  int count = pp.countval("most.roughness_file_name");
382  if (count > 1) {
383  AMREX_ALWAYS_ASSERT(count >= lev+1);
384  pp.query("most.roughness_file_name", fname, lev);
385  read_z0 = true;
386  } else if (count == 1) {
387  if (lev == 0) {
388  pp.query("most.roughness_file_name", fname);
389  } else {
390  // we will interpolate from the coarsest level
391  fname = "";
392  }
393  read_z0 = true;
394  }
395  // else use z0_const
396  }
397 
398  // Attributes for MFs and FABs
399  //--------------------------------------------------------
400  // Create a 2D ba, dm, & ghost cells
401  amrex::BoxArray ba = mf.boxArray();
402  amrex::BoxList bl2d = ba.boxList();
403  for (auto& b : bl2d) { b.setRange(2,0); }
404  amrex::BoxArray ba2d(std::move(bl2d));
405  const amrex::DistributionMapping& dm = mf.DistributionMap();
406  const int ncomp = 1;
407  amrex::IntVect ng{1,1,0};
408 
409  // Z0 heights FAB
410  //--------------------------------------------------------
411  z_0[lev].define(ba2d, dm, ncomp, ng);
412  z_0[lev].setVal(z0_const);
413  if (read_z0) {
414  read_custom_roughness(lev, fname);
415  }
416 
417  // 2D MFs for U*, T*, T_surf
418  //--------------------------------------------------------
419  u_star[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
420  u_star[lev]->setVal(1.E34);
421 
422  w_star[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
423  w_star[lev]->setVal(1.E34);
424 
425  t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
426  t_star[lev]->setVal(0.0); // default to neutral
427 
428  q_star[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
429  q_star[lev]->setVal(0.0); // default to dry
430 
431  olen[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
432  olen[lev]->setVal(1.E34);
433 
434  pblh[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
435  pblh[lev]->setVal(1.E34);
436 
437  t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
438  t_surf[lev]->setVal(default_land_surf_temp);
439 
440  q_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d, dm, ncomp, ng);
441  q_surf[lev]->setVal(default_land_surf_moist);
442 
443  // TODO: Do we want an enum struct for indexing?
444 
445  bool use_sst = (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]);
446  bool use_tsk = (!m_tsk_lev[lev].empty() && m_tsk_lev[lev][0]);
447  if (use_sst || use_tsk || m_has_lsm_tsurf) {
448  // Valid SST, TSK or LSM data; t_surf set before computing fluxes (avoids
449  // extended lambda capture) Note that land temp will be set from m_tsk_lev
450  // while sea temp will be set from m_sst_lev
452 
453  // Pathways in fill_tsurf_with_sst_and_tsk
454  amrex::Print() << "Using MOST with specified surface temperature ";
455  // NOTE: SST from the LOW file populates TSK in update_sst_tsk.
456  // So if we have TSK, it contains everything and has been
457  // sanity checked for valid SST values.
458  if (use_tsk) { m_ignore_sst = true; }
459  if (use_tsk) {
460  amrex::Print() << "(land: TSK, ";
461  } else {
462  amrex::Print() << "(land: T0, ";
463  }
464  if (use_tsk && !use_sst) {
465  amrex::Print() << "sea: TSK)" << std::endl;
466  } else {
467  amrex::Print() << "sea: SST)" << std::endl;
469  }
470  }
471  }
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:85
int lmask_min_reduce(amrex::iMultiFab &lmask, const int &nghost)
Definition: ERF_SurfaceLayer.H:618
amrex::Vector< std::string > m_lsm_data_name
Definition: ERF_SurfaceLayer.H:740
bool m_has_lsm_tsurf
Definition: ERF_SurfaceLayer.H:716
void read_custom_roughness(const int &lev, const std::string &fname)
Definition: ERF_SurfaceLayer.cpp:989
Here is the call graph for this function:

◆ read_custom_roughness()

void SurfaceLayer::read_custom_roughness ( const int &  lev,
const std::string &  fname 
)
991 {
992  // Read the file if we have it
993  if (!fname.empty()) {
994  // Only the ioproc reads the file
995  Gpu::HostVector<Real> m_x,m_y,m_z0;
996  if (ParallelDescriptor::IOProcessor()) {
997  Print()<<"Reading MOST roughness file at level " << lev << " : " << fname << std::endl;
998  std::ifstream file(fname);
999  Real value1,value2,value3;
1000  while(file>>value1>>value2>>value3){
1001  m_x.push_back(value1);
1002  m_y.push_back(value2);
1003  m_z0.push_back(value3);
1004  }
1005  file.close();
1006 
1007  AMREX_ALWAYS_ASSERT(m_x.size() == m_y.size());
1008  AMREX_ALWAYS_ASSERT(m_x.size() == m_z0.size());
1009  }
1010 
1011  // Broadcast the whole domain to every rank
1012  int ioproc = ParallelDescriptor::IOProcessorNumber();
1013  int nnode = m_x.size();
1014  ParallelDescriptor::Bcast(&nnode, 1, ioproc);
1015 
1016  if (!ParallelDescriptor::IOProcessor()) {
1017  m_x.resize(nnode);
1018  m_y.resize(nnode);
1019  m_z0.resize(nnode);
1020  }
1021  ParallelDescriptor::Bcast(m_x.data() , nnode, ioproc);
1022  ParallelDescriptor::Bcast(m_y.data() , nnode, ioproc);
1023  ParallelDescriptor::Bcast(m_z0.data(), nnode, ioproc);
1024 
1025  // Copy data to the GPU
1026  Gpu::DeviceVector<Real> d_x(nnode),d_y(nnode),d_z0(nnode);
1027  Gpu::copy(Gpu::hostToDevice, m_x.begin(), m_x.end(), d_x.begin());
1028  Gpu::copy(Gpu::hostToDevice, m_y.begin(), m_y.end(), d_y.begin());
1029  Gpu::copy(Gpu::hostToDevice, m_z0.begin(), m_z0.end(), d_z0.begin());
1030  Real* xp = d_x.data();
1031  Real* yp = d_y.data();
1032  Real* z0p = d_z0.data();
1033 
1034  // Each rank populates it's z_0[lev] MultiFab
1035  const int klo = m_geom[lev].Domain().smallEnd(2);
1036  for (MFIter mfi(z_0[lev]); mfi.isValid(); ++mfi)
1037  {
1038  Box gtbx = mfi.growntilebox();
1039 
1040  if (gtbx.smallEnd(2) != klo) { continue; }
1041 
1042  // Populate z_phys data
1043  Real tol = Real(1.0e-4);
1044  auto dx = m_geom[lev].CellSizeArray();
1045  auto ProbLoArr = m_geom[lev].ProbLoArray();
1046  int ilo = m_geom[lev].Domain().smallEnd(0);
1047  int jlo = m_geom[lev].Domain().smallEnd(1);
1048  int ihi = m_geom[lev].Domain().bigEnd(0);
1049  int jhi = m_geom[lev].Domain().bigEnd(1);
1050 
1051  Array4<Real> const& z0_arr = z_0[lev].array(mfi);
1052  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
1053  {
1054  // Clip indices for ghost-cells
1055  int ii = amrex::min(amrex::max(i,ilo),ihi);
1056  int jj = amrex::min(amrex::max(j,jlo),jhi);
1057 
1058  // Location of nodes
1059  Real x = ProbLoArr[0] + ii * dx[0];
1060  Real y = ProbLoArr[1] + jj * dx[1];
1061  int inode = ii + jj * (ihi-ilo+2); // stride is Nx+1
1062  if (std::sqrt(std::pow(x-xp[inode],2)+std::pow(y-yp[inode],2)) < tol) {
1063  z0_arr(i,j,klo) = z0p[inode];
1064  } else {
1065  // Unexpected list order, do brute force search
1066  Real z0loc = zero;
1067  bool found = false;
1068  for (int n=0; n<nnode; ++n) {
1069  Real delta=std::sqrt(std::pow(x-xp[n],2)+std::pow(y-yp[n],2));
1070  if (delta < tol) {
1071  found = true;
1072  z0loc = z0p[n];
1073  break;
1074  }
1075  }
1076  AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
1077  amrex::ignore_unused(found);
1078  z0_arr(i,j,klo) = z0loc;
1079  }
1080  });
1081  } // mfi
1082  } else {
1083  AMREX_ALWAYS_ASSERT(lev > 0);
1084 
1085  Print()<<"Interpolating MOST roughness at level " << lev << std::endl;
1086 
1087  // Create a BC mapper that uses FOEXTRAP at domain bndry
1088  Vector<int> bc_lo(3,ERFBCType::foextrap);
1089  Vector<int> bc_hi(3,ERFBCType::foextrap);
1090  Vector<BCRec> bcr; bcr.push_back(BCRec(bc_lo.data(),bc_hi.data()));
1091 
1092  // Create ref ratio
1093  IntVect ratio;
1094  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1095  ratio[idim] = m_geom[lev].Domain().length(idim) / m_geom[0].Domain().length(idim);
1096  }
1097 
1098  // Create interp object and interpolate from the coarsest grid
1099  MFInterpolater* interp = &mf_cell_cons_interp;
1100  interp->interp(z_0[0] , 0,
1101  z_0[lev], 0,
1102  1, z_0[lev].nGrowVect(),
1103  m_geom[0], m_geom[lev],
1104  m_geom[lev].Domain(),ratio,
1105  bcr, 0);
1106  }
1107 }
@ m_y
Definition: ERF_DataStruct.H:25
@ m_x
Definition: ERF_DataStruct.H:24
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
@ foextrap
Definition: ERF_IndexDefines.H:226

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
608 { q_surf[lev]->setVal(qsurf); }

◆ set_t_surf()

void SurfaceLayer::set_t_surf ( const int &  lev,
const amrex::Real  tsurf 
)
inline
605 { t_surf[lev]->setVal(tsurf); }

◆ update_fluxes()

void SurfaceLayer::update_fluxes ( const int &  lev,
const amrex::Real elapsed_time_since_start_low,
amrex::MultiFab &  cons_in,
const std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
const std::unique_ptr< amrex::MultiFab > &  walldist,
int  max_iters = 100 
)

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

Parameters
[in]levCurrent level
[in]max_itersmaximum iterations to use
18 {
19  // Update with SST/TSK data if we have a valid pointer
20  if (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]) {
21  fill_tsurf_with_sst_and_tsk(lev, elapsed_time_since_start_low);
22  }
23 
24  // Apply heating rate if needed
26  update_surf_temp(elapsed_time_since_start_low);
27  }
28 
29  // Update qsurf with qsat over sea
30  if (use_moisture) {
31  fill_qsurf_with_qsat(lev, cons_in, z_phys_nd);
32  }
33 
34  // Update land surface temp if we have a valid pointer
35  if (m_has_lsm_tsurf) { get_lsm_tsurf(lev); }
36 
37  // Fill interior ghost cells
38  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
39 
40  // Compute plane averages for all vars (regardless of flux type)
42 
43 
44  // ***************************************************************
45  // Iterate the fluxes if moeng type
46  // First iterate over land -- the only model for surface roughness
47  // over land is RoughCalcType::CONSTANT
48  // ***************************************************************
51  bool is_land = true;
52  // Do we have a constant flux for moisture over land?
53  bool cons_qflux = ( (moist_type == MoistCalcType::MOISTURE_FLUX) ||
57  surface_flux most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
58  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
59  } else {
60  amrex::Abort("Unknown value for rough_type_land");
61  }
64  surface_temp most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
65  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
66  } else {
67  amrex::Abort("Unknown value for rough_type_land");
68  }
69  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
73  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
74  } else {
75  amrex::Abort("Unknown value for rough_type_land");
76  }
77  } else {
78  amrex::Abort("Unknown value for theta_type");
79  }
80  } // MOENG -- LAND
81 
82  // ***************************************************************
83  // Iterate the fluxes if moeng type
84  // Next iterate over sea -- the models for surface roughness
85  // over sea are CHARNOCK, DONELAN, MODIFIED_CHARNOCK or WAVE_COUPLED
86  // ***************************************************************
89  bool is_land = false;
90  // NOTE: Do not allow default to adiabatic over sea (we have Qvs at surface)
91  // Do we have a constant flux for moisture over sea?
92  bool cons_qflux = (moist_type == MoistCalcType::MOISTURE_FLUX);
96  cnk_a, cnk_visc, cons_qflux);
97  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
100  depth, cons_qflux);
101  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
102  } else if (rough_type_sea == RoughCalcType::DONELAN) {
104  cons_qflux);
105  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
108  cons_qflux);
109  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
110  } else {
111  amrex::Abort("Unknown value for rough_type_sea");
112  }
113 
117  cnk_a, cnk_visc, cons_qflux);
118  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
121  depth, cons_qflux);
122  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
123  } else if (rough_type_sea == RoughCalcType::DONELAN) {
125  cons_qflux);
126  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
129  cons_qflux);
130  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
131  } else {
132  amrex::Abort("Unknown value for rough_type_sea");
133  }
134 
135  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
139  cnk_a, cnk_visc);
140  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
143  depth);
144  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
145  } else if (rough_type_sea == RoughCalcType::DONELAN) {
147  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
150  compute_fluxes(lev, max_iters, cons_in, most_flux, is_land);
151  } else {
152  amrex::Abort("Unknown value for rough_type_sea");
153  }
154  } else {
155  amrex::Abort("Unknown value for theta_type");
156  }
157 
158  } // MOENG -- SEA
159 
161  if (custom_rhosurf > 0) {
162  specified_rho_surf = true;
163  u_star[lev]->setVal(std::sqrt(custom_rhosurf) * custom_ustar);
164  t_star[lev]->setVal(custom_rhosurf * custom_tstar);
165  q_star[lev]->setVal(custom_rhosurf * custom_qstar);
166  } else {
167  u_star[lev]->setVal(custom_ustar);
168  t_star[lev]->setVal(custom_tstar);
169  q_star[lev]->setVal(custom_qstar);
170  }
171  }
172 
173  if (m_update_k_rans) {
174  const bool use_ref_theta = (theta_ref > 0);
176  const Real l_inv_Cmu2 = inv_Cmu2;
177  const int klo = m_geom[lev].Domain().smallEnd(2);
178  IntVect ng = u_star[lev]->nGrowVect(); ng[2] = 0;
179 
180  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi)
181  {
182  Box gpbx = mfi.tilebox(IntVect(0),ng);
183 
184  if (gpbx.smallEnd(2) != klo) { continue; }
185 
186  gpbx.makeSlab(2,klo);
187 
188  auto cons_arr = cons_in.array(mfi);
189  const auto& u_star_arr = u_star[lev]->const_array(mfi);
190  const auto& t_star_arr = t_star[lev]->const_array(mfi);
191  const auto& dist_arr = walldist->const_array(mfi);
192 
193  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
194  {
195  Real rho = cons_arr(i,j,k,Rho_comp);
196  if (t_star_arr(i,j,0) < -1e-8) {
197  // Only destabilizing buoyancy flux affects the boundary k
198  // tstar < 0 ==> B > 0
199  Real B = -CONST_GRAV * l_inv_theta0 * u_star_arr(i,j,0) * t_star_arr(i,j,0);
200  if (!use_ref_theta) {
201  B *= cons_arr(i,j,k,Rho_comp) /
202  cons_arr(i,j,k,RhoTheta_comp);
203  }
204 
205  // Axell & Liungman 2001, Eqn. 16
206  cons_arr(i,j,k,RhoKE_comp) = rho * l_inv_Cmu2 *
207  std::pow(
208  u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0)
209  + KAPPA * B * dist_arr(i,j,k),
210  two/three);
211  } else {
212  cons_arr(i,j,k,RhoKE_comp) = rho * l_inv_Cmu2 * u_star_arr(i,j,0) * u_star_arr(i,j,0);
213  }
214  });
215  }
216  }
217 }
constexpr amrex::Real three
Definition: ERF_Constants.H:9
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:30
constexpr amrex::Real two
Definition: ERF_Constants.H:8
const bool use_ref_theta
Definition: ERF_SetupDiff.H:16
const Real l_inv_theta0
Definition: ERF_SetupDiff.H:17
void compute_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:836
void update_surf_temp(const amrex::Real &time)
Definition: ERF_SurfaceLayer.H:566
void get_lsm_tsurf(const int &lev)
Definition: ERF_SurfaceLayer.cpp:858
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:815
void compute_fluxes(const int &lev, const int &max_iters, amrex::MultiFab &cons_in, const FluxIter &most_flux, bool is_land)
void fill_tsurf_with_sst_and_tsk(const int &lev, const amrex::Real &time)
Definition: ERF_SurfaceLayer.cpp:728
Definition: ERF_MOSTStress.H:189
Definition: ERF_MOSTStress.H:359
Definition: ERF_MOSTStress.H:278
Definition: ERF_MOSTStress.H:429
Definition: ERF_MOSTStress.H:140
Definition: ERF_MOSTStress.H:610
Definition: ERF_MOSTStress.H:840
Definition: ERF_MOSTStress.H:729
Definition: ERF_MOSTStress.H:948
Definition: ERF_MOSTStress.H:505
Definition: ERF_MOSTStress.H:1199
Definition: ERF_MOSTStress.H:1523
Definition: ERF_MOSTStress.H:1365
Definition: ERF_MOSTStress.H:1641
Definition: ERF_MOSTStress.H:1065
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
583  {
584  m_ma.update_field_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
585  }
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:269
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 
)
897 {
899  MYNNPBLH estimator;
900  compute_pblh(lev, vars, z_phys_cc, estimator, moisture_indices);
902  amrex::Error("YSU/MRF PBLH calc not implemented yet");
903  }
904 }
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)
Definition: ERF_PBLHeight.H:8

◆ update_sst_ptr()

void SurfaceLayer::update_sst_ptr ( const int  lev,
const int  itime,
amrex::MultiFab *  sst_ptr 
)
inline
638  {
639  m_sst_lev[lev][itime] = sst_ptr;
640  }

◆ update_surf_temp()

void SurfaceLayer::update_surf_temp ( const amrex::Real time)
inline
567  {
568  if (surf_heating_rate != 0) {
569  int nlevs = static_cast<int>(m_geom.size());
570  for (int lev = 0; lev < nlevs; lev++) {
571  t_surf[lev]->setVal(surf_temp + surf_heating_rate * time);
572  amrex::Print() << "Surface temp at t=" << time << ": "
573  << surf_temp + surf_heating_rate * time << std::endl;
574  }
575  }
576  }

◆ update_tsk_ptr()

void SurfaceLayer::update_tsk_ptr ( const int  lev,
const int  itime,
amrex::MultiFab *  tsk_ptr 
)
inline
642  {
643  m_tsk_lev[lev][itime] = tsk_ptr;
644  }

Member Data Documentation

◆ cnk_a

amrex::Real SurfaceLayer::cnk_a {amrex::Real(0.0185)}
private

Referenced by SurfaceLayer().

◆ cnk_visc

bool SurfaceLayer::cnk_visc {false}
private

Referenced by SurfaceLayer().

◆ custom_qstar

amrex::Real SurfaceLayer::custom_qstar {0}
private

Referenced by SurfaceLayer().

◆ custom_rhosurf

amrex::Real SurfaceLayer::custom_rhosurf {0}
private

Referenced by SurfaceLayer().

◆ custom_tstar

amrex::Real SurfaceLayer::custom_tstar {0}
private

Referenced by SurfaceLayer().

◆ custom_ustar

amrex::Real SurfaceLayer::custom_ustar {0}
private

Referenced by SurfaceLayer().

◆ default_land_surf_moist

amrex::Real SurfaceLayer::default_land_surf_moist {zero}
private

◆ default_land_surf_temp

amrex::Real SurfaceLayer::default_land_surf_temp {amrex::Real(300.)}
private

◆ depth

amrex::Real SurfaceLayer::depth {amrex::Real(30.0)}
private

Referenced by SurfaceLayer().

◆ flux_type

FluxCalcType SurfaceLayer::flux_type {FluxCalcType::MOENG}

◆ inv_Cmu2

amrex::Real SurfaceLayer::inv_Cmu2 = zero
private

Referenced by SurfaceLayer().

◆ m_Cd

amrex::Real SurfaceLayer::m_Cd = zero
private

Referenced by SurfaceLayer().

◆ m_Ch

amrex::Real SurfaceLayer::m_Ch = zero
private

Referenced by SurfaceLayer().

◆ m_Cq

amrex::Real SurfaceLayer::m_Cq = zero
private

Referenced by SurfaceLayer().

◆ m_eddyDiffs_lev

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

◆ m_final_low_time

amrex::Real SurfaceLayer::m_final_low_time
private

◆ m_geom

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

Referenced by update_surf_temp().

◆ m_has_lsm_tsurf

bool SurfaceLayer::m_has_lsm_tsurf = false
private

◆ m_Hwave_lev

amrex::Vector<amrex::MultiFab*> SurfaceLayer::m_Hwave_lev
private

◆ m_ignore_sst

bool SurfaceLayer::m_ignore_sst = false
private

◆ m_include_wstar

bool SurfaceLayer::m_include_wstar = false
private

Referenced by SurfaceLayer().

◆ m_lmask_lev

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

◆ m_low_time_interval

amrex::Real SurfaceLayer::m_low_time_interval
private

◆ m_lsm_data_lev

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

◆ m_lsm_data_name

amrex::Vector<std::string> SurfaceLayer::m_lsm_data_name
private

◆ m_lsm_flux_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceLayer::m_lsm_flux_lev
private

◆ m_lsm_flux_name

amrex::Vector<std::string> SurfaceLayer::m_lsm_flux_name
private

◆ m_lsm_tsurf_indx

int SurfaceLayer::m_lsm_tsurf_indx = -1
private

◆ m_Lwave_lev

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

◆ m_ma

MOSTAverage SurfaceLayer::m_ma
private

◆ m_rotate

bool SurfaceLayer::m_rotate = false
private

◆ m_sst_lev

amrex::Vector<amrex::Vector<amrex::MultiFab*> > SurfaceLayer::m_sst_lev
private

◆ m_start_low_time

amrex::Real SurfaceLayer::m_start_low_time
private

◆ m_tsk_lev

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

◆ m_update_k_rans

bool SurfaceLayer::m_update_k_rans = false
private

Referenced by SurfaceLayer().

◆ m_var_z0

bool SurfaceLayer::m_var_z0 {false}
private

◆ moist_type

MoistCalcType SurfaceLayer::moist_type {MoistCalcType::ADIABATIC}

Referenced by SurfaceLayer().

◆ olen

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

◆ pblh

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

◆ pblh_type

PBLHeightCalcType SurfaceLayer::pblh_type {PBLHeightCalcType::None}

Referenced by SurfaceLayer().

◆ q_star

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

◆ q_surf

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::q_surf
private

◆ rico_qsat_z0

amrex::Real SurfaceLayer::rico_qsat_z0 {amrex::Real(0.001)}
private

Referenced by SurfaceLayer().

◆ rico_theta_z0

amrex::Real SurfaceLayer::rico_theta_z0 {amrex::Real(298.0)}
private

Referenced by SurfaceLayer().

◆ rough_type_land

RoughCalcType SurfaceLayer::rough_type_land {RoughCalcType::CONSTANT}

Referenced by SurfaceLayer().

◆ rough_type_sea

RoughCalcType SurfaceLayer::rough_type_sea {RoughCalcType::CHARNOCK}

◆ specified_rho_surf

bool SurfaceLayer::specified_rho_surf {false}
private

◆ surf_heating_rate

amrex::Real SurfaceLayer::surf_heating_rate {0}
private

Referenced by SurfaceLayer(), and update_surf_temp().

◆ surf_moist

amrex::Real SurfaceLayer::surf_moist
private

Referenced by SurfaceLayer().

◆ surf_moist_flux

amrex::Real SurfaceLayer::surf_moist_flux {0}
private

Referenced by SurfaceLayer().

◆ surf_temp

amrex::Real SurfaceLayer::surf_temp
private

Referenced by SurfaceLayer(), and update_surf_temp().

◆ surf_temp_flux

amrex::Real SurfaceLayer::surf_temp_flux {0}
private

Referenced by SurfaceLayer().

◆ t_star

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::t_star
private

◆ t_surf

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::t_surf
private

◆ theta_ref

amrex::Real SurfaceLayer::theta_ref = zero
private

Referenced by SurfaceLayer().

◆ theta_type

◆ u_star

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

◆ use_moisture

bool SurfaceLayer::use_moisture
private

Referenced by SurfaceLayer().

◆ w_star

amrex::Vector<std::unique_ptr<amrex::MultiFab> > SurfaceLayer::w_star
private

◆ z0_const

amrex::Real SurfaceLayer::z0_const {amrex::Real(0.1)}
private

◆ z_0

amrex::Vector<amrex::MultiFab> SurfaceLayer::z_0
private

The documentation for this class was generated from the following files: