      SUBROUTINE CFFT( DATA, N, IFORW )
      DIMENSION DATA( N )
      ISIGN = -IFORW
      IP0 = 2
      IP3 = IP0 * N
      I3REV = 1

      DO I3 = 1, IP3, IP0
         IF ( I3 .LT. I3REV ) THEN
            TEMPR = DATA(I3)
            TEMPI = DATA(I3 + 1)
            DATA(I3) = DATA(I3REV)
            DATA(I3 + 1) = DATA(I3REV + 1)
            DATA(I3REV) = TEMPR
            DATA(I3REV + 1) = TEMPI
         ENDIF
         IP1 = IP3 / 2
30       IF( I3REV .LE. IP1 ) GOTO 50
            I3REV = I3REV - IP1
            IP1 = IP1 / 2
            IF(IP1 .GE. IP0) GOTO 30
50       I3REV = I3REV + IP1
      END DO

      IP1 = IP0

60    IF(IP1 .GE. IP3) GOTO 100
      IP2 = IP1 * 2
      THETA = 6.283185307 / FLOAT(ISIGN * IP2 / IP0)  
      SINTH = SIN(THETA / 2.)
      WSTPR = -2. * SINTH * SINTH
      WSTPI = SIN(THETA)
      WR = 1.
      WI = 0.
      DO I1 = 1, IP1, IP0
         DO I3 = I1, IP3, IP2
            I2A = I3
            I2B = I2A + IP1
            TEMPR = WR * DATA(I2B) - WI * DATA(I2B + 1)
            TEMPI = WR * DATA(I2B + 1) + WI * DATA(I2B)
            DATA(I2B) = DATA(I2A) - TEMPR
            DATA(I2B + 1) = DATA(I2A + 1) - TEMPI
            DATA(I2A) = DATA(I2A) + TEMPR
            DATA(I2A + 1) = DATA(I2A + 1) + TEMPI
         END DO

         TEMPR = WR
         WR = WR * WSTPR - WI * WSTPI + WR
         WI = WI * WSTPR + TEMPR * WSTPI + WI
      END DO

      IP1 = IP2
      GO TO 60

100   RETURN
      END
C**********************************************************************C
      SUBROUTINE CFFTSC( C, N )

C     SCALE COMPLEX FORWARD FFT RESULTS

C     M. PORTER 4 / 1 / 88

      COMPLEX C( * )

      DO I  = 1, N
         C(I)  = C(I) / N
      ENDDO

      RETURN
      END
