diff --git a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml
index c846de7cce..4cdf8b106d 100644
--- a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml
+++ b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml
@@ -41,7 +41,7 @@ diagnostics:
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
-  - short_name: [entr, detr, lmix, bgrad, strain, edt, evu]
+  - short_name: [entr, detr, lmix, bgrad, strain,] # edt, evu]
     period: 10mins
   - short_name: [rlut, rlutcs, rsut, rsutcs, clwvi, lwp, clivi, dsevi, clvi, prw, hurvi, husv]
     period: 10mins
diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl
index 42a4188d0a..a075adac1a 100644
--- a/src/cache/cloud_fraction.jl
+++ b/src/cache/cloud_fraction.jl
@@ -39,9 +39,23 @@ NVTX.@annotate function set_cloud_fraction!(
 )
     (; params) = p
     (; turbconv_model) = p.atmos
-    (; ᶜts, ᶜmixing_length, cloud_diagnostics_tuple) = p.precomputed
+    (; ᶜts, cloud_diagnostics_tuple) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
+    FT = eltype(p.params)
+    n = n_mass_flux_subdomains(turbconv_model)
     if isnothing(turbconv_model)
+        (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+        (;
+            ᶜlinear_buoygrad,
+            ᶜstrain_rate_norm,
+            ᶠu³⁰,
+            ᶠu³,
+            ᶜtke⁰,
+            ᶜentrʲs,
+            ᶜdetrʲs,
+            ᶠu³ʲs,
+            ᶜρa⁰,
+        ) = p.precomputed
         if p.atmos.call_cloud_diagnostics_per_stage isa
            CallCloudDiagnosticsPerStage
             (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
@@ -53,6 +67,39 @@ NVTX.@annotate function set_cloud_fraction!(
             @. ᶜgradᵥ_θ_liq_ice =
                 ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
         end
+        ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
+        @. ᶜprandtl_nvec =
+            turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
+        sfc_tke = Fields.level(ᶜtke⁰, 1)
+        z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+        ᶜz = Fields.coordinate_field(Y.c).z
+        ᶜdz = Fields.Δz_field(axes(Y.c))
+
+        ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
+        @. ᶜtke_exch = 0
+        for j in 1:n
+            ᶠu³ʲ = ᶠu³ʲs.:($j)
+            @. ᶜtke_exch +=
+                Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ *
+                (1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰)
+        end
+
+        ᶜmixing_length = @. lazy(master_mixing_length(
+            params,
+            ustar,
+            ᶜz,
+            z_sfc,
+            ᶜdz,
+            max(sfc_tke, eps(FT)),
+            ᶜlinear_buoygrad,
+            max(ᶜtke⁰, 0),
+            obukhov_length,
+            ᶜstrain_rate_norm,
+            ᶜprandtl_nvec,
+            ᶜtke_exch,
+            p.atmos.edmfx_model.scale_blending_method,
+        ))
         compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     end
     if moist_model isa EquilMoistModel
@@ -64,12 +111,12 @@ NVTX.@annotate function set_cloud_fraction!(
     else
         @. cloud_diagnostics_tuple = make_named_tuple(
             ifelse(
-                p.precomputed.ᶜspecific.q_liq + p.precomputed.ᶜspecific.q_ice > 0,
+                specific(Y.c.ρq_liq, Y.c.ρ) + specific(Y.c.ρq_ice, Y.c.ρ) > 0,
                 1,
                 0,
             ),
-            p.precomputed.ᶜspecific.q_liq,
-            p.precomputed.ᶜspecific.q_ice,
+            specific(Y.c.ρq_liq, Y.c.ρ),
+            specific(Y.c.ρq_ice, Y.c.ρ),
         )
     end
 end
@@ -81,13 +128,26 @@ NVTX.@annotate function set_cloud_fraction!(
 )
     SG_quad = qc.SG_quad
     (; params) = p
+    (; ustar, obukhov_length) = p.precomputed.sfc_conditions
 
     FT = eltype(params)
+    n = n_mass_flux_subdomains(turbconv_model)
     thermo_params = CAP.thermodynamics_params(params)
-    (; ᶜts, ᶜmixing_length, cloud_diagnostics_tuple) = p.precomputed
+    (; ᶜts, cloud_diagnostics_tuple) = p.precomputed
     (; turbconv_model) = p.atmos
 
     if isnothing(turbconv_model)
+        (;
+            ᶜlinear_buoygrad,
+            ᶜstrain_rate_norm,
+            ᶠu³⁰,
+            ᶠu³,
+            ᶜtke⁰,
+            ᶜentrʲs, 
+            ᶜdetrʲs, 
+            ᶠu³ʲs,
+            ᶜρa⁰,
+        ) = p.precomputed
         if p.atmos.call_cloud_diagnostics_per_stage isa
            CallCloudDiagnosticsPerStage
             (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
@@ -98,7 +158,47 @@ NVTX.@annotate function set_cloud_fraction!(
                 ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts)))
             @. ᶜgradᵥ_θ_liq_ice =
                 ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
+
+            @. ᶜprandtl_nvec = turbulent_prandtl_number(
+                params,
+                ᶜlinear_buoygrad,
+                ᶜstrain_rate_norm,
+            )
         end
+
+        ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
+        @. ᶜprandtl_nvec =
+            turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
+        sfc_tke = Fields.level(ᶜtke⁰, 1)
+        z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+        ᶜz = Fields.coordinate_field(Y.c).z
+        ᶜdz = Fields.Δz_field(axes(Y.c))
+
+        ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
+        @. ᶜtke_exch = 0
+        for j in 1:n
+            ᶠu³ʲ = ᶠu³ʲs.:($j)
+            @. ᶜtke_exch +=
+                Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ *
+                (1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰)
+        end
+
+        ᶜmixing_length = @. lazy(master_mixing_length(
+            params,
+            ustar,
+            ᶜz,
+            z_sfc,
+            ᶜdz,
+            max(sfc_tke, eps(FT)),
+            ᶜlinear_buoygrad,
+            max(ᶜtke⁰, 0),
+            obukhov_length,
+            ᶜstrain_rate_norm,
+            ᶜprandtl_nvec,
+            ᶜtke_exch,
+            p.atmos.edmfx_model.scale_blending_method,
+        ))
         compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     end
 
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index a34872b231..e4601cf8c8 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -61,7 +61,7 @@ NVTX.@annotate function set_diagnostic_edmfx_env_quantities_level!(
     local_geometry_halflevel,
     turbconv_model,
 )
-    @. u³⁰_halflevel = divide_by_ρa(
+    @. u³⁰_halflevel = specific(
         ρ_level * u³_halflevel -
         unrolled_dotproduct(ρaʲs_level, u³ʲs_halflevel),
         ρ_level,
@@ -950,7 +950,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     (; params) = p
     (; dt) = p
     (; ᶜp, ᶜu, ᶜts) = p.precomputed
-    (; q_tot) = p.precomputed.ᶜspecific
     (; ustar, obukhov_length) = p.precomputed.sfc_conditions
     (; ᶜtke⁰) = p.precomputed
     (;
@@ -1024,8 +1023,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     )
     @. ᶜmixing_length = ᶜmixing_length_tuple.master
 
-    @. ᶜK_u = eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length)
-    @. ᶜK_h = eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec)
+    ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length))
+    ᶜK_h = @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
 
     ρatke_flux_values = Fields.field_values(ρatke_flux)
     ρ_int_values = Fields.field_values(Fields.level(Y.c.ρ, 1))
@@ -1068,14 +1067,13 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
     microphys_0m_params = CAP.microphysics_0m_params(p.params)
     (; dt) = p
     (; ᶜts, ᶜSqₜᵖ⁰) = p.precomputed
-    (; q_tot) = p.precomputed.ᶜspecific
 
     # Environment precipitation sources (to be applied to grid mean)
     @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources(
         thermo_params,
         microphys_0m_params,
         dt,
-        q_tot,
+        specific(Y.c.ρq_tot, Y.c.ρ),
         ᶜts,
     )
     return nothing
diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl
index ad462c8473..5168512681 100644
--- a/src/cache/precipitation_precomputed_quantities.jl
+++ b/src/cache/precipitation_precomputed_quantities.jl
@@ -105,7 +105,6 @@ function set_precipitation_velocities!(
     precip_model::Microphysics2Moment,
 )
     (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
-    (; q_liq, q_ice, q_rai, q_sno) = p.precomputed.ᶜspecific
     (; ᶜΦ) = p.core
 
     cm1c = CAP.microphysics_cloud_params(p.params)
@@ -116,23 +115,23 @@ function set_precipitation_velocities!(
     # compute the precipitation terminal velocity [m/s]
     # TODO sedimentation of snow is based on the 1M scheme
     @. ᶜwnᵣ = getindex(
-        CM2.rain_terminal_velocity(cm2p.sb, cm2p.tv, q_rai, Y.c.ρ, Y.c.ρn_rai),
+        CM2.rain_terminal_velocity(cm2p.sb, cm2p.tv, specific(Y.c.ρq_rai, Y.c.ρ), Y.c.ρ, Y.c.ρn_rai),
         1,
     )
     @. ᶜwᵣ = getindex(
-        CM2.rain_terminal_velocity(cm2p.sb, cm2p.tv, q_rai, Y.c.ρ, Y.c.ρn_rai),
+        CM2.rain_terminal_velocity(cm2p.sb, cm2p.tv, specific(Y.c.ρq_rai, Y.c.ρ), Y.c.ρ, Y.c.ρn_rai),
         2,
     )
-    @. ᶜwₛ = CM1.terminal_velocity(cm1p.ps, cm1p.tv.snow, Y.c.ρ, q_sno)
+    @. ᶜwₛ = CM1.terminal_velocity(cm1p.ps, cm1p.tv.snow, Y.c.ρ, specific(Y.c.ρq_sno, Y.c.ρ))
     # compute sedimentation velocity for cloud condensate [m/s]
     # TODO sedimentation velocities of cloud condensates are based 
     # on the 1M scheme. Sedimentation velocity of cloud number concentration
     # is equal to that of the mass.
     @. ᶜwnₗ =
-        CMNe.terminal_velocity(cm1c.liquid, cm1c.Ch2022.rain, Y.c.ρ, q_liq)
-    @. ᶜwₗ = CMNe.terminal_velocity(cm1c.liquid, cm1c.Ch2022.rain, Y.c.ρ, q_liq)
+        CMNe.terminal_velocity(cm1c.liquid, cm1c.Ch2022.rain, Y.c.ρ, specific(Y.c.ρq_liq, Y.c.ρ))
+    @. ᶜwₗ = CMNe.terminal_velocity(cm1c.liquid, cm1c.Ch2022.rain, Y.c.ρ, specific(Y.c.ρq_liq, Y.c.ρ))
     @. ᶜwᵢ =
-        CMNe.terminal_velocity(cm1c.ice, cm1c.Ch2022.small_ice, Y.c.ρ, q_ice)
+        CMNe.terminal_velocity(cm1c.ice, cm1c.Ch2022.small_ice, Y.c.ρ, specific(Y.c.ρq_ice, Y.c.ρ))
 
     # compute their contributions to energy and total water advection
     @. ᶜwₜqₜ =
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index ccbe7136ca..6c5093b1ae 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -172,9 +172,7 @@ function precomputed_quantities(Y, atmos)
     advective_sgs_quantities =
         atmos.turbconv_model isa PrognosticEDMFX ?
         (;
-            ᶜmixing_length_tuple = similar(Y.c, MixingLength{FT}),
-            ᶜK_u = similar(Y.c, FT),
-            ᶜK_h = similar(Y.c, FT),
+            # ᶜmixing_length_tuple = similar(Y.c, MixingLength{FT}),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
             bdmr_l = similar(Y.c, BidiagonalMatrixRow{FT}),
             bdmr_r = similar(Y.c, BidiagonalMatrixRow{FT}),
@@ -194,7 +192,6 @@ function precomputed_quantities(Y, atmos)
         (;
             ᶜmixing_length_tuple = similar(Y.c, MixingLength{FT}),
             ᶜtke⁰ = similar(Y.c, FT),
-            ᶜK_u = similar(Y.c, FT),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
             ᶜK_h = similar(Y.c, FT),
         ) : (;)
@@ -398,34 +395,34 @@ function thermo_state(
     return get_ts(ρ, p, θ, e_int, q_tot, q_pt)
 end
 
-function thermo_vars(moisture_model, precip_model, specific, K, Φ)
-    energy_var = (; e_int = specific.e_tot - K - Φ)
+function thermo_vars(moisture_model, precip_model, ᶜY, K, Φ)
+    energy_var = (; e_int = specific(ᶜY.ρe_tot, ᶜY.ρ) - K - Φ)
     moisture_var = if moisture_model isa DryModel
         (;)
     elseif moisture_model isa EquilMoistModel
-        (; specific.q_tot)
+        (; q_tot = specific(ᶜY.ρq_tot, ᶜY.ρ))
     elseif moisture_model isa NonEquilMoistModel
         q_pt_args = (
-            specific.q_tot,
-            specific.q_liq + specific.q_rai,
-            specific.q_ice + specific.q_sno,
+            q_tot = specific(ᶜY.ρq_tot, ᶜY.ρ),
+            q_liq = specific(ᶜY.ρq_liq, ᶜY.ρ),
+            q_ice = specific(ᶜY.ρq_ice, ᶜY.ρ),
         )
         (; q_pt = TD.PhasePartition(q_pt_args...))
     end
     return (; energy_var..., moisture_var...)
 end
 
-ts_gs(thermo_params, moisture_model, precip_model, specific, K, Φ, ρ) =
+ts_gs(thermo_params, moisture_model, precip_model, ᶜY, K, Φ, ρ) =
     thermo_state(
         thermo_params;
-        thermo_vars(moisture_model, precip_model, specific, K, Φ)...,
+        thermo_vars(moisture_model, precip_model, ᶜY, K, Φ)...,
         ρ,
     )
 
-ts_sgs(thermo_params, moisture_model, precip_model, specific, K, Φ, p) =
+ts_sgs(thermo_params, moisture_model, precip_model, ᶜY, K, Φ, p) =
     thermo_state(
         thermo_params;
-        thermo_vars(moisture_model, precip_model, specific, K, Φ)...,
+        thermo_vars(moisture_model, precip_model, ᶜY, K, Φ)...,
         p,
     )
 
@@ -489,9 +486,9 @@ NVTX.@annotate function set_implicit_precomputed_quantities!(Y, p, t)
         # @. ᶜK += Y.c.sgs⁰.ρatke / Y.c.ρ
         # TODO: We should think more about these increments before we use them.
     end
-    @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, ᶜΦ, Y.c.ρ)
+    @. ᶜts = ts_gs(thermo_args..., Y.c, ᶜK, ᶜΦ, Y.c.ρ)
     @. ᶜp = TD.air_pressure(thermo_params, ᶜts)
-    @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot)
+    @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, specific(Y.c.ρe_tot, Y.c.ρ))
 
     if turbconv_model isa PrognosticEDMFX
         set_prognostic_edmf_precomputed_quantities_draft!(Y, p, ᶠuₕ³, t)
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index 3aa32fb8db..023e862409 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -30,15 +30,15 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!(
     end
 
     @. ᶜρa⁰ = ρa⁰(Y.c)
-    @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
-    @. ᶜmse⁰ = divide_by_ρa(
+    @. ᶜtke⁰ = specific(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
+    @. ᶜmse⁰ = specific(
         Y.c.ρ * (ᶜh_tot - ᶜK) - ρamse⁺(Y.c.sgsʲs),
         ᶜρa⁰,
         Y.c.ρ * (ᶜh_tot - ᶜK),
         Y.c.ρ,
         turbconv_model,
     )
-    @. ᶜq_tot⁰ = divide_by_ρa(
+    @. ᶜq_tot⁰ = specific(
         Y.c.ρq_tot - ρaq_tot⁺(Y.c.sgsʲs),
         ᶜρa⁰,
         Y.c.ρq_tot,
@@ -47,28 +47,28 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!(
     )
     if p.atmos.moisture_model isa NonEquilMoistModel &&
        p.atmos.precip_model isa Microphysics1Moment
-        @. ᶜq_liq⁰ = divide_by_ρa(
+        @. ᶜq_liq⁰ = specific(
             Y.c.ρq_liq - ρaq_liq⁺(Y.c.sgsʲs),
             ᶜρa⁰,
             Y.c.ρq_liq,
             Y.c.ρ,
             turbconv_model,
         )
-        @. ᶜq_ice⁰ = divide_by_ρa(
+        @. ᶜq_ice⁰ = specific(
             Y.c.ρq_ice - ρaq_ice⁺(Y.c.sgsʲs),
             ᶜρa⁰,
             Y.c.ρq_ice,
             Y.c.ρ,
             turbconv_model,
         )
-        @. ᶜq_rai⁰ = divide_by_ρa(
+        @. ᶜq_rai⁰ = specific(
             Y.c.ρq_rai - ρaq_rai⁺(Y.c.sgsʲs),
             ᶜρa⁰,
             Y.c.ρq_rai,
             Y.c.ρ,
             turbconv_model,
         )
-        @. ᶜq_sno⁰ = divide_by_ρa(
+        @. ᶜq_sno⁰ = specific(
             Y.c.ρq_sno - ρaq_sno⁺(Y.c.sgsʲs),
             ᶜρa⁰,
             Y.c.ρq_sno,
@@ -369,12 +369,10 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
 
     (; ᶜtke⁰, ᶜu, ᶜp, ᶜρa⁰, ᶠu³⁰, ᶜts⁰, ᶜq_tot⁰) = p.precomputed
     (;
-        ᶜmixing_length_tuple,
-        ᶜmixing_length,
+        # ᶜmixing_length_tuple,
+        # ᶜmixing_length,
         ᶜlinear_buoygrad,
         ᶜstrain_rate_norm,
-        ᶜK_u,
-        ᶜK_h,
         ρatke_flux,
     ) = p.precomputed
     (;
@@ -508,33 +506,28 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
     ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
     @. ᶜtke_exch = 0
     for j in 1:n
-        ᶠu³ʲ = ᶠu³ʲs.:($j)
         @. ᶜtke_exch +=
             Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ *
             (1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰)
     end
 
     sfc_tke = Fields.level(ᶜtke⁰, 1)
-    @. ᶜmixing_length_tuple = mixing_length(
-        p.params,
-        ustar,
-        ᶜz,
-        z_sfc,
-        ᶜdz,
-        max(sfc_tke, eps(FT)),
-        ᶜlinear_buoygrad,
-        max(ᶜtke⁰, 0),
-        obukhov_length,
-        ᶜstrain_rate_norm,
-        ᶜprandtl_nvec,
-        ᶜtke_exch,
-        p.atmos.edmfx_model.scale_blending_method,
-    )
-
-    @. ᶜmixing_length = ᶜmixing_length_tuple.master
-
-    @. ᶜK_u = eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length)
-    @. ᶜK_h = eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec)
+ 
+    # ᶜmixing_length = @. lazy(master_mixing_length(
+    #         p.params,
+    #         ustar,
+    #         ᶜz,
+    #         z_sfc,
+    #         ᶜdz,
+    #         max(sfc_tke, eps(FT)),
+    #         ᶜlinear_buoygrad,
+    #         max(ᶜtke⁰, 0),
+    #         obukhov_length,
+    #         ᶜstrain_rate_norm,
+    #         ᶜprandtl_nvec,
+    #         ᶜtke_exch,
+    #         p.atmos.edmfx_model.scale_blending_method,
+    #     ))
 
     ρatke_flux_values = Fields.field_values(ρatke_flux)
     ρa_sfc_values = Fields.field_values(Fields.level(ᶜρa⁰, 1)) # TODO: replace by surface value
diff --git a/src/diagnostics/edmfx_diagnostics.jl b/src/diagnostics/edmfx_diagnostics.jl
index ec6aab3b8c..51a45cf164 100644
--- a/src/diagnostics/edmfx_diagnostics.jl
+++ b/src/diagnostics/edmfx_diagnostics.jl
@@ -1287,10 +1287,20 @@ function compute_edt!(
     vert_diff::Nothing,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+
+    turbconv_params = CAP.turbconv_params(cache.params)   
+    (; ᶜmixing_length_tuple, ᶜtke⁰, ᶜlinear_buoygrad, ᶜstrain_rate_norm) = cache.precomputed  
+
+    ᶜprandtl_nvec = @. lazy(turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm))
+
+    ᶜmixing_length = ᶜmixing_length_tuple.master
+
+    ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length))
+    ᶜK_h = @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
     if isnothing(out)
-        return copy(cache.precomputed.ᶜK_h)
+        return copy(ᶜK_h)
     else
-        out .= cache.precomputed.ᶜK_h
+        out .= ᶜK_h
     end
 end
 
@@ -1342,10 +1352,15 @@ function compute_evu!(
     vert_diff::Nothing,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+    turbconv_params = CAP.turbconv_params(cache.params)   
+    (; ᶜmixing_length_tuple, ᶜtke⁰) = cache.precomputed  
+    ᶜmixing_length = ᶜmixing_length_tuple.master
+    ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length))
+
     if isnothing(out)
-        return copy(cache.precomputed.ᶜK_u)
+        return copy(ᶜK_u)
     else
-        out .= cache.precomputed.ᶜK_u
+        out .= ᶜK_u
     end
 end
 
diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl
index c08c92f49a..d3698a105c 100644
--- a/src/parameterized_tendencies/microphysics/cloud_condensate.jl
+++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl
@@ -29,13 +29,12 @@ function cloud_condensate_tendency!(
 )
     (; ᶜts) = p.precomputed
     (; params, dt) = p
-    (; q_rai, q_sno) = p.precomputed.ᶜspecific
     FT = eltype(params)
     thp = CAP.thermodynamics_params(params)
     cmc = CAP.microphysics_cloud_params(params)
 
-    @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, q_rai, dt)
-    @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, q_sno, dt)
+    @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, specific(Y.c.ρq_rai, Y.c.ρ), dt)
+    @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, specific(Y.c.ρq_sno, Y.c.ρ), dt)
 end
 
 function cloud_condensate_tendency!(
@@ -47,21 +46,20 @@ function cloud_condensate_tendency!(
 )
     (; ᶜts) = p.precomputed
     (; params, dt) = p
-    (; q_rai, q_sno, n_liq) = p.precomputed.ᶜspecific
     FT = eltype(params)
     thp = CAP.thermodynamics_params(params)
     cmc = CAP.microphysics_cloud_params(params)
 
-    @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, q_rai, dt)
-    @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, q_sno, dt)
+    @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, specific(Y.c.ρq_rai, Y.c.ρ), dt)
+    @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, specific(Y.c.ρq_sno, Y.c.ρ), dt)
 
     @. Yₜ.c.ρn_liq +=
         Y.c.ρ * aerosol_activation_sources(
             cmc.liquid,
             thp,
             ᶜts,
-            q_rai,
-            n_liq,
+            specific(Y.c.ρq_rai, Y.c.ρ),
+            specific(Y.c.ρn_liq, Y.c.ρ),
             cmc.N_cloud_liquid_droplets,
             dt,
         )
diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl
index dbc2079a2f..c5f3721020 100644
--- a/src/parameterized_tendencies/radiation/radiation.jl
+++ b/src/parameterized_tendencies/radiation/radiation.jl
@@ -440,7 +440,7 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS)
     @assert !(p.atmos.moisture_model isa DryModel)
 
     (; params) = p
-    (; ᶜspecific, ᶜts) = p.precomputed
+    (; ᶜts) = p.precomputed
     (; ᶜκρq, ∫_0_∞_κρq, ᶠ∫_0_z_κρq, isoline_z_ρ_q, ᶠradiation_flux) =
         p.radiation
     thermo_params = CAP.thermodynamics_params(params)
@@ -467,10 +467,10 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS)
     q_tot_isoline = FT(0.008)
     Operators.column_reduce!(
         (nt1, nt2) ->
-            abs(nt1.q_tot - q_tot_isoline) < abs(nt2.q_tot - q_tot_isoline) ?
+            abs(nt1.ρq_tot / nt1.ρ - q_tot_isoline) < abs(nt2.ρq_tot / nt2.ρ - q_tot_isoline) ?
             nt1 : nt2,
         isoline_z_ρ_q,
-        Base.broadcasted(NT ∘ tuple, ᶜz, Y.c.ρ, ᶜspecific.q_tot),
+        Base.broadcasted(NT ∘ tuple, ᶜz, Y.c.ρ, Y.c.ρq_tot),
     )
 
     zi = isoline_z_ρ_q.z
diff --git a/src/prognostic_equations/eddy_diffusion_closures.jl b/src/prognostic_equations/eddy_diffusion_closures.jl
index 7d87bdf8ef..bc7479743a 100644
--- a/src/prognostic_equations/eddy_diffusion_closures.jl
+++ b/src/prognostic_equations/eddy_diffusion_closures.jl
@@ -528,6 +528,41 @@ function mixing_length(
     return MixingLength{FT}(l_final, l_W, l_TKE, l_N, l_grid)
 end
 
+function master_mixing_length(    
+    params,
+    ustar,
+    ᶜz,
+    z_sfc,
+    ᶜdz,
+    sfc_tke,
+    ᶜN²_eff,
+    ᶜtke,
+    obukhov_length,
+    ᶜstrain_rate_norm,
+    ᶜPr,
+    ᶜtke_exch,
+    scale_blending_method,
+)
+    FT = eltype(params)
+    ᶜmixing_length_tuple = mixing_length(
+        params,
+        ustar,
+        ᶜz,
+        z_sfc,
+        ᶜdz,
+        max(sfc_tke, eps(FT)),
+        ᶜN²_eff,
+        max(ᶜtke, 0),
+        obukhov_length,
+        ᶜstrain_rate_norm,
+        ᶜPr,
+        ᶜtke_exch,
+        scale_blending_method,
+    )
+
+    return getproperty(ᶜmixing_length_tuple, :master)
+end
+
 """
     turbulent_prandtl_number(params, ᶜN²_eff, ᶜstrain_rate_norm)
 
@@ -701,8 +736,7 @@ Returns K_u in units of [m^2/s].
 function eddy_viscosity(turbconv_params, tke, mixing_length)
     FT = typeof(tke)
     c_m = CAP.tke_ed_coeff(turbconv_params)
-    K_u = c_m * mixing_length * sqrt(max(tke, FT(0)))
-    return K_u
+    return c_m * mixing_length * sqrt(max(tke, FT(0)))
 end
 
 """
@@ -729,6 +763,5 @@ Returns K_h in units of [m^2/s].
 """
 function eddy_diffusivity(turbconv_params, tke, mixing_length, prandtl_nvec)
     K_u = eddy_viscosity(turbconv_params, tke, mixing_length)
-    K_h = K_u / prandtl_nvec # prandtl_nvec is already bounded by eps_FT and Pr_max
-    return K_h
+    return K_u / prandtl_nvec
 end
diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
index 04ec3169de..a57190cfd0 100644
--- a/src/prognostic_equations/edmfx_sgs_flux.jl
+++ b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -14,7 +14,7 @@ function edmfx_sgs_mass_flux_tendency!(
 
     n = n_mass_flux_subdomains(turbconv_model)
     (; edmfx_sgsflux_upwinding) = p.atmos.numerics
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶠu³, ᶜh_tot) = p.precomputed
     (; ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = p.precomputed
     (; ᶜρa⁰, ᶜρ⁰, ᶠu³⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰) = p.precomputed
     if (
@@ -60,7 +60,7 @@ function edmfx_sgs_mass_flux_tendency!(
             for j in 1:n
                 @. ᶠu³_diff = ᶠu³ʲs.:($$j) - ᶠu³
                 @. ᶜa_scalar =
-                    (Y.c.sgsʲs.:($$j).q_tot - ᶜspecific.q_tot) *
+                    (Y.c.sgsʲs.:($$j).q_tot - specific(Y.c.ρq_tot, Y.c.ρ)) *
                     draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
                 vtt = vertical_transport(
                     ᶜρʲs.:($j),
@@ -72,7 +72,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 @. Yₜ.c.ρq_tot += vtt
             end
             @. ᶠu³_diff = ᶠu³⁰ - ᶠu³
-            @. ᶜa_scalar = (ᶜq_tot⁰ - ᶜspecific.q_tot) * draft_area(ᶜρa⁰, ᶜρ⁰)
+            @. ᶜa_scalar = (ᶜq_tot⁰ - specific(Y.c.ρq_tot, Y.c.ρ)) * draft_area(ᶜρa⁰, ᶜρ⁰)
             vtt = vertical_transport(
                 ᶜρ⁰,
                 ᶠu³_diff,
@@ -91,7 +91,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 @. ᶠu³_diff = ᶠu³ʲs.:($$j) - ᶠu³
 
                 @. ᶜa_scalar =
-                    (Y.c.sgsʲs.:($$j).q_liq - ᶜspecific.q_liq) *
+                    (Y.c.sgsʲs.:($$j).q_liq - specific(Y.c.ρq_liq, Y.c.ρ)) *
                     draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
                 vtt = vertical_transport(
                     ᶜρʲs.:($j),
@@ -103,7 +103,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 @. Yₜ.c.ρq_liq += vtt
 
                 @. ᶜa_scalar =
-                    (Y.c.sgsʲs.:($$j).q_ice - ᶜspecific.q_ice) *
+                    (Y.c.sgsʲs.:($$j).q_ice - specific(Y.c.ρq_ice, Y.c.ρ)) *
                     draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
                 vtt = vertical_transport(
                     ᶜρʲs.:($j),
@@ -115,7 +115,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 @. Yₜ.c.ρq_ice += vtt
 
                 @. ᶜa_scalar =
-                    (Y.c.sgsʲs.:($$j).q_rai - ᶜspecific.q_rai) *
+                    (Y.c.sgsʲs.:($$j).q_rai - specific(Y.c.ρq_rai, Y.c.ρ)) *
                     draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
                 vtt = vertical_transport(
                     ᶜρʲs.:($j),
@@ -127,7 +127,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 @. Yₜ.c.ρq_rai += vtt
 
                 @. ᶜa_scalar =
-                    (Y.c.sgsʲs.:($$j).q_sno - ᶜspecific.q_sno) *
+                    (Y.c.sgsʲs.:($$j).q_sno - specific(Y.c.ρq_sno, Y.c.ρ)) *
                     draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
                 vtt = vertical_transport(
                     ᶜρʲs.:($j),
@@ -140,7 +140,7 @@ function edmfx_sgs_mass_flux_tendency!(
             end
             @. ᶠu³_diff = ᶠu³⁰ - ᶠu³
 
-            @. ᶜa_scalar = (ᶜq_liq⁰ - ᶜspecific.q_liq) * draft_area(ᶜρa⁰, ᶜρ⁰)
+            @. ᶜa_scalar = (ᶜq_liq⁰ - specific(Y.c.ρq_liq, Y.c.ρ)) * draft_area(ᶜρa⁰, ᶜρ⁰)
             vtt = vertical_transport(
                 ᶜρ⁰,
                 ᶠu³_diff,
@@ -150,7 +150,7 @@ function edmfx_sgs_mass_flux_tendency!(
             )
             @. Yₜ.c.ρq_liq += vtt
 
-            @. ᶜa_scalar = (ᶜq_ice⁰ - ᶜspecific.q_ice) * draft_area(ᶜρa⁰, ᶜρ⁰)
+            @. ᶜa_scalar = (ᶜq_ice⁰ - specific(Y.c.ρq_ice, Y.c.ρ)) * draft_area(ᶜρa⁰, ᶜρ⁰)
             vtt = vertical_transport(
                 ᶜρ⁰,
                 ᶠu³_diff,
@@ -160,7 +160,7 @@ function edmfx_sgs_mass_flux_tendency!(
             )
             @. Yₜ.c.ρq_ice += vtt
 
-            @. ᶜa_scalar = (ᶜq_rai⁰ - ᶜspecific.q_rai) * draft_area(ᶜρa⁰, ᶜρ⁰)
+            @. ᶜa_scalar = (ᶜq_rai⁰ - specific(Y.c.ρq_rai, Y.c.ρ)) * draft_area(ᶜρa⁰, ᶜρ⁰)
             vtt = vertical_transport(
                 ᶜρ⁰,
                 ᶠu³_diff,
@@ -170,7 +170,7 @@ function edmfx_sgs_mass_flux_tendency!(
             )
             @. Yₜ.c.ρq_rai += vtt
 
-            @. ᶜa_scalar = (ᶜq_sno⁰ - ᶜspecific.q_sno) * draft_area(ᶜρa⁰, ᶜρ⁰)
+            @. ᶜa_scalar = (ᶜq_sno⁰ - specific(Y.c.ρq_sno, Y.c.ρ)) * draft_area(ᶜρa⁰, ᶜρ⁰)
             vtt = vertical_transport(
                 ᶜρ⁰,
                 ᶠu³_diff,
@@ -198,7 +198,7 @@ function edmfx_sgs_mass_flux_tendency!(
     a_max = CAP.max_area(turbconv_params)
     n = n_mass_flux_subdomains(turbconv_model)
     (; edmfx_sgsflux_upwinding) = p.atmos.numerics
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶠu³, ᶜh_tot) = p.precomputed
     (; ᶜρaʲs, ᶜρʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs) = p.precomputed
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
@@ -237,11 +237,11 @@ function edmfx_sgs_mass_flux_tendency!(
             for j in 1:n
                 @. ᶠu³_diff = ᶠu³ʲs.:($$j) - ᶠu³
                 # @. ᶜa_scalar =
-                #     (ᶜq_totʲs.:($$j) - ᶜspecific.q_tot) *
+                #     (ᶜq_totʲs.:($$j) - specific(Y.c.ρq_tot, Y.c.ρ)) *
                 #     draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j))
                 # TODO: remove this filter when mass flux is treated implicitly
                 @. ᶜa_scalar =
-                    (ᶜq_totʲs.:($$j) - ᶜspecific.q_tot) * min(
+                    (ᶜq_totʲs.:($$j) - specific(Y.c.ρq_tot, Y.c.ρ)) * min(
                         min(draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)), a_max),
                         FT(0.02) / max(
                             Geometry.WVector(
@@ -277,17 +277,57 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; ᶜρa⁰, ᶜu⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰, ᶜtke⁰, ᶜmixing_length) = p.precomputed
+    (; ᶜρa⁰, ᶜu⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰, ᶜtke⁰, ᶠu³ʲs, ᶠu³⁰, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜdetrʲs) = p.precomputed
     if (
         p.atmos.moisture_model isa NonEquilMoistModel &&
         p.atmos.precip_model isa Microphysics1Moment
     )
         (; ᶜq_liq⁰, ᶜq_ice⁰, ᶜq_rai⁰, ᶜq_sno⁰) = p.precomputed
     end
-    (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
+    (; ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
 
     if p.atmos.edmfx_model.sgs_diffusive_flux isa Val{true}
+        (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+        (; params) = p
+
+        n = n_mass_flux_subdomains(turbconv_model)
+
+        ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
+        @. ᶜprandtl_nvec =
+            turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
+        sfc_tke = Fields.level(ᶜtke⁰, 1)
+        z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+        ᶜz = Fields.coordinate_field(Y.c).z
+        ᶜdz = Fields.Δz_field(axes(Y.c))
+
+        ᶜ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
+
+        ᶜmixing_length = @. lazy(master_mixing_length(
+            p.params,
+            ustar,
+            ᶜz,
+            z_sfc,
+            ᶜdz,
+            max(sfc_tke, eps(FT)),
+            ᶜlinear_buoygrad,
+            max(ᶜtke⁰, 0),
+            obukhov_length,
+            ᶜstrain_rate_norm,
+            ᶜprandtl_nvec,
+            ᶜtke_exch,
+            p.atmos.edmfx_model.scale_blending_method,
+        ))
+
+        ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length))
+        ᶜK_h = @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
         ᶠρaK_h = p.scratch.ᶠtemp_scalar
         @. ᶠρaK_h = ᶠinterp(ᶜρa⁰) * ᶠinterp(ᶜK_h)
         ᶠρaK_u = p.scratch.ᶠtemp_scalar
@@ -378,7 +418,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; ᶜu, ᶜh_tot, ᶜspecific, ᶜtke⁰, ᶜmixing_length) = p.precomputed
+    (; ᶜu, ᶜh_tot, ᶜtke⁰, ᶜmixing_length) = p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
 
@@ -424,7 +464,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
                 bottom = Operators.SetValue(C3(FT(0))),
             )
             @. ᶜρχₜ_diffusion =
-                ᶜdivᵥ_ρq_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜspecific.q_tot)))
+                ᶜdivᵥ_ρq_tot(-(ᶠρaK_h * ᶠgradᵥ(specific(Y.c.ρq_tot, Y.c.ρ))))
             @. Yₜ.c.ρq_tot -= ᶜρχₜ_diffusion
             @. Yₜ.c.ρ -= ᶜρχₜ_diffusion
         end
diff --git a/src/prognostic_equations/edmfx_tke.jl b/src/prognostic_equations/edmfx_tke.jl
index f37c37028c..eeccdef72e 100644
--- a/src/prognostic_equations/edmfx_tke.jl
+++ b/src/prognostic_equations/edmfx_tke.jl
@@ -50,8 +50,16 @@ function edmfx_tke_tendency!(
 )
     n = n_mass_flux_subdomains(turbconv_model)
     (; ᶜturb_entrʲs, ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p.precomputed
-    (; ᶠu³⁰, ᶠu³, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰) = p.precomputed
-    (; ᶜK_u, ᶜK_h) = p.precomputed
+    (;
+        ᶠu³⁰,
+        ᶠu³,
+        ᶜstrain_rate_norm,
+        ᶜlinear_buoygrad,
+        ᶜtke⁰,
+    ) = p.precomputed
+    turbconv_params = CAP.turbconv_params(p.params)
+    FT = eltype(p.params)
+
     ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶜρa⁰ : Y.c.ρ
     nh_pressure3_buoyʲs =
         turbconv_model isa PrognosticEDMFX ?
@@ -74,6 +82,46 @@ function edmfx_tke_tendency!(
     end
 
     if use_prognostic_tke(turbconv_model)
+        (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+        (; params) = p
+
+        ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
+        @. ᶜprandtl_nvec =
+            turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
+        sfc_tke = Fields.level(ᶜtke⁰, 1)
+        z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+        ᶜz = Fields.coordinate_field(Y.c).z
+        ᶜdz = Fields.Δz_field(axes(Y.c))
+
+        ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
+        @. ᶜtke_exch = 0
+        for j in 1:n
+            ᶠu³ʲ = ᶠu³ʲs.:($j)
+            @. ᶜtke_exch +=
+                Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / max(ᶜρa⁰, eps(FT)) *
+                (1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰)
+        end
+        
+        ᶜmixing_length = @. lazy(master_mixing_length(
+            p.params,
+            ustar,
+            ᶜz,
+            z_sfc,
+            ᶜdz,
+            max(sfc_tke, eps(FT)),
+            ᶜlinear_buoygrad,
+            max(ᶜtke⁰, 0),
+            obukhov_length,
+            ᶜstrain_rate_norm,
+            ᶜprandtl_nvec,
+            ᶜtke_exch,
+            p.atmos.edmfx_model.scale_blending_method,
+        ))
+
+        ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length))
+        ᶜK_h = @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
+
         # shear production
         @. Yₜ.c.sgs⁰.ρatke += 2 * ᶜρa⁰ * ᶜK_u * ᶜstrain_rate_norm
         # buoyancy production
diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl
index ce559020ef..7a28c65ae8 100644
--- a/src/prognostic_equations/hyperdiffusion.jl
+++ b/src/prognostic_equations/hyperdiffusion.jl
@@ -242,8 +242,10 @@ NVTX.@annotate function prep_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     (; ᶜspecific) = p.precomputed
     (; ᶜ∇²specific_tracers) = p.hyperdiff
 
-    for χ_name in propertynames(ᶜ∇²specific_tracers)
-        @. ᶜ∇²specific_tracers.:($$χ_name) = wdivₕ(gradₕ(ᶜspecific.:($$χ_name)))
+    # TODO: double check this change
+    for ρχ_name in filter(is_tracer_var, propertynames(Y.c))
+        is_energy_var(ρχ_name) || continue
+        @. ᶜ∇²specific_tracers.:($$χ_name) = wdivₕ(gradₕ(specific(Y.c.:$$χ_name, Y.c.ρ)))
     end
 
     if turbconv_model isa PrognosticEDMFX
diff --git a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
index 892e11c7f6..d1438d6600 100644
--- a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
+++ b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
@@ -603,8 +603,50 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
     end
 
     if use_derivative(diffusion_flag)
+        turbconv_params = CAP.turbconv_params(params)
+        FT = eltype(params)
+        
+        (; ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰) = p.precomputed
+        (; ᶜρa⁰, ᶜdetrʲs, ᶠu³⁰, ᶠu³ʲs) = p.precomputed
+        (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+
+        ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
+        @. ᶜprandtl_nvec =
+            turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
+        sfc_tke = Fields.level(ᶜtke⁰, 1)
+        z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+        ᶜz = Fields.coordinate_field(Y.c).z
+        ᶜdz = Fields.Δz_field(axes(Y.c))
+
+        n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+        ᶜ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
+
+        ᶜmixing_length = @. lazy(master_mixing_length(
+            p.params,
+            ustar,
+            ᶜz,
+            z_sfc,
+            ᶜdz,
+            max(sfc_tke, eps(FT)),
+            ᶜlinear_buoygrad,
+            max(ᶜtke⁰, 0),
+            obukhov_length,
+            ᶜstrain_rate_norm,
+            ᶜprandtl_nvec,
+            ᶜtke_exch,
+            p.atmos.edmfx_model.scale_blending_method,
+        ))
+        ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length))
+        ᶜK_h = @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
+
         α_vert_diff_tracer = CAP.α_vert_diff_tracer(params)
-        (; ᶜK_h, ᶜK_u) = p.precomputed
         @. ᶜdiffusion_h_matrix =
             ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅
             ᶠgradᵥ_matrix()
diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl
index 485b21599c..4385a73891 100644
--- a/src/utils/variable_manipulations.jl
+++ b/src/utils/variable_manipulations.jl
@@ -1,67 +1,101 @@
+"""
+    specific(ρχ, ρ)
+    specific(ρaχ, ρa, ρχ, ρ, turbconv_model)
+
+Calculates the specific quantity `χ` (per unit mass) from a density-weighted
+quantity. This function uses multiple dispatch to select the appropriate
+calculation method based on the number of arguments.
+
+**Grid-Scale Method (2 arguments)**
+
+    specific(ρχ, ρ)
+
+Performs a direct division of the density-weighted quantity `ρχ` by the density `ρ`.
+This method is used for grid-mean quantities where the density `ρ` is well-defined
+and non-zero.
+
+**SGS Regularized Method (5 arguments)**
+
+    specific(ρaχ, ρa, ρχ, ρ, turbconv_model)
+
+Calculates the specific quantity `χ` for a subgrid-scale (SGS) component by
+dividing the density-area-weighted quantity `ρaχ` by the density-area
+product `ρa`.
+
+This method includes regularization to handle cases where the SGS area fraction
+(and thus `ρa`) is zero or vanishingly small. It performs a linear interpolation
+between the SGS specific quantity (`ρaχ / ρa`) and the grid-mean specific
+quantity (`ρχ / ρ`). The interpolation weight is computed by `sgs_weight_function`
+to ensure a smooth and numerically stable transition, preventing division by zero.
+Using this regularized version instead of directly computing `ρaχ / ρa` breaks the
+assumption of domain decomposition (sum of SGS domains equals GS) when the approximated 
+area fraction `a` is small.
+
+Arguments:
+- `ρχ`: The grid-mean density-weighted quantity (e.g., `ρe_tot`, `ρq_tot`).
+- `ρ`: The grid-mean density.
+- `ρaχ`: The density-area-weighted SGS quantity (e.g., `sgs.ρa * sgs.h_tot`).
+- `ρa`: The density-area product of the SGS component.
+- `ρχ_fallback`: The grid-mean density-weighted quantity used for the fallback value.
+- `ρ_fallback`: The grid-mean density used for the fallback value.
+- `turbconv_model`: The turbulence convection model, containing parameters for regularization (e.g., `a_half`).
+"""
+function specific(ρχ, ρ)
+    return ρχ / ρ 
+end
+
+function specific(ρaχ, ρa, ρχ, ρ, turbconv_model)
+    # TODO: Replace turbconv_model struct by parameters, and include a_half in 
+    # parameters, not in config
+    weight = sgs_weight_function(ρa / ρ, turbconv_model.a_half)
+    # If ρa is exactly zero, the weight function will be zero, causing the first
+    # term to be NaN (0 * ... / 0). The ifelse handles this case explicitly.
+    return ρa == 0 ? ρχ / ρ : weight * ρaχ / ρa + (1 - weight) * ρχ / ρ
+end
+
 """
     sgs_weight_function(a, a_half)
 
-Computes the weight of the SGS variables in the linear interpolation used in
-`divide_by_ρa`. This is a continuously differentiable and monotonically
-increasing function of `a` that is equal to 0 when `a ≤ 0`, is equal to 1 when
-`a ≥ 1`, is equal to `1 / 2` when `a = a_half`, grows very rapidly near
-`a = a_half`, and grows very slowly at all other values of `a`.  If `a_half` is
-sufficiently small, this function is essentially equal to 1 for all `a` more
-than a few times larger than `a_half` (up to floating-point precision).
-
-We will now provide a description of how this function was constructed. We need
-the function to be equal to 0 when `a ≤ 0` and equal to 1 when `a ≥ 1`. Since
-the function must also be continuously differentiable, its derivative at these
-values of `a` has to be 0. To obtain a function with these properties, we use a
-piecewise definition:
-    - For all `a < 0`, the function is equal to 0.
-    - For all `a > 1`, the function is equal to 1.
-    - For all `0 ≤ a ≤ 1`, the function is a sigmoid that connects the point
-      `(0, 0)` to the point `(1, 1)`, with a derivative of 0 at these points.
-Most well-known sigmoid functions connect the "points" `(-Inf, 0)` and
-`(Inf, 1)`, not `(0, 0)` and `(1, 1)`. To obtain the desired sigmoid curve, we
-begin with two simple sigmoid functions that go from `(-Inf, 0)` to `(Inf, 1)`
-at different rates. In this case, we use two `tanh` functions, scaled and
-translated so that they lie between 0 and 1:
-    - `fast_sigmoid(a) = (1 + tanh(a)) / 2` and
-    - `slow_sigmoid(a) = (1 + tanh(a / 2)) / 2`.
-Note that the second sigmoid is commonly called the "logistic" function. We then
-take the inverse of the sigmoid that grows more slowly, and we make that the
-input of the sigmoid that grows more quickly:
-    - `sigmoid(a) = fast_sigmoid(slow_sigmoid⁻¹(a)) =
-       (1 + tanh(2 * atanh(2 * a - 1))) / 2`.
-The resulting function goes from `(0, 0)` to `(1, 1)`, and, since the outer
-sigmoid grows more quickly, it has the same asymptotic behavior as the outer
-sigmoid, which means that its derivative at the boundary points is 0. If we had
-instead put the sigmoid that grows more slowly on the outside, the asymptotic
-behavior would come from the inverted inner sigmoid, which means that the
-derivative at the boundary points would be `Inf`.
-
-The sigmoid function we have constructed reaches `1 / 2` when `a = 1 / 2`. More
-generally, we need the weight function to reach `1 / 2` when `a` is some small
-value `a_half`. To achieve this, we replace the input to the sigmoid function
-with a smooth, monotonically increasing function that goes through `(0, 0)`,
-`(a_half, 1 / 2)`, and `(1, 1)`. The simplest option is the power function
-    - `power(a) = a^(-1 / log2(a_half))`.
-However, this function does not work well because, when `a_half < 1 / 2`, its
-derivative at `a = 0` is `Inf`, and its derivative at `a = 1` is some positive
-number. Making this power function the input to the sigmoid function causes the
-derivative of the sigmoid to become `Inf` at `a = 0` when `a_half < 1 / 4`, and
-it causes the sigmoid to grow too slowly from `1 / 2` to 1, only reaching 1 when
-`a` is significantly larger than `a_half`. In order to fix this, we transform
-the power function by replacing `a` with `1 - a`, `a_half` with `1 - a_half`,
-and `power(a)` with `1 - power(a)`, which gives us
-    - `power(a) = 1 - (1 - a)^(-1 / log2(1 - a_half))`.
-This transformed function works better because, when `a_half < 1 / 2`, its
-derivative at `a = 0` is some positive number, and its derivative at `a = 1` is
-0. When we make this the input to the sigmoid function, the result has a
-continuous derivative and is essentially equal to 1 for all `a` more than a few
-times larger than `a_half`. So, for all `0 ≤ a ≤ 1`, we define the weight
-function as
-    - `weight(a) = sigmoid(power(a)) =
-       (1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2`.
-"""
-sgs_weight_function(a, a_half) =
+Computes a smooth, monotonic weight function `w(a)` that ranges from 0 to 1.
+
+This function is used as the interpolation weight in the regularized `specific`
+function. It ensures a numerically stable and smooth transition between a subgrid-scale 
+(SGS) quantity and its grid-mean counterpart, especially when the SGS area fraction `a` 
+is small.
+
+**Key Properties:**
+- `w(a) = 0` for `a ≤ 0`.
+- `w(a) = 1` for `a ≥ 1`.
+- `w(a_half) = 0.5`.
+- The function is continuously differentiable, with derivatives equal to zero at
+  `a = 0` and `a = 1`, which ensures smooth blending.
+- The functions grows very rapidly near `a = a_half`, and grows very slowly at all other 
+  values of `a`.
+- For small `a_half`, the weight rapidly approaches 1 for values of `a` that are
+  a few times larger than `a_half`.
+
+**Construction Method:**
+The function is piecewise. For `a` between 0 and 1, it is a custom sigmoid curve
+constructed in two main steps to satisfy the key properties:
+1.  **Bounded Sigmoid Creation**: A base sigmoid is created that maps the interval
+    `(0, 1)` to `(0, 1)` with zero derivatives at the endpoints. This is achieved
+    by composing a standard `tanh` function with the inverse of a slower-growing
+    `tanh` function.
+2.  **Midpoint Control**: To ensure the function passes through the control point
+    `(a_half, 0.5)`, the input `a` is first transformed by a specially designed
+    power function (`1 - (1 - a)^k`) before being passed to the bounded sigmoid.
+    This transformation maps `a_half` to `0.5` while preserving differentiability 
+    at the boundaries.
+
+Arguments:
+- `a`: The input SGS area fraction (often approximated as `ρa / ρ`).
+- `a_half`: The value of `a` at which the weight function should be 0.5, controlling
+          the transition point of the sigmoid curve.
+
+Returns:
+- The computed weight, a value between 0 and 1.
+"""
+function sgs_weight_function(a, a_half)
     if a < 0
         zero(a)
     elseif a > 1
@@ -69,26 +103,6 @@ sgs_weight_function(a, a_half) =
     else
         (1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2
     end
-
-"""
-    divide_by_ρa(ρaχ, ρa, ρχ, ρ, turbconv_model)
-
-Computes `ρaχ / ρa`, regularizing the result to avoid issues when `a` is small.
-This is done by performing a linear interpolation from `ρaχ / ρa` to `ρχ / ρ`,
-using `sgs_weight_function(ρa / ρ, turbconv_model.a_half)` as the weight of
-`ρaχ / ρa` in the interpolation. Note that `ρa / ρ` is the "anelastic
-approximation" of `a`; we cannot directly use `a` to compute the weight because
-this function needs to be called before `a` has been computed. Also, note that
-using this function instead of directly computing `ρaχ / ρa` breaks the
-assumption of domain decomposition when the approximated `a` is small.
-"""
-function divide_by_ρa(ρaχ, ρa, ρχ, ρ, turbconv_model)
-    weight = sgs_weight_function(ρa / ρ, turbconv_model.a_half)
-    # If ρa = 0, we know that ρa / ρ = 0, which means that weight = 0. However,
-    # 0 * ρaχ / 0 = NaN, regardless of what ρaχ is, so the linear interpolation
-    # will always return NaN when ρa = 0. To avoid this problem, we need to add
-    # a special case for ρa = 0.
-    return ρa == 0 ? ρχ / ρ : weight * ρaχ / ρa + (1 - weight) * ρχ / ρ
 end
 
 # Helper functions for manipulating symbols in the generated functions:
@@ -121,7 +135,7 @@ end
 Converts every variable of the form `ρaχ` in the sub-grid-scale state `sgs` into
 the specific variable `χ` by dividing it by `ρa`. All other variables in `sgs`
 are omitted from the result. The division is computed as
-`divide_by_ρa(ρaχ, ρa, ρχ, ρ, turbconv_model)`, which is preferable to simply
+`specific(ρaχ, ρa, ρχ, ρ, turbconv_model)`, which is preferable to simply
 calling `ρaχ / ρa` because it avoids numerical issues that arise when `a` is
 small. The values of `ρ` and `ρχ` are taken from `gs`, but, when `ρχ` is not
 available in `gs` (e.g., when `χ` is a second moment variable like `tke`), its
@@ -136,7 +150,7 @@ value is assumed to be equal to the value of `ρaχ` in `sgs`.
         map(name -> remove_prefix(name, :ρa), relevant_sgs_names)
     relevant_gs_names = map(name -> Symbol(:ρ, name), specific_sgs_names)
     specific_sgs_values = map(
-        (sgs_name, gs_name) -> :(divide_by_ρa(
+        (sgs_name, gs_name) -> :(specific(
             sgs.$sgs_name,
             sgs.ρa,
             $(gs_name in gs_names ? :(gs.$gs_name) : :(sgs.$sgs_name)),
@@ -260,7 +274,7 @@ are computed from the tuples of subdomain densities and velocities `ρaʲs` and
 `u₃ʲs`. The division is computed using `divide_by_ρa` to avoid issues when `a⁺`
 is small.
 """
-u₃⁺(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) = divide_by_ρa(
+u₃⁺(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) = specific(
     unrolled_dotproduct(ρaʲs, u₃ʲs),
     reduce(+, ρaʲs),
     ρ * u₃,
@@ -278,7 +292,7 @@ are computed from the domain decomposition of the grid-scale quantities `ρw` an
 environment quantities. The division is computed using `divide_by_ρa` to avoid
 issues when `a⁰` is small.
 """
-u₃⁰(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) = divide_by_ρa(
+u₃⁰(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) = specific(
     ρ * u₃ - unrolled_dotproduct(ρaʲs, u₃ʲs),
     ρ - reduce(+, ρaʲs),
     ρ * u₃,