fill_halo_regions! for velocities in implicit free surface correctly

[glwagner]

cc @tomchor @simone-silvestri

[tomchor]

Thanks, @glwagner. I can confirm that if I run the following MWE things work as expected in this branch:

using Oceananigans
grid = RectilinearGrid(topology = (Bounded, Flat, Bounded), size = (4, 4), extent = (1, 1))

u₀ = 1
@inline u_func(z, t) = u₀

u_bcs = FieldBoundaryConditions(east = OpenBoundaryCondition(u_func), west = OpenBoundaryCondition(u_func))

model = HydrostaticFreeSurfaceModel(; grid, boundary_conditions = (; u = u_bcs,))
set!(model, u = u₀)

simulation = Simulation(model, Δt = 0.2*(minimum_xspacing(grid) / u₀), stop_time = 100)
run!(simulation)

And the end result is also what we would expect: a barotropic flow with u velocity u₀ = 1 everywhere:

julia> model.velocities.u
5×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 4×1×4 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 3×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Open, east: Open, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 11×1×10 OffsetArray(::Array{Float64, 3}, -2:8, 1:1, -2:7) with eltype Float64 with indices -2:8×1:1×-2:7
    └── max=1.0, min=1.0, mean=1.0

Although I understand that this does not mean results will be correct for non-trivial open boundary conditions.

[tomchor]

@glwagner merge?

[navidcy]

@tomchor if tests pass feel free to merge ;)

[navidcy]

@tomchor is this related to #3628? we should link the issue..