
      subroutine terp(xin,xat,nx,yout,yat,ny,log,spval,ier,iprnt)
      dimension xin(nx),xat(nx),yout(ny),yat(ny)
c
      ier = 0
      if (nx.le.0) then
        ier = 1
        if (iprnt.gt.0) write(6,"('>>> terp: bad nx=',i5)") nx
        return
      elseif (nx.eq.1) then
        do i=1,ny
          yout(i) = xin(1)
        enddo
        return
      endif
      do i=1,ny
        if (yat(i).lt.xat(1)) then
          k0 = 2
          k1 = 1
        elseif (yat(i).gt.xat(nx)) then
          k0 = nx-1
          k1 = nx
        else
          do ii=2,nx
            if (yat(i).ge.xat(ii-1).and.yat(i).le.xat(ii)) then
              k0 = ii-1
              k1 = ii
            endif
          enddo
        endif
c
c Linear or log interpolation:
c
        if (log.gt.0) then
          if (xin(k0).eq.spval.or.xin(k1).eq.spval.or.xin(k1).le.0.) 
     +      then
            yout(i) = spval
            if (iprnt.gt.0) then
              write(6,"('terp: cannot do exponential interp with ',
     +          ' xin(k0)=',e12.4,' xin(k1)=',e12.4,' k0,1=',2i4,
     +          ' (returning ',e12.4,')')") xin(k0),xin(k1),k0,k1,spval
            endif
          else
            exparg = (alog(xin(k1) / xin(k0)) / (xat(k1)-xat(k0))) *
     +               (yat(i)-xat(k0)) 
            yout(i) = xin(k0) * exp(exparg)
          endif
        else
          f1 = (yat(i)-xat(k0)) / (xat(k1)-xat(k0))
          yout(i) = f1 * xin(k1) + (1.-f1) * xin(k0)
        endif
      enddo
      return
      end
