      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
