diff --git a/src/Forcing.jl b/src/Forcing.jl
index 3c39e0ff36..f82f142a2a 100644
--- a/src/Forcing.jl
+++ b/src/Forcing.jl
@@ -50,8 +50,8 @@ function update(self::ForcingBase{ForcingStandard}, GMV::GridMeanVariables)
             # Apply large-scale subsidence tendencies
             H_cut = ccut_onesided(GMV.H.values, grid, k)
             q_tot_cut = ccut_onesided(GMV.QT.values, grid, k)
-            ∇H = ∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
-            ∇q_tot = ∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+            ∇H = c∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+            ∇q_tot = c∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
             GMV.H.tendencies[k] -= ∇H * self.subsidence[k]
             GMV.QT.tendencies[k] -= ∇q_tot * self.subsidence[k]
         end
@@ -79,8 +79,8 @@ function update(self::ForcingBase{ForcingDYCOMS_RF01}, GMV::GridMeanVariables)
     @inbounds for k in real_center_indicies(grid)
         H_cut = ccut_onesided(GMV.H.values, grid, k)
         q_tot_cut = ccut_onesided(GMV.QT.values, grid, k)
-        ∇H = ∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
-        ∇q_tot = ∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+        ∇H = c∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+        ∇q_tot = c∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
 
         GMV.QT.tendencies[k] += self.dqtdt[k]
         # Apply large-scale subsidence tendencies
@@ -140,8 +140,8 @@ function update(self::ForcingBase{ForcingLES}, GMV::GridMeanVariables)
 
             H_cut = ccut_onesided(GMV.H.values, grid, k)
             q_tot_cut = ccut_onesided(GMV.QT.values, grid, k)
-            ∇H = ∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
-            ∇q_tot = ∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+            ∇H = c∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+            ∇q_tot = c∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
             GMV_H_subsidence_k = -∇H * self.subsidence[k]
             GMV_QT_subsidence_k = -∇q_tot * self.subsidence[k]
         else
diff --git a/src/Operators.jl b/src/Operators.jl
index 61c8635074..091915c7da 100644
--- a/src/Operators.jl
+++ b/src/Operators.jl
@@ -35,26 +35,41 @@ end
 ∇f2c(f::SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[2] - bc.value) * grid.dzi
 ∇f2c(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
 
-function ∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
+function c∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
     if is_surface_center(grid, k)
-        return ∇_onesided(f_dual, grid, BottomBCTag(), bottom)
+        return c∇_onesided(f_dual, grid, BottomBCTag(), bottom)
     elseif is_toa_center(grid, k)
-        return ∇_onesided(f_dual, grid, TopBCTag(), top)
+        return c∇_onesided(f_dual, grid, TopBCTag(), top)
     else
-        return ∇_onesided(f_dual, grid, InteriorTag())
+        return c∇_onesided(f_dual, grid, InteriorTag())
     end
 end
-∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi
-∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2)
+c∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi
+c∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2)
 # TODO: this is a crud approximation, as we're specifying what should be the derivative
 # at the boundary, and we're taking this as the derivative at the first interior at the
 # top of the domain.
-∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value
-∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info
+c∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value
+c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info
 # TODO: this is a crud approximation, as we're specifying what should be the derivative
 # at the boundary, and we're taking this as the derivative at the first interior at the
 # top of the domain.
-∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
+c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
+
+function f∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
+    if is_surface_face(grid, k)
+        return f∇_onesided(f_dual, grid, BottomBCTag(), bottom)
+    elseif is_toa_face(grid, k)
+        return f∇_onesided(f_dual, grid, TopBCTag(), top)
+    else
+        return f∇_onesided(f_dual, grid, InteriorTag())
+    end
+end
+f∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi
+f∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2)
+f∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value
+f∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info
+f∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
 
 # Used when traversing cell faces
 
@@ -111,7 +126,7 @@ function upwind_advection_area(ρ0_half::Vector{Float64}, a_up::Vector{Float64},
     ρ_0_cut = ccut_onesided(ρ0_half, grid, k)
     a_up_cut = ccut_onesided(a_up, grid, k)
     m_cut = ρ_0_cut .* a_up_cut .* w_up_cut
-    ∇m = ∇_onesided(m_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+    ∇m = c∇_onesided(m_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
     # TODO: Why are we dividing by ρ0_half[k + 1]?
     return -∇m / ρ0_half[k + 1]
 end
@@ -124,6 +139,16 @@ function upwind_advection_velocity(ρ0::Vector{Float64}, a_up::Vector{Float64},
     return adv
 end
 
+function upwind_advection_velocity(ρ0::Vector{Float64}, a_up::Vector{Float64}, w_up::Vector{Float64}, grid, k)
+    a_dual = fdaul_onesided(a_up, grid, k)
+    ρ_0_dual = dual_faces(ρ0, grid, k)
+    w_up_dual = dual_faces(w_up, grid, k)
+    adv_dual = a_dual .* ρ_0_dual .* w_up_dual .* w_up_dual
+    ∇ρaw = f∇_onesided(adv_dual, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+    return ∇ρaw
+end
+
+
 function upwind_advection_scalar(
     ρ0_half::Vector{Float64},
     a_up::Vector{Float64},
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index 2a967f625e..cfcf6f5c93 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -1837,7 +1837,7 @@ function update_inversion(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, opti
 
         @inbounds for k in real_center_indicies(grid)
             ∇θ_liq_cut = ccut_onesided(GMV.H.values, grid, k)
-            ∇θ_liq = ∇_onesided(∇θ_liq_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+            ∇θ_liq = c∇_onesided(∇θ_liq_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
             if ∇θ_liq > ∇θ_liq_max
                 ∇θ_liq_max = ∇θ_liq
                 self.zi = grid.z[k]