Skip to content
Snippets Groups Projects
Commit 70fbeb95 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Do not use type snapshot in program `inst_eddies`

This was not useful since we do not need an array of snapshots in
`inst_eddies`. It is clearer without the integration in a snapshot
variable. The integration in a snapshot variable came from initial
commit, when there was only one program for detection of instantaneous
eddies and overlapping.
parent 60ab47a7
No related branches found
No related tags found
No related merge requests found
...@@ -11,7 +11,7 @@ program inst_eddies ...@@ -11,7 +11,7 @@ program inst_eddies
use numer_rec_95, only: indexx use numer_rec_95, only: indexx
use config_m, only: config, max_radius use config_m, only: config, max_radius
use derived_types, only: snapshot, shpc_slice_handler use derived_types, only: eddy, shpc_slice_handler
use input_ssh_m, only: input_ssh, nlon, nlat use input_ssh_m, only: input_ssh, nlon, nlat
use nearby_extr_m, only: nearby_extr use nearby_extr_m, only: nearby_extr
use set_all_extr_m, only: set_all_extr use set_all_extr_m, only: set_all_extr
...@@ -23,7 +23,20 @@ program inst_eddies ...@@ -23,7 +23,20 @@ program inst_eddies
implicit none implicit none
type(snapshot) s type(eddy), allocatable:: list(:) ! (number_extr) "Eddies" at a
! given date. In general, there may be both cyclonic and
! anticyclonic eddies in the list. Eddies include "eddies"
! without an outermost contour, that is, only the extremum is
! defined. The subscript value in list is the "identification
! number" of the extremum.
integer number_extr ! number of extrema
integer, allocatable:: extr_map(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is
! periodic in longitude, else (nlon, nlat). At a point of extremum
! SSH: identification number or this extremum. 0 at other points.
TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
integer iostat, i, l integer iostat, i, l
character(len = 200) iomsg character(len = 200) iomsg
...@@ -46,7 +59,7 @@ program inst_eddies ...@@ -46,7 +59,7 @@ program inst_eddies
integer:: date = 0, slice = 0 integer:: date = 0, slice = 0
namelist /dates_nml/ date, slice namelist /dates_nml/ date, slice
integer, allocatable:: sorted_extr(:) ! (s%number_extr) integer, allocatable:: sorted_extr(:) ! (number_extr)
! Sorted identifying number of extrema: first, minima sorted in ! Sorted identifying number of extrema: first, minima sorted in
! order of decreasing SSH, and second, maxima sorted in order of ! order of decreasing SSH, and second, maxima sorted in order of
! increasing SSH. ! increasing SSH.
...@@ -112,35 +125,35 @@ program inst_eddies ...@@ -112,35 +125,35 @@ program inst_eddies
write(unit, fmt = *) "CPU time for preamble, before computation:", t1 - t0, & write(unit, fmt = *) "CPU time for preamble, before computation:", t1 - t0, &
"s" "s"
t0 = t1 t0 = t1
call set_all_extr(s%extr_map, s%list, step, periodic, ssh) call set_all_extr(extr_map, list, step, periodic, ssh)
s%number_extr = size(s%list) number_extr = size(list)
allocate(sorted_extr(s%number_extr)) allocate(sorted_extr(number_extr))
forall (i = 1:s%number_extr) forall (i = 1:number_extr)
s%list(i)%out_cont%closed = .true. list(i)%out_cont%closed = .true.
! must be intialized to true because it is used in nearby_extr ! must be intialized to true because it is used in nearby_extr
end forall end forall
! Double sort of extrema on cyclonicity and SSH value: ! Double sort of extrema on cyclonicity and SSH value:
selection = argwhere(s%list%cyclone) selection = argwhere(list%cyclone)
n_cycl = size(selection) n_cycl = size(selection)
selection = selection(indexx(s%list(selection)%extr%ssh)) selection = selection(indexx(list(selection)%extr%ssh))
sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
selection = argwhere(.not. s%list%cyclone) selection = argwhere(.not. list%cyclone)
selection = selection(indexx(s%list(selection)%extr%ssh)) selection = selection(indexx(list(selection)%extr%ssh))
sorted_extr(n_cycl + 1:) = selection sorted_extr(n_cycl + 1:) = selection
! Done sorting ! Done sorting
do l = 1, s%number_extr do l = 1, number_extr
i = sorted_extr(l) i = sorted_extr(l)
! Define the geographical window around each eddy extremum: ! Define the geographical window around each eddy extremum:
llc = s%list(i)%extr%coord_proj - max_radius llc = list(i)%extr%coord_proj - max_radius
urc = s%list(i)%extr%coord_proj + max_radius urc = list(i)%extr%coord_proj + max_radius
llc(2) = max(llc(2), 1) llc(2) = max(llc(2), 1)
urc(2) = min(urc(2), nlat) urc(2) = min(urc(2), nlat)
...@@ -153,16 +166,16 @@ program inst_eddies ...@@ -153,16 +166,16 @@ program inst_eddies
! Done defining geographical window ! Done defining geographical window
outside_points & outside_points &
= nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), s%list, i) = nearby_extr(extr_map(llc(1):urc(1), llc(2):urc(2)), list, i)
! Shift the x coordinate of each point to a value close to the ! Shift the x coordinate of each point to a value close to the
! coordinate of the target extremum: ! coordinate of the target extremum:
if (periodic) outside_points(1, :) = outside_points(1, :) & if (periodic) outside_points(1, :) = outside_points(1, :) &
+ ceiling((llc(1) - outside_points(1, :)) / real(nlon)) * nlon + ceiling((llc(1) - outside_points(1, :)) / real(nlon)) * nlon
call set_contours(s%list(i)%out_cont, s%list(i)%speed_cont, & call set_contours(list(i)%out_cont, list(i)%speed_cont, &
s%list(i)%max_speed, s%list(i)%cyclone, s%list(i)%extr, & list(i)%max_speed, list(i)%cyclone, list(i)%extr, &
s%list(i)%innermost_level, step, ssh(llc(1):urc(1), llc(2):urc(2)), & list(i)%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)), & u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), &
llc, outside_points) llc, outside_points)
end do end do
...@@ -170,8 +183,8 @@ program inst_eddies ...@@ -170,8 +183,8 @@ program inst_eddies
call cpu_time(t1) call cpu_time(t1)
write(unit, fmt = *) "CPU time for computation, before output:", t1 - t0, "s" write(unit, fmt = *) "CPU time for computation, before output:", t1 - t0, "s"
t0 = t1 t0 = t1
print *, "inst_eddies: s%number_extr = ", s%number_extr print *, "inst_eddies: number_extr = ", number_extr
call write_snapshot(s%list, hshpc_cyclo, hshpc_anti, date) call write_snapshot(list, hshpc_cyclo, hshpc_anti, date)
CALL shpc_close(hshpc_cyclo) CALL shpc_close(hshpc_cyclo)
CALL shpc_close(hshpc_anti) CALL shpc_close(hshpc_anti)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment