WIP

WIP add pressure closure file cleanups in entrainment add entr_detr, incomplete move nh_pressure to closures move nh_pressure to closures update tables running with entr in closures running with entr in closures working with entr_detr apply formatter update tables resolve conflicts with yc/call_press_closure remove cleanup_covariance add entr_detr_model + merge turb entr + cleanup fix bug apply formatter
CliMA · Aug 10, 2021 · 3fa8046 · 3fa8046
1 parent 94a3f72
commit 3fa8046
Show file tree

Hide file tree

Showing 7 changed files with 221 additions and 333 deletions.
diff --git a/integration_tests/utils/generate_namelist.jl b/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

diff --git a/src/TurbulenceConvection.jl b/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
diff --git a/src/Turbulence_PrognosticTKE.jl b/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,103 +900,95 @@ 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]
-            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]
-
+            w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
+            a_upd = interpf2c(self.UpdVar.Area.values, grid, k, i)
+            if w_upd * a_upd > 0.0
                 # 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 = a_upd * (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.entrainment_ed_mf_sigma,
+                    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.frac_turb_entr[i, k] = 0.0
                 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
-            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
+            # end
 
             # pressure
             a_kfull = interpc2f(self.UpdVar.Area.values, grid, k, i)
@@ -1033,16 +1023,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 +1765,6 @@ function compute_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Ca
         self.UpdVar.H,
         self.UpdVar.QT,
     )
-    cleanup_covariance(self, GMV)
     return
 end
 
@@ -1872,49 +1851,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,

diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl
@@ -0,0 +1,85 @@
+#### Entrainment-Detrainment kernels
+"""
+    entr_detr(
+        param_set,
+        w_min,
+        β,
+        c_δ,
+        c_div,
+        c_μ,
+        c_μ0,
+        c_ed_mf,
+        χ_upd,
+        c_λ,
+        εδ_model,
+    )
+
+Returns the dynamic entrainment and detrainment rates,
+as well as the turbulent entrainment rate, following
+Cohen et al. (JAMES, 2020), given:
+ - `param_set`: parameter set
+ - `w_min`: minimum veritcal velocity
+ - `β`: sorting power for moist mixing
+ - `c_δ`: detrainment factor
+ - `c_div`: divergence factor for bubble case (zero otherwise)
+ - `c_μ`: logisitc function scale 
+ - `c_μ0`: logisitc function timescale
+ - `c_ed_mf`: ED-MF coeeficeicent for relating BL tke and velocity scales
+ - `χ_upd`: updraft mixing fraction
+ - `c_λ`: 
+ - `εδ_model`: updraft buoyancy
+ -- Where:
+ - `εδ_model.b_upd`: updraft buoyancy
+ - `εδ_model.w_upd`: updraft vertical velocity
+ - `εδ_model.a_upd`: updraft area fraction
+ - `εδ_model.a_env`: environment buoyancy
+ - `εδ_model.b_env`: environment vertical velocity
+ - `εδ_model.w_env`: environment area fraction
+ - `εδ_model.ql_up`: updraft liquid water
+ - `εδ_model.ql_env`: environment liquid water 
+ - `εδ_model.RH_upd`: updraft relative humidity
+ - `εδ_model.RH_env`: environment relative humidity
+ - `εδ_model.M`: updraft momentum
+ - `εδ_model.dMdz`: updraft momentum divergence
+ - `εδ_model.tke`: env TKE
+ - `εδ_model.N_up`: total number of updrafts
+ - `εδ_model.ρ`: referance density
+"""
+function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, c_ed_mf, χ_upd, c_λ, ε_t, εδ_model)
+
+    l = zeros(2)
+    c_ε = CPEDMF.c_ε(param_set)
+
+    if (εδ_model.ql_upd + εδ_model.ql_env) == 0.0
+        c_δ = 0.0
+    end
+
+    Δw = εδ_model.w_upd - εδ_model.w_env
+    if Δw < 0.0
+        Δw -= w_min
+    else
+        Δw += w_min
+    end
+
+    Δb = (εδ_model.b_upd - εδ_model.b_env)
+
+    D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model)
+
+    l[1] = c_λ * fabs(Δb / sqrt(εδ_model.tke + 1e-8))
+    l[2] = fabs(Δb / Δw)
+    λ = lamb_smooth_minimum(l, 0.1, 0.0005)
+
+    MdMdz = fmax(εδ_model.dMdz / fmax(εδ_model.M, 1e-12), 0.0)
+    MdMdz = fmax(-εδ_model.dMdz / fmax(εδ_model.M, 1e-12), 0.0)
+
+
+
+    # turbulent entrainment
+    k_ε = εδ_model.a_upd * ε_t * sqrt(fmax(εδ_model.tke, 0.0)) * εδ_model.R_up
+    ε_turb = (2.0 / εδ_model.R_up^2.0) * k_ε / (εδ_model.w_upd * εδ_model.a_upd)
+
+    ε_dyn = λ / Δw * (c_ε * D_ε + c_δ * M_ε) + MdMdz * c_div
+    δ_dyn = λ / Δw * (c_ε * D_δ + c_δ * M_δ) + MdMdz * c_div
+
+    return ε_dyn, δ_dyn, ε_turb, k_ε
+end