From 1a0f1ddf37cc4e50105a652e4900d300e191e31f Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Mon, 6 May 2019 18:31:18 +0200
Subject: [PATCH] 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 |  65 ++++++++++++-------------
 Tests/Input/uv_outermost.nc            | Bin 6428 -> 8048 bytes
 Tests/Input/uv_region_3.nc             | Bin 238272 -> 238340 bytes
 Tests/Input/uv_region_5.nc             | Bin 354180 -> 354252 bytes
 Tests/test_read_snapshot.f90           |   2 +
 Tests/test_successive_overlap.f90      |   2 +
 extraction_eddies.f90                  |   6 +--
 read_snapshot.f90                      |  10 ++--
 8 files changed, 43 insertions(+), 42 deletions(-)

diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index 421054a5..ba58f28a 100644
--- a/Documentation_texfol/documentation.tex
+++ b/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}
 
diff --git a/Tests/Input/uv_outermost.nc b/Tests/Input/uv_outermost.nc
index 8a93c0d74b491ffecfeb7e5264bd09f323ff05c4..beea86cfca7354fc72c63f7fda79dae0470b65cc 100644
GIT binary patch
delta 1743
zcmbPZ^uca|zelETex8DFVx@wDnS!CYm8qqbskwrYfuW_9LSAxFYF=V)s)BBrLa9z^
zdVVp8Ri;x0VwIM~=a-hG7Ukv_m+0jsb8QSzXWU%IxR8^B%|y>!&xnCxvOM>DCJ~Ow
z0zAK%M0h6i@Xlir(V4u4_djFK<TZS=87FM^;4fxmyfJx?Kse_e1_lNekQI{)gxV)B
z7CgfF1SA5oYIA{*A0v~v%H{*YF|74K(SZb`_77y}Gcb&X|7iFRWC)C=zkv*W1_lXS
G)&l^tFh1D;

delta 111
zcmexhH^*p#Kcmvd>FSJ|_cJcw<ls;;(KFREQJS32{hmqqz~l&?UrfSBCWr9OV-kKc
z`3vuVMwQ85_+~SDY@Wbh%*Z%lGK*k1=M)A82F6k#F?okj`{dVxM>ywzM9P4~<{d(Q
Oj7(y8Hgkx?umS)iV<OW4

diff --git a/Tests/Input/uv_region_3.nc b/Tests/Input/uv_region_3.nc
index d63cfb25fa8584f41cb4563b21073d0df41fbfa9..eb5f5793ee06b12de2673effef282d82a69400c9 100644
GIT binary patch
delta 195
zcmX?bl&|F&--KxQY~TDm1>eL<1qCw&LklYdODjV&1tSAPODl!E<f7EP#N1Q`-7<wz
zoznFDVi2oLrwqg@EsHNoP0!5Fi#OKGOXk`bqsh3ropAvt2b+nWxt<XN!)AT%8b-z~
zlh^S4X52A(1@C;u8IwKu{xe2Q_TZnzByBPI34b`_kI6Oy$GGGffPe*L!DI)a?#a^y
k=X1%kfP_FcZFUgSV`O^B(p(_iULec}#7x@@gqZ{S0mdUdsQ>@~

delta 125
zcmZoU#&_T---KvJ)s2fZ88@F~T*S%2p=6?Gs%N6KxtP0#k+ER%2cF-IC6nLq&Swmn
z+{5>uQDkxt|12h{FOxL{!Wp+rZV))eCAWuxfq}6UNNf%e)MDh4I{@OB0g1^0g58t(
ag%>kDIMBRBxP6N-BM>uf-y+N$&<_AQUML^{

diff --git a/Tests/Input/uv_region_5.nc b/Tests/Input/uv_region_5.nc
index dfb2cf2b78e8ab930555a63260e75ea3a8046549..11fd8a21410b91c43120920d6aff7f1b56633a33 100644
GIT binary patch
delta 199
zcmZoUFM8&@=!9tZB;Wiz1>eL<1qCw&LklZIV=EIg1tSAPODl!E<f7EP#N1Q`-7<wz
zoznFDVi2oLrwqg@EsHNoP0!5Fi#OHFOXk`bBhR?Gl5r6y2b+nWxt<XN!)9gf8b-!7
zlNa#(X527&4)1)%hRHU3|CnSHCfo4OVv=K+e1kun@xWvQf#aMf7#J8>Ko)Eg)Mn&7
n1LA=!nIzafdB4yS#xu<+!tE)-j6lo;#LPg<vOPtZRUrfbwEjK%

delta 131
zcmX?eUbN-B=!9rSnT?C&88;thT)@e}p=6?Gs%N6KIh(tNkuhcR3!dMM8Izyz&S$ik
z+`{*dN&3d*7XDdGvO6X#2!u0MOs)_(&e_1gz`$4vBqrAgbx*!5xPY?-BvJ+>HrEL0
dGcvX`uMuuvBg_cIOhC*G#4Ov_2(v1L002seDTV+5

diff --git a/Tests/test_read_snapshot.f90 b/Tests/test_read_snapshot.f90
index 9e3b9eaf..19f33097 100644
--- a/Tests/test_read_snapshot.f90
+++ b/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)
diff --git a/Tests/test_successive_overlap.f90 b/Tests/test_successive_overlap.f90
index a5cb8a80..6d42e716 100644
--- a/Tests/test_successive_overlap.f90
+++ b/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)
diff --git a/extraction_eddies.f90 b/extraction_eddies.f90
index b063fa76..b571c46a 100644
--- a/extraction_eddies.f90
+++ b/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
diff --git a/read_snapshot.f90 b/read_snapshot.f90
index 50b2736c..4438f1bc 100644
--- a/read_snapshot.f90
+++ b/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, &
-- 
GitLab