C Find autocorrelation to el Nino data using direct and FFT methods.
C
      INTEGER N
      PARAMETER(N=168)
      REAL EL(0:2*N-1),WSAVE(4*(2*N)+15),ACOV(0:N-1),A(2*N),B(2*N),
     *           ACOVR(0:N-1)
      COMPLEX CEL(0:2*N-1),CORR(0:2*N-1)
      LOGICAL EX
C
C needs N locations for el, 2N for acov, cel, corr
C and 4*(2N)+15 for wsave in complex case. N has maximum of 168.
      INQUIRE (FILE='ELNINO.DAT',EXIST=EX,ERR=1000)
      IF(.NOT.EX)GOTO 1000
      OPEN (UNIT=8,FILE='ELNINO.DAT',ERR=1000)
C
C Read data, find mean.
      SUM = 0.0
      DO 100 I = 0,N-1
         READ (8,*) EL(I)
         SUM = SUM + EL(I)
  100 CONTINUE
      DO 101 I = 0,N-1
         EL(I) = EL(I) - SUM/N
C Subtract mean, and add N zeros for either complex or real FFT usage
         CEL(I) = CMPLX(EL(I),0.0)
         EL(I+N) = 0.0
         CEL(I+N) = 0.0
  101 CONTINUE
C
C----------------------------------------------------------------------
C
C Direct calculation. Only sum as far as there is data. 
C Simple, but slow.
C      
      DO 110 J = 0,N-1
         ACOV(J) = 0.0
         DO 110 M = 0,N-1-J
            ACOV(J) = ACOV(J)+ EL(M) * EL(M+J)
  110 CONTINUE
C
C Write, scaled correlation
      WRITE (*,*) 
     * 'EX 11.6: AUTOCORRELATION (DIRECT) OUTPUT SUPRESSED'
CCCCC      WRITE (*,'(5E14.6)') ( (ACOV(I)/ACOV(0)), I = 0,N-1)
C----------------------------------------------------------------------
C
C FFT approach (complex).
C    Compute FFT of data of length 2N.
C    Compute square of magnitude of transform components and place  
C       in complex array as real parts.
C    Compute inverse transform, throwing away second half and
C       imaginary parts (which are zero), and multiply by length of 
C       sequence, 2N.  
      CALL CFFTI (2*N,WSAVE)
      CALL CFFTF (2*N,CEL,WSAVE)
C CFFTF returns unscaled transforms. Actual transforms are output
C divided by (2N).
      DO 120 I = 0,2*N-1
         CORR(I) = ABS(CEL(I) / (2*N)) **2 
  120 CONTINUE
C Since we compute transform times its conjugate, must divide by
C (2N) for each, i.e., (2N)**2.
      CALL CFFTB (2*N,CORR,WSAVE)
C
      DO 121 I = 0,N-1
         ACOV(I) = REAL(CORR(I))*(2*N)
  121 CONTINUE
C Autocovariance is inverse transform times sequence length, 2N.
C Normally, all the scaling would be done  only once
C by dividing by 2N. We've broken it up for exposition.
      WRITE (*,*) 
     *  'EX 11.6: AUTOCORRELATION (COMPLEX FFT) OUTPUT SUPRESSED'
CCCCC      WRITE (*,'(5E14.6)') ( (ACOV(I)/ACOV(0) ), I = 0,N-1)
C
C----------------------------------------------------------------------
C
C FFT approach (real).
C    Compute FFT of data of length 2N.
C    EZFFTF produces correctly scaled A's and B's so no extra scaling
C       needed to get transform.     
C    Compute array of square of each frequency component and place
C       in cosine array (A's) to be back transformed. Set B's to 0.
C    There are N A's, and N B's.
C    Note that care must be taken to compute magnitude correctly, 
C       0.5*(A(I)**2+B(I)**2) for I < N, twice that for I=N.
C    Compute back transform throwing away its second half.
C
      CALL EZFFTI (2*N,WSAVE)
      CALL EZFFTF (2*N,EL,AZERO,A,B,WSAVE)
      AZERO = AZERO*AZERO
C
      DO 150 I = 1,N
         IF(I.NE.N) THEN
            A(I) = (A(I)**2 + B(I)**2) / 2.0
         ELSE
            A(I) = (A(I)**2 + B(I)**2)
         ENDIF
         B(I) = 0.0
  150 CONTINUE
      CALL EZFFTB (2*N,ACOVR,AZERO,A,B,WSAVE)
      WRITE (*,*) 
     * 'EX 11.6: AUTOCORRELATION (REAL FFT) OUTPUT REDUCED'
      WRITE (*,'(5E14.6)') ( (ACOVR(I)/ACOVR(0) ), I = 0,19)
CCCCC      WRITE (*,'(5E14.6)') ( (ACOVR(I)/ACOVR(0) ), I = 0,N-1)
      WRITE (*,*)
      WRITE (*,*) 
     * 'REFERENCE RESULTS (EX 11.6 PARTIAL-LAST 4 LINES) FROM IBM PC/AT'
      WRITE(*,*)' 0.100000E+01  0.606740E+00  0.353857E+00  0.184303E+00 
     * -0.152914E-01'
      WRITE(*,*)'-0.219872E+00 -0.296218E+00 -0.291459E+00 -0.155056E+00 
     *  0.368081E-01'
      WRITE(*,*)' 0.173559E+00  0.266598E+00  0.304989E+00  0.201169E+00 
     *  0.178037E-01'
      WRITE(*,*)'-0.210367E+00 -0.377387E+00 -0.437960E+00 -0.460333E+00 
     * -0.425590E+00'
C
      STOP
 1000 WRITE (*,*) 'CANNOT FIND THE DATA FILE: ELNINO.DAT '
      END