(0.97.6) Move the computation of flux BC before the time-stepping

[simone-silvestri]

This PR is moving the computation of the boundary conditions just before the time-stepping of the prognostic variables, rather than at the end of the time-stepping procedure.

This is because, in Oceananigans, we add the flux bcs to the tendencies, which are computed at the end of a time step. However, if we need to couple a model, the exchange between components is typically computed at the end of the stepping of all the components, meaning that, if we do this, we do not add the updated fluxes to the tendency in an Oceananigans model; rather, the fluxes corresponding to the previous time-step.

This might become an issue when coupling on large time steps, where the fluxes lag behind and the model can overshoot.
This change will not affect oceananigans-only models, but will affect coupled models like the OceanSeaIceModel in ClimaOcean.jl.

[navidcy]

This is because, in Oceananigans, we add the flux bcs to the tendencies, which are computed at the end of a time step.

So, how was this working for multi-step timestepping schemes? Before this PR the flux bcs were added at the end of the Δt or at the end of each substep?

[simone-silvestri]

@navidcy This is the difference between main and this PR in a one-step method like AB2 blue boxes are computed in ClimaOcean, black boxes in Oceanananigans.
RK3 is a little more complex but still has only one blue box per complete Oceananigans timestep.

[simone-silvestri]

This is a bit oversemplified, but it is the RK3 timestepper

If we want to compute fluxes each substep we probably need a timestepper for the OceanSeaIceModel which would probably help a lot with stabiliy as we could substep the sea ice with the ocean