	PROGRAM proj

	
	DIMENSION grid(4096),store(8192)
	CHARACTER*64 File1,File2	

	WRITE(*,*) "Enter input file name"
	READ(*,*) File1
	WRITE(*,*) "ENter output file name"
	READ(*,*) File2
	WRITE(*,*) "ENTER nx value"
	READ(*,*) nxx
	WRITE(*,*) "ENter ny value"
	READ(*,*) nyx
	WRITE(*,*) "enter dx value"
	READ(*,*) dxx
	WRITE(*,*) "enter dy value"
	READ(*,*) dyx


	OPEN(UNIT=1,file=File1,status='unknown')
	OPEN(UNIT=8,file=File2,status='unknown')


 01	WRITE(*,*) "which function do you wish to perform?"
	WRITE(*,*) "1)upper contin"

	READ(*,*) n

	IF (n.gt.1) THEN
	WRITE(*,*) "not valid"
	GOTO 01
	ENDIF

	IF (n.eq.1) THEN
	 WRITE(*,*) "enter contin distance"
	 READ(*,*) dzx
	ENDIF

	READ(1,*) (grid(i), i=1,4096)
		
	CALL contin(grid,nxx,nyx,dxx,dyx,dzx,store)
	
	
	count=1
	DO j=1,nyx
	 DO i=1,nxx
	 WRITE(8,*) (i-1)*dxx, (j-1)*dyx, grid(count)
	 count=count+1
	 END DO
	END DO
	
	CLOSE(1)
	CLOSE(8,status='keep')
	END


      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


      subroutine fourn(data,nn,ndim,isign)
c
c  Replaces DATA by its NDIM-dimensional discrete Fourier transform, 
c  if ISIGN is input as 1.  NN is an integer array of length NDIM, 
c  containing the lengths of each dimension (number of complex values), 
c  which must all be powers of 2.  DATA is a real array of length twice 
c  the product of these lengths, in which the data are stored as in a 
c  multidimensional complex Fortran array.  If ISIGN is input as -1, 
c  DATA is replaced by its inverse transform times the product of the
c  lengths of all dimensions.  From Press, W.H., Flannery, B.P., 
c  Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical Recipes, 
c  Cambridge Univ. Press, p. 451-453.
c
      real*8 wr,wi,wpr,wpi,wtemp,theta
      dimension nn(ndim),data(*)
      ntot=1
      do 11 iidim=1,ndim
        ntot=ntot*nn(iidim)
11    continue
      nprev=1
      do 18 iidim=1,ndim
        n=nn(iidim)
        nrem=ntot/(n*nprev)
        ip1=2*nprev
        ip2=ip1*n
        ip3=ip2*nrem
        i2rev=1
        do 14 i2=1,ip2,ip1
          if(i2.lt.i2rev)then
            do 13 i1=i2,i2+ip1-2,2
              do 12 i3=i1,ip3,ip2
                i3rev=i2rev+i3-i2
                tempr=data(i3)
                tempi=data(i3+1)
                data(i3)=data(i3rev)
                data(i3+1)=data(i3rev+1)
                data(i3rev)=tempr
                data(i3rev+1)=tempi
12            continue
13          continue
          endif
          ibit=ip2/2
1         if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
            i2rev=i2rev-ibit
            ibit=ibit/2
          go to 1
          endif
          i2rev=i2rev+ibit
14      continue
        ifp1=ip1
2       if(ifp1.lt.ip2)then
          ifp2=2*ifp1
          theta=isign*6.28318530717959d0/(ifp2/ip1)
          wpr=-2.d0*dsin(0.5d0*theta)**2
          wpi=dsin(theta)
          wr=1.d0
          wi=0.d0
          do 17 i3=1,ifp1,ip1
            do 16 i1=i3,i3+ip1-2,2
              do 15 i2=i1,ip3,ifp2
                k1=i2
                k2=k1+ifp1
                tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                data(k2)=data(k1)-tempr
                data(k2+1)=data(k1+1)-tempi
                data(k1)=data(k1)+tempr
                data(k1+1)=data(k1+1)+tempi
15            continue
16          continue
            wtemp=wr
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wtemp*wpi+wi
17        continue
          ifp1=ifp2
        go to 2
        endif
        nprev=n*nprev
18    continue
      return
      end


      subroutine kvalue(i,j,nx,ny,dkx,dky,kx,ky)
c  Subroutine KVALUE finds the wavenumber coordinates of one 
c  element of a rectangular grid from subroutine FOURN.  
c
c  Input parameters:
c    i  - index in the ky direction.
c    j  - index in the kx direction.
c    nx - dimension of grid in ky direction (a power of two).
c    ny - dimension of grid in kx direction (a power of two).
c    dkx - sample interval in the kx direction.
c    dky - sample interval in the ky direction.
c
c  Output parameters:
c    kx - the wavenumber coordinate in the kx direction.
c    ky - the wavenumber coordinate in the ky direction.
c
      real kx,ky
      nyqx=nx/2+1
      nyqy=ny/2+1
      if(j.le.nyqx)then
         kx=(j-1)*dkx
         else
            kx=(j-nx-1)*dkx
            end if
      if(i.le.nyqy)then
         ky=(i-1)*dky
         else
            ky=(i-ny-1)*dky
            end if
      return
      end
