      SUBROUTINE RFFTI (N,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE RFFTI(N,WSAVE)
C
C       ****************************************************************
C
C     SUBROUTINE RFFTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     BOTH RFFTF AND RFFTB. 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. N MUST BE
C             EVEN
C
C     OUTPUT PARAMETER
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 2.5*N+15.
C             THE SAME WORK ARRAY CAN BE USED FOR BOTH RFFTF AND RFFTB
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 RFFTF OR RFFTB.
C
      DIMENSION       WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HRFFTI  ,7HFFTPACK,10HVERSION  1)
      NS2 = N/2
      NQM = (NS2-1)/2
      TPI = 8.*ATAN(1.)
      DT = TPI/FLOAT(N)
      DC = COS(DT)
      DS = SIN(DT)
      WSAVE(1) = DC
      WSAVE(NS2-1) = DS
      IF (NQM .LT. 2) GO TO 102
      DO 101 K=2,NQM
         KC = NS2-K
         WSAVE(K) = DC*WSAVE(K-1)-DS*WSAVE(KC+1)
         WSAVE(KC) = DS*WSAVE(K-1)+DC*WSAVE(KC+1)
  101 CONTINUE
  102 IW1 = NS2+1
      CALL CFFTI (NS2,WSAVE(IW1))
      RETURN
      END
      SUBROUTINE RFFTF (N,R,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE RFFTF(N,R,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE RFFTF COMPUTES THE FOURIER COEFFICIENTS OF A REAL
C     PERODIC SEQUENCE (FOURIER ANALYSIS). THE TRANSFORM IS DEFINED
C     BELOW AT OUTPUT PARAMETER R.
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE ARRAY R. N MUST BE EVEN AND THE METHOD
C             IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES.
C             N MAY CHANGE SO LONG AS DIFFERENT WORK ARRAYS ARE PROVIDED
C
C     R       A REAL ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE
C             TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 2.5*N+15
C             IN THE PROGRAM THAT CALLS RFFTF. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE RFFTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C             THE SAME WSAVE ARRAY CAN BE USED BY RFFTF AND RFFTB.
C
C
C     OUTPUT PARAMETERS
C
C     R       FOR K=2,...,N/2
C
C                  R(2*K-1)= THE SUM FROM I=1 TO I=N OF
C
C                       2.*R(I)*COS((K-1)*(I-1)*2*PI/N)
C
C                  R(2*K)= THE SUM FROM I=1 TO I=N OF
C
C                       2.*R(I)*SIN((K-1)*(I-1)*2*PI/N)
C
C             ALSO
C
C                  R(1)= THE SUM FROM I=1 TO I=N OF 2.*R(I)
C
C                  R(2)= THE SUM FROM I=1 TO I=N OF 2.*(-1)**(I+1)*R(I)
C
C      *****  NOTE
C                  THIS TRANSFORM IS UNNORMALIZED SINCE A CALL OF RFFTF
C                  FOLLOWED BY A CALL OF RFFTB WILL MULTIPLY THE INPUT
C                  SEQUENCE BY 2*N.
C
C     WSAVE   CONTAINS RESULTS WHICH MUST NOT BE DESTROYED BETWEEN
C             CALLS OF RFFTF OR RFFTB.
C
C
      DIMENSION       R(2,1)     ,WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HRFFTF  ,7HFFTPACK,10HVERSION  1)
      IF (N .GT. 2) GO TO 101
      R1 = 2.*(R(1,1)+R(2,1))
      R(2,1) = 2.*(R(1,1)-R(2,1))
      R(1,1) = R1
      RETURN
  101 IW1 = N/2+1
      CALL RFFTF1 (N,R,WSAVE(IW1),WSAVE)
      RETURN
      END
      SUBROUTINE RFFTF1 (N,X,XH,W)
      SAVE
      DIMENSION       X(2,1)     ,XH(2,1)    ,W(1)
      NS2 = N/2
      NS2P2 = NS2+2
      NQ = NS2/2
      IPAR = NS2-NQ-NQ
      NQM = NQ
      IF (IPAR .EQ. 0) NQM = NQM-1
      NQP = NQM+1
      CALL CFFTF (NS2,X,XH)
      IF (NQP .LT. 2) GO TO 107
      DO 101 K=2,NQP
         KC = NS2P2-K
         XH(1,KC) = X(2,K)+X(2,KC)
         XH(2,KC) = X(1,KC)-X(1,K)
  101 CONTINUE
      DO 102 K=2,NQP
         KC = NS2P2-K
         XH(1,K) = X(1,K)+X(1,KC)
         XH(2,K) = X(2,K)-X(2,KC)
  102 CONTINUE
      DO 103 K=2,NQP
         KC = NS2P2-K
         X(1,KC) = W(K-1)*XH(1,KC)+W(KC-1)*XH(2,KC)
         X(2,KC) = W(K-1)*XH(2,KC)-W(KC-1)*XH(1,KC)
  103 CONTINUE
      DO 104 K=2,NQP
         KC = NS2P2-K
         XH(1,KC) = X(1,KC)
         XH(2,KC) = X(2,KC)
  104 CONTINUE
      DO 105 K=2,NQP
         KC = NS2P2-K
         X(1,KC) = XH(1,K)-XH(1,KC)
         X(2,KC) = XH(2,K)-XH(2,KC)
  105 CONTINUE
      DO 106 K=2,NQP
         KC = NS2P2-K
         X(1,K) = XH(1,K)+XH(1,KC)
         X(2,K) = -XH(2,K)-XH(2,KC)
  106 CONTINUE
      IF (IPAR .NE. 0) GO TO 108
  107 X(1,NQP+1) = X(1,NQP+1)+X(1,NQP+1)
      X(2,NQP+1) = X(2,NQP+1)+X(2,NQP+1)
  108 XHOLD1 = X(2,1)+X(1,1)
      XHOLD2 = X(1,1)-X(2,1)
      X(2,1) = XHOLD2+XHOLD2
      X(1,1) = XHOLD1+XHOLD1
      RETURN
      END
      SUBROUTINE RFFTB (N,R,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE RFFTB(N,R,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE RFFTB COMPUTES THE REAL PERODIC SEQUENCE FROM ITS
C     FOURIER COEFFICIENTS (FOURIER SYNTHESIS). THE TRANSFORM IS DEFINED
C     BELOW AT OUTPUT PARAMETER R.
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE ARRAY R. N MUST BE EVEN AND THE METHOD
C             IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES.
C             N MAY CHANGE SO LONG AS DIFFERENT WORK ARRAYS ARE PROVIDED
C
C     R       A REAL ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE
C             TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 2.5*N+15
C             IN THE PROGRAM THAT CALLS RFFTB. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE RFFTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C             THE SAME WSAVE ARRAY CAN BE USED BY RFFTF AND RFFTB.
C
C
C     OUTPUT PARAMETERS
C
C     R       FOR I=1,...,N
C
C                  R(I)=X(1)+(-1)**(I+1)*X(2)
C
C                       PLUS THE SUM FROM K=2 TO K=N/2 OF
C
C                         2*R(2K-1)*COS((K-1)*(I-1)*2*PI/N)
C
C                        +2*R(2K)*SIN((K-1)*(I-1)*2*PI/N)
C
C      *****  NOTE
C                  THIS TRANSFORM IS UNNORMALIZED SINCE A CALL OF RFFTF
C                  FOLLOWED BY A CALL OF RFFTB WILL MULTIPLY THE INPUT
C                  SEQUENCE BY 2*N.
C
C     WSAVE   CONTAINS RESULTS WHICH MUST NOT BE DESTROYED BETWEEN
C             CALLS OF RFFTB OR RFFTF.
C
C
      DIMENSION       R(2,1)     ,WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HRFFTB  ,7HFFTPACK,10HVERSION  1)
      IF (N .GT. 2) GO TO 101
      R1 = R(1,1)+R(2,1)
      R(2,1) = R(1,1)-R(2,1)
      R(1,1) = R1
      RETURN
  101 IW1 = N/2+1
      CALL RFFTB1 (N,R,WSAVE(IW1),WSAVE)
      RETURN
      END
      SUBROUTINE RFFTB1 (N,X,XH,W)
      SAVE
      DIMENSION       X(2,1)     ,XH(2,1)    ,W(1)
      NS2 = N/2
      NS2P2 = NS2+2
      NQ = NS2/2
      IPAR = NS2-NQ-NQ
      NQM = NQ
      IF (IPAR .EQ. 0) NQM = NQM-1
      NQP = NQM+1
      XHOLD1 = X(1,1)-X(2,1)
      X(1,1) = X(2,1)+X(1,1)
      X(2,1) = XHOLD1
      IF (IPAR .NE. 0) GO TO 101
      X(1,NQP+1) = X(1,NQP+1)+X(1,NQP+1)
      X(2,NQP+1) = X(2,NQP+1)+X(2,NQP+1)
  101 IF (NQP .LT. 2) GO TO 108
      DO 102 K=2,NQP
         KC = NS2P2-K
         XH(1,K) = X(1,K)+X(1,KC)
         XH(2,K) = X(2,KC)-X(2,K)
  102 CONTINUE
      DO 103 K=2,NQP
         KC = NS2P2-K
         XH(1,KC) = X(1,K)-X(1,KC)
         XH(2,KC) = -X(2,K)-X(2,KC)
  103 CONTINUE
      DO 104 K=2,NQP
         KC = NS2P2-K
         X(1,KC) = XH(1,KC)
         X(2,KC) = XH(2,KC)
  104 CONTINUE
      DO 105 K=2,NQP
         KC = NS2P2-K
         XH(1,KC) = W(K-1)*X(1,KC)-W(KC-1)*X(2,KC)
         XH(2,KC) = W(K-1)*X(2,KC)+W(KC-1)*X(1,KC)
  105 CONTINUE
      DO 106 K=2,NQP
         KC = NS2P2-K
         X(1,K) = XH(1,K)-XH(2,KC)
         X(2,K) = XH(2,K)+XH(1,KC)
  106 CONTINUE
      DO 107 K=2,NQP
         KC = NS2P2-K
         X(1,KC) = XH(1,K)+XH(2,KC)
         X(2,KC) = XH(1,KC)-XH(2,K)
  107 CONTINUE
  108 CALL CFFTB (NS2,X,XH)
      RETURN
      END
      SUBROUTINE COSTI (N,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE COSTI(N,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE COSTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     SUBROUTINE COST. 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. N MUST BE
C             ODD AND GREATER THAN 1.
C
C     OUTPUT PARAMETER
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+12.
C             DIFFERENT WSAVE ARRAYS ARE REQUIRED FOR DIFFERENT VALUES
C             OF N. THE CONTENTS OF WSAVE MUST NOT BE CHANGED BETWEEN
C             CALLS OF COST.
C
      DIMENSION       WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCOSTI  ,7HFFTPACK,10HVERSION  1)
      IF (N .LE. 3) RETURN
      NM1 = N-1
      NS2 = NM1/2
      NS2M = NS2-1
      IW1 = NS2+1
      PI = 4.*ATAN(1.)
      DT = PI/FLOAT(NM1)
      DCS = COS(DT)
      DSS = SIN(DT)
      WSAVE(1) = DSS+DSS
      WCK = DCS+DCS
      IF (NS2M .LT. 2) GO TO 102
      DO 101 K=2,NS2M
         WSAVE(K) = DSS*WCK+DCS*WSAVE(K-1)
         WCK = DCS*WCK-DSS*WSAVE(K-1)
  101 CONTINUE
  102 CALL RFFTI (NM1,WSAVE(IW1))
      RETURN
      END
      SUBROUTINE COST (N,X,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE COST(N,X,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE COST COMPUTES THE DISCRETE FOURIER COSINE TRANSFORM
C     OF AN EVEN SEQUENCE X(I). THE TRANSFORM IS DEFINED BELOW AT OUTPUT
C     PARAMETER X.
C
C     COST IS THE UNNORMALIZED INVERSE OF ITSELF SINCE A CALL OF COST
C     FOLLOWED BY ANOTHER CALL OF COST WILL MULTIPLY THE INPUT SEQUENCE
C     X BY 8*(N+1). THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER X
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE COST MUST BE
C     INITIALIZED BY CALLING SUBROUTINE COSTI(N,WSAVE).
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE SEQUENCE X. N MUST BE ODD AND GREATER
C             THAN 1. THE METHOD IS MOST EFFICIENT WHEN N-1 IS A PRODUCT
C             OF SMALL PRIMES.
C
C     X       AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+15
C             IN THE PROGRAM THAT CALLS COST. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE COSTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C
C     OUTPUT PARAMETERS
C
C     X       FOR I=1,...,N
C
C                X(I)= 2*X(1)+2*(-1)**(I+1)*X(N)
C
C                  + THE SUM FROM K=2 TO K=N-1
C
C                    4*X(K)*COS((K-1)*(I-1)*PI/(N-1))
C
C                  A CALL OF COST FOLLOWED BY ANOTHER CALL OF
C                  COST WILL MULTIPLY THE SEQUENCE X BY 8(N+1)
C                  HENCE COST IS THE UNNORMALIZED INVERSE
C                  OF ITSELF.
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE
C             DESTROYED BETWEEN CALLS OF COST.
C
      DIMENSION       X(*)       ,WSAVE(*)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCOST   ,7HFFTPACK,10HVERSION  1)
      NM1 = N-1
      NP1 = N+1
      NS2 = NM1/2
      IW1 = NS2+1
      IF (N .GT. 3) GO TO 101
      X1PX3 = X(1)+X(3)
      TX13 = X1PX3+X1PX3
      FX2 = 4.*X(2)
      X(2) = 2.*(X(1)-X(3))
      X(1) = TX13+FX2
      X(3) = TX13-FX2
      RETURN
  101 C1 = X(1)-X(N)
      DO 102 K=2,NS2
         KC = NM1-K
         KS = NS2-K
         C1 = C1+WSAVE(KS+1)*(X(K)-X(KC+2))
  102 CONTINUE
      C1 = C1+C1
      X(1) = X(1)+X(N)
      DO 103 K=2,NS2
         KC = NP1-K
         T1 = X(K)+X(KC)
         T2 = WSAVE(K-1)*(X(K)-X(KC))
         X(K) = T1-T2
         X(KC) = T1+T2
  103 CONTINUE
      X(NS2+1) = X(NS2+1)+X(NS2+1)
      CALL RFFTF (NM1,X,WSAVE(IW1))
      CN = X(2)
      X(2) = C1
      DO 104 I=4,NM1,2
         X(I) = X(I)+X(I-2)
  104 CONTINUE
      X(N) = CN
      RETURN
      END
      SUBROUTINE SINTI (N,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE SINTI(N,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE SINTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     SUBROUTINE SINT. 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. N MUST BE
C             ODD.
C
C     OUTPUT PARAMETER
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+18.
C             DIFFERENT WSAVE ARRAYS ARE REQUIRED FOR DIFFERENT VALUES
C             OF N. THE CONTENTS OF WSAVE MUST NOT BE CHANGED BETWEEN
C             CALLS OF SINT.
C
      DIMENSION       WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HSINTI  ,7HFFTPACK,10HVERSION  1)
      IF (N .LE. 1) RETURN
      NP1 = N+1
      NS2 = NP1/2
      NS2M = NS2-1
      IW1 = NS2+1
      PI = 4.*ATAN(1.)
      DT = PI/FLOAT(NP1)
      DCS = COS(DT)
      DSS = SIN(DT)
      WSAVE(1) = DSS+DSS
      WCK = DCS+DCS
      IF (NS2M .LT. 2) GO TO 102
      DO 101 K=2,NS2M
         WSAVE(K) = DSS*WCK+DCS*WSAVE(K-1)
         WCK = DCS*WCK-DSS*WSAVE(K-1)
  101 CONTINUE
  102 CALL RFFTI (NP1,WSAVE(IW1))
      RETURN
      END
      SUBROUTINE SINT (N,X,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE SINT(N,X,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE SINT COMPUTES THE DISCRETE FOURIER SINE TRANSFORM
C     OF AN ODD SEQUENCE X(I). THE TRANSFORM IS DEFINED BELOW AT
C     OUTPUT PARAMETER X.
C
C     SINT IS THE UNNORMALIZED INVERSE OF ITSELF SINCE A CALL OF SINT
C     FOLLOWED BY ANOTHER CALL OF SINT WILL MULTIPLY THE INPUT SEQUENCE
C     X BY 8*(N+1).
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE SINT MUST BE
C     INITIALIZED BY CALLING SUBROUTINE SINTI(N,WSAVE).
C
C     INPUT PARAMETERS
C
C     N       N IS THE LENGTH OF X. N MUST BE ODD AND THE METHOD IS
C             MOST EFFICIENT WHEN N+1 IS A PRODUCT OF SMALL PRIMES
C
C     X       AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED
C
C                   ************IMPORTANT*************
C
C                   X MUST BE DIMENSIONED AT LEAST N+1
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3N+18
C             IN THE PROGRAM THAT CALLS SINT. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE SINTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C
C     OUTPUT PARAMETERS
C
C     X       FOR I=1,...,N
C
C                  X(I)= THE SUM FROM K=1 TO K=N
C
C                      4*X(K)*SIN(K*I*PI/(N+1))
C
C                  A CALL OF SINT FOLLOWED BY ANOTHER CALL OF
C                  SINT WILL MULTIPLY THE SEQUENCE X BY 8(N+1)
C                  HENCE SINT IS THE UNNORMALIZED INVERSE
C                  OF ITSELF.
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE
C             DESTROYED BETWEEN CALLS OF SINT.
C
      DIMENSION       X(1)       ,WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HSINT   ,7HFFTPACK,10HVERSION  1)
      IF (N .GT. 1) GO TO 101
      X(1) = 4.*X(1)
      RETURN
  101 NP1 = N+1
      NS2 = NP1/2
      NS2M = NS2-1
      IW1 = NS2+1
      X1 = X(1)
      X(1) = 0.
      DO 102 K=1,NS2M
         KC = NP1-K
         T1 = X1-X(KC)
         T2 = WSAVE(K)*(X1+X(KC))
         X1 = X(K+1)
         X(K+1) = T1+T2
         X(KC+1) = T2-T1
  102 CONTINUE
      X(NS2+1) = 4.*X1
      CALL RFFTF (NP1,X,WSAVE(IW1))
      X(1) = .5*X(1)
      DO 103 I=3,N,2
         X(I-1) = X(I+1)
         X(I) = X(I)+X(I-2)
  103 CONTINUE
      RETURN
      END
      SUBROUTINE COSQI (N,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE COSQI(N,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE COSQI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     BOTH COSQF AND COSQB. 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. N MUST BE
C             EVEN
C
C     OUTPUT PARAMETER
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3.5*N+15.
C             THE SAME WORK ARRAY CAN BE USED FOR BOTH COSQF AND COSQB
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 COSQF OR COSQB.
C
      DIMENSION       WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCOSQI  ,7HFFTPACK,10HVERSION  1)
      IW1 = N+1
      PIH = 2.*ATAN(1.)
      DT = PIH/FLOAT(N)
      DC = COS(DT)
      DS = SIN(DT)
      WSAVE(1) = DC
      WSK = DS
      DO 101 K=2,N
         WSAVE(K) = DC*WSAVE(K-1)-DS*WSK
         WSK = DS*WSAVE(K-1)+DC*WSK
  101 CONTINUE
      CALL RFFTI (N,WSAVE(IW1))
      RETURN
      END
      SUBROUTINE COSQF (N,X,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE COSQF(N,X,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE COSQF COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER
C     WAVE DATA. THAT IS , COSQF COMPUTES THE COEFFICIENTS IN A COSINE
C     SERIES REPRESENTATION WITH ONLY ODD WAVE NUMBERS. THE TRANSFORM
C     IS DEFINED BELOW AT OUTPUT PARAMETER X
C
C     COSQF IS THE UNNORMALIZED INVERSE OF COSQB SINCE A CALL OF COSQF
C     FOLLOWED BY A CALL OF COSQB WILL MULTIPLY THE INPUT SEQUENCE X
C     BY 8*N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE COSQF MUST BE
C     INITIALIZED BY CALLING SUBROUTINE COSQI(N,WSAVE).
C
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE ARRAY X. N MUST BE EVEN AND THE METHOD
C             IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES.
C
C     X       AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3.5N+15
C             IN THE PROGRAM THAT CALLS COSQF. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE COSQI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C
C     OUTPUT PARAMETERS
C
C     X       FOR I=1,...,N
C
C                  X(I)=2*X(1) PLUS THE SUM FROM K=2 TO K=N OF
C
C                     4.*X(K)*COS((2*I-1)*(K-1)*PI/(2*N))
C
C                  A CALL OF COSQF FOLLOWED BY A CALL OF
C                  COSQB WILL MULTIPLY THE SEQUENCE X BY 8*N
C                  THEREFORE COSQB IS THE UNNORMALIZED INVERSE
C                  OF COSQF.
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT
C             BE DESTROYED BETWEEN CALLS OF COSQF OR COSQB.
C
      DIMENSION       X(*)       ,WSAVE(*)
      DATA TSQ2/2.82842712474619/
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCOSQF  ,7HFFTPACK,10HVERSION  1)
      IF (N .GT. 2) GO TO 101
      TX = X(1)+X(1)
      TSQX = TSQ2*X(2)
      X(1) = TX+TSQX
      X(2) = TX-TSQX
      RETURN
  101 IW1 = N+1
      IW2 = N/2+IW1
      CALL COSQF1 (N,X,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE COSQF1 (N,X,W,W1,XH)
      SAVE
      DIMENSION       X(*)       ,XH(*)      ,W1(*)      ,W(*)
      NS2 = N/2
      NM1 = N-1
      NP2 = N+2
      DO 101 K=2,NS2
         KC = NP2-K
         XH(K) = X(K)+X(KC)
         XH(KC) = X(K)-X(KC)
  101 CONTINUE
      XH(NS2+1) = X(NS2+1)+X(NS2+1)
      DO 102 K=2,NS2
         KC = NP2-K
         X(K) = W(K-1)*XH(KC)+W(KC-1)*XH(K)
         X(KC) = W(K-1)*XH(K)-W(KC-1)*XH(KC)
  102 CONTINUE
      X(NS2+1) = W(NS2)*XH(NS2+1)
      CALL RFFTF (N,X,W1)
      XN = X(2)
      DO 103 I=3,NM1,2
         X(I-1) = X(I)+X(I+1)
         X(I) = X(I)-X(I+1)
  103 CONTINUE
      X(N) = XN
      RETURN
      END
      SUBROUTINE COSQB (N,X,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE COSQB(N,X,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE COSQB COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER
C     WAVE DATA. THAT IS , COSQB COMPUTES A SEQUENCE FROM ITS
C     REPRESENTATION IN TERMS OF A COSINE SERIES WITH ODD WAVE NUMBERS.
C     THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER X.
C
C     COSQB IS THE UNNORMALIZED INVERSE OF COSQF SINCE A CALL OF COSQB
C     FOLLOWED BY A CALL OF COSQF WILL MULTIPLY THE INPUT SEQUENCE X
C     BY 8*N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE COSQB MUST BE
C     INITIALIZED BY CALLING SUBROUTINE COSQI(N,WSAVE).
C
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE ARRAY X. N MUST BE EVEN AND THE METHOD
C             IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES.
C
C     X       AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3.5N+15
C             IN THE PROGRAM THAT CALLS COSQB. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE COSQI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C
C     OUTPUT PARAMETERS
C
C     X       FOR I=1,...,N
C
C                  X(I)= THE SUM FROM K=1 TO K=N OF
C
C                    4*X(K)*COS((2*K-1)*(I-1)*PI/(2*N))
C
C                  A CALL OF COSQB FOLLOWED BY A CALL OF
C                  COSQF WILL MULTIPLY THE SEQUENCE X BY 8*N
C                  THEREFORE COSQF IS THE UNNORMALIZED INVERSE
C                  OF COSQB.
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT
C             BE DESTROYED BETWEEN CALLS OF COSQB OR COSQF.
C
      DIMENSION       X(*)       ,WSAVE(*)
      DATA TSQ2/2.82842712474619/
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCOSQB  ,7HFFTPACK,10HVERSION  1)
      IF (N .GT. 2) GO TO 101
      X1 = 4.*(X(1)+X(2))
      X(2) = TSQ2*(X(1)-X(2))
      X(1) = X1
      RETURN
  101 IW1 = N+1
      IW2 = N/2+IW1
      CALL COSQB1 (N,X,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE COSQB1 (N,X,W,W1,XH)
      SAVE
      DIMENSION       X(*)       ,W(*)       ,W1(*)      ,XH(*)
         KC = NP2-K
         XH(K) = W(K-1)*X(KC)+W(KC-1)*X(K)
         XH(KC) = W(K-1)*X(K)-W(KC-1)*X(KC)
  102 CONTINUE
      X(NS2+1) = W(NS2)*(X(NS2+1)+X(NS2+1))
      DO 103 K=2,NS2
         KC = NP2-K
         X(K) = XH(K)+XH(KC)
         X(KC) = XH(K)-XH(KC)
  103 CONTINUE
      X(1) = X(1)+X(1)
      RETURN
      END
      SUBROUTINE SINQI (N,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE SINQI(N,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE SINQI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     BOTH SINQF AND SINQB. 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. N MUST BE
C             EVEN
C
C     OUTPUT PARAMETER
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3.5*N+15.
C             THE SAME WORK ARRAY CAN BE USED FOR BOTH SINQF AND SINQB
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 SINQF OR SINQB.
C
      DIMENSION       WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HSINQI  ,7HFFTPACK,10HVERSION  1)
      CALL COSQI (N,WSAVE)
      RETURN
      END
      SUBROUTINE SINQF (N,X,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE SINQF(N,X,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE SINQF COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER
C     WAVE DATA. THAT IS , SINQF COMPUTES THE COEFFICIENTS IN A SINE
C     SERIES REPRESENTATION WITH ONLY ODD WAVE NUMBERS. THE TRANSFORM
C     IS DEFINED BELOW AT OUTPUT PARAMETER X.
C
C     SINQB IS THE UNNORMALIZED INVERSE OF SINQF SINCE A CALL OF SINQF
C     FOLLOWED BY A CALL OF SINQB WILL MULTIPLY THE INPUT SEQUENCE X
C     BY 8*N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE SINQF MUST BE
C     INITIALIZED BY CALLING SUBROUTINE SINQI(N,WSAVE).
C
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE ARRAY X. N MUST BE EVEN AND THE METHOD
C             IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES.
C
C     X       AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3.5N+15
C             IN THE PROGRAM THAT CALLS SINQF. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE SINQI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C
C     OUTPUT PARAMETERS
C
C     X       FOR I=1,...,N
C
C                  X(I)=2*(-1)**(I+1)*X(N)
C
C                     + THE SUM FROM K=1 TO K=N-1 OF
C
C                     4*X(K)*SIN((2I-1)*K*PI/(2*N))
C
C                  A CALL OF SINQF FOLLOWED BY A CALL OF
C                  SINQB WILL MULTIPLY THE SEQUENCE X BY 8*N
C                  THEREFORE SINQB IS THE UNNORMALIZED INVERSE
C                  OF SINQF.
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT
C             BE DESTROYED BETWEEN CALLS OF SINQF OR SINQB.
C
      DIMENSION       X(*)       ,WSAVE(*)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HSINQF  ,7HFFTPACK,10HVERSION  1)
      NS2 = N/2
      DO 101 K=1,NS2
         KC = N-K
         XHOLD = X(K)
         X(K) = X(KC+1)
         X(KC+1) = XHOLD
  101 CONTINUE
      CALL COSQF (N,X,WSAVE)
      DO 102 K=2,N,2
         X(K) = -X(K)
  102 CONTINUE
      RETURN
      END
      SUBROUTINE SINQB (N,X,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE SINQB(N,X,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE SINQB COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER
C     WAVE DATA. THAT IS , SINQB COMPUTES A SEQUENCE FROM ITS
C     REPRESENTATION IN TERMS OF A SINE SERIES WITH ODD WAVE NUMBERS.
C     THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER X.
C
C     SINQF IS THE UNNORMALIZED INVERSE OF SINQB SINCE A CALL OF SINQB
C     FOLLOWED BY A CALL OF SINQF WILL MULTIPLY THE INPUT SEQUENCE X
C     BY 8*N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE SINQB MUST BE
C     INITIALIZED BY CALLING SUBROUTINE SINQI(N,WSAVE).
C
C
C     INPUT PARAMETERS
C
C     N       THE LENGTH OF THE ARRAY X. N MUST BE EVEN AND THE METHOD
C             IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES.
C
C     X       AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3.5N+15
C             IN THE PROGRAM THAT CALLS SINQB. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE SINQI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C
C     OUTPUT PARAMETERS
C
C     X       FOR I=1,...,N
C
C                  X(I)= THE SUM FROM K=1 TO K=N OF
C
C                    4*X(K)*SIN((2K-1)*I*PI/(2*N))
C
C                  A CALL OF SINQB FOLLOWED BY A CALL OF
C                  SINQF WILL MULTIPLY THE SEQUENCE X BY 8*N
C                  THEREFORE SINQF IS THE UNNORMALIZED INVERSE
C                  OF SINQB.
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT
C             BE DESTROYED BETWEEN CALLS OF SINQB OR SINQF.
C
      DIMENSION       X(*)       ,WSAVE(*)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HSINQB  ,7HFFTPACK,10HVERSION  1)
      NS2 = N/2
      DO 101 K=2,N,2
         X(K) = -X(K)
  101 CONTINUE
      CALL COSQB (N,X,WSAVE)
      DO 102 K=1,NS2
         KC = N-K
         XHOLD = X(K)
         X(K) = X(KC+1)
         X(KC+1) = XHOLD
  102 CONTINUE
      RETURN
      END
      SUBROUTINE CFFTI (N,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE CFFTI(N,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE CFFTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     BOTH CFFTF AND CFFTB. 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 CFFTF AND CFFTB
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 CFFTF OR CFFTB.
C
      DIMENSION       WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCFFTI  ,7HFFTPACK,10HVERSION  1)
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE CFFTI1 (N,WA,IFAC)
      SAVE
      DIMENSION       WA(*)      ,IFAC(*)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
      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.*ATAN(1.)
      ARG1 = TPI/FLOAT(N)
      DC = COS(ARG1)
      DS = SIN(ARG1)
      WA(1) = DC
      WA(2) = DS
      NT = N+N
      DO 108 I=4,NT,2
         WA(I-1) = DC*WA(I-3)-DS*WA(I-2)
         WA(I) = DS*WA(I-3)+DC*WA(I-2)
  108 CONTINUE
      RETURN
      END
      SUBROUTINE CFFTB (N,C,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE CFFTB(N,C,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE CFFTB COMPUTES THE BACKWARD COMPLEX DISCRETE FOURIER
C     TRANSFORM (THE FOURIER SYNTHESIS). EQUIVALENTLY , CFFTB COMPUTES
C     A COMPLEX PERIODIC SEQUENCE FROM ITS FOURIER COEFFICIENTS.
C     THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER C.
C
C     A CALL OF CFFTF FOLLOWED BY A CALL OF CFFTB WILL MULTIPLY THE
C     SEQUENCE BY N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE CFFTB MUST BE
C     INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE).
C
C     INPUT PARAMETERS
C
C
C     N      THE LENGTH OF THE COMPLEX SEQUENCE C. THE METHOD IS
C            MORE EFFICIENT WHEN N IS THE PRODUCT OF SMALL PRIMES.
C
C     C      A COMPLEX ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE
C
C     WSAVE   A REAL WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4N+15
C             IN THE PROGRAM THAT CALLS CFFTB. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C             THE SAME WSAVE ARRAY CAN BE USED BY CFFTF AND CFFTB.
C
C     OUTPUT PARAMETERS
C
C     C      FOR J=1,...,N
C
C                C(J)=THE SUM FROM K=1,...,N OF
C
C                      C(K)*EXP(I*J*K*2*PI/N)
C
C                            WHERE I=SQRT(-1)
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE
C             DESTROYED BETWEEN CALLS OF SUBROUTINE CFFTF OR CFFTB
      DIMENSION       C(1)       ,WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCFFTB  ,7HFFTPACK,10HVERSION  1)
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE CFFTB1 (N,C,CH,WA,IFAC)
      SAVE
      DIMENSION       CH(*)      ,C(*)       ,WA(*)      ,IFAC(*)
      NF = IFAC(2)
      L1 = 1
      DO 105 K1=1,NF
         IP = IFAC(K1+2)
         L2 = IP*L1
         IDO = N/L2
         IDOT = IDO+IDO
         IDL1 = IDOT*L1
         IF (IP .NE. 4) GO TO 101
         IX2 = L1+L1
         IX3 = IX2+L1
         CALL PASSB4 (IDOT,L1,IDL1,IX2,IX3,C,C,C,CH,CH,WA,WA,WA)
         GO TO 104
  101    IF (IP .NE. 2) GO TO 102
         CALL PASSB2 (IDOT,L1,IDL1,C,C,C,CH,CH,WA)
         GO TO 104
  102    IF (IP .NE. 3) GO TO 103
         IX2 = L1+L1
         CALL PASSB3 (IDOT,L1,IDL1,IX2,C,C,C,CH,CH,WA,WA)
         GO TO 104
  103    CALL PASSB (IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA)
  104    L1 = L2
  105 CONTINUE
      RETURN
      END
      SUBROUTINE PASSB2 (IDO,L1,IDL1,CC,C1,C2,CH,CH2,WA1)
      SAVE
      DIMENSION       CC(IDO,2,L1)           ,C1(IDO,L1,2)           ,
     1                C2(IDL1,2) ,CH(IDO,L1,2)           ,CH2(IDL1,2),
     2                WA1(L1,1)
      IDOT = IDO/2
      IF (IDO .LT. L1) GO TO 103
      DO 102 K=1,L1
         DO 101 I=1,IDO
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            CH(I,K,2) = CC(I,1,K)-CC(I,2,K)
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 I=1,IDO
         DO 104 K=1,L1
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            CH(I,K,2) = CC(I,1,K)-CC(I,2,K)
  104    CONTINUE
  105 CONTINUE
  106 DO 107 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  107 CONTINUE
      DO 108 K=1,L1
         C1(1,K,2) = CH(1,K,2)
         C1(2,K,2) = CH(2,K,2)
  108 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IF (IDOT .LT. L1) GO TO 111
      DO 110 K=1,L1
         DO 109 I=4,IDO,2
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)+WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)-
     1                    WA1(L1,I-2)*CH(I,K,2)
  109    CONTINUE
  110 CONTINUE
      RETURN
  111 DO 113 I=4,IDO,2
         DO 112 K=1,L1
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)+WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)-
     1                    WA1(L1,I-2)*CH(I,K,2)
  112    CONTINUE
  113 CONTINUE
      RETURN
      END
      SUBROUTINE PASSB3 (IDO,L1,IDL1,IX2,CC,C1,C2,CH,CH2,WA1,WA2)
      SAVE
      DIMENSION       CC(IDO,3,L1)           ,C1(IDO,L1,3)           ,
     1                C2(IDL1,3) ,CH(IDO,L1,3)           ,CH2(IDL1,3),
     2                WA1(L1,1)  ,WA2(IX2,1)
      DATA TAUR,TAUI /-.5,.866025403784439/
      IDOT = IDO/2
      IF (IDO .LT. L1) GO TO 103
      DO 102 K=1,L1
         DO 101 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
            CH(I,K,2) = CC(I,2,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)-CC(I,3,K)
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 I=1,IDO
         DO 104 K=1,L1
            CH(I,K,1) = CC(I,1,K)
            CH(I,K,2) = CC(I,2,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)-CC(I,3,K)
  104    CONTINUE
  105 CONTINUE
C
  106 DO 107 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)+CH2(IK,2)
         C2(IK,2) = CH2(IK,1)+TAUR*CH2(IK,2)
         C2(IK,3) = TAUI*CH2(IK,3)
  107 CONTINUE
      DO 108 IK=2,IDL1,2
         CH2(IK-1,2) = C2(IK-1,2)-C2(IK,3)
         CH2(IK-1,3) = C2(IK-1,2)+C2(IK,3)
  108 CONTINUE
      DO 109 IK=2,IDL1,2
         CH2(IK,2) = C2(IK,2)+C2(IK-1,3)
         CH2(IK,3) = C2(IK,2)-C2(IK-1,3)
  109 CONTINUE
      DO 111 J=2,3
         DO 110 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  110    CONTINUE
  111 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IF (IDOT-1 .LT. L1) GO TO 114
      DO 113 K=1,L1
         DO 112 I=4,IDO,2
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)+WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)-
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)+
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)-
     1                    WA2(IX2,I-2)*CH(I,K,3)
  112    CONTINUE
  113 CONTINUE
      RETURN
  114 DO 116 I=4,IDO,2
         DO 115 K=1,L1
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)+WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)-
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)+
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)-
     1                    WA2(IX2,I-2)*CH(I,K,3)
  115    CONTINUE
  116 CONTINUE
      RETURN
      END
      SUBROUTINE PASSB4 (IDO,L1,IDL1,IX2,IX3,CC,C1,C2,CH,CH2,WA1,WA2,
     1                   WA3)
      SAVE
      DIMENSION       CC(IDO,4,L1)           ,C1(IDO,L1,4)           ,
     1                C2(IDL1,4) ,CH(IDO,L1,4)           ,CH2(IDL1,4),
     2                WA1(L1,1)  ,WA2(IX2,1) ,WA3(IX3,1)
      IDOT = IDO/2
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 K=1,L1
         DO 101 I=2,IDO,2
            CH(I-1,K,4) = CC(I,4,K)-CC(I,2,K)
  101    CONTINUE
         DO 102 I=2,IDO,2
            CH(I,K,4) = CC(I-1,2,K)-CC(I-1,4,K)
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,2) = CC(I,1,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)+CC(I,4,K)
            CH(I,K,1) = CC(I,1,K)-CC(I,3,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 111
  106 DO 108 I=2,IDO,2
         DO 107 K=1,L1
            CH(I-1,K,4) = CC(I,4,K)-CC(I,2,K)
            CH(I,K,4) = CC(I-1,2,K)-CC(I-1,4,K)
  107    CONTINUE
  108 CONTINUE
      DO 110 I=1,IDO
         DO 109 K=1,L1
            CH(I,K,2) = CC(I,1,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)+CC(I,4,K)
            CH(I,K,1) = CC(I,1,K)-CC(I,3,K)
  109    CONTINUE
  110 CONTINUE
  111 DO 112 IK=1,IDL1
         C2(IK,1) = CH2(IK,2)+CH2(IK,3)
  112 CONTINUE
      DO 113 IK=1,IDL1
         CH2(IK,3) = CH2(IK,2)-CH2(IK,3)
  113 CONTINUE
      DO 114 IK=1,IDL1
         CH2(IK,2) = CH2(IK,1)+CH2(IK,4)
  114 CONTINUE
      DO 115 IK=1,IDL1
         CH2(IK,4) = CH2(IK,1)-CH2(IK,4)
  115 CONTINUE
      DO 117 J=2,4
         DO 116 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  116    CONTINUE
  117 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IF (IDOT .LT. L1) GO TO 120
      DO 119 K=1,L1
         DO 118 I=4,IDO,2
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)+WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)-
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)+
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)-
     1                    WA2(IX2,I-2)*CH(I,K,3)
            C1(I,K,4) = WA3(IX3-1,I-2)*CH(I,K,4)+
     1                  WA3(IX3,I-2)*CH(I-1,K,4)
            C1(I-1,K,4) = WA3(IX3-1,I-2)*CH(I-1,K,4)-
     1                    WA3(IX3,I-2)*CH(I,K,4)
  118    CONTINUE
  119 CONTINUE
      RETURN
  120 DO 122 I=4,IDO,2
         DO 121 K=1,L1
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)+WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)-
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)+
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)-
     1                    WA2(IX2,I-2)*CH(I,K,3)
            C1(I,K,4) = WA3(IX3-1,I-2)*CH(I,K,4)+
     1                  WA3(IX3,I-2)*CH(I-1,K,4)
            C1(I-1,K,4) = WA3(IX3-1,I-2)*CH(I-1,K,4)-
     1                    WA3(IX3,I-2)*CH(I,K,4)
  121    CONTINUE
  122 CONTINUE
      RETURN
      END
      SUBROUTINE PASSB (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
      SAVE
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
     1                C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP),
     2                CH2(IDL1,IP)
      IDOT = IDO/2
      NT = IP*IDL1
      IPP2 = IP+2
      IPPH = (IP+1)/2
      L1T = L1+L1
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 J=2,IPPH
         JC = IPP2-J
         DO 102 K=1,L1
            DO 101 I=1,IDO
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       CONTINUE
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 112
  106 DO 109 J=2,IPPH
         JC = IPP2-J
         DO 108 I=1,IDO
            DO 107 K=1,L1
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       CONTINUE
  108    CONTINUE
  109 CONTINUE
      DO 111 I=1,IDO
         DO 110 K=1,L1
            CH(I,K,1) = CC(I,1,K)
  110    CONTINUE
  111 CONTINUE
  112 DO 113 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  113 CONTINUE
      IDJ = 0
      DO 115 J=2,IPPH
         JC = IPP2-J
         IDJ = IDJ+IDL1
         DO 114 IK=1,IDL1
            C2(IK,J) = CH2(IK,1)+WA(IDJ-1)*CH2(IK,2)
            C2(IK,JC) = WA(IDJ)*CH2(IK,IP)
  114    CONTINUE
  115 CONTINUE
      DO 117 J=2,IPPH
         DO 116 IK=1,IDL1
            C2(IK,1) = C2(IK,1)+CH2(IK,J)
  116    CONTINUE
  117 CONTINUE
C
      IDL = 0
      DO 120 L=2,IPPH
         LC = IPP2-L
         IDL = IDL+IDL1
         IDLJ = IDL
         DO 119 J=3,IPPH
            JC = IPP2-J
            IDLJ = IDLJ+IDL
            IF (IDLJ .GT. NT) IDLJ = IDLJ-NT
            WAR = WA(IDLJ-1)
            WAI = WA(IDLJ)
            DO 118 IK=1,IDL1
               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
               C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC)
  118       CONTINUE
  119    CONTINUE
  120 CONTINUE
C
      DO 122 J=2,IPPH
         JC = IPP2-J
         DO 121 IK=2,IDL1,2
            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
  121    CONTINUE
  122 CONTINUE
      DO 124 J=2,IPPH
         JC = IPP2-J
         DO 123 IK=2,IDL1,2
            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  123    CONTINUE
  124 CONTINUE
C
      DO 126 J=2,IP
         DO 125 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  125    CONTINUE
  126 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IDJ = 0
      IF (IDOT .GT. L1) GO TO 130
      DO 129 J=2,IP
         IDJ = IDJ+L1T
         IDIJ = 0
         DO 128 I=4,IDO,2
            IDIJ = IDIJ+IDJ
            DO 127 K=1,L1
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  127       CONTINUE
  128    CONTINUE
  129 CONTINUE
      RETURN
  130 DO 134 J=2,IP
         IDJ = IDJ+L1T
         DO 133 K=1,L1
            IDIJ = 0
            DO 131 I=4,IDO,2
               IDIJ = IDIJ+IDJ
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
  131       CONTINUE
            IDIJ = 0
            DO 132 I=4,IDO,2
               IDIJ = IDIJ+IDJ
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  132       CONTINUE
  133    CONTINUE
  134 CONTINUE
      RETURN
      END
      SUBROUTINE CFFTF (N,C,WSAVE)
      SAVE
C     ******************************************************************
C
C     SUBROUTINE CFFTF(N,C,WSAVE)
C
C     ******************************************************************
C
C     SUBROUTINE CFFTF COMPUTES THE FORWARD COMPLEX DISCRETE FOURIER
C     TRANSFORM (THE FOURIER ANALYSIS). EQUIVALENTLY , CFFTF COMPUTES
C     THE FOURIER COEFFICIENTS OF A COMPLEX PERIODIC SEQUENCE.
C     THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER C.
C
C     THE TRANSFORM IS NOT NORMALIZED. TO OBTAIN A NORMALIZED TRANSFORM
C     THE OUTPUT MUST BE DIVIDED BY N. OTHERWISE A CALL OF CFFTF
C     FOLLOWED BY A CALL OF CFFTB WILL MULTIPLY THE SEQUENCE BY N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE CFFTF MUST BE
C     INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE).
C
C     INPUT PARAMETERS
C
C
C     N      THE LENGTH OF THE COMPLEX SEQUENCE C. THE METHOD IS
C            MORE EFFICIENT WHEN N IS THE PRODUCT OF SMALL PRIMES. N
C
C     C      A COMPLEX ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE
C
C     WSAVE   A REAL WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4N+15
C             IN THE PROGRAM THAT CALLS CFFTF. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C             THE SAME WSAVE ARRAY CAN BE USED BY CFFTF AND CFFTB.
C
C     OUTPUT PARAMETERS
C
C     C      FOR J=1,...,N
C
C                C(J)=THE SUM FROM K=1,...,N OF
C
C                      C(K)*EXP(-I*J*K*2*PI/N)
C
C                            WHERE I=SQRT(-1)
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE
C             DESTROYED BETWEEN CALLS OF SUBROUTINE CFFTF OR CFFTB
C
      DIMENSION       C(1)       ,WSAVE(1)
C
C THE FOLLOWING CALL IS FOR GATHERING LIBRARY STATISTICS AT NCAR
C
      CALL Q8QST4(4HXLIB,7HCFFTF  ,7HFFTPACK,10HVERSION  1)
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE CFFTF1 (N,C,CH,WA,IFAC)
      SAVE
      DIMENSION       CH(*)      ,C(*)       ,WA(*)      ,IFAC(*)
      NF = IFAC(2)
      L1 = 1
      DO 105 K1=1,NF
         IP = IFAC(K1+2)
         L2 = IP*L1
         IDO = N/L2
         IDOT = IDO+IDO
         IDL1 = IDOT*L1
         IF (IP .NE. 4) GO TO 101
         IX2 = L1+L1
         IX3 = IX2+L1
         CALL PASSF4 (IDOT,L1,IDL1,IX2,IX3,C,C,C,CH,CH,WA,WA,WA)
         GO TO 104
  101    IF (IP .NE. 2) GO TO 102
         CALL PASSF2 (IDOT,L1,IDL1,C,C,C,CH,CH,WA)
         GO TO 104
  102    IF (IP .NE. 3) GO TO 103
         IX2 = L1+L1
         CALL PASSF3 (IDOT,L1,IDL1,IX2,C,C,C,CH,CH,WA,WA)
         GO TO 104
  103    CALL PASSF (IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA)
  104    L1 = L2
  105 CONTINUE
      RETURN
      END
      SUBROUTINE PASSF2 (IDO,L1,IDL1,CC,C1,C2,CH,CH2,WA1)
      SAVE
      DIMENSION       CC(IDO,2,L1)           ,C1(IDO,L1,2)           ,
     1                C2(IDL1,2) ,CH(IDO,L1,2)           ,CH2(IDL1,2),
     2                WA1(L1,1)
      IDOT = IDO/2
      IF (IDO .LT. L1) GO TO 103
      DO 102 K=1,L1
         DO 101 I=1,IDO
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            CH(I,K,2) = CC(I,1,K)-CC(I,2,K)
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 I=1,IDO
         DO 104 K=1,L1
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            CH(I,K,2) = CC(I,1,K)-CC(I,2,K)
  104    CONTINUE
  105 CONTINUE
  106 DO 107 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  107 CONTINUE
      DO 108 K=1,L1
         C1(1,K,2) = CH(1,K,2)
         C1(2,K,2) = CH(2,K,2)
  108 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IF (IDOT .LT. L1) GO TO 111
      DO 110 K=1,L1
         DO 109 I=4,IDO,2
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)-WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)+
     1                    WA1(L1,I-2)*CH(I,K,2)
  109    CONTINUE
  110 CONTINUE
      RETURN
  111 DO 113 I=4,IDO,2
         DO 112 K=1,L1
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)-WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)+
     1                    WA1(L1,I-2)*CH(I,K,2)
  112    CONTINUE
  113 CONTINUE
      RETURN
      END
      SUBROUTINE PASSF3 (IDO,L1,IDL1,IX2,CC,C1,C2,CH,CH2,WA1,WA2)
      SAVE
      DIMENSION       CC(IDO,3,L1)           ,C1(IDO,L1,3)           ,
     1                C2(IDL1,3) ,CH(IDO,L1,3)           ,CH2(IDL1,3),
     2                WA1(L1,1)  ,WA2(IX2,1)
      DATA TAUR,TAUI /-.5,-.866025403784439/
      IDOT = IDO/2
      IF (IDO .LT. L1) GO TO 103
      DO 102 K=1,L1
         DO 101 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
            CH(I,K,2) = CC(I,2,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)-CC(I,3,K)
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 I=1,IDO
         DO 104 K=1,L1
            CH(I,K,1) = CC(I,1,K)
            CH(I,K,2) = CC(I,2,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)-CC(I,3,K)
  104    CONTINUE
  105 CONTINUE
C
  106 DO 107 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)+CH2(IK,2)
         C2(IK,2) = CH2(IK,1)+TAUR*CH2(IK,2)
         C2(IK,3) = TAUI*CH2(IK,3)
  107 CONTINUE
      DO 108 IK=2,IDL1,2
         CH2(IK-1,2) = C2(IK-1,2)-C2(IK,3)
         CH2(IK-1,3) = C2(IK-1,2)+C2(IK,3)
  108 CONTINUE
      DO 109 IK=2,IDL1,2
         CH2(IK,2) = C2(IK,2)+C2(IK-1,3)
         CH2(IK,3) = C2(IK,2)-C2(IK-1,3)
  109 CONTINUE
      DO 111 J=2,3
         DO 110 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  110    CONTINUE
  111 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IF (IDOT-1 .LT. L1) GO TO 114
      DO 113 K=1,L1
         DO 112 I=4,IDO,2
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)-WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)+
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)-
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)+
     1                    WA2(IX2,I-2)*CH(I,K,3)
  112    CONTINUE
  113 CONTINUE
      RETURN
  114 DO 116 I=4,IDO,2
         DO 115 K=1,L1
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)-WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)+
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)-
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)+
     1                    WA2(IX2,I-2)*CH(I,K,3)
  115    CONTINUE
  116 CONTINUE
      RETURN
      END
      SUBROUTINE PASSF4 (IDO,L1,IDL1,IX2,IX3,CC,C1,C2,CH,CH2,WA1,WA2,
     1                   WA3)
      SAVE
      DIMENSION       CC(IDO,4,L1)           ,C1(IDO,L1,4)           ,
     1                C2(IDL1,4) ,CH(IDO,L1,4)           ,CH2(IDL1,4),
     2                WA1(L1,1)  ,WA2(IX2,1) ,WA3(IX3,1)
      IDOT = IDO/2
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 K=1,L1
         DO 101 I=2,IDO,2
            CH(I-1,K,4) = CC(I,2,K)-CC(I,4,K)
  101    CONTINUE
         DO 102 I=2,IDO,2
            CH(I,K,4) = CC(I-1,4,K)-CC(I-1,2,K)
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,2) = CC(I,1,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)+CC(I,4,K)
            CH(I,K,1) = CC(I,1,K)-CC(I,3,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 111
  106 DO 108 I=2,IDO,2
         DO 107 K=1,L1
            CH(I,K,4) = CC(I-1,4,K)-CC(I-1,2,K)
            CH(I-1,K,4) = CC(I,2,K)-CC(I,4,K)
  107    CONTINUE
  108 CONTINUE
      DO 110 I=1,IDO
         DO 109 K=1,L1
            CH(I,K,2) = CC(I,1,K)+CC(I,3,K)
            CH(I,K,3) = CC(I,2,K)+CC(I,4,K)
            CH(I,K,1) = CC(I,1,K)-CC(I,3,K)
  109    CONTINUE
  110 CONTINUE
  111 DO 112 IK=1,IDL1
         C2(IK,1) = CH2(IK,2)+CH2(IK,3)
  112 CONTINUE
      DO 113 IK=1,IDL1
         CH2(IK,3) = CH2(IK,2)-CH2(IK,3)
  113 CONTINUE
      DO 114 IK=1,IDL1
         CH2(IK,2) = CH2(IK,1)+CH2(IK,4)
  114 CONTINUE
      DO 115 IK=1,IDL1
         CH2(IK,4) = CH2(IK,1)-CH2(IK,4)
  115 CONTINUE
      DO 117 J=2,4
         DO 116 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  116    CONTINUE
  117 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IF (IDOT .LT. L1) GO TO 120
      DO 119 K=1,L1
         DO 118 I=4,IDO,2
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)-WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)+
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)-
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)+
     1                    WA2(IX2,I-2)*CH(I,K,3)
            C1(I,K,4) = WA3(IX3-1,I-2)*CH(I,K,4)-
     1                  WA3(IX3,I-2)*CH(I-1,K,4)
            C1(I-1,K,4) = WA3(IX3-1,I-2)*CH(I-1,K,4)+
     1                    WA3(IX3,I-2)*CH(I,K,4)
  118    CONTINUE
  119 CONTINUE
      RETURN
  120 DO 122 I=4,IDO,2
         DO 121 K=1,L1
            C1(I,K,2) = WA1(L1-1,I-2)*CH(I,K,2)-WA1(L1,I-2)*CH(I-1,K,2)
            C1(I-1,K,2) = WA1(L1-1,I-2)*CH(I-1,K,2)+
     1                    WA1(L1,I-2)*CH(I,K,2)
            C1(I,K,3) = WA2(IX2-1,I-2)*CH(I,K,3)-
     1                  WA2(IX2,I-2)*CH(I-1,K,3)
            C1(I-1,K,3) = WA2(IX2-1,I-2)*CH(I-1,K,3)+
     1                    WA2(IX2,I-2)*CH(I,K,3)
            C1(I,K,4) = WA3(IX3-1,I-2)*CH(I,K,4)-
     1                  WA3(IX3,I-2)*CH(I-1,K,4)
            C1(I-1,K,4) = WA3(IX3-1,I-2)*CH(I-1,K,4)+
     1                    WA3(IX3,I-2)*CH(I,K,4)
  121    CONTINUE
  122 CONTINUE
      RETURN
      END
      SUBROUTINE PASSF (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
      SAVE
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
     1                C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP),
     2                CH2(IDL1,IP)
      IDOT = IDO/2
      NT = IP*IDL1
      IPP2 = IP+2
      IPPH = (IP+1)/2
      L1T = L1+L1
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 J=2,IPPH
         JC = IPP2-J
         DO 102 K=1,L1
            DO 101 I=1,IDO
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       CONTINUE
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 112
  106 DO 109 J=2,IPPH
         JC = IPP2-J
         DO 108 I=1,IDO
            DO 107 K=1,L1
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       CONTINUE
  108    CONTINUE
  109 CONTINUE
      DO 111 I=1,IDO
         DO 110 K=1,L1
            CH(I,K,1) = CC(I,1,K)
  110    CONTINUE
  111 CONTINUE
  112 DO 113 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  113 CONTINUE
      IDJ = 0
      DO 115 J=2,IPPH
         JC = IPP2-J
         IDJ = IDJ+IDL1
         DO 114 IK=1,IDL1
            C2(IK,J) = CH2(IK,1)+WA(IDJ-1)*CH2(IK,2)
            C2(IK,JC) = -WA(IDJ)*CH2(IK,IP)
  114    CONTINUE
  115 CONTINUE
      DO 117 J=2,IPPH
         DO 116 IK=1,IDL1
            C2(IK,1) = C2(IK,1)+CH2(IK,J)
  116    CONTINUE
  117 CONTINUE
C
      IDL = 0
      DO 120 L=2,IPPH
         LC = IPP2-L
         IDL = IDL+IDL1
         IDLJ = IDL
         DO 119 J=3,IPPH
            JC = IPP2-J
            IDLJ = IDLJ+IDL
            IF (IDLJ .GT. NT) IDLJ = IDLJ-NT
            WAR = WA(IDLJ-1)
            WAI = WA(IDLJ)
            DO 118 IK=1,IDL1
               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
               C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC)
  118       CONTINUE
  119    CONTINUE
  120 CONTINUE
C
      DO 122 J=2,IPPH
         JC = IPP2-J
         DO 121 IK=2,IDL1,2
            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
  121    CONTINUE
  122 CONTINUE
      DO 124 J=2,IPPH
         JC = IPP2-J
         DO 123 IK=2,IDL1,2
            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  123    CONTINUE
  124 CONTINUE
C
      DO 126 J=2,IP
         DO 125 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  125    CONTINUE
  126 CONTINUE
      IF (IDO .EQ. 2) RETURN
      IDJ = 0
      IF (IDOT .GT. L1) GO TO 130
      DO 129 J=2,IP
         IDJ = IDJ+L1T
         IDIJ = 0
         DO 128 I=4,IDO,2
            IDIJ = IDIJ+IDJ
            DO 127 K=1,L1
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  127       CONTINUE
  128    CONTINUE
  129 CONTINUE
      RETURN
  130 DO 134 J=2,IP
         IDJ = IDJ+L1T
         DO 133 K=1,L1
            IDIJ = 0
            DO 131 I=4,IDO,2
               IDIJ = IDIJ+IDJ
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
  131       CONTINUE
            IDIJ = 0
            DO 132 I=4,IDO,2
               IDIJ = IDIJ+IDJ
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  132       CONTINUE
  133    CONTINUE
  134 CONTINUE
      RETURN
      END
