      SUBROUTINE DEZFTI(N,WSAVE)   
C***BEGIN PROLOGUE  DEZFTI    
C***DATE WRITTEN   790601   (YYMMDD)    
C***REVISION DATE  860115   (YYMMDD)    
C***CATEGORY NO.  J1A1   
C***KEYWORDS  FOURIER TRANSFORM    
C***AUTHOR  SWARZTRAUBER, P. N., (NCAR) 
C***PURPOSE  Initialize DEZFTF and DEZFTB    
C***DESCRIPTION
C    
C  Subroutine DEZFTI initializes the array WSAVE which is used in
C  both DEZFTF and DEZFTB.  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 3*N+15.    
C          The same work array can be used for both DEZFTF and DEZFTB 
C          as long as N remains unchanged.  Different WSAVE arrays    
C          are required for different values of N.
C***REFERENCES  (NONE)   
C***ROUTINES CALLED  DEZFT1   
C***END PROLOGUE  DEZFTI 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WSAVE(*)
C***FIRST EXECUTABLE STATEMENT  DEZFTI  
      IF (N .EQ. 1) RETURN    
      CALL DEZFT1 (N,WSAVE(2*N+1),WSAVE(3*N+1))   
      RETURN   
      END 
      SUBROUTINE DEZFT1(N,WA,IFAC) 
C***BEGIN PROLOGUE  DEZFT1    
C***REFER TO  DEZFTI
C***ROUTINES CALLED  (NONE)   
C***END PROLOGUE  DEZFT1 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WA(*)      ,IFAC(*)    ,NTRYH(4) 
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
C***FIRST EXECUTABLE STATEMENT  DEZFT1  
      TPI = 8.0D0*ATAN(1.0D0)  
      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  
      ARGH = TPI/DBLE(N) 
      IS = 0   
      NFM1 = NF-1   
      L1 = 1   
      IF (NFM1 .EQ. 0) RETURN 
      DO 111 K1=1,NFM1   
         IP = IFAC(K1+2) 
         L2 = L1*IP 
         IDO = N/L2 
         IPM = IP-1 
         ARG1 = DBLE(L1)*ARGH 
         CH1 = 1.0D0   
         SH1 = 0.0D0   
         DCH1 = COS(ARG1)
         DSH1 = SIN(ARG1)
         DO 110 J=1,IPM  
            CH1H = DCH1*CH1-DSH1*SH1    
            SH1 = DCH1*SH1+DSH1*CH1
            CH1 = CH1H   
            I = IS+2
            WA(I-1) = CH1
            WA(I) = SH1  
            IF (IDO .LT. 5) GO TO 109   
            DO 108 II=5,IDO,2 
               I = I+2   
               WA(I-1) = CH1*WA(I-3)-SH1*WA(I-2)  
               WA(I) = CH1*WA(I-2)+SH1*WA(I-3)    
  108       CONTINUE
  109       IS = IS+IDO  
  110    CONTINUE   
         L1 = L2    
  111 CONTINUE 
      RETURN   
      END
      SUBROUTINE DRFFTI(N,WSAVE)    
C***BEGIN PROLOGUE  DRFFTI
C***DATE WRITTEN   790601   (YYMMDD)    
C***REVISION DATE  860115   (YYMMDD)    
C***CATEGORY NO.  J1A1   
C***KEYWORDS  FOURIER TRANSFORM    
C***AUTHOR  SWARZTRAUBER, P. N., (NCAR) 
C***PURPOSE  Initialize for DRFFTF and DRFFTB. 
C***DESCRIPTION
C    
C           From the book, "Numerical Methods and Software" by
C                D. Kahaner, C. Moler, S. Nash
C                Prentice Hall, 1989
C    
C  Subroutine DRFFTI initializes the array WSAVE which is used in 
C  both DRFFTF and DRFFTB.  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 DP work array which must be dimensioned at least 2*N+15.    
C          The same work array can be used for both DRFFTF and DRFFTB   
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 DRFFTF or DRFFTB. 
C    
C  *********************************************************************   
C  *                                                                   *   
C  *   SUBPROGRAM REVISION HISTORY                                     *   
C  *                                                                   *   
C  *   06/01/79  -  Original version by Paul Swarztrauber.             *   
C  *                Distributed by NCAR (ref. 1).                      *   
C  *   04/01/83  -  SLATEC Common Math Library Subcommittee.           *   
C  *                Modified to use SLATEC library source file format. *   
C  *                Distributed in the SLATEC library (ref. 2).        *   
C  *   01/15/86  -  Ron Boisvert, National Bureau of Standards.        *   
C  *                Modified to convert to portable Fortran 77.        *   
C  *                                                                   *   
C  *   The changes introduced in the most recent modification are      *   
C  *                                                                   *   
C  *   (a) Dummy array size declarations (1) changed to (*)            *   
C  *   (b) References to intrinsic function FLOAT changed to DP        *   
C  *   (c) Mathematical constants previously coded in DATA state-      *   
C  *       ments now computed at runtime using Fortran intrinsic       *   
C  *       functions.  The affected variables are                      *   
C  *                                                                   *   
C  *          PI      SQRT2   SQRT3   TAUR    TR11    TR12             *   
C  *          PIH     TSQRT2          TAUI    TI11    TI12             *   
C  *          TPI     HSQT2                                            *   
C  *                                                                   *   
C  *   References                                                      *   
C  *                                                                   *   
C  *   1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel         *   
C  *      Computations (G. Rodrigue, ed.), Academic Press, 1982,       *   
C  *      pp. 51-83.                                                   *   
C  *   2. B.L. Buzbee, The SLATEC Common Math Library, in Sources      *   
C  *      and Development of Mathematical Software (W. Cowell, ed.),   *   
C  *      Prentice-Hall, 1984, pp. 302-318.                            *   
C  *                                                                   *   
C  *********************************************************************   
C    
C***REFERENCES  (NONE)   
C***ROUTINES CALLED  DRFTI1   
C***END PROLOGUE  DRFFTI  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WSAVE(*)
C***FIRST EXECUTABLE STATEMENT  DRFFTI   
      IF (N .EQ. 1) RETURN    
      CALL DRFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))
      RETURN   
      END 
      SUBROUTINE DRFTI1(N,WA,IFAC) 
C***BEGIN PROLOGUE  DRFTI1    
C***REFER TO  DRFFTI 
C***ROUTINES CALLED  (NONE)   
C***END PROLOGUE  DRFTI1 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WA(*)      ,IFAC(*)    ,NTRYH(4) 
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
C***FIRST EXECUTABLE STATEMENT  DRFTI1  
      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.D0*ATAN(1.D0)  
      ARGH = TPI/N 
      IS = 0   
      NFM1 = NF-1   
      L1 = 1   
      IF (NFM1 .EQ. 0) RETURN 
      DO 110 K1=1,NFM1   
         IP = IFAC(K1+2) 
         LD = 0
         L2 = L1*IP 
         IDO = N/L2 
         IPM = IP-1 
         DO 109 J=1,IPM  
            LD = LD+L1   
            I = IS  
            ARGLD = LD*ARGH  
            FI = 0.D0 
            DO 108 II=3,IDO,2 
               I = I+2   
               FI = FI+1.D0
               ARG = FI*ARGLD 
               WA(I-1) = COS(ARG)  
               WA(I) = SIN(ARG)    
  108       CONTINUE
            IS = IS+IDO  
  109    CONTINUE   
         L1 = L2    
  110 CONTINUE 
      RETURN   
      END
