\documentstyle [11pt]{article} \oddsidemargin 0.1in \evensidemargin 0.1in \topmargin -0.5in \headsep .5in \textheight 9in \textwidth 6.5in \parskip 7pt \def\baselinestretch{1.7} \begin{document} \begin{center} \vspace*{5cm} {\large \bf Salinity-Driven Thermocline Transients in a Wind- and Thermohaline-Forced Isopycnic Coordinate Model of the North Atlantic} \\ \vspace{3cm} Rainer Bleck, Claes Rooth, Dingming Hu, and Linda T. Smith \\ \vspace{3cm} University of Miami \\ Rosenstiel School of Marine and Atmospheric Science \\ Miami, Florida 33149 \\ \vspace{3cm} FINAL VERSION \\ February 1992 \end{center} \newpage \begin{center} ABSTRACT \end{center} An isopycnic-coordinate oceanic circulation model formulated with the aim of simulating ther\-mo\-dynamically- and mechanically-driven flow in realistic basins is presented. Special emphasis is placed on the handling of diabatic surface processes and on thermocline ventilation. The model performance is illustrated by a 30-year spinup run with coarse horizontal resolution ($2^0$ mesh) in a domain with North Atlantic topography extending from $10^0$N to $60^0$N latitude. The vertical structure encompasses 10 isopycnic layers in steps of 0.2 $\sigma$ units, capped by a thermodynamically active mixed layer. From an initially isohaline state with isopycnals prescribed by zonally averaged climatology, the model is forced by seasonally varying wind stress, radiative and fresh-water fluxes, and by a thermal relaxation process at the surface. After a mechanical spinup time of about 15 years, a quasi-stationary pattern of mean circulation and annual variability ensues, characterized by pronounced subtropical mode-water formation and a gradual growth in the salinity contrast between the subtropics and the subpolar region. The effect of the fresh-water flux forcing on the ventilation of the thermocline is a key point of discussion. Finally, a low-viscosity experiment suggests that the thermohaline processes represented in the model are quite insensitive to dynamic noise development at the grid resolution limit. \newpage {\bf 1. INTRODUCTION} In contrast to the atmosphere, which is highly transparent to solar radiation and is therefore primarily driven by buoyancy input at its lower boundary, the world ocean is a rather opaque medium forced from above. Diabatic effects in its interior are mainly due to small-scale mixing. If approximated as diapycnic eddy diffusion, they demand diffusivity values of the order of $10^{-4} \mbox{m}^2 \mbox{s}^{-1}$. This implies that on decadal and shorter time scales, interior diabatic processes are likely to be significant on vertical scales of not more than a few hundred meters -- a prediction which is borne out by interleaving of trace property distributions in the deep basins, and by a high degree of lateral coherence of water mass characteristics observed when properties are mapped against potential density. Recognition of this quasi-adiabatic character of oceanic motion on fairly long time scales has encouraged the use of layer models in studying dynamic interaction modes, but a variety of technical difficulties has impeded their development into full-fledged oceanic general circulation models (OGCMs), the main one being the need to use computation-intensive methods for handling the movable intersections of coordinate surfaces with domain boundaries (Bleck and Boudra, 1981). Nevertheless, the intrinsic capacity of layer models to dynamically adapt their spatial vertical resolution to the occurrence of baroclinic zones, and the prospect of being able to handle long-term climate-oriented OGCM integrations in a framework where the diabatic exchange processes are totally under the experimenter's control, has encouraged the University of Miami modelling group to pursue the development of a practical isopycnic OGCM code. [A global isopycnic OGCM code featuring implicit time integration has been developed independently by Oberhuber (1990).] Our work has proceeded along two lines. The first one encompassed the development of an adiabatic, purely wind-forced model which resolves the various boundary-related difficulties, including the handling of variable bottom topography (Boudra {\em et al.}, 1988; Bleck and Smith, 1990; Smith {\em et al.}, 1990). The second line of work dealt with thermodynamic mixed-layer forcing and the associated treatment of inter-layer mass transfer in an idealized geometric setting (Bleck {\em et al.}, 1989), first with one and then with two prognostic thermodynamic variables linked by a nonlinear equation of state. We present here a synthesis of those two lines of development. Although a complete OGCM must include property and mass exchange mechanisms that operate between interior isopycnic layers, there exists a range of physical conditions and time scales for which the dominant diapycnal interactions are associated with the annual growth and retreat of the mixed layer. This applies to the ventilation of the permanent thermocline on sub-decadal time scales, and also to the impact of persistent upwelling processes on the state of the surface layer. Considering this, and the fact that the implementation of interior diapycnic processes is clearly separable (in this model class) from that of mixed-layer oriented processes, the present time is appropriate for the publication of a systematic model description, along with an account of a set of semi-realistic model experiments. The basic experiment described here is a 30-year run, undertaken with the dual purpose of illustrating the dynamics of surface-forced ventilation and of providing a reference case for future model experiments aimed at elucidating the significance of diapycnal processes in the interior. Since the latter experiments are intended to focus specifically on North Atlantic conditions, a low-resolution Atlantic domain was chosen for the present work. An exaggerated fresh-water exchange rate (evaporation minus precipitation) was adopted to test the model's ability to simulate dynamic regimes dominated by salinity contrasts. \newpage {\bf 2. MODEL OVERVIEW} In the present model, the oceanic interior is mapped into a set of isopycnal layers, giving it the performance characteristics of a ``stacked shallow water" model. A surface mixed layer of the Kraus and Turner (1967) type provides the linkage between externally imposed wind and thermohaline forces and the isopycnic grid domain. The annual forcing cycle causes the mixed layer to expand and shrink in the course of a model year, causing the mixed layer to exchange mass and other properties with interior layers and to thereby ventilate the latter. Contrary to some wind-driven isopycnic models which gain computational efficiency by keeping the lowermost coordinate layer at rest, thereby eliminating barotropic gravity waves, the present model permits motion in all layers and time-integrates the barotropic flow component in ``split-explicit" mode. With motion being allowed in the abyssal layers, much of our development effort has been directed at stabilizing the numerical solution in places where the flow interacts with topography. Another challenge was the development of transport algorithms that maintain their conservative properties, including positive-definiteness, in the presence of significant spatial and temporal variations in layer thickness. In this section we summarize the salient features of the model. While we have attempted to make the description as self-contained as possible, we will refer the reader for technical details to previous publications. These will be denoted by the letters BB (Bleck and Boudra, 1986), BHHK (Bleck, Hanson, Hu, and Kraus, 1989), BS (Bleck and Smith, 1990), and SBB (Smith, Boudra, and Bleck, 1990). Technical innovations too recent to be found in those papers will be presented in Appendices to the present paper. {\em 2.1 Governing equations.} The Miami model is a primitive-equation model containing 4 prognostic equations -- one for the horizontal velocity vector, a mass continuity or layer thickness tendency equation, and two conservation equations for salt and heat. The equations, written here for a generalized vertical coordinate ``$s$" (Bleck, 1978) to make them formally applicable to both the isopycnic domain and the surface mixed layer, are \begin{equation} \frac{\partial {\bf v}}{\partial t_s} + \nabla_s \frac{{\bf v}^2}{2} + (\zeta + f){\bf k} \times {\bf v} + \left( \dot{s} \frac{\partial p}{\partial s} \right) \frac{\partial {\bf v}}{\partial p} + \nabla_{\alpha} M = - g \frac{\partial \mbox{\boldmath $\tau$} }{\partial p} + \left( \frac{\partial p}{\partial s} \right)^{-1} \nabla_s \cdot \left( \nu \frac{\partial p}{\partial s} \,\nabla_s {\bf v} \right) \label{eqn:momentum} \end{equation} \begin{equation} \frac{\partial}{\partial t_s} \left( \frac{\partial p}{\partial s} \right) + \nabla_s \cdot \left( {\bf v} \frac{\partial p}{\partial s} \right) + \frac{\partial}{\partial s} \left( \dot{s} \frac{\partial p}{\partial s} \right) = 0 \label{eqn:continuity1} \end{equation} \begin{equation} \frac{\partial}{\partial t_s} \left( \frac{\partial p}{\partial s} T \right) + \nabla_s \cdot \left( {\bf v} \frac{\partial p}{\partial s} T \right) + \frac{\partial}{\partial s} \left( \dot{s} \frac{\partial p}{\partial s} T \right) = \nabla_s \cdot \left( \nu \frac{\partial p}{\partial s} \,\nabla_s T \right) + {\cal H}_T \label{eqn:transport1} \end{equation} Here, ${\bf v} = (u,v)$ is the horizontal velocity vector, $p$ is pressure, $T$ represents any one of the model's thermodynamic variables (temperature, salinity, buoyancy), $\alpha = \rho^{-1}$ is the specific volume, $\zeta = \partial v / \partial x_s - \partial u / \partial y_s$ is the relative vorticity, $M \equiv gz + p \alpha$ is the Montgomery potential, $gz$ is the geopotential, $f$ is the Coriolis parameter, ${\bf k}$ is the vertical unit vector, $\nu$ is a variable eddy viscosity coefficient, and {\boldmath $\tau$} is the wind- and/or bottom-drag induced shear stress vector. ${\cal H}_T$ represents the sum of the diabatic source terms acting on $T$. Subscripts indicate which variable is held constant during partial differentiation. Distances in $x,y$ direction (as well as their time derivatives $\dot{x} \equiv u$ and $\dot{y} \equiv v$) are measured in the projection onto a horizontal plane, a convention which eliminates metric terms related to the slope of the $s$ surface. After vertical integration over a coordinate layer bounded by two surfaces $s_{top}, s_{bot}$ (the ocean surface, the bottom of the mixed layer, the sea floor, and all isopycnic layer interfaces are regarded as $s$ surfaces here), the continuity equation (\ref{eqn:continuity1}) becomes a prognostic equation for the layer weight per unit area, $\Delta p = p_{bot} - p_{top}$: \begin{equation} \frac{\partial}{\partial t_s} \Delta p + \nabla_s \cdot ({\bf v} \Delta p) + \left( \dot{s} \frac{\partial p}{\partial s} \right)_{bot} - \left( \dot{s} \frac{\partial p}{\partial s} \right)_{top} = 0. \label{eqn:continuity} \end{equation} The expression $(\dot{s} \partial p / \partial s)$ represents the vertical mass flux -- taken to be positive if in $+p$ or downward direction -- across an $s$ surface. Multiplication of (\ref{eqn:momentum}) by $\partial p/\partial s$ and integration over the interval $(s_{top}, s_{bot})$, followed by division by $\Delta p / \Delta s$, changes the shear stress term in that equation into \[ \frac{g}{\Delta p} ( \mbox{\boldmath $\tau$}_{top} - \mbox{\boldmath $\tau$}_{bot} ) \] while the lateral momentum mixing term integrates to \begin{equation} (\Delta p)^{-1} \nabla_s \cdot (\nu \Delta p \, \nabla_s {\bf v}). \label{eqn:diffusionterm} \end{equation} All other terms in (\ref{eqn:momentum}) retain their formal appearance. Wind-induced stress is assumed to be zero at the bottom of the mixed layer, which therefore must not be allowed to become shallower than the Ekman layer. Bottom stress is prorated among coordinate layers in accordance with the assumption that it linearly decreases from -$\rho c_D |{\bf v}| {\bf v}$ to zero over the lowest $10\,m$ of the water column. ${\bf v}$ in this formula is the velocity vector averaged over the lowest $10\,cm$. The layer-integrated form of (\ref{eqn:transport1}) is \begin{equation} \frac{\partial}{\partial t} T \Delta p + \nabla_s \cdot ({\bf v} T \Delta p) + \left( \dot{s} \frac{\partial p}{\partial s} T \right)_{bot} - \left( \dot{s} \frac{\partial p}{\partial s} T \right)_{top} = \nabla_s \cdot (\nu \Delta p \,\nabla_s T) + {\cal H}_T. \label{eqn:transport} \end{equation} The above prognostic equations are complemented by diagnostic equations including the hydrostatic equation \begin{equation} \frac{\partial M}{\partial \alpha} = p, \label{eqn:hydrostat} \end{equation} an equation of state linking $T,S$ to $\alpha$, and a turbulence kinetic energy balance equation (outlined in Sec. 3a of BHHK) yielding the depth of the mixed layer. {\em 2.2 Horizontal and vertical grid.} We consider isotropy in horizontal grid resolution an important aspect of model numerics and therefore solve the equations on a regular mesh superimposed on a conformal (presently: Mercator) projection of the earth's surface. Thus, for a mesh size of $n$ degrees longitude the mesh size in meridional direction at latitude $\varphi$ is $(n\,\cos \varphi)$ degrees. The map-scale factor is a function of $\varphi$, but not of direction: small grid squares on the sphere project onto (nearly) square areas on the map. Velocity and mass field variables are staggered horizontally such that $u$ ($v$) grid points are positioned halfway between mass grid points in $x$ ($y$) direction. This staggering arrangement, known as the Arakawa C grid, is not considered optimal for coarse-mesh (noneddy-resolving) models (Arakawa and Lamb, 1977) and differs from the one employed in BHHK, but was chosen with future eddy-resolving applications in mind. Thermodynamic variables and ${\bf v}$ are treated as ``layer" variables which are vertically constant within layers but change discontinuously across layer interfaces. $p,z$, and $\dot{s} \partial p /\partial s$, on the other hand, are ``level" variables defined on interfaces. $M$ turns out to be a layer variable. {\em 2.3 Time differencing scheme.} In contrast to the model described in BHHK, which filters gravity waves using a rigid lid, the present model uses the split-explicit scheme described in BS to extend the time step in the baroclinic calculations. The barotropic component of the solution is advanced in time using the forward-backward scheme which uses the ``old" mass field to solve the continuity equation and then uses the ``new" pressure field in the momentum equations. A substantially longer time step is used in the leapfrog scheme advancing the baroclinic part of the solution. The two schemes interact once every baroclinic time step. Explicit integration of the barotropic mode was not only found to be faster than the rigid-lid scheme relying on a traditional iterative Poisson solver for the barotropic stream function, but also has the advantage of producing a sea-surface height field allowing direct comparison with, and possibly assimilation of, radar altimetry data. Liberal use is being made in the model of the time-splitting concept. For example, (\ref{eqn:transport}) is actually decomposed into 4 separate equations which in turn advect the $T$ field horizontally and vertically, diffuse it, and force it adiabatically. Vertical advection, being caused in part by cabbeling and by mixed-layer entrainment/detrainment, is itself a 2-step process. {\em 2.4 Isopycnal outcropping.} Outcropping of isopycnic coordinate layers is the Achilles' heel of isopycnal modeling. Our basic method for dealing with this process has not changed since BB. Briefly speaking, we assume each isopycnic layer to exist at every geographic location in the model, but allow it to become massless or dynamically ``invisible". Inflation or deflation of each layer, and thus the geographic location where the layer outcrops at a given time, is entirely controlled by the interplay between the continuity and momentum equations. Special algorithms maintaining positive-definiteness of layer thickness are required to make this scheme work; these are somewhat more expensive than conventional finite-difference methods for solving the continuity equation but do not require a special bookkeeping effort to keep track of layer outcrop lines or time-dependent lateral boundaries [the method chosen by Oberhuber (1990)]. We will return to this topic in Section 2.10. While the isopycnal layers in BB and SBB outcrop at the sea surface, they outcrop into the mixed layer in BHHK and in the present model -- a minor technical difference from the numerical point of view. Variable bottom topography requires that they also be allowed to intersect the ocean bottom. The bottom intersection problem is numerically similar to the surface outcrop problem. More on this subject can be found in BS. {\em 2.5 Surface forcing functions.} As in BHHK, the model is forced by $(i)$ wind stress; $(ii)$ mechanical stirring; $(iii)$ an imposed (i.e., model state-independent) thermal energy flux representing solar radiation; and $(iv)$ a relaxation boundary condition representing the combined effect of sensible and latent heat fluxes. $(iii)$ and $(iv)$ are presently designed to capture the meridional and temporal variability of the observed surface fluxes over the North Atlantic but are zonally averaged. Salinity, which was held constant in BHHK, is being forced by an imposed latitude-dependent (but time-independent) fresh-water flux. We selected this type of forcing, which differs from the commonly employed Newtonian relaxation toward observed surface salinity, for two reasons: $(a)$ A relaxation boundary condition is nonphysical because evaporation is virtually independent of surface salinity and precipitation is totally so. $(b)$ Imposed boundary fluxes are more likely to generate exotic salinity features; the model's reaction to such features provides valuable insight into its ability to simulate climates different from the present one. {\em 2.6 Mixed layer.} The thermodynamically active mixed layer and its modes of interaction with the interior isopycnic coordinate layers have been thoroughly documented in BHHK. However, the addition of a second thermodynamic variable has created new degrees of freedom in the detrainment algorithm. Alterations to the detrainment logic necessitated by the interplay between two independently acting thermodynamic variables are sketched below. As outlined in BHHK, discharge of mixed-layer water into the isopycnic interior during periods of mixed-layer retreat is accompanied by vertical buoyancy transfer within the model mixed layer. This buoyancy transfer is needed to adjust the density of the detrained water to the density of the receiving isopycnal layer; it is a truncation error in the sense that it is caused by finite grid resolution and that it can be suppressed to any degree desired by adding more layers to the model. To maintain static stability, the buoyancy transfer must always be directed upward, i.e., heat must be extracted from the column about to be detrained and added to the remaining, shallower mixed layer. To keep the latter from overheating, which would bias surface heat fluxes controlled by a relaxation boundary condition, a limit is placed on the mixed-layer heating rate. The price to be paid for this is an often substantially reduced detrainment rate: part of the ``fossil" mixed layer may for some length of time remain part of the model mixed layer, effectively extending the {\em physical} mixed layer whose depth in the Kraus-Turner model is given by the Monin-Obukhov length $L$. One obvious way to cap the vertical heat flux is to demand that the temperature rise in the model mixed layer (whose depth for reasons just explained may exceed $L$) be equal to the rate obtained by distributing the incoming thermal energy over the depth range $-L \leq z \leq 0$. This approach confines the numerically required heat transfer in the water column to water parcels below $-L$ and guarantees that the model ocean always displays the ``correct" sea surface temperature. The model catches up and restores the proper mixed-layer depth each time the mixed-layer density during springtime warming matches the density of an interior coordinate layer. As illustrated in BHHK, the exact moment for such detrainment ``bursts" comes when the mixed-layer buoyancy rises from $b_1$, a value just below the buoyancy $b_{k-1}$ of an interior layer, to a value $b'_1 \geq b_{k-1}$. An arbitrary amount of mixed-layer water can then be discharged into layers $k$ and $k$--1 without running the risk of overheating the mixed layer. We refer to these bursts as three-way detrainment events because they involve redistribution of mass among three layers, in contrast to the more frequent two-way detrainment events. It was shown in BHHK that apportioning the detrained column among layers $k$ and $k$--1 during a three-way detrainment event is a relatively easy task if buoyancy is a function of only one thermodynamic variable. In this special case the detrained column can be separated by internal heat transfer into subcolumns of buoyancies $b_k$ and $b_{k-1}$, with the assurance that detraining these into the interior will produce water of exactly the right density in the receiving layers. If both temperature and salinity are allowed to vary, ``cabbeling" (see Section 2.8) will be a factor and the above apportioning method will no longer be exact. The proper way of resolving the apportioning problem in the case of two thermodynamic variables is to specify that the heat required to raise the buoyancy in the upper part of the detrained column and in the remaining mixed layer from $b_1$ to $b_{k-1}$ and $b'_1$, respectively, be equal to the heat that needs to be extracted from layer $k$ to restore $b_k$ after the lower part of the detrained column (of buoyancy $b_1 > b_k$) has been mixed into it. A scheme that solves the above problem is described in Appendix E. This scheme has only recently been developed and tested. When we conducted the simulation experiment described in Section 3 we used the buoyancy-based apportioning method from BHHK with the following modifications. Water apportioned to layer $k$ by this method is tentatively mixed with layer-$k$ water. If the heat extracted from the mixture to restore its buoyancy to $b_k$ is not sufficient to raise the buoyancy of the remaining water to $b_{k-1}$, the deficit is taken from the mixed layer. If this causes the mixed-layer buoyancy to fall below $b_{k-1}$, the algorithm reverts from three-way to two-way detrainment, i.e., a small amount of water will be discharged into layer $k$ alone. Depending on the resulting mixed-layer buoyancy, three-way detrainment may still be feasible at the next time step; otherwise, two-way detrainment involving layer $k$--1 will take place. Safeguards are also put in place to check the direction of the buoyancy flux during two-way detrainment. Specifically, if cabbeling in the mixture of interior-layer water and detrained mixed-layer water requires a downward transfer of heat, and if this transfer makes the mixed layer heavier than the receiving layer, detrainment is postponed till the next time step. The effect of mixed-layer entrainment/detrainment on the momentum field is modeled by transferring momentum along with mass. Note that momentum is a {\em layer} variable in shallow-water models. {\em 2.7 Equation of state.} The equation of state used in BHHK is linear in salinity ($S$) and quadratic in temperature ($T$), which allows it to be easily inverted into an equation for $T$ as a function of $\rho$\footnote {Throughout this paper, $\rho \equiv \alpha^{-1}$ and $\sigma = \rho$--1000\,kg\,m$^{-3}$ denote potential density referred to sea level pressure.} and $S$. The main problem with this quadratic formula is that it renders the equation of state too nonlinear at moderate-to-high temperatures and too linear at low temperatures. This has the effect of exaggerating the cabbeling process at middle-to-low latitudes while suppressing it in subpolar regions. For this reason we have replaced the quadratic equation of state in BHHK by the $3^{rd}$-degree polynomial approximation proposed by Friedrich and Levitus (1972). Finding the roots of a $3^{rd}$-degree polynomial involves repeated calls to trigonometric functions to avoid complex-valued arithmetic and is therefore considerably more time-consuming than finding the roots of a quadratic equation. Execution time in our model, which uses the inverted equation of state in the cabbeling and the mixed-layer detrainment algorithm, increased by 12\% when the switchover took place. Of the three roots of the polynomial, the one guaranteed to be real lies near $180^0C$ and is therefore nonphysical. The other two roots, if they exist, form a pair straddling a temperature value in the range $-3^0$ to $-5^0C$ where $\rho$ for seawater reaches its maximum; the root above this value is the physical one and is used in computing $T(\rho,S)$. {\em 2.8 The Cabbeling algorithm.} ``Cabbeling" is the term used to describe the density increase resulting from the mixing of two water parcels of equal density but different temperature and salinity. The German verb {\em kabbeln}, which this term is derived from, indicates that such mixing usually produces a ``choppy" or fine-grain vertical velocity field. Cabbeling in a coarse-mesh ocean model is therefore a severely aliased process. In an isopycnic coordinate model, cabbeling generates a positive vertical velocity component $d \rho/d t$, i.e., a gradual transfer of mass from each coordinate layer to the one beneath it. In the Miami model the cabbeling process is treated as part of the more general problem of maintaining the reference density in coordinate layers where the $T,S$ fields evolve independently of each other. The numerical procedure adopted for this purpose is described in Appendix A. {\em 2.9 Thermodynamic equations.} As mentioned earlier, the model accomodates two independent thermodynamic variables. The prognostic equations for these variables (buoyancy $b$ and salinity $S$ in the mixed layer, temperature $T$ and $S$ in the interior) are all of the form (\ref{eqn:transport}). Horizontal advection (note that ``horizontal" translates into isopycnal everywhere in the model except in the mixed layer) is handled by the conservative, monotony-preserving, $3^{rd}$-order version of the Smolarkiewicz scheme (Smolarkiewicz and Clark, 1986; Smolarkiewicz and Grabowski, 1990). Details on the formulation of the turbulent mixing terms are given in Sections 2.12 and 2.13. There are presently two processes in the model -- cabbeling and mixed-layer entrain\-ment/\-detrainment -- that give rise to vertical mass fluxes $\dot{s} \partial p / \partial s$ and thereby allow thermodynamic variables to be transported across layer interfaces. {\em 2.10 Continuity equation.} The split-explicit time integration scheme requires that changes in layer thickness caused by the barotropic velocity component $\overline {\bf v}$ be evaluated separately from those caused by the baroclinic component ${\bf v}'$. This scheme has been described in considerable detail in BS. For the reader's convenience we recapitulate its major features in Appendix B. {\em 2.11 Momentum equations.} Even though the nonlinear terms in (\ref{eqn:momentum}) play no significant role in coarse-mesh basin-scale simulations such as the one reported here, the equations are written in their fully nonlinear form to accomodate future fine-mesh applications. As outlined in detail in Appendix B of BS, the Coriolis and nonlinear terms conserve potential vorticity and potential enstrophy except in the vicinity of zero layer thickness zones. The mixing terms are written in conservative form (see Sections 2.12 and 2.13) and use an eddy viscosity that reverts from a constant value to one proportional to the total deformation in regions of large horizontal shear (Smagorinsky, 1963). Specifically, the eddy viscosity used in isopycnal mixing of momentum is defined as \begin{equation} \nu = max \left( u_d \Delta x, \eta [(u_x-v_y)^2 + (v_x+u_y)^2]^{1/2} \Delta x^2 \right) \label{eqn:eddyvisc} \end{equation} where $\Delta x$ is the mesh size, subscripts denote partial derivatives, $u_d$ is a background diffusive velocity, and $\eta$ is a dimensionless factor presently set to 2.0. Near steep bottom topography, the horizontal pressure force $\nabla_{\rho} M$ is computed in non-centered form and thus does not require algorithms for extrapolating pressure to underground grid points. Details of the scheme, which basically infers a layer-depth weighted average of $\nabla_{\rho} M$ from points surrounding the grid point in question, are provided in Appendix A of BS. Lateral sidewalls are usually chosen to coincide with grid points of the {\em normal} velocity component to allow enforcement of the kinematic boundary condition. In the C grid, this implies that the tangential velocity grid point nearest a sidewall is displaced from the wall by half a grid distance. A ``mirror image" of this grid point on the opposite side of the wall is needed to solve the momentum equation at this location. Depending on how one defines the velocity at the mirror image point, one can run a C grid model with free-slip or no-slip lateral boundary conditions. The results reported in this paper have been obtained with no-slip conditions. As discussed elsewhere, we follow the convention of extending coordinate layers intersecting the ocean bottom as massless layers all the way to the coast. This convention, in conjunction with the layer depth-weighted lateral mixing formulation discussed in Section 2.13, makes it nearly impossible to communicate lateral sidewall drag to, say, grid points in the abyssal ocean next to the continental rise. A ``submerged sidewall" boundary condition has been developed to remedy this problem. Details can be found in Appendix C of BS. Since the momentum equations are written in advective rather than flux form, they can in principle be used to determine $u,v$ values at massless grid points. (Such values are needed to carry out finite-difference calculations at nearby non-massless grid points). However, it is safer in our experience to infer $u,v$ at massless points from non-massless points above or below in the column. All velocity values in the model are therefore replaced after each time step by vertical averages over a small depth interval centered on the layer in question. This procedure has no effect on the velocity values in locations where the layer thickness exceeds 10\,cm. {\em 2.12 Subgrid-scale mixing.} Let $T$ denote one of the prognostic variables in the model. Since the finite-difference analog of the product $\nu \partial T/ \partial x$ is numerically equivalent to the finite-difference analog of $u_d \partial T$ (where $\partial T$ is the $T$ difference between neighboring grid points, $u_d$ is a ``diffusive flux velocity" of magnitude $\nu /\partial x$, and $\partial x$ is the grid size), we have in BHHK and BS adopted the practice of specifying $u_d$ rather than $\nu$ as the parameter controlling isopycnal mixing. Setting $u_d$ constant throughout the model domain has the advantage of removing any geographic bias from the relative weight of diffusive and advective processes in shaping the character of the numerical solution. The $u_d$ value also provides a convenient measure of the degree to which diffusion overshadows advection. The isopycnic mixing term in the two thermodynamic equations is presently formulated with $u_d$ = 1\,cm\,s$^{-1}$ while a value of 2\,cm\,s$^{-1}$ is used for momentum mixing. Subgrid-scale turbulent flux terms in the prognostic differential equations formally result from Reynolds averaging of the various nonlinear flux terms. Since the term $\nabla \cdot ({\bf v} \Delta p)$ in (\ref{eqn:continuity}) is nonlinear and thus breaks down into \[ \overline{\nabla \cdot ({\bf v} \Delta p)} = \nabla \cdot (\overline{{\bf v}} \overline{\Delta p}) + \nabla \cdot (\overline{{\bf v}' \Delta p'}), \] upon Reynolds averaging, it would seem inconsistent not to add a ``thickness diffusion" term to (\ref{eqn:continuity}). Following the parameterization of subgrid-scale fluxes in the other model equations, we express the turbulent thickness diffusion term $\nabla \cdot (\overline{{\bf v}' \Delta p'})$ in the form $-\nabla \cdot (\nu \nabla \,\Delta p)$ where $\Delta p$ now represents a grid-scale average. Our numerical implementation of thickness diffusion is described in Appendix C. {\em 2.13 Transport and diffusion in the presence of massless grid points.} The appearance and disappearance of massless grid points requires special steps to assure conservation of dependent variables during isopycnal transport and mixing because conservative schemes require division by layer thickness. We will first discuss this problem and our adopted remedy in the context of grid-scale transport (``grid-scale" transport processes are those associated with resolved scales of motion, i.e., advection by the ${\bf v}$ field). The issue in this case is the interplay between the layer thickness tendency equation (\ref{eqn:continuity}) and the transport equation (\ref{eqn:transport}). The Flux Corrected Transport scheme used to solve (\ref{eqn:continuity}) (see Section 2.10) is capable of fine-tuning the horizontal mass fluxes ${\bf v} \Delta p$ such that a previously finite $\Delta p$ value is reduced to exactly zero during a model time step $\Delta t$. Unless the $T$ field is constant, solving (\ref{eqn:transport}) in that case will generally {\em not} result in a zero value of $T \Delta p$, rendering $T$ infinite. Equally challenging is the case where $\Delta p$ falls off to, say, 1\% of its previous value during $\Delta t$. Depending on the magnitude of $\nabla T$ and of the mass fluxes affecting the grid point in question, the product $T \Delta p$ obtained from (\ref{eqn:transport}), as well as the new $T$ value obtained by dividing $T \Delta p$ by $\Delta p$, may end up far outside their physical range in such a case. Our present method for dealing with this problem (which is not caused by roundoff errors but by the time and space discretization of the governing equations) is outlined in Appendix D. Vanishing layer thickness requires special measures not only in grid-scale transport algorithms but also in the formulation of subgrid-scale turbulent mixing terms. Appendix D addresses the issue of mixing as well. {\em 2.14 Time smoothing.} The leap frog time differencing scheme possesses a computational mode, i.e., the tendency to decouple the solution for even and odd time steps. We suppress this mode by applying an Asselin (1972) time filter, denoted hereafter by an overbar, on all prognostic variables. In the case of the layer thickness $\Delta p$, the smoothed value at time level $n$ is obtained from values at three consecutive time levels $n$--1,$n,n$+1 in the following manner: \begin{equation} \overline{\Delta p}^{(n)} = (1-2\gamma) \Delta p^{(n)} + \gamma \left( \overline{\Delta p}^{(n-1)} + \Delta p^{(n+1)} \right). \label{eqn:Asselin1} \end{equation} As shown by Asselin (1972), the weight $\gamma$ may be as large as 0.5, but a much smaller value is sufficient to suppress the computational mode. We presently use $\gamma = 2^{-6}$ for $\Delta p$ and for the thermodynamic variables, and $\gamma = 2^{-3}$ for $u,v$. Note that $\Delta p^{(n-1)}$ enters this formula in time-smoothed form; this eliminates the need to store smoothed and unsmoothed model variables simultaneously. The conservation properties of the time filter require a closer look. Since the leap frog scheme is only used on the ``baroclinic" part of the total $\Delta p$ field whose temporal changes do not affect the bottom pressure (see the discussion of the split-explicit time integration scheme in BS), it is only this part that needs to be subjected to the process indicated in (\ref{eqn:Asselin1}). Now if the sum of the unsmoothed $\Delta p$ values in each grid column is time-independent, one can infer from (\ref{eqn:Asselin1}) that the vertical sum of the $\overline{\Delta p}$ values will be time-independent as well. In other words, the time smoothing process conserves overall mass in the model. Mass transfer among individual layers, on the other hand, is unavoidable and is, in fact, the reason for using as small a $\gamma$ value as possible. In order to render the temperature and salinity time smoothing process conservative, $T$ and $S$ are smoothed according to \begin{equation} \overline{T}^{(n)} = \left( \overline{\Delta p}^{(n)} \right)^{-1} \left[ (1-2\gamma) T^{(n)} \Delta p^{(n)} + \gamma \left( \overline{T}^{(n-1)} \overline{\Delta p}^{(n-1)} + T^{(n+1)} \Delta p^{(n+1)} \right) \right]. \label{eqn:Asselin2} \end{equation} To accomodate the transition $\overline{\Delta p} \rightarrow 0$, (\ref{eqn:Asselin2}) is actually evaluated in the form \begin{equation} \overline{T}^{(n)} = \left( \overline{\Delta p}^{(n)} + \epsilon \right)^{-1} \left[ T^{(n)} \{ (1-2\gamma) \Delta p^{(n)} + \epsilon \} + \gamma \left( \overline{T}^{(n-1)} \overline{\Delta p}^{(n-1)} + T^{(n+1)} \Delta p^{(n+1)} \right) \right] \label{eqn:Asselin3} \end{equation} where $\epsilon$ represents a residual layer thickness presently set to 1\,mm. {\em 2.15 Diapycnal mixing.} At the present stage of model development, mixed-layer entrainment/detrainment and cabbeling (see Sections 2.6 and 2.8) are the only mechanisms capable of exchanging mass and other material properties among coordinate layers. Diapycnal mixing {\em per se} has not yet been implemented in the Atlantic basin version of the model. \newpage {\bf 3. A MECHANICAL-THERMOHALINE SPINUP EXPERIMENT} We proceed now to illustrate the dynamic characteristics of the Miami OGCM in a spinup experiment from an isohaline, zonally symmetric baroclinic initial state. While the surface boundary conditions imposed are somewhat realistic, this computation should not be viewed as an effort to simulate the present Atlantic climate, since no interior diapycnic mixing is imposed, nor are there buffer zones at the zonal boundaries to clamp the deep water stratification to climatological conditions. In addition, the fresh-water flux amplitude imposed at the surface is exaggerated in order to accelerate the thermohaline evolution in the model. It is thus the interplay of thermocline ventilation processes and the salinity field evolution with the wind-driven circulation, and the several associated surface and interior adjustment time scales that are illustrated and discussed. Results are presented for a basic model run of 30-year duration, with emphasis in the discussion on the slow evolution of the system state brought about by a gradual drift in the salinity and associated density fields, and for an experiment to test the model sensitivity to the imposed lateral viscosity. The second experiment uses as initial conditions the state at year 20 in the basic model run, and proceeds with all parameters unchanged except for a 50\% decrease in the lateral viscosity coefficients. {\em 3.1 Model configuration and forcing functions.} In an effort to resolve the thermocline ventilation process in an at least marginally satisfactory manner, we have increased the number of interior layers from 5, the configuration used in BHHK, to 10 covering the density range from 26.0 to 27.8 $\sigma$ units (26.0 - 27.8\,kg\,m$^{-3}$) in steps of 0.2 $\sigma$ units. The relatively coarse $2^0$ horizontal mesh in BHHK, on the other hand, has been retained for the sake of computational economy in this exploratory simulation. The model ocean basin extends from approximately $10^0$N to $60^0$N and is closed off by vertical sidewalls at those latitudes. Thus, the present model provides adequate room for the development of the subtropical and subpolar gyres, but the tropical current system and overflow processes from the Norwegian Sea are excluded. Limiting the ocean basin to the $10^0 - 60^0$ range allows us to use the atmospheric forcing functions developed in BHHK. Aside from the two zonally oriented sidewalls, the domain boundaries conform to the actual Atlantic coastline to the degree to which this is possible in a $2^0$ mesh. The bottom topography, which was taken from SBB, likewise can be called ``realistic" but is in fact too coarse to resolve important features like the Caribbean islands (see Fig.\,1 in SBB). As in BHHK, the initial mass distribution in the isopycnic layers is patterned after the annually and zonally averaged North Atlantic density distribution shown in Fig.\,112 of Levitus (1982). The mixed layer is given an initial uniform depth of 50\,m and a uniform density of $\sigma = 25.8\,\mbox{kg\,m}^{-3}$. Model integration starts at the time of spring equinox. The initial salinity is set to 35.0\,g/kg in all layers. Sea-surface latent and sensible heat fluxes are cast in the form of Newtonian relaxation boundary conditions. As outlined in BHHK, an attempt is made to simultaneously relax toward observed sea-surface temperature (SST) and observed heat fluxes. This double goal is achieved by defining an atmospheric pseudo-temperature field which, if used in the traditional bulk heat-flux formula, yields the desired heat flux {\em if} the model reproduces the correct SST. (We regard the alternative, relaxation toward the observed SST, as unattractive because such forcing renders the sensible and latent heat fluxes proportional, and thus overly sensitive, to differences between observed and model-generated SST. In fact, total compliance of the model with the observed SST would in such a forcing mode be accompanied by zero heat flux across the air-sea interface.) Since solar radiation is represented in the form of a separate, model state-independent (``imposed") flux condition, the pseudo-temperature so defined is considerably lower than the actual air temperature. Like the other forcing functions used in the model -- wind stress, short-wave radiation, and fresh-water flux --, the pseudo-air temperature incorporating the effects of latent and sensible heat fluxes varies with time and latitude but represents zonally averaged conditions. Specific expressions for wind stress, atmospheric pseudo-temperature and short-wave radiation are given in BHHK. We added a constant 12\,W\,m$^{-2}$ to the radiation forcing function to remove a bias necessitated by the idealized basin shape in BHHK. Fresh-water flux or ``$E$--$P$" (evaporation minus precipitation), the forcing function controlling the mixed-layer salinity, is patterned after Schmitt {\em et al.} (1989) but had to be modified considerably to achieve basin-wide salinity conservation. Details are as follows. After zonal averaging, the $E$--$P$ data shown in Fig.\,1 of Schmitt {\em et al.} (1989) can be represented in the latitude band $10^0 - 60^0$N to a fair degree of accuracy by the polynomial \[ b_0 - b_1 (y-y_0) \left[ 1 - b_0 \left( \frac{y-y_0}{15} \right)^4 \right] \] where $y$ is the distance from the equator on a Mercator projection, measured in units of $2^0$ longitude, $y_0$ = 21, $b_0$ = 0.45\,m/yr, and $b_1$ = 0.08\,m/yr. Averaged over the present model domain, this function produces a net freshwater loss of 0.7\,m/yr. In the actual North Atlantic, this deficit is compensated for by lateral fresh-water influx from rivers and adjacent ocean areas. To take river inflow into account, appropriate point sources could be added to the model. Influx from outside oceanic areas presents a greater problem; its effects are simulated in many models by relaxing water mass properties toward prescribed conditions along the artificial sidewalls. None of these refinements have been built into our model at this stage, giving us essentially a choice between tolerating a gradual salinity buildup during the model run and artificially increasing the amount of rainfall. We have opted for the second alternative. In order to allow salinity in the subtropical gyre to reach values in excess of 36\,g/kg within a reasonable length of time, we have chosen to add the extra fresh water to the subpolar gyre. This has been achieved by multiplying the above forcing function by a factor of two and adding fresh water at a rate of 1.4\,m/yr uniformly over the model grid. As shown in Fig.\,1, this procedure yields $E$--$P$ values slightly greater than 1\,m/yr in the subtropics, reasonably close to observed conditions, while changing $E$-$P$ in the extratropics from peak values near -33\,cm/yr to values as extreme as -2.07\,m/yr at $55^0$N. The sign change of $E$--$P$ occurs at $33.6^0$N, well to the south of the gyre boundary. The effect of the anomalous fresh-water influx on the high-latitude mixed-layer depth is one of the items to be discussed later. Ideally, the ``$E$" component of the $E$--$P$ forcing function should match the latent heat flux component of the thermodynamic forcing function. No attempt has been made to ensure consistency in this matter. \newpage {\bf 4. RESULTS AND DISCUSSION} The results from the basic 30-year model run are discussed first with emphasis on the evolution of the mean kinetic energy and its distribution between the mixed layer, the seasonally ventilated layers, and the essentially unventilated deep layers. This is done in the context of the large-scale stratification evolution, and of shifts in ventilation patterns between the layers. We proceed then to consider in some detail the regional aspects of the spinup transient, with emphasis on the evolving salinity contrast between the subtropical and subpolar basins, and on the relative effects of lateral (isopycnic) ventilation and seasonal convection. Finally, a low-viscosity run is described whose results are used to make two distinct points about model performance: (i) the large-scale characteristics of the property fields are not overly sensitive to the resulting computational noise, and (ii) the greater eastward penetration of the separated boundary current, with its enhanced salt advection, appears to strengthen the mode-water formation process. Due to a number of model idealizations concerning domain size, vertical and horizontal grid resolution, and forcing functions, the experiments reported here were not expected to produce a realistic set of climatic equilibrium fields. Our intent, rather, was to provide without undue computational effort a quasi-equilibrated set of motion and property fields which illustrate important aspects of basin-wide density field adjustments, and to demonstrate methods of field characterization tailored to isopycnic coordinate representation. {\em 4.1 Energetics and stratification evolution during the spinup phase.} As pointed out before, the optimal resolution properties of isopycnic coordinate models allow one to simulate stratified flow with a ``relatively" coarse vertical mesh. There is a corollary to this statement which is often overlooked: changes in the model state, particularly in water mass distribution, can be diagnosed from isopycnic model output with relative ease because of the small number of output parameters involved. This is illustrated in Figs.\,2 -- 4 which provide three different perspectives on the evolution of the dynamic and water mass structures in the system. They also illustrate the distinction between a mainly dynamic initial spinup regime and the mainly thermodynamic long-term adjustment phase in its development. Fig.\,2 shows the evolution of the basin-averaged kinetic energy, summed over all layers and split into contributions from the mixed layer, the seasonally ventilated layers 2--7, and the essentially unventilated layers 8-11. This characterization of the layers is supported by Fig.\,3, which shows the evolution of mean thickness anomalies in the individual layers. Since the total basin mass is conserved, the top curve (the mixed-layer mean thickness anomaly) is a mirror image of the sum of the remaining curves (layers 2--11). The seasonal mean thickness cycle in the interior layers is thus due to layer inflation associated with mixed-layer warming (and thinning) in spring and summer. It follows that the trend in the annual minima represents a gradual change in layer thickness outside the regions seasonally influenced by convection. These figures, taken together, illustrate several important aspects of the interplay of the thermohaline and mechanical adjustment processes: {\em a.} A rapid initial buildup of total kinetic energy occurs due to the initial process of geostrophic adjustment to the imposed density field. Over 60\% of the resulting kinetic energy resides in the four deepest layers, but most of this energy decays quasi-exponentially with a half time of approximately 4 years. By year 15 the kinetic energy is essentially confined to the surface and pycnocline regions. This rapid development of wind-driven gyre trapping by the pycnocline is typical of basin-wide spinup experiments that exclude (as this one does) significant bottom water formation and transient eddy activity. {\em b.} In Fig.\,3, during this same period, the seasonal inflation rate of layers 8-10 decreases, while layer 9 has gained mass. The reason for this gradual cessation of the seasonal deep-layer ventilation becomes evident when one considers the salinity evolution in two meridional sections, displayed in Fig.\,4. Under the initial isohaline conditions, cold-season thermal convection at high latitudes can reach quite deeply into the underlying water column. Because of the large (exaggerated with respect to natural conditions) meridional water vapor transport implied by the imposed fresh-water flux condition, an intensive halocline develops over the years in the cyclonically forced northern third of the model domain, until by year 20 the seasonal fluctuations seen in Fig.\,3 are barely noticeable in layer 7 and are absent below it. {\em c.} Layer 6 ($\sigma = 26.8\,\mbox{kg}\,\mbox{m}^{-3}$), as depicted in Figs.\,3 and 4, first loses mass because of the high-latitude convection process, while in the subtropics its depth distribution responds rapidly to the spinup of the subtropical gyre. By year 15, the negative trend in seasonal minimum thickness has come to a halt. Thereafter, layer 6 gradually gains mass as a distinct mode-water structure develops in response to the seasonal convection process. The bottom curve in Fig.\,3 gives the heat content of the model ocean as a function of time. It shows that during the first 15 years the ocean loses heat at a rate comparable to that encountered in BHHK. Instead of approaching a steady state, however, the model ocean is found to gain heat in later years. We surmise that the upturn is related to the freshening of the subpolar seasonal thermocline: as the isopycnic layers making up the seasonal thermocline cool to compensate for the salinity decrease (this adjustment is described in Section 2.8), the atmosphere during the mixed-layer entrainment season comes in contact with a slightly colder ocean from which it can extract less heat. {\em 4.2 Regional aspects of the spinup process.} The late-fall mixed-layer depths indicated in Fig.\,4 show little change between the three decadal states depicted, consistent with a mainly local character of the seasonal thermocline development. In contrast, the depth of wintertime convective penetration, as indicated by the mid-depth pycnocline positions, decreases with time in the subpolar domain and grows somewhat within the subtropical gyre. Note also that the subtropical pycnocline structure is only changing marginally in the third decade of the model integration, in spite of an ongoing drift in the salinity field. The persistent salinity gradient in the subtropical pycnocline is an inevitable consequence of the ongoing differential fresh-water input. In order to give a clear impression of the basinwide structure of these processes, a set of horizontal property maps is provided in Figs.\,5 -- 8. The decadal evolution of the ventilation depth is illustrated in Fig.\,5. Shown there is mixed-layer depth in March, approximately the time of its extreme vertical penetration. The impact of the high-latitude halocline evolution is seen in the disappearance of the layer-depth maximum at the northern domain boundary. The fact that increasing salinity (cf. Fig.\,4) allows the subtropical convection mode to dig deeper into the underlying layers, which still retain their original salinity, is also evident. The corresponding salinity field at the end of the third decade of integration and the mixed-layer circulation pattern are shown in Figs.\,6 and 7, respectively. Note the close association of the salinity front in the western basin with the confluence between the subtropical and subpolar circulations. Since the model is not eddy-resolving, the salinity contrast between the gyres will grow until a combination of lateral diffusion and Ekman drift effects in the frontal region balances the imposed meridional fresh-water flux. {\em 4.3 Mode-water development.} As a consequence of the rapid freshening of the high-latitude surface water, deep wintertime convection becomes concentrated along the southern flank of the separated western boundary current. This region is characterized by salt advection from the south which counteracts the tendency of the fresh-water input to impede production of dense water by wintertime cooling. The resulting tendency toward regionally concentrated mode-water formation is illustrated by the late-winter thickness patterns in layers 5 and 6 ($\sigma$ = 26.6 and 26.8\,kg\,m$^{-3}$) at the beginning and at the end of the third decade (see Fig.\,8). As already suggested in the discussion of Figs.\,3 and 4, a distinct mode-water signature is found mainly in layer 6, but a contribution of similar structure is also evident in layer 5 near the western boundary. Referring to Fig.\,4, we see that these layers occupy approximately the depth range characteristic of the subtropical mode water, the 18-degree water (Worthington, 1959) found in the western Atlantic. Its structure appears stable over the third decade of integration, but a gradual salinity increase in time suppresses the layer-5 contribution, that is, a drift takes place towards higher mean mode-water salinity. We will have an opportunity for further remarks on this trend in the next section, where a low-viscosity perturbation experiment is presented. {\em 4.4 A lowered-viscosity variation of the basic experiment.} The $2^0$ mesh employed in the present experiment is too coarse to resolve barotropically or baroclinically unstable eddies which in reality play a pivotal role in dissipating kinetic and available potential energy. Mixing terms of the form $\nu \nabla^2 u$ in the prognostic equations take over this task in noneddy-resolving ocean models such as the present one. If $\nu$ is chosen too small, the model can prevent unlimited energy growth only by creating its own, distorted version of the relevant hydrodynamic instabilities. Due to coarse grid resolution, these instabilities usually appear on the smallest resolved scale. A relevant question to ask is whether the resulting ``noise" in the model fields will cause numerical instabilities, and whether for this reason alone fine-tuning of $\nu$ is a critical aspect of model development. To shed light on this question, we have run a repeat integration of the final decade of the 30-year basic run, with identical parameters except for a lowering of all lateral mixing coefficients by half. Four parameters are affected: the $u_d$ value used for $T,S$ mixing, the two parameters $u_d, \eta$ appearing in (\ref{eqn:eddyvisc}), and the interface smoothing velocity $u_d$ in (\ref{eqn:thknssflux}). Fig.\,9, which picks up the time history from Fig.\,2 at year 15, shows the effect of this change on the mean kinetic energy and its seasonal variability. Energy levels are generally higher, but the model appears to be numerically stable. Fig.\,10 shows that the gross circulation characteristics remain virtually unchanged, except for a slight enhancement of the recirculation near the western boundary which is also reflected in the surface temperature patterns. This is an expected consequence of lowering the viscosity toward a value where the so-called Munk scale is marginally resolved. [To improve readability, we have removed two-grid interval noise from all horizontal fields by applying a 3-point (0.25--0.5--0.25) smoother in both grid directions. Layer interface traces in the cross sections, on the other hand, are unsmoothed.] The impact of the lowered viscosity on the formation and lateral spreading of mode water is more noticeable, as seen in the late-winter thickness patterns in Figs.\,11 and 12. The two panels in Fig.\,11 correspond to the layer-5 and layer-6 thickness maps shown in the lower part of Fig.\,8, while Fig.\,12 allows a comparison of the layer-7 thickness patterns at year 30 for both runs. The changes are most marked in layer 7 ($\sigma = 27.0\,\mbox{kg}\,\mbox{m}^{-3}$), where onset of mode-water formation is suggested in the longitude range from $40^0$W to $60^0$W. This is accompanied by a thinning of layer 6 in the same region, a shift in water mass distribution clearly seen in the vertical sections along $52^0$W shown in Fig.\,13. The physical process responsible for these changes is the enhanced eastward transport of saline subtropical water by the slightly less diffusive North Atlantic current. We note, finally, that the deep-layer kinetic energy seen in Fig.\,9 is continually increasing during the low-viscosity run, suggesting that near-bottom flows are being driven by the mid-depth adjustments in the density fields. A significant change in the layer-11 ($\sigma = 27.8\,\mbox{kg}\,\mbox{m}^{-3}$) circulation pattern is indeed indicated in Fig.\,14. Apparently a cyclonic circulation has been induced in the deep waters below. The precise interaction mode, if indeed a dominant one exists, is not apparent to us at the moment and thus needs further study. \newpage {\bf 5. SUMMARY} We have documented and tested a three-dimensional primitive-equation model which represents the essentially mature state of an isopycnic-coordinate OGCM if allowance is made for the absence of explicit cross-isopycnal mixing and sea-ice related processes. Special attention is paid to recently developed technical model aspects associated with inter-layer mass exchange due to surface forcing of temperature and salinity fields. The action of these processes is demonstrated in a mechanical-thermodynamic spinup integration from a baroclinic but isohaline initial state. By imposing a fresh-water flux at the surface which represents more than twice the present-day precipitation-evaporation difference, while seasonal winds and thermal boundary conditions were chosen to be representative of the current climate, we have produced within a single model run illustrations of a range of phenomena which can be expected to occur if major climate-related perturbations are imposed on the world ocean system. One must keep in mind that the coarse horizontal mesh excludes eddy-driven horizontal recirculations, and that the still fairly coarse density coordinate interval (corresponding to about $1^0$C in the main thermocline) limits representation of mode-water density changes to a relative thickness adjustment between generally two, and at most three isopycnic layers. It is also clear, on the other hand, that layer-thickness differences of tens of meters are well resolved on the gyre scale. While configured with North Atlantic geometry, the model state remains Atlantic-like only in the subtropical gyre, where the mechanical forcing imposes a density structure (baroclinicity) in the main pycnocline which is essentially equilibrated after two decades of integration. The salt advection by the western boundary current extension assures maintenance of the deep convection in the mode-water formation region in spite of the fact that the fresh-water influx at the surface is positive there. Meridional salinity gradients persist in the ventilated layers as a consequence of the continual salinity evolution. The importance of salt advection in the mode-water formation process is also borne out by the low-viscosity experiment. While the salinity in the isopycnal columns in Fig.\,13 appear at first glance nearly identical, closer inspection reveals a northward as well as downward displacement of the 36\,g/kg isohaline in the low-viscosity case. The sub-polar domain is, in contrast, more Pacific-like, with a pronounced halocline inhibiting deep convection there. While in the present experiment this is due to the negative buoyancy generation by the fresh-water input, it is also likely that without inclusion of the Arctic basin in the model domain the halocline state would not be avoided even with realistic values of the overall fresh-water fluxes. Important aspect of isopycnic OGCM performance are the clear depiction of mixed-layer depth, the ability to distinguish ventilated and unventilated thermocline regions (see BHHK), and the explicit modeling of long-term changes in the production rate of different water masses in response to transients in the atmospheric forcing. The gradual evolution of the wintertime mixed-layer depth seen in Figs.\,4 and 13, and of the layer thickness fields shown in Figs.\,8, 11, and 12, would be difficult to diagnose from mass fields produced by $z$-coordinate models, even with a 10\,m vertical mesh size. There remains a need to develop a deeper understanding of the consequences for model behavior of different approaches to system discretization, particularly in the interaction between the surface mixed layer and the interior layers considered by BHHK, and in lateral boundary layer dynamics. In particular, Chassignet and Gent (1991) have shown that jet separation processes at western basin boundaries behave quite differently in depth- and density-discretized models, at least with coarse vertical resolution. Regarding further model applications, the C-grid formulation was chosen for the horizontal mesh in order to allow straightforward migration to eddy-resolving mesh sizes. The model runs reported here did not include diapycnic mixing effects except in the context of cabbeling, but code for this purpose has been developed and has been tested with good success within the BHHK framework. Results will be reported soon. \vspace{1cm} \underline{Acknowledgements:} This work was supported by the National Science Foundation under Grant No. OCE-8812185. Model integrations were performed on the computers of the National Center for Atmospheric Research in Boulder, Colorado, operating under sponsorship of the National Science Foundation. \newpage \begin{center} APPENDIX A \\ {\bf Restoration of Reference Density in Isopycnic Coordinate Layers} \end{center} We presently make three passes through each grid column to bring the $T,S$ values at isopycnic grid points into compliance with the appropriate reference densities. The first pass, which affects all layers $k \geq 2$ ($k=1$ being the mixed layer) deals with water that is too light, a condition diagnosed by $T_k > T(\rho_k,S_k)$ where the function $T(\rho,S)$ is the inverted equation of state and $T_k, S_k$ are the actual temperature and salinity values found at grid point $\rho_k$. If this inequality holds, heat is redistributed vertically within the layer such that some water of the desired temperature $T(\rho_k,S_k)$ and some warmer water of temperature $T(\rho_{k-1},S_k)$ is created. The latter amount, usually a small fraction of the total mass $\Delta p_k$ in the layer, is given by \[ \epsilon = \Delta p_k \,\frac{T_k-T(\rho_k,S_k)}{T(\rho_{k-1},S_k)-T(\rho_k, S_k)} \] which reflects the requirement that total heat be conserved in the redistribution process. The amount $\epsilon$, in effect a numerical approximation to the cabbeling-related portion of $\dot{s} \partial p/ \partial s) \Delta t$ in (\ref{eqn:continuity}) and (\ref{eqn:transport}), is subsequently transferred to layer $k$--1. The second pass cycles through all isopycnic layers except the lowest one and remedies situations where $T_k < T(\rho_k,S_k)$. It results in the transfer of the amount \[ \epsilon = \Delta p_k \,\frac{T(\rho_k,S_k)-T_k}{T(\rho_k,S_k)-T(\rho_{k+1}, S_k)} \] to layer $k$+1. The third pass proceeds from the bottom on up. Its purpose is to lighten the water in the bottom layer which may end up being too heavy as a result of pass 2. Density reduction is achieved by entraining into the bottom layer an appropriate amount of water from the layer(s) above. Note that the salinity of water transferred out of a layer generally will not match the salinity of the receiving layer, causing additional density discrepancies down the line. The above algorithm should therefore be viewed as an iterative one; we find it necessary to invoke it at every model time step. \begin{center} APPENDIX B \\ {\bf Time Integration of the Continuity Equation} \end{center} The equation for the ``baroclinic" layer thickness tendency $\Delta p'$, reproduced here from BS (Eqn.\,10), is \[ \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'). \] (Note that cabbeling and mixed-layer entrainment/de\-train\-ment terms governing fluxes across layer interfaces are treated elsewhere.) To solve this equation numerically, we write it in the form \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 ) \label{eqn:continuity2} \end{equation} where ${\bf v} = \overline{\bf v} + {\bf v}'$. Time integration is by the Flux Corrected Transport scheme (Zalesak, 1979) which in its present implementation is a 4-step procedure. {\em Step 1.} The $\Delta p'$ field is advanced in time by low-order (diffusive) fluxes involving the total {\bf v} field. ``Low-order" in this case means forward-upstream, which in the C grid translates into flux expressions of the form \[ (u \Delta p')_{i-1/2}^{LO} = \left\{ \begin{array}{lll} u^{(mid)}_{i-1/2} \ \Delta p'^{(old)}_{i-1} & & \mbox{if $u_{i-1/2} > 0$} \\ u^{(mid)}_{i-1/2} \ \Delta p'^{(old)}_i & & \mbox{if $u_{i-1/2} < 0$} \end{array} \right. \] where $mid$ and $old$ refer to time levels in the leap frog scheme. (In writing this formula we have omitted scale factors needed to achieve mass conservation in non-uniform grids.) {\em Step 2.} The $\Delta p'$ field is further modified by so-called antidiffusive (high- minus low-order) fluxes, again formulated using the total velocity $\overline{\bf v} + {\bf v}'$. In a process central to the FCT scheme, the antidiffusive fluxes are clipped individually to avoid producing ``ripples" that are not already present in the low-order solution, and in particular negative $\Delta p'$ values. ``High-order" in the present case means $2^{nd}$-order space- and time-centered, i.e., \[ (u \Delta p')_{i-1/2}^{HI} = u^{(mid)}_{i-1/2} \ \frac{\Delta p'^{(mid)}_{i-1} + \Delta p'^{(mid)}_i}{2}. \] {\em Step 3.} High-order mass fluxes [even those based on the total velocity ${\bf v} = \overline{\bf v} + {\bf v}'$ if used in conjunction with the term on the right-hand side of (\ref{eqn:continuity2})] preserve bottom pressure. Any clipping of antidiffusive fluxes, i.e., multiplication of $({\bf v}^{HI} - {\bf v}^{LO})$ by a factor $\gamma < 1$, will cause bottom pressure changes. Hence, the fluxes $(1-\gamma)({\bf v}^{HI} - {\bf v}^{LO})$ discarded in the clipping process must be accounted for to restore bottom pressure. As discussed in Section 2 of BS, this is done by summing up the discarded fluxes vertically and prorating the sum among individual layers such that the resulting ``secondary" flux corrections do not interfere with the requirement $\Delta p' \geq 0$. Specifically, by setting \[ U = \sum_k (1-\gamma_k)({\bf v}^{HI} - {\bf v}^{LO})_k \] where $k$ is the layer index, we define the secondary flux correction in $x$ direction in layer $k$ to be \[ (u \Delta p')_{i-1/2,k} = \left\{ \begin{array}{lll} \frac{\Delta p'_{i-1,k}}{p_{b,i-1,k}} \,U_{i-1/2} & & \mbox{if $U_{i-1/2} > 0$} \\ \frac{\Delta p'_{i ,k}}{p_{b,i ,k}} \,U_{i-1/2} & & \mbox{if $U_{i-1/2} < 0$} \end{array} \right. \] The values $\Delta p'$ and $p_b$ appearing in this formula are those resulting from steps 1 and 2 above. It is easy to show that the fluxes so defined remove the effect of flux clipping on bottom pressure. {\em Step 4.} The function of the right-hand term in (\ref{eqn:continuity2}) is to keep $\sum \Delta p'_k$ ($k$ being the layer index) constant in each grid column. Since the term is of the form $const. \times \Delta p'_k$, its effect on the solution can be evaluated by multiplying the $\Delta p'$ values arrived at in steps 1 -- 3 by a constant reflecting the deviation of $\sum \Delta p'_k$ from the initially prescribed bottom pressure. Naturally, the constant will vary from grid column to grid column. \begin{center} APPENDIX C \\ {\bf Subgrid-scale mixing} \end{center} If a thickness diffusion term $-\nabla \cdot (\nu \nabla \,\Delta p)$ is added to the continuity equation, one must be careful not to modify the total water depth while diffusing $\Delta p$. In other words, the thickness diffusion term must integrate to zero in each grid column. This is most easily accomplished by smoothing not the thickness field but the pressure on interior coordinate interfaces (a definition which excludes the sea surface and bottom). This operation must be executed such that no two interfaces become intertwined. The details of the interface smoothing procedure are as follows. Let $p^k_i$ and $p^k_{i-1}$ be the pressure on the $k$th interface at two neighboring grid points $x_i, x_{i-1}$ and let $P_i$ and $P_{i-1}$ be the respective bottom pressures. The diffusive interface pressure flux, treated as positive if it leads to an increase of $p^k_i$ or a {\em lowering} of the interface at $x_i$, is then defined as \begin{equation} F^k_{i-1/2} = min \left( P_i-p_i, \,max \left[ p_{i-1}-P_{i-1}, \frac{u_d \Delta t}{\Delta x} (p_{i-1}-p_i) \right] \right). \label{eqn:thknssflux} \end{equation} (In writing this formula we have omitted scale factors needed to achieve mass conservation in unevenly spaced grids.) The thickness change of the layer sandwiched between interfaces $k$ and $k$--1 ($k$ increasing downward) at point $x_i$ during time interval $\Delta t$ is then given by \[ (F^{k-1}_{i+1/2} - F^{k-1}_{i-1/2}) - (F^k_{i+1/2} - F^k_{i-1/2}). \] The layer thickness change caused by interface pressure diffusion implies an isopycnic mass flux $({\bf v} \Delta p)_{diff}$ which for $u_d$ values of order 1\,cm\,s$^{-1}$ may be as important as other advection or diffusion processes in the model. The transport scheme we use to evaluate the effect of $({\bf v} \Delta p)_{diff}$ on the $T,S$ fields is once again the Smolarkiewicz scheme from which we omit the time-consuming antidiffusive steps normally giving the scheme its $3^{rd}$-order accuracy. The $u_d$ value presently used for thickness diffusion is 0.5\,cm\,s$^{-1}$. \begin{center} APPENDIX D \\ {\bf Transport and diffusion in the presence of massless grid points} \end{center} An ad-hoc way to prevent the term $\nabla_s \cdot ({\bf v} \,\Delta p)$ in (\ref{eqn:transport}) from producing unrealistic $T$ values during a model time step is to discard the new $T$ value and retain the old one whenever the new $\Delta p$ falls below some predetermined threshold. However, this approach introduces potentially dangerous bifurcations in the transport algorithm. To assure a gradual transition between nonzero and zero $\Delta p$ values, we replace the common finite-difference expression for the tendency term in (\ref{eqn:transport}), \[ \frac{T^{(new)} \Delta p^{(new)} - T^{(old)} \Delta p^{(old)}}{\Delta t} \] by \begin{equation} \frac{T^{(new)} (\Delta p^{(new)}+\epsilon) - T^{(old)} (\Delta p^{(old)}+\epsilon)}{\Delta t}. \label{eqn:tendency} \end{equation} The parameter $\epsilon$ is designed to be virtually zero unless a grid point loses more than 90\% of its mass during a single time step. If the mass loss exceeds this threshold, $\epsilon$ gradually increases to 10\% of $\Delta p^{(old)}$ as $\Delta p^{(new)}$ approaches zero. Specifically, we set \[ \epsilon = A + (\epsilon^2_1 + A^2)^{1/2} \] where $A = (0.1\,\Delta p^{(old)} - \Delta p^{(new)})/2$. The parameter $\epsilon_1$, which represents a layer thickness of order 1\,cm, assures that (\ref{eqn:tendency}) continues to work if $\Delta p^{(old)}$ and $\Delta p^{(new)}$ both approach zero. To handle turbulent mixing processes in a conservative manner, all eddy diffusion terms in the prognostic model equations must be weighted by $\Delta p$ in the manner shown in (\ref{eqn:diffusionterm}) and on the right-hand side of (\ref{eqn:transport}). To find an expression for $\Delta p$ in flux expressions like $-\nu \Delta p\, \nabla T$ which works well in case of strongly varying layer thickness we proceed as follows. Consider two adjacent grid points $x_i, x_{i-1}$ with associated layer thicknesses $\Delta p_i, \Delta p_{i-1}$. Let us treat $\Delta p$ as constant within a grid box centered on each grid point, i.e., assume that $\Delta p$ jumps from $\Delta p_{i-1}$ to $\Delta p_i$ at $x_{i-1/2}$. Assume the same stepwise behavior for $T$. Turbulent mixing across the interface at $x_{i-1/2}$ is now evaluated in two steps. In step 1, a mixed-water slab of width $2 \eta$ straddling the interface is created by mixing water found within a distance $\eta$ from the interface. The value of $T$ in the slab is \[ \frac{T_i \Delta p_{i-1} + T_{i-1} \Delta p_i} { \Delta p_{i-1} + \Delta p_i}. \] In step 2, the slab is split into the two original sub-slabs of width $\eta$, both of which are then re-absorbed into the grid box they originated from. Our intent is to formulate $-\nu \Delta p \,\nabla T$ such that it produces the same $T$ tendencies at $x_{i-1}$ and $x_i$ as the 2-step process sketched above. It is easy to show that this requires $\Delta p$ in $-\nu \Delta p \,\nabla T$ to be expressed as the {\em harmonic} mean $2/(\Delta p_{i-1}^{-1} + \Delta p_i^{-1})$ of the layer thicknesses involved. This prescription is adopted in all diffusion terms in the model. In order to smoothly approach the limit $\Delta p = 0$ [note that (\ref{eqn:diffusionterm}) requires division by $\Delta p$], we add a small $\epsilon$ to all $\Delta p$ values involved in the mixing term calculation. Information from non-massless grid points is thereby allowed to gradually infiltrate massless grid areas. $\epsilon$ is presently set to represent a residual layer thickness of 1\,mm. \begin{center} APPENDIX E \\ {\bf Detrainment into isopycnic layers in the case of variable temperature and salinity.} \end{center} The ability of the model to support variable $T,S$ fields complicates the mass transfer from the mixed layer to isopycnic coordinate layers, given the need to maintain the reference density in the receiving layers at all times. Our solution to this problem is described below. Let $z_1, T_1, S_1$, respectively, be the mixed-layer thickness, temperature, and salinity, and let $z_k, T_k, S_k$ be the corresponding quantities in the $k$-th isopycnic layer. Suppose $z_1$ exceeds the appropriate mixed-layer depth which in the Kraus-Turner model is given by the Monin-Obukhov length $L$. Our goal is to detrain as large a fraction of the excess amount $z_1-L$ as possible into layer $k$. Generally, this cannot be done without allowing some heat transfer from the detrained column to the remaining mixed layer. Let $\Delta T$ be the largest permissible mixed-layer temperature increase resulting from this transfer. (Our way to arrive at a plausible $\Delta T$ value is to argue that the thermal energy received during the last time step should have been spread over the depth range $L$ rather than $z_1$.) Let $z \ (\leq z_1-L)$ be the amount of mixed-layer water actually being detrained. Since mixed-layer water is more buoyant than water in layer $k$, mixing some of the former with the latter generally produces water too light to match the reference density $\sigma_k$. We correct this by extracting heat from the layer-$k$ mixture which we add to the remaining column $z_1-z$ with the intent of raising its temperature to $T_1+\Delta T$. The unknown in the problem is $z$. If the heat required to raise $T$ in the column $z_1-z$ by $\Delta T$ is taken out of the column $z$ at the beginning, the water to be detrained is cooled down by the amount $\Delta T (z_1-z)/z$. Adding this water to layer $k$ produces water of temperature \[ T_{new} = \frac{ \left( T_1 - \Delta T \frac{z_1-z}{z} \right) z + T_k z_k} {z+z_k} = \frac{(T_1+\Delta T) z - \Delta T z_1 + T_k z_k}{z+z_k} \] and salinity \[ S_{new} = \frac{ S_1 z + S_k z_k}{z+z_k}. \] We demand that \begin{equation} \sigma(T_{new},S_{new}) = \sigma_k. \label{eqn:ref_dens} \end{equation} Approximating the equation of state by a $3^{rd}$-order polynomial \begin{equation} \sigma(T,S) = c_1 + c_2 T + c_3 S + c_4 T^2 + c_5 TS + c_6 T^3 + c_7 T^2 S \label{eqn:state} \end{equation} (Friedrich and Levitus, 1972) and noting that $T_{new},S_{new}$ are of the form $(az+b)/(cz+d)$ and $(ez+f)/(cz+d)$, respectively, we find after some algebraic manipulation that (\ref{eqn:ref_dens}) is satisfied if $z$ is a root of the cubic equation \begin{equation} a_3 z^3 + a_2 z^2 + a_1 z + a_0 = 0 \label{eqn:z-eqn} \end{equation} where \begin{eqnarray*} a_0 & = & d^2[d(c_1-\sigma_k) + b\,c_2 + f\,c_3] + b[df\,c_5 + b(d\,c_4 + b\,c_6 + f\,c_7)] \\ a_1 & = & d[3cd(c_1-\sigma_k) + (2bc+ad)c_2 + (2cf+de)c_3] \\ & & + b[(2ad+bc)c_4 + 3ab\,c_6 + (2af+be)c_7] + [adf + b(de+cf)]c_5 \\ a_2 & = & c[3cd(c_1-\sigma_k) + (2ad+bc)c_2 + (2de+cf)c_3] \\ & & + a[(2bc+ad)c_4 + 3ab\,c_6 + (2be+af)c_7] + [bce + a(cf+de)]c_5 \\ a_3 &= & c^2[c(c_1-\sigma_k) + a\,c_2 + e\,c_3] + a[ce\,c_5 + a(c\,c_4 + a\,c_6 + e\,c_7)]. \end{eqnarray*} At low temperatures, where (\ref{eqn:state}) is most nonlinear, (\ref{eqn:z-eqn}) is occasionally found to have more than one positive root; however, only one of these usually falls into the range $0 \leq z \leq z_1-L$. If several roots lie in that range, the largest one is taken. If there is no positive root, no water is being detrained. If the smallest positive root exceeds $z_1-L$, all water below depth $L$ is detrained into layer $k$, provided the resulting mixed-layer temperature increase (which in this case is known to be less than $\Delta T$) is positive. \newpage \begin{center} REFERENCES \end{center} \parindent 0pt Arakawa, A., and V. R. Lamb, 1977: Computational design of the basic processes of the UCLA general circulation model, {\em Methods Comput. Phys., 17}, 174--265. Asselin, R., 1972: Frequency filter for time integrations. {\em Mon. Wea. Rev., 100}, 487--490. Bleck, R., 1978: Finite difference equations in generalized vertical coordinates. Part I: Total energy conservation. {\em Contrib. Atm. Phys., 51}, 360--372. ---, and D. B. Boudra, 1986: Wind-driven spin-up in eddy-resolving ocean models formulated in isopycnic and isobaric coordinates. {\em J. Geophys. Res., 91C}, 7611--7621. ---, H. P. Hanson, D. Hu, and E. B. Kraus, 1989: Mixed layer-thermocline interaction in a three-dimensional isopycnic coordinate model, {\em J. Phys. Oceanogr., 19}, 1417--1439. ---, and L. Smith, 1990: A wind-driven isopycnic coordinate model of the north and equatorial Atlantic Ocean. 1. Model development and supporting experiments. {\em J. Geophys. Res., 95C}, 3273--3285. Boudra, D. B., R. Bleck, and F. Schott, 1988: A numerical model of instabilities in the Florida Current. {\em J. Mar. Res., 46}, 715--751. Chassignet, E. P., and P. G. Gent, 1991: The influence of boundary conditions on mid-latitude jet separation in ocean numerical models. {\em J. Phys. Oceanogr., 21}, 1290--1299. Friedrich, H., and S. Levitus, 1972: An approximation to the equation of state for sea water, suitable for numerical ocean models. {\em J. Phys. Oceanogr., 2}, 514--517. Kraus, E. B., and J. S. Turner, 1967: A one-dimensional model of the seasonal thermocline: II. The general theory and its consequences. {\em Tellus, 19}, 98--106. Levitus, S., 1982: {\em Climatological Atlas of the World Ocean}. NOAA Professional Paper 13, 173 pp. Oberhuber, J., 1990: {\em Simulation of the Atlantic circulation with a coupled sea ice -- mixed layer -- isopycnal general circulation model.} Rep. No. 59, Max-Planck-Institut f\"ur Meteorologie, Hamburg, 86 pp. Schmitt, R. W., P. S. Bogden, and C. E. Dorman, 1989: Evaporation minus precipitation and density fluxes for the North Atlantic. {\em J. Phys. Oceanogr., 19}, 1208--1221. Smagorinsky, J. S., 1963: General circulation experiments with the primitive equations. I: The basic experiment. {\em Mon. Wea. Rev., 91}, 99--164. Smith, L. T., D. B. Boudra, and R. Bleck, 1990: A wind-driven isopycnic coordinate model of the north and equatorial Atlantic Ocean. 2. The Atlantic basin experiments. {\em J. Geophys. Res., 95C}, 13105--13128. Smolarkiewicz, P. K., and T. L. Clark 1986: The multidimensional positive definite advection transport algorithm: Further development and applications. {\em J. Comput. Phys., 67}, 396--438. ---, and W. W. Grabowski, 1990: The multidimensional positive definite advection transport algorithm: Nonoscillatory option. {\em J. Comput. Phys., 86}, 355--375. Worthington, L. V., 1959: The $18^0$ water in the Sargasso Sea. {\em Deep-Sea Research, 5}, 297--305. Zalesak, S., 1979: Fully multidimensional flux-corrected transport algorithms for fluids. {\em J. Comput. Phys., 31}, 335--362. \newpage \begin{center} FIGURE CAPTIONS \end{center} Fig.\,1. Evaporation minus precipitation (m/yr) versus latitude (deg). Solid: polynomial fit of Schmitt {\em et al.} (1989) data; dashed: modified curve yielding zero net $E$--$P$ flux over model domain. Fig.\,2. Basin-averaged kinetic energy ($10^2$J\,m$^{-2}$) as a function of time (yrs). Top curve: all layers; second curve: layers 2--7; third curve: mixed layer; bottom curve: layers 8--11. Fig.\,3. Top 11 curves: Basin-averaged thickness (m) of individual model layers as a function of time (yrs), plotted as deviations from long-term average. Bottom curve: deviation of total model heat content ($10^7$J\,m$^{-2}$) from long-term average. Fig.\,4. Vertical-meridional section through model grid along $52^0$W (top) and $36^0$W (bottom) in December of year 10 (left), year 20 (center) and year 30 (right). Solid lines: layer interfaces. Dashed: isopleths of salinity (g/kg). Coordinate labels: longitude and depth. Fig.\,5. Mixed-layer depth in March of year 10 (top), year 20 (center), and year 30 (bottom). Coordinates labeled in deg. latitude/longitude. Fig.\,6. Mixed layer salinity (g/kg) in March, year 30. Fig.\,7. Mass flux streamfunction for the late-winter mixed layer, year 30. Contour interval $2 \times 10^6 \mbox{m}^3 \mbox{s}^{-1}$. Fig.\,8. Thickness (m) of coordinate layers $\sigma$ = 26.6 kg\,m$^{-3}$ (left) and $\sigma$ = 26.8 kg\,m$^{-3}$ (right) in March of year 20 (top) and of year 30 (bottom). Fig.\,9. As in Fig.\,2, but with all lateral mixing coefficients reduced by 50\% at year 20. Fig.\,10. {\em Left:} Barotropic mass flux streamfunction in March, year 30. Contour interval $5 \times 10^6\mbox{m}^3\mbox{s}^{-1}$. {\em Right:} Mixed-layer temperature ($^0$C). Control run (top) and low-viscosity run (bottom). Fig.\,11. Thickness (m) of coordinate layers $\sigma$ = 26.6 kg\,m$^{-3}$ (left) and $\sigma$ = 26.8 kg\,m$^{-3}$ (right) in March of year 30, low-viscosity run. Fig.\,12. Thickness (m) of coordinate layer $\sigma$ = 27.0 kg\,m$^{-3}$ in March of year 30. Control run (left) and low-viscosity run (right). Fig.\,13. Vertical-meridional section through model grid along $52^0$W in March of year 30. Control run (top); low-viscosity run (bottom). Solid lines: layer interfaces. Dashed: isopleths of salinity (g/kg). Fig.\,14. Mass flux streamfunction for the bottom layer ($\sigma$ = 27.8 kg\,m$^{-3}$), March of year 30. Contour interval $0.5 \times 10^6\mbox{m}^3\mbox{s}^{-1}$. Control run (top) and low-viscosity run (bottom). \end{document}