diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index de6e795cba6..f7010e1955e 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -651,16 +651,14 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
         ),
     )
 
-    # TODO: This is not correct with topography or with grid stretching
+    # TODO: Currently the shear production only includes vertical gradients
     ᶠu⁰ = p.ᶠtemp_C123
     @. ᶠu⁰ = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³⁰)
-    ∇ᵥuvw = p.ᶜtemp_UVWxUVW
-    ∇ᵥuvw = Geometry.project.(Ref(Geometry.UVWAxis()), ᶜgradᵥ.(UVW.(ᶠu⁰)))
-    @. ᶜshear² = norm_sqr((∇ᵥuvw + adjoint(∇ᵥuvw)) / 2)
+    ᶜgradᵥ_ᶠuvw⁰ = p.ᶜtemp_UVWxUVW
+    ᶜgradᵥ_ᶠuvw⁰ =
+        Geometry.project.(Ref(Geometry.UVWAxis()), ᶜgradᵥ.(UVW.(ᶠu⁰)))
+    @. ᶜshear² = norm_sqr((ᶜgradᵥ_ᶠuvw⁰ + adjoint(ᶜgradᵥ_ᶠuvw⁰)) / 2)
 
-    # ct3_unit = p.ᶜtemp_CT3
-    # @. ct3_unit = CT3(Geometry.WVector(FT(1)), ᶜlg)
-    # @. ᶜshear² = norm_sqr(adjoint(CA.ᶜgradᵥ(ᶠu⁰)) * ct3_unit)
     ᶜprandtl_nvec = p.ᶜtemp_scalar
     @. ᶜprandtl_nvec = turbulent_prandtl_number(
         params,
diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl
index 5c9c6bf3eec..41a27e8d314 100644
--- a/src/cache/temporary_quantities.jl
+++ b/src/cache/temporary_quantities.jl
@@ -51,7 +51,10 @@ function temporary_quantities(atmos, center_space, face_space)
         ᶠtemp_CT12ʲs = Fields.Field(NTuple{n, CT12{FT}}, face_space), # ᶠω¹²ʲs
         ᶜtemp_C123 = Fields.Field(C123{FT}, center_space), # χ₁₂₃
         ᶠtemp_C123 = Fields.Field(C123{FT}, face_space), # χ₁₂₃
-        ᶜtemp_UVWxUVW = Fields.Field(typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'), center_space), # ∇ᵥuvw
+        ᶜtemp_UVWxUVW = Fields.Field(
+            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
+            center_space,
+        ), # ∇ᵥuvw
         # TODO: Remove this hack
         sfc_temp_C3 = Fields.Field(C3{FT}, Spaces.level(face_space, half)), # ρ_flux_χ
     )