update dry anelastic buoyancy to be consistent with moist (#1917)

erf-model · Oct 30, 2024 · 36e2d55 · 36e2d55
1 parent 4a4c136
commit 36e2d55
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 12 deletions.
diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H
@@ -11,23 +11,22 @@ buoyancy_dry_anelastic (int& i,
                         int& j,
                         int& k,
                         amrex::Real const& grav_gpu,
-                        amrex::Real const& rd_over_cp,
                         const amrex::Array4<const amrex::Real>& p0_arr,
                         const amrex::Array4<const amrex::Real>& r0_arr,
                         const amrex::Array4<const amrex::Real>& cell_data)
 {
-    amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
-    amrex::Real  t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
-    amrex::Real   t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp);
-    amrex::Real qplus  = (t_hi-t0_hi)/t0_hi;
+    // Note: this is the same term as the moist anelastic buoyancy when qv = qc = qt = 0
+    amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1);
+    amrex::Real theta_d_hi = cell_data(i,j,k  ,RhoTheta_comp)/r0_arr(i,j,k);
 
-    amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
-    amrex::Real  t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
-    amrex::Real   t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp);
-    amrex::Real qminus = (t_lo-t0_lo)/t0_lo;
+    amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi);
 
-    amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
-    return (-r0_q_avg * grav_gpu);
+    amrex::Real theta_d_0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)) / r0_arr(i,j,k-1);
+    amrex::Real theta_d_0_hi = getRhoThetagivenP(p0_arr(i,j,k  )) / r0_arr(i,j,k  );
+
+    amrex::Real theta_d_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi);
+
+    return (-grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface);
 }
 
 

diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp
@@ -73,7 +73,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
                     buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,
-                                                                   grav_gpu[2],rd_over_cp,
+                                                                   grav_gpu[2],
                                                                    p0_arr,r0_arr,cell_data);
                 });
             } else {