diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 9835b817ad..1e1301a682 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -267,12 +267,6 @@ function precomputed_quantities(Y, atmos)
     )
 end
 
-# Interpolates the third contravariant component of Y.c.uₕ to cell faces.
-function compute_ᶠuₕ³(ᶜuₕ, ᶜρ)
-    ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    return @. lazy(ᶠwinterp(ᶜρ * ᶜJ, CT3(ᶜuₕ)))
-end
-
 """
     set_velocity_at_surface!(Y, ᶠuₕ³, turbconv_model)
 
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index a932682049..fb7e2c7323 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -174,11 +174,11 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
     (; ᶜf³, ᶠf¹², ᶜΦ) = p.core
+    #TODO Test new compute_<> utility functions (set surface boundary condition) and rm from precomputed quantities
     (; ᶜu, ᶠu³, ᶜK) = p.precomputed
     (; edmfx_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing
     (; ᶜuʲs, ᶜKʲs, ᶠKᵥʲs) = n > 0 ? p.precomputed : all_nothing
     (; energy_upwinding, tracer_upwinding) = p.atmos.numerics
-    (; ᶜspecific) = p.precomputed
 
     ᶠu³⁰ =
         advect_tke ?
@@ -191,22 +191,20 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     ᶜtke⁰ = advect_tke ? p.precomputed.ᶜtke⁰ : nothing
     ᶜa_scalar = p.scratch.ᶜtemp_scalar
     ᶜω³ = p.scratch.ᶜtemp_CT3
-    ᶠω¹² = p.scratch.ᶠtemp_CT12
-    ᶠω¹²ʲs = p.scratch.ᶠtemp_CT12ʲs
 
     if point_type <: Geometry.Abstract3DPoint
-        @. ᶜω³ = curlₕ(Y.c.uₕ)
+        ᶜω³ = @. lazy(curlₕ(Y.c.uₕ))
     elseif point_type <: Geometry.Abstract2DPoint
-        @. ᶜω³ = zero(ᶜω³)
+        @. ᶜω³ = @. lazy(zero(ᶜω³))
     end
 
-    @. ᶠω¹² = ᶠcurlᵥ(Y.c.uₕ)
+    ᶠω¹² = @. lazy(ᶠcurlᵥ(Y.c.uₕ))
     for j in 1:n
-        @. ᶠω¹²ʲs.:($$j) = ᶠω¹²
+        ᶠω¹²ʲs.:($j) = @. lazy(ᶠω¹²)
     end
-    @. ᶠω¹² += CT12(curlₕ(Y.f.u₃))
+    ᶠω¹² += @. lazy(CT12(curlₕ(Y.f.u₃)))
     for j in 1:n
-        @. ᶠω¹²ʲs.:($$j) += CT12(curlₕ(Y.f.sgsʲs.:($$j).u₃))
+        @. ᶠω¹²ʲs.:($j) += CT12(curlₕ(Y.f.sgsʲs.:($j).u₃))
     end
     # Without the CT12(), the right-hand side would be a CT1 or CT2 in 2D space.
 
diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl
index be779d6418..58449493bb 100644
--- a/src/utils/utilities.jl
+++ b/src/utils/utilities.jl
@@ -47,6 +47,27 @@ determined by the filename.
 sort_files_by_time(files) =
     permute!(files, sortperm(time_from_filename.(files)))
 
+
+"""
+    Velocity manipulations
+"""
+function compute_ᶠuₕ³(Y)
+    return compute_ᶠuₕ³(Y.c.uₕ, Y.c.ρ)
+end
+
+function compute_ᶠuₕ³(ᶜuₕ, ᶜρ)
+    ᶜJ = Fields.local_geometry_field(ᶜρ).J
+    return @. lazy(ᶠwinterp(ᶜρ * ᶜJ, CT3(ᶜuₕ)))
+end
+
+function compute_ᶠu³(Y)
+    return compute_ᶠuₕ³(Y) + CT3(Y.f.u₃)
+end
+
+function compute_ᶜu(Y)
+    return C123(Y.c.uₕ) + ᶜinterp(C123(Y.f.u₃))
+end
+
 """
     κ .= compute_kinetic(uₕ::Field, uᵥ::Field)