c
c  --  PROGRAM BOUNDARY
c
c     A package of programs written by Rick Blakely, Bob Simpson, and
c  others (with inspiration from Lin Cordell, Tien Grauch, and Tom 
c  Hildenbrand) for locating edges of apparent magnetic or gravity 
c  sources.  See file BOUNDARY.HLP for further information, and put
c  proper location of BOUNDARY.HLP on open(10... statement.
c
c  *** MODIFICATIONS
c
c    Version 1.01 - (7/17/89) A bug was fixed that caused slight errors in 
c                 the location of maxspots.  Range of location error was 
c
c                       0 < ERROR < |max(dcol,drow)-sqrt(dcol**2+drow**2)| 
c
c                 with errors strongly skewed towards 0.  Magnitude of maxspots
c                 were not affected by this error.
c
c  Revisions re boundary.hlp file April 91
c
      write(6,100)
  100 format(/,' BOUNDARY - VERSION 1.01',//,
     &         ' See file BOUNDARY.HLP for information',//)
    1 write(6,101)
  101 format(/,' Operation [0-4]:',//,
     &         '     0.  Information',/,
     &         '     1.  Transform grid to pseudogravity',/,
     &         '     2.  Calculate horizontal gradient magnitude',/,
     &         '     3.  Locate maxima in grid',/,
     &         '     4.  Stop',/)
      ifunc=int_ask('  = ? ')
      if(ifunc.lt.0.or.ifunc.gt.4)call scold(*1)
      if(ifunc.eq.0)call info
      if(ifunc.eq.1)call pseudo
      if(ifunc.eq.2)call hgrad
      if(ifunc.eq.3)call maxima
      if(ifunc.eq.4)stop
      go to 1
      end
c-----------------------------------------------------------------------
      subroutine info
      character yn_ask*1,answer*1,line*75,flag*3,file1*87
c
c  --  !! Change following directory to location of BOUNDARY.HLP !!
c
c      open(10,file='\PF\HELP\BOUNDARY.HLP',status='old',
c     &      form='formatted',share='denywr')
      open(10,file='\pfhelp.tmp',form='formatted',status='old')
      read(10,102) line
  102 format(a)
      close(10)
      n=index(line,' ')
      file1=line(1:n-1)//'boundary.hlp'
      open(10,file=file1,form='formatted',status='old')
    2 read(10,100,end=3)flag,line
      if(flag.ne.'999')go to 2
    4 read(10,100,end=3)flag,line
  100 format(a3,a75)
      if(flag.eq.'999')then
         answer=yn_ask(' More info [y] ? ','y')
         if(answer.eq.'n')go to 3
         else
            write(6,101)flag,line
  101       format(1x,a3,a75)
            end if
      go to 4
    3 close(10)
      return
      end
c------------------------------------------------------------------------
      subroutine maxima
c
c  --  subroutine MAXIMA scans a standard grid along columns and rows.
c      Maxima are found, unprojected, and, if they meet certain criteria,
c      are sent to a binary xyz file.  Programmed by Rick Blakely from
c      an earlier algorithm by Bob Simpson.
c
c      common /label/ a(5000),b(5000),c(5000),x(50000),
c     &               y(50000),zmax(50000)
      dimension a(5000),b(5000),c(5000)
      real*8 sum,sum2,avg,signum,sigden,sig
      character id*56,pgm*8,yn_ask*1,answer*1,char_ask*1,downshift*1
c      data dval/0.1701412e39/,nmax/50000/
c      data nmax/50000/
c
c  --  Open the input file (standard grid) and read the header record
c
      write(6,102)
  102 format(/,' LOCATION OF MAXIMA WITHIN A GRID',//)
      call getfile(10,' Filename of input grid = ? ','in',
     &             'unformatted')
10    call rdheader(10,id,pgm,ncol,nrow,nz,col1,dcol,row1,drow,iproj,
     &              cmer,blat)
c
c  --  Fill arrays b and c with first two rows of grid
c
      call getrow(10,ncol,b)
      call getrow(10,ncol,c)
c
c  --  Initialize the projection routines
c
    8 write(6,105)
  105 format(
     & /,' Output is a binary xyz file.',/)
      answer=char_ask(
     &  ' x and y in cartesian or geodetic coordinates [g/c] ? ')
      answer=downshift(answer)
      if(answer.ne.'g'.and.answer.ne.'c')call scold(*8)
      if(answer.eq.'g')then
         if(iproj.eq.0)call getproj(iproj,cmer,blat)
         call projset(iproj,cmer,blat)
      else
         iproj=0
         end if
c
c  --  Determine level of significance
c
    6 write(6,100)
  100 format(
     & /, ' 1.  Significance levels = 1,2,3, and 4',
     & /, ' 2.       "         "    = 1',
     & /, ' 3.       "         "    = 2,3, and 4',
     & /, ' 4.       "         "    = 2',
     & /, ' 5.       "         "    = 3 and 4',
     & /, ' 6.       "         "    = 3',
     & /, ' 7.       "         "    = 4',/)
      leveltest=int_ask(' Pick one (1 through 7) ? ')
      if(leveltest.lt.1.or.leveltest.gt.7)call scold(*6)
c
c  --  Scan one row at a time
c
      close(12)
      open(12,status='scratch',form='unformatted')
      ntotal=0
      do 1 j=2,nrow-1
         do 2 i=1,ncol
            a(i)=b(i)
    2       b(i)=c(i)
         call getrow(10,ncol,c)
c
c  --  Call subroutine maximum at each grid intersection
c
         do 3 i=2,ncol-1
            if(b(i).ge.1.e38)go to 3
            call maximum(a(i-1),a(i),a(i+1),b(i-1),b(i),b(i+1),
     &                  c(i-1),c(i),c(i+1),dcol,drow,xoff,yoff,
     &                  ztest,level)
c
c  --  If significance level is satisfied, unproject the maxima and file away
c      for future use
c
            if(ipass(level,leveltest).eq.1)then
               ntotal=ntotal+1
c               if(ntotal.gt.nmax)then
c                  write(6,149)ntotal
c  149             format(/,' ***Number of maxima = ',i5,
c     &                     '...exceeds program limits')
c                  stop
c                   rewind(10)
c                   go to 10
c                  end if
               xmax=(i-1)*dcol+col1+xoff
               ymax=(j-1)*drow+row1+yoff
c               zmax(ntotal)=ztest
               zmax=ztest
               if(iproj.eq.0)then
c                  x(ntotal)=xmax
c                  y(ntotal)=ymax
                  x=xmax
                  y=ymax
               else
c                  call projinv(xmax,ymax,x(ntotal),y(ntotal))
                  call projinv(xmax,ymax,x,y)
               end if
               write(12) x,y,zmax
            end if
    3       continue
    1    continue
         rewind(12)
c
c  --  Calculate and report some statistics
c
      sum=0.d0
      sum2=0.d0
      zmaxmax=-1.e38
      zmaxmin= 1.e38
      do 4 i=1,ntotal
         read(12) x,y,zmax
c         zmaxmax=amax1(zmaxmax,zmax(i))
c         zmaxmin=amin1(zmaxmin,zmax(i))
c         sum=sum+zmax(i)
c    4    sum2=sum2+zmax(i)**2
         zmaxmax=amax1(zmaxmax,zmax)
         zmaxmin=amin1(zmaxmin,zmax)
         sum=sum+zmax
    4    sum2=sum2+zmax**2
      avg=sum/ntotal
      signum=float(ntotal)*sum2-sum**2
      sigden=float(ntotal)*float(ntotal-1)
      sig=sqrt(signum/sigden)
    7 write(6,101)ntotal,avg,sig,zmaxmax,zmaxmin
      rewind(12)
  101 format(
     & /,' The total set of maxima from this grid have the following',
     & /,' statistics:',//,
     &         ' Number = ',i7,/,
     &         ' Average = ',g11.4,/,
     &         ' Standard deviation = ',g11.4,/,
     &         ' Maximum = ',g11.4,/,
     &         ' Minimum = ',g11.4,/)
c
c  -- Determine the upper and lower thresholds
c
      write(6,104)
  104 format(/,' Specify range of maxima to be output...',/)
      thlow=float_ask(' Lower limit = ? ')
      thup=float_ask(' Upper limit = ? ')
      if(thup.le.thlow)call scold(*7)
      call getfile(11,' Filename for output list = ? ','out',
     &             'unformatted')
c
c  --  Send all maxima that lie between the thresholds to the ouput file
c
      nout=0
      do 5 i=1,ntotal
         read(12) x,y,zmax
c         if(zmax(i).ge.thlow.and.zmax(i).
         if(zmax.ge.thlow.and.zmax.
     &       lt.thup)then
            nout=nout+1
c            write(11)x(i),y(i),zmax(i)
            write(11) x,y,zmax
            end if
    5    continue
      write(6,103)nout
  103 format(/,i10,' records written',/)
      close(11)
c
c  --  Go back for new thresholds
c
      answer=yn_ask(
     &   ' Another output file with a new range [n] ? ','n')
      if(answer.eq.'y')go to 7
c
c  --  Terminate
c
      close(10)
      close(12)
      return
      end
c---------------------------------------------------------------------------
      subroutine maximum(a1,a2,a3,b1,b2,b3,c1,c2,c3,dcol,drow,xoff,
     &                   yoff,zmax,level)
c
c     Subroutine maximum determines if point b2 is a maximum with respect
c  to its neighbors, a1,a2,a3,b1,b3,c1,c2, and c3.
c
      xoff=0.
      yoff=0.
      zmax=-1.e38
      level=0
      diag=sqrt(dcol**2+drow**2)
      sinthe=drow/diag
      costhe=dcol/diag
      if(b2.gt.b1.and.b2.gt.b3)then
         level=level+1
         call interp(b1,b2,b3,dcol,dtest,ztest,ierr)
         if(ztest.gt.zmax)then
            zmax=ztest
            xoff=dtest
            yoff=0.
            end if
         end if
      if(b2.gt.a2.and.b2.gt.c2)then
         level=level+1
         call interp(a2,b2,c2,drow,dtest,ztest,ierr)
         if(ztest.gt.zmax)then
            zmax=ztest
            xoff=0.
            yoff=dtest
            end if
         end if
      if(b2.gt.a1.and.b2.gt.c3)then
         level=level+1
         call interp(a1,b2,c3,diag,dtest,ztest,ierr)
         if(ztest.gt.zmax)then
            zmax=ztest
            xoff=dtest*costhe
            yoff=dtest*sinthe
            end if
         end if
      if(b2.gt.a3.and.b2.gt.c1)then
         level=level+1
         call interp(a3,b2,c1,diag,dtest,ztest,ierr)
         if(ztest.gt.zmax)then
            zmax=ztest
            xoff=-dtest*costhe
            yoff=dtest*sinthe
            end if
         end if
      return
      end
c---------------------------------------------------------------------------
      subroutine interp(y1,y2,y3,d,xmax,ymax,ierr)
c
c     Subroutine interp fits a parabola through ordinants y1,y2 and y3, assumed
c  to be evenly spaced, and reports the x and y coordinate of the parabola's
c  maximum.  Note...y1 must be .lt. y2 and y3 must be .lt. y2.
c
      ierr=0
      a=0.5*(y1-2.*y2+y3)/(d*d)
      b=0.5*(y3-y1)/d
      c=y2
      if(a.eq.0.)then
         ierr=1
         return
         end if
      xmax=-0.5*b/a
      ymax=a*xmax**2+b*xmax+c
      return
      end
c-------------------------------------------------------------------------
      subroutine getfile(unit,question,inout,form)
c      Opens files...filenames not recognized by VAX will generate
c      error message and request a second try.
c      Sample call:
c           call getfile(21,' *Give input file:','in','formatted')
      integer unit
      character filename*50, status*8
      character*(*) question, inout, form
c
      print *, ' '
 10   write(*,101) question
 101  format(a,' ',$)
      read '(a50)', filename
c
      if(inout.eq.'in') then
        status='old'
      else
        status='unknown'
      endif
      if(form.eq.'formatted'.and.inout.eq.'out') then
        open(unit,file=filename,form=form,status=status,
     &   err=900)
        return
        end if
      if(inout.eq.'in')then
        open(unit,file=filename,form=form,status=status,
     &    share='denywr',err=900)
        else
        open(unit,file=filename,form=form,status=status,err=900)
        endif
        return
c
 900  print *,' ERROR IN OPENING ',filename
      print *,'   try again......'
      close(unit)
      go to 10
c
      end
c----------------------------------------------------------------------
      subroutine rdheader(unit,id,pgm,ncol,nrow,nz,
     & xo,dx,yo,dy,iproj,cm,bl)
      character id*56, pgm*8
      integer unit
      read(unit,err=20) id,pgm,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
c      print '(1h0,a,i2,2(a,f10.4))',
c     &    '  Proj=',iproj,'  cm=',cm,'  bl=',bl
      return
 20   rewind unit
      read(unit) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
c      print '(1h0,a)', '  Grid has no projection specs in header...'
      iproj=0
      cm=0.0
      bl=0.0
      return
      end
c-------------------------------------------------------------------------
      subroutine getrow(iounit,ncol,row)
      real row(ncol)
      read(iounit) dummy,row
      return
      end
c------------------------------------------------------------------------
      function yn_ask(request,default)
      character*1 yn_ask,default
      character*(*) request
    2 write(6,100)request
  100 format(a,' ',$)
      read(5,102,err=1)yn_ask
      if(yn_ask.eq.' ') yn_ask=default
  102 format(a1)
      if(yn_ask.eq.'Y')yn_ask='y'
      if(yn_ask.eq.'N')yn_ask='n'
      if(yn_ask.ne.'y'.and.yn_ask.ne.'n')call scold(*2)
      return
    1 call scold(*2)
      end
c-------------------------------------------------------------------------
      subroutine bell
      character*1 ding
      ding=char(007)
      print *, ding
      return
      end
c-----------------------------------------------------------------------
      subroutine scold(*)
      call bell
      write(6,100)
  100 format(/,' *** Unacceptable response...Try again')
      return 1
      end
c----------------------------------------------------------------------
      function ipass(level,leveltest)
      ipass=0
      if(level.eq.0)return
      go to (1,2,3,4,5,6,7)leveltest
    1 if(level.ge.1)ipass=1
      return
    2 if(level.eq.1)ipass=1
      return
    3 if(level.ge.2)ipass=1
      return
    4 if(level.eq.2)ipass=1
      return
    5 if(level.ge.3)ipass=1
      return
    6 if(level.eq.3)ipass=1
      return
    7 if(level.ge.4)ipass=1
      return
      end
c----------------------------------------------------------------------------
      subroutine pseudo
c
c     Calculates the pseudogravity transformation.  Modified by Rick Blakely
c  from Bob Simpson's GFILTER.
c
      real finc,fdec,minc,mdec,az
      character id*56,pgm*8
      common/gridid/ id,pgm
      common/gridspecs/ ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      common/origspecs/ ncol0,nrow0
c      complex grid(500000), fftgrid(500000)
      complex grid(10000), fftgrid(10000)
      common/grid/ grid
      common/fftgrid/ fftgrid
c      real kx(5000),ky(5000),kxsq(5000),kysq(5000)
      real kx(1000),ky(1000),kxsq(1000),kysq(1000)
      common/wavek/ kx,ky,kxsq,kysq
      data inunit/21/
c
      write(6,100)
  100 format(
     & /,' PSEUDOGRAVITY TRANSFORMATION',
     &//,' NOTE:  Automatic augmented dimensions are about 15% bigger.',
     & /,'        Shallow inclinations (<10?) may cause errors.',/)
c
c  --  Get grid to be filtered and read header
c
      call getfile(inunit,' Filename of input grid = ? ','in',
     &             'unformatted')
      call read_sfheader(inunit)
      if(ncol*nrow.gt.15000) then
        print *, ' GRID IS TOO BIG...ncol*nrow =',ncol*nrow
        print *, ' NCOL*NROW MUST BE .LE. 15000...'
        return
      endif
      ncol0=ncol
      nrow0=nrow
c
c  --  Read in grid, augment, and take fft...
      call get_newfft(inunit)
c
c  --  Calculate and store wavenumbers for fftfiltering...
c
 90   call wavenums(ncol,nrow,dx,dy)
c
c  --  Do the transformation
c
      finc=float_ask(' Field inclination (degrees) = ? ')
      fdec=float_ask(' Field declination = ? ')
      minc=float_ask(' Magnetization inclination = ? ')
      mdec=float_ask(' Magnetization declination = ? ')
      az=  float_ask
     &     (' Azimuth of grid relative to north (usually 0) = ? ')
         call get_units(ufactor)
         call psgrav(finc,fdec,minc,mdec,az,ufactor)
c
c  --  Inverse transform and output
c
         call output
c
c  --  Reset 'grid' to 'fftgrid' for further operations...
c
         do 185 ind=1,ncol*nrow
 185     grid(ind)=fftgrid(ind)
 900  return
      end
c-----------------------------------------------------------------------------
      subroutine get_newfft(inunit)
c
c  --  reads in grid, augments, and calculates the fftgrid
c
c      dimension row(5000)
c      complex crow(5000)
      dimension row(1000)
c      complex crow(1000)
c      dimension nn(2), work(5000)
      dimension nn(2), work(1000)
      character augans*1, yn_ask*1
      character id*56,pgm*8
      common/gridid/ id,pgm
      common/gridspecs/ ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      common/origspecs/ ncol0,nrow0
c      complex grid(500000), fftgrid(500000)
      complex grid(10000), fftgrid(10000)
      common/grid/ grid
      common/fftgrid/ fftgrid
      integer inunit
c      data outunit/22/
      ij(i,j,ncol)=ncol*(j-1)+i
c**** read in body of grid
      do 30 j=1,nrow
      call row_read(inunit,ncol,row)
      do 30 i=1,ncol
 30   grid(ij(i,j,ncol))=cmplx(row(i),0.)
      close(inunit)
c**** check grid for dvals...stop if found.
      if(ndval(2*ncol*nrow,grid).ne.0) then
       print *,' ERROR...GRID HAS DVALS...FILL THEM IN!'
       stop
      endif
c**** augment grid to speed fft
   40 augans=yn_ask(' Do you want automatic augmentation [y] ? ','y')
      if(augans.eq.'y') then
        call piknewdims(ncol0,nrow0,ncol,nrow,15.)
      else if(augans.eq.'n') then
        call fftdims(ncol0,nrow0)
      write(6,101)
  101 format(/,' Give aumented ncol and nrow ? '$)
        read *, ncol,nrow
      else
        goto 40
      endif
      write(6,1001) ncol,nrow
 1001 format(/,10x,'ncolp=',i8,'   nrowp=',i8,/)
      if(ncol*nrow.gt.10000) then
        print *, ' AUGMENTED GRID IS TOO BIG...ncol*nrow =',ncol*nrow
        print *, ' NCOL*NROW MUST BE .LE. 10000...'
        return
      endif
      call cmplxaug(ncol0,nrow0,ncol,nrow,+1)
c**** take fft of augmented grid.
      nn(1)=ncol
      nn(2)=nrow
      call fourt(grid,nn,2,-1,0,work)
c**** Transfer to fftgrid to save for future operations...
      do 50 ind=1,ncol*nrow
 50   fftgrid(ind)=grid(ind)
 90   return
      end
c-------------------------------------------------------------------------
      subroutine output
c
c  --  Takes inverse fft and outputs filtered standard file.
c
c      dimension row(5000)
c      dimension nn(2), work(5000)
      dimension row(1000)
      dimension nn(2), work(1000)
c      character ans*4
      character id*56,pgm*8
      common/gridid/ id,pgm
      common/gridspecs/ ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      common/origspecs/ ncol0,nrow0
c      complex grid(500000)
      complex grid(10000)
      common/grid/ grid
      integer outunit
      data outunit/22/
      ij(i,j,ncol)=ncol*(j-1)+i
      nn(1)=ncol
      nn(2)=nrow
      call fourt(grid,nn,2,+1,1,work)
      do 10 ind=1,ncol*nrow
 10   grid(ind)=grid(ind)/(ncol*nrow)
c**** deaugment
      call cmplxaug(ncol0,nrow0,ncol,nrow,-1)
c**** output grid
      call getfile(outunit,' Filename of output grid = ? ',
     &             'out','unformatted')
      write(6,101)
  101 format(/,' New ID (56 characters or less) = ? '$)
      read '(a56)', id
      pgm='filter  '
      call wrheader(outunit,id,pgm,ncol0,nrow0,nz,xo,dx,yo,dy,
     &   iproj,cm,bl)
      do 50 j=1,nrow0
      do 40 i=1,ncol0
 40   row(i)=real(grid(ij(i,j,ncol0)))
      call row_write(outunit,ncol0,row)
 50   continue
      close(outunit)
      print *,' '
      return
      end
c-----------------------------------------------------------------------
      complex function phase(x,y,z,kx,ky,k)
      real kx,ky,k
      phase=cmplx(+z,(x*kx+y*ky)/k)
      return
      end
c------------------------------------------------------------------------
      subroutine psgrav(finc,fdec,minc,mdec,az,ufactor)
c
c     Calculates pseudograv. f=field direction, m=mag direction
c  density contrast= 0.1 g/cc; mag= 0.01 emu.  y=north,x=east,z=dn,
c  or... az= angle of y from north
c
      real jay,kmtocm
      real finc,fdec,minc,mdec,fx,fy,fz,mx,my,mz
      real k
c      real kx(5000),ky(5000),kxsq(5000),kysq(5000)
      real kx(1000),ky(1000),kxsq(1000),kysq(1000)
      common/wavek/ kx,ky,kxsq,kysq
      complex phase,filt
c      complex grid(500000)
      complex grid(10000)
      common/grid/ grid
      common/gridspecs/ ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      data rho/0.1/,jay/.001/,gbig/6.670E-8/
      data kmtocm/1.0E5/,gamtogauss/1.0E-5/,galtomgal/1.0E+3/
      ij(i,j,ncol)=ncol*(j-1)+i
      call cmpnts(finc,fdec-az,fx,fy,fz)
      call cmpnts(minc,mdec-az,mx,my,mz)
      const=gbig*rho*kmtocm*gamtogauss*galtomgal/jay
      do 100 j=1,nrow
      do 100 i=1,ncol
      k=sqrt(kxsq(i)+kysq(j))
      filt=cmplx(0.,0.)
      if(i.eq.1.and.j.eq.1) go to 100
c      this avoids the case k=0.
      filt=1./(phase(fx,fy,fz,kx(i),ky(j),k)*
     1            phase(mx,my,mz,kx(i),ky(j),k))
      filt=const*filt/(k*ufactor)
c      [units specified]/ufactor = km
  100 grid(ij(i,j,ncol))=filt*grid(ij(i,j,ncol))
      return
      end
c-------------------------------------------------------------------------
      subroutine cmpnts(inc,dec,x,y,z)
      real inc,dec
      data rad/1.745329E-2/
      x=cos(inc*rad)*sin(dec*rad)
      y=cos(inc*rad)*cos(dec*rad)
      z=sin(inc*rad)
      return
      end
c----------------------------------------------------------------------
      subroutine get_units(ufactor)
      character units*8
    1 write(6,100)
  100 format(/,' Units of grid dx and dy (km,m,kft,ft,mi) = ? '$)
      read '(a8)', units
      if(units.eq.'km'.or.units.eq.'KM')then
         ufactor=1.0
         return
         end if
      if(units.eq.'m'.or.units.eq.'M')then
         ufactor=1000.
         return
         end if
      if(units.eq.'kft'.or.units.eq.'KFT')then
         ufactor=3.2802
         return
         end if
      if(units.eq.'ft'.or.units.eq.'FT')then
         ufactor=3.2802e3
         return
         end if
      if(units.eq.'mi'.or.units.eq.'MI')then
         ufactor=0.6214
         return
         end if
      call scold(*1)
      return
      end
c----------------------------------------------------------------------------
      subroutine wavenums(ncol,nrow,dx,dy)
c      calculate wavenumbers to be used in filters
c      real kx(5000),ky(5000),kxsq(5000),kysq(5000)
      real kx(1000),ky(1000),kxsq(1000),kysq(1000)
      common/wavek/ kx,ky,kxsq,kysq
      data pi/3.1415927/
      xscale=2.*pi/(ncol*dx)
      yscale=2.*pi/(nrow*dy)
      do 100 j=1,nrow
      jk=j-1
      if(jk.gt.nrow/2) jk=jk-nrow
      ky(j)=yscale*jk
 100  kysq(j)=ky(j)**2
      do 200 i=1,ncol
      ik=i-1
      if(ik.gt.ncol/2) ik=ik-ncol
      kx(i)=xscale*ik
 200  kxsq(i)=kx(i)**2
      return
      end
c--------------------------------------------------------------------------
      subroutine cmplxaug(ncol,nrow,ncolp,nrowp,idir)
c
c     Modified from a program written by R. Blakely.  Augments a 2d complex
c  grid stored in a linear array.  ncol,nrow = original dimensions;
c  ncolp,nrowp = augmented dimensions.  idir=+1 to augment, -1 to go back
c  to original size.  Grid is stored in common block grid row by row,
c  starting with bottom row.
c
      complex grid,first,last,diff
c      common/grid/ grid(500000)
      common/grid/ grid(10000)
c
      if(ncolp.lt.ncol.or.nrowp.lt.nrow) go to 900
      if(ncolp.eq.ncol.and.nrowp.eq.nrow) go to 910
c
c     Deaugmenting steps:
      if(idir.ne.-1) go to 200
      do 110 j=1,nrow
      do 110 i=1,ncol
 110  grid(ncol*(j-1)+i)=grid(ncolp*(j-1)+i)
      return
c
c     Augmenting steps:
 200  if(idir.ne.+1) go to 920
 210  do 220 j=nrow,1,-1
      do 220 i=ncol,1,-1
 220  grid(ncolp*(j-1)+i)=grid(ncol*(j-1)+i)
c     add to rows:
      if(ncolp.eq.ncol) go to 235
        do 230 j=1,nrow
        first=grid(ncolp*(j-1)+1)
        last=grid(ncolp*(j-1)+ncol)
        diff=(first-last)/float(ncolp-ncol+1)
        do 230 i=ncol+1,ncolp
 230    grid(ncolp*(j-1)+i)=last+diff*float(i-ncol)
c     add to cols:
 235  if(nrowp.eq.nrow) go to 250
        do 240 i=1,ncolp
        first=grid(i)
        last=grid(ncolp*(nrow-1)+i)
        diff=(first-last)/float(nrowp-nrow+1)
        do 240 j=nrow+1,nrowp
 240    grid(ncolp*(j-1)+i)=last+diff*float(j-nrow)
 250  return
c
 900  print *,' ***AUGMENTED GRID IS SMALLER THAN ORIGINAL!!!'
      return
 910  print *,' ***AUGMENTED GRID IS SAME SIZE AS ORIGINAL!!!'
      return
 920  print *,' idir in subroutine aug must be +1 or -1!!!'
      return
      end
c-------------------------------------------------------------------------
      subroutine piknewdims(ncol0,nrow0,ncol,nrow,pct)
c
c     Returns new dimensions approx 'pct' percent larger than originals.
c  New dimensions should be fast for fft.
c  first 2 cols from table III of Singleton's paper on fft....
c  (Singleton,R.C., 1969, An algotithm for computing the mixed radix fast
c  Fourier transform:  IEEE Trans., Audio and Electro-acoustics: v. AU-17,
c  p. 93-103.
c
      dimension nlist(104)
      data nlist/2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,
     1   32,36,40,45,48,50,54,60,64,72,75,80,81,90,96,100,108,120,
     2   125,128,135,144,150,160,162,180,192,200,216,225,240,
     3   243,250,256,270,288,300,320,324,360,375,384,400,405,432,450,
     4   480,486,500,512,540,576,600,625,640,648,675,720,729,750,
     5   768,800,810,864,900,960,972,1000,1024,1080,1125,1152,1200,
     6   1215,1250,1280,1296,1350,1440,1458,1500,1536,1600,1620,1728,
     7   1800,1875/
c
      ncol=(1.0+.01*pct)*ncol0
      nrow=(1.0+.01*pct)*nrow0
      do 50 i=1,104
      if(nlist(i).ge.ncol) ncol=nlist(i)
      if(nlist(i).eq.ncol) go to 55
 50   continue
 55   do 60 i=1,104
      if(nlist(i).ge.nrow) nrow=nlist(i)
      if(nlist(i).eq.nrow) go to 90
 60   continue
 90   if(ncol0.eq.1) ncol=1
      if(nrow0.eq.1) nrow=1
      return
      end
c------------------------------------------------------------------------------
      subroutine fftdims(ncol,nrow)
c     Prints out numbers less than 1024 with no prime factor greater than 5.
c  The FFT algorithm works fastest for these numbers.  Subroutine prints
c  out good numbers between min of (ncol,nrow) and max of (ncol+100,nrow+100).
c  first 2 cols from table III of Singleton's paper on fft....
c  (Singleton,R.C., 1969, An algotithm for computing the mixed radix fast
c  Fourier transform:  IEEE Trans., Audio and Electro-acoustics: v. AU-17,
c  p. 93-103.
c
      dimension nlist(104),ngood(100)
      data nlist/2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,
     1   32,36,40,45,48,50,54,60,64,72,75,80,81,90,96,100,108,120,
     2   125,128,135,144,150,160,162,180,192,200,216,225,240,
     3   243,250,256,270,288,300,320,324,360,375,384,400,405,432,450,
     4   480,486,500,512,540,576,600,625,640,648,675,720,729,750,
     5   768,800,810,864,900,960,972,1000,1024,1080,1125,1152,1200,
     6   1215,1250,1280,1296,1350,1440,1458,1500,1536,1600,1620,1728,
     7   1800,1875/
c
      min=min0(ncol,nrow)
      max=max0(ncol+100,nrow+100)
      if(min.eq.1) min=max0(ncol,nrow)
      n=0
      do 10 i=1,104
      if(nlist(i).lt.min) go to 10
      if(nlist(i).gt.max) go to 20
      n=n+1
      ngood(n)=nlist(i)
 10   continue
c
 20   print *,'    *Good fft dimensions:'
      write(6,101) (ngood(i),i=1,n)
 101  format('      ',14i5)
c
      return
      end
c-----------------------------------------------------------------------------
      subroutine read_sfheader(iunit)
      character id*56,pgm*8
      common/gridid/ id,pgm
      common/gridspecs/ ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      call rdheader(iunit,id,pgm,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl)
      call prheader(6,id,pgm,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl)
      return
      end
c------------------------------------------------------------------------------
      subroutine print_sfspecs
      character id*56,pgm*8
      common/gridid/id,pgm
      common/gridspecs/ ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      print *,'      id=',id
      print *,'      pgm=',pgm
      print *,'      nz=  ',nz
      print *,'      ncol=',ncol,'   xo=',xo,'   dx=',dx
      print *,'      nrow=',nrow,'   yo=',yo,'   dy=',dy
      print *,' '
      print *,'      iprj=',iproj,'   cm=',cm,'   bl=',bl
      return
      end
c---------------------------------------------------------------------------
      subroutine row_read(iunit,ncol,zrow)
      dimension zrow(ncol)
      read(iunit) dum,zrow
      return
      end
c----------------------------------------------------------------------------
      subroutine row_write(iunit,ncol,zrow)
      dimension zrow(ncol)
      dum=0.
      write(iunit) dum,zrow
      return
      end
c----------------------------------------------------------------------------
      subroutine crow_read(iunit,ncol,zrow)
      complex zrow(ncol)
      read(iunit) dum,zrow
      return
      end
c----------------------------------------------------------------------------
      subroutine crow_write(iunit,ncol,zrow)
      complex zrow(ncol)
      dum=0.
      write(iunit) dum,zrow
      return
      end
c--------------------------------------------------------------------------
      subroutine ask(ques)
      character*(*) ques
      character*40 string
      string=ques
c
c  --  This pads ques to left with blanks (or truncates).
c
      write(6,101) string
 101  format(5x,a40,$)
      return
      end        
c--------------------------------------------------------------------------
      subroutine hgrad
c
c     Calculates horizontal gradients...Modified slightly by Rick Blakely
c  from a program by Bob Simpson.
c
c      character*50 infil,outfil
      character id*56,pgm1*8,pgm2*8
c      dimension a(5000),b(5000),c(5000),grad(5000)
      dimension a(1000),b(1000),c(1000),grad(1000)
      data dcval/0.1701412e39/,dcvaltest/1.e37/,pgm2/'gradient'/
      rad=3.1415927/180.
c
c      print *,' *Mode: 1=magnitude of gradient'
c      print *,'        2=direction of gradient (0 to 360)'
c      print *,'        3=trend of gradient (0 to 180)'
c      print *,' *Give mode'
c      read *,mode
c
      write(6,101)
  101 format(/,' HORIZONTAL GRADIENT CALCULATION ',//)
      call getfile(10,' Filename of input grid = ? ','in',
     &             'unformatted')
      call getfile(11,' Filename of output grid = ? ','out',
     &             'unformatted')
c
      call rdheader(10,id,pgm1,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl)
      print *,id
      print *,pgm1
      print *,' ncol=',ncol,'     nrow=',nrow
      print *,' xo=  ',xo,'     yo=  ',yo
      print *,' dx=  ',dx,'     dy=  ',dy
      write(6,102)
  102 format(/,' New ID (56 characters or less) = ? '$)
      read '(a56)', id
      call wrheader(11,id,pgm2,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl)
c
      call getrow(10,ncol,b)
      call getrow(10,ncol,c)
c
      do 10 i=1,ncol
 10   grad(i)=dcval
      call putrow(11,ncol,grad)
c
      do 200 j=2,nrow-1
      do 120 i=1,ncol
      a(i)=b(i)
 120  b(i)=c(i)
      call getrow(10,ncol,c)
      do 40 i=1,ncol
 40   grad(i)=dcval
      do 100 i=2,ncol-1
      if(b(i-1).gt.dcvaltest.or.b(i+1).gt.dcvaltest) go to 100
      if(a(i).gt.dcvaltest.or.c(i).gt.dcvaltest) go to 100
      dzdx=(b(i+1)-b(i-1))/(2*dx)
      dzdy=(c(i)-a(i))/(2*dy)
      grad(i)=sqrt(dzdx**2+dzdy**2)
c      if(mode.eq.1) grad(i)=sqrt(dzdx**2+dzdy**2)
c      if(mode.eq.2) grad(i)=atan2(dzdy,dzdx)/rad
c      if(mode.eq.3) grad(i)=amod((atan2(dzdy,dzdx)+360.)
c     1  /rad,180.)
 100  continue
c
      call putrow(11,ncol,grad)
 200  continue
c
      do 510 i=1,ncol
 510  grad(i)=dcval
      call putrow(11,ncol,grad)
c
      close(10)
      close(11)
      return
      end
c------------------------------------------------------------------------------
      function char_ask(request)
      character*1 char_ask
      character*(*) request
    2 write(6,100)request
  100 format(a,' ',$)
      read(5,102,err=1)char_ask
  102 format(a1)
      return
    1 write(6,101)
  101 format(/,' **ERROR - try again')
      call bell
      go to 2
      end
c------------------------------------------------------------------------
      function downshift(lchar)
      character*1 downshift,lchar
      downshift=lchar
      i=ichar(lchar)
      if(i.ge.65.and.i.le.90)then
         downshift=char(i+32)
         end if
      return
      end
c------------------------------------------------------------------------------
      function float_ask(request)
      character*(*) request
    2 write(6,100)request
  100 format(a,' ',$)
      read(5,*,err=1)float_ask
      return
    1 write(6,101)
  101 format(/,' **ERROR - try again')
      call bell
      go to 2
      end
