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