[cig-commits] [commit] devel, master: sorts station list in output_solver.txt by epicentral distance (afc569e)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:29:35 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit afc569e9a357993e816c71a048e72c020ae211d5
Author: daniel peter <peterda at ethz.ch>
Date: Mon Sep 1 12:46:38 2014 +0200
sorts station list in output_solver.txt by epicentral distance
>---------------------------------------------------------------
afc569e9a357993e816c71a048e72c020ae211d5
src/specfem3D/locate_receivers.f90 | 135 ++++++++++++++++++++++++++++++++++---
1 file changed, 125 insertions(+), 10 deletions(-)
diff --git a/src/specfem3D/locate_receivers.f90 b/src/specfem3D/locate_receivers.f90
index 6fee35e..44e6945 100644
--- a/src/specfem3D/locate_receivers.f90
+++ b/src/specfem3D/locate_receivers.f90
@@ -158,6 +158,8 @@
double precision :: lat_min,lat_max,lon_min,lon_max
! search margin in degrees
double precision,parameter :: LAT_LON_MARGIN = 2.d0
+ ! sorting order
+ integer, allocatable, dimension(:) :: irec_dist_ordered
! get MPI starting time
time_start = wtime()
@@ -205,7 +207,7 @@
y_found(nrec), &
z_found(nrec), &
final_distance(nrec),stat=ier)
- if (ier /= 0 ) call exit_MPI(myrank,'Error allocating temporary receiver arrays')
+ if (ier /= 0) call exit_MPI(myrank,'Error allocating temporary receiver arrays')
! initializes search distances
final_distance(:) = HUGEVAL
@@ -332,13 +334,6 @@
epidist(irec) = acos(cos(theta)*cos(theta_source) + &
sin(theta)*sin(theta_source)*cos(phi-phi_source))*RADIANS_TO_DEGREES
- ! print some information about stations
- if (myrank == 0) then
- write(IMAIN,*) 'Station #',irec,': ',station_name(irec)(1:len_trim(station_name(irec))), &
- '.',network_name(irec)(1:len_trim(network_name(irec))), &
- ' epicentral distance: ',sngl(epidist(irec)),' degrees'
- endif
-
! record three components for each station
do iorientation = 1,3
@@ -494,6 +489,28 @@
if (USE_DISTANCE_CRITERION ) deallocate(xyz_midpoints)
+ ! print some information about stations
+ if (myrank == 0) then
+ ! sorts stations according to epicentral distances
+ allocate(irec_dist_ordered(nrec),stat=ier)
+ if (ier /= 0) call exit_MPI(myrank,'Error allocating temporary irec_dist_ordered array')
+
+ ! sorts array
+ call heap_sort_distances(nrec,epidist,irec_dist_ordered)
+
+ ! outputs info
+ write(IMAIN,*) 'Stations sorted by epicentral distance:'
+ do i = 1,nrec
+ irec = irec_dist_ordered(i)
+ write(IMAIN,*) 'Station #',irec,': ',station_name(irec)(1:len_trim(station_name(irec))), &
+ '.',network_name(irec)(1:len_trim(network_name(irec))), &
+ ' epicentral distance: ',sngl(epidist(irec)),' degrees'
+ enddo
+
+ deallocate(irec_dist_ordered)
+ endif
+
+
! create RECORDHEADERS file with usual format for normal-mode codes
if (myrank == 0) then
@@ -819,7 +836,7 @@
if (DISPLAY_DETAILS_STATIONS .or. final_distance(irec) > 0.01d0) then
write(IMAIN,*)
- write(IMAIN,*) 'station # ',irec,' ',station_name(irec),network_name(irec)
+ write(IMAIN,*) 'Station #',irec,': ',trim(station_name(irec))//'.'//trim(network_name(irec))
write(IMAIN,*) ' original latitude: ',sngl(stlat(irec))
write(IMAIN,*) ' original longitude: ',sngl(stlon(irec))
write(IMAIN,*) ' epicentral distance: ',sngl(epidist(irec))
@@ -831,7 +848,7 @@
! add warning if estimate is poor
! (usually means receiver outside the mesh given by the user)
if (final_distance(irec) > THRESHOLD_EXCLUDE_STATION) then
- write(IMAIN,*) 'station # ',irec,' ',station_name(irec),network_name(irec)
+ write(IMAIN,*) 'Station #',irec,': ',trim(station_name(irec))//'.'//trim(network_name(irec))
write(IMAIN,*) '*****************************************************************'
if (NCHUNKS_VAL == 6) then
write(IMAIN,*) '***** WARNING: receiver location estimate is poor, therefore receiver excluded *****'
@@ -972,3 +989,101 @@
end subroutine locate_receivers
+! sorting routine left here for inlining
+!
+! Implementation of a Heap Sort Routine
+! Input
+! n = Input
+! Length of arrays
+! X_in = Input
+! Vector to be sorted
+! dimension(n)
+! Y = Output
+! Sorted Indices of vector X
+!
+! Example:
+! D = [ 4.0 3.0 1.0 2.0 ] on Input
+! Y = [ 1 2 3 4 ] Computed Internally (in order)
+!
+! X = [ 1.0 2.0 3.0 4.0 ] Computed Internally
+! Y = [ 3 4 2 1 ] on Output
+!
+ subroutine heap_sort_distances(N, X_in, Y)
+
+ implicit none
+ integer, intent(in) :: N
+ double precision, dimension(N), intent(in) :: X_in
+ integer, dimension(N), intent(out) :: Y
+
+ ! local parameters
+ double precision, dimension(N) :: X
+ double precision :: tmp
+ integer :: itmp
+ integer :: i
+
+ do i = 1,N
+ Y(i) = i
+ X(i) = X_in(i)
+ enddo
+
+ ! checks if anything to do
+ if (N < 2) return
+
+ ! builds heap
+ do i = N/2, 1, -1
+ call heap_sort_siftdown(i, N)
+ enddo
+
+ ! sorts array
+ do i = N, 2, -1
+ ! swaps last and first entry in this section
+ tmp = X(1)
+ X(1) = X(i)
+ X(i) = tmp
+ itmp = Y(1)
+ Y(1) = Y(i)
+ Y(i) = itmp
+
+ call heap_sort_siftdown(1, i - 1)
+ enddo
+
+ contains
+
+ subroutine heap_sort_siftdown(start, bottom)
+
+ implicit none
+
+ integer, intent(in) :: start, bottom
+
+ ! local parameters
+ integer :: i, j
+ double precision :: xtmp
+ integer :: ytmp
+
+ i = start
+ xtmp = X(i)
+ ytmp = Y(i)
+
+ j = 2 * i
+ do while (j <= bottom)
+ ! chooses larger value first in this section
+ if (j < bottom) then
+ if (X(j) <= X(j+1)) j = j + 1
+ endif
+
+ ! checks if section already smaller than initial value
+ if (X(j) < xtmp) exit
+
+ X(i) = X(j)
+ Y(i) = Y(j)
+ i = j
+ j = 2 * i
+ enddo
+
+ X(i) = xtmp
+ Y(i) = ytmp
+
+ end subroutine heap_sort_siftdown
+
+ end subroutine heap_sort_distances
+
More information about the CIG-COMMITS
mailing list