fix dry default case.

erf-model · Oct 28, 2024 · 999d116 · 999d116
1 parent d6dedda
commit 999d116
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 3 deletions.
diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H
@@ -67,6 +67,31 @@ buoyancy_moist_anelastic (int& i,
     return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface);
 }
 
+AMREX_GPU_HOST
+AMREX_FORCE_INLINE
+amrex::Real
+buoyancy_dry_default (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(cell_data(i,j,k  ,Rho_comp), cell_data(i,j,k  ,RhoTheta_comp));
+    amrex::Real qplus  = (t_hi-t0_hi)/t0_hi;
+
+    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(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
+    amrex::Real qminus = (t_lo-t0_lo)/t0_lo;
+
+    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_GPU_HOST
 AMREX_FORCE_INLINE

diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp
@@ -150,9 +150,9 @@ 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,
-                                                                       p0_arr,r0_arr,cell_data);
+                        buoyancy_fab(i, j, k) = buoyancy_dry_default(i,j,k,
+                                                                     grav_gpu[2],rd_over_cp,
+                                                                     p0_arr,r0_arr,cell_data);
                     });
                 } // mfi
             } // buoyancy_type