%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Al1 %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4paper,11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
% amsfonts permet le \mathbb(C) attention, peu de lettres disponibles
\usepackage{amsfonts}
%%% xr est un .sty pour les references externes externaldocument
\usepackage{xr,eepic,psfig,pstricks,multirow}
\usepackage{floatflt}
\usepackage{fancyheadings}
\usepackage{footnpag}
\usepackage{rotating}
\usepackage{euscript}
\usepackage{tabularx,xspace}
%\usepackage{minitoc}
\usepackage{makeidx}
\usepackage[frenchb]{babel}
\usepackage[latin1]{inputenc}

\topmargin -1cm
\addtolength{\textheight}{4cm}
\tolerance=8000
\pretolerance=4000
\textwidth 15cm
\oddsidemargin 0.5cm

%
\newcommand{\cad}{c'est-\`a-dire\xspace}
\newcommand{\Minik}{Mini$\_$ker\xspace}
\newcommand{\bra}[1]{<{#1}\mid}
\newcommand{\ket}[1]{\mid{#1}>}
\newcommand{\mat}[1]{\underline{\underline{\mathcal{#1}}}}
\newcommand{\diag}[1]{{#1}^{d}}
\newcommand{\mean}[1]{\overline{{#1}}}
\newcommand{\ext}[2]{{\mid\, #1 >< #2\mid}}
\newcommand{\scal}[2]{{<\, #1 \, #2>}}
\newcommand{\scalm}[3]{{<\, #1 {#2}^{\dagger} \, #2 #3>}}
\newcommand{\qrayl}[3]{\frac{\scalm{#1}{#2}{#3}}{\scal{#1}{#3}}}
\newcommand{\trac}[1]{{\mathcal{T}r (#1)}}

\externaldocument[P-]{/home/al1/Cahiers/Gains/Propagateurs}
\externaldocument[A-]{/home/al1/Cahiers/Gains/Adjoint}
\externaldocument[SLT-]{/home/al1/Cahiers/Gains/SLT_Al1}


\makeindex
\begin{document}

\def\todec{{\overline {\overline {\Upsilon  }}}_{,dec} }
\def\tox{{\overline {\overline {\Upsilon  }}}_{,x} }
\def\dt{{\delta t}}
\def\edtaux{{exp (-{\delta t  \over \taux } ) }}
\def\deta{{\delta \eta } }
\def\dpi{{\delta \varpi } }
\def\Deta{{\Delta \eta } }
\def\Dt{{\delta t } }
\def\dfi{{\delta \varphi } }
\def\dmu{{\delta \mu } }
\def\Dfi{{\Delta \varphi } }
\def\dts{{\delta t\over 2 } }
\def\dtaf{{\overline {\overline {\partial _\eta f }}} }
\def\d(#1)G{{\overline {\overline {\partial_#1 G}}} }
\def\dA{{\overline {\overline {\partial A}}} }
\def\dtA{\partial_t A}
\def\VD{{\mathcal V \mathcal D}}


\title{Cahiers des algorithmes de Mini-ker\author{Al1, Août 2006}\date{Versions 1.02.xx}}
\maketitle 

\bigskip
{\it Ce manuel est une description de l'environnement
développé pour l'utilisation pratique de \Minik dans ses versions 1.02\footnote{à quelques
détails près, il suffit aussi à comprendre les principes des versions précédentes,
qui ne bénéficient pas de ce type de manuel.}. Il s'agit de
notes permettant de comprendre le code Mortran de \Minik\footnote{il est ainsi conseillé
d'avoir sous la main un ../cfs/exo$\_$o.tmp pour lire ces notes.}, et de pouvoir
le développer dans le même esprit.

\medskip
Le lecteur est renvoyé aux
différents documents concernant les bases théoriques du TEF\index{TEF}, de l'adjoint et des sensibilités,
cités dans les pages {\bf ZOOM}\index{ZOOM}\footnote{{\it http://www.lmd.jussieu.fr/Zoom}},
ainsi que ce qui concerne le filtre de Kalman.

\medskip
Le découpage des séquences de calcul du programme {\bf principal} de \Minik suit une logique
de présentation en même temps que l'ordre des instructions, excluant cependant des options
supplémentaires. A la lecture, il faut ainsi parcourir
plusieurs fois un exemple de {\bf exo$\_$o.tmp} en sautant des passages du code pour suivre
les explications du Cahier.}

\medskip
Enfin, la rédaction de ce cahier se faisant en parallèle avec la programmation de la version {\bf 102},
on notera le passage progressif de routiines de la {\bf mathlib} ancienne ou nouvelle convention
(ex. matmul -> omatmlt) LaPack.
\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Principes et structure du \Minik}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bigskip
Introduction générale aux choix de programmation.
\begin{itemize}
\item définition de langages dédiés dès qu'ils se dégagent
\item principe d'unicité (ou de non-duplication) des calculs;
\end{itemize}
%
L'écriture de langages pour tenter une remontée de la couche algorithmique est une
méthode systématique de la collaboration ZOOM. L'utilisation du précompilateur Mortran\footnote{SLAC, CERN
et RAMSÈS.}
offre les facilités nécessaires à l'écriture symbolique des descriptions de modèles
et des diverses options de calcul. Remarquons cependant que par rapport à ZOOM, cette
systématique en est restée à un niveau modeste dans \Minik, du fait de la simplicité
des algorithmes et de la gestion des matrices, qui sont rangées en un seul bloc par
matrice\footnote{un seul niveau d'emboitement est envisagé pour la version 200.}.

Le code Fortran généré utilise des séquences de calcul répétées dans des contextes différents,
pour le calcul de la trajectoire, des sensibilités et de l'adjoint. On a donc découpé ce code
final en séquences (des paquets d'instructions Mortran ou Fortran); ces séquences sont gérées
par CMZ\footnote{CERN Code Management using Zebra.}
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Résumé des calculs}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%====================================================
\subsection{Calculs de trajectoire}
%====================================================
Le principe de \Minik est celui du TEF-ZOOM\index{TEF}\index{ZOOM} en ne conservant qu'une cellule et
qu'un transfert (avec, ainsi, suppression de l'emboîtement des familles\footnote{du moins jusqu'à
la version 2.00.}).

Un système est ainsi défini par
%
\begin{eqnarray}
%\left{
\partial_t \eta (t) &=&  g(\eta(t),\varphi(t),t) \nonumber\\
\varphi(t) &=& f(\eta(t),\varphi(t),t)
%\right.
\label{tefbas}
\end{eqnarray}
%
duquel on dérive le système d'avance le long de l'espace tangent local :
%
\begin{equation}
\begin{bmatrix} \partial_t \deta \\ 0\end{bmatrix} \quad  = \quad
\begin{bmatrix} H  & B_b \\ C^{\dagger} & D\end{bmatrix} \quad 
\begin{bmatrix} \deta  \\ \dfi \end{bmatrix}\quad + \quad\begin{bmatrix} \Gamma  \\ \Omega \end{bmatrix}
\label{sltc1} 
\end{equation}
%
dans ce système, les $\delta x(\tau), \tau\in [t,t+\dt]$ sont des vecteurs évoluant dans le tube osculant
la trajectoire dans l'espace tangent local, valable de $t$ à $t+\dt$.
Les matrices sont les Jacobiennes\index{Jacobienne} de~(\ref{tefbas}) :
%
\begin{eqnarray}
A = \partial_{\eta} g &;& B_b = \partial_{\varphi} g; \nonumber\\
C^{\dagger} = \partial_{\eta} f &;& D = \partial_{\varphi}f 
\label{tefmat}
\end{eqnarray}
%
Le vecteur $\Gamma=\dt(g((t)) + \frac{\dt}{2}\partial_t g$, et $\Omega$ est en général
nul sauf cas particulier (source, transfert intégré) et est ici conservé par souci de
généralité de la forme de ce système -- cf calcul des sensibilités.

Après application du schéma de discrétisation au deuxième ordre en temps du TEF,
on aboutit au système algébrique suivant
%
\begin{equation}
\begin{bmatrix} A  & B \\ -C^{\dagger} & I-D\end{bmatrix} \quad 
\begin{bmatrix} \Deta  \\ \Dfi \end{bmatrix} \quad = \quad
\begin{bmatrix} \Gamma\\ \Omega \end{bmatrix}
\label{mini1} 
\end{equation}
%
avec les nouvelles matrices dépendant de $\dt$ et des Jacobiennes~: $A=I-\frac{\dt}{2}H$
et $B=-\frac{\dt}{2}B_b$.

La résolution suit les principes du TEF en appliquant une première étape d'élimination
des variables $\deta$, qui donne le système dit ``de couplage'' (entre composantes du transfert):
%
\begin{equation}
[\; I-D +C^{\dagger}A^{-1}B\;]\;\dfi = \Omega + C^{\dagger}A^{-1}\Gamma
\label{couplage} 
\end{equation}
%
C'est la matrice $\mathcal{C}=[\; I-D +C^{\dagger}A^{-1}B\;]$ qui est dite `` de couplage'', et
le vecteur $\dfi_{ins}= \Omega + C^{\dagger}A^{-1}\Gamma$ représente l'évolution
insensible du transfert (unique) de \Minik.

Après résolution du système linéaire~(\ref{couplage}), on a $\dfi$, et il reste à
calculer l'avance d'état~: $\deta = \Gamma + A^{-1}B\dfi$, ce qui finit de déterminer
l'avance du système en $dt$. L'état suivant~$\eta(t+\dt)=\eta(t)+\deta$ étant connu, on
détermine le nouveau vecteur des transferts par la résolution du système implicite
$\varphi(t) = f(\eta(t),\varphi(t),t)$, qui doit donner une solution peu éloignée de celle-ci
%
\begin{equation*}
\varphi(t+\dt)= \varphi(t)+\dfi\quad\sim\textit{ solution de } \varphi'= f(\eta(t+\dt),\varphi')
\end{equation*}
%
%====================================================
\subsection{Calculs de sensibilité}
%====================================================
Le système déterminant les sensibilités\index{sensibilité} aux conditions initiales est
%
\begin{equation}
\begin{bmatrix} \partial_t \Phi\\ \Psi \end{bmatrix} \quad = \quad
\begin{bmatrix}  H  & B_b \\ C^{\dagger} & D\end{bmatrix} \quad 
\begin{bmatrix} \Phi  \\ \Psi \end{bmatrix} (t,\tau) 
\label{sltsens} 
\end{equation}
%
et s'intègre ainsi comme le système d'avance~(\ref{sltc1}) -- avec la même matrice de couplage.
Les inconnues sont à présent les matrices $\Phi$ et$\Psi$, qui étendent alors les vecteurs
$\deta,\dfi$, avance $\Phi(0,0)=I$ matrice unité.

On profite de cette propriété pour résoudre ce système en même temps
que l'avance de trajectoire. Les différences ne portent que sur les termes sources,
\cad les vecteurs étendant $\Gamma, \Omega$.

L'application dans \Minik est alors très simple~:
%
\begin{itemize}
\item les jacobiennes\index{Jacobienne} sont calculées par dérivation formelle des modèles formulés
au début du code~({\bf principal});
\item  on passe des Jacobiennes au calcul des matrices de l'avance en prenant
en compte l'incrément temporel;
\item l'appel au sous-programme {\bf ker} effectue l'élimination des cellules~(\ref{couplage}),
\l'appel à LaPack (SGESV) pour la résolution de ce dernier système linéaire, mais complété
par un certain nombre choisi de colonnes de $\Phi,\Psi$.
\item L'avance de $\eta,\Phi$ est alors effectuée dans {\bf principal} au retour de {\bf ker};
\item il ne reste plus qu'à résoudre l'avance implicite de $\varphi$.
\item Quand ceci est fait, la matrice $D$ est
déterminée, $C^{\dagger}$ peut être calculée, ce qui suffit à déterminer $\Psi$.
\end{itemize}
%
En résumé, on intègre dans le même mouvement un ensemble de systèmes possédant les mêmes
Jacobiennes, seule la trajectoire nécessitant la résolution d'un système implicite.
Les différents systèmes ne diffèrent que par leurs vecteurs connus et leurs C.I.

En particulier, le système des sensibilités (cf cahier Effluents) mène à une variété
de sources et donc, vecteurs de type $\Gamma,\Omega$, comme nous l'expliciterons en temps utiles.
%
%====================================================
\subsection{Calculs dans la version 1.02}
%====================================================
La version \Minik 102 intègre en avant, en plus de la trajectoire et des sensibilités
aux variables du TEF\index{TEF}~:
%
\begin{itemize}
\item un calcul de l'avance d'un type de transfert nouveau, dit ``d'observables''
(ou ``probes'')~: vecteur $\mu$;
\item le calcul des sensibilités des variables du TEF à une liste de paramètres -- définis
par leurs symboles Fortran dans la formulation des modèles (formant un vecteur $\pi$ ordonné
par la liste);
\end{itemize}
%
Ce système a pour équation~:
%
\begin{equation}
\begin{bmatrix} \partial_t \dpi\\\partial_t \deta \\ \dfi \\ \dmu\end{bmatrix} \quad  = \quad
\begin{bmatrix} 0 & 0 & 0 & 0 \\ H_{\varpi} & H  & B_b &0 \\
                C^{\dagger}_{\varpi}&C^{\dagger} & D & 0 \\
                c_1 & c_2 & d_{\varphi\mu} & 0\end{bmatrix} \quad 
\begin{bmatrix} d\varpi \\ \deta  \\ \dfi \\ \dmu \end{bmatrix}\quad +
  \quad\begin{bmatrix} 0\\ \textit{termes}  \\ \textit{sources} \\0\end{bmatrix}
\label{sltc2} 
\end{equation}
%
dans lequel on peut faire les remarques suivantes~:
%
\begin{itemize}
\item les paramètres sont des cellules -- puisqu'il faut définir des C.I., d'équation
$\partial_t\varpi=0$;
\item les variables sont insensibles aux obs par définition; attention, dans le cas du filtrage
optimal de Kalman, il pourra y avoir bouclage par une matrice de gain $K$;
\item on a {\it a priori} exclu tout couplage instantané (par une matrice $D_{\mu}$) entre les
obs, ni envisagé de sources; ceci peut évidemment être changé;
\item le système a la même matrice de couplage que~(\ref{tefbas});
\end{itemize}
%
La conséquence sur le calcul reste standard : on a d'une part à faire avancer la trajectoire,
de l'autre à calculer des vecteurs connus supplémentaires ainsi que des systèmes linéaires
explicites supplémentaires. Pour la raison que le compilateur Fortran n'accepte pas les
dimensions de tableau nulles, on étend la dimension des matrices $H,C^{\dagger}$ et $D$,
l'objectif restant de ne pas stocker de blocs-matrices nuls.

Nous allons détailler les calculs à effectuer. Le système algébrique à intégrer
est à présent~:
%
\begin{equation}
\begin{bmatrix}    A  & B &0 \\
                C^{\dagger} & I-D & 0 \\
                 -c_2 & -d_{\varphi\mu} & I\end{bmatrix} \quad 
\begin{bmatrix} \deta  \\ \dfi \\ \dmu \end{bmatrix}\quad =
  \quad\begin{bmatrix} \Gamma  \\ \Omega \\0\end{bmatrix}
\label{mini2} 
\end{equation}
%
%-------------------------------------------------
\subsubsection{Calculs de la trajectoire}
%------------------------------------------------
%
Standards, puisque $\varpi$ n'avance pas et que les variables standards sont insensibles
aux obs. On avance les obs par~:
%
\begin{equation}
\Psi_{\mu} = c_2\deta + d_{\varphi\mu}\dfi
\label{sltc3} 
\end{equation}
%
et ceci dès que la détermination implicite de $\varphi$ a été effectuée et après calcul
des jacobiennes.
%-------------------------------------------------
\subsubsection{Calculs des sensibilités}
%------------------------------------------------
%
En remarquant que la sensibilité de $\varpi$ à $\varpi$ reste égale à $I$, écrivons-en le système complet des
sensibilités correspondant au SLTC (\ref{sltc2}) avant de le commenter~:
%
\begin{equation}
\begin{bmatrix} \partial_t [\,I\, ,\,0\,]\\ \partial_t [\Phi_{\varpi},\Phi] \\ \Psi_\varpi,\Psi \\ \Psi_{\varpi\mu},\Psi_{\mu}\end{bmatrix} \quad  = \quad
\begin{bmatrix} 0 & 0 & 0 & 0 \\ H_{\varpi} & H  & B_b &0 \\
                C^{\dagger}_{\varpi}&C^{\dagger} & D & 0 \\
                c_1 & c_2 & d_{\varphi\mu} & 0\end{bmatrix} \quad 
\begin{bmatrix}I,0 \\ \Phi_{\varpi},\Phi  \\ \Psi_{\varpi},\Psi \\ \Psi_{\varpi\mu},\Psi_{\mu} \end{bmatrix}\quad +
  \quad\begin{bmatrix} 0,0\\ \textit{termes}  \\ \textit{sources} \\0,0\end{bmatrix}
\label{sltc4} 
\end{equation}
%
dans lequel on a mis à part les colonnes relatives aux paramètres~($\varpi$).

On peut donc oublier la première ligne du système. La sous-matrice $H_{\varpi}$
n'intervient que dans la partie explicite ($\Gamma=g(t)\dt$) de ce système linéaire\index{sensibilité}.
On a donc simplement à résoudre le système étendu suivant~:
%
\begin{equation}
\begin{bmatrix}   A  & B &0 \\
                -C^{\dagger} & I-D & 0 \\
                -c_2 & -d_{\varphi\mu} & 0\end{bmatrix} \quad 
\begin{bmatrix}  \delta\Phi  \\ \delta\Psi  \\ \delta\Psi_{\mu} \end{bmatrix}\quad =
  \quad\begin{bmatrix}  \Gamma_{\Phi}  \\ \Omega_{\Psi} \\ 0\end{bmatrix}
\label{mini3} 
\end{equation}
%
commentaires~:
%
\begin{itemize}
\item la partie explicite de la première ligne est $\Gamma_{\Phi} = \dt[H\Phi+B_b\Psi]+$ sources
et pour les colonnes relatives à $\varpi$~:~ $\Gamma_{\varpi\Phi} = \dt[H_\varpi+H\Phi_\varpi+B_b\Psi_\varpi]$;
\item la détermination de $\Psi$ se fait par $(I-D)\Psi = C^{\dagger}\Phi +$ sources;
\item celle de $\Psi_\varpi$ par $(I-D)\Psi_\varpi = c1 + c_2\Phi_\varpi $;
\item le terme supplémentaire est donné par $\Psi_{\mu} =  c_2\Phi+d_{\varphi\mu}\Psi$
et  $\Psi_{_\varpi\mu} = c_1+ c_2\Phi_\varpi+d_{\varphi\mu}\Psi_\varpi$
 (cf~\ref{sltc4}), et ainsi la dernière ligne de~(\ref{mini3}) n'est pas utilisée~(on ne calcule pas les $d\Psi_\mu$);
\item les termes sources des vecteurs $\Gamma_{\Phi}$ er $\Omega_{\Psi}$ sont à définir colonne par colonne
en fonction du type de sensibilité choisi;
\item dans le cas des sensibilités à {\bf $n_I$} C.I., c'est $n_I$ colonnes de $\Phi(0,0)$ qui sont
mises à un coefficient 1;
\item le cas des sensibilités à une perturbation de $n_f$ composantes de transfert mène à au moins
deux cas standards déjà envisagés (cf Effluents), \cad un step avec source unité (qui ira dans une colonne
de $\Omega$), et celui d'un pulse -- qui consistera à remplir une colonne de $\Gamma$ par une colonne
de $(I-D)^{-1}B$ à l'instant du pulse; on pourra étendre à des perturbations supplémentaires (rampe, etc);
\end{itemize}
%
Par ailleurs, des raisons pratiques -- telle l'existence d'un transitoire initial -- il est souhaitable
de pouvoir choisir un premier instant d'appplication des sources ({\bf ... at 100;}).
Même en ce qui concerne la sensibilité
aux C.I., on peut être amené à se désintéresser du transitoire. Par cohérence, on donne la même possibilité
aux sensibilités paramétriques\index{paramètres}.

Enfin, on offre également la possibilité de ré-initialiser ces sensibilités à intervalles réguliers par
({\bf ... at 12 every 100;}). Cette option automatique peut toujours être effacée par une gestion
directe des flags de contrôles (variables logiques {\bf Zo$\_$p(.)}) dans la séquence {\bf ZSteer}).
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Les variables et symboles, rangements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On fait ici le lien entre le formulaire et les tableaux Fortran du calcul, avec la spécificité
de \Minik qui utilise un code symbolique dont certains symboles seront remplacés par
des expressions ou symboles différents en fonction des macros générées par le {\bf ZINIT}.
Pour simplifier, on reporte à la section de l'adjoint la description des éléments dépendants
de la sensibilité de la fonction de coût\index{fonction de coût} aux variables (adjointes),
à la commande et aux paramètres.

On commence par un tableau des variables de type état et transferts -- incluant donc les
extensions~:
%
\begin{flushleft}
\begin{tabular}{||p{1.4cm}||*{4}{c|}|}
\hline
{\it nom} & dim 1 & dim 2 & {\it signification} & data file\\
\hline
\hline
eta  & np & 1 & vecteur d'état & res.data\\
\hline
dot\_eta  & np & 1 & vel d'état & no data\\
\hline
dneta  & np & 1 & avance d'état & dres.data\\
\hline
Phi$\_$t  & np & nxp+1 & matrice des eta-sensibilités aux variables &sens.data \\
\hline
dPhi$\_$t  & np & nxp+1 & matrice des eta-sensibilités aux variables &no data \\
\hline
Phi$\_$p  & np & lp+1 & matrice des eta-sensibilités aux paramètres &sensp.data\\
\hline
dPhi$\_$p  & np & lp+1 & matrice des eta-sensibilités aux paramètres &no data\\
\hline
ff  & mp & 1 & vecteur de transfert & tr.data\\
\hline
ff(mp+1)  & mobs & 1 & vecteur des observables & obs.data\\
\hline
Psi$\_$t  & mp & nxp+1 & matrice des ff-sensibilités aux variables  & sigma.data\\
\hline
Psi$\_$p  & mp & lp+1 & matrice des ff-sensibilités aux paramètres & sigmap.data\\
\hline
Pso$\_$t  & mobs & nxp+1 & matrice des mu-sensibilités aux variables  & sigmo.data\\
\hline
Pso$\_$p  & mobs & lp+1 & matrice des mu-sensibilités aux paramètres & sigmop.data\\
\hline
\hline
\end{tabular}
\end{flushleft}
%
On remarque que les dimensions $+1$ ne sont là que pour éviter une erreur de compilation
avec une déclaration de dimension nulle. La solution implémentée est de déclarer que la
colonne 0 d'un tableau est ''equivalent'' à la dernière du précédent :
%
\begin{verbatim}
     real Phi_t(np,0:nxp),Psi_t(mp,0:nxp);
     real Phi_p(np,0:nxp),Psi_p(mp,0:nxp);
     equivalence (Phi_t(1,0),eta(1)),(Psi_t(1,0),ff(1));
     equivalence (Phi_p(1,0),Phi_t(1,nxp)),(Psi_p(1,0),Psi_t(1,nxp));
\end{verbatim}
%

Un ensemble de symboles Mortran\index{macros Mortran} sont dans les instructions de calcul des matrices  et vecteurs
de {\bf principal}, qui seront remplacées par des expressions ou des symboles différents selon
les directives générées à partir du {\bf ZINIT}. Les voici~:
%
\begin{flushleft}
\begin{tabular}{||p{4.0cm}||*{4}{c|}|}
\hline
{\it nom} & {\it signification} & origine du remplacement & selflag \\
\hline
\hline
deta$\_$tef[i]  & $g_i(\eta,\varphi)$ & fun: de set$\_$eta $<...>;$ &\\
\hline
deta$\_$tef[i]  & $g_i(\eta,\varphi)$ & fun: de set$\_$dwn$\_$eta $<...,>;$ & Grid1D\\
\hline
nde$\_$eta$\_$fun$\_$[.i](inode)  & $g_{i+nlo}(\eta,\varphi)$ & fun: de set$\_$node$\_$eta $<...,>;$ & Grid1D\\
\hline
deta$\_$tef(up[.i])  & $g_{i+o}(\eta,\varphi)$ & fun: de set$\_$up$\_$eta $<...,>;$ & Grid1D\\
%%                     & (avec un offset : o) &&\\
\hline
(/eta(np+[.j]))  &  composante de $\varpi$ & free$\_$parameters: liste & tjrs\\
\hline
\hline
dphi$\_$tef[i]  & $f_i(\eta,\varphi)$ & set$\_$phi $<...,$fun:...; ...$>;$ &\\
\hline
phi$\_$tef[i]  & $f_i(\eta,\varphi)$ & fun: de set$\_$dwn$\_$phi $<...,>;$ & Grid1D\\
\hline
nde$\_$ff$\_$fun$\_$[.i](inode)  & $f_{i+nlo}(\eta,\varphi)$ & fun: de set$\_$node$\_$phi $<...,>;$ & Grid1D\\
\hline
phi$\_$tef(up[.i])  & $f_{i+o}(\eta,\varphi)$ & fun: de set$\_$up$\_$phi $<...,>;$ & Grid1D\\
\hline
\hline
\end{tabular}
\end{flushleft}
%                       

\newpage
Les matrices du TEF\index{TEF}~:
%
\begin{flushleft}
\begin{tabular}{||p{1.9cm}||*{4}{c|}|}
\hline
{\it Jacobiennes} & sous-matrices & {\it signification} & TEF-mx\\
\hline
H  & H(1:np, 1:np) & $\partial_{\eta}g(\eta,\varphi)$ & $A=I-\frac{\dt}{2}H$\\
\hline
   & H(1:np, np+1:np+lp) & $\partial_{\varpi}g(\eta,\varphi,\varpi)$ & id\\
\hline
\hline
Bb  & Bb(1:np, 1:mp) & $\partial_{\varphi}g(\eta,\varphi)$ & $B = -\frac{\dt}{2}Bb$\\
\hline
\hline
Ct  & Ct(1:mp, 1:np) & $\partial_{\eta}f(\eta,\varphi)$ & Ct \\
\hline
    & Ct(1:mp, np+1:np+lp) & $\partial_{\varpi}f(\eta,\varphi,\varpi)$ & Ct \\
\hline
    & Ct(mp+1:mp+mobs, 1:mp) & $c_2=\partial_{\eta}f_{\mu}(\eta,\varphi)$ & Ct \\
\hline
    & Ct(mp+1:mp+mobs, np+1:np+lp) & $c_1=\partial_{\varpi}f_{\mu}(\eta,\varphi,\varpi)$ & Ct \\
\hline
\hline
D  & D(1:mp, 1:mp) & $\partial_{\varphi}f(\eta,\varphi)$ & D \\
\hline
   & D(mp+1:mp+mobs, 1:mp) & $\partial_{\varphi}f_{\mu}(\eta,\varphi,\varpi)$ & D \\
\hline
\hline
\end{tabular}
\end{flushleft}
%                       

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le découpage du code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Les séquences principales programmées par l'utilisateur sont au nombre de trois~:
%
\begin{itemize}
\item la séquence ZINIT est la principale, où le modèle mathématique est décrit, les variables
et paramètres initialisés en pseudo-Fortran et les options de calcul choisies;
\item  ZSTEER est exécutée en fin de boucle temporelle pour offrir des possibilités
de gestion (printout, données, calculs diagnostics particuliers, ...;
\item DIMETAPHI est la séquence de déclaration des paramètres de dimensionnement du modèle;
\end{itemize}
%
De plus, la programmation utilise de nombreuses séquences pour appliquer le principe de non-duplication;
la liste en est donnée ci-après. Cette option de découpage du code de calcul en sous-parties de programme
permet à l'utilisateur d'en importer dans sa source de code propre de façon à en modifier certaines
sans pour autant se détacher de l'évolution du reste du code.

Enfin, ce découpage du code est optionnel pour certaines parties de calculs spécifiques.
Ainsi, le code relatif au filtre de Kalman est activé par un flag logique ({\bf Kalman}, et n'est
généré en Fortran que si ce flag est activé (par {\it sel} ou {\it selmod}).
%
\subsection{Tableau des séquences}
\label{tabseq}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Par ordre de première apparition dans le programme principal, liste des séquences indépendantes
des flags logiques.

\begin{flushleft}
\begin{tabular}{||p{3.7cm}||*{2}{c|}|}
\hline
{\it nom} & {\it contenu} & {\it dans quels codes}\\
\hline
\hline
FTN$\_$opt  &implicit simple ou double précision & général\\
\hline
FTN$\_$OPT$\_$Principal   & restriction FTN option&principal \\
\hline
LO$\_$UNITS   & unité logiques FTN& partout\\
\hline
Definition$\_$Dimension & macros {\bf real}$\rightarrow$ {\bf double}& principal, subroutines\\
\hline
DimEtaPhi & paramètres de dimensionnement & principal, subroutines\\
             & du modèle &  et autres programmes\\
\hline
Definition$\_$Variable    &déclaration des variables& principal, subroutines\\
                          &et dimensionnements induits & \\
\hline
Derive$\_$Mac    & macros MTN dérivation formelle& principal\\
\hline
Set$\_$Eta$\_$Phi   & macros MTN descro modèle& principal\\
\hline
Print$\_$Hlp     & construction de Modele.hlp & principal\\
                 & et helps printout & \\
\hline
Z$\_$Pr          & macro facilité de printout& partout \\
\hline
ZINIT         &programmation du modèle & principal\\
\hline
ZCmd$\_$law      &loi de commande éventuelle & trajectoire, adjoint\\
\hline
matrice$\_$D     & calcul D-Loop& {\it ibidem}\\
\hline
D$\_$mx           & dérivation matrice D& trajectoire, adjoint \\
\hline
matrice$\_$ABCt   & dérivation des Jacobiennes & trajectoire, sensibilité \\
                  &  H, Bb et Ct               &  et adjoint\\
\hline
A$\_$Mx           & dérivation de Jacobienne H& {\it ibidem}\\
\hline
B$\_$Mx           & dérivation de Jacobienne Bb& {\it ibidem}\\
\hline
Ct$\_$Mx           & dérivation de Jacobienne Ct& {\it ibidem}\\
\hline
Gamma$\_$vec     & calcul vecteur $\Gamma$& traj., sensib. et adjoint\\
\hline
Omega$\_$vec     & calcul vecteur $\Omega$& {\it ibidem}\\
\hline
OmeGam$\_$sup    & analogues {\it ibidem}& extensions\\
\hline
Borel$\_$vecs    & vecteurs fonction de $\dt$& Borel sweep \\
\hline
ZSTEP         &(vide) avant incrémentations fin pas de temps&principal \\
\hline
ZSTEER         & gestion fin pas de temps&principal \\
\hline
prt$\_$mat        &print Jacobiennes et vecteurs &générale  \\
\hline
Cout$\_$Psi       & Fonction de coût terminal& direct et adjoint\\
\hline
Cout$\_$l         &intégrande fonction de coût  & trajectoire et adjoint\\
\hline
trans$\_$ABCD     & transposition des matrices du TEF& adjoint\\
\hline
\end{tabular}
\end{flushleft}
%

Il est toujours possible d'importer une séquence dans son {\bf patch zinproc} pour
en modifier le contenu, qui est alors prioritaire dans la génération du Fortan.
Ainsi de la séquence (102) ZSTEP, permettant d'introduire un critère sur {\bf dt},
sans incrémenter et en renvoyant à {\bf :ReDoStep:};
%
\subsection{Séquences supplémentaires}
Cette liste complémentaire concerne des parties du code générées sur activation
de flags logiques~:

\begin{flushleft}
\begin{tabular}{||p{4.0cm}||*{2}{c|}|}
\hline
{\it nom} & {\it contenu} & {\it options}\\
\hline
\hline
varcov$\_$update  &mise à jour matrice de variance-covariance & Obs\\
\hline
dPi$\_$HBbF  &sensibilité Jacobienne à un paramètre & dPi$\_$Aspha\\
\hline
dPi$\_$CtD  &{\it idem} & {\it idem} \\
\hline
dffk$\_$HBb  &{\it idem} & {\it idem}\\
\hline
SELF-code,IF=Kalman  &avance matrice $P$ et calcul de $K$ optimal & Kalman\\
\hline
\hline
\end{tabular}
\end{flushleft}
%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Les algorithmes}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------
\subsection{Calcul de la trajectoire}
%---------------------------------------------
%
Enchaînement des calculs :
\begin{itemize}
\item le ZINIT~: description du modèle et initialisations;
\item {\bf début de la boucle temporelle $\qquad\longleftarrow\bullet$};
\item La {\bf D-LOOP}~: résolution du système implicite $\varphi = f(\eta,\varphi)$;
\item calcul des matrices Jacobiennes $H,Bb,Ct$ et $D$ par dérivation symbolique;
\item  extensions ~: calcul des $\Psi$ (Phi$\_$t et Sens$\_$p);
\item initialisation éventuelle particulière de $\Phi$;
\item sortie des variables~(qui sont en phase après la D-loop);
\item calculs diagnostics (matrice d'avance d'état $A^{spha}$, fonction de coût $J$) et sauvegardes pour l'adjoint;
\item test de linéarité de l'avance de $\varphi$ ($f(t+\dt)\leftrightarrow f(t)+d\varphi$);
\item  (boucle sur $\dt$ si balayage de Borel$\rightarrow$ gain);
\item calcul des vecteurs et matrices du TEF (intervention de $\dt$);
\item préparation analogue pour le calcul des extensions;
\item {\bf appel à  ker pour la résolution};
\item séquence ZSTEP (pour gestion de $\dt$;
\item  écriture fichier de $\dt, \deta$;
\item avance $\eta$ et $t,\dt$ (et extensions);
\item passage par la séquence ZSTEER;
\item {\bf fin du pas de temps, boucle $\qquad\bullet\nearrow$};
\item fermeture fichier et arrêt.
\end{itemize}
%  ............................................
\subsubsection*{Le ZINIT}
%  ............................................
Cette séquence où l'utilisateur définit les modèles, fournit les conditions initiales,
et explicite les options de calcul ne présente de difficulté qu'en ce qui concerne
le langage, \cad le traitement par {\bf Mortran} du code entré.
Les explications en sont fournies dans un autre cahier~:
 ``Les macros Mortran''\footnote{en progression.}\index{macros Mortran}. Nous ne donnerons ici que les principes.

On considère le code suivant, qui déclare les modèles de cellule de l'exo Lorenz~:
%
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%
! Cell definition
!%%%%%%%%%%%%%%%%%%%%%%

set_eta
< var: eta_courant_L,
  fun: deta_conv = pi_prandtl*ff_T_czcx
                                  - pi_prandtl*eta_courant_L;
  var: eta_T_czcx,
  fun: deta_energy_zx = - ff_Psi_Tsz + pi_Rayleigh_ratio*ff_courant_L
                                                 - eta_T_czcx;
  var: eta_T_sz,
  fun: deta_energy_z = ff_Psi_Tczcx - pi_wave_nb_ratio*eta_T_sz;

>;
\end{verbatim}
%
On a ainsi défini trois variables d'état, chacune d'elles associée à un ``modèle'', \cad
un équation d'évolution $\partial_t\eta=g(\eta,\varphi)$. L'expression $g()$ sera dérivée
par rapport aux composantes de $\eta$ (matrice $H$) et $\varphi$ (matrice $B_b$).
Le calcul de ces dérivées est écrit en utilisant les symboles standards de variables
d'état et sa formule correspondante~:
%
\begin{verbatim}
!--------------------------------------------------------
! matrice    H(i,j) = dEta(i)/dEta(j)
!--------------------------------------------------------
       do j=1,np
       <
         ;mult[i]=1,np :  H([i],j) = f_d(deta_tef[i])(/eta(j));
       >;
\end{verbatim}
% 
En premier lieu, le symbole {\bf np} de la {\bf mult} est remplacée par ``3'', grâce à des macros
travaillant sur le code de {\bf dimetaphi} (explicite ou généré dans la 102).
Ensuite, la macro {\bf mult} génère trois lignes avec incrémentation de [i] de 1 à 3
(La {\bf set$\_$eta} a mis un compteur à la valeur 1).

Une des trois lignes est à présent
%
\begin{verbatim}

  H((2),j) = f_d(deta_tef(2))(/eta(j));

\end{verbatim}
% 
les macros {\bf var:} et {\bf fun:} on l'effet
%
\begin{itemize}
\item d'associer {\bf deta$\_$tef(2))} à la formule déclarée en deuxième position`:
%
\begin{verbatim}

  H((2),j) = f_d( - ff_Psi_Tsz 
                          + pi_Rayleigh_ratio*ff_courant_L
                                                - eta_T_czcx)
                                                             (/eta(j));

\end{verbatim}
% 
\item ensuite, chaque symbole déclaré {\bf var} dans le modèle sera remplacé par la composante
{\bf eta(.)} correspondante -- et {\bf ff(.)} pour les {\bf set$\_$phi};
%
\begin{verbatim}

  H((2),j) = f_d( -ff(2) +pi_Rayleigh_ratio*ff(3) -eta(2))(/eta(j));

\end{verbatim}
% 
\item il reste alors a effectuer la dérivation {\bf (/eta(j))}.
\end{itemize}
%
Dans la version historique du {\bf ZINIT}, on déclarait l'équivalence entre un nom
de formule (Phi$\_$TEF) et son expression avec un set$\_$f. Avec le langage plus avancé,
on a simplelement une étape supplémentaire faisant intervenir un troisième symbole
déclaré par {\bf fun: symbole = expression;}.

Tout le jeu des macros générées par les {\bf set, var et fun} s'explique par cet objectif
de remplacement, la
complexité apparente des macros tenant aux règles de parcours des symboles dans Mortan.
Bien sûr, la séquence décrite n'a pas lieu dans cet ordre, et chaque ligne [i] générée est d'abord
travaillée par des macros remplaçant [i] par un compteur, puis soumise aux macros générées
par le {\bf ZINIT}, avant le passage au compteur suivant, avec procédure d'arrêt, etc, ce qui fait
qu'un {\bf \&T3} de Mortan sur un calcul est un peu pénible à suivre -- ce n'est pas une raison
pour s'en priver.

De plus, les macros génèrent du code informatif ({\bf help code}) et appellent une routine
{\bf prohlp} qui a pour travail de conserver ce lien entre composante de vecteur et symbole-utilisateur
pour informer les sorties du code (et repérer les composantes des vecteurs et matrices
imprimées). Enfin, un {\bf modele.hlp} rassemble les formules du modèle.

\medskip
{\it On conçoit alors d'emblée l'importance de la définition d'un langage, puisque
non seulement il permet de séparer la formulation algorithmique du code des symboles du calcul,
mais de plus il permet de générer du code Fortran supplémentaire pour vérifier, conserver
et transmettre de l'information tout au long de la chaîne de calcul, dans les print$\_$out et
dans les fichiers résultats.
 Au moment où ces lignes sont écrites, seul Mortran offre
d'aussi vertigineuses potentialités\footnote{il est utile ici de citer un cadre du CNES
qualifiant Mortan~:\\ \emph{ ``c'est chiant, c'est vieux,..., on a pas trouvé mieux''} (dixit vip, Juillet 2006).}.}
%
%  ............................................
\subsubsection*{La $D$-Loop }
%  ............................................
La résolution de l'équation implicite des transferts $\varphi=f(\eta,\varphi)$ doit
être effectuée dès que l'état est connu, avant calcul des Jacobiennes\index{Jacobienne} -- le modèle
non linéaire faisant {\it a priori} dépendre ces matrices de $\varphi$ et de $\eta$.

C'est le cas dès la sortie de {\bf ZINIT} -- les C.I. étant données -- et après
le calcul d'avance d'état (après résolution par {\bf call ker} et $\eta = \eta +\deta$.
Il n'est donc pas étonnant que la boucle temporelle débute par cette séquence.
Par ailleurs, le calcul itératif de cette équation implicite doit éviter de possibles
singularités en $\varphi$ -- intervention de {\bf $\frac{1}{ff([i])}$} par exemple, ce qui
pourra exiger de fixer des valeurs initiales particulières à certaines composantes
de transfert, d'où l'algorithme suivant.

En initialisation de la variable $\varphi$~({tableau \bf ff}) lors de l'entrée
dans la boucle d'avance de la trajectoire, on conserve les valeurs
données par ZINIT.Par ailleurs, lors de la résolution du système implicite $\varphi=f(\eta,\varphi)$,
il est numériquement nécessaire que la solution reste indépendante de l'ordre
de calcul des composantes du vecteur $\varphi$. Ceci se fait par recopie de {\bf ff} dans
le tableau {\bf phiz} et calcul de {\bf ff(.) = f(eta,phiz}). La fonction $f()$ est
en correspondance symbolique avec {\bf Phi$\_$tef[.]}\footnote{cf dérivation formelle~\ref{deriv}.}.

La relaxation est ensuite effectuée par une méthode de Newton-Raphson, $\eta$ étant fixé~:
\begin{eqnarray}
\varphi_0 +\dfi &=& f(\eta,\varphi_0) + D\dfi\nonumber\\
\rightarrow \dfi &=& [I-D]^{-1}[f(\eta,\varphi_0) - \varphi_0] \nonumber
\end{eqnarray}
%
avec $D=\partial_{\varphi}f(\eta,\varphi)|_{\varphi_0}$. La variable vectorielle $\varphi_0$
est dans le tableau {\bf phiz} (routine {\bf newt}). Il faut donc calculer la matrice $D$ en
dérivant {\bf f($\varphi_z)$} en fonction de {\bf phiz}, d'où les macros MTN transformant
les formulations de {\bf Phi$\_$tef(ff)} en {\bf Phi$\_$tef(phiz)}, qui sont annulées en sortie de
la D-Loop (``retour à ff'').

 La relaxation est itérée jusqu'à
satisfaction du critère $f(\varphi)-\varphi_0 < TOLF$ ou  $\dfi < TOLX$, ou abandon après
{\bf NTRIAL} essais (avec diagnostic). En option {\bf Debug}, on a en plus une estimation
du conditionnement de $[I-D]$ et sortie de chaque itération.
%  ............................................
\subsubsection*{Fin du pas de temps et sorties des variables }
%  ............................................
La résolution du système implicite en $\varphi$ termine le calcul d'un pas de temps, puisqu'alors
les variables $\eta$ et $\varphi$ sont en phase. Les fichiers {\bf .data} sont alors abondés
d'une ligne pour chaque vecteur (avec la variable {\bf time} en tête~:
%
\begin{itemize}
\item {\bf ff} vecteur $\varphi$ et {\bf obs} vecteur $\mu$;
\item {\bf eta} vecteur $\eta$, suivi du logarithme du Wronskien~(cf SLT$\_$Al1 \ref{SLT-wrons});
\item {\bf Phi$\_$t} et le couple {\bf Psi$\_$t, PsiMu$\_$t} des sensibilités aux variables;
\item {\bf sensibilités paramétriques} dans les tableaux {\bf Sens$\_$p} et {\bf Sigma$\_$p};
\item {\bf dnEta}, la variation de l'état, correspond en sortie (avec la valeur de l'incrément
 {\bf dt}) à l'avance entre {\bf time} et {\bf time+dt};
\item {\bf autres incréments} de même que ci-dessus, {\bf dPhi$\_$t} et {\bf dSens$\_$par}
 sont sortis;
\item de manière provisoire pour vérifications, on sort aussi {\bf Phi$\_$t+dPhi$\_$t};
\end{itemize}
%
Tous les résultas sont écris dans des fichiers {\bf .data} assortis éventuellement d'un
suffixe différent ({\bf osuffix}) comme décrit dans le Manuel pratique.
%  ............................................
\subsubsection*{La dérivation des Jacobiennes }
\label{deriv}
%  ............................................
Un ``modèle'', associé à une variable, correspond à chaque ligne d'une Jacobienne. La formule
associée devra être dérivée par rapport à toutes les variables $\eta, \varphi, \mu, \varpi$
dans des matrices différentes. On distingue donc dans le code la génération d'une formule [i]
en correspondance avec le modèle d'une variable [i], et sa dérivation ($\backslash$var(j))
 qui est effectuée dans une boucle sur (j), comme explicité maintenant.

Le calcul des coefficients des matrices Jacobiennes\index{Jacobienne} est effectué dans
une boucle Fortran par lignes, chaque ligne correspondant á un ``modèle'' différent~:
%
\begin{verbatim}
!---------------------------------
! Matrice d(i,j) = dPhi(i)/dPhi(j)
!---------------------------------
        do j=1,mp
        < mult [i]=1,mp : d([i],j) = f_d(Phi_tef[i])(/ff(j));
        >;
;
\end{verbatim}
%
Une ligne est calculée par l'expansion du code {\bf mult} qui génère $i=1,mp$ instructions
Fortran. Le symbole {\bf Phi$\_$tef[i]} est ainsi remplacé d'abord par l'instanciation
{\bf Phi$\_$tef[1]}, ce dernier symbole lui-même remplacé par sa définition dans ZINIT,
\cad soit l'expression à droite du signe égal dans le cas d'une déclaration explicite
de type {\bf Phi$\_$tef[1] = expression FTN}, soit par un {\bf set$\_$phi} qui lui-même aura
généré une déclaration explicite de remplacement.

Enfin, la macro {\bf f$\_$d(S1)(/S2)} effectue la dérivation
de l'expression mathématique dans {\bf S1} par le symbole {\bf S2} - soit ici {\bf ff(j)}.
Remarquons à ce sujet que la dérivation d'une expression fonction de ff([i]) par une
variable indicée -- ici ff(j) -- est effectuée
en dérivant par la chaîne de caractères ff(), et le résultat multiplié par la fonction
implicite {\bf F$\_$delta([i],j)} nulle si $[i]\ne j$, égale à 1 sinon. Autrement dit,
lorsqu'on dérive une expression par ff(j) dans une boucle {\bf Do} portant sur $j$
-- c'est le cas de l'exemple ci-dessus, le code Fortran est identique, seule la valeur de
{\bf F$\_$delta} annule ou non cette dérivée.


\medskip {\it cas du Grid 1D}

On a d'une part les transferts {\bf dwn} et {\bf up} qui sont semblables aux standards~:
%
\begin{verbatim}
!---------------------------------
! Matrice d(i,j) = dPhi(i)/dPhi(j)
!---------------------------------
;
        if m_dwn .gt.0
        < do j=1,mp
          <  ;mult [i]=1,m_dwn :     d([.i],j) = F_D(Phi_tef[i])(/ff(j));
          >;
        >;
        do j=1,mp
        < do inode=1,m_node
          < ;mult [i]=1,m_mult :
              d(k2j_[.i](inode),j) = f_d(nde_ff_fun_[.i](inode))(/ff(j));
          >;
        >;
        if m_up .gt.0
        < do j=1,mp
          <  ;mult [i]=1,m_up :
                      d(joffsetup+[.i],j) = F_D(Phi_tef(up[.i]))(/ff(j));
          >;
        >;
;
\end{verbatim}
%
Pour les {\bf nodes}, on a  {\bf m$\_$node} séquences d'expressions identiques ne
différent que de la valeur de l'indice de boucle {\bf inode}. C'est encore la fonction
{\bf F$\_$delta(inode,j)} qui annulera ou pas la dérivée.

\medskip {\it les autres Jacobiennes et le passage aux matrices et vecteurs d'avance}

Les matrices Jacobiennes A,Bb,Ct sont calculées de manière semblable à D. Le calcul
du vecteur {\bf ff} aussi, mais sans dérivation.

Dans la mesure où ces calculs de Jacobiennes interviennent dans l'adjoint, leur code
est regroupé pour chacune dans une séquence (cf tableau~\ref{tabseq}).

Enfin pour finir avec les matrices du TEF\index{TEF}, les matrices $A=I-\frac{\dt}{2}H$ et $B=-\frac{\dt}{2}Bb$
sont effectuées dans le corps de {\bf principal}~:
%
\begin{verbatim}
!====================================================
! Complément de calculs dépendants de dt
!====================================================
! build a mx(init a I) from h : a = I -dt/2*h
       call scamat(-dt/2.,h,hdt2m,n,np,n+lp,np+lp);
       call diagmat(one,a,n,np+lp);
       call matadd(a,hdt2m,a,n,np,n+lp,np+lp);
;
! Build TEF B mx from Jacobian Bd
       call scamat(-dt/2.,Bb,B,n,n,m,m);
\end{verbatim}
%
Deux autres séquences concernent enfin le calcul des vecteurs $\Gamma$ et $\Omega$~:
%
\begin{verbatim}
!*KEEP,Gamma_vec.
!--------------------------------------------------
! vecteur gamma avec : (dt)eta = deta_tef(eta,ff,t)
!--------------------------------------------------
       mult[i]=1,np : gamma[i] = dt*(deta_tef[i]
                                     + 0.5*dt*(f_d(deta_tef[i])(/time))
                                    );

!*KEND.
!*KEEP,Omega_vec.
!---------------------------------------------
! vecteur omega avec : Phi = Phi_tef(eta,ff,t)
!---------------------------------------------
      ;mult[j]=1,mp : omega[j] = dt*f_d(Phi_tef[j])(/time);

!*KEND.
\end{verbatim}
%
où seule une dérivation par rapport à la variable {\bf time} est effectuée.

Ceci termine le calcul de l'ensemble des vecteurs (2) et matrices (4) nécessaires
au calcul de la trajectoire (de $t$ à $t +\dt$).
%  ............................................
\subsubsection*{La matrice $A^{spha}$}
%  ............................................
La matrice d'avance du système d'état -- hors TEF, \cad lorsque les transferts ont été
éliminés -- présente un intérêt théorique, car dans ce cas on se retrouve dans un cas
largement étudié de l'analyse des systèmes dynamiques. Cette matrice permet par exemple
le calcul de gains de rétroaction non-stationnaires (cf Effluents). Elle est également
nécessaire pour faire avancer le calcul de la matrice de var-covar de Kalman.

Lorsque la matrice $(I-D)^{-1}$ est régulière, la matrice dîte d'avance d'état est
$A^{\sharp}=H-B_b(I-D)^{-1}C^{\dagger}$ (tableau FTN Aspha)~:
%
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Calcul de la matrice d'avance dans l'espace de phase
!       d_t Eta = [ A + B.(I-D)-1.Ct ] Eta
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       call omatmlt(m,m,n,un_moins_d,mp,ct,mp+mobs,ctbuf,mp);
       call omatmlt(n,m,n,Bb,np,ctbuf,mp,aspha,np);
       call omatadd(n,n,h,np,aspha,np,aspha,np);
       if (Zout) write(loaspha,1001) time,((aspha(i,j),i=1,n),j=1,n);
! Calcul du Wronskien [d_t Log W = Tr(Aspha)]
       oTrace_adv = 0.; <i=1,n; [oTrace_adv :+aspha(i,i)]; >;
       oInt_Tradv = oInt_Tradv + dt*oTrace_adv;
! ------------------------------------------------------------------
\end{verbatim}
%
On voit que l'on en profite pour calculer le Wronskien, \cad le logarithme
du déterminant du propagateur, qui avance par la trace de $A^{\sharp}$ (cf SLT$\_$Al1 \ref{SLT-wrons}).
%  ............................................
\subsubsection*{Borel sweep et Boreleig-Boreigrun}
%  ............................................
On explique ensemble ces deux algorithmes parce qu'ils résolvent chacun, bien
qu'avec des méthodes orthogonales, le problème de la détermination
du gain, réponse et \index{Borel} effet de rétroaction stationnaire (cf Effluents).

Le principe par {\bf Borel sweep} est de déterminer numériquement un gain
en profitant de la formulation identique de l'avance par $\frac{\dt}{2}$ et
la transformation de Borel du système local (considéré alors comme autonome)
avec la variable $\tau \leftrightarrow \frac{\dt}{2}$.

On va alors balayer en $\dt$une réoslution par {\bf ker} et conserver
le gain $g_B(\tau)$ fourni par {\bf ker}. Pour limiter les effets
de ce balayage dans {\bf principal} . on renvoie à une routine~:
%
\begin{verbatim}

! --  Balayge de Borel en dTau
       call Bsweep(H,Bb,Ct,d,deta_B,deta_tB,dPhi_tB,dt,index_ff_gain,
            itau_max,tau_B_ini,tau_B_mult,time,istep,.false.);

\end{verbatim}
%
à qui on transmet les Jacobiennes\index{Jacobienne}. Il faut cependant aussi tenir compte
des vecteurs dépendants de $\dt$, dont on calcule et transmet la dérivée
temporelle~:
%
\begin{verbatim}
!*KEEP,borel_vecs.
!          vector computations for borel sweep
!--------------------------------------------------------
! sequence de calcul des parties independantes de dt
! des vecteurs gamma omega pour exportation B-sweep
!--------------------------------------------------------

!------------------------------------------------------
! pour vecteur gamma avec : (dt)eta = deta_tef(eta,ff,t)
!-----------------------------------------------------
       mult[i]=1,np : deta_B[i] = (deta_tef[i]);
       mult[i]=1,np : deta_tB[i] =f_d(deta_tef[i])(/time);
      ;mult[j]=1,mp : dPhi_tB[j] = f_d(Phi_tef[j])(/time);

!*KEND.
\end{verbatim}
%
Le principe dans {\bf ker} de la détermination du gain consiste à
éliminer toues les composantes de transferts sauf une (spécifiée par
le pointeur {\bf index$\_$ff$\_$gain}). La variable {\bf gain} donne
en sortie la valeur de $\tilde{g}_B(\tau)$ - cf {\bf ker}.
On a donc en sortie $\mathcal{B}[g_B(\tau)$, que l'on pourra
inverser par fit etc selon la fameuse horreur borelienne.

La routine {\bf Boreleig} (pour Borel par eigen-elements) calcule
en une seule fois les éléments propres des matrices $A^{\sharp}$ et $A^{\flat}$
(cf Effluents~(\ref{P-paire1})) pour permettre le calcul en off-line des éléments
de la fonction du temps physique (constituée de exp(), sin() et cos()).
{\bf Boreigrun} utilise les coefficients calulés par la Boreleig pour
en faire une série num érique en balayant la variable temporelle pour tracer
des graphes standards (mais rien n'empèche d'utiliser la fonction temporelle
lorsqu'elle n'est pas trop longue (Gnuplot est limité à quelques dizaines de
caractères).
Boregrun se trouve dans le patch {\bf gsrun}).

Enfin, on peut utiliser avantageusement
{\bf PAW} pour visualiser les fonctions statiques telles qu'elles sont
calulées le long d'une portion de trajectoire en utilisant la routine
(hors patch) {\bf gains.f} (voir directement le code FTN). On obtient une {\bf surf2z}, voire
un exemple d'utilisation dans le cahier ``Lorenz et gains''.
%
% ------------------------------------------------
\subsection{La routine oKer, solveur du système}
% ------------------------------------------------
%
Ayant en argument d'entrée les quatre matrices d'avance $A, B, C^{\dagger}, D$ du TEF et les vecteurs
connus $\Gamma, \Omega$, le sous-programme {\bf oker} effectue la résolution du système~(\ref{mini1})
en deux temps~: élimination des cellules menant au système dit des ``couplages'' entre les transferts,
résolution des transferts, et enfin ``descente'' pour résolution des états.
On a en sortie les variations $\deta$ et $d\phi$.

En réalité, {\bf oker} étend cette résolution à toutes les colonnes supplémentaires aux
deux vecteurs du TEF, soit {\bf nwp} colonnes supplémentaires~: ({$\Gamma_x, \Omega_x$}),
exploitant ainsi l'efficacité de la méthode de résolution de système linéaire de LaPack (routine SGESV).

Enfin, l'option de balayage par {\bf Borel-sweep}\index{Borel} utilise une colonne de plus, parce que le gain
lié à une composante $i$ de transfert (integer {\bf i$\_$gain}) s'obtient par résolution du système
de couplage pour le vecteur connu $\Gamma_g=\ket{1_i}=[0,0,...,1,0,...0]^{\dagger}$ avec $1$ en position
$i$. Le gain se trouve en position $i$ du vecteur solution.

Les autres arguments d'entrée ({\bf ZPrint, eta, etax}) ne servent qu'à l'impression, et on a en
sortie les deux arguments d'information LaPack ({\bf infoa, infores}).

La seule difficulté de {\bf oker} consiste ainsi au rangement des colonnes des vecteurs connus
et cette gestion du calcul optionnel du gain. Les matrices en entrée sont elles mêmes possiblement
des sous-matrices; on esquive cette difficulté en recopiant en entrée les blocs-matrices
utilisés dans des tableaux intermédiaires.

\bigskip
Les calculs et tableaux sont organisés comme suit~:
%
\begin{itemize}
\item ``vecteurs'' connus rassemblés dans la matrice {\bf bggx} dans l'ordre {\bf bw, gamma, gammx};
de cette manière, la ``montée'' déterminant le système des couplages fournit à la fois
la matrice $A^{-1}B$ et les variations découplées $A^{-1}\Gamma,\Gamma_x$ après l'appel
à un premier {\bf sgesv} (avec le flag LaPack {\bf infoa}). On a ainsi {\bf np+1+nwp} colonnes
de résultat;
\item il reste à multiplier par $C^{\dagger}$ pour obtenir la matrice de couplage d'une part,
les évolutions insensibles ($C^{\dagger}A^{-1}\Gamma,\Gamma_x$) d'autre part
%
\begin{verbatim}

   !  ** calcul de Ct*[A-1*B&gamma&gamx] = Ct*BGGx *****
   
   call omatmlt(mp,np,mp+1+nwp,ct,mp+mobs,bggx,np,ctbggx,mp);
   
   !  *********** calcul de matrice de couplage **
   !    couplage = I - (-Ct) A-1B - D
   
   CALL omatadd(mp,mp,c_unit,mp,ctbggx,mp,unplusctamb,mp);
   CALL omatsub(mp,mp,unplusctamb,mp,d,ldd,couplage,mp);
   !                                       ========
   
   if Zprint
   <     TITLE='>ker : Matrice de couplage';
         CALL printmat(title,couplage,mp,mp,mp,mp,2,2);
   >;
   
   !  *********** calcul de dphi_ins  ***********
   !    dphi_ins = Omega - (-Ct)A-1.Gamma
   
    CALL omatadd(mp,1+nwp,ctbggx(1,mp+1),mp,omegxx,mp,dphixx,mp);
  
   if Zprint
   <     TITLE='>ker : dphi_ins';
         CALL printmat(title,dphixx,mp,mp,1+nwp,1+nwp,0,2);
   >;
\end{verbatim}
%
\item on est à présent ``en haut de navette'' dans le système de couplage, avec un
``vecteur connu'' {\bf dphixx} de {\bf 1+nwp} colonnes.
Ici intervient l'option Borel-sweep, avec création du vecteur $\ket{1_i}$ rangé en {\bf 1+nwp+1}
%
\begin{verbatim}
   if (i_gain.EQ.0)
   < call sgesv(mp,nwp+1,couplage,mp,iPiv,dphixx,mp,infores);
     if infores.ne.0
     < print*,'*** SGESV(couplage,dphixx) error: ',infores;
     >;
   >else
   < call veczero(dphixx(1,nwp+2),mp);
     dphixx(i_gain,nwp+2)=one;
     call sgesv(mp,nwp+2,couplage,mp,iPiv,dphixx,mp,infores);
     if infores.ne.0
     < print*,'*** SGESV(couplage,dphixx+gain) error: ',infores;
     >;
     gain = dphixx(i_gain,nwp+2);  "en realite, =1/(1-g)"
   >;
   ;
\end{verbatim}
%
\item descente de navette avec calcul de $A^{-1}B\,[\dfi,\dfi_x]$
%
\begin{verbatim}
   !  ********************************************************
   !  *********** calcul de a-1b*dphi,x  ***********
     CALL omatmlt(np,mp,1+nwp,bggx,np,dphixx,mp,ambdphixx,np);
   !  ********************************************************
\end{verbatim}
%
\item et enfin calcul des avances $\deta,\deta_x$ du genre
%
\begin{verbatim}
   ! Calcul avance etat
   DO I=1,np
   < dneta(i)=agamdt(i)-ambdphixx(i,1);
     if (ZPRINT)
     write(*,9001)i,agamdt(i),ambdphixx(i,1),dneta(i),eta(i),
                                            EtaNam(i)(1:lenocc(EtaNam(i)));
   >;
   9001;format(2x,I3,4(1Pe12.4),1x,A);
\end{verbatim}
%
et des extensions, avec recopie dans les arguments de sortie;
\end{itemize}
%
La routine n'est pas parfaitement optimisée, elle fait un compromis entre efficacité
et clarté-robustesse aux développements en conservant quelques copies inutiles,
mais permettant des impressions de debugging.
% ......................................................
\subsubsection*{En sortie de oKer}
% ......................................................
%
On range le calcul de l'avance $\deta$ dans son {\bf dres.data}, et on avance Etat
%
\begin{verbatim}
! Calcul avance etat,x
! --------------------
call omatsub(np,1+nwp,bggx(1,mp+1),np,ambdphixx,np,dnetax,ldeta);

if Zprint
< print9000;print*,' ---- avance d Etat - State advance -----';
  print *,' i,  A-1gamdt(i), A-1Bdphi(i),   deta(i),     eta(i) Var name';
 DO I=1,np
  < write(*,9001)i,bggx(i,mp+1),ambdphixx(i,1),dnetax(i,1),etax(i,1),
                                           EtaNam(i)(1:lenocc(EtaNam(i)));
  >;
>;
\end{verbatim}
%
de plus, on utilise le vecteur $\dfi$ pour faire avancer les transferts 
à seule fin de contrôle des calculs de la trajectoire~:
%
\begin{verbatim}
!----------------------------------------------------------
! pour test linearite : ffl pour next step
     do i=1,m
     < ffl(i) = ff(i) + dphi(i);
     >;
\end{verbatim}
%
Ceci pour comparer l'avance linéaire\index{linéarité} des transferts aux transferts calculés par
l'algorithme du TEF\index{TEF} qui sont $\varphi = f(\eta,\varphi)$. On donne ici le test effectué en fin
de pas de temps~:
%
\begin{verbatim}
! ========================================================
! test linearite ffl - ff
! ========================================================
if (istep.gt.1)
< res=0.; <io=1,m; res = res +(ffl(io)-ff(io))**2/max(one,ff(io)*ff(io)); >;
  if (res .gt. TOL_FFL)
  < print*,'*** linearity pb: res > TOL_FFL at istep',istep,res,' > ',TOL_FFL;
    write(*,*)'  i,   ff(i)  , ff(i)-ffL(i), var Name';
    do io=1,m < write(*,'(i4,2(1pe12.5,1x),a)')
                                       io,ff(io),ff(io)-ffl(io),PhiNam(io); >;
    if (idxerror.eq.0.) idxerror = -1;
  >;
>;
\end{verbatim}
%
où on voit qu'il repose sur l'écart quadratique ``quasi-relatif'' entre les deux calculs de $\varphi$,
comparé au seuil {\bf TOL$\_$FFL} fixé soit dans ZINIT, soit à son défaut ($10^{-6}$), ceci
remarquons le indépendamment de la précision (simple ou double) des calculs\footnote{ce qui est la
méthode correcte pour confronter les deux types de calculs.}.
%
% ------------------------------------------------
\subsection{Calcul des sensibilités des variables}
% ------------------------------------------------
%
Comme on l'a expliqué dans le formulaire liminaire, l'ensemble de ces calculs
de sensibilité des variables $\eta, \varphi, \mu$ sont regroupés colonne par colonne
dans des tableaux en extension des vecteurs de variables;
% ............................................
\subsubsection*{Les extensions}
% ............................................
Description des équivalence et dimensions;

On va séparer l'explication des calculs des sensibilités aux variables et aux paramètres.
% ......................................................
\subsubsection*{Calcul des sensibilités à des variables}
% ......................................................
On a d'abord le calcul du $\Gamma_{\Phi}$\index{sensibilité}
%
\begin{verbatim}
! ---------------------------------------------------------
! extension Gamma_Phi_t = dt(H Phit + Bb Psit [+source])
     if nxp.gt.0
     < call omatmlt(n,n,nxp,H,np,Phi_t(1,1),np,gam_Phit(1,1),np);
       call omatmlt(n,m,nxp,Bb,np,Psi_t(1,1),mp+mobs,bobuf,np);
       call omatadd(n,nxp,gam_Phit(1,1),np,bobuf,np,gam_Phit(1,1),np);
       call oscamat(dt,n,nxp,gam_Phit(1,1),np,gam_Phit(1,1),np);
     >;
\end{verbatim}
%
Ensuite, le traitement du cas 4 (g/(1-g)) est effectué avant résolution et
avance~:
%
\begin{verbatim}
     if nxp.gt.0
     < jo=0;
       Do ko=1,nxp
       < if (istep.ge.ko_p_ini(ko).and.ko_p_typ(ko).eq.4) jo=1;
       >;
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!  Source pour propag. de transfert marche : gam_x  -> + |BbDm1>_k dt
!              de type g/1-g
! --------------------------------------------------------------------
       if jo.eq.1
       < call matmlt(Bb,un_moins_D,bobuf,n,np,m,mp,nxp,nxp);
         call scamat(dt,bobuf,bobuf,np,np,nxp,nxp);
         Do ko=1,nxp
         < if ko_P_typ(ko) .eq. 4
           < call oscamat(step_psi(ko),n,1,bobuf(1,ko_p_kvar(ko)),np
                                                ,bobuf(1,ko_p_kvar(ko)),np;
             call matadd(bobuf(1,ko_p_kvar(ko)),gam_Phit(1,ko)
                                                ,gam_Phit(1,ko),n,np,1,nxp);
           >;
         >;
       >;
     >;"End nxp"
! =================================
!      call oker
! =================================
! avance Phi_t
! ------------
   if (nxp.gt.0) call matadd(Phi_t,dPhi_t,Phi_t,np,np,nxp,nxp);
\end{verbatim}
%
Remarquer la multiplication par step\_psi de l'amplitude bobuf pour obtenir
la réponse effective à un step d'amplitude différengt de l'unité.
 Il reste à compléter la partie $\Psi$, avec les sources pulsées~:
%
\begin{verbatim}
! fin du pas de temps :
! ================================
!      D-Loop, Jacobiennes, puis :
! ================================
! Extension : calcul de (I-D) Psi_t = Ct Phi_t  [+source]
! ---------------------------------------------------------
     if nxp.gt.0
     <
! --- suivant indications des ko_p var, raz Phi_t, et sources sensib
        call pSteer("in:" np,nxp, Bb,np, un_moins_d,mp,
                          ZPrint,lp,mobs,EtaNam,PhiNam,ParNam,
                                                        "out:" Phi_t,np);
! --- calcul de Ct Phi_t [+source]
       call omatmlt(m,n,nxp,Ct,mp+mobs,Phi_t(1,1),np,Psi_t(1,1),mp+mobs);
! ---
!    -- add perpetual unit source for type 3 pert
       Do ko=1,nxp
       < if ko_P_typ(ko).eq.3 .and. istep.ge.ko_p_ini(ko)
         < call oscamat(step_psi(ko),m,1,s_unitm(1,ko_p_kvar(ko)),mp
                                                              ,T_tmp1(1,1),mp);
           call omatadd(m,1,T_tmp1(1,1),mp
                                     ,Psi_t(1,ko),mp+mobs,Psi_t(1,ko),mp+mobs);
         >;
       >;
! --- Calcul Dm1Psit
    call omatmlt(m,m,nxp,un_moins_d,mp,Psi_t(1,1),mp+mobs,Psi_t(1,1),mp+mobs);
;
     >;"End nxp"
\end{verbatim}
%
après quoi on complète par une multiplication par la $C^{\dagger}$ avant
de finir avec $(I-D)^{-1}$. Le buffer T\_tmp1 ne sert qu'à prendre une amplitude de step
définie dans ZINIT au lieu de l'unité.
%
Autre perturbation, autre source, l'ensemble étant géré dans la routine {\bf PSteer},
avec par exemple~:
%
\begin{verbatim}
!     --- source pour perturbation-marche
       if io.eq.1
       < print*,' --- (re) setting step(s) on transfer at:',time,istep,' -----';
         print*,'  -- for 1/(1-g) response-type -------------------------';
         print*,'  -  (setting Unit Source: Dm1 Psi_t = CtPhit + I)   ---';
         Do ko=1,nxp
         <
!    -- RAZ Phi_t, Psi_t si nouveau depart Tr-Heaviside
           if ko_P_typ(ko).eq.3
           < if Zko_p(ko)  " RAZ Phi_t, Psi_t "
             < call veczero(Phi_t(1,ko),np); call veczero(Psi_t(1,ko),mp);
              if ZPrint
              < print*,'   for Tr #:',ko_p_kvar(ko),' on:',ko,' Psit colon ---';
              >;
             >;
             if ko_P_typ(ko).eq.3 .and. istep.ge.ko_p_ini(ko)
             < call matadd(s_unitm(1,ko_p_kvar(ko))
                                          ,Psi_t(1,ko),Psi_t(1,ko),m,mp,1,nxp);
             >;
           >;"end 3"
         >;
       >;
\end{verbatim}
%
et même final que pour la source précédente.

On ne décrit pas ici la partie remise à zéro et initialisation, qui sont effectuées
en fonction du time-flag suivant~:
%
\begin{verbatim}
       Do ko=1,nxp
       < Zko_p(ko)=(istep.eq.ko_p_ini(ko)
                     .or. mod(istep-ko_p_ini(ko)+1,ko_p_mod(ko)).eq.1);
       >;
\end{verbatim}
%
qui est {\bf .True.} pour les {\bf istep = ini [mod]}, correspondant aux {\bf at every},
qui vaut en particulier {\bf .False.} pour la valeur de -1 donnée à {\bf ko$\_$p$\_$mod}
(\cad en l'absence de {\bf every}.
Il faut cependant traiter à part le cas du premier pas de boucle (\cad absence de
{\bf at} et donc de {\bf every} qui doit conduire à une action au premier pas.
%  ............................................
\subsubsection{Calcul des sensibilités aux paramètres}
%  ............................................
%
On a deux calculs de ``matrice connue'' à effectuer, simplifiés par
le fait que la sous-matrice $A_{\varpi}=I$\index{sensibilité}.

La fin du pas de temps résoud l'équation~(cf~\ref{sltc4})
%
\begin{equation}
(I-D)\Psi_{\varpi} = C^{\dagger}\Phi_{\varpi}+C^{\dagger} _{\varpi}
\end{equation}
%
de la façon suivante~:
%
\begin{verbatim}
! *** Partie parametres : calcul de (I-D) Psi_p = Ct Phi_p + Ct_pi
! ----------------------------------------------------------------
     if lp.gt.0
     < Zko_p(nxp+1)=(istep.eq.ko_p_ini(nxp+1)
                      .or. mod(istep-ko_p_ini(nxp+1)+1,ko_p_mod(nxp+1)).eq.1);
       if Zko_p(nxp+1)
       < call veczero(Phi_p(1,1),np*lp);
            print*,' ---  resetting sensitivity to Parameters at:',time,istep;
       >;
       call omatmlt(m,n,lp,Ct,mp+mobs,Phi_p(1,1),np,Psi_p(1,1),mp+mobs);
       call omatadd(m,lp,Ct(1,np+1),mp+mobs,Psi_p(1,1),mp+mobs
                                                           ,Psi_p(1,1),mp+mobs);
!  -- Dm1Psip
       call omatmlt(m,m,lp,un_moins_d,mp,Psi_p(1,1),mp+mobs,Psi_p(1,1),mp+mobs);
! --- partie Obs : Psi_mu =  c1 + c2 Phi_p + d_fmu Psi_p
       if mobs.gt.0
       < call omatmlt(mobs,n,lp,Ct(mp+1,1),mp+mobs,Phi_p(1,1),np
                                                        ,Psi_p(mp+1,1),mp+mobs);
         call omatmlt(mobs,m,lp,D(mp+1,1),mp+mobs,Psi_p(1,1),mp+mobs
                                      ,Ctbuf,mobs);
         call omatadd(mobs,lp,Ctbuf,mobs,Psi_p(mp+1,1),mp+mobs
                                                        ,Psi_p(mp+1,1),mp+mobs);
         call omatadd(mobs,lp,Ct(mp+1,np+1),mp+mobs
                                  ,Psi_p(mp+1,1),mp+mobs,Psi_p(mp+1,1),mp+mobs);
       >;
     >;"End lp"
! ---------------------------------------------------------
! **************   Ecriture resultats sur fichier  ********
\end{verbatim}
%
Où on voit que c'est les {\bf ko$\_$p$\_$ (nxp+1)} qui gèrent les {\bf at . every .} de l'ensemble
des paramètres (qui sont donc synchronisés contrairement aux sensibilités aux variables).

La partie calcul de $\Gamma_{\varpi}$ se calcule par~(cf~\ref{sltc4})
%
\begin{equation}
\Gamma_{\varpi} = [H\Phi_{\varpi} + H_{\varpi} + B_b \Psi_{\varpi}] dt
\end{equation}
%
avec sa programmation~:
%
\begin{verbatim}
! --- partie parametres ---
! -------------------------
! extension Omega du Psi_p
     if lp.gt.0
     < call veczero(omeg_Psip(1,1),mp*lp);  "Omega vec null"

! partie Gam_p = [H Phi_p + H_pi + Bb Psi_p] dt
! ---------------------------------------------
       call omatmlt(n,n,lp,H,np,Phi_p(1,1),np,gam_Phip(1,1),np);
       call omatmlt(n,m,lp,Bb,np,Psi_p(1,1),mp+mobs,bobuf,np);
       call omatadd(n,lp,gam_Phip(1,1),np,bobuf,np,gam_Phip(1,1),np);
       call omatadd(n,lp,gam_Phip(1,1),np,H(1,np+1),np,gam_Phip(1,1),np);
       call oscamat(dt,n,lp,gam_Phip(1,1),np,gam_Phip(1,1),np);
     >;"End lp"
\end{verbatim}
%
qui fixe aussi à zéro la matrice $\Omega_{\varpi}$.

La résolution est prise en compte avec les extensions aux matrices $A$ et $C^{\dagger}$.
La routine {\bf oker} voit en effet toutes les extensions~:
%
\begin{verbatim}
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! -------------------------------------------------------------------
   CALL oker("in"   n,m,nxp+lp,  A,np,B,np,Ct,mp+mobs,D,mp+mobs,
                                 Gamma,np,Omega,mp,Eta,np,ff,mp,
                            lp,mobs,EtaNam,PhiNam,ParNam,zprint,
                 "out" dnEta,dPhi, "in" 0, "out" gain,infoa,infores);
! --------------------******-----------------------------------------
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\end{verbatim}
%
et traite, en plus des vecteurs $\Gamma, \Omega$ du TEF\index{TEF}, {\bf nxp+lp}
colonnes, correspondant aux {\bf nxp} sensibilités aux variables demandées
et aux {\bf lp} sensibilités aux paramètres. La liste des symboles
est rangées dans le tableau {\bf ParNam(1:lp)} lors du traitement dans {\bf ZINIT}
de la {\bf Free$\_$parameters} liste.

{\it Remarque~: Le remplacement de {\bf ker} par {\bf oker} correrspond à l'adoption
de la norme LaPack des dimensions de matrices (Mx,LDMx), \cad la matrice suivie
du nombre des lignes décarées dans {\bf dimension Mx(LDMx,*)}};
%
% ======================================================
\subsection*{L'extension au vecteur des observables}
% ======================================================
%
Avec la {\bf version 102} apparait un nouveau vecteur $\mu$ dît des observables\index{observables},
en supposant qu'il
s'agit de transferts n'ayant aucune rétroaction sur les variables du TEF (matrices $B_b$ nuules).
On suppose de plus (au moins provisoirement) que les composantes sont indépendantes ($D_{\mu}=0$).

Tout ceci ne modifie donc en rien la résolution, mais on a à calculer en plus une série
de {\bf mobs} fonction explicites en plus, à la fin d'un pas de temps~:
%
\begin{verbatim}
! ---   Calcul des Observables-probes ---------------------
! ---------------------------------------------------------
;mult [j]=1,mobs : ff(mp+[.j]) = (Phi_tef(mp+[.j]));

\end{verbatim}
%
on peut ainsi constater que l'on place ce vecteur à la suite de {\bf ff(.)}, bien que les
fichiers {\bf tr et obs.data} restent séparés par raison de commodité.

Bien sûr, ces obervables sont sensibles aux variables du TEF et par la même, directement ou non,
aux {\bf Free-parameters}. On a donc besoin de sous-matrices Jacobiennes définies dans~(\ref{sltc4}),
programmées\index{Jacobienne} comme suit~:
%
\begin{verbatim}
 mult [i]=1,mobs : d(mp+[.i],j) = f_d(Phi_tef(mp+[.i]))(/ff(j));              <<< D_mx

!------------------------------------
! matrices Ct(i,j) = dFi(i) / dEta(j)
!------------------------------------
        do j=1,np
        < mult[i]=1,mp : ct([.i],j) = f_d(Phi_tef[i])(/eta(j));
        >;
!  partie Obs / probes
       if mobs.gt.0
       < do j=1,np
         < mult[i]=1,mobs : ct(mp+[.i],j) = f_d(Phi_tef(mp+[.i]))(/eta(j));   <<< Ct 2
         >;
       >;"End obs"
! partie parametres
! -----------------
       if lp.gt.0
       < ;mult[j]=1,lp: mult [i]=1,mp:
           Ct([.i],np+[.j]) = F_D(Phi_tef[i])(/eta(np+[.j]));
!  partie Obs
         if mobs.gt.0
         < ;mult[j]=1,lp: mult [i]=1,mobs:
           Ct(mp+[.i],np+[.j]) = F_D(Phi_tef(mp+[.i]))(/eta(np+[.j]));        <<< Ct 1
         >;
       >;
\end{verbatim}
%

On a ainsi deux séries de colonnes à calculer, qui sont rangées dans les tableaux en extension
du nombre de lignes de {\bf Psi$\_$t} et {\bf Psi$\_$p}~:
%
\begin{itemize}
\item La partie sensible aux variables
\begin{verbatim}
! --- partie Obs : Psi_mu =  c2 Phi_t + d_fmu Psi_t
       if mobs.gt.0
       < call omatmlt(mobs,n,nxp,Ct(mp+1,1),mp+mobs,Phi_t(1,1),np  "c2 Phit"
                                                  ,Psi_t(mp+1,1),mp+mobs);
         call omatmlt(mobs,n,nxp,D(mp+1,1),mp+mobs,Psi_t(1,1),mp+mobs
                                      ,Ctbuf,mobs);
         call omatadd(mobs,nxp,Ctbuf,mobs,Psi_t(mp+1,1),mp+mobs
                                                  ,Psi_t(mp+1,1),mp+mobs);
       >;
     >;"End nxp"
\end{verbatim}
%
\item La partie sensible aux paramètres
\begin{verbatim}
! --- partie Obs : Psi_mu =  c1 + c2 Phi_p + d_fmu Psi_p

       if mobs.gt.0
       < call omatmlt(mobs,n,lp,Ct(mp+1,1),mp+mobs,Phi_p(1,1),np
                                                        ,Psi_p(mp+1,1),mp+mobs);
         call omatmlt(mobs,m,lp,D(mp+1,1),mp+mobs,Psi_p(1,1),mp+mobs
                                      ,Ctbuf,mobs);
         call omatadd(mobs,lp,Ctbuf,mobs,Psi_p(mp+1,1),mp+mobs
                                                        ,Psi_p(mp+1,1),mp+mobs);
         call omatadd(mobs,lp,Ct(mp+1,np+1),mp+mobs
                                  ,Psi_p(mp+1,1),mp+mobs,Psi_p(mp+1,1),mp+mobs);
       >;
     >;"End lp"
\end{verbatim}
%
\end{itemize}
%
Enfin, une remise à zéro des sensibilités suit la même logique de déroulement pour
les options {\bf at [i] every [m];}.
%
% ======================================================
\subsection*{Sensibilité des jacobiennes à un paramètre}
% ======================================================
%
Sur activation du selflag {\bf dPI$\_$aspha}, le code de {\bf principal}
calcule la sensibilité\index{sensibilité} des Jacobiennes à un paramètre. Ces dérivées élementaires
sont alors disponibles pour permettre différents calculs de sensibilité dans
{\bf ZSTEER}.

Il faut noter que la dérivation partielle par rapport
au paramètre\index{paramètres} ne peut se faire pour une expression qui composerait un coefficient
de matrice, masi qui n'est pas disponible -- puisque celles-ci sont juste calculées
numériquement. Il est donc
nécessaire de repartir des formules du modèle (Deta$\_$tef et Phi$\_$tef) et d'en effectuer
une double dérivation symbolique.

Pour fixer un objectif et comprendre l'utilité du code suivant, On donne l'exemple
d'utilisation pour le calcul de la dérivée partielle
de $A^{\sharp}$, avec les formules sont les suivantes~:
%
\begin{eqnarray}
  d_{\varpi} A^{\sharp} =&& d_{\varpi}H + d_{\varphi}H \, d_{\varpi} \varphi \nonumber\\
                 &+& \left\{ d_{\varphi}B \, d_{\varpi} \varphi + d_{\varpi}B \right\} (I-D)^{-1} C^{\dagger}\\
                 &+& B \left\{ (I-D)^{-1} d_{\varpi}C^{\dagger} - (I-D)^{-1} d_{\varpi}D\, (I-D)^{-1} C^{\dagger} \right\}\nonumber
\label{dpiaspha}\end{eqnarray}

Et pour suivre le code nécessaire, on donne l'exemple de dérivation double du calcul
de la dérivée de $H$ et $B_b$
par rapport au paramètre (de symbole {\bf po$\_$asp} -- on aussi besoin de la dérivée
de $\varphi$~:
%
\begin{verbatim}
!*KEEP,dPi_HBbF.
       ;mult [j]=1,mp : d_pi_ff[j]    = F_D(Phi_tef[j])(/po_asp);
       do j=1,np
       < mult[i]=1,np : d_pi_H([i],j) = F_D(F_D(deta_tef[i])(/eta(j)))(/po_asp);
       >;
       do j=1,mp
       < mult[i]=1,np :  d_pi_Bb([i],j)= F_D(F_D(deta_tef[i])(/ff(j)))(/po_asp);
       >;
\end{verbatim}
%
Ces dérivées élémentaires sont combinées pour construire la sensibilité de aspha~:
%
\begin{verbatim}
+SELF,IF=DPi_Aspha. **********************************************
! ----------------------
! Calcul de d_pi_aspha .
! ----------------------------------------------------------------
!  D_p M = D_p[H] + D_f[H] D_p[f]
!                 + { D_f[B] D_p[f] + D_p[B] } (I-D)^ Ct
!                 + B { (I-D)^ D_p[Ct] - (I-D)^ D_p[D] (I-D)^ Ct };
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SEQ,dPi_HBbF.       ================> d_p ff et celles de H et Bb, puis :
       call matcopy(d_pi_H,d_pi_aspha,n,np,n,np);
       call matmlt(d_pi_Bb,Ctbuf,To,n,np,m,mp,n,np);
       call matadd(d_pi_aspha,To,d_pi_aspha,n,np,n,np);
+SEQ,dPi_CtD.          ==============> dérivations doubles de Ct et D, suivies de :
       call matmlt(d_pi_D,Ctbuf,T_tmp1,m,mp,m,mp,n,np);
       call matsub(d_pi_Ct,T_tmp1,T_tmp2,m,mp,n,np);
       call matmlt(un_moins_d,T_tmp2,T_tmp1,m,mp,m,mp,n,np);
       call matmlt(Bb,T_tmp1,To,n,np,m,mp,n,np);
       call matadd(d_pi_aspha,To,d_pi_aspha,n,np,n,np);
       call veczero(T_tmp2,n*m);
\end{verbatim}
%
Les calculs de la deuxième ligne (\ref{dpiaspha}) sont effectués en dérivant $H$ et $B_b$
par $\varphi_k$ dans une boucle sur $k$~:
%
\begin{verbatim}
   Do ko=1,m
   <
+SEQ,dffk_HBb. dérivation de H et Bb par rapport à ff(ko) :
     call scamat(d_pi_ff(ko),d_ffk_H,To,n,np,n,np);
     call matadd(d_pi_aspha,To,d_pi_aspha,n,np,n,np);
     call scamat(d_pi_ff(ko),d_ffk_Bb,T_tmp1,n,np,m,mp);
     call matadd(T_tmp2,T_tmp1,T_tmp2,n,np,m,mp);
   >;"End ko loop"

    call matmlt(T_tmp2,Ctbuf,To,n,np,m,mp,n,np);
    call matadd(d_pi_aspha,To,d_pi_aspha,n,np,n,np);
    if Zprint
    < title=' d_pi_aspha';
      CALL printmat(title,d_pi_aspha,n,np,n,np,0,0);
    >;
+SELF.
\end{verbatim}
%
Ces calculs étant effectués dans le corps de principal, on n'a plus qu'à y
ajouter son propre post-traitement, qui ici concerne ici comme illustration
le calcul de la sensibilité des valeurs
propres de $A^{\sharp}$ au paramètre, selon une idée de stb.
Les calculs sont les suivants~:
%
\begin{equation}
\scal{f_i}{e_i} d_{\varpi}\lambda_i = \bra{f_i}d_{\varpi} M \ket{e_i}
\end{equation}
%
avec $f_i$ et $e_i$ les couples de vecteurs propres à gauche et à droite de $M$.
On a une formule équivalente pour les valeurs singulières\index{Valeurs propres}
\index{Valeurs singulières}.

Ceci se fait dans {\bf ZSTEER}\footnote{noter que le passage
en double précision remplacera les {\bf call SGEEVX} en {\bf call DGEEVX}
automatiquement par activation du {\bf sel Double}.}. Noter au passage que la partie
la plus laborieuse est ici le rangement des éléments propores déivrés par LaPack !
Le codage est de nature pédagogique, et ne recherche pas à faire court.
%
\begin{verbatim}
+SELF,IF=dPi_aspha.

!!=============  CALCUL Valeurs et Vecteurs propres================
!            <fi|ei> d_pi Li = <fi| d_pi_aspha |ei>
!! author Stephane BLANCO
!!=================================================================
      call matcopy(aspha,a1a,n,n,n,n);
      call sgeevx( 'N', 'V', 'V', 'N', n, a1a, n,
                     WR, WI, VL, n, VR, n, ILO, IHI, SCALE, ABNRM,
                     RCONDE, RCONDV, WORK,LoWORK,IoWORK,INFO );

5100; format (501(1x,e12.5));
5101; format (501(1x,e14.7));
      do jjk=1,n
      < tmp_cmp(jjk)=CMPLX(WR(jjk),WI(jjk));
      >;
     if (Zoutbef) < write(70,5101)time,tmp_cmp;>;

      do ival=1,n
      <
       IF (WI(ival).eq.zero)
         <
          do i_coord=1,n
          <
           vec_propre_l_reel(i_coord,ival)=VL(i_coord,ival);
           vec_propre_l_im(i_coord,ival)=zero;
           vec_propre_r_reel(i_coord,ival)=VR(i_coord,ival);
           vec_propre_r_im(i_coord,ival)=zero;
          >;
         >else
         <
          IF (WI(ival).gt.zero)
          <
          do i_coord=1,n
           <
           vec_propre_l_reel(i_coord,ival+1)=VL(i_coord,ival);
           vec_propre_l_im(i_coord,ival+1)=VL(i_coord,ival+1);
           vec_propre_r_reel(i_coord,ival)=VR(i_coord,ival);
           vec_propre_r_im(i_coord,ival)=VR(i_coord,ival+1);
           >;
          >else
          <
          do i_coord=1,n
           <
           vec_propre_l_reel(i_coord,ival-1)=VL(i_coord,ival-1);
           vec_propre_l_im(i_coord,ival-1)=-VL(i_coord,ival);
           vec_propre_r_reel(i_coord,ival)=VR(i_coord,ival-1);
           vec_propre_r_im(i_coord,ival)=-VR(i_coord,ival);
           >;
          >;
         >; "end IF (WI(ival).eq.zero"
      >; "end ival"

      do ival=1,n
      <
       do i_coord=1,n
       <
         comp_vec_propre_l(i_coord,ival)=
         CMPLX(vec_propre_l_reel(i_coord,ival),vec_propre_l_im(i_coord,ival));
         comp_vec_propre_r(i_coord,ival)=
         CMPLX(vec_propre_r_reel(i_coord,ival),vec_propre_r_im(i_coord,ival));
       >
      >
!!=============  CALCUL sensibilites valeurs propres ================!!
     do i_val=1,n
      <
        do i_coord=1,n
        <
         comp_vec_l(i_coord)=comp_vec_propre_l(i_coord,i_val);
         comp_vec_r(i_coord)=comp_vec_propre_r(i_coord,i_val);
        >;
        call vecvec_comp(comp_vec_l,comp_vec_r,comp_lr,n,n);
        call matvec_comp(d_aspha_dpi,comp_vec_r,comp_tmp_vec,n,n,n,n);
        call vecvec_comp(comp_vec_l,comp_tmp_vec,comp_lMr,n,n);

      tmp_cmp(i_val)=comp_lMr/comp_lr;
      >;
     if (Zoutbef)
     < write(71,5101)time,tmp_cmp;
     >;

+SELF.
\end{verbatim}
%
% ========================================================
\subsection*{Autres calculs dépendants des Jacobiennes}
% ========================================================
%---------------------------------------------
\subsubsection{Intervention de mesures d'observables}
%---------------------------------------------
% ......................................................................
\subsubsection*{Fonction critère et rencontre avec les points de mesure}
% ......................................................................
% .................................................................................
\subsubsection*{Avance de la méthode de Newton pour la détermination de paramètres}
% .................................................................................
%..................................
\subsubsection*{Filtrage de Kalman}
%..................................
%  ............................................
% ========================================================
\subsection*{Processeurs peignes pour les modèles maillés}
% ========================================================
On illustre directement la méthode avec l'exemple suivant, tiré du modèle
d'urbanysme de Gus\index{fonctions peigne}. La fonction calcule l'écart quadratique de la variable {\bf util2(.)}
qui est un transfert\footnote{vari\_eta(...) est la même fonction pour une variable d'état.}.
On précise donc la variable du premier n\oe ud et le nombre de mailles~:
%
\begin{verbatim}
set_Probe
<  eqn: u_quad = vari_phi(n_node,util2(inode:1));
>;
\end{verbatim}
%
Une macro transforme cette fonction en point d'entrée Fortran~:
%
\begin{verbatim}
   ff(mp+1)=((VARI_NODE(n_node,ff(k2j_5(1)),j,FF(1),M_MULT)))
\end{verbatim}
%
avec les arguments Fortran permettant à la routine de calculer les positions
des variables dans le vecteur {\bf ff(.)}. Ceci se fait en comptant M\_MULT positions
(ou le double en double precision) des adresses de taleaux.

La dérivation demande d'une part une macro~:
%
\begin{verbatim}
&'(VARI_NODE(#,#,#,FF(1),M_MULT))(/FF(#))'='(DVARI_NODE(#1,#2,#4,FF(1),M_MULT))'
\end{verbatim}
%
avec appel au point d'entrée {\bf dvari\_node} qui, connaissant la vraiable de dérivation
{\bf ff(\# 4)} fournit la dérivée. Remarquons l'absence de règle de dérivation (/eta(.)) qui
résulte en une dérivée nulle.
\newpage
Pour les fonctions avec plus d'arguments, une autre procédure est utilisée en déclarant
explicitement le type des variables\index{macros Mortran}, exemple avec un calcul de moyenne pondérée~:
%
\begin{verbatim}
! Moyenne ponderee
&'AVRG_POIDS(#,PHI:#,ETA:#)'='AVRG_TWO(#1,#2,#3,j,FF(1),ETA(1),M_MULT,N_MULT)'
&'AVRG_POIDS(#,ETA:#,PHI:#)'='AVRG_TWO(#1,#3,#2,j,FF(1),ETA(1),M_MULT,N_MULT)'
&'AVRG_POIDS(#,PHI:#,PHI:#)'='AVRG_TWO(#1,#2,#3,j,FF(1),FF(1),M_MULT,M_MULT)'
&'AVRG_POIDS(#,ETA:#,ETA:#)'='AVRG_TWO(#1,#2,#3,j,ETA(1),ETA(1),N_MULT,N_MULT)'
! poids type eta * ff
&'(AVRG_POIDS(#,PHI:#,ETA:#))(/FF(#))'
='(DAVRG_TWO1(#1,#2,#3,#4,FF(1),ETA(1),M_MULT,N_MULT))'
&'(AVRG_POIDS(#,ETA:#,PHI:#))(/FF(#))'
='(DAVRG_TWO1(#1,#3,#2,#4,FF(1),ETA(1),M_MULT,N_MULT))'
&'(AVRG_POIDS(#,PHI:#,ETA:#))(/ETA(#))'
='(DAVRG_TWO2(#1,#2,#3,#4,FF(1),ETA(1),M_MULT,N_MULT))'
&'(AVRG_POIDS(#,ETA:#,PHI:#))(/ETA(#))'
='(DAVRG_TWO2(#1,#3,#2,#4,FF(1),ETA(1),M_MULT,N_MULT))'
! poids type ff * ff
&'(AVRG_POIDS(#,PHI:#,PHI:#))(/FF(#))'
='(DAVRG_TWO1(#1,#2,#3,#4,FF(1),FF(1),M_MULT,M_MULT)
+DAVRG_TWO2(#1,#2,#3,#4,FF(1),FF(1),M_MULT,M_MULT))'
&'(AVRG_POIDS(#,ETA:#,ETA:#))(/ETA(#))'
='(DAVRG_TWO1(#1,#2,#3,#4,ETA(1),ETA(1),N_MULT,N_MULT)
+DAVRG_TWO2(#1,#2,#3,#4,ETA(1),ETA(1),N_MULT,N_MULT))'
\end{verbatim}
%
La récupération des variables et leur dérivées suivent la même logique que dans
le cas plus simple.

Au 27/02/07, on a en stock les fonctions de calcul de moyenne simple ou pondérée,
de variance en accord\index{moyenne}\index{variance}, et un exemple de fonction ad'hoc calculant
un indice de Gini, avec dérivées par accroissements finis.
\newpage
% ------------------------------
\subsection{Calcul de l'adjoint}
% ------------------------------
%
On obtient le système adjoint à partir de~(\ref{sltc2})
%
\begin{equation}
\begin{bmatrix} -\partial_t v_{\varpi}\\ -\partial_t v \\ w \\ w_{\mu}\end{bmatrix} \quad  = \quad
\begin{bmatrix} 0 & H_{\varpi}^{\dagger}& C_{\varpi}  & c_1^{\dagger} \\
                0 & H^{\dagger}  &C& c_2^{\dagger}  \\
                0 & B_b^{\dagger} & D^{\dagger} & d_{\varphi\mu}^{\dagger} \\
                0 &0& 0 & 0\end{bmatrix} \quad 
\begin{bmatrix} v_{\varpi} \\ v  \\ w \\ w_\mu \end{bmatrix}
\quad + \quad\begin{bmatrix} l_{\varpi}\\ l_{\eta}  \\ l_{\varphi} \\ l_\mu\end{bmatrix}
\label{adj} 
\end{equation}
%
système dans lequel $v_\varpi(\tau)$ représente la sensibilité\index{sensibilité} (indirecte) de la fonction
de coût $J(T)$ à une modification des paramètres\index{paramètres} (à $t=\tau < T$), et ainsi des autres variables
adjointes. Les sources $l_x$
sont donnés par les dérivées partielles de l'intégrande {\bf cout$\_$l} (sensibilités directes).

Bien sùr, la sensibilité (indirecte) $w_\mu$ aux probes est découplée du système, puisque
les autres variables n'en dépendent pas par construction; par contre, $J$ pourra explicitement être fonction
de ces observables ($l_\mu =\partial_\mu\; \{(\mu(t_i)-M_i)^2\}$ dans le cas classique de la minimisation
de l'écart quadratique obs-mesures).

 Bref, on a en gros une situation où à nouveau le système couplé est le même
que sans paramètres ni probes, et des simplifications de calcul de ces composantes.
%
%.....................................
\subsubsection*{Trajectoire en arrière}
%.....................................
L'intégration à la TEF\index{TEF} du système~(\ref{adj}) donne
%
\begin{equation}
\begin{bmatrix}  -I-\frac{dt}{2} H^{\dagger}  &-\frac{dt}{2}C  \\
                 -B_b^{\dagger} & I-D^{\dagger} \end{bmatrix}
\begin{bmatrix}   \delta v \\ \delta w \end{bmatrix} \quad  = \quad
\begin{bmatrix}   \Gamma_v \\ \Omega_w \end{bmatrix} \quad  + \quad
\begin{bmatrix} l_{\eta}\delta t  \\ l_{\varphi}\end{bmatrix}
\label{adj2} 
\end{equation}
%
pour la partie couplée du système~(cf (\ref{A-fadjint}) du cahier Adjoint), avec
$\Omega_w=0$, et ou la contribution explicite est $\Gamma_v = \delta t[ H^{\dagger} v + C w ]$.

La première étape consiste à récupérer les Jacobiennes\index{Jacobienne}. Le principe d'économie consiste en un
compromis~: on sauve la trajectoire avec à la fois les vecteurs $\eta$ et $\varphi$ (mais pas
$\varphi_\varpi$), pour éviter la D-loop (A). On récupère également la commande lorsqu'une
loi est calculée en avant (B).

 Les Jacobiennes sont calculées classiquement (C), nous verrons
le calcul concernant la commande, les paramètres et les probes ultérieurement.

On passe des Jacobiennes aux matrices adjointes du TEF par transposition ({\bf call otransmat}, (D))
et avec $A^{\dagger}= -I -\frac{dt}{2}H^{\dagger}$ et c'est la $C=C^{\dagger\dagger}$ qui est
ici multipliée par $-\frac{dt}{2}$~:
%
\begin{verbatim}
! **************  Restitution fin pas retour **************

     time=time_sauve(istep);
     dt = time-timebef;
     <i=1,np; eta(i) = eta_sauve(i,istep); >;                 (A)
     <j=1,mp; ff(j) = ff_sauve(j,istep); >;
;
     if Zcommand
     < <i=1,np; ux_com(i) = ux_sauve(i,istep); >;             (B)
       <j=1,mp; uy_com(j) = uy_sauve(j,istep); >;
     >;

if ZPRINT< z_pr:timebef,time,dt,eta; z_pr:ff; z_pr:cout_J; >;

! Fonction de cout  (dt<0) par trapèzes
  coutstep = - (coutlbef+cout_l)*0.5*dt;
  cout_J = cout_J + coutstep;
;
! **************   Appel de la resolution  ****************
! ---------------- ( time backward ) ----------------------
+SEQ,D_Mx.
+SEQ,Matrice_ABCt.                                            (C)
! ====== calculs dependants de dt
! build adj a mx(init a I) from h : a = -I -dt/2*h
       call oscamat(-dt/2.,n,n+lp,h,np,hdt2m,np);
       call diagmat(-one,a,n,np+lp);
       call omatadd(n,n+lp,a,np,hdt2m,np,a,np);
! Build adj B mx from Bd  (identique)
       call omatcopy(n,m,Bb,np,B,np);

+SEQ,trans_abcd. et calculs vecteurs                          (D)

! -------------------------------------------------------------
  CALL oker("in" n,m,0,a_c,np+lp,ct_c,np+lp,b_c,mp,d_c,mp,
                 gam_adj,np,omeg_adj,mp,v_adj,np,w_adj,mp,    (E)
            lp,mobs,EtaNam,PhiNam,ParNam,zprint,
           "out" dv_adj,dphi, "in:" 0, "out:"gain,infoa,infores);
! -------------------------------------------------------------
! No error report as there were no error during the direct run
;
! incrementations
! -------------------------------------------------           (F)
  <i=1,np+lp;  v_adj(i) = v_adj(i) + dv_adj(i); >;
   coutlbef = cout_l;
   timebef = time;
   ZPRINT = mod(istep-1,modzprint).eq.0 .or.istep.eq.2
                                             .or.istep.eq.nstep;
   zprint=zprint.and.zprcall;
   Zout = mod(istep-1,modzout).eq.0;
   Zout = Zout .or. istep.eq.1 .or.(istep+1).eq.1;

  >;"end time loop"
! --------------------------------------------------- end back loop
\end{verbatim}
%
La résolution se fait par appel de {\bf oker} (E) avec échange des matrices
transposées par rapport aux calculs en avant, et on incrémente l'état adjoint
(ici avec {\bf nzp} nul, \cad  sans les extensions du direct, qui consisteraient en un calcul
de plusieurs fonctions de coût à la fois\footnote{extension possible mais non encore envisagée.}).

Le parcourt de la trajectoire s'effectue à rebrousse-temps avec {\bf dt} négatif~:
%
\begin{verbatim}
!---------------------------------------------------- back time loop
  DO istep = nstep,0,-1
  <
    IF ZPRINT< print9000; print *,'back_istep',istep,time; print9000; >;

! ---------------------------------------------------------
!  ********** calcul du w debut de istep avec D et B bef **
!             ( I-Dt) w = Bt v + l_y
! -----------------------------[sauf a l'entree]-----------
\end{verbatim}
%
avec pour le pas de temps initial utilisation des matrices $D^{\dagger}$ et $Bb^{\dagger}$
du dernier pas en avant -- alors que le calcul en avant devait commencer par calculer
la D-loop et $C^{\dagger}$. La source en provenance de l'intégrande de $J$ ne contribue qu'à
partir du pas de temps suivant, \cad le premier pas de l'intégration.

 Comme en avant, on commence par effectuer les calculs de la fin
d'un incrément (ou des valeurs initiales)~:
%
\begin{verbatim}
      CALL matvec(b_c,v_adj,b_c_v_adj,m,mp,n,np);

! -------------------------------------------------
! fin calcul w fin bef : ajout l_y et resolution
! -------------------------------------------------

    CALL matvec(un_moins_d_c,b_c_v_adj,w_adj,m,mp,m,mp);
;
    >;"end if istep courant"

! **************  Ecriture des valeurs a timebef **********
              o o o 

      if istep.eq.0 < goto :exitB:; >;
;
\end{verbatim}
%
%.............................................................
\subsubsection*{Sensibilité aux paramètres et vecteurs connus}
%.............................................................
%
Comme en avant, la sensibilité\index{sensibilité} aux paramètres de la {\bf free$\_$parameters} liste
est explicite\index{paramètres} une fois connue l'avance des autres~(cf~\ref{adj})~:
%
\begin{equation}
-\delta v_\varpi = \delta t \left[ H^{\dagger}_\varpi v + C w + l_\varpi \right] + \frac{dt}{2}H^{\dagger}_\varpi \delta v
+ \frac{dt}{2}C_\varpi w
\end{equation}
%
la seule difficulté réside ainsi dans la prise des bons signes et conventions, qui ont été
choisies pour l'appel à {\bf ker}, avec le résultat~:
%
\begin{verbatim}
! partie parametres : calcul direct de p_adj / v_adj et w_adj
! -------------------------------------------------------------
! d p_adj = A' dv_adj + C' dw_adj - gam_adj'   (a1a used as working space)
      if lp.gt.0
      < nw=np+1;
        call matvec(a_c(nw,1),dv_adj,dv_adj(nw),lp,np+lp,n,np);
        call matvec(ct_c(nw,1),dphi,a1a,lp,np+lp,m,mp);
        call vecadd(dv_adj(nw),a1a,dv_adj(nw),lp);
        call vecsub(dv_adj(nw),gam_adj(nw),dv_adj(nw),lp);
      >;
;
! -------------------------------------------------
! incrementations
! -------------------------------------------------
  <i=1,np+lp;  v_adj(i) = v_adj(i) + dv_adj(i); >;
\end{verbatim}
%
Enfin, les vecteurs connus (genre $\Gamma, \Omega$) sont auparavent déterminés
classiquement, ce qui fait que nous regroupons ici ces calculs~:
%
\begin{verbatim}
! calcul des vecteurs Gam_adj et omeg_adj
! Gam = dt ( l_x + Ht v + C w  )
! -----------------------------------------
     CALL matvec(H_c,v_adj,gam_adj1,n,np+lp,n,np);            (A)
     CALL matvec(ct_c,w_adj,gam_adj2,n,np+lp,m,mp);
     CALL vecadd(gam_adj1,gam_adj2,gam_adj,n);
! Omeg = 0
     <j=1,mp; omeg_adj(j) = 0.; >;

! ajout source l_x                                            (B)
     <i=1,n; gam_adj(i) = gam_adj(i) + f_d(cout_l)(/eta(i)); >

! partie parametres :
! gam_pi_adj = dt ( l'_x + Ht' v + C' w  )
! -----------------------------------------                   (C)
     if lp.ne.0
     < nw=np+1;
       CALL matvec(H_c(nw,1),v_adj,gam_adj1(nw),lp,np+lp,n,np);
       CALL matvec(ct_c(nw,1),w_adj,gam_adj2(nw),lp,np+lp,m,mp);
       CALL vecadd(gam_adj1(nw),gam_adj2(nw),gam_adj(nw),lp);
!    addition de l'x
       ;mult[i]=1,lp : gam_adj(np+[.i]) =
                       gam_adj(np+[.i]) + f_d(cout_l)(/eta(np+[.i]));
"---">;"end parametres ------------------"
;
     CALL scavec(dt,gam_adj,gam_adj,n+lp);                    (D)
\end{verbatim}
%
On a successivement~:
%
\begin{itemize}
\item la partie linéaire parallèle à $H\eta+B_b\varphi$ (A);
\item à laquelle on ajoute la contribution explicite de $J$~: $\partial_\eta J$ (B);
\item les paramètres ont leur incidence sur $v_{adj}$ et $w_{adj}$ données par
      l'extensions aux matrices (C);
\item et l'intervention du pas de temps (D);
\end{itemize}
%
%
%.........................................
\subsubsection*{Sensibilité à la commande}
%.........................................
%
La grande différence avec le direct est le calcul de sensibilité de $J$
à la commande\index{loi de commande}, qui est traitée à part, cf formule~(\ref{duuj2})~(la~\ref{A-duuj2} du cahier Adjoint),
alors que l'on aurait pu la considérer comme étant spécifiée par des composantes
de transferts -- leur différence ici étant le choix d'absence de matrice $D$
et la simplicité de leur application (une loi distincte par composante d'état
et de transfert)~:
%
\begin{eqnarray}
\nabla_{ux}J(\tau,T) &=& \bra{v(\tau)} B^u (\tau) \quad\,+\, \bra{l_{ux}(\eta(\tau),h^x(\tau))} \nonumber \\
\nabla_{uy}J(\tau,T) &=& \bra{w(\tau)} D^u (\tau) \quad+\, \bra{l_{uy}(\eta(\tau),h^y(\tau))}
\label{duuj2}
\end{eqnarray}
%
On a donc incidence de Jacobiennes\index{Jacobienne} du modèle par rapport aux commandes à
déterminer, d'où la programmation avec la partie commande via états (E)
suivie de commande via transferts (F)~:
%
\begin{verbatim}
!  calcul de Grad_ue J = <vB'| + l_ux
!  calcul de Grad_uf J = <wD'| + l_uy :
! -------------------------------------
      if Zcommand
      <                                                         (E)
        ;mult[i]=1,np: GradueJ[i] = f_d( v_adj[i]*(deta_tef[i])
                                        + cout_l )(/ux_com[i]);
        ;mult[j]=1,mp: GradufJ[j] = w_adj[j]*f_d(Phi_tef[j])(/uy_com[j])
                                   + f_d(cout_l)(/uy_com[j]);
      >;                                                        (F)
;
    >;"end if istep courant"
! **************  Ecriture des valeurs a timebef **********
\end{verbatim}
%
Il reste l'imputation de la fonction de coût\index{fonction de coût} terminale {\bf coût-Psi} que l'on
calcule donc en tête de l'adjoint (le code se trouve dans la séquence
{\bf Cout$\_$Psi})
%
\begin{verbatim}
!       Initialisation v(T) -> pour zinit-adjoint
! -----------------------------------------------------------------
!  |v(T)> = |d_eta cout_Psi> + Sum_j d_ff(j) cout_Psi |C(I-D')-1>_j
! -----------------------------------------------------------------
! - partie dependance de eta
  <i=1,n; v_adj(i) = f_d(cout_Psi)(/eta(i)) >;                  (G)

! - partie dependance de ff (Ctbuf=CtDm1(T) est disponible)
  do j=1,m
  < aaa = f_d(cout_Psi)(/ff(j));
    if abs(aaa).ge.1.E-06
    < <i=1,n; v_adj(i) = v_adj(i) + aaa*Ctbuf(j,i) >;           (H)
    >;
  >;
  if abs(aaa).ge.1.E-06 <Print*,' *** cout_Psi dependant de ff(T) pulse ***';>;

! - partie parametres
  if lp.gt.0                                                    (I)
  < ;mult[i]=1,lp : v_adj(np+[.i]) = f_d(cout_Psi)(/eta(np+[.i]));
  >;
\end{verbatim}
%
On remarque en plus de la contribution de l'état final (G) une
éventuelle contribution de la valeur d'un tranfert final (H) -- considéré
alors comme générant un pulse sur l'état comme il est à présent admis
(cf cahier Effluents).

Enfin des paramètres\index{paramètres} pouvant apparaître explicitement dans le finale,
on a ce terme (I) -- dans lequel rappelons que {\bf eta(np+[.i])} sera
remplacé tour à tour ({\bf mult[i]}) par Mortran\index{macros Mortran} par chaque
symbole des paramètres de la {\bf Free} liste.
%
%......................................
\subsubsection*{Traitement des probes}
%......................................
En suivant cette fois l'ordre séquentiel\index{observables}, on calcule d'abord $w_\mu$~:
%
\begin{verbatim}}
!  ********** calcul du w_mu debut de istep avec D et B bef **
!              w_mu =  l_mu
! ---------------------------------------------------------
    if mobs.gt.0
    < if istep .lt. nstep
      < <j=1,mobs; w_adj(mp+j) = f_d(cout_l)(/ff(mp+j)); >;
      >else
      < <j=1,mobs; w_adj(mp+j) = f_d(cout_Psi)(/ff(mp+j)); >;
      >;
    >;
\end{verbatim}
%
qui fournit un term additionnel (J) au calcul de $w$ à la fin
du pas de temps~:
%
\begin{verbatim}}
! ---------------------------------------------------------
!  ********** calcul du w debut de istep avec D et B bef **
!             ( I-Dt) w = Bt v + d'_ffmu w_mu + l_y
!                                  (J)

      CALL matvec(b_c,v_adj,b_c_v_adj,m,mp,n,np);
      CALL matvec(d_c(1,mp+1),w_adj(mp+1),d_c_w_mu,m,mp,mobs,mobs);  (J)

    if mobs.gt.0
    < <j=1,mp; b_c_v_adj(j) = b_c_v_adj(j) + d_c_w_mu(j); >;         (J)
    >;

    if istep .lt. nstep
    < <j=1,mp; b_c_v_adj(j) = b_c_v_adj(j) + f_d(cout_l)(/ff(j)); >;
    >else
    < <j=1,mp; b_c_v_adj(j) = b_c_v_adj(j) + f_d(cout_Psi)(/ff(j)); >;
    >;

    CALL matvec(un_moins_d_c,b_c_v_adj,w_adj,m,mp,m,mp);
\end{verbatim}
%
On a aussi des contributions au calcul des $\Gamma$ (K)~:
%
\begin{verbatim}}
! ---------------------------------------------------
! calcul des vecteurs Gam_adj (et omeg_adj  nul)
! Gam_adj    = dt ( l_x + Ht v + C w + c_2 w_mu )
! Gam_pi_adj = dt ( l_pi + Ht' v + C' w + c_1 w_pi )
! ---------------------------------------------------

     CALL matvec(H_c,v_adj,gam_adj1,n+lp,np+lp,n,np);
     CALL matvec(ct_c,w_adj,gam_adj2,n+lp,np+lp,m+mobs,mp+mpbs);     (K)
     CALL vecadd(gam_adj1,gam_adj2,gam_adj,n);

     <j=1,mp; omeg_adj(j) = 0.; >;

! ajout sources l_x, l_pi
     <i=1,n; gam_adj(i) = gam_adj(i) + f_d(cout_l)(/eta(i)); >
     if lp.gt.0
     < ;mult[i]=1,lp : gam_adj(np+[.i]) =                            (K)
                       gam_adj(np+[.i]) + f_d(cout_l)(/eta(np+[.i]));
     >;
\end{verbatim}
%

%......................................
\subsubsection*{Calcul du Hamiltonien}
%......................................
%
Toujours dans le cahiers Adjoint, on a une formulation de la fonction dîte
du Hamiltonien\index{Hamiltonien} qui doit être calculée à chaque instant~:
%
\begin{eqnarray}
H(\eta,\varphi,v,w,h^x,h^y) = \scal{v\,}{ g(\eta,\varphi ,h^x)} 
+\scal{w\,}{ f(\eta,\varphi ,h^y)} + l(\eta,\varphi,h)
\end{eqnarray}
%
son intérêt rappelons le concernant son extrêmalité en $u$ considérée comme
simple variavle d'une part, et de l'autre lors de la recherche d'optimalité
du temps terminal ($H(T)$ devant alors s'annuler.
%
\begin{verbatim}
!  ********** Fonctions dependant de <v| et <w| (t*********
! -----------------------------[sauf a l'entree]-----------
    if istep .lt. nstep
    <
! -********************************
!   Calcul du Hamiltonien courant :
! -********************************
!  H(eta,phi,v,w,ux,uy) =
!                   <v| g(eta,ff,ux)> + <w| f(eta,ff,uy) + l(eta,phi,u)
      Hamilton = 0.;
      mult [i]=1,np : Hamilton = Hamilton + v_adj[i]*(deta_tef[i]);
      mult [j]=1,mp : Hamilton = Hamilton + w_adj[j]*(Phi_tef[j]);
      Hamilton = Hamilton + cout_l;
\end{verbatim}
%
la variable {\bf cout$\_$l} étant l'intégrande de $J$ donnée dans {\bf Zinit}
par une {\bf f$\_$set}.
%..................................................
\subsubsection*{Variables et fichiers de l'adjoint}
%..................................................
%
%
\begin{flushleft}
\begin{tabular}{||p{1.7cm}||*{4}{c|}|}
\hline
{\it nom} & dim 1 &  & {\it signification} & data file\\
\hline
\hline
v$\_$adj  & np & 1 & vecteur adjoint d'état & vadj.data\\
\hline
v$\_$adj  & +lp & 1 & vecteur adjoint paramètres & gradpj.data\\
\hline
dv$\_$adj  & np+lp & 1 & avance  & no data\\
\hline
GradueJ  & np & 1 & sensibilités à la commande d'état  & graduxj.data\\
\hline
w$\_$adj  & mp & 1 & vecteur adjoint de transfert & wadj.data\\
\hline
w$\_$adj(+)  & +mobs & 1 & partie observables & gradmuj.data\\
\hline
GradufJ  & mp & 1 & sensibilités à la commande de transferts  & graduyj.data\\
\hline
\hline
\end{tabular}
\end{flushleft}
%
%
% -----------------------------------------------
\subsection{Recherche d'extremum de fonctionnelle}
% -----------------------------------------------
\clearpage
% ------------------------------------------
\section{de quelques analyses dynamiques}
% ------------------------------------------
\subsection{Vecteurs de Floquet}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L'idée serait ici, un brouillon avant codage, d'avoir une cohérence entre
le calcul des Floquets et des gains et réponses du système périodique.
Les algorithmes suivants sont pensés pour leur passage à MK2 et à sa navette.
%
\paragraph{recherche d'agorithme efficace}
%-----------------------------------------
On part des remarques des {\bf Effluents}, on prolonge jusqu'à l'explicite\index{Floquet!(vecteurs de)}.
L'idée est d'éviter le plus possible les variables complexes, en profitant
de leur rangement par couple de colonnes $\{\ket{x_r}, \ket{x_i}\}$ pour tout
ce qui leur applique des produits de matrices réelles. On a alors une partie
du calcul qui avance comme les sensibilités, et on post-traite pour
appliquer les corrections dûes aux produits complexes. {\it Dans le exos de Céline et Pablo,
on s'embêtait pas et on calculait $\Psi^F(t)=\Phi(t,0)Zexp(-\Lambda t)$ en complexes; mais
d'une part, on avance les deux parties de chaque paire, l'imag nulle des réels etc.
Par ailleurs, du fait de la précision numérique, on obtenait des V.F. erratiques
sur lesquels on fermait les yeux --- garbage out...}

On suppose que $\Psi^r(t)=\Psi^F(t)$ --- \cad formées des mêmes colonnes --- et
$\Psi^r(t+\tau)=\Phi(t+\tau,t)\Psi^r(t)$, ce calcul de $\delta\Psi^r$ séparant
les deux étant effectué dans $\bf oker$. On a alors
%
\begin{eqnarray}
\text{\`a }t:\;\Phi^r&=&\Phi^F\rightarrow\text{MK}\rightarrow\Phi^r(t+\tau)=\Phi(t+\tau,t)\Phi^F\nonumber\\
&&\text{d'o\`u :}\;\delta\Phi^r\;;\nonumber\\
\Phi^F(t+\tau) &=& \Phi(t+\tau,t)\Phi(t,0)Ze^{-\Lambda(t+\tau)} \nonumber\\
 &=& \Phi(t+\tau,t)\Phi^F(t)e^{-\Lambda\tau} \nonumber\\
 &=& \Phi^r(t)\;e^{-\Lambda\tau} \nonumber
\label{psiminik2}\end{eqnarray}
%
d'où il suit que 
%
\begin{equation}
\boxed{
\Phi^F=\Phi^r \exp(-\Lambda\tau).
}
\label{dFlo}\end{equation}
%
On a ainsi une correction avec matrice diagonale de coeffs complexes ou non.
On a autant de colonnes à avancer que celles de $\Phi_t$, d'où l'idee de l'initialiser
avec les bonnes parties de $Z$, l'avancer classiquement et appliquer~(\ref{psiminik2})
vers ZSTEER.
Il y a ainsi à construire une gestion des colonnes et corrections, dont la mise au
point suit ce paragraphe.
%
\paragraph{Algorithme de construction des pointeurs et tableau des v.p.}
%----------------------------------------------------------------------
%
On part des sorties de LaPack au moment du stockage de $Z$ à la première période.
On suppose pour fixer les idées $np=7$, avec une suite de vp réelle, complexe, réelle, complexe,
réelle. Comme remarque supplémentaire, on doit tenir compte de l'existence possible
d'un multiplicateur négatif, correspondant à un V.F. réel corrigé par un facteur de type
$\cos(\omega t)exp(-\lambda t$, c'est pourquoi on conserve sous forme complexe
tous les exposants de Floquets pour intégredr ce cas à la correction du cas réel.
%
\begin{verbatim}
+SELF,IF=Flo_vec.
      complex Flo_Z(np,np),Flo_sigma(np),ovpF(np),wox;
      real Flo_Period,aor,bor,orvec(np),oivec(np),ww,dtby2;
      integer ko_colF(2*np);

+SELF,IF=Flo_vec.
! --------------------------------------------------------
!                r + - r + - r    (rgt en sortie de Lapack)
!   exemple WR : * * * * * * *    n=7
!           WI : 0 *-* 0 *-* 0        => =2n (Lapack donne toutes les paires)
!                1 2 3 4 5 6 7
!
!  rgt des v.p r+is                          Il doit etre de meme dim que Phi_F
!        vpF  : r 0 r i 0 r i r 0   = 2n    mais en n*cmplx
!  ptr colonnes de Phi^F et yes/no correction   ************
!       kcolF : 1,no 2,yes 4,no 5,yes 7,no  <= 2n  (ici 10)
! ! --------------------------------------------------------  
! ozvf dit r1 pour colonne relle, R2,I2 pour complexe etc
character*2 ozvf(np);
;
! calcul des exposants de Floquet
! -------------------------------
OPEN (50, File='flovp.init', status='OLD');
      READ (50,1009) Flo_Period, (Flo_sigma(i),i=1,np);
      CLOSE(50);
;
! et initialisation de Phi_t avec Z
! ---------------------------------
OPEN (50, File='floVP.init', status='OLD');
      READ (50,1009) ww, ((Flo_Z(i,j),i=1,np),j=1,np);
      CLOSE(50);
;
  <i=1,2*n; ko_colF(i)=0; >;
   i=1;"ptr sur une premiere v.p. des paires"
   k=1;"info VP (une info par paire)"
  While i.le.n
  < ovpF(i)= log(Flo_Sigma(i))/Flo_Period;

    if Imag(Flo_Sigma(i)).eq.0
    < ko_colF(k)=i; ko_colF(k+1)=-1"no";
      ozvf(i)='r'//charit2(i);
      <j=1,n;Phi_t(j,i)=Real(Flo_Z(j,i));>;
     [i:+1];
    >else
    < ko_colF(k)=i; ko_colF(k+1)=1"yes";
      ozvf(i)='R'//charit2(i); ozvf(i+1)='I'//charit2(i);
      <j=1,n;Phi_t(j,i)=Real(Flo_Z(j,i));Phi_t(j,i+1)=Imag(Flo_Z(j,i));>;
     [i:+2];
    >;
   [k:+2]
  >;"EndDo"
  print9000;
  Z_pr/Floquet/: ((j,Phi_t(i,j),i=1,np),j=1,np);
  Z_pr: Flo_period, Flo_Sigma;
  Z_pr/mult de Floquet/: ovpF;

\end{verbatim}
%
\paragraph{Calculs correctifs}
%-----------------------------
%
On ne fait évidemment avancer qu'un seul Floquet complexe de chaque paire.
À faire vers ZStep (on utilise les dPhi et Phi-(t))~:
%
\begin{verbatim}
!##################
+SELF, IF=FLO_VEC.
!##################
!-! Epreuve d'algorithme a rentrer dans le code standard
!-! -------------------------------------------------------------
!-! Phi_t = Phi_t exp(-Lambda.dt) ou Ph-It(t) exp(-Lambda. t)
!-! à choisir en fonction des multiplicateurs.
!-! -------------------------------------------------------------
!-  Do io=1,2*n,2
!-  <
!-   ko=ko_colF(io); "ptr 1st col"
!-    if ko_colF(io+1).eq.1                     "cas complexe"
!-    <
!-! ------------------------------------------------------------
!-!    Calcul du facteur correctif
!-!    ---------------------------
!-     wox = cmplx(ovpF(ko),ovpF(ko+1));
!-     wox = exp(-wox*dt);
!-;
!-!    Correction du cas complexe
!-!    --------------------------
!-     aor = real(wox); bor = imag(wox);
!-     call omatcopy(n,1,Phi_t(1,ko),np,orvec,np);
!-     call omatcopy(n,1,Phi_t(1,ko+1),np,oivec,np);
!-     Do jo=1,np
!-     < Phi_t(jo,ko)   = aor*orvec(jo) - bor*oivec(jo);
!-       Phi_t(jo,ko+1) = bor*orvec(jo) + aor*oivec(jo);
!-     >;"EndDo"
!-!    ----------------------------
!-!     Correction du cas reel ----
!-!    ----------------------------
!-    >else                                        "cas reel"
!-    < aor = exp(-dt*ovpF(ko));
!-      Do jo=1,np
!-      < Phi_t(jo,ko) = aor*Phi_t(jo,ko);
!-      >;"EndDo"
!-;
!-    >;"EndIf"
!-;
!-  >;"end Do"
! Epreuve d'algorithme a rentrer dans le code standard
! -------------------------------------------------------------
! Phi_t = Phi_t exp(-Lambda.dt) ou Ph-It(t) exp(-Lambda. t)
! à choisir en fonction des multiplicateurs.
! Version integrant un multiplicateur reel negatif
! -------------------------------------------------------------
  Do io=1,2*n,2
  <
   ko=ko_colF(io); "ptr 1st col"
!    Calcul du facteur correctif
!    ---------------------------
   wox = exp(-ovpF(ko)*dt);
   aor = real(wox); bor = imag(wox);
    if ko_colF(io+1).eq.1                     "cas complexe"
    <
!    Correction du cas complexe
!    --------------------------
     call omatcopy(n,1,Phi_t(1,ko),np,orvec,np);
     call omatcopy(n,1,Phi_t(1,ko+1),np,oivec,np);
     Do jo=1,np
     < Phi_t(jo,ko)   = aor*orvec(jo) - bor*oivec(jo);
       Phi_t(jo,ko+1) = bor*orvec(jo) + aor*oivec(jo);
     >;"EndDo"
!    ----------------------------
!     Correction du cas reel ----
!    ----------------------------
    >elseif ko_colF(io+1).eq.-1                       "cas reel"
    < 
      Do jo=1,np
      < Phi_t(jo,ko) = aor*Phi_t(jo,ko);
      >;"EndDo"
;
    >;"EndIf"
;
  >;"end Do"
\end{verbatim}
%
On recalcule les coefficients à partir des vp Floquet pour l'avance des vecteurs, parce que l'incrément
de temps peut changer. On devrait alors avoir successivement
%
\begin{enumerate}
\item à la détermination des Floquets : stockage Z et $\sigma$ en complexe (dimension np);
\item init~: calcul des multiplicateurs de Floquet; et rangement kCol(1:<2n) et vpF (dim np);
\item initialisation de Phi\_t avec Z ;
\item à chaque fin du pas de calcul : coefficients dépendant du dt et corrections d'avance complexe
ou réelle;
\item orvec et oivec sont deux vecteurs colonnes de dim np;
\end{enumerate}
%
avec cette solution, rien à changer dans les calculs d'avance. Essais validés avec ClimSI, reste
à faire avec des Floquet complexes.
%
\paragraph{Analyse de Floquet de l'Effet-en-retour}
%-----------------------------------------------
%
ici pour l'instant~:

pour la partie avance d'une période
%
\begin{verbatim}
! Calcul de Phi_t et Psi_t en avant (nxp=np)
! ------------------------------------------
 Sensy_to_var
 <         var: T_ocean_eta   , pert: init;
           var: T_tropo_eta   , pert: init;
           var: T_strato_eta  , pert: init;
           var: Q_tropo_eta   , pert: init;
           var: CO2_conc_eta   , pert: init;
! Gain de rétroaction
           var: T_WV_phi      , pert: step_eff;
 >;
! moderateur sensib a xCO2 :
!***  Phi_t(5,5)=3.30d-04; "XXX"
\end{dverbatim}
%
pour la partie avance d'une période (II) varrho_m
%
\begin{verbatim}
! Calcul de Phi_t et Psi_t en avant (nxp=np)
! ------------------------------------------
 Sensy_to_Eta
 < "avance de sigma^F : pour chaque rep-eff"
!                il faut en prendre une et initialiser la colonne avec sig_T
           var: T_ocean_eta   , pert: init;   on a un pb de notation alors ?
           var: T_tropo_eta   , pert: init;
           var: T_strato_eta  , pert: init;
           var: Q_tropo_eta   , pert: init;
           var: CO2_conc_eta  , pert: init;  celle 
 >;
! Rep-eff rétroaction
 Sensy_to_Phi
 <           var: T_WV_phi      , pert: step_eff;
 >;
! moderateur sensib a xCO2 :
!***  Phi_t(5,5)=3.30d-04; "XXX" à présent dans minikq.car :
       step_psi(6)=3.30d-04;

! copie de sigma(T) dans les colonnes de Phi_t
! copie de Z dans Phi_t nxn
\end{verbatim}
%
Il faudrait une notation pour reserver la place des Floquet --- et pour $\sigma^F$ ? 

En profiter pour ajouter une source symbolique attachée à une rep-eff ? Il faudrait
lister les sources dans une ko\_p\_source(k)=nom de la fonction définie par fun\_set source=f(t),
calculer ces sources par une mult[]  et les attacher à step\_psi(ko)=source;

{\bf Finalement, hors flo-resp, on a nxp total, nyp=nyp2=0 ou le nbre des phi sensy; et pour
flo-resp, nxp total, nyp2 de phi-sensy, nyp=2*nyp (on doit avoir np=nxp-nyp).}

Les macros 102 de sensy-to-var corrigées pour Floquet-effets-en-retour~:
%
\begin{verbatim}
! ----------------------------------------------------------------------------
! Sensitivity definition
! Le principe des defauts est que l'instruction
!  if (mod(istep+1,ko_p_mod(1)).eq.ko_p_ini(1)+1) Zko_p(1)=.true.;
! done un toujours .false. [et mod(x,-1) donne tjrs 0 ]
&'$PERT:INIT;'='ko_p_typ(#C-1)=1;' &'$PERT:STEP_EFF;'='ko_p_typ(#C-1)=4;'
&'$PERT:STEP_RESP;'='ko_p_typ(#C-1)=3;' &'$PERT:PULSE;'='ko_p_typ(#C-1)=2;'
&'$AT#;'='ko_p_ini(#C-1)=#1;ko_p_mod(#C-1)=-1;'
&'$AT#EVERY#;'='ko_p_ini(#C-1)=#1;ko_p_mod(#C-1)=#2;'
&'ETA(#)$VT$'='ko_p_tvar(#C)=1;ko_p_kvar(#C)=#1;'
&'FF(#)$VT$'=' ko_p_tvar(#C)=2;ko_p_kvar(#C)=#1;'
&'$PRT:INIT;'='ko_p_typ(nxp-nyp+#C-1)=1;'
&'$PRT:STEP_EFF;'='ko_p_typ(nxp-nyp+#C-1)=4;'
&'$PRT:STEP_RESP;'='ko_p_typ(nxp-nyp+#C-1)=3;'
&'$PRT:PULSE;'='ko_p_typ(nxp-nyp+#C-1)=2;'
&'$AS#;'='ko_p_ini(nxp-nyp+#C-1)=#1;ko_p_mod(nxp-nyp+#C-1)=-1;'
&'$AS#EVERY#;'='ko_p_ini(nxp-nyp+#C-1)=#1;ko_p_mod(nxp-nyp+#C-1)=#2;'
&'ETA(#)$VTF$'='ko_p_tvar(nxp-nyp+#C)=1;ko_p_kvar(nxp-nyp+#C)=#1;'
&'FF(#)$VTF$'=' ko_p_tvar(nxp-nyp+#C)=2;ko_p_kvar(nxp-nyp+#C)=#1;'
! on pourrait compter les etats et transferts et remplir nyp (a voir avec l unsecte)
&';SENSY_TO_VAR<'='Call veczero(PHI_T(1,1),np*nxp);WRITE(*,*);
WRITE(*,*)'' --------------- Informing on Sensitivity to variables ----------'';
;WRITE(*,*)''  Var-name,              Perturbation-type,    column in Phit Mx'';
#L8#S1#L #S6#N; &''VAR:##,''=''##U6##S6;##1$VT$;,##1$''
 &'',##$PERT:##;''='' ko_p_ini(##C)=1;ko_p_mod(##C)=-1;
 WRITE(*,''''(2A24,I3)'''')''''##1'''',''''##2'''',##U6##C;$PERT:##2;##A1##S6;''
 &'',##$PERT:## AT''=''
 WRITE(*,''''(2A24,I3)'''')''''##1'''',''''##2'''',##U6##C;$PERT:##2;##A1##S6;
$AT'' &''$8$''=''##U6;##N;
 WRITE(*,*)'''' ----------------------------------------------------'''';
 IF(nxp.NE.##C-1)
WRITE(*,*)''''**WARNING : Mismatched definition / parameter nxp='''' ,nxp,##C-1;
 &''''####K$8$'''' &''''####KVAR:####,'''' &''''####K,####$PERT:####;''''
jo=0;ko=0;
DO io=1,nxp;IF (ko_p_typ(io).eq.1.and.ko_p_tvar(io).eq.1.and.ko_p_ini(io).eq.1)
THEN;Phi_t(ko_p_kvar(io),io)=1.;ko=ko+1;ENDIF;
IF (ko_p_tvar(io).eq.1.and.ko_p_typ(io).ne.1)THEN;
write(*,*)''''***ERROR : Cannot set Pert_type'''',ko_p_typ(io),
'''' on State variable:'''',EtaNam(ko_p_kvar(io));jo=1;ENDIF;
IF (ko_p_tvar(io).ne.1.and.ko_p_typ(io).eq.1)THEN;
write(*,*)''''***ERROR : Cannot set C.I. sensy type for Transfer variable:''''
,PhiNam(ko_p_kvar(io));jo=1;ENDIF;ENDDO;
+SELF,IF=Flo_resp.   --------------------------------------------------------
IF(ko.NE.N)THEN;WRITE(*,*)''''***ERROR all State Variable required in sensy
to Var block'''';jo=1;
! copie de Z dans Phit et de sigmaT dans son double
;####* === Initialisations: Z et sigma_T dans PhiT(1:n,1:2n) ============;;
! *** A FAIRE ***
+SELF.
IF (jo.ne.0)write(loerr,*)'''' *** ERROR in ZINIT Sensy_to_var'''',char(7);'''
! ------------------- new commands added ---------------------------------------
&';SENSY_TO_ETA<'=';WRITE(*,*);
WRITE(*,*)'' --------------- Informing on Sensitivity to Eta-variables ------'';
;WRITE(*,*)''  Var-name,              Perturbation-type,    column in Phit Mx'';
#L8#S1#L #S6#N; &''VAR:##,''=''##U6##S6;##1$VT$;,##1$''
 &'',##$PERT:##;''='' ko_p_ini(##C)=1;ko_p_mod(##C)=-1;
 WRITE(*,''''(2A24,I3)'''')''''##1'''',''''##2'''',##U6##C;$PERT:##2;##A1##S6;''
 &'',##$PERT:## AT''=''
 WRITE(*,''''(2A24,I3)'''')''''##1'''',''''##2'''',##U6##C;$PERT:##2;##A1##S6;
$AT'' &''$8$''=''##U6;##N;
 WRITE(*,*)'''' ----------------------------------------------------'''';
 &''''####K$8$'''' &''''####KVAR:####,'''' &''''####K,####$PERT:####;'''' jo=0;
Call veczero(PHI_T(1,1),np*(nxp-nyp));                                             ****
DO io=1,##C-1;IF (ko_p_typ(io).eq.1.and.ko_p_tvar(io).eq.1.and.ko_p_ini(io)
.eq.1)Phi_t(ko_p_kvar(io),io)=1.;
IF (ko_p_tvar(io).eq.1.and.ko_p_typ(io).ne.1)THEN;
write(*,*)''''***ERROR Cannot set Pert_type'''',ko_p_typ(io),
'''' on State variable:'''',EtaNam(ko_p_kvar(io));jo=1;ENDIF;
IF (ko_p_tvar(io).ne.1)THEN;write(*,*)''''***ERROR Transfer variable in sensy_to
_Eta block:'''';jo=1;ENDIF;ENDDO;
+SELF,IF=Flo_resp.   --------------------------------------------------------
IF(##C-1.NE.NP)THEN;WRITE(*,*)''''***ERROR all State Variable required in sensy
 to Eta block'''';jo=1;
! copie de Z dans Phit et de sigmaT dans son double
;####* === Initialiser : Z dans PhiT(1:n,1:n) ============;;
! *** A FAIRE ***
+SELF.
IF (jo.ne.0)write(loerr,*)'''' *** ERROR in ZINIT Sensy_to_Eta'''',char(7);'''
! ------------------------------------------------------------------------------
&';SENSY_TO_PHI<'=';WRITE(*,*);
WRITE(*,*)'' --------------- Informing on Sensitivity to Phi-variables ------'';
;WRITE(*,*)''  Var-name,              Perturbation-type,    column in Phit Mx'';
#L8#S1#L #S6#N; &''VAR:##,''=''##U6##S6;##1$VTF$;,##1$''
 &'',##$PERT:##;''='' ko_p_ini(nxp-nyp+##C)=1;ko_p_mod(nxp-nyp+##C)=-1;
 WRITE(*,''''(2A24,I3)'''')''''##1'''',''''##2'''',nxp-nyp+##U6##C;$PRT:##2;
##A1##S6;''
 &'',##$PERT:## AT''=''
 WRITE(*,''''(2A24,I3)'''')''''##1'''',''''##2'''',nxp-nyp+##U6##C;$PRT:##2;
##A1##S6;$AS'' &''$8$''=''##U6;##N;
 WRITE(*,*)'''' ----------------------------------------------------'''';
 &''''####K$8$'''' &''''####KVAR:####,'''' &''''####K,####$PRT:####;'''' jo=0;
Call veczero(PHI_T(1,nxp-nyp+1),np*nyp);
DO io=nxp-nyp+1,nxp-nyp2;IF (ko_p_tvar(io).eq.1)THEN;
write(*,*)''''***ERROR Eta variable in sensy_to_Phi block: ''''
,EtaNam(ko_p_kvar(io));jo=1;ENDIF;
IF (ko_p_tvar(io).ne.1.and.ko_p_typ(io).eq.1)THEN;
write(*,*)''''***ERROR Cannot set C.I. pert type for Transfer variable:''''
,PhiNam(ko_p_kvar(io));jo=1;ENDIF;ENDDO;
IF (jo.ne.0)write(loerr,*)''''*** ERROR in ZINIT Sensy_to_Phi'''',char(7);'''
! ----------------------------------------------------------------------------
;


  --------------- Informing on Sensitivity to variables ----------
   Var-name,              Perturbation-type,    column in Phit Mx
             T_ocean_eta                    init  1
             T_tropo_eta                    init  2
             Q_tropo_eta                    init  3
                T_WV_phi               step_resp  4
            T_Planck_phi               step_resp  5
              T_Evap_phi               step_resp  6
              T_sens_phi               step_resp  7
            T_FlxCMO_phi               step_resp  8
                xCO2_phi               step_resp  9
  ----------------------------------------------------
\end{verbatim}
%
On décide que pour le calcul des vecterus de Floquet, on place $Z$ dans $\Phi_t(t=0)$ \cad pas
de chgt à sensy\_to\_var$<>$.
On a déjà les rep-eff qui y sont demandées, de taille nyp, on double cette dimension pour y
adjoindre la partie $\overline{\sigma(t)}$ derrière~: nxp=[np]+2*nyp2, avec nyp=2*nyp2.
dans ce cas de figure, le bloc sensy\_to\_var$<>$ ne serait pas à modifier, on ajouterait seulement
les initialisations, et on ajouterait aux dimensions nyp et nyp2.
pSteer devra se limiter à nxp-nyp2.

Macros dodimetaphi~:
%
\begin{verbatim}
! --------- al1 15/09/09 chgt #L -------
                           | ici pound sign pas codé par Latex latin1
&';SENSY_TO_VAR<'='#LS#S1#L#S6 &''##K;IF NXP.gt.0<''
 &'';VAR:##,##;''=''##U6##A1##S6;''
 &'';$S$;''='' &''''####K;$S$;'''' &''''####K;VAR:####,####;''''
+SELF,IF=-Flo_resp.   --------------------------------------------------------
 &'''';parameter(nxp=0,''''='''';parameter(nxp=##U6##C,'''';'';'
+SELF,IF=Flo_resp.    --------------------------------------------------------
 &'''';parameter(nxp=0,nyp=0''''=
'''';parameter(nxp=##U6##C+nyp,nyp=##C-np'''';'';'
+SELF.
&';SENSY_TO_ETA<'='#LS#S1#L#S6 &''##K;IF NXP.gt.0<''
 &'';VAR:##,##;''=''##U6##A1##S6;''
 &'';$S$;''='' &''''####K;$S$;'''' &''''####K;VAR:####,####;''''
+SELF,IF=-Flo_resp.   --------------------------------------------------------
 &'''';parameter(nxp=####,''''='''' &''''''''########K;parameter(nxp=########,''
'''''';parameter(nxp=####1+##U6##C,'''';'';'
+SELF,IF=Flo_resp.    --------------------------------------------------------
 &'''';parameter(nxp=0,''''='''';parameter(nxp=##U6##C+nyp,''''
;character*2 ozvf(np);'';'
+SELF.
&';SENSY_TO_PHI<'='#LS#S1#L#S6 &''##K;IF NXP.gt.0<''
 &'';VAR:##,##;''=''##U6##A1##S6;''
 &'';$S$;''='' &''''####K;$S$;'''' &''''####K;VAR:####,####;''''
+SELF,IF=-Flo_resp.   --------------------------------------------------------
 &'''';parameter(nxp=####,nyp=0''''=''''
;parameter(nxp=####1+##U6##C,nyp=##C'''';'';'
+SELF,IF=Flo_resp.    --------------------------------------------------------
 &'''';parameter(nxp=####,nyp=0,nyp2=0''''=''''
;parameter(nxp=####1,nyp=2*##U6##C,nyp2=##C'''';'';'
+SELF.
! ----------------------------------------------------------------------------
Il vaut mieux calculer eta et ff indépendamment *** cf ci-dessus
et en accord :  (il n'y a que nxp-nyp2 variables sensibles) mais nxp colonnes
(mais les modifs sont OK pour sensy-to-var quand même)

! extension pour Phi_t et Psi_t (propagateurs TEF)
   if nxp.gt.0
   < OPEN(loPhit,FILE='sens'//osuffix,STATUS='UNKNOWN');
      write(loPhit,1011)' Sensy_to:',
+SELF,IF=-Flo_resp.
       (EtaNam(np*(ko_p_tvar(jo)-1)+ko_p_kvar(jo))
                           ,(charit2(io+2+np*(jo-1)),io=1,np-1),jo=1,nxp)
+SELF,IF=Flo_resp.
       (ozvf(jo)         ,(charit2(io+2+np*(jo-1)),io=1,np-1),jo=1,nxp-nyp)
                            ,('$'//EtaNam(np*(ko_p_tvar(jo)-1)+ko_p_kvar(jo))
                 ,(charit2(io+2+np*(jo-1)),io=1,np-1),jo=nxp-nyp+1,nxp-nyp2)
                             ,(EtaNam(np*(ko_p_tvar(jo-nyp2)-1)+ko_p_kvar(jo-nyp2))
                       ,(charit2(io+2+np*(jo-1)),io=1,np-1),jo=nxp-nyp2+1,nxp)
+SELF.
      ;write(loPhit,1011)' time \\  of:',((EtaNam(io),io=1,np),jo=1,nxp)
     ;OPEN(loPsit,FILE='sigma'//osuffix,STATUS='UNKNOWN');
      write(loPsit,1011)' Sensy_to:',
+SELF,IF=-Flo_resp.
       (EtaNam(np*(ko_p_tvar(jo)-1)+ko_p_kvar(jo))
                             ,(charit2(io+2+mp*(jo-1)),io=1,mp-1),jo=1,nxp)
+SELF,IF=Flo_resp.
       (ozvf(jo)           ,(charit2(io+2+mp*(jo-1)),io=1,mp-1),jo=1,nxp-nyp)
                             ,('$'//EtaNam(np*(ko_p_tvar(jo)-1)+ko_p_kvar(jo))
                   ,(charit2(io+2+mp*(jo-1)),io=1,mp-1),jo=nxp-nyp+1,nxp-nyp2)
                        ,(EtaNam(np*(ko_p_tvar(jo-nyp2)-1)+ko_p_kvar(jo-nyp2))
                   ,(charit2(io+2+mp*(jo-nyp-1)),io=1,mp-1),jo=nxp-nyp2+1,nxp)
+SELF.
      ;write(loPsit,1011)' time \\  of:',((PhiNam(io),io=1,mp),jo=1,nxp);
    ;if mobs.gt.0
     < OPEN(loPsot,FILE='sigmo'//osuffix,STATUS='UNKNOWN');
      write(loPsot,1011)' Sensy_to:',
+SELF,IF=-Flo_resp.
     (EtaNam(np*(ko_p_tvar(jo)-1)+ko_p_kvar(jo))
                          ,(charit2(io+2+mobs*(jo-1)),io=1,mobs-1),jo=1,nxp)
+SELF,IF=Flo_resp.
     (ozvf(jo)           ,(charit2(io+2+mobs*(jo-1)),io=1,mobs-1),jo=1,nxp-nyp)
                              ,('$'//EtaNam(np*(ko_p_tvar(jo)-1)+ko_p_kvar(jo))
                ,(charit2(io+2+mobs*(jo-1)),io=1,mobs-1),jo=nxp-nyp+1,nxp-nyp2)
                         ,(EtaNam(np*(ko_p_tvar(jo-nyp2)-1)+ko_p_kvar(jo-nyp2))
                ,(charit2(io+2+mobs*(jo-nyp-1)),io=1,mobs-1),jo=nxp-nyp2+1,nxp)
+SELF.
      ;write(loPsot,1011)' time \\  of:',((PhiNam(mp+io),io=1,mobs),jo=1,nxp)
;
     >;
   >;
\end{verbatim}
%

On a en plus les tableaux de pointeurs et de calcul à la fois pour Floquet et 2.
Dans le cas où on ne s'intéresse qu'à une partie des Floquet, Floquet2 reste OK,
mais il y aura un travail spécifique dans un Floquet2 pour extraire l'information
tronquée des pointeurs --- si c'est les premiers, facile ...

{\bf question subs~:} Que se passe-t-il avec des variables constantes
de la liste des états ? Peut-on dans ce cas restreindre la dim pour le calcul complet ?


Pour les notations, on reprend les noms de variable avec un \$ en préfixe (\$T\_WV\_Eta),
et on crée un tableau character*2 ozvf(np) pour les Floquets, qui disent si r\# (reel)
ou R\# I\# (cmplx).

On aurait alors pour un calcul de gains :
%
\begin{verbatim}
 Sensy_to_Eta
 < 
           var: T_ocean_eta   , pert: init;
           var: T_tropo_eta   , pert: init;
           var: T_strato_eta  , pert: init;
           var: Q_tropo_eta   , pert: init;
           var: CO2_conc_eta  , pert: init; 
 >;
! Rep-eff rétroaction
 Sensy_to_Phi
 <         var: T_WV_phi      , pert: step_eff;
 >;
! moderateur sensib a xCO2 :
!***  Phi_t(5,5)=3.30d-04; "XXX"
       step_psi(6)=3.30d-04;
\end{verbatim}
%
en gardant sensy-to-var pour le standard (les macros tueront celles qu'il faut).
Rq~: faudra dans ce cas tuer les sources ?
\bigskip

En résumé\index{Floquet!effet-en-retour}., considérant l'une des formes dynamiques possibles de~(cf~\ref{sensg})~:
%
\begin{equation}
\left \lbrace
\begin{array}{lll}
s_g(t,\tau) &=& \int_{\tau}^{t} \Phi (t,s)\ket{b(s)} \, ds\nonumber\\
\sigma_g(t,\tau) &=& \bra{c^{\dagger}(t)}\int_{\tau}^{t} {\Phi (t,s)\ket{b(s)}}\, ds
\label{sensgencore}
\end{array}
\right.
\end{equation}
%
%
on espère qu'une décomposition du propagateur {\it à la Floquet} permettra d'affiner le
regard sur les dynamiqes plus fines sous-jacentes au modèle. On a deux décompositions utiles
en pratique --- une seule sera {\it in fine} retenue~:
%
\begin{eqnarray}
\Phi(t,\tau) &=& \sum_i e^{\lambda_i(t-\tau)}\ext{\psi^F_i(t)}{\psi^{-F}_i(\tau)}\nonumber\\
             &=& \sum_i \ext{y_i(t)}{y^{-1}_i(\tau)}\nonumber
\label{sensgencorep}
\end{eqnarray}
%
où nous avons ajouté celle en $Y^F$ liée aux Floquets par~$\ket{\psi_i ^F(t)}=\ket{y_i(t)}\exp(-\lambda_i t)$. 

On doit de plus remplacer les amortissements simples mailles de Floquet par une forme plus complexe, et
pour se faire, déterminer les passages à $t+T$, que l'on renouvelle ici en $Y^F$.

En revenant à la décomposition $Y(t) = \Phi(t,0)Z$ (de \ref{adjes}), avec les propriétés\\
$Y(t+T)=Y(t)\exp(\Lambda T),\quad Y(0)=Z\quad$ et $\Phi(t,\tau)=Y(t)Y^{-1}(\tau)$\footnote{puisque $Y$
est la matrice des solutions fondamentales ayant $Z$ comme C.I.}, on peut écrire~:
%
\begin{eqnarray}
s_i^g(t+T) &=& \int_0 ^T{\sum_i e^{\lambda_i T}\ket{y_i(t)}\scal{y^{-1}_i(\tau)}{b(\tau)}}\,d\tau\nonumber\\
       &&+ \int_T ^{T+t}{ Y(t+T)Y^{-1}(\tau)\ket{b(\tau)}}\,d\tau\nonumber
\end{eqnarray}
%
et après changement dans la deuxième intégrale et en résolvant les périodicités~:
%
\begin{eqnarray}
s_i^g(t+T) &=& s_i^g(t) + e^{\lambda_i T}\overline{\sigma_i} (t)\nonumber\\
\overline{\sigma_i}(t) &=& \ket{y_i(t)}\int_0 ^T{ \scal{y^{-1}_i(\tau)}{b(\tau)}}\,d\tau\nonumber
\end{eqnarray}
%
qui mène directement en amont de la formule~(\ref{fTF2}).

 La décomposition précédente est celle qui permet
de passer aux EDP sans dériver de matrice, et conduit de plus au système déjà TEFisé~($s,\sigma$).
On a ainsi deux déterminations à mener~: i) une réponse de base
(ou maille) de $t=0$ à la période; ii) une fonction d'amortisssement d'une période à
l'autre~:$\overline{\sigma_i}(t)$.

On utilise la forme intégrale du vecteur des sensibilités pou dériver les EDP,
comme ici en $Y$~:
%
\begin{equation}
\left \lbrace
\begin{array}{lll}
\partial_t s_i^g (t) &=& A(t)s_i^g(t) + \ket{y^F_i(t)}\scal{y^{-1}_i(t)}{b(t)}\nonumber\\
\sigma_g(t,\tau) &=& \bra{c^{\dagger}(t)}s_i^g
\label{miniF3}
\end{array}
\right.
\end{equation}
%

La forme intégrale de $\overline{\sigma_i}(t)$ est déterminée par
 $\partial_t\overline{\sigma_i}(t)=A\overline{\sigma_i}(t)$,
sans source donc; mais il faut encore déterminer $\overline{\sigma_i}(0)$, soit~:
%
\begin{eqnarray}
\overline{\sigma_i}(0) &=& \ket{y_i(0)}\int_0 ^T{ \scal{y^{-1}_i(\tau)}{b(\tau)}}\,d\tau\nonumber\\
              &=& e^{-\lambda_i T}\ket{y_i(T)}\int_0 ^T{ \scal{y^{-1}_i(\tau)}{b(\tau)}}\,d\tau\nonumber\\
              &=& e^{-\lambda_i T} s^g_i(T)
\end{eqnarray}
%
On commence donc par déterminer par intégration de~(\ref{miniF3}) les sensibilités-réponses mailles
de Floquet en avant de $0$ à $T$. On refait un run pour déterminer $\overline{\sigma_i}(t)$
en partant de $\overline{\sigma_i}(0)=e^{-\lambda_i T} s^g_i(T)$; on a alors tous les élements d'analyse
possible ...
\smallskip
Bref, il y a un peu d'essais concrets à faire avant de trouver ce que l'on cherche, et je compte sur la
nocturne du prochain Zircon pour y contribuer\footnote{des dangers du Floquage à Jussieu, on sera mieux
à Tarbes.}. Ceci dit, on devrait s'orienter vers un vrai calcul des vecteurs tefisés, dans la
mesure où les termes en $-\lambda\ket{\psi}$ régularisent par essence les divergences--contractions
du propagateurs {\it à la roduct}. Ce qui signifie pour \Minik un calcul de type {\bf oker} en complexes, qui serait
effectué après (ou en parallèle) à la résolution standard. L'avance deu $\Psi$ complet peut pauser
problème pour les systèmes très contractants --- pour lesquels on ne peut déterminer les
éléments de Floquet ayant déjà numériquement disparu à $T$ --- et il faut probablement se rabattre
sur une avance conjointe des lignes de $\Psi^{-F}$. Pour revenir aux systèmes rapidement contractants,
qui vont exhiber des multiplicateurs de Floquet très faibles menant à une impossibilité numérique
de calculer les V.P. correspondants, que faire ?

\bigskip
{\it quid du Wronskien ? fait le 23/10/2006, cf~(\ref{wronskF}).}


\medskip{\it il y a également l'examen dans ce contexte de la notion de ``cycle'' climatologique
(de l'eau, du carbone, de l'énergie mécanique atmosphérique ...}

En entendant par là qu'il est légitime de considérer un modèle qui, en suivant une molécule finalement recyclée
par exemple, sous entend de manière certes très implicite l'idée de récurrence, voire de périodicité.
La question est comment formuler\footnote{cf l'article P.Humbert et al, chap.~6~: On the relative Humidity
of the Atmosphere (note du 13 mai 2009). } un tel modèle ?

\printindex
%
\tableofcontents
\end{document}
\end

