subroutine contin(grid,nx,ny,dx,dy,dz,store)
c
c Subroutine CONTIN upward continues gridded potential field
c data using the following steps: (1) Fourier transform the field,
c (2) multiply by the continuation filter, and (3) inverse Fourier
c transform the product. Field values are specified on a
c rectangular grid with x and y axes directed north and east,
c respectively; north is arbitrary. Z axis is down. Requires
c subroutines FOURN and KVALUE.
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 a power of two.)
c grid - a singly dimensioned real array containing the
c two-dimensional potential field. 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 in the
c calling program.
c dx - sample interval in the x direction.
c dy - sample interval in the y direction.
c dz - continuation distance,in same units as dx and dy. Should
c be greater than zero for upward continuation.
c
c Output parameters:
c grid - upon output, grid contains the upward-continued
c potential field with same orientation as above.
c
dimension grid(nx*ny),store(2*nx*ny),nn(2)
real kx,ky,k
complex cgrid,cmplx
data pi/3.14159265/
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)=grid(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)
k=sqrt(kx**2+ky**2)
cont=exp(-k*dz)
cgrid=cmplx(store(2*ij-1),store(2*ij))*cont
store(2*ij-1)=real(cgrid)
20 store(2*ij)=aimag(cgrid)
call fourn(store,nn,ndim,+1)
do 30 j=1,nx
do 30 i=1,ny
ij=index(i,j,ny)
30 grid(ij)=store(2*ij-1)/(nx*ny)
return
end