29.6.3 Time stepping algorithm

The focus in this subsection is on time and depth discretization, with Figure 29.3 summarizing the following algorithm. For purposes of brevity, the horizontal spatial discretization discussed in the previous subsection will not be exposed. Also, discrete baroclinic times and time steps will be denoted by the Greek $\tau$ and $\Delta\tau$, respectively, whereas the barotropic analogs will use the Latin t and $\Delta t$.

The basic idea is to split the velocity at an arbitrary depth level k and baroclinic time $\tau' = \tau + \Delta \tau$ into two components

 

$$u_{k}(\tau') B_{km}(\tau) \, u_{m}(\tau') + (\delta_{km} - B_{km}(\tau) ) \, u_{m}(\tau')$$ =

$$\equiv \hat{u}_{k}(\tau,\tau') + \overline{u}(\tau,\tau').$$ (29.86)

The following ``baroclinicity operator'' is used to affect this split

$$B_{km}(\tau) = \delta_{km} - D(\tau)^{-1} \, h^{u}_{m}(\tau),$$ (29.87)

where $\delta_{km}$ is the Kronecker delta, summation over the repeated vertical level index m is implied, and

$$D(\tau) = \sum_{k} h^{u}_{k}(\tau) = D_{o} + \eta^{u}(\tau)$$ (29.88)

is the ocean depth at baroclinic time $\tau$ over a column of velocity points, with D_o the resting ocean depth.

The introduction of two baroclinic time labels to equation (29.87) is necessitated by the freedom afforded the ocean depth to change in time. This property is in contrast to the case with a rigid lid, or the fixed volume free surface models of Killworth et al. (1991) and Dukowicz and Smith (1994).

Equation (29.87) is an identity which is valid for any baroclinic times $\tau'$ and $\tau$. Its utility depends on the ability to render a relatively clean and stable split between the fast and slow dynamics. The form of the baroclinicity operator $B_{km}(\tau)$ is motivated by an attempt to perform such a split. That is, it projects out the approximate baroclinic portion of a field, where the projection is based on the distribution of cell thicknesses at time $\tau$.

If the split introduced in equation (29.87) is successful, the baroclinic velocity $\hat{u}_{k}(\tau,\tau')$ will evolve on a slow time scale $\Delta\tau$. In turn, the barotropic velocity $\overline{u}(\tau,\tau')$ will evolve on the fast time scale $\Delta t = 2 \, N^{-1} \, \Delta \tau$, with N determined by the ratio of external to internal gravity wave speeds ( $N \approx 100$ for climate models). The method therefore proceeds by separately updating $\hat{u}_{k}(\tau,\tau')$ and $\overline{u}(\tau,\tau')$ by exploiting the time scale split. Upon doing so, the right hand side of the identity (29.87) will be specified, hence allowing for an update of the full velocity field $u_{k}(\tau+\Delta \tau)$. The following discussion details these ideas.

  

Figure 29.3: Schematic of the split-explicit time stepping scheme. Time increases to the right. The baroclinic time steps are denoted by $\tau -\Delta \tau , \tau , \tau + \Delta \tau$, and $\tau +2\Delta \tau$. The curved line represents a baroclinic leap-frog time step, and the smaller barotropic time steps $N \, \Delta t = 2 \, \Delta \tau$ are denoted by the zig-zag line. First, a baroclinic leap-frog time step updates the baroclinic mode to $\tau +\Delta \tau$. Then, using the vertically integrated forcing $G(\tau )$ (equation 29.86) computed at baroclinic time step $\tau$, a forward-backward time stepping scheme integrates the surface height and vertically integrated velocity from $\tau$ to $\tau +2\Delta \tau$ using Nbarotropic time steps of length $\Delta t$, while keeping $G(\tau )$and the ocean depth $D(\tau )$ fixed. Time averaging the barotropic fields over the N+1 time steps (endpoints included) centers the vertically integrated velocity and free surface height at baroclinic time step $\tau +\Delta \tau$. Note that for accurate centering, Nis set to an even integer.

By construction, evolution of the baroclinic mode is unaffected by vertically independent forces, such as those from surface pressure gradients. Therefore, it is sufficient to update the ``primed'' velocity

$$u_{k}'(\tau + \Delta \tau) = u_{k}^{R}(\tau - \Delta \tau) + 2 \Delta \tau \, [f \, v_{k}(\tau) + \tilde{G}_{k}(\tau)],$$ (29.89)

which represents a temporal discretization of the full momentum equation (29.83), yet without the surface pressure gradient. The lagged velocity $u_{k}^{R}(\tau - \Delta \tau)$ is a Robert time filtered version of the full velocity field. A weak form of such filtering has been found sufficient to suppress a splitting between the two branches of the leap-frog (e.g., Haltiner and Williams 1980). It is also useful to dampen any fast dynamics which may partially leak through the baroclinicity operator due to the imperfect baroclinic/barotropic split discussed in Section 5.1.1. The baroclinic piece of the primed velocity $u_{k}'(\tau + \Delta \tau)$ is equivalent to the updated baroclinic velocity

$$\hat{u}_{k}(\tau, \tau + \Delta \tau) = B(\tau)_{km} \, u_{m}'(\tau + \Delta \tau),$$ (29.90)

thus specifying one half of the identity (29.87). To update the barotropic velocity $\overline{u}(\tau,\tau+\Delta \tau)$, we employ a forward-backward time stepping scheme (e.g., Haltiner and Williams 1980). For this scheme, the surface height is time stepped using the small barotropic time step $\Delta t$

 

$$\eta^{*}(t_{n+1}) = \eta^{*}(t_{n}) - \Delta t \, [\nabla_{h} \cdot {\bf U}^{*}(t_{n}) - q_{w}(\tau) ],$$ (29.91)

with $t_{n} = n\,\Delta t$ the barotropic time, and $n \in [0,N]$. The asterisk is used to denote intermediate values of the barotropic fields, each of which are updated on the barotropic time steps. For stability purposes, it is important to take the initial condition $\eta^{*}(0)$ as the time average $\overline{\eta}^{*}$ from the previous barotropic integration (time averaging is defined by equation (29.95) discussed below). Note that it is assumed that the fresh water flux q_w is constant over the small barotropic time steps, since the hydrological fluxes are typically updated at a period no shorter than the baroclinic time step.

Having the surface height $\eta^{*}$ updated to the new barotropic time step allows for an update of the transport

 

$$U^{*}(t_{n}) + f \, \Delta t \, \left( \frac{V^{*}(t_{n})+V^{*}(t_{n+1})}{2} \right)$$ U^*(t_n+1) =

$$\Delta t \, [G(\tau) - D(\tau) \, \partial_{x} \, \tilde{p}^{*}_{s}(t_{n+1})],$$ + (29.92)

where $\tilde{p}^{*}_{s}(t_{n+1}) = g \, \eta^{*}(t_{n+1}) \, \rho(\tau) / \rho_{o}$ is the surface pressure normalized by the Boussinesq density. Both the ocean depth $D(\tau )$ and depth integrated forcing $G(\tau )$ are assumed to evolve on the baroclinic time scale, and so are held constant over the extent of the barotropic time steps from $\tau$ to $\tau +2\Delta \tau$. The Coriolis force is computed using a Crank-Nicholson semi-implicit time stepping scheme (e.g., Haltiner and Williams 1980).

After N barotropic time steps, the vertical transport and surface height are time averaged to produce

  

$$\overline{U}^{*}(\tau+\Delta \tau) \frac{1}{N+1} \, \sum_{n=0}^{N} \, U^{*}(t_{n})$$ = (29.93)

$$\overline{\eta}^{*}(\tau+\Delta \tau) \frac{1}{N+1} \, \sum_{n=0}^{N} \, \eta^{*}(t_{n}).$$ = (29.94)

The time averaged fields are centered on the baroclinic time step $\tau +\Delta \tau$, so long as N is an even integer. Equating the barotropic velocity $\overline{u}(\tau,\tau+\Delta \tau)$ to $\overline{U}^{*}(\tau+\Delta \tau)/D(\tau)$ allows for the full velocity field $u(\tau+\Delta \tau)$ to be updated according to equation (29.87). The time averaged field $\overline{\eta}^{*}(\tau+\Delta \tau)$ is used to initialize the next suite of barotropic integrations from $\tau +\Delta \tau$ to $\tau+3\Delta \tau$.

To begin the next baroclinic time step, it is necessary to determine the updated surface height $\eta(\tau+\Delta \tau)$ and vertically integrated velocity $U(\tau + \Delta \tau)$. A natural choice is to equate these fields to the time averages of the intermediate fields $\eta^{*}$ and U^* given by equations (29.94) and (29.95). Unfortunately, this choice does not lead to a self-consistent and conservative algorithm. The reason is that it is essential to maintain consistency between the updated ocean depth, full velocity, barotropic velocity, all while maintaining basic conservation properties. The following approach has been found to be sufficient for these purposes.

As discussed in Section 29.6.9, since the tracer concentration is time stepped with a leap-frog scheme, quasi-conservation of tracer results if the surface height is similarly time stepped

$$\eta(\tau+\Delta \tau) = \overline{\eta}^{*}(\tau-\Delta \t... ...\Delta \tau \, [\nabla_{h} \cdot {\bf U}(\tau) - q_{w}(\tau)].$$ (29.95)

To maintain stability and smoothness of the solutions, it has been found necessary to use the lagged surface height $\overline{\eta}^{*}(\tau-\Delta \tau)$ rather than a more traditional Robert filtered height. As so defined, $\eta(\tau+\Delta \tau)$ is used to update thicknesses of the surface grid cells $h^{t}_{k}(\tau + \Delta \tau)$ and $h^{u}_{k}(\tau + \Delta \tau)$. Given these new thicknesses, the updated transport $U(\tau + \Delta \tau)$ can be diagnosed from the known full velocity $u_{k}(\tau+\Delta \tau)$through

$$U(\tau + \Delta \tau) = \sum_{k} h^{u}_{k}(\tau + \Delta \tau) \, u_{k}(\tau + \Delta \tau),$$ (29.96)

thus ensuring self-consistency between the updated full velocity and the updated vertically integrated velocity. That is, with time dependent thicknesses, $U(\tau + \Delta \tau)$ differs from the time averaged field $\overline{U}^{*}(\tau+\Delta \tau)$ because of the changing level thicknesses. To complete the time stepping, the updated barotropic velocity is diagnosed through

$$\overline{u}(\tau + \Delta \tau) = U(\tau + \Delta \tau)/D(\tau + \Delta \tau),$$ (29.97)

where $D(\tau + \Delta \tau) = \sum h^{u}_{k}(\tau + \Delta \tau)$ is the new depth of a velocity cell column.