Simplify test_max_speed_contour_ssh. The idea is to make unit tests:

call only one procedure per test, and read all the required
data. (Also: there was much duplication of code between
test_set_outermost_contour and test_max_speed_contour_ssh.)

GNUmakefile

@@ -9,7 +9,7 @@ src_test_local_extrema = test_local_extrema.f local_extrema.f
 
 src_test_set_outermost_contour = good_contour.f test_set_outermost_contour.f derived_types.f set_outermost_contour.f outermost_possible_level.f spherical_polygon_area.f local_extrema.f
 
-src_test_max_speed_contour_ssh = good_contour.f test_max_speed_contour_ssh.f max_speed_contour_ssh.f set_outermost_contour.f outermost_possible_level.f derived_types.f inside.f spherical_polygon_area.f local_extrema.f
+src_test_max_speed_contour_ssh = test_max_speed_contour_ssh.f max_speed_contour_ssh.f inside.f
 
 ##src_test_get_snapshot = test_get_snapshot.f get_snapshot.f dispatch_snapshot.f write_eddy.f send_snapshot.f receive_snapshot.f local_extrema.f get_eddy.f outermost_possible_level.f set_outermost_contour.f max_speed_contour_ssh.f good_contour.f spherical_polygon_area.f mean_speed.f inside.f
 

Sources/Tests/test_max_speed_contour_ssh.f

@@ -2,16 +2,10 @@ program test_max_speed_contour_ssh
 
   use, intrinsic:: ISO_FORTRAN_ENV
 
-  use derived_types, only: eddy
-  use local_extrema_m, only: local_extrema, construct_mask_center
   use max_speed_contour_ssh_m, only: max_speed_contour_ssh
   use netcdf, only: nf90_nowrite
   use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var
-  use nr_util, only: pi
-  use set_outermost_contour_m, only: set_outermost_contour
-  use shapelib, only: shpt_polygon, shpfileobject, ftdouble, shpclose
-  use shapelib_03, only: shp_create_03, dbf_add_field_03, &
-       shp_append_simple_object_03, dbf_write_attribute_03
+  use nr_util, only: assert
 
   implicit none
 
@@ -19,59 +13,55 @@ program test_max_speed_contour_ssh
   integer:: ilon_urc = 49, ilat_urc = 231 ! upper right corner
   integer n_lon, n_lat
 
+  integer length, status
+  character(len = :), allocatable:: adt_file, velocity_file
   integer ncid, varid
-  real, allocatable:: longitude(:) ! (n_lon) in rad
-  real, allocatable:: latitude(:) ! (n_lat) in rad
   real, allocatable:: ssh(:, :) ! (n_lon, n_lat) sea-surface height, in m
   real, allocatable:: u(:, :), v(:, :) ! (n_lon, n_lat) wind, in m s-1
-  integer, allocatable:: extr_map(:, :) ! (n_lon, n_lat) map of extrema
-
-  integer, allocatable:: ind_extr(:, :) ! (2, n_extr)
-  ! indices in the two dimensions of extrema
+  integer:: ind_extr(2) = [19, 11] ! indices in the two dimensions of
+                                   ! the extremum
+  real:: ssh_outermost = 0.460744 ! in m
+  logical:: cyclone = .true.
+  real level
 
-  real, allocatable:: innermost_level(:) ! (n_extr)
-  ! level of innermost contour, for each extremum
+  namelist /main_nml/ ilon_llc, ilat_llc, ilon_urc, ilat_urc, ind_extr, &
+       ssh_outermost, cyclone
 
-  logical, allocatable:: cyclone(:) ! (n_extr)
+  !----------------------------------------------------------------
 
-  type(eddy) e
-  real, parameter:: step = 0.25 / 180. * pi ! in rad
-  TYPE(shpfileobject) shphandle
-  integer field_number, shape_number
-  real level
-  integer:: i = 6 ! identifying number of the target extremum
+  call assert(COMMAND_ARGUMENT_COUNT() == 2, &
+       "Required arguments: ADT-file velocity-file")
 
-  namelist /main_nml/ ilon_llc, ilat_llc, ilon_urc, ilat_urc, i
+  call get_command_argument(number = 1, length = length, status = status)
+  allocate(character(len = length):: adt_file)
+  call get_command_argument(1, adt_file)
 
-  !----------------------------------------------------------------
+  call get_command_argument(number = 2, length = length, status = status)
+  allocate(character(len = length):: velocity_file)
+  call get_command_argument(2, velocity_file)
 
   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)
+  
   n_lon = ilon_urc - ilon_llc + 1
   n_lat = ilat_urc - ilat_llc + 1
-  allocate(longitude(n_lon), latitude(n_lat))
-  allocate(ssh(n_lon, n_lat), extr_map(n_lon, n_lat))
-  allocate(u(n_lon, n_lat), v(n_lon, n_lat))
-
-  call nf95_open("h_2006_01_01.nc", nf90_nowrite, ncid)
-  
-  call nf95_inq_varid(ncid, "lon", varid)
-  call nf95_get_var(ncid, varid, longitude, start = [ilon_llc])
-  longitude = longitude / 180. * pi ! degrees east to rad
+  allocate(ssh(n_lon, n_lat))
 
-  call nf95_inq_varid(ncid, "lat", varid)
-  call nf95_get_var(ncid, varid, latitude, start = [ilat_llc])
-  latitude = latitude / 180. * pi ! degrees north to rad
+  print *, "Reading from ", adt_file, "..."
+  call nf95_open(adt_file, nf90_nowrite, ncid)
 
   call nf95_inq_varid(ncid, "adt", varid)
   call nf95_get_var(ncid, varid, ssh, start = [ilon_llc, ilat_llc, 1])
 
   call nf95_close(ncid)
-  
-  call nf95_open("uv_2006_01_01.nc", nf90_nowrite, ncid)
-  
+
+  allocate(u(n_lon, n_lat), v(n_lon, n_lat))
+
+  print *, "Reading from ", velocity_file, "..."
+  call nf95_open(velocity_file, nf90_nowrite, ncid)
+
   call nf95_inq_varid(ncid, "u", varid)
   call nf95_get_var(ncid, varid, u, start = [ilon_llc, ilat_llc, 1])
 
@@ -80,33 +70,9 @@ program test_max_speed_contour_ssh
 
   call nf95_close(ncid)
 
-  call construct_mask_center
-  call local_extrema(ssh, extr_map, ind_extr, innermost_level, cyclone)
+  level = max_speed_contour_ssh(ssh, u, v, ind_extr, ssh_outermost, cyclone)
 
-  e%coord_extr = [longitude(ind_extr(1, i)), latitude(ind_extr(2, i))]
-  e%cyclone = cyclone(i)
-  call set_outermost_contour(e, ind_extr(:, i), innermost_level(i), extr_map, &
-       ssh, corner = [longitude(1), latitude(1)], step = [step, step], &
-       noise_around = .false.)
-  
-  print *, "e%outermost_contour%closed = ", e%outermost_contour%closed
-  print *, "Radius of disk of equal area: ", &
-       sqrt(e%outermost_contour%area / pi) / 1e3, "km"
-
-  call shp_create_03("test_max_speed_contour_ssh", shpt_polygon, shphandle)
-  call dbf_add_field_03(field_number, shphandle, 'level', ftdouble, &
-       nwidth = 13, ndecimals = 6)
-  call shp_append_simple_object_03(shape_number, shphandle, shpt_polygon, &
-       e%outermost_contour%points / pi * 180.)
-  call dbf_write_attribute_03(shphandle, shape_number, ifield = 0, &
-       fieldvalue = e%outermost_contour%ssh)
-  CALL shpclose(shphandle)
-  print *, 'Created shapefile "test_max_speed_contour_ssh".'
-
-  level = max_speed_contour_ssh(ssh, u, v, ind_extr(:, i), &
-       e%outermost_contour%ssh, e%cyclone)
-
-  if (level == ssh(ind_extr(1, i), ind_extr(2, i))) then
+  if (level == ssh(ind_extr(1), ind_extr(2))) then
      print *, "No max speed contour because outermost contour is within " &
           // "one grid point of extremum."
   else

Sources/Tests/test_set_outermost_contour.f

@@ -85,7 +85,7 @@ program test_set_outermost_contour
   call set_outermost_contour(e, ind_extr(:, i), innermost_level(i), extr_map, &
        ssh, corner = [longitude(1), latitude(1)], step = [step, step], &
        noise_around = noise_around)
-  
+
   print *, "e%outermost_contour%closed = ", e%outermost_contour%closed
   print *, "Radius of disk of equal area: ", &
        sqrt(e%outermost_contour%area / pi) / 1e3, "km"

depend.mk

@@ -6,7 +6,7 @@ send_snapshot.o : derived_types.o
 test_good_contour.o : good_contour.o 
 test_inside.o : inside.o 
 test_local_extrema.o : local_extrema.o 
-test_max_speed_contour_ssh.o : outermost_possible_level.o set_outermost_contour.o max_speed_contour_ssh.o local_extrema.o derived_types.o 
+test_max_speed_contour_ssh.o : max_speed_contour_ssh.o 
 test_mean_speed.o : mean_speed.o 
 test_set_outermost_contour.o : set_outermost_contour.o local_extrema.o derived_types.o 
 write_eddy.o : derived_types.o