Move call to read_field_indices out of read_snapshot. In program

overlap, we will have to read a set of three shapefiles, extremum,
outermost_contour and max_speed_contour, containing several
snapshots. We do not want to re-read field indices for each snapshot.

Rename u, v to ugos, vgos in NetCDF files.

Documentation_texfol/documentation.tex

@@ -579,33 +579,6 @@ ne ferme pas les contours passant par le bord. J'avais essayé, pour
 inférieure mais cela pose des problèmes de cohérence dans
 l'algorithme.
 
-\subsubsection{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 ? Il
-faudrait tester après l'appel de \verb+get_1_outerm+ si le cadre le
-plus extérieur touche le bord du domaine.
-
 \subsubsection{good\_contour}
 
 Trouve, si elle existe, une ligne de niveau à un niveau donné,
@@ -658,6 +631,31 @@ Il manque des tests pour lesquels direction vaut 1 et 4.
 \subsubsection{set\_all\_outerm}
 \label{sec:set_all_outerm}
 
+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 ? Il
+faudrait tester après l'appel de \verb+get_1_outerm+ si le cadre le
+plus extérieur touche le bord du domaine.
+
 Notons ici :
 \begin{align*}
   & u = \mathtt{urc(1)} \\
@@ -735,14 +733,13 @@ 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
+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.
+weight à u et v, et même u, v à plusieurs dates. Cela remettrait en
+question toute l'organisation des programmes.
 
 \subsubsection{spher\_polyline\_area}
 

Tests/test_read_snapshot.f90

@@ -9,6 +9,7 @@ program test_read_snapshot
 
   use derived_types, only: snapshot
   use init_shapefiles_m, only: init_shapefiles
+  use read_field_indices_m, only: read_field_indices
   use read_snapshot_m, only: read_snapshot
   use write_eddy_m, only: write_eddy
   use write_extr_map_m, only: write_extr_map
@@ -53,6 +54,7 @@ program test_read_snapshot
        hshp = hshp_outermost)
   call shp_open_03("max_speed_contour_old", pszaccess = "rb", &
        hshp = hshp_max_speed)
+  call read_field_indices(hshp_extremum, hshp_outermost, hshp_max_speed)
   call read_snapshot(s, k, hshp_extremum, hshp_outermost, hshp_max_speed, &
        corner * deg_to_rad, step * deg_to_rad, nlon, nlat, periodic, dist_lim)
   CALL shpclose(hshp_extremum)

Tests/test_successive_overlap.f90

@@ -9,6 +9,7 @@ program test_successive_overlap
   use shapelib_03, only: shp_open_03
 
   use derived_types, only: snapshot
+  use read_field_indices_m, only: read_field_indices
   use read_snapshot_m, only: read_snapshot
   use successive_overlap_m, only: successive_overlap
 
@@ -52,6 +53,7 @@ program test_successive_overlap
        hshp = hshp_outermost)
   call shp_open_03("max_speed_contour", pszaccess = "rb", &
        hshp = hshp_max_speed)
+  call read_field_indices(hshp_extremum, hshp_outermost, hshp_max_speed)
   call read_snapshot(flow(1), k, hshp_extremum, hshp_outermost, &
        hshp_max_speed, corner * deg_to_rad, step * deg_to_rad, nlon, nlat, &
        periodic, dist_lim)

extraction_eddies.f90

@@ -23,9 +23,9 @@ program extraction_eddies
   implicit none
 
   type(snapshot) s
-  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
+  TYPE(shpfileobject) hshp_extremum ! shapefile extremum
+  TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour
+  TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour
   integer i, k
 
   real:: min_amp = 0.001

read_snapshot.f90

@@ -7,6 +7,8 @@ contains
   subroutine read_snapshot(s, k, hshp_extremum, hshp_outermost, &
        hshp_max_speed, corner, step, nlon, nlat, periodic, dist_lim)
 
+    ! Note: read_field_indices must have been called before read_snapshot.
+
     ! Libraries:
     use contour_531, only: convert_to_ind 
     use nr_util, only: assert, deg_to_rad
@@ -15,17 +17,16 @@ contains
 
     use derived_types, only: snapshot, eddy
     use read_eddy_m, only: read_eddy
-    use read_field_indices_m, only: read_field_indices
 
     type(snapshot), intent(out):: s ! completely defined
     integer, intent(out):: k ! date index
-    TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum_1
+    TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum
 
     TYPE(shpfileobject), intent(inout):: hshp_outermost
-    ! shapefile outermost_contour_1
+    ! shapefile outermost_contour
 
     TYPE(shpfileobject), intent(inout):: hshp_max_speed
-    ! shapefile max_speed_contour_1
+    ! shapefile max_speed_contour
 
     real, intent(in):: corner(:) ! (2) longitude and latitude of the
     ! corner of the whole grid, in rad
@@ -48,7 +49,6 @@ contains
 
     call shp_get_info_03(hshp_extremum, n_entities = s%number_vis_extr)
     allocate(s%list_vis(s%number_vis_extr))
-    call read_field_indices(hshp_extremum, hshp_outermost, hshp_max_speed)
 
     ! The first shape gives the date index:
     call read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &