      SUBROUTINE SCFT (A, MODE, M, ND, MS, S)
c>> 1997-03-31 SCFT Krogh  Increased KEDIM, more sine table checks.
C>> 1994-10-20 SCFT Krogh  Change CERROR to simplify conversion to C.
C>> 1994-10-20 SCFT Krogh  Changes to use M77CON
C>> 1994-04-19 SCFT CLL Edited to make DP & SP files similar.
c>> 1989-06-16 FTK Fix error message on MODE.
c>> 1989-06-05 WVS Change length of MODE from *(ND) to *(*).
c>> 1989-05-07 FTK & CLL
c>> 1989-04-21 FTK & CLL
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c
c     This subroutine computes complex Fourier transforms in up to 6
c     dimensions using the Cooley-Tukey fast Fourier transform.
c
c     Variables in the calling sequence have the following types.
      REAL             A(*), S(*)
      INTEGER  ND, M(ND), MS
      CHARACTER MODE*(*)
c
c     Programmed by Fred T. Krogh at the Jet Propulsion Laboratory,
c     Pasadena, Calif.   August 1, 1969.
c
c     Values for A, MODE, M, ND, and MS must be specified before
c     calling the subroutine.
c
c     In describing the usage the following notation is used
c     N(K) = 2**M(K)
c     WK = EXP(2*PI*I/N(K)), where I = sqrt(-1) and PI = 3.14159...
c     MA = M(1) + M(2) + ... + M(ND)
c     NA = 2**MA
c     MX = MAX(M(1), M(2), ..., M(ND))
c     NX = 2**MX
c
c The usage is as follows
c
c A() on input is an array of function values, X, if one is doing
c   Fourier analysis,
c   and is an array of Fourier coefficients, C, if one is
c   doing Fourier synthesis.  On output, A contains C if doing analysis,
c   and contains X if doing synthesis.  In our description here, we
c   assume that A has been declared by the user as a complex array with
c   dimension A(N(1),N(2),...,N(ND)). (Of course this is not necessary
c   as long as one declares enough storage for A, and keeps in mind that
c   the array is treated in this subroutine as if it had been specified
c   in this way.)  It is assumed that the imaginary part of a complex
c   number is stored in the cell immediately following its real part.
c
c MODE  A character variable of length at least ND.
c   The k-th character must be 'A' or 'a' to select Analysis in the
c   k-th dimension, or 'S' or 's' to select Synthesis in the
c   k-th dimension.
c   The subroutine sets  A(j1+1, j2+1, ..., jND+1) to the sum
c  over 0 .le. k1 .le. N(1)-1, 0 .le. k2 .le. N(2)-1, ...,
c  0 .le. kND .le. N(ND)-1 of A(k1+1, k2+1, ..., kND+1) * T(1,j1,k1) *
c  T(2,j2,k2) * ... * T(ND,jND,kND), 0 .le. j1 .le. N(1)-1, ...,
c  0 .le. jnd .le. N(ND)-1
c  with T(L,j,k) =
c        WL ** (j*k)               if MODE(L:L) = 'S' or 's'
c        (1/N(L)) * WL ** (-j*k)   if MODE(L:L) = 'A' or 'a'
c
c M() is a vector used to indicate N(k) = 2**M(k), the number of
c  complex points in the k-th variable.  M must be such that MX < 21.
c  If M(L)=0, no change is made with respect to the L-th dimension.
c
c ND is the dimension of (i.e. the number of subscripts in) the
c    array A.  ND must satisfy 1 .le. ND .le. 6
c
c MS gives the state of the sine table.  If MS > 0, there are NT =
c    2 ** (MS-2) good entries in the sine table.  On the initial call,
c    MS must be set to 0, and when the return is made, it will be set
c    to MX, which is the value of MS required for computing a
c    transform of size NX.  If MS = -1, the sine table will be computed
c    as for the case MS = 0, and then a return to the user will be made
c    with MS set as before, but no transform will be computed.  This
c    option is useful if the user would like access to the sine table
c    before computing the FFT.
c
c S is a vector, S(j) = sin(pi*j/(2*NT)), j = 1, 2, ..., NT-1, where
c  NT is defined in the description of MS above.  S is computed by the
c  subroutine if MX .gt. MS.  (If S is altered, set MS=0 so that S
c  is recomputed.)
c     ------------------------------------------------------------------
c                Notes on COMMON, PARAMETER's, and local variables
C
c     NDMAX = the maximum value for ND, and MAXM = the maximum
c     permitted for M(1), ..., M(ND)
c     The dimension of KE must be at least as large as MAXM.
c     The named common CSFFTC is used for communication between this
c     subroutine and the subroutine SFFT which computes a one
c     dimensional complex Fourier transform and computes the sine table.
c     The use of the variables in CSFFTC is contained in the listing
c     of SFFT.
c
c     NF(1) = 2, NF(K+1) = 2**(1+M(1)+...+M(K)), K = 1,..., ND
c     ------------------------------------------------------------------
c--S replaces "?": ?CFT, ?FFT, C?FFTC
c     Both versions use IERM1
c     ------------------------------------------------------------------
      INTEGER MAXM, NDMAX
      PARAMETER (NDMAX = 6)

      INTEGER I, II, IIC, IR, IRC
      INTEGER J, JL
      INTEGER K, KEDIM, KSM
      INTEGER L
      INTEGER MA, MMAX, MSI
      INTEGER NDD, NF(NDMAX+1), NTOT2

      CHARACTER*13  CERROR

      REAL             FN, ONE, SPI4

      PARAMETER (ONE = 1.E0)
      PARAMETER (SPI4 = .70710 67811 86547 52440 08443 62104 8490E0)

c Common variables
      PARAMETER (KEDIM=30)
      LOGICAL NEEDST
      INTEGER MT, NT, MM, KS, ILAST, KE(KEDIM), KEE(KEDIM+1)
c Note that KEE(1) is equivalent to ILAST.
      EQUIVALENCE (KE(1), KEE(2))
      COMMON /CSFFTC/ NEEDST, MT, NT, MM, KS, ILAST, KE
      SAVE /CSFFTC/
      PARAMETER (MAXM = KEDIM)
      DATA CERROR / 'MODE(K:K) =  ' /
c     ------------------------------------------------------------------
c
      NDD = ND
      IF ((NDD .LE. 0) .OR. (NDD .GT. NDMAX)) THEN
c                               Fatal error, default is to stop in IERM1
         CALL IERM1 ('SCFT', 1, 2, 'BAD ND', 'ND', ND, '.')
         MS = -2
         RETURN
      END IF
      NF(1) = 2
      MA = 0
      MMAX = 0
      FN = ONE
      DO 10 K = 1, NDD
         MM = M(K)
         MMAX = MAX(MM, MMAX)
         MA = MA + MM
         I = 2**MM
         IF ((MM .LT. 0) .OR. (MM .GT. MAXM)) THEN
c                               Fatal error, default is to stop in IERM1
         CALL IERM1('SCFT',2,2,'Need 0 .le. M(K) .le. 30','M(K)',MM,'.')
         MS = -2
         RETURN
         END IF
         NF(K+1) = NF(K) * I
         if (MODE(K:K) .eq. 'A' .or. MODE(K:K) .eq. 'a') then
            FN = FN / I
         else if (MODE(K:K) .ne. 'S' .and. MODE(K:K) .ne. 's') then
            CERROR(13:13) = MODE(K:K)
            CALL IERM1 ('SCFT',3,2,CERROR, 'for K =',K,'.')
            MS = -2
            RETURN
         end if
   10 CONTINUE
      MSI = MS
      NEEDST = MMAX .GT. MSI
      if (.NOT. NEEDST) then
c  Check internal parameters to catch certain user errors.
         if (MT .lt. KEDIM) then
            if (MMAX .le. MT + 2) then
c              Skip sine table computation if all appears O.K.
               if (MT .le. 0) go to 20
               if (abs(S(NT/2) - SPI4) .le. 1.E-7) go to 20
            end if
         end if
         NEEDST = .true.
         call ERMSG('DRFT1', 3, 1,
     1      'Invalid sine table (re)computed', '.')
      end if
      MS = MMAX
      MT = MMAX - 2
      CALL SFFT (A, A, S)
      IF (MSI .EQ. -1) RETURN
c                   All setup for computation now
   20 NTOT2 = NF(NDD+1)
      IF (FN .NE. ONE) THEN
         DO 30 I = 1, NTOT2
            A(I) = A(I) * FN
   30    CONTINUE
      END IF
c     Beginning of loop for computing multiple sum
      DO 70 K = 1, NDD
         if(MODE(K:K) .eq. 'A' .or. MODE(K:K) .eq. 'a') then
            IRC = 1
            IIC = 0
         ELSE
            IRC = 0
            IIC = 1
         END IF
         MM = M(K)
         KS = NF(K)
         KSM = KS - 1
         ILAST = NF(K+1)
         DO 40 L = 1, MM
            KEE(L+1) = KEE(L) / 2
   40    CONTINUE
         DO 60 J = 1, NTOT2, ILAST
            JL = J + KSM
            DO 50 I = J, JL, 2
               IR = I + IRC
               II = I + IIC
               CALL SFFT (A(IR), A(II), S)
   50       CONTINUE
   60    CONTINUE
   70 CONTINUE
      RETURN
c        End of loop for computing multiple sum
c
      END
