Move entr and detr to separate file

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.4297134939930601e-01
-best_mse["updraft_area"] = 1.9829008493073268e+03
-best_mse["updraft_w"] = 3.1072236638537851e+02
-best_mse["updraft_qt"] = 1.2698381273033208e+01
-best_mse["updraft_thetal"] = 2.7686390354792074e+01
-best_mse["u_mean"] = 8.7998547277817863e+01
-best_mse["tke_mean"] = 6.7319699431828644e+02
-best_mse["temperature_mean"] = 1.3960321509272764e-04
-best_mse["ql_mean"] = 2.9854035802671558e+02
-best_mse["thetal_mean"] = 1.4200441599983529e-04
-best_mse["Hvar_mean"] = 1.3755159256021079e+03
-best_mse["QTvar_mean"] = 2.6566877259844244e+02
+best_mse["qt_mean"] = 3.4420193413687117e-01
+best_mse["updraft_area"] = 1.9818317322584239e+03
+best_mse["updraft_w"] = 3.1101869507553300e+02
+best_mse["updraft_qt"] = 1.2779005319987142e+01
+best_mse["updraft_thetal"] = 2.7687174775467490e+01
+best_mse["u_mean"] = 8.7998547277817835e+01
+best_mse["tke_mean"] = 6.7532125088372368e+02
+best_mse["temperature_mean"] = 1.3976231155597154e-04
+best_mse["ql_mean"] = 2.9406894196734129e+02
+best_mse["thetal_mean"] = 1.4216021128644291e-04
+best_mse["Hvar_mean"] = 1.4219930941047571e+03
+best_mse["QTvar_mean"] = 2.7979236616322686e+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"] = 9.5080504754167433e-02
-best_mse["updraft_area"] = 7.2534663873051477e+02
-best_mse["updraft_w"] = 2.6358909128208378e+01
-best_mse["updraft_qt"] = 4.0334302720334296e+00
-best_mse["updraft_thetal"] = 2.1547995304586468e+01
-best_mse["v_mean"] = 6.5231550206742128e+01
-best_mse["u_mean"] = 5.3292889729286458e+01
-best_mse["tke_mean"] = 3.8860563196392924e+01
-best_mse["temperature_mean"] = 3.8436817065445734e-05
-best_mse["ql_mean"] = 3.7878522086312807e+00
-best_mse["thetal_mean"] = 3.9063754163395865e-05
-best_mse["Hvar_mean"] = 7.2639655951863034e+01
-best_mse["QTvar_mean"] = 2.2435640414468629e+01
+best_mse["qt_mean"] = 9.5377282321874049e-02
+best_mse["updraft_area"] = 7.2522420468798543e+02
+best_mse["updraft_w"] = 2.6282677670021386e+01
+best_mse["updraft_qt"] = 4.0341139540222635e+00
+best_mse["updraft_thetal"] = 2.1548164483602367e+01
+best_mse["v_mean"] = 6.5247595142753596e+01
+best_mse["u_mean"] = 5.3292989796462209e+01
+best_mse["tke_mean"] = 3.8891129116688681e+01
+best_mse["temperature_mean"] = 3.8547176978568602e-05
+best_mse["ql_mean"] = 3.8134281154981489e+00
+best_mse["thetal_mean"] = 3.9162071364386570e-05
+best_mse["Hvar_mean"] = 7.5510623226760771e+01
+best_mse["QTvar_mean"] = 2.2962565404185188e+01
 
 @testset "Bomex" begin
     println("Running Bomex...")

integration_tests/Nieuwstadt.jl

@@ -11,14 +11,14 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["updraft_area"] = 5.9412545314262150e+02
-best_mse["updraft_w"] = 2.6952552348430249e+01
-best_mse["updraft_thetal"] = 3.0475271481866773e+01
-best_mse["u_mean"] = 1.5247657003047263e+02
-best_mse["tke_mean"] = 7.3600391357561151e+01
-best_mse["temperature_mean"] = 1.1971868506884125e-05
-best_mse["thetal_mean"] = 1.2117924470979204e-05
-best_mse["Hvar_mean"] = 1.8622261767371862e+02
+best_mse["updraft_area"] = 5.9412544437504687e+02
+best_mse["updraft_w"] = 2.6952552450364237e+01
+best_mse["updraft_thetal"] = 3.0475271481865867e+01
+best_mse["u_mean"] = 1.5247656989308092e+02
+best_mse["tke_mean"] = 7.3600391623055913e+01
+best_mse["temperature_mean"] = 1.1971868653164508e-05
+best_mse["thetal_mean"] = 1.2117924617884567e-05
+best_mse["Hvar_mean"] = 1.8622261701198534e+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.2888329142054360e-01
-best_mse["updraft_area"] = 1.5804026921690113e+03
-best_mse["updraft_w"] = 3.0412249337738064e+02
-best_mse["updraft_qt"] = 1.8041619418228048e+01
-best_mse["updraft_thetal"] = 6.3565634303932690e+01
-best_mse["v_mean"] = 1.0637419570945200e+02
-best_mse["u_mean"] = 1.1444385357013196e+02
-best_mse["tke_mean"] = 8.6459012010773438e+02
-best_mse["temperature_mean"] = 1.2577596756983134e-04
-best_mse["ql_mean"] = 5.0938207504953844e+01
-best_mse["thetal_mean"] = 1.2604124901901089e-04
-best_mse["Hvar_mean"] = 1.4598483199237567e+03
-best_mse["QTvar_mean"] = 7.0850350571186209e+02
+best_mse["qt_mean"] = 4.2745465446436787e-01
+best_mse["updraft_area"] = 1.5810580213880007e+03
+best_mse["updraft_w"] = 3.0293830804824944e+02
+best_mse["updraft_qt"] = 1.8036840808532091e+01
+best_mse["updraft_thetal"] = 6.3565923109974378e+01
+best_mse["v_mean"] = 1.0637284782337001e+02
+best_mse["u_mean"] = 1.1444354015737491e+02
+best_mse["tke_mean"] = 8.6553999307149547e+02
+best_mse["temperature_mean"] = 1.2540517259642333e-04
+best_mse["ql_mean"] = 5.1325163153016852e+01
+best_mse["thetal_mean"] = 1.2561593164473264e-04
+best_mse["Hvar_mean"] = 1.4739277424672016e+03
+best_mse["QTvar_mean"] = 7.1366176286414179e+02
 
 @testset "Rico" begin
     println("Running Rico...")

integration_tests/Soares.jl

@@ -11,16 +11,16 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 1.9368730351926020e-01
-best_mse["updraft_area"] = 4.8515130524256091e+02
-best_mse["updraft_w"] = 3.1260485555848359e+01
-best_mse["updraft_qt"] = 1.1137263871861222e+01
-best_mse["updraft_thetal"] = 2.2395183773350311e+01
-best_mse["u_mean"] = 7.2111916707127898e+02
-best_mse["tke_mean"] = 4.7882808390238374e+01
-best_mse["temperature_mean"] = 1.3818642271008120e-05
-best_mse["thetal_mean"] = 1.4383770217873919e-05
-best_mse["Hvar_mean"] = 2.3235278609470686e+02
+best_mse["qt_mean"] = 1.9368730493269354e-01
+best_mse["updraft_area"] = 4.8515131142491634e+02
+best_mse["updraft_w"] = 3.1260485803642322e+01
+best_mse["updraft_qt"] = 1.1137263873078020e+01
+best_mse["updraft_thetal"] = 2.2395183773363378e+01
+best_mse["u_mean"] = 7.2111916666642423e+02
+best_mse["tke_mean"] = 4.7882808859338290e+01
+best_mse["temperature_mean"] = 1.3818642368852120e-05
+best_mse["thetal_mean"] = 1.4383770315305301e-05
+best_mse["Hvar_mean"] = 2.3235278783484819e+02
 
 @testset "Soares" begin
     println("Running Soares...")

integration_tests/utils/generate_namelist.jl

@@ -57,7 +57,7 @@ function default_namelist(case_name::String)
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["max_area"] = 0.9
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"] = 0.13
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["detrainment_factor"] = 0.51
-    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.4
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.0
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["turbulent_entrainment_factor"] = 0.015
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_ed_mf_sigma"] = 50.0
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_smin_tke_coeff"] = 0.3
@@ -345,7 +345,7 @@ function DryBubble(namelist_defaults)
     namelist["meta"]["simname"] = "DryBubble"
     namelist["meta"]["casename"] = "DryBubble"
 
-    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment"] = "moisture_deficit_div"
+    namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.4
 
     return namelist
 end

src/TurbulenceConvection.jl

@@ -73,5 +73,7 @@ include("forcing_functions.jl")
 include("Surface.jl")
 include("surface_functions.jl")
 include("closures/perturbation_pressure.jl")
+include("closures/entr_detr.jl")
+include("closures/nondimensional_exchange_functions.jl")
 
 end

src/Turbulence_PrognosticTKE.jl

@@ -146,8 +146,6 @@ function io(self::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats, TS::TimeStepping)
     write_profile(Stats, "horiz_K_eddy", mean_horiz_K_eddy[cinterior])
     write_profile(Stats, "entrainment_sc", mean_entr_sc[cinterior])
     write_profile(Stats, "detrainment_sc", mean_detr_sc[cinterior])
-    write_profile(Stats, "sorting_function", mean_sorting_function[cinterior])
-    write_profile(Stats, "b_mix", mean_b_mix[cinterior])
     write_profile(Stats, "nh_pressure", mean_nh_pressure[cinterior])
     write_profile(Stats, "nh_pressure_adv", mean_nh_pressure_adv[cinterior])
     write_profile(Stats, "nh_pressure_drag", mean_nh_pressure_drag[cinterior])
@@ -902,102 +900,71 @@ function get_env_covar_from_GMV(
 end
 
 function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase)
-    ret = entr_struct()
-    input = entr_in_struct()
+    εδ_model = entr_detr_model()
     sa = eos_struct()
     quadrature_order = 3
     grid = get_grid(self)
     param_set = parameter_set(self)
 
     upd_cloud_diagnostics(self.UpdVar, reference_state(self))
 
-    input.wstar = self.wstar
-
-    input.dz = get_grid(self).dz
-    input.zbl = compute_zbl_qt_grad(GMV)
-    input.sort_pow = self.sorting_power
-    input.c_det = self.detrainment_factor
-    input.c_div = self.entrainment_Mdiv_factor
-    input.c_mu = self.entrainment_sigma
-    input.c_mu0 = self.entrainment_scale
-    input.c_ed_mf = self.entrainment_ed_mf_sigma
-    input.chi_upd = self.updraft_mixing_frac
-    input.tke_coef = self.entrainment_smin_tke_coeff
-
     # TODO: merge this with grid.cinterior
     kc_surf = kc_surface(grid)
     kc_toa = kc_top_of_atmos(grid)
     cinterior = kc_surf:kc_toa
 
     ref_state = reference_state(self)
-    tau = get_mixing_tau(self.zi, self.wstar)
-
-    ae = 1 .- self.UpdVar.Area.bulkvalues
     l = zeros(2)
 
     @inbounds for k in real_center_indicies(grid)
         @inbounds for i in xrange(self.n_updrafts)
-            alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])), 1)
-            avals = self.UpdVar.Area.values[i, cinterior]
-            input.zi = self.UpdVar.cloud_base[i]
             # entrainment
-            input.buoy_ed_flux = self.EnvVar.TKE.buoy[k]
+            w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
             if self.UpdVar.Area.values[i, k] > 0.0
-                input.b_upd = self.UpdVar.B.values[i, k]
-                input.w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
-                input.z = grid.z_half[k]
-                input.a_upd = self.UpdVar.Area.values[i, k]
-                input.a_env = 1.0 - self.UpdVar.Area.bulkvalues[k]
-                input.tke = self.EnvVar.TKE.values[k]
-                input.ql_env = self.EnvVar.QL.values[k]
-                input.b_env = self.EnvVar.B.values[k]
-                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.RH_upd = self.UpdVar.RH.values[i, k]
-                input.RH_env = self.EnvVar.RH.values[k]
-
                 # compute dMdz at half levels
                 gmv_w_k = interpf2c(GMV.W.values, grid, k)
                 gmv_w_km = interpf2c(GMV.W.values, grid, k - 1)
                 upd_w_km = interpf2c(self.UpdVar.W.values, grid, k - 1, i)
-                input.M = input.a_upd * (input.w_upd - gmv_w_k)
+                M = self.UpdVar.Area.values[i, k] * (w_upd - gmv_w_k)
                 Mm = self.UpdVar.Area.values[i, k - 1] * (upd_w_km - gmv_w_km)
-                input.dMdz = (input.M - Mm) * grid.dzi
-                input.tke = self.EnvVar.TKE.values[k]
-
-                ret = self.entr_detr_fp(param_set, input)
-                self.entr_sc[i, k] = ret.entr_sc
-                self.detr_sc[i, k] = ret.detr_sc
-                self.sorting_function[i, k] = ret.sorting_function
-                self.b_mix[i, k] = ret.b_mix
-
+                dMdz = (M - Mm) * grid.dzi
+                w_min = 0.001
+
+                εδ_model.b_upd = self.UpdVar.B.values[i, k]
+                εδ_model.b_env = self.EnvVar.B.values[k]
+                εδ_model.w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
+                εδ_model.w_env = interpf2c(self.EnvVar.W.values, grid, k)
+                εδ_model.a_upd = self.UpdVar.Area.values[i, k]
+                εδ_model.a_env = self.EnvVar.Area.values[k]
+                εδ_model.ql_upd = self.UpdVar.QL.values[i, k]
+                εδ_model.ql_env = self.EnvVar.QL.values[k]
+                εδ_model.RH_upd = self.UpdVar.RH.values[i, k]
+                εδ_model.RH_env = self.EnvVar.RH.values[k]
+                εδ_model.M = M
+                εδ_model.dMdz = dMdz
+                εδ_model.tke = self.EnvVar.TKE.values[k]
+                εδ_model.n_up = self.n_updrafts
+                εδ_model.ρ = ref_state.rho0[k]
+                εδ_model.R_up = self.pressure_plume_spacing[i]
+
+                self.entr_sc[i, k], self.detr_sc[i, k], self.frac_turb_entr[i, k], self.horiz_K_eddy[i, k] = entr_detr(
+                    param_set,
+                    w_min,
+                    self.sorting_power,
+                    self.detrainment_factor,
+                    self.entrainment_Mdiv_factor,
+                    self.entrainment_sigma,
+                    self.entrainment_scale,
+                    self.updraft_mixing_frac,
+                    self.entrainment_smin_tke_coeff,
+                    self.turbulent_entrainment_factor,
+                    εδ_model,
+                )
             else
                 self.entr_sc[i, k] = 0.0
                 self.detr_sc[i, k] = 0.0
-                self.sorting_function[i, k] = 0.0
-                self.b_mix[i, k] = self.EnvVar.B.values[k]
-            end
-
-            # horizontal eddy diffusivities
-            if self.UpdVar.Area.values[i, k] > 0.0
-                self.horiz_K_eddy[i, k] =
-                    self.UpdVar.Area.values[i, k] *
-                    self.turbulent_entrainment_factor *
-                    sqrt(fmax(self.EnvVar.TKE.values[k], 0.0)) *
-                    self.pressure_plume_spacing[i]
-            else
-                self.horiz_K_eddy[i, k] = 0.0
-            end
-
-            # turbulent entrainment
-            wu_half = interp2pt(self.UpdVar.W.values[i, k], self.UpdVar.W.values[i, k - 1])
-            if self.UpdVar.Area.values[i, k] * wu_half > 0.0
-                self.frac_turb_entr[i, k] =
-                    (2.0 / self.pressure_plume_spacing[i]^2.0) * self.horiz_K_eddy[i, k] / wu_half /
-                    self.UpdVar.Area.values[i, k]
-            else
                 self.frac_turb_entr[i, k] = 0.0
+                self.horiz_K_eddy[i, k] = 0.0
             end
 
             # pressure
@@ -1033,16 +1000,6 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
     return
 end
 
-function compute_zbl_qt_grad(GMV::GridMeanVariables)
-    grid = get_grid(GMV)
-    # computes inversion height as z with max gradient of qt
-    z∇q_tot = map(real_center_indicies(grid)) do k
-        (grid.z_half[k], abs(∇c2f(GMV.QT.values, grid, k)))
-    end
-    k_star = argmax(last.(z∇q_tot))
-    return first.(z∇q_tot)[k_star]
-end
-
 function compute_pressure_plume_spacing(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase)
 
     @inbounds for i in xrange(self.n_updrafts)
@@ -1785,7 +1742,6 @@ function compute_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Ca
         self.UpdVar.H,
         self.UpdVar.QT,
     )
-    cleanup_covariance(self, GMV)
     return
 end
 
@@ -1872,49 +1828,6 @@ function initialize_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables,
     return
 end
 
-function cleanup_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables)
-
-    tmp_eps = 1e-18
-    grid = get_grid(self)
-
-    @inbounds for k in real_center_indicies(grid)
-        if GMV.TKE.values[k] < tmp_eps
-            GMV.TKE.values[k] = 0.0
-        end
-        if GMV.Hvar.values[k] < tmp_eps
-            GMV.Hvar.values[k] = 0.0
-        end
-        if GMV.QTvar.values[k] < tmp_eps
-            GMV.QTvar.values[k] = 0.0
-        end
-        if fabs(GMV.HQTcov.values[k]) < tmp_eps
-            GMV.HQTcov.values[k] = 0.0
-        end
-        if fabs(GMV.W_third_m.values[k]) < tmp_eps
-            GMV.W_third_m.values[k] = 0.0
-        end
-        if fabs(GMV.H_third_m.values[k]) < tmp_eps
-            GMV.H_third_m.values[k] = 0.0
-        end
-        if fabs(GMV.QT_third_m.values[k]) < tmp_eps
-            GMV.QT_third_m.values[k] = 0.0
-        end
-        if self.EnvVar.Hvar.values[k] < tmp_eps
-            self.EnvVar.Hvar.values[k] = 0.0
-        end
-        if self.EnvVar.TKE.values[k] < tmp_eps
-            self.EnvVar.TKE.values[k] = 0.0
-        end
-        if self.EnvVar.QTvar.values[k] < tmp_eps
-            self.EnvVar.QTvar.values[k] = 0.0
-        end
-        if fabs(self.EnvVar.HQTcov.values[k]) < tmp_eps
-            self.EnvVar.HQTcov.values[k] = 0.0
-        end
-    end
-end
-
-
 function compute_covariance_shear(
     self::EDMF_PrognosticTKE,
     GMV::GridMeanVariables,