diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index 9317072e43c7debd7640b5835ec2a0075c793a41..de094c8cd58b348d4e71c8430840f82ffbd894e3 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -14,9 +14,8 @@
\usepackage[np]{numprint}
-\usepackage{hyperref} \hypersetup{pdftitle={Documentation du programme
- en Fortran de détection de tourbillons}, pdfauthor={Lionel Guez},
- hypertexnames=false}
+\usepackage{hyperref} \hypersetup{pdftitle={Détection de tourbillons},
+ pdfauthor={Lionel Guez}, hypertexnames=false}
\usepackage{algorithm}
@@ -40,7 +39,7 @@
\author{Lionel GUEZ}
-\title{Documentation du programme en Fortran de détection de
+\title{Documentation des programmes en Fortran de détection de
tourbillons}
%---------------
@@ -196,8 +195,7 @@ rester sur le même point), allant jusqu'au premier point du maillage
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 \verb+outermost_possible_level+.
+sauf le dernier) varie de façon monotone.
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
@@ -283,8 +281,7 @@ 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 \verb+get_1_outerm+, \verb+local_extrema+, inside et
-\verb+outermost_possible_level+.
+impactées sont \verb+get_1_outerm+, \verb+local_extrema+ et inside.
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
@@ -390,7 +387,7 @@ 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 \verb+set_max_speed+.
+par \verb+nearby_extr+.
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
@@ -589,12 +586,9 @@ 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.
+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}
@@ -1024,8 +1018,9 @@ indicés par le numéro de processus donc deux processus ne peuvent
\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.
+On suppose que les extremums sont sur les n\oe{}uds d'un maillage
+rectangulaire uniforme. On suppose que ce maillage est le même à
+toutes les dates.
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
diff --git a/GNUmakefile b/GNUmakefile
index add7920fb1cf6e6873a032b67b28bcaa93d3fb2b..7cdeaa271b7e21281d1cbe813970f449a81d88c3 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -1,19 +1,14 @@
# This is a makefile for GNU make.
-makefile_dir = .
-include ${general_compiler_options_dir}/settings.mk
-
# 1. Source files and libraries
-VPATH += ${makefile_dir}/Tests
-
-src_test_local_extrema = test_local_extrema.f local_extrema.f
+src_test_local_extrema = test_local_extrema.f local_extrema.f write_extr_map.f
-src_test_get_1_outerm = good_contour.f test_get_1_outerm.f derived_types.f get_1_outerm.f outermost_possible_level.f spher_polyline_area.f
+src_test_get_1_outerm = good_contour.f test_get_1_outerm.f derived_types.f get_1_outerm.f spher_polyline_area.f
src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f good_contour.f max_speed_contour_ssh.f mean_speed.f spher_polyline_area.f inside_4.f
-src_extraction_eddies = extraction_eddies.f write_eddy.f local_extrema.f set_max_speed.f outermost_possible_level.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spher_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.f
+src_extraction_eddies = extraction_eddies.f write_eddy.f local_extrema.f set_max_speed.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spher_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.f
src_test_set_all_outerm = test_set_all_outerm.f derived_types.f set_all_outerm.f local_extrema.f get_1_outerm.f good_contour.f spher_polyline_area.f nearby_extr.f get_var.f
@@ -23,7 +18,7 @@ src_test_spher_polygon_area = spher_polygon_area.f test_spher_polygon_area.f sph
src_test_read_eddy = test_read_eddy.f derived_types.f init_shapefiles.f read_eddy.f write_eddy.f read_field_indices.f
-src_test_read_snapshot = test_read_snapshot.f derived_types.f init_shapefiles.f read_snapshot.f write_eddy.f read_eddy.f read_field_indices.f
+src_test_read_snapshot = test_read_snapshot.f derived_types.f init_shapefiles.f read_snapshot.f write_eddy.f read_eddy.f read_field_indices.f write_extr_map.f
src_test_successive_overlap = test_successive_overlap.f derived_types.f successive_overlap.f read_snapshot.f spher_polygon_area.f spher_polyline_area.f weight.f read_eddy.f read_field_indices.f
@@ -35,6 +30,10 @@ sources := $(sort ${src_test_local_extrema} ${src_test_get_1_outerm} ${src_test_
lib_list = GPC_F contour_531 numer_rec_95 GPC shapelib_03 netcdf95 geometry jumble netcdff fortrangis shp fortranc nr_util
+makefile_dir = .
+include ${general_compiler_options_dir}/settings.mk
+VPATH += ${makefile_dir}/Tests
+
# 2. Objects and executable files
obj_test_local_extrema := $(src_test_local_extrema:.f=.o)
diff --git a/Tests/read_snapshot.f b/Tests/read_snapshot.f
index ca1795a6164c39b7f87b46b15981fc314a39c291..1bc38b8fbb96d5794d628bf6479508608bc7f2ea 100644
--- a/Tests/read_snapshot.f
+++ b/Tests/read_snapshot.f
@@ -5,17 +5,18 @@ module read_snapshot_m
contains
subroutine read_snapshot(s, k, hshp_extremum, hshp_outermost, &
- hshp_max_speed, corner, nlon, nlat)
+ hshp_max_speed, corner, step, nlon, nlat, periodic)
! Libraries:
use contour_531, only: convert_to_ind
use nr_util, only: assert, deg_to_rad
use shapelib, only: shpfileobject
use shapelib_03, only: shp_get_info_03
-
+
use derived_types, only: snapshot, eddy
use read_eddy_m, only: read_eddy
use read_field_indices_m, only: read_field_indices
+ use successive_overlap_m, only: dist_lim
type(snapshot), intent(out):: s ! completely defined
integer, intent(out):: k ! date index
@@ -23,17 +24,19 @@ contains
TYPE(shpfileobject), intent(inout):: hshp_outermost
! shapefile outermost_contour_1
-
+
TYPE(shpfileobject), intent(inout):: hshp_max_speed
! shapefile max_speed_contour_1
real, intent(in):: corner(:) ! (2) longitude and latitude of the
- ! corner of the whole grid, in rad
+ ! corner of the whole grid, in rad
+ real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
integer, intent(in):: nlon, nlat
+ logical, intent(in):: periodic ! grid is periodic in longitude
! Local:
- integer ishape
+ integer ishape, copy
integer k2 ! date index
integer i ! eddy index
type(eddy) e
@@ -48,7 +51,7 @@ contains
call read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
ishape = 0)
s%list_vis(i) = e
-
+
do ishape = 1, s%number_vis_extr - 1
call read_eddy(e, k2, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
ishape)
@@ -58,22 +61,28 @@ contains
allocate(s%ind_extr(2, s%number_vis_extr))
forall (i = 1:s%number_vis_extr) s%ind_extr(:, i) &
- = nint(convert_to_ind(s%list_vis(i)%coord_extr, corner, &
- step = [0.25, 0.25] * deg_to_rad))
+ = nint(convert_to_ind(s%list_vis(i)%coord_extr, corner, step))
- allocate(s%extr_map(nlon, nlat))
- s%extr_map = 0
+ s%number_eddies = s%number_vis_extr
+
+ ! Define s%extr_map from s%ind_extr:
+ copy = merge(dist_lim, 0, periodic)
+ allocate(s%extr_map(1 - copy:nlon + copy, nlat))
+ s%extr_map(1:nlon, :) = 0
+
! Since we are not sure that there is no duplicate tuple in
! s%ind_extr, we must use a do construct and not a forall:
do i = 1, s%number_vis_extr
s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) = i
end do
- call assert(count(s%extr_map /= 0) == s%number_vis_extr, &
+ call assert(count(s%extr_map(1:nlon, :) /= 0) == s%number_vis_extr, &
"read_snapshot: one eddy per grid point")
- s%number_eddies = s%number_vis_extr
+ ! Extend in longitude:
+ s%extr_map(1 - copy:0, :) = s%extr_map(nlon + 1 - copy:nlon, :)
+ s%extr_map(nlon + 1:nlon + copy, :) = s%extr_map(1:copy, :)
end subroutine read_snapshot
diff --git a/Tests/short_tests.json b/Tests/short_tests.json
index 7c2cab5a07523733068b8047014afbe3c6048fba..6a526fb2b213f25c3424759dc0b545bd57649ecd 100644
--- a/Tests/short_tests.json
+++ b/Tests/short_tests.json
@@ -293,22 +293,8 @@
"title" : "Nearby_extr",
"required":
[
- [
- "$tests_old_dir/Extraction_eddies_region_1_noise/extremum.shp",
- "extremum.shp"
- ],
- [
- "$tests_old_dir/Extraction_eddies_region_1_noise/extremum.shx",
- "extremum.shx"
- ],
- [
- "$tests_old_dir/Extraction_eddies_region_1_noise/extremum.dbf",
- "extremum.dbf"
- ],
- [
- "$tests_old_dir/Local_extrema/test_local_extrema.nc",
- "extr_map.nc"
- ]
+ "$tests_old_dir/Extraction_eddies_region_1_noise/extremum.*",
+ "$tests_old_dir/Local_extrema/extr_map.nc"
]
},
{
diff --git a/Tests/test_local_extrema.f b/Tests/test_local_extrema.f
index fe893ffe073894957a1b7726c7ab0ec5571bd495..b63d0b1820c50c886d3f44630441c25ed2340596 100644
--- a/Tests/test_local_extrema.f
+++ b/Tests/test_local_extrema.f
@@ -7,18 +7,17 @@ program test_local_extrema
! Libraries:
use jumble, only: new_unit, get_command_arg_dyn
- use netcdf, only: nf90_nowrite, nf90_clobber, NF90_FLOAT, NF90_INT
+ use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
- nf95_put_var, nf95_create, nf95_def_dim, nf95_def_var, nf95_put_att, &
- nf95_enddef, nf95_gw_var, nf95_get_att
+ nf95_gw_var, nf95_get_att
use local_extrema_m, only: local_extrema
+ use write_extr_map_m, only: write_extr_map
implicit none
integer nlon, nlat
- integer ncid, varid, dimid_lat, dimid_lon
- integer varid_lat, varid_lon, varid_extr_map
+ integer ncid, varid
integer i, unit
character(len = :), allocatable:: filename
real, allocatable:: longitude(:) ! (nlon) in degrees
@@ -95,35 +94,11 @@ program test_local_extrema
close(unit)
print *, 'Created file "test_local_extrema.csv".'
- call nf95_create("test_local_extrema.nc", NF90_CLOBBER, ncid)
-
- call nf95_def_dim(ncid, "lat", nlat, dimid_lat)
- call nf95_def_dim(ncid, "lon", nlon + 2 * copy, dimid_lon)
-
- call nf95_def_var(ncid, "lat", NF90_FLOAT, dimid_lat, varid_lat)
- call nf95_put_att(ncid, varid_lat, "standard_name", "latitude")
- call nf95_put_att(ncid, varid_lat, "units", "degrees_north")
-
- call nf95_def_var(ncid, "lon", NF90_FLOAT, dimid_lon, varid_lon)
- call nf95_put_att(ncid, varid_lon, "standard_name", "longitude")
- call nf95_put_att(ncid, varid_lon, "units", "degrees_east")
-
- call nf95_def_var(ncid, "extr_map", NF90_INT, [dimid_lon, dimid_lat], &
- varid_extr_map)
-
- call nf95_enddef(ncid)
-
- if (periodic) then
- call nf95_put_var(ncid, varid_lon, &
- [longitude(nlon) - 360., longitude, longitude(1) + 360.])
- else
- call nf95_put_var(ncid, varid_lon, longitude)
- end if
-
- call nf95_put_var(ncid, varid_lat, latitude)
- call nf95_put_var(ncid, varid_extr_map, extr_map)
-
- call nf95_close(ncid)
- print *, 'Created file "test_local_extrema.nc".'
+ if (periodic) then
+ call write_extr_map([longitude(nlon) - 360., longitude, longitude(1) &
+ + 360.], latitude, extr_map)
+ else
+ call write_extr_map(longitude, latitude, extr_map)
+ end if
end program test_local_extrema
diff --git a/Tests/test_read_snapshot.f b/Tests/test_read_snapshot.f
index 9ef0336ae6f08ce0d52db9f96e2dfa94b51c0703..917eef36a8c302368ea5ae93863464d8415d0e99 100644
--- a/Tests/test_read_snapshot.f
+++ b/Tests/test_read_snapshot.f
@@ -3,14 +3,16 @@ program test_read_snapshot
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
- use nr_util, only: deg_to_rad
+ use nr_util, only: deg_to_rad, twopi, pi, arth, assert
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 successive_overlap_m, only: dist_lim
use write_eddy_m, only: write_eddy
+ use write_extr_map_m, only: write_extr_map
implicit none
@@ -18,13 +20,15 @@ program test_read_snapshot
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
+ integer i, k, copy
real:: corner(2) = [0.125, - 59.875]
! longitude and latitude of the corner of the whole grid, in degrees
+ real:: step(2) = [0.25, 0.25] ! longitude and latitude steps, in degrees
+ logical periodic ! grid is periodic in longitude
integer:: nlon = 120, nlat = 120
- namelist /main_nml/ corner, nlon, nlat
+ namelist /main_nml/ corner, step, nlon, nlat
!-------------------------------------------------------------------------
@@ -33,13 +37,20 @@ program test_read_snapshot
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
+ ! As we are requiring the grid spacing to be uniform, the value of
+ ! "periodic" must be consistent with the values of step(1) and nlon:
+ periodic = nint(twopi / step(1)) == nlon
+ print *, "periodic = ", periodic
+ if (periodic) call assert(2 * dist_lim * step(1) < pi, &
+ "test_read_snapshot dist_lim")
+
call shp_open_03("extremum_old", pszaccess = "rb", hshp = hshp_extremum)
call shp_open_03("outermost_contour_old", pszaccess = "rb", &
hshp = hshp_outermost)
call shp_open_03("max_speed_contour_old", pszaccess = "rb", &
hshp = hshp_max_speed)
call read_snapshot(s, k, hshp_extremum, hshp_outermost, hshp_max_speed, &
- corner * deg_to_rad, nlon, nlat)
+ corner * deg_to_rad, step * deg_to_rad, nlon, nlat, periodic)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
CALL shpclose(hshp_max_speed)
@@ -66,4 +77,8 @@ program test_read_snapshot
print *, s%ind_extr(:, i)
end do
+ copy = merge(dist_lim, 0, periodic)
+ call write_extr_map(arth(corner(1) - copy * step(1), step(1), nlon + 2 &
+ * copy), arth(corner(2), step(2), nlat), s%extr_map)
+
end program test_read_snapshot
diff --git a/Tests/test_successive_overlap.f b/Tests/test_successive_overlap.f
index 147808c78123e9a1aa9643f92d70cc81ad26f6ae..3c3c47097b29280271200aab4e3ca47d2186408d 100644
--- a/Tests/test_successive_overlap.f
+++ b/Tests/test_successive_overlap.f
@@ -4,13 +4,13 @@ program test_successive_overlap
! Libraries:
use jumble, only: new_unit
- use nr_util, only: deg_to_rad
+ use nr_util, only: deg_to_rad, twopi, assert, pi
use shapelib, only: shpfileobject, shpclose
use shapelib_03, only: shp_open_03
use derived_types, only: snapshot
use read_snapshot_m, only: read_snapshot
- use successive_overlap_m, only: successive_overlap
+ use successive_overlap_m, only: successive_overlap, dist_lim
implicit none
@@ -23,8 +23,10 @@ program test_successive_overlap
real:: corner(2) = [huge(0.), huge(0.)]
! longitude and latitude of the corner of the whole grid, in degrees
+ real:: step(2) = [0.25, 0.25] ! longitude and latitude steps, in degrees
+ logical periodic ! grid is periodic in longitude
integer:: nlon = - 1, nlat = - 1
- namelist /main_nml/ corner, nlon, nlat
+ namelist /main_nml/ corner, step, nlon, nlat
!-------------------------------------------------------------------------
@@ -33,13 +35,21 @@ program test_successive_overlap
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
+ ! As we are requiring the grid spacing to be uniform, the value of
+ ! "periodic" must be consistent with the values of step(1) and nlon:
+ periodic = nint(twopi / step(1)) == nlon
+ print *, "periodic = ", periodic
+ if (periodic) call assert(2 * dist_lim * step(1) < pi, &
+ "test_successive_overlap dist_lim")
+
call shp_open_03("extremum", pszaccess = "rb", hshp = hshp_extremum)
call shp_open_03("outermost_contour", pszaccess = "rb", &
hshp = hshp_outermost)
call shp_open_03("max_speed_contour", pszaccess = "rb", &
hshp = hshp_max_speed)
call read_snapshot(flow(1), k, hshp_extremum, hshp_outermost, &
- hshp_max_speed, corner * deg_to_rad, nlon, nlat)
+ hshp_max_speed, corner * deg_to_rad, step * deg_to_rad, nlon, nlat, &
+ periodic)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
CALL shpclose(hshp_max_speed)
diff --git a/Tests/write_extr_map.f b/Tests/write_extr_map.f
new file mode 100644
index 0000000000000000000000000000000000000000..5ae076b77c6bd9e4ec06567e32dd51eab78110fb
--- /dev/null
+++ b/Tests/write_extr_map.f
@@ -0,0 +1,49 @@
+module write_extr_map_m
+
+ implicit none
+
+contains
+
+ subroutine write_extr_map(longitude, latitude, extr_map)
+
+ use netcdf, only: nf90_clobber, NF90_FLOAT, NF90_INT
+ use netcdf95, only: nf95_close, nf95_put_var, nf95_create, nf95_def_dim, &
+ nf95_def_var, nf95_put_att, nf95_enddef
+
+ real, intent(in):: longitude(:) ! in degrees
+ real, intent(in):: latitude(:) ! in degrees
+ integer, intent(in):: extr_map(:, :) ! map of extrema
+
+ ! Local:
+ integer ncid, varid_lat, varid_lon, varid_extr_map, dimid_lat, dimid_lon
+
+ !--------------------------------------------------------------------------
+
+ call nf95_create("extr_map.nc", NF90_CLOBBER, ncid)
+
+ call nf95_def_dim(ncid, "lat", size(latitude), dimid_lat)
+ call nf95_def_dim(ncid, "lon", size(longitude), dimid_lon)
+
+ call nf95_def_var(ncid, "lat", NF90_FLOAT, dimid_lat, varid_lat)
+ call nf95_put_att(ncid, varid_lat, "standard_name", "latitude")
+ call nf95_put_att(ncid, varid_lat, "units", "degrees_north")
+
+ call nf95_def_var(ncid, "lon", NF90_FLOAT, dimid_lon, varid_lon)
+ call nf95_put_att(ncid, varid_lon, "standard_name", "longitude")
+ call nf95_put_att(ncid, varid_lon, "units", "degrees_east")
+
+ call nf95_def_var(ncid, "extr_map", NF90_INT, [dimid_lon, dimid_lat], &
+ varid_extr_map)
+
+ call nf95_enddef(ncid)
+
+ call nf95_put_var(ncid, varid_lon, longitude)
+ call nf95_put_var(ncid, varid_lat, latitude)
+ call nf95_put_var(ncid, varid_extr_map, extr_map)
+
+ call nf95_close(ncid)
+ print *, 'Created file "extr_map.nc".'
+
+ end subroutine write_extr_map
+
+end module write_extr_map_m
diff --git a/depend.mk b/depend.mk
index efff70c0f45bb29088d11c7d467f52374046c013..48aa00ee2d75bcb259112d47aff18d2eda1f5425 100644
--- a/depend.mk
+++ b/depend.mk
@@ -3,17 +3,17 @@ get_1_outerm.o : spher_polyline_area.o good_contour.o derived_types.o
max_speed_contour_ssh.o : derived_types.o
nearby_extr.o : derived_types.o
read_eddy.o : read_field_indices.o derived_types.o
-read_snapshot.o : read_field_indices.o read_eddy.o derived_types.o
+read_snapshot.o : successive_overlap.o read_field_indices.o read_eddy.o derived_types.o
set_all_outerm.o : local_extrema.o nearby_extr.o get_1_outerm.o derived_types.o
set_max_speed.o : spher_polyline_area.o mean_speed.o max_speed_contour_ssh.o inside_4.o good_contour.o derived_types.o
spher_polygon_area.o : spher_polyline_area.o
successive_overlap.o : weight.o spher_polyline_area.o spher_polygon_area.o derived_types.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
-test_local_extrema.o : local_extrema.o
+test_local_extrema.o : write_extr_map.o local_extrema.o
test_max_speed_contour_ssh.o : max_speed_contour_ssh.o get_var.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 : write_eddy.o read_snapshot.o init_shapefiles.o derived_types.o
+test_read_snapshot.o : write_extr_map.o write_eddy.o successive_overlap.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_spher_polygon_area.o : spher_polygon_area.o
diff --git a/derived_types.f b/derived_types.f
index 8ae356d6afe72f49a1fbdf9ffc2981fdd17967df..9170cb4b293519cf05d503d018da2486bbdfc637 100644
--- a/derived_types.f
+++ b/derived_types.f
@@ -46,10 +46,10 @@ module derived_types
integer number_vis_extr ! number of visible extrema
integer, allocatable:: extr_map(:, :)
- ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is
- ! periodic in longitude, else (nlon, nlat). At a point of
- ! extremum SSH: identification number or this extremum. 0 at
- ! other points.
+ ! (1 - copy:nlon + copy, nlat) if the grid is periodic in
+ ! longitude, else (nlon, nlat). copy can be max_radius(1) or
+ ! disy_lim, depending on the program. At a point of extremum SSH:
+ ! identification number or this extremum. 0 at other points.
integer, allocatable:: ind_extr(:, :) ! (2, number_vis_extr)
! List of coordinates of ssh extrema in the global grid, in index
diff --git a/extraction_eddies.f b/extraction_eddies.f
index 3a16b527cbf171fe49f058fb4cd2c9d05bd0b79d..53d7c896a6424143d591488734a4d6bce8b23d52 100644
--- a/extraction_eddies.f
+++ b/extraction_eddies.f
@@ -112,7 +112,8 @@ program extraction_eddies
! "periodic" must be consistent with the values of step(1) and nlon:
periodic = nint(twopi / step(1)) == nlon
print *, "periodic = ", periodic
- call assert(2 * max_radius(1) * step(1) < pi, "extraction_eddies max_radius")
+ if (periodic) call assert(2 * max_radius(1) * step(1) < pi, &
+ "extraction_eddies max_radius")
copy = merge(max_radius(1), 0, periodic)
allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), &
diff --git a/get_snapshot.f b/get_snapshot.f
index 2b514572b5ade6f6d474d0aab04497498d397933..2e97aca65385065b129a81d61cee0cfa46726bb5 100644
--- a/get_snapshot.f
+++ b/get_snapshot.f
@@ -4,7 +4,8 @@ module get_snapshot_m
contains
- subroutine get_snapshot(s, m, n_proc, k_end, max_delta, nlon, nlat, k, corner)
+ subroutine get_snapshot(s, m, n_proc, k_end, max_delta, nlon, nlat, k, &
+ corner, step, periodic)
use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
@@ -32,6 +33,11 @@ contains
real, intent(in):: corner(:) ! (2) longitude and latitude of the
! corner of the whole grid, in rad
+ real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
+ integer, intent(in):: nlon, nlat
+
+ logical, intent(in):: periodic ! grid is periodic in longitude
+
! Local:
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
@@ -47,13 +53,11 @@ contains
hshp = hshp_outermost)
call shp_open_03("max_speed_contour", pszaccess = "rb", &
hshp = hshp_max_speed)
- call read_snapshot(s, k2, hshp_extremum, hshp_outermost, hshp_max_speed, &
- corner, nlon, nlat)
+ call read_snapshot(s, k2, hshp_extremum, hshp_outermost, &
+ hshp_max_speed, corner, step, nlon, nlat, periodic)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
CALL shpclose(hshp_max_speed)
-
- s%number_eddies = s%number_vis_extr
else
call receive_snapshot(s, source = m + 1, tag = k)
end if
diff --git a/local_extrema.f b/local_extrema.f
index aba4ecd60f9637de71b3325c886854390ff937ab..c0fc4b791b92f0a9e8286465891e20464e74e13d 100644
--- a/local_extrema.f
+++ b/local_extrema.f
@@ -24,7 +24,7 @@ contains
integer, intent(out):: extr_map(:, :) ! (m, n) Map of
! identification numbers of extrema. Identification numbers are
- ! indices in the second dimension of ind_extr. During a call to
+ ! subscripts in the second dimension of ind_extr. During a call to
! local_extrema, we use the array extr_map to temporarily record
! whether the extremum is a maximum or a minimum: extr_map is
! positive for a maximum, negative for a minimum, 0
@@ -34,7 +34,9 @@ contains
! equal to extr_map(2, :).
integer, intent(out), allocatable:: ind_extr(:, :) ! (2, n_extr)
- ! Indices in the two dimensions of extrema of field, relative to my_lbound.
+ ! List of couple of subscripts of extrema of field, in the two
+ ! dimensions of the field, considering that the lower bounds of
+ ! field is given by my_lbound.
real, intent(out), allocatable:: innermost_level(:) ! (n_extr)
! Level of innermost contour, for each extremum. By construction,
@@ -102,10 +104,11 @@ contains
end do
end if
- ! Define ind_extr, relative to (1, 1). We will adjust it later
- ! relative to my_lbound. We do not call pack_indices because we
- ! already know the number of extrema so the following algorithm is
- ! faster (stops at n_extr) and uses less memory.
+ ! Define ind_extr, considering that the lower bounds of field are
+ ! 1 and 1. We will adjust it later for lower bounds my_lbound. We
+ ! do not call pack_indices because we already know the number of
+ ! extrema so the following algorithm is faster (stops at n_extr)
+ ! and uses less memory.
allocate(ind_extr(2, n_extr))
@@ -131,6 +134,9 @@ contains
! {abs(extr_map(ind_extr(1, k), ind_extr(2, k))) == k}
end do
+ ! Assertion: \forall k, ind_extr(1, k) \in \{2, \dots, m - 1\} et
+ ! ind_extr(2, k) \in \{2, \dots, n - 1\}
+
allocate(innermost_level(n_extr), local_min(n_extr))
forall (k = 1:n_extr)
innermost_level(k) = innermost_map(ind_extr(1, k), ind_extr(2, k))
diff --git a/set_all_outerm.f b/set_all_outerm.f
index 01a8d14e2c51c839e5fac95022c776b19f94cc9b..dcff12737c0c8772f6c3b7c14d38b96a2e78cac9 100644
--- a/set_all_outerm.f
+++ b/set_all_outerm.f
@@ -88,6 +88,11 @@ contains
call local_extrema(s%extr_map(1 - copy:nlon + copy, :), s%ind_extr, &
innermost_level, cyclone, ssh(1 - copy:nlon + copy, :), periodic, &
my_lbound = [1 - copy, 1])
+ ! Assertion: \forall k, s%ind_extr(1, k) \in \{2 - copy, \dots,
+ ! nlon + copy - 1\} et s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
+
+ ! Assertion: \forall k, s%ind_extr(1, k) \in \{1, \dots, nlon\} et
+ ! s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
if (periodic) then
! Extend in longitude:
@@ -102,6 +107,9 @@ contains
forall (i = 1:s%number_vis_extr)
s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
+ ! (Even when periodic, this is within the original NetCDF grid,
+ ! that is the grid without duplicated longitudes.)
+
s%list_vis(i)%ssh_extr = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
s%list_vis(i)%cyclone = cyclone(i)
s%list_vis(i)%interpolated = .false.
diff --git a/successive_overlap.f b/successive_overlap.f
index ff710dba50b7f4d4dc1ad0d67f953ea438d3521c..63ba00777d60e37a913e8dfc3286da0258316e6f 100644
--- a/successive_overlap.f
+++ b/successive_overlap.f
@@ -2,6 +2,10 @@ module successive_overlap_m
implicit none
+ integer, parameter:: dist_lim = 12
+ ! We look for an overlapping eddy at dist_lim (in grid points) of
+ ! the first extremum.
+
contains
subroutine successive_overlap(unit_edgelist, nlon, nlat, flow, j, k)
@@ -30,10 +34,6 @@ contains
type(polyline) polyline_1, polyline_2
type(polygon) res_pol
- integer, parameter:: dist_lim = 12
- ! We look for an overlapping eddy at dist_lim (in grid points) of
- ! the first extremum.
-
!-----------------------------------------------------------------------
do i1 = 1, flow(j - 1)%number_vis_extr