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