test_overlap.f90

program test_overlap

  use, intrinsic:: ISO_FORTRAN_ENV

  ! Libraries:
  use jumble, only: get_command_arg_dyn, read_column, new_unit
  use mpi_f08, only: mpi_init, mpi_finalize, MPI_Comm_rank, MPI_Comm_world, &
       MPI_Comm_size, mpi_abort
  use nr_util, only: deg_to_rad, assert
  use shapelib_03, only: dbf_read_attribute_03

  use derived_types, only: snapshot, shpc
  use overlap_m, only: overlap
  use read_snapshot_m, only: read_snapshot
  use shpc_close_m, only: shpc_close
  use shpc_create_m, only: shpc_create
  use shpc_open_m, only: shpc_open
  use unit_edge_m, only: open_edge_file, unit_edge
  use write_overlap_m, only: init_interpolated_eddy

  implicit none

  character(len = :), allocatable:: shpc_dir
  integer k1
  integer:: k_test_1 = 20454, k_test_2 = 20455
  integer unit, i, copy, rank, n_proc
  integer, allocatable:: ishape_last(:)
  type(snapshot), allocatable:: flow(:) ! (max_delta + 1)
  TYPE(shpc) hshp

  real:: corner_deg(2) = [huge(0.), huge(0.)], corner(2)
  ! longitude and latitude of the corner of the whole grid, in degrees
  ! and in rad

  real:: step_deg(2) = [0.25, 0.25], step(2)
  ! longitude and latitude steps, in degrees and in rad

  logical periodic ! grid is periodic in longitude
  integer:: nlon = - 1, nlat = - 1, max_delta = 1

  integer:: dist_lim = 12
  ! We look for an overlapping eddy at dist_lim (in grid points) of
  ! the first extremum.

  namelist /grid_nml/ corner_deg, step_deg, nlon, nlat
  namelist /main_nml/ dist_lim, max_delta, k_test_1, k_test_2

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

  call mpi_init
  call MPI_Comm_rank(MPI_Comm_world, rank)
  call mpi_comm_size(mpi_comm_world, n_proc)
  
  if (n_proc /= 1) then
     if (rank == 0) print *, "test_overlap: 1 process only"
     call mpi_abort(MPI_Comm_world, errorcode = 1)
  end if

  call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")
  call new_unit(unit)
  open(unit, file = shpc_dir // "/grid_nml.txt", status = "old", &
       action = "read", position = "rewind")
  read(unit, nml = grid_nml)
  close(unit)
  write(unit = error_unit, nml = main_nml)
  write(unit = error_unit, fmt = *) "Enter namelist main_nml."
  read(unit = *, nml = main_nml)
  write(unit = *, nml = main_nml)

  ! As we are requiring the grid spacing to be uniform, the value of
  ! "periodic" may be deduced from the values of step_deg(1) and nlon:
  periodic = nint(360. / step_deg(1)) == nlon
  print *, "periodic = ", periodic
  if (periodic) call assert(2 * dist_lim * step_deg(1) < 180., &
       "test_overlap dist_lim")
  copy = merge(dist_lim, 0, periodic)
  
  corner = corner_deg * deg_to_rad
  step = step_deg * deg_to_rad
  allocate(flow(max_delta + 1))
  call read_column(ishape_last, &
       file = trim(shpc_dir) // "/ishape_last.txt")
  call shpc_open(hshp, trim(shpc_dir), rank = 0)
  call dbf_read_attribute_03(k1, hshp%extremum, hshp%extr_date, ishape = 0)
  call read_snapshot(flow(1), k_test_1, hshp, corner, step, nlon, nlat, copy, &
       k1, ishape_last)
  call read_snapshot(flow(max_delta + 1), k_test_2, hshp, corner, step, nlon, &
       nlat, copy, k1, ishape_last)
  CALL shpc_close(hshp)
  print *, "Enter flow(1)%list_vis%delta_out (array with ", &
       flow(1)%number_vis_extr, " values):"
  read *, flow(1)%list_vis%delta_out
  print *, "Enter flow(max_delta + 1)%list_vis%delta_in (array with ", &
       flow(max_delta + 1)%number_vis_extr, " values)):"
  read *, flow(max_delta + 1)%list_vis%delta_in
  flow(2:max_delta)%number_eddies = 20000
  call open_edge_file(rank = 0)

  ! Title lines:
  write(unit_edge, fmt = "(1x, a)") '"predecessor date subscript" ' &
       // '"predecessor eddy subscript" "successor date subscript" ' &
       // '"successor eddy subscript"'
  write(unit_edge, fmt = *) "k1 i1 k2 i2 weight"

  call shpc_create(hshp, shpc_dir = "SHPC")
  call init_interpolated_eddy
  call overlap(flow, nlon, nlat, periodic, dist_lim, hshp, k = k_test_2, &
       delta = max_delta, j = max_delta + 1)
  close(unit_edge)
  print *, 'Created file "edgelist.csv".'
  CALL shpc_close(hshp)
  print *, 'Created shapefiles for interpolated eddies in SHPC.'

  print *, k_test_1, ":"
  print *, "Valid isolated eddies:"
  do i = 1, flow(1)%number_vis_extr
     if (flow(1)%list_vis(i)%valid .and. flow(1)%list_vis(i)%delta_out &
          == huge(0)) write(unit = *, fmt = "(i0, 1x)", advance = "no") i
  end do
  print *

  print *
  print *, k_test_2, ":"
  print *, "Valid isolated eddies:"

  do i = 1, flow(max_delta + 1)%number_vis_extr
     if (flow(max_delta + 1)%list_vis(i)%valid &
          .and. flow(max_delta + 1)%list_vis(i)%delta_in == huge(0)) &
          write(unit = *, fmt = "(i0, 1x)", advance = "no") i
  end do
  
  print *
  call mpi_finalize

end program test_overlap