wip

CliMA · Apr 6, 2024 · 73c9f59 · 73c9f59
1 parent abdb213
commit 73c9f59
Show file tree

Hide file tree

Showing 4 changed files with 42 additions and 38 deletions.
diff --git a/config/model_configs/prognostic_edmfx_bomex_column.yml b/config/model_configs/prognostic_edmfx_bomex_column.yml
@@ -10,8 +10,8 @@ implicit_sgs_advection: false
 approximate_linear_solve_iters: 2
 max_newton_iters_ode: 3
 edmfx_upwinding: first_order
-edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_entr_model: "PiGroups"
+edmfx_detr_model: "PiGroups"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true

diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -182,7 +182,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
         ᶜK_h,
         ρatke_flux,
     ) = p.precomputed
-    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p.precomputed
+    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs, ᶜtke⁰) = p.precomputed
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
 
     ᶜz = Fields.coordinate_field(Y.c).z
@@ -200,14 +200,14 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             z_sfc,
             ᶜp,
             Y.c.ρ,
-            buoyancy_flux,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
             TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
             ᶜphysical_buoyancy(params, Y.c.ρ, ᶜρʲs.:($$j)),
             get_physical_w(ᶜu, ᶜlg),
             TD.relative_humidity(thermo_params, ᶜts⁰),
             FT(0),
+            max(ᶜtke⁰, 0),
             p.atmos.edmfx_entr_model,
         )
         @. ᶜentrʲs.:($$j) = limit_entrainment(
@@ -224,7 +224,6 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             z_sfc,
             ᶜp,
             Y.c.ρ,
-            buoyancy_flux,
             Y.c.sgsʲs.:($$j).ρa,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
@@ -236,6 +235,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             ᶜentrʲs.:($$j),
             ᶜvert_div,
             ᶜmassflux_vert_div,
+            max(ᶜtke⁰, 0),
             p.atmos.edmfx_detr_model,
         )
         @. ᶜdetrʲs.:($$j) = limit_detrainment(

diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -2,9 +2,6 @@
 ##### EDMF entrainment detrainment
 #####
 
-# return a harmonic mean of (a, 1-a)
-hm_limiter(a) = 2 * a * (1 - a)
-
 """
    Return entrainment rate [1/s].
 
@@ -25,14 +22,14 @@ function entrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
     ᶜbuoyʲ::FT,
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜtke⁰::FT,
     ::NoEntrainment,
 ) where {FT}
     return FT(0)
@@ -44,38 +41,34 @@ function entrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
     ᶜbuoyʲ::FT,
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜtke⁰::FT,
     ::PiGroupsEntrainment,
 ) where {FT}
     if ᶜaʲ <= FT(0)
         return FT(0)
     else
         g = CAP.grav(params)
-
-        # pressure scale height (height where pressure drops by 1/e)
         ref_H = ᶜp / (ᶜρ * g)
-        # convective velocity
-        w_star = get_wstar(buoy_flux_surface)
 
         # non-dimensional pi-groups
-        # TODO - using Π₁ blows things up
         Π₁ =
             (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
-            ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
+            ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 100
+        Π₂ = ᶜtke⁰ / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 2
         Π₃ = sqrt(ᶜaʲ)
         Π₄ = ᶜRHʲ - ᶜRH⁰
-        Π₆ = (ᶜz - z_sfc) / ref_H
+        Π₅ = (ᶜz - z_sfc) / ref_H
         entr =
-            abs(ᶜwʲ) / (ᶜz - z_sfc) * (
-                -4.013288 - 0.000968 * Π₁ + 0.356974 * Π₃ - 0.403124 * Π₄ +
-                1.503261 * Π₆
+            abs(ᶜwʲ - ᶜw⁰) / (ᶜz - z_sfc) * (
+                -0.32332 + 4.79372 * Π₁ + 3.98108 * Π₂ - 21.64173 * Π₃ + 18.395 * Π₄ +
+                1.12799 * Π₅
             )
 
         return entr
@@ -88,14 +81,14 @@ function entrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
     ᶜbuoyʲ::FT,
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜtke⁰::FT,
     ::GeneralizedEntrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -146,6 +139,7 @@ function detrainment(
     ᶜentr::FT,
     ᶜvert_div::FT,
     ᶜmassflux_vert_div::FT,
+    ᶜtke⁰::FT,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -157,7 +151,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
@@ -169,34 +162,28 @@ function detrainment(
     ᶜentr::FT,
     ᶜvert_div::FT,
     ᶜmassflux_vert_div::FT,
+    ᶜtke⁰::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
 
     if ᶜaʲ <= FT(0)
         return FT(0)
     else
         g = CAP.grav(params)
-        turbconv_params = CAP.turbconv_params(params)
-        ᶜaʲ_max = CAP.max_area(turbconv_params)
-
-        max_area_limiter = FT(0.1) * exp(-10 * (ᶜaʲ_max - ᶜaʲ))
-        # pressure scale height (height where pressure drops by 1/e)
         ref_H = ᶜp / (ᶜρ * g)
-        # convective velocity
-        w_star = get_wstar(buoy_flux_surface)
 
         # non-dimensional pi-groups
-        # TODO - using Π₁ blows things up
         Π₁ =
             (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
-            ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
+            ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 100
+        Π₂ = ᶜtke⁰ / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 2
         Π₃ = sqrt(ᶜaʲ)
         Π₄ = ᶜRHʲ - ᶜRH⁰
-        Π₆ = (ᶜz - z_sfc) / ref_H
+        Π₅ = (ᶜz - z_sfc) / ref_H
         detr =
-            abs(ᶜwʲ) * (
-                3.535208 + 0.598496 * Π₁ + 1.583348 * Π₃ + 0.046275 * Π₄ -
-                0.344836 * Π₆ + max_area_limiter
+            -min(ᶜmassflux_vert_div, 0) / ᶜρaʲ * (
+                0.3410 - 0.56153 * Π₁ - 0.53411 * Π₂ + 6.01925 * Π₃ - 1.47516 * Π₄ -
+                3.85788 * Π₅
             )
         return detr
     end
@@ -208,7 +195,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
@@ -220,6 +206,7 @@ function detrainment(
     ᶜentr::FT,
     ᶜvert_div::FT,
     ᶜmassflux_vert_div::FT,
+    ᶜtke⁰::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -258,7 +245,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,

diff --git a/toml/prognostic_edmfx_bomex.toml b/toml/prognostic_edmfx_bomex.toml
@@ -41,4 +41,22 @@ value = 1
 value = 0
 
 [pressure_normalmode_drag_coeff]
-value = 50.0
+value = 20.0
+
+[pressure_normalmode_buoy_coeff1]
+value = 0.857
+
+[mixing_length_eddy_viscosity_coefficient]
+value = 0.041
+
+[mixing_length_diss_coeff]
+value = 0.174
+
+[mixing_length_static_stab_coeff]
+value = 0.292
+
+[mixing_length_tke_surf_scale]
+value = 9.609
+
+[mixing_length_Prandtl_number_0]
+value = 0.515