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