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