c  program regrid
c  input is a grid with a constant increment in x and y,
c  (dx,dy are not necessarily equal)
c  output is a grid with a specified increment, and optionally
c  covering a specified area.
c  the interpolater is a 1-dimensional cubic spline applied first
c  in the y direction, then in x.
c  maximum grid dimension 1500x1500 (tested ok, 2500x2500 may
c  have a problem.)
c
c      Mike Webring
c    Geophysics branch, USGS
c    Lakewood, Colo.
c
        common /regrd/ buf(50000),z(2630)
        common /spl15/ wrk(13150)
        real dxy2(2),xoyo(4)
	character*56 ifile,jfile,kfile
        character id*56,id2*56,p*8
	data kfile/'regrid.tmp'/,inp/5/,p2/'regrid'/,
     1 dval/0.1701412e39/,isw,jsw,ksw/12,13,14/
c
        call askin
        call pfinit('regrid')
        ifile=' '
400     call askc('file to be interpolated',ifile,ierr)
        if(ierr.eq.-2) stop
        call gopen(isw,ifile,'old','read',ierr)
        if(ierr.ne.0) then
          print*,'Error - try again'
          go to 400
        endif
        call gheader('r',isw,id,nc,nr,xo,dx,yo,dy,ierr)
        call gpgm(p)
        write (*, 210)id,p,nc,nr,xo,yo,dx,dy
210	format(/,1x,a56,'  pgm:',a8,/,' ncol=',i5,' nrow=',i5,/,
     1 ' xo,yo=',1p2g15.5,'   dx,dy=',2g15.5)
c
c        jfile=' '
        jfile=ifile(1:index(ifile,'.'))//'rgd'
401      call askc('output file',jfile,ierr)
        if(ierr.eq.-2) then
          close(isw)
          go to 400
        endif
        call gopen(jsw,jfile,'new','write',ierr)
        if(ierr.ne.0) then
          print*,'Error - try again'
          go to 401
        endif
        id2=id
        call askalt
402     call askc('title',id2,ierr)
        if(ierr.eq.-2) then
          close(jsw)
          go to 401
        endif
        call askalt
        dxy2(1)=dx
        dxy2(2)=dy
403     call askf4a('new dx,dy',dxy2,2,ierr)
        if(ierr.eq.-2) go to 402
        dx2=dxy2(1)
        dy2=dxy2(2)
	if(dx*dx2.lt.0.0) stop ' dx,dx2 have differing signs'
	if(dy*dy2.lt.0.0) stop ' dy,dy2 have differing signs'
        ians=0
404     call aski4l('change location or areal coverage?',ians,ierr)
        if(ierr.eq.-2) go to 403
        if(ians.eq.1) then
          xoyo(1)=xo
          xoyo(2)=yo
          xoyo(3)=nc
          xoyo(4)=nr
405       call askalt
          call askf4a('enter new  x origin, y origin, ncol, nrow',xoyo,
     1                 4,ierr)
          call askalt
          if(ierr.eq.-2) go to 404
          xo2=xoyo(1)
          yo2=xoyo(2)
          nc2=xoyo(3)
          nr2=xoyo(4)
        else
	  xo2=xo
	  yo2=yo
	  nc2=int((float(nc-1)*dx)/dx2+1.00001)
	  nr2=int((float(nr-1)*dy)/dy2+1.00001)
	  write ( *,*) ' new grid will be '
	  write (*,*)  nc2,' columns by ',nr2,' rows'
	endif
        if(nc2.gt.2630 .or. nr2.gt.2630
     1 .or. nc.gt.2630 .or. nr.gt.2630) then
	  write ( *,*) '  input col=',nc, '   input row=',nr
	  write ( *,*) ' output col=',nc2,'  output row=',nr2
          write (*,*) ' in&out x&y max dimension is 2630'
	  stop
	endif
        call gheader('w',jsw,id2,nc2,nr2,xo2,dx2,yo2,dy2,ierr)
c
	m = int(float(nc)/20.0+1.00001)
	m20=m*20
	n = int(float(nr)/20.0+1.00001)
	n2=int(float(nr2)/20.0+1.00001)
	n=max(n,n2)
	n20=n*20
	open(unit=ksw,file=kfile,access='direct',form='unformatted',
     1 status='unknown',recl=1600)
	call rio(buf,m,n,1,1,m20,-1,isw,ksw)
c
	ifill=0
	ic=1
	do 100 ibuf=1,m
c  assemble buffer of 20 columns
	loc=ibuf
	call rio(buf,n,0,loc,m,n20,1,isw,ksw)
c  spline in y
	iof=0
	do 255 ir=1,20
	if(ifill.eq.0) then
	  do 260 i=1,nr
	    z(i)=buf(i+iof)
260     continue
	  call splint(nr,yo,dy,z,nr2,yo2,dy2,buf(1+iof))
	  if(ic.eq.nc) ifill=1
	else
	  do 270 i=1,n20
	    buf(i+iof)=dval
270     continue
	endif
	ic=ic+1
	iof=iof+n20
255     continue
c  output to same addresses
	call rio(buf,n,0,loc,m,n20,0,isw,ksw)
100	continue
c
c  spline in x, output final grid rows
	ir=1
	loc=1
	do 280 ibuf=1,n2
	  call rio(buf,m,0,loc,1,m20,1,isw,ksw)
	  iof=1
	  do 290 j=1,20
	    call splint(nc,xo,dx,buf(iof),nc2,xo2,dx2,z)
	    call rowio(nc2,z,0,jsw,jsw,iend)
	    if(ir.eq.nr2) go to 999
	    iof=iof+m20
	    ir=ir+1
290     continue
	  loc=loc+m
280     continue
c
999	continue
	close(isw)
	close(jsw)
	close(ksw)
	stop
	end
c******************************************************************************
	subroutine rio(z[huge],nblk,nbuf,loc2,idloc,lenz,iop,isw,ksw)
c iop=-1 read grid, pad with dval's, and output with implied iop=1 
c    =0 output 20x20 tranposed blocks
c    =1 input  20x20 block
c buffer z is dimensioned lenz x 20
c random access is nblk x nbuf, where 'loc' increments by 'idl'
	dimension hdr(16),w(400),z(1)
	logical more
	data dval/.1701412e39/
	loc=loc2
	assign 999 to iret
c
	if(iop)1,2,3
1	assign 99 to iret
	more=.true.
	itot=0
	rewind(isw)
	read(isw) hdr,nc,nr,nz,xo,dx,yo,dy
	do 100 ibuf=1,nbuf
	do 105 i=1,lenz*20
	  z(i)=dval
105     continue
	ireadz=1
	do 106 irow=1,20
	  if(more) then
	    call rowio(nc,z(ireadz),-1,isw,isw,iend)
	    ireadz=ireadz+lenz
	    if(iend.eq.0) itot=itot+1
	    if(iend.ne.0 .or. itot.ge.nr) more=.false.
	  endif
106     continue
	  go to 2
99	  continue
100	continue
	if(itot.ne.nr) write(*,*) ' found only ',itot,' of ',nr,' rows'
	return
c
c  extract 20x20 block, tranpose, and output
2	do 107 iblk=1,nblk
	iws=1
	izs=(iblk-1)*20+1
	do 108 j=1,20
	  iw=iws
	  iz=izs
	  do 109 i=1,20
	    w(iw)=z(iz)
	    iw=iw+20
	    iz=iz+1
109     continue
	  iws=iws+1
	  izs=izs+lenz
108     continue
	write(ksw,rec=loc) w
	loc=loc+idloc
107     continue
	go to iret,(99,999)
c
c  insert 20x20 block into buffer
3	do 113 iblk=1,nblk
	read(ksw ,rec=loc) w
	iw=1
	izs=(iblk-1)*20+1
	do 114 j=1,20
	  iz=izs
	  do 111 i=1,20
	    z(iz)=w(iw)
	    iw=iw+1
	    iz=iz+1
111     continue
	izs=izs+lenz
114     continue
	loc=loc+idloc
113     continue
c
999	return
	end
c******************************************************************************
	subroutine splint(nz,txo,tdx,z[huge], nz2,txo2,tdx2,zint[huge])
c  spline interpolation on defined intervals, no extrapolation
	dimension d(2),z(nz),zint(nz2)
        common /spl15/ a(2630),b(2630),c(2630),
     1 w1(2630),w2(2630)
	data tol/.001/, d/0.,0./, dval,test/0.1701412e39,1.0e29/
c
	if(tdx.eq.0.0 .or. tdx2.eq.0.0) go to 999
	do 110 i=1,nz2
	  zint(i)=dval
110     continue
	xo = txo*sign(1.0,tdx)
	dx = abs(tdx)
	xo2= txo2*sign(1.0,tdx2)
	dx2= abs(tdx2)
	dxtol=amin1(dx*tol,dx2*tol)
	x1c=xo-dx
	x2c=xo2-dx2
	xend2 = float(nz2)*dx2 + x2c
c
	ie=0
	xs = -test
	xe = -test
c  cycle the output array
	do 100 ioutp = 1, nz2
	xint = float(ioutp)*dx2 + x2c
	if (xint.lt.xs) go to 100
	if (xint.ge.xe) then
c  get start and ending indices of next input segment
5	next=ie+1
	if(next.ge.nz) go to 999
	do 120 is = next, nz
	  if(z(is).lt.test) then
	    do 130 ie = is+1, nz
	      if(z(ie).gt.test) go to 10
130     continue
	    ie = nz+1
	    go to 10
	  endif
120     continue
	go to 999
10	ie=ie-1
	ni=ie-is+1
	xs = float(is)*dx + x1c - dxtol
	if (xs.gt.xend2) go to 999
	xe = float(ie)*dx + x1c + dxtol
	if (xe.lt.xint) go to 5
	if ( ni.lt.2 ) go to 100
	  idata=is
	  x1 = xs - dxtol
	  x2 = xs + dx + dxtol
	  if (ni.gt.2) then
	    call spln1(ni,dx,xdum,z(is),a,b,c,0,d,w1,w2)
	  else
	    a(1)=(z(ie)-z(is))/dx
	    b(1)=0.0
	    c(1)=0.0
	  endif
	endif
c
c  find input interval that contains xint
	if(xint.ge.x1 .and. xint.lt.x2) go to 50
	idata1 = idata + 1
	if(xint.ge.x2) then
	  if(idata.eq.ie) go to 100
	  do 140 idata = idata1, ie-1
	    x1 = float(idata)*dx + x1c
	    x2 = x1 + dx + dxtol
	    x1 = x1 - dxtol
	    if(xint.ge.x1 .and. xint.lt.x2) go to 50
140     continue
	else
	  if(idata.eq.is) go to 100
	  write(*,*)'splint: backtracking'
	  idata1 = idata - 1
	  do 150 idata = idata1, is, -1
	    x1 = float(idata)*dx + x1c
	    x2 = x1 + dx + dxtol
	    x1 = x1 - dxtol
	    if(xint.ge.x1 .and. xint.lt.x2) go to 50
150     continue
	endif
	go to 100
c
c  cubic interpolation
50	xd=xint-x1
	j=(idata-is)+1
	zint(ioutp)=z(idata)+((c(j)*xd+b(j))*xd+a(j))*xd
c
100	continue
999	return
	end
c******************************************************************************
	subroutine spln1(m,h,x,y[huge],a,b,c,t,d,p,s)
c	coded by W. L. Anderson, U.S.G.S.
	dimension x(m),y(m),a(m),b(m),c(m),d(2),p(m),s(m)
	real mul
	integer t
	if(t.lt.0 .or. t.gt.1 .or. h.lt.0. .or. m.lt.3) go to 999
	n=m-1
	if(t.eq.0) go to 20
c  1ST DERIV BOUNDS GIVEN
	ne=n-1
	if(h)999,11,1
c EQUAL SPACING H>0  T=1
1	hh=3./h
	do 2 i=1,ne
	b(i)=4.
	c(i)=1.
	a(i)=1.
2	p(i)=hh*(y(i+2)-y(i))
	p(1)=p(1)-d(1)
	p(ne)=p(ne)-d(2)
c SOLN OF TRIDIAGONAL MATRIX ORDER NE
3	c(1)=c(1)/b(1)
	p(1)=p(1)/b(1)
	do 4 i=2,ne
	mul=1./(b(i)-a(i)*c(i-1))
	c(i)=mul*c(i)
4	p(i)=mul*(p(i)-a(i)*p(i-1))
c OBTAIN SPLINE COEF
	a(ne+t)=p(ne)
	i=ne-1
5	a(i+t)=p(i)-c(i)*a(i+t+1)
	i=i-1
	if(i.ge.1) go to 5
	if(t.eq.0) go to 6
	a(1)=d(1)
	a(m)=d(2)
6	if(h.eq.0.) go to 14
	hh=1./h
	do 7 i=1,n
	i1=i+1
	mul=hh*(y(i1)-y(i))
	b(i)=hh*(3.*mul-(a(i1)+2.*a(i)))
7	c(i)=hh*hh*(-2.*mul+a(i1)+a(i))
	return
c UNEQUAL SPACING H=0  T=1
11	do 12 i=1,n
	i1=i+1
12	s(i1)=x(i1)-x(i)
	do 13 i=1,ne
	i1=i+1
	i2=i+2
	b(i)=2.*(s(i1)+s(i2))
	c(i)=s(i1)
	a(i)=s(i2)
	sq1=s(i1)*s(i1)
	sq2=s(i2)*s(i2)
13	p(i)=3.*(sq1*(y(i2)-y(i1))+sq2*(y(i1)-y(i)))/(s(i1)*s(i2))
c       B(I)=2.*(S(I+1)+S(I+2))
c       C(I)=S(I+1)
c       A(I)=S(I+2)
c13     P(I)=3.*(S(I+1)**2*(Y(I+2)-Y(I+1))+S(I+2)**2*(Y(I+1)-Y(I)))/
c       1  (S(I+1)*S(I+2))
	p(1)=p(1)-s(3)*d(1)
	p(ne)=p(ne)-s(n)*d(2)
	go to 3
14	do 15 i=1,n
	i1=i+1
	hh=1./s(i1)
	hs=hh*hh
	mul=(y(i1)-y(i))*hs
	b(i)=3.*mul-(a(i1)+2.*a(i))*hh
15	c(i)=-2.*mul*hh+(a(i1)+a(i))*hs
	return
c  2ND DERIV BOUNDS GIVEN
20	ne=n+1
	if(h)999,31,21
c  EQUAL SPACING H>0  T=0
21	hh=3./h
	do 22 i=2,n
	b(i)=4.
	c(i)=1.
	a(i)=1.
22	p(i)=hh*(y(i+1)-y(i-1))
	b(1)=2.
	b(ne)=2.
	c(1)=1.
	c(ne)=1.
	a(ne)=1.
	p(1)=hh*(y(2)-y(1))-.5*h*d(1)
	p(ne)=hh*(y(m)-y(n))+.5*h*d(2)
	go to 3
c  UNEQUAL SPACING H=0  T=0
31	do 32 i=1,n
	i1=i+1
32	s(i1)=x(i1)-x(i)
	n1=n-1
	do 33 i=1,n1
	i1=i+1
	i2=i+2
	b(i1)=2.*(s(i1)+s(i2))
	c(i1)=s(i1)
	a(i1)=s(i2)
	sq1=s(i1)*s(i1)
	sq2=s(i2)*s(i2)
33	p(i1)=3.*(sq1*(y(i2)-y(i1))+sq2*(y(i1)-y(i)))/(s(i1)*s(i2))
c       B(I+1)=2.*(S(I+1)+S(I+2))
c       C(I+1)=S(I+1)
c       A(I+1)=S(I+2)
c33     P(I+1)=3.*(S(I+1)**2*(Y(I+2)-Y(I+1))+S(I+2)**2*(Y(I+1)-Y(I)))/
c       1 (S(I+1)*S(I+2))
	b(1)=2.
	b(ne)=2.
	c(1)=1.
	c(ne)=1.
	a(ne)=1.
	p(1)=3.*(y(2)-y(1))/s(2)-.5*s(2)*d(1)
	p(ne)=3.*(y(m)-y(n))/s(m)+.5*s(m)*d(2)
	go to 3
999	m=-abs(m)
	return
	end
c******************************************************************************
	subroutine rowio(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
c******************************************************************************
        function noyes()
	character inp*8
	icnt = 0
	noyes = -1
1	inp = ' '
	read (5,2,end=5,err=5) inp
2	format(a8)
	do 20 n = 1,8
	  if (inp(n:n).ne.' ') then
	    if (inp(n:n).eq.'y' .or. inp(n:n).eq.'Y') noyes = 1
	    if (inp(n:n).eq.'n' .or. inp(n:n).eq.'N') noyes = 0
	    if (noyes.gt.-1) return
	    go to 5
	  endif
20      continue
5	if (icnt.lt.3) then
	  icnt = icnt + 1
	  write (*,6)
6	  format(' y or n:'$)
	  go to 1
	endif
	write(*,*)' count exceeded, answering no'
	noyes = 0
	return
	end


