--- manual/s_examples/baroclinic_gyre/fourlayer.tex	2001/10/24 19:43:07	1.4
+++ manual/s_examples/baroclinic_gyre/fourlayer.tex	2003/07/30 13:42:53	1.14
@@ -1,8 +1,9 @@
-% $Header: /home/ubuntu/mnt/e9_copy/manual/s_examples/baroclinic_gyre/fourlayer.tex,v 1.4 2001/10/24 19:43:07 cnh Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/manual/s_examples/baroclinic_gyre/fourlayer.tex,v 1.14 2003/07/30 13:42:53 edhill Exp $
 % $Name:  $
 
-\section{Example: Four layer Baroclinic Ocean Gyre In Spherical Coordinates}
-\label{sec:eg-fourlayer}
+\section{Four Layer Baroclinic Ocean Gyre In Spherical Coordinates}
+\label{www:tutorials}
+\label{sect:eg-fourlayer}
 
 \bodytext{bgcolor="#FFFFFFFF"}
 
@@ -19,17 +20,18 @@
 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{sec:eg-baro}
-ilustrated how to configure the code for a single layer 
-simulation in a cartesian grid. In this example a similar physical problem
+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 simliar
+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}.
 \\
@@ -37,13 +39,13 @@
 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_x$, that varies sinusoidally
-in the north-south direction. Topologically the simulated 
+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:fcori}
+\label{EQ:eg-fourlayer-fcori}
 f(\varphi) = 2 \Omega \sin( \varphi )
 \end{equation}
  
@@ -54,21 +56,21 @@
 
 \begin{equation}
 \label{EQ:taux}
-\tau_x(\varphi) = \tau_{0}\sin(\pi \frac{\varphi}{L_{\varphi}})
+\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:simulation_config}
-summarises the configuration simulated.
-In contrast to the example in section \ref{sec:eg-baro}, the 
-current experiment simulates a spherical polar domain. However, as indicated
+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 orthoganal coordinate $(x,y,z)$. For this experiment description 
-of this document the local orthogonal model coordinate $(x,y,z)$ is synonomous 
-with the spherical polar coordinate shown in figure 
+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}
 \\
 
@@ -82,20 +84,20 @@
 linear
 
 \begin{equation}
-\label{EQ:linear1_eos}
+\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:linear1_eos_pert}
+\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 synonomous with
+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
@@ -110,33 +112,36 @@
 \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.
-In the four-layer case an initial temperature stratification is 
+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:simulation_config}
+\label{FIG:eg-fourlayer-simulation_config}
 \end{figure}
 
 \subsection{Equations solved}
-
-The implicit free surface form of the 
-pressure equation described in Marshall et. al \cite{Marshall97a} is 
-employed. 
-A horizontal laplacian operator $\nabla_{h}^2$ provides viscous
-dissipation. The wind-stress momentum input is added to the momentum equation
-for the ``zonal flow'', $u$. Other terms in the model
-are explicitly switched off for this experiement configuration (see section
-\ref{SEC:code_config} ). This yields an active set of equations in 
+\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:model_equations}
+\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}
+\cal{F}_{\lambda}
 \\
 \frac{Dv}{Dt} + fu + 
   \frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \varphi} - 
@@ -144,47 +149,72 @@
 & = &
 0
 \\
-\frac{\partial \eta}{\partial t} + \frac{\partial H \hat{u}}{\partial \lambda} +
-\frac{\partial H \hat{v}}{\partial \varphi}
+\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} |_{s} & = & \frac{\tau_{x}}{\rho_{0}\Delta z_{s}}
+{\cal F}_{\lambda} |_{s} & = & \frac{\tau_{\lambda}}{\rho_{0}\Delta z_{s}}
 \\
-{\cal F} |_{i} & = & 0
+{\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\hat{u}$ and $H\hat{v}$ are the components of the term
-integrated in eqaution \ref{eq:free-surface}, as descirbed in section
+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 suffices ${s},{i}$ indicate surface and interior of the domain.
-The forcing $\cal F$ is only applied at the surface.
 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}$, over the water column.
+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 model is configured in hydrostatic form.  The domain is discretised with 
+ 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 initialised from the values of
+variables $x$ and $y$ are initialized from the values of
 $\lambda$, $\varphi$, $\Delta \lambda$ and $\Delta \varphi$ in
 radians according to
 
@@ -195,7 +225,7 @@
 
 The procedure for generating a set of internal grid variables from a
 spherical polar grid specification is discussed in section 
-\ref{sec:spatial_discrete_horizontal_grid}.
+\ref{sect:spatial_discrete_horizontal_grid}.
 
 \noindent\fbox{ \begin{minipage}{5.5in}
 {\em S/R INI\_SPHERICAL\_POLAR\_GRID} ({\em
@@ -216,46 +246,60 @@
 
 
 
-As described in \ref{sec:tracer_equations}, the time evolution of potential 
+As described in \ref{sect:tracer_equations}, the time evolution of potential 
 temperature, 
-$\theta$, equation is solved prognostically.
+$\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.
+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 dissipation coefficient, $A_{h}$, is set to $400 m s^{-1}$.
-This value is chosen to yield a Munk layer width \cite{Adcroft_thesis},
+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:munk_layer}
+\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$, ensuring that the frictional 
+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 \cite{Adcroft_thesis}
+parameter to the horizontal Laplacian friction
 
 \begin{eqnarray}
-\label{EQ:laplacian_stability}
+\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 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:laplacian_stability_z}
+\label{EQ:eg-fourlayer-laplacian_stability_z}
 S_{l} = 4 \frac{A_{z} \delta t}{{\Delta z}^2}
 \end{eqnarray}
 
@@ -266,10 +310,9 @@
 \\
 
 \noindent The numerical stability for inertial oscillations
-\cite{Adcroft_thesis} 
 
 \begin{eqnarray}
-\label{EQ:inertial_stability}
+\label{EQ:eg-fourlayer-inertial_stability}
 S_{i} = f^{2} {\delta t}^2
 \end{eqnarray}
 
@@ -277,35 +320,36 @@
 limit for stability.
 \\
 
-\noindent The advective CFL \cite{Adcroft_thesis} for a extreme maximum 
+\noindent The advective CFL for a extreme maximum 
 horizontal flow
 speed of $ | \vec{u} | = 2 ms^{-1}$
 
 \begin{eqnarray}
-\label{EQ:cfl_stability}
-S_{a} = \frac{| \vec{u} | \delta t}{ \Delta x}
+\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 
-\cite{Adcroft_thesis}
+\noindent The stability parameter for internal gravity waves
+propagating at $2~{\rm m}~{\rm s}^{-1}$ 
 
 \begin{eqnarray}
-\label{EQ:igw_stability}
+\label{EQ:eg-fourlayer-igw_stability}
 S_{c} = \frac{c_{g} \delta t}{ \Delta x}
 \end{eqnarray}
 
-\noindent evaluates to $5 \times 10^{-2}$. This is well below the linear
+\noindent evaluates to $\approx 5 \times 10^{-2}$. This is well below the linear
 stability limit of 0.25.
   
 \subsection{Code Configuration}
-\label{SEC:code_config}
+\label{www:tutorials}
+\label{SEC:eg_fourl_code_config}
 
 The model configuration for this experiment resides under the 
-directory {\it verification/exp1/}.  The experiment files 
+directory {\it verification/exp2/}.  The experiment files 
 \begin{itemize}
 \item {\it input/data}
 \item {\it input/data.pkg}
@@ -317,10 +361,11 @@
 \item {\it code/SIZE.h}. 
 \end{itemize}
 contain the code customisations and parameter settings for this 
-experiements. Below we describe the customisations
+experiments. 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
@@ -334,17 +379,17 @@
 the initial and reference values of potential temperature at each model
 level in units of $^{\circ}$C.
 The entries are ordered from surface to depth. For each
-depth level the inital and reference profiles will be uniform in
+depth level the initial and reference profiles will be uniform in
 $x$ and $y$. The values specified here are read into the
 variable 
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/OK.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/OK.htm> \end{rawhtml}
 tRef
 \begin{rawhtml} </A>\end{rawhtml}
 } 
 in the model code, by procedure 
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
 }.
@@ -357,7 +402,7 @@
 %% \file{ini_parms}
 \newcommand{\VARtref}{
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/OK.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/OK.htm> \end{rawhtml}
 tRef
 \begin{rawhtml} </A>\end{rawhtml}
 } 
@@ -372,7 +417,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/98.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/98.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -380,27 +425,29 @@
 
 \item Line 6, 
 \begin{verbatim} viscAz=1.E-2, \end{verbatim} 
-this line sets the vertical laplacian dissipation coefficient to
+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 
 {\bf 
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/ZQ.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/ZQ.htm> \end{rawhtml}
 viscAz
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
 }
 and is copied into model general vertical coordinate variable 
 {\bf 
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/PF.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/PF.htm> \end{rawhtml}
 viscAr
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+}. At each time step, the viscous term contribution to the momentum equations
+is calculated in routine
+{\it S/R CALC\_DIFFUSIVITY}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -408,7 +455,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/53.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/53.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -422,16 +469,16 @@
 for this operator are specified later.
 The variable 
 {\bf 
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/SI.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/SI.htm> \end{rawhtml}
 viscAh
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and applied in routines {\it CALC\_MOM\_RHS} and {\it CALC\_GW}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -439,7 +486,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/60.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/60.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -450,7 +497,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/58.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/58.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -465,16 +512,17 @@
 $\frac{\partial v}{\partial x}$=0 along boundaries in $x$.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/UT.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/UT.htm> \end{rawhtml}
 no\_slip\_sides
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and the boundary condition is evaluated in routine
+{\it S/R CALC\_MOM\_RHS}.
 
 
 \fbox{
@@ -483,7 +531,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/60.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/60.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -497,16 +545,16 @@
 e.g. $u=v=0$ at $z=-H$, where $H$ is the local depth of the domain.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/UK.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/UK.htm> \end{rawhtml}
 no\_slip\_bottom
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and is applied in the routine {\it S/R CALC\_MOM\_RHS}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -514,7 +562,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/60.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/60.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -529,23 +577,23 @@
 all boundaries.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/RC.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/RC.htm> \end{rawhtml}
 diffKhT
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and used in routine {\it S/R CALC\_GT}.
 
 \fbox{ \begin{minipage}{5.0in}
 {\it S/R CALC\_GT}({\it calc\_gt.F})
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/57.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/57.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -559,29 +607,29 @@
 operator is $\frac{\partial}{\partial z}$ = 0 on all boundaries.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/ZT.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/ZT.htm> \end{rawhtml}
 diffKzT
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
 }.
 It is copied into model general vertical coordinate variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/PD.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/PD.htm> \end{rawhtml}
 diffKrT
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} which is used in routine {\it S/R CALC\_DIFFUSIVITY}.
 
 \fbox{ \begin{minipage}{5.0in}
 {\it S/R CALC\_DIFFUSIVITY}({\it calc\_diffusivity.F})
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/53.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/53.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -596,16 +644,16 @@
 to $2 \times 10^{-4}\,{\rm degrees}^{-1}$
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/ZV.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/ZV.htm> \end{rawhtml}
 tAlpha 
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+}. The routine {\it S/R FIND\_RHO} makes use of {\bf tAlpha}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -613,7 +661,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/79.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/79.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -625,16 +673,17 @@
 This line selects the linear form of the equation of state.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/WV.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/WV.htm> \end{rawhtml}
 eosType
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+}. The values of {\bf eosType} sets which formula in routine
+{\it FIND\_RHO} is used to calculate density.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -642,7 +691,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/79.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/79.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -655,21 +704,23 @@
 \end{verbatim}
 This line requests that the simulation be performed in a 
 spherical polar coordinate system. It affects the interpretation of
-grid inoput parameters, for exampl {\bf delX} and {\bf delY} and
-causes the grid generation routines to initialise an internal grid based
+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
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/10T.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/10T.htm> \end{rawhtml}
 usingSphericalPolarGrid
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+}. 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 {\it INI\_SPEHRICAL\_POLAR\_GRID}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -677,7 +728,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/97.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/97.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -689,22 +740,22 @@
 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 initialisation of the coriolis force.
+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
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/110.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/110.htm> \end{rawhtml}
 phiMin
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and is used in routine {\it INI\_SPEHRICAL\_POLAR\_GRID}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -712,7 +763,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/97.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/97.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -725,16 +776,16 @@
 in the discrete grid to $1^{\circ}$ in longitude.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/10Z.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/10Z.htm> \end{rawhtml}
 delX
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and is used in routine {\it INI\_SPEHRICAL\_POLAR\_GRID}. 
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -742,7 +793,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/97.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/97.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -755,16 +806,16 @@
 in the discrete grid to $1^{\circ}$ in latitude.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/UB.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/UB.htm> \end{rawhtml}
 delY   
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+} and is used in routine {\it INI\_SPEHRICAL\_POLAR\_GRID}. 
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -772,7 +823,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/97.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/97.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -786,23 +837,23 @@
 is $2\,{\rm km}$.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/10W.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/10W.htm> \end{rawhtml}
 delZ
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
 }.
 It is copied into the internal
 model coordinate variable 
 {\bf 
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/10Y.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/10Y.htm> \end{rawhtml}
 delR
 \begin{rawhtml} </A>\end{rawhtml}
-}. 
+} which is used in routine {\it INI\_VERTICAL\_GRID}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -810,7 +861,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/100.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/100.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -832,16 +883,16 @@
 bathymetry file.
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/179.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/179.htm> \end{rawhtml}
 bathyFile
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+}. The bathymetry file is read in the routine {\it INI\_DEPTHS}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -849,7 +900,7 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/88.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/88.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
@@ -860,7 +911,7 @@
 zonalWindFile='windx.sin_y'
 \end{verbatim}
 This line specifies the name of the file from which the x-direction
-surface wind stress is read. This file is also a two-dimensional
+(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 
@@ -868,16 +919,17 @@
 file.  
 The variable
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/names/13W.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/names/13W.htm> \end{rawhtml}
 zonalWindFile
 \begin{rawhtml} </A>\end{rawhtml}
 }
 is read in the routine
 {\it
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/94.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/94.htm> \end{rawhtml}
 INI\_PARMS
 \begin{rawhtml} </A>\end{rawhtml}
-}.
+}.  The wind-stress file is read in the routine 
+{\it EXTERNAL\_FIELDS\_LOAD}.
 
 \fbox{
 \begin{minipage}{5.0in}
@@ -885,16 +937,14 @@
 \end{minipage}
 }
 {\bf
-\begin{rawhtml} <A href=../../../code_reference/vdb/code/75.htm> \end{rawhtml}
+\begin{rawhtml} <A href=../code_reference/vdb/code/75.htm> \end{rawhtml}
 goto code
 \begin{rawhtml} </A>\end{rawhtml}
 }
 
 \end{itemize}
 
-\noindent other lines in the file {\it input/data} are standard values
-that are described in the MITgcm Getting Started and MITgcm Parameters
-notes.
+\noindent other lines in the file {\it input/data} are standard values.
 
 \begin{rawhtml}<PRE>\end{rawhtml}
 \begin{small}
@@ -903,37 +953,45 @@
 \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}$.
-Although $\tau_{x}$ is only a function of $y$n in this experiment
+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. The included matlab program {\it input/gendata.m} gives a complete
+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
-$0m$ or $-2000\,{\rm m}$, corresponding respectively to a wall or to deep
+$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
 
@@ -960,17 +1018,20 @@
 \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}
@@ -983,22 +1044,26 @@
 \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 the Getting Started
-Guide \cite{MITgcm_Getting_Started}.
+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/exp1/ }
+{\it verification/exp2/ }
 
 \subsubsection{Running the Experiment}
+\label{www:tutorials}
 
  To run the experiment
 
@@ -1015,9 +1080,9 @@
 % pwd
 \end{verbatim}
 
- You shold see a response on the screen ending in
+ You should see a response on the screen ending in
 
-{\it verification/exp1/input }
+{\it verification/exp2/input }
 
 
 \item Run the genmake script to create the experiment {\it Makefile}
