      Program HURRLDF
c
c  Program to analyze hurricane wind data using the De Haan procedure
c
c  Although this program is specific to the Milepost data sets,
c  it can easily be adapted to other data sets.  The main routine
c  that needs to be modified is the GETDAT routine (which actually
c  reads the data).  You can also modify the desired reccurence
c  intervals (change NUMR and the values specifed for R).
c
c     MAXPT   = Maximum number of observations
c     MINPT   = Minimum number of observations required for analysis
c     NUMR    = Number of recurrence intervals
c     MAXIT = Maximum number of iterations
c
      PARAMETER(MAXPT=10000)
      PARAMETER(MINPT=25)
      PARAMETER(NUMR=4)
      PARAMETER(MAXIT=50)
C
      REAL XTEMP(MAXPT)
      REAL R(NUMR),XR(NUMR),LDFCTR(NUMR)
      REAL MXR(NUMR), MLF(NUMR)
      REAL LAMBDA,XMAX
      REAL MC, MSDC, MXMAX
      INTEGER MTHRSH,MNXEED
C
      CHARACTER*255 IFILE
      CHARACTER*40  MPOST
      CHARACTER*40  IFORMT
C
      LOGICAL IERROR
      LOGICAL IDEBUG
C
      DATA R/50, 2000, 10000, 100000/
      DATA IERROR/.FALSE./
      DATA IDEBUG/.FALSE./
C
C  UNCOMMENT FOLLOWING LINE TO ACTIVATE SOME DEBUG CODE
C
CCCCC IDEBUG=.TRUE.
C
C  Step 1: Open the input/output files
C
      CALL OPENFI(IFILE,NIN,NOUT1,NOUT2,IDEBUG,IERROR)
      IF (IDEBUG)THEN
         WRITE(*,*)'NIN,NOUT1,NOUT2,IERROR=',NIN,NOUT1,NOUT2,IERROR
      END IF
      IF(IERROR)GOTO9999
C
C  Step 2: Read the data
C
C          Note that this routine needs to be modifed if you
C          are not using the Milepost data
C
C          MPOST is the Milepost identifier
C          URATE is the estimated mean annual rate of occurrence
C                of hurricanes/tropical storms for this milepost
C          XTEMP is the maximum observed wind speed (regardless of
c                direction) for the ith simulation
C
      CALL GETDAT(NIN,MAXPT,MINPT,XTEMP,NSTRMS,MPOST,NCID,URATE,
     1            IERROR)
      IF (IDEBUG)THEN
         WRITE(*,*)'NSTRMS,NCID,URATE=',NSTRMS,NCID,URATE
         WRITE(*,*)'XTEMP(1) = ',XTEMP(1)
         WRITE(*,*)'XTEMP(NSTRMS) = ',XTEMP(NSTRMS)
      END IF
      IF(IERROR)GOTO9999
C
C  Step 3: Write some initial information
C
      WRITE(NOUT1,10) MPOST(1:NCID)
      WRITE(NOUT1,20) URATE
      WRITE(NOUT1,30) NSTRMS
      WRITE(NOUT1,40)
      WRITE(NOUT1,'(A)') '      '
   10 FORMAT('MILEPOST:                     ',A4)
   20 FORMAT('URATE:                        ',F7.4)
   30 FORMAT('Number of simulated storms:   ',I4)
   40 FORMAT(1X)
C
C  Step 4: Start loop for thresholds here
C
C          1. Sort the data
C          2. Specify the initial threshold so that top MINPT.
C             observations selected (ties will be included).
C          3. U is the highest threshold considered.
C
      CALL SORT(XTEMP,NSTRMS,XTEMP)
C
      NUSE = 25
      NSTART = NSTRMS - NUSE + 1
      XCUTOFF = XTEMP(NSTART)
      DO110I=NSTART-1,1,-1
        IF(XTEMP(I).EQ.XCUTOFF)THEN
          NSTART = I
        ELSE
          GOTO119
        ENDIF
 110  CONTINUE
 119  CONTINUE
C
      NFIRST=NSTART-1
      U = XTEMP(NFIRST) + 1.0
      IFORMT=' '
      IFORMT(1:36)='(I2,1X,I4,1X,2(F7.3,1X),  (F9.3,1X))'
      WRITE(IFORMT(25:26),'(I2)') 2*NUMR+1
      WRITE(NOUT2,'(A4)') MPOST(1:MIN(4,NCID))
C
C  Step 5: Main computational loop
C
C          Create MAXIT sets corresponding the thresholds
C          U, U-1, U-2, ..., U-MAXIT.
C
      DO 1000 ITER=1,MAXIT
C
C        Step 5a: Calculate for the next threshold
C
         U = U - 1.0
         NPTS=0
         DO1109I=1,NSTRMS
            IF(XTEMP(I).GT.U)THEN
               NSTART=I
               GOTO1119
            ENDIF
 1109    CONTINUE
 1119    CONTINUE
         NFIRST=NSTART-1
         IF(NFIRST.LT.1)NFIRST=1
         IF (IDEBUG)THEN
            WRITE(*,*)'NFIRST,U = ',NFIRST,U
         ENDIF
C
C        Step 5b: UTEMP is the threshold for this iteration
C                 LAMBDA is the rate of threshold crossings per year.
C
         UTEMP = XTEMP(NFIRST)
         NPTS = NSTRMS - NFIRST + 1
         NABOVE = NPTS - 1
         LAMBDA=REAL(NABOVE)*URATE/REAL(NSTRMS)
C
         WRITE(NOUT1,*) 'U, UTEMP = ',U,UTEMP
         WRITE(NOUT1,*) 'NUMBER OF POINTS ABOVE THRESHOLD = ',NABOVE
         WRITE(NOUT1,*) 'LAMBDA                           = ',LAMBDA
         WRITE(NOUT1,40)
C
C        Step 5c: Check for insufficient points (< 6)
C
         IF(NPTS.LE.5)THEN
            MC=-99.999
            MSDC=-99.999
            MNXEED=0
            MTHRSH=0
            GOTO1000
         ENDIF
C
C        Step 5d: Analyze the data using the De Haan procedure
C
C                 The De Haan routine returns the tail length
C                 parameter "C" of the generalized Pareto
C                 distribution and SDC is the standard deviation
C                 of the tail length parameter.
C
C                 A is the scale parameter of the generalized
C                 Pareto distribution (see equation 5.3 on
C                 page 108 of De Haan  "Extreme Value Statistics",
C                 complete reference given in DEHAAN subroutine).
C                 XR is the wind speed with an R-year mean
C                 recurrence interval.
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
         IF (IDEBUG)THEN
            WRITE(*,*) 'C,SDC,A = ',C,SDC,A
         ENDIF
         WRITE(NOUT1,*) 'USING DE HAAN PROCEDURE:'
         WRITE(NOUT1,*) '     A = ',A
         WRITE(NOUT1,*) '     C = ',C
         WRITE(NOUT1,*) '  SD C = ',SDC
         DO1410I=1,NUMR
            XR(I)=UTEMP-A*(1.0 - (LAMBDA*R(I))**C)/C
            WRITE(NOUT1,1411)R(I),XR(I)
 1410    CONTINUE
 1411    FORMAT('R, XR   = ',2F12.5)
C
C        Step 5e: Based on the just-determined values of A, C, and
C                 XR(50), calculate
C
C                 XMAX   = the maximum possible wind speed for this
C                          probability distribution
C                 LDFCTR = the corresponding load factor
C
         IF(C .LT. 0.0) THEN
           XMAX = UTEMP - (A/C)
           LDFCTR(NUMR) = (XMAX/XR(1))**2
         ELSE
           XMAX = 999.
           LDFCTR(NUMR) = 999.
         ENDIF
         DO 1503 IK=2,NUMR
            LDFCTR(IK-1) = (XR(IK)/XR(1))**2
 1503    CONTINUE
         WRITE(NOUT1,1818) XMAX
         WRITE(NOUT1,1819) (LDFCTR(IK),IK=1,NUMR)
         WRITE(NOUT1,'(A)') '      '
         WRITE(NOUT1,'(A)') '      '
 1818    FORMAT(' XMAX   = ',F12.5)
 1819    FORMAT(' Load Factors = ',4(F10.5,1X))
         MTHRSH=INT(U+0.5)
         MC=C
         MNXEED=NABOVE
         MSDC=SDC
         DO1450JJ=1,NUMR
            MXR(JJ)=XR(JJ)
 1450    CONTINUE
         MXMAX=XMAX
         DO 1513 JJ=1,NUMR
            MLF(JJ)=LDFCTR(JJ)
 1513    CONTINUE
C
C        Step 5f: Write results to file
C
         WRITE(NOUT2,IFORMT)MTHRSH,MNXEED,MC,MSDC,
     1                      (MXR(KK),KK=1,NUMR),MXMAX,
     1                      (MLF(IK),IK=1,NUMR)
 1000 CONTINUE
C
 9999 CONTINUE
      STOP
      END
