Redefine `outside_points` in projection space

In main program units `inst_eddies`, `examine_eddy` and
`test_get_1_outerm`, redefine `outside_points` as coordinates in
projection space, avoiding a round trip to longitude and latitude. In
`examine_eddy`, we therefore store `coord_proj` instead of coord.

Inst_eddies/Tests/Input/Region_1/outside_points.csv

-"degrees east" "degrees north"
-"longitude" "latitude"
-8.375 -36.125
-6.125 -35.375
-8.125 -34.625
-10.875 -34.625
-5.875 -33.125
+ x y
+          13           8
+          24           8
+           4          14
+          14           2
+           5           5

Inst_eddies/Tests/examine_eddy.f90

@@ -20,15 +20,13 @@ program examine_eddy
   ! The main difficulty here was to get outside_points. We need in
   ! outside_points all the extrema of opposite sign (including those
   ! not associated to an eddy) plus the extrema of same sign
-  ! associated to an eddy. Note that coord and outside_points may not
-  ! be exactly equal to what they were in inst_eddies because there is
-  ! a multiplication of coord by rad_to_deg when writing the shapefile
-  ! and a multiplication by deg_to_rad when reading it.
+  ! associated to an eddy. Note that e%extr%coord may not be exactly
+  ! equal to what it was in inst_eddies because there is a
+  ! multiplication by rad_to_deg when writing the shapefile and a
+  ! multiplication by deg_to_rad when reading it.
 
   ! Libraries:
-  use contour_531, only: convert_to_ind
-  use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, rad_to_deg, &
-       twopi, new_unit
+  use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, new_unit
   use shapelib, only: shpfileobject, shpclose
   use shapelib_03, only: dbf_read_attribute_03, dbf_get_field_index_03, &
        shp_read_point, shp_open_03
@@ -47,7 +45,7 @@ program examine_eddy
   implicit none
 
   integer:: date = huge(0), eddy_i_target = huge(0)
-  integer ishape_first, d0, ishape, i, l, n1, n2, number_extr, unit
+  integer ishape_first, d0, ishape, i, l, n1, n2, number_extr, unit, nlon
   logical:: cyclone = .true. ! orientation of target eddy
   namelist /main_nml/ cyclone, date, eddy_i_target
   TYPE(shpfileobject) hshp
@@ -71,26 +69,24 @@ program examine_eddy
   ! (nlon, nlat). Wind, in m s-1.
 
   type(snapshot) s
+  real coord(2) ! longitude and latitude, in rad
 
-  real, allocatable:: coord(:, :)
-  ! (2, number_extr) longitude and latitude, in rad
+  integer, allocatable:: coord_proj(:, :)
+  ! (2, number_extr) coordinates in projection space
 
   integer, allocatable:: selection(:)
   ! identifying numbers of a selection of eddies
 
   logical, allocatable:: in_window(:)
 
-  real, allocatable:: outside_points(:, :) ! (2, :) longitude and
-  ! latitude, in rad, of all the significant extrema, except the
-  ! target extremum
+  integer, allocatable:: outside_points(:, :) ! (2, :)
+  ! coordinates in projection space of all the significant extrema,
+  ! except the target extremum
 
   TYPE(shpc_slice_handler) hshpc
   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
-
   !------------------------------------------------------------------------
 
   call config
@@ -122,15 +118,16 @@ program examine_eddy
   end if
 
   number_extr = ishape_last(date) - ishape_first + 1
-  allocate(coord(2, number_extr), in_window(number_extr))
+  allocate(coord_proj(2, number_extr), in_window(number_extr))
 
   do ishape = ishape_first, ishape_last(date)
-     call shp_read_point(coord(:, ishape - ishape_first + 1), hshp, ishape)
+     call shp_read_point(coord, hshp, ishape)
+     coord_proj(:, ishape - ishape_first + 1) = nint((coord * deg_to_rad &
+          - corner) / step + 1.)
   end do
 
-  coord = coord * deg_to_rad
-  e%extr%coord = coord(:, eddy_i_target)
-  e%extr%coord_proj = nint(convert_to_ind(e%extr%coord, corner, step))
+  e%extr%coord = (coord_proj(:, eddy_i_target) - 1) * step + corner
+  e%extr%coord_proj = coord_proj(:, eddy_i_target)
   call dbf_get_field_index_03(hshp, "ssh", ifield)
   call dbf_read_attribute_03(e%extr%ssh, hshp, ifield, &
        ishape = ishape_first + eddy_i_target - 1)
@@ -149,18 +146,18 @@ program examine_eddy
      urc(1) = min(urc(1), size(ssh, 1))
   end if
 
-  corner_window = corner + (llc - 1) * step
-
   ! Done defining geographical window
 
+  if (periodic) nlon = size(ssh, 1) - 2 * max_radius(1)
+
   do i = 1, number_extr
      if (i /= eddy_i_target) then
         ! Shift the longitude to a value close to the longitude of the
         ! target extremum:
-        if (periodic) coord(1, i) = coord(1, i) &
-             + ceiling((corner_window(1) - coord(1, i)) / twopi) * twopi
+        if (periodic) coord_proj(1, i) = coord_proj(1, i) &
+             + ceiling((llc(1) - coord_proj(1, i)) / real(nlon)) * nlon
 
-        in_window(i) = all(abs(nint((coord(:, i) - e%extr%coord) / step)) &
+        in_window(i) = all(abs(coord_proj(:, i) - e%extr%coord_proj) &
              < max_radius)
      else
         in_window(i) = .false.
@@ -174,12 +171,12 @@ program examine_eddy
   selection = pack(selection, s%list(selection)%cyclone .neqv. cyclone)
   n1 = size(selection)
   allocate(outside_points(2, n1 + n2))
-  forall (i = 1:n1) outside_points(:, i) = s%list(selection(i))%extr%coord
+  forall (i = 1:n1) outside_points(:, i) = s%list(selection(i))%extr%coord_proj
   l = n1 + 1
 
   do i = 1, number_extr
      if (in_window(i)) then
-        outside_points(:, l) = coord(:, i)
+        outside_points(:, l) = coord_proj(:, i)
         l = l + 1
      end if
   end do
@@ -190,7 +187,7 @@ program examine_eddy
   call set_contours(e%out_cont, e%speed_cont, e%max_speed, cyclone, e%extr, &
        e%innermost_level, 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)), llc, &
-       convert_to_ind(outside_points, corner, step))
+       outside_points)
   call close_cont_list
   call shpc_create(hshpc, shpc_dir = "SHPC", cyclone = cyclone, slice = 0, &
        grid_lon_lat = .true.)
@@ -199,7 +196,6 @@ program examine_eddy
   CALL shpc_close(hshpc)
   print *, "llc = ", llc
   print *, "urc = ", urc
-  print *, "corner_window (in degrees):", corner_window * rad_to_deg
   print *, "Number of extrema of opposite orientation which must be outside ", &
        "outermost contour:", n1
   print *, "Number of (valid) eddies of same orientation which must be ", &
@@ -207,11 +203,10 @@ program examine_eddy
 
   call new_unit(unit)
   open(unit, file = "outside_points.csv", status = "replace", action = "write")
-  write(unit, fmt = *) '"degrees east" "degrees north"'
-  write(unit, fmt = *) "longitude latitude"
+  write(unit, fmt = *) "x y"
 
   do i = 1, n1 + n2
-     write(unit, fmt = *) outside_points(:, i) * rad_to_deg
+     write(unit, fmt = *) outside_points(:, i)
   end do
 
   close(unit)

Inst_eddies/Tests/test_get_1_outerm.f90

@@ -6,8 +6,8 @@ program test_get_1_outerm
   ! cont_list and cont_list_proj.
 
   ! Libraries:
-  use contour_531, only: polyline, convert_to_ind
-  use jumble, only: csvread, deg_to_rad, rad_to_deg
+  use contour_531, only: polyline
+  use jumble, only: csvread, rad_to_deg
   use shapelib, only: shpfileobject, ftdouble, shpt_polygon, shpclose
   use shapelib_03, only: shp_create_03, dbf_add_field_03, &
        shp_append_object_03, dbf_write_attribute_03
@@ -47,9 +47,9 @@ program test_get_1_outerm
   TYPE(shpc_slice_handler) hshpc
   type(eddy) e
 
-  real, allocatable:: outside_points(:, :) ! (2, :) longitude and
-  ! latitude, in rad, of all the significant extrema, except the
-  ! target extremum
+  integer, allocatable:: outside_points(:, :) ! (2, :)
+  ! coordinates in projection space of all the significant extrema,
+  ! except the target extremum
 
   logical:: grid_lon_lat = .true.
   integer, parameter:: n_max_cont = 31
@@ -89,8 +89,8 @@ program test_get_1_outerm
   inquire(file = "outside_points.csv", exist = exist)
 
   if (exist) then
-     call csvread("outside_points.csv", outside_points, first_r = 3)
-     outside_points = transpose(outside_points) * deg_to_rad
+     call csvread("outside_points.csv", outside_points, first_r = 2)
+     outside_points = transpose(outside_points)
   else
      print *, '"outside_points.csv" not found. Assuming no outside points...'
      allocate(outside_points(2, 0))
@@ -99,8 +99,8 @@ program test_get_1_outerm
   e%extr%coord_proj = ind_targ_extr
   e%extr%coord = corner + (ind_targ_extr - 1) * step
   call get_1_outerm(e%out_cont, cont_list, cont_list_proj, n_cont, cyclone, &
-       real(ind_targ_extr), innermost_level, convert_to_ind(outside_points, &
-       corner, step), ssh, corner = [1., 1.], step = step)
+       real(ind_targ_extr), innermost_level, real(outside_points), ssh, &
+       corner = [1., 1.], step = step)
 
   if (e%out_cont%closed) then
      e%extr%ssh = ssh(ind_targ_extr(1), ind_targ_extr(2))

Inst_eddies/inst_eddies.f90

@@ -7,13 +7,12 @@ program inst_eddies
   use, intrinsic:: ISO_FORTRAN_ENV, only: ERROR_UNIT
 
   ! Libraries:
-  use contour_531, only: convert_to_ind, convert_to_reg_coord
-  use jumble, only: assert, new_unit, argwhere, twopi
+  use jumble, only: assert, new_unit, argwhere
   use numer_rec_95, only: indexx
 
   use config_m, only: config
   use derived_types, only: snapshot, shpc_slice_handler
-  use input_ssh_m, only: input_ssh, max_radius, corner
+  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
@@ -59,12 +58,9 @@ program inst_eddies
   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
+  integer, allocatable:: outside_points(:, :) ! (2, :)
+  ! coordinates in projection space of all the significant extrema,
+  ! except the target extremum
 
   !--------------------------------------------------------------
 
@@ -137,7 +133,7 @@ program inst_eddies
 
   ! Done sorting
 
-  if (.not. periodic) nlon = size(ssh, 1)
+  nlon = size(ssh, 1) - 2 * merge(max_radius(1), 0, periodic)
   nlat = size(ssh, 2)
 
   do l = 1, s%number_extr
@@ -156,25 +152,21 @@ program inst_eddies
         urc(1) = min(urc(1), nlon)
      end if
 
-     corner_window = corner + (llc - 1) * step
-
      ! Done defining geographical window
 
      outside_points &
-          = convert_to_reg_coord(nearby_extr(s%extr_map(llc(1):urc(1), &
-          llc(2):urc(2)), s%list, i), corner, step)
+          = 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:
+     ! Shift the x coordinate of each point to a value close to the
+     ! coordinate of the target extremum:
      if (periodic) outside_points(1, :) = outside_points(1, :) &
-          + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
-          * twopi
+          + ceiling((llc(1) - outside_points(1, :)) / real(nlon)) * nlon
 
      call set_contours(s%list(i)%out_cont, s%list(i)%speed_cont, &
           s%list(i)%max_speed, s%list(i)%cyclone, s%list(i)%extr, &
           s%list(i)%innermost_level, 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)), &
-          llc, convert_to_ind(outside_points, corner, step))
+          llc, outside_points)
   end do
 
   call cpu_time(t1)

Inst_eddies/set_contours.f90

@@ -46,7 +46,7 @@ contains
     ! corresponding to ssh(1, 1). A priori, this is not the corner of
     ! the grid of the NetCDF file.
 
-    real, intent(in):: outside_points_proj(:, :) ! (2, :)
+    integer, intent(in):: outside_points_proj(:, :) ! (2, :)
     ! coordinates in projection space of all the significant extrema,
     ! except the target extremum
 
@@ -94,8 +94,8 @@ contains
          innermost_level, abs(extr%ssh - innermost_level) < min_amp)
 
     call get_1_outerm(out_cont, cont_list, cont_list_proj, n_cont, cyclone, &
-         real(extr%coord_proj), innermost_level_2, outside_points_proj, ssh, &
-         real(corner), step)
+         real(extr%coord_proj), innermost_level_2, real(outside_points_proj), &
+         ssh, real(corner), step)
 
     ! Done with outermost contour, now let us take care of the
     ! max-speed contour.