Skip to content
Snippets Groups Projects
Commit 0148e104 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Make the program test_set_all_outerm more to the point: testing the

procedure set_all_outerm. So do not define and output null max-speed
contours.

Bug fix in program test_set_all_outerm: take into account missing
values in ADT.

Make plot_snapshot.py work even if there is no max-speed contour
shapefile. So we can use plot_snapshot.py on results of program
test_set_all_outerm.
parent c0eb232f
No related branches found
No related tags found
No related merge requests found
......@@ -13,7 +13,7 @@ src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f go
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 set_max_speed.f outermost_possible_level.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spherical_polygon_area.f mean_speed.f inside_4.f set_all_outerm.f
src_test_set_all_outerm = test_set_all_outerm.f derived_types.f dispatch_snapshot.f set_all_outerm.f write_eddy.f send_snapshot.f local_extrema.f get_1_outerm.f good_contour.f spherical_polygon_area.f
src_test_set_all_outerm = test_set_all_outerm.f derived_types.f set_all_outerm.f local_extrema.f get_1_outerm.f good_contour.f spherical_polygon_area.f
sources := $(shell cat ${makefile_dir}/file_list)
......@@ -55,7 +55,7 @@ test_get_snapshot: ${obj_test_get_snapshot}
test_set_all_outerm: ${obj_test_set_all_outerm}
depend ${makefile_dir}/depend.mk:
makedepf90 -free -Wmissing -Wconfused $(addprefix -I, ${VPATH}) -nosrc $(addprefix -u , ${lib_list} shapelib netcdf) ${sources} >${makefile_dir}/depend.mk
makedepf90 -free -Wmissing -Wconfused $(addprefix -I, ${VPATH}) -nosrc $(addprefix -u , ${lib_list} shapelib netcdf intrinsic) ${sources} >${makefile_dir}/depend.mk
TAGS: ${sources}
ctags -e --language-force=fortran $^
......
......@@ -16,6 +16,7 @@ import matplotlib.pyplot as plt
import netCDF4
from mpl_toolkits import basemap
import math
import itertools
parser = argparse.ArgumentParser()
parser.add_argument("-v", "--velocity", help = "plot velocity field",
......@@ -40,11 +41,20 @@ m.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), latlon = True, marker = "+",
reader_extr = shapefile.Reader("extremum_1")
reader_outer = shapefile.Reader("outermost_contour_1")
reader_m_s = shapefile.Reader("max_speed_contour_1")
try:
reader_m_s = shapefile.Reader("max_speed_contour_1")
except shapefile.ShapefileException:
print("Shapefile max_speed_contour_1 not found. "
"Max-speed contours will not be plotted.")
w = shapefile.Writer(shapeType = shapefile.NULL)
m_s_iterShapes = itertools.repeat(w)
else:
m_s_iterShapes = reader_m_s.iterShapes()
for shape_rec_extr, shape_outer, shape_m_s \
in zip(reader_extr.iterShapeRecords(), reader_outer.iterShapes(),
reader_m_s.iterShapes()):
m_s_iterShapes):
points = shape_rec_extr.shape.points[0]
if shape_rec_extr.record[4] == 0:
......
......@@ -2,16 +2,16 @@ program test_set_all_outerm
use contour_531, only: null_polyline
use derived_types, only: snapshot
use dispatch_snapshot_m, only: dispatch_snapshot
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, &
find_coord, nf95_inquire_dimension
find_coord, nf95_inquire_dimension, nf95_get_att
use nr_util, only: pi
use set_all_outerm_m, only: set_all_outerm
use shapelib, only: shpt_point, shpt_polygon, shpfileobject, ftdouble, &
shpclose, ftinteger
use shapelib_03, only: shp_create_03, dbf_add_field_03
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_simple_object_03
implicit none
......@@ -33,9 +33,11 @@ program test_set_all_outerm
! grid points
real, parameter:: deg_over_rad = pi / 180.
real, parameter:: rad_over_deg = 180. / pi
integer ncid, varid, dimid, i
real, allocatable:: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
real Fill_Value
TYPE(shpfileobject) hshp_extremum
! shapefile extremum_$m. The fields in the DBF file are, in that
......@@ -46,15 +48,12 @@ program test_set_all_outerm
! shapefile outermost_contour_$m. The fields in the DBF file are,
! in that order: area, ssh, date index, eddy index, twice.
TYPE(shpfileobject) hshp_max_speed
! shapefile x_speed_contour_$m. The fields in the DBF file are, in
! that order: area, ssh, spped, date index, eddy index.
integer ifield, unit_isolated, unit_number_eddies
integer ifield, unit_isolated, ishape
!--------------------------------------------------------------
call get_command_arg_dyn(1, filename)
call get_command_arg_dyn(1, filename, &
"Required argument: NetCDF file containing adt")
print *, "min_amp = ? "
read *, min_amp
......@@ -75,6 +74,8 @@ program test_set_all_outerm
allocate(ssh(nlon, nlat))
call nf95_inq_varid(ncid, "adt", varid)
call nf95_get_var(ncid, varid, ssh)
call nf95_get_att(ncid, varid, "_FillValue", Fill_Value)
where (ssh == Fill_Value) ssh = huge(0.)
call nf95_close(ncid)
......@@ -83,15 +84,6 @@ program test_set_all_outerm
= [longitude(2) - longitude(1), latitude(2) - latitude(1)] &
* deg_over_rad, ssh = ssh)
do i = 1, s%number_vis_eddies
s%list_vis(i)%max_speed_contour%polyline = null_polyline()
s%list_vis(i)%max_speed_contour%ssh = s%list_vis(i)%ssh_extremum
s%list_vis(i)%max_speed_contour%area = 0.
s%list_vis(i)%max_speed = 0.
end do
s%number_eddies = s%number_vis_eddies
call shp_create_03("extremum_1", shpt_point, hshp_extremum)
call dbf_add_field_03(ifield, hshp_extremum, 'ssh', ftdouble, nwidth = 13, &
ndecimals = 6)
......@@ -118,44 +110,58 @@ program test_set_all_outerm
call dbf_add_field_03(ifield, hshp_outermost, 'twice', ftinteger, &
nwidth = 1, ndecimals = 0)
call shp_create_03("max_speed_contour_1", shpt_polygon, hshp_max_speed)
call dbf_add_field_03(ifield, hshp_max_speed, 'area', ftdouble, nwidth = 20, &
ndecimals = 6)
call dbf_add_field_03(ifield, hshp_max_speed, 'ssh', ftdouble, nwidth = 13, &
ndecimals = 6)
call dbf_add_field_03(ifield, hshp_max_speed, 'speed', &
ftdouble, nwidth = 13, ndecimals = 6)
call dbf_add_field_03(ifield, hshp_max_speed, 'date_index', ftinteger, &
nwidth = 4, ndecimals = 0)
call dbf_add_field_03(ifield, hshp_max_speed, 'eddy_index', ftinteger, &
nwidth = 5, ndecimals = 0)
call new_unit(unit_isolated)
open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
action = "write")
write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
call new_unit(unit_number_eddies)
open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
action = "write")
write(unit_number_eddies, fmt = *) '"date index" ' &
// '"number of visible eddies" "number of interpolated eddies"'
! (title line)
do i = 1, s%number_vis_eddies
call shp_append_point_03(ishape, hshp_extremum, &
s%list_vis(i)%coord_extr * rad_over_deg)
call dbf_write_attribute_03(hshp_extremum, ishape, 0, &
s%list_vis(i)%ssh_extremum)
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)%suff_amp))
if (s%list_vis(i)%interpolated) then
call shp_append_null_03(ishape, hshp_outermost)
else
if (s%list_vis(i)%outermost_contour%n_points == 0) then
call shp_append_null_03(ishape, hshp_outermost)
else
call shp_append_simple_object_03(ishape, hshp_outermost, &
shpt_polygon, &
s%list_vis(i)%outermost_contour%points * rad_over_deg)
end if
end if
call dbf_write_attribute_03(hshp_outermost, ishape, 0, &
s%list_vis(i)%outermost_contour%area)
call dbf_write_attribute_03(hshp_outermost, ishape, 1, &
s%list_vis(i)%outermost_contour%ssh)
call dbf_write_attribute_03(hshp_outermost, ishape, 2, 1)
call dbf_write_attribute_03(hshp_outermost, ishape, 3, i)
call dbf_write_attribute_03(hshp_outermost, ishape, 4, &
merge(1, 0, s%list_vis(i)%twice))
if (s%list_vis(i)%delta_in == huge(0) .and. s%list_vis(i)%delta_out &
== huge(0)) write(unit_isolated, fmt = *) 1, i
end do
call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
max_delta = 4)
print *, "s%number_vis_eddies = ", s%number_vis_eddies
CALL shpclose(hshp_extremum)
print *, 'Created shapefile "extremum_1".'
CALL shpclose(hshp_outermost)
print *, 'Created shapefile "outermost_contour_1".'
CALL shpclose(hshp_max_speed)
print *, 'Created shapefile "max_speed_contour_1".'
close(unit_isolated)
print *, 'Created "isolated_nodes_1.csv".'
close(unit_number_eddies)
print *, 'Created "number_eddies_1.csv".'
print *, "Average number of points per outermost contour: ", &
sum(s%list_vis%outermost_contour%n_points) / real(s%number_vis_eddies)
......
......@@ -145,12 +145,12 @@
"input": "1e-3\n"
},
{
"args": ["$compil_prod_dir/test_inside_4", "$input_dir/outermost_5"],
"args": ["$compil_prod_dir/test_inside_4", "$input_dir/outermost_eddy_5"],
"title": "Inside_4_true",
"stdin": "$stdin_dir/inside_4_true_nml.txt"
},
{
"args": ["$compil_prod_dir/test_inside_4", "$input_dir/outermost_5"],
"args": ["$compil_prod_dir/test_inside_4", "$input_dir/outermost_eddy_5"],
"title": "Inside_4_false",
"stdin": "$stdin_dir/inside_4_false_nml.txt"
},
......
......@@ -11,7 +11,7 @@ test_inside_4.o : inside_4.o
test_local_extrema.o : local_extrema.o
test_max_speed_contour_ssh.o : max_speed_contour_ssh.o
test_mean_speed.o : mean_speed.o
test_set_all_outerm.o : set_all_outerm.o dispatch_snapshot.o derived_types.o
test_set_all_outerm.o : set_all_outerm.o derived_types.o
test_set_max_speed.o : set_max_speed.o derived_types.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
write_eddy.o : derived_types.o
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment