diff --git a/GNUmakefile b/GNUmakefile
index 1efd84781a64a5d0d896bc357773c6ce800e4a5c..e29763e0fcd41ffd7affc9a1b6d88f107f1004df 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -11,9 +11,9 @@ src_test_get_1_outerm = good_contour.f test_get_1_outerm.f derived_types.f get_1
src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f good_contour.f max_speed_contour_ssh.f mean_speed.f spherical_polyline_area.f inside_4.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 set_max_speed.f outermost_possible_level.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spherical_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.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 set_max_speed.f outermost_possible_level.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spherical_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.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_polyline_area.f nearby_extr.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_polyline_area.f nearby_extr.f get_var.f
src_test_weight = test_weight.f weight.f derived_types.f
diff --git a/Tests/Input/Region_4/max_speed_contour_1.dbf b/Tests/Input/Region_4/max_speed_contour_1.dbf
index 358fe194c4f9640a020cb14479ef9faf9450ff4f..9ef3accdd7e36c10b06042935fd5d489cfd44354 100644
Binary files a/Tests/Input/Region_4/max_speed_contour_1.dbf and b/Tests/Input/Region_4/max_speed_contour_1.dbf differ
diff --git a/Tests/Input/Region_4/outermost_contour_1.dbf b/Tests/Input/Region_4/outermost_contour_1.dbf
index 28f9ea35bfdcecb073a814db71a54985cd2b2e2a..ed036b0bdab4a325fb5c8a2bd9bc4fcc61593a2f 100644
Binary files a/Tests/Input/Region_4/outermost_contour_1.dbf and b/Tests/Input/Region_4/outermost_contour_1.dbf differ
diff --git a/Tests/examine_ssh_values.py b/Tests/examine_ssh_values.py
index 5d9065fed0fa9334da29ba37d7f27a774ff14a76..b562ae6e7998d005889f0c6aee482b5e3f38bf28 100755
--- a/Tests/examine_ssh_values.py
+++ b/Tests/examine_ssh_values.py
@@ -4,8 +4,9 @@ import netCDF4
import matplotlib.pyplot as plt
import numpy as np
import math
+import sys
-with netCDF4.Dataset("h.nc") as f:
+with netCDF4.Dataset(sys.argv[1]) as f:
ssh = f.variables["adt"][:].squeeze()
longitude = f.variables["lon"][:]
latitude = f.variables["lat"][:]
@@ -16,7 +17,8 @@ lat_center = float(input("latitude of center = ? "))
# Base 0 indices, assuming regular grid::
i_center = round(float((lon_center - longitude[0])
/ (longitude[1] - longitude[0])))
-j_center = round(float((lat_center - latitude[0]) / (latitude[1] - latitude[0])))
+j_center = round(float((lat_center - latitude[0])
+ / (latitude[1] - latitude[0])))
# (float to convert from numpy.float64, so that round produces an int)
level = float(input("level = ? "))
@@ -32,6 +34,7 @@ plt.plot(longitude[i_center], latitude[j_center], "o", markersize = 10,
for i in np.arange(i_center - delta,i_center + delta + 1):
for j in np.arange(j_center - delta,j_center + delta + 1):
- plt.annotate(ssh[j,i], (longitude[i], latitude[j]))
+ if 0 <= i < longitude.size and 0 <= j < latitude.size:
+ plt.annotate(ssh[j,i], (longitude[i], latitude[j]))
plt.show()
diff --git a/Tests/read_eddy.f b/Tests/read_eddy.f
index 93c67d4aa1b2d3b95f1139667174a60cba332b53..8bacd743afdf89bf179833e1b19982d5acb1f525 100644
--- a/Tests/read_eddy.f
+++ b/Tests/read_eddy.f
@@ -12,7 +12,7 @@ contains
! Libraries:
use contour_531, only: null_polyline
use gpc_f, only: shpobj2pol, polygon
- use nr_util, only: deg_to_rad
+ use nr_util, only: deg_to_rad, pi
use shapelib, only: shpfileobject, shpobject, dbfreadattribute, &
shpdestroyobject
use shapelib_03, only: shp_read_object_03, dbf_get_field_index_03
@@ -20,8 +20,8 @@ contains
use derived_types, only: eddy
use read_field_indices_m, only: ifield_extr_ssh, ifield_extr_date, &
ifield_extr_eddy_index, ifield_extr_interp, ifield_extr_cycl, &
- ifield_extr_valid, ifield_extr_speed, ifield_out_area, &
- ifield_out_ssh, ifield_out_radius4, ifield_max_speed_area, &
+ ifield_extr_valid, ifield_extr_speed, ifield_out_r_eq_area, &
+ ifield_out_ssh, ifield_out_radius4, ifield_max_speed_r_eq_area, &
ifield_max_speed_ssh
type(eddy), intent(out):: e
@@ -75,8 +75,15 @@ contains
e%coord_extr = [psobject%padfx(1), psobject%padfy(1)] * deg_to_rad
call shpdestroyobject(psobject)
- call dbfreadattribute(hshp_outermost, ishape, ifield_out_area, attr = attr)
- e%out_cont%area = attr * 1e6 ! km2 to m2
+ call dbfreadattribute(hshp_outermost, ishape, ifield_out_r_eq_area, &
+ attr = attr)
+
+ if (attr >= 0._c_double) then
+ e%out_cont%area = pi * real(attr)**2 * 1e6 ! km2 to m2
+ else
+ e%out_cont%area = - 1e6
+ end if
+ ! (no merge statement because of the different real kinds)
call dbfreadattribute(hshp_outermost, ishape, ifield_out_ssh, attr = attr)
e%out_cont%ssh = attr
@@ -95,9 +102,15 @@ contains
e%out_cont%polyline%points = e%out_cont%polyline%points * deg_to_rad
end if
- call dbfreadattribute(hshp_max_speed, ishape, ifield_max_speed_area, &
- attr = attr)
- e%speed_cont%area = attr * 1e6 ! km2 to m2
+ call dbfreadattribute(hshp_max_speed, ishape, &
+ ifield_max_speed_r_eq_area, attr = attr)
+
+ if (attr >= 0._c_double) then
+ e%speed_cont%area = pi * real(attr)**2 * 1e6 ! km2 to m2
+ else
+ e%speed_cont%area = - 1e6
+ end if
+ ! (no merge statement because of the different real kinds)
call dbfreadattribute(hshp_max_speed, ishape, ifield_max_speed_ssh, &
attr = attr)
diff --git a/Tests/read_field_indices.f b/Tests/read_field_indices.f
index adb76f7353ede1ec8e52c886672abf72bb8940fe..ce4fe9e481555aee945ff4d4e207fcefdff76b8c 100644
--- a/Tests/read_field_indices.f
+++ b/Tests/read_field_indices.f
@@ -5,8 +5,9 @@ module read_field_indices_m
integer, protected:: ifield_extr_ssh, ifield_extr_date, &
ifield_extr_eddy_index, ifield_extr_interp, ifield_extr_cycl, &
ifield_extr_valid, ifield_extr_speed
- integer, protected:: ifield_out_area, ifield_out_ssh, ifield_out_radius4
- integer, protected:: ifield_max_speed_area, ifield_max_speed_ssh
+ integer, protected:: ifield_out_r_eq_area, ifield_out_ssh, &
+ ifield_out_radius4
+ integer, protected:: ifield_max_speed_r_eq_area, ifield_max_speed_ssh
contains
@@ -35,11 +36,13 @@ contains
call dbf_get_field_index_03(hshp_extremum, "valid", ifield_extr_valid)
call dbf_get_field_index_03(hshp_extremum, "speed", ifield_extr_speed)
- call dbf_get_field_index_03(hshp_outermost, "area", ifield_out_area)
+ call dbf_get_field_index_03(hshp_outermost, "r_eq_area", &
+ ifield_out_r_eq_area)
call dbf_get_field_index_03(hshp_outermost, "ssh", ifield_out_ssh)
call dbf_get_field_index_03(hshp_outermost, "radius4", ifield_out_radius4)
- call dbf_get_field_index_03(hshp_max_speed, "area", ifield_max_speed_area)
+ call dbf_get_field_index_03(hshp_max_speed, "r_eq_area", &
+ ifield_max_speed_r_eq_area)
call dbf_get_field_index_03(hshp_max_speed, "ssh", ifield_max_speed_ssh)
end subroutine read_field_indices
diff --git a/Tests/test_set_all_outerm.f b/Tests/test_set_all_outerm.f
index 4789fc254f122f27ab1fba33b43fa055601d31cf..4d42e010e47359811b8d8d06c6881c3bf81ec758 100644
--- a/Tests/test_set_all_outerm.f
+++ b/Tests/test_set_all_outerm.f
@@ -1,18 +1,21 @@
program test_set_all_outerm
+ use, intrinsic:: ISO_FORTRAN_ENV
+
! Libraries:
use contour_531, only: null_polyline
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, nf95_get_att
- use nr_util, only: pi, assert, deg_to_rad, rad_to_deg
+ 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 derived_types, only: snapshot
+ use get_var_m, only: get_var
use set_all_outerm_m, only: set_all_outerm
implicit none
@@ -24,37 +27,34 @@ program test_set_all_outerm
! in m
character(len = :), allocatable:: filename
- integer nlon, nlat
+ integer nlon, nlat, copy
real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
- integer, parameter:: max_radius(2) = [20, 20]
+ integer:: max_radius(2) = [20, 20]
! maximum radius of an eddy in longitude and latitude, in number of
! grid points
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
- ! order: ssh, date index, eddy index, interpolated, cyclone, suff
- ! amplitude.
-
- TYPE(shpfileobject) hshp_outermost
- ! shapefile outermost_contour_$m. The fields in the DBF file are,
- ! in that order: area, ssh, date index, eddy index.
+ real, allocatable:: ssh(:, :)
+ ! (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, unit_isolated, ishape
real step(2) ! longitude and latitude steps, in rad
+ logical periodic ! grid is periodic in longitude
+ namelist /main_nml/ min_amp, max_radius
!--------------------------------------------------------------
call get_command_arg_dyn(1, filename, &
"Required argument: NetCDF file containing adt")
- print *, "min_amp = ? "
- read *, min_amp
- print *, "min_amp = ", min_amp
+ 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)
! Read ssh:
@@ -70,23 +70,27 @@ program test_set_all_outerm
call nf95_get_var(ncid, varid, lat_min)
call nf95_get_var(ncid, varid, lat_max, start = [nlat])
- 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)
-
call assert(lon_max > lon_min .and. lat_max > lat_min, &
"test_set_all_outerm coordinate order")
step = [(lon_max - lon_min) / (nlon - 1), (lat_max - lat_min) / (nlat - 1)] &
* deg_to_rad
+
+ ! As we are requiring the grid spacing to be uniform, the value of
+ ! "periodic" must be consistent with the values of step(1) and nlon:
+ periodic = nint(twopi / step(1)) == nlon
+ print *, "periodic = ", periodic
call assert(2 * max_radius(1) * step(1) < pi, &
"test_set_all_outerm max_radius")
- call set_all_outerm(s, min_amp, max_radius, &
- corner = [lon_min, lat_min] * deg_to_rad, step = step, &
- periodic = .false., ssh = ssh)
+
+ 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", &
+ new_fill_value = huge(0.))
+
+ call nf95_close(ncid)
+
+ call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, &
+ corner = [lon_min, lat_min] * deg_to_rad)
call shp_create_03("extremum_1", shpt_point, hshp_extremum)
call dbf_add_field_03(ifield, hshp_extremum, 'ssh', ftdouble, nwidth = 13, &
@@ -103,8 +107,8 @@ program test_set_all_outerm
nwidth = 1, ndecimals = 0)
call shp_create_03("outermost_contour_1", shpt_polygon, hshp_outermost)
- call dbf_add_field_03(ifield, hshp_outermost, 'area', ftdouble, nwidth = 12, &
- ndecimals = 4)
+ 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, &
@@ -142,8 +146,13 @@ program test_set_all_outerm
end if
end if
- call dbf_write_attribute_03(hshp_outermost, ishape, 0, &
- s%list_vis(i)%out_cont%area / 1e6)
+ 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)
diff --git a/Tests/test_set_max_speed.f b/Tests/test_set_max_speed.f
index b8b587b09fd06d5a125ee32c604463047286e1a4..698a09f01010a1491ad0322789efd852c91916c9 100644
--- a/Tests/test_set_max_speed.f
+++ b/Tests/test_set_max_speed.f
@@ -8,7 +8,7 @@ program test_set_max_speed
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
find_coord, nf95_inquire_dimension
- use nr_util, only: deg_to_rad, rad_to_deg
+ use nr_util, only: deg_to_rad, rad_to_deg, pi
use shapelib, only: shpt_polygon, shpfileobject, ftdouble, shpclose, &
shpobject, dbfreadattribute
use shapelib_03, only: shp_open_03, shp_create_03, dbf_add_field_03, &
@@ -119,7 +119,7 @@ program test_set_max_speed
call shp_create_03("test_set_max_speed", shpt_polygon, hshp)
call dbf_add_field_03(ifield, hshp, 'ssh', ftdouble, nwidth = 13, &
ndecimals = 6)
- call dbf_add_field_03(ifield, hshp, 'area', ftdouble, nwidth = 12, &
+ call dbf_add_field_03(ifield, hshp, 'r_eq_area', ftdouble, nwidth = 10, &
ndecimals = 4)
call dbf_add_field_03(ifield, hshp, 'speed', ftdouble, nwidth = 13, &
ndecimals = 6)
@@ -127,7 +127,14 @@ program test_set_max_speed
call shp_append_object_03(ishape, hshp, shpt_polygon, &
e%speed_cont%points * rad_to_deg)
call dbf_write_attribute_03(hshp, ishape, 0, e%speed_cont%ssh)
- call dbf_write_attribute_03(hshp, ishape, 1, e%speed_cont%area / 1e6)
+
+ if (e%speed_cont%area >= 0) then
+ call dbf_write_attribute_03(hshp, ishape, 1, &
+ sqrt(e%speed_cont%area / 1e6 / pi))
+ else
+ call dbf_write_attribute_03(hshp, ishape, 1, - 100.)
+ end if
+
call dbf_write_attribute_03(hshp, ishape, 2, e%max_speed)
CALL shpclose(hshp)
diff --git a/Tests/tests.json b/Tests/tests.json
index 28450706f702394408873279e0f1e46a7d76e616..508399c2da2d93d82f3b721d18ba850de2977871 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -133,7 +133,7 @@
"args": ["$compil_prod_dir/test_set_all_outerm",
"$input_dir/h_region_3.nc"],
"title": "Set_all_outerm",
- "input": "0.001\n"
+ "input": "&main_nml min_amp = 0.001/\n"
},
{
"args": "$compil_prod_dir/test_get_snapshot",
@@ -359,5 +359,12 @@
"title" : "Local_extrema_periodic_2",
"input" : "t",
"description": "The data is really global so there is no discontinuity in the field and the longitude grid assuming periodicity is regular."
+ },
+ {
+ "args": ["$compil_prod_dir/test_set_all_outerm",
+ "$input_dir/h_2006_01_01_coarse.nc"],
+ "title": "Set_all_outerm_periodic",
+ "input": "&main_nml MAX_RADIUS = 4, 4/\n",
+ "description": "test_set_all_outerm with periodicity."
}
]
diff --git a/depend.mk b/depend.mk
index b2f70db5b6c804c51214cf26926a005f5723c809..f8a9e17951d32e48e03ec1e589d9d0f8d1a6b385 100644
--- a/depend.mk
+++ b/depend.mk
@@ -1,6 +1,6 @@
dispatch_snapshot.o : write_eddy.o send_snapshot.o derived_types.o
get_1_outerm.o : spherical_polyline_area.o good_contour.o derived_types.o
-get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o nearby_extr.o derived_types.o
+get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o nearby_extr.o get_var.o derived_types.o
nearby_extr.o : derived_types.o
read_eddy.o : read_field_indices.o derived_types.o
read_snapshot.o : read_field_indices.o read_eddy.o derived_types.o
@@ -16,7 +16,7 @@ test_local_extrema.o : local_extrema.o
test_nearby_extr.o : nearby_extr.o derived_types.o
test_read_eddy.o : write_eddy.o read_field_indices.o read_eddy.o init_shapefiles.o derived_types.o
test_read_snapshot.o : read_snapshot.o init_shapefiles.o dispatch_snapshot.o derived_types.o
-test_set_all_outerm.o : set_all_outerm.o derived_types.o
+test_set_all_outerm.o : set_all_outerm.o get_var.o derived_types.o
test_set_max_speed.o : set_max_speed.o derived_types.o
test_spherical_polygon_area.o : spherical_polygon_area.o
test_successive_overlap.o : successive_overlap.o read_snapshot.o derived_types.o
diff --git a/derived_types.f b/derived_types.f
index 8e1cc019f8aee57aa48e0cac74883597326a857a..91f54a58450eb40490777d88658fc372ea1baf2c 100644
--- a/derived_types.f
+++ b/derived_types.f
@@ -19,7 +19,8 @@ module derived_types
type(ssh_contour) speed_cont ! contour with maximum average azimuthal speed
real max_speed ! average of azimuthal speed on speed_cont, in m s-1
- logical valid ! valid (passing all requirements) out_cont found
+ logical valid ! valid out_cont found: not a null contour, and
+ ! with sufficient area
logical interpolated
diff --git a/get_snapshot.f b/get_snapshot.f
index 8f693f3b0109585cdeda17e49f9f51ad32283b6d..9ca95b5e9037953ff495a70d67000b60f301e01d 100644
--- a/get_snapshot.f
+++ b/get_snapshot.f
@@ -15,6 +15,7 @@ contains
use netcdf95, only: nf95_open, nf95_close
use derived_types, only: snapshot, null_ssh_contour, missing_speed
+ use get_var_m, only: get_var
use nearby_extr_m, only: nearby_extr
use receive_snapshot_m, only: receive_snapshot
use set_max_speed_m, only: set_max_speed
@@ -65,21 +66,22 @@ contains
! Read ssh, u and v at date k:
call nf95_open("h.nc", nf90_nowrite, ncid)
- call get_var(ssh, name = "adt", new_fill_value = huge(0.))
+ call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, 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(u, name = "u", &
+ call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
- call get_var(v, name = "v", &
+ call get_var(periodic, max_radius(1), v, ncid, nlon, k, 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.)
call nf95_close(ncid)
- call set_all_outerm(s, min_amp, max_radius, corner, step, periodic, ssh)
+ call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, corner)
! Done with outermost contours, now let us take care of the
! max-speed contours.
@@ -94,8 +96,9 @@ contains
urc = ceiling(convert_to_ind(maxval( &
s%list_vis(i)%out_cont%points, dim = 2), corner, step))
- if (.not. periodic) urc = min(urc, [nlon, nlat])
- ! (should have no effect except because of roundup error)
+ ! Should have no effect except because of roundup error:
+ urc(2) = min(urc(2), nlat)
+ if (.not. periodic) urc(1) = min(urc(1), nlon)
call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
@@ -115,42 +118,6 @@ contains
call receive_snapshot(s, source = m + 1, tag = k)
end if
- contains
-
- subroutine get_var(values, name, new_fill_value)
-
- ! Libraries:
- use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_get_att
-
- real, intent(out):: values(:, :)
- ! (1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1), 0,
- ! periodic), nlat) ssh, u or v
-
- character(len = *), intent(in):: name ! of NetCDF variable
- real, intent(in):: new_fill_value
-
- ! Local:
- integer varid
- real Fill_Value
-
- !-------------------------------------------------------------------------
-
- call nf95_inq_varid(ncid, name, varid)
- call nf95_get_var(ncid, varid, values(1:nlon, :), start = [1, 1, k])
- call nf95_get_att(ncid, varid, "_FillValue", Fill_Value)
-
- ! Change the missing value:
- where (values(1:nlon, :) == Fill_Value) values(1:nlon, :) = new_fill_value
-
- if (periodic) then
- ! Extend in longitude:
- values(1 - max_radius(1):0, :) &
- = values(nlon + 1 - max_radius(1):nlon, :)
- values(nlon + 1:nlon + max_radius(1), :) = values(1:max_radius(1), :)
- end if
-
- end subroutine get_var
-
end subroutine get_snapshot
end module get_snapshot_m
diff --git a/get_var.f b/get_var.f
new file mode 100644
index 0000000000000000000000000000000000000000..57452beb5611a399dd7b6f06e1bd38ef7c9c7b96
--- /dev/null
+++ b/get_var.f
@@ -0,0 +1,52 @@
+module get_var_m
+
+ implicit none
+
+contains
+
+ subroutine get_var(periodic, max_rad_lon, values, ncid, nlon, k, name, &
+ new_fill_value)
+
+ ! Read a NetCDF variable, change the missing value and extend it
+ ! in longitude.
+
+ ! Libraries:
+ use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_get_att
+
+ logical, intent(in):: periodic ! grid is periodic in longitude
+
+ integer, intent(in):: max_rad_lon ! maximum radius of an eddy in
+ ! longitude, in number of grid points
+
+ real, intent(out):: values(1 - merge(max_rad_lon, 0, periodic):, :)
+ ! (1 - merge(max_rad_lon, 0, periodic):nlon + merge(max_rad_lon, 0,
+ ! periodic), nlat) ssh, u or v
+
+ 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
+
+ ! Local:
+ integer varid
+ real Fill_Value
+
+ !-------------------------------------------------------------------------
+
+ call nf95_inq_varid(ncid, name, varid)
+ call nf95_get_var(ncid, varid, values(1:nlon, :), start = [1, 1, k])
+ call nf95_get_att(ncid, varid, "_FillValue", Fill_Value)
+
+ ! Change the missing value:
+ where (values(1:nlon, :) == Fill_Value) values(1:nlon, :) = new_fill_value
+
+ if (periodic) then
+ ! Extend in longitude:
+ values(1 - max_rad_lon:0, :) &
+ = values(nlon + 1 - max_rad_lon:nlon, :)
+ values(nlon + 1:nlon + max_rad_lon, :) = values(1:max_rad_lon, :)
+ end if
+
+ end subroutine get_var
+
+end module get_var_m
diff --git a/init_shapefiles.f b/init_shapefiles.f
index 24d6cb5bdf0f7a093fcaff16c585e34d77d164ae..30a0ef63daa022f66802d5d114fbc06eecfbec16 100644
--- a/init_shapefiles.f
+++ b/init_shapefiles.f
@@ -5,9 +5,9 @@ module init_shapefiles_m
integer, protected:: ifield_extr_ssh, ifield_extr_date, &
ifield_extr_eddy_index, ifield_extr_interp, ifield_extr_cycl, &
ifield_extr_valid, ifield_extr_speed
- integer, protected:: ifield_out_area, ifield_out_ssh, ifield_out_date, &
- ifield_out_eddy_index, ifield_out_radius4
- integer, protected:: ifield_max_speed_area, ifield_max_speed_ssh, &
+ integer, protected:: ifield_out_r_eq_area, ifield_out_ssh, &
+ ifield_out_date, ifield_out_eddy_index, ifield_out_radius4
+ integer, protected:: ifield_max_speed_r_eq_area, ifield_max_speed_ssh, &
ifield_max_speed_date, ifield_max_speed_eddy_index
contains
@@ -45,8 +45,8 @@ contains
nwidth = 13, ndecimals = 6)
call shp_create_03("outermost_contour_1", shpt_polygon, hshp_outermost)
- call dbf_add_field_03(ifield_out_area, hshp_outermost, 'area', ftdouble, &
- nwidth = 12, ndecimals = 4)
+ call dbf_add_field_03(ifield_out_r_eq_area, hshp_outermost, &
+ 'r_eq_area', ftdouble, nwidth = 10, ndecimals = 4)
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', &
@@ -57,8 +57,8 @@ contains
ftinteger, nwidth = 2, ndecimals = 0)
call shp_create_03("max_speed_contour_1", shpt_polygon, hshp_max_speed)
- call dbf_add_field_03(ifield_max_speed_area, hshp_max_speed, 'area', &
- ftdouble, nwidth = 12, ndecimals = 4)
+ call dbf_add_field_03(ifield_max_speed_r_eq_area, hshp_max_speed, &
+ 'r_eq_area', ftdouble, nwidth = 10, ndecimals = 4)
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', &
diff --git a/set_all_outerm.f b/set_all_outerm.f
index 893d14d23bb45f7616b8a94c2b8d7fe080a4857c..d8c29307c0424a7434e3321dee5287ec1c07f86f 100644
--- a/set_all_outerm.f
+++ b/set_all_outerm.f
@@ -4,7 +4,7 @@ module set_all_outerm_m
contains
- subroutine set_all_outerm(s, min_amp, max_radius, corner, step, periodic, ssh)
+ subroutine set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, corner)
! Extraction of eddies: find extrema and set all outermost
! contours in the snapshot "s". Not a function because "s" is not
@@ -30,16 +30,16 @@ contains
integer, intent(in):: max_radius(:) ! (2) maximum radius of an
! eddy in longitude and latitude, in number of grid points
- real, intent(in):: corner(:) ! (2) longitude and latitude of the
- ! corner of the global grid, in rad
-
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
logical, intent(in):: periodic ! grid is periodic in longitude
- real, intent(in):: ssh(:, :)
+ real, intent(in):: ssh(1 - merge(max_radius(1), 0, periodic):, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
! in longitude, else (nlon, nlat). Sea-surface height, in m.
+ real, intent(in):: corner(:) ! (2) longitude and latitude of the
+ ! corner of the global grid, in rad
+
! Local:
integer nlon, nlat
@@ -125,9 +125,12 @@ contains
llc = s%ind_extr(:, i) - max_radius
urc = s%ind_extr(:, i) + max_radius
+ llc(2) = max(llc(2), 1)
+ urc(2) = min(urc(2), nlat)
+
if (.not. periodic) then
- llc = max(llc, 1)
- urc = min(urc, [nlon, nlat])
+ llc(1) = max(llc(1), 1)
+ urc(1) = min(urc(1), nlon)
end if
corner_window = corner + (llc - 1) * step
diff --git a/write_eddy.f b/write_eddy.f
index feacd9ae2753571e82c230be7b2629d9e3df0421..bdd95b0c7a7e4c10e7adf7f6b898371a0e4b8e35 100644
--- a/write_eddy.f
+++ b/write_eddy.f
@@ -11,10 +11,11 @@ contains
use derived_types, only: eddy, missing_speed
use init_shapefiles_m, only: ifield_extr_ssh, ifield_extr_date, &
ifield_extr_eddy_index, ifield_extr_interp, ifield_extr_cycl, &
- ifield_extr_valid, ifield_extr_speed, ifield_out_area, &
+ ifield_extr_valid, ifield_extr_speed, ifield_out_r_eq_area, &
ifield_out_ssh, ifield_out_date, ifield_out_eddy_index, &
- ifield_out_radius4, ifield_max_speed_area, ifield_max_speed_ssh, &
- ifield_max_speed_date, ifield_max_speed_eddy_index
+ ifield_out_radius4, ifield_max_speed_r_eq_area, &
+ ifield_max_speed_ssh, ifield_max_speed_date, &
+ ifield_max_speed_eddy_index
use nr_util, only: pi
use shapelib, only: shpfileobject, shpt_polygon
use shapelib_03, only: shp_append_point_03, dbf_write_attribute_03, &
@@ -80,8 +81,14 @@ contains
end if
end if
- call dbf_write_attribute_03(hshp_outermost, ishape, ifield_out_area, &
- e%out_cont%area / 1e6)
+ if (e%out_cont%area >= 0) then
+ call dbf_write_attribute_03(hshp_outermost, ishape, &
+ ifield_out_r_eq_area, sqrt(e%out_cont%area / 1e6 / pi))
+ else
+ call dbf_write_attribute_03(hshp_outermost, ishape, &
+ ifield_out_r_eq_area, - 100.)
+ end if
+
call dbf_write_attribute_03(hshp_outermost, ishape, ifield_out_ssh, &
e%out_cont%ssh)
call dbf_write_attribute_03(hshp_outermost, ishape, ifield_out_date, k)
@@ -90,8 +97,14 @@ contains
call dbf_write_attribute_03(hshp_outermost, ishape, ifield_out_radius4, &
e%radius4)
- call dbf_write_attribute_03(hshp_max_speed, ishape, ifield_max_speed_area, &
- e%speed_cont%area / 1e6)
+ if (e%speed_cont%area >= 0) then
+ call dbf_write_attribute_03(hshp_max_speed, ishape, &
+ ifield_max_speed_r_eq_area, sqrt(e%speed_cont%area / 1e6 / pi))
+ else
+ call dbf_write_attribute_03(hshp_max_speed, ishape, &
+ ifield_max_speed_r_eq_area, - 100.)
+ end if
+
call dbf_write_attribute_03(hshp_max_speed, ishape, ifield_max_speed_ssh, &
e%speed_cont%ssh)
call dbf_write_attribute_03(hshp_max_speed, ishape, ifield_max_speed_date, &