subroutine cfft2d(a,wk,n,iop)
c
c    compute the 2-d fft transform! Must call cffti(n,wk) to initiate !
c    iop=-1: fft
c    iop=1:  inverse fft
c
c     n  <= 512

      double complex a(n,n),wk(1),b(512)
      if (iop.ne.1) then
      do 1 i=1,n
         call cfftf(n,a(1,i),wk)
 1    continue
      do 2 i=1,n
         do 5 j=1,n
            b(j)=a(i,j)
 5       continue
         call cfftf(n,b,wk)
         do 6 j=1,n
            a(i,j)=b(j)
6        continue
 2    continue
      else
      do 10 i=1,n
         call cfftb(n,a(1,i),wk)
 10   continue
      do 20 i=1,n
         do 50 j=1,n
            b(j)=a(i,j)
 50      continue
         call cfftb(n,b,wk)
         do 60 j=1,n
            a(i,j)=b(j)
 60      continue
 20   continue
      endif
      return
      end