\documentstyle[A4,11pt]{article}

\begin{document}

\title{Transport de l'\'energie dans le mod\`ele du LMD}
\author{ Fr\'ed\'eric Hourdin}
\maketitle

\bibliographystyle{abbrv}

% cross-references
\def\sc#1{Section~\ref{sc:#1}}
\def\an#1{Appendix~\ref{an:#1}}
\def\ch#1{Chapter~\ref{ch:#1}}
\def\fig#1{Fig.~\ref{fg:#1}}
\def\fg#1{Fig.~\ref{fg:#1}}
\def\figs#1{Figs.~\ref{fg:#1}}
\def\eq#1{Eq.~\ref{eq:#1}}
\def\eqs#1{Eqs.~\ref{eq:#1}}
\def\tb#1{Table~\ref{tb:#1}}


\def\localdir{.}
\def\rcp{\kappa}
\def\uabs{\tilde{u}_{abs}}
\def\err{\epsilon}
\def\s{ {\dep{\frac{p}{\p0}}}^{\rcp} }
\def\dsig{\dz \sigma}
\def\ps{p_s}
\def\psk{{\ps}^\kappa}
\def\ucov{\tilde{u}}
\def\vcov{\tilde{v}}
\def\ucont{\tilde{\ucov}}
\def\vcont{\tilde{\vcov}}
\def\cu{c_u}
\def\cv{c_v}
\def\h{\theta}
\def\pext{\tilde{p}_s}
\def\fext{f}
\def\dx{\delta_x}
\def\dy{\delta_y}
\def\dz{\delta_z}
\def\K{\frac{1}{2}
\left( \av{\ucov \ucont}{X} + \av{\vcov \vcont}{Y} \right)}
\def\Z{\frac{\dx \vcov - \dy \ucov + \fext}{\av{\pext}{X,Y}}}
\def\Zm{\frac{- \dy \ucov + \fext}{\av{\pext}{Y}}}

\newcommand{\glob}[1]{ \left< #1 \right> }
\newcommand{\dt}[1]{\frac{\partial #1 }{\partial t} }
\newcommand{\av}[2]{{\overline{#1}}^{ #2 }}
\newcommand{\avg}[1]{\left< #1 \right>}
% \newcommand{\emis}[2]{E \left( #1 , #2 \right) }
\newcommand{\der}[2]{\frac{\partial #1 }{\partial #2} }
\newcommand{\dep}[1]{\left( #1 \right) }
\newcommand{\depb}[1]{\left[ #1 \right] }
\newcommand{\resol}[3]{ {\em #1 longitudes $\times$  #2 latitudes $\times$ #3  niveaux verticaux}}
\def\t{\theta}
\def\cd{C_D}

\section{Rappel des \'equations anciennes du MCG du LMD} 


La coordonn\'ee  verticale du mod\`ele est la pression normalis\'ee
par sa valeur \`a la surface: $\sigma=p/\ps$.
On utilise en fait $\sigma$ aux niveaux inter-couches
et $s=\sigma^\kappa$ au milieu des couches.
On note $X$ et $Y$ les coordonn\'ees horizontales:
$X$ (resp. $Y$) est une fonction biunivoque de la longitude $\lambda$
(resp. de la latitude $\phi$).
Les variables scalaires
(hentalpie potentielle $\h = c_p T/\psk$, g\'eopotentiel $\Phi$
et pression de surface $\ps$) sont \'evalu\'ees aux points
correspondant \`a des couples de valeurs enti\`eres $(X,Y)=(i,j)$.
Les variables dynamiques sont d\'ecal\'ees par rapport aux variables
scalaires en utilisant une grille $C$ dans la d\'efinition de
Arakawa~\cite{Arak:77}: le vent zonal est calcul\'e
aux points $(X,Y)=(i+1/2,j)$ et  le vent
m\'eridien aux points $(X,Y)=(i,j+1/2)$.
La disposition des variables sur la grille est illustr\'ee sur la
Figure~\ref{fg:grille}.
\input{\localdir/fggrille}

On utilise en fait les composantes covariantes
($\ucov$ et $\vcov$) et contravariantes ($\ucont$ et $\vcont$)
du vent d\'efinies par
\begin{equation}
\begin{array}{llllllllll}
\ucov = \cu u & \mbox{et} & \ucont = u / \cu & \mbox{avec} & 
\cu = a \cos{\phi} \left( d\lambda/dX \right)  \\
\vcov = \cv v & \mbox{et} & \vcont = v / \cv & \mbox{avec} &
\cv = a \left( d\phi / dY \right)
\end{array}
\end{equation}
%
o\`u $u$ et $v$ sont les composantes physiques du vecteur vent
horizontal.
On introduit \'egalement:
%
\paragraph{la pression extensive:}
$\pext$ (pression au sol multipli\'ee
par l'aire de la maille).
%
\paragraph{les trois composantes du flux de masse:}
\begin{equation}
U=\av{\pext}{X} \ucont ,\  V= \av{\pext}{Y} \vcont \  \mbox{et} \  
W= \pext \dot{\sigma}
\ \mbox{avec}\ \dot{\sigma}=\frac{d\sigma}{dt}
\end{equation}
%
\paragraph{le facteur de Coriolis multipli\'e par l'aire de la maille:}
$\fext=2\Omega \sin{\phi} \cu \cv$\\
o\`u $\Omega$ est la vitesse de rotation de la plan\`ete.
%
\paragraph{la vorticit\'e potentielle absolue:}
\begin{equation}
Z=\Z
\end{equation}
%
\paragraph{l'\'energie cin\'etique}
\begin{equation}
K=\K
\end{equation}\\
%
La notation $\delta_X$ signifie simplement qu'on
effectue la diff\'erence entre deux points cons\'ecutifs
suivant la direction $X$.
La notation $\av{a}{X}$ signifie qu'on prend la moyenne arithm\'etique
de la quantit\'e $a$ suivant la direction $X$.
Les \'equations discr\'etis\'ees sont \'ecrites sous la forme
suivante:
\paragraph{\'equations du mouvement:}
\begin{equation} \label{eq:u1}
\dt{\ucov} -
\av{Z}{Y} \av{V}{X,Y} 
+ \dx \left( \Phi + K \right)
+s \av{\h}{X} \dx \psk
+\frac{\av{\av{W}{X} \dz \ucov}{Z}}{\av{\pext}{X}\dz \sigma}
=S_{\ucov}
\end{equation}
\begin{equation} \label{eq:v1}
\dt{\vcov} + \av{Z}{X} \av{U}{X,Y} + \dy \left( \Phi + K \right)
+s \av{\h}{Y} \dy \psk
+\frac{\av{\av{W}{Y} \dz \vcov}{Z}}{\av{\pext}{Y}\dz \sigma}
=S_{\vcov}
\end{equation}
%
\paragraph{\'equation thermodynamique:}
%
\begin{equation}
\label{eq:thermo}
\dt{\dep{\pext \h}}
+\dx \dep{\av{\h}{X}U} +\dy \dep{\av{\h}{Y}V}
+\frac{\dz \dep{\av{\h}{Z} W}}{\dz \sigma}=S_\h
\end{equation}
%
\paragraph{\'equation hydrostatique:}
\begin{equation}
\label{eq:hydrostat}
\dz \Phi=-\psk  \av{\h}{z} \dz s
\end{equation}
with $\av{\h}{z}_1=\h_1$ and $\dz s_1=1-s_1$.
%
\paragraph{\'equations de continuit\'e:}
%
\begin{equation}
\label{eq:cont}
\dt{\pext} + \dx U+ \dy V +\frac{\dz W}{\dz\sigma}=0
\end{equation}
%
On a not\'e $S$ les termes sources dans les diff\'erentes \'equations.

\paragraph{Remarque:} on utilise ici la formulation originale, sans
contrainte sur la conservation du moment cin\'etique. Cette formulation
est plus proche d'une conservation num\'erique de l'\'energie.

\section{Numerical conservation of the total energy}

The specific energy is the sum of the specific kinetic energy $K$,
internal energy $I=C_v T$ and potential energy $\Phi$.

\subsection{Relationship between potential and internal energy}

One important property of the hydrostatic approximation is that the
total energy content of a given air column in term of internal energy
on one hand and potential energy on the other are proportional
\begin{equation}
\int_0^\infty \Phi \rho dz =\int_0^\infty RT \rho dz
\end{equation}
This relationship is no more straitforward in the numerical formulation.
However the following finite diference relationship
\begin{equation}
\label{eq:ener1}
\sum_{l=1}^N\Phi_l\dz\sigma_l -\Phi_s = \sum_{l=1}^N RT_l\dz\sigma_l = 
\psk\sum_{l=1}^N \kappa \h_k s_k\dz\sigma_k
\end{equation}
can be imposed by chosing a particular relationship between the location
of the $\sigma$ and $s$ levels as is shown now.

From the discretized hydrostatic equation \eq{hydrostat}, we can write
\begin{equation}
\Phi_l=\Phi_s+\psk \sum_{k=1}^l \depb{\av{\h}{z}}_k \dz s_k
\end{equation}
with the convention at the lower boundary that
\begin{equation}
\depb{\av{\h}{z}}_1=\h_1
\end{equation}
The gravitational energy of a complete column of air is thus given by
\begin{eqnarray}
\sum_{l=1}^N\Phi_l\dz\sigma_l -\Phi_s
&=&\psk\sum_{l=1}^N\sum_{1\le k \le l} \depb{\av{\h}{z}}_k
\dz s_k\dz\sigma_l\\
&=&\psk\sum_{k=1}^N\sum_{k\le l \le N} \depb{\av{\h}{z}}_k
\dz s_k\dz\sigma_l\\
&=&\psk\sum_{k=1}^N\depb{\av{\h}{z}}_k\dz s_k \sigma_k\\
&=&\psk\sum_{k=1}^N\h_k\depb{\av{\dz s_k \sigma_k}{z}}_k
\end{eqnarray}
where
\begin{equation}
\depb{\av{\sigma\dz s}{z}}_1=\frac{\sigma_2 \dz s_2}{2}+\dz s_1
\end{equation}
and
\begin{equation}
\depb{\av{\sigma\dz s}{z}}_N=\frac{\sigma_{N-1} \dz s_{N-1}}{2}
\end{equation}
The \eq{ener1} can hence be satisfied by just imposing that
 $\av{\sigma\dz s}{z}=\kappa s \dz\sigma$.


\subsection{Numerical conservation}

We will now derive the time evolution of the total energy content as
obtained from the finite diference formulation of the primitive
equations.
We game is to show that, for adiabatic flows we have
\begin{equation}
\dt{\sum_\sigma\dz \sigma\avg{\pext\dep{V^2/2+C_v T+\Phi}}}=0
\end{equation}
where $\sum_\sigma$ stands for the summation over all the vertical
column of air and $\avg{}$ accounts for the summation over all the grid
points (the area of the mesh is included in the extensif pressure
$\pext$).

\paragraph{Preliminary remark:}
as usually done in finite difference analysis, the
numerical conservation is only concerned with partial and not temporal
derivative. The time derivative are suppose to be exact.\\

The time derivative of the kinetic energy alone can be computed as
\begin{equation}
\dt{K}=\frac{1}{2}
\dt{\depb{\av{\ucov\ucont}{x}}+\depb{\av{\vcov\vcont}{y}}}
\end{equation}
\begin{equation}
\frac{1}{2}\dt{\av{\ucov\ucont}{x}}=\av{\ucont\dt{\ucov}}{x}
\end{equation}
When integrated over the total planet, the classical symmetry property
of the averaging operator ($\avg{a\av{b}{x}}=\avg{b\av{a}{x}}$) can be applied
and one find that
\begin{eqnarray}
\avg{\pext\dt{K}}
&=&\avg{\av{\pext}{x}\ucont\dt{\ucov}}+\avg{\av{\pext}{y}\vcont\dt{\vcov}}\\
&=&\avg{U\dt{\ucov}}+\avg{V\dt{\vcov}}\label{eq:avK}
\end{eqnarray}

\paragraph{The vorticity term.}
By choosing a different discretization of the vorticity term
\cite{Sado:75,Sado:75a,Sado:80} the term
\begin{equation}
\avg{U\depb{Z V}}-\avg{V\depb{Z U}}
\end{equation}
would cancel exactly. This is not the case in the present formulation
which is based on enstrophy conservation but this term will be
forgotten hereafter.

\paragraph{Other terms.}
The time evolution of the kinetic energy then
reduces to
\begin{eqnarray}
\nonumber
\avg{\pext\dt{K}}
&=& -\avg{
 U\dx \dep{\Phi+K}
+sU\av{\h}{X} \dx \psk
+U\frac{\av{\av{W}{X} \dz \ucov}{Z}}{\av{\pext}{X}\dz \sigma}
}\\
&&-\avg{
V \dy \dep{\Phi+K}
+s V \av{\h}{Y} \dy \psk
+V\frac{\av{\av{W}{Y} \dz \vcov}{Z}}{\av{\pext}{Y}\dz \sigma}
}
\end{eqnarray}
By using the antisymmetry property of the partial differential operator
$\avg{a\dx b}=-\avg{b\dx a}$ (integration by part)
\begin{eqnarray}
\nonumber\dz \sigma \avg{\pext\dt{K}}\label{eq:ener2}
&=& \dz\sigma \avg{\dep{\Phi+K} \dep{\dx  U +\dy V} }\\
&&+s \dz\sigma \avg{\psk \depb{\dx\dep{U\av{\h}{x}}+\dy\dep{V\av{\h}{y}}}}\\
&&-\avg{\av{\ucont\av{W}{x}\dz\ucov}{z}}-\avg{\av{\vcont\av{W}{y}\dz\vcov}{z}}
\end{eqnarray}
which can be transformed using the continuity and energy equation
\begin{eqnarray}
\dz \sigma \avg{\pext\dt{K}}\nonumber
&=& -\dz\sigma \avg{\dep{\Phi+K} \dep{\dt{\pext}+
\der{W}{\sigma} }}\\\nonumber
&&-s \dz\sigma \avg{\psk
\depb{\dt{\dep{\pext\h}} +\frac{\dz\dep{W\av{\h}{z}}}{\dz\sigma} }}\\
&&-\avg{\av{\ucont\av{W}{x}\dz\ucov}{z}}-\avg{\av{\vcont\av{W}{y}\dz\vcov}{z}}
\end{eqnarray}
or
\begin{eqnarray}
\sum_\sigma\dz \sigma \avg{\dt{\dep{\pext K}}}\nonumber
&=& -\dz\sigma \avg{
\Phi \depb{ \dt{\pext}+\frac{\dz W}{\dz\sigma} }
+s \psk
\depb{  \dt{\dep{\pext\h}} +\frac{\dz \dep{W\av{\h}{z}}}{\dz \sigma}  }
}\\
\label{eq:ener3}
&&-\avg{\av{\ucont\av{W}{x}\dz\ucov}{z}+\av{\vcont\av{W}{y}\dz\vcov}{z}
-K\dz W}
\end{eqnarray}
When integrated vertically, the vertical means in the first two
terms of the last line disapear
and it can then be shown easily that the terms of the last line exactly
cancel. 
This is no more the case indeed when using the new numerical formulation
which includes a constraint on the numerical conservation of angular momentum
\cite{Hour:92e}. However, one has to remenmber that the energy is not
exactly conserved in both cases because of the discretization of the
vorticity term.

The hydrostatic relationship,
\begin{eqnarray}\nonumber
\sum_\sigma\Phi\dz W&=&-\sum_\sigma W \dz\Phi\\
&&=\sum_\sigma W\psk  \av{\h}{z} \dz s
\end{eqnarray}
leads to the cancellation of two other terms in the first line of
\eq{ener3}:
\begin{equation}
\sum_\sigma\dz \sigma \avg{
\Phi\dz W +s \psk \frac{\dz \dep{W\av{\h}{z}}}{\dz \sigma}
}=0
\end{equation}

Finally, by applying \eq{ener1}, the evolution of the total  kinetic energy can
be written as
\begin{eqnarray}\nonumber
\sum_\sigma\dz \sigma\avg{\dt{\dep{\pext K}}}
&=&-\avg{\Phi_s\dt{\pext}}\\\nonumber
&&-\sum_\sigma\dz\sigma\avg{\kappa\psk s\h\dt{\pext}+s\psk\dt{\dep{\pext\h}}}
\end{eqnarray}

This term is exactly the oposite of the time derivative of the sum of the total
internal plus potential energy (or total "total potential" energy)
since \eq{ener1} insures that
\begin{equation}
\dt{\sum_\sigma \dz\sigma\pext\dep{C_v T + \Phi}}
=\Phi_s\dt{\ps}+\sum_\sigma s \dz\sigma\dt{\dep{\pext\psk \h}}
\end{equation}
which leads finally to the numerical conservation of the total energy:
\begin{equation}
\dt { \depb{\sum_\sigma\dz \sigma\avg{
K+C_v T + \Phi}}}
=0
\end{equation}
Once again, for the particular numerical formulation under
consideration, this is only true for non rotational flows.

\subsection{Horizontal energy budget}

The main part of those derivations are concerned with vertical
integrations over a column of air. Therefore, the previous analysis can
be used in order to study the horizontal transport of energy by the
horizontal circulation. 
\eq{avK} introduces some numerical error in the energy budget if applied
to a given air column.
The only other part where horizontal averages or derivatives play a role
is for the derivation of \eq{ener2} based on the use of the
anti-symmetry property of the finite-difference derivative. In fact,
this property hides an integration by part:
\begin{equation}
a_1 b_1 - a_0 b_0 = \frac{\dep{a_1 - a_0} \dep{b_1+b_0}}{2}+
\frac{\dep{a_1 + a_0} \dep{b_1-b_0}}{2}
\end{equation}
or
\begin{equation}
\label{eq:parder}
\dx{ab}=\av{a}{X}\dx b+\dx a \av{b}{X}
\end{equation}
where the term on the left hand side disapears when sumed on the
total domain.
When sumed on a limited area, the same relation can be applied
introducing the boundary terms which are fundamentaly fluxes for our
particular application. 
Therefore, for each time the antisymmetry property is used, a boudary
term has to be introudced.
The anti-symmerty property is applied for the term of the type $U\dx{\Phi+K}$
and $s U \av{\h}{X}\dx\psk$ (same thing with $V$ instead of $U$) so that
the time evolution of the total energy of a given air column
can be computed as
\begin{equation}
\dt{\depb{\sum_\sigma\dz\sigma \pext\dep{K+C_v T+\Phi}}}=\\
\sum_\sigma \dz\sigma \depb{
\dx \dep{U\av{\Phi+K+s\h\psk}{X}}+\dy \dep{V\av{\Phi+K+s\h\psk}{Y}}
}
\end{equation}
where the term $U\av{s\h\psk}{X}$ accounts for both the advection of the
internal energy and the work of the pressure forces on the boundaries.

\paragraph{Warning:}
this does not give the exact time evolution of the column of
air as computed by the model since \eq{parder} is not exactly the one which
is used in the derivation of \eq{ener2}.

\bibliography{/d2/hourdin/tex/biblio/fred}

\end{document}
