-
Lionel GUEZ authored
Use `SHP_triplet` instead of Snapshot for the name of the directory containing the three shapefiles extremum, `outermost_contour` and `max_speed_contour`. This is meant as a general name, analogous to a new file type. Generally, the shapefiles may contain eddies at different dates, not just a snapshot.
Lionel GUEZ authoredUse `SHP_triplet` instead of Snapshot for the name of the directory containing the three shapefiles extremum, `outermost_contour` and `max_speed_contour`. This is meant as a general name, analogous to a new file type. Generally, the shapefiles may contain eddies at different dates, not just a snapshot.
test_read_snapshot.f90 3.09 KiB
program test_read_snapshot
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
use nr_util, only: deg_to_rad, arth, assert
use shapelib, only: shpfileobject, shpclose
use shapelib_03, only: shp_open_03
use derived_types, only: snapshot
use init_shapefiles_m, only: init_shapefiles
use read_field_indices_m, only: read_field_indices
use read_snapshot_m, only: read_snapshot
use write_eddy_m, only: write_eddy
use write_extr_map_m, only: write_extr_map
implicit none
type(snapshot) s
TYPE(shpfileobject) hshp_extremum ! shapefile extremum
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour
TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour
integer i, k, copy
real:: corner(2) = [0.125, - 59.875]
! longitude and latitude of the corner of the whole grid, in degrees
real:: step(2) = [0.25, 0.25] ! longitude and latitude steps, in degrees
logical periodic ! grid is periodic in longitude
integer:: nlon = 120, nlat = 120
integer:: dist_lim = 12
! We look for an overlapping eddy at dist_lim (in grid points) of
! the first extremum.
namelist /main_nml/ corner, step, nlon, nlat, dist_lim
!-------------------------------------------------------------------------
write(unit = error_unit, nml = main_nml)
write(unit = error_unit, fmt = *) "Enter namelist main_nml."
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
! 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(360. / step(1)) == nlon
print *, "periodic = ", periodic
if (periodic) call assert(2 * dist_lim * step(1) < 180., &
"test_read_snapshot dist_lim")
call shp_open_03(hshp_extremum, "SHP_triplet_old/extremum", pszaccess = "rb")
call shp_open_03(hshp_outermost, "SHP_triplet_old/outermost_contour", &
pszaccess = "rb")
call shp_open_03(hshp_max_speed, "SHP_triplet_old/max_speed_contour", &
pszaccess = "rb")
call read_field_indices(hshp_extremum, hshp_outermost, hshp_max_speed)
call read_snapshot(s, k, hshp_extremum, hshp_outermost, hshp_max_speed, &
corner * deg_to_rad, step * deg_to_rad, nlon, nlat, periodic, dist_lim)
CALL shpclose(hshp_extremum)
CALL shpclose(hshp_outermost)
CALL shpclose(hshp_max_speed)
call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
! Write snapshot:
do i = 1, s%number_vis_extr
call write_eddy(s%list_vis(i), hshp_extremum, hshp_outermost, &
hshp_max_speed, k, i)
end do
print *, "Number of extrema:", s%number_vis_extr
CALL shpclose(hshp_extremum)
print *, 'Created shapefile "extremum".'
CALL shpclose(hshp_outermost)
print *, 'Created shapefile "outermost_contour".'
CALL shpclose(hshp_max_speed)
print *, 'Created shapefile "max_speed_contour".'
print *, "s%ind_extr:"
do i = 1, s%number_vis_extr
print *, s%ind_extr(:, i)
end do
copy = merge(dist_lim, 0, periodic)
call write_extr_map(arth(corner(1) - copy * step(1), step(1), nlon + 2 &
* copy), arth(corner(2), step(2), nlat), s%extr_map)
end program test_read_snapshot