`interpolate!` fails silently on `Flat` topologies with explicitly specified coordinate in `Flat` directions

This is perhaps an obscure problem as it happens on `Flat` topologies where those flat topologies have specified coordinate values. Here's an MWE: ```julia using Oceananigans using Oceananigans.Fields: interpolate!, interpolate arch = CPU() z = (-1, 0) # Source: 4 cells, default x = y = 0 source_grid = RectilinearGrid(arch; size=4, z=z, topology=(Flat, Flat, Bounded)) source = CenterField(source_grid) interior(source, 1, 1, :) .= [1.0, 2.0, 3.0, 4.0] # case A: target with omitted x and y on an XYFlatGrid target_A_grid = RectilinearGrid(arch; size=40, z=z, topology=(Flat, Flat, Bounded)) target_A = CenterField(target_A_grid) interpolate!(target_A, source) # case B: target with explicitly specified x and y on an XYFlatGrid target_B_grid = RectilinearGrid(arch; size=40, x=-144.9, y=50.1, z=z, topology=(Flat, Flat, Bounded)) target_B = CenterField(target_B_grid) interpolate!(target_B, source) # case C: scalar interpolate per cell (1-tuple path; correct dispatch) target_C = CenterField(target_B_grid) z_targets = collect(znodes(target_B_grid, Center())) target_C_vals = [interpolate((z,), source) for z in z_targets] interior(target_C, 1, 1, :) .= target_C_vals @info "case A (target x=y=0) interior[1:5] = ", vec(interior(target_A))[1:5] @info "case B (target x=-144.9,y=50.1) interior[1:5] = ", vec(interior(target_B))[1:5] @info "case C (scalar loop, same target as B) [1:5] = ", vec(interior(target_C))[1:5] ``` which gives ```julia [ Info: ("case A (target x=y=0) interior[1:5] = ", [0.5499999999999998, 0.6499999999999999, 0.75, 0.8500000000000001, 0.9500000000000002]) [ Info: ("case B (target x=-144.9,y=50.1) interior[1:5] = ", [6.9145708468793e-311, 6.9145708468793e-311, 6.9145708468793e-311, 6.9145708468793e-311, 6.9145708468793e-311]) [ Info: ("case C (scalar loop, same target as B) [1:5] = ", [0.5499999999999998, 0.6499999999999999, 0.75, 0.8500000000000001, 0.9500000000000002]) ``` , noting that case B is wrong. This is because the `_fractional_indices` dispatch is wrong. For an `XYFlatGrid` source the dispatch is https://github.com/CliMA/Oceananigans.jl/blob/10c0f41db2245e3c753fc51248eb9d81b653505b/src/Fields/interpolate.jl#L206 which calls https://github.com/CliMA/Oceananigans.jl/blob/10c0f41db2245e3c753fc51248eb9d81b653505b/src/Fields/interpolate.jl#L210-L215 since `at_node` is a 3-tuple `(-144.9, 50.1, some_z)` but not `(nothing, nothing, some_z)`. I suppose there are two solutions to this. 1. We have to first address the question of whether we want to allow nonzero coordinates in `Flat` directions. Maybe there are some cases where the nonzero coordinates come into importance. Note that in the OS Papa example in Numerical Earth `x` and `y` directions are provided with nonzero values while they are `Flat` ```julia grid = RectilinearGrid(size = 200, x = λ★, y = φ★, z = (-400, 0), topology = (Flat, Flat, Bounded)) ``` (taken from https://github.com/NumericalEarth/NumericalEarth.jl/blob/fbee506ef71541e00745a6c4345035402d72671a/examples/single_column_os_papa_simulation.jl#L34-L38) If we think this is not the correct way to set up `grid`, we can throw an error at grid construction level. 2. We add additional `_fractional_indices` dispatch like ```julia @inline function _fractional_indices((x, y, z), grid, ::Nothing, ::Nothing, ℓz) kk = fractional_z_index(z, (nothing, nothing, ℓz), grid) return FractionalIndices(nothing, nothing, kk) end ``` However, having coordinates specified even on `Flat` directions could be useful for interpolation or computing properties like beta effect or albedo. So we should preserve this functionality. So we should go with option 2.