Fix bomex MOST. (erf-model#1627)

Co-authored-by: Aaron Lattanzi <[email protected]> Co-authored-by: Ann Almgren <[email protected]>
wiersema1 · May 24, 2024 · 120f71f · 120f71f
1 parent 517e56c
commit 120f71f
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 8 deletions.
diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H
@@ -1198,7 +1198,7 @@ struct custom_flux
         je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
-        dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz);
+        dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) - moflux*rho/eta*deltaz;
 
         return moflux;
     }
@@ -1238,7 +1238,7 @@ struct custom_flux
 
         amrex::Real ustar   = u_star_arr(ic,jc,zlo);
         amrex::Real tstar   = t_star_arr(ic,jc,zlo);
-        amrex::Real moflux  = (std::abs(tstar) > eps) ? -tstar*ustar : 0.0;
+        amrex::Real moflux  = (std::abs(tstar) > eps) ? -rho * tstar : 0.0;
         amrex::Real deltaz  = dz * (zlo - k);
 
         if (exp_most) {
@@ -1254,10 +1254,7 @@ struct custom_flux
             je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
             eta   = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
             eta   = amrex::max(eta,eta_eps);
-            // Note: Kh = eta/rho
-            //      hfx = -Kh dT/dz  ==> +ve hfx corresponds to heating from the surface
-            // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel
-            dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz);
+            dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) - moflux*rho/eta*deltaz;
         }
 
         return moflux;

diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp
@@ -191,10 +191,11 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22,
           // Calculate SFS quantities
           // - dissipation
           Real Ce;
-          if ((l_C_e_wall > 0) && (k==0))
+          if ((l_C_e_wall > 0) && (k==0)) {
               Ce = l_C_e_wall;
-          else
+          } else {
               Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf;
+          }
           diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length;
           // - heat flux
           //   (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will