move nh_pressure to closures

move nh_pressure to closures

update tables

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...")

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...")

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...")

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...")

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...")

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...")

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...")

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)

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

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