ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
440  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
441  do {
442  ustar = u_star_arr(i,j,k);
443  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);
444  if (w_star_arr) {
445  // update w* and Umagmean
446  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
447  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
448  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
449  umm = std::max(umm, WSMIN);
450  }
451  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
452  zeta = mdata.zref / Olen;
453  psi_m = sfuns.calc_psi_m(zeta);
454  psi_h = sfuns.calc_psi_h(zeta);
455  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0_arr(i,j,k)) - psi_m);
456  ++iter;
457  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
458  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
459 
460  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h) /
461  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
462  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
463  olen_arr(i,j,k) = Olen;
464  }
465 
466 private:
469  const amrex::Real tol = 1.0e-5;
470  const amrex::Real WSMIN = 0.1; // minimum wind speed
471 };
472 
473 
474 /**
475  * Surface flux with charnock roughness
476  */
478 {
479  surface_flux_charnock (amrex::Real zref,
480  amrex::Real flux,
481  amrex::Real cnk_a,
482  bool cnk_visc)
483  {
484  mdata.zref = zref;
485  mdata.surf_temp_flux = flux;
486  mdata.Cnk_a = cnk_a;
487  mdata.visc = cnk_visc;
488  }
489 
490  AMREX_GPU_DEVICE
491  AMREX_FORCE_INLINE
492  void
493  iterate_flux (const int& i,
494  const int& j,
495  const int& k,
496  const int& max_iters,
497  const amrex::Array4<amrex::Real>& z0_arr,
498  const amrex::Array4<const amrex::Real>& umm_arr,
499  const amrex::Array4<const amrex::Real>& tm_arr,
500  const amrex::Array4<const amrex::Real>& tvm_arr,
501  const amrex::Array4<const amrex::Real>& qvm_arr,
502  const amrex::Array4<amrex::Real>& u_star_arr,
503  const amrex::Array4<amrex::Real>& w_star_arr,
504  const amrex::Array4<amrex::Real>& t_star_arr,
505  const amrex::Array4<amrex::Real>& q_star_arr,
506  const amrex::Array4<amrex::Real>& t_surf_arr,
507  const amrex::Array4<amrex::Real>& olen_arr,
508  const amrex::Array4<amrex::Real>& pblh_arr,
509  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
510  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
511  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
512  {
513  int iter = 0;
514  amrex::Real ustar = 0.0;
515  amrex::Real wstar = 0.0;
516  amrex::Real tflux = 0.0;
517  amrex::Real z0 = 0.0;
518  amrex::Real zeta = 0.0;
519  amrex::Real psi_m = 0.0;
520  amrex::Real psi_h = 0.0;
521  amrex::Real Olen = 0.0;
522  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
523  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
524  do {
525  ustar = u_star_arr(i,j,k);
526  if (mdata.Cnk_a > 0) {
527  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
528  if (mdata.visc) {
529  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, 0.05);
530  }
531  } else {
532  z0 = COARE3_roughness(mdata.zref, umm, ustar);
533  }
534  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);
535  if (w_star_arr) {
536  // update w* and Umagmean
537  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
538  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
539  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
540  umm = std::max(umm, WSMIN);
541  }
542  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
543  zeta = mdata.zref / Olen;
544  psi_m = sfuns.calc_psi_m(zeta);
545  psi_h = sfuns.calc_psi_h(zeta);
546  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
547  ++iter;
548  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
549  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
550 
551  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
552  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
553  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
554  olen_arr(i,j,k) = Olen;
555  z0_arr(i,j,k) = z0;
556  }
557 
558 private:
561  const amrex::Real tol = 1.0e-5;
562  const amrex::Real WSMIN = 0.1; // minimum wind speed
563 };
564 
565 
566 /**
567  * Surface flux with modified charnock roughness
568  */
570 {
571  surface_flux_mod_charnock (amrex::Real zref,
572  amrex::Real flux,
573  amrex::Real depth)
574  {
575  mdata.zref = zref;
576  mdata.surf_temp_flux = flux;
577  mdata.Cnk_d = depth;
578  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
579  }
580 
581  AMREX_GPU_DEVICE
582  AMREX_FORCE_INLINE
583  void
584  iterate_flux (const int& i,
585  const int& j,
586  const int& k,
587  const int& max_iters,
588  const amrex::Array4<amrex::Real>& z0_arr,
589  const amrex::Array4<const amrex::Real>& umm_arr,
590  const amrex::Array4<const amrex::Real>& tm_arr,
591  const amrex::Array4<const amrex::Real>& tvm_arr,
592  const amrex::Array4<const amrex::Real>& qvm_arr,
593  const amrex::Array4<amrex::Real>& u_star_arr,
594  const amrex::Array4<amrex::Real>& w_star_arr,
595  const amrex::Array4<amrex::Real>& t_star_arr,
596  const amrex::Array4<amrex::Real>& q_star_arr,
597  const amrex::Array4<amrex::Real>& t_surf_arr,
598  const amrex::Array4<amrex::Real>& olen_arr,
599  const amrex::Array4<amrex::Real>& pblh_arr,
600  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
601  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
602  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
603  {
604  int iter = 0;
605  amrex::Real ustar = 0.0;
606  amrex::Real wstar = 0.0;
607  amrex::Real tflux = 0.0;
608  amrex::Real z0 = 0.0;
609  amrex::Real zeta = 0.0;
610  amrex::Real psi_m = 0.0;
611  amrex::Real psi_h = 0.0;
612  amrex::Real Olen = 0.0;
613  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
614  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
615  do {
616  ustar = u_star_arr(i,j,k);
617  z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) );
618  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);
619  if (w_star_arr) {
620  // update w* and Umagmean
621  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
622  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
623  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
624  umm = std::max(umm, WSMIN);
625  }
626  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
627  zeta = mdata.zref / Olen;
628  psi_m = sfuns.calc_psi_m(zeta);
629  psi_h = sfuns.calc_psi_h(zeta);
630  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
631  ++iter;
632  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
633  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
634 
635  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
636  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
637  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
638  olen_arr(i,j,k) = Olen;
639  z0_arr(i,j,k) = z0;
640  }
641 
642 private:
645  const amrex::Real tol = 1.0e-5;
646  const amrex::Real WSMIN = 0.1; // minimum wind speed
647 };
648 
649 
650 /**
651  * Surface flux with donelan roughness
652  */
654 {
655  surface_flux_donelan (amrex::Real zref,
656  amrex::Real flux)
657  {
658  mdata.zref = zref;
659  mdata.surf_temp_flux = flux;
660  }
661 
662  AMREX_GPU_DEVICE
663  AMREX_FORCE_INLINE
664  void
665  iterate_flux (const int& i,
666  const int& j,
667  const int& k,
668  const int& max_iters,
669  const amrex::Array4<amrex::Real>& z0_arr,
670  const amrex::Array4<const amrex::Real>& umm_arr,
671  const amrex::Array4<const amrex::Real>& tm_arr,
672  const amrex::Array4<const amrex::Real>& tvm_arr,
673  const amrex::Array4<const amrex::Real>& qvm_arr,
674  const amrex::Array4<amrex::Real>& u_star_arr,
675  const amrex::Array4<amrex::Real>& w_star_arr,
676  const amrex::Array4<amrex::Real>& t_star_arr,
677  const amrex::Array4<amrex::Real>& q_star_arr,
678  const amrex::Array4<amrex::Real>& t_surf_arr,
679  const amrex::Array4<amrex::Real>& olen_arr,
680  const amrex::Array4<amrex::Real>& pblh_arr,
681  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
682  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
683  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
684  {
685  int iter = 0;
686  amrex::Real ustar = 0.0;
687  amrex::Real wstar = 0.0;
688  amrex::Real tflux = 0.0;
689  amrex::Real z0 = 0.0;
690  amrex::Real zeta = 0.0;
691  amrex::Real psi_m = 0.0;
692  amrex::Real psi_h = 0.0;
693  amrex::Real Olen = 0.0;
694  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
695  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
696  do {
697  ustar = u_star_arr(i,j,k);
698  z0 = Donelan_roughness(ustar);
699  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);
700  if (w_star_arr) {
701  // update w* and Umagmean
702  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
703  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
704  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
705  umm = std::max(umm, WSMIN);
706  }
707  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
708  zeta = mdata.zref / Olen;
709  psi_m = sfuns.calc_psi_m(zeta);
710  psi_h = sfuns.calc_psi_h(zeta);
711  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
712  ++iter;
713  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
714  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
715 
716  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
717  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
718  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
719  olen_arr(i,j,k) = Olen;
720  z0_arr(i,j,k) = z0;
721  }
722 
723 private:
726  const amrex::Real tol = 1.0e-5;
727  const amrex::Real WSMIN = 0.1; // minimum wind speed
728 };
729 
730 
731 /**
732  * Surface flux with wave-coupled roughness
733  */
735 {
736  surface_flux_wave_coupled (amrex::Real zref,
737  amrex::Real flux)
738  {
739  mdata.zref = zref;
740  mdata.surf_temp_flux = flux;
741  }
742 
743  AMREX_GPU_DEVICE
744  AMREX_FORCE_INLINE
745  void
746  iterate_flux (const int& i,
747  const int& j,
748  const int& k,
749  const int& max_iters,
750  const amrex::Array4<amrex::Real>& z0_arr,
751  const amrex::Array4<const amrex::Real>& umm_arr,
752  const amrex::Array4<const amrex::Real>& tm_arr,
753  const amrex::Array4<const amrex::Real>& tvm_arr,
754  const amrex::Array4<const amrex::Real>& qvm_arr,
755  const amrex::Array4<amrex::Real>& u_star_arr,
756  const amrex::Array4<amrex::Real>& w_star_arr,
757  const amrex::Array4<amrex::Real>& t_star_arr,
758  const amrex::Array4<amrex::Real>& q_star_arr,
759  const amrex::Array4<amrex::Real>& t_surf_arr,
760  const amrex::Array4<amrex::Real>& olen_arr,
761  const amrex::Array4<amrex::Real>& pblh_arr,
762  const amrex::Array4<amrex::Real>& Hwave_arr,
763  const amrex::Array4<amrex::Real>& Lwave_arr,
764  const amrex::Array4<amrex::Real>& eta_arr) const
765  {
766  int iter = 0;
767  amrex::Real ustar = 0.0;
768  amrex::Real wstar = 0.0;
769  amrex::Real tflux = 0.0;
770  amrex::Real z0 = 0.0;
771  amrex::Real zeta = 0.0;
772  amrex::Real psi_m = 0.0;
773  amrex::Real psi_h = 0.0;
774  amrex::Real Olen = 0.0;
775  int ie, je;
776  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
777  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
778  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
779  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
780  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
781  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
782  do {
783  ustar = u_star_arr(i,j,k);
784  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 )
785  + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
786  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);
787  if (w_star_arr) {
788  // update w* and Umagmean
789  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
790  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
791  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
792  umm = std::max(umm, WSMIN);
793  }
794  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
795  zeta = mdata.zref / Olen;
796  psi_m = sfuns.calc_psi_m(zeta);
797  psi_h = sfuns.calc_psi_h(zeta);
798  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
799  ++iter;
800  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
801  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
802 
803  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) /
804  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
805  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
806  olen_arr(i,j,k) = Olen;
807  z0_arr(i,j,k) = z0;
808  }
809 
810 private:
813  const amrex::Real tol = 1.0e-5;
814  const amrex::Real eps = 1e-15;
815  const amrex::Real z0_eps = 1.0e-6;
816  const amrex::Real z0_max = 0.1;
817  const amrex::Real WSMIN = 0.1; // minimum wind speed
818 };
819 
820 
821 /**
822  * Surface temperature with constant roughness
823  */
825 {
826  surface_temp (amrex::Real zref,
827  amrex::Real flux)
828  {
829  mdata.zref = zref;
830  mdata.surf_temp_flux = flux;
831  }
832 
833  AMREX_GPU_DEVICE
834  AMREX_FORCE_INLINE
835  void
836  iterate_flux (const int& i,
837  const int& j,
838  const int& k,
839  const int& max_iters,
840  const amrex::Array4<const amrex::Real>& z0_arr,
841  const amrex::Array4<const amrex::Real>& umm_arr,
842  const amrex::Array4<const amrex::Real>& tm_arr,
843  const amrex::Array4<const amrex::Real>& tvm_arr,
844  const amrex::Array4<const amrex::Real>& qvm_arr,
845  const amrex::Array4<amrex::Real>& u_star_arr,
846  const amrex::Array4<amrex::Real>& w_star_arr,
847  const amrex::Array4<amrex::Real>& t_star_arr,
848  const amrex::Array4<amrex::Real>& q_star_arr,
849  const amrex::Array4<amrex::Real>& t_surf_arr,
850  const amrex::Array4<amrex::Real>& olen_arr,
851  const amrex::Array4<amrex::Real>& pblh_arr,
852  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
853  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
854  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
855  {
856  int iter = 0;
857  amrex::Real ustar = 0.0;
858  amrex::Real wstar = 0.0;
859  amrex::Real tflux = 0.0;
860  amrex::Real zeta = 0.0;
861  amrex::Real psi_m = 0.0;
862  amrex::Real psi_h = 0.0;
863  amrex::Real Olen = 0.0;
864  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
865  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
866  do {
867  ustar = u_star_arr(i,j,k);
868  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
869  (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h); // <w'T'>
870  tflux *= (1 + 0.61*qvm_arr(i,j,k));
871  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
872  if (w_star_arr) {
873  // update w* and Umagmean
874  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
875  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
876  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
877  umm = std::max(umm, WSMIN);
878  }
879  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
880  zeta = mdata.zref / Olen;
881  psi_m = sfuns.calc_psi_m(zeta);
882  psi_h = sfuns.calc_psi_h(zeta);
883  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0_arr(i,j,k)) - psi_m);
884  ++iter;
885  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
886  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
887 
888  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
889  (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h);
890  olen_arr(i,j,k) = Olen;
891  }
892 
893 private:
896  const amrex::Real tol = 1.0e-5;
897  const amrex::Real WSMIN = 0.1; // minimum wind speed
898 };
899 
900 
901 /**
902  * Surface temperature with charnock roughness
903  */
905 {
906  surface_temp_charnock (amrex::Real zref,
907  amrex::Real flux,
908  amrex::Real cnk_a,
909  bool cnk_visc)
910  {
911  mdata.zref = zref;
912  mdata.surf_temp_flux = flux;
913  mdata.Cnk_a = cnk_a;
914  mdata.visc = cnk_visc;
915  }
916 
917  AMREX_GPU_DEVICE
918  AMREX_FORCE_INLINE
919  void
920  iterate_flux (const int& i,
921  const int& j,
922  const int& k,
923  const int& max_iters,
924  const amrex::Array4<amrex::Real>& z0_arr,
925  const amrex::Array4<const amrex::Real>& umm_arr,
926  const amrex::Array4<const amrex::Real>& tm_arr,
927  const amrex::Array4<const amrex::Real>& tvm_arr,
928  const amrex::Array4<const amrex::Real>& qvm_arr,
929  const amrex::Array4<amrex::Real>& u_star_arr,
930  const amrex::Array4<amrex::Real>& w_star_arr,
931  const amrex::Array4<amrex::Real>& t_star_arr,
932  const amrex::Array4<amrex::Real>& q_star_arr,
933  const amrex::Array4<amrex::Real>& t_surf_arr,
934  const amrex::Array4<amrex::Real>& olen_arr,
935  const amrex::Array4<amrex::Real>& pblh_arr,
936  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
937  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
938  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
939  {
940  int iter = 0;
941  amrex::Real ustar = 0.0;
942  amrex::Real wstar = 0.0;
943  amrex::Real z0 = 0.0;
944  amrex::Real tflux = 0.0;
945  amrex::Real zeta = 0.0;
946  amrex::Real psi_m = 0.0;
947  amrex::Real psi_h = 0.0;
948  amrex::Real Olen = 0.0;
949  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
950  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
951  do {
952  ustar = u_star_arr(i,j,k);
953  if (mdata.Cnk_a > 0) {
954  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
955  if (mdata.visc) {
956  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, 0.05);
957  }
958  } else {
959  z0 = COARE3_roughness(mdata.zref, umm, ustar);
960  }
961  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
962  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
963  tflux *= (1 + 0.61*qvm_arr(i,j,k));
964  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
965  if (w_star_arr) {
966  // update w* and Umagmean
967  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
968  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
969  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
970  umm = std::max(umm, WSMIN);
971  }
972  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
973  zeta = mdata.zref / Olen;
974  psi_m = sfuns.calc_psi_m(zeta);
975  psi_h = sfuns.calc_psi_h(zeta);
976  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
977  ++iter;
978  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
979  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
980 
981  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
982  (std::log(mdata.zref / z0) - psi_h);
983  olen_arr(i,j,k) = Olen;
984  z0_arr(i,j,k) = z0;
985  }
986 
987 private:
990  const amrex::Real tol = 1.0e-5;
991  const amrex::Real WSMIN = 0.1; // minimum wind speed
992 };
993 
994 
995 /**
996  * Surface temperature with modified charnock roughness
997  */
999 {
1000  surface_temp_mod_charnock (amrex::Real zref,
1001  amrex::Real flux,
1002  amrex::Real depth)
1003  {
1004  mdata.zref = zref;
1005  mdata.surf_temp_flux = flux;
1006  mdata.Cnk_d = depth;
1007  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
1008  }
1009 
1010  AMREX_GPU_DEVICE
1011  AMREX_FORCE_INLINE
1012  void
1013  iterate_flux (const int& i,
1014  const int& j,
1015  const int& k,
1016  const int& max_iters,
1017  const amrex::Array4<amrex::Real>& z0_arr,
1018  const amrex::Array4<const amrex::Real>& umm_arr,
1019  const amrex::Array4<const amrex::Real>& tm_arr,
1020  const amrex::Array4<const amrex::Real>& tvm_arr,
1021  const amrex::Array4<const amrex::Real>& qvm_arr,
1022  const amrex::Array4<amrex::Real>& u_star_arr,
1023  const amrex::Array4<amrex::Real>& w_star_arr,
1024  const amrex::Array4<amrex::Real>& t_star_arr,
1025  const amrex::Array4<amrex::Real>& q_star_arr,
1026  const amrex::Array4<amrex::Real>& t_surf_arr,
1027  const amrex::Array4<amrex::Real>& olen_arr,
1028  const amrex::Array4<amrex::Real>& pblh_arr,
1029  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1030  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1031  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1032  {
1033  int iter = 0;
1034  amrex::Real ustar = 0.0;
1035  amrex::Real wstar = 0.0;
1036  amrex::Real z0 = 0.0;
1037  amrex::Real tflux = 0.0;
1038  amrex::Real zeta = 0.0;
1039  amrex::Real psi_m = 0.0;
1040  amrex::Real psi_h = 0.0;
1041  amrex::Real Olen = 0.0;
1042  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1043  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
1044  do {
1045  ustar = u_star_arr(i,j,k);
1046  z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) );
1047  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1048  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
1049  tflux *= (1 + 0.61*qvm_arr(i,j,k));
1050  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
1051  if (w_star_arr) {
1052  // update w* and Umagmean
1053  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1054  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1055  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1056  umm = std::max(umm, WSMIN);
1057  }
1058  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1059  zeta = mdata.zref / Olen;
1060  psi_m = sfuns.calc_psi_m(zeta);
1061  psi_h = sfuns.calc_psi_h(zeta);
1062  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
1063  ++iter;
1064  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1065  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
1066 
1067  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1068  (std::log(mdata.zref / z0) - psi_h);
1069  olen_arr(i,j,k) = Olen;
1070  z0_arr(i,j,k) = z0;
1071  }
1072 
1073 private:
1076  const amrex::Real tol = 1.0e-5;
1077  const amrex::Real WSMIN = 0.1; // minimum wind speed
1078 };
1079 
1080 
1081 /**
1082  * Surface temperature with donelan roughness
1083  */
1085 {
1086  surface_temp_donelan (amrex::Real zref,
1087  amrex::Real flux)
1088  {
1089  mdata.zref = zref;
1090  mdata.surf_temp_flux = flux;
1091  }
1092 
1093  AMREX_GPU_DEVICE
1094  AMREX_FORCE_INLINE
1095  void
1096  iterate_flux (const int& i,
1097  const int& j,
1098  const int& k,
1099  const int& max_iters,
1100  const amrex::Array4<amrex::Real>& z0_arr,
1101  const amrex::Array4<const amrex::Real>& umm_arr,
1102  const amrex::Array4<const amrex::Real>& tm_arr,
1103  const amrex::Array4<const amrex::Real>& tvm_arr,
1104  const amrex::Array4<const amrex::Real>& qvm_arr,
1105  const amrex::Array4<amrex::Real>& u_star_arr,
1106  const amrex::Array4<amrex::Real>& w_star_arr,
1107  const amrex::Array4<amrex::Real>& t_star_arr,
1108  const amrex::Array4<amrex::Real>& q_star_arr,
1109  const amrex::Array4<amrex::Real>& t_surf_arr,
1110  const amrex::Array4<amrex::Real>& olen_arr,
1111  const amrex::Array4<amrex::Real>& pblh_arr,
1112  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1113  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1114  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1115  {
1116  int iter = 0;
1117  amrex::Real ustar = 0.0;
1118  amrex::Real wstar = 0.0;
1119  amrex::Real z0 = 0.0;
1120  amrex::Real tflux = 0.0;
1121  amrex::Real zeta = 0.0;
1122  amrex::Real psi_m = 0.0;
1123  amrex::Real psi_h = 0.0;
1124  amrex::Real Olen = 0.0;
1125  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1126  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
1127  do {
1128  ustar = u_star_arr(i,j,k);
1129  z0 = Donelan_roughness(ustar);
1130  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1131  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
1132  tflux *= (1 + 0.61*qvm_arr(i,j,k));
1133  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
1134  if (w_star_arr) {
1135  // update w* and Umagmean
1136  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1137  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1138  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1139  umm = std::max(umm, WSMIN);
1140  }
1141  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1142  zeta = mdata.zref / Olen;
1143  psi_m = sfuns.calc_psi_m(zeta);
1144  psi_h = sfuns.calc_psi_h(zeta);
1145  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
1146  ++iter;
1147  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1148  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
1149 
1150  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1151  (std::log(mdata.zref / z0) - psi_h);
1152  olen_arr(i,j,k) = Olen;
1153  z0_arr(i,j,k) = z0;
1154  }
1155 
1156 private:
1159  const amrex::Real tol = 1.0e-5;
1160  const amrex::Real WSMIN = 0.1; // minimum wind speed
1161 };
1162 
1163 
1164 /**
1165  * Surface temperature with wave-coupled roughness
1166  */
1168 {
1169  surface_temp_wave_coupled (amrex::Real zref,
1170  amrex::Real flux)
1171  {
1172  mdata.zref = zref;
1173  mdata.surf_temp_flux = flux;
1174  }
1175 
1176  AMREX_GPU_DEVICE
1177  AMREX_FORCE_INLINE
1178  void
1179  iterate_flux (const int& i,
1180  const int& j,
1181  const int& k,
1182  const int& max_iters,
1183  const amrex::Array4<amrex::Real>& z0_arr,
1184  const amrex::Array4<const amrex::Real>& umm_arr,
1185  const amrex::Array4<const amrex::Real>& tm_arr,
1186  const amrex::Array4<const amrex::Real>& tvm_arr,
1187  const amrex::Array4<const amrex::Real>& qvm_arr,
1188  const amrex::Array4<amrex::Real>& u_star_arr,
1189  const amrex::Array4<amrex::Real>& w_star_arr,
1190  const amrex::Array4<amrex::Real>& t_star_arr,
1191  const amrex::Array4<amrex::Real>& q_star_arr,
1192  const amrex::Array4<amrex::Real>& t_surf_arr,
1193  const amrex::Array4<amrex::Real>& olen_arr,
1194  const amrex::Array4<amrex::Real>& pblh_arr,
1195  const amrex::Array4<amrex::Real>& Hwave_arr,
1196  const amrex::Array4<amrex::Real>& Lwave_arr,
1197  const amrex::Array4<amrex::Real>& eta_arr) const
1198  {
1199  int iter = 0;
1200  amrex::Real ustar = 0.0;
1201  amrex::Real wstar = 0.0;
1202  amrex::Real z0 = 0.0;
1203  amrex::Real tflux = 0.0;
1204  amrex::Real zeta = 0.0;
1205  amrex::Real psi_m = 0.0;
1206  amrex::Real psi_h = 0.0;
1207  amrex::Real Olen = 0.0;
1208  int ie, je;
1209  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1210  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1211  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1212  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1213  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1214  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(mdata.zref / z0_arr(i,j,k));
1215  do {
1216  ustar = u_star_arr(i,j,k);
1217  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 )
1218  + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
1219  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1220  (std::log(mdata.zref / z0) - psi_h); // <w'T'>
1221  tflux *= (1 + 0.61*qvm_arr(i,j,k));
1222  tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= <w'Tv'>
1223  if (w_star_arr) {
1224  // update w* and Umagmean
1225  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1226  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1227  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1228  umm = std::max(umm, WSMIN);
1229  }
1230  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1231  zeta = mdata.zref / Olen;
1232  psi_m = sfuns.calc_psi_m(zeta);
1233  psi_h = sfuns.calc_psi_h(zeta);
1234  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(mdata.zref / z0) - psi_m);
1235  ++iter;
1236  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1237  AMREX_ASSERT_WITH_MESSAGE(iter < max_iters, "Maximum number of MOST iterations reached.");
1238 
1239  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1240  (std::log(mdata.zref / z0) - psi_h);
1241  olen_arr(i,j,k) = Olen;
1242  z0_arr(i,j,k) = z0;
1243  }
1244 
1245 private:
1248  const amrex::Real tol = 1.0e-5;
1249  const amrex::Real eps = 1e-15;
1250  const amrex::Real z0_eps = 1.0e-6;
1251  const amrex::Real z0_max = 0.1;
1252  const amrex::Real WSMIN = 0.1; // minimum wind speed
1253 };
1254 
1255 
1256 /**
1257  * Moeng flux formulation
1258  */
1260 {
1261  moeng_flux (int l_zlo)
1262  : zlo(l_zlo) {}
1263 
1264 
1265  AMREX_GPU_DEVICE
1266  AMREX_FORCE_INLINE
1267  amrex::Real
1268  compute_q_flux (const int& i,
1269  const int& j,
1270  const int& k,
1271  const int& n,
1272  const int& icomp,
1273  const amrex::Real& dz,
1274  const amrex::Real& dz1,
1275  const bool& exp_most,
1276  const amrex::Array4<const amrex::Real>& eta_arr,
1277  const amrex::Array4<const amrex::Real>& cons_arr,
1278  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1279  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1280  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1281  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
1282  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1283  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
1284  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1285  const amrex::Array4<amrex::Real>& dest_arr) const
1286  {
1287  amrex::Real rho, eta;
1288 
1289  int ic, jc;
1290  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1291  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1292  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1293  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1294 
1295  rho = cons_arr(ic,jc,zlo,Rho_comp);
1296 
1297  // TODO: Integrate MOST with moisture and MOENG FLUX type
1298  amrex::Real deltaz = dz * (zlo - k);
1299 
1300  // NOTE: this is rho*<q'w'> = -K dqdz
1301  amrex::Real moflux = 0.0;
1302 
1303  if (exp_most) {
1304  // surface gradient equal to gradient at first zface
1305  amrex::Real rqvgrad = (cons_arr(ic,jc,zlo+1,RhoQ1_comp) - cons_arr(ic,jc,zlo ,RhoQ1_comp)) / (0.5*(dz+dz1));
1306  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
1307  } else {
1308  int ie, je;
1309  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1310  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1311  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1312  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1313  eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
1314  eta = amrex::max(eta,eta_eps);
1315  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1316  }
1317 
1318  return moflux;
1319  }
1320 
1321  AMREX_GPU_DEVICE
1322  AMREX_FORCE_INLINE
1323  amrex::Real
1324  compute_t_flux (const int& i,
1325  const int& j,
1326  const int& k,
1327  const int& n,
1328  const int& icomp,
1329  const amrex::Real& dz,
1330  const amrex::Real& dz1,
1331  const bool& exp_most,
1332  const amrex::Array4<const amrex::Real>& eta_arr,
1333  const amrex::Array4<const amrex::Real>& cons_arr,
1334  const amrex::Array4<const amrex::Real>& velx_arr,
1335  const amrex::Array4<const amrex::Real>& vely_arr,
1336  const amrex::Array4<const amrex::Real>& umm_arr,
1337  const amrex::Array4<const amrex::Real>& tm_arr,
1338  const amrex::Array4<const amrex::Real>& u_star_arr,
1339  const amrex::Array4<const amrex::Real>& t_star_arr,
1340  const amrex::Array4<const amrex::Real>& t_surf_arr,
1341  const amrex::Array4<amrex::Real>& dest_arr) const
1342  {
1343  amrex::Real velx, vely, rho, theta, eta;
1344  int ix, jx, iy, jy, ic, jc;
1345 
1346  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
1347  jx = j < lbound(velx_arr).y ? lbound(velx_arr).y : j;
1348  ix = ix > ubound(velx_arr).x-1 ? ubound(velx_arr).x-1 : ix;
1349  jx = jx > ubound(velx_arr).y ? ubound(velx_arr).y : jx;
1350 
1351  iy = i < lbound(vely_arr).x ? lbound(vely_arr).x : i;
1352  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1353  iy = iy > ubound(vely_arr).x ? ubound(vely_arr).x : iy;
1354  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1355 
1356  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1357  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1358  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1359  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1360 
1361  velx = 0.5 *( velx_arr(ix,jx,zlo) + velx_arr(ix+1,jx ,zlo) );
1362  vely = 0.5 *( vely_arr(iy,jy,zlo) + vely_arr(iy ,jy+1,zlo) );
1363  rho = cons_arr(ic,jc,zlo,Rho_comp);
1364  theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho;
1365 
1366  amrex::Real theta_mean = tm_arr(ic,jc,zlo);
1367  amrex::Real wsp_mean = umm_arr(ic,jc,zlo);
1368  amrex::Real ustar = u_star_arr(ic,jc,zlo);
1369  amrex::Real tstar = t_star_arr(ic,jc,zlo);
1370  amrex::Real theta_surf = t_surf_arr(ic,jc,zlo);
1371 
1372  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1373  amrex::Real num1 = wsp * (theta_mean-theta_surf);
1374  amrex::Real num2 = wsp_mean * (theta-theta_mean);
1375  amrex::Real deltaz = dz * (zlo - k);
1376 
1377  wsp_mean = std::max(wsp_mean, WSMIN);
1378 
1379  // NOTE: this is rho*<T'w'> = -K dTdz
1380  amrex::Real moflux = (std::abs(tstar) > eps) ?
1381  -rho*tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0;
1382 
1383  if (exp_most) {
1384  // surface gradient equal to gradient at first zface
1385  amrex::Real rthetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1));
1386  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
1387  } else {
1388  int ie, je;
1389  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1390  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1391  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1392  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1393  eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
1394  eta = amrex::max(eta,eta_eps);
1395  // Note: Kh = eta/rho
1396  // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface
1397  // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel
1398  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1399  }
1400 
1401  return moflux;
1402  }
1403 
1404  AMREX_GPU_DEVICE
1405  AMREX_FORCE_INLINE
1406  amrex::Real
1407  compute_u_flux (const int& i,
1408  const int& j,
1409  const int& k,
1410  const int& icomp,
1411  const amrex::Real& dz,
1412  const amrex::Real& dz1,
1413  const bool& exp_most,
1414  const amrex::Array4<const amrex::Real>& eta_arr,
1415  const amrex::Array4<const amrex::Real>& cons_arr,
1416  const amrex::Array4<const amrex::Real>& velx_arr,
1417  const amrex::Array4<const amrex::Real>& vely_arr,
1418  const amrex::Array4<const amrex::Real>& umm_arr,
1419  const amrex::Array4<const amrex::Real>& um_arr,
1420  const amrex::Array4<const amrex::Real>& u_star_arr,
1421  const amrex::Array4<amrex::Real>& dest_arr) const
1422  {
1423  amrex::Real velx, vely, rho, eta;
1424  int jy, ic, jc;
1425 
1426  int iylo = i <= lbound(vely_arr).x ? lbound(vely_arr).x : i-1;
1427  int iyhi = i > ubound(vely_arr).x ? ubound(vely_arr).x : i;
1428 
1429  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1430  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1431 
1432  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
1433  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1434  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1435  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1436 
1437  velx = velx_arr(i,j,zlo);
1438  vely = 0.25*( vely_arr(iyhi,jy,zlo)+vely_arr(iyhi,jy+1,zlo)
1439  + vely_arr(iylo,jy,zlo)+vely_arr(iylo,jy+1,zlo) );
1440  rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
1441  + cons_arr(ic ,jc,zlo,Rho_comp) );
1442 
1443  amrex::Real umean = um_arr(i,j,zlo);
1444  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) );
1445  amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
1446 
1447  // Note: The surface mean shear stress is decomposed into tau_xz by
1448  // multiplying the modeled shear stress (rho*ustar^2) with
1449  // a factor of umean/wsp_mean for directionality; this factor
1450  // modifies the denominator from what is in Moeng 1984.
1451  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1452  amrex::Real num1 = wsp * umean;
1453  amrex::Real num2 = wsp_mean * (velx-umean);
1454  amrex::Real deltaz = dz * (zlo - k);
1455 
1456  wsp_mean = std::max(wsp_mean, WSMIN);
1457 
1458  // NOTE: this is rho*<u'w'> = -K dudz
1459  amrex::Real stressx = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
1460 
1461  if (exp_most) {
1462  // surface gradient equal to gradient at first zface
1463  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
1464  dest_arr(i,j,k,icomp) = velx - ugrad * deltaz;
1465  } else {
1466  int ie, je;
1467  ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
1468  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1469  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1470  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1471  eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
1472  + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) );
1473  eta = amrex::max(eta,eta_eps);
1474  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressx/eta*deltaz;
1475  }
1476 
1477  return stressx;
1478  }
1479 
1480  AMREX_GPU_DEVICE
1481  AMREX_FORCE_INLINE
1482  amrex::Real
1483  compute_v_flux (const int& i,
1484  const int& j,
1485  const int& k,
1486  const int& icomp,
1487  const amrex::Real& dz,
1488  const amrex::Real& dz1,
1489  const bool& exp_most,
1490  const amrex::Array4<const amrex::Real>& eta_arr,
1491  const amrex::Array4<const amrex::Real>& cons_arr,
1492  const amrex::Array4<const amrex::Real>& velx_arr,
1493  const amrex::Array4<const amrex::Real>& vely_arr,
1494  const amrex::Array4<const amrex::Real>& umm_arr,
1495  const amrex::Array4<const amrex::Real>& vm_arr,
1496  const amrex::Array4<const amrex::Real>& u_star_arr,
1497  const amrex::Array4<amrex::Real>& dest_arr) const
1498  {
1499  amrex::Real velx, vely, rho, eta;
1500  int ix, ic, jc;
1501 
1502  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
1503  ix = ix > ubound(velx_arr).x ? ubound(velx_arr).x : ix;
1504 
1505  int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
1506  int jxhi = j > ubound(velx_arr).y ? ubound(velx_arr).y : j;
1507 
1508  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1509  jc = j < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
1510  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1511  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1512 
1513  velx = 0.25*( velx_arr(ix,jxhi,zlo)+velx_arr(ix+1,jxhi,zlo)
1514  + velx_arr(ix,jxlo,zlo)+velx_arr(ix+1,jxlo,zlo) );
1515  vely = vely_arr(i,j,zlo);
1516  rho = 0.5*( cons_arr(ic,jc-1,zlo,Rho_comp)
1517  + cons_arr(ic,jc ,zlo,Rho_comp) );
1518 
1519  amrex::Real vmean = vm_arr(i,j,zlo);
1520  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) );
1521  amrex::Real ustar = 0.5 * ( u_star_arr(ic,jc-1,zlo) + u_star_arr(ic,jc,zlo) );
1522 
1523  // Note: The surface mean shear stress is decomposed into tau_yz by
1524  // multiplying the modeled shear stress (rho*ustar^2) with
1525  // a factor of vmean/wsp_mean for directionality; this factor
1526  // modifies the denominator from what is in Moeng 1984.
1527  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1528  amrex::Real num1 = wsp * vmean;
1529  amrex::Real num2 = wsp_mean * (vely-vmean);
1530  amrex::Real deltaz = dz * (zlo - k);
1531 
1532  wsp_mean = std::max(wsp_mean, WSMIN);
1533 
1534  // NOTE: this is rho*<v'w'> = -K dvdz
1535  amrex::Real stressy = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
1536 
1537  if (exp_most) {
1538  // surface gradient equal to gradient at first zface
1539  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
1540  dest_arr(i,j,k,icomp) = vely - vgrad * deltaz;
1541  } else {
1542  int ie, je;
1543  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1544  je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
1545  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1546  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1547  eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
1548  + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) );
1549  eta = amrex::max(eta,eta_eps);
1550  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressy/eta*deltaz;
1551  }
1552 
1553  return stressy;
1554  }
1555 
1556 private:
1557  int zlo;
1558  const amrex::Real eps = 1e-15;
1559  const amrex::Real eta_eps = 1e-8;
1560  const amrex::Real WSMIN = 0.1; // minimum wind speed
1561 };
1562 
1563 
1564 /**
1565  * Donelan flux formulation
1566  */
1568 {
1569  donelan_flux (int l_zlo)
1570  : zlo(l_zlo) {}
1571 
1572 
1573  AMREX_GPU_DEVICE
1574  AMREX_FORCE_INLINE
1575  amrex::Real
1576  compute_q_flux (const int& i,
1577  const int& j,
1578  const int& k,
1579  const int& n,
1580  const int& icomp,
1581  const amrex::Real& dz,
1582  const amrex::Real& dz1,
1583  const bool& exp_most,
1584  const amrex::Array4<const amrex::Real>& eta_arr,
1585  const amrex::Array4<const amrex::Real>& cons_arr,
1586  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1587  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1588  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1589  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
1590  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1591  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
1592  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1593  const amrex::Array4<amrex::Real>& dest_arr) const
1594  {
1595  amrex::Real rho, eta;
1596 
1597  int ic, jc;
1598  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1599  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1600  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1601  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1602 
1603  rho = cons_arr(ic,jc,zlo,Rho_comp);
1604 
1605  // TODO: Integrate MOST with moisture and DONELAN FLUX type
1606  amrex::Real deltaz = dz * (zlo - k);
1607 
1608  // NOTE: this is rho*<q'w'> = -K dqdz
1609  amrex::Real moflux = 0.0;
1610 
1611  if (exp_most) {
1612  // surface gradient equal to gradient at first zface
1613  amrex::Real rqvgrad = (cons_arr(ic,jc,zlo+1,RhoQ1_comp) - cons_arr(ic,jc,zlo ,RhoQ1_comp)) / (0.5*(dz+dz1));
1614  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
1615  } else {
1616  int ie, je;
1617  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1618  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1619  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1620  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1621  eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
1622  eta = amrex::max(eta,eta_eps);
1623  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1624  }
1625 
1626  return moflux;
1627  }
1628 
1629  AMREX_GPU_DEVICE
1630  AMREX_FORCE_INLINE
1631  amrex::Real
1632  compute_t_flux (const int& i,
1633  const int& j,
1634  const int& k,
1635  const int& n,
1636  const int& icomp,
1637  const amrex::Real& dz,
1638  const amrex::Real& dz1,
1639  const bool& exp_most,
1640  const amrex::Array4<const amrex::Real>& eta_arr,
1641  const amrex::Array4<const amrex::Real>& cons_arr,
1642  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1643  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1644  const amrex::Array4<const amrex::Real>& umm_arr,
1645  const amrex::Array4<const amrex::Real>& tm_arr,
1646  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1647  const amrex::Array4<const amrex::Real>& /*t_star_arr*/,
1648  const amrex::Array4<const amrex::Real>& t_surf_arr,
1649  const amrex::Array4<amrex::Real>& dest_arr) const
1650  {
1651  amrex::Real rho, eta;
1652 
1653  int ic, jc;
1654  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1655  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1656  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1657  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1658 
1659  rho = cons_arr(ic,jc,zlo,Rho_comp);
1660 
1661  amrex::Real Ch = 0.0012;
1662  amrex::Real wsp_mean = umm_arr(ic,jc,zlo);
1663  amrex::Real theta_surf = t_surf_arr(ic,jc,zlo);
1664  amrex::Real theta_mean = tm_arr(ic,jc,zlo);
1665  amrex::Real deltaz = dz * (zlo - k);
1666 
1667  // NOTE: this is rho*<T'w'> = -K dTdz
1668  amrex::Real moflux = -rho * Ch * wsp_mean * (theta_mean - theta_surf);
1669 
1670  if (exp_most) {
1671  // surface gradient equal to gradient at first zface
1672  amrex::Real rthetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1));
1673  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
1674  } else {
1675  int ie, je;
1676  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1677  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1678  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1679  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1680  eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
1681  eta = amrex::max(eta,eta_eps);
1682  // Note: Kh = eta/rho
1683  // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface
1684  // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel
1685  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1686  }
1687 
1688  return moflux;
1689  }
1690 
1691  AMREX_GPU_DEVICE
1692  AMREX_FORCE_INLINE
1693  amrex::Real
1694  compute_u_flux (const int& i,
1695  const int& j,
1696  const int& k,
1697  const int& icomp,
1698  const amrex::Real& dz,
1699  const amrex::Real& dz1,
1700  const bool& exp_most,
1701  const amrex::Array4<const amrex::Real>& eta_arr,
1702  const amrex::Array4<const amrex::Real>& cons_arr,
1703  const amrex::Array4<const amrex::Real>& velx_arr,
1704  const amrex::Array4<const amrex::Real>& vely_arr,
1705  const amrex::Array4<const amrex::Real>& umm_arr,
1706  const amrex::Array4<const amrex::Real>& /*um_arr*/,
1707  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1708  const amrex::Array4<amrex::Real>& dest_arr) const
1709  {
1710  amrex::Real velx, vely, rho, eta;
1711  int jy, ic, jc;
1712 
1713  int iylo = i <= lbound(vely_arr).x ? lbound(vely_arr).x : i-1;
1714  int iyhi = i > ubound(vely_arr).x ? ubound(vely_arr).x : i;
1715 
1716  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1717  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1718 
1719  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
1720  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1721  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1722  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1723 
1724  velx = velx_arr(i,j,zlo);
1725  vely = 0.25*( vely_arr(iyhi,jy,zlo)+vely_arr(iyhi,jy+1,zlo)
1726  + vely_arr(iylo,jy,zlo)+vely_arr(iylo,jy+1,zlo) );
1727  rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
1728  + cons_arr(ic ,jc,zlo,Rho_comp) );
1729 
1730  amrex::Real Cd = 0.001;
1731  const amrex::Real c = 7e-5;
1732  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1733  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) );
1734  if (wsp_mean <= 5.0) {
1735  Cd = 0.001;
1736  } else if (wsp_mean < 25.0 && wsp_mean > 5.0) {
1737  Cd = 0.001 + c * (wsp_mean - 5.0);
1738  } else {
1739  Cd = 0.0024;
1740  }
1741  amrex::Real deltaz = dz * (zlo - k);
1742 
1743  // NOTE: this is rho*<u'w'> = -K dudz
1744  amrex::Real stressx = -rho * Cd * velx * wsp;
1745 
1746  if (exp_most) {
1747  // surface gradient equal to gradient at first zface
1748  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
1749  dest_arr(i,j,k,icomp) = velx - ugrad * deltaz;
1750  } else {
1751  int ie, je;
1752  ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
1753  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1754  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1755  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1756  eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
1757  + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) );
1758  eta = amrex::max(eta,eta_eps);
1759  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressx/eta*deltaz;
1760  }
1761 
1762  return stressx;
1763  }
1764 
1765  AMREX_GPU_DEVICE
1766  AMREX_FORCE_INLINE
1767  amrex::Real
1768  compute_v_flux (const int& i,
1769  const int& j,
1770  const int& k,
1771  const int& icomp,
1772  const amrex::Real& dz,
1773  const amrex::Real& dz1,
1774  const bool& exp_most,
1775  const amrex::Array4<const amrex::Real>& eta_arr,
1776  const amrex::Array4<const amrex::Real>& cons_arr,
1777  const amrex::Array4<const amrex::Real>& velx_arr,
1778  const amrex::Array4<const amrex::Real>& vely_arr,
1779  const amrex::Array4<const amrex::Real>& umm_arr,
1780  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
1781  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1782  const amrex::Array4<amrex::Real>& dest_arr) const
1783  {
1784  amrex::Real velx, vely, rho, eta;
1785  int ix, ic, jc;
1786 
1787  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
1788  ix = ix > ubound(velx_arr).x ? ubound(velx_arr).x : ix;
1789 
1790  int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
1791  int jxhi = j > ubound(velx_arr).y ? ubound(velx_arr).y : j;
1792 
1793  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1794  jc = j < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
1795  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1796  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1797 
1798  velx = 0.25*( velx_arr(ix,jxhi,zlo)+velx_arr(ix+1,jxhi,zlo)
1799  + velx_arr(ix,jxlo,zlo)+velx_arr(ix+1,jxlo,zlo) );
1800  vely = vely_arr(i,j,zlo);
1801  rho = 0.5*( cons_arr(ic,jc-1,zlo,Rho_comp)
1802  + cons_arr(ic,jc ,zlo,Rho_comp) );
1803 
1804  amrex::Real Cd = 0.001;
1805  const amrex::Real c = 7e-5;
1806  amrex::Real wsp = sqrt(velx*velx+vely*vely);
1807  amrex::Real wsp_mean = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) );
1808  if (wsp_mean <= 5.0) {
1809  Cd = 0.001;
1810  } else if (wsp_mean < 25.0 && wsp_mean > 5.0) {
1811  Cd = 0.001 + c * (wsp_mean - 5.0);
1812  } else {
1813  Cd = 0.0024;
1814  }
1815  amrex::Real deltaz = dz * (zlo - k);
1816 
1817  // NOTE: this is rho*<v'w'> = -K dvdz
1818  amrex::Real stressy = -rho * Cd * vely * wsp;
1819 
1820  if (exp_most) {
1821  // surface gradient equal to gradient at first zface
1822  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
1823  dest_arr(i,j,k,icomp) = vely - vgrad * deltaz;
1824  } else {
1825  int ie, je;
1826  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1827  je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
1828  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1829  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1830  eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
1831  + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) );
1832  eta = amrex::max(eta,eta_eps);
1833  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressy/eta*deltaz;
1834  }
1835 
1836  return stressy;
1837  }
1838 
1839 private:
1840  int zlo;
1841  const amrex::Real eta_eps = 1e-8;
1842 };
1843 
1844 
1845 /**
1846  * Custom flux formulation
1847  */
1849 {
1850  custom_flux (int l_zlo)
1851  : zlo(l_zlo) {}
1852 
1853 
1854  AMREX_GPU_DEVICE
1855  AMREX_FORCE_INLINE
1856  amrex::Real
1857  compute_q_flux (const int& i,
1858  const int& j,
1859  const int& k,
1860  const int& n,
1861  const int& icomp,
1862  const amrex::Real& dz,
1863  const amrex::Real& dz1,
1864  const bool& exp_most,
1865  const amrex::Array4<const amrex::Real>& eta_arr,
1866  const amrex::Array4<const amrex::Real>& cons_arr,
1867  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1868  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1869  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1870  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
1871  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1872  const amrex::Array4<const amrex::Real>& q_star_arr,
1873  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1874  const amrex::Array4<amrex::Real>& dest_arr) const
1875  {
1876  amrex::Real rho, eta;
1877 
1878  int ic, jc;
1879  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1880  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1881  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1882  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1883 
1884  rho = cons_arr(ic,jc,zlo,Rho_comp);
1885 
1886  amrex::Real qstar = q_star_arr(ic,jc,zlo);
1887  amrex::Real deltaz = dz * (zlo - k);
1888 
1889  // NOTE: this is rho*<q'w'> = -K dqdz
1890  amrex::Real moflux = (std::abs(qstar) > eps) ? rho * qstar : 0.0;
1891 
1892  if (exp_most) {
1893  // surface gradient equal to gradient at first zface
1894  amrex::Real rqvgrad = ( cons_arr(ic,jc,zlo+1,RhoQ1_comp) - cons_arr(ic,jc,zlo ,RhoQ1_comp) ) / (0.5*(dz+dz1));
1895  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
1896  } else {
1897  int ie, je;
1898  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1899  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1900  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1901  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1902  eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
1903  eta = amrex::max(eta,eta_eps);
1904  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1905  }
1906 
1907  return moflux;
1908  }
1909 
1910  AMREX_GPU_DEVICE
1911  AMREX_FORCE_INLINE
1912  amrex::Real
1913  compute_t_flux (const int& i,
1914  const int& j,
1915  const int& k,
1916  const int& n,
1917  const int& icomp,
1918  const amrex::Real& dz,
1919  const amrex::Real& dz1,
1920  const bool& exp_most,
1921  const amrex::Array4<const amrex::Real>& eta_arr,
1922  const amrex::Array4<const amrex::Real>& cons_arr,
1923  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
1924  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
1925  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1926  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
1927  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
1928  const amrex::Array4<const amrex::Real>& t_star_arr,
1929  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
1930  const amrex::Array4<amrex::Real>& dest_arr) const
1931  {
1932  amrex::Real rho, eta;
1933 
1934  int ic, jc;
1935  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
1936  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1937  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1938  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1939 
1940  rho = cons_arr(ic,jc,zlo,Rho_comp);
1941 
1942  amrex::Real tstar = t_star_arr(ic,jc,zlo);
1943  amrex::Real deltaz = dz * (zlo - k);
1944 
1945  // NOTE: this is rho*<T'w'> = -K dTdz
1946  amrex::Real moflux = (std::abs(tstar) > eps) ? rho * tstar : 0.0;
1947 
1948  if (exp_most) {
1949  // surface gradient equal to gradient at first zface
1950  amrex::Real rthetagrad = ( cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo ,RhoTheta_comp) ) / (0.5*(dz+dz1));
1951  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
1952  } else {
1953  int ie, je;
1954  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1955  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1956  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1957  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1958  eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
1959  eta = amrex::max(eta,eta_eps);
1960  dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) + moflux*rho/eta*deltaz;
1961  }
1962 
1963  return moflux;
1964  }
1965 
1966  AMREX_GPU_DEVICE
1967  AMREX_FORCE_INLINE
1968  amrex::Real
1969  compute_u_flux (const int& i,
1970  const int& j,
1971  const int& k,
1972  const int& icomp,
1973  const amrex::Real& dz,
1974  const amrex::Real& dz1,
1975  const bool& exp_most,
1976  const amrex::Array4<const amrex::Real>& eta_arr,
1977  const amrex::Array4<const amrex::Real>& cons_arr,
1978  const amrex::Array4<const amrex::Real>& velx_arr,
1979  const amrex::Array4<const amrex::Real>& vely_arr,
1980  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
1981  const amrex::Array4<const amrex::Real>& /*um_arr*/,
1982  const amrex::Array4<const amrex::Real>& u_star_arr,
1983  const amrex::Array4<amrex::Real>& dest_arr) const
1984  {
1985  amrex::Real velx, vely, rho, eta;
1986  int jy, ic, jc;
1987 
1988  int iylo = i <= lbound(vely_arr).x ? lbound(vely_arr).x : i-1;
1989  int iyhi = i > ubound(vely_arr).x ? ubound(vely_arr).x : i;
1990 
1991  jy = j < lbound(vely_arr).y ? lbound(vely_arr).y : j;
1992  jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
1993 
1994  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
1995  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
1996  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
1997  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
1998 
1999  velx = velx_arr(i,j,zlo);
2000  vely = 0.25*( vely_arr(iyhi,jy,zlo)+vely_arr(iyhi,jy+1,zlo)
2001  + vely_arr(iylo,jy,zlo)+vely_arr(iylo,jy+1,zlo) );
2002  rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
2003  + cons_arr(ic ,jc,zlo,Rho_comp) );
2004 
2005  amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
2006  amrex::Real wsp = sqrt(velx*velx+vely*vely);
2007  amrex::Real deltaz = dz * (zlo - k);
2008 
2009  // NOTE: this is rho*<u'w'> = -K dudz
2010  amrex::Real stressx = -rho * ustar * ustar * velx / wsp;
2011 
2012  if (exp_most) {
2013  // surface gradient equal to gradient at first zface
2014  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
2015  dest_arr(i,j,k,icomp) = velx - ugrad * deltaz;
2016  } else {
2017  int ie, je;
2018  ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
2019  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
2020  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
2021  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
2022  eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
2023  + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) );
2024  eta = amrex::max(eta,eta_eps);
2025  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressx/eta*deltaz;
2026  }
2027 
2028  return stressx;
2029  }
2030 
2031  AMREX_GPU_DEVICE
2032  AMREX_FORCE_INLINE
2033  amrex::Real
2034  compute_v_flux (const int& i,
2035  const int& j,
2036  const int& k,
2037  const int& icomp,
2038  const amrex::Real& dz,
2039  const amrex::Real& dz1,
2040  const bool& exp_most,
2041  const amrex::Array4<const amrex::Real>& eta_arr,
2042  const amrex::Array4<const amrex::Real>& cons_arr,
2043  const amrex::Array4<const amrex::Real>& velx_arr,
2044  const amrex::Array4<const amrex::Real>& vely_arr,
2045  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2046  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2047  const amrex::Array4<const amrex::Real>& u_star_arr,
2048  const amrex::Array4<amrex::Real>& dest_arr) const
2049  {
2050  amrex::Real velx, vely, rho, eta;
2051  int ix, ic, jc;
2052 
2053  ix = i < lbound(velx_arr).x ? lbound(velx_arr).x : i;
2054  ix = ix > ubound(velx_arr).x ? ubound(velx_arr).x : ix;
2055 
2056  int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
2057  int jxhi = j > ubound(velx_arr).y ? ubound(velx_arr).y : j;
2058 
2059  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
2060  jc = j < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
2061  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2062  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2063 
2064  velx = 0.25*( velx_arr(ix,jxhi,zlo)+velx_arr(ix+1,jxhi,zlo)
2065  + velx_arr(ix,jxlo,zlo)+velx_arr(ix+1,jxlo,zlo) );
2066  vely = vely_arr(i,j,zlo);
2067  rho = 0.5*( cons_arr(ic,jc-1,zlo,Rho_comp)
2068  + cons_arr(ic,jc ,zlo,Rho_comp) );
2069 
2070  amrex::Real ustar = 0.5 * ( u_star_arr(ic,jc-1,zlo) + u_star_arr(ic,jc,zlo) );
2071  amrex::Real wsp = sqrt(velx*velx+vely*vely);
2072  amrex::Real deltaz = dz * (zlo - k);
2073 
2074  // NOTE: this is rho*<v'w'> = -K dvdz
2075  amrex::Real stressy = -rho * ustar * ustar * vely / wsp;
2076 
2077  if (exp_most) {
2078  // surface gradient equal to gradient at first zface
2079  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
2080  dest_arr(i,j,k,icomp) = vely - vgrad * deltaz;
2081  } else {
2082  int ie, je;
2083  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
2084  je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
2085  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
2086  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
2087  eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
2088  + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) );
2089  eta = amrex::max(eta,eta_eps);
2090  dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) + stressy/eta*deltaz;
2091  }
2092 
2093  return stressy;
2094  }
2095 
2096 private:
2097  int zlo;
2098  const amrex::Real eps = 1e-15;
2099  const amrex::Real eta_eps = 1e-8;
2100 };
2101 
2102 /**
2103  * Rotate flux formulation
2104  */
2106 {
2107  rotate_flux (int l_zlo)
2108  : zlo(l_zlo) {}
2109 
2110 
2111  AMREX_GPU_DEVICE
2112  AMREX_FORCE_INLINE
2113  amrex::Real
2114  compute_q_flux (const int& i,
2115  const int& j,
2116  const int& k,
2117  const int& n,
2118  const int& icomp,
2119  const amrex::Real& dz,
2120  const amrex::Real& dz1,
2121  const bool& /*exp_most*/,
2122  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2123  const amrex::Array4<const amrex::Real>& cons_arr,
2124  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2125  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2126  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2127  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
2128  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2129  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
2130  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
2131  const amrex::Array4<amrex::Real>& dest_arr) const
2132  {
2133  int ic, jc;
2134  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
2135  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
2136  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2137  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2138 
2139  // surface gradient equal to gradient at first zface
2140  amrex::Real deltaz = dz * (zlo - k);
2141  amrex::Real rqvgrad = ( cons_arr(ic,jc,zlo+1,RhoQ1_comp)
2142  - cons_arr(ic,jc,zlo ,RhoQ1_comp) ) / (0.5*(dz+dz1));
2143  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoQ1_comp) - rqvgrad * deltaz;
2144 
2145  // NOTE: this is the total stress
2146  amrex::Real rho = cons_arr(ic,jc,zlo,Rho_comp);
2147  amrex::Real qstar = 0.0; //q_star_arr(ic,jc,zlo);
2148  amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0;
2149 
2150  return moflux;
2151  }
2152 
2153  AMREX_GPU_DEVICE
2154  AMREX_FORCE_INLINE
2155  amrex::Real
2156  compute_t_flux (const int& i,
2157  const int& j,
2158  const int& k,
2159  const int& n,
2160  const int& icomp,
2161  const amrex::Real& dz,
2162  const amrex::Real& dz1,
2163  const bool& /*exp_most*/,
2164  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2165  const amrex::Array4<const amrex::Real>& cons_arr,
2166  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2167  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2168  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2169  const amrex::Array4<const amrex::Real>& tm_arr,
2170  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2171  const amrex::Array4<const amrex::Real>& t_star_arr,
2172  const amrex::Array4<const amrex::Real>& t_surf_arr,
2173  const amrex::Array4<amrex::Real>& dest_arr) const
2174  {
2175  int ic, jc;
2176  ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
2177  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
2178  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2179  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2180 
2181  // surface gradient equal to gradient at first zface
2182  amrex::Real deltaz = dz * (zlo - k);
2183  amrex::Real rthetagrad = ( cons_arr(ic,jc,zlo+1,RhoTheta_comp)
2184  - cons_arr(ic,jc,zlo ,RhoTheta_comp) ) / (0.5*(dz+dz1));
2185  dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rthetagrad * deltaz;
2186 
2187  // NOTE: this is the total stress
2188  amrex::Real theta_mean = tm_arr(ic,jc,zlo);
2189  amrex::Real theta_surf = t_surf_arr(ic,jc,zlo);
2190  amrex::Real rho = cons_arr(ic,jc,zlo,Rho_comp);
2191  amrex::Real tstar = t_star_arr(ic,jc,zlo);
2192  amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar * (theta_mean-theta_surf) : 0.0;
2193 
2194  return moflux;
2195  }
2196 
2197  AMREX_GPU_DEVICE
2198  AMREX_FORCE_INLINE
2199  amrex::Real
2200  compute_u_flux (const int& i,
2201  const int& j,
2202  const int& k,
2203  const int& icomp,
2204  const amrex::Real& dz,
2205  const amrex::Real& dz1,
2206  const bool& /*exp_most*/,
2207  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2208  const amrex::Array4<const amrex::Real>& cons_arr,
2209  const amrex::Array4<const amrex::Real>& velx_arr,
2210  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2211  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2212  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2213  const amrex::Array4<const amrex::Real>& u_star_arr,
2214  const amrex::Array4<amrex::Real>& dest_arr) const
2215  {
2216  int ic, jc;
2217  ic = i < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
2218  jc = j < lbound(cons_arr).y ? lbound(cons_arr).y : j;
2219  ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
2220  jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;
2221 
2222  // surface gradient equal to gradient at first zface
2223  amrex::Real deltaz = dz * (zlo - k);
2224  amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx_arr(i,j,zlo)) / (0.5*(dz+dz1));
2225  dest_arr(i,j,k,icomp) = velx_arr(i,j,zlo) - ugrad * deltaz;
2226 
2227  // NOTE: this is the total stress
2228  amrex::Real rho = 0.5 *( cons_arr(ic-1,jc,zlo,Rho_comp)
2229  + cons_arr(ic ,jc,zlo,Rho_comp) );
2230  amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
2231  amrex::Real stressx = -rho * ustar * ustar;
2232 
2233  return stressx;
2234  }
2235 
2236  AMREX_GPU_DEVICE
2237  AMREX_FORCE_INLINE
2238  amrex::Real
2239  compute_v_flux (const int& i,
2240  const int& j,
2241  const int& k,
2242  const int& icomp,
2243  const amrex::Real& dz,
2244  const amrex::Real& dz1,
2245  const bool& /*exp_most*/,
2246  const amrex::Array4<const amrex::Real>& /*eta_arr*/,
2247  const amrex::Array4<const amrex::Real>& /*cons_arr*/,
2248  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2249  const amrex::Array4<const amrex::Real>& vely_arr,
2250  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2251  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2252  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2253  const amrex::Array4<amrex::Real>& dest_arr) const
2254  {
2255  // surface gradient equal to gradient at first zface
2256  amrex::Real deltaz = dz * (zlo - k);
2257  amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely_arr(i,j,zlo)) / (0.5*(dz+dz1));
2258  dest_arr(i,j,k,icomp) = vely_arr(i,j,zlo) - vgrad * deltaz;
2259 
2260  // NOTE: this is the total stress
2261  amrex::Real stressy = 0.0;
2262 
2263  return stressy;
2264  }
2265 
2266 private:
2267  int zlo;
2268  const amrex::Real eps = 1e-15;
2269 };
2270 #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:168
@ Q_v
Definition: ERF_IndexDefines.H:171
@ Mom_v
Definition: ERF_IndexDefines.H:167
@ theta
Definition: ERF_MM5.H:20
@ rho
Definition: ERF_Kessler.H:22
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:1849
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:2034
int zlo
Definition: ERF_MOSTStress.H:2097
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:1969
const amrex::Real eta_eps
Definition: ERF_MOSTStress.H:2099
custom_flux(int l_zlo)
Definition: ERF_MOSTStress.H:1850
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:1913
const amrex::Real eps
Definition: ERF_MOSTStress.H:2098
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:1857
Definition: ERF_MOSTStress.H:1568
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:1694
donelan_flux(int l_zlo)
Definition: ERF_MOSTStress.H:1569
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:1632
const amrex::Real eta_eps
Definition: ERF_MOSTStress.H:1841
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:1768
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:1576
int zlo
Definition: ERF_MOSTStress.H:1840
Definition: ERF_MOSTStress.H:1260
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:1407
moeng_flux(int l_zlo)
Definition: ERF_MOSTStress.H:1261
int zlo
Definition: ERF_MOSTStress.H:1557
const amrex::Real eps
Definition: ERF_MOSTStress.H:1558
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:1268
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:1483
const amrex::Real eta_eps
Definition: ERF_MOSTStress.H:1559
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1560
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:1324
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:2106
int zlo
Definition: ERF_MOSTStress.H:2267
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:2114
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:2156
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:2200
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:2239
const amrex::Real eps
Definition: ERF_MOSTStress.H:2268
rotate_flux(int l_zlo)
Definition: ERF_MOSTStress.H:2107
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:478
const amrex::Real tol
Definition: ERF_MOSTStress.H:561
similarity_funs sfuns
Definition: ERF_MOSTStress.H:560
most_data mdata
Definition: ERF_MOSTStress.H:559
surface_flux_charnock(amrex::Real zref, amrex::Real flux, amrex::Real cnk_a, bool cnk_visc)
Definition: ERF_MOSTStress.H:479
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:493
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:562
Definition: ERF_MOSTStress.H:654
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:665
const amrex::Real tol
Definition: ERF_MOSTStress.H:726
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:727
most_data mdata
Definition: ERF_MOSTStress.H:724
similarity_funs sfuns
Definition: ERF_MOSTStress.H:725
surface_flux_donelan(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:655
Definition: ERF_MOSTStress.H:570
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:646
most_data mdata
Definition: ERF_MOSTStress.H:643
surface_flux_mod_charnock(amrex::Real zref, amrex::Real flux, amrex::Real depth)
Definition: ERF_MOSTStress.H:571
const amrex::Real tol
Definition: ERF_MOSTStress.H:645
similarity_funs sfuns
Definition: ERF_MOSTStress.H:644
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:584
Definition: ERF_MOSTStress.H:735
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:817
surface_flux_wave_coupled(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:736
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:746
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:816
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:815
most_data mdata
Definition: ERF_MOSTStress.H:811
similarity_funs sfuns
Definition: ERF_MOSTStress.H:812
const amrex::Real tol
Definition: ERF_MOSTStress.H:813
const amrex::Real eps
Definition: ERF_MOSTStress.H:814
Definition: ERF_MOSTStress.H:399
similarity_funs sfuns
Definition: ERF_MOSTStress.H:468
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:470
most_data mdata
Definition: ERF_MOSTStress.H:467
surface_flux(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:400
const amrex::Real tol
Definition: ERF_MOSTStress.H:469
Definition: ERF_MOSTStress.H:905
most_data mdata
Definition: ERF_MOSTStress.H:988
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:920
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:991
const amrex::Real tol
Definition: ERF_MOSTStress.H:990
surface_temp_charnock(amrex::Real zref, amrex::Real flux, amrex::Real cnk_a, bool cnk_visc)
Definition: ERF_MOSTStress.H:906
similarity_funs sfuns
Definition: ERF_MOSTStress.H:989
Definition: ERF_MOSTStress.H:1085
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1158
most_data mdata
Definition: ERF_MOSTStress.H:1157
surface_temp_donelan(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:1086
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:1096
const amrex::Real tol
Definition: ERF_MOSTStress.H:1159
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1160
Definition: ERF_MOSTStress.H:999
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1075
surface_temp_mod_charnock(amrex::Real zref, amrex::Real flux, amrex::Real depth)
Definition: ERF_MOSTStress.H:1000
const amrex::Real tol
Definition: ERF_MOSTStress.H:1076
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1077
most_data mdata
Definition: ERF_MOSTStress.H:1074
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:1013
Definition: ERF_MOSTStress.H:1168
const amrex::Real eps
Definition: ERF_MOSTStress.H:1249
const amrex::Real tol
Definition: ERF_MOSTStress.H:1248
most_data mdata
Definition: ERF_MOSTStress.H:1246
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:1179
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:1250
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1247
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:1251
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1252
surface_temp_wave_coupled(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:1169
Definition: ERF_MOSTStress.H:825
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:836
similarity_funs sfuns
Definition: ERF_MOSTStress.H:895
surface_temp(amrex::Real zref, amrex::Real flux)
Definition: ERF_MOSTStress.H:826
const amrex::Real tol
Definition: ERF_MOSTStress.H:896
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:897
most_data mdata
Definition: ERF_MOSTStress.H:894