Problem with boundary values in .core.grids._horizontal_gradient

[lesommer]

numpy.gradient compute the spatial derivatives with a second order scheme except near boundaries where the derivatives are evaluated with a first order scheme. oocgcm.core.grids._horizontal_gradient wraps numpy.gradient for xarray.DataArray with dask.

in order for the result not to depend on the chunk size, we need to set depth and use dask.array.map_overlap.

but it is not clear how to specify that the calculation should be first order near the boundary when the chunk is close to the boundary. numpy.gradient sees the ghost cell and therefore applies the second order derivative...

the current implementation explicitely sets a np.nan boundary condition to avoid errors.

this could be a question for dask mailing list.