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

#include <ERF_ABLMost.H>

Collaboration diagram for ABLMost:

Public Types

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

Public Member Functions

 ABLMost (const amrex::Vector< amrex::Geometry > &geom, bool &use_exp_most, bool &use_rot_most, 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::Vector< std::unique_ptr< amrex::MultiFab >> &z_phys_nd, amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab >>> &sst_lev, amrex::Vector< amrex::Vector< std::unique_ptr< amrex::iMultiFab >>> &lmask_lev, amrex::Vector< amrex::Vector< amrex::MultiFab * >> lsm_data, amrex::Vector< amrex::Vector< amrex::MultiFab * >> lsm_flux, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Hwave, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Lwave, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &eddyDiffs, amrex::Real start_bdy_time=0.0, amrex::Real bdy_time_interval=0.0)
 
void update_fluxes (const int &lev, const amrex::Real &time, int max_iters=25)
 
template<typename FluxIter >
void compute_fluxes (const int &lev, const int &max_iters, const FluxIter &most_flux, bool is_land)
 
void impose_most_bcs (const int &lev, const amrex::Vector< amrex::MultiFab * > &mfs, amrex::MultiFab *xxmom_flux, amrex::MultiFab *yymom_flux, amrex::MultiFab *zzmom_flux, amrex::MultiFab *xymom_flux, amrex::MultiFab *yxmom_flux, amrex::MultiFab *xzmom_flux, amrex::MultiFab *zxmom_flux, amrex::MultiFab *yzmom_flux, amrex::MultiFab *zymom_flux, amrex::MultiFab *xheat_flux, amrex::MultiFab *yheat_flux, amrex::MultiFab *zheat_flux, amrex::MultiFab *xqv_flux, amrex::MultiFab *yqv_flux, amrex::MultiFab *zqv_flux, amrex::MultiFab *z_phys)
 
template<typename FluxCalc >
void compute_most_bcs (const int &lev, const amrex::Vector< amrex::MultiFab * > &mfs, amrex::MultiFab *xxmom_flux, amrex::MultiFab *yymom_flux, amrex::MultiFab *zzmom_flux, amrex::MultiFab *xymom_flux, amrex::MultiFab *yxmom_flux, amrex::MultiFab *xzmom_flux, amrex::MultiFab *zxmom_flux, amrex::MultiFab *yzmom_flux, amrex::MultiFab *zymom_flux, amrex::MultiFab *xheat_flux, amrex::MultiFab *yheat_flux, amrex::MultiFab *zheat_flux, amrex::MultiFab *xqv_flux, amrex::MultiFab *yqv_flux, amrex::MultiFab *zqv_flux, amrex::MultiFab *z_phys, const FluxCalc &flux_comp)
 
void time_interp_sst (const int &lev, const amrex::Real &time)
 
void get_lsm_tsurf (const int &lev)
 
void update_pblh (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars, amrex::MultiFab *z_phys_cc, const int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp)
 
template<typename PBLHeightEstimator >
void compute_pblh (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars, amrex::MultiFab *z_phys_cc, const PBLHeightEstimator &est, const int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp)
 
void read_custom_roughness (const int &lev, const std::string &fname)
 
void update_surf_temp (const amrex::Real &time)
 
void update_mac_ptrs (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_old, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Theta_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qv_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qr_prim)
 
const amrex::MultiFab * get_u_star (const int &lev)
 
const amrex::MultiFab * get_w_star (const int &lev)
 
const amrex::MultiFab * get_t_star (const int &lev)
 
const amrex::MultiFab * get_q_star (const int &lev)
 
const amrex::MultiFab * get_olen (const int &lev)
 
const amrex::MultiFab * get_pblh (const int &lev)
 
const amrex::MultiFab * get_mac_avg (const int &lev, int comp)
 
const amrex::MultiFab * get_t_surf (const int &lev)
 
amrex::Real get_zref ()
 
const amrex::FArrayBox * get_z0 (const int &lev)
 
bool have_variable_sea_roughness ()
 
const amrex::iMultiFab * get_lmask (const int &lev)
 
int lmask_min_reduce (amrex::iMultiFab &lmask, const int &nghost)
 
template<typename FluxCalc >
void compute_most_bcs (const int &lev, const Vector< MultiFab * > &mfs, MultiFab *xxmom_flux, MultiFab *yymom_flux, MultiFab *zzmom_flux, MultiFab *xymom_flux, MultiFab *yxmom_flux, MultiFab *xzmom_flux, MultiFab *zxmom_flux, MultiFab *yzmom_flux, MultiFab *zymom_flux, MultiFab *xheat_flux, MultiFab *yheat_flux, MultiFab *zheat_flux, MultiFab *xqv_flux, MultiFab *yqv_flux, MultiFab *zqv_flux, MultiFab *z_phys, const FluxCalc &flux_comp)
 
template<typename PBLHeightEstimator >
void compute_pblh (const int &lev, Vector< Vector< MultiFab >> &vars, MultiFab *z_phys_cc, const PBLHeightEstimator &est, int RhoQv_comp, int RhoQc_comp, int RhoQr_comp)
 

Public Attributes

FluxCalcType flux_type {FluxCalcType::MOENG}
 
ThetaCalcType theta_type {ThetaCalcType::ADIABATIC}
 
RoughCalcType rough_type_land {RoughCalcType::CONSTANT}
 
RoughCalcType rough_type_sea {RoughCalcType::CHARNOCK}
 
PBLHeightCalcType pblh_type {PBLHeightCalcType::None}
 

Private Attributes

bool use_moisture
 
bool m_exp_most = false
 
bool m_rotate = false
 
bool m_include_wstar = false
 
amrex::Real z0_const {0.1}
 
amrex::Real surf_temp
 
amrex::Real surf_heating_rate {0}
 
amrex::Real surf_temp_flux {0}
 
amrex::Real custom_ustar {0}
 
amrex::Real custom_tstar {0}
 
amrex::Real custom_qstar {0}
 
amrex::Real cnk_a {0.0185}
 
bool cnk_visc {false}
 
amrex::Real depth {30.0}
 
amrex::Real m_start_bdy_time
 
amrex::Real m_bdy_time_interval
 
amrex::Vector< amrex::Geometry > m_geom
 
amrex::Vector< amrex::FArrayBox > z_0
 
bool m_var_z0 {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< amrex::Vector< amrex::MultiFab * > > m_sst_lev
 
amrex::Vector< amrex::Vector< amrex::iMultiFab * > > m_lmask_lev
 
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_data_lev
 
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_flux_lev
 
amrex::Vector< amrex::MultiFab * > m_Hwave_lev
 
amrex::Vector< amrex::MultiFab * > m_Lwave_lev
 
amrex::Vector< amrex::MultiFab * > m_eddyDiffs_lev
 

Detailed Description

Monin-Obukhov surface layer profile

van der Laan, P., Kelly, M. C., & Sørensen, N. N. (2017). A new k-epsilon model consistent with Monin-Obukhov similarity theory. Wind Energy, 20(3), 479–489. https://doi.org/10.1002/we.2017

Consistent with Dyer (1974) formulation from page 57, Chapter 2, Modeling the vertical ABL structure in Modelling of Atmospheric Flow Fields, Demetri P Lalas and Corrado F Ratto, January 1996, https://doi.org/10.1142/2975.

Member Enumeration Documentation

◆ FluxCalcType

enum ABLMost::FluxCalcType
strong
Enumerator
MOENG 

Moeng functional form.

DONELAN 

Donelan functional form.

CUSTOM 

Custom constant flux functional form.

ROTATE 

Terrain rotation flux functional form.

474  {
475  MOENG = 0, ///< Moeng functional form
476  DONELAN, ///< Donelan functional form
477  CUSTOM, ///< Custom constant flux functional form
478  ROTATE ///< Terrain rotation flux functional form
479  };

◆ PBLHeightCalcType

Enumerator
None 
MYNN25 
YSU 
495  {
496  None, MYNN25, YSU
497  };

◆ RoughCalcType

Enumerator
CONSTANT 

Constant z0.

CHARNOCK 
MODIFIED_CHARNOCK 
DONELAN 
WAVE_COUPLED 
487  {
488  CONSTANT = 0, ///< Constant z0
489  CHARNOCK,
490  MODIFIED_CHARNOCK,
491  DONELAN,
492  WAVE_COUPLED
493  };

◆ ThetaCalcType

Enumerator
ADIABATIC 
HEAT_FLUX 

Heat-flux specified.

SURFACE_TEMPERATURE 

Surface temperature specified.

481  {
482  ADIABATIC = 0,
483  HEAT_FLUX, ///< Heat-flux specified
484  SURFACE_TEMPERATURE ///< Surface temperature specified
485  };

Constructor & Destructor Documentation

◆ ABLMost()

ABLMost::ABLMost ( const amrex::Vector< amrex::Geometry > &  geom,
bool &  use_exp_most,
bool &  use_rot_most,
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::Vector< std::unique_ptr< amrex::MultiFab >> &  z_phys_nd,
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab >>> &  sst_lev,
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::iMultiFab >>> &  lmask_lev,
amrex::Vector< amrex::Vector< amrex::MultiFab * >>  lsm_data,
amrex::Vector< amrex::Vector< amrex::MultiFab * >>  lsm_flux,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  Hwave,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  Lwave,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  eddyDiffs,
amrex::Real  start_bdy_time = 0.0,
amrex::Real  bdy_time_interval = 0.0 
)
inlineexplicit
52  : m_exp_most(use_exp_most),
53  m_rotate(use_rot_most),
54  m_start_bdy_time(start_bdy_time),
55  m_bdy_time_interval(bdy_time_interval),
56  m_geom(geom),
57  m_ma(geom,vars_old,Theta_prim,Qv_prim,Qr_prim,z_phys_nd)
58  {
59  // We have a moisture model if Qv_prim is a valid pointer
60  use_moisture = (Qv_prim[0].get());
61 
62  // Get roughness
63  amrex::ParmParse pp("erf");
64  pp.query("most.z0", z0_const);
65 
66  // Specify how to compute the flux
67  if (use_rot_most) {
69  } else {
70  std::string flux_string{"moeng"};
71  pp.query("most.flux_type", flux_string);
72  if (flux_string == "donelan") {
74  } else if (flux_string == "moeng") {
76  } else if (flux_string == "custom") {
78  } else {
79  amrex::Abort("Undefined MOST flux type!");
80  }
81  }
82 
83  // Include w* to handle free convection (Beljaars 1995, QJRMS)
84  pp.query("most.include_wstar", m_include_wstar);
85 
86  std::string pblh_string{"none"};
87  pp.query("most.pblh_calc", pblh_string);
88  if (pblh_string == "none") {
90  } else if (pblh_string == "MYNN2.5") {
92  } else if (pblh_string == "YSU") {
94  } else {
95  amrex::Abort("Undefined PBLH calc type!");
96  }
97 
98  // Get surface temperature
99  auto erf_st = pp.query("most.surf_temp", surf_temp);
100 
101  // Custom type user must specify the fluxes
104  pp.query("most.ustar", custom_ustar);
105  pp.query("most.tstar", custom_tstar);
106  pp.query("most.qstar", custom_qstar);
107  if (custom_qstar != 0) {
108  AMREX_ASSERT_WITH_MESSAGE(use_moisture, "Specified custom MOST qv flux without moisture model!");
109  }
110  amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
111  << custom_ustar << " "
112  << custom_tstar << " "
113  << custom_qstar << std::endl;
114 
115  // Specify surface temperature or surface flux
116  } else {
117  if (erf_st) {
119  pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
120  surf_heating_rate = surf_heating_rate / 3600.0; // [K/s]
121  if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
122  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
123  }
124  } else {
125  pp.query("most.surf_temp_flux", surf_temp_flux);
126  if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
127  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
128  }
129  if (std::abs(surf_temp_flux) > std::numeric_limits<amrex::Real>::epsilon()) {
131  } else {
133  }
134  }
135  }
136 
137  // Specify how to compute the surface flux over land (if there is any)
138  std::string rough_land_string{"constant"};
139  pp.query("most.roughness_type_land", rough_land_string);
140  if (rough_land_string == "constant") {
142  } else {
143  amrex::Abort("Undefined MOST roughness type for land!");
144  }
145 
146  // Specify how to compute the surface flux over sea (if there is any)
147  std::string rough_sea_string{"charnock"};
148  pp.query("most.roughness_type_sea", rough_sea_string);
149  if (rough_sea_string == "charnock") {
151  pp.query("most.charnock_constant",cnk_a);
152  pp.query("most.charnock_viscosity",cnk_visc);
153  if (cnk_a > 0) {
154  amrex::Print() << "If there is water, Charnock relation with C_a=" << cnk_a
155  << (cnk_visc? " and viscosity" : "")
156  << " will be used"
157  << std::endl;
158  } else {
159  amrex::Print() << "If there is water, Charnock relation with variable Charnock parameter (COARE3.0)"
160  << (cnk_visc? " and viscosity" : "")
161  << " will be used"
162  << std::endl;
163  }
164  } else if (rough_sea_string == "coare3.0") {
166  amrex::Print() << "If there is water, Charnock relation with variable Charnock parameter (COARE3.0)"
167  << (cnk_visc? " and viscosity" : "")
168  << " will be used"
169  << std::endl;
170  cnk_a = -1;
171  } else if (rough_sea_string == "donelan") {
173  } else if (rough_sea_string == "modified_charnock") {
175  pp.query("most.modified_charnock_depth",depth);
176  } else if (rough_sea_string == "wave_coupled") {
178  } else if (rough_sea_string == "constant") {
180  } else {
181  amrex::Abort("Undefined MOST roughness type for sea!");
182  }
183 
184  // Check if there is a user-specified roughness file to be read
185  std::string fname;
186  bool read_z0 = false;
187  if ( (flux_type == FluxCalcType::MOENG) ||
189  read_z0 = pp.query("most.roughness_file_name", fname);
190  }
191 
192  // Do we have a time-varying surface roughness that needs to be saved?
193  const int nghost = 0; // ghost cells not included
194  const int level = 0;
195  int lmask_min = lmask_min_reduce(*lmask_lev[level][0].get(), nghost);
196  amrex::ParallelDescriptor::ReduceIntMin(lmask_min);
197 
198  m_var_z0 = (lmask_min < 1) & (rough_type_sea != RoughCalcType::CONSTANT);
199  if (m_var_z0) {
200  amrex::Print() << "Variable sea roughness (type " << rough_sea_string << ")" << std::endl;
201  }
202 
203  // Size the MOST params for all levels
204  int nlevs = m_geom.size();
205  z_0.resize(nlevs);
206  u_star.resize(nlevs);
207  w_star.resize(nlevs);
208  t_star.resize(nlevs);
209  q_star.resize(nlevs);
210  t_surf.resize(nlevs);
211  olen.resize(nlevs);
212  pblh.resize(nlevs);
213 
214  // Get pointers to SST and LANDMASK data
215  m_sst_lev.resize(nlevs);
216  m_lmask_lev.resize(nlevs);
217 
218  for (int lev(0); lev<nlevs; ++lev) {
219  int nt_tot_sst = sst_lev[lev].size();
220  m_sst_lev[lev].resize(nt_tot_sst);
221  for (int nt(0); nt<nt_tot_sst; ++nt) {
222  m_sst_lev[lev][nt] = sst_lev[lev][nt].get();
223  }
224  int nt_tot_lmask = lmask_lev[lev].size();
225  m_lmask_lev[lev].resize(nt_tot_lmask);
226  for (int nt(0); nt<nt_tot_lmask; ++nt) {
227  m_lmask_lev[lev][nt] = lmask_lev[lev][nt].get();
228  }
229  } // lev
230 
231  // Get pointers to LSM data and Fluxes
232  m_lsm_data_lev.resize(nlevs);
233  m_lsm_flux_lev.resize(nlevs);
234  for (int lev(0); lev<nlevs; ++lev) {
235  int nvar = lsm_data[lev].size();
236  m_lsm_data_lev[lev].resize(nvar);
237  m_lsm_flux_lev[lev].resize(nvar);
238  for (int n(0); n<nvar; ++n) {
239  m_lsm_data_lev[lev][n] = lsm_data[lev][n];
240  m_lsm_flux_lev[lev][n] = lsm_flux[lev][n];
241  }
242  } // lev
243 
244  // Get pointers to wave data
245  m_Hwave_lev.resize(nlevs);
246  m_Lwave_lev.resize(nlevs);
247  m_eddyDiffs_lev.resize(nlevs);
248  for (int lev(0); lev<nlevs; ++lev) {
249  m_Hwave_lev[lev] = Hwave[lev].get();
250  m_Lwave_lev[lev] = Lwave[lev].get();
251  m_eddyDiffs_lev[lev] = eddyDiffs[lev].get();
252  }
253 
254  for (int lev = 0; lev < nlevs; lev++) {
255  // Attributes for MFs and FABs
256  //--------------------------------------------------------
257  auto& mf = vars_old[lev][Vars::cons];
258  // Create a 2D ba, dm, & ghost cells
259  amrex::BoxArray ba = mf.boxArray();
260  amrex::BoxList bl2d = ba.boxList();
261  for (auto& b : bl2d) {
262  b.setRange(2,0);
263  }
264  amrex::BoxArray ba2d(std::move(bl2d));
265  const amrex::DistributionMapping& dm = mf.DistributionMap();
266  const int ncomp = 1;
267  amrex::IntVect ng = mf.nGrowVect(); ng[2]=0;
268 
269  // Z0 heights FAB
270  //--------------------------------------------------------
271  amrex::Arena* Arena_Used = amrex::The_Arena();
272 #ifdef AMREX_USE_GPU
273  Arena_Used = amrex::The_Pinned_Arena();
274 #endif
275  amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng);
276  bx.setSmall(2,0);
277  bx.setBig(2,0);
278  z_0[lev].resize(bx,1,Arena_Used);
279  z_0[lev].setVal<amrex::RunOn::Device>(z0_const);
280  if (read_z0) read_custom_roughness(lev,fname);
281 
282  // 2D MFs for U*, T*, T_surf
283  //--------------------------------------------------------
284  u_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
285  u_star[lev]->setVal(1.E34);
286 
287  w_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
288  w_star[lev]->setVal(1.E34);
289 
290  t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
291  t_star[lev]->setVal(0.0); // default to neutral
292 
293  q_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
294  q_star[lev]->setVal(0.0); // default to dry
295 
296  olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
297  olen[lev]->setVal(1.E34);
298 
299  pblh[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
300  pblh[lev]->setVal(1.E34);
301 
302  t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
303 
304  // TODO: Do we want lsm_data to have theta at 0 index always?
305  // Do we want an external enum struct for indexing?
306  if (m_sst_lev[lev][0] || m_lsm_data_lev[lev][0]) {
307  // Valid SST or LSM data; t_surf set before computing fluxes (avoids extended lambda capture)
309  } else if (erf_st) {
310  t_surf[lev]->setVal(surf_temp);
311  } else {
312  t_surf[lev]->setVal(0.0);
313  }
314  }// lev
315  }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:219
amrex::Real m_bdy_time_interval
Definition: ERF_ABLMost.H:521
amrex::Vector< amrex::MultiFab * > m_eddyDiffs_lev
Definition: ERF_ABLMost.H:541
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_star
Definition: ERF_ABLMost.H:530
amrex::Vector< std::unique_ptr< amrex::MultiFab > > olen
Definition: ERF_ABLMost.H:531
amrex::Vector< std::unique_ptr< amrex::MultiFab > > w_star
Definition: ERF_ABLMost.H:528
amrex::Vector< amrex::FArrayBox > z_0
Definition: ERF_ABLMost.H:523
amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_ABLMost.H:522
amrex::Real custom_ustar
Definition: ERF_ABLMost.H:514
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_surf
Definition: ERF_ABLMost.H:533
amrex::Real custom_qstar
Definition: ERF_ABLMost.H:516
amrex::Real depth
Definition: ERF_ABLMost.H:519
amrex::Real surf_heating_rate
Definition: ERF_ABLMost.H:512
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_flux_lev
Definition: ERF_ABLMost.H:538
MOSTAverage m_ma
Definition: ERF_ABLMost.H:526
@ SURFACE_TEMPERATURE
Surface temperature specified.
@ HEAT_FLUX
Heat-flux specified.
int lmask_min_reduce(amrex::iMultiFab &lmask, const int &nghost)
Definition: ERF_ABLMost.H:455
bool cnk_visc
Definition: ERF_ABLMost.H:518
amrex::Real m_start_bdy_time
Definition: ERF_ABLMost.H:520
amrex::Vector< std::unique_ptr< amrex::MultiFab > > pblh
Definition: ERF_ABLMost.H:532
amrex::Vector< amrex::MultiFab * > m_Hwave_lev
Definition: ERF_ABLMost.H:539
void read_custom_roughness(const int &lev, const std::string &fname)
Definition: ERF_ABLMost.cpp:592
amrex::Real surf_temp
Definition: ERF_ABLMost.H:511
amrex::Real cnk_a
Definition: ERF_ABLMost.H:517
ThetaCalcType theta_type
Definition: ERF_ABLMost.H:500
PBLHeightCalcType pblh_type
Definition: ERF_ABLMost.H:503
@ MOENG
Moeng functional form.
@ CUSTOM
Custom constant flux functional form.
@ ROTATE
Terrain rotation flux functional form.
@ DONELAN
Donelan functional form.
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_data_lev
Definition: ERF_ABLMost.H:537
amrex::Real surf_temp_flux
Definition: ERF_ABLMost.H:513
amrex::Vector< amrex::Vector< amrex::iMultiFab * > > m_lmask_lev
Definition: ERF_ABLMost.H:536
RoughCalcType rough_type_land
Definition: ERF_ABLMost.H:501
bool m_rotate
Definition: ERF_ABLMost.H:508
amrex::Vector< amrex::MultiFab * > m_Lwave_lev
Definition: ERF_ABLMost.H:540
bool m_var_z0
Definition: ERF_ABLMost.H:524
RoughCalcType rough_type_sea
Definition: ERF_ABLMost.H:502
bool use_moisture
Definition: ERF_ABLMost.H:506
bool m_include_wstar
Definition: ERF_ABLMost.H:509
bool m_exp_most
Definition: ERF_ABLMost.H:507
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_star
Definition: ERF_ABLMost.H:529
amrex::Vector< std::unique_ptr< amrex::MultiFab > > u_star
Definition: ERF_ABLMost.H:527
amrex::Real z0_const
Definition: ERF_ABLMost.H:510
FluxCalcType flux_type
Definition: ERF_ABLMost.H:499
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_sst_lev
Definition: ERF_ABLMost.H:535
@ CONSTANT
Constant z0.
amrex::Real custom_tstar
Definition: ERF_ABLMost.H:515
@ cons
Definition: ERF_IndexDefines.H:129
Here is the call graph for this function:

Member Function Documentation

◆ compute_fluxes()

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

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

Parameters
[in]levCurrent level
[in]max_itersmaximum iterations to use
[in]most_fluxstructure to iteratively compute ustar and tstar
149 {
150  // Pointers to the computed averages
151  const auto *const tm_ptr = m_ma.get_average(lev,2); // potential temperature
152  const auto *const qvm_ptr = m_ma.get_average(lev,3); // water vapor mixing ratio
153  const auto *const tvm_ptr = m_ma.get_average(lev,4); // virtual potential temperature
154  const auto *const umm_ptr = m_ma.get_average(lev,5); // horizontal velocity magnitude
155 
156  for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
157  {
158  Box gtbx = mfi.growntilebox();
159 
160  auto u_star_arr = u_star[lev]->array(mfi);
161  auto t_star_arr = t_star[lev]->array(mfi);
162  auto q_star_arr = q_star[lev]->array(mfi);
163  auto t_surf_arr = t_surf[lev]->array(mfi);
164  auto olen_arr = olen[lev]->array(mfi);
165 
166  const auto tm_arr = tm_ptr->array(mfi);
167  const auto tvm_arr = tvm_ptr->array(mfi);
168  const auto qvm_arr = qvm_ptr->array(mfi);
169  const auto umm_arr = umm_ptr->array(mfi);
170  const auto z0_arr = z_0[lev].array();
171 
172  // PBL height if we need to calculate wstar for the Beljaars correction
173  // TODO: can/should we apply this in LES mode?
174  const auto w_star_arr = (m_include_wstar) ? w_star[lev].get()->array(mfi) : Array4<Real> {};
175  const auto pblh_arr = (m_include_wstar) ? pblh[lev].get()->array(mfi) : Array4<Real> {};
176 
177  // Wave properties if they exist
178  const auto Hwave_arr = (m_Hwave_lev[lev]) ? m_Hwave_lev[lev]->array(mfi) : Array4<Real> {};
179  const auto Lwave_arr = (m_Lwave_lev[lev]) ? m_Lwave_lev[lev]->array(mfi) : Array4<Real> {};
180  const auto eta_arr = (m_eddyDiffs_lev[lev]) ? m_eddyDiffs_lev[lev]->array(mfi) : Array4<Real> {};
181 
182  // Land mask array if it exists
183  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
184  Array4<int> {};
185 
186  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
187  {
188  if (( is_land && lmask_arr(i,j,k) == 1) ||
189  (!is_land && lmask_arr(i,j,k) == 0))
190  {
191  most_flux.iterate_flux(i, j, k, max_iters,
192  z0_arr, umm_arr, tm_arr, tvm_arr, qvm_arr,
193  u_star_arr, w_star_arr, // to be updated
194  t_star_arr, q_star_arr, // to be updated
195  t_surf_arr, olen_arr, // to be updated
196  pblh_arr, Hwave_arr, Lwave_arr, eta_arr);
197  }
198  });
199  }
200 }
const amrex::MultiFab * get_average(int lev, int comp) const
Definition: ERF_MOSTAverage.H:92

◆ compute_most_bcs() [1/2]

template<typename FluxCalc >
void ABLMost::compute_most_bcs ( const int &  lev,
const amrex::Vector< amrex::MultiFab * > &  mfs,
amrex::MultiFab *  xxmom_flux,
amrex::MultiFab *  yymom_flux,
amrex::MultiFab *  zzmom_flux,
amrex::MultiFab *  xymom_flux,
amrex::MultiFab *  yxmom_flux,
amrex::MultiFab *  xzmom_flux,
amrex::MultiFab *  zxmom_flux,
amrex::MultiFab *  yzmom_flux,
amrex::MultiFab *  zymom_flux,
amrex::MultiFab *  xheat_flux,
amrex::MultiFab *  yheat_flux,
amrex::MultiFab *  zheat_flux,
amrex::MultiFab *  xqv_flux,
amrex::MultiFab *  yqv_flux,
amrex::MultiFab *  zqv_flux,
amrex::MultiFab *  z_phys,
const FluxCalc &  flux_comp 
)

◆ compute_most_bcs() [2/2]

template<typename FluxCalc >
void ABLMost::compute_most_bcs ( const int &  lev,
const Vector< MultiFab * > &  mfs,
MultiFab *  xxmom_flux,
MultiFab *  yymom_flux,
MultiFab *  zzmom_flux,
MultiFab *  xymom_flux,
MultiFab *  yxmom_flux,
MultiFab *  xzmom_flux,
MultiFab *  zxmom_flux,
MultiFab *  yzmom_flux,
MultiFab *  zymom_flux,
MultiFab *  xheat_flux,
MultiFab *  yheat_flux,
MultiFab *  zheat_flux,
MultiFab *  xqv_flux,
MultiFab *  yqv_flux,
MultiFab *  zqv_flux,
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
306 {
307  const int klo = 0;
308  const int icomp = 0;
309  const auto& dxInv = m_geom[lev].InvCellSizeArray();
310  const auto& dz_no_terrain = m_geom[lev].CellSize(2);
311  for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
312  {
313  // TODO: No LSM lateral ghost cells, should this change?
314  // Valid CC box
315  Box vbx = mfi.validbox(); vbx.makeSlab(2,klo-1);
316 
317  Box vbxx = surroundingNodes(vbx,0);
318  Box vbxy = surroundingNodes(vbx,1);
319 
320  // Expose for GPU
321  bool exp_most = m_exp_most;
322  bool rot_most = m_rotate;
323 
324  // Get field arrays
325  const auto cons_arr = mfs[Vars::cons]->array(mfi);
326  const auto velx_arr = mfs[Vars::xvel]->array(mfi);
327  const auto vely_arr = mfs[Vars::yvel]->array(mfi);
328  const auto velz_arr = mfs[Vars::zvel]->array(mfi);
329 
330  // Explicit MOST vars
331  auto t13_arr = (m_exp_most) ? xzmom_flux->array(mfi) : Array4<Real>{};
332  auto t31_arr = (m_exp_most && zxmom_flux) ? zxmom_flux->array(mfi) : Array4<Real>{};
333 
334  auto t23_arr = (m_exp_most) ? yzmom_flux->array(mfi) : Array4<Real>{};
335  auto t32_arr = (m_exp_most && zymom_flux) ? zymom_flux->array(mfi) : Array4<Real>{};
336 
337  auto hfx3_arr = (m_exp_most) ? zheat_flux->array(mfi) : Array4<Real>{};
338  auto qfx3_arr = (m_exp_most && zqv_flux) ? zqv_flux->array(mfi) : Array4<Real>{};
339 
340  // Rotated MOST vars
341  auto t11_arr = (m_rotate) ? xxmom_flux->array(mfi) : Array4<Real>{};
342  auto t22_arr = (m_rotate) ? yymom_flux->array(mfi) : Array4<Real>{};
343  auto t33_arr = (m_rotate) ? zzmom_flux->array(mfi) : Array4<Real>{};
344  auto t12_arr = (m_rotate) ? xymom_flux->array(mfi) : Array4<Real>{};
345  auto t21_arr = (m_rotate) ? yxmom_flux->array(mfi) : Array4<Real>{};
346 
347  auto hfx1_arr = (m_rotate) ? xheat_flux->array(mfi) : Array4<Real>{};
348  auto hfx2_arr = (m_rotate) ? yheat_flux->array(mfi) : Array4<Real>{};
349  auto qfx1_arr = (m_rotate && xqv_flux) ? xqv_flux->array(mfi) : Array4<Real>{};
350  auto qfx2_arr = (m_rotate && yqv_flux) ? yqv_flux->array(mfi) : Array4<Real>{};
351 
352  // Viscosity and terrain
353  const auto eta_arr = (!m_exp_most) ? m_eddyDiffs_lev[lev]->array(mfi) : Array4<const Real>{};
354  const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};
355 
356  // Get average arrays
357  const auto *const u_mean = m_ma.get_average(lev,0);
358  const auto *const v_mean = m_ma.get_average(lev,1);
359  const auto *const t_mean = m_ma.get_average(lev,2);
360  const auto *const q_mean = m_ma.get_average(lev,3);
361  const auto *const u_mag_mean = m_ma.get_average(lev,5);
362 
363  const auto um_arr = u_mean->array(mfi);
364  const auto vm_arr = v_mean->array(mfi);
365  const auto tm_arr = t_mean->array(mfi);
366  const auto qm_arr = q_mean->array(mfi);
367  const auto umm_arr = u_mag_mean->array(mfi);
368 
369  // Get derived arrays
370  const auto u_star_arr = u_star[lev]->array(mfi);
371  const auto t_star_arr = t_star[lev]->array(mfi);
372  const auto q_star_arr = q_star[lev]->array(mfi);
373  const auto t_surf_arr = t_surf[lev]->array(mfi);
374 
375  // Get LSM fluxes
376  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
377  Array4<int> {};
378  auto lsm_flux_arr = (m_lsm_flux_lev[lev][0]) ? m_lsm_flux_lev[lev][0]->array(mfi) :
379  Array4<Real> {};
380 
381  for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx)
382  {
383  const Box& bx = (*mfs[var_idx])[mfi].box();
384  auto dest_arr = (*mfs[var_idx])[mfi].array();
385 
386  if (var_idx == Vars::cons) {
387  Box b2d = bx;
388  b2d.setBig(2,klo-1);
389  int n = RhoTheta_comp;
390  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
391  {
392  Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
393  Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain;
394  Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr,
395  cons_arr, velx_arr, vely_arr,
396  umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr,
397  dest_arr);
398 
399  if (rot_most) {
400  rotate_scalar_flux(i, j, klo, Tflux, dxInv, zphys_arr,
401  hfx1_arr, hfx2_arr, hfx3_arr);
402  } else {
403  int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
404  if (is_land && lsm_flux_arr && vbx.contains(i,j,k)) {
405  lsm_flux_arr(i,j,klo) = Tflux;
406  }
407  else if ((k == klo-1) && vbx.contains(i,j,k) && exp_most) {
408  hfx3_arr(i,j,klo) = Tflux;
409  }
410  }
411  });
412 
413  // TODO: Generalize MOST q flux
415  n = RhoQ1_comp;
416  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
417  {
418  Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
419  Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain;
420  Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr,
421  cons_arr, velx_arr, vely_arr,
422  umm_arr, qm_arr, u_star_arr, q_star_arr, t_surf_arr,
423  dest_arr);
424 
425  if (rot_most) {
426  rotate_scalar_flux(i, j, klo, Qflux, dxInv, zphys_arr,
427  qfx1_arr, qfx2_arr, qfx3_arr);
428  } else {
429  if ((k == klo-1) && vbx.contains(i,j,k) && exp_most) {
430  qfx3_arr(i,j,klo) = Qflux;
431  }
432  }
433  });
434  }
435 
436  } else if (var_idx == Vars::xvel) {
437 
438  Box xb2d = surroundingNodes(bx,0);
439  xb2d.setBig(2,klo-1);
440 
441  ParallelFor(xb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
442  {
443  Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
444  Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain;
445  Real stressx = flux_comp.compute_u_flux(i, j, k, icomp, dz, dz1, exp_most, eta_arr,
446  cons_arr, velx_arr, vely_arr,
447  umm_arr, um_arr, u_star_arr,
448  dest_arr);
449 
450  if (rot_most) {
451  rotate_stress_tensor(i, j, klo, stressx, dxInv, zphys_arr,
452  velx_arr, vely_arr, velz_arr,
453  t11_arr, t22_arr, t33_arr,
454  t12_arr, t21_arr,
455  t13_arr, t31_arr,
456  t23_arr, t32_arr);
457  } else {
458  if ((k == klo-1) && vbxx.contains(i,j,k) && exp_most) {
459  t13_arr(i,j,klo) = stressx;
460  if (t31_arr) t31_arr(i,j,klo) = stressx;
461  }
462  }
463  });
464 
465  } else if (var_idx == Vars::yvel) {
466 
467  Box yb2d = surroundingNodes(bx,1);
468  yb2d.setBig(2,klo-1);
469 
470  ParallelFor(yb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
471  {
472  Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
473  Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain;
474  Real stressy = flux_comp.compute_v_flux(i, j, k, icomp, dz, dz1, exp_most, eta_arr,
475  cons_arr, velx_arr, vely_arr,
476  umm_arr, vm_arr, u_star_arr,
477  dest_arr);
478 
479  // NOTE: One stress rotation for ALL the stress components
480  if (!rot_most) {
481  if ((k == klo-1) && vbxy.contains(i,j,k) && exp_most) {
482  t23_arr(i,j,klo) = stressy;
483  if (t32_arr) t32_arr(i,j,klo) = stressy;
484  }
485  }
486  });
487  }
488  } // var_idx
489  } // mfiter
490 }
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
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:454
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:475
@ xvel
Definition: ERF_IndexDefines.H:130
@ zvel
Definition: ERF_IndexDefines.H:132
@ NumTypes
Definition: ERF_IndexDefines.H:133
@ yvel
Definition: ERF_IndexDefines.H:131
Here is the call graph for this function:

◆ compute_pblh() [1/2]

template<typename PBLHeightEstimator >
void ABLMost::compute_pblh ( const int &  lev,
amrex::Vector< amrex::Vector< amrex::MultiFab >> &  vars,
amrex::MultiFab *  z_phys_cc,
const PBLHeightEstimator &  est,
const int  RhoQv_comp,
const int  RhoQc_comp,
const int  RhoQr_comp 
)

◆ compute_pblh() [2/2]

template<typename PBLHeightEstimator >
void ABLMost::compute_pblh ( const int &  lev,
Vector< Vector< MultiFab >> &  vars,
MultiFab *  z_phys_cc,
const PBLHeightEstimator &  est,
int  RhoQv_comp,
int  RhoQc_comp,
int  RhoQr_comp 
)
585 {
586  est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
587  vars[lev][Vars::cons],m_lmask_lev[lev][0],
588  RhoQv_comp, RhoQc_comp, RhoQr_comp);
589 }

◆ get_lmask()

const amrex::iMultiFab* ABLMost::get_lmask ( const int &  lev)
inline
452 { return m_lmask_lev[lev][0]; }

◆ get_lsm_tsurf()

void ABLMost::get_lsm_tsurf ( const int &  lev)
529 {
530  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
531  {
532  Box gtbx = mfi.growntilebox();
533 
534  // TODO: LSM does not carry lateral ghost cells.
535  // This copies the valid box into the ghost cells.
536  // Fillboundary is called after this to pick up the
537  // interior ghost and periodic directions. Is there
538  // a better approach?
539  Box vbx = mfi.validbox();
540  int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
541  int j_lo = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1);
542 
543  auto t_surf_arr = t_surf[lev]->array(mfi);
544  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
545  Array4<int> {};
546  const auto lsm_arr = m_lsm_data_lev[lev][0]->const_array(mfi);
547 
548  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
549  {
550  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
551  if (is_land) {
552  int li = amrex::min(amrex::max(i, i_lo), i_hi);
553  int lj = amrex::min(amrex::max(j, j_lo), j_hi);
554  t_surf_arr(i,j,k) = lsm_arr(li,lj,k);
555  }
556  });
557  }
558 }

◆ get_mac_avg()

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

◆ get_olen()

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

◆ get_pblh()

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

◆ get_q_star()

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

◆ get_t_star()

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

◆ get_t_surf()

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

◆ get_u_star()

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

◆ get_w_star()

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

◆ get_z0()

const amrex::FArrayBox* ABLMost::get_z0 ( const int &  lev)
inline
447 { return &z_0[lev]; }

◆ get_zref()

amrex::Real ABLMost::get_zref ( )
inline
444 { return m_ma.get_zref(); }
amrex::Real get_zref() const
Definition: ERF_MOSTAverage.H:95
Here is the call graph for this function:

◆ have_variable_sea_roughness()

bool ABLMost::have_variable_sea_roughness ( )
inline
449 { return m_var_z0; }

◆ impose_most_bcs()

void ABLMost::impose_most_bcs ( const int &  lev,
const amrex::Vector< amrex::MultiFab * > &  mfs,
amrex::MultiFab *  xxmom_flux,
amrex::MultiFab *  yymom_flux,
amrex::MultiFab *  zzmom_flux,
amrex::MultiFab *  xymom_flux,
amrex::MultiFab *  yxmom_flux,
amrex::MultiFab *  xzmom_flux,
amrex::MultiFab *  zxmom_flux,
amrex::MultiFab *  yzmom_flux,
amrex::MultiFab *  zymom_flux,
amrex::MultiFab *  xheat_flux,
amrex::MultiFab *  yheat_flux,
amrex::MultiFab *  zheat_flux,
amrex::MultiFab *  xqv_flux,
amrex::MultiFab *  yqv_flux,
amrex::MultiFab *  zqv_flux,
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
226 {
227  const int klo = 0;
229  moeng_flux flux_comp(klo);
230  compute_most_bcs(lev, mfs,
231  xxmom_flux,
232  yymom_flux,
233  zzmom_flux,
234  xymom_flux, yxmom_flux,
235  xzmom_flux, zxmom_flux,
236  yzmom_flux, zymom_flux,
237  xheat_flux, yheat_flux, zheat_flux,
238  xqv_flux, yqv_flux, zqv_flux,
239  z_phys, flux_comp);
240  } else if (flux_type == FluxCalcType::DONELAN) {
241  donelan_flux flux_comp(klo);
242  compute_most_bcs(lev, mfs,
243  xxmom_flux,
244  yymom_flux,
245  zzmom_flux,
246  xymom_flux, yxmom_flux,
247  xzmom_flux, zxmom_flux,
248  yzmom_flux, zymom_flux,
249  xheat_flux, yheat_flux, zheat_flux,
250  xqv_flux, yqv_flux, zqv_flux,
251  z_phys, flux_comp);
252  } else if (flux_type == FluxCalcType::ROTATE) {
253  rotate_flux flux_comp(klo);
254  compute_most_bcs(lev, mfs,
255  xxmom_flux,
256  yymom_flux,
257  zzmom_flux,
258  xymom_flux, yxmom_flux,
259  xzmom_flux, zxmom_flux,
260  yzmom_flux, zymom_flux,
261  xheat_flux, yheat_flux, zheat_flux,
262  xqv_flux, yqv_flux, zqv_flux,
263  z_phys, flux_comp);
264  } else {
265  custom_flux flux_comp(klo);
266  compute_most_bcs(lev, mfs,
267  xxmom_flux,
268  yymom_flux,
269  zzmom_flux,
270  xymom_flux, yxmom_flux,
271  xzmom_flux, zxmom_flux,
272  yzmom_flux, zymom_flux,
273  xheat_flux, yheat_flux, zheat_flux,
274  xqv_flux, yqv_flux, zqv_flux,
275  z_phys, flux_comp);
276  }
277 }
void compute_most_bcs(const int &lev, const amrex::Vector< amrex::MultiFab * > &mfs, amrex::MultiFab *xxmom_flux, amrex::MultiFab *yymom_flux, amrex::MultiFab *zzmom_flux, amrex::MultiFab *xymom_flux, amrex::MultiFab *yxmom_flux, amrex::MultiFab *xzmom_flux, amrex::MultiFab *zxmom_flux, amrex::MultiFab *yzmom_flux, amrex::MultiFab *zymom_flux, amrex::MultiFab *xheat_flux, amrex::MultiFab *yheat_flux, amrex::MultiFab *zheat_flux, amrex::MultiFab *xqv_flux, amrex::MultiFab *yqv_flux, amrex::MultiFab *zqv_flux, amrex::MultiFab *z_phys, const FluxCalc &flux_comp)
Definition: ERF_MOSTStress.H:1848
Definition: ERF_MOSTStress.H:1567
Definition: ERF_MOSTStress.H:1259
Definition: ERF_MOSTStress.H:2105

◆ lmask_min_reduce()

int ABLMost::lmask_min_reduce ( amrex::iMultiFab &  lmask,
const int &  nghost 
)
inline
456  {
457  int lmask_min = amrex::ReduceMin(lmask, nghost,
458  [=] AMREX_GPU_HOST_DEVICE (amrex::Box const& bx, amrex::Array4<int const> const& lm_arr) -> int
459  {
460  int locmin = std::numeric_limits<int>::max();
461  const auto lo = lbound(bx);
462  const auto hi = ubound(bx);
463  for (int j = lo.y; j <= hi.y; ++j) {
464  for (int i = lo.x; i <= hi.x; ++i) {
465  locmin = std::min(locmin, lm_arr(i,j,0));
466  }
467  }
468  return locmin;
469  });
470 
471  return lmask_min;
472  }

Referenced by ABLMost().

Here is the caller graph for this function:

◆ read_custom_roughness()

void ABLMost::read_custom_roughness ( const int &  lev,
const std::string &  fname 
)
594 {
595  // Read the file if we are on the coarsest level
596  if (lev==0) {
597  // Only the ioproc reads the file and broadcasts
598  if (ParallelDescriptor::IOProcessor()) {
599  Print()<<"Reading MOST roughness file: "<< fname << std::endl;
600  std::ifstream file(fname);
601  Gpu::HostVector<Real> m_x,m_y,m_z0;
602  Real value1,value2,value3;
603  while(file>>value1>>value2>>value3){
604  m_x.push_back(value1);
605  m_y.push_back(value2);
606  m_z0.push_back(value3);
607  }
608  file.close();
609 
610  // Copy data to the GPU
611  int nnode = m_x.size();
612  Gpu::DeviceVector<Real> d_x(nnode),d_y(nnode),d_z0(nnode);
613  Gpu::copy(Gpu::hostToDevice, m_x.begin(), m_x.end(), d_x.begin());
614  Gpu::copy(Gpu::hostToDevice, m_y.begin(), m_y.end(), d_y.begin());
615  Gpu::copy(Gpu::hostToDevice, m_z0.begin(), m_z0.end(), d_z0.begin());
616  Real* xp = d_x.data();
617  Real* yp = d_y.data();
618  Real* z0p = d_z0.data();
619 
620  // Populate z_phys data
621  Real tol = 1.0e-4;
622  auto dx = m_geom[lev].CellSizeArray();
623  auto ProbLoArr = m_geom[lev].ProbLoArray();
624  int ilo = m_geom[lev].Domain().smallEnd(0);
625  int jlo = m_geom[lev].Domain().smallEnd(1);
626  int klo = 0;
627  int ihi = m_geom[lev].Domain().bigEnd(0);
628  int jhi = m_geom[lev].Domain().bigEnd(1);
629 
630  // Grown box with no z range
631  Box xybx = z_0[lev].box();
632  xybx.setRange(2,0);
633 
634  Array4<Real> const& z0_arr = z_0[lev].array();
635  ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
636  {
637  // Clip indices for ghost-cells
638  int ii = amrex::min(amrex::max(i,ilo),ihi);
639  int jj = amrex::min(amrex::max(j,jlo),jhi);
640 
641  // Location of nodes
642  Real x = ProbLoArr[0] + ii * dx[0];
643  Real y = ProbLoArr[1] + jj * dx[1];
644  int inode = ii + jj * (ihi-ilo+2); // stride is Nx+1
645  if (std::sqrt(std::pow(x-xp[inode],2)+std::pow(y-yp[inode],2)) < tol) {
646  z0_arr(i,j,klo) = z0p[inode];
647  } else {
648  // Unexpected list order, do brute force search
649  Real z0loc = 0.0;
650  bool found = false;
651  for (int n=0; n<nnode; ++n) {
652  Real delta=std::sqrt(std::pow(x-xp[n],2)+std::pow(y-yp[n],2));
653  if (delta < tol) {
654  found = true;
655  z0loc = z0p[n];
656  break;
657  }
658  }
659  AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
660  amrex::ignore_unused(found);
661  z0_arr(i,j,klo) = z0loc;
662  }
663  });
664  } // Is ioproc
665 
666  int ioproc = ParallelDescriptor::IOProcessorNumber();
667  ParallelDescriptor::Barrier();
668  ParallelDescriptor::Bcast(z_0[lev].dataPtr(),z_0[lev].box().numPts(),ioproc);
669  } else {
670  // Create a BC mapper that uses FOEXTRAP at domain bndry
671  Vector<int> bc_lo(3,ERFBCType::foextrap);
672  Vector<int> bc_hi(3,ERFBCType::foextrap);
673  Vector<BCRec> bcr; bcr.push_back(BCRec(bc_lo.data(),bc_hi.data()));
674 
675  // Create ref ratio
676  IntVect ratio;
677  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
678  ratio[idim] = m_geom[lev].Domain().length(idim) / m_geom[0].Domain().length(idim);
679  }
680 
681  // Create interp object and interpolate from the coarsest grid
682  Interpolater* interp = &cell_cons_interp;
683  interp->interp(z_0[0] , 0,
684  z_0[lev], 0,
685  1, z_0[lev].box(),
686  ratio, m_geom[0], m_geom[lev],
687  bcr, 0, 0, RunOn::Gpu);
688  }
689 }
@ foextrap
Definition: ERF_IndexDefines.H:179

Referenced by ABLMost().

Here is the caller graph for this function:

◆ time_interp_sst()

void ABLMost::time_interp_sst ( const int &  lev,
const amrex::Real &  time 
)
495 {
496  // Time interpolation
497  Real dT = m_bdy_time_interval;
498  Real time_since_start = time - m_start_bdy_time;
499  int n_time = static_cast<int>( time_since_start / dT);
500  Real alpha = (time_since_start - n_time * dT) / dT;
501  AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
502  Real oma = 1.0 - alpha;
503  AMREX_ALWAYS_ASSERT( (n_time >= 0) && (n_time < (m_sst_lev[lev].size()-1)));
504 
505  // Populate t_surf
506  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
507  {
508  Box gtbx = mfi.growntilebox();
509 
510  auto t_surf_arr = t_surf[lev]->array(mfi);
511  const auto sst_hi_arr = m_sst_lev[lev][n_time+1]->const_array(mfi);
512  const auto sst_lo_arr = m_sst_lev[lev][n_time ]->const_array(mfi);
513  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
514  Array4<int> {};
515 
516  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
517  {
518  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
519  if (!is_land) {
520  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
521  + alpha * sst_hi_arr(i,j,k);
522  }
523  });
524  }
525 }

◆ update_fluxes()

void ABLMost::update_fluxes ( const int &  lev,
const amrex::Real &  time,
int  max_iters = 25 
)

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

Parameters
[in]levCurrent level
[in]max_itersmaximum iterations to use
15 {
16  // Update SST data if we have a valid pointer
17  if (m_sst_lev[lev][0]) time_interp_sst(lev, time);
18 
19  // TODO: we want 0 index to always be theta?
20  // Update land surface temp if we have a valid pointer
21  if (m_lsm_data_lev[lev][0]) get_lsm_tsurf(lev);
22 
23  // Fill interior ghost cells
24  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
25 
26  // Compute plane averages for all vars (regardless of flux type)
28 
29  // ***************************************************************
30  // Iterate the fluxes if moeng type
31  // First iterate over land -- the only model for surface roughness
32  // over land is RoughCalcType::CONSTANT
33  // ***************************************************************
36  bool is_land = true;
40  compute_fluxes(lev, max_iters, most_flux, is_land);
41  } else {
42  amrex::Abort("Unknown value for rough_type_land");
43  }
45  update_surf_temp(time);
48  compute_fluxes(lev, max_iters, most_flux, is_land);
49  } else {
50  amrex::Abort("Unknown value for rough_type_land");
51  }
52  } else if (theta_type == ThetaCalcType::ADIABATIC) {
54  adiabatic most_flux(m_ma.get_zref(), surf_temp_flux);
55  compute_fluxes(lev, max_iters, most_flux, is_land);
56  } else {
57  amrex::Abort("Unknown value for rough_type_land");
58  }
59  } else {
60  amrex::Abort("Unknown value for theta_type");
61  }
62  } // MOENG -- LAND
63 
64  // ***************************************************************
65  // Iterate the fluxes if moeng type
66  // Next iterate over sea -- the models for surface roughness
67  // over sea are CHARNOCK, DONELAN, MODIFIED_CHARNOCK or WAVE_COUPLED
68  // ***************************************************************
71  bool is_land = false;
75  compute_fluxes(lev, max_iters, most_flux, is_land);
78  compute_fluxes(lev, max_iters, most_flux, is_land);
79  } else if (rough_type_sea == RoughCalcType::DONELAN) {
81  compute_fluxes(lev, max_iters, most_flux, is_land);
84  compute_fluxes(lev, max_iters, most_flux, is_land);
85  } else {
86  amrex::Abort("Unknown value for rough_type_sea");
87  }
88 
90  update_surf_temp(time);
93  compute_fluxes(lev, max_iters, most_flux, is_land);
96  compute_fluxes(lev, max_iters, most_flux, is_land);
97  } else if (rough_type_sea == RoughCalcType::DONELAN) {
99  compute_fluxes(lev, max_iters, most_flux, is_land);
102  compute_fluxes(lev, max_iters, most_flux, is_land);
103  } else {
104  amrex::Abort("Unknown value for rough_type_sea");
105  }
106 
107  } else if (theta_type == ThetaCalcType::ADIABATIC) {
110  compute_fluxes(lev, max_iters, most_flux, is_land);
113  compute_fluxes(lev, max_iters, most_flux, is_land);
114  } else if (rough_type_sea == RoughCalcType::DONELAN) {
116  compute_fluxes(lev, max_iters, most_flux, is_land);
119  compute_fluxes(lev, max_iters, most_flux, is_land);
120  } else {
121  amrex::Abort("Unknown value for rough_type_sea");
122  }
123  } else {
124  amrex::Abort("Unknown value for theta_type");
125  }
126 
127  } // MOENG -- SEA
128 
130  u_star[lev]->setVal(custom_ustar);
131  t_star[lev]->setVal(custom_tstar);
132  q_star[lev]->setVal(custom_qstar);
133  }
134 }
void compute_fluxes(const int &lev, const int &max_iters, const FluxIter &most_flux, bool is_land)
Definition: ERF_ABLMost.cpp:145
void time_interp_sst(const int &lev, const amrex::Real &time)
Definition: ERF_ABLMost.cpp:493
void update_surf_temp(const amrex::Real &time)
Definition: ERF_ABLMost.H:396
void get_lsm_tsurf(const int &lev)
Definition: ERF_ABLMost.cpp:528
void compute_averages(int lev)
Definition: ERF_MOSTAverage.cpp:664
Definition: ERF_MOSTStress.H:136
Definition: ERF_MOSTStress.H:270
Definition: ERF_MOSTStress.H:207
Definition: ERF_MOSTStress.H:330
Definition: ERF_MOSTStress.H:90
Definition: ERF_MOSTStress.H:477
Definition: ERF_MOSTStress.H:653
Definition: ERF_MOSTStress.H:569
Definition: ERF_MOSTStress.H:734
Definition: ERF_MOSTStress.H:399
Definition: ERF_MOSTStress.H:904
Definition: ERF_MOSTStress.H:1084
Definition: ERF_MOSTStress.H:998
Definition: ERF_MOSTStress.H:1167
Definition: ERF_MOSTStress.H:824

◆ update_mac_ptrs()

void ABLMost::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
417  { m_ma.update_field_ptrs(lev,vars_old,Theta_prim,Qv_prim,Qr_prim); }
void update_field_ptrs(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:228
Here is the call graph for this function:

◆ update_pblh()

void ABLMost::update_pblh ( const int &  lev,
amrex::Vector< amrex::Vector< amrex::MultiFab >> &  vars,
amrex::MultiFab *  z_phys_cc,
const int  RhoQv_comp,
const int  RhoQc_comp,
const int  RhoQr_comp 
)
567 {
569  MYNNPBLH estimator;
570  compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQc_comp, RhoQr_comp);
571  } else if (pblh_type == PBLHeightCalcType::YSU) {
572  amrex::Error("YSU PBLH calc not implemented yet");
573  }
574 }
void compute_pblh(const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars, amrex::MultiFab *z_phys_cc, const PBLHeightEstimator &est, const int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp)
Definition: ERF_PBLHeight.H:8

◆ update_surf_temp()

void ABLMost::update_surf_temp ( const amrex::Real &  time)
inline
397  {
398  if (surf_heating_rate != 0) {
399  int nlevs = m_geom.size();
400  for (int lev = 0; lev < nlevs; lev++)
401  {
402  t_surf[lev]->setVal(surf_temp + surf_heating_rate*time);
403  amrex::Print() << "Surface temp at t=" << time
404  << ": "
405  << surf_temp + surf_heating_rate*time
406  << std::endl;
407  }
408  }
409  }

Member Data Documentation

◆ cnk_a

amrex::Real ABLMost::cnk_a {0.0185}
private

Referenced by ABLMost().

◆ cnk_visc

bool ABLMost::cnk_visc {false}
private

Referenced by ABLMost().

◆ custom_qstar

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

Referenced by ABLMost().

◆ custom_tstar

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

Referenced by ABLMost().

◆ custom_ustar

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

Referenced by ABLMost().

◆ depth

amrex::Real ABLMost::depth {30.0}
private

Referenced by ABLMost().

◆ flux_type

FluxCalcType ABLMost::flux_type {FluxCalcType::MOENG}

Referenced by ABLMost().

◆ m_bdy_time_interval

amrex::Real ABLMost::m_bdy_time_interval
private

◆ m_eddyDiffs_lev

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

Referenced by ABLMost().

◆ m_exp_most

bool ABLMost::m_exp_most = false
private

◆ m_geom

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

Referenced by ABLMost(), and update_surf_temp().

◆ m_Hwave_lev

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

Referenced by ABLMost().

◆ m_include_wstar

bool ABLMost::m_include_wstar = false
private

Referenced by ABLMost().

◆ m_lmask_lev

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

Referenced by ABLMost(), and get_lmask().

◆ m_lsm_data_lev

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

Referenced by ABLMost().

◆ m_lsm_flux_lev

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

Referenced by ABLMost().

◆ m_Lwave_lev

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

Referenced by ABLMost().

◆ m_ma

MOSTAverage ABLMost::m_ma
private

◆ m_rotate

bool ABLMost::m_rotate = false
private

◆ m_sst_lev

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

Referenced by ABLMost().

◆ m_start_bdy_time

amrex::Real ABLMost::m_start_bdy_time
private

◆ m_var_z0

bool ABLMost::m_var_z0 {false}
private

◆ olen

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

Referenced by ABLMost(), and get_olen().

◆ pblh

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

Referenced by ABLMost(), and get_pblh().

◆ pblh_type

Referenced by ABLMost().

◆ q_star

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

Referenced by ABLMost(), and get_q_star().

◆ rough_type_land

RoughCalcType ABLMost::rough_type_land {RoughCalcType::CONSTANT}

Referenced by ABLMost().

◆ rough_type_sea

RoughCalcType ABLMost::rough_type_sea {RoughCalcType::CHARNOCK}

Referenced by ABLMost().

◆ surf_heating_rate

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

Referenced by ABLMost(), and update_surf_temp().

◆ surf_temp

amrex::Real ABLMost::surf_temp
private

Referenced by ABLMost(), and update_surf_temp().

◆ surf_temp_flux

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

Referenced by ABLMost().

◆ t_star

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

Referenced by ABLMost(), and get_t_star().

◆ t_surf

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

◆ theta_type

ThetaCalcType ABLMost::theta_type {ThetaCalcType::ADIABATIC}

Referenced by ABLMost().

◆ u_star

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

Referenced by ABLMost(), and get_u_star().

◆ use_moisture

bool ABLMost::use_moisture
private

Referenced by ABLMost().

◆ w_star

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

Referenced by ABLMost(), and get_w_star().

◆ z0_const

amrex::Real ABLMost::z0_const {0.1}
private

Referenced by ABLMost().

◆ z_0

amrex::Vector<amrex::FArrayBox> ABLMost::z_0
private

Referenced by ABLMost(), and get_z0().


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