From 6c00adf7d85273ad83884e1c6e50d7d7d79992dd Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Fri, 2 Jun 2023 12:07:01 +0200
Subject: [PATCH] Create program `examine_eddy`
---
Common/CMakeLists.txt | 2 +
Inst_eddies/CMakeLists.txt | 4 +
Inst_eddies/Tests/CMakeLists.txt | 11 +-
Inst_eddies/Tests/examine_eddy.f90 | 186 +++++++++++++++++++++++++++++
4 files changed, 202 insertions(+), 1 deletion(-)
create mode 100644 Inst_eddies/Tests/examine_eddy.f90
diff --git a/Common/CMakeLists.txt b/Common/CMakeLists.txt
index 06e267a0..732cdcce 100644
--- a/Common/CMakeLists.txt
+++ b/Common/CMakeLists.txt
@@ -12,6 +12,8 @@ target_sources(inst_eddies PRIVATE spher_polyline_area.f90
get_slice_dir.f90 read_field_indices.f90)
target_sources(test_write_null PRIVATE derived_types.f90 shpc_open.f90
shpc_close.f90 read_field_indices.f90 get_slice_dir.f90)
+target_sources(examine_eddy PRIVATE derived_types.f90 get_slice_dir.f90
+ spher_polyline_area.f90 shpc_close.f90)
if(MPI_Fortran_HAVE_F08_MODULE)
target_sources(test_overlap PRIVATE derived_types.f90 spher_polyline_area.f90
diff --git a/Inst_eddies/CMakeLists.txt b/Inst_eddies/CMakeLists.txt
index cd13b6df..d6ca08b6 100644
--- a/Inst_eddies/CMakeLists.txt
+++ b/Inst_eddies/CMakeLists.txt
@@ -30,3 +30,7 @@ target_sources(test_set_max_speed PRIVATE set_max_speed.F90 good_contour.f90
ccw_orient.f90 write_eddy.f90 complete_ssh.f90)
target_sources(test_write_null PRIVATE shpc_create.f90)
target_sources(test_complete_ssh PRIVATE complete_ssh.f90)
+target_sources(examine_eddy PRIVATE config.f90 input_ssh.f90 get_var.f90
+ shpc_create.f90 set_all_extr.f90 local_extrema.f90 write_eddy.f90
+ set_contours.f90 get_1_outerm.f90 ccw_orient.f90 good_contour.f90
+ set_max_speed.F90 complete_ssh.f90 mean_speed.f90)
diff --git a/Inst_eddies/Tests/CMakeLists.txt b/Inst_eddies/Tests/CMakeLists.txt
index 2324f7b6..734af5e1 100644
--- a/Inst_eddies/Tests/CMakeLists.txt
+++ b/Inst_eddies/Tests/CMakeLists.txt
@@ -46,8 +46,17 @@ target_link_libraries(test_complete_ssh PRIVATE Jumble::jumble)
# inst_eddies
target_sources(inst_eddies PRIVATE cont_list.f90)
+# examine_eddies
+add_executable(examine_eddy examine_eddy.f90 cont_list.f90
+ ${PROJECT_SOURCE_DIR}/Overlap/read_eddy.f90)
+target_link_libraries(examine_eddy PRIVATE Contour_531::contour_531
+ Jumble::jumble Shapelib_03::shapelib_03 NetCDF95::netcdf95 gpc_f
+ Geometry::geometry Numer_Rec_95::numer_rec_95)
+target_compile_definitions(examine_eddy PRIVATE $<$<CONFIG:Debug>:DEBUG>)
+
foreach(my_target IN ITEMS test_get_1_outerm test_good_contour test_mean_speed
- test_nearby_extr test_local_extrema test_set_max_speed test_complete_ssh)
+ test_nearby_extr test_local_extrema test_set_max_speed test_complete_ssh
+ examine_eddy)
set_target_properties(${my_target} PROPERTIES Fortran_MODULE_DIRECTORY
${CMAKE_CURRENT_BINARY_DIR}/${my_target}_modules)
target_include_directories(${my_target} PRIVATE
diff --git a/Inst_eddies/Tests/examine_eddy.f90 b/Inst_eddies/Tests/examine_eddy.f90
new file mode 100644
index 00000000..f25fca9d
--- /dev/null
+++ b/Inst_eddies/Tests/examine_eddy.f90
@@ -0,0 +1,186 @@
+program examine_eddy
+
+ ! Libraries:
+ use contour_531, only: convert_to_ind
+ use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, pi, &
+ rad_to_deg, twopi
+ use shapelib, only: shpfileobject, shpclose
+ use shapelib_03, only: dbf_read_attribute_03, dbf_get_field_index_03, &
+ shp_read_point, shp_open_03
+
+ use config_m, only: config, min_radius
+ use cont_list_m, only: create_cont_list, close_cont_list
+ use derived_types, only: snapshot, shpc_slice_handler, eddy, missing_ssh
+ use input_ssh_m, only: input_ssh, max_radius
+ use read_eddy_m, only: read_eddy
+ use shpc_create_m, only: shpc_create
+ use set_all_extr_m, only: set_all_extr
+ use set_contours_m, only: set_contours
+ use shpc_close_m, only: shpc_close
+ use write_eddy_m, only: write_eddy
+
+ implicit none
+
+ integer:: date = huge(0), eddy_i_target = huge(0)
+ integer ishape_first, d0, ishape, i, l, n1, n2, number_extr
+ logical:: cyclone = .true. ! orientation of target eddy
+ namelist /main_nml/ cyclone, date, eddy_i_target
+ TYPE(shpfileobject) hshp
+ character(len = :), allocatable:: slice_dir, slice_o_dir
+ type(eddy) e
+ integer extr_date ! field index in DBF file
+
+ integer, allocatable:: ishape_last(:) ! (d0:)
+ ! shape index (0-based) in the collection of shapefiles of the last
+ ! shape at a given date index
+
+ real corner(2)
+ ! longitude and latitude of the corner of the whole grid, in rad
+
+ real step(2) ! longitude and latitude steps, in rad
+ logical periodic ! grid is periodic in longitude
+
+ 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.
+
+ type(snapshot) s
+
+ real, allocatable:: coord(:, :)
+ ! (2, number_extr) longitude and latitude, in rad
+
+ integer, allocatable:: selection(:)
+ ! identifying numbers of a selection of eddies
+
+ logical, allocatable:: in_window(:)
+
+ real, allocatable:: outside_points(:, :) ! (2, :) longitude and
+ ! latitude, in rad, of all the significant extrema, except the
+ ! target extremum
+
+ TYPE(shpc_slice_handler) hshpc
+ integer llc(2) ! indices in global grid of lower left corner
+ integer urc(2) ! indices in global grid of upper right corner
+
+ real corner_window(2) ! longitude and latitude of the window around each
+ ! extremum, in rad
+
+ !------------------------------------------------------------------------
+
+ call config
+ call input_ssh(corner, step, periodic, ssh, u, v)
+ call get_command_arg_dyn(3, slice_dir, &
+ "Required 3rd argument: input-slice-directory")
+ write(unit = *, nml = main_nml)
+ print *, "examine_eddy: Enter main_nml:"
+ read(unit = *, nml = main_nml)
+ write(unit = *, nml = main_nml)
+
+ if (cyclone) then
+ slice_o_dir = slice_dir // "/Cyclones"
+ else
+ slice_o_dir = slice_dir // "/Anticyclones"
+ end if
+
+ call shp_open_03(hshp, pszshapefile = slice_o_dir // "/extremum", &
+ pszaccess = "rb")
+ call dbf_get_field_index_03(hshp, "date", extr_date)
+ call dbf_read_attribute_03(d0, hshp, extr_date, ishape = 0)
+ call read_column(ishape_last, file = slice_o_dir // "/ishape_last.txt", &
+ my_lbound = d0)
+
+ if (date == d0) then
+ ishape_first = 0
+ else
+ ishape_first = ishape_last(date - 1) + 1
+ end if
+
+ number_extr = ishape_last(date) - ishape_first + 1
+ allocate(coord(2, number_extr), in_window(number_extr))
+
+ do ishape = ishape_first, ishape_last(date)
+ call shp_read_point(coord(:, ishape - ishape_first + 1), hshp, ishape)
+ end do
+
+ coord = coord * deg_to_rad
+ e%extr%coord = coord(:, eddy_i_target)
+ e%extr%coord_proj = nint(convert_to_ind(e%extr%coord, corner, step))
+ e%extr%ssh = missing_ssh ! because we do not care about this
+ e%cyclone = cyclone
+
+ ! Define the geographical window around the target eddy extremum:
+
+ llc = e%extr%coord_proj - max_radius
+ urc = e%extr%coord_proj + max_radius
+
+ llc(2) = max(llc(2), 1)
+ urc(2) = min(urc(2), size(ssh, 2))
+
+ if (.not. periodic) then
+ llc(1) = max(llc(1), 1)
+ urc(1) = min(urc(1), size(ssh, 1))
+ end if
+
+ corner_window = corner + (llc - 1) * step
+
+ ! Done defining geographical window
+
+ do i = 1, number_extr
+ if (i /= eddy_i_target) then
+ ! Shift the longitude to a value close to the longitude of the
+ ! target extremum:
+ if (periodic) coord(1, i) = coord(1, i) &
+ + ceiling((corner_window(1) - coord(1, i)) / twopi) * twopi
+
+ in_window(i) = all(abs(nint((coord(:, i) - e%extr%coord) / step)) &
+ < max_radius)
+ else
+ in_window(i) = .false.
+ end if
+ end do
+
+ n2 = count(in_window)
+ call set_all_extr(s, step, periodic = .false., &
+ ssh = ssh(llc(1):urc(1), llc(2):urc(2)), corner = corner_window)
+ selection = pack(s%extr_map, s%extr_map /= 0)
+ selection = pack(selection, s%list(selection)%cyclone .neqv. cyclone)
+ n1 = size(selection)
+ allocate(outside_points(2, n1 + n2))
+ forall (i = 1:n1) outside_points(:, i) = s%list(selection(i))%extr%coord
+ l = n1 + 1
+
+ do i = 1, number_extr
+ if (in_window(i)) then
+ outside_points(:, l) = coord(:, i)
+ l = l + 1
+ end if
+ end do
+
+ i = s%extr_map(e%extr%coord_proj(1) - llc(1) + 1, &
+ e%extr%coord_proj(2) - llc(2) + 1)
+ e%innermost_level = s%list(i)%innermost_level
+ call create_cont_list
+ call set_contours(e%out_cont, e%speed_cont, e%max_speed, cyclone, e%extr, &
+ e%innermost_level, step, ssh(llc(1):urc(1), llc(2):urc(2)), &
+ u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), &
+ corner_window, pi * (min_radius * 1e3)**2, outside_points)
+ call close_cont_list
+ call shpc_create(hshpc, shpc_dir = "SHPC", cyclone = cyclone, slice = 0, &
+ grid_lon_lat = .true.)
+ call write_eddy(e, hshpc, date, 1)
+ write(hshpc%unit, fmt = *) 0
+ CALL shpc_close(hshpc)
+ print *, "llc = ", llc
+ print *, "urc = ", urc
+ print *, "corner_window (in degrees):", corner_window * rad_to_deg
+ print *, "n1 = ", n1
+ print *, "n2 = ", n2
+ print *, "outside_points = ", outside_points
+ print *, "e%innermost_level = ", e%innermost_level
+ print *, "Created shapefile cont_list and shapefile collection SHPC."
+
+end program examine_eddy
--
GitLab