diff --git a/FindNetCDF_Fortran.cmake b/FindNetCDF_Fortran.cmake
index b13c608869c59def407737a79955386f23ab8ec0..10ee8c8169a8eaddd160d432bbcce9b0b9c7e2ec 100644
--- a/FindNetCDF_Fortran.cmake
+++ b/FindNetCDF_Fortran.cmake
@@ -5,22 +5,25 @@ if(TARGET NetCDF_Fortran::netcdff)
else()
# Find NetCDF dependency:
- option(use_find_netcdf_module "Use the find module for NetCDF")
-
- unset(extraArgs)
+ option(FIND_PACKAGE_PREFER_MODULE_netCDF
+ "Use directly the find module for NetCDF")
if(${CMAKE_FIND_PACKAGE_NAME}_FIND_QUIETLY)
- list(APPEND extraArgs QUIET)
+ set(maybe_quiet QUIET)
endif()
if(${CMAKE_FIND_PACKAGE_NAME}_FIND_REQUIRED)
- list(APPEND extraArgs REQUIRED)
+ set(maybe_required REQUIRED)
endif()
- if(use_find_netcdf_module)
- find_package(netCDF ${extraArgs})
+ if(FIND_PACKAGE_PREFER_MODULE_netCDF)
+ find_package(netCDF ${maybe_quiet} ${maybe_required})
else()
- find_package(netCDF CONFIG ${extraArgs})
+ find_package(netCDF CONFIG ${maybe_quiet})
+
+ if(NOT netCDF_FOUND)
+ find_package(netCDF ${maybe_quiet} ${maybe_required})
+ endif()
endif()
#-
diff --git a/Inst_eddies/Analysis/eddy_dump.py b/Inst_eddies/Analysis/eddy_dump.py
index 46868e22a05b9e7ce30e7f28e86f227eddd1247c..1eb2c95914c7cffedfb9b13737f9dd65c13ff4aa 100755
--- a/Inst_eddies/Analysis/eddy_dump.py
+++ b/Inst_eddies/Analysis/eddy_dump.py
@@ -15,25 +15,31 @@ import util_eddies
parser = argparse.ArgumentParser()
parser.add_argument("directory", help = "containing the three shapefiles")
+parser.add_argument("-i", "--ishape", type = int)
args = parser.parse_args()
-reply = input("days_1950, eddy_index = ? ").split(",")
-days_1950 = int(reply[0])
-eddy_index = int(reply[1])
-
handler = util_eddies.open_shpc(args.directory)
-# Find ishape:
-if handler["ishape_last"] is None:
- print("ishape_last not found. Iterating through shapes...")
-
- for ishape, rec in enumerate(handler["readers"]["extremum"].iterRecords()):
- if rec["days_1950"] == days_1950 and rec["eddy_index"] == eddy_index:
- break
- else:
- sys.exit("Eddy with this date and index not found")
+if args.ishape:
+ ishape = args.ishape
else:
- ishape = util_eddies.comp_ishape(handler, days_1950, eddy_index)
+ reply = input("days_1950, eddy_index = ? ").split(",")
+ days_1950 = int(reply[0])
+ eddy_index = int(reply[1])
+
+ # Find ishape:
+ if handler["ishape_last"] is None:
+ print("ishape_last not found. Iterating through shapes...")
+
+ for ishape, rec \
+ in enumerate(handler["readers"]["extremum"].iterRecords()):
+ if rec["days_1950"] == days_1950 \
+ and rec["eddy_index"] == eddy_index:
+ break
+ else:
+ sys.exit("Eddy with this date and index not found")
+ else:
+ ishape = util_eddies.comp_ishape(handler, days_1950, eddy_index)
for layer in ["extremum", "outermost_contour", "max_speed_contour", "centroid"]:
if layer in handler["readers"]:
diff --git a/Inst_eddies/Analysis/plot_eddy_contours.py b/Inst_eddies/Analysis/plot_eddy_contours.py
index 3ea6cef02ed43306d4ee30d85f99a3ae8997c789..13c2fbe14839e051a86408046737ef7d8fe162ac 100755
--- a/Inst_eddies/Analysis/plot_eddy_contours.py
+++ b/Inst_eddies/Analysis/plot_eddy_contours.py
@@ -131,15 +131,10 @@ def snapshot(ax, ishape_list, handler, *, dashed = False, light = False,
if not light: lines[0].set_marker("o")
if dashed: lines[0].set_linestyle("dashed")
-def plot_grid_bb(shpc_dir, ax):
+def plot_grid_bb(grid_nml, ax):
"""Grid bounding box."""
-
- file = path.join(shpc_dir, "grid_nml.txt")
- try:
- grid_nml = f90nml.read(file)["grid_nml"]
- except FileNotFoundError:
- print("grid_nml.txt not found. Will not plot bounding box.")
- else:
+
+ if grid_nml is not None:
step = grid_nml["STEP_DEG"]
rect = patches.Rectangle(grid_nml["corner_deg"],
(grid_nml["nlon"] - 1) * step[0],
@@ -153,10 +148,6 @@ if __name__ == "__main__":
import netCDF4
parser = argparse.ArgumentParser()
- parser.add_argument("-v", "--velocity", help = "plot velocity field",
- action = "store_true")
- parser.add_argument("-s", "--scale", default = 20, type = float,
- help = "scale of arrows for the velocity field")
parser.add_argument("-d", "--date", type = int,
help = "date, as days since 1950-1-1")
parser.add_argument("-g", "--grid", help = "plot grid",
@@ -175,17 +166,24 @@ if __name__ == "__main__":
help = "Save file to specified format")
args = parser.parse_args()
- if args.grid or args.velocity:
- with netCDF4.Dataset("h.nc") as f:
- if "lon" in f.variables:
- lon = "lon"
- lat = "lat"
- else:
- lon = "longitude"
- lat = "latitude"
-
- longitude = f[lon][:]
- latitude = f[lat][:]
+ if args.grid or args.window is None:
+ file = path.join(args.shpc_dir[0], "grid_nml.txt")
+
+ try:
+ grid_nml = f90nml.read(file)["grid_nml"]
+ except FileNotFoundError:
+ print("grid_nml.txt not found. Will not plot bounding box.")
+ grid_nml = None
+
+ if args.grid:
+ width = (grid_nml["NLON"] - 1) * grid_nml["STEP_DEG"][0]
+ height = (grid_nml["NLat"] - 1) * grid_nml["STEP_DEG"][1]
+ longitude = np.linspace(grid_nml["CORNER_DEG"][0],
+ grid_nml["CORNER_DEG"][0] + width,
+ grid_nml["NLON"])
+ latitude = np.linspace(grid_nml["CORNER_DEG"][1],
+ grid_nml["CORNER_DEG"][1] + height,
+ grid_nml["NLat"])
if args.window is not None:
llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = args.window
@@ -193,14 +191,14 @@ if __name__ == "__main__":
if urcrnrlon - llcrnrlon > 360:
sys.exit("bad values of urcrnrlon and llcrnrlon")
- if args.grid or args.velocity:
+ if args.grid:
longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
# (in [llcrnrlon, llcrnrlon + 2 pi[)
lon_mask = longitude <= urcrnrlon
lat_mask = np.logical_and(latitude >= llcrnrlat,
latitude <= urcrnrlat)
- elif args.grid or args.velocity:
+ elif args.grid:
lon_mask = np.ones(len(longitude), dtype = bool)
lat_mask = np.ones(len(latitude), dtype = bool)
@@ -216,19 +214,7 @@ if __name__ == "__main__":
ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
marker = "+", color = "gray", linestyle = "None")
- if args.window is None: plot_grid_bb(args.shpc_dir[0], ax)
-
- if args.velocity:
- with netCDF4.Dataset("uv.nc") as f:
- quiver_return = ax.quiver(longitude[lon_mask],
- latitude[lat_mask],
- f["ugos"][0, lat_mask][:, lon_mask],
- f["vgos"][0, lat_mask][:, lon_mask],
- scale = args.scale,
- scale_units = "width",
- transform = src_crs)
- plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
- coordinates = "figure")
+ if args.window is None: plot_grid_bb(grid_nml, ax)
for shpc_dir in args.shpc_dir:
handler = util_eddies.open_shpc(shpc_dir)
diff --git a/Inst_eddies/Analysis/plot_velocity.py b/Inst_eddies/Analysis/plot_velocity.py
new file mode 100755
index 0000000000000000000000000000000000000000..f0eb401b6624a32181392391d7e8beed21716693
--- /dev/null
+++ b/Inst_eddies/Analysis/plot_velocity.py
@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import cartopy.crs as ccrs
+import sys
+import matplotlib.pyplot as plt
+import argparse
+import netCDF4
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-s", "--scale", default = 20, type = float,
+ help = "scale of arrows for the velocity field")
+parser.add_argument("-w", "--window", help = "choose a limited plot window",
+ type = float, nargs = 4,
+ metavar = ("llcrnrlon", "llcrnrlat", "urcrnrlon",
+ "urcrnrlat"))
+parser.add_argument("--save", metavar = "format",
+ help = "Save file to specified format")
+parser.add_argument("-u", "--undefined", action = "store_true",
+ help = "plot points where velocity is not defined")
+parser.add_argument("input_file", help = "NetCDF file containing velocity")
+args = parser.parse_args()
+
+with netCDF4.Dataset(args.input_file) as f:
+ if "lon" in f.variables:
+ lon = "lon"
+ lat = "lat"
+ else:
+ lon = "longitude"
+ lat = "latitude"
+
+ longitude = f[lon][:]
+ latitude = f[lat][:]
+
+if args.window is not None:
+ llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = args.window
+
+ if urcrnrlon - llcrnrlon > 360:
+ sys.exit("bad values of urcrnrlon and llcrnrlon")
+
+ longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
+ # (in [llcrnrlon, llcrnrlon + 2 pi[)
+
+ lon_mask = longitude <= urcrnrlon
+ lat_mask = np.logical_and(latitude >= llcrnrlat,
+ latitude <= urcrnrlat)
+else:
+ lon_mask = np.ones(len(longitude), dtype = bool)
+ lat_mask = np.ones(len(latitude), dtype = bool)
+
+fig = plt.figure()
+src_crs = ccrs.PlateCarree()
+projection = ccrs.PlateCarree()
+##projection = ccrs.NorthPolarStereo()
+ax = plt.axes(projection = projection)
+
+with netCDF4.Dataset(args.input_file) as f:
+ if args.undefined:
+ undef_velocity = np.logical_or(f["ugos"][0].mask, f["vgos"][0].mask)
+ lon_2d, lat_2d = np.meshgrid(longitude[lon_mask],
+ latitude[lat_mask])
+ ax.plot(lon_2d[undef_velocity].reshape(-1),
+ lat_2d[undef_velocity].reshape(-1), transform = src_crs,
+ marker = "*", color = "violet", linestyle = "None")
+ else:
+ quiver_return = ax.quiver(longitude[lon_mask],
+ latitude[lat_mask],
+ f["ugos"][0, lat_mask][:, lon_mask],
+ f["vgos"][0, lat_mask][:, lon_mask],
+ scale = args.scale,
+ scale_units = "width",
+ transform = src_crs)
+ plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
+ coordinates = "figure")
+
+ax.gridlines(draw_labels = True)
+ax.coastlines()
+
+if args.save:
+ plt.savefig(f"plot_velocity.{args.save}")
+ print(f'Created "plot_velocity.{args.save}".')
+else:
+ plt.show()
diff --git a/Inst_eddies/Analysis/tests.json b/Inst_eddies/Analysis/tests.json
index 91951186ee71f1c75e59e1d67c07ef8db99919ed..b1e8ef7bb2551ac5c3a2d0ab8577910c3d295a32 100644
--- a/Inst_eddies/Analysis/tests.json
+++ b/Inst_eddies/Analysis/tests.json
@@ -59,18 +59,12 @@
"description": "With window option."
},
{
- "title": "Plot_eddy_contours_vel",
+ "title": "Plot_velocity",
"command":
[
- "$src_dir/Inst_eddies/Analysis/plot_eddy_contours.py",
- "$tests_old_dir/Extraction_eddies_region_4/SHPC_anti",
- "--velocity", "--save=png"
- ],
- "description": "With velocity option.",
- "required":
- [
- ["$src_dir/Inst_eddies/Tests/Input/h_region_4.nc", "h.nc"],
- ["$src_dir/Inst_eddies/Tests/Input/uv_region_4.nc", "uv.nc"]
+ "$src_dir/Inst_eddies/Analysis/plot_velocity.py",
+ "$src_dir/Inst_eddies/Tests/Input/uv_region_4.nc",
+ "--save=png"
]
},
{
diff --git a/Inst_eddies/Tests/CMakeLists.txt b/Inst_eddies/Tests/CMakeLists.txt
index e6cbcc95ed53cea12770bd4648f9fae21fcac8af..98e9b65eb4b1de0ba652a92c2579513dabc563e4 100644
--- a/Inst_eddies/Tests/CMakeLists.txt
+++ b/Inst_eddies/Tests/CMakeLists.txt
@@ -100,9 +100,22 @@ add_executable(test_set_max_speed
${PROJECT_SOURCE_DIR}/Common/read_field_indices.F90
${PROJECT_SOURCE_DIR}/Common/read_eddy.f90
${PROJECT_SOURCE_DIR}/Common/shpc_create.f90
- ${PROJECT_SOURCE_DIR}/Common/write_eddy.f90)
+ ${PROJECT_SOURCE_DIR}/Common/write_eddy.f90 get_var.f90)
target_link_libraries(test_set_max_speed Geometry::geometry
Numer_Rec_95::numer_rec_95 NetCDF95::netcdf95 Shapelib_03::shapelib_03
Contour_531::contour_531 Jumble::jumble NR_util::nr_util
NetCDF_Fortran::netcdff gpc_f)
+
+foreach(my_target IN ITEMS test_get_1_outerm test_set_all_outerm
+ test_good_contour test_inside_4 test_mean_speed
+ test_max_speed_contour_ssh test_nearby_extr test_local_extrema
+ test_set_max_speed)
+
+ set_target_properties(${my_target} PROPERTIES Fortran_MODULE_DIRECTORY
+ ${CMAKE_CURRENT_BINARY_DIR}/${my_target}_modules)
+
+ target_include_directories(${my_target} PRIVATE
+ ${PROJECT_BINARY_DIR}/${my_target}_modules>)
+
+endforeach()
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.dbf b/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..b43c318c495d62b98feb7fce2bf58317b59ef65e
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.dbf differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.shp b/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.shp
new file mode 100644
index 0000000000000000000000000000000000000000..29f893d5a6185f7407890b5a0eabd27760ade83d
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.shp differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.shx b/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.shx
new file mode 100644
index 0000000000000000000000000000000000000000..56a03bc6e6c9eb4566fd10213ddca20c25d155d6
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/extremum.shx differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.dbf b/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..0b01d3563aafcb726abfede3022bba1bf413721c
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.dbf differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.shp b/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.shp
new file mode 100644
index 0000000000000000000000000000000000000000..e08e311a72141da8c2709642b36b47313d475bf2
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.shp differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.shx b/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.shx
new file mode 100644
index 0000000000000000000000000000000000000000..d2645682a6962fc3a459c6fd5204dd2874cc72cd
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/max_speed_contour.shx differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/orientation.txt b/Inst_eddies/Tests/Input/SHPC_degenerate/orientation.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0e187bc2016925068243c8bbc4ca6eb8cfc73275
--- /dev/null
+++ b/Inst_eddies/Tests/Input/SHPC_degenerate/orientation.txt
@@ -0,0 +1 @@
+anticyclones
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.dbf b/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..5f453ce9304d2838b28515529408b2163c7cb6cd
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.dbf differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.shp b/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.shp
new file mode 100644
index 0000000000000000000000000000000000000000..05284bb77a0052d2b8ce72fa24bf2bedb004d638
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.shp differ
diff --git a/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.shx b/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.shx
new file mode 100644
index 0000000000000000000000000000000000000000..ecfb70a9f8ce65b456159a48e23dbbc3acf5e23c
Binary files /dev/null and b/Inst_eddies/Tests/Input/SHPC_degenerate/outermost_contour.shx differ
diff --git a/Inst_eddies/Tests/Input/degenerated_SSH.nc b/Inst_eddies/Tests/Input/degenerated_SSH.nc
new file mode 100644
index 0000000000000000000000000000000000000000..eaa5d3f79d863b71da0d5c8fde2f86ae01245013
Binary files /dev/null and b/Inst_eddies/Tests/Input/degenerated_SSH.nc differ
diff --git a/Inst_eddies/Tests/short_tests.json b/Inst_eddies/Tests/short_tests.json
index 5d72490281bacfc1593b594ab014397cd60a821c..1ddb6c3fa94386b0d41f47993fee606138e711db 100644
--- a/Inst_eddies/Tests/short_tests.json
+++ b/Inst_eddies/Tests/short_tests.json
@@ -270,6 +270,27 @@
]
]
},
+ {
+ "required":
+ [
+ ["$src_dir/Inst_eddies/Tests/Input/degenerated_SSH.nc", "h.nc"],
+ ["$src_dir/Inst_eddies/Tests/Input/degenerated_SSH.nc", "uv.nc"],
+ [
+ "$src_dir/Inst_eddies/Tests/Input/empty_outside_points.csv",
+ "outside_points.csv"
+ ]
+ ],
+ "title": "Set_max_speed_degenerate",
+ "description": "The SSH value of the extremum is exactly equal to the SSH value of the point next to the extremum, eastward. This points is selected as the point with maximum azimuthal speed. Then good_contour returns a null polyline.",
+ "commands":
+ [
+ ["mkdir", "SHPC"],
+ [
+ "$build_dir/Inst_eddies/test_set_max_speed",
+ "$src_dir/Inst_eddies/Tests/Input/SHPC_degenerate"
+ ]
+ ]
+ },
{
"create_file": ["main_nml.txt", "&main_nml min_amp = 0./\n"],
"input": "20454\n",
diff --git a/Inst_eddies/Tests/test_set_max_speed.f90 b/Inst_eddies/Tests/test_set_max_speed.f90
index 6e57a972d7d4c98ee871a01d6f6bb5efbb96e782..0df3250d40f4007f6cd03f1e8bb1ad54594f06af 100644
--- a/Inst_eddies/Tests/test_set_max_speed.f90
+++ b/Inst_eddies/Tests/test_set_max_speed.f90
@@ -1,6 +1,8 @@
program test_set_max_speed
use, intrinsic:: ISO_FORTRAN_ENV
+ use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
+
! Libraries:
use jumble, only: new_unit, count_lines, get_command_arg_dyn
@@ -10,6 +12,7 @@ program test_set_max_speed
use nr_util, only: deg_to_rad
use derived_types, only: eddy, shpc
+ use get_var_m, only: get_var
use read_eddy_m, only: read_eddy
use set_max_speed_m, only: set_max_speed
use shpc_close_m, only: shpc_close
@@ -57,15 +60,16 @@ program test_set_max_speed
print *, "Reading from uv.nc..."
call nf95_open("uv.nc", nf90_nowrite, ncid)
-
allocate(u(nlon, nlat), v(nlon, nlat))
-
- call nf95_inq_varid(ncid, "ugos", varid)
- call nf95_get_var(ncid, varid, u)
-
- call nf95_inq_varid(ncid, "vgos", varid)
- call nf95_get_var(ncid, varid, v)
-
+ call get_var(periodic = .false., max_rad_lon = 0, values = u, ncid = ncid, &
+ nlon = nlon, name = "ugos", &
+ new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+ call get_var(periodic = .false., max_rad_lon = 0, values = v, ncid = ncid, &
+ nlon = nlon, name = "vgos", &
+ 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
+ ! ssh of max-speed contours.)
call nf95_close(ncid)
print *, "Reading from shapefiles..."
diff --git a/Inst_eddies/max_speed_contour_ssh.f90 b/Inst_eddies/max_speed_contour_ssh.f90
index 695c49fddc3e176db547f1e28bfc96691fe5ea29..8b70fcefe3de65fe7cea13a994cd70ed67d779b6 100644
--- a/Inst_eddies/max_speed_contour_ssh.f90
+++ b/Inst_eddies/max_speed_contour_ssh.f90
@@ -33,7 +33,7 @@ contains
!----------------------------------------------------------------------
- ! Assuming that ascending indices in u and v mean correspond to
+ ! Assuming that ascending indices in u and v correspond to
! ascending longitude and latitude:
v_azim(east, :) = v(ind_extr(1) + 1:ind_extr(1) + radius - 1, &
ind_extr(2))
@@ -50,7 +50,12 @@ contains
! outermost contour as max-speed contour.
max_speed_contour_ssh = missing_ssh
else
- l = maxloc(abs(sum(v_azim, dim = 1) / 4.), dim = 1)
+ if (radius >= 3) then
+ l = maxloc(abs(sum(v_azim, dim = 1) / 4.), dim = 1)
+ else
+ l = 1
+ end if
+
direction = maxloc(abs(v_azim(:, l)), dim = 1)
select case (direction)
diff --git a/Inst_eddies/mean_speed.f90 b/Inst_eddies/mean_speed.f90
index 795e4040bf93f2761c48e3e5bd6f8d9c0c834597..6fb0c4411b5e37c995973d1f3eba9b171aa3089b 100644
--- a/Inst_eddies/mean_speed.f90
+++ b/Inst_eddies/mean_speed.f90
@@ -16,6 +16,7 @@ contains
! Libraries:
use contour_531, only: polyline
use numer_rec_95, only: bilinear_interp2_reg
+ use nr_util, only: assert
real, intent(in):: u(:, :), v(:, :) ! velocity
type(polyline), intent(in):: p ! should be closed
@@ -41,6 +42,8 @@ contains
!-------------------------------------------------------------------
+ call assert(p%n_points >= 2, &
+ "mean_speed: requires non-null, closed contour")
ni = p%n_points - 1
call bilinear_interp2_reg(corner, step, u, v, p%points(:, :ni), ui, vi)
diff --git a/Inst_eddies/set_max_speed.f90 b/Inst_eddies/set_max_speed.f90
index cea64de059a7286ccb9d6fd038a8d0faccb1f842..16d0776f87945bab75b2e413dbc5c4363ac5eb22 100644
--- a/Inst_eddies/set_max_speed.f90
+++ b/Inst_eddies/set_max_speed.f90
@@ -84,30 +84,41 @@ contains
e%speed_cont%ssh = max_speed_contour_ssh(ssh, u, v, ind_targ_extr, &
e%radius4)
-
+
if (e%speed_cont%ssh == missing_ssh) then
e%speed_cont = null_ssh_contour()
e%max_speed = speed_outerm
else
e%speed_cont%polyline = good_contour(corner, step, ssh, &
e%speed_cont%ssh, e%coord_extr, outside_points)
- e%max_speed = mean_speed(u, v, e%speed_cont%polyline, e%coord_extr, &
- corner, step)
- if (IEEE_IS_NAN(speed_outerm)) then
- ! This may happen when the outermost contour is near land.
- ! Stick to the contour coming from max_speed_contour_ssh.
- e%speed_cont%area = spher_polyline_area(e%speed_cont%polyline)
+ if (e%speed_cont%n_points == 0) then
+ ! This may happen if there is the exact same value of ssh
+ ! next to the extremum and this value is
+ ! e%speed_cont%ssh.
+ e%speed_cont%ssh = missing_ssh
+ e%speed_cont%area = - 1e6
+ e%max_speed = speed_outerm
else
- ! Note the following test should raise an invalid exception
- ! if e%max_speed is NaN.
- if (abs(speed_outerm) > abs(e%max_speed)) then
- ! Abandon the contour coming from max_speed_contour_ssh:
- e%speed_cont = null_ssh_contour()
- e%max_speed = speed_outerm
- else
+ e%max_speed = mean_speed(u, v, e%speed_cont%polyline, &
+ e%coord_extr, corner, step)
+
+ if (IEEE_IS_NAN(speed_outerm)) then
+ ! This may happen when the outermost contour is near land.
! Stick to the contour coming from max_speed_contour_ssh.
e%speed_cont%area = spher_polyline_area(e%speed_cont%polyline)
+ else
+ ! Note the following test should raise an invalid exception
+ ! if e%max_speed is NaN.
+ if (abs(speed_outerm) > abs(e%max_speed)) then
+ ! Abandon the contour coming from max_speed_contour_ssh:
+ e%speed_cont = null_ssh_contour()
+ e%max_speed = speed_outerm
+ else
+ ! Stick to the contour coming from max_speed_contour_ssh.
+ e%speed_cont%area &
+ = spher_polyline_area(e%speed_cont%polyline)
+ end if
end if
end if
end if
diff --git a/Overlap/eddy_graph.f90 b/Overlap/eddy_graph.f90
index 2697068cb223f701fb06b4803e59f1766ad10be0..028b375183b6d3bddd7d07b7d0c5aadd0e67fd70 100644
--- a/Overlap/eddy_graph.f90
+++ b/Overlap/eddy_graph.f90
@@ -93,7 +93,7 @@ program eddy_graph
copy = merge(dist_lim, 0, periodic)
call read_column(ishape_last, &
file = trim(shpc_dir) // "/ishape_last.txt", last = n_dates)
- n_dates = size(ishape_last)
+ if (n_dates == huge(0)) n_dates = size(ishape_last)
call assert(n_dates >= max_delta + 1, &
"eddy_graph: n_dates should be >= max_delta + 1")
call assert(n_proc <= n_dates / (max_delta + 1), &