       subroutine blii2d(x,y,f,xi,yi,fi,fxi,fyi,fxyi)
	  
c***********************************************************************
c      Barycentric Linear Interpolation 2D 
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Gambelas, 15/06/2011 at 14:00
c***********************************************************************

       real*8 f(2,2)
       real*8 x(2),y(2)
       real*8 px(2),py(2)
       real*8 fi,fxi,fyi,fxyi
       real*8 xi,x1,x2
       real*8 yi,y1,y2
       real*8 a12,a21,a22
       real*8 pxy

       x1 = x(1)
       x2 = x(2)
       
       y1 = y(1)
       y2 = y(2)

       pxy = (x2-x1)*(y2-y1)
       
       a12 = (f(1,1)-f(1,2))/pxy
       a21 = (f(1,1)-f(2,1))/pxy
       a22 = (f(2,2)-f(1,1))/pxy
       
       px(1) = xi - x1
       px(2) = xi - x2
       
       py(1) = yi - y1
       py(2) = yi - y2
       
       fi = f(1,1)+a12*px(1)*py(2)+a21*px(2)*py(1)+a22*px(1)*py(1)
     
       fxi  = a12*py(2) + (a21+a22)*py(1)
       fyi  = (a12+a22)*px(1) + a21*px(2)
       fxyi = a12 + a22 + a21
       
c***********************************************************************
c     Back to main:
c***********************************************************************
       
       return 
       
       end

