Fix thetav for PBL schemes. (#2016)

erf-model · Dec 9, 2024 · 7849906 · 7849906
1 parent add69a0
commit 7849906
Show file tree

Hide file tree

Showing 13 changed files with 64 additions and 33 deletions.
diff --git a/Source/BoundaryConditions/ERF_ABLMost.H b/Source/BoundaryConditions/ERF_ABLMost.H
@@ -368,15 +368,19 @@ public:
     update_pblh (const int& lev,
                  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
                  amrex::MultiFab* z_phys_cc,
-                 const int RhoQv_comp, const int RhoQr_comp);
+                 const int RhoQv_comp,
+                 const int RhoQc_comp,
+                 const int RhoQr_comp);
 
     template <typename PBLHeightEstimator>
     void
     compute_pblh (const int& lev,
                   amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
                   amrex::MultiFab* z_phys_cc,
                   const PBLHeightEstimator& est,
-                  const int RhoQv_comp, const int RhoQr_comp);
+                  const int RhoQv_comp,
+                  const int RhoQc_comp,
+                  const int RhoQr_comp);
 
     void
     read_custom_roughness (const int& lev,

diff --git a/Source/BoundaryConditions/ERF_ABLMost.cpp b/Source/BoundaryConditions/ERF_ABLMost.cpp
@@ -561,11 +561,13 @@ void
 ABLMost::update_pblh (const int& lev,
                       Vector<Vector<MultiFab>>& vars,
                       MultiFab* z_phys_cc,
-                      int RhoQv_comp, int RhoQr_comp)
+                      int RhoQv_comp,
+                      int RhoQc_comp,
+                      int RhoQr_comp)
 {
     if (pblh_type == PBLHeightCalcType::MYNN25) {
         MYNNPBLH estimator;
-        compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQr_comp);
+        compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQc_comp, RhoQr_comp);
     } else if (pblh_type == PBLHeightCalcType::YSU) {
         amrex::Error("YSU PBLH calc not implemented yet");
     }
@@ -577,11 +579,13 @@ ABLMost::compute_pblh (const int& lev,
                        Vector<Vector<MultiFab>>& vars,
                        MultiFab* z_phys_cc,
                        const PBLHeightEstimator& est,
-                       int RhoQv_comp, int RhoQr_comp)
+                       int RhoQv_comp,
+                       int RhoQc_comp,
+                       int RhoQr_comp)
 {
     est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
                      vars[lev][Vars::cons],m_lmask_lev[lev][0],
-                     RhoQv_comp, RhoQr_comp);
+                     RhoQv_comp, RhoQc_comp, RhoQr_comp);
 }
 
 void

diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H
@@ -107,17 +107,22 @@ struct SolverChoice {
         pp.query_enum_case_insensitive("moisture_model",moisture_type);
         if (moisture_type == MoistureType::SAM) {
             RhoQv_comp = RhoQ1_comp;
+            RhoQc_comp = RhoQ2_comp;
             RhoQr_comp = RhoQ4_comp;
         } else if (moisture_type == MoistureType::SAM_NoIce) {
             RhoQv_comp = RhoQ1_comp;
+            RhoQc_comp = RhoQ2_comp;
             RhoQr_comp = RhoQ4_comp;
         } else if (moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
             RhoQv_comp = RhoQ1_comp;
+            RhoQc_comp = RhoQ2_comp;
         } else if (moisture_type == MoistureType::Kessler) {
             RhoQv_comp = RhoQ1_comp;
+            RhoQc_comp = RhoQ2_comp;
             RhoQr_comp = RhoQ3_comp;
         } else if (moisture_type == MoistureType::Kessler_NoRain) {
             RhoQv_comp = RhoQ1_comp;
+            RhoQc_comp = RhoQ2_comp;
         }
 
         // TODO: should we set default for dry??
@@ -637,6 +642,7 @@ struct SolverChoice {
     bool do_precip {true};
     bool use_moist_background {false};
     int RhoQv_comp {-1};
+    int RhoQc_comp {-1};
 
     // This component will be model-dependent:
     // if a model with no rain, this will stay -1

diff --git a/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp b/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp
@@ -489,7 +489,8 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
         ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity,
                                  geom, turbChoice, most, use_moisture,
                                  level, bc_ptr, vert_only, z_phys_nd,
-                                 solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
+                                 solverChoice.RhoQv_comp, solverChoice.RhoQc_comp,
+                                 solverChoice.RhoQr_comp);
     } else if (turbChoice.pbl_type == PBLType::YSU) {
         ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity,
                               geom, turbChoice, most, use_moisture,

diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp
@@ -713,6 +713,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
         auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;
 
         const int rhoqv_comp = solverChoice.RhoQv_comp;
+        const int rhoqc_comp = solverChoice.RhoQc_comp;
         const int rhoqr_comp = solverChoice.RhoQr_comp;
 
         ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -728,7 +729,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
             cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
                                                                   mu_turb,cellSizeInv,domain,
                                                                   pbl_mynn_B1_l,tm_arr(i,j,0),
-                                                                  rhoqv_comp, rhoqr_comp,
+                                                                  rhoqv_comp, rhoqc_comp, rhoqr_comp,
                                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi,

diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp
@@ -824,6 +824,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
         auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;
 
         const int rhoqv_comp = solverChoice.RhoQv_comp;
+        const int rhoqc_comp = solverChoice.RhoQc_comp;
         const int rhoqr_comp = solverChoice.RhoQr_comp;
 
         ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -839,7 +840,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
             cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
                                                                   mu_turb,cellSizeInv,domain,
                                                                   pbl_mynn_B1_l,tm_arr(i,j,0),
-                                                                  rhoqv_comp, rhoqr_comp,
+                                                                  rhoqv_comp, rhoqc_comp, rhoqr_comp,
                                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi,

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -1061,7 +1061,9 @@ ERF::InitData_post ()
                 // we don't want to call update_fluxes multiple times because
                 // it will change u* and theta* from their previous values
                 m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get(),
-                                    solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
+                                    solverChoice.RhoQv_comp,
+                                    solverChoice.RhoQc_comp,
+                                    solverChoice.RhoQr_comp);
                 m_most->update_fluxes(lev, time);
             }
         }

diff --git a/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp b/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp
@@ -21,6 +21,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
                           bool /*vert_only*/,
                           const std::unique_ptr<MultiFab>& z_phys_nd,
                           const int RhoQv_comp,
+                          const int RhoQc_comp,
                           const int RhoQr_comp)
 {
     const bool use_terrain = (z_phys_nd != nullptr);
@@ -140,7 +141,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
                                           u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                           v_ext_dir_on_zlo, v_ext_dir_on_zhi,
                                           dthetadz, dudz, dvdz,
-                                          RhoQv_comp, RhoQr_comp, use_most);
+                                          RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most);
 
             // Spatially varying MOST
             Real theta0 = tm_arr(i,j,0);

diff --git a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp
@@ -204,13 +204,14 @@ ComputeDiffusivityYSU (const MultiFab& xvel,
                     constexpr Real prandtl_max = 4.0;
                     Real dthetadz, dudz, dvdz;
                     const int RhoQv_comp = -1;
+                    const int RhoQc_comp = -1;
                     const int RhoQr_comp = -1;
                     ComputeVerticalDerivativesPBL(i, j, k,
                                                   uvel, vvel, cell_data, izmin, izmax, 1.0/dz_terrain,
                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi,
-                                                  dthetadz, dudz, dvdz, RhoQv_comp, RhoQr_comp);
+                                                  dthetadz, dudz, dvdz, RhoQv_comp, RhoQc_comp, RhoQr_comp);
                     const Real shear_squared = dudz*dudz + dvdz*dvdz + 1.0e-9; // 1.0e-9 from WRF to avoid divide by zero
                     const Real theta = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
                     Real richardson = CONST_GRAV / theta * dthetadz / shear_squared;

diff --git a/Source/PBL/ERF_PBLHeight.H b/Source/PBL/ERF_PBLHeight.H
@@ -25,6 +25,7 @@ struct MYNNPBLH {
                       const amrex::MultiFab& cons,
                       const amrex::iMultiFab* lmask,
                       const int RhoQv_comp,
+                      const int RhoQc_comp,
                       const int RhoQr_comp) const
     {
 #if 0
@@ -35,7 +36,7 @@ struct MYNNPBLH {
         auto thetav_min = amrex::ReduceToPlane<amrex::ReduceOpMin,amrex::Real>(dir, bxlow, cons,
             [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real
             {
-                return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQr_comp);
+                return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQc_comp,RhoQr_comp);
             });
 #endif
 
@@ -88,7 +89,7 @@ struct MYNNPBLH {
                     int jj = amrex::max(amrex::min(j,jmax),jmin);
 
                     if (zphys_arr(ii,jj,k) < thetamin_height) {
-                        amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQr_comp);
+                        amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
                         if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
                     }
                 });
@@ -104,8 +105,8 @@ struct MYNNPBLH {
                         //
                         // Find PBL height based on thetav increase (best for CBLs)
                         //
-                        amrex::Real thv  = Thetav(i,j,k  ,cons_arr,RhoQv_comp,RhoQr_comp);
-                        amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQr_comp);
+                        amrex::Real thv  = Thetav(i,j,k  ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
+                        amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
 
                         int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
                         if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
@@ -164,7 +165,7 @@ struct MYNNPBLH {
 
                 ParallelFor(gtbxlow, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
-                    amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQr_comp);
+                    amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
                     if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
                 });
 
@@ -176,8 +177,8 @@ struct MYNNPBLH {
                         //
                         // Find PBL height based on thetav increase (best for CBLs)
                         //
-                        amrex::Real thv  = Thetav(i,j,k  ,cons_arr,RhoQv_comp,RhoQr_comp);
-                        amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQr_comp);
+                        amrex::Real thv  = Thetav(i,j,k  ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
+                        amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
 
                         int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
                         if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)

diff --git a/Source/PBL/ERF_PBLModels.H b/Source/PBL/ERF_PBLModels.H
@@ -33,7 +33,9 @@ ComputeDiffusivityMYNN25 (const amrex::MultiFab& xvel,
                           const amrex::BCRec* bc_ptr,
                           bool /*vert_only*/,
                           const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
-                          const int RhoQv_comp, const int RhoQr_comp);
+                          const int RhoQv_comp,
+                          const int RhoQc_comp,
+                          const int RhoQr_comp);
 
 /**
  * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
@@ -94,23 +96,24 @@ ComputeVerticalDerivativesPBL (int i, int j, int k,
                                amrex::Real& dudz,
                                amrex::Real& dvdz,
                                const int RhoQv_comp,
+                               const int RhoQc_comp,
                                const int RhoQr_comp,
                                const bool use_most=true)
 {
     if ( k==izmax && c_ext_dir_on_zhi ) {
-        dthetadz = (1.0/3.0)*(       -Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp)
-                              - 3.0 * Thetav(i,j,k  ,cell_data,RhoQv_comp,RhoQr_comp)
-                              + 4.0 * Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
+        dthetadz = (1.0/3.0)*(       -Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+                              - 3.0 * Thetav(i,j,k  ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+                              + 4.0 * Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
     } else if ( k==izmin && c_ext_dir_on_zlo ) {
-        dthetadz = (1.0/3.0)*(        Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp)
-                              + 3.0 * Thetav(i,j,k  ,cell_data,RhoQv_comp,RhoQr_comp)
-                              - 4.0 * Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
+        dthetadz = (1.0/3.0)*(        Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+                              + 3.0 * Thetav(i,j,k  ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+                              - 4.0 * Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
     } else if ( k==izmin && use_most ) {
-        dthetadz = (  Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp)
-                    - Thetav(i,j,k  ,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
+        dthetadz = (  Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+                    - Thetav(i,j,k  ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
     } else {
-        dthetadz = 0.5*(  Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp)
-                        - Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
+        dthetadz = 0.5*(  Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+                        - Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
     }
 
     if ( k==izmax && u_ext_dir_on_zhi ) {
@@ -165,6 +168,7 @@ ComputeQKESourceTerms (int i, int j, int k,
                        amrex::Real pbl_mynn_B1_l,
                        const amrex::Real theta_mean,
                        const int RhoQv_comp,
+                       const int RhoQc_comp,
                        const int RhoQr_comp,
                        bool c_ext_dir_on_zlo,
                        bool c_ext_dir_on_zhi,
@@ -195,7 +199,7 @@ ComputeQKESourceTerms (int i, int j, int k,
                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi,
                                   dthetadz, dudz, dvdz,
-                                  RhoQv_comp, RhoQr_comp, use_most);
+                                  RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most);
 
     // Note: Transport terms due to turbulence and pressure are included when
     // DiffusionSrcForState_* is called from ERF_slow_rhs_post.

diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp
@@ -110,7 +110,9 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/)
             //       Reassign the field ptrs for MAC avg computation.
             m_most->update_mac_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
             m_most->update_pblh(lev, vars_old, z_phys_cc[lev].get(),
-                                solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
+                                solverChoice.RhoQv_comp,
+                                solverChoice.RhoQc_comp,
+                                solverChoice.RhoQr_comp);
             m_most->update_fluxes(lev, time);
         }
     }

diff --git a/Source/Utils/ERF_Thetav.H b/Source/Utils/ERF_Thetav.H
@@ -10,15 +10,18 @@ amrex::Real
 Thetav (int i, int j, int k,
         const amrex::Array4<const amrex::Real>& cell_data,
         const int RhoQv_comp,
+        const int RhoQc_comp,
         const int RhoQr_comp)
 {
     amrex::Real thetav = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
 
     if (RhoQr_comp > 0) {
         thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQv_comp) / cell_data(i,j,k,Rho_comp)
+                              - cell_data(i,j,k,RhoQc_comp) / cell_data(i,j,k,Rho_comp)
                               - cell_data(i,j,k,RhoQr_comp) / cell_data(i,j,k,Rho_comp));
     } else if (RhoQv_comp > 0) {
-        thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQ1_comp) / cell_data(i,j,k,Rho_comp));
+        thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQ1_comp) / cell_data(i,j,k,Rho_comp)
+                              - cell_data(i,j,k,RhoQc_comp) / cell_data(i,j,k,Rho_comp));
     }
 
     return thetav;