SUBROUTINE TIME(X,N)
C
C     PURPOSE--THIS SUBROUTINE PERFORMS A TIME SERIES ANALYSIS
C              ON THE DATA IN THE INPUT VECTOR X. 
C              THE ANALYSIS CONSISTS OF THE FOLLOWING--
C              1) A PLOT OF AUTOCORRELATION VERSUS LAG NUMBER;
C              2) A TEST FOR WHITE NOISE (ASSUMING NORMALITY);
C              3) A 'PILOT' SPECTRUM; AND
C              4) 4 OTHER ESTIMATED SPECTRA--EACH BASED
C                 ON A DIFFERING BANDWIDTH.
C
C              IN ORDER THAT THE RESULTS OF THE TIME SERIES ANALYSIS
C              BE VALID AND PROPERLY INTERPRETED, THE INPUT DATA
C              IN X SHOULD BE EQUI-SPACED IN TIME 
C              (OR WHATEVER VARIABLE CORRESPONDS TO TIME).
C
C              THE HORIZONTAL AXIS OF THE SPECTRA PRODUCED
C              BY THIS SUBROUTINE IS FREQUENCY.
C              THIS FREQUENCY IS MEASURED IN UNITS OF
C              CYCLES PER 'DATA POINT' OR, MORE PRECISELY, IN
C              CYCLES PER UNIT TIME WHERE
C              'UNIT TIME' IS DEFINED AS THE
C              ELAPSED TIME BETWEEN ADJACENT OBSERVATIONS.
C              THE RANGE OF THE FREQUENCY AXIS IS 0.0 TO 0.5.
C
C     INPUT ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                               (UNSORTED) OBSERVATIONS.
C                      N      = THE INTEGER NUMBER OF OBSERVATIONS
C                               IN THE VECTOR X.
C     OUTPUT--7 TO 11 PAGES (DEPENDING ON
C             THE INPUT SAMPLE SIZE) OF 
C             AUTOMATIC PRINTOUT--
C             1) A PLOT OF AUTOCORRELATION VERSUS LAG NUMBER;
C                THIS PLOT MAY TAKE AS LITTLE AS 1
C                OR AS MANY AS 5 PAGES
C                (THE EXACT NUMBER DEPENDING ON
C                THE INPUT SAMPLE SIZE N);
C             2) A TEST FOR WHITE NOISE (ASSUMING NORMALITY);
C             3) A 'PILOT' SPECTRUM; AND
C             4) AN ESTIMATED SPECTRUM BASED ON A 
C                BANDWIDTH DERIVED FROM THE DATA SET;
C             5) AN ESTIMATED SPECTRUM BASED ON A 
C                BANDWIDTH ONLY 1/2 AS WIDE AS IN 4;
C             6) AN ESTIMATED SPECTRUM BASED ON A 
C                BANDWIDTH ONLY 1/4 AS WIDE AS IN 4;
C             7) AN ESTIMATED SPECTRUM BASED ON A 
C                BANDWIDTH ONLY 1/8 AS WIDE AS IN 4;
C     PRINTING--YES.
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C                 --THE SAMPLE SIZE N MUST BE GREATER
C                   THAN OR EQUAL TO 3. 
C     OTHER DATAPAC   SUBROUTINES NEEDED--PLOTC0, PLOTSP, AND CHSPPF. 
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--THE 'FAST FOURIER TRANSFORM' IS NOT USED
C              IN THIS VERSION OF TIME, BUT WILL BE
C              IMPLEMENTED IN A FUTURE VERSION.
C            --THE USUAL MAXIMUM NUMBER OF LAGS
C              FOR WHICH THE AUTOCORRELATION IS
C              COMPUTED IS N/4 WHERE N IS
C              THE SAMPLE SIZE (LENGTH OF THE
C              DATA RECORD IN THE VECTOR X).
C              THIS RULE IS OVERRIDDEN IN
C              LARGE DATA SETS AND IS REPLACED
C              BY THE RULE THAT THE MAXIMUM
C              NUMBER OF LAGS = 500
C              (WHICH CORRESPONDS TO AN 
C              AUTOCORRELATION PLOT COVERING
C              5 COMPUTER PAGES).
C              IF MORE LAGS ARE DESIRED,
C              CHANGE THE VALUE OF THE
C              VARIABLE     MAXLAG
C              WITHIN THIS SUBROUTINE
C              FROM 500 TO WHATEVER DESIRED,
C              AND ALSO CHANGE THE DIMENSION OF
C              THE VECTOR R FROM ITS PRESENT 500 TO HOWEVER 
C              MANY LAGS ARE DESIRED.
C            --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED
C              TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME,
C              THEN THE FREQUENCY AXIS OF THE RESULTING
C              SPECTRA WOULD BE IN UNITS OF HERTZ 
C              (= CYCLES PER SECOND).
C            --THE FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE
C              IN THE DATA OF INFINITE (= 1/(0.0))
C              LENGTH OR PERIOD.
C              THE FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE
C              IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS.
C            --ANY EQUI-SPACED TIME SERIES ANALYSIS IS
C              INTRINSICALLY LIMITED TO DETECTING FREQUENCIES
C              NO LARGER THAN 0.5 CYCLES PER DATA POINT;
C              THIS CORRESPONDS TO THE FACT THAT THE
C              SMALLEST DETECTABLE CYCLE IN THE DATA
C              IS 2 DATA POINTS PER CYCLE.
C     REFERENCES--JENKINS AND WATTS, ESPECIALLY PAGE 290.
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     UPDATED         --FEBRUARY  1977. 
C
C---------------------------------------------------------------------
C
      DIMENSION X(1)
      DIMENSION R(500)
      DIMENSION S(125)
      DIMENSION PSSQ(6),PMSQ(6),PS(6),P(5),L(4)
      DATA PI/3.14159265358979/
C
      IPR=6
      ILOWER=3
      MAXLAG=500
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.ILOWER)GOTO50
      HOLD=X(1)
      DO65I=2,N
      IF(X(I).NE.HOLD)GOTO90
   65 CONTINUE
      WRITE(IPR, 9)HOLD
      RETURN
   50 WRITE(IPR,17)ILOWER
      WRITE(IPR,47)N
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE TIME   SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 96H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 TIME   SUBROUTINE IS OUTSIDE THE ALLOWABLE (,I6,11H,INFINITY) ,
     1 14HINTERVAL *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      AN=N
C
C     COMPUTE THE SAMPLE MEAN 
C
      SUM=0.0
      DO100I=1,N
      SUM=SUM+X(I)
  100 CONTINUE
      XBAR=SUM/AN
C
C     COMPUTE THE SAMPLE VARIANCE AND THE SUM OF SQUARED DEVIATIONS
C
      SUM=0.0
      DO200I=1,N
      SUM=SUM+(X(I)-XBAR)*(X(I)-XBAR)
  200 CONTINUE
      SSQ=SUM
      VARB=SSQ/AN
      VAR=SSQ/(AN-1.0)
      SD=SQRT(VAR)
C
C     COMPUTE THE SAMPLE AUTOCORRELATIONS
C     REFERENCE--JENKINS AND WATTS, PAGES 290 AND 259 (7.1.6)
C
      KMAX=N/4
      IF(N.LE.32)KMAX=N/2
      IF(N.LE.16)KMAX=N
      IF(KMAX.GT.MAXLAG)KMAX=MAXLAG
      DO300K=1,KMAX 
      SUM=0.0
      NMK=N-K
      DO400I=1,NMK
      J=I+K
      SUM=SUM+(X(I)-XBAR)*(X(J)-XBAR)
  400 CONTINUE
      R(K)=SUM/SSQ
  300 CONTINUE
C
C     PLOT THE SAMPLE AUTOCORRELATIONS
C
      CALL PLOTCO(R,KMAX)
      WRITE(IPR,1705)
      WRITE(IPR,1710)N,KMAX
C
C     DO A WHITE NOISE ANALYSIS
C
      SDR=1.0/SQRT(AN)
      R975=1.96*SDR 
      IF(R975.GT.1.0)R975=1.0 
      R025=-R975
      NUMOUT=0
      DO410K=1,KMAX 
      ABSR=R(K)
      IF(ABSR.LT.0.0)ABSR=-ABSR
      IF(ABSR.GT.R975)NUMOUT=NUMOUT+1
  410 CONTINUE
      PEROUT=FLOAT(NUMOUT)/FLOAT(KMAX)
      PEROUT=100.0*PEROUT
      WRITE(IPR,999)
      WRITE(IPR,420)R025,R975 
      WRITE(IPR,430)
      WRITE(IPR,440)NUMOUT,KMAX,PEROUT
      DO415I=1,5
      WRITE(IPR,999)
  415 CONTINUE
      WRITE(IPR,1715)N
      WRITE(IPR,1720)XBAR
      WRITE(IPR,1725)VAR
      WRITE(IPR,1730)SD
      WRITE(IPR,1735)VARB
C
C     COMPUTE THE PILOT SPECTRUM FOR THE REDUCED (2**J) SAMPLE
C     REFERENCE--JENKINS AND WATTS, PAGE 288
C
      DO450I=1,20
      NDIV=N/(2**I) 
      IF(NDIV.EQ.0)I2=I-1
      IF(NDIV.EQ.0)GOTO460
  450 CONTINUE
  460 IF(7.LT.I2)I2=7
      N2=2**I2
      AN2=N2
      DO500K=1,I2
      SUM=0.0
      IMIN=2**K
      JMAX=IMIN/2
      DO600I=IMIN,N2,IMIN
      SUM1=0.0
      SUM2=0.0
      DO700J=1,JMAX 
      IARG1=I+J-JMAX
      IARG2=IARG1-JMAX
      SUM1=SUM1+X(IARG1)
      SUM2=SUM2+X(IARG2)
  700 CONTINUE
      SUM=SUM+(SUM1-SUM2)*(SUM1-SUM2)
  600 CONTINUE
      PSSQ(K)=SUM/FLOAT(IMIN) 
      PMSQ(K)=PSSQ(K)/AN2
      PS(K)=FLOAT(2*IMIN)*PMSQ(K)
      PS(K)=PS(K)/VARB
  500 CONTINUE
C
C     FORM THE PILOT SPECTRUM PLOT
C
      DO900I=1,I2
      IREV=I2-I+1
      JMIN=(120/(2**I))+1
      IF(I.EQ.I2)JMIN=1
      JMAX=120/(2**(I-1))
      DO1000J=JMIN,JMAX
      S(J)=PS(I)
 1000 CONTINUE
  900 CONTINUE
      CALL PLOTSP(S,120,0)
      WRITE(IPR,999)
      WRITE(IPR,1750)
C
C     DEFINE 4 LAG WINDOW TRUNCATION POINTS
C     REFERENCE--JENKINS AND WATTS, PAGES 290 AND 260
C
      P(1)=.2
      P(2)=.1
      P(3)=.05
      P(4)=.025
      P(5)=.01
      LMAX=0
      DO1100I=1,5
      DO1200K=1,KMAX
      KREV=KMAX-K+1 
      RK=R(KREV)
      IF(RK.LT.0.0)RK=-RK
      IF(RK.GE.P(I))LMAX=KREV 
      IF(RK.GE.P(I))GOTO1350
 1200 CONTINUE
 1100 CONTINUE
      IF(LMAX.NE.0)GOTO1350
      RMAX=ABS(R(1))
      DO1300K=1,KMAX
      RK=R(K)
      IF(RK.LT.0.0)RK=-RK
      IF(RK.GE.RMAX)LMAX=K
      IF(RK.GE.RMAX)RMAX=RK
 1300 CONTINUE
 1350 ALMAX=LMAX
      L(1)=(3.0/2.0)*ALMAX
      IF(L(1).LE.32)LMAX=32
      IF(L(1).LE.32)ALMAX=32.0
      IF(L(1).LE.32)L(1)=32
      IF(L(1).GE.KMAX)LMAX=KMAX
      IF(L(1).GE.KMAX)ALMAX=KMAX
      IF(L(1).GE.KMAX)L(1)=KMAX
      L(2)=(ALMAX/2.0)+0.1
      L(3)=(ALMAX/4.0)+0.1
      L(4)=(ALMAX/8.0)+0.1
      IF(L(4).GE.3)NUMSP=4
      IF(L(4).GE.3)GOTO1380
      IF(L(3).GE.3)NUMSP=3
      IF(L(3).GE.3)GOTO1380
      IF(L(2).GE.3)NUMSP=2
      IF(L(2).GE.3)GOTO1380
      IF(L(1).GE.3)NUMSP=1
      IF(L(1).GE.3)GOTO1380
      WRITE(IPR,1390)N
      RETURN
 1380 CONTINUE
C
C     COMPUTE THE 4 SPECTRUM ESTIMATES
C     REFERENCE--JENKINS AND WATTS, PAGES 260 AND 244
C
C     COMPUTE BANDWIDTHS
C     REFERENCE--JENKINS AND WATTS, PAGES 257 AND 252
C
C     COMPUTE DEGREES OF FREEDOM FOR THE SPECTAL DENSITY ESTIMATE AT INDIVIDUAL 
C     FREQUENCIES
C     REFERENCE--JENKINS AND WATTS, PAGES 254 AND 252
C
C     COMPUTE 95 PERCENT CONFIDENCE INTERVAL LENGTHS FOR THE LOG SPECTRAL
C     DENSITY ESTIMATES
C     REFERENCE--JENKINS AND WATTS, PAGES 255 AND 252
C
C     WRITE OUT THE 4 SPECTRUM PLOTS
C
      DO1400I=1,NUMSP
      AL=L(I)
      LM1=L(I)-1
      DO1500LLP1=1,121
      LL=LLP1-1
      ALL=LL
      SUM=0.0
      DO1600K=1,LM1 
      AK=K
      ARG1=PI*AK/AL 
      ARG2=PI*AK*ALL/120.0
      WK=0.0
      IF(K.LE.L(I))WK=0.5*(1.0+COS(ARG1))
      SUM=SUM+R(K)*WK*COS(ARG2)
 1600 CONTINUE
      SUM=2.0+4.0*SUM
      S(LLP1)=SUM
 1500 CONTINUE
      BW=(4.0/3.0)/FLOAT(L(I))
      DF=(8.0/3.0)*AN/FLOAT(L(I))
      IDF=DF+0.5
      CALL PLOTSP(S,121,IDF)
      DFROUN=IDF
      WRITE(IPR,999)
      WRITE(IPR,1760)L(I),BW,DFROUN
 1400 CONTINUE
C
  420 FORMAT(1H ,103HUNDER THE NULL HYPOTHESIS OF WHITE NOISE (AND NORMA
     1LITY), A 2-SIDED 95 PERCENT ACCEPTANCE INTERVAL IS (,F6.4,1H,,F6.4
     1,1H))
  430 FORMAT(1H ,125HUNDER THE NULL HYPOTHESIS, ONLY 5 PERCENT (ON THE A
     1VERAGE) OF THE OBSERVED AUTOCORRELATIONS SHOULD FALL OUTSIDE THIS
     1INTERVAL)
  440 FORMAT(1H ,20HIT IS OBSERVED THAT ,I5,12H OUT OF THE ,I5,11H (THAT
     1 IS, ,F5.1,72H PERCENT) OF THE COMPUTED AUTOCORRELATIONS FALL OUTS
     1IDE OF THIS INTERVAL)
 1390 FORMAT(1H ,45HDUE TO THE SMALL NUMBER OF OBSERVATIONS (N = ,I6,59H
     1), THERE ARE NOT ENOUGH LAGS TO PRODUCE A RELIABLE SPECTRUM)
 1705 FORMAT(1H ,30X, 63HAUTOCORRELATION PLOT--PLOT OF SAMPLE AUTOCORREL
     1ATION VERSUS LAG)
 1710 FORMAT(1H ,10X,29HTHE NUMBER OF OBSERVATIONS = ,I6,10X,56HTHE NUMB
     1ER OF COMPUTED (AND PLOTTED) AUTOCORRELATIONS = ,I6)
 1715 FORMAT(1H ,18HTHE SAMPLE SIZE = ,I6)
 1720 FORMAT(1H ,18HTHE SAMPLE MEAN = ,E15.8)
 1725 FORMAT(1H ,22HTHE SAMPLE VARIANCE = ,E15.8) 
 1730 FORMAT(1H ,32HTHE SAMPLE STANDARD DEVIATION = ,E15.8) 
 1735 FORMAT(1H ,29HTHE BIASED SAMPLE VARIANCE = ,E15.8)
 1750 FORMAT(1H ,50X,14HPILOT SPECTRUM) 
 1760 FORMAT(1H ,17HNUMBER OF LAGS = ,I5,10X,11HBANDWIDTH =,F10.3,10X,21
     1HDEGREES OF FREEDOM = ,F10.3)
  999 FORMAT(1H )
C
      RETURN
      END