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 <limits>
#include "ERF_ImmersedForcing.H"
#include "ERF_NumericalDiffusion.H"
#include "ERF_PlaneAverage.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_SrcHeaders.H"
#include "ERF_Utils.H"
Include dependency graph for ERF_MakeMomSources.cpp:

Functions

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

Function Documentation

◆ make_mom_sources()

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

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

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