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
}
 
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, amrex::Real start_time, amrex::Real stop_time, amrex::Real bdy_time_interval=0.0)
 
void make_SurfaceLayer_at_level (const int &lev, int nlevs, const amrex::Vector< amrex::MultiFab * > &mfv, std::unique_ptr< amrex::MultiFab > &Theta_prim, std::unique_ptr< amrex::MultiFab > &Qv_prim, std::unique_ptr< amrex::MultiFab > &Qr_prim, std::unique_ptr< amrex::MultiFab > &z_phys_nd, amrex::MultiFab *Hwave, amrex::MultiFab *Lwave, amrex::MultiFab *eddyDiffs, amrex::Vector< amrex::MultiFab * > lsm_data, amrex::Vector< 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 &time, const amrex::MultiFab &cons_in, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, int max_iters=100)
 
template<typename FluxIter >
void compute_fluxes (const int &lev, const int &max_iters, 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=1.0, const amrex::Real zscale=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)
 
template<typename FluxCalc >
void compute_SurfaceLayer_bcs (const int &lev, amrex::Vector< const amrex::MultiFab * > mfs, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Tau_lev, amrex::MultiFab *xheat_flux, amrex::MultiFab *yheat_flux, amrex::MultiFab *zheat_flux, amrex::MultiFab *xqv_flux, amrex::MultiFab *yqv_flux, amrex::MultiFab *zqv_flux, const amrex::MultiFab *z_phys, const FluxCalc &flux_comp)
 
void fill_tsurf_with_sst_and_tsk (const int &lev, const amrex::Real &time)
 
void 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 FluxCalc >
void compute_SurfaceLayer_bcs (const int &lev, Vector< const MultiFab * > mfs, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab *xheat_flux, MultiFab *yheat_flux, MultiFab *zheat_flux, MultiFab *xqv_flux, MultiFab *yqv_flux, MultiFab *zqv_flux, const MultiFab *z_phys, const FluxCalc &flux_comp)
 
template<typename PBLHeightEstimator >
void compute_pblh (const int &lev, Vector< Vector< MultiFab >> &vars, MultiFab *z_phys_cc, const PBLHeightEstimator &est, 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_time
 
amrex::Real m_stop_time
 
amrex::Real m_bdy_time_interval
 
bool m_include_wstar = false
 
amrex::Real z0_const {0.1}
 
amrex::Real default_land_surf_temp {300.}
 
amrex::Real surf_temp
 
amrex::Real surf_heating_rate {0}
 
amrex::Real surf_temp_flux {0}
 
amrex::Real default_land_surf_moist {0.}
 
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 {0.0185}
 
bool cnk_visc {false}
 
amrex::Real depth {30.0}
 
amrex::Vector< amrex::MultiFab > z_0
 
bool m_var_z0 {false}
 
bool use_moisture
 
bool m_has_lsm_tsurf = false
 
int m_lsm_tsurf_indx = -1
 
amrex::Real m_Cd = 0.0
 
amrex::Real m_Ch = 0.0
 
amrex::Real m_Cq = 0.0
 
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
 

Detailed Description

Abstraction layer for different surface layer schemes (e.g. MOST, Cd)

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

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

Member Enumeration Documentation

◆ FluxCalcType

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.

600  {
601  MOENG = 0, ///< Moeng functional form
602  DONELAN, ///< Donelan functional form
603  CUSTOM, ///< Custom constant flux functional form
604  BULK_COEFF, ///< Bulk transfer coefficient functional form
605  ROTATE ///< Terrain rotation flux functional form
606  };

◆ MoistCalcType

Enumerator
ADIABATIC 
MOISTURE_FLUX 

Qv-flux specified.

SURFACE_MOISTURE 

Surface Qv specified.

614  {
615  ADIABATIC = 0,
616  MOISTURE_FLUX, ///< Qv-flux specified
617  SURFACE_MOISTURE ///< Surface Qv specified
618  };

◆ PBLHeightCalcType

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

◆ RoughCalcType

Enumerator
CONSTANT 

Constant z0.

CHARNOCK 
MODIFIED_CHARNOCK 
DONELAN 
WAVE_COUPLED 
620  {
621  CONSTANT = 0, ///< Constant z0
622  CHARNOCK,
623  MODIFIED_CHARNOCK,
624  DONELAN,
625  WAVE_COUPLED
626  };

◆ ThetaCalcType

Enumerator
ADIABATIC 
HEAT_FLUX 

Heat-flux specified.

SURFACE_TEMPERATURE 

Surface temperature specified.

608  {
609  ADIABATIC = 0,
610  HEAT_FLUX, ///< Heat-flux specified
611  SURFACE_TEMPERATURE ///< Surface temperature specified
612  };

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,
amrex::Real  start_time,
amrex::Real  stop_time,
amrex::Real  bdy_time_interval = 0.0 
)
inlineexplicit
45  : m_geom(geom),
46  m_rotate(use_rot_surface_flux),
47  m_start_time(start_time),
48  m_stop_time(stop_time),
49  m_bdy_time_interval(bdy_time_interval),
50  m_ma(geom, (z_phys_nd[0] != nullptr), a_pp_prefix, a_mesh_type, a_terrain_type)
51  {
52  // We have a moisture model if Qv_prim is a valid pointer
53  use_moisture = (Qv_prim[0].get());
54 
55  // Get roughness
56  amrex::ParmParse pp("erf");
57  pp.query("most.z0", z0_const);
58 
59  // Specify how to compute the flux
60  if (use_rot_surface_flux) {
62  } else {
63  std::string flux_string_in;
64  std::string flux_string{"moeng"};
65  auto read_flux = pp.query("surface_layer.flux_type", flux_string_in);
66  if (read_flux) {
67  flux_string = amrex::toLower(flux_string_in);
68  }
69  if (flux_string == "donelan") {
71  } else if (flux_string == "moeng") {
73  } else if (flux_string == "bulk_coeff") {
75  } else if (flux_string == "custom") {
77  } else {
78  amrex::Abort("Undefined MOST flux type!");
79  }
80  }
81 
82  // Include w* to handle free convection (Beljaars 1995, QJRMS)
83  pp.query("most.include_wstar", m_include_wstar);
84 
85  std::string pblh_string_in;
86  std::string pblh_string{"none"};
87  auto read_pblh = pp.query("most.pblh_calc", pblh_string_in);
88  if (read_pblh) {
89  pblh_string = amrex::toLower(pblh_string_in);
90  }
91  if (pblh_string == "none") {
93  } else if (pblh_string == "mynn25") {
95  } else if (pblh_string == "mynnedmf") {
97  } else if (pblh_string == "ysu") {
99  } else if (pblh_string == "mrf") {
101  } else {
102  amrex::Abort("Undefined PBLH calc type!");
103  }
104 
105  // Get surface temperature
106  auto erf_st = pp.query("most.surf_temp", surf_temp);
107  if (erf_st) { default_land_surf_temp = surf_temp; }
108 
109  // Get surface moisture
110  bool erf_sq = false;
111  if (use_moisture) { erf_sq = pp.query("most.surf_moist", surf_moist); }
112  if (erf_sq) { default_land_surf_moist = surf_moist; }
113 
114  // Custom type user must specify the fluxes
118  pp.get("most.ustar", custom_ustar);
119  pp.get("most.tstar", custom_tstar);
120  pp.get("most.qstar", custom_qstar);
121  pp.query("most.rhosurf", custom_rhosurf);
122  if (custom_qstar != 0) {
123  AMREX_ASSERT_WITH_MESSAGE(use_moisture,
124  "Specified custom MOST qv flux without moisture model!");
125  }
126  amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
127  << custom_ustar << " " << custom_tstar << " "
128  << custom_qstar << std::endl;
129 
130  // Bulk transfer coefficient (must specify coeffs and surface values)
131  } else if (flux_type == FluxCalcType::BULK_COEFF) {
132  pp.get("most.Cd", m_Cd);
133  pp.get("most.Ch", m_Ch);
134  pp.get("most.Cq", m_Cq);
135  pp.get("most.surf_temp", default_land_surf_temp);
136  pp.get("most.surf_moist", default_land_surf_moist);
137  amrex::Print() << "Using specified Cd, Ch, Cq for MOST = "
138  << m_Cd << " " << m_Ch << " "
139  << m_Cq << std::endl;
140 
141  // Specify surface temperature/moisture or surface flux
142  } else {
143  if (erf_st) {
145  pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
146  surf_heating_rate = surf_heating_rate / 3600.0; // [K/s]
147  if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
148  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
149  }
150  } else {
151  pp.query("most.surf_temp_flux", surf_temp_flux);
152 
153  if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
154  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
155  }
156  if (std::abs(surf_temp_flux) >
159  } else {
161  }
162  }
163 
164  if (erf_sq) {
166  } else {
167  pp.query("most.surf_moist_flux", surf_moist_flux);
168  if (std::abs(surf_moist_flux) >
171  } else {
173  }
174  }
175  }
176 
177  // Make sure the inputs file doesn't try to use most.roughness_type
178  std::string bogus_input;
179  if (pp.query("most.roughness_type", bogus_input) > 0) {
180  amrex::Abort("most.roughness_type is deprecated; use "
181  "most.roughness_type_land and/or most.roughness_type_sea");
182  }
183 
184  // Specify how to compute the surface flux over land (if there is any)
185  std::string rough_land_string_in;
186  std::string rough_land_string{"constant"};
187  auto read_rough_land =
188  pp.query("most.roughness_type_land", rough_land_string_in);
189  if (read_rough_land) {
190  rough_land_string = amrex::toLower(rough_land_string_in);
191  }
192  if (rough_land_string == "constant") {
194  } else {
195  amrex::Abort("Undefined MOST roughness type for land!");
196  }
197 
198  // Specify how to compute the surface flux over sea (if there is any)
199  std::string rough_sea_string_in;
200  std::string rough_sea_string{"charnock"};
201  auto read_rough_sea = pp.query("most.roughness_type_sea", rough_sea_string_in);
202  if (read_rough_sea) {
203  rough_sea_string = amrex::toLower(rough_sea_string_in);
204  }
205  if (rough_sea_string == "charnock") {
207  pp.query("most.charnock_constant", cnk_a);
208  pp.query("most.charnock_viscosity", cnk_visc);
209  if (cnk_a > 0) {
210  amrex::Print() << "If there is water, Charnock relation with C_a="
211  << cnk_a << (cnk_visc ? " and viscosity" : "")
212  << " will be used" << std::endl;
213  } else {
214  amrex::Print() << "If there is water, Charnock relation with variable "
215  "Charnock parameter (COARE3.0)"
216  << (cnk_visc ? " and viscosity" : "") << " will be used"
217  << std::endl;
218  }
219  } else if (rough_sea_string == "coare3.0") {
221  amrex::Print() << "If there is water, Charnock relation with variable "
222  "Charnock parameter (COARE3.0)"
223  << (cnk_visc ? " and viscosity" : "") << " will be used"
224  << std::endl;
225  cnk_a = -1;
226  } else if (rough_sea_string == "donelan") {
228  } else if (rough_sea_string == "modified_charnock") {
230  pp.query("most.modified_charnock_depth", depth);
231  } else if (rough_sea_string == "wave_coupled") {
233  } else if (rough_sea_string == "constant") {
235  } else {
236  amrex::Abort("Undefined MOST roughness type for sea!");
237  }
238 
239  // use skin temperature instead of sea-surface temperature
240  // (wrfinput data may have lower resolution SST data)
241  pp.query("most.ignore_sst", m_ignore_sst);
242  } // constructor
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:233
ThetaCalcType theta_type
Definition: ERF_SurfaceLayer.H:631
bool m_include_wstar
Definition: ERF_SurfaceLayer.H:645
bool m_rotate
Definition: ERF_SurfaceLayer.H:640
PBLHeightCalcType pblh_type
Definition: ERF_SurfaceLayer.H:635
bool use_moisture
Definition: ERF_SurfaceLayer.H:665
amrex::Real m_Cq
Definition: ERF_SurfaceLayer.H:671
RoughCalcType rough_type_land
Definition: ERF_SurfaceLayer.H:633
amrex::Real z0_const
Definition: ERF_SurfaceLayer.H:646
amrex::Real cnk_a
Definition: ERF_SurfaceLayer.H:659
amrex::Real m_Ch
Definition: ERF_SurfaceLayer.H:670
amrex::Real surf_temp
Definition: ERF_SurfaceLayer.H:648
amrex::Real surf_moist_flux
Definition: ERF_SurfaceLayer.H:653
RoughCalcType rough_type_sea
Definition: ERF_SurfaceLayer.H:634
amrex::Real surf_moist
Definition: ERF_SurfaceLayer.H:652
bool m_ignore_sst
Definition: ERF_SurfaceLayer.H:673
amrex::Real m_bdy_time_interval
Definition: ERF_SurfaceLayer.H:643
amrex::Real m_stop_time
Definition: ERF_SurfaceLayer.H:642
amrex::Real custom_qstar
Definition: ERF_SurfaceLayer.H:656
amrex::Real custom_rhosurf
Definition: ERF_SurfaceLayer.H:657
amrex::Real m_start_time
Definition: ERF_SurfaceLayer.H:641
@ 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:661
amrex::Real default_land_surf_moist
Definition: ERF_SurfaceLayer.H:651
amrex::Real surf_temp_flux
Definition: ERF_SurfaceLayer.H:650
amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_SurfaceLayer.H:639
amrex::Real custom_tstar
Definition: ERF_SurfaceLayer.H:655
bool cnk_visc
Definition: ERF_SurfaceLayer.H:660
amrex::Real surf_heating_rate
Definition: ERF_SurfaceLayer.H:649
FluxCalcType flux_type
Definition: ERF_SurfaceLayer.H:630
MoistCalcType moist_type
Definition: ERF_SurfaceLayer.H:632
amrex::Real custom_ustar
Definition: ERF_SurfaceLayer.H:654
amrex::Real m_Cd
Definition: ERF_SurfaceLayer.H:669
amrex::Real default_land_surf_temp
Definition: ERF_SurfaceLayer.H:647
@ SURFACE_TEMPERATURE
Surface temperature specified.
@ HEAT_FLUX
Heat-flux specified.
MOSTAverage m_ma
Definition: ERF_SurfaceLayer.H:675
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
Here is the call graph for this function:

Member Function Documentation

◆ compute_fluxes()

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

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

Parameters
[in]levCurrent level
[in]max_itersmaximum iterations to use
[in]most_fluxstructure to iteratively compute ustar and tstar
186 {
187  // Pointers to the computed averages
188  const auto *const tm_ptr = m_ma.get_average(lev,2); // potential temperature
189  const auto *const qvm_ptr = m_ma.get_average(lev,3); // water vapor mixing ratio
190  const auto *const tvm_ptr = m_ma.get_average(lev,4); // virtual potential temperature
191  const auto *const umm_ptr = m_ma.get_average(lev,5); // horizontal velocity magnitude
192  const auto *const zref_ptr = m_ma.get_zref(lev); // reference height
193 
194  for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
195  {
196  Box gtbx = mfi.growntilebox();
197 
198  auto u_star_arr = u_star[lev]->array(mfi);
199  auto t_star_arr = t_star[lev]->array(mfi);
200  auto q_star_arr = q_star[lev]->array(mfi);
201  auto t_surf_arr = t_surf[lev]->array(mfi);
202  auto q_surf_arr = q_surf[lev]->array(mfi);
203  auto olen_arr = olen[lev]->array(mfi);
204 
205  const auto tm_arr = tm_ptr->array(mfi);
206  const auto tvm_arr = tvm_ptr->array(mfi);
207  const auto qvm_arr = qvm_ptr->array(mfi);
208  const auto umm_arr = umm_ptr->array(mfi);
209  const auto zref_arr = zref_ptr->array(mfi);
210  const auto z0_arr = z_0[lev].array(mfi);
211 
212  // PBL height if we need to calculate wstar for the Beljaars correction
213  // TODO: can/should we apply this in LES mode?
214  const auto w_star_arr = (m_include_wstar) ? w_star[lev].get()->array(mfi) : Array4<Real> {};
215  const auto pblh_arr = (m_include_wstar) ? pblh[lev].get()->array(mfi) : Array4<Real> {};
216 
217  // Wave properties if they exist
218  const auto Hwave_arr = (m_Hwave_lev[lev]) ? m_Hwave_lev[lev]->array(mfi) : Array4<Real> {};
219  const auto Lwave_arr = (m_Lwave_lev[lev]) ? m_Lwave_lev[lev]->array(mfi) : Array4<Real> {};
220  const auto eta_arr = (m_eddyDiffs_lev[lev]) ? m_eddyDiffs_lev[lev]->array(mfi) : Array4<Real> {};
221 
222  // Land mask array if it exists
223  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
224  Array4<int> {};
225 
226  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
227  {
228  if (( is_land && lmask_arr(i,j,k) == 1) ||
229  (!is_land && lmask_arr(i,j,k) == 0))
230  {
231  most_flux.iterate_flux(i, j, k, max_iters,
232  zref_arr, // set in most average
233  z0_arr, // updated if(!is_land)
234  umm_arr, tm_arr, tvm_arr, qvm_arr,
235  u_star_arr, // updated
236  w_star_arr, // updated if(m_include_wstar)
237  t_star_arr, q_star_arr, // updated
238  t_surf_arr, q_surf_arr, olen_arr, // updated
239  pblh_arr, // updated if(m_include_wstar)
240  Hwave_arr, Lwave_arr, eta_arr);
241  }
242  });
243  }
244 }
amrex::MultiFab * get_zref(const int &lev) const
Definition: ERF_MOSTAverage.H:105
const amrex::MultiFab * get_average(const int &lev, const int &comp) const
Definition: ERF_MOSTAverage.H:102
amrex::Vector< amrex::Vector< amrex::iMultiFab * > > m_lmask_lev
Definition: ERF_SurfaceLayer.H:687
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_surf
Definition: ERF_SurfaceLayer.H:682
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_star
Definition: ERF_SurfaceLayer.H:679
amrex::Vector< amrex::MultiFab * > m_Lwave_lev
Definition: ERF_SurfaceLayer.H:693
amrex::Vector< amrex::MultiFab > z_0
Definition: ERF_SurfaceLayer.H:662
amrex::Vector< std::unique_ptr< amrex::MultiFab > > w_star
Definition: ERF_SurfaceLayer.H:677
amrex::Vector< std::unique_ptr< amrex::MultiFab > > u_star
Definition: ERF_SurfaceLayer.H:676
amrex::Vector< std::unique_ptr< amrex::MultiFab > > q_surf
Definition: ERF_SurfaceLayer.H:683
amrex::Vector< amrex::MultiFab * > m_Hwave_lev
Definition: ERF_SurfaceLayer.H:692
amrex::Vector< std::unique_ptr< amrex::MultiFab > > t_star
Definition: ERF_SurfaceLayer.H:678
amrex::Vector< amrex::MultiFab * > m_eddyDiffs_lev
Definition: ERF_SurfaceLayer.H:694
amrex::Vector< std::unique_ptr< amrex::MultiFab > > olen
Definition: ERF_SurfaceLayer.H:680
amrex::Vector< std::unique_ptr< amrex::MultiFab > > pblh
Definition: ERF_SurfaceLayer.H:681

◆ 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 
)
709 {
710  est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
711  vars[lev][Vars::cons],m_lmask_lev[lev][0],
712  moisture_indices);
713 }
@ cons
Definition: ERF_IndexDefines.H:140

◆ 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
323 {
324  bool rotate = m_rotate;
325  const int klo = m_geom[lev].Domain().smallEnd(2);
326  const auto& dxInv = m_geom[lev].InvCellSizeArray();
327  for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
328  {
329  // Get field arrays
330  const auto cons_arr = mfs[Vars::cons]->array(mfi);
331  const auto velx_arr = mfs[Vars::xvel]->array(mfi);
332  const auto vely_arr = mfs[Vars::yvel]->array(mfi);
333  const auto velz_arr = mfs[Vars::zvel]->array(mfi);
334 
335  // Diffusive stress vars
336  auto t13_arr = Tau_lev[TauType::tau13]->array(mfi);
337  auto t31_arr = (Tau_lev[TauType::tau31]) ? Tau_lev[TauType::tau31]->array(mfi) : Array4<Real>{};
338 
339  auto t23_arr = Tau_lev[TauType::tau23]->array(mfi);
340  auto t32_arr = (Tau_lev[TauType::tau32]) ? Tau_lev[TauType::tau32]->array(mfi) : Array4<Real>{};
341 
342  auto hfx3_arr = zheat_flux->array(mfi);
343  auto qfx3_arr = (zqv_flux) ? zqv_flux->array(mfi) : Array4<Real>{};
344 
345  // Rotated stress vars
346  auto t11_arr = (m_rotate) ? Tau_lev[TauType::tau11]->array(mfi) : Array4<Real>{};
347  auto t22_arr = (m_rotate) ? Tau_lev[TauType::tau22]->array(mfi) : Array4<Real>{};
348  auto t33_arr = (m_rotate) ? Tau_lev[TauType::tau33]->array(mfi) : Array4<Real>{};
349  auto t12_arr = (m_rotate) ? Tau_lev[TauType::tau12]->array(mfi) : Array4<Real>{};
350  auto t21_arr = (m_rotate) ? Tau_lev[TauType::tau21]->array(mfi) : Array4<Real>{};
351 
352  auto hfx1_arr = (m_rotate) ? xheat_flux->array(mfi) : Array4<Real>{};
353  auto hfx2_arr = (m_rotate) ? yheat_flux->array(mfi) : Array4<Real>{};
354  auto qfx1_arr = (m_rotate && xqv_flux) ? xqv_flux->array(mfi) : Array4<Real>{};
355  auto qfx2_arr = (m_rotate && yqv_flux) ? yqv_flux->array(mfi) : Array4<Real>{};
356 
357  // Terrain
358  const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};
359 
360  // Get average arrays
361  const auto *const u_mean = m_ma.get_average(lev,0);
362  const auto *const v_mean = m_ma.get_average(lev,1);
363  const auto *const t_mean = m_ma.get_average(lev,2);
364  const auto *const q_mean = m_ma.get_average(lev,3);
365  const auto *const u_mag_mean = m_ma.get_average(lev,5);
366 
367  const auto um_arr = u_mean->array(mfi);
368  const auto vm_arr = v_mean->array(mfi);
369  const auto tm_arr = t_mean->array(mfi);
370  const auto qm_arr = q_mean->array(mfi);
371  const auto umm_arr = u_mag_mean->array(mfi);
372 
373  // Get derived arrays
374  const auto u_star_arr = u_star[lev]->array(mfi);
375  const auto t_star_arr = t_star[lev]->array(mfi);
376  const auto q_star_arr = q_star[lev]->array(mfi);
377  const auto t_surf_arr = t_surf[lev]->array(mfi);
378  const auto q_surf_arr = q_surf[lev]->array(mfi);
379 
380  // Get LSM fluxes
381  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
382  Array4<int> {};
383  auto lsm_t_flux_arr = Array4<Real> {};
384  auto lsm_q_flux_arr = Array4<Real> {};
385  auto lsm_tau13_arr = Array4<Real> {};
386  auto lsm_tau23_arr = Array4<Real> {};
387  for (int n(0); n<m_lsm_flux_lev[lev].size(); ++n) {
388  if (toLower(m_lsm_flux_name[n]) == "t_flux") { lsm_t_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
389  if (toLower(m_lsm_flux_name[n]) == "q_flux") { lsm_q_flux_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
390  if (toLower(m_lsm_flux_name[n]) == "tau13") { lsm_tau13_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
391  if (toLower(m_lsm_flux_name[n]) == "tau23") { lsm_tau23_arr = m_lsm_flux_lev[lev][n]->array(mfi); }
392  }
393 
394 
395  // Rho*Theta flux
396  //============================================================================
397  Box bx = mfi.tilebox();
398  if (bx.smallEnd(2) != klo) { continue; }
399  bx.makeSlab(2,klo);
400  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
401  {
402  // Valid theta flux from LSM and over land
403  Real Tflux;
404  int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
405  if (lsm_t_flux_arr && is_land) {
406  Tflux = lsm_t_flux_arr(i,j,k);
407  } else {
408  Tflux = flux_comp.compute_t_flux(i, j, k,
409  cons_arr, velx_arr, vely_arr,
410  umm_arr, tm_arr, u_star_arr,
411  t_star_arr, t_surf_arr);
412 
413  }
414 
415  // Do scalar flux rotations?
416  if (rotate) {
417  rotate_scalar_flux(i, j, k, Tflux, dxInv, zphys_arr,
418  hfx1_arr, hfx2_arr, hfx3_arr);
419  } else {
420  hfx3_arr(i,j,klo) = Tflux;
421  }
422  });
423 
424  // Rho*Qv flux
425  //============================================================================
426  if (use_moisture) {
427  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
428  {
429  // Valid qv flux from LSM and over land
430  Real Qflux;
431  int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
432  if (lsm_q_flux_arr && is_land) {
433  Qflux = lsm_q_flux_arr(i,j,k);
434  } else {
435  Qflux = flux_comp.compute_q_flux(i, j, k,
436  cons_arr, velx_arr, vely_arr,
437  umm_arr, qm_arr, u_star_arr,
438  q_star_arr, q_surf_arr);
439  }
440 
441  // Do scalar flux rotations?
442  if (rotate) {
443  rotate_scalar_flux(i, j, k, Qflux, dxInv, zphys_arr,
444  qfx1_arr, qfx2_arr, qfx3_arr);
445  } else {
446  qfx3_arr(i,j,k) = Qflux;
447  }
448  });
449  } // custom
450 
451  if (!rotate) {
452  // Rho*u flux
453  //============================================================================
454  Box bxx = surroundingNodes(bx,0);
455  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
456  {
457  // Valid tau13 from LSM and over land
458  Real stressx;
459  int is_land_hi = (lmask_arr) ? lmask_arr(i ,j,klo) : 1;
460  int is_land_lo = (lmask_arr) ? lmask_arr(i-1,j,klo) : 1;
461  if (lsm_tau13_arr && (is_land_hi || is_land_lo)) {
462  stressx = 0.;
463  if (!is_land_hi || !is_land_lo) {
464  stressx += 0.5 * flux_comp.compute_u_flux(i, j, k,
465  cons_arr, velx_arr, vely_arr,
466  umm_arr, um_arr, u_star_arr);
467  }
468  if (is_land_hi) {
469  stressx += 0.5 * lsm_tau13_arr(i ,j,k);
470  }
471  if (is_land_lo) {
472  stressx += 0.5 * lsm_tau13_arr(i-1,j,k);
473  }
474  } else {
475  stressx = flux_comp.compute_u_flux(i, j, k,
476  cons_arr, velx_arr, vely_arr,
477  umm_arr, um_arr, u_star_arr);
478  }
479 
480  t13_arr(i,j,k) = stressx;
481  if (t31_arr) { t31_arr(i,j,k) = stressx; }
482  });
483 
484  // Rho*v flux
485  //============================================================================
486  Box bxy = surroundingNodes(bx,1);
487  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
488  {
489  // Valid tau13 from LSM and over land
490  Real stressy;
491  int is_land_hi = (lmask_arr) ? lmask_arr(i,j ,klo) : 1;
492  int is_land_lo = (lmask_arr) ? lmask_arr(i,j-1,klo) : 1;
493  if (lsm_tau23_arr && (is_land_hi || is_land_lo)) {
494  stressy = 0.;
495  if (!is_land_hi || !is_land_lo) {
496  stressy += 0.5 * flux_comp.compute_v_flux(i, j, k,
497  cons_arr, velx_arr, vely_arr,
498  umm_arr, vm_arr, u_star_arr);
499  }
500  if (is_land_hi) {
501  stressy += 0.5 * lsm_tau23_arr(i,j ,k);
502  }
503  if (is_land_lo) {
504  stressy += 0.5 * lsm_tau23_arr(i,j-1,k);
505  }
506  } else {
507  stressy = flux_comp.compute_v_flux(i, j, k,
508  cons_arr, velx_arr, vely_arr,
509  umm_arr, vm_arr, u_star_arr);
510  }
511 
512  t23_arr(i,j,k) = stressy;
513  if (t32_arr) { t32_arr(i,j,k) = stressy; }
514  });
515  } else {
516  // All fluxes with rotation
517  //============================================================================
518  Box bxxy = convert(bx, IntVect(1,1,0));
519  ParallelFor(bxxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
520  {
521  Real stresst = flux_comp.compute_u_flux(i, j, k,
522  cons_arr, velx_arr, vely_arr,
523  umm_arr, um_arr, u_star_arr);
524  rotate_stress_tensor(i, j, k, stresst, dxInv, zphys_arr,
525  velx_arr, vely_arr, velz_arr,
526  t11_arr, t22_arr, t33_arr,
527  t12_arr, t21_arr,
528  t13_arr, t31_arr,
529  t23_arr, t32_arr);
530  });
531  }
532  } // mfiter
533 }
@ 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::Real Real
Definition: ERF_ShocInterface.H:19
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:609
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:630
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_flux_lev
Definition: ERF_SurfaceLayer.H:689
amrex::Vector< std::string > m_lsm_flux_name
Definition: ERF_SurfaceLayer.H:691
@ xvel
Definition: ERF_IndexDefines.H:141
@ zvel
Definition: ERF_IndexDefines.H:143
@ yvel
Definition: ERF_IndexDefines.H:142
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 
)
620 {
621  // NOTE: We have already tested a moisture model exists
622 
623  // Populate q_surf with qsat over water
624  auto dz = m_geom[lev].CellSize(2);
625  for (MFIter mfi(*q_surf[lev]); mfi.isValid(); ++mfi)
626  {
627  Box gtbx = mfi.growntilebox();
628 
629  auto t_surf_arr = t_surf[lev]->array(mfi);
630  auto q_surf_arr = q_surf[lev]->array(mfi);
631  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
632  Array4<int> {};
633  const auto cons_arr = cons_in.const_array(mfi);
634  const auto z_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) :
635  Array4<const Real> {};
636 
637  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
638  {
639  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
640  if (!is_land) {
641  auto deltaZ = (z_arr) ? Compute_Zrel_AtCellCenter(i,j,k,z_arr) :
642  0.5*dz;
643  auto Rho = cons_arr(i,j,k,Rho_comp);
644  auto RTh = cons_arr(i,j,k,RhoTheta_comp);
645  auto Qv = cons_arr(i,j,k,RhoQ1_comp) / Rho;
646  auto P_cc = getPgivenRTh(RTh, Qv);
647  P_cc += Rho*CONST_GRAV*deltaZ;
648  P_cc *= 0.01;
649  erf_qsatw(t_surf_arr(i,j,k), P_cc, q_surf_arr(i,j,k));
650  }
651  });
652  }
653  q_surf[lev]->FillBoundary(m_geom[lev].periodicity());
654 }
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:166
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:387
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 
)
538 {
539  int n_times_in_sst = m_sst_lev[lev].size();
540 
541  int n_time_lo, n_time_hi;
542  Real alpha;
543 
544  if (n_times_in_sst > 1) {
545  // Time interpolation
547  int n_time = static_cast<int>( elapsed_time / dT);
548  n_time_lo = n_time;
549  n_time_hi = n_time+1;
550  alpha = (elapsed_time - n_time * dT) / dT;
551  if ((elapsed_time == m_stop_time-m_start_time) && (alpha==0)) {
552  // stop time coincides with final lowinp slice -- don't try to
553  // interpolate from the following time slice
554  n_time -= 1;
555  n_time_lo -= 1;
556  n_time_hi -= 1;
557  alpha = 1.0;
558  }
559  AMREX_ALWAYS_ASSERT( (n_time >= 0) && (n_time < (m_sst_lev[lev].size()-1)));
560  } else {
561  n_time_lo = 0;
562  n_time_hi = 0;
563  alpha = 1.0;
564  }
565  AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
566 
567  Real oma = 1.0 - alpha;
568 
569  // Define a default land surface temperature if we don't read in tsk
571 
572  bool use_tsk = (m_tsk_lev[lev][0]);
573  bool ignore_sst = m_ignore_sst;
574 
575  // Populate t_surf
576  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
577  {
578  Box gtbx = mfi.growntilebox();
579 
580  auto t_surf_arr = t_surf[lev]->array(mfi);
581 
582  const auto sst_lo_arr = m_sst_lev[lev][n_time_lo]->const_array(mfi);
583  const auto sst_hi_arr = m_sst_lev[lev][n_time_hi]->const_array(mfi);
584 
585  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
586  Array4<int> {};
587 
588  if (use_tsk) {
589  const auto tsk_arr = m_tsk_lev[lev][n_time_lo]->const_array(mfi);
590  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
591  {
592  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
593  if (!is_land && !ignore_sst) {
594  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
595  + alpha * sst_hi_arr(i,j,k);
596  } else {
597  t_surf_arr(i,j,k) = tsk_arr(i,j,k);
598  }
599  });
600  } else {
601  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
602  {
603  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
604  if (!is_land) {
605  t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
606  + alpha * sst_hi_arr(i,j,k);
607  } else {
608  t_surf_arr(i,j,k) = lst;
609  }
610  });
611  }
612  }
613  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
614 }
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_sst_lev
Definition: ERF_SurfaceLayer.H:685
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_tsk_lev
Definition: ERF_SurfaceLayer.H:686

◆ get_lmask()

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

◆ get_lsm_tsurf()

void SurfaceLayer::get_lsm_tsurf ( const int &  lev)
658 {
659  for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
660  {
661  Box gtbx = mfi.growntilebox();
662 
663  // NOTE: LSM does not carry lateral ghost cells.
664  // This copies the valid box into the ghost cells.
665  // Fillboundary is called after this to pick up the
666  // interior ghost and periodic directions.
667  Box vbx = mfi.validbox();
668  int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
669  int j_lo = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1);
670 
671  auto t_surf_arr = t_surf[lev]->array(mfi);
672  auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) :
673  Array4<int> {};
674  const auto lsm_arr = m_lsm_data_lev[lev][m_lsm_tsurf_indx]->const_array(mfi);
675 
676  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
677  {
678  int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1;
679  if (is_land) {
680  int li = amrex::min(amrex::max(i, i_lo), i_hi);
681  int lj = amrex::min(amrex::max(j, j_lo), j_hi);
682  t_surf_arr(i,j,k) = lsm_arr(li,lj,k);
683  }
684  });
685  }
686 }
int m_lsm_tsurf_indx
Definition: ERF_SurfaceLayer.H:667
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_lsm_data_lev
Definition: ERF_SurfaceLayer.H:688

◆ get_mac_avg()

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

◆ get_olen()

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

◆ get_pblh()

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

◆ get_q_star()

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

◆ get_q_surf()

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

◆ get_t_star()

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

◆ get_t_surf()

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

◆ get_u_star()

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

◆ get_w_star()

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

◆ get_z0()

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

◆ get_zref()

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

◆ 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
265 {
267  moeng_flux flux_comp;
268  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
269  xheat_flux, yheat_flux, zheat_flux,
270  xqv_flux, yqv_flux, zqv_flux,
271  z_phys, flux_comp);
272  } else if (flux_type == FluxCalcType::DONELAN) {
273  donelan_flux flux_comp;
274  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
275  xheat_flux, yheat_flux, zheat_flux,
276  xqv_flux, yqv_flux, zqv_flux,
277  z_phys, flux_comp);
278  } else if (flux_type == FluxCalcType::ROTATE) {
279  rotate_flux flux_comp;
280  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
281  xheat_flux, yheat_flux, zheat_flux,
282  xqv_flux, yqv_flux, zqv_flux,
283  z_phys, flux_comp);
284  } else if (flux_type == FluxCalcType::BULK_COEFF) {
285  bulk_coeff_flux flux_comp(m_Cd, m_Ch, m_Cq);
286  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
287  xheat_flux, yheat_flux, zheat_flux,
288  xqv_flux, yqv_flux, zqv_flux,
289  z_phys, flux_comp);
290  } else if (flux_type == FluxCalcType::CUSTOM) {
291  custom_flux flux_comp(specified_rho_surf);
292  compute_SurfaceLayer_bcs(lev, mfs, Tau_lev,
293  xheat_flux, yheat_flux, zheat_flux,
294  xqv_flux, yqv_flux, zqv_flux,
295  z_phys, flux_comp);
296  } else {
297  amrex::Abort("Unknown surface layer flux calculation type");
298  }
299 }
bool specified_rho_surf
Definition: ERF_SurfaceLayer.H:658
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:2180
Definition: ERF_MOSTStress.H:2063
Definition: ERF_MOSTStress.H:1931
Definition: ERF_MOSTStress.H:1768
Definition: ERF_MOSTStress.H:2307

◆ 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 = 1.0,
const amrex::Real  zscale = 700.0 
)
721 {
722  Print() << "Initializing TKE from surface layer ustar on level " << lev << std::endl;
723 
724  constexpr Real small = 0.01;
725 
726  for (MFIter mfi(cons); mfi.isValid(); ++mfi)
727  {
728  Box gtbx = mfi.tilebox();
729 
730  auto const& u_star_arr = u_star[lev]->const_array(mfi);
731  auto const& z_phys_arr = z_phys_nd->const_array(mfi);
732 
733  auto const& cons_arr = cons.array(mfi);
734 
735  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
736  {
737  Real rho = cons_arr(i, j, k, Rho_comp);
738  Real ust = u_star_arr(i, j, 0);
739  Real tke0 = tkefac * ust * ust; // surface value
740  Real zagl = Compute_Zrel_AtCellCenter(i, j, k, z_phys_arr);
741 
742  // linearly tapering profile -- following WRF, approximate top of
743  // PBL as ustar * zscale
744  cons_arr(i, j, k, RhoKE_comp) = rho * tke0 * std::max(
745  (ust * zscale - zagl) / (std::max(ust, small) * zscale),
746  small);
747  });
748  }
749 }
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
@ rho
Definition: ERF_Kessler.H:22
Here is the call graph for this function:

◆ lmask_min_reduce()

int SurfaceLayer::lmask_min_reduce ( amrex::iMultiFab &  lmask,
const int &  nghost 
)
inline
574  {
575  int lmask_min = amrex::ReduceMin(lmask, nghost, [=] AMREX_GPU_HOST_DEVICE(
576  amrex::Box const& bx, amrex::Array4<int const> const& lm_arr) -> int
577  {
578  int locmin = std::numeric_limits<int>::max();
579  const auto lo = lbound(bx);
580  const auto hi = ubound(bx);
581  for (int j = lo.y; j <= hi.y; ++j) {
582  for (int i = lo.x; i <= hi.x; ++i) {
583  locmin = std::min(locmin, lm_arr(i, j, 0));
584  }
585  }
586  return locmin;
587  });
588 
589  return lmask_min;
590  }

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

◆ read_custom_roughness()

void SurfaceLayer::read_custom_roughness ( const int &  lev,
const std::string &  fname 
)
755 {
756  // Read the file if we have it
757  if (!fname.empty()) {
758  // Only the ioproc reads the file
759  Gpu::HostVector<Real> m_x,m_y,m_z0;
760  if (ParallelDescriptor::IOProcessor()) {
761  Print()<<"Reading MOST roughness file at level " << lev << " : " << fname << std::endl;
762  std::ifstream file(fname);
763  Real value1,value2,value3;
764  while(file>>value1>>value2>>value3){
765  m_x.push_back(value1);
766  m_y.push_back(value2);
767  m_z0.push_back(value3);
768  }
769  file.close();
770 
771  AMREX_ALWAYS_ASSERT(m_x.size() == m_y.size());
772  AMREX_ALWAYS_ASSERT(m_x.size() == m_z0.size());
773  }
774 
775  // Broadcast the whole domain to every rank
776  int ioproc = ParallelDescriptor::IOProcessorNumber();
777  int nnode = m_x.size();
778  ParallelDescriptor::Bcast(&nnode, 1, ioproc);
779 
780  if (!ParallelDescriptor::IOProcessor()) {
781  m_x.resize(nnode);
782  m_y.resize(nnode);
783  m_z0.resize(nnode);
784  }
785  ParallelDescriptor::Bcast(m_x.data() , nnode, ioproc);
786  ParallelDescriptor::Bcast(m_y.data() , nnode, ioproc);
787  ParallelDescriptor::Bcast(m_z0.data(), nnode, ioproc);
788 
789  // Copy data to the GPU
790  Gpu::DeviceVector<Real> d_x(nnode),d_y(nnode),d_z0(nnode);
791  Gpu::copy(Gpu::hostToDevice, m_x.begin(), m_x.end(), d_x.begin());
792  Gpu::copy(Gpu::hostToDevice, m_y.begin(), m_y.end(), d_y.begin());
793  Gpu::copy(Gpu::hostToDevice, m_z0.begin(), m_z0.end(), d_z0.begin());
794  Real* xp = d_x.data();
795  Real* yp = d_y.data();
796  Real* z0p = d_z0.data();
797 
798  // Each rank populates it's z_0[lev] MultiFab
799  for (MFIter mfi(z_0[lev]); mfi.isValid(); ++mfi)
800  {
801  Box gtbx = mfi.growntilebox();
802 
803  // Populate z_phys data
804  Real tol = 1.0e-4;
805  auto dx = m_geom[lev].CellSizeArray();
806  auto ProbLoArr = m_geom[lev].ProbLoArray();
807  int ilo = m_geom[lev].Domain().smallEnd(0);
808  int jlo = m_geom[lev].Domain().smallEnd(1);
809  int klo = 0;
810  int ihi = m_geom[lev].Domain().bigEnd(0);
811  int jhi = m_geom[lev].Domain().bigEnd(1);
812 
813  Array4<Real> const& z0_arr = z_0[lev].array(mfi);
814  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
815  {
816  // Clip indices for ghost-cells
817  int ii = amrex::min(amrex::max(i,ilo),ihi);
818  int jj = amrex::min(amrex::max(j,jlo),jhi);
819 
820  // Location of nodes
821  Real x = ProbLoArr[0] + ii * dx[0];
822  Real y = ProbLoArr[1] + jj * dx[1];
823  int inode = ii + jj * (ihi-ilo+2); // stride is Nx+1
824  if (std::sqrt(std::pow(x-xp[inode],2)+std::pow(y-yp[inode],2)) < tol) {
825  z0_arr(i,j,klo) = z0p[inode];
826  } else {
827  // Unexpected list order, do brute force search
828  Real z0loc = 0.0;
829  bool found = false;
830  for (int n=0; n<nnode; ++n) {
831  Real delta=std::sqrt(std::pow(x-xp[n],2)+std::pow(y-yp[n],2));
832  if (delta < tol) {
833  found = true;
834  z0loc = z0p[n];
835  break;
836  }
837  }
838  AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
839  amrex::ignore_unused(found);
840  z0_arr(i,j,klo) = z0loc;
841  }
842  });
843  } // mfi
844  } else {
845  AMREX_ALWAYS_ASSERT(lev > 0);
846 
847  Print()<<"Interpolating MOST roughness at level " << lev << std::endl;
848 
849  // Create a BC mapper that uses FOEXTRAP at domain bndry
850  Vector<int> bc_lo(3,ERFBCType::foextrap);
851  Vector<int> bc_hi(3,ERFBCType::foextrap);
852  Vector<BCRec> bcr; bcr.push_back(BCRec(bc_lo.data(),bc_hi.data()));
853 
854  // Create ref ratio
855  IntVect ratio;
856  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
857  ratio[idim] = m_geom[lev].Domain().length(idim) / m_geom[0].Domain().length(idim);
858  }
859 
860  // Create interp object and interpolate from the coarsest grid
861  MFInterpolater* interp = &mf_cell_cons_interp;
862  interp->interp(z_0[0] , 0,
863  z_0[lev], 0,
864  1, z_0[lev].nGrowVect(),
865  m_geom[0], m_geom[lev],
866  m_geom[lev].Domain(),ratio,
867  bcr, 0);
868  }
869 }
@ m_y
Definition: ERF_DataStruct.H:25
@ m_x
Definition: ERF_DataStruct.H:24
@ foextrap
Definition: ERF_IndexDefines.H:208

Referenced by make_SurfaceLayer_at_level().

Here is the caller graph for this function:

◆ set_q_surf()

void SurfaceLayer::set_q_surf ( const int &  lev,
const amrex::Real  qsurf 
)
inline
562 { q_surf[lev]->setVal(qsurf); }

◆ set_t_surf()

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

◆ update_fluxes()

void SurfaceLayer::update_fluxes ( const int &  lev,
const amrex::Real time,
const amrex::MultiFab &  cons_in,
const std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
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
17 {
18  // Update with SST/TSK data if we have a valid pointer
19  if (m_sst_lev[lev][0]) {
20  fill_tsurf_with_sst_and_tsk(lev, time);
21  }
22 
23  // Apply heating rate if needed
25  update_surf_temp(time);
26  }
27 
28  // Update qsurf with qsat over sea
29  if (use_moisture) {
30  fill_qsurf_with_qsat(lev, cons_in, z_phys_nd);
31  }
32 
33  // Update land surface temp if we have a valid pointer
34  if (m_has_lsm_tsurf) { get_lsm_tsurf(lev); }
35 
36  // Fill interior ghost cells
37  t_surf[lev]->FillBoundary(m_geom[lev].periodicity());
38 
39  // Compute plane averages for all vars (regardless of flux type)
41 
42 
43  // ***************************************************************
44  // Iterate the fluxes if moeng type
45  // First iterate over land -- the only model for surface roughness
46  // over land is RoughCalcType::CONSTANT
47  // ***************************************************************
50  bool is_land = true;
51  // Do we have a constant flux for moisture over land?
52  bool cons_qflux = ( (moist_type == MoistCalcType::MOISTURE_FLUX) ||
56  surface_flux most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
57  compute_fluxes(lev, max_iters, most_flux, is_land);
58  } else {
59  amrex::Abort("Unknown value for rough_type_land");
60  }
63  surface_temp most_flux(surf_temp_flux, surf_moist_flux, cons_qflux);
64  compute_fluxes(lev, max_iters, most_flux, is_land);
65  } else {
66  amrex::Abort("Unknown value for rough_type_land");
67  }
68  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
72  compute_fluxes(lev, max_iters, most_flux, is_land);
73  } else {
74  amrex::Abort("Unknown value for rough_type_land");
75  }
76  } else {
77  amrex::Abort("Unknown value for theta_type");
78  }
79  } // MOENG -- LAND
80 
81  // ***************************************************************
82  // Iterate the fluxes if moeng type
83  // Next iterate over sea -- the models for surface roughness
84  // over sea are CHARNOCK, DONELAN, MODIFIED_CHARNOCK or WAVE_COUPLED
85  // ***************************************************************
88  bool is_land = false;
89  // NOTE: Do not allow default to adiabatic over sea (we have Qvs at surface)
90  // Do we have a constant flux for moisture over sea?
91  bool cons_qflux = (moist_type == MoistCalcType::MOISTURE_FLUX);
95  cnk_a, cnk_visc, cons_qflux);
96  compute_fluxes(lev, max_iters, most_flux, is_land);
99  depth, cons_qflux);
100  compute_fluxes(lev, max_iters, most_flux, is_land);
101  } else if (rough_type_sea == RoughCalcType::DONELAN) {
103  cons_qflux);
104  compute_fluxes(lev, max_iters, most_flux, is_land);
107  cons_qflux);
108  compute_fluxes(lev, max_iters, most_flux, is_land);
109  } else {
110  amrex::Abort("Unknown value for rough_type_sea");
111  }
112 
116  cnk_a, cnk_visc, cons_qflux);
117  compute_fluxes(lev, max_iters, most_flux, is_land);
120  depth, cons_qflux);
121  compute_fluxes(lev, max_iters, most_flux, is_land);
122  } else if (rough_type_sea == RoughCalcType::DONELAN) {
124  cons_qflux);
125  compute_fluxes(lev, max_iters, most_flux, is_land);
128  cons_qflux);
129  compute_fluxes(lev, max_iters, most_flux, is_land);
130  } else {
131  amrex::Abort("Unknown value for rough_type_sea");
132  }
133 
134  } else if ((theta_type == ThetaCalcType::ADIABATIC) &&
138  cnk_a, cnk_visc);
139  compute_fluxes(lev, max_iters, most_flux, is_land);
142  depth);
143  compute_fluxes(lev, max_iters, most_flux, is_land);
144  } else if (rough_type_sea == RoughCalcType::DONELAN) {
146  compute_fluxes(lev, max_iters, most_flux, is_land);
149  compute_fluxes(lev, max_iters, most_flux, is_land);
150  } else {
151  amrex::Abort("Unknown value for rough_type_sea");
152  }
153  } else {
154  amrex::Abort("Unknown value for theta_type");
155  }
156 
157  } // MOENG -- SEA
158 
160  if (custom_rhosurf > 0) {
161  specified_rho_surf = true;
162  u_star[lev]->setVal(std::sqrt(custom_rhosurf) * custom_ustar);
163  t_star[lev]->setVal(custom_rhosurf * custom_tstar);
164  q_star[lev]->setVal(custom_rhosurf * custom_qstar);
165  } else {
166  u_star[lev]->setVal(custom_ustar);
167  t_star[lev]->setVal(custom_tstar);
168  q_star[lev]->setVal(custom_qstar);
169  }
170  }
171 }
void compute_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:732
void update_surf_temp(const amrex::Real &time)
Definition: ERF_SurfaceLayer.H:520
void get_lsm_tsurf(const int &lev)
Definition: ERF_SurfaceLayer.cpp:657
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:617
void compute_fluxes(const int &lev, const int &max_iters, const FluxIter &most_flux, bool is_land)
Definition: ERF_SurfaceLayer.cpp:182
void fill_tsurf_with_sst_and_tsk(const int &lev, const amrex::Real &time)
Definition: ERF_SurfaceLayer.cpp:536
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

◆ 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
537  {
538  m_ma.update_field_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
539  }
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:265
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 
)
693 {
695  MYNNPBLH estimator;
696  compute_pblh(lev, vars, z_phys_cc, estimator, moisture_indices);
698  amrex::Error("YSU/MRF PBLH calc not implemented yet");
699  }
700 }
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
592  {
593  m_sst_lev[lev][itime] = sst_ptr;
594  }

◆ update_surf_temp()

void SurfaceLayer::update_surf_temp ( const amrex::Real time)
inline
521  {
522  if (surf_heating_rate != 0) {
523  int nlevs = static_cast<int>(m_geom.size());
524  for (int lev = 0; lev < nlevs; lev++) {
525  t_surf[lev]->setVal(surf_temp + surf_heating_rate * time);
526  amrex::Print() << "Surface temp at t=" << time << ": "
527  << surf_temp + surf_heating_rate * time << std::endl;
528  }
529  }
530  }

◆ update_tsk_ptr()

void SurfaceLayer::update_tsk_ptr ( const int  lev,
const int  itime,
amrex::MultiFab *  tsk_ptr 
)
inline
596  {
597  m_tsk_lev[lev][itime] = tsk_ptr;
598  }

Member Data Documentation

◆ cnk_a

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

Referenced by SurfaceLayer().

◆ cnk_visc

bool SurfaceLayer::cnk_visc {false}
private

Referenced by SurfaceLayer().

◆ custom_qstar

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

Referenced by SurfaceLayer().

◆ custom_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 {0.}
private

◆ default_land_surf_temp

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

◆ depth

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

Referenced by SurfaceLayer().

◆ flux_type

FluxCalcType SurfaceLayer::flux_type {FluxCalcType::MOENG}

◆ m_bdy_time_interval

amrex::Real SurfaceLayer::m_bdy_time_interval
private

◆ m_Cd

amrex::Real SurfaceLayer::m_Cd = 0.0
private

Referenced by SurfaceLayer().

◆ m_Ch

amrex::Real SurfaceLayer::m_Ch = 0.0
private

Referenced by SurfaceLayer().

◆ m_Cq

amrex::Real SurfaceLayer::m_Cq = 0.0
private

Referenced by SurfaceLayer().

◆ m_eddyDiffs_lev

amrex::Vector<amrex::MultiFab*> SurfaceLayer::m_eddyDiffs_lev
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_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_time

amrex::Real SurfaceLayer::m_start_time
private

◆ m_stop_time

amrex::Real SurfaceLayer::m_stop_time
private

◆ m_tsk_lev

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

◆ 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

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

◆ z_0

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

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