s_examples/tracer_adjsens/co2sens.tex

% $Header: /u/u0/gcmpack/manual/part3/case_studies/carbon_outgassing_sensitivity/co2sens.tex,v 1.5 2001/11/13 19:01:42 adcroft Exp $
% $Name:  $

\section{Centennial Time Scale Tracer Injection}
\label{sect:eg-simple-tracer}

\bodytext{bgcolor="#FFFFFFFF"}

%\begin{center} 
%{\Large \bf Using MITgcm to Look at Centennial Time Scale
%Sensitivities}
%
%\vspace*{4mm}
%
%\vspace*{3mm}
%{\large May 2001}
%\end{center}

\subsection{Introduction}

This document describes the fourth example MITgcm experiment.
This example illustrates the use of
the MITgcm to perform sensitivity analysis in a
large scale ocean circulation simulation.

\subsection{Overview}

This example experiment demonstrates using the MITgcm to simulate
the planetary ocean circulation. The simulation is configured
with realistic geography and bathymetry on a
$4^{\circ} \times 4^{\circ}$ spherical polar grid.
Twenty vertical layers are used in the vertical, ranging in thickness
from $50\,{\rm m}$ at the surface to $815\,{\rm m}$ at depth,
giving a maximum model depth of $6\,{\rm km}$.
At this resolution, the configuration
can be integrated forward for thousands of years on a single 
processor desktop computer.
\\

The model is forced with climatological wind stress data and surface
flux data from Da Silva \cite{DaSilva94}. Climatological data
from Levitus \cite{Levitus94} is used to initialize the model hydrography.
Levitus data is also used throughout the calculation
to derive air-sea fluxes of heat at the ocean surface.
These fluxes are combined with climatological estimates of
surface heat flux and fresh water, resulting in a mixed boundary
condition of the style described in Haney \cite{Haney}.
Altogether, this yields the following forcing applied
in the model surface layer.

\begin{eqnarray}
\label{EQ:eg-simple-tracer-global_forcing}
\label{EQ:eg-simple-tracer-global_forcing_fu}
{\cal F}_{u} & = & \frac{\tau_{x}}{\rho_{0} \Delta z_{s}}
\\
\label{EQ:eg-simple-tracer-global_forcing_fv}
{\cal F}_{v} & = & \frac{\tau_{y}}{\rho_{0} \Delta z_{s}}
\\
\label{EQ:eg-simple-tracer-global_forcing_ft}
{\cal F}_{\theta} & = & - \lambda_{\theta} ( \theta - \theta^{\ast} ) 
 - \frac{1}{C_{p} \rho_{0} \Delta z_{s}}{\cal Q}
\\
\label{EQ:eg-simple-tracer-global_forcing_fs}
{\cal F}_{s} & = & - \lambda_{s} ( S - S^{\ast} ) 
 + \frac{S_{0}}{\Delta z_{s}}({\cal E} - {\cal P} - {\cal R})
\end{eqnarray}

\noindent where ${\cal F}_{u}$, ${\cal F}_{v}$, ${\cal F}_{\theta}$,
${\cal F}_{s}$ are the forcing terms in the zonal and meridional
momentum and in the potential temperature and salinity 
equations respectively.
The term $\Delta z_{s}$ represents the top ocean layer thickness.
It is used in conjunction with the reference density, $\rho_{0}$
(here set to $999.8\,{\rm kg\,m^{-3}}$), the
reference salinity, $S_{0}$ (here set to 35ppt),
and a specific heat capacity $C_{p}$ to convert
wind-stress fluxes given in ${\rm N}\,m^{-2}$, 
\\


The configuration is illustrated in figure \ref{simulation_config}.


\subsection{Discrete Numerical Configuration}


 The model is configured in hydrostatic form.  The domain is discretised with 
a uniform grid spacing in latitude and longitude of
 $\Delta x=\Delta y=4^{\circ}$, so 
that there are ninety grid cells in the $x$ and forty in the 
$y$ direction (Arctic polar regions are not
included in this experiment). Vertically the 
model is configured with twenty layers with the following thicknesses
$\Delta z_{1} = 50\,{\rm m},\,
 \Delta z_{2} = 50\,{\rm m},\,
 \Delta z_{3} = 55\,{\rm m},\,
 \Delta z_{4} = 60\,{\rm m},\,
 \Delta z_{5} = 65\,{\rm m},\,
$
$
 \Delta z_{6}~=~70\,{\rm m},\,
 \Delta z_{7}~=~80\,{\rm m},\,
 \Delta z_{8}~=95\,{\rm m},\,
 \Delta z_{9}=120\,{\rm m},\,
 \Delta z_{10}=155\,{\rm m},\,
$
$
 \Delta z_{11}=200\,{\rm m},\,
 \Delta z_{12}=260\,{\rm m},\,
 \Delta z_{13}=320\,{\rm m},\,
 \Delta z_{14}=400\,{\rm m},\,
 \Delta z_{15}=480\,{\rm m},\,
$
$
 \Delta z_{16}=570\,{\rm m},\,
 \Delta z_{17}=655\,{\rm m},\,
 \Delta z_{18}=725\,{\rm m},\,
 \Delta z_{19}=775\,{\rm m},\,
 \Delta z_{20}=815\,{\rm m}
$ (here the numeric subscript indicates the model level index number, ${\tt k}$).
The implicit free surface form of the pressure equation described in Marshall et. al 
\cite{marshall:97a} is employed. A Laplacian operator, $\nabla^2$, provides viscous
dissipation. Thermal and haline diffusion is also represented by a Laplacian operator.
\\

Wind-stress momentum inputs are added to the momentum equations for both
the zonal flow, $u$ and the meridional flow $v$, according to equations 
(\ref{EQ:eg-simple-tracer-global_forcing_fu}) and (\ref{EQ:eg-simple-tracer-global_forcing_fv}).
Thermodynamic forcing inputs are added to the equations for
potential temperature, $\theta$, and salinity, $S$, according to equations 
(\ref{EQ:eg-simple-tracer-global_forcing_ft}) and (\ref{EQ:eg-simple-tracer-global_forcing_fs}).
This produces a set of equations solved in this configuration as follows:
% {\fracktur}


\begin{eqnarray}
\label{EQ:eg-simple-tracer-model_equations}
\frac{Du}{Dt} - fv + 
  \frac{1}{\rho}\frac{\partial p^{'}}{\partial x} - 
  A_{h}\nabla_{h}^2u - A_{z}\frac{\partial^{2}u}{\partial z^{2}} 
& = &
{\cal F}_{u}
\\
\frac{Dv}{Dt} + fu + 
  \frac{1}{\rho}\frac{\partial p^{'}}{\partial y} - 
  A_{h}\nabla_{h}^2v - A_{z}\frac{\partial^{2}v}{\partial z^{2}} 
& = &
{\cal F}_{v}
\\
\frac{\partial \eta}{\partial t} + \nabla_{h}\cdot \vec{u}
&=&
0
\\
\frac{D\theta}{Dt} -
 K_{h}\nabla_{h}^2\theta  - \Gamma(K_{z})\frac{\partial^{2}\theta}{\partial z^{2}} 
& = &
{\cal F}_{\theta}
\\
\frac{D s}{Dt} -
 K_{h}\nabla_{h}^2 s  - \Gamma(K_{z})\frac{\partial^{2} s}{\partial z^{2}} 
& = &
{\cal F}_{s}
\\
g\rho_{0} \eta + \int^{0}_{-z}\rho^{'} dz & = & p^{'}
\\
\end{eqnarray}

\noindent where $u$ and $v$ are the $x$ and $y$ components of the
flow vector $\vec{u}$. The suffices ${s},{i}$ indicate surface and
interior model levels respectively. As described in
MITgcm Numerical Solution Procedure \ref{chap:discretization}, the time 
evolution of potential temperature, $\theta$, equation is solved prognostically.
The total pressure, $p$, is diagnosed by summing pressure due to surface 
elevation $\eta$ and the hydrostatic pressure.
\\

\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:95},

\begin{eqnarray}
\label{EQ:eg-simple-tracer-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 
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:95}

\begin{eqnarray}
\label{EQ:eg-simple-tracer-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. 
\\

\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-simple-tracer-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
\cite{adcroft:95} 

\begin{eqnarray}
\label{EQ:eg-simple-tracer-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 \cite{adcroft:95} for a extreme maximum 
horizontal flow
speed of $ | \vec{u} | = 2 ms^{-1}$

\begin{eqnarray}
\label{EQ:eg-simple-tracer-cfl_stability}
S_{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:95}

\begin{eqnarray}
\label{EQ:eg-simple-tracer-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
stability limit of 0.25.
  
\subsection{Code Configuration}
\label{SEC:code_config}

The model configuration for this experiment resides under the 
directory {\it verification/exp1/}.  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 customizations and parameter settings for this 
experiments. Below we describe the customizations
to these files associated with this experiment.

\subsubsection{File {\it input/data}}

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}$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$.

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


\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. This variable is copied into
model general vertical coordinate variable {\bf viscAr}.

\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.

\item Lines 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$.

\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.

\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.

\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.

\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}$

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

\item Line 18,
\begin{verbatim}
eosType='LINEAR'
\end{verbatim}
This line selects the linear form of the equation of state.

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


\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.

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

\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.

\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.

\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.

\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 {\bf delZ} is copied into the internal
model coordinate variable {\bf delR}

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

\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.


\item Line 50,
\begin{verbatim}
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
($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.  

\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.

\begin{small}
% \input{part3/case_studies/carbon_outgassing_sensitivity/input/data}
\end{small}

\subsubsection{File {\it input/data.pkg}}

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

\subsubsection{File {\it input/eedata}}

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

\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
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
code for creating the {\it input/windx.sin\_y} file.

\subsubsection{File {\it input/topog.box}}


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
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}}

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{code/SIZE.h}
\end{small}

\subsubsection{File {\it code/CPP\_OPTIONS.h}}

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


\subsubsection{File {\it code/CPP\_EEOPTIONS.h}}

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

\subsubsection{Other Files }

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}.
\item {\it model/src/ini\_spherical\_polar\_grid.F}
\item {\it model/src/ini\_parms.F},
\item {\it input/windx.sin\_y},
\end{itemize}
contain the code customizations and parameter settings for this 
experiments. Below we describe the customizations
to these files associated with this experiment.
1	% $Header: /u/u0/gcmpack/manual/part3/case_studies/carbon_outgassing_sensitivity/co2sens.tex,v 1.5 2001/11/13 19:01:42 adcroft Exp $
2	% $Name: $
3
4	\section{Centennial Time Scale Tracer Injection}
5	\label{sect:eg-simple-tracer}
6
7	\bodytext{bgcolor="#FFFFFFFF"}
8
9	%\begin{center}
10	%{\Large \bf Using MITgcm to Look at Centennial Time Scale
11	%Sensitivities}
12	%
13	%\vspace*{4mm}
14	%
15	%\vspace*{3mm}
16	%{\large May 2001}
17	%\end{center}
18
19	\subsection{Introduction}
20
21	This document describes the fourth example MITgcm experiment.
22	This example illustrates the use of
23	the MITgcm to perform sensitivity analysis in a
24	large scale ocean circulation simulation.
25
26	\subsection{Overview}
27
28	This example experiment demonstrates using the MITgcm to simulate
29	the planetary ocean circulation. The simulation is configured
30	with realistic geography and bathymetry on a
31	$4^{\circ} \times 4^{\circ}$ spherical polar grid.
32	Twenty vertical layers are used in the vertical, ranging in thickness
33	from $50\,{\rm m}$ at the surface to $815\,{\rm m}$ at depth,
34	giving a maximum model depth of $6\,{\rm km}$.
35	At this resolution, the configuration
36	can be integrated forward for thousands of years on a single
37	processor desktop computer.
38	\\
39
40	The model is forced with climatological wind stress data and surface
41	flux data from Da Silva \cite{DaSilva94}. Climatological data
42	from Levitus \cite{Levitus94} is used to initialize the model hydrography.
43	Levitus data is also used throughout the calculation
44	to derive air-sea fluxes of heat at the ocean surface.
45	These fluxes are combined with climatological estimates of
46	surface heat flux and fresh water, resulting in a mixed boundary
47	condition of the style described in Haney \cite{Haney}.
48	Altogether, this yields the following forcing applied
49	in the model surface layer.
50
51	\begin{eqnarray}
52	\label{EQ:eg-simple-tracer-global_forcing}
53	\label{EQ:eg-simple-tracer-global_forcing_fu}
54	{\cal F}_{u} & = & \frac{\tau_{x}}{\rho_{0} \Delta z_{s}}
55	\\
56	\label{EQ:eg-simple-tracer-global_forcing_fv}
57	{\cal F}_{v} & = & \frac{\tau_{y}}{\rho_{0} \Delta z_{s}}
58	\\
59	\label{EQ:eg-simple-tracer-global_forcing_ft}
60	{\cal F}_{\theta} & = & - \lambda_{\theta} ( \theta - \theta^{\ast} )
61	- \frac{1}{C_{p} \rho_{0} \Delta z_{s}}{\cal Q}
62	\\
63	\label{EQ:eg-simple-tracer-global_forcing_fs}
64	{\cal F}_{s} & = & - \lambda_{s} ( S - S^{\ast} )
65	+ \frac{S_{0}}{\Delta z_{s}}({\cal E} - {\cal P} - {\cal R})
66	\end{eqnarray}
67
68	\noindent where ${\cal F}_{u}$, ${\cal F}_{v}$, ${\cal F}_{\theta}$,
69	${\cal F}_{s}$ are the forcing terms in the zonal and meridional
70	momentum and in the potential temperature and salinity
71	equations respectively.
72	The term $\Delta z_{s}$ represents the top ocean layer thickness.
73	It is used in conjunction with the reference density, $\rho_{0}$
74	(here set to $999.8\,{\rm kg\,m^{-3}}$), the
75	reference salinity, $S_{0}$ (here set to 35ppt),
76	and a specific heat capacity $C_{p}$ to convert
77	wind-stress fluxes given in ${\rm N}\,m^{-2}$,
78	\\
79
80
81	The configuration is illustrated in figure \ref{simulation_config}.
82
83
84	\subsection{Discrete Numerical Configuration}
85
86
87	The model is configured in hydrostatic form. The domain is discretised with
88	a uniform grid spacing in latitude and longitude of
89	$\Delta x=\Delta y=4^{\circ}$, so
90	that there are ninety grid cells in the $x$ and forty in the
91	$y$ direction (Arctic polar regions are not
92	included in this experiment). Vertically the
93	model is configured with twenty layers with the following thicknesses
94	$\Delta z_{1} = 50\,{\rm m},\,
95	\Delta z_{2} = 50\,{\rm m},\,
96	\Delta z_{3} = 55\,{\rm m},\,
97	\Delta z_{4} = 60\,{\rm m},\,
98	\Delta z_{5} = 65\,{\rm m},\,
99	$
100	$
101	\Delta z_{6}~=~70\,{\rm m},\,
102	\Delta z_{7}~=~80\,{\rm m},\,
103	\Delta z_{8}~=95\,{\rm m},\,
104	\Delta z_{9}=120\,{\rm m},\,
105	\Delta z_{10}=155\,{\rm m},\,
106	$
107	$
108	\Delta z_{11}=200\,{\rm m},\,
109	\Delta z_{12}=260\,{\rm m},\,
110	\Delta z_{13}=320\,{\rm m},\,
111	\Delta z_{14}=400\,{\rm m},\,
112	\Delta z_{15}=480\,{\rm m},\,
113	$
114	$
115	\Delta z_{16}=570\,{\rm m},\,
116	\Delta z_{17}=655\,{\rm m},\,
117	\Delta z_{18}=725\,{\rm m},\,
118	\Delta z_{19}=775\,{\rm m},\,
119	\Delta z_{20}=815\,{\rm m}
120	$ (here the numeric subscript indicates the model level index number, ${\tt k}$).
121	The implicit free surface form of the pressure equation described in Marshall et. al
122	\cite{marshall:97a} is employed. A Laplacian operator, $\nabla^2$, provides viscous
123	dissipation. Thermal and haline diffusion is also represented by a Laplacian operator.
124	\\
125
126	Wind-stress momentum inputs are added to the momentum equations for both
127	the zonal flow, $u$ and the meridional flow $v$, according to equations
128	(\ref{EQ:eg-simple-tracer-global_forcing_fu}) and (\ref{EQ:eg-simple-tracer-global_forcing_fv}).
129	Thermodynamic forcing inputs are added to the equations for
130	potential temperature, $\theta$, and salinity, $S$, according to equations
131	(\ref{EQ:eg-simple-tracer-global_forcing_ft}) and (\ref{EQ:eg-simple-tracer-global_forcing_fs}).
132	This produces a set of equations solved in this configuration as follows:
133	% {\fracktur}
134
135
136	\begin{eqnarray}
137	\label{EQ:eg-simple-tracer-model_equations}
138	\frac{Du}{Dt} - fv +
139	\frac{1}{\rho}\frac{\partial p^{'}}{\partial x} -
140	A_{h}\nabla_{h}^2u - A_{z}\frac{\partial^{2}u}{\partial z^{2}}
141	& = &
142	{\cal F}_{u}
143	\\
144	\frac{Dv}{Dt} + fu +
145	\frac{1}{\rho}\frac{\partial p^{'}}{\partial y} -
146	A_{h}\nabla_{h}^2v - A_{z}\frac{\partial^{2}v}{\partial z^{2}}
147	& = &
148	{\cal F}_{v}
149	\\
150	\frac{\partial \eta}{\partial t} + \nabla_{h}\cdot \vec{u}
151	&=&
152	0
153	\\
154	\frac{D\theta}{Dt} -
155	K_{h}\nabla_{h}^2\theta - \Gamma(K_{z})\frac{\partial^{2}\theta}{\partial z^{2}}
156	& = &
157	{\cal F}_{\theta}
158	\\
159	\frac{D s}{Dt} -
160	K_{h}\nabla_{h}^2 s - \Gamma(K_{z})\frac{\partial^{2} s}{\partial z^{2}}
161	& = &
162	{\cal F}_{s}
163	\\
164	g\rho_{0} \eta + \int^{0}_{-z}\rho^{'} dz & = & p^{'}
165	\\
166	\end{eqnarray}
167
168	\noindent where $u$ and $v$ are the $x$ and $y$ components of the
169	flow vector $\vec{u}$. The suffices ${s},{i}$ indicate surface and
170	interior model levels respectively. As described in
171	MITgcm Numerical Solution Procedure \ref{chap:discretization}, the time
172	evolution of potential temperature, $\theta$, equation is solved prognostically.
173	The total pressure, $p$, is diagnosed by summing pressure due to surface
174	elevation $\eta$ and the hydrostatic pressure.
175	\\
176
177	\subsubsection{Numerical Stability Criteria}
178
179	The Laplacian dissipation coefficient, $A_{h}$, is set to $400 m s^{-1}$.
180	This value is chosen to yield a Munk layer width \cite{adcroft:95},
181
182	\begin{eqnarray}
183	\label{EQ:eg-simple-tracer-munk_layer}
184	M_{w} = \pi ( \frac { A_{h} }{ \beta } )^{\frac{1}{3}}
185	\end{eqnarray}
186
187	\noindent of $\approx 100$km. This is greater than the model
188	resolution in mid-latitudes $\Delta x$, ensuring that the frictional
189	boundary layer is well resolved.
190	\\
191
192	\noindent The model is stepped forward with a
193	time step $\delta t=1200$secs. With this time step the stability
194	parameter to the horizontal Laplacian friction \cite{adcroft:95}
195
196	\begin{eqnarray}
197	\label{EQ:eg-simple-tracer-laplacian_stability}
198	S_{l} = 4 \frac{A_{h} \delta t}{{\Delta x}^2}
199	\end{eqnarray}
200
201	\noindent evaluates to 0.012, which is well below the 0.3 upper limit
202	for stability.
203	\\
204
205	\noindent The vertical dissipation coefficient, $A_{z}$, is set to
206	$1\times10^{-2} {\rm m}^2{\rm s}^{-1}$. The associated stability limit
207
208	\begin{eqnarray}
209	\label{EQ:eg-simple-tracer-laplacian_stability_z}
210	S_{l} = 4 \frac{A_{z} \delta t}{{\Delta z}^2}
211	\end{eqnarray}
212
213	\noindent evaluates to $4.8 \times 10^{-5}$ which is again well below
214	the upper limit.
215	The values of $A_{h}$ and $A_{z}$ are also used for the horizontal ($K_{h}$)
216	and vertical ($K_{z}$) diffusion coefficients for temperature respectively.
217	\\
218
219	\noindent The numerical stability for inertial oscillations
220	\cite{adcroft:95}
221
222	\begin{eqnarray}
223	\label{EQ:eg-simple-tracer-inertial_stability}
224	S_{i} = f^{2} {\delta t}^2
225	\end{eqnarray}
226
227	\noindent evaluates to $0.0144$, which is well below the $0.5$ upper
228	limit for stability.
229	\\
230
231	\noindent The advective CFL \cite{adcroft:95} for a extreme maximum
232	horizontal flow
233	speed of $ \| \vec{u} \| = 2 ms^{-1}$
234
235	\begin{eqnarray}
236	\label{EQ:eg-simple-tracer-cfl_stability}
237	S_{a} = \frac{\| \vec{u} \| \delta t}{ \Delta x}
238	\end{eqnarray}
239
240	\noindent evaluates to $5 \times 10^{-2}$. This is well below the stability
241	limit of 0.5.
242	\\
243
244	\noindent The stability parameter for internal gravity waves
245	\cite{adcroft:95}
246
247	\begin{eqnarray}
248	\label{EQ:eg-simple-tracer-igw_stability}
249	S_{c} = \frac{c_{g} \delta t}{ \Delta x}
250	\end{eqnarray}
251
252	\noindent evaluates to $5 \times 10^{-2}$. This is well below the linear
253	stability limit of 0.25.
254
255	\subsection{Code Configuration}
256	\label{SEC:code_config}
257
258	The model configuration for this experiment resides under the
259	directory {\it verification/exp1/}. The experiment files
260	\begin{itemize}
261	\item {\it input/data}
262	\item {\it input/data.pkg}
263	\item {\it input/eedata},
264	\item {\it input/windx.sin\_y},
265	\item {\it input/topog.box},
266	\item {\it code/CPP\_EEOPTIONS.h}
267	\item {\it code/CPP\_OPTIONS.h},
268	\item {\it code/SIZE.h}.
269	\end{itemize}
270	contain the code customizations and parameter settings for this
271	experiments. Below we describe the customizations
272	to these files associated with this experiment.
273
274	\subsubsection{File {\it input/data}}
275
276	This file, reproduced completely below, specifies the main parameters
277	for the experiment. The parameters that are significant for this configuration
278	are
279
280	\begin{itemize}
281
282	\item Line 4,
283	\begin{verbatim} tRef=20.,10.,8.,6., \end{verbatim}
284	this line sets
285	the initial and reference values of potential temperature at each model
286	level in units of $^{\circ}$C.
287	The entries are ordered from surface to depth. For each
288	depth level the initial and reference profiles will be uniform in
289	$x$ and $y$.
290
291	\fbox{
292	\begin{minipage}{5.0in}
293	{\it S/R INI\_THETA}({\it ini\_theta.F})
294	\end{minipage}
295	}
296
297
298	\item Line 6,
299	\begin{verbatim} viscAz=1.E-2, \end{verbatim}
300	this line sets the vertical Laplacian dissipation coefficient to
301	$1 \times 10^{-2} {\rm m^{2}s^{-1}}$. Boundary conditions
302	for this operator are specified later. This variable is copied into
303	model general vertical coordinate variable {\bf viscAr}.
304
305	\fbox{
306	\begin{minipage}{5.0in}
307	{\it S/R CALC\_DIFFUSIVITY}({\it calc\_diffusivity.F})
308	\end{minipage}
309	}
310
311	\item Line 7,
312	\begin{verbatim}
313	viscAh=4.E2,
314	\end{verbatim}
315	this line sets the horizontal Laplacian frictional dissipation coefficient to
316	$1 \times 10^{-2} {\rm m^{2}s^{-1}}$. Boundary conditions
317	for this operator are specified later.
318
319	\item Lines 8,
320	\begin{verbatim}
321	no_slip_sides=.FALSE.
322	\end{verbatim}
323	this line selects a free-slip lateral boundary condition for
324	the horizontal Laplacian friction operator
325	e.g. $\frac{\partial u}{\partial y}$=0 along boundaries in $y$ and
326	$\frac{\partial v}{\partial x}$=0 along boundaries in $x$.
327
328	\item Lines 9,
329	\begin{verbatim}
330	no_slip_bottom=.TRUE.
331	\end{verbatim}
332	this line selects a no-slip boundary condition for bottom
333	boundary condition in the vertical Laplacian friction operator
334	e.g. $u=v=0$ at $z=-H$, where $H$ is the local depth of the domain.
335
336	\item Line 10,
337	\begin{verbatim}
338	diffKhT=4.E2,
339	\end{verbatim}
340	this line sets the horizontal diffusion coefficient for temperature
341	to $400\,{\rm m^{2}s^{-1}}$. The boundary condition on this
342	operator is $\frac{\partial}{\partial x}=\frac{\partial}{\partial y}=0$ at
343	all boundaries.
344
345	\item Line 11,
346	\begin{verbatim}
347	diffKzT=1.E-2,
348	\end{verbatim}
349	this line sets the vertical diffusion coefficient for temperature
350	to $10^{-2}\,{\rm m^{2}s^{-1}}$. The boundary condition on this
351	operator is $\frac{\partial}{\partial z}$ = 0 on all boundaries.
352
353	\item Line 13,
354	\begin{verbatim}
355	tAlpha=2.E-4,
356	\end{verbatim}
357	This line sets the thermal expansion coefficient for the fluid
358	to $2 \times 10^{-4}\,{\rm degrees}^{-1}$
359
360	\fbox{
361	\begin{minipage}{5.0in}
362	{\it S/R FIND\_RHO}({\it find\_rho.F})
363	\end{minipage}
364	}
365
366	\item Line 18,
367	\begin{verbatim}
368	eosType='LINEAR'
369	\end{verbatim}
370	This line selects the linear form of the equation of state.
371
372	\fbox{
373	\begin{minipage}{5.0in}
374	{\it S/R FIND\_RHO}({\it find\_rho.F})
375	\end{minipage}
376	}
377
378
379
380	\item Line 40,
381	\begin{verbatim}
382	usingSphericalPolarGrid=.TRUE.,
383	\end{verbatim}
384	This line requests that the simulation be performed in a
385	spherical polar coordinate system. It affects the interpretation of
386	grid input parameters, for example {\bf delX} and {\bf delY} and
387	causes the grid generation routines to initialize an internal grid based
388	on spherical polar geometry.
389
390	\fbox{
391	\begin{minipage}{5.0in}
392	{\it S/R INI\_SPEHRICAL\_POLAR\_GRID}({\it ini\_spherical\_polar\_grid.F})
393	\end{minipage}
394	}
395
396	\item Line 41,
397	\begin{verbatim}
398	phiMin=0.,
399	\end{verbatim}
400	This line sets the southern boundary of the modeled
401	domain to $0^{\circ}$ latitude. This value affects both the
402	generation of the locally orthogonal grid that the model
403	uses internally and affects the initialization of the coriolis force.
404	Note - it is not required to set
405	a longitude boundary, since the absolute longitude does
406	not alter the kernel equation discretisation.
407
408	\item Line 42,
409	\begin{verbatim}
410	delX=60*1.,
411	\end{verbatim}
412	This line sets the horizontal grid spacing between each y-coordinate line
413	in the discrete grid to $1^{\circ}$ in longitude.
414
415	\item Line 43,
416	\begin{verbatim}
417	delY=60*1.,
418	\end{verbatim}
419	This line sets the horizontal grid spacing between each y-coordinate line
420	in the discrete grid to $1^{\circ}$ in latitude.
421
422	\item Line 44,
423	\begin{verbatim}
424	delZ=500.,500.,500.,500.,
425	\end{verbatim}
426	This line sets the vertical grid spacing between each z-coordinate line
427	in the discrete grid to $500\,{\rm m}$, so that the total model depth
428	is $2\,{\rm km}$. The variable {\bf delZ} is copied into the internal
429	model coordinate variable {\bf delR}
430
431	\fbox{
432	\begin{minipage}{5.0in}
433	{\it S/R INI\_VERTICAL\_GRID}({\it ini\_vertical\_grid.F})
434	\end{minipage}
435	}
436
437	\item Line 47,
438	\begin{verbatim}
439	bathyFile='topog.box'
440	\end{verbatim}
441	This line specifies the name of the file from which the domain
442	bathymetry is read. This file is a two-dimensional ($x,y$) map of
443	depths. This file is assumed to contain 64-bit binary numbers
444	giving the depth of the model at each grid cell, ordered with the x
445	coordinate varying fastest. The points are ordered from low coordinate
446	to high coordinate for both axes. The units and orientation of the
447	depths in this file are the same as used in the MITgcm code. In this
448	experiment, a depth of $0m$ indicates a solid wall and a depth
449	of $-2000m$ indicates open ocean. The matlab program
450	{\it input/gendata.m} shows an example of how to generate a
451	bathymetry file.
452
453
454	\item Line 50,
455	\begin{verbatim}
456	zonalWindFile='windx.sin_y'
457	\end{verbatim}
458	This line specifies the name of the file from which the x-direction
459	surface wind stress is read. This file is also a two-dimensional
460	($x,y$) map and is enumerated and formatted in the same manner as the
461	bathymetry file. The matlab program {\it input/gendata.m} includes example
462	code to generate a valid
463	{\bf zonalWindFile}
464	file.
465
466	\end{itemize}
467
468	\noindent other lines in the file {\it input/data} are standard values
469	that are described in the MITgcm Getting Started and MITgcm Parameters
470	notes.
471
472	\begin{small}
473	% \input{part3/case_studies/carbon_outgassing_sensitivity/input/data}
474	\end{small}
475
476	\subsubsection{File {\it input/data.pkg}}
477
478	This file uses standard default values and does not contain
479	customizations for this experiment.
480
481	\subsubsection{File {\it input/eedata}}
482
483	This file uses standard default values and does not contain
484	customizations for this experiment.
485
486	\subsubsection{File {\it input/windx.sin\_y}}
487
488	The {\it input/windx.sin\_y} file specifies a two-dimensional ($x,y$)
489	map of wind stress ,$\tau_{x}$, values. The units used are $Nm^{-2}$.
490	Although $\tau_{x}$ is only a function of $y$n in this experiment
491	this file must still define a complete two-dimensional map in order
492	to be compatible with the standard code for loading forcing fields
493	in MITgcm. The included matlab program {\it input/gendata.m} gives a complete
494	code for creating the {\it input/windx.sin\_y} file.
495
496	\subsubsection{File {\it input/topog.box}}
497
498
499	The {\it input/topog.box} file specifies a two-dimensional ($x,y$)
500	map of depth values. For this experiment values are either
501	$0m$ or $-2000\,{\rm m}$, corresponding respectively to a wall or to deep
502	ocean. The file contains a raw binary stream of data that is enumerated
503	in the same way as standard MITgcm two-dimensional, horizontal arrays.
504	The included matlab program {\it input/gendata.m} gives a complete
505	code for creating the {\it input/topog.box} file.
506
507	\subsubsection{File {\it code/SIZE.h}}
508
509	Two lines are customized in this file for the current experiment
510
511	\begin{itemize}
512
513	\item Line 39,
514	\begin{verbatim} sNx=60, \end{verbatim} this line sets
515	the lateral domain extent in grid points for the
516	axis aligned with the x-coordinate.
517
518	\item Line 40,
519	\begin{verbatim} sNy=60, \end{verbatim} this line sets
520	the lateral domain extent in grid points for the
521	axis aligned with the y-coordinate.
522
523	\item Line 49,
524	\begin{verbatim} Nr=4, \end{verbatim} this line sets
525	the vertical domain extent in grid points.
526
527	\end{itemize}
528
529	\begin{small}
530	% \include{code/SIZE.h}
531	\end{small}
532
533	\subsubsection{File {\it code/CPP\_OPTIONS.h}}
534
535	This file uses standard default values and does not contain
536	customizations for this experiment.
537
538
539	\subsubsection{File {\it code/CPP\_EEOPTIONS.h}}
540
541	This file uses standard default values and does not contain
542	customizations for this experiment.
543
544	\subsubsection{Other Files }
545
546	Other files relevant to this experiment are
547	\begin{itemize}
548	\item {\it model/src/ini\_cori.F}. This file initializes the model
549	coriolis variables {\bf fCorU}.
550	\item {\it model/src/ini\_spherical\_polar\_grid.F}
551	\item {\it model/src/ini\_parms.F},
552	\item {\it input/windx.sin\_y},
553	\end{itemize}
554	contain the code customizations and parameter settings for this
555	experiments. Below we describe the customizations
556	to these files associated with this experiment.