From 0993858f70dde264383ebc8ff3a3d062bd7b0674 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 30 Oct 2024 09:07:15 -0700 Subject: [PATCH] update dry anelastic buoyancy to be consistent with moist --- Source/SourceTerms/ERF_buoyancy_utils.H | 21 ++++++++++----------- Source/SourceTerms/ERF_make_buoyancy.cpp | 2 +- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H index 91f2a8ea0..c5d5b75be 100644 --- 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& p0_arr, const amrex::Array4& r0_arr, const amrex::Array4& 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 index 11a3743a9..528c82edf 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -73,7 +73,7 @@ void make_buoyancy (Vector& 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 {