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