c                 program rgrd
c
c   linear interpolater with constant gradients along cell edges.
c   or one-dimensional spline in both directions
c
       subroutine rgrd (buf,igrid,speed,ierr)
       dimension buf(1)
       character*1 speed
       character*8 kfile
       character*64 title
	common /igrid/ nc,   nr,  xo,  dx,  yo,  dy
	common /ogrid/ nc2, nr2, xo2, dx2, yo2, dy2
       common/scalef/iwc0,jwc0,nxhpix,nyhpix,xpix,ypix,xscr,yscr
       logical linear
       data dval/1.e38/,kfile/'rgrd.tmp'/

         ierr = 0
         rewind(igrid)
         read(igrid,err=90,iostat=icheck) title,nc,nr,nz,xo,dx,yo,dy

	  xo2 = xo
	  yo2 = yo
	  nc2 = xscr 
	  nr2 = yscr
         dx2=((nc-1)*dx)/(xscr-1)
         dy2=((nr-1)*dy)/(yscr-1)
         open(igrid+1,status='scratch',form='unformatted',
     &        err=90,iostat=icheck)
         write(igrid+1,err=90,iostat=icheck) title,nc2,nr2,nz,xo2,
     &   dx2,yo2,dy2

c  adjust if close to integer subdivisor
c  call intdx( dx, dx2 )
c  call intdx( dy, dy2 )

         isw=igrid
         jsw=igrid+1
         ksw=igrid+2

	if ( nc2 .gt. 800  .or.  nr2 .gt. 800 .or.
     1    nc  .gt. 800  .or.  nr  .gt. 800.or.speed.eq.'f' ) then
          linear = .true.
       else
          linear = .false.
       endif

c  interpolate

       if(linear) then
	  n2 =  1 + nc
	  n3 = n2 + nc
	  n4 = n3 + nc
	  call lin2d( isw, jsw, dval, buf(1), buf(n2), buf(n3),
     1             buf(n4), ierr )
       else

	open(unit=ksw,file=kfile,access='direct',form='unformatted',
     1 status='unknown',recl=1600,err=90,iostat=icheck)
	call spln2d( isw, jsw, ksw, dval,buf,ierr )
       endif

	return
90     print *,'error no. ',icheck,' in subroutine rgrd'
       ierr = 1
       close(ksw,status='delete')
       return
	end
c***********************************************************************
	subroutine intdx( dx, dx2 )
c  assume the user was shooting for integral subdivision
c  if dx is an integral multiple of dx2 within .001
	if ( abs( dx * dx2 ) .lt. 1e-30 ) return
	t1   = abs ( dx/dx2 )
	t2   = abs ( dx2/dx )
	t    = amax1 ( t1, t2 )
	tint = aint( t + .001 )
	even = abs(  t - tint )
	if ( even .le. 0.001 ) then
	  if ( dx2 .lt. dx ) then
	    dx2 = dx / tint	
	    else
	    dx2 = dx * tint
	  endif
	endif
	return
	end
c***********************************************************************
	subroutine lin2d( isw, jsw, dval, b1, b2, flag, r, ierr )
c  grid linear interpolation assuming constant edge gradients
c  on four sides of the input cell.
c  f(ix+dx,iy+dy) = f(ix,iy) + a*dx + b*dy + c*dx*dy.
c  general placement of input and output grids.
c
	common /igrid/ nc,  nr,  xo,  dx,  yo,  dy
	common /ogrid/ nc2, nr2, xo2, dx2, yo2, dy2
       common/scalef/iwc0,jwc0,nxhpix,nyhpix,xpix,ypix,xscr,yscr
       character*1 ctype
       common/captur/icapt,ctype
	dimension b1(nc), b2(nc), flag(nc), r(nc2)
c
c  switch internal coordinates to right cartesian
	if ( dx  .lt. 0 )  xo = -xo
	if ( dy  .lt. 0 )  yo = -yo
	if ( dx2 .lt. 0 ) xo2 = -xo2
	if ( dy2 .lt. 0 ) yo2 = -yo2
	dx  = abs( dx  )
	dx2 = abs( dx2 )
	dy  = abs( dy  )
	dy2 = abs( dy2 )
c
c  initialize output row
	do i = 1, nc2
	  r(i) = dval
	enddo
c
c  register input and output grids
c
	rewind isw
	read( isw )
	irow  = 0
	nout  = 0
	ywork = yo2
c	
c  output dval rows until output catches up to input.
	do while ( ywork .lt. yo )
	  write( jsw,err=900,iostat=icheck ) ywork, r
	  nout  = nout + 1
	  yout  = yo2  + dy2 * float( nout - 1 )
	  ywork = yout + dy2
	  if ( nout .eq. nr2 ) then
	    print *, ' output blank grid: yend, ynext=', yout, ywork
           ierr = 1
	    return
	  endif
	enddo
c
c  read input until yl & yu bracket ywork
	yu = -1.e35
	do while ( ywork .gt. yu )
	  read( isw, end = 777 ) xx, b1
	  irow = irow + 1
	  yl   = yo + dy * float( irow - 1 )
	  yu   = yl + dy
	enddo
	read( isw ) xx, b2
	irow = irow + 1
	yu   = yo + dy * float( irow - 1 )
c
	test = ( yu - yl ) / dy - 1.0
	if ( abs( test ) .gt. 0.1 ) go to 800
	if ( ywork .gt. yu   .or.   ywork .lt. yl ) go to 820
	call setflg( nc, b1, b2, flag )
c
c  interpolation loop
c
	istop = 0
100	ywork = yo2 + dy2 * float( nout )
c
c       if ( ywork .gt. yu ) then 
c*******THIS STATEMENT CHANGED BY R.G.
        if ( ywork .gt. yu .and. nout .lt. (nr2-1) ) then 
c  update buffer
110	  do i = 1, nc
	    b1(i) = b2(i)
	  enddo
	  read( isw, end=120, iostat=ios ) xx, b2
120	  irow = irow + 1
	  yu   = yo + dy * float( irow - 1 )
	  yl   = yu - dy
	  if ( ios .eq. 0 ) then
	    if ( ywork .gt. yu ) go to 110
	    call setflg( nc, b1, b2, flag )
	    else
	    istop = 1
	  endif
	endif
c
	if ( istop .eq. 0 ) then
c  interpolate row
	  call lin1d( dval, yl, ywork, b1, b2, flag, r )
	  else
c  output blank rows to end of output grid
	  do i = 1, nc2
	    r(i) = dval
	  enddo
	  do while ( 1 .eq. 1 )
	    write (jsw,err=900,iostat=icheck) ywork, r
	    nout  = nout + 1
           if(nr2.gt.99) then
             call perout(nout,nr2,1,1,'REGRIDDED')
           endif
	    ywork = yo2 + dy2 * float( nout )
	    if ( nout .eq. nr2 ) return
	  enddo
	endif
c
c  output row
	write (jsw,err=900,iostat=icheck) ywork, r
	nout = nout + 1
           if(nr2.gt.99) then
             call perout(nout,nr2,1,1,'REGRIDDED')
           endif
	if ( nout .lt. nr2 ) go to 100
	if ( nout .gt. nr2 ) go to 890
	return
c
777	print *, ' input eof - check grid origins'
       ierr = 1
	return
800	print *, 'lin2d: buffer init,  test =', test
820	print *, 'lin2d: irow, yl, yu' , irow, yl, yu
890	print *, ' nout, ywork', nout, ywork
       ierr = 1
900    print *,'error no. ',icheck,' writing output in subroutine rgrd'
       ierr = 1
	return
	end
c***********************************************************************
	subroutine lin1d( dval, yl, ywork, b1, b2, flag, r )
c  linear interpolate profile
	common /igrid/ nc,  nr,  xo,  dx,  yo,  dy
	common /ogrid/ nc2, nr2, xo2, dx2, yo2, dy2
	dimension      b1(nc), b2(nc), flag(nc), r(nc2)
	do i = 1, nc2
	  r(i) = dval
	enddo
c
	rdx  = 1.0 / dx
	dely = ( ywork - yl ) / dy
	ix   = int( ( xo2 - xo ) / dx ) - 1
	if ( ix .le. 0  ) ix = 1
	if ( ix .ge. nc ) go to 999
	xu   = xo + dx  * float( ix )
	xl   = xu - dx
c
	do 100 i2 = 1, nc2
	  x2 = xo2 + dx2 * float( i2 - 1 )
c
	  if ( x2 .lt. xl ) go to 100
c10	  if ( x2 .gt. xu ) then
c********THIS STATEMENT CHANGED BY R.G.
10	  if ( x2 .gt. xu .and. i2 .ne. nc2 ) then
11	    ix = ix + 1
	    xu = xo + dx * float( ix )
	    xl = xu - dx
	    if ( ix .ge. nc ) go to 999
	    if ( x2 .gt. xu ) go to 11
	  endif
c
	if ( flag(ix) .eq. 0.  ) then
	  delx  = ( x2 - xl ) * rdx
	  ixu   =   ix + 1
	  a     =   b1(ixu) - b1(ix)
	  b     =   b2(ix)  - b1(ix)
	  c     =   b2(ixu) - b2(ix) - a
	  r(i2) =   b1(ix)  +  delx * a + dely * ( b + delx * c )
	endif

100	continue
c
999	return
	end
c***********************************************************************
	subroutine setflg( nc, b1, b2, flag )
	dimension b1(nc), b2(nc), flag(nc)
	do i = 1, nc
	  flag(i) = 0.
	enddo
c
	if ( b1(1) .gt. 1.e29  .or.  b2(1) .gt. 1.e29 ) flag(1) = 1.
	do i = 2, nc
	  if ( b1(i) .gt. 1.e29  .or.  b2(i) .gt. 1.e29 ) then
	    flag(i)   = 1.
	    flag(i-1) = 1.
	  endif
	enddo
	return
	end
c***********************************************************************
	subroutine spln2d( isw, jsw, ksw, dval,buf, ierr )
       dimension z(800),buf(1)
	common /spl15/ wrk(4000)
	common /igrid/  nc,  nr,  xo,  dx,  yo,  dy
	common /ogrid/ nc2, nr2, xo2, dx2, yo2, dy2
	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
       numspl = 20 * m + 20 * n2
	call rio( buf, m, n, 1, 1, m20, -1, isw, ksw,ierr )
       if(ierr.ne.0) return
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,ierr)
       if(ierr.ne.0) return
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),ierr)
         if(ierr.ne.0) return
         call perout(ic,numspl,1,1,'REGRIDDED')
	  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,ierr)
       if(ierr.ne.0) return
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,ierr)
         if(ierr.ne.0) return
	  iof=1
	  do 290 j=1,20
	    call splint(nc,xo,dx,buf(iof),nc2,xo2,dx2,z,ierr)
           if(ierr.ne.0) return
           call perout(ic+ir,numspl,1,1,'REGRIDDED')
	    call rrowio(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
       call perout(numspl,numspl,1,1,'REGRIDDED')
	close(ksw,status='delete')
	return
	end
c******************************************************************************
	subroutine rio(z,nblk,nbuf,loc2,idloc,lenz,iop,isw,ksw,ierr)
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/1.0e38/
	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 rrowio(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) then
         write(*,*) ' found only ',itot,' of ',nr,' rows'
         ierr = 1
       endif
	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,err=900,iostat=icheck) 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
900    print *,'error no. ',icheck,' writing output in subroutine rgrd'
       ierr = 1
	return
	end
c******************************************************************************
	subroutine splint(nz,txo,tdx,z, nz2,txo2,tdx2,zint,ierr)
c  spline interpolation on defined intervals, no extrapolation
	dimension d(2),z(nz),zint(nz2)
	common /spl15/ a(800),b(800),c(800),
     1 w1(800),w2(800)
	data tol/.001/, d/0.,0./, dval,test/1.0e38,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
c	if (xint.ge.xe) then
c********THIS STATEMENT CHANGED BY R.G.
	if (xint.ge.xe.and.ioutp.ne.nz2) 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
c********THIS STATEMENT INSERTED BY R.G.
         if(ioutp.eq.nz2.and.(abs(xint-x2).le.000001)) go to 50
	  if(idata.eq.ie) go to 100
	  do 140 k = idata1, ie-1
           idata = k
	    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
        idata = ie - 1
	else
	  if(idata.eq.is) go to 100
	  write(*,*)'splint: backtracking'
         print *,' run again using linear interpolation'
         print *,' by pressing ALT F5 or ALT F6'
         ierr = 1
         return
	  idata1 = idata - 1
	  do 150 k = idata1, is, -1
           idata = k
	    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
        idata = is 
	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,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 rrowio(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
