subroutine glayer(rho,nx,ny,dx,dy,z1,z2,store)
c
c  Subroutine GLAYER calculates the vertical gravitational 
c  attraction on a two-dimensional grid caused by a two-
c  dimensional density confined to a horizontal layer.  The 
c  following steps are involved:  (1) Fourier transform the 
c  density, (2) multiply by the earth filter, and (3) inverse 
c  Fourier transform the product.  Density is specified on a 
c  rectangular grid with x and y axes directed north and east,
c  respectively.  Z axis is down.  Requires subroutines FOURN, 
c  KVALUE, and GFILT.  
c
c  Input parameters:
c    nx - number of elements in the south-to-north direction.  
c    ny - number of elements in the west-to-east direction.  
c         (NOTE: Both nx and ny must be powers of two.
c    rho - a singly dimensioned real array containing the 
c         two-dimensional density, in kg/(m**3).  Elements 
c         should be in order of west to east, then south to 
c         north (i.e., element 1 is the southwest corner, 
c         element ny is the southeast corner, element 
c         (nx-1)*ny+1 is the northwest corner, and element ny*nx
c         is the 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, in km.
c    dy - sample interval in the y direction, in km. 
c    z1 - depth to top of layer, in km.  Must be > 0.
c    z2 - depth to bottom of layer, in km.  Must be > z1.
c
c  Output parameters:
c    rho - upon output, rho will contain the gravity anomaly,
c         in mGal, with same orientation as above.
c
      complex crho,cmplx
      real kx,ky,km2m
      dimension rho(nx*ny),store(2*nx*ny),nn(2)
      data pi/3.14159265/,si2mg/1.e5/,km2m/1.e3/
      index(i,j,ncol)=(j-1)*ncol+i
      nn(1)=ny
      nn(2)=nx
      ndim=2
      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)=rho(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)
            crho=cmplx(store(2*ij-1),store(2*ij))
            crho=crho*gfilt(kx,ky,z1,z2)
            store(2*ij-1)=real(crho)
   20       store(2*ij)=aimag(crho)
      call fourn(store,nn,ndim,+1)
      do 30 j=1,nx
         do 30 i=1,ny
            ij=index(i,j,ny)
   30       rho(ij)=store(2*ij-1)*si2mg*km2m/(nx*ny)
      return
      end