Make `cont_list` a polyline

Make `cont_list` a polyline and `cont_list_proj` an SSH contour,
instead of vice-versa. The motivation is clarity in `get_1_outerm`: we
only define `cont_list` at the end.

Inst_eddies/Tests/test_get_1_outerm.f90

@@ -49,22 +49,23 @@ program test_get_1_outerm
   logical:: grid_lon_lat = .true.
   integer, parameter:: n_max_cont = 31
 
-  type(ssh_contour) cont_list(n_max_cont)
+  type(polyline) cont_list(n_max_cont)
   ! Contour list, with points given by longitude and latitude, in
-  ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
-  ! the outermost contour is element with subscript n_cont. The
-  ! contours are in monotonic order of ssh. For i <= n_cont,
-  ! cont_list(i)%area has its default initialization value and the
-  ! order of points in cont_list(i)%polyline (clockwise or
+  ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1, the
+  ! outermost contour is element with subscript n_cont. For i <=
+  ! n_cont, the order of points in cont_list(i) (clockwise or
   ! counter-clockwise) is not specified.
 
-  type(polyline) cont_list_proj(n_max_cont)
-  ! Contour list, with points given by projection
-  ! coordinates. Defined only for subscripts 1:n_cont. The contour
-  ! with a given subscript in cont_list_proj corresponds to the
-  ! contour with the same subscript in cont_list. For i <= n_cont,
-  ! the order of points in cont_list_proj(i) (clockwise or
-  ! counter-clockwise) is not specified.
+  type(ssh_contour) cont_list_proj(n_max_cont)
+  ! Contour list, with points given by projection coordinates. Defined
+  ! only for subscripts 1:n_cont. If n_cont >= 1, the outermost
+  ! contour is element with subscript n_cont. The contours are in
+  ! monotonic order of ssh. The contour with a given subscript in
+  ! cont_list_proj corresponds to the contour with the same subscript
+  ! in cont_list. For i <= n_cont, cont_list_proj(i)%area has its
+  ! default initialization value and the order of points in
+  ! cont_list_proj(i) (clockwise or counter-clockwise) is not
+  ! specified.
 
   integer n_cont ! number of good contours found and stored in cont_list
   TYPE(shpfileobject) hshp_cont_list
@@ -117,7 +118,7 @@ program test_get_1_outerm
         call shp_append_object_03(ishape, hshp_cont_list, shpt_polygon, &
              cont_list(i)%points * rad_to_deg)
         call dbf_write_attribute_03(hshp_cont_list, ishape, ifield, &
-             cont_list(i)%ssh)
+             cont_list_proj(i)%ssh)
      end do
 
      CALL shpclose(hshp_cont_list)

Inst_eddies/Tests/test_set_max_speed.f90

@@ -40,22 +40,23 @@ program test_set_max_speed
 
   integer, parameter:: n_max_cont = 31
 
-  type(ssh_contour) cont_list(n_max_cont)
+  type(polyline) cont_list(n_max_cont)
   ! Contour list, with points given by longitude and latitude, in
-  ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
-  ! the outermost contour is element with subscript n_cont. The
-  ! contours are in monotonic order of ssh. For i <= n_cont,
-  ! cont_list(i)%area has its default initialization value and the
-  ! order of points in cont_list(i)%polyline (clockwise or
+  ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1, the
+  ! outermost contour is element with subscript n_cont. For i <=
+  ! n_cont, the order of points in cont_list(i) (clockwise or
   ! counter-clockwise) is not specified.
 
-  type(polyline) cont_list_proj(n_max_cont)
-  ! Contour list, with points given by projection
-  ! coordinates. Defined only for subscripts 1:n_cont. The contour
-  ! with a given subscript in cont_list_proj corresponds to the
-  ! contour with the same subscript in cont_list. For i <= n_cont,
-  ! the order of points in cont_list_proj(i) (clockwise or
-  ! counter-clockwise) is not specified.
+  type(ssh_contour) cont_list_proj(n_max_cont)
+  ! Contour list, with points given by projection coordinates. Defined
+  ! only for subscripts 1:n_cont. If n_cont >= 1, the outermost
+  ! contour is element with subscript n_cont. The contours are in
+  ! monotonic order of ssh. The contour with a given subscript in
+  ! cont_list_proj corresponds to the contour with the same subscript
+  ! in cont_list. For i <= n_cont, cont_list_proj(i)%area has its
+  ! default initialization value and the order of points in
+  ! cont_list_proj(i) (clockwise or counter-clockwise) is not
+  ! specified.
 
   TYPE(shpfileobject) hshp
   type(polygon) p
@@ -102,13 +103,13 @@ program test_set_max_speed
         call shp_read_pol(p, hshp, ishape)
 
         if (p%nparts == 0) then
-           cont_list(i)%polyline = null_polyline()
+           cont_list(i) = null_polyline()
         else
-           cont_list(i)%polyline = p%part(1)
+           cont_list(i) = p%part(1)
            cont_list(i)%points = cont_list(i)%points * deg_to_rad
         end if
 
-        call dbf_read_attribute_03(cont_list(i)%ssh, hshp, ifield = 0, &
+        call dbf_read_attribute_03(cont_list_proj(i)%ssh, hshp, ifield = 0, &
              ishape = ishape)
         cont_list_proj(i)%n_points = cont_list(i)%n_points
         cont_list_proj(i)%closed = cont_list(i)%closed

Inst_eddies/get_1_outerm.f90

@@ -21,12 +21,12 @@ contains
     ! the outermost contour does not have sufficient area. If an
     ! outermost contour is found then its level is farther from the
     ! extremum than innermost_level_2. The outermost contour is
-    ! duplicated in the array cont_list. In cont_list, the contours,
-    ! including the outermost contour, are in monotonic order of ssh:
-    ! ascending if the extremum is a minimum, descending if the
-    ! extremum is a maximum. The area is not computed for contours in
-    ! cont_list. On return, the array
-    ! abs(ediff1d(cont_list(:n_cont)%ssh)) is in descending
+    ! duplicated in the array cont_list. In cont_list_proj, the
+    ! contours, including the outermost contour, are in monotonic
+    ! order of ssh: ascending if the extremum is a minimum, descending
+    ! if the extremum is a maximum. The area is not computed for
+    ! contours in cont_list_proj. On return, the array
+    ! abs(ediff1d(cont_list_proj(:n_cont)%ssh)) is in descending
     ! order. (See notes for a proof.) On return, n_cont == 1 means we
     ! have only found one good contour at innermost_level_2 and it has
     ! sufficient area.
@@ -55,22 +55,22 @@ contains
     type(ssh_contour), intent(out):: out_cont
     ! points given by longitude, latitude, in rad
 
-    type(ssh_contour), intent(out):: cont_list(:) ! (n_max_cont)
+    type(polyline), intent(out):: cont_list(:) ! (n_max_cont)
     ! Contour list, with points given by longitude and latitude, in
     ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
-    ! the outermost contour is element with subscript n_cont. The
-    ! contours are in monotonic order of ssh. For i <= n_cont,
-    ! cont_list(i)%area has its default initialization value and the
-    ! order of points in cont_list(i)%polyline (clockwise or
+    ! the outermost contour is element with subscript n_cont. For i <=
+    ! n_cont, the order of points in cont_list(i) (clockwise or
     ! counter-clockwise) is not specified.
 
-    type(polyline), intent(out):: cont_list_proj(:) ! (n_max_cont)
+    type(ssh_contour), intent(out):: cont_list_proj(:) ! (n_max_cont)
     ! Contour list, with points given by projection
     ! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
     ! 1, the outermost contour is element with subscript n_cont. The
-    ! contour with a given subscript in cont_list_proj corresponds to
-    ! the contour with the same subscript in cont_list. For i <=
-    ! n_cont, the order of points in cont_list_proj(i) (clockwise or
+    ! contours are in monotonic order of ssh. The contour with a given
+    ! subscript in cont_list_proj corresponds to the contour with the
+    ! same subscript in cont_list. For i <= n_cont,
+    ! cont_list_proj(i)%area has its default initialization value and
+    ! the order of points in cont_list_proj(i) (clockwise or
     ! counter-clockwise) is not specified.
 
     integer, intent(out):: n_cont
@@ -111,23 +111,23 @@ contains
     corner_proj = (corner - corner_whole) / step + 1.
     coord_extr_proj = (coord_extr - corner_whole) / step + 1.
     outside_points_proj = convert_to_ind(outside_points, corner_whole, step)
-    n_max_cont = size(cont_list)
+    n_max_cont = size(cont_list_proj)
     n_cont = 0
-    cont_list_proj(1) = good_contour(corner_proj, ssh, innermost_level_2, &
-         coord_extr_proj, outside_points_proj)
+    cont_list_proj(1)%polyline = good_contour(corner_proj, ssh, &
+         innermost_level_2, coord_extr_proj, outside_points_proj)
 
     test_n_points: if (.not. cont_list_proj(1)%closed) then
        ! cont_list_proj(1)%n_points == 0
        out_cont = null_ssh_contour()
     else test_n_points
        n_cont = 1
-       cont_list(1)%ssh = innermost_level_2
+       cont_list_proj(1)%ssh = innermost_level_2
        level_good = innermost_level_2
        mask = ssh /= huge(0.)
        level_try = merge(maxval(ssh, mask), minval(ssh, mask), cyclone)
        call assert(merge(level_try > level_good, level_try < level_good, &
             cyclone), "get_1_outerm level_try")
-       cont_list_proj(2) = good_contour(corner_proj, ssh, level_try, &
+       cont_list_proj(2)%polyline = good_contour(corner_proj, ssh, level_try, &
             coord_extr_proj, outside_points_proj)
 
        if (cont_list_proj(2)%closed) then
@@ -136,24 +136,23 @@ contains
           ! everywhere around the path of contour. So it should
           ! probably never happen.
           n_cont = 2
-          cont_list(2)%ssh = level_try
+          cont_list_proj(2)%ssh = level_try
        else
           ! Assertion: {no good contour at maxval or minval}
           level_bad = level_try
 
           do while (abs(level_bad - level_good) > precision)
              level_try = (level_good + level_bad) / 2.
-             cont_list_proj(n_cont + 1) = good_contour(corner_proj, ssh, &
-                  level_try, coord_extr_proj, outside_points_proj)
+             cont_list_proj(n_cont + 1)%polyline = good_contour(corner_proj, &
+                  ssh, level_try, coord_extr_proj, outside_points_proj)
 
              if (cont_list_proj(n_cont + 1)%closed) then
                 ! cont_list_proj(n_cont + 1)%n_points /= 0
-                cont_list(n_cont + 1)%ssh = level_try
+                cont_list_proj(n_cont + 1)%ssh = level_try
 
                 if (n_cont == n_max_cont - 1) then
                    ! Replace the previous good contour found
                    cont_list_proj(n_cont) = cont_list_proj(n_max_cont)
-                   cont_list(n_cont)%ssh = cont_list(n_max_cont)%ssh
                 else
                    n_cont = n_cont + 1
                 end if
@@ -163,7 +162,7 @@ contains
                 level_bad = level_try
              end if
 
-             ! Assertion: cont_list_proj(n_cont) ==
+             ! Assertion: cont_list_proj(n_cont)%polyline ==
              ! good_contour(. . ., level_good) and
              ! cont_list_proj(n_cont)%closed and no good contour at
              ! level_bad
@@ -171,7 +170,8 @@ contains
           ! Assertion: n_cont <= n_max_cont - 1
        end if
 
-       out_cont%polyline = convert_to_reg_coord(cont_list_proj(n_cont), &
+       out_cont%polyline &
+            = convert_to_reg_coord(cont_list_proj(n_cont)%polyline, &
             corner_whole, step)
        out_cont%area = spher_polyline_area(out_cont%polyline)
 
@@ -180,10 +180,11 @@ contains
           n_cont = 0
        else
           ! n_cont >= 1
-          forall (i = 1:n_cont - 1) cont_list(i)%polyline &
-               = convert_to_reg_coord(cont_list_proj(i), corner_whole, step)
-          cont_list(n_cont)%polyline = out_cont%polyline
-          out_cont%ssh = cont_list(n_cont)%ssh
+          forall (i = 1:n_cont - 1) cont_list(i) &
+               = convert_to_reg_coord(cont_list_proj(i)%polyline, &
+               corner_whole, step)
+          cont_list(n_cont) = out_cont%polyline
+          out_cont%ssh = cont_list_proj(n_cont)%ssh
           call ccw_orient(out_cont)
        end if
     end if test_n_points

Inst_eddies/set_contours.f90

@@ -64,20 +64,21 @@ contains
 
     integer, parameter:: n_max_cont = 31 ! must be >= 3
 
-    type(ssh_contour) cont_list(n_max_cont)
+    type(polyline) cont_list(n_max_cont)
     ! Contour list, with points given by longitude and latitude, in
     ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
-    ! the outermost contour is element with subscript n_cont. The
-    ! contours are in monotonic order of ssh. For i <= n_cont,
-    ! cont_list(i)%area has its default initialization value and the
-    ! order of points in cont_list(i)%polyline (clockwise or
+    ! the outermost contour is element with subscript n_cont. For i <=
+    ! n_cont, the order of points in cont_list(i) (clockwise or
     ! counter-clockwise) is not specified.
 
-    type(polyline) cont_list_proj(n_max_cont)
+    type(ssh_contour) cont_list_proj(n_max_cont)
     ! Contour list, with points given by projection
-    ! coordinates. Defined only for subscripts 1:n_cont. The contour
-    ! with a given subscript in cont_list_proj corresponds to the
-    ! contour with the same subscript in cont_list. For i <= n_cont,
+    ! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
+    ! 1, the outermost contour is element with subscript n_cont. The
+    ! contours are in monotonic order of ssh. The contour with a given
+    ! subscript in cont_list_proj corresponds to the contour with the
+    ! same subscript in cont_list. For i <= n_cont,
+    ! cont_list_proj(i)%area has its default initialization value and
     ! the order of points in cont_list_proj(i) (clockwise or
     ! counter-clockwise) is not specified.
 

Inst_eddies/set_max_speed.f90

@@ -42,20 +42,21 @@ contains
     ! speed_cont is not a null ssh contour then max_speed is the
     ! speed on speed_cont.
 
-    type(ssh_contour), intent(inout):: cont_list(:) ! (n_max_cont)
+    type(polyline), intent(inout):: cont_list(:) ! (n_max_cont)
     ! Contour list, with points given by longitude and latitude, in
-    ! radians. Defined only for subscripts 1:n_cont. On entry, the
-    ! outermost contour is element with subscript n_cont. The contours
-    ! between 1 and i_outer are in monotonic order of
-    ! ssh. cont_list%area has its default initialization value and the
-    ! order of points in the polyline components (clockwise or
+    ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
+    ! the outermost contour is element with subscript n_cont. For i <=
+    ! n_cont, the order of points in cont_list(i) (clockwise or
     ! counter-clockwise) is not specified.
 
-    type(polyline), intent(inout):: cont_list_proj(:) ! (n_max_cont)
+    type(ssh_contour), intent(inout):: cont_list_proj(:) ! (n_max_cont)
     ! Contour list, with points given by projection
-    ! coordinates. Defined only for subscripts 1:n_cont. The contour
-    ! with a given subscript in cont_list_proj corresponds to the
-    ! contour with the same subscript in cont_list. For i <= n_cont,
+    ! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
+    ! 1, the outermost contour is element with subscript n_cont. The
+    ! contours are in monotonic order of ssh. The contour with a given
+    ! subscript in cont_list_proj corresponds to the contour with the
+    ! same subscript in cont_list. For i <= n_cont,
+    ! cont_list_proj(i)%area has its default initialization value and
     ! the order of points in cont_list_proj(i) (clockwise or
     ! counter-clockwise) is not specified.
 
@@ -89,13 +90,13 @@ contains
     corner_proj = (corner - corner_whole) / step + 1.
     extr_coord_proj = (extr_coord - corner_whole) / step + 1.
     i_outer = n_cont
-    if (n_cont >= 2) call complete_ssh(cont_list%ssh, n_cont)
+    if (n_cont >= 2) call complete_ssh(cont_list_proj%ssh, n_cont)
 
     ! Now find the contours associated to the new values of SSH:
     do i = i_outer + 1, n_cont
-       cont_list_proj(i) = good_contour(corner_proj, ssh, cont_list(i)%ssh, &
-            extr_coord_proj)
-       cont_list(i)%polyline = convert_to_reg_coord(cont_list_proj(i), &
+       cont_list_proj(i)%polyline = good_contour(corner_proj, ssh, &
+            cont_list_proj(i)%ssh, extr_coord_proj)
+       cont_list(i) = convert_to_reg_coord(cont_list_proj(i)%polyline, &
             corner_whole, step)
     end do
 
@@ -103,8 +104,8 @@ contains
 
     do i = 1, n_cont
        if (cont_list(i)%closed) then
-          speed(i) = mean_speed(u, v, cont_list(i)%polyline, &
-               cont_list_proj(i)%points, extr_coord, corner, step)
+          speed(i) = mean_speed(u, v, cont_list(i), cont_list_proj(i)%points, &
+               extr_coord, corner, step)
        else
           ! It is possible that good_contour returned a null polyline
           ! if there is a missing value of SSH inside the outermost
@@ -118,7 +119,7 @@ contains
           call shp_append_object_03(ishape, hshp_cont_list, shpt_polygon, &
                cont_list(i)%points * rad_to_deg)
           call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_ssh, &
-               fieldvalue = cont_list(i)%ssh)
+               fieldvalue = cont_list_proj(i)%ssh)
 
           if (ieee_is_nan(speed(i))) then
              call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_speed, &
@@ -145,7 +146,8 @@ contains
           ! Maximum speed on the outermost contour
           speed_cont = null_ssh_contour()
        else
-          speed_cont = cont_list(i)
+          speed_cont%polyline = cont_list(i)
+          speed_cont%ssh = cont_list_proj(i)%ssh
           speed_cont%area = spher_polyline_area(speed_cont%polyline)
           call ccw_orient(speed_cont)
        end if