test_mean_speed.f90

program test_mean_speed

  use, intrinsic:: ISO_FORTRAN_ENV

  ! Libraries:
  use gpc_f, only: shp_read_pol, polygon
  use jumble, only: get_command_arg_dyn
  use netcdf, only: nf90_nowrite
  use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
       nf95_find_coord, nf95_inquire_dimension
  use jumble, only: deg_to_rad, assert
  use shapelib, only: shpclose, shpfileobject
  use shapelib_03, only: shp_open_03
  
  use mean_speed_m, only: mean_speed

  implicit none

  real corner(2)
  real:: center(2) = [9.625, - 33.875] ! longitude and latitude, in degrees
  integer ncid, varid, dimid
  integer:: ishape = 0
  integer nlon, nlat
  real lon_max, lat_max ! longitude and latitude, in degrees
  real step(2) ! longitude and latitude steps, in degrees
  real, allocatable:: u(:, :), v(:, :) ! (nlon, nlat) wind, in m s-1
  TYPE(shpfileobject) hshp
  type(polygon) p
  character(len = :), allocatable:: contour_filename
  real min_lon_p ! minimum longitude of p, in degrees

  namelist /main_nml/ center, ishape

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

  call get_command_arg_dyn(1, contour_filename, "Required argument: shapefile")
  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)
  
  call nf95_open("uv.nc", nf90_nowrite, ncid)
  
  call nf95_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 nf95_get_var(ncid, varid, lon_max, start = [nlon])

  call nf95_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))
  call nf95_get_var(ncid, varid, lat_max, start = [nlat])

  call assert(lon_max > corner(1) .and. lat_max > corner(2), &
       "test_mean_speed coordinate order")
  step = [(lon_max - corner(1)) / (nlon - 1), &
       (lat_max - corner(2)) / (nlat - 1)]

  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)

  call shp_open_03(hshp, contour_filename, "rb")
  call shp_read_pol(p, hshp, ishape)
  CALL shpclose(hshp)

  ! Shift the longitude of the polyline so that it is in the grid:
  min_lon_p = minval(p%part(1)%points(1, :))
  p%part(1)%points(1, :) = p%part(1)%points(1, :) &
       + ceiling((corner(1) - min_lon_p) / 360.) * 360.
  
  if (corner(1) > center(1) .or. center(1) > lon_max &
       .or. minval(p%part(1)%points(1, :)) > center(1) &
       .or. center(1) > maxval(p%part(1)%points(1, :))) then
     print *, "corner(1) = ", corner(1)
     print *, "center(1) = ", center(1)
     print *, "lon_max = ", lon_max
     print *, "minval(p%part(1)%points(1, :)) = ", &
          minval(p%part(1)%points(1, :))
     print *, "maxval(p%part(1)%points(1, :)) = ", &
          maxval(p%part(1)%points(1, :))
     print *, "test_mean_speed: modulo 2 pi needed"
     stop 1
  end if
  
  p%part(1)%points = p%part(1)%points * deg_to_rad
  print *, "mean azimuthal speed = ", mean_speed(u, v, p%part(1), &
       center * deg_to_rad, corner * deg_to_rad, step  * deg_to_rad), "m s-1"

end program test_mean_speed