From 2b64791c1672a5b963cbaff38e3e4a6ce3b97e56 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 14 Mar 2019 17:55:47 +0100
Subject: [PATCH] Write number of days since 1950-1-1 as date_index in
shapefiles (instead of date_index = 1). This is the value read from the
NetCDF files. We have to increase the width of field date_index to 5 digits
(4 digits would only allow approximately 30 years since 1950).
In procedure read_snapshot, check that all input eddies have the same
date index, and return it as an additional argument.
Assume that the input NetCDF files contain a single date.
In program test_set_all_outerm, do not write file
"isolated_nodes_1.csv". Since s is defined by set_all_outerm only, all
the eddies are necessarily isolated.
---
Tests/read_snapshot.f | 24 +++++++++++++-----------
Tests/test_max_speed_contour_ssh.f | 6 ++----
Tests/test_read_snapshot.f | 6 +++---
Tests/test_set_all_outerm.f | 22 ++++++----------------
Tests/test_successive_overlap.f | 6 +++---
extraction_eddies.f | 25 ++++++++++++++++++-------
get_snapshot.f | 3 ++-
get_var.f | 5 ++---
init_shapefiles.f | 6 +++---
9 files changed, 52 insertions(+), 51 deletions(-)
diff --git a/Tests/read_snapshot.f b/Tests/read_snapshot.f
index aefc3bc6..ca1795a6 100644
--- a/Tests/read_snapshot.f
+++ b/Tests/read_snapshot.f
@@ -4,8 +4,8 @@ module read_snapshot_m
contains
- subroutine read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
- corner, nlon, nlat)
+ subroutine read_snapshot(s, k, hshp_extremum, hshp_outermost, &
+ hshp_max_speed, corner, nlon, nlat)
! Libraries:
use contour_531, only: convert_to_ind
@@ -18,6 +18,7 @@ contains
use read_field_indices_m, only: read_field_indices
type(snapshot), intent(out):: s ! completely defined
+ integer, intent(out):: k ! date index
TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum_1
TYPE(shpfileobject), intent(inout):: hshp_outermost
@@ -32,13 +33,8 @@ contains
integer, intent(in):: nlon, nlat
! Local:
-
- integer ishape
-
- integer k
- ! Date index, not used. We could check that it has the same value
- ! for all input shapes.
-
+ integer ishape
+ integer k2 ! date index
integer i ! eddy index
type(eddy) e
@@ -48,9 +44,15 @@ contains
allocate(s%list_vis(s%number_vis_extr))
call read_field_indices(hshp_extremum, hshp_outermost, hshp_max_speed)
- do ishape = 0, s%number_vis_extr - 1
- call read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ ! The first shape gives the date index:
+ call read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ ishape = 0)
+ s%list_vis(i) = e
+
+ do ishape = 1, s%number_vis_extr - 1
+ call read_eddy(e, k2, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
ishape)
+ call assert(k2 == k, "read_snapshot: date index")
s%list_vis(i) = e
end do
diff --git a/Tests/test_max_speed_contour_ssh.f b/Tests/test_max_speed_contour_ssh.f
index 08ade81c..aa3669e2 100644
--- a/Tests/test_max_speed_contour_ssh.f
+++ b/Tests/test_max_speed_contour_ssh.f
@@ -54,11 +54,9 @@ program test_max_speed_contour_ssh
print *, "Reading from ", velocity_file, "..."
call nf95_open(velocity_file, nf90_nowrite, ncid)
call get_var(periodic = .false., max_rad_lon = 0, values = u, ncid = ncid, &
- nlon = nlon, k = 1, name = "u", &
- new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+ nlon = nlon, name = "u", new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
call get_var(periodic = .false., max_rad_lon = 0, values = v, ncid = ncid, &
- nlon = nlon, k = 1, name = "v", &
- new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+ nlon = nlon, name = "v", new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
call nf95_close(ncid)
print *, "level = ", max_speed_contour_ssh(ssh(:, :, 1), u, v, ind_extr, &
diff --git a/Tests/test_read_snapshot.f b/Tests/test_read_snapshot.f
index 2c6232b8..9ef0336a 100644
--- a/Tests/test_read_snapshot.f
+++ b/Tests/test_read_snapshot.f
@@ -18,7 +18,7 @@ program test_read_snapshot
TYPE(shpfileobject) hshp_extremum ! shapefile extremum
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour
TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour
- integer i
+ integer i, k
real:: corner(2) = [0.125, - 59.875]
! longitude and latitude of the corner of the whole grid, in degrees
@@ -38,7 +38,7 @@ program test_read_snapshot
hshp = hshp_outermost)
call shp_open_03("max_speed_contour_old", pszaccess = "rb", &
hshp = hshp_max_speed)
- call read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ call read_snapshot(s, k, hshp_extremum, hshp_outermost, hshp_max_speed, &
corner * deg_to_rad, nlon, nlat)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
@@ -48,7 +48,7 @@ program test_read_snapshot
! Write snapshot:
do i = 1, s%number_vis_extr
- call write_eddy(s%list_vis(i), 1, i, hshp_extremum, hshp_outermost, &
+ call write_eddy(s%list_vis(i), k, i, hshp_extremum, hshp_outermost, &
hshp_max_speed)
end do
diff --git a/Tests/test_set_all_outerm.f b/Tests/test_set_all_outerm.f
index 5d241b6b..32e89ab8 100644
--- a/Tests/test_set_all_outerm.f
+++ b/Tests/test_set_all_outerm.f
@@ -44,7 +44,7 @@ program test_set_all_outerm
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
- integer ifield, unit_isolated, ishape
+ integer ifield, ishape
real step(2) ! longitude and latitude steps, in rad
logical periodic ! grid is periodic in longitude
namelist /main_nml/ min_amp, max_radius, min_radius
@@ -87,7 +87,7 @@ program test_set_all_outerm
copy = merge(max_radius(1), 0, periodic)
allocate(ssh(1 - copy:nlon + copy, nlat))
- call get_var(periodic, max_radius(1), ssh, ncid, nlon, k = 1, name = "adt", &
+ call get_var(periodic, max_radius(1), ssh, ncid, nlon, name = "adt", &
new_fill_value = huge(0.))
call nf95_close(ncid)
@@ -96,7 +96,7 @@ program test_set_all_outerm
corner = [lon_min, lat_min] * deg_to_rad, &
min_area = pi * (min_radius * 1e3)**2)
- call shp_create_03("extremum_1", shpt_point, hshp_extremum)
+ call shp_create_03("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, &
@@ -110,7 +110,7 @@ program test_set_all_outerm
call dbf_add_field_03(ifield, hshp_extremum, 'valid', ftinteger, &
nwidth = 1, ndecimals = 0)
- call shp_create_03("outermost_contour_1", shpt_polygon, hshp_outermost)
+ call shp_create_03("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, &
@@ -120,11 +120,6 @@ program test_set_all_outerm
call dbf_add_field_03(ifield, hshp_outermost, '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
-
do i = 1, s%number_vis_extr
call shp_append_point_03(ishape, hshp_extremum, &
s%list_vis(i)%coord_extr * rad_to_deg)
@@ -161,19 +156,14 @@ program test_set_all_outerm
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)
-
- 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
print *, "s%number_vis_extr = ", s%number_vis_extr
CALL shpclose(hshp_extremum)
- print *, 'Created shapefile "extremum_1".'
+ print *, 'Created shapefile "extremum".'
CALL shpclose(hshp_outermost)
- print *, 'Created shapefile "outermost_contour_1".'
- close(unit_isolated)
- print *, 'Created "isolated_nodes_1.csv".'
+ print *, 'Created shapefile "outermost_contour".'
print *, "Average number of points per outermost contour: ", &
sum(s%list_vis%out_cont%n_points) / real(s%number_vis_extr)
diff --git a/Tests/test_successive_overlap.f b/Tests/test_successive_overlap.f
index 0fefd796..147808c7 100644
--- a/Tests/test_successive_overlap.f
+++ b/Tests/test_successive_overlap.f
@@ -14,7 +14,7 @@ program test_successive_overlap
implicit none
- integer unit_edgelist, i
+ integer unit_edgelist, i, k
type(snapshot) flow(2)
TYPE(shpfileobject) hshp_extremum ! shapefile extremum
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour
@@ -38,8 +38,8 @@ program test_successive_overlap
hshp = hshp_outermost)
call shp_open_03("max_speed_contour", pszaccess = "rb", &
hshp = hshp_max_speed)
- call read_snapshot(flow(1), hshp_extremum, hshp_outermost, hshp_max_speed, &
- corner * deg_to_rad, nlon, nlat)
+ call read_snapshot(flow(1), k, hshp_extremum, hshp_outermost, &
+ hshp_max_speed, corner * deg_to_rad, nlon, nlat)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
CALL shpclose(hshp_max_speed)
diff --git a/extraction_eddies.f b/extraction_eddies.f
index 70e650c8..3a16b527 100644
--- a/extraction_eddies.f
+++ b/extraction_eddies.f
@@ -8,7 +8,7 @@ program extraction_eddies
use contour_531, only: convert_to_ind
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
- nf95_get_var, nf95_close
+ nf95_get_var, nf95_close, nf95_get_att
use nr_util, only: assert, deg_to_rad, twopi, pi
use shapelib, only: shpfileobject, shpclose
@@ -39,6 +39,7 @@ program extraction_eddies
real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
real step(2) ! longitude and latitude steps, in rad
logical periodic ! grid is periodic in longitude
+ real time
integer:: max_radius(2) = [20, 20]
! maximum radius of an eddy in longitude and latitude, in number of
@@ -67,6 +68,8 @@ program extraction_eddies
! latitude, in rad, of all the significant extrema, except the
! target extremum
+ character(len = 30) time_unit
+
namelist /main_nml/ min_amp, max_radius, min_radius
!--------------------------------------------------------------
@@ -81,6 +84,15 @@ program extraction_eddies
call nf95_open("h.nc", nf90_nowrite, ncid)
+ call find_coord(ncid, varid = varid, std_name = "time")
+ call nf95_get_att(ncid, varid, name = "units", values = time_unit)
+ call assert(time_unit == "days since 1950-01-01 00:00:00" &
+ .or. time_unit == "days since 1950-1-1 00:00:00", "main: bad time unit")
+ ! We are assuming there is a single date in the input file:
+ call nf95_get_var(ncid, varid, time)
+ k = nint(time)
+ call assert(abs(time - k) < 0.1, "main: bad time value")
+
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, lon_min)
@@ -105,21 +117,20 @@ program extraction_eddies
copy = merge(max_radius(1), 0, periodic)
allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), &
v(1 - copy:nlon + copy, nlat))
- k = 1
corner = [lon_min, lat_min] * deg_to_rad
- ! Read ssh, u and v at date k:
+ ! Read ssh, u and v:
- call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, name = "adt", &
+ call get_var(periodic, max_radius(1), ssh, ncid, nlon, name = "adt", &
new_fill_value = huge(0.))
! (We cannot keep the original fill value because Contour_531
! works with an upper limit for valid values.)
call nf95_close(ncid)
call nf95_open("uv.nc", nf90_nowrite, ncid)
- call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
+ call get_var(periodic, max_radius(1), u, ncid, nlon, name = "u", &
new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
- call get_var(periodic, max_radius(1), v, ncid, nlon, k, name = "v", &
+ call get_var(periodic, max_radius(1), v, ncid, nlon, name = "v", &
new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
! (We will need quiet NaNs rather the original fill values for u and
! v when we compute the max-speed contours and when we search the
@@ -172,7 +183,7 @@ program extraction_eddies
! Write snapshot:
do i = 1, s%number_vis_extr
- call write_eddy(s%list_vis(i), 1, i, hshp_extremum, hshp_outermost, &
+ call write_eddy(s%list_vis(i), k, i, hshp_extremum, hshp_outermost, &
hshp_max_speed)
end do
diff --git a/get_snapshot.f b/get_snapshot.f
index d8b6fc91..2b514572 100644
--- a/get_snapshot.f
+++ b/get_snapshot.f
@@ -36,6 +36,7 @@ contains
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
+ integer k2
!--------------------------------------------------------------
@@ -46,7 +47,7 @@ contains
hshp = hshp_outermost)
call shp_open_03("max_speed_contour", pszaccess = "rb", &
hshp = hshp_max_speed)
- call read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ call read_snapshot(s, k2, hshp_extremum, hshp_outermost, hshp_max_speed, &
corner, nlon, nlat)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
diff --git a/get_var.f b/get_var.f
index 67fe2382..6c5a5730 100644
--- a/get_var.f
+++ b/get_var.f
@@ -4,7 +4,7 @@ module get_var_m
contains
- subroutine get_var(periodic, max_rad_lon, values, ncid, nlon, k, name, &
+ subroutine get_var(periodic, max_rad_lon, values, ncid, nlon, name, &
new_fill_value)
! Read a NetCDF variable, change the missing value and extend it
@@ -25,7 +25,6 @@ contains
! periodic.
integer, intent(in):: ncid, nlon
- integer, intent(in):: k ! date index
character(len = *), intent(in):: name ! of NetCDF variable
real, intent(in):: new_fill_value
@@ -36,7 +35,7 @@ contains
!-------------------------------------------------------------------------
call nf95_inq_varid(ncid, name, varid)
- call nf95_get_var(ncid, varid, values(1:nlon, :), start = [1, 1, k])
+ call nf95_get_var(ncid, varid, values(1:nlon, :))
call nf95_get_att(ncid, varid, "_FillValue", Fill_Value)
! Change the missing value:
diff --git a/init_shapefiles.f b/init_shapefiles.f
index 4aae6472..c262ac46 100644
--- a/init_shapefiles.f
+++ b/init_shapefiles.f
@@ -32,7 +32,7 @@ contains
call dbf_add_field_03(ifield_extr_ssh, hshp_extremum, 'ssh', ftdouble, &
nwidth = 13, ndecimals = 6)
call dbf_add_field_03(ifield_extr_date, hshp_extremum, 'date_index', &
- ftinteger, nwidth = 4, ndecimals = 0)
+ ftinteger, nwidth = 5, ndecimals = 0)
call dbf_add_field_03(ifield_extr_eddy_index, hshp_extremum, 'eddy_index', &
ftinteger, nwidth = 5, ndecimals = 0)
call dbf_add_field_03(ifield_extr_interp, hshp_extremum, 'interpolat', &
@@ -50,7 +50,7 @@ contains
call dbf_add_field_03(ifield_out_ssh, hshp_outermost, 'ssh', ftdouble, &
nwidth = 13, ndecimals = 6)
call dbf_add_field_03(ifield_out_date, hshp_outermost, 'date_index', &
- ftinteger, nwidth = 4, ndecimals = 0)
+ ftinteger, nwidth = 5, ndecimals = 0)
call dbf_add_field_03(ifield_out_eddy_index, hshp_outermost, 'eddy_index', &
ftinteger, nwidth = 5, ndecimals = 0)
call dbf_add_field_03(ifield_out_radius4, hshp_outermost, 'radius4', &
@@ -62,7 +62,7 @@ contains
call dbf_add_field_03(ifield_max_speed_ssh, hshp_max_speed, 'ssh', &
ftdouble, nwidth = 13, ndecimals = 6)
call dbf_add_field_03(ifield_max_speed_date, hshp_max_speed, 'date_index', &
- ftinteger, nwidth = 4, ndecimals = 0)
+ ftinteger, nwidth = 5, ndecimals = 0)
call dbf_add_field_03(ifield_max_speed_eddy_index, hshp_max_speed, &
'eddy_index', ftinteger, nwidth = 5, ndecimals = 0)
--
GitLab