ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MOSTStress.H
Go to the documentation of this file.
1 #ifndef ERF_MOSTStress_H
2 #define ERF_MOSTStress_H
3 
4 #include <ERF_Constants.H>
5 #include <ERF_IndexDefines.H>
6 #include <ERF_MOSTRoughness.H>
7 #include <ERF_Wstar.H>
8 
9 /**
10  * Structure of plain old data relevant to MOST BCs
11  */
12 struct most_data
13 {
14 public:
15  amrex::Real zref{10.0}; ///< Reference height (m)
16  amrex::Real z0_const{0.1}; ///< Roughness height -- default constant value(m)
17  amrex::Real kappa{KAPPA}; ///< von Karman constant
18  amrex::Real gravity{CONST_GRAV}; ///< Acceleration due to gravity (m/s^2)
19  amrex::Real surf_temp_flux{0.0}; ///< Heat flux
20 
21  amrex::Real Cnk_a{0.0185}; ///< Standard Charnock constant https://doi.org/10.1175/JAMC-D-17-0137.1
22  amrex::Real Cnk_b1{1.0/30.0}; ///< Modified Charnock Eq (4) https://doi.org/10.1175/JAMC-D-17-0137.1
23  amrex::Real Cnk_b2{1260.0}; ///< Modified Charnock Eq (4) https://doi.org/10.1175/JAMC-D-17-0137.1
24  amrex::Real Cnk_d{30.0}; ///< Modified Charnock Eq (4) https://doi.org/10.1175/JAMC-D-17-0137.1
25  amrex::Real Cnk_b;
26  bool visc{false}; ///< Use viscous Charnock formulation
27 
28  const amrex::Real Bjr_beta = 1.2; // Empirical parameter from Beljaars 1995 QJRMS
29 };
30 
31 
32 /**
33  * Structure of similarity functions for Moeng formulation
34  */
36 {
37  AMREX_GPU_HOST_DEVICE
38  AMREX_FORCE_INLINE
39  amrex::Real
40  calc_psi_m (amrex::Real zeta) const
41  {
42  if (zeta > 0) {
43  return -beta_m * zeta;
44  } else {
45  amrex::Real x = std::sqrt(std::sqrt(1.0 - gamma_m * zeta));
46  return 2.0 * std::log(0.5 * (1.0 + x)) + log(0.5 * (1.0 + x * x)) -
47  2.0 * std::atan(x) + PIoTwo;
48  }
49  }
50 
51  AMREX_GPU_HOST_DEVICE
52  AMREX_FORCE_INLINE
53  amrex::Real
54  calc_psi_h (amrex::Real zeta) const
55  {
56  if (zeta > 0) {
57  return -beta_h * zeta;
58  } else {
59  amrex::Real x = std::sqrt(1.0 - gamma_h * zeta);
60  return 2.0 * std::log(0.5 * (1.0 + x));
61  }
62  }
63 
64 private:
65  amrex::Real beta_m{5.0}; ///< Constants from Dyer, BLM, 1974
66  amrex::Real beta_h{5.0}; ///< https://doi.org/10.1007/BF00240838
67  amrex::Real gamma_m{16.0};
68  amrex::Real gamma_h{16.0};
69 };
70 
71 
72 /**
73  * Empirical kinematic viscosity [m2/s] formula from Andreas (1989) CRREL Rep.
74  * 89-11, valid between -173 and 277 deg C.
75  */
76 AMREX_GPU_DEVICE
77 AMREX_FORCE_INLINE
78 amrex::Real
79 air_viscosity (amrex::Real T_degK)
80 {
81  amrex::Real TC = T_degK - 273.15;
82  return 1.326e-5*(1. + 6.542e-3*TC + 8.301e-6*TC*TC - 4.84e-9*TC*TC*TC);
83 }
84 
85 
86 /**
87  * Adiabatic with constant roughness
88  */
89 struct adiabatic
90 {
91  adiabatic (amrex::Real zref,
92  amrex::Real flux)
93  {
94  mdata.zref = zref;
95  mdata.surf_temp_flux = flux;
96  }
97 
98  AMREX_GPU_DEVICE
99  AMREX_FORCE_INLINE
100  void
101  iterate_flux (const int& i,
102  const int& j,
103  const int& k,
104  const int& /*max_iters*/,
105  const amrex::Array4<const amrex::Real>& z0_arr,
106  const amrex::Array4<const amrex::Real>& umm_arr,
107  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
108  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
109  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
110  const amrex::Array4<amrex::Real>& u_star_arr,
111  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
112  const amrex::Array4<amrex::Real>& t_star_arr,
113  const amrex::Array4<amrex::Real>& /*q_star_arr*/,
114  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
115  const amrex::Array4<amrex::Real>& olen_arr,
116  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
117  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
118  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
119  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
120  {
121  u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k));
122  t_star_arr(i,j,k) = 0.0;
123  olen_arr(i,j,k) = 1.0e16;
124  }
125 
126 private:
129 };
130 
131 
132 /**
133  * Adiabatic with charnock roughness
134  */
136 {
137  adiabatic_charnock (amrex::Real zref,
138  amrex::Real flux,
139  amrex::Real cnk_a,
140  bool cnk_visc)
141  {
142  mdata.zref = zref;
143  mdata.surf_temp_flux = flux;
144  mdata.Cnk_a = cnk_a;
145  mdata.visc = cnk_visc;
146  }
147 
148  AMREX_GPU_DEVICE
149  AMREX_FORCE_INLINE
150  void
151  iterate_flux (const int& i,
152  const int& j,
153  const int& k,
154  const int& max_iters,
155  const amrex::Array4<amrex::Real>& z0_arr,
156  const amrex::Array4<const amrex::Real>& umm_arr,
157  const amrex::Array4<const amrex::Real>& tm_arr,
158  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
159  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
160  const amrex::Array4<amrex::Real>& u_star_arr,
161  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
162  const amrex::Array4<amrex::Real>& t_star_arr,
163  const amrex::Array4<amrex::Real>& /*q_star_arr*/,
164  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
165  const amrex::Array4<amrex::Real>& olen_arr,
166  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
167  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
168  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
169  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
170  {
171  int iter = 0;
172  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
173  amrex::Real ustar = 0.0;
174  amrex::Real z0 = 0.0;
175  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
176  do {
177  ustar = u_star_arr(i,j,k);
178  if (mdata.Cnk_a > 0) {
179  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
180  if (mdata.visc) {
181  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, 0.05);
182  }
183  } else {
184  z0 = COARE3_roughness(mdata.zref, umm, ustar);
185  }
186  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0);
187  ++iter;
188  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
189 
190  t_star_arr(i,j,k) = 0.0;
191  olen_arr(i,j,k) = 1.0e16;
192  z0_arr(i,j,k) = z0;
193  }
194 
195 private:
198  const amrex::Real tol = 1.0e-5;
199  const amrex::Real WSMIN = 0.1; // minimum wind speed
200 };
201 
202 
203 /**
204  * Adiabatic with modified charnock roughness
205  */
207 {
208  adiabatic_mod_charnock (amrex::Real zref,
209  amrex::Real flux,
210  amrex::Real depth)
211  {
212  mdata.zref = zref;
213  mdata.surf_temp_flux = flux;
214  mdata.Cnk_d = depth;
215  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
216  }
217 
218  AMREX_GPU_DEVICE
219  AMREX_FORCE_INLINE
220  void
221  iterate_flux (const int& i,
222  const int& j,
223  const int& k,
224  const int& max_iters,
225  const amrex::Array4<amrex::Real>& z0_arr,
226  const amrex::Array4<const amrex::Real>& umm_arr,
227  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
228  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
229  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
230  const amrex::Array4<amrex::Real>& u_star_arr,
231  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
232  const amrex::Array4<amrex::Real>& t_star_arr,
233  const amrex::Array4<amrex::Real>& /*q_star_arr*/,
234  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
235  const amrex::Array4<amrex::Real>& olen_arr,
236  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
237  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
238  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
239  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
240  {
241  int iter = 0;
242  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
243  amrex::Real ustar = 0.0;
244  amrex::Real z0 = 0.0;
245  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
246  do {
247  ustar = u_star_arr(i,j,k);
248  z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) );
249  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0);
250  ++iter;
251  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
252 
253  t_star_arr(i,j,k) = 0.0;
254  olen_arr(i,j,k) = 1.0e16;
255  z0_arr(i,j,k) = z0;
256  }
257 
258 private:
261  const amrex::Real tol = 1.0e-5;
262  const amrex::Real WSMIN = 0.1; // minimum wind speed
263 };
264 
265 
266 /**
267  * Adiabatic with Donelan roughness
268  */
270 {
271  adiabatic_donelan (amrex::Real zref,
272  amrex::Real flux)
273  {
274  mdata.zref = zref;
275  mdata.surf_temp_flux = flux;
276  }
277 
278  AMREX_GPU_DEVICE
279  AMREX_FORCE_INLINE
280  void
281  iterate_flux (const int& i,
282  const int& j,
283  const int& k,
284  const int& max_iters,
285  const amrex::Array4<amrex::Real>& z0_arr,
286  const amrex::Array4<const amrex::Real>& umm_arr,
287  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
288  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
289  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
290  const amrex::Array4<amrex::Real>& u_star_arr,
291  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
292  const amrex::Array4<amrex::Real>& t_star_arr,
293  const amrex::Array4<amrex::Real>& /*q_star_arr*/,
294  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
295  const amrex::Array4<amrex::Real>& olen_arr,
296  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
297  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
298  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
299  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
300  {
301  int iter = 0;
302  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
303  amrex::Real ustar = 0.0;
304  amrex::Real z0 = 0.0;
305  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
306  do {
307  ustar = u_star_arr(i,j,k);
308  z0 = Donelan_roughness(ustar);
309  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0);
310  ++iter;
311  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
312 
313  t_star_arr(i,j,k) = 0.0;
314  olen_arr(i,j,k) = 1.0e16;
315  z0_arr(i,j,k) = z0;
316  }
317 
318 private:
321  const amrex::Real tol = 1.0e-5;
322  const amrex::Real WSMIN = 0.1; // minimum wind speed
323 };
324 
325 
326 /**
327  * Adiabatic with wave-coupled roughness
328  */
330 {
331  adiabatic_wave_coupled (amrex::Real zref,
332  amrex::Real flux)
333  {
334  mdata.zref = zref;
335  mdata.surf_temp_flux = flux;
336  }
337 
338  AMREX_GPU_DEVICE
339  AMREX_FORCE_INLINE
340  void
341  iterate_flux (const int& i,
342  const int& j,
343  const int& k,
344  const int& max_iters,
345  const amrex::Array4<amrex::Real>& z0_arr,
346  const amrex::Array4<const amrex::Real>& umm_arr,
347  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
348  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
349  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
350  const amrex::Array4<amrex::Real>& u_star_arr,
351  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
352  const amrex::Array4<amrex::Real>& t_star_arr,
353  const amrex::Array4<amrex::Real>& /*q_star_arr*/,
354  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
355  const amrex::Array4<amrex::Real>& olen_arr,
356  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
357  const amrex::Array4<amrex::Real>& Hwave_arr,
358  const amrex::Array4<amrex::Real>& Lwave_arr,
359  const amrex::Array4<amrex::Real>& eta_arr) const
360  {
361  int iter = 0;
362  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
363  amrex::Real ustar = 0.0;
364  amrex::Real z0 = 0.0;
365  int ie, je;
366  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
367  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
368  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
369  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
370  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
371  do {
372  ustar = u_star_arr(i,j,k);
373  z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
374  + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
375  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0);
376  ++iter;
377  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
378 
379  t_star_arr(i,j,k) = 0.0;
380  olen_arr(i,j,k) = 1.0e16;
381  z0_arr(i,j,k) = z0;
382  }
383 
384 private:
387  const amrex::Real tol = 1.0e-5;
388  const amrex::Real eps = 1e-15;
389  const amrex::Real z0_eps = 1.0e-6;
390  const amrex::Real z0_max = 0.1;
391  const amrex::Real WSMIN = 0.1; // minimum wind speed
392 };
393 
394 
395 /**
396  * Surface flux with constant roughness
397  */
399 {
400  surface_flux (amrex::Real zref,
401  amrex::Real flux)
402  {
403  mdata.zref = zref;
404  mdata.surf_temp_flux = flux;
405  }
406 
407  AMREX_GPU_DEVICE
408  AMREX_FORCE_INLINE
409  void
410  iterate_flux (const int& i,
411  const int& j,
412  const int& k,
413  const int& max_iters,
414  const amrex::Array4<const amrex::Real>& z0_arr,
415  const amrex::Array4<const amrex::Real>& umm_arr,
416  const amrex::Array4<const amrex::Real>& tm_arr,
417  const amrex::Array4<const amrex::Real>& tvm_arr,
418  const amrex::Array4<const amrex::Real>& qvm_arr,
419  const amrex::Array4<amrex::Real>& u_star_arr,
420  const amrex::Array4<amrex::Real>& w_star_arr,
421  const amrex::Array4<amrex::Real>& t_star_arr,
422  const amrex::Array4<amrex::Real>& q_star_arr,
423  const amrex::Array4<amrex::Real>& t_surf_arr,
424  const amrex::Array4<amrex::Real>& olen_arr,
425  const amrex::Array4<amrex::Real>& pblh_arr,
426  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
427  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
428  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
429  {
430  int iter = 0;
431  amrex::Real ustar = 0.0;
432  amrex::Real wstar = 0.0;
433  amrex::Real tflux = 0.0;
434  amrex::Real zeta = 0.0;
435  amrex::Real psi_m = 0.0;
436  amrex::Real psi_h = 0.0;
437  amrex::Real Olen = 0.0;
438  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
439  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
440  do {
441  ustar = u_star_arr(i,j,k);
442  tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k);
443  if (w_star_arr) {
444  // update w* and Umagmean
445  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
446  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
447  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
448  umm = std::max(umm, WSMIN);
449  }
450  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
451  zeta = mdata.zref / Olen;
452  psi_m = sfuns.calc_psi_m(zeta);
453  psi_h = sfuns.calc_psi_h(zeta);
454  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0_arr(i,j,k)) - psi_m);
455  ++iter;
456  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
457  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
458 
459  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h) /
460  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
461  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
462  olen_arr(i,j,k) = Olen;
463  }
464 
465 private:
468  const amrex::Real tol = 1.0e-5;
469  const amrex::Real WSMIN = 0.1; // minimum wind speed
470 };
471 
472 
473 /**
474  * Surface flux with charnock roughness
475  */
477 {
478  surface_flux_charnock (amrex::Real zref,
479  amrex::Real flux,
480  amrex::Real cnk_a,
481  bool cnk_visc)
482  {
483  mdata.zref = zref;
484  mdata.surf_temp_flux = flux;
485  mdata.Cnk_a = cnk_a;
486  mdata.visc = cnk_visc;
487  }
488 
489  AMREX_GPU_DEVICE
490  AMREX_FORCE_INLINE
491  void
492  iterate_flux (const int& i,
493  const int& j,
494  const int& k,
495  const int& max_iters,
496  const amrex::Array4<amrex::Real>& z0_arr,
497  const amrex::Array4<const amrex::Real>& umm_arr,
498  const amrex::Array4<const amrex::Real>& tm_arr,
499  const amrex::Array4<const amrex::Real>& tvm_arr,
500  const amrex::Array4<const amrex::Real>& qvm_arr,
501  const amrex::Array4<amrex::Real>& u_star_arr,
502  const amrex::Array4<amrex::Real>& w_star_arr,
503  const amrex::Array4<amrex::Real>& t_star_arr,
504  const amrex::Array4<amrex::Real>& q_star_arr,
505  const amrex::Array4<amrex::Real>& t_surf_arr,
506  const amrex::Array4<amrex::Real>& olen_arr,
507  const amrex::Array4<amrex::Real>& pblh_arr,
508  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
509  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
510  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
511  {
512  int iter = 0;
513  amrex::Real ustar = 0.0;
514  amrex::Real wstar = 0.0;
515  amrex::Real tflux = 0.0;
516  amrex::Real z0 = 0.0;
517  amrex::Real zeta = 0.0;
518  amrex::Real psi_m = 0.0;
519  amrex::Real psi_h = 0.0;
520  amrex::Real Olen = 0.0;
521  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
522  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
523  do {
524  ustar = u_star_arr(i,j,k);
525  if (mdata.Cnk_a > 0) {
526  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
527  if (mdata.visc) {
528  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, 0.05);
529  }
530  } else {
531  z0 = COARE3_roughness(mdata.zref, umm, ustar);
532  }
533  tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k);
534  if (w_star_arr) {
535  // update w* and Umagmean
536  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
537  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
538  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
539  umm = std::max(umm, WSMIN);
540  }
541  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
542  zeta = mdata.zref / Olen;
543  psi_m = sfuns.calc_psi_m(zeta);
544  psi_h = sfuns.calc_psi_h(zeta);
545  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
546  ++iter;
547  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
548  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
549 
550  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
551  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
552  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
553  olen_arr(i,j,k) = Olen;
554  z0_arr(i,j,k) = z0;
555  }
556 
557 private:
560  const amrex::Real tol = 1.0e-5;
561  const amrex::Real WSMIN = 0.1; // minimum wind speed
562 };
563 
564 
565 /**
566  * Surface flux with modified charnock roughness
567  */
569 {
570  surface_flux_mod_charnock (amrex::Real zref,
571  amrex::Real flux,
572  amrex::Real depth)
573  {
574  mdata.zref = zref;
575  mdata.surf_temp_flux = flux;
576  mdata.Cnk_d = depth;
577  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
578  }
579 
580  AMREX_GPU_DEVICE
581  AMREX_FORCE_INLINE
582  void
583  iterate_flux (const int& i,
584  const int& j,
585  const int& k,
586  const int& max_iters,
587  const amrex::Array4<amrex::Real>& z0_arr,
588  const amrex::Array4<const amrex::Real>& umm_arr,
589  const amrex::Array4<const amrex::Real>& tm_arr,
590  const amrex::Array4<const amrex::Real>& tvm_arr,
591  const amrex::Array4<const amrex::Real>& qvm_arr,
592  const amrex::Array4<amrex::Real>& u_star_arr,
593  const amrex::Array4<amrex::Real>& w_star_arr,
594  const amrex::Array4<amrex::Real>& t_star_arr,
595  const amrex::Array4<amrex::Real>& q_star_arr,
596  const amrex::Array4<amrex::Real>& t_surf_arr,
597  const amrex::Array4<amrex::Real>& olen_arr,
598  const amrex::Array4<amrex::Real>& pblh_arr,
599  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
600  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
601  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
602  {
603  int iter = 0;
604  amrex::Real ustar = 0.0;
605  amrex::Real wstar = 0.0;
606  amrex::Real tflux = 0.0;
607  amrex::Real z0 = 0.0;
608  amrex::Real zeta = 0.0;
609  amrex::Real psi_m = 0.0;
610  amrex::Real psi_h = 0.0;
611  amrex::Real Olen = 0.0;
612  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
613  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
614  do {
615  ustar = u_star_arr(i,j,k);
616  z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) );
617  tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k);
618  if (w_star_arr) {
619  // update w* and Umagmean
620  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
621  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
622  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
623  umm = std::max(umm, WSMIN);
624  }
625  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
626  zeta = mdata.zref / Olen;
627  psi_m = sfuns.calc_psi_m(zeta);
628  psi_h = sfuns.calc_psi_h(zeta);
629  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
630  ++iter;
631  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
632  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
633 
634  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
635  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
636  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
637  olen_arr(i,j,k) = Olen;
638  z0_arr(i,j,k) = z0;
639  }
640 
641 private:
644  const amrex::Real tol = 1.0e-5;
645  const amrex::Real WSMIN = 0.1; // minimum wind speed
646 };
647 
648 
649 /**
650  * Surface flux with donelan roughness
651  */
653 {
654  surface_flux_donelan (amrex::Real zref,
655  amrex::Real flux)
656  {
657  mdata.zref = zref;
658  mdata.surf_temp_flux = flux;
659  }
660 
661  AMREX_GPU_DEVICE
662  AMREX_FORCE_INLINE
663  void
664  iterate_flux (const int& i,
665  const int& j,
666  const int& k,
667  const int& max_iters,
668  const amrex::Array4<amrex::Real>& z0_arr,
669  const amrex::Array4<const amrex::Real>& umm_arr,
670  const amrex::Array4<const amrex::Real>& tm_arr,
671  const amrex::Array4<const amrex::Real>& tvm_arr,
672  const amrex::Array4<const amrex::Real>& qvm_arr,
673  const amrex::Array4<amrex::Real>& u_star_arr,
674  const amrex::Array4<amrex::Real>& w_star_arr,
675  const amrex::Array4<amrex::Real>& t_star_arr,
676  const amrex::Array4<amrex::Real>& q_star_arr,
677  const amrex::Array4<amrex::Real>& t_surf_arr,
678  const amrex::Array4<amrex::Real>& olen_arr,
679  const amrex::Array4<amrex::Real>& pblh_arr,
680  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
681  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
682  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
683  {
684  int iter = 0;
685  amrex::Real ustar = 0.0;
686  amrex::Real wstar = 0.0;
687  amrex::Real tflux = 0.0;
688  amrex::Real z0 = 0.0;
689  amrex::Real zeta = 0.0;
690  amrex::Real psi_m = 0.0;
691  amrex::Real psi_h = 0.0;
692  amrex::Real Olen = 0.0;
693  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
694  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
695  do {
696  ustar = u_star_arr(i,j,k);
697  z0 = Donelan_roughness(ustar);
698  tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k);
699  if (w_star_arr) {
700  // update w* and Umagmean
701  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
702  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
703  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
704  umm = std::max(umm, WSMIN);
705  }
706  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
707  zeta = mdata.zref / Olen;
708  psi_m = sfuns.calc_psi_m(zeta);
709  psi_h = sfuns.calc_psi_h(zeta);
710  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
711  ++iter;
712  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
713  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
714 
715  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
716  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
717  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
718  olen_arr(i,j,k) = Olen;
719  z0_arr(i,j,k) = z0;
720  }
721 
722 private:
725  const amrex::Real tol = 1.0e-5;
726  const amrex::Real WSMIN = 0.1; // minimum wind speed
727 };
728 
729 
730 /**
731  * Surface flux with wave-coupled roughness
732  */
734 {
735  surface_flux_wave_coupled (amrex::Real zref,
736  amrex::Real flux)
737  {
738  mdata.zref = zref;
739  mdata.surf_temp_flux = flux;
740  }
741 
742  AMREX_GPU_DEVICE
743  AMREX_FORCE_INLINE
744  void
745  iterate_flux (const int& i,
746  const int& j,
747  const int& k,
748  const int& max_iters,
749  const amrex::Array4<amrex::Real>& z0_arr,
750  const amrex::Array4<const amrex::Real>& umm_arr,
751  const amrex::Array4<const amrex::Real>& tm_arr,
752  const amrex::Array4<const amrex::Real>& tvm_arr,
753  const amrex::Array4<const amrex::Real>& qvm_arr,
754  const amrex::Array4<amrex::Real>& u_star_arr,
755  const amrex::Array4<amrex::Real>& w_star_arr,
756  const amrex::Array4<amrex::Real>& t_star_arr,
757  const amrex::Array4<amrex::Real>& q_star_arr,
758  const amrex::Array4<amrex::Real>& t_surf_arr,
759  const amrex::Array4<amrex::Real>& olen_arr,
760  const amrex::Array4<amrex::Real>& pblh_arr,
761  const amrex::Array4<amrex::Real>& Hwave_arr,
762  const amrex::Array4<amrex::Real>& Lwave_arr,
763  const amrex::Array4<amrex::Real>& eta_arr) const
764  {
765  int iter = 0;
766  amrex::Real ustar = 0.0;
767  amrex::Real wstar = 0.0;
768  amrex::Real tflux = 0.0;
769  amrex::Real z0 = 0.0;
770  amrex::Real zeta = 0.0;
771  amrex::Real psi_m = 0.0;
772  amrex::Real psi_h = 0.0;
773  amrex::Real Olen = 0.0;
774  int ie, je;
775  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
776  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
777  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
778  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
779  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
780  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
781  do {
782  ustar = u_star_arr(i,j,k);
783  z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
784  + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
785  tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k);
786  if (w_star_arr) {
787  // update w* and Umagmean
788  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
789  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
790  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
791  umm = std::max(umm, WSMIN);
792  }
793  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
794  zeta = mdata.zref / Olen;
795  psi_m = sfuns.calc_psi_m(zeta);
796  psi_h = sfuns.calc_psi_h(zeta);
797  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
798  ++iter;
799  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
800  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
801 
802  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
803  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
804  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
805  olen_arr(i,j,k) = Olen;
806  z0_arr(i,j,k) = z0;
807  }
808 
809 private:
812  const amrex::Real tol = 1.0e-5;
813  const amrex::Real eps = 1e-15;
814  const amrex::Real z0_eps = 1.0e-6;
815  const amrex::Real z0_max = 0.1;
816  const amrex::Real WSMIN = 0.1; // minimum wind speed
817 };
818 
819 
820 /**
821  * Surface temperature with constant roughness
822  */
824 {
825  surface_temp (amrex::Real zref,
826  amrex::Real flux)
827  {
828  mdata.zref = zref;
829  mdata.surf_temp_flux = flux;
830  }
831 
832  AMREX_GPU_DEVICE
833  AMREX_FORCE_INLINE
834  void
835  iterate_flux (const int& i,
836  const int& j,
837  const int& k,
838  const int& max_iters,
839  const amrex::Array4<const amrex::Real>& z0_arr,
840  const amrex::Array4<const amrex::Real>& umm_arr,
841  const amrex::Array4<const amrex::Real>& tm_arr,
842  const amrex::Array4<const amrex::Real>& tvm_arr,
843  const amrex::Array4<const amrex::Real>& qvm_arr,
844  const amrex::Array4<amrex::Real>& u_star_arr,
845  const amrex::Array4<amrex::Real>& w_star_arr,
846  const amrex::Array4<amrex::Real>& t_star_arr,
847  const amrex::Array4<amrex::Real>& q_star_arr,
848  const amrex::Array4<amrex::Real>& t_surf_arr,
849  const amrex::Array4<amrex::Real>& olen_arr,
850  const amrex::Array4<amrex::Real>& pblh_arr,
851  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
852  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
853  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
854  {
855  int iter = 0;
856  amrex::Real ustar = 0.0;
857  amrex::Real wstar = 0.0;
858  amrex::Real tflux = 0.0;
859  amrex::Real zeta = 0.0;
860  amrex::Real psi_m = 0.0;
861  amrex::Real psi_h = 0.0;
862  amrex::Real Olen = 0.0;
863  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
864  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
865  do {
866  ustar = u_star_arr(i,j,k);
867  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
868  (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h); // <w'T'>
869  tflux *= (1 + 0.61*qvm_arr(i,j,k));
870  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
871  if (w_star_arr) {
872  // update w* and Umagmean
873  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
874  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
875  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
876  umm = std::max(umm, WSMIN);
877  }
878  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
879  zeta = mdata.zref / Olen;
880  psi_m = sfuns.calc_psi_m(zeta);
881  psi_h = sfuns.calc_psi_h(zeta);
882  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0_arr(i,j,k)) - psi_m);
883  ++iter;
884  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
885  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
886 
887  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
888  (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h);
889  olen_arr(i,j,k) = Olen;
890  }
891 
892 private:
895  const amrex::Real tol = 1.0e-5;
896  const amrex::Real WSMIN = 0.1; // minimum wind speed
897 };
898 
899 
900 /**
901  * Surface temperature with charnock roughness
902  */
904 {
905  surface_temp_charnock (amrex::Real zref,
906  amrex::Real flux,
907  amrex::Real cnk_a,
908  bool cnk_visc)
909  {
910  mdata.zref = zref;
911  mdata.surf_temp_flux = flux;
912  mdata.Cnk_a = cnk_a;
913  mdata.visc = cnk_visc;
914  }
915 
916  AMREX_GPU_DEVICE
917  AMREX_FORCE_INLINE
918  void
919  iterate_flux (const int& i,
920  const int& j,
921  const int& k,
922  const int& max_iters,
923  const amrex::Array4<amrex::Real>& z0_arr,
924  const amrex::Array4<const amrex::Real>& umm_arr,
925  const amrex::Array4<const amrex::Real>& tm_arr,
926  const amrex::Array4<const amrex::Real>& tvm_arr,
927  const amrex::Array4<const amrex::Real>& qvm_arr,
928  const amrex::Array4<amrex::Real>& u_star_arr,
929  const amrex::Array4<amrex::Real>& w_star_arr,
930  const amrex::Array4<amrex::Real>& t_star_arr,
931  const amrex::Array4<amrex::Real>& q_star_arr,
932  const amrex::Array4<amrex::Real>& t_surf_arr,
933  const amrex::Array4<amrex::Real>& olen_arr,
934  const amrex::Array4<amrex::Real>& pblh_arr,
935  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
936  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
937  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
938  {
939  int iter = 0;
940  amrex::Real ustar = 0.0;
941  amrex::Real wstar = 0.0;
942  amrex::Real z0 = 0.0;
943  amrex::Real tflux = 0.0;
944  amrex::Real zeta = 0.0;
945  amrex::Real psi_m = 0.0;
946  amrex::Real psi_h = 0.0;
947  amrex::Real Olen = 0.0;
948  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
949  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
950  do {
951  ustar = u_star_arr(i,j,k);
952  if (mdata.Cnk_a > 0) {
953  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
954  if (mdata.visc) {
955  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, 0.05);
956  }
957  } else {
958  z0 = COARE3_roughness(mdata.zref, umm, ustar);
959  }
960  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
961  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
962  tflux *= (1 + 0.61*qvm_arr(i,j,k));
963  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
964  if (w_star_arr) {
965  // update w* and Umagmean
966  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
967  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
968  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
969  umm = std::max(umm, WSMIN);
970  }
971  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
972  zeta = mdata.zref / Olen;
973  psi_m = sfuns.calc_psi_m(zeta);
974  psi_h = sfuns.calc_psi_h(zeta);
975  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
976  ++iter;
977  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
978  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
979 
980  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
981  (std::log(mdata.zref / z0) - psi_h);
982  olen_arr(i,j,k) = Olen;
983  z0_arr(i,j,k) = z0;
984  }
985 
986 private:
989  const amrex::Real tol = 1.0e-5;
990  const amrex::Real WSMIN = 0.1; // minimum wind speed
991 };
992 
993 
994 /**
995  * Surface temperature with modified charnock roughness
996  */
998 {
999  surface_temp_mod_charnock (amrex::Real zref,
1000  amrex::Real flux,
1001  amrex::Real depth)
1002  {
1003  mdata.zref = zref;
1004  mdata.surf_temp_flux = flux;
1005  mdata.Cnk_d = depth;
1006  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
1007  }
1008 
1009  AMREX_GPU_DEVICE
1010  AMREX_FORCE_INLINE
1011  void
1012  iterate_flux (const int& i,
1013  const int& j,
1014  const int& k,
1015  const int& max_iters,
1016  const amrex::Array4<amrex::Real>& z0_arr,
1017  const amrex::Array4<const amrex::Real>& umm_arr,
1018  const amrex::Array4<const amrex::Real>& tm_arr,
1019  const amrex::Array4<const amrex::Real>& tvm_arr,
1020  const amrex::Array4<const amrex::Real>& qvm_arr,
1021  const amrex::Array4<amrex::Real>& u_star_arr,
1022  const amrex::Array4<amrex::Real>& w_star_arr,
1023  const amrex::Array4<amrex::Real>& t_star_arr,
1024  const amrex::Array4<amrex::Real>& q_star_arr,
1025  const amrex::Array4<amrex::Real>& t_surf_arr,
1026  const amrex::Array4<amrex::Real>& olen_arr,
1027  const amrex::Array4<amrex::Real>& pblh_arr,
1028  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1029  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1030  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1031  {
1032  int iter = 0;
1033  amrex::Real ustar = 0.0;
1034  amrex::Real wstar = 0.0;
1035  amrex::Real z0 = 0.0;
1036  amrex::Real tflux = 0.0;
1037  amrex::Real zeta = 0.0;
1038  amrex::Real psi_m = 0.0;
1039  amrex::Real psi_h = 0.0;
1040  amrex::Real Olen = 0.0;
1041  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1042  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
1043  do {
1044  ustar = u_star_arr(i,j,k);
1045  z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) );
1046  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1047  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
1048  tflux *= (1 + 0.61*qvm_arr(i,j,k));
1049  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
1050  if (w_star_arr) {
1051  // update w* and Umagmean
1052  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1053  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1054  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1055  umm = std::max(umm, WSMIN);
1056  }
1057  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1058  zeta = mdata.zref / Olen;
1059  psi_m = sfuns.calc_psi_m(zeta);
1060  psi_h = sfuns.calc_psi_h(zeta);
1061  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
1062  ++iter;
1063  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1064  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
1065 
1066  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1067  (std::log(mdata.zref / z0) - psi_h);
1068  olen_arr(i,j,k) = Olen;
1069  z0_arr(i,j,k) = z0;
1070  }
1071 
1072 private:
1075  const amrex::Real tol = 1.0e-5;
1076  const amrex::Real WSMIN = 0.1; // minimum wind speed
1077 };
1078 
1079 
1080 /**
1081  * Surface temperature with donelan roughness
1082  */
1084 {
1085  surface_temp_donelan (amrex::Real zref,
1086  amrex::Real flux)
1087  {
1088  mdata.zref = zref;
1089  mdata.surf_temp_flux = flux;
1090  }
1091 
1092  AMREX_GPU_DEVICE
1093  AMREX_FORCE_INLINE
1094  void
1095  iterate_flux (const int& i,
1096  const int& j,
1097  const int& k,
1098  const int& max_iters,
1099  const amrex::Array4<amrex::Real>& z0_arr,
1100  const amrex::Array4<const amrex::Real>& umm_arr,
1101  const amrex::Array4<const amrex::Real>& tm_arr,
1102  const amrex::Array4<const amrex::Real>& tvm_arr,
1103  const amrex::Array4<const amrex::Real>& qvm_arr,
1104  const amrex::Array4<amrex::Real>& u_star_arr,
1105  const amrex::Array4<amrex::Real>& w_star_arr,
1106  const amrex::Array4<amrex::Real>& t_star_arr,
1107  const amrex::Array4<amrex::Real>& q_star_arr,
1108  const amrex::Array4<amrex::Real>& t_surf_arr,
1109  const amrex::Array4<amrex::Real>& olen_arr,
1110  const amrex::Array4<amrex::Real>& pblh_arr,
1111  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1112  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1113  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1114  {
1115  int iter = 0;
1116  amrex::Real ustar = 0.0;
1117  amrex::Real wstar = 0.0;
1118  amrex::Real z0 = 0.0;
1119  amrex::Real tflux = 0.0;
1120  amrex::Real zeta = 0.0;
1121  amrex::Real psi_m = 0.0;
1122  amrex::Real psi_h = 0.0;
1123  amrex::Real Olen = 0.0;
1124  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1125  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
1126  do {
1127  ustar = u_star_arr(i,j,k);
1128  z0 = Donelan_roughness(ustar);
1129  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1130  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
1131  tflux *= (1 + 0.61*qvm_arr(i,j,k));
1132  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
1133  if (w_star_arr) {
1134  // update w* and Umagmean
1135  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1136  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1137  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1138  umm = std::max(umm, WSMIN);
1139  }
1140  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1141  zeta = mdata.zref / Olen;
1142  psi_m = sfuns.calc_psi_m(zeta);
1143  psi_h = sfuns.calc_psi_h(zeta);
1144  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
1145  ++iter;
1146  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1147  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
1148 
1149  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1150  (std::log(mdata.zref / z0) - psi_h);
1151  olen_arr(i,j,k) = Olen;
1152  z0_arr(i,j,k) = z0;
1153  }
1154 
1155 private:
1158  const amrex::Real tol = 1.0e-5;
1159  const amrex::Real WSMIN = 0.1; // minimum wind speed
1160 };
1161 
1162 
1163 /**
1164  * Surface temperature with wave-coupled roughness
1165  */
1167 {
1168  surface_temp_wave_coupled (amrex::Real zref,
1169  amrex::Real flux)
1170  {
1171  mdata.zref = zref;
1172  mdata.surf_temp_flux = flux;
1173  }
1174 
1175  AMREX_GPU_DEVICE
1176  AMREX_FORCE_INLINE
1177  void
1178  iterate_flux (const int& i,
1179  const int& j,
1180  const int& k,
1181  const int& max_iters,
1182  const amrex::Array4<amrex::Real>& z0_arr,
1183  const amrex::Array4<const amrex::Real>& umm_arr,
1184  const amrex::Array4<const amrex::Real>& tm_arr,
1185  const amrex::Array4<const amrex::Real>& tvm_arr,
1186  const amrex::Array4<const amrex::Real>& qvm_arr,
1187  const amrex::Array4<amrex::Real>& u_star_arr,
1188  const amrex::Array4<amrex::Real>& w_star_arr,
1189  const amrex::Array4<amrex::Real>& t_star_arr,
1190  const amrex::Array4<amrex::Real>& q_star_arr,
1191  const amrex::Array4<amrex::Real>& t_surf_arr,
1192  const amrex::Array4<amrex::Real>& olen_arr,
1193  const amrex::Array4<amrex::Real>& pblh_arr,
1194  const amrex::Array4<amrex::Real>& Hwave_arr,
1195  const amrex::Array4<amrex::Real>& Lwave_arr,
1196  const amrex::Array4<amrex::Real>& eta_arr) const
1197  {
1198  int iter = 0;
1199  amrex::Real ustar = 0.0;
1200  amrex::Real wstar = 0.0;
1201  amrex::Real z0 = 0.0;
1202  amrex::Real tflux = 0.0;
1203  amrex::Real zeta = 0.0;
1204  amrex::Real psi_m = 0.0;
1205  amrex::Real psi_h = 0.0;
1206  amrex::Real Olen = 0.0;
1207  int ie, je;
1208  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1209  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1210  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1211  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1212  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1213  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
1214  do {
1215  ustar = u_star_arr(i,j,k);
1216  z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
1217  + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
1218  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1219  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
1220  tflux *= (1 + 0.61*qvm_arr(i,j,k));
1221  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
1222  if (w_star_arr) {
1223  // update w* and Umagmean
1224  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1225  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1226  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1227  umm = std::max(umm, WSMIN);
1228  }
1229  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1230  zeta = mdata.zref / Olen;
1231  psi_m = sfuns.calc_psi_m(zeta);
1232  psi_h = sfuns.calc_psi_h(zeta);
1233  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
1234  ++iter;
1235  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1236  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
1237 
1238  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1239  (std::log(mdata.zref / z0) - psi_h);
1240  olen_arr(i,j,k) = Olen;
1241  z0_arr(i,j,k) = z0;
1242  }
1243 
1244 private:
1247  const amrex::Real tol = 1.0e-5;
1248  const amrex::Real eps = 1e-15;
1249  const amrex::Real z0_eps = 1.0e-6;
1250  const amrex::Real z0_max = 0.1;
1251  const amrex::Real WSMIN = 0.1; // minimum wind speed
1252 };
1253 
1254 
1255 /**
1256  * Moeng flux formulation
1257  */
1259 {
1260  moeng_flux (int l_zlo)
1261  : zlo(l_zlo) {}
1262 
1263 
1264  AMREX_GPU_DEVICE
1265  AMREX_FORCE_INLINE
1266  amrex::Real
1267  compute_q_flux (const int& i,
1268  const int& j,
1269  const int& k,
1270  const int& n,
1271  const int& icomp,
1272  const amrex::Real& dz,
1273  const amrex::Real& dz1,
1274  const bool& exp_most,
1275  const amrex::Array4<const amrex::Real>& eta_arr,
1276  const amrex::Array4<const amrex::Real>& cons_arr,
1277  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1278  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1279  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1280  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
1281  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1282  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
1283  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1284  const amrex::Array4<amrex::Real>& dest_arr) const
1285  {
1286  amrex::Real rho, eta;
1287 
1288  int ic, jc;
1289  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1290  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1291  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1292  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1293 
1294  rho = cons_arr(ic,jc,zlo,Rho_comp);
1295 
1296  // TODO: Integrate MOST with moisture and MOENG FLUX type
1297  amrex::Real deltaz = dz * (zlo - k);
1298 
1299  // NOTE: this is rho*<q'w'> = -K dqdz
1300  amrex::Real moflux = 0.0;
1301 
1302  if (exp_most) {
1303  // surface gradient equal to gradient at first zface
1304  amrex::Real rqvgrad = (cons_arr(ic,jc,zlo+1,RhoQ1_comp) - cons_arr(ic,jc,zlo ,RhoQ1_comp)) / (0.5*(dz+dz1));
1305  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
1306  } else {
1307  int ie, je;
1308  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1309  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1310  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1311  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1312  eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
1313  eta = amrex::max(eta,eta_eps);
1314  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1315  }
1316 
1317  return moflux;
1318  }
1319 
1320  AMREX_GPU_DEVICE
1321  AMREX_FORCE_INLINE
1322  amrex::Real
1323  compute_t_flux (const int& i,
1324  const int& j,
1325  const int& k,
1326  const int& n,
1327  const int& icomp,
1328  const amrex::Real& dz,
1329  const amrex::Real& dz1,
1330  const bool& exp_most,
1331  const amrex::Array4<const amrex::Real>& eta_arr,
1332  const amrex::Array4<const amrex::Real>& cons_arr,
1333  const amrex::Array4<const amrex::Real>& velx_arr,
1334  const amrex::Array4<const amrex::Real>& vely_arr,
1335  const amrex::Array4<const amrex::Real>& umm_arr,
1336  const amrex::Array4<const amrex::Real>& tm_arr,
1337  const amrex::Array4<const amrex::Real>& u_star_arr,
1338  const amrex::Array4<const amrex::Real>& t_star_arr,
1339  const amrex::Array4<const amrex::Real>& t_surf_arr,
1340  const amrex::Array4<amrex::Real>& dest_arr) const
1341  {
1342  amrex::Real velx, vely, rho, theta, eta;
1343  int ix, jx, iy, jy, ic, jc;
1344 
1345  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
1346  jx = j < lbound(velx_arr).y ? lbound(velx_arr).y : j;
1347  ix = ix > ubound(velx_arr).x-1 ? ubound(velx_arr).x-1 : ix;
1348  jx = jx > ubound(velx_arr).y ? ubound(velx_arr).y : jx;
1349 
1350  iy = i < lbound(vely_arr).x ? lbound(vely_arr).x : i;
1351  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1352  iy = iy > ubound(vely_arr).x ? ubound(vely_arr).x : iy;
1353  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1354 
1355  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1356  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1357  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1358  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1359 
1360  velx = 0.5 *( velx_arr(ix,jx,zlo) + velx_arr(ix+1,jx ,zlo) );
1361  vely = 0.5 *( vely_arr(iy,jy,zlo) + vely_arr(iy ,jy+1,zlo) );
1362  rho = cons_arr(ic,jc,zlo,Rho_comp);
1363  theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho;
1364 
1365  amrex::Real theta_mean = tm_arr(ic,jc,zlo);
1366  amrex::Real wsp_mean = umm_arr(ic,jc,zlo);
1367  amrex::Real ustar = u_star_arr(ic,jc,zlo);
1368  amrex::Real tstar = t_star_arr(ic,jc,zlo);
1369  amrex::Real theta_surf = t_surf_arr(ic,jc,zlo);
1370 
1371  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1372  amrex::Real num1 = wsp * (theta_mean-theta_surf);
1373  amrex::Real num2 = wsp_mean * (theta-theta_mean);
1374  amrex::Real deltaz = dz * (zlo - k);
1375 
1376  wsp_mean = std::max(wsp_mean, WSMIN);
1377 
1378  // NOTE: this is rho*<T'w'> = -K dTdz
1379  amrex::Real moflux = (std::abs(tstar) > eps) ?
1380  -rho*tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0;
1381 
1382  if (exp_most) {
1383  // surface gradient equal to gradient at first zface
1384  amrex::Real rthetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1));
1385  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
1386  } else {
1387  int ie, je;
1388  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1389  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1390  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1391  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1392  eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
1393  eta = amrex::max(eta,eta_eps);
1394  // Note: Kh = eta/rho
1395  // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface
1396  // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel
1397  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1398  }
1399 
1400  return moflux;
1401  }
1402 
1403  AMREX_GPU_DEVICE
1404  AMREX_FORCE_INLINE
1405  amrex::Real
1406  compute_u_flux (const int& i,
1407  const int& j,
1408  const int& k,
1409  const int& icomp,
1410  const amrex::Real& dz,
1411  const amrex::Real& dz1,
1412  const bool& exp_most,
1413  const amrex::Array4<const amrex::Real>& eta_arr,
1414  const amrex::Array4<const amrex::Real>& cons_arr,
1415  const amrex::Array4<const amrex::Real>& velx_arr,
1416  const amrex::Array4<const amrex::Real>& vely_arr,
1417  const amrex::Array4<const amrex::Real>& umm_arr,
1418  const amrex::Array4<const amrex::Real>& um_arr,
1419  const amrex::Array4<const amrex::Real>& u_star_arr,
1420  const amrex::Array4<amrex::Real>& dest_arr) const
1421  {
1422  amrex::Real velx, vely, rho, eta;
1423  int jy, ic, jc;
1424 
1425  int iylo = i <= lbound(vely_arr).x ? lbound(vely_arr).x : i-1;
1426  int iyhi = i > ubound(vely_arr).x ? ubound(vely_arr).x : i;
1427 
1428  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1429  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1430 
1431  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
1432  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1433  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1434  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1435 
1436  velx = velx_arr(i,j,zlo);
1437  vely = 0.25*( vely_arr(iyhi,jy,zlo)+vely_arr(iyhi,jy+1,zlo)
1438  + vely_arr(iylo,jy,zlo)+vely_arr(iylo,jy+1,zlo) );
1439  rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
1440  + cons_arr(ic ,jc,zlo,Rho_comp) );
1441 
1442  amrex::Real umean = um_arr(i,j,zlo);
1443  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) );
1444  amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
1445 
1446  // Note: The surface mean shear stress is decomposed into tau_xz by
1447  // multiplying the modeled shear stress (rho*ustar^2) with
1448  // a factor of umean/wsp_mean for directionality; this factor
1449  // modifies the denominator from what is in Moeng 1984.
1450  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1451  amrex::Real num1 = wsp * umean;
1452  amrex::Real num2 = wsp_mean * (velx-umean);
1453  amrex::Real deltaz = dz * (zlo - k);
1454 
1455  wsp_mean = std::max(wsp_mean, WSMIN);
1456 
1457  // NOTE: this is rho*<u'w'> = -K dudz
1458  amrex::Real stressx = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
1459 
1460  if (exp_most) {
1461  // surface gradient equal to gradient at first zface
1462  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
1463  dest_arr(i,j,k,icomp) = velx - ugrad * deltaz;
1464  } else {
1465  int ie, je;
1466  ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
1467  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1468  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1469  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1470  eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
1471  + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) );
1472  eta = amrex::max(eta,eta_eps);
1473  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressx/eta*deltaz;
1474  }
1475 
1476  return stressx;
1477  }
1478 
1479  AMREX_GPU_DEVICE
1480  AMREX_FORCE_INLINE
1481  amrex::Real
1482  compute_v_flux (const int& i,
1483  const int& j,
1484  const int& k,
1485  const int& icomp,
1486  const amrex::Real& dz,
1487  const amrex::Real& dz1,
1488  const bool& exp_most,
1489  const amrex::Array4<const amrex::Real>& eta_arr,
1490  const amrex::Array4<const amrex::Real>& cons_arr,
1491  const amrex::Array4<const amrex::Real>& velx_arr,
1492  const amrex::Array4<const amrex::Real>& vely_arr,
1493  const amrex::Array4<const amrex::Real>& umm_arr,
1494  const amrex::Array4<const amrex::Real>& vm_arr,
1495  const amrex::Array4<const amrex::Real>& u_star_arr,
1496  const amrex::Array4<amrex::Real>& dest_arr) const
1497  {
1498  amrex::Real velx, vely, rho, eta;
1499  int ix, ic, jc;
1500 
1501  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
1502  ix = ix > ubound(velx_arr).x ? ubound(velx_arr).x : ix;
1503 
1504  int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
1505  int jxhi = j > ubound(velx_arr).y ? ubound(velx_arr).y : j;
1506 
1507  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1508  jc = j < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
1509  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1510  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1511 
1512  velx = 0.25*( velx_arr(ix,jxhi,zlo)+velx_arr(ix+1,jxhi,zlo)
1513  + velx_arr(ix,jxlo,zlo)+velx_arr(ix+1,jxlo,zlo) );
1514  vely = vely_arr(i,j,zlo);
1515  rho = 0.5*( cons_arr(ic,jc-1,zlo,Rho_comp)
1516  + cons_arr(ic,jc ,zlo,Rho_comp) );
1517 
1518  amrex::Real vmean = vm_arr(i,j,zlo);
1519  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) );
1520  amrex::Real ustar = 0.5 * ( u_star_arr(ic,jc-1,zlo) + u_star_arr(ic,jc,zlo) );
1521 
1522  // Note: The surface mean shear stress is decomposed into tau_yz by
1523  // multiplying the modeled shear stress (rho*ustar^2) with
1524  // a factor of vmean/wsp_mean for directionality; this factor
1525  // modifies the denominator from what is in Moeng 1984.
1526  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1527  amrex::Real num1 = wsp * vmean;
1528  amrex::Real num2 = wsp_mean * (vely-vmean);
1529  amrex::Real deltaz = dz * (zlo - k);
1530 
1531  wsp_mean = std::max(wsp_mean, WSMIN);
1532 
1533  // NOTE: this is rho*<v'w'> = -K dvdz
1534  amrex::Real stressy = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
1535 
1536  if (exp_most) {
1537  // surface gradient equal to gradient at first zface
1538  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
1539  dest_arr(i,j,k,icomp) = vely - vgrad * deltaz;
1540  } else {
1541  int ie, je;
1542  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1543  je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
1544  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1545  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1546  eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
1547  + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) );
1548  eta = amrex::max(eta,eta_eps);
1549  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressy/eta*deltaz;
1550  }
1551 
1552  return stressy;
1553  }
1554 
1555 private:
1556  int zlo;
1557  const amrex::Real eps = 1e-15;
1558  const amrex::Real eta_eps = 1e-8;
1559  const amrex::Real WSMIN = 0.1; // minimum wind speed
1560 };
1561 
1562 
1563 /**
1564  * Donelan flux formulation
1565  */
1567 {
1568  donelan_flux (int l_zlo)
1569  : zlo(l_zlo) {}
1570 
1571 
1572  AMREX_GPU_DEVICE
1573  AMREX_FORCE_INLINE
1574  amrex::Real
1575  compute_q_flux (const int& i,
1576  const int& j,
1577  const int& k,
1578  const int& n,
1579  const int& icomp,
1580  const amrex::Real& dz,
1581  const amrex::Real& dz1,
1582  const bool& exp_most,
1583  const amrex::Array4<const amrex::Real>& eta_arr,
1584  const amrex::Array4<const amrex::Real>& cons_arr,
1585  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1586  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1587  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1588  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
1589  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1590  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
1591  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1592  const amrex::Array4<amrex::Real>& dest_arr) const
1593  {
1594  amrex::Real rho, eta;
1595 
1596  int ic, jc;
1597  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1598  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1599  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1600  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1601 
1602  rho = cons_arr(ic,jc,zlo,Rho_comp);
1603 
1604  // TODO: Integrate MOST with moisture and DONELAN FLUX type
1605  amrex::Real deltaz = dz * (zlo - k);
1606 
1607  // NOTE: this is rho*<q'w'> = -K dqdz
1608  amrex::Real moflux = 0.0;
1609 
1610  if (exp_most) {
1611  // surface gradient equal to gradient at first zface
1612  amrex::Real rqvgrad = (cons_arr(ic,jc,zlo+1,RhoQ1_comp) - cons_arr(ic,jc,zlo ,RhoQ1_comp)) / (0.5*(dz+dz1));
1613  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
1614  } else {
1615  int ie, je;
1616  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1617  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1618  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1619  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1620  eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
1621  eta = amrex::max(eta,eta_eps);
1622  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1623  }
1624 
1625  return moflux;
1626  }
1627 
1628  AMREX_GPU_DEVICE
1629  AMREX_FORCE_INLINE
1630  amrex::Real
1631  compute_t_flux (const int& i,
1632  const int& j,
1633  const int& k,
1634  const int& n,
1635  const int& icomp,
1636  const amrex::Real& dz,
1637  const amrex::Real& dz1,
1638  const bool& exp_most,
1639  const amrex::Array4<const amrex::Real>& eta_arr,
1640  const amrex::Array4<const amrex::Real>& cons_arr,
1641  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1642  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1643  const amrex::Array4<const amrex::Real>& umm_arr,
1644  const amrex::Array4<const amrex::Real>& tm_arr,
1645  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1646  const amrex::Array4<const amrex::Real>& /*t_star_arr*/,
1647  const amrex::Array4<const amrex::Real>& t_surf_arr,
1648  const amrex::Array4<amrex::Real>& dest_arr) const
1649  {
1650  amrex::Real rho, eta;
1651 
1652  int ic, jc;
1653  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1654  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1655  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1656  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1657 
1658  rho = cons_arr(ic,jc,zlo,Rho_comp);
1659 
1660  amrex::Real Ch = 0.0012;
1661  amrex::Real wsp_mean = umm_arr(ic,jc,zlo);
1662  amrex::Real theta_surf = t_surf_arr(ic,jc,zlo);
1663  amrex::Real theta_mean = tm_arr(ic,jc,zlo);
1664  amrex::Real deltaz = dz * (zlo - k);
1665 
1666  // NOTE: this is rho*<T'w'> = -K dTdz
1667  amrex::Real moflux = -rho * Ch * wsp_mean * (theta_mean - theta_surf);
1668 
1669  if (exp_most) {
1670  // surface gradient equal to gradient at first zface
1671  amrex::Real rthetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1));
1672  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
1673  } else {
1674  int ie, je;
1675  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1676  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1677  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1678  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1679  eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
1680  eta = amrex::max(eta,eta_eps);
1681  // Note: Kh = eta/rho
1682  // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface
1683  // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel
1684  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1685  }
1686 
1687  return moflux;
1688  }
1689 
1690  AMREX_GPU_DEVICE
1691  AMREX_FORCE_INLINE
1692  amrex::Real
1693  compute_u_flux (const int& i,
1694  const int& j,
1695  const int& k,
1696  const int& icomp,
1697  const amrex::Real& dz,
1698  const amrex::Real& dz1,
1699  const bool& exp_most,
1700  const amrex::Array4<const amrex::Real>& eta_arr,
1701  const amrex::Array4<const amrex::Real>& cons_arr,
1702  const amrex::Array4<const amrex::Real>& velx_arr,
1703  const amrex::Array4<const amrex::Real>& vely_arr,
1704  const amrex::Array4<const amrex::Real>& umm_arr,
1705  const amrex::Array4<const amrex::Real>& /*um_arr*/,
1706  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1707  const amrex::Array4<amrex::Real>& dest_arr) const
1708  {
1709  amrex::Real velx, vely, rho, eta;
1710  int jy, ic, jc;
1711 
1712  int iylo = i <= lbound(vely_arr).x ? lbound(vely_arr).x : i-1;
1713  int iyhi = i > ubound(vely_arr).x ? ubound(vely_arr).x : i;
1714 
1715  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1716  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1717 
1718  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
1719  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1720  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1721  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1722 
1723  velx = velx_arr(i,j,zlo);
1724  vely = 0.25*( vely_arr(iyhi,jy,zlo)+vely_arr(iyhi,jy+1,zlo)
1725  + vely_arr(iylo,jy,zlo)+vely_arr(iylo,jy+1,zlo) );
1726  rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
1727  + cons_arr(ic ,jc,zlo,Rho_comp) );
1728 
1729  amrex::Real Cd = 0.001;
1730  const amrex::Real c = 7e-5;
1731  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1732  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) );
1733  if (wsp_mean <= 5.0) {
1734  Cd = 0.001;
1735  } else if (wsp_mean < 25.0 && wsp_mean > 5.0) {
1736  Cd = 0.001 + c * (wsp_mean - 5.0);
1737  } else {
1738  Cd = 0.0024;
1739  }
1740  amrex::Real deltaz = dz * (zlo - k);
1741 
1742  // NOTE: this is rho*<u'w'> = -K dudz
1743  amrex::Real stressx = -rho * Cd * velx * wsp;
1744 
1745  if (exp_most) {
1746  // surface gradient equal to gradient at first zface
1747  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
1748  dest_arr(i,j,k,icomp) = velx - ugrad * deltaz;
1749  } else {
1750  int ie, je;
1751  ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
1752  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1753  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1754  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1755  eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
1756  + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) );
1757  eta = amrex::max(eta,eta_eps);
1758  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressx/eta*deltaz;
1759  }
1760 
1761  return stressx;
1762  }
1763 
1764  AMREX_GPU_DEVICE
1765  AMREX_FORCE_INLINE
1766  amrex::Real
1767  compute_v_flux (const int& i,
1768  const int& j,
1769  const int& k,
1770  const int& icomp,
1771  const amrex::Real& dz,
1772  const amrex::Real& dz1,
1773  const bool& exp_most,
1774  const amrex::Array4<const amrex::Real>& eta_arr,
1775  const amrex::Array4<const amrex::Real>& cons_arr,
1776  const amrex::Array4<const amrex::Real>& velx_arr,
1777  const amrex::Array4<const amrex::Real>& vely_arr,
1778  const amrex::Array4<const amrex::Real>& umm_arr,
1779  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
1780  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1781  const amrex::Array4<amrex::Real>& dest_arr) const
1782  {
1783  amrex::Real velx, vely, rho, eta;
1784  int ix, ic, jc;
1785 
1786  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
1787  ix = ix > ubound(velx_arr).x ? ubound(velx_arr).x : ix;
1788 
1789  int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
1790  int jxhi = j > ubound(velx_arr).y ? ubound(velx_arr).y : j;
1791 
1792  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1793  jc = j < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
1794  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1795  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1796 
1797  velx = 0.25*( velx_arr(ix,jxhi,zlo)+velx_arr(ix+1,jxhi,zlo)
1798  + velx_arr(ix,jxlo,zlo)+velx_arr(ix+1,jxlo,zlo) );
1799  vely = vely_arr(i,j,zlo);
1800  rho = 0.5*( cons_arr(ic,jc-1,zlo,Rho_comp)
1801  + cons_arr(ic,jc ,zlo,Rho_comp) );
1802 
1803  amrex::Real Cd = 0.001;
1804  const amrex::Real c = 7e-5;
1805  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1806  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) );
1807  if (wsp_mean <= 5.0) {
1808  Cd = 0.001;
1809  } else if (wsp_mean < 25.0 && wsp_mean > 5.0) {
1810  Cd = 0.001 + c * (wsp_mean - 5.0);
1811  } else {
1812  Cd = 0.0024;
1813  }
1814  amrex::Real deltaz = dz * (zlo - k);
1815 
1816  // NOTE: this is rho*<v'w'> = -K dvdz
1817  amrex::Real stressy = -rho * Cd * vely * wsp;
1818 
1819  if (exp_most) {
1820  // surface gradient equal to gradient at first zface
1821  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
1822  dest_arr(i,j,k,icomp) = vely - vgrad * deltaz;
1823  } else {
1824  int ie, je;
1825  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1826  je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
1827  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1828  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1829  eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
1830  + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) );
1831  eta = amrex::max(eta,eta_eps);
1832  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressy/eta*deltaz;
1833  }
1834 
1835  return stressy;
1836  }
1837 
1838 private:
1839  int zlo;
1840  const amrex::Real eta_eps = 1e-8;
1841 };
1842 
1843 
1844 /**
1845  * Custom flux formulation
1846  */
1848 {
1849  custom_flux (int l_zlo)
1850  : zlo(l_zlo) {}
1851 
1852 
1853  AMREX_GPU_DEVICE
1854  AMREX_FORCE_INLINE
1855  amrex::Real
1856  compute_q_flux (const int& i,
1857  const int& j,
1858  const int& k,
1859  const int& n,
1860  const int& icomp,
1861  const amrex::Real& dz,
1862  const amrex::Real& dz1,
1863  const bool& exp_most,
1864  const amrex::Array4<const amrex::Real>& eta_arr,
1865  const amrex::Array4<const amrex::Real>& cons_arr,
1866  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1867  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1868  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1869  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
1870  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1871  const amrex::Array4<const amrex::Real>& q_star_arr,
1872  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1873  const amrex::Array4<amrex::Real>& dest_arr) const
1874  {
1875  amrex::Real rho, eta;
1876 
1877  int ic, jc;
1878  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1879  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1880  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1881  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1882 
1883  rho = cons_arr(ic,jc,zlo,Rho_comp);
1884 
1885  amrex::Real qstar = q_star_arr(ic,jc,zlo);
1886  amrex::Real deltaz = dz * (zlo - k);
1887 
1888  // NOTE: this is rho*<q'w'> = -K dqdz
1889  amrex::Real moflux = (std::abs(qstar) > eps) ? rho * qstar : 0.0;
1890 
1891  if (exp_most) {
1892  // surface gradient equal to gradient at first zface
1893  amrex::Real rqvgrad = ( cons_arr(ic,jc,zlo+1,RhoQ1_comp) - cons_arr(ic,jc,zlo ,RhoQ1_comp) ) / (0.5*(dz+dz1));
1894  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
1895  } else {
1896  int ie, je;
1897  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1898  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1899  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1900  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1901  eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
1902  eta = amrex::max(eta,eta_eps);
1903  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1904  }
1905 
1906  return moflux;
1907  }
1908 
1909  AMREX_GPU_DEVICE
1910  AMREX_FORCE_INLINE
1911  amrex::Real
1912  compute_t_flux (const int& i,
1913  const int& j,
1914  const int& k,
1915  const int& n,
1916  const int& icomp,
1917  const amrex::Real& dz,
1918  const amrex::Real& dz1,
1919  const bool& exp_most,
1920  const amrex::Array4<const amrex::Real>& eta_arr,
1921  const amrex::Array4<const amrex::Real>& cons_arr,
1922  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1923  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1924  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1925  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
1926  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1927  const amrex::Array4<const amrex::Real>& t_star_arr,
1928  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1929  const amrex::Array4<amrex::Real>& dest_arr) const
1930  {
1931  amrex::Real rho, eta;
1932 
1933  int ic, jc;
1934  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1935  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1936  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1937  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1938 
1939  rho = cons_arr(ic,jc,zlo,Rho_comp);
1940 
1941  amrex::Real tstar = t_star_arr(ic,jc,zlo);
1942  amrex::Real deltaz = dz * (zlo - k);
1943 
1944  // NOTE: this is rho*<T'w'> = -K dTdz
1945  amrex::Real moflux = (std::abs(tstar) > eps) ? rho * tstar : 0.0;
1946 
1947  if (exp_most) {
1948  // surface gradient equal to gradient at first zface
1949  amrex::Real rthetagrad = ( cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo ,RhoTheta_comp) ) / (0.5*(dz+dz1));
1950  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
1951  } else {
1952  int ie, je;
1953  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1954  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1955  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1956  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1957  eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
1958  eta = amrex::max(eta,eta_eps);
1959  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1960  }
1961 
1962  return moflux;
1963  }
1964 
1965  AMREX_GPU_DEVICE
1966  AMREX_FORCE_INLINE
1967  amrex::Real
1968  compute_u_flux (const int& i,
1969  const int& j,
1970  const int& k,
1971  const int& icomp,
1972  const amrex::Real& dz,
1973  const amrex::Real& dz1,
1974  const bool& exp_most,
1975  const amrex::Array4<const amrex::Real>& eta_arr,
1976  const amrex::Array4<const amrex::Real>& cons_arr,
1977  const amrex::Array4<const amrex::Real>& velx_arr,
1978  const amrex::Array4<const amrex::Real>& vely_arr,
1979  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1980  const amrex::Array4<const amrex::Real>& /*um_arr*/,
1981  const amrex::Array4<const amrex::Real>& u_star_arr,
1982  const amrex::Array4<amrex::Real>& dest_arr) const
1983  {
1984  amrex::Real velx, vely, rho, eta;
1985  int jy, ic, jc;
1986 
1987  int iylo = i <= lbound(vely_arr).x ? lbound(vely_arr).x : i-1;
1988  int iyhi = i > ubound(vely_arr).x ? ubound(vely_arr).x : i;
1989 
1990  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1991  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1992 
1993  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
1994  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1995  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1996  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1997 
1998  velx = velx_arr(i,j,zlo);
1999  vely = 0.25*( vely_arr(iyhi,jy,zlo)+vely_arr(iyhi,jy+1,zlo)
2000  + vely_arr(iylo,jy,zlo)+vely_arr(iylo,jy+1,zlo) );
2001  rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
2002  + cons_arr(ic ,jc,zlo,Rho_comp) );
2003 
2004  amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
2005  amrex::Real wsp = sqrt(velx*velx+vely*vely);
2006  amrex::Real deltaz = dz * (zlo - k);
2007 
2008  // NOTE: this is rho*<u'w'> = -K dudz
2009  amrex::Real stressx = -rho * ustar * ustar * velx / wsp;
2010 
2011  if (exp_most) {
2012  // surface gradient equal to gradient at first zface
2013  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
2014  dest_arr(i,j,k,icomp) = velx - ugrad * deltaz;
2015  } else {
2016  int ie, je;
2017  ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
2018  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
2019  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
2020  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
2021  eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
2022  + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) );
2023  eta = amrex::max(eta,eta_eps);
2024  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressx/eta*deltaz;
2025  }
2026 
2027  return stressx;
2028  }
2029 
2030  AMREX_GPU_DEVICE
2031  AMREX_FORCE_INLINE
2032  amrex::Real
2033  compute_v_flux (const int& i,
2034  const int& j,
2035  const int& k,
2036  const int& icomp,
2037  const amrex::Real& dz,
2038  const amrex::Real& dz1,
2039  const bool& exp_most,
2040  const amrex::Array4<const amrex::Real>& eta_arr,
2041  const amrex::Array4<const amrex::Real>& cons_arr,
2042  const amrex::Array4<const amrex::Real>& velx_arr,
2043  const amrex::Array4<const amrex::Real>& vely_arr,
2044  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2045  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2046  const amrex::Array4<const amrex::Real>& u_star_arr,
2047  const amrex::Array4<amrex::Real>& dest_arr) const
2048  {
2049  amrex::Real velx, vely, rho, eta;
2050  int ix, ic, jc;
2051 
2052  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
2053  ix = ix > ubound(velx_arr).x ? ubound(velx_arr).x : ix;
2054 
2055  int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
2056  int jxhi = j > ubound(velx_arr).y ? ubound(velx_arr).y : j;
2057 
2058  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
2059  jc = j < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
2060  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2061  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2062 
2063  velx = 0.25*( velx_arr(ix,jxhi,zlo)+velx_arr(ix+1,jxhi,zlo)
2064  + velx_arr(ix,jxlo,zlo)+velx_arr(ix+1,jxlo,zlo) );
2065  vely = vely_arr(i,j,zlo);
2066  rho = 0.5*( cons_arr(ic,jc-1,zlo,Rho_comp)
2067  + cons_arr(ic,jc ,zlo,Rho_comp) );
2068 
2069  amrex::Real ustar = 0.5 * ( u_star_arr(ic,jc-1,zlo) + u_star_arr(ic,jc,zlo) );
2070  amrex::Real wsp = sqrt(velx*velx+vely*vely);
2071  amrex::Real deltaz = dz * (zlo - k);
2072 
2073  // NOTE: this is rho*<v'w'> = -K dvdz
2074  amrex::Real stressy = -rho * ustar * ustar * vely / wsp;
2075 
2076  if (exp_most) {
2077  // surface gradient equal to gradient at first zface
2078  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
2079  dest_arr(i,j,k,icomp) = vely - vgrad * deltaz;
2080  } else {
2081  int ie, je;
2082  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
2083  je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
2084  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
2085  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
2086  eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
2087  + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) );
2088  eta = amrex::max(eta,eta_eps);
2089  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressy/eta*deltaz;
2090  }
2091 
2092  return stressy;
2093  }
2094 
2095 private:
2096  int zlo;
2097  const amrex::Real eps = 1e-15;
2098  const amrex::Real eta_eps = 1e-8;
2099 };
2100 
2101 /**
2102  * Rotate flux formulation
2103  */
2105 {
2106  rotate_flux (int l_zlo)
2107  : zlo(l_zlo) {}
2108 
2109 
2110  AMREX_GPU_DEVICE
2111  AMREX_FORCE_INLINE
2112  amrex::Real
2113  compute_q_flux (const int& i,
2114  const int& j,
2115  const int& k,
2116  const int& n,
2117  const int& icomp,
2118  const amrex::Real& dz,
2119  const amrex::Real& dz1,
2120  const bool& /*exp_most*/,
2121  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2122  const amrex::Array4<const amrex::Real>& cons_arr,
2123  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2124  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2125  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2126  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
2127  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2128  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
2129  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
2130  const amrex::Array4<amrex::Real>& dest_arr) const
2131  {
2132  int ic, jc;
2133  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
2134  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
2135  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2136  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2137 
2138  // surface gradient equal to gradient at first zface
2139  amrex::Real deltaz = dz * (zlo - k);
2140  amrex::Real rqvgrad = ( cons_arr(ic,jc,zlo+1,RhoQ1_comp)
2141  - cons_arr(ic,jc,zlo ,RhoQ1_comp) ) / (0.5*(dz+dz1));
2142  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
2143 
2144  // NOTE: this is the total stress
2145  amrex::Real rho = cons_arr(ic,jc,zlo,Rho_comp);
2146  amrex::Real qstar = 0.0; //q_star_arr(ic,jc,zlo);
2147  amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0;
2148 
2149  return moflux;
2150  }
2151 
2152  AMREX_GPU_DEVICE
2153  AMREX_FORCE_INLINE
2154  amrex::Real
2155  compute_t_flux (const int& i,
2156  const int& j,
2157  const int& k,
2158  const int& n,
2159  const int& icomp,
2160  const amrex::Real& dz,
2161  const amrex::Real& dz1,
2162  const bool& /*exp_most*/,
2163  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2164  const amrex::Array4<const amrex::Real>& cons_arr,
2165  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2166  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2167  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2168  const amrex::Array4<const amrex::Real>& tm_arr,
2169  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2170  const amrex::Array4<const amrex::Real>& t_star_arr,
2171  const amrex::Array4<const amrex::Real>& t_surf_arr,
2172  const amrex::Array4<amrex::Real>& dest_arr) const
2173  {
2174  int ic, jc;
2175  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
2176  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
2177  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2178  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2179 
2180  // surface gradient equal to gradient at first zface
2181  amrex::Real deltaz = dz * (zlo - k);
2182  amrex::Real rthetagrad = ( cons_arr(ic,jc,zlo+1,RhoTheta_comp)
2183  - cons_arr(ic,jc,zlo ,RhoTheta_comp) ) / (0.5*(dz+dz1));
2184  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
2185 
2186  // NOTE: this is the total stress
2187  amrex::Real theta_mean = tm_arr(ic,jc,zlo);
2188  amrex::Real theta_surf = t_surf_arr(ic,jc,zlo);
2189  amrex::Real rho = cons_arr(ic,jc,zlo,Rho_comp);
2190  amrex::Real tstar = t_star_arr(ic,jc,zlo);
2191  amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar * (theta_mean-theta_surf) : 0.0;
2192 
2193  return moflux;
2194  }
2195 
2196  AMREX_GPU_DEVICE
2197  AMREX_FORCE_INLINE
2198  amrex::Real
2199  compute_u_flux (const int& i,
2200  const int& j,
2201  const int& k,
2202  const int& icomp,
2203  const amrex::Real& dz,
2204  const amrex::Real& dz1,
2205  const bool& /*exp_most*/,
2206  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2207  const amrex::Array4<const amrex::Real>& cons_arr,
2208  const amrex::Array4<const amrex::Real>& velx_arr,
2209  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2210  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2211  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2212  const amrex::Array4<const amrex::Real>& u_star_arr,
2213  const amrex::Array4<amrex::Real>& dest_arr) const
2214  {
2215  int ic, jc;
2216  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
2217  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
2218  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2219  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2220 
2221  // surface gradient equal to gradient at first zface
2222  amrex::Real deltaz = dz * (zlo - k);
2223  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx_arr(i,j,zlo)) / (0.5*(dz+dz1));
2224  dest_arr(i,j,k,icomp) = velx_arr(i,j,zlo) - ugrad * deltaz;
2225 
2226  // NOTE: this is the total stress
2227  amrex::Real rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
2228  + cons_arr(ic ,jc,zlo,Rho_comp) );
2229  amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
2230  amrex::Real stressx = -rho * ustar * ustar;
2231 
2232  return stressx;
2233  }
2234 
2235  AMREX_GPU_DEVICE
2236  AMREX_FORCE_INLINE
2237  amrex::Real
2238  compute_v_flux (const int& i,
2239  const int& j,
2240  const int& k,
2241  const int& icomp,
2242  const amrex::Real& dz,
2243  const amrex::Real& dz1,
2244  const bool& /*exp_most*/,
2245  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2246  const amrex::Array4<const amrex::Real>& /*cons_arr*/,
2247  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2248  const amrex::Array4<const amrex::Real>& vely_arr,
2249  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2250  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2251  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2252  const amrex::Array4<amrex::Real>& dest_arr) const
2253  {
2254  // surface gradient equal to gradient at first zface
2255  amrex::Real deltaz = dz * (zlo - k);
2256  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely_arr(i,j,zlo)) / (0.5*(dz+dz1));
2257  dest_arr(i,j,k,icomp) = vely_arr(i,j,zlo) - vgrad * deltaz;
2258 
2259  // NOTE: this is the total stress
2260  amrex::Real stressy = 0.0;
2261 
2262  return stressy;
2263  }
2264 
2265 private:
2266  int zlo;
2267  const amrex::Real eps = 1e-15;
2268 };
2269 #endif
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:7
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Donelan_roughness(amrex::Real ustar)
Definition: ERF_MOSTRoughness.H:29
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real COARE3_roughness(amrex::Real zref, amrex::Real umm, amrex::Real ustar)
Definition: ERF_MOSTRoughness.H:9
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real air_viscosity(amrex::Real T_degK)
Definition: ERF_MOSTStress.H:79
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=0.0)
Definition: ERF_Wstar.H:13
@ Theta_v
Definition: ERF_IndexDefines.H:157
@ Q_v
Definition: ERF_IndexDefines.H:160
@ Mom_v
Definition: ERF_IndexDefines.H:156
@ theta
Definition: ERF_MM5.H:20
@ rho
Definition: ERF_Kessler.H:30
Definition: ERF_MOSTStress.H:136
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:199
adiabatic_charnock(amrex::Real zref, amrex::Real flux, amrex::Real cnk_a, bool cnk_visc)
Definition: ERF_MOSTStress.H:137
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< 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 > &, 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 > &, 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_MOSTStress.H:151
most_data mdata
Definition: ERF_MOSTStress.H:196
const amrex::Real tol
Definition: ERF_MOSTStress.H:198
similarity_funs sfuns
Definition: ERF_MOSTStress.H:197
Definition: ERF_MOSTStress.H:270
similarity_funs sfuns
Definition: ERF_MOSTStress.H:320
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:322
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< 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 > &, 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_MOSTStress.H:281
adiabatic_donelan(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:271
most_data mdata
Definition: ERF_MOSTStress.H:319
const amrex::Real tol
Definition: ERF_MOSTStress.H:321
Definition: ERF_MOSTStress.H:207
adiabatic_mod_charnock(amrex::Real zref, amrex::Real flux, amrex::Real depth)
Definition: ERF_MOSTStress.H:208
similarity_funs sfuns
Definition: ERF_MOSTStress.H:260
const amrex::Real tol
Definition: ERF_MOSTStress.H:261
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:262
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< 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 > &, 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_MOSTStress.H:221
most_data mdata
Definition: ERF_MOSTStress.H:259
Definition: ERF_MOSTStress.H:330
const amrex::Real eps
Definition: ERF_MOSTStress.H:388
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< 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 > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &Hwave_arr, const amrex::Array4< amrex::Real > &Lwave_arr, const amrex::Array4< amrex::Real > &eta_arr) const
Definition: ERF_MOSTStress.H:341
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:389
most_data mdata
Definition: ERF_MOSTStress.H:385
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:391
adiabatic_wave_coupled(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:331
similarity_funs sfuns
Definition: ERF_MOSTStress.H:386
const amrex::Real tol
Definition: ERF_MOSTStress.H:387
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:390
Definition: ERF_MOSTStress.H:90
adiabatic(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:91
similarity_funs sfuns
Definition: ERF_MOSTStress.H:128
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 > &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 > &, 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_MOSTStress.H:101
most_data mdata
Definition: ERF_MOSTStress.H:127
Definition: ERF_MOSTStress.H:1848
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_v_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:2033
int zlo
Definition: ERF_MOSTStress.H:2096
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_u_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1968
const amrex::Real eta_eps
Definition: ERF_MOSTStress.H:2098
custom_flux(int l_zlo)
Definition: ERF_MOSTStress.H:1849
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_t_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &t_star_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1912
const amrex::Real eps
Definition: ERF_MOSTStress.H:2097
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_q_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &q_star_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1856
Definition: ERF_MOSTStress.H:1567
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_u_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1693
donelan_flux(int l_zlo)
Definition: ERF_MOSTStress.H:1568
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_t_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &t_surf_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1631
const amrex::Real eta_eps
Definition: ERF_MOSTStress.H:1840
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_v_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1767
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_q_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1575
int zlo
Definition: ERF_MOSTStress.H:1839
Definition: ERF_MOSTStress.H:1259
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_u_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &um_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1406
moeng_flux(int l_zlo)
Definition: ERF_MOSTStress.H:1260
int zlo
Definition: ERF_MOSTStress.H:1556
const amrex::Real eps
Definition: ERF_MOSTStress.H:1557
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_q_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1267
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_v_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &vm_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1482
const amrex::Real eta_eps
Definition: ERF_MOSTStress.H:1558
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1559
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_t_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &exp_most, const amrex::Array4< const amrex::Real > &eta_arr, 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 > &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< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:1323
Definition: ERF_MOSTStress.H:13
amrex::Real zref
Reference height (m)
Definition: ERF_MOSTStress.H:15
amrex::Real Cnk_b2
Modified Charnock Eq (4) https://doi.org/10.1175/JAMC-D-17-0137.1.
Definition: ERF_MOSTStress.H:23
amrex::Real Cnk_b
Definition: ERF_MOSTStress.H:25
amrex::Real Cnk_d
Modified Charnock Eq (4) https://doi.org/10.1175/JAMC-D-17-0137.1.
Definition: ERF_MOSTStress.H:24
amrex::Real kappa
von Karman constant
Definition: ERF_MOSTStress.H:17
amrex::Real gravity
Acceleration due to gravity (m/s^2)
Definition: ERF_MOSTStress.H:18
amrex::Real Cnk_a
Standard Charnock constant https://doi.org/10.1175/JAMC-D-17-0137.1.
Definition: ERF_MOSTStress.H:21
const amrex::Real Bjr_beta
Definition: ERF_MOSTStress.H:28
amrex::Real Cnk_b1
Modified Charnock Eq (4) https://doi.org/10.1175/JAMC-D-17-0137.1.
Definition: ERF_MOSTStress.H:22
amrex::Real z0_const
Roughness height – default constant value(m)
Definition: ERF_MOSTStress.H:16
bool visc
Use viscous Charnock formulation.
Definition: ERF_MOSTStress.H:26
amrex::Real surf_temp_flux
Heat flux.
Definition: ERF_MOSTStress.H:19
Definition: ERF_MOSTStress.H:2105
int zlo
Definition: ERF_MOSTStress.H:2266
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_q_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:2113
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_t_flux(const int &i, const int &j, const int &k, const int &n, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &t_star_arr, const amrex::Array4< const amrex::Real > &t_surf_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:2155
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_u_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &cons_arr, const amrex::Array4< const amrex::Real > &velx_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:2199
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_v_flux(const int &i, const int &j, const int &k, const int &icomp, const amrex::Real &dz, const amrex::Real &dz1, const bool &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &vely_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &dest_arr) const
Definition: ERF_MOSTStress.H:2238
const amrex::Real eps
Definition: ERF_MOSTStress.H:2267
rotate_flux(int l_zlo)
Definition: ERF_MOSTStress.H:2106
Definition: ERF_MOSTStress.H:36
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:40
amrex::Real beta_m
Constants from Dyer, BLM, 1974.
Definition: ERF_MOSTStress.H:65
amrex::Real beta_h
https://doi.org/10.1007/BF00240838
Definition: ERF_MOSTStress.H:66
amrex::Real gamma_h
Definition: ERF_MOSTStress.H:68
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:54
amrex::Real gamma_m
Definition: ERF_MOSTStress.H:67
Definition: ERF_MOSTStress.H:477
const amrex::Real tol
Definition: ERF_MOSTStress.H:560
similarity_funs sfuns
Definition: ERF_MOSTStress.H:559
most_data mdata
Definition: ERF_MOSTStress.H:558
surface_flux_charnock(amrex::Real zref, amrex::Real flux, amrex::Real cnk_a, bool cnk_visc)
Definition: ERF_MOSTStress.H:478
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< 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 > &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_MOSTStress.H:492
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:561
Definition: ERF_MOSTStress.H:653
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< 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 > &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_MOSTStress.H:664
const amrex::Real tol
Definition: ERF_MOSTStress.H:725
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:726
most_data mdata
Definition: ERF_MOSTStress.H:723
similarity_funs sfuns
Definition: ERF_MOSTStress.H:724
surface_flux_donelan(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:654
Definition: ERF_MOSTStress.H:569
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:645
most_data mdata
Definition: ERF_MOSTStress.H:642
surface_flux_mod_charnock(amrex::Real zref, amrex::Real flux, amrex::Real depth)
Definition: ERF_MOSTStress.H:570
const amrex::Real tol
Definition: ERF_MOSTStress.H:644
similarity_funs sfuns
Definition: ERF_MOSTStress.H:643
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< 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 > &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_MOSTStress.H:583
Definition: ERF_MOSTStress.H:734
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:816
surface_flux_wave_coupled(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:735
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< 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 > &olen_arr, const amrex::Array4< amrex::Real > &pblh_arr, const amrex::Array4< amrex::Real > &Hwave_arr, const amrex::Array4< amrex::Real > &Lwave_arr, const amrex::Array4< amrex::Real > &eta_arr) const
Definition: ERF_MOSTStress.H:745
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:815
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:814
most_data mdata
Definition: ERF_MOSTStress.H:810
similarity_funs sfuns
Definition: ERF_MOSTStress.H:811
const amrex::Real tol
Definition: ERF_MOSTStress.H:812
const amrex::Real eps
Definition: ERF_MOSTStress.H:813
Definition: ERF_MOSTStress.H:399
similarity_funs sfuns
Definition: ERF_MOSTStress.H:467
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 > &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 > &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_MOSTStress.H:410
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:469
most_data mdata
Definition: ERF_MOSTStress.H:466
surface_flux(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:400
const amrex::Real tol
Definition: ERF_MOSTStress.H:468
Definition: ERF_MOSTStress.H:904
most_data mdata
Definition: ERF_MOSTStress.H:987
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< 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 > &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_MOSTStress.H:919
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:990
const amrex::Real tol
Definition: ERF_MOSTStress.H:989
surface_temp_charnock(amrex::Real zref, amrex::Real flux, amrex::Real cnk_a, bool cnk_visc)
Definition: ERF_MOSTStress.H:905
similarity_funs sfuns
Definition: ERF_MOSTStress.H:988
Definition: ERF_MOSTStress.H:1084
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1157
most_data mdata
Definition: ERF_MOSTStress.H:1156
surface_temp_donelan(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:1085
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< 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 > &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_MOSTStress.H:1095
const amrex::Real tol
Definition: ERF_MOSTStress.H:1158
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1159
Definition: ERF_MOSTStress.H:998
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1074
surface_temp_mod_charnock(amrex::Real zref, amrex::Real flux, amrex::Real depth)
Definition: ERF_MOSTStress.H:999
const amrex::Real tol
Definition: ERF_MOSTStress.H:1075
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1076
most_data mdata
Definition: ERF_MOSTStress.H:1073
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< 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 > &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_MOSTStress.H:1012
Definition: ERF_MOSTStress.H:1167
const amrex::Real eps
Definition: ERF_MOSTStress.H:1248
const amrex::Real tol
Definition: ERF_MOSTStress.H:1247
most_data mdata
Definition: ERF_MOSTStress.H:1245
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< 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 > &olen_arr, const amrex::Array4< amrex::Real > &pblh_arr, const amrex::Array4< amrex::Real > &Hwave_arr, const amrex::Array4< amrex::Real > &Lwave_arr, const amrex::Array4< amrex::Real > &eta_arr) const
Definition: ERF_MOSTStress.H:1178
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:1249
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1246
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:1250
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1251
surface_temp_wave_coupled(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:1168
Definition: ERF_MOSTStress.H:824
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 > &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 > &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_MOSTStress.H:835
similarity_funs sfuns
Definition: ERF_MOSTStress.H:894
surface_temp(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:825
const amrex::Real tol
Definition: ERF_MOSTStress.H:895
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:896
most_data mdata
Definition: ERF_MOSTStress.H:893