SUBROUTINE EXTREM(X,N)
C
C     PURPOSE--THIS SUBROUTINE PERFORMS AN EXTREME VALUE ANALYSIS
C              ON THE DATA IN THE INPUT VECTOR X. 
C              THIS ANALYSIS CONSISTS OF DETERMINING THAT PARTICULAR
C              EXTREME VALUE TYPE 1 OR EXTREME VALUE TYPE 2 DISTRIBUTION
C              WHICH BEST FITS THE DATA SET.
C              THE GOODNESS OF FIT CRITERION IS THE MAXIMUM PROBABILITY
C              PLOT CORRELATION COEFFICIENT CRITERION.
C              AFTER THE BEST-FIT DISTRIBUTION IS DETERMINED,
C              ESTIMATES ARE COMPUTED AND PRINTED OUT FOR THE
C              LOCATION AND SCALE PARAMETERS.
C              TWO PROBABILITY PLOTS ARE ALSO PRINTED OUT-- 
C              THE BEST-FIT TYPE 2 PROBABILITY PLOT
C              (IF THE BEST FIT WAS IN FACT A TYPE 2),
C              AND THE TYPE 1 PROBABILITY PLOT.
C              PREDICTED EXTREMES FOR VARIOUS RETURN PERIODS ARE
C              ALSO COMPUTED AND PRINTED OUT.
C     INPUT ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                               (UNSORTED OR SORTED) OBSERVATIONS.
C                      N      = THE INTEGER NUMBER OF OBSERVATIONS
C                               IN THE VECTOR X.
C     OUTPUT--6 PAGES OF AUTOMATIC PRINTOUT.
C     PRINTING--YES.
C     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N
C                   FOR THIS SUBROUTINE IS 7500.
C     OTHER DATAPAC   SUBROUTINES NEEDED--SORT, UNIMED, EV1PLT,
C                                         EV2PLT, PLOT.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, ALOG.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     REFERENCES--FILLIBEN (1972), 'TECHNIQUES FOR TAIL LENGTH
C                 ANALYSIS', PROCEEDINGS OF THE EIGHTEENTH
C                 CONFERENCE ON THE DESIGN OF EXPERIMENTS IN
C                 ARMY RESEARCH AND TESTING, PAGES 425-450. 
C               --FILLIBEN, 'THE PERCENT POINT FUNCTION',
C                 UNPUBLISHED MANUSCRIPT.
C               --JOHNSON AND KOTZ (1970), CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS-1, 1970, PAGES 272-295.
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         --DECEMBER  1974. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --MAY       1976. 
C
C---------------------------------------------------------------------
C
      CHARACTER*4 BLANK,ALPHAM,ALPHAA,ALPHAX
      CHARACTER*4 ALPHAI,ALPHAN,ALPHAF,ALPHAT,ALPHAY
      CHARACTER*4 ALPHAG,EQUAL
C
      CHARACTER*4 IFLAG1
      CHARACTER*4 IFLAG2
      CHARACTER*4 IFLAG3
C
      DIMENSION W(3000)
      DIMENSION X(1)
      DIMENSION Y(7500),Z(7500)
      DIMENSION GAMTAB(50),CORR(50)
      DIMENSION YI(50),YS(50),T(50)
      DIMENSION IFLAG1(50),IFLAG2(50),IFLAG3(50)
CCCCC DIMENSION C(10)
      DIMENSION AM(50)
      DIMENSION SCRAT(50)
C
      DIMENSION AINDEX(50)
CCCCC DIMENSION P0(10)
      DIMENSION H(60,2)
      COMMON /BLOCK2/ WS(15000)
      EQUIVALENCE (Y(1),WS(1)),(Z(1),WS(7501))
      DATA BLANK,ALPHAM,ALPHAA,ALPHAX/' ','M','A','X'/
      DATA ALPHAI,ALPHAN,ALPHAF,ALPHAT,ALPHAY/'I','N','F','T','Y'/
      DATA ALPHAG,EQUAL/'G','='/
      DATA GAMTAB(1),GAMTAB(2),GAMTAB(3),GAMTAB(4),GAMTAB(5),
     1GAMTAB(6),GAMTAB(7),GAMTAB(8),GAMTAB(9),GAMTAB(10),
     1GAMTAB(11),GAMTAB(12),GAMTAB(13),GAMTAB(14),GAMTAB(15),
     1GAMTAB(16),GAMTAB(17),GAMTAB(18),GAMTAB(19),GAMTAB(20),
     1GAMTAB(21),GAMTAB(22),GAMTAB(23),GAMTAB(24),GAMTAB(25)
     1/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,
     113.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25./
      DATA GAMTAB(26),GAMTAB(27),GAMTAB(28),GAMTAB(29),GAMTAB(30),
     1GAMTAB(31),GAMTAB(32),GAMTAB(33),GAMTAB(34),GAMTAB(35),
     1GAMTAB(36),GAMTAB(37),GAMTAB(38),GAMTAB(39),GAMTAB(40),
     1GAMTAB(41),GAMTAB(42)
     1/30.,35.,40.,45.,50.,60.,70.,80.,90.,100.,150.,200.,250.,
     1350.,500.,750.,1000./
CCCCC DATA C(1),C(2),C(3),C(4),C(5),C(6),C(7),C(8),C(9),C(10)
CCCCC1/60.,75.,100.,150.,250.,500.,1000.,10000.,100000.,1000000./
CCCCC DATA P0(1),P0(2),P0(3),P0(4),P0(5),P0(6),P0(7),P0(8),P0(9),P0(10)
CCCCC1/.0,.5,.75,.9,.95,.975,.99,.999,.9999,.99999/
      DATA T(1),T(2),T(3),T(4),T(5),T(6),T(7),T(8),T(9),T(10),
     1T(11),T(12),T(13),T(14),T(15),T(16),T(17),T(18),T(19),T(20),
     1T(21),T(22),T(23),T(24),T(25)
     1/10.18011,3.39672,2.47043,2.14609,1.98712,1.89429,1.83394,
     11.79175,1.76069,1.73691,1.71814,1.70297,1.69045,1.67996,
     11.67103,1.66335,1.65667,1.65082,1.64564,1.64102,1.63689,
     11.63316,1.62979,1.62672,1.62391/
      DATA T(26),T(27),T(28),T(29),T(30),
     1T(31),T(32),T(33),T(34),T(35),T(36),T(37),T(38),T(39),T(40),
     1T(41),T(42),T(43)
     1/1.61287,1.60516,1.59947,1.59510,1.59164,1.58651,1.58289,
     11.58019,1.57811,1.57645,1.57152,1.56908,1.56763,1.56666,
     11.56546,1.56377,1.56330,1.56187/
      DATA AINDEX(1),AINDEX(2),AINDEX(3),AINDEX(4),AINDEX(5),
     1AINDEX(6),AINDEX(7),AINDEX(8),AINDEX(9),AINDEX(10),
     1AINDEX(11),AINDEX(12),AINDEX(13),AINDEX(14),AINDEX(15),
     1AINDEX(16),AINDEX(17),AINDEX(18),AINDEX(19),AINDEX(20),
     1AINDEX(21),AINDEX(22),AINDEX(23),AINDEX(24),AINDEX(25)
     1/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,
     113.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25./
      DATA AINDEX(26),AINDEX(27),AINDEX(28),AINDEX(29),AINDEX(30),
     1AINDEX(31),AINDEX(32),AINDEX(33),AINDEX(34),AINDEX(35),
     1AINDEX(36),AINDEX(37),AINDEX(38),AINDEX(39),AINDEX(40),
     1AINDEX(41),AINDEX(42),AINDEX(43),AINDEX(44),AINDEX(45),
     1AINDEX(46),AINDEX(47),AINDEX(48),AINDEX(49),AINDEX(50)
     1/26.,27.,28.,29.,30.,31.,32.,33.,34.,35.,36.,37.,38., 
     139.,40.,41.,42.,43.,44.,45.,46.,47.,48.,49.,50./
C
      IPR=6
      IUPPER=7500
      NUMDIS=43
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1.OR.N.GT.IUPPER)GOTO50
      IF(N.EQ.1)GOTO55
      HOLD=X(1)
      DO60I=2,N
      IF(X(I).NE.HOLD)GOTO90
   60 CONTINUE
      WRITE(IPR, 9)HOLD
      RETURN
   50 WRITE(IPR,17)IUPPER
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,18) 
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE EXTREM SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 EXTREM SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL *
     1****)
   18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME
     1NT TO THE EXTREM SUBROUTINE HAS THE VALUE 1 *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      AN=N
C
C     COMPUTE THE SAMPLE MINIMUM AND SAMPLE MAXIMUM
C
      XMIN=X(1)
      XMAX=X(1)
      DO140I=2,N
      IF(X(I).LT.XMIN)XMIN=X(I)
      IF(X(I).GT.XMAX)XMAX=X(I)
  140 CONTINUE
C
C     COMPUTE THE PROB PLOT CORRELATION COEFFICIENTS FOR THE VARIOUS VALUES
C     OF GAMMA
C
      CALL SORT(X,N,Y)
      CALL UNIMED(N,Z)
C
      DO100IDIS=1,NUMDIS
      IF(IDIS.EQ.NUMDIS)GOTO150
      A=GAMTAB(IDIS)
      DO110I=1,N
      W(I)=(-ALOG(Z(I)))**(-1.0/A)
  110 CONTINUE
      GOTO170
  150 DO160I=1,N
      W(I)=-ALOG(ALOG(1.0/Z(I)))
  160 CONTINUE
C
  170 SUM1=0.0
      SUM2=0.0
      DO200I=1,N
      SUM1=SUM1+Y(I)
      SUM2=SUM2+W(I)
  200 CONTINUE
      YBAR=SUM1/AN
      WBAR=SUM2/AN
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      DO300I=1,N
      SUM2=SUM2+(Y(I)-YBAR)*(W(I)-WBAR) 
      SUM1=SUM1+(Y(I)-YBAR)*(Y(I)-YBAR) 
      SUM3=SUM3+(W(I)-WBAR)*(W(I)-WBAR) 
  300 CONTINUE
      SY=SQRT(SUM1/(AN-1.0))
      CC=SUM2/SQRT(SUM3*SUM1) 
      YSLOPE=SUM2/SUM3
      YINT=YBAR-YSLOPE*WBAR
      CORR(IDIS)=CC 
      YI(IDIS)=YINT 
      YS(IDIS)=YSLOPE
  100 CONTINUE
C
C     DETERMINE THAT DISTRIBUTION WITH THE MAX PROB PLOT CORR COEFFICIENT
C
      IDISMX=1
      CORRMX=CORR(1)
      DO400IDIS=1,NUMDIS
      IF(CORR(IDIS).GT.CORRMX)IDISMX=IDIS
      IF(CORR(IDIS).GT.CORRMX)CORRMX=CORR(IDIS)
  400 CONTINUE
      DO500IDIS=1,NUMDIS
      IFLAG1(IDIS)=BLANK
      IFLAG2(IDIS)=BLANK
      IFLAG3(IDIS)=BLANK
      IF(IDIS.EQ.IDISMX)GOTO550
      GOTO500
  550 IFLAG1(IDIS)=ALPHAM
      IFLAG2(IDIS)=ALPHAA
      IFLAG3(IDIS)=ALPHAX
  500 CONTINUE
C
C     WRITE OUT THE TABLE OF PROB PLOT CORR COEFFICIENTS FOR VARIOUS GAMMA
C
      WRITE(IPR,998)
      WRITE(IPR,305)
      WRITE(IPR,999)
      WRITE(IPR,310)N
      WRITE(IPR,311)YBAR
      WRITE(IPR,312)SY
      WRITE(IPR,313)XMIN
      WRITE(IPR,314)XMAX
      WRITE(IPR,999)
      WRITE(IPR,323)
      WRITE(IPR,324)
      WRITE(IPR,325)
      WRITE(IPR,999)
C
      NUMDM1=NUMDIS-1
      IF(NUMDM1.LT.1)GOTO850
      DO800I=1,NUMDM1
      WRITE(IPR,805)GAMTAB(I),CORR(I),IFLAG1(I),IFLAG2(I),IFLAG3(I),
     1YI(I),YS(I),T(I)
  800 CONTINUE
  850 I=NUMDIS
      WRITE(IPR,806)ALPHAI,ALPHAN,ALPHAF,ALPHAI,ALPHAN,ALPHAI,
     1ALPHAT,ALPHAY,CORR(I),IFLAG1(I),IFLAG2(I),IFLAG3(I),
     1YI(I),YS(I),T(I)
C
C     PLOT THE PROB PLOT CORR COEFFICIENT VERSUS GAMMA VALUE INDEX
C
      CALL PLOT(CORR,AINDEX,NUMDIS)
      WRITE(IPR,810)ALPHAG,ALPHAA,ALPHAM,ALPHAM,ALPHAA,EQUAL,
     1GAMTAB(1),GAMTAB(12),GAMTAB(23),GAMTAB(34), 
     1ALPHAI,ALPHAN,ALPHAF,ALPHAI,ALPHAN,ALPHAI,ALPHAT,ALPHAY
      WRITE(IPR,999)
      WRITE(IPR,812)
      WRITE(IPR,813)
C
C     IF THE OPTIMAL GAMMA IS FINITE, PLOT OUT THE EXTREME VALUE
C     TYPE 2 PROBABILITY PLOT FOR THE OPTIMAL VALUE
C     OF GAMMA.
C
      IF(IDISMX.LT.NUMDIS)CALL EV2PLT(X,N,GAMTAB(IDISMX))
C
C     PLOT OUT AN EXTREME VALUE TYPE 1 PROBABILITY PLOT
C
      CALL EV1PLT(X,N)
C
C     FORM THE VARIOUS RETURN PERIOD VALUES
C
 1650 K=0 
      DO2100I=1,4
      DO2200J=1,9
      K=K+1
      AM(K)=J*(10**(I-1))
 2200 CONTINUE
 2100 CONTINUE
      K=K+1
      AM(K)=10000.
      K=K+1
      AM(K)=50000.
      K=K+1
      AM(K)=100000. 
      K=K+1
      AM(K)=500000. 
      K=K+1
      AM(K)=1000000.
      K=K+1
      AM(K)=N
      NUMAM=K
      CALL SORT(AM,NUMAM,SCRAT)
      DO2300I=1,NUMAM
      AM(I)=SCRAT(I)
 2300 CONTINUE
C
C     IF THE OPTIMAL GAMMA IS FINITE, COMPUTE THE 
C     PREDICTED EXTREME (= F(1-(1/M)) FOR VARIOUS RETURN PERIODS M
C     FOR THE OPTIMAL EXTREME VALUE TYPE 2 DISTRIBUTION.
C
      IF(IDISMX.EQ.NUMDIS)GOTO2450
      A=GAMTAB(IDISMX)
      YINT=YI(IDISMX)
      YSLOPE=YS(IDISMX)
      DO2400I=2,NUMAM
      R=1.0/AM(I)
      P=1.0-R
      ARG=-ALOG(P)
      IF(ARG.LE.0.0)GOTO2400
      H(I,1)=YINT+YSLOPE*(ARG**(-1.0/A))
 2400 CONTINUE
C
C     COMPUTE THE PREDICTED EXTREME (= F(1-(1/M)) FOR VARIOUS RETURN
C     PERIODS M FOR THE EXTREME VALUE TYPE 1 DISTRIBUTION.
C
 2450 YINT=YI(NUMDIS)
      YSLOPE=YS(NUMDIS)
      DO2500I=2,NUMAM
      R=1.0/AM(I)
      P=1.0-R
      ARG=-ALOG(P)
      IF(ARG.LE.0.0)GOTO2500
      H(I,2)=YINT+YSLOPE*(-ALOG(ARG))
 2500 CONTINUE
C
C     WRITE OUT THE PAGE WITH THE RETURN PERIODS AND THE PREDICTED EXTREMES
C     FOR THE 2 DISTRIBUTIONS--OPTIMAL EXTREME VALUE TYPE 2, AND EXTREME
C     VALUE TYPE 1. 
C
      WRITE(IPR,998)
      IF(IDISMX.EQ.NUMDIS)GOTO2750
      WRITE(IPR,2602)
      WRITE(IPR,2604)
      WRITE(IPR,2606)
      WRITE(IPR,2608)
      WRITE(IPR,2610)GAMTAB(IDISMX)
      WRITE(IPR,999)
      DO2700I=2,NUMAM
      WRITE(IPR,2705)AM(I),H(I,1),H(I,2)
      J=I-1
      JSKIP=J-5*(J/5)
      IF(JSKIP.EQ.0)WRITE(IPR,999)
 2700 CONTINUE
      RETURN
C
 2750 WRITE(IPR,2802)
      WRITE(IPR,2804)
      WRITE(IPR,2806)
      WRITE(IPR,2808)
      WRITE(IPR,999)
      DO2900I=2,NUMAM
      WRITE(IPR,2705)AM(I),H(I,2)
      J=I-1
      JSKIP=J-5*(J/5)
      IF(JSKIP.EQ.0)WRITE(IPR,999)
 2900 CONTINUE
C
  998 FORMAT(1H1)
  999 FORMAT(1H )
  305 FORMAT(1H ,40X,22HEXTREME VALUE ANALYSIS)
  310 FORMAT(1H ,37X,20HTHE SAMPLE SIZE N = ,I7)
  311 FORMAT(1H ,34X,18HTHE SAMPLE MEAN = ,F14.7) 
  312 FORMAT(1H ,28X,32HTHE SAMPLE STANDARD DEVIATION = ,F14.7)
  313 FORMAT(1H ,32X,21HTHE SAMPLE MINIMUM = ,F14.7)
  314 FORMAT(1H ,32X,21HTHE SAMPLE MAXIMUM = ,F14.7)
  323 FORMAT(1H ,85H     EXTREME VALUE      PROBABILITY PLOT     LOCATIO
     1N         SCALE       TAIL LENGTH)
  324 FORMAT(1H ,83H  TYPE 2 TAIL LENGTH      CORRELATION        ESTIMAT
     1E        ESTIMATE       MEASURE)
  325 FORMAT(1H ,37H   PARAMETER (GAMMA)      COEFFICIENT)
  805 FORMAT(1H ,3X,F10.2,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5) 
  806 FORMAT(1H ,5X,8A1,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5)
  810 FORMAT(1H ,12X,5A1,1X,A1,F14.7,11X,F14.7,11X,F14.7,11X,F14.7,
     115X,8A1)
  812 FORMAT(1H ,96HTHE ABOVE IS A PLOT OF THE 46 PROBABILITY PLOT CORRE
     1LATION COEFFICIENTS (FROM THE PREVIOUS PAGE))
  813 FORMAT(1H ,16X,41HVERSUS THE 46 EXTREME VALUE DISTRIBUTIONS)
 2602 FORMAT(1H ,43H   RETURN PERIOD     PREDICTED EXTREME WIND,
     1 27H     PREDICTED EXTREME WIND)
 2604 FORMAT(1H ,43H    (IN YEARS)          BASED ON OPTIMAL   ,
     1 20H            BASED ON)
 2606 FORMAT(1H ,42H                      EXTREME VALUE TYPE 2,
     1 27H       EXTREME VALUE TYPE 1)
 2608 FORMAT(1H ,43H                          DISTRIBUTION     ,
     1 22H          DISTRIBUTION)
 2610 FORMAT(1H ,30H                     (GAMMA = ,F12.5,1H))
 2705 FORMAT(1H ,2X,F9.1,13X,F10.2,17X,F10.2)
 2802 FORMAT(1H ,43H   RETURN PERIOD     PREDICTED EXTREME WIND)
 2804 FORMAT(1H ,36H    (IN YEARS)              BASED ON)
 2806 FORMAT(1H ,42H                      EXTREME VALUE TYPE 1)
 2808 FORMAT(1H ,38H                          DISTRIBUTION) 
C
      RETURN
      END