Mask ReducedFields accounting the immersed column

[simone-silvestri]

From a MWE similar to the one in #3899 (adding some more reduced directions):

using Oceananigans
using GLMakie

grid = RectilinearGrid(size = (10, 10, 10), extent = (10, 10, 10))
bottom(x, y) = - rand() * 5 # Between -5 and 0
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom))

f1 = Field{Center, Center, Nothing}(grid)
f2 = Field{Center, Nothing, Center}(grid)
f3 = Field{Nothing, Center, Center}(grid)
set!(f1, (x, y) -> rand())
set!(f2, (x, y) -> rand())
set!(f3, (x, y) -> rand())

fig = Figure(); ax1 = Axis(fig[1, 1]); ax2 = Axis(fig[1, 2]); ax3 = Axis(fig[1, 3])
heatmap!(ax1, f1)
heatmap!(ax2, f2)
heatmap!(ax3, f3)

Should we add some tests?

closes #3899