Fix mismatches. (#2102)

erf-model · Feb 5, 2025 · aaad4ba · aaad4ba
1 parent d91d64b
commit aaad4ba
Show file tree

Hide file tree

Showing 5 changed files with 14 additions and 10 deletions.
diff --git a/Source/Initialization/ERF_Init1D.cpp b/Source/Initialization/ERF_Init1D.cpp
@@ -218,7 +218,7 @@ ERF::erf_enforce_hse (int lev,
 
                 pres_arr(i,j,klo) = p_0 - hz * rho_arr(i,j,klo) * l_gravity;
                   pi_arr(i,j,klo) = getExnergivenP(pres_arr(i,j,klo), rdOcp);
-                  th_arr(i,j,klo) =getRhoThetagivenP(pres_arr(i,j,klo)) / rho_arr(i,j,klo);
+                  th_arr(i,j,klo) = getRhoThetagivenP(pres_arr(i,j,klo)) / rho_arr(i,j,klo);
 
                 //
                 // Set ghost cell with dz and rho at boundary

diff --git a/Source/Initialization/ERF_InitFromInputSounding.cpp b/Source/Initialization/ERF_InitFromInputSounding.cpp
@@ -275,7 +275,7 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx,
         r_hse_arr (i,j,k) = rho_k * (1.0 + qv_k);
         p_hse_arr (i,j,k) = getPgivenRTh(rhoTh_k, qv_k);
         pi_hse_arr(i,j,k) = getExnergivenRTh(rhoTh_k, l_rdOcp);
-        th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k)) / r_hse_arr(i,j,k);
+        th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k), qv_k) / rho_k;
 
         // TODO: we should be setting this to the number of ghost cells of base_state[lev]
         //       instead of hard-wiring it here!

diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp
@@ -1162,7 +1162,7 @@ init_base_state_from_metgrid (const bool use_moisture,
             r_hse_arr(i,j,k) *= (1.0 + Qv);
 
             pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
-            th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k)) / r_hse_arr(i,j,k);
+            th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k), Qv) / new_data(i,j,k,Rho_comp);
         });
 
         // FOEXTRAP hse arrays to fill ghost cells. FillBoundary is

diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp
@@ -622,7 +622,7 @@ init_base_state_from_wrfinput (const Box& domain,
             r_hse_arr(i,j,k)  = Rhse_Sum;
             p_hse_arr(i,j,k)  = Ptot;
             pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
-            th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k)) / r_hse_arr(i,j,k);
+            th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k), Qv) / cons_arr(ii,jj,kk,Rho_comp);
         });
     }
 }

diff --git a/Source/SourceTerms/ERF_BuoyancyUtils.H b/Source/SourceTerms/ERF_BuoyancyUtils.H
@@ -34,18 +34,22 @@ buoyancy_moist_anelastic (int& i, int& j, int& k,
                           const amrex::Array4<const amrex::Real>& th0_arr,
                           const amrex::Array4<const amrex::Real>& cell_data)
 {
-    amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1);
-    amrex::Real qv_lo      = cell_data(i,j,k-1,RhoQ1_comp)   /r0_arr(i,j,k-1);
-    amrex::Real qc_lo      = cell_data(i,j,k-1,RhoQ2_comp)   /r0_arr(i,j,k-1);
+    amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
+    amrex::Real qv_lo      = cell_data(i,j,k-1,RhoQ1_comp)   /cell_data(i,j,k-1,Rho_comp);
+    amrex::Real qc_lo      = cell_data(i,j,k-1,RhoQ2_comp)   /cell_data(i,j,k-1,Rho_comp);
     amrex::Real qt_lo      = qv_lo + qc_lo;
     amrex::Real theta_v_lo = theta_d_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo);
 
-    amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k);
-    amrex::Real qv_hi      = cell_data(i,j,k,RhoQ1_comp)   /r0_arr(i,j,k);
-    amrex::Real qc_hi      = cell_data(i,j,k,RhoQ2_comp)   /r0_arr(i,j,k);
+    amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/cell_data(i,j,k,Rho_comp);
+    amrex::Real qv_hi      = cell_data(i,j,k,RhoQ1_comp)   /cell_data(i,j,k,Rho_comp);
+    amrex::Real qc_hi      = cell_data(i,j,k,RhoQ2_comp)   /cell_data(i,j,k,Rho_comp);
     amrex::Real qt_hi      = qv_hi + qc_hi;
     amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi);
 
+    // NOTE: The following has three errors that depend upon storing a qv0 base state
+    //       1. theta_v0_wface needs to incorporate qv0
+    //       2. rho0_wface should be rhod0_wface (rho0 includes (1 + qv0) from balance w/ moisture)
+    //       3. another term with qv' = (qv - qv0) should be included
     amrex::Real  theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi);
     amrex::Real theta_v0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1));
     amrex::Real     rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));