Move the loop on eddies out of `set_all_contours`

Move the loop on eddies out of `set_all_contours`, up into the main
program unit. The motivation is to be able to easily study a
particular eddy, having access to `cont_list` and corresponding speed
values.

Inst_eddies/inst_eddies.f90

@@ -7,9 +7,10 @@ program inst_eddies
   use, intrinsic:: ISO_FORTRAN_ENV
 
   ! Libraries:
-  use jumble, only: assert, new_unit
+  use jumble, only: assert, new_unit, pi, argwhere
+  use numer_rec_95, only: indexx
 
-  use config_m, only: config
+  use config_m, only: config, min_radius
   use derived_types, only: snapshot, shpc_slice_handler
   use input_ssh_m, only: input_ssh
   use set_all_extr_m, only: set_all_extr
@@ -23,7 +24,7 @@ program inst_eddies
 
   type(snapshot) s
   TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
-  integer iostat
+  integer iostat, i, l
   character(len = 200) iomsg
   integer unit
 
@@ -47,6 +48,18 @@ program inst_eddies
   integer:: date = 0, slice = 0
   namelist /dates_nml/ date, slice
 
+  integer, allocatable:: sorted_extr(:) ! (s%number_extr)
+  ! Sorted identifying number of extrema: first, minima sorted in
+  ! order of decreasing SSH, and second, maxima sorted in order of
+  ! increasing SSH.
+
+  real min_area ! minimum area of an outermost contour, in m2
+
+  integer, allocatable:: selection(:)
+  ! identifying numbers of a selection of extrema
+
+  integer n_cycl ! number of cyclones
+
   !--------------------------------------------------------------
 
   call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
@@ -96,7 +109,33 @@ program inst_eddies
        "s"
   t0 = t1
   call set_all_extr(s, step, periodic, ssh, corner)
-  call set_all_contours(s, step, periodic, ssh, u, v, corner)
+
+  allocate(sorted_extr(s%number_extr))
+  min_area = pi * (min_radius * 1e3)**2
+
+  forall (i = 1:s%number_extr)
+     s%list(i)%out_cont%closed = .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(s%list%cyclone)
+  n_cycl = size(selection)
+  selection = selection(indexx(s%list(selection)%extr%ssh))
+  sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
+
+  selection = argwhere(.not. s%list%cyclone)
+  selection = selection(indexx(s%list(selection)%extr%ssh))
+  sorted_extr(n_cycl + 1:) = selection
+
+  ! Done sorting
+
+  loop_extr: do l = 1, s%number_extr
+     i = sorted_extr(l)
+     call set_all_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
+  end do loop_extr
+
   call cpu_time(t1)
   write(unit, fmt = *) "CPU time for computation, before output:", t1 - t0, "s"
   t0 = t1

Inst_eddies/set_all_contours.f90

@@ -4,16 +4,15 @@ module set_all_contours_m
 
 contains
 
-  subroutine set_all_contours(s, step, periodic, ssh, u, v, corner)
+  subroutine set_all_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
 
     ! This procedure sets all contours in the snapshot "s".
 
     ! Libraries:
     use contour_531, only: convert_to_ind
-    use jumble, only: argwhere, twopi, pi
-    use numer_rec_95, only: indexx
+    use jumble, only: twopi
 
-    use config_m, only: min_amp, min_radius
+    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
@@ -37,21 +36,13 @@ contains
     real, intent(in):: corner(:) ! (2) longitude and latitude of the
     ! corner of the global grid, in rad
 
+    real, intent(in):: min_area ! minimum area of an outermost contour, in m2
+    integer, intent(in):: i
+
     ! Local:
 
-    real min_area ! minimum area of an outermost contour, in m2
     integer nlon, nlat
     real innermost_level_2
-    integer i, l
-    integer n_cycl ! number of cyclones
-
-    integer, allocatable:: sorted_extr(:) ! (s%number_extr)
-    ! Sorted identifying number of extrema: first, minima sorted in
-    ! order of decreasing SSH, and second, maxima sorted in order of
-    ! increasing SSH.
-
-    integer, allocatable:: selection(:)
-    ! identifying numbers of a selection of extrema
 
     ! Window around each extremum:
 
@@ -83,88 +74,64 @@ contains
 
     if (.not. periodic) nlon = size(ssh, 1)
     nlat = size(ssh, 2)
-    allocate(sorted_extr(s%number_extr))
-    min_area = pi * (min_radius * 1e3)**2
 
-    forall (i = 1:s%number_extr)
-       s%list(i)%out_cont%closed = .true.
-       ! must be intialized to true because it is used in nearby_extr
-    end forall
+    associate (e => s%list(i))
+      ! Define the geographical window around each eddy extremum:
 
-    ! Double sort of extrema on cyclonicity and SSH value:
+      llc = e%extr%coord_proj - max_radius
+      urc = e%extr%coord_proj + max_radius
 
-    selection = argwhere(s%list%cyclone)
-    n_cycl = size(selection)
-    selection = selection(indexx(s%list(selection)%extr%ssh))
-    sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
+      llc(2) = max(llc(2), 1)
+      urc(2) = min(urc(2), nlat)
 
-    selection = argwhere(.not. s%list%cyclone)
-    selection = selection(indexx(s%list(selection)%extr%ssh))
-    sorted_extr(n_cycl + 1:) = selection
+      if (.not. periodic) then
+         llc(1) = max(llc(1), 1)
+         urc(1) = min(urc(1), nlon)
+      end if
 
-    ! Done sorting
+      corner_window = corner + (llc - 1) * step
 
-    loop_extr: do l = 1, s%number_extr
-       i = sorted_extr(l)
+      ! 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)
 
-       associate (e => s%list(i))
-         ! Define the geographical window around each eddy extremum:
+      outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
+           llc(2):urc(2)), s%list, i)
 
-         llc = e%extr%coord_proj - max_radius
-         urc = e%extr%coord_proj + max_radius
+      ! 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
 
-         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
+      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)
 
-         ! 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)
+      ! Done with outermost contour, now let us take care of the
+      ! max-speed contour.
+      if (e%out_cont%closed) then
+         ! Restrict the field to the outermost contour:
 
-         outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
-              llc(2):urc(2)), s%list, i)
+         llc = floor(convert_to_ind(minval(e%out_cont%points, dim = 2), &
+              corner, step))
 
-         ! 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
+         urc = ceiling(convert_to_ind(maxval(e%out_cont%points, dim = 2), &
+              corner, step))
 
-         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)
-
-         ! Done with outermost contour, now let us take care of the
-         ! max-speed contour.
-         if (e%out_cont%closed) then
-            ! Restrict the field to the outermost contour:
-
-            llc = floor(convert_to_ind(minval(e%out_cont%points, dim = 2), &
-                 corner, step))
-
-            urc = ceiling(convert_to_ind(maxval(e%out_cont%points, dim = 2), &
-                 corner, step))
-
-            ! Should have no effect except because of roundup error:
-            urc(2) = min(urc(2), nlat)
-            if (.not. periodic) urc(1) = min(urc(1), nlon)
+         ! Should have no effect except because of roundup error:
+         urc(2) = min(urc(2), nlat)
+         if (.not. periodic) urc(1) = min(urc(1), nlon)
 
-            corner_window = corner + (llc - 1) * step
+         corner_window = corner + (llc - 1) * step
 
-            call set_max_speed(e, cont_list, n_cont, &
-                 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, step)
-         end if
-       end associate
-    end do loop_extr
+         call set_max_speed(e, cont_list, n_cont, &
+              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, step)
+      end if
+    end associate
 
   end subroutine set_all_contours