diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index b343c827d835b9b7abc3200b13410a99f6bfc5c4..eaf07eaba0736a2deffe42b888daf727aa1c3958 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -50,88 +50,39 @@
 \maketitle
 \tableofcontents
 
-\section{Divers}
+\section{Extraction des tourbillons}
+
+\subsection{Divers}
 
 Chaque ligne de niveau de la hauteur correspond à une ligne de courant
 de la vitesse.
 
-La zone intérieure à la ligne de vitesse maximale est le c\oe{}ur du
-tourbillon. On cherche les recouvrements entre ces c\oe{}urs, l'idée
-étant que les extremums à deux dates correspondent au même tourbillon
-si les c\oe{}urs se recouvrent.
-
-But du programme : produire une base de données de tourbillons et un
-graphe. Les sommets du graphe sont les tourbillons. Les arcs relient
-les tourbillons dont les c\oe{}urs se recouvrent à deux dates
-successives.
-
-Le graphe est orienté acyclique. L'ordre topologique est tout trouvé :
-numérotation par dates successives. Le graphe non orienté associé est
-une forêt.
-
-Pour comparer les tourbillons à deux dates, il semble coûteux de
-n'utiliser que les listes de tourbillons aux deux dates. Cela implique
-une double boucle sur les listes de tourbillons, à l'intérieur d'une
-boucle sur les dates, soit potentiellement
-$8000 \times 8000 \times 5000$ comparaisons. Pour éviter cette
-comparaison aveugle, on peut conserver l'information sur la position
-des tourbillons dans des cartes aux différentes dates. La conservation
-en mémoire vive à 5000 dates semble aussi trop lourde. Il faut se
-restreindre aux cinq (max\_delta + 1) dates nécessaires à la
-profondeur de vue.
-
-Lorsqu'on compare des tourbillons à distance temporelle delta, on veut
-tester s'ils n'ont pas de prédécesseur ou successeur visible à une
-date intermédiaire. Il faut donc, si on crée un arc entre deux dates
-successives, ou delta arcs à distance delta (par interpolation de
-tourbillons invisibles), enregistrer ce fait en prévision du passage
-au delta supérieur. Mais d'un autre côté, on ne veut pas que le fait
-de mettre des arcs à l'étape delta ait une importance pour la suite de
-l'étape delta. (Et accessoirement, on ne veut pas que l'ordre dans
-lequel on passe les tourbillons à l'étape delta ait une importance.)
-On est donc obligé d'enregistrer, pour chaque tourbillon, non
-seulement s'il a des prédécesseurs visibles et s'il a des successeurs
-visibles, mais aussi à quelle distance temporelle.
-
-On ne cherche des superpositions qu'entres tourbillons visibles, en
-excluant les tourbillons interpolés.
-
-Possibilité de parallélisme dans l'extraction des tourbillons mais au
-prix d'une conservation de ssh, uv, v en mémoire vive à plusieurs
-dates.
-
-On pourrait écrire l'algorithme principal en explicitant les
-get\_snpashot, en explicitant le cas m < n\_proc et en supprimant la
-variable k\_end\_main\_loop. Ce serait mieux du point de vue de la
-performance, mais je pense que ce serait finalement moins clair.
-
-Cf. figure (\ref{fig:call_graph}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{call_graph}
-  \caption{Arbre des appels du programme de détection de tourbillons.}
-  \label{fig:call_graph}
-\end{figure}
-
-Cas pathologiques. out\_cont est vide si et seulement s'il n'y a pas
-de bon contour à innermost\_level. Le champ speed de extremum.dbf vaut
-missing\_speed si et seulement si radius4 == 0 ou speed\_outerm dans
-set\_max\_speed est NaN. Si radius4 $\ge$ 2 et speed dans extremum.dbf
-vaut missing\_speed alors speed\_outerm dans set\_max\_speed est NaN,
-e\%speed\_cont est non vide dans max\_speed\_contour.shp et
-mean\_speed sur e\%speed\_cont a donné un NaN. En d'autres termes, les
-deux appels à mean\_speed ont tous les deux produit un
-NaN. (e\%speed\_cont est non vide dans max\_speed\_contour.shp et
-speed dans extremum.dbf vaut missing\_speed) si et seulement si
-(radius4 $\ge$ 2 et speed dans extremum.dbf vaut missing\_speed).
-
-Cas où les composantes out\_cont ou speed\_cont d'un tourbillon sont
-vides, c'est-à-dire égales à null\_ssh\_contour(). On peut distinguer
-les cas du tableau (\ref{tab:null_ssh_contour}).
+But du programme : produire une base de données de tourbillons. Le
+programme ne traite qu'une date.
+
+Cas pathologiques. \verb+out_cont+ est vide si et seulement s'il n'y a
+pas de bon contour à \verb+innermost_level+. Le champ speed de
+extremum.dbf vaut \verb+missing_speed+ si et seulement si radius4 == 0
+ou \verb+speed_outerm+ dans \verb+set_max_speed+ est NaN. Si radius4
+$\ge$ 2 et speed dans extremum.dbf vaut \verb+missing_speed+ alors
+\verb+speed_outerm+ dans \verb+set_max_speed+ est NaN,
+\verb+e%speed_cont+ est
+non vide dans \verb+max_speed_contour.shp+ et \verb+mean_speed+ sur
+\verb+e%speed_cont+ a donné un NaN. En d'autres termes, les deux
+appels à \verb+mean_speed+ ont tous les deux produit un
+NaN. (\verb+e%speed_cont+ est
+non vide dans \verb+max_speed_contour.shp+ et speed dans extremum.dbf
+vaut \verb+missing_speed+) si et seulement si (radius4 $\ge$ 2 et
+speed dans extremum.dbf vaut \verb+missing_speed+).
+
+Cas où les composantes \verb+out_cont+ ou \verb+speed_cont+ d'un
+tourbillon sont vides, c'est-à-dire égales à
+\verb+null_ssh_contour()+. On peut distinguer les cas du tableau
+(\ref{tab:null_ssh_contour}).
 \begin{table}[htbp]
   \centering
   \begin{tabular}{lllll|l}
-    out\_cont & valid & radius4 & max\_speed
+    \verb+out_cont+ & valid & radius4 & \verb+max_speed+
     & fraction & note \\
     \hline
     null & F & 0 & \np{e4} & 5 \% & 1 \\
@@ -152,11 +103,11 @@ les cas du tableau (\ref{tab:null_ssh_contour}).
 \end{table}
 Sur la région 5, sur 260 extremums, je trouve les nombres de cas
 suivants pour les 5 lignes du tableau (\ref{tab:null_ssh_contour}) :
-12, 15, 7, 122, 54. Au total : 210 cas où les composantes out\_cont ou
-speed\_cont d'un tourbillon sont vides, sur 260 tourbillons, soit 81
+12, 15, 7, 122, 54. Au total : 210 cas où les composantes \verb+out_cont+ ou
+\verb+speed_cont+ d'un tourbillon sont vides, sur 260 tourbillons, soit 81
 \% des tourbillons.
 
-\section{Calcul de la vitesse azimutale par rapport à l'extremum}
+\subsection{Calcul de la vitesse azimutale par rapport à l'extremum}
 
 Notons $C$ la position de l'extremum et $M$ un point où l'on veut
 calculer la vitesse azimutale, $\vec k_C$ et $\vec k$ les vecteurs
@@ -225,8 +176,7 @@ distance entre l'extremum et $M$, la vitesse azimutale est :
   V_\theta = \frac{x v - y u}{\rho}
 \end{equation*}
 
-
-\section{Bons contours}
+\subsection{Bons contours}
 
 Soit un extremum de SSH. Un \og bon contour\fg{} de cet extremum, est
 une ligne de niveau de SSH fermée contenant cet extremum et ne
@@ -247,7 +197,7 @@ strictement extérieur au contour.  Je crois, sans l'avoir démontré,
 que, pour chacun des quatre chemins de sortie, la suite des valeurs de
 SSH aux points rencontrés à l'intérieur du contour (tous les points
 sauf le dernier) varie de façon monotone. J'utilise cette hypothèse
-dans outermost\_possible\_level.
+dans \verb+outermost_possible_level+.
 
 Admettons les propriétés ci-dessus. Considérons l'un des quatre
 chemins de sortie. Indiçons, mettons avec $l$ variant de 1 (pour
@@ -304,7 +254,7 @@ contour) forcément les points voisins de l'extremum. Donc le chemin de
 sortie de ce bon contour le plus éloigné serait de longueur supérieure
 à 3.
 
-\section{Problèmes de dégénérescence}
+\subsection{Problèmes de dégénérescence}
 
 Possibilité que deux points adajacents aient exactement la même valeur
 de SSH. J. Tierny mentionne ce type de problème dans son
@@ -333,8 +283,8 @@ faudrait faire ces changements jusque dans les procédures de recherche
 de contour.
 
 Si j'admets la possibilité de dégénérescence alors les procédures
-impactées sont get\_1\_outerm, local\_extrema, inside et
-outermost\_possible\_level.
+impactées sont \verb+get_1_outerm+, \verb+local_extrema+, inside et
+\verb+outermost_possible_level+.
 
 Je cherche un bon contour en ne spécifiant qu'un seul point, qui doit
 être intérieur au bon contour, et un ensemble de points qui doivent
@@ -348,28 +298,25 @@ contour autour des deux. Il faut donc n'en retenir qu'un des deux
 comme valide. Je peux prendre arbitrairement le premier rencontré dans
 l'ordre de parcours du tableau.
 
-Il faut prévoir dans get\_1\_outerm la possibilité de ne pas
+Il faut prévoir dans \verb+get_1_outerm+ la possibilité de ne pas
 trouver de bon contour. La logique est impactée. \`A strictement
-parler, la valeur de innermost\_level ne permet plus de prévoir si
+parler, la valeur de \verb+innermost_level+ ne permet plus de prévoir si
 l'amplitude va être suffisante, et on ne sait plus quelle valeur
-intiale de level\_good prendre dans get\_1\_outerm. Je garde
-néanmoins innermost\_level comme valeur intiale de level\_good. La
+intiale de \verb+level_good+ prendre dans \verb+get_1_outerm+. Je garde
+néanmoins \verb+innermost_level+ comme valeur intiale de \verb+level_good+. La
 conséquence est que je risque de ne pas détecter de contour autour
 d'un extremum dégénéré alors qu'on pourrait peut-être en trouver un en
-prenant une valeur initiale de level\_good plus proche de
+prenant une valeur initiale de \verb+level_good+ plus proche de
 l'extremum. On peut aussi imaginer le cas bizarre de trois maximums
 (par exemple) égaux adjacents. Dans ce cas, l'algorithme choisi
 conserve deux de ces maximaux et il n'y a pas de bon contour possible,
 alors qu'on pourrait en principe en trouver. On se désintéresse de ces
 cas pathologiques.
 
-\section{Entrées et sorties}
+\subsection{Sorties}
 
 Format shapefile. Un type de géométrie par fichier donc il faut
-séparer les contours et les positions des extremums. Par ailleurs, on
-n'enregistre pas de contours pour les tourbillons interpolés. Les
-tourbillons interpolés sont écrits avec les autres, en utilisant une
-forme NULL pour le polygone correspondant à un contour interpolé.
+séparer les contours et les positions des extremums.
 \begin{itemize}
 \item \verb+extremum_$m.shp+ : points
 \item \verb+extremum_$m.dbf+
@@ -378,7 +325,7 @@ forme NULL pour le polygone correspondant à un contour interpolé.
   de vitesse sur le contour de vitesse maximale
 \item \verb+outermost_contour_$m.shp+ : polygones
 \item \verb+outermost_contour_$m.dbf+ : aire, valeur de SSH, indice de
-  date, indice de tourbillon à cette date, twice, radius4
+  date, indice de tourbillon à cette date, radius4
 \item \verb+max_speed_contour_$m.shp+ : polygones
 \item \verb+max_speed_contour_$m.dbf+
   : aire, valeur de SSH, indice de date, indice de tourbillon à cette
@@ -400,10 +347,7 @@ L'indice de date et l'indice de tourbillon dans
 sont en principe inutiles puisqu'ils sont aussi dans
 \verb+extremum_$m.dbf+ et que les tourbillons sont stockés dans le
 même ordre dans les trois shapefiles. La redondance me semble
-préférable pour la solidité et le confort d'utilisation. La
-méta-donnée logique \og interpolé \fg{} dans \verb+extremum_$m.dbf+
-est nécessaire parce qu'il est possible de ne pas trouver de contour
-extérieur autour d'un extremum détecté.
+préférable pour la solidité et le confort d'utilisation.
 
 Les réels dans un fichier shp sont stockés en double précision donc
 les polygones occupent chacun : 16 B $\times$ nombre moyen de côtés
@@ -417,56 +361,9 @@ coup sûr les polygones pour 16 ans dans un seul fichier, même en
 séparant cyclones et anticyclones. Faire un fichier par an ?  Séparer
 en outre cyclones et anticyclones ?
 
-Nous stockons le graphe au format liste d'arcs (en texte),
-c'est-à-dire, pour chaque sommet du graphe, une ligne :
-\begin{verbatim}
-indice date, indice tourbillon, indice date, indice tourbillon, poids
-\end{verbatim}
-Fichier
-\verb+edgelist_$m.txt+. L'estimation de la taille du fichier, en
-comptant 10 caractères par entier, 13 caractères par réel, \np{5e7}
-sommets et deux arcs par sommet en moyenne, est de 6 GiB.
-
-Problème des sommets isolés. L'information est déjà présente, de façon
-implicite, dans les fichiers évoqués ci-dessus. En effet, on a le
-nombre de tourbillons à chaque date d'après la liste des identifiants
-de tourbillons dans les shapefiles, et on a la liste des tourbillons
-non isolés d'après la liste des arcs. Mais la liste des sommets isolés
-n'est donc pas très facilement trouvée. On ajoute par conséquent les
-fichiers de sortie suivants :
-\begin{description}
-\item[number\_eddies\_\$m.csv] Une colonne indice de date, une colonne
-  nombre de tourbillons visibles, une colonne nombre de tourbillons
-  interpolés.
-\item[isolated\_nodes\_\$m.txt] Sur chaque ligne : indice de date, indice
-  de tourbillon. Peut être lu comme une simple liste
-  d'adjacence. Remarque : par construction, les tourbillons interpolés
-  ne sont jamais isolés.
-\end{description}
-
-Les entrées-sorties sont dans : l'algorithme principal directement
-(lecture de longitude et latitude), get\_snapshot (lecture de SSH, u,
-v), successive\_overlap (écriture dans
-\verb+edgelist_$m.txt+, arcs entre k\_begin(m) et k\_end\_main\_loop(m)),
-non\_successive\_overlap directement (écriture dans
-\verb+edgelist_$m.txt+),
-write\_eddy (écriture dans extremum\_\$m, outermost\_contour\_\$m,
-max\_speed\_contour\_\$m), dispatch\_snapshot directement (écriture
-dans isolated\_nodes\_\$m, number\_eddies\_\$m). Les processus
-écrivent dans des fichiers indicés par le numéro de processus donc
-deux processus ne peuvent écrire dans le même fichier. Cf. figure
-(\ref{fig:input_output}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{input_output}
-  \caption{Entrées et sorties du programme en Fortran.}
-  \label{fig:input_output}
-\end{figure}
-
-\section{Déclarations pour l'algorithme principal}
+\subsection{Structures de données}
 
-On suppose que le maillage est le même à toutes les dates et qu'il est
-régulier.
+On suppose que le maillage est régulier.
 
 Type dérivé polyline : défini dans \verb+Contour_531+. Mémoire
 utilisée par un scalaire de ce type :
@@ -493,19 +390,16 @@ calculer une vitesse.
 
 Nous avons besoin que la composante extr\_map de snapshot soit étendue
 en longitude dans le cas d'un domaine périodique pour son utilisation
-par set\_max\_speed et successive\_overlap.
+par \verb+set_max_speed+.
 
-Les tailles ci-dessous sont calculées pour max\_delta = 4, 16 ans,
-soit environ 5000 dates, 8000 tourbillons par date, 30 côtés par
+Les tailles ci-dessous sont calculées pour 8000 tourbillons, 30 côtés par
 polygone (sans différence entre contour extérieur et contour de
 vitesse maximale), 720 latitudes par 1440 longitudes.
 \begin{description}
-\item[max\_delta] Scalaire entier. Intervalle maximal d'indices de
-  dates auxquelles on cherche des recouvrements. Normalement 4.
 \item[longitude, latitude] Vecteurs de réels
 \item[ssh, u, v] Tableaux de réels, deux dimensions (nlon,
   nlat). Environ 4 MiB chacun.
-\item[eddy] La comparaison entre out\_cont\%ssh et ssh\_extr donne la
+\item[eddy] La comparaison entre \verb+out_cont+\%ssh et ssh\_extr donne la
   distinction anticyclone -- cyclone mais il est quand même plus
   confortable d'avoir cette information directement dans les
   sorties. Pour cela, j'ajoute un champ cyclone. Mémoire utilisée par
@@ -539,10 +433,6 @@ soit environ \np{0.5} KiB.
   soit environ 4 MiB.
 \item[snapshot\%ind\_extr] Taille : $2 n_{ve}$ mots, soit environ
   \np{0.06} MiB.
-\item[flow] : Vecteur de type snapshot, de taille max\_delta + 1,
-  chaque élément correspondant à une date. Environ 40 MiB.
-\item[dist\_lim] Scalaire entier constant, = 3° / résolution.
-\item[weight\_delta] Scalaire réel.
 \end{description}
 On a :
 \begin{align*}
@@ -563,147 +453,724 @@ global. J'ai donc choisi donc de remplacer ces champs par un champ
 ind\_extr dans snapshot. Et j'ajoute à la place des champs longitude,
 latitude dans le type eddy.
 
-i désigne un indice de tourbillon, j un indice de position dans la
-fenêtre temporelle (entre 1 et max\_delta + 1) et k un indice de
-date. Cf. figure (\ref{fig:window}) et algorithme
-(\ref{alg:main_sequential}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics{window}
-  \caption{Indice de position dans la fenêtre glissante et indice de
-    date.}
-  \label{fig:window}
-\end{figure}
-\begin{algorithm}[htbp]
-  \begin{algorithmic}
-    \STATE \COMMENT{\{max\_delta $\ge 1$; n\_dates $\ge$ max\_delta +
-      1\}} 
+\subsection{Résumé du problème de géométrie algorithmique}
 
-    \STATE \COMMENT{prologue}
+\begin{algorithmic}
+  \STATE recherche de tous les extremums locaux
 
-    \STATE entrer(longitude, latitude)
-    \STATE \COMMENT{longitude et latitude dans l'ordre croissant} 
-    \FOR{k = 1 \TO max\_delta + 1}
+  \FOR{chaque extremum local}
 
-    \STATE appel de get\_snapshot(flow(k), k)
+  \STATE recherche du \og bon contour le plus extérieur\fg{}
+  \ENDFOR
+\end{algorithmic}
 
-    \ENDFOR
-    \FOR{k = 2 \TO max\_delta + 1}
+\og Bon contour\fg{} : contient l'extremum considéré et aucun
+autre. \og Le plus extérieur\fg{} : plus grande différence de niveau
+avec l'extremum considéré, ou encore amplitude maximale.
 
-    \STATE appel de successive\_overlap(k, k, flow)
-    \ENDFOR
-    \FOR{delta = 2 \TO max\_delta}
+Recherche du bon contour le plus extérieur : par dichotomie.
+\begin{algorithmic}
+  \STATE $a_1 = 0$
 
-    \FOR{k = delta + 1 \TO max\_delta + 1}
+  \STATE $a_2$ = max ou min de SSH sur le domaine - SSH de l'extremum cible
 
-    \STATE appel de non\_successive\_overlap(k, k, delta, flow)
+  \WHILE{$a_2 - a_1 >$ précision}
 
-    \ENDFOR
-    \ENDFOR
+  \STATE recherche du bon contour à l'amplitude
+  $\frac{a_1 + a_2}{2}$
 
-    \STATE \COMMENT{boucle principale}
+  \IF{bon contour trouvé}
 
-    \FOR{k = max\_delta + 2 \TO n\_dates}
+  \STATE $a_1 = \frac{a_1 + a_2}{2}$
 
-    \STATE appel de dispatch\_snapshot(flow(1), k - max\_delta - 1)
-    \STATE flow(:max\_delta) = flow(2:)
+  \ELSE
 
-    \STATE appel de get\_snapshot(flow(max\_delta + 1), k)
+  \STATE $a_2 = \frac{a_1 + a_2}{2}$
 
-    \STATE appel de successive\_overlap(max\_delta + 1, k, flow)
+  \ENDIF
 
-    \FOR{delta = 2 \TO max\_delta}
+  \ENDWHILE
+\end{algorithmic}
 
-    \STATE appel de non\_successive\_overlap(max\_delta + 1, k, delta,
-    flow)
+Recherche d'un bon contour à une amplitude donnée. On cherche tous les
+contours à cette amplitude. Pour chaque contour, pour chaque extremum,
+on regarde si l'extremum est à l'intérieur du contour (test de point
+intérieur à polygone).
 
-    \ENDFOR
-    \ENDFOR
-    \STATE \COMMENT{épilogue}
+Les procédures géométriques utilisées :
+\begin{itemize}
+\item polygon\_area\_2d (bibliothèque Geometry)
+\item polygon\_contains\_point (bibliothèque Geometry)
+\item find\_contours\_reg\_grid (bibliothèque Contour\_531)
+\end{itemize}
+travaillent dans le plan
+et non sur la sphère. Ces procédures sont appelées par :
+\begin{itemize}
+\item spherical\_polyline\_area
+\item inside\_4
+\item good\_contour
+\end{itemize}
+Dans le plan, par deux points passe un unique segment, et la donnée
+d'un triangle est équivalente à la donnée de l'ensemble de ses trois
+sommets. Sur la sphère, par deux points passent deux arcs de grand
+cercle, et plusieurs triangles peuvent être associés à un ensemble de
+trois sommets. Par ailleurs, les procédures géométriques dans le plan
+ne peuvent pas admettre en arguments d'entrée des longitudes définies
+modulo $2 \pi$. Pour que les valeurs de longitudes en entrée des
+procédures géométriques soient contraintes sans ambiguïté, nous
+imposons que tout domaine de longitudes en entrée des procédures
+géométriques soit de longueur strictement inférieure à $\pi$.
+\begin{equation}
+  \label{eq:length_pi}
+  \forall (\lambda_1, \lambda_2) \in \mathbb{R}^2, \forall k \in \mathbb{Z},
+  |\lambda_2 - \lambda_1| < \pi
+  \Rightarrow |\lambda_2 + 2 k \pi - \lambda_1| \ge \pi
+\end{equation}
+(Toute modulation de $2 \pi$ d'une longitude viole la contrainte.) Les
+procédures qui appellent les procédures géométriques doivent se
+charger d'ajouter le multiple de $2 \pi$ pour ramener les longitudes
+dans un intervalle de longueur strictement inférieure à $\pi$. La
+contrainte de domaine de longitudes s'applique aux procédures
+\verb+get_1_outerm+ et \verb+set_max_speed+.
 
-    \FOR{j = 1 \TO max\_delta + 1}
+\subsection{Sous-algorithmes}
 
-    \STATE appel de dispatch\_snapshot(flow(j), n\_dates - max\_delta - 1
-    + j)
-    \ENDFOR
-  \end{algorithmic}
-  \caption{Algorithme principal, séquentiel}
-  \label{alg:main_sequential}
-\end{algorithm}
+\subsubsection{get\_1\_outerm}
 
-\section{Décomposition de domaine}
+Au lieu d'initialiser l'amplitude maximale à 4 m plus du bruit, comme
+dans le programme pour Matlab :
+\begin{algorithmic}
+   \STATE appel de random\_number(level\_try)
 
-Cf. figure (\ref{fig:processes}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{processes}
-  \caption{Raccordement des processus. Les comparaisons représentées
-    ici en couleurs sont faites par le processus m. En rouge dans la
-    boucle principale, en bleu, vert et marron dans l'épilogue.}
-  \label{fig:processes}
-\end{figure}
-Il faut, pour m < n\_proc :
-\begin{equation}
-  \label{eq:stitching}
-  \mathtt{k\_end}(m)
-  = \mathtt{k\_begin}(m + 1) + \mathtt{max\_delta} - 1
-\end{equation}
+   \STATE level\_try = e\%ssh\_extr $\pm$ [4 + (2 * level\_try -
+   1) * precision] \COMMENT{un peu de bruit sur l'amplitude, autour de
+     4}
+\end{algorithmic}
+et de supposer qu'il n'y a pas de bon contour avec cette amplitude, on
+cherche une amplitude maximale possible. Dans le cas le plus simple où
+noise\_around est faux, on purrait parcourir les points dans une
+direction jusqu'à ce que le sens de variation de SSH change. Dans mes
+tests, le nombre d'itérations dans la recherche par dichotomie est
+typiquement une dizaine, pour une précision de \np{0.5} mm. En
+considérant simplement l'extremum de ssh sur la fenêtre, la taille de
+l'intervalle initial [\verb+innermost_level+, level\_try] est plus
+grande. Sur mes tests, elle est plus grande typiquement d'un facteur 2
+environ seulement, donc le nombre d'itérations change peu.
 
-k\_end - k\_begin + 1 appels à get\_snapshot, pour les dates k\_begin
-à k\_end. k\_end - k\_begin + 1 appels à dispatch\_snapshot, pour les
-dates k\_begin à k\_end. Si m > 1 alors max\_delta envois. Si m <
-n\_proc alors max\_delta réceptions.
+Pour l'amplitude minimale, si je prends simplement l'amplitude
+correspondant à \verb+innermost_level+ alors l'algorithme échoue lorsque
+l'extremum est à un point du bord et que le point à \verb+innermost_level+
+est au bord. Dans ce cas, good\_contour échoue parce que Contour\_531
+ne ferme pas les contours passant par le bord. J'avais essayé, pour
+éviter le problème, de prendre une amplitude minimale légèrement
+inférieure mais cela pose des problèmes de cohérence dans
+l'algorithme.
 
-Nombre d'appels à get\_snapshot avec lecture. Si m < n\_proc : k\_end -
-k\_begin - max\_delta + 1. Si m = n\_proc : k\_end - k\_begin + 1.
+\subsubsection{get\_snapshot}
 
-successive\_overlap est appelé pour tous les indices k compris entre
-k\_begin + 1 et k\_end\_main\_loop, seulement. Nombre d'appels à
-successive\_overlap : k\_end\_main\_loop - k\_begin. Donc si m <
-n\_proc : k\_end - k\_begin - max\_delta + 1 appels. Si m = n\_proc : k\_end
-- k\_begin appels.
+Pour chaque tourbillon du résultat, le contour le plus extérieur doit
+être fermé, doit contenir l'extremum de SSH du tourbillon considéré et
+ne doit contenir aucun autre extremum de SSH avec une amplitude
+significative. Amplitude significative : supérieure à un certain
+minimum, entre l'extremum de SSH et le contour le plus
+lointain.
 
-Nombre d'appels à non\_successive\_overlap dans le prologue :
-\begin{equation*}
-  \sum_{\delta = 2} ^{\mathtt{max\_delta}}
-  (\mathtt{max\_delta} - \delta +1)
-  = \frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}
-\end{equation*}
-Dans la boucle principale : (max\_delta - 1)(k\_end\_main\_loop -
-k\_begin - max\_delta). Dans l'épilogue, si m < n\_proc :
-$\frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}$. Total si m
-< n\_proc : (max\_delta - 1)(k\_end - k\_begin - max\_delta +
-1). Total si m = n\_proc : (max\_delta - 1)(k\_end - k\_begin -
-$\frac{\mathtt{max\_delta}}{2}$).
+Dans une première version qui n'envisageait pas d'amplitude minimum,
+j'avais une procédure get\_eddy qui définissait un tourbillon, et qui
+trouvait donc le contour le plus extérieur et le contour de maximum de
+vitesse. Mais avec la possibilité d'une amplitude minimale, il faut
+d'abord se préoccuper de tous les contours extérieurs, qui ne sont
+plus indépendants les uns des autres, pour ne pas chercher inutilement
+des contours de maximum de vitesse. Dans le nouvel algorithme, il n'y
+a donc plus de procédure get\_eddy.
 
-Si k\_end - k\_begin était le même pour m = n\_proc et m <
-n\_proc, le processus m=n\_proc aurait des plus grands nombres
-d'appels aux trois procédures. Notons $\Delta k$ = k\_end - k\_begin,
-considéré comme une fonction de m. Si, pour m < n\_proc :
-\begin{equation}
-  \label{eq:charge}
-  \Delta k(m) \approx \Delta k(\mathtt{n\_proc}) + \mathtt{max\_delta}
-\end{equation}
-alors la charge devrait être à peu près équilibrée.
+Si le domaine considéré n'est pas global, comment traiter un
+tourbillon au bord du domaine ? La taille du contour le plus extérieur
+peut être artificiellement diminuée. (Ce n'est pas le même cas qu'un
+tourbillon près d'une côte.) C'est-à-dire qu'on pourrait trouver un
+contour plus extérieur avec un domaine élargi. Il pourrait être
+intéressant d'éliminer un tel tourbillon. Comment le détecter ? Le
+fait d'aller jusqu'au bord du domaine dans \verb+outermost_possible_level+
+n'est pas un critère suffisant, parce que \verb+outermost_possible_level+
+ne teste que la direction est. Il faudrait tester après l'appel de
+\verb+get_1_outerm+ si le cadre le plus extérieur touche le bord du
+domaine.
 
-On peut prendre, pour tout m :
-\begin{equation*}
-  \mathtt{k\_begin}(m)
-  = 1 + E[(m - 1) \mathtt{n\_dates} / \mathtt{n\_proc}]
-\end{equation*}
-(On prend la partie entière au lieu de nint pour pouvoir démontrer des
+\subsubsection{good\_contour}
+
+Trouve, si elle existe, une ligne de niveau à un niveau donné,
+contenant un point donné et ne contenant pas un ensemble de points
+donnés.  Le résultat est un polygone à zéro point s'il n'existe pas de
+ligne de niveau satisfaisante. Correspond à une partie de la fonction
+extraction\_eddies dans le programme pour Matlab.
+
+\subsubsection{local\_extrema}
+
+a(ind\_extr(1, k), ind\_extr(2, k)) est un extremum de a par rapport
+aux huit points voisins. extr\_map vaut 0 en un point où il n'y a pas
+d'extremum, le numéro de l'extremum dans ind\_extr en un point
+extremum. Accélération possible, au prix de tableaux supplémentaires,
+en utilisant le fait que deux maximums ne peuvent être côte à côte et
+deux minimums ne peuvent être côte à côte.
+
+J'avais d'abord codé le sens de l'extremum (maximum ou minimum) par le
+signe dans extr\_map. Mais c'est plus clair et plus simple sans ce
+codage. J'ajoute donc un argument vecteur cyclone à \verb+local_extrema+,
+qui prend moins de place qu'une carte supplémentaire. Pour le sens de
+l'extremum, je n'ai pas besoin que l'information soit sur une carte.
+
+Cf. aussi fontion \verb+skimage.feature.peak_local_max+ dans le paquet
+Python scikit-image. La fonction applique
+\verb+scipy.ndimage.maximum_filter+ et compare le résultat au tableau
+de départ. \verb+scipy.ndimage.maximum_filter+ applique
+\verb+scipy.ndimage.maximum_filter1d+ dans chaque dimension. Cf. aussi
+Douglas (1996 k0968).
+
+Cf. figure (\ref{fig:periodicity}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics{periodicity}
+  \caption{Cas de périodicité dans local\_extrema. Le tableau
+    représenté est extr\_map. Pour définir la case (2, j), il faut
+    avoir copié la case (m - 1, j - 1) (qui peut être non nulle). Pour
+    définir la case (m - 1, j), il faut avoir copié la case (2, j)
+    (qui peut être non nulle). (La case (2, j - 1) doit aussi avoir
+    été copiée mais cela a été fait au passage sur la colonne
+    précédente.)}
+  \label{fig:periodicity}
+\end{figure}
+Cas de périodicité : field(i + m - 2, :) = field(i, :).
+
+\subsubsection{max\_speed\_contour\_ssh}
+
+Il manque des tests pour lesquels direction vaut 1 et 4.
+
+\subsubsection{set\_all\_outerm}
+
+Notons ici :
+\begin{align*}
+  & u = \mathtt{urc(1)} \\
+  & l = \mathtt{llc(1)} \\
+  & s = \mathtt{step(1)} \\
+  & e = \mathtt{s\%list\_vis(i)\%coord\_extr(1)} \\
+  & w = \mathtt{corner\_window(1)} \\
+  & m = \mathtt{max\_radius(1)}
+\end{align*}
+Pour que la contrainte de domaine de longitudes dans \verb+get_1_outerm+
+soit satisfaite, il faut que, dans set\_all\_outerm, au moment de
+l'appel à \verb+get_1_outerm+ :
+\begin{equation}
+  \label{eq:max_radius}
+  (u - l) s < \pi  
+\end{equation}
+et que $e$ et outside\_points(1, :) soient dans l'intervalle
+$[w, w + (u - l) s]$. Pour que l'inégalité (\ref{eq:max_radius}) soit
+satisfaite, il suffit que :
+\begin{equation*}
+  2 m s < \pi
+\end{equation*}
+ce que nous imposons en amont de l'appel de set\_all\_outerm. Si
+periodic est vrai alors :
+\begin{equation*}
+  [w, w + (u - l) s] = [e - m s, e + m s]
+\end{equation*}
+donc $e$ appartient bien à cet intervalle. Si periodic est faux,
+$[w, w + (u - l) s]$ peut être plus réduit que $[e - m s, e + m s]$
+mais nous savons aussi que $e$ est dans l'intervalle [corner(1),
+corner(1) + (nlon - 1) s], donc $e$ est bien toujours dans
+$[w, w + (u - l) s]$.
+
+Il reste le problème de outside\_points(1, :). outside\_points(1, :)
+contient des coordonnées prises dans s\%list\_vis\%coord\_extr(1). Si
+le domaine est périodique et si l'extremum cible est proche en
+longitude de corner(1) par exemple alors outside\_points(1, :) peut
+contenir des longitudes proches de corner(1) + (nlon - 1) s. Il faut
+leur retrancher $2 \pi$. Nous posons le problème plus généralement de
+la façon suivante. Soit $\lambda \in \mathbb{R}$. On cherche $k \in
+\mathbb{Z}$ tel que :
+\begin{equation*}
+  \lambda + 2 k \pi \in [w, w + (u - l) s]
+\end{equation*}
+Supposons que $k$ existe (c'est le cas si $\lambda$ est une longitude
+retournée par nearby\_extr). Alors $k$ est unique (cf. propriété
+(\ref{eq:length_pi})). On a :
+\begin{equation*}
+  \lambda + 2 k \pi \in [w, w + 2 \pi[
+\end{equation*}
+C'est-à-dire :
+\begin{equation*}
+  k - 1 < \frac{w - \lambda}{2 \pi} \le k
+\end{equation*}
+Ou encore :
+\begin{equation*}
+  k = \left\lceil \frac{w - \lambda}{2 \pi} \right\rceil
+\end{equation*}
+
+
+\subsubsection{set\_max\_speed}
+
+Notons $d$ la distance dans l'espace des indices entre l'extremum
+cible et le countour extérieur. Par définition de la partie entière,
+les points à distance $E(d)$ de l'extremum cible sont à l'intérieur du
+contour. Certains points à distance $E(d) + 1$ peuvent être à
+l'extérieur. D'où : radius = $E(d) + 1$.
+
+Dans le cas où radius vaut 1, je calcule \verb+max_speed+. Ce calcul ne sera
+utile que dans weight si le tourbillon a une intersection avec avec un
+tourbillon à une autre date. Donc, pour certains tourbillons, ce
+calcul est inutile. Difficile d'éviter ces calculs inutiles. Il
+faudrait les reporter à weight mais alors il faudrait avoir accès dans
+weight à u et v, c'est-à-dire qu'il faudrait conserver u, v à
+plusieurs dates en dehors de \verb+get_snpashot+. Cela remettrait en
+question tout l'algorithme.
+
+\subsubsection{spherical\_polyline\_area}
+
+La procédure de calcul des contours donne une suite de sommets. Les
+courbes reliant ces sommets ne sont en fait pas imposés
+intrinsèquement par cette procédure. Par exemple, il n'est ni plus ni
+moins juste de supposer que ces courbes sont des portions de grands
+cercles ou de supposer qu'elles donnent des segments en projection
+plate carrée. Néanmoins, dans le reste du programme, pour décider si
+un contour entoure un extremum, ou pour calculer l'intersection de
+deux tourbillons, nous supposons que les contours sont des polygones
+en $(\lambda, \phi)$ (c'est-à-dire que les courbes reliant les sommets
+donnent des segments en projection plate carrée).
+
+Un domaine $S$ quelconque à la surface de la Terre a pour aire :
+\begin{equation*}
+  a^2  \iint_S \cos \phi\ \ud \phi\ \ud \lambda
+  = a^2  \iint_S \ud \lambda\ \ud(\sin \phi)
+\end{equation*}
+Nous avons une procédure qui calcule la surface $\iint \ud x\ \ud y$
+d'un polygone en $(x, y)$. Nous calculons la surface du polygone en
+$(\lambda, \sin \phi)$ passant par les sommets du
+contour. C'est-à-dire que nous supposons que les courbes reliant les
+somments donnent des segments en projection conforme cylindrique.
+
+Rémi fait un autre calcul : la surface du polygone en
+$(\lambda \cos \phi, \phi)$ passant par les sommets du contour. Le
+changement de variable :
+\begin{align*}
+  & x = \lambda \cos \phi \\
+  & y = \phi
+\end{align*}
+a pour jacobien $\cos \phi$ donc on a bien :
+\begin{equation*}
+  \iint \cos \phi\ \ud \phi\ \ud \lambda = \iint \ud x\ \ud y
+\end{equation*}
+Les deux surfaces obtenues, celle du polygone en
+$(\lambda, \sin \phi)$ et celle du polygone en
+$(\lambda \cos \phi, \phi)$ sont très proches si les sommets sont très
+proches. C'est le cas puisque les sommets sont séparés au maximum de
+\np{0.25}° en longitude et latitude. On peut aussi le voir en
+considérant la formule de calcul de surface du polygone :
+\begin{align*}
+  S & = \frac{1}{2} \sum_i x_i (y_{i + 1} - y_{i - 1}) \\
+    & =
+      \frac{1}{2} \sum_i \lambda_i \cos \phi_i (\phi_{i + 1} - \phi_{i - 1}) \\
+    & \approx \frac{1}{2} \sum_i \lambda_i (\sin \phi_{i + 1} - \sin \phi_{i - 1})
+\end{align*}
+si les valeurs successives des $\phi_i$ sont proches.
+
+\subsubsection{write\_eddy}
+
+On pourrait passer en arguments d'entrée les numéros des
+champs dans les fichiers DBF mais cela ferait 13
+arguments.
+
+\subsection{Tests}
+
+Pour faire des tests : données au 1\ier{} janvier 2006. Tests sur
+différents domaines. Cf. figure (\ref{fig:regions}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{regions}
+  \caption{Régions pour les tests.}
+  \label{fig:regions}
+\end{figure}
+Régions 1, 2 et 3 : 29 $\times$ 17, 57 $\times$ 33; 120 $\times$ 120,
+contenant respectivement 8, 27 et 216 extremums.
+
+Pour la région 1, coin en :
+\begin{equation*}
+  \lambda = \np{5.125}^\circ, \phi = -\np{36.375}^\circ
+\end{equation*}
+soit $(\np[rad]{0.0894}, - \np[rad]{0.6349})$. Indices base 1 du coin
+: 21, 215. Pas de grille : \np{0.25}°, soit \np[rad]{0.0044}. 5
+minimums et 3 maximums dans ce cadre. Cf. figure (\ref{fig:extrema}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{extrema}
+  \caption{Région 1, lignes de niveau de SSH (sans signification
+    particulière). Maximums : points rouges, minimums : points bleus.}
+  \label{fig:extrema}
+\end{figure}
+Minimum cible pour le test \verb+Get_1_outerm+ en :
+\begin{equation*}
+  \lambda = \np{9.625}^\circ, \phi = - \np{33.875}^\circ
+\end{equation*}
+Indices base 1 de ce minimum : 39, 225. SSH en ce mimimum :
+\np{0.2759}. Extremum numéro 6 dans ce cadre. Contour le plus
+intérieur : SSH = \np{0.293300003} m. Contour le plus extérieur : SSH
+= \np{0.460744} m. Contour de vitesse maximale : SSH = \np{0.3815}
+m. Cf. figures (\ref{fig:test_get_eddy}) et (\ref{fig:test_get_snapshot}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{test_get_eddy}
+  \caption{Test Set\_max\_speed. Projection stéréographique centrée
+    sur un minimum de SSH. Champ de vitesse, contour de SSH le plus
+    extérieur en bleu, contour de SSH de vitesse maximale en vert. Le
+    minimum de SSH est marqué par un point.}
+  \label{fig:test_get_eddy}
+\end{figure}
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{test_get_snapshot}
+  \caption{Test Get\_snapshot\_region\_1. Région 1, sans minimum
+    d'amplitude. Contours extérieurs et contours de vitesse
+    maximale. En rouge les anticyclones, en bleu les cyclones. Croix :
+    contour extérieur sans contour de maximum de vitesse
+    distinct. Cercles : contour extérieur et contour de maximum de
+    vitesse distinct. Les contours extérieurs 1, 2 et 8 passent près
+    de l'extremum et il n'y a donc pas de contour de vitesse maximale
+    correspondant. Les contours extérieurs 1 et 7 touchent le bord et
+    sont donc peut-être artificiellement sous-estimés.}
+  \label{fig:test_get_snapshot}
+\end{figure}
+
+Comparaison sur la région 1 avec le programme Matlab, pour les
+tourbillons 4 et 5. Les coordonnées des extremums sont bien les
+mêmes. Les SSH sur les contours extérieurs sont pratiquement les mêmes
+: différence de l'ordre de \np{0.1} mm. Du coup les contours
+extérieurs sont quasi-identiques. Les rayons des contours extérieurs
+sont presque les mêmes aussi : j'ai \np{87.7} km au lieu de \np{87.5}
+km et \np{90.80} km au lieu de \np{90.81} km. Les SSH sur les contours
+de vitesse max diffèrent de 2 cm environ. Cf. figure
+(\ref{fig:comparaison_Matlab}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{comparaison_201}
+  \includegraphics[width=\textwidth]{comparaison_202}
+  \caption{Comparaison avec le programme Matlab. Pour chacun des deux
+    tourbillons, le contour sans marqueur est à la valeur de SSH
+    donnée par le programme Matlab.}
+  \label{fig:comparaison_Matlab}
+\end{figure}
+J'obtiens comme vitesse maximale : $\np[m s^{-1}]{0.29}$ au lieu de
+$\np[m s^{-1}]{0.30}$ et $\np[m s^{-1}]{0.23}$ au lieu de
+$\np[m s^{-1}]{0.26}$.
+
+Test sur la région 3, sans minimum d'amplitude. 14 extremums sans
+contour extérieur (valid faux), dont 5 sur les bords et 11
+dégénérés. Avec minimum d'amplitude : 34 extremums pour lesquels
+valid est faux. Je vérifie que, pour chaque extremum, si valid
+est faux dans le cas sans minimum, valid reste faux dans le cas
+avec minimum. Il y a donc 20 extremums qui passent à valid faux à
+cause du minimum d'amplitude. Je vérifie que, pour ces 20 extremums,
+le contour extérieur dans le cas sans minimum est le même que dans le
+cas avec minimum. Entre le cas sans minimum et le cas avec minimum,
+seulement trois contours extérieurs grossissent (je le vois en
+comparant les fichiers outermost\_contour\_1.dbf) : dans le cas avec
+minimum, ils englobent un extremum d'amplitude insuffisante. Dans le
+cas avec minimum, 48 contours extérieurs sont calculés deux fois.
+
+Sur la région 3, avec un mm de minimum d'amplitude, le contour
+extérieur 168 a un diamètre de 27 points, environ 7 degrés en
+longitude. \`A 40° S, donc un diamètre de 600 km environ. On ne peut
+donc pas prendre max\_radius = 12. Ayant fait le test avec max\_radius
+= 12, je constate aussi que le temps d'exécution est très sensible à
+max\_radius : un facteur \np{2.7} environ entre max\_radius = 12 et
+max\_radius = 20. Le temps d'exécution semble donc proportionnel à
+max\_radius$^2$, ce qui peut se comprendre.
+
+\section{Recouvrements}
+
+\subsection{Divers}
+
+La zone intérieure à la ligne de vitesse maximale est le c\oe{}ur du
+tourbillon. On cherche les recouvrements entre ces c\oe{}urs, l'idée
+étant que les extremums à deux dates correspondent au même tourbillon
+si les c\oe{}urs se recouvrent.
+
+But du programme : produire un graphe. Les sommets du graphe sont les
+tourbillons. Les arcs relient les tourbillons dont les c\oe{}urs se
+recouvrent à deux dates successives.
+
+Le graphe est orienté acyclique. L'ordre topologique est tout trouvé :
+numérotation par dates successives. Le graphe non orienté associé est
+une forêt.
+
+Pour comparer les tourbillons à deux dates, il semble coûteux de
+n'utiliser que les listes de tourbillons aux deux dates. Cela implique
+une double boucle sur les listes de tourbillons, à l'intérieur d'une
+boucle sur les dates, soit potentiellement
+$8000 \times 8000 \times 5000$ comparaisons. Pour éviter cette
+comparaison aveugle, on peut conserver l'information sur la position
+des tourbillons dans des cartes aux différentes dates. La conservation
+en mémoire vive à 5000 dates semble aussi trop lourde. Il faut se
+restreindre aux cinq (\verb+max_delta+ + 1) dates nécessaires à la
+profondeur de vue.
+
+Lorsqu'on compare des tourbillons à distance temporelle delta (entre 1
+et \verb+max_delta+), on veut tester s'ils n'ont pas de prédécesseur ou
+successeur visible à une date intermédiaire. Il faut donc, si on crée
+un arc entre deux dates successives, ou delta arcs à distance delta
+(par interpolation de tourbillons invisibles), enregistrer ce fait en
+prévision du passage au delta supérieur. Mais d'un autre côté, on ne
+veut pas que le fait de mettre des arcs à l'étape delta ait une
+importance pour la suite de l'étape delta. (Et accessoirement, on ne
+veut pas que l'ordre dans lequel on passe les tourbillons à l'étape
+delta ait une importance.)  On est donc obligé d'enregistrer, pour
+chaque tourbillon, non seulement s'il a des prédécesseurs visibles et
+s'il a des successeurs visibles, mais aussi à quelle distance
+temporelle.
+
+On ne cherche des superpositions qu'entres tourbillons visibles, en
+excluant les tourbillons interpolés.
+
+On pourrait écrire l'algorithme principal en explicitant les
+\verb+get_snpashot+, en explicitant le cas m < \verb+n_proc+ et en supprimant la
+variable \verb+k_end_main_loop+. Ce serait mieux du point de vue de la
+performance, mais je pense que ce serait finalement moins clair.
+
+Cf. figure (\ref{fig:call_graph}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{call_graph}
+  \caption{Arbre des appels du programme de recherche de recouvrements.}
+  \label{fig:call_graph}
+\end{figure}
+
+Les procédures géométriques utilisées :
+\begin{itemize}
+\item polygon\_area\_2d (bibliothèque Geometry)
+\item gpc\_polygon\_clip\_f (bibliothèque GPC\_F)
+\end{itemize}
+Ces procédures sont appelées par :
+\begin{itemize}
+\item spherical\_polyline\_area
+\item \verb+successive_overlap+
+\end{itemize}
+Cf. discussion autour de l'équation~(\ref{eq:length_pi}).
+
+\subsection{Entrées et sorties}
+
+On n'enregistre pas de contours pour les tourbillons interpolés. Les
+tourbillons interpolés sont écrits avec les autres, en utilisant une
+forme NULL pour le polygone correspondant à un contour interpolé.
+
+La méta-donnée logique \og interpolé \fg{} dans
+\verb+extremum_$m.dbf+ est nécessaire parce qu'il est possible de ne
+pas trouver de contour extérieur autour d'un extremum détecté.
+
+Nous stockons le graphe au format liste d'arcs (en texte),
+c'est-à-dire, pour chaque sommet du graphe, une ligne :
+\begin{verbatim}
+indice date, indice tourbillon, indice date, indice tourbillon, poids
+\end{verbatim}
+Fichier
+\verb+edgelist_$m.txt+. L'estimation de la taille du fichier, en
+comptant 10 caractères par entier, 13 caractères par réel, \np{5e7}
+sommets et deux arcs par sommet en moyenne, est de 6 GiB.
+
+Problème des sommets isolés. L'information est déjà présente, de façon
+implicite, dans les fichiers évoqués ci-dessus. En effet, on a le
+nombre de tourbillons à chaque date d'après la liste des identifiants
+de tourbillons dans les shapefiles, et on a la liste des tourbillons
+non isolés d'après la liste des arcs. Mais la liste des sommets isolés
+n'est donc pas très facilement trouvée. On ajoute par conséquent les
+fichiers de sortie suivants :
+\begin{description}
+\item[number\_eddies\_\$m.csv] Une colonne indice de date, une colonne
+  nombre de tourbillons visibles, une colonne nombre de tourbillons
+  interpolés.
+\item[isolated\_nodes\_\$m.txt] Sur chaque ligne : indice de date, indice
+  de tourbillon. Peut être lu comme une simple liste
+  d'adjacence. Remarque : par construction, les tourbillons interpolés
+  ne sont jamais isolés.
+\end{description}
+
+Les entrées-sorties sont dans : l'algorithme principal directement
+(lecture de longitude et latitude), \verb+get_snapshot+ (lecture de
+SSH, u, v), \verb+successive_overlap+ (écriture dans
+\verb+edgelist_$m.txt+, arcs entre \verb+k_begin+(m) et
+\verb+k_end_main_loop+(m)), \verb+non_successive_overlap+ directement
+(écriture dans \verb+edgelist_$m.txt+), \verb+write_eddy+ (écriture
+dans \verb+extremum_$m+, \verb+outermost_contour_$m+,
+\verb+max_speed_contour_$m+), \verb+dispatch_snapshot+ directement
+(écriture dans \verb+isolated_nodes_$m+,
+\verb+number_eddies_$m+). Les processus écrivent dans des fichiers
+indicés par le numéro de processus donc deux processus ne peuvent
+écrire dans le même fichier. Cf. figure (\ref{fig:input_output}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{input_output}
+  \caption{Entrées et sorties du programme de recouvrement.}
+  \label{fig:input_output}
+\end{figure}
+
+\subsection{Déclarations pour l'algorithme principal}
+
+On suppose que le maillage est le même à toutes les dates et qu'il est
+régulier.
+
+Nous avons besoin que la composante extr\_map de snapshot soit étendue
+en longitude dans le cas d'un domaine périodique pour son utilisation
+par \verb+successive_overlap+.
+
+Les tailles ci-dessous sont calculées pour \verb+max_delta+ = 4 et une
+durée de 16 ans, soit environ 5000 dates.
+\begin{description}
+\item[max\_delta] Scalaire entier. Intervalle maximal d'indices de
+  dates auxquelles on cherche des recouvrements. Normalement 4.
+\item[flow] : Vecteur de type snapshot, de taille \verb+max_delta+ + 1,
+  chaque élément correspondant à une date. Environ 40 MiB.
+\item[dist\_lim] Scalaire entier constant, = 3° / résolution.
+\item[weight\_delta] Scalaire réel.
+\end{description}
+
+i désigne un indice de tourbillon, j un indice de position dans la
+fenêtre temporelle (entre 1 et \verb+max_delta+ + 1) et k un indice de
+date. Cf. figure (\ref{fig:window}) et algorithme
+(\ref{alg:main_sequential}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics{window}
+  \caption{Indice de position dans la fenêtre glissante et indice de
+    date.}
+  \label{fig:window}
+\end{figure}
+\begin{algorithm}[htbp]
+  \begin{algorithmic}
+    \STATE \COMMENT{\{max\_delta $\ge 1$; n\_dates $\ge$ max\_delta +
+      1\}} 
+
+    \STATE \COMMENT{prologue}
+
+    \STATE entrer(longitude, latitude)
+    \STATE \COMMENT{longitude et latitude dans l'ordre croissant} 
+    \FOR{k = 1 \TO max\_delta + 1}
+
+    \STATE appel de get\_snapshot(flow(k), k)
+
+    \ENDFOR
+    \FOR{k = 2 \TO max\_delta + 1}
+
+    \STATE appel de successive\_overlap(k, k, flow)
+    \ENDFOR
+    \FOR{delta = 2 \TO max\_delta}
+
+    \FOR{k = delta + 1 \TO max\_delta + 1}
+
+    \STATE appel de non\_successive\_overlap(k, k, delta, flow)
+
+    \ENDFOR
+    \ENDFOR
+
+    \STATE \COMMENT{boucle principale}
+
+    \FOR{k = max\_delta + 2 \TO n\_dates}
+
+    \STATE appel de dispatch\_snapshot(flow(1), k - max\_delta - 1)
+    \STATE flow(:max\_delta) = flow(2:)
+
+    \STATE appel de get\_snapshot(flow(max\_delta + 1), k)
+
+    \STATE appel de successive\_overlap(max\_delta + 1, k, flow)
+
+    \FOR{delta = 2 \TO max\_delta}
+
+    \STATE appel de non\_successive\_overlap(max\_delta + 1, k, delta,
+    flow)
+
+    \ENDFOR
+    \ENDFOR
+    \STATE \COMMENT{épilogue}
+
+    \FOR{j = 1 \TO max\_delta + 1}
+
+    \STATE appel de dispatch\_snapshot(flow(j), n\_dates - max\_delta - 1
+    + j)
+    \ENDFOR
+  \end{algorithmic}
+  \caption{Algorithme principal, séquentiel}
+  \label{alg:main_sequential}
+\end{algorithm}
+
+\subsection{Décomposition de domaine}
+
+Cf. figure (\ref{fig:processes}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=\textwidth]{processes}
+  \caption{Raccordement des processus. Les comparaisons représentées
+    ici en couleurs sont faites par le processus m. En rouge dans la
+    boucle principale, en bleu, vert et marron dans l'épilogue.}
+  \label{fig:processes}
+\end{figure}
+Il faut, pour m < \verb+n_proc+ :
+\begin{equation}
+  \label{eq:stitching}
+  \mathtt{k\_end}(m)
+  = \mathtt{k\_begin}(m + 1) + \mathtt{max\_delta} - 1
+\end{equation}
+
+k\_end - \verb+k_begin+ + 1 appels à \verb+get_snapshot+, pour les dates \verb+k_begin+
+à k\_end. k\_end - \verb+k_begin+ + 1 appels à \verb+dispatch_snapshot+, pour les
+dates \verb+k_begin+ à k\_end. Si m > 1 alors \verb+max_delta+ envois. Si m <
+\verb+n_proc+ alors \verb+max_delta+ réceptions.
+
+Nombre d'appels à \verb+get_snapshot+ avec lecture. Si m < \verb+n_proc+ : k\_end -
+\verb+k_begin+ - \verb+max_delta+ + 1. Si m = \verb+n_proc+ : k\_end - \verb+k_begin+ + 1.
+
+\verb+successive_overlap+ est appelé pour tous les indices k compris entre
+\verb+k_begin+ + 1 et \verb+k_end_main_loop+, seulement. Nombre d'appels à
+\verb+successive_overlap+ : \verb+k_end_main_loop+ - \verb+k_begin+. Donc si m <
+\verb+n_proc+ : k\_end - \verb+k_begin+ - \verb+max_delta+ + 1 appels. Si m = \verb+n_proc+ : k\_end
+- \verb+k_begin+ appels.
+
+Nombre d'appels à \verb+non_successive_overlap+ dans le prologue :
+\begin{equation*}
+  \sum_{\delta = 2} ^{\mathtt{max\_delta}}
+  (\mathtt{max\_delta} - \delta +1)
+  = \frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}
+\end{equation*}
+Dans la boucle principale : (\verb+max_delta+ - 1)(\verb+k_end_main_loop+ -
+\verb+k_begin+ - \verb+max_delta+). Dans l'épilogue, si m < \verb+n_proc+ :
+$\frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}$. Total si m
+< \verb+n_proc+ : (\verb+max_delta+ - 1)(k\_end - \verb+k_begin+ - \verb+max_delta+ +
+1). Total si m = \verb+n_proc+ : (\verb+max_delta+ - 1)(k\_end - \verb+k_begin+ -
+$\frac{\mathtt{max\_delta}}{2}$).
+
+Si k\_end - \verb+k_begin+ était le même pour m = \verb+n_proc+ et m <
+\verb+n_proc+, le processus m=\verb+n_proc+ aurait des plus grands nombres
+d'appels aux trois procédures. Notons $\Delta k$ = k\_end - \verb+k_begin+,
+considéré comme une fonction de m. Si, pour m < \verb+n_proc+ :
+\begin{equation}
+  \label{eq:charge}
+  \Delta k(m) \approx \Delta k(\mathtt{n\_proc}) + \mathtt{max\_delta}
+\end{equation}
+alors la charge devrait être à peu près équilibrée.
+
+On peut prendre, pour tout m :
+\begin{equation*}
+  \mathtt{k\_begin}(m)
+  = 1 + E[(m - 1) \mathtt{n\_dates} / \mathtt{n\_proc}]
+\end{equation*}
+(On prend la partie entière au lieu de nint pour pouvoir démontrer des
 encadrements.)  D'où, avec la contrainte (\ref{eq:stitching}), pour m
-< n\_proc :
+< \verb+n_proc+ :
 \begin{equation*}
   \mathtt{k\_end}(m)
   = E(m\ \mathtt{n\_dates} / \mathtt{n\_proc}) + \mathtt{max\_delta}
 \end{equation*}
 Alors on peut montrer que $\Delta k(m) - \Delta k(\mathtt{n\_proc})$
-vaut max\_delta ou max\_delta - 1. On a donc l'approximation
+vaut \verb+max_delta+ ou \verb+max_delta+ - 1. On a donc l'approximation
 (\ref{eq:charge}). Chaque processus fait aussi le même nombre de
 lectures de ssh, u, v.
 
-Pour n\_proc $\ne 1$ et m < n\_proc, le processus m + 1 ne doit pas
+Pour \verb+n_proc+ $\ne 1$ et m < \verb+n_proc+, le processus m + 1 ne doit pas
 relire une date déjà lue par le processus m. On doit donc avoir :
 \begin{equation}
   \label{eq:min_delta_k}
@@ -727,19 +1194,20 @@ Il suffit donc que :
 \end{equation*}
 
 L'équation (\ref{eq:stitching}) et l'inégalité (\ref{eq:min_delta_k})
-impliquent que, pour m < n\_proc :
+impliquent que, pour m < \verb+n_proc+ :
 \begin{equation}
   \label{eq:delta_k}
   \mathtt{k\_begin} + \mathtt{max\_delta}
   \le \mathtt{k\_end} - \mathtt{max\_delta}
 \end{equation}
-donc les appels de get\_snapshot dans le prologue font des
-lectures. Les appels de get\_snapshot dans l'épilogue font des
+donc les appels de \verb+get_snapshot+ dans le prologue font des
+lectures. Les appels de \verb+get_snapshot+ dans l'épilogue font des
 réceptions. L'inégalité (\ref{eq:delta_k}) implique en outre que, pour
-m < n\_proc, la boucle principale fait au moins une itération. Si m <
-n\_proc : tous les appels à get\_snapshot de la boucle principale sauf
-le dernier font des lectures ; le dernier appel à get\_snapshot de la
-boucle principale fait une réception.
+m < \verb+n_proc+, la boucle principale fait au moins une
+itération. Si m < \verb+n_proc+ : tous les appels à
+\verb+get_snapshot+ de la boucle principale sauf le dernier font des
+lectures ; le dernier appel à \verb+get_snapshot+ de la boucle
+principale fait une réception.
 
 Si $m \ge 2$ alors :
 \begin{equation*}
@@ -751,7 +1219,7 @@ Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
 \begin{table}[htbp]
   \centering
   \begin{tabular}{llll}
-    m & k\_begin & k\_end & k\_end\_main\_loop \\
+    m & \verb+k_begin+ & k\_end & \verb+k_end_main_loop+ \\
     \hline
     1 & 1 & 9 & 6 \\
     2 & 6 & 10 & 10
@@ -765,7 +1233,7 @@ Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
 \begin{table}[htbp]
   \centering
   \begin{tabular}{llll}
-    m & k\_begin & k\_end & k\_end\_main\_loop \\
+    m & \verb+k_begin+ & k\_end & \verb+k_end_main_loop+ \\
     \hline
     1 & 1 & 9 & 6 \\
     2 & 6 & 14 & 11 \\
@@ -784,16 +1252,17 @@ Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
   \label{fig:m3}
 \end{figure}
 
-Le dispatch\_snapshot dans la dernière boucle sur j n'est pas
-forcément une écriture. Exemple : n\_proc = 2, max\_delta = 4,
+Le \verb+dispatch_snapshot+ dans la dernière boucle sur j n'est pas
+forcément une écriture. Exemple : \verb+n_proc+ = 2, \verb+max_delta+ = 4,
 n\_dates = 10, m = 2.
 
-Si m = 1 alors write\_eddy est appelé pour les dates k\_begin à
-k\_end. Si m > 1 alors write\_eddy est appelé pour les dates
-k\_begin + max\_delta à k\_end. Avec l'équation (\ref{eq:stitching}),
-une date donnée est bien écrite par un seul processus.
+Si m = 1 alors \verb+write_eddy+ est appelé pour les dates
+\verb+k_begin+ à k\_end. Si m > 1 alors \verb+write_eddy+ est appelé
+pour les dates \verb+k_begin+ + \verb+max_delta+ à k\_end. Avec
+l'équation (\ref{eq:stitching}), une date donnée est bien écrite par
+un seul processus.
 
-\section{Algorithme principal, parallèle}
+\subsection{Algorithme principal, parallèle}
 
 \begin{algorithmic}[1]
   \STATE \COMMENT{\{max\_delta $\ge 1$; n\_dates $\ge$ max\_delta + 1;
@@ -891,237 +1360,44 @@ une date donnée est bien écrite par un seul processus.
   
   \FOR{j = 1 \TO max\_delta + 1}
 
-  \STATE appel de dispatch\_snapshot(flow(j), k\_end - max\_delta - 1
-  + j)
-
-  \ENDFOR
-\end{algorithmic}
-
-\section{Résumé du problème de géométrie algorithmique}
-
-\begin{algorithmic}
-  \STATE recherche de tous les extremums locaux
-
-  \FOR{chaque extremum local}
-
-  \STATE recherche du \og bon contour le plus extérieur\fg{}
-  \ENDFOR
-\end{algorithmic}
-
-\og Bon contour\fg{} : contient l'extremum considéré et aucun
-autre. \og Le plus extérieur\fg{} : plus grande différence de niveau
-avec l'extremum considéré, ou encore amplitude maximale.
-
-Recherche du bon contour le plus extérieur : par dichotomie.
-\begin{algorithmic}
-  \STATE $a_1 = 0$
-
-  \STATE $a_2$ = max ou min de SSH sur le domaine - SSH de l'extremum cible
-
-  \WHILE{$a_2 - a_1 >$ précision}
-
-  \STATE recherche du bon contour à l'amplitude
-  $\frac{a_1 + a_2}{2}$
-
-  \IF{bon contour trouvé}
-
-  \STATE $a_1 = \frac{a_1 + a_2}{2}$
-
-  \ELSE
-
-  \STATE $a_2 = \frac{a_1 + a_2}{2}$
-
-  \ENDIF
-
-  \ENDWHILE
-\end{algorithmic}
-
-Recherche d'un bon contour à une amplitude donnée. On cherche tous les
-contours à cette amplitude. Pour chaque contour, pour chaque extremum,
-on regarde si l'extremum est à l'intérieur du contour (test de point
-intérieur à polygone).
-
-Les procédures géométriques utilisées :
-\begin{itemize}
-\item polygon\_area\_2d (bibliothèque Geometry)
-\item polygon\_contains\_point (bibliothèque Geometry)
-\item find\_contours\_reg\_grid (bibliothèque Contour\_531)
-\item gpc\_polygon\_clip\_f (bibliothèque GPC\_F)
-\end{itemize}
-travaillent dans le plan
-et non sur la sphère. Ces procédures sont appelées par :
-\begin{itemize}
-\item spherical\_polyline\_area
-\item inside\_4
-\item good\_contour
-\item successive\_overlap
-\end{itemize}
-Dans le plan, par deux points passe un unique
-segment, et la donnée d'un triangle est équivalente à la donnée de
-l'ensemble de ses trois sommets. Sur la sphère, par deux points
-passent deux arcs de grand cercle, et plusieurs triangles peuvent être
-associés à un ensemble de trois sommets. Par ailleurs, les procédures
-géométriques dans le plan ne peuvent pas admettre en arguments
-d'entrée des longitudes définies modulo $2 \pi$. Pour que les valeurs
-de longitudes en entrée des procédures géométriques soient contraintes
-sans ambiguïté, nous imposons que tout domaine de longitudes en entrée
-des procédures géométriques soit de longueur strictement inférieure à
-$\pi$.
-\begin{equation}
-  \label{eq:length_pi}
-  \forall (\lambda_1, \lambda_2) \in \mathbb{R}^2, \forall k \in \mathbb{Z},
-  |\lambda_2 - \lambda_1| < \pi
-  \Rightarrow |\lambda_2 + 2 k \pi - \lambda_1| \ge \pi
-\end{equation}
-(Toute modulation de $2 \pi$ d'une longitude viole la contrainte.) Les
-procédures qui appellent les procédures géométriques doivent se
-charger d'ajouter le multiple de $2 \pi$ pour ramener les longitudes
-dans un intervalle de longueur strictement inférieure à $\pi$. La
-contrainte de domaine de longitudes s'applique aux procédures
-get\_1\_outerm, set\_max\_speed et successive\_overlap.
-
-\section{Sous-algorithmes}
-
-\subsection{get\_1\_outerm}
-
-Au lieu d'initialiser l'amplitude maximale à 4 m plus du bruit, comme
-dans le programme pour Matlab :
-\begin{algorithmic}
-   \STATE appel de random\_number(level\_try)
-
-   \STATE level\_try = e\%ssh\_extr $\pm$ [4 + (2 * level\_try -
-   1) * precision] \COMMENT{un peu de bruit sur l'amplitude, autour de
-     4}
-\end{algorithmic}
-et de supposer qu'il n'y a pas de bon contour avec cette amplitude, on
-cherche une amplitude maximale possible. Dans le cas le plus simple où
-noise\_around est faux, on purrait parcourir les points dans une
-direction jusqu'à ce que le sens de variation de SSH change. Dans mes
-tests, le nombre d'itérations dans la recherche par dichotomie est
-typiquement une dizaine, pour une précision de \np{0.5} mm. En
-considérant simplement l'extremum de ssh sur la fenêtre, la taille de
-l'intervalle initial [innermost\_level, level\_try] est plus
-grande. Sur mes tests, elle est plus grande typiquement d'un facteur 2
-environ seulement, donc le nombre d'itérations change peu.
-
-Pour l'amplitude minimale, si je prends simplement l'amplitude
-correspondant à innermost\_level alors l'algorithme échoue lorsque
-l'extremum est à un point du bord et que le point à innermost\_level
-est au bord. Dans ce cas, good\_contour échoue parce que Contour\_531
-ne ferme pas les contours passant par le bord. J'avais essayé, pour
-éviter le problème, de prendre une amplitude minimale légèrement
-inférieure mais cela pose des problèmes de cohérence dans
-l'algorithme.
-
-\subsection{get\_snapshot}
-
-Pour chaque tourbillon du résultat, le contour le plus extérieur doit
-être fermé, doit contenir l'extremum de SSH du tourbillon considéré et
-ne doit contenir aucun autre extremum de SSH avec une amplitude
-significative. Amplitude significative : supérieure à un certain
-minimum, entre l'extremum de SSH et le contour le plus
-lointain.
-
-Dans une première version qui n'envisageait pas d'amplitude minimum,
-j'avais une procédure get\_eddy qui définissait un tourbillon, et qui
-trouvait donc le contour le plus extérieur et le contour de maximum de
-vitesse. Mais avec la possibilité d'une amplitude minimale, il faut
-d'abord se préoccuper de tous les contours extérieurs, qui ne sont
-plus indépendants les uns des autres, pour ne pas chercher inutilement
-des contours de maximum de vitesse. Dans le nouvel algorithme, il n'y
-a donc plus de procédure get\_eddy.
-
-Si le domaine considéré n'est pas global, comment traiter un
-tourbillon au bord du domaine ? La taille du contour le plus extérieur
-peut être artificiellement diminuée. (Ce n'est pas le même cas qu'un
-tourbillon près d'une côte.) C'est-à-dire qu'on pourrait trouver un
-contour plus extérieur avec un domaine élargi. Il pourrait être
-intéressant d'éliminer un tel tourbillon. Comment le détecter ? Le
-fait d'aller jusqu'au bord du domaine dans outermost\_possible\_level
-n'est pas un critère suffisant, parce que outermost\_possible\_level
-ne teste que la direction est. Il faudrait tester après l'appel de
-get\_1\_outerm si le cadre le plus extérieur touche le bord du
-domaine.
-
-\subsection{good\_contour}
-
-Trouve, si elle existe, une ligne de niveau à un niveau donné,
-contenant un point donné et ne contenant pas un ensemble de points
-donnés.  Le résultat est un polygone à zéro point s'il n'existe pas de
-ligne de niveau satisfaisante. Correspond à une partie de la fonction
-extraction\_eddies dans le programme pour Matlab.
-
-\subsection{interpolate\_eddy}
-
-Tous les arguments sont des données. e1, e2 : scalaires de type
-eddy. j1, j2, indices de position dans la fenêtre correspondant à e1,
-e2. j : indice de position dans la fenêtre auquel on veut
-interpoler. k, i : indices de date et de tourbillon. Cf. algorithme
-(\ref{alg:interpolate_eddy}).
-\begin{algorithm}[htbp]
-  \begin{algorithmic}
-    \STATE e = interpolation entre e1 et e2
-
-    \STATE appel de write\_eddy(e, k, i)
-  \end{algorithmic}
-  \caption{Sous-algorithme interpolate\_eddy(e1, e2, j1, j2, j, k, i)}
-  \label{alg:interpolate_eddy}
-\end{algorithm}
-On peut mettre le champ valid à faux pour un tourbillon interpolé,
-en cohérence avec un champ out\_cont vide. Il faut alors un
-champ interpolated pour distinguer les tourbillons interpolés de ceux
-pour lesquels un extremum a été détecté mais sans obtenir de contour
-extérieur.
-
-\subsection{local\_extrema}
-
-a(ind\_extr(1, k), ind\_extr(2, k)) est un extremum de a par rapport
-aux huit points voisins. extr\_map vaut 0 en un point où il n'y a pas
-d'extremum, le numéro de l'extremum dans ind\_extr en un point
-extremum. Accélération possible, au prix de tableaux supplémentaires,
-en utilisant le fait que deux maximums ne peuvent être côte à côte et
-deux minimums ne peuvent être côte à côte.
-
-J'avais d'abord codé le sens de l'extremum (maximum ou minimum) par le
-signe dans extr\_map. Mais c'est plus clair et plus simple sans ce
-codage. J'ajoute donc un argument vecteur cyclone à local\_extrema,
-qui prend moins de place qu'une carte supplémentaire. Pour le sens de
-l'extremum, je n'ai pas besoin que l'information soit sur une carte.
-
-Cf. aussi fontion \verb+skimage.feature.peak_local_max+ dans le paquet
-Python scikit-image. La fonction applique
-\verb+scipy.ndimage.maximum_filter+ et compare le résultat au tableau
-de départ. \verb+scipy.ndimage.maximum_filter+ applique
-\verb+scipy.ndimage.maximum_filter1d+ dans chaque dimension. Cf. aussi
-Douglas (1996 k0968).
-
-Cf. figure (\ref{fig:periodicity}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics{periodicity}
-  \caption{Cas de périodicité dans local\_extrema. Le tableau
-    représenté est extr\_map. Pour définir la case (2, j), il faut
-    avoir copié la case (m - 1, j - 1) (qui peut être non nulle). Pour
-    définir la case (m - 1, j), il faut avoir copié la case (2, j)
-    (qui peut être non nulle). (La case (2, j - 1) doit aussi avoir
-    été copiée mais cela a été fait au passage sur la colonne
-    précédente.)}
-  \label{fig:periodicity}
-\end{figure}
-Cas de périodicité : field(i + m - 2, :) = field(i, :).
-
-\subsection{max\_speed\_contour\_ssh}
+  \STATE appel de dispatch\_snapshot(flow(j), k\_end - max\_delta - 1
+  + j)
 
-Il manque des tests pour lesquels direction vaut 1 et 4.
+  \ENDFOR
+\end{algorithmic}
+
+\subsection{Sous-algorithmes}
+
+\subsubsection{interpolate\_eddy}
+
+Tous les arguments sont des données. e1, e2 : scalaires de type
+eddy. j1, j2, indices de position dans la fenêtre correspondant à e1,
+e2. j : indice de position dans la fenêtre auquel on veut
+interpoler. k, i : indices de date et de tourbillon. Cf. algorithme
+(\ref{alg:interpolate_eddy}).
+\begin{algorithm}[htbp]
+  \begin{algorithmic}
+    \STATE e = interpolation entre e1 et e2
 
-\subsection{non\_successive\_overlap}
+    \STATE appel de write\_eddy(e, k, i)
+  \end{algorithmic}
+  \caption{Sous-algorithme interpolate\_eddy(e1, e2, j1, j2, j, k, i)}
+  \label{alg:interpolate_eddy}
+\end{algorithm}
+On peut mettre le champ valid à faux pour un tourbillon interpolé,
+en cohérence avec un champ \verb+out_cont+ vide. Il faut alors un
+champ interpolated pour distinguer les tourbillons interpolés de ceux
+pour lesquels un extremum a été détecté mais sans obtenir de contour
+extérieur.
+
+\subsubsection{non\_successive\_overlap}
 
 j, k, delta $\in \{2, \dots, \mathtt{max\_delta}\}$ : scalaires
 entiers, données. flow : donnée-résultat. Cherche les arcs entre
 flow(j - delta) et flow(j), correspondant aux dates k - delta et
 k. Écrit ces arcs dans \verb+graph.txt+, écrit les tourbillons
 interpolés dans les shapefiles via interpolate\_eddy et
-write\_eddy. Met à jour flow(j - delta)\%list\_vis\%delta\_out,
+\verb+write_eddy+. Met à jour flow(j - delta)\%list\_vis\%delta\_out,
 flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
 (\ref{alg:non_successive_overlapping}).
 \begin{algorithm}[htbp]
@@ -1200,152 +1476,20 @@ flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
   \label{alg:non_successive_overlapping}
 \end{algorithm}
 
-\subsection{set\_all\_outerm}
-
-Notons ici :
-\begin{align*}
-  & u = \mathtt{urc(1)} \\
-  & l = \mathtt{llc(1)} \\
-  & s = \mathtt{step(1)} \\
-  & e = \mathtt{s\%list\_vis(i)\%coord\_extr(1)} \\
-  & w = \mathtt{corner\_window(1)} \\
-  & m = \mathtt{max\_radius(1)}
-\end{align*}
-Pour que la contrainte de domaine de longitudes dans get\_1\_outerm
-soit satisfaite, il faut que, dans set\_all\_outerm, au moment de
-l'appel à get\_1\_outerm :
-\begin{equation}
-  \label{eq:max_radius}
-  (u - l) s < \pi  
-\end{equation}
-et que $e$ et outside\_points(1, :) soient dans l'intervalle
-$[w, w + (u - l) s]$. Pour que l'inégalité (\ref{eq:max_radius}) soit
-satisfaite, il suffit que :
-\begin{equation*}
-  2 m s < \pi
-\end{equation*}
-ce que nous imposons en amont de l'appel de set\_all\_outerm. Si
-periodic est vrai alors :
-\begin{equation*}
-  [w, w + (u - l) s] = [e - m s, e + m s]
-\end{equation*}
-donc $e$ appartient bien à cet intervalle. Si periodic est faux,
-$[w, w + (u - l) s]$ peut être plus réduit que $[e - m s, e + m s]$
-mais nous savons aussi que $e$ est dans l'intervalle [corner(1),
-corner(1) + (nlon - 1) s], donc $e$ est bien toujours dans
-$[w, w + (u - l) s]$.
-
-Il reste le problème de outside\_points(1, :). outside\_points(1, :)
-contient des coordonnées prises dans s\%list\_vis\%coord\_extr(1). Si
-le domaine est périodique et si l'extremum cible est proche en
-longitude de corner(1) par exemple alors outside\_points(1, :) peut
-contenir des longitudes proches de corner(1) + (nlon - 1) s. Il faut
-leur retrancher $2 \pi$. Nous posons le problème plus généralement de
-la façon suivante. Soit $\lambda \in \mathbb{R}$. On cherche $k \in
-\mathbb{Z}$ tel que :
-\begin{equation*}
-  \lambda + 2 k \pi \in [w, w + (u - l) s]
-\end{equation*}
-Supposons que $k$ existe (c'est le cas si $\lambda$ est une longitude
-retournée par nearby\_extr). Alors $k$ est unique (cf. propriété
-(\ref{eq:length_pi})). On a :
-\begin{equation*}
-  \lambda + 2 k \pi \in [w, w + 2 \pi[
-\end{equation*}
-C'est-à-dire :
-\begin{equation*}
-  k - 1 < \frac{w - \lambda}{2 \pi} \le k
-\end{equation*}
-Ou encore :
-\begin{equation*}
-  k = \left\lceil \frac{w - \lambda}{2 \pi} \right\rceil
-\end{equation*}
-
-
-\subsection{set\_max\_speed}
-
-Notons $d$ la distance dans l'espace des indices entre l'extremum
-cible et le countour extérieur. Par définition de la partie entière,
-les points à distance $E(d)$ de l'extremum cible sont à l'intérieur du
-contour. Certains points à distance $E(d) + 1$ peuvent être à
-l'extérieur. D'où : radius = $E(d) + 1$.
-
-Dans le cas où radius vaut 1, je calcule max\_speed. Ce calcul ne sera
-utile que dans weight si le tourbillon a une intersection avec avec un
-tourbillon à une autre date. Donc, pour certains tourbillons, ce
-calcul est inutile. Difficile d'éviter ces calculs inutiles. Il
-faudrait les reporter à weight mais alors il faudrait avoir accès dans
-weight à u et v, c'est-à-dire qu'il faudrait conserver u, v à
-plusieurs dates en dehors de get\_snpashot. Cela remettrait en
-question tout l'algorithme.
-
-\subsection{spherical\_polyline\_area}
-
-La procédure de calcul des contours donne une suite de sommets. Les
-courbes reliant ces sommets ne sont en fait pas imposés
-intrinsèquement par cette procédure. Par exemple, il n'est ni plus ni
-moins juste de supposer que ces courbes sont des portions de grands
-cercles ou de supposer qu'elles donnent des segments en projection
-plate carrée. Néanmoins, dans le reste du programme, pour décider si
-un contour entoure un extremum, ou pour calculer l'intersection de
-deux tourbillons, nous supposons que les contours sont des polygones
-en $(\lambda, \phi)$ (c'est-à-dire que les courbes reliant les sommets
-donnent des segments en projection plate carrée).
-
-Un domaine $S$ quelconque à la surface de la Terre a pour aire :
-\begin{equation*}
-  a^2  \iint_S \cos \phi\ \ud \phi\ \ud \lambda
-  = a^2  \iint_S \ud \lambda\ \ud(\sin \phi)
-\end{equation*}
-Nous avons une procédure qui calcule la surface $\iint \ud x\ \ud y$
-d'un polygone en $(x, y)$. Nous calculons la surface du polygone en
-$(\lambda, \sin \phi)$ passant par les sommets du
-contour. C'est-à-dire que nous supposons que les courbes reliant les
-somments donnent des segments en projection conforme cylindrique.
-
-Rémi fait un autre calcul : la surface du polygone en
-$(\lambda \cos \phi, \phi)$ passant par les sommets du contour. Le
-changement de variable :
-\begin{align*}
-  & x = \lambda \cos \phi \\
-  & y = \phi
-\end{align*}
-a pour jacobien $\cos \phi$ donc on a bien :
-\begin{equation*}
-  \iint \cos \phi\ \ud \phi\ \ud \lambda = \iint \ud x\ \ud y
-\end{equation*}
-Les deux surfaces obtenues, celle du polygone en
-$(\lambda, \sin \phi)$ et celle du polygone en
-$(\lambda \cos \phi, \phi)$ sont très proches si les sommets sont très
-proches. C'est le cas puisque les sommets sont séparés au maximum de
-\np{0.25}° en longitude et latitude. On peut aussi le voir en
-considérant la formule de calcul de surface du polygone :
-\begin{align*}
-  S & = \frac{1}{2} \sum_i x_i (y_{i + 1} - y_{i - 1}) \\
-    & =
-      \frac{1}{2} \sum_i \lambda_i \cos \phi_i (\phi_{i + 1} - \phi_{i - 1}) \\
-    & \approx \frac{1}{2} \sum_i \lambda_i (\sin \phi_{i + 1} - \sin \phi_{i - 1})
-\end{align*}
-si les valeurs successives des $\phi_i$ sont proches.
-
-\subsection{successive\_overlap}
+\subsubsection{successive\_overlap}
 
-j, k, scalaires entiers, données. flow : donnée-résultat. Cherche les
-arcs entre flow(j - 1) et flow(j), correspondant aux dates k - 1 et
-k. Écrit ces arcs dans \verb+graph.txt+. Met à jour flow(j -
-1)\%list\_vis\%delta\_out et
-flow(j)\%list\_vis\%delta\_in. Accélération possible en ne comparant à
-i1 que les trois tourbillons les plus proches (distance entre les
-extremums), comme dans la fonction \verb+Eddy_Intersection+ du
-programme pour Matlab.
+Accélération possible en ne comparant à i1 que les trois tourbillons
+les plus proches (distance entre les extremums), comme dans la
+fonction \verb+Eddy_Intersection+ du programme pour Matlab.
 
 On ne peut éviter la copie de polylines due à l'alternative contour
 extérieur -- contour de vitesse maximale. Pour ne pas faire la copie
 dans polygon\_1 et polygon\_2, il faudrait faire des copies de
-out\_cont dans speed\_cont dans get\_snapshot, on ne gagnerait rien en
-quantité de copie, et on perdrait en mémoire vive occupée.
+\verb+out_cont+ dans \verb+speed_cont+ dans \verb+get_snapshot+, on ne
+gagnerait rien en quantité de copie, et on perdrait en mémoire vive
+occupée.
 
-\subsection{weight}
+\subsubsection{weight}
 
 D'autant plus proche de 0 que les tourbillons sont
 ressemblants. Comment calculer ce poids ? Comment calculer le poids
@@ -1355,115 +1499,4 @@ intervenir les tourbillons interpolés le poids associé aux tourbillons
 visibles entre lesquels on interpole. Normalement, ce choix doit
 donner un poids plus faible à ces arcs.
 
-\subsection{write\_eddy}
-
-On pourrait passer en arguments d'entrée les numéros des
-champs dans les fichiers DBF mais cela ferait 13
-arguments.
-
-\section{Tests}
-
-Pour faire des tests : données au 1\ier{} janvier 2006. Tests sur
-différents domaines. Cf. figure (\ref{fig:regions}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{regions}
-  \caption{Régions pour les tests.}
-  \label{fig:regions}
-\end{figure}
-Régions 1, 2 et 3 : 29 $\times$ 17, 57 $\times$ 33; 120 $\times$ 120,
-contenant respectivement 8, 27 et 216 extremums.
-
-Pour la région 1, coin en :
-\begin{equation*}
-  \lambda = \np{5.125}^\circ, \phi = -\np{36.375}^\circ
-\end{equation*}
-soit $(\np[rad]{0.0894}, - \np[rad]{0.6349})$. Indices base 1 du coin
-: 21, 215. Pas de grille : \np{0.25}°, soit \np[rad]{0.0044}. 5
-minimums et 3 maximums dans ce cadre. Cf. figure (\ref{fig:extrema}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{extrema}
-  \caption{Région 1, lignes de niveau de SSH (sans signification
-    particulière). Maximums : points rouges, minimums : points bleus.}
-  \label{fig:extrema}
-\end{figure}
-Minimum cible pour le test \verb+Get_1_outerm+ en :
-\begin{equation*}
-  \lambda = \np{9.625}^\circ, \phi = - \np{33.875}^\circ
-\end{equation*}
-Indices base 1 de ce minimum : 39, 225. SSH en ce mimimum :
-\np{0.2759}. Extremum numéro 6 dans ce cadre. Contour le plus
-intérieur : SSH = \np{0.293300003} m. Contour le plus extérieur : SSH
-= \np{0.460744} m. Contour de vitesse maximale : SSH = \np{0.3815}
-m. Cf. figures (\ref{fig:test_get_eddy}) et (\ref{fig:test_get_snapshot}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{test_get_eddy}
-  \caption{Test Set\_max\_speed. Projection stéréographique centrée
-    sur un minimum de SSH. Champ de vitesse, contour de SSH le plus
-    extérieur en bleu, contour de SSH de vitesse maximale en vert. Le
-    minimum de SSH est marqué par un point.}
-  \label{fig:test_get_eddy}
-\end{figure}
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{test_get_snapshot}
-  \caption{Test Get\_snapshot\_region\_1. Région 1, sans minimum
-    d'amplitude. Contours extérieurs et contours de vitesse
-    maximale. En rouge les anticyclones, en bleu les cyclones. Croix :
-    contour extérieur sans contour de maximum de vitesse
-    distinct. Cercles : contour extérieur et contour de maximum de
-    vitesse distinct. Les contours extérieurs 1, 2 et 8 passent près
-    de l'extremum et il n'y a donc pas de contour de vitesse maximale
-    correspondant. Les contours extérieurs 1 et 7 touchent le bord et
-    sont donc peut-être artificiellement sous-estimés.}
-  \label{fig:test_get_snapshot}
-\end{figure}
-
-Comparaison sur la région 1 avec le programme Matlab, pour les
-tourbillons 4 et 5. Les coordonnées des extremums sont bien les
-mêmes. Les SSH sur les contours extérieurs sont pratiquement les mêmes
-: différence de l'ordre de \np{0.1} mm. Du coup les contours
-extérieurs sont quasi-identiques. Les rayons des contours extérieurs
-sont presque les mêmes aussi : j'ai \np{87.7} km au lieu de \np{87.5}
-km et \np{90.80} km au lieu de \np{90.81} km. Les SSH sur les contours
-de vitesse max diffèrent de 2 cm environ. Cf. figure
-(\ref{fig:comparaison_Matlab}).
-\begin{figure}[htbp]
-  \centering
-  \includegraphics[width=\textwidth]{comparaison_201}
-  \includegraphics[width=\textwidth]{comparaison_202}
-  \caption{Comparaison avec le programme Matlab. Pour chacun des deux
-    tourbillons, le contour sans marqueur est à la valeur de SSH
-    donnée par le programme Matlab.}
-  \label{fig:comparaison_Matlab}
-\end{figure}
-J'obtiens comme vitesse maximale : $\np[m s^{-1}]{0.29}$ au lieu de
-$\np[m s^{-1}]{0.30}$ et $\np[m s^{-1}]{0.23}$ au lieu de
-$\np[m s^{-1}]{0.26}$.
-
-Test sur la région 3, sans minimum d'amplitude. 14 extremums sans
-contour extérieur (valid faux), dont 5 sur les bords et 11
-dégénérés. Avec minimum d'amplitude : 34 extremums pour lesquels
-valid est faux. Je vérifie que, pour chaque extremum, si valid
-est faux dans le cas sans minimum, valid reste faux dans le cas
-avec minimum. Il y a donc 20 extremums qui passent à valid faux à
-cause du minimum d'amplitude. Je vérifie que, pour ces 20 extremums,
-le contour extérieur dans le cas sans minimum est le même que dans le
-cas avec minimum. Entre le cas sans minimum et le cas avec minimum,
-seulement trois contours extérieurs grossissent (je le vois en
-comparant les fichiers outermost\_contour\_1.dbf) : dans le cas avec
-minimum, ils englobent un extremum d'amplitude insuffisante. Dans le
-cas avec minimum, 48 contours extérieurs sont calculés deux fois.
-
-Sur la région 3, avec un mm de minimum d'amplitude, le contour
-extérieur 168 a un diamètre de 27 points, environ 7 degrés en
-longitude. \`A 40° S, donc un diamètre de 600 km environ. On ne peut
-donc pas prendre max\_radius = 12. Ayant fait le test avec max\_radius
-= 12, je constate aussi que le temps d'exécution est très sensible à
-max\_radius : un facteur \np{2.7} environ entre max\_radius = 12 et
-max\_radius = 20. Le temps d'exécution semble donc proportionnel à
-max\_radius$^2$, ce qui peut se comprendre.
-
 \end{document}
diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f
index 63291acc266b420da209b6a34c75309ae2be5000..a5443bf4f9255b00ec707bcb5ad78d6009685c9b 100644
--- a/Tests/test_get_snapshot.f
+++ b/Tests/test_get_snapshot.f
@@ -5,17 +5,16 @@ program test_get_snapshot
   use, intrinsic:: ISO_FORTRAN_ENV
 
   ! Libraries:
-  use jumble, only: new_unit
   use netcdf, only: nf90_nowrite
   use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
        nf95_get_var, nf95_close
   use nr_util, only: assert, deg_to_rad, twopi, pi
   use shapelib, only: shpfileobject, shpclose
-  
+
   use derived_types, only: snapshot
-  use dispatch_snapshot_m, only: dispatch_snapshot
   use get_snapshot_m, only: get_snapshot
   use init_shapefiles_m, only: init_shapefiles
+  use write_eddy_m, only: write_eddy
 
   implicit none
 
@@ -23,8 +22,7 @@ program test_get_snapshot
   TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
   TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
   TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_$m
-
-  integer unit_isolated, unit_number_eddies
+  integer i
 
   real:: min_amp = 0.
   ! minimum amplitude of ssh, between extremum and outermost contour,
@@ -32,7 +30,7 @@ program test_get_snapshot
 
   real:: min_radius = 20.
   ! minimum radius of the equal-area disk for an outermost contour, in km
-  
+
   integer ncid, dimid, varid, nlon, nlat
   real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
   real step(2) ! longitude and latitude steps, in rad
@@ -67,7 +65,7 @@ program test_get_snapshot
   call nf95_get_var(ncid, varid, lat_max, start = [nlat])
 
   call nf95_close(ncid)
-  
+
   call assert(lon_max > lon_min .and. lat_max > lat_min, &
        "test_get_snapshot coordinate order")
   step = [(lon_max - lon_min) / (nlon - 1), (lat_max - lat_min) / (nlat - 1)]  &
@@ -85,32 +83,19 @@ program test_get_snapshot
        step = step, periodic = periodic)
 
   call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
-  
-  call new_unit(unit_isolated)
-  open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
-       action = "write")
-  write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
-
-  call new_unit(unit_number_eddies)
-  open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
-       action = "write")
-  write(unit_number_eddies, fmt = *) '"date index" ' &
-       // '"number of visible extrema" "number of interpolated eddies"'
-  ! (title line)
-
-  call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
-       unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
-       max_delta = 4)
 
+  ! Write snapshot:
+  do i = 1, s%number_vis_extr
+     call write_eddy(s%list_vis(i), 1, i, hshp_extremum, hshp_outermost, &
+          hshp_max_speed)
+  end do
+
+  print *, "Number of extrema:", s%number_vis_extr
   CALL shpclose(hshp_extremum)
   print *, 'Created shapefile "extremum_1".'
   CALL shpclose(hshp_outermost)
   print *, 'Created shapefile "outermost_contour_1".'
   CALL shpclose(hshp_max_speed)
   print *, 'Created shapefile "max_speed_contour_1".'
-  close(unit_isolated)
-  print *, 'Created "isolated_nodes_1.csv".'
-  close(unit_number_eddies)
-  print *, 'Created "number_eddies_1.csv".'
 
 end program test_get_snapshot
diff --git a/Tests/test_read_snapshot.f b/Tests/test_read_snapshot.f
index 52c27cd83f021e51408ec141765b2fad913f4464..857753c6db0f5fde97cdcb7cef2d10a53368ea3d 100644
--- a/Tests/test_read_snapshot.f
+++ b/Tests/test_read_snapshot.f
@@ -2,22 +2,23 @@ program test_read_snapshot
 
   use, intrinsic:: ISO_FORTRAN_ENV
 
-  use derived_types, only: snapshot
-  use dispatch_snapshot_m, only: dispatch_snapshot
-  use init_shapefiles_m, only: init_shapefiles
-  use jumble, only: new_unit
+  ! Libraries:
   use nr_util, only: deg_to_rad
-  use read_snapshot_m, only: read_snapshot
   use shapelib, only: shpfileobject, shpclose
   use shapelib_03, only: shp_open_03
 
+  use derived_types, only: snapshot
+  use init_shapefiles_m, only: init_shapefiles
+  use read_snapshot_m, only: read_snapshot
+  use write_eddy_m, only: write_eddy
+
   implicit none
 
   type(snapshot) s
   TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
   TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
   TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
-  integer unit_isolated, unit_number_eddies, i
+  integer i
 
   real:: corner(2) = [0.125, - 59.875]
   ! longitude and latitude of the corner of the whole grid, in degrees
@@ -45,32 +46,19 @@ program test_read_snapshot
 
   call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
 
-  call new_unit(unit_isolated)
-  open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
-       action = "write")
-  write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
-
-  call new_unit(unit_number_eddies)
-  open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
-       action = "write")
-  write(unit_number_eddies, fmt = *) '"date index" ' &
-       // '"number of visible eddies" "number of interpolated eddies"'
-  ! (title line)
-
-  call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
-       unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
-       max_delta = 4)
+  ! Write snapshot:
+  do i = 1, s%number_vis_extr
+     call write_eddy(s%list_vis(i), 1, i, hshp_extremum, hshp_outermost, &
+          hshp_max_speed)
+  end do
 
+  print *, "Number of extrema:", s%number_vis_extr
   CALL shpclose(hshp_extremum)
   print *, 'Created shapefile "extremum_1".'
   CALL shpclose(hshp_outermost)
   print *, 'Created shapefile "outermost_contour_1".'
   CALL shpclose(hshp_max_speed)
   print *, 'Created shapefile "max_speed_contour_1".'
-  close(unit_isolated)
-  print *, 'Created "isolated_nodes_1.csv".'
-  close(unit_number_eddies)
-  print *, 'Created "number_eddies_1.csv".'
 
   print *, "s%ind_extr:"
 
diff --git a/depend.mk b/depend.mk
index f8a9e17951d32e48e03ec1e589d9d0f8d1a6b385..67e5337c7080c3bd3709772f0247c8691e36b78e 100644
--- a/depend.mk
+++ b/depend.mk
@@ -1,4 +1,4 @@
-dispatch_snapshot.o : write_eddy.o send_snapshot.o derived_types.o 
+dispatch_snapshot.o : send_snapshot.o derived_types.o 
 get_1_outerm.o : spherical_polyline_area.o good_contour.o derived_types.o 
 get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o nearby_extr.o get_var.o derived_types.o 
 nearby_extr.o : derived_types.o 
@@ -11,11 +11,11 @@ set_max_speed.o : spherical_polyline_area.o mean_speed.o max_speed_contour_ssh.o
 spherical_polygon_area.o : spherical_polyline_area.o 
 successive_overlap.o : weight.o spherical_polyline_area.o spherical_polygon_area.o derived_types.o 
 test_get_1_outerm.o : get_1_outerm.o derived_types.o 
-test_get_snapshot.o : init_shapefiles.o get_snapshot.o dispatch_snapshot.o derived_types.o 
+test_get_snapshot.o : write_eddy.o init_shapefiles.o get_snapshot.o derived_types.o 
 test_local_extrema.o : local_extrema.o 
 test_nearby_extr.o : nearby_extr.o derived_types.o 
 test_read_eddy.o : write_eddy.o read_field_indices.o read_eddy.o init_shapefiles.o derived_types.o 
-test_read_snapshot.o : read_snapshot.o init_shapefiles.o dispatch_snapshot.o derived_types.o 
+test_read_snapshot.o : write_eddy.o read_snapshot.o init_shapefiles.o derived_types.o 
 test_set_all_outerm.o : set_all_outerm.o get_var.o derived_types.o 
 test_set_max_speed.o : set_max_speed.o derived_types.o 
 test_spherical_polygon_area.o : spherical_polygon_area.o 
diff --git a/dispatch_snapshot.f b/dispatch_snapshot.f
index a4004a01554f3e136bfaa9cd80dd58399395d61e..759ffaec2c58b061410d28d71638029f86d47cc3 100644
--- a/dispatch_snapshot.f
+++ b/dispatch_snapshot.f
@@ -4,27 +4,14 @@ module dispatch_snapshot_m
 
 contains
 
-  subroutine dispatch_snapshot(s, hshp_extremum, hshp_outermost, &
-       hshp_max_speed, unit_isolated, unit_number_eddies, k, m, k_begin, &
-       max_delta)
-
-    ! Library:
-    use shapelib, only: shpfileobject
+  subroutine dispatch_snapshot(s, unit_isolated, unit_number_eddies, k, m, &
+       k_begin, max_delta)
 
     use derived_types, only: snapshot
     use send_snapshot_m, only: send_snapshot
-    use write_eddy_m, only: write_eddy
 
     type(snapshot), intent(in):: s
 
-    TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum_$m
-
-    TYPE(shpfileobject), intent(inout):: hshp_outermost
-    ! shapefile outermost_contour_$m
-
-    TYPE(shpfileobject), intent(inout):: hshp_max_speed
-    ! shapefile x_speed_contour_$m
-
     integer, intent(in):: unit_isolated
     ! logical unit for file isolated_nodes_$m.csv
 
@@ -40,11 +27,7 @@ contains
     !------------------------------------------------------------------
 
     if (m == 1 .or. k >= k_begin + max_delta) then
-       ! Write snapshot:
-
        do i = 1, s%number_vis_extr
-          call write_eddy(s%list_vis(i), k, i, hshp_extremum, hshp_outermost, &
-               hshp_max_speed)
           if (s%list_vis(i)%delta_in == huge(0) .and. s%list_vis(i)%delta_out &
                == huge(0)) write(unit_isolated, fmt = *) k, i
        end do