C Find autocorrelation to el Nino data using direct and FFT methods.
C
      REAL EL(0:335),WSAVE(1360),ACOV(0:335),A(336),B(336),ACOVR(0:335)
      COMPLEX CEL(0:335),CORR(0:335)
      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.
      N = 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) '
      WRITE (*,*) ( (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) '
      WRITE (*,*) ( (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       and multiply by length of sequence, 2N.  
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) '
      WRITE (*,*) ( (ACOVR(I)* (2*N)/ACOVR(0) ), I = 0,N-1)
      WRITE (*,*)
      WRITE (*,*) 
     * 'REFERENCE RESULTS (EX 11.6 PARTIAL-LAST 4 LINES) FROM IBM PC/AT'
      WRITE (*,*)'  -0.109820        2.48239        2.09944        4.468
     *69        2.07501'
      WRITE (*,*)'    1.78581       0.224745       0.702122E-01   0.1963
     *69       -1.50187'
      WRITE (*,*)'   0.668857      -0.102865        1.70797        1.289
     *04        1.95623'
      WRITE (*,*)'    1.13909        1.09402       0.570615'
C
      STOP
 1000 WRITE (*,*) 'CANNOT FIND THE DATA FILE: ELNINO.DAT '
      END