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