
implicit none

 real :: x(3), y(3), z(3)
 real :: xx, yy, zz
 integer :: i,j


  x(1) = 0; y(1) = 0; z(1) =  sin(x(1)+y(1))
  x(2) = 1; y(2) = 0; z(2) =  sin(x(2)+y(2))
  x(3) = 0; y(3) = 1; z(3) =  sin(x(3)+y(3))

   
 do i=1, 3
  do j=1, 3
   xx = 0.5*(i-1); yy=0.5*(j-1)
    call interpol2d( x, y, z, xx, yy, zz )

    write(*,*) xx, yy, zz
  end do
 end do
end


subroutine interpol2d(x, y, z, xx, yy, zz )
implicit none
 real :: x(3), y(3), z(3)
 real :: xx, yy, zz

 real ::  u1, u2, u3, v1, v2, v3
 real ::  nx, ny, nz, d

  u1 = x(2)-x(1); u2=y(2)-y(1); u3=z(2)-z(1)
  v1 = x(3)-x(1); v2=y(3)-y(1); v3=z(3)-z(1)

 call outer_product(  u1, u2, u3, v1, v2, v3, nx, ny, nz )
 call innner_product( nx, ny, nz, x(1), y(1), z(1), d )
!    eq :   nx * x + ny* y + nz * z =  d

 if( nz /= 0.0 ) then
  zz =  - (nx * xx + ny * yy - d)/nz
 else
  zz = 0
  write(0,*) ' Can not interporate!!! '
  stop
 end if

end



subroutine outer_product( u1, u2, u3, v1, v2, v3, w1, w2, w3)
!  w = u x v
  implicit none
  real :: u1, u2, u3, v1, v2, v3, w1, w2, w3

   w1 = u2*v3-u3*v2
   w2 = u3*v1-u1*v3
   w3 = u1*v2-v1*u2

end

subroutine innner_product( u1, u2, u3, v1, v2, v3, s )
! s = (u, v)
  implicit none
  real :: u1, u2, u3, v1, v2, v3, s

  s = u1*v1+u2*v2+u3*v3
end
