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