add prognostic tke equation for prognostic edmf

config/model_configs/prognostic_edmfx_bomex_box.yml

@@ -11,7 +11,7 @@ edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
-prognostic_tke: false
+prognostic_tke: true
 moist: "equil"
 config: "box"
 hyperdiff: "true"
@@ -24,78 +24,78 @@ y_elem: 2
 z_elem: 60
 z_stretch: false
 perturb_initstate: false
-dt: "5secs"
+dt: "2secs"
 t_end: "6hours"
-dt_save_to_disk: "10mins"
+dt_save_to_disk: "2mins"
 toml: [toml/prognostic_edmfx_box.toml]
 diagnostics:
   - short_name: ts
-    period: 10mins
+    period: 1hours
   - short_name: ta
-    period: 10mins
+    period: 1hours
   - short_name: thetaa
-    period: 10mins
+    period: 1hours
   - short_name: ha
-    period: 10mins
+    period: 1hours
   - short_name: pfull
-    period: 10mins
+    period: 1hours
   - short_name: rhoa
-    period: 10mins
+    period: 1hours
   - short_name: ua
-    period: 10mins
+    period: 1hours
   - short_name: va
-    period: 10mins
+    period: 1hours
   - short_name: wa
-    period: 10mins
+    period: 1hours
   - short_name: hur
-    period: 10mins
+    period: 1hours
   - short_name: hus
-    period: 10mins
+    period: 1hours
   - short_name: clw
-    period: 10mins
+    period: 1hours
   - short_name: cli
-    period: 10mins
+    period: 1hours
   - short_name: hussfc
-    period: 10mins
+    period: 1hours
   - short_name: evspsbl
-    period: 10mins
+    period: 1hours
   - short_name: arup
-    period: 10mins
+    period: 1hours
   - short_name: waup
-    period: 10mins
+    period: 1hours
   - short_name: taup
-    period: 10mins
+    period: 1hours
   - short_name: thetaaup
-    period: 10mins
+    period: 1hours
   - short_name: haup
-    period: 10mins
+    period: 1hours
   - short_name: husup
-    period: 10mins
+    period: 1hours
   - short_name: hurup
-    period: 10mins
+    period: 1hours
   - short_name: clwup
-    period: 10mins
+    period: 1hours
   - short_name: cliup
-    period: 10mins
+    period: 1hours
   - short_name: rhoaen
-    period: 10mins
+    period: 1hours
   - short_name: waen
-    period: 10mins
+    period: 1hours
   - short_name: taen
-    period: 10mins
+    period: 1hours
   - short_name: thetaaen
-    period: 10mins
+    period: 1hours
   - short_name: haen
-    period: 10mins
+    period: 1hours
   - short_name: husen
-    period: 10mins
+    period: 1hours
   - short_name: huren
-    period: 10mins
+    period: 1hours
   - short_name: clwen
-    period: 10mins
+    period: 1hours
   - short_name: clien
-    period: 10mins
+    period: 1hours
   - short_name: tke
-    period: 10mins
+    period: 1hours
   - short_name: lmix
-    period: 10mins
+    period: 1hours

config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml

@@ -11,7 +11,7 @@ edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
-prognostic_tke: false
+prognostic_tke: true
 moist: equil
 config: box
 hyperdiff: "true"
@@ -23,7 +23,7 @@ x_elem: 2
 y_elem: 2
 z_elem: 30
 z_stretch: false
-dt: 5secs
+dt: 2secs
 t_end: 4hours
-dt_save_to_disk: 10mins
+dt_save_to_disk: 2mins
 toml: [toml/prognostic_edmfx_box.toml]

config/model_configs/prognostic_edmfx_gabls_box.yml

@@ -9,7 +9,7 @@ edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
-prognostic_tke: false
+prognostic_tke: true
 moist: "equil"
 config: "box"
 hyperdiff: "true"

config/model_configs/prognostic_edmfx_rico_box.yml

@@ -10,7 +10,7 @@ edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
-prognostic_tke: false
+prognostic_tke: true
 moist: "equil"
 precip_model: "0M"
 config: "box"

config/model_configs/prognostic_edmfx_trmm_box.yml

@@ -9,7 +9,7 @@ edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
-prognostic_tke: false
+prognostic_tke: true
 moist: equil
 apply_limiter: false
 precip_model: "0M"

src/cache/precomputed_quantities.jl

@@ -78,6 +78,7 @@ function precomputed_quantities(Y, atmos)
             ᶜmixing_length = similar(Y.c, FT),
             ᶜK_u = similar(Y.c, FT),
             ᶜK_h = similar(Y.c, FT),
+            ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
             ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}),
             ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
             ᶜKʲs = similar(Y.c, NTuple{n, FT}),

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -22,9 +22,9 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     n = n_mass_flux_subdomains(turbconv_model)
     thermo_args = (thermo_params, energy_form, moisture_model)
 
-    (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p
+    (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref, ᶜu) = p
     (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
-    (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p
+    (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h, ρatke_flux) = p
     (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p
     (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
 
@@ -226,9 +226,7 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
         @. ᶜtke_exch +=
             Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ * (
                 1 / 2 *
-                (
-                    get_physical_w(ᶜuʲs.:($$j), ᶜlg) - get_physical_w(ᶜu⁰, ᶜlg)
-                )^2 - ᶜtke⁰
+                norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰
             )
     end
 
@@ -239,9 +237,9 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
         ᶜz,
         z_sfc,
         ᶜdz,
-        sfc_tke,
+        max(sfc_tke, 0),
         ᶜlinear_buoygrad,
-        ᶜtke⁰,
+        max(ᶜtke⁰, 0),
         obukhov_length,
         ᶜstrain_rate_norm,
         ᶜprandtl_nvec,
@@ -251,8 +249,25 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     turbconv_params = CAP.turbconv_params(params)
     c_m = TCP.tke_ed_coeff(turbconv_params)
     @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0))
-    # TODO: add Prantdl number
     @. ᶜK_h = ᶜK_u / ᶜprandtl_nvec
 
+    ρatke_flux_values = Fields.field_values(ρatke_flux)
+    ρ_int_values = Fields.field_values(Fields.level(ᶜρa⁰, 1))
+    u_int_values = Fields.field_values(Fields.level(ᶜu, 1))
+    ustar_values = Fields.field_values(ustar)
+    int_local_geometry_values =
+        Fields.field_values(Fields.level(Fields.local_geometry_field(Y.c), 1))
+    sfc_local_geometry_values = Fields.field_values(
+        Fields.level(Fields.local_geometry_field(Y.f), half),
+    )
+    @. ρatke_flux_values = surface_flux_tke(
+        turbconv_params,
+        ρ_int_values,
+        u_int_values,
+        ustar_values,
+        int_local_geometry_values,
+        sfc_local_geometry_values,
+    )
+
     return nothing
 end

src/prognostic_equations/edmfx_closures.jl

@@ -194,10 +194,10 @@ function mixing_length(
     # kz scale (surface layer)
     if obukhov_length < 0.0 #unstable
         l_W =
-            vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m) *
+            vkc * (ᶜz - z_sfc) / (max(sqrt(sfc_tke / ustar / ustar), eps(FT)) * c_m) *
             min((1 - 100 * (ᶜz - z_sfc) / obukhov_length)^FT(0.2), 1 / vkc)
     else # neutral or stable
-        l_W = vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m)
+        l_W = vkc * (ᶜz - z_sfc) / (max(sqrt(sfc_tke / ustar / ustar), eps(FT)) * c_m)
     end
 
     # compute l_TKE - the production-dissipation balanced length scale
@@ -235,15 +235,19 @@ function mixing_length(
         l_smag = c_smag * ᶜdz
     end
 
+    l_smag = FT(10)
     # add limiters
-    l = SA.SVector(
-        (l_N < eps(FT) || l_N > l_z) ? l_z : l_N,
-        (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE,
-        (l_W < eps(FT) || l_W > l_z) ? l_z : l_W,
-    )
+    # l = SA.SVector(
+    #     (l_N < eps(FT) || l_N > l_z) ? l_z : l_N,
+    #     (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE,
+    #     (l_W < eps(FT) || l_W > l_z) ? l_z : l_W,
+    # )
+    # l_limited = lamb_smooth_minimum(l, smin_ub, smin_rm)
+    l = SA.SVector(l_N, l_TKE, l_W)
+    l_limited = max(l_smag, min(lamb_smooth_minimum(l, smin_ub, smin_rm), l_z))
     # get soft minimum
     # TODO: limit it with l_smag
-    return lamb_smooth_minimum(l, smin_ub, smin_rm)
+    return l_limited
 end
 
 """

src/prognostic_equations/edmfx_entr_detr.jl

@@ -106,7 +106,10 @@ function entrainment(
     ::ConstantCoefficientEntrainment,
 ) where {FT}
     entr_coeff = CAP.entr_coeff(params)
-    entr = min(entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc), 1 / dt)
+    min_area_limiter = FT(2e-3) * exp(-12 * ᶜaʲ)
+    #entr = min(entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc), 1 / dt)
+    entr =
+        min(abs(ᶜwʲ) * (entr_coeff / (ᶜz - z_sfc) + min_area_limiter), 1 / dt)
     return entr
 end
 
@@ -258,11 +261,13 @@ function detrainment(
 ) where {FT}
     detr_coeff = CAP.detr_coeff(params)
     detr_buoy_coeff = CAP.detr_buoy_coeff(params)
+    max_area_limiter = FT(0.1) * exp(-10 * (0.6 - ᶜaʲ))
     detr = min(
         abs(ᶜwʲ) * (
             detr_coeff +
             detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-            max(eps(FT), (ᶜwʲ - ᶜw⁰) * (ᶜwʲ - ᶜw⁰))
+            max(eps(FT), (ᶜwʲ - ᶜw⁰) * (ᶜwʲ - ᶜw⁰)) +
+            max_area_limiter
         ),
         1 / dt,
     )

src/prognostic_equations/edmfx_sgs_flux.jl

@@ -169,7 +169,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     if p.atmos.edmfx_sgs_diffusive_flux
         # energy
         ᶠρaK_h = p.ᶠtemp_scalar
-        @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_h[colidx])
+        @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_h[colidx]), FT(0))
 
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
@@ -192,7 +192,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
 
         # momentum
         ᶠρaK_u = p.ᶠtemp_scalar
-        @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_u[colidx])
+        @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_u[colidx]), FT(0))
         ᶠstrain_rate = p.ᶠtemp_UVWxUVW
         compute_strain_rate_face!(ᶠstrain_rate[colidx], ᶜu⁰[colidx])
         @. Yₜ.c.uₕ[colidx] -= C12(

src/prognostic_equations/edmfx_tke.jl

@@ -10,7 +10,7 @@ function edmfx_tke_tendency!(
     p,
     t,
     colidx,
-    turbconv_model::DiagnosticEDMFX,
+    turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
 
     FT = Spaces.undertype(axes(Y.c))
@@ -21,12 +21,13 @@ function edmfx_tke_tendency!(
     (; ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p
     (; ᶠu³⁰, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰, ᶜmixing_length) = p
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p
+    ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : Y.c.ρ
     ᶠgradᵥ = Operators.GradientC2F()
 
     ᶠρaK_u = p.ᶠtemp_scalar
     if use_prognostic_tke(turbconv_model)
         # turbulent transport (diffusive flux)
-        @. ᶠρaK_u[colidx] = ᶠinterp(Y.c.ρ[colidx]) * ᶠinterp(ᶜK_u[colidx])
+        @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_u[colidx])
         # boundary condition for the diffusive flux
         ᶜdivᵥ_ρatke = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
@@ -36,15 +37,15 @@ function edmfx_tke_tendency!(
             ᶜdivᵥ_ρatke(-(ᶠρaK_u[colidx] * ᶠgradᵥ(ᶜtke⁰[colidx])))
         # shear production
         @. Yₜ.c.sgs⁰.ρatke[colidx] +=
-            2 * Y.c.ρ[colidx] * ᶜK_u[colidx] * ᶜstrain_rate_norm[colidx]
+            2 * ᶜρa⁰[colidx] * ᶜK_u[colidx] * ᶜstrain_rate_norm[colidx]
         # buoyancy production
         @. Yₜ.c.sgs⁰.ρatke[colidx] -=
-            Y.c.ρ[colidx] * ᶜK_h[colidx] * ᶜlinear_buoygrad[colidx]
+            ᶜρa⁰[colidx] * ᶜK_h[colidx] * ᶜlinear_buoygrad[colidx]
         # entrainment and detraiment
         # using ᶜu⁰ and local geometry results in allocation
         for j in 1:n
             @. Yₜ.c.sgs⁰.ρatke[colidx] +=
-                Y.c.ρ[colidx] * (
+                ᶜρa⁰[colidx] * (
                     ᶜentrʲs.:($$j)[colidx] * 1 / 2 * norm_sqr(
                         ᶜinterp(ᶠu³⁰[colidx]) - ᶜinterp(ᶠu³ʲs.:($$j)[colidx]),
                     ) - ᶜdetrʲs.:($$j)[colidx] * ᶜtke⁰[colidx]
@@ -53,7 +54,7 @@ function edmfx_tke_tendency!(
         # pressure work
         # dissipation
         @. Yₜ.c.sgs⁰.ρatke[colidx] -=
-            Y.c.ρ[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) /
+            ᶜρa⁰[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) /
             ᶜmixing_length[colidx]
     end