subroutine mlayer(mag,nx,ny,dx,dy,z1,z2,mi,md,fi,fd,store)
c
c  Subroutine MLAYER calculates the total-field anomaly on a two-
c  dimensional grid due to a horizontal layer with two-
c  dimensional magnetization.  The following steps are involved:
c  (1) Fourier transform the magnetization, (2) multiply by the  
c  earth filter, and (3) inverse Fourier transform the product.   
c  Magnetization is specified on a rectangular grid with x and y  
c  axes directed north and east, respectively.  Z axis is down.  
c  Distance units irrelevant but must be consistent.  Requires 
c  subroutines FOURN, DIRCOS, KVALUE, and MFILT.  
c 
c  Input parameters:
c    nx - number of elements in the sout_to_north direction.  
c    ny - number of elements in the west-to-east direction.  
c         (NOTE:  both nx and ny must be a power of two.)
c    mag - a singly dimensioned real array containing the 
c         two-dimensional magnetization (in A/m).  Elements should 
c         be in order of west to east, then south to north (i.e., 
c         element 1 is southwest corner, element ny is 
c         southeast corner, element (nx-1)*ny+1 is northwest 
c         corner, and element ny*nx is northeast corner. 
c    store - a singly dimensioned real array used internally.  
c         It should be dimensioned at least 2*nx*ny. 
c    dx - sample interval in the x direction.
c    dy - sample interval in the y direction. 
c    z1 - depth to top of layer.  Must be > 0.
c    z2 - depth to bottom of layer.  Must be > z1.
c    mi - inclination of magnetization, in degrees positive below 
c         horizontal.
c    md - declination of magnetization, in degrees east of north.
c    fi - inclination of regional field.
c    fd - declination of regional field.
c
c  Output parameters:
c    mag - upon output, mag contains the total-field anomaly
c         (in nT) with same orientation as above.
c
      complex cmag,mfilt,cmplx
      real mag,mi,md,mx,my,mz,kx,ky
      dimension mag(nx*ny),store(2*nx*ny),nn(2)
      data pi/3.14159265/,t2nt/1.e9/
      index(i,j,ncol)=(j-1)*ncol+i
      nn(1)=ny
      nn(2)=nx
      ndim=2
      call dircos(mi,md,0.,mx,my,mz)
      call dircos(fi,fd,0.,fx,fy,fz)
      dkx=2.*pi/(nx*dx)
      dky=2.*pi/(ny*dy)
      do 10 j=1,nx
         do 10 i=1,ny
            ij=index(i,j,ny)
            store(2*ij-1)=mag(ij)
   10       store(2*ij)=0.
      call fourn(store,nn,ndim,-1)
      do 20 j=1,nx
         do 20 i=1,ny
            ij=index(i,j,ny)
            call kvalue(i,j,nx,ny,dkx,dky,kx,ky)
            cmag=cmplx(store(2*ij-1),store(2*ij))
            cmag=cmag*mfilt(kx,ky,mx,my,mz,fx,fy,fz,z1,z2)
            store(2*ij-1)=real(cmag)
   20       store(2*ij)=aimag(cmag)
      call fourn(store,nn,ndim,+1)
      do 30 j=1,nx
         do 30 i=1,ny
            ij=index(i,j,ny)
   30       mag(ij)=store(2*ij-1)*t2nt/(nx*ny)
      return
      end