diff --git a/docs/src/operators.md b/docs/src/operators.md
index 217a46e460..a61e2460d0 100644
--- a/docs/src/operators.md
+++ b/docs/src/operators.md
@@ -91,8 +91,8 @@ CurlC2F
 
 ```@docs
 SetBoundaryOperator
-FirstOrderOneSided
-ThirdOrderOneSided
+OneSided1stOrder
+OneSided3rdOrder
 ```
 
 ## Finite difference boundary conditions
diff --git a/examples/column/bb_fct_advection.jl b/examples/column/bb_fct_advection.jl
index 746ed3792e..33af59f976 100644
--- a/examples/column/bb_fct_advection.jl
+++ b/examples/column/bb_fct_advection.jl
@@ -48,12 +48,12 @@ function tendency!(yₜ, y, parameters, t)
         top = Operators.Extrapolate(),
     )
     upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.ThirdOrderOneSided(),
-        top = Operators.ThirdOrderOneSided(),
+        bottom = Operators.OneSided3rdOrder(),
+        top = Operators.OneSided3rdOrder(),
     )
     FCTBorisBook = Operators.FCTBorisBook(
-        bottom = Operators.FirstOrderOneSided(),
-        top = Operators.FirstOrderOneSided(),
+        bottom = Operators.OneSided1stOrder(),
+        top = Operators.OneSided1stOrder(),
     )
     @. yₜ.q =
         -divf2c(
diff --git a/examples/column/fct_advection.jl b/examples/column/fct_advection.jl
index 8b047367f0..86495ed5a5 100644
--- a/examples/column/fct_advection.jl
+++ b/examples/column/fct_advection.jl
@@ -45,8 +45,8 @@ function f!(dydt, y, parameters, t)
         top = Operators.Extrapolate(),
     )
     third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.ThirdOrderOneSided(),
-        top = Operators.ThirdOrderOneSided(),
+        bottom = Operators.OneSided3rdOrder(),
+        top = Operators.OneSided3rdOrder(),
     )
     divf2c = Operators.DivergenceF2C(
         bottom = Operators.SetValue(Geometry.WVector(FT(0.0))),
diff --git a/examples/column/zalesak_fct_advection.jl b/examples/column/zalesak_fct_advection.jl
index efc20fbf44..5ec508c8ac 100644
--- a/examples/column/zalesak_fct_advection.jl
+++ b/examples/column/zalesak_fct_advection.jl
@@ -48,12 +48,12 @@ function tendency!(yₜ, y, parameters, t)
         top = Operators.Extrapolate(),
     )
     upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.ThirdOrderOneSided(),
-        top = Operators.ThirdOrderOneSided(),
+        bottom = Operators.OneSided3rdOrder(),
+        top = Operators.OneSided3rdOrder(),
     )
     FCTZalesak = Operators.FCTZalesak(
-        bottom = Operators.FirstOrderOneSided(),
-        top = Operators.FirstOrderOneSided(),
+        bottom = Operators.OneSided1stOrder(),
+        top = Operators.OneSided1stOrder(),
     )
     @. yₜ.q =
         -divf2c(
diff --git a/examples/hybrid/box/limiters_advection.jl b/examples/hybrid/box/limiters_advection.jl
index daaa8155a7..046ed58cd2 100644
--- a/examples/hybrid/box/limiters_advection.jl
+++ b/examples/hybrid/box/limiters_advection.jl
@@ -179,8 +179,8 @@ function vertical_tendency!(yₜ, y, cache, t)
         bottom = Operators.SetValue(Geometry.Contravariant3Vector(FT(0))),
     )
     upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.ThirdOrderOneSided(),
-        top = Operators.ThirdOrderOneSided(),
+        bottom = Operators.OneSided3rdOrder(),
+        top = Operators.OneSided3rdOrder(),
     )
     ax12 = (Geometry.Covariant12Axis(),)
     ax3 = (Geometry.Covariant3Axis(),)
diff --git a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
index 6b3d1107fe..70a68337c4 100644
--- a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
+++ b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
@@ -172,8 +172,8 @@ function rhs_invariant!(dY, Y, _, t)
     )
     # 1.c) vertical upwinding
     third_order_upwind_c2f = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.ThirdOrderOneSided(),
-        top = Operators.ThirdOrderOneSided(),
+        bottom = Operators.OneSided3rdOrder(),
+        top = Operators.OneSided3rdOrder(),
     )
     # we want the total u³ at the boundary to be zero: we can either constrain
     # both to be zero, or allow one to be non-zero and set the other to be its
diff --git a/examples/hybrid/plane/topo_schar_nh.jl b/examples/hybrid/plane/topo_schar_nh.jl
index d4277ef37a..64181f523e 100644
--- a/examples/hybrid/plane/topo_schar_nh.jl
+++ b/examples/hybrid/plane/topo_schar_nh.jl
@@ -230,8 +230,8 @@ function rhs_invariant!(dY, Y, _, t)
 
     f_upwind_product1 = Operators.UpwindBiasedProductC2F()
     f_upwind_product3 = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.FirstOrderOneSided(),
-        top = Operators.FirstOrderOneSided(),
+        bottom = Operators.OneSided1stOrder(),
+        top = Operators.OneSided1stOrder(),
     )
 
     dρ .= 0 .* cρ
diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl
index db995ae2f8..3697ba667e 100644
--- a/examples/hybrid/sphere/deformation_flow.jl
+++ b/examples/hybrid/sphere/deformation_flow.jl
@@ -67,16 +67,16 @@ const upwind1 = Operators.UpwindBiasedProductC2F(
     top = Operators.Extrapolate(),
 )
 const upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
-    bottom = Operators.ThirdOrderOneSided(),
-    top = Operators.ThirdOrderOneSided(),
+    bottom = Operators.OneSided3rdOrder(),
+    top = Operators.OneSided3rdOrder(),
 )
 const FCTZalesak = Operators.FCTZalesak(
-    bottom = Operators.FirstOrderOneSided(),
-    top = Operators.FirstOrderOneSided(),
+    bottom = Operators.OneSided1stOrder(),
+    top = Operators.OneSided1stOrder(),
 )
 const FCTBorisBook = Operators.FCTBorisBook(
-    bottom = Operators.FirstOrderOneSided(),
-    top = Operators.FirstOrderOneSided(),
+    bottom = Operators.OneSided1stOrder(),
+    top = Operators.OneSided1stOrder(),
 )
 
 # Reference pressure and density
diff --git a/examples/hybrid/sphere/hadley_circulation.jl b/examples/hybrid/sphere/hadley_circulation.jl
index 949b2bf9c5..597ad6c4bb 100644
--- a/examples/hybrid/sphere/hadley_circulation.jl
+++ b/examples/hybrid/sphere/hadley_circulation.jl
@@ -94,12 +94,12 @@ function tendency!(yₜ, y, parameters, t)
         top = Operators.Extrapolate(),
     )
     upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
-        bottom = Operators.ThirdOrderOneSided(),
-        top = Operators.ThirdOrderOneSided(),
+        bottom = Operators.OneSided3rdOrder(),
+        top = Operators.OneSided3rdOrder(),
     )
     FCTZalesak = Operators.FCTZalesak(
-        bottom = Operators.FirstOrderOneSided(),
-        top = Operators.FirstOrderOneSided(),
+        bottom = Operators.OneSided1stOrder(),
+        top = Operators.OneSided1stOrder(),
     )
     hdiv = Operators.Divergence()
     hwdiv = Operators.WeakDivergence()
diff --git a/examples/hybrid/staggered_nonhydrostatic_model.jl b/examples/hybrid/staggered_nonhydrostatic_model.jl
index 302de7696f..1879104818 100644
--- a/examples/hybrid/staggered_nonhydrostatic_model.jl
+++ b/examples/hybrid/staggered_nonhydrostatic_model.jl
@@ -52,8 +52,8 @@ const ᶜFC = Operators.FluxCorrectionC2C(
 )
 const ᶠupwind_product1 = Operators.UpwindBiasedProductC2F()
 const ᶠupwind_product3 = Operators.Upwind3rdOrderBiasedProductC2F(
-    bottom = Operators.ThirdOrderOneSided(),
-    top = Operators.ThirdOrderOneSided(),
+    bottom = Operators.OneSided3rdOrder(),
+    top = Operators.OneSided3rdOrder(),
 )
 
 const ᶜinterp_stencil = Operators.Operator2Stencil(ᶜinterp)
diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl
index 973d59b842..0066f87b3d 100644
--- a/src/Operators/Operators.jl
+++ b/src/Operators/Operators.jl
@@ -20,7 +20,9 @@ using ..RecursiveApply
 include("common.jl")
 include("spectralelement.jl")
 include("numericalflux.jl")
-include("finitedifference.jl")
+include("finitedifference/finitedifference.jl")
+include("finitedifference/upwinding.jl")
+include("finitedifference/deprecated.jl")
 include("stencilcoefs.jl")
 include("operator2stencil.jl")
 include("pointwisestencil.jl")
diff --git a/src/Operators/finitedifference/deprecated.jl b/src/Operators/finitedifference/deprecated.jl
new file mode 100644
index 0000000000..1c847425e4
--- /dev/null
+++ b/src/Operators/finitedifference/deprecated.jl
@@ -0,0 +1,534 @@
+# to be deprecated in future
+
+# these were renamed
+const LeftBiasedC2F = LeftBiased1stOrderC2F
+const RightBiasedC2F = RightBiased1stOrderC2F
+const FirstOrderOneSided = OneSided1stOrder
+const ThirdOrderOneSided = OneSided3rdOrder
+
+# F2C aren't used
+
+"""
+    L = LeftBiasedF2C(;boundaries)
+    L.(x)
+
+Interpolate a face-value field to a center-valued field from the left.
+```math
+L(x)[i+\\tfrac{1}{2}] = x[i]
+```
+
+Only the left boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+L(x)[1] = x_0
+```
+"""
+struct LeftBiasedF2C{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+LeftBiasedF2C(; kwargs...) = LeftBiasedF2C(NamedTuple(kwargs))
+
+return_space(::LeftBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) =
+    Spaces.CenterFiniteDifferenceSpace(space)
+return_space(::LeftBiasedF2C, space::Spaces.FaceExtrudedFiniteDifferenceSpace) =
+    Spaces.CenterExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::LeftBiasedF2C, arg) = ((-half, -half),)
+Base.@propagate_inbounds stencil_interior(
+    ::LeftBiasedF2C,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+) = getidx(space, arg, loc, idx - half, hidx)
+left_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiasedF2C,
+    ::AbstractBoundaryCondition,
+    arg,
+) = left_idx(space)
+right_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiasedF2C,
+    ::AbstractBoundaryCondition,
+    arg,
+) = right_idx(space)
+
+left_interior_idx(space::AbstractSpace, ::LeftBiasedF2C, ::SetValue, arg) =
+    left_idx(space) + 1
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::LeftBiasedF2C,
+    bc::SetValue,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == left_center_boundary_idx(space)
+    getidx(space, bc.val, loc, nothing, hidx)
+end
+
+
+
+
+"""
+    L = LeftBiased3rdOrderF2C(;boundaries)
+    L.(x)
+
+Interpolate a face-value field to a center-valued field from the left, using a 3rd-order reconstruction.
+```math
+L(x)[i+\\tfrac{1}{2}] =  \\left(-2 x[i-1] + 10 x[i] + 4 x[i+1] \\right) / 12
+```
+
+Only the left boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+L(x)[1] = x_0
+```
+"""
+struct LeftBiased3rdOrderF2C{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+LeftBiased3rdOrderF2C(; kwargs...) = LeftBiased3rdOrderF2C(NamedTuple(kwargs))
+
+return_space(::LeftBiased3rdOrderF2C, space::Spaces.FaceFiniteDifferenceSpace) =
+    Spaces.CenterFiniteDifferenceSpace(space)
+return_space(
+    ::LeftBiased3rdOrderF2C,
+    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::LeftBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),)
+Base.@propagate_inbounds stencil_interior(
+    ::LeftBiased3rdOrderF2C,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+) =
+    (
+        -2 * getidx(space, arg, loc, idx - 1 - half, hidx) +
+        10 * getidx(space, arg, loc, idx - half, hidx) +
+        4 * getidx(space, arg, loc, idx + half, hidx)
+    ) / 12
+
+left_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiased3rdOrderF2C,
+    ::AbstractBoundaryCondition,
+    arg,
+) = left_idx(space) + 1
+right_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiased3rdOrderF2C,
+    ::AbstractBoundaryCondition,
+    arg,
+) = right_idx(space)
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::LeftBiased3rdOrderF2C,
+    bc::SetValue,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == left_center_boundary_idx(space)
+    getidx(space, bc.val, loc, nothing, hidx)
+end
+
+
+"""
+    R = RightBiasedF2C(;boundaries)
+    R.(x)
+
+Interpolate a face-valued field to a center-valued field from the right.
+```math
+R(x)[i] = x[i+\\tfrac{1}{2}]
+```
+
+Only the right boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+R(x)[n+\\tfrac{1}{2}] = x_0
+```
+"""
+struct RightBiasedF2C{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+RightBiasedF2C(; kwargs...) = RightBiasedF2C(NamedTuple(kwargs))
+
+return_space(::RightBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) =
+    Spaces.CenterFiniteDifferenceSpace(space)
+return_space(
+    ::RightBiasedF2C,
+    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::RightBiasedF2C, arg) = ((half, half),)
+Base.@propagate_inbounds stencil_interior(
+    ::RightBiasedF2C,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+) = getidx(space, arg, loc, idx + half, hidx)
+
+left_interior_idx(
+    space::AbstractSpace,
+    ::RightBiasedF2C,
+    ::AbstractBoundaryCondition,
+    arg,
+) = left_idx(space)
+right_interior_idx(
+    space::AbstractSpace,
+    ::RightBiasedF2C,
+    ::AbstractBoundaryCondition,
+    arg,
+) = right_idx(space)
+
+right_interior_idx(space::AbstractSpace, ::RightBiasedF2C, ::SetValue, arg) =
+    right_idx(space) - 1
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::RightBiasedF2C,
+    bc::SetValue,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == right_center_boundary_idx(space)
+    getidx(space, bc.val, loc, nothing, hidx)
+end
+
+
+"""
+    R = RightBiased3rdOrderF2C(;boundaries)
+    R.(x)
+
+Interpolate a face-valued field to a center-valued field from the right, using a 3rd-order reconstruction.
+```math
+R(x)[i] = \\left(4 x[i] + 10 x[i+1] -2 x[i+2]  \\right) / 12
+```
+
+Only the right boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+R(x)[n+\\tfrac{1}{2}] = x_0
+```
+"""
+struct RightBiased3rdOrderF2C{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+RightBiased3rdOrderF2C(; kwargs...) = RightBiased3rdOrderF2C(NamedTuple(kwargs))
+
+return_space(
+    ::RightBiased3rdOrderF2C,
+    space::Spaces.FaceFiniteDifferenceSpace,
+) = Spaces.CenterFiniteDifferenceSpace(space)
+return_space(
+    ::RightBiased3rdOrderF2C,
+    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::RightBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),)
+Base.@propagate_inbounds stencil_interior(
+    ::RightBiased3rdOrderF2C,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+) =
+    (
+        4 * getidx(space, arg, loc, idx - half, hidx) +
+        10 * getidx(space, arg, loc, idx + half, hidx) -
+        2 * getidx(space, arg, loc, idx + half + 1, hidx)
+    ) / 12
+
+boundary_width(::RightBiased3rdOrderF2C, ::SetValue) = 1
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::RightBiased3rdOrderF2C,
+    bc::SetValue,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == right_center_boundary_idx(space)
+    getidx(space, bc.val, loc, nothing, hidx)
+end
+
+
+"""
+    U = UpwindBiasedProductC2F(;boundaries)
+    U.(v, x)
+
+Compute the product of the face-valued vector field `v` and a center-valued
+field `x` at cell faces by upwinding `x` according to the direction of `v`.
+
+More precisely, it is computed based on the sign of the 3rd contravariant
+component, and it returns a `Contravariant3Vector`:
+```math
+U(\\boldsymbol{v},x)[i] = \\begin{cases}
+  v^3[i] x[i-\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] > 0 \\\\
+  v^3[i] x[i+\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] < 0
+  \\end{cases}
+```
+where ``\\boldsymbol{e}_3`` is the 3rd covariant basis vector.
+
+Supported boundary conditions are:
+- [`SetValue(x₀)`](@ref): set the value of `x` to be `x₀` in a hypothetical
+  ghost cell on the other side of the boundary. On the left boundary the stencil
+  is
+  ```math
+  U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = \\begin{cases}
+    v^3[\\tfrac{1}{2}] x_0  \\boldsymbol{e}_3 \\textrm{, if }  v^3[\\tfrac{1}{2}] > 0 \\\\
+    v^3[\\tfrac{1}{2}] x[1] \\boldsymbol{e}_3 \\textrm{, if }  v^3[\\tfrac{1}{2}] < 0
+    \\end{cases}
+  ```
+- [`Extrapolate()`](@ref): set the value of `x` to be the same as the closest
+  interior point. On the left boundary, the stencil is
+  ```math
+  U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = U(\\boldsymbol{v},x)[1 + \\tfrac{1}{2}]
+  ```
+"""
+struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator
+    bcs::BCS
+end
+UpwindBiasedProductC2F(; kwargs...) = UpwindBiasedProductC2F(NamedTuple(kwargs))
+
+return_eltype(::UpwindBiasedProductC2F, V, A) =
+    Geometry.Contravariant3Vector{eltype(eltype(V))}
+
+return_space(
+    ::UpwindBiasedProductC2F,
+    velocity_space::Spaces.FaceFiniteDifferenceSpace,
+    arg_space::Spaces.CenterFiniteDifferenceSpace,
+) = velocity_space
+return_space(
+    ::UpwindBiasedProductC2F,
+    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = velocity_space
+
+
+stencil_interior_width(::UpwindBiasedProductC2F, velocity, arg) =
+    ((0, 0), (-half, half))
+
+Base.@propagate_inbounds function stencil_interior(
+    ::UpwindBiasedProductC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    u = stencil_interior(
+        Upwind1stOrderC2F(),
+        loc,
+        space,
+        idx,
+        hidx,
+        Contravariant3Vector(v³),
+        arg,
+    )
+    return Contravariant3Vector(v³ * u)
+end
+
+boundary_width(::UpwindBiasedProductC2F, ::AbstractBoundaryCondition) = 1
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::UpwindBiasedProductC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    u = stencil_left_boundary(
+        Upwind1stOrderC2F(),
+        bc,
+        loc,
+        space,
+        idx,
+        hidx,
+        Contravariant3Vector(v³),
+        arg,
+    )
+    return Contravariant3Vector(v³ * u)
+end
+
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::UpwindBiasedProductC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    u = stencil_right_boundary(
+        Upwind1stOrderC2F(),
+        bc,
+        loc,
+        space,
+        idx,
+        hidx,
+        Contravariant3Vector(v³),
+        arg,
+    )
+    return Contravariant3Vector(v³ * u)
+end
+
+"""
+    U = Upwind3rdOrderBiasedProductC2F(;boundaries)
+    U.(v, x)
+
+Compute the product of a face-valued vector field `v` and a center-valued field
+`x` at cell faces by upwinding `x`, to third-order of accuracy, according to `v`
+```math
+U(v,x)[i] = \\begin{cases}
+  v[i] \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12  \\textrm{, if } v[i] > 0 \\\\
+  v[i] \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}]  \\right) / 12  \\textrm{, if } v[i] < 0
+  \\end{cases}
+```
+This stencil is based on [WickerSkamarock2002](@cite), eq. 4(a).
+
+Supported boundary conditions are:
+- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind scheme to compute `x` on the left boundary,
+  and the first-order upwind scheme to compute `x` on the right boundary.
+and the third-order upwind reconstruction to compute `x` on the right boundary.
+
+!!! note
+    These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a [`DivergenceF2C`](@ref) operator, with a [`SetValue`](@ref) boundary.
+"""
+struct Upwind3rdOrderBiasedProductC2F{BCS} <: AdvectionOperator
+    bcs::BCS
+end
+Upwind3rdOrderBiasedProductC2F(; kwargs...) =
+    Upwind3rdOrderBiasedProductC2F(NamedTuple(kwargs))
+
+return_eltype(::Upwind3rdOrderBiasedProductC2F, V, A) =
+    Geometry.Contravariant3Vector{eltype(eltype(V))}
+
+return_space(
+    ::Upwind3rdOrderBiasedProductC2F,
+    velocity_space::Spaces.FaceFiniteDifferenceSpace,
+    arg_space::Spaces.CenterFiniteDifferenceSpace,
+) = velocity_space
+return_space(
+    ::Upwind3rdOrderBiasedProductC2F,
+    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = velocity_space
+
+
+stencil_interior_width(::Upwind3rdOrderBiasedProductC2F, velocity, arg) =
+    ((0, 0), (-half - 1, half + 1))
+
+Base.@propagate_inbounds function stencil_interior(
+    ::Upwind3rdOrderBiasedProductC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    u = stencil_interior(
+        Upwind3rdOrderC2F(),
+        loc,
+        space,
+        idx,
+        hidx,
+        Contravariant3Vector(v³),
+        arg,
+    )
+    return Contravariant3Vector(v³ * u)
+end
+
+boundary_width(::Upwind3rdOrderBiasedProductC2F, ::AbstractBoundaryCondition) =
+    2
+
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::Upwind3rdOrderBiasedProductC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    u = stencil_left_boundary(
+        Upwind3rdOrderC2F(),
+        bc,
+        loc,
+        space,
+        idx,
+        hidx,
+        Contravariant3Vector(v³),
+        arg,
+    )
+    return Contravariant3Vector(v³ * u)
+end
+
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::Upwind3rdOrderBiasedProductC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    u = stencil_right_boundary(
+        Upwind3rdOrderC2F(),
+        bc,
+        loc,
+        space,
+        idx,
+        hidx,
+        Contravariant3Vector(v³),
+        arg,
+    )
+    return Contravariant3Vector(v³ * u)
+end
diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference/finitedifference.jl
similarity index 76%
rename from src/Operators/finitedifference.jl
rename to src/Operators/finitedifference/finitedifference.jl
index a15cb9cf7e..4fc428a298 100644
--- a/src/Operators/finitedifference.jl
+++ b/src/Operators/finitedifference/finitedifference.jl
@@ -136,18 +136,19 @@ Set the value at the boundary to be the same as the closest interior point.
 struct Extrapolate <: AbstractBoundaryCondition end
 
 """
-    FirstOrderOneSided()
+    OneSided1stOrder()
 
 Use a first-order up/down-wind scheme to compute the value at the boundary.
 """
-struct FirstOrderOneSided <: AbstractBoundaryCondition end
+struct OneSided1stOrder <: AbstractBoundaryCondition end
 
 """
-    ThirdOrderOneSided()
+    OneSided3rdOrder()
 
 Use a third-order up/down-wind scheme to compute the value at the boundary.
 """
-struct ThirdOrderOneSided <: AbstractBoundaryCondition end
+struct OneSided3rdOrder <: AbstractBoundaryCondition end
+
 
 abstract type Location end
 abstract type Boundary <: Location end
@@ -514,1038 +515,225 @@ Base.@propagate_inbounds function stencil_right_boundary(
     a⁻
 end
 
-"""
-    L = LeftBiasedC2F(;boundaries)
-    L.(x)
 
-Interpolate a center-value field to a face-valued field from the left.
-```math
-L(x)[i] = x[i-\\tfrac{1}{2}]
-```
 
-Only the left boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-L(x)[\\tfrac{1}{2}] = x_0
-```
-"""
-struct LeftBiasedC2F{BCS} <: InterpolationOperator
-    bcs::BCS
-end
-LeftBiasedC2F(; kwargs...) = LeftBiasedC2F(NamedTuple(kwargs))
 
-return_space(::LeftBiasedC2F, space::Spaces.CenterFiniteDifferenceSpace) =
-    Spaces.FaceFiniteDifferenceSpace(space)
-return_space(
-    ::LeftBiasedC2F,
-    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
 
-stencil_interior_width(::LeftBiasedC2F, arg) = ((-half, -half),)
-Base.@propagate_inbounds stencil_interior(
-    ::LeftBiasedC2F,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-) = getidx(space, arg, loc, idx - half, hidx)
 
-left_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiasedC2F,
-    ::AbstractBoundaryCondition,
-    arg,
-) = left_idx(space) + 1
-right_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiasedC2F,
-    ::AbstractBoundaryCondition,
-    arg,
-) = right_idx(space)
 
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::LeftBiasedC2F,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-)
-    @assert idx == left_face_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
+abstract type WeightedInterpolationOperator <: InterpolationOperator end
+# TODO: this is not in general correct and the return type
+# should be based on the component operator types (rdiv, rmul) but we don't have a good way
+# of creating ex. one(field_type) for complex fields for inference
+return_eltype(::WeightedInterpolationOperator, weights, arg) = eltype(arg)
 
 """
-    L = LeftBiasedF2C(;boundaries)
-    L.(x)
+    WI = WeightedInterpolateF2C(; boundaries)
+    WI.(w, x)
 
-Interpolate a face-value field to a center-valued field from the left.
+Interpolate a face-valued field `x` to centers, weighted by a face-valued field
+`w`, using the stencil
 ```math
-L(x)[i+\\tfrac{1}{2}] = x[i]
+WI(w, x)[i] = \\frac{
+        w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] +  w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}])
+    }{
+        w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}]
+    }
 ```
 
-Only the left boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-L(x)[1] = x_0
-```
+No boundary conditions are required (or supported)
 """
-struct LeftBiasedF2C{BCS} <: InterpolationOperator
+struct WeightedInterpolateF2C{BCS} <: WeightedInterpolationOperator
     bcs::BCS
 end
-LeftBiasedF2C(; kwargs...) = LeftBiasedF2C(NamedTuple(kwargs))
-
-return_space(::LeftBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) =
-    Spaces.CenterFiniteDifferenceSpace(space)
-return_space(::LeftBiasedF2C, space::Spaces.FaceExtrudedFiniteDifferenceSpace) =
-    Spaces.CenterExtrudedFiniteDifferenceSpace(space)
+WeightedInterpolateF2C(; kwargs...) = WeightedInterpolateF2C(NamedTuple(kwargs))
 
-stencil_interior_width(::LeftBiasedF2C, arg) = ((-half, -half),)
-Base.@propagate_inbounds stencil_interior(
-    ::LeftBiasedF2C,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-) = getidx(space, arg, loc, idx - half, hidx)
-left_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiasedF2C,
-    ::AbstractBoundaryCondition,
-    arg,
-) = left_idx(space)
-right_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiasedF2C,
-    ::AbstractBoundaryCondition,
-    arg,
-) = right_idx(space)
+return_space(
+    ::WeightedInterpolateF2C,
+    weight_space::Spaces.FaceFiniteDifferenceSpace,
+    arg_space::Spaces.FaceFiniteDifferenceSpace,
+) = Spaces.CenterFiniteDifferenceSpace(arg_space)
+return_space(
+    ::WeightedInterpolateF2C,
+    weight_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+    arg_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+) = Spaces.CenterExtrudedFiniteDifferenceSpace(arg_space)
 
-left_interior_idx(space::AbstractSpace, ::LeftBiasedF2C, ::SetValue, arg) =
-    left_idx(space) + 1
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::LeftBiasedF2C,
-    bc::SetValue,
+stencil_interior_width(::WeightedInterpolateF2C, weight, arg) =
+    ((-half, half), (-half, half))
+Base.@propagate_inbounds function stencil_interior(
+    ::WeightedInterpolateF2C,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
 )
-    @assert idx == left_center_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
+    w⁺ = getidx(space, weight, loc, idx + half, hidx)
+    w⁻ = getidx(space, weight, loc, idx - half, hidx)
+    a⁺ = getidx(space, arg, loc, idx + half, hidx)
+    a⁻ = getidx(space, arg, loc, idx - half, hidx)
+    RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻))
 end
 
+boundary_width(::WeightedInterpolateF2C, ::AbstractBoundaryCondition) = 0
+
 """
-    L = LeftBiased3rdOrderC2F(;boundaries)
-    L.(x)
+    WI = WeightedInterpolateC2F(; boundaries)
+    WI.(w, x)
 
-Interpolate a center-value field to a face-valued field from the left, using a 3rd-order reconstruction.
+Interpolate a center-valued field `x` to faces, weighted by a center-valued field
+`w`, using the stencil
 ```math
-L(x)[i] =  \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12
+WI(w, x)[i] = \\frac{
+    w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] +  w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}])
+}{
+    w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}]
+}
 ```
 
-Only the left boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-L(x)[\\tfrac{1}{2}] = x_0
-```
+Supported boundary conditions are:
+
+- [`SetValue(val)`](@ref): set the value at the boundary face to be `val`.
+- [`SetGradient`](@ref): set the value at the boundary such that the gradient is `val`.
+- [`Extrapolate`](@ref): use the closest interior point as the boundary value.
+
+These have the same stencil as in [`InterpolateC2F`](@ref).
 """
-struct LeftBiased3rdOrderC2F{BCS} <: InterpolationOperator
+struct WeightedInterpolateC2F{BCS} <: WeightedInterpolationOperator
     bcs::BCS
 end
-LeftBiased3rdOrderC2F(; kwargs...) = LeftBiased3rdOrderC2F(NamedTuple(kwargs))
+WeightedInterpolateC2F(; kwargs...) = WeightedInterpolateC2F(NamedTuple(kwargs))
 
 return_space(
-    ::LeftBiased3rdOrderC2F,
-    space::Spaces.CenterFiniteDifferenceSpace,
-) = Spaces.FaceFiniteDifferenceSpace(space)
+    ::WeightedInterpolateC2F,
+    weight_space::Spaces.CenterFiniteDifferenceSpace,
+    arg_space::Spaces.CenterFiniteDifferenceSpace,
+) = Spaces.FaceFiniteDifferenceSpace(arg_space)
 return_space(
-    ::LeftBiased3rdOrderC2F,
-    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
+    ::WeightedInterpolateC2F,
+    weight_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = Spaces.FaceExtrudedFiniteDifferenceSpace(arg_space)
 
-stencil_interior_width(::LeftBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),)
-Base.@propagate_inbounds stencil_interior(
-    ::LeftBiased3rdOrderC2F,
+stencil_interior_width(::WeightedInterpolateC2F, weight, arg) =
+    ((-half, half), (-half, half))
+Base.@propagate_inbounds function stencil_interior(
+    ::WeightedInterpolateC2F,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
-) =
-    (
-        -2 * getidx(space, arg, loc, idx - 1 - half, hidx) +
-        10 * getidx(space, arg, loc, idx - half, hidx) +
-        4 * getidx(space, arg, loc, idx + half, hidx)
-    ) / 12
-
-left_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiased3rdOrderC2F,
-    ::AbstractBoundaryCondition,
-    arg,
-) = left_idx(space) + 2
-right_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiased3rdOrderC2F,
-    ::AbstractBoundaryCondition,
-    arg,
-) = right_idx(space) - 1
+)
+    w⁺ = getidx(space, weight, loc, idx + half, hidx)
+    w⁻ = getidx(space, weight, loc, idx - half, hidx)
+    a⁺ = getidx(space, arg, loc, idx + half, hidx)
+    a⁻ = getidx(space, arg, loc, idx - half, hidx)
+    RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻))
+end
 
+boundary_width(::WeightedInterpolateC2F, ::AbstractBoundaryCondition) = 1
 Base.@propagate_inbounds function stencil_left_boundary(
-    ::LeftBiased3rdOrderC2F,
+    ::WeightedInterpolateC2F,
     bc::SetValue,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
 )
     @assert idx == left_face_boundary_idx(space)
     getidx(space, bc.val, loc, nothing, hidx)
 end
-
-"""
-    L = LeftBiased3rdOrderF2C(;boundaries)
-    L.(x)
-
-Interpolate a face-value field to a center-valued field from the left, using a 3rd-order reconstruction.
-```math
-L(x)[i+\\tfrac{1}{2}] =  \\left(-2 x[i-1] + 10 x[i] + 4 x[i+1] \\right) / 12
-```
-
-Only the left boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-L(x)[1] = x_0
-```
-"""
-struct LeftBiased3rdOrderF2C{BCS} <: InterpolationOperator
-    bcs::BCS
-end
-LeftBiased3rdOrderF2C(; kwargs...) = LeftBiased3rdOrderF2C(NamedTuple(kwargs))
-
-return_space(::LeftBiased3rdOrderF2C, space::Spaces.FaceFiniteDifferenceSpace) =
-    Spaces.CenterFiniteDifferenceSpace(space)
-return_space(
-    ::LeftBiased3rdOrderF2C,
-    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)
-
-stencil_interior_width(::LeftBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),)
-Base.@propagate_inbounds stencil_interior(
-    ::LeftBiased3rdOrderF2C,
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::WeightedInterpolateC2F,
+    bc::SetValue,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
-) =
-    (
-        -2 * getidx(space, arg, loc, idx - 1 - half, hidx) +
-        10 * getidx(space, arg, loc, idx - half, hidx) +
-        4 * getidx(space, arg, loc, idx + half, hidx)
-    ) / 12
-
-left_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiased3rdOrderF2C,
-    ::AbstractBoundaryCondition,
-    arg,
-) = left_idx(space) + 1
-right_interior_idx(
-    space::AbstractSpace,
-    ::LeftBiased3rdOrderF2C,
-    ::AbstractBoundaryCondition,
-    arg,
-) = right_idx(space)
+)
+    @assert idx == right_face_boundary_idx(space)
+    getidx(space, bc.val, loc, nothing, hidx)
+end
 
 Base.@propagate_inbounds function stencil_left_boundary(
-    ::LeftBiased3rdOrderF2C,
-    bc::SetValue,
+    ::WeightedInterpolateC2F,
+    bc::SetGradient,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
 )
-    @assert idx == left_center_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
-
-"""
-    R = RightBiasedC2F(;boundaries)
-    R.(x)
-
-Interpolate a center-valued field to a face-valued field from the right.
-```math
-R(x)[i] = x[i+\\tfrac{1}{2}]
-```
-
-Only the right boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-R(x)[n+\\tfrac{1}{2}] = x_0
-```
-"""
-struct RightBiasedC2F{BCS} <: InterpolationOperator
-    bcs::BCS
+    @assert idx == left_face_boundary_idx(space)
+    a⁺ = getidx(space, arg, loc, idx + half, hidx)
+    v₃ = Geometry.covariant3(
+        getidx(space, bc.val, loc, nothing, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    a⁺ ⊟ RecursiveApply.rdiv(v₃, 2)
 end
-RightBiasedC2F(; kwargs...) = RightBiasedC2F(NamedTuple(kwargs))
-
-return_space(::RightBiasedC2F, space::Spaces.CenterFiniteDifferenceSpace) =
-    Spaces.FaceFiniteDifferenceSpace(space)
-return_space(
-    ::RightBiasedC2F,
-    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
-
-stencil_interior_width(::RightBiasedC2F, arg) = ((half, half),)
-Base.@propagate_inbounds stencil_interior(
-    ::RightBiasedC2F,
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::WeightedInterpolateC2F,
+    bc::SetGradient,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
-) = getidx(space, arg, loc, idx + half, hidx)
+)
+    @assert idx == right_face_boundary_idx(space)
+    a⁻ = getidx(space, arg, loc, idx - half, hidx)
+    v₃ = Geometry.covariant3(
+        getidx(space, bc.val, loc, nothing, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    a⁻ ⊞ RecursiveApply.rdiv(v₃, 2)
+end
 
-left_interior_idx(
-    space::AbstractSpace,
-    ::RightBiasedC2F,
-    ::AbstractBoundaryCondition,
-    arg,
-) = left_idx(space)
-right_interior_idx(
-    space::AbstractSpace,
-    ::RightBiasedC2F,
-    ::AbstractBoundaryCondition,
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::WeightedInterpolateC2F,
+    bc::Extrapolate,
+    loc,
+    space,
+    idx,
+    hidx,
+    weight,
     arg,
-) = right_idx(space) - 1
-
+)
+    @assert idx == left_face_boundary_idx(space)
+    a⁺ = getidx(space, arg, loc, idx + half, hidx)
+    a⁺
+end
 Base.@propagate_inbounds function stencil_right_boundary(
-    ::RightBiasedC2F,
-    bc::SetValue,
+    ::WeightedInterpolateC2F,
+    bc::Extrapolate,
     loc,
     space,
     idx,
     hidx,
+    weight,
     arg,
 )
     @assert idx == right_face_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
+    a⁻ = getidx(space, arg, loc, idx - half, hidx)
+    a⁻
 end
 
-"""
-    R = RightBiasedF2C(;boundaries)
-    R.(x)
-
-Interpolate a face-valued field to a center-valued field from the right.
-```math
-R(x)[i] = x[i+\\tfrac{1}{2}]
-```
 
-Only the right boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-R(x)[n+\\tfrac{1}{2}] = x_0
-```
-"""
-struct RightBiasedF2C{BCS} <: InterpolationOperator
-    bcs::BCS
-end
-RightBiasedF2C(; kwargs...) = RightBiasedF2C(NamedTuple(kwargs))
 
-return_space(::RightBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) =
-    Spaces.CenterFiniteDifferenceSpace(space)
-return_space(
-    ::RightBiasedF2C,
-    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)
+abstract type AdvectionOperator <: FiniteDifferenceOperator end
+return_eltype(::AdvectionOperator, velocity, arg) = eltype(arg)
 
-stencil_interior_width(::RightBiasedF2C, arg) = ((half, half),)
-Base.@propagate_inbounds stencil_interior(
-    ::RightBiasedF2C,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-) = getidx(space, arg, loc, idx + half, hidx)
-
-left_interior_idx(
-    space::AbstractSpace,
-    ::RightBiasedF2C,
-    ::AbstractBoundaryCondition,
-    arg,
-) = left_idx(space)
-right_interior_idx(
-    space::AbstractSpace,
-    ::RightBiasedF2C,
-    ::AbstractBoundaryCondition,
-    arg,
-) = right_idx(space)
-
-right_interior_idx(space::AbstractSpace, ::RightBiasedF2C, ::SetValue, arg) =
-    right_idx(space) - 1
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::RightBiasedF2C,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-)
-    @assert idx == right_center_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
-
-
-"""
-    R = RightBiased3rdOrderC2F(;boundaries)
-    R.(x)
-
-Interpolate a center-valued field to a face-valued field from the right, using a 3rd-order reconstruction.
-```math
-R(x)[i] = \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}]  \\right) / 12
-```
-
-Only the right boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-R(x)[n+\\tfrac{1}{2}] = x_0
-```
-"""
-struct RightBiased3rdOrderC2F{BCS} <: InterpolationOperator
-    bcs::BCS
-end
-RightBiased3rdOrderC2F(; kwargs...) = RightBiased3rdOrderC2F(NamedTuple(kwargs))
-
-return_space(
-    ::RightBiased3rdOrderC2F,
-    space::Spaces.CenterFiniteDifferenceSpace,
-) = Spaces.FaceFiniteDifferenceSpace(space)
-return_space(
-    ::RightBiased3rdOrderC2F,
-    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
-
-stencil_interior_width(::RightBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),)
-Base.@propagate_inbounds stencil_interior(
-    ::RightBiased3rdOrderC2F,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-) =
-    (
-        4 * getidx(space, arg, loc, idx - half, hidx) +
-        10 * getidx(space, arg, loc, idx + half, hidx) -
-        2 * getidx(space, arg, loc, idx + half + 1, hidx)
-    ) / 12
-
-boundary_width(::RightBiased3rdOrderC2F, ::SetValue) = 1
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::RightBiased3rdOrderC2F,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-)
-    @assert idx == right_face_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
-
-"""
-    R = RightBiased3rdOrderF2C(;boundaries)
-    R.(x)
-
-Interpolate a face-valued field to a center-valued field from the right, using a 3rd-order reconstruction.
-```math
-R(x)[i] = \\left(4 x[i] + 10 x[i+1] -2 x[i+2]  \\right) / 12
-```
-
-Only the right boundary condition should be set. Currently supported is:
-- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
-```math
-R(x)[n+\\tfrac{1}{2}] = x_0
-```
-"""
-struct RightBiased3rdOrderF2C{BCS} <: InterpolationOperator
-    bcs::BCS
-end
-RightBiased3rdOrderF2C(; kwargs...) = RightBiased3rdOrderF2C(NamedTuple(kwargs))
-
-return_space(
-    ::RightBiased3rdOrderF2C,
-    space::Spaces.FaceFiniteDifferenceSpace,
-) = Spaces.CenterFiniteDifferenceSpace(space)
-return_space(
-    ::RightBiased3rdOrderF2C,
-    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)
-
-stencil_interior_width(::RightBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),)
-Base.@propagate_inbounds stencil_interior(
-    ::RightBiased3rdOrderF2C,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-) =
-    (
-        4 * getidx(space, arg, loc, idx - half, hidx) +
-        10 * getidx(space, arg, loc, idx + half, hidx) -
-        2 * getidx(space, arg, loc, idx + half + 1, hidx)
-    ) / 12
-
-boundary_width(::RightBiased3rdOrderF2C, ::SetValue) = 1
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::RightBiased3rdOrderF2C,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    arg,
-)
-    @assert idx == right_center_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
-
-abstract type WeightedInterpolationOperator <: InterpolationOperator end
-# TODO: this is not in general correct and the return type
-# should be based on the component operator types (rdiv, rmul) but we don't have a good way
-# of creating ex. one(field_type) for complex fields for inference
-return_eltype(::WeightedInterpolationOperator, weights, arg) = eltype(arg)
-
-"""
-    WI = WeightedInterpolateF2C(; boundaries)
-    WI.(w, x)
-
-Interpolate a face-valued field `x` to centers, weighted by a face-valued field
-`w`, using the stencil
-```math
-WI(w, x)[i] = \\frac{
-        w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] +  w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}])
-    }{
-        w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}]
-    }
-```
-
-No boundary conditions are required (or supported)
-"""
-struct WeightedInterpolateF2C{BCS} <: WeightedInterpolationOperator
-    bcs::BCS
-end
-WeightedInterpolateF2C(; kwargs...) = WeightedInterpolateF2C(NamedTuple(kwargs))
-
-return_space(
-    ::WeightedInterpolateF2C,
-    weight_space::Spaces.FaceFiniteDifferenceSpace,
-    arg_space::Spaces.FaceFiniteDifferenceSpace,
-) = Spaces.CenterFiniteDifferenceSpace(arg_space)
-return_space(
-    ::WeightedInterpolateF2C,
-    weight_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-    arg_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-) = Spaces.CenterExtrudedFiniteDifferenceSpace(arg_space)
-
-stencil_interior_width(::WeightedInterpolateF2C, weight, arg) =
-    ((-half, half), (-half, half))
-Base.@propagate_inbounds function stencil_interior(
-    ::WeightedInterpolateF2C,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    w⁺ = getidx(space, weight, loc, idx + half, hidx)
-    w⁻ = getidx(space, weight, loc, idx - half, hidx)
-    a⁺ = getidx(space, arg, loc, idx + half, hidx)
-    a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻))
-end
-
-boundary_width(::WeightedInterpolateF2C, ::AbstractBoundaryCondition) = 0
-
-"""
-    WI = WeightedInterpolateC2F(; boundaries)
-    WI.(w, x)
-
-Interpolate a center-valued field `x` to faces, weighted by a center-valued field
-`w`, using the stencil
-```math
-WI(w, x)[i] = \\frac{
-    w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] +  w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}])
-}{
-    w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}]
-}
-```
-
-Supported boundary conditions are:
-
-- [`SetValue(val)`](@ref): set the value at the boundary face to be `val`.
-- [`SetGradient`](@ref): set the value at the boundary such that the gradient is `val`.
-- [`Extrapolate`](@ref): use the closest interior point as the boundary value.
-
-These have the same stencil as in [`InterpolateC2F`](@ref).
-"""
-struct WeightedInterpolateC2F{BCS} <: WeightedInterpolationOperator
-    bcs::BCS
-end
-WeightedInterpolateC2F(; kwargs...) = WeightedInterpolateC2F(NamedTuple(kwargs))
-
-return_space(
-    ::WeightedInterpolateC2F,
-    weight_space::Spaces.CenterFiniteDifferenceSpace,
-    arg_space::Spaces.CenterFiniteDifferenceSpace,
-) = Spaces.FaceFiniteDifferenceSpace(arg_space)
-return_space(
-    ::WeightedInterpolateC2F,
-    weight_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = Spaces.FaceExtrudedFiniteDifferenceSpace(arg_space)
-
-stencil_interior_width(::WeightedInterpolateC2F, weight, arg) =
-    ((-half, half), (-half, half))
-Base.@propagate_inbounds function stencil_interior(
-    ::WeightedInterpolateC2F,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    w⁺ = getidx(space, weight, loc, idx + half, hidx)
-    w⁻ = getidx(space, weight, loc, idx - half, hidx)
-    a⁺ = getidx(space, arg, loc, idx + half, hidx)
-    a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻))
-end
-
-boundary_width(::WeightedInterpolateC2F, ::AbstractBoundaryCondition) = 1
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::WeightedInterpolateC2F,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    @assert idx == left_face_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::WeightedInterpolateC2F,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    @assert idx == right_face_boundary_idx(space)
-    getidx(space, bc.val, loc, nothing, hidx)
-end
-
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::WeightedInterpolateC2F,
-    bc::SetGradient,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    @assert idx == left_face_boundary_idx(space)
-    a⁺ = getidx(space, arg, loc, idx + half, hidx)
-    v₃ = Geometry.covariant3(
-        getidx(space, bc.val, loc, nothing, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    a⁺ ⊟ RecursiveApply.rdiv(v₃, 2)
-end
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::WeightedInterpolateC2F,
-    bc::SetGradient,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    @assert idx == right_face_boundary_idx(space)
-    a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    v₃ = Geometry.covariant3(
-        getidx(space, bc.val, loc, nothing, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    a⁻ ⊞ RecursiveApply.rdiv(v₃, 2)
-end
-
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::WeightedInterpolateC2F,
-    bc::Extrapolate,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    @assert idx == left_face_boundary_idx(space)
-    a⁺ = getidx(space, arg, loc, idx + half, hidx)
-    a⁺
-end
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::WeightedInterpolateC2F,
-    bc::Extrapolate,
-    loc,
-    space,
-    idx,
-    hidx,
-    weight,
-    arg,
-)
-    @assert idx == right_face_boundary_idx(space)
-    a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    a⁻
-end
-
-
-abstract type AdvectionOperator <: FiniteDifferenceOperator end
-return_eltype(::AdvectionOperator, velocity, arg) = eltype(arg)
-
-"""
-    U = UpwindBiasedProductC2F(;boundaries)
-    U.(v, x)
-
-Compute the product of the face-valued vector field `v` and a center-valued
-field `x` at cell faces by upwinding `x` according to the direction of `v`.
-
-More precisely, it is computed based on the sign of the 3rd contravariant
-component, and it returns a `Contravariant3Vector`:
-```math
-U(\\boldsymbol{v},x)[i] = \\begin{cases}
-  v^3[i] x[i-\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] > 0 \\\\
-  v^3[i] x[i+\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] < 0
-  \\end{cases}
-```
-where ``\\boldsymbol{e}_3`` is the 3rd covariant basis vector.
-
-Supported boundary conditions are:
-- [`SetValue(x₀)`](@ref): set the value of `x` to be `x₀` in a hypothetical
-  ghost cell on the other side of the boundary. On the left boundary the stencil
-  is
-  ```math
-  U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = \\begin{cases}
-    v^3[\\tfrac{1}{2}] x_0  \\boldsymbol{e}_3 \\textrm{, if }  v^3[\\tfrac{1}{2}] > 0 \\\\
-    v^3[\\tfrac{1}{2}] x[1] \\boldsymbol{e}_3 \\textrm{, if }  v^3[\\tfrac{1}{2}] < 0
-    \\end{cases}
-  ```
-- [`Extrapolate()`](@ref): set the value of `x` to be the same as the closest
-  interior point. On the left boundary, the stencil is
-  ```math
-  U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = U(\\boldsymbol{v},x)[1 + \\tfrac{1}{2}]
-  ```
-"""
-struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator
-    bcs::BCS
-end
-UpwindBiasedProductC2F(; kwargs...) = UpwindBiasedProductC2F(NamedTuple(kwargs))
-
-return_eltype(::UpwindBiasedProductC2F, V, A) =
-    Geometry.Contravariant3Vector{eltype(eltype(V))}
-
-return_space(
-    ::UpwindBiasedProductC2F,
-    velocity_space::Spaces.FaceFiniteDifferenceSpace,
-    arg_space::Spaces.CenterFiniteDifferenceSpace,
-) = velocity_space
-return_space(
-    ::UpwindBiasedProductC2F,
-    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = velocity_space
-
-function upwind_biased_product(v, a⁻, a⁺)
-    RecursiveApply.rdiv(
-        ((v ⊞ RecursiveApply.rmap(abs, v)) ⊠ a⁻) ⊞
-        ((v ⊟ RecursiveApply.rmap(abs, v)) ⊠ a⁺),
-        2,
-    )
-end
-
-stencil_interior_width(::UpwindBiasedProductC2F, velocity, arg) =
-    ((0, 0), (-half, half))
-
-Base.@propagate_inbounds function stencil_interior(
-    ::UpwindBiasedProductC2F,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg)
-    a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg)
-    vᶠ = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, a⁻, a⁺))
-end
-
-boundary_width(::UpwindBiasedProductC2F, ::AbstractBoundaryCondition) = 1
-
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::UpwindBiasedProductC2F,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx == left_face_boundary_idx(space)
-    aᴸᴮ = getidx(space, bc.val, loc, nothing, hidx)
-    a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg)
-    vᶠ = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, aᴸᴮ, a⁺))
-end
-
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::UpwindBiasedProductC2F,
-    bc::SetValue,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx == right_face_boundary_idx(space)
-    a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg)
-    aᴿᴮ = getidx(space, bc.val, loc, nothing, hidx)
-    vᶠ = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, a⁻, aᴿᴮ))
-end
-
-Base.@propagate_inbounds function stencil_left_boundary(
-    op::UpwindBiasedProductC2F,
-    ::Extrapolate,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx == left_face_boundary_idx(space)
-    stencil_interior(op, loc, space, idx + 1, hidx, velocity, arg)
-end
-
-Base.@propagate_inbounds function stencil_right_boundary(
-    op::UpwindBiasedProductC2F,
-    ::Extrapolate,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx == right_face_boundary_idx(space)
-    stencil_interior(op, loc, space, idx - 1, hidx, velocity, arg)
-end
-
-"""
-    U = Upwind3rdOrderBiasedProductC2F(;boundaries)
-    U.(v, x)
-
-Compute the product of a face-valued vector field `v` and a center-valued field
-`x` at cell faces by upwinding `x`, to third-order of accuracy, according to `v`
-```math
-U(v,x)[i] = \\begin{cases}
-  v[i] \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12  \\textrm{, if } v[i] > 0 \\\\
-  v[i] \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}]  \\right) / 12  \\textrm{, if } v[i] < 0
-  \\end{cases}
-```
-This stencil is based on [WickerSkamarock2002](@cite), eq. 4(a).
-
-Supported boundary conditions are:
-- [`FirstOrderOneSided(x₀)`](@ref): uses the first-order downwind scheme to compute `x` on the left boundary,
-  and the first-order upwind scheme to compute `x` on the right boundary.
-- [`ThirdOrderOneSided(x₀)`](@ref): uses the third-order downwind reconstruction to compute `x` on the left boundary,
-and the third-order upwind reconstruction to compute `x` on the right boundary.
-
-!!! note
-    These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a [`DivergenceF2C`](@ref) operator, with a [`SetValue`](@ref) boundary.
-"""
-struct Upwind3rdOrderBiasedProductC2F{BCS} <: AdvectionOperator
-    bcs::BCS
-end
-Upwind3rdOrderBiasedProductC2F(; kwargs...) =
-    Upwind3rdOrderBiasedProductC2F(NamedTuple(kwargs))
-
-return_eltype(::Upwind3rdOrderBiasedProductC2F, V, A) =
-    Geometry.Contravariant3Vector{eltype(eltype(V))}
-
-return_space(
-    ::Upwind3rdOrderBiasedProductC2F,
-    velocity_space::Spaces.FaceFiniteDifferenceSpace,
-    arg_space::Spaces.CenterFiniteDifferenceSpace,
-) = velocity_space
-return_space(
-    ::Upwind3rdOrderBiasedProductC2F,
-    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
-    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
-) = velocity_space
-
-function upwind_3rdorder_biased_product(v, a⁻, a⁻⁻, a⁺, a⁺⁺)
-    RecursiveApply.rdiv(
-        (v ⊠ (7 ⊠ (a⁺ + a⁻) ⊟ (a⁺⁺ + a⁻⁻))) ⊟
-        (RecursiveApply.rmap(abs, v) ⊠ (3 ⊠ (a⁺ - a⁻) ⊟ (a⁺⁺ - a⁻⁻))),
-        12,
-    )
-end
-
-stencil_interior_width(::Upwind3rdOrderBiasedProductC2F, velocity, arg) =
-    ((0, 0), (-half - 1, half + 1))
-
-Base.@propagate_inbounds function stencil_interior(
-    ::Upwind3rdOrderBiasedProductC2F,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    a⁻⁻ = getidx(space, arg, loc, idx - half - 1, hidx)
-    a⁺ = getidx(space, arg, loc, idx + half, hidx)
-    a⁺⁺ = getidx(space, arg, loc, idx + half + 1, hidx)
-    vᶠ = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    return Geometry.Contravariant3Vector(
-        upwind_3rdorder_biased_product(vᶠ, a⁻, a⁻⁻, a⁺, a⁺⁺),
-    )
-end
-
-boundary_width(::Upwind3rdOrderBiasedProductC2F, ::AbstractBoundaryCondition) =
-    2
-
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::Upwind3rdOrderBiasedProductC2F,
-    bc::FirstOrderOneSided,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx <= left_face_boundary_idx(space) + 1
-    v = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg)
-    a⁺ = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
-    return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺))
-end
-
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::Upwind3rdOrderBiasedProductC2F,
-    bc::FirstOrderOneSided,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx >= right_face_boundary_idx(space) - 1
-    v = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    a⁻ = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
-    a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg)
-    return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺))
-
-end
-
-Base.@propagate_inbounds function stencil_left_boundary(
-    ::Upwind3rdOrderBiasedProductC2F,
-    bc::ThirdOrderOneSided,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx <= left_face_boundary_idx(space) + 1
-
-    vᶠ = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    a = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
-
-    return Geometry.Contravariant3Vector(vᶠ * a)
-end
-
-Base.@propagate_inbounds function stencil_right_boundary(
-    ::Upwind3rdOrderBiasedProductC2F,
-    bc::ThirdOrderOneSided,
-    loc,
-    space,
-    idx,
-    hidx,
-    velocity,
-    arg,
-)
-    @assert idx <= right_face_boundary_idx(space) - 1
-
-    vᶠ = Geometry.contravariant3(
-        getidx(space, velocity, loc, idx, hidx),
-        Geometry.LocalGeometry(space, idx, hidx),
-    )
-    a = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
-
-    return Geometry.Contravariant3Vector(vᶠ * a)
-end
 
 """
     U = FCTBorisBook(;boundaries)
@@ -1564,7 +752,7 @@ where ``s[i] = +1`` if  `` v[i] \\geq 0`` and ``s[i] = -1`` if  `` v[i] \\leq 0`
 This formulation is based on [BorisBook1973](@cite), as reported in [durran2010](@cite) section 5.4.1.
 
 Supported boundary conditions are:
-- [`FirstOrderOneSided(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary.
+- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary.
 
 !!! note
     Similar to the [`Upwind3rdOrderBiasedProductC2F`](@ref) operator, these boundary conditions do not define the value at the actual boundary faces,
@@ -1645,7 +833,7 @@ boundary_width(::FCTBorisBook, ::AbstractBoundaryCondition) = 2
 
 Base.@propagate_inbounds function stencil_left_boundary(
     ::FCTBorisBook,
-    bc::FirstOrderOneSided,
+    bc::OneSided1stOrder,
     loc,
     space,
     idx,
@@ -1664,7 +852,7 @@ end
 
 Base.@propagate_inbounds function stencil_right_boundary(
     ::FCTBorisBook,
-    bc::FirstOrderOneSided,
+    bc::OneSided1stOrder,
     loc,
     space,
     idx,
@@ -1701,7 +889,7 @@ This stencil is based on [zalesak1979fully](@cite), as reported in [durran2010](
 the corrected antidiffusive flux.
 
 Supported boundary conditions are:
-- [`FirstOrderOneSided(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary.
+- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary.
 
 !!! note
     Similar to the [`Upwind3rdOrderBiasedProductC2F`](@ref) operator, these boundary conditions do not define
@@ -1838,7 +1026,7 @@ boundary_width(::FCTZalesak, ::AbstractBoundaryCondition) = 2
 
 Base.@propagate_inbounds function stencil_left_boundary(
     ::FCTZalesak,
-    bc::FirstOrderOneSided,
+    bc::OneSided1stOrder,
     loc,
     space,
     idx,
@@ -1854,7 +1042,7 @@ end
 
 Base.@propagate_inbounds function stencil_right_boundary(
     ::FCTZalesak,
-    bc::FirstOrderOneSided,
+    bc::OneSided1stOrder,
     loc,
     space,
     idx,
diff --git a/src/Operators/finitedifference/upwinding.jl b/src/Operators/finitedifference/upwinding.jl
new file mode 100644
index 0000000000..8ee1f57535
--- /dev/null
+++ b/src/Operators/finitedifference/upwinding.jl
@@ -0,0 +1,543 @@
+"""
+    L = LeftBiased1stOrderC2F(;boundaries)
+    L.(x)
+
+Interpolate a center-value field to a face-valued field from the left.
+```math
+L(x)[i] = x[i-\\tfrac{1}{2}]
+```
+
+Only the left boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+L(x)[\\tfrac{1}{2}] = x_0
+```
+"""
+struct LeftBiased1stOrderC2F{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+LeftBiased1stOrderC2F(; kwargs...) = LeftBiased1stOrderC2F(NamedTuple(kwargs))
+
+return_space(
+    ::LeftBiased1stOrderC2F,
+    space::Spaces.CenterFiniteDifferenceSpace,
+) = Spaces.FaceFiniteDifferenceSpace(space)
+return_space(
+    ::LeftBiased1stOrderC2F,
+    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::LeftBiased1stOrderC2F, arg) = ((-half, -half),)
+Base.@propagate_inbounds stencil_interior(
+    ::LeftBiased1stOrderC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+) = getidx(space, arg, loc, idx - half, hidx)
+
+left_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiased1stOrderC2F,
+    ::AbstractBoundaryCondition,
+    arg,
+) = left_idx(space) + 1
+right_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiased1stOrderC2F,
+    ::AbstractBoundaryCondition,
+    arg,
+) = right_idx(space)
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::LeftBiased1stOrderC2F,
+    bc::SetValue,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == left_face_boundary_idx(space)
+    getidx(space, bc.val, loc, nothing, hidx)
+end
+
+
+
+"""
+    L = LeftBiased3rdOrderC2F(;boundaries)
+    L.(x)
+
+Interpolate a center-value field to a face-valued field from the left, using a 3rd-order reconstruction.
+```math
+L(x)[i] =  \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12
+```
+
+Only the left boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+L(x)[\\tfrac{1}{2}] = x_0
+```
+"""
+struct LeftBiased3rdOrderC2F{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+LeftBiased3rdOrderC2F(; kwargs...) = LeftBiased3rdOrderC2F(NamedTuple(kwargs))
+
+return_space(
+    ::LeftBiased3rdOrderC2F,
+    space::Spaces.CenterFiniteDifferenceSpace,
+) = Spaces.FaceFiniteDifferenceSpace(space)
+return_space(
+    ::LeftBiased3rdOrderC2F,
+    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::LeftBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),)
+Base.@propagate_inbounds function stencil_interior(
+    ::LeftBiased3rdOrderC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    FT = Spaces.undertype(space)
+
+    FT(4 / 12) ⊠ getidx(space, arg, loc, idx + half, hidx) ⊞
+    FT(10 / 12) ⊠ getidx(space, arg, loc, idx - half, hidx) ⊟
+    FT(2 / 12) ⊠ getidx(space, arg, loc, idx - half - 1, hidx)
+end
+
+left_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiased3rdOrderC2F,
+    ::AbstractBoundaryCondition,
+    arg,
+) = left_idx(space) + 2
+right_interior_idx(
+    space::AbstractSpace,
+    ::LeftBiased3rdOrderC2F,
+    ::AbstractBoundaryCondition,
+    arg,
+) = right_idx(space) - 1
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::LeftBiased3rdOrderC2F,
+    bc::OneSided1stOrder,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == left_face_boundary_idx(space) + 1
+    stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+end
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::LeftBiased3rdOrderC2F,
+    bc::OneSided3rdOrder,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == left_face_boundary_idx(space) + 1
+    stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
+end
+
+
+"""
+    R = RightBiased1stOrderC2F(;boundaries)
+    R.(x)
+
+Interpolate a center-valued field to a face-valued field from the right.
+```math
+R(x)[i] = x[i+\\tfrac{1}{2}]
+```
+
+Only the right boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+R(x)[n+\\tfrac{1}{2}] = x_0
+```
+"""
+struct RightBiased1stOrderC2F{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+RightBiased1stOrderC2F(; kwargs...) = RightBiased1stOrderC2F(NamedTuple(kwargs))
+
+return_space(
+    ::RightBiased1stOrderC2F,
+    space::Spaces.CenterFiniteDifferenceSpace,
+) = Spaces.FaceFiniteDifferenceSpace(space)
+return_space(
+    ::RightBiased1stOrderC2F,
+    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::RightBiased1stOrderC2F, arg) = ((half, half),)
+Base.@propagate_inbounds stencil_interior(
+    ::RightBiased1stOrderC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+) = getidx(space, arg, loc, idx + half, hidx)
+
+left_interior_idx(
+    space::AbstractSpace,
+    ::RightBiased1stOrderC2F,
+    ::AbstractBoundaryCondition,
+    arg,
+) = left_idx(space)
+right_interior_idx(
+    space::AbstractSpace,
+    ::RightBiased1stOrderC2F,
+    ::AbstractBoundaryCondition,
+    arg,
+) = right_idx(space) - 1
+
+
+"""
+    R = RightBiased3rdOrderC2F(;boundaries)
+    R.(x)
+
+Interpolate a center-valued field to a face-valued field from the right, using a 3rd-order reconstruction.
+```math
+R(x)[i] = \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}]  \\right) / 12
+```
+
+Only the right boundary condition should be set. Currently supported is:
+- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
+```math
+R(x)[n+\\tfrac{1}{2}] = x_0
+```
+"""
+struct RightBiased3rdOrderC2F{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+RightBiased3rdOrderC2F(; kwargs...) = RightBiased3rdOrderC2F(NamedTuple(kwargs))
+
+return_space(
+    ::RightBiased3rdOrderC2F,
+    space::Spaces.CenterFiniteDifferenceSpace,
+) = Spaces.FaceFiniteDifferenceSpace(space)
+return_space(
+    ::RightBiased3rdOrderC2F,
+    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)
+
+stencil_interior_width(::RightBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),)
+Base.@propagate_inbounds function stencil_interior(
+    ::RightBiased3rdOrderC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    FT = Spaces.undertype(space)
+
+    FT(4 / 12) ⊠ getidx(space, arg, loc, idx - half, hidx) ⊞
+    FT(10 / 12) ⊠ getidx(space, arg, loc, idx + half, hidx) ⊟
+    FT(2 / 12) ⊠ getidx(space, arg, loc, idx + half + 1, hidx)
+end
+
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::RightBiased3rdOrderC2F,
+    bc::OneSided1stOrder,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == right_face_boundary_idx(space) - 1
+    stencil_interior(RightBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+end
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::RightBiased3rdOrderC2F,
+    bc::OneSided3rdOrder,
+    loc,
+    space,
+    idx,
+    hidx,
+    arg,
+)
+    @assert idx == right_face_boundary_idx(space) - 1
+    stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
+end
+
+"""
+    U = Upwind1stOrderC2F(;boundaries)
+    U.(v, x)
+
+Interpolate the center-valued field `x` to cell faces by 1st-order upwinding `x`
+according to the direction of `v`.
+
+More precisely, it is computed based on the sign of the 3rd contravariant
+component, and it returns a `Contravariant3Vector`:
+```math
+U(\\boldsymbol{v},x)[i] = \\begin{cases}
+  x[i-\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] > 0 \\\\
+  x[i+\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] < 0
+  \\end{cases}
+```
+where ``\\boldsymbol{e}_3`` is the 3rd covariant basis vector.
+
+Supported boundary conditions are:
+- [`SetValue(x₀)`](@ref): set the value of `x` to be `x₀` in a hypothetical
+  ghost cell on the other side of the boundary. On the left boundary the stencil
+  is
+  ```math
+  U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = \\begin{cases}
+    x_0  \\boldsymbol{e}_3 \\textrm{, if }  v^3[\\tfrac{1}{2}] > 0 \\\\
+    x[1] \\boldsymbol{e}_3 \\textrm{, if }  v^3[\\tfrac{1}{2}] < 0
+    \\end{cases}
+  ```
+- [`Extrapolate()`](@ref): set the value of `x` to be the same as the closest
+  interior point. On the left boundary, the stencil is
+  ```math
+  U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = U(\\boldsymbol{v},x)[1 + \\tfrac{1}{2}]
+  ```
+"""
+struct Upwind1stOrderC2F{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+Upwind1stOrderC2F(; kwargs...) = Upwind1stOrderC2F(NamedTuple(kwargs))
+
+return_eltype(::Upwind1stOrderC2F, v, arg) = eltype(arg)
+
+return_space(
+    ::Upwind1stOrderC2F,
+    velocity_space::Spaces.FaceFiniteDifferenceSpace,
+    arg_space::Spaces.CenterFiniteDifferenceSpace,
+) = velocity_space
+return_space(
+    ::Upwind1stOrderC2F,
+    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = velocity_space
+
+stencil_interior_width(::Upwind1stOrderC2F, velocity, arg) =
+    ((0, 0), (-half, half))
+
+
+Base.@propagate_inbounds function stencil_interior(
+    ::Upwind1stOrderC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    if signbit(v³)
+        # v³ < 0        
+        stencil_interior(RightBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+    else
+        # v³ > 0
+        stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+    end
+end
+
+boundary_width(::Upwind1stOrderC2F, ::AbstractBoundaryCondition) = 1
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::Upwind1stOrderC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    @assert idx == left_face_boundary_idx(space)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    if signbit(v³)
+        # v³ < 0        
+        stencil_interior(RightBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+    else
+        # v³ > 0
+        stencil_left_boundary(
+            LeftBiased1stOrderC2F(),
+            bc,
+            loc,
+            space,
+            idx,
+            hidx,
+            arg,
+        )
+    end
+end
+
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::Upwind1stOrderC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    @assert idx == right_face_boundary_idx(space)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    if signbit(v³)
+        # v³ < 0        
+        stencil_right_boundary(
+            RightBiased1stOrderC2F(),
+            bc,
+            loc,
+            space,
+            idx,
+            hidx,
+            arg,
+        )
+    else
+        # v³ > 0
+        stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+    end
+end
+
+
+"""
+    U = Upwind3rdOrderC2F(;boundaries)
+    U.(v, x)
+
+Interpolate the center-valued field `x` to cell faces by 3rd-order upwinding `x`
+according to the direction of `v`.
+    
+Compute the product of a face-valued vector field `v` and a center-valued field
+`x` at cell faces by upwinding `x`, to third-order of accuracy, according to `v`
+```math
+U(v,x)[i] = \\begin{cases}
+  \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12  \\textrm{, if } v[i] > 0 \\\\
+  \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}]  \\right) / 12  \\textrm{, if } v[i] < 0
+  \\end{cases}
+```
+This stencil is based on [WickerSkamarock2002](@cite), eq. 4(a).
+
+Supported boundary conditions are:
+- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind scheme to compute `x` on the left boundary,
+  and the first-order upwind scheme to compute `x` on the right boundary.
+- [`OneSided3rdOrder(x₀)`](@ref): uses the third-order downwind reconstruction to compute `x` on the left boundary,
+and the third-order upwind reconstruction to compute `x` on the right boundary.
+
+!!! note
+    These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a [`DivergenceF2C`](@ref) operator, with a [`SetValue`](@ref) boundary.
+"""
+struct Upwind3rdOrderC2F{BCS} <: InterpolationOperator
+    bcs::BCS
+end
+Upwind3rdOrderC2F(; kwargs...) = Upwind3rdOrderC2F(NamedTuple(kwargs))
+
+return_eltype(::Upwind3rdOrderC2F, V, arg) = eltype(arg)
+
+return_space(
+    ::Upwind3rdOrderC2F,
+    velocity_space::Spaces.FaceFiniteDifferenceSpace,
+    arg_space::Spaces.CenterFiniteDifferenceSpace,
+) = velocity_space
+return_space(
+    ::Upwind3rdOrderC2F,
+    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+    arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+) = velocity_space
+
+stencil_interior_width(::Upwind3rdOrderC2F, velocity, arg) =
+    ((0, 0), (-half - 1, half + 1))
+
+Base.@propagate_inbounds function stencil_interior(
+    ::Upwind3rdOrderC2F,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    if signbit(v³)
+        stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
+    else
+        stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
+    end
+end
+
+boundary_width(::Upwind3rdOrderC2F, ::AbstractBoundaryCondition) = 2
+
+Base.@propagate_inbounds function stencil_left_boundary(
+    ::Upwind3rdOrderC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    if signbit(v³)
+        stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg)
+    else
+        stencil_left_boundary(
+            LeftBiased1stOrderC2F(),
+            bc,
+            loc,
+            space,
+            idx,
+            hidx,
+            arg,
+        )
+    end
+end
+Base.@propagate_inbounds function stencil_right_boundary(
+    ::Upwind3rdOrderC2F,
+    bc::AbstractBoundaryCondition,
+    loc,
+    space,
+    idx,
+    hidx,
+    velocity,
+    arg,
+)
+    v³ = Geometry.contravariant3(
+        getidx(space, velocity, loc, idx, hidx),
+        Geometry.LocalGeometry(space, idx, hidx),
+    )
+    if signbit(v³)
+        stencil_right_boundary(
+            RightBiased3rdOrderC2F(),
+            bc,
+            loc,
+            space,
+            idx,
+            hidx,
+            arg,
+        )
+    else
+        stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg)
+    end
+end
diff --git a/test/Operators/finitedifference/column.jl b/test/Operators/finitedifference/column.jl
index 98d0ef6a24..c99b6f0cd6 100644
--- a/test/Operators/finitedifference/column.jl
+++ b/test/Operators/finitedifference/column.jl
@@ -731,7 +731,7 @@ end
     @test conv_adv_wc[3] ≈ 2 atol = 0.1
 end
 
-@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with FirstOrderOneSided + DivergenceF2C SetValue BCs, constant w" begin
+@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided1stOrder + DivergenceF2C SetValue BCs, constant w" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -761,8 +761,8 @@ end
             s = sin.(centers)
 
             third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
-                bottom = Operators.FirstOrderOneSided(),
-                top = Operators.FirstOrderOneSided(),
+                bottom = Operators.OneSided1stOrder(),
+                top = Operators.OneSided1stOrder(),
             )
 
             divf2c = Operators.DivergenceF2C(
@@ -792,7 +792,7 @@ end
     end
 end
 
-@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with ThirdOrderOneSided + DivergenceF2C SetValue BCs, varying sign w" begin
+@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder + DivergenceF2C SetValue BCs, varying sign w" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -821,8 +821,8 @@ end
             c = sin.(centers)
 
             third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
-                bottom = Operators.ThirdOrderOneSided(),
-                top = Operators.ThirdOrderOneSided(),
+                bottom = Operators.OneSided3rdOrder(),
+                top = Operators.OneSided3rdOrder(),
             )
 
             divf2c = Operators.DivergenceF2C(
@@ -903,7 +903,7 @@ end
     @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2]
 end
 
-@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with FirstOrderOneSided BCs" begin
+@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided1stOrder BCs" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -938,8 +938,8 @@ end
                 top = Operators.Extrapolate(),
             )
             third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
-                bottom = Operators.FirstOrderOneSided(),
-                top = Operators.FirstOrderOneSided(),
+                bottom = Operators.OneSided1stOrder(),
+                top = Operators.OneSided1stOrder(),
             )
 
             divf2c = Operators.DivergenceF2C(
@@ -973,7 +973,7 @@ end
     end
 end
 
-@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with ThirdOrderOneSided BCs" begin
+@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder BCs" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -1006,8 +1006,8 @@ end
                 top = Operators.Extrapolate(),
             )
             third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
-                bottom = Operators.ThirdOrderOneSided(),
-                top = Operators.ThirdOrderOneSided(),
+                bottom = Operators.OneSided3rdOrder(),
+                top = Operators.OneSided3rdOrder(),
             )
 
             divf2c = Operators.DivergenceF2C(
diff --git a/test/Operators/finitedifference/column_benchmark_utils.jl b/test/Operators/finitedifference/column_benchmark_utils.jl
index b9467127f3..bf78459d37 100644
--- a/test/Operators/finitedifference/column_benchmark_utils.jl
+++ b/test/Operators/finitedifference/column_benchmark_utils.jl
@@ -187,18 +187,18 @@ function set_vec_value_bcs(c)
 end
 
 function set_upwind_biased_bcs(c)
-    return (;bottom = Operators.FirstOrderOneSided(),
-             top = Operators.FirstOrderOneSided())
+    return (;bottom = Operators.OneSided1stOrder(),
+             top = Operators.OneSided1stOrder())
 end
 
 function set_upwind_biased_bcs(c)
-    return (;bottom = Operators.ThirdOrderOneSided(),
-             top = Operators.ThirdOrderOneSided())
+    return (;bottom = Operators.OneSided3rdOrder(),
+             top = Operators.OneSided3rdOrder())
 end
 
 function set_upwind_biased_3_bcs(c)
-    return (;bottom = Operators.ThirdOrderOneSided(),
-             top = Operators.ThirdOrderOneSided())
+    return (;bottom = Operators.OneSided3rdOrder(),
+             top = Operators.OneSided3rdOrder())
 end
 
 function set_top_value_bc(c)
@@ -438,7 +438,7 @@ function test_results(t_ave)
     @test t_ave[(:no_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 2.926*μs*buffer
     @test t_ave[(:no_h_space, op_UpwindBiasedProductC2F!, :SetValue, :SetValue)] < 697.341*ns*buffer
     @test t_ave[(:no_h_space, op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate)] < 659.267*ns*buffer
-    @test t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 4.483*μs*buffer
+    @test t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 4.483*μs*buffer
     @test t_ave[(:no_h_space, op_divgrad_CC!, :SetValue, :SetValue, :none)] < 1.607*μs*buffer
     @test t_ave[(:no_h_space, op_divgrad_FF!, :none, :SetDivergence, :SetDivergence)] < 1.529*μs*buffer
     @test t_ave[(:no_h_space, op_div_interp_CC!, :SetValue, :SetValue, :none)] < 1.510*μs*buffer
@@ -468,7 +468,7 @@ function test_results(t_ave)
     @test t_ave[(:has_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 3.159*μs*buffer
     @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :SetValue, :SetValue)] < 5.197*μs*buffer
     @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate)] < 5.304*μs*buffer
-    @test t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 14.304*μs*buffer
+    @test t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 14.304*μs*buffer
     @test t_ave[(:has_h_space, op_divgrad_CC!, :SetValue, :SetValue, :none)] < 8.593*μs*buffer
     @test t_ave[(:has_h_space, op_divgrad_FF!, :none, :SetDivergence, :SetDivergence)] < 8.597*μs*buffer
     @test t_ave[(:has_h_space, op_div_interp_CC!, :SetValue, :SetValue, :none)] < 1.735*μs*buffer
@@ -479,7 +479,7 @@ function test_results(t_ave)
     # Broken tests
     @test_broken t_ave[(:no_h_space, op_CurlC2F!, :SetCurl, :SetCurl)] < 500
     @test_broken t_ave[(:no_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 500
-    @test_broken t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 500
+    @test_broken t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 500
     @test_broken t_ave[(:no_h_space, op_divgrad_uₕ!, :none, :SetValue, :Extrapolate)] < 500
     @test_broken t_ave[(:no_h_space, op_divgrad_uₕ!, :none, :SetValue, :SetValue)] < 500 # different with/without h_space
     @test_broken t_ave[(:has_h_space, op_DivergenceF2C!, :none)] < 500
@@ -489,7 +489,7 @@ function test_results(t_ave)
     @test_broken t_ave[(:has_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 500
     @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :SetValue, :SetValue)] < 800
     @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate)] < 800
-    @test_broken t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 500
+    @test_broken t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 500
     @test_broken t_ave[(:has_h_space, op_divgrad_CC!, :SetValue, :SetValue, :none)] < 500
     @test_broken t_ave[(:has_h_space, op_divgrad_FF!, :none, :SetDivergence, :SetDivergence)] < 500
     @test t_ave[(:has_h_space, op_div_interp_CC!, :SetValue, :SetValue, :none)] < 800