c  prep.for
c
c  To extend a grid, tilt, 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 FFTFIL likes.  The program checks that
c  nrow and ncol will not be increased by FFTFIL.
c
c  Lin Cordell Jan 1990.
c  Revised June 1991.
c
        character*56 ifile,ofile
c      dimension g(2048),g1(2048),g2(2048),id(14),pgm(2),ncnr(2)
      dimension g(2048),g1(2048),g2(2048),ncnr(2)
      character id*56, prjnam*80
      dimension cmbl(2)
      data lnri/16/
      data ddval/1.e29/
	pi=3.141597
        ifile=' '
        call askin
        call pfinit('prep')
400     call askc('Enter input grid name (NO DVALS !)',ifile,ierr)
        if(ierr.lt.0) stop
        call gopen(10,ifile,'old','read',ierr)
        if(ierr.ne.0) then
          print*,'Error - try again'
          go to 400
        endif
        call gheader('r',10,id,nc,nr,xo,dx,yo,dy,ierr)
        print*,id
        if(ierr.ne.0) stop 'error reading header'
        call fftdims(nc,nr,ncnr(1),ncnr(2))
401     call aski4a('Enter (new) ncol and nrow',ncnr,2,ierr)
        if(ierr.eq.-2) then
          close(10)
          go to 400
        endif
        ncol=ncnr(1)
        nrow=ncnr(2)
      if(nc.gt.ncol.or.nr.gt.nrow)then
        print*,'Too small.'
        stop
      endif
c  Check to make sure FFTFIL will use these unchanged.
c  Added 19 June 91.
      kr=11
      iw=6
      ir=5
      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
      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+2*nadd.
c  m=no. of rows, l=no. from 9-16, k=interger.
      n1=n1+2*nadd
  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 n2=n2+2*nadd
      if(n1.gt.2048.or.n2.gt.2048) then
        print*,'Too big.'
        stop
      endif
      nxa=nxa+2*nadd
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
      itest=n1-nx
      if(itest.ne.0) then
        print*,'Select another nrow; FFTFIL would change this one.'
        stop
      endif
c  Tilt grid.
      write(iw,275)nri
  275 format(' blocking rows =',i3)
        call grow('r',10,1,g,nc,ierr)
        do 9 i=1,nc
9       if(g(i).gt.ddval) stop 'dvals found'
	a=g(1)
	b=g(nc)
	do 10 j2=2,nr
        call grow('r',10,j2,g,nc,ierr)
        do 10 i=1,nc
10      if(g(i).gt.ddval) stop 'dvals found'
        c=g(1)
	d=g(nc)
	alpha=(c-b)/float(2*(nc-1))
	beta=(b-c)/float(2*(nr-1))
	print*, ' alpha = ',alpha,', beta = ',beta
c        rewind(10)
        call gheader('r',10,id,nc,nr,xo,dx,yo,dy,ierr)
        call mpsnam('r',idum,prjnam,ncharp,ierr)
        if(prjnam.eq.'not specified') then
          cmbl(1)=0.
          cmbl(2)=0.
        else
          call mpsref('r',idum,cmbl,nref,ierr)
        endif
      open(14,file='prep.rec',status='new',err=280,form='formatted')
      go to 281
  280 ians=1
      call aski4l('ok to overwrite prep.rec?',ians,ierr)
      if(ians.eq.0) stop 'user exit'
      open(14,file='prep.rec',status='unknown',form='formatted')
  281   write(14,300)ofile,alpha,beta,nc,nr,prjnam,cmbl
300     format(a50/2e16.8,2i5/a50/2e16.8)
      close(14)
c  Output extended file.
        ofile=' '
        ofile=ifile(1:index(ifile,'.'))//'prp'
301     call askc('Enter output grid name',ofile,ierr)
        if(ierr.eq.-2) then
          close(10)
          go to 401
        endif
        call gopen(11,ofile,'new','write',ierr)
        if(ierr.ne.0) then
          print*,'Error -try again'
          go to 301
        endif
        call gheader('w',11,id,ncol,nrow,xo,dx,yo,dy,ierr)
c  First row (saved).
	j2=1
        call grow('r',10,j2,g,nc,ierr)
      do 211 j1=1,nc
	x=float(j1-1)
	g(j1)=g(j1)+alpha*x
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 grow('w',11,j2,g,ncol,ierr)
c  Rows down to nr.
	do 204 j2=2,nr
	y=float(j2-1)
        call grow('r',10,j2,g,nc,ierr)
	do 202 j1=1,nc
	x=float(j1-1)
202	g(j1)=g(j1)+alpha*x+beta*y
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 grow('w',11,j2,g,ncol,ierr)
204	continue
c  Last row.
	do 205 j1=1,ncol
205	g2(j1)=g(j1)
c  Rows nr to nrow.
c  Cosine taper.
	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 grow('w',11,j2,g,ncol,ierr)
207	continue
	go to 999
900	print*, 'Cant handle it.'
	close(10)
999	close(11)
	stop
	end
c **************************************************************
      subroutine fftdims(ncol, nrow, ncout, nrout)
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)
      ncout=-1
      nrout=-1
      n = 0
      do 10 i = 1, 54
      if (nlist(i) .lt. ncol) goto 10
      if (nlist(i) .gt. ncol+ncol) 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)
      do 11 i=1,n
      if(ngood(n-i+1).ge.(ncol+20)) ncout=ngood(n-i+1)
   11 continue
      n = 0
      do 30 i = 1, 54
      if (nlist(i) .lt. nrow) goto 30
      if (nlist(i) .gt. nrow+nrow) 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)
      do 31 i=1,n
      if(ngood(n-i+1).ge.(nrow+20)) nrout=ngood(n-i+1)
   31 continue
      write(unit=6, fmt=901) ncout, nrout
      return 
      end

