ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MakeMomSources.cpp File Reference
#include <AMReX_MultiFab.H>
#include <AMReX_ArrayLim.H>
#include <AMReX_BCRec.H>
#include <AMReX_TableData.H>
#include <AMReX_GpuContainers.H>
#include "ERF_ImmersedForcing.H"
#include "ERF_NumericalDiffusion.H"
#include "ERF_PlaneAverage.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_SrcHeaders.H"
#include "ERF_Utils.H"
Include dependency graph for ERF_MakeMomSources.cpp:

Functions

void make_mom_sources (double time_d, double dt, const Vector< MultiFab > &S_data, const MultiFab *z_phys_nd, const MultiFab *z_phys_cc, Vector< Real > &stretched_dz_h, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &wvel, MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, const MultiFab &base_state, MultiFab *forest_drag, MultiFab *terrain_blank, MultiFab *terrain_blank_xface, MultiFab *terrain_blank_yface, MultiFab *terrain_blank_zface, MultiFab *cosPhi_mf, MultiFab *sinPhi_mf, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &, const Real *dptr_u_geos, const Real *dptr_v_geos, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const amrex::Real *d_sinesq_at_lev, const amrex::Real *d_sinesq_stag_at_lev, const Vector< Real * > d_sponge_ptrs_at_lev, const Vector< MultiFab > *forecast_state_at_lev, const InputSoundingData &input_sounding_data, const LargeScaleForcingData &lsf_data, std::unique_ptr< amrex::MultiFab > &lsf_tendencies, const eb_ &ebfact, bool is_slow_step)
 

Function Documentation

◆ make_mom_sources()

void make_mom_sources ( double  time_d,
double  dt,
const Vector< MultiFab > &  S_data,
const MultiFab *  z_phys_nd,
const MultiFab *  z_phys_cc,
Vector< Real > &  stretched_dz_h,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  wvel,
MultiFab &  xmom_src,
MultiFab &  ymom_src,
MultiFab &  zmom_src,
const MultiFab &  base_state,
MultiFab *  forest_drag,
MultiFab *  terrain_blank,
MultiFab *  terrain_blank_xface,
MultiFab *  terrain_blank_yface,
MultiFab *  terrain_blank_zface,
MultiFab *  cosPhi_mf,
MultiFab *  sinPhi_mf,
const Geometry  geom,
const SolverChoice &  solverChoice,
Vector< std::unique_ptr< MultiFab >> &  ,
const Real *  dptr_u_geos,
const Real *  dptr_v_geos,
const Real *  dptr_wbar_sub,
const Vector< Real * >  d_rayleigh_ptrs_at_lev,
const amrex::Real *  d_sinesq_at_lev,
const amrex::Real *  d_sinesq_stag_at_lev,
const Vector< Real * >  d_sponge_ptrs_at_lev,
const Vector< MultiFab > *  forecast_state_at_lev,
const InputSoundingData &  input_sounding_data,
const LargeScaleForcingData &  lsf_data,
std::unique_ptr< amrex::MultiFab > &  lsf_tendencies,
const eb_ &  ebfact,
bool  is_slow_step 
)

Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum.

Parameters
[in]timecurrent time
[in]dtcurrent slow or fast timestep size
[in]S_datacurrent solution
[in]xvelx-component of velocity
[in]yvely-component of velocity
[in]xmom_srcsource terms for x-momentum
[in]ymom_srcsource terms for y-momentum
[in]zmom_srcsource terms for z-momentum
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]mapfacmap factors
[in]dptr_u_geoscustom geostrophic wind profile
[in]dptr_v_geoscustom geostrophic wind profile
[in]dptr_wbar_subsubsidence source term
[in]d_rayleigh_ptrs_at_levVector of {strength of Rayleigh damping, reference value for xvel/yvel/zvel/theta} used to define Rayleigh damping
[in]d_sinesq_at_levsin( (pi/2) (z-z_t)/(damping depth)) at cell centers
[in]d_sinesq_stag_at_levsin( (pi/2) (z-z_t)/(damping depth)) at z-faces
74 {
75  BL_PROFILE_REGION("erf_make_mom_sources()");
76 
77  Real time = static_cast<Real>(time_d);
78 
79  Box domain(geom.Domain());
80  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
81 
82  // Initialize sources to zero each time we may use them
83  xmom_src.setVal(0.0);
84  ymom_src.setVal(0.0);
85  zmom_src.setVal(0.0);
86 
87  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
88 
89  // flags to apply certain source terms in substep call only
90  bool use_Rayleigh_fast_uv = ( (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit) ||
91  (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastImplicit) );
92  bool use_Rayleigh_fast_w = (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit);
93  bool use_canopy_fast = solverChoice.forest_substep;
94  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
95 
96  // *****************************************************************************
97  // Define source term for all three components of momenta from
98  // 1. Coriolis forcing for (xmom,ymom,zmom)
99  // 2. Rayleigh damping for (xmom,ymom,zmom)
100  // 3. Constant / height-dependent geostrophic forcing
101  // 4. Subsidence
102  // 5. Nudging towards input sounding data
103  // 6. Numerical diffusion for (xmom,ymom,zmom)
104  // 7. Sponge
105  // 8. Forest canopy
106  // 9a. Immersed forcing for terrain
107  // 9b. Immersed forcing for buildings
108  // 10. Constant mass flux
109  // *****************************************************************************
110  // NOTE: buoyancy is now computed in a separate routine - it should not appear here
111  // *****************************************************************************
112  //const bool l_use_ndiff = solverChoice.use_num_diff;
113 
114  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
115  if (solverChoice.do_forest_drag) {
116  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
117  }
118  }
119 
120 
121  // *****************************************************************************
122  // Data for Coriolis forcing
123  // *****************************************************************************
124  auto use_coriolis = solverChoice.use_coriolis;
125  auto coriolis_factor = solverChoice.coriolis_factor;
126  auto cosphi = solverChoice.cosphi;
127  auto sinphi = solverChoice.sinphi;
128  auto var_coriolis = solverChoice.variable_coriolis;
129 
130  // *****************************************************************************
131  // Flag for Geostrophic forcing
132  // *****************************************************************************
133  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
134  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
135 
136  // *****************************************************************************
137  // Data for Rayleigh damping
138  // *****************************************************************************
139  auto rayleigh_damp_U = solverChoice.dampingChoice.rayleigh_damp_U;
140  auto rayleigh_damp_V = solverChoice.dampingChoice.rayleigh_damp_V;
141  auto rayleigh_damp_W = solverChoice.dampingChoice.rayleigh_damp_W;
142 
143  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
144  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
145  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
146 
147  // *****************************************************************************
148  // Data for constant mass flux
149  // *****************************************************************************
150  bool enforce_massflux_x = (solverChoice.const_massflux_u != 0);
151  bool enforce_massflux_y = (solverChoice.const_massflux_v != 0);
152  Real U_target = solverChoice.const_massflux_u;
153  Real V_target = solverChoice.const_massflux_v;
154  int massflux_klo = solverChoice.massflux_klo;
155  int massflux_khi = solverChoice.massflux_khi;
156 
157  // These will be updated by integrating through the planar average profiles
158  Real rhoUA_target{0};
159  Real rhoVA_target{0};
160  Real rhoUA{0};
161  Real rhoVA{0};
162 
163  // *****************************************************************************
164  // Planar averages for subsidence, nudging, or constant mass flux
165  // *****************************************************************************
166  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
167  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
168 
169  if (is_slow_step && (dptr_wbar_sub ||
170  (solverChoice.nudging_from_input_sounding &&
171  (solverChoice.large_scale_forcing || solverChoice.nudging_u)) ||
172  enforce_massflux_x || enforce_massflux_y))
173  {
174  // The plane averaging operates at fixed z rather than fixed physical
175  // height. It is valid for a variable-dz fitted mesh only when that
176  // mesh has been explicitly validated as horizontally flat.
177  const bool planar_averages_valid =
178  solverChoice.mesh_type != MeshType::VariableDz ||
179  solverChoice.flat_terrain;
180  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
181  planar_averages_valid,
182  "Fixed-index planar averages require a non-variable-dz mesh or "
183  "a validated flat fitted mesh (erf.flat_terrain = true).");
184 
185  const int offset = 1;
186  const int u_offset = 1;
187  const int v_offset = 1;
188 
189  //
190  // We use the alias here to control ncomp inside the PlaneAverage
191  //
192  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
193 
194  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = offset;
195  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
196  r_ave.compute_averages(ZDir(), r_ave.field());
197 
198  int ncell = r_ave.ncell_line();
199  Gpu::HostVector< Real> r_plane_h(ncell);
200  Gpu::DeviceVector< Real> r_plane_d(ncell);
201 
202  r_ave.line_average(Rho_comp, r_plane_h);
203 
204  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
205 
206  Real* dptr_r = r_plane_d.data();
207 
208  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
209  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
210 
211  dptr_r_plane = r_plane_tab.table();
212  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
213  {
214  dptr_r_plane(k-offset) = dptr_r[k];
215  });
216 
217  // U and V momentum
218  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
219  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
220 
221  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
222  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
223 
224  u_ave.compute_averages(ZDir(), u_ave.field());
225  v_ave.compute_averages(ZDir(), v_ave.field());
226 
227  int u_ncell = u_ave.ncell_line();
228  int v_ncell = v_ave.ncell_line();
229  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
230  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
231 
232  u_ave.line_average(0, u_plane_h);
233  v_ave.line_average(0, v_plane_h);
234 
235  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
236  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
237 
238  Real* dptr_u = u_plane_d.data();
239  Real* dptr_v = v_plane_d.data();
240 
241  Box udomain = domain; udomain.grow(2,ng_u[2]);
242  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
243  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
244  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
245 
246  dptr_u_plane = u_plane_tab.table();
247  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
248  {
249  dptr_u_plane(k-u_offset) = dptr_u[k];
250  });
251 
252  dptr_v_plane = v_plane_tab.table();
253  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
254  {
255  dptr_v_plane(k-v_offset) = dptr_v[k];
256  });
257 
258  // sum in z for massflux adjustment
259  if (enforce_massflux_x || enforce_massflux_y) {
260  Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
261  Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
262 
263  if (solverChoice.mesh_type == MeshType::ConstantDz) {
264  // note: massflux_khi corresponds to unstaggered indices in this case
265  rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
266  u_plane_h.begin() + u_offset + massflux_khi+1, zero);
267  rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
268  v_plane_h.begin() + v_offset + massflux_khi+1, zero);
269  rhoUA_target = std::accumulate(r_plane_h.begin() + offset + massflux_klo,
270  r_plane_h.begin() + offset + massflux_khi+1, zero);
271  rhoVA_target = rhoUA_target;
272 
273  rhoUA *= geom.CellSize(2) * Ly;
274  rhoVA *= geom.CellSize(2) * Lx;
275  rhoUA_target *= geom.CellSize(2) * Ly;
276  rhoVA_target *= geom.CellSize(2) * Lx;
277 
278  } else if (solverChoice.mesh_type == MeshType::StretchedDz) {
279  // note: massflux_khi corresponds to staggered indices in this case
280  for (int k=massflux_klo; k < massflux_khi; ++k) {
281  rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
282  rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
283  rhoUA_target += r_plane_h[k + offset] * stretched_dz_h[k];
284  }
285  rhoVA_target = rhoUA_target;
286 
287  rhoUA *= Ly;
288  rhoVA *= Lx;
289  rhoUA_target *= Ly;
290  rhoVA_target *= Lx;
291  }
292 
293  // at this point, this is integrated rho*dA
294  rhoUA_target *= U_target;
295  rhoVA_target *= V_target;
296 
297  Print() << "Integrated mass flux : " << rhoUA << " " << rhoVA
298  << " (target: " << rhoUA_target << " " << rhoVA_target << ")"
299  << std::endl;
300  }
301  }
302 
303  // *****************************************************************************
304  // Add all the other forcings
305  // *****************************************************************************
306  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
307  {
308  Box tbx = mfi.nodaltilebox(0);
309  Box tby = mfi.nodaltilebox(1);
310  Box tbz = mfi.nodaltilebox(2);
311  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
312 
313  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
314  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
315  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
316  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
317 
318  const Array4<const Real>& u = xvel.array(mfi);
319  const Array4<const Real>& v = yvel.array(mfi);
320  const Array4<const Real>& w = wvel.array(mfi);
321 
322  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
323  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
324  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
325 
326  const Array4<const Real>& r0 = r_hse.const_array(mfi);
327 
328  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
329  Array4<const Real>{};
330  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
331  Array4<const Real>{};
332  const Array4<const Real>& t_blank_xface_arr = (terrain_blank_xface) ? terrain_blank_xface->const_array(mfi) :
333  Array4<const Real>{};
334  const Array4<const Real>& t_blank_yface_arr = (terrain_blank_yface) ? terrain_blank_yface->const_array(mfi) :
335  Array4<const Real>{};
336  const Array4<const Real>& t_blank_zface_arr = (terrain_blank_zface) ? terrain_blank_zface->const_array(mfi) :
337  Array4<const Real>{};
338 
339  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
340  Array4<const Real>{};
341  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
342  Array4<const Real>{};
343 
344  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
345  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
346 
347 
348  // *****************************************************************************
349  // 1. Add Coriolis forcing (this assumes east is +x, north is +y)
350  // *****************************************************************************
351  if (use_coriolis && is_slow_step) {
352  if(solverChoice.init_type == InitType::HindCast) {
353  const Array4<const Real>& latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
354  ParallelFor(tbx, tby, tbz,
355  [=] AMREX_GPU_DEVICE (int i, int j, int k)
356  {
357  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
358  Real v_loc = fourth * (v(i,j+1,k) + v(i,j,k) + v(i-1,j+1,k) + v(i-1,j,k));
359  Real w_loc = fourth * (w(i,j,k+1) + w(i,j,k) + w(i-1,j,k+1) + w(i-1,j,k));
360  Real latitude = latlon_arr(i,j,k,0);
361  Real sphi_loc = std::sin(latitude*PI/Real(180.0));
362  Real cphi_loc = std::cos(latitude*PI/Real(180.0));
363  xmom_src_arr(i, j, k) += coriolis_factor * rho_on_u_face * (v_loc * sphi_loc - w_loc * cphi_loc);
364  },
365  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
366  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
367  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j-1,k) + u(i,j-1,k));
368  Real latitude = latlon_arr(i,j,k,0);
369  Real sphi_loc = std::sin(latitude*PI/Real(180.0));
370  ymom_src_arr(i, j, k) += -coriolis_factor * rho_on_v_face * u_loc * sphi_loc;
371  },
372  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
373  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
374  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j,k-1) + u(i,j,k-1));
375  Real latitude = latlon_arr(i,j,k,0);
376  Real cphi_loc = std::cos(latitude*PI/Real(180.0));
377  zmom_src_arr(i, j, k) += coriolis_factor * rho_on_w_face * u_loc * cphi_loc;
378  });
379  }
380  else if (var_coriolis && (sinPhi_mf) && (cosPhi_mf)) {
381  ParallelFor(tbx, tby, tbz,
382  [=] AMREX_GPU_DEVICE (int i, int j, int k)
383  {
384  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
385  Real v_loc = fourth * (v(i,j+1,k) + v(i,j,k) + v(i-1,j+1,k) + v(i-1,j,k));
386  Real w_loc = fourth * (w(i,j,k+1) + w(i,j,k) + w(i-1,j,k+1) + w(i-1,j,k));
387  Real sphi_loc = myhalf * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
388  Real cphi_loc = myhalf * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
389  xmom_src_arr(i, j, k) += coriolis_factor * rho_on_u_face * (v_loc * sphi_loc - w_loc * cphi_loc);
390  },
391  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
392  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
393  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j-1,k) + u(i,j-1,k));
394  Real sphi_loc = myhalf * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
395  ymom_src_arr(i, j, k) += -coriolis_factor * rho_on_v_face * u_loc * sphi_loc;
396  },
397  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
398  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
399  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j,k-1) + u(i,j,k-1));
400  zmom_src_arr(i, j, k) += coriolis_factor * rho_on_w_face * u_loc * cphi_arr(i,j,0);
401  });
402  } else {
403  if (solverChoice.terrain_type == TerrainType::EB) {
404  Array4<const Real> u_volfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
405  Array4<const Real> v_volfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
406  Array4<const Real> w_volfrac = (ebfact.get_w_const_factory())->getVolFrac().const_array(mfi);
407  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
408  Real v_loc = 0.0;
409  Real w_loc = 0.0;
410  Real v_vol = v_volfrac(i,j+1,k) + v_volfrac(i,j,k) + v_volfrac(i-1,j+1,k) + v_volfrac(i-1,j,k);
411  Real w_vol = w_volfrac(i,j,k+1) + w_volfrac(i,j,k) + w_volfrac(i-1,j,k+1) + w_volfrac(i-1,j,k);
412  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
413  if (v_vol > 0.0) {
414  v_loc = ( v_volfrac(i ,j+1,k) * v(i ,j+1,k) + v_volfrac(i ,j,k) * v(i ,j,k)
415  + v_volfrac(i-1,j+1,k) * v(i-1,j+1,k) + v_volfrac(i-1,j,k) * v(i-1,j,k)) / v_vol;
416  }
417  if (w_vol > 0.0) {
418  w_loc = ( w_volfrac(i ,j,k+1) * w(i ,j,k+1) + w_volfrac(i ,j,k) * w(i ,j,k)
419  + w_volfrac(i-1,j,k+1) * w(i-1,j,k+1) + w_volfrac(i-1,j,k) * w(i-1,j,k)) / w_vol;
420  }
421  xmom_src_arr(i, j, k) += coriolis_factor * rho_on_u_face * (v_loc * sinphi - w_loc * cosphi);
422  });
423  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
424  Real u_loc = 0.0;
425  Real u_vol = u_volfrac(i+1,j,k) + u_volfrac(i,j,k) + u_volfrac(i+1,j-1,k) + u_volfrac(i,j-1,k);
426  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
427  if (u_vol > 0.0) {
428  u_loc = ( u_volfrac(i+1,j ,k) * u(i+1,j ,k) + u_volfrac(i,j ,k) * u(i,j ,k)
429  + u_volfrac(i+1,j-1,k) * u(i+1,j-1,k) + u_volfrac(i,j-1,k) * u(i,j-1,k)) / u_vol;
430  }
431  ymom_src_arr(i, j, k) += -coriolis_factor * rho_on_v_face * u_loc * sinphi;
432  });
433  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
434  Real u_loc = 0.0;
435  Real u_vol = u_volfrac(i+1,j,k) + u_volfrac(i,j,k) + u_volfrac(i+1,j,k-1) + u_volfrac(i,j,k-1);
436  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
437  if (u_vol > 0.0) {
438  u_loc = ( u_volfrac(i+1,j,k ) * u(i+1,j,k ) + u_volfrac(i,j,k) * u(i,j,k )
439  + u_volfrac(i+1,j,k-1) * u(i+1,j,k-1) + u_volfrac(i,j,k-1) * u(i,j,k-1)) / u_vol;
440  }
441  zmom_src_arr(i, j, k) += coriolis_factor * rho_on_w_face * u_loc * cosphi;
442  });
443  } else {
444  ParallelFor(tbx, tby, tbz,
445  [=] AMREX_GPU_DEVICE (int i, int j, int k)
446  {
447  Real v_loc = fourth * (v(i,j+1,k) + v(i,j,k) + v(i-1,j+1,k) + v(i-1,j,k));
448  Real w_loc = fourth * (w(i,j,k+1) + w(i,j,k) + w(i-1,j,k+1) + w(i-1,j,k));
449  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
450  xmom_src_arr(i, j, k) += rho_on_u_face * ( coriolis_factor * (v_loc * sinphi - w_loc * cosphi) );
451  },
452  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
453  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
454  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j-1,k) + u(i,j-1,k));
455  ymom_src_arr(i, j, k) += rho_on_v_face * ( -coriolis_factor * u_loc * sinphi );
456  },
457  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
458  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j,k-1) + u(i,j,k-1));
459  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
460  zmom_src_arr(i, j, k) += rho_on_w_face * ( coriolis_factor * u_loc * cosphi );
461  });
462  }
463  } // var_coriolis
464  } // use_coriolis
465 
466  // *****************************************************************************
467  // 2. Add RAYLEIGH damping
468  // *****************************************************************************
469  Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef;
470 
471  if ( (is_slow_step && !use_Rayleigh_fast_uv) || (!is_slow_step && use_Rayleigh_fast_uv)) {
472  if (rayleigh_damp_U) {
473  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
474  {
475  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
476  Real uu = rho_u(i,j,k) / rho_on_u_face;
477  Real sinesq = d_sinesq_at_lev[k];
478  xmom_src_arr(i, j, k) -= dampcoef*sinesq * (uu - ubar[k]) * rho_on_u_face;
479  });
480  }
481 
482  if (rayleigh_damp_V) {
483  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
484  {
485  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
486  Real vv = rho_v(i,j,k) / rho_on_v_face;
487  Real sinesq = d_sinesq_at_lev[k];
488  ymom_src_arr(i, j, k) -= dampcoef*sinesq * (vv - vbar[k]) * rho_on_v_face;
489  });
490  }
491  } // fast or slow step
492 
493  if ( (is_slow_step && !use_Rayleigh_fast_w) || (!is_slow_step && use_Rayleigh_fast_w)) {
494  if (rayleigh_damp_W) {
495  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
496  {
497  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
498  Real ww = rho_w(i,j,k) / rho_on_w_face;
499  Real sinesq = d_sinesq_stag_at_lev[k];
500  zmom_src_arr(i, j, k) -= dampcoef*sinesq * (ww - wbar[k]) * rho_on_w_face;
501  });
502  }
503  } // fast or slow step
504 
505  // *****************************************************************************
506  // 3a. Add constant GEOSTROPHIC forcing
507  // *****************************************************************************
508  if (is_slow_step) {
509  ParallelFor(tbx, tby, tbz,
510  [=] AMREX_GPU_DEVICE (int i, int j, int k)
511  {
512  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
513  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
514  },
515  [=] AMREX_GPU_DEVICE (int i, int j, int k)
516  {
517  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
518  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
519  },
520  [=] AMREX_GPU_DEVICE (int i, int j, int k)
521  {
522  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
523  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
524  });
525  }
526 
527  // *****************************************************************************
528  // 3b. Add height-dependent GEOSTROPHIC forcing
529  // *****************************************************************************
530  if (geo_wind_profile && is_slow_step) {
531  ParallelFor(tbx, tby,
532  [=] AMREX_GPU_DEVICE (int i, int j, int k)
533  {
534  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
535  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
536  },
537  [=] AMREX_GPU_DEVICE (int i, int j, int k)
538  {
539  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
540  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
541  });
542  } // geo_wind_profile
543 
544  // *****************************************************************************
545  // 4. Add custom SUBSIDENCE terms
546  // *****************************************************************************
547  if (solverChoice.custom_w_subsidence && is_slow_step && solverChoice.do_mom_advection) {
548  if (solverChoice.custom_forcing_prim_vars) {
549  const int nr = Rho_comp;
550  ParallelFor(tbx, tby,
551  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
552  {
553  Real dzInv = myhalf*dxInv[2];
554  if (z_nd_arr) {
555  Real z_xf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
556  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
557  Real z_xf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
558  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
559  dzInv = one / (z_xf_hi - z_xf_lo);
560  }
561  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
562  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
563  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
564  Real wbar_xf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
565  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
566  },
567  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
568  {
569  Real dzInv = myhalf*dxInv[2];
570  if (z_nd_arr) {
571  Real z_yf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
572  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
573  Real z_yf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
574  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
575  dzInv = one / (z_yf_hi - z_yf_lo);
576  }
577  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
578  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
579  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
580  Real wbar_yf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
581  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
582  });
583  } else {
584  ParallelFor(tbx, tby,
585  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
586  {
587  Real dzInv = myhalf*dxInv[2];
588  if (z_nd_arr) {
589  Real z_xf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
590  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
591  Real z_xf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
592  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
593  dzInv = one / (z_xf_hi - z_xf_lo);
594  }
595  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
596  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
597  Real wbar_xf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
598  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
599  },
600  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
601  {
602  Real dzInv = myhalf*dxInv[2];
603  if (z_nd_arr) {
604  Real z_yf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
605  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
606  Real z_yf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
607  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
608  dzInv = one / (z_yf_hi - z_yf_lo);
609  }
610  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
611  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
612  Real wbar_yf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
613  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
614  });
615  }
616  }
617 
618  if (solverChoice.large_scale_forcing && is_slow_step)
619  {
620  // subsidence terms for U and V
621  auto lsf_arr = lsf_tendencies->const_array(mfi);
622 
623  const int kmin = domain.smallEnd(2) + 1; // minimum k for vertical subsidence
624  const int kmax = domain.bigEnd(2) - 1; // maximum k for vertical subsidence
625 
626  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
627  {
628  if (k >= kmin && k <= kmax) {
629  int k1, k2;
630  if (lsf_arr(i, j, k, 2) >= 0.0)
631  {
632  k1 = k;
633  k2 = k-1;
634  } else {
635  k1 = k+1;
636  k2 = k;
637  }
638 
639  // one-sided difference over the single cell k1->k2
640  Real dzInv = dxInv[2];
641  if (z_nd_arr) {
642  Real z_uf_1 = fourth * ( z_nd_arr(i,j,k1 ) + z_nd_arr(i,j+1,k1 )
643  + z_nd_arr(i,j,k1+1) + z_nd_arr(i,j+1,k1+1) );
644  Real z_uf_2 = fourth * ( z_nd_arr(i,j,k2 ) + z_nd_arr(i,j+1,k2 )
645  + z_nd_arr(i,j,k2+1) + z_nd_arr(i,j+1,k2+1) );
646  dzInv = one / (z_uf_1 - z_uf_2);
647  }
648  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
649  amrex::Real utend = -(dzInv * lsf_arr(i, j, k, 2)) * ( u(i, j, k1) - u(i, j, k2) );
650  xmom_src_arr(i, j, k) += utend * rho_on_u_face;
651  }
652  });
653 
654  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
655  {
656  if (k >= kmin && k <= kmax) {
657  int k1, k2;
658  if (lsf_arr(i, j, k, 2) >= 0.0)
659  {
660  k1 = k;
661  k2 = k-1;
662  } else {
663  k1 = k+1;
664  k2 = k;
665  }
666 
667  // one-sided difference over the single cell k1->k2
668  Real dzInv = dxInv[2];
669  if (z_nd_arr) {
670  Real z_vf_1 = fourth * ( z_nd_arr(i,j,k1 ) + z_nd_arr(i+1,j,k1 )
671  + z_nd_arr(i,j,k1+1) + z_nd_arr(i+1,j,k1+1) );
672  Real z_vf_2 = fourth * ( z_nd_arr(i,j,k2 ) + z_nd_arr(i+1,j,k2 )
673  + z_nd_arr(i,j,k2+1) + z_nd_arr(i+1,j,k2+1) );
674  dzInv = one / (z_vf_1 - z_vf_2);
675  }
676  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
677  amrex::Real vtend = -(dzInv * lsf_arr(i, j, k, 2)) * ( v(i, j, k1) - v(i, j, k2) );
678  ymom_src_arr(i, j, k) += vtend * rho_on_v_face;
679  }
680  });
681  }
682 
683  // *************************************************************************************
684  // 5. Add nudging towards value specified in input sounding
685  // NOTE: nudging_u only gates sounding-based wind nudging (the !large_scale_forcing
686  // branch below); LSF-based wind nudging (towards uls/vls) is controlled solely by
687  // large_scale_forcing, per the NOTE in Inputs.rst. Likewise, nudging_u_z1/z2 only
688  // restrict the sounding-based branch.
689  // *************************************************************************************
690  if (solverChoice.nudging_from_input_sounding && is_slow_step &&
691  (solverChoice.large_scale_forcing || solverChoice.nudging_u))
692  {
693  const bool l_lsf = solverChoice.large_scale_forcing;
694  const Real u_z1 = solverChoice.nudging_u_z1;
695  const Real u_z2 = solverChoice.nudging_u_z2;
696  Real uv_coeff_n = 1.0;
697  Real uv_coeff_np1 = 0.0;
698  Real tau = Real(1.0) / input_sounding_data.tau_nudging;
699  const Real *u_nudge_n, *u_nudge_np1, *v_nudge_n, *v_nudge_np1;
700  if (!l_lsf)
701  {
702  int itime_n = 0;
703  int itime_np1 = 0;
704  int n_sounding_times = input_sounding_data.input_sounding_time.size();
705 
706  for (int nt = 1; nt < n_sounding_times; nt++) {
707  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
708  }
709  if (itime_n == n_sounding_times-1) {
710  itime_np1 = itime_n;
711  } else {
712  itime_np1 = itime_n+1;
713  uv_coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
714  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
715  uv_coeff_n = Real(1.0) - uv_coeff_np1;
716  }
717  u_nudge_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr() + 1;
718  u_nudge_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr() + 1;
719  v_nudge_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr() + 1;
720  v_nudge_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr() + 1;
721  } else {
722  int itime_curr = 0;
723  int itime_next = 0;
724  uv_coeff_n = 1.0;
725  uv_coeff_np1 = 0.0;
726  tau = 1.0 / lsf_data.tau_lsf; // only applies to u,v LSF nudging
727 
728  lsf_data.get_forcing_time_coeffs(time, itime_curr, itime_next, uv_coeff_n, uv_coeff_np1);
729  u_nudge_n = lsf_data.u_int_lsf_d[itime_curr].dataPtr();
730  u_nudge_np1 = lsf_data.u_int_lsf_d[itime_next].dataPtr();
731  v_nudge_n = lsf_data.v_int_lsf_d[itime_curr].dataPtr();
732  v_nudge_np1 = lsf_data.v_int_lsf_d[itime_next].dataPtr();
733  }
734 
735  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
736  {
737  Real z = z_cc_arr(i,j,k);
738  if (l_lsf || (z >= u_z1 && z <= u_z2)) {
739  Real unudge = -((dptr_u_plane(k)/dptr_r_plane(k)) - (uv_coeff_n*u_nudge_n[k] + uv_coeff_np1*u_nudge_np1[k]));
740  xmom_src_arr(i, j, k) += tau * unudge * dptr_r_plane(k);
741  }
742  });
743 
744  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
745  {
746  Real z = z_cc_arr(i,j,k);
747  if (l_lsf || (z >= u_z1 && z <= u_z2)) {
748  Real vnudge = -((dptr_v_plane(k)/dptr_r_plane(k)) - (uv_coeff_n*v_nudge_n[k] + uv_coeff_np1*v_nudge_np1[k]));
749  ymom_src_arr(i, j, k) += tau * vnudge * dptr_r_plane(k);
750  }
751  });
752  }
753 
754  // *****************************************************************************
755  // 6. Add NUMERICAL DIFFUSION terms
756  // *****************************************************************************
757 #if 0
758  if (l_use_ndiff) {
759  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
760  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
761  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
762  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
763  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
764  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
765  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
766  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
767  }
768 #endif
769 
770  // *****************************************************************************
771  // 7. Add SPONGING
772  // *****************************************************************************
773  if (is_slow_step) {
774  if (solverChoice.spongeChoice.sponge_type == SpongeType::Input_Sponge)
775  {
776  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
777  z_cc_arr, xmom_src_arr, ymom_src_arr,
778  rho_u, rho_v, d_sponge_ptrs_at_lev);
779  }
780  else if (solverChoice.spongeChoice.sponge_type == SpongeType::Standard)
781  {
782  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
783  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
784  r0, z_nd_arr, z_cc_arr);
785  }
786 
787  if(solverChoice.init_type == InitType::HindCast and solverChoice.hindcast_lateral_forcing){
788 
789  const Array4<const Real>& rho_u_forecast_state = (*forecast_state_at_lev)[IntVars::xmom].array(mfi);
790  const Array4<const Real>& rho_v_forecast_state = (*forecast_state_at_lev)[IntVars::ymom].array(mfi);
791  const Array4<const Real>& rho_w_forecast_state = (*forecast_state_at_lev)[IntVars::zmom].array(mfi);
792  const Array4<const Real>& cons_forecast_state = (*forecast_state_at_lev)[IntVars::cons].array(mfi);
793  ApplyBndryForcing_Forecast(solverChoice, geom, tbx, tby, tbz, z_nd_arr,
794  xmom_src_arr, ymom_src_arr, zmom_src_arr,
795  rho_u, rho_v, rho_w,
796  rho_u_forecast_state, rho_v_forecast_state, rho_w_forecast_state,
797  cons_forecast_state);
798  }
799  }
800 
801  // *****************************************************************************
802  // 8. Add CANOPY source terms
803  // *****************************************************************************
804  if (solverChoice.do_forest_drag &&
805  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
806  ParallelFor(tbx, tby, tbz,
807  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
808  {
809  const Real ux = u(i, j, k);
810  const Real uy = fourth * ( v(i, j , k ) + v(i-1, j , k )
811  + v(i, j+1, k ) + v(i-1, j+1, k ) );
812  const Real uz = fourth * ( w(i, j , k ) + w(i-1, j , k )
813  + w(i, j , k+1) + w(i-1, j , k+1) );
814  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
815  const Real f_drag = myhalf * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
816  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
817  },
818  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
819  {
820  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
821  + u(i+1, j , k ) + u(i+1, j-1, k ) );
822  const Real uy = v(i, j, k);
823  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
824  + w(i , j , k+1) + w(i , j-1, k+1) );
825  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
826  const Real f_drag = myhalf * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
827  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
828  },
829  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
830  {
831  const amrex::Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
832  + u(i , j , k-1) + u(i+1, j , k-1) );
833  const amrex::Real uy = fourth * ( v(i , j , k ) + v(i , j+1, k )
834  + v(i , j , k-1) + v(i , j+1, k-1) );
835  const amrex::Real uz = w(i, j, k);
836  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
837  const Real f_drag = myhalf * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
838  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
839  });
840  }
841  // *****************************************************************************
842  // 9a. Add immersed source terms for terrain
843  // *****************************************************************************
844  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
845  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
846 
847  ImmersedForcingTerrain_Xmom(tbx, u, v, w, cell_data, t_blank_arr, t_blank_xface_arr,
848  z_cc_arr, xmom_src_arr, geom, solverChoice, dt);
849  ImmersedForcingTerrain_Ymom(tby, u, v, w, cell_data, t_blank_arr, t_blank_yface_arr,
850  z_cc_arr, ymom_src_arr, geom, solverChoice, dt);
851  ImmersedForcingTerrain_Zmom(tbz, u, v, w, cell_data, t_blank_arr, t_blank_zface_arr,
852  z_cc_arr, zmom_src_arr, geom, solverChoice, dt);
853  }
854 
855  // *****************************************************************************
856  // 9b. Add immersed source terms for buildings
857  // *****************************************************************************
858  // geometric properties
859  const Real* dx_arr = geom.CellSize();
860  const Real delta_xy = std::sqrt(dx_arr[0] * dx_arr[1]);
861  if ((solverChoice.buildings_type == BuildingsType::ImmersedForcing) &&
862  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)) &&
863  (delta_xy <= 50.0)) { // only apply immersed forcing when grid spacing is less than 50m
864 
865  ImmersedForcingBuildings_Xmom(tbx, u, v, w, cell_data, t_blank_arr, t_blank_xface_arr,
866  z_cc_arr, xmom_src_arr, geom, solverChoice, dt);
867  ImmersedForcingBuildings_Ymom(tby, u, v, w, cell_data, t_blank_arr, t_blank_yface_arr,
868  z_cc_arr, ymom_src_arr, geom, solverChoice, dt);
869  ImmersedForcingBuildings_Zmom(tbz, u, v, w, cell_data, t_blank_arr, t_blank_zface_arr,
870  z_cc_arr, zmom_src_arr, geom, solverChoice, dt);
871  }
872 
873  // *****************************************************************************
874  // 10. Enforce constant mass flux
875  // *****************************************************************************
876  if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
877  Real tau_inv = one / solverChoice.const_massflux_tau;
878 
879  ParallelFor(tbx, tby,
880  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
881  xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
882  },
883  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
884  ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
885  });
886  }
887 
888  } // mfi
889 }
ApplyBndryForcing_Forecast
void ApplyBndryForcing_Forecast(const SolverChoice &solverChoice, const Geometry geom, const Box &tbx, const Box &tby, const Box &tbz, const Array4< const Real > &z_phys_nd, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &rho_w, const Array4< const Real > &rho_u_initial_state, const Array4< const Real > &rho_v_initial_state, const Array4< const Real > &rho_w_initial_state, const Array4< const Real > &cons_initial_state)
Definition: ERF_ApplyBndryForcing_Forecast.cpp:8
ApplySpongeZoneBCsForMom
void ApplySpongeZoneBCsForMom(const SpongeChoice &spongeChoice, const Geometry geom, const Box &tbx, const Box &tby, const Box &tbz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &rho_w, const Array4< const Real > &r0, const Array4< const Real > &z_phys_nd, const Array4< const Real > &z_phys_cc)
Apply sponge zone damping to momentum variables.
Definition: ERF_ApplySpongeZoneBCs.cpp:211
ApplySpongeZoneBCsForMom_ReadFromFile
void ApplySpongeZoneBCsForMom_ReadFromFile(const SpongeChoice &spongeChoice, const Geometry geom, const Box &tbx, const Box &tby, const Array4< const Real > &cell_data, const Array4< const Real > &z_phys_cc, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Vector< Real * > d_sponge_ptrs_at_lev)
Definition: ERF_ApplySpongeZoneBCs_ReadFromFile.cpp:23
RayleighDampingType::FastExplicit
@ FastExplicit
RayleighDampingType::FastImplicit
@ FastImplicit
ubar
@ ubar
Definition: ERF_DataStruct.H:194
wbar
@ wbar
Definition: ERF_DataStruct.H:194
vbar
@ vbar
Definition: ERF_DataStruct.H:194
Coord::z
@ z
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:55
ImmersedForcingBuildings_Ymom
void ImmersedForcingBuildings_Ymom(const Box &tby, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_yface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
Definition: ERF_ImmersedForcing.cpp:562
ImmersedForcingTerrain_Ymom
void ImmersedForcingTerrain_Ymom(const Box &tby, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_yface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
Definition: ERF_ImmersedForcing.cpp:194
ImmersedForcingBuildings_Xmom
void ImmersedForcingBuildings_Xmom(const Box &tbx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_xface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
Definition: ERF_ImmersedForcing.cpp:363
ImmersedForcingTerrain_Xmom
void ImmersedForcingTerrain_Xmom(const Box &tbx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_xface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
Definition: ERF_ImmersedForcing.cpp:79
ImmersedForcingBuildings_Zmom
void ImmersedForcingBuildings_Zmom(const Box &tbz, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_zface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &zmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
Definition: ERF_ImmersedForcing.cpp:760
ImmersedForcingTerrain_Zmom
void ImmersedForcingTerrain_Zmom(const Box &tbz, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_zface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &zmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
Definition: ERF_ImmersedForcing.cpp:308
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:42
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
PI
constexpr amrex::Real PI
Definition: ERF_NumericalConstants.H:46
NumericalDiffusion_Ymom
void NumericalDiffusion_Ymom(const Box &bx, const double dt, const Real num_diff_coeff, const Array4< const Real > &prim_data, const Array4< const Real > &cell_data, const Array4< Real > &rhs, const Array4< const Real > &mfx_arr, const Array4< const Real > &mfy_arr)
Definition: ERF_NumericalDiffusion.cpp:151
NumericalDiffusion_Xmom
void NumericalDiffusion_Xmom(const Box &bx, const double dt, const Real num_diff_coeff, const Array4< const Real > &prim_data, const Array4< const Real > &cell_data, const Array4< Real > &rhs, const Array4< const Real > &mfx_arr, const Array4< const Real > &mfy_arr)
Definition: ERF_NumericalDiffusion.cpp:85
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:32
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
PlaneAverage
Definition: ERF_PlaneAverage.H:14
eb_::get_w_const_factory
eb_aux_ const * get_w_const_factory() const noexcept
Return the ERF auxiliary z-face EB factory.
Definition: ERF_EB.H:133
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Return the ERF auxiliary y-face EB factory.
Definition: ERF_EB.H:131
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Return the ERF auxiliary x-face EB factory.
Definition: ERF_EB.H:129
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:79
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:237
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:235
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:238
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:236
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:47
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:218
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:217
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:219
WSM6SedCellScratch::ww
@ ww
Definition: ERF_AdvanceWSM6.cpp:107
module_mp_morr_two_moment::k1
real(c_double), private k1
Definition: ERF_module_mp_morr_two_moment.F90:213
mp_wdm6::r0
real(kind=kind_phys), parameter, private r0
Definition: ERF_module_mp_wdm6.F90:75
DampingChoice::rayleigh_damping_type
RayleighDampingType rayleigh_damping_type
Selected Rayleigh damping time-integration treatment.
Definition: ERF_DampingStruct.H:111
DampingChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Rayleigh damping inverse time scale [1/s].
Definition: ERF_DampingStruct.H:98
DampingChoice::rayleigh_damp_U
bool rayleigh_damp_U
Whether Rayleigh damping is applied to x-momentum.
Definition: ERF_DampingStruct.H:94
DampingChoice::rayleigh_damp_V
bool rayleigh_damp_V
Whether Rayleigh damping is applied to y-momentum.
Definition: ERF_DampingStruct.H:95
DampingChoice::rayleigh_damp_W
bool rayleigh_damp_W
Whether Rayleigh damping is applied to vertical momentum.
Definition: ERF_DampingStruct.H:96
InputSoundingData::input_sounding_time
amrex::Vector< amrex::Real > input_sounding_time
Physical time associated with each input sounding file.
Definition: ERF_InputSoundingData.H:496
InputSoundingData::V_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Device y velocity profiles.
Definition: ERF_InputSoundingData.H:509
InputSoundingData::U_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Device x velocity profiles.
Definition: ERF_InputSoundingData.H:509
InputSoundingData::tau_nudging
amrex::Real tau_nudging
Nudging time scale for relaxing toward input sounding data.
Definition: ERF_InputSoundingData.H:493
LargeScaleForcingData::get_forcing_time_coeffs
void get_forcing_time_coeffs(const amrex::Real &time, int &curr, int &next, amrex::Real &coeff_curr, amrex::Real &coeff_next) const
Determine interpolation coefficients to apply tendencies for the given time.
Definition: ERF_LargeScaleForcingData.H:335
LargeScaleForcingData::v_int_lsf_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > v_int_lsf_d
Definition: ERF_LargeScaleForcingData.H:457
LargeScaleForcingData::tau_lsf
amrex::Real tau_lsf
Definition: ERF_LargeScaleForcingData.H:407
LargeScaleForcingData::u_int_lsf_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > u_int_lsf_d
Definition: ERF_LargeScaleForcingData.H:457
SolverChoice::const_massflux_v
amrex::Real const_massflux_v
Target constant mass flux in the y direction.
Definition: ERF_DataStruct.H:2624
SolverChoice::init_type
static InitType init_type
Initial-condition source selected for the run.
Definition: ERF_DataStruct.H:2218
SolverChoice::coriolis_factor
amrex::Real coriolis_factor
Twice the planetary rotation rate used for Coriolis forcing.
Definition: ERF_DataStruct.H:2353
SolverChoice::flat_terrain
static bool flat_terrain
Permit fixed-index planar averages after validating a flat fitted mesh.
Definition: ERF_DataStruct.H:2227
SolverChoice::variable_coriolis
bool variable_coriolis
Whether spatially varying Coriolis forcing is enabled.
Definition: ERF_DataStruct.H:2570
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2572
SolverChoice::const_massflux_u
amrex::Real const_massflux_u
Target constant mass flux in the x direction.
Definition: ERF_DataStruct.H:2623
SolverChoice::forest_substep
bool forest_substep
Whether canopy source terms are applied only during substeps.
Definition: ERF_DataStruct.H:2314
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Numerical diffusion coefficient after input scaling.
Definition: ERF_DataStruct.H:2521
SolverChoice::hindcast_lateral_forcing
bool hindcast_lateral_forcing
Whether hindcast lateral forcing is enabled.
Definition: ERF_DataStruct.H:2633
SolverChoice::nudging_u_z2
amrex::Real nudging_u_z2
Definition: ERF_DataStruct.H:2377
SolverChoice::use_coriolis
bool use_coriolis
Whether Coriolis forcing is enabled.
Definition: ERF_DataStruct.H:2307
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:2236
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Sponge-layer options.
Definition: ERF_DataStruct.H:2251
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Whether a geostrophic wind profile has been configured.
Definition: ERF_DataStruct.H:2568
SolverChoice::dampingChoice
DampingChoice dampingChoice
Damping-related options.
Definition: ERF_DataStruct.H:2250
SolverChoice::nudging_u_z1
amrex::Real nudging_u_z1
Definition: ERF_DataStruct.H:2376
SolverChoice::do_forest_drag
bool do_forest_drag
Whether forest canopy drag is enabled.
Definition: ERF_DataStruct.H:2615
SolverChoice::massflux_khi
int massflux_khi
Upper vertical index for constant-mass-flux forcing.
Definition: ERF_DataStruct.H:2629
SolverChoice::large_scale_forcing
bool large_scale_forcing
Definition: ERF_DataStruct.H:2387
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:2224
SolverChoice::cosphi
amrex::Real cosphi
Cosine of the latitude used for Coriolis forcing.
Definition: ERF_DataStruct.H:2354
SolverChoice::buildings_type
static BuildingsType buildings_type
Building representation.
Definition: ERF_DataStruct.H:2230
SolverChoice::massflux_klo
int massflux_klo
Lower vertical index for constant-mass-flux forcing.
Definition: ERF_DataStruct.H:2628
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Applied geostrophic-wind forcing vector.
Definition: ERF_DataStruct.H:2566
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Whether custom vertical subsidence is enabled.
Definition: ERF_DataStruct.H:2360
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:2313
SolverChoice::const_massflux_tau
amrex::Real const_massflux_tau
Relaxation time scale for constant-mass-flux forcing.
Definition: ERF_DataStruct.H:2625
SolverChoice::sinphi
amrex::Real sinphi
Sine of the latitude used for Coriolis forcing.
Definition: ERF_DataStruct.H:2355
SolverChoice::do_mom_advection
bool do_mom_advection
Whether custom vertical subsidence is applied to momentum.
Definition: ERF_DataStruct.H:2362
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Whether custom forcing operates on primitive variables.
Definition: ERF_DataStruct.H:2364
SolverChoice::nudging_u
bool nudging_u
Definition: ERF_DataStruct.H:2383
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Whether solution fields are nudged toward input sounding data.
Definition: ERF_DataStruct.H:2370
SpongeChoice::sponge_type
static SpongeType sponge_type
Selected sponge damping model.
Definition: ERF_SpongeStruct.H:100
Here is the call graph for this function: