      SUBROUTINE DCFT2D(N,F,LDF,W,FORWD)
C***BEGIN PROLOGUE  DCFT2D
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*16)        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 DCFT2D. For example,
C                             if you declare F by either
C                               COMPLEX*16 F(0:15,0:20)  OR  F(16,21) then
C                             set LDF=16. You must have LDF >= N, NOT CHECKED.  
C     W: (COMPLEX*16)         Array for internal use as work storage.
C                             Must be dimensioned by calling program to be
C                             at least 6N+15 D. P. words or 3N+8 COMPLEX*16 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*16)         Forward or reverse transformed input matrix.
C                          Output is unscaled, that is, a call to DCFT2D
C                          with FORWD=.TRUE. followed by a call to DCFT2D
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.0D0
C
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED DCFFTI, DCFFTF, DCFFTB
C***END PROLOGUE  DCFT2D
      COMPLEX*16 F(0:LDF-1,0:*), W(0:*)
      LOGICAL FORWD
C
C  NOTE:  DOUBLE PRECISION COMPLEX IS NOT STANDARD FORTRAN.
C         SOME COMPILER VENDORS MAY FAIL TO IMPLEMENT IT
C         OR MAY IMPLEMENT IN A DIFFERENT MANNER.  IF YOUR
C         COMPILER DOES NOT ACCEPT THE COMPLEX*16 STATEMENT
C         WE RECOMMEND THAT YOU CHECK WITH YOUR VENDOR AND
C         MODIFY THE ABOVE ACCORDINGLY.
C
C***FIRST EXECUTABLE STATEMENT  DCFT2D
C   Find transform of each row, then of each column
C   First row transforms
      CALL DCFFTI(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 DCFFTF(N,W,W(N))
            ELSE
               CALL DCFFTB(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 DCFFTF or DCFFTB by passing F(0,J). 
C         Similarly DCFFTF or DCFFTB places results back there.
            IF(FORWD) THEN
               CALL DCFFTF(N, F(0,J), W(N))
            ELSE
               CALL DCFFTB(N, F(0,J), W(N))
            ENDIF
   20 CONTINUE
      RETURN
      END   

