Merge pull request #3772 from CliMA/zs/edmf_nh_buoy

move virtual mass term to prognostic edmf momentum equation

Author: szy21

config/default_configs/default_config.yml

@@ -325,8 +325,9 @@ implicit_sgs_entr_detr:
   help: "Whether to treat the subgrid-scale entrainment and detrainment tendency implicitly [`false` (default), `true`]. Setting it to true only works if implicit_sgs_advection is set to true."
   value: false
 implicit_sgs_nh_pressure:
-  help: "Whether to treat the subgrid-scale nonhydrostatic pressure closure implicitly [`false` (default), `true`]. Setting it to true only works if implicit_sgs_advection is set to true.
-  This flag only controls whether the drag term in the pressure closure is treated implicitly. The buoyancy term is always treated explicitly."
+  help: "Whether to treat the subgrid-scale nonhydrostatic pressure closure implicitly. Setting it to true only works if 
+  `implicit_sgs_advection` is set to true. This flag only controls whether the drag term in the pressure closure is treated implicitly. 
+  The buoyancy term is always treated implicitly. [`false` (default), `true`]"
   value: false
 implicit_sgs_mass_flux:
   help: "Whether to treat the subgrid-scale mass flux tendency implicitly or explicitly in grid-mean equations. Currently updraft only with Jacobian terms 0. [`false` (default), `true`]. Setting it to true only works if both implicit_sgs_advection and implicit_diffusion are set to true."
@@ -338,7 +339,8 @@ edmfx_filter:
   help: "If set to true, it switches on the relaxation of negative velocity in EDMFX.  [`true`, `false` (default)]"
   value: false
 edmfx_nh_pressure:
-  help: "If set to true, it switches on EDMFX pressure drag closure.  [`true`, `false` (default)]"
+  help: "If set to true, it switches on EDMFX pressure drag closure.  For prognostic EDMF, this only controls the drag term 
+  in the pressure closure. The buoyancy term is always applied. [`true`, `false` (default)]"
   value: false
 edmfx_entr_model:
   help: "EDMFX entrainment closure.  [`nothing` (default), `PiGroups`, `Generalized`, `GeneralizedHarmonics`]"

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-233
+234
 
 # **README**
 #
@@ -21,34 +21,46 @@
 
 #=
 
+234
+- Move the virtual mass term in pressure closure to prognostic edmf momentum equation
+  (This only affects prognostic edmf)
+
 233
 - modify the derivative of p and \rho with respect to qt
 
 232
 - Use lazy broadcasting instead of temp scalar in implicit solver for kappa_m vars,
   which fixes a bug that the temp scalar is updated before it is reused.
+  (This only affects prognostic edmf)
 
 231
 - Add mass flux derivatives with respect to grid-mean u_3
+  (This only affects prognostic edmf)
 
 230
 - Add u_{3,m} (updraft) Jacobians to updraft MSE, rho*a, and q_tot prognostic equations. Move sgs ∂ᶠu₃ʲ derivatives to BlockLowerTriangularSolve.
+  (This only affects prognostic edmf)
 
 229
 - Remove derivatives with respect to grid mean rho in edmf implicit solver
+  (This only affects prognostic edmf)
 
 228
 - Only treat the drag term in edmf pressure closure implicitly
+  (This only affects prognostic edmf)
 
 227
 - Move nonhydrostatic pressure drag calculation to implicit precomputed quantities
+  (This only affects prognostic edmf)
 
 226
 - Add updraft rho*a and u_{3,m} jacobian terms
+  (This only affects prognostic edmf)
 
 225
 - Move nonhydrostatic pressure drag calculation to precomputed quantities and 
   remove one reproducibility job
+  (This only affects prognostic edmf)
 
 224
 - Machine precision differences due to https://github.com/CliMA/ClimaCore.jl/pull/2232

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -110,12 +110,9 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft!(
 
     n = n_mass_flux_subdomains(turbconv_model)
     thermo_params = CAP.thermodynamics_params(p.params)
-    turbconv_params = CAP.turbconv_params(p.params)
 
     (; ᶜΦ,) = p.core
-    (; ᶜspecific, ᶜp, ᶜh_tot, ᶜK) = p.precomputed
-    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
-    (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
+    (; ᶜp, ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
 
     for j in 1:n
         ᶜuʲ = ᶜuʲs.:($j)
@@ -468,15 +465,14 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
             dt,
         )
 
-        # nonhydrostatic pressure closure buoyancy term
-        if p.atmos.edmfx_model.nh_pressure isa Val{true}
-            @. ᶠnh_pressure₃_buoyʲs.:($$j) = ᶠupdraft_nh_pressure_buoyancy(
-                params,
-                ᶠbuoyancy(ᶠinterp(Y.c.ρ), ᶠinterp(ᶜρʲs.:($$j)), ᶠgradᵥ_ᶜΦ),
-            )
-        else
-            @. ᶠnh_pressure₃_buoyʲs.:($$j) = C3(0)
-        end
+        # The buoyancy term in the nonhydrostatic pressure closure is always applied
+        # for prognostic edmf. The tendency is combined with the buoyancy term in the 
+        # updraft momentum equation in `edmfx_sgs_vertical_advection_tendency!`. This
+        # term is still calculated here as it is used explicitly in the TKE equation.
+        @. ᶠnh_pressure₃_buoyʲs.:($$j) = ᶠupdraft_nh_pressure_buoyancy(
+            params,
+            ᶠbuoyancy(ᶠinterp(Y.c.ρ), ᶠinterp(ᶜρʲs.:($$j)), ᶠgradᵥ_ᶜΦ),
+        )
     end
 
     (; ᶜgradᵥ_θ_virt⁰, ᶜgradᵥ_q_tot⁰, ᶜgradᵥ_θ_liq_ice⁰) = p.precomputed

src/prognostic_equations/advection.jl

@@ -201,6 +201,8 @@ function edmfx_sgs_vertical_advection_tendency!(
     (; ᶠu³ʲs, ᶠKᵥʲs, ᶜρʲs) = p.precomputed
     (; ᶠgradᵥ_ᶜΦ) = p.core
 
+    turbconv_params = CAP.turbconv_params(params)
+    α_b = CAP.pressure_normalmode_buoy_coeff1(turbconv_params)
     ᶠz = Fields.coordinate_field(Y.f).z
     ᶜa_scalar = p.scratch.ᶜtemp_scalar
     ᶜu₃ʲ = p.scratch.ᶜtemp_C3
@@ -215,9 +217,10 @@ function edmfx_sgs_vertical_advection_tendency!(
         )
         # For the updraft u_3 equation, we assume the grid-mean to be hydrostatic
         # and calcuate the buoyancy term relative to the grid-mean density.
+        # We also include the buoyancy term in the nonhydrostatic pressure closure here.
         @. Yₜ.f.sgsʲs.:($$j).u₃ -=
-            (ᶠinterp(ᶜρʲs.:($$j) - Y.c.ρ) * ᶠgradᵥ_ᶜΦ) / ᶠinterp(ᶜρʲs.:($$j)) +
-            ᶠgradᵥ(ᶜKᵥʲ)
+            (1 - α_b) * (ᶠinterp(ᶜρʲs.:($$j) - Y.c.ρ) * ᶠgradᵥ_ᶜΦ) /
+            ᶠinterp(ᶜρʲs.:($$j)) + ᶠgradᵥ(ᶜKᵥʲ)
 
         # buoyancy term in mse equation
         @. Yₜ.c.sgsʲs.:($$j).mse +=

src/prognostic_equations/edmfx_closures.jl

@@ -88,21 +88,6 @@ function ᶠupdraft_nh_pressure_drag(params, ᶠlg, ᶠu3ʲ, ᶠu3⁰, scale_hei
            max(scale_height, H_up_min)
 end
 
-edmfx_nh_pressure_buoyancy_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
-function edmfx_nh_pressure_buoyancy_tendency!(
-    Yₜ,
-    Y,
-    p,
-    t,
-    turbconv_model::PrognosticEDMFX,
-)
-    (; ᶠnh_pressure₃_buoyʲs) = p.precomputed
-    n = n_mass_flux_subdomains(turbconv_model)
-    for j in 1:n
-        @. Yₜ.f.sgsʲs.:($$j).u₃ -= ᶠnh_pressure₃_buoyʲs.:($$j)
-    end
-end
-
 edmfx_nh_pressure_drag_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
 function edmfx_nh_pressure_drag_tendency!(
     Yₜ,

src/prognostic_equations/implicit/implicit_solver.jl

@@ -935,6 +935,9 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
                     Δcp_v * TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1))
                 ) / abs2(TD.cp_m(thermo_params, ᶜtsʲs.:(1)))
 
+            turbconv_params = CAP.turbconv_params(params)
+            α_b = CAP.pressure_normalmode_buoy_coeff1(turbconv_params)
+
             ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
                 matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)]
             @. ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
@@ -1085,7 +1088,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
                 matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)]
             @. ∂ᶠu₃ʲ_err_∂ᶜq_totʲ =
                 dtγ * DiagonalMatrixRow(
-                    ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
+                    (1 - α_b) * ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) /
+                    (ᶠinterp(ᶜρʲs.:(1)))^2,
                 ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
                     (ᶜρʲs.:(1))^2 / ᶜp * (
                         ᶜkappa_mʲ / (ᶜkappa_mʲ + 1) * ∂e_int_∂q_tot +
@@ -1100,7 +1104,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
                 matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)]
             @. ∂ᶠu₃ʲ_err_∂ᶜmseʲ =
                 dtγ * DiagonalMatrixRow(
-                    ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
+                    (1 - α_b) * ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) /
+                    (ᶠinterp(ᶜρʲs.:(1)))^2,
                 ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
                     ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
                 )

src/prognostic_equations/remaining_tendency.jl

@@ -118,7 +118,6 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     if p.atmos.sgs_mf_mode == Explicit()
         edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end
-    edmfx_nh_pressure_buoyancy_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     if p.atmos.sgs_nh_pressure_mode == Explicit()
         edmfx_nh_pressure_drag_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end