test_local_extrema.f90

program test_local_extrema

  ! Reads ADT from a NetCDF file and writes the list of extrema to a
  ! CSV file. Also creates a NetCDF file containing the map of extrema.

  ! Note: the grid spacing does not need to be uniform.

  ! Libraries:
  use jumble, only: new_unit, get_command_arg_dyn
  use netcdf, only: nf90_nowrite
  use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
       nf95_gw_var, nf95_get_att, find_coord

  use local_extrema_m, only: local_extrema
  use write_extr_map_m, only: write_extr_map

  implicit none

  integer nlon, nlat
  integer ncid, varid
  integer i, unit
  character(len = :), allocatable:: filename
  real, allocatable:: longitude(:) ! (nlon) in degrees
  real, allocatable:: latitude(:) ! (nlat) in degrees

  real, allocatable:: ssh(:, :) ! (1 - copy:nlon + copy, nlat)
  ! sea-surface height, in m
  
  integer, allocatable:: extr_map(:, :) ! (1 - copy:nlon + copy, nlat)
  ! map of extrema

  integer, allocatable:: ind_extr(:, :) ! (2, n_extr)
  ! indices in the two dimensions of extrema

  real, allocatable:: innermost_level(:) ! (n_extr)
  ! level of innermost contour, for each extremum

  logical, allocatable:: cyclone(:) ! (n_extr)
  real Fill_Value
  logical periodic ! in longitude
  integer copy

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

  call get_command_arg_dyn(1, filename, message = "Required argument: ADT file")
  print *, "periodic = ? "
  read *, periodic
  print *, "periodic = ", periodic
  print *, "Reading from ", filename, "..."
  call nf95_open(filename, nf90_nowrite, ncid)
  
  call find_coord(ncid, std_name = "longitude", varid = varid)
  call nf95_gw_var(ncid, varid, longitude)

  call find_coord(ncid, std_name = "latitude", varid = varid)
  call nf95_gw_var(ncid, varid, latitude)

  nlon = size(longitude)
  nlat = size(latitude)
  copy = merge(1, 0, periodic)
  allocate(ssh(1 - copy:nlon + copy, nlat), &
       extr_map(1 - copy:nlon + copy, nlat))

  call nf95_inq_varid(ncid, "adt", varid)
  call nf95_get_var(ncid, varid, ssh(1:nlon, :))
  call nf95_get_att(ncid, varid, "_FillValue", Fill_Value)
  where (ssh(1:nlon, :) == Fill_Value) ssh(1:nlon, :) = huge(0.)

  if (periodic) then
     ssh(0, :) = ssh(nlon, :)
     ssh(nlon + 1, :) = ssh(1, :)
  end if

  call nf95_close(ncid)

  call local_extrema(extr_map, ind_extr, innermost_level, cyclone, ssh, &
       periodic, my_lbound = [1 - copy, 1])

  call new_unit(unit)
  open(unit, file = "test_local_extrema.csv", status = "replace", &
       action = "write")

  ! Title lines:
  write(unit, fmt = *) '"longitude index" "latitude index"'
  write(unit, fmt = *) '"" "" "" "degrees east" "degrees north" "m" "m"'
  write(unit, fmt = *) 'ilon ilat cyclone longitude latitude ssh ', &
       'innermost_level'

  do i = 1, size(ind_extr, 2)
     write(unit, fmt = *) ind_extr(1, i), ind_extr(2, i), cyclone(i), &
          longitude(ind_extr(1, i)), latitude(ind_extr(2, i)), &
          ssh(ind_extr(1, i), ind_extr(2, i)), innermost_level(i)
  end do
  close(unit)
  print *, 'Created file "test_local_extrema.csv".'

    if (periodic) then
       call write_extr_map(extr_map, &
            [longitude(nlon) - 360., longitude, longitude(1) + 360.], latitude)
    else
       call write_extr_map(extr_map, longitude, latitude)
    end if

end program test_local_extrema