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