Surface pressure gradient (dynspg.F90)

!-----------------------------------------------------------------------
!namdyn_spg    !   surface pressure gradient   (CPP key only)
!-----------------------------------------------------------------------
!                          !  explicit free surface                     ("key_dynspg_exp")
!                          !  filtered free surface                     ("key_dynspg_flt")
!                          !  split-explicit free surface               ("key_dynspg_ts")

$$

Options are defined through the namdyn_spg namelist variables. The surface pressure gradient term is related to the representation of the free surface (§2.2). The main distinction is between the fixed volume case (linear free surface) and the variable volume case (nonlinear free surface, key_ vvl is defined). In the linear free surface case (§2.2.2) the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case (§2.2.2). With both linear and nonlinear free surface, external gravity waves are allowed in the equations, which imposes a very small time step when an explicit time stepping is used. Two methods are proposed to allow a longer time step for the three-dimensional equations: the filtered free surface, which is a modification of the continuous equations (see ()), and the split-explicit free surface described below. The extra term introduced in the filtered method is calculated implicitly, so that the update of the next velocities is done in module dynspg_flt.F90 and not in dynnxt.F90.

The form of the surface pressure gradient term depends on how the user wants to handle the fast external gravity waves that are a solution of the analytical equation (§2.2). Three formulations are available, all controlled by a CPP key (ln_dynspg_xxx): an explicit formulation which requires a small time step ; a filtered free surface formulation which allows a larger time step by adding a filtering term into the momentum equation ; and a split-explicit free surface formulation, described below, which also allows a larger time step.

The extra term introduced in the filtered method is calculated implicitly, so that a solver is used to compute it. As a consequence the update of the $next$ velocities is done in module dynspg_flt.F90 and not in dynnxt.F90.

Explicit free surface (key_ dynspg_exp)

In the explicit free surface formulation (key_ dynspg_exp defined), the model time step is chosen to be small enough to resolve the external gravity waves (typically a few tens of seconds). The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time), is thus simply given by :

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

Note that in the non-linear free surface case ($i.e.$ key_ vvl defined), the surface pressure gradient is already included in the momentum tendency through the level thickness variation allowed in the computation of the hydrostatic pressure gradient. Thus, nothing is done in the dynspg_exp.F90 module.

Split-Explicit free surface (key_ dynspg_ts)

!-----------------------------------------------------------------------
&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
!-----------------------------------------------------------------------
   ln_bt_fw      =    .TRUE.           !  Forward integration of barotropic equations
   ln_bt_av      =    .TRUE.           !  Time filtering of barotropic variables
   ln_bt_nn_auto =    .TRUE.           !  Set nn_baro automatically to be just below
                                       !  a user defined maximum courant number (rn_bt_cmax)
   nn_baro       =    30               !  Number of iterations of barotropic mode
                                       !  during rn_rdt seconds. Only used if ln_bt_nn_auto=F
   rn_bt_cmax    =    0.8              !  Maximum courant number allowed if ln_bt_nn_auto=T
   nn_bt_flt     =    1                !  Time filter choice
                                       !  = 0 None
                                       !  = 1 Boxcar over   nn_baro barotropic steps
                                       !  = 2 Boxcar over 2*nn_baro     "        "
/

The split-explicit free surface formulation used in NEMO (key_ dynspg_ts defined), also called the time-splitting formulation, follows the one proposed by Shchepetkin and McWilliams [2005]. The general idea is to solve the free surface equation and the associated barotropic velocity equations with a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables (Fig. 6.2). The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through the nn_baro namelist parameter as: $\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (ln_bt_nn_auto=true) considering that the stability of the barotropic system is essentially controled by external waves propagation. Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry. Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than rn_bt_cmax.

The barotropic mode solves the following equations:

$$\begin{equation}\frac{\partial {\rm\overline{{\bf U}}_h} }{\parti... ...( {H+\eta } \right) \; {\rm {\bf\overline{U}}}_h } \right]+P-E \end{equation}$$

where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. The third term on the right hand side of (6.25a) represents the bottom stress (see section §10.4), explicitly accounted for at each barotropic iteration. Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm detailed in Shchepetkin and McWilliams [2005]. AB3-AM4 coefficients used in NEMO follow the second-order accurate, "multi-purpose" stability compromise as defined in Shchepetkin and McWilliams [2008] (see their figure 12, lower left).

Figure 6.2: Schematic of the split-explicit time stepping scheme for the external and internal modes. Time increases to the right. In this particular exemple, a boxcar averaging window over $nn\_baro$ barotropic time steps is used ( $nn\_bt\_flt=1$) and $nn\_baro=5$. Internal mode time steps (which are also the model time steps) are denoted by $t-\rdt$, $t$ and $t+\rdt$. Variables with $k$ superscript refer to instantaneous barotropic variables, $< >$ and $« »$ operator refer to time filtered variables using respectively primary (red vertical bars) and secondary weights (blue vertical bars). The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged transports to advect tracers. a) Forward time integration: ln_bt_fw=true, ln_bt_av=true. b) Centred time integration: ln_bt_fw=false, ln_bt_av=true. c) Forward time integration with no time filtering (POM-like scheme): ln_bt_fw=true, ln_bt_av=false.

In the default case (ln_bt_fw=true), the external mode is integrated between now and after baroclinic time-steps (Fig. 6.2a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic quantities (ln_bt_av=true). In that case, the integration is extended slightly beyond after time step to provide time filtered quantities. These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, asselin filtering is not applied to barotropic quantities.
Alternatively, one can choose to integrate barotropic equations starting from before time step (ln_bt_fw=false). Although more computationaly expensive ( nn_baro additional iterations are indeed necessary), the baroclinic to barotropic forcing term given at now time step become centred in the middle of the integration window. It can easily be shown that this property removes part of splitting errors between modes, which increases the overall numerical robustness.

As far as tracer conservation is concerned, barotropic velocities used to advect tracers must also be updated at now time step. This implies to change the traditional order of computations in NEMO: most of momentum trends (including the barotropic mode calculation) updated first, tracers' after. This de facto makes semi-implicit hydrostatic pressure gradient (see section §6.4.5) and time splitting not compatible. Advective barotropic velocities are obtained by using a secondary set of filtering weights, uniquely defined from the filter coefficients used for the time averaging (Shchepetkin and McWilliams [2005]). Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to obtain exact conservation.

One can eventually choose to feedback instantaneous values by not using any time filter (ln_bt_av=false). In that case, external mode equations are continuous in time, ie they are not re-initialized when starting a new sub-stepping sequence. This is the method used so far in the POM model, the stability being maintained by refreshing at (almost) each barotropic time step advection and horizontal diffusion terms. Since the latter terms have not been added in NEMO for computational efficiency, removing time filtering is not recommended except for debugging purposes. This may be used for instance to appreciate the damping effect of the standard formulation on external gravity waves in idealized or weakly non-linear cases. Although the damping is lower than for the filtered free surface, it is still significant as shown by Levier et al. [2007] in the case of an analytical barotropic Kelvin wave.

Filtered free surface (key_ dynspg_flt)

The filtered formulation follows the Roullet and Madec [2000] implementation. The extra term introduced in the equations (see §2.2.2) is solved implicitly. The elliptic solvers available in the code are documented in §15.

Note that in the linear free surface formulation (key_ vvl not defined), the ocean depth is time-independent and so is the matrix to be inverted. It is computed once and for all and applies to all ocean time steps.