ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MOSTStress.H
Go to the documentation of this file.
1 #ifndef ERF_MOSTStress_H
2 #define ERF_MOSTStress_H
3 
4 #include <ERF_Constants.H>
5 #include <ERF_IndexDefines.H>
6 #include <ERF_MOSTRoughness.H>
7 #include <ERF_Wstar.H>
8 
9 /**
10  * Structure of plain old data relevant to MOST BCs
11  */
12 struct most_data
13 {
14 public:
15  amrex::Real z0_const{amrex::Real(0.1)}; ///< Roughness height -- default constant value(m)
16  amrex::Real kappa{KAPPA}; ///< von Karman constant
17  amrex::Real gravity{CONST_GRAV}; ///< Acceleration due to gravity (m/s^2)
18  amrex::Real surf_temp_flux{zero}; ///< Heat flux TODO: decide whether this is <θ'w'> or <θv'w'> under moist conditions
19  amrex::Real surf_moist_flux{zero}; ///< Moisture flux
20 
21  amrex::Real Cnk_a{amrex::Real(0.0185)}; ///< Standard Charnock constant https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137.1)
22  amrex::Real Cnk_b1{one/amrex::Real(30.0)}; ///< Modified Charnock Eq (4) https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137.1)
23  amrex::Real Cnk_b2{amrex::Real(1260.0)}; ///< Modified Charnock Eq (4) https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137.1)
24  amrex::Real Cnk_d{amrex::Real(30.0)}; ///< Modified Charnock Eq (4) https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137.1)
26  bool visc{false}; ///< Use viscous Charnock formulation
27 
31 
32  const amrex::Real Bjr_beta = amrex::Real(1.2); // Empirical parameter from Beljaars 1995 QJRMS
33 };
34 
35 
36 /**
37  * Structure of similarity functions for Moeng formulation
38  */
40 {
41  //
42  // Similarity functions from Jimenez (2012)
43  //
44  /**
45  * Compute the Jimenez momentum stability correction.
46  *
47  * @param[in] zeta stability coordinate z/L
48  */
49  AMREX_GPU_HOST_DEVICE
50  AMREX_FORCE_INLINE
52  calc_psi_m2 (amrex::Real zeta) const
53  {
54  if (zeta > 0) {
55  amrex::Real x = std::pow(one + std::pow(zeta, amrex::Real(2.5)), one/amrex::Real(2.5));
56  return ( -amrex::Real(6.1)*std::log(zeta + x) );
57  } else {
58  amrex::Real x = std::pow(one - amrex::Real(16.0)*zeta, fourth);
59  amrex::Real psi_k_m = two * std::log(myhalf * (one + x)) + std::log(myhalf * (one + x * x)) -
60  two * std::atan(x) + PIoTwo;
61  amrex::Real y = std::pow(one - amrex::Real(10.0)*zeta, one/three);
62  amrex::Real psi_c_m = (three/two)*std::log((y*y + y + one)/three)
63  - std::sqrt(three)*std::atan((two*y + one)/std::sqrt(three))
64  + PI/std::sqrt(three);
65  return ( (psi_k_m + zeta*zeta*psi_c_m) / (one + zeta*zeta) );
66  }
67  }
68 
69  /**
70  * Compute the Jimenez heat stability correction.
71  *
72  * @param[in] zeta stability coordinate z/L
73  */
74  AMREX_GPU_HOST_DEVICE
75  AMREX_FORCE_INLINE
77  calc_psi_h2 (amrex::Real zeta) const
78  {
79  if (zeta > 0) {
80  amrex::Real x = std::pow(one + std::pow(zeta, amrex::Real(1.1)), one/amrex::Real(1.1));
81  return ( -amrex::Real(5.3)*std::log(zeta + x) );
82  } else {
83  amrex::Real x = std::sqrt(one - amrex::Real(16.0)*zeta);
84  amrex::Real psi_k_h = two * std::log(myhalf * (one + x));
85  amrex::Real y = std::pow(one - amrex::Real(34.0)*zeta, one/three);
86  amrex::Real psi_c_h = (three/two)*std::log((y*y + y + one)/three)
87  - std::sqrt(three)*std::atan((two*y + one)/std::sqrt(three))
88  + PI/std::sqrt(three);
89  return ( (psi_k_h + zeta*zeta*psi_c_h) / (one + zeta*zeta) );
90  }
91  }
92 
93 
94  //
95  // Similarity functions from Businger & Dyer (1966)
96  //
97  /**
98  * Compute the Businger-Dyer momentum stability correction.
99  *
100  * @param[in] zeta stability coordinate z/L
101  */
102  AMREX_GPU_HOST_DEVICE
103  AMREX_FORCE_INLINE
105  calc_psi_m (amrex::Real zeta) const
106  {
107  if (zeta > 0) {
108  return -beta_m * zeta;
109  } else {
110  amrex::Real x = std::sqrt(std::sqrt(one - gamma_m * zeta));
111  return two * std::log(myhalf * (one + x)) + std::log(myhalf * (one + x * x)) -
112  two * std::atan(x) + PIoTwo;
113  }
114  }
115 
116  /**
117  * Compute the Businger-Dyer heat stability correction.
118  *
119  * @param[in] zeta stability coordinate z/L
120  */
121  AMREX_GPU_HOST_DEVICE
122  AMREX_FORCE_INLINE
124  calc_psi_h (amrex::Real zeta) const
125  {
126  if (zeta > 0) {
127  return -beta_h * zeta;
128  } else {
129  amrex::Real x = std::sqrt(one - gamma_h * zeta);
130  return two * std::log(myhalf * (one + x));
131  }
132  }
133 
134 private:
135  amrex::Real beta_m{amrex::Real(5.0)}; ///< Constants from Dyer, BLM, 1974
136  amrex::Real beta_h{amrex::Real(5.0)}; ///< https://doi.org/amrex::Real(10.1007)/BF00240838
139 };
140 
141 
142 /**
143  * Empirical kinematic viscosity [m2/s] formula from Andreas (1989) CRREL Rep.
144  * 89-11, valid between -173 and 277 deg C.
145  *
146  * @param[in] T_degK air temperature in Kelvin
147  */
148 AMREX_GPU_DEVICE
149 AMREX_FORCE_INLINE
152 {
153  amrex::Real TC = T_degK - amrex::Real(273.15);
154  return amrex::Real(1.326e-5)*(one + amrex::Real(6.542e-3)*TC + amrex::Real(8.301e-6)*TC*TC - amrex::Real(4.84e-9)*TC*TC*TC);
155 }
156 
157 
158 /**
159  * Adiabatic with constant roughness
160  */
161 struct adiabatic
162 {
163  /**
164  * Construct adiabatic constant-roughness MOST data.
165  *
166  * @param[in] Tflux prescribed heat flux
167  * @param[in] Qvflux prescribed moisture flux
168  */
170  amrex::Real Qvflux)
171  {
172  mdata.surf_temp_flux = Tflux;
173  mdata.surf_moist_flux = Qvflux;
174  }
175 
176  /**
177  * Iterate MOST state for one surface point.
178  *
179  * @par Calling sequence
180  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
181  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
182  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
183  * eta_arr.
184  */
185  AMREX_GPU_DEVICE
186  AMREX_FORCE_INLINE
187  void
188  iterate_flux (const int& i,
189  const int& j,
190  const int& k,
191  const int& /*max_iters*/,
192  const amrex::Array4<const amrex::Real>& zref_arr,
193  const amrex::Array4<const amrex::Real>& z0_arr,
194  const amrex::Array4<const amrex::Real>& umm_arr,
195  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
196  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
197  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
198  const amrex::Array4<amrex::Real>& u_star_arr,
199  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
200  const amrex::Array4<amrex::Real>& t_star_arr,
201  const amrex::Array4<amrex::Real>& q_star_arr,
202  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
203  const amrex::Array4<amrex::Real>& /*q_surf_arr*/,
204  const amrex::Array4<amrex::Real>& olen_arr,
205  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
206  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
207  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
208  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
209  {
210  olen_arr(i,j,k) = bogus_large_value;
211  u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(zref_arr(i,j,k) / z0_arr(i,j,k));
212  t_star_arr(i,j,k) = zero;
213  q_star_arr(i,j,k) = zero;
214  }
215 
216 private:
219 };
220 
221 
222 /**
223  * Adiabatic with charnock roughness
224  */
226 {
227  /**
228  * Construct adiabatic Charnock-roughness MOST data.
229  *
230  * @param[in] Tflux prescribed heat flux
231  * @param[in] Qvflux prescribed moisture flux
232  * @param[in] cnk_a Charnock parameter
233  * @param[in] cnk_visc whether to include viscous roughness
234  */
236  amrex::Real Qvflux,
237  amrex::Real cnk_a,
238  bool cnk_visc)
239  {
240  mdata.surf_temp_flux = Tflux;
241  mdata.surf_moist_flux = Qvflux;
242  mdata.Cnk_a = cnk_a;
243  mdata.visc = cnk_visc;
244  }
245 
246  /**
247  * Iterate MOST state for one surface point.
248  *
249  * @par Calling sequence
250  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
251  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
252  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
253  * eta_arr.
254  */
255  AMREX_GPU_DEVICE
256  AMREX_FORCE_INLINE
257  void
258  iterate_flux (const int& i,
259  const int& j,
260  const int& k,
261  const int& max_iters,
262  const amrex::Array4<const amrex::Real>& zref_arr,
263  const amrex::Array4<amrex::Real>& z0_arr,
264  const amrex::Array4<const amrex::Real>& umm_arr,
265  const amrex::Array4<const amrex::Real>& tm_arr,
266  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
267  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
268  const amrex::Array4<amrex::Real>& u_star_arr,
269  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
270  const amrex::Array4<amrex::Real>& t_star_arr,
271  const amrex::Array4<amrex::Real>& q_star_arr,
272  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
273  const amrex::Array4<amrex::Real>& /*q_surf_arr*/,
274  const amrex::Array4<amrex::Real>& olen_arr,
275  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
276  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
277  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
278  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
279  {
280  amrex::Real ustar = zero;
281  amrex::Real zref = zref_arr(i,j,k);
282  amrex::Real z0 = z0_arr(i,j,k);
283  amrex::Real z0_old = z0;
284  amrex::Real psi_m = zero;
285  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
286  amrex::Real C = std::log(zref / z0);
287 
288  // Fixed point iteration on roughness
289  int iter_z = 0;
290  do {
291  // Transfer curr to old
292  z0_old = z0;
293 
294  // Update
295  C = std::log(zref / z0_old);
296  ustar = mdata.kappa * umm / (C - psi_m);
297  if (mdata.Cnk_a > 0) {
298  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
299  if (mdata.visc) {
300  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, amrex::Real(0.05));
301  }
302  } else {
303  z0 = COARE3_roughness(zref, umm, ustar);
304  }
305 
306  ++iter_z;
307  } while ( (std::abs(z0 - z0_old) > tol_z) && (iter_z <= max_iters) );
308  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter_z < max_iters,
309  "Maximum number of MOST roughness iterations reached.");
310  C = std::log(zref / z0);
311 
312  // Populate the stored MOST arrays
313  z0_arr(i,j,k) = z0;
314  olen_arr(i,j,k) = bogus_large_value;
315  u_star_arr(i,j,k) = mdata.kappa * umm / (C - psi_m);
316  t_star_arr(i,j,k) = zero;
317  q_star_arr(i,j,k) = zero;
318  }
319 
320 private:
323 #ifdef AMREX_USE_FLOAT
324  const amrex::Real tol_z = amrex::Real(1.0e-6);
325 #else
326  const amrex::Real tol_z = amrex::Real(1.0e-10);
327 #endif
328  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
329 };
330 
331 
332 /**
333  * Adiabatic with modified charnock roughness
334  */
336 {
337  /**
338  * Construct adiabatic modified-Charnock MOST data.
339  *
340  * @param[in] Tflux prescribed heat flux
341  * @param[in] Qvflux prescribed moisture flux
342  * @param[in] depth water depth for modified Charnock roughness
343  */
345  amrex::Real Qvflux,
346  amrex::Real depth)
347  {
348  mdata.surf_temp_flux = Tflux;
349  mdata.surf_moist_flux = Qvflux;
350  mdata.Cnk_d = depth;
351  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
352  }
353 
354  /**
355  * Iterate MOST state for one surface point.
356  *
357  * @par Calling sequence
358  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
359  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
360  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
361  * eta_arr.
362  */
363  AMREX_GPU_DEVICE
364  AMREX_FORCE_INLINE
365  void
366  iterate_flux (const int& i,
367  const int& j,
368  const int& k,
369  const int& max_iters,
370  const amrex::Array4<const amrex::Real>& zref_arr,
371  const amrex::Array4<amrex::Real>& z0_arr,
372  const amrex::Array4<const amrex::Real>& umm_arr,
373  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
374  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
375  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
376  const amrex::Array4<amrex::Real>& u_star_arr,
377  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
378  const amrex::Array4<amrex::Real>& t_star_arr,
379  const amrex::Array4<amrex::Real>& q_star_arr,
380  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
381  const amrex::Array4<amrex::Real>& /*q_surf_arr*/,
382  const amrex::Array4<amrex::Real>& olen_arr,
383  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
384  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
385  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
386  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
387  {
388  amrex::Real ustar = zero;
389  amrex::Real zref = zref_arr(i,j,k);
390  amrex::Real z0 = z0_arr(i,j,k);
391  amrex::Real z0_old = z0;
392  amrex::Real psi_m = zero;
393  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
394  amrex::Real C = std::log(zref / z0);
395 
396  // Fixed point iteration on roughness
397  int iter_z = 0;
398  do {
399  // Transfer curr to old
400  z0_old = z0;
401 
402  // Update
403  C = std::log(zref / z0_old);
404  ustar = mdata.kappa * umm / (C - psi_m);
405  z0 = std::exp( (amrex::Real(2.7)*ustar - amrex::Real(1.8)/mdata.Cnk_b) / (ustar + amrex::Real(0.17)/mdata.Cnk_b) );
406 
407  ++iter_z;
408  } while ( (std::abs(z0 - z0_old) > tol_z) && (iter_z <= max_iters) );
409  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter_z < max_iters,
410  "Maximum number of MOST roughness iterations reached.");
411  C = std::log(zref / z0);
412 
413  // Populate the stored MOST arrays
414  z0_arr(i,j,k) = z0;
415  olen_arr(i,j,k) = bogus_large_value;
416  u_star_arr(i,j,k) = mdata.kappa * umm / (C - psi_m);
417  t_star_arr(i,j,k) = zero;
418  q_star_arr(i,j,k) = zero;
419  }
420 
421 private:
424 #ifdef AMREX_USE_FLOAT
425  const amrex::Real tol_z = amrex::Real(1.0e-6);
426 #else
427  const amrex::Real tol_z = amrex::Real(1.0e-10);
428 #endif
429  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
430 };
431 
432 
433 /**
434  * Adiabatic with Donelan roughness
435  */
437 {
438  /**
439  * Construct adiabatic Donelan-roughness MOST data.
440  *
441  * @param[in] Tflux prescribed heat flux
442  * @param[in] Qvflux prescribed moisture flux
443  */
445  amrex::Real Qvflux)
446  {
447  mdata.surf_temp_flux = Tflux;
448  mdata.surf_moist_flux = Qvflux;
449  }
450 
451  /**
452  * Iterate MOST state for one surface point.
453  *
454  * @par Calling sequence
455  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
456  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
457  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
458  * eta_arr.
459  */
460  AMREX_GPU_DEVICE
461  AMREX_FORCE_INLINE
462  void
463  iterate_flux (const int& i,
464  const int& j,
465  const int& k,
466  const int& max_iters,
467  const amrex::Array4<const amrex::Real>& zref_arr,
468  const amrex::Array4<amrex::Real>& z0_arr,
469  const amrex::Array4<const amrex::Real>& umm_arr,
470  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
471  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
472  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
473  const amrex::Array4<amrex::Real>& u_star_arr,
474  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
475  const amrex::Array4<amrex::Real>& t_star_arr,
476  const amrex::Array4<amrex::Real>& q_star_arr,
477  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
478  const amrex::Array4<amrex::Real>& /*q_surf_arr*/,
479  const amrex::Array4<amrex::Real>& olen_arr,
480  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
481  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
482  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
483  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
484  {
485  int iter = 0;
486  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
487  amrex::Real ustar = zero;
488  amrex::Real zref = zref_arr(i,j,k);
489  amrex::Real z0 = zero;
490  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
491  if (u_star_arr(i,j,k) == bogus_large_value) {
492  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
493  }
494  do {
495  ustar = u_star_arr(i,j,k);
496  z0 = Donelan_roughness(ustar);
497  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0);
498  ++iter;
499  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
500  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
501  "Maximum number of MOST iterations reached.");
502 
503  // Populate the stored MOST arrays
504  z0_arr(i,j,k) = z0;
505  olen_arr(i,j,k) = bogus_large_value;
506  t_star_arr(i,j,k) = zero;
507  q_star_arr(i,j,k) = zero;
508  }
509 
510 private:
513  const amrex::Real tol = amrex::Real(1.0e-5);
514  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
515 };
516 
517 
518 /**
519  * Adiabatic with wave-coupled roughness
520  */
522 {
523  /**
524  * Construct adiabatic wave-coupled roughness MOST data.
525  *
526  * @param[in] Tflux prescribed heat flux
527  * @param[in] Qvflux prescribed moisture flux
528  */
530  amrex::Real Qvflux)
531  {
532  mdata.surf_temp_flux = Tflux;
533  mdata.surf_moist_flux = Qvflux;
534  }
535 
536  /**
537  * Iterate MOST state for one surface point.
538  *
539  * @par Calling sequence
540  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
541  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
542  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
543  * eta_arr.
544  */
545  AMREX_GPU_DEVICE
546  AMREX_FORCE_INLINE
547  void
548  iterate_flux (const int& i,
549  const int& j,
550  const int& k,
551  const int& max_iters,
552  const amrex::Array4<const amrex::Real>& zref_arr,
553  const amrex::Array4<amrex::Real>& z0_arr,
554  const amrex::Array4<const amrex::Real>& umm_arr,
555  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
556  const amrex::Array4<const amrex::Real>& /*tvm_arr*/,
557  const amrex::Array4<const amrex::Real>& /*qvm_arr*/,
558  const amrex::Array4<amrex::Real>& u_star_arr,
559  const amrex::Array4<amrex::Real>& /*w_star_arr*/,
560  const amrex::Array4<amrex::Real>& t_star_arr,
561  const amrex::Array4<amrex::Real>& q_star_arr,
562  const amrex::Array4<amrex::Real>& /*t_surf_arr*/,
563  const amrex::Array4<amrex::Real>& /*q_surf_arr*/,
564  const amrex::Array4<amrex::Real>& olen_arr,
565  const amrex::Array4<amrex::Real>& /*pblh_arr*/,
566  const amrex::Array4<amrex::Real>& Hwave_arr,
567  const amrex::Array4<amrex::Real>& Lwave_arr,
568  const amrex::Array4<amrex::Real>& eta_arr) const
569  {
570  int iter = 0;
571  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
572  amrex::Real ustar = zero;
573  amrex::Real zref = zref_arr(i,j,k);
574  amrex::Real z0 = zero;
575  int ie, je;
576  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
577  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
578  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
579  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
580  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
581  do {
582  ustar = u_star_arr(i,j,k);
583  z0 = std::min( std::max(amrex::Real(1200.0) * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), amrex::Real(4.5) )
584  + amrex::Real(0.11) * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
585  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0);
586  ++iter;
587  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
588  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
589  "Maximum number of MOST iterations reached.");
590 
591  // Populate the stored MOST arrays
592  z0_arr(i,j,k) = z0;
593  olen_arr(i,j,k) = bogus_large_value;
594  t_star_arr(i,j,k) = zero;
595  q_star_arr(i,j,k) = zero;
596  }
597 
598 private:
601  const amrex::Real tol = amrex::Real(1.0e-5);
602 #ifdef AMREX_USE_FLOAT
603  const amrex::Real eps = amrex::Real(1e-8);
604 #else
605  const amrex::Real eps = amrex::Real(1e-15);
606 #endif
607  const amrex::Real z0_eps = amrex::Real(1.0e-6);
609  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
610 };
611 
612 
613 /**
614  * Surface flux with constant roughness
615  */
617 {
618  /**
619  * Construct specified-flux, constant-roughness MOST data.
620  *
621  * @param[in] Tflux prescribed heat flux
622  * @param[in] Qvflux prescribed moisture flux
623  * @param[in] cons_qflux whether the moisture flux is specified
624  */
626  amrex::Real Qvflux,
627  bool cons_qflux)
628  {
629  mdata.surf_temp_flux = Tflux;
630  mdata.surf_moist_flux = Qvflux;
631  spec_qflux = cons_qflux;
632  }
633 
634  /**
635  * Iterate MOST state for one surface point.
636  *
637  * @par Calling sequence
638  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
639  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
640  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
641  * eta_arr.
642  */
643  AMREX_GPU_DEVICE
644  AMREX_FORCE_INLINE
645  void
646  iterate_flux (const int& i,
647  const int& j,
648  const int& k,
649  const int& max_iters,
650  const amrex::Array4<const amrex::Real>& zref_arr,
651  const amrex::Array4<const amrex::Real>& z0_arr,
652  const amrex::Array4<const amrex::Real>& umm_arr,
653  const amrex::Array4<const amrex::Real>& tm_arr,
654  const amrex::Array4<const amrex::Real>& tvm_arr,
655  const amrex::Array4<const amrex::Real>& qvm_arr,
656  const amrex::Array4<amrex::Real>& u_star_arr,
657  const amrex::Array4<amrex::Real>& w_star_arr,
658  const amrex::Array4<amrex::Real>& t_star_arr,
659  const amrex::Array4<amrex::Real>& q_star_arr,
660  const amrex::Array4<amrex::Real>& t_surf_arr,
661  const amrex::Array4<amrex::Real>& q_surf_arr,
662  const amrex::Array4<amrex::Real>& olen_arr,
663  const amrex::Array4<amrex::Real>& pblh_arr,
664  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
665  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
666  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
667  {
668  int iter = 0;
669  amrex::Real ustar = zero;
670  amrex::Real wstar = zero;
671  amrex::Real tflux = zero;
672  amrex::Real qflux = zero;
673  amrex::Real zeta = zero;
674  amrex::Real psi_m = zero;
675  amrex::Real psi_h = zero;
676  amrex::Real Olen = zero;
677  amrex::Real zref = zref_arr(i,j,k);
678  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
679  if (u_star_arr(i,j,k) == bogus_large_value) {
680  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
681  } else {
682  Olen = olen_arr(i,j,k);
683  zeta = zref / Olen;
684  psi_m = sfuns.calc_psi_m(zeta);
685  psi_h = sfuns.calc_psi_h(zeta);
686  }
687  do {
688  ustar = u_star_arr(i,j,k);
689  qflux = (spec_qflux) ? mdata.surf_moist_flux :
690  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
691  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
692  tflux = mdata.surf_temp_flux*(one + epsv*qvm_arr(i,j,k)) + qflux*epsv*tm_arr(i,j,k);
693  if (w_star_arr) {
694  // update w* and Umagmean
695  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
696  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
697  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
698  umm = std::max(umm, WSMIN);
699  }
700  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
701  zeta = zref / Olen;
702  psi_m = sfuns.calc_psi_m(zeta);
703  psi_h = sfuns.calc_psi_h(zeta);
704  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0_arr(i,j,k)) - psi_m);
705  ++iter;
706  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
707  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
708  "Maximum number of MOST iterations reached.");
709 
710  // Populate the stored MOST arrays
711  olen_arr(i,j,k) = Olen;
712  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
713  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
714  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
715  if (spec_qflux) {
716  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
717  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
718  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
719  } else {
720  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
721  (std::log(zref / z0_arr(i,j,k)) - psi_h);
722  }
723  }
724 
725 private:
729  const amrex::Real tol = amrex::Real(1.0e-5);
730  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
731 };
732 
733 
734 /**
735  * Surface flux with charnock roughness
736  */
738 {
739  /**
740  * Construct specified-flux, Charnock-roughness MOST data.
741  *
742  * @param[in] Tflux prescribed heat flux
743  * @param[in] Qvflux prescribed moisture flux
744  * @param[in] cnk_a Charnock parameter
745  * @param[in] cnk_visc whether to include viscous roughness
746  * @param[in] cons_qflux whether the moisture flux is specified
747  */
749  amrex::Real Qvflux,
750  amrex::Real cnk_a,
751  bool cnk_visc,
752  bool cons_qflux)
753  {
754  mdata.surf_temp_flux = Tflux;
755  mdata.surf_moist_flux = Qvflux;
756  mdata.Cnk_a = cnk_a;
757  mdata.visc = cnk_visc;
758  spec_qflux = cons_qflux;
759  }
760 
761  /**
762  * Iterate MOST state for one surface point.
763  *
764  * @par Calling sequence
765  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
766  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
767  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
768  * eta_arr.
769  */
770  AMREX_GPU_DEVICE
771  AMREX_FORCE_INLINE
772  void
773  iterate_flux (const int& i,
774  const int& j,
775  const int& k,
776  const int& max_iters,
777  const amrex::Array4<const amrex::Real>& zref_arr,
778  const amrex::Array4<amrex::Real>& z0_arr,
779  const amrex::Array4<const amrex::Real>& umm_arr,
780  const amrex::Array4<const amrex::Real>& tm_arr,
781  const amrex::Array4<const amrex::Real>& tvm_arr,
782  const amrex::Array4<const amrex::Real>& qvm_arr,
783  const amrex::Array4<amrex::Real>& u_star_arr,
784  const amrex::Array4<amrex::Real>& w_star_arr,
785  const amrex::Array4<amrex::Real>& t_star_arr,
786  const amrex::Array4<amrex::Real>& q_star_arr,
787  const amrex::Array4<amrex::Real>& t_surf_arr,
788  const amrex::Array4<amrex::Real>& q_surf_arr,
789  const amrex::Array4<amrex::Real>& olen_arr,
790  const amrex::Array4<amrex::Real>& pblh_arr,
791  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
792  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
793  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
794  {
795  int iter = 0;
796  amrex::Real ustar = zero;
797  amrex::Real wstar = zero;
798  amrex::Real tflux = zero;
799  amrex::Real qflux = zero;
800  amrex::Real z0 = zero;
801  amrex::Real zeta = zero;
802  amrex::Real psi_m = zero;
803  amrex::Real psi_h = zero;
804  amrex::Real Olen = zero;
805  amrex::Real zref = zref_arr(i,j,k);
806  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
807  if (u_star_arr(i,j,k) == bogus_large_value) {
808  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
809  } else {
810  Olen = olen_arr(i,j,k);
811  zeta = zref / Olen;
812  psi_m = sfuns.calc_psi_m(zeta);
813  psi_h = sfuns.calc_psi_h(zeta);
814  }
815  do {
816  ustar = u_star_arr(i,j,k);
817  if (mdata.Cnk_a > 0) {
818  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
819  if (mdata.visc) {
820  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, amrex::Real(0.05));
821  }
822  } else {
823  z0 = COARE3_roughness(zref, umm, ustar);
824  }
825  qflux = (spec_qflux) ? mdata.surf_moist_flux :
826  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
827  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
828  tflux = mdata.surf_temp_flux*(one + epsv*qvm_arr(i,j,k)) + qflux*epsv*tm_arr(i,j,k);
829  if (w_star_arr) {
830  // update w* and Umagmean
831  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
832  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
833  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
834  umm = std::max(umm, WSMIN);
835  }
836  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
837  zeta = zref / Olen;
838  psi_m = sfuns.calc_psi_m(zeta);
839  psi_h = sfuns.calc_psi_h(zeta);
840  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0) - psi_m);
841  ++iter;
842  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
843  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
844  "Maximum number of MOST iterations reached.");
845 
846  // Populate the stored MOST arrays
847  z0_arr(i,j,k) = z0;
848  olen_arr(i,j,k) = Olen;
849  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(zref / z0) - psi_h) /
850  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
851  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
852  if (spec_qflux) {
853  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
854  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
855  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
856  } else {
857  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
858  (std::log(zref / z0_arr(i,j,k)) - psi_h);
859  }
860  }
861 
862 private:
866  const amrex::Real tol = amrex::Real(1.0e-5);
867  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
868 };
869 
870 
871 /**
872  * Surface flux with modified charnock roughness
873  */
875 {
876  /**
877  * Construct specified-flux, modified-Charnock MOST data.
878  *
879  * @param[in] Tflux prescribed heat flux
880  * @param[in] Qvflux prescribed moisture flux
881  * @param[in] depth water depth for modified Charnock roughness
882  * @param[in] cons_qflux whether the moisture flux is specified
883  */
885  amrex::Real Qvflux,
886  amrex::Real depth,
887  bool cons_qflux)
888  {
889  mdata.surf_temp_flux = Tflux;
890  mdata.surf_moist_flux = Qvflux;
891  mdata.Cnk_d = depth;
892  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
893  spec_qflux = cons_qflux;
894  }
895 
896  /**
897  * Iterate MOST state for one surface point.
898  *
899  * @par Calling sequence
900  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
901  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
902  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
903  * eta_arr.
904  */
905  AMREX_GPU_DEVICE
906  AMREX_FORCE_INLINE
907  void
908  iterate_flux (const int& i,
909  const int& j,
910  const int& k,
911  const int& max_iters,
912  const amrex::Array4<const amrex::Real>& zref_arr,
913  const amrex::Array4<amrex::Real>& z0_arr,
914  const amrex::Array4<const amrex::Real>& umm_arr,
915  const amrex::Array4<const amrex::Real>& tm_arr,
916  const amrex::Array4<const amrex::Real>& tvm_arr,
917  const amrex::Array4<const amrex::Real>& qvm_arr,
918  const amrex::Array4<amrex::Real>& u_star_arr,
919  const amrex::Array4<amrex::Real>& w_star_arr,
920  const amrex::Array4<amrex::Real>& t_star_arr,
921  const amrex::Array4<amrex::Real>& q_star_arr,
922  const amrex::Array4<amrex::Real>& t_surf_arr,
923  const amrex::Array4<amrex::Real>& q_surf_arr,
924  const amrex::Array4<amrex::Real>& olen_arr,
925  const amrex::Array4<amrex::Real>& pblh_arr,
926  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
927  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
928  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
929  {
930  int iter = 0;
931  amrex::Real ustar = zero;
932  amrex::Real wstar = zero;
933  amrex::Real tflux = zero;
934  amrex::Real qflux = zero;
935  amrex::Real z0 = zero;
936  amrex::Real zeta = zero;
937  amrex::Real psi_m = zero;
938  amrex::Real psi_h = zero;
939  amrex::Real Olen = zero;
940  amrex::Real zref = zref_arr(i,j,k);
941  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
942  if (u_star_arr(i,j,k) == bogus_large_value) {
943  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
944  } else {
945  Olen = olen_arr(i,j,k);
946  zeta = zref / Olen;
947  psi_m = sfuns.calc_psi_m(zeta);
948  psi_h = sfuns.calc_psi_h(zeta);
949  }
950  do {
951  ustar = u_star_arr(i,j,k);
952  z0 = std::exp( (amrex::Real(2.7)*ustar - amrex::Real(1.8)/mdata.Cnk_b) / (ustar + amrex::Real(0.17)/mdata.Cnk_b) );
953  qflux = (spec_qflux) ? mdata.surf_moist_flux :
954  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
955  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
956  tflux = mdata.surf_temp_flux*(one + epsv*qvm_arr(i,j,k)) + qflux*epsv*tm_arr(i,j,k);
957  if (w_star_arr) {
958  // update w* and Umagmean
959  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
960  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
961  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
962  umm = std::max(umm, WSMIN);
963  }
964  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
965  zeta = zref / Olen;
966  psi_m = sfuns.calc_psi_m(zeta);
967  psi_h = sfuns.calc_psi_h(zeta);
968  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0) - psi_m);
969  ++iter;
970  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
971  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
972  "Maximum number of MOST iterations reached.");
973 
974  // Populate the stored MOST arrays
975  z0_arr(i,j,k) = z0;
976  olen_arr(i,j,k) = Olen;
977  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(zref / z0) - psi_h) /
978  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
979  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
980  if (spec_qflux) {
981  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
982  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
983  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
984  } else {
985  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
986  (std::log(zref / z0_arr(i,j,k)) - psi_h);
987  }
988  }
989 
990 private:
994  const amrex::Real tol = amrex::Real(1.0e-5);
995  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
996 };
997 
998 
999 /**
1000  * Surface flux with donelan roughness
1001  */
1003 {
1004  /**
1005  * Construct specified-flux, Donelan-roughness MOST data.
1006  *
1007  * @param[in] Tflux prescribed heat flux
1008  * @param[in] Qvflux prescribed moisture flux
1009  * @param[in] cons_qflux whether the moisture flux is specified
1010  */
1012  amrex::Real Qvflux,
1013  bool cons_qflux)
1014  {
1015  mdata.surf_temp_flux = Tflux;
1016  mdata.surf_moist_flux = Qvflux;
1017  spec_qflux = cons_qflux;
1018  }
1019 
1020  /**
1021  * Iterate MOST state for one surface point.
1022  *
1023  * @par Calling sequence
1024  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1025  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1026  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1027  * eta_arr.
1028  */
1029  AMREX_GPU_DEVICE
1030  AMREX_FORCE_INLINE
1031  void
1032  iterate_flux (const int& i,
1033  const int& j,
1034  const int& k,
1035  const int& max_iters,
1036  const amrex::Array4<const amrex::Real>& zref_arr,
1037  const amrex::Array4<amrex::Real>& z0_arr,
1038  const amrex::Array4<const amrex::Real>& umm_arr,
1039  const amrex::Array4<const amrex::Real>& tm_arr,
1040  const amrex::Array4<const amrex::Real>& tvm_arr,
1041  const amrex::Array4<const amrex::Real>& qvm_arr,
1042  const amrex::Array4<amrex::Real>& u_star_arr,
1043  const amrex::Array4<amrex::Real>& w_star_arr,
1044  const amrex::Array4<amrex::Real>& t_star_arr,
1045  const amrex::Array4<amrex::Real>& q_star_arr,
1046  const amrex::Array4<amrex::Real>& t_surf_arr,
1047  const amrex::Array4<amrex::Real>& q_surf_arr,
1048  const amrex::Array4<amrex::Real>& olen_arr,
1049  const amrex::Array4<amrex::Real>& pblh_arr,
1050  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1051  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1052  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1053  {
1054  int iter = 0;
1055  amrex::Real ustar = zero;
1056  amrex::Real wstar = zero;
1057  amrex::Real tflux = zero;
1058  amrex::Real qflux = zero;
1059  amrex::Real z0 = zero;
1060  amrex::Real zeta = zero;
1061  amrex::Real psi_m = zero;
1062  amrex::Real psi_h = zero;
1063  amrex::Real Olen = zero;
1064  amrex::Real zref = zref_arr(i,j,k);
1065  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1066  if (u_star_arr(i,j,k) == bogus_large_value) {
1067  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
1068  } else {
1069  Olen = olen_arr(i,j,k);
1070  zeta = zref / Olen;
1071  psi_m = sfuns.calc_psi_m(zeta);
1072  psi_h = sfuns.calc_psi_h(zeta);
1073  }
1074  do {
1075  ustar = u_star_arr(i,j,k);
1076  z0 = Donelan_roughness(ustar);
1077  qflux = (spec_qflux) ? mdata.surf_moist_flux :
1078  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
1079  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
1080  tflux = mdata.surf_temp_flux*(one + epsv*qvm_arr(i,j,k)) + qflux*epsv*tm_arr(i,j,k);
1081  if (w_star_arr) {
1082  // update w* and Umagmean
1083  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1084  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1085  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1086  umm = std::max(umm, WSMIN);
1087  }
1088  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1089  zeta = zref / Olen;
1090  psi_m = sfuns.calc_psi_m(zeta);
1091  psi_h = sfuns.calc_psi_h(zeta);
1092  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0) - psi_m);
1093  ++iter;
1094  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1095  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
1096  "Maximum number of MOST iterations reached.");
1097 
1098  // Populate the stored MOST arrays
1099  z0_arr(i,j,k) = z0;
1100  olen_arr(i,j,k) = Olen;
1101  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(zref / z0) - psi_h) /
1102  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
1103  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
1104  if (spec_qflux) {
1105  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
1106  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
1107  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
1108  } else {
1109  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
1110  (std::log(zref / z0_arr(i,j,k)) - psi_h);
1111  }
1112  }
1113 
1114 private:
1118  const amrex::Real tol = amrex::Real(1.0e-5);
1119  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
1120 };
1121 
1122 
1123 /**
1124  * Surface flux with wave-coupled roughness
1125  */
1127 {
1128  /**
1129  * Construct specified-flux, wave-coupled roughness MOST data.
1130  *
1131  * @param[in] Tflux prescribed heat flux
1132  * @param[in] Qvflux prescribed moisture flux
1133  * @param[in] cons_qflux whether the moisture flux is specified
1134  */
1136  amrex::Real Qvflux,
1137  bool cons_qflux)
1138  {
1139  mdata.surf_temp_flux = Tflux;
1140  mdata.surf_moist_flux = Qvflux;
1141  spec_qflux = cons_qflux;
1142  }
1143 
1144  /**
1145  * Iterate MOST state for one surface point.
1146  *
1147  * @par Calling sequence
1148  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1149  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1150  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1151  * eta_arr.
1152  */
1153  AMREX_GPU_DEVICE
1154  AMREX_FORCE_INLINE
1155  void
1156  iterate_flux (const int& i,
1157  const int& j,
1158  const int& k,
1159  const int& max_iters,
1160  const amrex::Array4<const amrex::Real>& zref_arr,
1161  const amrex::Array4<amrex::Real>& z0_arr,
1162  const amrex::Array4<const amrex::Real>& umm_arr,
1163  const amrex::Array4<const amrex::Real>& tm_arr,
1164  const amrex::Array4<const amrex::Real>& tvm_arr,
1165  const amrex::Array4<const amrex::Real>& qvm_arr,
1166  const amrex::Array4<amrex::Real>& u_star_arr,
1167  const amrex::Array4<amrex::Real>& w_star_arr,
1168  const amrex::Array4<amrex::Real>& t_star_arr,
1169  const amrex::Array4<amrex::Real>& q_star_arr,
1170  const amrex::Array4<amrex::Real>& t_surf_arr,
1171  const amrex::Array4<amrex::Real>& q_surf_arr,
1172  const amrex::Array4<amrex::Real>& olen_arr,
1173  const amrex::Array4<amrex::Real>& pblh_arr,
1174  const amrex::Array4<amrex::Real>& Hwave_arr,
1175  const amrex::Array4<amrex::Real>& Lwave_arr,
1176  const amrex::Array4<amrex::Real>& eta_arr) const
1177  {
1178  int iter = 0;
1179  amrex::Real ustar = zero;
1180  amrex::Real wstar = zero;
1181  amrex::Real tflux = zero;
1182  amrex::Real qflux = zero;
1183  amrex::Real z0 = zero;
1184  amrex::Real zeta = zero;
1185  amrex::Real psi_m = zero;
1186  amrex::Real psi_h = zero;
1187  amrex::Real Olen = zero;
1188  amrex::Real zref = zref_arr(i,j,k);
1189  int ie, je;
1190  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1191  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1192  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1193  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1194  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1195  if (u_star_arr(i,j,k) == bogus_large_value) {
1196  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
1197  } else {
1198  Olen = olen_arr(i,j,k);
1199  zeta = zref / Olen;
1200  psi_m = sfuns.calc_psi_m(zeta);
1201  psi_h = sfuns.calc_psi_h(zeta);
1202  }
1203  do {
1204  ustar = u_star_arr(i,j,k);
1205  z0 = std::min( std::max(amrex::Real(1200.0) * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), amrex::Real(4.5) )
1206  + amrex::Real(0.11) * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
1207  qflux = (spec_qflux) ? mdata.surf_moist_flux :
1208  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
1209  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
1210  tflux = mdata.surf_temp_flux*(one + epsv*qvm_arr(i,j,k)) + qflux*epsv*tm_arr(i,j,k);
1211  if (w_star_arr) {
1212  // update w* and Umagmean
1213  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1214  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1215  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1216  umm = std::max(umm, WSMIN);
1217  }
1218  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1219  zeta = zref / Olen;
1220  psi_m = sfuns.calc_psi_m(zeta);
1221  psi_h = sfuns.calc_psi_h(zeta);
1222  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0) - psi_m);
1223  ++iter;
1224  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1225  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
1226  "Maximum number of MOST iterations reached.");
1227 
1228  // Populate the stored MOST arrays
1229  z0_arr(i,j,k) = z0;
1230  olen_arr(i,j,k) = Olen;
1231  t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(zref / z0) - psi_h) /
1232  (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k);
1233  t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k);
1234  if (spec_qflux) {
1235  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
1236  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
1237  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
1238  } else {
1239  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
1240  (std::log(zref / z0_arr(i,j,k)) - psi_h);
1241  }
1242  }
1243 
1244 private:
1248  const amrex::Real tol = amrex::Real(1.0e-5);
1249 #ifdef AMREX_USE_FLOAT
1250  const amrex::Real eps = amrex::Real(1e-8);
1251 #else
1252  const amrex::Real eps = amrex::Real(1e-15);
1253 #endif
1254  const amrex::Real z0_eps = amrex::Real(1.0e-6);
1256  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
1257 };
1258 
1259 
1260 /**
1261  * Surface temperature with constant roughness
1262  */
1264 {
1265  /**
1266  * Construct specified-temperature, constant-roughness MOST data.
1267  *
1268  * @param[in] Tflux prescribed heat flux
1269  * @param[in] Qvflux prescribed moisture flux
1270  * @param[in] cons_qflux whether the moisture flux is specified
1271  */
1273  amrex::Real Qvflux,
1274  bool cons_qflux)
1275  {
1276  mdata.surf_temp_flux = Tflux;
1277  mdata.surf_moist_flux = Qvflux;
1278  spec_qflux = cons_qflux;
1279  }
1280 
1281  /**
1282  * Iterate MOST state for one surface point.
1283  *
1284  * @par Calling sequence
1285  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1286  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1287  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1288  * eta_arr.
1289  */
1290  AMREX_GPU_DEVICE
1291  AMREX_FORCE_INLINE
1292  void
1293  iterate_flux (const int& i,
1294  const int& j,
1295  const int& k,
1296  const int& max_iters,
1297  const amrex::Array4<const amrex::Real>& zref_arr,
1298  const amrex::Array4<const amrex::Real>& z0_arr,
1299  const amrex::Array4<const amrex::Real>& umm_arr,
1300  const amrex::Array4<const amrex::Real>& tm_arr,
1301  const amrex::Array4<const amrex::Real>& tvm_arr,
1302  const amrex::Array4<const amrex::Real>& qvm_arr,
1303  const amrex::Array4<amrex::Real>& u_star_arr,
1304  const amrex::Array4<amrex::Real>& w_star_arr,
1305  const amrex::Array4<amrex::Real>& t_star_arr,
1306  const amrex::Array4<amrex::Real>& q_star_arr,
1307  const amrex::Array4<amrex::Real>& t_surf_arr,
1308  const amrex::Array4<amrex::Real>& q_surf_arr,
1309  const amrex::Array4<amrex::Real>& olen_arr,
1310  const amrex::Array4<amrex::Real>& pblh_arr,
1311  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1312  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1313  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1314  {
1315  amrex::Real Rib = zero;
1316  amrex::Real zeta = zero;
1317  amrex::Real zeta_old = zero;
1318  amrex::Real psi_m = zero;
1319  amrex::Real psi_h = zero;
1320  amrex::Real num = zero;
1321  amrex::Real den = zero;
1322  amrex::Real zref = zref_arr(i,j,k);
1323  amrex::Real z0 = z0_arr(i,j,k);
1324  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1325  amrex::Real C = std::log(zref / z0);
1326 
1327  // First iteration we assume neutral (L -> inf)
1328  if (u_star_arr(i,j,k) == bogus_large_value) { olen_arr(i,j,k) = amrex::Real(1.0e3); }
1329  zeta = zref / olen_arr(i,j,k);
1330 
1331  // Water vapor in atmos and surface
1332  amrex::Real qv_s, qv_a;
1333  if (q_surf_arr(i,j,k) > zero) {
1334  qv_s = q_surf_arr(i,j,k);
1335  } else {
1336  // First iteration and no qv_surf was specified
1337  // Use mean since there will be no flux
1338  qv_s = qvm_arr(i,j,k);
1339  }
1340  qv_a = qvm_arr(i,j,k);
1341 
1342  // update w* and Umagmean from Beljaars (1995)
1343  if (w_star_arr) {
1344  // NOTE: Thv flux is lagged, similar to WRF
1345  psi_m = sfuns.calc_psi_m2(zeta);
1346  psi_h = sfuns.calc_psi_h2(zeta);
1347  amrex::Real ustar = mdata.kappa * umm / (C - psi_m);
1348  amrex::Real tstar = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / (C - psi_h);
1350  -ustar * mdata.kappa * (qv_a - qv_s) / (C - psi_h);
1351  amrex::Real tflux = -ustar*tstar*(one + epsv*qvm_arr(i,j,k)) + epsv*tm_arr(i,j,k)*qflux;
1352  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1353  amrex::Real wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1354  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1355  umm = std::max(umm, WSMIN);
1356  }
1357 
1358  // Bulk Richardson number w/ moisture
1359  amrex::Real thv_s = t_surf_arr(i,j,k) * (one + epsv*qv_s);
1360  amrex::Real thv_a = tm_arr(i,j,k) * (one + epsv*qv_a);
1361  Rib = ( (mdata.gravity * zref) / tm_arr(i,j,k) ) *
1362  ( (thv_a - thv_s) / (umm * umm) );
1363  Rib = std::min(std::max(Rib,-amrex::Real(4.0)),amrex::Real(4.0));
1364 
1365  // Fixed point iteration on zeta
1366  int iter = 0;
1367  do {
1368  // Transfer curr to old
1369  zeta_old = zeta;
1370 
1371  // Stability functions
1372  psi_m = sfuns.calc_psi_m2(zeta_old);
1373  psi_h = sfuns.calc_psi_h2(zeta_old);
1374 
1375  // Limiting
1376  num = std::max(C - psi_m, amrex::Real(1.0));
1377  den = std::max(C - psi_h, amrex::Real(1.0));
1378 
1379  // Update with under relaxation
1380  zeta = (one - alpha) * zeta_old + alpha * Rib * num * num / den;
1381 
1382  ++iter;
1383  } while ( (std::abs(zeta - zeta_old) > tol) && (iter <= max_iters) );
1384  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
1385  "Maximum number of MOST iterations reached.");
1386 
1387  // Populate the stored MOST arrays
1388  olen_arr(i,j,k) = zref / zeta;
1389  u_star_arr(i,j,k) = mdata.kappa * umm / (C - psi_m);
1390  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / (C - psi_h);
1391  if (spec_qflux) {
1392  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (C - psi_h) /
1393  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
1394  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
1395  } else {
1396  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) / (C - psi_h);
1397  }
1398  }
1399 
1400 private:
1404  const amrex::Real tol = amrex::Real(1.0e-3);
1406  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
1407 };
1408 
1409 
1410 /**
1411  * Surface temperature with charnock roughness
1412  */
1414 {
1415  /**
1416  * Construct specified-temperature, Charnock-roughness MOST data.
1417  *
1418  * @param[in] Tflux prescribed heat flux
1419  * @param[in] Qvflux prescribed moisture flux
1420  * @param[in] cnk_a Charnock parameter
1421  * @param[in] cnk_visc whether to include viscous roughness
1422  * @param[in] cons_qflux whether the moisture flux is specified
1423  */
1425  amrex::Real Qvflux,
1426  amrex::Real cnk_a,
1427  bool cnk_visc,
1428  bool cons_qflux)
1429  {
1430  mdata.surf_temp_flux = Tflux;
1431  mdata.surf_moist_flux = Qvflux;
1432  mdata.Cnk_a = cnk_a;
1433  mdata.visc = cnk_visc;
1434  spec_qflux = cons_qflux;
1435  }
1436 
1437  /**
1438  * Iterate MOST state for one surface point.
1439  *
1440  * @par Calling sequence
1441  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1442  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1443  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1444  * eta_arr.
1445  */
1446  AMREX_GPU_DEVICE
1447  AMREX_FORCE_INLINE
1448  void
1449  iterate_flux (const int& i,
1450  const int& j,
1451  const int& k,
1452  const int& max_iters,
1453  const amrex::Array4<const amrex::Real>& zref_arr,
1454  const amrex::Array4<amrex::Real>& z0_arr,
1455  const amrex::Array4<const amrex::Real>& umm_arr,
1456  const amrex::Array4<const amrex::Real>& tm_arr,
1457  const amrex::Array4<const amrex::Real>& tvm_arr,
1458  const amrex::Array4<const amrex::Real>& qvm_arr,
1459  const amrex::Array4<amrex::Real>& u_star_arr,
1460  const amrex::Array4<amrex::Real>& w_star_arr,
1461  const amrex::Array4<amrex::Real>& t_star_arr,
1462  const amrex::Array4<amrex::Real>& q_star_arr,
1463  const amrex::Array4<amrex::Real>& t_surf_arr,
1464  const amrex::Array4<amrex::Real>& q_surf_arr,
1465  const amrex::Array4<amrex::Real>& olen_arr,
1466  const amrex::Array4<amrex::Real>& pblh_arr,
1467  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1468  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1469  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1470  {
1471  amrex::Real Rib = zero;
1472  amrex::Real zeta = zero;
1473  amrex::Real zeta_old = zero;
1474  amrex::Real psi_m = zero;
1475  amrex::Real psi_h = zero;
1476  amrex::Real num = zero;
1477  amrex::Real den = zero;
1478  amrex::Real z0 = z0_arr(i,j,k);
1479  amrex::Real z0_old = z0;
1480  amrex::Real ustar = zero;
1481  amrex::Real zref = zref_arr(i,j,k);
1482  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1483  amrex::Real C = std::log(zref / z0);
1484 
1485  // First iteration we assume neutral (L -> inf)
1486  if (u_star_arr(i,j,k) == bogus_large_value) {
1487  olen_arr(i,j,k) = amrex::Real(1.0e3);
1488  }
1489  zeta = zref / olen_arr(i,j,k);
1490 
1491  // Water vapor in atmos and surface
1492  amrex::Real qv_s, qv_a;
1493  if (q_surf_arr(i,j,k) > zero) {
1494  qv_s = q_surf_arr(i,j,k);
1495  } else {
1496  // First iteration and no qv_surf was specified
1497  // Use mean since there will be no flux
1498  qv_s = qvm_arr(i,j,k);
1499  }
1500  qv_a = qvm_arr(i,j,k);
1501 
1502  // update w* and Umagmean from Beljaars (1995)
1503  if (w_star_arr) {
1504  // NOTE: Thv flux is lagged, similar to WRF
1505  psi_m = sfuns.calc_psi_m2(zeta);
1506  psi_h = sfuns.calc_psi_h2(zeta);
1507  ustar = mdata.kappa * umm / (C - psi_m);
1508  amrex::Real tstar = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / (C - psi_h);
1510  -ustar * mdata.kappa * (qv_a - qv_s) / (C - psi_h);
1511  amrex::Real tflux = -ustar*tstar*(one + epsv*qvm_arr(i,j,k)) + epsv*tm_arr(i,j,k)*qflux;
1512  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1513  amrex::Real wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1514  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1515  umm = std::max(umm, WSMIN);
1516  }
1517 
1518  // Bulk Richardson number w/ moisture
1519  amrex::Real thv_s = t_surf_arr(i,j,k) * (one + epsv*qv_s);
1520  amrex::Real thv_a = tm_arr(i,j,k) * (one + epsv*qv_a);
1521  Rib = ( (mdata.gravity * zref) / tm_arr(i,j,k) ) *
1522  ( (thv_a - thv_s) / (umm * umm) );
1523  Rib = std::min(std::max(Rib,-amrex::Real(4.0)),amrex::Real(4.0));
1524 
1525  // Fixed point iteration on zeta
1526  int iter = 0;
1527  do {
1528  // Transfer curr to old
1529  zeta_old = zeta;
1530 
1531  // Stability functions
1532  psi_m = sfuns.calc_psi_m2(zeta_old);
1533  psi_h = sfuns.calc_psi_h2(zeta_old);
1534 
1535  // Fixed point iteration on roughness
1536  int iter_z = 0;
1537  do {
1538  // Transfer curr to old
1539  z0_old = z0;
1540 
1541  // Update
1542  C = std::log(zref / z0_old);
1543  ustar = mdata.kappa * umm / (C - psi_m);
1544  if (mdata.Cnk_a > 0) {
1545  z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar;
1546  if (mdata.visc) {
1547  z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, amrex::Real(0.05));
1548  }
1549  } else {
1550  z0 = COARE3_roughness(zref, umm, ustar);
1551  }
1552 
1553  ++iter_z;
1554  } while ( (std::abs(z0 - z0_old) > tol_z) && (iter_z <= max_iters) );
1555  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter_z < max_iters,
1556  "Maximum number of MOST roughness iterations reached.");
1557  C = std::log(zref / z0);
1558 
1559  // Limiting
1560  num = std::max(C - psi_m, amrex::Real(1.0));
1561  den = std::max(C - psi_h, amrex::Real(1.0));
1562 
1563  // Update with under relaxation
1564  zeta = (amrex::Real(1.0) - alpha) * zeta_old + alpha * Rib * num * num / den;
1565 
1566  ++iter;
1567  } while ( (std::abs(zeta - zeta_old) > tol) && (iter <= max_iters) );
1568  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
1569  "Maximum number of MOST iterations reached.");
1570 
1571  // Populate the stored MOST arrays
1572  z0_arr(i,j,k) = z0;
1573  olen_arr(i,j,k) = zref / zeta;
1574  u_star_arr(i,j,k) = mdata.kappa * umm / (C - psi_m);
1575  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / (C - psi_h);
1576  if (spec_qflux) {
1577  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (C - psi_h) /
1578  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
1579  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
1580  } else {
1581  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) / (C - psi_h);
1582  }
1583  }
1584 
1585 private:
1589  const amrex::Real tol = amrex::Real(1.0e-3);
1590 #ifdef AMREX_USE_FLOAT
1591  const amrex::Real tol_z = amrex::Real(1.0e-6);
1592 #else
1593  const amrex::Real tol_z = amrex::Real(1.0e-10);
1594 #endif
1596  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
1597 };
1598 
1599 
1600 /**
1601  * Surface temperature with modified charnock roughness
1602  */
1604 {
1605  /**
1606  * Construct specified-temperature, modified-Charnock MOST data.
1607  *
1608  * @param[in] Tflux prescribed heat flux
1609  * @param[in] Qvflux prescribed moisture flux
1610  * @param[in] depth water depth for modified Charnock roughness
1611  * @param[in] cons_qflux whether the moisture flux is specified
1612  */
1614  amrex::Real Qvflux,
1615  amrex::Real depth,
1616  bool cons_qflux)
1617  {
1618  mdata.surf_temp_flux = Tflux;
1619  mdata.surf_moist_flux = Qvflux;
1620  mdata.Cnk_d = depth;
1621  mdata.Cnk_b = mdata.Cnk_b1 * std::log(mdata.Cnk_b2 / mdata.Cnk_d);
1622  spec_qflux = cons_qflux;
1623  }
1624 
1625  /**
1626  * Iterate MOST state for one surface point.
1627  *
1628  * @par Calling sequence
1629  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1630  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1631  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1632  * eta_arr.
1633  */
1634  AMREX_GPU_DEVICE
1635  AMREX_FORCE_INLINE
1636  void
1637  iterate_flux (const int& i,
1638  const int& j,
1639  const int& k,
1640  const int& max_iters,
1641  const amrex::Array4<const amrex::Real>& zref_arr,
1642  const amrex::Array4<amrex::Real>& z0_arr,
1643  const amrex::Array4<const amrex::Real>& umm_arr,
1644  const amrex::Array4<const amrex::Real>& tm_arr,
1645  const amrex::Array4<const amrex::Real>& tvm_arr,
1646  const amrex::Array4<const amrex::Real>& qvm_arr,
1647  const amrex::Array4<amrex::Real>& u_star_arr,
1648  const amrex::Array4<amrex::Real>& w_star_arr,
1649  const amrex::Array4<amrex::Real>& t_star_arr,
1650  const amrex::Array4<amrex::Real>& q_star_arr,
1651  const amrex::Array4<amrex::Real>& t_surf_arr,
1652  const amrex::Array4<amrex::Real>& q_surf_arr,
1653  const amrex::Array4<amrex::Real>& olen_arr,
1654  const amrex::Array4<amrex::Real>& pblh_arr,
1655  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1656  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1657  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1658  {
1659  amrex::Real Rib = zero;
1660  amrex::Real zeta = zero;
1661  amrex::Real zeta_old = zero;
1662  amrex::Real psi_m = zero;
1663  amrex::Real psi_h = zero;
1664  amrex::Real num = zero;
1665  amrex::Real den = zero;
1666  amrex::Real z0 = z0_arr(i,j,k);
1667  amrex::Real z0_old = z0;
1668  amrex::Real ustar = zero;
1669  amrex::Real zref = zref_arr(i,j,k);
1670  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1671  amrex::Real C = std::log(zref / z0);
1672 
1673  // First iteration we assume neutral (L -> inf)
1674  if (u_star_arr(i,j,k) == bogus_large_value) {
1675  olen_arr(i,j,k) = amrex::Real(1000);
1676  }
1677  zeta = zref / olen_arr(i,j,k);
1678 
1679  // Water vapor in atmos and surface
1680  amrex::Real qv_s, qv_a;
1681  if (q_surf_arr(i,j,k) > zero) {
1682  qv_s = q_surf_arr(i,j,k);
1683  } else {
1684  // First iteration and no qv_surf was specified
1685  // Use mean since there will be no flux
1686  qv_s = qvm_arr(i,j,k);
1687  }
1688  qv_a = qvm_arr(i,j,k);
1689 
1690  // update w* and Umagmean from Beljaars (1995)
1691  if (w_star_arr) {
1692  // NOTE: Thv flux is lagged, similar to WRF
1693  psi_m = sfuns.calc_psi_m2(zeta);
1694  psi_h = sfuns.calc_psi_h2(zeta);
1695  ustar = mdata.kappa * umm / (C - psi_m);
1696  amrex::Real tstar = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / (C - psi_h);
1698  -ustar * mdata.kappa * (qv_a - qv_s) / (C - psi_h);
1699  amrex::Real tflux = -ustar*tstar*(one + epsv*qvm_arr(i,j,k)) + epsv*tm_arr(i,j,k)*qflux;
1700  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1701  amrex::Real wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1702  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1703  umm = std::max(umm, WSMIN);
1704  }
1705 
1706  // Bulk Richardson number w/ moisture
1707  amrex::Real thv_s = t_surf_arr(i,j,k) * (one + epsv*qv_s);
1708  amrex::Real thv_a = tm_arr(i,j,k) * (one + epsv*qv_a);
1709  Rib = ( (mdata.gravity * zref) / tm_arr(i,j,k) ) *
1710  ( (thv_a - thv_s) / (umm * umm) );
1711  Rib = std::min(std::max(Rib,-amrex::Real(4.0)),amrex::Real(4.0));
1712 
1713  // Fixed point iteration on zeta
1714  int iter = 0;
1715  do {
1716  // Transfer curr to old
1717  zeta_old = zeta;
1718 
1719  // Stability functions
1720  psi_m = sfuns.calc_psi_m2(zeta_old);
1721  psi_h = sfuns.calc_psi_h2(zeta_old);
1722 
1723  // Fixed point iteration on roughness
1724  int iter_z = 0;
1725  do {
1726  // Transfer curr to old
1727  z0_old = z0;
1728 
1729  // Update
1730  C = std::log(zref / z0_old);
1731  ustar = mdata.kappa * umm / (C - psi_m);
1732  z0 = std::exp( (amrex::Real(2.7)*ustar - amrex::Real(1.8)/mdata.Cnk_b) / (ustar + amrex::Real(0.17)/mdata.Cnk_b) );
1733 
1734  ++iter_z;
1735  } while ( (std::abs(z0 - z0_old) > tol_z) && (iter_z <= max_iters) );
1736  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter_z < max_iters,
1737  "Maximum number of MOST roughness iterations reached.");
1738  C = std::log(zref / z0);
1739 
1740  // Limiting
1741  num = std::max(C - psi_m, amrex::Real(1.0));
1742  den = std::max(C - psi_h, amrex::Real(1.0));
1743 
1744  // Update with under relaxation
1745  zeta = (amrex::Real(1.0) - alpha) * zeta_old + alpha * Rib * num * num / den;
1746 
1747  ++iter;
1748  } while ( (std::abs(zeta - zeta_old) > tol) && (iter <= max_iters) );
1749  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
1750  "Maximum number of MOST iterations reached.");
1751 
1752  // Populate the stored MOST arrays
1753  z0_arr(i,j,k) = z0;
1754  olen_arr(i,j,k) = zref / zeta;
1755  u_star_arr(i,j,k) = mdata.kappa * umm / (C - psi_m);
1756  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / (C - psi_h);
1757  if (spec_qflux) {
1758  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (C - psi_h) /
1759  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
1760  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
1761  } else {
1762  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) / (C - psi_h);
1763  }
1764  }
1765 
1766 private:
1770  const amrex::Real tol = amrex::Real(1.0e-3);
1771 #ifdef AMREX_USE_FLOAT
1772  const amrex::Real tol_z = amrex::Real(1.0e-6);
1773 #else
1774  const amrex::Real tol_z = amrex::Real(1.0e-10);
1775 #endif
1777  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
1778 };
1779 
1780 
1781 /**
1782  * Surface temperature with donelan roughness
1783  */
1785 {
1786  /**
1787  * Construct specified-temperature, Donelan-roughness MOST data.
1788  *
1789  * @param[in] Tflux prescribed heat flux
1790  * @param[in] Qvflux prescribed moisture flux
1791  * @param[in] cons_qflux whether the moisture flux is specified
1792  */
1794  amrex::Real Qvflux,
1795  bool cons_qflux)
1796  {
1797  mdata.surf_temp_flux = Tflux;
1798  mdata.surf_moist_flux = Qvflux;
1799  spec_qflux = cons_qflux;
1800  }
1801 
1802  /**
1803  * Iterate MOST state for one surface point.
1804  *
1805  * @par Calling sequence
1806  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1807  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1808  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1809  * eta_arr.
1810  */
1811  AMREX_GPU_DEVICE
1812  AMREX_FORCE_INLINE
1813  void
1814  iterate_flux (const int& i,
1815  const int& j,
1816  const int& k,
1817  const int& max_iters,
1818  const amrex::Array4<const amrex::Real>& zref_arr,
1819  const amrex::Array4<amrex::Real>& z0_arr,
1820  const amrex::Array4<const amrex::Real>& umm_arr,
1821  const amrex::Array4<const amrex::Real>& tm_arr,
1822  const amrex::Array4<const amrex::Real>& tvm_arr,
1823  const amrex::Array4<const amrex::Real>& qvm_arr,
1824  const amrex::Array4<amrex::Real>& u_star_arr,
1825  const amrex::Array4<amrex::Real>& w_star_arr,
1826  const amrex::Array4<amrex::Real>& t_star_arr,
1827  const amrex::Array4<amrex::Real>& q_star_arr,
1828  const amrex::Array4<amrex::Real>& t_surf_arr,
1829  const amrex::Array4<amrex::Real>& q_surf_arr,
1830  const amrex::Array4<amrex::Real>& olen_arr,
1831  const amrex::Array4<amrex::Real>& pblh_arr,
1832  const amrex::Array4<amrex::Real>& /*Hwave_arr*/,
1833  const amrex::Array4<amrex::Real>& /*Lwave_arr*/,
1834  const amrex::Array4<amrex::Real>& /*eta_arr*/) const
1835  {
1836  int iter = 0;
1837  amrex::Real ustar = zero;
1838  amrex::Real wstar = zero;
1839  amrex::Real z0 = zero;
1840  amrex::Real tflux = zero;
1841  amrex::Real qflux = zero;
1842  amrex::Real zeta = zero;
1843  amrex::Real psi_m = zero;
1844  amrex::Real psi_h = zero;
1845  amrex::Real Olen = zero;
1846  amrex::Real Oleno = zero;
1847  amrex::Real zref = zref_arr(i,j,k);
1848  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1849  if (u_star_arr(i,j,k) == bogus_large_value) {
1850  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
1851  } else {
1852  Olen = olen_arr(i,j,k);
1853  Oleno = Olen;
1854  zeta = zref / Olen;
1855  psi_m = sfuns.calc_psi_m(zeta);
1856  psi_h = sfuns.calc_psi_h(zeta);
1857  }
1858  do {
1859  ustar = u_star_arr(i,j,k);
1860  z0 = Donelan_roughness(ustar);
1861  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
1862  (std::log(zref / z0) - psi_h); // <w'T'>
1863  tflux *= (one + epsv*qvm_arr(i,j,k));
1864  qflux = (spec_qflux) ? mdata.surf_moist_flux :
1865  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
1866  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
1867  tflux += epsv * tm_arr(i,j,k) * qflux; // ~= <w'Tv'>
1868  if (w_star_arr) {
1869  // update w* and Umagmean
1870  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
1871  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
1872  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
1873  umm = std::max(umm, WSMIN);
1874  }
1875  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
1876  if ( (((Olen >= zero) && (Oleno <= zero)) ||
1877  ((Olen <= zero) && (Oleno >= zero))) &&
1878  std::fabs(Olen) + std::fabs(Oleno) < amrex::Real(1.0)) {
1879  Olen = myhalf * (Olen + Oleno);
1880  }
1881  Oleno = Olen;
1882  zeta = zref / Olen;
1883  psi_m = sfuns.calc_psi_m(zeta);
1884  psi_h = sfuns.calc_psi_h(zeta);
1885  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0) - psi_m);
1886  ++iter;
1887  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
1888  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
1889  "Maximum number of MOST iterations reached.");
1890 
1891  // Populate the stored MOST arrays
1892  z0_arr(i,j,k) = z0;
1893  olen_arr(i,j,k) = Olen;
1894  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
1895  (std::log(zref / z0) - psi_h);
1896  if (spec_qflux) {
1897  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
1898  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
1899  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
1900  } else {
1901  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
1902  (std::log(zref / z0_arr(i,j,k)) - psi_h);
1903  }
1904  }
1905 
1906 private:
1910  const amrex::Real tol = amrex::Real(1.0e-5);
1911  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
1912 };
1913 
1914 
1915 /**
1916  * Surface temperature with wave-coupled roughness
1917  */
1919 {
1920  /**
1921  * Construct specified-temperature, wave-coupled roughness MOST data.
1922  *
1923  * @param[in] Tflux prescribed heat flux
1924  * @param[in] Qvflux prescribed moisture flux
1925  * @param[in] cons_qflux whether the moisture flux is specified
1926  */
1928  amrex::Real Qvflux,
1929  bool cons_qflux)
1930  {
1931  mdata.surf_temp_flux = Tflux;
1932  mdata.surf_moist_flux = Qvflux;
1933  spec_qflux = cons_qflux;
1934  }
1935 
1936  /**
1937  * Iterate MOST state for one surface point.
1938  *
1939  * @par Calling sequence
1940  * i, j, k, max_iters, zref_arr, z0_arr, umm_arr, tm_arr, tvm_arr,
1941  * qvm_arr, u_star_arr, w_star_arr, t_star_arr, q_star_arr,
1942  * t_surf_arr, q_surf_arr, olen_arr, pblh_arr, Hwave_arr, Lwave_arr,
1943  * eta_arr.
1944  */
1945  AMREX_GPU_DEVICE
1946  AMREX_FORCE_INLINE
1947  void
1948  iterate_flux (const int& i,
1949  const int& j,
1950  const int& k,
1951  const int& max_iters,
1952  const amrex::Array4<const amrex::Real>& zref_arr,
1953  const amrex::Array4<amrex::Real>& z0_arr,
1954  const amrex::Array4<const amrex::Real>& umm_arr,
1955  const amrex::Array4<const amrex::Real>& tm_arr,
1956  const amrex::Array4<const amrex::Real>& tvm_arr,
1957  const amrex::Array4<const amrex::Real>& qvm_arr,
1958  const amrex::Array4<amrex::Real>& u_star_arr,
1959  const amrex::Array4<amrex::Real>& w_star_arr,
1960  const amrex::Array4<amrex::Real>& t_star_arr,
1961  const amrex::Array4<amrex::Real>& q_star_arr,
1962  const amrex::Array4<amrex::Real>& t_surf_arr,
1963  const amrex::Array4<amrex::Real>& q_surf_arr,
1964  const amrex::Array4<amrex::Real>& olen_arr,
1965  const amrex::Array4<amrex::Real>& pblh_arr,
1966  const amrex::Array4<amrex::Real>& Hwave_arr,
1967  const amrex::Array4<amrex::Real>& Lwave_arr,
1968  const amrex::Array4<amrex::Real>& eta_arr) const
1969  {
1970  int iter = 0;
1971  amrex::Real ustar = zero;
1972  amrex::Real wstar = zero;
1973  amrex::Real z0 = zero;
1974  amrex::Real tflux = zero;
1975  amrex::Real qflux = zero;
1976  amrex::Real zeta = zero;
1977  amrex::Real psi_m = zero;
1978  amrex::Real psi_h = zero;
1979  amrex::Real Olen = zero;
1980  amrex::Real Oleno = zero;
1981  amrex::Real zref = zref_arr(i,j,k);
1982  int ie, je;
1983  ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i;
1984  je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j;
1985  ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
1986  je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
1987  amrex::Real umm = std::max(umm_arr(i,j,k), WSMIN);
1988  if (u_star_arr(i,j,k) == bogus_large_value) {
1989  u_star_arr(i,j,k) = mdata.kappa * umm / std::log(zref / z0_arr(i,j,k));
1990  } else {
1991  Olen = olen_arr(i,j,k);
1992  Oleno = Olen;
1993  zeta = zref / Olen;
1994  psi_m = sfuns.calc_psi_m(zeta);
1995  psi_h = sfuns.calc_psi_h(zeta);
1996  }
1997  do {
1998  ustar = u_star_arr(i,j,k);
1999  z0 = std::min( std::max(amrex::Real(1200.0) * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), amrex::Real(4.5) )
2000  + amrex::Real(0.11) * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
2001  tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
2002  (std::log(zref / z0) - psi_h); // <w'T'>
2003  tflux *= (one + epsv*qvm_arr(i,j,k));
2004  qflux = (spec_qflux) ? mdata.surf_moist_flux :
2005  -(qvm_arr(i,j,k) - q_surf_arr(i,j,k)) * ustar * mdata.kappa /
2006  (std::log(zref / z0_arr(i,j,k)) - psi_h); // <w'Qv'>
2007  tflux += epsv * tm_arr(i,j,k) * qflux; // ~= <w'Tv'>
2008  if (w_star_arr) {
2009  // update w* and Umagmean
2010  w_star_arr(i,j,k) = calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,k));
2011  wstar = mdata.Bjr_beta * w_star_arr(i,j,k);
2012  umm = std::sqrt(umm_arr(i,j,k)*umm_arr(i,j,k) + wstar*wstar);
2013  umm = std::max(umm, WSMIN);
2014  }
2015  Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux);
2016  if ( (((Olen >= zero) && (Oleno <= zero)) ||
2017  ((Olen <= zero) && (Oleno >= zero))) &&
2018  std::fabs(Olen) + std::fabs(Oleno) < amrex::Real(1.0)) {
2019  Olen = myhalf * (Olen + Oleno);
2020  }
2021  Oleno = Olen;
2022  zeta = zref / Olen;
2023  psi_m = sfuns.calc_psi_m(zeta);
2024  psi_h = sfuns.calc_psi_h(zeta);
2025  u_star_arr(i,j,k) = mdata.kappa * umm / (std::log(zref / z0) - psi_m);
2026  ++iter;
2027  } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
2028  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(iter < max_iters,
2029  "Maximum number of MOST iterations reached.");
2030 
2031  // Populate the stored MOST arrays
2032  z0_arr(i,j,k) = z0;
2033  olen_arr(i,j,k) = Olen;
2034  t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) /
2035  (std::log(zref / z0) - psi_h);
2036  if (spec_qflux) {
2037  q_surf_arr(i,j,k) = mdata.surf_moist_flux * (std::log(zref / z0_arr(i,j,k)) - psi_h) /
2038  (u_star_arr(i,j,k) * mdata.kappa) + qvm_arr(i,j,k);
2039  q_star_arr(i,j,k) = -mdata.surf_moist_flux / u_star_arr(i,j,k);
2040  } else {
2041  q_star_arr(i,j,k) = mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) /
2042  (std::log(zref / z0_arr(i,j,k)) - psi_h);
2043  }
2044  }
2045 
2046 private:
2050  const amrex::Real tol = amrex::Real(1.0e-5);
2051 #ifdef AMREX_USE_FLOAT
2052  const amrex::Real eps = amrex::Real(1e-8);
2053 #else
2054  const amrex::Real eps = amrex::Real(1e-15);
2055 #endif
2056  const amrex::Real z0_eps = amrex::Real(1.0e-6);
2058  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
2059 };
2060 
2061 
2062 /**
2063  * Moeng flux formulation
2064  */
2066 {
2067  /**
2068  * Construct the Moeng flux calculator.
2069  */
2071 
2072  /**
2073  * Compute moisture flux for one surface point.
2074  *
2075  * @par Calling sequence
2076  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, qvm_arr,
2077  * u_star_arr, q_star_arr, q_surf_arr.
2078  */
2079  AMREX_GPU_DEVICE
2080  AMREX_FORCE_INLINE
2081  amrex::Real
2082  compute_q_flux (const int& i,
2083  const int& j,
2084  const int& k,
2085  const amrex::Array4<const amrex::Real>& cons_arr,
2086  const amrex::Array4<const amrex::Real>& velx_arr,
2087  const amrex::Array4<const amrex::Real>& vely_arr,
2088  const amrex::Array4<const amrex::Real>& umm_arr,
2089  const amrex::Array4<const amrex::Real>& qvm_arr,
2090  const amrex::Array4<const amrex::Real>& u_star_arr,
2091  const amrex::Array4<const amrex::Real>& q_star_arr,
2092  const amrex::Array4<const amrex::Real>& q_surf_arr) const
2093  {
2094  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2095  amrex::Real qv = cons_arr(i,j,k,RhoQ1_comp) / rho;
2096  amrex::Real velx = myhalf * ( velx_arr(i,j,k) + velx_arr(i+1,j ,k) );
2097  amrex::Real vely = myhalf * ( vely_arr(i,j,k) + vely_arr(i ,j+1,k) );
2098 
2099  amrex::Real qv_mean = qvm_arr(i,j,0);
2100  amrex::Real ustar = u_star_arr(i,j,0);
2101  amrex::Real qstar = q_star_arr(i,j,0);
2102  amrex::Real qv_surf = q_surf_arr(i,j,0);
2103  amrex::Real wsp_mean = umm_arr(i,j,0);
2104  wsp_mean = std::max(wsp_mean, WSMIN);
2105 
2106  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2107  amrex::Real num1 = wsp * (qv_mean-qv_surf);
2108  amrex::Real num2 = wsp_mean * (qv-qv_mean);
2109 
2110  // NOTE: this is rho*<Qv'w'> = -K dQvdz
2111  amrex::Real moflux = (std::abs(qstar) > eps) ?
2112  -rho*qstar*ustar*(num1+num2)/((qv_mean-qv_surf)*wsp_mean) : zero;
2113 
2114  return moflux;
2115  }
2116 
2117  /**
2118  * Compute heat flux for one surface point.
2119  *
2120  * @par Calling sequence
2121  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, tm_arr,
2122  * u_star_arr, t_star_arr, t_surf_arr.
2123  */
2124  AMREX_GPU_DEVICE
2125  AMREX_FORCE_INLINE
2126  amrex::Real
2127  compute_t_flux (const int& i,
2128  const int& j,
2129  const int& k,
2130  const amrex::Array4<const amrex::Real>& cons_arr,
2131  const amrex::Array4<const amrex::Real>& velx_arr,
2132  const amrex::Array4<const amrex::Real>& vely_arr,
2133  const amrex::Array4<const amrex::Real>& umm_arr,
2134  const amrex::Array4<const amrex::Real>& tm_arr,
2135  const amrex::Array4<const amrex::Real>& u_star_arr,
2136  const amrex::Array4<const amrex::Real>& t_star_arr,
2137  const amrex::Array4<const amrex::Real>& t_surf_arr) const
2138  {
2139  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2140  amrex::Real theta = cons_arr(i,j,k,RhoTheta_comp) / rho;
2141  amrex::Real velx = myhalf * ( velx_arr(i,j,k) + velx_arr(i+1,j ,k) );
2142  amrex::Real vely = myhalf * ( vely_arr(i,j,k) + vely_arr(i ,j+1,k) );
2143 
2144  amrex::Real theta_mean = tm_arr(i,j,0);
2145  amrex::Real ustar = u_star_arr(i,j,0);
2146  amrex::Real tstar = t_star_arr(i,j,0);
2147  amrex::Real theta_surf = t_surf_arr(i,j,0);
2148  amrex::Real wsp_mean = umm_arr(i,j,0);
2149  wsp_mean = std::max(wsp_mean, WSMIN);
2150 
2151  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2152  amrex::Real num1 = wsp * (theta_mean-theta_surf);
2153  amrex::Real num2 = wsp_mean * (theta-theta_mean);
2154 
2155  // NOTE: this is rho*<T'w'> = -K dTdz
2156  amrex::Real moflux = (std::abs(tstar) > eps) ?
2157  -rho*tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : zero;
2158 
2159  return moflux;
2160  }
2161 
2162  /**
2163  * Compute x-momentum flux for one surface face.
2164  *
2165  * @par Calling sequence
2166  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, um_arr,
2167  * u_star_arr.
2168  */
2169  AMREX_GPU_DEVICE
2170  AMREX_FORCE_INLINE
2171  amrex::Real
2172  compute_u_flux (const int& i,
2173  const int& j,
2174  const int& k,
2175  const amrex::Array4<const amrex::Real>& cons_arr,
2176  const amrex::Array4<const amrex::Real>& velx_arr,
2177  const amrex::Array4<const amrex::Real>& vely_arr,
2178  const amrex::Array4<const amrex::Real>& umm_arr,
2179  const amrex::Array4<const amrex::Real>& um_arr,
2180  const amrex::Array4<const amrex::Real>& u_star_arr) const
2181  {
2182  amrex::Real velx = velx_arr(i,j,k);
2183  amrex::Real vely = fourth * ( vely_arr(i ,j,k) + vely_arr(i ,j+1,k)
2184  + vely_arr(i-1,j,k) + vely_arr(i-1,j+1,k) );
2185  amrex::Real rho = myhalf * ( cons_arr(i-1,j,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2186 
2187  amrex::Real umean = um_arr(i,j,0);
2188  amrex::Real ustar = myhalf * ( u_star_arr(i-1,j,0) + u_star_arr(i,j,0) );
2189  amrex::Real wsp_mean = myhalf * ( umm_arr(i-1,j,0) + umm_arr(i,j,0) );
2190  wsp_mean = std::max(wsp_mean, WSMIN);
2191 
2192  // Note: The surface mean shear stress is decomposed into tau_xz by
2193  // multiplying the modeled shear stress (rho*ustar^2) with
2194  // a factor of umean/wsp_mean for directionality; this factor
2195  // modifies the denominator from what is in Moeng amrex::Real(1984.)
2196  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2197  amrex::Real num1 = wsp * umean;
2198  amrex::Real num2 = wsp_mean * (velx-umean);
2199 
2200  // NOTE: this is rho*<u'w'> = -K dudz
2201  amrex::Real stressx = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
2202 
2203  return stressx;
2204  }
2205 
2206  /**
2207  * Compute y-momentum flux for one surface face.
2208  *
2209  * @par Calling sequence
2210  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr,
2211  * u_star_arr.
2212  */
2213  AMREX_GPU_DEVICE
2214  AMREX_FORCE_INLINE
2215  amrex::Real
2216  compute_v_flux (const int& i,
2217  const int& j,
2218  const int& k,
2219  const amrex::Array4<const amrex::Real>& cons_arr,
2220  const amrex::Array4<const amrex::Real>& velx_arr,
2221  const amrex::Array4<const amrex::Real>& vely_arr,
2222  const amrex::Array4<const amrex::Real>& umm_arr,
2223  const amrex::Array4<const amrex::Real>& vm_arr,
2224  const amrex::Array4<const amrex::Real>& u_star_arr) const
2225  {
2226  amrex::Real velx = fourth * ( velx_arr(i,j ,k) + velx_arr(i+1,j ,k)
2227  + velx_arr(i,j-1,k) + velx_arr(i+1,j-1,k) );
2228  amrex::Real vely = vely_arr(i,j,k);
2229  amrex::Real rho = myhalf * ( cons_arr(i,j-1,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2230 
2231  amrex::Real vmean = vm_arr(i,j,0);
2232  amrex::Real ustar = myhalf * ( u_star_arr(i,j-1,0) + u_star_arr(i,j,0) );
2233  amrex::Real wsp_mean = myhalf * ( umm_arr(i,j-1,0) + umm_arr(i,j,0) );
2234  wsp_mean = std::max(wsp_mean, WSMIN);
2235 
2236  // Note: The surface mean shear stress is decomposed into tau_yz by
2237  // multiplying the modeled shear stress (rho*ustar^2) with
2238  // a factor of vmean/wsp_mean for directionality; this factor
2239  // modifies the denominator from what is in Moeng amrex::Real(1984.)
2240  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2241  amrex::Real num1 = wsp * vmean;
2242  amrex::Real num2 = wsp_mean * (vely-vmean);
2243 
2244  // NOTE: this is rho*<v'w'> = -K dvdz
2245  amrex::Real stressy = -rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
2246 
2247  return stressy;
2248  }
2249 
2250 private:
2251 #ifdef AMREX_USE_FLOAT
2252  const amrex::Real eps = amrex::Real(1e-8);
2253 #else
2254  const amrex::Real eps = amrex::Real(1e-15);
2255 #endif
2256  const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed
2257 };
2258 
2259 
2260 /**
2261  * Donelan flux formulation
2262  */
2264 {
2265  /**
2266  * Construct the Donelan flux calculator.
2267  */
2269 
2270  /**
2271  * Compute moisture flux for one surface point.
2272  *
2273  * @par Calling sequence
2274  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, qvm_arr,
2275  * u_star_arr, q_star_arr, q_surf_arr.
2276  */
2277  AMREX_GPU_DEVICE
2278  AMREX_FORCE_INLINE
2279  amrex::Real
2280  compute_q_flux (const int& /*i*/,
2281  const int& /*j*/,
2282  const int& /*k*/,
2283  const amrex::Array4<const amrex::Real>& /*cons_arr*/,
2284  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2285  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2286  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2287  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
2288  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2289  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
2290  const amrex::Array4<const amrex::Real>& /*q_surf_arr*/) const
2291  {
2292  // NOTE: this is rho*<q'w'> = -K dqdz
2293  amrex::Real moflux = zero;
2294 
2295  return moflux;
2296  }
2297 
2298  /**
2299  * Compute heat flux for one surface point.
2300  *
2301  * @par Calling sequence
2302  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, tm_arr,
2303  * u_star_arr, t_star_arr, t_surf_arr.
2304  */
2305  AMREX_GPU_DEVICE
2306  AMREX_FORCE_INLINE
2307  amrex::Real
2308  compute_t_flux (const int& i,
2309  const int& j,
2310  const int& k,
2311  const amrex::Array4<const amrex::Real>& cons_arr,
2312  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2313  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2314  const amrex::Array4<const amrex::Real>& umm_arr,
2315  const amrex::Array4<const amrex::Real>& tm_arr,
2316  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2317  const amrex::Array4<const amrex::Real>& /*t_star_arr*/,
2318  const amrex::Array4<const amrex::Real>& t_surf_arr) const
2319  {
2320  amrex::Real Ch = amrex::Real(0.0012);
2321  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2322 
2323  amrex::Real theta_surf = t_surf_arr(i,j,0);
2324  amrex::Real theta_mean = tm_arr(i,j,0);
2325  amrex::Real wsp_mean = umm_arr(i,j,0);
2326 
2327  // NOTE: this is rho*<T'w'> = -K dTdz
2328  amrex::Real moflux = -rho * Ch * wsp_mean * (theta_mean - theta_surf);
2329 
2330  return moflux;
2331  }
2332 
2333  /**
2334  * Compute x-momentum flux for one surface face.
2335  *
2336  * @par Calling sequence
2337  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, um_arr,
2338  * u_star_arr.
2339  */
2340  AMREX_GPU_DEVICE
2341  AMREX_FORCE_INLINE
2342  amrex::Real
2343  compute_u_flux (const int& i,
2344  const int& j,
2345  const int& k,
2346  const amrex::Array4<const amrex::Real>& cons_arr,
2347  const amrex::Array4<const amrex::Real>& velx_arr,
2348  const amrex::Array4<const amrex::Real>& vely_arr,
2349  const amrex::Array4<const amrex::Real>& umm_arr,
2350  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2351  const amrex::Array4<const amrex::Real>& /*u_star_arr*/) const
2352  {
2353  amrex::Real velx = velx_arr(i,j,k);
2354  amrex::Real vely = fourth * ( vely_arr(i ,j,k) + vely_arr(i ,j+1,k)
2355  + vely_arr(i-1,j,k) + vely_arr(i-1,j+1,k) );
2356  amrex::Real rho = myhalf * ( cons_arr(i-1,j,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2357 
2358  amrex::Real Cd = amrex::Real(0.001);
2359  const amrex::Real c = amrex::Real(7e-5);
2360  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2361  amrex::Real wsp_mean = myhalf * ( umm_arr(i-1,j,0) + umm_arr(i,j,0) );
2362  if (wsp_mean <= amrex::Real(5.0)) {
2363  Cd = amrex::Real(0.001);
2364  } else if (wsp_mean < amrex::Real(25.0) && wsp_mean > amrex::Real(5.0)) {
2365  Cd = amrex::Real(0.001) + c * (wsp_mean - amrex::Real(5.0));
2366  } else {
2367  Cd = amrex::Real(0.0024);
2368  }
2369 
2370  // NOTE: this is rho*<u'w'> = -K dudz
2371  amrex::Real stressx = -rho * Cd * velx * wsp;
2372 
2373  return stressx;
2374  }
2375 
2376  /**
2377  * Compute y-momentum flux for one surface face.
2378  *
2379  * @par Calling sequence
2380  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr,
2381  * u_star_arr.
2382  */
2383  AMREX_GPU_DEVICE
2384  AMREX_FORCE_INLINE
2385  amrex::Real
2386  compute_v_flux (const int& i,
2387  const int& j,
2388  const int& k,
2389  const amrex::Array4<const amrex::Real>& cons_arr,
2390  const amrex::Array4<const amrex::Real>& velx_arr,
2391  const amrex::Array4<const amrex::Real>& vely_arr,
2392  const amrex::Array4<const amrex::Real>& umm_arr,
2393  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2394  const amrex::Array4<const amrex::Real>& /*u_star_arr*/) const
2395  {
2396  amrex::Real velx = fourth * ( velx_arr(i,j ,k) + velx_arr(i+1,j ,k)
2397  + velx_arr(i,j-1,k) + velx_arr(i+1,j-1,k) );
2398  amrex::Real vely = vely_arr(i,j,k);
2399  amrex::Real rho = myhalf * ( cons_arr(i,j-1,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2400 
2401  amrex::Real Cd = amrex::Real(0.001);
2402  const amrex::Real c = amrex::Real(7e-5);
2403  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2404  amrex::Real wsp_mean = myhalf * ( umm_arr(i,j-1,0) + umm_arr(i,j,0) );
2405  if (wsp_mean <= amrex::Real(5.0)) {
2406  Cd = amrex::Real(0.001);
2407  } else if (wsp_mean < amrex::Real(25.0) && wsp_mean > amrex::Real(5.0)) {
2408  Cd = amrex::Real(0.001) + c * (wsp_mean - amrex::Real(5.0));
2409  } else {
2410  Cd = amrex::Real(0.0024);
2411  }
2412 
2413  // NOTE: this is rho*<v'w'> = -K dvdz
2414  amrex::Real stressy = -rho * Cd * vely * wsp;
2415 
2416  return stressy;
2417  }
2418 
2419 private:
2420 
2421 };
2422 
2423 
2424 /**
2425  * Custom flux formulation
2426  */
2428 {
2429  /**
2430  * Construct the custom flux calculator.
2431  *
2432  * @param[in] specified_rho_surf whether the supplied custom fluxes include density
2433  */
2434  custom_flux (bool specified_rho_surf)
2435  : fluxes_include_rho(specified_rho_surf)
2436  {}
2437 
2438  /**
2439  * Compute moisture flux for one surface point.
2440  *
2441  * @par Calling sequence
2442  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, qvm_arr,
2443  * u_star_arr, q_star_arr, q_surf_arr.
2444  */
2445  AMREX_GPU_DEVICE
2446  AMREX_FORCE_INLINE
2447  amrex::Real
2448  compute_q_flux (const int& i,
2449  const int& j,
2450  const int& k,
2451  const amrex::Array4<const amrex::Real>& cons_arr,
2452  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2453  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2454  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2455  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
2456  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2457  const amrex::Array4<const amrex::Real>& q_star_arr,
2458  const amrex::Array4<const amrex::Real>& /*q_surf_arr*/) const
2459  {
2460  amrex::Real rho = (fluxes_include_rho) ? one : cons_arr(i,j,k,Rho_comp);
2461  amrex::Real qstar = q_star_arr(i,j,0);
2462 
2463  // NOTE: this is rho*<q'w'> = -K dqdz
2464  amrex::Real moflux = (std::abs(qstar) > eps) ? rho * qstar : zero;
2465 
2466  return moflux;
2467  }
2468 
2469  /**
2470  * Compute heat flux for one surface point.
2471  *
2472  * @par Calling sequence
2473  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, tm_arr,
2474  * u_star_arr, t_star_arr, t_surf_arr.
2475  */
2476  AMREX_GPU_DEVICE
2477  AMREX_FORCE_INLINE
2478  amrex::Real
2479  compute_t_flux (const int& i,
2480  const int& j,
2481  const int& k,
2482  const amrex::Array4<const amrex::Real>& cons_arr,
2483  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2484  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2485  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2486  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
2487  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2488  const amrex::Array4<const amrex::Real>& t_star_arr,
2489  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/) const
2490  {
2491  amrex::Real rho = (fluxes_include_rho) ? one : cons_arr(i,j,k,Rho_comp);
2492  amrex::Real tstar = t_star_arr(i,j,0);
2493 
2494  // NOTE: this is rho*<T'w'> = -K dTdz
2495  amrex::Real moflux = (std::abs(tstar) > eps) ? rho * tstar : zero;
2496 
2497  return moflux;
2498  }
2499 
2500  /**
2501  * Compute x-momentum flux for one surface face.
2502  *
2503  * @par Calling sequence
2504  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, um_arr,
2505  * u_star_arr.
2506  */
2507  AMREX_GPU_DEVICE
2508  AMREX_FORCE_INLINE
2509  amrex::Real
2510  compute_u_flux (const int& i,
2511  const int& j,
2512  const int& k,
2513  const amrex::Array4<const amrex::Real>& cons_arr,
2514  const amrex::Array4<const amrex::Real>& velx_arr,
2515  const amrex::Array4<const amrex::Real>& vely_arr,
2516  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2517  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2518  const amrex::Array4<const amrex::Real>& u_star_arr) const
2519  {
2520  amrex::Real velx = velx_arr(i,j,k);
2521  amrex::Real vely = fourth * ( vely_arr(i ,j,k) + vely_arr(i ,j+1,k)
2522  + vely_arr(i-1,j,k) + vely_arr(i-1,j+1,k) );
2523  amrex::Real rho = (fluxes_include_rho) ? one : myhalf * ( cons_arr(i-1,j,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2524 
2525  amrex::Real ustar = myhalf * ( u_star_arr(i-1,j,0) + u_star_arr(i,j,0) );
2526  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2527 
2528  // NOTE: this is rho*<u'w'> = -K dudz
2529  amrex::Real stressx = -rho * ustar * ustar * velx / wsp;
2530 
2531  return stressx;
2532  }
2533 
2534  /**
2535  * Compute y-momentum flux for one surface face.
2536  *
2537  * @par Calling sequence
2538  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr,
2539  * u_star_arr.
2540  */
2541  AMREX_GPU_DEVICE
2542  AMREX_FORCE_INLINE
2543  amrex::Real
2544  compute_v_flux (const int& i,
2545  const int& j,
2546  const int& k,
2547  const amrex::Array4<const amrex::Real>& cons_arr,
2548  const amrex::Array4<const amrex::Real>& velx_arr,
2549  const amrex::Array4<const amrex::Real>& vely_arr,
2550  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2551  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2552  const amrex::Array4<const amrex::Real>& u_star_arr) const
2553  {
2554  amrex::Real velx = fourth * ( velx_arr(i,j ,k) + velx_arr(i+1,j ,k)
2555  + velx_arr(i,j-1,k) + velx_arr(i+1,j-1,k) );
2556  amrex::Real vely = vely_arr(i,j,k);
2557  amrex::Real rho = (fluxes_include_rho) ? one : myhalf * ( cons_arr(i,j-1,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2558 
2559  amrex::Real ustar = myhalf * ( u_star_arr(i,j-1,0) + u_star_arr(i,j,0) );
2560  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2561 
2562  // NOTE: this is rho*<v'w'> = -K dvdz
2563  amrex::Real stressy = -rho * ustar * ustar * vely / wsp;
2564 
2565  return stressy;
2566  }
2567 
2568 private:
2569 #ifdef AMREX_USE_FLOAT
2570  const amrex::Real eps = amrex::Real(1e-8);
2571 #else
2572  const amrex::Real eps = amrex::Real(1e-15);
2573 #endif
2574  const bool fluxes_include_rho{false};
2575 };
2576 
2577 
2578 /**
2579  * Bulk coefficient flux formulation
2580  */
2582 {
2583  /**
2584  * Construct the bulk-coefficient flux calculator.
2585  *
2586  * @param[in] m_Cd drag coefficient
2587  * @param[in] m_Ch heat-transfer coefficient
2588  * @param[in] m_Cq moisture-transfer coefficient
2589  */
2591  amrex::Real m_Ch,
2592  amrex::Real m_Cq)
2593  {
2594  mdata.Cd = m_Cd;
2595  mdata.Ch = m_Ch;
2596  mdata.Cq = m_Cq;
2597  }
2598 
2599  /**
2600  * Compute moisture flux for one surface point.
2601  *
2602  * @par Calling sequence
2603  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, qvm_arr,
2604  * u_star_arr, q_star_arr, q_surf_arr.
2605  */
2606  AMREX_GPU_DEVICE
2607  AMREX_FORCE_INLINE
2608  amrex::Real
2609  compute_q_flux (const int& i,
2610  const int& j,
2611  const int& k,
2612  const amrex::Array4<const amrex::Real>& cons_arr,
2613  const amrex::Array4<const amrex::Real>& velx_arr,
2614  const amrex::Array4<const amrex::Real>& vely_arr,
2615  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2616  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
2617  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2618  const amrex::Array4<const amrex::Real>& /*q_star_arr*/,
2619  const amrex::Array4<const amrex::Real>& q_surf_arr) const
2620  {
2621  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2622  amrex::Real qv = cons_arr(i,j,k,RhoQ1_comp)/cons_arr(i,j,k,Rho_comp);
2623  amrex::Real qvsurf = q_surf_arr(i,j,0);
2624  amrex::Real velx = myhalf * ( velx_arr(i,j,k) + velx_arr(i+1,j ,k) );
2625  amrex::Real vely = myhalf * ( vely_arr(i,j,k) + vely_arr(i ,j+1,k) );
2626  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2627 
2628  // NOTE: this is rho*<q'w'> = -K dqdz
2629  amrex::Real moflux = rho * mdata.Cq * wsp * (qvsurf - qv);
2630 
2631  return moflux;
2632  }
2633 
2634  /**
2635  * Compute heat flux for one surface point.
2636  *
2637  * @par Calling sequence
2638  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, tm_arr,
2639  * u_star_arr, t_star_arr, t_surf_arr.
2640  */
2641  AMREX_GPU_DEVICE
2642  AMREX_FORCE_INLINE
2643  amrex::Real
2644  compute_t_flux (const int& i,
2645  const int& j,
2646  const int& k,
2647  const amrex::Array4<const amrex::Real>& cons_arr,
2648  const amrex::Array4<const amrex::Real>& velx_arr,
2649  const amrex::Array4<const amrex::Real>& vely_arr,
2650  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2651  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
2652  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2653  const amrex::Array4<const amrex::Real>& /*t_star_arr*/,
2654  const amrex::Array4<const amrex::Real>& t_surf_arr) const
2655  {
2656  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2657  amrex::Real th = cons_arr(i,j,k,RhoTheta_comp)/cons_arr(i,j,k,Rho_comp);
2658  amrex::Real thsurf = t_surf_arr(i,j,0);
2659  amrex::Real velx = myhalf * ( velx_arr(i,j,k) + velx_arr(i+1,j ,k) );
2660  amrex::Real vely = myhalf * ( vely_arr(i,j,k) + vely_arr(i ,j+1,k) );
2661  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2662 
2663  // NOTE: this is rho*<T'w'> = -K dTdz
2664  amrex::Real moflux = rho * mdata.Ch * wsp * (thsurf - th);
2665 
2666  return moflux;
2667  }
2668 
2669  /**
2670  * Compute x-momentum flux for one surface face.
2671  *
2672  * @par Calling sequence
2673  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, um_arr,
2674  * u_star_arr.
2675  */
2676  AMREX_GPU_DEVICE
2677  AMREX_FORCE_INLINE
2678  amrex::Real
2679  compute_u_flux (const int& i,
2680  const int& j,
2681  const int& k,
2682  const amrex::Array4<const amrex::Real>& cons_arr,
2683  const amrex::Array4<const amrex::Real>& velx_arr,
2684  const amrex::Array4<const amrex::Real>& vely_arr,
2685  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2686  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2687  const amrex::Array4<const amrex::Real>& /*u_star_arr*/) const
2688  {
2689  amrex::Real rho = myhalf * ( cons_arr(i-1,j,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2690  amrex::Real velx = velx_arr(i,j,k);
2691  amrex::Real vely = fourth * ( vely_arr(i ,j,k) + vely_arr(i ,j+1,k)
2692  + vely_arr(i-1,j,k) + vely_arr(i-1,j+1,k) );
2693  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2694 
2695  // NOTE: this is rho*<u'w'> = -K dudz
2696  // NOTE: tau_tot = rho * Cd * wsp^2, multiply by u/wsp to get tau_13
2697  amrex::Real stressx = -rho * mdata.Cd * wsp * velx;
2698 
2699  return stressx;
2700  }
2701 
2702  /**
2703  * Compute y-momentum flux for one surface face.
2704  *
2705  * @par Calling sequence
2706  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr,
2707  * u_star_arr.
2708  */
2709  AMREX_GPU_DEVICE
2710  AMREX_FORCE_INLINE
2711  amrex::Real
2712  compute_v_flux (const int& i,
2713  const int& j,
2714  const int& k,
2715  const amrex::Array4<const amrex::Real>& cons_arr,
2716  const amrex::Array4<const amrex::Real>& velx_arr,
2717  const amrex::Array4<const amrex::Real>& vely_arr,
2718  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2719  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2720  const amrex::Array4<const amrex::Real>& /*u_star_arr*/) const
2721  {
2722  amrex::Real rho = myhalf * ( cons_arr(i,j-1,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2723  amrex::Real velx = fourth * ( velx_arr(i,j ,k) + velx_arr(i+1,j ,k)
2724  + velx_arr(i,j-1,k) + velx_arr(i+1,j-1,k) );
2725  amrex::Real vely = vely_arr(i,j,k);
2726  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2727 
2728  // NOTE: this is rho*<v'w'> = -K dvdz
2729  // NOTE: tau_tot = rho * Cd * wsp^2, multiply by v/wsp to get tau_13
2730  amrex::Real stressy = -rho * mdata.Cd * wsp * vely;
2731 
2732  return stressy;
2733  }
2734 
2735 private:
2737 };
2738 
2739 
2740 /**
2741  * RICO flux formulation
2742  */
2744 {
2745  /**
2746  * Construct the RICO flux calculator.
2747  *
2748  * @param[in] l_theta_z0 surface reference potential temperature
2749  * @param[in] l_qsat_z0 surface saturation specific humidity
2750  */
2751  rico_flux (amrex::Real l_theta_z0, amrex::Real l_qsat_z0)
2752  : theta_z0{l_theta_z0}, qsat_z0{l_qsat_z0} {}
2753 
2754 
2755  /**
2756  * Compute moisture flux for one surface point.
2757  *
2758  * @par Calling sequence
2759  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, qvm_arr,
2760  * u_star_arr, q_star_arr, q_surf_arr.
2761  */
2762  AMREX_GPU_DEVICE
2763  AMREX_FORCE_INLINE
2764  amrex::Real
2765  compute_q_flux (const int& i,
2766  const int& j,
2767  const int& k,
2768  const amrex::Array4<const amrex::Real>& cons_arr,
2769  const amrex::Array4<const amrex::Real>& velx_arr,
2770  const amrex::Array4<const amrex::Real>& vely_arr,
2771  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2772  const amrex::Array4<const amrex::Real>& /*qm_arr*/,
2773  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2774  const amrex::Array4<const amrex::Real>& q_star_arr,
2775  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/) const
2776  {
2777  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2778  amrex::Real qstar = q_star_arr(i,j,0);
2779  amrex::Real q = cons_arr(i, j, k, RhoQ1_comp) / rho;
2780 
2781  amrex::Real velx = myhalf * (velx_arr(i,j,k) + velx_arr(i+1,j,k));
2782  amrex::Real vely = myhalf * (vely_arr(i,j,k) + vely_arr(i,j+1,k));
2783  amrex::Real wsp = std::sqrt(velx*velx + vely*vely);
2784 
2785  // NOTE: this is rho*<q'w'> = -K dqdz
2786  amrex::Real moflux = (std::abs(qstar) > eps) ? - rho * qstar * wsp * (q - qsat_z0) : zero;
2787 
2788  return moflux;
2789  }
2790 
2791  /**
2792  * Compute heat flux for one surface point.
2793  *
2794  * @par Calling sequence
2795  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, tm_arr,
2796  * u_star_arr, t_star_arr, t_surf_arr.
2797  */
2798  AMREX_GPU_DEVICE
2799  AMREX_FORCE_INLINE
2800  amrex::Real
2801  compute_t_flux (const int& i,
2802  const int& j,
2803  const int& k,
2804  const amrex::Array4<const amrex::Real>& cons_arr,
2805  const amrex::Array4<const amrex::Real>& velx_arr,
2806  const amrex::Array4<const amrex::Real>& vely_arr,
2807  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2808  const amrex::Array4<const amrex::Real>& /*tm_arr*/,
2809  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2810  const amrex::Array4<const amrex::Real>& t_star_arr,
2811  const amrex::Array4<const amrex::Real>& /*t_surf_arr*/) const
2812  {
2813  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2814  amrex::Real tstar = t_star_arr(i, j, 0);
2815  amrex::Real theta = cons_arr(i, j, k, RhoTheta_comp) / rho;
2816 
2817  amrex::Real velx = myhalf * (velx_arr(i,j,k) + velx_arr(i+1,j,k));
2818  amrex::Real vely = myhalf * (vely_arr(i,j,k) + vely_arr(i,j+1,k));
2819  amrex::Real wsp = std::sqrt(velx*velx + vely*vely);
2820 
2821  // NOTE: this is rho*<T'w'> = -K dTdz
2822  amrex::Real moflux = (std::abs(tstar) > eps) ? - rho * tstar * wsp * (theta - theta_z0) : zero;
2823 
2824  return moflux;
2825  }
2826 
2827  /**
2828  * Compute x-momentum flux for one surface face.
2829  *
2830  * @par Calling sequence
2831  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, um_arr,
2832  * u_star_arr.
2833  */
2834  AMREX_GPU_DEVICE
2835  AMREX_FORCE_INLINE
2836  amrex::Real
2837  compute_u_flux (const int& i,
2838  const int& j,
2839  const int& k,
2840  const amrex::Array4<const amrex::Real>& cons_arr,
2841  const amrex::Array4<const amrex::Real>& velx_arr,
2842  const amrex::Array4<const amrex::Real>& vely_arr,
2843  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2844  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2845  const amrex::Array4<const amrex::Real>& u_star_arr) const
2846  {
2847  amrex::Real velx = velx_arr(i,j,k);
2848  amrex::Real vely = fourth * ( vely_arr(i ,j,k) + vely_arr(i ,j+1,k)
2849  + vely_arr(i-1,j,k) + vely_arr(i-1,j+1,k) );
2850  amrex::Real rho = myhalf * ( cons_arr(i-1,j,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2851 
2852  amrex::Real ustar = myhalf * ( u_star_arr(i-1,j,0) + u_star_arr(i,j,0) );
2853  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2854 
2855  // NOTE: this is rho*<u'w'> = -K dudz
2856  amrex::Real stressx = -rho * ustar * wsp * velx;
2857 
2858  return stressx;
2859  }
2860 
2861  /**
2862  * Compute y-momentum flux for one surface face.
2863  *
2864  * @par Calling sequence
2865  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr,
2866  * u_star_arr.
2867  */
2868  AMREX_GPU_DEVICE
2869  AMREX_FORCE_INLINE
2870  amrex::Real
2871  compute_v_flux (const int& i,
2872  const int& j,
2873  const int& k,
2874  const amrex::Array4<const amrex::Real>& cons_arr,
2875  const amrex::Array4<const amrex::Real>& velx_arr,
2876  const amrex::Array4<const amrex::Real>& vely_arr,
2877  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2878  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
2879  const amrex::Array4<const amrex::Real>& u_star_arr) const
2880  {
2881  amrex::Real velx = fourth * ( velx_arr(i,j ,k) + velx_arr(i+1,j ,k)
2882  + velx_arr(i,j-1,k) + velx_arr(i+1,j-1,k) );
2883  amrex::Real vely = vely_arr(i,j,k);
2884  amrex::Real rho = myhalf * ( cons_arr(i,j-1,k,Rho_comp) + cons_arr(i,j,k,Rho_comp) );
2885 
2886  amrex::Real ustar = myhalf * ( u_star_arr(i,j-1,0) + u_star_arr(i,j,0) );
2887  amrex::Real wsp = std::sqrt(velx*velx+vely*vely);
2888 
2889  // NOTE: this is rho*<v'w'> = -K dvdz
2890  amrex::Real stressy = -rho * ustar * wsp * vely;
2891 
2892  return stressy;
2893  }
2894 
2895 private:
2896 #ifdef AMREX_USE_FLOAT
2897  const amrex::Real eps = amrex::Real(1e-8);
2898 #else
2899  const amrex::Real eps = amrex::Real(1e-15);
2900 #endif
2903 };
2904 
2905 
2906 /**
2907  * Rotate flux formulation
2908  */
2910 {
2911  /**
2912  * Construct the rotated-flux calculator.
2913  */
2915 
2916  /**
2917  * Compute moisture flux for one surface point.
2918  *
2919  * @par Calling sequence
2920  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, qvm_arr,
2921  * u_star_arr, q_star_arr, q_surf_arr.
2922  */
2923  AMREX_GPU_DEVICE
2924  AMREX_FORCE_INLINE
2925  amrex::Real
2926  compute_q_flux (const int& i,
2927  const int& j,
2928  const int& k,
2929  const amrex::Array4<const amrex::Real>& cons_arr,
2930  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2931  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2932  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2933  const amrex::Array4<const amrex::Real>& qvm_arr,
2934  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2935  const amrex::Array4<const amrex::Real>& q_star_arr,
2936  const amrex::Array4<const amrex::Real>& q_surf_arr) const
2937  {
2938  // NOTE: this is the total stress
2939  amrex::Real qv_mean = qvm_arr(i,j,0);
2940  amrex::Real qv_surf = q_surf_arr(i,j,0);
2941  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2942  amrex::Real qstar = q_star_arr(i,j,0);
2943  amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar * (qv_mean-qv_surf): zero;
2944 
2945  return moflux;
2946  }
2947 
2948  /**
2949  * Compute heat flux for one surface point.
2950  *
2951  * @par Calling sequence
2952  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, tm_arr,
2953  * u_star_arr, t_star_arr, t_surf_arr.
2954  */
2955  AMREX_GPU_DEVICE
2956  AMREX_FORCE_INLINE
2957  amrex::Real
2958  compute_t_flux (const int& i,
2959  const int& j,
2960  const int& k,
2961  const amrex::Array4<const amrex::Real>& cons_arr,
2962  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2963  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2964  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2965  const amrex::Array4<const amrex::Real>& tm_arr,
2966  const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
2967  const amrex::Array4<const amrex::Real>& t_star_arr,
2968  const amrex::Array4<const amrex::Real>& t_surf_arr) const
2969  {
2970  // NOTE: this is the total stress
2971  amrex::Real theta_mean = tm_arr(i,j,0);
2972  amrex::Real theta_surf = t_surf_arr(i,j,0);
2973  amrex::Real rho = cons_arr(i,j,k,Rho_comp);
2974  amrex::Real tstar = t_star_arr(i,j,0);
2975  amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar * (theta_mean-theta_surf) : zero;
2976 
2977  return moflux;
2978  }
2979 
2980  /**
2981  * Compute x-momentum flux for one surface face.
2982  *
2983  * @par Calling sequence
2984  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, um_arr,
2985  * u_star_arr.
2986  */
2987  AMREX_GPU_DEVICE
2988  AMREX_FORCE_INLINE
2989  amrex::Real
2990  compute_u_flux (const int& i,
2991  const int& j,
2992  const int& k,
2993  const amrex::Array4<const amrex::Real>& cons_arr,
2994  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
2995  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
2996  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
2997  const amrex::Array4<const amrex::Real>& /*um_arr*/,
2998  const amrex::Array4<const amrex::Real>& u_star_arr) const
2999  {
3000  // NOTE: this is the total stress
3001  amrex::Real rho = myhalf *( cons_arr(i-1,j,k,Rho_comp) + cons_arr(i ,j,k,Rho_comp) );
3002  amrex::Real ustar = myhalf * ( u_star_arr(i-1,j,0) + u_star_arr(i,j,0) );
3003  amrex::Real stressx = -rho * ustar * ustar;
3004 
3005  return stressx;
3006  }
3007 
3008  /**
3009  * Compute y-momentum flux for one surface face.
3010  *
3011  * @par Calling sequence
3012  * i, j, k, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr,
3013  * u_star_arr.
3014  */
3015  AMREX_GPU_DEVICE
3016  AMREX_FORCE_INLINE
3017  amrex::Real
3018  compute_v_flux (const int& /*i*/,
3019  const int& /*j*/,
3020  const int& /*k*/,
3021  const amrex::Array4<const amrex::Real>& /*cons_arr*/,
3022  const amrex::Array4<const amrex::Real>& /*velx_arr*/,
3023  const amrex::Array4<const amrex::Real>& /*vely_arr*/,
3024  const amrex::Array4<const amrex::Real>& /*umm_arr*/,
3025  const amrex::Array4<const amrex::Real>& /*vm_arr*/,
3026  const amrex::Array4<const amrex::Real>& /*u_star_arr*/) const
3027  {
3028  // NOTE: this is the total stress
3029  amrex::Real stressy = zero;
3030 
3031  return stressy;
3032  }
3033 
3034 private:
3035 #ifdef AMREX_USE_FLOAT
3036  const amrex::Real eps = amrex::Real(1e-8);
3037 #else
3038  const amrex::Real eps = amrex::Real(1e-15);
3039 #endif
3040 };
3041 #endif
constexpr amrex::Real epsv
Definition: ERF_Constants.H:53
constexpr amrex::Real three
Definition: ERF_Constants.H:11
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:63
constexpr amrex::Real two
Definition: ERF_Constants.H:10
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:43
@ num
Definition: ERF_DataStruct.H:27
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
rho
Definition: ERF_InitCustomPert_Bubble.H:107
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Donelan_roughness(amrex::Real ustar)
Definition: ERF_MOSTRoughness.H:35
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real COARE3_roughness(amrex::Real zref, amrex::Real umm, amrex::Real ustar)
Definition: ERF_MOSTRoughness.H:13
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real air_viscosity(amrex::Real T_degK)
Definition: ERF_MOSTStress.H:151
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real calc_wstar(const amrex::Real &ust, const amrex::Real &tst, const amrex::Real &qst, const amrex::Real &pblh, const amrex::Real &th, const amrex::Real &thv, const amrex::Real &qv=amrex::Real(0))
Definition: ERF_Wstar.H:13
@ Mom_v
Definition: ERF_IndexDefines.H:211
@ theta
Definition: ERF_SLM.H:20
@ qv
Definition: ERF_Kessler.H:30
@ q
Definition: ERF_WSM6.H:184
@ den
Definition: ERF_AdvanceWSM6.cpp:109
Definition: ERF_MOSTStress.H:226
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:328
most_data mdata
Definition: ERF_MOSTStress.H:321
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:258
const amrex::Real tol_z
Definition: ERF_MOSTStress.H:326
similarity_funs sfuns
Definition: ERF_MOSTStress.H:322
adiabatic_charnock(amrex::Real Tflux, amrex::Real Qvflux, amrex::Real cnk_a, bool cnk_visc)
Definition: ERF_MOSTStress.H:235
Definition: ERF_MOSTStress.H:437
similarity_funs sfuns
Definition: ERF_MOSTStress.H:512
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:514
most_data mdata
Definition: ERF_MOSTStress.H:511
adiabatic_donelan(amrex::Real Tflux, amrex::Real Qvflux)
Definition: ERF_MOSTStress.H:444
const amrex::Real tol
Definition: ERF_MOSTStress.H:513
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:463
Definition: ERF_MOSTStress.H:336
similarity_funs sfuns
Definition: ERF_MOSTStress.H:423
const amrex::Real tol_z
Definition: ERF_MOSTStress.H:427
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:429
adiabatic_mod_charnock(amrex::Real Tflux, amrex::Real Qvflux, amrex::Real depth)
Definition: ERF_MOSTStress.H:344
most_data mdata
Definition: ERF_MOSTStress.H:422
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:366
Definition: ERF_MOSTStress.H:522
const amrex::Real eps
Definition: ERF_MOSTStress.H:605
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:548
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:607
most_data mdata
Definition: ERF_MOSTStress.H:599
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:609
similarity_funs sfuns
Definition: ERF_MOSTStress.H:600
const amrex::Real tol
Definition: ERF_MOSTStress.H:601
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:608
adiabatic_wave_coupled(amrex::Real Tflux, amrex::Real Qvflux)
Definition: ERF_MOSTStress.H:529
Definition: ERF_MOSTStress.H:162
adiabatic(amrex::Real Tflux, amrex::Real Qvflux)
Definition: ERF_MOSTStress.H:169
similarity_funs sfuns
Definition: ERF_MOSTStress.H:218
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< const amrex::Real > &z0_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &t_star_arr, const amrex::Array4< amrex::Real > &q_star_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &) const
Definition: ERF_MOSTStress.H:188
most_data mdata
Definition: ERF_MOSTStress.H:217
Definition: ERF_MOSTStress.H:2582
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 > &) const
Definition: ERF_MOSTStress.H:2679
bulk_coeff_flux(amrex::Real m_Cd, amrex::Real m_Ch, amrex::Real m_Cq)
Definition: ERF_MOSTStress.H:2590
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 > &, 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_surf_arr) const
Definition: ERF_MOSTStress.H:2609
most_data mdata
Definition: ERF_MOSTStress.H:2736
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 > &, 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_surf_arr) const
Definition: ERF_MOSTStress.H:2644
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 > &) const
Definition: ERF_MOSTStress.H:2712
Definition: ERF_MOSTStress.H:2428
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:2510
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:2448
custom_flux(bool specified_rho_surf)
Definition: ERF_MOSTStress.H:2434
const bool fluxes_include_rho
Definition: ERF_MOSTStress.H:2574
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:2479
const amrex::Real eps
Definition: ERF_MOSTStress.H:2572
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:2544
Definition: ERF_MOSTStress.H:2264
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:2343
donelan_flux()
Definition: ERF_MOSTStress.H:2268
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:2386
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:2308
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:2280
Definition: ERF_MOSTStress.H:2066
moeng_flux()
Definition: ERF_MOSTStress.H:2070
const amrex::Real eps
Definition: ERF_MOSTStress.H:2254
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:2172
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:2256
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:2127
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:2082
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:2216
Definition: ERF_MOSTStress.H:13
amrex::Real surf_moist_flux
Moisture flux.
Definition: ERF_MOSTStress.H:19
amrex::Real Cnk_b2
Modified Charnock Eq (4) https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137....
Definition: ERF_MOSTStress.H:23
amrex::Real Cnk_b
Definition: ERF_MOSTStress.H:25
amrex::Real Ch
Definition: ERF_MOSTStress.H:29
amrex::Real Cnk_d
Modified Charnock Eq (4) https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137....
Definition: ERF_MOSTStress.H:24
amrex::Real kappa
von Karman constant
Definition: ERF_MOSTStress.H:16
amrex::Real gravity
Acceleration due to gravity (m/s^2)
Definition: ERF_MOSTStress.H:17
amrex::Real Cnk_a
Standard Charnock constant https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137....
Definition: ERF_MOSTStress.H:21
amrex::Real Cd
Definition: ERF_MOSTStress.H:28
amrex::Real Cq
Definition: ERF_MOSTStress.H:30
const amrex::Real Bjr_beta
Definition: ERF_MOSTStress.H:32
amrex::Real Cnk_b1
Modified Charnock Eq (4) https://doi.org/amrex::Real(10.1175)/JAMC-D-17-amrex::Real(0137....
Definition: ERF_MOSTStress.H:22
amrex::Real z0_const
Roughness height – default constant value(m)
Definition: ERF_MOSTStress.H:15
bool visc
Use viscous Charnock formulation.
Definition: ERF_MOSTStress.H:26
amrex::Real surf_temp_flux
Heat flux TODO: decide whether this is <θ'w'> or <θv'w'> under moist conditions.
Definition: ERF_MOSTStress.H:18
Definition: ERF_MOSTStress.H:2744
rico_flux(amrex::Real l_theta_z0, amrex::Real l_qsat_z0)
Definition: ERF_MOSTStress.H:2751
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 > &, 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:2765
amrex::Real qsat_z0
Definition: ERF_MOSTStress.H:2902
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 > &, 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:2801
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:2837
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:2871
const amrex::Real eps
Definition: ERF_MOSTStress.H:2899
amrex::Real theta_z0
Definition: ERF_MOSTStress.H:2901
Definition: ERF_MOSTStress.H:2910
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:3018
rotate_flux()
Definition: ERF_MOSTStress.H:2914
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:2990
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:2958
const amrex::Real eps
Definition: ERF_MOSTStress.H:3038
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:2926
Definition: ERF_MOSTStress.H:40
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:105
amrex::Real beta_m
Constants from Dyer, BLM, 1974.
Definition: ERF_MOSTStress.H:135
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h2(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:77
amrex::Real beta_h
https://doi.org/amrex::Real(10.1007)/BF00240838
Definition: ERF_MOSTStress.H:136
amrex::Real gamma_h
Definition: ERF_MOSTStress.H:138
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:124
amrex::Real gamma_m
Definition: ERF_MOSTStress.H:137
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m2(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:52
Definition: ERF_MOSTStress.H:738
const amrex::Real tol
Definition: ERF_MOSTStress.H:866
surface_flux_charnock(amrex::Real Tflux, amrex::Real Qvflux, amrex::Real cnk_a, bool cnk_visc, bool cons_qflux)
Definition: ERF_MOSTStress.H:748
similarity_funs sfuns
Definition: ERF_MOSTStress.H:865
most_data mdata
Definition: ERF_MOSTStress.H:863
bool spec_qflux
Definition: ERF_MOSTStress.H:864
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:773
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:867
Definition: ERF_MOSTStress.H:1003
bool spec_qflux
Definition: ERF_MOSTStress.H:1116
const amrex::Real tol
Definition: ERF_MOSTStress.H:1118
surface_flux_donelan(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Definition: ERF_MOSTStress.H:1011
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1119
most_data mdata
Definition: ERF_MOSTStress.H:1115
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1117
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:1032
Definition: ERF_MOSTStress.H:875
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:908
surface_flux_mod_charnock(amrex::Real Tflux, amrex::Real Qvflux, amrex::Real depth, bool cons_qflux)
Definition: ERF_MOSTStress.H:884
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:995
most_data mdata
Definition: ERF_MOSTStress.H:991
bool spec_qflux
Definition: ERF_MOSTStress.H:992
const amrex::Real tol
Definition: ERF_MOSTStress.H:994
similarity_funs sfuns
Definition: ERF_MOSTStress.H:993
Definition: ERF_MOSTStress.H:1127
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1256
surface_flux_wave_coupled(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Definition: ERF_MOSTStress.H:1135
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:1255
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:1254
most_data mdata
Definition: ERF_MOSTStress.H:1245
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1247
const amrex::Real tol
Definition: ERF_MOSTStress.H:1248
bool spec_qflux
Definition: ERF_MOSTStress.H:1246
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, 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:1156
const amrex::Real eps
Definition: ERF_MOSTStress.H:1252
Definition: ERF_MOSTStress.H:617
similarity_funs sfuns
Definition: ERF_MOSTStress.H:728
surface_flux(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Definition: ERF_MOSTStress.H:625
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< const amrex::Real > &z0_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &tvm_arr, const amrex::Array4< const amrex::Real > &qvm_arr, const amrex::Array4< amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &w_star_arr, const amrex::Array4< amrex::Real > &t_star_arr, const amrex::Array4< amrex::Real > &q_star_arr, const amrex::Array4< amrex::Real > &t_surf_arr, const amrex::Array4< amrex::Real > &q_surf_arr, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &pblh_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &) const
Definition: ERF_MOSTStress.H:646
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:730
most_data mdata
Definition: ERF_MOSTStress.H:726
const amrex::Real tol
Definition: ERF_MOSTStress.H:729
bool spec_qflux
Definition: ERF_MOSTStress.H:727
Definition: ERF_MOSTStress.H:1414
most_data mdata
Definition: ERF_MOSTStress.H:1586
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:1449
bool spec_qflux
Definition: ERF_MOSTStress.H:1587
surface_temp_charnock(amrex::Real Tflux, amrex::Real Qvflux, amrex::Real cnk_a, bool cnk_visc, bool cons_qflux)
Definition: ERF_MOSTStress.H:1424
const amrex::Real tol_z
Definition: ERF_MOSTStress.H:1593
const amrex::Real alpha
Definition: ERF_MOSTStress.H:1595
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1596
const amrex::Real tol
Definition: ERF_MOSTStress.H:1589
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1588
Definition: ERF_MOSTStress.H:1785
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1909
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:1814
most_data mdata
Definition: ERF_MOSTStress.H:1907
const amrex::Real tol
Definition: ERF_MOSTStress.H:1910
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1911
bool spec_qflux
Definition: ERF_MOSTStress.H:1908
surface_temp_donelan(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Definition: ERF_MOSTStress.H:1793
Definition: ERF_MOSTStress.H:1604
surface_temp_mod_charnock(amrex::Real Tflux, amrex::Real Qvflux, amrex::Real depth, bool cons_qflux)
Definition: ERF_MOSTStress.H:1613
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1769
const amrex::Real tol
Definition: ERF_MOSTStress.H:1770
const amrex::Real tol_z
Definition: ERF_MOSTStress.H:1774
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1777
most_data mdata
Definition: ERF_MOSTStress.H:1767
bool spec_qflux
Definition: ERF_MOSTStress.H:1768
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:1637
const amrex::Real alpha
Definition: ERF_MOSTStress.H:1776
Definition: ERF_MOSTStress.H:1919
const amrex::Real eps
Definition: ERF_MOSTStress.H:2054
surface_temp_wave_coupled(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Definition: ERF_MOSTStress.H:1927
const amrex::Real tol
Definition: ERF_MOSTStress.H:2050
most_data mdata
Definition: ERF_MOSTStress.H:2047
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< 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:1948
const amrex::Real z0_eps
Definition: ERF_MOSTStress.H:2056
similarity_funs sfuns
Definition: ERF_MOSTStress.H:2049
bool spec_qflux
Definition: ERF_MOSTStress.H:2048
const amrex::Real z0_max
Definition: ERF_MOSTStress.H:2057
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:2058
Definition: ERF_MOSTStress.H:1264
surface_temp(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Definition: ERF_MOSTStress.H:1272
similarity_funs sfuns
Definition: ERF_MOSTStress.H:1403
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void iterate_flux(const int &i, const int &j, const int &k, const int &max_iters, const amrex::Array4< const amrex::Real > &zref_arr, const amrex::Array4< const amrex::Real > &z0_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &tvm_arr, const amrex::Array4< const amrex::Real > &qvm_arr, const amrex::Array4< amrex::Real > &u_star_arr, const amrex::Array4< amrex::Real > &w_star_arr, const amrex::Array4< amrex::Real > &t_star_arr, const amrex::Array4< amrex::Real > &q_star_arr, const amrex::Array4< amrex::Real > &t_surf_arr, const amrex::Array4< amrex::Real > &q_surf_arr, const amrex::Array4< amrex::Real > &olen_arr, const amrex::Array4< amrex::Real > &pblh_arr, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &, const amrex::Array4< amrex::Real > &) const
Definition: ERF_MOSTStress.H:1293
const amrex::Real tol
Definition: ERF_MOSTStress.H:1404
bool spec_qflux
Definition: ERF_MOSTStress.H:1402
const amrex::Real alpha
Definition: ERF_MOSTStress.H:1405
const amrex::Real WSMIN
Definition: ERF_MOSTStress.H:1406
most_data mdata
Definition: ERF_MOSTStress.H:1401