subroutine interpCoeff16(x,y,xp,yp,intco)

      implicit none

      real,dimension(4,4)::a,x,y,intco
      real::xp,yp
      real,dimension(4)::b,ypp
      integer::i,j,k

      do i=1,4
      do j=1,4
         a(i,j)=1.
	 do k=1,4
	   if(k.ne.i)then
	     a(i,j)=a(i,j)*(xp-x(k,j))/(x(i,j)-x(k,j))
           endif
         enddo
      enddo
      enddo

      do j=1,4
        ypp(j)=0.
	do i=1,4
	  ypp(j)=ypp(j)+a(i,j)*y(i,j)
        enddo
      enddo

      do j=1,4
        b(j)=1.
	do k=1,4
	  if(k.ne.j)then
	    b(j)=b(j)*(yp-ypp(k))/(ypp(j)-ypp(k))
          endif
        enddo
      enddo

      do i=1,4
      do j=1,4
        intco(i,j)=a(i,j)*b(j)
      enddo
      enddo

      end subroutine interpCoeff16