
\def\scald#1{\scal{#1}}
\def\tn{n}
\def\tnp1{{n+1}}
\def\scaln#1{{\scal{#1}}^{\tn}}
\def\scalnp1#1{{\scal{#1}}^{\tnp1}}
\def\transp#1{{#1}^T}
\def\Mn{M^{\tn}}
\def\Mnp1{M^{\tnp1}}
\def\cn{c^{\tn}}
\def\cnp{c^{\tnp1}}
\def\cnim1{c_{i-1}^{\tn}}
\def\cni{c_i^{\tn}}
\def\cnj{c_j^{\tn}}
\def\cnip{c_{i+1}^{\tn}}
\def\cnpi{c_i^{\tnp1}}
\def\cnpip{c_{i+1}^{\tnp1}}
\def\cnpim1{c_{i-1}^{\tnp1}}
\def\cnpj{c_j^{\tnp1}}
\def\cstn{{c^*}^{\tn}}
\def\cstnp1{{c^*}^{\tnp1}}
\def\gradJn{{\nabla}_{\cn}J}
\def\gradJnp1{{\nabla}_{\cnp}J}
\def\E{{\cal E}}
\def\F{{\cal F}}
\def\inv#1{{\dep{#1}}^{-1}}

\section{Symétrie des modèles numériques\label{sec:numeric}}

\subsection{Mise en {\oe}uvre dans LMDZ}

Les développements présentés ci-dessus ont été mis en {\oe}uvre dans LMDZ.

Les intégrations à rebours dans le temps sont évidemment effectuées
uniquement en mode débranché. Il n'est en aucun cas question ici
d'inverser la météorologie elle-même.

\def\Fc{F_{\mbox{conv}}}
\def\Fd{F_{\mbox{diff}}}

On obtient finalement le mode rétro à partir du mode direct débranché dont
les équations complètes s'écrivent
\begin{equation}
\label{eq:direct}
\frac{\partial c}{\partial t}+\V.\grad\ c+\lambda c
+\frac{1}{\rho}\frac{\partial \dep{\Fd+\Fc}}{\partial z}=\sigma
\end{equation}
\begin{equation}\label{eq:conserv}
\dt{\rho}+\div{\rho \V}=0
\end{equation}
\begin{equation}
\Fd=-{K_z\rho\frac{\partial c}{\partial z}}
\end{equation}
\begin{equation}
{\Fd}_{|_{\surf}}=\ssource
\end{equation}
\begin{equation}
{\Fd}_{|_{\infty}}=0 
\end{equation}
\begin{eqnarray}
\Fc&=&
\fu\cu-\fd\cd-(\fu-\fd)c\label{eq:mf3}
\end{eqnarray}
\begin{eqnarray}
\frac{\partial \fu\cu}{\partial z}&=&\eu c -\du \cu\label{eq:up}\label{eq:mf1}\\
-\frac{\partial \fd\cd}{\partial z}&=&\ed c -\dd \cd\label{eq:down}\label{eq:mf2}\\
\frac{\partial \fu}{\partial z}&=&\eu-\du\label{eq:contup}\\
-\frac{\partial \fd}{\partial z}&=&\ed-\dd\label{eq:contdown}
\end{eqnarray}
en relisant à rebours
les archives météorologiques et en remplaçant
$\rho\V$, $K_z$, $\eu$, $\du$, $\ed$, $\dd$
par
$-\rho\V$, $K_z$, $\dd$, $\ed$, $\du$, $\eu$.


\subsection{Illustration numérique}

On montre ici une illustration et un test numérique de la symétrie
du transport atmosphérique en utilisant à nouveau le contexte de la
campagne ETEX-1.
\begin{figure}
\centerline{\includegraphics[height=10.cm]{\local/tr01.eps}
            \includegraphics[height=7.7cm]{\local/legend.eps}
              \includegraphics[height=10.cm]{\local/tr06.eps}}
\caption{Expérience ETEX-1 : cartes de la concentration en surface de
PMCH (en ng/m$^3$) sur une période de 36 heures suivant l'injection
au temps $t_0$.
Le pas de temps pour l'injection et les sorties est de 3 heures
(le panache à $t_0+12$h correspond par exemple en fait à une moyenne entre
$t_0+12$h et  $t_0+15$h).
Le trait discontinu correspond au transect montré sur la figure
suivante.\label{fg:etex}}
\end{figure}

\begin{figure*}
\centerline{\includegraphics[height=5.6cm]{\local/tr01z.eps}
            \includegraphics[height=5.3cm]{\local/legendz.eps}
              \includegraphics[height=5.6cm]{\local/tr06z.eps}}
\caption{
Coupe verticale, dans un repère longitude-pression (hPa) le long du transect
représenté sur la figure précédente (avec les mêmes
unités et conventions).\label{fg:etexz}}
\end{figure*}

La colonne de gauche de la \fig{etex}
montre l'évolution temporelle simulée de la concentration
de surface pour une émission de 
340~kg de PMCH à Monterfil, entre $t_0$ et $t_0+3$~h.
La source est repérée par un cercle sur les figures.
24 heures après l'émission, le panache atteint la station allemande D05,
repérée par des carrés sur les mêmes figures.
La concentration maximum est obtenue à cette station à $t_0$+36h.
Dans la colonne de droite, on montre les résultats d'une simulation
à rebours, pour laquelle les 340~kg de PMCH sont injectés à la 
station D05 de façon régulière entre $t_0+39$~h et $t_0+36$~h.
La réciprocité du transport est illustrée par le fait que la même concentration
est observée dans le carré 36 heures après l'émission pour la simulation
directe (carte en bas à gauche de la \fig{etex} montrant la moyenne de la concentration
directe entre $t_0+36$ et $t_0+39$~h) et dans le cercle au moment de l'injection
réelle pour
la simulation à rebours (panneau en haut à droite montrant la
rétro-concentration moyenne entre $t_0+3$~h et $t_0$).
La \fig{etexz} montre les coupes verticales correspondantes le long d'un
transect reliant la source et la station.
On voit sur la figure que le panache est advecté plus rapidement dans sa
partie supérieure. C'est la combinaison du transport horizontal par
un vent cisaillé et du mélange vertical  par la turbulence qui explique
la dispersion rapide du panache et du rétro-panache dans un plan horizontal.


\begin{figure}
\centerline{
\begin{tabular}{cc}
\includegraphics[width=7.6cm]{\local/FIGURES/D05.eps}
\includegraphics[width=8.0cm]{\local/FIGURES/fmts.eps}
\end{tabular}
}
\caption{Test numérique de la symétrie du transport atmosphérique sur
le cas ETEX-1.\label{fg:retex}}
A gauche :
Evolution temporelle de la concentration de PMCH observée
et simulée  à la station D05 avec les schémas de Van Leer I et de Godunov ainsi que
les résultats obtenus à partir de séries de rétro-simulations
effectuées à partir de la station.

A droite :
FMTs pour les 11 stations retenues pour les études d'intercomparaison
de modèles.
On montre les FMTs mesurant l'écart entre simulation directe et observation
pour le schéma  I de Van Leer et pour le schéma de Godunov, ainsi que les
FMTs mesurant l'écart entre simulation directe et reconstitution inverse
pour ces deux schémas.
%\label{fg:FMTs}}
%\label{fg:station}}
\end{figure}

La partie de gauche de la
\fig{retex} montre la comparaison entre les estimations directe et
rétro de l'évolution temporelle de la concentration de PMCH
pour la station D05 ainsi que les observations.
Ici, le PMCH est injecté de façon uniforme entre
$t_0$ et $t_0+12$h comme dans la réalité et on considère la moyenne
des concentrations sur 3 heures.
De façon symétrique, on considère la moyenne sur les 12 heures d'injection
de rétro-panaches émis toutes les 3 heures.
Avec le schéma I de Van Leer, les simulations directe (carrés noirs)
et rétro (signes {\bf +}) diffèrent, mais la différence est moindre qu'entre
chacune des simulations et les observations.
Les petites différences entre estimations directe et rétro proviennent
en fait pour l'essentiel de la violation de la symétrie du transport
par les limiteurs de pentes introduits dans le schéma de Van Leer pour garantir
positivité et monotonie.
Avec le schéma de Godunov (carrés blancs et étoiles) en effet, les simulations
directe et rétro sont presque confondues à la précision de la figure.

Ce comportement est confirmé par le calcul des FMTs pour les 11
stations privilégiées
pour les analyses ETEX (se reporter à la \sec{validation} pour plus de détails 
sur l'analyse des simulations directes).
La simulation directe avec le schéma I de Van Leer comparé
aux observations (carrés noirs) montre un FMT moyen
de 40$\%$ environ alors que, pour le même schéma,
les FMTs mesurant l'écart entre estimations directe et
rétro pour le schéma I de Van Leer dépassent toujours 75$\%$ (cercles noirs).
La différence entre estimation directe et rétro est bien plus faible
quand on utilise le schéma de Godunov
(ronds blancs, avec un FMT moyen de 98,2$\%$).
Enfin, la symétrie est quasiment exacte quand on inverse l'ordre d'appel
aux différents modules de transport dans le calcul rétro (signes {\bf +},
FMT moyenne de 99,5\%) comme on l'explique plus loin.

On étudie ci-dessous systématiquement la symétrie des algorithmes
utilisés dans LMDZ. Dans le modèle, l'advection, la diffusion turbulente,
la convection et la décroissance radioactive sont appelées de façon
séquentielle. L'intégration du modèle peut donc être vue comme une succession
de pas de transition entre un champ de concentration $\cn$ et
$\cnp$, un pas de temps étant une succession particulière de tels pas.
Nous allons considérer la symétrie de chaque pas individuel en comparant
l'algorithme de rétro-transport obtenu à partir de transformations
simples sur le modèle direct (changement de signe et permutation des flux de
masse pour la convection par exemple) et l'adjoint du code numérique.

\subsection{Modèle adjoint}


Il faut donc, pour discuter en détail ces aspects numériques, introduire
l'adjoint d'un code numérique.
La définition mathématique de l'opérateur adjoint est la suivante.
Soient $\E$ et $\F$ deux espaces de Hilbert munis de deux produits scalaires
$\scal{\ ,\ }_\E$ et $\scal{\ ,\ }_\F$, et soit
$L$ un opérateur de $\E$ dans $\F$, l'adjoint 
$L^*$ de $L$ est un opérateur de $\F$ dans $\E$ défini par la relation 
\begin{equation}
\scal{L\phi,\psi}_\F=\scal{\phi,L^*\psi}_\E
\end{equation}
pour tout vecteur $\phi$ de $\E$ et tout vecteur $\psi$ de $\F$.

Etant donné un produit scalaire discrétisé au pas $n$,
$\scaln{\phi,\psi}= \transp{\phi}\Mn\psi$,
où $\Mn$ est une matrice symétrique définie positive, l'adjoint
du modèle numérique direct permet de calculer
l'évolution à rebours dans le temps
du gradient  ${\nabla}_{c}J$ d'une fonction objective quelconque $J$
par rapport aux variables d'état du modèle
\cite[cf. e.~g.][]{Tala:87}.
Ceci se montre simplement par comparaison des expressions du gradient
aux pas 
$\tn$ et $\tnp1$
\begin{eqnarray}
dJ &=& \transp{\dep{\gradJn}}\Mn d\cn \label{eq:graddef1}\\
   &=& \transp{\dep{\gradJnp1}}\Mnp1 d\cnp\label{eq:graddef2}
\end{eqnarray}
\def\inv#1{{\dep{#1}}^{-1}}
Après introduction du modèle linéaire tangent $L$ (le modèle lui-même
pour les algorithmes linéaires), qui décrit l'évolution
directe d'une perturbation $dc$ de $c$
\begin{equation}\label{eq:tmp2}
d\cnp=L d\cn
\end{equation}
et en prenant les transposées, on obtient
$\gradJn=L^*\gradJnp1$, où
\begin{equation}\label{eq:recmatr}
L^*=\inv{\Mn}\transp{L}\Mnp1
\end{equation}
est bien l'adjoint de $L$ par rapport aux produits scalaires $\Mn$ et
$\Mnp1$.


Pour le produit scalaire canonique, $M=I$, on trouve le résultat classique~:
 $L^*=\transp{L}$.

 Le produit scalaire pondéré par la masse de l'air correspond à
$M=\mbox{diag}(m_i)$ où
$m_i$ est la masse d'air dans la maille $i$.
Si on prend comme convention d'indiçage
$\psi_i=\sum_j L_{i,j} \phi_j$,
on trouve, pour les éléments des matrices directe et adjointe, la relation
\begin{equation}\label{eq:reccomp}
m_i^n L^*_{i,j} = m_j^{n+1} L_{j,i}
\end{equation}

\def\MES{\cout}
Si le modèle de transport est linéaire, pour tout pas de temps $n$ entre la
source et la détection, une mesure linéaire de $c$ peut être évaluée
comme 
$\MES=\scaln{{\nabla}_{\cn}\MES ,\cn}$, qui, pour
${c^*}^n=\gradJn$, est un équivalent numérique de l'\eq{dtmesure}.

\subsection{Vérification de la symétrie temporelle des algorithmes}

La symétrie du code numérique est donc équivalente à
l'identité entre le mode rétro-transport $R$ du modèle direct $L$ (obtenu
pour LMDZ par la transformation
$(\V,\eu,\du,\ed,\dd)\rightarrow (-\V,\dd,\ed,\du,\eu)$)
et l'opérateur adjoint du modèle numérique direct.
La symétrie sera donc testée  au moyen de la relation
(\ref{eq:recmatr})
(ou (\ref{eq:reccomp})) avec $L^*=R$.


Pour des paramétrisations linéaires qui ne modifient pas la densité
de l'air -- c'est le cas par exemple de la diffusion
turbulente, de la convection, ou de l'advection par un champ de vent non divergent --
c'est à dire dans le cas où
$\Mn=\Mnp1=M$,
il est plus commode de tester une version légèrement modifiée des relations
(\ref{eq:recmatr}) ou (\ref{eq:reccomp}) en écrivant le modèle
sous une forme intégrale comme
\begin{equation}\label{eq:tmp0}
M\dep{\cnp-\cn}=A\cn
\end{equation}
où $L=I+M^{-1}A$. D'après \eqb{recmatr},
 $L^*=I+M^{-1}A^T$ de sorte que le modèle adjoint s'écrit
 \begin{equation}
 M\dep{\cstnp1-\cstn}=A^T\cstn
 \end{equation}
 Pour un modèle de la forme (\ref{eq:tmp0}),
 la symétrie sera donc assurée dans le monde numérique si la matrice 
 $B$ associée au mode rétro-transport est égale à la transposée de la
 matrice du modèle direct, 
 $B=A^T$.


 La même relation $B=A^T$ garantit en fait plus généralement
 la symétrie du schéma
\begin{equation}\label{eq:tmp4}
M\dep{\cnp-\cn}=A \depb{\gamma\cn+\dep{1-\gamma}\cnp}
\end{equation}
qui couvre les cas des schémas temporels 
explicite ($\gamma=1$), implicite ($\gamma=0$)
et de Crank-Nicholson ($\gamma=1/2$).
Après avoir transformé le schéma généralisé en
\begin{equation}\label{eq:tmp3}
\cnp={\depb{I-\dep{1-\gamma}M^{-1}A}}^{-1}\depb{I+\gamma M^{-1}A}\cn
\end{equation}
l'application de l'\eq{recmatr} conduit en effet à
\begin{eqnarray}
\cstn
&=& M^{-1}\depb{I+\gamma A^TM^{-1}}{\depb{I-\dep{1-\gamma}A^TM^{-1}}}^{-1}M\ \cstnp1 \\
&=& \depb{I+\gamma M^{-1}A^T} {\depb{I-\dep{1-\gamma}M^{-1}A^T}}^{-1} \cstnp1
\end{eqnarray}
Or, les matrices $[I+aQ]$ et ${[I+bQ]}^{-1}$ commutent pour n'importe quelle
matrice $Q$ et scalaires $a$ et $b$. On obtient ainsi
\begin{eqnarray}
\cstn
&=& {\depb{I-\dep{1-\gamma}M^{-1}A^T}}^{-1} \depb{I+\gamma M^{-1}A^T}
\cstnp1
\end{eqnarray}
qui est de la forme (\ref{eq:tmp3}), $A$ étant remplacée par $A^T$.



\subsection{Symétrie du schéma d'advection}

\def\dt{\delta_t}

Les courbes relatives au schéma de Godunov sur la \fig{retex} suggèrent que
le schéma de Godunov est symétrique, ce que nous démontrons ci-dessous.

Considérons dans un premier temps le cas de l'advection unidimensionnelle
avec un champ 
de vent non divergent $u>0$, une grille régulière de pas $\delta x$
et un pas de temps $\delta t$.
Si on note $\alpha=u\delta t/\delta x$ le nombre de Courant,
l'évolution de la concentration du traceur pour la maille $i$ est
donnée dans le schéma de Godunov par
\begin{equation}
c_i^{n+1}-c_i^n=\alpha \dep{c_{i-1}^n-c_i^n}
\end{equation}
qui est de la forme (\ref{eq:tmp0}) avec $M=I$.
Quand on l'applique au rétro-traceur $c^*$ avec le vent $-u$, le schéma amont
s'écrit
\begin{equation}
{c^*}_i^{n}-{c^*}_i^{n+1}
=\alpha \dep{{c^*}_{i+1}^{n+1}-{c^*}_i^{n+1}}
\end{equation}
La matrice du second schéma est bien la transposée du schéma direct
ce qui établit l'équivalence entre modèle de rétro-transport et modèle
adjoint.

\def\mnp{{m^{n+1}}}
\def\mn{{m^{n}}}
\def\mnpi{{{m_i}^{n+1}}}
\def\mni{{{m_i}^{n}}}
% \def\mnpi{{{m_{i+1}}^{n}}}
\def\mnpip{{{m_{i+1}}^{n+1}}}

La symétrie
du transport est également vérifiée pour un champ de vent non divergent.
Si on note $U_{i+\dem}$ le transfert de masse entre les mailles
$i$ et $i+1$ et les instants $t$ et $t+\delta t$ 
et si on suppose -- pour fixer les idées -- que ce transfert est positif,
le schéma amont s'écrit
\begin{equation}\label{eq:up1}
\mnpi\cnpi-\mni \cni=U_{i-\dem}\cnim1-U_{i+\dem}\cni
\end{equation}
avec
\begin{equation}\label{eq:up2}
\mnpi-\mni =U_{i-\dem}-U_{i+\dem}
\end{equation}
L'\eq{up1} est de la forme (\ref{eq:tmp2}), $L$ étant donné par
\begin{equation}
L_{i,i}=\frac{\mni-U_{i+\dem}}{\mnpi}\mbox{  et  }
L_{i,i-1}=\frac{U_{i-\dem}}{\mnpi}
\end{equation}  
Pour le rétro-transport, le signe de $U$ et les indices doivent être
inversés. Le rétro-transport est associé à une matrice $R$
donnée par
\begin{equation}
R_{i,i}=\frac{\mnpi-U_{i-\dem}}{\mni}=\frac{\mni-U_{i+\dem}}{\mni}
=\frac{\mnpi}{\mni}L_{i,i}
\end{equation}
et
\begin{equation}
R_{i,i+1}=\frac{U_{i+\dem}}{\mni}=\frac{\mnpip}{\mni}L_{i+1,i}
\end{equation}
de sorte que, suivant (\ref{eq:reccomp}), $R=L^*$.

Comme on l'a expliqué plus haut (\sec{volumesfinis}),
le passage en dimension 3 est effectué dans le modèle au moyen
d'un calcul successif dans les 3 directions, successivement en
$x$, $y$, $z$, $y$ et $x$ avec un pas de temps deux fois plus petit dans
les deux directions horizontales.
Dans chaque direction, 
l'\eq{up1} pour le traceur et 
l'\eq{up2} pour l'air sont intégrées simultanément.
Ceci assure la symétrie du schéma tridimensionnel.


On voit donc une fois de plus que le schéma amont de Godunov présente
beaucoup de bonnes propriétés (positivité, monotonie, diminution de la 
variation totale, linéarité, et maintenant
symétrie) mais au prix d'une diffusion numérique importante.
On a déjà vu dans la \sec{Godunov} que le schéma de Godunov peut se récrire
\begin{equation}\label{eq:godu2}
c_i^{n+1}-c_i^n=
\alpha \frac{c_{i-1}^n-c_{i+1}^n}{2}
+\alpha \frac{c_{i-1}^n-2c_i^n+c_{i+1}^n}{2}
\end{equation}
comme la somme d'un schéma d'advection centré du second ordre et d'un 
terme de diffusion numérique avec une diffusivité
$\alpha(\delta x)^2/(2\delta t)=u\delta x/2$.
En appliquant la même transformation au rétro-transport, on obtient
\begin{equation}\label{eq:godur}
c_i^{*\ n}-c_i^{*\ n+1}
=
-\alpha \frac{{c^*}_{i-1}^{n+1}-{c^*}_{i+1}^{n+1}}{2}
+\alpha \frac{{c^*}_{i-1}^{n+1}-2{c^*}_i^{n+1}+{c^*}_{i+1}^{n+1}}{2}
\end{equation}
qui est à nouveau la somme d'un schéma d'advection centré et d'un terme
de diffusion. Comme pour la diffusion turbulente, la diffusion numérique
est la même dans les schémas direct et rétro.
A noter également que le schéma d'ordre 2 linéaire obtenu en ignorant le
dernier terme des Eqs~\eqb{godu2} et \eqb{godur}, est également symétrique.

Les schémas en volumes finis plus sophistiqués, comme le schéma I
de Van Leer, peuvent souvent aussi être décrits comme la somme d'un schéma
centré et d'un terme de  diffusion, introduit pour éviter les oscillations
numériques ou la dispersion. Cependant, ce terme de diffusion dépend alors
en général de la concentration du traceur (elle est plus active là où le
traceur montre de brusques variations de concentration), rompant la linéarité
et donc la symétrie du schéma.

Le cas particulier
du schéma de Van Leer et les implications pour l'utilisation pour
l'inversion du transport sont discutées plus loin.


\subsection{Diffusion turbulente}

Pour la diffusion turbulente, on utilise dans le modèle LMDZ
un schéma implicite en temps et centré sur la verticale~:
\def\kzd{\tilde{K}}
\begin{equation}\label{eq:Kdiffd}
m_i\dep{\cnpi-\cni}=\kzd_{i+\dem}\dep{\cnpip-\cnpi}-\kzd_{i-\dem}\dep{\cnpi-\cnpim1}
\end{equation}
($m$ n'est pas affecté par le mélange)
où $i$ est à présent l'indice vertical et où $\kzd_{i+\dem}$ est une estimation de
$K_z \rho\delta_x\delta_y\delta_t/\delta_z$ à
l'interface entre les couches $i$ et $i+1$.
L'\eq{Kdiffd} est de la forme (\ref{eq:tmp4}), $A$ étant donnée par
\begin{eqnarray}
A_{i,i-1}&=&\kzd_{i-\dem}\\
A_{i,i}  &=&-\kzd_{i+\dem}-\kzd_{i-\dem}\\
A_{i,j}&=&0 \mbox{\ \ \ pour\ \ \ }|i-j|>1
\end{eqnarray}
Nous avons montré plus haut sur des bases physiques que la diffusion rétro
était identique à la diffusion directe.
On vérifie aussi que $A=A^T$, ce qui
montre que la forme adjointe et la forme rétro-transport 
du schéma (\ref{eq:Kdiffd}) sont les mêmes.


\subsection{Modèles de convection en flux de masse}

Les schémas en flux de masse utilisés dans le modèle sont
également symétriques, mais c'est un peu plus fastidieux à démontrer.
Comme précédemment, on va se restreindre
au cas
d'une ascendance compensée par une subsidence dans l'environnement.

Comme on l'a déjà expliqué au \ch{lmdzt},
les équations de continuité pour l'air et le traceur dans l'ascendance
(resp. Eqs. \ref{eq:contup}
et \ref{eq:mf1}) sont discrétisées comme
\begin{equation}\label{eq:contupd}
E_i+F_{i-\dem}=D_i +F_{i+\dem}
\end{equation}
et
\begin{equation}
E_i\cn_i+F_{i-\dem}\cu_{i-1}=\cu_i\dep{D_i+F_{i+\dem}}
\end{equation}
où $E_i\simeq \eu\delta z \delta t$ et $D_i\simeq \du\delta z \delta t$
sont les entraînement et détraînement vers et depuis l'ascendance pour
la couche $i$ durant le pas de temps $\delta t$,
et $F_{i+\dem}\simeq \fu\delta t$ est le transfert de masse entre les couches
$i$ et $i+1$,
avec les conditions supplémentaire que $F$ est nul en haut et en bas
de la colonne convective.

L'évolution temporelle de la concentration grande échelle de traceur
restreinte au transport convectif,  obtenue en combinant
les Eqs.~\ref{eq:direct} et \ref{eq:mf3} avec $\V=0$, $\lambda=0$, $\Fd=0$ et
$\fd=0$, qui s'écrit
\begin{equation}
\rho\frac{\partial c}{\partial t}=
-\frac{\partial \fu\dep{\cu-c}}{\partial z}
\end{equation}
est discrétisée sous la forme
\begin{eqnarray}
m_i \cnp_i - m_i \cn_i &=& F_{i-\dem} \cu_{i-1}-F_{i+\dem}\cu_{i}
                      +F_{i+\dem}\cnip-F_{i-\dem}\cni\\
&=& D_i \cu_i -E_i \cn_i
                      +F_{i+\dem}\cnip-F_{i-\dem}\cni
\end{eqnarray}

Ce schéma est de la forme (\ref{eq:tmp0}).
Puisque le modèle est linéaire, les éléments de la matrice $A$ peuvent
être obtenus en calculant
$m_i\cnp_i-m_i\cn_i$ dans la couche  $i$ pour un traceur injecté dans
la couche $j$ ($c_k^n=\delta_{k,j}$).

Si on calcule d'abord les concentrations dans l'ascendance, pour cette
injection particulière, on trouve, pour $k<j$, ${\cu}_k=0$.
Pour la couche d'injection, la concentration dans l'ascendance est donnée par
\begin{equation}
\cu_j\dep{D_j+F_{j+\dem}}=E_j\cnj
\end{equation}

Pour $k>j$
\begin{equation}
\cu_k\dep{D_k+F_{k+\dem}}=\cu_{k-1}F_{k-\dem}
\end{equation}
et, après itération,
\begin{equation}
\cu_i
=\prod_{k=j+1}^i \frac{F_{k-\dem}}{D_k+F_{k+\dem}}\cu_j
=\frac{\prod_{k=j+1}^i F_{k-\dem}}{\prod_{k=j}^{i} \dep{D_k+F_{k+\dem}}}E_j\cnj
\end{equation}
Le traceur est finalement détraîné dans l'environnement, donnant pour un
détraînement dans la couche $i$ avec $i>j$, 
\begin{equation}
A_{i,j}=D_i\cu_i=
\frac{\prod_{k=j+1}^i F_{k-\dem}}{\prod_{k=j}^{i} \dep{D_k+F_{k+\dem}}}
E_j D_i
\end{equation}
Dans la couche d'injection $i=j$, un calcul similaire donne
\begin{equation}
A_{j,j}=-F_{j-\dem}-E_j+\frac{D_j E_j}{\dep{D_jF_{j+\dem}}}
\end{equation}
La concentration grande échelle est également modifiée par la
subsidence juste en dessous de la couche d'injection
($m_{j-1}\cnp_{j-1}=F_{j-\dem}$)
d'où l'on tire
\begin{equation}
A_{j-1,j}=F_{j-\dem} \mbox{\ \ \  et \ \ \ }
A_{i,j}=0 \mbox{\ \ \ pour \ \ \ } i<j-1 \mbox{\ .}
\end{equation}

Le schéma pour le rétro-transport est obtenu en inversant le rôle de
$D$  et $E$ et en changeant le sens vertical de la propagation des
indices.
Le schéma rétro est donc également  de la forme
(\ref{eq:tmp0}), la matrice  $A$ étant remplacée par la matrice $B$ dont les
éléments sont obtenus en considérant
une rétro-injection dans la couche $j$, ce qui donne
\begin{eqnarray}
B_{i,j}&=&
\frac{\prod_{k=i}^{j-1} F_{k+\dem}}{\prod_{k=i}^{j} \dep{E_k+F_{k-\dem}}}
D_j E_i \mbox{\ \ \ pour \ \ \ }i<j\\
B_{j,j}&=& -F_{j+\dem}-D_j+\frac{E_j D_j}{\dep{E_j+F_{j-\dem}}}\\
B_{j+1,j}&=&F_{j+\dem}\mbox{\ \ \ pour \ \ \ et \ \ \ }
B_{i,j}=0 \mbox{\ \ \ pour \ \ \ } i>j+1 \mbox{\ .}
\end{eqnarray}

En utilisant l'équation de continuité
(\ref{eq:contupd}), on vérifie que
$B=A^T$ ce qui montre que les transports convectifs
adjoint et rétro sont identiques.


\subsection{Puits linéaires}

Les puits linéaires conservent aussi la masse d'air et peuvent être 
discrétisés au moyen de la matrice diagonale $A=\mbox{diag}(m_i \lambda_i)$
qui assure évidemment la symétrie.

\subsection{Symétrie du modèle complet}

Dans le modèle direct, l'advection de grande échelle, la diffusion 
turbulente et le transport convectif sont traités de façon séquentielle.
Si on veut vraiment obtenir la symétrie numérique avec le schéma de
Godunov, il faut inverser complètement l'ordre de cette séquence dans
le modèle. On ne le fait pas en standard dans le mode rétro de LMDZ
qui utilise de toutes façons le schéma I de Van Leer, non symétrique.

C'est la non inversion de
cet ordre des opérateurs qui explique la différence résiduelle
entre calcul direct et rétro dans le cas où on utilise le schéma
de Godunov (cercles blancs sur la partie de droite de la \fig{retex}).
 Comme on l'a déjà dit, si on inverse cette séquence dans le calcul
 rétro
 (signes ${\bf +}$), le FMT moyen pour les 11 stations passe
 de 98,2 à 99,5\%.

 



