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_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 (Real time, Real, 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 *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 MultiFab *surface_state_at_lev, InputSoundingData &input_sounding_data, bool is_slow_step)
 

Function Documentation

◆ make_mom_sources()

void make_mom_sources ( Real  time,
Real  ,
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 *  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 MultiFab *  surface_state_at_lev,
InputSoundingData &  input_sounding_data,
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
68 {
69  BL_PROFILE_REGION("erf_make_mom_sources()");
70 
71  Box domain(geom.Domain());
72  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
73 
74  // Initialize sources to zero each time we may use them
75  xmom_src.setVal(0.0);
76  ymom_src.setVal(0.0);
77  zmom_src.setVal(0.0);
78 
79  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
80 
81  // flags to apply certain source terms in substep call only
82  bool use_Rayleigh_fast_uv = ( (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit) ||
83  (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastImplicit) );
84  bool use_Rayleigh_fast_w = (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit);
85  bool use_canopy_fast = solverChoice.forest_substep;
86  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
87 
88  // *****************************************************************************
89  // Define source term for all three components of momenta from
90  // 1. Coriolis forcing for (xmom,ymom,zmom)
91  // 2. Rayleigh damping for (xmom,ymom,zmom)
92  // 3. Constant / height-dependent geostrophic forcing
93  // 4. Subsidence
94  // 5. Nudging towards input sounding data
95  // 6. Numerical diffusion for (xmom,ymom,zmom)
96  // 7. Sponge
97  // 8. Forest canopy
98  // 9a. Immersed forcing for terrain
99  // 9b. Immersed forcing for buildings
100  // 10. Constant mass flux
101  // *****************************************************************************
102  // NOTE: buoyancy is now computed in a separate routine - it should not appear here
103  // *****************************************************************************
104  //const bool l_use_ndiff = solverChoice.use_num_diff;
105 
106  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
107  if (solverChoice.do_forest_drag) {
108  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
109  }
110  }
111 
112 
113  // *****************************************************************************
114  // Data for Coriolis forcing
115  // *****************************************************************************
116  auto use_coriolis = solverChoice.use_coriolis;
117  auto coriolis_factor = solverChoice.coriolis_factor;
118  auto cosphi = solverChoice.cosphi;
119  auto sinphi = solverChoice.sinphi;
120  auto var_coriolis = solverChoice.variable_coriolis;
121 
122  // *****************************************************************************
123  // Flag for Geostrophic forcing
124  // *****************************************************************************
125  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
126  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
127 
128  // *****************************************************************************
129  // Data for Rayleigh damping
130  // *****************************************************************************
131  auto rayleigh_damp_U = solverChoice.dampingChoice.rayleigh_damp_U;
132  auto rayleigh_damp_V = solverChoice.dampingChoice.rayleigh_damp_V;
133  auto rayleigh_damp_W = solverChoice.dampingChoice.rayleigh_damp_W;
134 
135  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
136  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
137  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
138 
139  // *****************************************************************************
140  // Data for constant mass flux
141  // *****************************************************************************
142  bool enforce_massflux_x = (solverChoice.const_massflux_u != 0);
143  bool enforce_massflux_y = (solverChoice.const_massflux_v != 0);
144  Real U_target = solverChoice.const_massflux_u;
145  Real V_target = solverChoice.const_massflux_v;
146  int massflux_klo = solverChoice.massflux_klo;
147  int massflux_khi = solverChoice.massflux_khi;
148 
149  // These will be updated by integrating through the planar average profiles
150  Real rhoUA_target{0};
151  Real rhoVA_target{0};
152  Real rhoUA{0};
153  Real rhoVA{0};
154 
155  // *****************************************************************************
156  // Planar averages for subsidence, nudging, or constant mass flux
157  // *****************************************************************************
158  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
159  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
160 
161  if (is_slow_step && (dptr_wbar_sub || solverChoice.nudging_from_input_sounding ||
162  enforce_massflux_x || enforce_massflux_y))
163  {
164  const int offset = 1;
165  const int u_offset = 1;
166  const int v_offset = 1;
167 
168  //
169  // We use the alias here to control ncomp inside the PlaneAverage
170  //
171  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
172 
173  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = offset;
174  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
175  r_ave.compute_averages(ZDir(), r_ave.field());
176 
177  int ncell = r_ave.ncell_line();
178  Gpu::HostVector< Real> r_plane_h(ncell);
179  Gpu::DeviceVector< Real> r_plane_d(ncell);
180 
181  r_ave.line_average(Rho_comp, r_plane_h);
182 
183  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
184 
185  Real* dptr_r = r_plane_d.data();
186 
187  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
188  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
189 
190  dptr_r_plane = r_plane_tab.table();
191  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
192  {
193  dptr_r_plane(k-offset) = dptr_r[k];
194  });
195 
196  // U and V momentum
197  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
198  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
199 
200  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
201  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
202 
203  u_ave.compute_averages(ZDir(), u_ave.field());
204  v_ave.compute_averages(ZDir(), v_ave.field());
205 
206  int u_ncell = u_ave.ncell_line();
207  int v_ncell = v_ave.ncell_line();
208  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
209  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
210 
211  u_ave.line_average(0, u_plane_h);
212  v_ave.line_average(0, v_plane_h);
213 
214  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
215  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
216 
217  Real* dptr_u = u_plane_d.data();
218  Real* dptr_v = v_plane_d.data();
219 
220  Box udomain = domain; udomain.grow(2,ng_u[2]);
221  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
222  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
223  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
224 
225  dptr_u_plane = u_plane_tab.table();
226  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
227  {
228  dptr_u_plane(k-u_offset) = dptr_u[k];
229  });
230 
231  dptr_v_plane = v_plane_tab.table();
232  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
233  {
234  dptr_v_plane(k-v_offset) = dptr_v[k];
235  });
236 
237  // sum in z for massflux adjustment
238  if (enforce_massflux_x || enforce_massflux_y) {
239  Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
240  Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
241 
242  if (solverChoice.mesh_type == MeshType::ConstantDz) {
243  // note: massflux_khi corresponds to unstaggered indices in this case
244  rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
245  u_plane_h.begin() + u_offset + massflux_khi+1, 0.0);
246  rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
247  v_plane_h.begin() + v_offset + massflux_khi+1, 0.0);
248  rhoUA_target = std::accumulate(r_plane_h.begin() + offset + massflux_klo,
249  r_plane_h.begin() + offset + massflux_khi+1, 0.0);
250  rhoVA_target = rhoUA_target;
251 
252  rhoUA *= geom.CellSize(2) * Ly;
253  rhoVA *= geom.CellSize(2) * Lx;
254  rhoUA_target *= geom.CellSize(2) * Ly;
255  rhoVA_target *= geom.CellSize(2) * Lx;
256 
257  } else if (solverChoice.mesh_type == MeshType::StretchedDz) {
258  // note: massflux_khi corresponds to staggered indices in this case
259  for (int k=massflux_klo; k < massflux_khi; ++k) {
260  rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
261  rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
262  rhoUA_target += r_plane_h[k + offset] * stretched_dz_h[k];
263  }
264  rhoVA_target = rhoUA_target;
265 
266  rhoUA *= Ly;
267  rhoVA *= Lx;
268  rhoUA_target *= Ly;
269  rhoVA_target *= Lx;
270  }
271 
272  // at this point, this is integrated rho*dA
273  rhoUA_target *= U_target;
274  rhoVA_target *= V_target;
275 
276  Print() << "Integrated mass flux : " << rhoUA << " " << rhoVA
277  << " (target: " << rhoUA_target << " " << rhoVA_target << ")"
278  << std::endl;
279  }
280  }
281 
282  // *****************************************************************************
283  // Add all the other forcings
284  // *****************************************************************************
285  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
286  {
287  Box tbx = mfi.nodaltilebox(0);
288  Box tby = mfi.nodaltilebox(1);
289  Box tbz = mfi.nodaltilebox(2);
290  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
291 
292  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
293  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
294  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
295  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
296 
297  const Array4<const Real>& u = xvel.array(mfi);
298  const Array4<const Real>& v = yvel.array(mfi);
299  const Array4<const Real>& w = wvel.array(mfi);
300 
301  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
302  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
303  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
304 
305  const Array4<const Real>& r0 = r_hse.const_array(mfi);
306 
307  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
308  Array4<const Real>{};
309  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
310  Array4<const Real>{};
311 
312  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
313  Array4<const Real>{};
314  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
315  Array4<const Real>{};
316 
317  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
318  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
319 
320 
321  // *****************************************************************************
322  // 1. Add CORIOLIS forcing (this assumes east is +x, north is +y)
323  // *****************************************************************************
324  if (use_coriolis && is_slow_step) {
325  if(solverChoice.init_type == InitType::HindCast) {
326  const Array4<const Real>& latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
327  ParallelFor(tbx, tby, tbz,
328  [=] AMREX_GPU_DEVICE (int i, int j, int k)
329  {
330  Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k));
331  Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k));
332  Real latitude = latlon_arr(i,j,k,0);
333  Real sphi_loc = std::sin(latitude*PI/180.0);
334  Real cphi_loc = std::cos(latitude*PI/180.0);
335  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
336  },
337  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
338  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k));
339  Real latitude = latlon_arr(i,j,k,0);
340  Real sphi_loc = std::sin(latitude*PI/180.0);
341  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
342  },
343  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
344  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1));
345  Real latitude = latlon_arr(i,j,k,0);
346  Real cphi_loc = std::cos(latitude*PI/180.0);
347  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_loc;
348  });
349  }
350  else if (var_coriolis && (sinPhi_mf) && (cosPhi_mf)) {
351  ParallelFor(tbx, tby, tbz,
352  [=] AMREX_GPU_DEVICE (int i, int j, int k)
353  {
354  Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k));
355  Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i-1,j,k+1) + rho_w(i-1,j,k));
356  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
357  Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
358  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
359  },
360  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
361  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k));
362  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
363  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
364  },
365  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
366  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1));
367  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
368  });
369  } else {
370  ParallelFor(tbx, tby, tbz,
371  [=] AMREX_GPU_DEVICE (int i, int j, int k)
372  {
373  Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k));
374  Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i-1,j,k+1) + rho_w(i-1,j,k));
375  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
376  },
377  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
378  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k));
379  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
380  },
381  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
382  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1));
383  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
384  });
385  } // var_coriolis
386  } // use_coriolis
387 
388  // *****************************************************************************
389  // 2. Add RAYLEIGH damping
390  // *****************************************************************************
391  Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef;
392 
393  if ( (is_slow_step && !use_Rayleigh_fast_uv) || (!is_slow_step && use_Rayleigh_fast_uv)) {
394  if (rayleigh_damp_U) {
395  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
396  {
397  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
398  Real uu = rho_u(i,j,k) / rho_on_u_face;
399  Real sinesq = d_sinesq_at_lev[k];
400  xmom_src_arr(i, j, k) -= dampcoef*sinesq * (uu - ubar[k]) * rho_on_u_face;
401  });
402  }
403 
404  if (rayleigh_damp_V) {
405  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
406  {
407  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
408  Real vv = rho_v(i,j,k) / rho_on_v_face;
409  Real sinesq = d_sinesq_at_lev[k];
410  ymom_src_arr(i, j, k) -= dampcoef*sinesq * (vv - vbar[k]) * rho_on_v_face;
411  });
412  }
413  } // fast or slow step
414 
415  if ( (is_slow_step && !use_Rayleigh_fast_w) || (!is_slow_step && use_Rayleigh_fast_w)) {
416  if (rayleigh_damp_W) {
417  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
418  {
419  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
420  Real ww = rho_w(i,j,k) / rho_on_w_face;
421  Real sinesq = d_sinesq_stag_at_lev[k];
422  zmom_src_arr(i, j, k) -= dampcoef*sinesq * (ww - wbar[k]) * rho_on_w_face;
423  });
424  }
425  } // fast or slow step
426 
427  // *****************************************************************************
428  // 3a. Add constant GEOSTROPHIC forcing
429  // *****************************************************************************
430  if (is_slow_step) {
431  ParallelFor(tbx, tby, tbz,
432  [=] AMREX_GPU_DEVICE (int i, int j, int k)
433  {
434  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
435  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
436  },
437  [=] AMREX_GPU_DEVICE (int i, int j, int k)
438  {
439  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
440  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
441  },
442  [=] AMREX_GPU_DEVICE (int i, int j, int k)
443  {
444  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
445  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
446  });
447  }
448 
449  // *****************************************************************************
450  // 3b. Add height-dependent GEOSTROPHIC forcing
451  // *****************************************************************************
452  if (geo_wind_profile && is_slow_step) {
453  ParallelFor(tbx, tby,
454  [=] AMREX_GPU_DEVICE (int i, int j, int k)
455  {
456  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
457  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
458  },
459  [=] AMREX_GPU_DEVICE (int i, int j, int k)
460  {
461  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
462  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
463  });
464  } // geo_wind_profile
465 
466  // *****************************************************************************
467  // 4. Add custom SUBSIDENCE terms
468  // *****************************************************************************
469  if (solverChoice.custom_w_subsidence && is_slow_step && solverChoice.do_mom_advection) {
470  if (solverChoice.custom_forcing_prim_vars) {
471  const int nr = Rho_comp;
472  ParallelFor(tbx, tby,
473  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
474  {
475  Real dzInv = 0.5*dxInv[2];
476  if (z_nd_arr) {
477  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
478  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
479  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
480  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
481  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
482  }
483  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
484  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
485  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
486  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
487  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
488  },
489  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
490  {
491  Real dzInv = 0.5*dxInv[2];
492  if (z_nd_arr) {
493  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
494  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
495  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
496  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
497  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
498  }
499  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
500  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
501  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
502  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
503  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
504  });
505  } else {
506  ParallelFor(tbx, tby,
507  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
508  {
509  Real dzInv = 0.5*dxInv[2];
510  if (z_nd_arr) {
511  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
512  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
513  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
514  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
515  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
516  }
517  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
518  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
519  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
520  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
521  },
522  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
523  {
524  Real dzInv = 0.5*dxInv[2];
525  if (z_nd_arr) {
526  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
527  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
528  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
529  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
530  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
531  }
532  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
533  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
534  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
535  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
536  });
537  }
538  }
539 
540  // *************************************************************************************
541  // 5. Add nudging towards value specified in input sounding
542  // *************************************************************************************
543  if (solverChoice.nudging_from_input_sounding && is_slow_step)
544  {
545  int itime_n = 0;
546  int itime_np1 = 0;
547  Real coeff_n = Real(1.0);
548  Real coeff_np1 = Real(0.0);
549 
550  Real tau_inv = Real(1.0) / input_sounding_data.tau_nudging;
551 
552  int n_sounding_times = input_sounding_data.input_sounding_time.size();
553 
554  for (int nt = 1; nt < n_sounding_times; nt++) {
555  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
556  }
557  if (itime_n == n_sounding_times-1) {
558  itime_np1 = itime_n;
559  } else {
560  itime_np1 = itime_n+1;
561  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
562  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
563  coeff_n = Real(1.0) - coeff_np1;
564  }
565 
566  int nr = Rho_comp;
567 
568  const Real* u_inp_sound_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr();
569  const Real* u_inp_sound_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr();
570  const Real* v_inp_sound_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr();
571  const Real* v_inp_sound_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr();
572  ParallelFor(tbx, tby,
573  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
574  {
575  Real nudge_u = (coeff_n*u_inp_sound_n[k] + coeff_np1*u_inp_sound_np1[k]) - (dptr_u_plane(k)/dptr_r_plane(k));
576  nudge_u *= tau_inv;
577  xmom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_u;
578  },
579  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
580  {
581  Real nudge_v = (coeff_n*v_inp_sound_n[k] + coeff_np1*v_inp_sound_np1[k]) - (dptr_v_plane(k)/dptr_r_plane(k));
582  nudge_v *= tau_inv;
583  ymom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_v;
584  });
585  }
586 
587  // *****************************************************************************
588  // 6. Add NUMERICAL DIFFUSION terms
589  // *****************************************************************************
590 #if 0
591  if (l_use_ndiff) {
592  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
593  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
594  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
595  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
596  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
597  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
598  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
599  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
600  }
601 #endif
602 
603  // *****************************************************************************
604  // 7. Add SPONGING
605  // *****************************************************************************
606  if (is_slow_step) {
607  if (solverChoice.spongeChoice.sponge_type == "input_sponge")
608  {
609  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
610  z_cc_arr, xmom_src_arr, ymom_src_arr,
611  rho_u, rho_v, d_sponge_ptrs_at_lev);
612  }
613  else
614  {
615  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
616  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
617  r0, z_nd_arr, z_cc_arr);
618  }
619 
620  if(solverChoice.init_type == InitType::HindCast and solverChoice.hindcast_lateral_forcing){
621 
622  const Array4<const Real>& rho_u_forecast_state = (*forecast_state_at_lev)[IntVars::xmom].array(mfi);
623  const Array4<const Real>& rho_v_forecast_state = (*forecast_state_at_lev)[IntVars::ymom].array(mfi);
624  const Array4<const Real>& rho_w_forecast_state = (*forecast_state_at_lev)[IntVars::zmom].array(mfi);
625  const Array4<const Real>& cons_forecast_state = (*forecast_state_at_lev)[IntVars::cons].array(mfi);
626  ApplyBndryForcing_Forecast(solverChoice, geom, tbx, tby, tbz, z_nd_arr,
627  xmom_src_arr, ymom_src_arr, zmom_src_arr,
628  rho_u, rho_v, rho_w,
629  rho_u_forecast_state, rho_v_forecast_state, rho_w_forecast_state,
630  cons_forecast_state);
631  }
632  if(solverChoice.init_type == InitType::HindCast and solverChoice.hindcast_surface_bcs) {
633  const Array4<const Real>& surface_state_arr = (*surface_state_at_lev).array(mfi);
634  ApplySurfaceTreatment_BulkCoeff_Mom(tbx, tby,
635  xmom_src_arr, ymom_src_arr,
636  rho_u, rho_v,
637  cell_data, z_nd_arr,
638  surface_state_arr);
639  }
640  }
641 
642  // *****************************************************************************
643  // 8. Add CANOPY source terms
644  // *****************************************************************************
645  if (solverChoice.do_forest_drag &&
646  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
647  ParallelFor(tbx, tby, tbz,
648  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
649  {
650  const Real ux = u(i, j, k);
651  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
652  + v(i, j+1, k ) + v(i-1, j+1, k ) );
653  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
654  + w(i, j , k+1) + w(i-1, j , k+1) );
655  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
656  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
657  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
658  },
659  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
660  {
661  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
662  + u(i+1, j , k ) + u(i+1, j-1, k ) );
663  const Real uy = v(i, j, k);
664  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
665  + w(i , j , k+1) + w(i , j-1, k+1) );
666  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
667  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
668  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
669  },
670  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
671  {
672  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
673  + u(i , j , k-1) + u(i+1, j , k-1) );
674  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
675  + v(i , j , k-1) + v(i , j+1, k-1) );
676  const amrex::Real uz = w(i, j, k);
677  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
678  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
679  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
680  });
681  }
682  // *****************************************************************************
683  // 9a. Add immersed source terms for terrain
684  // *****************************************************************************
685  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
686  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
687  // geometric properties
688  const Real* dx_arr = geom.CellSize();
689  const Real dx_x = dx_arr[0];
690  const Real dx_y = dx_arr[1];
691  const Real dx_z = dx_arr[2];
692 
693  const Real alpha_m = solverChoice.if_Cd_momentum;
694  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
695  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
696  const Real U_s = 1.0; // unit velocity scale
697 
698  // MOST parameters
699  similarity_funs sfuns;
700  const Real ggg = CONST_GRAV;
701  const Real kappa = KAPPA;
702  const Real z0 = solverChoice.if_z0;
703  const Real tflux_in = solverChoice.if_surf_temp_flux;
704  const Real Olen_in = solverChoice.if_Olen_in;
705  const bool l_use_most = solverChoice.if_use_most;
706 
707  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
708  {
709  const Real ux = u(i, j, k);
710  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
711  + v(i, j+1, k ) + v(i-1, j+1, k ) );
712  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
713  + w(i, j , k+1) + w(i-1, j , k+1) );
714  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
715  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
716  const Real t_blank_above = 0.5 * (t_blank_arr(i, j, k+1) + t_blank_arr(i-1, j, k+1));
717  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
718  const Real rho_xface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
719 
720  if ((t_blank > 0 && (t_blank_above == 0.0)) && l_use_most) { // force to MOST value
721  // calculate tangential velocity one cell above
722  const Real ux2r = u(i, j, k+1) ;
723  const Real uy2r = 0.25 * ( v(i, j , k+1) + v(i-1, j , k+1)
724  + v(i, j+1, k+1) + v(i-1, j+1, k+1) ) ;
725  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
726 
727  // MOST
728  const Real theta_xface = (0.5 * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i-1,j,k, RhoTheta_comp))) / rho_xface;
729  const Real rho_xface_below = 0.5 * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp) );
730  const Real theta_xface_below = (0.5 * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i-1,j,k-1, RhoTheta_comp))) / rho_xface_below;
731  const Real theta_surf = theta_xface_below;
732 
733  Real psi_m = 0.0;
734  Real psi_h = 0.0;
735  Real ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
736  Real tflux = (tflux_in != 1e-8) ? tflux_in : -(theta_xface - theta_surf) * ustar * kappa / (std::log(1.5 * dx_z / z0) - psi_h);
737  Real Olen = (Olen_in != 1e-8) ? Olen_in : -ustar * ustar * ustar * theta_xface / (kappa * ggg * tflux + tiny);
738  Real zeta = 1.5 * dx_z / Olen;
739 
740  // similarity functions
741  psi_m = sfuns.calc_psi_m(zeta);
742  psi_h = sfuns.calc_psi_h(zeta);
743  ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m);
744 
745  // prevent some unphysical math
746  if (!(ustar > 0.0 && !std::isnan(ustar))) { ustar = 0.0; }
747  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
748  if (psi_m > std::log(0.5 * dx_z / z0)) { psi_m = std::log(0.5 * dx_z / z0); }
749 
750  // determine target velocity
751  const Real uTarget = ustar / kappa * (std::log(0.5 * dx_z / z0) - psi_m);
752  Real uxTarget = uTarget * ux2r / (tiny + h_windspeed2r);
753  const Real bc_forcing_x = -(uxTarget - ux); // BC forcing pushes nonrelative velocity toward target velocity
754  xmom_src_arr(i, j, k) -= (1-t_blank) * rho_xface * CdM * U_s * bc_forcing_x; // if Vf low, force more strongly to MOST. If high, less forcing.
755  } else {
756  xmom_src_arr(i, j, k) -= t_blank * rho_xface * CdM * ux * windspeed;
757  }
758  });
759  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
760  {
761  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
762  + u(i+1, j , k ) + u(i+1, j-1, k ) );
763  const Real uy = v(i, j, k);
764  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
765  + w(i , j , k+1) + w(i , j-1, k+1) );
766  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
767  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
768  const Real t_blank_above = 0.5 * (t_blank_arr(i, j, k+1) + t_blank_arr(i, j-1, k+1));
769  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
770  const Real rho_yface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
771 
772  if ((t_blank > 0 && (t_blank_above == 0.0)) && l_use_most) { // force to MOST value
773  // calculate tangential velocity one cell above
774  const Real ux2r = 0.25 * ( u(i , j , k+1) + u(i , j-1, k+1)
775  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
776  const Real uy2r = v(i, j, k+1) ;
777  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
778 
779  // MOST
780  const Real theta_yface = (0.5 * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k, RhoTheta_comp))) / rho_yface;
781  const Real rho_yface_below = 0.5 * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
782  const Real theta_yface_below = (0.5 * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i,j-1,k-1, RhoTheta_comp))) / rho_yface_below;
783  const Real theta_surf = theta_yface_below;
784 
785  Real psi_m = 0.0;
786  Real psi_h = 0.0;
787  Real ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
788  Real tflux = (tflux_in != 1e-8) ? tflux_in : -(theta_yface - theta_surf) * ustar * kappa / (std::log(1.5 * dx_z / z0) - psi_h);
789  Real Olen = (Olen_in != 1e-8) ? Olen_in : -ustar * ustar * ustar * theta_yface / (kappa * ggg * tflux + tiny);
790  Real zeta = 1.5 * dx_z / Olen;
791 
792  // similarity functions
793  psi_m = sfuns.calc_psi_m(zeta);
794  psi_h = sfuns.calc_psi_h(zeta);
795  ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m);
796 
797  // prevent some unphysical math
798  if (!(ustar > 0.0 && !std::isnan(ustar))) { ustar = 0.0; }
799  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
800  if (psi_m > std::log(0.5 * dx_z / z0)) { psi_m = std::log(0.5 * dx_z / z0); }
801 
802  // determine target velocity
803  const Real uTarget = ustar / kappa * (std::log(0.5 * dx_z / z0) - psi_m);
804  Real uyTarget = uTarget * uy2r / (tiny + h_windspeed2r);
805  const Real bc_forcing_y = -(uyTarget - uy); // BC forcing pushes nonrelative velocity toward target velocity
806  ymom_src_arr(i, j, k) -= (1 - t_blank) * rho_yface * CdM * U_s * bc_forcing_y; // if Vf low, force more strongly to MOST. If high, less forcing.
807  } else {
808  ymom_src_arr(i, j, k) -= t_blank * rho_yface * CdM * uy * windspeed;
809  }
810  });
811  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
812  {
813  const Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
814  + u(i , j , k-1) + u(i+1, j , k-1) );
815  const Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
816  + v(i , j , k-1) + v(i , j+1, k-1) );
817  const Real uz = w(i, j, k);
818  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
819  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
820  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
821  const Real rho_zface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
822  zmom_src_arr(i, j, k) -= t_blank * rho_zface * CdM * uz * windspeed;
823  });
824  }
825 
826  // *****************************************************************************
827  // 9b. Add immersed source terms for buildings
828  // *****************************************************************************
829  if ((solverChoice.buildings_type == BuildingsType::ImmersedForcing ) &&
830  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)))
831  {
832  // geometric properties
833  const Real* dx_arr = geom.CellSize();
834  const Real dx_x = dx_arr[0];
835  const Real dx_y = dx_arr[1];
836 
837  const Real alpha_m = solverChoice.if_Cd_momentum;
838  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
839  const Real min_t_blank = 0.005; // threshold for where immersed forcing acts
840 
841  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
842  {
843  const Real ux = u(i, j, k);
844  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
845  + v(i, j+1, k ) + v(i-1, j+1, k ) );
846  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
847  + w(i, j , k+1) + w(i-1, j , k+1) );
848  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
849 
850  Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
851  if (t_blank < min_t_blank) { t_blank = 0.0; }
852  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
853  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
854  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
855  const Real rho_xface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
856  xmom_src_arr(i, j, k) -= t_blank * rho_xface * CdM * ux * windspeed;
857  });
858  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
859  {
860  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
861  + u(i+1, j , k ) + u(i+1, j-1, k ) );
862  const Real uy = v(i, j, k);
863  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
864  + w(i , j , k+1) + w(i , j-1, k+1) );
865  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
866 
867  Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
868  if (t_blank < min_t_blank) { t_blank = 0.0; }
869  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
870  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
871  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
872  const Real rho_yface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
873  ymom_src_arr(i, j, k) -= t_blank * rho_yface * CdM * uy * windspeed;
874  });
875  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
876  {
877  const Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
878  + u(i , j , k-1) + u(i+1, j , k-1) );
879  const Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
880  + v(i , j , k-1) + v(i , j+1, k-1) );
881  const Real uz = w(i, j, k);
882  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
883 
884  Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
885  if (t_blank < min_t_blank) { t_blank = 0.0; }
886  const Real dx_z = (z_nd_arr) ? (z_nd_arr(i,j,k) - z_nd_arr(i,j,k-1)) : dx_arr[2]; // ASW double check
887  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
888  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
889  const Real rho_zface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
890  zmom_src_arr(i, j, k) -= t_blank * rho_zface * CdM * uz * windspeed;
891  });
892  }
893 
894  // *****************************************************************************
895  // 10. Enforce constant mass flux
896  // *****************************************************************************
897  if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
898  Real tau_inv = Real(1.0) / solverChoice.const_massflux_tau;
899 
900  ParallelFor(tbx, tby,
901  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
902  xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
903  },
904  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
905  ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
906  });
907  }
908 
909  } // mfi
910 }
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)
Definition: ERF_ApplySpongeZoneBCs.cpp:118
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:8
ApplySurfaceTreatment_BulkCoeff_Mom
void ApplySurfaceTreatment_BulkCoeff_Mom(const Box &tbx, const Box &tby, 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 Array4< const Real > &cons_state, const Array4< const Real > &z_phys_nd, const Array4< const Real > &surface_state_arr)
Definition: ERF_ApplySurfaceTreatment_BulkCoeff.cpp:8
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
RayleighDampingType::FastExplicit
@ FastExplicit
RayleighDampingType::FastImplicit
@ FastImplicit
ubar
@ ubar
Definition: ERF_DataStruct.H:98
wbar
@ wbar
Definition: ERF_DataStruct.H:98
vbar
@ vbar
Definition: ERF_DataStruct.H:98
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:16
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+0.5) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
NumericalDiffusion_Ymom
void NumericalDiffusion_Ymom(const Box &bx, const Real 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 Real 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
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
PlaneAverage
Definition: ERF_PlaneAverage.H:14
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:63
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:176
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:174
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:177
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:175
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:45
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:157
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:156
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:158
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
DampingChoice::rayleigh_damp_V
bool rayleigh_damp_V
Definition: ERF_DampingStruct.H:85
DampingChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Definition: ERF_DampingStruct.H:88
DampingChoice::rayleigh_damp_W
bool rayleigh_damp_W
Definition: ERF_DampingStruct.H:86
DampingChoice::rayleigh_damping_type
RayleighDampingType rayleigh_damping_type
Definition: ERF_DampingStruct.H:101
DampingChoice::rayleigh_damp_U
bool rayleigh_damp_U
Definition: ERF_DampingStruct.H:84
InputSoundingData::V_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Definition: ERF_InputSoundingData.H:407
InputSoundingData::input_sounding_time
amrex::Vector< amrex::Real > input_sounding_time
Definition: ERF_InputSoundingData.H:394
InputSoundingData::tau_nudging
amrex::Real tau_nudging
Definition: ERF_InputSoundingData.H:391
InputSoundingData::U_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Definition: ERF_InputSoundingData.H:407
SolverChoice::do_mom_advection
bool do_mom_advection
Definition: ERF_DataStruct.H:1161
SolverChoice::coriolis_factor
amrex::Real coriolis_factor
Definition: ERF_DataStruct.H:1152
SolverChoice::mesh_type
static MeshType mesh_type
Definition: ERF_DataStruct.H:1068
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Definition: ERF_DataStruct.H:1132
SolverChoice::if_use_most
bool if_use_most
Definition: ERF_DataStruct.H:1136
SolverChoice::dampingChoice
DampingChoice dampingChoice
Definition: ERF_DataStruct.H:1078
SolverChoice::const_massflux_v
amrex::Real const_massflux_v
Definition: ERF_DataStruct.H:1230
SolverChoice::if_z0
amrex::Real if_z0
Definition: ERF_DataStruct.H:1131
SolverChoice::cosphi
amrex::Real cosphi
Definition: ERF_DataStruct.H:1153
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:1202
SolverChoice::hindcast_lateral_forcing
bool hindcast_lateral_forcing
Definition: ERF_DataStruct.H:1239
SolverChoice::massflux_klo
int massflux_klo
Definition: ERF_DataStruct.H:1234
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Definition: ERF_DataStruct.H:1159
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:1169
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:1124
SolverChoice::sinphi
amrex::Real sinphi
Definition: ERF_DataStruct.H:1154
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:1204
SolverChoice::const_massflux_u
amrex::Real const_massflux_u
Definition: ERF_DataStruct.H:1229
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Definition: ERF_DataStruct.H:1135
SolverChoice::use_coriolis
bool use_coriolis
Definition: ERF_DataStruct.H:1118
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:1188
SolverChoice::variable_coriolis
bool variable_coriolis
Definition: ERF_DataStruct.H:1206
SolverChoice::if_Cd_momentum
amrex::Real if_Cd_momentum
Definition: ERF_DataStruct.H:1129
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:1163
SolverChoice::buildings_type
static BuildingsType buildings_type
Definition: ERF_DataStruct.H:1059
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:1056
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:1079
SolverChoice::init_type
static InitType init_type
Definition: ERF_DataStruct.H:1050
SolverChoice::hindcast_surface_bcs
bool hindcast_surface_bcs
Definition: ERF_DataStruct.H:1240
SolverChoice::do_forest_drag
bool do_forest_drag
Definition: ERF_DataStruct.H:1226
SolverChoice::const_massflux_tau
amrex::Real const_massflux_tau
Definition: ERF_DataStruct.H:1231
SolverChoice::massflux_khi
int massflux_khi
Definition: ERF_DataStruct.H:1235
SolverChoice::forest_substep
bool forest_substep
Definition: ERF_DataStruct.H:1125
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:1208
SpongeChoice::sponge_type
std::string sponge_type
Definition: ERF_SpongeStruct.H:58
similarity_funs
Definition: ERF_MOSTStress.H:40
similarity_funs::calc_psi_m
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:90
similarity_funs::calc_psi_h
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:104
Here is the call graph for this function: