A Parallel Navier-Stokes Ocean Model: Formulation and Documentation
Daniel Jamous, Chris Hill, Alistair Adcroft, John Marshall
December 1997
Massachusetts Institute of Technology
CONTENTS...................................................................................................................... 3
INTRODUCTION.............................................................................................................. 4
Chapter 1 Continuous
formulation........................................................... 6
1.1 The continuous
equations of the model ocean..................................................................................................... 6
1.2 Boundary conditions................................................................................................................................................ 7
1.3 Coordinate systems.................................................................................................................................................. 9
1.4 Hydrostatic,
non-hydrostatic and quasi-hydrostatic forms............................................................................. 11
1.5 Finding the
pressure field.................................................................................................................................... 12
1.5.1 Depth-averaged
pressure................................................................................................................................ 12
1.5.2 Surface,
hydrostatic and non-hydrostatic pressure.................................................................................... 13
1.6 Summary of solution
strategies.......................................................................................................................... 15
1.7 Forcing and
dissipation......................................................................................................................................... 16
1.7.1 Forcing................................................................................................................................................................ 16
1.7.2 Dissipation......................................................................................................................................................... 17
CHAPTER 2 Discrete
formulation............................................................... 19
2.1 Temporal
discretization........................................................................................................................................ 19
2.2 Spatial
discretization............................................................................................................................................. 21
2.2.1 Volume
quantities............................................................................................................................................. 22
2.2.2 Face
quantities.................................................................................................................................................. 24
2.2.3 Coriolis terms..................................................................................................................................................... 25
2.2.4 Pressure
gradient force.................................................................................................................................... 25
2.2.5 Forcing and
dissipation................................................................................................................................... 25
2.2.6 The C-D scheme................................................................................................................................................ 26
2.2.7 Conservation
properties.................................................................................................................................. 26
2.3 The Elliptic problem.............................................................................................................................................. 27
2.3.1 Discrete
formulation......................................................................................................................................... 27
2.3.2 Preconditioned
conjugate-gradient solution method.................................................................................. 30
CHAPTER 3 PARALLEL
IMPLEMENTATION........................................................ 36
2.1 Enumeration of
arrays........................................................................................................................................... 36
Chapter 4 Procedures and variables.................................................. 38
4.1 Procedures.............................................................................................................................................................. 38
4.1.1 Flow chart........................................................................................................................................................... 38
4.1.2 Alphabetical
list of procedures and 'include' files....................................................................................... 40
4.2 Variables.................................................................................................................................................................. 40
Chapter 5 Getting
started: practical considerations............. 41
5.1 Accessing the code
- directory structure.......................................................................................................... 41
5.2 Configuring the
model and setting up a run..................................................................................................... 41
5.2.1 domain,
geometry.............................................................................................................................................. 42
5.2.2 Equation of
state............................................................................................................................................... 42
5.2.3 Hydrostatic,
quasi-hydrostatic, or non-hydrostatic regime....................................................................... 43
5.2.4 Momentum
equations...................................................................................................................................... 43
5.2.5 Diagnosis of
pressure...................................................................................................................................... 44
5.2.6 Tracer
equations (temperature and salinity)................................................................................................. 45
5.2.7 Run setup........................................................................................................................................................... 46
5.3 Compiling, running,
and restarting the code.................................................................................................... 47
5.3.1 compilation......................................................................................................................................................... 47
5.3.2 execution............................................................................................................................................................ 47
5.3.3 output and
restart............................................................................................................................................. 48
5.4 Examples :
description of test cases.................................................................................................................... 49
5.4.1 Two-layer
ocean in a rectangular basin driven by a steady wind-stress................................................. 49
5.4.2 Neutral ocean
in a doubly periodic domain driven by buoyancy forcing................................................ 50
5.4.3 Exponentially
stratified ocean on a wind-driven channel........................................................................... 50
5.4.4 Global ocean
driven by realistic wind and buoyancy forcing.................................................................... 50
5.4.5 Model
configuration for the test cases......................................................................................................... 51
APPENDICES................................................................................................................ 53
References.............................................................................................................. 54
This manual aims to provide the reader with the
information necessary to carry out numerical experiments using the M.I.T GCM.
The manual provides a comprehensive description of the numerical paradigm
employed by the M.I.T. GCM and a description of the associated program code.
Examples showing the use of the code in small process studies and for
basin/planetary scale simulations are presented.
Any
practical numerical rendition involves a series of compromises. Our guiding principle has been to devise
methods which, as far as is possible, are competitive across a large range of
scales from
, the depth of the fluid, right up to horizontal resolutions
, coarser than the Rossby radius of deformation. We did not
want to make the hydrostatic assumption a-priori,
since it precludes the study of many interesting small-scale phenomenon. Rather it was important to us that our
approach be also well-suited to the convective, non-hydrostatic limit. We therefore adopt pressure (or height) as a
vertical coordinate and employ a ‘finite volume’ approach, in which property
fluxes are defined normal to the faces that define the volumes, leading to a
very natural and robust discrete analogue of ‘divergence’. In the special case that the volumes are of
regular shape, the arrangement of the model variables in the horizontal reduces
to a ‘C’ grid, using the nomenclature of Arakawa and Lamb (1977), and so
carries with it well-documented strengths and weaknesses in the treatment of
gravity, inertial and Rossby wave modes (for a recent discussion see Dukowicz,
1995).
The
basis for the numerical procedure employed in the M.I.T GCM is the
incompressible Navier-Stokes equations. The model integrates forward the
Boussinesq approximation to the full Navier-Stokes equation
governing fluid motion. The model’s numerical integration method is arrived at
by utilising the “pressure method” [Harlow and Welch, 1965; Williams, 1969;
Potter, 1976]. In this scheme the continuity relation 
is used to give an
elliptic equation for the pressure field. Starting from these basis equations
the model is formulated to allow three different solution strategies. In the
non-hydrostatic - NH - mode of
integration, the basis equations are stepped forward without further
approximation - in particular full non-hydrostatic dynamics are retained. The
hydrostatic primitive equations - HPE
- mode of the model integrates an approximated form of the full equations,
which assumes exact hydrostatic balance. The quasi-hydrostatic - QH - mode of integration admits
additional Coriolis terms, which augment HPE
with an analytically exact angular momentum principle.
The
‘kernel’ algorithm solves the incompressible Navier Stokes equations on the
sphere and can be employed as an atmospheric model or as an ocean model in a
geometry as complicated as that of the ocean basins with irregular coastlines
and islands. The algorithm builds on ideas developed in the computational fluid
community. The numerical challenge is
to ensure that the evolving velocity field remains non-divergent; a ‘pressure
correction’ to the velocity field is used to guarantee non-divergence. The correction step is equivalent to, and is
solved as, a Poisson equation for the pressure field with Neumann boundary
conditions. This Poisson inversion is
the most computationally demanding part of the algorithm and, because it is not
localised in space, presents the biggest challenge in mapping on to parallel
computers since it demands ‘communication’ across the grid to the boundary, and
hence between processors. We have
devised a method for solving this Poisson equation which exploits knowledge of
the dynamics (by separating the pressure field in to hydrostatic,
non-hydrostatic and surface pressure components), the geometry of the ocean
basin (they are much wider than they are deep) and naturally lends itself to
data-parallel implementation. The model
will, of course, also run on a single processor.
Much of
what is here has been already described in the two reference papers on the
model published in the 'Journal of
Geophysical Research' [Marshall et al., 1997a,b]. The manual is organized
as follows: in Chapter 1, we describe the continuous formulation of the model
as well as the different solution strategies employed to find the pressure
field for the HPE, QH, and NH modes. This chapter corresponds essentially to sections 3 and 4
of Marshall et al., 1997a. Chapter 2 describes the temporal and spatial
discretizations used in the model and the numerical approach used to solve the
Elliptic problem for the pressure field. This chapter corresponds essentially to
sections 3 and 4 of Marshall et al., 1997b. Chapter 3 is not written yet and is
meant to describe the parallel implementation of the model code in Fortran 90
and High Performance Fortran (HPF). A flow chart and an alphabetical list of
the procedures and variables are given in Chapter 4. In Chapter 5, the most
practical section of the manual, we review what you need to know to run an
experiment with the M.I.T. GCM (accessing the code, configuration, compilation,
description of output files etc…). We provide some examples that have been
published in the literature, ranging from the convective scale up to the
planetary/global scale. These examples should help the user to get more
familiar with the model.
It would be nice to include an history of the
model in the introduction!
The
state of the ocean at any time is characterized by the distribution of currents
v, potential temperature T, salinity S, pressure p and density
r. The equations that govern the evolution of these
fields, obtained by applying the laws of classical mechanics and thermodynamics
to a Boussinesq incompressible fluid are, using height as the vertical
coordinate:
(1.1.1)
(1.1.2)
(1.1.3)
(1.1.4)
(1.1.5)
where
(1.1.6)
is the
velocity in the lateral and vertical direction respectively. In (1.1.1), 
is the pressure-gradient force, where
(1.1.7)
is the
dynamic pressure expressed in terms of dp, the deviation of the pressure from that of a
resting, hydrostatically balanced ocean
(1.1.8)
where,
is a
constant reference density. 
represents the advective, Coriolis, buoyancy, forcing and
dissipative terms in the momentum equations, i.e. in vector form:
(1.1.9)
In
(1.1.9), W is the angular speed of rotation of the Earth, g is gravity and Fv and Dv
are the forcing and dissipative terms, respectively. The explicit expression of
Gv, component by
component, is given in spherical coordinates in section 1.3.
The
tendencies of temperature and salt are given by
(1.1.10)
(1.1.91)
where 
, 
, 
, 
represent the forcing and dissipative terms in the
temperature and salt equations respectively.
Unlike
the prognostic variables u, v, w,
T and S, the pressure field must be obtained diagnostically. Taking the
divergence of (1.1.1) and using the continuity equation (1.1.2), leads to a
three-dimensional Elliptic equation for the pressure:
(1.1.12)
For a
given field of F, (1.1.12) must be inverted for p subject to appropriate choice of boundary conditions. This method is usually called The Pressure Method [Harlow and Welch, 1965; Williams, 1969;
Potter, 1976]. In the hydrostatic primitive equations case (HPE), the 3-d problem, (1.1.12),
degenerates in to a 2-d inversion:- see section 1.5.
The
configuration of the ocean basin is defined by its depth H(l,f) and allows arbitrary specification of the coastline,
bottom topography and connectedness - see Fig.1. We apply the condition of no normal flow through all solid
boundaries - the coasts and the bottom:
(1.2.1)
where n is a vector of unit length normal to
the boundary. The kinematic condition (1.2.1) is also applied to the vertical
velocity at the surface of the ocean in both free-surface and rigid-lid models.
No-slip (
) or slip (
) conditions are employed on the tangential component of
velocity 
at all solid boundaries, depending on the form chosen for the
dissipative terms in the momentum equations - see section 1.7.
Eq.(1.2.1)
implies, making use of (1.1.1), that:
(1.2.2)
presenting
inhomogeneous Neumann boundary conditions to the Elliptic problem
(1.1.12). As shown, for example, by Williams (1969), one can exploit classical
3D potential theory and, by introducing an appropriately chosen d-function sheet of ‘source-charge’, replace the
inhomogenous boundary condition on pressure by a homogeneous one. The source term
F in (1.1.12) is the divergence of the vector 
. By simultaneously
setting 
and 
on the boundary the following self-consistent but simpler homogenised
Elliptic problem is obtained:
(1.2.3)
where 
is a modified 
such that 
. As is implied by (1.2.2) the modified
boundary condition becomes:
(1.2.4)
At the
ocean bottom and side the diffusive flux of heat and salt is set to zero:
(1.2.5)
where 
is a ‘diffusion’
coefficient normal to the boundary.
Include Fig.1 here: schematic diagram of an ocean
basin showing the irregular geometry: coastlines and islands.
In this
section, the explicit expression of 
in spherical coordinates is given. u, v,
w the velocity components in the
zonal, meridional and vertical direction respectively, are given by (see Fig.2)
:
(1.3.1)
Here f is the latitude, l the longitude, r the radial distance of the particle from the center of the earth,
W is the angular speed of rotation of the Earth and
D/Dt is the total derivative.
Include
Fig.2 here
Fig.2. The spherical polar velocities (u,v,w), the latitude is f and the longitude l.
are inertial,
Coriolis, metric, gravitational, and forcing/ dissipation terms in the zonal,
meridional and vertical directions defined by:
(1.3.2)
(1.3.3)
(1.3.4)
where 
are non-conservative
forces acting on the fluid. The precise
forms of the 
’s which are available in the model are described in section
1.7.
The
tendency of the temperature and salt are given by:
(1.3.5)
(1.3.6)
where 
, 
, 
, 
are the forcing/dissipation terms in the T and S
equations respectively; they are described in section 1.7.
In the
above the ‘grad’ (Ñ) and ‘div’ (Ñ.) operators are defined by, in spherical
coordinates:
(1.3.7)
(1.3.8)
A
consistent Cartesian version of the model can also readily be posed - it does
not carry metric terms and cosf Coriolis terms - see Appendix.
In the
non-hydrostatic model all terms in equations (1.3.2 ® 1.3.4) are retained. A three dimensional elliptic equation (1.2.3) must be solved with
boundary conditions (1.2.4). It is
important to note that use of the full NH
does not admit any new ‘fast’ waves in to the system - the incompressible
condition (1.1.2) has already filtered out acoustic modes. It does, however, ensure that the gravity
waves are treated accurately with an exact dispersion relation.
NH has the following energy equation:
(1.4.1)
where 
is the
three-dimensional velocity vector, Q
is the buoyancy forcing, F
represents the forcing/dissipation terms in the momentum equations, and 
. Note that the
pressure work term 
vanishes when integrated over the ocean basin if
all bounding surfaces, including the upper one, are assumed rigid. NH has a complete angular momentum
principle - see White and Bromley, 1995; Marshall et.al.1997a.
In HPE all the underlined terms in Eqs.
(1.3.2 ® 1.3.4) are neglected and ‘r’ is replaced by ‘a’,
the mean radius of the earth. The 3-d
Elliptic problem reduces to a 2-dimensional one since once the pressure is
known at one level (we choose this level to be the surface) then it can be
computed at all other levels from the hydrostatic relation. An energy equation analogous to (1.4.1) is
obtained except that the contribution of
to the kinetic energy is absent and, on the rhs, only the
horizontal components of the frictional forces do work. In QH only the terms underlined twice in Eqs. (1.3.2 ® 1.3.4) are neglected and, simultaneously, the
shallow atmosphere approximation is relaxed.
Thus all the metric terms must be retained and the full variation
of the radial position of a particle monitored. QH has good energetic credentials - they are the same as for HPE.
Importantly, however, it has the same angular momentum principle as NH - see Marshall et.al., 1997a. Again
in QH, the 3-d elliptic problem is
reduced to a 2-d one. Strict balance
between gravity and vertical pressure gradients is not imposed, however, since
the 2Wucosf Coriolis term plays a role in balancing ‘g’ in (1.3.4).
If the
ocean had a flat bottom, then (1.2.3) could readily be solved by projecting p on to the eigenfunctions of the 
operator with
boundary condition (1.2.4) applied at the (flat) upper and lower boundaries (see
Marshall et al., 1997a). Here such a modal approach can not be employed because
the geometry can be as complex as that of an ocean basin. However, the modal
approach points to the advantage of separating out, as far as is possible, the
depth-averaged pressure field.
For an
arbitrary function y we can define its vertical average 
as:
(1.5.1)
where H=H(l,f) is the
local depth.
The
vertically-averaged gradient operator is then:
(1.5.2)
The
vertical integral of the horizontal component of (1.1.1) is then, using the
rule (1.5.2) :
(1.5.3)
Since
there can be no net convergence of mass over the water column, we have,
assuming a rigid lid at the surface (see next chapter and appendix for the free
surface case):
(1.5.4)
At this
point, and following Bryan (1969), ocean modelers often introduce a stream
function for the depth-averaged flow. However, as argued by Dukowicz et al.
(1993) (see also Marshall et al., 1997a), it is advantageous to couch the
inversion problem in terms of pressure rather than a stream function. Applying
the horizontal gradient operator to (1.5.3) and making use of (1.5.4), we
obtain our desired equation for the depth-averaged pressure:
(1.5.5)
where
(1.5.6)
It is
now clear that if the ocean has a flat bottom (H=constant) then 
and equation (1.5.5)
does not have any pressure-dependent terms on the right-hand side and can be
solved unambiguously for 
. But if the depth of the basin varies from point to point
one cannot solve for the depth-averaged pressure without knowledge of 
. We can, however, make progress by separating the pressure
field in to hydrostatic, non-hydrostatic and surface pressure parts.
Add something about free-surface case here.
Let us
write the pressure p as a sum of three terms:
(1.5.7)
The
first term,
, only varies in the horizontal and is independent of depth.
The second term is the hydrostatic pressure defined in terms of the weight of
water in a vertical column above the depth z,
(1.5.8a)
where
(1.5.8b)
Note
that (1.5.8) is a generalized statement of hydrostatic balance, balancing
vertical pressure gradients with gravity and, in QH, also Coriolis and metric terms.
It
should be noted that the pressure at the rigid lid,
, has a contribution from both
and
(by definition,
, at z=0). In the
hydrostatic limit
= 0 and
is the surface pressure.
By
substituting Eq. (1.5.7) in (1.5.5) an equation for
results:
(1.5.9)
where 
is given by:
(1.5.10)
In HPE and QH the (doubly) underlined terms on the rhs of (1.5.9), which
depend on the non-hydrostatic pressure, are set to zero and
is found by solving:
(1.5.11)
The
vertical velocity is obtained through integration of the continuity equation
vertically:
(1.5.12)
In NH, instead, we first solve (1.5.11)
to give a provisional solution for
and then find
from the Elliptic
equation obtained by substituting (1.5.7) in to (1.1.1), and noting, as before,
that 
at each point in the fluid:
(1.5.13a)
where
(1.5.13b)
Eq.(1.5.13)
is solved with the boundary conditions:
(1.5.14)
If the
flow is ‘close’ to hydrostatic balance then the 3-d inversion converges rapidly
because 
is then only a small
correction to the hydrostatic pressure field.
The
solution 
to (1.5.13) and (1.5.14) does not vanish at the upper
surface, and so refines the pressure at z=0; it is in this sense that the 
obtained from
(1.5.11) is provisional. In the
interior of the fluid non-hydrostatic pressure gradients,
, drive motion:- w
is found by prognostic integration of the vertical velocity equation:
(1.5.15a)
where
(1.5.15b)
Note
that in (1.5.15), the vertical gradient of the hydrostatic pressure has been
canceled out with 
, Eq.(1.5.8b), rendering it suitable for prognostic integration.
The
method of solution employed in the HPE,
QH and NH models are summarized in Fig.3 below.
Fig.3. Outline of
hydrostatic (HPE), quasi-hydrostatic
(QH) and non-hydrostatic (NH) algorithms.
There
is no penalty in implementing QH
over HPE except, of course, some
complication that goes with the inclusion of cosf Coriolis terms and the relaxation of the shallow
atmosphere approximation. But this
leads to negligible increase in computation.
In NH, in contrast, one
additional elliptic equation - a three-dimensional one - must be inverted. However we show that this ‘overhead’ of the NH model is essentially negligible in
the hydrostatic limit (as the non-hydrostatic parameter N®0 - see appendix II) resulting in a non-hydrostatic
algorithm that, in the hydrostatic limit, is as computationally economic as the
HPEs.
In the momentum equations, the
forcing comes from the action of the wind at the sea surface and is represented
by a term 
such that:
(1.7.1)
where X is the stress such that:
(1.7.2)
where 
and 
are the zonal and
meridional components of the wind stress at the sea surface. Similarly,
buoyancy forcing at the sea surface induce forcing terms in the temperature and
salinity equations:
(1.7.3)
where 
is the specific heat
of the ocean and Q the heat flux such
that:
(1.7.4)
where 
is the air-sea heat
flux. Similarly:
(1.7.5)
where e is the freshwater flux such that:
(1.7.6)
where 
is a (constant)
reference salinity and E and P are the evaporation and precipitation
rates, respectively. Often, the heat and freshwater fluxes at the sea surface
are parameterized by restoring the surface temperature and salinity to observed
values. In this case, terms of the form 
, 
, where 
and 
are relaxation
time-scale constants, are added to, or replace, the surface terms in equations
(1.7.3) and (1.7.5).
The
forms of friction available in the model for momentum are Laplacian and
biharmonic frictions:
(1.7.7)
where 
and 
are (constant)
horizontal and vertical viscosity coefficients and 
is the horizontal
coefficient for biharmonic friction. These coefficients are the same for all
velocity components.
The
dissipation terms for the temperature and salinity equations have a similar
form except that the diffusion tensor can be non diagonal and have varying
coefficients (see appendix on the parameterization of eddy fluxes of heat and
salt).
(1.7.8)
where 
is the diffusion
tensor and 
the horizontal
coefficient for biharmonic diffusion. In the simplest case where the
subgrid-scale fluxes of heat and salt are parameterized with constant
horizontal and vertical diffusion coefficients, 
is reduced to a
diagonal matrix with constant coefficients:
(1.7.9)
where 
and 
are the horizontal
and vertical diffusion coefficients. Different values for temperature and
salinity can be specified. More physically-based parameterizations, such as
Gent-McWilliams - GM - or Green-Stone - GS -, are variants on (1.7.8), (1.7.9)
and are described in the appendix.
The
continuum equations employed are discretized using a finite volume strategy. In
this approach the domain to be simulated is divided into cells and the
evolution of the fluid is simulated, by stepping forward discrete forms of the
Navier-Stokes equations integrated over these cells. In the limit that the
cells have a regular cuboid structure the model is identical to a finite
difference model.
The variables are stepped
forward in time using a quasi-second-order Adams-Bashforth time-stepping
scheme. The pressure field, which
ensures that evolving currents remain non-divergent, is found by inversion of
Elliptic operators. In HPE and QH a 2-d Elliptic problem must be inverted; in NH the elliptic problem is 3-dimensional. A major objective is to make this 3-D inversion - the ‘overhead’
of NH - efficient and hence
non-hydrostatic modeling affordable.
The pressure field is separated in to surface pressure,
, hydrostatic pressure, 
and non-hydrostatic
pressure, 
and the component
parts found sequentially. In the 3-D
inversion for 
a preconditioner is
used which, in the hydrostatic limit, is an exact integral of the Elliptic
operator, and so leads to an algorithm that seamlessly moves from
non-hydrostatic to hydrostatic limits.
Thus, in the hydrostatic limit, the NH
is ‘fast’, competitive with the fastest ocean climate models in use today based
on the HPEs. But as the resolution is increased the model
dynamics asymptotes smoothly to the Navier Stokes equations.
We
write (1.1.1®1.1.4) in semi-discrete
form to second order in the time step Dt in which, as yet, only time is discretized, and using
the notation of chapter 1 (see Fig.3):
(2.1.1)
(2.1.2)
(2.1.3)
(2.1.4)
In
(2.1.1), we employ a ‘tracer’ parameter ‘q’
which takes on the value zero in the HPEs
and QH, and the value unity in NH.
Given
that the 
’s are known at time level n, an accurate explicit form for the time-stepping which involves
just two time levels is the Adams-Bashforth method (AB2) which makes use of
time-levels n and n-1 thus:
(2.1.5)
AB2 is
a linear extrapolation in time to a point that is just, by an amount e, on
the n+1 side of the mid-point
. AB2 has the
advantage of being quasi-second-order in time. A computational mode exists but
is damped by the e factor. Furthermore it can be implemented by
evaluating the G’s only once and
storing them for use on the next time-step. The code is written in a
sufficiently flexible manner that other time-stepping schemes could be readily
implemented.
The
time-step is limited by the following condition (see for example Potter, 1976;
p.69):
if 
where v is the fastest propagation velocity anywhere on the mesh of size D. Typically we set e to a value of 0.1.
It should be noted that if e=0 then AB2 is unstable in the inviscid case.
The
prognostic and diagnostic steps of our calculation in HPE, QH and NH are outlined in the schematic
diagram below - Fig.4. We now go on to
describe the spatial discretization employed and, in section 2.3, our Elliptic
inversion procedure.
Fig.4. Outline of the HPE, QH and NH algorithms.
The ocean is carved up in to a
large number of ‘volumes’ which can be called ‘zones’ or ‘cells’ - see
Fig.5. They take on a regular shape
over the interior of the ocean but can be modified in shape when they abut a
solid boundary - see Fig.5. We
associate tracer quantities with these cells; the cells have a volume 
and six faces 
, where ‘J’
denotes a generic tracer (such as T
and S). Except where they abut a solid boundary, the faces must be chosen
to coincide with our (orthogonal) coordinate system. If ‘x’, ‘y’ and ‘z’ are three orthogonal coordinates, increasing (nominally)
eastward, northward and upward - see Fig.5 - the faces normal to the ‘x’ axis have area 
, faces normal to the ‘y’
axis area 
, and faces normal to the ‘z’ axis area 
. Velocity components
are always normal to the faces. The geometry of the computational domain, and
the coordinate system employed within it (whether, Cartesian, spherical-polar,
cylindrical etc), is set up by prescribing the volumes and surface areas of the
faces of all the cells of which it is comprised.
Fig.5 The faces of
the cells are coincident with three orthogonal coordinate axes, sketched here
for a Cartesian grid (except where the cells abut a solid boundary). Velocities
are ‘face’ quantities defined normal to the faces of the cells; T, S
and p are ‘volume’ quantities.
In the
terminology of Arakawa and Lamb (1977), the grid used here is the 'C' grid.
The volume or cell quantities, p, r, T and S are defined as volume
averages over the cells. Prognostic
equations for these quantities are found using a `finite volume’ approach i.e.
by integration of the continuous equations over cells making use of Gauss’
theorem. For example, applying Gauss’ theorem to the continuity equation
(1.1.2) over a cell it becomes:
(2.2.1)
where
and similarly for the dy and dz terms.
The velocities are defined perpendicular
to the faces of the cells so ensuring that divergence and gradient operators
can always be naturally represented by symmetric operators. The divergence of
the flux of a quantity over the cell is:
(2.2.2)
where 
is the flux with
components defined normal to the faces and 
is the volume of the
cell. The presence of a solid boundary
is indicated by setting the appropriate flux normal to the boundary to zero.
Let us now consider the
discrete evaluation of the G’s for
the volume quantities required for integration of (2.1.4). The flux of, say, the temperature T through a ‘u-face’ of a cell is:
where 
, the temperature on the face, is given by :
the average of the temperature
in the two cells to which the face is common. The divergence of the flux of T over a cell,
, required in the evaluation of 
, for example, is:
(2.2.3)
where the 
is the spatial
operator in the x direction (and analogously
for y and z). It can be shown that
this simple average ensures that second moments of T are conserved i.e. the identity 
is mimicked by the
discrete form (2.2.3) if 
is chosen to be the
arithmetic mean. In (2.2.2) and (2.2.3), 
is the volume of the
cell in question.
It is clear from (2.2.3) that
velocities only appear as averages over faces of the cells and so we choose
only to define them there. To deduce
appropriate discrete forms of the momentum equations, however, we associate
volumes with these face quantities - 
- and again use Gauss’ theorem. Let us first consider the special case in which these volumes are
not shaved. Then the advecting terms
that make up the 
’s in Eqs. (1.3.2®1.3.4)
have the form (using continuity equation (1.1.2):
(2.2.4)
where ¾ is the spatial-averaging operator and
(2.2.5)
the average of the volumes of the cells on either side of
the face in question.
Exactly analogous expressions
are written down for 
and 
in
which the 
average of the ‘area times
velocity’ terms is replaced by, respectively, 
and 
.The forms (2.2.4, 2.2.5) guarantees:
-conservation
of the first moment because it is written in flux form,
-conservation
of the second moment because the fluxing velocities satisfy a non-divergence
condition:
exactly analogous to (2.2.1).
Velocity components are
staggered in space and so the evaluation of Coriolis terms in (1.3.2®1.3.4) involves spatial
averaging:
(2.2.6)
where 
. The above form
ensures energy conservation; see, for example, Arakawa and Lamb (1977).
Spatial
averaging of Coriolis has undesirable consequences - at coarse resolution the
w-field can be prone to grid-scale noise.
Considerable attempts have been made to ameliorate this problem - see
section 2.2.6 where the 'C-D' scheme, devised by Alistair Adcroft, is described.
The three components of the pressure gradient force are:
(2.2.7)
The form (2.2.7) ensures that
the discrete operator representing 
, which must be inverted to find the pressure field, is a symmetric
operator facilitating the use of conjugate-gradient methods - see section
2.3. It also leads to an appropriate
discrete analogue of energy conservation - see section 2.2.7 and, in
particular, Adcroft et al.(1997).
‘Laplacian’ diffusion terms are represented for cell
quantities thus:
(2.2.8a)
and for ‘face’ quantities:
(2.2.8b)
where
is
given by (2.2.5).
Higher order dissipative forms
are obtained by repeated application of the 
operators defined above.
Include form of discrete forcing terms here.
The spatial averaging required
in the evaluation of Coriolis terms on the ‘C’ grid is a significant
disadvantage when the model is used at resolutions which are coarse relative to
the Rossby radius of deformation. A
checkerboard mode in the horizontal divergence field can readily be excited and
implicit treatment of Coriolis terms is significantly more complicated than on
unstaggered grids (such as the ‘B’ grid, for example). As shown in Adcroft (1995), however, both of
these difficulties can be significantly ameliorated by use of a C-D grid in
which a ‘D’ grid is used to step forward velocity components used in the
evaluation of Coriolis terms.
Include the discrete equations here.
The scheme has been
particularly useful in coarse-resolution studies of ocean circulation in which
the horizontal resolution is many times the Rossby radius of deformation.
The above discrete forms give
the model excellent integral conservation properties. They conserve integrals of mass, heat, salt, variance of
temperature and salinity and total energy, when effects of mixing and transfer
across boundaries are absent. The
discrete analogue of kinetic energy is:
These integral properties are
not compromised by the use of shaved cells provided care is taken in the definition
of control volumes for variables adjacent to the boundary (see appendix). The model does not conserve enstrophy,
however.
Include
section on shaved cells (maybe in appendix?)
We now formulate the discrete
analogues of the 2-d and 3-d elliptic problems for the pressure which
guarantees that the velocity fields evolving from n to n+1 according to
(2.1.1) and (2.1.2) satisfy the non-divergence condition (2.1.3) at step n+1.
The parallel implementation of the preconditioned conjugate-gradient
methods used to invert them are discussed in section 2.3.2.
Two-dimensional
Because of the ‘finite volume’
approach adopted in our treatment of volume quantities, in which velocities are
defined normal to the faces of the cells, the divergence operator has a natural
form leading to a local Elliptic operator which has a five-point stencil. By setting q=0 in the momentum equations (2.1.1) and summing them over the
whole depth of the ocean, invoking the continuity equation (2.2.1) and applying
boundary conditions (1.2.1), the following equation for 
results:
(2.3.1)
where
(2.3.2)
Here 
is the discrete
analogue of 
, a vertical integral over the whole depth of the ocean; we
sum over the vertical faces of the cells each of depth Dz|k making up the
column of ocean. The ‘divh’
and ‘gradh’ operators are horizontal components of (2.2.2) and
(2.2.7), respectively.
Note that on the rhs of (2.3.1) we have retained a
term,
which ensures that 
as the
model steps forward. This ‘relaxation’ technique is often used in solving the
Navier-Stokes equations- (see, for example, Williams;1969). In the present context it obviates the need
to step forward barotropic equations for
separately. In our
method only the prognostic equations for interior velocities, Eq (4.1.1), are
stepped forward and the horizontal divergence of the depth-integrated
horizontal velocities at each time-step evaluated and used to modify the
source-function to the 2-d Elliptic problem accordingly, so ‘tying together’
the velocity and the pressure field.
The Elliptic problem
Eqs.(2.3.1) and (2.3.2) can be written in the concise matrix notation:
(2.3.3)
where 
is a
symmetric, positive-definite matrix
comprised of
and 
(matrix
representations of the ‘div’ and ‘grad’ operators), 
is a column vector of
surface pressure elements and 
a column vector
containing the elements of the rhs of (2.3.1).
The system can thus be solved using a standard conjugate-gradient
method, preconditioned for efficient solution on a parallel computer - see
section 2.3.2.
In summary, then, because (i)
the divergence constraint on the evolving velocities is applied in integral
form employing Gauss’ theorem and (ii) the kinematic boundary condition v.n=0 (where n is a unit vector normal to a solid boundary) is applied at the face
of a cell, then the form of (2.3.3) is unchanged, even if the face of a cell which abuts the solid boundary is inclined
to coordinate surfaces. This
property enables one to rather easily ‘shave’ zones which abut boundaries to
improve the representation of flow over topography (see appendix).
Three-dimensional
In non-hydrostatic calculations
a three-dimensional elliptic equation must also be inverted for 
to ensure that the local
divergence vanishes. The appropriate
discrete form can be deduced in a manner that exactly parallels that which was
used to deduce (2.3.1). Taking Ñh.(2.1.1)+
, invoking (2.1.3) and boundary conditions (1.2.1), the
resulting elliptic equation can be written:
(2.3.4)
where 
, like 
, is a
symmetric,
positive-definite matrix representing the discrete representation of 
, but now in three dimensions. Here 
and 
are (1´N) column vectors containing the source function (the
discrete form of Eq.(1.5.13)) and non-hydrostatic pressure, in each of the 
cells in to which
the fluid has been carved. The column
vectors have singly subscripted elements
in which the elements are first enumerated in each vertical
column, and only then in the horizontal directions thus:
where i is an index increasing eastwards, j is an index increasing southwards and k increases downwards.
Although huge, of size 
, 
has a particularly
simple form which can be exploited to devise a highly efficient method of
inversion.
has 7 diagonals representing the coupling in the three space
dimensions. In any particular row of
the matrix, the three leading diagonals are the coefficients multiplying the
pressure in the same vertical column of
the fluid; the four diagonals in the wings are the coefficients multiplying
pressures in cells in the same horizontal
plane. The inner three diagonals
can usefully be blocked and arranged as shown below:
(2.3.5)
Here each block 
is a tri-diagonal
matrix representing 
and e is a
diagonal matrix representing 
. D and e are matrices of size 
if there are
cells
in each column of the ocean:- there are 
such blocks, one for
each column of the ocean.
One further property of 
should be noted. If the vertical dimension of a cell is much
smaller than its horizontal extent (as is often the case because the ocean is
much shallower than it is wide) then the 
»
and 
is dominated by the
blocks D along its diagonal. The elements e are smaller than those of D
by an amount 
. Moreover, the
tracer parameter ‘q’ always appears
as a multiplier of the e matrix -
see (2.3.5); thus when q is set to
zero (corresponding to the hydrostatic limit), 
is comprised only of the blocks D and so can readily be inverted.
Thus if we ‘precondition’ (2.3.4) by premultiplying it by a matrix which
is comprised of the inverse of these blocks, that preconditioner will be an
exact inverse of 
in the hydrostatic
limit. These properties of 
will be exploited in
our chosen method of solution.
Many standard references to
preconditioned conjugate gradient methods exist - see, for example, Press et
al; 1988 and references therein. We
briefly outline the ‘nub’ of the method here, emphasising the use we make of
preconditioners.
Our problem is to find p given A and f (see (2.3.3) and (2.3.4)), where:
(2.3.6)
and A is a
symmetric, positive-definite matrix.
In serious modelling
applications the size of A is too
large for direct methods to be useful in 3-dimensions and often in 2-dimensions
too. Since A does not change in time, it could be inverted once and
stored. However, although A is sparse, its inverse is dense and
so operating with its inverse would involve N2
multiplications; an unrealistic task given that typically N >> 100000. So we
adopt an iterative procedure, preconditioned conjugate gradient iteration,
which exploits the sparseness of A
and its diagonal dominance. The
procedure involves repeated multiplication of the iterative solution by A and by another sparse matrix K, an approximate inverse of A; K
is called the preconditioner.
The algorithm can be understood thus:
Let us pre-multiply Eq.(2.3.6)
by a (carefully chosen) matrix K,
which is an approximate inverse of A,
so that 
. Then (2.3.6) can be
written:
where 
. To the extent that
|C| ®0, the above suggests the
following iterative scheme - 'i' is the
iteration step:
where
is called the search
direction and
is the residual
vector. The r's and b's can be
deduced from the previous iteration using the relations:
The above iterative procedure
can be accelerated by choosing those search directions that minimize the
magnitude of the residual vector as measured by the inner product 
. Our chosen method,
the conjugate gradient method, takes a form much like the simple scheme
outlined above, but the search directions are optimized. It can be written:
(2.3.7)
Here [ , ]
is the inner product of two vectors.
The conjugate-gradient
procedure involves, at each iteration, multiplication of vectors by A
and K and the evaluation of two
global sums to form the inner product.
The cost of the preconditioner is one operation by it per iteration.
We have designed
preconditioners, K, for both the 2-d
and 3-d problems, so that (i) it can be easily stored (ii) the number of
operations required when multiplying by it, is as small as possible and (iii)
it is a good approximation to A-1
so the iterative procedure converges rapidly.
Because the true inverse is dense, our choice will be a compromise
between these, sometimes conflicting, ideal properties.
Two-dimensional
At each horizontal location (2.3.3) can be written:
where the superscripts denote
the center, west, east, north and south locations in the operator’s five point
stencil. The leading diagonal of A2d will be comprised of the
and the four non-zero
off-diagonals will be comprised of the 
. Because of the
dominance of the leading diagonal, to a first approximation:
One of
the simplest forms that could be chosen for 
, then, is to suppose that it is a diagonal matrix comprised
of elements equal to the reciprocal of the corresponding elements along the
leading diagonal of A. However, this can be improved by
constructing an approximate local inverse as follows.
At the
five points in the stencil surrounding ‘c’ then, to the same level of
approximation, we may write:
where 
denotes the point to
the west etc.
Hence we can make a better
approximation to 
thus:
To arrive at a symmetric preconditioner we replace 
etc
thus:
giving, finally, our local preconditioner for use in
(2.3.7):
(2.3.8)
Although simple, (2.3.8) is
highly effective, reducing the number of iterations required for convergence
by, typically, a factor of four.
Three-dimensional
After considerable
experimentation we have chosen a block-diagonal preconditioner, a matrix whose
diagonal is made up of the inverse of the tri-diagonal matrices D defined above:
(2.3.9)
Evaluation of the inner products to compute a and b in
(2.3.7), involves forming the vector x,
where
or, since 
, where D is the matrix of diagonal elements of
(2.3.5).
If our ocean model has many
levels, then D-1 will be
dense and so storing and multiplying by D-1
is also demanding of resources.
Instead, we exploit the fact that D
is tri-diagonal, and use ‘LU’ decomposition to solve the preconditioning
equations for x in the form:
We write D=LU where L
is a lower-triangular matrix and U
an upper triangular matrix which have the form:
, 
and 
It is easy to see, that if 
we have a system of
simple equations:
and so on.
Finding the inverse of D, or solving problem Dx=r
is equivalent to solving the two sets of equations:
(2.3.10)
first for z and then for x. This is straightforward because of the
triangular structure of L and U.
Importantly the number of operations required to solve (2.3.10) for x scales as NZ compared to NZ2 if we had used the inverse of D directly.
The resulting preconditioner
was found to be a good compromise; its use significantly reduces the number of
iterations required to find a solution to (2.3.6) because it is an acceptable
approximation to the inverse of A,
it is sparse and, most importantly, its application requires no communication
across the network in the data-parallel implementation of the algorithm.
Not complete yet
The current model implementation uses a
conservative form of Fortran 90. The array syntax notation feature of Fortran
90 is heavily used, however, the more elaborate features of Fortran 90 such as
overloading of operators, modules, array valued functions are avoided (the
performance of these constructs varies widely from one compiler to another).
Apart from the array notation, the only new syntax elements that experienced
Fortran 77 or C programmers need to understand to work with the model, are the
Fortran 90 SHIFT primitives.
In the model implementation the
grid cells are mapped to regular three-dimensional arrays. As shown in Fig.6
the index i is used to index the x-axis, j to index the y-axis and
k to index the z-axis. The rule that coordinate (i=1,j=1,k=1) of any variable must fall within
the western most, southern most and upper most cell and be as far toward the
west, south and top of that cell is used throughout the code. This is illustrated
in Fig.6.
Fig 3 .The discrete domain is mapped to multi-dimensional arrays using an
enumeration scheme that follows the two rules illustrated. Both the major
prognostic variables and the intermediate terms abide by this convention.
See below a flow chart of the code. Some procedures are
not mentioned (diagnostic procedures for example).
MAIN
|----CONTROL
|
|----SET_DEFAULTS Set default values for model parameters
| |----DATE Record date and time of run
| |----MACHINE Record machine identifying string
|
|----READ_DATA Read file specifying model parameters
| |----READ_TOPOGRAPHY Read
land/water map
|
|----INITIALISE
| |----INITIALISE_MASKS Set land/water masks arrays
| |----INITIALISE_GRIDDING Set up arrays
expressing geometry
| | |----INITIALISE_CARTESIAN_GRID
| | |----INITIALISE_SPHERICAL_POLAR_GRID
| |----INITIALISE_OPERATORS Set up differencing operators
| |----INITIALISE_FIELDS Set initial values of model variables
| | |----INITIALISE_U
| | |----INITIALISE_V
| | |----INITIALISE_S
| | |----INITIALISE_T
| |----INITIALISE_CORIOLIS Set Coriolis term operator
| |----INITIALISE_METRIC Set
metric term operator
| |----CG3DAI Define
3d laplacian "A matrix" operator
| |----CG3DPI Define
3d "A matrix" preconditioner
| |----CG2DAI Define
2d laplacian "A matrix" operator
| |----CG2DPI Define
2d "A matrix" preconditioner
|
|----CHECK_INPUT Validate input data
|
|----PRINT_HEADING Print out run data summary page
|
|----MODEL Perform
integration
|----RHO_WGRID_CALC Evaluate initial density
|
MAIN
LOOP =====>
|
|----FORWARD Control explicit stepping of u, v, and w
| |----TERMS0 Set up temporary arrays
| |----TERMS1 needed throughout stepping
| |----GU_CALCULATE Calculate explicit dU/dt
| | |----EXTERNAL_FORCING_U Add forcing
| |----GV_CALCULATE Calculate
explicit dV/dt
| | |----EXTERNAL_FORCING_V Add forcing
| |----GW_CALCULATE Calculate explicit dW/dt
|
|----PFIND
Diagnose pressure field
| |----PFIND_HSTATIC Find hydrostatic pressure
| |----PFIND_SURFACE Find surface pressure
| | |----CG2D 2d preconditioned conjgrad
| | | |----CG2DMA Multiply by 2d laplacian
| | | |----CG2DPC Precondition
| |----PFIND_3D Find remaining pressure
| | |----PFIND_RHS3D Form right hand side
| | |----CG3D 3d preconditioned conjgrad
| | | |----CG3DMA Multiply by 3d laplacian
| | | |----CG3DPC Precondition
|
|----UPDATE_VELOCITIES
|
|----EPARAM_SET_MIXING Computes 3d diffusion tensor
|
|----INC_TRACER Step forward temperarture
| |----EXTERNAL_FORCING_TEMPERATURE Add
forcing
|
|----INC_TRACER Step forward salinity
| |----EXTERNAL_FORCING_SALINITY Add forcing
|
|----CADJTRACE Mix unstable water column
|
|----RHO_WGRID_CALC Evaluate density field
|
TO MAIN LOOP <=====|----OUT Dump out model fields
Include
the list here with a brief description (indicate the 'parent' procedure).
Include an alphabetical list of variables here with the following
information: type (integer, real, character etc…), units (for real variables),
brief description, first occurrence in the code. Use table.
Notation : Fortran variable names
of the code are written in italic.
Subroutine and file names are written in CAPITAL LETTERS.
leave for now
This section is meant to review
the parameters that need to be set up before one can run the model. Unless
noted otherwise, these parameters are Fortran variables declared in file
PARAMS.h and set up in routine SET_DEFAULTS. The default values can be changed
by means of a namelist. The list that follows is not exhaustive. For example,
parameters related to the forcing files are expected to be user specific and
are not mentioned here. Others are left to the appendices. However, we have
attempted to include all the ones that will need to be set up no matter what
configuration you use. For the parameters not described here, look into file
PARAMS.h and routine SET_DEFAULTS. In what follows, the parameters are grouped
into categories related to the geometry and the equations solved in the model.
We start with some constants
describing the state of the ocean. The real variables g, Cp, RoNil represent respectively gravity (in
m s-2), the specific heat (in J K-1 kg-1), and
the reference density (in kg m-3) of the ocean.
The vertical coordinate used in
the model is actually pressure so that the same code can be used to simulate
fluid motions in the atmosphere. In the namelist, variables depending on the
choice of vertical coordinate (such as vertical viscosity and diffusivity
coefficients) have their vertical length expressed in meters. Later in the
code, this vertical length is converted in pressure unit (this involves
typically a simple multiplication or division by RoNil*g).
-
dimensions
The dimensions of the domain Nx, Ny,
Nz, are declared and set in file
SIZE.h. The integer variables L, M, N
in the namelist must be set to the values of Nx, Ny, Nz, respectively.
-
grid
The character variable gridStyle controls the type of grid. Two
versions are currently available. For a Cartesian grid, set gridStyle to 'Cartesian'. Then set up
the horizontal grid spacing dM (in
m). For a spherical polar grid, set gridStyle
to 'Spherical Polar'. Then specify the resolution along latitude and longitude dPhi and dTheta (in degrees) and the minimum latitude phiMin (in degrees). Whatever grid you choose, also specify the 1D
array representing the vertical layers thickness deltaz (in m) in the namelist. Note that in the code, the layer
thicknesses are expressed in pressure units (Pa) and are represented by the variable
delps (equal to deltaz multiplied by RoNil*g).
-
topography
The
model code applies without modification to enclosed, periodic (x or y) and
double periodic (x and y) domains. The geometry of the model is controlled by
the topography file. A real 3D array pmask
is defined having the value 1 on water and 0 on land. Periodicity is assumed by
default and is suppressed by defining zero open interface areas for cells at
the limits of the computational domain.
Note: The present elliptic
solver procedures CG2D and CG3D assume that they are solving a single global
problem. This means that all cells with non-zero volumes must interconnect,
either directly or indirectly, through the surface. If disconnected basins are
defined, the results are unpredictable.
Initial reference temperature
and salinity profiles need to be provided to compute the perturbation density
and are represented respectively by the 1D arrays ths (in oC) and ssppt
(in ppt). Two forms of equation of state are available in the baseline
implementation and are controlled by the character variable equationofState. If set to 'linear',
then specify the constant expansion coefficients tAlpha (in K-1) and sBeta (in ppt-1) (through the variables alphaTemperature and alphaSalt in file PARAMS.h). If you set equationofState to 'polynomial', density
is evaluated from a third order polynomial expansion of the full state
equation. The polynomial coefficients are read from a file called
POLY3D.COEFFS. Any other arbitrary formulae for deriving density can readily be
plugged into the model.
The logical variable nonHydrostatic controls whether or not
you want the model to run in a nonhydrostatic regime (NH) as described in the first two chapters of this manual. If nonHydrostatic is set to .TRUE., the
model will step forward for the vertical velocity and solve a 3D elliptic
equation for the pressure. Note that the vertical velocity can still be
diagnosed from the continuity equation if the logical variable stepW is set to .FALSE. In the
nonhydrostatic regime, you also need to set up the logical variable buoyancy that controls whether or not
you want to include the buoyancy term in the vertical momentum equation. The
distinction between hydrostatic (HPE)
and quasi-hydrostatic (QH) regime is
controlled by the logical variable verticalCoriolis
(see next section, momentum equations).
Unless noted otherwise, what
follows is not applicable to the vertical velocity w if the variable nonHydrostatic
has been set to .FALSE. (HPE or QH regime).
The logical variable momentumStepping controls whether or not
you want to step forward for the velocity components u, v and w. If set to .TRUE., the following
parameters have to be specified.
-
advection terms
If you want to include these
terms in the momentum equations, set the logical variable momentumAdvection to .TRUE.. A centered second-order advection
scheme is then used.
-
metric terms
If you want to include these
terms in the momentum equations, set the logical variable metricTerms to .TRUE.. Note that all the metric terms are
included. This choice is to be made only if the grid is spherical polar (the
variable is automatically set to .FALSE. if the grid is Cartesian).
-
Coriolis terms
If you want to include these
terms in the momentum equations, set the logical variable Coriolis to .TRUE.. If this is the case, also set the real variable
raja and the logical variable verticalCoriolis. The former is a number
between 0 and 1 that controls the "C-D scheme" (see chapter 2). If raja is set to a value greater than 0.,
horizontal velocity components stepped forward on a 'D' grid are used in the
evaluation of the Coriolis terms with a weight given by raja. If raja is set to
0., the Coriolis terms are estimated using the 'C' grid velocity components
only. The variable verticalCoriolis
controls whether or not you want to include the cos (latitude) terms in the
zonal and vertical momentum equations (I
think the name verticalCoriolis is misleading. The additional terms appear both
in the zonal and vertical momentum equations. I would suggest to replace the
name verticalCoriolis by cosphiCoriolis). If verticalCoriolis is set to .TRUE. and nonHydrostatic to .FALSE., the model is in the quasi-hydrostatic (QH) regime. If the grid is Cartesian,
you have to specify the Coriolis parameter dependence on latitude with the
character variable CoriolisPlane. If
this variable is set to 'betaplane', set the reference Coriolis parameter fcori (in s-1) and its
derivative with respect to latitude beta
(in m-1 s-1) (Note
: the Beta plane formulation is not available yet). If CoriolisPlane is set to 'fplane', only fcori needs to be set. Note that you can set up a non-rotating
problem on a Cartesian grid by simply setting fcori to 0.
-
dissipation terms
If you want to include these
terms in the momentum equations, set the logical variable momentumDiffusion to .TRUE.. If this is the case, specify the
horizontal and vertical viscosity coefficients CAPUV and RKPUV respectively
in the namelist (in m2 s-1). Note that in the code, these
variables are called a2MomXY and a2MomZ and that the unit of the latter
is Pa2.s-1 (multiplication by [RoNil*g]2). If you also want to include biharmonic
diffusion terms in the momentum equations, set the logical variable momBhDiffusion to .TRUE.. Then, set the
value of the horizontal biharmonic diffusion coefficient C4PUV (in m4 s-1) in the namelist. Note that
in the code, this variable is called a4MomXY.
-
forcing terms
This section only applies to
the horizontal velocity components.
If you want to include these
terms in the momentum equations, set the logical variable momentumForcing to .TRUE.. Then, additional terms will be added to
the internally induced tendency terms calculated by the model. It is assumed
that you have wind data available (or that you have prescribed an analytical
forcing) and a routine that estimates the forcing terms at each time step.
-
boundary conditions
Slip or no-slip conditions on
side walls and at the bottom and top are controlled by logical variables freeslipSide, freeSlipBottom and freeSlipTop,
respectively. If these variables are set to .TRUE., free slip boundary
conditions are applied (otherwise no-slip conditions are applied).
The logical variable pressureStepping controls whether or not
you want to solve for the pressure field. If you do, the following parameters
need to be specified :
-
hydrostatic pressure
The logical variable seperatePh controls whether or not you
want to calculate the hydrostatic pressure anomaly. If you do, then specify the
surface boundary condition with the character variable phSurfBC. Two options are available : the hydrostatic pressure
anomaly is set to zero at the surface (phSurfBC='zero')
or is set to the density of the first layer times half the layer depth (phSurfBC='topBoxPhIntegralofRho').
-
surface pressure
The logical variable seperatePs controls whether or not you
want to solve for the surface pressure. If you do, the real constant variables freeSurfFac, psGam, uDivAlpha, and gBaro need to be set up and control
whether or not you want to have a rigid lid or an implicit free surface. For a
rigid lid, set freeSurfFac, psGam, uDivAlpha, and gBaro to
0., 1., 0., and g (gravity)
respectively. To control the operation of the 2D iterative solver routine CG2D,
specify the values of the variables toler2d,
max2dIt, freqCheckToler2d, divhMax,
divhMin. Finally, the logical
variable stericEffect controls
whether or not you want to include a steric anomaly in surface elevation. If
you do, then specify the value of variable rhorefSurf,
reference density for surface layers.
-
non hydrostatic pressure
This part is relevant only if
the variable nonHydrostatic has been
set to .TRUE. (see section 5.2.3). Then specify the values of the variables toler3d, max3dIt, freqCheckToler3d,
divMax, divMin that control the operation of the 3D iterative solver
routine CG3D. Furthermore, the logical variable findD2pnhDZ2 controls whether or not you want to solve the
Dirichlet problem for the vertical gradient of nonhydrostatic pressure.
This section also applies to
other tracers (if necessary). The logical variable tempStepping (saltStepping)
controls whether or not you want to step forward for temperature (salinity). If
you do, it is assumed that you have data defining initial conditions (or you
can specify initial vertical profiles through the variables ths and ssppt, see section 5.2.3). Moreover, the following parameters need
to be specified.
-
advection terms
If you want to include these
terms in the temperature (salinity) equation, set the logical variable tempAdvection (saltAdvection) to .TRUE.. If this is the case, choose the tracer
advection scheme with the character variable TradvecScheme. Two options are available : 'upstream' or
'centered'. (Note : in the current
version, TradvecScheme is not declared in PARAMS.h but is a local variable of
routine INC_TRACER). In the current version, the advection scheme is common
to all tracers.
-
dissipation terms
If you want to include
diffusion terms in the temperature (salinity) equation, set the logical
variable tempDiffusion (saltDiffusion) to .TRUE.. Then, choose
the tracer diffusion scheme with the character variable TrdifScheme. Several options are available. Here, we only consider
the case of constant horizontal and vertical diffusivity coefficients. For
other cases, see appendix to know what parameters need to be set up. So, if TrdifScheme is set to 'constant', the
variables CAPTH (CAPS) and RKPTH (RKPS) in the namelist need to be set and represent the
horizontal and vertical diffusivity coefficients, respectively (in m2 s-1). Note that in the code, the
variables are named a2TempXY (a2SaltXY) and a2TempZ (a2SaltZ) and
that the latter is expressed in Pa2 s-1. If you also want
to include biharmonic diffusion terms in the temperature (salinity) equation,
set the logical variable tempBhdiffusion
(saltBhdiffusion) to .TRUE.. Then,
set the value of the horizontal biharmonic diffusion coefficient C4PTH (C4PS) in the namelist (in m4 s-1). Note that in the code, the variable is named a4TempXY (a4SaltXY).
-
forcing terms
If you want to include these
terms in the momentum equations, set the logical variable tempForcing (saltForcing)
to .TRUE.. Then, additional terms will be added to the internally induced
tendency terms calculated by the model. It is assumed that you have air-sea
heat and freshwater fluxes data available (or that you have prescribed an
analytical forcing) and a routine that estimates the forcing terms at each time
step.
-
convective adjustment
The logical variable ConvectiveAdjustment controls whether or
not you want the static instabilities to be dealt with by a convective
adjustment scheme. If you do, specify the value of variable CadjFreq, frequency (in s) with which
the convective adjustment scheme is applied. For details about the scheme, see
appendix.
-
run duration
The starting time, ending time
and time step of an integration need to be set up with the integer variables IPUSTTTMP, IENDSECTMP and IDELT in
the namelist (in s). Note that in the code, these parameters are real variables
named startTime, endTime and delt. If the
starting time is non-zero, the model will try to read in a set of checkpoint
files from a previous integration (see section 5.3 for a description of the
checkpoint files).
-
forcing controls
Certain parameters related to
the reading of the forcing files might need to be set up (for example after a
restart). In general, this will depend on the particular user configuration.
Therefore, this is not extended any further.
-
output
Real constant variables
defining frequencies (in s) with which output files are written on disk need to
be set up. dStateFreq is the
frequency with which the state of the model is saved. dPUFreq is the frequency with which pick-up files (suitable to
restart an integration) are saved. dCkptFreq
is the frequency with which "rolling" checkpoint files are saved. See
section 5.3.3 for a description of the output files.
Type
make wherever you are. There is a good chance it will work! To set the
precision, look into utility r8 (in directory ProgUtils). This utility is
called by the makefile to "preprocess" the procedures of the code
converting all the 'REAL' statements into 'REAL*4' or 'REAL*8'. The makefile
will create an executable named MITGCM. Use the UNIX command 'size' to know how
much memory is required to execute the code on a single processor.
First, we review the input
files needed to run the model :
-
namelist (file named DATA),
-
topography (file named PMASK.SUN.B),
-
polynomial coefficients (file named POLY3D.COEFFS) if the
full equation of state is used,
and if necessary :
-
initial temperature and salinity,
-
wind data,
-
heat and freshwater fluxes data.
If you want to restart the
code, pick-up files are also needed (see next section). To run the code, type
MITGCM < DATA > OUTPUT. The building of scripts to make the submission
and execution of the code more automatic is left to the user at this point.
Depending upon the values of
the various frequencies mentioned in section 5.2.7, output files describing the
state of the model will be written on disk. The name of these files is written
in the following way : "XXX.0000nIter", where XXX is a capital letter
or a group of capital letters indicating the quantity being saved (for example
U for the zonal velocity) and nIter is a chain of figures indicating the instant
(in number of time steps) at which saving occurs. The maximum value for nIter
is 109 and 0s will fill in
the space after XXX when nIter is less than its maximum value as shown above.
Below follows a list of the output files. Unless noted otherwise, the different
quantities are expressed in SI units.
U : zonal
component of velocity vector
V : meridional
component of velocity vector
W : vertical
component of velocity vector (in Pa s-1)
T : potential
temperature (in oC)
S : salinity
(in ppt)
PS : surface
pressure (in m)
PH : hydrostatic
pressure (in m)
PNH : nonhydrostatic
pressure (in m)
PTOTAL : total
pressure (PS+PH+PNH) (in m)
GTNM1 : temperature
tendency at time step n-1 (in oC s-1)
GSNM1 : salinity
tendency at time step n-1 (in ppt s-1)
GUNM1 : zonal velocity
tendency (except pressure term) at
time step n-1
GVNM1 : meridional
velocity tendency (except pressure term) at
time step n-1
GWNM1 : vertical
velocity tendency (except pressure term) at
time step n-1 (in Pa s-2)
PNM1 : total
pressure at time step n-1
UAJA : 'D'
grid zonal velocity (used in C-D scheme)
VAJA : 'D'
grid meridional velocity (used in C-D scheme)
UBXBYNM1 : interpolated
'D' grid zonal velocity (from 'C' grid component) at time step n-1
VBYBXNM1 : interpolated 'D'
grid meridional velocity (from 'C' grid
component) at time step n-1
These are the files needed to
restart an integration (pick-up files) and their output is controlled by the
variable dPUFreq (see section 5.2.7).
Of course, if you are stepping forward other quantities (other tracer, ice…),
they will need to be output along with the files above. The variable dStateFreq controls the output of the
quantities U, V, W, T, S, PTOTAL only, from which most diagnostics can be made.
Rolling checkpoint files, the output of which are controlled by the variable dCkptFreq, are the same as the ones
above (pick-up files) but are named differently. Their name end with the suffix
'ckptA' or 'ckptB' instead of '0000nIter'. They can be used to restart the
model but are overwritten every other time they are output to save disk space
during long integrations.
First, we describe briefly each
case. Refer to the references for more detail.
This case corresponds to the
numerical experiments of Holland and Lin (1975). In this study, the authors
explore the generation of mesoscale eddies and their contribution to the
oceanic circulation. Their model has two homogeneous layers on a simple closed
rectangular basin 1000 km ´ 1000
km with a flat bottom at 5 km depth, and is forced by a steady wind-stress that
varies in a simple sinusoidal manner with latitude and has an amplitude of 0.1
Pa. The horizontal resolution is 20 km on a b-plane.
For sufficiently small viscosity values, the model, spun up from rest, exhibits
mesoscale eddies and eventually reaches a statistical steady state. Analysis of
the energy budget shows that the eddies are generated by baroclinic
instability.
The configuration given in
5.4.5 corresponds to the preliminary experiment described in Holland and Lin
(1975). In our height-coordinate model, quasi reduced-gravity conditions can be
configurated by prescribing an initial temperature profile mimicking the
reduced-gravity used in Holland and Lin set-up.
This case corresponds to the
numerical experiment of Jones and Marshall (1993) done with the present code.
In this study, the authors explore open ocean deep convection. The model has
dimensions 32 km ´ 32 km ´ 2 km and 19 levels on the
vertical. Periodic boundary conditions are applied on the x and y directions.
The horizontal and vertical resolutions are 250 m and 100 m, respectively. The
model is initialized with a homogeneous, resting ocean governed by a linear
equation of state without salt. Constant cooling corresponding to a heat loss
of 800 W m-2 is induced over a disk of radius 8 km at the center of
the domain and over a depth of 200 m from the surface. The model is run in the
nonhydrostatic mode and convective plumes develop. These plumes efficiently mix
the water column and generate a dense chimney of fluid, which subsequently
breaks up through the mechanism of baroclinic instability to form spinning
cones.
The configuration given
hereafter corresponds to the reference simulation of Jones and Marshall (1993).
This case corresponds to the
wind-driven channel experiment of Visbeck et. al. (1997) done with the present
code. In this study, the authors compare eddy-resolving 3D simulations with 2D
simulations in which the eddies are parameterized in order to test and assess
different schemes for specifying eddy transfer coefficients. The model has
dimensions 1500 km ´ 500 km
´ 4.5 km on a f-plane. The grid spacing is 10 km on
the horizontal and varies on the vertical from 25 m in the upper layers to 400
m at depth (20 vertical levels). Periodic conditions are applied at the
meridional boundaries. The model is initialized with an exponential temperature
profile of scale depth 900 m and is forced by a cosine wind stress of maximum
strength 0.2 Pa at the channel center. After six years or so, a steady state is
reached where the input of potential energy by the wind is balanced by release
of potential energy due to baroclinic eddies.
This case corresponds to
an-going study of the large-scale global ocean circulation. The spherical
domain extends from 80oS to 80oN in latitude and a
realistic topography is used. The horizontal resolution is 4 o and
the vertical grid spacing varies from 25 m in the upper layers to 500 m at
depth (20 vertical levels). The full equation of state is employed. The model
is initialized with the Levitus climatology and forced by climatological
forcing, namely the ….
The parameter values are
summarized on Table 5.1.
Table 5.1
Model configuration for the test cases. "-" indicate not
applicable or that the model results are not sensitive to the values chosen.
|
|
2-layer
|
convection
|
Channel
|
Global
4 deg
|
|
Nx
Ny
Nz-1
|
100
100
2
|
132
132
19
|
150
50
20
|
90
40
20
|
|
gridStyle
|
'cartesian'
|
'cartesian'
|
'cartesian'
|
'spherical
polar'
|
|
dm (km)
|
20
|
0.25
|
10
|
-
|
|
dphi
dTheta
phiMin (deg)
|
-
|
-
|
-
|
4
4
-80
|
|
deltaz (min)
(max)
(m)
|
1000
4000
|
100
|
25
400
|
25
500
|
|
Geometry
Topography
|
rectangular
box
flat
bottom
|
squareboxd-periodic
flat
bottom
|
channel
s-periodic
flat
bottom
|
sphere
realistic
topo
|
|
EquationofState
ths (oC)
ssppt (ppt)
tAlpha (K-1)
sBeta (ppt-1)
|
'linear'
20 -
10
0
?
0
|
'linear'
?
0
2 ´ 10-4
0
|
'linear'
exponential
?
2 ´ 10-4
?
|
'polynomial'
0
0
-
-
|
|
nonHydrostatic
|
.FALSE.
|
.TRUE.
|
.FALSE.
|
.FALSE.
|
|
momentumStepping
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
momentumAdvection
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
metricTerms
|
-
|
-
|
-
|
.TRUE.
|
|
Coriolis
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
VerticalCoriolis
|
.FALSE.
|
.FALSE.
|
.FALSE.
|
?
|
|
raja
|
0
|
0
|
0
|
?
|
|
CoriolisPlane
|
'betaplane'
|
'fPlane'
|
'fPlane'
|
-
|
|
fcori (s-1)
|
7.3 ´ 10-5
|
10-4
|
10-4
|
-
|
|
beta (m-1 s-1)
|
2 ´ 10-11
|
-
|
-
|
-
|
|
MomentumDiffusion
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
CAPUV
RKPUV (m2 s-1)
|
330
0
|
5
0.2
|
?
10-3
|
?
?
|
|
MomBhDiffusion
|
.FALSE.
|
.FALSE.
|
?
|
?
|
|
C4PUV (m4 s-1)
|
-
|
-
|
?
|
?
|
|
MomentumForcing
|
.TRUE.
cosine
wind
0.1
Pa
|
.FALSE.
|
.TRUE.
cosine
wind
0.2
Pa
|
.TRUE.
wind
data
|
|
FreeslipSide
FreeslipBottom
FreeslipTop
|
.TRUE.
.TRUE.
.TRUE.
|
?
?
?
|
?
?
?
|
?
?
?
|
|
PressureStepping
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
SeparatePh
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
phSurfBC
|
'zero'
|
'zero'
|
'zero'
|
'zero'
|
|
SeparatePs
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
rigid
lid
|
yes
|
yes
|
yes
|
no
|
|
toler2d
max2dIt
FreqCheckToler2d
divHmax
divHmin
|
|
|
|
|
|
stericEffect
|
.FALSE.
|
.FALSE.
|
.FALSE.
|
?
|
|
rhorefsurf
|
-
|
-
|
-
|
?
|
|
toler3d
max3dIt
FreqCheckToler3d
divmax
divmin
|
-
|
|
-
|
-
|
|
findD2pnhDZ2
|
-
|
?
|
-
|
-
|
|
Initial
condition
|
see ths
|
see ths
|
see ths
|
Levitus
data
|
|
TempStepping
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
TempAdvection
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
TradvecScheme
|
'centered'
|
'centered'
|
'centered'
|
'centered'
|
|
TempDiffusion
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
.TRUE.
|
|
TrdifScheme
|
'constant'
|
'constant'
|
'constant'
|
'GM'
|
|
CAPTH
RKPTH (m2 s-1)
|
?
?
|
5
0.2
|
0
10-5
|
-
|
|
tempBhdiffusion
|
.FALSE.
|
.FALSE.
|
.TRUE.
|
.FALSE.
|
|
C4PTH (m4s-1)
|
-
|
-
|
1010
|
-
|
|
tempForcing
|
.FALSE.
|
.TRUE.
constant
cooling
800
W m-2
|
.FALSE.
|
.TRUE.
heat
and
freshwater
fluxes
data
|
|
convectiveAdjustment
|
?
|
?
|
?
|
.TRUE.
|
|
cadjFreq
|
?
|
?
|
?
|
?
|
|
Idelt (s)
|
1200
|
60
|
?
|
3600
|
Not
written yet
Adcroft, A. (1995) Numerical
algorithms for use in a dynamical model of the ocean. Ph.D. thesis, Imperial College,
London.
Adcroft, A., Hill, C., and
Marshall, J. (1997) Representation of topography by shaved cells in a height
coordinate ocean model. Mon. Weather Rev.,
125, 2293-2315.
Arakawa, A., and Lamb, V.
(1977) Computational design of the basic dynamical processes of the UCLA
General Circulation Model. Methods
Comput. Phys., 17, 174-267.
Bryan, K. (1969) A numerical
model for the study of the circulation of the world ocean. J. Comput. Phys., 4,
347-376.
Dukowicz, J. K., Smith, R. D.,
and Malone, R. C. (1993) A reformulation and implementation of the
Bryan-Cox-Semtner ocean model on the connection machine. J. Atmos. Oceanic Technol., 10,
195-208.
Dukowicz, J. K. (1995) I don't have this one!
Harlow, F. H., and Welch, J. E.
(1965) Numerical calculation of time-dependent viscous incompressible flow of
fluid with free surface. Phys. Fluids,
8, 2182-2189.
Holland, W. R., and Lin, L. B.
(1975) On the generation of mesoscale eddies and their contribution to the
oceanic general circulation I. A preliminary numerical experiment. J. Phys. Oceanogr., 5, 642-657.
Jones, H., and Marshall, J.
(1993) Convection with rotation in a neutral ocean: a study of open-ocean deep
convection. J. Phys. Oceanogr., 23, 1009-1039.
Marshall, J., Hill, C.,
Perelman, L., and Adcroft, A. (1997a) Hydrostatic, quasi-hydrostatic, and
nonhydrostatic ocean modeling. J.
Geophys. Res., 102, 5733-5752.
Marshall, J., Adcroft, A.,
Hill, C., Perelman, L., and Heisey, C. (1997b) A finite volume, incompressible
Navier Stokes model for studies of the ocean on parallel computers. J. Geophys. Res., 102, 5753-5766.
Potter, D. (1976) Computational Physics, John Wiley, New
York.
Press, W. H., Flannery, B. P.,
Teukolsky, S. A., and Vettering, W. T. (1986) Numerical Recipes, Cambridge University Press, New York.
Visbeck, M., Marshall, J.,
Haine, T., and Spall, M. (1997) Specification of eddy transfer coefficients in
coarse-resolution ocean circulation models. J.
Phys. Oceanogr., 27, 381-402.
White, A. A., and Bromley, R.
A. (1995) Dynamically consistent, quasi-hydrostatic equations for global models
with a complete representation of the Coriolis force. Q. J. R. Meteorol. Soc., 121,
399-418.
Williams, G. P. (1969)
Numerical integrations of the three-dimensional Navier-Stokes equations for
incompressible flow. J. Fluid Mech., 37, 727-750.
has five diagonals corresponding to the coupling of the
central point with surrounding points along the four ‘arms’ of the horizontal 
; the search directions on consecutive iterations are
conjugate to one another.