s_examples/held_suarez_cs/held_suarez_cs.tex

% $Header: /u/gcmpack/manual/part3/case_studies/held_suarez_cs/held_suarez_cs.tex,v 1.6 2006/06/27 19:08:22 molod Exp $
% $Name:  $

\section[Held-Suarez Atmosphere MITgcm Example]{Held-Suarez atmospheric simulation
on cube-sphere grid with 32 square cube faces.}
\label{www:tutorials}
\label{sect:eg-hs}
\begin{rawhtml}
<!-- CMIREDIR:eg-hs: -->
\end{rawhtml}

\bodytext{bgcolor="#FFFFFFFF"}

%\begin{center} 
%{\Large \bf Using MITgcm to Simulate Global Climatological Ocean Circulation
%At Four Degree Resolution with Asynchronous Time Stepping}
%
%\vspace*{4mm}
%
%\vspace*{3mm}
%{\large May 2001}
%\end{center}

This example illustrates the use of the MITgcm as an Atmospheric GCM,
using simple \cite{held-suar:94} forcing 
to simulate Atmospheric Dynamics on global scale.
The set-up use the rescaled pressure coordinate ($p^*$)\cite[]{adcroft:04a}
in the vertical direction, with 20 equaly-spaced levels, and
the conformal cube-sphere grid (C32) \cite[]{adcroft:04b}.
The files for this experiment can be found in the verification directory
under tutorial\_held\_suarez\_cs.

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

This example demonstrates using the MITgcm to simulate
the planetary atmospheric circulation, with flat orography
and simplified forcing.
In particular, only dry air processes are considered and 
radiation effects are represented by a simple newtownien cooling,
Thus this example does not rely on any particular atmospheric 
physics package.
This kind of simplified atmospheric simulation has been widely 
used in GFD-type experiments and in intercomparison projects of 
AGCM dynamical cores \cite[]{held-suar:94}.

The horizontal grid is obtain from the projection of a uniform gridded cube 
to the sphere. Each of the 6 faces has the same resolution, with 
$32 \times 32$ grid points. The equator line coincide with a grid line
and crosses, right in the midle, 4 of the 6 faces, leaving 2 faces 
for the Northern and Southern polar regions.
This curvilinear grid requires the use of the 2nd generation exchange
topology ({\it pkg/exch2}) to connect tile and face edges,
but without any limitation on the number of processors.

The use of the $p^*$ coordinate with 20 equally spaced levels
($20 \times 50\,{\rm mb}$, from $p^*=1000,{\rm mb}$ to $0$ at the 
top of the atmosphere) follows the choice of \cite{held-suar:94}. 
Note that without topography, the $p^*$ coordinate and the normalized 
pressure coordinate ($\sigma_p$) coincide exactly.
No viscosity and zero diffusion are used here, but
a $8^th$ order \cite{Shapiro_70} filter is applied to both momentum and
potential temperature, to remove selectively grid scale noise.
Apart from the horizontal grid, this experiment is made very similar to
the grid-point model case used in \cite{held-suar:94} study.

At this resolution, the configuration can be integrated forward 
for many years on a single processor desktop computer.
\\

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

The model is forced by relaxation to a radiative equilibrium temperature from
\cite{held-suar:94}.
A linear frictional drag (Rayleigh damping) is applied in the lower 
part of the atmosphere and account from surface friction and momentum
dissipation in the boundary layer.
Altogether, this yields the following forcing 
\cite[from][]{held-suar:94} that is applied to the fluid:

\begin{eqnarray}
\label{EQ:eg-hs-global_forcing}
\label{EQ:eg-hs-global_forcing_fv}
\vec{{\cal F}_\mathbf{v}} & = & -k_\mathbf{v}(p)\vec{\mathbf{v}}_h
\\
\label{EQ:eg-hs-global_forcing_ft}
{\cal F}_{\theta} & = & -k_{\theta}(\varphi,p)[\theta-\theta_{eq}(\varphi,p)]
\end{eqnarray}

\noindent where $\vec{\cal F}_\mathbf{v}$, ${\cal F}_{\theta}$,
are the forcing terms in the zonal and meridional
momentum and in the potential temperature equations respectively.
The term $k_\mathbf{v}$ in equation (\ref{EQ:eg-hs-global_forcing_fv}) applies a
Rayleigh damping that is active within the planetary boundary layer. 
It is defined so as to decay as pressure decreases according to
\begin{eqnarray*}
\label{EQ:eg-hs-define_kv}
k_\mathbf{v} & = & k_{f}~\max[0,~(p^*/P^{0}_{s}-\sigma_{b})/(1-\sigma_{b})]
\\
\sigma_{b} & = & 0.7 ~~{\rm and}~~
k_{f}  =  1/86400 ~{\rm s}^{-1}
\end{eqnarray*}

where $p^*$ is the pressure level of the cell center 
and $P^{0}_{s}$ is the pressure at the base of the atmospheric column,
which is constant and uniform here ($= 10^5 {\rm Pa}$), in the absence 
of topography.

The Equilibrium temperature $\theta_{eq}$ and relaxation time scale $k_{\theta}$ 
are set to:
\begin{eqnarray}
\label{EQ:eg-hs-define_Teq}
\theta_{eq}(\varphi,p^*) & = & \max \{ 200.K (P^{0}_{s}/p^*)^\kappa,\\
\nonumber
& & \hspace{8mm} 315.K - \Delta T_y~\sin^2(\varphi) 
  - \Delta \theta_z \cos^2(\varphi) \log(p^*/P^{0}_s) \}
\\
\label{EQ:eg-hs-define_kT}
k_{\theta}(\varphi,p^*) & = &
k_a + (k_s -k_a)~\cos^4(\varphi)~\max[0,(p^*/P^{0}_{s}-\sigma_{b})/(1-\sigma_{b})]
\end{eqnarray}
with:
\begin{eqnarray*}
 \Delta T_y = 60.K & k_a = 1/(40 \cdot 86400) ~{\rm s}^{-1}\\
\Delta \theta_z = 10.K & k_s = 1/(4 \cdot 86400) ~{\rm s}^{-1}
\end{eqnarray*}

Initial conditions correspond to a resting state with horizontally uniform 
stratified fluid. The initial temperature profile is simply the 
horizontally average of the radiative equilibrium temperature.

\subsection{Set-up description}
%\subsection{Discrete Numerical Configuration}
\label{www:tutorials}

The model is configured in hydrostatic form, using non-boussinesq
$p^*$ coordinate.
The vertical resolution is uniform, $\Delta p^* = 50.10^2 Pa$,
with 20 levels, from $p^*=10^5 Pa$ to $0$ at the top.
The domain is discretised using C32 cube-sphere grid \cite[]{adcroft:04b}
that cover the whole sphere with a relatively uniform grid-spacing.
The resolution at the equator or along the Greenwitch meridian
is similar to the $128 \times 64$ equaly spaced longitude-latitude grid,
but requires $25\%$ less grid points.
Grid spacing and grid-point location are not computed by the model but
read from files.

The vector-invariant form of the momentum equation (see section
\ref{sect:vect-inv_momentum_equations}) is used so that no explicit 
metrics are necessary.

Applying the vector-invariant discretization to the
atmospheric equations \ref{eq:atmos-prime}, and adding the 
forcing term 
(\ref{EQ:eg-hs-global_forcing_fv}, \ref{EQ:eg-hs-global_forcing_ft})
on the right-hand-side,
leads to the set of equations that are solved in this configuration:

%the The set of equations solved here is der
%Wind-stress forcing is added to the momentum equations for both
%the zonal flow, $u$ and the meridional flow $v$, according to equations 
%(\ref{EQ:eg-hs-global_forcing_fv}) and (\ref{EQ:eg-hs-global_forcing_fv}).
%Thermodynamic forcing inputs are added to the equations for
%potential temperature, $\theta$, and salinity, $S$, according to equations 
%(\ref{EQ:eg-hs-global_forcing_ft}) and (\ref{EQ:eg-hs-global_forcing_fs}).

\begin{eqnarray}
\label{EQ:eg-hs-model_equations}
\frac{\partial \vec{\mathbf{v}}_h}{\partial t}
+(f + \zeta)\hat{\mathbf{k}} \times \vec{\mathbf{v}}_h
%+\mathbf{\nabla }_{p} ({\rm KE})
+\mathbf{\nabla }_{p} (\mbox{\sc ke})
+ \omega \frac{\partial \vec{\mathbf{v}}_h }{\partial p}
+\mathbf{\nabla }_p \Phi ^{\prime } 
&=& 
% \vec{{\cal F}_v} = 
-k_\mathbf{v}\vec{\mathbf{v}}_h
\\
\frac{\partial \Phi ^{\prime }}{\partial p} 
+\frac{\partial \Pi }{\partial p}\theta ^{\prime } &=&0
\\
\mathbf{\nabla }_{p}\cdot \vec{\mathbf{v}}_h+\frac{\partial \omega }{
\partial p} &=&0
\\
\frac{\partial \theta }{\partial t} 
+ \mathbf{\nabla }_{p}\cdot (\theta \vec{\mathbf{v}}_h)
+ \frac{\partial (\theta \omega)}{\partial p}
%= \frac{\mathcal{Q}}{\Pi } 
&=& -k_{\theta}[\theta-\theta_{eq}]
\end{eqnarray}

%\begin{equation}
%\partial_t \vec{v} + ( 2\vec{\Omega} + \vec{\zeta}) \wedge \vec{v}
%- b \hat{r}
%+ \vec{\nabla} B = \vec{\nabla} \cdot \vec{\bf \tau}
%\end{equation}
%{\cal F}_{\theta} & = & -k_{\theta}(\varphi,p)[\theta-\theta_{eq}(\varphi,p)]

\noindent where $\vec{\mathbf{v}}_h$ and $\omega = \frac{Dp}{Dt}$ 
are the horizontal velocity vector and the vertical velocity in pressure coordinate,
$\zeta$ is the relative vorticity and $f$ the Coriolis parameter, 
$\hat{\mathbf{k}}$ is the vertical unity vector, 
{\sc ke} is the kinetic energy, $\Phi$ is the geopotential
and $\Pi$ the Exner function 
($\Pi = C_p (p/p_c)^\kappa ~{\rm with}~ p_c = 10^5 Pa$).
Variables marked with ' corresponds to anomaly from 
the resting, uniformly stratified state.

As described in MITgcm Numerical Solution Procedure \ref{chap:discretization}, 
the continuity equation is integrated vertically, to give a prognostic
equation for the surface pressure $p_s$:
\begin{equation}
\frac{\partial p_s}{\partial t} + \nabla_{h}\cdot \int_{0}^{p_s} \vec{\mathbf{v}}_h dp
= 0
\end{equation}

The implicit free surface form of the pressure equation described in 
\cite{marshall:97a} is employed to solve for $p_s$; 
Integrating vertically the hydrostatic balance 
gives the geopotential $\Phi'$ and allow to step forward the momentum equation
\ref{EQ:eg-hs-model_equations}.
The potential temperature, $\theta$, is stepped forward using the 
new velocity field ({\it staggered time-step}, section 
\ref{sect:adams-bashforth-staggered}).
\\

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

\noindent The numerical stability for inertial oscillations
\cite{adcroft:95} 

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

\noindent evaluates to $4.\times10^{-3}$ at the poles, 
for $f=2\Omega\sin(\pi / 2) =1.45\times10^{-4}~{\rm s}^{-1}$, 
which is well below the $S_{i} < 1$ upper limit for stability.
\\

\noindent The advective CFL \cite{adcroft:95} 
for a extreme maximum horizontal flow speed of $ | \vec{u} | = 90. {\rm m/s}$~ 
and the smallest horizontal grid spacing $ \Delta x = 1.1\times10^5 {\rm m}$~:

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

\noindent evaluates to $0.37$, which is close to the stability 
limit of 0.5.
\\

\noindent The stability parameter for internal gravity waves propagating
with a maximum speed of $c_{g}=100~{\rm m/s}$
\cite{adcroft:95}

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

\noindent evaluates to $4 \times 10^{-1}$. This is close to the linear
stability limit of 0.5.
  
\subsection{Experiment Configuration}
\label{www:tutorials}
\label{SEC:eg-hs_examp_exp_config}

The model configuration for this experiment resides under the 
directory {\it verification/tutorial\_held\_suarez\_cs}.  The experiment files 
\begin{itemize}
\item {\it input/data}
\item {\it input/data.pkg}
\item {\it input/eedata},
\item {\it input/data.shap},
\item {\it code/packages.conf},
\item {\it code/CPP\_OPTIONS.h},
\item {\it code/SIZE.h},
\item {\it code/DIAGNOSTICS\_SIZE.h},
\item {\it code/external\_forcing.F},
\end{itemize}
contain the code customizations and parameter settings for these
experiments. Below we describe the customizations
to these files associated with this experiment.

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

\input{part3/case_studies/held_suarez_cs/inp_data}

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

\input{part3/case_studies/held_suarez_cs/inp_data.pkg}

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

This file uses standard default values except line 6:
\begin{verbatim}
 useCubedSphereExchange=.TRUE.,
\end{verbatim}
This line selects the cubed-sphere specific exchanges to
to connect tiles and faces edges.

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

\input{part3/case_studies/held_suarez_cs/inp_data.shap}

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

Four lines are customized in this file for the current experiment

\begin{itemize}

\item Line 39, 
\begin{verbatim} sNx=32, \end{verbatim}
sets the lateral domain extent in grid points along the x-direction,
for 1 face.

\item Line 40,
\begin{verbatim} sNy=32, \end{verbatim} 
sets the lateral domain extent in grid points along the y-direction,
for 1 face.

\item Line 43,
\begin{verbatim} nSx=6, \end{verbatim} 
sets the number of tiles in the x-direction, for the model domain
decomposition. In this simple case (one processor and 1 tile per 
face), this number corresponds to the total number of faces.

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

\end{itemize}

%\begin{small}
%\input{part3/case_studies/held_suarez_cs/code/SIZE.h}
%\end{small}

\subsubsection{File {\it code/packages.conf}}
\label{www:tutorials}

\input{part3/case_studies/held_suarez_cs/cod_packages.conf}

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

This file uses standard default except for Line 40\\
({\it diff CPP\_OPTIONS.h ../../../model/inc}):
\begin{verbatim}
#define NONLIN_FRSURF 
\end{verbatim}
This line allow to use the non-linear free-surface part of the code,
which is required for the $p^*$ coordinate formulation.

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

Other files relevant to this experiment are
\begin{itemize}
\item {\it code/external\_forcing.F}
\item {\it input/grid\_cs32.face00[n].bin}, with $n=1,2,3,4,5,6$
\end{itemize}
contain the code customisations and binary input files for this 
experiments. Below we describe the customisations
to these files associated with this experiment.\\

The file {\it code/external\_forcing.F} contains 4 subroutines
that calculate the forcing terms (Right-Hand side term) in the
momentum equation (\ref{EQ:eg-hs-global_forcing_fv}, 
{\it S/R EXTERNAL\_FORCING\_U} and {\it EXTERNAL\_FORCING\_V})
and in the potential temperature equation 
(\ref{EQ:eg-hs-global_forcing_ft}, {\it S/R  EXTERNAL\_FORCING\_T}). 
The water-vapour forcing subroutine ({\it S/R EXTERNAL\_FORCING\_S})
is left empty for this experiment.\\

The grid-files {\it input/grid\_cs32.face00[n].bin}, with $n=1,2,3,4,5,6$,
are binary files (direct-access, big-endian 64.bits real) that 
contains all the cubed-sphere grid lengths, areas and grid-point
positions, with one file per face.
Each file contains 18 2-D arrays (dimension $33 \times 33$) that correspond
to the model variables:
{\it 
XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG AngleCS AngleSN 
}
(see {\it GRID.h} file)

1	molod	1.7	% $Header: /u/gcmpack/manual/part3/case_studies/held_suarez_cs/held_suarez_cs.tex,v 1.6 2006/06/27 19:08:22 molod Exp $
2	jmc	1.1	% $Name: $
3
4			\section[Held-Suarez Atmosphere MITgcm Example]{Held-Suarez atmospheric simulation
5			on cube-sphere grid with 32 square cube faces.}
6			\label{www:tutorials}
7			\label{sect:eg-hs}
8			\begin{rawhtml}
9			<!-- CMIREDIR:eg-hs: -->
10			\end{rawhtml}
11
12			\bodytext{bgcolor="#FFFFFFFF"}
13
14			%\begin{center}
15			%{\Large \bf Using MITgcm to Simulate Global Climatological Ocean Circulation
16			%At Four Degree Resolution with Asynchronous Time Stepping}
17			%
18			%\vspace*{4mm}
19			%
20			%\vspace*{3mm}
21			%{\large May 2001}
22			%\end{center}
23
24			This example illustrates the use of the MITgcm as an Atmospheric GCM,
25	jmc	1.2	using simple \cite{held-suar:94} forcing
26	jmc	1.1	to simulate Atmospheric Dynamics on global scale.
27			The set-up use the rescaled pressure coordinate ($p^*$)\cite[]{adcroft:04a}
28			in the vertical direction, with 20 equaly-spaced levels, and
29			the conformal cube-sphere grid (C32) \cite[]{adcroft:04b}.
30	molod	1.6	The files for this experiment can be found in the verification directory
31	molod	1.7	under tutorial\_held\_suarez\_cs.
32	jmc	1.1
33			\subsection{Overview}
34			\label{www:tutorials}
35
36			This example demonstrates using the MITgcm to simulate
37			the planetary atmospheric circulation, with flat orography
38			and simplified forcing.
39			In particular, only dry air processes are considered and
40			radiation effects are represented by a simple newtownien cooling,
41	jmc	1.4	Thus this example does not rely on any particular atmospheric
42	jmc	1.1	physics package.
43			This kind of simplified atmospheric simulation has been widely
44			used in GFD-type experiments and in intercomparison projects of
45			AGCM dynamical cores \cite[]{held-suar:94}.
46
47			The horizontal grid is obtain from the projection of a uniform gridded cube
48			to the sphere. Each of the 6 faces has the same resolution, with
49			$32 \times 32$ grid points. The equator line coincide with a grid line
50			and crosses, right in the midle, 4 of the 6 faces, leaving 2 faces
51			for the Northern and Southern polar regions.
52			This curvilinear grid requires the use of the 2nd generation exchange
53			topology ({\it pkg/exch2}) to connect tile and face edges,
54			but without any limitation on the number of processors.
55
56			The use of the $p^*$ coordinate with 20 equally spaced levels
57			($20 \times 50\,{\rm mb}$, from $p^*=1000,{\rm mb}$ to $0$ at the
58	jmc	1.2	top of the atmosphere) follows the choice of \cite{held-suar:94}.
59			Note that without topography, the $p^*$ coordinate and the normalized
60			pressure coordinate ($\sigma_p$) coincide exactly.
61			No viscosity and zero diffusion are used here, but
62			a $8^th$ order \cite{Shapiro_70} filter is applied to both momentum and
63			potential temperature, to remove selectively grid scale noise.
64			Apart from the horizontal grid, this experiment is made very similar to
65			the grid-point model case used in \cite{held-suar:94} study.
66	jmc	1.1
67			At this resolution, the configuration can be integrated forward
68			for many years on a single processor desktop computer.
69			\\
70
71	jmc	1.2	\subsection{Forcing}
72			\label{www:tutorials}
73
74	jmc	1.1	The model is forced by relaxation to a radiative equilibrium temperature from
75			\cite{held-suar:94}.
76			A linear frictional drag (Rayleigh damping) is applied in the lower
77			part of the atmosphere and account from surface friction and momentum
78			dissipation in the boundary layer.
79			Altogether, this yields the following forcing
80			\cite[from][]{held-suar:94} that is applied to the fluid:
81
82			\begin{eqnarray}
83			\label{EQ:eg-hs-global_forcing}
84	jmc	1.2	\label{EQ:eg-hs-global_forcing_fv}
85	jmc	1.4	\vec{{\cal F}_\mathbf{v}} & = & -k_\mathbf{v}(p)\vec{\mathbf{v}}_h
86	jmc	1.1	\\
87			\label{EQ:eg-hs-global_forcing_ft}
88	jmc	1.4	{\cal F}_{\theta} & = & -k_{\theta}(\varphi,p)[\theta-\theta_{eq}(\varphi,p)]
89	jmc	1.1	\end{eqnarray}
90
91	jmc	1.4	\noindent where $\vec{\cal F}_\mathbf{v}$, ${\cal F}_{\theta}$,
92	jmc	1.1	are the forcing terms in the zonal and meridional
93			momentum and in the potential temperature equations respectively.
94	jmc	1.4	The term $k_\mathbf{v}$ in equation (\ref{EQ:eg-hs-global_forcing_fv}) applies a
95	jmc	1.2	Rayleigh damping that is active within the planetary boundary layer.
96			It is defined so as to decay as pressure decreases according to
97			\begin{eqnarray*}
98	jmc	1.1	\label{EQ:eg-hs-define_kv}
99	jmc	1.4	k_\mathbf{v} & = & k_{f}~\max[0,~(p^*/P^{0}_{s}-\sigma_{b})/(1-\sigma_{b})]
100	jmc	1.1	\\
101	jmc	1.2	\sigma_{b} & = & 0.7 ~~{\rm and}~~
102			k_{f} = 1/86400 ~{\rm s}^{-1}
103			\end{eqnarray*}
104	jmc	1.1
105			where $p^*$ is the pressure level of the cell center
106	jmc	1.2	and $P^{0}_{s}$ is the pressure at the base of the atmospheric column,
107	jmc	1.4	which is constant and uniform here ($= 10^5 {\rm Pa}$), in the absence
108	jmc	1.2	of topography.
109	jmc	1.1
110			The Equilibrium temperature $\theta_{eq}$ and relaxation time scale $k_{\theta}$
111			are set to:
112			\begin{eqnarray}
113	jmc	1.2	\label{EQ:eg-hs-define_Teq}
114	jmc	1.4	\theta_{eq}(\varphi,p^) & = & \max \{ 200.K (P^{0}_{s}/p^)^\kappa,\\
115	jmc	1.2	\nonumber
116	jmc	1.4	& & \hspace{8mm} 315.K - \Delta T_y~\sin^2(\varphi)
117			- \Delta \theta_z \cos^2(\varphi) \log(p^*/P^{0}_s) \}
118	jmc	1.2	\\
119	jmc	1.1	\label{EQ:eg-hs-define_kT}
120	jmc	1.4	k_{\theta}(\varphi,p^*) & = &
121			k_a + (k_s -k_a)~\cos^4(\varphi)~\max[0,(p^*/P^{0}_{s}-\sigma_{b})/(1-\sigma_{b})]
122	jmc	1.1	\end{eqnarray}
123	jmc	1.2	with:
124			\begin{eqnarray*}
125			\Delta T_y = 60.K & k_a = 1/(40 \cdot 86400) ~{\rm s}^{-1}\\
126			\Delta \theta_z = 10.K & k_s = 1/(4 \cdot 86400) ~{\rm s}^{-1}
127			\end{eqnarray*}
128	jmc	1.1
129			Initial conditions correspond to a resting state with horizontally uniform
130			stratified fluid. The initial temperature profile is simply the
131			horizontally average of the radiative equilibrium temperature.
132
133	jmc	1.3	\subsection{Set-up description}
134			%\subsection{Discrete Numerical Configuration}
135	jmc	1.1	\label{www:tutorials}
136
137	jmc	1.2	The model is configured in hydrostatic form, using non-boussinesq
138			$p^*$ coordinate.
139			The vertical resolution is uniform, $\Delta p^* = 50.10^2 Pa$,
140			with 20 levels, from $p^*=10^5 Pa$ to $0$ at the top.
141			The domain is discretised using C32 cube-sphere grid \cite[]{adcroft:04b}
142			that cover the whole sphere with a relatively uniform grid-spacing.
143			The resolution at the equator or along the Greenwitch meridian
144			is similar to the $128 \times 64$ equaly spaced longitude-latitude grid,
145			but requires $25\%$ less grid points.
146			Grid spacing and grid-point location are not computed by the model but
147			read from files.
148
149			The vector-invariant form of the momentum equation (see section
150			\ref{sect:vect-inv_momentum_equations}) is used so that no explicit
151			metrics are necessary.
152
153			Applying the vector-invariant discretization to the
154			atmospheric equations \ref{eq:atmos-prime}, and adding the
155			forcing term
156			(\ref{EQ:eg-hs-global_forcing_fv}, \ref{EQ:eg-hs-global_forcing_ft})
157			on the right-hand-side,
158			leads to the set of equations that are solved in this configuration:
159
160			%the The set of equations solved here is der
161			%Wind-stress forcing is added to the momentum equations for both
162			%the zonal flow, $u$ and the meridional flow $v$, according to equations
163			%(\ref{EQ:eg-hs-global_forcing_fv}) and (\ref{EQ:eg-hs-global_forcing_fv}).
164			%Thermodynamic forcing inputs are added to the equations for
165			%potential temperature, $\theta$, and salinity, $S$, according to equations
166			%(\ref{EQ:eg-hs-global_forcing_ft}) and (\ref{EQ:eg-hs-global_forcing_fs}).
167	jmc	1.1
168			\begin{eqnarray}
169			\label{EQ:eg-hs-model_equations}
170	jmc	1.2	\frac{\partial \vec{\mathbf{v}}_h}{\partial t}
171			+(f + \zeta)\hat{\mathbf{k}} \times \vec{\mathbf{v}}_h
172			%+\mathbf{\nabla }_{p} ({\rm KE})
173			+\mathbf{\nabla }_{p} (\mbox{\sc ke})
174			+ \omega \frac{\partial \vec{\mathbf{v}}_h }{\partial p}
175	jmc	1.5	+\mathbf{\nabla }_p \Phi ^{\prime }
176	jmc	1.2	&=&
177			% \vec{{\cal F}_v} =
178	jmc	1.4	-k_\mathbf{v}\vec{\mathbf{v}}_h
179	jmc	1.2	\\
180	jmc	1.5	\frac{\partial \Phi ^{\prime }}{\partial p}
181	jmc	1.2	+\frac{\partial \Pi }{\partial p}\theta ^{\prime } &=&0
182			\\
183			\mathbf{\nabla }_{p}\cdot \vec{\mathbf{v}}_h+\frac{\partial \omega }{
184			\partial p} &=&0
185			\\
186			\frac{\partial \theta }{\partial t}
187			+ \mathbf{\nabla }_{p}\cdot (\theta \vec{\mathbf{v}}_h)
188			+ \frac{\partial (\theta \omega)}{\partial p}
189			%= \frac{\mathcal{Q}}{\Pi }
190			&=& -k_{\theta}[\theta-\theta_{eq}]
191			\end{eqnarray}
192
193			%\begin{equation}
194			%\partial_t \vec{v} + ( 2\vec{\Omega} + \vec{\zeta}) \wedge \vec{v}
195			%- b \hat{r}
196			%+ \vec{\nabla} B = \vec{\nabla} \cdot \vec{\bf \tau}
197			%\end{equation}
198	jmc	1.4	%{\cal F}_{\theta} & = & -k_{\theta}(\varphi,p)[\theta-\theta_{eq}(\varphi,p)]
199	jmc	1.2
200	jmc	1.4	\noindent where $\vec{\mathbf{v}}_h$ and $\omega = \frac{Dp}{Dt}$
201			are the horizontal velocity vector and the vertical velocity in pressure coordinate,
202	jmc	1.2	$\zeta$ is the relative vorticity and $f$ the Coriolis parameter,
203			$\hat{\mathbf{k}}$ is the vertical unity vector,
204			{\sc ke} is the kinetic energy, $\Phi$ is the geopotential
205			and $\Pi$ the Exner function
206			($\Pi = C_p (p/p_c)^\kappa ~{\rm with}~ p_c = 10^5 Pa$).
207	jmc	1.4	Variables marked with ' corresponds to anomaly from
208	jmc	1.2	the resting, uniformly stratified state.
209
210			As described in MITgcm Numerical Solution Procedure \ref{chap:discretization},
211			the continuity equation is integrated vertically, to give a prognostic
212			equation for the surface pressure $p_s$:
213			\begin{equation}
214	jmc	1.4	\frac{\partial p_s}{\partial t} + \nabla_{h}\cdot \int_{0}^{p_s} \vec{\mathbf{v}}_h dp
215	jmc	1.2	= 0
216			\end{equation}
217
218			The implicit free surface form of the pressure equation described in
219			\cite{marshall:97a} is employed to solve for $p_s$;
220			Integrating vertically the hydrostatic balance
221	jmc	1.4	gives the geopotential $\Phi'$ and allow to step forward the momentum equation
222	jmc	1.2	\ref{EQ:eg-hs-model_equations}.
223			The potential temperature, $\theta$, is stepped forward using the
224			new velocity field ({\it staggered time-step}, section
225			\ref{sect:adams-bashforth-staggered}).
226	jmc	1.1	\\
227
228			\subsubsection{Numerical Stability Criteria}
229			\label{www:tutorials}
230
231			\noindent The numerical stability for inertial oscillations
232			\cite{adcroft:95}
233
234			\begin{eqnarray}
235			\label{EQ:eg-hs-inertial_stability}
236	jmc	1.3	S_{i} = f^{2} {\Delta t}^2
237	jmc	1.1	\end{eqnarray}
238
239	jmc	1.3	\noindent evaluates to $4.\times10^{-3}$ at the poles,
240			for $f=2\Omega\sin(\pi / 2) =1.45\times10^{-4}~{\rm s}^{-1}$,
241			which is well below the $S_{i} < 1$ upper limit for stability.
242	jmc	1.1	\\
243
244	jmc	1.3	\noindent The advective CFL \cite{adcroft:95}
245			for a extreme maximum horizontal flow speed of $ \| \vec{u} \| = 90. {\rm m/s}$~
246			and the smallest horizontal grid spacing $ \Delta x = 1.1\times10^5 {\rm m}$~:
247	jmc	1.1
248			\begin{eqnarray}
249			\label{EQ:eg-hs-cfl_stability}
250	jmc	1.3	S_{a} = \frac{\| \vec{u} \| \Delta t}{ \Delta x}
251	jmc	1.1	\end{eqnarray}
252
253	jmc	1.3	\noindent evaluates to $0.37$, which is close to the stability
254	jmc	1.1	limit of 0.5.
255			\\
256
257			\noindent The stability parameter for internal gravity waves propagating
258	jmc	1.3	with a maximum speed of $c_{g}=100~{\rm m/s}$
259	jmc	1.1	\cite{adcroft:95}
260
261			\begin{eqnarray}
262			\label{EQ:eg-hs-gfl_stability}
263	jmc	1.3	S_{c} = \frac{c_{g} \Delta t}{ \Delta x}
264	jmc	1.1	\end{eqnarray}
265
266	jmc	1.3	\noindent evaluates to $4 \times 10^{-1}$. This is close to the linear
267	jmc	1.1	stability limit of 0.5.
268
269			\subsection{Experiment Configuration}
270			\label{www:tutorials}
271			\label{SEC:eg-hs_examp_exp_config}
272
273			The model configuration for this experiment resides under the
274	jmc	1.4	directory {\it verification/tutorial\_held\_suarez\_cs}. The experiment files
275	jmc	1.1	\begin{itemize}
276			\item {\it input/data}
277			\item {\it input/data.pkg}
278			\item {\it input/eedata},
279	jmc	1.3	\item {\it input/data.shap},
280			\item {\it code/packages.conf},
281	jmc	1.1	\item {\it code/CPP\_OPTIONS.h},
282	jmc	1.3	\item {\it code/SIZE.h},
283			\item {\it code/DIAGNOSTICS\_SIZE.h},
284			\item {\it code/external\_forcing.F},
285	jmc	1.1	\end{itemize}
286			contain the code customizations and parameter settings for these
287			experiments. Below we describe the customizations
288			to these files associated with this experiment.
289
290			\subsubsection{File {\it input/data}}
291			\label{www:tutorials}
292
293	jmc	1.3	\input{part3/case_studies/held_suarez_cs/inp_data}
294	jmc	1.1
295			\subsubsection{File {\it input/data.pkg}}
296			\label{www:tutorials}
297
298	jmc	1.3	\input{part3/case_studies/held_suarez_cs/inp_data.pkg}
299
300	jmc	1.1	\subsubsection{File {\it input/eedata}}
301			\label{www:tutorials}
302
303	jmc	1.3	This file uses standard default values except line 6:
304			\begin{verbatim}
305			useCubedSphereExchange=.TRUE.,
306			\end{verbatim}
307			This line selects the cubed-sphere specific exchanges to
308			to connect tiles and faces edges.
309	jmc	1.1
310	jmc	1.3	\subsubsection{File {\it input/data.shap}}
311	jmc	1.1	\label{www:tutorials}
312
313	jmc	1.3	\input{part3/case_studies/held_suarez_cs/inp_data.shap}
314	jmc	1.1
315			\subsubsection{File {\it code/SIZE.h}}
316			\label{www:tutorials}
317
318	jmc	1.3	Four lines are customized in this file for the current experiment
319	jmc	1.1
320			\begin{itemize}
321
322			\item Line 39,
323	jmc	1.3	\begin{verbatim} sNx=32, \end{verbatim}
324	jmc	1.4	sets the lateral domain extent in grid points along the x-direction,
325	jmc	1.3	for 1 face.
326
327			\item Line 40,
328			\begin{verbatim} sNy=32, \end{verbatim}
329	jmc	1.4	sets the lateral domain extent in grid points along the y-direction,
330	jmc	1.3	for 1 face.
331
332			\item Line 43,
333			\begin{verbatim} nSx=6, \end{verbatim}
334	jmc	1.4	sets the number of tiles in the x-direction, for the model domain
335	jmc	1.3	decomposition. In this simple case (one processor and 1 tile per
336	jmc	1.4	face), this number corresponds to the total number of faces.
337	jmc	1.1
338			\item Line 49,
339	jmc	1.3	\begin{verbatim} Nr=20, \end{verbatim}
340			sets the vertical domain extent in grid points.
341	jmc	1.1
342			\end{itemize}
343
344	jmc	1.3	%\begin{small}
345			%\input{part3/case_studies/held_suarez_cs/code/SIZE.h}
346			%\end{small}
347	jmc	1.1
348	jmc	1.3	\subsubsection{File {\it code/packages.conf}}
349	jmc	1.1	\label{www:tutorials}
350
351	jmc	1.3	\input{part3/case_studies/held_suarez_cs/cod_packages.conf}
352	jmc	1.1
353	jmc	1.3	\subsubsection{File {\it code/CPP\_OPTIONS.h}}
354	jmc	1.1	\label{www:tutorials}
355
356	jmc	1.3	This file uses standard default except for Line 40\\
357			({\it diff CPP\_OPTIONS.h ../../../model/inc}):
358			\begin{verbatim}
359			#define NONLIN_FRSURF
360			\end{verbatim}
361			This line allow to use the non-linear free-surface part of the code,
362			which is required for the $p^*$ coordinate formulation.
363	jmc	1.1
364			\subsubsection{Other Files }
365			\label{www:tutorials}
366
367			Other files relevant to this experiment are
368			\begin{itemize}
369	jmc	1.3	\item {\it code/external\_forcing.F}
370			\item {\it input/grid\_cs32.face00[n].bin}, with $n=1,2,3,4,5,6$
371	jmc	1.1	\end{itemize}
372	jmc	1.3	contain the code customisations and binary input files for this
373	jmc	1.1	experiments. Below we describe the customisations
374	jmc	1.3	to these files associated with this experiment.\\
375
376			The file {\it code/external\_forcing.F} contains 4 subroutines
377			that calculate the forcing terms (Right-Hand side term) in the
378			momentum equation (\ref{EQ:eg-hs-global_forcing_fv},
379			{\it S/R EXTERNAL\_FORCING\_U} and {\it EXTERNAL\_FORCING\_V})
380			and in the potential temperature equation
381			(\ref{EQ:eg-hs-global_forcing_ft}, {\it S/R EXTERNAL\_FORCING\_T}).
382			The water-vapour forcing subroutine ({\it S/R EXTERNAL\_FORCING\_S})
383			is left empty for this experiment.\\
384
385			The grid-files {\it input/grid\_cs32.face00[n].bin}, with $n=1,2,3,4,5,6$,
386			are binary files (direct-access, big-endian 64.bits real) that
387			contains all the cubed-sphere grid lengths, areas and grid-point
388			positions, with one file per face.
389	jmc	1.4	Each file contains 18 2-D arrays (dimension $33 \times 33$) that correspond
390	jmc	1.3	to the model variables:
391			{\it
392			XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG AngleCS AngleSN
393			}
394			(see {\it GRID.h} file)
395