Interpolation and Coupling with Oceananigans

[skygering]

The current interpolation and points of ocean/atmosphere velocity/constants don't match with oceananigans. Oceananigans uses a c-grid. See image below to see what that means. Basically, the vector values are stored on half grid lines and constants are stored in the center of the grid. We will need to change that in order to match easily with Oceananigans. If the grid is Nx by Ny grid cells then the ocean.u field should be a [Nx+1, Ny] matrix, the ocean.v field should be a [Nx, Ny+1] matrix, and the ocean.temp field should be a [Nx, Ny] field.

[skygering]

Note that our current method matches with the MATLAB version of Subzero, so changing this would be a divergence from existing behavior and would invalidate our tests on this issue.

[skygering]

Also need to check that the stress from ice to ocean matches up with stress from ocean to ice and that we are calculating atmospheric stress correctly. Furthermore, the ocean stress needs to be on the same grids as the u and v velocities as mentioned above. This will require some refactoring.

[sdbrenner]

I want to revive discussion on this issue. I spent a little bit of time trying to figure out how to account for this. My current work-around is to interpolate the velocity fields from Oceananigans onto the Subzero grid when passing the velocities back and forth. This is made easy by the interpolate function in Oceananigans.Fields:

xi = iceSimulation.model.grid.x0:iceSimulation.model.grid.Δx:iceSimulation.model.grid.xf
yi = iceSimulation.model.grid.y0:iceSimulation.model.grid.Δy:iceSimulation.model.grid.yf
iceSimulation.model.ocean.u .= [interpolate((x,y,zoC[nz]),ocnModel.velocities.u) for x in xi, y in yi]
iceSimulation.model.ocean.v .= [interpolate((x,y,zoC[nz]),ocnModel.velocities.v) for x in xi, y in yi]

(though there may be a smarter approach)

Also note that Oceananigans also has built-in functions to interpolate onto either the Faces or Centers of ocean grid cells, so it's easy to put the velocities onto the same grid as each other; e.g.:

uccc = @at (Center, Center, Center) ocnModel.velocities.u
vccc = @at (Center, Center, Center) ocnModel.velocities.v

(but I don't think that the corners of the cell, where the Subzero grid is defined, is an option).

However, this doesn't solve the issue of going the other way: using the ocean stresses, τx and τy computed from Subzero as boundary conditions in Oceananigans. For now, these probably also need to be interpolated back to the Oceananigans grid. I think that interpolating to locations Center, Center is okay though, as the boundary condition can be defined as: Field{Center, Center, Nothing}(ocnGrid). But the interpolate function in Oceananigans.Fields works specifically on Oceananigans Fields, so we need a different interpolation approach.

I wonder if it would be useful to redefine the Subzero grid so that the ocean fields match with the Center, Center locations from Oceananigans, rather than the corners, to simplify the process (especially for the coupling fields, τx,τy).