35.1.8.6 Linear numerical stability for the RM98 operator

For linear stability, it is sufficient, and conservative, to consider the numerical stability issues raised by Cox (1987) and Griff ies et al. (1998) isoneutral diffusion papers. In particular, one requires

$$\frac{(B \, \vert \nabla_{h}^{2} \, S_{x} \vert) \, \Delta t} { \Delta x \, \Delta z} \le \frac{1}{4}$$ (35.59)

$$\frac{(B \, \vert\nabla_{h}^{2} \, S_{y}\vert) \, \Delta t} { \Delta y \, \Delta z} \le \frac{1}{4}.$$ (35.60)

For a conservative estimate of the stability, introduce the dimensionless grid factor

$$\delta = \frac{\Delta z \; (\Delta a)^{3} } {4 \, B \, \Delta t},$$ (35.61)

where $\Delta a$ is the minimum horizontal grid spacing over the extent of the model domain, and $\Delta z$ is the minimum vertical spacing. Note that $\Delta a$ should take into account the convergence of the meridions. With this notation, the stability constraint takes the form

$$\le \delta$$ |S_x| (35.62)

$$\le \delta.$$ |S_y| (35.63)

As with the isoneutral diffusion scheme, these constraints place limits on the maximum isoneutral slope that can be realized without introducing some form of tapering to the biharmonic coefficient B. When either component of the slope vector has a magnitude larger than $\delta$, then the tapering of dm_taper (section 34.1.9.1) or gkw_taper (Section 34.1.9.2; the default) is employed.

An example is useful. Consider a typical mid-latitude channel model with

$$\Delta a \cos(\pi/4) \, (.25^{\circ})$$ = (35.64)

$$\Delta z$$ = 10m (35.65)

$$\Delta t$$ = 1000 secs (35.66)

B = 10¹⁹ cm⁴/sec. (35.67)

Linear stability says that tapering of the biharmonic coefficient must occur when the slope is larger than

$$\delta = 75/100.$$ (35.68)

In the ocean, this is a rather steep slope, and so should not place a serious constraint on the regions for which the biharmonic dissipation acts with the full diffusivity. Note that the use of B = 10¹⁹ cm⁴/sec implies a dissipation time scale of grid sized tracer anomalies (see Section 33.4 for details of how this damping time is derived)

$$\tau = (\Delta a / 2)^{4} / B$$ (35.69)

equal to 1day. Such a time scale corresponds to the time needed to damp out a wave with wavenumber $k=\pi/\Delta a$, which is the highest wavenumber available to the grid.

Note that for most model grids, the huge maximum slope of 75/100 is indistinguishable from the more modest 1/100 value commonly used with redi_diffusion or gent_mcwilliams. The reason is that the grid aspect ratio is typically on the order of 1/1000 in the models. Therefore, the model grid essentially equates a slope $\ge 1/100$ to an infinite slope. Hence, it makes little difference whether one uses the same maximum slope for the biharmonic_rm option as for the redi_diffusion or gent_mcwilliams options, or if one allows the maximum slope for the biharmonic_rm option to be larger. For coding simplicity and consistency between the different schemes, MOM uses the same maximum slope slmx for all of the isoneutral mixing schemes.