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"] = 3.0998954346529128e-01
-best_mse["updraft_area"] = 1.9456011695010488e+03
-best_mse["updraft_w"] = 3.1996277691153108e+02
-best_mse["updraft_qt"] = 1.2645380500158543e+01
-best_mse["updraft_thetal"] = 2.7684653720274000e+01
-best_mse["u_mean"] = 8.7998547277817863e+01
-best_mse["tke_mean"] = 6.2536497518678300e+02
-best_mse["temperature_mean"] = 1.3316949170699788e-04
-best_mse["ql_mean"] = 2.5314704151953435e+02
-best_mse["thetal_mean"] = 1.3123118572261322e-04
-best_mse["Hvar_mean"] = 1.0770693519961310e+03
-best_mse["QTvar_mean"] = 1.9003025309800736e+02
+best_mse["qt_mean"] = 3.4761388491700224e-01
+best_mse["updraft_area"] = 1.9655632954763928e+03
+best_mse["updraft_w"] = 3.2377093752385053e+02
+best_mse["updraft_qt"] = 1.2452570404535185e+01
+best_mse["updraft_thetal"] = 2.7685659900036992e+01
+best_mse["u_mean"] = 8.7998547277817906e+01
+best_mse["tke_mean"] = 5.9501141635391093e+02
+best_mse["temperature_mean"] = 1.4210607813065063e-04
+best_mse["ql_mean"] = 2.1890556031513017e+02
+best_mse["thetal_mean"] = 1.3953340654098429e-04
+best_mse["Hvar_mean"] = 1.1148198912643536e+03
+best_mse["QTvar_mean"] = 3.0703058500395360e+02
 
 @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"] = 1.1113745694252045e-01
-best_mse["updraft_area"] = 7.2538205557139236e+02
-best_mse["updraft_w"] = 2.5943804836195852e+01
-best_mse["updraft_qt"] = 4.0075397466347100e+00
-best_mse["updraft_thetal"] = 2.1549266568293614e+01
-best_mse["v_mean"] = 6.5506192741355221e+01
-best_mse["u_mean"] = 5.3295866692805944e+01
-best_mse["tke_mean"] = 3.9596900219869937e+01
-best_mse["temperature_mean"] = 4.6114814306316929e-05
-best_mse["ql_mean"] = 6.7039710795893939e+00
-best_mse["thetal_mean"] = 4.6514684598969401e-05
-best_mse["Hvar_mean"] = 2.5387994922244559e+02
-best_mse["QTvar_mean"] = 4.5843881812484575e+01
+best_mse["qt_mean"] = 1.1068396401839847e-01
+best_mse["updraft_area"] = 7.2515940450436165e+02
+best_mse["updraft_w"] = 2.5856202193762641e+01
+best_mse["updraft_qt"] = 4.0061911576818243e+00
+best_mse["updraft_thetal"] = 2.1549386071253458e+01
+best_mse["v_mean"] = 6.5501952940732622e+01
+best_mse["u_mean"] = 5.3295928431425537e+01
+best_mse["tke_mean"] = 3.9582886707067345e+01
+best_mse["temperature_mean"] = 4.5842727523070819e-05
+best_mse["ql_mean"] = 6.5924727808279036e+00
+best_mse["thetal_mean"] = 4.6222694206747763e-05
+best_mse["Hvar_mean"] = 2.5900864672158934e+02
+best_mse["QTvar_mean"] = 4.6280994372607324e+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"] = 1.5987989356173508e-02
-best_mse["ql_mean"] = 5.9262542639081852e+00
-best_mse["updraft_area"] = 2.3887443993964089e+02
-best_mse["updraft_w"] = 4.3317357122037770e+00
-best_mse["updraft_qt"] = 1.1983341056077728e+00
-best_mse["updraft_thetal"] = 1.2739627176096105e+01
-best_mse["v_mean"] = 3.9736878160027324e+01
-best_mse["u_mean"] = 3.7037227811674803e+01
-best_mse["tke_mean"] = 1.4955313626928097e+01
-best_mse["temperature_mean"] = 1.6122516172306611e-05
-best_mse["thetal_mean"] = 1.7630641025752975e-05
-best_mse["Hvar_mean"] = 1.2452024821778659e+04
-best_mse["QTvar_mean"] = 7.6626261862502452e+02
-
+best_mse["qt_mean"] = 1.5987989288888820e-02
+best_mse["ql_mean"] = 5.9262500323990217e+00
+best_mse["updraft_area"] = 2.3887443639936907e+02
+best_mse["updraft_w"] = 4.3317426538007693e+00
+best_mse["updraft_qt"] = 1.1983341214387264e+00
+best_mse["updraft_thetal"] = 1.2739627177179765e+01
+best_mse["v_mean"] = 3.9736878190256498e+01
+best_mse["u_mean"] = 3.7037227839730683e+01
+best_mse["tke_mean"] = 1.4955295527017412e+01
+best_mse["temperature_mean"] = 1.6122517253407575e-05
+best_mse["thetal_mean"] = 1.7630642054888154e-05
+best_mse["Hvar_mean"] = 1.2452017842800817e+04
+best_mse["QTvar_mean"] = 7.6626218978684426e+02
 
 @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"] = 5.9123316810068854e+02
-best_mse["updraft_w"] = 2.6743480109271673e+01
-best_mse["updraft_thetal"] = 3.0475272195292995e+01
-best_mse["u_mean"] = 1.5245784141331271e+02
-best_mse["tke_mean"] = 7.3707701211117509e+01
-best_mse["temperature_mean"] = 1.2199756418560522e-05
-best_mse["thetal_mean"] = 1.2066937466844642e-05
-best_mse["Hvar_mean"] = 1.9140385882308689e+02
+best_mse["updraft_area"] = 5.9103576007452034e+02
+best_mse["updraft_w"] = 2.6741082264407293e+01
+best_mse["updraft_thetal"] = 3.0475272117346940e+01
+best_mse["u_mean"] = 1.5245873874565478e+02
+best_mse["tke_mean"] = 7.3707728452173882e+01
+best_mse["temperature_mean"] = 1.2199914281992930e-05
+best_mse["thetal_mean"] = 1.2067109520894632e-05
+best_mse["Hvar_mean"] = 1.9140690462510639e+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.2382033011491649e-01
-best_mse["updraft_area"] = 1.6696925513643459e+03
-best_mse["updraft_w"] = 2.8274002139289530e+02
-best_mse["updraft_qt"] = 1.7817151504225862e+01
-best_mse["updraft_thetal"] = 6.3573789244725759e+01
-best_mse["v_mean"] = 1.0633863265842136e+02
-best_mse["u_mean"] = 1.1434504431562398e+02
-best_mse["tke_mean"] = 8.6963627090840509e+02
-best_mse["temperature_mean"] = 1.2597722889848628e-04
-best_mse["ql_mean"] = 8.8699692429062537e+01
-best_mse["thetal_mean"] = 1.1916582780639058e-04
-best_mse["Hvar_mean"] = 2.0710735541431495e+03
-best_mse["QTvar_mean"] = 8.8908792485557092e+02
+best_mse["qt_mean"] = 4.2359420114784679e-01
+best_mse["updraft_area"] = 1.6722742675049749e+03
+best_mse["updraft_w"] = 2.8300484466673169e+02
+best_mse["updraft_qt"] = 1.7904453430368505e+01
+best_mse["updraft_thetal"] = 6.3573804437983789e+01
+best_mse["v_mean"] = 1.0634096403180666e+02
+best_mse["u_mean"] = 1.1434656097951681e+02
+best_mse["tke_mean"] = 8.6824438228445820e+02
+best_mse["temperature_mean"] = 1.2576042829475203e-04
+best_mse["ql_mean"] = 8.9402683009257942e+01
+best_mse["thetal_mean"] = 1.1896802334339157e-04
+best_mse["Hvar_mean"] = 2.0228054123757493e+03
+best_mse["QTvar_mean"] = 8.7748179858993274e+02
 
 @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"] = 1.9338238569368907e-01
-best_mse["updraft_area"] = 4.8852480428600916e+02
-best_mse["updraft_w"] = 3.1250380743373622e+01
-best_mse["updraft_qt"] = 1.1136647638388030e+01
-best_mse["updraft_thetal"] = 2.2395194056858262e+01
-best_mse["u_mean"] = 7.2109941452760484e+02
-best_mse["tke_mean"] = 4.7919320017479706e+01
-best_mse["temperature_mean"] = 1.4374232234123472e-05
-best_mse["thetal_mean"] = 1.4379965951428025e-05
-best_mse["Hvar_mean"] = 2.3305430832031414e+02
-
+best_mse["qt_mean"] = 1.9338238625419904e-01
+best_mse["updraft_area"] = 4.8852467389339699e+02
+best_mse["updraft_w"] = 3.1250501928635639e+01
+best_mse["updraft_qt"] = 1.1136647539128006e+01
+best_mse["updraft_thetal"] = 2.2395194054335001e+01
+best_mse["u_mean"] = 7.2109940121868794e+02
+best_mse["tke_mean"] = 4.7919324815294324e+01
+best_mse["temperature_mean"] = 1.4374232147002715e-05
+best_mse["thetal_mean"] = 1.4379965874656967e-05
+best_mse["Hvar_mean"] = 2.3305425470797826e+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"] = 1.7430143247852277e+00
-best_mse["updraft_area"] = 3.2544587990881948e+04
-best_mse["updraft_w"] = 9.2113862144403333e+02
-best_mse["updraft_qt"] = 3.0303277939652101e+01
-best_mse["updraft_thetal"] = 1.1108868468322090e+02
-best_mse["v_mean"] = 2.9291515854682558e+02
-best_mse["u_mean"] = 1.6874269993592438e+03
-best_mse["tke_mean"] = 1.5877278080585038e+03
-best_mse["temperature_mean"] = 6.3087898742812351e-04
-best_mse["ql_mean"] = 9.4714318856310877e+02
-best_mse["thetal_mean"] = 4.1952114677068139e-04
-best_mse["Hvar_mean"] = 8.1299695367399718e+03
-best_mse["QTvar_mean"] = 2.7048653627007479e+03
+best_mse["qt_mean"] = 1.7430285968163504e+00
+best_mse["updraft_area"] = 3.2545401197014067e+04
+best_mse["updraft_w"] = 9.2113782722003907e+02
+best_mse["updraft_qt"] = 3.0306781659879494e+01
+best_mse["updraft_thetal"] = 1.1108867837608678e+02
+best_mse["v_mean"] = 2.9291577336310036e+02
+best_mse["u_mean"] = 1.6874237352911189e+03
+best_mse["tke_mean"] = 1.5857027462501208e+03
+best_mse["temperature_mean"] = 6.3088676948654218e-04
+best_mse["ql_mean"] = 9.4722581733362244e+02
+best_mse["thetal_mean"] = 4.1952296308447006e-04
+best_mse["Hvar_mean"] = 8.1275664794363738e+03
+best_mse["QTvar_mean"] = 2.7043370821860535e+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