diff --git a/Project.toml b/Project.toml
index 39ea9cd692..0f0bc54fd2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
diff --git a/src/Forcing.jl b/src/Forcing.jl
index b429187510..9638e0f06a 100644
--- a/src/Forcing.jl
+++ b/src/Forcing.jl
@@ -58,8 +58,8 @@ function update(self::ForcingBase{ForcingStandard}, GMV::GridMeanVariables)
     if self.apply_subsidence
         @inbounds for k in real_center_indicies(self.Gr)
             # Apply large-scale subsidence tendencies
-            GMV.H.tendencies[k] -= (GMV.H.values[k + 1] - GMV.H.values[k]) * self.Gr.dzi * self.subsidence[k]
-            GMV.QT.tendencies[k] -= (GMV.QT.values[k + 1] - GMV.QT.values[k]) * self.Gr.dzi * self.subsidence[k]
+            GMV.H.tendencies[k] -= ∇_upwind(GMV.H.values, self.Gr, k) * self.subsidence[k]
+            GMV.QT.tendencies[k] -= ∇_upwind(GMV.QT.values, self.Gr, k) * self.subsidence[k]
         end
     end
     if self.apply_coriolis
@@ -85,8 +85,8 @@ function update(self::ForcingBase{ForcingDYCOMS_RF01}, GMV::GridMeanVariables)
     @inbounds for k in real_center_indicies(self.Gr)
         GMV.QT.tendencies[k] += self.dqtdt[k]
         # Apply large-scale subsidence tendencies
-        GMV.H.tendencies[k] -= (GMV.H.values[k + 1] - GMV.H.values[k]) * self.Gr.dzi * self.subsidence[k]
-        GMV.QT.tendencies[k] -= (GMV.QT.values[k + 1] - GMV.QT.values[k]) * self.Gr.dzi * self.subsidence[k]
+        GMV.H.tendencies[k] -= ∇_upwind(GMV.H.values, self.Gr, k) * self.subsidence[k]
+        GMV.QT.tendencies[k] -= ∇_upwind(GMV.QT.values, self.Gr, k) * self.subsidence[k]
     end
 
     if self.apply_coriolis
diff --git a/src/Operators.jl b/src/Operators.jl
new file mode 100644
index 0000000000..c2d6946b8b
--- /dev/null
+++ b/src/Operators.jl
@@ -0,0 +1,111 @@
+
+abstract type AbstractBCTag end
+struct BottomBCTag <: AbstractBCTag end
+struct TopBCTag <: AbstractBCTag end
+struct InteriorTag <: AbstractBCTag end
+
+struct NoBCGivenError end
+struct SetValue{FT}
+    value::FT
+end
+struct SetGradient{FT}
+    value::FT
+end
+
+function ∇f2c(f, grid::Grid, k::Int)
+    return (f[k] - f[k - 1]) * grid.dzi
+end
+
+function ∇c2f(f, grid::Grid, k::Int)
+    return (f[k + 1] - f[k]) * grid.dzi
+end
+
+function ∇_upwind(f, grid::Grid, k::Int)
+    return (f[k + 1] - f[k]) * grid.dzi
+end
+
+function interpc2f(f, grid::Grid, k::Int)
+    return 0.5 * (f[k + 1] + f[k])
+end
+
+#####
+##### ∇(center data)
+#####
+
+c∇(f, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) = c∇(cut(f, k), grid, k; bottom, top)
+
+function c∇(f_cut::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
+    gw = grid.gw
+    if k == gw # NOTE: 0-based indexing (k = 3 for 1-based indexing) bottom boundary
+        return c∇(f_cut, grid, BottomBCTag(), bottom)
+    elseif k == grid.nzg - gw - 1 # NOTE: 0-based indexing
+        return c∇(f_cut, grid, TopBCTag(), top)
+    else
+        return c∇(f_cut, grid, InteriorTag())
+    end
+end
+c∇(f::SVector, grid::Grid, ::AbstractBCTag, ::NoBCGivenError) = error("No BC given")
+function c∇(f::SVector, grid::Grid, ::InteriorTag)
+    @assert length(f) == 3
+    f_dual⁺ = SVector(f[2], f[3])
+    f_dual⁻ = SVector(f[1], f[2])
+    return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
+end
+function c∇(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue)
+    @assert length(f) == 3
+    # 2fb = cg+ci => cg = 2fb-ci
+    f_dual⁺ = SVector(f[2], 2 * bc.value - f[2])
+    f_dual⁻ = SVector(f[1], f[2])
+    return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
+end
+function c∇(f::SVector, grid::Grid, ::BottomBCTag, bc::SetValue)
+    @assert length(f) == 3
+    # 2fb = cg+ci => cg = 2fb-ci
+    f_dual⁺ = SVector(f[2], f[3])
+    f_dual⁻ = SVector(2 * bc.value - f[2], f[2])
+    return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
+end
+function c∇(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient)
+    @assert length(f) == 3
+    f_dual⁺ = SVector(f[2], f[3])
+    f_dual⁻ = SVector(f[1], f[2])
+    return (bc.value + ∇_staggered(f_dual⁻, grid)) / 2
+end
+function c∇(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient)
+    @assert length(f) == 3
+    f_dual⁺ = SVector(f[2], f[3])
+    f_dual⁻ = SVector(f[1], f[2])
+    return (∇_staggered(f_dual⁺, grid) + bc.value) / 2
+end
+
+#####
+##### Generic functions
+#####
+
+function ∇_staggered(f::SVector, grid::Grid)
+    @assert length(f) == 2
+    return (f[2] - f[1]) * grid.dzi
+end
+
+# A 3-point field stencil
+function cut(f::AbstractVector, k::Int)
+    return SVector(f[k - 1], f[k], f[k + 1])
+end
+
+# A 3-point field stencil for updraft variables
+function cut(f::AbstractMatrix, k::Int, i_up::Int)
+    return SVector(f[i_up, k - 1], f[i_up, k], f[i_up, k + 1])
+end
+
+function ∇_collocated(f::SVector, grid::Grid)
+    @assert length(f) == 3
+    return (f[3] - f[1]) * grid.dzi
+end
+
+# TODO: use this implementation
+# function ∇_collocated(f::SVector, grid::Grid)
+#     @assert length(f) == 3
+#     f_dual⁺ = SVector(f[2], f[3])
+#     f_dual⁻ = SVector(f[1], f[2])
+#     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
+# end
diff --git a/src/Turbulence.jl b/src/Turbulence.jl
index ceb918620b..7843ccb6c1 100644
--- a/src/Turbulence.jl
+++ b/src/Turbulence.jl
@@ -27,7 +27,7 @@ end
 # Update the diagnosis of the inversion height, using the maximum temperature gradient method
 function update_inversion(self::ParameterizationBase, GMV::GridMeanVariables, option)
     theta_rho = center_field(self.Gr)
-    maxgrad = 0.0
+    ∇θ_liq_max = 0.0
     k_fi = first_center(self.Gr)
 
     @inbounds for k in real_center_indicies(self.Gr)
@@ -46,9 +46,9 @@ function update_inversion(self::ParameterizationBase, GMV::GridMeanVariables, op
     elseif option == "thetal_maxgrad"
 
         @inbounds for k in real_center_indicies(self.Gr)
-            grad = (GMV.THL.values[k + 1] - GMV.THL.values[k]) * self.Gr.dzi
-            if grad > maxgrad
-                maxgrad = grad
+            ∇θ_liq = ∇_upwind(GMV.THL.values, self.Gr, k)
+            if ∇θ_liq > ∇θ_liq_max
+                ∇θ_liq_max = ∇θ_liq
                 self.zi = self.Gr.z[k]
             end
         end
diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl
index 778f99d7e1..ad021dd299 100644
--- a/src/TurbulenceConvection.jl
+++ b/src/TurbulenceConvection.jl
@@ -2,6 +2,7 @@ module TurbulenceConvection
 
 import Dierckx
 using StatsBase
+using StaticArrays
 import Statistics
 using PoissonRandom: pois_rand
 import LambertW
@@ -33,6 +34,7 @@ include("python_primitives.jl")
 include("parameters.jl")
 include("Grid.jl")
 include("Fields.jl")
+include("Operators.jl")
 include("NetCDFIO.jl")
 include("ReferenceState.jl")
 include("types.jl")
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index 4090398697..197eb712f7 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -410,10 +410,12 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
         end
 
         # Buoyancy-shear-subdomain exchange-dissipation TKE equilibrium scale
+        U_cut = cut(GMV.U.values, k)
+        V_cut = cut(GMV.V.values, k)
         shear2 =
-            pow((GMV.U.values[k + 1] - GMV.U.values[k - 1]) * 0.5 * dzi, 2) +
-            pow((GMV.V.values[k + 1] - GMV.V.values[k - 1]) * 0.5 * dzi, 2) +
-            pow((self.EnvVar.W.values[k] - self.EnvVar.W.values[k - 1]) * dzi, 2)
+            pow(∇_collocated(U_cut, grid), 2) +
+            pow(∇_collocated(V_cut, grid), 2) +
+            pow(∇f2c(self.EnvVar.W.values, grid, k), 2)
         qt_dry = self.EnvThermo.qt_dry[k]
         th_dry = self.EnvThermo.th_dry[k]
         t_cloudy = self.EnvThermo.t_cloudy[k]