From 24c9eac08efa66cc56d7d751422dbcd8b2efcc8a Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ipsl.fr>
Date: Wed, 11 Sep 2024 18:50:38 +0200
Subject: [PATCH] Create function `invert_proj`

Preparing to generalize to a non-uniform longitude-latitude grid.
---
 Inst_eddies/CMakeLists.txt              | 10 +++---
 Inst_eddies/Tests/examine_eddy.f90      |  3 +-
 Inst_eddies/Tests/test_get_1_outerm.f90 |  5 +--
 Inst_eddies/get_1_outerm.f90            |  9 +++--
 Inst_eddies/invert_proj.f90             | 44 +++++++++++++++++++++++++
 Inst_eddies/set_all_extr.f90            |  5 +--
 Inst_eddies/set_max_speed.f90           |  5 ++-
 7 files changed, 64 insertions(+), 17 deletions(-)
 create mode 100644 Inst_eddies/invert_proj.f90

diff --git a/Inst_eddies/CMakeLists.txt b/Inst_eddies/CMakeLists.txt
index fff8c050..a654d80e 100644
--- a/Inst_eddies/CMakeLists.txt
+++ b/Inst_eddies/CMakeLists.txt
@@ -2,7 +2,8 @@
 add_executable(inst_eddies inst_eddies.f90 local_extrema.f90 set_max_speed.f90
   get_1_outerm.f90 good_contour.f90 mean_speed.f90 set_contours.f90
   nearby_extr.f90 config.f90 input_ssh.f90 shpc_create.f90 write_snapshot.f90
-  ccw_orient.f90 write_eddy.f90 set_all_extr.f90 complete_ssh.f90)
+  ccw_orient.f90 write_eddy.f90 set_all_extr.f90 complete_ssh.f90
+  invert_proj.f90)
 target_link_libraries(inst_eddies PRIVATE Contour_531::contour_531
   Geometry::geometry NetCDF95::netcdf95 Shapelib_03::shapelib_03
   Numer_Rec_95::numer_rec_95 Jumble::jumble NetCDF_Fortran::netcdff)
@@ -18,17 +19,18 @@ configure_file(inst_eddies_Aviso.py . COPYONLY FILE_PERMISSIONS
 # Tests:
 add_subdirectory(Tests)
 target_sources(test_get_1_outerm PRIVATE get_1_outerm.f90 good_contour.f90
-  config.f90 input_ssh.f90 shpc_create.f90 write_eddy.f90 ccw_orient.f90)
+  config.f90 input_ssh.f90 shpc_create.f90 write_eddy.f90 ccw_orient.f90
+  invert_proj.f90)
 target_sources(test_good_contour PRIVATE good_contour.f90)
 target_sources(test_mean_speed PRIVATE mean_speed.f90 input_ssh.f90 config.f90)
 target_sources(test_nearby_extr PRIVATE nearby_extr.f90)
 target_sources(test_local_extrema PRIVATE local_extrema.f90)
 target_sources(test_set_max_speed PRIVATE set_max_speed.f90 good_contour.f90
   mean_speed.f90 config.f90 input_ssh.f90 shpc_create.f90 ccw_orient.f90
-  write_eddy.f90 complete_ssh.f90)
+  write_eddy.f90 complete_ssh.f90 invert_proj.f90)
 target_sources(test_write_null PRIVATE shpc_create.f90)
 target_sources(test_complete_ssh PRIVATE complete_ssh.f90)
 target_sources(examine_eddy PRIVATE config.f90 input_ssh.f90 shpc_create.f90
   set_all_extr.f90 local_extrema.f90 write_eddy.f90 set_contours.f90
   get_1_outerm.f90 ccw_orient.f90 good_contour.f90 set_max_speed.f90
-  complete_ssh.f90 mean_speed.f90)
+  complete_ssh.f90 mean_speed.f90 invert_proj.f90)
diff --git a/Inst_eddies/Tests/examine_eddy.f90 b/Inst_eddies/Tests/examine_eddy.f90
index 526a2909..03bca16e 100644
--- a/Inst_eddies/Tests/examine_eddy.f90
+++ b/Inst_eddies/Tests/examine_eddy.f90
@@ -36,6 +36,7 @@ program examine_eddy
   use cont_list_m, only: create_cont_list, close_cont_list
   use derived_types, only: snapshot, shpc_slice_handler, eddy
   use input_ssh_m, only: input_ssh, corner_whole, step, nlon, uniform_lon_lat
+  use invert_proj_m, only: invert_proj
   use read_eddy_m, only: read_eddy
   use shpc_create_m, only: shpc_create
   use set_all_extr_m, only: set_all_extr
@@ -155,7 +156,7 @@ program examine_eddy
      end if
   end do
 
-  e%extr%coord = (coord_proj(:, eddy_i_target) - 1) * step + corner_whole
+  e%extr%coord = invert_proj(coord_proj(:, eddy_i_target))
   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, &
diff --git a/Inst_eddies/Tests/test_get_1_outerm.f90 b/Inst_eddies/Tests/test_get_1_outerm.f90
index 64151470..26279587 100644
--- a/Inst_eddies/Tests/test_get_1_outerm.f90
+++ b/Inst_eddies/Tests/test_get_1_outerm.f90
@@ -16,7 +16,8 @@ program test_get_1_outerm
   use derived_types, only: eddy, shpc_slice_handler, null_ssh_contour, &
        missing_speed, ssh_contour
   use get_1_outerm_m, only: get_1_outerm
-  use input_ssh_m, only: input_ssh, corner_whole, step, uniform_lon_lat
+  use input_ssh_m, only: input_ssh, uniform_lon_lat
+  use invert_proj_m, only: invert_proj
   use shpc_close_m, only: shpc_close
   use shpc_create_m, only: shpc_create
   use write_eddy_m, only: write_eddy
@@ -95,7 +96,7 @@ program test_get_1_outerm
   end if
 
   e%extr%coord_proj = ind_targ_extr
-  e%extr%coord = corner_whole + (ind_targ_extr - 1) * step
+  e%extr%coord = invert_proj(ind_targ_extr)
   call get_1_outerm(e%out_cont, cont_list, cont_list_proj, cont_list_ssh, &
        n_cont, cyclone, real(ind_targ_extr), innermost_level, &
        real(outside_points), ssh, corner = [1., 1.])
diff --git a/Inst_eddies/get_1_outerm.f90 b/Inst_eddies/get_1_outerm.f90
index 699c1fe6..3520a08c 100644
--- a/Inst_eddies/get_1_outerm.f90
+++ b/Inst_eddies/get_1_outerm.f90
@@ -42,14 +42,14 @@ contains
     ! automatic shifting by a mulitple of 360 degrees.
 
     ! Libraries:
-    use contour_531, only: polyline, convert_to_reg_coord
+    use contour_531, only: polyline
     use jumble, only: assert
 
     use ccw_orient_m, only: ccw_orient
     use config_m, only: min_area
     use derived_types, only: ssh_contour, null_ssh_contour
     use good_contour_m, only: good_contour
-    use input_ssh_m, only: corner_whole, step
+    use invert_proj_m, only: invert_proj
     use spher_polyline_area_m, only: spher_polyline_area
 
     type(ssh_contour), intent(out):: out_cont
@@ -168,8 +168,7 @@ contains
           ! Assertion: n_cont <= n_max_cont - 1
        end if
 
-       out_cont%polyline &
-            = convert_to_reg_coord(cont_list_proj(n_cont), corner_whole, step)
+       out_cont%polyline = invert_proj(cont_list_proj(n_cont))
        out_cont%area = spher_polyline_area(out_cont%polyline)
 
        if (abs(out_cont%area) < min_area) then
@@ -178,7 +177,7 @@ contains
        else
           ! n_cont >= 1
           forall (i = 1:n_cont - 1) cont_list(i) &
-               = convert_to_reg_coord(cont_list_proj(i), corner_whole, step)
+               = invert_proj(cont_list_proj(i))
           cont_list(n_cont) = out_cont%polyline
           out_cont%ssh = cont_list_ssh(n_cont)
 
diff --git a/Inst_eddies/invert_proj.f90 b/Inst_eddies/invert_proj.f90
new file mode 100644
index 00000000..aa0bcdcc
--- /dev/null
+++ b/Inst_eddies/invert_proj.f90
@@ -0,0 +1,44 @@
+module invert_proj_m
+
+  use input_ssh_m, only: corner_whole, step
+
+  implicit none
+
+  interface invert_proj
+     ! Invert projection: convert from projection coordinates x, y to
+     ! longitude and latitude.
+
+     module procedure invert_proj_int, invert_proj_pol
+  end interface invert_proj
+
+  private
+  public invert_proj
+
+contains
+
+  pure function invert_proj_int(ind)
+
+    real invert_proj_int(2)
+    integer, intent(in):: ind(:) ! (2)
+
+    !-------------------------------------------------------
+
+    invert_proj_int = corner_whole + (ind - 1.) * step
+
+  end function invert_proj_int
+
+  !**************************************************************************
+
+  pure type(polyline) function invert_proj_pol(ind)
+
+    use contour_531, only: polyline, convert_to_reg_coord
+
+    type(polyline), intent(in):: ind
+
+    !-------------------------------------------------------
+
+    invert_proj_pol = convert_to_reg_coord(ind, corner_whole, step)
+
+  end function invert_proj_pol
+
+end module invert_proj_m
diff --git a/Inst_eddies/set_all_extr.f90 b/Inst_eddies/set_all_extr.f90
index 1689d7e7..0f951d44 100644
--- a/Inst_eddies/set_all_extr.f90
+++ b/Inst_eddies/set_all_extr.f90
@@ -10,7 +10,8 @@ contains
 
     use config_m, only: max_radius
     use derived_types, only: eddy
-    use input_ssh_m, only: corner_whole, step, nlon, nlat
+    use input_ssh_m, only: nlon, nlat
+    use invert_proj_m, only: invert_proj
     use local_extrema_m, only: local_extrema
 
     integer, allocatable, intent(out):: extr_map(:, :)
@@ -74,7 +75,7 @@ contains
     list%innermost_level = innermost_level
 
     forall (i = 1:number_extr)
-       list(i)%extr%coord = corner_whole + (ind_extr(:, i) - 1) * step
+       list(i)%extr%coord = invert_proj(ind_extr(:, i))
        ! (Even when periodic, this is within the original NetCDF grid,
        ! that is the grid without duplicated longitudes.)
 
diff --git a/Inst_eddies/set_max_speed.f90 b/Inst_eddies/set_max_speed.f90
index d7048759..9ae3593f 100644
--- a/Inst_eddies/set_max_speed.f90
+++ b/Inst_eddies/set_max_speed.f90
@@ -29,7 +29,7 @@ contains
     use complete_ssh_m, only: complete_ssh
     use derived_types, only: eddy, null_ssh_contour, ssh_contour, missing_speed
     use good_contour_m, only: good_contour
-    use input_ssh_m, only: corner_whole, step
+    use invert_proj_m, only: invert_proj
     use mean_speed_m, only: mean_speed
     use spher_polyline_area_m, only: spher_polyline_area
     use cont_list_m, only: write_cont_list, hshp_cont_list, &
@@ -99,8 +99,7 @@ contains
     do i = i_outer + 1, n_cont
        cont_list_proj(i) = good_contour(corner, ssh, cont_list_ssh(i), &
             extr_coord_proj)
-       cont_list(i) = convert_to_reg_coord(cont_list_proj(i), corner_whole, &
-            step)
+       cont_list(i) = invert_proj(cont_list_proj(i))
     end do
 
     allocate(speed(n_cont))
-- 
GitLab