ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_EBMOSTStress.H
Go to the documentation of this file.
1 /**
2  * \file ERF_EBMOSTStress.H
3  * \brief Defines EB Monin-Obukhov surface-layer flux functors.
4  */
5 #ifndef ERF_EBMOSTStress_H
6 #define ERF_EBMOSTStress_H
7 
8 #include <ERF_Constants.H>
9 #include <ERF_IndexDefines.H>
10 
11 
12 //! EB surface-layer model for adiabatic constant-roughness fluxes.
14 {
15  /**
16  * \brief Construct an adiabatic EB flux functor.
17  * \param Tflux Prescribed surface temperature flux.
18  * \param Qvflux Prescribed surface moisture flux.
19  */
21  amrex::Real Qvflux)
22  {
23  mdata.surf_temp_flux = Tflux;
24  mdata.surf_moist_flux = Qvflux;
25  }
26 
27  //! Populate neutral MOST flux scales for an adiabatic EB surface.
28  AMREX_GPU_DEVICE
29  AMREX_FORCE_INLINE
30  void
31  iterate_flux (const int& i,
32  const int& j,
33  const int& k,
34  const int& /*max_iters*/,
35  const amrex::Array4<const amrex::Real>& zref_arr,
36  const amrex::Array4<const amrex::Real>& z0_arr,
37  const amrex::Array4<const amrex::Real>& umm_arr,
38  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
39  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
40  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
41  const amrex::Array4<amrex::Real>& u_star_arr,
42  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
43  const amrex::Array4<amrex::Real>& t_star_arr,
44  const amrex::Array4<amrex::Real>& q_star_arr,
45  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
46  const amrex::Array4<amrex::Real>& /*q_surf_arr*/,
47  const amrex::Array4<amrex::Real>& olen_arr,
48  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
49  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
50  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
51  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
52  {
53  olen_arr(i,j,k) = bogus_large_value;
54  u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,0) / std::log(zref_arr(i,j,0) / z0_arr(i,j,k));
55  t_star_arr(i,j,k) = zero;
56  q_star_arr(i,j,k) = zero;
57  }
58 
59 private:
62 };
63 
64 
65 //! EB surface-layer model with prescribed surface temperature and constant roughness.
67 {
68  /**
69  * \brief Construct a prescribed-temperature EB flux functor.
70  * \param Tflux Prescribed surface temperature flux.
71  * \param Qvflux Prescribed surface moisture flux.
72  * \param cons_qflux Whether the moisture flux is prescribed directly.
73  */
75  amrex::Real Qvflux,
76  bool cons_qflux)
77  {
78  mdata.surf_temp_flux = Tflux;
79  mdata.surf_moist_flux = Qvflux;
80  spec_qflux = cons_qflux;
81  }
82 
83  //! Iterate MOST stability functions and update surface flux scales.
84  AMREX_GPU_DEVICE
85  AMREX_FORCE_INLINE
86  void
87  iterate_flux (const int& i,
88  const int& j,
89  const int& k,
90  const int& max_iters,
91  const amrex::Array4<const amrex::Real>& zref_arr,
92  const amrex::Array4<const amrex::Real>& z0_arr,
93  const amrex::Array4<const amrex::Real>& umm_arr,
94  const amrex::Array4<const amrex::Real>& tm_arr,
95  const amrex::Array4<const amrex::Real>& tvm_arr,
96  const amrex::Array4<const amrex::Real>& qvm_arr,
97  const amrex::Array4<amrex::Real>& u_star_arr,
98  const amrex::Array4<amrex::Real>& w_star_arr,
99  const amrex::Array4<amrex::Real>& t_star_arr,
100  const amrex::Array4<amrex::Real>& q_star_arr,
101  const amrex::Array4<amrex::Real>& t_surf_arr,
102  const amrex::Array4<amrex::Real>& q_surf_arr,
103  const amrex::Array4<amrex::Real>& olen_arr,
104  const amrex::Array4<amrex::Real>& pblh_arr,
105  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
106  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
107  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
108  {
109  amrex::Real Rib = zero;
110  amrex::Real zeta = zero;
111  amrex::Real zeta_old = zero;
112  amrex::Real psi_m = zero;
113  amrex::Real psi_h = zero;
114  amrex::Real num = zero;
115  amrex::Real den = zero;
116  amrex::Real zref = zref_arr(i,j,0);
117  amrex::Real z0 = z0_arr(i,j,k);
118  amrex::Real umm = std::max(umm_arr(i,j,0), WSMIN);
119  amrex::Real C = std::log(zref / z0);
120 
121  // First iteration we assume neutral (L -> inf)
122  if (u_star_arr(i,j,k) == bogus_large_value) {
123  olen_arr(i,j,k) = amrex::Real(1000);
124  }
125  zeta = zref / olen_arr(i,j,k);
126 
127  // Water vapor in atmos and surface
128  amrex::Real qv_s, qv_a;
129  if (q_surf_arr(i,j,k) > zero) {
130  qv_s = q_surf_arr(i,j,k);
131  } else {
132  // First iteration and no qv_surf was specified
133  // Use mean since there will be no flux
134  qv_s = qvm_arr(i,j,0);
135  }
136  qv_a = qvm_arr(i,j,0);
137 
138  // update w* and Umagmean from Beljaars (1995)
139  if (w_star_arr) {
140  // NOTE: Thv flux is lagged, similar to WRF
141  psi_m = sfuns.calc_psi_m2(zeta);
142  psi_h = sfuns.calc_psi_h2(zeta);
143  amrex::Real ustar = mdata.kappa * umm / (C - psi_m);
144  amrex::Real tstar = mdata.kappa * (tm_arr(i,j,0) - t_surf_arr(i,j,k)) / (C - psi_h);
146  -ustar * mdata.kappa * (qv_a - qv_s) / (C - psi_h);
147  amrex::Real tflux = -ustar*tstar*(one + epsv*qvm_arr(i,j,0)) + epsv*tm_arr(i,j,0)*qflux;
148  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,0));
149  amrex::Real wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
150  umm = std::sqrt(umm_arr(i,j,0)*umm_arr(i,j,0) + wstar*wstar);
151  umm = std::max(umm, WSMIN);
152  }
153 
154  // Bulk Richardson number w/ moisture
155  amrex::Real thv_s = t_surf_arr(i,j,k) * (one + epsv*qv_s);
156  amrex::Real thv_a = tm_arr(i,j,0) * (one + epsv*qv_a);
157  Rib = ( (mdata.gravity * zref) / tm_arr(i,j,0) ) *
158  ( (thv_a - thv_s) / (umm * umm) );
159  Rib = std::min(std::max(Rib,-amrex::Real(4.0)),amrex::Real(4.0));
160 
161  // Fixed point iteration on zeta
162  int iter = 0;
163  do {
164  // Transfer curr to old
165  zeta_old = zeta;
166 
167  // Stability functions
168  psi_m = sfuns.calc_psi_m2(zeta_old);
169  psi_h = sfuns.calc_psi_h2(zeta_old);
170 
171  // Limiting
172  num = std::max(C - psi_m, amrex::Real(1.0));
173  den = std::max(C - psi_h, amrex::Real(1.0));
174 
175  // Update with under relaxation
176  zeta = (one - alpha) * zeta_old + alpha * Rib * num * num / den;
177 
178  ++iter;
179  } while ( (std::abs(zeta - zeta_old) > tol) && (iter <= max_iters) );
180  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
181  "Maximum number of MOST iterations reached.");
182 
183  // Populate the stored MOST arrays
184  olen_arr(i,j,k) = zref / zeta;
185  u_star_arr(i,j,k) = mdata.kappa * umm / (C - psi_m);
186  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,0) - t_surf_arr(i,j,k)) / (C - psi_h);
187  if (spec_qflux) {
188  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (C - psi_h) /
189  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,0);
190  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
191  } else {
192  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,0) - q_surf_arr(i,j,k)) / (C - psi_h);
193  }
194  }
195 
196 private:
200  const amrex::Real tol = amrex::Real(1.0e-3);
202  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
203 };
204 
205 //! EB surface-layer model with prescribed surface fluxes and constant roughness.
207 {
208  /**
209  * \brief Construct a prescribed-flux EB functor.
210  * \param Tflux Prescribed surface temperature flux.
211  * \param Qvflux Prescribed surface moisture flux.
212  * \param cons_qflux Whether the moisture flux is prescribed directly.
213  */
215  amrex::Real Qvflux,
216  bool cons_qflux)
217  {
218  mdata.surf_temp_flux = Tflux;
219  mdata.surf_moist_flux = Qvflux;
220  spec_qflux = cons_qflux;
221  }
222 
223  //! Iterate MOST stability functions and update implied surface values.
224  AMREX_GPU_DEVICE
225  AMREX_FORCE_INLINE
226  void
227  iterate_flux (const int& i,
228  const int& j,
229  const int& k,
230  const int& max_iters,
231  const amrex::Array4<const amrex::Real>& zref_arr,
232  const amrex::Array4<const amrex::Real>& z0_arr,
233  const amrex::Array4<const amrex::Real>& umm_arr,
234  const amrex::Array4<const amrex::Real>& tm_arr,
235  const amrex::Array4<const amrex::Real>& tvm_arr,
236  const amrex::Array4<const amrex::Real>& qvm_arr,
237  const amrex::Array4<amrex::Real>& u_star_arr,
238  const amrex::Array4<amrex::Real>& w_star_arr,
239  const amrex::Array4<amrex::Real>& t_star_arr,
240  const amrex::Array4<amrex::Real>& q_star_arr,
241  const amrex::Array4<amrex::Real>& t_surf_arr,
242  const amrex::Array4<amrex::Real>& q_surf_arr,
243  const amrex::Array4<amrex::Real>& olen_arr,
244  const amrex::Array4<amrex::Real>& pblh_arr,
245  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
246  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
247  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
248  {
249  int iter = 0;
250  amrex::Real ustar = zero;
251  amrex::Real wstar = zero;
252  amrex::Real tflux = zero;
253  amrex::Real qflux = zero;
254  amrex::Real zeta = zero;
255  amrex::Real psi_m = zero;
256  amrex::Real psi_h = zero;
257  amrex::Real Olen = zero;
258  amrex::Real zref = zref_arr(i,j,0);
259  amrex::Real umm = std::max(umm_arr(i,j,0), WSMIN);
260  if (u_star_arr(i,j,k) == bogus_large_value) {
261  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
262  } else {
263  Olen = olen_arr(i,j,k);
264  zeta = zref / Olen;
265  psi_m = sfuns.calc_psi_m(zeta);
266  psi_h = sfuns.calc_psi_h(zeta);
267  }
268  do {
269  ustar = u_star_arr(i,j,k);
270  qflux = (spec_qflux) ? mdata.surf_moist_flux :
271  -(qvm_arr(i,j,0) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
272  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
273  tflux = mdata.surf_temp_flux*(one + epsv*qvm_arr(i,j,0)) + qflux*epsv*tm_arr(i,j,0);
274  if (w_star_arr) {
275  // update w* and Umagmean
276  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,0));
277  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
278  umm = std::sqrt(umm_arr(i,j,0)*umm_arr(i,j,0) + wstar*wstar);
279  umm = std::max(umm, WSMIN);
280  }
281  Olen = -ustar * ustar * ustar * tvm_arr(i,j,0) / (mdata.kappa * mdata.gravity * tflux);
282  zeta = zref / Olen;
283  psi_m = sfuns.calc_psi_m(zeta);
284  psi_h = sfuns.calc_psi_h(zeta);
285  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0_arr(i,j,k)) - psi_m);
286  ++iter;
287  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
288  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
289  "Maximum number of MOST iterations reached.");
290 
291  // Populate the stored MOST arrays
292  olen_arr(i,j,k) = Olen;
293  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
294  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,0);
295  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
296  if (spec_qflux) {
297  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
298  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,0);
299  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
300  } else {
301  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,0) - q_surf_arr(i,j,k)) /
302  (std::log(zref / z0_arr(i,j,k)) - psi_h);
303  }
304  }
305 
306 private:
310  const amrex::Real tol = amrex::Real(1.0e-5);
311  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
312 };
313 
314 //! EB implementation of the Moeng surface-flux formulation.
316 {
317  //! Construct a Moeng EB flux functor.
319 
320  //! Compute the EB moisture flux at a cell-centered cut-cell surface.
321  AMREX_GPU_DEVICE
322  AMREX_FORCE_INLINE
324  compute_q_flux (const int& i,
325  const int& j,
326  const int& k,
327  const amrex::Array4<const amrex::Real>& cons_arr,
328  const amrex::Array4<const amrex::Real>& velx_arr,
329  const amrex::Array4<const amrex::Real>& vely_arr,
330  const amrex::Array4<const amrex::Real>& umm_arr,
331  const amrex::Array4<const amrex::Real>& qvm_arr,
332  const amrex::Array4<const amrex::Real>& u_star_arr,
333  const amrex::Array4<const amrex::Real>& q_star_arr,
334  const amrex::Array4<const amrex::Real>& q_surf_arr,
335  const amrex::Array4<const amrex::Real>& u_vfrac_arr,
336  const amrex::Array4<const amrex::Real>& v_vfrac_arr) const
337  {
338  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
339  amrex::Real qv = cons_arr(i,j,k,RhoQ1_comp) / rho;
340 
341  // Volume-weighted average of x-face velocities to cell center
342  amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
343  amrex::Real velx = (u_vfrac_sum > eps) ?
344  (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
345  / u_vfrac_sum : zero;
346 
347  // Volume-weighted average of y-face velocities to cell center
348  amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
349  amrex::Real vely = (v_vfrac_sum > eps) ?
350  (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
351  / v_vfrac_sum : zero;
352 
353  amrex::Real qv_mean = qvm_arr(i,j,0);
354  amrex::Real ustar = u_star_arr(i,j,k);
355  amrex::Real qstar = q_star_arr(i,j,k);
356  amrex::Real qv_surf = q_surf_arr(i,j,k);
357  amrex::Real wsp_mean = umm_arr(i,j,0);
358  wsp_mean = std::max(wsp_mean, WSMIN);
359 
360  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
361  amrex::Real num1 = wsp * (qv_mean-qv_surf);
362  amrex::Real num2 = wsp_mean * (qv-qv_mean);
363 
364  // NOTE: this is rho*<Qv'w'> = -K dQvdz
365  amrex::Real moflux = (std::abs(qstar) > eps) ?
366  -rho*qstar*ustar*(num1+num2)/((qv_mean-qv_surf)*wsp_mean) : zero;
367 
368  return moflux;
369  }
370 
371  //! Compute the EB temperature flux at a cell-centered cut-cell surface.
372  AMREX_GPU_DEVICE
373  AMREX_FORCE_INLINE
375  compute_t_flux (const int& i,
376  const int& j,
377  const int& k,
378  const amrex::Array4<const amrex::Real>& cons_arr,
379  const amrex::Array4<const amrex::Real>& velx_arr,
380  const amrex::Array4<const amrex::Real>& vely_arr,
381  const amrex::Array4<const amrex::Real>& velz_arr,
382  const amrex::Array4<const amrex::Real>& umm_arr,
383  const amrex::Array4<const amrex::Real>& tm_arr,
384  const amrex::Array4<const amrex::Real>& u_star_arr,
385  const amrex::Array4<const amrex::Real>& t_star_arr,
386  const amrex::Array4<const amrex::Real>& t_surf_arr,
387  const amrex::Array4<const amrex::Real>& u_vfrac_arr,
388  const amrex::Array4<const amrex::Real>& v_vfrac_arr,
389  const amrex::Array4<const amrex::Real>& w_vfrac_arr,
390  const amrex::Array4<const amrex::Real>& bnorm_arr) const
391  {
392  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
393  amrex::Real theta = cons_arr(i,j,k,RhoTheta_comp) / rho;
394 
395  // Volume-weighted average of x-face velocities to cell center
396  amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
397  amrex::Real velx = (u_vfrac_sum > eps) ?
398  (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
399  / u_vfrac_sum : zero;
400 
401  // Volume-weighted average of y-face velocities to cell center
402  amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
403  amrex::Real vely = (v_vfrac_sum > eps) ?
404  (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
405  / v_vfrac_sum : zero;
406 
407  // Volume-weighted average of z-face velocities to cell center
408  amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1);
409  amrex::Real velz = (w_vfrac_sum > eps) ?
410  (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1))
411  / w_vfrac_sum : zero;
412 
413  // Get boundary normal components
414  amrex::Real nx = bnorm_arr(i,j,k,0);
415  amrex::Real ny = bnorm_arr(i,j,k,1);
416  amrex::Real nz = bnorm_arr(i,j,k,2);
417 
418  // Project velocity onto tangent plane (remove normal component)
419  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
420  amrex::Real velx_tangent = velx - v_dot_n * nx;
421  amrex::Real vely_tangent = vely - v_dot_n * ny;
422 
423  amrex::Real theta_mean = tm_arr(i,j,0);
424  amrex::Real ustar = u_star_arr(i,j,k);
425  amrex::Real tstar = t_star_arr(i,j,k);
426  amrex::Real theta_surf = t_surf_arr(i,j,k);
427  amrex::Real wsp_mean = umm_arr(i,j,0);
428  wsp_mean = std::max(wsp_mean, WSMIN);
429 
430  // Use tangential velocity magnitude instead of Cartesian
431  amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent);
432  amrex::Real num1 = wsp * (theta_mean-theta_surf);
433  amrex::Real num2 = wsp_mean * (theta-theta_mean);
434 
435  // NOTE: this is rho*<T'w'> = -K dTdz
436  amrex::Real moflux = (std::abs(tstar) > eps) ?
437  -rho*tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : zero;
438 
439  return moflux;
440  }
441 
442  //! Compute the EB x-momentum stress on a staggered face.
443  AMREX_GPU_DEVICE
444  AMREX_FORCE_INLINE
447  int j,
448  int k,
449  const amrex::Array4<const amrex::Real>& cons_arr,
450  const amrex::Array4<const amrex::Real>& velx_arr,
451  const amrex::Array4<const amrex::Real>& vely_arr,
452  const amrex::Array4<const amrex::Real>& velz_arr,
453  const amrex::Array4<const amrex::Real>& umm_arr,
454  const amrex::Array4<const amrex::Real>& um_arr,
455  const amrex::Array4<const amrex::Real>& u_star_arr,
456  const amrex::Array4<const amrex::Real>& u_vfrac_arr,
457  const amrex::Array4<const amrex::Real>& v_vfrac_arr,
458  const amrex::Array4<const amrex::Real>& w_vfrac_arr,
459  const amrex::Array4<const amrex::Real>& cc_vfrac_arr,
460  const amrex::Array4<const amrex::EBCellFlag>& cc_flag_arr,
461  const amrex::Array4<const amrex::Real>& bnorm_arr,
462  int idir = 0) const
463  {
464  amrex::Real velx, vely, rho, ustar, wsp_mean;
465  amrex::Real velx_tangent, vely_tangent;
466 
467  if (idir == 0) {
468  // x-face: average to x-face
469  velx = velx_arr(i,j,k);
470 
471  // Volume-weighted average of y-face velocities to x-face
472  amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k) +
473  v_vfrac_arr(i-1,j,k) + v_vfrac_arr(i-1,j+1,k);
474  vely = (v_vfrac_sum > eps) ?
475  (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k) +
476  vely_arr(i-1,j,k) * v_vfrac_arr(i-1,j,k) + vely_arr(i-1,j+1,k) * v_vfrac_arr(i-1,j+1,k))
477  / v_vfrac_sum : zero;
478 
479  // Volume-weighted average of z-face velocities to x-face
480  amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
481  w_vfrac_arr(i-1,j,k) + w_vfrac_arr(i-1,j,k+1);
482  amrex::Real velz = (w_vfrac_sum > eps) ?
483  (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
484  velz_arr(i-1,j,k) * w_vfrac_arr(i-1,j,k) + velz_arr(i-1,j,k+1) * w_vfrac_arr(i-1,j,k+1))
485  / w_vfrac_sum : zero;
486 
487  // Get boundary normal at x-face (already at correct staggered location)
488  amrex::Real nx = bnorm_arr(i,j,k,0);
489  amrex::Real ny = bnorm_arr(i,j,k,1);
490  amrex::Real nz = bnorm_arr(i,j,k,2);
491 
492  // Project velocity onto tangent plane
493  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
494  velx_tangent = velx - v_dot_n * nx;
495  vely_tangent = vely - v_dot_n * ny;
496 
497  // Volume-weighted average of cell-centered density to x-face
498  amrex::Real cc_vfrac_sum = cc_vfrac_arr(i-1,j,k) + cc_vfrac_arr(i,j,k);
499  rho = (cc_vfrac_sum > eps) ?
500  (cons_arr(i-1,j,k,Rho_comp) * cc_vfrac_arr(i-1,j,k) + cons_arr(i,j,k,Rho_comp) * cc_vfrac_arr(i,j,k))
501  / cc_vfrac_sum : zero;
502 
503  // Average cell-centered u_star and wsp_mean to x-face, using only valid (SingleValued) cells
504  bool low_valid = cc_flag_arr(i-1,j,k).isSingleValued();
505  bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
506 
507  if (low_valid && high_valid) {
508  ustar = myhalf * (u_star_arr(i-1,j,k) + u_star_arr(i,j,k));
509  wsp_mean = myhalf * (umm_arr(i-1,j,0) + umm_arr(i,j,0));
510  } else if (low_valid) {
511  ustar = u_star_arr(i-1,j,k);
512  wsp_mean = umm_arr(i-1,j,0);
513  } else if (high_valid) {
514  ustar = u_star_arr(i,j,k);
515  wsp_mean = umm_arr(i,j,0);
516  } else {
517  ustar = zero;
518  wsp_mean = WSMIN;
519  }
520 
521  } else if (idir == 1) {
522  // y-face: average to y-face
523  vely = vely_arr(i,j,k);
524 
525  // Volume-weighted average of x-face velocities to y-face
526  amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k) +
527  u_vfrac_arr(i,j-1,k) + u_vfrac_arr(i+1,j-1,k);
528  velx = (u_vfrac_sum > eps) ?
529  (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k) +
530  velx_arr(i,j-1,k) * u_vfrac_arr(i,j-1,k) + velx_arr(i+1,j-1,k) * u_vfrac_arr(i+1,j-1,k))
531  / u_vfrac_sum : zero;
532 
533  // Volume-weighted average of z-face velocities to y-face
534  amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
535  w_vfrac_arr(i,j-1,k) + w_vfrac_arr(i,j-1,k+1);
536  amrex::Real velz = (w_vfrac_sum > eps) ?
537  (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
538  velz_arr(i,j-1,k) * w_vfrac_arr(i,j-1,k) + velz_arr(i,j-1,k+1) * w_vfrac_arr(i,j-1,k+1))
539  / w_vfrac_sum : zero;
540 
541  // Get boundary normal at y-face (already at correct staggered location)
542  amrex::Real nx = bnorm_arr(i,j,k,0);
543  amrex::Real ny = bnorm_arr(i,j,k,1);
544  amrex::Real nz = bnorm_arr(i,j,k,2);
545 
546  // Project velocity onto tangent plane
547  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
548  velx_tangent = velx - v_dot_n * nx;
549  vely_tangent = vely - v_dot_n * ny;
550 
551  // Volume-weighted average of cell-centered density to y-face
552  amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j-1,k) + cc_vfrac_arr(i,j,k);
553  rho = (cc_vfrac_sum > eps) ?
554  (cons_arr(i,j-1,k,Rho_comp) * cc_vfrac_arr(i,j-1,k) + cons_arr(i,j,k,Rho_comp) * cc_vfrac_arr(i,j,k))
555  / cc_vfrac_sum : zero;
556 
557  // Average cell-centered u_star and wsp_mean to y-face, using only valid (SingleValued) cells
558  bool low_valid = cc_flag_arr(i,j-1,k).isSingleValued();
559  bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
560 
561  if (low_valid && high_valid) {
562  ustar = myhalf * (u_star_arr(i,j-1,k) + u_star_arr(i,j,k));
563  wsp_mean = myhalf * (umm_arr(i,j-1,0) + umm_arr(i,j,0));
564  } else if (low_valid) {
565  ustar = u_star_arr(i,j-1,k);
566  wsp_mean = umm_arr(i,j-1,0);
567  } else if (high_valid) {
568  ustar = u_star_arr(i,j,k);
569  wsp_mean = umm_arr(i,j,0);
570  } else {
571  ustar = zero;
572  wsp_mean = WSMIN;
573  }
574 
575  } else {
576  // z-face: average to z-face
577  // Volume-weighted average of x-face velocities to z-face
578  amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k-1) + u_vfrac_arr(i+1,j,k-1) +
579  u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
580  velx = (u_vfrac_sum > eps) ?
581  (velx_arr(i,j,k-1) * u_vfrac_arr(i,j,k-1) + velx_arr(i+1,j,k-1) * u_vfrac_arr(i+1,j,k-1) +
582  velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
583  / u_vfrac_sum : zero;
584 
585  // Volume-weighted average of y-face velocities to z-face
586  amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k-1) + v_vfrac_arr(i,j+1,k-1) +
587  v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
588  vely = (v_vfrac_sum > eps) ?
589  (vely_arr(i,j,k-1) * v_vfrac_arr(i,j,k-1) + vely_arr(i,j+1,k-1) * v_vfrac_arr(i,j+1,k-1) +
590  vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
591  / v_vfrac_sum : zero;
592 
593  // z-velocity is already at z-face
594  amrex::Real velz = velz_arr(i,j,k);
595 
596  // Get boundary normal at z-face (already at correct staggered location)
597  amrex::Real nx = bnorm_arr(i,j,k,0);
598  amrex::Real ny = bnorm_arr(i,j,k,1);
599  amrex::Real nz = bnorm_arr(i,j,k,2);
600 
601  // Project velocity onto tangent plane
602  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
603  velx_tangent = velx - v_dot_n * nx;
604  vely_tangent = vely - v_dot_n * ny;
605 
606  // Volume-weighted average of cell-centered density to z-face
607  amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j,k-1) + cc_vfrac_arr(i,j,k);
608  rho = (cc_vfrac_sum > eps) ?
609  (cons_arr(i,j,k-1,Rho_comp) * cc_vfrac_arr(i,j,k-1) + cons_arr(i,j,k,Rho_comp) * cc_vfrac_arr(i,j,k))
610  / cc_vfrac_sum : zero;
611 
612  // Average cell-centered u_star and wsp_mean to z-face, using only valid (SingleValued) cells
613  bool low_valid = cc_flag_arr(i,j,k-1).isSingleValued();
614  bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
615 
616  if (low_valid && high_valid) {
617  ustar = myhalf * (u_star_arr(i,j,k-1) + u_star_arr(i,j,k));
618  wsp_mean = myhalf * (umm_arr(i,j,0) + umm_arr(i,j,0));
619  } else if (low_valid) {
620  ustar = u_star_arr(i,j,k-1);
621  wsp_mean = umm_arr(i,j,0);
622  } else if (high_valid) {
623  ustar = u_star_arr(i,j,k);
624  wsp_mean = umm_arr(i,j,0);
625  } else {
626  ustar = zero;
627  wsp_mean = WSMIN;
628  }
629  }
630 
631  wsp_mean = std::max(wsp_mean, WSMIN);
632  amrex::Real umean = um_arr(i,j,0);
633 
634  // Note: The surface mean shear stress is decomposed into tau_xz by
635  // multiplying the modeled shear stress (rho*ustar^2) with
636  // a factor of umean/wsp_mean for directionality; this factor
637  // modifies the denominator from what is in Moeng amrex::Real(1984.)
638  amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent);
639  amrex::Real num1 = wsp * umean;
640  amrex::Real num2 = wsp_mean * (velx_tangent-umean);
641 
642  // NOTE: this is rho*<u'w'> = -K dudz
643  amrex::Real stressx = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
644 
645  return stressx;
646  }
647 
648  //! Compute the EB y-momentum stress on a staggered face.
649  AMREX_GPU_DEVICE
650  AMREX_FORCE_INLINE
653  int j,
654  int k,
655  const amrex::Array4<const amrex::Real>& cons_arr,
656  const amrex::Array4<const amrex::Real>& velx_arr,
657  const amrex::Array4<const amrex::Real>& vely_arr,
658  const amrex::Array4<const amrex::Real>& velz_arr,
659  const amrex::Array4<const amrex::Real>& umm_arr,
660  const amrex::Array4<const amrex::Real>& vm_arr,
661  const amrex::Array4<const amrex::Real>& u_star_arr,
662  const amrex::Array4<const amrex::Real>& u_vfrac_arr,
663  const amrex::Array4<const amrex::Real>& v_vfrac_arr,
664  const amrex::Array4<const amrex::Real>& w_vfrac_arr,
665  const amrex::Array4<const amrex::Real>& cc_vfrac_arr,
666  const amrex::Array4<const amrex::EBCellFlag>& cc_flag_arr,
667  const amrex::Array4<const amrex::Real>& bnorm_arr,
668  int idir = 0) const
669  {
670  amrex::Real velx, vely, rho, ustar, wsp_mean;
671  amrex::Real velx_tangent, vely_tangent;
672 
673  if (idir == 0) {
674  // x-face: average from cells (i-1) and (i)
675  // Volume-weighted average of y-face velocities to x-face
676  amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k) +
677  v_vfrac_arr(i-1,j,k) + v_vfrac_arr(i-1,j+1,k);
678  vely = (v_vfrac_sum > eps) ?
679  (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k) +
680  vely_arr(i-1,j,k) * v_vfrac_arr(i-1,j,k) + vely_arr(i-1,j+1,k) * v_vfrac_arr(i-1,j+1,k))
681  / v_vfrac_sum : zero;
682 
683  velx = velx_arr(i,j,k);
684 
685  // Volume-weighted average of z-face velocities to x-face
686  amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
687  w_vfrac_arr(i-1,j,k) + w_vfrac_arr(i-1,j,k+1);
688  amrex::Real velz = (w_vfrac_sum > eps) ?
689  (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
690  velz_arr(i-1,j,k) * w_vfrac_arr(i-1,j,k) + velz_arr(i-1,j,k+1) * w_vfrac_arr(i-1,j,k+1))
691  / w_vfrac_sum : zero;
692 
693  // Get boundary normal at x-face (already at correct staggered location)
694  amrex::Real nx = bnorm_arr(i,j,k,0);
695  amrex::Real ny = bnorm_arr(i,j,k,1);
696  amrex::Real nz = bnorm_arr(i,j,k,2);
697 
698  // Project velocity onto tangent plane
699  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
700  velx_tangent = velx - v_dot_n * nx;
701  vely_tangent = vely - v_dot_n * ny;
702 
703  // Volume-weighted average of cell-centered density to x-face
704  amrex::Real cc_vfrac_sum = cc_vfrac_arr(i-1,j,k) + cc_vfrac_arr(i,j,k);
705  rho = (cc_vfrac_sum > eps) ?
706  (cons_arr(i-1,j,k,Rho_comp) * cc_vfrac_arr(i-1,j,k) + cons_arr(i,j,k,Rho_comp) * cc_vfrac_arr(i,j,k))
707  / cc_vfrac_sum : zero;
708 
709  // Average cell-centered u_star and wsp_mean to x-face, using only valid (SingleValued) cells
710  bool low_valid = cc_flag_arr(i-1,j,k).isSingleValued();
711  bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
712 
713  if (low_valid && high_valid) {
714  ustar = myhalf * (u_star_arr(i-1,j,k) + u_star_arr(i,j,k));
715  wsp_mean = myhalf * (umm_arr(i-1,j,0) + umm_arr(i,j,0));
716  } else if (low_valid) {
717  ustar = u_star_arr(i-1,j,k);
718  wsp_mean = umm_arr(i-1,j,0);
719  } else if (high_valid) {
720  ustar = u_star_arr(i,j,k);
721  wsp_mean = umm_arr(i,j,0);
722  } else {
723  ustar = zero;
724  wsp_mean = WSMIN;
725  }
726 
727  } else if (idir == 1) {
728  // y-face: average from cells (j-1) and (j)
729  // Volume-weighted average of x-face velocities to y-face
730  amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k) +
731  u_vfrac_arr(i,j-1,k) + u_vfrac_arr(i+1,j-1,k);
732  velx = (u_vfrac_sum > eps) ?
733  (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k) +
734  velx_arr(i,j-1,k) * u_vfrac_arr(i,j-1,k) + velx_arr(i+1,j-1,k) * u_vfrac_arr(i+1,j-1,k))
735  / u_vfrac_sum : zero;
736 
737  vely = vely_arr(i,j,k);
738 
739  // Volume-weighted average of z-face velocities to y-face
740  amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
741  w_vfrac_arr(i,j-1,k) + w_vfrac_arr(i,j-1,k+1);
742  amrex::Real velz = (w_vfrac_sum > eps) ?
743  (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
744  velz_arr(i,j-1,k) * w_vfrac_arr(i,j-1,k) + velz_arr(i,j-1,k+1) * w_vfrac_arr(i,j-1,k+1))
745  / w_vfrac_sum : zero;
746 
747  // Get boundary normal at y-face (already at correct staggered location)
748  amrex::Real nx = bnorm_arr(i,j,k,0);
749  amrex::Real ny = bnorm_arr(i,j,k,1);
750  amrex::Real nz = bnorm_arr(i,j,k,2);
751 
752  // Project velocity onto tangent plane
753  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
754  velx_tangent = velx - v_dot_n * nx;
755  vely_tangent = vely - v_dot_n * ny;
756 
757  // Volume-weighted average of cell-centered density to y-face
758  amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j-1,k) + cc_vfrac_arr(i,j,k);
759  rho = (cc_vfrac_sum > eps) ?
760  (cons_arr(i,j-1,k,Rho_comp) * cc_vfrac_arr(i,j-1,k) + cons_arr(i,j,k,Rho_comp) * cc_vfrac_arr(i,j,k))
761  / cc_vfrac_sum : zero;
762 
763  // Average cell-centered u_star and wsp_mean to y-face, using only valid (SingleValued) cells
764  bool low_valid = cc_flag_arr(i,j-1,k).isSingleValued();
765  bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
766 
767  if (low_valid && high_valid) {
768  ustar = myhalf * (u_star_arr(i,j-1,k) + u_star_arr(i,j,k));
769  wsp_mean = myhalf * (umm_arr(i,j-1,0) + umm_arr(i,j,0));
770  } else if (low_valid) {
771  ustar = u_star_arr(i,j-1,k);
772  wsp_mean = umm_arr(i,j-1,0);
773  } else if (high_valid) {
774  ustar = u_star_arr(i,j,k);
775  wsp_mean = umm_arr(i,j,0);
776  } else {
777  ustar = zero;
778  wsp_mean = WSMIN;
779  }
780 
781  } else {
782  // z-face: average from cells (k-1) and (k)
783  // Volume-weighted average of x-face velocities to z-face
784  amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k-1) + u_vfrac_arr(i+1,j,k-1) +
785  u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
786  velx = (u_vfrac_sum > eps) ?
787  (velx_arr(i,j,k-1) * u_vfrac_arr(i,j,k-1) + velx_arr(i+1,j,k-1) * u_vfrac_arr(i+1,j,k-1) +
788  velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
789  / u_vfrac_sum : zero;
790 
791  // Volume-weighted average of y-face velocities to z-face
792  amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k-1) + v_vfrac_arr(i,j+1,k-1) +
793  v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
794  vely = (v_vfrac_sum > eps) ?
795  (vely_arr(i,j,k-1) * v_vfrac_arr(i,j,k-1) + vely_arr(i,j+1,k-1) * v_vfrac_arr(i,j+1,k-1) +
796  vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
797  / v_vfrac_sum : zero;
798 
799  // z-velocity is already at z-face
800  amrex::Real velz = velz_arr(i,j,k);
801 
802  // Get boundary normal at z-face (already at correct staggered location)
803  amrex::Real nx = bnorm_arr(i,j,k,0);
804  amrex::Real ny = bnorm_arr(i,j,k,1);
805  amrex::Real nz = bnorm_arr(i,j,k,2);
806 
807  // Project velocity onto tangent plane
808  amrex::Real v_dot_n = velx*nx + vely*ny + velz*nz;
809  velx_tangent = velx - v_dot_n * nx;
810  vely_tangent = vely - v_dot_n * ny;
811 
812  // Volume-weighted average of cell-centered density to z-face
813  amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j,k-1) + cc_vfrac_arr(i,j,k);
814  rho = (cc_vfrac_sum > eps) ?
815  (cons_arr(i,j,k-1,Rho_comp) * cc_vfrac_arr(i,j,k-1) + cons_arr(i,j,k,Rho_comp) * cc_vfrac_arr(i,j,k))
816  / cc_vfrac_sum : zero;
817 
818  // Average cell-centered u_star and wsp_mean to z-face, using only valid (SingleValued) cells
819  bool low_valid = cc_flag_arr(i,j,k-1).isSingleValued();
820  bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
821 
822  if (low_valid && high_valid) {
823  ustar = myhalf * (u_star_arr(i,j,k-1) + u_star_arr(i,j,k));
824  wsp_mean = myhalf * (umm_arr(i,j,0) + umm_arr(i,j,0));
825  } else if (low_valid) {
826  ustar = u_star_arr(i,j,k-1);
827  wsp_mean = umm_arr(i,j,0);
828  } else if (high_valid) {
829  ustar = u_star_arr(i,j,k);
830  wsp_mean = umm_arr(i,j,0);
831  } else {
832  ustar = zero;
833  wsp_mean = WSMIN;
834  }
835  }
836 
837  wsp_mean = std::max(wsp_mean, WSMIN);
838  amrex::Real vmean = vm_arr(i,j,0);
839 
840  // Note: The surface mean shear stress is decomposed into tau_yz by
841  // multiplying the modeled shear stress (rho*ustar^2) with
842  // a factor of vmean/wsp_mean for directionality; this factor
843  // modifies the denominator from what is in Moeng amrex::Real(1984.)
844  amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent);
845  amrex::Real num1 = wsp * vmean;
846  amrex::Real num2 = wsp_mean * (vely_tangent-vmean);
847 
848  // NOTE: this is rho*<v'w'> = -K dvdz
849  amrex::Real stressy = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
850 
851  return stressy;
852  }
853 
854 private:
855 #ifdef AMREX_USE_FLOAT
856  const amrex::Real eps = amrex::Real(1e-6);
857 #else
858  const amrex::Real eps = amrex::Real(1e-12);
859 #endif
860  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
861 };
862 #endif
constexpr amrex::Real epsv
Definition: ERF_Constants.H:53
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
@ num
Definition: ERF_DataStruct.H:27
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
rho
Definition: ERF_InitCustomPert_Bubble.H:107
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real calc_wstar(const amrex::Real &ust, const amrex::Real &tst, const amrex::Real &qst, const amrex::Real &pblh, const amrex::Real &th, const amrex::Real &thv, const amrex::Real &qv=amrex::Real(0))
Definition: ERF_Wstar.H:13
@ theta
Definition: ERF_SLM.H:20
@ qv
Definition: ERF_Kessler.H:30
@ den
Definition: ERF_AdvanceWSM6.cpp:109
EB surface-layer model for adiabatic constant-roughness fluxes.
Definition: ERF_EBMOSTStress.H:14
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< const amrex::Real > &z0_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &t_star_arr, const amrex::Array4< amrex::Real > &q_star_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &) const
Populate neutral MOST flux scales for an adiabatic EB surface.
Definition: ERF_EBMOSTStress.H:31
adiabatic_eb(amrex::Real Tflux, amrex::Real Qvflux)
Construct an adiabatic EB flux functor.
Definition: ERF_EBMOSTStress.H:20
most_data mdata
Definition: ERF_EBMOSTStress.H:60
similarity_funs sfuns
Definition: ERF_EBMOSTStress.H:61
EB implementation of the Moeng surface-flux formulation.
Definition: ERF_EBMOSTStress.H:316
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_u_flux(int i, int j, int k, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &velx_arr, const amrex::Array4< const amrex::Real > &vely_arr, const amrex::Array4< const amrex::Real > &velz_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &um_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr, const amrex::Array4< const amrex::Real > &w_vfrac_arr, const amrex::Array4< const amrex::Real > &cc_vfrac_arr, const amrex::Array4< const amrex::EBCellFlag > &cc_flag_arr, const amrex::Array4< const amrex::Real > &bnorm_arr, int idir=0) const
Compute the EB x-momentum stress on a staggered face.
Definition: ERF_EBMOSTStress.H:446
const amrex::Real WSMIN
Definition: ERF_EBMOSTStress.H:860
moeng_flux_eb()
Construct a Moeng EB flux functor.
Definition: ERF_EBMOSTStress.H:318
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_t_flux(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &velx_arr, const amrex::Array4< const amrex::Real > &vely_arr, const amrex::Array4< const amrex::Real > &velz_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< const amrex::Real > &t_star_arr, const amrex::Array4< const amrex::Real > &t_surf_arr, const amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr, const amrex::Array4< const amrex::Real > &w_vfrac_arr, const amrex::Array4< const amrex::Real > &bnorm_arr) const
Compute the EB temperature flux at a cell-centered cut-cell surface.
Definition: ERF_EBMOSTStress.H:375
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_v_flux(int i, int j, int k, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &velx_arr, const amrex::Array4< const amrex::Real > &vely_arr, const amrex::Array4< const amrex::Real > &velz_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &vm_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr, const amrex::Array4< const amrex::Real > &w_vfrac_arr, const amrex::Array4< const amrex::Real > &cc_vfrac_arr, const amrex::Array4< const amrex::EBCellFlag > &cc_flag_arr, const amrex::Array4< const amrex::Real > &bnorm_arr, int idir=0) const
Compute the EB y-momentum stress on a staggered face.
Definition: ERF_EBMOSTStress.H:652
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_q_flux(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &velx_arr, const amrex::Array4< const amrex::Real > &vely_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &qvm_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< const amrex::Real > &q_star_arr, const amrex::Array4< const amrex::Real > &q_surf_arr, const amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr) const
Compute the EB moisture flux at a cell-centered cut-cell surface.
Definition: ERF_EBMOSTStress.H:324
const amrex::Real eps
Definition: ERF_EBMOSTStress.H:858
Definition: ERF_MOSTStress.H:13
amrex::Real surf_moist_flux
Moisture flux.
Definition: ERF_MOSTStress.H:19
amrex::Real kappa
von Karman constant
Definition: ERF_MOSTStress.H:16
amrex::Real gravity
Acceleration due to gravity (m/s^2)
Definition: ERF_MOSTStress.H:17
const amrex::Real Bjr_beta
Definition: ERF_MOSTStress.H:32
amrex::Real surf_temp_flux
Heat flux TODO: decide whether this is <θ'w'> or <θv'w'> under moist conditions.
Definition: ERF_MOSTStress.H:18
Definition: ERF_MOSTStress.H:40
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:105
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h2(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:77
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:124
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m2(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:52
EB surface-layer model with prescribed surface fluxes and constant roughness.
Definition: ERF_EBMOSTStress.H:207
const amrex::Real WSMIN
Definition: ERF_EBMOSTStress.H:311
similarity_funs sfuns
Definition: ERF_EBMOSTStress.H:309
const amrex::Real tol
Definition: ERF_EBMOSTStress.H:310
most_data mdata
Definition: ERF_EBMOSTStress.H:307
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< const amrex::Real > &z0_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &tvm_arr, const amrex::Array4< const amrex::Real > &qvm_arr, const amrex::Array4< amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &w_star_arr, const amrex::Array4< amrex::Real > &t_star_arr, const amrex::Array4< amrex::Real > &q_star_arr, const amrex::Array4< amrex::Real > &t_surf_arr, const amrex::Array4< amrex::Real > &q_surf_arr, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &pblh_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &) const
Iterate MOST stability functions and update implied surface values.
Definition: ERF_EBMOSTStress.H:227
bool spec_qflux
Definition: ERF_EBMOSTStress.H:308
surface_flux_eb(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Construct a prescribed-flux EB functor.
Definition: ERF_EBMOSTStress.H:214
EB surface-layer model with prescribed surface temperature and constant roughness.
Definition: ERF_EBMOSTStress.H:67
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< const amrex::Real > &z0_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &tvm_arr, const amrex::Array4< const amrex::Real > &qvm_arr, const amrex::Array4< amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &w_star_arr, const amrex::Array4< amrex::Real > &t_star_arr, const amrex::Array4< amrex::Real > &q_star_arr, const amrex::Array4< amrex::Real > &t_surf_arr, const amrex::Array4< amrex::Real > &q_surf_arr, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &pblh_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &) const
Iterate MOST stability functions and update surface flux scales.
Definition: ERF_EBMOSTStress.H:87
const amrex::Real tol
Definition: ERF_EBMOSTStress.H:200
const amrex::Real alpha
Definition: ERF_EBMOSTStress.H:201
const amrex::Real WSMIN
Definition: ERF_EBMOSTStress.H:202
similarity_funs sfuns
Definition: ERF_EBMOSTStress.H:199
surface_temp_eb(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Construct a prescribed-temperature EB flux functor.
Definition: ERF_EBMOSTStress.H:74
bool spec_qflux
Definition: ERF_EBMOSTStress.H:198
most_data mdata
Definition: ERF_EBMOSTStress.H:197