Add Brunt-Vaisala freq into Smag model

Docs/sphinx_doc/MOST.rst

@@ -141,6 +141,7 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in
    .. math::
 
       (\rho \theta)_{z} = \frac{(\rho \theta)_{i,j,1} - (\rho \theta)_{i,j,0}}{\Delta z}
+
       (\rho \theta)_{i,j,-n} = (\rho \theta)_{i,j,0} - (\rho \theta)_{z} n \Delta z .
 
    Finally, it must be noted that complex terrain will modify the surface normal and tangent vectors. Consequently, the MOST implementation with terrain will require local vector rotations. While the ERF dycore accounts for

Exec/DevTests/Bomex/input_SAM

@@ -19,6 +19,7 @@ erf.most.tstar  = 8.0e-3 # theta flux
 erf.most.qstar  = 5.2e-5 # qv    flux
 erf.most.z0     = 0.1
 erf.most.zref   = 20.0
+erf.use_explicit_most = true
 
 # NOTE: This should have a qv grad too (use hoextrapcc?!)
 zhi.type = "SlipWall"
@@ -58,7 +59,7 @@ erf.dryscal_vert_adv_type    = WENOZ5
 erf.moistscal_horiz_adv_type = WENOZ5
 erf.moistscal_vert_adv_type  = WENOZ5
 
-erf.moisture_model  = "SAM"
+erf.moisture_model  = "Kessler"
 
 erf.molec_diff_type = "None"
 

Source/BoundaryConditions/ABLMost.H

@@ -252,6 +252,7 @@ public:
                      amrex::MultiFab* xzmom_flux, amrex::MultiFab* zxmom_flux,
                      amrex::MultiFab* yzmom_flux, amrex::MultiFab* zymom_flux,
                      amrex::MultiFab* heat_flux,
+                     amrex::MultiFab* qv_flux,
                      amrex::MultiFab* z_phys);
 
     template<typename FluxCalc>
@@ -261,6 +262,7 @@ public:
                       amrex::MultiFab* xzmom_flux, amrex::MultiFab* zxmom_flux,
                       amrex::MultiFab* yzmom_flux, amrex::MultiFab* zymom_flux,
                       amrex::MultiFab* heat_flux,
+                      amrex::MultiFab* qv_flux,
                       amrex::MultiFab* z_phys,
                       const amrex::Real& dz_no_terrain,
                       const FluxCalc& flux_comp);

Source/BoundaryConditions/ABLMost.cpp

@@ -191,6 +191,7 @@ ABLMost::impose_most_bcs (const int& lev,
                           MultiFab* xzmom_flux, MultiFab* zxmom_flux,
                           MultiFab* yzmom_flux, MultiFab* zymom_flux,
                           MultiFab* heat_flux,
+                          MultiFab* qv_flux,
                           MultiFab* z_phys)
 {
     const int klo = 0;
@@ -200,20 +201,23 @@ ABLMost::impose_most_bcs (const int& lev,
                          xzmom_flux, zxmom_flux,
                          yzmom_flux, zymom_flux,
                          heat_flux,
+                         qv_flux,
                          z_phys, m_geom[lev].CellSize(2), flux_comp);
     } else if (flux_type == FluxCalcType::DONELAN) {
         donelan_flux flux_comp(klo);
         compute_most_bcs(lev, mfs,
                          xzmom_flux, zxmom_flux,
                          yzmom_flux, zymom_flux,
                          heat_flux,
+                         qv_flux,
                          z_phys, m_geom[lev].CellSize(2), flux_comp);
     } else {
         custom_flux flux_comp(klo);
         compute_most_bcs(lev, mfs,
                          xzmom_flux, zxmom_flux,
                          yzmom_flux, zymom_flux,
                          heat_flux,
+                         qv_flux,
                          z_phys, m_geom[lev].CellSize(2), flux_comp);
     }
 }
@@ -234,6 +238,7 @@ ABLMost::compute_most_bcs (const int& lev,
                            MultiFab* xzmom_flux, MultiFab* zxmom_flux,
                            MultiFab* yzmom_flux, MultiFab* zymom_flux,
                            MultiFab* heat_flux,
+                           MultiFab* qv_flux,
                            MultiFab* z_phys,
                            const Real& dz_no_terrain,
                            const FluxCalc& flux_comp)
@@ -242,7 +247,6 @@ ABLMost::compute_most_bcs (const int& lev,
     const int icomp = 0;
     for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
     {
-        // TODO: No LSM lateral ghost cells, should this change?
         // Valid CC box
         Box vbx = mfi.validbox(); vbx.makeSlab(2,klo-1);
 
@@ -262,6 +266,7 @@ ABLMost::compute_most_bcs (const int& lev,
         auto t31_arr = (zxmom_flux && m_exp_most) ? zxmom_flux->array(mfi) : Array4<Real>{};
         auto t32_arr = (zymom_flux && m_exp_most) ? zymom_flux->array(mfi) : Array4<Real>{};
         auto hfx_arr = (m_exp_most) ? heat_flux->array(mfi) : Array4<Real>{};
+        auto qfx_arr = (m_exp_most && use_moisture ) ? qv_flux->array(mfi) : Array4<Real>{};
 
         const auto  eta_arr  = m_eddyDiffs_lev[lev]->array(mfi);
 
@@ -292,6 +297,11 @@ ABLMost::compute_most_bcs (const int& lev,
         auto lsm_flux_arr = (m_lsm_flux_lev[lev][0]) ? m_lsm_flux_lev[lev][0]->array(mfi) :
                                                        Array4<Real> {};
 
+        // NOTE: The returned fluxes (tflux/qflux/stressx/stressy) are the diffusive fluxes
+        //       defined as F_d = -K \partial \phi / \partial z. This is consistent with the
+        //       computation of the diffusion source term as well as the buoyancy source for
+        //       turbulence modeling (e.g. Deardorff)
+
         for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx)
         {
             const Box& bx = (*mfs[var_idx])[mfi].box();
@@ -305,7 +315,7 @@ ABLMost::compute_most_bcs (const int& lev,
                 {
                     Real dz  = (zphys_arr) ? ( zphys_arr(i,j,klo  ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
                     Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo  ) ) : dz_no_terrain;
-                    Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr,
+                    Real tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr,
                                                           cons_arr, velx_arr, vely_arr,
                                                           umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr,
                                                           dest_arr);
@@ -314,10 +324,10 @@ ABLMost::compute_most_bcs (const int& lev,
                     // TODO: make sure not to double-count surface heat flux if using a LSM
                     int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
                     if (is_land && lsm_flux_arr && vbx.contains(i,j,k)) {
-                        lsm_flux_arr(i,j,klo) = Tflux;
+                        lsm_flux_arr(i,j,klo) = tflux;
                     }
                     else if ((k == klo-1) && vbx.contains(i,j,k) && exp_most) {
-                        hfx_arr(i,j,klo-1) = Tflux;
+                        hfx_arr(i,j,klo) = tflux;
                     }
                 });
 
@@ -328,11 +338,14 @@ ABLMost::compute_most_bcs (const int& lev,
                     {
                         Real dz  = (zphys_arr) ? ( zphys_arr(i,j,klo  ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
                         Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo  ) ) : dz_no_terrain;
-                        Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr,
+                        Real qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr,
                                                               cons_arr, velx_arr, vely_arr,
                                                               umm_arr, qm_arr, u_star_arr, q_star_arr, t_surf_arr,
                                                               dest_arr);
-                        amrex::ignore_unused(Qflux);
+
+                        if ((k == klo-1) && vbx.contains(i,j,k) && exp_most) {
+                            qfx_arr(i,j,klo) = qflux;
+                        }
                     });
                 }
 
@@ -350,10 +363,9 @@ ABLMost::compute_most_bcs (const int& lev,
                                                             umm_arr, um_arr, u_star_arr,
                                                             dest_arr);
                     if ((k == klo-1) && vbxx.contains(i,j,k) && exp_most) {
-                        t13_arr(i,j,klo) = -stressx;
-                        if (t31_arr) t31_arr(i,j,klo) = -stressx;
+                        t13_arr(i,j,klo) = stressx;
+                        if (t31_arr) t31_arr(i,j,klo) = stressx;
                     }
-                    amrex::ignore_unused(stressx);
                 });
 
             } else if (var_idx == Vars::yvel) {
@@ -370,10 +382,9 @@ ABLMost::compute_most_bcs (const int& lev,
                                                             umm_arr, vm_arr, u_star_arr,
                                                             dest_arr);
                     if ((k == klo-1) && vbxy.contains(i,j,k) && exp_most) {
-                        t23_arr(i,j,klo) = -stressy;
-                        if (t32_arr) t32_arr(i,j,klo) = -stressy;
+                        t23_arr(i,j,klo) = stressy;
+                        if (t32_arr) t32_arr(i,j,klo) = stressy;
                     }
-                    amrex::ignore_unused(stressy);
                 });
             }
         } // var_idx

Source/BoundaryConditions/ERF_FillPatch.cpp

@@ -420,6 +420,7 @@ ERF::FillIntermediatePatch (int lev, Real time,
                                 Tau13_lev[lev].get(), Tau31_lev[lev].get(),
                                 Tau23_lev[lev].get(), Tau32_lev[lev].get(),
                                 SFS_hfx3_lev[lev].get(),
+                                SFS_qfx3_lev[lev].get(),
                                 z_phys_nd[lev].get());
     }