\documentstyle [11pt]{article} \oddsidemargin 0.5in \evensidemargin 0.5in \topmargin -0.5in \headsep .5in \textheight 8.8in \textwidth 6in \parskip 7pt \def\baselinestretch {1.7} \begin{document} \vspace* {3cm} \begin{center} \bf A wind-driven isopycnic coordinate model of the North and Equatorial Atlantic Ocean. I: Model development and supporting experiments \rm \vspace {3cm} by \vspace {1cm} Rainer Bleck and Linda T. Smith \vspace {3cm} Rosenstiel School of Marine and Atmospheric Science \\ University of Miami \\ 4600 Rickenbacker Causeway \\ Miami, Florida 33149 \vspace {1cm} \newpage \bf Abstract \rm \end{center} Numerical approximations to the dynamic equations are given which allow basin-size ocean circulation models formulated in isopycnic coordinates to accomodate variable bottom topography and irregular coastlines. Emphasis is placed on computational economy through the use of a split-explicit time integration scheme, on the proper formulation of the advection and Coriolis terms in the momentum equation in case of strongly varying layer thickness, and on the correct estimation of the horizontal pressure gradient force in grid boxes truncated by steep bottom slopes. The algorithms are tested in a series of two- and three-layer double-gyre experiments. In cases of steady forcing leading to a steady circulation we are able to reproduce the motionless final state in coordinate layers that are below the direct influence of the wind forcing. This includes layers intersecting topographic obstacles. A long-term (25-- to 30--$yr$) vacillation tentatively associated with the outcropping of isopycnals along the edge of the cyclonic gyre in a steadily forced model is documented. The paper forms the basis for a subsequent one in which circulation features obtained in a realistic North Atlantic setting will be discussed. \newpage \bf 1. Introductory Remarks \rm Numerical modeling of geophysical fluid systems using electronic computers began approximately four decades ago. Until the middle to late 1960s, ocean model development lagged significantly behind that of atmospheric models. A number of reasons (apart from societal needs) can be cited for this lag, among them the inherently greater complexity of circulation systems confined to closed basins, a highly nonlinear equation of state for seawater, and the lack of three-dimensional synoptic observations for ocean model initialization and verification. Furthermore, the computing power required to resolve the relevant hydrodynamic instability processes in grid-point space is far greater in the ocean than in the atmosphere because of the much smaller radius of deformation. The above factors, taken together, imply that a major commitment of institutional resources was required in the past to mount a credible oceanic modeling effort. Thus oceanic modeling has not become a ``cottage industry" like atmospheric modeling. Fiscal constraints often prohibited the luxury of a ``duplication of efforts", and thus steered the oceanographic community away from model diversity. Only now, with supercomputers being within easy reach of university researchers, can such a diversity be achieved and should, in fact, be encouraged. Our present modeling effort should be viewed in this context. This is the first of two papers describing our efforts to simulate the wind-driven circulation in the North and Equatorial Atlantic using a multilayer, isopycnic coordinate model. The present paper deals largely with numerical problems encountered while incorporating variable bottom topography into an isopycnic coordinate model. We also discuss in this paper certain long-term ($\sim$30 years) fluctuations observed in the strength of wind-forced gyres in coarse-mesh (i.e., noneddy resolving) model simulations. The subject of the second paper will be the actual results obtained from the Atlantic model. In isopycnic coordinate representation, the ocean is viewed as a stack of immiscible layers, each of which is characterized by a constant value of density and is governed by dynamic equations resembling the shallow-water equations. The layers interact through hydrostatically transmitted pressure forces. The analytic models of {\em Welander} [1966] and {\em Parsons} [1969], in which the vertical structure of the ocean is represented in terms of layers of constant density, can be considered as prototypes of the model employed in this study. {\em Huang} [1986] has recently extended the work of Welander and Parsons in numerical models with idealized forcing and basin geometry. The ventilated thermocline model of {\em Luyten et al.} [1983], which seeks to describe the vertical structure of a subtropical gyre, also represents that structure by means of constant-density layers. In contrast to multilevel models of stratified fluids formulated in conventional Cartesian coordinates, a layered (Lagrangian) model employing density as a vertical coordinate is capable of reducing the vertical truncation error by concentrating coordinate surfaces in regions characterized by large vertical and horizontal density gradients. In addition, the isopycnic framework conforms to the accepted view that large-scale oceanic mixing processes take place primarily along surfaces of constant potential density [{\em Montgomery}, 1938]. In the particular case of a purely wind-driven calculation, one important advantage of an isopycnic coordinate model lies in its ability to retain the buoyancy contrast between water masses regardless of the length of time for which the equations are integrated. These advantages are offset by numerical difficulties encountered where interior density surfaces come into contact with either the ocean surface or the ocean bottom, both of which are treated as coordinate surfaces. The quasi-isopycnic model of {\em Bleck and Boudra} [1981] avoids such coordinate intersections by locally reverting to a Cartesian coordinate representation wherever minimum isopycnic grid-point separation in the vertical cannot be maintained. By contrast, the pure isopycnic model described by {\em Bleck and Boudra} [1986] depicts outcrop lines as boundaries between regions of zero and nonzero isopycnic layer thickness and solves the mass continuity equation in such a way that the ``drying up" of a layer (tantamount to zero vertical grid-point separation) does not cause numerical instability. One particular algorithm which we have found to be well suited to this purpose is the flux corrected transport (FCT) algorithm [{\em Zalesak}, 1979]. We have chosen to build our present model upon the ``pure" isopycnic framework, primarily because we judge the inconvenience of employing the FCT algorithm in the mass continuity equation to be minor compared with the inconvenience of carrying a density advection equation, which the quasi-isopycnic model requires. As outlined by {\em Bleck and Boudra} [1986], the chief disadvantage of the pure isopycnic approach is that vanishing layer thickness prohibits the use of quadratic-conservative finite-difference operators. We will propose one possible solution to this problem. Some of the stability measures to be described below may appear rather arbitrary to readers unfamiliar with quasi-Lagrangian vertical coordinate modeling; however, in designing these measures, we have been careful not to jeopardize numerical consistency, defined here as the property of a finite-difference equation to approach the correct differential equation as the mesh size in space and time goes to zero. Until now, we have not attempted to incorporate separate temperature and salinity advection equations and an appropriate equation of state into the North Atlantic version of our isopycnic coordinate model. (These constraints are being removed in a parallel modeling effort; see {\em Bleck et al.}, [1989].) We do not see any conceptual problems with incorporating more realistic thermodynamics into the model, however. All diabatic processes can ultimately be reduced to prescriptions of the vertical velocity component $ d \alpha /dt $ ($ \alpha $ being specific volume), which can easily be added to the present set of equations. On the other hand, we do not claim that there will be no numerical difficulties associated with sharply rising isopycnals in the vicinity of deep convection zones. While the present model is ``rigorously" isopycnic, diabatic processes of the type just mentioned may be treated more appropriately in the quasi-isopycnic, hybrid coordinate framework employed by {\em Bleck and Boudra} [1981]. A numerical model is worthless if it is fraught with stability problems or is computationally inefficient. Most of our development efforts naturally were aimed at overcoming these difficulties. Another, somewhat more subtle design criterion which we chose to adopt is that the model should reproduce the total spindown of coordinate layers inaccessible to direct wind forcing in situations where both the external forcing and the ensuing circulation are steady [{\em Welander}, 1966; {\em Anderson and Killworth}, 1977]. We were able to reach an acceptable level of numerical efficiency through the use of a split-explicit time integration scheme. We also found a way to reliably diagnose the horizontal pressure force in the immediate vicinity of steep bottom slopes, a problem that has received considerable attention in the meteorological literature [e.g., {\em Mesinger}, 1982]. However, combining these two schemes in a way which does not interfere with the spindown process proved to be a more difficult task than we had anticipated. The paper is organized as follows. In section 2 we will discuss the model equations, including the split-explicit time differencing scheme. In section 3 we will present results of double-gyre spinup-spindown experiments conducted with a two- and three-layer version of the model set in a rectangular ocean basin. Emphasis will be on the effect of bottom topography on the spindown in situations where the topography rises above the lowest coordinate layer and thereby creates zones of zero layer thickness. We will also discuss a long-term vacillation found in the three-layer model. Section 4 will contain a summary. Technical details of the pressure gradient formula, quadratic conservation properties of the finite-difference equations, and the implementation of lateral stresses along topographic obstacles will be discussed in Appendices A--C. \vspace {1cm} \bf 2. Model equations \rm The layer-integrated nonlinear momentum equations in $ x,y,\alpha $ coordinates ($ \alpha $ being the specific volume) under adiabatic flow conditions are \begin{equation} \frac{\partial u}{\partial t} + \frac{\partial}{\partial x} \frac{{\bf v}^2}{2} - (\zeta + f)v = - \frac{\partial M}{\partial x} + \frac{1}{\Delta p} [g \Delta \tau_x + \nabla \cdot (\nu \Delta p \: \nabla u)] \label{eq:A-1a} \end{equation} \begin{equation} \frac{\partial v}{\partial t} + \frac{\partial}{\partial y} \frac{{\bf v}^2}{2} + (\zeta + f)u = - \frac{\partial M}{\partial y} + \frac{1}{\Delta p} [g \Delta \tau_y + \nabla \cdot (\nu \Delta p \: \nabla v)] \label{eq:A-1b} \end{equation} where $ \zeta = (\partial v / \partial x)_{\alpha} - (\partial u / \partial y)_{\alpha} $ is the relative vorticity, $ M \equiv gz + p \alpha $ is the Bernoulli function or Montgomery potential, $ \Delta p $ is the pressure thickness of the isopycnic layer, $ (\Delta \tau_x, \Delta \tau_y) $ are the wind and bottom drag-related stress differences in $ x $ and $ y $ direction between the top and bottom of the coordinate layer, and $ \nu $ is an eddy viscosity coefficient. All horizontal derivatives are evaluated at constant $ \alpha $, and the local time derivative applies to a fixed point in $ x,y,\alpha $ (as opposed to $ x,y,z $) space. Distances in $ x $ and $ y $ direction are measured in the projection onto a horizontal plane. No metric terms related to the slope of the $ \alpha $ surface need to be neglected in deriving (\ref{eq:A-1a}) and (\ref{eq:A-1b}). Under adiabatic flow conditions the thermodynamic equation is satisfied by setting the vertical velocity component $ d \alpha /dt $ to zero. The nonlinear continuity equation, written in terms of a generalized material vertical coordinate $ \; s \; $ [{\em Bleck}, 1978], is \begin{equation} \frac{\partial}{\partial t_{s}} \left(\frac{\partial p}{\partial s} \right) + \frac{\partial}{\partial x_{s}} \left(u \frac{\partial p}{\partial s} \right) + \frac{\partial}{\partial y_{s}} \left(v \frac{\partial p}{\partial s} \right) = 0 \label{eq:A-2} \end{equation} where $ \partial / \partial t_{s} $ , $ \partial / \partial x_{s} $, and $ \partial / \partial y_{s} $ denote differentiations with $ s $ held constant. (The usefulness of introducing a generalized coordinate will become apparent shortly.) Since barotropic gravity waves typically propagate through an incompressible ocean at least 20 times faster than other signals (the incompressibility condition serves to filter out sound waves), numerical efficiency can be greatly increased by removing these waves from the time-consuming layer-by-layer integration of the prognostic equations. This can be done in a number of ways: by keeping the lowest model layer at rest, by topping the ocean with a ``rigid lid", or by separating barotropic gravity waves from the slower wave modes supported by the stratified medium and advancing them in time in a highly simplified single-layer model. The rigid lid approximation [{\em Bryan}, 1969] requires solution of a Poisson equation after each ``baroclinic" time step. In the special case of a rectangular basin with bottom topography varying in only one coordinate direction, this equation can be solved efficiently by a direct (i.e., noniterative) method based on the fast Fourier transform (FFT). Iterative methods based on successive overrelaxation (SOR) perform under less restrictive conditions but require considerable insight into the precision requirements of the solution. We have found that in applications involving irregularly shaped basins and $x$- and $y$-dependent bottom topography, an explicit integration of the barotropic component of motion in split-explicit mode [{\em Gadd}, 1978] is more economical than a rigid lid scheme utilizing SOR, even if the latter is fully vectorized. Specifically, we determined that in the case of 1900 grid points (the number of points in our Atlantic model) an explicit time integration of the barotropic equations over a single ``baroclinic" time step requires as much time on a CRAY X-MP computer as 14 iterations in the vectorized SOR scheme. Since a realistic iteration count for a slowly evolving model state is about twice that number, the SOR scheme was found to be roughly twice as expensive as the explicit scheme. Our model tests have also shown that the baroclinic calculation consumes approximately the same amount of computer time per layer as does the barotropic calculation. Thus as the number of layers increases (to five, in the case of our Atlantic experiments), the percentage of total time consumed by the barotropic calculation decreases. While the experiments reported below would have allowed the use of an FFT-based rigid-lid scheme because they are based on idealized basin configurations, the fact that these are supporting experiments for the North Atlantic model made it seem advisable to utilize the split-explicit method exclusively. The split-explicit scheme requires that the prognostic variables be separated into their barotropic (i.e., depth-independent) and baroclinic components. The decomposition of the velocity components is straightforward. We write \[ u = u' + \overline{u} \hspace{2cm} \overline{u'} = 0 \]\[ v = v' + \overline{v} \hspace{2cm} \overline{v'} = 0 \] where the overbars indicate a mass-weighted vertical average. In decomposing the pressure field, we take note of the fact that layer thickness perturbations caused by barotropic flux divergences are proportional to the layer thickness itself. Rather than treating the pressure field as the sum of a barotropic and a baroclinic component, we therefore opt for a breakdown into multiplicative factors. We set \[ p = p'(1 + \eta) \] where the nondimensional variable $ \eta $, representing the barotropic component of the pressure field, is treated as depth-independent, while the baroclinic pressure field component $ p' $ has the property that its bottom value $ p'_b $ is time-independent. Since $ \eta \ll 1 $, we adopt the convention of approximating $(1 + \eta)$ by 1 wherever this expression appears as a factor. The barotropic and baroclinic velocity components then satisfy \[ \overline{u} = \frac{1}{p'_b} \int_{0}^{p'_b} u \: dp \hspace{2cm} \int_{0}^{p'_b} u' \: dp = 0 \]\[ \overline{v} = \frac{1}{p'_b} \int_{0}^{p'_b} v \: dp \hspace{2cm} \int_{0}^{p'_b} v' \: dp = 0 \] where the pressure at the ocean surface is set to zero. Equations for $ (\overline{u},\overline{v}) $ can be obtained, at least in principle, by converting (\ref{eq:A-1a}) and (\ref{eq:A-1b}) into layer-integrated flux form and summing up the resulting equations over all layers. Equations for $ (u',v') $ in each layer are then obtained by subtracting from (\ref{eq:A-1a}) and (\ref{eq:A-1b}) the barotropic equations so defined. We follow a more pragmatic approach. We write the equations for $ (\overline{u},\overline{v}) $ in the form \begin{equation} \frac{\partial \overline{u}}{\partial t} - f \overline{v} = - \alpha_{0} \frac{\partial p'_b \eta}{\partial x} + \frac{\partial \overline{u}^{*}}{\partial t} \label{eq:A-3a} \end{equation} \begin{equation} \frac{\partial \overline{v}}{\partial t} + f \overline{u} = - \alpha_{0} \frac{\partial p'_b \eta}{\partial y} + \frac{\partial \overline{v}^{*}}{\partial t} \label{eq:A-3b} \end{equation} where $ \partial (\overline{u}^{*}, \overline{v}^{*})/ \partial t $ comprise all those terms appearing in (\ref{eq:A-1a}) and (\ref{eq:A-1b}) that are not explicitly accounted for in (\ref{eq:A-3a}) and (\ref{eq:A-3b}), namely, the vertically integrated inertial, friction, surface stress, and bottom stress terms. Subtracting (\ref{eq:A-3a}), (\ref{eq:A-3b}) from (\ref{eq:A-1a}), (\ref{eq:A-1b}) respectively yields \begin{equation} \frac{\partial u'}{\partial t} + \frac{\partial}{\partial x} \frac{{\bf v}^2}{2} - (\zeta + f)v' - \zeta \overline{v} = - \frac{\partial M'}{\partial x} + \frac{1}{\Delta p'} [g \Delta \tau_x + \nabla \cdot (\nu \Delta p' \: \nabla u)] - \frac{\partial \overline{u}^{*}}{\partial t} \label{eq:A-4a} \end{equation} \begin{equation} \frac{\partial v'}{\partial t} + \frac{\partial}{\partial y} \frac{{\bf v}^2}{2} + (\zeta + f)u' + \zeta \overline{u} = - \frac{\partial M'}{\partial y} + \frac{1}{\Delta p'} [g \Delta \tau_y + \nabla \cdot (\nu \Delta p' \: \nabla v)] - \frac{\partial \overline{v}^{*}}{\partial t} \label{eq:A-4b} \end{equation} where $ M' = M - p'_b \eta \: \alpha_{0} $. The role of $ \partial (\overline{u}^{*}, \overline{v}^{*})/ \partial t $ in (\ref{eq:A-4a}) and (\ref{eq:A-4b}) is to assure that $ \int (u',v') dp $ remains zero despite the appearance of several terms in (\ref{eq:A-4a}) and (\ref{eq:A-4b}) which generally do not integrate to zero in the vertical. In other words, the vector $ \partial (\overline{u}^{*}, \overline{v}^{*})/ \partial t $ represents the rate at which the remaining terms in (\ref{eq:A-4a}) and (\ref{eq:A-4b}) acting on $ (u',v') $ generate a barotropic velocity component. We account for the presence of $ \partial (\overline{u}^{*}, \overline{v}^{*})/ \partial t $ in solving (\ref{eq:A-3a})--(\ref{eq:A-4b}) as follows. We start out at time $ t $ with a three-dimensional $ (u',v') $ field that satisfies $ \int (u',v') dp = 0 $, and solve (\ref{eq:A-4a}) and (\ref{eq:A-4b}) in all model layers for $ (u' + \overline{u}^{*}, v' + \overline{v}^{*}) $ at $ t + \Delta t $. The condition that $ (u',v') $ must integrate to zero at all times is then used to separate $ (u',v') $ from $ (\overline{u}^{*}, \overline{v}^{*}) $ at time $ t + \Delta t $. Having thus isolated the average value of $ \partial (\overline{u}^{*}, \overline{v}^{*}) / \partial t $ over the time interval $ \Delta t $, we feed this quantity, divided by $ \Delta t $, into (\ref{eq:A-3a}) and (\ref{eq:A-3b}) as a forcing term as we advance the barotropic flow over many small time steps from $ t $ to $ t + \Delta t $. We now turn our attention to the continuity equation (\ref{eq:A-2}). After splitting the dependent variables into their barotropic and baroclinic parts (note that $ \partial p / \partial s = (1 + \eta) \partial p' / \partial s $), and approximating $ (1 + \eta) $ by 1 as mentioned before, we obtain \begin{equation} \left(\frac{\partial p'}{\partial s} \right) \frac{\partial \eta}{\partial t} + \frac{\partial}{\partial t_{s}} \left( \frac{\partial p'}{\partial s} \right) + \nabla \cdot \left( \overline{\bf v} \frac{\partial p'}{\partial s} \right) + \nabla \cdot \left( {\bf v}' \frac{\partial p'}{\partial s} \right) = 0 \label{eq:A-5} \end{equation} By defining both the top and the bottom of the ocean as $ s $-surfaces we can interchange temporal and spatial differentiation on the one hand with vertical integration over the total water column on the other hand. Since $ p'_b $ is by definition time-independent and the vertical integral over $ u' $ and $ v' $ is zero, terms 2 and 4 in (\ref{eq:A-5}) vanish upon vertical integration. The remaining terms form the desired predictive equation for $ \eta $: \begin{equation} \frac{\partial p'_b \eta}{\partial t} + \nabla \cdot ( \overline{\bf v} p'_b) = 0 \label{eq:A-6} \end{equation} Tendency equations for the thickness $ \Delta p' $ of individual coordinate layers are obtained by integrating (\ref{eq:A-5}) over each layer separately. By substituting from (\ref{eq:A-6}), we obtain \begin{equation} \frac{\partial}{\partial t} \Delta p' + \nabla \cdot ( {\bf v}' \Delta p') = \frac{\Delta p'}{p'_b} \nabla \cdot ( \overline{\bf v} p'_b ) - \nabla \cdot ( \overline{\bf v} \Delta p') \label{eq:A-7} \end{equation} Note that the three flux divergence terms in (\ref{eq:A-7}) cancel exactly upon summation over all layers, as they should if $ p'_b $ is to remain time-independent. It follows that (\ref{eq:A-4a}), (\ref{eq:A-4b}), and (\ref{eq:A-7}) do not transmit barotropic gravity waves and can therefore be solved using a comparatively long, ``baroclinic" time step. One prognostic variable in the split-explicit model which is only indirectly available from the rigid-lid scheme is the sea surface elevation. It is given by $ M/g \; = \; (M' + p'_b \eta \: \alpha_{0})/g $ in layer 1. We will make use of this field when we discuss North Atlantic simulations in Part II of this paper. As mentioned in the introduction, (\ref{eq:A-7}) is solved in the model by the FCT method [{\em Zalesak}, 1979; {\em Bleck and Boudra}, 1986]. The purpose of this scheme is to modify the finite-difference approximations of $ {\bf v}' \Delta p' $ in (\ref{eq:A-7}) to prevent overshooting and undershooting (particularly the occurrence of negative values) in the $ \Delta p' $ field. Note that we apply the FCT scheme only to one of the three flux divergence terms in (\ref{eq:A-7}), $ \nabla \cdot ({\bf v}' \Delta p') $. The other two divergence terms appearing on the right-hand side of (\ref{eq:A-7}) act as small source terms; their effect on the $ \Delta p' $ field is taken into account separately. While these two terms are small and seemingly unimportant, they play a crucial role in maintaining the consistency of the split-explicit scheme. Since the flux-limiting process is carried out in each coordinate layer separately, there is no guarantee that the resulting mass flux divergences add up to zero in the vertical, which they must to conserve $ p'_b $. To prevent the FCT scheme from inadvertently changing $ p'_b $, the amounts removed by the flux-limiting algorithm are summed up in each grid column. The resulting values are treated as a barotropic correction to the three-dimensional mass flux field, their purpose being the restoration of the correct bottom pressure. By distributing these corrective fluxes uniformly over the ``upstream" grid column (upstream from the point of view of the corrective fluxes), they can be incorporated into the flux divergence calculation without endangering the positive-definiteness of the layer thickness field. The barotropic equations are time-integrated with the forward-backward scheme which permits a time step almost twice as large as the one used in leapfrog integrations. The scheme employs a forward (Euler) step in the linearized barotropic continuity equation (\ref{eq:A-6}), \[ p'_b \eta^{n+1} = p'_b \eta^{n} - \Delta t [\frac{\partial}{\partial x} (\overline{u} p'_b) + \frac{\partial}{\partial y} (\overline{v} p'_b) ]^{n} \] (the superscript $ n $ here indicates the time level), while the pressure force term in the momentum equations is subsequently evaluated in backward mode. The Coriolis term is treated in de facto leapfrog mode by switching the temporal order in which the two momentum equations are solved, and by always using the ``latest" available value of $ f \overline{u}, f \overline{v} $: \vspace{11 pt} n even \[ \overline{u}^{n+1} = \overline{u}^{n} + \Delta t [ - \alpha_{0} \frac{\partial p'_b \eta^{n+1}}{\partial x} + f \overline{v}^{n}] \]\[ \overline{v}^{n+1} = \overline{v}^{n} + \Delta t [ - \alpha_{0} \frac{\partial p'_b \eta^{n+1}}{\partial y} - f \overline{u}^{n+1}] \] \vspace{11 pt} n odd \[ \overline{v}^{n+1} = \overline{v}^{n} + \Delta t [ - \alpha_{0} \frac{\partial p'_b \eta^{n+1}}{\partial y} - f \overline{u}^{n}] \]\[ \overline{u}^{n+1} = \overline{u}^{n} + \Delta t [ - \alpha_{0} \frac{\partial p'_b \eta^{n+1}}{\partial x} + f \overline{v}^{n+1}] \] Since viscosity in (\ref{eq:A-4a}) and (\ref{eq:A-4b}) is weighted by layer thickness, the large spatial variability of $ \Delta p' $ may jeopardize linear numerical stability. To prevent this, we use a $ 10 \times 10^4 $ Pa cutoff value (equivalent to a 10-m layer thickness) for $ \Delta p' $ in the denominator of the viscous stress term in (\ref{eq:A-4a}) and (\ref{eq:A-4b}). Coordinate interface stresses arising from wind forcing and bottom drag are evaluated by assuming fixed values for both the top and the bottom boundary layer depth. (These values are presently chosen as $ 100 \times 10^4 $ and $ 10 \times 10^4 $ Pa, or 100 m and 10 m, respectively.) The stress at the sea surface is extraneously prescribed (see the following section), while the bottom stress is parameterized by the standard bulk method as \[ \mbox{\boldmath $\tau$} = - c_D \: \rho_0 \: |{\bf v}| \: {\bf v} \] where $ {\bf v} $ is the near-bottom velocity and $ c_D = 3 \times 10^{-3} $. The calculation of the interfacial stress differences $ (\Delta \tau_x, \Delta \tau_y) $ in (\ref{eq:A-1a}) and (\ref{eq:A-1b}) assumes a linear stress profile in both boundary layers and is carried out in such a way that the stress exerted on a zero-thickness layer located within a boundary layer defaults to $ \partial \mbox{\boldmath $\tau$} / \partial z $. This assures that the velocity field in massless layers is subjected to the same combination of forces as the velocity field in adjoining mass-containing layers. We neglect the momentum advection terms in the horizontal momentum equations (\ref{eq:A-1a}) and (\ref{eq:A-1b}), an acceptable omission for basin-scale experiments such as those reported here. However, the full form of the continuity equation (\ref{eq:A-2}) is used, as is typically done in isopycnic coordinate modeling where large horizontal variations in layer thickness preclude linearization of this equation. (The predominance of the role of this latter nonlinearity over that associated with the inertial terms, in the case of basin-scale studies, has been emphasized by {\em Huang} [1986, 1987].) Even though conservation of enstrophy or potential enstrophy is not strictly an issue if the nonlinear terms are removed from the momentum equations, we found that the model reacts unfavorably to a strongly varying bottom slope unless the Coriolis terms $ (-vf, uf) $ are formulated like the nonlinear potential-enstrophy conserving operators discussed in Appendix B. This formulation entails replacing the velocities $ (u,v) $ by the mass fluxes $ (u \Delta p, v \Delta p) $, and the Coriolis parameter $ f $ by the linear potential vorticity $ f / \Delta p $. If simpler formulations are used, the model produces unrealistic solutions or becomes unstable due to false energy production by the Coriolis terms. Simultaneous conservation of finite-difference potential enstrophy in individual layers and in the total column is achieved by formulating the linearized finite-difference analogs of the terms $ [- (\zeta + f)v' - \zeta \overline{v} $, $ (\zeta + f)u' + \zeta \overline{u}] $ in (\ref{eq:A-4a}) and ( \ref{eq:A-4b}) as explicit differences between the finite-difference analogs of $ [- (f/\Delta p) \: v \Delta p $, $ (f/\Delta p) \: u \Delta p] $ and $ [- (f/p'_b) \: \overline{v} p'_b $, $ (f/p'_b) \: \overline{u} p'_b] $. Use of the latter expressions in (\ref{eq:A-3a}) and (\ref{eq:A-3b}) guarantees that the sum of the barotropic and baroclinic equations has the proper finite-difference form. During our initial experiments with the actual topography of the North Atlantic, we discovered that further modifications to the finite-difference equations are needed to assure numerical stability. These modifications are discussed in Appendix C. Although the relatively simple bottom topography in the experiments discussed below may not call for these extra measures, we considered it nonetheless prudent to incorporate them. \vspace {1cm} \bf 3. Supporting experiments \rm In this section we will demonstrate the numerical properties of the isopycnic model in the simplified case of a rectangular ocean basin with an eastward-sloping bottom near the western basin boundary. We use the terms ``low-cut" and ``high-cut" topography to indicate the extent to which the slope extends up the western sidewall. In the high-cut case, the slope ascends from 4000 m to a depth of 200 m over a horizontal distance of 2000 km, whereas in the low-cut case it ascends to 1550 m over a distance of 1400 km. We will present results of coarse-resolution (i.e., noneddy resolving), double-gyre, wind-driven spinup experiments with two and three model layers. The rectangular ocean basin in these experiments is dimensioned 4400 km east-west by 6000 km north-south, with uniform grid spacing $\Delta x = \Delta y =$ 200 km. The vertical structure for the simplest model configuration consists of two layers whose initial thicknesses are 1000 and 3000 m (1000 and 3000 $\times\,10^4$ Pa, to be precise) in the flat-bottom region. The specific volume values in these layers are 0.974 and 0.973 $\times\,10^{-3}\mbox{m}^3\mbox{kg}^{-1}$, respectively, equivalent to potential density values 26.6 and 27.4. For the three-layer model configuration, the corresponding values are: layer thickness 300, 700, and 3000 m; specific volume 0.975, 0.974, and 0.973$\,\times\,10^{-3}\mbox{m}^3 \mbox{kg}^{-1}$; potential density 26.0, 27.0, and 27.7. The Coriolis parameter varies as $ f=f_0+(y-Y/2){\beta} $, where $ y $ is the north-south coordinate, $ Y $ is the north-south basin size, $ f_0= 0.83 \times 10^{-4}\mbox{s}^{-1} $, and $ \beta = 2 \times 10^{-11} \mbox{m}^{-1}\mbox{s}^{-1} $. The model is driven by an east-west surface wind stress pattern of the form \[ \tau_x =\tau_0 \: \cos \left[ 2 \pi \left( \frac{y}{Y} - \frac{1}{2} \right) \right] \] where $ \tau_0 / \rho_0 = 1.0 \times 10^{-4} \mbox{m}^2 \mbox{s}^{-2} $. The horizontal eddy viscosity coefficient $ \nu $ is set equal to $ 10^5 \mbox{m}^2 \mbox{s}^{-1} $. No-slip boundary conditions are imposed on all lateral boundaries; the implementation of these conditions along submerged sidewalls is discussed in Appendix C. The low-cut and high-cut bottom topography configurations are illustrated in Figure 1. {\em Two-layer model, low-cut topography}. The bottom topography in this simplest of our cases does not reach high enough to intersect the layer interface. Since the bottom layer is well beyond the reach of direct wind forcing, the end result to be expected in the absence of time-dependent forcing or internally generated hydrodynamic instabilities is a totally quiet bottom layer [{\em Anderson and Killworth}, 1977]. This is borne out in Figure 2 by the fact that the barotropic mass-flux stream function gives no indication of the presence of the sloping bottom. (We refrain from showing stream function plots in individual model layers because a plot of the layer-1 stream function is identical to Figure 2, while a plot of the layer-2 stream function is an empty frame.) {\em Two-layer model, high-cut topography}. In this case the bottom topography does intersect the layer interface. Westward of the intersection line, where the bottom layer is massless, the model is expected to behave like a single-layer model, i.e., the bottom topography there should distort the barotropic stream function so as to allow the flow to follow isolines of $ f / \Delta p $. (A discussion of the underlying principles can be found, for example, in {\em Holland} [1967].) Eastward of the intersection line, a bottom layer exists. As in the previous experiment, this layer should spin down completely and thus should insulate the stream function east of the intersection line from the effects of variable water depth. The resulting ``mixture" of one-and two-layer model behavior is shown in Figure 3. In order to give an impression of the time required for the complete spinup-spindown process, we show in Figures 4 and 5 time traces of layer kinetic energy for the two-layer, high-cut topography experiment. The difference between the two figures is in the pressure gradient formula used (see Appendix A). Figure 5 shows that the kinetic energy in the bottom layer is virtually zero after about 10 years of integration in the run based on pressure gradient scheme 2. Figure 4, on the other hand, indicates that a nonzero level of kinetic energy is being maintained against dissipative forces in the bottom layer if pressure gradient scheme 1 is used. The explanation for this behavior is as follows. The wedge formed by the bottom layer in the slope zone is too thin at some grid points to allow in situ calculation of the pressure gradient force. As outlined in Appendix A, scheme 1 then infers the pressure gradient from the layer above. Since the top layer is not motionless, it is subject to nonzero pressure forces which are hereby transmitted to the bottom layer. Scheme 2, which does not go outside the layer to estimate the pressure gradient in grid boxes truncated by bottom topography, avoids this problem. {\em Three-layer model, flat bottom}. By choosing a rather shallow top layer in the 3-layer runs, we allow Ekman suction to pull layer 2 to the surface in the cyclonic gyre. ({\em Huang} [1986] refers to this as the supercritical state.) The outcropping interface between layers 1 and 2 delineates a surface front skirting the cyclonic gyre. Even though the model is nominally noneddy resolving, the baroclinic structure thus created does appear to lead to frontal instabilities which cause a general unsteadiness in the kinetic energy trace (Figure 6) and allow energy to be transmitted to the bottom layer. (Note that a steady interface can transmit torques but not kinetic energy.) There is also a prominent long-term ($\sim$25 years) vacillation in the model energetics, both with respect to kinetic energy (shown) and potential energy (not shown). We are presently not equipped to analyze the connection, if one exists, between the relatively short-term energy fluctuations and the long-term vacillation visible in Figure 6. It is fairly straightforward, however, to link the 25-year cusps in the kinetic energy curve to a particular event in the evolution of the flow field. To illustrate the event, we show in Figure 7 snapshots of layer-1 thickness (plotted as depth of the layer-1/layer-2 interface) at three consecutive times a few years before, at the height of, and a few years after the energy peak near year 38 (see Figure 6). The item of interest is the isolated body of layer-1 water seen in Figure $7a$ near the northern basin boundary. This body can be traced back clockwise to the eastern sector of the cyclonic gyre from which it emerged several years earlier as an identifiable entity. Its source region is characterized by considerable long-term variability in the position of the surface front. The cusp in the kinetic energy curve occurs when the body of layer-1 water is being sucked into the western boundary current. Figure $7b$ shows that the body is hidden from view during that phase. Whether it emerges intact from the boundary current and continues its counterclockwise motion around the cyclonic gyre is difficult to assess. Figure $7c$ does show a perturbation in the layer thickness isopleths in the southern sector of the anticyclonic gyre at a time when energy returns to its normal level. However, a tracer advection experiment would have to be carried out to establish a clear connection between this rather diffuse thickness anomaly and the far more distinct pool of layer-1 water seen earlier along the northern boundary. Such an experiment could also be used to establish whether remnants of this disturbance, as it comes around full circle, trigger the formation of a new one in the eastern part of the cyclonic gyre, or simply a rejuvenation of the old one. Experiments of this type are planned. {\em Three-layer model, high-cut topography}. Time traces of layer kinetic energy for this experiment are shown in Figure 8. Once again, there is a marked 25-year cycle in the energy curves, except that the cusps mark minima, rather than maxima as seen in the flat-bottom case. Figure 6 also differs from Figure 8 in the energy level of the lowest layer. Since most of the kinetic energy in our noneddy-resolving model is found in the western boundary currents, we attribute the differences between Figures 6 and 8 to the small water depth at the western edge. As in the flat-bottom experiment, a cusp in the energy curve occurs every time a body of layer-1 water originating in the eastern part of the cyclonic gyre and migrating anticlockwise around it gets sucked into the southward flowing boundary current. This sequence of events is shown in Figure 9 which is rather similar to Figure 7. We do not pretend to understand why the seemingly undisturbed circulation state shown in Figures $7b$ and $9b$ is associated with an energy maximum in one case and with a minimum in the other. The effect of the small water depth in the high-cut experiment on the energy level in the bottom layer seems straightforward by comparison. More analysis will be required to fully explain these phenomena. As expected, the barotropic mass-flux stream functions for the 3-layer flat-bottom and high-cut experiments are similar to those shown in Figures 2 and 3, respectively, and will therefore not be shown here. There are slight differences in gyre strength which correspond to variations in the kinetic energy level seen in Figures 5, 6, and 8. We wish to conclude by pointing out that vacillations of the type discussed above also occur in the North Atlantic version of our model, where their period is close to 30 years. In fact, we found them there first, which led to the simplified 3-layer experiments described above. These vacillations, incidentally, are not sensitive to the presence of the nonlinear advection terms in the momentum equations. \vspace {1cm} \bf 4. Summary \rm This is the first of two papers dealing with the development of a model of the North and Equatorial Atlantic Ocean framed in isopycnic coordinates. The forthcoming second paper will deal with circulation features obtained by driving the model with an observed wind stress distribution in a basin with coastline geometry and bathymetry appropriate to an Atlantic simulation. The emphasis in the present paper is on numerical algorithms for solving the dynamic equations in isopycnic coordinates in the presence of steep bottom topography which is allowed to intersect layer interfaces. The model contains no thermodynamic equation, other than the statement that isopycnic layer interfaces are material surfaces. The remaining prognostic equations may be viewed as the shallow-water equations applied to a stack of individual fluid layers. The layers interact not through interface friction (although this can be added if desired), but through hydrostatically transmitted pressure torques. Due to computer time limitations, the model is not eddy resolving at present. A large part of the paper is taken up by the documentation of the split-explicit time integration scheme which frees the equations solved in individual layers from the severe time-step constraint posed by barotropic gravity waves. Another aspect which is thoroughly documented because of its apparent impact on the stability of simulations with strongly varying layer thickness is the formulation of the advection and Coriolis terms in the momentum equation. Finally, we give broad room to the problem of estimating the horizontal pressure force in grid boxes truncated by steep bottom slopes. Among the numerous trial experiments that paved the way for the North Atlantic model we selected a few that demonstrate (1) our ability to achieve near-complete cessation of motion in layers not directly influenced by the wind stress (particularly layers that intersect steep bottom slopes), and (2) the tendency toward long-term vacillations in model versions configured to show isopycnal outcropping along the fringe of the cyclonic gyre. These experiments go beyond the scope of previous multi-layer wind-driven model studies which either confine topographic variations to the lowest layer [e.g., {\em Thompson and Schmitz}, 1989] or assume a priori the existence of a steady state [{\em Huang}, 1986, 1987]. The generation mechanism for the instabilities associated with both short-term and long-term energy fluctuations in our nominally noneddy-resolving 3-layer simulations is not fully understood at this time and deserves further study. \newpage \begin{center} APPENDIX A \\ \bf Determination of the horizontal pressure force on steep bottom slopes \rm \end{center} In isopycnic coordinates the horizontal acceleration $ \alpha \nabla_{z} p $ is given by the isopycnic gradient of Montgomery potential $ \nabla_{\alpha} M $, where $ M = gz + \alpha p $. Problems arise if the acceleration in truncated grid boxes, like the shaded area in Figure 10, is to be computed. Since we assume the $k$th layer ($ k $ being the layer index) to extend in deflated or ``massless" state up the slope past the intersection point, we could, in principle, express $ (\partial M / \partial x)_{\alpha} $ at $ (x_2-x_1)/2 $ by $ [M_{k}(x_2) - M_{k}(x_1)] / (x_2 - x_1) $ with $ M_{k}(x_1) = gz_b(x_1) + \alpha_{k} p_b(x_1) $ (where the subscript $ b $ denotes bottom values). Though numerically consistent, this approach is not satisfactory because it results in a nonzero pressure force even if the fluid is at rest, a state characterized by horizontal isopycnals in the fluid interior. To reduce this error, we need to construct an ``underground" estimate of $ M $ by downward integration of the hydrostatic equation from the ocean bottom to an extrapolated interface pressure $ p_{k+1/2}^{(extr.)} (x_1) $. Since the layer variable $ M $ is piecewise constant in the vertical with discontinuities at layer interfaces, the information needed for such a downward integration includes the interface pressures of all underground surfaces down to $ k + \frac{1}{2} $. A number of schemes, employing either sideways or vertical extrapolation of the $ p $ field, can be formulated to accomplish this task. Vertical extrapolation schemes, which include virtually all methods ever invented for expressing $ \alpha \nabla_{z} p $ as the sum of $ \nabla_{s} \phi $ and $ \alpha \nabla_{s} p $ (where $ \nabla_{s} $ denotes a gradient evaluated on a model coordinate surface) are found to work reasonably well in atmospheric models because of the comparatively gentle slope of terrestrial mountains. However, we have not discovered a vertical extrapolation scheme that works well near steep oceanic bottom topography. We therefore limit our attention to two schemes which rely predominantly or even exclusively on sideways extrapolation. They are discussed below. \underline{Scheme 1}. This scheme infers the pressure gradient in ``marginal" (i.e., severely truncated) grid boxes from the $ M $ gradient in a slab of finite thickness which touches, but does not intersect, the bottom. A finite pressure thickness $ \epsilon $, as indicated by cross-hatching in Figure 10, must be assigned to this slab to eliminate temporal discontinuities in $ \nabla_{\alpha} M $ which otherwise would occur whenever a coordinate interface sweeps past a bottom grid point. A zero-order extrapolation of the $ M $ gradient downward from the cross-hatched area in Figure 10 implies that all coordinate surfaces in the column below are treated as isobaric. Viewed from this perspective, the above scheme belongs to the class that infers underground $ M $ values from a horizontally extrapolated $ p $ field. Some aspects of vertical extrapolation are retained, however. Consider a situation where the mass field in the model changes in such a way that the coordinate surface labeled $ m + \frac{1}{2} $ in Figure 10 approaches the ground. As long as $ p_b - p_{m+1/2} $ (subscript $ b $ denoting the bottom value) exceeds $ \epsilon $, $ \partial M / \partial x $ in the truncated grid boxes is determined from $ M $ values in layer $ m $. Once the interface $ m + \frac{1}{2} $ touches bottom at $ x_1 $, our scheme stipulates that the $ \partial M / \partial x $ value in layer $ m $ is set equal to $ \partial M / \partial x $ computed in layer $ m+1 $. From this time on, the pressure force in the $ m $-th layer has a value we would obtain by treating $ p_{m+1/2}(x_1) $ as equal to $ p_{m+1/2} (x_2) $.During the intervening stage, when $ 0 < p_b - p_{m+1/2} < \epsilon $, the gradual shift in the $ M $ gradient calculation from layer $ m $ to $ m+1 $ can be associated conceptually with a gradual migration of $ p_{m+1/2} $ from its original above-bottom value $ p_b - \epsilon $ to $ p_{m+1/2}(x_2) $. It is during this stage that Scheme 1 retains a vertical extrapolation aspect. In the interest of computational economy we only extend the vertical averaging operation over two layers, namely, the one in the truncated box under consideration and the layer above it. Specifically, we evaluate the pressure force in x direction in the k-th layer by \[ \frac{1}{\epsilon} \int_{max(p_b - \epsilon, p_{k+3/2})}^{p_b} \frac{\partial M}{\partial x} dp \] where $ p_b = min[p_b(x_1),p_b(x_2)] $. Note that we divide the integral by $ \epsilon $ regardless of the actual integration limits. This has the effect of reducing the pressure force to zero in grid boxes which are in, or are approaching, massless state. \underline{Scheme 2}. This scheme, which is based on a suggestion by A. Tafferner (personal communication,1988), relies exclusively on sideways extrapolation within each isopycnic layer. The underlying idea is to infer the $ M $ gradient in ``marginal" (severely truncated) grid boxes from the $ M $ gradient in neighboring grid boxes exhibiting sufficient layer thickness. Our implementation of Scheme 2 is as follows. Suppose we want to compute the pressure gradient force in $ x $ direction at a $ u $ velocity point $ i,j $. (Relative to this point, mass field variables have fractional indices in $ i $ direction.) The first step is to form a weighted average of the quantity $ \partial M / \partial x $ based on its value at the four nearest grid points $(i',j')$ = $(i\pm1,j)$, $(i,j\pm1)$, \[ \left( \frac{\partial M}{\partial x} \right)_{av} = \left[ \sum w_{i',j'} (\frac{\partial M}{\partial x})_{i',j'} \right] / \left[ \sum w_{i',j'} \right]. \] The four weights, defined as \[ w_{i\pm1,j} = \min\,(\epsilon, \Delta p_{i \pm 1/2,j}, \Delta p_{i \pm 3/2,j}) \]\[ w_{i,j\pm1} = \min\,(\epsilon, \Delta p_{i-1/2,j\pm1}, \Delta p_{i+1/2,j\pm1}) \] are designed to assure that only $ M $ gradients from grid boxes containing a minimal amount of mass are used in the average. $ \epsilon_1 $ is a small positive quantity added to the denominator to prevent division by zero. The intent is to substitute the horizontally interpolated quantity $ (\partial M / \partial x)_{av} $ for the locally computed value $ (\partial M / \partial x)_{ij} $ if the grid box $ i,j $ is severely truncated, i.e, if the minimum of the two layer thickness values $ \Delta p_{i \pm 1/2,j} $ is too close to zero to allow the pressure force to be given by $ (\partial M / \partial x)_{ij} $. In order to avoid temporal discontinuities, we do not force the algorithm to make a clearcut choice between $ (\partial M / \partial x)_{ij} $ and $ (\partial M / \partial x)_{av} $ at every point in space and time. Instead, we define the finally chosen pressure force value as a weighted average of the two alternate expressions, i.e., \[ \left( \frac{\partial M}{\partial x} \right)_{final} = \left( \frac{\partial M}{\partial x} \right)_{av} + \frac{\min\,(\epsilon, \Delta p_{i-1/2,j}, \Delta p_{i+1/2,j})}{\epsilon} \left[ \left( \frac{\partial M}{\partial x} \right)_{ij} - \left( \frac{\partial M}{\partial x} \right)_{av} \right] \] Note that the above scheme has the seeming defect of reverting to $ (\partial M / \partial x)_{ij} $ if the four weights $ w_{i\pm1,j},w_{i,j\pm1} $ and \[ \min\,(\Delta p_{i-1/2,j}, \Delta p_{i+1/2,j}) \] are all zero. However, the $ u $ point in question is embedded in a cluster of massless grid points in this case; accelerations computed there are therefore inconsequential. In this context we should mention one other measure which has been found to improve the spatial coherence of velocity data at massless and near-massless grid points. After the $ u $ and $ v $ field has been advanced over one time step, each grid point value $ u_k $ and $ v_k $ ($ k $ being the vertical grid index) is replaced by a vertical average taken over a 5-m interval centered at $ (p_{k+1/2} + p_{k-1/2})/2 $. Since $ u $ and $ v $ are assumed constant in the vertical within each layer, this averaging operation has no effect on the velocity field except where the layer thickness drops below 5 m. Massless velocity points in particular are hereby assigned velocity values found in nonmassless layers above or below. \newpage \begin{center} APPENDIX B \\ \bf Conservation properties of nonlinear terms \rm \end{center} As mentioned in the Introduction, an isopycnic coordinate model resembles a stack of shallow-water models. A question worth pursuing is whether finite difference techniques developed for conserving quadratic properties in shallow-water models can be incorporated into multi-layer models. Particularly relevant are the laws governing conservation of layer potential vorticity and layer potential enstrophy. {\em Sadourny} [1975] has shown that the nonlinear terms in the shallow-water momentum equations will conserve potential enstrophy if the advection and Coriolis terms are combined according to \begin{equation} {\bf v} \cdot \nabla {\bf v} + f {\bf k} \times {\bf v} = \nabla \frac{{\bf v}^2}{2} + (\zeta+f) {\bf k} \times {\bf v} \label{eq:B-1} \end{equation} and if the components of $ (\zeta+f) {\bf k} \times {\bf v} $ are written in finite difference form as \begin{equation} -v(\zeta+f) = - \overline{V}^{xy} \overline{Q}^y \hspace{2cm} \label{eq:B-2a} \end{equation} \begin{equation} u(\zeta+f) = \overline{U}^{xy} \overline{Q}^x \hspace{2cm} \label{eq:B-2b} \end{equation} Here, $ (U,V) = (u \Delta \overline{p}^x, v \Delta \overline{p}^y) $ are the horizontal mass flux components, $ Q = (\delta_xv - \delta_y u + f)/ \Delta \overline{p}^{xy} $ is the potential vorticity, and $ \Delta p $ is the layer thickness. Overbar and $ \delta $ operator are defined here in the standard fashion as horizontal two-point averaging and differentiation operators, respectively. Variables are staggered horizontally in such a way that $ u $ and $ v $ points lie midway in $ x $ and $ y $ direction, respectively, between $ \Delta p $ points. This arrangement, which places potential vorticity points in the center of the grid boxes formed by $ \Delta p $ points, is commonly referred to as the C grid [{\em Arakawa and Lamb}, 1977]. (We note in passing that {\em Sadourny} [1975] also gives kinetic energy conserving variants of the above nonlinear operators: \[ -v(\zeta+f) = - \overline{\overline{V}^x Q}^y \]\[ u(\zeta+f) = \overline{\overline{U}^y Q}^x \] Operators which conserve both potential enstrophy and kinetic energy have been developed by {\em Arakawa and Lamb} [1981].) An important point to note with regard to potential vorticity conservation is that the ``ordinary" vorticity equation is the flux form of the potential vorticity equation. This is to say that an area integral of the shallow-water vorticity equation \begin{equation} \frac{\partial}{\partial t} (\zeta+f) + \nabla \cdot [(\zeta+f){\bf v}] = 0 \label{eq:B-3} \end{equation} if taken over a closed domain, yields not only the well-known conservation principle for total absolute vorticity, but expresses conservation of mass-integrated potential vorticity as well. To build this conservation law into the finite-difference momentum equations requires no specific numerical measures other than the assurance that the finite-difference analog of the term $ \nabla {\bf v}^2/2 $ in (\ref{eq:B-1}) vanishes upon application of the finite-difference curl operator. In particular, the finite-difference form of the term $ (\zeta+f) {\bf k} \times {\bf v} $ is irrelevant as long as one is merely concerned with potential vorticity, as opposed to potential enstrophy, conservation. Stated differently, conservation of layer-integrated potential vorticity is virtually ``automatic" as long as the nonlinear momentum advection terms are formulated along the lines of (\ref{eq:B-1}). Aside from benefitting the long-term stability of a shallow-water model integration, {\em Sadourny's} [1975] potential enstrophy conserving operators (\ref{eq:B-2a}) and (\ref{eq:B-2b}) have the desirable property of allowing reduction of the flux form (\ref{eq:B-3}) of the potential vorticity equation to the advective form. With the help of the shallow-water continuity equation \[ \frac{\partial}{\partial t} \Delta p + \delta_x U + \delta_y V = 0 \] the latter is found to be \begin{equation} \frac{\partial Q}{\partial t} + (\Delta \overline{p}^{xy})^{-1} \left( \overline{\overline{U}^{xy} \delta_x Q}^x + \overline{\overline{V}^{xy} \delta_y Q}^y \right) = 0 \label{eq:B-4} \end{equation} Note the appearance of terms $ \delta_x Q, \delta_y Q $: these imply that the potential-enstrophy conserving operators (\ref{eq:B-2a}) and (\ref{eq:B-2b}) in the momentum equations keep the potential vorticity field constant in regions where it is spatially uniform to begin with. An arbitrarily chosen finite difference rendition of the flux form (\ref{eq:B-3}) of the potential vorticity equation will not necessarily offer this feature, i.e., it will conserve potential vorticity only in the layer-integrated sense. The appearance of the variable $ Q $ in (\ref{eq:B-2a}) and (\ref{eq:B-2b}) indicates that potential enstrophy conserving operators can only be used if the minimum thickness of individual coordinate layers is nonzero at all times. In a general-purpose circulation model allowing isopycnals to intersect the sea surface and the bottom topography, this condition is clearly not satisfied. Fortunately, this does not close the door on all attempts to use potential-enstrophy conserving finite difference operators. There are various ways to ``shield" the finite-difference equations in the oceanic interior, where layer thickness remains finite, from the detrimental effects of zero layer thickness. We will describe here one possible modification to the enstrophy-conserving scheme (\ref{eq:B-2a}) and (\ref{eq:B-2b}) which allows it to revert smoothly to a more robust finite difference expression wherever the denominator in the expression for $ Q $ is about to reach zero. The numerical overflow occurring in the calculation of $ Q $ is only one aspect of the problem posed by vanishing layer thickness. Another aspect, and perhaps a physically more relevant one, is that the layer integral of squared, mass-weighted potential vorticity (i.e., the potential enstrophy in the layer) may be infinite if layer depth varies gradually between zero and nonzero values. An infinite reservoir of potential enstrophy clearly makes potential enstrophy conservation useless as a guarantor of long-term numerical stability. We may also expect the unbounded and poorly resolved gradient of $ Q $ in the transition zone between zero and non-zero layer thickness to push the truncation and phase errors in the $ Q $ advection equation (\ref{eq:B-4}) far beyond acceptable limits. In light of these problems, it seems inadvisable to rely on (\ref{eq:B-2a}) and (\ref{eq:B-2b}) when solving the momentum equation in a zone where the layer depth gradually approaches zero. The proper fall-back position is to require that the terms (\ref{eq:B-2a}) and (\ref{eq:B-2b}) revert to a form of $ (\zeta+f) {\bf k} \times {\bf v} $ not involving layer thickness. Our procedure for accomplishing this is as follows. Consider the task of evaluating the term $ - \overline{V}^{xy} \overline{Q}^y $ at a $ u $ velocity point. The required data and their grid locations relative to the $ u $ point are shown in Figure 11. (Data requirements for evaluating $ \overline{U}^{xy} \overline{Q}^x $ at $ v $ points are analogous.) The ``normal" situation, in which it would be appropriate to enforce the potential enstrophy conservation law, is one where all six $ \Delta p $ values shown in Figure 11 are of comparable magnitude and therefore contribute equally to the calculation of the two $ Q $ and the four $ V $ values involved. Let us now study the ``anomalous" situation where some $ \Delta p $ values are significantly larger than others. For symmetry reasons it is sufficient to limit this discussion to the relative magnitudes of the four points labeled $ \Delta p_1 ... \Delta p_4 $ and to assume that $ \Delta p_5, \Delta p_6 $ are zero. {\em Case 1}. $ Three of the four $ \Delta p $ values are large and of comparable magnitude. In this case at least one of the two values $ \Delta p_3, \Delta p_4 $ will belong to the group of large values, guaranteeing that the two $ Q $ values in Figure 11 will be of comparable magnitude. In combination with what amounts to a mass-weighted average of the four surrounding $ v $ values, a reliable estimate of $ -v(\zeta+f) $ at the central location should result. {\em Case 2}. $ Two of the four $ \Delta p $ values are large and of comparable magnitude. If $ \Delta p_3 $ or $ \Delta p_4 $ belongs to this group (case $2a$), the situation is essentially the same as described in case 1. On the other hand, if $ \Delta p_1 $ and $ \Delta p_2 $ are the two large values (case $2b$), the central U point in Figure 11 should be considered too close to the edge of the zero thickness domain to allow enforcement of the potential enstrophy conservation law. We will return to this case after discussing the somewhat simpler situation of case 3. {\em Case 3}. One of the four $ \Delta p $ values greatly exceeds the other three. If either $ \Delta p_3 $ or $ \Delta p_4 $ is the dominant value, the situation is essentially the same as in case 1. No special steps are required then. The problematic case is the one where an outlying $ \Delta p $ value, say, $ \Delta p_1 $, is dominant. In this case three $ V $ values will be negligible compared with the fourth one as they represent conditions at virtually ``massless" grid points. (Note the term ``virtual". If $ \Delta p_3 $ and $ \Delta p_4 $ were exactly zero, there would be no point in solving the $ u $ equation at the central grid point.) This is definitely a situation where the expression $ - \overline{V}^{xy} \overline{Q}^x $ should default to $ -(\zeta+f)v $. Since only one of the four $ V $ values represents motion of a sizeable amount of fluid, this grid point should be singled out for determining the magnitude of the Coriolis term. One way to achieve the goal of giving it full weight (rather than 25\% of the weight as stipulated by $ \overline{V}^{xy} $) is to replace the denominator $ \Delta \overline{p}^{xy} $ at both $ Q $ points by $ \Delta p_1/8 $. We will refer to this algorithm as the ``one-eighth" rule. Returning now to case $2b$, we note that there are now two massless $ V $ points while the other two $ V $ values (the ones in the upper half of Figure 11) will be finite and of comparable magnitude. Application of the one-eighth rule would in this case overestimate the Coriolis force by as much as a factor of 2. This defect can be remedied by using $ (\Delta p_1 + \Delta p_2)/8 $ in the denominator of the two $ Q $ values involved. Note that this modification of the one-eighth rule is compatible with the procedure chosen in case 3. We have streamlined the one-eighth rule as follows. Before computing the $ Q $ field we determine in each grid square the sum of the two largest $ \Delta p $ corner values; we shall denote this sum by $ \Pi_{ij} $ ( $ i,j $ being the grid indices). The denominator in $ Q_{ij} $ is then given the value \begin{equation} \max \left\{ \begin{array}{l} \Delta \overline{p}^{xy}_{ij} \\ \frac{1}{8} \max\,(\Pi_{ij},\Pi_{i-1,j},\Pi_{i+1,j},\Pi_{i,j-1},\Pi_{i,j+1}) \end{array} \right\} \label{eq:B-5} \end{equation} If both alternatives are zero, all possible $ \overline{V}^{xy} \overline{Q}^y $ and $ \overline{U}^{xy} \overline{Q}^x $ combinations involving $ Q_{ij} $ will be comprised of zero $ U $s and $ V $s; in this case there is no need to solve the momentum equations in the vicinity of grid point $ i,j $ because we are too far from grid points containing actual mass. An arbitrary nonzero value in the denominator of $ Q $ may be used in this case to prevent overflow. The one-eighth rule guarantees that layer potential enstrophy will be conserved as long as $ \Delta p $ within octagonal clusters of 12 grid points varies by no more than a factor of four. This can be used to prove numerical consistency of the scheme (\ref{eq:B-5}): As grid resolution increases and the grid points in Figure 11 migrate toward each other, the spatial variation in the $ \Delta p $ field will eventually fall below the factor-of-four threshold. The scheme then reverts to (\ref{eq:B-2a}) and (\ref{eq:B-2b}), whose numerical consistency properties are obvious. While we have tried to take into account as many configurations of the $ \Delta p $ field as possible, the one-eighth rule does not guarantee optimal results under all conditions. However, it does satisfy two important criteria: it avoids (1) complicated decision-making algorithms and (2) temporal discontinuities in the $ Q $ field. \newpage \begin{center} APPENDIX C \\ \bf Kinematic effects of steep bottom slopes \rm \end{center} The continental shelf in an oceanic basin inhibits lateral displacement of abyssal water in a manner similar to a vertical sidewall. In a model that allows all isopycnic layers to extend to the actual shoreline (where the water depth may only be a few tens of meters) the kinematic boundary conditions posed at the coast clearly are relevant only for near-surface water. Some additional boundary conditions expressing the blocking and lateral drag effect of steep bottom slopes are therefore required. One of the most serious problems encountered in the absence of blocking constraints is the occasional extrusion of abyssal water into grid boxes located on the shelf. Depending on the three-dimensional configuration of the bottom topography, katabatic accelerations resulting from these extrusions can lead to numerical instabilities; these appear to be related in part to the need for spatial averaging of the velocity field in the Coriolis term. The blocking technique we have adopted is reminiscent of the Geophysical Fluid Dynamics Laboratory (GFDL) model ``step" mountain. Rather than assuming that the water depth varies linearly between adjacent mass field grid points, we treat the bottom within each mass field grid box as flat. Thus a ``sill" is defined on the boundary between two adjacent grid boxes that have different water depth. The depth of the flow over the sill is then given by the minimum (as opposed to the average) of the water depths in the upstream and downstream box. The pressure on interior interfaces at $ u $ and $ v $ points, however, is still given by $ \overline{p}^x $ and $ \overline{p}^y $, respectively, provided these values do not exceed the sill depth as defined above. To summarize, the depth of the $ k $th interface at, say, a $ u $ velocity point labeled $ i,j $ is given by \[ \min\,(\frac{p_{i-1/2,j,k} + p_{i+1/2,j,k}}{2}, p_{i-1/2,j,bot}, p_{i+1/2,j,bot}) \] where the index $ i $ varies in $ x $ direction and the subscript ``bot" denotes bottom pressure. A lateral drag condition (free-slip or no-slip) is applied to flow tangent to submerged orographic obstacles. Recall that lateral drag is incorporated in the momentum equations (\ref{eq:A-1a}) and (\ref{eq:A-1b}) through the eddy momentum flux $ \nu \Delta p \: \nabla {\bf v} $ . Consider the exchange of $ u $-momentum in $ y $ direction between two grid boxes, $ (i,j) $ and $ (i,j-1) $ (abbreviated here as $ j $ and $ j-1 $). If point $ j $ is a boundary point, that is, if a sidewall is present at $ j - \frac{1}{2} $, boundary conditions are invoked by setting $ (\partial u / \partial y) \Delta y $ at $ j-\frac{1}{2} $ to $ u_j - \sigma u_j $, where $ \sigma = 1 $ for free-slip and $ \sigma = -1 $ for no-slip conditions. If a ``partial" sidewall is present in the form of a seafloor rise between $ j $ and $ j-1 $, the above drag condition is generalized by writing \[ \frac{\partial u}{\partial y} \Delta y = u_j - [w u_{j-1} + (1-w) \sigma u_j] \] where the aspect ratio $ 1-w $, $ 0 \leq w \leq 1 $ indicates the extent to which the seafloor elevation at $ j-1 $ creates a sidewall if viewed from grid box $ j $. The vectorized code for implementation of sidewall friction has been found to account for roughly 6\% of the total computation time for a 2-layer model experiment. \newpage {\em Acknowledgments}. Support for this work by the National Science Foundation under grants OCE-8600593 and OCE-8812185 is gratefully acknowledged. Most computations were carried out at the National Center for Atmospheric Research, sponsored by the National Science Foundation. Computer time was also made available in 1985 by the Geophysical Fluid Dynamics Laboratory, Princeton. We thank Douglas B. Boudra (deceased) and Claes Rooth for their continuing interest in this work and for many helpful discussions. \newpage \begin{center} \bf REFERENCES \rm \end{center} \parindent 0pt Anderson, D. L. T., and P. D. Killworth, Spin-up of a stratified ocean, with topography. {\em Deep-Sea Res., 24}, 709--732, 1977. Arakawa, A., and V. R. Lamb, Computational design of the basic processes of the UCLA general circulation model. {\em Methods of Computational Physics, 17}, Academic Press, New York, 174 -- 265, 1977. ---, and ---, A potential enstrophy- and energy-conserving scheme for the shallow water equations. {\em Mon. Wea. Rev., 109}, 18 -- 36, 1981. Bleck, R., Finite difference equations in generalized vertical coordinates. Part I: Total energy conservation. {\em Contrib. Atm. Phys., 51}, 360 -- 372, 1978. ---, and D.B. Boudra, Initial testing of a numerical ocean circulation model using a hybrid (quasi-isopycnic) vertical coordinate. {\em J. Phys. Oceanogr., 11}, 755 -- 770, 1981. ---, and ---, Wind-driven spin-up in eddy-resolving ocean models formulated in isopycnic and isobaric coordinates. {\em J. Geophys. Res., 91C}, 7611 -- 7621, 1986. ---, H. P. Hanson, D. Hu, and E. B. Kraus, Mixed layer/thermocline interaction in a three-dimensional isopycnal coordinate model. {\em J. Phys. Oceanogr.}, in press (Sept. issue), 1989. Bryan, K., A numerical method for the study of the circulation of the world ocean. {\em J. Comput. Phys., 4}, 347 -- 376, 1969. Gadd, A. J., A split-explicit integration scheme for numerical weather prediction. {\em Quart. J. R. Met. Soc., 104}, 569 -- 582, 1978. Holland, W. R., On the wind-driven circulation in an ocean with bottom topography. {\em Tellus, 19}, 582 -- 599, 1967. Huang, R. X., Numerical simulation of wind-driven circulation in a subtropical/ subpolar basin. {\em J. Phys. Oceanogr., 16}, 1636 -- 1650, 1986. ---, A three-layer model for wind-driven circulation in a subtropical/subpolar basin. Part I: Model formulation and the subcritical state. {\em J. Phys. Oceanogr., 17}, 664 -- 678, 1987. Luyten, J. R., J. Pedlosky, and H. Stommel, The ventilated thermocline. {\em J. Phys. Oceanogr., 13}, 292 -- 309, 1983. Mesinger, F., On the convergence and error problems of the calculation of the pressure gradient force in sigma coordinate models. {\em Geophys. Astrophys. Fluid Dyn., 19}, 105 -- 117, 1982. Montgomery, R. B., Circulation in upper layers of southern North Atlantic deduced with use of isentropic analysis. {\em Papers in Physical Oceanography and Meteorology, VI}, Woods Hole Oceanographic Institution, Woods Hole, Mass., 55 pp., 1938. Parsons, A. T., A two-layer model of Gulf Stream separation. {\em J. Fluid Mech., 39}, 511 -- 528, 1969. Sadourny, R., The dynamics of finite-difference models of the shallow-water equations. {\em J. Atmos. Sci., 32}, 680 -- 689, 1975. Thompson, J. D., and W. J. Schmitz, Jr., A limited-area model of the Gulf Stream: Design, initial experiments, and model-data intercomparison. {\em J. Phys. Oceanogr., 19}, 791 -- 814, 1989. Welander, P., A two-layer frictional model of wind-driven motion in a rectangular oceanic basin. {\em Tellus, 18}, 54 -- 62, 1966. Zalesak, S., Fully multidimensional flux-corrected transport algorithms for fluids. {\em J. Comput. Phys., 31}, 335 -- 362, 1979. \newpage \bf Figure Captions \rm Fig. 1. Basin configuration and dimensions. Dashed line indicates ``low-cut" shelf configuration. Fig. 2. Barotropic mass-flux stream function obtained in 2-layer, low-cut topography experiment. Contour interval: 4$\: \times \:10^6m^3 s^{-1} $. Solid (dashed) lines denote clockwise (counterclockwise) flow. Fig. 3. Barotropic mass-flux stream function obtained in 2-layer, high-cut topography experiment. Contour interval: 4$\: \times \:10^6m^3 s^{-1} $. Fig. 4. Total kinetic energy (top curve), layer-1 and layer-2 kinetic energy (lower 2 curves) versus time in 2-layer, high-cut experiment using pressure gradient scheme 1 (see Appendix A). Fig. 5. As in Fig. 4, but using pressure gradient scheme 2. Fig. 6. Total kinetic energy (top curve), layer-1, layer-2 and layer-3 kinetic energy (lower 3 curves) versus time in 3-layer, flat-bottom experiment. Vertical lines indicate times shown in Fig. 7. Fig. 7. Depth of the layer-1 lower interface at three times marked by vertical lines in Fig. 6. Solid (dashed) lines denote interface depth greater (less) than the $ 300 \: m $ depth at time zero. Contour interval $ 50 \: m $. Fig. 8. As in Fig. 6, but for the 3-layer, high-cut topography experiment. Vertical lines indicate times shown in Fig. 9. Fig. 9. As in Fig. 7, but for the three times marked by vertical lines in Fig. 8. Fig. 10. Vertical section showing isopycnic layers near sloping bottom. Fig. 11. Grid variables used in evaluating the Coriolis term $fv$ at the central $ U $ point. \end{document}