subroutine mtopo(mag,ztop,nx,ny,dx,dy,mi,md,fi,fd,nstop,err, & store,cstore) c c Subroutine MTOPO calculates the total-field anomaly on a two- c dimensional grid due to an infinite half-space with uneven c top surface and two-dimensional magnetization. Method c according to Parker (Geophys. J., v. 31, p. 447-455, 1972). c Magnetization is specified on a rectangular grid with x and y c axes directed north and east, respectively. Z axis is down c and anomaly is calculated at z=0; topographic grid should be c arranged accordingly. Units of distance irrelevant but must be c consistent. Requires subroutines FOURN, DIRCOS, KVALUE, and FAC. 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 mag - a singly dimensioned real array containing the two- c dimensional magnetization (in A/m). Elements should c be in order of west to east, then south to north c (i.e., 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 ztop - a singly dimensioned real array containing the c topography of the upper surface, in same units as dx c and dy. Grid layout same as mag. Note: z axis is c positive down and anomaly is calculated at z=0. Array c ztop is modified by subroutine. c store - a singly dimensioned real array used internally. c Should be dimensioned at least 2*nx*ny. c cstore - a singly dimensioned complex array used internally. c Should be dimensioned at least nx*ny. c dx - sample interval in the x direction. c dy - sample interval in the y direction. c mi - inclination of magnetization, in degrees positive below c horizontal. c md - declination of magnetization, in degrees east of north. c fi - inclination of regional field. c fd - declination of regional field. c nstop - maximum number of iterations to try. c err - convergence criterion. Iterations stop when the c contribution of the last term to the summation is less c than err times the contribution of all previous terms. c c Output parameters: c mag - upon output, mag contains the total-field anomaly c (in nT) with same grid orientation as above. c dimension mag(nx*ny),ztop(nx*ny),store(2*nx*ny),cstore(nx*ny), & nn(2) real mag,mi,md,mx,my,mz,kx,ky,k complex cmplx,cstore,cstep,thetam,thetaf data pi/3.14159265/,t2nt/1.e9/,cm/1.e-7/ index(i,j,ncol)=(j-1)*ncol+i nn(1)=ny nn(2)=nx ndim=2 call dircos(mi,md,0.,mx,my,mz) call dircos(fi,fd,0.,fx,fy,fz) dkx=2.*pi/(nx*dx) dky=2.*pi/(ny*dy) ztpmax=-1.e20 ztpmin= 1.e20 do 1 j=1,nx do 1 i=1,ny ij=index(i,j,ny) ztpmax=amax1(ztpmax,ztop(ij)) 1 ztpmin=amin1(ztpmin,ztop(ij)) ztpmed=ztpmin+(ztpmax-ztpmin)/2. do 2 j=1,nx do 2 i=1,ny ij=index(i,j,ny) ztop(ij)=ztop(ij)-ztpmed 2 cstore(ij)=0. write(*,100) 100 format(/,' Ratio = contribution of Nth term divided by',/, & ' contribution of 0 through N-1 terms',//, & ' N | Ratio',/, & ' -----|---------') n=-1 3 n=n+1 do 4 j=1,nx do 4 i=1,ny ij=index(i,j,ny) store(2*ij-1)=mag(ij)*ztop(ij)**n 4 store(2*ij)=0. call fourn(store,nn,ndim,-1) abnew=0. abold=0. do 5 j=1,nx do 5 i=1,ny ij=index(i,j,ny) call kvalue(i,j,nx,ny,dkx,dky,kx,ky) k=sqrt(kx**2+ky**2) arg=((-k)**n)*exp(-k*ztpmed)/fac(n) cstep=arg*cmplx(store(2*ij-1),store(2*ij)) abnew=abnew+cabs(cstep) abold=abold+cabs(cstore(ij)) 5 cstore(ij)=cstore(ij)+cstep if(n.eq.0.)go to 3 ratio=abnew/abold write(*,101)n,ratio 101 format(1x,i5,g10.3) if(ratio.gt.err.and.n.lt.nstop)go to 3 do 6 j=1,nx do 6 i=1,ny ij=index(i,j,ny) if(ij.eq.1)then store(2*ij-1)=0. store(2*ij)=0. else call kvalue(i,j,nx,ny,dkx,dky,kx,ky) k=sqrt(kx**2+ky**2) thetam=cmplx(mz,(mx*kx+my*ky)/k) thetaf=cmplx(fz,(fx*kx+fy*ky)/k) cstore(ij)=2.*pi*cm*thetam*thetaf*cstore(ij) store(2*ij-1)=real(cstore(ij)) store(2*ij)=aimag(cstore(ij)) end if 6 continue call fourn(store,nn,ndim,+1) do 7 j=1,nx do 7 i=1,ny ij=index(i,j,ny) 7 mag(ij)=store(2*ij-1)*t2nt/(nx*ny) return end