Create program `examine_eddy`

Common/CMakeLists.txt

@@ -12,6 +12,8 @@ target_sources(inst_eddies PRIVATE spher_polyline_area.f90
   get_slice_dir.f90 read_field_indices.f90)
 target_sources(test_write_null PRIVATE derived_types.f90 shpc_open.f90
   shpc_close.f90 read_field_indices.f90 get_slice_dir.f90)
+target_sources(examine_eddy PRIVATE derived_types.f90 get_slice_dir.f90
+  spher_polyline_area.f90 shpc_close.f90)
 
 if(MPI_Fortran_HAVE_F08_MODULE)
   target_sources(test_overlap PRIVATE derived_types.f90 spher_polyline_area.f90

Inst_eddies/CMakeLists.txt

@@ -30,3 +30,7 @@ target_sources(test_set_max_speed PRIVATE set_max_speed.F90 good_contour.f90
   ccw_orient.f90 write_eddy.f90 complete_ssh.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 get_var.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)

Inst_eddies/Tests/CMakeLists.txt

@@ -46,8 +46,17 @@ target_link_libraries(test_complete_ssh PRIVATE Jumble::jumble)
 # inst_eddies
 target_sources(inst_eddies PRIVATE cont_list.f90)
 
+# examine_eddies
+add_executable(examine_eddy examine_eddy.f90 cont_list.f90
+  ${PROJECT_SOURCE_DIR}/Overlap/read_eddy.f90)
+target_link_libraries(examine_eddy PRIVATE Contour_531::contour_531
+  Jumble::jumble Shapelib_03::shapelib_03 NetCDF95::netcdf95 gpc_f
+  Geometry::geometry Numer_Rec_95::numer_rec_95)
+target_compile_definitions(examine_eddy PRIVATE $<$<CONFIG:Debug>:DEBUG>)
+
 foreach(my_target IN ITEMS test_get_1_outerm test_good_contour test_mean_speed
-    test_nearby_extr test_local_extrema test_set_max_speed test_complete_ssh)
+    test_nearby_extr test_local_extrema test_set_max_speed test_complete_ssh
+    examine_eddy)
   set_target_properties(${my_target} PROPERTIES Fortran_MODULE_DIRECTORY
     ${CMAKE_CURRENT_BINARY_DIR}/${my_target}_modules)
   target_include_directories(${my_target} PRIVATE

Inst_eddies/Tests/examine_eddy.f90

+program examine_eddy
+
+  ! Libraries:
+  use contour_531, only: convert_to_ind
+  use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, pi, &
+       rad_to_deg, twopi
+  use shapelib, only: shpfileobject, shpclose
+  use shapelib_03, only: dbf_read_attribute_03, dbf_get_field_index_03, &
+       shp_read_point, shp_open_03
+
+  use config_m, only: config, min_radius
+  use cont_list_m, only: create_cont_list, close_cont_list
+  use derived_types, only: snapshot, shpc_slice_handler, eddy, missing_ssh
+  use input_ssh_m, only: input_ssh, max_radius
+  use read_eddy_m, only: read_eddy
+  use shpc_create_m, only: shpc_create
+  use set_all_extr_m, only: set_all_extr
+  use set_contours_m, only: set_contours
+  use shpc_close_m, only: shpc_close
+  use write_eddy_m, only: write_eddy
+
+  implicit none
+
+  integer:: date = huge(0), eddy_i_target = huge(0)
+  integer ishape_first, d0, ishape, i, l, n1, n2, number_extr
+  logical:: cyclone = .true. ! orientation of target eddy
+  namelist /main_nml/ cyclone, date, eddy_i_target
+  TYPE(shpfileobject) hshp
+  character(len = :), allocatable:: slice_dir, slice_o_dir
+  type(eddy) e
+  integer extr_date ! field index in DBF file
+
+  integer, allocatable:: ishape_last(:) ! (d0:)
+  ! shape index (0-based) in the collection of shapefiles of the last
+  ! shape at a given date index
+
+  real corner(2)
+  ! longitude and latitude of the corner of the whole grid, in rad
+
+  real step(2) ! longitude and latitude steps, in rad
+  logical periodic ! grid is periodic in longitude
+
+  real, allocatable:: ssh(:, :)
+  ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
+  ! in longitude, else (nlon, nlat). Sea-surface height, in m.
+
+  real, allocatable:: u(:, :), v(:, :)
+  ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+  ! (nlon, nlat). Wind, in m s-1.
+
+  type(snapshot) s
+
+  real, allocatable:: coord(:, :)
+  ! (2, number_extr) longitude and latitude, in rad
+
+  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
+
+  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
+  call input_ssh(corner, step, periodic, ssh, u, v)
+  call get_command_arg_dyn(3, slice_dir, &
+       "Required 3rd argument: input-slice-directory")
+  write(unit = *, nml = main_nml)
+  print *, "examine_eddy: Enter main_nml:"
+  read(unit = *, nml = main_nml)
+  write(unit = *, nml = main_nml)
+
+  if (cyclone) then
+     slice_o_dir = slice_dir // "/Cyclones"
+  else
+     slice_o_dir = slice_dir // "/Anticyclones"
+  end if
+
+  call shp_open_03(hshp, pszshapefile = slice_o_dir // "/extremum", &
+       pszaccess = "rb")
+  call dbf_get_field_index_03(hshp, "date", extr_date)
+  call dbf_read_attribute_03(d0, hshp, extr_date, ishape = 0)
+  call read_column(ishape_last, file = slice_o_dir // "/ishape_last.txt", &
+       my_lbound = d0)
+
+  if (date == d0) then
+     ishape_first = 0
+  else
+     ishape_first = ishape_last(date - 1) + 1
+  end if
+
+  number_extr = ishape_last(date) - ishape_first + 1
+  allocate(coord(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)
+  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%ssh = missing_ssh ! because we do not care about this
+  e%cyclone = cyclone
+
+  ! Define the geographical window around the target eddy extremum:
+
+  llc = e%extr%coord_proj - max_radius
+  urc = e%extr%coord_proj + max_radius
+
+  llc(2) = max(llc(2), 1)
+  urc(2) = min(urc(2), size(ssh, 2))
+
+  if (.not. periodic) then
+     llc(1) = max(llc(1), 1)
+     urc(1) = min(urc(1), size(ssh, 1))
+  end if
+
+  corner_window = corner + (llc - 1) * step
+
+  ! Done defining geographical window
+
+  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
+
+        in_window(i) = all(abs(nint((coord(:, i) - e%extr%coord) / step)) &
+             < max_radius)
+     else
+        in_window(i) = .false.
+     end if
+  end do
+
+  n2 = count(in_window)
+  call set_all_extr(s, step, periodic = .false., &
+       ssh = ssh(llc(1):urc(1), llc(2):urc(2)), corner = corner_window)
+  selection = pack(s%extr_map, s%extr_map /= 0)
+  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
+  l = n1 + 1
+
+  do i = 1, number_extr
+     if (in_window(i)) then
+        outside_points(:, l) = coord(:, i)
+        l = l + 1
+     end if
+  end do
+
+  i = s%extr_map(e%extr%coord_proj(1) - llc(1) + 1, &
+       e%extr%coord_proj(2) - llc(2) + 1)
+  e%innermost_level = s%list(i)%innermost_level
+  call create_cont_list
+  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)), &
+       corner_window, pi * (min_radius * 1e3)**2, outside_points)
+  call close_cont_list
+  call shpc_create(hshpc, shpc_dir = "SHPC", cyclone = cyclone, slice = 0, &
+       grid_lon_lat = .true.)
+  call write_eddy(e, hshpc, date, 1)
+  write(hshpc%unit, fmt = *) 0
+  CALL shpc_close(hshpc)
+  print *, "llc = ", llc
+  print *, "urc = ", urc
+  print *, "corner_window (in degrees):", corner_window * rad_to_deg
+  print *, "n1 = ", n1
+  print *, "n2 = ", n2
+  print *, "outside_points = ", outside_points
+  print *, "e%innermost_level = ", e%innermost_level
+  print *, "Created shapefile cont_list and shapefile collection SHPC."
+
+end program examine_eddy