subroutine intrp (itab, xtab, ytab, ntab, x, y, ierr)
c
c Interpolator using CODIM1 algorithm which is admittedly
c obscure  but works well.
c
c       itab - a label between 1 and 20 to identify the table and its
c              most recent search index
c       xtab - array of length ntab containing the x-values
c       ytab -                                     y-values
c       ntab - number of x,y pairs in the table
c          x - independent for which a y-value is desired
c          y - returned interpolated (or extrapolated) value
c       ierr - =0 for ok, -1 for extrapolation
c
        double precision xtab(ntab), ytab(ntab), x, y
        integer itab, ierr, index
        double precision t(4), u(4)
c       integer savind
c       data savind/-1/
c
c----- Only 1 pt in table
        if (ntab .eq. 1) then
                y = ytab(1)
                ierr = 0
                return
        endif
c
c-----
c Locate search index
        call srch (itab, x, xtab, ntab, index, ierr)
c       if (index .eq. savind) go to 2000
c       savind = index
c
c-----
c Set interpolator index flags
        i1 = 2
        i2 = 3
        iload = max0 (index-2, 1)
c
        if (ntab .gt. 2) then
                if (index.eq.   2) i2 = 4
c
                if (index.eq.ntab) i1 = 1
        endif
c
        if (index.gt.2 .and. index.lt.ntab) then
                i1 = 1
                i2 = 4
        endif
c-----
c Load interpolation arrays
        do 1000 i = i1, i2
                j = iload + (i-i1)
                t(i) = xtab (j)
                u(i) = ytab (j)
1000    continue
c
c-----
c Get interpolated value
2000    call codim1 (x, t, u, i1, i2, y)
        return
        end
c
c--------------------------------------------------------------
c
        subroutine srch (itab, x, xtab, ntab, index, ierr)
c
c Search table of x-values to bracket the desired interpolant, x
c
c The returned search index will be:
c       2 - if extrapolation below the table is required
c    ntab -                  above
c   index - points to value just above x in the table if bounded.
c
c The index is saved as a starting point for subsequent entries
c in an array indexed through 'itab' which serves to label the
c set of saved search indices. Itab may be between 1 and 20.
c
c       itab - The table identifier (1-20)
c          x - The value for which an index is desired
c       xtab - The table containing the x-values (array of length ntab)
c       ntab - number of elements in the table
c      index - returned index into the table (points just above x)
c       ierr - 0 for ok, -1 for extrapolation
c
c       Modified to remove entry points.  3/20/86 Valdes
c
        integer ntab, index, init
        double precision xtab(ntab), x
c
c common for subroutines intrp0 and intrpi
c
        common /insvcm/ insave(20)
c
c initialize
        data init/0/
c
c-----
c Initialize
        if (init.eq.0) then
            do 1110 i = 1, 20
1110            insave(i) = 0
            init = 1
        endif
c
c Determine direction of table, ascending or descending
        idir = sign (1.0d0, xtab(ntab) - xtab(1))
c
c-----
c Reset error flag
        ierr = 0
c
c-----
c Check for previous insaved index
        last = insave(itab)
        if (last .eq. 0 .or. last .gt. ntab) then
c
c-----
c no previous entry
                isrch = 1
c check for extrapolation
                if ((x-xtab(   1)) * idir .lt. 0.0d0) go to 2000
                if ((x-xtab(ntab)) * idir .gt. 0.0d0) go to 2100
        else
c
c-----
c previous entry left a valid index
                isrch = last
c
c check for still wihin bounds - difference from above should be opposite
c                                sign of difference from below
c
                if ((xtab(last)-x) * (xtab(last-1)-x) .lt. 0.0d0) then
                        index = last
                        return
                endif
        endif
c
c -----
c Begin searching - first determine direction
c
c This change made because x = xtab(1) was considered extrapolation.
c       if ((x - xtab(isrch)) * idir .gt. 0.0d0) then
        if ((x - xtab(isrch)) * idir .ge. 0.0d0) then
c forward
                do 1100 i = isrch+1, ntab
                        if ((x-xtab(i)) * idir .gt. 0.0d0) go to 1100
                        go to 1500
1100            continue
c fall thru implies extrapolation required at high end
                go to 2100
        else
c
c-----
c negative direction search
                do 1200 i = isrch-1,1,-1
                        if ((x-xtab(i)) * idir .lt. 0.0d0) go to 1200
                        go to 1400
1200            continue
c fall through implies extrapolation at low end
                go to 2000
        endif
c
c-----
c point has been bounded
1400    index = i + 1
        go to 3000
1500    index = i
        go to 3000
c
c-----
c extrapolations
2000    index = 2
        ierr = -1
        go to 3000
2100    index = ntab
        ierr = -1
        go to 3000
c
c-----
c insave index
3000    insave(itab) = index
        end
c
c------
c Subroutine to reset saved index
        subroutine intrp0 (itab)
        integer itab
        common /insvcm/ insave(20)
c
        insave(itab) = 0
        end
c
c-----
c Subroutine to return current index
        subroutine intrpi (itab, ind)
        integer itab, ind
        common /insvcm/ insave(20)
c
        ind = insave(itab)
        end
c
c-------------------------------------------------------------------
c
        subroutine codim1 (x, t, u, i1, i2, y)
c
c this subroutine performs an interposlation in a fashion
c not really understandable, but it works well.
c
c       x - input independent variable
c       t - array of 4 table independents surrounding x if possible
c       u - array of 4 table dependents corresponding to the t array
c
c  i1, i2 - indicators as follows:
c
c           i1 = 1, i2 = 4  :  4 pts available in t and u arrays
c           i1 = 1, i2 = 3  :  3 pts available (x near right edge of table)
c           i1 = 2, i2 = 4  :                  (x near left  edge of table)
c           i1 = 2, i2 = 3  :  2 pts available
c           i1 = 3, i3 = 3  :  1 pt  available
c
c       y - output interpolated (or extrapolated) dependent value
c
        double precision t(4), u(4), x, y
        integer i1, i2
        double precision s, v, z, a1, a2, a3, c1, c2, c3, a4, c4, c5, c6
	double precision e1, e2, p1, p2, slope1, slope2, al, bt, xe

c
c variable xk affects the extrapolation procedure. a value of -1.0
c appears to be a reliable value.
c
        data xk/-1.0d0/
c
        v = x
c the following code is extracted from an original source
c
      a2=v-t(2)
      al=a2/(t(3)-t(2))
      s=al*u(3)+(1.-al)*u(2)
      if(i1.gt.1.and.i2.lt.4)goto1530
      a3=v-t(3)
      if(i1.gt.1)goto1185
1180  a1=v-t(1)
      c1=a2/(t(1)-t(2))*a3/(t(1)-t(3))
      c2=a1/(t(2)-t(1))*a3/(t(2)-t(3))
      c3=a1/(t(3)-t(1))*a2/(t(3)-t(2))
      p1=c1*u(1)+c2*u(2)+c3*u(3)
      if(i2.lt.4)goto1400
1185  a4=v-t(4)
      c4=a3/(t(2)-t(3))*a4/(t(2)-t(4))
      c5=a2/(t(3)-t(2))*a4/(t(3)-t(4))
      c6=a2/(t(4)-t(2))*a3/(t(4)-t(3))
      p2=c4*u(2)+c5*u(3)+c6*u(4)
      if(i1.eq.1)goto1500
1200  if(xk.lt.0.)goto1230
      xe=xk
      goto1260
1230  slope1=abs((u(4)-u(3))/(t(4)-t(3)))
      slope2=abs((u(3)-u(2))/(t(3)-t(2)))
      xe=1.0d0
      if(slope1+slope2.ne.0.)xe=1.-abs(slope1-slope2)/(slope1+slope2)
1260  p1=s+xe*(p2-s)
      goto1500
1400  if(xk.lt.0.)goto1430
      xe=xk
      goto1460
1430  slope1=abs((u(2)-u(1))/(t(2)-t(1)))
      slope2=abs((u(3)-u(2))/(t(3)-t(2)))
      xe=1.0d0
      if(slope1+slope2.ne.0.)xe=1.-abs(slope1-slope2)/(slope1+slope2)
1460  p2=s+xe*(p1-s)
1500  e1=abs(p1-s)
      e2=abs(p2-s)
      if(e1+e2.gt.0.)goto1560
1530  z=s
      goto1700
1560  bt=(e1*al)/(e1*al+(1.-al)*e2)
      z=bt*p2+(1.-bt)*p1
c
1700    y = z
        return
        end
c
c----------------------------------------------------------------------
c
        subroutine lintrp (itab, xtab, ytab, ntab, x, y, ierr)
c
c Linear interpolator with last index save
c
c Arguments are identical to INTRP, and uses the same index search
c scheme so that values for ITAB should not clash with calls
c to INTRP and LINTRP.
c
        double precision xtab(ntab), ytab(ntab), x , y
        integer itab, ierr
c
c----- Only 1 pt in table
        if (ntab .eq. 1) then
                y = ytab (1)
                ierr = 0
                return
        endif
c
c-----locate search index
        call srch (itab, x, xtab, ntab, index, ierr)
c
c----- index points just above x
        y = ytab(index-1) + (x - xtab(index-1)) * 
     1   (ytab(index) - ytab(index-1)) / (xtab(index) - xtab(index-1))
c
        return
        end