ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MakeSources.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_SrcHeaders.H>
#include <ERF_TI_slow_headers.H>
#include <ERF_MOSTStress.H>
Include dependency graph for ERF_MakeSources.cpp:

Functions

void make_sources (int level, int, double dt, double time_d, const Vector< MultiFab > &S_data, const MultiFab &S_prim, MultiFab &source, const MultiFab &base_state, const MultiFab *z_phys_cc, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, const MultiFab *qheating_rates, MultiFab *terrain_blank, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &mapfac, const MultiFab *rhotheta_src, const MultiFab *rhoqt_src, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const Real *d_sinesq_at_lev, InputSoundingData &input_sounding_data, TurbulentPerturbation &turbPert, bool is_slow_step)
 

Function Documentation

◆ make_sources()

void make_sources ( int  level,
int  ,
double  dt,
double  time_d,
const Vector< MultiFab > &  S_data,
const MultiFab &  S_prim,
MultiFab &  source,
const MultiFab &  base_state,
const MultiFab *  z_phys_cc,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  zvel,
const MultiFab *  qheating_rates,
MultiFab *  terrain_blank,
const Geometry  geom,
const SolverChoice &  solverChoice,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
const MultiFab *  rhotheta_src,
const MultiFab *  rhoqt_src,
const Real *  dptr_wbar_sub,
const Vector< Real * >  d_rayleigh_ptrs_at_lev,
const Real *  d_sinesq_at_lev,
InputSoundingData &  input_sounding_data,
TurbulentPerturbation &  turbPert,
bool  is_slow_step 
)

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

Parameters
[in]levellevel of resolution
[in]nrkwhich RK stage
[in]dtslow time step
[in]S_datacurrent solution
[in]S_primprimitive variables (i.e. conserved variables divided by density)
[in]sourcesource terms for conserved variables
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]mapfacmap factors
[in]dptr_rhotheta_srccustom temperature source term
[in]dptr_rhoqt_srccustom moisture source term
[in]dptr_wbar_subsubsidence source term
[in]d_rayleigh_ptrs_at_levVector of {strength of Rayleigh damping, reference value of theta} used to define Rayleigh damping
[in]d_sinesq_at_levsin( (pi/2) (z-z_t)/(damping depth)) at cell centers
59 {
60  BL_PROFILE_REGION("erf_make_sources()");
61 
62  Real time = static_cast<Real>(time_d);
63 
64  // *****************************************************************************
65  // Initialize source to zero since we re-compute it every RK stage
66  // *****************************************************************************
67  if (is_slow_step) {
68  source.setVal(0.);
69  } else {
70  source.setVal(0.0,Rho_comp,2);
71  }
72 
73  const bool l_use_ndiff = solverChoice.use_num_diff;
74 
75  TurbChoice tc = solverChoice.turbChoice[level];
76  const bool l_use_KE = tc.use_tke;
77  const bool l_diff_KE = tc.diffuse_tke_3D;
78 
79  const Box& domain = geom.Domain();
80 
81  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
82  const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
83 
84  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
85  MultiFab th_hse (base_state, make_alias, BaseState::th0_comp , 1);
86  MultiFab qv_hse (base_state, make_alias, BaseState::qv0_comp , 1);
87 
88  Real* thetabar = d_rayleigh_ptrs_at_lev[Rayleigh::thetabar];
89 
90  // flags to apply certain source terms in substep call only
91  bool use_Rayleigh_fast = ( (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit) ||
92  (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastImplicit) );
93  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
94 
95  // flag for a moisture model
96  bool has_moisture = (solverChoice.moisture_type != MoistureType::None);
97 
98  // *****************************************************************************
99  // Planar averages for subsidence terms
100  // *****************************************************************************
101  Table1D<Real> dptr_r_plane, dptr_t_plane, dptr_qv_plane, dptr_qc_plane;
102  TableData<Real, 1> r_plane_tab, t_plane_tab, qv_plane_tab, qc_plane_tab;
103  bool compute_averages = ( is_slow_step && (dptr_wbar_sub || solverChoice.nudging_from_input_sounding) );
104 
105  if (compute_averages)
106  {
107  // The plane averaging operates at fixed z not fixed height so is not correct for variable dz
108  AMREX_ALWAYS_ASSERT(solverChoice.mesh_type != MeshType::VariableDz);
109 
110  //
111  // The call to "compute_averages" currently does all the components in one call
112  // We can then extract each component separately with the "line_average" call
113  //
114  // We need just one ghost cell in the vertical
115  //
116  IntVect ng_c(S_data[IntVars::cons].nGrowVect()); ng_c[2] = 1;
117  //
118  // With no moisture we only (rho) and (rho theta); with moisture we also do qv and qc
119  // We use the alias here to control ncomp inside the PlaneAverage
120  //
121  int ncomp = (!has_moisture) ? 2 : RhoQ2_comp+1;
122  MultiFab cons(S_data[IntVars::cons], make_alias, 0, ncomp);
123 
124  PlaneAverage cons_ave(&cons, geom, solverChoice.ave_plane, ng_c);
125  cons_ave.compute_averages(ZDir(), cons_ave.field());
126 
127  int ncell = cons_ave.ncell_line();
128 
129  Gpu::HostVector< Real> r_plane_h(ncell);
130  Gpu::DeviceVector< Real> r_plane_d(ncell);
131 
132  Gpu::HostVector< Real> t_plane_h(ncell);
133  Gpu::DeviceVector< Real> t_plane_d(ncell);
134 
135  cons_ave.line_average(Rho_comp , r_plane_h);
136  cons_ave.line_average(RhoTheta_comp, t_plane_h);
137 
138  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
139  Gpu::copyAsync(Gpu::hostToDevice, t_plane_h.begin(), t_plane_h.end(), t_plane_d.begin());
140 
141  Real* dptr_r = r_plane_d.data();
142  Real* dptr_t = t_plane_d.data();
143 
144  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
145  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
146  t_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
147 
148  int offset = ng_c[2];
149 
150  dptr_r_plane = r_plane_tab.table();
151  dptr_t_plane = t_plane_tab.table();
152  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
153  {
154  dptr_r_plane(k-offset) = dptr_r[k];
155  dptr_t_plane(k-offset) = dptr_t[k];
156  });
157 
158  if (has_moisture)
159  {
160  Gpu::HostVector< Real> qv_plane_h(ncell), qc_plane_h(ncell);
161  Gpu::DeviceVector<Real> qv_plane_d(ncell), qc_plane_d(ncell);
162 
163  // Water vapor
164  cons_ave.line_average(RhoQ1_comp, qv_plane_h);
165  Gpu::copyAsync(Gpu::hostToDevice, qv_plane_h.begin(), qv_plane_h.end(), qv_plane_d.begin());
166 
167  // Cloud water
168  cons_ave.line_average(RhoQ2_comp, qc_plane_h);
169  Gpu::copyAsync(Gpu::hostToDevice, qc_plane_h.begin(), qc_plane_h.end(), qc_plane_d.begin());
170 
171  Real* dptr_qv = qv_plane_d.data();
172  Real* dptr_qc = qc_plane_d.data();
173 
174  qv_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
175  qc_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
176 
177  dptr_qv_plane = qv_plane_tab.table();
178  dptr_qc_plane = qc_plane_tab.table();
179  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
180  {
181  dptr_qv_plane(k-offset) = dptr_qv[k];
182  dptr_qc_plane(k-offset) = dptr_qc[k];
183  });
184  }
185  }
186 
187  // *****************************************************************************
188  // Radiation flux vector for four stream approximation
189  // *****************************************************************************
190  // NOTE: The fluxes live on w-faces
191  int klo = domain.smallEnd(0);
192  int khi = domain.bigEnd(2);
193  int nk = khi - klo + 2;
194  Gpu::DeviceVector<Real> radiation_flux(nk,zero);
195  Gpu::DeviceVector<Real> q_integral(nk,zero);
196  Real* rad_flux = radiation_flux.data();
197  Real* q_int = q_integral.data();
198 
199  // *****************************************************************************
200  // Define source term for cell-centered conserved variables, from
201  // one user-defined source terms for (rho theta) and (rho q_t)
202  // two radiation for (rho theta)
203  // three Rayleigh damping for (rho theta)
204  // Real(4.) custom forcing for (rho theta) and (rho Q1)
205  // Real(5.) custom subsidence for (rho theta) and (rho Q1)
206  // Real(6.) numerical diffusion for (rho theta)
207  // Real(7.) sponging
208  // Real(8.) turbulent perturbation
209  // Real(9.) nudging towards input sounding values (only for theta)
210  // 10a. Immersed forcing for terrain
211  // 10b. Immersed forcing for buildings
212  // Real(11.) Four stream radiation source for (rho theta)
213  // *****************************************************************************
214 
215  // ***********************************************************************************************
216  // Add remaining source terms
217  // ***********************************************************************************************
218 #ifdef _OPENMP
219 #pragma omp parallel if (Gpu::notInLaunchRegion())
220 #endif
221  {
222  for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
223  {
224  Box bx = mfi.tilebox();
225 
226  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
227  const Array4<const Real>& cell_prim = S_prim.array(mfi);
228  const Array4<Real> & cell_src = source.array(mfi);
229 
230  const Array4<const Real>& r0 = r_hse.const_array(mfi);
231  const Array4<const Real>& th0 = th_hse.const_array(mfi);
232  const Array4<const Real>& qv0 = qv_hse.const_array(mfi);
233 
234  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
235 
236  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
237  Array4<const Real>{};
238 
239 
240  // *************************************************************************************
241  // two Add radiation source terms to (rho theta)
242  // *************************************************************************************
243  if (solverChoice.rad_type != RadiationType::None && is_slow_step) {
244  auto const& qheating_arr = qheating_rates->const_array(mfi);
245  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
246  {
247  // Short-wavelength and long-wavelength radiation source terms
248  cell_src(i,j,k,RhoTheta_comp) += cell_data(i,j,k,Rho_comp) * ( qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1) );
249  });
250  }
251 
252 
253  // *************************************************************************************
254  // three Add Rayleigh damping for (rho theta)
255  // *************************************************************************************
256  Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef;
257 
258  if (solverChoice.dampingChoice.rayleigh_damp_T) {
259  if ((is_slow_step && !use_Rayleigh_fast) || (!is_slow_step && use_Rayleigh_fast)) {
260  int n = RhoTheta_comp;
261  int nr = Rho_comp;
262  int np = PrimTheta_comp;
263 
264  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
265  {
266  Real theta = cell_prim(i,j,k,np);
267  Real sinesq = d_sinesq_at_lev[k];
268  cell_src(i, j, k, n) -= dampcoef*sinesq * (theta - thetabar[k]) * cell_data(i,j,k,nr);
269  });
270  }
271  }
272 
273  // *************************************************************************************
274  // Real(4.) Add custom forcing for (rho theta)
275  // *************************************************************************************
276  if (solverChoice.custom_rhotheta_forcing && is_slow_step) {
277  const int n = RhoTheta_comp;
278  auto const& rhotheta_src_arr = rhotheta_src->const_array(mfi);
279  if (solverChoice.spatial_rhotheta_forcing)
280  {
281  if (solverChoice.custom_forcing_prim_vars) {
282  const int nr = Rho_comp;
283  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
284  {
285  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhotheta_src_arr(i, j, k);
286  });
287  } else {
288  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
289  {
290  cell_src(i, j, k, n) += rhotheta_src_arr(i, j, k);
291  });
292  }
293  } else {
294  if (solverChoice.custom_forcing_prim_vars) {
295  const int nr = Rho_comp;
296  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
297  {
298  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhotheta_src_arr(0, 0, k);
299  });
300  } else {
301  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
302  {
303  cell_src(i, j, k, n) += rhotheta_src_arr(0, 0, k);
304  });
305  }
306  }
307  }
308 
309  // *************************************************************************************
310  // Real(4.) Add custom forcing for RhoQ1
311  // *************************************************************************************
312  if (solverChoice.custom_moisture_forcing && is_slow_step) {
313  const int n = RhoQ1_comp;
314  auto const& rhoqt_src_arr = rhoqt_src->const_array(mfi);
315  if (solverChoice.spatial_moisture_forcing)
316  {
317  if (solverChoice.custom_forcing_prim_vars) {
318  const int nr = Rho_comp;
319  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
320  {
321  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhoqt_src_arr(i, j, k);
322  });
323  } else {
324  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
325  {
326  cell_src(i, j, k, n) += rhoqt_src_arr(i, j, k);
327  });
328  }
329  } else {
330  if (solverChoice.custom_forcing_prim_vars) {
331  const int nr = Rho_comp;
332  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
333  {
334  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhoqt_src_arr(0, 0, k);
335  });
336  } else {
337  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
338  {
339  cell_src(i, j, k, n) += rhoqt_src_arr(0, 0, k);
340  });
341  }
342  }
343  }
344 
345  // *************************************************************************************
346  // Real(5.) Add custom subsidence for (rho theta)
347  // *************************************************************************************
348  if (solverChoice.custom_w_subsidence && is_slow_step && solverChoice.do_theta_advection) {
349  const int n = RhoTheta_comp;
350  if (solverChoice.custom_forcing_prim_vars) {
351  const int nr = Rho_comp;
352  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
353  {
354  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
355  Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
356  Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
357  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
358  cell_src(i, j, k, n) -= cell_data(i,j,k,nr) * wbar_cc * (T_hi - T_lo) * dzInv;
359  });
360  } else {
361  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
362  {
363  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
364  Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
365  Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
366  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
367  cell_src(i, j, k, n) -= wbar_cc * (T_hi - T_lo) * dzInv;
368  });
369  }
370  }
371 
372  // *************************************************************************************
373  // Real(5.) Add custom subsidence for RhoQ1 and RhoQ2
374  // *************************************************************************************
375  if (solverChoice.custom_w_subsidence && (solverChoice.moisture_type != MoistureType::None) && is_slow_step) {
376  const int nv = RhoQ1_comp;
377  if (solverChoice.custom_forcing_prim_vars) {
378  const int nr = Rho_comp;
379  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
380  {
381  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
382  Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
383  Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
384  Real Qc_hi = dptr_qc_plane(k+1) / dptr_r_plane(k+1);
385  Real Qc_lo = dptr_qc_plane(k-1) / dptr_r_plane(k-1);
386  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
387  cell_src(i, j, k, nv ) -= cell_data(i,j,k,nr) * wbar_cc * (Qv_hi - Qv_lo) * dzInv;
388  cell_src(i, j, k, nv+1) -= cell_data(i,j,k,nr) * wbar_cc * (Qc_hi - Qc_lo) * dzInv;
389  });
390  } else {
391  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
392  {
393  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
394  Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
395  Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
396  Real Qc_hi = dptr_qc_plane(k+1) / dptr_r_plane(k+1);
397  Real Qc_lo = dptr_qc_plane(k-1) / dptr_r_plane(k-1);
398  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
399  cell_src(i, j, k, nv ) -= wbar_cc * (Qv_hi - Qv_lo) * dzInv;
400  cell_src(i, j, k, nv+1) -= wbar_cc * (Qc_hi - Qc_lo) * dzInv;
401  });
402  }
403  }
404 
405  // *************************************************************************************
406  // Real(6.) Add numerical diffusion for rho and (rho theta)
407  // *************************************************************************************
408  if (l_use_ndiff && is_slow_step)
409  {
410  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
411  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
412 
413  // Rho is a special case
414  NumericalDiffusion_Scal(bx, 0, 1, dt, solverChoice.num_diff_coeff,
415  cell_data, cell_data, cell_src, mf_mx, mf_my);
416 
417  // Other scalars proceed as normal
418  NumericalDiffusion_Scal(bx, 1, 1, dt, solverChoice.num_diff_coeff,
419  cell_prim, cell_data, cell_src, mf_mx, mf_my);
420 
421 
422  if (l_use_KE && l_diff_KE) {
423  NumericalDiffusion_Scal(bx, RhoKE_comp, 1, dt, solverChoice.num_diff_coeff,
424  cell_prim, cell_data, cell_src, mf_mx, mf_my);
425  }
426 
427  NumericalDiffusion_Scal(bx, RhoScalar_comp, NSCALARS, dt, solverChoice.num_diff_coeff,
428  cell_prim, cell_data, cell_src, mf_mx, mf_my);
429  }
430 
431  // *************************************************************************************
432  // Real(7.) Add sponging
433  // *************************************************************************************
434  if ( (solverChoice.spongeChoice.sponge_type != SpongeType::None) && is_slow_step) {
435  const int n_qstate = S_data[IntVars::cons].nComp() - (NDRY + NSCALARS);
436  ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data, r0, th0, qv0, z_cc_arr, n_qstate);
437  }
438 
439  // *************************************************************************************
440  // Real(8.) Add perturbation
441  // *************************************************************************************
442  if (solverChoice.use_source_perturbation(level) && is_slow_step) {
443  auto m_ixtype = S_data[IntVars::cons].boxArray().ixType(); // Conserved term
444  const amrex::Array4<const amrex::Real>& pert_cell = turbPert.pb_cell[level].const_array(mfi);
445  turbPert.apply_tpi(level, bx, RhoTheta_comp, m_ixtype, cell_src, pert_cell); // Applied as source term
446  }
447 
448  // *************************************************************************************
449  // Real(9.) Add nudging towards value specified in input sounding
450  // *************************************************************************************
451  if (solverChoice.nudging_from_input_sounding && is_slow_step)
452  {
453  int itime_n = 0;
454  int itime_np1 = 0;
455  Real coeff_n = one;
456  Real coeff_np1 = zero;
457 
458  Real tau_inv = one / input_sounding_data.tau_nudging;
459 
460  int n_sounding_times = input_sounding_data.input_sounding_time.size();
461 
462  for (int nt = 1; nt < n_sounding_times; nt++) {
463  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
464  }
465  if (itime_n == n_sounding_times-1) {
466  itime_np1 = itime_n;
467  } else {
468  itime_np1 = itime_n+1;
469  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
470  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
471  coeff_n = one - coeff_np1;
472  }
473 
474  const Real* theta_inp_sound_n = input_sounding_data.theta_inp_sound_d[itime_n].dataPtr();
475  const Real* theta_inp_sound_np1 = input_sounding_data.theta_inp_sound_d[itime_np1].dataPtr();
476 
477  const int n = RhoTheta_comp;
478  const int nr = Rho_comp;
479 
480  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
481  {
482  Real nudge = (coeff_n*theta_inp_sound_n[k] + coeff_np1*theta_inp_sound_np1[k]) - (dptr_t_plane(k)/dptr_r_plane(k));
483  nudge *= tau_inv;
484  cell_src(i, j, k, n) += cell_data(i, j, k, nr) * nudge;
485  });
486  }
487 
488  // *************************************************************************************
489  // 10a. Add immersed source terms for terrain
490  // *************************************************************************************
491  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
492  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)))
493  {
494  const Array4<const Real>& u = xvel.array(mfi);
495  const Array4<const Real>& v = yvel.array(mfi);
496 
497  // geometric properties
498  const Real* dx_arr = geom.CellSize();
499  const Real dx_x = dx_arr[0];
500  const Real dx_y = dx_arr[1];
501 
502  const Real alpha_h = solverChoice.if_Cd_scalar;
503  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
504  const Real U_s = one; // unit velocity scale
505 
506  // MOST parameters
507  similarity_funs sfuns;
508  const Real ggg = CONST_GRAV;
509  const Real kappa = KAPPA;
510  const Real z0 = solverChoice.if_z0;
511  const Real tflux = solverChoice.if_surf_temp_flux;
512  const Real init_surf_temp = solverChoice.if_init_surf_temp;
513 
514  // Note this has been converted to K / s when it was read in;
515  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
516 
517  const Real Olen_in = solverChoice.if_Olen_in;
518 
519  ParallelFor(bx, [=]
520  AMREX_GPU_DEVICE(int i, int j, int k) noexcept
521  {
522  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
523  const Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
524 
525  const Real t_blank = t_blank_arr(i, j, k);
526  const Real t_blank_above = t_blank_arr(i, j, k+1);
527  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
528  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
529  const Real h_windspeed2r = std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r);
530 
531  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
532  const Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
533 
534  // SURFACE TEMP AND HEATING/COOLING RATE
535  if (init_surf_temp > zero) {
536  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
537  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
538  const Real bc_forcing_rt_srf = -(cell_data(i,j,k-1,Rho_comp) * surf_temp - cell_data(i,j,k-1,RhoTheta_comp));
539  cell_src(i, j, k-1, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf; // k-1
540  }
541  }
542 
543  // SURFACE HEAT FLUX
544  if (tflux != Real(1e-8)){
545  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
546  Real psi_m = zero;
547  Real psi_h = zero;
548  Real psi_h_neighbor = zero;
549  Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
550  const Real Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
551  const Real zeta = (myhalf) * dx_z / Olen;
552  const Real zeta_neighbor = (Real(1.5)) * dx_z / Olen;
553 
554  // similarity functions
555  psi_m = sfuns.calc_psi_m(zeta);
556  psi_h = sfuns.calc_psi_h(zeta);
557  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
558  ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
559 
560  // prevent some unphysical math
561  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
562  if (!(ustar < two && !std::isnan(ustar))) { ustar = two; }
563  if (psi_h_neighbor > std::log(Real(1.5) * dx_z / z0)) { psi_h_neighbor = std::log(Real(1.5) * dx_z / z0); }
564  if (psi_h > std::log(myhalf * dx_z / z0)) { psi_h = std::log(myhalf * dx_z / z0); }
565 
566  // We do not know the actual temperature so use cell above
567  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
568  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
569  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
570 
571  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
572  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
573  }
574  }
575 
576  // OBUKHOV LENGTH
577  if (Olen_in != Real(1e-8)){
578  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
579  const Real Olen = Olen_in;
580  const Real zeta = (myhalf) * dx_z / Olen;
581  const Real zeta_neighbor = (Real(1.5)) * dx_z / Olen;
582 
583  // similarity functions
584  const Real psi_m = sfuns.calc_psi_m(zeta);
585  const Real psi_h = sfuns.calc_psi_h(zeta);
586  const Real psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
587  const Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
588 
589  // We do not know the actual temperature so use cell above
590  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
591  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
592  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
593 
594  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
595  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
596  }
597  }
598 
599  });
600  }
601 
602  // *************************************************************************************
603  // 10b. Add immersed source terms for buildings
604  // *************************************************************************************
605  // geometric properties
606  const Real* dx_arr = geom.CellSize();
607  const Real dx_x = dx_arr[0];
608  const Real dx_y = dx_arr[1];
609  const Real delta_xy = std::pow(dx_x*dx_y, myhalf);
610  if ((solverChoice.buildings_type == BuildingsType::ImmersedForcing ) &&
611  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)) &&
612  (delta_xy <= 50.0)) // only apply immersed forcing when grid spacing is less than 50m
613  {
614  const Array4<const Real>& u = xvel.array(mfi);
615  const Array4<const Real>& v = yvel.array(mfi);
616  const Array4<const Real>& w = zvel.array(mfi);
617 
618  const Real alpha_h = solverChoice.if_Cd_scalar;
619  const Real U_s = one; // unit velocity scale
620  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
621  const Real min_t_blank = Real(1.e-4);
622 
623  // MOST parameters
624  similarity_funs sfuns;
625  const Real ggg = CONST_GRAV;
626  const Real kappa = KAPPA;
627  const Real z0 = solverChoice.if_z0;
628  const Real tflux = solverChoice.if_surf_temp_flux;
629  const Real init_surf_temp = solverChoice.if_init_surf_temp;
630  const Real surf_heating_rate = solverChoice.if_surf_heating_rate; // Note this has been converted to K / s when it was read in;
631  const Real Olen_in = solverChoice.if_Olen_in;
632 
633  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
634  {
635  Real t_blank = t_blank_arr(i, j, k);
636  Real t_blank_below = t_blank_arr(i, j, k-1);
637  Real t_blank_above = t_blank_arr(i, j, k+1);
638  Real t_blank_north = t_blank_arr(i , j+1, k);
639  Real t_blank_south = t_blank_arr(i , j-1, k);
640  Real t_blank_east = t_blank_arr(i+1, j , k);
641  Real t_blank_west = t_blank_arr(i-1, j , k);
642  if (t_blank < min_t_blank) { t_blank = zero; } // deal with situations where very small volfrac exist
643  if (t_blank_below < min_t_blank) { t_blank_below = zero; }
644  if (t_blank_north < min_t_blank) { t_blank_north = zero; }
645  if (t_blank_south < min_t_blank) { t_blank_south = zero; }
646  if (t_blank_east < min_t_blank) { t_blank_east = zero; }
647  if (t_blank_west < min_t_blank) { t_blank_west = zero; }
648 
649  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
650  Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
651 
652  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
653  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
654  const Real h_windspeed2r = std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r);
655 
656  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
657  Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
658 
659  // SURFACE TEMP AND HEATING/COOLING RATE
660  if (init_surf_temp > zero) {
661  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
662  if (t_blank > 0 && (t_blank_above == zero) && (t_blank_below == one)) { // building roof
663  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
664  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
665 
666  } else if (((t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) ||
667  (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero) ||
668  (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero) ||
669  (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero))) {
670  // this should enter for just building walls
671  // walls are currently separated to allow for flexibility in the future to heat walls differently
672 
673  // south face
674  if ((t_blank < t_blank_north) && (t_blank_north == one)) {
675  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
676  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
677  }
678 
679  // north face
680  if ((t_blank < t_blank_south) && (t_blank_south == one)) {
681  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
682  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
683  }
684 
685  // west face
686  if ((t_blank < t_blank_east) && (t_blank_east == one)) {
687  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
688  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
689  }
690 
691  // east face
692  if ((t_blank < t_blank_west) && (t_blank_west == one)) {
693  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
694  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
695  }
696 
697  }
698  }
699 
700  // SURFACE HEAT FLUX
701  if (tflux != Real(1.e-8)){
702  if (t_blank > zero && (t_blank_above == zero)) { // building roof
703  Real psi_m = zero;
704  Real psi_h = zero;
705  Real psi_h_neighbor = zero;
706  Real ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
707  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
708 
709  for (int iter = 0; iter < 2; ++iter) {
710  if (iter > 0) { Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny); }
711  Real zeta = (myhalf) * dx_z / Olen;
712  Real zeta_neighbor = (1.5) * dx_z / Olen;
713 
714  // similarity functions
715  psi_m = sfuns.calc_psi_m(zeta);
716  psi_h = sfuns.calc_psi_h(zeta);
717  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
718  ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
719  }
720 
721  // prevent some unphysical math
722  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
723  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
724  if (psi_h_neighbor > std::log(1.5 * dx_z / z0)) { psi_h_neighbor = std::log(1.5 * dx_z / z0); }
725  if (psi_h > std::log(myhalf * dx_z / z0)) { psi_h = std::log(myhalf * dx_z / z0); }
726 
727  // We do not know the actual temperature so use cell above
728  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
729  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * dx_z / z0) - psi_h_neighbor);
730  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
731 
732  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
733  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
734 
735  } else if (((t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) ||
736  (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero) ||
737  (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero) ||
738  (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero))) { // this should enter for just building walls
739 
740  Real ux_cellaway = zero;
741  Real uy_cellaway = zero;
742  Real uz_cellaway = zero;
743  Real u1 = zero;
744  Real u2 = zero;
745  Real delta = zero;
746 
747  // south face
748  if (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero) {
749  ux_cellaway = myhalf * (u(i ,j-1,k) + u(i+1,j-1,k ));
750  uz_cellaway = myhalf * (w(i ,j-1,k) + w(i ,j-1,k+1));
751  u1 = ux_cellaway;
752  u2 = uz_cellaway;
753  delta = dx_y;
754 
755  // MOST
756  theta_neighbor = cell_data(i,j-1,k,RhoTheta_comp) / cell_data(i,j-1,k,Rho_comp);
757  }
758 
759  // north face
760  if (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero) {
761  ux_cellaway = myhalf * (u(i ,j+1,k) + u(i+1,j+1,k ));
762  uz_cellaway = myhalf * (w(i ,j+1,k) + w(i ,j+1,k+1));
763  u1 = ux_cellaway;
764  u2 = uz_cellaway;
765  delta = dx_y;
766 
767  // MOST
768  theta_neighbor = cell_data(i,j+1,k,RhoTheta_comp) / cell_data(i,j+1,k,Rho_comp);
769  }
770 
771  // west face
772  if (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero) {
773  uy_cellaway = myhalf * (u(i-1,j ,k) + u(i-1,j+1,k ));
774  uz_cellaway = myhalf * (w(i-1,j ,k) + w(i-1,j ,k+1));
775  u1 = uy_cellaway;
776  u2 = uz_cellaway;
777  delta = dx_x;
778 
779  // MOST
780  theta_neighbor = cell_data(i-1,j,k,RhoTheta_comp) / cell_data(i-1,j,k,Rho_comp);
781  }
782 
783  // east face
784  if (t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) {
785  uy_cellaway = myhalf * (u(i+1,j ,k) + u(i+1,j+1,k ));
786  uz_cellaway = myhalf * (w(i+1,j ,k) + w(i+1,j ,k+1));
787  u1 = uy_cellaway;
788  u2 = uz_cellaway;
789  delta = dx_x;
790 
791  // MOST
792  theta_neighbor = cell_data(i+1,j,k,RhoTheta_comp) / cell_data(i+1,j,k,Rho_comp);
793  }
794 
795  Real tan_wspd = std::sqrt(u1 * u1 + u2 * u2);
796 
797  Real psi_m = zero;
798  Real psi_h = zero;
799  Real psi_h_neighbor = zero;
800  Real ustar = tan_wspd * kappa / (std::log(1.5 * delta / z0) - psi_m);
801  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
802 
803  for (int iter = 0; iter < 2; ++iter) {
804  if (iter > 0) { Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny); }
805  Real zeta = (myhalf) * delta / Olen;
806  Real zeta_neighbor = (1.5) * delta / Olen;
807 
808  // similarity functions
809  psi_m = sfuns.calc_psi_m(zeta);
810  psi_h = sfuns.calc_psi_h(zeta);
811  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
812  ustar = tan_wspd * kappa / (std::log((1.5) * delta / z0) - psi_m);
813  }
814 
815  // prevent some unphysical math
816  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
817  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
818  if (psi_h_neighbor > std::log(1.5 * delta / z0)) { psi_h_neighbor = std::log(1.5 * delta / z0); }
819  if (psi_h > std::log(myhalf * delta / z0)) { psi_h = std::log(myhalf * delta / z0); }
820 
821  // We do not know the actual temperature so use cell above
822  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
823  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * delta / z0) - psi_h_neighbor);
824  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * delta / z0) - psi_h);
825 
826  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
827  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
828  }
829  }
830  });
831  }
832 
833  // *************************************************************************************
834  // Real(11.) Add 4 stream radiation src to RhoTheta
835  // *************************************************************************************
836  if (solverChoice.four_stream_radiation && has_moisture && is_slow_step)
837  {
838  AMREX_ALWAYS_ASSERT((bx.smallEnd(2) == klo) && (bx.bigEnd(2) == khi));
839  Real D = Real(3.75e-6); // [s^-1]
840  Real F0 = 70; // [W/m^2]
841  Real F1 = 22; // [W/m^2]
842  Real krad = 85; // [m^2 kg^-1]
843  Real qt_i = Real(0.008);
844 
845  Box xybx = makeSlab(bx,2,klo);
846  ParallelFor(xybx, [=]
847  AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
848  {
849  // Inclusive scan at w-faces for the Q integral (also find "i" values)
850  q_int[0] = zero;
851  Real zi = myhalf * (z_cc_arr(i,j,khi) + z_cc_arr(i,j,khi-1));
852  Real rhoi = myhalf * (cell_data(i,j,khi,Rho_comp) + cell_data(i,j,khi-1,Rho_comp));
853  for (int k(klo+1); k<=khi+1; ++k) {
854  int lk = k - klo;
855  // Average to w-faces when looping w-faces
856  Real dz = (z_cc_arr) ? myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-2)) : dx[2];
857  q_int[lk] = q_int[lk-1] + krad * cell_data(i,j,k-1,Rho_comp) * cell_data(i,j,k-1,RhoQ2_comp) * dz;
858  Real qt_hi = cell_data(i,j,k ,RhoQ1_comp) + cell_data(i,j,k ,RhoQ2_comp);
859  Real qt_lo = cell_data(i,j,k-1,RhoQ1_comp) + cell_data(i,j,k-1,RhoQ2_comp);
860  if ( (qt_lo > qt_i) && (qt_hi < qt_i) ) {
861  zi = myhalf * (z_cc_arr(i,j,k) + z_cc_arr(i,j,k-1));
862  rhoi = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp));
863  }
864  }
865 
866  // Decompose the integral to get the fluxes at w-faces
867  Real q_int_inf = q_int[khi+1];
868  for (int k(klo); k<=khi+1; ++k) {
869  int lk = k - klo;
870  Real z = myhalf * (z_cc_arr(i,j,k) + z_cc_arr(i,j,k-1));
871  rad_flux[lk] = F1*std::exp(-q_int[lk]) + F0*std::exp(-(q_int_inf - q_int[lk]));
872  if (z > zi) {
873  rad_flux[lk] += rhoi * Cp_d * D * ( std::pow(z-zi,Real(4.)/three)/Real(4.) + zi*std::pow(z-zi,one/three) ) ;
874  }
875  }
876 
877  // Compute the radiative heating source
878  for (int k(klo); k<=khi; ++k) {
879  int lk = k - klo;
880  // Average to w-faces when looping CC
881  Real dzInv = (z_cc_arr) ? one/ (myhalf * (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1))) : dxInv[2];
882  // NOTE: Fnet = Up - Dn (all fluxes are up here)
883  // dT/dt = dF/dz * (1/(-rho*Cp))
884  Real dTdt = (rad_flux[lk+1] - rad_flux[lk]) * dzInv / (-cell_data(i,j,k,Rho_comp)*Cp_d);
885  Real qv = cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp);
886  Real iexner = one/getExnergivenRTh(cell_data(i,j,k,RhoTheta_comp), RdoCp, qv);
887  // Convert dT/dt to dTheta/dt and multiply rho
888  cell_src(i,j,k,RhoTheta_comp) += cell_data(i,j,k,Rho_comp) * dTdt * iexner;
889  }
890  });
891  }
892  } // mfi
893  } // OMP
894 }
ApplySpongeZoneBCsForCC
void ApplySpongeZoneBCsForCC(const SpongeChoice &spongeChoice, const Geometry geom, const Box &bx, const Array4< Real > &cell_rhs, const Array4< const Real > &cell_data, const Array4< const Real > &r0, const Array4< const Real > &th0, const Array4< const Real > &qv0, const Array4< const Real > &z_phys_cc, int n_qstate)
Definition: ERF_ApplySpongeZoneBCs.cpp:7
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:63
Cp_d
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:49
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
RdoCp
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:54
RayleighDampingType::FastExplicit
@ FastExplicit
RayleighDampingType::FastImplicit
@ FastImplicit
thetabar
@ thetabar
Definition: ERF_DataStruct.H:152
m_y
@ m_y
Definition: ERF_DataStruct.H:28
m_x
@ m_x
Definition: ERF_DataStruct.H:27
Coord::z
@ z
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
getExnergivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:156
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
NDRY
#define NDRY
Definition: ERF_IndexDefines.H:13
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
PrimTheta_comp
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:55
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
NumericalDiffusion_Scal
void NumericalDiffusion_Scal(const Box &bx, const int start_comp, const int num_comp, 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:18
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:31
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
PlaneAverage
Definition: ERF_PlaneAverage.H:14
BaseState::qv0_comp
@ qv0_comp
Definition: ERF_IndexDefines.H:77
BaseState::th0_comp
@ th0_comp
Definition: ERF_IndexDefines.H:76
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:73
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:194
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:20
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:30
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:46
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:177
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:176
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:179
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:178
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
WSM6SedNodeScratch::zi
@ zi
Definition: ERF_AdvanceWSM6.cpp:133
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
module_mp_morr_two_moment::rhoi
real(c_double), private rhoi
Definition: ERF_module_mp_morr_two_moment.F90:188
mp_wsm6::r0
real(kind=kind_phys), parameter, private r0
Definition: ERF_module_mp_wsm6.F90:21
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_T
bool rayleigh_damp_T
Whether Rayleigh damping is applied to potential temperature.
Definition: ERF_DampingStruct.H:97
InputSoundingData::input_sounding_time
amrex::Vector< amrex::Real > input_sounding_time
Physical time associated with each input sounding file.
Definition: ERF_InputSoundingData.H:462
InputSoundingData::theta_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > theta_inp_sound_d
Device potential temperature profiles [K].
Definition: ERF_InputSoundingData.H:475
InputSoundingData::tau_nudging
amrex::Real tau_nudging
Nudging time scale for relaxing toward input sounding data.
Definition: ERF_InputSoundingData.H:459
SolverChoice::turbChoice
amrex::Vector< TurbChoice > turbChoice
Turbulence options for each AMR level.
Definition: ERF_DataStruct.H:1393
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:1604
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Input Obukhov length for immersed-forcing MOST [m].
Definition: ERF_DataStruct.H:1462
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:1618
SolverChoice::if_z0
amrex::Real if_z0
Immersed-forcing roughness length [m].
Definition: ERF_DataStruct.H:1458
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Numerical diffusion coefficient after input scaling.
Definition: ERF_DataStruct.H:1601
SolverChoice::do_theta_advection
bool do_theta_advection
Whether custom vertical subsidence is applied to rho-theta.
Definition: ERF_DataStruct.H:1493
SolverChoice::spatial_moisture_forcing
bool spatial_moisture_forcing
Whether spatially varying moisture forcing is enabled.
Definition: ERF_DataStruct.H:1498
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:1377
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Sponge-layer options.
Definition: ERF_DataStruct.H:1392
SolverChoice::four_stream_radiation
bool four_stream_radiation
Whether the four-stream radiation approximation is enabled.
Definition: ERF_DataStruct.H:1445
SolverChoice::dampingChoice
DampingChoice dampingChoice
Damping-related options.
Definition: ERF_DataStruct.H:1391
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:1368
SolverChoice::if_Cd_scalar
amrex::Real if_Cd_scalar
Immersed-forcing drag coefficient for scalars.
Definition: ERF_DataStruct.H:1453
SolverChoice::buildings_type
static BuildingsType buildings_type
Building representation.
Definition: ERF_DataStruct.H:1371
SolverChoice::custom_rhotheta_forcing
bool custom_rhotheta_forcing
Whether custom rho-theta forcing is enabled.
Definition: ERF_DataStruct.H:1490
SolverChoice::if_init_surf_temp
amrex::Real if_init_surf_temp
Initial immersed-forcing surface temperature [K].
Definition: ERF_DataStruct.H:1460
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Immersed-forcing surface temperature flux [K m/s].
Definition: ERF_DataStruct.H:1459
SolverChoice::spatial_rhotheta_forcing
bool spatial_rhotheta_forcing
Whether spatially varying rho-theta forcing is enabled.
Definition: ERF_DataStruct.H:1497
SolverChoice::use_source_perturbation
bool use_source_perturbation(int lev) const
Query whether source-term turbulent perturbations are enabled on a level.
Definition: ERF_DataStruct.H:1547
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Whether custom vertical subsidence is enabled.
Definition: ERF_DataStruct.H:1492
SolverChoice::custom_moisture_forcing
bool custom_moisture_forcing
Whether custom moisture forcing is enabled.
Definition: ERF_DataStruct.H:1491
SolverChoice::if_surf_heating_rate
amrex::Real if_surf_heating_rate
Immersed-forcing surface heating rate [K/hr].
Definition: ERF_DataStruct.H:1461
SolverChoice::rad_type
RadiationType rad_type
Radiation model.
Definition: ERF_DataStruct.H:1608
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:1448
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Whether custom forcing operates on primitive variables.
Definition: ERF_DataStruct.H:1496
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Whether solution fields are nudged toward input sounding data.
Definition: ERF_DataStruct.H:1502
SolverChoice::use_num_diff
bool use_num_diff
Whether sixth-order numerical diffusion is enabled.
Definition: ERF_DataStruct.H:1600
SpongeChoice::sponge_type
static SpongeType sponge_type
Selected sponge damping model.
Definition: ERF_SpongeStruct.H:100
TurbChoice
Definition: ERF_TurbStruct.H:114
TurbChoice::diffuse_tke_3D
bool diffuse_tke_3D
Whether three-dimensional numerical diffusion is applied to TKE/QKE.
Definition: ERF_TurbStruct.H:731
TurbChoice::use_tke
bool use_tke
Whether any TKE or QKE closure is active.
Definition: ERF_TurbStruct.H:666
TurbulentPerturbation::pb_cell
amrex::Vector< amrex::MultiFab > pb_cell
Per-cell perturbation amplitude storage.
Definition: ERF_TurbPertStruct.H:755
TurbulentPerturbation::apply_tpi
void apply_tpi(const int &lev, const amrex::Box &vbx, const int &comp, const amrex::IndexType &m_ixtype, const amrex::Array4< amrex::Real > &src_arr, const amrex::Array4< amrex::Real const > &pert_cell)
Apply stored turbulent perturbations to a source or state array.
Definition: ERF_TurbPertStruct.H:401
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:105
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:124
Here is the call graph for this function: