Hydrostatic pressure gradient (dynhpg.F90)

!-----------------------------------------------------------------------
&namdyn_hpg    !   Hydrostatic pressure gradient option
!-----------------------------------------------------------------------
   ln_hpg_zco  = .false.   !  z-coordinate - full steps
   ln_hpg_zps  = .true.    !  z-coordinate - partial steps (interpolation)
   ln_hpg_sco  = .false.   !  s-coordinate (standard jacobian formulation)
   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to isf
   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
                                 !           centered      time scheme  (F)
/

Options are defined through the namdyn_hpg namelist variables. The key distinction between the different algorithms used for the hydrostatic pressure gradient is the vertical coordinate used, since HPG is a horizontal pressure gradient, $i.e.$ computed along geopotential surfaces. As a result, any tilt of the surface of the computational levels will require a specific treatment to compute the hydrostatic pressure gradient.

The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, $i.e.$ the density appearing in its expression is centred in time (now $\rho$), or a semi-implcit scheme. At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied.

$z$-coordinate with full step (ln_dynhpg_zco=true)

The hydrostatic pressure can be obtained by integrating the hydrostatic equation vertically from the surface. However, the pressure is large at great depth while its horizontal gradient is several orders of magnitude smaller. This may lead to large truncation errors in the pressure gradient terms. Thus, the two horizontal components of the hydrostatic pressure gradient are computed directly as follows:

for $k=km$ (surface layer, $jk=1$ in the code)

$$\left\{ \begin{aligned}\left. \delta _{i+1/2} \left[ p^h \right... ..._{3w} \rho \right] \right\vert _{k=km} \end{aligned} \right.$$

for $1<k<km$ (interior layer)

$$\left\{ \begin{aligned}\left. \delta _{i+1/2} \left[ p^h \right... ... {\rho}^{k+1/2} \right] \right\vert _{k} \end{aligned} \right.$$

Note that the $1/2$ factor in (6.18) is adequate because of the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface level ($z=0$). Note also that in case of variable volume level (key_ vvl defined), the surface pressure gradient is included in (6.18) and (6.19) through the space and time variations of the vertical scale factor $e_{3w}$.

$z$-coordinate with partial step (ln_dynhpg_zps=true)

With partial bottom cells, tracers in horizontally adjacent cells generally live at different depths. Before taking horizontal gradients between these tracer points, a linear interpolation is used to approximate the deeper tracer as if it actually lived at the depth of the shallower tracer point.

Apart from this modification, the horizontal hydrostatic pressure gradient evaluated in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case. As explained in detail in section §5.9, the nonlinearity of pressure effects in the equation of state is such that it is better to interpolate temperature and salinity vertically before computing the density. Horizontal gradients of temperature and salinity are needed for the TRA modules, which is the reason why the horizontal gradients of density at the deepest model level are computed in module zpsdhe.F90 located in the TRA directory and described in §5.9.

$s$- and $z$-$s$-coordinates

Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers ($e.g.$, Shchepetkin and McWilliams [2005], Song [1998]). A number of different pressure gradient options are coded but the ROMS-like, density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.

$\bullet$ Traditional coding (see for example Madec et al. [1996]: (ln_dynhpg_sco=true)

$$\left\{ \begin{aligned}- \frac{1} {\rho_o e_{1u}} \; \delta ... ...}} \; \delta _{j+1/2} \left[ z_t \right] \end{aligned} \right.$$

Where the first term is the pressure gradient along coordinates, computed as in (6.18) - (6.19), and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point ($e_{3w}$).

$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (ln_dynhpg_prj=true)

$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) [Shchepetkin and McWilliams, 2005] (ln_dynhpg_djc=true) (currently disabled; under development)

Note that expression (6.20) is commonly used when the variable volume formulation is activated (key_ vvl) because in that case, even with a flat bottom, the coordinate surfaces are not horizontal but follow the free surface [Levier et al., 2007]. The pressure jacobian scheme (ln_dynhpg_prj=true) is available as an improved option to ln_dynhpg_sco=true when key_ vvl is active. The pressure Jacobian scheme uses a constrained cubic spline to reconstruct the density profile across the water column. This method maintains the monotonicity between the density nodes The pressure can be calculated by analytical integration of the density profile and a pressure Jacobian method is used to solve the horizontal pressure gradient. This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme.

Ice shelf cavity

Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and the pressure gradient due to the ocean load. If cavities are present (ln_isfcav = true) these two terms can be calculated by setting ln_dynhpg_isf = true. No other scheme is working with ice shelves.

$\bullet$ The main hypothesis to compute the ice shelf load is that the ice shelf is in isostatic equilibrium. The top pressure is computed integrating a reference density profile (prescribed as density of a water at 34.4 PSU and -1.9$^{\circ}C$) from the sea surface to the ice shelf base, which corresponds to the load of the water column in which the ice shelf is floatting. This top pressure is constant over time. A detailed description of this method is described in Losch [2008].

$\bullet$ The ocean load is computed using the expression (6.20) described in 6.4.3. A treatment of the top and bottom partial cells similar to the one described in 6.4.2 is done to reduce the residual circulation generated by the top partial cell.

Time-scheme (ln_dynhpg_imp= true/false)

The default time differencing scheme used for the horizontal pressure gradient is a leapfrog scheme and therefore the density used in all discrete expressions given above is the now density, computed from the now temperature and salinity. In some specific cases (usually high resolution simulations over an ocean domain which includes weakly stratified regions) the physical phenomenon that controls the time-step is internal gravity waves (IGWs). A semi-implicit scheme for doubling the stability limit associated with IGWs can be used [Maltrud et al., 1998, Brown and Campana, 1978]. It involves the evaluation of the hydrostatic pressure gradient as an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ ($i.e.$ before, now and after time-steps), rather than at the central time level $t$ only, as in the standard leapfrog scheme.

$\bullet$ leapfrog scheme (ln_dynhpg_imp=true):

$$\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; -\frac{1}{\rho _o e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right]$$ (6.20)

$\bullet$ semi-implicit scheme (ln_dynhpg_imp=true):

$$\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; -\frac{1}{4 \r... ...o e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2 p_h^t +p_h^{t-\rdt} \right]$$ (6.21)

The semi-implicit time scheme (6.22) is made possible without significant additional computation since the density can be updated to time level $t+\rdt$ before computing the horizontal hydrostatic pressure gradient. It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using (6.22) compared to that using the standard leapfrog scheme (6.21). Note that (6.22) is equivalent to applying a time filter to the pressure gradient to eliminate high frequency IGWs. Obviously, when using (6.22), the doubling of the time-step is achievable only if no other factors control the time-step, such as the stability limits associated with advection or diffusion.

In practice, the semi-implicit scheme is used when ln_dynhpg_imp=true. In this case, we choose to apply the time filter to temperature and salinity used in the equation of state, instead of applying it to the hydrostatic pressure or to the density, so that no additional storage array has to be defined. The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows:

$$\rho^t = \rho( \widetilde{T},\widetilde {S},z_t) \quad \widetilde{X} = 1 / 4 \left( X^{t+\rdt} +2 X^t + X^{t-\rdt} \right)$$ (6.22)

Note that in the semi-implicit case, it is necessary to save the filtered density, an extra three-dimensional field, in the restart file to restart the model with exact reproducibility. This option is controlled by nn_dynhpg_rst, a namelist parameter.