-
Lionel GUEZ authoredLionel GUEZ authored
examine_eddy.f90 6.25 KiB
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 (in degrees): "
do i = 1, n1 + n2
print *, outside_points(:, i) * rad_to_deg
end do
print *, "e%innermost_level = ", e%innermost_level
print *, "Created shapefile cont_list and shapefile collection SHPC."
end program examine_eddy