c   PMAG3D
c
c   computes the magnetic anomaly due to a magnetic layer.
c   The layer can be specified by up to three standard
c   grid files: the top surface, the bottom surface, and
c   the magnetization.  Alternatively up to two of these
c   grids may be replaced by constants.  Each
c   input grid must be complete (no dvals), and its
c   column dimension must be less than 2049.
c
c   link with SFFTMG and SFOURT
c
c   Jeff Phillips April 1993
c
      dimension a1(2,2048,16),work(4096)
      character*80 ifile,ofilet,ofileb,ofilem
      character*1 ans,drive
      common/mane/id(14),pgm(2),dx,dy,nc,nr
      common/top/topscale,tconst,tmin,tmax
      common/bottom/botscale,bconst,bmin,bmax
      common/other/cmag,sincl,sdecl,fincl,fdecl,nstop,err2,zh
      common/fast/drive,ofilet,ofileb,ofilem
      data iw/6/
      data dval/1.e37/
      data eps/.1e-10/
      print'(a\)',' Enter the (RAM) drive letter for temporary files: '
      read(5,1) drive
c get level of observation surface
      print*,'Enter the survey elevation in units of dx and dy;'
      print'(a\)','   remember z is positive down: '
      read(5,*) zh
c get top surface
      topscale=0.
      tconst=0.
      print'(a\)',' Is upper surface flat ? '
      read(5,1) ans
      if(ans.eq.'n'.or.ans.eq.'N') then
c upper surface is specified as a grid
        print'(a\)',' Standard file containing gridded upper surface: '
        read(5,1)ifile
1       format(a)
        open(10,file=ifile,status='old',form='unformatted')
        read(10)id,pgm,nc,nr,nz,xo,dx,yo,dy
        nctest=nc
        nrtest=nr
c get tmin,tmax and test for dvals
        tmin=1.e37
        tmax=-1.e37
        do 20 j=1,nr
        call rowio(work,yr,nz,nc,10,1)
        do 20 i=1,nc
        if(work(i).ge.dval) stop 'dvals found'
        if(work(i).lt.tmin) tmin=work(i)
        if(work(i).gt.tmax) tmax=work(i)
   20   continue
c get scale factor for upper surface
        print*,'grid dx/dy = ',dx,dy
        print*,'surface min/max = ',tmin,tmax
        print*,'What single number should these values be multiplied by'
        print*,'  so that they increase vertically down and have the'
        print'(a\)','   same units as dx?: '
        read(5,*) topscale
        tmin=tmin*topscale
        tmax=tmax*topscale
c get dimensions of expanded upper surface grid
        call check_dims(nc,nr,ncol,nrow)
      else
c get constant upper surface level
11      print'(a\)',' Elevation: '
        read(5,*) tconst
        if(tconst.lt.(zh+eps)) then
          print*,'Warning: upper surface lies above the survey'
          print*,'  elevation!  Try again...'
          go to 11
        endif
        tmin=tconst
        tmax=tconst
      endif
c get bottom surface
      botscale=0.
      bconst=0.
      print'(a\)',' Is bottom surface flat ? '
      read(5,1) ans
      if(ans.eq.'n'.or.ans.eq.'N') then
c bottom surface is specified as a grid
        print'(a\)',' Standard file containing gridded bottom surface: '
        read(5,1)ifile
        open(11,file=ifile,status='old',form='unformatted')
        read(11)id,pgm,nc,nr,nz,xo,dx,yo,dy
c either test for registration with upper surface
        if(topscale.ne.0.) then
          if(nc.ne.nctest.or.nr.ne.nrtest) stop 'grids do not overlay'
        else
          nctest=nc
          nrtest=nr
        endif
c get bmin,bmax and test for dvals
        bmin=1.e37
        bmax=-1.e37
        do 21 j=1,nr
        call rowio(work,yr,nz,nc,11,1)
        do 21 i=1,nc
        if(work(i).ge.dval) stop 'dvals found'
        if(work(i).lt.bmin) bmin=work(i)
        if(work(i).gt.bmax) bmax=work(i)
   21   continue
c get scale factor for bottom surface
        print*,'grid dx/dy = ',dx,dy
        print*,'surface min/max = ',bmin,bmax
        print*,'What single number should these values be multiplied by'
        print*,'  so that they increase vertically down and have the'
        print'(a\)','   same units as dx?:'
        read(5,*) botscale
        bmin=bmin*botscale
        bmax=bmax*botscale
c get dimensions of expanded bottom surface grid
        if(topscale.eq.0.0) call check_dims(nc,nr,ncol,nrow)
      else
c get a constant bottom surface level
22      print'(a\)',' Elevation: '
        read(5,*) bconst
        if(bconst.lt.tmax) then
          print*,'Warning: the layer has negative thickness at some'
          print*,'  point!  Try again...'
          go to 22
        endif
        bmin=bconst
        bmax=bconst
      endif
c check surfaces for negative thickness or sources above observations
      if(topscale.ne.0) then
        rewind(10)
        open(12,form='unformatted',status='scratch')
        read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
        write(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
        tmax=-1.e37
        tmin=1.e37
      else
        tmax=tconst
        tmin=tconst
      endif
      if(botscale.ne.0) then
        rewind(11)
        open(13,form='unformatted',status='scratch')
        read(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
        write(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
        bmax=-1.e37
        bmin=1.e37
      else
        bmax=bconst
        bmin=bconst
      endif
      nerr1=0
      nerr2=0
      do j=1,nr
      if(topscale.ne.0) then
        call rowio(work,yr,nz,nc,10,1)
        do i=1,nc
          work(i)=work(i)*topscale
          if(work(i).le.(zh+eps)) then
            nerr2=nerr2+1
            work(i)=zh+eps
          endif
        enddo
        if(botscale.ne.0) then
          call rowio(work(2049),yr,nz,nc,11,1)
          do i=1,nc
            work(2048+i)=work(2048+i)*botscale
            if(work(i).gt.work(2048+i)) then
              nerr1=nerr1+1
              work(2048+i)=work(i)
            endif
            if(work(2048+i).lt.bmin) bmin=work(2048+i)
            if(work(2048+i).gt.bmax) bmax=work(2048+i)
          enddo
          call rowio(work(2049),yr,nz,nc,13,2)
        else
          do i=1,nc
            if(work(i).gt.bconst) then
              nerr2=nerr2+1
              work(i)=bconst
            endif
          enddo
        endif
        do i=1,nc
          if(work(i).lt.tmin) tmin=work(i)
          if(work(i).gt.tmax) tmax=work(i)
        enddo
        call rowio(work,yr,nz,nc,12,2)
      else
        if(botscale.ne.0) then
          call rowio(work,yr,nz,nc,11,1)
          do i=1,nc
            work(i)=work(i)*botscale
            if(work(i).lt.tconst) then
              nerr1=nerr1+1
              work(i)=tconst
            endif
            if(work(i).lt.bmin) bmin=work(i)
            if(work(i).gt.bmax) bmax=work(i)
          enddo
          call rowio(work,yr,nz,nc,13,2)
        endif
      endif
      enddo
      close(10)
      close(11)
      if(nerr2.ne.0) then
        print*,'Warning: the layer rises above the survey elevation'
        print*,'  at ',nerr2,' points.  The top of the layer will be'
        print*,'  truncated accordingly.'
      endif
      if(nerr1.ne.0) then
        print*,'Warning: the layer has negative thickness at ',nerr1
        print*,'  points.  The body will be assigned zero thickness'
        print*,'  at these locations.'
      endif
c extend the grids
      if(topscale.ne.0.0) then
        ofilet=drive//':\pmag.50'
        call prep(12,ofilet,topscale,ncol,nrow,a1,work,work(2049))
        print *,'Extended upper surface grid written to ',ofilet
      endif
      if(botscale.ne.0.0) then
        ofileb=drive//':\pmag.51'
        call prep(13,ofileb,botscale,ncol,nrow,a1,work,work(2049))
        print *,'Extended lower surface grid written to ',ofileb
      endif
c get magnetization
      cmag=0.
      print'(a\)',' Use a constant magnetization? '
      read(5,1) ans
      if(ans.eq.'n'.or.ans.eq.'N') then
c the magnetization is specified as a grid
        print'(a\)',' Standard file containing gridded magnetization: '
        read(5,1)ifile
        open(10,file=ifile,status='old',form='unformatted')
        read(10)id,pgm,nc,nr,nz,xo,dx,yo,dy
        if(topscale.ne.0.0.or.botscale.ne.0.0) then
          if(nc.ne.nctest.or.nr.ne.nrtest) stop 'grids do not overlay'
        else
          call check_dims(nc,nr,ncol,nrow)
        endif
        ofilem=drive//':\pmag.52'
        call prep(10,ofilem,1.0,ncol,nrow,a1,work,work(2049))
        print *,'Extended magnetization grid written to ',ofilem
        close(10)
      else
c the magnetization is a constant
        if(topscale.eq.0.0.and.botscale.eq.0.0) then
          stop 'at least one grid is required'
        endif
        print'(a\)',' What is the value of the magnetization: '
        read(5,*) cmag
      endif
c get other parameters
      print*
      print*,'  sincl = Inclination of magnetization in degrees'
      print*,'  sdecl = Declination of magnetization in degrees'
      print*,'  fincl = Inclination of regional field in degrees'
      print*,'  fdecl = Declination of regional field in degrees'
      print*,'  nstop = Limit on number of summations (eg. 20)'
      print*,'  err2  = Convergence criterion for summations (eg. 0.001
     1)'
      print*
      print'(a\)',' Enter sincl,sdecl,fincl,fdecl,nstop,err2 :'
      read(5,*) sincl,sdecl,fincl,fdecl,nstop,err2
      if(nstop.gt.31) nstop=31
c get output file and set up calculations
      print'(a\)',' Standard file to recieve gridded anomaly:'
      read(5,1) ifile
      open(10,file=ifile,form='unformatted',status='unknown')
      pgm='pmag3d'
      write(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
      call get_dims(nrow,ncol,nx,ny,nri,id2,nxa)
      write(iw,275) nri
  275 format(' blocking factor for rows =',i3)
      call fmag3d(a1,work,nri,nrow,ncol,id2)
      stop
      end
c **************************************************************
      subroutine check_dims(nc,nr,ncol,nrow)
        print*,nc,nr
10      call fftdims(nc,nr)
        print*, 'Enter (new) ncol and nrow:'
        read*, ncol,nrow
c  Check to make sure SFFTMG will accept these.
c  Added 19 June 91.
c        kr=11
c        iw=6
c        ir=5
c        nadd=0
c  note that x & y have been switched from normal (usgs) grid specficati
c  ...convention used is x-north y-east and z-down
        n1=nrow
        n2=ncol
        call get_dims(n1,n2,nx,ny,nri,id2,nxa)
        itest=n1-nx
        if(itest.ne.0) then
          print*,'Select another nrow; SFFTMG cannot handle this one.'
          go to 10
        endif
        if(nc.gt.ncol.or.nr.gt.nrow)then
          print*,'New ncol, nrow are too small.'
          go to 10
        endif
        return
        end

c  read header record of extended input file.
c  180 open(kr,status='old',form='unformatted',file='n:\mfilt.50')
c  note that x & y have been switched from normal (usgs) grid specficati
c  ...convention used is x-north y-east and z-down
c      read(kr) id,pgm,n2,n1,nz,yo,dy,xo,dx
c      pgm(1)='mfin'
c      pgm(2)='it  '
c      call get_dims(n1,n2,nx,ny,nri,id2,nxa)
c      write(iw,275) nri
c  275 format(' blocking factor for rows =',i3)
c      write(14) nri,n1,n2,id2
c      close(14)
c  input parameters for fft are set: call main program.
c      call mfftfil(a1,work,nri,n1,n2,id2)
c      stop
c      end
c **************************************************************
      subroutine prep(iunit,ofile,scale,ncol,nrow,g,g1,g2)
c
c  Open grid and cosine taper to reduce FFT wraparound.
c  The new nrow, ncol should be chosen from the list of rich-in-
c  -factors-of-2 numbers that SFFTMG likes.
c
c  Lin Cordell Jan 1990.
c  Revised June 1991.
c  Modified by Jeff Phillips Feb 1992.
c
        character*80 ofile
      dimension g(1),g1(1),g2(1),id(14),pgm(2)
c      index(i,j,nc) = ((j - 1) * nc) + i
        pi=3.1415927
c  Tilt grid.
c      call edge(alpha,beta,gamma,nc,nr,10)
c      close(10)
c        print*, ' alpha = ',alpha,', beta = ',beta,' ,gamma = ',gamma
c        open(10,file=ifile,status='old',form='unformatted')
        rewind(iunit)
        read(iunit)id,pgm,nc,nr,nz,xo,dx,yo,dy
c  Output extended file.
      open(14,file=ofile,status='unknown',form='unformatted')
	pgm(1)='prep'
	pgm(2)='    '
        write(14)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
c      open(14,file='mfilt.6',status='unknown',form='unformatted')
c        write(14)id,pgm,nc,nr,nz,xo,dx,yo,dy
c        write(14)alpha,beta,gamma,ncol,nrow
c  First row (saved).
        j2=1
        call rowio(g,yr,nz,nc,iunit,1)
        y=float(j2)
      do 211 j1=1,nc
        x=float(j1)
c        g(j1)=(g(j1)-alpha*x-beta*y-gamma)*scale
      g(j1)=g(j1)*scale
211   g1(j1)=g(j1)
c  Cosine taper right side
        s=0.5*(g(nc)+g(1))
	r=0.5*(g(nc)-g(1))
	lim=ncol-nc
	rlamb=float(lim+1)
	arg=pi/rlamb
	do 201 k=1,lim
	rk=float(k)
	ge=s+r*cos(arg*rk)
	g(nc+k)=ge
201     g1(nc+k)=ge
        call rowio(g,yr,nz,ncol,14,2)
c  Rows down to nr
        do 204 j2=2,nr
        y=float(j2)
        call rowio(g,yr,nz,nc,iunit,1)
	do 202 j1=1,nc
        x=float(j1)
c202     g(j1)=g(j1)-alpha*x-beta*y-gamma
202     g(j1)=g(j1)*scale
c  Cosine taper right side.
        s=0.5*(g(nc)+g(1))
	r=0.5*(g(nc)-g(1))
	do 203 k=1,lim
	rk=float(k)
203     g(nc+k)=s+r*cos(arg*rk)
        call rowio(g,yr,nz,ncol,14,2)
204	continue
c  Last row
c  Extend right side
c  Rows nr+1 to nrow.
        do 205 j1=1,ncol
205     g2(j1)=g(j1)
        lim=nrow-nr
	rlamb=float(lim+1)
	arg=pi/rlamb
	do 207 k=1,lim
	yr=float(nr+k-1)*dx+yo
	rk=float(k)
	do 206 j1=1,ncol
	s=0.5*(g2(j1)+g1(j1))
	r=0.5*(g2(j1)-g1(j1))
206     g(j1)=s+r*cos(arg*rk)
        call rowio(g,yr,nz,ncol,14,2)
207	continue
	go to 999
900     print*, "Can't handle it."
        close(iunit)
999     close(14)
        return
	end
c **************************************************************
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
c **************************************************************
      subroutine fftdims(ncol, nrow)
c
c     suggests new dimensions for fourier transform
c
c     Jeff Phillips 1992
c
      dimension nlist(54), ngood(54)
      data nlist / 9, 10, 11, 12, 13, 16
     &, 18, 20, 22, 24, 26, 32
     &, 36, 40, 44, 48, 52, 56, 64
     &, 72, 80, 88, 96, 104, 112, 128
     &, 144, 160, 176, 192, 208, 224, 256
     &, 288, 320, 352, 384, 416, 448, 512
     &, 576, 640, 704, 768, 832, 896, 1024
     &, 1152, 1280, 1408, 1536, 1664, 1792, 2048 /
c
      write(unit=6, fmt=901) ncol, nrow
  901 format(' ncol = ',i5,' nrow = ',i5)
      n = 0
      do 10 i = 1, 54
      if (nlist(i) .lt. ncol) goto 10
      if (nlist(i) .gt. ncol+ncol/2) goto 20
      n = n + 1
      ngood(n) = nlist(i)
   10 continue
   20 write(unit=6, fmt=902)
  902 format(' Suggested new values:')
      write(unit=6, fmt=101) (ngood(i),i = 1, n)
  101 format('  ncol = ',14i5)
      n = 0
      do 30 i = 1, 54
      if (nlist(i) .lt. nrow) goto 30
      if (nlist(i) .gt. nrow+nrow/2) goto 40
      n = n + 1
      ngood(n) = nlist(i)
   30 continue
   40   write(unit=6, fmt=102) (ngood(i),i = 1, n)
  102 format('  nrow = ',14i5)
      return 
      end
c **************************************************************
      subroutine fmag3d(a1,work[huge],nri,nrow,ncol,id2)
      dimension a1(2,ncol,nri),work(id2)
      dimension fact(31),s(2048),ztop(2048),zbot(2048)
      real mx,my,mz
      character*80 ofilet,ofileb,ofilem
      character*1 drive
      common/top/topscale,tconst,tmin,tmax
      common/bottom/botscale,bconst,bmin,bmax
      common/other/cmag,sincl,sdecl,fincl,fdecl,nstop,err2,zh
      common/mane/id(14),pgm(2),dx,dy,nc,nr
      common/fast/drive,ofilet,ofileb,ofilem
      data twopi / 6.2831853 /
      data conv / .01745329 /
c      data epslon / .1e-10 /
      data iunita/12/,iunitb/13/,iunitc/14/,iunitd/15/
c initialize some parameters
      do 22 i=1,nstop
      fact(i)=fac(i,ierr)
   22 continue
      ntotal=nrow*ncol
      tmed = tmin + ((tmax - tmin) / 2.)
      bmed = bmin + ((bmax - bmin) / 2.)
      write(6,108) tmin, tmax, tmed, bmin, bmax, bmed
  108 format(/,40h THE HIGH, LOW, AND MEDIAN POINTS OF THE
     &16h TOPOGRAPHY ARE ,//,3g11.4,///,
     &47h THE HIGH, LOW, AND MEDIAN POINTS OF THE BOTTOM,
     &13h SURFACE ARE ,//,3g11.4/)
      mx = cos(sincl * conv) * cos(sdecl * conv)
      my = cos(sincl * conv) * sin(sdecl * conv)
      mz = sin(sincl * conv)
      fx = cos(fincl * conv) * cos(fdecl * conv)
      fy = cos(fincl * conv) * sin(fdecl * conv)
      fz = sin(fincl * conv)
      write(6,107) mx, my, mz, fx, fy, fz
  107 format(/,45h THE COMPONENTS OF MAGNETIZATION AND REGIONAL
     &13h FIELD ARE...,//,5x,5hmx = ,f6.3,5x,5hmy = ,f6.3,5x,5hmz = 
     &,f6.3,/,5x,5hfx = ,f6.3,5x,5hfy = ,f6.3,5x,5hfz = ,f6.3/)
      aa1 = (mx * fx) - (mz * fz)
      a2 = (my * fy) - (mz * fz)
      a3 = (mx * fy) + (my * fx)
      b1 = (- (mx * fz)) - (mz * fx)
      b2 = (- (my * fz)) - (mz * fy)
c  do the n=0 term
        nc2=2*ncol
c   zero out the imaginary part
        do 30 j=1,ncol
        if(cmag.ne.0.0) a1(1,j,1)=cmag
   30   a1(2,j,1)=0.0
        open(iunita,access='direct',status='unknown',
     1  form='unformatted',file=drive//':\slave1.tmp',recl=nc2*4)
      if(cmag.eq.0.0) then
        open(11,file=ofilem,form='unformatted',status='old')
        read(11) id
      endif
c read magnetization grid and write to iunita
      do 450 jr=1,nrow
      if(cmag.eq.0.0) read(11)dum,(a1(1,i,1),i=1,ncol)
      call swrda(iunita,jr,a1,nc2)
  450 continue
      rnr=twopi/(float(nrow)*dy)
      rnc=twopi/(float(ncol)*dx)
c fourier transform the magnetization
  980 isign=-1
      open(iunitb,access='direct',status='unknown',
     1 form='unformatted',file=drive//':\slave2.tmp',recl=nc2*4)
      call sfftmg(iunita,iunitb,ncol,nrow,a1,nri,isign,work)
      ncp1=ncol+1
      nnhr=(nrow/2)+1
      nnhc=(ncol/2)+1
c   loop to operate on nnhr rows of the fourier transform a1
      open(iunitc,access='direct',status='unknown',
     1  form='unformatted',file=drive//':\slave3.tmp',recl=nc2*4)
      do 1010 jr=1,nnhr
      jj=jr-1
      x=float(jj)*rnr
      xsq=x*x
c   read a row of the fourier transform a1
      call srdda(iunita,jr,a1,nc2)
c   loop to operate on ncol elements in the row
      do 1090 i=1,ncol
      ii=i-1
      if(i.gt.nnhc) ii=-(ncp1-i)
      y=float(ii)*rnc
      xy=sqrt(xsq+y*y)
      argb=-xy*(bmed-zh)
      argt=-xy*(tmed-zh)
      if(argb.lt.-88.028) argb=-88.028
      if(argt.lt.-88.028) argt=-88.028
      argb=exp(argb)-exp(argt)
      a1(1,i,1)=a1(1,i,1)*argb
      a1(2,i,1)=a1(2,i,1)*argb
 1090 continue
c  save the filtered output to iunitc
      call swrda(iunitc,jr,a1,nc2)
 1010 continue
 1110 continue
      if(topscale.eq.0.0.and.botscale.eq.0.0) go to 5000
c get ready for summation
      write(6,110) err2
  110 format(/,25h PERFORMANCE OF SUMMATION,/,17h TRYING TO REACH ,f5.3,
     &10h OF STOTAL/)
      write(6,101)
  101 format(/,4x,1hN,4x,6hSTOTAL,5x,5hSLAST,4x,6hSRATIO/)
c begin summations
      n=0
   14 n=n+1
      open(iunitd,access='direct',status='unknown',
     1  form='unformatted',file=drive//':\slave4.tmp',recl=nc2*4)
c      nfact=fact(n)
      if(topscale.eq.0.0) go to 13
c      close(11)
      close(9)
      open(9,file=ofilet,form='unformatted',status='old')
      read(9) id
      if(cmag.eq.0.0) then
        rewind(11)
        read(11) id
      endif
      do 21 j=1,nrow
      call rowio(ztop,dlt,1,ncol,9,1)
      if(cmag.eq.0.0) then
        call rowio(s,dlt,1,ncol,11,1)
        do i=1,ncol
        a1(1,i,1)=s(i)*(ztop(i)-tmed)**n
        a1(2,i,1)=0.0
        enddo
      else
        do i=1,ncol
        a1(1,i,1)=cmag*(ztop(i)-tmed)**n
        a1(2,i,1)=0.0
        enddo
      endif
      call swrda(iunita,j,a1,nc2)
   21 continue
      call sfftmg(iunita,iunitb,ncol,nrow,a1,nri,isign,work)
c   loop to operate on nnhr rows of the fourier transform a1
      do 2010 jr=1,nnhr
      jj=jr-1
      x=float(jj)*rnr
      xsq=x*x
c   read a row of the fourier transform a1
      call srdda(iunita,jr,a1,nc2)
c   loop to operate on ncol elements in the row
      do 2090 i=1,ncol
      ii=i-1
      if(i.gt.nnhc) ii=-(ncp1-i)
      y=float(ii)*rnc
      xy=sqrt(xsq+y*y)
      argt=-xy*(tmed-zh)
      if(argt.lt.-88.028) argt=-88.028
      argt=(-xy)**n*exp(argt)/fact(n)
      a1(1,i,1)=-a1(1,i,1)*argt
      a1(2,i,1)=-a1(2,i,1)*argt
 2090 continue
c  save the filtered output to iunitd
      call swrda(iunitd,jr,a1,nc2)
 2010 continue
c      close(11)
   13 if(botscale.eq.0.0) go to 16
      close(9)
      open(9,file=ofileb,form='unformatted',status='old')
      read(9) id
      if(cmag.eq.0.0) then
        rewind(11)
        read(11) id
      endif
      do 31 j=1,nrow
      call rowio(zbot,dlt,1,ncol,9,1)
      if(cmag.eq.0.0) then
        call rowio(s,dlt,1,ncol,11,1)
        do i=1,ncol
        a1(1,i,1)=s(i)*(zbot(i)-bmed)**n
        a1(2,i,1)=0.0
        enddo
      else
        do i=1,ncol
        a1(1,i,1)=cmag*(zbot(i)-bmed)**n
        a1(2,i,1)=0.0
        enddo
      endif
      call swrda(iunita,j,a1,nc2)
   31 continue
      call sfftmg(iunita,iunitb,ncol,nrow,a1,nri,isign,work)
c   loop to operate on nnhr rows of the fourier transform a1
      do 3010 jr=1,nnhr
      jj=jr-1
      x=float(jj)*rnr
      xsq=x*x
c   read a row of the fourier transform a1
      call srdda(iunita,jr,a1,nc2)
c   loop to operate on ncol elements in the row
      do 3090 i=1,ncol
      ii=i-1
      if(i.gt.nnhc) ii=-(ncp1-i)
      y=float(ii)*rnc
      xy=sqrt(xsq+y*y)
      argb=-xy*(bmed-zh)
      if(argb.lt.-88.028) argb=-88.028
      argb=(-xy)**n*exp(argb)/fact(n)
      a1(1,i,1)=a1(1,i,1)*argb
      a1(2,i,1)=a1(2,i,1)*argb
 3090 continue
c  save the filtered output to iunitd
      if(topscale.ne.0.0) then
        call srdda(iunitd,jr,a1(1,1,2),nc2)
        do i=1,ncol
        a1(1,i,1)=a1(1,i,1)+a1(1,i,2)
        a1(2,i,1)=a1(2,i,1)+a1(2,i,2)
        enddo
      endif
      call swrda(iunitd,jr,a1,nc2)
 3010 continue
   16 continue
      slast=0.
      stotal=0.
      do j=1,nnhr
      call srdda(iunitc,j,a1,nc2)
      call srdda(iunitd,j,a1(1,1,2),nc2)
      do i=1,ncol
      slast=slast+sqrt(a1(1,i,2)**2+a1(2,i,2)**2)
      stotal=stotal+sqrt(a1(1,i,1)**2+a1(2,i,1)**2)
      a1(1,i,1)=a1(1,i,1)+a1(1,i,2)
      a1(2,i,1)=a1(2,i,1)+a1(2,i,2)
      enddo
      call swrda(iunitc,j,a1,nc2)
      enddo
      slast=slast/ntotal
      stotal=stotal/ntotal
      stest=err2*stotal
      sratio=slast/(stotal+.1e-20)
      write(6,104) n,stotal,slast,sratio
  104 format(i5,3g10.3)
      if(n.lt.nstop.and.slast.ge.stest) go to 14
 5000 continue
      do 4010 jr=1,nnhr
      jj=jr-1
      x=float(jj)*rnr
      xsq=x*x
      call srdda(iunitc,jr,a1,nc2)
      do 8 i=1,ncol
      ii=i-1
      if(i.gt.nnhc)ii=-(ncp1-i)
      y=float(ii)*rnc
      ysq=y*y
      xy=xsq+ysq
      if(xy.eq.0.0) go to 8
      phir=(aa1*xsq+a2*ysq+a3*x*y)/xy
      phii=(b1*x+b2*y)/sqrt(xy)
      temp=a1(1,i,1)
      a1(1,i,1)=temp*phir-a1(2,i,1)*phii
      a1(2,i,1)=temp*phii+a1(2,i,1)*phir
    8 continue
      call swrda(iunita,jr,a1,nc2)
      if(jr.eq.1) go to 4010
      jmr=nrow-jr+2
      if(jr.eq.jmr) go to 4110
      a1(2,1,1)=-a1(2,1,1)
      iwr=ncp1
      do 4100 i=2,nnhc
      iwr=iwr-1
      x=a1(2,i,1)
      y=a1(1,i,1)
      a1(2,i,1)=-a1(2,iwr,1)
      a1(1,i,1)=a1(1,iwr,1)
      a1(1,iwr,1)=y
      a1(2,iwr,1)=-x
 4100 continue
      call swrda(iunita,jmr,a1,nc2)
 4010 continue
 4110 continue
      isign=1
      call sfftmg(iunita,iunitb,ncol,nrow,a1,nri,isign,work)
      area=1.e5*twopi/float(ntotal)
      dum=0.0
      do 1210 jr=1,nr
      call srdda(iunita,jr,a1,nc2)
      do 1180 i=1,nc
 1180 a1(1,i,1)=a1(1,i,1)*area
 1210 write(10) dum,(a1(1,i,1),i=1,nc)
      return
      end
c***********************************************************************
c  Function fac finds the factorial of n
      function fac(n, ier)
      ier = 0
      if (n .gt. 1) goto 10
      if (n .ge. 0) goto 20
      ier = 1
   20 fac = 1.
      return 
   10 fac = n
      fac2 = fac
   30 fac2 = fac2 - 1.
      if (fac2 .eq. 1.) return 
      fac = fac * fac2
      goto 30
      end
c***********************************************************************
      subroutine srdbin(no,dat[huge],n)
c       subroutine 'srdbin' reads standard grids (binary).
c***********************************************************************
      dimension dat(n)
      read(no)dum,dat
      return
      end
c***********************************************************************
      subroutine swrbin(no,dat[huge],n)
c       subroutine 'swrbin' writes standard grids (binary).
c***********************************************************************
      dimension dat(n)
      data dum/0.0/
      write(no)dum,dat
      return
      end
c***********************************************************************
       subroutine srdda(no,ipos,dat[huge],n)
c       subroutine 'srrda' reads keyed sequential files.
c***********************************************************************
      dimension dat(n)
      read(no,rec=ipos)dat
      return
      end
c***********************************************************************
      subroutine swrda(no,ipos,dat[huge],n)
c       subroutine 'swrda' writes keyed sequential files.
c***********************************************************************
      dimension dat(n)
      write(no,rec=ipos)dat
      return
      end
c***********************************************************************
      subroutine get_dims(n1,n2,nx,ny,nri,id2,nxa)
      data lnri/16/,iw/6/
c
c  nri = blocking factor for rows (l)
c  nxa = number of rows added (n1-nx)
c  id2 = 2*max(n1,n2)
c
      nx=n1
      ny=n2
      l=lnri+1
      lnri21=lnri/2+1
c  set no. of rows for fft: need m=l*2**k, m.gt. or .eq. nx.
c  m=no. of rows, l=no. from 9-16, k=integer.
  190 l=l-1
      if(l.lt.lnri21) go to 260
      mr=n1/l+0.0000001
      k=1
      idiv=2
  200 iquot=mr/idiv+0.0000001
      if(iquot.lt.idiv) go to 210
      k=k+1
      mr=iquot
      go to 200
  210 k=k+1
      m=l*2**k
  220 mtest=l*2**(k-1)
      if(mtest.lt.n1) go to 230
      k=k-1
      m=mtest
      if(k.eq.0) go to 230
      go to 220
  230 lnxa=m-n1
      if(l.ne.lnri) go to 250
      nxa16=lnxa
  240 nri=l
      nxa=lnxa
      go to 190
  250 if(lnxa.ge.nxa) go to 190
      go to 240
  260 n1=n1+nxa
c  check to see if row block size of 16 will be more efficient
      n116=n1-nxa+nxa16
      ntest=0.9*n116
      if(ntest.gt.n1.or.n116.gt.2048) go to 270
      n1=n1-nxa+nxa16
      nxa=nxa16
      nri=16
  270 continue
      if(n1.gt.2048.or.n2.gt.2048) go to 320
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
      return
  320 write(iw,330)nx,ny,n1,n2
  330 format(' #no. of extended rows or columns exceeds 2048:'/
     1' input no. of rows and columns='2i4/
     2' no. of rows and columns required for fft =',2i4,/)
      stop
      end
      function ask_int(request)
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(1x,a,$)
      read(unit=5, fmt=*, err=1) ask_int
      return 
    1 write(unit=6, fmt=101) 
  101 format(/,20h **ERROR - try again)
      call bell
      goto 2
      end
c
      subroutine bell()
      character ding*1
      ding = char(7)
      write(unit=*, fmt=*) ding
      return 
      end
