extrem(3f) - [M_datapac:STATISTICS] determine whether a type 1 or
type 2 extreme value distribution better fits a given data set
SUBROUTINE EXTREM(X,N)
extrem(3f) performs an extreme value analysis on the data in the
input vector x.
this analysis consists of determining that particular extreme value
type 1 or extreme value type 2 distribution which best fits the
data set.
the goodness of fit criterion is the maximum probability plot
correlation coefficient criterion.
after the best-fit distribution is determined, estimates are computed
and printed out for the location and scale parameters.
two probability plots are also printed out-- the best-fit type 2
probability plot (if the best fit was in fact a type 2), and the type
1 probability plot.
predicted extremes for various return periods are also computed and
printed out.
X description of parameter
Y description of parameter
Sample program:
program demo_extrem
use M_datapac, only : extrem
implicit none
! call extrem(x,y)
end program demo_extrem
Results:
The original DATAPAC library was written by James Filliben of the Statistical
Engineering Division, National Institute of Standards and Technology.
John Urban, 2022.05.31
CC0-1.0
fcdf(3f) - [M_datapac:CUMULATIVE_DISTRIBUTION] compute the F cumulative distribution
function
SUBROUTINE FCDF(X,Nu1,Nu2,Cdf)
REAL(kind=wp) :: X
INTEGER :: Nu1
INTEGER :: Nu2
REAL(kind=wp) :: Cdf
FCDF(3f) computes the cumulative distribution function value for the F
distribution with integer degrees of freedom parameters = NU1 and NU2.
This distribution is defined for all non-negative X. The probability
density function is given in the references below.
X The value at which the cumulative distribution function is to
be evaluated. X should be non-negative.
NU1 The integer degrees of freedom for the numerator of the F
ratio. NU1 should be positive.
NU2 The integer degrees of freedom for the denominator of the F
ratio. NU2 should be positive.
CDF The cumulative distribution function value for the F distribution
Sample program:
program demo_fcdf
use M_datapac, only : fcdf
implicit none
! call fcdf(x,y)
end program demo_fcdf
Results:
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
John Urban, 2022.05.31
CC0-1.0
fourie(3f) - [M_datapac:ANALYSIS] perform a Fourier analysis of a data set
SUBROUTINE FOURIE(X,N)
REAL(kind=wp),intent(in) :: X(:)
INTEGER :: N
FOURIE(3f) performs a Fourier analysis of the data in the input vector
X. The analysis consists of the following--
1. computing (and printing)
(for each of the harmonic frequencies
1/n, 2/n, 3/n, ..., 1/2)
the corresponding fourier coefficients,
the amplitude, the phase,
the contribution to the total variance,
and the relative contribution to the total
variance.
2. plotting out a fourier line spectrum =
the periodogram = the plot of relative
contribution to total variance
(at each fourier frequency) versus
the fourier frequency.
In order that the results of the Fourier analysis be valid and properly
interpreted, the input data in X should be equi-spaced in time (or
whatever variable corresponds to time).
The horizontal axis of the spectra produced by fourie(3f) is frequency.
This frequency is measured in units of cycles per 'data point' or,
more precisely, in cycles per unit time where 'unit time' is defined
as the elapsed time between adjacent observations.
The range of the frequency axis is 0.0 to 0.5.
Fourier analysis differs from spectral analysis (as, for example,
produced by the datapac TIMESE(3f) subroutine) in that a Fourier
analysis does no smoothing on the spectral estimates whereas a spectral
analysis does smooth the spectral estimates. The net result is that
the spectral estimates obtained from a Fourier analysis are almost
always more variable than those obtained in a spectral analysis.
The practical conclusion is that when the data analyst has a choice
of whether to perform a Fourier analysis or a spectral analysis,
the spectral analysis should almost always be preferred.
the maximum number of Fourier frequencies for which the Fourier
coefficients is computed (and listed) is N/2 where N is the sample
size (length of the data record in the vector X). This rule is
overridden (for listing purposes only) in large data sets and is
replaced by the rule that the maximum number of lags listed = 800
(which corresponds to an 8-page listing of Fourier coefficients.
If more pages are desired, change the value of the variable MAXPAG
within this subroutine from 8 to whatever is desired.
If the input observations in X are considered to have been collected
1 second apart in time, then the frequency axis of the resulting
spectra would be in units of Hertz (= cycles per second).
The frequency of 0.0 corresponds to a cycle in the data of infinite
(= 1/(0.0)) length or period. the frequency of 0.5 corresponds to
a cycle in the data of length = 1/(0.5) = 2 data points.
Any equi-spaced fourier analysis is intrinsically limited to detecting
frequencies no larger than 0.5 cycles per data point; this corresponds
to the fact that the smallest detectable cycle in the data is 2 data
points per cycle.
X The vector of (unsorted) observations.
N The integer number of observations in the vector X.
The maximum allowable value of N for this subroutine is 15000.
The sample size N must be greater than or equal to 3.
2 to 10 pages (depending on the input sample size) of automatic
printout--
1. a listing of the amplitude, phase, contribution to the total
variance, and relative contribution to the total variance for
each of the fourier frequencies (1/n, 2/n, 3/n, ..., 1/2).
this listing may take as little as 1 page or as many as n/100
pages (the exact number depending on the input sample size n).
this listing is terminated after at most 8 computer pages.
if more pages are desired, change the value of the variable
maxpag within this subroutine from 8 to whatever desired.
2. a plot of the relative contribution to the total variance versus
frequency.
Sample program:
program demo_fourie
use M_datapac, only : fourie
implicit none
! call fourie(x,y)
end program demo_fourie
Results:
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
John Urban, 2022.05.31
CC0-1.0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | dimension(:) | :: | X | |||
| integer | :: | N |
| Type | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|
| real | :: | WS(15000) |
| Type | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|
| real | :: | WS(15000) |
SUBROUTINE EXTREM(X,N) REAL(kind=wp) :: a, aindex, am, an, arg, cc, corr, corrmx, gamtab, h, hold, p, r, scrat, sum1, sum2, sum3, sy, t, w REAL(kind=wp) :: wbar, WS, X, xmax, xmin, Y, ybar, yi, yint, ys, yslope, Z INTEGER :: i, idis, idismx, iupper, j, jskip, k, N, numam, numdis, numdm1 ! ! INPUT ARGUMENTS--X = THE VECTOR OF ! (UNSORTED OR SORTED) OBSERVATIONS. ! N = THE INTEGER NUMBER OF OBSERVATIONS ! IN THE VECTOR X. ! OUTPUT--6 pages OF AUTOMATIC PRINTOUT. ! PRINTING--YES. ! RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N ! FOR THIS SUBROUTINE IS 7500. ! OTHER DATAPAC SUBROUTINES NEEDED--SORT, UNIMED, EV1PLT, ! EV2PLT, PLOT. ! FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, LOG. ! MODE OF INTERNAL OPERATIONS--. ! ORIGINAL VERSION--JUNE 1972. ! UPDATED --DECEMBER 1974. ! UPDATED --NOVEMBER 1975. ! UPDATED --MAY 1976. ! !--------------------------------------------------------------------- ! CHARACTER(len=4) :: blank , alpham , alphaa , alphax CHARACTER(len=4) :: alphai , alphan , alphaf , alphat , alphay CHARACTER(len=4) :: alphag , equal ! CHARACTER(len=4) :: iflag1 CHARACTER(len=4) :: iflag2 CHARACTER(len=4) :: iflag3 ! DIMENSION w(3000) DIMENSION X(:) DIMENSION Y(7500) , Z(7500) DIMENSION gamtab(50) , corr(50) DIMENSION yi(50) , ys(50) , t(50) DIMENSION iflag1(50) , iflag2(50) , iflag3(50) DIMENSION am(50) DIMENSION scrat(50) ! DIMENSION aindex(50) DIMENSION h(60,2) COMMON /BLOCK2_real32/ WS(15000) EQUIVALENCE (Y(1),WS(1)) EQUIVALENCE (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) , & & gamtab(6) , gamtab(7) , gamtab(8) , gamtab(9) , gamtab(10) , & & gamtab(11) , gamtab(12) , gamtab(13) , gamtab(14) , & & gamtab(15) , gamtab(16) , gamtab(17) , gamtab(18) , & & gamtab(19) , gamtab(20) , gamtab(21) , gamtab(22) , & & gamtab(23) , gamtab(24) , gamtab(25)/1.0_wp , 2.0_wp , 3.0_wp , 4.0_wp , 5.0_wp ,& & 6.0_wp , 7.0_wp , 8.0_wp , 9.0_wp , 10.0_wp , 11.0_wp , 12.0_wp , 13.0_wp , 14.0_wp , 15.0_wp , 16.0_wp ,& & 17.0_wp , 18.0_wp , 19.0_wp , 20.0_wp , 21.0_wp , 22.0_wp , 23.0_wp , 24.0_wp , 25.0_wp/ DATA gamtab(26) , gamtab(27) , gamtab(28) , gamtab(29) , & & gamtab(30) , gamtab(31) , gamtab(32) , gamtab(33) , & & gamtab(34) , gamtab(35) , gamtab(36) , gamtab(37) , & & gamtab(38) , gamtab(39) , gamtab(40) , gamtab(41) , & & gamtab(42)/30.0_wp , 35.0_wp , 40.0_wp , 45.0_wp , 50.0_wp , 60.0_wp , 70.0_wp , 80.0_wp , & & 90.0_wp , 100.0_wp , 150.0_wp , 200.0_wp , 250.0_wp , 350.0_wp , 500.0_wp , 750.0_wp , 1000.0_wp/ !CCCC DATA C(1),C(2),C(3),C(4),C(5),C(6),C(7),C(8),C(9),C(10) !CCCC1/60.0_wp,75.0_wp,100.0_wp,150.0_wp,250.0_wp,500.0_wp,1000.0_wp,10000.0_wp,100000.0_wp,1000000.0_wp/ !CCCC DATA P0(1),P0(2),P0(3),P0(4),P0(5),P0(6),P0(7),P0(8),P0(9),P0(10) !CCCC1/0.0_wp,0.5_wp,0.75_wp,0.9_wp,0.95_wp,0.975_wp,0.99_wp,0.999_wp,0.9999_wp,0.99999_wp/ DATA t(1) , t(2) , t(3) , t(4) , t(5) , t(6) , t(7) , t(8) , & & t(9) , t(10) , t(11) , t(12) , t(13) , t(14) , t(15) , & & t(16) , t(17) , t(18) , t(19) , t(20) , t(21) , t(22) , & & t(23) , t(24) , t(25)/10.18011_wp , 3.39672_wp , 2.47043_wp , & & 2.14609_wp , 1.98712_wp , 1.89429_wp , 1.83394_wp , 1.79175_wp , 1.76069_wp , & & 1.73691_wp , 1.71814_wp , 1.70297_wp , 1.69045_wp , 1.67996_wp , 1.67103_wp , & & 1.66335_wp , 1.65667_wp , 1.65082_wp , 1.64564_wp , 1.64102_wp , 1.63689_wp , & & 1.63316_wp , 1.62979_wp , 1.62672_wp , 1.62391_wp/ DATA t(26) , t(27) , t(28) , t(29) , t(30) , t(31) , t(32) , & & t(33) , t(34) , t(35) , t(36) , t(37) , t(38) , t(39) , & & t(40) , t(41) , t(42) , t(43)/1.61287 , 1.60516 , 1.59947 , & & 1.59510_wp , 1.59164_wp , 1.58651_wp , 1.58289_wp , 1.58019_wp , 1.57811_wp , & & 1.57645_wp , 1.57152_wp , 1.56908_wp , 1.56763_wp , 1.56666_wp , 1.56546_wp , & & 1.56377_wp , 1.56330_wp , 1.56187_wp/ DATA aindex(1) , aindex(2) , aindex(3) , aindex(4) , aindex(5) , & & aindex(6) , aindex(7) , aindex(8) , aindex(9) , aindex(10) , & & aindex(11) , aindex(12) , aindex(13) , aindex(14) , & & aindex(15) , aindex(16) , aindex(17) , aindex(18) , & & aindex(19) , aindex(20) , aindex(21) , aindex(22) , & & aindex(23) , aindex(24) , aindex(25)/1.0_wp , 2.0_wp , 3.0_wp , 4.0_wp , 5.0_wp ,& & 6.0_wp , 7.0_wp , 8.0_wp , 9.0_wp , 10.0_wp , 11.0_wp , 12.0_wp , 13.0_wp , 14.0_wp , 15.0_wp , 16.0_wp ,& & 17.0_wp , 18.0_wp , 19.0_wp , 20.0_wp , 21.0_wp , 22.0_wp , 23.0_wp , 24.0_wp , 25.0_wp/ DATA aindex(26) , aindex(27) , aindex(28) , aindex(29) , & & aindex(30) , aindex(31) , aindex(32) , aindex(33) , & & aindex(34) , aindex(35) , aindex(36) , aindex(37) , & & aindex(38) , aindex(39) , aindex(40) , aindex(41) , & & aindex(42) , aindex(43) , aindex(44) , aindex(45) , & & aindex(46) , aindex(47) , aindex(48) , aindex(49) , & & aindex(50)/26.0_wp , 27.0_wp , 28.0_wp , 29.0_wp , 30.0_wp , 31.0_wp , 32.0_wp , 33.0_wp , & & 34.0_wp , 35.0_wp , 36.0_wp , 37.0_wp , 38.0_wp , 39.0_wp , 40.0_wp , 41.0_wp , 42.0_wp , 43.0_wp , & & 44.0_wp , 45.0_wp , 46.0_wp , 47.0_wp , 48.0_wp , 49.0_wp , 50.0_wp/ ! iupper = 7500 numdis = 43 ! ! CHECK THE INPUT ARGUMENTS FOR ERRORS ! IF ( N<1 .OR. N>iupper ) THEN WRITE (G_IO,99001) iupper 99001 FORMAT (' ', & &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE EXTREM SUBROU& &TINE IS OUTSIDE THE ALLOWABLE (1,',I0,') INTERVAL *****') WRITE (G_IO,99002) N 99002 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****') RETURN ELSE IF ( N==1 ) THEN WRITE (G_IO,99003) 99003 FORMAT (' ', & &'***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUMENT TO THE EXTR& &EM SUBROUTINE HAS THE VALUE 1 *****') RETURN ELSE hold = X(1) DO i = 2 , N IF ( X(i)/=hold ) GOTO 50 ENDDO WRITE (G_IO,99004) hold 99004 FORMAT (' ', & &'***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUMENT (A VECTOR) & &TO THE EXTREM SUBROUTINE HAS ALL ELEMENTS = ',E15.8,' *****') RETURN ENDIF ! !-----START POINT----------------------------------------------------- ! 50 an = N ! ! COMPUTE THE SAMPLE MINIMUM AND SAMPLE MAXIMUM ! xmin = X(1) xmax = X(1) DO i = 2 , N IF ( X(i)<xmin ) xmin = X(i) IF ( X(i)>xmax ) xmax = X(i) ENDDO ! ! COMPUTE THE PROB PLOT CORRELATION COEFFICIENTS FOR THE VARIOUS VALUES ! OF GAMMA ! CALL SORT(X,N,Y) CALL UNIMED(N,Z) ! DO idis = 1 , numdis IF ( idis==numdis ) THEN DO i = 1 , N w(i) = -LOG(LOG(1.0_wp/Z(i))) ENDDO ELSE a = gamtab(idis) DO i = 1 , N w(i) = (-LOG(Z(i)))**(-1.0_wp/a) ENDDO ENDIF ! sum1 = 0.0_wp sum2 = 0.0_wp DO i = 1 , N sum1 = sum1 + Y(i) sum2 = sum2 + w(i) ENDDO ybar = sum1/an wbar = sum2/an sum1 = 0.0_wp sum2 = 0.0_wp sum3 = 0.0_wp DO i = 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) ENDDO sy = SQRT(sum1/(an-1.0_wp)) cc = sum2/SQRT(sum3*sum1) yslope = sum2/sum3 yint = ybar - yslope*wbar corr(idis) = cc yi(idis) = yint ys(idis) = yslope ENDDO ! ! DETERMINE THAT DISTRIBUTION WITH THE MAX PROB PLOT CORR COEFFICIENT ! idismx = 1 corrmx = corr(1) DO idis = 1 , numdis IF ( corr(idis)>corrmx ) idismx = idis IF ( corr(idis)>corrmx ) corrmx = corr(idis) ENDDO DO idis = 1 , numdis iflag1(idis) = blank iflag2(idis) = blank iflag3(idis) = blank IF ( idis==idismx ) THEN iflag1(idis) = alpham iflag2(idis) = alphaa iflag3(idis) = alphax ENDIF ENDDO ! ! WRITE OUT THE TABLE OF PROB PLOT CORR COEFFICIENTS FOR VARIOUS GAMMA ! WRITE (G_IO,99028) WRITE (G_IO,99005) 99005 FORMAT (' ',40X,'EXTREME VALUE ANALYSIS') WRITE (G_IO,99029) WRITE (G_IO,99006) N 99006 FORMAT (' ',37X,'THE SAMPLE SIZE N = ',I0) WRITE (G_IO,99007) ybar 99007 FORMAT (' ',34X,'THE SAMPLE MEAN = ',F14.7) WRITE (G_IO,99008) sy 99008 FORMAT (' ',28X,'THE SAMPLE STANDARD DEVIATION = ',F14.7) WRITE (G_IO,99009) xmin 99009 FORMAT (' ',32X,'THE SAMPLE MINIMUM = ',F14.7) WRITE (G_IO,99010) xmax 99010 FORMAT (' ',32X,'THE SAMPLE MAXIMUM = ',F14.7) WRITE (G_IO,99029) WRITE (G_IO,99011) 99011 FORMAT (' ', & &' EXTREME VALUE PROBABILITY PLOT LOCATION SCA& &LE TAIL LENGTH') WRITE (G_IO,99012) 99012 FORMAT (' ', & &' TYPE 2 TAIL LENGTH CORRELATION ESTIMATE ESTI& &MATE MEASURE') WRITE (G_IO,99013) 99013 FORMAT (' ',' PARAMETER (GAMMA) COEFFICIENT') WRITE (G_IO,99029) ! numdm1 = numdis - 1 IF ( numdm1>=1 ) THEN DO i = 1 , numdm1 WRITE (G_IO,99014) gamtab(i) , corr(i) , iflag1(i) , & & iflag2(i) , iflag3(i) , yi(i) , ys(i) ,& & t(i) 99014 FORMAT (' ',3X,F10.2,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7, & & 3X,F10.5) ENDDO ENDIF i = numdis WRITE (G_IO,99015) alphai , alphan , alphaf , alphai , alphan , & & alphai , alphat , alphay , corr(i) , & & iflag1(i) , iflag2(i) , iflag3(i) , yi(i) , & & ys(i) , t(i) 99015 FORMAT (' ',5X,8A1,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5) ! ! PLOT THE PROB PLOT CORR COEFFICIENT VERSUS GAMMA VALUE INDEX ! CALL PLOT(corr,aindex,numdis) WRITE (G_IO,99016) alphag , alphaa , alpham , alpham , alphaa , & & equal , gamtab(1) , gamtab(12) , gamtab(23) ,& & gamtab(34) , alphai , alphan , alphaf , & & alphai , alphan , alphai , alphat , alphay 99016 FORMAT (' ',12X,5A1,1X,A1,F14.7,11X,F14.7,11X,F14.7,11X,F14.7, & & 15X,8A1) WRITE (G_IO,99029) WRITE (G_IO,99017) 99017 FORMAT (' ', & &'THE ABOVE IS A PLOT OF THE 46 PROBABILITY PLOT CORRELATION COEFFI& &CIENTS (FROM THE PREVIOUS page)') WRITE (G_IO,99018) 99018 FORMAT (' ',16X,'VERSUS THE 46 EXTREME VALUE DISTRIBUTIONS') ! ! IF THE OPTIMAL GAMMA IS FINITE, PLOT OUT THE EXTREME VALUE ! TYPE 2 PROBABILITY PLOT FOR THE OPTIMAL VALUE ! OF GAMMA. ! IF ( idismx<numdis ) CALL EV2PLT(X,N,gamtab(idismx)) ! ! PLOT OUT AN EXTREME VALUE TYPE 1 PROBABILITY PLOT ! CALL EV1PLT(X,N) ! ! FORM THE VARIOUS RETURN PERIOD VALUES ! k = 0 DO i = 1 , 4 DO j = 1 , 9 k = k + 1 am(k) = j*(10**(i-1)) ENDDO ENDDO k = k + 1 am(k) = 10000.0_wp k = k + 1 am(k) = 50000.0_wp k = k + 1 am(k) = 100000.0_wp k = k + 1 am(k) = 500000.0_wp k = k + 1 am(k) = 1000000.0_wp k = k + 1 am(k) = N numam = k CALL SORT(am,numam,scrat) DO i = 1 , numam am(i) = scrat(i) ENDDO ! ! IF THE OPTIMAL GAMMA IS FINITE, COMPUTE THE ! PREDICTED EXTREME (= F(1-(1/M)) FOR VARIOUS RETURN PERIODS M ! FOR THE OPTIMAL EXTREME VALUE TYPE 2 DISTRIBUTION. ! IF ( idismx/=numdis ) THEN a = gamtab(idismx) yint = yi(idismx) yslope = ys(idismx) DO i = 2 , numam r = 1.0_wp/am(i) p = 1.0_wp - r arg = -LOG(p) IF ( arg>0.0_wp ) h(i,1) = yint + yslope*(arg**(-1.0_wp/a)) ENDDO ENDIF ! ! COMPUTE THE PREDICTED EXTREME (= F(1-(1/M)) FOR VARIOUS RETURN ! PERIODS M FOR THE EXTREME VALUE TYPE 1 DISTRIBUTION. ! yint = yi(numdis) yslope = ys(numdis) DO i = 2 , numam r = 1.0_wp/am(i) p = 1.0_wp - r arg = -LOG(p) IF ( arg>0.0_wp ) h(i,2) = yint + yslope*(-LOG(arg)) ENDDO ! ! WRITE OUT THE page WITH THE RETURN PERIODS AND THE PREDICTED EXTREMES ! FOR THE 2 DISTRIBUTIONS--OPTIMAL EXTREME VALUE TYPE 2, AND EXTREME ! VALUE TYPE 1. ! WRITE (G_IO,99028) IF ( idismx==numdis ) THEN ! WRITE (G_IO,99019) 99019 FORMAT (' ',' RETURN PERIOD PREDICTED EXTREME WIND') WRITE (G_IO,99020) 99020 FORMAT (' ',' (IN YEARS) BASED ON') WRITE (G_IO,99021) 99021 FORMAT (' ',' EXTREME VALUE TYPE 1') WRITE (G_IO,99022) 99022 FORMAT (' ',' DISTRIBUTION') WRITE (G_IO,99029) DO i = 2 , numam WRITE (G_IO,99030) am(i) , h(i,2) j = i - 1 jskip = j - 5*(j/5) IF ( jskip==0 ) WRITE (G_IO,99029) ENDDO GOTO 99999 ENDIF ENDIF WRITE (G_IO,99023) 99023 FORMAT (' ',' RETURN PERIOD PREDICTED EXTREME WIND', & & ' PREDICTED EXTREME WIND') WRITE (G_IO,99024) 99024 FORMAT (' ',' (IN YEARS) BASED ON OPTIMAL ', & & ' BASED ON') WRITE (G_IO,99025) 99025 FORMAT (' ',' EXTREME VALUE TYPE 2', & & ' EXTREME VALUE TYPE 1') WRITE (G_IO,99026) 99026 FORMAT (' ',' DISTRIBUTION ', & & ' DISTRIBUTION') WRITE (G_IO,99027) gamtab(idismx) 99027 FORMAT (' ',' (GAMMA = ',F12.5,')') WRITE (G_IO,99029) DO i = 2 , numam WRITE (G_IO,99030) am(i) , h(i,1) , h(i,2) j = i - 1 jskip = j - 5*(j/5) IF ( jskip==0 ) WRITE (G_IO,99029) ENDDO RETURN ! 99028 FORMAT ('1') 99029 FORMAT (' ') 99030 FORMAT (' ',2X,F9.1,13X,F10.2,17X,F10.2) ! 99999 END SUBROUTINE EXTREM