Implicit free surface implementation in `NonhydrostaticModel`

[glwagner]

This issue documents the development needed to implement an implicit free surface formulation in NonhydrostaticModel, which is prototyped by #3968 and #4166

The implicit free surface algorithm boils down to a modified pressure solve which invokes a Robin boundary condition on pressure, as opposed to the homogeneous Neumann boundary condition that's used for a rigid lid formulation. In addition to implementing the modified boundary condition and pressure solver, we also update the free surface elevation. Here's a summary of the modified algorithm:

1. first solve_for_pressure!, incorporating the current Δt and also gravitational acceleration here

2. fill_halo_regions! on pressure, using the Robin boundary condition, which depends on Δt and g. The modified boundary condition is needed to compute the vertical pressure gradient across the top boundary. This is protoyped here

3. Call make_pressure_correction!(model, Δt) as usual.

4. Call update_free_surface!(model, Δt) which effectively executes this line

The todo list is therefore

• Implement RobinBoundaryCondition (CombinationBoundaryCondition to complement Gradient and Value boundary conditions #4490)
• Develop the specific RobinBoundaryCondition needed for pressure with the implicit free surface formulation.
  • Note: filling halo regions on pressure will then require Δt and g
• (Modify NonhydrostaticModel to have a free surface property, [as in original PR(https://github.com/shriyafruitwala/Oceananigans.jl/blob/e3a02e1ce5cfdb742335e07b510c9b1a27e1fb19/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl#L223-L226))
• Implement a modified FourierTridigonalPoissonSolver where diagonals can be updated each time step dependent on the passed parameters Δt and g

[glwagner]

cc @joernc, @shriyafruitwala