Move some computation to main program unit

Move the computation of the search window for the outermost contour
and `outside_points` from procedure `set_contours` to main program
unit. The motivation is to have a clearer interface of procedure
`set_contours`, without `s` and `i`.

Inst_eddies/inst_eddies.f90

@@ -7,12 +7,13 @@ program inst_eddies
   use, intrinsic:: ISO_FORTRAN_ENV
 
   ! Libraries:
-  use jumble, only: assert, new_unit, pi, argwhere
+  use jumble, only: assert, new_unit, pi, argwhere, twopi
   use numer_rec_95, only: indexx
 
   use config_m, only: config, min_radius
   use derived_types, only: snapshot, shpc_slice_handler
-  use input_ssh_m, only: input_ssh
+  use input_ssh_m, only: input_ssh, max_radius
+  use nearby_extr_m, only: nearby_extr
   use set_all_extr_m, only: set_all_extr
   use set_contours_m, only: set_contours
   use shpc_close_m, only: shpc_close
@@ -24,7 +25,7 @@ program inst_eddies
 
   type(snapshot) s
   TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
-  integer iostat, i, l
+  integer iostat, i, l, nlon, nlat
   character(len = 200) iomsg
   integer unit
 
@@ -59,6 +60,16 @@ program inst_eddies
   ! identifying numbers of a selection of extrema
 
   integer n_cycl ! number of cyclones
+  integer llc(2) ! indices in global grid of lower left corner
+  integer urc(2) ! indices in global grid of upper right corner
+
+  real corner_window(2) ! longitude and latitude of the window around each
+  ! extremum, in rad
+
+  real, allocatable:: outside_points(:, :) ! (2, :) longitude and
+  ! latitude, in rad, of all the significant extrema, except the
+  ! target extremum
+
 
   !--------------------------------------------------------------
 
@@ -133,7 +144,38 @@ program inst_eddies
 
   loop_extr: do l = 1, s%number_extr
      i = sorted_extr(l)
-     call set_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
+     if (.not. periodic) nlon = size(ssh, 1)
+     nlat = size(ssh, 2)
+
+     associate (e => s%list(i))
+       ! Define the geographical window around each eddy extremum:
+
+       llc = e%extr%coord_proj - max_radius
+       urc = e%extr%coord_proj + max_radius
+
+       llc(2) = max(llc(2), 1)
+       urc(2) = min(urc(2), nlat)
+
+       if (.not. periodic) then
+          llc(1) = max(llc(1), 1)
+          urc(1) = min(urc(1), nlon)
+       end if
+
+       corner_window = corner + (llc - 1) * step
+     end associate
+
+     outside_points = nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
+          s%list, i)
+
+     ! Shift the longitude of each point to a value close to the
+     ! longitude of the target extremum:
+     if (periodic) outside_points(1, :) = outside_points(1, :) &
+          + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
+          * twopi
+
+     call set_contours(s, step, ssh(llc(1):urc(1), llc(2):urc(2)), &
+          u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), &
+          corner_window, min_area, i, outside_points)
   end do loop_extr
 
   call cpu_time(t1)

Inst_eddies/set_contours.f90

@@ -4,19 +4,17 @@ module set_contours_m
 
 contains
 
-  subroutine set_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
+  subroutine set_contours(s, step, ssh, u, v, corner, min_area, i, &
+       outside_points)
 
     ! This procedure sets all contours in the snapshot "s".
 
     ! Libraries:
     use contour_531, only: convert_to_ind
-    use jumble, only: twopi
 
     use config_m, only: min_amp
     use derived_types, only: snapshot, null_ssh_contour, ssh_contour
     use get_1_outerm_m, only: get_1_outerm
-    use input_ssh_m, only: max_radius
-    use nearby_extr_m, only: nearby_extr
     use set_max_speed_m, only: set_max_speed
 
     type(snapshot), intent(inout):: s
@@ -25,9 +23,8 @@ contains
     ! s%list%innermost_level, s%extr_map should be defined on entry.
 
     real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
-    logical, intent(in):: periodic ! grid is periodic in longitude
 
-    real, intent(in), dimension(merge(1 - max_radius(1), 1, periodic):, :):: &
+    real, intent(in), dimension(:, :):: &
          ssh, u, v
     ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is
     ! periodic in longitude, else (nlon, nlat). Sea-surface height, in
@@ -39,6 +36,10 @@ contains
     real, intent(in):: min_area ! minimum area of an outermost contour, in m2
     integer, intent(in):: i
 
+    real, intent(in):: outside_points(:, :) ! (2, :) longitude and
+    ! latitude, in rad, of all the significant extrema, except the
+    ! target extremum
+
     ! Local:
 
     integer nlon, nlat
@@ -52,10 +53,6 @@ contains
     real corner_window(2) ! longitude and latitude of the window around each
     ! extremum, in rad
 
-    real, allocatable:: outside_points(:, :) ! (2, :) longitude and
-    ! latitude, in rad, of all the significant extrema, except the
-    ! target extremum
-
     integer, parameter:: n_max_cont = 31 ! must be >= 3
 
     type(ssh_contour) cont_list(n_max_cont)
@@ -72,42 +69,18 @@ contains
 
     !--------------------------------------------------------------
 
-    if (.not. periodic) nlon = size(ssh, 1)
+    nlon = size(ssh, 1)
     nlat = size(ssh, 2)
 
     associate (e => s%list(i))
-      ! Define the geographical window around each eddy extremum:
-
-      llc = e%extr%coord_proj - max_radius
-      urc = e%extr%coord_proj + max_radius
-
-      llc(2) = max(llc(2), 1)
-      urc(2) = min(urc(2), nlat)
-
-      if (.not. periodic) then
-         llc(1) = max(llc(1), 1)
-         urc(1) = min(urc(1), nlon)
-      end if
-
-      corner_window = corner + (llc - 1) * step
-
       ! Only look for good contours with amplitudes >= min_amp:
       innermost_level_2 = merge(e%extr%ssh + merge(min_amp, - min_amp, &
            e%cyclone), e%innermost_level, &
            abs(e%extr%ssh - e%innermost_level) < min_amp)
 
-      outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
-           llc(2):urc(2)), s%list, i)
-
-      ! Shift the longitude of each point to a value close to the
-      ! longitude of the target extremum:
-      if (periodic) outside_points(1, :) = outside_points(1, :) &
-           + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
-           * twopi
-
       call get_1_outerm(e%out_cont, cont_list, n_cont, e%cyclone, &
            e%extr%coord, innermost_level_2, outside_points, &
-           ssh(llc(1):urc(1), llc(2):urc(2)), corner_window, step, min_area)
+           ssh, corner, step, min_area)
 
       ! Done with outermost contour, now let us take care of the
       ! max-speed contour.
@@ -122,7 +95,7 @@ contains
 
          ! Should have no effect except because of roundup error:
          urc(2) = min(urc(2), nlat)
-         if (.not. periodic) urc(1) = min(urc(1), nlon)
+         urc(1) = min(urc(1), nlon)
 
          corner_window = corner + (llc - 1) * step