Program MAIN c c Program to analyze wind data using the De Haan procedure c PARAMETER(MAXLNS=3000) PARAMETER(MAXPT=10000) PARAMETER(MINPT=50) PARAMETER(MINSIZ=5) PARAMETER(MAXR=10) PARAMETER(MAXTHR=25) C REAL XTEMP(MAXPT) INTEGER INDX(MAXPT) REAL R(MAXR),XR(MAXR) REAL MD(MAXTHR),MSDD(MAXTHR) INTEGER MTD(MAXTHR),MND(MAXTHR) REAL MXM(MAXTHR),MLF(MAXTHR) REAL MR(MAXTHR,MAXR) REAL VAL(13*MAXLNS) C REAL LAMBDA,XMAX,LDFCTR C CHARACTER*255 IFILE CHARACTER*80 IFORMT CHARACTER*27 INAME C LOGICAL IDEBUG LOGICAL IERROR C C STEP 1: Open the input/output files C Specify the interval size (to control correlation) C Specify the mean recurrence intervals C IDEBUG=.FALSE. CALL OPENFI(IFILE,NIN,NOUT1,NOUT2,MAXR,NUMR,R,INTSIZ, 1 IDEBUG,IERROR) IF(IERROR)GOTO9999 IPNTR=0 DO100I=1,NUMR IF(ABS(R(I) - 50.0).LE.1.0)IPNTR=I 100 CONTINUE IF(IPNTR.EQ.0)THEN NUMR=NUMR+1 IF(NUMR.GT.MAXR)NUMR=MAXR IPNTR=NUMR ENDIF C C STEP 2: Read in the daily wind data. C MAXLIN=13*MAXLNS CALL GETDAT(NIN,MAXLIN,MINPT,VAL,NDAYS,INAME,NYRS,IDEBUG,IERROR) IF(IERROR)GOTO9999 WRITE(NOUT1,*) 'Number of years of wind data: ',NYRS WRITE(NOUT1,'(A1)') ' ' WRITE(NOUT1,'(A1)') ' ' C C STEP 3: From the daily wind data, calculate the maximum windspeed C for each interval, then sort all the maximas. C CALL INTRVL(VAL,NDAYS,XTEMP,NINTS,INTSIZ,INDX) CALL SORT(XTEMP,NINTS,XTEMP) C C STEP 4: Start the loop for various threshold. The first C threshold is chosen so that the 15 largest values are used. C C If there are ties for the 15th higest value, all of these C points are included. C C U defines the threshold (all values above U are used). C C STEP 4A: Determine the intial value for U. C NUSE = 15 NSTART = NINTS - NUSE + 1 XCUTOFF = XTEMP(NSTART) DO910I=NSTART-1,1,-1 IF(XTEMP(I).EQ.XCUTOFF)THEN NSTART = I ELSE GOTO919 ENDIF 910 CONTINUE 919 CONTINUE C NFIRST=NSTART-1 U = XTEMP(NFIRST) C C STEP 4B: Perform MAXTHR iterations. These correspond to C the thresholds U, U-1, U-2, ... , U-24. C C UTEMP = threshold for current iteration C LAMBDA = rate of threshold crossing per year C DO 1000 ITER=1,MAXTHR C IF(ITER.GT.1) U = U - 1.0 NPTS=0 DO1109I=1,NINTS IF(XTEMP(I).GT.U)THEN NSTART=I GOTO1119 ENDIF 1109 CONTINUE 1119 CONTINUE NFIRST=NSTART-1 IF(NFIRST.LT.1)NFIRST=1 C UTEMP = XTEMP(NFIRST) NPTS = NINTS - NFIRST + 1 NABOVE = NPTS - 1 LAMBDA=REAL(NABOVE)/REAL(NYRS) C WRITE(NOUT1,*) 'U, UTEMP = ',U,UTEMP WRITE(NOUT1,*) 'NUMBER OF POINTS ABOVE THRESHOLD = ',NABOVE WRITE(NOUT1,*) 'LAMBDA = ',LAMBDA WRITE(NOUT1,'(A1)') ' ' WRITE(NOUT1,'(A1)') ' ' C C STEP 4C: Check for minimum number of points C IF(NPTS.LE.MINSIZ)THEN MD(ITER)=-99.999 MSDD(ITER)=-99.999 MND(ITER)=0 MTD(ITER)=0 GOTO1000 ENDIF C C STEP 4D: Analyze the data with the De Haan procedure C C The De Haan procedure estimates the tail length C parameter "C", the standard deviation of C, and C scale parameter "A" of the generalized Pareto C distribution. Details of the computation are given C in the reference below. C C Once GPD parameters are estimated, compute the C wind speeds (XR) with specified mean recurrence C intervals (R). Also calculate the nominal ultimate C wind speed for the reverse Weibull distribution C (XMAX) and the corresponding load factor. C CALL DEHAAN(XTEMP(NFIRST),NPTS,C,SDC,AMN1) IF(C.GE.0.0)THEN RHO1=1.0 ELSE RHO1=1.0/(1.0-C) ENDIF A = UTEMP*AMN1/RHO1 WRITE(NOUT1,*) 'USING DE HAAN PROCEDURE:' WRITE(NOUT1,*) ' A = ',A WRITE(NOUT1,*) ' C = ',C WRITE(NOUT1,*) ' SD C = ',SDC DO1410I=1,NUMR IF(IDEBUG)THEN WRITE(*,*) 'I, R(I) = ',I,R(I) ENDIF XR(I)=UTEMP-A*(1.0 - (LAMBDA*R(I))**C)/C WRITE(NOUT1,1411) R(I),XR(I) 1411 FORMAT(' R, XR = ',2F12.5) 1410 CONTINUE C IF(C .LT. 0.0) THEN XMAX = UTEMP - (A/C) LDFCTR = (XMAX/XR(IPNTR))**2 ELSE XMAX = 999.0 LDFCTR = 999.0 ENDIF WRITE(NOUT1,1818) XMAX WRITE(NOUT1,1819) LDFCTR WRITE(NOUT1,'(A)') ' ' WRITE(NOUT1,'(A)') ' ' 1818 FORMAT(' XMAX = ',F12.5) 1819 FORMAT(' LDFCTR = ',F12.5) MTD(ITER)=INT(U+0.5) MD(ITER)=C MND(ITER)=NPTS MSDD(ITER)=SDC DO 8012 JJ=1,NUMR MR(ITER,JJ)=XR(JJ) 8012 CONTINUE MXM(ITER)=XMAX MLF(ITER)=LDFCTR 1499 CONTINUE C 1000 CONTINUE C C STEP 5: Write the results to an ASCII files (for additional C processing in other programs) C IFORMT=' ' IFORMT='(I2,1X,I4,1X,F7.3,1X,F7.3,1X, (F7.3,1X),F8.3,1X,F8.3)' WRITE(IFORMT(30:31),'(I2)')NUMR WRITE(NOUT2,'(A27)') INAME(1:27) DO 5010 JJ=1,MAXTHR WRITE(NOUT2,IFORMT)MTD(JJ),MND(JJ),MD(JJ),MSDD(JJ), 1 (MR(JJ,KK),KK=1,NUMR),MXM(JJ),MLF(JJ) 5010 CONTINUE GOTO9999 C 9010 CONTINUE WRITE(*,*) 'ERROR READING PARAMETER FILE' GOTO9999 C 9020 CONTINUE WRITE(*,*) 'ERROR READING DATA FILE' GOTO9999 C 9999 CONTINUE STOP END SUBROUTINE OPENFI(IFILE,NIN,NOUT1,NOUT2,MAXR,NUMR,R,INTSIZ, 1 IDEBUG,IERROR) C C This subroutine is used to prompt for and open the input C file. C C Input Arguemnts: C C IFILE = the name of the input file C MAXR = maximum number of mean recurrence intervals C IDEBUG = debug switch C C Output Arguemnts: C C NIN = the unit number for the input file C NOUT1 = the unit number for the first output file C NOUT2 = the unit number for the second output file C NUMR = number of mean recurrence intervals C R = array to store means recurrence intervals C INTSIZ = interval size C IERROR = error flag C REAL R(*) C LOGICAL IERROR LOGICAL IDEBUG LOGICAL IPERD C CHARACTER*40 ID CHARACTER*255 IFILE,IFILE2,IFILE3 CHARACTER*1 IA2 C C Step 0: Initialize C IFILE=' ' IFILE2=' ' IFILE3=' ' NCFIL1=0 NCFIL2=0 NIN=11 NOUT1=12 NOUT2=13 IERROR=.FALSE. C C Step 1: Prompt for the file name and open the input file C C Determine number of characters in file name and C the position of the last "." in the name. C WRITE(*,*) 'Name of file containing wind data:' READ(*,'(A255)',END=9010,ERR=9020) IFILE C DO100J=255,1,-1 IF(IFILE(J:J).NE.' ')THEN NCFIL1=J GOTO119 ENDIF 100 CONTINUE GOTO9020 119 CONTINUE C IPERD=.FALSE. NCFIL2=NCFIL1 DO200J=NCFIL1,1,-1 IF(IFILE(J:J).EQ.'.')THEN NCFIL2=J IPERD=.TRUE. GOTO219 ENDIF 200 CONTINUE 219 CONTINUE C IF(IDEBUG)THEN WRITE(*,*) 'NCFIL1,NCFIL2,IPERD=',NCFIL1,NCFIL2,IPERD ENDIF OPEN(NIN,FILE=IFILE,ERR=9030) C C Step 2: Open the output files C IFILE2(1:NCFIL2)=IFILE(1:NCFIL2) IF(IPERD)THEN IFILE2(NCFIL2+1:NCFIL2+3)='out' ELSE IFILE2(NCFIL2+1:NCFIL2+4)='.out' ENDIF OPEN(NOUT1,FILE=IFILE2,ERR=9040) C IFILE3(1:NCFIL2)=IFILE(1:NCFIL2) IF(IPERD)THEN IFILE3(NCFIL2+1:NCFIL2+3)='dhn' ELSE IFILE3(NCFIL2+1:NCFIL2+4)='.dhn' ENDIF OPEN(NOUT2,FILE=IFILE3,ERR=9050) C IF(IDEBUG)THEN WRITE(*,*) 'AFTER ALL FILES OPENED: IERROR=',IERROR ENDIF C C Now prompt for: C 1) Interval to use (typically 4 or 8) C 2) Mean recurrence intervals C WRITE(*,*) WRITE(*,*) 'The daily wind data may have highly correlated.' WRITE(*,*) 'To minimize this, we split the data into intervals' WRITE(*,*) 'and find the maximums of these intervals.' WRITE(*,*) WRITE(*,*) 'Enter the size of the interval (4 or 8 days ', 1 'recommended):' READ(*,*) INTSIZ IF(INTSIZ.LT.1)INTSIZ=4 IF(INTSIZ.GT.16)INTSIZ=16 C WRITE(*,*) WRITE(*,*) 'Specify the mean recurrence intervals to use:' WRITE(*,*) ' 0 - use default (25, 50, 100, 500, 1000, 2000)' WRITE(*,*) ' -1 - use (50, 500, 2000, 10000, 100000)' WRITE(*,*) ' -2 - specify which intervals to use' READ(*,*) ITEMP IF(ITEMP.EQ.0)THEN NUMR=6 R(1)=25. R(2)=50. R(3)=100. R(4)=500. R(5)=1000. R(6)=2000. ELSEIF(ITEMP.EQ.-1)THEN NUMR=5 R(1)=50. R(2)=500. R(3)=2000. R(4)=10000. R(5)=100000. ELSE NUMR=0 600 CONTINUE WRITE(*,*) 'Enter mean recurrence interval ',NUMR+1, 1 ' (0 to quit):' READ(*,*)ITEMP IF(ITEMP.GT.0)THEN NUMR=NUMR+1 R(NUMR)=REAL(ITEMP) IF(NUMR.GE.MAXR)GOTO9999 GOTO600 ELSE IF(NUMR.LT.1)THEN NUMR=6 R(1)=25. R(2)=50. R(3)=100. R(4)=500. R(5)=1000. R(6)=2000. ELSE GOTO9999 ENDIF ENDIF ENDIF GOTO9999 C 9010 CONTINUE WRITE(*,*) 'No file name entered' IERROR=.TRUE. GOTO9999 C 9020 CONTINUE WRITE(*,*) 'Error entering the file name' IERROR=.TRUE. GOTO9999 C 9030 CONTINUE WRITE(*,*) 'Error opening the input data file' IERROR=.TRUE. GOTO9999 C 9040 CONTINUE WRITE(*,*) 'Error opening the first output file' GOTO9999 C 9050 CONTINUE WRITE(*,*) 'Error opening the second output file' GOTO9999 C 9999 CONTINUE C IF(IDEBUG)THEN WRITE(*,*) 'INTSIZ = ',INTSIZ WRITE(*,*) 'MAXR, NUMR = ',MAXR,NUMR DO9003I=1,NUMR WRITE(*,*) 'I, R(I) = ',I,R(I) 9003 CONTINUE ENDIF C RETURN END SUBROUTINE GETDAT(NIN,MAXPT,MINPT,Y,N,INAME,NYRS,IDEBUG,IERROR) C C This subroutine is used to read the data. The default version C is used to read the Milepost data. For other wind data, C make the appropriate modifications to this routine. C C Input Arguemnts: C C NIN = unit number to read data from (file previosly C opened in OPENFI) C MAXPT = maximum number of points that can be read C MINPT = minimum number of points required for the analysis C C Output Arguemnts: C C Y = the data array that will contain the wind speeds C N = the number of data values actually read C NYRS = the number of years in the data file C INAME = Name of site (a city) C IERROR = error flag C REAL Y(MAXPT) C LOGICAL IERROR LOGICAL IDEBUG C CHARACTER*27 INAME CHARACTER*80 IA C PARAMETER (MAXLNS=3000) C C Step 0: Initialize C IA=' ' IERROR=.FALSE. C C Step 1: Read the city C READ(NIN,'(A80)',END=9010,ERR=9020)IA INAME(1:27)=IA(1:27) READ(NIN,'(A80)',END=9010,ERR=9020)IA C C Step 2: Read the daily wind data and determine the total number C of years and the total number of days in the sample C IFRST = 0 N = 0 DO 100 I=1,MAXLNS READ(NIN,'(A80)',END=119) IA IF(IA(5:5) .EQ. '/') THEN READ(IA(1:4),'(I4)') IYR IF(IFRST .EQ. 0) IFRST = IYR ILAST = IYR GOTO 100 ELSE DO 628 J = 1,13 NF = (J-1)*6 + 4 NL = NF + 2 IF(IA(NF:NL) .EQ. ' ') GOTO 100 N = N+1 IF(N.GT.MAXPT)THEN WRITE(*,*) 'Maximum number of points (',MAXPT, 1 ') exceeded.' IERROR=.TRUE. GOTO9999 ENDIF READ(IA(NF:NL),'(F3.0)') Y(N) 628 CONTINUE ENDIF 100 CONTINUE 119 CONTINUE IF(N.LT.MINPT)THEN WRITE(*,*) 'Minimum number of points (',MINPT,' not read.' IERROR=.TRUE. GOTO9999 ENDIF NYRS = ILAST-IFRST+1 C GOTO9999 C 9010 CONTINUE WRITE(*,*) 'End-of-File encountered while reading the header' WRITE(*,*) 'information.' IERROR=.TRUE. GOTO9999 C 9020 CONTINUE WRITE(*,*) 'Error while reading the header information.' IERROR=.TRUE. GOTO9999 C 9030 CONTINUE IF(N.LT.MINPT)THEN WRITE(*,*) 'End-of-File encountered before minimum' WRITE(*,*) 'number of observations (',MINPT,') read.' IERROR=.TRUE. ELSE WRITE(*,*) N,' wind speeds read.' IERROR=.FALSE. ENDIF GOTO9999 C 9040 CONTINUE WRITE(*,*) 'Error reading the wind speeds data.' WRITE(*,*) N,' records were successfully read.' IERROR=.TRUE. GOTO9999 C 9999 CONTINUE C IF(IDEBUG)THEN WRITE(*,*) 'N, NYRS = ',N,NYRS DO9003I=1,N WRITE(18,*) 'I, Y(I) = ',I,Y(I) 9003 CONTINUE ENDIF C RETURN END SUBROUTINE INTRVL(DWINDS,NUMDAYS,PWINDS,NUMINTS,INTSIZ,INDX) C C This subroutine extracts from the daily wind data DWINDS a reduced C list PWINDS containing the maximum windspeed in each of NUMINTS C intervals of length INTSIZ days. To reduce the correlation of C these data, the subroutine checks to see if maxima from adjacent C intervals are within CORSIZ days of each other, and changes the C lower of these two (using a method described below) if the C possibility of correlation exists. C REAL DWINDS(*),PWINDS(*) INTEGER INDX(*) C INTEGER CORSIZ LOGICAL CORFLAG C CORSIZ = INTSIZ/2 - 1 NUMINTS = NUMDAYS/INTSIZ IF(MOD(NUMDAYS,INTSIZ) .NE. 0) NUMINTS = NUMINTS+1 C C Find the max windspeed in each interval. C DO 1000 J=1,NUMINTS IFIRST = INTSIZ*(J-1) + 1 ILAST = INTSIZ*J IF(ILAST .GT. NUMDAYS) ILAST = NUMDAYS INDX(J) = MAXPTR(DWINDS,IFIRST,ILAST) PWINDS(J) = DWINDS(INDX(J)) 1000 CONTINUE C C Check for correlations between maxima of adjacent intervals. More C precisely, check for "backward" correlations, ie., correlations C between the maximum of the current interval and that of the C previous interval. C DO 2000 I=2,NUMINTS DIST = INDX(I) - INDX(I-1) IF(DIST .LE. CORSIZ) THEN C C If a correlation does exist, we keep the larger of the two and C replace the smaller one with the maximum of all data in the C relevant interval that are more than CORSIZ days away from the C kept maximum. C IF(PWINDS(I) .LE. PWINDS(I-1)) THEN C C In this situation, the maximum of the current interval is not C larger than the maximum of the previous interval. Therefore, we C replace the current maximum. C NEWFIRST = INDX(I-1)+CORSIZ+1 NEWLAST = I*INTSIZ IF(NEWLAST .GT. NUMDAYS) NEWLAST = NUMDAYS IF(NEWFIRST .GT. NEWLAST) NEWFIRST = (I-1)*INTSIZ + 1 INDX(I) = MAXPTR(DWINDS,NEWFIRST,NEWLAST) PWINDS(I) = DWINDS(INDX(I)) ELSE C C In this situation, the current maximum is larger than the C previous one. Therefore, we replace the previous maximum and C then check that this does not introduce a correlation with C *its* previous maximum. C DO 5454 II = I-1,1,-1 NEWFIRST = (II-1)*INTSIZ+1 NEWLAST = INDX(II+1)-CORSIZ-1 INDX(II) = MAXPTR(DWINDS,NEWFIRST,NEWLAST) PWINDS(II) = DWINDS(INDX(II)) IF(II .EQ. 1) GOTO 5455 IF(INDX(II)-INDX(II-1) .GT. CORSIZ) GOTO 5455 5454 CONTINUE 5455 CONTINUE ENDIF ENDIF 2000 CONTINUE C RETURN END INTEGER FUNCTION MAXPTR(VALUES,JF,JL) C C A function that returns an index MAXPTR such that C VALUES(MAXPTR) .GE. VALUES(J) for all J in the interval (JF,JL). C REAL VALUES(*) C VMAX = VALUES(JF) MAXPTR = JF DO 1111 J = JF,JL IF(VALUES(J) .LT. VMAX) GOTO 1111 VMAX = VALUES(J) MAXPTR = J 1111 CONTINUE RETURN END SUBROUTINE DEHAAN(X,N,GAMMA,SD,ANM1) C C SUBROUTINE IMPLEMENTING THE DE HAAN-DEKKER MOMENT-BASED EXTREME C VALUE INDEX ESTIMATOR AS DOCUMENTED IN C "Extreme Value Theory and Applications", Edited by C Galambos, Lechner, and Simiu, pp. 93-122, Kluwer Academic C Publishers, Boston, 1994. C PARAMETER(MAXPT=2000) C REAL GAMNUM,GAMDEN,GAMMA REAL X(*) C C THE MAIN LOOP C COMPUTE THE DE HAAN-DEKKER INDEX C (CALLED GAMMA BELOW) C FOR THE K HIGHEST ORDER STATISTICS, ITERATING ON K. C NI=N KK = NI - 1 C C GAMNUM AND GAMDEN ARE THE MN(1) AND MN(2) ON PAGE 100 OF C THE REFERENCE CITED ABOVE. C GAMNUM=0. GAMDEN=0. C DO 50 J=1,KK C JM1=J-1 TERM1=LOG(X(NI-JM1))-LOG(X(NI-KK)) GAMNUM=GAMNUM+TERM1 GAMDEN=GAMDEN+TERM1*TERM1 C 50 CONTINUE C GAMNUM=GAMNUM/FLOAT(KK) GAMDEN=GAMDEN/FLOAT(KK) ANM1=GAMNUM ANM2=GAMDEN C TERM1=ANM1**2/ANM2 GAMMA=ANM1+1.0-0.5*(1.0/(1.0-TERM1)) C C COMPUTE THE STANDARD DEVIATION OF C. C IF(GAMMA.GE.0.)THEN SD=SQRT((1.0 + GAMMA*GAMMA)/REAL(KK)) ELSE TERM1 = (1.0 - GAMMA)*(1.0 - GAMMA)*(1.0 - 2.0*GAMMA) TERM2 = 4.0 - 8.0*(1.0 - 2.0*GAMMA)/(1.0 - 3.0*GAMMA) TERM3 = (5.0-11.0*GAMMA)*(1.0-2.0*GAMMA)/ * ((1.0-3.0*GAMMA)*(1.0 - 4.0*GAMMA)) SD = SQRT(TERM1*(TERM2 + TERM3)/REAL(KK)) ENDIF C RETURN END SUBROUTINE SORT(X,N,Y) C C PURPOSE--THIS SUBROUTINE SORTS (IN ASCENDING ORDER) C THE N ELEMENTS OF THE SINGLE PRECISION VECTOR X C AND PUTS THE RESULTING N SORTED VALUES INTO THE C SINGLE PRECISION VECTOR Y. C INPUT ARGUMENTS--X = THE SINGLE PRECISION VECTOR OF C OBSERVATIONS TO BE SORTED. C --N = THE INTEGER NUMBER OF OBSERVATIONS C IN THE VECTOR X. C OUTPUT ARGUMENTS--Y = THE SINGLE PRECISION VECTOR C INTO WHICH THE SORTED DATA VALUES C FROM X WILL BE PLACED. C OUTPUT--THE SINGLE PRECISION VECTOR Y C CONTAINING THE SORTED C (IN ASCENDING ORDER) VALUES C OF THE SINGLE PRECISION VECTOR X. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THE DIMENSIONS OF THE VECTORS IL AND IU C (DEFINED AND USED INTERNALLY WITHIN C THIS SUBROUTINE) DICTATE THE MAXIMUM C ALLOWABLE VALUE OF N FOR THIS SUBROUTINE. C IF IL AND IU EACH HAVE DIMENSION K, C THEN N MAY NOT EXCEED 2**(K+1) - 1. C FOR THIS SUBROUTINE AS WRITTEN, THE DIMENSIONS C OF IL AND IU HAVE BEEN SET TO 36, C THUS THE MAXIMUM ALLOWABLE VALUE OF N IS C APPROXIMATELY 137 BILLION. C SINCE THIS EXCEEDS THE MAXIMUM ALLOWABLE C VALUE FOR AN INTEGER VARIABLE IN MANY COMPUTERS, C AND SINCE A SORT OF 137 BILLION ELEMENTS C IS PRESENTLY IMPRACTICAL AND UNLIKELY, C THEN THERE IS NO PRACTICAL RESTRICTION C ON THE MAXIMUM VALUE OF N FOR THIS SUBROUTINE. C (IN LIGHT OF THE ABOVE, NO CHECK OF THE C UPPER LIMIT OF N HAS BEEN INCORPORATED C INTO THIS SUBROUTINE.) C OTHER DATAPAC SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--NONE. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--THE SMALLEST ELEMENT OF THE VECTOR X C WILL BE PLACED IN THE FIRST POSITION C OF THE VECTOR Y, C THE SECOND SMALLEST ELEMENT IN THE VECTOR X C WILL BE PLACED IN THE SECOND POSITION C OF THE VECTOR Y, ETC. C COMMENT--THE INPUT VECTOR X REMAINS UNALTERED. C COMMENT--IF THE ANALYST DESIRES A SORT 'IN PLACE', C THIS IS DONE BY HAVING THE SAME C OUTPUT VECTOR AS INPUT VECTOR IN THE CALLING SEQUENCE. C THUS, FOR EXAMPLE, THE CALLING SEQUENCE C CALL SORT(X,N,X) C IS ALLOWABLE AND WILL RESULT IN C THE DESIRED 'IN-PLACE' SORT. C COMMENT--THE SORTING ALGORTHM USED HEREIN C IS THE BINARY SORT. C THIS ALGORTHIM IS EXTREMELY FAST AS THE C FOLLOWING TIME TRIALS INDICATE. C THESE TIME TRIALS WERE CARRIED OUT ON THE C UNIVAC 1108 EXEC 8 SYSTEM AT NBS C IN AUGUST OF 1974. C BY WAY OF COMPARISON, THE TIME TRIAL VALUES C FOR THE EASY-TO-PROGRAM BUT EXTREMELY C INEFFICIENT BUBBLE SORT ALGORITHM HAVE C ALSO BEEN INCLUDED-- C NUMBER OF RANDOM BINARY SORT BUBBLE SORT C NUMBERS SORTED C N = 10 .002 SEC .002 SEC C N = 100 .011 SEC .045 SEC C N = 1000 .141 SEC 4.332 SEC C N = 3000 .476 SEC 37.683 SEC C N = 10000 1.887 SEC NOT COMPUTED C REFERENCES--CACM MARCH 1969, PAGE 186 (BINARY SORT ALGORITHM C BY RICHARD C. SINGLETON). C --CACM JANUARY 1970, PAGE 54. C --CACM OCTOBER 1970, PAGE 624. C --JACM JANUARY 1961, PAGE 41. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-2315 C ORIGINAL VERSION--JUNE 1972. C UPDATED --NOVEMBER 1975. C C--------------------------------------------------------------------- C DIMENSION X(1),Y(1) DIMENSION IU(36),IL(36) C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1)GOTO50 IF(N.EQ.1)GOTO55 HOLD=X(1) DO60I=2,N IF(X(I).NE.HOLD)GOTO90 60 CONTINUE WRITE(IPR, 9)HOLD DO61I=1,N Y(I)=X(I) 61 CONTINUE RETURN 50 WRITE(IPR,15) WRITE(IPR,47)N RETURN 55 WRITE(IPR,18) Y(1)=X(1) RETURN 90 CONTINUE 9 FORMAT(1H ,108H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT (A VECTOR) TO THE SORT SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 SORT SUBROUTINE IS NON-POSITIVE *****) 18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME 1NT TO THE SORT SUBROUTINE HAS THE VALUE 1 *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C C COPY THE VECTOR X INTO THE VECTOR Y DO100I=1,N Y(I)=X(I) 100 CONTINUE C C CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED C NM1=N-1 DO200I=1,NM1 IP1=I+1 IF(Y(I).LE.Y(IP1))GOTO200 GOTO250 200 CONTINUE RETURN 250 M=1 I=1 J=N 305 IF(I.GE.J)GOTO370 310 K=I MID=(I+J)/2 AMED=Y(MID) IF(Y(I).LE.AMED)GOTO320 Y(MID)=Y(I) Y(I)=AMED AMED=Y(MID) 320 L=J IF(Y(J).GE.AMED)GOTO340 Y(MID)=Y(J) Y(J)=AMED AMED=Y(MID) IF(Y(I).LE.AMED)GOTO340 Y(MID)=Y(I) Y(I)=AMED AMED=Y(MID) GOTO340 330 Y(L)=Y(K) Y(K)=TT 340 L=L-1 IF(Y(L).GT.AMED)GOTO340 TT=Y(L) 350 K=K+1 IF(Y(K).LT.AMED)GOTO350 IF(K.LE.L)GOTO330 LMI=L-I JMK=J-K IF(LMI.LE.JMK)GOTO360 IL(M)=I IU(M)=L I=K M=M+1 GOTO380 360 IL(M)=K IU(M)=J J=L M=M+1 GOTO380 370 M=M-1 IF(M.EQ.0)RETURN I=IL(M) J=IU(M) 380 JMI=J-I IF(JMI.GE.11)GOTO310 IF(I.EQ.1)GOTO305 I=I-1 390 I=I+1 IF(I.EQ.J)GOTO370 AMED=Y(I+1) IF(Y(I).LE.AMED)GOTO390 K=I 395 Y(K+1)=Y(K) K=K-1 IF(AMED.LT.Y(K))GOTO395 Y(K+1)=AMED GOTO390 END