diff --git a/Inst_eddies/get_1_outerm.f90 b/Inst_eddies/get_1_outerm.f90 index 7f4c56fcc98a9fa4249f17a428c26fcd9d086e44..e89e3bdaa137d92f6cb59575277b06a30b978b03 100644 --- a/Inst_eddies/get_1_outerm.f90 +++ b/Inst_eddies/get_1_outerm.f90 @@ -88,6 +88,7 @@ contains ! longitude and latitude steps, in rad ! Local: + real level_try, level_good, level_bad ! in m real, parameter:: precision = 5e-4 ! in m logical mask(size(ssh, 1), size(ssh, 2)) @@ -96,6 +97,15 @@ contains real corner_proj(2), coord_extr_proj(2) real outside_points_proj(2, size(outside_points, 2)) + type(polyline) cont_list_proj(size(cont_list)) ! (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 + ! counter-clockwise) is not specified. + !----------------------------------------------------------------- corner_proj = (corner - corner_whole) / step + 1. @@ -103,11 +113,12 @@ contains outside_points_proj = convert_to_ind(outside_points, corner_whole, step) n_max_cont = size(cont_list) n_cont = 0 - cont_list(1)%polyline = convert_to_reg_coord(good_contour(corner_proj, & - ssh, innermost_level_2, coord_extr_proj, outside_points_proj), & + cont_list_proj(1) = good_contour(corner_proj, ssh, innermost_level_2, & + coord_extr_proj, outside_points_proj) + cont_list(1)%polyline = convert_to_reg_coord(cont_list_proj(1), & corner_whole, step) - test_n_points: if (.not. cont_list(1)%closed) then + test_n_points: if (.not. cont_list_proj(1)%closed) then out_cont = null_ssh_contour() else test_n_points n_cont = 1 @@ -117,11 +128,12 @@ contains 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(2)%polyline = convert_to_reg_coord(good_contour(corner_proj, & - ssh, level_try, coord_extr_proj, outside_points_proj), & + cont_list_proj(2) = good_contour(corner_proj, ssh, level_try, & + coord_extr_proj, outside_points_proj) + cont_list(2)%polyline = convert_to_reg_coord(cont_list_proj(2), & corner_whole, step) - if (cont_list(2)%closed) then + if (cont_list_proj(2)%closed) then ! This can only happen if there are degenerate extrema ! everywhere around the path of contour. So it should ! probably never happen. @@ -133,16 +145,18 @@ contains 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(n_cont + 1)%polyline & - = convert_to_reg_coord(good_contour(corner_proj, ssh, & - level_try, coord_extr_proj, outside_points_proj), & + = convert_to_reg_coord(cont_list_proj(n_cont + 1), & corner_whole, step) - if (cont_list(n_cont + 1)%closed) then + if (cont_list_proj(n_cont + 1)%closed) then cont_list(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) = cont_list(n_max_cont) else n_cont = n_cont + 1 @@ -153,9 +167,9 @@ contains level_bad = level_try end if - ! Assertion: cont_list(n_cont)%polyline == + ! Assertion: cont_list_proj(n_cont) == ! good_contour(. . ., level_good) and - ! cont_list(n_cont)%closed and no good contour at + ! cont_list_proj(n_cont)%closed and no good contour at ! level_bad end do ! Assertion: n_cont <= n_max_cont - 1 diff --git a/Inst_eddies/set_max_speed.f90 b/Inst_eddies/set_max_speed.f90 index 29c4047bf47e5c38fe163e0419616f67ceb02713..2e37b5ce2a6bbbe77c384b1d6bf4500cbe81bfbd 100644 --- a/Inst_eddies/set_max_speed.f90 +++ b/Inst_eddies/set_max_speed.f90 @@ -22,7 +22,7 @@ contains use jumble, only: rad_to_deg use shapelib, only: shpt_polygon use shapelib_03, only: shp_append_object_03, dbf_write_attribute_03 - use contour_531, only: convert_to_reg_coord + use contour_531, only: convert_to_reg_coord, polyline use ccw_orient_m, only: ccw_orient use complete_ssh_m, only: complete_ssh @@ -77,6 +77,14 @@ contains integer i, i_outer, ishape real corner_proj(2), extr_coord_proj(2) + type(polyline) cont_list_proj(size(cont_list)) ! (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. + !--------------------------------------------------------------- corner_proj = (corner - corner_whole) / step + 1. @@ -86,8 +94,10 @@ contains ! Now find the contours associated to the new values of SSH: do i = i_outer + 1, n_cont - cont_list(i)%polyline = convert_to_reg_coord(good_contour(corner_proj, & - ssh, cont_list(i)%ssh, extr_coord_proj), corner_whole, step) + 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), & + corner_whole, step) end do allocate(speed(n_cont))