5 #ifndef ERF_EBMOSTStress_H
6 #define ERF_EBMOSTStress_H
35 const amrex::Array4<const amrex::Real>& zref_arr,
36 const amrex::Array4<const amrex::Real>& z0_arr,
37 const amrex::Array4<const amrex::Real>& umm_arr,
38 const amrex::Array4<const amrex::Real>& ,
39 const amrex::Array4<const amrex::Real>& ,
40 const amrex::Array4<const amrex::Real>& ,
41 const amrex::Array4<amrex::Real>& u_star_arr,
42 const amrex::Array4<amrex::Real>& ,
43 const amrex::Array4<amrex::Real>& t_star_arr,
44 const amrex::Array4<amrex::Real>& q_star_arr,
45 const amrex::Array4<amrex::Real>& ,
46 const amrex::Array4<amrex::Real>& ,
47 const amrex::Array4<amrex::Real>& olen_arr,
48 const amrex::Array4<amrex::Real>& ,
49 const amrex::Array4<amrex::Real>& ,
50 const amrex::Array4<amrex::Real>& ,
51 const amrex::Array4<amrex::Real>& )
const
54 u_star_arr(i,j,k) =
mdata.
kappa * umm_arr(i,j,0) / std::log(zref_arr(i,j,0) / z0_arr(i,j,k));
55 t_star_arr(i,j,k) =
zero;
56 q_star_arr(i,j,k) =
zero;
91 const amrex::Array4<const amrex::Real>& zref_arr,
92 const amrex::Array4<const amrex::Real>& z0_arr,
93 const amrex::Array4<const amrex::Real>& umm_arr,
94 const amrex::Array4<const amrex::Real>& tm_arr,
95 const amrex::Array4<const amrex::Real>& tvm_arr,
96 const amrex::Array4<const amrex::Real>& qvm_arr,
97 const amrex::Array4<amrex::Real>& u_star_arr,
98 const amrex::Array4<amrex::Real>& w_star_arr,
99 const amrex::Array4<amrex::Real>& t_star_arr,
100 const amrex::Array4<amrex::Real>& q_star_arr,
101 const amrex::Array4<amrex::Real>& t_surf_arr,
102 const amrex::Array4<amrex::Real>& q_surf_arr,
103 const amrex::Array4<amrex::Real>& olen_arr,
104 const amrex::Array4<amrex::Real>& pblh_arr,
105 const amrex::Array4<amrex::Real>& ,
106 const amrex::Array4<amrex::Real>& ,
107 const amrex::Array4<amrex::Real>& )
const
125 zeta = zref / olen_arr(i,j,k);
129 if (q_surf_arr(i,j,k) >
zero) {
130 qv_s = q_surf_arr(i,j,k);
134 qv_s = qvm_arr(i,j,0);
136 qv_a = qvm_arr(i,j,0);
146 -ustar *
mdata.
kappa * (qv_a - qv_s) / (C - psi_h);
148 w_star_arr(i,j,k) =
calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,0));
150 umm = std::sqrt(umm_arr(i,j,0)*umm_arr(i,j,0) + wstar*wstar);
151 umm = std::max(umm,
WSMIN);
158 ( (thv_a - thv_s) / (umm * umm) );
179 }
while ( (std::abs(zeta - zeta_old) >
tol) && (iter <= max_iters) );
181 "Maximum number of MOST iterations reached.");
184 olen_arr(i,j,k) = zref / zeta;
185 u_star_arr(i,j,k) =
mdata.
kappa * umm / (C - psi_m);
186 t_star_arr(i,j,k) =
mdata.
kappa * (tm_arr(i,j,0) - t_surf_arr(i,j,k)) / (C - psi_h);
189 (u_star_arr(i,j,k) *
mdata.
kappa) + qvm_arr(i,j,0);
192 q_star_arr(i,j,k) =
mdata.
kappa * (qvm_arr(i,j,0) - q_surf_arr(i,j,k)) / (C - psi_h);
230 const int& max_iters,
231 const amrex::Array4<const amrex::Real>& zref_arr,
232 const amrex::Array4<const amrex::Real>& z0_arr,
233 const amrex::Array4<const amrex::Real>& umm_arr,
234 const amrex::Array4<const amrex::Real>& tm_arr,
235 const amrex::Array4<const amrex::Real>& tvm_arr,
236 const amrex::Array4<const amrex::Real>& qvm_arr,
237 const amrex::Array4<amrex::Real>& u_star_arr,
238 const amrex::Array4<amrex::Real>& w_star_arr,
239 const amrex::Array4<amrex::Real>& t_star_arr,
240 const amrex::Array4<amrex::Real>& q_star_arr,
241 const amrex::Array4<amrex::Real>& t_surf_arr,
242 const amrex::Array4<amrex::Real>& q_surf_arr,
243 const amrex::Array4<amrex::Real>& olen_arr,
244 const amrex::Array4<amrex::Real>& pblh_arr,
245 const amrex::Array4<amrex::Real>& ,
246 const amrex::Array4<amrex::Real>& ,
247 const amrex::Array4<amrex::Real>& )
const
261 u_star_arr(i,j,k) =
mdata.
kappa * umm / std::log(zref / z0_arr(i,j,k));
263 Olen = olen_arr(i,j,k);
269 ustar = u_star_arr(i,j,k);
271 -(qvm_arr(i,j,0) - q_surf_arr(i,j,k)) * ustar *
mdata.
kappa /
272 (std::log(zref / z0_arr(i,j,k)) - psi_h);
276 w_star_arr(i,j,k) =
calc_wstar(tflux, pblh_arr(i,j,k), tvm_arr(i,j,0));
278 umm = std::sqrt(umm_arr(i,j,0)*umm_arr(i,j,0) + wstar*wstar);
279 umm = std::max(umm,
WSMIN);
285 u_star_arr(i,j,k) =
mdata.
kappa * umm / (std::log(zref / z0_arr(i,j,k)) - psi_m);
287 }
while ((std::abs(u_star_arr(i,j,k) - ustar) >
tol) && iter <= max_iters);
289 "Maximum number of MOST iterations reached.");
292 olen_arr(i,j,k) = Olen;
294 (u_star_arr(i,j,k) *
mdata.
kappa) + tm_arr(i,j,0);
298 (u_star_arr(i,j,k) *
mdata.
kappa) + qvm_arr(i,j,0);
301 q_star_arr(i,j,k) =
mdata.
kappa * (qvm_arr(i,j,0) - q_surf_arr(i,j,k)) /
302 (std::log(zref / z0_arr(i,j,k)) - psi_h);
327 const amrex::Array4<const amrex::Real>& cons_arr,
328 const amrex::Array4<const amrex::Real>& velx_arr,
329 const amrex::Array4<const amrex::Real>& vely_arr,
330 const amrex::Array4<const amrex::Real>& umm_arr,
331 const amrex::Array4<const amrex::Real>& qvm_arr,
332 const amrex::Array4<const amrex::Real>& u_star_arr,
333 const amrex::Array4<const amrex::Real>& q_star_arr,
334 const amrex::Array4<const amrex::Real>& q_surf_arr,
335 const amrex::Array4<const amrex::Real>& u_vfrac_arr,
336 const amrex::Array4<const amrex::Real>& v_vfrac_arr)
const
342 amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
344 (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
345 / u_vfrac_sum :
zero;
348 amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
350 (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
351 / v_vfrac_sum :
zero;
358 wsp_mean = std::max(wsp_mean,
WSMIN);
366 -
rho*qstar*ustar*(num1+num2)/((qv_mean-qv_surf)*wsp_mean) :
zero;
378 const amrex::Array4<const amrex::Real>& cons_arr,
379 const amrex::Array4<const amrex::Real>& velx_arr,
380 const amrex::Array4<const amrex::Real>& vely_arr,
381 const amrex::Array4<const amrex::Real>& velz_arr,
382 const amrex::Array4<const amrex::Real>& umm_arr,
383 const amrex::Array4<const amrex::Real>& tm_arr,
384 const amrex::Array4<const amrex::Real>& u_star_arr,
385 const amrex::Array4<const amrex::Real>& t_star_arr,
386 const amrex::Array4<const amrex::Real>& t_surf_arr,
387 const amrex::Array4<const amrex::Real>& u_vfrac_arr,
388 const amrex::Array4<const amrex::Real>& v_vfrac_arr,
389 const amrex::Array4<const amrex::Real>& w_vfrac_arr,
390 const amrex::Array4<const amrex::Real>& bnorm_arr)
const
396 amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
398 (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
399 / u_vfrac_sum :
zero;
402 amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
404 (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
405 / v_vfrac_sum :
zero;
408 amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1);
410 (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1))
411 / w_vfrac_sum :
zero;
428 wsp_mean = std::max(wsp_mean,
WSMIN);
431 amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent);
437 -
rho*tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) :
zero;
449 const amrex::Array4<const amrex::Real>& cons_arr,
450 const amrex::Array4<const amrex::Real>& velx_arr,
451 const amrex::Array4<const amrex::Real>& vely_arr,
452 const amrex::Array4<const amrex::Real>& velz_arr,
453 const amrex::Array4<const amrex::Real>& umm_arr,
454 const amrex::Array4<const amrex::Real>& um_arr,
455 const amrex::Array4<const amrex::Real>& u_star_arr,
456 const amrex::Array4<const amrex::Real>& u_vfrac_arr,
457 const amrex::Array4<const amrex::Real>& v_vfrac_arr,
458 const amrex::Array4<const amrex::Real>& w_vfrac_arr,
459 const amrex::Array4<const amrex::Real>& cc_vfrac_arr,
460 const amrex::Array4<const amrex::EBCellFlag>& cc_flag_arr,
461 const amrex::Array4<const amrex::Real>& bnorm_arr,
469 velx = velx_arr(i,j,k);
472 amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k) +
473 v_vfrac_arr(i-1,j,k) + v_vfrac_arr(i-1,j+1,k);
474 vely = (v_vfrac_sum >
eps) ?
475 (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k) +
476 vely_arr(i-1,j,k) * v_vfrac_arr(i-1,j,k) + vely_arr(i-1,j+1,k) * v_vfrac_arr(i-1,j+1,k))
477 / v_vfrac_sum :
zero;
480 amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
481 w_vfrac_arr(i-1,j,k) + w_vfrac_arr(i-1,j,k+1);
483 (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
484 velz_arr(i-1,j,k) * w_vfrac_arr(i-1,j,k) + velz_arr(i-1,j,k+1) * w_vfrac_arr(i-1,j,k+1))
485 / w_vfrac_sum :
zero;
494 velx_tangent = velx - v_dot_n * nx;
495 vely_tangent = vely - v_dot_n * ny;
498 amrex::Real cc_vfrac_sum = cc_vfrac_arr(i-1,j,k) + cc_vfrac_arr(i,j,k);
499 rho = (cc_vfrac_sum >
eps) ?
500 (cons_arr(i-1,j,k,
Rho_comp) * cc_vfrac_arr(i-1,j,k) + cons_arr(i,j,k,
Rho_comp) * cc_vfrac_arr(i,j,k))
501 / cc_vfrac_sum :
zero;
504 bool low_valid = cc_flag_arr(i-1,j,k).isSingleValued();
505 bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
507 if (low_valid && high_valid) {
508 ustar =
myhalf * (u_star_arr(i-1,j,k) + u_star_arr(i,j,k));
509 wsp_mean =
myhalf * (umm_arr(i-1,j,0) + umm_arr(i,j,0));
510 }
else if (low_valid) {
511 ustar = u_star_arr(i-1,j,k);
512 wsp_mean = umm_arr(i-1,j,0);
513 }
else if (high_valid) {
514 ustar = u_star_arr(i,j,k);
515 wsp_mean = umm_arr(i,j,0);
521 }
else if (idir == 1) {
523 vely = vely_arr(i,j,k);
526 amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k) +
527 u_vfrac_arr(i,j-1,k) + u_vfrac_arr(i+1,j-1,k);
528 velx = (u_vfrac_sum >
eps) ?
529 (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k) +
530 velx_arr(i,j-1,k) * u_vfrac_arr(i,j-1,k) + velx_arr(i+1,j-1,k) * u_vfrac_arr(i+1,j-1,k))
531 / u_vfrac_sum :
zero;
534 amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
535 w_vfrac_arr(i,j-1,k) + w_vfrac_arr(i,j-1,k+1);
537 (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
538 velz_arr(i,j-1,k) * w_vfrac_arr(i,j-1,k) + velz_arr(i,j-1,k+1) * w_vfrac_arr(i,j-1,k+1))
539 / w_vfrac_sum :
zero;
548 velx_tangent = velx - v_dot_n * nx;
549 vely_tangent = vely - v_dot_n * ny;
552 amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j-1,k) + cc_vfrac_arr(i,j,k);
553 rho = (cc_vfrac_sum >
eps) ?
554 (cons_arr(i,j-1,k,
Rho_comp) * cc_vfrac_arr(i,j-1,k) + cons_arr(i,j,k,
Rho_comp) * cc_vfrac_arr(i,j,k))
555 / cc_vfrac_sum :
zero;
558 bool low_valid = cc_flag_arr(i,j-1,k).isSingleValued();
559 bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
561 if (low_valid && high_valid) {
562 ustar =
myhalf * (u_star_arr(i,j-1,k) + u_star_arr(i,j,k));
563 wsp_mean =
myhalf * (umm_arr(i,j-1,0) + umm_arr(i,j,0));
564 }
else if (low_valid) {
565 ustar = u_star_arr(i,j-1,k);
566 wsp_mean = umm_arr(i,j-1,0);
567 }
else if (high_valid) {
568 ustar = u_star_arr(i,j,k);
569 wsp_mean = umm_arr(i,j,0);
578 amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k-1) + u_vfrac_arr(i+1,j,k-1) +
579 u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
580 velx = (u_vfrac_sum >
eps) ?
581 (velx_arr(i,j,k-1) * u_vfrac_arr(i,j,k-1) + velx_arr(i+1,j,k-1) * u_vfrac_arr(i+1,j,k-1) +
582 velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
583 / u_vfrac_sum :
zero;
586 amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k-1) + v_vfrac_arr(i,j+1,k-1) +
587 v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
588 vely = (v_vfrac_sum >
eps) ?
589 (vely_arr(i,j,k-1) * v_vfrac_arr(i,j,k-1) + vely_arr(i,j+1,k-1) * v_vfrac_arr(i,j+1,k-1) +
590 vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
591 / v_vfrac_sum :
zero;
603 velx_tangent = velx - v_dot_n * nx;
604 vely_tangent = vely - v_dot_n * ny;
607 amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j,k-1) + cc_vfrac_arr(i,j,k);
608 rho = (cc_vfrac_sum >
eps) ?
609 (cons_arr(i,j,k-1,
Rho_comp) * cc_vfrac_arr(i,j,k-1) + cons_arr(i,j,k,
Rho_comp) * cc_vfrac_arr(i,j,k))
610 / cc_vfrac_sum :
zero;
613 bool low_valid = cc_flag_arr(i,j,k-1).isSingleValued();
614 bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
616 if (low_valid && high_valid) {
617 ustar =
myhalf * (u_star_arr(i,j,k-1) + u_star_arr(i,j,k));
618 wsp_mean =
myhalf * (umm_arr(i,j,0) + umm_arr(i,j,0));
619 }
else if (low_valid) {
620 ustar = u_star_arr(i,j,k-1);
621 wsp_mean = umm_arr(i,j,0);
622 }
else if (high_valid) {
623 ustar = u_star_arr(i,j,k);
624 wsp_mean = umm_arr(i,j,0);
631 wsp_mean = std::max(wsp_mean,
WSMIN);
638 amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent);
640 amrex::Real num2 = wsp_mean * (velx_tangent-umean);
643 amrex::Real stressx = -
rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
655 const amrex::Array4<const amrex::Real>& cons_arr,
656 const amrex::Array4<const amrex::Real>& velx_arr,
657 const amrex::Array4<const amrex::Real>& vely_arr,
658 const amrex::Array4<const amrex::Real>& velz_arr,
659 const amrex::Array4<const amrex::Real>& umm_arr,
660 const amrex::Array4<const amrex::Real>& vm_arr,
661 const amrex::Array4<const amrex::Real>& u_star_arr,
662 const amrex::Array4<const amrex::Real>& u_vfrac_arr,
663 const amrex::Array4<const amrex::Real>& v_vfrac_arr,
664 const amrex::Array4<const amrex::Real>& w_vfrac_arr,
665 const amrex::Array4<const amrex::Real>& cc_vfrac_arr,
666 const amrex::Array4<const amrex::EBCellFlag>& cc_flag_arr,
667 const amrex::Array4<const amrex::Real>& bnorm_arr,
676 amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k) +
677 v_vfrac_arr(i-1,j,k) + v_vfrac_arr(i-1,j+1,k);
678 vely = (v_vfrac_sum >
eps) ?
679 (vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k) +
680 vely_arr(i-1,j,k) * v_vfrac_arr(i-1,j,k) + vely_arr(i-1,j+1,k) * v_vfrac_arr(i-1,j+1,k))
681 / v_vfrac_sum :
zero;
683 velx = velx_arr(i,j,k);
686 amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
687 w_vfrac_arr(i-1,j,k) + w_vfrac_arr(i-1,j,k+1);
689 (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
690 velz_arr(i-1,j,k) * w_vfrac_arr(i-1,j,k) + velz_arr(i-1,j,k+1) * w_vfrac_arr(i-1,j,k+1))
691 / w_vfrac_sum :
zero;
700 velx_tangent = velx - v_dot_n * nx;
701 vely_tangent = vely - v_dot_n * ny;
704 amrex::Real cc_vfrac_sum = cc_vfrac_arr(i-1,j,k) + cc_vfrac_arr(i,j,k);
705 rho = (cc_vfrac_sum >
eps) ?
706 (cons_arr(i-1,j,k,
Rho_comp) * cc_vfrac_arr(i-1,j,k) + cons_arr(i,j,k,
Rho_comp) * cc_vfrac_arr(i,j,k))
707 / cc_vfrac_sum :
zero;
710 bool low_valid = cc_flag_arr(i-1,j,k).isSingleValued();
711 bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
713 if (low_valid && high_valid) {
714 ustar =
myhalf * (u_star_arr(i-1,j,k) + u_star_arr(i,j,k));
715 wsp_mean =
myhalf * (umm_arr(i-1,j,0) + umm_arr(i,j,0));
716 }
else if (low_valid) {
717 ustar = u_star_arr(i-1,j,k);
718 wsp_mean = umm_arr(i-1,j,0);
719 }
else if (high_valid) {
720 ustar = u_star_arr(i,j,k);
721 wsp_mean = umm_arr(i,j,0);
727 }
else if (idir == 1) {
730 amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k) +
731 u_vfrac_arr(i,j-1,k) + u_vfrac_arr(i+1,j-1,k);
732 velx = (u_vfrac_sum >
eps) ?
733 (velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k) +
734 velx_arr(i,j-1,k) * u_vfrac_arr(i,j-1,k) + velx_arr(i+1,j-1,k) * u_vfrac_arr(i+1,j-1,k))
735 / u_vfrac_sum :
zero;
737 vely = vely_arr(i,j,k);
740 amrex::Real w_vfrac_sum = w_vfrac_arr(i,j,k) + w_vfrac_arr(i,j,k+1) +
741 w_vfrac_arr(i,j-1,k) + w_vfrac_arr(i,j-1,k+1);
743 (velz_arr(i,j,k) * w_vfrac_arr(i,j,k) + velz_arr(i,j,k+1) * w_vfrac_arr(i,j,k+1) +
744 velz_arr(i,j-1,k) * w_vfrac_arr(i,j-1,k) + velz_arr(i,j-1,k+1) * w_vfrac_arr(i,j-1,k+1))
745 / w_vfrac_sum :
zero;
754 velx_tangent = velx - v_dot_n * nx;
755 vely_tangent = vely - v_dot_n * ny;
758 amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j-1,k) + cc_vfrac_arr(i,j,k);
759 rho = (cc_vfrac_sum >
eps) ?
760 (cons_arr(i,j-1,k,
Rho_comp) * cc_vfrac_arr(i,j-1,k) + cons_arr(i,j,k,
Rho_comp) * cc_vfrac_arr(i,j,k))
761 / cc_vfrac_sum :
zero;
764 bool low_valid = cc_flag_arr(i,j-1,k).isSingleValued();
765 bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
767 if (low_valid && high_valid) {
768 ustar =
myhalf * (u_star_arr(i,j-1,k) + u_star_arr(i,j,k));
769 wsp_mean =
myhalf * (umm_arr(i,j-1,0) + umm_arr(i,j,0));
770 }
else if (low_valid) {
771 ustar = u_star_arr(i,j-1,k);
772 wsp_mean = umm_arr(i,j-1,0);
773 }
else if (high_valid) {
774 ustar = u_star_arr(i,j,k);
775 wsp_mean = umm_arr(i,j,0);
784 amrex::Real u_vfrac_sum = u_vfrac_arr(i,j,k-1) + u_vfrac_arr(i+1,j,k-1) +
785 u_vfrac_arr(i,j,k) + u_vfrac_arr(i+1,j,k);
786 velx = (u_vfrac_sum >
eps) ?
787 (velx_arr(i,j,k-1) * u_vfrac_arr(i,j,k-1) + velx_arr(i+1,j,k-1) * u_vfrac_arr(i+1,j,k-1) +
788 velx_arr(i,j,k) * u_vfrac_arr(i,j,k) + velx_arr(i+1,j,k) * u_vfrac_arr(i+1,j,k))
789 / u_vfrac_sum :
zero;
792 amrex::Real v_vfrac_sum = v_vfrac_arr(i,j,k-1) + v_vfrac_arr(i,j+1,k-1) +
793 v_vfrac_arr(i,j,k) + v_vfrac_arr(i,j+1,k);
794 vely = (v_vfrac_sum >
eps) ?
795 (vely_arr(i,j,k-1) * v_vfrac_arr(i,j,k-1) + vely_arr(i,j+1,k-1) * v_vfrac_arr(i,j+1,k-1) +
796 vely_arr(i,j,k) * v_vfrac_arr(i,j,k) + vely_arr(i,j+1,k) * v_vfrac_arr(i,j+1,k))
797 / v_vfrac_sum :
zero;
809 velx_tangent = velx - v_dot_n * nx;
810 vely_tangent = vely - v_dot_n * ny;
813 amrex::Real cc_vfrac_sum = cc_vfrac_arr(i,j,k-1) + cc_vfrac_arr(i,j,k);
814 rho = (cc_vfrac_sum >
eps) ?
815 (cons_arr(i,j,k-1,
Rho_comp) * cc_vfrac_arr(i,j,k-1) + cons_arr(i,j,k,
Rho_comp) * cc_vfrac_arr(i,j,k))
816 / cc_vfrac_sum :
zero;
819 bool low_valid = cc_flag_arr(i,j,k-1).isSingleValued();
820 bool high_valid = cc_flag_arr(i,j,k).isSingleValued();
822 if (low_valid && high_valid) {
823 ustar =
myhalf * (u_star_arr(i,j,k-1) + u_star_arr(i,j,k));
824 wsp_mean =
myhalf * (umm_arr(i,j,0) + umm_arr(i,j,0));
825 }
else if (low_valid) {
826 ustar = u_star_arr(i,j,k-1);
827 wsp_mean = umm_arr(i,j,0);
828 }
else if (high_valid) {
829 ustar = u_star_arr(i,j,k);
830 wsp_mean = umm_arr(i,j,0);
837 wsp_mean = std::max(wsp_mean,
WSMIN);
844 amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent);
846 amrex::Real num2 = wsp_mean * (vely_tangent-vmean);
849 amrex::Real stressy = -
rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
855 #ifdef AMREX_USE_FLOAT
constexpr amrex::Real epsv
Definition: ERF_Constants.H:53
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
@ num
Definition: ERF_DataStruct.H:27
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
rho
Definition: ERF_InitCustomPert_Bubble.H:107
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real calc_wstar(const amrex::Real &ust, const amrex::Real &tst, const amrex::Real &qst, const amrex::Real &pblh, const amrex::Real &th, const amrex::Real &thv, const amrex::Real &qv=amrex::Real(0))
Definition: ERF_Wstar.H:13
@ theta
Definition: ERF_SLM.H:20
@ qv
Definition: ERF_Kessler.H:30
@ den
Definition: ERF_AdvanceWSM6.cpp:109
EB surface-layer model for adiabatic constant-roughness fluxes.
Definition: ERF_EBMOSTStress.H:14
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
Populate neutral MOST flux scales for an adiabatic EB surface.
Definition: ERF_EBMOSTStress.H:31
adiabatic_eb(amrex::Real Tflux, amrex::Real Qvflux)
Construct an adiabatic EB flux functor.
Definition: ERF_EBMOSTStress.H:20
most_data mdata
Definition: ERF_EBMOSTStress.H:60
similarity_funs sfuns
Definition: ERF_EBMOSTStress.H:61
EB implementation of the Moeng surface-flux formulation.
Definition: ERF_EBMOSTStress.H:316
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_u_flux(int i, int j, 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 > &velz_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 amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr, const amrex::Array4< const amrex::Real > &w_vfrac_arr, const amrex::Array4< const amrex::Real > &cc_vfrac_arr, const amrex::Array4< const amrex::EBCellFlag > &cc_flag_arr, const amrex::Array4< const amrex::Real > &bnorm_arr, int idir=0) const
Compute the EB x-momentum stress on a staggered face.
Definition: ERF_EBMOSTStress.H:446
const amrex::Real WSMIN
Definition: ERF_EBMOSTStress.H:860
moeng_flux_eb()
Construct a Moeng EB flux functor.
Definition: ERF_EBMOSTStress.H:318
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 > &velz_arr, const amrex::Array4< const amrex::Real > &umm_arr, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::Array4< const amrex::Real > &u_star_arr, const amrex::Array4< const amrex::Real > &t_star_arr, const amrex::Array4< const amrex::Real > &t_surf_arr, const amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr, const amrex::Array4< const amrex::Real > &w_vfrac_arr, const amrex::Array4< const amrex::Real > &bnorm_arr) const
Compute the EB temperature flux at a cell-centered cut-cell surface.
Definition: ERF_EBMOSTStress.H:375
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real compute_v_flux(int i, int j, 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 > &velz_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 amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr, const amrex::Array4< const amrex::Real > &w_vfrac_arr, const amrex::Array4< const amrex::Real > &cc_vfrac_arr, const amrex::Array4< const amrex::EBCellFlag > &cc_flag_arr, const amrex::Array4< const amrex::Real > &bnorm_arr, int idir=0) const
Compute the EB y-momentum stress on a staggered face.
Definition: ERF_EBMOSTStress.H:652
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 amrex::Array4< const amrex::Real > &u_vfrac_arr, const amrex::Array4< const amrex::Real > &v_vfrac_arr) const
Compute the EB moisture flux at a cell-centered cut-cell surface.
Definition: ERF_EBMOSTStress.H:324
const amrex::Real eps
Definition: ERF_EBMOSTStress.H:858
Definition: ERF_MOSTStress.H:13
amrex::Real surf_moist_flux
Moisture flux.
Definition: ERF_MOSTStress.H:19
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
const amrex::Real Bjr_beta
Definition: ERF_MOSTStress.H:32
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:40
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:105
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h2(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:77
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:124
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m2(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:52
EB surface-layer model with prescribed surface fluxes and constant roughness.
Definition: ERF_EBMOSTStress.H:207
const amrex::Real WSMIN
Definition: ERF_EBMOSTStress.H:311
similarity_funs sfuns
Definition: ERF_EBMOSTStress.H:309
const amrex::Real tol
Definition: ERF_EBMOSTStress.H:310
most_data mdata
Definition: ERF_EBMOSTStress.H:307
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
Iterate MOST stability functions and update implied surface values.
Definition: ERF_EBMOSTStress.H:227
bool spec_qflux
Definition: ERF_EBMOSTStress.H:308
surface_flux_eb(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Construct a prescribed-flux EB functor.
Definition: ERF_EBMOSTStress.H:214
EB surface-layer model with prescribed surface temperature and constant roughness.
Definition: ERF_EBMOSTStress.H:67
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
Iterate MOST stability functions and update surface flux scales.
Definition: ERF_EBMOSTStress.H:87
const amrex::Real tol
Definition: ERF_EBMOSTStress.H:200
const amrex::Real alpha
Definition: ERF_EBMOSTStress.H:201
const amrex::Real WSMIN
Definition: ERF_EBMOSTStress.H:202
similarity_funs sfuns
Definition: ERF_EBMOSTStress.H:199
surface_temp_eb(amrex::Real Tflux, amrex::Real Qvflux, bool cons_qflux)
Construct a prescribed-temperature EB flux functor.
Definition: ERF_EBMOSTStress.H:74
bool spec_qflux
Definition: ERF_EBMOSTStress.H:198
most_data mdata
Definition: ERF_EBMOSTStress.H:197