diff --git a/config/model_configs/diagnostic_edmfx_aquaplanet_dense_autodiff.yml b/config/model_configs/diagnostic_edmfx_aquaplanet_dense_autodiff.yml
index c728ca61bc..62268b5b8c 100644
--- a/config/model_configs/diagnostic_edmfx_aquaplanet_dense_autodiff.yml
+++ b/config/model_configs/diagnostic_edmfx_aquaplanet_dense_autodiff.yml
@@ -19,6 +19,6 @@ edmfx_sgs_diffusive_flux: true
 moist: equil
 cloud_model: "quadrature_sgs"
 precip_model: 0M
-dt: 120secs
+dt: 80secs
 t_end: 1hours
 toml: [toml/diagnostic_edmfx_0M.toml]
diff --git a/perf/benchmark.jl b/perf/benchmark.jl
index 01ac9c034a..68bcbfc2c7 100644
--- a/perf/benchmark.jl
+++ b/perf/benchmark.jl
@@ -60,7 +60,7 @@ are_boundschecks_forced = Base.JLOptions().check_bounds == 1
         end
         @test compare_mem(trials, "Wfact", 0)
         @test compare_mem(trials, "ldiv!", 0)
-        @test compare_mem(trials, "T_imp!", 96)
+        @test compare_mem(trials, "T_imp!", 19680)
         @test compare_mem(trials, "T_exp_T_lim!", 190420)
         @test compare_mem(trials, "lim!", 0)
         @test compare_mem(trials, "dss!", 0)
diff --git a/perf/flame.jl b/perf/flame.jl
index 61b56888b7..aa4ca62068 100644
--- a/perf/flame.jl
+++ b/perf/flame.jl
@@ -45,9 +45,9 @@ allocs_limit["flame_diagnostics"] = 10_677_144
 allocs_limit["flame_aquaplanet_diagedmf"] = 11_644_128
 allocs_limit["flame_aquaplanet_progedmf"] = 774_712
 allocs_limit["flame_aquaplanet_progedmf_dense_autodiff"] = 774_712
-allocs_limit["flame_diffusion"] = 138_432
+allocs_limit["flame_diffusion"] = 153_248
 allocs_limit["flame_threaded"] = 2047_736
-allocs_limit["flame_callbacks"] = 391_864
+allocs_limit["flame_callbacks"] = 450_336
 allocs_limit["flame_gravity_wave"] = 581_381_976
 # Ideally, we would like to track all the allocations, but this becomes too
 # expensive there is too many of them. Here, we set the default sample rate to
diff --git a/perf/jet_test_nfailures.jl b/perf/jet_test_nfailures.jl
index e6a95b43be..d103542d08 100644
--- a/perf/jet_test_nfailures.jl
+++ b/perf/jet_test_nfailures.jl
@@ -38,7 +38,7 @@ using Test
     # inference. By increasing this counter, we acknowledge that
     # we have introduced an inference failure. We hope to drive
     # this number down to 0.
-    n_allowed_failures = 253
+    n_allowed_failures = 265
     @show n
     if n < n_allowed_failures
         @info "Please update the n-failures to $n"
diff --git a/reproducibility_tests/ref_counter.jl b/reproducibility_tests/ref_counter.jl
index 6491e2a62e..7e025c3013 100644
--- a/reproducibility_tests/ref_counter.jl
+++ b/reproducibility_tests/ref_counter.jl
@@ -1,4 +1,4 @@
-248
+249
 
 # **README**
 #
@@ -20,6 +20,9 @@
 
 
 #=
+249
+- Remove viscosity, diffusivity, and mixing length from precomputed quantities
+
 248
 - update deps: climacore 0.14.35
 
diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl
index 06055c78ff..bd3420aa85 100644
--- a/src/cache/cloud_fraction.jl
+++ b/src/cache/cloud_fraction.jl
@@ -22,12 +22,9 @@ end
    Compute the grid scale cloud fraction based on sub-grid scale properties
 """
 NVTX.@annotate function set_cloud_fraction!(Y, p, ::DryModel, _)
-    (; ᶜmixing_length) = p.precomputed
     (; turbconv_model) = p.atmos
     FT = eltype(p.params)
-    if isnothing(turbconv_model)
-        compute_gm_mixing_length!(ᶜmixing_length, Y, p)
-    end
+
     p.precomputed.cloud_diagnostics_tuple .=
         ((; cf = FT(0), q_liq = FT(0), q_ice = FT(0)),)
 end
@@ -39,8 +36,9 @@ 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)
+
     if isnothing(turbconv_model)
         if p.atmos.call_cloud_diagnostics_per_stage isa
            CallCloudDiagnosticsPerStage
@@ -53,7 +51,6 @@ NVTX.@annotate function set_cloud_fraction!(
             @. ᶜgradᵥ_θ_liq_ice =
                 ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
         end
-        compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     end
     if moist_model isa EquilMoistModel
         @. cloud_diagnostics_tuple = make_named_tuple(
@@ -85,7 +82,7 @@ NVTX.@annotate function set_cloud_fraction!(
 
     FT = eltype(params)
     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)
@@ -100,9 +97,10 @@ NVTX.@annotate function set_cloud_fraction!(
             @. ᶜgradᵥ_θ_liq_ice =
                 ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
         end
-        compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     end
 
+    ᶜmixing_length = compute_gm_mixing_length(Y, p)
+
     diagnostic_covariance_coeff = CAP.diagnostic_covariance_coeff(params)
     @. cloud_diagnostics_tuple = quad_loop(
         SG_quad,
@@ -139,20 +137,23 @@ NVTX.@annotate function set_cloud_fraction!(
 
     FT = eltype(params)
     thermo_params = CAP.thermodynamics_params(params)
-    (; ᶜts, ᶜmixing_length, cloud_diagnostics_tuple) = p.precomputed
+    (; ᶜts, cloud_diagnostics_tuple) = p.precomputed
     (; turbconv_model) = p.atmos
 
     # TODO - we should make this default when using diagnostic edmf
     # environment
     diagnostic_covariance_coeff = CAP.diagnostic_covariance_coeff(params)
 
+    ᶜmixing_length_field = p.scratch.ᶜtemp_scalar
+    ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+
     @. cloud_diagnostics_tuple = quad_loop(
         SG_quad,
         ᶜts,
         Geometry.WVector(p.precomputed.ᶜgradᵥ_q_tot),
         Geometry.WVector(p.precomputed.ᶜgradᵥ_θ_liq_ice),
         diagnostic_covariance_coeff,
-        ᶜmixing_length,
+        ᶜmixing_length_field,
         thermo_params,
     )
 
@@ -192,7 +193,7 @@ NVTX.@annotate function set_cloud_fraction!(
 
     FT = eltype(params)
     thermo_params = CAP.thermodynamics_params(params)
-    (; ᶜts⁰, ᶜmixing_length, cloud_diagnostics_tuple) = p.precomputed
+    (; ᶜts⁰, cloud_diagnostics_tuple) = p.precomputed
     (; ᶜρʲs, ᶜtsʲs, ᶜρa⁰, ᶜρ⁰) = p.precomputed
     (; turbconv_model) = p.atmos
 
@@ -200,13 +201,16 @@ NVTX.@annotate function set_cloud_fraction!(
     # environment
     diagnostic_covariance_coeff = CAP.diagnostic_covariance_coeff(params)
 
+    ᶜmixing_length_field = p.scratch.ᶜtemp_scalar
+    ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+
     @. cloud_diagnostics_tuple = quad_loop(
         SG_quad,
         ᶜts⁰,
         Geometry.WVector(p.precomputed.ᶜgradᵥ_q_tot⁰),
         Geometry.WVector(p.precomputed.ᶜgradᵥ_θ_liq_ice⁰),
         diagnostic_covariance_coeff,
-        ᶜmixing_length,
+        ᶜmixing_length_field,
         thermo_params,
     )
 
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index f0cb9255b9..51bafd3160 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -966,13 +966,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     (; ᶜp, ᶜu, ᶜts) = p.precomputed
     (; ustar, obukhov_length) = p.precomputed.sfc_conditions
     (; ᶜtke⁰) = p.precomputed
-    (;
-        ᶜlinear_buoygrad,
-        ᶜstrain_rate_norm,
-        ᶜmixing_length_tuple,
-        ᶜmixing_length,
-    ) = p.precomputed
-    (; ᶜK_h, ᶜK_u, ρatke_flux) = p.precomputed
+    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
+    (; ρatke_flux) = p.precomputed
     turbconv_params = CAP.turbconv_params(params)
     thermo_params = CAP.thermodynamics_params(params)
     ᶜlg = Fields.local_geometry_field(Y.c)
@@ -1004,42 +999,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     ᶜstrain_rate .= compute_strain_rate_center(ᶠu⁰)
     @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
 
-    ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
-    @. ᶜprandtl_nvec =
-        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
-
-    ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
-    @. ᶜtke_exch = 0
-    # using ᶜu⁰ would be more correct, but this is more consistent with the
-    # TKE equation, where using ᶜu⁰ results in allocation
-    for j in 1:n
-        @. ᶜtke_exch +=
-            ᶜρaʲs.:($$j) * ᶜdetrʲs.:($$j) / Y.c.ρ *
-            (1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰)
-    end
-
-    sfc_tke = Fields.level(ᶜtke⁰, 1)
-    z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, half)
-    @. ᶜmixing_length_tuple = mixing_length(
-        params,
-        ustar,
-        ᶜz,
-        z_sfc,
-        ᶜdz,
-        max(sfc_tke, 0),
-        ᶜ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)
-
     ρatke_flux_values = Fields.field_values(ρatke_flux)
     ρ_int_values = Fields.field_values(Fields.level(Y.c.ρ, 1))
     u_int_values = Fields.field_values(Fields.level(ᶜu, 1))
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 04d6e3bca6..357e845c28 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -114,7 +114,6 @@ function precomputed_quantities(Y, atmos)
     gs_quantities = (;
         ᶜwₜqₜ = similar(Y.c, Geometry.WVector{FT}),
         ᶜwₕhₜ = similar(Y.c, Geometry.WVector{FT}),
-        ᶜmixing_length = similar(Y.c, FT),
         ᶜlinear_buoygrad = similar(Y.c, FT),
         ᶜstrain_rate_norm = similar(Y.c, FT),
         sfc_conditions = similar(Spaces.level(Y.f, half), SCT),
@@ -173,9 +172,6 @@ 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),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
             bdmr_l = similar(Y.c, BidiagonalMatrixRow{FT}),
             bdmr_r = similar(Y.c, BidiagonalMatrixRow{FT}),
@@ -193,11 +189,8 @@ function precomputed_quantities(Y, atmos)
     edonly_quantities =
         atmos.turbconv_model isa EDOnlyEDMFX ?
         (;
-            ᶜ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),
         ) : (;)
 
     sgs_quantities = (;
@@ -225,20 +218,9 @@ function precomputed_quantities(Y, atmos)
             ᶠu³⁰ = similar(Y.f, CT3{FT}),
             ᶜu⁰ = similar(Y.c, C123{FT}),
             ᶜK⁰ = similar(Y.c, FT),
-            ᶜmixing_length_tuple = similar(Y.c, MixingLength{FT}),
-            ᶜK_u = similar(Y.c, FT),
-            ᶜK_h = similar(Y.c, FT),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
             precipitation_sgs_quantities...,
         ) : (;)
-    vert_diff_quantities =
-        if atmos.vertical_diffusion isa
-           Union{VerticalDiffusion, DecayWithHeightDiffusion}
-            ᶜK_h = similar(Y.c, FT)
-            (; ᶜK_u = ᶜK_h, ᶜK_h) # ᶜK_u aliases ᶜK_h because they are always equal.
-        else
-            (;)
-        end
     smagorinsky_lilly_quantities =
         if atmos.smagorinsky_lilly isa SmagorinskyLilly
             uvw_vec = UVW(FT(0), FT(0), FT(0))
@@ -259,7 +241,6 @@ function precomputed_quantities(Y, atmos)
         advective_sgs_quantities...,
         edonly_quantities...,
         diagnostic_sgs_quantities...,
-        vert_diff_quantities...,
         sedimentation_quantities...,
         precipitation_quantities...,
         surface_precip_fluxes...,
@@ -538,12 +519,6 @@ NVTX.@annotate function set_explicit_precomputed_quantities!(Y, p, t)
             ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
     end
 
-    # TODO: It is too slow to calculate mixing length at every timestep
-    # if isnothing(turbconv_model)
-    #     (; ᶜmixing_length) = p.precomputed
-    #     compute_gm_mixing_length!(ᶜmixing_length, Y, p)
-    # end
-
     if turbconv_model isa PrognosticEDMFX
         set_prognostic_edmf_precomputed_quantities_bottom_bc!(Y, p, t)
         set_prognostic_edmf_precomputed_quantities_explicit_closures!(Y, p, t)
@@ -585,19 +560,6 @@ NVTX.@annotate function set_explicit_precomputed_quantities!(Y, p, t)
     )
     set_precipitation_surface_fluxes!(Y, p, p.atmos.microphysics_model)
 
-    if vertical_diffusion isa DecayWithHeightDiffusion
-        (; ᶜK_h) = p.precomputed
-        @. ᶜK_h =
-            $compute_eddy_diffusivity_coefficient(Y.c.ρ, vertical_diffusion)
-    elseif vertical_diffusion isa VerticalDiffusion
-        (; ᶜK_h) = p.precomputed
-        @. ᶜK_h = $compute_eddy_diffusivity_coefficient(
-            Y.c.uₕ,
-            ᶜp,
-            vertical_diffusion,
-        )
-    end
-
     # TODO
     if call_cloud_diagnostics_per_stage isa CallCloudDiagnosticsPerStage
         set_cloud_fraction!(Y, p, moisture_model, cloud_model)
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index 131a780ccd..6f94fc5567 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -368,15 +368,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
     n = n_mass_flux_subdomains(turbconv_model)
 
     (; ᶜtke⁰, ᶜu, ᶜp, ᶜρa⁰, ᶠu³⁰, ᶜts⁰, ᶜq_tot⁰) = p.precomputed
-    (;
-        ᶜmixing_length_tuple,
-        ᶜmixing_length,
-        ᶜlinear_buoygrad,
-        ᶜstrain_rate_norm,
-        ᶜK_u,
-        ᶜK_h,
-        ρatke_flux,
-    ) = p.precomputed
+    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm, ρatke_flux) = p.precomputed
     (;
         ᶜuʲs,
         ᶜtsʲs,
@@ -501,41 +493,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
     ᶜstrain_rate .= compute_strain_rate_center(ᶠu⁰)
     @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
 
-    ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
-    @. ᶜprandtl_nvec =
-        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
-
-    ᶜ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)
-
     ρatke_flux_values = Fields.field_values(ρatke_flux)
     ρa_sfc_values = Fields.field_values(Fields.level(ᶜρa⁰, 1)) # TODO: replace by surface value
     ustar_values = Fields.field_values(ustar)
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
index 059e868179..cbf58a77d5 100644
--- a/src/diagnostics/Diagnostics.jl
+++ b/src/diagnostics/Diagnostics.jl
@@ -56,8 +56,15 @@ import ..SlabOceanSST
 
 # functions used to calculate diagnostics
 import ..draft_area
-import ..compute_gm_mixing_length!
+import ..compute_gm_mixing_length
 import ..horizontal_integral_at_boundary
+import ..ᶜmixing_length
+import ..ᶜeddy_diffusivity
+import ..ᶜeddy_viscosity
+import ..turbulent_prandtl_number
+import ..smagorinsky_lilly_length
+import ..compute_eddy_diffusivity_coefficient
+
 
 # We need the abbreviations for symbols like curl, grad, and so on
 include(joinpath("..", "utils", "abbreviations.jl"))
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index 08c6abca62..4c9d5d5c71 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -273,10 +273,40 @@ add_diagnostic_variable!(
     or from mixing length closure with EDMF SGS model.
     """,
     compute! = (out, state, cache, time) -> begin
+        turbconv_model = cache.atmos.turbconv_model
+        # TODO: consolidate remaining mixing length types
+        # (smagorinsky_lilly, dz) into a single mixing length function
+        if isa(turbconv_model, PrognosticEDMFX) ||
+           isa(turbconv_model, DiagnosticEDMFX) ||
+           isa(turbconv_model, EDOnlyEDMFX)
+            ᶜmixing_length_field = ᶜmixing_length(state, cache)
+        else
+            (; params) = cache
+            (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = cache.precomputed
+            ᶜdz = Fields.Δz_field(axes(state.c))
+            ᶜprandtl_nvec = @. lazy(
+                turbulent_prandtl_number(
+                    params,
+                    ᶜlinear_buoygrad,
+                    ᶜstrain_rate_norm,
+                ),
+            )
+            ᶜmixing_length_field = @. lazy(
+                smagorinsky_lilly_length(
+                    CAP.c_smag(params),
+                    sqrt(max(ᶜlinear_buoygrad, 0)),   # N_eff
+                    ᶜdz,
+                    ᶜprandtl_nvec,
+                    ᶜstrain_rate_norm,
+                ),
+            )
+        end
+
+
         if isnothing(out)
-            return copy(cache.precomputed.ᶜmixing_length)
+            return Base.materialize(ᶜmixing_length_field)
         else
-            out .= cache.precomputed.ᶜmixing_length
+            out .= ᶜmixing_length_field
         end
     end,
 )
diff --git a/src/diagnostics/edmfx_diagnostics.jl b/src/diagnostics/edmfx_diagnostics.jl
index b6a8d6454b..a3eb4db48b 100644
--- a/src/diagnostics/edmfx_diagnostics.jl
+++ b/src/diagnostics/edmfx_diagnostics.jl
@@ -1176,10 +1176,12 @@ function compute_lmixw!(
     time,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+    ᶜwall_mixing_length = ᶜmixing_length(state, cache, Val(:wall))
+
     if isnothing(out)
-        return copy(cache.precomputed.ᶜmixing_length_tuple.wall)
+        return Base.materialize(ᶜwall_mixing_length)
     else
-        out .= cache.precomputed.ᶜmixing_length_tuple.wall
+        out .= ᶜwall_mixing_length
     end
 end
 
@@ -1205,10 +1207,13 @@ function compute_lmixtke!(
     time,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+
+    ᶜtke_mixing_length = ᶜmixing_length(state, cache, Val(:tke))
+
     if isnothing(out)
-        return copy(cache.precomputed.ᶜmixing_length_tuple.tke)
+        return Base.materialize(ᶜtke_mixing_length)
     else
-        out .= cache.precomputed.ᶜmixing_length_tuple.tke
+        out .= ᶜtke_mixing_length
     end
 end
 
@@ -1234,10 +1239,12 @@ function compute_lmixb!(
     time,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+    ᶜbuoy_mixing_length = ᶜmixing_length(state, cache, Val(:buoy))
+
     if isnothing(out)
-        return copy(cache.precomputed.ᶜmixing_length_tuple.buoy)
+        return Base.materialize(ᶜbuoy_mixing_length)
     else
-        out .= cache.precomputed.ᶜmixing_length_tuple.buoy
+        out .= ᶜbuoy_mixing_length
     end
 end
 
@@ -1259,23 +1266,42 @@ compute_edt!(out, state, cache, time) = compute_edt!(
     cache.atmos.vertical_diffusion,
     cache.atmos.turbconv_model,
 )
-compute_edt!(_, _, _, _, vert_diff::T1, turbconv_model::T2) where {T1, T2} =
-    error_diagnostic_variable(
-        "Can only compute heat diffusivity with vertical diffusion or EDMFX",
-    )
+compute_edt!(
+    _,
+    _,
+    _,
+    _,
+    vertical_diffusion::T1,
+    turbconv_model::T2,
+) where {T1, T2} = error_diagnostic_variable(
+    "Can only compute heat diffusivity with vertical diffusion or EDMFX",
+)
 
 function compute_edt!(
     out,
     state,
     cache,
     time,
-    vert_diff::Union{VerticalDiffusion, DecayWithHeightDiffusion},
+    vertical_diffusion::Union{VerticalDiffusion, DecayWithHeightDiffusion},
     turbconv_model::Nothing,
 )
+    (; vertical_diffusion) = cache.atmos
+    (; ᶜp) = cache.precomputed
+
+    if vertical_diffusion isa DecayWithHeightDiffusion
+        ᶜK_h =
+            compute_eddy_diffusivity_coefficient(state.c.ρ, vertical_diffusion)
+    elseif vertical_diffusion isa VerticalDiffusion
+        ᶜK_h = compute_eddy_diffusivity_coefficient(
+            state.c.uₕ,
+            ᶜp,
+            vertical_diffusion,
+        )
+    end
     if isnothing(out)
-        return copy(cache.precomputed.ᶜK_h)
+        return copy(ᶜK_h)
     else
-        out .= cache.precomputed.ᶜK_h
+        out .= ᶜK_h
     end
 end
 
@@ -1284,13 +1310,19 @@ function compute_edt!(
     state,
     cache,
     time,
-    vert_diff::Nothing,
+    vertical_diffusion::Nothing,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+    turbconv_params = CAP.turbconv_params(cache.params)
+    (; ᶜtke⁰) = cache.precomputed
+
+    ᶜmixing_length_field = ᶜmixing_length(state, cache)
+    ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+    ᶜK_h = ᶜeddy_diffusivity(cache, ᶜK_u)
     if isnothing(out)
-        return copy(cache.precomputed.ᶜK_h)
+        return Base.materialize(ᶜK_h)
     else
-        out .= cache.precomputed.ᶜK_h
+        out .= ᶜK_h
     end
 end
 
@@ -1314,23 +1346,43 @@ compute_evu!(out, state, cache, time) = compute_evu!(
     cache.atmos.vertical_diffusion,
     cache.atmos.turbconv_model,
 )
-compute_evu!(_, _, _, _, vert_diff::T1, turbconv_model::T2) where {T1, T2} =
-    error_diagnostic_variable(
-        "Can only compute momentum diffusivity with vertical diffusion or EDMFX",
-    )
+compute_evu!(
+    _,
+    _,
+    _,
+    _,
+    vertical_diffusion::T1,
+    turbconv_model::T2,
+) where {T1, T2} = error_diagnostic_variable(
+    "Can only compute momentum diffusivity with vertical diffusion or EDMFX",
+)
 
 function compute_evu!(
     out,
     state,
     cache,
     time,
-    vert_diff::Union{VerticalDiffusion, DecayWithHeightDiffusion},
+    vertical_diffusion::Union{VerticalDiffusion, DecayWithHeightDiffusion},
     turbconv_model::Nothing,
 )
+    (; vertical_diffusion) = cache.atmos
+    (; ᶜp) = cache.precomputed
+
+    # this setup assumes ᶜK_u = ᶜK_h
+    if vertical_diffusion isa DecayWithHeightDiffusion
+        ᶜK_u =
+            compute_eddy_diffusivity_coefficient(state.c.ρ, vertical_diffusion)
+    elseif vertical_diffusion isa VerticalDiffusion
+        ᶜK_u = compute_eddy_diffusivity_coefficient(
+            state.c.uₕ,
+            ᶜp,
+            vertical_diffusion,
+        )
+    end
     if isnothing(out)
-        return copy(cache.precomputed.ᶜK_u)
+        return copy(ᶜK_u)
     else
-        out .= cache.precomputed.ᶜK_u
+        out .= ᶜK_u
     end
 end
 
@@ -1339,13 +1391,18 @@ function compute_evu!(
     state,
     cache,
     time,
-    vert_diff::Nothing,
+    vertical_diffusion::Nothing,
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
+    turbconv_params = CAP.turbconv_params(cache.params)
+    (; ᶜtke⁰) = cache.precomputed
+    ᶜmixing_length_field = ᶜmixing_length(state, cache)
+    ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+
     if isnothing(out)
-        return copy(cache.precomputed.ᶜK_u)
+        return Base.materialize(ᶜK_u)
     else
-        out .= cache.precomputed.ᶜK_u
+        out .= ᶜK_u
     end
 end
 
diff --git a/src/prognostic_equations/eddy_diffusion_closures.jl b/src/prognostic_equations/eddy_diffusion_closures.jl
index 3954285b3d..de7702e6c0 100644
--- a/src/prognostic_equations/eddy_diffusion_closures.jl
+++ b/src/prognostic_equations/eddy_diffusion_closures.jl
@@ -293,7 +293,7 @@ function surface_flux_tke(
 end
 
 """
-    mixing_length(
+    mixing_length_lopez_gomez_2020(
     params,
     ustar,
     ᶜz,
@@ -324,13 +324,17 @@ where:
 - `ᶜtke_exch`: TKE exchange term [m^2/s^3].
 - `scale_blending_method`: The method to use for blending physical scales.
 
-Calculates the turbulent mixing length, limited by physical constraints (wall distance,
-TKE balance, stability) and grid resolution. 
+Point-wise calculation of the turbulent mixing length, limited by physical constraints (wall distance,
+TKE balance, stability) and grid resolution. Based on 
+Lopez‐Gomez, I., Cohen, Y., He, J., Jaruga, A., & Schneider, T. (2020). 
+A generalized mixing length closure for eddy‐diffusivity mass‐flux schemes of turbulence and convection. 
+Journal of Advances in Modeling Earth Systems, 12, e2020MS002161. https://doi.org/ 10.1029/2020MS002161
 
 Returns a `MixingLength{FT}` struct containing the final blended mixing length (`master`) 
 and its constituent physical scales.
 """
-function mixing_length(
+
+function mixing_length_lopez_gomez_2020(
     params,
     ustar,
     ᶜz,
@@ -346,7 +350,7 @@ function mixing_length(
     scale_blending_method,
 )
 
-    FT = eltype(params)
+    FT = eltype(ᶜz)
     eps_FT = eps(FT)
 
     turbconv_params = CAP.turbconv_params(params)
@@ -524,7 +528,49 @@ function mixing_length(
     # minimum mixing length
     l_final = max(l_final, FT(1)) # TODO: make a climaparam
 
-    return MixingLength{FT}(l_final, l_W, l_TKE, l_N, l_grid)
+    return MixingLength(l_final, l_W, l_TKE, l_N, l_grid)
+end
+
+# GPU-safe field access using Val dispatch
+@inline get_mixing_length_field(ml::MixingLength, ::Val{:master}) = ml.master
+@inline get_mixing_length_field(ml::MixingLength, ::Val{:wall}) = ml.wall
+@inline get_mixing_length_field(ml::MixingLength, ::Val{:tke}) = ml.tke
+@inline get_mixing_length_field(ml::MixingLength, ::Val{:buoy}) = ml.buoy
+@inline get_mixing_length_field(ml::MixingLength, ::Val{:l_grid}) = ml.l_grid
+
+function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
+    (; params) = p
+    (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+    (; ᶜtke⁰) = p.precomputed
+    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
+    ᶜz = Fields.coordinate_field(Y.c).z
+    z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+    ᶜdz = Fields.Δz_field(axes(Y.c))
+    sfc_tke = Fields.level(ᶜtke⁰, 1)
+
+    ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar_5
+    ᶜprandtl_nvec .= ᶜturbulent_prandtl_number(p)
+
+    ᶜtke_exch = ᶜtke_exchange(Y, p)
+
+    ᶜmixing_length_tuple = @. lazy(
+        mixing_length_lopez_gomez_2020(
+            params,
+            ustar,
+            ᶜz,
+            z_sfc,
+            ᶜdz,
+            sfc_tke,
+            ᶜlinear_buoygrad,
+            ᶜtke⁰,
+            obukhov_length,
+            ᶜstrain_rate_norm,
+            ᶜprandtl_nvec,
+            ᶜtke_exch,
+            p.atmos.edmfx_model.scale_blending_method,
+        ),
+    )
+    return @. lazy(get_mixing_length_field(ᶜmixing_length_tuple, property))
 end
 
 """
@@ -612,6 +658,70 @@ function turbulent_prandtl_number(params, ᶜN²_eff, ᶜstrain_rate_norm)
     return min(max(prandtl_nvec, eps_FT), Pr_max)
 end
 
+
+function ᶜturbulent_prandtl_number(p)
+    (; params) = p
+    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
+    return @. lazy(
+        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm),
+    )
+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 ? p.precomputed.ᶜρa⁰ : Y.c.ρ
+
+
+    if p.atmos.turbconv_model isa PrognosticEDMFX
+        (; ᶜdetrʲs, ᶜtke⁰, ᶠ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, ᶜtke⁰, ᶠ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,
@@ -719,44 +829,33 @@ function lamb_smooth_minimum(l, smoothness_param, λ_floor)
 end
 
 """
-    eddy_viscosity(turbconv_params, tke, mixing_length)
+    ᶜeddy_viscosity(turbconv_params, tke, mixing_length)
 
 Calculates the eddy viscosity (K_u) for momentum based on the turbulent
 kinetic energy (TKE) and the mixing length.
 
 Returns K_u in units of [m^2/s].
 """
-function eddy_viscosity(turbconv_params, tke, mixing_length)
-    FT = typeof(tke)
+function ᶜeddy_viscosity(turbconv_params, tke, mixing_length)
     c_m = CAP.tke_ed_coeff(turbconv_params)
-    K_u = c_m * mixing_length * sqrt(max(tke, FT(0)))
-    return K_u
+    return @. lazy(c_m * mixing_length * sqrt(max(tke, 0)))
 end
 
 """
-    eddy_diffusivity(K_u, prandtl_nvec)
+    ᶜeddy_diffusivity(turbconv_params, tke, mixing_length, prandtl_nvec)
 
-Calculates the eddy diffusivity (K_h) for scalars given the
-eddy viscosity (K_u) and the turbulent Prandtl number.
+Calculates the eddy diffusivity (K_h) for scalars given turbulent kinetic energy (TKE),
+the mixing length, and the turbulent Prandtl number.
 
 Returns K_h in units of [m^2/s].
 """
-
-function eddy_diffusivity(K_u, prandtl_nvec)
-    return K_u / prandtl_nvec # prandtl_nvec is already bounded by eps_FT and Pr_max
+function ᶜeddy_diffusivity(turbconv_params, tke, mixing_length, prandtl_nvec)
+    K_u = ᶜeddy_viscosity(turbconv_params, tke, mixing_length)
+    return K_u / prandtl_nvec
 end
 
 
-"""
-    eddy_diffusivity(turbconv_params, tke, mixing_length, prandtl_nvec)
-
-Calculates the eddy diffusivity (K_h) for scalars given turbulent kinetic energy (TKE),
-the mixing length, and the turbulent Prandtl number.
-
-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
+function ᶜeddy_diffusivity(p, K_u)
+    ᶜprandtl_nvec = ᶜturbulent_prandtl_number(p)
+    return @. lazy(K_u / ᶜprandtl_nvec) # prandtl_nvec is already bounded by eps_FT and Pr_max
 end
diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
index 5fced40005..733878f177 100644
--- a/src/prognostic_equations/edmfx_sgs_flux.jl
+++ b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -248,10 +248,6 @@ function edmfx_sgs_mass_flux_tendency!(
         ᶜa_scalar = p.scratch.ᶜtemp_scalar
         for j in 1:n
             @. ᶠu³_diff = ᶠu³ʲs.:($$j) - ᶠu³
-            # @. ᶜa_scalar =
-            #     (ᶜmseʲs.:($$j) + ᶜKʲs.:($$j) - ᶜh_tot) *
-            #     draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j))
-            # TODO: remove this filter when mass flux is treated implicitly
             @. ᶜa_scalar =
                 (ᶜmseʲs.:($$j) + ᶜKʲs.:($$j) - ᶜh_tot) * min(
                     min(draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)), a_max),
@@ -345,17 +341,22 @@ 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⁰) = p.precomputed
     if (
         p.atmos.moisture_model isa NonEquilMoistModel &&
         p.atmos.microphysics_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}
+
+        ᶜmixing_length_field = p.scratch.ᶜtemp_scalar_2
+        ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+        ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+        ᶜK_h = ᶜeddy_diffusivity(p, ᶜK_u)
         ᶠρaK_h = p.scratch.ᶠtemp_scalar
         @. ᶠρaK_h = ᶠinterp(ᶜρa⁰) * ᶠinterp(ᶜK_h)
         ᶠρaK_u = p.scratch.ᶠtemp_scalar
@@ -382,7 +383,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
                         turbconv_params,
                         Y.c.sgs⁰.ρatke,
                         ᶜtke⁰,
-                        ᶜmixing_length,
+                        ᶜmixing_length_field,
                     ),
                     Y.c.sgs⁰.ρatke / float(dt),
                 )
@@ -450,11 +451,17 @@ 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, ᶜtke⁰, ᶜmixing_length) = p.precomputed
-    (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
+    (; ᶜu, ᶜh_tot, ᶜtke⁰) = p.precomputed
+    (; ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
 
     if p.atmos.edmfx_model.sgs_diffusive_flux isa Val{true}
+
+        ᶜmixing_length_field = p.scratch.ᶜtemp_scalar_2
+        ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+        ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+        ᶜK_h = ᶜeddy_diffusivity(p, ᶜK_u)
+
         ᶠρaK_h = p.scratch.ᶠtemp_scalar
         @. ᶠρaK_h = ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h)
         ᶠρaK_u = p.scratch.ᶠtemp_scalar
@@ -482,7 +489,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
                         turbconv_params,
                         Y.c.sgs⁰.ρatke,
                         ᶜtke⁰,
-                        ᶜmixing_length,
+                        ᶜmixing_length_field,
                     ),
                     Y.c.sgs⁰.ρatke / float(dt),
                 )
diff --git a/src/prognostic_equations/edmfx_tke.jl b/src/prognostic_equations/edmfx_tke.jl
index ae394835a9..e0091ee576 100644
--- a/src/prognostic_equations/edmfx_tke.jl
+++ b/src/prognostic_equations/edmfx_tke.jl
@@ -26,8 +26,13 @@
 edmfx_tke_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
 
 function edmfx_tke_tendency!(Yₜ, Y, p, t, turbconv_model::EDOnlyEDMFX)
-    (; ᶜstrain_rate_norm, ᶜlinear_buoygrad) = p.precomputed
-    (; ᶜK_u, ᶜK_h) = p.precomputed
+    (; ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰) = p.precomputed
+    turbconv_params = CAP.turbconv_params(p.params)
+
+    ᶜmixing_length_field = p.scratch.ᶜtemp_scalar
+    ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+    ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+    ᶜK_h = ᶜeddy_diffusivity(p, ᶜK_u)
 
     # shear production
     @. Yₜ.c.sgs⁰.ρatke += 2 * Y.c.ρ * ᶜK_u * ᶜstrain_rate_norm
@@ -45,7 +50,9 @@ 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
+    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 ?
@@ -68,6 +75,12 @@ function edmfx_tke_tendency!(
     end
 
     if use_prognostic_tke(turbconv_model)
+
+        ᶜmixing_length_field = p.scratch.ᶜtemp_scalar_2
+        ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+        ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+        ᶜK_h = ᶜeddy_diffusivity(p, ᶜK_u)
+
         # shear production
         @. Yₜ.c.sgs⁰.ρatke += 2 * ᶜρa⁰ * ᶜK_u * ᶜstrain_rate_norm
         # buoyancy production
diff --git a/src/prognostic_equations/gm_sgs_closures.jl b/src/prognostic_equations/gm_sgs_closures.jl
index 9a5b4b15b8..6460d2d5ea 100644
--- a/src/prognostic_equations/gm_sgs_closures.jl
+++ b/src/prognostic_equations/gm_sgs_closures.jl
@@ -36,7 +36,7 @@ function smagorinsky_lilly_length(c_smag, N_eff, dz, Pr, ϵ_st)
 end
 
 """
-    compute_gm_mixing_length!(ᶜmixing_length, Y, p)
+    compute_gm_mixing_length(Y, p)
 
 Computes the grid-mean subgrid-scale (SGS) mixing length using the
 Smagorinsky-Lilly formulation and stores it in `ᶜmixing_length`.
@@ -61,7 +61,7 @@ Modifies `ᶜmixing_length` in place. Also modifies fields in `p.precomputed`
 (like `ᶜlinear_buoygrad`, `ᶜstrain_rate_norm`) and uses `p.scratch` for
 intermediate calculations.
 """
-NVTX.@annotate function compute_gm_mixing_length!(ᶜmixing_length, Y, p)
+NVTX.@annotate function compute_gm_mixing_length(Y, p)
     (; params) = p
     thermo_params = CAP.thermodynamics_params(params)
 
@@ -92,11 +92,13 @@ NVTX.@annotate function compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     @. ᶜprandtl_nvec =
         turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
 
-    @. ᶜmixing_length = smagorinsky_lilly_length(
-        CAP.c_smag(params),
-        sqrt(max(ᶜlinear_buoygrad, 0)),   # N_eff
-        ᶜdz,
-        ᶜprandtl_nvec,
-        ᶜstrain_rate_norm,
+    return @. lazy(
+        smagorinsky_lilly_length(
+            CAP.c_smag(params),
+            sqrt(max(ᶜlinear_buoygrad, 0)),   # N_eff
+            ᶜdz,
+            ᶜprandtl_nvec,
+            ᶜstrain_rate_norm,
+        ),
     )
 end
diff --git a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
index 872a7185ea..16fe44db66 100644
--- a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
+++ b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
@@ -544,8 +544,30 @@ 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)
+        (; vertical_diffusion) = p.atmos
+        (; ᶜp) = p.precomputed
+        if vertical_diffusion isa DecayWithHeightDiffusion
+            ᶜK_h =
+                compute_eddy_diffusivity_coefficient(Y.c.ρ, vertical_diffusion)
+            ᶜK_u = ᶜK_h
+        elseif vertical_diffusion isa VerticalDiffusion
+            ᶜK_h = compute_eddy_diffusivity_coefficient(
+                Y.c.uₕ,
+                ᶜp,
+                vertical_diffusion,
+            )
+            ᶜK_u = ᶜK_h
+        else
+            (; ᶜtke⁰,) = p.precomputed
+            ᶜmixing_length_field = p.scratch.ᶜtemp_scalar_3
+            ᶜmixing_length_field .= ᶜmixing_length(Y, p)
+            ᶜK_u = ᶜeddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+            ᶜK_h = ᶜeddy_diffusivity(p, ᶜK_u)
+        end
+
         α_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()
@@ -613,12 +635,14 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
             turbconv_params = CAP.turbconv_params(params)
             c_d = CAP.tke_diss_coeff(turbconv_params)
             (; dt) = p
-            (; ᶜtke⁰, ᶜmixing_length) = p.precomputed
+            (; ᶜtke⁰) = p.precomputed
             ᶜρa⁰ =
                 p.atmos.turbconv_model isa PrognosticEDMFX ?
                 p.precomputed.ᶜρa⁰ : ᶜρ
             ᶜρatke⁰ = Y.c.sgs⁰.ρatke
 
+            ᶜmixing_length_field = ᶜmixing_length(Y, p)
+
             @inline tke_dissipation_rate_tendency(tke⁰, mixing_length) =
                 tke⁰ >= 0 ? c_d * sqrt(tke⁰) / mixing_length : 1 / float(dt)
             @inline ∂tke_dissipation_rate_tendency_∂tke⁰(tke⁰, mixing_length) =
@@ -627,8 +651,10 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
 
             ᶜdissipation_matrix_diagonal = p.scratch.ᶜtemp_scalar
             @. ᶜdissipation_matrix_diagonal =
-                ᶜρatke⁰ *
-                ∂tke_dissipation_rate_tendency_∂tke⁰(ᶜtke⁰, ᶜmixing_length)
+                ᶜρatke⁰ * ∂tke_dissipation_rate_tendency_∂tke⁰(
+                    ᶜtke⁰,
+                    ᶜmixing_length_field,
+                )
 
             ∂ᶜρatke⁰_err_∂ᶜρ = matrix[@name(c.sgs⁰.ρatke), @name(c.ρ)]
             ∂ᶜρatke⁰_err_∂ᶜρatke⁰ =
@@ -644,7 +670,10 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
                         ᶜdiffusion_u_matrix -
                         DiagonalMatrixRow(ᶜdissipation_matrix_diagonal)
                     ) ⋅ DiagonalMatrixRow(1 / ᶜρa⁰) - DiagonalMatrixRow(
-                        tke_dissipation_rate_tendency(ᶜtke⁰, ᶜmixing_length),
+                        tke_dissipation_rate_tendency(
+                            ᶜtke⁰,
+                            ᶜmixing_length_field,
+                        ),
                     )
                 ) - (I,)
         end
diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
index bca3cd65c2..98aa10a7ca 100644
--- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
+++ b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
@@ -81,16 +81,24 @@ function vertical_diffusion_boundary_layer_tendency!(
     ::Union{VerticalDiffusion, DecayWithHeightDiffusion},
 )
     FT = eltype(Y)
+    (; vertical_diffusion) = p.atmos
     α_vert_diff_tracer = CAP.α_vert_diff_tracer(p.params)
-    (; ᶜu, ᶜh_tot, ᶜspecific, ᶜK_u, ᶜK_h) = p.precomputed
+    (; ᶜu, ᶜh_tot, ᶜspecific, ᶜp) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
 
+    if vertical_diffusion isa DecayWithHeightDiffusion
+        ᶜK_h = compute_eddy_diffusivity_coefficient(Y.c.ρ, vertical_diffusion)
+    elseif vertical_diffusion isa VerticalDiffusion
+        ᶜK_h =
+            compute_eddy_diffusivity_coefficient(Y.c.uₕ, ᶜp, vertical_diffusion)
+    end
+
     if !disable_momentum_vertical_diffusion(p.atmos.vertical_diffusion)
         ᶠstrain_rate = p.scratch.ᶠtemp_UVWxUVW
         ᶠstrain_rate .= compute_strain_rate_face(ᶜu)
         @. Yₜ.c.uₕ -= C12(
-            ᶜdivᵥ(-2 * ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_u) * ᶠstrain_rate) / Y.c.ρ,
-        )
+            ᶜdivᵥ(-2 * ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠstrain_rate) / Y.c.ρ,
+        ) # assumes ᶜK_u = ᶜK_h
     end
 
     ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
diff --git a/src/solver/types.jl b/src/solver/types.jl
index 3e2fbb8388..eb38db31a5 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -346,6 +346,11 @@ struct MixingLength{FT}
     l_grid::FT
 end
 
+function MixingLength(master, wall, tke, buoy, l_grid)
+    return MixingLength(promote(master, wall, tke, buoy, l_grid)...)
+end
+
+
 abstract type AbstractEDMF end
 
 """