H. \underline{Numerical and Coupling Improvements}

a. Third-Order Upwind Algorithm

Third-order upwind schemes and variants of this scheme are discussed fully in Leonard (1979). A somewhat more complicated scheme has been put into the MOM code by the FRAM group, see Farrow and Stevens (1995). The third-order upwind scheme is switched on by setting the ifdef option upwind3.

The finite-difference expression for the advection of tracers is given by

\begin{aligned}ADV (i, j, k) = & − \frac{1}{dxt_{i} cst_{j}} (\overline{u}_{E} T_{E}^{∗} − \overline{u}_{W} T_{W}^{∗}) − \frac{1}{dyt_{j}} (\overline{v}_{N} T_{N}^{∗} − \overline{v}_{S} T_{S}^{∗}) \\ & − \frac{1}{dzt_{k}} (w_{k−1} T_{T}^{∗} − w_{k} T_{B}^{∗}) \end{aligned}

where \overline{u}_{W} , \overline{u}_{E}, \overline{v}_{S}, \overline{v}_{N} are the T-grid cell face centered velocity components.

\begin{aligned}\overline{u}_{E}(i) & = \frac{u_{i,j}dyu_{j} + u_{i,j−1}dyu_{j−1}}{2dyt_{j}} \overline{u}_{W} (i) = \overline{u}_{E}(i − 1) \\ \overline{v}_{N}(j) & = \frac{(v_{i,j}dxu_{i} + v_{i−1,j}dxu_{i−1})csu_{j}}{2dxt_{i}cst_{j}} \\ \overline{v}_{S}(j) & = \frac{(v_{i,j−1}dxu_{i} + v_{i−1,j−1}dxu_{i−1})csu_{j−1}}{2dxt_{i}cst_{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 MOM, these are taken as the arithmetic average of the tracer concentrations in adjacent cells, e.g.,

T_{E}^{∗} = \frac{1}{2}(T_{i+1,j} + T_{i,j})

In the third-order upwind formulation, the T^{∗}s are replaced by an upwind direction biased three-point interpolation to the interface. In the zonal direction they are:

T_{E}^{∗}(i) = \{ \begin{matrix}γT_{i−1} + β^{+}T_{i} + α^{+}T_{i+1} \overline{u}_{E} > 0 \\ β^{−}T_{i} + α^{−}T_{i+1} + δT_{i+2} \overline{u}_{E} < 0 \end{matrix}

T_{W}^{∗} (i) = T_{E}^{∗}(i − 1)

H-1 where the coefficients α, β, γ, δ are given by

\begin{aligned}α^{+} & = \frac{dxt_{i} (2dxt_{i} + dxt_{i−1})}{(dxt_{i} + dxt_{i+1}) (dxt_{i+1} + dxt_{i−1} + 2dxt_{i})} \\ α^{−} & = \frac{dxt_{i}(2dxt_{i+1} + dxt_{i+2})}{(dxt_{i} + dxt_{i+1}) (dxt_{i+1} + dxt_{i+2})} \\ β^{+} & = \frac{dxt_{i+1} (2dxt_{i} + dxt_{i−1})}{(dxt_{i} + dxt_{i+1}) (dxt_{i} + dxt_{i−1})} \\ β^{−} & = \frac{dxt_{i+1}(2dxt_{i+1} + dxt_{i+2})}{(dxt_{i} + dxt_{i+1}) (dxt_{i} + dxt_{i+2} + 2dxt_{i+1})} \\ γ & = \frac{−dxt_{i} ∗ dxt_{i+1}}{(dxt_{i} + dxt_{i−1}) (dxt_{i+1} + dxt_{i−1} + 2dxt_{i})} \\ δ & = \frac{−dxt_{i} ∗ dxt_{i+1}}{(dxt_{i+1} + dxt_{i+2}) (dxt_{i} + dxt_{i+2} + 2dxt_{i+1})} \end{aligned}

The interface tracer values T_{N}^{∗}, T_{S}^{∗}, T_{T}^{∗}, andT_{B}^{∗} are obtained in an analogous manner.

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–1741.

Leonard, B.P., 1979: A stable and accurate convective modeling procedure based on quadratic upstream interpolation. Comp. Meth. in Appl. Eng., 19, 59–98.

b. Single Vertical Velocity

The second improvement concerns the vertical velocity computation. In the original scheme, two independent vertical velocity computations, one in the momentum and one in the tracer equations, are performed, resulting in two different vertical velocity fields, with the former typically much noisier than the latter. In the new method, the vertical velocity is computed only for the tracer equations. Consequently, this velocity is then area-averaged to obtain the vertical velocity for the momentum equations, eliminating the second computation. The related advection operator in the momentum equations is also modified accordingly to maintain the flux and kinetic energy conservation properties. This new treatment of the vertical velocity produces a unique and smooth vertical velocity field with rather modest effects on global distributions. Further details can be found in Webb (1995), where the new scheme is described in Section 7.

Webb, D.J., 1995: The vertical advection of momentum in Bryan-Cox-Semtner ocean general circulation models. J. Phys. Oceanogr., 25, 3186–3195.

H-2 c. North Pole Implementation

This section describes a consistent way to treat the north pole in B-grid models. The traditional finite difference approach breaks down because the zonal grid spacing goes to zero and there is no set of values to the north of the north pole. The approach described here treats the north pole as a circular grid point of radius dyt_{jmt} which carries values for all tracers and the transport streamfunction at its center. All grid points to the south of the pole are arranged in a standard B-grid, with internal velocity points offset one half grid spacing to the north and east from tracer and transport streamfunction points. There is no velocity point which corresponds to the tracer point at the pole. This option is activated by specifying the ifdef option northpole.

The momentum equations may be solved at the row just to the south of the north pole (jmt-1) using the standard solution procedure outlined in Section C. It is assumed here that a Laplacian frictional parameterization is used, some considerations may be necessary in order not to overflow arrays when using biharmonic or higher frictional parameterizations. The value of the internal mode velocities at j=jmt is not defined, but this poses no problem because the area of the interface between the jmt-1 and jmt velocity grid boxes is zero. All terms in the finite difference formulation which multiply u_{i,jmt} or v_{i,jmt} are also multiplied by cst_{jmt} (the cosine at the latitude of the the j=jmt T-point) which, for the north pole, is zero.

The solution of the tracer equation at the north pole requires some special treatment. The finite difference form of the tracer evolution equation at the circular grid may be written (neglecting friction) as

T^{n+1} = T^{n−1} + 2Δt[ ∑_{i=2}^{imt−1}(v_{i,jmt−1}dxu_{i} + v_{i−1,jmt−1}dxu_{i−1})

(T_{i,jmt} + T_{i,jmt−1})csu_{jmt−1}/(dyt_{jmt}dxt_{i})].

This is simply the sum of all meridional fluxes into the circular grid cell from the regular grid to the south. All terms on the right hand side are known. Once solved, the single new tracer value at the north pole can be filled into the j=jmt row of the tracer array so that the meridional fluxes through the interfaces with the interior grid can be calculated during the next time step.

The transport streamfunction is located at the same points as the tracers, so this will also require some special treatment at the pole. We take a similar approach by vertically averaging the internal mode zonal momentum equation and computing its line integral along the velocity points at row j=jmt-1. The contribution due to the surface pressure gradient identically vanishes, leaving

∫ \overline{u}_{t}dλ = ∫ \overline{u′}_{t}dλ.

H-3 where \overline{u′_{t}} is the vertical average of the internal mode tendency and λ is longitude. Using the definition of the transport streamfunction, this relation may be rewritten as

− ∫ \frac{1}{aH} Ψ̇_{θ}dλ = ∫ \overline{u′}_{t}dλ,

where a is the mean radius of the earth and H is the model depth at tracer points.

Hence, the barotropic circulation tendency around the pole determines the value of the streamfunction tendency (Ψ̇) at the pole, relative to the average interior value. This equation may then be used to solve for the time rate of change of the transport streamfunction at the north pole (Ψ̇_{pole}) as

∑_{i=2}^{imt−1} \frac{2dxt_{i}}{H_{i,jmt−1} + H_{i−1,jmt−1}} Ψ̇_{pole} =

∑_{i=2}^{imt−1} \frac{2dxt_{i}}{H_{i,jmt−1} + H_{i−1,jmt−1}} Ψ̇_{i,jmt−1} − ∑_{i=2}^{imt−1} dxt_{i}dyu_{jmt−1}\overline{u′_{t}}.

The forcing term on the right hand side is known and an estimate of the streamfunction tendency at j=jmt-1 is obtained during the standard iterative solution for the streamfunction in the interior. Once the estimate for the tendency at the pole is found, the usual iterative solution procedure is resumed. As for the tracer arrays, the single values for the streamfunction tendency and the streamfunction, once solved, must be stored in the appropriate array so the interior solution can proceed.

REMARK: The present implementation of the northpole ifdef option assumes equal grid intervals in the zonal direction when computing fluxes into the north pole gridbox.

d. Tapering of Horizontal and Isopycnal Mixing Coefficients

Our previous experience with global models showed that the diffusive computational stability at high latitudes, rather than advective stability, is usually the limiting factor for the model time step. Consequently, in the present model, we allow the horizontal eddy viscosity coefficients A_{MH}, A_{I}, and A_{IT D} to vary in the meridional direction. In particular, they are reduced at high latitudes. This treatment necessitates no changes in the scalar tracer equations. In contrast, the vector momentum equation requires additional horizontal friction terms so that pure solid body rotation does not produce any stress on the fluid. Consequently, we modified these viscous terms to follow Wajsowicz (1993). Because the model uses the leap-frog time integration scheme, the viscous/diffusive computational stability limits in the zonal direction in which the grid spacing becomes smaller at high latitudes are given by

A_{MH} \begin{matrix}Δt & U \\ (Δx cosϕ) & 2 \end{matrix} ≤ \frac{1}{4}, A_{I} \begin{matrix}Δt & T \\ (Δx cosϕ) & 2 \end{matrix} ≤ \frac{1}{4}.

H-4 In the present discussion, we consider only this direction for simplicity, because it is the most restrictive. In the above equations, Δt_{U} and Δt_{T} are the maximum momentum and tracer time steps, respectively, Δx is the longitudinal grid spacing, and ϕ is latitude. The model determines the allowable A_{MH}, A_{I}, and A_{IT D} from these stability equations, so that all these coefficients are kept at their full values as poleward as possible. In addition, to prevent excessively small A_{MH}, A_{I}, and A_{IT D} which can produce instabilities due to lack of enough diffusion, we restrict the minimum values to be 20% and 6% of their specified values in the 3×3 and 2×2 models, respectively. In the 3×3 model, A_{MH} tapering starts at 81.9^{∘}N, and it is kept constant at the 20% level starting at 87.3^{∘}N. A_{I} and A_{IT D} satisfy the above constraint at all latitudes. In the 2×2 model, A_{MH}, and A_{I} and A_{IT D} taperings start at 84.6^{∘}N and 85.2^{∘}N, respectively. They are kept constant at the 6% level north of 88.8^{∘}N. No tapering is applied in the Southern Hemisphere in both model configurations. This option is activated by setting the ifdef option taperhrzvd.

Wajsowicz, R. C., 1993: A consistent formulation of the anisotropic stress tensor for use in models of the large scale ocean circulation. J. Comput. Phys., 105, 333-338.

e. K. Bryan Acceleration Technique

The model time integration is performed in two stages. First, the acceleration technique of Bryan (1984) is applied in which the momentum and tracer equations have variable time steps, i. e.,

\begin{matrix}α[\frac{û^{n+1} − û^{n−1}}{2Δt} ] = F^{n−1,n}, \\ γ(z) & [\frac{[θ, S]^{n+1} − [θ, S]^{n−1}}{2Δt} ] = G^{n−1,n}. \end{matrix}

The temporal discretization is done using the leap-frog scheme, and n − 1, n, and n + 1 represent the previous, current, and next time levels, respectively. Also, û is the horizontal baroclinic velocity vector, Δt is the fundamental (surface tracer) time step, θ is potential temperature, and S is salinity. F and G represent the rest of the terms in these equations, involving quantities from both the n − 1 and n time levels. α and γ, which can vary in the vertical (z), are the time scale, or acceleration, factors. In the case of α > γ, the momentum time step is smaller than the tracer time step by a factor of γ/α.

The model calender and all the forcing fields are based on the surface tracer time step (γ^{−1}Δt) by setting γ = 1 above 1058 m. The depth and timesteps given here and below apply to the 3x3 model version. With this choice, the annual cycles of heat and salt forcing are consistently applied. Here, Δt = 440 min was chosen, well within the advective computational stability limit. The smaller current speeds at depth are taken advantage of by setting γ = 0.2 between 1058 and 2502 m and γ = 0.1 between 2502 and 5000 m depth. Consequently, the deep ocean tracer time step is about 3 days; ten times larger than at the surface. To illustrate a disadvantage of this depth acceleration, consider a vertical heat flux, Q, between two adjacent model layers, k and k + 1. Integration of the equation above

H-5 leads to a change in the total heat content of the two layers by an amount

ΔH = 2Δt Q (γ_{k}^{−1} − γ_{k+1}^{−1}) .

Conservation of heat requires ΔH = 0, which is violated whenever γ_{k} ≠ γ_{k+1}. Therefore, vertical discontinuities in γ are restricted to depths where the fluxes are thought to be small enough that tracers are well conserved.

In the momentum equations, horizontal diffusion determines computational stability. At high latitudes this factor limits the baroclinic velocity time step to about 44 min, α = 10, at all depths. To be consistent a time step of α^{−1}Δt is also used in the barotropic streamfunction equation. With this time step, the surface wind forcing is not consistently applied, because the annual cycle harmonic is applied 10 times per momentum year. As a result, some distortion of the seasonal cycle of ocean currents is expected. Analogous problems are not serious in the tracer fields, because their depth acceleration (γ ≠ 1) is restricted to depths where the influence of the seasonal cycle on tracers is negligible. The acceleration technique is implemented by setting the values of dtts, dtuv, dtsf, and dtxcel in the model input data.

In the second stage of the time integration, the accelerated results were used as initial conditions for a synchronous time integration in which equal time steps and no depth acceleration were applied, i. e., Δt = 44 min and α = γ = 1. This integration needs to be long enough to remedy any erroneous behavior spawned by violating conservation laws and by forcing mismatches in the acceleration stage, see Danabasoglu et al. (1996).

Bryan, K., 1984: Accelerating the convergence to equilibrium of ocean-climate models. J. Phys. Oceanogr., 14, 666–673.

Danabasoglu, G., J.C. McWilliams, and W.G. Large, 1996: Approach to equilibrium in accelerated global oceanic models. J. Climate, in press.

f. Ocean Surface Slope

The components of the ocean surface slope are computed exactly the same way as the diagnostic surface pressure gradient in the original MOM code. These slope components are accumulated at each timestep and only time-averaged fields over the nysnc steps in a coupling interval are passed to the flux coupler. The surface slope is nondimensionalized from the surface pressure gradient by dividing by the gravitational acceleration. This calculation is activated by setting the ifdef option sflxpvmicetilt.

g. Sea-Ice Formation and Melting

This improvement is needed when the ocean model is coupled to a sea-ice model, and is activated by the ifdef option sflxiceflx. At the end of each timestep, n, after surface fluxes and model transports have been applied, the upper kmxice model layers, where kmxice is a model input parameter, are scanned for temperatures below the freezing point,

H-6 T_{f}. At present T_{f} is set to the constant value of −1.8^{∘}C. However, the freezing point is a function of salinity and a more accurate equation for T_{f}, that will be implemented in the near future, is

\begin{aligned}T_{f} = & − .0575 S_{K} + 0.001710523 S_{K}^{3/2} \\ & − 0.0002154996 S_{K}^{2} , \end{aligned}

at the local salinity, S_{K}. At each layer the potential mass per unit area of ice formation (potice > 0) or ice melt (potice < 0) is computed as

potice (k) = MAX(\frac{(ϱC_{p})_{o}}{Λ_{f}} H_{k} (T_{f} − T_{k}) , qice (k − 1)) ,

where (ϱC_{p})_{o} = 4.1 × 10^{6} j/m^{3}/K, is the product of ocean density and heat capacity, and Λ_{f} = 334000 j/Kg is the latent heat of fusion. Any ice that forms at depth is assumed to float towards the surface and this ice flux, qice(k) (defined positive downwards), so qice(k) ≤ 0 is accumulated bottom to top as

qice (k) = ∑_{kmxice}^{k} −potice (k) ,

where qice(kmxice−1) = 0, assumes no ice formation below the kmxice layer. As ice floats to the surface, it can either partially or completely melt in upper layers whose temperatures are above freezing.

At each layer the temperature and salinity are both adjusted in accord with the ice formed or melted in the layer:

\begin{aligned}T_{k} & = T_{k} + \frac{Λ_{f}}{H_{k}(ϱC_{p})_{o}} potice(k) \\ S_{k} & = S_{k} + \frac{(S_{o} − S_{i})}{ϱ_{o} H_{k}} potice(k) . \end{aligned}

The ice flux is accumulated at each location over the number of timesteps in a coupling interval, nsync, as

aqice(i, j) = ∑_{n=1}^{nsync} qice(k = 1)_{n}.

If, at the end of a vertical scan, the upper layer temperature T_{1} remains greater than freezing, (qice(k = 1) = 0), then the excess heat melts previously formed ice and aqice(i, j), T_{1} and S_{1} are adjusted accordingly. Thus over the nsync timesteps, ice is assumed to remain where it was formed. After nsync timesteps, the accumulated ice flux aqice(i, j) can be passed through the Flux Coupler to the sea-ice model as an equivalent downward heat flux (W/m^{2}).

qflux(i, j) = \frac{−Λ_{f}}{nsync Δt} aqice(i, j)

where aqice includes any adjustment due to melting previously formed ice in the last timestep, and Δt is the model timestep.

H-7

Thus, ocean temperature and salinity adjustments corresponding to the total ice formed (qflux > 0) are made prior to the ice being passed to the sea-ice model. However, warm surface model temperatures result in (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, temperature and salinity adjustments are delayed until the appropriate heat and freshwater fluxes are received back from the sea-ice model through the Flux Coupler. Note, if sufficient ice is present in the absence of all other fluxes, the heat flux will be equal to the cooling needed to make T_{1} = T_{f}, when applied over the next nsync timesteps.

H-8