Coriolis and Advection: vector invariant form

!-----------------------------------------------------------------------
&namdyn_adv    !   formulation of the momentum advection
!-----------------------------------------------------------------------
   ln_dynadv_vec = .true.  !  vector form (T) or flux form (F)
   nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction
   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection
/

The vector invariant form of the momentum equations (ln_dynhpg_vec = true) is the one most often used in applications of the NEMO ocean model. The flux form option (ln_dynhpg_vec =false) (see next section) has been present since version $2$. Options are defined through the namdyn_adv namelist variables. Coriolis and momentum advection terms are evaluated using a leapfrog scheme, $i.e.$ the velocity appearing in these expressions is centred in time (now velocity). At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following Chap.8.

Vorticity term (dynvor.F90)

!-----------------------------------------------------------------------
&namdyn_vor    !   option of physics/algorithm (not control by CPP keys)
!-----------------------------------------------------------------------
   ln_dynvor_ene = .false. !  enstrophy conserving scheme
   ln_dynvor_ens = .false. !  energy conserving scheme
   ln_dynvor_mix = .false. !  mixed scheme
   ln_dynvor_een = .true.  !  energy & enstrophy scheme
   ln_dynvor_een_old = .false.  !  energy & enstrophy scheme - original formulation
/

Options are defined through the namdyn_vor namelist variables. Four discretisations of the vorticity term (ln_dynvor_xxx=true) are available: conserving potential enstrophy of horizontally non-divergent flow (ENS scheme) ; conserving horizontal kinetic energy (ENE scheme) ; conserving potential enstrophy for the relative vorticity term and horizontal kinetic energy for the planetary vorticity term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy (EEN scheme) (see Appendix C.5). In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of vorticity term with analytical equations (ln_dynvor_con=true). The vorticity terms are all computed in dedicated routines that can be found in the dynvor.F90 module.

Enstrophy conserving scheme (ln_dynvor_ens=true)

In the enstrophy conserving case (ENS scheme), the discrete formulation of the vorticity term provides a global conservation of the enstrophy ( $[ (\zeta +f ) / e_{3f} ]^2$ in $s$-coordinates) for a horizontally non-divergent flow ($i.e.$ $\chi$=0), but does not conserve the total kinetic energy. It is given by:

$$\left\{ \begin{aligned}{+\frac{1}{e_{1u} } } & {\overline {\lef... ...e_{2u} e_{3u}\;u} \right)}} }^{ i+1/2, j} \end{aligned} \right.$$

Energy conserving scheme (ln_dynvor_ene=true)

The kinetic energy conserving scheme (ENE scheme) conserves the global kinetic energy but not the global enstrophy. It is given by:

$$\left\{ \begin{aligned}{+\frac{1}{e_{1u}}\; {\overline {\left( ... ... e_{3u}\;u} \right)} ^{ j+1/2}} }^{ i} } \end{aligned} \right.$$

Mixed energy/enstrophy conserving scheme (ln_dynvor_mix=true)

For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used. It consists of the ENS scheme (C.13) for the relative vorticity term, and of the ENE scheme (6.6) applied to the planetary vorticity term.

$$\left\{ { \begin{aligned}{+\frac{1}{e_{1u} }\; {\overline {\lef... ...;u} \right)} ^{ j+1/2}} }^{ i} } \hfill \end{aligned} } \right.$$

Energy and enstrophy conserving scheme (ln_dynvor_een=true)

In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$ averages of the velocity allows for the presence of grid point oscillation structures that will be invisible to the operator. These structures are computational modes that will be at least partly damped by the momentum diffusion operator ($i.e.$ the subgrid-scale advection), but not by the resolved advection term. The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field. Such noise would result in more noise in the vertical velocity field, an undesirable feature. This is a well-known characteristic of $C$-grid discretization where $u$ and $v$ are located at different grid points, a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid.

A very nice solution to the problem of double averaging was proposed by Arakawa and Hsu [1990]. The idea is to get rid of the double averaging by considering triad combinations of vorticity. It is noteworthy that this solution is conceptually quite similar to the one proposed by [Griffies et al., 1998] for the discretization of the iso-neutral diffusion operator (see App.C).

The Arakawa and Hsu [1990] vorticity advection scheme for a single layer is modified for spherical coordinates as described by Arakawa and Lamb [1981] to obtain the EEN scheme. First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:

$$q = \frac{\zeta +f} {e_{3f} }$$ (6.8)

where the relative vorticity is defined by (6.1), the Coriolis parameter is given by $f=2 \Omega \;\sin \varphi _f$ and the layer thickness at $f$-points is:

$$e_{3f} = \overline{\overline {e_{3t} }} ^{ i+1/2,j+1/2}$$ (6.9)

Figure 6.1: Triads used in the energy and enstrophy conserving scheme (een) for $u$-component (upper panel) and $v$-component (lower panel).

A key point in (6.9) is how the averaging in the i- and j- directions is made. It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks (ln_dynvor_een_old = false), or just by $4$ (ln_dynvor_een_old = true). The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and extends by continuity the value of $e_{3f}$ into the land areas. This case introduces a sub-grid-scale topography at f-points (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) that tends to reinforce the topostrophy of the flow ($i.e.$ the tendency of the flow to follow the isobaths) [Penduff et al., 2007].

Next, the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as the following triad combinations of the neighbouring potential vorticities defined at f-points (Fig. 6.1):

$$_i^j \mathbb{Q}^{i_p}_{j_p} = \frac{1}{12} \left( q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p} \right)$$ (6.10)

where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$.

Finally, the vorticity terms are represented as:

$$\left\{ { \begin{aligned}+q e_3 v &\equiv +\frac{1}{e_{1u} ... ...e_{3u} \;u \right)^{i+i_p}_{j+1/2-j_p} \end{aligned} } \right.$$

This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow ($i.e.$ $\chi$=0) (see Appendix C.5). Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of the noise in the vertical velocity field [Le Sommer et al., 2009]. Furthermore, used in combination with a partial steps representation of bottom topography, it improves the interaction between current and topography, leading to a larger topostrophy of the flow [Penduff et al., 2007, Barnier et al., 2006].

Kinetic Energy Gradient term (dynkeg.F90)

As demonstrated in Appendix C, there is a single discrete formulation of the kinetic energy gradient term that, together with the formulation chosen for the vertical advection (see below), conserves the total kinetic energy:

$$\left\{ \begin{aligned}-\frac{1}{2 \; e_{1u} } & \delta _{i+1... ...{u^2}^{ i} + \overline{v^2}^{ j}} \right] \end{aligned} \right.$$

Vertical advection term (dynzad.F90)

The discrete formulation of the vertical advection, together with the formulation chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic energy. Indeed, the change of KE due to the vertical advection is exactly balanced by the change of KE due to the gradient of KE (see Appendix C).

$$\left\{ \begin{aligned}-\frac{1} {e_{1u} e_{2u} e_{3u}} & \o... ...\;\delta _{k+1/2} \left[ u \right] }^{ k} \end{aligned} \right.$$

When ln_dynzad_zts = true, a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. This option can be useful when the value of the timestep is limited by vertical advection [Lemarié et al., 2015]. Note that in this case, a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, an option which is only available with a TVD scheme (see ln_traadv_tvd_zts in §5.1.2).