Global Ocean Modeling on Massively Parallel Computers

1. Abstract
2. Project Objectives
3. Recent Accomplishments
4. Future Plans
5. References

Abstract

Very high resolution simulations of the ocean currents and distributions of temperature and salinity in the three-dimensional global ocean have been carried out on the 1024-node CM-5 at Los Alamos. With 1280 points in longitude (0.28deg.) and 896 variably spaced points in latitude (0.17deg.), the model has a horizontal resolution of 31 km at the equator decreasing to 7 km at 78deg. latitude. Mesoscale ocean eddies and narrow western boundary currents like the Gulf Stream are well resolved. Thanks to the improved model formulation and numerical methods, it is possible to integrate one year of simulated time at this very high spatial resolution in only 28 hours on 256 nodes of the CM-5, including output of about 50 gigabytes of data every simulated year. Comparisons between simulated and observed data are very encouraging.

This work is sponsored by the DOE CHAMMP and HPCC programs.

Project Objectives

Recent Accomplishments

(1) The Semtner-Chervin Model on the Connection Machine: We began with a version of the Bryan-Cox-Semtner ocean model developed by Albert Semtner at the Naval Postgraduate School in Monterey, CA, and Robert Chervin at the National Center for Atmospheric Research in Boulder, CO [Semtner and Chervin, 1988]. That code is highly vectorized and multi-tasked to run on multi-processor Cray XMP and YMP computers. First, the model was rewritten in data-parallel FORTRAN 90 for the CM-2 Connection Machine, and although it was completely reorganized to accommodate a new data structure more appropriate for parallel architectures, all of the basic numerical algorithms were left unchanged. We found that the performance of this original version (for a high-resolution problem: 1/2deg. x 1/2deg. with 20 vertical levels) on the full 64K CM-2 (with 2048 FPUs) was approximately equivalent to the performance of the Semtner-Chervin code on a 4-processor XMP. However, it was apparent that certain parts of the code performed poorly on the CM because the numerical algorithms did not parallelize well. For historical reasons, dating back to the original development of the model by Kirk Bryan in the late 1960's, the equations for the barotropic mode, which describe the vertically averaged flow, were expressed in terms of a volume-transport streamfunction. The boundary conditions on the streamfunction involve line integrals around islands, which lead to computational inefficiencies on parallel machines since they require nonlocal communication. In the original CM-2 version of the code, the barotropic solver consumed ~75% of the total run time. The line integrals also lead to inefficient vectorization on vector machines, and for this reason Semtner and Chervin were forced to simplify the topography to include only three islands in their 1/2deg. degree global ocean simulations.

(2) The Surface Pressure Formulation: To overcome this problem we reformulated the model by changing from the streamfunction (SF) to a surface-pressure (SP) formulation of the barotropic equations. The SP method involves local boundary conditions at coastal points, and hence does not require any line integrals around islands. This method has several other advantages over the SF method: (a) any number of islands can be included in the computational grid without a performance penalty; (b) the model can handle steep gradients in the bottom topography (unlike the SF formulation where the depth field has to be smoothed in order to prevent numerical instabilities); and (c) the surface pressure is now computed directly, and since surface pressure is approximately proportional to surface elevation above mean sea-level, this simplifies assimilation of altimetric data from satellite observations into the model.

Both SF and SP formulations are based on the rigid-lid approximation, which eliminates fast surface waves but leads to an elliptic Poisson-like equation for the streamfunction or surface pressure field. The Semtner-Chervin model uses a simple Jacobi relaxation method to solve this elliptic equation. We replaced this with a much more powerful conjugate gradient (CG) solver. However, in order to use a standard CG algorithm, which requires a symmetric operator matrix, we used an approximate factorization technique [Dukowicz and Dvinsky, 1992] to split off the Coriolis terms which are antisymmetric. This splitting maintains the second-order accuracy of the time discretization scheme, but leaves a symmetric operator to which the standard CG algorithm can be applied. We also developed a new parallelizable preconditioning method in order to accelerate convergence in the solver. Together, these improvements produced better than a factor of two improvement in run time over the original method on the CM-2 in a 1/2deg. global model, even though the mesh for the SP method was unsmoothed and contained 80 islands, while the mesh for the SF method had to be smoothed and contained only 3 islands. These changes to the model, including the SP formulation, the improvements in the elliptic solver, and the implementation and performance of the model on the CM-2, have been documented in two publications [Smith et al., 1992; Dukowicz et al., 1993].

(3) The Implicit Free-Surface Formulation: The SP formulation substantially reduces the time spent solving the barotropic equations (especially if many islands are included in the computational grid), but in spite of the improvements in the elliptic solver it still takes typically half of the time spent in the full code. Furthermore, it was found that the SP formulation required the use of a 9-point operator in the elliptic equation for reasons of energetic consistency (unlike the SF formulation where a 5-point operator could be used). The disadvantage of the 9-point operator is that it contains a "checkerboard" null space which can allow unwanted checkerboard patterns to appear in the solution for the surface pressure field. In the course of investigating this phenomenon, we realized that we could both eliminate the checkerboard null space and greatly reduce the time spent in the elliptic solver if we abandoned the rigid-lid approximation in favor of a free surface. The fast surface waves may be handled instead using an implicit time discretization scheme which leads to an elliptic operator for the surface pressure field that is virtually identical to the operator in the rigid-lid SP formulation, except that it is much better conditioned, and results in much more rapid convergence in the solver.

The implicit free-surface (FS) method retains the advantages of the SP formulation, and in addition to alleviating checkerboarding and substantially improving the computational efficiency, it has the following advantages: (a) it improves accuracy in the computation of long-wavelength Rossby waves; (b) it provides robustness by selectively damping computational modes; and (c) the sea-surface elevation is exactly proportional to the surface pressure (instead of only approximately, in the limit of long time scales, as in the rigid-lid formulation). The FS formulation is computationally much more efficient than either of the rigid-lid formulations. The solution of the barotropic mode with the FS method is an order of magnitude faster than the SP method if the same convergence criterion is used, and takes up less than 25% of the time spent in the full code on the CM. Of course, these new numerical methods should also improve performance of the Cray codes. We have implemented the implicit free-surface method as an option in both the Semtner-Chervin code and the widely-distributed Modular Ocean Model (MOM) developed at the Geophysical Fluid Dynamics laboratory (GFDL) in Princeton, NJ. The formulation, implementation, and testing of the implicit free-surface method is documented in Dukowicz and Smith [1994].

(4) Other Numerical Improvements: Another improvement was the implementation of a well-known technique for doubling the stability limit associated with baroclinic gravity waves [Brown and Campana, 1978], which involves time-averaging the hydrostatic pressure gradient. This method has apparently been tried in other versions of Bryan-Cox models without much success and is generally not used. However, it did work in our model, and with it we have achieved a factor of two increase in the time step. We eliminated other factors, besides gravity waves, that might be controlling the time step by performing a detailed linear stability analysis of the advection and mixing terms. Furthermore, our new barotropic formulations do not suffer from the topographic instabilities which limit the time step in the streamfunction formulation [Killworth, 1987]. Thus we determined that gravity waves alone control our time step, and this allows the pressure-averaging technique to be effective. As an example, using this technique in combination with the other improvements discussed above, our 1/2deg. 20-level global model can be integrated with a 90-minute time step, which compares to a time step of 15 minutes used by Semtner and Chervin in their 1/2deg. runs.

(5) New Horizontal Grid: The original model uses the standard polar longitude-latitude grid which has a singularity at the poles, making it difficult to include the Arctic basin. The customary way to avoid the severe time step restriction associated with thin cell-widths at the pole is to apply ad-hoc spatial filtering to the prognostic fields. On parallel machines such filtering is expensive, because it involves long-distance communication, and it results in a load-balance problem, since the filter is applied only at high latitudes. We have extended the model to allow for general orthogonal coordinates for the horizontal grid, which allows us to include the Arctic using a grid with a displaced pole. This maintains finite cell widths everywhere and makes filtering unnecessary. Long runs at ~1deg. resolution have been made with this grid.

(6) The Parallel Ocean Program (POP): Our ocean model, originally written for the CM-2 and CM-200 Connection Machine systems, has since been completely rewritten to make it more modular and easily portable to other machine architectures, including both data-parallel and message-passing parallel architectures, workstation clusters, and more conventional vector or serial computers. Since we wanted to retain the simplicity and compactness of array-syntax, the source code is written in Fortran 90, but is set up to be easily translated into Fortran 77 if the target machine does not support this. To avoid problems that might arise to due to differences in Fortran 90 compilers, and to simplify the translation to Fortran 77, only the simplest Fortran 90 features have been used in the code, namely, array syntax and masked array operations (WHERE statement).

The main portability issues have to do with the parallel implementations: how is the data distributed across processors, and how is the inter-processor communication handled, especially given the variety of different message-passing protocols? To deal with these problems, POP is structured so that all communication occurs in a few low-level subroutines, called stencil and global-reduction routines. All the finite-difference operators in the code can be reduced to a handful of local stencils (e.g., a 5-point stencil involving a gridpoint and its nearest-neighbors to the north, south, east, and west). The only other communication is in the form of global SUM, MAX, and MIN routines (used in the conjugate-gradient solver and for computing various diagnostic quantities). The idea is that these few routines may be rewritten to optimize communication on a particular architecture, while the remainder of the code is machine-independent and easily portable. To date, POP has been successfully ported and run on many different systems, including: Vector / Serial: Cray YMP and C90, Sun workstation; Data Parallel: CM-200, CM-5, MasPar; Message-Passing: CM-5, Intel Paragon, Cray T3D, IBM workstation cluster.

(7) High Resolution Global Ocean Simulations: The culmination of our three-year CHAMMP project is a set of very high resolution global ocean simulations carried out with POP on the CM-5. The model uses a Mercator grid covering the global ocean from 78S to 78N with 1280 grid points in longitude, 896 points in latitude, and 20 vertical levels. The resolution varies from 31 km at the equator to 6.5 km at the highest latitudes, with an average latitudinal resolution of 1/6deg., sufficient to resolve mesoscale eddies and western boundary currents quite well. This work is being done in collaboration with Albert Semtner, who provided us with the grid topography, the initial temperature and salinity fields from his 1/4deg.-average run, and with surface temperature, salinity, and wind stress fields to drive the model, all interpolated to the model grid. The buoyancy forcing consists of restoring to monthly Levitus climatology with a one-month timescale at the surface and with a 3-year timescale at depth (< 2 km) in selected regions (latitudes >70N, <70S, and in the Mediterranean outflow region). A 9-year integration is currently underway using daily ECMWF winds covering the period from January 1985 to December 1993.

With the improved numerics described above, we can integrate the model with a 30-minute time step at this resolution. The run time, including output of 50 GB per simulated year, is 18 hours per simulated year on one-half of the Los Alamos 1024-node CM-5 or 28 hours per year on one-quarter of the machine. With these run times, multi-decade simulations are feasible. While this is not nearly long enough to equilibrate the deep ocean (which would require >1000 year simulations), it is long enough to study important features of the upper ocean circulation such as narrow western-boundary currents and the effects of mesoscale eddies.

Future Plans

It is now well established that the oceans play a fundamental role in climate variability on these long timescales. In order to accurately model such variability it is necessary to run simulations for hundreds or even thousands of years: the time required for the deep ocean to equilibrate to changes in surface forcing. Eddy-resolving models with spatial resolutions of order 10 to 30 km are not yet capable of such a task. For example, integrating the 1/6deg. model discussed below for 1000 years would take more than two full years of cpu time on 512 nodes of the CM-5. One possibility for overcoming this problem is to develop a fully implicit ocean GCM capable of taking much longer time steps, which could then be efficiently integrated for thousands of years. Another possibility is to find an accurate way of parameterizing the effects of the eddies in coarser resolution models. Other key small-scale processes such as convective overturning are unlikely ever to be resolved in direct numerical simulations of global ocean circulation, and improved subgrid-scale parameterizations of these processes also need to be developed and tested.

Given these considerations, our research over the next few years will concentrate in four areas: (A) continuation of our ongoing series of high resolution (~1/6deg. and possibly higher), multi-decade global ocean simulations to study the role of mesoscale eddies and other features of the upper ocean circulation with further refinements in surface forcing; (B) the implementation and testing of parameterizations of subgrid-scale phenomena and surface forcing that are relevant to climate studies, involving relatively long (100-200 year) simulations at moderate resolution (~1deg.) which will also be used to investigate inter-decadal variability; (C) the development and implementation of a sea-ice model compatible with the numerics and general-coordinate grid capability of POP; and (D) the development of an implicit 3-D primitive equation model to provide the capability for very long (>1000 year) integrations at moderate resolution.

(A) High Resolution Simulations: We will conduct additional high-resolution, multi-decade experiments on the Los Alamos CM-5, like the 1/6deg. simulation currently underway, with improved representations of surface fluxes, subgrid-scale parameterizations, and possibly also increased resolution. The model output will be analyzed and compared in detail with available experimental observations, in particular: satellite altimetry data from Geosat and Topex/Poseidon which will provide information on the distribution and magnitude of surface-height variability and eddy kinetic energy; and data on temperature, salinity, and currents being gathered by the World Ocean Circulation Experiment (WOCE). Our work will complement and support the WOCE program, whose underlying scientific goal is "to improve global ocean models to advance a climate prediction capability extending to decadal timescales" [U.S. WOCE Implementation Plan, 1992].

We will carry out extensive analyses of these 1/6deg. runs, looking at quantities such as transports of the major currents; time-mean and eddy transports of momentum, heat and salt; seasonal and annual heat transports in the various ocean basins; barotropic and meridional transport streamfunctions; and surface height and slope variability.

The results of this analysis will undoubtedly provide insight into aspects of the model physics or forcing that should be improved. This will aid in the design of further simulations at 1/6deg.. We already have several improvements planned. In the simulation now underway, we have switched from monthly-mean winds to daily winds from the same ECMWF analysis. With monthly winds, the model gives fairly realistic surface-height variability in regions containing strong currents, but the variability predicted in quiescent regions of the ocean is weaker than expected from satellite measurements. In simulations at 1/2deg. resolution, Chervin [1994] has recently shown that the use of daily winds produces a more realistic level of this "background" variability.

We also intend to test improved treatments of surface buoyancy forcing in the next series of high resolution experiments. In particular, we plan to replace the current restoring boundary condition on temperature with a Haney-Han type flux boundary condition, using results from the recent analysis of Barnier et al. [1994], in which heat fluxes were estimated from a 3-year seasonal climatology of ECMWF surface atmospheric fields. In this formulation, the heat flux includes a climatological term plus a correction term which varies with model SST and is designed to account for the feedback from the ocean to an atmosphere that remains in a climatological state.

To investigate the issue of model convergence as a function of resolution, the model output will be compared with Semtner and Chervin's 1/2deg. and 1/4deg.-average runs, and detailed comparisons with satellite altimetry data will be made to estimate the minimum resolution required to produce the correct spectrum of mesoscale eddies. Depending on the results of these analyses, it may be necessary to design experiments with higher horizontal and/or vertical resolution. At ~1/10deg., a 3-5 year integration is feasible on the CM-5.

(B) Testing Model Physics Parameterizations: Ocean simulations designed to study the long-term variability of the global thermohaline circulation require an accurate representation of the processes which carry heat, salt and passive tracers from the surface to deeper levels. Deep water formation occurs primarily by one of two types of processes: deep sinking in polar regions due to cooling and/or increased salinity at the surface leading to unstable stratification and convective overturning; and intermediate water formation due to convection and wind-induced downwelling at mid-latitudes involving in mixing and transport along isopycnal surfaces. The first type of process occurs on very small scales, and cannot be directly simulated in any global ocean model, so it must be parameterized. In the second process, mesoscale eddies play an important role, because they mix water properties along isopycnals and because eddy generation by baroclinic instability acts to smooth out lateral gradients in stratification. These effects must also be accurately parameterized, at least in the non-eddy-resolving models that can be used in long-term climate studies. We plan to implement some recently developed parameterizations of these processes in POP, and to test them by performing relatively long (~100 year) integrations at moderate (~2/3deg.) resolution. Integrations at other resolutions, including the high-resolution runs discussed above, may also be carried out to test the scalability of the parameterizations. We will also use the output of the high-resolution simulations to test the validity of eddy parameterizations used at lower resolution.

(1) Parameterization of Mesoscale Eddies: Recently, Gent and McWilliams have proposed a new parameterization designed to represent the transport of tracers by mesoscale eddies in non-eddy-resolving models [Gent and McWilliams,1990; Gent et al., 1994]. In Eulerian models such as Bryan-Cox-Semtner, this parameterization has two components. The first is a rotation of the lateral diffusion terms so that they only mix tracers along surfaces of constant potential density (isopycnals); this replaces the physically unjustifiable horizontal mixing terms commonly used. The second is to add a term to the flux velocities for tracer advection to account for the transport due to unresolved eddies. This "eddy-induced transport velocity" is modeled as adiabatic down-gradient diffusion of the thickness between neighboring isopycnal surfaces. It therefore acts to smooth out lateral gradients in stratification and thereby mimics baroclinic instability. Recent studies at low resolution [Danabasoglu et al. 1994, Boning et al. 1994] show that this parameterization appears to dramatically improve features of the ocean circulation and eliminate problems that have traditionally plagued the standard formulation with horizontal mixing: it produces a sharper, more realistic main thermocline and a cooler abyssal ocean; it eliminates unphysical upwelling in the Gulf stream region (the so-called "Veronis effect") resulting in a meridionally-expanded overturning streamfunction in the North Atlantic that strengthens the poleward and surface heat fluxes; and it drastically reduces the number of convective-adjustment events, confining them to regions where deep convection is known to occur. We plan to collaborate with Gent and McWilliams to implement and test their parameterization in POP. They will work with us on the design and analysis of simulations at moderate and high resolution, which will focus on effects of the parameterization on the upper-ocean circulation that can be investigated with relatively short-duration integrations.

(2) Parameterization of Deep Convection: It is believed that, at least in the North Atlantic, deep water formation due to unstable stratification takes place in convective "chimneys" tens of kilometers in diameter, within which many small-scale plumes develop that carry surface water to depth and in the process entrain surrounding water. Depending on conditions, the effect of these plumes may be to simply homogenize the water column, or it may actually cause an exchange of surface and deep water properties. Current models typically parameterize such deep water formation by vertical mixing of water properties in unstable regions, using either convective adjustment or an implicit vertical mixing scheme. These are most likely inadequate for modeling convective overturning under the broad range of conditions present in the real ocean. We plan to implement in POP the convective parameterization developed by Paluszkiewicz, Denbo, and Skyllingstad (PNL). This model was designed to parameterize the effects of plume convection and was validated by comparison to direct numerical simulation of plumes using a Large Eddy Simulation model.

(C) Development of a Sea-Ice Model for POP: A sea-ice model is an essential component of a climate modeling system, because sea ice can significantly affect the large scale transfer of energy between the atmosphere and the oceans. For instance, the formation of ice reduces the transfer of heat between the ocean and the atmosphere, thereby inhibiting the normal moderating effect of the ocean. Furthermore, its high surface albedo causes a large fraction of the solar energy to be reflected back to space. Under certain conditions these effects might result in a positive feedback causing an even colder atmosphere and further ice growth. In addition, the melting and freezing of sea ice affects the sea-surface salt flux and thereby provides a primary driving source for the thermohaline circulation. These effects are crucial for properly modeling the climate.

We will formulate, develop, and incorporate an improved sea-ice model in POP. The stress will be on improving existing models and on developing discretization methods that are suitable for massively parallel computers and compatible with the general-coordinate capability of POP. Because of this capability, POP is ready to accept a sea-ice model without the problem of a mesh singularity in the Arctic, and without the need for filtering at high latitudes which is common in other codes. This work will be done in collaboration with Drs. Steve Piacsek and Ruth Preller of Naval Research Laboratory / Stennis Space Center (NRL/SSC), both of whom are experienced in sea-ice modeling and are familiar with the extensive sea-ice data available at NRL/SSC for model validation.

(D) Development of an Implicit 3D Primitive Equation Model: In the last several years it has become increasingly evident, based on both GCM simulations and observational evidence, that the thermohaline circulation may undergo natural internal variability on timescales ranging from decades to millennia [Weaver and Hughes, 1992]. In order to study the thermohaline circulation, its bifurcations to multiple solutions, and the long-term response of the ocean to external forcing, routine integrations for ~1000 years and longer will be necessary. As discussed above, POP is currently not a suitable tool for such long-term computations. Its explicit treatment of the baroclinic modes means that it is constrained by internal gravity-wave stability limits, which in practice implies a time step limit of 3 to 4 hours at ~1deg. resolution, and 1/2 to 1 hour at ~1/6deg..

Long-term ocean simulations typically employ very coarse resolutions (3 - 4deg.) for reasons of computational economy. However, it is not clear that such coarse resolutions are sufficiently accurate, because important spatial scales such as the width of the western boundary currents are not resolved. Furthermore, although there is some promise, it is not yet clear that turbulence parameterizations exist that apply at such coarse resolutions. We would therefore like to be able to do such calculations at finer resolutions (e.g., < 1deg.). This implies, however, that time steps much longer than those used presently are needed for the computation to remain affordable, and that new solution techniques are needed to overstep the current stability limits, which are set by internal gravity time scales (minutes to hours), the inertial time scale (two hours), as well as advective time scales (days). In general we would like to eliminate all these limits, but, depending on the resolution, the elimination of just the gravity and inertial limits will produce significant gains.

A method is being developed that we expect will overcome external and internal gravity wave and inertial stability limits. If advective stability limits can also be successfully overcome, then the remaining stability limit is that associated with Rossby waves, which is of the order of several days for the barotropic mode and order one month for the baroclinic modes. This is a great improvement over the present capabilities of POP, whose practical limit is currently set by the inertial limit of about two hours even at coarse resolution, and it therefore should permit economical long-term integrations.

References

Boning, C. W., F. O. Bryan, and W. R. Holland, 1994: An Overlooked Problem in Model Simulations of the Thermohaline Circulation and Heat Transport in the Atlantic Ocean, submitted to J. Climate.

Brown, J. A. Jr., and K. A. Campana, 1978: An Economical Time-Differencing System for Numerical Weather Prediction, Mon. Wea. Rev. 106, 1125-1136.

Chervin, R. M., 1994: Dependence of Global Ocean Circulation on Wind Stress Forcing, Fifth Symposium on Global Change Studies, Nashville, Tennessee, January 23-28, 1994.

Danabasoglu, G., J. C. McWilliams, and P. R. Gent, 1994: The Role of Mesoscale Tracer Transports in the Global Ocean Circulation, Science, 264, 1123.

Dukowicz, J. K. and A. S. Dvinsky, 1992: Approximate Factorization as a High Order Splitting for the Implicit Incompressible Flow Equations, J. Comput. Phys., 102, 336-347.

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. Atmospheric & Oceanic Tech., 10, 195-208.

Dukowicz, J. K. and R. D. Smith, 1994: Implicit Free-Surface Method for the Bryan-Cox-Semtner Ocean Model, J. Geophys. Res., 99, 7991-8014.

Gent, P. R., and J. C. McWilliams, 1990: Isopycnal Mixing in Ocean Circulation Models, J. Phys. Oceanogr., 20, 150-155.

Gent, P. R., J. Willebrand, T. J. McDougal, and J. C. McWilliams, 1994: Parameterizing Eddy-Induced Tracer Transports in Ocean Circulation Models, submitted to J. Phys. Oceanogr.

Killworth, P. D., 1987: Topographical Instabilities in Level Model OGCMs, Ocean Modeling.

Semtner, A. J. Jr. and R. M. Chervin, 1988: A simulation of the Global Ocean Circulation with Resolved Eddies, J. Geophys. Res., 93, 15502-15522.

Semtner, A. J. Jr. and R. M. Chervin, 1992: Ocean General Circulation from a Global Eddy-Resolving Model, J. Geophys. Res., 97, 5493-5550.

Smith, R. D., J. K. Dukowicz and R. C. Malone, 1992: Parallel Ocean General Circulation Modeling, Physica D., 60, 38-61.

Weaver, A., J. and T. M. C. Hughes, 1992: Stability and Variability of the Thermohaline Circulation and its Link to Climate; to appear in Trends in Physical Oceanography , Council of Scientific Research Integrations, Research Trends Series, Trivandrum, India, April, 1992.