1. Cover

The Parallel Ocean Program (POP) Reference Manual

Ocean Component of the Community Climate System Model (CCSM) and Community Earth System Model (CESM) [1]

R. Smith [2], P. Jones [2], B. Briegleb [3], F. Bryan [3], G. Danabasoglu [3], J. Dennis [3], J. Dukowicz [2], C. Eden [4], B. Fox-Kemper [5], P. Gent [3], M. Hecht [2], S. Jayne [6], M. Jochum [7], W. Large [3], K. Lindsay [3], M. Maltrud [2], N. Norton [2], S. Peacock [3], M. Vertenstein [3], S. Yeager [3]
[1]CCSM is used throughout this document, though applies to both model names
[2](1, 2, 3, 4, 5, 6) Los Alamos National Laboratory (LANL)
[3](1, 2, 3, 4, 5, 6, 7, 8, 9, 10) National Center for Atmospheric Research (NCAR)
[4]IFM-GEOMAR, Univ. Kiel
[5]Brown University, Rhode Island
[6]Woods Hole Oceanographic Institution (WHOI)
[7]Neils Bohr Institute, University of Copenhagen

2. Introduction

This manual describes the details of the numerical methods and discretization used in the Parallel Ocean Program (POP), a level-coordinate ocean general circulation model that solves the three-dimensional primitive equations for ocean dynamics. It is designed for users who want an in-depth description of the numerics in the model, including details of the computational grid, the space and time discretization of the hydrodynamical core and subgridscale parameterizations, and other features of the model.

The POP model is a descendant of the Bryan-Cox-Semtner class of models (Semtner 1986). In the early 1990’s it was written for the Connection Machine by Rick Smith, John Dukowicz and Bob Malone (Smith, Dukowicz, and Malone 1992). An implicit free surface formulation and other numerical improvements were added by Dukowicz and Smith (1994). Later, the capability for general orthogonal coordinates for the horizontal mesh was implemented (Smith, Kortas, and Meltz 1995).

Since then many new features and physics packages have been added by various people, including the authors of this document. In addition, the software infrastructure has continued to evolve to run on all high performance architectures. In 2001, POP was officially adopted as the ocean component of the CCSM based at NCAR. Substantial effort at both LANL and NCAR has gone into adding various features to meet the needs of the CCSM coupled model. This manual includes descriptions of several of the features and options used in the POP 2.0.1 release (Jan 2004), the ocean model configuration of the Spring 2002 release of CCSM2, the Spring 2004 release of CCSM3, and the 2010 release of CCSM4.

Acknowledgments: The work at LANL has been supported by the U.S. Department of Energy’s Office of Science through the CHAMMP Program and later the Climate Change Prediction Program. The work at NCAR has been supported by the National Science Foundation.

Comments, typos or other errors should be reported using the official POP web site at http://climate.lanl.gov

3. The Primitive Equations in General Coordinates

Ocean dynamics are described by the 3-D primitive equations for a thin stratified fluid using the hydrostatic and Boussinesq approximations. Before deriving the equations in general coordinates, we first present, as a reference point, the continuous equations in spherical polar coordinates with vertical z -coordinate. These are standard in Bryan-Cox models (Semtner 1986; Pacanowski and Griffies 2000).

Momentum equations:

(1)\frac{\partial}{\partial t}u + {\cal L} (u) - (uv\tan \phi)/a - fv =
-\frac{1}{\rho_{_0} a\cos \phi} \frac{\partial p}{\partial\lambda} +
{\cal F}_{Hx} (u,v) + {\cal F}_V (u)

(2)\frac{\partial}{\partial t}v + {\cal L} (v) + (u^2\tan \phi)/a + fu =
-\frac{1}{\rho_{_0}a} \frac{\partial p}{\partial\phi} + {\cal F}_{Hy}
(u,v) + {\cal F}_V (v)

(3){\cal L}(\alpha) = \frac{1}{a\cos\phi} \left[\frac{\partial}{
\partial\lambda} (u\alpha) + \frac{\partial}{\partial\phi}
(\cos\phi v\alpha)\right] + \frac{\partial}{\partial z} (w\alpha)

(4){\cal F}_{Hx} (u,v) = A_M\left\{\nabla^2 u + u(1 - \tan^2\phi)/a^2 -
\frac{2\sin\phi
}{a^2\cos^2\phi} \frac{\partial v }{\partial\lambda}\right\}

(5){\cal F}_{Hy} (u,v) = A_M\left\{\nabla^2v + v (1 - \tan^2\phi)/a^2 +
\frac{2\sin\phi
}{a^2\cos^2\phi} \frac{\partial u}{\partial\lambda}\right\}

(6)\nabla^2\alpha = \frac{1}{a^2\cos^2\phi} \frac{\partial^2 \alpha}{
\partial\lambda^2} +
\frac{1}{a^2\cos\phi} \frac{\partial}{\partial\phi}
\left(\cos\phi \frac{\partial \alpha}{\partial\phi}\right)

(7){\cal F}_V (\alpha) = \frac{\partial}{\partial z}\mu \frac{\partial
}{\partial z}\alpha

Continuity equation:

(8){\cal L} (1) = 0

Hydrostatic equation:

(9)\frac{\partial p}{\partial z} = -\rho g

Equation of state:

(10)\rho = \rho (\Theta, S, p)  \rightarrow  \rho (\Theta, S, z)

Tracer transport:

(11)\frac{\partial}{\partial t} \varphi + {\cal L} (\varphi)
= {\cal D}_H (\varphi) + {\cal D}_V (\varphi)

(12){\cal D}_H (\varphi) = A_H \nabla^2 \varphi

(13){\cal D}_V (\varphi) = \frac{\partial}{\partial z} \kappa \frac{\partial
}{\partial z} \varphi,

where \lambda, \phi, z=r-a are longitude, latitude, and depth relative to mean sea level r=a; g is the acceleration due to gravity, f = 2\Omega\sin\phi is the Coriolis parameter, and \rho_{_0} is the background density of seawater. The prognostic variables in these equations are the eastward and northward velocity components (u, v), the vertical velocity w, the pressure p, the density \rho, and the potential temperature \Theta and salinity S. In (11) \varphi represents \Theta, S or a passive tracer. The pressure dependence of the equation of state is usually approximated to be a function of depth only (see Sec. Equation of State). A_H and A_M are the coefficients (here assumed to be spatially constant) for horizontal diffusion and viscosity, respectively, and \kappa and \mu are the corresponding vertical mixing coefficients which typically depend on the local state and mixing parameterization (see Chapter ref:`ch-subgrid). The third terms on the left-hand side in (1) and (2) are metric terms due to the convective derivatives in d{\bf u}/dt acting on the unit vectors in the \lambda, \phi directions, and the second and third terms in brackets in (4) and (5) ensure that no stresses are generated due to solid-body rotation (Williams 1972). The forcing terms due to wind stress and heat and fresh water fluxes are applied as surface boundary conditions to the friction and diffusive terms {\cal F}_V and {\cal D}_V. The bottom and lateral boundary conditions applied in POP (and in most other Bryan-Cox models) are no-flux for tracers (zero tracer gradient normal to boundaries) and no-slip for velocities (both components of velocity zero on bottom and lateral boundaries).

To derive the primitive equations in general coordinates, consider the transformation from Cartesian coordinates (\xi_1, \xi_2, \xi_3 with origin at the center of the Earth) to general horizontal coordinates (q_x, q_y, z), where q_x and q_y are arbitrary curvilinear coordinates in the horizontal directions, and z=r-a is again the vertical coordinate normal to the surface of the sphere. The actual distances along the curvilinear coordinates are denoted x and y, which typically lie along the circumpolar (longitude-like) and azimuthal (latitude-like) coordinate lines, respectively, on dipole grids with two arbitrarily located poles (see Sec. Global Orthogonal Grids). The differential length element ds is given by

(14)ds^2 = d\xi^2_1 + d\xi^2_2 + d\xi^2_3
= h_{ij}^2 dq_i dq_j + dz^2

(15)h_{ij}^2 = \frac{\partial \xi_k}{\partial q_i}  \frac{\partial\xi_k}{\partial q_j},

(where i,j = x,y and repeated indices are summed). The metric coefficients h_{ij} depend on the local curvature of the coordinates. Differential lengths in the z direction are assumed independent of x and y, so no metric coefficients involving z appear. Further restricting ourselves to orthogonal grids, the cross terms vanish, and we have

(16)\begin{aligned}
h_i \equiv h_{ii},& \quad h_{ij} = 0 (i\neq j),& ds_i=h_idq_i.
\end{aligned}

For the purpose of constructing horizontal finite-difference operators corresponding to the various terms in the primitive equations, define:

\Delta_i \equiv ds_i \nonumber

(17)\delta_i \equiv \frac{\partial}{\partial s_i} =
\frac{1}{h_i}\frac{\partial}{\partial q_i},

where \Delta_i and \delta_i can be interpreted as either infinitesimal or finite differences and their associated derivatives. Formulas for the basic horizontal operators (gradient, divergence, curl) can be found in standard textbooks (Arfken 1970). The gradient is

(18)\begin{aligned}
\nabla\psi &=& \hat{\bf x} \frac{1}{h_x}\frac{\partial\psi}{\partial
q_x} +
\hat{\bf y} \frac{1}{h_y}  \frac{\partial\psi}{\partial q_y}
\nonumber \\
&=& \hat{\bf x} \delta_x \psi + \hat{\bf y}\delta_y \psi.
\end{aligned}

where \hat{\bf x} and \hat{\bf y} are unit vectors in the x,y directions. The horizontal divergence is:

(19)\begin{aligned}
\nabla \cdot {\bf u} &=& \frac{1}{h_xh_y} \frac{\partial}{\partial q_x}
(h_y u_x) + \frac{1}{h_xh_y} \frac{\partial}{\partial q_y} (h_xu_y)
\nonumber \\
&=& \frac{1}{\Delta_y} \delta_x (\Delta_y u_x) + \frac{1}{\Delta_x}
\delta_y (\Delta_x u_y),
\end{aligned}

where u_x and u_y are the velocity components along the x and y directions. The advection operator (3) is similar:

(20){\cal L}(\alpha) = \frac{1}{\Delta_y} \delta_x (\Delta_y u_x \alpha)
+ \frac{1}{\Delta_x}
\delta_y (\Delta_x u_y \alpha) + \delta_z (w\alpha).

The vertical component of the curl operator is

(21)\begin{aligned}
\hat{\bf z} \cdot \nabla \times {\bf u} &=& \frac{1}{h_xh_y}
\frac{\partial}{\partial q_x} (h_yu_y) - \frac{1}{h_xh_y}
\frac{\partial}{\partial q_y} (h_xu_x)
\nonumber \\
&=& \frac{1}{\Delta_y} \delta_x (\Delta_y u_y) - \frac{1}{\Delta_x}
\delta_y (\Delta_x u_x).
\end{aligned}

Laplacian-type operators, which appear in the viscous and diffusive terms, have the form

(22)\begin{aligned}
\nabla \cdot G\nabla\psi
= \frac{1}{\Delta_y} \delta_x (\Delta_y G\delta_x \psi) +
\frac{1}{\Delta_x} \delta_y (\Delta_x G\delta_y \psi).
\end{aligned}

where G is an arbitrary scalar function of x and y. Note that these operators have been expressed in terms of the differences and derivatives \Delta_i and \delta_i, and hence there is no explicit dependence on the new coordinates q_i or the metric coefficients h_i. In the discrete operators the same is true: it is not necessary to have an analytic transformation with metric coefficients describing the new coordinate system, it is only necessary to know the location of the discrete grid points and the distances between neighboring grid points along the coordinate directions.

The other horizontal finite difference operators appearing in the primitive equations can also be derived in general coordinates. The Coriolis terms are simply given by

(23)2{\Omega} \times {\bf u} = -\hat{\bf x} fu_y + \hat{\bf y} fu_x.

The metric momentum advection terms corresponding to the third terms on the left in (1) and (2) are given by Haltiner and Williams (1980) (p. 442):

(24)(uv\tan\phi)/a \rightarrow u_xu_yk_y - u^2_y k_x

(25)(u^2\tan\phi)/a \rightarrow u_xu_yk_x - u^2_x k_y

(26)k_x \equiv \frac{1}{h_xh_y} \frac{\partial}{\partial q_x} h_y =
\frac{1}{\Delta_y} \delta_x \Delta_y

(27)k_y \equiv \frac{1}{h_xh_y} \frac{\partial}{\partial q_y} h_x =
\frac{1}{\Delta_x} \delta_y \Delta_x

Note that these revert to the standard forms (left of arrows) in spherical polar coordinates, where h_x = a\cos\phi, h_y = a, u_x = u and u_y = v. The metric terms in the viscous operators (second and third terms on the right in (4) and (5) require a more careful treatment. These terms were derived by Williams (1972) in spherical coordinates, by applying the thin-shell approximation (r\to a) to the viscous terms expressed as the divergence of a stress tensor whose components are linearly proportional to the components of the rate-of-strain tensor. This form is transversely isotropic and ensures that for solid rotation the fluid is stress-free. The general coordinate versions of these terms are derived in Smith, Kortas, and Meltz (1995). The results are

(28)\begin{aligned}
{\cal F}_{Hx} (u_x,u_y) &=& A_M \big\{\nabla^2 u_x - u_x (\delta_x k_x +
\delta_y k_y + 2 k^2_x + 2k^2_y)
\nonumber \\
&&+ u_y (\delta_x k_y - \delta_y k_x) + 2 k_y (\delta_x u_y) -
2 k_x (\delta_y u_y)\big\}
\end{aligned}

The formula for {\cal F}_{Hy}(u_x,u_y) is the same with x and y interchanged everywhere on the r.h.s. It is straightforward to show that these also reduce to the correct form in the spherical polar limit (4), (5). The above forms assume a spatially constant viscosity A_M. More terms appear if A_M is allowed to vary spatially. Wajsowicz (1993) derives the extra terms for spherical polar coordinates. In general orthogonal coordinates they take the form:

(29)\begin{aligned}
{\cal F}_{Hx} (u_x,u_y) &=& \nabla\cdot A_M\nabla u_x
- u_x\big\{\delta_x A_M k_x +
\delta_y A_M k_y + 2 A_M (k^2_x + k^2_y)\big\}
\nonumber \\
&&+ u_y (\delta_x A_M k_y - \delta_y A_M k_x) \\
&&+ (2A_M k_y +\delta_y A_M) (\delta_x u_y) -
(2A_M k_x  +\delta_x A_M)(\delta_y u_y). \nonumber
\end{aligned}

The formula for {\cal F}_{Hy}(u_x,u_y) is again the same with x and y interchanged everywhere.

The general coordinate forms of the anisotropic and biharmonic viscous operators are given in Sec. Horizontal Viscosity below, and the Gent-McWilliams and biharmonic forms of the horizontal diffusion terms are given in Sec. ”ref:sec-horiz-diff.

4. Spatial Discretization

4.1. Discrete Horizontal and Vertical Grids

The placement of model variables on the horizontal B-grid is illustrated in Figure3.1. The solid lines enclose a “T-cell” and the hatched lines enclose a “U-cell”. Scalars (T,S,p,\rho) are located at “T-points” (circles) at the centers of T-cells, and horizontal vectors (u_x,u_y) are located at “U-points” (diamonds) at the corners of T-cells. The indexing for points (i,j) in the logically-rectangular 2-D horizontal grid is such that i increases in the x direction (eastward for spherical polar coordinates), and j increases in the y direction (northward for spherical polar coordinates). A U-point with logical indices (i,j) lies to the upper right (\sim northeast) of the T-point with the same indices (i,j). The index for the vertical dimension k increases with depth, although the vertical coordinate z, measured from the mean surface level z=0, decreases with depth.

../_images/FigHorzStagger.jpg

Figure 3.1: The staggered horizontal B-grid. The x-coordinate grid index i increases to the right (generally eastward), and the y-coordinate index j increases upward (generally northward). Solid lines enclose a T-cell, hatched lines a U-cell. The quantities labeled HTN, HTE, HUW, HUS, as well as the model prognostic variables (\Theta,S,p,u_x,u_y) at the locations shown all have grid indices (i,j).

When the horizontal grid is generated, the latitude and longitude of each U-point and the distances HTN and HTE (see Figure3.1 ) along the coordinates between adjacent U-points are first constructed. Then the latitude and longitude of T-points are computed as the straight average of the latitude and longitude of the four surrounding U-points, and the along-coordinate distances HUW (HUS) between adjacent T-points are computed as the straight average of the four surrounding values of HTE (HTN). Thus T-points are located exactly in the middle of the T-cell, but because the grid spacing in either direction may be non-uniform, the U-points are not located exactly in the middle of the U-cell.

In addition to the grid spacings HTN, HTE, HUS, HUW, several other lengths and areas are also used in the code. These are defined as follows (see Figure3.1 ):

(30)\begin{aligned}
\mathrm{DXU}_{i,j}&=&[\mathrm{HTN}_{i,j}+\mathrm{HTN}_{i+1,j}]/2
\nonumber \\
\mathrm{DYU}_{i,j}&=&[\mathrm{HTE}_{i,j}+\mathrm{HTE}_{i,j+1}]/2
\nonumber \\
\mathrm{DXT}_{i,j}&=&[\mathrm{HTN}_{i,j}+\mathrm{HTN}_{i,j-1}]/2
\nonumber \\
\mathrm{DYT}_{i,j}&=&[\mathrm{HTE}_{i,j}+\mathrm{HTE}_{i-1,j}]/2
\nonumber \\
\mathrm{UAREA}_{i,j}&=&\mathrm{DXU}_{i,j}\mathrm{DYU}_{i,j}
\nonumber \\
\mathrm{TAREA}_{i,j}&=&\mathrm{DXT}_{i,j}\mathrm{DYT}_{i,j}
\end{aligned}

DXU and DYU are the grid lengths centered on U-points, and DXT and DYT are centered on T-points. TAREA and UAREA are the horizontal areas of T-cells and U-cells, respectively.

The construction of the semi-analytic dipole and tripole grids commonly used in POP is described in detail in Smith, Kortas, and Meltz (1995), and is briefly reviewed in Chapter Grids. These grids are based on an underlying orthogonal curvilinear coordinate system with the one or two singularities in the northern hemisphere displaced into land masses (typically North America, Asia, or Greenland). The equator is usually retained as a grid line and the southern hemisphere is a standard spherical polar grid with the southern grid pole located at the true South Pole. These grids are topologically equivalent to a cylinder (periodic in x but not in y), and therefore can be mapped onto a logically-rectangular 2-D array (i,j) which is cyclic in i. Tripole grids require additional communication along the northern boundary of the grid in order to “sew up” the grid along the line between the two northern grid poles. The grids are constructed off line and a file is generated which contains the following 2-D fields: ULAT,ULONG,HTN,HTE,HUS,HUW,ANGLE, where \mathrm{ULAT}_{i,j} and \mathrm{ULONG}_{i,j} are the true latitude and longitude of U-points, and \mathrm{ANGLE}_{i,j} is the angle between the x-direction and true east at the U-point (i,j).

In addition to constructing the grids, the fields ULAT, ULONG and ANGLE are used to interpolate the wind stress fields from a latitude-longitude grid to the model grid if needed. ULAT is also used to compute the Coriolis parameter f at each model grid point.

Figure 3.2 is a diagram of the full three-dimensional T-cell, showing the location of the vertical velocities w and w^U, which advect tracers and momentum, respectively. Note that the vertical velocities w are located in the middle of the top and bottom faces of the T-cell, while the horizontal velocities are located at the midpoints of the vertical edges.

../_images/FigTcell3d.jpg

Figure 3.2: The 3-D T-cell, showing the location of the vertical velocity in T-columns and U-columns.

Since POP is a z-level model, the depth of each point (i,j,k) is independent of its horizontal location (unless partial bottom cells are used, see Sec. Partial Bottom Cells). The vertical discretization is illustrated in Figure3.3. The discrete index k increases from the surface (k=1) to the deepest level (k=km). The thickness of cells at level k is dz_k. T-points are located exactly in the middle of each level, but because the vertical grid may be non-uniform (dz_k\ne dz_{k+1}), the interfaces where the vertical velocities w lie are not exactly halfway between the T-points. The vertical distances between T-points dzw_k are just the average dzw_k=0.5(dz_k+dz_{k+1}), except at the surface where dzw_0=0.5dz_1. The depth of a T-point at level k is zt_k, and zw_k is the depth of the bottom of cells at level k. Note that while the coordinate z is positive upward, zt_k and zw_k are positive depths. Vertical profiles of dz_k are usually generated offline and read in by the code, but there is an option for generating this profile internally. Usually dz_k is small in the upper ocean and increases with depth according to a smooth analytic function describing the thickness as a function of depth. This is necessary in order to maintain the formal second-order accuracy of the vertical discretization; if the vertical spacing changes suddenly the scheme reverts to first order accuracy (Smith, Kortas, and Meltz 1995).

../_images/FigVertGrid.jpg

Figure 3.3: The vertical grid.

The topography is defined in the T-cells, which are completely filled with either land or ocean (except when optional partial bottom cells are used; see Sec. Partial Bottom Cells). Thus U-points lie exactly on the lateral boundaries between land and ocean, and w points lie exactly on the ocean floor. These boundary velocities are by default set to zero for no-slip boundary conditions, though can be modified by parameterizations like the overflow parameterization (Sec. Overflow Parameterization). Vertical velocities w^U along the rims in the stair-step topography may also be nonzero (see the discussion of velocity boundary conditions in Sec. Momentum Advection). The topography is determined by the 2-D integer field \mathrm{KMT}_{i,j} which gives the number of open ocean points in each vertical column of T-cells. The KMT field is usually generated offline and read in from a file in the code. Thus 0\le \mathrm{KMT}\le km, and \mathrm{KMT}=0 indicates a surface land point. In some situations the ocean depth in a column of U-points is needed, and this is defined by the field \mathrm{KMU}_{i,j}, which is just the minimum of the four surrounding values of KMT:

\mathrm{KMU}_{i,j} = \mathrm{min}\{\mathrm{KMT}_{i,j},
\mathrm{KMT}_{i-1,j},\mathrm{KMT}_{i,j-1},\mathrm{KMT}_{i-1,j-1}\}

The depths of columns of ocean T-points and U-points are given, respectively, by:

(31)\begin{aligned}
H_T&=& zw(\mathrm{KMT})\;,\nonumber\\
H_U &=& zw(\mathrm{KMU})\;.
\end{aligned}

With partial bottom cells the depth of the deepest ocean cell in each column has variable thickness, and the above formulas are modified accordingly (see Sec. Partial Bottom Cells).

4.2. Finite-difference operators

The exact finite-difference versions of the differential operators can be easily derived for the various types of staggered horizontal grids A,B,C,D,E (Arakawa and Lamb 1977) given only the forms of the fundamental operators: divergence, gradient, and curl for that type of mesh. POP employs a B-grid (scalars at cell centers, vectors at cell corners) while some OGCM’s use a C-grid (scalars at cell centers, vector components normal to cell faces). We will use standard notation (Semtner 1986) for finite-difference derivatives and averages:

(32)\begin{aligned}
\delta_x\psi &=& \left[\psi\left(x + \Delta_x/2\right) -
\psi\left(x - \Delta_x/2\right)\right]\big/\Delta_x
\\
\overline{\psi}^x &=& \left[\psi\left(x + \Delta_x/2\right) +
\psi \left(x - \Delta_x/2\right)\right]\big/2\;,
\end{aligned}

with similar definitions for differences and averages in the y and z directions. These formulas strictly apply for uniform grid spacing; where, for example, if \psi is a tracer located at T-points, then \psi\left(x + \Delta_x/2\right) is located on the east face of a T-cell. For nonuniform grid spacing, the above definitions should be interpreted such that variables lie exactly at T- or U-cell centers and faces, as appropriate.

The fundamental operators on C-grids have the same form as (18) -(22). On B-grids the derivatives involve transverse averaging, and the fundamental operators are given by:

(33)\begin{aligned}
\nabla\psi          & =  \hat{\bf x} \delta_x \overline\psi^y + \hat{\bf y}  \delta_y \overline\psi^x  \\
\nabla \cdot{\bf u} & =  \frac{1}{\Delta_y} \delta_x \overline{\Delta_y u_x}^y + \frac{1}{\Delta_x} \delta_y \overline{\Delta_x u_y}^x \\
\hat{\bf z} \cdot \nabla \times {\bf u} & = \frac{1}{\Delta_y} \delta_x \overline{\Delta_y u_y}^y - \frac{1}{\Delta_x} \delta_y \overline{\Delta_x u_x}^x \\
\nabla\cdot G\nabla \psi                & = \frac{1}{\Delta_y} \delta_x  \overline{\big[\Delta_y G\delta_x \overline{\psi}^y \big]}^y
                                          + \frac{1}{\Delta_x} \delta_y  \overline{\big[\Delta_x G\delta_y \overline{\psi}^x \big]}^x\;.
\end{aligned}

The gradient is located at U-points and the divergence, curl and Laplacian are located at T-points. In the Laplacian operator G must also be defined at U-points. The factors \Delta_x, \Delta_y inside the difference operators \delta_x, \delta_y are located at U-points and are given by DXU, DYU, respectively, while the factors 1/\Delta_x, 1/\Delta_y outside the difference operators, as well as similar factors in the denominators of the difference operators \delta_x, \delta_y, are evaluated at T-points. For example, the first term on the r.h.s. of the divergence (33) at the T-point (i,j) is given by

\begin{aligned}
&&0.5[\mathrm{DYU}_{i,j}(u_x)_{i,j}
+\mathrm{DYU}_{i,j-1}(u_x)_{i,j-1} \nonumber \\
&&-\mathrm{DYU}_{i-1,j}(u_x)_{i-1,j}
-\mathrm{DYU}_{i-1,j-1}(u_x)_{i-1,j-1}]/\mathrm{TAREA}_{i,j}\end{aligned}

In POP (and in other Bryan-Cox models which use a B-grid formulation) all viscous and diffusive terms are given in terms of an approximate C-grid discretization in order to ensure they will damp checkerboard oscillations on the scale of the grid spacing (see Secs. Horizontal Tracer Diffusion, Horizontal Tracer Diffusion, and Horizontal Viscosity).

4.3. Discrete Tracer Transport Equations

The discrete tracer transport equations are:

(34)\frac{\partial}{\partial t}(1+\xi) \varphi + {\cal L}_{T} (\varphi)
= {\cal D}_H (\varphi) + {\cal D}_V (\varphi) + {\cal F}_{_W}(\varphi)

where {\cal L}_T is the advection operator in T-cells, and {\cal D}_H, {\cal D}_V are the horizontal and vertical diffusion operators, respectively. The factor (1+\xi) is associated with the change in volume of the surface layer due to undulations of the free surface, and {\cal F}_{_W}(\varphi) is the change in tracer concentration associated with the freshwater flux. \xi and {\cal F}_{_W}(\varphi) are given by

(35)\xi = \frac{\delta_{k1}}{dz_1}\,\eta

(36){\cal F}_{_W}(\varphi) = \frac{\delta_{k1}}{dz_1}\,
\sum_{m} q^{(m)}_{_{W}}\varphi^{(m)}_{_W}

where \delta_{k1} is the Kronecker delta, equal to 1 for k=1 and zero otherwise. \eta is the displacement of the free surface relative to z=0, q^{(m)}_{_{W}} is the freshwater flux per unit area associated with a specific source (labeled m) of freshwater. Thus, q{_{W}}=\sum_{m} q^{(m)}_{_{W}} = P - E + R - F_{ice} + M_{ice} is the total freshwater flux per unit area associated with precipitation P, evaporation E, river runoff R, freezing F_{ice} and melting M_{ice} of sea ice, and \varphi^{(m)}_{_W} is the tracer concentration in the freshwater associated with source m. The change in volume of the surface layer due to the freshwater flux is discussed in Sec.:ref:sec-linear-free-surface, and the natural boundary conditions for tracers associated with freshwater flux are discussed in Sec. Variable-Thickness Surface Layer. The boundary conditions on tracers are no-flux normal to bottom and lateral boundaries, unless modified by parameterizations like the overflow parameterization (see Sec. Overflow Parameterization).

4.3.1. Tracer Advection

POP has multiple options for tracer advection. Here we describe a basic second-order centered advection scheme as an example of a discretization. Other advection schemes are described in Chapter Advection. In the standard second-order scheme, the advection operator is given by:

(37){\cal L}_{T}(\varphi) = \frac{1}{\Delta_y} \delta_x
(\overline{\Delta_y u_x}^y \overline{\varphi}^x)
+ \frac{1}{\Delta_x} \delta_y
(\overline{\Delta_x u_y}^x \overline{\varphi}^y)
+ \delta_z(w \overline{\varphi}^z)

Again, \Delta_x and \Delta_y inside the difference operator are located at U-points, and the mass fluxes \overline{\Delta_y u_x}^y, \overline{\Delta_x u_y}^x are located on the lateral faces of T-cells. \overline{\varphi}^x and \overline{\varphi}^y are also located on the lateral faces of T-cells, while \overline{\varphi}^z is located on the top and bottom faces of T-cells. At the surface, \overline{\varphi}^z is set equal to zero because there is no advection of tracers across the surface. The vertical velocity w at T-points is determined from the solution of the continuity equation

(38)\frac{1}{\Delta_y} \delta_x
(\overline{\Delta_y u_x}^y)
+ \frac{1}{\Delta_x} \delta_y
(\overline{\Delta_x u_y}^x)
+ \delta_z (w)=0

which is integrated in a column of T-cells downward from the top with the boundary conditions:

(39)w=\frac{\partial}{\partial t}\eta + q_{_{W}}\;\;\;\mathrm{at}\;\;\;z=0\;

(40)w=0\;\;\;\mathrm{at}\;\;\;z=-H_T\;.

The boundary condition (39) is discussed in more detail in Sec.:ref:sec-linear-free-surface.

When integrating from the top down using (38) and (39), the bottom boundary condition (40) will only be satisfied to the extent that the solution of the elliptic system of barotropic equations for the surface height \eta has converged exactly (see Sec. Barotropic Equations). In practice, this exact convergence is not achieved, so the bottom velocity in T-columns is set to zero to ensure no tracers are fluxed through the bottom. This amounts to allowing a very small divergence of velocity in the bottom cell.

In versions of POP prior to 1.4.3, the volume of the surface cells is assumed constant, and no account is taken of the change in volume of the surface cells when the free surface height changes. Thus \xi = 0 and q_{_{W}} = 0 in (34), (35) and (36), and the freshwater flux is approximated as a virtual salinity flux and imposed as a boundary condition on the vertical diffusion operator. In addition, the tracers are advected through the surface in this formulation using (37) with the vertical velocity given by (39) and \overline{\varphi}^z = \varphi_1 at the surface. One problem with this approximation is that the advective flux of tracers through the surface is not zero in global average. The globally integrated vertical mass flux vanishes, but the integrated tracer flux does not. In practice, we have found that the residual surface tracer fluxes associated with this are usually small, but in some situations they may be non-negligible. (Note: the global mean residual surface tracer fluxes are standard diagnostic model output in the earlier versions of POP without the variable thickness surface layer.) Because the residual surface flux is nonzero, the global mean tracers are not conserved. In POP version 1.4.3, the surface layer is based on (34), (35) and (36), and global mean tracers (in particular total salt) are conserved exactly. The variable-thickness surface layer is detailed in Sec. Variable-Thickness Surface Layer.

4.3.2. Horizontal Tracer Diffusion

POP has three options for horizontal tracer diffusion: 1) horizontal Laplacian diffusion, 2) horizontal biharmonic diffusion, 3) the Gent-McWilliams parameterization, which includes along-isopycnal tracer diffusion and tracer advection with an additional eddy-induced transport velocity. All of these are implemented for a spatially varying diffusivity. The first of these options (Laplacian diffusion) is described here, the other two options are discussed in Sec. Horizontal Tracer Diffusion. The discrete horizontal Laplacian diffusion operator is given by:

(41){\cal D}_H(\varphi) =
  \frac{1}{\Delta_y} \delta_x
(\overline{A_H}^x\Delta_y \delta_x \varphi )
+ \frac{1}{\Delta_x} \delta_y
(\overline{A_H}^y\Delta_x \delta_y \varphi)

Note that no lateral averaging is involved as in (33). Thus the Laplacian is approximated as a five-point stencil as would be used on a C-grid. The factors \Delta_x, \Delta_y inside the difference operator are given by HTN, HTE, respectively (see Figure3.1 ), and A_H, defined at T-points, is averaged across the T-cell faces. As mentioned above, all the horizontal diffusive operators in POP that would normally involve nine-point operators (including the GM parameterization and the horizontal friction operators), are approximated by five-point C-grid operators in order to ensure that they damp checkerboard noise on the grid scale. B-grid Laplacian-type operators like (33) have a checkerboard null space, i.e., they yield zero when applied to a +/- checkerboard field, and thus cannot damp noise of this character (see Sec. Nullspace Removal). The only Laplacian-type operator which uses a B-grid discretization is the elliptic operator in the implicit barotropic system. There the B-grid discretization is required in order to maintain energetic consistency (Smith, Dukowicz, and Malone 1992; Dukowicz, Smith, and Malone 1993), see also Secs. Energetic Consistency and Barotropic Equations). The boundary conditions on the diffusive operator (41) are that tracer gradients \delta_x \varphi, \delta_y \varphi are zero normal to lateral boundaries.

4.3.3. Vertical Tracer Diffusion

The spatial discretization of the vertical diffusion operator is given by:

(42)\begin{aligned}
{\cal D}_V(\varphi) & = & \delta_z(\kappa\delta_z\varphi)
\nonumber \\
  & = &
\frac{1}{dz_k} \left(
 \frac{\kappa_{k-\frac{1}{2}}(\varphi_{k-1}-\varphi_{k})}
{dz_{k-\frac{1}{2}}} -\frac{\kappa_{k+\frac{1}{2}}
(\varphi_{k}-\varphi_{k+1})}{dz_{k+\frac{1}{2}}} \right).
\end{aligned}

where \varphi_{k} is a tracer at level k, and \kappa_{k-\frac{1}{2}}, \kappa_{k+\frac{1}{2}} are evaluated on the top and bottom faces, respectively, of the T-cell at level k, and dz_{k-\frac{1}{2}}=dzw_{k-1}, dz_{k+\frac{1}{2}}=dzw_{k}. The boundary conditions at the top and bottom of the column are

(43)\begin{aligned}
\kappa\delta_z\varphi \rightarrow Q_{\varphi} &
\mathrm{at} & z=0
\nonumber \\
\kappa\delta_z\varphi \rightarrow 0 &
\mathrm{at} & z=-H_T
\end{aligned}

where Q_{\varphi} is the surface flux of tracer \varphi (e.g., heat flux for temperature and equivalent salt flux associated with freshwater flux for salinity). The modifications to this discretization when partial bottom cells are used is described in Sec. Partial Bottom Cells. The diffusive term may either be evaluated explicitly or implicitly. The implicit treatment is described in Sec. Implicit Vertical Diffusion. With explicit mixing, a convective adjustment routine may also be used to more efficiently mix tracers when the column is unstable (see Sec. Convection). Various subgrid-scale parameterizations for the vertical diffusivity are discussed in Sec. Vertical Mixing.

4.4. Discrete Momentum Equations.

The momentum equations discretized on the B-grid are given by:

(44)\frac{\partial}{\partial t}u_x + {\cal L}_{U} (u_x)
+ u_xu_yk_y - u^2_y k_x -fu_y =
-\frac{1}{\rho_{_0}} \delta_x\overline{p}^y
+ {\cal F}_{Hx} (u_x,u_y) + {\cal F}_V (u_x)

(45)\frac{\partial}{\partial t}u_y + {\cal L}_{U} (u_y)
+ u_xu_yk_x - u^2_x k_y +fu_x =
-\frac{1}{\rho_{_0}}\; \delta_y\overline{p}^x
+ {\cal F}_{Hy} (u_y,u_x) + {\cal F}_V (u_y)

In these equations no account has been taken of the change in volume of the surface layer due to undulations of the free surface. Therefore, no terms involving \xi appear as in the tracer transport equation (34). The justification for this is that the global mean momentum, unlike the global mean tracers, is not conserved in the absence of forcing, so there is less motivation to correct for momentum nonconservation due to surface height fluctuations. Furthermore, the error introduced is typically small compared to the uncertainty in the applied wind stress.

Note: Currently the code is in cgs units and it is assumed that \rho_{_0}=1.0g cm^{-3}, so it never explicitly appears. If the Boussinesq correction (Sec. Boussinesq Correction) is used, then this factor is already taken into account, because the factor r(p) in (229) is normalized such that the pressure gradient should be divided by \rho_{_0} = 1.0.

4.4.1. Momentum Advection

The nonlinear momentum advection term is discretized as:

(46){\cal L}_{U}(\alpha) = \frac{1}{\Delta_y} \delta_x
\big[\overline{(\overline{\Delta_y u_x}^y)}^{xy}
\overline{\alpha}^x\big]
+ \frac{1}{\Delta_x} \delta_y
\big[\overline{(\overline{\Delta_x u_y}^x)}^{xy}
\overline{\alpha}^y\big]
+ \delta_z (w^U\overline{\alpha}^z)\;.

This is a second-order centered advection scheme and is currently the only option available in POP for momentum advection. It has the property that global mean kinetic energy is conserved by advection. Momentum is conserved in the interior by advection, but not on the boundaries (see Sec. Energetic Consistency on energetic consistency). The mass fluxes in the operator include an extra average in both horizontal directions, denoted \overline{(..)}^{xy}. Both the horizontal and vertical mass fluxes in a U-cell are the average of the four surrounding T-cell mass fluxes. This is illustrated in Figure3.4 for the mass flux in the x-direction. As a consequence, the vertical velocity in a U-cell is exactly the area-weighted average of the four surrounding T-cell vertical velocities w. This averaging is necessary in order to maintain the energetic balance between the global mean work done by the pressure gradient and the change in gravitational potential energy (Smith, Dukowicz, and Malone 1992), see also Sec. Energetic Consistency on energetic consistency. An additional advantage of this flux averaging in U-cells is that it substantially reduces noise in the vertical velocity field compared to other approaches (Webb 1995). At the bottom of a column of ocean U-cells w^U is not necessarily zero, because it is the weighted average of the surrounding w’s, some of which may be nonzero if it is a “rim” point where k=\mathrm{KMU} but at least one of the surrounding T-points has k<\mathrm{KMT}. It can be shown that the value of w^U at these points approximates the boundary condition for tangential flow along the sloping bottom w = -{\bf u}\cdot\nabla H (Semtner 1986). Momentum is advected through the surface using the vertical velocity (39) averaged to U-points and with \overline{\alpha}^z = \alpha_1 (where \alpha = u_x or u_y) in (46).

../_images/FigTadvect.jpg

Figure 3.4: Advective mass fluxes through lateral faces of T-cells (narrow solid lines) and U-cells (thick solid line). U-cell fluxes are the average of the four surrounding T-cell fluxes. A similar averaging applies to the vertical fluxes.

4.4.2. Metric Terms

The metric terms (third and fourth terms on the l.h.s. in (45)) are constructed using a simple C-grid discretization of the metric factors, as in (26) and (27). Specifically, k_x and k_y at U-points are:

(47)\begin{aligned}
\mathrm{KXU}_{i,j}&=&[\mathrm{HUW}_{i+1,j}-\mathrm{HUW}_{i,j}]/\mathrm{UAREA}_{i,j}
\nonumber \\
\mathrm{KYU}_{i,j}&=&[\mathrm{HUS}_{i,j+1}-\mathrm{HUS}_{i,j}]/\mathrm{UAREA}_{i,j}
\end{aligned}

4.4.3. Horizontal Friction

Several options for horizontal friction are available in the code. In this section a simple spatial discretization of the Laplacian-type formulation with a spatially varying viscous coefficient as given by (29) is presented. The biharmonic friction operator obtained by applying (28) twice is described in Sec. Biharmonic Horizontal Viscosity. More sophisticated formulations of the viscosity based on a functional discretization of the friction operator formulated as the divergence of a viscous stress that is linearly related to the components of the strain-rate tensor are described in Sections Anisotropic Horizontal Viscosity and Smagorinsky Nonlinear Viscous Coefficients. These include an anisotropic formulation of the viscosity and the use of Smagorinksy-type nonlinear viscous coefficients.

The Laplacian horizontal friction terms (28) and (29) are constructed from a C-grid discretization of both the Laplacian and the metric terms. The discrete Laplacian terms are given by:

(48)\begin{aligned}
&&\nabla\cdot A_M \nabla{\bf u} =
  \frac{1}{\Delta_y} \delta_x
(\overline{A_M}^x\Delta_y \delta_x {\bf u} )
+ \frac{1}{\Delta_x} \delta_y
(\overline{A_M}^y\Delta_x \delta_y {\bf u})
 \end{aligned}

where A_M, defined at U-points, is averaged across cell faces inside the divergence. Terms proportional to k_x and k_y in (28) and (29) are evaluated using (47). Terms involving derivatives of k_x and k_y are evaluated with k_x, k_y defined at T-points and averaged along U-cell faces. For example, the term \delta_x A_M k_x is evaluated as:

\begin{aligned}
&&\delta_x \overline{A_M}^x\overline{k_x}^y =
 \nonumber \\
&&0.25\{[(A_M)_{i+1,j}+(A_M)_{i,j}][\mathrm{KXT}_{i+1,j+1}+\mathrm{KXT}_{i+1,j}]
 \nonumber \\
&&-[(A_M)_{i,j}+(A_M)_{i-1,j}][\mathrm{KXT}_{i,j+1}+\mathrm{KXT}_{i,j}]\}
/\mathrm{DXU}_{i,j}\end{aligned}

where A_M is averaged across the U-cell faces, and k_x, k_y at T-points are given by

(49)\begin{aligned}
\mathrm{KXT}_{i,j}&=&[\mathrm{HTE}_{i,j}-\mathrm{HTE}_{i-1,j}]/\mathrm{TAREA}_{i,j}
\nonumber \\
\mathrm{KYT}_{i,j}&=&[\mathrm{HTN}_{i,j}-\mathrm{HTN}_{i,j-1}]/\mathrm{TAREA}_{i,j}
\end{aligned}

Finally, in those terms in (28),:eq:horiz-fric.2 that involve single derivatives of the velocities or viscosities (e.g. \delta_x u_y or, \delta_x A_M) the derivatives are evaluated as differences across the cell without using the central value. For example, \delta_x u_y at point (i,j) is evaluated as:

\begin{aligned}
[(u_y)_{i+1,j}-(u_y)_{i-1,j}]/[\mathrm{HTN}_{i,j}+\mathrm{HTN}_{i+1,j}]\end{aligned}

The no-slip boundary conditions are implemented by simply setting u_x = u_y = 0 on all lateral boundary points. Modifications to the operator when partial bottom cells are used are described in Sec. Partial Bottom Cells.

4.4.4. Vertical Friction

The spatial discretization of the vertical friction terms is essentially identical to that of the vertical diffusion described in Sec. Vertical Tracer Diffusion. The spatial discretization of both {\cal F}_V (u_x) and {\cal F}_V (u_y) is identical to (42) with \kappa replaced by the vertical viscosity \mu and \varphi replaced by one of the velocity components u_x or u_y. Modifications for partial bottom cells are discussed in Sec. Partial Bottom Cells. The boundary conditions at the top and bottom of the column of U-points are

(50)\begin{aligned}
\mu\delta_z (u_x,u_y) \rightarrow (\tau_x,\tau_y) &
\mathrm{at} & z=0
\nonumber \\
\mu\delta_z (u_x,u_y) \rightarrow c|{\bf u}|(u_x,u_y) &
\mathrm{at} & z=-H_U
\end{aligned}

where (\tau_x,\tau_y) are the components of the surface wind stress along the coordinate directions, and a quadratic drag term is applied at the bottom of the column. The dimensionless constant c is typically chosen to be of order 10^{-3}. The semi-implicit treatment of these terms is described in Sec. Semi-Implicit Treatment of Coriolis and Vertical Friction Terms. Various subgrid-scale parameterizations for the vertical viscosity are described in Sec. Vertical Mixing.

4.5. Energetic Consistency

The energetic balances in the free-surface formulation are described in detail in Dukowicz and Smith (1994). However, this reference predates the implementation of a variable-surface thickness layer (Sec. Variable-Thickness Surface Layer) and partial bottom cells (Sec. Partial Bottom Cells), both of which influence the energetic balances.

5. Time Discretization

The POP model uses a three-time-level second-order-accurate modified leapfrog scheme for stepping forward in time. It is modified in the sense that some terms are evaluated semi-implicitly, and of the terms that are treated explicitly, only the advection operators are actually evaluated at the central time level, as in a pure leapfrog scheme. The diffusive terms are evaluated using a forward step. The reason for this is that the centered advection scheme is unstable for forward steps, and the diffusive scheme is unstable for leapfrog steps.

5.1. Filtering Timesteps

Leapfrog schemes can develop computational noise due to the partial decoupling of even and odd timesteps. In a pure leapfrog scheme they are completely decoupled and the solutions on the even and odd steps can evolve independently, leading to 2\Delta t oscillations in time. There are several methods to damp the leapfrog computational mode, two of which are currently implemented in POP.

One is to occasionally take a forward step or an Euler forward-backward step, sometimes called a Matsuno timestep (Haltiner and Williams 1980). The Matsuno step is a forward predictor step followed by a “backward” step which is essentially a repeat of the forward step but using the predicted prognostic variables from the first pass to evaluate all terms except the time-tendency term (see Sec. Tracer Transport Equations). The Matsuno step is more expensive than a forward step, but it is stable for advection.

The other method is to occasionally perform an averaging of the solution at three successive time levels to the two intermediate times, back up half a timestep and proceed. The later procedure is referred to as an “averaging timestep” (Dukowicz and Smith 1994) and is the recommended method for eliminating the leapfrog computational mode. The leapfrog scheme generates two different “trajectories” of the solution, one corresponding to even and one to odd steps. The advantage of the averaging step is that it places the solution on the average trajectory, whereas the forward and Matsuno steps select only one trajectory, corresponding to either the even or the odd solution. Experience has shown that some model configurations are not stable using Matsuno filtering timesteps, and this is especially true with the variable-thickness surface layer (Sec. Variable-Thickness Surface Layer).

On the very first time step of a spinup from rest a forward step is taken to avoid immediately exciting a leapfrog computational mode (this feature is hardwired into the code). In the presentation below, the time discretization of various terms will first be presented for a regular leapfrog step, then the discretization for the forward and backward steps will be given.

5.2. Tracer Transport Equations

Labeling the three time levels on a given step as n-1, n, and n+1, the tracer transport equation (34) is discretized in time as follows:

(51)(1+\xi^{n+1})\varphi^{n+1}-(1+\xi^{n-1})\varphi^{n-1} = \\
\tau [- {\cal L}^{(n)}_T(\varphi^n) + {\cal D}_H(\varphi^{n-1})
+{\cal D}_V(\overline{\varphi}^{\lambda})
+ {\cal F}^n_{_W}]

\overline{\varphi}^{\lambda} =
\lambda\varphi^{(n+1)} + (1-\lambda)\varphi^{(n-1)}

where \Delta t is the timestep and \tau = 2\Delta t. {\cal F}^n_{_W} is given by (36) with q_{_{W}} and \varphi evaluated at time n, and \xi^n is given by (35) with \eta evaluated at time n. The superscript (n) on the advection operator indicates that the advective mass fluxes are evaluated using the time n velocities. The vertical diffusion term may be evaluated either explicitly or semi-implicitly. In the explicit case \lambda=0. For semi-implicit diffusion \lambda\ge 1/2 is required for stability, and the code is usually run with \lambda=1. The surface forcing (43), applied as a boundary condition on the vertical diffusion operator, is evaluated at time n for both explicit and semi-implicit mixing. The modifications to the tracer transport equations with implicit vertical mixing are described in Sec. Implicit Vertical Diffusion.

With pressure averaging (see Sec. Pressure Averaging) the potential temperature and salinity at the new time are needed to evaluate the pressure gradient at the new time. This is required in the baroclinic momentum equations before the barotropic equations have been solved for the surface height at the new time. Therefore, (51) cannot be used to predict the new tracers because \xi^{n+1} is not yet known. Instead an approximation is made to predict the new temperature and salinity which are then used to evaluate the pressure gradient at the new time. After the barotropic equations have been solved, the new potential temperature and salinity are corrected so that they satisfy (51) exactly. The equations for predicting and correcting the tracers at the new time are obtained as follows. The l.h.s. of (51), which involves the unknown surface height \eta^{n+1}, is approximated as:

\begin{aligned}
 & (1+\xi^{n+1})\varphi^{n+1}-(1+\xi^{n-1})\varphi^{n-1}
    = (\varphi^{n+1}-\varphi^{n-1}) + & \nonumber \\
& \frac{1}{2}\big(\xi^{n+1}+\xi^{n-1}\big)
(\varphi^{n+1}-\varphi^{n-1})
+ (\xi^{n+1}-\xi^{n-1})\frac{1}{2}
\big(\varphi^{n+1}+\varphi^{n-1}\big) &
\nonumber \\
& \approx (1+\xi^n)(\varphi^{n+1}-\varphi^{n-1})
+ 2(\xi^n - \xi^{n-1})\varphi^n & \end{aligned}

With this approximation, the equations for predicting and correcting \Theta and S are given by:

Predictor:

(52)(1+\xi^n)(\widehat{\varphi}^{n+1}-\varphi^{n-1})
= - 2(\xi^n-\xi^{n-1})\varphi^n + \tau F

Corrector:

(53)(1+\xi^{n+1})\varphi^{n+1} =
(1+\xi^n)\widehat{\varphi}^{n+1}
+(\xi^n-\xi^{n-1})(2\varphi^n-\varphi^{n-1})

where \widehat{\varphi}^{n+1} is the predicted tracer at the new time, and F represents all terms in brackets on the r.h.s. of (51). Note that these equations are used only to predict and correct \Theta and S; all passive tracers are updated directly using (51).

The time discretization for the two-time-level forward and backward steps is given by

Forward step:

(54)\begin{aligned}
(1+\xi^{*})\varphi^{*}-(1+\xi^{n})\varphi^{n} & = &
\tau [- {\cal L}^{(n)}_T(\varphi^n) + {\cal D}_H(\varphi^{n})
+{\cal D}_V(\overline{\varphi}^{\lambda})
+ {\cal F}^n_{_W}] \nonumber \\
\overline{\varphi}^{\lambda} & = &
\lambda\varphi^{*} + (1-\lambda)\varphi^{n}
\end{aligned}

Backward step:

(55)\begin{aligned}
(1+\xi^{n+1})\varphi^{n+1}-(1+\xi^{n})\varphi^{n} & = &
\tau [- {\cal L}^{(*)}_T(\varphi^*) + {\cal D}_H(\varphi^{*})
+{\cal D}_V(\overline{\varphi}^{\lambda})
+ {\cal F}^n_{_W}] \nonumber \\
\overline{\varphi}^{\lambda} & = &
\lambda\varphi^{n+1} + (1-\lambda)\varphi^{n}
\end{aligned}

Here \tau = \Delta t and \varphi^*, the predicted tracer at the new time from the forward step, is used to evaluate the r.h.s. in the backward step. \xi^{*} is given by (35) with \eta=\eta^{*}. {\cal L}^{(*)}_T is the advection operator evaluated using the predicted velocities from the forward step (see Sec. Baroclinic Momentum Equations). For a forward step only, \varphi^*=\varphi^{n+1}. The surface forcing applied to {\cal D}_V, as well as the freshwater tracer flux {\cal F}_{_W}, are evaluated at time n in both forward and backward steps.

5.2.1. Tracer Acceleration

The acceleration technique of Bryan (1984) can be used to more quickly spin up the model to a steady state by modifying the timestep for the momentum equation, or tracers in the deep ocean. The basic method is to modify the time step \tau as follows in the baroclinic and barotropic momentum equations: (65) and (83)

\begin{aligned}
\tau \rightarrow \tau^{(u)}_k = \tau^{(u)}/\alpha\end{aligned}

and in the tracer transport equation (51)

\begin{aligned}
\tau \rightarrow \tau^{(\varphi)}_k = \tau^{(\varphi)}/\gamma_k\end{aligned}

where superscripts denote timesteps for the momentum and continuity (u), and tracers (\varphi). \alpha and \gamma_k are acceleration factors specified in the namelist model input: dtuxcel = \alpha^{-1}, dttxcel(k) = (\gamma_k)^{-1}, in addition to the tracer timestep. Note that \tau^{(u)} is not a model input.

There are some points to make about the use of acceleration. First, a disadvantage of depth-dependent tracer acceleration is that it leads to non-conservative advective and diffusive fluxes when \gamma_k \neq \gamma_{k+1}. For this reason, it is recommended that profiles of \gamma_k be smooth functions of depth, and that the largest vertical discontinuities in \gamma_k be restricted to depths where fluxes are small enough that tracers are well conserved. The depth-dependent tracer timestep must also be accounted for in the convective adjustment and vertical mixing routines. These changes are described in Secs. Convection, Tridiagonal Solver, and Semi-Implicit Treatment of Coriolis and Vertical Friction Terms. Danabasoglu (2004) shows some detrimental effects of the tracer non-conservation with this type of acceleration method, including the model approaching an incorrect steady state solution. Therefore, it is recommended that the depth-dependent tracer acceleration not be used.

Second, if \alpha\neq \gamma_1, then the surface forcing for momentum does not occur at the correct time because the model calender and forcing fields are based on the surface tracer time step. For this reason, if momentum acceleration, \alpha\neq 1, is used, then a subsequent synchronous integration is required to get the momentum and tracer fields back in step. This acceleration technique is recommended by Danabasoglu (2004), who used a value of \alpha = 5. Third, if the variable thickness surface layer option is used with virtual salt flux boundary conditions, see Sec. Variable-Thickness Surface Layer, Danabasoglu (2004) shows the momentum acceleration technique is stable, despite the fact that the time stepping discretizations in the tracer transport and barotropic continuity equations are inconsistent.

5.2.2. Implicit Vertical Diffusion

The tracer equations with implicit vertical diffusion involve the solution of a tridiagonal system in each vertical column of grid points. This is a relatively easy thing to do because the system does not involve any coupling with neighboring points in the horizontal direction. The equations solved in the code are

(56)\begin{aligned}
[1+\xi^{n+1}-\lambda \tau {\sf A}(\kappa)]\Delta\varphi & = &
-(\xi^{n+1}-\xi^{n-1})\varphi^{n-1} \nonumber \\
 & + & \tau [- {\cal L}_T^{(n)}(\varphi^n) +
{\cal D}_H(\varphi^{n-1}) + {\sf A}(\kappa)\varphi^{n-1}
+ {\cal F}^n_{_W}]
\nonumber \\
{\sf A}(\kappa) = \delta_z \kappa \delta_z,
&& \Delta\varphi\equiv\varphi^{n+1}-\varphi^{n-1}
\end{aligned}

where \tau=2\Delta t. For explicit diffusion \lambda=0, so the term in brackets on the right corresponds to the exact r.h.s. in the explicit case. The system is solved for the change in tracer \Delta\varphi, subject to the no-flux boundary conditions

(57)\begin{aligned}
\delta_z\Delta\varphi=0 \mathrm{\ at \  both \ } z=0
\mathrm{\ and \ } z=-H_T\;.
\end{aligned}

Note that because the surface forcing and freshwater tracer fluxes are evaluated at time n, they are entirely contained in the r.h.s. of (56) and hence are not directly imposed as a boundary condition on the operator.

The predictor and corrector steps for updating \Theta and S with pressure averaging (Sec. Pressure Averaging) are given by:

Predictor:

(58)\begin{aligned}
[1+\xi^{n}-\lambda \tau {\sf A}(\kappa)]\Delta\varphi & = &
-2(\xi^{n}-\xi^{n-1})\varphi^n  + \tau F \nonumber \\
\Delta\varphi & \equiv & \widehat{\varphi}^{n+1}-\varphi^{n-1}
\end{aligned}

Corrector:

(59)\begin{aligned}
[1+\xi^{n+1}-\lambda \tau {\sf A}(\kappa)]\Delta\varphi & = &
(\xi^n-\xi^{n-1})(2\varphi^n-\varphi^{n-1})
-(\xi^{n+1}-\xi^n)\widehat{\varphi}^{n+1} \nonumber \\
\Delta\varphi & \equiv & \varphi^{n+1}-\widehat{\varphi}^{n+1}
 \end{aligned}

where F represents all terms in brackets on the right-hand side of (56). These reduce to (52) and (53) in the limit \lambda = 0. Note that both the predictor and corrector steps involve the solution of a tridiagonal system.

In forward and backward timesteps the implicit equations have the form:

Forward step:

(60)\begin{aligned}
[1+\xi^{*}-\lambda \tau {\sf A}(\kappa)]\Delta\varphi & = &
-(\xi^*-\xi^n)\varphi^n \nonumber \\
 & + & \tau [- {\cal L}_{T}^{(n)}(\varphi^n)
+ {\cal D}_H(\varphi^n) + {\sf A}(\kappa)\varphi^{n}
+ {\cal F}^n_{_W}] \nonumber \\
\Delta\varphi & \equiv & \varphi^{*}-\varphi^{n}
\end{aligned}

Backward step:

(61)\begin{aligned}
[1+\xi^{n+1}-\lambda \tau {\sf A}(\kappa)]\Delta\varphi & = &
-(\xi^{n+1}-\xi^n)\varphi^n \nonumber \\
& + & \tau [- {\cal L}_{T}^{(*)}(\varphi^*)
+ {\cal D}_H(\varphi^*) + {\sf A}(\kappa)\varphi^{n}
+ {\cal F}^n_{_W}] \nonumber \\
\Delta\varphi & \equiv & \varphi^{n+1}-\varphi^{n}
\end{aligned}

where \tau = \Delta t.

5.2.3. Tridiagonal Solver

For the tridiagonal solution of (56), (58), (59), (60) or (61), a new algorithm is used (Schopf and Loughe 1995) which is more accurate and stable than traditional methods. These equations have the form:

(62)(1+\xi)X_k - \lambda\frac{\tau_k}{dz_k}\big[
\frac{\kappa_{k-\frac{1}{2}}}{dz_{k-\frac{1}{2}}}
(X_{k-1} - X_{k})
- \frac{\kappa_{k+\frac{1}{2}}}{dz_{k+\frac{1}{2}}}
(X_{k} - X_{k+1})\big] = R_k

where X_k = \Delta\varphi at level k, R_k is the r.h.s., and a depth-dependent timestep \Delta t_k (\tau_k = 2\Delta t_k for leapfrog and \tau_k=\Delta t_k for forward/backward timesteps) is used with tracer acceleration (see Sec. Tracer Acceleration). Here \xi=\xi^{n+1}, except in the predictor step (58) where \xi=\xi^{n}. (62) can be rewritten in the form of a tridiagonal system:

(63)\begin{aligned}
-A_k X_{k+1} + C_k X_k & - & A_{k-1} X_{k-1} = R'_k \\
A_k & = & \lambda\frac{\kappa_{k+\frac{1}{2}}}{dz_{k+\frac{1}{2}}} \nonumber \\
C_k & = & h'_k + A_k + A_{k-1} \nonumber \\
R'_k & = & h_k R_k \nonumber \\
h_k & = & \frac{dz_k}{\tau_k} \nonumber \\
h'_k & = & h_k(1+\xi) \nonumber
\end{aligned}

The algorithm for the solution of this system involves a loop over vertical levels to determine the coefficients:

(64)\begin{aligned}
E_k & = & A_k/D_k \\
F_k & = & (R'_k + A_{k-1} F_{k-1}) / D_k \nonumber \\
& \mathrm{where} & \nonumber \\
B_k & = & A_k ( h'_k + B_{k-1})/ D_k \nonumber \\
D_k & = & h'_k + A_k + B_{k-1} \nonumber
\end{aligned}

The loop begins at k=1 with A_0 = B_0 = F_0 =0. This is followed by another vertical loop to determine the solution by back substitution:

\begin{aligned}
X_k = E_k X_{k+1} + F_k \end{aligned}

This loop begins at the bottom with X_{k+1}=0 when k=\mathrm{KMT}.

5.3. Splitting of Barotropic and Baroclinic Modes

The barotropic mode of the primitive equations supports fast gravity waves with speeds of \sqrt{gH}\sim 200 m s:math:^{-1}. If resolved numerically, these waves impose a severe restriction on the model timestep. However, they have little effect on the dynamics, especially on timescales longer than a day or so. To overcome this severe limitation on the timestep, the barotropic mode is split off and solved as a separate 2-D system (see Sec. Barotropic Equations).

The barotropic equations are taken to be the vertically integrated momentum and continuity equations. The true barotropic mode is not exactly isolated by the vertically integrated system, except in the limit of a flat-bottom topography, a rigid lid, and a depth-independent buoyancy frequency. Nevertheless, it suffices in practice to isolate and treat the vertically integrated system. Then, provided advective and diffusive CFL limits do not control the timestep, the baroclinic equations can be integrated with a timestep that is controlled by the gravity-wave speed of the first internal mode, which is of order 2 m s^{-1}, two orders of magnitude smaller than the barotropic wave speed. The procedure for solving the split barotropic-baroclinic system is as follows.

1) First the momentum equations are solved, without including the surface pressure gradient, for an auxiliary velocity {\bf u}'. This is what the momentum at the new time would be in the absence of the surface pressure gradient, which is depth-independent and hence determined only by the solution of the vertically integrated system. The time discretization of the resulting “baroclinic momentum equations” written in vector form is:

(65)\begin{aligned}
\frac{{\bf u}'-{\bf u}^{n-1}}{\tau} + \widetilde{\cal L}_U({\bf u}^n)
-f\hat{\bf z}\times\overline{{\bf u}}^{\alpha\gamma} =
-\frac{1}{\rho_o}\nabla\overline{p_h} + {\bf{\cal F}}_{H}({\bf u}^{n-1})
+ {\bf{\cal F}}_{V}(\overline{{\bf u}}^{\lambda})
\end{aligned}

where \tau = 2\Delta t, \widetilde{\cal L}_U represents the advection operator plus metric terms acting on both components of velocity, and {\bf{\cal F}}_{H} and {\bf{\cal F}}_{V} are now horizontal vectors. The overbars indicate various averages over the three time levels for the semi-implicit treatment of the Coriolis, hydrostatic pressure gradient, and vertical mixing terms. The velocities are averaged as:

(66)\overline{{\bf u}}^{\alpha\gamma} =
\alpha{\bf u}' + \gamma{\bf u}^n + (1-\alpha-\gamma){\bf u}^{n-1}

(67)\overline{\bf u}^{\lambda} =
\lambda{\bf u}'+(1-\lambda){\bf u}^{n-1}

where \alpha,\gamma,\lambda are coefficients used to vary the time-centering of the velocities. The averaging for the hydrostatic pressure \overline{p_h} is discussed in Sec. Hydrostatic Pressure. To maintain an accurate time discretization of geostrophic balance, it is important that, in the averaging over the three time levels, the velocity and pressure are centered at the same time, i.e. if they are centered at time (n), then the coefficients for the variables at times (n+1) and (n-1) must be equal.

2) Next the vertical average of {\bf u}' is subtracted. The result is the baroclinic velocity:

(68)\begin{aligned}
\tilde{\bf u}'_k = {\bf u}'_k -
\frac{1}{H_U}\sum_{k'=1}^{km} dz_{k'}\;{\bf u}'_{k'}
\end{aligned}

3) Finally, the barotropic system is solved for the barotropic velocity at the new time {\bf U}^{n+1} (see Sec. Barotropic Equations), and this is added to the normalized auxiliary velocities to obtain the full velocities at the new time:

(69)\begin{aligned}
{\bf u}^{n+1}_k = \tilde{\bf u}'_k + {\bf U}^{n+1}
\end{aligned}

The barotropic velocity is defined by

(70)\begin{aligned}
{\bf U} &=& \frac{1}{H+\eta}\int^{\eta}_{\!-H}\!dz\;{\bf u}(z)
\nonumber\\
&\approx &\frac{1}{H}\int^{0}_{\!-H}\!dz\;{\bf u}(z)
= \frac{1}{H_U}\sum_{k=1}^{km} dz_k\;{\bf u}_{k}\;.
\end{aligned}

where \eta is the displacement of the free surface relative to z=0. As discussed in Secs. Barotropic Equations and Variable-Thickness Surface Layer, in the current version of the POP model we assume \eta\ll dz_1 and approximate the barotropic velocity as the vertical integral from z=-H to z=0 in (70).

5.4. Baroclinic Momentum Equations

The time discretization of the baroclinic momentum equations for leapfrog steps is given by (65). For forward and backward steps, the time discretization is

Forward steps:

(71)\begin{aligned}
\frac{{\bf u}'^*-{\bf u}^{n}}{\tau} + \widetilde{\cal L}_U({\bf u}^n)
-f\hat{\bf z}\times\overline{{\bf u}}^{\theta} & = &
-\frac{1}{\rho_o}\nabla p^n_h + {\bf{\cal F}}_{H}({\bf u}^{n})
+ {\bf{\cal F}}_{V}(\overline{{\bf u}}^{\lambda'})
\nonumber \\
\overline{{\bf u}}^{\theta} & = &
\theta{\bf u}'^* + (1-\theta){\bf u}^{n}
\nonumber \\
\overline{{\bf u}}^{\lambda'} & = &
\lambda'{\bf u}'^* + (1-\lambda'){\bf u}^{n}
\end{aligned}

Backward steps:

(72)\begin{aligned}
\frac{{\bf u}'-{\bf u}^{n}}{\tau} + \widetilde{\cal L}_U({\bf u}'^*)
-f\hat{\bf z}\times\overline{{\bf u}}^{\theta} & = &
-\frac{1}{\rho_o}\nabla p^*_h + {\bf{\cal F}}_{H}({\bf u}'^*)
+ {\bf{\cal F}}_{V}(\overline{{\bf u}}^{\lambda'})
\nonumber \\
\overline{{\bf u}}^{\theta} & = &
\theta{\bf u}' + (1-\theta){\bf u}^{n}
\nonumber \\
\overline{{\bf u}}^{\lambda'} & = &
\lambda'{\bf u}' + (1-\lambda'){\bf u}^{n}
\end{aligned}

where \tau = \Delta t, and p^*_h is the predicted hydrostatic pressure from the forward step. {\bf u}' is the unnormalized momentum, as in (65) and (68), and {\bf u}'^* is the same quantity predicted by the forward step only. \theta = 0.5 is hardwired into the code. This choice was made (rather than choosing \theta = 0 where both Coriolis and pressure gradient would be centered and time n) because the forward step is unstable if the Coriolis terms is treated explicitly.

5.4.1. Pressure Averaging

The method of Brown and Campana (1978) is used for the semi-implicit treatment of the hydrostatic pressure gradient on leapfrog timesteps, where the averaging of the pressure over the three time levels in (65) is given by

(73)\begin{aligned}
\frac{1}{\rho_o}\nabla\overline{p_h} =
\frac{1}{\rho_o}\nabla(\frac{1}{4}p^{n+1}
+\frac{1}{2}p^{n}+\frac{1}{4}p^{n-1})
\end{aligned}

This choice of coefficients allows for a factor of two increase in the CFL limit associated with internal gravity waves, and hence a factor of two increase in timestep if internal gravity waves are the controlling factor (see Maltrud et al. (1998), (1998), for a simple proof of this property).

The new pressure at time (n+1) is obtained by updating the tracer transport equations for the temperature \Theta and salinity S before the baroclinic momentum equations are solved. Then the density and hydrostatic pressure at the new time can be diagnosed from the equation of state and the hydrostatic equation. The convective adjustment is done afterwards and this modifies the new \Theta and S and hence the new pressure. So, in this case the new pressure before convective adjustment is used in (73). With implicit vertical mixing the vertical loops over tracer and momentum are separated, and the implicit vertical diffusion of tracers is done after the tracer loop and before the momentum loop, so in that case the exact pressure at the new time is used in (73). Equation (73) applies to leapfrog timesteps. On forward steps the time n pressure is used, and on backward steps the predicted pressure from the first pass is used (see (71), (72)).

5.4.2. Semi-Implicit Treatment of Coriolis and Vertical Friction Terms

A semi-implicit treatment of the Coriolis terms can allow a somewhat longer timestep due to filtering of inertial waves and barotropic Rossby waves, but the main motivation in POP is to damp the Rossby-wave computational mode which appears in the implicit free-surface formulation of the barotropic equations, as shown in Sec. Barotropic Equations and Dukowicz and Smith (1994). Since the barotropic equations are the exact vertical average of the momentum equations, the Coriolis terms in the baroclinic equations must also be treated semi-implicitly. While it is possible to run POP with explicit treatment of the Coriolis terms, we strongly recommend against this because of the above-mentioned computational mode. The following discussion assumes the Coriolis terms are treated semi-implicitly.

The simultaneous semi-implicit treatment of both Coriolis and vertical mixing terms leads to a coupled system where both components of velocity must be solved for simultaneously. To avoid this, we employ an operator splitting which maintains the second-order accuracy of the time discretization. To illustrate this splitting, we write the momentum equations in the matrix form

(74)\begin{aligned}

{\bf u}' - {\bf u}^{n-1}
+ \tau {\sf B}\overline{{\bf u}}^{\alpha\gamma}
&=\tau {\bf F} + \tau {\sf A}(\mu)\overline{{\bf u}}^{\lambda} \\
\tau = 2\Delta t,\
{\bf u} =
\left( \begin{array}{l}
         u_x \\
         u_y \end{array}\right),
&\  {\sf B} =
\left( \begin{array}{lr}
         0 & -f \\
         f & 0   \end{array}\right),
\ {\sf A}(\mu) = \delta_z\mu\delta_z \nonumber\end{aligned}

where {\bf u} is the velocity vector organized in block form and {\sf B} is a 2\times2 matrix in the space of the two velocity components. {\bf F} represents all other terms in the momentum equation that are treated explicitly (advection, metric, pressure gradient and horizontal diffusion terms). The time-averaging of the velocities \overline{{\bf u}}^{\alpha\gamma} and \overline{{\bf u}}^{\lambda} are given by (66) and (67). Equation (74) can be rewritten as

[{\bf I} - \lambda\tau{\sf A}(\mu) + \alpha\tau{\sf B}]
\Delta{\bf u} = \tau \mathbf{F'} \nonumber

\mathbf{F'} \equiv
\{ {\bf F} - {\sf B}[\gamma{\bf u}^n + (1-\gamma){\bf u}^{n-1}]
+{\sf A}(\mu){\bf u}^{n-1}\} \nonumber

(75)\Delta{\bf u} \equiv {\bf u}' - {\bf u}^{n-1}

where {\bf I} is the identity matrix. The operator splitting is given by

(76)\begin{aligned}
[{\bf I} - \lambda\tau{\sf A}(\mu) + \alpha\tau{\sf B}]
= [{\bf I} + \alpha\tau{\sf B}]
  [{\bf I} - \lambda\tau{\sf A}(\mu) ]
+ {\cal O}(\tau^2)
\end{aligned}

and thus second-order accuracy in time is maintained by dropping the {\cal O}(\tau^2) terms and solving the simpler system

(77)\begin{aligned}
[{\bf I} - \lambda\tau{\sf A}(\mu) ]\Delta{\bf u}
= ({\bf I} + \alpha\tau{\sf B})^{-1}\tau{\bf F}'
\end{aligned}

The r.h.s. of (77) can be evaluated analytically using

\begin{aligned}
({\bf I} + \alpha\tau{\sf B})^{-1}
=\frac{1}{1+(\alpha\tau f)^2}({\bf I} - \alpha\tau{\sf B})
\;\;.\end{aligned}

Note again that the surface forcing is contained only on the r.h.s. of (77). Furthermore, the quadratic bottom drag term (50) is evaluated at the lagged time (n-1) and also appears only on the r.h.s. By lagging the bottom drag term, the tridiagonal systems for u_x and u_y are linearized and decoupled, which greatly facilitates their solution. Again, \lambda=0 is the explicit vertical mixing case. For implicit vertical mixing, \lambda > 0.5 is required, and the code is typically run with \lambda = 1. In the barotropic equations (Sec. Barotropic Equations) we choose \alpha=\gamma=1/3, which is hardwired in the code, and it is clear from (66) that the Coriolis terms are centered at time n.

On forward and backward steps the splitting is given by

Forward step:

(78)\begin{aligned}
{\bf u}'^* - {\bf u}^{n} + \tau {\sf B}\overline{{\bf u}}^{\theta}
 & = \tau {\bf F} + \tau {\sf A}(\mu)\overline{{\bf u}}^{\lambda'}
 \nonumber \\
 & \mathrm{or} \nonumber \\
 [{\bf I} - \lambda'\tau{\sf A}(\mu) + \theta\tau{\sf B}]
\Delta{\bf u} & = \tau [{\bf F}^n - {\sf B}{\bf u}^n
+{\sf A}(\mu){\bf u}^{n}]
\nonumber \\
\Delta{\bf u} & \equiv {\bf u}'^* - {\bf u}^{n}
\end{aligned}

Backward step:

(79)\begin{aligned}
{\bf u}' - {\bf u}^{n} + \tau {\sf B}\overline{{\bf u}}^{\theta}
& = \tau {\bf F} + \tau {\sf A}(\mu)\overline{{\bf u}}^{\lambda'}
\nonumber \\
 & \mathrm{or} \nonumber \\
[{\bf I} - \lambda'\tau{\sf A}(\mu) + \theta\tau{\sf B}]
\Delta{\bf u} & = \tau [ {\bf F}^* - {\sf B}{\bf u}^n
+{\sf A}(\mu){\bf u}^{n}] \nonumber \\
\Delta{\bf u} & \equiv {\bf u}' - {\bf u}^{n}
\end{aligned}

where \tau = \Delta t, {\bf F}^n represents the advection and hydrostatic pressure gradient terms evaluated at time n, and {\bf F}^* represents the same terms evaluated using the predicted variables from the forward step. In (78) and (79) \overline{{\bf u}}^{\lambda'} is defined as in (71) and (72), respectively. Now employing the operator splitting

(80)\begin{aligned}
[{\bf I} - \lambda'\tau{\sf A}(\mu) + \theta\tau{\sf B}]
= [{\bf I} + \theta\tau{\sf B}]
  [{\bf I} - \lambda'\tau{\sf A}(\mu) ]
+ {\cal O}(\tau^2)
\end{aligned}

the tridiagonal systems analogous to (77) are:

Forward step:

[{\bf I} - \lambda'\tau{\sf A}(\mu) ]\Delta{\bf u}
= ({\bf I} + \theta\tau{\sf B})^{-1}\tau[{\bf F}^n
- {\sf B}{\bf u}^n +{\sf A}(\mu){\bf u}^{n}]
\nonumber

(81)\Delta{\bf u} \equiv {\bf u}'^* - {\bf u}^{n}

Backward step:

[{\bf I} - \lambda'\tau{\sf A}(\mu) ]\Delta{\bf u}
= ({\bf I} + \theta\tau{\sf B})^{-1}\tau[{\bf F}^*
- {\sf B}{\bf u}^n +{\sf A}(\mu){\bf u}^{n}]
\nonumber

(82)\Delta{\bf u} \equiv {\bf u}' - {\bf u}^{n}

Again, \theta=0.5 is hardwired into the code.

5.5. Barotropic Equations

POP uses the implicit free-surface formulation of the barotropic equations developed by Dukowicz and Smith (1994), and this formulation is presented here. Other possible options are the rigid-lid streamfunction approach (Bryan 1969), the rigid-lid surface pressure approach (Smith, Dukowicz, and Malone 1992; Dukowicz, Smith, and Malone 1993), and the explicit free-surface method (Killworth et al. 1991), which involves subcycling the barotropic mode with a smaller timestep than that used in the baroclinic equations.

The prognostic equation for the barotropic velocity, defined by (70), is obtained by vertically integrating the momentum and continuity equations. The barotropic momentum equation in block form, analogous to (74), is given by:

(83)\begin{aligned}
&&{\bf U}^{n+1} - {\bf U}^{n-1}
+ \tau {\sf B}\overline{\bf U}^{\alpha\gamma}
=
- \tau g\nabla \overline{\eta}^{\alpha\gamma}
+ \tau {\bf F}_B
\nonumber \\
&&\overline{\bf U}^{\alpha\gamma} =
\alpha{\bf U}^{n+1} + \gamma{\bf U}^{n} + (1-\alpha-\gamma){\bf U}^{n-1}
\nonumber \\
&&\overline{\eta}^{\alpha\gamma} =
\alpha\eta^{n+1} + \gamma\eta^{n} + (1-\alpha-\gamma)\eta^{n-1}
\end{aligned}

where \eta = p_s/\rho_{_0} g is the surface height, and {\bf F}_B is the vertical integral of all terms besides the time-tendency, Coriolis, and surface pressure gradient terms in the momentum equation.

5.5.1. Linear Free Surface Model

A prognostic equation for the free surface height \eta is obtained by vertically integrating the continuity equation

(84)\int^{\eta}_{-H}dz(\nabla\cdot{\bf u}+\frac{\partial w}{\partial z}) =
\frac{\partial}{\partial t}\eta + \nabla\cdot(H+\eta){\bf U} - q_{_{W}} = 0,

where we have used the surface boundary condition on the vertical velocity:

(85)w(\eta)=d\eta/dt-q_{_{W}} = \frac{\partial}{\partial t}\eta
+ {\bf u}(\eta)\cdot\nabla\eta - q_{_{W}},

and {\bf U} is the vertically averaged horizontal velocity:

(86){\bf U}=\frac{1}{H+\eta}\int^{\eta}_{-H}dz{\bf u}(z)

where w(\eta) and {\bf u}(\eta) are the vertical and horizontal velocities at the surface, and q_{_W} is freshwater flux. This result is derived using Leibnitz’s Theorem to interchange the order of integration and differentiation. A difficulty with this form of the barotropic continuity equation is that, in the implicit time discretization of the barotropic equations, the term proportional to \eta inside the divergence in (84) leads to a nonlinear elliptic system, and standard solution methods such as conjugate gradient algorithms cannot be directly applied to it. To avoid this, POP uses a linearized form of the barotropic continuity equation which is derived as follows. Integrating the continuity equation over depth as before, but modifying the boundary condition (85) by dropping the term involving \nabla\eta (which can be shown to be of order |\eta|/dz_1 relative to the time tendency term), we find:

(87)\int^{\eta}_{-H}dz(\nabla\cdot{\bf u}+\frac{\partial
w}{\partial z}) =
\frac{\partial}{\partial t}\eta + \nabla\cdot H{\bf U}
+ \int^{\eta}_{0}dz\nabla\cdot{\bf u} - q_{_{W}} = 0

(88)w(\eta)= \frac{\partial}{\partial t}\eta - q_{_{W}}

(89){\bf U}=\frac{1}{H}\int^{0}_{-H}dz{\bf u}(z)

To obtain the linearized barotropic continuity equation we drop the term \int^{\eta}_{0}\!dz\;\nabla\cdot{\bf u} in (87), which corresponds to neglecting the horizontal mass flux between z=0 and z=\eta:

(90)\begin{aligned}
\frac{\partial}{\partial t}\eta + \nabla\cdot H{\bf U} = q_{_{W}}
\;\;.
 \end{aligned}

This derivation makes it clear that in the advection operators the horizontal mass flux between z=0 and z=\eta must be neglected in order to be consistent with (90). So there are four ingredients to the linear free surface model: 1) the barotropic continuity equation is given by (90) instead of (84); 2) the barotropic velocity is given by (89) instead of (86); 3) the vertical velocity at the surface, which is used to integrate the continuity equation from the top down, is given by (88) instead of (85); and 4) the horizontal mass fluxes between z=0 and z=\eta should be neglected in the advection operators and when integrating the continuity equation to find the vertical velocities. In the discrete equations this means that the horizontal mass fluxes in the surface cells are proportional to the full cell height dz_1, rather than dz_1+\eta (see Sec. Variable-Thickness Surface Layer). This linear approximation is valid provided the surface displacement is small compared to the depth of the full ocean: |\eta|\ll H, and in the discrete equations the surface displacement must also be small compared to the depth of the surface layer: |\eta|\ll dz_1.

With a non-zero freshwater flux q_{_{W}} the mean volume of the ocean is not constant, even though the velocity field is divergence free. Integrating (90) over horizontal area, we find

(91)\begin{aligned}
\frac{d}{dt}\int\!dV = \int\!da\; \frac{\partial}{\partial t}\eta
=\int\!da\; q_{_{W}}\;
\end{aligned}

where da = dx dy is the horizontal area element. So the total volume of the ocean will change unless the surface-integrated freshwater flux vanishes.

5.5.2. Time Discretization of the Barotropic Equations

In the implicit free-surface formulation the continuity equation (90) is discretized in time using a forward step as follows:

\frac{\eta^{n+1}-\eta^{n}}{\Delta t}
+ \nabla\cdot H\overline{\bf U}^{\alpha'} = q^n_{_W}
\nonumber

(92)\overline{\bf U}^{\alpha'}
=\alpha'{\bf U}^{n+1} + (1-\alpha'){\bf U}^{n}

where \alpha' is a coefficient used to vary the time centering of the velocity in the continuity equation. The barotropic equations support three types of linear waves: two gravity waves (one in each horizontal direction) and one Rossby wave. In a pure leapfrog discretization, there would be three computational modes, one associated with each of these waves. By choosing the forward discretization (92), one computational mode is eliminated, leaving one gravity-like and one Rossby-like computational mode. The Rossby computational mode is damped by the discretization scheme if the Coriolis terms are treated implicitly, and the gravity-wave computational mode is strongly damped if \alpha' is close to 1. In Dukowicz and Smith (1994) it was shown that the optimal set of time-centering coefficients which maximally damps the computational modes, minimally damps the physical modes, and minimally distorts the phase speed of the physical modes is given by the parameter set:

(93)\alpha=\gamma=\frac{1}{3},\alpha'=1.

Thus in (83)

\overline{\eta}^{\alpha\gamma} = \frac{1}{3}
(\eta^{n+1}+\eta^{n}+\eta^{n-1}) \nonumber

and similarly for \overline{\bf u}^{\alpha\gamma}. These choices are hardwired into the code. The physical gravity waves are damped at small space and time scales in this implicit scheme, but the physical Rossby waves are essentially unaffected.

By inserting the barotropic momentum equation (83) into the continuity equation (92) we obtain an elliptic equation for the surface height at the new time \eta^{n+1}. However, due to the presence of the Coriolis terms, the resulting elliptic operator is not symmetric, making it much more difficult to invert, because standard solvers such as conjugate gradient methods require symmetric positive-definite linear operators. This problem can be overcome by using an operator splitting technique that maintains the second-order accuracy of the time discretization scheme (Smith, Dukowicz, and Malone 1992; Dukowicz, Smith, and Malone 1993). Defining an auxiliary velocity

(94)\hat{\bf U}\equiv{\bf U}^{n+1}
+ \alpha\tau g\nabla(\eta^{n+1}-\eta^{n-1})

equation (83) can be written:

(95)\begin{aligned}
({\bf I} + \alpha\tau{\sf B})(\hat{\bf U} - {\bf U}^{n-1})
& =\tau\{{\bf F}_B -
{\sf B}[\gamma{\bf U}^n + (1-\gamma){\bf U}^{n-1}] \nonumber \\
& - g\nabla[\gamma\eta^n + (1-\gamma)\eta^{n-1}]\}
+{\cal O}(f\tau^3)
\end{aligned}

Dropping the {\cal O}(f\tau^3) terms (which are the same order as the time truncation error in this second-order scheme), and using the continuity equation (92) with \alpha'=1, we arrive at the elliptic system

(96){\big [} \nabla\cdot H\nabla - \frac{1}{g\alpha\tau\Delta t}
\big] \eta^{n+1} =
\nabla\cdot H\big[\frac{\hat{\bf U}}{g\alpha\tau}
+\nabla\eta^{n-1} \big] - \frac{\eta^{n}}{g\alpha\tau\Delta t}
-\frac{q^n_{_W}}{g\alpha\tau}

The procedure is to first solve (95) for \hat{\bf U}, then solve (96) for \eta^{n+1} (POP actually solves for the surface pressure p_s = \rho_{_0} g\eta), and finally use (94) to obtain {\bf U}^{n+1}. This system is solved in POP on leapfrog timesteps with \alpha=1/3 using a diagonally preconditioned conjugate gradient algorithm described in Sec. Conjugate Gradient Algorithm. Note that the terms dropped in (95) are only {\cal O}(\tau^3) if the timestep is small compared to the inertial period 1/f. We therefore recommend that the model timestep be no greater than about two hours. The scheme is stable for longer timesteps, but the barotropic mode will start to become inaccurate as the timestep is increased above two hours. The divergence on the r.h.s. in (96) is discretized with the correct B-grid discretization, i.e., the quantity in brackets is transversely averaged as in (33). The Laplacian-like operator on the l.h.s. in (96) is a nine-point stencil with the correct B-grid discretization:

(97)\nabla\cdot H\nabla\eta =
\frac{1}{\Delta_y} \delta_x
\overline{\big[\Delta_y H_U\delta_x \overline{\eta}^y \big]}^y
+ \frac{1}{\Delta_x} \delta_y
\overline{\big[\Delta_x H_U\delta_y \overline{\eta}^x \big]}^x.

Unlike the friction and diffusion operators, this operator cannot be approximated as a five-point stencil, because doing so would violate the energetic balance between pressure work and change in potential energy (see Sec. Energetic Consistency). Since the nine-point B-grid operator in (97) has a checkerboard null space (see Sec. Nullspace Removal), the solution of (96) is prone to have local patches of checkerboard noise. Strictly speaking, the operator on the l.h.s. of (96) has no checkerboard null space due to the presence of the extra diagonal term (second term in brackets). However, if the solution in some region is in a near steady state, this diagonal term is cancelled by the next to last term on the r.h.s., and checkerboard noise may appear. This is particularly true in isolated bays where the solution is only weakly coupled to the interior. On the other hand, the checkerboard noise has little effect on the dynamics because only the gradient of surface height enters the barotropic momentum equation, and the B-grid gradient operator (33) does not see a checkerboard field. The only way a checkerboard surface height field can affect the dynamics is through the vertical velocity at the surface (39), which depends on the change in surface height, and is used as the surface boundary condition to integrate the continuity equation to find the advection velocities. However, experience has so far shown that this does not lead to serious problems with the model simulations.

It should be noted that the solution for surface height will only satisfy the continuity equation (92) to the extent that the solution of (96) has converged in the iterative solution of the elliptic solver. As discussed in Sec. Tracer Advection, this can lead to a small non-divergent mass flux in the bottom-most ocean cell when the 3-D continuity equation is integrated from the top down with w given by (88) at the surface, as discussed in Sec. Tracer Advection. In practice, we suggest that the convergence criterion for the iterative solver be chosen so that the global mean balance between pressure work and change in potential energy is accurate to within about three significant digits (see Sec. Energetic Consistency on energetic consistency).

The time discretization for forward and backwards steps, corresponding to (83) and (92) are given by:

Forward step:

{\bf U}^* - {\bf U}^{n}
+ \tau {\sf B}\overline{\bf U}^{\theta} =
- \tau g\nabla \overline{\eta}^{\theta}
+ \tau {\bf F}_B^n \nonumber

\frac{\eta^*-\eta^{n}}{\Delta t}
+ \nabla\cdot H{\bf U}^* = q^n_{_W} \nonumber

\overline{\bf U}^{\theta} =\theta{\bf U}^* +
(1-\theta){\bf U}^{n} \nonumber

(98)\overline\eta^{\theta} =\theta\eta^* + (1-\theta)\eta^{n}

Backward step:

{\bf U}^{n+1} - {\bf U}^{n}
+ \tau {\sf B}\overline{\bf U}^{\theta} =
- \tau g\nabla \overline{\eta}^{\theta}
+ \tau {\bf F}^*_B \nonumber

\frac{\eta^{n+1}-\eta^{n}}{\Delta t}
+ \nabla\cdot H{\bf U}^{n+1} = q^n_{_W}
\nonumber

\overline{\bf U}^{\theta} =\theta{\bf U}^{n+1} +
(1-\theta){\bf U}^{n} \nonumber

(99)\overline\eta^{\theta}
=\theta\eta^{n+1} + (1-\theta)\eta^{n}

where we have assumed \alpha' = 1 as in (92) and (93). {\bf U}^* and \eta^* are the predicted barotropic velocity and surface height from the forward step. {\bf F}^n_B contains all terms other than the Coriolis, surface pressure gradient and time-tendency terms in the barotropic momentum equation, all evaluated at time n, and {\bf F}^*_B contains the same terms but evaluated using the prognostic variables predicted by the forward step. The operator splittings for the forward and backward steps, analogous to (94) and (95) are:

Forward step:

\hat{\bf U}\equiv{\bf U}^*
+ \theta\tau g\nabla(\eta^*-\eta^{n}) \nonumber

(100)({\bf I} + \theta\tau{\sf B})(\hat{\bf U} - {\bf U}^{n})
=\tau \{{\bf F}^n_B - {\sf B}{\bf u}^n - g\nabla\eta^n \}
+{\cal O}(f\tau^3)

Backward step:

\hat{\bf U}\equiv{\bf U}^{n+1}
+ \theta\tau g\nabla(\eta^{n+1}-\eta^*)
\nonumber \\

(101)({\bf I} + \theta\tau{\sf B})(\hat{\bf U} - {\bf U}^{n})
=\tau\{{\bf F}^*_B - {\sf B}{\bf u}^n
- g[\theta\nabla\eta^* +(1-\theta)\nabla\eta^n]\}
+{\cal O}(f\tau^3)

Finally, the elliptic equations for the forward and backward steps analogous to (96) are given by:

Forward step:

(102){\big [} \nabla\cdot H\nabla - \frac{1}{g\theta\tau\Delta t}
\big] \eta^* =
\nabla\cdot H\big[\frac{\hat{\bf U}}{g\theta\tau}
+\nabla\eta^{n} \big] - \frac{\eta^{n}}{g\theta\tau\Delta t}
-\frac{q^n_{_W}}{g\theta\tau}

Backward step:

(103){\big [} \nabla\cdot H\nabla - \frac{1}{g\theta\tau\Delta t}
\big] \eta^{n+1} =
\nabla\cdot H\big[\frac{\hat{\bf U}}{g\theta\tau}
+\nabla\eta^* \big] - \frac{\eta^{n}}{g\theta\tau\Delta t}
-\frac{q^n_{_W}}{g\theta\tau}

5.5.3. Conjugate Gradient Algorithm

The most efficient method we have found for solving the elliptic system of barotropic equations (96) is to use a standard Preconditioned Conjugate Gradient solver. The algorithm consists of the following steps to solve the system Ax = b given by (96) multiplied by the T-cell area. A is the symmetric positive-definite Laplacian type operator on the l.h.s. in (96) (it is important to note that if (96) is not multiplied by the T-cell area then the operator is not symmetric), b is the r.h.s., and x is the solution. Lower-case Greek and Roman letters are used, respectively for scalars and 2-D arrays, and (x,y) denotes a dot product: (x,y) = \sum_{ij} x_{ij} y_{ij}.

PCG Algorithm:

Initial guess: x_o

Compute initial residual r_o = b-Ax_o

Set \beta_o=1, s_o=0

For k=1,2,...,k_{max} do

\begin{aligned}
r'_{k-1} &=& Zr_{k-1} \nonumber \\
\beta_k &=& (r_{k-1},r'_{k-1}) \nonumber   \\
s_k &=& r'_{k-1} + (\beta_k/\beta_{k-1})s_{k-1} \nonumber   \\
s'_k &=& A s_k   \\
\alpha_k &=&\beta_k/(s_k,s'_k) \nonumber   \\
x_k &=& x_{k-1} + \alpha_k s_k  \nonumber  \\
r_k &=& r_{k-1} - \alpha_k s'_k
\end{aligned}

Exit if converged:

(r_k,r_k)^{\frac{1}{2}} < \epsilon\overline{a}

End do

Here Z is a preconditioning matrix. This is usually taken to be a diagonal preconditioner Z = C_0^{-1}, where C_0 is the diagonal matrix composed of the central coefficient in the nine-point stencil corresponding to the operator in (96). The code is set up to use a more sophisticated preconditioner consisting of a nine-point operator which is a local approximate inverse of the true operator (Smith, Dukowicz, and Malone 1992). This preconditioner is expensive to compute but is time independent and can be computed offline and read in by the code. To improve the initial guess, x_o is linearly extrapolated in time from the solutions at the two previous timesteps. In the convergence criterion (r_k,r_k)^{\frac{1}{2}} < \epsilon\,\overline{a}, \overline{a} is the rms T-cell area. This normalization is somewhat arbitrary, but is chosen so that \epsilon has the same dimensions as the operator (96) before it is multiplied by the T-cell area averaged over surface ocean points. As discussed in Sec. Time Discretization of the Barotropic Equations, we suggest that \epsilon be chosen such that the global mean balance between pressure work and change in potential energy is accurate to within about three significant digits. Typically this is achieved for values of \epsilon between 10^{-12} and 10^{-13}.

The conjugate gradient solver outlined above has poor scalability on many distributed memory computer systems. The fundamental limit to scalability is the need for three dot products, each requiring a global summation operation and adding an inherent sychronization. A variant of the standard conjugate gradient solution method, known as the Chronopoulos-Gear algorithm (D’Azevedo, Eijkhout, and Romine 1993), combines the separate global sums into a single vector global sum, which reduces the latency to one third. This option is now the default in POP, though the previous PCG algorithm is still supported.

6. Grids

6.1. Global Orthogonal Grids

As discussed in Sec. Spatial Discretization, the POP model is designed to handle any orthogonal curvilinear coordinate system for the horizontal grid. The only information the code needs locally at each grid point is the longitude and latitude of velocity points and the distances to neighboring points along the two coordinate directions. This information is computed offline and read in from a file. However, there are restrictions on the global connectivity of the grid. In the current model release, it is assumed that the computational grid can be mapped onto a 2-D domain. When such global grids are mapped onto the sphere they necessarily involve two singularities, analogous to the north and south poles on a standard polar grid. In these ``dipole” grids, originally developed independently by two groups (Madec and Imbard 1996; Smith, Kortas, and Meltz 1995), the two singularities can be placed inside land masses (e.g., North America, Greenland or Asia) to construct global grids. POP also has the ability to use what we call “tripole” grids, which were originally developed by Murray (1996). In these grids, there are two “poles” in the northern hemisphere (usually placed in North America and Asia). These grids can also be mapped onto a cylinder but special boundary communications are required in order to “sew up” the line between the two northern poles. In both dipole and tripole grids, the meshes are constructed so that there is a smooth transition moving north from the Equator, from the southern hemisphere polar grid (usually Mercator grid) to the northern “distorted” grid. Various factors enter into the design of a given grid, and different choices (pole locations, dipole vs. tripole) may be preferred for different problems. In global grids the Equator is typically chosen as a grid line so that the numerical dispersion relation for the B-grid linear Kelvin waves is accurate. At resolutions coarser than about \frac{1}{2}^\circ, the meridional spacing near the Equator is enhanced (to about \frac{1}{2}^\circ) in order to more accurately resolve equatorial waves and currents. Note: Software for generation of global dipole and tripole grids is not released with the code, but will eventually be made public.

6.1.1. Dipole Grids

The simplest orthogonal grid with two arbitrarily located poles can be constructed analytically in one of three ways:

  1. Electrostatic dipole charge distribution on a sphere: solve Poisson’s Equation on the sphere with two electric charges of opposite sign at the chosen pole points; the electric field lines and the lines of constant electrostatic potential form an orthogonal coordinate system;
  2. Conformal map method (Murray 1996): project a spherical polar grid onto a tangent plane at its equator from the point on the sphere opposite the tangent point, then reproject back onto a larger sphere tangent to the plane at the same point; the new pole points on the larger sphere will be shifted toward its equator and the tangent point;
  3. Geometric method (Smith, Kortas, and Meltz 1995): choose two axes, one which intersects the surface of the sphere at the two pole points and one located at the intersection of the two planes tangent to the pole points, then the two families of planes that contain these axes intersect the surface of the sphere in two sets of circles which form the coordinates of an orthogonal coordinate system.

All these methods produce equivalent orthogonal grids. The only problem with these simple analytic grids is that if one pole is located in Antarctica and the other is displaced away from the true North Pole, then the Equator will not be a grid line. To get around this problem, grids can be constructed by an iterative procedure (Madec and Imbard 1996; Smith, Kortas, and Meltz 1995), starting at the equator of a southern hemisphere polar grid, where the northern pole point is gradually shifted away from the true North Pole as successive latitude-like circles are constructed. In this manner it is possible to construct a semi-analytic dipole grid which smoothly matches onto the southern polar grid, and in which local changes in grid spacing are everywhere small and smooth. Composite grids that have sudden changes in grid spacing normal to matching boundaries suffer a loss of formal second-order accuracy of the spatial discretization scheme (Smith, Kortas, and Meltz 1995).

An example of a semi-analytic dipole grid is shown in Figure5.1 . In this grid the northern pole is in North America and the southern hemisphere is a Mercator grid with a pole at the true South Pole. Note that the meridional grid spacing near the Equator has been enhanced as discussed above. The grid spacing in the northern hemisphere is chosen so as to minimize the global mean aspect ratio, attempting to make each cell as nearly “square” as possible (as in a true Mercator grid). In both hemispheres the grid terminates in a circle enclosing the pole, and anything within this circle is excluded (e.g., Hudson Bay in Figure5.1.

../_images/FigDipoleGrid.jpg

Figure 5.1 A dipole grid.

6.1.2. Tripole Grids

Analytic grids can also be constructed which involve more than two singularities using the different methods described above for dipole grids. Various types of multipole grids are discussed by Murray (1996). One variation that is particularly useful for ocean models is shown in Figure 11 of Murray (1996). It has two singularities in the northern hemisphere and smoothly matches onto a southern hemisphere Mercator grid. A “tripole” grid of this type is shown in Figure5.2 . Note that the grid spacing in the Arctic is much more uniform and the cell aspect ratios are much closer to one than in the dipole grid in Figure5.1 . This particular grid was constructed for the French OPA model (Madec, Imbard, and Levy 1998), and was made available, courtesy of G. Madec, for testing in POP. Like the dipole grid, the tripole grid is mapped onto a cylinder (a logical 2-D array that is periodic in the x direction). The difference is that in the dipole grid the northernmost (``upper”) boundary of the grid corresponds to the circle surrounding the pole (e.g. the open circle containing Hudson bay in Figure5.1 ), whereas in the tripole grid the two northern poles both lie on the upper boundary, one is located at the two upper corners of the logical grid (which correspond to the same point because it is periodic) and the other at the midpoint of the upper boundary. When the 2-D cylindrical array is mapped onto the sphere (as in Figure5.2 ) there is a chopped pole in the southern hemisphere, but in the northern hemisphere there is a “cut” in the grid along the line between the poles, and flow across this line from one side of the Arctic to the other requires non-local communication along the upper boundary of the logical grid: it must be folded back along itself at its mid point in order to make neighboring cells in the physical grid be adjacent in the logical grid. The indexing for this nonlocal communication is different for quantities located at different places (i.e. U-points or T-points) on the grid. POP does assume that in the logical grid the upper boundary lies on the northern edge of the last row of T-cells, so that U-points lie exactly on the upper grid boundary, along the line between the two poles. Furthermore, it assumes that the middle pole point lies at a U-point exactly halfway across the top of the grid.

../_images/FigTripoleGrid.jpg

Figure 5.2: A tripole grid.

6.2. Vertical Grid Variants

The current vertical grid in POP is a fixed Eulerian grid with depth as the vertical coordinate. However, two options are available to mitigate some problems encountered with such a grid. A variable-thickness surface layer can be used to more accurately treat changes in sea surface height and fluxes of tracers at the surface. Partial bottom cells can be used to more accurately represent bottom topography.

6.2.1. Variable-Thickness Surface Layer

The tracer transport equations have the form

(104)\frac{\partial}{\partial t}\varphi
+ \nabla\cdot{\bf u}\varphi
+ \frac{\partial}{\partial z}w\varphi
= \nabla\cdot {\bf F} + \frac{\partial}{\partial z} F_V

where \nabla\cdot {\bf F} = {\cal D}_H and \frac{\partial}{\partial z} F_V = {\cal D}_V are, respectively, the horizontal and vertical diffusion operators. Integrating this equation over the model surface level (from -h_1 to \eta, where h_1 = dz_1 is the depth of the upper level) we find:

(105)\begin{aligned}
\frac{\partial}{\partial t}\int^{\eta}_{-h_1}dz\varphi
+ \nabla\cdot\int^{0}_{-h_1}dz{\bf u}\varphi
-w(h_1)\varphi(h_1) = \nonumber \\
q_{_{W}}\varphi(\eta) + \nabla\cdot\int^{0}_{-h_1}dz{\bf F}
+F_V(\eta) - F_V(h_1)
\end{aligned}

In deriving (105) the same approximations used in the linear free surface model (Sec. Linear Free Surface Model) were employed: the boundary condition (88) was applied, and the advective and diffusive horizontal fluxes between z=0 and z=\eta were set to zero. Note that there is no advection of tracers through the surface with the vertical velocity w(\eta).

Some care must be taken to specify the tracer fluxes through the air-sea interface. The total tracer flux seen by the ocean model at the surface is given by

(106)\begin{aligned}
F_T = q_{_{W}}\varphi(\eta) + F_V(\eta)\;,
\end{aligned}

where \varphi(\eta) is the tracer concentration at the sea surface and q_{_{W}}\varphi(\eta) represents the advection of tracers in the ocean relative to the sea surface due to the freshwater flux. The atmosphere (or sea ice model) sees a flux of tracers into the ocean given by

(107)\begin{aligned}
Q_T = q_{_W}\varphi_{_W} + Q_\varphi\;,
\end{aligned}

where \varphi_{_W} is the tracer concentration in the freshwater being added or removed from the ocean, (for simplicity we assume here that locally there will be only one source of freshwater, and omit the sum over different types freshwater sources – see (36)), and Q_{\varphi} is any additional flux of tracers (e.g., for temperature in the open ocean it corresponds to the sum of sensible, latent, shortwave and longwave heat fluxes). In reality there is a very thin boundary layer where the tracer concentration varies continuously and \varphi(\eta)=\varphi_{_W} at the interface. The models do not resolve this boundary layer and so in general \varphi(\eta)\ne\varphi_{_W}. However, the total flux must be conserved across the boundary layer, so that F_T = Q_T, and this can be used to eliminate \varphi(\eta) in terms of \varphi_{_W}. Substituting (107) for (106) in (105), we obtain

(108)\begin{aligned}
\frac{\partial}{\partial t}(h_1+\eta)\varphi
+ h_1\nabla\cdot{\bf u}_1\varphi_1
-w_1\varphi(h_1)
= \nonumber \\ h_1\nabla\cdot {\bf F}_1
+Q_{\varphi} - F_V(h_1) + q_{_{W}}\varphi_{_W}
\end{aligned}

where now the flux Q_\varphi is applied as the surface boundary condition on {\cal D}_V. The integrals over the surface layer have been evaluated by assuming the model variables are constant within the surface layer. w_1 is the velocity at the bottom of the surface layer and w_1\varphi(h_1) represents tracer advection through the bottom of the surface layer. Dividing (108) by h_1, we arrive at the tracer transport equation (34).

It remains to specify the freshwater tracer concentrations \varphi_{_W}. In the case of salinity the default choice in the code is \varphi_{_W} = S_{_W} = 0, that is, the freshwater flux q_{_{W}} = P - E + R - F_{ice} + M_{ice} is assumed to have zero salinity (where P, E, R, F_{ice} and M_{ice} are the freshwater fluxes associated with precipitation, evaporation, river runoff, freezing and melting of sea ice). In the case of potential temperature, the default is to assume the freshwater has the same temperature as the model surface layer, so that \varphi_{_W} = \Theta_{_W} = \Theta_1. However, in general the tracer concentration in the freshwater may vary depending on its source (precipitation, evaporation, etc.), and the default assumption that these are all given by a single value \varphi_{_W} as in (107) may be inadequate. When coupling to atmospheric or sea-ice models a more accurate accounting of the fluxes may be needed to ensure they are conserved between component models.

Combining (108) with the transport equations in the subsurface levels and integrating over volume, we find

(109)\begin{aligned}
\frac{d}{dt}\int\!dV\;\varphi = \int\!da\; (Q_{\varphi} +
q_{_{W}}\varphi_{_W})\;.
\end{aligned}

Thus tracers are conserved in the absence of surface fluxes. This is also true of the time-discretized equations in the model. In particular, if there are no surface fluxes of salinity, Q_S = 0, and if the freshwater flux has zero salinity (\varphi_{_W} = S_{_W} = 0), then total salt is conserved.

The time discretization of (108) or (34) is given by (51) for leapfrog steps and by (54) and (55) for forward and backward (Matsuno) timesteps. All of these discrete equations conserve global mean tracers exactly in the sense of (109).

There is also a small error in the time discretization associated with the fact that we have adopted a three-time-level leapfrog scheme for the tracer transport equations (51), but a two-time-level discretization for the barotropic continuity equation (92). Note that in the limit of a constant tracer \varphi=1 the continuous transport equation (104), in the absence of surface forcing, reduces to the continuity equation \nabla\cdot{\bf u}+\frac{\partial
w}{\partial z}=0, and the vertically integrated transport equations also reduce to the barotropic continuity equation (90). This does not hold for the discrete equations because of the inconsistency between the two- and three-time-level schemes. The consequence of this is that a tracer which is initially constant in space may evolve to have some small spatial variations. However, this error is second-order in time and is typically extremely small. In a test problem involving the Goldsborough-Stommel circulation (Griffies et al. 2001), which involves only freshwater forcing at the surface, the temperature field, which is initially constant, deviates from its initial value by only 1 part in 10^7 in a 100-year integration, and furthermore, the error does not accumulate in time. Therefore, we feel that this error is small enough that it will not cause problems in long-term climate simulations.

Finally, note that the time discretization error when momentum acceleration is used, see Sec. Tracer Acceleration, may lead to a computational instability when the variable thickness surface layer option is used with explicit fresh water fluxes. Because of this possibility, we advise using the variable thickness surface layer option with virtual salt fluxes if momentum acceleration, with \alpha \ne 1, is being used.

6.2.1.1. Nullspace Removal

The operator for the implicit solution of the barotropic equations for the surface height \eta in (96) contains the Laplacian-like term (97). On the B-grid this nine-point Laplacian-like operator contains two null eigenvectors: a constant field and a global checkerboard (+/- or black/white) field (a checkerboard field has the value of +1 in black cells and -1 in white cells). The Laplacian operator on the B-grid annihilates both these types of null fields (i.e., they have zero eigenvalue). On the C-grid, the global checkerboard is not a null eigenvector; for this reason all the horizontal diffusive operators in POP have an approximate C-grid discretization, allowing grid-scale checkerboard noise to be damped by these operators. However, the presence of a checkerboard component in the surface height field does not significantly affect the dynamical solution in most cases. The dominant effect of the surface height is through its contribution to the surface pressure gradient, and on the B-grid the gradient operator also does not see the checkerboard field. There is only a small indirect effect on the vertical velocities in T-columns because they are determined by vertically integrating the continuity equation with the surface boundary condition (39), which depends on \eta.

In the implicit free surface formulation, the extra term on the left in (96) is a diagonal term associated with the free surface (the \partial_t\,\eta term in the barotropic continuity equation). Without this diagonal term, as in the case in the rigid-lid surface pressure formulation (Smith, Dukowicz, and Malone 1992), the null eigenvectors are not removed by the operator, and can creep into the iterative solution. The diagonal term shifts the eigenvalue of both constant and checkerboard fields away from zero so that they are damped. Therefore, in most cases it is not necessary to explicitly remove constant or checkerboard fields from the solution of the implicit free-surface system. However, experience has shown that the solution sometimes does develop large but localized regions where checkerboard noise appears in the surface height. The exact nature of such regional checkerboarding is not completely understood. It can even occur in the rigid-lid streamfunction formulation using an approximate five-point operator for the barotropic operator which has no global checkerboard nullspace. Note that a five-point approximation is allowed in this case but not in the surface pressure or free-surface formulations, because the five-point approximation to the nine-point B-grid operator violates energetic consistency (Smith, Dukowicz, and Malone 1992). Experience has also shown that this regional checkerboarding does not affect the dynamical solution, again because it is almost completely invisible to the gradient operator. However, with the variable-thickness surface layer (Sec. Variable-Thickness Surface Layer), such regional checkerboarding can appear in the solution and then slowly accumulate over time, eventually leading to large checkerboard components that can affect the solution, either by contributing a large roundoff error to the pressure gradient, or by becoming so large that in some cells the surface height descends below the bottom of the surface layer.

When complicated lateral boundaries are included in the computational mesh, other global null eigenvectors besides the constant and checkerboard fields can also exist. These occur when the domain contains quasi-isolated bays that are only connected to the rest of the ocean through openings of one or two cell faces. These regions can have their own constant and checkerboard nullspaces. There are two types of quasi-isolated bays that need to be considered: 1) bays that are connected to the rest of the ocean at the surface by only one face of a T-cell, and 2) bays that are connected by two adjacent faces of the same T-cell, so that the mouth of the bay is across the diagonal of the T-cell. If either of these types of bay is in the grid, a checkerboard nullspace component of the solution can accumulate in the bay. In the first type of bay, a constant nullspace can also accumulate. In these smaller areas, the problems mentioned above (leading to roundoff error or complete emptying of surface cells) accumulate more quickly in time. Normally, bays of the first type do not appear, because it is the usual practice during construction of a grid to eliminate all points that are connected to the open ocean by only one cell face (such points have no influence on the flow except through diffusion across the open face). But when the variable-thickness surface layer is used, it is necessary to remove all isolated bays of both types in order for the model to run stably. Once these bays are eliminated from the grid, it is still necessary to remove the global constant and checkerboard nullspaces associated with the full domain, to prevent them from slowly accumulating over time. The null fields are projected out of the solution subject to two constraints: 1) the resulting field has no checkerboard component; and 2) the projection leaves the global mean sea level unchanged. Denoting the null eigenvectors {\bf P}_{const} and {\bf P}_{check}, which are normalized to 1 and \pm1, respectively, the nullspaces can be projected out with the following operation (here {\bf X} is the surface height field considered as a one-dimensional vector, and {\bf X}' is the same field after the projection):

(110)\begin{aligned}
{\bf X}'& =& {\bf X} + a\,{\bf P}_{const} - b\,{\bf P}_{check} \nonumber \\
a& =& \frac{X_{check}V_{check}}{N\,V_{const} - M\,V_{check}} \nonumber \\
b& =& \frac{X_{check}V_{const}}{N\,V_{const} - M\,V_{check}} \nonumber \\
X_{check}& =& \sum \, X\cdot{\bf P}_{check}\;, \nonumber \\
N& =& \sum \, {\bf P}_{const}\cdot{\bf P}_{const}\;, \nonumber \\
M& =& \sum \, {\bf P}_{check}\cdot{\bf P}_{check}\;, \nonumber \\
V_{check}& =& \sum\,d{\bf a}\cdot{\bf P}_{check}\;, \nonumber \\
V_{const}& =& \sum\,d{\bf a}\cdot{\bf P}_{const}\;,
\end{aligned}

where d{\bf a} is the one-dimensional vector of cell areas, and the sums are over all surface cells. It is straightforward to show that \sum \,{\bf X}'\cdot{\bf P}_{check} = 0 and \sum \,{\bf X}'\cdot d{\bf a} = \sum \,{\bf X}\cdot d{\bf a}. The latter constraint ensures the projection does not change the global mean sea level. This projection is applied once each timestep immediately after the barotropic equations are solved for \eta.

6.2.2. Partial Bottom Cells

To better represent the bottom topography, we employ the idea of partial bottom cells introduced by Adcroft, Hill, and Marshall (1997) and later developed for the MOM model by Pacanowski and Gnanadesikan (1998), for which the thickness of bottom cells is allowed to vary in space. When we place tracer or pressure points at the center of partial bottom cells, as shown in Figure5.3 , the finite difference for horizontal derivatives between two adjacent points introduces an error because these points are located at different vertical positions. More specifically, when we compute the horizontal pressure gradient, we have

(111)\begin{aligned}
     \left.\frac{\partial p}{\partial x}\,\right\vert_{z=-h(x)}
     = \frac{\partial}{\partial x}\left [\, p\,\vert_{z=-h(x)}\right ]
     - g\rho\,\vert_{z=-h(x)}\, h(x)\,,
     \end{aligned}

where h(x) can be regarded as the position of tracer points. Without the second term on the right hand side of (111), the pressure gradient does not vanish even for \rho=\rho(z) for which we expect p_x\vert_{z=-h(x)}=0. Even though we include the correction term, p_x vanishes only when \rho(z) varies linearly in depth, because p is computed from the hydrostatic equation p_z=-\rho g by using the trapezoidal rule. Likewise, when we compute the horizontal diffusion of tracers without taking the effects of different vertical positions of tracer points into account, we also introduce spurious diffusion errors. To avoid these errors, the interpolation between two vertical points has been used in z-coordinate ocean models. More specifically, after finding {\overline T}_{i,k} by linear interpolation between T_{i,k} and T_{i,k-1},

(112){\overline T}_{i,k} =
        T_{i,k} + \frac{\min({h_z}_{i,k},{h_z}_{i+1,k})-{h_z}_{i,k}}
                       {{h_z}_{i,k}-{h_z}_{i,k-1}}
     \Big (T_{i,k} - T_{i,k-1}\Big ),

where h_z is the height of T-cells, the horizontal gradient is computed as \Delta_x {\overline T}= T_{i+1,k}-
{\overline T}_{i,k}, as shown in Figure5.4 . This approach seems to be correct but it does not guarantee the negative definiteness of tracer variance: for example, for the Laplacian horizontal mixing,

(113)\sum_{i,k}\,
     \left (\Delta_x^2 {\overline T}_{i,k}\right )\,T_{i,k}\,,

is not sign-definite! A correct approach is to work with the functional which can be written as

(114){\cal G}[T]=-\frac{1}{2}\sum\frac{h_yh_z}{h_x}\Big (\Delta_x
     {\overline T}\Big)^2\,,

where h_x and h_y are the horizontal grid sizes of T-cells in the x and y directions, respectively. Taking the functional derivative leads to the following discretization for the horizontal Laplacian diffusion:

\begin{aligned}
%        \frac{\delta {\cal G}}{\delta T} =
        dV_T\,D[T]= \Delta_x\left [\,\left(\frac{h_yh_z}{h_x}\right )\,
        {\overline \kappa}^x\, \Delta_x \overline T\,\,\right ]
        + \Delta_x\left [\left (\frac{h_y}{h_x}\right )
        \Delta_z\left (\gamma h_z{\overline \kappa}^x\Delta_x
        \overline T\,\right )\,\right ]\nonumber\end{aligned}

(115)\begin{aligned}
     \hspace{2cm}
     +\Delta_y\left [\,\left(\frac{h_xh_z}{h_y}\right )\,
        {\overline \kappa}^y\, \Delta_y\overline T\,\,\right ]
        + \Delta_y\left [\left (\frac{h_x}{h_y}\right )
        \Delta_z\left (\gamma h_z{\overline \kappa}^y\Delta_y
        \overline T\,\right )\,\right ]\,.
        \end{aligned}

../_images/FigPBC-PT.jpg

Figure 5.3: Pressure and tracer points in the x-z plane above partially-filled bottom cells.

../_images/FigPBCInterp.jpg

Figure 5.4: Interpolating tracer points to the same level before taking horizontal derivatives.

Since the correct discretization even for the simple Laplacian mixing is too complicated and thus expensive to employ, we have made the partial-bottom-cell implementation in POP in the following ways: (1) we leave the pressure points at the same z-levels which introduces no pressure error and guarantees energetic consistency. Even though some pressure points are underneath the ocean floor, the pressure is computed by extrapolating tracers, therefore density, to pressure points. (2) when we compute the tracer gradients in diffusion operators, we make no interpolation for tracers. This approach introduces spurious diffusion errors but guarantees at least the sign-definite tracer variances, which has been given up in other OGCMs. Our test results show the error introduced by the second approach is so small that we can avoid any complication involving the interpolation.

6.2.2.1. Discretization

After reading the height of bottom T-cells, POP uses two additional three-dimensional arrays for the thickness of T-cells, {\rm DZT}_{i,j,k}, and the thickness of U-cells defined as

(116)\mathrm{DZU}_{i,j,k} =
   \min(\mathrm{DZT}_{i,j  ,k}, \mathrm{DZT}_{i+1,j  ,k},
        \mathrm{DZT}_{i,j+1,k}, \mathrm{DZT}_{i+1,j+1,k}).

Due to these variable cell heights, the spatial discretizations described in Sec. Spatial Discretization need to be modified as shown in the following.

Tracer Advection:

(117)\begin{aligned}
     {\cal L}_T (\varphi) & = \frac{1}{H^T \Delta_y^T}\delta_x
     \left(\overline{H^U \Delta_y^U \,u_x}^y\overline{\varphi}^x\right)
         \nonumber \\
     & + \frac{1}{H^T\Delta_x^T }\delta_y
     \left(\overline{H^U\Delta_x^U \,u_y}^x\overline{\varphi}^y\right )
     +\delta_z\left (w\overline{\varphi}^z\right ),
     \end{aligned}

where \Delta_x^{(U,T)}={\rm DX(U,T)}, \Delta_y^{(U,T)}={\rm DY(U,T)}, H^{(U,T)}={\rm DZ(U,T)}, and the vertical velocity w at bottom T-cells is found as

(118)w_{k+1}-w_{k}=\frac{1}{\Delta_y^T}\,
     \delta_x\left(\overline{\Delta_y^{U}H^U u_x}^y\right )+
     \frac{1}{\Delta_x^T}\,
     \delta_y\left(\overline{\Delta_x^{U}H^U u_y}^x\right )\,.

Laplacian Horizontal Tracer Diffusion:

(119){\cal D}_H(\varphi) = \frac{1}{H^T\Delta_y^T}\delta_x\left(
     \overline{A_H}^x\Delta_y\Delta_z^i\delta_x\varphi\right )
     +\frac{1}{H^T\Delta_x^T}\delta_y\left(
     \overline{A_H}^y\Delta_x\Delta_z^j\delta_y\varphi\right )\,,

where \Delta_x=HTN, \Delta_y=HTE, and

(120)\Delta_z^i=\min({\rm DZT}_{i,j,k},{\rm DZT}_{i+1,j,k}),\quad
     \Delta_z^j=\min({\rm DZT}_{i,j,k},{\rm DZT}_{i,j+1,k}).

Vertical Tracer Diffusion:

(121){\cal D}_V(\varphi) = \frac{1}{H^T}\left [
     \frac{\kappa_{k+\frac{1}{2}}(\varphi_{k+1}-\varphi_k)}
             {({H^T}_{k+1}+{H^T}_{k})/2} -
     \frac{\kappa_{k-\frac{1}{2}}(\varphi_k-\varphi_{k-1})}
             {(H^T_{k}+H^T_{k-1})/2}\right ].

Momentum Advection:

(122)\begin{aligned}
     {\cal L}_U (\alpha) & = \frac{1}{H^U \Delta_y^U}\delta_x
     \left(\overline{\overline{H^U \Delta_y^U \,u_x}^y}^{xy}
     \overline{\alpha}^x\right ) \nonumber \\
     & +\frac{1}{H^U\Delta_x^U }\delta_y
     \left(\overline{\overline{H^U\Delta_x^U \,u_y}^x}^{xy}
     \overline{\alpha}^y\right )
     +\delta_z\left (w\overline{\alpha}^z\right ).
     \end{aligned}

Laplacian Horizontal Friction:

(123)\begin{aligned}
     {\cal D}_H({\bf u}) & = \frac{1}{H^U\Delta_y^U}\delta_x\left(
     \overline{A_M}^x\Delta_y^i\Delta_z^i\delta_x{\bf u}\right )
     +\frac{1}{H^U \Delta_x^U}\delta_y\left(
     \overline{A_M}^y\Delta_x^j\Delta_z^j\delta_y{\bf u}\right )
        \nonumber \\
     & - \frac{A_M}{H^U\Delta_x^U\Delta_y^U}
     \left[\,\overline{\frac{\Delta_y^i (H^U-\Delta_z^i)}{\Delta_x}^x}
     +\overline{\frac{\Delta_x^j (H^U-\Delta_z^i)}{\Delta_y}^y}
     \right ]{\bf u},
     \end{aligned}

where \Delta_x=HTN, \Delta_y=HTE, \Delta_x^j=HUS, \Delta_y^i=HUW, and

(124)\Delta_z^i=\min({\rm DZU}_{i,j,k},{\rm DZU}_{i+1,j,k}),\quad
     \Delta_z^j=\min({\rm DZU}_{i,j,k},{\rm DZU}_{i,j+1,k}).

The last term of (123) is the drag due to the no-slip boundary condition at lateral boundaries. Metric terms in the momentum advection and horizontal friction remain unaffected by partial bottom cells.

Vertical Friction:

(125){\cal D}_V({\bf u}) = \frac{1}{H^U}\left [
     \frac{\kappa_{k+\frac{1}{2}}({\bf u}_{k+1}-{\bf u}_k)}
             {(H^U_{k+1}+H^U_{k})/2}-
     \frac{\kappa_{k-\frac{1}{2}}({\bf u}_k-{\bf u}_{k-1})}
             {(H^U_{k}+H^U_{k-1})/2}\right ].

Hydrostatic Pressure: Since we leave the pressure points at the middle of full cells, the hydrostatic pressure is computed in the same way as for the full-cell case.

Biharmonic Horizontal Diffusion and Friction: The same operations for the Laplacian diffusion and friction given by (119) and (123) for \varphi and {\bf u}, respectively, are applied twice.

Isopycnal Diffusion:

(126)\begin{aligned}
     {\cal D}_H (\varphi) & =
     \frac{1}{H^T\Delta_y^T}\delta_x
     \left[{\Delta_y}\left (\overline{H^T{\overline\kappa}^{z}}^x
     \delta_x\varphi -
      \overline{H^T \overline{\kappa S_x\delta_z\varphi}^{z}}^x
         \right )\right ]  \\
     & - \frac{1}{\Delta_x^T \Delta_y^T}
     \delta_z\left[\frac{1}{\overline {H^T}^z}\left (\overline{H^T
     \overline{{\Delta_x^j}{\Delta_y}
     \kappa  S_x \delta_x\varphi}^x}^z - \overline{H^T\overline{
     {\Delta_x^j}{\Delta_y} \kappa S_x^2}^x}^z
     \delta_z \varphi \right )\right ]
     \nonumber \\
     & + \frac{1}{H^T\Delta_x^T}\delta_y
     \left[{\Delta_x}\left (\overline{H^T{\overline\kappa}^{z}}^y
     \delta_y\varphi -
      \overline{H^T\overline {\kappa S_y
     \delta_z\varphi}^{z}}^y\right )\right ]
     \nonumber \\
     & - \frac{1}{\Delta_x^T \Delta_y^T}
     \delta_z\left[\frac{1}{\overline {H^T}^z}\left (\overline{H^T
     \overline{{\Delta_x}{\Delta_y^i}
     \kappa  S_y \delta_y\varphi}^y}^z - \overline{H^T\overline{
     {\Delta_x}{\Delta_y^i} \kappa S_y^2}^y}^z
     \delta_z \varphi \right )\right ],
     \nonumber\end{aligned}

where \Delta_x=HTN, \Delta_y=HTE, \Delta_x^j=HUS, \Delta_y^i=HUW, S_x=\delta_x\rho/\delta_z\rho, and S_y=\delta_y\rho/\delta_z\rho. In (126), for example, in the (x,z)-plane, \kappa in the two quarter cells in either the east or west direction of T-cells right next to any partial cell is reduced to

(127)\kappa^e_{i,j,k}=\frac{H^T_{i+1,j,k}}{dz_k}\, \kappa_{i,j,k},\quad
     \kappa^w_{i,j,k}=\frac{H^T_{i-1,j,k}}{dz_k}\, \kappa_{i,j,k}\,,

where dz_k is the full-cell height. Similar reductions are made for \kappa in the (y,z)-plane, too. Also \kappa in the lower half of T-cells right above any partial cells and \kappa in the upper half of any partial cells are modified to

(128)\kappa^{bottom}_{i,j,k-1}=\frac{H^T_{i,j,k}}{dz_k}\kappa_{i,j,k-1},\quad
     \kappa^{top}_{i,j,k}=\frac{{H^T}_{i,j,k}}{dz_k}\kappa_{i,j,k}\,.\quad

Remember that \kappa at the lower half of bottom T-cells is set to zero.

Anisotropic Viscosity:

(129)\begin{aligned}
     {\cal D}_H^x({\bf u}) &= \frac{1}{\Delta_y^U}\delta_x
     (\overline{\Delta_y\gamma\sigma_{11}}^x)+\frac{1}{\Delta_x^U}\delta_y
     (\overline{\Delta_x\gamma\sigma_{12}}^y) \nonumber \\
       &-\overline{\Delta_x k_1 \overline{\Delta_y \gamma \sigma_{22}}^x}^x
     +\overline{\Delta_y k_2 \overline{\Delta_x \gamma \sigma_{12}}^y}^y,
     \end{aligned}

(130)\begin{aligned}
     {\cal D}_H^y({\bf u}) &= \frac{1}{\Delta_x^U}\delta_y
     (\overline{\Delta_x\gamma\sigma_{22}}^y)+\frac{1}{\Delta_y^U}\delta_x
     (\overline{\Delta_y\gamma\sigma_{12}}^x) \nonumber \\
       &-\overline{\Delta_y k_2 \overline{\Delta_x \gamma \sigma_{11}}^y}^y
     +\overline{\Delta_x k_1 \overline{\Delta_y \gamma \sigma_{12}}^x}^x,
     \end{aligned}

where \Delta_x=HTN, \Delta_y=HTE, \gamma=1 at four quarter cells at the centered U-cell and \gamma’s are defined in two quarter cells adjacent to the U-cell in the east, west, north and south directions as

(131)\begin{aligned}
     \gamma^e=\frac{\min(H^U_{i,j,k},H^U_{i+1,j,k})}{H^U_{i,j,k}},
       &\gamma^w=\frac{\min(H^U_{i,j,k},H^U_{i-1,j,k})}{H^U_{i,j,k}},
      \\
     \gamma^n=\frac{\min(H^U_{i,j,k},H^U_{i,j+1,k})}{H^U_{i,j,k}},
       &\gamma^s=\frac{\min(H^U_{i,j,k},H^U_{i,j-1,k})}{H^U_{i,j,k}}.
     \end{aligned}

Due to no-slip boundary conditions imposed at lateral boundaries, the velocity gradients to compute the stress tensor \sigma are re-defined as

(132)\begin{aligned}
   (\delta_x u)^e = \frac{1}{\Delta_x}
                    \left(\gamma^e u_{i+1,j,k}- u_{i,j,k}\right),
  &(\delta_x u)^w = \frac{1}{\Delta_x}
                    \left(u_{i,j,k}- \gamma^w u_{i-1,j,k}\right),
     \\
   (\delta_x u)^n = \frac{1}{\Delta_y}
                    \left(\gamma^n u_{i,j+1,k}- u_{i,j,k}\right),
  &(\delta_x u)^s = \frac{1}{\Delta_y}
                    \left(u_{i,j,k}- \gamma^s u_{i,j-1,k}\right),
     \\
   (k_1 {\bf u})^e = k_1^e
           \frac{{\bf u}_{i,j,k}+ \gamma^e {\bf u}_{i+1,j,k}}{2},
  &(k_1 {\bf u})^w = k_1^w
           \frac{{\bf u}_{i,j,k}+ \gamma^w {\bf u}_{i-1,j,k}}{2},
    \\
   (k_2 {\bf u})^n = k_2^n
           \frac{{\bf u}_{i,j,k}+ \gamma^n {\bf u}_{i,j+1,k}}{2},
  &(k_2 {\bf u})^s=k_2^s
           \frac{{\bf u}_{i,j,k}+ \gamma^s {\bf u}_{i,j-1,k}}{2},
    \end{aligned}

where k_1^e, k_1^w, k_2^n and k_2^s are the metric terms defined on the east, west, north and south faces.

7. Advection

POP supports a few advection schemes, primarily for tracers where conservation and monotonicity requirements are stronger. While basic centered schemes for tracers and advection are described in Chapter Spatial Discretization, this chapter describes other advection schemes implemented in POP.

7.1. Third-Order Upwind Advection

Third-order upwind schemes and variants of this scheme were introduced by Leonard (1979a). Our implementation of third-order upwinding simply involves the substitution of three-point interpolation for the two-point centered average in the determination of the tracer value at cell faces, and remains in leapfrog temporal implementation.

It should be noted that improved temporal implementations have been explored in Farrow and Stevens (1995), where they apply a predictor-corrector sequence, and in Holland, Bryan, and Chow (1998), where in addition to documenting our implementation they also describe a reorganization of terms into two groups, one of which is entirely centered in space, the other having a dissipative form.

The finite-difference expression for the advection of tracers is

(133)\begin{aligned}
ADV(i,j,k) = &- ( {u}_E~T^*_E - {u}_W~T^*_W)/\mathrm{DXT} \nonumber \\
             &- ( {v}_N~T^*_N -  {v}_S~T^*_S)/\mathrm{DYT} \nonumber \\
             &- (w_{k}~T^*_T - w_{k+1}~T^*_B)/dz
 \end{aligned}

where :math:` {u}_W, {u}_E, {v}_S, {v}_N` are the T-grid cell face centered velocity components and the directions W, E, S and N are with respect to the logical coordinates.

(134)\begin{aligned}
{u}_E(i) & = (u_{i,j  } \mathrm{DYU}_{i,j  } +
              u_{i,j-1} \mathrm{DYU}_{i,j-1})/
              (2\mathrm{DXT}_{i,j}) \nonumber \\
{u}_W(i) & = {u}_E(i-1) \nonumber \\
{v}_N(j) & = (v_{i  ,j}\mathrm{DXU}_{i  ,j} +
              v_{i-1,j}\mathrm{DXU}_{i-1,j})/
              (2\mathrm{DXT}_{i,j}) \nonumber \\
{v}_S(j) & = (v_{i  ,j-1}\mathrm{DXU}_{i  ,j} +
              v_{i-1,j-1}\mathrm{DXU}_{i-1,j})/
              (2\mathrm{DXT}_{i,j})
 \end{aligned}

and the T^* quantities are the cell face centered tracer concentration at the east, west, north, south, top, and bottom faces. In the standard second-order scheme of POP, these are taken as the arithmetic average of the tracer concentrations in adjacent cells, e.g.,

(135)\begin{aligned}
T^*_E = \frac{1}{2} (T_{i+1, j} + T _{i,j})
 \end{aligned}

The third point included in the three-point interpolation is taken to be in the upwind direction,

(136)T^*_E(i) =
\begin{cases}
 \gamma T_{i-1}~+ \beta ^+ T_i + \alpha ^+ T_{i+1} & {u}_E > 0 \\
 \beta ^- T_i + \alpha ^- T_{i+1} + \delta T_{i+2} & {u}_E < 0
\end{cases}

T^*_W (i) = T^* _E (i-1)

where the coefficients \alpha, \beta, \gamma, \delta are given by

(137)\begin{aligned}
\alpha^+ &= \frac{dx_i(2dx_i + dx_{i-1})}
                 {(dx_i + dx_{i+1})(dx_{i+1} + dx_{i-1} + 2dx _i)}
            \nonumber \\
\alpha^- &= \frac{dx_i(2dx_{i+1} + dx_{i+2})}
                 {(dx_i + dx_{i+1})(dx_{i+1} + dx_{i+2})}
            \nonumber \\
\beta^+ &= \frac{dx_{i+1}(2dx_i + dx_{i-1})}
                {(dx_i + dx_{i+1})(dx_i + dx_{i-1})}
            \nonumber \\
\beta^- &= \frac{dx_{i+1}(2dx_{i+1} + dx_{i+2})}
                {(dx_i + dx_{i+1})(dx_i + dx_{i+2} + 2dx_{i+1})}
            \nonumber \\
\gamma &= \frac{ -dx_i dx_{i+1}}
               {(dx_i + dx _{i-1})(dx_{i+1} + dx_{i-1} + 2dx_i)}
            \nonumber \\
\delta &= \frac{ -dx_i dx _{i+1}}
               {(dx_{i+1} + dx_{i+2})(dx_i + dx_{i+2} + 2dx_{i+1})}
 \end{aligned}

where dx=\mathrm{DXT} and the j-index is suppressed for clarity. The interface tracer values T^*_N, T^*_S, T^*_T, and T^*_B are obtained in an analogous manner.

7.2. Flux-Limited Lax-Wendroff Algorithm

A flux-limited version of the second-order Lax-Wendroff advection scheme has been introduced into POP. The scheme is designed to reduce, but not completely eliminate, the production and amplification of artificial extrema. The artificial extrema are not completely eliminated because the flux limiter is not a true multidimensional limiter. The scheme is constructed by combining a one-dimensional flux-limited scheme with a dimensional splitting technique. We describe in order: time stepping, the non-flux limited version of the one dimensional scheme, the flux limiter, the dimensional splitting technique, and how boundary conditions are handled.

The advection schemes currently available in POP, centered and third-order upwind, use leap-frog time stepping. That is, when computing the tracer tendency between time levels n-1 and n+1 the advective operator {\cal L}_{T} is applied to the tracer \phi at time level n. However, flux-limited advection schemes are generally designed for forward-in-time time stepping. So this scheme applies the advective operator to the tracer \phi time level n-1. Note that velocities are still from time level n and thus are time-centered.

The advection equation in one dimension, in advective form, is

(138)\phi_{t} + u \phi_{x} = 0,

where u is the component of the velocity in the dimension under consideration. We rewrite (138) as

(139)\phi_{t} + \left( u \phi \right)_{x} = u_{x} \phi.

While we assume that \nabla \cdot \vec{U} = 0, we do not assume the velocity to be divergence-free in one dimension. Because of this, (138) and (139) do not necessarily conserve total tracer content. Conservation of total tracer content will later be recovered with the splitting scheme described below.

For our description of the non-flux limited one-dimensional scheme, we temporarily change notation from that used elsewhere in this manual. We consider a one dimensional grid where the jth interval of the grid is denoted I_{j}. Its width is \Delta x_{j}, its midpoint is x_{j}, and its left and right endpoints are x_{j-1/2} and x_{j+1/2} respectively. Thus, we have \Delta x_{j} = x_{j+1/2} -
x_{j-1/2}. Additionally, we define \Delta x_{j+1/2} \equiv \left(
\Delta x_{j} + \Delta x_{j+1} \right) / 2 = x_{j+1} - x_{j}. The advection scheme under consideration discretizes (139) as

(140)\phi^{n+1}_{j} = \phi^{n-1}_{j}
   - \frac{2 \Delta t}{\Delta x_{j}}
      \left( F_{j+1/2} - F_{j-1/2} \right)
   + \frac{2 \Delta t}{\Delta x_{j}}
      \left( u_{j+1/2} - u_{j-1/2} \right) \phi^{n-1}_{j},

where super- and sub-scripts indicate time and space indices respectively and \Delta t is the time step. F_{j+1/2}, the advective flux of \phi through x_{j+1/2}, is given by

F_{j+1/2} \equiv u_{j+1/2} \phi_{j+1/2},

where \phi_{j+1/2} is an approximation to the average value of \phi at x_{j+1/2} over t^{n-1} \leq t \leq t^{n+1}. Since \phi_{j+1/2} does not represent \phi at a particular time level, we do not denote a time level in a superscript. \phi_{j+1/2} is constructed by

  1. approximating \phi\left(x,t^{n-1}\right) with a linear fit p(x),
  2. advecting p(x) exactly with velocity u_{j+1/2},
  3. taking \phi_{j+1/2} to be the average value of the advected p(x) at x_{j+1/2} over t^{n-1} \leq t \leq t^{n+1}.

The linear fit p(x) is constructed by equating p(x)’s average values over intervals I_{j+k} near x_{j+1/2} to \phi^{n-1}_{j+k} for k \in {0,1}. That is, p(x) satisfies the constraints

(141)\frac{1}{\Delta x_{j+k}} \int\limits_{I_{j+k}} p(x) \,\mathrm{d} x
   = \phi^{n-1}_{j+k}.

When the constructed p(x) is advected and averaged at x=x_{j+1/2} over t^{n-1} \leq t \leq t^{n+1}, the result is

(142)\begin{aligned}
\phi_{j+1/2} & = \phi^{n-1}_{j} + \frac{1}{2}
   \left( \Delta x_{j} - 2 u_{j+1/2} \Delta t \right)
   \frac{\phi^{n-1}_{j+1} - \phi^{n-1}_{j}}{\Delta x_{j+1/2}}
 \\
& = \phi^{n-1}_{j+1} - \frac{1}{2}
   \left( \Delta x_{j+1} + 2 u_{j+1/2} \Delta t \right)
   \frac{\phi^{n-1}_{j+1} - \phi^{n-1}_{j}}{\Delta x_{j+1/2}}.
\end{aligned}

This is the familiar Lax-Wendroff scheme in a form that is appropriate for a stretched grid. We present the different forms of the same expression in anticipation of the application of flux limiters.

The above description assumes that the velocity is specified at interval boundaries. This leads to a natural implementation in the multidimensional setting when normal velocities are specified at midpoints of tracer cell faces, as is done with Arakawa C-grids (Arakawa and Lamb 1977). However, POP utilizes an Arakawa B-grid, where velocities are specified at tracer cell corners. Since the one-dimensional scheme use the face velocity for the tracer face value construction, e.g. (142), it is necessary to compute a face velocity from the B-grid corner velocities. The normal velocity for the east face of the i,j tracer cell is computed as

(143)\mathrm{U}^{\mathrm{east}}_{i,j} = 0.5
   \left( \mathrm{DYU}_{i,j} \mathrm{U}_{i,j}
          + \mathrm{DYU}_{i,j-1} \mathrm{U}_{i,j-1} \right)
   / \mathrm{HTE}_{i,j},

where \mathrm{DYU}_{i,j} = \left(\mathrm{HTE}_{i,j} + \mathrm{HTE}_{i,j+1}\right)/2, \mathrm{HTE}_{i,j} is the length of the east face, and \mathrm{U}_{i,j} and \mathrm{U}_{i,j-1} are the zonal velocities at the north-east and south-east corners respectively. We note that (143) requires adjustment if POP is configured with partial bottom cells and that a formula analogous to (143) holds for \mathrm{V}^{\mathrm{north}}_{i,j}.

In one dimension, we employ a flux limiter that is based on the work of Hundsdorfer and Trompert (1994), noting that the flux limiter presented there is in fact equivalent to the ULTIMATE flux limiter of Leonard (1979b). Solutions of (138) are constant along characteristics, which implies that local extrema are not created or amplified. The purpose of the flux limiter is to constrain the discretization to have this property as well. This is done by modifying, as necessary, the advective flux F_{j+1/2}. In order to simplify the presentation, we initially assume u>0. We write F_{j+1/2} as

(144)\begin{aligned}
F_{j+1/2} &= u_{j+1/2} \phi_{j+1/2} \nonumber \\
          &= u_{j+1/2} \left( \phi^{n-1}_{j} + \psi_{j+1/2}
      \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) \right),
\end{aligned}

where \psi_{j+1/2} is to be determined. A priori, we will impose 0 \leq \psi_{j+1/2} \leq 1 in order to avoid extrapolation in the construction of \phi_{j+1/2}. We seek further constraints on \psi_{j+1/2} that will ensure

(145)\min \left( \phi^{n-1}_{j-1}, \phi^{n-1}_{j} \right)
\leq \phi^{n+1}_{j} \leq
\max \left( \phi^{n-1}_{j-1}, \phi^{n-1}_{j} \right).

Note that constraints on \psi_{j+1/2} are ensuring extrema constraints on the interval that is upwind of x_{j+1/2} and that these extrema constraints are expressed in terms of upwind values. Given this form for F_{j+1/2}, (140) becomes

(146)\begin{split}
\phi^{n+1}_{j} & = \phi^{n-1}_{j} - \frac{2 \Delta t}{\Delta x_{j}} \Big[
   \mathbin{\phantom{-}} u_{j-1/2} \phi^{n-1}_{j\phantom{-1}}
      + u_{j+1/2} \psi_{j+1/2}
        \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) \\
& \mathrel{\phantom{=}}
  \phantom{\phi^{n-1}_{j} - \frac{2 \Delta t}{\Delta x_{j}} \Big[}
   - u_{j-1/2} \phi^{n-1}_{j-1} - u_{j-1/2} \psi_{j-1/2}
      \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right) \Big]
\end{split}

We now consider two very specific cases and then the more general third case:

: \phi^{n-1}_{j} = \phi^{n-1}_{j-1}

In this case, (146) reduces to

\phi^{n+1}_{j} = \phi^{n-1}_{j} -
                 \frac{2u_{j+1/2}\Delta t}{\Delta x_{j}}
   \psi_{j+1/2} \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right), \nonumber

and (145) reduces to \phi^{n+1}_{j} = \phi^{n-1}_{j}. Together, these imply \psi_{j+1/2}=0.

: \phi^{n-1}_{j} = \phi^{n-1}_{j+1}

In this case, (146) reduces to

\phi^{n+1}_{j} = \phi^{n-1}_{j} -
                 \frac{2u_{j-1/2}\Delta t}{\Delta x_{j}}
   \left( 1 - \psi_{j-1/2} \right)
   \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right), \nonumber

so there is no constraint on \psi_{j+1/2}.

: \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right) \neq 0

Letting \theta_{j} = \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right) \big/ \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right), (146) becomes

\begin{split}
\phi^{n+1}_{j} & = \bigg[ 1 - \frac{2 \Delta t}{\Delta x_{j}}
   \left( u_{j-1/2} \left( 1 - \psi_{j-1/2} \right)
          + u_{j+1/2} \psi_{j+1/2} / \theta_{j} \right)
   \bigg] \phi^{n-1}_{j} \\
& \mathrel{\phantom{=}}\phantom{\bigg[ 1}
   + \frac{2 \Delta t}{\Delta x_{j}}
   \left( u_{j-1/2} \left( 1 - \psi_{j-1/2} \right)
          + u_{j+1/2} \psi_{j+1/2} / \theta_{j} \right)
   \phi^{n-1}_{j-1}. \nonumber
\end{split}

This expresses \phi^{n+1}_{j} as a linear combination of \phi^{n-1}_{j} and \phi^{n-1}_{j-1} where the weights sum to 1. Given such a representation, (145) is equivalent to requiring the weights to be between 0 and 1. Thus, (145) is equivalent to

0 \leq \frac{2 \Delta t}{\Delta x_{j}}
   \left( u_{j-1/2} \left( 1 - \psi_{j-1/2} \right)
          + u_{j+1/2} \psi_{j+1/2} / \theta_{j} \right)
  \leq 1,

which in turn is equivalent to

(147)0 \leq \frac{u_{j-1/2}}{u_{j+1/2}}
       \left( 1 - \psi_{j-1/2} \right) + \psi_{j+1/2} / \theta_{j}
  \leq \frac{\Delta x_{j}}{2 u_{j+1/2} \Delta t},

where we have used the assumption u_{j+1/2}>0. We now suppose that

(148)\psi_{j+1/2} = \max \left( 0,
   \min \left( 1, \psi^{0}_{j+1/2},
               \mu_{j} \theta_{j} \right) \right),

where \mu_{j} \ge 0 is to be determined, and \psi^{0}_{j+1/2} is the value of \psi_{j+1/2} that corresponds to the non-flux limited advection scheme. Then

0 \leq \psi_{j+1/2} \leq 1 \mbox{ and }
0 \leq \psi_{j+1/2} / \theta_{j} \leq \mu_j.

It follows that

0 \leq \frac{u_{j-1/2}}{u_{j+1/2}} \left( 1 - \psi_{j-1/2} \right)
       + \psi_{j+1/2} / \theta_{j}
  \leq \frac{u_{j-1/2}}{u_{j+1/2}} + \mu_{j}.

Comparing this to (147), we see that (147), and thus (145), will be ensured if

\mu_{j} \leq \frac{\Delta x_{j}}{2 u_{j+1/2} \Delta t}
   - \frac{u_{j-1/2}}{u_{j+1/2}}
   = \frac{\Delta x_{j} - 2 u_{j-1/2} \Delta t}{2 u_{j+1/2} \Delta t}.

In the above presentation we assumed that u>0. The velocities used in the analysis were u_{j+1/2} and the velocity in the upwind direction, u_{j-1/2}. With this in mind, the remaining cases to consider are u<0, u_{j+1/2}>0 and u_{j-1/2}<0, and u_{j+1/2}<0 and u_{j+3/2}>0. The last velocity sign pair is merely a shift by +1 of the j indices from the second pair, so it is not necessary to separately consider this case. We note that the converging velocity pair u_{j+1/2}<0 and u_{j-1/2}>0 does not arise in our analysis because when u_{j+1/2}<0, the analysis depends on u_{j+3/2}.

We now consider the case where u<0. In order to maintain the notion that smaller values of \psi_{j+1/2} correspond to greater diffusion, (144) is replaced with

(149)F_{j+1/2} = u_{j+1/2} \phi_{j+1/2}
   = u_{j+1/2} \left( \phi^{n-1}_{j+1} - \psi_{j+1/2}
      \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) \right).

In order to maintain the notion that constraints on \psi_{j+1/2} are ensuring extrema constraints on the upwind interval, (145) is replaced with

(150)\min \left( \phi^{n-1}_{j+2}, \phi^{n-1}_{j+1} \right)
   \leq \phi^{n+1}_{j+1} \leq
\max \left( \phi^{n-1}_{j+2}, \phi^{n-1}_{j+1} \right).

Assuming \psi_{j+1/2} = \max \left( 0,  \min \left( 1, \psi^{0}_{j+1/2},  \mu_{j+1} / \theta_{j+1} \right) \right), we find, with similar algebra, that (150) will be satisfied if

\mu_{j+1} \leq
   \frac{\Delta x_{j+1} + 2 u_{j+3/2} \Delta t}{-2 u_{j+1/2} \Delta t}.

We finally consider the case where u_{j+1/2}>0 and u_{j-1/2}<0. In this case, we use (144) for F_{j+1/2} and (149) for F_{j-1/2}. Then (140) becomes

\begin{split}
\phi^{n+1}_{j} & = \phi^{n-1}_{j} - \frac{2 \Delta t}{\Delta x_{j}} \Big[
   \mathbin{\phantom{+}}
   u_{j+1/2} \psi_{j+1/2} \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) \\
& \mathrel{\phantom{=}}
  \phantom{\phi^{n-1}_{j} - \frac{2 \Delta t}{\Delta x_{j}} \Big[}
   + u_{j-1/2} \psi_{j-1/2} \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right)
   \Big].
\end{split}

Since x_{j} is upwind of both x_{j+1/2} and x_{j-1/2}, we use constraints on both \psi_{j+1/2} and \psi_{j-1/2} to ensure extrema constraints on \phi^{n+1}_{j}. Since there are no points upwind of x_{j}, the extrema constraints reduce to \phi^{n+1}_{j}=\phi^{n-1}_{j}. This can only be satisfied if

(151)0 = u_{j+1/2} \psi_{j+1/2} \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right)
  + u_{j-1/2} \psi_{j-1/2} \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right)

We now consider two very specific cases and then the more general third case:

: \phi^{n-1}_{j} = \phi^{n-1}_{j-1}

Then we must have \psi_{j+1/2}=0. The value of \psi_{j-1} is irrelevant.

: \phi^{n-1}_{j+1} = \phi^{n-1}_{j}

Then we must have \psi_{j-1/2}=0. The value of \psi_{j+1} is irrelevant.

: \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) \left( \phi^{n-1}_{j} - \phi^{n-1}_{j-1} \right) \neq 0

Then (151) can be divided by \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right) to obtain

(152)0 = u_{j+1/2} \psi_{j+1/2} + u_{j-1/2} \psi_{j-1/2} \theta_{j}.

If \theta_{j}<0, the only solution that satisfies 0 \leq \psi \leq1 is \psi_{j+1/2} = \psi_{j-1/2} = 0. If \theta_{j}>0, we choose the largest values of \psi_{j-1/2} and \psi_{j+1/2} that satisfy (152), \psi_{j-1/2} \le \max \left( 0, \min \left(1, \psi^{0}_{j-1/2} \right) \right), and \psi_{j+1/2} \le \max \left( 0, \min \left(1, \psi^{0}_{j+1/2} \right) \right). We do this by setting

\begin{aligned}
\psi^{*}_{j+1/2} & = \max
   \left( 0, \min \left( 1, \psi^{0}_{j+1/2} \right) \right), \\
\psi_{j-1/2} & = \max
   \left( 0, \min \left( 1, \psi^{0}_{j-1/2},
          -u_{j+1/2} / u_{j-1/2}
          \cdot \psi^{*}_{j+1/2} / \theta_{j} \right) \right), \\
\psi_{j+1/2} & = -u_{j-1/2} / u_{j+1/2} \cdot \psi_{j-1/2} \theta_{j}.\end{aligned}

The dimensional splitting technique follows the technique described in section 2.17.4 of Adcroft et al. (2005). In our implementation of this technique, (140) and (144) are applied sequentially in the z, x, and y directions. Each application of this pair of equations uses the results from the previous application, except for the last term of (140), where we use \phi^{n-1} in each iteration. Using \phi^{n-1} for this term in each iteration leads to conservation of total tracer content, assuming as we do that the velocity field is divergence free. Our choice of the ordering of dimensions is based on simplifying the implementation. We tested alternating the dimension order in simulations of the Stommel gyre test problem of Hecht, Holland, and Rasch (1995) and did not find a meaningful benefit, so we do not vary the dimension order in the three-dimensional setting.

The flux limiter extends the discrete domain of dependence of Lax-Wendroff in the upstream direction, potentially beyond lateral boundaries. For example, consider the situation where we are computing the ’zonal’ flux at x_{j+1/2}, u_{j+1/2}>0, and there is a solid boundary at x_{j-1/2}. Using (144) for F_{j+1/2} and u_{j-1/2}=F_{j-1/2}=0 in (140) yields

\phi^{n+1}_{j} = \phi^{n-1}_{j}
   - \frac{2 \Delta t}{\Delta x_{j}} \psi_{j+1/2} u_{j+1/2}
      \left( \phi^{n-1}_{j+1} - \phi^{n-1}_{j} \right).

Since x_{j} is upwind of x_{j+1/2}, we use constraints on \psi_{j+1/2} to ensure extrema constraints on \phi^{n+1}_{j}. Since there are no points upwind of x_{j} in the direction under consideration, the extrema constraints reduce to \phi^{n+1}_{j}=\phi^{n-1}_{j}. This can only be satisfied if \psi_{j+1/2}=0. Examining (148), this can be achieved by using \theta_{j}=0, which is equivalent to using a fictitious \phi^{n-1}_{j-1}=\phi^{n-1}_{j}. In other words, we use a zero Neumann boundary condition at no-flow boundaries.

8. Sub-grid Scale Mixing Parameterizations

8.1. Horizontal Tracer Diffusion

8.1.1. Laplacian Horizontal Diffusivity

The spatial discretization of the standard Laplacian horizontal tracer diffusion with a spatially varying tracer diffusivity A_{_H} is described in Sec. Horizontal Tracer Diffusion.

8.1.2. Biharmonic Horizontal Diffusivity

Biharmonic horizontal tracer diffusion is implemented by applying the Laplacian operator (41) twice. Specifically,

(153){\cal D}_H^{(4)}(\varphi) = {\cal D}_H^{(2)}
(A_H^{(4)}{\cal D}_H^{(2)}(\varphi))

where the superscripts (4) and (2) refer to biharmonic and harmonic (Laplacian) operators, and {\cal D}_H^{(2)} is given by (41) with A_H=1. A_H^{(4)} is the biharmonic diffusivity, which may be spatially varying, and should be negative for positive-definite diffusion. Here the biharmonic diffusivity is sandwiched between the two applications of the Laplacian operator. An alternate approach (not implemented in the code), following Griffies and Hallberg (2000), is to include the square-root of -A_H inside the second derivative of each {\cal D}_H^{(2)} operator.

8.1.3. The Gent-McWilliams Parameterization

The transport equation of tracer \varphi is given by

(154)\frac{\partial}{\partial t}\varphi +
  ({\bf u} + {\bf u}^*)\cdot \nabla\varphi +
  (w+w^*)\frac{\partial}{\partial z}\varphi =
  R(\varphi) + {\cal D}_V(\varphi),

where the bolus velocity induced by mesoscale eddies is parameterized, from Gent and McWilliams (1990), as

(155){\bf u}^* =\left(\nu \frac{\nabla \rho}{\rho_z}\right)_z, \qquad
     w^* = -\nabla\cdot \left(
     \nu \frac{\nabla \rho}{\rho_z}\right ),

where \nu is a thickness diffusivity and subscripts x,y,z on \rho and tracers \varphi denote partial derivatives with respect to those variables (this convention will be followed below). The Redi isoneutral diffusion operator (Redi 1982) for small slope can be written as

(156)R(\varphi) = \nabla_3 \cdot (\mathrm{\bf K} \cdot \nabla_3 \varphi),

where the subscript 3 indicates the three-dimensional gradient or divergence, i.e., \nabla_3 = (\nabla, \frac{\partial}{\partial z}). The symmetric tensor K is defined as

(157){\rm \bf K} = \kappa_I \left (
  \begin{array}{ccc}
  1 & 0 & -{\rho_x/\rho_z} \\
  0 & 1 & -{\rho_y/\rho_z} \\
  -{\rho_x/\rho_z} & -{\rho_y/\rho_z} &
       {\vert \nabla \rho\vert^2/\rho_z^2}
  \end{array}\right ),

This tensor describes along-isopycnal diffusion that is isotropic in the two horizontal dimensions. The general anisotropic form of (157) is given in Smith (1999). The isopycnal diffusivity \kappa_I is in general a function of space and time, and a parameterization for variable \kappa_I will be described at the end of this section. In POP, we write the bolus velocity in the skew-flux form (Griffies 1998):

(158){\bf u}_3^*\cdot \nabla_3\varphi
  =   \nabla_3\cdot ({\bf u}_3^*\varphi)
  = - \nabla_3\cdot ({\rm \bf B}\cdot \nabla_3\varphi),

where we have used :math:` nabla_3cdot {bf u}_3^*=0`. The subscript 3 on the velocity indicates the three-dimensional velocity, i.e., {\bf u}_3^* = ({\bf u}^*, w^*). The antisymmetric tensor B is given by

(159){\rm \bf B}=\nu \left (
  \begin{array}{ccc}
  0 & 0 & {\rho_x/\rho_z} \\
  0 & 0 & {\rho_y/\rho_z} \\
  -{\rho_x/\rho_z} & -{\rho_y/\rho_z} & 0
  \end{array}\right ).

By substituting (159), the transport equation (154) can be written

(160)\frac{\partial}{\partial t}\varphi + {\bf u}\cdot\nabla\varphi +
  w\frac{\partial}{\partial z}\varphi =
  \overline R(\varphi) + {\cal D}_V(\varphi), \qquad
  \overline R(\varphi) = \nabla_3\cdot(\mathrm{\bf K+B})\cdot
  \nabla_3\varphi.

With a nonlinear equation of state, the isopycnal slopes (L_x,L_y) = (-\rho_x /\rho_z, -\rho_y /\rho_z) that appear in (157) and (159) should be replaced with slopes along the local neutral surfaces (McDougall 1987), which are given by:

(161)\begin{aligned}
L_x =  -\frac{\alpha_p\Theta_x - \beta_p S_x}
             {\alpha_p\Theta_z - \beta_p S_z}
\end{aligned}

with a similar definition for L_y. Here \alpha_p = -\partial\, \rho_p / \partial \Theta and \beta_p = \partial\, \rho_p / \partial S, where \Theta and S are the model potential temperature and salinity, and \rho_p is the potential density referenced to the local pressure (or depth).

8.1.3.1. Discretization of the Diffusion Operator

The functional formalism for deriving positive-definite diffusive operators in ocean models was originally used to implement the Gent-McWilliams parameterization in POP. Later, a collaboration between GFDL and LANL led to an improved implementation that is described in Griffies et al. (1998). The main change in the new formulation is to discretize the neutral slopes (161) in a particular way that avoids exciting a nonlinear instability.

The discrete functional formalism is based on the property of the continuous dissipation operator (160) that it can be expressed as the functional derivative of the integral of a positive-definite functional {\cal G}:

(162)\begin{aligned}
{\cal G}[\varphi] &= \frac{1}{2}\int dV
   (\nabla_3\varphi)\cdot{\bf K}\cdot(\nabla_3\varphi) \nonumber \\
   &= \frac{1}{2} \int dV \kappa_I
       \left[(\varphi_x + L_x\varphi_z)^2
           + (\varphi_y + L_y\varphi_z)^2\right]
    \\
R(\varphi) &= -\frac{\delta}{\delta \varphi} {\cal G}[\varphi]
\end{aligned}

where \delta denotes a functional derivative. That is, {\cal G} is a functional whose derivative with respect to a tracer field at a given point in space yields the isoneutral diffusive operator R at that point. Note that, as a consequence of the assumption of isotropic diffusion (157), the 3-D functional splits into two terms which are effectively 2-D, one in the x-z plane and one in the y-z plane.

The procedure for deriving the discrete operator is to first discretize the functional and then take its derivatives with respect to the tracer field at a given point on the computational grid; this yields the positive-definite diffusion operator at that point. In order to discretize the functional on the POP grid, we need to define all the quantities appearing in {\cal G} within each cell. Because gradients involve differences across cell faces, it is convenient to subdivide each T-cell into four subcells, as shown in Figure7.1 , which shows a central T-cell, its four interior subcells, and the eight surrounding subcells that contribute to the friction operator when the discrete functional derivative is taken with respect to the tracer value at the central point. A unique value of the tracer gradients and isoneutral slopes can be assigned to each subcell, and from these the complete discrete functional {\cal G} can be constructed and summed over the grid. Associated with each subcell is the set of its three nearest points (known as the “triad” of points associated with that subcell), namely, the central point of the full cell to which it belongs, and the two nearest points in the horizontal and vertical directions. For example, the triad associated with the subcell in Figure7.1  labeled (o,ne) is the set of points labeled (o), (n), (e). All the derivatives of the tracer fields needed to construct the functional in a given subcell only involve the tracer values at its triad points. For example, the x and z derivatives of \varphi in the subcell labeled (o,ne), are given by \varphi_x^{(o,ne)}=(\varphi^{(e)}-\varphi^{(o)})/h_x and \varphi_z^{(o,ne)}=(\varphi^{(n)}-\varphi^{(o)})/h_z, where h_x and h_z are defined in each subcell as the distances between the central point (o) and the other two triad points in the x and z directions, respectively. Similarly, the derivatives in the subcell labeled (o,nw) are given by \varphi_x^{(o,nw)}=(\varphi^{(o)}-\varphi^{(w)})/h_x and \varphi_z^{(o,nw)}=(\varphi^{(n)}-\varphi^{(o)})/h_z. The diffusion coefficient is assumed to have a constant value across all four subcells within a given full cell, but can vary from full-cell to full cell.

Using these definitions of the quantities in each subcell appearing in (162), it can be written in the discrete form:

(163){\cal G}[\varphi] = \frac{1}{2}\sum_{ij}\sum_{n=1}^4
     (\kappa_I)_{ij} V_{ijn} [(\varphi_x + L_x\varphi_z)^2
                            + (\varphi_y + L_y\varphi_z)^2]

where subscripts ij run over all full cells, and subscripts n over subcells within a given full cell. V_{ijn} = \frac{1}{4}h_x h_y dz(k) is the subcell volume (where h_y is the full cell width centered at the same point as h_x, and note that the height of each subcell is taken to be half its full-cell height dz(k) .)

The diffusion operator is then given by the derivative of the discrete functional with respect to the velocity at a given point (i',j'):

(164)R_{i'j'} = -\frac{1}{V^T_{i'j'}}
            \frac{\partial}{\partial\varphi_{i'j'}}
           G[\varphi]

where V^T_{i'j'}= A^T_{i'j'}\,dz(k) full-cell volume of cell (i,j,k), and A^T_{i'j'} is the full-cell area. (Note: in the code A^T=\mathrm{TAREA}, which is not exactly equal to the sum of the four subcell areas given above, however, the difference is small and is absorbed in an overall normalization factor.) When taking the derivative with respect to \varphi at the central point, only those subcells in the full functional whose triad points contain the central point will contribute to the derivative, and these are just the 12 subcells shown in Figure7.1 . It is straightforward to take the derivative in (164) in these 12 subcells and assemble the results into a compact form for the diffusion operator. The result is:

../_images/FigGMSubCells.jpg

Figure 7.1: The discretization of the x-z part of the isoneutral diffusion functional is accomplished by subdividing each T-cell (shown bounded by thick lines) into four subcells. The central tracer point (T-point) is denoted by “(o)” and the surrounding T-points by “(n)”, “(s)”, “(ne)”, etc. (this is for notation only, points “(n)” and “(s)” are displaced in the vertical direction.) The four subcells surrounding the central point are labeled “(o,ne)”, “(o,nw)”, etc., and the subcells surrounding the T-cell to the east are labeled “(e,nw)”, etc. Within each subcell a unique value is defined for all quantities (e.g., the gradients and slope factors appearing in (162) that are needed to evaluate the functional. Only the twelve subcells shown contribute to the friction operator at the central point (see text). (Note: this figure shows a case where the fundamental grid is constructed with U-points in the exact middle of T-cells in a grid with variable horizontal spacing. More commonly T-points are located at the exact center of U-cells. This choice is made during grid generation and is only determined in the code through the input grid file.) This figure can also be used for constructing the anisotropic dissipation functional (see Sec. Anisotropic Horizontal Viscosity) where the central point represents a U-point.

(165)\begin{aligned}
V^{T}_{ij}R_{ij}(\varphi) &=
    \overline{(V/h_x)\kappa_I(\varphi_x + L_x\varphi_z)}^e -
    \overline{(V/h_x)\kappa_I(\varphi_x + L_x\varphi_z)}^w
    \nonumber \\
  &+\overline{(V/h_z)\kappa_I(L_x\varphi_x + L_x^2\varphi_z)}^n -
    \overline{(V/h_z)\kappa_I(L_x\varphi_x + L_x^2\varphi_z)}^s
    \nonumber \\
  &+ \ \mathrm{similar\ terms\ with}\ x\ \mathrm{replaced\ by}
       \ y\ \mathrm{everywhere.}
\end{aligned}

Here overbars with superscripts e, w, etc., denote sums over the four subcells surrounding the east, west, etc., faces of the central full cell. In computing the slopes L_x and L_y given by (161) it is important that the coefficients \alpha_p and \beta_p be evaluated at the pressure corresponding to the level of the central point of the triad in a given subcell (Griffies et al. 1998); this avoids exciting the nonlinear instability mentioned above. No-flux boundary conditions on the tracers are implemented by setting \kappa_I = 0 in all subcells adjacent to lateral and bottom boundaries. There is no way to derive the bolus transport term in (160) and (159) from a functional like the diffusive term, but the bolus term in (160) involves terms proportional to the off-diagonal terms in the diffusion tensor, and these are discretized in the same way. Thus we find for the full skew flux (\overline R(\varphi) in (160))

(166)\begin{aligned}
V^{T}_{ij}\overline{R}_{ij}(\varphi) &=
    \overline{(V/h_x)[\kappa_I\varphi_x + (\kappa_I-\nu)L_x\varphi_z]}^e
    \nonumber \\
 &- \overline{(V/h_x)[\kappa_I\varphi_x + (\kappa_I-\nu)L_x\varphi_z]}^w
    \nonumber \\
 &+ \overline{(V/h_z)[(\kappa_I+\nu)L_x\varphi_x +
                       \kappa_I   L_x^2\varphi_z]}^n
    \nonumber \\
 &- \overline{(V/h_z)[(\kappa_I+\nu)L_x\varphi_x +
                       \kappa_I L_x^2\varphi_z)}^s
    \nonumber \\
 &+ \ \mathrm{similar\ terms\ with}\ x\ \mathrm{replaced\ by}
                                   \ y\ \mathrm{everywhere.}
\end{aligned}

If the tracer and thickness diffusion coefficients are chosen equal (\kappa_I = \nu), then the second terms in brackets in the first line of (166) vanish and need not be computed. (Note: in the code, the above formulae are evaluated using the slope functions \mathrm{SLX} = -(h_x/h_z)L_x, \mathrm{SLY} = -(h_y/h_z)L_y, and the differences \Delta_x\varphi = h_x\varphi_x, etc.) Finally, it should be noted that the last terms in (166) involving L_x^2 (and L_y^2) are separated and treated implicitly, because the effective viscous coefficient can be large enough to violate explicit CFL conditions (Cox 1987). These terms can be written

(167)\begin{aligned}
  \overline{(V/h_z)\kappa_I L_x^2\varphi_z}^n
- \overline{(V/h_z)\kappa_I L_x^2\varphi_z}^s
= \delta_z(\kappa' \delta_z\varphi),
\nonumber \\
 \mathrm{where\ } \kappa' = (A^T)^{-1}
  \overline{(V/h_z)\kappa_I (L_x^2+ L_y^2)}^n
\end{aligned}

and the last line gives \kappa' on the top of the T-cell. \kappa' is then added to the vertical diffusivity \kappa in (42) and the vertical diffusion is solved implicitly. For this reason implicit vertical mixing must be used with the Gent-McWilliams parameterization. One of the design features of the diffusive operator is that with a linear equation of state it should be zero when applied to \varphi = \rho. This condition is violated when the above term is treated implicitly while the remaining terms are treated explicitly, but in practice this leads to a very small error and is ignored.

Two tapering functions are often applied to the GM coefficient \nu and the Redi mixing coefficient \kappa_I. The first is for physical reasons and reduces the GM coefficient in the upper boundary layer, see Large et al. (1997). Here the ocean is not nearly adiabatic, and so the GM parameterization does not apply; instead the isopycnals tend to become more vertical because of the strong vertical mixing in the upper boundary layer. The tapering occurs where the depth is less than the local isopycnal slope multiplied by the local Rossby radius L_R defined by L_R=\vert f\vert/c_0 with c_0=2m/s and L_R is bounded by 15\ {\rm km}\le L_R\le 100\ {\rm km}. Two forms for the factor are in the code, with the second being faster because there is no sine function.

tanh option:

(168)f_1 = \frac{1}{2} \Big[1+
 \sin\Big(\pi (\min(1, -z/(\vert{\bf L}\vert L_R))-0.5)\Big)\Big]

else:

(169)\begin{aligned}
f_1 = \frac{1}{2} +
  2&\Big(min(1, -z/(\vert{\bf L}\vert L_R))-0.5\Big)\times \nonumber \\
   &\Big(1 - \vert min(1, -z/(\vert{\bf L}\vert L_R))-0.5 \vert \Big)
 \end{aligned}

where |{\bf L}| = (L_x^2+L_y^2)^{\frac{1}{2}}. The second function ensures numerical stability when the isopycnals become too steep. There are four options to reduce the coefficients \kappa_I and \nu for \vert{\bf L}\vert \le S_M^{(\kappa_I,\nu)}, where the maximum slope allowed can be different for S_M^{\kappa_I} and S_M^{\nu}. The default values are S_M^{\kappa_I} =S_M^{\nu} = 0.01. The first option was described in Danabasoglu and McWilliams (1995) as

(170)f_2 = \frac{1}{2}
  \Big [1-\tanh (10\vert{\bf L}\vert/ S_M^{(\kappa_I,\nu)}-4)\Big ]

for \vert{\bf L}\vert > S_M^{(\kappa_I,\nu)}, \kappa_I=\nu=0. The second option is a faster version without a tanh function:

(171)f_2 = \frac{1}{2}
\Big [1-\Big(2.5\vert{\bf L}\vert/ S_M^{(\kappa_I,\nu)}-1\Big)
\Big( 4 - \vert 10\vert{\bf L}\vert/ S_M^{(\kappa_I,\nu)}-4\vert\Big)\Big]

with f_2 = 1 for \vert{\bf L}\vert \le 0.2 S_M^{(\kappa_I,\nu)} and f_2 = 0 for \vert{\bf L}\vert \ge 0.6 S_M^{(\kappa_I,\nu)}.

The third option is the clipping option described in Cox (1987) where the slopes are set to the maximum if they are diagnosed as larger. In this case, mixing occurs along the maximum allowed slope. This implies some cross isopycnal mixing when the slopes are steep. The default values are S_M^{\kappa_I} =S_M^{\nu} = 0.01.

The fourth option is described in Gerdes, Koberle, and Willebrand (1991) and the factor is

(172)f_2 = (S_M^{(\kappa_I,\nu)}/\vert{\bf L}\vert)^2

when the slope is larger than the maximum slope. This option and the first two options retain the adiabatic nature of the parameterization. When \kappa_I=\nu and S_M^{\kappa_I} =S_M^{\nu}, the code runs much faster because the terms proportional to \kappa_I-\nu disappear. Because CCSM uses the Near Surface Eddy Flux Paramaterization described below, no f_1 function is needed. For f_2, CCSM uses the second option, but the maximum slopes allowed are much larger than the default values with S_M^{\kappa_I} =S_M^{\nu} = 0.3.

8.1.3.2. Variable GM Diffusivity

Two options for specifying a variable diffusivity are available. The first follows Visbeck et al. (1997), where it is assumed that

(173)\kappa_I = \alpha\frac{L_E^2}{T_E}

with a typical value of the constant \alpha as 0.13. In (173), the baroclinic length scale L_E and the Eady time scale T_E are defined as

L_E = \min\left(c/\vert f\vert,\sqrt{c/2\beta}\right ),

\frac{1}{T_E}=\frac{\max(\vert f\vert, \sqrt{2c\beta})}{\sqrt{Ri}},

Ri=-\frac{g}{\rho_{_0} H}\int_{z_2}^{z_1}
        \frac{\rho_z}{\vert {\bf u}_z\vert^2}{\rm d}z,

where \beta is the meridional gradient of the Coriolis parameter and c is the first baroclinic wave speed given by c^2=-gH\int_{z_2}^{z_1}\rho_z\,{\rm d} z/\rho_{_0}, z_1=-50{\rm m} and z_2=-1000{\rm m}. L_E is also bounded below by the smaller horizontal side of each grid box. In the code, \kappa_I is set to be bounded by 3\times10^6\le\kappa_I\le 2\times 10^7 cm^2/s and \kappa_I is set to the lower bounding value where the model depth is less than -z_1.

A second option is meant to parameterize submesoscale eddies. This option varies the thickness, \nu, and isopycnal diffusivity, \kappa_I, using a simplified version of the closure of Eden and Greatbatch (2008) (note that thickness diffusivity and isopycnal diffusivity are here assumed to be identical). The thickness diffusivity is given by \kappa_I=c L^2 \sigma, where \sigma denotes an inverse eddy time scale and c a constant parameter of order one. L denotes an eddy length scale, which is taken as the minimum of the local Rossby radius L_r and the Rhines scale L_{Rhi}. The latter is estimated from model variables as L_{Rhi}=\sigma /\beta, where \beta=\partial_y f, while L_r  is given by L_r = min \left[ c_r/|f|, \sqrt{c_r/2\beta}\right] where c_r denotes the first baroclinic Rossby wave speed, which is calculated approximately by c_r \approx \int_{-h}^0 N/\pi dz, where h denotes the local water depth. The inverse eddy time scale \sigma is given by \sigma=|\nabla_h \bar b|/N, which is, by the thermal wind relation, for mid-latitudes identical to the Eady growth rate \sigma=|f| Ri^{-1/2}, where Ri=N^{2}/(\partial_z u^2 +\partial_z v^2 ) denotes the local Richardson number. To prevent a singularity for \sigma for N \to 0 we use \sigma = |f| (Ri + \gamma)^{-1/2} with \gamma>0 which acts effectively as an upper limit for \sigma and consequently for \kappa_I. To prevent a singularity for the time scale \sigma at the equator, f is replaced by \max(|f|, \sqrt{2 \beta c_r}) when computing \sigma. The default values for the parameters of the closure are set to \gamma=200 and c=2. More details about the closure can be found in Eden and Greatbatch (2008) and Eden, Jochum, and Danabasoglu (2009).

8.1.3.3. Near Surface Eddy Flux Parameterization

We have implemented a simplified version of the near-boundary eddy flux parameterization of Ferrari et al. (2008), ((2008), referred to as FMCD08) for the surface boundary in POP. We refer to this scheme as the Near Surface Eddy Flux (NSEF) parameterization. Here, we present a brief summary of its implementation in POP and refer to Danabasoglu, Ferrari, and McWilliams (2008) for further details.

As indicated earlier, the Gent and McWilliams (1990) ((1990), GM90) parameterization was developed for the quasi-adiabatic interior, and it is not valid near the boundaries. Therefore, the usual practice is to taper the effects of parameterized eddy fluxes as the surface (or any other boundary) is approached (see equations (168) and (169)). For example, both \kappa_I and \nu are tapered to zero when the ocean surface is approached using the f_1 taper function in addition to the f_2 large slope taper function using

(174)\kappa_I^* = f_1 f_2 \kappa_I \quad \mathrm{and} \quad
\nu^* = f_1 f_2 \nu;

the asterisk denotes the tapered values of these diffusivities (Large et al. 1997). This near boundary attenuation of the eddy effects is not physically consistent with observational evidence, particularly near the surface where the diabatic mesoscale fluxes may dominate the mixing. FMCD08 includes the effects of these upper-ocean diabatic fluxes and eliminates the ad-hoc, near-surface taper functions.

The NSEF parameterization requires estimation of two vertical length scales: the boundary layer depth (BLD) and the transition layer thickness (TLT). The BLD can exhibit rapid fluctuations during which mesoscale eddies are unlikely to change. In contrast, the mixed layer depth (MLD) represents a lower-frequency, i.e. time-filtered, envelope of the surface boundary layer region. In particular, the mixed layer in our simulations records the maximum depth of the boundary layer after sustained deep mixing events, and it can be considered a measure of the deepest penetration of turbulent mixing in the recent past. Therefore, in our present implementation, we use the MLD instead of the BLD. Following Large et al. (1997), we define the MLD as the shallowest depth where the local, interpolated buoyancy gradient matches the maximum buoyancy gradient between the surface and any discrete depth within that water column. One pass of a grid-scale, five-point spatial filter is applied to the MLD to eliminate any small scale features. For simplicity, we choose to apply no filtering in time.

FMCD08 defines the transition layer as the layer containing all isopycnals within an averaging area and time interval that are intermittently exposed to strong turbulent mixing. Thus, TLT is defined by the range of isopycnals that can be lifted into the boundary layer by subgrid-scale eddy heaving and/or the subgrid-scale (sub-synoptic) variations of the boundary layer. For this purpose, we first calculate a depth D, using the same equation as employed in the f_1 taper evaluation above, i.e.,

(175)D = L_R \vert {\bf L} \vert.

Here, L_R is the Rossby deformation radius, representing the preferred horizontal length scale of the baroclinic eddies, and {\bf L} is the isopycnal slope. For simplicity, the former is determined from L_R = c_0 / \vert f \vert, subject to an additional restriction of 15 km \le \ L_R \ \le 100 km. Here, c_0 \ = 2 ms^{-1} is a typical value for the first baroclinic wave speed. With this prescription, L_R is constant at 100 km equatorwards of 8^\circ latitude and no other equation is used for the equatorial deformation radius. We compute D at each grid point in a vertical column below MLD and search for the shallowest depth where D does not reach MLD. This is equivalent to finding the depth d where,

(176){\rm MLD} < d - D

is satisfied for the first time. TLT is simply obtained from

(177){\rm TLT} = d - {\rm MLD}.

Because we do not impose any maximum slope restrictions in these computations, TLT can become large or even reach the ocean bottom when the column is still very weakly stratified below the mixed layer. Conversely, TLT can vanish if \vert {\bf L} \vert = 0. The sum of MLD and TLT then defines the diabatic layer depth (DLD) over which the upper-ocean eddy fluxes depart from their interior formulas.

The eddy-induced vector streamfunction can be split into its mixed layer, {\bf \Psi}_{\rm ML}, and transition layer, {\bf \Psi}_{\rm TL}, expressions as follows

(178){\bf \Psi}_\mathrm{ML} = \frac{\eta - z}{\eta + \mathrm{MLD}}
                         {\bf \Psi}_\mathrm{o}
  \ {\rm for}\  -\mathrm{MLD} \le z \le \eta

and

(179)\begin{aligned}
{\bf \Psi}_\mathrm{TL} =
     {\left( \frac{z + \mathrm{MLD}}
                  {\mathrm{TLT}} \right)}^2 {\bf \Phi}
   + \left( \frac{\eta - z}{\eta + \mathrm{MLD}}\right){\bf \Psi}_\mathrm{o} \\
 \mathrm{for}\ -\mathrm{DLD} \le z < -\mathrm{MLD}, \nonumber \end{aligned}

respectively. The eddy-induced velocity is obtained using \nabla \times {\bf \Psi}. In the above equations z is the vertical coordinate, positive upward, and \eta is the displacement of the free surface relative to z=0. The two functions {\bf \Psi}_\mathrm{o} and {\bf \Phi} are chosen such that {\bf \Psi} and its vertical derivative are continuous across the base of MLD and the base of TLT. These constraints then yield

(180){\bf \Psi}_\mathrm{o} = \frac{\eta + \mathrm{MLD}}
                             {2(\eta + \mathrm{MLD}) + \mathrm{TLT}}
                        \bigl( 2{\bf \Psi}_I + \mathrm{TLT}\
                                \partial_z {\bf \Psi}_I \bigr)

and

(181){\bf \Phi} = -\frac{\mathrm{TLT}}
                   {2(\eta + \mathrm{MLD}) + \mathrm{TLT}}
              \bigl({\bf \Psi}_I + (\eta + \mathrm{DLD})
                    \partial_z{\bf \Psi}_I \bigr),

where subscript z denotes partial differentiation. In (180) and (181), {\bf \Psi}_I is the interior eddy-induced streamfunction at the base of the transition layer given by the GM90 parameterization, i.e., {\bf \Psi}_I = {\bf \Psi}_\mathrm{GM}(z=-\mathrm{DLD}) where

(182){\bf \Psi}_\mathrm{GM} = - \nu\frac{{\bf z} \times \nabla_h\rho}
                                   {\partial_z \rho}.

In (182) \rho is the local potential density and \nabla_h is the horizontal gradient operator.

In (178), {\bf \Psi}_{\rm o} represents the vector streamfunction value at the base of MLD. We note that {\bf \Psi}_\mathrm{ML} is linear within MLD, going to zero at the ocean surface. This implies constant eddy-induced velocities with no vertical shear within MLD. These eddy-induced velocities, however, must develop a shear in the transition layer to match the interior values. For this purpose, (179) represents the simplest choice for {\bf \Psi}_\mathrm{TL}, which is parabolic, or equivalently linear in z for the eddy-induced velocities.

In our present implementation we neglect any variations of the free surface height in the above equations by setting \eta = 0, because all the horizontal fluxes between z=0 and z=\eta are already neglected in POP due to the linearization of the barotropic continuity equation. The existing discrete implementation of the isopycnal transport parameterization readily subdivides a vertical grid cell into a top and a bottom half. We naturally take advantage of this doubled vertical grid resolution and use the depths of the mid-points of the top and bottom halves as reference depths to determine if a particular half is within MLD, TLT, or the interior. We then evaluate separate {\bf \Psi}_I’s for the appropriate top and bottom halves of the vertical grid cells that straddle the base of the transition layer. These {\bf \Psi}_I’s are also used to compute \partial_z {\bf \Psi}_I along the transition layer - interior interface.

With {\bf \Psi}_I and \partial_z {\bf \Psi}_I now available at z=-{\rm DLD}, {\bf \Psi}_{\rm ML} and {\bf \Psi}_{\rm TL} can be evaluated using (178), (179), (180), and (181). The interior eddy-induced streamfunction values below the base of the transition layer are then determined from (182).

The mesoscale eddy fluxes still mix tracers along isopycnal surfaces in the ocean interior as represented by equation (156). Within the mixed and transition layers, the parameterization for diffusive tracer flux (see FMCD08 and Danabasoglu, Ferrari, and McWilliams (2008)) is given by

(183){\bf F} (\varphi) \ = \ - c(z) \ \kappa_H \nabla_h \varphi \ - \
[1 \ - \ c(z)] \ {\rm {\bf K}} \cdot \nabla \varphi.

In (183), \kappa_H represents the untapered horizontal eddy diffusivity and \nabla_h is the horizontal gradient operator. Of course, the rate of mixing within both the interior and boundary layer remain the same, so \kappa_H has the same value as the interior \kappa_I present in K. The vertical function c(z) is defined by

\nonumber
c(z) =
\begin{cases}
1, & \mathrm{for} \ -\mathrm{MLD}<z\le \eta \\
(z + \mathrm{MLD} + \mathrm{TLT})/\mathrm{TLT},
        &\mathrm{for} -\mathrm{TLT}-\mathrm{MLD}<z \le -\mathrm{MLD}.
\end{cases}

Note that we drop the vertical component of {\bf F} (\varphi) within the mixed layer, because it is found to generate negligible flux divergences at resolutions used in climate models. With the NSEF parameterization, we do not need the f_1 taper any longer, i.e.,

\kappa_I^* \ = \ f_2 \ \kappa_I \ \ \ {\rm and} \ \ \
\nu^* \ = \ f_2 \ \nu,

In addition, f_2 is only applied in the interior below the transition layer. We note that further sensitivity experiments that also use \nu^*=\nu in the ocean interior, i.e., no f_2, lead to numerical instabilities with our present parameter choices. This is the default option used in CCSM4.

8.1.3.4. Vertically-varying Isopycnal and Thickness Diffusivity Coefficients

We provide a vertical dependence of the isopycnal and thickness diffusivity coefficients which vary with the stratification, N^{2}, following Ferreira, Marshall, and Heimbach (2005) and Ferreira and Marshall (2006). The formulation produces diffusivities that are large in the upper ocean and small in the abyss. We regard use of an N^{2} dependency as a convenient way of introducing a surface intensification of the diffusivity coefficients. A brief summary of its implementation in POP2 is presented here and further details are given in Danabasoglu and Marshall (2007).

The vertical variation of \kappa_I and \nu is specified using

(184)\kappa_I \ = \ {\biggl [} \frac{ N^2}{N_{ref}^2 } {\biggr ]} \
{\bigl [} \kappa_I {\bigr ]}_{ref},

and

(185)\nu \ = \ {\biggl [} \frac{ N^2}{N_{ref}^2 } {\biggr ]} \
{\bigl [} \nu {\bigr ]}_{ref},

where N is the local buoyancy frequency computed at the vertical level interfaces, and {\bigl [} \kappa_I {\bigr ]}_{ref} and {\bigl [} \nu {\bigr ]}_{ref} are the constant reference values of \kappa_I and \nu within the surface diabatic layer (SDL). N_{ref} is the reference buoyancy frequency obtained just below the SDL, provided that N^2 > 0 there. Otherwise N_{ref}^2 is the first stable N^2 below the SDL. The ratio N^2 / N_{ref}^2 is set to 1 for all shallower depths, implying no vertical variation of \kappa_I and \nu, particularly within the SDL. Between the depth at which N^2 = N_{ref}^2 and the ocean bottom, we also ensure that

(186)N_{min} \ \le \ \frac{N^2}{N_{ref}^2} \ \le \ 1.0,

where N_{min}(>0) is a specified lower limit. Equation (186) also implies that in statically unstable regions, i.e. N^2 < 0, N^2 / N_{ref}^2 = N_{min}. This approach relies on the enhancement of vertical mixing coefficients to alleviate any local static instability.

The specification of the SDL depends on some other parameterization choices. The internal model default for the SDL is the first model level thickness. If the K-Profile Parameterization (???) is used, then the SDL is set to the diagnosed boundary layer depth. When the NSEF scheme is also specified, the diabatic layer depth (DLD) becomes the SDL. If one uses neither scheme, then the default internal value likely needs to be modified to reflect the bulk mixed layer depth.

The N^2 / N_{ref}^2 profiles are evaluated once per day at the beginning of a day. We avoid any vertical averaging of the profiles to preserve their extrema. Instead, N^2 / N_{ref}^2 computed at the top face of a grid cell is simply used for that grid cell. The resulting diffusivities from (184) and (185) are subject to the usual diffusive numerical stability criterion. In addition, the zero eddy flux boundary condition at the ocean bottom is imposed by setting the diffusivities to zero in the bottom halves of the grid boxes. Thus, the effective diffusivities for these bottom grid cells are further reduced by one half.

One can choose to use (185) only for \nu while using a constant \kappa_I. However, given that there are no strong arguments for setting \kappa_I different from \nu, we use (184) and (185) for both \kappa_I and \nu with the same reference diffusivity values. We note that, at least for the case of constant diffusivities, the theoretical study of Dukowicz and Smith (1997) supports specifying equal values. In the standard nominal \times 1^\circ and \times 3^\circ configurations, we set {\bigl [} \kappa_I {\bigr ]}_{ref} = {\bigl [} \nu {\bigr ]}_{ref} = \kappa_H = 3000 m^2 s^{-1} and {\bigl [} \kappa_I {\bigr ]}_{ref} = {\bigl [} \nu {\bigr ]}_{ref} = \kappa_H = 4000 m^2 s^{-1}, respectively. The default for N_{min} is N_{min}=0.1. Thus, by about 2000-m depth, these reference values are reduced to 300 and 400 m^2 s^{-1}, respectively. We note that \kappa_H is used if the NSEF option is on (default). In addition to the NSEF parameterization, this prescription for vertical variations of the diffusivity coefficients is used in CCSM4.

8.1.4. Parameterization of Mixed Layer Eddies

In principle, the effects of all scales smaller than the model grid scale should be parameterized with subgridscale models. Historically, focus has been on the mesoscale (\approx100 km, e.g. Gent and McWilliams 1990) and finescale (\approx10 m, e.g. ???) parameterizations. However, recent interest has focused on the submesoscale (\approx1 km) physics, because of its ubiquitous restratifying effect on the mixed layer.

Recent realistic simulations (Oschlies 2002; Capet, Campos, and Paiva 2008; ???; ???; Capet et al. 2008), idealized simulations (Haine and Marshall 1998; Thomas 2008; Mahadevan, Tandon, and Ferrari 2009), and theory (Thomas and Lee 2005; Thomas 2005; Boccaletti, Ferrari, and Fox-Kemper 2007; Thomas, Tandon, and Mahadevan 2008; Thomas and Ferrari 2008) have studied and quantified some effects of the submesoscale that may potentially affect the larger scales resolved by climate models. The parameterization implemented in POP follows the theory for one important submesoscale effect–restratification by mixed layer eddies–proposed by Fox-Kemper, Ferrari, and Hallberg (2008 hereafter FFH). The FFH parameterization is validated against idealized submesoscale-resolving simulations (Fox-Kemper and Ferrari 2008) and has been implemented and documented in a number of global climate models (Fox-Kemper et al. 2008; Fox-Kemper et al. 2010). These papers have substantially more detail than the short presentation herein and the reader is encouraged to seek them out.

Mixed Layer Eddies (MLEs) are the finite-amplitude result of ageostropic baroclinic instabilities of the density fronts common throughout the ocean mixed layer (Rudnick and Ferrari 1999; Ferrari and Rudnick 2000; Hosegood, Gregg, and Alford 2006; Boccaletti, Ferrari, and Fox-Kemper 2007). The MLE horizontal scale is near the deformation radius of this mixed layer, which is O(1 km), so these eddies will not be routinely resolved in global- or basin-scale simulations for some time. Nonetheless they do have effects that are important to be included in coarse-resolution models. The primary effect of these eddies is to extract potential energy by restratifying the mean density profile of the upper ocean.

The parameterized restratification is carried out by an eddy-induced overturning streamfunction ({\bf \Psi}), which can be thought of as producing an eddy-induced or quasi-Stokes velocity field ({\bf u}^\ast=\nabla\times{\bf \Psi}). Advection by the eddy-induced velocity provides the eddy fluxes of tracers, including the buoyancy skew flux (\overline{{\bf u'}b'}={\bf\Psi}\times\nabla{\bar{b}}). Buoyancy is the negative density anomaly rescaled to have units of acceleration b\equiv g (\rho_0-\rho)/\rho_0, as appropriate for a Boussinesq fluid. Overlines are used to represent the fields in a coarse-resolution model, that is, one not resolving the submesoscale eddies. The MLE fluxes are assumed to occur in addition to resolved or otherwise parameterized fluxes.

The FFH parameterization is given by

(187)\begin{aligned}

&{\bf \Psi}_0=C_e \frac{H^2 \nabla {\overline b}^z\times{\hat{\bf z}}}{|f|}\mu(z),&\\
&\!\!\!\!\!\!\mu(z)=\max\left\{0,\left[1-\left(\frac{2 z}{H}+1\right)^2\right]\!\!\!\left[1+\frac{5}{21}\left(\frac{2 z}{H}+1\right)^2\right]\right\},&\nonumber\end{aligned}

where H is mixed layer depth, f is the Coriolis parameter, and {\hat{\bf z}} is the unit vertical vector. The overline with subscript z on \nabla {\overline b}^z is understood to be the depth-average of \nabla{\bar{b}} over the mixed layer. The efficiency coefficient C_e is approximately 0.06 (Fox-Kemper, Ferrari, and Hallberg 2008). Note that the parameterization is only nonzero within the mixed layer.

An adaptation to (187) that is suitable and justified in a global coarse-resolution model (Fox-Kemper et al. 2010) is:

(188){\bf \Psi}=C_e \frac{\Delta s}{L_f}\frac{H^2 \nabla {\overline b}^z\times
{\hat{\bf z}}}{\sqrt{f^2+\tau^{-2}}}\mu(z).

\Delta s is an appropriate local gridscale of the coarse model, L_f is an estimate of the typical local width of mixed layer fronts, and \tau is a timescale for mixing momentum vertically across the mixed layer (\approx1-10 days). Hosegood, Gregg, and Alford (2006) suggest L_f is close to the mixed layer deformation radius N H/f, where N is the buoyancy frequency of the mixed layer stratification. For stability, a lower cutoff is used: L_f=\max(N H/|f|, L_{f,min}) where L_{f,min} is chosen in the 1 to 5km range. Likewise, if \Delta s is very large, an upper limit near 1^\circ should be imposed. See Fox-Kemper et al. (2010) for details.

8.1.4.1. Details of Implementation

Coarse Resolution Rescaling: The FFH parameterization (187) depends on the horizontal density gradient. Coarser models have weaker gradients than finer, and sparser observations have weaker gradients than denser. The \Delta s/L_f factor in (188) accounts for the expected dependence of horizontal density gradients on resolution. The vertical buoyancy flux given by (187) scales as

{\overline{w'b'}}\equiv{\bf \Psi}\times\nabla{\bar{b}}\approx{\bf \Psi}
\times\nabla{\overline b}^z\propto \frac{H^2 |\nabla {\overline b}^z\times
{\hat{\bf z}}|^2}{|f|}.

One would like the same vertical buoyancy flux regardless of model resolution, and the \Delta s/L_f factor is a good approximation (Fox-Kemper et al. 2010). Note that the horizontal fluxes are not necessarily accurately rescaled by this method, and the overturning streamfunction may become very large. However, these submesoscale eddies are so small that they generally have negligible horizontal fluxes, and the mixed layer vertical stratification so weak that the parameterized horizontal fluxes are small even when the streamfunction is large. These spurious effects are consequences of the resolution rescaling that should be inconsequential, although the user should keep their spurious nature in mind.

Timestep Limitation: The overturning streamfunction is thus rescaled as \Delta s/L_f. The eddy-induced velocity will go as u^\ast \propto \Delta s/L_f, so the timestep is limited as though the gridscale were L_f instead of \Delta s. In practice u^\ast is usually smaller than other velocities in the model–consistent with the small MLE horizontal fluxes. However, it is possible for the FFH parameterization to require somewhat smaller timesteps in certain conditions on certain grids, especially with fine vertical resolution where w^\ast may determine the timestep.

Steps are taken in POP to ameliorate this timestep limitation for small L_f. 1) A minimum value of L_f, L_{f,min}, is used (usually in the one to five km range). 2) \min[\Delta s, 1^\circ]/\max[L_f,L_{f,min}] constrains the scaleup in very coarse resolution models. As one might expect, the high latitudes where |f| is large are more strongly affected by the choice of cutoff. The mixed layer frontal scale is often obseved to be smaller than 5km. However, the buoyancy frequency, N, in the mixed layer is highly sensitive to other parameterization choices (e.g. maximum diffusivity of mixing in KPP), thus this cutoff is required and should be adjusted along with other model parameters. Generally, the cutoff value should be as small as can be safely integrated. In test cases, a 1km cutoff in POP allows a ML depth nearly 20% shallower in high latitudes than a 5km cutoff.

High-Resolution Usage: The \Delta s/L_f scaling automatically handles regional variations of eddy scales in a high-resolution model. If L_f is resolved in some regions – e.g. where the mixed layer is particularly deep after deep convection – and not in other regions, the scaling adjusts. A just-resolved front (where L_f=\Delta s) has no scale-up and insufficient resolution to form any MLEs. A resolved front with MLEs permitted but not resolved (L_f=4 \Delta s) is boosted by the parameterization, and a well resolved feature (L_f=20 \Delta s) has negligible effects from the parameterization.

Tracers Other than Buoyancy: The rescaling relies on {\overline{w'b'}} \propto |\nabla {\overline b}^z
\times{\hat{\bf z}}|^2, which is appropriate for the scaling of density fronts of the near-surface ocean (???). The FFH parameterization is also used for other tracers, and these fronts will lend similar scaling behavior to those tracers as well. Thus, the rescaling remains appropriate.

Equatorial Regularization: The division by |f| in the scaling for (187) for {\bf\Psi}_0 means that this form cannot be used at the equator. Physically, Boccaletti, Ferrari, and Fox-Kemper (2007) and Fox-Kemper, Ferrari, and Hallberg (2008) note that MLEs are largely geostrophic and thus do not restratify in the same way as the equator is approached. Under typical midlatitude situations, the growth of instabilities rivals the timescales of mixing events and the eddy fluxes are only intermittently interrupted. However, as {\bf\Psi}_0 grows near the equator, frictional constraints enter. The FFH parameterization (187) is therefore adapted in (188) to converge to a friction slumping rate (Ferrari and Young 1997) as the equator is approached. The frictional time constant \tau should be on the order of the timescale between mixing events driven by atmospheric synoptic variability–a few days to a few weeks.

Adjustable Parameters: The scaling (187) has the simulation-based parameter C_e\approx0.06. This constant is an efficiency factor of mixed layer eddies that is apparently constant across different simulations and should not be changed. However, the parameters L_f and \tau are not presently known from observations, theory, or simulations, and therefore they may be tuned to reduce model bias.

Discretization: The discretization of advection by the FFH parameterization is treated identically to that of the advection by the Gent and McWilliams (1990) parameterization. However, since there is no division by N^2 in the FFH parameterization (187), eq:parammod, tapering and slope-limiting are not required.

Regions of Impact: The MLE parameterization is proportional to mixed layer depth squared and horizontal buoyancy gradient. Thus, regions where these parameters are large (e.g., western boundary currents, deep convection sites) will have large effects from the FFH parameterization. However, the parameterization is not intended to be turned off regionally, so it may be active in other regions as well. For example, regions with shallow mixed layers often have weak mixing and weak convection and thus weak vertical tracer transport, so {\overline{w'b'}} from the FFH parameterization may be comparable to these processes.

8.2. Horizontal Viscosity

8.2.1. Laplacian Horizontal Viscosity

The spatial discretization of the standard Laplacian horizontal friction terms with a spatially varying viscosity A_{_M} given by (29) is described in Sec. Horizontal Friction.

8.2.2. Biharmonic Horizontal Viscosity

In analogy to the biharmonic diffusion operator (153), the biharmonic viscosity is constructed by applying the Laplacian-like operator (28) twice. Specifically,

(189)\begin{aligned}
{\cal F}_{Hx}^{(4)}(u_x,u_y) =
{\cal F}_{Hx}^{(2)}(A_M^{(4)}{\cal F}_{Hx}^{(2)}(u_x,u_y),
A_M^{(4)}{\cal F}_{Hy}^{(2)}(u_x,u_y))
\end{aligned}

with a similar expression for {\cal F}_{Hy}^{(4)}. The superscripts (4) and (2) refer to biharmonic and harmonic operators, and {\cal F}_{Hx}^{(2)} is given by (28) with A_M=1. A_M^{(4)} is the biharmonic viscosity, which may be spatially varying, and should be negative for positive-definite dissipation of kinetic energy. Here the biharmonic viscosity is sandwiched between the two applications of the harmonic operator. An alternate approach (not implemented in the code) is to include the square-root of -A_M inside the second derivative of each harmonic operator (Griffies and Hallberg 2000).

8.2.3. Anisotropic Horizontal Viscosity

An anisotropic formulation of the friction terms in the momentum equation was first introduced into an ocean GCM by Large et al. (2001). This was implemented in a model with a spherical polar grid, and their formulation was specifically tied to that coordinate system. A more general formulation of anisotropic viscosity that can be applied in any general orthogonal coordinate system was developed by Smith and McWilliams (2003) ((2003), hereafter SM), and is implemented in the POP code. The friction operator is formulated as the divergence of a viscous stress tensor which is linearly related to the velocity gradients. The general anisotropic formulation of the stress involves four independent viscous coefficients, as well as a unit vector \hat{\bf n} that specifies a preferred horizontal direction which breaks the transverse isotropy. The coefficients and the direction vector may vary in both space and time. A modification of the original formulation has been implemented by Jochum et al. (2008) that changes the values, but not the structure of viscosity. Equations (195) and (197) are modified from the SM version to allow for a better representation of equatorial and high latitude boundary currents.

SM discuss two different reduced forms involving only two viscous coefficients, A and B. In both of these forms, the friction operator in Cartesian coordinates has the following approximate form if the x-coordinate is aligned with \hat{\bf n}:

(190)\begin{aligned}
F_{x} &= A\partial^{2}_{x}u+B\partial^{2}_{y}u \nonumber \\
F_{y} &= B\partial^{2}_{x}v+A\partial^{2}_{y}v.
\end{aligned}

Thus, the friction acts to diffuse momentum along the direction \hat{\bf n} with viscosity A, and perpendicular to \hat{\bf n} with viscosity B; this is true for the general operator, even when the coordinates are not aligned with \hat{\bf n}.

The three independent elements \sigma_{11}, \sigma_{22} and \sigma_{12} of the symmetric transverse stress tensor \sigma_{ij} are linearly related to the elements of the rate-of-strain tensor \dot{e}_{ij} (where the subscripts 1 and 2 refer to the x and y coordinates, respectively). In the first form discussed by SM, which is currently implemented in the code, the stress and strain-rate are related by

(191)\begin{aligned}

\left( \begin{array}{c}
       {\bf\sigma}_{11} \\
       {\bf\sigma}_{22} \\
       {\bf\sigma}_{12} \end{array}\right)
& =
\left( \begin{array}{rrr}
       A & -B & 0 \\
      -B &  A & 0 \\
       0 &  0 & B \end{array} \right)
\left( \begin{array}{c}
        \dot{e}_{11} \\
        \dot{e}_{22} \\
       2\dot{e}_{12} \end{array} \right) \\
& + (A-B)n_{1}n_{2}
\left( \begin{array}{ccc}
      -2n_1 n_2 &      2n_1 n_2 & n_1^2 - n_2^2 \\
       2n_1 n_2 &     -2n_1 n_2 & n_2^2 - n_1^2 \\
  n_1^2 - n_2^2 & n_2^2 - n_1^2 & 2n_1 n_2
\end{array} \right)
\left( \begin{array}{c}
        \dot{e}_{11} \\
        \dot{e}_{22} \\
       2\dot{e}_{12} \end{array} \right) \nonumber\end{aligned}

where n_1 and n_2 are the components of \hat{\bf n} along the x and y coordinates, respectively. In the second two-coefficient form, the relation between the stress and strain-rate is also given by (191), except that the first matrix inside the brackets on the right is replaced by

\begin{aligned}
\left( \begin{array}{rrr}
       {\scriptstyle\frac{1}{2}}(A'+B')
& -{\scriptstyle\frac{1}{2}}(A'+B')
& 0 \\
      -{\scriptstyle\frac{1}{2}}(A'+B')
&  {\scriptstyle\frac{1}{2}}(A'+B')
& 0 \\
       0  & 0 & B' \end{array} \right).
\nonumber\end{aligned}

The second form is actually independent of the horizontal divergence, and can be written in the more compact form:

(192)\begin{aligned}

\left( \begin{array}{c}
       \sigma_{_T} \\
       \sigma_{_S} \end{array}\right)
=
\Biggl[
& \left( \begin{array}{cc}
       A+B & 0   \\
         0 & 2B \end{array} \right) \\
& + 2(A-B)n_{1}n_{2}
\left( \begin{array}{cc}
            -2n_1 n_2 & n_1^2 - n_2^2   \\
       n_1^2 - n_2^2  & 2n_1 n_2  \end{array} \right)
\Biggr]
\left( \begin{array}{c}
       \dot{e}_{_T} \\
       \dot{e}_{_S} \end{array}\right). \nonumber \end{aligned}

where \sigma_{_T} = \sigma_{11}-\sigma_{22}, \sigma_{_S} = 2\sigma_{12}, \dot{e}_{_T} = \dot{e}_{11}-\dot{e}_{22} and \dot{e}_{_S} = 2\dot{e}_{12}. The two forms are equivalent up to terms proportional to the horizontal velocity divergence \nabla\cdot{\bf u} = \dot{e}_{11}+\dot{e}_{22} which is very small compared to typical velocity gradients in geophysical flows, as discussed by SM. The first form was essentially constructed so that it exactly produces the Cartesian friction operator (190), whereas the second form was derived as the limit \nabla\cdot{\bf u}\rightarrow 0 of the general four-coefficient form of the stress. For this reason, the second form has a more sound physical basis, however, the results using the two forms should in most cases be very similar.

In both forms the stress is invariant under a rotation of \hat{\bf n} by 90^\circ. Currently there are three options in the code for the field of unit vectors \hat{\bf n}: (1) aligned along the flow direction (\hat{\bf n} = {\bf u}/|{\bf u}|); (2) aligned with due east (or, equivalently, with any of the four cardinal directions); and (3) aligned with the local grid coordinates (in polar coordinates this is equivalent to (2), but it is different in general orthogonal coordinate systems such as the displaced-pole grids described in Sec. Global Orthogonal Grids). With any of these choices, the isotropic limit is obtained by equating the two coefficients A=B. The first form reduces to the original anisotropic formulation of Large et al. (2001) when \hat{\bf n} is oriented eastward.

The spatial discretization of the anisotropic friction operator is described in detail by SM, and will only be briefly reviewed here. The basic idea of this approach is to take advantage of the fact that the friction operator can be written as the functional derivative of the area-integrated energy dissipation rate. Using the no-slip lateral boundary conditions on the velocity, it can be shown that the domain-averaged dissipation of kinetic energy due to friction is given by

\begin{aligned}
\int da\, {\bf u}\cdot{\bf F}
 = \int da\, {\cal D}\end{aligned}

where {\bf F} is the friction vector and {\cal D} is the energy dissipation rate:

\begin{aligned}
{\cal D} = \frac{1}{2}\,\dot{e}:\sigma \:.\end{aligned}

To ensure positive-definite dissipation of kinetic energy, the viscous coefficients must be chosen to satisfy {\cal D}\ge 0. It can be shown that in the first form this leads to the constraint A\ge B\ge 0, whereas the second form has the less restrictive constraint A\ge 0,\,B\ge 0. The i^{\mathrm{th}} component of the friction is given by

\begin{aligned}
F_i=-\frac{\delta}{\delta u_i}
\int dV\: {\cal D}\end{aligned}

where \frac{\delta}{\delta u_i} denotes a functional derivative with respect to the i^{\mathrm{th}} component of the velocity. The procedure for deriving the discrete operator is to first discretize the functional, making sure that the discrete form is positive definite, and then take its derivatives with respect to the velocity components at a given point on the computational grid; this yields a friction operator at that point that is guaranteed to dissipate kinetic energy. As was done in the case of the GM operator described in Sec. The Gent-McWilliams Parameterization, the method of discretizing the functional is to subdivide each horizontal cell into four subcells (see Figure7.1. , with the central point interpreted as a U-point). This allows all quantities (such as viscous coefficients, gradients of velocities, etc.) required to construct the dissipation rate to be assigned a unique value in each subcell. In this manner the components of the strain-rate tensor are computed in each subcell, and from these the stresses are computed using (191). The specific details of the derivation are given by SM, and we will only quote the final result here. The discrete friction operator has the form

(193)\begin{aligned}
F_x = &\frac{1}{V}[
\overline{h_2\sigma_{11}}^e - \overline{h_2\sigma_{11}}^w
+ \overline{h_1\sigma_{12}}^n - \overline{h_1\sigma_{12}}^s \nonumber \\
&+\frac{1}{2}\,(h_2^n k_2^n \overline{h_1\sigma_{12}}^n
+h_2^s k_2^s \overline{h_1\sigma_{12}}^s) \nonumber \\
&-\frac{1}{2}\,(h_1^e k_1^e \overline{h_2\sigma_{22}}^e
+h_1^w k_1^w \overline{h_2\sigma_{22}}^w)]
\nonumber \\
\nonumber \\
F_y = &\frac{1}{V}[
\overline{h_1\sigma_{22}}^n - \overline{h_1\sigma_{22}}^s
+ \overline{h_2\sigma_{12}}^e - \overline{h_2\sigma_{12}}^w \nonumber \\
&+\frac{1}{2}\,(h_1^e k_1^e \overline{h_2\sigma_{12}}^e
+h_1^w k_1^w \overline{h_2\sigma_{12}}^w) \nonumber \\
&-\frac{1}{2}\,(h_2^n k_2^n \overline{h_1\sigma_{11}}^n
+h_2^s k_2^s \overline{h_1\sigma_{11}}^s)]
\end{aligned}

where F_x and F_y are the components of the friction operator along the two general orthogonal coordinate directions x and y. The factors h_1 and h_2 are distances between grid points along the two coordinate directions, and k_1 and k_2 are metric factors proportional to the derivatives of h_2 and h_1 in the x and y directions, respectively. The overbars denote an average over the four subcells surrounding a given face (denoted by the superscript n, s, e, or w, for the north, south, east and west faces), and V is the full cell area. For a more detailed description of the quantities in (193), see the appendix of SM.

The viscous coefficients are assumed to have a constant value within a given full cell, i,e., all four subcells within each full cell have the same values of the coefficients A and B. There are currently three options in the code for the form of the viscous coefficients: (1) constant values of A and B in both space and time, (2) spatially varying but constant in time; and (3) Smagorinsky-type nonlinear coefficients that depend on the local deformation rate and hence vary in both space and time. In options (2) and (3), the viscosities are tapered if they exceed one-half the viscous CFL limit. Specifically, both A and B are tapered as follows:

(194)\begin{aligned}
A &= \min(A,\frac{1}{2}A_{\mathrm{cfl}}) \nonumber \\
A_{\mathrm{cfl}} &= \frac{(dx)^2 + (dy)^2}{4\,dt}
\end{aligned}

where dt is the model momentum timestep, and dx and dy are the grid spacings in the two coordinate directions. For option (1), A and B are simply set to different, constant values. The model does not enforce the limit, if they exceed A_{\mathrm{cfl}}; only a warning message is printed. If A < B, the integration is terminated with an error exit.

The spatial dependence of the coefficients is computed following Jochum et al. (2008), which also demonstrates their importance for climate. The two primary design criteria are to have viscosity at values appropriate for the parameterization of missing meso-scale eddies wherever possible and to use other values only where required by the numerics. The parallel coefficient A_1 is given by

(195)A_1 = \max \biggl \{ {c_3 \beta dx^3 \exp(-p(x)^2),
      A_{eddy}} {\biggr \}},

The first part represents a numerical, Munk viscous western boundary constraint with

(196)p(x) = L_M^{-1}  ~ \max (0,x-x_N)

where p(x) causes A_1 to fall off as fast as possible away from the western boundaries. Here, x_N is the zonal coordinate of the N^{th} grid point east of the nearest western boundary and L_M is a length scale. Viscosity is not similarly increased near zonal and eastern boundaries because doing so does not reduce numerical noise. Also, c_3 (\approx 0.2) is a dimensionless coefficient, and \beta is the meridional gradient of the Coriolis parameter. A_{eddy} is a physical, lower bound set to account for all the missing mesoscale eddy activity. The perpendicular coefficient B_1 is constructed as

(197)B_1 = \max {\biggl \{} c_3\ \beta \ dx^3 \ exp(-p(x)^2),\ B_{eddy}
\ {\bigl (} 1 + c_2 (1 - cos [2\phi^\prime] ) {\bigr )} {\biggr \}}.

The first part represents again the Munk viscous western boundary constraint. In the second part, \phi^\prime = 90^\circ
\ min (\vert \phi \vert, \phi_I ) / \phi_I < 90^\circ and B_{eddy} represents a physical lower bound. This form allows B_1 to be small at the equator and increase poleward for latitude between \pm \phi_I, using the dimensionless coefficient c_2. A preferred option is to set (1+2c_2)  =  A_{eddy}/B_{eddy}, so that B_1 becomes equal to A_1 poleward of \phi_I. In most applications, however, B_{eddy} is set equal to A_{eddy}.

To obtain the final distributions of A and B, the viscosities are tapered if they exceed half of the viscous CFL limit:

(198)\begin{aligned}
A &= \min (A_1, \frac{1}{2} A_{\mathrm{cfl}}), \nonumber \\
B &= \min (B_1, \frac{1}{2} A_{\mathrm{cfl}}).
\end{aligned}

8.2.3.1. Smagorinsky Nonlinear Viscous Coefficients

Jochum et al. (2008) shows that Smagorinsky-type nonlinear viscosity leads to excessive viscosity along boundaries and the equator. However, the code is still set up to optionally allow for Smagorinksy-type nonlinear dependence of the coefficients on the deformation rate (Smagorinsky 1993) if the anisotropic functional discretization described in Sec. Anisotropic Horizontal Viscosity is used. This option is not available with the simple Laplacian-like friction operator given by (29). An isotropic form with the same Smagorinsky-type nonlinear viscosity can be obtained by simply setting the two anisotropic viscous coefficients A and B equal to one another.

The deformation rate \dot{D} is proportional to the norm of the strain-rate tensor:

(199)\begin{aligned}
\frac{1}{2}\dot{D}^2 = \|\dot{e}\|^2
= \dot{e}_{11}^2 + \dot{e}_{22}^2 + 2\dot{e}_{12}^2\approx
\frac{1}{2}(\dot{e}_{_T}^2+\dot{e}_{_S}^2)
\end{aligned}

This is particularly easy to evaluate because the strain-rate tensor is already computed in the code. Specifically, the nonlinear coefficients A and B are given by

(200)\begin{aligned}
A &= \min\left[\max\left(C_{_{A}} \dot{D}\,
ds^2,V_{_{A}}\,ds\right),\: \frac{1}{2}A_{\mathrm{cfl}}\right]
\nonumber \\
B &= \min\left[\max\left(C_{_{B}} \dot{D}\,
ds^2,V_{_{B}}\,ds\right),\:\frac{1}{2}A_{\mathrm{cfl}}\right]
\end{aligned}

where ds=\mathrm{min}(dx,dy) is the minimum grid spacing in the two coordinate directions. C_{_{A}} and C_{_{B}} are dimensionless coefficients of order 1, and V_{_{A}} and V_{_{B}} are velocity scales associated with the grid Reynolds number which determine a minimum background viscosity in regions where the nonlinear viscosity is too small to control grid-point noise. Typically V_{_{A}} and V_{_{B}} are order 1 cm s^{-1}. A_{\mathrm{cfl}} is the maximum viscosity allowed by the viscous CFL limit, see (194). The advantage of the nonlinear viscosities (200) is that they selectively apply the friction operator in regions of strong shear.

In moderate resolution configurations of POP, the Smagorinsky coefficients have been modified in a similar manner to the anisotropic coefficients in Large et al. (2001) described above. Again the reason is to control numerical noise in the solutions. Equation (200) has been modified such that

(201)\begin{aligned}
A &= \min \left[ \max \left( C_A F_A(\phi)\dot{D}ds^2, A_1 \right),
     \ \frac{1}{2} A_\mathrm{cfl} \right], \nonumber \\
B &= \min \left[ \max \left( C_B F_B(\phi)\dot{D}ds^2, B_1 \right),
     \ \frac{1}{2} A_\mathrm{cfl} \right],
\end{aligned}

with

(202)\begin{aligned}
F_A (\phi) & = & 1, \nonumber \\
F_B (\phi) & = & 0.02 \ \ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ {\rm for} \
               \vert \phi \vert \le 20^\circ, \cr
           & = & 1 \ - \ 0.98 \ exp( - (\phi - 20)^2/98) \ \ {\rm for} \
               \vert \phi \vert  > 20^\circ.
\end{aligned}

Here, ds={\mathrm{min}}(dx,dy) is the minimum grid spacing in the two coordinate directions. C_A and C_B are dimensionless coefficients of order 1 whose latitudinal dependencies are given by F_A and F_B, respectively. The positive-definite dissipation of kinetic energy is enforced by the last equation in (198). Where these nonlinear viscosities are too small, A_1 and B_1, with the same formulations as in the previous section, provide lower limits to control any grid-point noise.

8.3. Vertical Mixing

Vertical mixing is treated using the spatial discretizations presented in sections Vertical Tracer Diffusion and Vertical Friction. The implicit time discretizations are presented in Secs. Implicit Vertical Diffusion, Tridiagonal Solver, and Semi-Implicit Treatment of Coriolis and Vertical Friction Terms. POP contains three different parameterizations for computing the vertical diffusivity \kappa (VDC in the code) and viscosity \mu (VVC in the code): constant coefficients, a parameterization based on the gradient Richardson number by Pacanowski and Philander (1981) and the K-profile parameterization (KPP) of (???).

8.3.1. Convection

Convective instability can be handled in two ways within the code. If implicit vertical mixing is used, convection is generally treated by assigning a very large diffusivity and viscosity when the density profile between two cells in the vertical direction is statically unstable. The diffusivity and viscosity in each case is determined at run time through namelist inputs.

If explicit vertical mixing is used, convection is usually treated using convective adjustment. During convective adjustment, multiple passes through the water column are performed to check for stability and to adjust tracer values. Each pass in turn consists of two sweeps. The first sweep checks odd-valued vertical levels; the second checks even-valued levels. At each level k, the stability between level k and level k+1 is checked (see below for stability criterion). If the two cells are found to be unstably-stratified, the tracer in each level is adjusted using

(203)\phi_k = \phi_{k+1} =  \frac{{(dz_k \phi_k + dz_{k+1}\phi_{k+1})}}
                            {{2dz_{k+\frac{1}{2}}}}.

If tracer acceleration is being used, the thicknesses dz are adjusted by the acceleration factor

(204)\begin{aligned}
dz_k^* &= dz_k/\gamma_k \\
dz_{k+\frac{1}{2}}* &= \frac{1}{2}(dz_k^* + dz_{k+1}^*).
 \end{aligned}

Instability can be determined using one of two methods. In both the Richardson number parameterization and KPP, the Richardson number is computed and the column is unstable when the Richardson number is less than zero. The calculation of Richardson number differs slightly between the two and is shown in each section below.

If the Richardson number is not computed, stability is determined by adiabatically displacing a water mass from the current vertical level k to the level k+1 just below. The density of the parcel after this displacement, denoted \rho^*_k, is computed by calling the equation of state with the temperature and salinity from level k, but evaluating the equation of state at level k+1. If the density after the displacement is greater than the actual density at the level k+1, the column is unstable.

8.3.2. Constant Vertical Viscosity

If this option is chosen, the vertical diffusivity VDC and vertical viscosity VVC are simply constant at all levels and for all tracers. The scalar constant values are determined at run time through namelist input.

8.3.3. Richardson Number Dependent Mixing

In this parameterization, the vertical diffusivity and viscosity are functions of the Richardson number,

(205)\begin{aligned}
Ri & = -g \frac{{(\rho^*_k - \rho_{k+1})}}
       {{\frac{1}{2}(dz(k) + dz(k+1))}}
       {\left(\frac{{\partial V}}{{\partial z}}\right)^{-2}} \nonumber \\
   & = \frac{{-gdz_{k+\frac{1}{2}}(\rho_k^* - \rho_{k+1})}}
          {{((u_k - u_{k+1})^2 + (v_k - v_{k+1})^2 + \epsilon)}}
 \end{aligned}

where the velocities are evaluated at tracer points and \epsilon is simply a small number to avoid dividing by zero.

The particular functional forms for the diffusivity \kappa (VDC) and viscosity \mu (VVC) used in the code are

(206)\begin{aligned}
\kappa &= \kappa_w + (\mu_w +
                 Rich\_mix/(1 + 5Ri)^2)/(1 + 5Ri) \\
\mu &= \mu_w + Rich\_mix/(1 + 5Ri)^2
 \end{aligned}

where the background values \kappa_w and \mu_w and Richardson mixing coefficient Rich\_mix are determined at run time through namelist input. Note that the Richardson number used for the viscosity \mu has been averaged from tracer points to velocity points.

8.3.4. The KPP Boundary-Layer Parameterization

The full KPP parameterization is detailed in (???), hereafter LMD, so only a general overview plus specifics of the current POP implementation and later modifications (see also Appendix A of Danabasoglu et al. (2006)) are described here. The scheme provides all the coefficients required to compute the vertical mixing contributions to both {\cal F}_V (\alpha) in (7), and {\cal D}_V (\varphi) in (11). As such, it should not be used in conjunction with any other parameterization of vertical mixing. Implicit vertical mixing is required to solve (7) and a modified version of (13); namely,

(207){\cal D}_V (\varphi) = \frac{\partial}{\partial z} \kappa
    \left( \frac{\partial}{\partial z}\varphi - \gamma_{\varphi}\right),

where the additional non-local transport term, \gamma_{\varphi} is non-zero only in convective (unstable) forcing, and above a diagnosed boundary layer depth, h. In the ocean interior below h, the viscosity, \mu_{IN}, and diffusivity, \kappa_{IN}, represent very different physics than their boundary layer, z < -h, counterparts, denoted \mu_B and \kappa_B, respectively.

The inputs to the scheme on the tracer grid are the surface forcing and vertical profiles of temperature and salinity. The specific forcings are the surface friction velocity, u^*, both the solar, B_{sol}, and non-solar, B_{ns}, buoyancy fluxes, and the kinematic surface tracer fluxes, \overline{w \varphi}_o. From the equation of state and its buoyancy form, b(\Theta,S,p), derived profiles are the thermal and haline Boussinesq coefficients, \alpha_T and \beta_S, the local buoyancy difference at each interface, \Delta b_{k+.5}, and the excess buoyancy of first level water when moved down to each grid level,

(208)\Delta B_k = b(\Theta_1, S_1, -z_k) ~-~ b(\Theta_k, S_k, -z_k) ~.

The local vertical shear at each interface is squared before being averaged onto the tracer grid, as Sh^2_{k+.5}. In contrast, the square of velocity differences with the first level on the tracer grid are taken as the maximum from the four neighboring velocity points:

(209)|\Delta V|^2_k = {\max}_{\mathrm{(over\ 4\ velocity\ points)}}
                 \left|(u_1,v_1) - (u_k,v_k)\right|^2.

For vertical resolutions of 1m or less, excessive sensitivity can be avoided if \Theta_1,  S_1, u_1 and v_1 in the above equations are replaced by averages over the upper 10% of the boundary layer (LMD). A local gradient Richardson number, Ri_g, buoyancy frequency, N, and density ratio, R_{\rho}, on the tracer grid are computed as

(210)\begin{aligned}
N^2  &= \Delta b_{k+\frac{1}{2}}/(z_k - z_{k+1}) \nonumber \\
Ri_g &= N^2/Sh^2_{k+\frac{1}{2}}                 \nonumber \\
R_{\rho} &= \alpha_T \partial_z\Theta/(\beta_S \partial_z S).
\end{aligned}

It is then a five step process to complete the algorithm, as described below. First, the interior mixing coefficients are computed at all model interfaces, on the tracer grid, as if there were no boundary layer scheme. These are the sum of individual viscosities and diffusivities corresponding to a number (currently up to four) of different physical processes:

(211)\begin{aligned}
\mu_{IN}    &= \mu_w + \mu_s + \mu_c  \nonumber \\
\kappa_{IN} &= \kappa_w + \kappa_s + \kappa_c + \kappa_d
\end{aligned}

The first coefficients on the right hand side are background values. With the exception of the surface mixed layer and the bottom boundary layer, the parameterizations of which are documented below, the ocean is mostly adiabatic. However, even there diapycnal mixing happens due to internal waves and other mechanisms. In POP this is represented by the background diffusivity \kappa_w and the background viscosity \mu_w. Their vertical variation has the general form

(212)\begin{aligned}
\kappa_w &= vdc1 + vdc2 \tan^{-1} ((|z| - dpth) linv), \nonumber \\
\mu_w    &= Pr_w \kappa_w
\end{aligned}

where vdc1 equals the vertical diffusivity at |z|=D, vdc2 is the amplitude of variation, linv is the inverse length scale of the transition region, and dpth is the depth where diffusivity equals vdc1. The form allows for an increase in diffusivity with depth, as a crude parameterization of the observed increase in deep mixing over rough topography. Pr_w is only poorly constrained by observations and thus set to a constant value of 10. \kappa_w, too, used to be only weakly constrained by observations, but since the release of CCSM3 in 2004 it became clear that there is not only a latitudinal structure to \kappa_w, but also that variations to it can have a dramatic impact on climate (Jochum and Potemra 2008). Thus, the value of \kappa_w has been changed from constant values of typically 0.1 cm^2/s to 0.17 cm^2/s (bckgrnd_vdc1 + bckgrnd_vdc_eq = 0.16 + 0.01 = 0.17) mostly everywhere. The regions with different values are the Banda Sea (bckgrnd_vdc_ban = 1.0 cm^2/s), the latitude bands around 30^{\circ}N/S (bckgrnd_vdc1 + bckgrnd_vdc_eq + bckgrnd_vdc_psim = 0.3 cm^2/s), and the equator (bckgrnd_vdc_eq = 0.01 cm^2/s). The theories and observations that led to this new structure, as well as their impact on climate, are discussed in Jochum (2009).

The viscosity and diffusivity associated with shear instability mixing are equal (Pr_s = 1) and parameterized as a function of Ri_g:

(213)\mu_s = \kappa_s =  \kappa_s^0 \left[1-
                    \left(\frac{Ri_g}{Ri_0}\right)^2 \right]^3,
     \quad  0 < Ri_g  <  Ri_0  .

For an unstable profile with negative Ri_g the coefficients remain constant at \kappa_s^0 (rich_mix (cm^2/s) in namelist), and are zero for all Ri_g \ge Ri_0. This function falls most rapidly near Ri_g = 0.4 Ri_0, where it approximates the onset of shear instability. In this neighborhood rapid changes in Ri_g can cause instabilities to develop in the vertical, but these are largely controlled by vertically smoothing Ri_g profiles. A 1-2-1 smoother is can be applied repeated a specified number of times.

Convective instability in the interior ocean is relieved by setting the mixing coefficients, \mu_c and \kappa_c, to large values whenever the density profile is unstable, N<0. Otherwise, they are set to zero.

Second, double diffusion processes have the potential to significantly enhance diffusivities, with R_{\rho} the governing parameter. Their effects on viscosity are ignored in (211), because they are not well known and because the large value of \mu_w, relative to \kappa_w in (211) means that they are unlikely to be dominant. In the salt fingering regime (destabilizing salinity profile and 1 < R_{\rho} < R_{\rho}^0 = 2.55), the diffusivity is a fit to observational estimates (St.Laurent and Schmitt 1999)

(214)\kappa_d = \kappa_d^0 \left[1-\frac{{R_{\rho}-1}}{{R_{\rho}^0 -1}}\right]^3.

The values of \kappa_d^0 are internally set to 1 cm^2/s for salt and 0.7 cm^2/s for heat. Diffusive convective instability occurs where the temperature is destabilizing and 0 < R_{\rho} < 1. For temperature

(215)\begin{aligned}
\kappa_d = \mathrm{VISCM} \times .909
     \exp\left(4.6\exp\left[-.54\left(R_{\rho}^{-1} - 1\right)\right]\right),
\end{aligned}

where VISCM is molecular viscosity. Multiplying this diffusivity by a

(216)\mathrm{factor} =
\begin{cases}
\left(1.85 - 0.85R_{\rho}^{-1}\right) R_{\rho} & 0.5\leq R_{\rho}<1 \\
0.15  R_{\rho}                                 & R_{\rho} < 0.5,
\end{cases}

gives diffusivities for other tracers, including salt.

Third, the diagnostic boundary layer depth, h, is determined on the tracer grid. A bulk Richardson number relative to the surface is defined at each vertical level as

(217)Ri_b = \frac{ -z_k \Delta B_k }{|\Delta V|^2_k  + V_t^2}

LMD gives the rationale and expression for the shear contribution from unresolved turbulence, V_t / d. The buoyancy frequency in this expression is taken from (210) at the interface below the boundary layer depth. Furthermore, the factor which accounts for smoothing of the buoyancy profile at the entrainment depth is an empirical function of this frequency, and ranges linearly from 2.1 at zero to 1.7 at 0.002 s^{-1} and all higher frequencies. The boundary layer depth is equated to the shallowest depth where Ri_b equals an internally specified critical value, Ri_{cr}=0.3 (LMD), using a quadratic, rather than linear, interpolation whenever possible. With stabilizing surface forcing, options exist (namelist flag lcheckekmo = .true.) to limit h to be no greater than either the Ekman depth (0.7 u^{*}/ f), or the Monin-Obukhov depth (L_{mo} = u^{*3}/ (vonk~~B_f)), where vonk is the von Karman constant, and the surface buoyancy flux, B_f, is B_{ns} plus a fraction of B_{sol} (LMD). The grid level immediately below h is denoted as k_b.

Fourth, the boundary layer mixing coefficients are computed on the tracer grid and replace the interior values from step 1, for 1.5 \leq k \leq k_b-\frac{1}{2}. The analytic expression is

(218)\begin{aligned}
\mu_B(\sigma)    &= vonk \  w_{\alpha}hG(\sigma) \nonumber \\
\kappa_B(\sigma) &= vonk \  w_{\varphi}hG(\sigma),
\end{aligned}

where \sigma = -z / h varies from 0 to 1 over the boundary layer. The turbulent velocity scales, w_{\alpha} and w_{\varphi} are usually proportional to u^*, but become proportional to the convective velocity scale as u^* \rightarrow 0 in convective forcing (see LMD for details). The shape function is a cubic polynomial whose coefficients are chosen such that G(0) = 0, fluxes vary linearly near the surface, and interior (excluding the convective contributions) and boundary layer coefficients, and their first vertical derivatives, are continuous at z=-h. An inherent bias to shallow boundary layers is ameliorated by making the coefficients at the k_b-\frac{1}{2} interface linear combinations of the interior values at this interface, and the boundary layer values at both the interface and at the nearest higher level, k_b-1 (LMD). This feature is less important when quadratic interpolation is used to find the boundary layer depth.

In convective (unstable) forcing situations the nonlocal term is non-zero, and evaluated as

(219)\gamma_{\varphi} = C_{\varphi}
                   \frac{\overline{w\varphi_o}}{u^* h},

where the constant C_{\varphi} is prescribed as in LMD.

Fifth and finally, the viscosity is averaged from the tracer grid to velocity points.

8.3.5. Tidally Driven Mixing

The tidally driven mixing parameterization implemented in CCSM POP is described in Jayne (2009).  There are multiple parameterizations of mixing phenomena in the model (as well as numerical mixing), however here we are only concerned with the vertical mixing in the interior away from the surface boundary layer. The tidally driven vertical mixing parameterization of Laurent, Simmons, and Jayne (2002) was implemented in a similar fashion as in Simmons et al. (2004). This mixing is added to the interior mixing part of the KPP vertical mixing scheme as another component, \kappa_t, in equation (211):

(220)\kappa_t = \kappa_w + \frac{\Gamma\varepsilon}{N^2} =
\kappa_w + \frac{\Gamma q E_f(x,y)F(z)}{\rho N^2}

We note that \kappa_t essentially replaces \kappa_w in equation (211), because \kappa_w is already included in \kappa_t.

The parameterization for the vertical diffusivity is related to the turbulent dissipation of internal wave dissipation using the Osborn (1980) relation and depends on the stratification, N. The mixing efficiency of turbulence is set by \Gamma and is taken to be the canonical value of \Gamma = 0.2 (Osborn 1980). The tidal dissipation efficiency is given by the parameter, q = \frac{1}{3}, and represents the part of the internal wave energy flux, E(x,y), that is estimated to be dissipated locally (Laurent and Garrett 2002). The rest of the internal wave energy (1-q = 2/3) is presumed to radiate to the far field and contribute to the background internal wave field (Garrett and Munk 1975). The vertical structure function, F(z), models the distribution of the turbulent mixing in the vertical, and is implemented as a simple exponential decaying upwards away from the bottom, with a vertical scale of \zeta = 500m (Laurent and Nash 2004):

(221)F(z) = \frac{e^{(z-H)/\zeta}}{\zeta (1-e^{-H/\zeta})}

so it is normalized such that vertical integral over the water column is unity (Laurent, Simmons, and Jayne 2002). In the POP code implementation, the vertical exponential function is summed numerically for the normalization.

The constant background value of \kappa_w is an attempt to model the low level of mixing in the ocean in areas away from topography (Gregg 1987). Exactly what contributes the mechanical energy required to maintain order 0.1\; \mathrm{cm}^2\mathrm{s}^{-1} background diffusivity is less clear. However, it is the order of magnitude expected from the self-interaction of the background Garrett-Munk internal wave spectrum (Gregg 1987; Gregg 1989; Hibiya, Nagasawa, and Niwa 2006).

The vertical viscosity is calculated from the vertical diffusivity assuming Pr_w = 10, as is typically done in other parameterizations (e.g. KPP). In the limit of N \rightarrow 0 (or becoming negative), both the vertical diffusivity and viscosity are capped at 100\; \mathrm{cm}^2\mathrm{s}^{-1}. This is a departure from the original Simmons et al. (2004) who instead impose a lower limit on N^2 of 10^{-8}\mathrm{s}^{-2} (additionally convective mixing was used to limit negative stratification regions). In addition, both the diffusivity and viscosity coefficients within the two bottom-most grid cells are checked and modified, if necessary, to ensure that they do not diminish with depth near the bottom topography. Such decreasing-with-depth diffusivities result from large N^2 in a few isolated grid points, likely due to advection errors. Modifications of these diffusivity coefficients, e.g., capping the diffusivities, change the energy consistency of the scheme, and so should be chosen carefully.  One other change from the original Simmons et al. (2004) implementation was that the parameterization was applied everywhere, whereas they arbitrarily limited it to depths greater than 1000m.  The conversion of tidal energy into internal waves that occurs in the deep ocean also occurs along the shelf break, so there is no reason to limit it to the abyssal ocean.

The internal wave energy map, E(x,y), is derived in a similar manner as in Jayne and Laurent (2001) from a barotropic model of the tides utilizing a parameterization of the conversion of barotropic tidal energy into internal waves. It is read from a user-supplied input file. The essential goal of the parameterization is to represent the momentum exchange between the barotropic tides and the unrepresented internal waves induced by the tidal flow over rough topography in a stratified ocean. In the parameterization of Jayne and Laurent (2001), the conversion of barotropic tidal kinetic energy is given by:

(222)E(x,y) = \frac{1}{2}\rho_0kh^2N\mathbf{u}^2

for the energy flux per unit area, E, where (k,h) are the wavenumber and amplitude that characterize the bathymetry, and \mathbf{u} is the barotropic tidal velocity vector. The topographic roughness, h^2, is derived from the high-resolution bathymetry ETOPO2v2 (Smith and Sandwell 1997; NGDC 2006)) as the root-mean-square of the topography over a 50km smoothing radius, and k is a free parameter set as k = 2\pi/125 \mathrm{km}^{-1}. It should be emphasized that internal wave conversion equation is a scale relation, and not a precise specification of internal tide energy-flux. In the barotropic tidal model, the value of k was tuned to give the best fit to the observed tides.

9. Other Model Physics

9.1. Equation of State

POP requires an equation of state in the form

(223)\rho = \rho(\Theta,S,p),

to relate density, \rho, to the prognostic variables \Theta and S; see (10). To avoid a nonlinear integration of the hydrostatic equation (9), the pressure, p, in the equation of state is approximately evaluated as a time-independent function of depth p_0(z) by means of the equation

(224)\begin{aligned}
p_{_0}(z) = 0.059808[\mathrm{exp}(-0.025z)-1]
+0.100766z + 2.28405\times 10^{-7}z^2,
\end{aligned}

where pressure p is in bars and depth z is in meters. This formula is derived in Dukowicz (2001) from the latest Levitus global mean climatology. At a number of places in the code partial derivatives of \rho with respect to \Theta and S may also be necessary.

There are four options in the code for the equation of state:

(A) The density may be obtained from the formula of Jackett and McDougall (1995).

(225)\rho = \frac{\rho(\Theta,S,0)}{1 - p/K(\Theta,S,p)},

The numerator, \rho(\Theta,S,0), is a 15-term equation in powers of S and \Theta with coefficients given by the UNESCO international standard equation of state (Fofonoff and Millard 1983). The secant bulk modulus, K(\Theta,S,p), is a 26-term equation in powers of \Theta, S and p, obtained from a least-squares fit to the expression for the secant bulk modulus as a function of the in situ temperature (Fofonoff and Millard 1983). This density equation is valid in the range 0 \le S \le 42 psu, -2^\circ \le \Theta \le 40^{\circ}C, and 0 \le p \le 1000 bar. The total of 41 terms in this equation of state makes it the most expensive of the four available options.

(B) A cubic polynomial fit to the UNESCO international standard equation of state (Fofonoff and Millard 1983) which has the form

(226)\rho(\Theta,S,z_k) = \rho_k(\Theta-\Theta_{ref}(k),S-S_{ref}(k)),

where \Theta_{ref}(k) and S_{ref}(k) are reference values of the potential temperature and salinity for each model level. The set of nine coefficients in the cubic polynomial for each level are pre-calculated based on a least-squares fit to the international standard equation of state as prescribed by Bryan and Cox (1972). The range of validity of this equation is depth dependent and is specified during the fitting procedure to encompass typically observed oceanographic conditions. This results in a fairly narrow range of validity at depth which may not be appropriate for certain simulations or for evaluating density resulting from large adiabatic displacements. The reduction from 41 terms to 9, however, represents a considerable cost savings over the full UNESCO equation of state.

(C) A 25-coefficient approximation of the Feistel and Hagen (1995) equation of state (which is more accurate relative to laboratory data than the UNESCO form), derived by McDougall et al. (2003).

(227)\rho(\Theta,S,p) = \frac{P_1(\Theta,S,p)}{P_2(\Theta,S,p)},

where P_1 is a 12 term polynomial and P_2 is a 13-term polynomial. This density equation is valid in the range 0 \le S \le 40 psu, -2^\circ \le \Theta \le 33^{\circ}C, at the surface, diminishing to 30 \le S \le 40 psu, -2^\circ \le \Theta \le 12^{\circ}C at 550 bar. However, the authors report that the equation is well behaved in the range 0 \le S \le 50 psu, -10^\circ \le \Theta \le 50^{\circ}C, and 0 \le p \le 1000 bar This equation of state is intermediate in cost between the UNESCO form (A) and the cubic polynomial form (B). It is currently the preferred option for CCSM integrations. An updated form of this equation of state Jackett et al. (2006) has been implemented in a future version of POP and should be available soon.

  1. A simple linear equation of state given by

(228)\rho=\rho_{_0} - 2.5\time 10^{-4} \Theta + 7.6\times 10^{-4} S\;,

where \rho is in gm/cm^3, \Theta is in ^\circC, and S is in practical salinity units (psu). The \rho_{_0} term is not included in the code for the full density with the linear EOS.

9.1.1. Boussinesq Correction

The pressure gradient terms in the momentum equations (44) and (45) are an approximation of \rho^{-1}\nabla p, introduced as part of the Boussinesq approximation. Furthermore, the conversion of pressure to depth in the equation of state by means of (224) can have significant dynamic consequences through its effect on the pressure gradient (Dewar et al. 1998). Both of these errors can be greatly reduced by transforming the density as follows (Dukowicz 2001):

(229)\rho = r(p)\rho^{*},

where \rho^{*} is termed the thermobaric density, and r(p) is a nondimensional function of pressure that extracts the pressure-dependent part of the adiabatic compressibility from the density along the global-mean Levitus climatology:

(230)\begin{aligned}
r(p) = 1.02819 - 2.93161\times 10^{-4}
\mathrm{exp}(-0.05p) + 4.4004\times 10^{-5}p\;,
\end{aligned}

where p is in bars. This leads to the definition of an associated thermobaric pressure p^*:

(231)p^*(p) = \int^p_0\frac{dp'}{r(p')},

such that the hydrostatic equation (9) becomes

(232)\frac{\partial p^*}{\partial z} = -\rho^* g,

and the pressure gradient force is transformed into

(233)\frac{1}{\rho}\nabla p = \frac{1}{\rho^*}\nabla p^*
\approx \frac{1}{\rho_{_0}}\nabla p^*.

The effective equation of state in terms of these new variables becomes

(234)\rho^* = \rho(\Theta,S,p(p^*))/r(p(p^*)) = \rho^*(\Theta,S,p^*).

The advantage of this transformation is that the effective adiabatic compressibility associated with the equation of state is now at least an order of magnitude lower. This means that the variation of density with depth is much lower and that fluctuations of pressure in the equation of state have a much smaller effect on density. Thus, the errors associated with the Boussinesq approximation in the pressure gradient force and the linearization of the equation of state in the transformed variables are at least an order of magnitude smaller than without the transformation, as explained in Dukowicz (2001).

The pressure gradient force is approximated by \rho_{_0}^{-1}\nabla p^*, and the hydrostatic equation becomes

(235)\frac{\partial p^*}{\partial z} =
-\frac{\rho(\Theta,S,p_{_0}(z))g}
{r(p_{_0}(z))}.

This implies that the pressure variable handled in POP is the thermobaric pressure p^*, not p. The pressure enters only in the pressure gradient and in the equation of state and effects nothing else. Should the true pressure ever be required it is easily obtainable from p^* using the relation (231). Given this, we shall drop the notation p^* and henceforth interpret p to imply the thermobaric pressure.

9.1.2. Hydrostatic Pressure

The pressure at depth z is obtained by integrating the hydrostatic equation from z=0 (the “hydrostatic pressure” p_h) and adding the contribution from the surface pressure p_s associated with undulations of the free surface:

(236)\begin{aligned}
p(x,y,z) &=& p_s(x,y) + p_h(x,y,z) \nonumber\\
p_h(x,y,z) &=& \int^{0}_{z} dz' g\rho(x,y,z')
\end{aligned}

In the code only the horizontal pressure gradients are needed, and these are evaluated by integrating the horizontal gradients of density in (236). The discrete formulas for the horizontal pressure gradients at level k are:

(237)\begin{aligned}
&&\delta_x\overline{p_k}^y = \delta_x\overline{p_s}^y
+ g\sum_{m=1}^{k}\frac{1}{2}[\delta_x\overline{\rho^{*}_{m-1}}^y
+\delta_x\overline{\rho^{*}_m}^y]dz_{m-\frac{1}{2}}
\nonumber\\
&&\delta_y\overline{p_k}^x = \delta_y\overline{p_s}^x
+ g\sum_{m=1}^{k}\frac{1}{2}[\delta_y\overline{\rho^{*}_{m-1}}^x
+\delta_y\overline{\rho^{*}_m}^x]dz_{m-\frac{1}{2}}
\end{aligned}

where dz_{m-\frac{1}{2}} = dzw_{m-1}, and dz_{-\frac{1}{2}}=0.5dz_{1}. Here \rho^{*}_m = \rho_m / r(p_{_0}zt_m) if the Boussinesq correction is used, otherwise \rho^{*}_m =\rho_m, where \rho_m is the density at level m, and \rho_{m}=\rho_1 when m=0.

9.1.3. Expansion Diagnostic

This is a feature that is currently not in the code but will be added later. The expansion diagnostic is supposed to account for the change in global mean sea level due to two effects: 1) the change in volume associated with the net accumulated freshwater flux into the ocean (from precipitation, evaporation, melting and freezing of sea ice, and river runoff); and 2) the change in volume associated with steric expansion due to changes in global mean density. If the model employs virtual salinity fluxes, either from restoring to climatological surface salinity or from the conversion of actual freshwater fluxes to virtual salinity fluxes (e.g., when the model is run without using the natural boundary conditions for freshwater flux as discussed in Sec. Variable-Thickness Surface Layer) then the expansion effects due to 1) and 2) are not easily separated: a virtual salinity flux will change the salinity and hence the density through the equation of state, on the other hand, a virtual salinity flux is usually associated with an actual freshwater flux which changes the volume as well, and it is not clear how to cleanly separate these volume and density changes. However, in the steric effect the change in density is primarily associated with thermal expansion due to heating or cooling, rather than changes in salinity. Hence, we expect that computing the steric effect due to density changes without including volume changes due to freshwater flux provides a useful diagnostic. Therefore, in the code the change in surface elevation due to the steric effect can be computed assuming there is no change in total volume. The change in surface elevation due to steric expansion is computed as:

(238)\begin{aligned}
\langle\eta^n\rangle = \langle H\rangle \left( \frac{\langle\rho^{_0}\rangle}
{\langle\rho^{n}\rangle} - 1\right)
\end{aligned}

where \langle H\rangle, the ratio of ocean volume to surface area, is the mean ocean depth, \langle\rho^{_0}\rangle is the initial global mean density, \langle\rho^{n}\rangle is the global mean density at the n^{\mathrm{th}} model timestep computed assuming there is no change in total volume, and \langle\eta^n\rangle is the estimated change in sea-surface elevation at the n^{\mathrm{th}} timestep due to steric expansion only. If natural boundary conditions for freshwater flux are employed, then the actual volume of the ocean in the model does change (see Sec. Variable-Thickness Surface Layer), and the corresponding change in mean sea level is diagnosed as the actual global-mean sea level in the model. This change can be added to the change due to steric expansion to obtain an approximate total change in global mean sea level due to the combined effects of 1) and 2).

9.2. Sea-Ice Formation and Melting

Air-sea heat fluxes, when applied over a time interval, can produce regions of subfreezing model upper-layer temperatures. This situation can be alleviated by adjusting the model \Theta and S (and the first layer thickness, if the variable surface thickness formulation is used with freshwater fluxes) that is associated with frazil sea-ice formation. The choice of when to perform such adjustments is judiciously made to ensure tracer conservation with minimum adjustments. So, in coupled simulations with averaging as the time mixing option, such ice formation time steps are performed at both the coupling time step and the time step just prior to it so that both branches of the leap-frog time steps are brought to the freezing temperature which will then be used over the next coupling interval. If Matsuno time mixing option is selected, checking for ice formation only at the coupling time steps is adequate because the first time step of the next coupling interval is Matsuno and the unadjusted branch of the leap-frog time steps will be ignored. The following algorithm describes the procedure for the virtual salt flux formulation, default for the present integrations.

After \Theta and S are updated, the upper kmxice model layers (denoted by index k) are scanned for \Theta below the freezing point, \Theta_f. Here, kmxice is a model namelist input in ice_nml and, at present, \Theta_f is not a function of local salinity S_k, but \Theta_f= - 1.8^\circ C. Although we describe the most general case here, in practice we do recommend setting kmxice=1, because unless a monotone advection scheme is in place, model advective errors could lead to \Theta far below freezing in some deep locations and the scheme will form too much ice.

At each layer, the potential mass per unit area of ice melt ({\tt POTICE}_k < 0) or ice formation ({\tt POTICE}_k > 0) is computed as

(239){\tt POTICE}_k = \max \left[ \frac{\rho_{sw} c_p}{L_f}
 \Delta z_k (\Theta_f - \Theta_k), {\tt QICE}_{k+1} \right],

where \rho_{sw} and c_p are the density and heat capacity of sea water, respectively, L_f is the latent heat of fusion, \Delta z_k is the layer thickness, and \Theta_k is the local potential temperature. Any ice that forms at depth is assumed to float towards the surface and this ice flux, {\tt QICE}_k (defined positive downwards so {\tt QICE}_k \leq 0), is accumulated bottom to top as

(240){\tt QICE}_k = \sum_{kmxice}^k -{\tt POTICE}_k ,

assuming no ice formation below the kmxice layer, i.e. {\tt QICE}_{k=kmxice+1}=0. As ice floats to the surface, it can either partially or completely melt in upper layers whose temperatures are above freezing.

At each layer, \Theta and S are both adjusted in accord with the ice formed or melted in the layer:

(241)\begin{aligned}

\Theta_k &= \Theta_k +
            \frac{L_f}{\Delta z_k \rho_{sw} c_p} {\tt POTICE}_k, \\
S_k      &= S_k + \frac{(S_o - S_i)}{\rho_{sw} \Delta z_k}{\tt POTICE}_k.
\end{aligned}

For S, this is equivalent to replacing a volume of formed ice at salinity S_i with an identical volume of water at salinity S_o. Here S_i and S_o represent sea ice and ocean reference salinities, respectively, and they are set to constant values in order to ensure salt conservation. In addition, if POP is coupled to a sea-ice model, both must use the same value for S_i to ensure conservation. Note that in CCSM2, S_i=0, but in CCSM3 S_i=4 psu.

The ice flux is accumulated at each location over the number of ice formation time steps, N, during a coupling interval as

(242){\tt AQICE} = \sum_1^N weight_1 \ {\tt QICE}_{k=1},

where weight_1 is either 1/2 or 1 depending on whether it is an averaging time step or not.

If, at the end of a vertical scan, the surface temperature \Theta_1 remains greater than freezing, ({\tt QICE}_{k=1} =0), then the excess heat melts previously formed ice, and AQICE, \Theta_1, and S_1 are adjusted accordingly. Thus, over a coupling interval, ice is assumed to remain where it was formed. In coupled simulations, the accumulated ice during a coupling interval is passed to an ice model via the flux coupler as an equivalent downward heat flux, QFLUX in W m^{-2}:

(243){\tt QFLUX} = \frac{-L_f}{\Delta t^*} \ weight_2 \ {\tt AQICE}

where AQICE includes any adjustment due to melting previously formed ice over the last timestep, \Delta t^* = 1 day is the coupling interval, and weight_2 is either 1/2 or 1 depending on averaging or Matsuno time steps, respectively. In the last two equations, the presence of weights assures that \Theta and S budgets will be conserved when the averaging option is selected where, after averaging, actual \Theta and S changes become 1/2 of what is implied by the ice formation fluxes.

Thus, \Theta and S adjustments for the total ice formed ({\tt QFLUX} > 0) are made prior to the ice being passed to the sea-ice model. However, warm surface model temperatures result in {\tt QFLUX} < 0, which is a potential to melt ice in the sea-ice model. Since the ocean does not know if sufficient ice is present at a given location, \Theta and S adjustments are delayed until the appropriate heat and freshwater fluxes are received back from the sea-ice model through the coupler. If sufficient ice is present in the absence of all other fluxes, the heat flux will be equal to the cooling needed to make \Theta_1 = \Theta_f, when applied over the next coupling interval.

The ice formation option can also be activated in ocean-only (uncoupled) simulations subject to observational or reanalysis based air-sea heat fluxes. In this case, the accumulated ice is used internally for local ice formation and melt, and AQICE is saved in restart files for exact continuation. If the Matsuno time step mixing option is chosen, the ice time steps are set to be the time steps just before a mixing step by default. For the averaging time step option, the default is to form ice every time step.

The POP ice formation subroutine contains a modified algorithm when the freshwater flux formulation is chosen. Here, the volume of ice formed at layer k is replaced with an equal volume of water (because the thickness of the below surface layers cannot change) from the layer above with salinity S_{k-1}. The surface layer S can change both through exchanges of S with lower layers and change of thickness due to formed ice sent to the ice model. Because there still remain some consistencies and scientific issues with ice formation and melt when freshwater fluxes are used, this algorithm has not yet been tested, and is not recommended to be used at the present time.

9.3. Freshwater Balancing over Marginal Seas

Unlike the surface heat flux and sea surface temperature, there are no appreciable feedbacks between the surface freshwater flux and S. This is particularly so in isolated marginal sea regions where the freshwater fluxes can produce unphysical S values throughout the water column. The situation can be especially severe in coupled integrations when the marginal sea regions receive river runoff fluxes. In uncoupled integrations, one obvious remedy is to use restoring freshwater or salt fluxes instead in these regions where the surface S is restored to some climatological distribution with a relatively small restoring time scale. In this section, we describe a more novel approach, particularly designed for coupled integrations. Here, the amount of excess or deficit of freshwater flux over a marginal sea is transported to or from its associated active-ocean region, thus implicitly connecting marginal seas with active ocean and providing a means for marginal sea runoff to eventually discharge into the open ocean. This process assures that the volume-mean S stays constant throughout the integration within each marginal sea, eliminating any unphysical values in S.

When balancing is requested, the active-ocean regions corresponding to each marginal sea, ms, are determined at the beginning of each integration. For this computation, a longitude, a latitude (both in degrees), and a distribution active-ocean area size (in cm^2) are provided for each region in an input file. Using the longitude and latitude values as starting points, a search window is created to find the active-ocean points associated with each ms. This iterative process continues until the cumulative active-ocean grid areas match the input distribution size. At each iteration, the search window is increased by 1^\circ on all four sides. If the hard-coded values for either the maximum number of iterations or the maximum number of active-ocean grid points per ms are exceeded, POP will stop with an error message.

The final product of this initialization procedure is a global array of area fractions for each ms, F_{i,j}^{ms}. In F, the fractions are computed simply as the ratios of grid areas to the total computed distribution area. It is ensured that the sum of these area fractions for each ms is unity. The distribution regions for different ms can coincide. Note that, the zero elements of F represent the ocean points outside of the designated regions.

In fully coupled integrations in which the surface fluxes stay constant until the next coupling, balancing is performed once per coupling interval, immediately following the arrival of the new fluxes. If the fluxes do change during a coupling interval (e.g. partially-coupled option), then balancing is done every time step. For each ms, a transport term, T, is evaluated in Kg s^{-1} of freshwater:

(244)\begin{aligned}

T^{ms} = \sum_{ms} \bigl( \bigl[ & \max(0,{\tt QFLUX}_{i,j})c_q + \\
  & (E_{i,j} + P_{i,j} + R_{i,j} + M_{i,j}) + F_{i,j}^S c_s \bigr]
  \Delta x_{i,j} \Delta y_{i,j} \bigr), \nonumber\end{aligned}

where only the fluxes over the marginal seas contribute to the sum which is performed over each ms individually. In the above equation, {\tt QFLUX} > 0 is the frazil ice formation in W m^{-2} and E, P, R, and M represent the evaporation, precipitation, river runoff, and ice melt freshwater fluxes from the coupler in Kg m^{-2} s^{-1}. F^S is the salt flux due to ice melt in Kg of salt m^{-2} s^{-1}. F^S was zero in the CCSM2 because S_i=0, but is nonzero in the CCSM3 because S_i=4 psu in the later version. Finally, \Delta x and \Delta y are the zonal and meridional grid spacings (in m) centered at \Theta points, respectively, and c_q and c_s represent unit conversion factors given by

(245)c_q = - \frac{1}{L_f} \frac{S_o - S_i}{S_o}, \quad
c_s = - \frac{1}{S_o} \frac{\rho_{fw}}{\rho_{sw}}.

T^{ms} represents the excess (or deficit) freshwater flux for marginal sea ms that needs to be transported to or from its associated active-ocean region. How much of T^{ms} is transported to or from an active-ocean grid point is determined by

(246)\mathrm{MSTF}_{i,j} = \frac{ F_{i,j}^{ms} T^{ms} }
                           { \Delta x_{i,j} \Delta y_{i,j}}.

At these points, the surface freshwater or salt flux is modified as

(247)SF_{i,j}^{new} = SF_{i,j}^{old} + c_t \mathrm{MSTF}_{i,j},

where c_t=-S_o/\rho_{fw} or c_t = 1 depending on if POP is forced with virtual salt or freshwater flux boundary conditions, respectively.

In marginal seas where there is no frazil ice formation, the surface freshwater flux is set to zero. Otherwise, the surface flux is simply computed as -\max(0, {\tt QFLUX}_{i,j}) c_q  c_t to undo the S adjustment that has been done in the previous time steps.

Because this algorithm has not yet been used with the freshwater flux formulation, its use is not recommended with that formulation.

9.4. Overflow Parameterization

A new overflow parameterization (OFP) for deep channel overflows (e.g., Denmark Strait and Faroe Bank Channel) and continental shelf overflows (e.g., Weddell and Ross Sea) has been developed and implemented in POP. Although this OFP is based on the Price and Yang (1998) Marginal Sea Boundary Condition (MSBC) scheme, there are substantial differences between the two. The OFP consists of two parts: i) calculation of the overflow properties based on the evolving ocean model state, and ii) modifications of the model boundary conditions, equations, and bottom topography to incorporate the overflow volume and tracer fluxes as well as to ensure volume and tracer conservations. The first part uses the same exchange and entrainment formulas as in the MSBC, but all the necessary fields are obtained prognostically. Furthermore, the marginal seas providing the overflow source waters (e.g., the Nordic Sea) are a part of the prognostic model domain rather than just some marginal sea boundary conditions as in the MSBC and the inflow into these marginal seas is accomplished by the resolved flow in contrast with a parameterized inflow in the MSBC. The second part of the OFP is entirely new, particularly in its treatment of the baroclinic and barotropic momentum and continuity equations to conserve volume.

In POP, overflows can be activated in either interactive or diagnostic mode, where the latter computes the overflow without affecting the model solution. The OFP is used to represent the Denmark Strait, Faroe Bank Channel, Ross Sea, and Weddell Sea overflows. When the overflow module is used either in active or passive mode, the model requires an input file to specify all the details of each overflow region (Briegleb, Danabasoglu, and Large 2010). We note that the OFP is not used to represent the Mediterranean overflow through the Strait of Gibraltar. Instead, we use an open Gibraltar Strait and a sufficiently deep model topography on the Atlantic side of the strait, so that the overflow can convect as deeply as needed locally. This configuration avoids the excess entrainment associated with staircase topography (Wu, Danabasoglu, and Large 2007).

The OFP and its implementation in POP are detailed in Briegleb, Danabasoglu, and Large (2010) and Danabasoglu, Large, and Briegleb (2010), so only a general overview is described here. A schematic of the Nordic Sea overflows and its parameterization are shown in Figure8.1  to help with the following description. For the same purpose, Figure8.2  gives the bottom topography of the model in the vicinity of the Nordic Sea overflows, including the lateral extents of the tracer averaging regions and the pre-specified locations of the possible injection sites for the product waters.

../_images/FigOverflow.jpg

Figure 8.1: A schematic of the Nordic Sea overflows. T, S, \rho, and M represent potential temperature, salinity, density, and volume transport, respectively. The subscripts i, s, e, and p indicate inflow, source, entrainment, and product water properties, respectively. d is depth, x_{ssb} is the distance from the sill to the shelf break, \alpha is the maximum slope of the continental shelf near the shelf break, C_d is the bottom drag coefficient, W_s is the width of the channel at the sill depth, and h_u is the upstream thickness of the source water. The question marks denote the prespecified product water injection locations. The thick, short arrows indicate flow directions. The green box shows the raised bottom topography. The other boxes represent the regions whose T and S are used to compute the necessary densities.

../_images/FigNAtlTopo.jpg

Figure 8.2: Bottom topography as represented in the model in the vicinity of the Denmark Strait (DMS) and Faroe Bank Channel (FBC) overflows. The colors indicate the model vertical levels. The corresponding depths are given above the color bar. The boxed regions denoted by I, E, and S indicate the Inflow, Entrainment (thin box), and Source regions in the horizontal, respectively, whose potential temperature and salinity properties are used to compute the necessary densities. The source and entrainment box edges at which the respective water properties and transports are imposed as side boundary conditions in the model are indicated by the black arrows, showing directions corresponding to flows out of the model domain. The white lines denoted by P show the prespecified product water injection locations into the model domain. All product water sites have the same injection direction as denoted by the white arrows drawn at only a few of the sites for clarity.

9.4.1. Calculation of Overflow Properties

We use the subscripts i, s, e, and p to denote inflow, source, entrainment, and product water properties, respectively. In addition to the equations given below, the overflow dynamics depend on the geography of the straits and channels under consideration. Here, \varphi is the latitude, h_u is the upstream source thickness, W_s is the width of the strait, d_s and d_e are the sill depth and the depth of the shelf break where entrainment occurs, respectively, x_{ssb} is the distance from the channel exit to the shelf break, \alpha is the maximum bottom slope near the shelf break, and C_d is the bottom drag coefficient.

The overflows are driven by the density difference between the source density, \rho_s, and the open ocean inflow density, \rho_i, as expressed by the source reduced gravity

(248)g_s^\prime = \frac{\rho_s - \rho_i}{\rho_0} \  g

where g is the gravitational acceleration and \rho_0 = 1027 kg m^{-3} is a reference density. \rho_i and \rho_s are evaluated at the sill depth using

(249)\begin{aligned}

\rho_{i} &= \rho (T_i, S_i, d_s), \nonumber \\
\rho_{s} &= \rho (T_s, S_s, d_s). \end{aligned}

Here, T and S represent volume-average potential temperature and salinity, respectively, calculated within the corresponding horizontal regions shown in Figure8.2 . There is no averaging in the vertical as only one vertical level is used.

As long as g_s^\prime > 0 source overflow transport will occur. Assuming that the strait or channel width is greater than the radius of deformation and the inflow is not geometrically constrained, following Whitehead, Leetmaa, and Knox (1974), this source overflow transport, M_s, is obtained using the expression for rotating, hydraulically controlled maximum geostropic flow through a strait

(250)M_s = \frac{g_s^\prime h_u^2}{2f}.

It is important to note that (250) uses h_u (< d_s) which is different than the source thickness above the sill as the source waters exit the strait, denoted as h_s. Again following Whitehead, Leetmaa, and Knox (1974), h_s is calculated from h_u as

(251)h_s = \frac{2}{3}h_u.

Assuming a rectangular cross sectional area with height h_s and width W_s, i.e., A_s=h_s W_s, at the strait exit, an associated source speed, U_s, can be evaluated

(252)U_s = \frac{M_s}{A_s}.

As the source water flows down along the continental slope, it accelerates, spreads in width, and therefore thins. It also entrains ambient waters in the open ocean, largely near the shelf-slope break. These processes are parameterized using the end-point model of an entraining, rotating density current developed by Price and Baringer (1994). The entrainment is driven by the entrainment reduced gravity given by

(253)g_e^\prime = \frac{\rho_s^\prime - \rho_e}{\rho_0} \ g,

where

\begin{aligned}
\rho_{s}^\prime &= \rho (T_s, S_s, d_e), \nonumber \\
\rho_{e}        &= \rho (T_e, S_e, d_e),\end{aligned}

with \rho_s^\prime and \rho_e the source and entrainment region densities, respectively, both computed at the entrainment depth d_e, using volume-average T and S for the corresponding lateral regions shown in Figure8.2 . As long as g_e^\prime > 0 entrainment will occur.

At the shelf break, the flow is assumed to have a characteristic speed governed by geostrophic balance with the reduced gravity flow

(254)U_{geo} = \frac{g_e^\prime \alpha}{f}.

The average flow speed between the channel exit and the shelf break point is then given by

(255)U_{avg} = 0.5 \ (U_s + U_{geo}).

The spreading width is assumed to increase linearly with distance from the sill such that when the source water reaches the shelf break it has width W(x_{ssb}) given by

(256)W(x_{ssb}) = W_s + 2 K_{geo} x_{ssb}

with the Ekman number, K_{geo}, specified by the ratio of bottom drag to Coriolis force over downslope flow

(257)K_{geo} = \frac{C_d U_{avg}}{\frac{1}{2} f (h_s + h_{geo})}.

During the lateral spread, the overflow thickness decreases by volume conservation. At the shelf break, the thickness is given by

(258)h_{geo} = \frac{U_s h_s W_s}{U_{geo} W(x_{ssb})}.

Using (257), (256) and (258) can be solved simultaneously for h_{geo}.

A geostrophic Froude number, F_{geo}, for the entrainment mixing at the shelf-break is defined as

(259)F_{geo} = U_{geo} / \sqrt{g_e^\prime h_{geo}}

from which an entrainment parameter, \vartheta, representing the ratio of entrained to product water volume transports, can be evaluated as

(260)\vartheta = \frac{M_e}{M_p} = 1 - F_{geo}^{-2/3}.

If g_e^\prime \le 0 or F_{geo} \le 1, \vartheta is set to 0. Given M_s and \vartheta, the entrainment volume transport is

(261)M_e = M_s \frac{\vartheta}{1 - \vartheta},

and using volume conservation, the product volume transport is then

(262)M_p = M_s + M_e.

Finally, from tracer conservation, the product water potential temperature and salinity are calculated using

(263)\begin{aligned}

T_p &= T_s (1-\vartheta) + T_e \vartheta, \nonumber \\
S_p &= S_s (1-\vartheta) + S_e \vartheta.\end{aligned}

We note that any additional tracers that the model carries, e.g., the ideal age and CFCs, use the same averaging procedure and conservation relations as in T and S.

9.4.2. Model Modifications

As indicated by the green box in Figure8.1 , we raise the discrete model bottom topography by one vertical level at the sill depth for each parameterized overflow. This popped-up topography occurs at a minimum of three consecutive lateral grid points per overflow along the source water injection site and effectively reduces flow over the sill, replacing it with the parameterized source flow M_s. These raised vertical levels represent the only topographic modifications done when the overflow parameterization is activated compared to the original topography used in simulations without the OFP. The source flow velocity and the associated tracer fluxes are imposed as boundary conditions on the source marginal sea side of the raised level, i.e., at the sill depth level d_s, as denoted by the blue box in Figure8.1 . Similarly, the entrainment flow velocity and tracer fluxes are specified as side boundary conditions on the open ocean side (Atlantic Ocean in the figure) at the entrainment depth d_e as shown by the brown box in Figure8.1 . We note that both the source and entrainment flows are directed into the topography as indicated by the arrows.

We pre-specify a product water pathway for each overflow, consisting of N_P number of possible injection sites as shown by the white lines in Figure8.2 . To determine where to inject the product water along this path, we compare the product water density, \rho_p^n = \rho (T_p, S_p, d_p^n), with the open ocean ambient water density evaluated at the same depth d_p^n, using the volume-average T and S from the immediate downstream of the nth product water site. Assuming that the product water injection site depths monotonically increase with n, i.e., the N_Pth site representing the deepest one, our search starts with the (N_P-1)th site and ends at the first one. If \rho_p^n exceeds the ambient density at level n, the product water is inserted as a lateral boundary condition at the next deepest site n+1. This process ensures that the product water goes to the deepest possible site. The injection occurs at the shallowest site if \rho_p^n is never larger than the ambient density at any of the sites. As shown in Figure8.1 , the product water is directed out of the bottom topography. In effect, the flow between source to entrainment and then entrainment to product occurs below the model bottom topography through the parameterization. We remark that the ambient ocean density along the product path is not necessarily monotonically increasing, because the product path extends laterally as it deepens.

The specification of the overflow volume transports as side wall velocity boundary conditions for the source, entrainment, and product sites requires careful modifications of the baroclinic and barotropic equations to ensure local mass (volume) conservation from which the global conservation follows automatically. Here, we repeat some equations for clarity. With baroclinic and barotropic split, the total velocity is written as

(264){\bf u} = {\bf \tilde u}^\prime + {\bf U}

where {\bf \tilde u}^\prime and {\bf U} are the baroclinic and barotropic velocities, respectively. They are defined by

(265){\bf U} = \frac{1}{H + \eta} \int_{-H}^\eta {\bf u}~ dz, \quad
\frac{1}{H+\eta} \int_{-H}^\eta {\bf \tilde u}^\prime ~ dz =  0.

We assume that the injected side wall overflow velocities represent the total velocity {\bf u}. For those U-grid columns where these injections occur, H must be extended downwards to the bottom of the vertical grid level where the overflow velocities are actually specified. This new depth is given by

(266)H^\prime = H + \Sigma \Delta z

where \Delta z denotes the vertical level thicknesses and the summation is done for the side wall heights from H down to the bottom of the injection level. To avoid complications arising from moving product water injection locations, we change H to H^\prime at all N_P product water sites. If the injection level is several vertical grid levels below H, the thicknesses of the levels in between are included in this summation. These in-between levels are assumed to have {\bf u}=0. Thus, we replace H with H^\prime in all model equations where appropriate, including any relevant vertical integrals. For example, (265) becomes

(267){\bf U} = \frac{1}{H^\prime + \eta} \int_{-H^\prime}^\eta {\bf u}~dz, \quad
\frac{1}{H^\prime+\eta} \int_{-H^\prime}^\eta {\bf \tilde u}^\prime ~dz = 0.

As described earlier, POP first solves the momentum equations, without including the surface pressure gradient, for an auxiliary velocity {\bf u}^{\prime}, giving us {\bf u}^{\prime} between H \le z \le \eta. The relationship between this intermediate velocity and {\bf \tilde u}^\prime is

(268){\bf \tilde u}^\prime = {\bf u}^{\prime} -
      \frac{1}{H^\prime + \eta}\int_{-H^\prime}^\eta{\bf u}^{\prime} ~ dz
    = {\bf u}^{\prime} - {\overline {{\bf u}^{\prime}}}.

Our assumption of {\bf u}=0 at the side walls between H and any overflow injection levels leads to {\bf \tilde u}^\prime = -{\bf U} at these in-between levels. Similarly, because {\bf u} = {\bf u}_{OVF}, we have {\bf \tilde u}^\prime = {\bf u}_{OVF} - {\bf U} at all the injection locations. Here, {\bf u} = {\bf u}_{OVF} is a generic variable used to indicate the overflow velocities. So, knowing {\bf U} (see below), {\bf u}^{\prime} between -H \le z \le \eta, and {\bf \tilde u}^\prime between -H^\prime \le z < H, we use (267) and (268) to obtain an equation for {\overline {{\bf u}^{\prime}}} for the entire column

(269){\overline {{\bf u}^{\prime}}} =
\frac{1}{H} \left[ \int_{-H}^\eta {\bf u}^{\prime} ~ dz +
\int_{-H^\prime}^{-H} {\bf \tilde u}^\prime ~ dz \right].

We then evaluate {\bf \tilde u}^\prime using (268) for -H \le z \le \eta. This process results in local volume conservation, consistent with the imposed overflow velocities, for the T-grid columns to the immediate upstream (to the right in Figure8.1 ) of the source and to the immediate downstream (to the left in Figure8.1 ) of the entrainment and product injection locations.

These velocity modifications obviously impact the horizontal velocity divergences, thus destroying the local volume conservation, at the T-grid columns immediately downstream (to the left in Figure8.1 ) of the source and immediately upstream (to the right in Figure8.1. ) of the entrainment and product injection sites. For example, the T-grid column above the raised sill has a non-zero vertical velocity at the bottom of the column. The magnitudes of these imbalances are equal to the volume transports associated with the source, entrainment, and product waters for each affected column. Therefore, we modify the vertically integrated continuity equation for these columns to account for these transports and thus enforce continuity, viz.,

(270)\frac{\partial \eta}{\partial t} + \nabla \cdot (H^\prime + \eta) {\bf U}
-q_w - \delta_{OVF} \frac{M_{OVF}}{\Delta A} = 0

where \delta_{OVF} is 1 for the affected columns and 0 elsewhere, M_{OVF} represents either M_s, M_e, or M_p, and \Delta A is the appropriate surface area of the affected T-grid columns. (270) together with the vertically integrated momentum equation are then used to compute both \eta and {\bf U}. It is important to reiterate here that the last term on the left-hand side of (270) accounts for the changes in the horizontal velocity divergences within a specified T-grid column and it should not be viewed as mass (volume) injections / extractions from either the surface or the ocean bottom.

The changes in the tracer equations are relatively straightforward compared to the momentum equations. Specifically, the advection algorithms are modified to incorporate the non-homogeneous flux boundary conditions at the injection sites. By construction, this process conserves all tracers. For further details, see Briegleb, Danabasoglu, and Large (2010).

9.5. Passive Tracers

Passive tracers are defined as tracers that do not impact the physical solution, but are carried along with the flow with perhaps some internal transformations that do not feed back on the prognostic tracer. POP contains a tracer infrastructure that accommodates such tracers. Two examples are described here that provide diagnostics for the age of ocean water and long term transport.

9.5.1. Ideal Age

Ideal age is a tracer used to estimate ventilation timescales in the ocean (Thiele and Sarmiento 1990; England 1995). Ideal age satisfies the same tracer transport equation as potential temperature and salinity, e.g., equations (11) and (34), but with different surface boundary conditions. In particular, ideal age is set to zero in the model surface layer at each time step. In startup runs, ideal age is initialized to zero everywhere.

Each water parcel ages by one year for each year since last contact with the surface ocean. In practice, ideal age is incremented by the model time step at each time step. The only way that indefinite growth of ideal age is limited is by contact with the surface, or by mixing with waters that have recently contacted the surface. The oldest ages will be found in waters which have been removed from the surface for the longest period of time. Ideal age equilibrates to its steady state solution on a timescale of thousands of years. The ocean model carries ideal age as a default tracer, and it is reported in units of years.

9.5.2. Chlorofluorocarbons

Chlorofluorocarbons (CFCs) are non-toxic chemicals containing chlorine, carbon, and fluorine atoms. They were developed in the early 1930s and widely used as refrigerants, in aerosol applications, and as solvents. Their concentration increased in the atmosphere until it was recognized that they contribute to ozone depletion in the stratosphere, and they began to be phased out by the Montreal Protocol in the late twentieth century. CFCs enter the ocean via air-sea gas exchange with the atmosphere. Both CFC-11 and CFC-12 are inert in seawater, and have been widely measured throughout the world ocean in the late 20th century and early 21st century. CFCs provide a useful diagnostic for timescales and pathways of transport in ocean models (Dutay et al. 2002) because they can be directly compared against observations.

In POP, surface CFC fluxes are calculated from 1931.5 through 2008.5, largely following OCMIP 2 protocols (OCMIP 2000). However, instead of using the atmospheric fields specified by OCMIP, atmospheric datasets for 10m winds, atmospheric pressure and ice fraction are used in computing the air-sea CFC flux. Values specified for the atmospheric CFC partial pressure (pCFC) history up to 1998.5 are from OCMIP 2. Values from 1999.5 on are from Walker, Weiss, and Salameh (2009). Values were extended from 1999.5 to 2008.5 by converting from SIO-1998 scale to SIO-1993 scale in order to be consistent with pre-1998 values. Both CFC-11 and CFC-12 are reported in units of \mathrm{fmol}/\mathrm{cm}^3 (which is within a few percent of pmol/kg).

10. Forcing

Most of the surface forcing is implemented as a surface flux and applied as the flux across the top face in vertical mixing routines. The remaining details pertain to which fields are required and the options are covered in the POP User’s Guide. Here we describe the solar radiative forcing and some conservation details that require more description.

10.1. Penetration of Solar Radiation

POP contains several options for allowing solar radiation to penetrate the water column.

10.1.1. Jerlov Water Type

One form of this absorption allows the user to specify a Jerlov water type. This assumption for solar absorption has the form of a double exponential, with the fraction of surface shortwave solar radiation, SWFRAC, that reaches a depth, -z, given by :

(271)SWFRAC = coef \exp(\frac{z}{depth1}) + (1-coef)\exp(\frac{z}{depth2}).

The first term on the right-hand side represents the rapid absorption of the longer wavelengths (reds), so the extinction depth1 is small (0.35 \rightarrow 1.4m) and a fraction (0.58 \leq coef \leq 0.78) of the radiation is in this band. The extinction depth2 for the shorter wavelengths (blues) is much greater (7.9m \rightarrow 23m). These parameters are set according to a specific Jerlov water type, and are varied neither spatially, nor in time. An important restriction is that no solar radiation passes through the bottom of the model ocean, lest it be lost to the system. Also, in order to avoid numerical underflows associated with vanishingly small exponentials, SWFRAC=0 for depths deeper than 200m. Jerlov water type Ib is typical.

10.1.2. Absorption Based on Chlorophyll

Another option in POP activates a chlorophyll dependent solar absorption that is based on the parameterization of Ohlmann (2003). Its form is a double exponential, with the transmission from surface to depth, -z, given by:

(272)Tr =  A_1 \exp(B_1 z) + A_2 \exp(B_2 z)

When this transmission is multiplied by the net shortwave absorption of the entire column; i.e. down shortwave at surface times one minus ocean albedo, there results the net flux transmitted to depth, -z, for a given chlorophyll amount. The A,B coefficients are chlorophyll dependent. The effects of solar zenith angle and cloud amount are small, and are ignored in the present implementation. Chlorophyll is assumed to be uniformly distributed in the column. As with the Jerlov-based transmission, no solar radiation passes through the bottom of the model ocean. Chlorophyll amounts range from 0.01~mg~m^{-3} for pure ocean water, to 10~mg~m^{-3} in chlorophyll blooms. Transmissions are estimated to be accurate to 0.01.

For this option, the user must specify the chlorophyll amount in an input data file. This file contains twelve monthly chlorophyll amounts in mg~m^{-3} on the ocean horizontal grid. The monthly chlorophyll distributions were provided by Ohlmann, based on a limited set of SeaWIFS data spanning the period September 1997 to November 2001. Regional values range from 0.05~mg~m^{-3} typical of subtropical oceans far from continents, to mid-latitude, coastal and equatorial values from 0.3 to 0.6 mg~m^{-3}. A few coastal regions have bloom values of 10~mg~m^{-3}.

Transmissions vary spatially and are updated once a month, with a repeating annual cycle. Chlorophyll based transmissions are nearly the same as the standard Jerlov water type Ib transmissions for a chlorophyll amount of 0.2~mg~m^{-3}. Generally, the subtropical oceans away from continents have more transmissive surface layers than the Jerlov case, while midlatitude, coastal and equatorial surface layers are less transmissive.

10.1.3. Diurnal Cycle

Because POP is often run in situations where it receives only daily shortwave fluxes, there are three options for specifying how the net shortwave flux is distributed in time to mimic a diurnal cycle. The first option keeps the spatially varying flux constant for the duration of the coupling interval.

The second option distributes the flux across a 12 hour window centered at noon. Inside this window, the flux is proportional to a sine term. The flux is zero outside of this window. The distribution is normalized to conserve the integral of the shortwave flux over the coupling interval. This option, whose impact is described in Danabasoglu et al. (2006), works only when the coupling frequency is once per day. Note that the sine distribution curve is independent of longitude, latitude, and time of year.

In the third option, the flux is proportional to the cosine of the solar zenith angle. The distribution is normalized to conserve the integral of the shortwave flux over the coupling interval. This option can be used with any coupling frequency. The solar zenith angle is determined from longitude, latitude, time of year, and the solar declination angle, which is determined from orbital parameters passed from the flux coupler. The zenith angle and declination angle computations make use of CCSM shared code.

10.2. Ice Runoff

In CCSM4, the ocean model receives an ice runoff field, I, in addition to the liquid water runoff, R. To conserve heat in the coupled system, the associated phase change is accounted for in POP by adding the latent heat, -L_fI, to the total ocean surface heat flux. The freshwater flux due to total runoff is then R+I.

10.3. Tracer Fluxes and Averaging Time Steps

In coupled applications it is essential that an accurate accounting of tracer fluxes through the sea surface be maintained in order to allow for conservative exchange within the full coupled system. The use of a leapfrog scheme, with its even and odd branches, presents some complexity in this respect, and the use of averaging steps to suppress the leapfrog computational mode (Section Filtering Timesteps) must also be taken into account.

We consider our accounting problem over a coupling interval, defined by successive times of information exchange between ocean model and coupler. It is sufficient to consider the use of a single averaging step within that interval, but with the restriction that none of the three time levels which are to be averaged fall outside of the coupling interval. We note that this implies that the last time step prior to a coupling time step is always a leapfrog time step.

../_images/FigAvgSteps.jpg

Figure 9.1: Illustration of a coupling interval containing a single averaging step, with cpl identifying times at which information is exchanged with the coupler (time advances to the right). Steps taken before averaging are represented with dashed lines, steps taken after averaging are represented with solid lines.

The illustration presented in Figure9.1  includes five leapfrog steps within the coupling interval from time n to time (n+4 \frac{1}{2}). To streamline notation here we consider only the influence of the surface flux F on the leapfrog evolution of \varphi, omitting the reference point n from the time index and the geometric factors which would multiply the surface flux. Considering only the tendency due to surface fluxes then the tracer evolution of the two leapfrog branches begins with

(273)\begin{aligned}

\varphi^{1} &= \varphi^{-1}  + 2 \Delta t F^{0}  \nonumber \\
\varphi^{2} &= \varphi^{0}  + 2 \Delta t F^{1}\end{aligned}

An averaging step is then applied to time levels 0, 1 and 2 in order to suppress the leapfrog computational mode (Section Filtering Timesteps):

(274)\begin{aligned}

\varphi^{\frac{1}{2}} &= ( \varphi^{1} + \varphi^{0}) / 2  \nonumber \\
\varphi^{\frac{3}{2}} &= ( \varphi^{2} + \varphi^{1}) / 2\end{aligned}

We resume ordinary leapfrog timestepping, still using information exchanged at the beginning of the coupling interval, until we have produced \varphi^{4 \frac{1}{2}}, the tracer at the end of the coupling interval:

(275)\begin{aligned}

\varphi^{2 \frac{1}{2}} &= \varphi^{  \frac{1}{2}}  +
                           2\Delta t F^{1 \frac{1}{2}}  \nonumber \\
\varphi^{3 \frac{1}{2}} &= \varphi^{1 \frac{1}{2}}  +
                           2\Delta t F^{2 \frac{1}{2}}  \nonumber \\
\varphi^{4 \frac{1}{2}} &= \varphi^{2 \frac{1}{2}}  +
                           2\Delta t F^{3 \frac{1}{2}}\end{aligned}

Summing (273) and (275), and making use of (274), the change in the average of the leapfrog branches over the coupling interval may be expressed in terms of the fluxes which span the coupling interval, as

(276)\frac{1}{2} (\varphi^{4 \frac{1}{2}} + \varphi^{3 \frac{1}{2}}) -
\frac{1}{2} (\varphi^{0} + \varphi^{-1}) =
\Delta t  ( F^{0} +\frac{1}{2}F^{1} + F^{1 \frac{1}{2}} +
                  F^{2 \frac{1}{2}} + F^{3 \frac{1}{2}} )

We accurately account for the fluxes passing through the ocean surface based on this expression, with the exception of the very first coupling interval of a run which starts with a forward step, as in initial startup from climatology (we permit this one-time departure from strict conservation).

Under current practice, at the time this document is being prepared, the only temporal variation which may occur within the fluxes spanning a coupling interval is due to the use of an idealized diurnal cycle (Section Penetration of Solar Radiation), and due to the latent heat associated with the formation or melting of ice. The idealized diurnal cycle is conservative and used only when coupling is daily. It does not change the amount of incoming SW, but only its temporal distribution. Ice formation is only considered at the time of coupling, at which point both branches of the leap-frog time steps are modified, as explained in Section Sea-Ice Formation and Melting, resulting in the addition of the QFLUX term of (243) to the right hand side of (276).

If another method of controlling the leapfrog computational mode is to be used then the accounting of fluxes must be reconsidered. With use of Matsuno steps an exact accounting is straightforward. An exact accounting of fluxes is more problematic with use of a Robert-Asselin time filter, as discussed in Leclair and Madec (2009).

11. Parallel Implementation

POP has been designed to run effectively on a wide variety of high performance computing architectures. It has also been designed for portability, so that no changes in the code itself are required to compile and run on a new architecture. This has been achieved by adhering to language standards and by abstracting communication functions to support multiple messaging or other communication paradigms. No preprocessing directives are used, except to maintain some differences between the CCSM version and the standalone POP, primarily due to coupling and CCSM shared code issues. In this chapter, details of the domain decomposition algorithms, the communication architecture and parallel I/O are described. Instructions on configuring POP inputs for performance are included in the POP User’s Guide.

11.1. Domain Decomposition

In order to compute efficiently on parallel architectures, work must be divided among each distributed computing node. POP, like many other climate or fluid dynamics models, achieves this parallelism through domain decomposition, dividing the gridded domain into pieces that each node can compute relatively independently. To increase this independence, ghost cells or halo regions are defined so that a node will have some information about state variables on neighboring nodes (see Figure10.1 ). This enables the calculation on each node to proceed without communicating with other nodes until these regions need to be updated. In POP, using a halo width of two grid cells will ensure that the entire baroclinic tendencies can be computed without any inter-node communication.

../_images/FigGhostCells.jpg

Figure 10.1: An illustration of a ghost cell or halo region, showing the halo for a particular block based on neighboring blocks.

Originally, POP decomposed the domain into single blocks per node. From POP 2.0 onward, POP has used a more flexible decomposition. In the new decomposition, the horizontal domain continues to be divided into blocks, but multiple blocks can be distributed to each node, as shown in Figure10.2a , Figure10.2b , Figure10.2c , Figure10.2d . This structure enables the optimization of block size for cache performance. Multiple blocks on a node also provides a mechanism for hybrid parallelism, in which a threaded parallelism (e.g. OpenMP) can be used to distribute work within a multi-processor shared memory node, while message passing (e.g. MPI) is still used to communicate between distributed nodes. Furthermore, different numbers of blocks can be assigned to each node to balance the computational load. Finally, this scheme permits the elimination of blocks which contain only land points. Such a structure provides additional optimization tuning flexibility. The use of small blocks can eliminating more land, provide better load balancing and better performance on cache-based microprocessors. However, because each block has a halo region, blocks that are too small have a high surface-to-volume ratio and communication and synchronization issues become dominant. Generally, block sizes of 20-40 grid points in each direction appear to provide a reasonable balance, though each computer is different.

../_images/FigDomainAll.jpg

Figure 10.2a: The POP subblocking decomposition, showing the full domain with shaded land cells

../_images/FigDomainCart.jpg

Figure 10.2b: The POP subblocking decomposition, showing a traditional single block Cartesian

../_images/FigDomainSubblock.jpg

Figure 10.2c: The POP subblocking decomposition, showing a subblock decomposition

../_images/FigDomainBalance.jpg

Figure 10.2d: The POP subblocking decomposition, showing a distribution of blocks in which land block have been eliminated and the computational load has been balanced.

POP supports three different schemes for distributing blocks among processors. All of them eliminate land blocks, but each represents options for balancing the number of blocks per node.

11.1.1. Cartesian Distribution

The first option for block distribution is a simple Cartesian distribution. In this case, the sublocks in a domain are distributed according to their logical position in the domain with a two-dimensional rectangular decomposition, similar to that shown in Figure10.2b , except that the single blocks in that figure would correspond to the set of subblocks associated with those large portions of the domain. No load balancing occurs and some nodes may be left with no work to do if all blocks are land.

11.1.2. Rake Distribution

A second distribution scheme uses a simple rake algorithm to balance the work per node. The algorithm starts with the Cartesian distribution described above. Then an average number of blocks per node is computed. In each dimension, the algorithm starts with the first block in that dimension and determines whether the number of blocks exceeds the mean and if so, blocks are passed on to the next node in line. This can be visualized as a rake passing over each node and dragging excess work into the next available “hole”. In POP, additional constraints are placed on block movements, so that if a block is moved once, it is not likely to be moved again and another block is chosen for the move. This attempts to keep neighboring blocks as local as possible. Because the rake is performed in each logical direction separately, there are instances where the algorithm can actually result in a worse load imbalance as blocks get raked into a corner. Also, if there are relatively few blocks per node, load balancing is not very effective.

11.1.3. Space-filling Curve Distribution

A third distribution method is a partitioning algorithm based on space-filling curves (SFC). This algorithm draws a curve through the set of blocks in a way that maintains locality and can more effectively balance both the computational load and the communications costs. Details can be found in Dennis (2007). While the space-filling curve algorithm is supported for all resolutions, it is particularly effective at reducing the cost of high resolution simulations. Use of the SFC partitioning places restrictions on the size sub-blocks used by POP. In particular the size of the sub-blocks nxBlock and nyBlock must be selected under the constraint that nbx = nxGlobal/nxBlock and nby = nyGlobal/nxBlock are integers such that nbx, nby \in 2^n 3^m 5^p where n, m, and p are integers and nxGlobal and nyGlobal refer to the global size of the horizontal computational grid.

As mentioned above, in addition to improving the computational load balance, SFC partitioning also has a tendency to improve the communication load balance of POP (Dennis 2007). In particular it randomizes network traffic on machines with fixed interconnect topologies, resulting in decreased network congestion and communication time.

11.2. Parallel I/O

In a parallel implementation, the lack of a parallel I/O subsystem creates a serial bottleneck that limits performance. Originally, POP created a parallel I/O system by gathering a single horizontal slice of a 3-D array onto one of the I/O tasks and that task wrote the slice into a specified record of a binary direct-access file. This system was not always robust, due to local machine I/O issues (e.g. system I/O buffers being overwritten by multiple writers) and because this limits the number of I/O tasks and requires each task to keep a full global slice of the data. Finally, this implementation did not support netCDF.

A new parallel I/O library, PIO, has been developed as a collaborative effort by NCAR/CISL, DOE/SciDAC and NCAR/CSEG. PIO was initially designed to allow 0.1^\circ POP to execute and write history and restart files on Blue Gene/L using less than 256 MB per MPI task due to better memory management, relaxing the requirement for retaining a full global slice on each I/O task.

Since that initial prototype version, PIO has developed into a general purpose parallel I/O library that currently supports netCDF (serial), pnetCDF and MPI_IO and has been implemented throughout the entire CCSM system. PIO is a software interface layer designed to encapsulate the complexities of parallel IO and to make it easier to replace the lower level software backend. PIO calls are collective, an MPI communicator is set in a call to PIO_init and all tasks associated with that communicator must participate in all subsequent calls to PIO.

One of the key features of PIO is that it takes the model’s decomposition and redistributes it to an I/O “friendly” decomposition on the requested number of I/O tasks. In using the PIO library for netCDF or pnetCDF I/O, the user must specify the number of iotasks to be used, the stride or number of tasks between iotasks and whether the I/O will use the serial netCDF or pnetCDF library. By varying the number of iotasks, the user can easily reduce the serial I/O memory bottleneck (by increasing the number of iotasks), even with the use of serial netCDF.

12. References

Adcroft, A., J.-M. Campin, P. Heimbach, C. Hill, and J. Marshall. 2005. MITgcm Manual. http://mitgcm.org/pelican/online\_documents/manual.html.

Adcroft, A., C. Hill, and J. Marshall. 1997. “Representation of Topography by Shaved Cells in a Height Coordinate Ocean Model.” Mon. Wea. Rev. 125: 2293–2315.

Arakawa, A., and V.R. Lamb. 1977. “Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model.” Methods of Computational Physics 17: 173.

Arfken, G. 1970. Mathematical Methods for Physicists. Academic, New York.

Boccaletti, Giulio, Raffaele Ferrari, and Baylor Fox-Kemper. 2007. “Mixed Layer Instabilities and Restratification.” J. Phys. Oceanogr. 37: 2228–50. doi:10.1175/JPO3101.1.

Briegleb, B.P., G. Danabasoglu, and W.G. Large. 2010. An over Flow Parameterization for the Ocean Component of the Community Climate System Model. NCAR Tech. Note NCAR/TN-xxxx + STR. National Center for Atmospheric Research.

Brown, J.A., and K. A. Campana. 1978. “An Economical Time-Differencing Scheme for Numerical Weather Prediction.” Mon. Weather Rev. 106: 1125–36.

Bryan, K. 1969. “A Numerical Method for the Study of the World Ocean.” J. Comp. Phys. 4: 347–76.

———. 1984. “Accelerating the Convergence to Equilibrium of Ocean-Climate Models.” J. Phys. Oceanogr. 14: 666–73.

Capet, X., E. J. Campos, and A. M. Paiva. 2008. “Submesoscale Activity over the Argentinian Shelf.” Geophysical Research Letters 35. doi:ARTN L15605.

Capet, X., J. C. McWilliams, M. J. Molemaker, and A. F. Shchepetkin. 2008. “Mesoscale to Submesoscale Transition in the California Current System. Part III: Energy Balance and Flux.” J. Phys. Oceanogr. 38: 2256–69. doi:DOI 10.1175/2008JPO3810.1.

Cox, M.D. 1987. “Isopycnal Diffusion in a Z-Coordinate Ocean Model.” Ocean Modelling 74: 1–5.

Danabasoglu, G. 2004. “A Comparison of Global Ocean General Circulation Model Solutions Obtained with Synchronous and Accelerated Integration Methods.” Ocean Modelling 7: 323–41.

Danabasoglu, G., and J. Marshall. 2007. “Effects of Vertical Variations of Thickness Diffusivity in an Ocean General Circulation Model.” Ocean Modelling 18: 122–41, doi:10.1016/j.ocemod.2007.03.006.

Danabasoglu, G., and J. C. McWilliams. 1995. “Sensitivity of the Global Ocean Circulation to Parameterizations of Mesoscale Tracer Transports.” J. Climate 8: 2967–87.

Danabasoglu, G., R. Ferrari, and J.C. McWilliams. 2008. “Sensitivity of an Ocean General Circulation Model to a Parameterization of Near-Surface Eddy Fluxes.” J. Climate 21: 1192–1208, doi:10.1175/2007JCLI1508.1.

Danabasoglu, G., W.G. Large, and B.P. Briegleb. 2010. “Climate Impacts of Parameterized Nordic Sea Overflows.” Submitted.

Danabasoglu, G., W.G. Large, J.J. Tribbia, P.R. Gent, B.P. Briegleb, and J.C. McWilliams. 2006. “Diurnal Coupling in the Tropical Oceans of CCSM3.” J. Climate 19: 2347–65.

Dennis, J. M. 2007. “Inverse Space-Filling Curve Partitioning of a Global Ocean Model.” In Proceedings of IEEE International Parallel and Distributed Processing Symposium.

Dewar, W., Y. Hsueh, T. J. McDougall, and D. Yuan. 1998. “Calculation of Pressure in Ocean Simulations.” J. Phys. Oceanogr. 28: 577–88.

Dukowicz, J. K., R. D. Smith, and R. C. Malone. 1993. “A Reformulation and Implementation of the Bryan-Cox-Semtner Ocean Model on the Connection Machine.” J. Atmos. Ocean Tech. 14: 294–317.

Dukowicz, J.K. 2001. “Reduction of Density and Pressure Gradient Errors in Ocean Simulations.” J. Phys. Oceanogr. 31: 1915–21.

Dukowicz, J.K., and R.D. Smith. 1994. “Implicit Free-Surface Formulation of the Bryan-Cox-Semtner Ocean Model.” J. Geophys. Res. 99: 7991–8014.

———. 1997. “Stochastic Theory of Compressible Turbulent Fluid Transport.” Phys. Fluids 9: 3523–29.

Dutay, J.-C., J. Bullister, S.Doney, J. Orr, R. Naijar, K. Caldeira, J.-M. Campin, et al. 2002. “Evaluation of Ocean Model Ventilation with CFC-11: Comparison of 13 Global Ocean Models.” Ocean Modelling 4: 89–120.

D’Azevedo, E.F., V.L. Eijkhout, and C.H. Romine. 1993. LAPACK Working Note 56. Conjugate Gradient Algorithms with Reduced Synchronization Overhead on Distributed Memory Multiprocessors. Tech. Rep. CS-93-185. Knoxville, TN: Computer Science Department, University of Tennessee.

Eden, C., and R.J. Greatbatch. 2008. “Towards a Mesoscale Eddy Closure.” Ocean Modelling 20: 223–39.

Eden, C., M. Jochum, and G. Danabasoglu. 2009. “Effects of Different Closures for Thickness Diffusivity.” Ocean Modelling 26: 47–59.

England, M. 1995. “The Age of Water and Ventilation Timescales in a Global Ocean Model.” J. of Phys. Oceanogr. 25: 2756–77.

Farrow, D.E., and D.P. Stevens. 1995. “A New Tracer Advection Scheme for Bryan and Cox Type Ocean General Circulation Models.” J. Phys. Oceanogr. 25: 1731–41.

Feistel, R., and E. Hagen. 1995. “On the GIBBS Thermodynamic Potential of Seawater.” Prog. in Oceanogr. 36: 249–327.

Ferrari, R, and D L Rudnick. 2000. “Thermohaline Variability in the Upper Ocean.” Journal of Geophysical Research-Oceans 105: 16857–83.

Ferrari, R, and W R Young. 1997. “On the Development of Thermohaline Correlations as a Result of Nonlinear Diffusive Parameterizations.” Journal of Marine Research 55: 1069–1101.

Ferrari, R., J.C. McWilliams, V.M. Canuto, and M. Dubovikov. 2008. “Parameterization of Eddy Fluxes Near Oceanic Boundaries.” J. Climate 21: 2770–89, doi:10.1175/2007JCLI1510.1.

Ferreira, D., and J. Marshall. 2006. “Formulation and Implementation of a ‘Residual-Mean’ Ocean Circulation Model.” Ocean Modelling 13: 86–107.

Ferreira, D., J. Marshall, and P. Heimbach. 2005. “Estimating Eddy Stresses by Fitting Dynamics to Observations Using a Residual-Mean Ocean Circulation Model and Its Adjoint.” J. Phys. Oceanogr. 35: 1891–1910.

Fofonoff, N.P., and R. C. Millard. 1983. “Algorithms for Computation of Fundamental Properties of Seawater.” UNESCO Technical Papers in Marine Science 44: 53.

Fox-Kemper, B., G. Danabasoglu, R. Ferrari, and R. W. Hallberg. 2008. “Parameterizing Submesoscale Physics in Global Climate Models.” CLIVAR Exchanges 13: 3–5.

Fox-Kemper, B., G. Danabasoglu, R. Ferrari, S. M. Griffies, R. W. Hallberg, M. M. Holland, S. Peacock, and B. L. Samuels. 2010. “Parameterization of Mixed Layer Eddies. III: Global Implementation and Impact on Ocean Climate Simulations.”

Fox-Kemper, Baylor, and Raffaele Ferrari. 2008. “Parameterization of Mixed Layer Eddies. Part II: Prognosis and Impact.” J. Phys. Oceanogr. 38: 1166–79. doi:10.1175/2007JPO3788.1.

Fox-Kemper, Baylor, Raffaele Ferrari, and Robert Hallberg. 2008. “Parameterization of Mixed Layer Eddies. Part I: Theory and Diagnosis.” J. Phys. Oceanogr. 38: 1145–65. doi:10.1175/2007JPO3792.1.

Garrett, C.J.R., and W. H. Munk. 1975. “Space-Time Scales of Internal Waves: A Progress Report.” J. Geophys. Res. 80: 291–97.

Gent, P.R., and J. C. McWilliams. 1990. “Isopycnal Mixing in Ocean Circulation Models.” J. Phys. Oceanogr. 20: 150–55.

Gerdes, R., C. Koberle, and J. Willebrand. 1991. “The Influence of Numerical Advection Schemes on the Results of Ocean General Circulation Models.” Climate Dynamics 5: 211–26.

Gregg, M.C. 1987. “Diapynal Mixing in the Thermocline: A Review.” J. Geophys. Res. 92: 5249–86.

———. 1989. “Scaling Turbulent Dissipation in the Thermocline.” J. Geophys. Res. 94: 9686–98.

Griffies, S.M. 1998. “The Gent-McWilliams Skew Flux.” J. Phys. Oceanogr. 28: 831–41.

Griffies, S.M., and R. W. Hallberg. 2000. “Biharmonic Friction with a Smagorinsky-Like Viscosity for Use in Large-Scale Eddy-Permitting Ocean Models.” Monthly Weather Review 128: 2935–46.

Griffies, S.M., A. Gnanadesikan, R. C. Pacanowski, V. Larichev, J. K. Dukowicz, and R. D. Smith. 1998. “Isoneutral Diffusion in a z-Coordinate Ocean Model.” J. Phys. Oceanogr. 28: 805–30.

Griffies, S.M., R. C. Pacanowski, M. Schmidt, and V. Balaji. 2001. “Tracer Conservation with an Explicit Free Surface Method for Z-Coordinate Ocean Models.” Monthly Weather Review 129: 1081–98.

Haine, T. W. N., and J. C. Marshall. 1998. “Gravitational, Symmetric and Baroclinic Instability of the Ocean Mixed Layer.” J. Phys. Oceanogr. 28: 634–58.

Haltiner, G.J., and R.T. Williams. 1980. Numerical Prediction and Dynamic Meteorology. 2nd ed. John Wiley; Sons, New York.

Hecht, M.W., W.R. Holland, and P.J. Rasch. 1995. “Upwind-Weighted Advection Schemes for Ocean Tracer Transport: An Evaluation in a Passive Tracer Context.” J. Geophys. Res. 100: 20763–78.

Hibiya, T., M. Nagasawa, and Y. Niwa. 2006. “Global Mapping of Diaplycnal Diffusivity in the Deep Ocean Based on the Results of Expendable Current Profiler (XCP) Surveys.” Geophys. Res. Lett. 33: L03611, doi:10.1029/2005GL025218.

Holland, W.R., F.O. Bryan, and J.C. Chow. 1998. “Application of a Third-Order Upwind Scheme in the NCAR Ocean Model.” J. Climate 11: 1487–93.

Hosegood, P., M. C. Gregg, and M. H. Alford. 2006. “Sub-Mesoscale Lateral Density Structure in the Oceanic Surface Mixed Layer.” Geophys. Res. Lett. 33: doi:10.1029/2006GL026797.

Hundsdorfer, W., and R.A. Trompert. 1994. “Method of Lines and Direct Discretization: A Comparison for Linear Advection.” Appl. Numer. Math. 13: 469–90.

Jackett, D.R., and T. J. McDougall. 1995. “Minimal Adjustment of Hydrographic Profiles to Achieve Static Stability.” J. of Atmos. and Oceanic Tech. 12: 381–89.

Jackett, D.R., T.J. McDougall, R. Feistel, D.G. Wright, and S. Griffies. 2006. “Algorithms for Density, Potential Temperature, Conservative Temperature, and the Freezing Temperature of Seawater.” J. of Atmos. and Oceanic Tech. 23: 1709–28.

Jayne, S.R. 2009. “The Impact of Abyssal Mixing Parameterizations in an Ocean General Circulation Model.” J. Phys. Oceanogr. 39: 1756–75.

Jayne, S.R., and L. C. St. Laurent. 2001. “Parameterizing Tidal Dissipation over Rough Topography.” Geophys. Res. Lett. 28: 811–14, doi:10.1029/2000GL012044.

Jochum, M. 2009. “Simulated Climate Impacts of Latitudinal Variations in Diapycnal Diffusivity.” J. Geophys. Res. 114: C01010, doi:10.1029/2008JC005030.

Jochum, M., and J. Potemra. 2008. “Sensitivity of Tropical Rainfall to Banda Sea Diffusivity in the Community Climate System Model.” J. Climate 21: 6445–54.

Jochum, M., G. Danabasoglu, M. Holland, Y.-O. Kwon, and W.G. Large. 2008. “Ocean Viscosity and Climate.” J. Geophys. Res. 113: C06017, doi:10.1029/2007JC004515.

Killworth, P.D., D. Stainforth, D.J. Webb, and S.M. Paterson. 1991. “The Development of a Free-Surface Bryan-Cox-Semtner Ocean Model.” J. Phys. Oceanogr. 21: 1333–48.

Large, W.G., G. Danabasoglu, S.C. Doney, and J.C. McWilliams. 1997. “Sensitivity to Surface Forcing and Boundary Layer Mixing in the NCAR CSM Ocean Model: Annual-Mean Climatology.” J. Phys. Oceanogr. 27: 2418–47.

Large, W.G., G. Danabasoglu, J.C. McWilliams, P.R. Gent, and F.O. Bryan. 2001. “Equatorial Circulation of a Global Ocean Climate Model with Anisotropic Horizontal Viscosity.” J. Phys. Oceanogr. 31: 518–36.

Laurent, L.C. St., and C. J. R. Garrett. 2002. “The Role of Internal Tides in Mixing the Deep Ocean.” J. Phys. Oceanogr. 32: 2882–99.

Laurent, L.C. St., and J. D. Nash. 2004. “An Examination of the Radiative and Dissipative Properties of Deep Ocean Internal Tides.” Deep-Sea Res. II 51: 3029–42, doi:10.1016/j.dsr2.2004.09.008.

Laurent, L.C. St., H. L. Simmons, and S. R. Jayne. 2002. “Estimates of Tidally Driven Enhanced Mixing in the Deep Ocean.” Geophys. Res. Lett. 29: 2106, doi:10.1029/2002GL015633.

Leclair, M., and G. Madec. 2009. “A Conservative Leapfrog Time Stepping Method.” Ocean Modelling, doi 10.1016/j.ocemod.2009.06.006.

Leonard, B.P. 1979a. “A Stable and Accurate Convective Modeling Procedure Based on Quadratic Upstream Interpolation.” Comp. Meth. in Appl. Eng. 19: 59–98.

———. 1979b. “The ULTIMATE Conservative Difference Scheme Applied to Unsteady One-Dimensional Advection.” Comp. Meth. in Appl. Mech. Eng. 88: 17–74.

Madec, G., and M. Imbard. 1996. “A Global Ocean Mesh to Overcome the North Pole Singularity.” Climate Dyn. 12: 381–88.

Madec, G., M. Imbard, and C. Levy. 1998. OPA 8.1 Ocean General Circulation Model Reference Manual. Note du Pole der modelisation number 11. France: Institute Pierre-Simon Laplace (IPSL).

Mahadevan, A., A. Tandon, and R. Ferrari. 2009. “Rapid Changes in Mixed Layer Stratification Driven by Submesoscale Instabilities and Winds.”

Maltrud, M., R.D. Smith, A.J. Semtner, and R.C. Malone. 1998. “Global Eddy-Resolving Ocean Simulations Driven by 1985-1995 Atmospheric Winds.” J. Geophys. Res. 103: 825–30.

McDougall, T.J. 1987. “Neutral Surfaces.” J. Phys. Oceanogr. 17: 1950–64.

McDougall, T.J., D.R. Jackett, D.G. Wright, and R. Feistel. 2003. “Accurate and Computationally Efficient Algorithms for Potential Temperature and Density of Seawater.” J. of Atmos. and Oceanic Tech. 20: 730–41.

Murray, R.J. 1996. “Explicit Generation of Orthogonal Grids for Ocean Models.” J. Comp. Phys. 126: 251–73.

NGDC. 2006. Two-Minute Gridded Global Relief Data (ETOPO2v2). National Geophysical Data Center. National Oceanic; Atmospheric Administration, U.S. Department of Commerce.

OCMIP. 2000. “Ocean Carbon Model Intercomparison Project.” http://www.ipsl.jussieu.fr/OCMIP/phase2/simulations/CFC/HOWTO-CFC.html.

Ohlmann, J.C. 2003. “Ocean Radiant Heating in Climate Models.” J. Climate 16: 1337–51.

Osborn, T.R. 1980. “Estimates of the Local Rate of Vertical Diffusion from Dissipation Measurements.” J. Phys. Oceanogr. 10: 83–89.

Oschlies, A. 2002. “Improved Representation of Upper-Ocean Dynamics and Mixed Layer Depths in a Model of the North Atlantic on Switching from Eddy-Permitting to Eddy-Resolving Grid Resolution.” J. Phys. Oceanogr. 32: 2277–98.

Pacanowski, R.C., and A. Gnanadesikan. 1998. “Transient Response in a z-Level Ocean Model That Resolves Topography with Partial Cells.” J. Phys. Oceanogr. 126: 3248–70.

Pacanowski, R.C., and S.M. Griffies. 2000. The MOM3 Manual. Princeton, NJ, USA: Geophysical Fluid Dynamics Laboratory/NOAA.

Pacanowski, R.C., and S.G.H. Philander. 1981. “Parameterization of Vertical Mixing in Numerical Models of the Tropical Oceans.” J. Phys. Oceanogr. 11: 1443–51.

Price, J., and M.O. Baringer. 1994. “Outflows and Deep Water Production by Marginal Seas.” Prog. Oceanogr. 33: 161–200.

Price, J., and J. Yang. 1998. “Marginal Sea Overflows for Climate Simulations.” In Ocean Modeling and Parameterization, edited by E. P. Chassignet and J. Verron, 155–70. Kluwer Academic.

Redi, M.H. 1982. “Oceanic Isopycnal Mixing by Coordinate Rotation.” J. Phys. Oceanogr. 12: 1154–58.

Rudnick, D L, and R Ferrari. 1999. “Compensation of Horizontal Temperature and Salinity Gradients in the Ocean Mixed Layer.” Science 283: 526–29.

Schopf, P.S., and A. Loughe. 1995. “A Reduced-Gravity Isopycnal Ocean Model: Hindcasts of El Nino.” Monthly Weather Review 123: 2839–63.

Semtner, A.J. 1986. “Finite-Difference Formulation of a World Ocean Model.” In Advanced Physical Oceanographic Numerical Modelling, edited by J.J. O’Brien, 187–202. Norwell, Mass.: D. Reidel.

Simmons, H.L., S. R. Jayne, L. C. St. Laurent, and A. J. Weaver. 2004. “Tidally Driven Mixing in a Numerical Model of the Ocean General Circulation.” Ocean Modelling 6: 245–63, doi:10.1016/S1463–5003(03)00011–8.

Smagorinsky, J. 1993. “Some Historical Remarks on the Use of Nonlinear Viscosities.” In Large Eddy Simulation of Complex Engineering and Geophysical Flows, edited by B. Galperin and S.A. Orszag. Cambridge University Press.

Smith, R.D. 1999. “The Primitive Equations in the Stochastic Theory of Adiabatic Stratified Turbulence.” J. Phys. Oceanogr. 29: 1865–80.

Smith, R.D., and J. C. McWilliams. 2003. “Anisotropic Horizontal Viscosity for Ocean Models.” Ocean Modelling 5: 129–56.

Smith, R.D., J.K. Dukowicz, and R.C. Malone. 1992. “Parallel Ocean General Circulation Modeling.” Physica D. 60: 38–61.

Smith, R.D., S. Kortas, and B. Meltz. 1995. Curvilinear Coordinates for Global Ocean Models. LA-UR-95-1146. Los Alamos National Laboratory.

Smith, W.H.F., and D. T. Sandwell. 1997. “Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings.” Science 277: 1956–62, doi:10.1126/science.277.5334.1956.

St.Laurent, L., and R.W. Schmitt. 1999. “The Contribution of Salt Fingers to Vertical Mixing in the North Atlantic Tracer Release Experiment.” J. Phys. Oceanogr. 29: 1404–24.

Thiele, G., and J. Sarmiento. 1990. “Tracer Dating and Ocean Ventilation.” J. Geophys. Res. 95: 9377–91.

Thomas, L N, A Tandon, and A Mahadevan. 2008. “Submesoscale Processes and Dynamics.” In Ocean Modeling in an Eddying Regime, edited by M. Hecht and H. Hasumi, 177:17–38. AGU Geophysical Monograph Series.

Thomas, L.N. 2005. “Destruction of Potential Vorticity by Winds.” J. Phys. Oceanogr. 35: 2457–66.

———. 2008. “Formation of Intrathermocline Eddies at Ocean Fronts by Wind-Driven Destruction of Potential Vorticity.” Dyn. Atmos. Ocean. 45: 252–73.

Thomas, L.N., and C. M. Lee. 2005. “Intensification of Ocean Fronts by down-Front Winds.” J. Phys. Oceanogr. 35: 1086–1102.

Thomas, Leif N., and Raffaele Ferrari. 2008. “Friction, Frontogenesis, and the Stratification of the Surface Mixed Layer.” J. Phys. Oceanogr. 38: 2501–18. doi:DOI 10.1175/2008JPO3797.1.

Visbeck, M., J. Marshall, T. Haine, and M. Spall. 1997. “Specification of Eddy Transfer Coefficients in Coarse-Resolution Ocean Circulation Models.” J. Phys. Oceanogr. 27: 381–402.

Wajsowicz, R.C. 1993. “A Consistent Formulation of the Anisotropic Stress Tensor for Use in Models of the Large-Scale Ocean Circulation.” J. Comp. Phys. 10: 333–38.

Walker, S.J., R.F. Weiss, and P.K. Salameh. 2009. “Reconstructed Histories of the Annual Mean Atmospheric Mole Fractions for the Halocarbons CFC-11, CFC-12, CFC-113 and Carbon Tetrachloride.” http://bluemoon.ucsd.edu/pub/cfchist/.

Webb, D.J. 1995. “The Vertical Advection of Momentum in Bryan-Cox-Semtner Ocean General Circulation Models.” J. Phys. Oceanogr. 25: 3186–95.

Whitehead, J.A., A. Leetmaa, and R.A. Knox. 1974. “Rotating Hydraulics of Strait and Sill Flows.” Geophys. Fluid Dyn. 6: 101–25.

Williams, G.P. 1972. “Friction Term Formulation and Convective Instability in a Shallow Atmosphere.” J. Atmos. Sci. 29: 870–76.

Wu, W., G. Danabasoglu, and W.G. Large. 2007. “On the Effects of Parameterized Mediterranean Overflow on North Atlantic Ocean Circulation and Climate.” Ocean Modelling 19: 31–52, doi:10.1016/j.ocemod.2007.06.003.