Use higher level procedures

Use higher level procedures in main program unit `test_set_all_outerm`.

Inst_eddies/Tests/CMakeLists.txt

@@ -17,7 +17,10 @@ add_executable(test_set_all_outerm
   ${CMAKE_SOURCE_DIR}/derived_types.f90 set_all_outerm.f90
   local_extrema.f90 get_1_outerm.f90 good_contour.f90
   ${CMAKE_SOURCE_DIR}/spher_polyline_area.f90 nearby_extr.f90
-  get_var.f90 ${CMAKE_CURRENT_LIST_DIR}/test_set_all_outerm.f90 config.f90)
+  get_var.f90 ${CMAKE_CURRENT_LIST_DIR}/test_set_all_outerm.f90
+  config.f90 ${CMAKE_SOURCE_DIR}/shp_tr_create.f90
+  ${CMAKE_SOURCE_DIR}/write_eddy.f90
+  ${CMAKE_SOURCE_DIR}/shp_tr_close.f90)
 
 target_link_libraries(test_set_all_outerm geometry NetCDF95::netcdf95
   numer_rec_95 shapelib_03 contour_531 jumble nr_util

Inst_eddies/Tests/test_set_all_outerm.f90

@@ -8,15 +8,14 @@ program test_set_all_outerm
   use netcdf95, only: nf95_open, nf95_close, nf95_get_var, find_coord, &
        nf95_inquire_dimension
   use nr_util, only: pi, assert, deg_to_rad, rad_to_deg, twopi
-  use shapelib, only: shpt_point, shpt_polygon, shpfileobject, ftdouble, &
-       shpclose, ftinteger
-  use shapelib_03, only: shp_create_03, dbf_add_field_03, shp_append_point_03, &
-       dbf_write_attribute_03, shp_append_null_03, shp_append_object_03
   
   use config_m, only: config, max_radius_deg, min_radius
-  use derived_types, only: snapshot
+  use derived_types, only: snapshot, shp_tr, null_ssh_contour, missing_speed
   use get_var_m, only: get_var
   use set_all_outerm_m, only: set_all_outerm
+  use shp_tr_close_m, only: shp_tr_close
+  use shp_tr_create_m, only: shp_tr_create
+  use write_eddy_m, only: write_eddy
 
   implicit none
 
@@ -35,11 +34,9 @@ program test_set_all_outerm
   ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
   ! in longitude, else (nlon, nlat). Sea-surface height, in m.
 
-  TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
-  TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
-  integer ifield, ishape
   real step_deg(2), step(2) ! longitude and latitude steps, in degrees and rad
   logical periodic ! grid is periodic in longitude
+  TYPE(shp_tr) hshp
 
   !--------------------------------------------------------------
 
@@ -86,75 +83,22 @@ program test_set_all_outerm
        corner = [lon_min, lat_min] * deg_to_rad, &
        min_area = pi * (min_radius * 1e3)**2)
 
-  call shp_create_03("SHP_triplet/extremum", shpt_point, hshp_extremum)
-  call dbf_add_field_03(ifield, hshp_extremum, 'ssh', ftdouble, nwidth = 13, &
-       ndecimals = 6)
-  call dbf_add_field_03(ifield, hshp_extremum, 'date_index', ftinteger, &
-       nwidth = 4, ndecimals = 0)
-  call dbf_add_field_03(ifield, hshp_extremum, 'eddy_index', ftinteger, &
-       nwidth = 5, ndecimals = 0)
-  call dbf_add_field_03(ifield, hshp_extremum, 'interpolat', ftinteger, &
-       nwidth = 1, ndecimals = 0)
-  call dbf_add_field_03(ifield, hshp_extremum, 'cyclone', ftinteger, &
-       nwidth = 1, ndecimals = 0)
-  call dbf_add_field_03(ifield, hshp_extremum, 'valid', ftinteger, &
-       nwidth = 1, ndecimals = 0)
-
-  call shp_create_03("SHP_triplet/outermost_contour", shpt_polygon, &
-       hshp_outermost)
-  call dbf_add_field_03(ifield, hshp_outermost, 'r_eq_area', ftdouble, &
-       nwidth = 10, ndecimals = 4)
-  call dbf_add_field_03(ifield, hshp_outermost, 'ssh', ftdouble, nwidth = 13, &
-       ndecimals = 6)
-  call dbf_add_field_03(ifield, hshp_outermost, 'date_index', ftinteger, &
-       nwidth = 4, ndecimals = 0)
-  call dbf_add_field_03(ifield, hshp_outermost, 'eddy_index', ftinteger, &
-       nwidth = 5, ndecimals = 0)
+  do i = 1, s%number_vis_extr
+     s%list_vis(i)%speed_cont = null_ssh_contour()
+     s%list_vis(i)%max_speed = missing_speed
+     s%list_vis(i)%radius4 = 0
+  end do
+
+  call shp_tr_create(hshp, shp_tr_dir = "SHP_triplet")
 
   do i = 1, s%number_vis_extr
-     call shp_append_point_03(ishape, hshp_extremum, &
-          s%list_vis(i)%coord_extr * rad_to_deg)
-     call dbf_write_attribute_03(hshp_extremum, ishape, 0, &
-          s%list_vis(i)%ssh_extr)
-     call dbf_write_attribute_03(hshp_extremum, ishape, 1, 1)
-     call dbf_write_attribute_03(hshp_extremum, ishape, 2, i)
-     call dbf_write_attribute_03(hshp_extremum, ishape, 3, &
-          merge(1, 0, s%list_vis(i)%interpolated))
-     call dbf_write_attribute_03(hshp_extremum, ishape, 4, &
-          merge(1, 0, s%list_vis(i)%cyclone))
-     call dbf_write_attribute_03(hshp_extremum, ishape, 5, &
-          merge(1, 0, s%list_vis(i)%valid))
-
-     if (s%list_vis(i)%interpolated) then
-        call shp_append_null_03(ishape, hshp_outermost)
-     else
-        if (s%list_vis(i)%out_cont%n_points == 0) then
-           call shp_append_null_03(ishape, hshp_outermost)
-        else
-           call shp_append_object_03(ishape, hshp_outermost, shpt_polygon, &
-                s%list_vis(i)%out_cont%points * rad_to_deg)
-        end if
-     end if
-
-     if (s%list_vis(i)%out_cont%area >= 0) then
-        call dbf_write_attribute_03(hshp_outermost, ishape, 0, &
-             sqrt(s%list_vis(i)%out_cont%area / 1e6 / pi))
-     else
-        call dbf_write_attribute_03(hshp_outermost, ishape, 0, - 100.)
-     end if
-      
-     call dbf_write_attribute_03(hshp_outermost, ishape, 1, &
-          s%list_vis(i)%out_cont%ssh)
-     call dbf_write_attribute_03(hshp_outermost, ishape, 2, 1)
-     call dbf_write_attribute_03(hshp_outermost, ishape, 3, i)
+     call write_eddy(s%list_vis(i), hshp, k = 1, i = i)
   end do
 
   print *, "s%number_vis_extr = ", s%number_vis_extr
 
-  CALL shpclose(hshp_extremum)
-  print *, 'Created shapefile "SHP_triplet/extremum".'
-  CALL shpclose(hshp_outermost)
-  print *, 'Created shapefile "SHP_triplet/outermost_contour".'
+  CALL shp_tr_close(hshp)
+  print *, 'Created shapefiles in SHP_triplet.'
   print *, "Average number of points per outermost contour: ", &
        sum(s%list_vis%out_cont%n_points) / real(s%number_vis_extr)