program examine_eddy
! For a particular target instantaneous eddy, identified by its
! orientation, date and intra-date eddy index, the program reads in
! an SHPC the coordinates of all the extrema at the target date, and
! the SSH of the extremum of the target instantaneous eddy. The
! program also reads SSH, u, v at the target date.
! The program writes:
! -- the shapefiles cont_list and cont_list_proj for the target
! -- instantaneous eddy;
! -- an SHPC containing a single instantaneous eddy, that is the
! target instantaneous eddy;
! -- outside_points for the target instantaneous eddy;
! -- innermost_level for the target instantaneous eddy.
! The main difficulty here was to get outside_points. We need in
! outside_points all the extrema of opposite sign (including those
! not associated to an eddy) plus the extrema of same sign
! associated to an eddy. Note that e%extr%coord may not be exactly
! equal to what it was in inst_eddies because there is a
! multiplication by rad_to_deg when writing the shapefile and a
! multiplication by deg_to_rad when reading it.
! Libraries:
use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, new_unit, &
assert
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, max_radius
use cont_list_m, only: create_cont_list, close_cont_list
use derived_types, only: snapshot, shpc_slice_handler, eddy
use input_ssh_m, only: input_ssh, corner_whole, step, nlon, uniform_lon_lat
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, unit, n_arg
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 ifield ! 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
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 coord(2) ! longitude and latitude, in rad
integer, allocatable:: coord_proj(:, :)
! (2, number_extr) coordinates in projection space
integer, allocatable:: selection(:)
! identifying numbers of a selection of eddies
logical, allocatable:: in_window(:)
integer, allocatable:: outside_points(:, :) ! (2, :)
! coordinates in projection space 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
character(len = :), allocatable:: ssh_fname, u_fname, v_fname ! file names
!------------------------------------------------------------------------
n_arg = command_argument_count()
call assert(n_arg >= 2 .and. n_arg <= 4, &
"Usage: examine_eddy input-slice-directory h-file [u-file [v-file]]")
call get_command_arg_dyn(1, slice_dir)
call get_command_arg_dyn(2, ssh_fname)
if (n_arg == 2) then
u_fname = ssh_fname
v_fname = ssh_fname
else
call get_command_arg_dyn(3, u_fname)
if (n_arg == 3) then
v_fname = u_fname
else
call get_command_arg_dyn(4, v_fname)
end if
end if
call config
call input_ssh(periodic, ssh, u, v, ssh_fname, u_fname, v_fname)
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
if (uniform_lon_lat) then
call shp_open_03(hshp, pszshapefile = slice_o_dir // "/extremum", &
pszaccess = "rb")
else
call shp_open_03(hshp, pszshapefile = slice_o_dir // "/extr_proj", &
pszaccess = "rb")
end if
call dbf_get_field_index_03(hshp, "date", ifield)
call dbf_read_attribute_03(d0, hshp, ifield, 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_proj(2, number_extr), in_window(number_extr))
do ishape = ishape_first, ishape_last(date)
call shp_read_point(coord, hshp, ishape)
if (uniform_lon_lat) then
coord_proj(:, ishape - ishape_first + 1) = nint((coord * deg_to_rad &
- corner_whole) / step + 1.)
else
coord_proj(:, ishape - ishape_first + 1) = nint(coord)
end if
end do
e%extr%coord = (coord_proj(:, eddy_i_target) - 1) * step + corner_whole
e%extr%coord_proj = coord_proj(:, eddy_i_target)
call dbf_get_field_index_03(hshp, "ssh", ifield)
call dbf_read_attribute_03(e%extr%ssh, hshp, ifield, &
ishape = ishape_first + eddy_i_target - 1)
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
! 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_proj(1, i) = coord_proj(1, i) &
+ ceiling((llc(1) - coord_proj(1, i)) / real(nlon)) * nlon
in_window(i) = all(abs(coord_proj(:, i) - e%extr%coord_proj) &
< max_radius)
else
in_window(i) = .false.
end if
end do
n2 = count(in_window)
call set_all_extr(s%extr_map, s%list, periodic, ssh)
s%number_extr = size(s%list)
selection = pack(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
s%extr_map(llc(1):urc(1), llc(2):urc(2)) /= 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_proj
l = n1 + 1
do i = 1, number_extr
if (in_window(i)) then
outside_points(:, l) = coord_proj(:, i)
l = l + 1
end if
end do
i = s%extr_map(e%extr%coord_proj(1), e%extr%coord_proj(2))
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, 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)), llc, &
outside_points)
call close_cont_list
call shpc_create(hshpc, shpc_dir = "SHPC", cyclone = cyclone, slice = 0, &
with_proj = .false.)
call write_eddy(e, hshpc, date, 1)
write(hshpc%unit, fmt = *) 0
CALL shpc_close(hshpc)
print *, "llc = ", llc
print *, "urc = ", urc
print *, "Number of extrema of opposite orientation which must be outside ", &
"outermost contour:", n1
print *, "Number of (valid) eddies of same orientation which must be ", &
"outside outermost contour:", n2
call new_unit(unit)
open(unit, file = "outside_points.csv", status = "replace", action = "write")
write(unit, fmt = *) "x y"
do i = 1, n1 + n2
write(unit, fmt = *) outside_points(:, i)
end do
close(unit)
print *, "e%innermost_level = ", e%innermost_level
print *, "e%extr%coord_proj = ", e%extr%coord_proj
print *, "Created outside_points.csv, shapefiles cont_list and ", &
"cont_list_proj and shapefile collection SHPC."
end program examine_eddy