Change `MutableVerticalDiscretization` to `GeneralizedVerticalDiscretization` and implement horizontal derivatives

[glwagner]

I think that MutableVerticalDiscretization is misnamed a bit. Instead it should be called GeneralizedVerticalDiscretization, because this is a general (non-orthogonal) vertical coordinate implementation.

Also, horizontal derivatives are currently incorrect when using MutableVerticalDiscretization. For example from the CM1 documentation:

In the above, $\sigma$ is the generalized vertical coordinate (what we call $r$).

In Oceananigans, metric terms are implemented in a non-generic way for the hydrostatic pressure gradient and z-star specifically, eg

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/z_star_vertical_spacing.jl

Lines 219 to 220 in d695969

	@inline grid_slope_contribution_x(i, j, k, grid::MutableGridOfSomeKind, buoyancy, ::ZStarCoordinate, model_fields) =
	ℑxᶠᵃᵃ(i, j, k, grid, buoyancy_perturbationᶜᶜᶜ, buoyancy.formulation, model_fields) * ∂x_z(i, j, k, grid)

the non-generic implementation has a few issues:

• While the dynamics are correct, the diagnostics are wrong; more generally anywhere horizontal derivatives are invoked are wrong
• It produces less code-reuse if we are to implement other coordinates other than z-star.

I propose correctly implementing horizontal derivatives for GeneralizeVerticalDiscretization. We will leave the "difference operators" unchanged, and only change operators like ∂x.

There also seems to be confusion about the notation. For example we have an operator ∂x_z. But if z is height, this is 0. It should read ∂x_r, so that we can interpret z as height.

@simone-silvestri

[simone-silvestri]

Yeah, we can add the horizontal derivatives and keep the differences as is. I think the horizontal derivative is only used for the hydrostatic pressure and the closures in the dynamics?