s_examples/baroclinic_gyre/fourlayer.tex

% $Header: /u/gcmpack/manual/part3/case_studies/fourlayer_gyre/fourlayer.tex,v 1.19 2004/10/16 03:40:13 edhill Exp $
% $Name:  $

\section[Baroclinic Gyre MITgcm Example]{Four Layer Baroclinic Ocean Gyre In Spherical Coordinates}
\label{www:tutorials}
\label{sect:eg-fourlayer}
\begin{rawhtml}
<!-- CMIREDIR:eg-fourlayer: -->
\end{rawhtml}

\bodytext{bgcolor="#FFFFFFFF"}

%\begin{center} 
%{\Large \bf Using MITgcm to Simulate a Baroclinic Ocean Gyre In Spherical
%Polar Coordinates}
%
%\vspace*{4mm}
%
%\vspace*{3mm}
%{\large May 2001}
%\end{center}

This document describes an example experiment using MITgcm
to simulate a baroclinic ocean gyre in spherical
polar coordinates. The barotropic
example experiment in section \ref{sect:eg-baro}
illustrated how to configure the code for a single layer 
simulation in a Cartesian grid. In this example a similar physical problem
is simulated, but the code is now configured
for four layers and in a spherical polar coordinate system.

\subsection{Overview}
\label{www:tutorials}

This example experiment demonstrates using the MITgcm to simulate
a baroclinic, wind-forced, ocean gyre circulation. The experiment 
is a numerical rendition of the gyre circulation problem similar
to the problems described analytically by Stommel in 1966 
\cite{Stommel66} and numerically in Holland et. al \cite{Holland75}.
\\

In this experiment the model is configured to represent a mid-latitude 
enclosed sector of fluid on a sphere, $60^{\circ} \times 60^{\circ}$ in 
lateral extent. The fluid is $2$~km deep and is forced
by a constant in time zonal wind stress, $\tau_{\lambda}$, that varies 
sinusoidally in the north-south direction. Topologically the simulated 
domain is a sector on a sphere and the coriolis parameter, $f$, is defined 
according to latitude, $\varphi$

\begin{equation}
\label{EQ:eg-fourlayer-fcori}
f(\varphi) = 2 \Omega \sin( \varphi )
\end{equation}
 
\noindent with the rotation rate, $\Omega$ set to $\frac{2 \pi}{86400s}$.
\\
  
 The sinusoidal wind-stress variations are defined according to 

\begin{equation}
\label{EQ:taux}
\tau_{\lambda}(\varphi) = \tau_{0}\sin(\pi \frac{\varphi}{L_{\varphi}})
\end{equation}
 
\noindent where $L_{\varphi}$ is the lateral domain extent ($60^{\circ}$) and 
$\tau_0$ is set to $0.1N m^{-2}$. 
\\

Figure \ref{FIG:eg-fourlayer-simulation_config}
summarizes the configuration simulated.
In contrast to the example in section \ref{sect:eg-baro}, the 
current experiment simulates a spherical polar domain. As indicated
by the axes in the lower left of the figure the model code works internally
in a locally orthogonal coordinate $(x,y,z)$. For this experiment description 
the local orthogonal model coordinate $(x,y,z)$ is synonymous 
with the coordinates $(\lambda,\varphi,r)$ shown in figure
\ref{fig:spherical-polar-coord}
\\

The experiment has four levels in the vertical, each of equal thickness,
$\Delta z = 500$~m. Initially the fluid is stratified with a reference
potential temperature profile,
$\theta_{250}=20^{\circ}$~C,
$\theta_{750}=10^{\circ}$~C,
$\theta_{1250}=8^{\circ}$~C,
$\theta_{1750}=6^{\circ}$~C. The equation of state used in this experiment is 
linear

\begin{equation}
\label{EQ:eg-fourlayer-linear1_eos}
\rho = \rho_{0} ( 1 - \alpha_{\theta}\theta^{'} )
\end{equation}

\noindent which is implemented in the model as a density anomaly equation

\begin{equation}
\label{EQ:eg-fourlayer-linear1_eos_pert}
\rho^{'} = -\rho_{0}\alpha_{\theta}\theta^{'}
\end{equation}

\noindent with $\rho_{0}=999.8\,{\rm kg\,m}^{-3}$ and 
$\alpha_{\theta}=2\times10^{-4}\,{\rm degrees}^{-1} $. Integrated forward in
this configuration the model state variable {\bf theta} is equivalent to
either in-situ temperature, $T$, or potential temperature, $\theta$. For 
consistency with later examples, in which the equation of state is
non-linear, we use $\theta$ to represent temperature here. This is
the quantity that is carried in the model core equations.

\begin{figure}
%% \begin{center}
%%  \resizebox{7.5in}{5.5in}{
%%    \includegraphics*[0.2in,0.7in][10.5in,10.5in]
%%    {part3/case_studies/fourlayer_gyre/simulation_config.eps} }
%% \end{center}
\centerline{
  \scalefig{.95}
  \epsfbox{part3/case_studies/fourlayer_gyre/simulation_config.eps}
}
\caption{Schematic of simulation domain and wind-stress forcing function 
for the four-layer gyre numerical experiment. The domain is enclosed by solid
walls at $0^{\circ}$~E, $60^{\circ}$~E, $0^{\circ}$~N and $60^{\circ}$~N.
An initial stratification is 
imposed by setting the potential temperature, $\theta$, in each layer.
The vertical spacing, $\Delta z$, is constant and equal to $500$m.
}
\label{FIG:eg-fourlayer-simulation_config}
\end{figure}

\subsection{Equations solved}
\label{www:tutorials}
For this problem
the implicit free surface, {\bf HPE} (see section \ref{sect:hydrostatic_and_quasi-hydrostatic_forms}) form of the 
equations described in Marshall et. al \cite{marshall:97a} are
employed. The flow is three-dimensional with just temperature, $\theta$, as 
an active tracer.  The equation of state is linear.
A horizontal Laplacian operator $\nabla_{h}^2$ provides viscous
dissipation and provides a diffusive sub-grid scale closure for the 
temperature equation. A wind-stress momentum forcing is added to the momentum 
equation for the zonal flow, $u$. Other terms in the model
are explicitly switched off for this experiment configuration (see section
\ref{SEC:eg_fourl_code_config} ). This yields an active set of equations
solved in this configuration, written in spherical polar coordinates as 
follows

\begin{eqnarray}
\label{EQ:eg-fourlayer-model_equations}
\frac{Du}{Dt} - fv + 
  \frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \lambda} - 
  A_{h}\nabla_{h}^2u - A_{z}\frac{\partial^{2}u}{\partial z^{2}} 
& = &
\cal{F}_{\lambda}
\\
\frac{Dv}{Dt} + fu + 
  \frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \varphi} - 
  A_{h}\nabla_{h}^2v - A_{z}\frac{\partial^{2}v}{\partial z^{2}} 
& = &
0
\\
\frac{\partial \eta}{\partial t} + \frac{\partial H \widehat{u}}{\partial \lambda} +
\frac{\partial H \widehat{v}}{\partial \varphi}
&=&
0
\label{eq:fourl_example_continuity}
\\
\frac{D\theta}{Dt} -
 K_{h}\nabla_{h}^2\theta  - K_{z}\frac{\partial^{2}\theta}{\partial z^{2}} 
& = &
0
\label{eq:eg_fourl_theta}
\\
p^{\prime} & = &
g\rho_{0} \eta + \int^{0}_{-z}\rho^{\prime} dz
\\
\rho^{\prime} & = &- \alpha_{\theta}\rho_{0}\theta^{\prime}
\\
{\cal F}_{\lambda} |_{s} & = & \frac{\tau_{\lambda}}{\rho_{0}\Delta z_{s}}
\\
{\cal F}_{\lambda} |_{i} & = & 0
\end{eqnarray}

\noindent where $u$ and $v$ are the components of the horizontal
flow vector $\vec{u}$ on the sphere ($u=\dot{\lambda},v=\dot{\varphi}$).
The terms $H\widehat{u}$ and $H\widehat{v}$ are the components of the vertical
integral term given in equation \ref{eq:free-surface} and
explained in more detail in section \ref{sect:pressure-method-linear-backward}.
However, for the problem presented here, the continuity relation (equation
\ref{eq:fourl_example_continuity}) differs from the general form given
in section \ref{sect:pressure-method-linear-backward},
equation \ref{eq:linear-free-surface=P-E+R}, because the source terms
${\cal P}-{\cal E}+{\cal R}$ 
are all $0$.

The pressure field, $p^{\prime}$, is separated into a barotropic part
due to variations in sea-surface height, $\eta$, and a hydrostatic
part due to variations in density, $\rho^{\prime}$, integrated
through the water column.

The suffices ${s},{i}$ indicate surface layer and the interior of the domain.
The windstress forcing, ${\cal F}_{\lambda}$, is applied in the surface layer 
by a source term in the zonal momentum equation. In the ocean interior
this term is zero.

In the momentum equations
lateral and vertical boundary conditions for the $\nabla_{h}^{2}$
and $\frac{\partial^{2}}{\partial z^{2}}$ operators are specified
when the numerical simulation is run - see section 
\ref{SEC:eg_fourl_code_config}. For temperature
the boundary condition is ``zero-flux'' 
e.g. $\frac{\partial \theta}{\partial \varphi}=
\frac{\partial \theta}{\partial \lambda}=\frac{\partial \theta}{\partial z}=0$.


\subsection{Discrete Numerical Configuration}
\label{www:tutorials}

 The domain is discretised with 
a uniform grid spacing in latitude and longitude
 $\Delta \lambda=\Delta \varphi=1^{\circ}$, so 
that there are sixty grid cells in the zonal and meridional directions. 
Vertically the 
model is configured with four layers with constant depth, 
$\Delta z$, of $500$~m. The internal, locally orthogonal, model coordinate 
variables $x$ and $y$ are initialized from the values of
$\lambda$, $\varphi$, $\Delta \lambda$ and $\Delta \varphi$ in
radians according to

\begin{eqnarray}
x=r\cos(\varphi)\lambda,~\Delta x & = &r\cos(\varphi)\Delta \lambda \\
y=r\varphi,~\Delta y &= &r\Delta \varphi
\end{eqnarray}

The procedure for generating a set of internal grid variables from a
spherical polar grid specification is discussed in section 
\ref{sect:spatial_discrete_horizontal_grid}.

\noindent\fbox{ \begin{minipage}{5.5in}
{\em S/R INI\_SPHERICAL\_POLAR\_GRID} ({\em
model/src/ini\_spherical\_polar\_grid.F})

$A_c$, $A_\zeta$, $A_w$, $A_s$: {\bf rAc}, {\bf rAz}, {\bf rAw}, {\bf rAs}
({\em GRID.h})

$\Delta x_g$, $\Delta y_g$: {\bf DXg}, {\bf DYg} ({\em GRID.h})

$\Delta x_c$, $\Delta y_c$: {\bf DXc}, {\bf DYc} ({\em GRID.h})

$\Delta x_f$, $\Delta y_f$: {\bf DXf}, {\bf DYf} ({\em GRID.h})

$\Delta x_v$, $\Delta y_u$: {\bf DXv}, {\bf DYu} ({\em GRID.h})

\end{minipage} }\\


As described in \ref{sect:tracer_equations}, the time evolution of potential 
temperature, 
$\theta$, (equation \ref{eq:eg_fourl_theta})
is evaluated prognostically. The centered second-order scheme with
Adams-Bashforth time stepping described in section 
\ref{sect:tracer_equations_abII} is used to step forward the temperature 
equation. Prognostic terms in
the momentum equations are solved using flux form as
described in section \ref{sect:flux-form_momentum_eqautions}.
The pressure forces that drive the fluid motions, (
$\frac{\partial p^{'}}{\partial \lambda}$ and $\frac{\partial p^{'}}{\partial \varphi}$), are found by summing pressure due to surface 
elevation $\eta$ and the hydrostatic pressure. The hydrostatic part of the 
pressure is diagnosed explicitly by integrating density. The sea-surface
height, $\eta$, is diagnosed using an implicit scheme. The pressure
field solution method is described in sections
\ref{sect:pressure-method-linear-backward} and 
\ref{sect:finding_the_pressure_field}.

\subsubsection{Numerical Stability Criteria}
\label{www:tutorials}

The Laplacian viscosity coefficient, $A_{h}$, is set to $400 m s^{-1}$.
This value is chosen to yield a Munk layer width,

\begin{eqnarray}
\label{EQ:eg-fourlayer-munk_layer}
M_{w} = \pi ( \frac { A_{h} }{ \beta } )^{\frac{1}{3}}
\end{eqnarray}

\noindent  of $\approx 100$km. This is greater than the model
resolution in mid-latitudes 
$\Delta x=r \cos(\varphi) \Delta \lambda \approx 80~{\rm km}$ at
$\varphi=45^{\circ}$, ensuring that the frictional 
boundary layer is well resolved.
\\

\noindent The model is stepped forward with a 
time step $\delta t=1200$secs. With this time step the stability 
parameter to the horizontal Laplacian friction

\begin{eqnarray}
\label{EQ:eg-fourlayer-laplacian_stability}
S_{l} = 4 \frac{A_{h} \delta t}{{\Delta x}^2}
\end{eqnarray}

\noindent evaluates to 0.012, which is well below the 0.3 upper limit
for stability for this term under ABII time-stepping.
\\

\noindent The vertical dissipation coefficient, $A_{z}$, is set to 
$1\times10^{-2} {\rm m}^2{\rm s}^{-1}$. The associated stability limit

\begin{eqnarray}
\label{EQ:eg-fourlayer-laplacian_stability_z}
S_{l} = 4 \frac{A_{z} \delta t}{{\Delta z}^2}
\end{eqnarray}

\noindent evaluates to $4.8 \times 10^{-5}$ which is again well below
the upper limit.
The values of $A_{h}$ and $A_{z}$ are also used for the horizontal ($K_{h}$) 
and vertical ($K_{z}$) diffusion coefficients for temperature respectively.
\\

\noindent The numerical stability for inertial oscillations

\begin{eqnarray}
\label{EQ:eg-fourlayer-inertial_stability}
S_{i} = f^{2} {\delta t}^2
\end{eqnarray}

\noindent evaluates to $0.0144$, which is well below the $0.5$ upper 
limit for stability.
\\

\noindent The advective CFL for a extreme maximum 
horizontal flow
speed of $ | \vec{u} | = 2 ms^{-1}$

\begin{eqnarray}
\label{EQ:eg-fourlayer-cfl_stability}
C_{a} = \frac{| \vec{u} | \delta t}{ \Delta x}
\end{eqnarray}

\noindent evaluates to $5 \times 10^{-2}$. This is well below the stability 
limit of 0.5.
\\

\noindent The stability parameter for internal gravity waves
propagating at $2~{\rm m}~{\rm s}^{-1}$ 

\begin{eqnarray}
\label{EQ:eg-fourlayer-igw_stability}
S_{c} = \frac{c_{g} \delta t}{ \Delta x}
\end{eqnarray}

\noindent evaluates to $\approx 5 \times 10^{-2}$. This is well below the linear
stability limit of 0.25.
  
\subsection{Code Configuration}
\label{www:tutorials}
\label{SEC:eg_fourl_code_config}

The model configuration for this experiment resides under the 
directory {\it verification/exp2/}.  The experiment files 
\begin{itemize}
\item {\it input/data}
\item {\it input/data.pkg}
\item {\it input/eedata},
\item {\it input/windx.sin\_y},
\item {\it input/topog.box},
\item {\it code/CPP\_EEOPTIONS.h}
\item {\it code/CPP\_OPTIONS.h},
\item {\it code/SIZE.h}. 
\end{itemize}
contain the code customisations and parameter settings for this
experiment. Below we describe the customisations to these files
associated with this experiment.

\subsubsection{File {\it input/data}}
\label{www:tutorials}

This file, reproduced completely below, specifies the main parameters 
for the experiment. The parameters that are significant for this configuration
are

\begin{itemize}

\item Line 4, 
\begin{verbatim} tRef=20.,10.,8.,6., \end{verbatim} 
this line sets the initial and reference values of potential
temperature at each model level in units of $^{\circ}\mathrm{C}$.  The entries
are ordered from surface to depth. For each depth level the initial
and reference profiles will be uniform in $x$ and $y$. The values
specified here are read into the variable \varlink{tRef}{tRef} in the
model code, by procedure \filelink{INI\_PARMS}{model-src-ini_parms.F}

\fbox{
  \begin{minipage}{5.0in}
    {\it S/R INI\_THETA}({\it ini\_theta.F})
  \end{minipage}
}
\filelink{ini\_theta.F}{model-src-ini_theta.F}

\item Line 6, 
\begin{verbatim} viscAz=1.E-2, \end{verbatim} 
this line sets the vertical Laplacian dissipation coefficient to $1
\times 10^{-2} {\rm m^{2}s^{-1}}$. Boundary conditions for this
operator are specified later.  The variable \varlink{viscAz}{viscAz}
is read in the routine \filelink{ini\_parms.F}{model-src-ini_parms.F}
and is copied into model general vertical coordinate variable
\varlink{viscAr}{viscAr} At each time step, the viscous term
contribution to the momentum equations is calculated in routine
\varlink{CALC\_DIFFUSIVITY}{CALC_DIFFUSIVITY}

\fbox{
\begin{minipage}{5.0in}
{\it S/R CALC\_DIFFUSIVITY}({\it calc\_diffusivity.F})
\end{minipage}
}

\item Line 7, 
\begin{verbatim}
viscAh=4.E2,
\end{verbatim} 
  this line sets the horizontal laplacian frictional dissipation
  coefficient to $1 \times 10^{-2} {\rm m^{2}s^{-1}}$. Boundary
  conditions for this operator are specified later.  The variable
  \varlink{viscAh}{viscAh} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS} and applied in routines
  \varlink{CALC\_MOM\_RHS}{CALC_MOM_RHS} and
  \varlink{CALC\_GW}{CALC_GW}.

\fbox{
  \begin{minipage}{5.0in}
    {\it S/R CALC\_MOM\_RHS}({\it calc\_mom\_rhs.F})
  \end{minipage}
}
\fbox{
  \begin{minipage}{5.0in}
    {\it S/R CALC\_GW}({\it calc\_gw.F})
  \end{minipage}
}

\item Line 8,
\begin{verbatim}
no_slip_sides=.FALSE.
\end{verbatim}
  this line selects a free-slip lateral boundary condition for the
  horizontal laplacian friction operator e.g. $\frac{\partial
    u}{\partial y}$=0 along boundaries in $y$ and $\frac{\partial
    v}{\partial x}$=0 along boundaries in $x$.  The variable
  \varlink{no\_slip\_sides}{no_slip_sides} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS} and the boundary condition is
  evaluated in routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R CALC\_MOM\_RHS}({\it calc\_mom\_rhs.F})
    \end{minipage}
  }
  \filelink{calc\_mom\_rhs.F}{calc_mom_rhs.F}
  
\item Lines 9,
\begin{verbatim}
no_slip_bottom=.TRUE.
\end{verbatim}
  this line selects a no-slip boundary condition for bottom boundary
  condition in the vertical laplacian friction operator e.g. $u=v=0$
  at $z=-H$, where $H$ is the local depth of the domain.  The variable
  \varlink{no\_slip\_bottom}{no\_slip\_bottom} is read in the routine
  \filelink{INI\_PARMS}{model-src-ini_parms.F} and is applied in the
  routine \varlink{CALC\_MOM\_RHS}{CALC_MOM_RHS}.

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R CALC\_MOM\_RHS}({\it calc\_mom\_rhs.F})
    \end{minipage}
  }
  \filelink{calc\_mom\_rhs.F}{calc_mom_rhs.F}

\item Line 10,
\begin{verbatim}
diffKhT=4.E2,
\end{verbatim}
  this line sets the horizontal diffusion coefficient for temperature
  to $400\,{\rm m^{2}s^{-1}}$. The boundary condition on this operator
  is $\frac{\partial}{\partial x}=\frac{\partial}{\partial y}=0$ at
  all boundaries.  The variable \varlink{diffKhT}{diffKhT} is read in
  the routine \varlink{INI\_PARMS}{INI_PARMS} and used in routine
  \varlink{CALC\_GT}{CALC_GT}.

  \fbox{ \begin{minipage}{5.0in}
      {\it S/R CALC\_GT}({\it calc\_gt.F})
    \end{minipage}
  }
  \filelink{calc\_gt.F}{model-src-calc_gt.F}

\item Line 11,
\begin{verbatim}
diffKzT=1.E-2,
\end{verbatim}
  this line sets the vertical diffusion coefficient for temperature to
  $10^{-2}\,{\rm m^{2}s^{-1}}$. The boundary condition on this
  operator is $\frac{\partial}{\partial z}$ = 0 on all boundaries.
  The variable \varlink{diffKzT}{diffKzT} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS}. It is copied into model general
  vertical coordinate variable \varlink{diffKrT}{diffKrT} which is
  used in routine \varlink{CALC\_DIFFUSIVITY}{CALC_DIFFUSIVITY}.

  \fbox{ \begin{minipage}{5.0in}
      {\it S/R CALC\_DIFFUSIVITY}({\it calc\_diffusivity.F})
    \end{minipage}
  }
  \filelink{calc\_diffusivity.F}{model-src-calc_diffusivity.F}

\item Line 13,
\begin{verbatim}
tAlpha=2.E-4,
\end{verbatim}
  This line sets the thermal expansion coefficient for the fluid to $2
  \times 10^{-4}\,{\rm degrees}^{-1}$ The variable
  \varlink{tAlpha}{tAlpha} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS}. The routine
  \varlink{FIND\_RHO}{FIND\_RHO} makes use of {\bf tAlpha}.

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R FIND\_RHO}({\it find\_rho.F})
    \end{minipage}
  }
  \filelink{find\_rho.F}{model-src-find_rho.F}

\item Line 18,
\begin{verbatim}
eosType='LINEAR'
\end{verbatim}
  This line selects the linear form of the equation of state.  The
  variable \varlink{eosType}{eosType} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS}. The values of {\bf eosType} sets
  which formula in routine {\it FIND\_RHO} is used to calculate
  density.

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R FIND\_RHO}({\it find\_rho.F})
    \end{minipage}
  }
  \filelink{find\_rho.F}{model-src-find_rho.F}

\item Line 40,
\begin{verbatim}
usingSphericalPolarGrid=.TRUE.,
\end{verbatim}
  This line requests that the simulation be performed in a spherical
  polar coordinate system. It affects the interpretation of grid input
  parameters, for example {\bf delX} and {\bf delY} and causes the
  grid generation routines to initialize an internal grid based on
  spherical polar geometry.  The variable
  \varlink{usingSphericalPolarGrid}{usingSphericalPolarGrid} is read
  in the routine \varlink{INI\_PARMS}{INI_PARMS}. When set to {\bf
    .TRUE.} the settings of {\bf delX} and {\bf delY} are taken to be
  in degrees. These values are used in the routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R INI\_SPEHRICAL\_POLAR\_GRID}({\it ini\_spherical\_polar\_grid.F})
    \end{minipage}
  }
  \filelink{ini\_spherical\_polar\_grid.F}{model-src-ini_spherical_polar_grid.F}

\item Line 41,
\begin{verbatim}
phiMin=0.,
\end{verbatim}
  This line sets the southern boundary of the modeled domain to
  $0^{\circ}$ latitude. This value affects both the generation of the
  locally orthogonal grid that the model uses internally and affects
  the initialization of the coriolis force.  Note - it is not required
  to set a longitude boundary, since the absolute longitude does not
  alter the kernel equation discretisation.  The variable
  \varlink{phiMin}{phiMin} is read in the
  routine \varlink{INI\_PARMS}{INI_PARMS} and is used in routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R INI\_SPEHRICAL\_POLAR\_GRID}({\it ini\_spherical\_polar\_grid.F})
    \end{minipage}
  }
  \filelink{ini\_spherical\_polar\_grid.F}{model-src-ini_spherical_polar_grid.F}

\item Line 42,
\begin{verbatim}
delX=60*1.,
\end{verbatim}
  This line sets the horizontal grid spacing between each y-coordinate
  line in the discrete grid to $1^{\circ}$ in longitude.  The variable
  \varlink{delX}{delX} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS} and is used in routine
 
  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R INI\_SPEHRICAL\_POLAR\_GRID}({\it ini\_spherical\_polar\_grid.F})
    \end{minipage}
  }
  \filelink{ini\_spherical\_polar\_grid.F}{model-src-ini_spherical_polar_grid.F}

\item Line 43,
\begin{verbatim}
delY=60*1.,
\end{verbatim}
  This line sets the horizontal grid spacing between each y-coordinate
  line in the discrete grid to $1^{\circ}$ in latitude.  The variable
  \varlink{delY}{delY} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS} and is used in routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R INI\_SPEHRICAL\_POLAR\_GRID}({\it ini\_spherical\_polar\_grid.F})
    \end{minipage}
  }
  \filelink{ini\_spherical\_polar\_grid.F}{model-src-ini_spherical_polar_grid.F}

\item Line 44,
\begin{verbatim}
delZ=500.,500.,500.,500.,
\end{verbatim}
  This line sets the vertical grid spacing between each z-coordinate
  line in the discrete grid to $500\,{\rm m}$, so that the total model
  depth is $2\,{\rm km}$.  The variable \varlink{delZ}{delZ} is read
  in the routine \varlink{INI\_PARMS}{INI_PARMS}.  It is copied into
  the internal model coordinate variable \varlink{delR}{delR} which is
  used in routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R INI\_VERTICAL\_GRID}({\it ini\_vertical\_grid.F})
    \end{minipage}
  }
  \filelink{ini\_vertical\_grid.F}{model-src-ini_vertical_grid.F}

\item Line 47,
\begin{verbatim}
bathyFile='topog.box'
\end{verbatim}
  This line specifies the name of the file from which the domain
  bathymetry is read. This file is a two-dimensional ($x,y$) map of
  depths. This file is assumed to contain 64-bit binary numbers giving
  the depth of the model at each grid cell, ordered with the x
  coordinate varying fastest. The points are ordered from low
  coordinate to high coordinate for both axes. The units and
  orientation of the depths in this file are the same as used in the
  MITgcm code. In this experiment, a depth of $0m$ indicates a solid
  wall and a depth of $-2000m$ indicates open ocean. The matlab
  program {\it input/gendata.m} shows an example of how to generate a
  bathymetry file.  The variable \varlink{bathyFile}{bathyFile} is
  read in the routine \varlink{INI\_PARMS}{INI_PARMS}.  The bathymetry
  file is read in the routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R INI\_DEPTHS}({\it ini\_depths.F})
    \end{minipage}
  }
  \filelink{ini\_depths.F}{model-src-ini_depths.F}

\item Line 50,
\begin{verbatim}
zonalWindFile='windx.sin_y'
\end{verbatim}
  This line specifies the name of the file from which the x-direction
  (zonal) surface wind stress is read. This file is also a
  two-dimensional ($x,y$) map and is enumerated and formatted in the
  same manner as the bathymetry file. The matlab program {\it
    input/gendata.m} includes example code to generate a valid {\bf
    zonalWindFile} file.  The variable
  \varlink{zonalWindFile}{zonalWindFile} is read in the routine
  \varlink{INI\_PARMS}{INI_PARMS}.  The wind-stress file is read in
  the routine

  \fbox{
    \begin{minipage}{5.0in}
      {\it S/R EXTERNAL\_FIELDS\_LOAD}({\it external\_fields\_load.F})
    \end{minipage}
  }
  \filelink{external\_fields\_load.F}{model-src-external_fields_load.F}

\end{itemize}

\noindent other lines in the file {\it input/data} are standard values.

\begin{rawhtml}<PRE>\end{rawhtml}
\begin{small}
\input{part3/case_studies/fourlayer_gyre/input/data}
\end{small}
\begin{rawhtml}</PRE>\end{rawhtml}

\subsubsection{File {\it input/data.pkg}}
\label{www:tutorials}

This file uses standard default values and does not contain
customisations for this experiment.

\subsubsection{File {\it input/eedata}}
\label{www:tutorials}

This file uses standard default values and does not contain
customisations for this experiment.

\subsubsection{File {\it input/windx.sin\_y}}
\label{www:tutorials}

The {\it input/windx.sin\_y} file specifies a two-dimensional ($x,y$)
map of wind stress ,$\tau_{x}$, values. The units used are $Nm^{-2}$
(the default for MITgcm).  Although $\tau_{x}$ is only a function of
latitude, $y$, in this experiment this file must still define a
complete two-dimensional map in order to be compatible with the
standard code for loading forcing fields in MITgcm (routine {\it
  EXTERNAL\_FIELDS\_LOAD}.  The included matlab program {\it
  input/gendata.m} gives a complete code for creating the {\it
  input/windx.sin\_y} file.

\subsubsection{File {\it input/topog.box}}
\label{www:tutorials}


The {\it input/topog.box} file specifies a two-dimensional ($x,y$) 
map of depth values. For this experiment values are either
$0~{\rm m}$ or $-2000\,{\rm m}$, corresponding respectively to a wall or to deep
ocean. The file contains a raw binary stream of data that is enumerated
in the same way as standard MITgcm two-dimensional, horizontal arrays.
The included matlab program {\it input/gendata.m} gives a complete
code for creating the {\it input/topog.box} file.

\subsubsection{File {\it code/SIZE.h}}
\label{www:tutorials}

Two lines are customized in this file for the current experiment

\begin{itemize}

\item Line 39, 
\begin{verbatim} sNx=60, \end{verbatim} this line sets
the lateral domain extent in grid points for the
axis aligned with the x-coordinate.

\item Line 40, 
\begin{verbatim} sNy=60, \end{verbatim} this line sets
the lateral domain extent in grid points for the
axis aligned with the y-coordinate.

\item Line 49, 
\begin{verbatim} Nr=4,   \end{verbatim} this line sets
the vertical domain extent in grid points.

\end{itemize}

\begin{small}
\include{part3/case_studies/fourlayer_gyre/code/SIZE.h}
\end{small}

\subsubsection{File {\it code/CPP\_OPTIONS.h}}
\label{www:tutorials}

This file uses standard default values and does not contain
customisations for this experiment.


\subsubsection{File {\it code/CPP\_EEOPTIONS.h}}
\label{www:tutorials}

This file uses standard default values and does not contain
customisations for this experiment.

\subsubsection{Other Files }
\label{www:tutorials}

Other files relevant to this experiment are
\begin{itemize}
\item {\it model/src/ini\_cori.F}. This file initializes the model
coriolis variables {\bf fCorU} and {\bf fCorV}.
\item {\it model/src/ini\_spherical\_polar\_grid.F} This file
initializes the model grid discretisation variables {\bf
dxF, dyF, dxG, dyG, dxC, dyC}.
\item {\it model/src/ini\_parms.F}.
\end{itemize}

\subsection{Running The Example}
\label{www:tutorials}
\label{SEC:running_the_example}

\subsubsection{Code Download}
\label{www:tutorials}

 In order to run the examples you must first download the code distribution.
Instructions for downloading the code can be found in section
\ref{sect:obtainingCode}.

\subsubsection{Experiment Location}
\label{www:tutorials}

 This example experiments is located under the release sub-directory

\vspace{5mm}
{\it verification/exp2/ }

\subsubsection{Running the Experiment}
\label{www:tutorials}

 To run the experiment

\begin{enumerate}
\item Set the current directory to {\it input/ }

\begin{verbatim}
% cd input
\end{verbatim}

\item Verify that current directory is now correct

\begin{verbatim}
% pwd
\end{verbatim}

 You should see a response on the screen ending in

{\it verification/exp2/input }


\item Run the genmake script to create the experiment {\it Makefile}

\begin{verbatim}
% ../../../tools/genmake -mods=../code
\end{verbatim}

\item Create a list of header file dependencies in {\it Makefile}

\begin{verbatim}
% make depend
\end{verbatim}

\item Build the executable file.

\begin{verbatim}
% make
\end{verbatim}

\item Run the {\it mitgcmuv} executable

\begin{verbatim}
% ./mitgcmuv
\end{verbatim}

\end{enumerate}

