SUBROUTINE CFFT2D(N,F,LDF,W,FORWD)
C***BEGIN PROLOGUE  CFFT2D
C***DATE WRITTEN   870811   (YYMMDD)
C***REVISION DATE  870811   (YYMMDD)
C***CATEGORY NO.  J1B
C***KEYWORDS  TWO DIMENSIONAL FOURIER TRANSFORM, FFT
C***AUTHOR  KAHANER, DAVID K., (NBS) 
C***PURPOSE  Two dimensional complex fast Fourier transform.
C***DESCRIPTION
C           From the book, "Numerical Methods and Software" by
C                D. Kahaner, C. Moler, S. Nash
C                Prentice Hall, 1988
C       Two dimensional fast Fourier transform, forward or backward
C          of complex N*N matrix F.
C
C   Input: 
C     N: (INTEGER)        Number of rows and columns in the matrix F to be
C                            transformed. You must set N > 0, NOT CHECKED.
C     F: (COMPLEX)        Array of N*N complex values to be transformed.
C                            This array is overwritten on output.
C     LDF: (INTEGER)       Leading (first) dimension of the complex array F 
C                             in the subroutine that calls CFFT2D. For example,
C                             if you declare F by either
C                               COMPLEX F(0:15,0:20)  OR  F(16,21) then
C                             set LDF=16. You must have LDF >= N, NOT CHECKED.  
C     W: (COMPLEX)         Array for internal use as work storage.
C                             Must be dimensioned by calling program to be
C                             at least 6N+15 REAL words or 3N+8 COMPLEX words.
C     FORWD: (LOGICAL)     Direction of transform. Set to .TRUE. for
C                            forward transform, set to .FALSE. for backward
C                            transform.
C
C
C   Output:
C     F: (COMPLEX)         Forward or reverse transformed input matrix.
C                          Output is unscaled, that is, a call to CFFT2D
C                          with FORWD=.TRUE. followed by a call to CFFT2D
C                          with FORWD=.FALSE. returns original data 
C                          multiplied by N*N.
C
C
C   Remark:
C     For some applications it is desirable to have the transform scaled so
C        the center of the N by N frequency square corresponds to zero
C        frequency. The user can do this replacing the original input data
C        F(I,J) by F(I,J)*(-1.)**(I+J),  I,J =0,...,N-1.
C
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED CFFTI, CFFTF, CFFTB
C***END PROLOGUE  CFFT2D
      COMPLEX F(0:LDF-1,0:*), W(0:*)
      LOGICAL FORWD
C***FIRST EXECUTABLE STATEMENT  CFFT2D
C   Find transform of each row, then of each column
C   First row transforms
      CALL CFFTI(N,W(N))
      DO 10 I=0,N-1
C         Place row in beginning of array W
            DO 5 J=0,N-1
               W(J) = F(I,J)
    5       CONTINUE
C         Compute fft of row
            IF(FORWD) THEN
               CALL CFFTF(N,W,W(N))
            ELSE
               CALL CFFTB(N,W,W(N))
            ENDIF
C         Copy back to row of f
            DO 6 J=0,N-1
               F(I,J)=W(J)
    6       CONTINUE          
   10 CONTINUE
C   Column transforms
      DO 20 J=0,N-1
C         Pass column of F to CFFTF or CFFTB by passing F(0,J). 
C         Similarly CFFTF or CFFTB places results back there.
            IF(FORWD) THEN
               CALL CFFTF(N, F(0,J), W(N))
            ELSE
               CALL CFFTB(N, F(0,J), W(N))
            ENDIF
   20 CONTINUE
      RETURN
      END