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