Use bcs with operators in upwind_advection_velocity

src/Operators.jl

@@ -56,6 +56,21 @@ c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] -
 # top of the domain.
 c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
 
+function f∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
+    if is_surface_face(grid, k)
+        return f∇_onesided(f_dual, grid, BottomBCTag(), bottom)
+    elseif is_toa_face(grid, k)
+        return f∇_onesided(f_dual, grid, TopBCTag(), top)
+    else
+        return f∇_onesided(f_dual, grid, InteriorTag())
+    end
+end
+f∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi
+f∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2)
+f∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value
+f∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info
+f∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
+
 # Used when traversing cell faces
 
 interpc2f(f, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) =
@@ -116,14 +131,16 @@ function upwind_advection_area(ρ0_half::Vector{Float64}, a_up::Vector{Float64},
     return -∇m / ρ0_half[k + 1]
 end
 
-
 function upwind_advection_velocity(ρ0::Vector{Float64}, a_up::Vector{Float64}, w_up::Vector{Float64}, grid, k)
-    a_k = interp2pt(a_up[k], a_up[k + 1])
-    a_km = interp2pt(a_up[k - 1], a_up[k])
-    adv = (ρ0[k] * a_k * w_up[k] * w_up[k] * grid.dzi - ρ0[k - 1] * a_km * w_up[k - 1] * w_up[k - 1] * grid.dzi)
-    return adv
+    a_dual = fdaul_onesided(a_up, grid, k)
+    ρ_0_dual = dual_faces(ρ0, grid, k)
+    w_up_dual = dual_faces(w_up, grid, k)
+    adv_dual = a_dual .* ρ_0_dual .* w_up_dual .* w_up_dual
+    ∇ρaw = f∇_onesided(adv_dual, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
+    return ∇ρaw
 end
 
+
 function upwind_advection_scalar(
     ρ0_half::Vector{Float64},
     a_up::Vector{Float64},