diff --git a/integration_tests/ARM_SGP.jl b/integration_tests/ARM_SGP.jl
index d759c57fbc..0996673563 100644
--- a/integration_tests/ARM_SGP.jl
+++ b/integration_tests/ARM_SGP.jl
@@ -11,18 +11,18 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 6.4559638590248924e-01
-best_mse["updraft_area"] = 2.9732339500518958e+02
-best_mse["updraft_w"] = 1.0894427841756977e+00
-best_mse["updraft_qt"] = 5.2188511467832939e+01
-best_mse["updraft_thetal"] = 6.9082849512451588e+01
-best_mse["u_mean"] = 8.7994360629094146e+01
-best_mse["tke_mean"] = 3.0426947022991828e+00
-best_mse["temperature_mean"] = 2.9081377572297745e-04
-best_mse["ql_mean"] = 0.0000000000000000e+00
-best_mse["thetal_mean"] = 2.8896514406812518e-04
-best_mse["Hvar_mean"] = 8.9154045834618110e+00
-best_mse["QTvar_mean"] = 4.3180594810766799e+01
+best_mse["qt_mean"] = 6.4706501357843726e-01
+best_mse["updraft_area"] = 2.9291358318960744e+02
+best_mse["updraft_w"] = 9.9252910692095353e-01
+best_mse["updraft_qt"] = 5.2976988305713625e+01
+best_mse["updraft_thetal"] = 6.9069465033893721e+01
+best_mse["u_mean"] = 8.7994360629094245e+01
+best_mse["tke_mean"] = 2.9571812916570166e+00
+best_mse["temperature_mean"] = 2.9042903519224878e-04
+best_mse["ql_mean"] = 4.2347733720323376e-02
+best_mse["thetal_mean"] = 2.8846582418377522e-04
+best_mse["Hvar_mean"] = 7.0891467145334675e+00
+best_mse["QTvar_mean"] = 3.0196125753072817e+01
 
 @testset "ARM_SGP" begin
     println("Running ARM_SGP...")
diff --git a/integration_tests/Bomex.jl b/integration_tests/Bomex.jl
index 4915efa8bf..9d33dd6fd6 100644
--- a/integration_tests/Bomex.jl
+++ b/integration_tests/Bomex.jl
@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 8.9644597356208541e-02
-best_mse["updraft_area"] = 6.9327817412026047e+02
-best_mse["updraft_w"] = 8.9057089858539797e+01
-best_mse["updraft_qt"] = 5.9179688063528078e+00
-best_mse["updraft_thetal"] = 2.2933264810554988e+01
-best_mse["v_mean"] = 1.2348440297549568e+02
-best_mse["u_mean"] = 5.3484778922733959e+01
-best_mse["tke_mean"] = 3.2847767205446033e+01
-best_mse["temperature_mean"] = 3.4874387520736297e-05
-best_mse["ql_mean"] = 2.1344864368686135e+01
-best_mse["thetal_mean"] = 3.5327188090292184e-05
-best_mse["Hvar_mean"] = 1.7623971358198961e+02
-best_mse["QTvar_mean"] = 4.0583986764371645e+01
+best_mse["qt_mean"] = 8.9967689186130878e-02
+best_mse["updraft_area"] = 6.9301645004274201e+02
+best_mse["updraft_w"] = 8.9039019651975735e+01
+best_mse["updraft_qt"] = 5.9168523443194641e+00
+best_mse["updraft_thetal"] = 2.2933370324850653e+01
+best_mse["v_mean"] = 1.2348563017334499e+02
+best_mse["u_mean"] = 5.3484900051947236e+01
+best_mse["tke_mean"] = 3.2906542599846873e+01
+best_mse["temperature_mean"] = 3.4859944892702567e-05
+best_mse["ql_mean"] = 2.2542367939656252e+01
+best_mse["thetal_mean"] = 3.5354307559496704e-05
+best_mse["Hvar_mean"] = 1.8023023675915451e+02
+best_mse["QTvar_mean"] = 4.0992883352799481e+01
 
 @testset "Bomex" begin
     println("Running Bomex...")
diff --git a/integration_tests/DYCOMS_RF01.jl b/integration_tests/DYCOMS_RF01.jl
index da6bb4836e..fd8f91568b 100644
--- a/integration_tests/DYCOMS_RF01.jl
+++ b/integration_tests/DYCOMS_RF01.jl
@@ -11,20 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 6.3055740870523944e-02
-best_mse["ql_mean"] = 1.0898102515993909e+01
-best_mse["updraft_area"] = 2.1643245563605532e+02
-best_mse["updraft_w"] = 3.7742918972750585e+00
-best_mse["updraft_qt"] = 1.4833369346056364e+00
-best_mse["updraft_thetal"] = 1.2735585448251884e+01
-best_mse["v_mean"] = 4.0028336968928720e+01
-best_mse["u_mean"] = 3.5748548802210514e+01
-best_mse["tke_mean"] = 1.3334607271849251e+01
-best_mse["temperature_mean"] = 1.3282908741777973e-05
-best_mse["thetal_mean"] = 1.7266062255939275e-05
-best_mse["Hvar_mean"] = 8.5678621532669626e+04
-best_mse["QTvar_mean"] = 6.5174106643485084e+03
-
+best_mse["qt_mean"] = 6.3055738374403369e-02
+best_mse["ql_mean"] = 1.0898094895741076e+01
+best_mse["updraft_area"] = 2.1643245239719272e+02
+best_mse["updraft_w"] = 3.7742979507011780e+00
+best_mse["updraft_qt"] = 1.4833369559877259e+00
+best_mse["updraft_thetal"] = 1.2735585449339219e+01
+best_mse["v_mean"] = 4.0028336999479734e+01
+best_mse["u_mean"] = 3.5748548829194654e+01
+best_mse["tke_mean"] = 1.3334591135952765e+01
+best_mse["temperature_mean"] = 1.3282907861220565e-05
+best_mse["thetal_mean"] = 1.7266061266968340e-05
+best_mse["Hvar_mean"] = 8.5678573512520568e+04
+best_mse["QTvar_mean"] = 6.5174070171062540e+03
 
 @testset "DYCOMS_RF01" begin
     println("Running DYCOMS_RF01...")
diff --git a/integration_tests/Nieuwstadt.jl b/integration_tests/Nieuwstadt.jl
index 0bd8832b76..9227e1e0ab 100644
--- a/integration_tests/Nieuwstadt.jl
+++ b/integration_tests/Nieuwstadt.jl
@@ -11,14 +11,14 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["updraft_area"] = 6.7123367331397571e+02
-best_mse["updraft_w"] = 4.3115995379486961e+01
-best_mse["updraft_thetal"] = 2.9769609162505365e+01
-best_mse["u_mean"] = 7.6376436377614436e+02
-best_mse["tke_mean"] = 7.2281909958288438e+01
-best_mse["temperature_mean"] = 1.5773181513042552e-05
-best_mse["thetal_mean"] = 1.5665330726140129e-05
-best_mse["Hvar_mean"] = 2.0155355014688541e+02
+best_mse["updraft_area"] = 6.6806536316930988e+02
+best_mse["updraft_w"] = 4.3096286654742102e+01
+best_mse["updraft_thetal"] = 2.9769609085248476e+01
+best_mse["u_mean"] = 7.6375379428605470e+02
+best_mse["tke_mean"] = 7.2286201690254956e+01
+best_mse["temperature_mean"] = 1.5776736605110006e-05
+best_mse["thetal_mean"] = 1.5668943844822401e-05
+best_mse["Hvar_mean"] = 2.0155012796131999e+02
 
 @testset "Nieuwstadt" begin
     println("Running Nieuwstadt...")
diff --git a/integration_tests/Rico.jl b/integration_tests/Rico.jl
index 44ae53b235..e52abcd3f2 100644
--- a/integration_tests/Rico.jl
+++ b/integration_tests/Rico.jl
@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 4.9523533702786315e-01
-best_mse["updraft_area"] = 1.7419532961183722e+03
-best_mse["updraft_w"] = 5.3893010779565452e+02
-best_mse["updraft_qt"] = 2.4530572533274270e+01
-best_mse["updraft_thetal"] = 6.5509509302767754e+01
-best_mse["v_mean"] = 1.0615326468150907e+02
-best_mse["u_mean"] = 1.1393191702720974e+02
-best_mse["tke_mean"] = 8.4090034387266689e+02
-best_mse["temperature_mean"] = 1.4274314152822176e-04
-best_mse["ql_mean"] = 2.1441653002617767e+02
-best_mse["thetal_mean"] = 1.3609606969128344e-04
-best_mse["Hvar_mean"] = 3.7305257502890240e+03
-best_mse["QTvar_mean"] = 1.2241261056877631e+03
+best_mse["qt_mean"] = 4.9368437132590526e-01
+best_mse["updraft_area"] = 1.7440556056666451e+03
+best_mse["updraft_w"] = 5.3170885681379491e+02
+best_mse["updraft_qt"] = 2.4531006750397474e+01
+best_mse["updraft_thetal"] = 6.5509500570416122e+01
+best_mse["v_mean"] = 1.0615615061647286e+02
+best_mse["u_mean"] = 1.1393423823185431e+02
+best_mse["tke_mean"] = 8.3936753844031784e+02
+best_mse["temperature_mean"] = 1.4250540709147947e-04
+best_mse["ql_mean"] = 2.1055525812660557e+02
+best_mse["thetal_mean"] = 1.3576737175490100e-04
+best_mse["Hvar_mean"] = 4.2976683292044045e+03
+best_mse["QTvar_mean"] = 1.3538355966380366e+03
 
 @testset "Rico" begin
     println("Running Rico...")
diff --git a/integration_tests/Soares.jl b/integration_tests/Soares.jl
index 07ac51b047..469d0c2069 100644
--- a/integration_tests/Soares.jl
+++ b/integration_tests/Soares.jl
@@ -11,17 +11,16 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 3.2147414017234566e-01
-best_mse["updraft_area"] = 7.5523687950242140e+02
-best_mse["updraft_w"] = 3.5564080086151499e+01
-best_mse["updraft_qt"] = 1.0702812362369221e+01
-best_mse["updraft_thetal"] = 2.1534108878970788e+01
-best_mse["u_mean"] = 3.9385484811475408e+03
-best_mse["tke_mean"] = 5.2022934994314660e+01
-best_mse["temperature_mean"] = 1.9333588833166558e-05
-best_mse["thetal_mean"] = 1.9486304719984318e-05
-best_mse["Hvar_mean"] = 2.6039266998435176e+02
-
+best_mse["qt_mean"] = 3.2147414156385279e-01
+best_mse["updraft_area"] = 7.5523663990367277e+02
+best_mse["updraft_w"] = 3.5564218002606921e+01
+best_mse["updraft_qt"] = 1.0702812274245508e+01
+best_mse["updraft_thetal"] = 2.1534108876550594e+01
+best_mse["u_mean"] = 3.9385483882997114e+03
+best_mse["tke_mean"] = 5.2022941689385171e+01
+best_mse["temperature_mean"] = 1.9333589218733757e-05
+best_mse["thetal_mean"] = 1.9486305081638964e-05
+best_mse["Hvar_mean"] = 2.6039260879907039e+02
 
 @testset "Soares" begin
     println("Running Soares...")
diff --git a/integration_tests/TRMM_LBA.jl b/integration_tests/TRMM_LBA.jl
index 692fa0eee8..58f6012878 100644
--- a/integration_tests/TRMM_LBA.jl
+++ b/integration_tests/TRMM_LBA.jl
@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 3.0609461714234958e+00
-best_mse["updraft_area"] = 3.5992222800156262e+04
-best_mse["updraft_w"] = 1.0100451430716917e+03
-best_mse["updraft_qt"] = 2.7507151910663865e+01
-best_mse["updraft_thetal"] = 1.1074070908480367e+02
-best_mse["v_mean"] = 2.9438922611828133e+02
-best_mse["u_mean"] = 1.6906200771006561e+03
-best_mse["tke_mean"] = 2.5854592261914117e+03
-best_mse["temperature_mean"] = 9.2883217582418853e-04
-best_mse["ql_mean"] = 2.1864394985446243e+03
-best_mse["thetal_mean"] = 6.4455936975381660e-04
-best_mse["Hvar_mean"] = 1.2394474700136290e+04
-best_mse["QTvar_mean"] = 9.2129571715440597e+03
+best_mse["qt_mean"] = 3.0609725711974320e+00
+best_mse["updraft_area"] = 3.5994430618299564e+04
+best_mse["updraft_w"] = 1.0100729940018912e+03
+best_mse["updraft_qt"] = 2.7508554326878087e+01
+best_mse["updraft_thetal"] = 1.1074070149524526e+02
+best_mse["v_mean"] = 2.9438910525490803e+02
+best_mse["u_mean"] = 1.6906163382368013e+03
+best_mse["tke_mean"] = 2.5824652671530885e+03
+best_mse["temperature_mean"] = 9.2886791818677325e-04
+best_mse["ql_mean"] = 2.1883827035908516e+03
+best_mse["thetal_mean"] = 6.4459500424083574e-04
+best_mse["Hvar_mean"] = 1.2398403742209141e+04
+best_mse["QTvar_mean"] = 9.2234342640415034e+03
 
 @testset "TRMM_LBA" begin
     println("Running TRMM_LBA...")
diff --git a/src/Operators.jl b/src/Operators.jl
index f8e2aa1d7c..fc4f6be316 100644
--- a/src/Operators.jl
+++ b/src/Operators.jl
@@ -16,6 +16,9 @@ struct Extrapolate end
 function ∇f2c(f, grid::Grid, k::Int)
     return (f[k] - f[k - 1]) * grid.dzi
 end
+function ∇f2c(f, grid::Grid, k::Int, i::Int)
+    return (f[i, k] - f[i, k - 1]) * grid.dzi
+end
 function ∇f2c(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
     if is_surface_face(grid, k - 1)
         return ∇f2c(f_dual, grid, BottomBCTag(), bottom)
diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl
index 834879f476..c7426749f3 100644
--- a/src/TurbulenceConvection.jl
+++ b/src/TurbulenceConvection.jl
@@ -68,5 +68,6 @@ include("Radiation.jl")
 include("forcing_functions.jl")
 include("Surface.jl")
 include("surface_functions.jl")
+include("closures/perturbation_pressure.jl")
 
 end
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index a77db67af7..58bd0159af 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -24,8 +24,6 @@ function initialize_io(self::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats)
     add_profile(Stats, "nh_pressure_drag")
     add_profile(Stats, "nh_pressure_b")
     add_profile(Stats, "asp_ratio")
-    add_profile(Stats, "b_coeff")
-
     add_profile(Stats, "horiz_K_eddy")
     add_profile(Stats, "sorting_function")
     add_profile(Stats, "b_mix")
@@ -99,8 +97,6 @@ function io(self::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats, TS::TimeStepping)
     mean_nh_pressure_drag = center_field(grid)
     mean_nh_pressure_b = center_field(grid)
     mean_asp_ratio = center_field(grid)
-    mean_b_coeff = center_field(grid)
-
     mean_detr_sc = center_field(grid)
     massflux = face_field(grid)
     mf_h = face_field(grid)
@@ -135,8 +131,6 @@ function io(self::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats, TS::TimeStepping)
                     self.UpdVar.Area.values[i, k] * self.nh_pressure_drag[i, k] / self.UpdVar.Area.bulkvalues[k]
                 mean_asp_ratio[k] +=
                     self.UpdVar.Area.values[i, k] * self.asp_ratio[i, k] / self.UpdVar.Area.bulkvalues[k]
-                mean_b_coeff[k] += self.UpdVar.Area.values[i, k] * self.b_coeff[i, k] / self.UpdVar.Area.bulkvalues[k]
-
                 mean_frac_turb_entr[k] +=
                     self.UpdVar.Area.values[i, k] * self.frac_turb_entr[i, k] / self.UpdVar.Area.bulkvalues[k]
                 mean_horiz_K_eddy[k] +=
@@ -159,8 +153,6 @@ function io(self::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats, TS::TimeStepping)
     write_profile(Stats, "nh_pressure_drag", mean_nh_pressure_drag[cinterior])
     write_profile(Stats, "nh_pressure_b", mean_nh_pressure_b[cinterior])
     write_profile(Stats, "asp_ratio", mean_asp_ratio[cinterior])
-    write_profile(Stats, "b_coeff", mean_b_coeff[cinterior])
-
     write_profile(Stats, "massflux", massflux[cinterior])
     write_profile(Stats, "massflux_h", mf_h[cinterior])
     write_profile(Stats, "massflux_qt", mf_qt[cinterior])
@@ -931,9 +923,6 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
     input.c_ed_mf = self.entrainment_ed_mf_sigma
     input.chi_upd = self.updraft_mixing_frac
     input.tke_coef = self.entrainment_smin_tke_coeff
-    ret_b = pressure_buoy_struct()
-    ret_w = pressure_drag_struct()
-    input_p = pressure_in_struct()
 
     # TODO: merge this with grid.cinterior
     kc_surf = kc_surface(grid)
@@ -948,10 +937,8 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
 
     @inbounds for k in real_center_indicies(grid)
         @inbounds for i in xrange(self.n_updrafts)
-            input_p.updraft_top = self.UpdVar.updraft_top[i]
             alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])), 1)
             avals = self.UpdVar.Area.values[i, cinterior]
-            input_p.a_med = Statistics.median(avals[1:alen])
             input.zi = self.UpdVar.cloud_base[i]
             # entrainment
             input.buoy_ed_flux = self.EnvVar.TKE.buoy[k]
@@ -967,7 +954,6 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
                 input.w_env = interpf2c(self.EnvVar.W.values, grid, k)
                 input.ql_up = self.UpdVar.QL.values[i, k]
                 input.ql_env = self.EnvVar.QL.values[k]
-                input.nh_pressure = self.nh_pressure[i, k]
                 input.RH_upd = self.UpdVar.RH.values[i, k]
                 input.RH_env = self.EnvVar.RH.values[k]
 
@@ -1015,55 +1001,32 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
             end
 
             # pressure
-            input_p.a_kfull = interpc2f(self.UpdVar.Area.values, grid, k, i)
-            if input_p.a_kfull >= self.minimum_area
-                input_p.dzi = grid.dzi
-
-                input_p.b_kfull = interpc2f(self.UpdVar.B.values, grid, k, i)
-                input_p.rho0_kfull = ref_state.rho0[k]
-                input_p.alpha1 = self.pressure_normalmode_buoy_coeff1
-                input_p.alpha2 = self.pressure_normalmode_buoy_coeff2
-                input_p.beta1 = self.pressure_normalmode_adv_coeff
-                input_p.beta2 = self.pressure_normalmode_drag_coeff
-                input_p.w_kfull = self.UpdVar.W.values[i, k]
-                input_p.w_kmfull = self.UpdVar.W.values[i, k - 1]
-                input_p.w_kenv = self.EnvVar.W.values[k]
-
-                if self.asp_label == "z_dependent"
-                    input_p.asp_ratio = input_p.updraft_top / 2.0 / sqrt(input_p.a_kfull) / input_p.rd
-                elseif self.asp_label == "median"
-                    input_p.asp_ratio = input_p.updraft_top / 2.0 / sqrt(input_p.a_med) / input_p.rd
-                elseif self.asp_label == "const"
-                    input_p.asp_ratio = 1.0
-                end
-
-                if input_p.a_kfull > 0.0
-                    ret_b = self.pressure_func_buoy(input_p)
-                    ret_w = self.pressure_func_drag(input_p)
-                    self.nh_pressure_b[i, k] = ret_b.nh_pressure_b
-                    self.nh_pressure_adv[i, k] = ret_w.nh_pressure_adv
-                    self.nh_pressure_drag[i, k] = ret_w.nh_pressure_drag
-
-                    self.b_coeff[i, k] = ret_b.b_coeff
-                    self.asp_ratio[i, k] = input_p.asp_ratio
-
-                else
-                    self.nh_pressure_b[i, k] = 0.0
-                    self.nh_pressure_adv[i, k] = 0.0
-                    self.nh_pressure_drag[i, k] = 0.0
-
-                    self.b_coeff[i, k] = 0.0
-                    self.asp_ratio[i, k] = 0.0
-                end
+            a_kfull = interpc2f(self.UpdVar.Area.values, grid, k, i)
+            if a_kfull > 0.0
+                b_kfull = interpc2f(self.UpdVar.B.values, grid, k, i)
+                ∇w_up = ∇f2c(self.UpdVar.W.values, grid, k, i)
+                asp_ratio = 1.0
+                self.nh_pressure_b[i, k], self.nh_pressure_adv[i, k], self.nh_pressure_drag[i, k] =
+                    perturbation_pressure(
+                        self.UpdVar.updraft_top[i],
+                        500.0,
+                        a_kfull,
+                        b_kfull,
+                        ref_state.rho0[k],
+                        self.pressure_normalmode_buoy_coeff1,
+                        self.pressure_normalmode_buoy_coeff2,
+                        self.pressure_normalmode_adv_coeff,
+                        self.pressure_normalmode_drag_coeff,
+                        self.UpdVar.W.values[i, k],
+                        ∇w_up,
+                        self.EnvVar.W.values[k],
+                        asp_ratio,
+                    )
             else
                 self.nh_pressure_b[i, k] = 0.0
                 self.nh_pressure_adv[i, k] = 0.0
                 self.nh_pressure_drag[i, k] = 0.0
-
-                self.b_coeff[i, k] = 0.0
-                self.asp_ratio[i, k] = 0.0
             end
-
             self.nh_pressure[i, k] = self.nh_pressure_b[i, k] + self.nh_pressure_adv[i, k] + self.nh_pressure_drag[i, k]
         end
     end
diff --git a/src/closures/perturbation_pressure.jl b/src/closures/perturbation_pressure.jl
new file mode 100644
index 0000000000..93392e950e
--- /dev/null
+++ b/src/closures/perturbation_pressure.jl
@@ -0,0 +1,54 @@
+"""
+    perturbation_pressure(
+                updraft_top,
+                a_up,
+                b_up,
+                ρ0_k,
+                α₁,
+                α₂,
+                β₁,
+                β₂,
+                w_up,
+                ∇w_up,
+                w_en,
+                asp_ratio
+                )
+Returns the value of perturbation pressure gradient
+for updraft i following He et al. (JAMES, 2020), given:
+ - `updraft_top,`, the heigh of the updraft in the previous timestep
+ - `a_up,`, updraft area
+ - `b_up,`, updraft buoyancy
+ - `ρ0_k,`, reference denbsity
+ - `α₁,`, pressure closure free parameter
+ - `α₂,`, pressure closure free parameter
+ - `β₁,`, pressure closure free parameter
+ - `β₂,`, pressure closure free parameter
+ - `w_up,`, updraft vertical velocity
+ - `∇w_up,`, updraft divergence of vertical velocity
+ - `w_en,`, environment vertical velocity
+ - `asp_ratio`, the aspect radio of the updraft
+"""
+function perturbation_pressure(
+    updraft_top,
+    min_updraft_top,
+    a_up,
+    b_up,
+    ρ0_k,
+    α₁,
+    α₂,
+    β₁,
+    β₂,
+    w_up,
+    ∇w_up,
+    w_en,
+    asp_ratio,
+)
+
+    nh_press_buoy = -α₁ / (1 + α₂ * asp_ratio^2) * ρ0_k * a_up * b_up
+
+    nh_pressure_adv = ρ0_k * a_up * β₁ * w_up * ∇w_up
+    # drag as w_dif and account for downdrafts
+    nh_pressure_drag = -1.0 * ρ0_k * a_up * β₂ * (w_up - w_en) * fabs(w_up - w_en) / fmax(updraft_top, min_updraft_top)
+
+    return nh_press_buoy, nh_pressure_adv, nh_pressure_drag
+end;
diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl
index f6e02d5090..27754264c4 100644
--- a/src/turbulence_functions.jl
+++ b/src/turbulence_functions.jl
@@ -142,39 +142,6 @@ function entr_detr_env_moisture_deficit_div(param_set::APS, entr_in::entr_struct
     return _ret
 end
 
-
-function pressure_normalmode_buoy(press_in::pressure_in_struct)
-    _ret = pressure_buoy_struct()
-
-    _ret.b_coeff = press_in.alpha1 / (1 + press_in.alpha2 * press_in.asp_ratio^2)
-    _ret.nh_pressure_b = -1.0 * press_in.rho0_kfull * press_in.a_kfull * press_in.b_kfull * _ret.b_coeff
-
-    return _ret
-end
-
-function pressure_normalmode_drag(press_in::pressure_in_struct)
-    _ret = pressure_drag_struct()
-
-    _ret.nh_pressure_adv =
-        press_in.rho0_kfull *
-        press_in.a_kfull *
-        press_in.beta1 *
-        press_in.w_kfull *
-        (press_in.w_kfull - press_in.w_kmfull) *
-        press_in.dzi
-
-    # drag as w_dif and account for downdrafts
-    _ret.nh_pressure_drag =
-        -1.0 *
-        press_in.rho0_kfull *
-        press_in.a_kfull *
-        press_in.beta2 *
-        (press_in.w_kfull - press_in.w_kenv) *
-        fabs(press_in.w_kfull - press_in.w_kenv) / fmax(press_in.updraft_top, 500.0)
-
-    return _ret
-end
-
 # convective velocity scale
 function get_wstar(bflux, zi)
     return cbrt(fmax(bflux * zi, 0.0))
diff --git a/src/types.jl b/src/types.jl
index 6e6f32480c..29ee0adf17 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,32 +1,3 @@
-Base.@kwdef mutable struct pressure_drag_struct
-    nh_pressure_adv::Float64 = 0
-    nh_pressure_drag::Float64 = 0
-end
-
-Base.@kwdef mutable struct pressure_buoy_struct
-    b_coeff::Float64 = 0
-    nh_pressure_b::Float64 = 0
-end
-
-Base.@kwdef mutable struct pressure_in_struct
-    updraft_top::Float64 = 0
-    asp_label::Float64 = 0
-    a_med::Float64 = 0
-    a_kfull::Float64 = 0
-    b_kfull::Float64 = 0
-    rho0_kfull::Float64 = 0
-    alpha1::Float64 = 0
-    alpha2::Float64 = 0
-    beta1::Float64 = 0
-    beta2::Float64 = 0
-    w_kfull::Float64 = 0
-    w_kmfull::Float64 = 0
-    w_kenv::Float64 = 0
-    dzi::Float64 = 0
-    drag_sign::Float64 = 0
-    asp_ratio::Float64 = 0
-end
-
 Base.@kwdef mutable struct entr_struct
     entr_sc::Float64 = 0
     detr_sc::Float64 = 0
@@ -753,9 +724,7 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
     n_updrafts::Int
     use_const_plume_spacing::Bool
     entr_detr_fp::Function
-    pressure_func_buoy::Function
     drag_sign::Int
-    pressure_func_drag::Function
     asp_label
     extrapolate_buoyancy::Bool
     surface_area::Float64
@@ -799,7 +768,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
     nh_pressure_adv::A2
     nh_pressure_drag::A2
     asp_ratio::A2
-    b_coeff::A2
     m::A2
     mixing_length::A1
     horiz_K_eddy::A2
@@ -864,8 +832,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
             error("Turbulence--EDMF_PrognosticTKE: Entrainment rate namelist option is not recognized")
         end
 
-        pressure_func_buoy = pressure_normalmode_buoy
-
         pressure_func_drag_str = parse_namelist(
             namelist,
             "turbulence",
@@ -875,12 +841,10 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
             valid_options = ["normalmode", "normalmode_signdf"],
         )
 
-        drag_sign = false
-        pressure_func_drag = if pressure_func_drag_str == "normalmode"
-            pressure_normalmode_drag
+        if pressure_func_drag_str == "normalmode"
+            drag_sign = false
         elseif pressure_func_drag_str == "normalmode_signdf"
             drag_sign = true
-            pressure_normalmode_drag
         end
 
         asp_label = parse_namelist(
@@ -994,7 +958,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
         nh_pressure_adv = center_field(Gr, n_updrafts)
         nh_pressure_drag = center_field(Gr, n_updrafts)
         asp_ratio = center_field(Gr, n_updrafts)
-        b_coeff = center_field(Gr, n_updrafts)
 
         # Mass flux
         m = face_field(Gr, n_updrafts)
@@ -1058,9 +1021,7 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
             n_updrafts,
             use_const_plume_spacing,
             entr_detr_fp,
-            pressure_func_buoy,
             drag_sign,
-            pressure_func_drag,
             asp_label,
             extrapolate_buoyancy,
             surface_area,
@@ -1104,7 +1065,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
             nh_pressure_adv,
             nh_pressure_drag,
             asp_ratio,
-            b_coeff,
             m,
             mixing_length,
             horiz_K_eddy,