C Find autocorrelation to el Nino data using direct and FFT methods.
C
      INTEGER N
      PARAMETER(N=168)
      DOUBLE PRECISION EL(0:2*N-1),WSAVE(4*(2*N)+15),ACOV(0:N-1),
     *           A(2*N),B(2*N),ACOVR(0:2*N-1),SUM,AZERO
      COMPLEX*16 CEL(0:2*N-1),CORR(0:2*N-1)
      LOGICAL EX
C
C  NOTE:  DOUBLE PRECISION COMPLEX IS NOT STANDARD FORTRAN.
C         SOME COMPILER VENDORS MAY FAIL TO IMPLEMENT IT
C         OR MAY IMPLEMENT IN A DIFFERENT MANNER.  IF YOUR
C         COMPILER DOES NOT ACCEPT THE COMPLEX*16 STATEMENT
C         WE RECOMMEND THAT YOU CHECK WITH YOUR VENDOR AND
C         MODIFY THE ABOVE ACCORDINGLY.
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='DELNINO.DAT',EXIST=EX,ERR=1000)
      IF(.NOT.EX)GOTO 1000
      OPEN (UNIT=8,FILE='DELNINO.DAT',ERR=1000)
C
C Read data, find mean.
      SUM = 0.0D0
      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 d.p. FFT usage
         CEL(I) = CMPLX(EL(I),0.0D0)
         EL(I+N) = 0.0D0
         CEL(I+N) = 0.0D0
  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.0D0
         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 (*,'(5D14.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 d.p. 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 DCFFTI (2*N,WSAVE)
      CALL DCFFTF (2*N,CEL,WSAVE)
C DCFFTF 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 DCFFTB (2*N,CORR,WSAVE)
C
      DO 121 I = 0,N-1
         ACOV(I) = DBLE(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 (*,'(5D14.6)') ( (ACOV(I)/ACOV(0) ), I = 0,N-1)
C
C----------------------------------------------------------------------
C
C FFT approach (d.p.).
C    Compute FFT of data of length 2N.
C    DEZFTF 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 DEZFTI (2*N,WSAVE)
      CALL DEZFTF (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.0D0
         ELSE
            A(I) = (A(I)**2 + B(I)**2)
         ENDIF
         B(I) = 0.0D0
  150 CONTINUE
      CALL DEZFTB (2*N,ACOVR,AZERO,A,B,WSAVE)
      WRITE (*,*) 
     * 'EX 11.6: AUTOCORRELATION (DBLE FFT) OUTPUT REDUCED'
      WRITE (*,'(1X,3D18.10)') ( (ACOVR(I)/ACOVR(0) ), I = 0,19)
CCCCC      WRITE (*,'(5D14.6)') ( (ACOVR(I)/ACOVR(0) ), I = 0,N-1)
      WRITE (*,*)
      WRITE (*,*) 
     * 'REFERENCE RESULTS (EX 11.6 PARTIAL-LAST 7 LINES) FROM IBM PC/AT'
      WRITE (*,*)    
     *'  0.1000000000E+01  0.6067398062E+00  0.3538571074E+00'
      WRITE (*,*)
     *'  0.1843027185E+00 -0.1529114172E-01 -0.2198717698E+00'
      WRITE (*,*)
     *' -0.2962181307E+00 -0.2914580887E+00 -0.1550562043E+00'
      WRITE (*,*)
     *'  0.3680803130E-01  0.1735588664E+00  0.2665974779E+00'
      WRITE (*,*)
     *'  0.3049896546E+00  0.2011691615E+00  0.1780362880E-01'
      WRITE (*,*)
     *' -0.2103672514E+00 -0.3773868663E+00 -0.4379599788E+00'
      WRITE (*,*)
     *' -0.4603327597E+00 -0.4255892823E+00'
C
      STOP
 1000 WRITE (*,*) 'CANNOT FIND THE DATA FILE: DELNINO.DAT '
      END