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	cnh	1.6	% $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	cnh	1.2	% $Name: $
3	adcroft	1.1
4	cnh	1.6	\section{Centennial Time Scale Tracer Injection}
5			\label{sect:eg-simple-tracer}
6	adcroft	1.1
7			\bodytext{bgcolor="#FFFFFFFF"}
8
9			%\begin{center}
10	cnh	1.2	%{\Large \bf Using MITgcm to Look at Centennial Time Scale
11	adcroft	1.1	%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	cnh	1.2	This example illustrates the use of
23			the MITgcm to perform sensitivity analysis in a
24	adcroft	1.1	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	cnh	1.2	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	adcroft	1.1	Levitus data is also used throughout the calculation
44			to derive air-sea fluxes of heat at the ocean surface.
45	cnh	1.2	These fluxes are combined with climatological estimates of
46	adcroft	1.1	surface heat flux and fresh water, resulting in a mixed boundary
47	cnh	1.2	condition of the style described in Haney \cite{Haney}.
48	adcroft	1.1	Altogether, this yields the following forcing applied
49			in the model surface layer.
50
51			\begin{eqnarray}
52	cnh	1.6	\label{EQ:eg-simple-tracer-global_forcing}
53			\label{EQ:eg-simple-tracer-global_forcing_fu}
54	adcroft	1.1	{\cal F}_{u} & = & \frac{\tau_{x}}{\rho_{0} \Delta z_{s}}
55			\\
56	cnh	1.6	\label{EQ:eg-simple-tracer-global_forcing_fv}
57	adcroft	1.1	{\cal F}_{v} & = & \frac{\tau_{y}}{\rho_{0} \Delta z_{s}}
58			\\
59	cnh	1.6	\label{EQ:eg-simple-tracer-global_forcing_ft}
60	adcroft	1.1	{\cal F}_{\theta} & = & - \lambda_{\theta} ( \theta - \theta^{\ast} )
61			- \frac{1}{C_{p} \rho_{0} \Delta z_{s}}{\cal Q}
62			\\
63	cnh	1.6	\label{EQ:eg-simple-tracer-global_forcing_fs}
64	adcroft	1.1	{\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	adcroft	1.5	\cite{marshall:97a} is employed. A Laplacian operator, $\nabla^2$, provides viscous
123	cnh	1.2	dissipation. Thermal and haline diffusion is also represented by a Laplacian operator.
124	adcroft	1.1	\\
125
126			Wind-stress momentum inputs are added to the momentum equations for both
127	cnh	1.2	the zonal flow, $u$ and the meridional flow $v$, according to equations
128	cnh	1.6	(\ref{EQ:eg-simple-tracer-global_forcing_fu}) and (\ref{EQ:eg-simple-tracer-global_forcing_fv}).
129	adcroft	1.1	Thermodynamic forcing inputs are added to the equations for
130			potential temperature, $\theta$, and salinity, $S$, according to equations
131	cnh	1.6	(\ref{EQ:eg-simple-tracer-global_forcing_ft}) and (\ref{EQ:eg-simple-tracer-global_forcing_fs}).
132	adcroft	1.1	This produces a set of equations solved in this configuration as follows:
133			% {\fracktur}
134
135
136			\begin{eqnarray}
137	cnh	1.6	\label{EQ:eg-simple-tracer-model_equations}
138	adcroft	1.1	\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	adcroft	1.4	MITgcm Numerical Solution Procedure \ref{chap:discretization}, the time
172	adcroft	1.1	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	cnh	1.2	The Laplacian dissipation coefficient, $A_{h}$, is set to $400 m s^{-1}$.
180	adcroft	1.3	This value is chosen to yield a Munk layer width \cite{adcroft:95},
181	adcroft	1.1
182			\begin{eqnarray}
183	cnh	1.6	\label{EQ:eg-simple-tracer-munk_layer}
184	adcroft	1.1	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	adcroft	1.3	parameter to the horizontal Laplacian friction \cite{adcroft:95}
195	adcroft	1.1
196			\begin{eqnarray}
197	cnh	1.6	\label{EQ:eg-simple-tracer-laplacian_stability}
198	adcroft	1.1	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	cnh	1.6	\label{EQ:eg-simple-tracer-laplacian_stability_z}
210	adcroft	1.1	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	adcroft	1.3	\cite{adcroft:95}
221	adcroft	1.1
222			\begin{eqnarray}
223	cnh	1.6	\label{EQ:eg-simple-tracer-inertial_stability}
224	adcroft	1.1	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	adcroft	1.3	\noindent The advective CFL \cite{adcroft:95} for a extreme maximum
232	adcroft	1.1	horizontal flow
233			speed of $ \| \vec{u} \| = 2 ms^{-1}$
234
235			\begin{eqnarray}
236	cnh	1.6	\label{EQ:eg-simple-tracer-cfl_stability}
237	adcroft	1.1	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	cnh	1.2	\noindent The stability parameter for internal gravity waves
245	adcroft	1.3	\cite{adcroft:95}
246	adcroft	1.1
247			\begin{eqnarray}
248	cnh	1.6	\label{EQ:eg-simple-tracer-igw_stability}
249	adcroft	1.1	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	cnh	1.2	contain the code customizations and parameter settings for this
271			experiments. Below we describe the customizations
272	adcroft	1.1	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	cnh	1.2	depth level the initial and reference profiles will be uniform in
289	adcroft	1.1	$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	cnh	1.2	this line sets the vertical Laplacian dissipation coefficient to
301	adcroft	1.1	$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	cnh	1.2	this line sets the horizontal Laplacian frictional dissipation coefficient to
316	adcroft	1.1	$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	cnh	1.2	the horizontal Laplacian friction operator
325	adcroft	1.1	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	cnh	1.2	boundary condition in the vertical Laplacian friction operator
334	adcroft	1.1	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	cnh	1.2	grid input parameters, for example {\bf delX} and {\bf delY} and
387			causes the grid generation routines to initialize an internal grid based
388	adcroft	1.1	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	cnh	1.2	uses internally and affects the initialization of the coriolis force.
404	adcroft	1.1	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	cnh	1.2	customizations for this experiment.
480	adcroft	1.1
481			\subsubsection{File {\it input/eedata}}
482
483			This file uses standard default values and does not contain
484	cnh	1.2	customizations for this experiment.
485	adcroft	1.1
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	cnh	1.2	customizations for this experiment.
537	adcroft	1.1
538
539			\subsubsection{File {\it code/CPP\_EEOPTIONS.h}}
540
541			This file uses standard default values and does not contain
542	cnh	1.2	customizations for this experiment.
543	adcroft	1.1
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	cnh	1.2	contain the code customizations and parameter settings for this
555			experiments. Below we describe the customizations
556	adcroft	1.1	to these files associated with this experiment.