In this section, we compare the analytic vorticity budget with the equivalent discretized vorticity budget for a Cgrid shallow water (SW) model and explain why the two budgets do not match. We then give results for the discretized vorticity budget and discuss the implications in terms of modelling of wind driven gyres in presence of steplike coastlines.
We consider the shallow water equations
(4.4) 
(4.5) 
Eqns. (4.14.2) or (4.24.3) can be
discretized in different ways. To simplify the discussion, we leave the
time derivative being continuous, and restrict ourselves to the Cgrid.
A useful general form of the SW equation is the following:
(4.9) 
(4.11) 
(4.12) 
(4.13) 
It is always possible to approximate the vorticity budget at the model boundary by using offcentered derivatives and interpolating some of the variables to the boundary. The model numerics, however, make no use of variable values found by such an interpolation and therefore, a vorticity budget calculated in this way must be considered distinct from the model vorticity budget. Such a budget might misrepresent the contribution of the different terms of the discretized equations of the model, especially if the error introduced by the coastline discretization is of lower order than are the truncation errors of the model. For this reason, we prefer to use the model vorticity budget. We note also that the truncation errors in the model vorticity budget are larger than the truncation errors in solving the shallow water equations, since vorticity is a higher order variable.
Figure 4.2 compares the rotated and nonrotated basin cases. The integrand ( i.e. the local ) is plotted as a function of distance around the basin perimeter and the position of the steps is evident from the abrupt jumps in the integrand value. When summed along the perimeter, is nonzero and is larger for the rotated basin experiment compared to the nonrotated basin experiment. Starting from this observation, we are interested in quantifying the importance of the steps over the global vorticity budget. First, as increased resolution leads to more steps, and due to the singular behavior of close to steps, it is no longer obvious that converges to zero with increasing resolution. From this point of view, is probably very sensitive to the formulation of the advective terms in (4.10), and extension in (4.1,4.3). Second, we want to investigate whether the overall circulation is sensitive to the presence of an extra term in the global vorticity budget, as can be a source or a sink term, depending on its sign.
We are interested in applying different formulations for the advectionCoriolis terms since we noted that was generally nonzero for the single gyre Munk problem, with the integrand being particularly large at steps. The two advective numerical schemes that we consider are the conventional formulation (based on Eq. 4.1) and the potential enstrophy conserving formulation (based on Eq. 4.3).
In addition to testing for sensitivity to the choice of advective schemes, we also consider different formulations of the stress tensor. That the overall circulation is primarily sensitive to the formulation of the stress tensor is the main result of AM, who found that the  formulation gave better results than the conventional formulation. We refer to [Gent, 1993] and [Shchepetkin and O'Brien, 1996] for a more complete discussion on appropriate viscous stress tensor formulations for the shallow water equations and we limit ourselves to the two stress tensor formulations used by AM. Below, we review these two formulations. Thus, we are interested in testing four combinations of two advective and two diffusive formulations. Table 4.1 summarizes these four different combinations.
With respect to the advectionCoriolis terms, we compare the conventional formulation to the potential enstrophy conserving formulation of [Sadourny, 1975] . For the conventional formulation, $\Cor$ and $\Phi$ are given by
4.15 
and for the potential enstrophy conserving formulation, $\Cor$ and $\Phi$ are given by
4.16 
Both formulations ensure a second order accuracy to the discretized SW equations. For the conventional formulation, changes are made to incorporate the boundary conditions at second order of accuracy, by using offcentered differencings close to the boundary. No boundary condition for the vorticity is required. However, since the enstrophy conserving scheme explicitly uses the vorticity, this formulation requires that vorticity be specified at boundary points. We choose to set the relative vorticity to zero along the model boundary, which is consistent with the freeslip boundary condition along straight walls. Also, contrary to the conventional formulation, no offcentered differencing is needed at the boundary for the computation of $\Cor$.
The two numerical formulations that we consider for the viscous terms are the divergencevorticity tensor formulation of [Madec et al, 1991] and the conventional fivepoint Laplacian. For the latter,
4.17 
As with the conventional advection formulation, changes are made here to incorporate the boundary conditions at second order accuracy, by using offcentered differencings. Another technical remark concerns the treatment of velocity points close to tips of land. For those points (the $u$ and $v$ points of Fig.4.3, the tip of the land is half a grid cell away. Let us focus on the $u$point. The problem is to evaluate the Laplacian of $u$ at this point. A fivepoint Laplacian requires knowledge of $\partial u/\partial y$ in the center of the northern and southern sides of the cell surrounding the $u$point, and $\partial u/\partial x$ on the eastern and western sides. The problem lies with $\partial u/\partial y$ on the northern side. The usual treatment would have
4.18 
which simplifies to
4.19 
because of the impermeability condition which sets $u_{i,j+1}$ to zero. Alternatively, one might take impermeability to imply that the tip is a stagnation point, in which case an offcentered differencing leads to
4.20 
A third logical possibility would be to apply the free slip condition at the tip to conclude that
4.21 
We choose the latter (4.21), in order to let the ``fluid'' slip as much as possible along the walls since the first two conditions (4.19,4.20) tend to slow down the boundary currents. A more accurate formulation of the boundary condition close to the steps can be derived using a finite volume formulation, which treats the northern viscous flux as a mean between (4.19) and (4.21). However, this would slow down the boundary current due to the use of (4.19). In addition, more accurate treatment of the steps have limited value as the steps are artificial.
The divergencevorticity (  ) form of the stress tensor leads to the following form for the Laplacians
4.22 
where $\delta$ is the divergence expressed at the $h$location (center of the cell). This formulation is more general in the sense that there is no adjustment of the formulation at the boundary. Another remark concerns the case of straight walls. In that particular case, there is no difference between the  stress tensor formulation and the traditional formulation. The difference is in the treatment of steps.
To illustrate this, we consider a comparison between the two stress tensor formulations for a forward step along a northflowing western boundary current. Choose $(i,j)$ so that, in Figure~\ref{cgrid3}, the $\zeta$point right at the tip of the land corner would have $(i,j+1)$ indices (see Fig.~A.1 for indices arrangement). Thus, the viscous terms under the  formulation are
4.23 
A final remark is that the divergence part of the viscous forces cancels out in the vorticity equation. Therefore, in the discretized vorticity equation, the  formulation leads to a viscous term that takes the form of the fivepoint Laplacian of the vorticity. This is not true of the conventional formulation.

By studying , we want to address several issues related to the accuracy of the different combinations of the advection and diffusion formulations and their influence on the strength of the overall circulation. Firstly, a major requirement is that, whatever the geometry of the basin, should converge to zero as the resolution goes to infinity. This test allows us to rank the performances of the model for the different combinations of advective and diffusive schemes. Of particular interest will be the importance of the advective formulation. A second concern is to assess whether the size of the artificial source or sink of vorticity due to influences the overall strength of the gyres. A third concern relates to the general accuracy of model vorticity budgets.
To address these issues, we make use of the conceptual experiment proposed by AM, in which a single gyre Munk circulation is computed in rotated and nonrotated square basins. In both cases, all parameters and forcing are unchanged except for the discretized coastline. The four combinations (A, B, C, D) of numerical formulations we propose to test are detailed in Table 4.1. One remark concerns the nonrotated basin results. There, since the conventional and  stress tensor formulations are identical, the results for the B combination are identical to the results for A. The same applies for the C and D cases.
We reproduce the results of AM in Figure 4.4. This figure shows the elevation fields for the A and B cases and for no rotation and a small rotation angle of 3.4$^o$. Clearly, the A case shows circulation patterns collapsing as the number of steps along the walls increases whereas, for the B case, the circulation is quite similar to the original nonrotated circulation. The results for C are not shown but are very similar to the results for A. The results for D show a small increase in the strength of the gyre compared to A, but the original overall circulation of AB with no rotation is not recovered (not shown).
Figure 4.5a shows the kinetic energy as a function of resolution for the various combinations and for a rotation angle of $3.4^o$. Only the B combination converges to nonrotated solutions. The A and C results are almost identical, but appear to converge to a kinetic energy that is reduced by over a factor of 2 compared to the nonrotated cases. For the D combination, kinetic energy decreases and then tends to slightly increase with increasing resolution and is generally much lower than for AB with no rotation or B with rotation.
As mentioned, the first consistency test related to the vorticity budget is to verify that converges to zero with increasing resolution. For all rotation angles considered and for the B combination, this statement appears to be true. For the other combinations (A, C, D), such is not the case, at least for certain angles. For instance, tends to increase or stay constant for the A, C and D combinations at $3.4^o$ (Fig. 4.5b). For the D case, increases dramatically with increasing resolution so much that becomes larger than the wind input. Associated with this is a reverse (negative) viscous flux. This behavior may have consequences on the stability of the model. Although no obvious numerical instabilities occurred for a rotation angle of $3.4^o$, numerical instabilities cause the model to crash for other angles, for example at $30^o$. It seems plausible that this behaviour is associated with the large (and opposing) advective and diffusive fluxes of vorticity near the model perimeter. In any event, it seems reasonable to conclude that the D combination is inappropriate. This implies that the  viscous formulation performs well only when used is conjunction with the enstrophy conserving advection. This finding complements that of AM. For the A and C combinations (Fig. 4.5b), does not converge toward zero with increasing resolution. Hence, these two combinations seem inappropriate, even if the resulting solutions are always stable.
We now address the issue of possible correlation between and the kinetic energy. Given that inertial runaway (the inability of simple models of the ocean to converge to a reasonable statistical mean solution as the eddy viscosity is decreased to the real value of the viscosity found in water) appears to be related to ``difficulties'' in balancing the global vorticity budget [Pedlosky, 1996] , it seems reasonable to ask whether the sign of is correlated with an indicator of the overall strength of the gyre, such as total kinetic energy. For example, when is negative, it adds to the wind input of vorticity and one might expect a stronger gyre to result. Some evidence that this may be the case is found by comparing B and D, which share the same formulation of the viscous terms. Figure 4.5b shows that is positive and larger for D than is the case for combination B. Thus the total wind plus advective input of vorticity to the basin is stronger in case B. As might have been anticipated, B shows a more energetic circulation (Fig.~\ref{ene4}a). It is also interesting to see whether there is any correlation between kinetic energy and the sign/strength of for a given formulation of the numerics. We restrict this discussion to the use of the B combination. From figures 4.6 and 4.7, which show the kinetic energy and advective/wind vorticity input ratio for a range of resolution and rotation angles, there does not appear to be any striking correlation. In particular, if we focus on the region of negative values of (i.e., for a case where has the same sign as the wind input), the kinetic energy for this region is not larger than the kinetic energy at the same resolution but for an opposite angle (in fact, the kinetic is slightly lower). Presumably the added advective flux in this region is locally balanced by the viscous terms, so that processes analogous to those thought to be responsible for inertial runaway do not lead to an increase in the overall strength of the gyre.
To conclude this section, we investigate the general accuracy of model vorticity budgets with respect to using the B combination, only, since this combination is the only one showing a robust convergence to zero with increasing resolution. As should ideally not be present in the vorticity budget, the viscous flux, $\fov$, can be either underestimated or overestimated (which modifies the local balance at the wall and therefore the strength of the gyre) and can be viewed as an error. From Figure 4.8 and for the range of resolution we used, varies between 5\% and 50\% of the wind input. The order of the convergence for with increasing resolution is fairly close to unity or slightly lower for all positive angles. For negative angles, we did not compute the convergence order because goes through a minimum (Figure 4.7) and had not asymptoted to an uniform convergence order at the highest resolutions we considered. A noteworthy point is that the effect of increasing the rotation angle (introducing more steps) seems to decrease the convergence order (1/2 at $20^{o}$). Paradoxically, however, the convergence order increases again to reach unity for $45^{o}$, the rotation angle at which the number of steps is maximum. In fact, at this angle even shows a negative offset compared to the $0^{o}$ angle.
Except for effects related to steplike boundaries, that the convergence order is unity follows directly from the order of discretization of the vorticity. Since the vorticity is one order higher a variable than is velocity, and since the velocity is computed at second order accuracy, it follows that the vorticity is at best accurate to first order. Therefore, can be considered an explicit first order (at best) error in the vorticity budget. For the B combination, we observe that the convergence order for varies between 1/2 and unity, depending on the rotation angle. In the 1/2 order case, errors (or discrepancies) vary between 25% (high resolution) and 50% (coarse resolution) and, in the first order case, they vary between 6% and 22%. The errors are much larger for the other combinations and can reach 100%. This implies that the accuracy of computing vorticity budgets from primitive equations models is fairly low, especially in absence of attention to the numerics. These errors may also vary a lot considerably with the discretized domain geometry.
Get the PS or PDF version here