c  PREP2.FOR
c
c  2D version of PREP.  Detrends and extrapolates 1-row grid profile for FFT.
c  Constants stored in prep2.rec file, to be used by program DE_PREP2 to trim
c  the profile back and (optionally) restore the linear trend.
c
c  Extrapolation is between ncol and new-ncol.  FFT is fastest if new-ncol
c  is a power of 2, e.g. 64,128,256,...2048.  2048 is max allowed.
c
c  Lin Cordell
c  22 October 92
c
      character*50 ifile,ofile
      dimension f(2000),id(14),pgm(2)
c      data dval/0.170412e+39/
      dval2=1.e20
C      print*,'Enter input 1-row-grid profile file name:'
C      print 1
C1     format(' *'$)
C      read(5,100)ifile
C100   format(a50)
        ifile=' '
        call askin
400     call askc('Enter input 1-row-grid profile file name',ifile,ierr)
        if(ierr.lt.0) stop
C      print*,'Enter output file name:'
C      print 1
C      read(5,100)ofile
        ofile=' '
        ofile=ifile(1:index(ifile,'.'))//'prp'
401     call askc('Enter output grid name',ofile,ierr)
        if(ierr.eq.-2) go to 400
C      open(10,file=ifile,form='unformatted',status='old')
C      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
        call gopen(10,ifile,'old','read',ierr)
        if(ierr.ne.0) stop 'Problem opening file'
        call gheader('r',10,id,ncol,nrow,x0,dx,y0,dy,ierr)
c      if(ncol.gt.2048.or.nz.ne.1.or.dx.eq.0.e0)then
      if(ncol.gt.2048.or.dx.eq.0.e0)then
        print*,'Too big.'
        stop
      endif
C      open(11,file=ofile,form='unformatted',status='new')
        call gopen(11,ofile,'new','write',ierr)
C      print*,'Enter new id:'
C      print 1
C      read(5,203)(id(i),i=1,14)
C203   format(14a4)
C201   print*,'Enter new ncol:'
C      print 1
C      read*,nncol
201     nncol=ncol
        call aski4('Enter new ncol',nncol,ierr)
        if(ierr.eq.-2) go to 401
      if(nncol.lt.ncol.or.nncol.gt.2048)then
        print*,'Cant handle it.'
        go to 201
      endif
C      write(11)id,pgm,nncol,nrow,nz,x0,dx,y0,dy
        call gheader('w',11,id,nncol,nrow,x0,dx,y0,dy,ierr)
C      call rowio(f,y,1,ncol,10,1)
        call grow('r',10,1,f,ncol,ierr)
c  Detrend
      sx=0.
      sf=0.
      sxx=0.
      sfx=0.
      rn=ncol
      do 102 i=1,ncol
      if(f(i).gt.dval2)then
        print*,'Dvals not allowed'
        stop
      endif
      x=i*dx
      sx=sx+x
      sf=sf+f(i)
      sxx=sxx+x*x
      sfx=sfx+f(i)*x
102   continue
      den=rn*sxx-sx*sx
      if(den.eq.0.)then
        print*,'Jacobian = 0, no solution possible'
        go to 99
      endif
      a=(sf*sxx-sx*sfx)/den
      b=(rn*sfx-sf*sx)/den
      print*,'A =',a,', B =',b
      open(12,file='prep2.rec',status='new',err=280,form='formatted')
      go to 281
  280 ians=1
      call aski4l('ok to overwrite prep2.rec?',ians,ierr)
      if(ians.eq.0) stop 'user exit'
      open(12,file='prep.rec',status='unknown',form='formatted')
  281 write(12,*)a,b,ncol
      close(12)
      do 103 i=1,ncol
      ri=i
      f(i)=f(i)-a-b*ri*dx
103   continue
c  Extrapolate with sine taper.
      if(nncol.eq.ncol)go to 105
c  Sine taper.
      pi=3.1415927
      arg=pi/(2.*(nncol-ncol+1))
      r=f(ncol)-f(1)
      do 104 i=1,nncol-ncol
104   f(ncol+i)=f(ncol)-r*sin(i*arg)
c  Smooth the near-extrapolated section with moving average.
      if(nncol-ncol.gt.4)then
        f(ncol-3)=(f(ncol-4)+f(ncol-3)+f(ncol-2))/3.
        f(nncol-1)=(f(nncol-2)+f(nncol-1)+f(nncol))/3.
        do 52 i=ncol-2,nncol-2
        s=0
        do 51 j=-2,2
        s=s+f(i+j)
51      continue
        f(i)=.2*s
52      continue
      endif
C105   call rowio(f,y,1,nncol,11,2)
105     call grow('w',11,1,f,nncol,ierr)
      go to 99
90    print*,'Problem opening file',ifile
99    stop
      end
      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
      subroutine splnon(f,h,p,n)
c
c  finds the horizontal derivative for irregularly-spaced data.
c
      dimension f(n),h(n),p(n),c(1000),d(1000)
      c(1)=.5
      d(1)=1.5*(f(2)-f(1))/h(2)
      n1=n-1
      do 10 i=2,n1
      dnm=2.*(h(i)+h(i+1))-h(i+1)*c(i-1)
      c(i)=h(i)/dnm
10    d(i)=(3.*(h(i)*(f(i+1)-f(i))/h(i+1)+h(i+1)*(f(i)-f(i-1))/h(i))-
     1     h(i+1)*d(i-1))/dnm
      p(n)=(3.*(f(n)-f(n1))/h(n)-d(n1))/(2.-c(n1))
      do 9 i=1,n1
      k=n-i
9     p(k)=d(k)-c(k)*p(k+1)
      return
      end
