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