From c8d0c9cc6ea2651d101c970a0c3a1a2d117646dd Mon Sep 17 00:00:00 2001
From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com>
Date: Sat, 23 Mar 2024 22:48:14 -0700
Subject: [PATCH] use scale height in nh pressure

---
 src/cache/precomputed_quantities.jl        |  1 -
 src/prognostic_equations/edmfx_closures.jl | 40 +++-------------------
 toml/prognostic_edmfx_bomex_box.toml       |  3 ++
 toml/prognostic_edmfx_box.toml             |  3 ++
 toml/prognostic_edmfx_dycoms_rf01_box.toml |  3 ++
 5 files changed, 14 insertions(+), 36 deletions(-)

diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 43cc2dcf50..34b6f88b40 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -74,7 +74,6 @@ function precomputed_quantities(Y, atmos)
             ᶜK_u = similar(Y.c, FT),
             ᶜK_h = similar(Y.c, FT),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
-            ᶜupdraft_top = similar(Fields.level(Y.c, 1), FT),
             ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}),
             ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
             ᶜKʲs = similar(Y.c, NTuple{n, FT}),
diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl
index e603a78f14..55df94579d 100644
--- a/src/prognostic_equations/edmfx_closures.jl
+++ b/src/prognostic_equations/edmfx_closures.jl
@@ -71,9 +71,8 @@ function ᶠupdraft_nh_pressure(
     ᶠbuoyʲ,
     ᶠu3ʲ,
     ᶠu3⁰,
-    updraft_top,
+    plume_height,
 )
-
     if !nh_pressure_flag
         return zero(ᶠu3ʲ)
     else
@@ -86,12 +85,11 @@ function ᶠupdraft_nh_pressure(
         # Independence of aspect ratio hardcoded: α₂_asp_ratio² = FT(0)
 
         H_up_min = CAP.min_updraft_top(turbconv_params)
-        plume_scale_height = max(updraft_top, H_up_min)
 
         # We also used to have advection term here: α_a * w_up * div_w_up
         return α_b * ᶠbuoyʲ +
                α_d * (ᶠu3ʲ - ᶠu3⁰) * CC.Geometry._norm(ᶠu3ʲ - ᶠu3⁰, ᶠlg) /
-               plume_scale_height
+               max(plume_height, H_up_min)
     end
 end
 
@@ -108,40 +106,12 @@ function edmfx_nh_pressure_tendency!(
     n = n_mass_flux_subdomains(turbconv_model)
     (; params) = p
     (; ᶠgradᵥ_ᶜΦ) = p.core
-    (; ᶜρʲs, ᶠnh_pressure₃ʲs, ᶠu₃⁰, ᶜupdraft_top) = p.precomputed
-    FT = eltype(Y)
-    ᶜz = Fields.coordinate_field(Y.c).z
-    z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
-    z_sfc_data = Fields.field_values(z_sfc[colidx])
+    (; ᶜρʲs, ᶠnh_pressure₃ʲs, ᶠu₃⁰) = p.precomputed
     ᶠlg = Fields.local_geometry_field(Y.f)
 
-    turbconv_params = CAP.turbconv_params(params)
-    a_min = CAP.min_area(turbconv_params)
+    scale_height = CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
 
     for j in 1:n
-
-        # look for updraft top
-        ᶜupdraft_top_data = Fields.field_values(ᶜupdraft_top[colidx])
-        @. ᶜupdraft_top_data = FT(0)
-        for level in 1:Spaces.nlevels(axes(ᶜz))
-            ρaʲ_lev = Fields.field_values(
-                Spaces.level(Y.c.sgsʲs.:($j).ρa[colidx], level),
-            )
-            ρʲ_lev =
-                Fields.field_values(Spaces.level(ᶜρʲs.:($j)[colidx], level))
-            ᶜz_lev = Fields.field_values(Spaces.level(ᶜz[colidx], level))
-            @. ᶜupdraft_top_data = ifelse(
-                draft_area(ρaʲ_lev, ρʲ_lev) > a_min,
-                ᶜz_lev,
-                ᶜupdraft_top_data,
-            )
-        end
-        @. ᶜupdraft_top_data = ᶜupdraft_top_data - z_sfc_data
-
-        # There's only one updraft_top per column, so it's
-        # safe to use at cell centers and cell faces, correct?:
-        ᶠupdraft_top = Fields.Field(ᶜupdraft_top_data, axes(ᶠu₃⁰[colidx]))
-
         @. ᶠnh_pressure₃ʲs.:($$j)[colidx] = ᶠupdraft_nh_pressure(
             params,
             p.atmos.edmfx_nh_pressure,
@@ -153,7 +123,7 @@ function edmfx_nh_pressure_tendency!(
             ),
             Y.f.sgsʲs.:($$j).u₃[colidx],
             ᶠu₃⁰[colidx],
-            ᶠupdraft_top,
+            scale_height,
         )
 
         @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -= ᶠnh_pressure₃ʲs.:($$j)[colidx]
diff --git a/toml/prognostic_edmfx_bomex_box.toml b/toml/prognostic_edmfx_bomex_box.toml
index 20ffe2d1af..5be8788fc4 100644
--- a/toml/prognostic_edmfx_bomex_box.toml
+++ b/toml/prognostic_edmfx_bomex_box.toml
@@ -36,3 +36,6 @@ value = 1
 
 [max_area_limiter_scale]
 value = 0
+
+[pressure_normalmode_drag_coeff]
+value = 40.0
diff --git a/toml/prognostic_edmfx_box.toml b/toml/prognostic_edmfx_box.toml
index 7341d979b5..654a543e81 100644
--- a/toml/prognostic_edmfx_box.toml
+++ b/toml/prognostic_edmfx_box.toml
@@ -33,3 +33,6 @@ value = 1
 
 [max_area_limiter_scale]
 value = 0
+
+[pressure_normalmode_drag_coeff]
+value = 40.0
diff --git a/toml/prognostic_edmfx_dycoms_rf01_box.toml b/toml/prognostic_edmfx_dycoms_rf01_box.toml
index 05b6c0cb91..023ace6560 100644
--- a/toml/prognostic_edmfx_dycoms_rf01_box.toml
+++ b/toml/prognostic_edmfx_dycoms_rf01_box.toml
@@ -36,3 +36,6 @@ value = 1
 
 [max_area_limiter_scale]
 value = 0.01
+
+[pressure_normalmode_drag_coeff]
+value = 40.0