--- manual/s_examples/baroclinic_gyre/fourlayer.tex 2001/10/24 15:21:27 1.3
+++ manual/s_examples/baroclinic_gyre/fourlayer.tex 2001/10/25 12:06:56 1.8
@@ -1,4 +1,4 @@
-% $Header: /home/ubuntu/mnt/e9_copy/manual/s_examples/baroclinic_gyre/fourlayer.tex,v 1.3 2001/10/24 15:21:27 cnh Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/manual/s_examples/baroclinic_gyre/fourlayer.tex,v 1.8 2001/10/25 12:06:56 cnh Exp $
% $Name: $
\section{Example: Four layer Baroclinic Ocean Gyre In Spherical Coordinates}
@@ -37,8 +37,8 @@
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$
@@ -54,7 +54,7 @@
\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
@@ -64,11 +64,11 @@
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
+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
+the local orthogonal model coordinate $(x,y,z)$ is synonomous
+with the coordinates $(\lambda,\varphi,r)$ shown in figure
\ref{fig:spherical-polar-coord}
\\
@@ -95,7 +95,7 @@
\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,7 +110,7 @@
\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.
}
@@ -118,68 +118,97 @@
\end{figure}
\subsection{Equations solved}
-
-The implicit free surface form of the
-pressure equation described in Marshall et. al \cite{Marshall97a} is
-employed.
+For this problem
+the implicit free surface, {\bf HPE} (see section \ref{sec:hydrostatic_and_quasi-hydrostatic_forms}) form of the
+equations described in Marshall et. al \cite{Marshall97a} 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. The wind-stress momentum input is added to the momentum equation
-for the ``zonal flow'', $u$. Other terms in the model
+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 experiement configuration (see section
-\ref{SEC:code_config} ). This yields an active set of equations in
+\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}
\frac{Du}{Dt} - fv +
- \frac{1}{\rho}\frac{\partial p^{'}}{\partial \lambda} -
+ \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^{'}}{\partial \varphi} -
+ \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 \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^{'} & = &
-g\rho_{0} \eta + \int^{0}_{-z}\rho^{'} dz
+p^{\prime} & = &
+g\rho_{0} \eta + \int^{0}_{-z}\rho^{\prime} dz
\\
-\rho^{'} & = &- \alpha_{\theta}\rho_{0}\theta^{'}
+\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 suffices ${s},{i}$ indicate surface and interior of the domain.
-The forcing $\cal F$ is only applied at the surface.
-The pressure field $p^{'}$ is separated into a barotropic part
+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^{'}$, 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}
- 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 a four layers with constant depth,
+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
$\lambda$, $\varphi$, $\Delta \lambda$ and $\Delta \varphi$ in
@@ -215,15 +244,26 @@
As described in \ref{sec: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{sec: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{sec: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{sec:finding_the_pressure_field}.
\subsubsection{Numerical Stability Criteria}
-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}
@@ -231,13 +271,15 @@
\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}
@@ -245,7 +287,7 @@
\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
@@ -263,7 +305,6 @@
\\
\noindent The numerical stability for inertial oscillations
-\cite{Adcroft_thesis}
\begin{eqnarray}
\label{EQ:inertial_stability}
@@ -274,35 +315,35 @@
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}
+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
+propogating at $2~{\rm m}~{\rm s}^{-1}$
\begin{eqnarray}
\label{EQ: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{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}
@@ -397,7 +438,9 @@
\begin{rawhtml} \end{rawhtml}
viscAr
\begin{rawhtml} \end{rawhtml}
-}.
+}. At each time step, the viscous term contribution to the momentum eqautions
+is calculated in routine
+{\it S/R CALC\_DIFFUSIVITY}.
\fbox{
\begin{minipage}{5.0in}
@@ -428,7 +471,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+} and applied in routines {\it CALC\_MOM\_RHS} and {\it CALC\_GW}.
\fbox{
\begin{minipage}{5.0in}
@@ -471,7 +514,8 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+} and the boundary condition is evaluated in routine
+{\it S/R CALC\_MOM\_RHS}.
\fbox{
@@ -503,7 +547,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+} and is applied in the routine {\it S/R CALC\_MOM\_RHS}.
\fbox{
\begin{minipage}{5.0in}
@@ -535,7 +579,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \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})
@@ -571,7 +615,7 @@
\begin{rawhtml} \end{rawhtml}
diffKrT
\begin{rawhtml} \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})
@@ -602,7 +646,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+}. The routine {\it S/R FIND\_RHO} makes use of {\bf tAlpha}.
\fbox{
\begin{minipage}{5.0in}
@@ -631,7 +675,8 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \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}
@@ -666,7 +711,9 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \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}
@@ -701,7 +748,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+} and is used in routine {\it INI\_SPEHRICAL\_POLAR\_GRID}.
\fbox{
\begin{minipage}{5.0in}
@@ -731,7 +778,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+} and is used in routine {\it INI\_SPEHRICAL\_POLAR\_GRID}.
\fbox{
\begin{minipage}{5.0in}
@@ -761,7 +808,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+} and is used in routine {\it INI\_SPEHRICAL\_POLAR\_GRID}.
\fbox{
\begin{minipage}{5.0in}
@@ -799,7 +846,7 @@
\begin{rawhtml} \end{rawhtml}
delR
\begin{rawhtml} \end{rawhtml}
-}.
+} which is used in routine {\it INI\_VERTICAL\_GRID}.
\fbox{
\begin{minipage}{5.0in}
@@ -838,7 +885,7 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+}. The bathymetry file is read in the routine {\it INI\_DEPTHS}.
\fbox{
\begin{minipage}{5.0in}
@@ -857,7 +904,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
@@ -874,7 +921,8 @@
\begin{rawhtml} \end{rawhtml}
INI\_PARMS
\begin{rawhtml} \end{rawhtml}
-}.
+}. The wind-stress file is read in the routine
+{\it EXTERNAL\_FIELDS\_LOAD}.
\fbox{
\begin{minipage}{5.0in}
@@ -889,9 +937,7 @@
\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}\end{rawhtml}
\begin{small}
@@ -912,11 +958,14 @@
\subsubsection{File {\it input/windx.sin\_y}}
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 latituted, $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}}
@@ -924,7 +973,7 @@
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
@@ -985,15 +1034,15 @@
\subsubsection{Code Download}
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}
This example experiments is located under the release sub-directory
\vspace{5mm}
-{\it verification/exp1/ }
+{\it verification/exp2/ }
\subsubsection{Running the Experiment}
@@ -1014,7 +1063,7 @@
You shold 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}