test_set_max_speed.f90

program test_set_max_speed

  use, intrinsic:: ISO_FORTRAN_ENV
  use, intrinsic:: ISO_c_binding

  ! Libraries:
  use jumble, only: new_unit, count_lines
  use netcdf, only: nf90_nowrite
  use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
       find_coord, nf95_inquire_dimension
  use nr_util, only: deg_to_rad, rad_to_deg, pi
  use shapelib, only: shpt_polygon, shpfileobject, ftdouble, shpclose, &
       shpobject, dbfreadattribute
  use shapelib_03, only: shp_open_03, shp_create_03, dbf_add_field_03, &
       shp_append_object_03, dbf_write_attribute_03, shp_read_object_03

  use derived_types, only: eddy
  use set_max_speed_m, only: set_max_speed

  implicit none

  integer ncid, varid, dimid
  integer nlon, nlat, unit, n, l
  real, allocatable:: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
  real, allocatable:: u(:, :), v(:, :) ! (nlon, nlat) wind, in m s-1
  type(eddy) e

  integer:: ind_targ_extr(2) = [5, 4]
  ! indices of the target extremum, relative to corner

  real corner(2)
  TYPE(shpfileobject) hshp
  TYPE(shpobject) psobject
  integer ifield, ishape
  logical:: cyclone = .true.
  real:: ssh_extremum = 0.2759
  real:: coord_extr(2) = [9.625, - 33.875]
  real(kind = c_double) outermost_contour_ssh

  real, allocatable:: outside_points(:, :) ! (2, :) longitude and
  ! latitude, in rad, of all the extrema except the target extremum

  namelist /main_nml/ cyclone, ssh_extremum, coord_extr, ind_targ_extr

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

  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)

  print *, "Reading from h.nc..."
  call nf95_open("h.nc", nf90_nowrite, ncid)

  call find_coord(ncid, dimid = dimid, varid = varid, std_name = "longitude")
  call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
  call nf95_get_var(ncid, varid, corner(1))

  call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
  call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
  call nf95_get_var(ncid, varid, corner(2))

  allocate(ssh(nlon, nlat))
  call nf95_inq_varid(ncid, "adt", varid)
  call nf95_get_var(ncid, varid, ssh)

  call nf95_close(ncid)

  print *, "Reading from uv.nc..."
  call nf95_open("uv.nc", nf90_nowrite, ncid)

  allocate(u(nlon, nlat), v(nlon, nlat))

  call nf95_inq_varid(ncid, "ugos", varid)
  call nf95_get_var(ncid, varid, u)

  call nf95_inq_varid(ncid, "vgos", varid)
  call nf95_get_var(ncid, varid, v)

  call nf95_close(ncid)

  print *, "Reading from shapefile test_get_1_outerm..."
  call shp_open_03(hshp, "test_get_1_outerm", "rb")
  call shp_read_object_03(hshp, 0, psobject)
  call dbfreadattribute(hshp, 0, 0, outermost_contour_ssh)
  CALL shpclose(hshp)

  e%out_cont%ssh = outermost_contour_ssh
  e%cyclone = cyclone
  e%ssh_extr = ssh_extremum
  e%coord_extr = coord_extr * deg_to_rad
  e%out_cont%n_points = psobject%nvertices
  e%out_cont%closed = .true.
  allocate(e%out_cont%points(2, e%out_cont%n_points))
  e%out_cont%points(1, :) = psobject%padfx * deg_to_rad
  e%out_cont%points(2, :) = psobject%padfy * deg_to_rad

  print *, "Reading from outside_points.csv..."
  call new_unit(unit)
  open(unit, file = "outside_points.csv", status = "old", action = "read", &
       position = "rewind")
  call count_lines(unit, n)
  n = n - 2 ! 2 title lines
  allocate(outside_points(2, n))
  rewind(unit)
  read(unit, fmt = *)
  read(unit, fmt = *)

  do l = 1, n
     read(unit, fmt = *) outside_points(:, l)
  end do
  
  outside_points = outside_points * deg_to_rad
  close(unit)
  
  call set_max_speed(e, ind_targ_extr, outside_points, ssh, u, v, &
       corner * deg_to_rad, step = [0.25, 0.25]  * deg_to_rad)

  call shp_create_03("test_set_max_speed", shpt_polygon, hshp)
  call dbf_add_field_03(ifield, hshp, 'ssh', ftdouble, nwidth = 13, &
       ndecimals = 6)
  call dbf_add_field_03(ifield, hshp, 'r_eq_area', ftdouble, nwidth = 10, &
       ndecimals = 4)
  call dbf_add_field_03(ifield, hshp, 'speed', ftdouble, nwidth = 13, &
       ndecimals = 6)

  call shp_append_object_03(ishape, hshp, shpt_polygon, &
       e%speed_cont%points * rad_to_deg)
  call dbf_write_attribute_03(hshp, ishape, 0, e%speed_cont%ssh)

  if (e%speed_cont%area >= 0) then
     call dbf_write_attribute_03(hshp, ishape, 1, &
          sqrt(e%speed_cont%area / 1e6 / pi))
  else
     call dbf_write_attribute_03(hshp, ishape, 1, - 100.)
  end if
     
  call dbf_write_attribute_03(hshp, ishape, 2, e%max_speed)

  CALL shpclose(hshp)
  print *, 'Created shapefile "test_set_max_speed".'

end program test_set_max_speed