Extract into new procedure `set_all_extr`

Extract definition of extrema from procedure `set_all_outerm` into new
procedure `set_all_extr`. So we are reverting part of commit
f5f4298d. The drawback, for now, is that `innermost_level` has to
transit in the main program unit from `set_all_extr` to
`set_all_outerm`.

Inst_eddies/CMakeLists.txt

@@ -3,7 +3,7 @@ add_executable(inst_eddies inst_eddies.f90 local_extrema.f90 set_max_speed.f90
   get_1_outerm.f90 max_speed_contour_ssh.f90 good_contour.f90 mean_speed.f90
   inside_4.f90 set_all_outerm.f90 nearby_extr.f90 get_var.f90 config.f90
   input_ssh.f90 shpc_create.f90 write_snapshot.f90 ccw_orient.f90
-  write_eddy.f90)
+  write_eddy.f90 set_all_extr.f90)
 target_link_libraries(inst_eddies PRIVATE Contour_531::contour_531
   Geometry::geometry NetCDF95::netcdf95 Shapelib_03::shapelib_03
   Numer_Rec_95::numer_rec_95 Jumble::jumble NetCDF_Fortran::netcdff)
@@ -24,7 +24,7 @@ target_sources(test_get_1_outerm PRIVATE get_1_outerm.f90 good_contour.f90
 target_sources(test_set_all_outerm PRIVATE set_all_outerm.f90 local_extrema.f90
   get_1_outerm.f90 good_contour.f90 nearby_extr.f90 get_var.f90 config.f90
   input_ssh.f90 shpc_create.f90 write_snapshot.f90 ccw_orient.f90
-  write_eddy.f90)
+  write_eddy.f90 set_all_extr.f90)
 target_sources(test_good_contour PRIVATE good_contour.f90)
 target_sources(test_inside_4 PRIVATE inside_4.f90)
 target_sources(test_mean_speed PRIVATE mean_speed.f90)

Inst_eddies/Tests/test_set_all_outerm.f90

@@ -7,6 +7,7 @@ program test_set_all_outerm
   use derived_types, only: snapshot, shpc_slice_handler, null_ssh_contour, &
        missing_speed
   use input_ssh_m, only: input_ssh
+  use set_all_extr_m, only: set_all_extr
   use set_all_outerm_m, only: set_all_outerm
   use shpc_close_m, only: shpc_close
   use shpc_create_m, only: shpc_create
@@ -33,11 +34,17 @@ program test_set_all_outerm
   logical periodic ! grid is periodic in longitude
   TYPE(shpc_slice_handler) hshp_cyclo, hshp_anti
 
+  real, allocatable:: innermost_level(:) ! (s%number_extr)
+  ! SSH level of the innermost contour around each extremum, in
+  ! m. By construction, innermost_level < extremum for a maximum, >
+  ! extremum for a minimum.
+
   !--------------------------------------------------------------
 
   call config
   call input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v)
-  call set_all_outerm(s, step, periodic, ssh, corner, &
+  call set_all_extr(s, innermost_level, step, periodic, ssh, corner)
+  call set_all_outerm(s, step, periodic, ssh, corner, innermost_level, &
        min_area = pi * (min_radius * 1e3)**2)
 
   do i = 1, s%number_extr

Inst_eddies/inst_eddies.f90

@@ -12,6 +12,7 @@ program inst_eddies
        shpc_slice_handler
   use input_ssh_m, only: input_ssh
   use nearby_extr_m, only: nearby_extr
+  use set_all_extr_m, only: set_all_extr
   use set_all_outerm_m, only: set_all_outerm
   use set_max_speed_m, only: set_max_speed
   use shpc_close_m, only: shpc_close
@@ -56,6 +57,11 @@ program inst_eddies
   ! latitude, in rad, of all the significant extrema, except the
   ! target extremum
 
+  real, allocatable:: innermost_level(:) ! (s%number_extr)
+  ! SSH level of the innermost contour around each extremum, in
+  ! m. By construction, innermost_level < extremum for a maximum, >
+  ! extremum for a minimum.
+
   logical exist
   real t0, t1 ! CPU times, in s
   integer:: date = 0, slice = 0
@@ -100,7 +106,8 @@ program inst_eddies
   write(unit, fmt = *) "CPU time for preamble, before computation:", t1 - t0, &
        "s"
   t0 = t1
-  call set_all_outerm(s, step, periodic, ssh, corner, &
+  call set_all_extr(s, innermost_level, step, periodic, ssh, corner)
+  call set_all_outerm(s, step, periodic, ssh, corner, innermost_level, &
        min_area = pi * (min_radius * 1e3)**2)
 
   ! Done with outermost contours, now let us take care of the

Inst_eddies/set_all_extr.f90

+module set_all_extr_m
+
+  implicit none
+
+contains
+
+  subroutine set_all_extr(s, innermost_level, step, periodic, ssh, corner)
+
+    ! This procedures finds all the extrema of ssh. Not a function
+    ! because "s" is not completely defined.
+
+    use derived_types, only: snapshot
+    use input_ssh_m, only: max_radius
+    use local_extrema_m, only: local_extrema
+
+    type(snapshot), intent(out):: s
+    ! Define only s%extr_map, s%ind_extr, s%number_extr,
+    ! s%list%coord_extr, s%list%ssh_extr, s%list%cyclone
+
+    real, allocatable, intent(out):: innermost_level(:) ! (s%number_extr)
+    ! SSH level of the innermost contour around each extremum, in
+    ! m. By construction, innermost_level < extremum for a maximum, >
+    ! extremum for a minimum.
+
+    real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
+    logical, intent(in):: periodic ! grid is periodic in longitude
+
+    real, intent(in):: ssh(1 - merge(max_radius(1), 0, periodic):, :)
+    ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
+    ! in longitude, else (nlon, nlat). Sea-surface height, in m.
+
+    real, intent(in):: corner(:) ! (2) longitude and latitude of the
+    ! corner of the global grid, in rad
+
+    ! Local:
+
+    integer nlon, nlat
+
+    logical, allocatable:: cyclone(:) ! (s%number_extr)
+    integer i, copy
+
+    !--------------------------------------------------------------
+
+    copy = merge(max_radius(1), 0, periodic)
+    nlon = size(ssh, 1) - 2 * copy
+    nlat = size(ssh, 2)
+    allocate(s%extr_map(1 - copy:nlon + copy, nlat))
+    copy = merge(1, 0, periodic)
+    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\} and s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
+
+    ! Assertion: \forall k, s%ind_extr(1, k) \in \{1, \dots, nlon\} and
+    ! s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
+
+    if (periodic) then
+       ! Extend in longitude:
+       s%extr_map(1 - max_radius(1):- 1, :) &
+            = s%extr_map(nlon + 1 - max_radius(1):nlon - 1, :)
+       s%extr_map(nlon + 2:nlon + max_radius(1), :) &
+            = s%extr_map(2:max_radius(1), :)
+    end if
+
+    s%number_extr = size(s%ind_extr, 2)
+    allocate(s%list(s%number_extr))
+
+    forall (i = 1:s%number_extr)
+       s%list(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(i)%ssh_extr = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
+       s%list(i)%cyclone = cyclone(i)
+    end forall
+
+  end subroutine set_all_extr
+
+end module set_all_extr_m

Inst_eddies/set_all_outerm.f90

@@ -4,11 +4,11 @@ module set_all_outerm_m
 
 contains
 
-  subroutine set_all_outerm(s, step, periodic, ssh, corner, min_area)
+  subroutine set_all_outerm(s, step, periodic, ssh, corner, innermost_level, &
+       min_area)
 
-    ! Extraction of eddies: find extrema and set all outermost
-    ! contours in the snapshot "s". Not a function because "s" is not
-    ! completely defined.
+    ! This procedure sets all outermost contours in the snapshot
+    ! "s".
 
     ! Libraries:
     use jumble, only: argwhere
@@ -20,11 +20,9 @@ contains
     use get_1_outerm_m, only: get_1_outerm
     use input_ssh_m, only: max_radius
     use nearby_extr_m, only: nearby_extr
-    use local_extrema_m, only: local_extrema
 
-    type(snapshot), intent(out):: s
-    ! Specifically: define everything in s except
-    ! s%list%speed_cont, s%list%max_speed and s%list%radius4.
+    type(snapshot), intent(inout):: s
+    ! Define s%list%valid, s%list%out_cont
 
     real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
     logical, intent(in):: periodic ! grid is periodic in longitude
@@ -36,21 +34,19 @@ contains
     real, intent(in):: corner(:) ! (2) longitude and latitude of the
     ! corner of the global grid, in rad
 
+    real, intent(in):: innermost_level(:) ! (s%number_extr)
+    ! SSH level of the innermost contour around each extremum, in
+    ! m. By construction, innermost_level < extremum for a maximum, >
+    ! extremum for a minimum.
+
     real, intent(in):: min_area
     ! minimum area of an outermost contour, in m2
-    
+
     ! Local:
 
     integer nlon, nlat
-
-    real, allocatable:: innermost_level(:) ! (s%number_extr)
-    ! SSH level of the innermost contour around each extremum, in
-    ! m. By construction, innermost_level < extremum for a maximum, >
-    ! extremum for a minimum.
-
     real innermost_level_2
-    logical, allocatable:: cyclone(:) ! (s%number_extr)
-    integer i, l, copy
+    integer i, l
     integer n_cycl ! number of cyclones
 
     integer, allocatable:: sorted_extr(:) ! (s%number_extr)
@@ -75,51 +71,24 @@ contains
 
     !--------------------------------------------------------------
 
-    copy = merge(max_radius(1), 0, periodic)
-    nlon = size(ssh, 1) - 2 * copy
+    if (.not. periodic) nlon = size(ssh, 1)
     nlat = size(ssh, 2)
-    allocate(s%extr_map(1 - copy:nlon + copy, nlat))
-    copy = merge(1, 0, periodic)
-    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\} and s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
-    
-    ! Assertion: \forall k, s%ind_extr(1, k) \in \{1, \dots, nlon\} and
-    ! s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
-
-    if (periodic) then
-       ! Extend in longitude:
-       s%extr_map(1 - max_radius(1):- 1, :) &
-            = s%extr_map(nlon + 1 - max_radius(1):nlon - 1, :)
-       s%extr_map(nlon + 2:nlon + max_radius(1), :) &
-            = s%extr_map(2:max_radius(1), :)
-    end if
-
-    s%number_extr = size(s%ind_extr, 2)
-    allocate(s%list(s%number_extr), sorted_extr(s%number_extr))
 
-    forall (i = 1:s%number_extr)
-       s%list(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(i)%ssh_extr = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
-       s%list(i)%cyclone = cyclone(i)
+    allocate(sorted_extr(s%number_extr))
 
+    forall (i = 1:s%number_extr)
        s%list(i)%valid = .true.
        ! must be intialized to true because it is used in nearby_extr
     end forall
 
     ! Double sort of extrema on cyclonicity and SSH value:
 
-    selection = argwhere(cyclone)
+    selection = argwhere(s%list%cyclone)
     n_cycl = size(selection)
     selection = selection(indexx(s%list(selection)%ssh_extr))
     sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
 
-    selection = argwhere(.not. cyclone)
+    selection = argwhere(.not. s%list%cyclone)
     selection = selection(indexx(s%list(selection)%ssh_extr))
     sorted_extr(n_cycl + 1:) = selection