\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 $kj$ \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 \ \ \ } ij+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\%.