c  prep3.for
c
c  Extends a grid 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.  Undo the extension
c  with program de_prep3.
c
c  prep:  Lin Cordell Jan 1990.
c         Revised June 1991.
c  prep3: Jeff Phillips February 1996
c
        character*56 ifile,ofile,id
        character prjnam*80
c      dimension id(14),pgm(2),ncnr(2),cmbl(2)
      dimension ncnr(2),cmbl(2)
      data lnri/16/
        pi=3.1415927
1       format(a)
        ifile=' '
        call askin
        call pfinit('prep3')
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) stop 'Too small'
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  open and check input file
      write(iw,275)nri
  275 format(' blocking rows =',i3)
      open(14,file='prep4.rec',status='new',err=280,form='formatted')
      go to 281
  280 ians=1
c      read(14,*) nc1,nr1
c      if(nc1.eq.nc.and.nr1.eq.nr) go to 282
c      rewind(14)
c      call bell
      call aski4l('ok to overwrite prep4.rec?',ians,ierr)
      if(ians.eq.0) stop 'user exit'
      open(14,file='prep4.rec',status='unknown',form='formatted')
  281 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
      write(14,300)nc,nr,prjnam,cmbl
  300 format(2i5/a50/2e16.8)
  282 close(14)

c  get output file name
        ofile=' '
        ofile=ifile(1:index(ifile,'.'))//'pr4'
  301   call askc('Enter output grid name',ofile,ierr)
        if(ierr.eq.-2) then
          close(10)
          go to 401
        endif
c
        call gopen(18,ofile,'new','write',ierr)
        if(ierr.ne.0) then
          print*,'Error - try again'
          go to 301
        endif
        call gopen(11,ofile,'scratch','readwrite',ierr)
        call gheader('w',11,id,ncol,nr,xo,dx,yo,dy,ierr)
c  extend grid to right
      print*,'extending grid to the right'
      call extend(10,11,nc,nr,nz,nr,ncol)
      rewind(11)
c  transpose grid
      print*,'transposing extended grid'
      call transpos(11,12)
      close(11)
      rewind(12)
        call gheader('r',12,id,nr,ncol,yo,dy,xo,dx,ierr)
        call gopen(14,ofile,'scratch','readwrite',ierr)
        call gheader('w',14,id,nrow,ncol,yo,dy,xo,dx,ierr)
c  extend grid to right
      print*,'extending transposed grid to the right'
      call extend(12,14,nr,ncol,nz,nc,nrow)
      close(12)
        call gopen(15,ofile,'scratch','readwrite',ierr)
c  smooth first extention
      print*,'smoothing across first extension'
      call smooth(14,15,nc)
      close(14)
      rewind(15)
c  transpose back
      print*,'transposing back'
      call transpos(15,16)
      close(15)
        call gopen(17,ofile,'scratch','readwrite',ierr)
c  smooth second extension
      print*,'smoothing across second extension'
      call smooth(16,17,nr)
      close(16)
c      rewind(17)
      print*,'centering data'
      call center(17,18,nc,nr)
      close(17)
      close(18)
        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)
      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
      return 
      end
c***********************************************************************
      subroutine extend(in,io,nc,nr,nz,nrow,ncol)
c extend grid to the right using prediction filters
c   in:   input unit number
c   io:   output unit number
c   nc:   columns in input grid
c   nr:   rows in input grid
c   nz:   1
c   nrow: last row to be used in designing the prediction filter
c   ncol: columns in the extended grid
      dimension we(16),wesav(16)
      dimension g(2048)
c  First row.
        call grow('r',in,1,g,nc,ierr)
      call predfilt(nc,ncol,g,16,we)
        do 11 i=1,16
   11   wesav(i)=we(i)
c  Rows down to nrow.
        do 204 j2=2,nrow
        call grow('r',in,j2,g,nc,ierr)
      call predfilt(nc,ncol,g,16,we)
        do 12 i=1,16
   12   wesav(i)=wesav(i)+we(i)
204     continue
      do 13 i=1,16
   13 wesav(i)=wesav(i)/float(nrow)
      rewind(in)
      read(in) id
      do 14 j=1,nr
        call grow('r',in,j,g,nc,ierr)
      call filter(nc,ncol,g,16,wesav)
        call grow('w',io,j,g,ncol,ierr)
   14 continue
      return
      end
c *****************************
      subroutine predfilt(n,lc,a,nco,we)
c determine prediction filter coefs we from profile a
      dimension a(lc),we(nco)
      dimension aa(2048),v(2048)
      real*8 ap,xind
      data ddval/1.e30/
        pi=3.1415927
      if(lc.gt.2048) stop 'too big'
      do 20 i=1,n
      if(a(i).gt.ddval) stop 'dvals found'
      aa(i)=a(i)
   20 v(i)=a(i)
      we(1)=1.
      do 23 j=2,nco
      ap=0.
      xind=0.
      do 21 i=j,n
      ap=ap+aa(i)*aa(i)+v(i-j+1)*v(i-j+1)
   21 xind=xind+aa(i)*v(i-j+1)
      if(ap.eq.0.) stop 'error: ap=0'
      rc=-2.*xind/ap
      do 22 i=j,n
      temp=aa(i)
      aa(i)=aa(i)+rc*v(i-j+1)
   22 v(i-j+1)=v(i-j+1)+rc*temp
      we(j)=0.
      jh=(j+1)/2
      do 23 i=1,jh
      k=j-i+1
      temp=we(k)+rc*we(i)
      we(i)=we(i)+rc*we(k)
   23 we(k)=temp
      return
c *****************************
      entry filter(n,lc,a,nco,we)
c extend profile
      lp=n+1
      do 27 i=1,n
   27 aa(i)=a(n-i+1)
      do 28 i=lp,lc
      a(i)=0.
      aa(i)=0.
      do 28 j=2,nco
      a(i)=a(i)-we(j)*a(i-j+1)
   28 aa(i)=aa(i)-we(j)*aa(i-j+1)
      do 29 i=lp,lc
      w=pi*float(i-lp)/float(lc-lp)
      w=(cos(w)+1.)/2.
      a(i)=w*a(i)+(1.-w)*aa(lc-i+lp)
   29 continue
      return
      end
c *****************************
      subroutine transpos(in,io)
      dimension buf(50000)
c      dimension id(14)
      character id*56
        call gheader('r',in,id,nc,nr,xo,dx,yo,dy,ierr)
        call gopen(io,'dummy','scratch','readwrite',ierr)
        call gheader('w',io,id,nr,nc,yo,dy,xo,dx,ierr)
      n=int((float(nr) / 20.) + .999)
      m=int((float(nc) / 20.) + .999)
      n10 = n*20
      m10 = m*20
      nm = n10 * m10
      open(13,access='direct',status='scratch',recl=1600)
      call rio(nc,buf,m,n,1,1,m10,-1,in)
      ic=1
      do 10 ibuf = 1, m
      loc = ibuf
      call rio(nr, buf, n, 0, loc, m, n10, 1, in)
      ndx = 1
      do 20 i = 1, 20
      call rowio2(nr, buf(ndx), 0, io, io, iend)
      if (ic .eq. nc) goto 99
      ic = ic + 1
   20 ndx = ndx + n10
   10 continue
   99 continue
      close(13)
      return
      end
c *****************************
      subroutine rio(na, z[huge], nblk, nbuf, loc2, idl, n10, iop, in)
      dimension w(400)
      dimension z(1), is(2)
cv    data dval/'37777677777'o/
      data dval / 1.0e+38 /
      loc = loc2
      assign 999 to ir
      if (iop) 1, 2, 3
    1 iflag = 1
      assign 99 to ir
      is(1) = na + 1
      is(2) = 1
      assign 10 to ig
      if (n10 .eq. na) goto 15
      assign 12 to ig
   15 do 100 ibuf = 1, nbuf
      ii = 1
      do 10 irow = 1, 20
      goto (11, 12), iflag
   11 call rowio2(na, z(ii), -1, in, in, ie)
      if (ie .eq. 1) goto 14
      ii = ii + n10
      goto ig
   14 iflag = 2
   12 kk = (irow - 1) * n10
      do 13 k = kk + is(iflag), kk + n10
   13 z(k) = dval
   10 continue
      goto 2
   99 continue
  100 continue
      loc1 = loc - 1
c
      return 
    2 do 30 iblk = 1, nblk
      lw = 1
      ls = ((iblk - 1) * 20) + 1
      do 21 j = 1, 20
      l2 = ls
      l = lw
      do 20 i = 1, 20
      w(l) = z(l2)
      l = l + 20
   20 l2 = l2 + 1
      lw = lw + 1
   21 ls = ls + n10
      write(13, rec=loc) w
   30 loc = loc + idl
c
      goto ir
    3 do 40 iblk = 1, nblk
      read(13, rec=loc) w
      loc = loc + idl
      l = 1
      ls = ((iblk - 1) * 20) + 1
      do 51 j = 1, 20
      l2 = ls
      do 50 i = 1, 20
      z(l2) = w(l)
      l = l + 1
   50 l2 = l2 + 1
   51 ls = ls + n10
   40 continue
  999 return 
      end
c *****************************
      subroutine rowio2(n, z[huge], iop, idev, jdev, iend)
c  WHERE IOP<0 READ; IOP=0 WRITE; IOP>0 READ&WRITE
      dimension z(n)
      iend = 0
      if (iop) 1, 2, 1
    1 read(idev, end=10) xo, z
      if (iop) 9, 9, 2
    2 write(jdev) xo, z
    9 return
   10 iend = 1
      return
      end
c *****************************
      subroutine smooth(in,io,nr)
c smooth across rows
c   in:   input unit number
c   io:   output unit number
c   nr:   number of rows to skip
      dimension g(2048),h(2048)
c      dimension id(14)
      character id*56
        pi=3.1415927
      rewind(in)
        call gheader('r',in,id,nc,nrow,xo,dx,yo,dy,ierr)
        call gheader('w',io,id,nc,nrow,xo,dx,yo,dy,ierr)
      do 10 j=1,nr
        call grow('r',in,j,g,nc,ierr)
        call grow('w',io,j,g,nc,ierr)
  10  continue
      do 30 j=nr+1,nrow
c new
        call grow('r',in,j,g,nc,ierr)
      k=min(j-nr,nrow-j+1)
      do 20 i=1,k
      h(i)=g(i)
      do 19 ii=1,k
      kk=mod(nc+i-ii,nc)
      if(kk.eq.0) kk=nc
   19 h(i)=h(i)+(k-ii+1)*(g(i+ii)+g(kk))/float(k+1)
   20 h(i)=h(i)/float(k+1)
c
      do 22 i=k+1,nc-k
      h(i)=g(i)
      do 21 ii=1,k
  21  h(i)=h(i)+(k-ii+1)*(g(i+ii)+g(i-ii))/float(k+1)
  22  h(i)=h(i)/float(k+1)
c
      do 25 i=nc-k+1,nc
      h(i)=g(i)
      do 24 ii=1,k
      kk=mod(i+ii,nc)
      if(kk.eq.0) kk=nc
   24 h(i)=h(i)+(k-ii+1)*(g(i-ii)+g(kk))/float(k+1)
   25 h(i)=h(i)/float(k+1)

   30  call grow('w',io,j,h,nc,ierr)
      return
      end
c *****************************
      subroutine center(in,io,nc,nr)
c smooth across rows
c   in:   input unit number
c   io:   output unit number
c   nc:   original number of columns
c   nr:   original number of rows
      dimension g(2048)
c      dimension id(14)
      character id*56
c        pi=3.1415927
      rewind(in)
        call gheader('r',in,id,ncol,nrow,xo,dx,yo,dy,ierr)
        nx=(ncol-nc)/2
        ny=(nrow-nr)/2
        xo=xo-nx*dx
        yo=yo-ny*dy
        call gheader('w',io,id,ncol,nrow,xo,dx,yo,dy,ierr)
      do 10 j=1,nrow-ny
        call grow('r',in,j,g,ncol,ierr)
  10  continue
      do 30 j=nrow-ny+1,nrow
        call grow('r',in,j,g,ncol,ierr)
        do 25 k=1,nx
        temp=g(ncol)
        do 20 i=ncol,2,-1
        g(i)=g(i-1)
  20  continue
        g(1)=temp
  25  continue
        call grow('w',io,j,g,ncol,ierr)
  30  continue
      rewind(in)
        call gheader('r',in,id,ncol,nrow,xo,dx,yo,dy,ierr)
      do 50 j=1,nrow-ny
        call grow('r',in,j,g,ncol,ierr)
        do 45 k=1,nx
        temp=g(ncol)
        do 40 i=ncol,2,-1
        g(i)=g(i-1)
  40  continue
        g(1)=temp
  45  continue
        call grow('w',io,j,g,ncol,ierr)
  50  continue
      return
      end
c *****************************
      subroutine bell
      character*1 ding
      ding=char(007)
      print *, ding
      return
      end
