	function polint (xa, ya, nx, ny, j, x)
	dimension xa(nx), ya(nx, ny)

	ns= 0
        do i= 1, nx-1 
    	  if ((x .lt. xa(i+1)) .and. (x .gt. xa(i))) ns= i
 	enddo

	if (ns .eq. 0) then
	  polint= 0.0	
	  return
	endif
        
	t1= ((ya(ns+1,j)-ya(ns,j))/(xa(ns+1)-xa(ns)))*x
	t2= ((ya(ns+1,j) - ya(ns,j))/(xa(ns+1)-xa(ns)))*xa(ns)
	y= t1 + (ya(ns,j) - t2)

	polint= y

	end
