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