       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 
