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, InputSoundingData &input_sounding_data, 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,
InputSoundingData &  input_sounding_data,
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 not fixed height so is not correct for variable dz
175  AMREX_ALWAYS_ASSERT(solverChoice.mesh_type != MeshType::VariableDz);
176 
177  const int offset = 1;
178  const int u_offset = 1;
179  const int v_offset = 1;
180 
181  //
182  // We use the alias here to control ncomp inside the PlaneAverage
183  //
184  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
185 
186  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = offset;
187  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
188  r_ave.compute_averages(ZDir(), r_ave.field());
189 
190  int ncell = r_ave.ncell_line();
191  Gpu::HostVector< Real> r_plane_h(ncell);
192  Gpu::DeviceVector< Real> r_plane_d(ncell);
193 
194  r_ave.line_average(Rho_comp, r_plane_h);
195 
196  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
197 
198  Real* dptr_r = r_plane_d.data();
199 
200  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
201  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
202 
203  dptr_r_plane = r_plane_tab.table();
204  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
205  {
206  dptr_r_plane(k-offset) = dptr_r[k];
207  });
208 
209  // U and V momentum
210  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
211  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
212 
213  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
214  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
215 
216  u_ave.compute_averages(ZDir(), u_ave.field());
217  v_ave.compute_averages(ZDir(), v_ave.field());
218 
219  int u_ncell = u_ave.ncell_line();
220  int v_ncell = v_ave.ncell_line();
221  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
222  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
223 
224  u_ave.line_average(0, u_plane_h);
225  v_ave.line_average(0, v_plane_h);
226 
227  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
228  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
229 
230  Real* dptr_u = u_plane_d.data();
231  Real* dptr_v = v_plane_d.data();
232 
233  Box udomain = domain; udomain.grow(2,ng_u[2]);
234  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
235  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
236  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
237 
238  dptr_u_plane = u_plane_tab.table();
239  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
240  {
241  dptr_u_plane(k-u_offset) = dptr_u[k];
242  });
243 
244  dptr_v_plane = v_plane_tab.table();
245  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
246  {
247  dptr_v_plane(k-v_offset) = dptr_v[k];
248  });
249 
250  // sum in z for massflux adjustment
251  if (enforce_massflux_x || enforce_massflux_y) {
252  Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
253  Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
254 
255  if (solverChoice.mesh_type == MeshType::ConstantDz) {
256  // note: massflux_khi corresponds to unstaggered indices in this case
257  rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
258  u_plane_h.begin() + u_offset + massflux_khi+1, zero);
259  rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
260  v_plane_h.begin() + v_offset + massflux_khi+1, zero);
261  rhoUA_target = std::accumulate(r_plane_h.begin() + offset + massflux_klo,
262  r_plane_h.begin() + offset + massflux_khi+1, zero);
263  rhoVA_target = rhoUA_target;
264 
265  rhoUA *= geom.CellSize(2) * Ly;
266  rhoVA *= geom.CellSize(2) * Lx;
267  rhoUA_target *= geom.CellSize(2) * Ly;
268  rhoVA_target *= geom.CellSize(2) * Lx;
269 
270  } else if (solverChoice.mesh_type == MeshType::StretchedDz) {
271  // note: massflux_khi corresponds to staggered indices in this case
272  for (int k=massflux_klo; k < massflux_khi; ++k) {
273  rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
274  rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
275  rhoUA_target += r_plane_h[k + offset] * stretched_dz_h[k];
276  }
277  rhoVA_target = rhoUA_target;
278 
279  rhoUA *= Ly;
280  rhoVA *= Lx;
281  rhoUA_target *= Ly;
282  rhoVA_target *= Lx;
283  }
284 
285  // at this point, this is integrated rho*dA
286  rhoUA_target *= U_target;
287  rhoVA_target *= V_target;
288 
289  Print() << "Integrated mass flux : " << rhoUA << " " << rhoVA
290  << " (target: " << rhoUA_target << " " << rhoVA_target << ")"
291  << std::endl;
292  }
293  }
294 
295  // *****************************************************************************
296  // Add all the other forcings
297  // *****************************************************************************
298  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
299  {
300  Box tbx = mfi.nodaltilebox(0);
301  Box tby = mfi.nodaltilebox(1);
302  Box tbz = mfi.nodaltilebox(2);
303  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
304 
305  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
306  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
307  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
308  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
309 
310  const Array4<const Real>& u = xvel.array(mfi);
311  const Array4<const Real>& v = yvel.array(mfi);
312  const Array4<const Real>& w = wvel.array(mfi);
313 
314  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
315  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
316  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
317 
318  const Array4<const Real>& r0 = r_hse.const_array(mfi);
319 
320  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
321  Array4<const Real>{};
322  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
323  Array4<const Real>{};
324  const Array4<const Real>& t_blank_xface_arr = (terrain_blank_xface) ? terrain_blank_xface->const_array(mfi) :
325  Array4<const Real>{};
326  const Array4<const Real>& t_blank_yface_arr = (terrain_blank_yface) ? terrain_blank_yface->const_array(mfi) :
327  Array4<const Real>{};
328  const Array4<const Real>& t_blank_zface_arr = (terrain_blank_zface) ? terrain_blank_zface->const_array(mfi) :
329  Array4<const Real>{};
330 
331  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
332  Array4<const Real>{};
333  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
334  Array4<const Real>{};
335 
336  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
337  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
338 
339 
340  // *****************************************************************************
341  // 1. Add Coriolis forcing (this assumes east is +x, north is +y)
342  // *****************************************************************************
343  if (use_coriolis && is_slow_step) {
344  if(solverChoice.init_type == InitType::HindCast) {
345  const Array4<const Real>& latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
346  ParallelFor(tbx, tby, tbz,
347  [=] AMREX_GPU_DEVICE (int i, int j, int k)
348  {
349  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
350  Real v_loc = fourth * (v(i,j+1,k) + v(i,j,k) + v(i-1,j+1,k) + v(i-1,j,k));
351  Real w_loc = fourth * (w(i,j,k+1) + w(i,j,k) + w(i-1,j,k+1) + w(i-1,j,k));
352  Real latitude = latlon_arr(i,j,k,0);
353  Real sphi_loc = std::sin(latitude*PI/Real(180.0));
354  Real cphi_loc = std::cos(latitude*PI/Real(180.0));
355  xmom_src_arr(i, j, k) += coriolis_factor * rho_on_u_face * (v_loc * sphi_loc - w_loc * cphi_loc);
356  },
357  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
358  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
359  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j-1,k) + u(i,j-1,k));
360  Real latitude = latlon_arr(i,j,k,0);
361  Real sphi_loc = std::sin(latitude*PI/Real(180.0));
362  ymom_src_arr(i, j, k) += -coriolis_factor * rho_on_v_face * u_loc * sphi_loc;
363  },
364  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
365  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
366  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j,k-1) + u(i,j,k-1));
367  Real latitude = latlon_arr(i,j,k,0);
368  Real cphi_loc = std::cos(latitude*PI/Real(180.0));
369  zmom_src_arr(i, j, k) += coriolis_factor * rho_on_w_face * u_loc * cphi_loc;
370  });
371  }
372  else if (var_coriolis && (sinPhi_mf) && (cosPhi_mf)) {
373  ParallelFor(tbx, tby, tbz,
374  [=] AMREX_GPU_DEVICE (int i, int j, int k)
375  {
376  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
377  Real v_loc = fourth * (v(i,j+1,k) + v(i,j,k) + v(i-1,j+1,k) + v(i-1,j,k));
378  Real w_loc = fourth * (w(i,j,k+1) + w(i,j,k) + w(i-1,j,k+1) + w(i-1,j,k));
379  Real sphi_loc = myhalf * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
380  Real cphi_loc = myhalf * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
381  xmom_src_arr(i, j, k) += coriolis_factor * rho_on_u_face * (v_loc * sphi_loc - w_loc * cphi_loc);
382  },
383  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
384  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
385  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j-1,k) + u(i,j-1,k));
386  Real sphi_loc = myhalf * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
387  ymom_src_arr(i, j, k) += -coriolis_factor * rho_on_v_face * u_loc * sphi_loc;
388  },
389  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
390  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
391  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j,k-1) + u(i,j,k-1));
392  zmom_src_arr(i, j, k) += coriolis_factor * rho_on_w_face * u_loc * cphi_arr(i,j,0);
393  });
394  } else {
395  if (solverChoice.terrain_type == TerrainType::EB) {
396  Array4<const Real> u_volfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
397  Array4<const Real> v_volfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
398  Array4<const Real> w_volfrac = (ebfact.get_w_const_factory())->getVolFrac().const_array(mfi);
399  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
400  Real v_loc = 0.0;
401  Real w_loc = 0.0;
402  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);
403  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);
404  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
405  if (v_vol > 0.0) {
406  v_loc = ( v_volfrac(i ,j+1,k) * v(i ,j+1,k) + v_volfrac(i ,j,k) * v(i ,j,k)
407  + 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;
408  }
409  if (w_vol > 0.0) {
410  w_loc = ( w_volfrac(i ,j,k+1) * w(i ,j,k+1) + w_volfrac(i ,j,k) * w(i ,j,k)
411  + 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;
412  }
413  xmom_src_arr(i, j, k) += coriolis_factor * rho_on_u_face * (v_loc * sinphi - w_loc * cosphi);
414  });
415  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
416  Real u_loc = 0.0;
417  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);
418  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
419  if (u_vol > 0.0) {
420  u_loc = ( u_volfrac(i+1,j ,k) * u(i+1,j ,k) + u_volfrac(i,j ,k) * u(i,j ,k)
421  + 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;
422  }
423  ymom_src_arr(i, j, k) += -coriolis_factor * rho_on_v_face * u_loc * sinphi;
424  });
425  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
426  Real u_loc = 0.0;
427  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);
428  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
429  if (u_vol > 0.0) {
430  u_loc = ( u_volfrac(i+1,j,k ) * u(i+1,j,k ) + u_volfrac(i,j,k) * u(i,j,k )
431  + 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;
432  }
433  zmom_src_arr(i, j, k) += coriolis_factor * rho_on_w_face * u_loc * cosphi;
434  });
435  } else {
436  ParallelFor(tbx, tby, tbz,
437  [=] AMREX_GPU_DEVICE (int i, int j, int k)
438  {
439  Real v_loc = fourth * (v(i,j+1,k) + v(i,j,k) + v(i-1,j+1,k) + v(i-1,j,k));
440  Real w_loc = fourth * (w(i,j,k+1) + w(i,j,k) + w(i-1,j,k+1) + w(i-1,j,k));
441  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
442  xmom_src_arr(i, j, k) += rho_on_u_face * ( coriolis_factor * (v_loc * sinphi - w_loc * cosphi) );
443  },
444  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
445  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
446  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j-1,k) + u(i,j-1,k));
447  ymom_src_arr(i, j, k) += rho_on_v_face * ( -coriolis_factor * u_loc * sinphi );
448  },
449  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
450  Real u_loc = fourth * (u(i+1,j,k) + u(i,j,k) + u(i+1,j,k-1) + u(i,j,k-1));
451  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
452  zmom_src_arr(i, j, k) += rho_on_w_face * ( coriolis_factor * u_loc * cosphi );
453  });
454  }
455  } // var_coriolis
456  } // use_coriolis
457 
458  // *****************************************************************************
459  // 2. Add RAYLEIGH damping
460  // *****************************************************************************
461  Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef;
462 
463  if ( (is_slow_step && !use_Rayleigh_fast_uv) || (!is_slow_step && use_Rayleigh_fast_uv)) {
464  if (rayleigh_damp_U) {
465  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
466  {
467  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
468  Real uu = rho_u(i,j,k) / rho_on_u_face;
469  Real sinesq = d_sinesq_at_lev[k];
470  xmom_src_arr(i, j, k) -= dampcoef*sinesq * (uu - ubar[k]) * rho_on_u_face;
471  });
472  }
473 
474  if (rayleigh_damp_V) {
475  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
476  {
477  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
478  Real vv = rho_v(i,j,k) / rho_on_v_face;
479  Real sinesq = d_sinesq_at_lev[k];
480  ymom_src_arr(i, j, k) -= dampcoef*sinesq * (vv - vbar[k]) * rho_on_v_face;
481  });
482  }
483  } // fast or slow step
484 
485  if ( (is_slow_step && !use_Rayleigh_fast_w) || (!is_slow_step && use_Rayleigh_fast_w)) {
486  if (rayleigh_damp_W) {
487  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
488  {
489  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
490  Real ww = rho_w(i,j,k) / rho_on_w_face;
491  Real sinesq = d_sinesq_stag_at_lev[k];
492  zmom_src_arr(i, j, k) -= dampcoef*sinesq * (ww - wbar[k]) * rho_on_w_face;
493  });
494  }
495  } // fast or slow step
496 
497  // *****************************************************************************
498  // 3a. Add constant GEOSTROPHIC forcing
499  // *****************************************************************************
500  if (is_slow_step) {
501  ParallelFor(tbx, tby, tbz,
502  [=] AMREX_GPU_DEVICE (int i, int j, int k)
503  {
504  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
505  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
506  },
507  [=] AMREX_GPU_DEVICE (int i, int j, int k)
508  {
509  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
510  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
511  },
512  [=] AMREX_GPU_DEVICE (int i, int j, int k)
513  {
514  Real rho_on_w_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
515  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
516  });
517  }
518 
519  // *****************************************************************************
520  // 3b. Add height-dependent GEOSTROPHIC forcing
521  // *****************************************************************************
522  if (geo_wind_profile && is_slow_step) {
523  ParallelFor(tbx, tby,
524  [=] AMREX_GPU_DEVICE (int i, int j, int k)
525  {
526  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
527  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
528  },
529  [=] AMREX_GPU_DEVICE (int i, int j, int k)
530  {
531  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
532  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
533  });
534  } // geo_wind_profile
535 
536  // *****************************************************************************
537  // 4. Add custom SUBSIDENCE terms
538  // *****************************************************************************
539  if (solverChoice.custom_w_subsidence && is_slow_step && solverChoice.do_mom_advection) {
540  if (solverChoice.custom_forcing_prim_vars) {
541  const int nr = Rho_comp;
542  ParallelFor(tbx, tby,
543  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
544  {
545  Real dzInv = myhalf*dxInv[2];
546  if (z_nd_arr) {
547  Real z_xf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
548  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
549  Real z_xf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
550  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
551  dzInv = one / (z_xf_hi - z_xf_lo);
552  }
553  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
554  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
555  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
556  Real wbar_xf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
557  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
558  },
559  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
560  {
561  Real dzInv = myhalf*dxInv[2];
562  if (z_nd_arr) {
563  Real z_yf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
564  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
565  Real z_yf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
566  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
567  dzInv = one / (z_yf_hi - z_yf_lo);
568  }
569  Real rho_on_v_face = myhalf * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
570  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
571  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
572  Real wbar_yf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
573  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
574  });
575  } else {
576  ParallelFor(tbx, tby,
577  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
578  {
579  Real dzInv = myhalf*dxInv[2];
580  if (z_nd_arr) {
581  Real z_xf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
582  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
583  Real z_xf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
584  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
585  dzInv = one / (z_xf_hi - z_xf_lo);
586  }
587  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
588  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
589  Real wbar_xf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
590  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
591  },
592  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
593  {
594  Real dzInv = myhalf*dxInv[2];
595  if (z_nd_arr) {
596  Real z_yf_lo = fourth * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
597  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
598  Real z_yf_hi = fourth * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
599  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
600  dzInv = one / (z_yf_hi - z_yf_lo);
601  }
602  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
603  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
604  Real wbar_yf = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
605  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
606  });
607  }
608  }
609 
610  if (solverChoice.large_scale_forcing && is_slow_step)
611  {
612  // subsidence terms for U and V
613  auto lsf_arr = lsf_tendencies->const_array(mfi);
614 
615  const int kmin = domain.smallEnd(2) + 1; // minimum k for vertical subsidence
616  const int kmax = domain.bigEnd(2) - 1; // maximum k for vertical subsidence
617 
618  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
619  {
620  if (k >= kmin && k <= kmax) {
621  int k1, k2;
622  if (lsf_arr(i, j, k, 2) >= 0.0)
623  {
624  k1 = k;
625  k2 = k-1;
626  } else {
627  k1 = k+1;
628  k2 = k;
629  }
630 
631  // one-sided difference over the single cell k1->k2
632  Real dzInv = dxInv[2];
633  if (z_nd_arr) {
634  Real z_uf_1 = fourth * ( z_nd_arr(i,j,k1 ) + z_nd_arr(i,j+1,k1 )
635  + z_nd_arr(i,j,k1+1) + z_nd_arr(i,j+1,k1+1) );
636  Real z_uf_2 = fourth * ( z_nd_arr(i,j,k2 ) + z_nd_arr(i,j+1,k2 )
637  + z_nd_arr(i,j,k2+1) + z_nd_arr(i,j+1,k2+1) );
638  dzInv = one / (z_uf_1 - z_uf_2);
639  }
640  Real rho_on_u_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
641  amrex::Real utend = -(dzInv * lsf_arr(i, j, k, 2)) * ( u(i, j, k1) - u(i, j, k2) );
642  xmom_src_arr(i, j, k) += utend * rho_on_u_face;
643  }
644  });
645 
646  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
647  {
648  if (k >= kmin && k <= kmax) {
649  int k1, k2;
650  if (lsf_arr(i, j, k, 2) >= 0.0)
651  {
652  k1 = k;
653  k2 = k-1;
654  } else {
655  k1 = k+1;
656  k2 = k;
657  }
658 
659  // one-sided difference over the single cell k1->k2
660  Real dzInv = dxInv[2];
661  if (z_nd_arr) {
662  Real z_vf_1 = fourth * ( z_nd_arr(i,j,k1 ) + z_nd_arr(i+1,j,k1 )
663  + z_nd_arr(i,j,k1+1) + z_nd_arr(i+1,j,k1+1) );
664  Real z_vf_2 = fourth * ( z_nd_arr(i,j,k2 ) + z_nd_arr(i+1,j,k2 )
665  + z_nd_arr(i,j,k2+1) + z_nd_arr(i+1,j,k2+1) );
666  dzInv = one / (z_vf_1 - z_vf_2);
667  }
668  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
669  amrex::Real vtend = -(dzInv * lsf_arr(i, j, k, 2)) * ( v(i, j, k1) - v(i, j, k2) );
670  ymom_src_arr(i, j, k) += vtend * rho_on_v_face;
671  }
672  });
673  }
674 
675  // *************************************************************************************
676  // 5. Add nudging towards value specified in input sounding
677  // NOTE: nudging_u only gates sounding-based wind nudging (the !large_scale_forcing
678  // branch below); LSF-based wind nudging (towards uls/vls) is controlled solely by
679  // large_scale_forcing, per the NOTE in Inputs.rst. Likewise, nudging_u_z1/z2 only
680  // restrict the sounding-based branch.
681  // *************************************************************************************
682  if (solverChoice.nudging_from_input_sounding && is_slow_step &&
683  (solverChoice.large_scale_forcing || solverChoice.nudging_u))
684  {
685  const bool l_lsf = solverChoice.large_scale_forcing;
686  const Real u_z1 = solverChoice.nudging_u_z1;
687  const Real u_z2 = solverChoice.nudging_u_z2;
688  Real uv_coeff_n = 1.0;
689  Real uv_coeff_np1 = 0.0;
690  Real tau = Real(1.0) / input_sounding_data.tau_nudging;
691  Real* u_nudge_n, *u_nudge_np1, *v_nudge_n, *v_nudge_np1;
692  if (!l_lsf)
693  {
694  int itime_n = 0;
695  int itime_np1 = 0;
696  int n_sounding_times = input_sounding_data.input_sounding_time.size();
697 
698  for (int nt = 1; nt < n_sounding_times; nt++) {
699  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
700  }
701  if (itime_n == n_sounding_times-1) {
702  itime_np1 = itime_n;
703  } else {
704  itime_np1 = itime_n+1;
705  uv_coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
706  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
707  uv_coeff_n = Real(1.0) - uv_coeff_np1;
708  }
709  u_nudge_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr() + 1;
710  u_nudge_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr() + 1;
711  v_nudge_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr() + 1;
712  v_nudge_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr() + 1;
713  } else {
714  int itime_curr = 0;
715  int itime_next = 0;
716  uv_coeff_n = 1.0;
717  uv_coeff_np1 = 0.0;
718  tau = 1.0 / lsf_data.tau_lsf; // only applies to u,v LSF nudging
719 
720  lsf_data.get_forcing_time_coeffs(time, itime_curr, itime_next, uv_coeff_n, uv_coeff_np1);
721  u_nudge_n = lsf_data.u_int_lsf_d[itime_curr].dataPtr();
722  u_nudge_np1 = lsf_data.u_int_lsf_d[itime_next].dataPtr();
723  v_nudge_n = lsf_data.v_int_lsf_d[itime_curr].dataPtr();
724  v_nudge_np1 = lsf_data.v_int_lsf_d[itime_next].dataPtr();
725  }
726 
727  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
728  {
729  Real z = z_cc_arr(i,j,k);
730  if (l_lsf || (z >= u_z1 && z <= u_z2)) {
731  Real unudge = -((dptr_u_plane(k)/dptr_r_plane(k)) - (uv_coeff_n*u_nudge_n[k] + uv_coeff_np1*u_nudge_np1[k]));
732  xmom_src_arr(i, j, k) += tau * unudge * dptr_r_plane(k);
733  }
734  });
735 
736  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
737  {
738  Real z = z_cc_arr(i,j,k);
739  if (l_lsf || (z >= u_z1 && z <= u_z2)) {
740  Real vnudge = -((dptr_v_plane(k)/dptr_r_plane(k)) - (uv_coeff_n*v_nudge_n[k] + uv_coeff_np1*v_nudge_np1[k]));
741  ymom_src_arr(i, j, k) += tau * vnudge * dptr_r_plane(k);
742  }
743  });
744  }
745 
746  // *****************************************************************************
747  // 6. Add NUMERICAL DIFFUSION terms
748  // *****************************************************************************
749 #if 0
750  if (l_use_ndiff) {
751  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
752  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
753  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
754  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
755  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
756  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
757  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
758  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
759  }
760 #endif
761 
762  // *****************************************************************************
763  // 7. Add SPONGING
764  // *****************************************************************************
765  if (is_slow_step) {
766  if (solverChoice.spongeChoice.sponge_type == SpongeType::Input_Sponge)
767  {
768  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
769  z_cc_arr, xmom_src_arr, ymom_src_arr,
770  rho_u, rho_v, d_sponge_ptrs_at_lev);
771  }
772  else if (solverChoice.spongeChoice.sponge_type == SpongeType::Standard)
773  {
774  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
775  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
776  r0, z_nd_arr, z_cc_arr);
777  }
778 
779  if(solverChoice.init_type == InitType::HindCast and solverChoice.hindcast_lateral_forcing){
780 
781  const Array4<const Real>& rho_u_forecast_state = (*forecast_state_at_lev)[IntVars::xmom].array(mfi);
782  const Array4<const Real>& rho_v_forecast_state = (*forecast_state_at_lev)[IntVars::ymom].array(mfi);
783  const Array4<const Real>& rho_w_forecast_state = (*forecast_state_at_lev)[IntVars::zmom].array(mfi);
784  const Array4<const Real>& cons_forecast_state = (*forecast_state_at_lev)[IntVars::cons].array(mfi);
785  ApplyBndryForcing_Forecast(solverChoice, geom, tbx, tby, tbz, z_nd_arr,
786  xmom_src_arr, ymom_src_arr, zmom_src_arr,
787  rho_u, rho_v, rho_w,
788  rho_u_forecast_state, rho_v_forecast_state, rho_w_forecast_state,
789  cons_forecast_state);
790  }
791  }
792 
793  // *****************************************************************************
794  // 8. Add CANOPY source terms
795  // *****************************************************************************
796  if (solverChoice.do_forest_drag &&
797  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
798  ParallelFor(tbx, tby, tbz,
799  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
800  {
801  const Real ux = u(i, j, k);
802  const Real uy = fourth * ( v(i, j , k ) + v(i-1, j , k )
803  + v(i, j+1, k ) + v(i-1, j+1, k ) );
804  const Real uz = fourth * ( w(i, j , k ) + w(i-1, j , k )
805  + w(i, j , k+1) + w(i-1, j , k+1) );
806  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
807  const Real f_drag = myhalf * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
808  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
809  },
810  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
811  {
812  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
813  + u(i+1, j , k ) + u(i+1, j-1, k ) );
814  const Real uy = v(i, j, k);
815  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
816  + w(i , j , k+1) + w(i , j-1, k+1) );
817  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
818  const Real f_drag = myhalf * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
819  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
820  },
821  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
822  {
823  const amrex::Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
824  + u(i , j , k-1) + u(i+1, j , k-1) );
825  const amrex::Real uy = fourth * ( v(i , j , k ) + v(i , j+1, k )
826  + v(i , j , k-1) + v(i , j+1, k-1) );
827  const amrex::Real uz = w(i, j, k);
828  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
829  const Real f_drag = myhalf * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
830  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
831  });
832  }
833  // *****************************************************************************
834  // 9a. Add immersed source terms for terrain
835  // *****************************************************************************
836  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
837  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
838 
839  ImmersedForcingTerrain_Xmom(tbx, u, v, w, cell_data, t_blank_arr, t_blank_xface_arr,
840  z_cc_arr, xmom_src_arr, geom, solverChoice, dt);
841  ImmersedForcingTerrain_Ymom(tby, u, v, w, cell_data, t_blank_arr, t_blank_yface_arr,
842  z_cc_arr, ymom_src_arr, geom, solverChoice, dt);
843  ImmersedForcingTerrain_Zmom(tbz, u, v, w, cell_data, t_blank_arr, t_blank_zface_arr,
844  z_cc_arr, zmom_src_arr, geom, solverChoice, dt);
845  }
846 
847  // *****************************************************************************
848  // 9b. Add immersed source terms for buildings
849  // *****************************************************************************
850  // geometric properties
851  const Real* dx_arr = geom.CellSize();
852  const Real delta_xy = std::sqrt(dx_arr[0] * dx_arr[1]);
853  if ((solverChoice.buildings_type == BuildingsType::ImmersedForcing) &&
854  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)) &&
855  (delta_xy <= 50.0)) { // only apply immersed forcing when grid spacing is less than 50m
856 
857  ImmersedForcingBuildings_Xmom(tbx, u, v, w, cell_data, t_blank_arr, t_blank_xface_arr,
858  z_cc_arr, xmom_src_arr, geom, solverChoice, dt);
859  ImmersedForcingBuildings_Ymom(tby, u, v, w, cell_data, t_blank_arr, t_blank_yface_arr,
860  z_cc_arr, ymom_src_arr, geom, solverChoice, dt);
861  ImmersedForcingBuildings_Zmom(tbz, u, v, w, cell_data, t_blank_arr, t_blank_zface_arr,
862  z_cc_arr, zmom_src_arr, geom, solverChoice, dt);
863  }
864 
865  // *****************************************************************************
866  // 10. Enforce constant mass flux
867  // *****************************************************************************
868  if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
869  Real tau_inv = one / solverChoice.const_massflux_tau;
870 
871  ParallelFor(tbx, tby,
872  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
873  xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
874  },
875  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
876  ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
877  });
878  }
879 
880  } // mfi
881 }
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:201
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
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
RayleighDampingType::FastExplicit
@ FastExplicit
RayleighDampingType::FastImplicit
@ FastImplicit
ubar
@ ubar
Definition: ERF_DataStruct.H:153
wbar
@ wbar
Definition: ERF_DataStruct.H:153
vbar
@ vbar
Definition: ERF_DataStruct.H:153
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:498
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:177
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:335
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:70
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:661
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:283
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
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);})
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:31
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:123
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:121
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:119
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:76
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:234
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:232
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:235
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:233
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:46
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
WSM6SedCellScratch::ww
@ ww
Definition: ERF_AdvanceWSM6.cpp:105
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::v_int_lsf_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > v_int_lsf_d
Definition: ERF_LargeScaleForcingData.H:451
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)
Determine interpolation coefficients to apply tendencies for the given time.
Definition: ERF_LargeScaleForcingData.H:329
LargeScaleForcingData::tau_lsf
amrex::Real tau_lsf
Definition: ERF_LargeScaleForcingData.H:401
LargeScaleForcingData::u_int_lsf_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > u_int_lsf_d
Definition: ERF_LargeScaleForcingData.H:451
SolverChoice::const_massflux_v
amrex::Real const_massflux_v
Target constant mass flux in the y direction.
Definition: ERF_DataStruct.H:2183
SolverChoice::init_type
static InitType init_type
Initial-condition source selected for the run.
Definition: ERF_DataStruct.H:1833
SolverChoice::coriolis_factor
amrex::Real coriolis_factor
Twice the planetary rotation rate used for Coriolis forcing.
Definition: ERF_DataStruct.H:1956
SolverChoice::variable_coriolis
bool variable_coriolis
Whether spatially varying Coriolis forcing is enabled.
Definition: ERF_DataStruct.H:2139
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2141
SolverChoice::const_massflux_u
amrex::Real const_massflux_u
Target constant mass flux in the x direction.
Definition: ERF_DataStruct.H:2182
SolverChoice::forest_substep
bool forest_substep
Whether canopy source terms are applied only during substeps.
Definition: ERF_DataStruct.H:1920
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Numerical diffusion coefficient after input scaling.
Definition: ERF_DataStruct.H:2121
SolverChoice::hindcast_lateral_forcing
bool hindcast_lateral_forcing
Whether hindcast lateral forcing is enabled.
Definition: ERF_DataStruct.H:2192
SolverChoice::nudging_u_z2
amrex::Real nudging_u_z2
Definition: ERF_DataStruct.H:1976
SolverChoice::use_coriolis
bool use_coriolis
Whether Coriolis forcing is enabled.
Definition: ERF_DataStruct.H:1913
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:1848
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Sponge-layer options.
Definition: ERF_DataStruct.H:1863
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Whether a geostrophic wind profile has been configured.
Definition: ERF_DataStruct.H:2137
SolverChoice::dampingChoice
DampingChoice dampingChoice
Damping-related options.
Definition: ERF_DataStruct.H:1862
SolverChoice::nudging_u_z1
amrex::Real nudging_u_z1
Definition: ERF_DataStruct.H:1975
SolverChoice::do_forest_drag
bool do_forest_drag
Whether forest canopy drag is enabled.
Definition: ERF_DataStruct.H:2174
SolverChoice::massflux_khi
int massflux_khi
Upper vertical index for constant-mass-flux forcing.
Definition: ERF_DataStruct.H:2188
SolverChoice::large_scale_forcing
bool large_scale_forcing
Definition: ERF_DataStruct.H:1986
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:1839
SolverChoice::cosphi
amrex::Real cosphi
Cosine of the latitude used for Coriolis forcing.
Definition: ERF_DataStruct.H:1957
SolverChoice::buildings_type
static BuildingsType buildings_type
Building representation.
Definition: ERF_DataStruct.H:1842
SolverChoice::massflux_klo
int massflux_klo
Lower vertical index for constant-mass-flux forcing.
Definition: ERF_DataStruct.H:2187
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Applied geostrophic-wind forcing vector.
Definition: ERF_DataStruct.H:2135
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Whether custom vertical subsidence is enabled.
Definition: ERF_DataStruct.H:1963
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:1919
SolverChoice::const_massflux_tau
amrex::Real const_massflux_tau
Relaxation time scale for constant-mass-flux forcing.
Definition: ERF_DataStruct.H:2184
SolverChoice::sinphi
amrex::Real sinphi
Sine of the latitude used for Coriolis forcing.
Definition: ERF_DataStruct.H:1958
SolverChoice::do_mom_advection
bool do_mom_advection
Whether custom vertical subsidence is applied to momentum.
Definition: ERF_DataStruct.H:1965
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Whether custom forcing operates on primitive variables.
Definition: ERF_DataStruct.H:1967
SolverChoice::nudging_u
bool nudging_u
Definition: ERF_DataStruct.H:1982
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Whether solution fields are nudged toward input sounding data.
Definition: ERF_DataStruct.H:1973
SpongeChoice::sponge_type
static SpongeType sponge_type
Selected sponge damping model.
Definition: ERF_SpongeStruct.H:100
Here is the call graph for this function: