      subroutine sfftmg(unita,unitb,n,m,f[huge],l,isign,work)
c     performs 2-dimensional fast fourier transform using disk
c     arguments -
c       unita - unit number of open random access disk file containing
c               one complex data row per record.
c               this file is used for input to fftmg and output
c               of the transform from fftmg.
c       unitb - unit number of open random access work file
c               the same as the file of unita.
c       n     - number of y columns = length of x rows.
c       m     - number of x rows = length of y columns,
c               must be of the form l*2**k (k integer).
c       f     - complex work array for core manipulations.
c       l     - number of x rows which can be held in f.
c       isign - +1 or -1 for forward or reverse transform.
c       work  - work array for in-core transform routine sfourt.
c               work must be complex, of length max(n,l).
c   developed and coded by ray watts--usgs,denver
c   note: randin,randou,inset, and outset are entry points
c   in subroutine sdummy
c***********************************************************************
      integer unita,unitb,output,ndim(2)
      complex f(n,l),twidle,twidlf,w,work(1)
      nfourt=2*n*l
      nwork=2*n
      if(l.gt.n)nwork=2*l
c     **********************************************************
c     *    the number of disk accesses (reads+writes) used by  *
c     *    this program is m*(4*log2(m/l)+3).                  *
c     **********************************************************
c     --------------------------------------------------------
c     check arguments
      if(l*(m/l).ne.m)stop 'l not a factor of m in fftmg'
      ml=m/l
      mt=1
   10 if(mt-ml)20,40,30
   20 mt=2*mt
      go to 10
   30 stop 'm/l not a power of 2 in fftmg'
   40 continue
c     ------------------------------------------------------------
c     -------------------------------------------------------------
c     start processing
c     skip decimation if l .eq. m
      if(l.eq.m)go to 80
c     set starting input and output units
      input=unita
      output=unitb
c     set starting blocksize for decimation in input domain
      iblock=m
c     -------------------------------------------------------
c     decimate this blocksize
c     set up sub-block size
   50 isub=iblock/2
c     set up for random i/0
c      call inset(input)
c      call outset(output)
c     for each block of this size
      do 60 iblst=1,m,iblock
c     copy even members into first sub-block, odd members into second
      do 60 i=1,isub 
c      call randin(n,f,2*(i-1)+iblst)
c      call randou(n,f,i-1+iblst)
c      call randin(n,f,2*i-1+iblst)
c   60 call randou(n,f,i-1+iblst+isub)
        read (input,rec=2*(i-1)+iblst) (f(k,1),k=1,n)
        write (output,rec=i-1+iblst) (f(k,1),k=1,n)
        read(input,rec=2*i-1+iblst) (f(k,1),k=1,n)
60      write (output,rec=i-1+iblst+isub) (f(k,1),k=1,n)
c     swap input/output functions
      inputt=input
      input=output
      output=inputt
c     next smaller blocksize
      iblock=isub
c     loop if not done
      if(iblock.gt.l)go to 50
c     decimation complete
c     ---------------------------------------------------------
c     ---------------------------------------------------------
c     ensure that data set is back on unita
      if(input.eq.unita)go to 80
c      call inset(input)
c      call outset(output)
      do 70 i=1,m
c      call randin(n,f,i)
       read (input,rec=i) (f(k,1),k=1,n)
c   70 call randou(n,f,i)
   70  write (output,rec=i) (f(k,1),k=1,n)
   80 continue
c     --------------------------------------------------------
c     start transform.
c     set block size
      iblock=l
c     set up for random i/o 
c      call inset(unita)
c      call outset(unita)
c     run through the blocks
      do 100 istart=1,m,iblock
      into=1
c     run through the records in a block
      do 90 irecrd=istart,istart+iblock-1
c     read the record
c      call randin(n,f(1,into),irecrd)
        read(unita,rec=irecrd) (f(k,into),k=1,n)
c     proceed with next record in block
   90 into=into+1
c     fourier transform in both directions
      ndim(1)=n
      ndim(2)=l  
      call sfourt(f,ndim,2,isign,1,work,nfourt,nwork) 
c     ******************************************************
c     write the block back onto disk
      ifrom=1
      do 100 irecrd=istart,istart+iblock-1
c      call randou(n,f(1,ifrom),irecrd)
        write (unita,rec=irecrd) (f(k,ifrom),k=1,n)
  100 ifrom=ifrom+1
c     if l.eq.m, the transform is now complete.
      if(l.eq.m)return
c     done with first-step processing
c     ---------------------------------------------------------
c     finish the processing by pairs of rows, applying twiddle factors
c     set the twiddle generator
      twidlf=cexp(cmplx(0.,3.141593*isign/(iblock)))
c     run through the blocks by pairs
  110 do 140 istart=1,m,iblock*2
c     start the twiddle factor
      twidle=(1.,0.)
c     scan the block pair
      do 130 irecrd=istart,istart+iblock-1
      jrecrd=irecrd+iblock
c      call randin(n,f,irecrd)
        read(unita,rec=irecrd) (f(k,1),k=1,n)
c      call randin(n,f(1,2),jrecrd)
        read(unita,rec=jrecrd) (f(k,2),k=1,n)
c     combine, using twiddle factor
      do 120 i=1,n
      w=twidle*f(i,2)
      f(i,2)=f(i,1)-w
  120 f(i,1)=f(i,1)+w
c     write back onto disk
c      call randou(n,f,irecrd)
        write (unita,rec=irecrd)(f(k,1),k=1,n)
c      call randou(n,f(1,2),jrecrd)
        write (unita,rec=jrecrd)(f(k,2),k=1,n)
c     increment the twiddle factor
  130 twidle=twidle*twidlf
c     go on to next block pair
  140 continue
c     increase the block size, change twiddle increment
      iblock=2*iblock
      twidlf=csqrt(twidlf)
      if(isign*aimag(twidlf).lt.0)twidlf=-twidlf
c     check if done
      if(iblock.lt.m)go to 110
c     ----------------------------------------------------------
c     done processing
      end
c***************************************************************
        subroutine sdummy
c       random access i/o routine
c***************************************************************
        integer out
        complex f(n),g(n)
c       entry to receive input setup
        entry inset(iin)
        in=iin
        return
c       entry to receive output setup
        entry outset(iout)
        out=iout
        return
c       entry to do input
        entry randin(n,f[huge],nrec)
        read (in,rec=nrec) f
        return
c       entry to do output
        entry randout(n,g[huge],nrec)
        write (out,rec=nrec) g
        return
        end




