program aggregate_gradient_chk_v3
c *******************************************************************
c * *
c * Program to take results from the maximum elevation gradient *
c * checker program for a number of storms and aggregate those *
c * results to output the element number, the number of storms *
c * that involve that element, an average of the gradient and also *
c * the x and y coordinates of the point in the middle of the *
c * element. This is for use with reading in the file into SMS as *
c * a scatterpoint so one can see the elements that have gradient *
c * warnings. *
c * *
c * The program writes the element numbers out in reverse numerical *
c * order of the number of storms involved for each element (so *
c * the first element listed has the most storms involved). *
c * *
c * The input file is named input_file.in and the format is: *
c * Line 1: number of files to open and aggregate *
c * Lines 2-#files+1: filename to open *
c * *
c * This program also reads in the grid file using the standard *
c * name of fort.14. It produces two output files, agg_ele.out *
c * which is the main output file and storm_name.out which lists *
c * the elements and storm numbers. Further modification of this *
c * program is expected to have a better way to output the storm *
c * names. *
c * *
c * Crystal W. Fulcher *
c * v3 - 5/25/2011 *
c * *
c *******************************************************************
c Set up variable names
character*60 filename
c character*60 storm(1201000,1395)
real*8 avg_x,avg_y,sum_x,sum_y
real*8, allocatable :: gradient(:,:)
real*8, allocatable :: grad_sum(:)
real*8, allocatable :: x(:)
real*8, allocatable :: y(:)
real*8 grad_avg
integer, allocatable :: nele_node(:,:)
integer, allocatable :: nele_aggregate(:)
integer, allocatable :: nele_sort(:)
c Open the files used in this program
open(10,file='input_files.in')
open(11,file='agg_ele.out')
open(13,file='storm_name.out')
open(14,file='fort.14')
c Read in the number of storm files to read in and allocate the variable
read(10,*) nstorms
allocate ( nele_aggregate(nstorms) )
c Read in the grid file header lines and allocate variables
read(14,*)
read(14,*) nele,nnodes
allocate ( nele_node(nele,3) )
allocate ( gradient(nstorms,nele) )
allocate ( grad_sum(nele) )
allocate ( nele_sort(nele) )
allocate ( x(nnodes) )
allocate ( y(nnodes) )
c Read in node table from grid (fort.14) file
do i=1,nnodes
read(14,*) n,x(i),y(i),z
enddo
c Read in element table from grid (fort.14) file, zero out variables and set sort variable
c to equal the element number
do i=1,nele
read(14,*) n,ndum,nele_node(i,1),nele_node(i,2),nele_node(i,3)
nele_aggregate(i)=0
grad_sum(i)=0
nele_sort(i)=i
enddo
c Loop to read in each storm file and then increment the number for each storm that contains
c an element, also add the gradient to the gradient sum for a final averaging in the end
c This loop also writes out the element and the filename to a separate file, this is something I need
c to work on further and perfect
do i=1,nstorms
read(10,10,end=20) filename
open(12,file=filename)
do j=1,10000
read(12,*,end=5) num_ele,gradient(i,num_ele)
nele_aggregate(num_ele)=nele_aggregate(num_ele)+1
grad_sum(num_ele)=grad_sum(num_ele)+gradient(i,num_ele)
write(13,15) num_ele,filename
enddo
5 continue
close(12)
enddo
10 format(a60)
15 format(i8,2x,a60)
20 continue
c Call the sorting routine to organize the element numbers in the order of how many storms were
c involved for each element, this will order from 1 to ....
CALL SORT3(nele,nele_aggregate,nele_sort,grad_sum)
c Zero out the sum_x and sum_y for finding the middle of the element
sum_x=0.d0
sum_y=0.d0
c Prepare to write out the element number in reverse order of how many storms the element had a gradient
c warning in (so the element with the largest number of storms involved is first). This loop also
c calculates the gradient average and then writes out the element number, the number of storms that
c have gradient warnings for that element, the average gradient for the element, the midpoint of the element
c x-coordinate and the midpoint of the element y-coordinate. If there were no storms with gradient warnings
c for an element, it is not written out.
do i=nele,1,-1
if (nele_aggregate(i).eq.0) then
goto 30
endif
grad_avg=grad_sum(i)/nele_aggregate(i)
do j=1,3
sum_x=sum_x+x(nele_node(nele_sort(i),j))
sum_y=sum_y+y(nele_node(nele_sort(i),j))
enddo
avg_x=sum_x/3.d0
avg_y=sum_y/3.d0
sum_x=0.d0
sum_y=0.d0
write(11,35) nele_sort(i),nele_aggregate(i),grad_avg,avg_x,
& avg_y
30 continue
enddo
35 format(i10,2x,i4,2x,f12.10,2x,f13.9,2x,f12.9)
write(*,*) 'it has finished the last loop'
END
SUBROUTINE SORT3(N,RA,RB,RC)
c Sorts an array RA of length N into ascending numerical order using the Heapsort algorithm.
c while making the corresponding rearrangement of the array RB and RC.
c
c This subroutine is slightly modified from Numerical Recipes
INTEGER RA(N),RB(N),RRA,RRB
REAL*8 RC(N), RRC
L=N/2+1
IR=N
c The index L will be decremented from its inital value down to 1 during the "hiring" (heap
c creation phase. Once it reaches 1, the index IR will be decremented from its initial value down
c to 1 during the "retirement-and-promotion" (heap selection) phase.
10 CONTINUE
IF (L.GT.1) THEN
L=L-1
RRA=RA(L)
RRB=RB(L)
RRC=RC(L)
ELSE
RRA=RA(IR)
RRB=RB(IR)
RRC=RC(IR)
RA(IR)=RA(1)
RB(IR)=RB(1)
RC(IR)=RC(1)
IR=IR-1
IF (IR.EQ.1) THEN
RA(1)=RRA
RB(1)=RRB
RC(1)=RRC
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF (J.LE.IR) THEN
IF (J.LT.IR) THEN
IF (RA(J).LT.RA(J+1)) J=J+1
ENDIF
IF (RRA.LT.RA(J)) THEN
RA(I)=RA(J)
RB(I)=RB(J)
RC(I)=RC(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
RA(I)=RRA
RB(I)=RRB
RC(I)=RRC
GO TO 10
END