Merge pull request #4088 from CliMA/zs/rm_tke_function

remove tke_exchange and specific_tke

Author: szy21

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-273
+274
 
 # **README**
 #
@@ -20,6 +20,9 @@
 
 
 #=
+274
+- Remove unused calculation of TKE exchange in mixing length
+
 273
 - Updated deps
 

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -505,7 +505,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
         ),
     )
 
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, Y.c.ρ, turbconv_model))
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
 
     for i in 2:Spaces.nlevels(axes(Y.c))
         ρ_level = Fields.field_values(Fields.level(Y.c.ρ, i))

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -24,9 +24,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!(
     (; ᶜp, ᶜK) = p.precomputed
     (; ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰) = p.precomputed
 
-    ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
-
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
     set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
     set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
     # @. ᶜK⁰ += ᶜtke⁰
@@ -346,7 +344,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
     ᶜlg = Fields.local_geometry_field(Y.c)
     ᶠlg = Fields.local_geometry_field(Y.f)
     ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
 
     ᶜvert_div = p.scratch.ᶜtemp_scalar
     ᶜmassflux_vert_div = p.scratch.ᶜtemp_scalar_2

src/diagnostics/Diagnostics.jl

@@ -67,7 +67,7 @@ import ..turbulent_prandtl_number
 import ..smagorinsky_lilly_length
 import ..ᶜcompute_eddy_diffusivity_coefficient
 import ..ρa⁰
-import ..specific_tke
+import ..specific
 import ..ᶜspecific_env_value
 
 

src/diagnostics/edmfx_diagnostics.jl

@@ -1398,7 +1398,7 @@ function compute_tke!(
     end
 
     ᶜtke = @. lazy(
-        specific_tke(state.c.ρ, state.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model),
+        specific(state.c.sgs⁰.ρatke, state.c.ρ),
     )
     if isnothing(out)
         return Base.materialize(ᶜtke)
@@ -1574,7 +1574,7 @@ function compute_edt!(
         turbconv_model isa PrognosticEDMFX ?
         (@. lazy(ρa⁰(state.c.ρ, state.c.sgsʲs, turbconv_model))) : state.c.ρ
     ᶜtke⁰ = @. lazy(
-        specific_tke(state.c.ρ, state.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model),
+        specific(state.c.sgs⁰.ρatke, state.c.ρ),
     )
     ᶜmixing_length_field = ᶜmixing_length(state, cache)
     ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field))
@@ -1663,7 +1663,7 @@ function compute_evu!(
         turbconv_model isa PrognosticEDMFX ?
         (@. lazy(ρa⁰(state.c.ρ, state.c.sgsʲs, turbconv_model))) : state.c.ρ
     ᶜtke⁰ = @. lazy(
-        specific_tke(state.c.ρ, state.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model),
+        specific(state.c.sgs⁰.ρatke, state.c.ρ),
     )
     ᶜmixing_length_field = ᶜmixing_length(state, cache)
     ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field))

src/prognostic_equations/advection.jl

@@ -231,7 +231,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     end
     ᶜtke⁰ =
         advect_tke ?
-        (@. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))) :
+        (@. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))) :
         nothing
     ᶜa_scalar = p.scratch.ᶜtemp_scalar
     ᶜω³ = p.scratch.ᶜtemp_CT3

src/prognostic_equations/eddy_diffusion_closures.jl

@@ -512,19 +512,13 @@ function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
     ᶜdz = Fields.Δz_field(axes(Y.c))
 
-    turbconv_model = p.atmos.turbconv_model
-    ᶜρa⁰ =
-        turbconv_model isa PrognosticEDMFX ?
-        (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
     sfc_tke = Fields.level(ᶜtke⁰, 1)
 
     ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar_5
     @. ᶜprandtl_nvec =
         turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
 
-    ᶜtke_exch = ᶜtke_exchange(Y, p)
-
     ᶜmixing_length_tuple = @. lazy(
         mixing_length_lopez_gomez_2020(
             params,
@@ -629,62 +623,6 @@ function turbulent_prandtl_number(params, ᶜN²_eff, ᶜstrain_rate_norm)
     return min(max(prandtl_nvec, eps_FT), Pr_max)
 end
 
-
-"""
-    ᶜtke_exchange(Y, p)
-
-Calculates the turbulent kinetic energy (TKE) exchange tendency between the
-environment and updrafts due to detrainment.
-
-Arguments:
-- `Y`: The prognostic state vector.
-- `p`: Cache
-
-Returns:
-- The TKE exchange tendency term [m²/s³].
-"""
-function ᶜtke_exchange(Y, p)
-    (; turbconv_model) = p.atmos
-    n = n_mass_flux_subdomains(turbconv_model)
-    ᶜρa⁰ =
-        p.atmos.turbconv_model isa PrognosticEDMFX ?
-        (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
-
-    if p.atmos.turbconv_model isa PrognosticEDMFX
-        (; ᶜdetrʲs, ᶠu³⁰, ᶠu³ʲs) = p.precomputed
-        ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
-        @. ᶜtke_exch = 0
-        for j in 1:n
-            @. ᶜtke_exch +=
-                Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ * (
-                    1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) -
-                    ᶜtke⁰
-                )
-        end
-
-        return ᶜtke_exch
-    elseif p.atmos.turbconv_model isa DiagnosticEDMFX
-        (; ᶜdetrʲs, ᶠu³⁰, ᶠu³ʲs, ᶜρaʲs) = p.precomputed
-        ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
-        @. ᶜtke_exch = 0
-        for j in 1:n
-            @. ᶜtke_exch +=
-                ᶜρaʲs.:($$j) * ᶜdetrʲs.:($$j) / ᶜρa⁰ * (
-                    1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) -
-                    ᶜtke⁰
-                )
-        end
-
-        return ᶜtke_exch
-        # ED only or none-EDMF model does not have updrafts (or detrainment),
-        # so tke exchange is 0
-    else
-        return 0
-    end
-
-end
-
 """
     blend_scales(
         method::AbstractScaleBlending,

src/prognostic_equations/edmfx_sgs_flux.jl

@@ -394,7 +394,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
     ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
 
     if p.atmos.edmfx_model.sgs_diffusive_flux isa Val{true}
 
@@ -535,7 +535,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; ᶜu, ᶜts) = p.precomputed
     (; ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, Y.c.ρ, turbconv_model))
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
 
     if p.atmos.edmfx_model.sgs_diffusive_flux isa Val{true}
 

src/prognostic_equations/edmfx_tke.jl

@@ -29,8 +29,7 @@ function edmfx_tke_tendency!(Yₜ, Y, p, t, turbconv_model::EDOnlyEDMFX)
     (; params) = p
     (; ᶜstrain_rate_norm, ᶜlinear_buoygrad) = p.precomputed
     turbconv_params = CAP.turbconv_params(p.params)
-    ᶜρa⁰ = Y.c.ρ # EDOnly
-    ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
+    ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
     ᶜmixing_length_field = p.scratch.ᶜtemp_scalar
     ᶜmixing_length_field .= ᶜmixing_length(Y, p)
     ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field))
@@ -57,20 +56,13 @@ function edmfx_tke_tendency!(
     turbconv_params = CAP.turbconv_params(p.params)
     FT = eltype(p.params)
     thermo_params = CAP.thermodynamics_params(p.params)
-    ᶜρa⁰ =
-        turbconv_model isa PrognosticEDMFX ?
-        (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ
 
     if use_prognostic_tke(turbconv_model)
         (; ᶜρʲs) = p.precomputed
         ᶠz = Fields.coordinate_field(Y.f).z
         ᶜmixing_length_field = p.scratch.ᶜtemp_scalar_2
         ᶜmixing_length_field .= ᶜmixing_length(Y, p)
-        ᶜρa⁰ =
-            turbconv_model isa PrognosticEDMFX ?
-            (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ
-        ᶜtke⁰ =
-            @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
+        ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
         ᶜK_u = @. lazy(
             eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field),
         )
@@ -98,6 +90,7 @@ function edmfx_tke_tendency!(
         end
         # Note: Adding the following tendency breaks bm_aquaplanet_progedmf_dense_autodiff
         if turbconv_model isa PrognosticEDMFX
+            ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
             (; ᶜts⁰) = p.precomputed
             ᶜρ⁰ = @. lazy(TD.air_density(thermo_params, ᶜts⁰))
             @. Yₜ.c.sgs⁰.ρatke -=

src/prognostic_equations/hyperdiffusion.jl

@@ -114,10 +114,7 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     @. ᶜ∇²specific_energy = wdivₕ(gradₕ(specific(Y.c.ρe_tot, Y.c.ρ) + ᶜp / Y.c.ρ - ᶜh_ref))
 
     if diffuse_tke
-        ᶜρa⁰ =
-            turbconv_model isa PrognosticEDMFX ?
-            (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ
-        ᶜtke⁰ = @. lazy(specific_tke(Y.c.ρ, Y.c.sgs⁰.ρatke, ᶜρa⁰, turbconv_model))
+        ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
         (; ᶜ∇²tke⁰) = p.hyperdiff
         @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜtke⁰))
     end