Accurate calculation of pressure forces on cells defined by steeply sloping coordinates

The NEMO ocean model has a fairly generic s-level vertical co-ordinate option that allows the model levels to follow either the ocean bathymetry resolved by the model or a smoothed ``envelope'' version of it. Sufficiently accurate calculation of the net horizontal pressure forces on the resulting steeply sloping cells is well known to be difficult and many schemes have attempted to address this issue. We focus here on two types of schemes that calculate these forces by integrating around the cell perimeters (or faces). The first type calculates the pressure forces on the faces of each velocity cell as first proposed by Lin (1997). The second type is a re-arrangement of the Song (1998) density Jacobian scheme derived by Shchepetkin & McWilliams (2003). In these schemes the density and the height along coordinate surfaces are approximated by local piecewise continuous polynomial reconstructions and the integrals are calculated by direct integration of the polynomials or by quadrature methods. Higher order schemes are obtained by using higher order polynomials. In some cases these polynomials are constrained to avoid over-fitting of potentially noisy fields. We have implemented these schemes in NEMO treating model variables as point values rather than cell-mean values. The implementation for the second type follows previously published work closely whilst that for the first type differs somewhat from previous work. We have also implemented a ``pre-processing'' step that calculates cubic reconstructions of the density profile at the deepest point in the stencil used to calculate the net horizontal forces. The ``pre-processing'' step subtracts this interpolated profile from all the other profiles in the stencil or uses it to assist the interpolation on surfaces of constant s. We present results for these schemes in a standard sea-mount test case with steep bathymetry and horizontal density surfaces in which the currents should be identically zero. We show that both types of schemes using appropriate quadratic or cubic reconstructions give currents that are about one order of magnitude smaller than those given by second order accurate schemes. This is the case even if the fields are initialised as cell-mean values rather than point values. Using ``pre-processing'' steps, one can obtain smaller initial velocities but small velocities are not retained in longer integrations.