This is a work-in-progress PR that aims to reformulate the timestepping scheme of the HydrostaticFreeSurfaceModel to
make sure that three numerical properties are satisfied:

• tracers are globally conserved (the integral of a tracer field remains constant in the absence of forcings)
• tracers are locally conserved (an initially constant tracer remains constant in the absence of forcings)
• The free surface is evolved consistently with the grid thickness.

This is easier to achieve with an explicit discretization of the free surface (ExplicitFreeSurface and SplitExplicitFreeSurface) rather than an implicit free surface discretization and with RK3 rather than AB2.
However, this PR reformulates much of the timestepping algorithm, so it will be a lengthy WIP that will need to be thoroughly validated for all combinations of barotropic and baroclinic timestepping schemes.

More changes are performed to correctly abstract the timestepping schemes, which, at the moment, are tailored specifically to the NonhydrostaticModel (for example, the timesteppers include a compute_pressure_correction! or a correct_velocities function, which are meaningful only for a non-hydrostatic Navier Stokes solver)

Edit: another objective I have is to remove any timestepping from the update_state! function. We should be able to call update_state! multiple times in a row without problems and timestepping inside update_state! negates this possibility.

TODO:

• Remove timestepping from update_state! In another PR
• Clean up the time-stepping code
• Introduce transport weights for transport averaging in the SplitExplicitFreeSurface
• Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the Explicit free surface
• Implement a conservative, low-storage Runge Kutta method (of arbitrary order) coupled with the SplitExplicitFreeSurface
• Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the (PCG) ImplicitFreeSurface
• Ensure that the conservation is satisfied on DistributedGrids
• Implement a semi - conservative method for AB2 In another PR
• Address precision errors that tend to destroy conservation for constant fields In another PR