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.7 2006/06/27 20:47:05 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}
\begin{center}
(in directory: {\it verification/tutorial\_held\_suarez\_cs/})
\end{center}

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