diff --git a/Inst_eddies/CMakeLists.txt b/Inst_eddies/CMakeLists.txt
index e01c9920a5bf1ec94438c401f17139db9c049ae6..0bf64ba27076980ae067e14cc9c1d0bbab8f3492 100644
--- a/Inst_eddies/CMakeLists.txt
+++ b/Inst_eddies/CMakeLists.txt
@@ -2,7 +2,7 @@
add_executable(inst_eddies inst_eddies.f90 local_extrema.f90 set_max_speed.f90
get_1_outerm.f90 max_speed_contour_ssh.f90 good_contour.f90 mean_speed.f90
inside_4.f90 set_all_outerm.f90 nearby_extr.f90 get_var.f90 write_aux.f90
- config.f90)
+ config.f90 input_ssh.f90)
target_link_libraries(inst_eddies PRIVATE Contour_531::contour_531
Geometry::geometry NetCDF95::netcdf95 Shapelib_03::shapelib_03
Numer_Rec_95::numer_rec_95 Jumble::jumble NetCDF_Fortran::netcdff)
@@ -17,10 +17,11 @@ file(COPY inst_eddies_Aviso.py DESTINATION . FILE_PERMISSIONS
# Tests:
add_subdirectory(Tests)
-target_sources(test_get_1_outerm PRIVATE get_1_outerm.f90 good_contour.f90)
+target_sources(test_get_1_outerm PRIVATE get_1_outerm.f90 good_contour.f90
+ config.f90 input_ssh.f90 get_var.f90)
target_sources(test_set_all_outerm PRIVATE set_all_outerm.f90
local_extrema.f90 get_1_outerm.f90 good_contour.f90 nearby_extr.f90
- get_var.f90 config.f90)
+ get_var.f90 config.f90 input_ssh.f90)
target_sources(test_good_contour PRIVATE good_contour.f90)
target_sources(test_inside_4 PRIVATE inside_4.f90)
target_sources(test_mean_speed PRIVATE mean_speed.f90)
@@ -29,4 +30,5 @@ target_sources(test_max_speed_contour_ssh PRIVATE max_speed_contour_ssh.f90
target_sources(test_nearby_extr PRIVATE nearby_extr.f90)
target_sources(test_local_extrema PRIVATE local_extrema.f90)
target_sources(test_set_max_speed PRIVATE set_max_speed.f90 good_contour.f90
- max_speed_contour_ssh.f90 mean_speed.f90 inside_4.f90 get_var.f90)
+ max_speed_contour_ssh.f90 mean_speed.f90 inside_4.f90 get_var.f90
+ config.f90 input_ssh.f90)
diff --git a/Inst_eddies/Tests/long_tests.json b/Inst_eddies/Tests/long_tests.json
index 27b2aff6e449d90f7f1b5d87035e622dd9f958a2..c24c09bdd8bbda2571547af026c69933a77986d0 100644
--- a/Inst_eddies/Tests/long_tests.json
+++ b/Inst_eddies/Tests/long_tests.json
@@ -3,13 +3,15 @@
"commands":
[
["mkdir", "SHPC_cyclo", "SHPC_anti"],
- [
- "$build_dir/Inst_eddies/Tests/test_set_all_outerm",
- "$src_dir/Inst_eddies/Tests/Input/h_region_3.nc"
- ]
+ "$build_dir/Inst_eddies/Tests/test_set_all_outerm"
],
"title": "Set_all_outerm",
- "create_file": ["main_nml.txt", "&main_nml min_radius = 20./\n"]
+ "create_file": ["main_nml.txt", "&main_nml min_radius = 20./\n"],
+ "symlink":
+ [
+ ["$src_dir/Inst_eddies/Tests/Input/h_region_3.nc", "h.nc"],
+ ["$src_dir/Inst_eddies/Tests/Input/uv_region_3.nc", "uv.nc"]
+ ]
},
{
"commands":
diff --git a/Inst_eddies/Tests/short_tests.json b/Inst_eddies/Tests/short_tests.json
index 12833e5d2b60eb395b3d60c70d84ceb5af96707a..b694ce088bd02409957be37d403cc3841aa46f5a 100644
--- a/Inst_eddies/Tests/short_tests.json
+++ b/Inst_eddies/Tests/short_tests.json
@@ -1,12 +1,17 @@
[
{
"input": "&main_nml /\n20454\n",
+ "create_file": ["config_nml.txt", "&config_nml /\n"],
"title": "Get_1_outerm",
"symlink": [
[
"$src_dir/Inst_eddies/Tests/Input/Region_1/huv_region_1_2006_01_01.nc",
"h.nc"
],
+ [
+ "$src_dir/Inst_eddies/Tests/Input/Region_1/huv_region_1_2006_01_01.nc",
+ "uv.nc"
+ ],
[
"$src_dir/Inst_eddies/Tests/Input/Region_1/outside_points_get_1_outerm.csv",
"outside_points.csv"
@@ -19,12 +24,17 @@
"description":
"Assume insufficient amplitude for extrema 2 and 8. Even if extremum 8 does not have a particularly small amplitude, we just want to see the target extremum, 6, grow its outermost contour compared to the case without minimum amplitude.",
"input": "&main_nml /\n20454\n",
+ "create_file": ["config_nml.txt", "&config_nml /\n"],
"title": "Get_1_outerm_noise_2_8",
"symlink": [
[
"$src_dir/Inst_eddies/Tests/Input/Region_1/huv_region_1_2006_01_01.nc",
"h.nc"
],
+ [
+ "$src_dir/Inst_eddies/Tests/Input/Region_1/huv_region_1_2006_01_01.nc",
+ "uv.nc"
+ ],
"$src_dir/Inst_eddies/Tests/Input/Region_1/outside_points.csv"
],
"commands":
@@ -33,12 +43,17 @@
{
"description": "Negative value for extremum 2.",
"input": "&main_nml /\n20454\n",
+ "create_file": ["config_nml.txt", "&config_nml /\n"],
"title": "Get_1_outerm_noise_2",
"symlink": [
[
"$src_dir/Inst_eddies/Tests/Input/Region_1/huv_region_1_2006_01_01.nc",
"h.nc"
],
+ [
+ "$src_dir/Inst_eddies/Tests/Input/Region_1/huv_region_1_2006_01_01.nc",
+ "uv.nc"
+ ],
[
"$src_dir/Inst_eddies/Tests/Input/Region_1/outside_points_get_1_outerm_noise_2.csv",
"outside_points.csv"
@@ -222,6 +237,7 @@
]
},
{
+ "create_file": ["config_nml.txt", "&config_nml /\n"],
"symlink":
[
[
@@ -248,6 +264,7 @@
]
},
{
+ "create_file": ["config_nml.txt", "&config_nml /\n"],
"symlink":
[
[
@@ -271,6 +288,7 @@
]
},
{
+ "create_file": ["config_nml.txt", "&config_nml /\n"],
"symlink":
[
["$src_dir/Inst_eddies/Tests/Input/degenerated_SSH.nc", "h.nc"],
@@ -402,9 +420,14 @@
"commands":
[
["mkdir", "SHPC_cyclo", "SHPC_anti"],
- [
- "$build_dir/Inst_eddies/Tests/test_set_all_outerm",
- "$src_dir/Inst_eddies/Tests/Input/h_2006_01_01_coarse.nc"
+ "$build_dir/Inst_eddies/Tests/test_set_all_outerm"
+ ],
+ "symlink":
+ [
+ ["$src_dir/Inst_eddies/Tests/Input/h_2006_01_01_coarse.nc", "h.nc"],
+ [
+ "$src_dir/Inst_eddies/Tests/Input/uv_2006_01_01_coarse.nc",
+ "uv.nc"
]
]
},
diff --git a/Inst_eddies/Tests/test_get_1_outerm.f90 b/Inst_eddies/Tests/test_get_1_outerm.f90
index e4f7c764a4827a737dbf8355fb3db5cc55027867..0e6f6026b92dae16eb26a8a1ab8ce9718884d251 100644
--- a/Inst_eddies/Tests/test_get_1_outerm.f90
+++ b/Inst_eddies/Tests/test_get_1_outerm.f90
@@ -4,25 +4,34 @@ program test_get_1_outerm
! Libraries:
use jumble, only: new_unit, count_lines
- use netcdf, only: nf90_nowrite
- use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
- nf95_gw_var
use jumble, only: deg_to_rad
+ use config_m, only: config
use derived_types, only: eddy, shpc_slice_handler, null_ssh_contour, &
missing_speed
use get_1_outerm_m, only: get_1_outerm
+ use input_ssh_m, only: input_ssh
use shpc_close_m, only: shpc_close
use shpc_create_m, only: shpc_create
use write_eddy_m, only: write_eddy
implicit none
- integer ncid, varid, unit, n, l, date
- real longitude(2) ! in rad
- real latitude(2) ! in rad
- real step(2) ! in rad
- real, allocatable:: ssh(:, :, :) ! (:, :, 1) sea-surface height, in m
+ integer unit, n, l, date, nlon, nlat
+
+ real corner_deg(2), corner(2)
+ ! longitude and latitude of the corner of the whole grid, in degrees
+ ! and in rad
+
+ real step_deg(2), step(2) ! longitude and latitude steps, in degrees and rad
+
+ 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.
+
+ real, allocatable:: u(:, :), v(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+ ! (nlon, nlat) wind, in m s-1
integer:: ind_targ_extr(2) = [19, 11]
! indices in the two dimensions of target extremum
@@ -30,6 +39,7 @@ program test_get_1_outerm
real:: innermost_level = 0.2933
! level of innermost contour, for each extremum
+ logical periodic ! grid is periodic in longitude
logical:: cyclone = .true.
TYPE(shpc_slice_handler) hshp
type(eddy) e
@@ -46,26 +56,11 @@ program test_get_1_outerm
write(unit = error_unit, fmt = *) "Enter namelist main_nml."
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
-
- print *, 'Reading from "h.nc"...'
- call nf95_open("h.nc", nf90_nowrite, ncid)
-
- call nf95_inq_varid(ncid, "lon", varid)
- call nf95_get_var(ncid, varid, longitude)
- longitude = longitude * deg_to_rad
-
- call nf95_inq_varid(ncid, "lat", varid)
- call nf95_get_var(ncid, varid, latitude)
- latitude = latitude * deg_to_rad
-
- call nf95_inq_varid(ncid, "adt", varid)
- call nf95_gw_var(ncid, varid, ssh)
-
- call nf95_close(ncid)
+ call config(exist = .false.)
+ call input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v, corner_deg, &
+ step_deg)
print *, "date = ?"
read *, date
- step = [longitude(2) - longitude(1), latitude(2) - latitude(1)]
-
print *, "Reading from outside_points.csv..."
call new_unit(unit)
open(unit, file = "outside_points.csv", status = "old", action = "read", &
@@ -83,15 +78,13 @@ program test_get_1_outerm
outside_points = outside_points * deg_to_rad
close(unit)
- e%coord_extr = [longitude(1), latitude(1)] + (ind_targ_extr - 1) * step
+ e%coord_extr = corner + (ind_targ_extr - 1) * step
- e%out_cont = get_1_outerm(cyclone, &
- coord_extr = e%coord_extr, &
- innermost_level = innermost_level, outside_points = outside_points, &
- ssh = ssh(:, :, 1), corner = [longitude(1), latitude(1)], step = step)
+ e%out_cont = get_1_outerm(cyclone, e%coord_extr, innermost_level, &
+ outside_points, ssh, corner, step)
if (e%out_cont%closed) then
- e%ssh_extr = ssh(ind_targ_extr(1), ind_targ_extr(2), 1)
+ e%ssh_extr = ssh(ind_targ_extr(1), ind_targ_extr(2))
e%cyclone = cyclone
e%speed_cont = null_ssh_contour()
e%max_speed = missing_speed
diff --git a/Inst_eddies/Tests/test_set_all_outerm.f90 b/Inst_eddies/Tests/test_set_all_outerm.f90
index d2319d02c7ecb859106922626de527728e35e23d..f208eba09bc7fc9a2ff43db03f5b8e07a5c7124b 100644
--- a/Inst_eddies/Tests/test_set_all_outerm.f90
+++ b/Inst_eddies/Tests/test_set_all_outerm.f90
@@ -3,16 +3,12 @@ program test_set_all_outerm
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
- use jumble, only: get_command_arg_dyn
- use netcdf, only: nf90_nowrite
- use netcdf95, only: nf95_open, nf95_close, nf95_get_var, find_coord, &
- nf95_inquire_dimension
- use jumble, only: pi, assert, deg_to_rad, twopi
+ use jumble, only: pi
- use config_m, only: config, max_radius_deg, min_radius
+ use config_m, only: config, min_radius
use derived_types, only: snapshot, shpc_slice_handler, null_ssh_contour, &
missing_speed
- use get_var_m, only: get_var
+ use input_ssh_m, only: input_ssh
use set_all_outerm_m, only: set_all_outerm
use shpc_close_m, only: shpc_close
use shpc_create_m, only: shpc_create
@@ -21,66 +17,31 @@ program test_set_all_outerm
implicit none
type(snapshot) s
- character(len = :), allocatable:: filename
- integer nlon, nlat, copy
- real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
+ integer nlon, nlat
+ integer i, n_cyclo, n_anti
- integer max_radius(2)
- ! maximum radius of an eddy in longitude and latitude, in number of
- ! grid points
-
- integer ncid, varid, dimid, i, n_cyclo, n_anti
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.
+ real, allocatable:: u(:, :), v(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+ ! (nlon, nlat) wind, in m s-1
+
+ real corner_deg(2), corner(2)
+ ! longitude and latitude of the corner of the whole grid, in degrees
+ ! and in rad
+
real step_deg(2), step(2) ! longitude and latitude steps, in degrees and rad
logical periodic ! grid is periodic in longitude
TYPE(shpc_slice_handler) hshp_cyclo, hshp_anti
!--------------------------------------------------------------
- call get_command_arg_dyn(1, filename, &
- "Required argument: NetCDF file containing adt")
call config(exist = .false.)
-
- ! Read ssh:
-
- call nf95_open(filename, nf90_nowrite, ncid)
-
- 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)
- call nf95_get_var(ncid, varid, lon_max, start = [nlon])
-
- call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
- call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
- call nf95_get_var(ncid, varid, lat_min)
- call nf95_get_var(ncid, varid, lat_max, start = [nlat])
-
- call assert(lon_max > lon_min .and. lat_max > lat_min, &
- "test_set_all_outerm coordinate order")
- step_deg = [(lon_max - lon_min) / (nlon - 1), &
- (lat_max - lat_min) / (nlat - 1)]
- step = step_deg * deg_to_rad
- max_radius = nint(max_radius_deg / step_deg)
-
- ! As we are requiring the grid spacing to be uniform, the value of
- ! "periodic" may be deduced from 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")
-
- copy = merge(max_radius(1), 0, periodic)
- allocate(ssh(1 - copy:nlon + copy, nlat))
- call get_var(periodic, max_radius(1), ssh, ncid, nlon, name = "adt", &
- new_fill_value = huge(0.))
-
- call nf95_close(ncid)
-
- call set_all_outerm(s, max_radius, step, periodic, ssh, &
- corner = [lon_min, lat_min] * deg_to_rad, &
+ call input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v, corner_deg, &
+ step_deg)
+ call set_all_outerm(s, step, periodic, ssh, corner, &
min_area = pi * (min_radius * 1e3)**2)
do i = 1, s%number_extr
diff --git a/Inst_eddies/Tests/test_set_max_speed.f90 b/Inst_eddies/Tests/test_set_max_speed.f90
index 9a214d49ac01d201a42baeb4d7e5fa03ec4ab49c..8c1733265aa6cf59e02c7ca165dcc56fc84748a6 100644
--- a/Inst_eddies/Tests/test_set_max_speed.f90
+++ b/Inst_eddies/Tests/test_set_max_speed.f90
@@ -1,17 +1,13 @@
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, deg_to_rad
- use netcdf, only: nf90_nowrite
- use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
- find_coord, nf95_inquire_dimension
+ use config_m, only: config
use derived_types, only: eddy, shpc_slice_handler
- use get_var_m, only: get_var
+ use input_ssh_m, only: input_ssh
use read_eddy_m, only: read_eddy
use set_max_speed_m, only: set_max_speed
use shpc_close_m, only: shpc_close
@@ -21,12 +17,17 @@ program test_set_max_speed
implicit none
- integer ncid, varid, dimid
+ logical periodic ! grid is periodic in longitude
integer nlon, nlat, unit, n, l, d, i
real, allocatable:: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
real, allocatable:: u(:, :), v(:, :) ! (nlon, nlat) wind, in m s-1
type(eddy) e
- real corner(2), step(2) ! in rad
+
+ real corner_deg(2), corner(2)
+ ! longitude and latitude of the corner of the whole grid, in degrees
+ ! and in rad
+
+ real step_deg(2), step(2) ! longitude and latitude steps, in degrees and rad
TYPE(shpc_slice_handler) hshp
character(len = :), allocatable:: shpc_dir
@@ -36,40 +37,9 @@ program test_set_max_speed
!----------------------------------------------------------------
call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")
- print *, "Reading from h.nc..."
- call nf95_open("h.nc", nf90_nowrite, ncid)
-
- ! Read corner:
-
- 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, corner(1))
-
- call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
- call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
- call nf95_get_var(ncid, varid, corner(2))
-
- corner = corner * deg_to_rad
-
- allocate(ssh(nlon, nlat))
- call nf95_inq_varid(ncid, "adt", varid)
- call nf95_get_var(ncid, varid, ssh)
-
- call nf95_close(ncid)
-
- print *, "Reading from uv.nc..."
- call nf95_open("uv.nc", nf90_nowrite, ncid)
- allocate(u(nlon, nlat), v(nlon, nlat))
- 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)
+ call config(exist = .true.)
+ call input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v, corner_deg, &
+ step_deg)
print *, "Reading from shapefiles..."
call shpc_open(hshp, shpc_dir, pszaccess = "rb")
diff --git a/Inst_eddies/input_ssh.f90 b/Inst_eddies/input_ssh.f90
new file mode 100644
index 0000000000000000000000000000000000000000..7c2234bc6ed8bc5f391ca00cff13da2c264403fc
--- /dev/null
+++ b/Inst_eddies/input_ssh.f90
@@ -0,0 +1,110 @@
+module input_ssh_m
+
+ implicit none
+
+ integer, protected:: max_radius(2)
+ ! maximum radius of an eddy in longitude and latitude, in number of
+ ! grid points
+
+contains
+
+ subroutine input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v, &
+ corner_deg, step_deg)
+
+ use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
+
+ ! Libraries:
+ use jumble, only: assert, deg_to_rad, twopi
+ use netcdf, only: nf90_nowrite
+ use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
+ nf95_get_var, nf95_close
+
+ use config_m, only: max_radius_deg
+ use get_var_m, only: get_var
+
+ real, intent(out):: corner(:) ! (2)
+ ! longitude and latitude of the corner of the whole grid, in rad
+
+ real, intent(out):: step(:) ! (2) longitude and latitude steps, in rad
+
+ integer, intent(out):: nlon, nlat
+ ! size of ssh array in input NetCDF, assuming no repeated point if
+ ! the grid is global
+
+ logical, intent(out):: periodic ! grid is periodic in longitude
+
+ real, allocatable, intent(out):: ssh(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+ ! (nlon, nlat) sea-surface height, in m
+
+ real, allocatable, intent(out):: u(:, :), v(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+ ! (nlon, nlat) wind, in m s-1
+
+ real, intent(out):: corner_deg(:) ! (2)
+ ! longitude and latitude of the corner of the whole grid, in degrees
+
+ real, intent(out):: step_deg(:) ! (2)
+ ! longitude and latitude steps, in degrees
+
+ ! Local:
+
+ integer ncid, dimid, varid
+ real lon_max, lat_max ! longitude and latitude, in degrees
+ integer copy
+
+ !----------------------------------------------------------------------
+
+ call nf95_open("h.nc", nf90_nowrite, ncid)
+
+ 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, corner_deg(1))
+ call nf95_get_var(ncid, varid, lon_max, start = [nlon])
+
+ call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
+ call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
+ call nf95_get_var(ncid, varid, corner_deg(2))
+ call nf95_get_var(ncid, varid, lat_max, start = [nlat])
+
+ call assert(lon_max > corner_deg(1) .and. lat_max > corner_deg(2), &
+ "inst_eddies coordinate order")
+ step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
+ (lat_max - corner_deg(2)) / (nlat - 1)]
+ step = step_deg * deg_to_rad
+
+ ! As we are requiring the grid spacing to be uniform, the value of
+ ! "periodic" may be deduced from the values of step(1) and nlon:
+ periodic = nint(twopi / step(1)) == nlon
+ print *, "inst_eddies: periodic = ", periodic
+
+ call assert(2 * max_radius_deg(1) < 180., "inst_eddies max_radius_deg")
+ ! (Let us require this even if not periodic. This is clearer.)
+
+ max_radius = nint(max_radius_deg / step_deg)
+ 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))
+ corner = corner_deg * deg_to_rad
+
+ ! Read ssh, u and v:
+
+ 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, name = "ugos", &
+ new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+ call get_var(periodic, max_radius(1), v, ncid, 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)
+
+ end subroutine input_ssh
+
+end module input_ssh_m
diff --git a/Inst_eddies/inst_eddies.f90 b/Inst_eddies/inst_eddies.f90
index 88d61dada2378a40361680fe5efdc8fcdec80f0c..751967650cf8084381235c6d7466a7f6d6a52f14 100644
--- a/Inst_eddies/inst_eddies.f90
+++ b/Inst_eddies/inst_eddies.f90
@@ -1,19 +1,16 @@
program inst_eddies
use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
- ieee_support_nan, ieee_value, IEEE_QUIET_NAN
+ ieee_support_nan
! Libraries:
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
- use jumble, only: assert, deg_to_rad, twopi, pi, new_unit
+ use jumble, only: assert, twopi, pi, new_unit
- use config_m, only: config, max_radius_deg, min_radius
+ use config_m, only: config, min_radius
use derived_types, only: snapshot, null_ssh_contour, missing_speed, &
shpc_slice_handler
- use get_var_m, only: get_var
+ use input_ssh_m, only: input_ssh
use nearby_extr_m, only: nearby_extr
use set_all_outerm_m, only: set_all_outerm
use set_max_speed_m, only: set_max_speed
@@ -28,21 +25,15 @@ program inst_eddies
type(snapshot) s
TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
integer i, date, n_cyclo, n_anti
-
- integer ncid, dimid, varid, copy, unit
+ integer unit
integer nlon, nlat
! size of ssh array in input NetCDF, assuming no repeated point if
! the grid is global
- real lon_max, lat_max ! longitude and latitude, in degrees
real step_deg(2), step(2) ! longitude and latitude steps, in degrees and rad
logical periodic ! grid is periodic in longitude
- integer max_radius(2)
- ! maximum radius of an eddy in longitude and latitude, in number of
- ! grid points
-
real, allocatable:: ssh(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
! (nlon, nlat) sea-surface height, in m
@@ -90,60 +81,12 @@ program inst_eddies
print *, "date = ?"
read *, date
- call nf95_open("h.nc", nf90_nowrite, ncid)
-
- 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, corner_deg(1))
- call nf95_get_var(ncid, varid, lon_max, start = [nlon])
-
- call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
- call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
- call nf95_get_var(ncid, varid, corner_deg(2))
- call nf95_get_var(ncid, varid, lat_max, start = [nlat])
-
- call assert(lon_max > corner_deg(1) .and. lat_max > corner_deg(2), &
- "inst_eddies coordinate order")
- step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
- (lat_max - corner_deg(2)) / (nlat - 1)]
- step = step_deg * deg_to_rad
-
- ! As we are requiring the grid spacing to be uniform, the value of
- ! "periodic" may be deduced from the values of step(1) and nlon:
- periodic = nint(twopi / step(1)) == nlon
- print *, "inst_eddies: periodic = ", periodic
-
- call assert(2 * max_radius_deg(1) < 180., "inst_eddies max_radius_deg")
- ! (Let us require this even if not periodic. This is clearer.)
-
- max_radius = nint(max_radius_deg / step_deg)
- 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))
- corner = corner_deg * deg_to_rad
-
- ! Read ssh, u and v:
-
- 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, name = "ugos", &
- new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
- call get_var(periodic, max_radius(1), v, ncid, 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)
-
+ call input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v, corner_deg, &
+ step_deg)
call cpu_time(t1)
write(unit, fmt = *) "CPU time for input, before computation:", t1 - t0, "s"
t0 = t1
- call set_all_outerm(s, max_radius, step, periodic, ssh, corner, &
+ call set_all_outerm(s, step, periodic, ssh, corner, &
min_area = pi * (min_radius * 1e3)**2)
! Done with outermost contours, now let us take care of the
diff --git a/Inst_eddies/set_all_outerm.f90 b/Inst_eddies/set_all_outerm.f90
index fcb89c89f051a29e003a2c9b03aacbd198b46b12..40b95e8de1caba35409be417ba7f0cac78e8cbd9 100644
--- a/Inst_eddies/set_all_outerm.f90
+++ b/Inst_eddies/set_all_outerm.f90
@@ -4,8 +4,7 @@ module set_all_outerm_m
contains
- subroutine set_all_outerm(s, max_radius, step, periodic, ssh, corner, &
- min_area)
+ subroutine set_all_outerm(s, step, periodic, ssh, corner, min_area)
! Extraction of eddies: find extrema and set all outermost
! contours in the snapshot "s". Not a function because "s" is not
@@ -19,6 +18,7 @@ contains
use config_m, only: min_amp
use derived_types, only: snapshot
use get_1_outerm_m, only: get_1_outerm
+ use input_ssh_m, only: max_radius
use nearby_extr_m, only: nearby_extr
use local_extrema_m, only: local_extrema
@@ -26,9 +26,6 @@ contains
! Specifically: define everything in s except
! s%list%speed_cont, s%list%max_speed and s%list%radius4.
- integer, intent(in):: max_radius(:) ! (2) maximum radius of an
- ! eddy in longitude and latitude, in number of grid points
-
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
logical, intent(in):: periodic ! grid is periodic in longitude