This is interesting. I like the concept of reusing the scheme kwarg. But I think this isn't sustainable given the future planned changes to the advection scheme API, right? Ie an "advection scheme" is going to be more comprehensive. Also, it doesn't make sense unless we are using vector invariant.

Another possibility is to add another kwarg.

Yes, I think we should probably put the high-order reconstructions in Operators

We can also do

advection = UpwindBiasedFifthOrder()
coriolis = HydrostaticSphericalCoriolis(scheme = advection)

Yes, I think we should probably put the high-order reconstructions in Operators

When we give an advection scheme to scheme, is there an assumption that we are using vector invariant? I'm still a little confused by the logic of the API. Maybe clarifying here will help?

The advection-like Coriolis scheme was a non-sequitur because f is a very regular field, so upwinding it was just decreasing performance without a significative increase in quality of the simulation.

On the other hand, upwinding u is very much discouraged because the energy builds up rapidly (by upwinding the velocity the divergence of the reconstructed tangential velocity is not a direct interpolation of the divergence of the original velocity, which is a necessary condition to maintain the algorithm stable).

The only thing I can think to increase the order of velocity interpolation in the Coriolis force is to use a centered high-order scheme to interpolate velocity, but that would not help with the noise since a centered scheme is dispersive in nature.

I converted this PR to implement a WetPointCoriolisScheme (described in Numerical boundary layers and spurious residual flows).
This is just a simple addition to an enstrophy conserving scheme where edge ("dry") points are neglected in the interpolation of the velocity in the tangential direction.

A comparison of the output of this scheme in a global 1 degree setup will follow

additionally a very small but significant bug has been corrected in the GM tapering, where, if we had bz == 0 the tapering would become Inf.

This is very rare occasion but it was happening with some particular initial conditions

Comparison with the default HydrostaticSphericalCoriolis() on main (left) and HydrostaticSphericalCoriolis(scheme = WetPointEnstrophyConservingScheme()) (right)
The figure below is the zonal velocity (u, top) and the vertical velocity (w, bottom) at the surface averaged over a year after 20 years of evolution from ECCO interpolated initial conditions.
Note the strikingly lower noise emanating from the boundaries (especially in the w velocity) using the new coriolis scheme

Also perhaps merge main because I think 0.77.3 is the current version.