test_overlap.f90

program test_overlap

  ! Libraries:
  use jumble, only: get_command_arg_dyn, read_opcol, new_unit, ediff1d, assert
  use shapelib_03, only: dbf_read_attribute_03

  use config_graph_m, only: config_graph, copy, max_delta, cyclone
  use derived_types, only: snapshot, shpc_slice_handler, shpc_slice_meta
  use overlap_m, only: overlap
  use read_grid_m, only: read_grid, uniform_lon_lat
  use read_snapshot_m, only: read_snapshot
  use shpc_close_m, only: shpc_close
  use shpc_open_m, only: shpc_open
  use unit_edge_m, only: open_edge_file, unit_edge

  implicit none

  character(len = :), allocatable:: shpc_dir
  integer:: k_test_1 = 0, k_test_2 = 1, i_slice = 0
  integer unit, i, n_dates
  type(snapshot), allocatable:: flow(:) ! (max_delta + 1)
  TYPE(shpc_slice_handler) hshp
  type(shpc_slice_meta) ssm
  integer e_overestim ! over-estimation of the number of eddies at each date

  namelist /main_nml/ k_test_1, k_test_2, i_slice

  !-------------------------------------------------------------------------

  call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")
  call read_grid(shpc_dir, rank = 0)
  call config_graph(rank = 0)

  ! main_nml:
  write(unit = *, nml = main_nml)
  print *, "Enter namelist main_nml."
  read(unit = *, nml = main_nml)
  write(unit = *, nml = main_nml)

  call shpc_open(hshp, trim(shpc_dir), cyclone, i_slice, &
       with_proj = .not. uniform_lon_lat, pszaccess = "rb")
  call dbf_read_attribute_03(ssm%d0, hshp%extremum, hshp%extr_date, ishape = 0)
  call read_opcol(ssm%ishape_last, hshp%unit, my_lbound = ssm%d0)
  n_dates = size(ssm%ishape_last)
  call assert(ssm%d0 <= [k_test_1, k_test_2] .and. [k_test_1, k_test_2] &
       < ssm%d0 + n_dates, "test_overlap k_test_1, k_test_2")
  e_overestim = maxval([ssm%ishape_last(ssm%d0) + 1, ediff1d(ssm%ishape_last)])
  call new_unit(unit)
  open(unit, file = "e_overestim.txt", status = "replace", action = "write")
  write(unit, fmt = *) e_overestim
  close(unit)
  allocate(flow(max_delta + 1))
  call read_snapshot(flow(1), [hshp], [ssm], k_test_1, copy)
  call read_snapshot(flow(max_delta + 1), [hshp], [ssm], k_test_2, copy)
  CALL shpc_close(hshp)
  print *, "Enter flow(1)%list%delta_out (array with ", &
       flow(1)%number_extr, " values):"
  read *, flow(1)%list%delta_out
  print *, "Enter flow(max_delta + 1)%list%delta_in (array with ", &
       flow(max_delta + 1)%number_extr, " values)):"
  read *, flow(max_delta + 1)%list%delta_in
  call open_edge_file(rank = 0)
  call overlap(flow, e_overestim, k = k_test_2, delta = max_delta, &
       j = max_delta + 1)
  close(unit_edge)
  print *, 'Created file "edgelist.csv".'

  print *, k_test_1, ":"
  print *, "Isolated eddies:"
  do i = 1, flow(1)%number_extr
     if (flow(1)%list(i)%delta_out == huge(0)) &
          write(unit = *, fmt = "(i0, 1x)", advance = "no") i
  end do
  print *

  print *
  print *, k_test_2, ":"
  print *, "Isolated eddies:"

  do i = 1, flow(max_delta + 1)%number_extr
     if (flow(max_delta + 1)%list(i)%delta_in == huge(0)) &
          write(unit = *, fmt = "(i0, 1x)", advance = "no") i
  end do

  print *

end program test_overlap