\section{The MIT ocean-atmosphere model}
\label{s:model}

The MIT ocean-atmosphere model is built around a versatile incompressible
Navier-Stokes solver developed at MIT for the numerical study and analysis
of rotating fluids.  A hydro-dynamical kernel, designed to be sensitive to
the transition from non-hydrostatic to hydrostatic dynamics, lies at the
heart of the model, providing computational efficiency across a broad
spectrum of simulation regimes.  A review of the model's utility for
research in a range of problems related to ocean and climate dynamics can
be found in Marshall et.al.\cite{Marshall97a}.  Implementations of the
model targeted to the MPP CM-5 and to CRAY vector supercomputers are being
actively used by researchers at
MIT\cite{Klinger96,Speer95,Visbeck96,Jones96,Haine96,Menem96,men97}.

The discrete algorithm uses a finite volume\cite{adc96}, 
time-stepping formulation
which is rich in data parallelism and can handle the intricate, order one
variations in geometry of ocean basins.  The CM-5 code achieves sustained
floating point throughput of over 10 GFlop/s on a 512 node CM-5.  Daily
production runs of the model sustain between 500 MFlop/s and 1 GFlop/s on a
32 node CM-5, and between 100 MFlop/s and 200 MFlop/s on a CRAY Y-MP.

\subsection{The numerical scheme}

The numerical procedure used in the model is described in detail in
Marshall et.al.\cite{Marshall97b} The scheme is a variant on the theme set
out by Harlow and Welch\cite{Harlow65}, in which a pressure correction to
the velocity field is used to guarantee non-divergence.

The numerical kernel is a discrete form of the Boussinesq incompressible Navier
Stokes fluid equations: it can be written in semi-discrete form to second 
order in time $\Delta t$ thus:\\

\begin{equation}\label{eq-Vt}
\frac{\underline{v}^{n+1} - \underline{v}^{n}}{\Delta t}=
\underline{G}_{v}^{n+\frac{1}{2}} - \nabla p^{n+\frac{1}{2}} 
\end{equation}

\begin{equation}\label{eq-Continuity}
\nabla .\underline{v}^{n+1} = 0
\end{equation}

Equations \ref{eq-Vt} and \ref{eq-Continuity} describe how the fluid evolves
($\underline{v}=(\underline{v}_h,w)$ is the velocity in the horizontal and
the vertical) in response to forcing due to
$\underline{G}$ (representing inertial, Coriolis, metric, gravitational,
and forcing/dissipation terms) and to the pressure gradient $\nabla p$.
The same kernel algorithm is used by us to study both the 
atmosphere and the ocean - the anelastic approximation as described in
Brugge et.al.\cite{Brugge91} and in Miller\cite{Miller74b} is used
to render the same model suitable for the study of both fluids.

The pressure field required in Equations \ref{eq-Vt} and
\ref{eq-Continuity} is found by separating $p$ into hydrostatic, surface and
non-hydrostatic parts.  Vertical variations in $p$ can be computed
hydrostatically yielding $p_{hy}$ using only information in each column of
the fluid.  Horizontal variations in $p_{s}$ are diagnosed by solving an
elliptic equation which ensures non-divergent depth integrated flow:

\begin{equation}\label{eq-Ps}
\nabla_{h}. H\nabla_{h} {p_{s}} =
\nabla_{h}.\overline{\underline{G}_{v_{h}}^{n+\frac{1}{2}}}^H -
                        \nabla_{h}.\overline{\nabla_{h} {p_{hy}}^\prime}^H
\end{equation}

\noindent
( $\stackrel{\overline{\;\;\;\;\;\;}}{\;\;\;\;\;\;}\stackrel{H}{\;}$
indicates the vertical integral over the local depth, $H$, of the domain
and the subscript $_{h}$ denotes horizontal). In the non-hydrostatic limit
a three-dimensional elliptic equation is inverted for the non-hydrostatic
pressure. 
In this study we consider the hydrostatic case involving the
two-dimensional problem, Equation \ref{eq-Ps}.
More details of the numerical procedure can be found in Marshall
et.al\cite{Marshall97a,Marshall97b}.

\subsection{Program structure}

Figure {\ref{fig:FLOWDIAG}{\it } shows the high-level program structure of
the MIT model. Almost all of the execution time is spent in the
time-stepping loop, which is divided into two parts -- {\bf Forward\_Step}
and {\bf Inversion\_Step}.  The {\bf Forward\_Step} comprises a series of
nearest-neighbor computations involving two time levels of the model state
variables and their time derivatives. This stage computes explicit time
derivative terms for the velocity and tracer equations, steps them forward
in time and updates the density field through an equation of state.  There
is considerable latitude for different orderings of the computations within
{\bf Forward\_Step}, and, as discussed in Shaw et.al.\cite{Arvind95}, this
stage contains an abundance of both data parallelism and functional
parallelism.  The {\bf Inversion\_Step} requires global communication to
solve a set of simultaneous equations to ensure that the flow remains 
non-divergent.

\begin{fig}{fig:FLOWDIAG}{High-level structure of the MIT model. The 
model iterates repeatedly over a time-stepping loop, comprised of two main
blocks, {\bf Forward\_Step} and {\bf Inversion\_Step}.}
\begin{center}
\fbox{
\begin{minipage}{4in}
\sffamily
\begin{tabbing}
\hspace{.15cm} \= \hspace{.15cm} \= \hspace{.15cm} \= \hspace{3cm} \= \kill
 BEGIN \\
 \> Initialize \\
 \> \> Define topography, initialize flow field, tracer distributions\\
 \> FOR each time step DO \\
 \> \> {\sffamily\bfseries Forward\_Step} \\
 \> \> \> Calculate velocity time derivatives  \\
 \> \> \> Calculate tracer field time derivatives \\
 \> \> \> Step forward state variables \\
 \> \> \> Update density field \\
 \> \> {\sffamily\bfseries Inversion\_Step} \\
 \> \> \> Solve elliptic problem to ensure non-divergent flow\\
 \> END FOR \\
 END
\end{tabbing}
\end{minipage}
}
\end{center}
\end{fig}
