      SUBROUTINE DCFFTI(N,WSAVE)
C***BEGIN PROLOGUE  DCFFTI
C***DATE WRITTEN   790601   (YYMMDD)
C***REVISION DATE  860115   (YYMMDD)
C***CATEGORY NO.  J1A2
C***KEYWORDS  FOURIER TRANSFORM
C***AUTHOR  SWARZTRAUBER, P. N., (NCAR)
C***PURPOSE  Initialize for DCFFTF and DCFFTB.
C***DESCRIPTION
C           From the book, "Numerical Methods and Software" by
C                D. Kahaner, C. Moler, S. Nash
C                Prentice Hall, 1988
C
C  Subroutine DCFFTI initializes the array WSAVE which is used in
C  both DCFFTF and DCFFTB.  The prime factorization of N together with
C  a tabulation of the trigonometric functions are computed and
C  stored in WSAVE.
C
C  Input Parameter
C
C  N       the length of the sequence to be transformed
C
C  Output Parameter
C
C  WSAVE   a work array which must be dimensioned at least 4*N+15.
C          The same work array can be used for both DCFFTF and DCFFTB
C          as long as N remains unchanged.  Different WSAVE arrays
C          are required for different values of N.  The contents of
C          WSAVE must not be changed between calls of DCFFTF or DCFFTB.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DCFTI1
C***END PROLOGUE  DCFFTI
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WSAVE(*)
C***FIRST EXECUTABLE STATEMENT  DCFFTI
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL DCFTI1 (N,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE DCFTI1(N,WA,IFAC)
C***BEGIN PROLOGUE  DCFTI1
C***REFER TO  DCFFTI
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  DCFTI1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WA(*)      ,IFAC(*)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
C***FIRST EXECUTABLE STATEMENT  DCFTI1
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J-4) 102,102,103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR) 101,105,101
  105 NF = NF+1
      IFAC(NF+2) = NTRY
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         IFAC(IB+2) = IFAC(IB+1)
  106 CONTINUE
      IFAC(3) = 2
  107 IF (NL .NE. 1) GO TO 104
      IFAC(1) = N
      IFAC(2) = NF
      TPI = 8.0D0*ATAN(1.0D0)
      ARGH = TPI/DBLE(N)
      I = 2
      L1 = 1
      DO 110 K1=1,NF
         IP = IFAC(K1+2)
         LD = 0
         L2 = L1*IP
         IDO = N/L2
         IDOT = IDO+IDO+2
         IPM = IP-1
         DO 109 J=1,IPM
            I1 = I
            WA(I-1) = 1.0D0
            WA(I) = 0.0D0
            LD = LD+L1
            FI = 0.0D0
            ARGLD = DBLE(LD)*ARGH
            DO 108 II=4,IDOT,2
               I = I+2
               FI = FI+1.0D0
               ARG = FI*ARGLD
               WA(I-1) = COS(ARG)
               WA(I) = SIN(ARG)
  108       CONTINUE
            IF (IP .LE. 5) GO TO 109
            WA(I1-1) = WA(I-1)
            WA(I1) = WA(I)
  109    CONTINUE
         L1 = L2
  110 CONTINUE
      RETURN
      END
