GFDL Modular Ocean Model (MOM)

Design philosophy

The first numerical implementation of the primitive equations was performed at the Geophysical Fluid Dynamics Laboratory (GFDL) by Bryan (1969) and Cox (1984). The original objective of this development was to build a tool for numerical modeling of the large-scale (basin-wide or global) ocean circulation using then-current algorithms based upon decomposition of the ocean into rectangular boxes of variable size, and the use of low-order finite difference approximations.

By default, the equations are solved in spherical coordinates, and the box boundaries follow latitude circles (), meridians () and geopotential surfaces (z). The geopotential level formulation is a simple and convenient concept that has been used for a variety of applications; however, the primary applications of MOM have been long-term integrations (years to decades and centuries) of the basin-wide wind driven and thermohaline circulation in realistic geometries. A detailed model description and user's guide for the most recent release is available (Pacanowski, 1996) and is being updated continually. Other versions of this family of models are the free sea surface code (Killworth et al., 1991) and the Semtner and Chervin code (POCM and POP; see McClean et al., 1997)

System of equations

For geopotential coordinates, , and the derivatives in the unforced, inviscid primitive equations become the vertical velocity w. Furthermore, the hydrostatic equation assumes its familiar form,

and the pressure gradient terms in the unforced, inviscid primitive equations reduce to

This is because z is orthogonal to the horizontal coordinates x and y.

The hydrostatic primitive equations in MOM are formulated in spherical coordinates, which are (for completeness)

where P is the dynamic pressure (), and D and F symbolically denote dissipation and forcing terms, respectively.

Kinematic boundary conditions in the vertical are

at the surface, and

at the bottom. In order to save computer time, the nonlinear equation of state is implemented for each model level separately by fitting a third-order polynomial approximation of the full equation of state (Bryan and Cox, 1972).

Depth-integrated flow

The two-dimensional circulation (often referred to as the barotropic or external mode) requires a different treatment than the baroclinic or internal mode. Most applications of MOM still use the rigid lid option, wherein the vertically integrated flow is assumed to be nondivergent, and the streamfunction has to be found from an elliptic equation. The mass transport streamfunction () is defined as

where U and V are the depth-averaged velocities. The time-varying streamfunction is determined by first taking the curl of the depth-integrated equations. The resulting vorticity equation is then time-stepped. The depth-integrated velocities can be computed after the elliptic problem is solved, subject to the appropriate sidewall boundary conditions. Several solution algorithms are offered for the elliptic equation: simple 5 or 9 point successive over-relaxation (SOR) schemes, and a conjugate gradient solver. Islands require special treatment, and their number should be chosen as small as possible.

Recently, two additional strategies for the external mode have been implemented; one is based on an implicit free surface (Dukowicz and Smith, 1994), the other on an elliptic equation for dynamic pressure P rather than streamfunction (Dukowicz et al., 1993). While the explicit free surface version (Killworth et al., 1991) requires time splitting (with many small time steps for the external mode centered around the long time-steps for the internal modes), both implicit methods have advantages for treating islands and are potentially more efficient in high-resolution applications with many islands and on massively parallel computer systems.

Spatial discretization, grids and topography

MOM uses a discretization based on the so-called "box-concept", wherein the ocean is subdivided into a number of rectangular boxes. The physical variables are arranged on the Arakawa "B" grid with tracers in the center of the boxes, horizontal velocities at the corners (centered in the vertical) and vertical velocity centered at the bottom and top surfaces of the box. Figure 1 shows two illustrations of vertical discretization in an area with strong variation in water depth. The first example uses levels having equidistant geopotential spacing; and the second has increased resolution near the surface for an improved representation of upper thermocline processes. The latter is typical of most applications (e.g., Cox 1985; Bryan and Holland, 1989; FRAM, 1991; Semtner and Chervin, 1992; Lehmann, 1995).

Figure 1: Vertical disrecretization of a tanh-shaped topogrphy of a geopotential (z) coordinate model and 30 levels: (a) equidistant grid spacing; (b) typical typical discretization with higher resolution near the surface.

Note in figure 1 that the exact form of the topography changes with vertical resolution and placement of the grid boxes, especially in regions with very gentle and very steep topographic slopes. This leads to two separate constraints on the choice of the levels, the representation of topography and the representation of stratification. Some flexibility exists in choosing the levels, although a grid spacing that varies according to a smooth analytic function has some formal advantages (Treguier et al., 1996). However, it is usually not possible to find an equally suitable discretization for both topography and stratification.

In practice, this vertical discretization may lead to a large number of inactive grid cells, which both increases the computer memory requirements and the number of needless operations on vector computers. Note also that this choice of vertical coordinate requires a lateral boundary condition at the topographic mask. MOM adopts a no-slip condition, which is most "natural" for the "B" grid, as the velocity points are located on the boundaries. Unfortunately, this leads to a systematic numerical error associated with baroclinic leakage of energy. In general, optimal representation of topography is achieved with (Winton et al., 1998)

while much smaller values lead to inconsistencies (Hughes, 1995).

Although spherical coordinates are hard-wired into the code, some flexibility can be achieved with the "rotated grid" configuration, for which the poles are moved from their geographical position in such a way that the convergence of the meridians no longer poses a problem in high latitudes (see, e.g., Gerdes and Köberle, 1995).

Finally, the "B" grid requires special attention during the specification of the topography. The numerical grid may be inadequate if one-grid-cell bays or passages exist which are "non-advective" due to the placement of the variables on the "B" grid. This is, for example, crucial in narrow straits and trenches. As Redler and Böning (1997) point out, the deep circulation in a geopotential coordinate model on the "B" grid may depend crucially on artificially widened fracture zones if these are not adequately resolved.

Semi-discrete equations

We define the averaging operators as

and the derivative operators as

The negative sign in the last relation reflects the fact that the vertical index k increases in the negative z-direction. We also introduce the notation for a vertical sum. The semi-discrete equations for MOM are then:

Here, a modified averaging scheme for the vertical momentum advection term (Webb, 1995) has been added to improve the grid dispersion properties of the code

The vertical velocity is computed as

Finally, the barotropic vorticity is

and the barotropic velocities from the streamfunction are

where is the minimum of the four surrounding depth values. Note that on a "B" grid, averaging of the prognostic variables is necessary for both and .

Time-stepping

Time integration in MOM is done with the traditional leapfrog scheme with an occasional Euler forward step to eliminate the computational mode. A sophisticated time managing subroutine in MOM synchronizes model time and calendar time, such that leap years can be taken into account when, for example, daily forcing datasets are used.

Over the years, a number of applications with specific needs have led to several model extensions, which are usually included in the code as independently specifiable options. The most important are described in the following paragraphs.

MOM comes with a number of global datasets for wind stress (Hellerman and Rosenstein, 1983), air temperature climatology (Oort, 1983), ocean potential temperature and salinity climatology (Levitus, 1982), and bathymetry. All the necessary interpolation routines are also provided and facilitate coupled atmosphere-ocean simulations for climate studies.

The treatment of artificial domain boundaries is a concern for all non-global simulations. Often, a simple closed boundary will be sufficient, if placed along gyre boundaries where inter-gyre exchange is weak. In addition, the assumption of a "sponge layer", where water mass transformations take place, has been used extensively.

In instances of strong advective exchange, the assumption of no normal flow will be unrealistic. Therefore, a method is required to successfully allow for an inflow (with prescribed water mass characteristics) and an outflow dependent on internal model dynamics. Stevens (1990) devised such a method, based on Orlanski's (1976) radiation condition, for any quantity b:

where and are the advective and phase velocities normal to the boundary, respectively, and is an inverse relaxation time scale to the reference value . It was first applied to the northern boundary in the Southern Ocean model FRAM (Stevens, 1991). Since then, these methods have been used successfully by, e.g., Döscher and Redler (1997), and Redler and Böning (1997) in the Community Modeling Effort (CME). They have now become an option in the MOM package.

In an explicit formulation, the time-step is determined by the smallest grid spacing. For a model domain that extends to high latitudes, the convergence of coordinate lines in a spherical system causes an increasingly severe constraint on the size of the time-step. A method to relax this limitation is to use Fourier filtering in high latitudes. This method eliminates the shortest wave lengths which are contaminated most by effects of the violation of the CFL criterion.

The semi-implicit treatment of the Coriolis term may be advantageous for coarse resolution simulations which extend into the polar regions. In this case, the time-step is dictated by inertial waves which have periods close to 12 hours near the pole. A longer time-step is possible with implicit treatment of the Coriolis term. Note, however, that this procedure also damps Rossby and topographic waves.

For studies that investigate steady state solutions only, computer time can be saved by so-called asynchronous integration, a term describing the use of different time steps for tracers on one side and momentum on the other (i.e., the artificial acceleration of the tracer advection-diffusion equation relative to the momentum equations). Based on the fact that time scales for evolution of tracer and velocity fields are often an order of magnitude different, it was hypothesized (Bryan, 1984) that a steady state can be reached after different integration periods for different variables. Of course, this assumes that only one steady state exists for a given choice of model parameters and forcing, an assumption that is not always justified.

There is an extensive list of standard diagnostics built into the code. Output of large-scale integral quantities like the global energy cycle, gyre components, meridional overturning and tracer budgets can be chosen. Other tools include the diagnostic computation of sea surface height for rigid lid runs, and the integration of particle trajectories. Plotting and other graphic visualization, however, is left to the user. While plotting is not part of MOM, a NetCDF option provides an interface to a wide variety of utilities for graphic representation of model results.

The model is written in f77 and uses cgs units throughout. Individual options can be specified via cpp options. Versions for parallel computers exist (see, e.g., McClean et al., 1997) but may not be freely available. MOM code and documentation can be obtained under

http://www.gfdl.gov/

on the World Wide Web.

Concluding remarks

Finally, a few words about the efficiency, accuracy, flexibility and handling of MOM. Typical applications have more than 50 percent land points, a large overhead in memory and computations. Consequently, a sophisticated mechanism has been developed to run the model even on machines with limited memory: the computations are performed sequentially for a number of zonal-vertical slabs, which are read in from and written out to physical (or virtual) disk. This general concept is very prone to errors, especially for beginners, who often get confused with the placement of data in storage.

The representation of topography as a series of steps has consequences for several aspects of the ocean circulation, affecting both the wind-driven and thermohaline circulation. For example, the solution for the barotropic mode converges only slowly with increasing horizontal resolution. Also, downslope flow is poorly represented as alternating advective and convective events, with huge artificial mixing. There are a number of recent attempts to overcome the problems of the staircase topography in geopotential coordinate models. One is the implementation of "shaved grid cells" at the bottom, thus introducing slopes in the lower-most grid boxes (Adcroft et al., 1997). Other attempts add some form of a bottom boundary layer sub-model to the code: Dietrich et al. (1987) use a thin shell model to include the effects of boundary-parallel flow into their model. Beckmann and Döscher's (1997) approach redirects advective and diffusive fluxes in the bottom-most layer from horizontal into the boundary-parallel direction. A similar approach was proposed by Campin and Goosse (1999). Killworth and Edwards (1998) employ a variable thickness single-layer concept (much akin to surface mixed layer models) to construct a sub-model that predicts height and circulation within the bottom boundary layer. A constant thickness slab layer model that explicitly includes the bottom slope is being tested by Gnanadesikan (1999). Systematic testing will be necessary to evaluate the applicability and to quantify the improvements by all these schemes.

Other geopotential coordinate models include:

Unlike MOM, these are all on the "C" grid and differ further in numerics, parameterizations and coding.