35.1.8.8 Discretization details for the RM98 operator

Discretization of the RM98 skew-fluxes directly parallels that of the GM skew-fluxes discussed in Section 34.1.6.1. Just as for the Laplacian GM scheme, the skew-flux approach using the Griff ies et al. (1998) algorithm guarantees that the biharmonic scheme will not alter the tracer variance. Additionally, the biharmonic tracer fluxes are much less noisy than the fluxes resulting from computing the eddy-velocity (see discussion in Section 34.1.6.1). However, when applied to the biharmonic scheme, the Griff ies et al. (1998) algorithm brings a relatively heavy computational load. The algorithm described by RM98 is likely cheaper, but it computes eddy-velocities and so has not been implemented at this time.

The main difference between the Laplacian GM scheme and the biharmonic scheme is the extension of the stencil due to the Laplacian operator acting on the slopes. In order to maintain numerical stability in steep sloped regions, each slope comprising the elements of the Laplacian is weighted by its corresponding diffusivity. In the case of constant diffusivity and small slopes, this discretization reduces to that prescribed in the continuum. All flux components are integrated explicitly in time, since there is no diagonal contribution to the vertical tracer flux component.

To get started, recall that the off-diagonal term in the discretized small angle x-component of the Redi diffusive flux is given by

$${-F^{off-diag}_{i,k,j} = {1 \over 4 } \sum_{kr=0}^{1} \Delta^{(i,k-1+kr,j)}_{(i,k,j)} }$$

$$\sum_{ip=0}^{1} Ax^{(i+ip,k)}_{(i,k \vert i+ip,k-1+kr)} \; Sx^{(i+ip,k)}_{(i,k \vert i+ip,k-1+kr)} \delta_{z}T_{i+ip,k-1+kr},$$ (35.70)

where the terms in this equation are defined in Section 34.1.5.1. This form for the flux is taken for the discretization of the RM98 operator, with the diffusivity times the slope becoming the horizontal Laplacian of the biharmonic diffusivity times the slope

$${-F^{x-biharm}_{i,k,j} = {1 \over 4 } \sum_{kr=0}^{1} \Delta^{(i,k-1+kr,j)}_{(i,k,j)} }$$

$$\sum_{ip=0}^{1} \; \nabla^{2}_{h} \left( B x^{(i+ip,k)}_{(i,k \ve... ...} \, Sx^{(i+ip,k)}_{(i,k \vert i+ip,k-1+kr)} \right) \delta_{z}T_{i+ip,k-1+kr}.$$ (35.71)

Again, the inclusion of B x^(i+ip,k)_{(i,k | i+ip,k-1+kr)} inside the Laplacian is necessary for maintaining numerical stability in regions of steep slopes. Otherwise, for example, $\nabla^{2}_{h} Sx^{(i+ip,k)}_{(i,k \vert i+ip,k-1+kr)}$ could be large due to the presence of a steep slope, say Sx^(i+1+ip,k)_{(i+1,k | i+1+ip,k-1+kr)}. With B x^(i+ip,k)_{(i,k | i+ip,k-1+kr)} outside the Laplacian, the product $B x^{(i+ip,k)}_{(i,k \vert i+ip,k-1+kr)}\nabla^{2}_{h} Sx^{(i+ip,k)}_{(i,k \vert i+ip,k-1+kr)}$ will then be larger than allowed by numerical stability.

The meridional flux is given likewise by

$${-F^{y-biharm}_{i,k,j} = {1 \over 4 \cos\phi^{T}_{j} } \sum_{kr=0}^{1} \Delta^{(i,k-1+kr,j)}_{(i,k,j)} }$$

$$\sum_{jq=0}^{1} \nabla_{h}^{2} \left( B y^{(k,j+jq)}_{(k,j \vert ... ..., Sy^{(k,j+jq)}_{(k,j \vert k-1+kr,j+jq)} \right) \; \delta_{z}T_{k-1+kr,j+jq}.$$ (35.72)

The vertical flux is given by

$${-F^{z-biharm}_{i,k,j} = -{1 \over 4 dxt_{i} \; dht_{i,k,j} } \; \sum_{ip=0}^{1} dxu_{i-1+ip} }$$

$$\sum_{kr=0}^{1} \Delta^{(i-1+ip,k,j)}_{(i-1+ip,k+kr,j)}$$

$$\times \; \nabla_{h}^{2} \left( B x^{(i,k+kr)}_{(i-1+ip,k+kr \ver... ...)} Sx^{(i,k+kr)}_{(i-1+ip,k+kr \vert i,k)} \right) \delta_{x}T_{i-1+ip,k+kr} \;$$

$${1 \over 4 \cos\phi^{T}_{j} dyt_{j} \; dht_{i,k,j} } \; \sum_{jq=... ...}_{j-1+jq} dyu_{j-1+jq} \sum_{kr=0}^{1} \Delta^{(i,k,j-1+jq)}_{(i,k+kr,j-1+jq)}$$ -

$$\times \; \nabla_{h}^{2} \left( B y^{(k+kr,j)}_{(k+kr,j-1+jq \ver... ...} Sy^{(k+kr,j)}_{(k+kr,j-1+jq \vert k,j)} \right) \; \delta_{y}T_{k+kr,j-1+jq}.$$ (35.73)

In each of the above expressions, the horizontal Laplacian operator is given by

$$\nabla^{2}_{h} \alpha \frac{1}{dxt_{i} \, (\cos\phi^{T}_{j})^{2}} \left( \frac{ \alpha_... ...} - \alpha_{i}}{dxu_{i}} - \frac{ \alpha_{i} - \alpha_{i-1}}{dxu_{i-1}} \right)$$ =

$$\frac{1}{dyt_{j} \, \cos\phi^{T}_{j}} \left( \frac{ \cos\phi^{U}_... ...}} - \frac{\cos\phi^{U}_{j-1} (\alpha_{j} - \alpha_{j-1})} {dyu_{j-1}} \right).$$ + (35.74)