 SUBROUTINE triangle(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn,surf,xm,ym,zm)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    GEO/triangle  *
!
!    NOM
! triangle(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn,surf,xm,ym,zm)
!    DESCRIPTION  
!*
! Calcul du vecteur normal, de la surface et du barycentre d un triangle
!
! Calculate the normal vector, area and barycentre of a triangle
!*
!
!    ENTREES / INPUT
!    x1,y1,z1,x2,y2,z2,x3,y3,z3 : 3 sommets du triangle / 3 vertices of the triangle
!
!    SORTIES / OUTPUT
!    xn,yn,zn,surf,xm,ym,zm: vecteur normal, surface et barycentre du triangle 
!                          / normal vector, area and barycentre of the triangle
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  IMPLICIT NONE
!.2-----  Declaration
  real(kind=kind(0.d0)), intent(in)   :: x1,y1,z1,x2,y2,z2,x3,y3,z3
!! coordonnees des sommets du triangle / coordinates of the vertices of the triangle
  real(kind=kind(0.d0)), intent(out)   :: xn,yn,zn 
!! vecteur normal du triangle / normal vector of the triangle
  real(kind=kind(0.d0)), intent(out)   :: surf !! Surface du triangle / Surface of the triangle
  real(kind=kind(0.d0)), intent(out)   :: xm,ym,zm
!! coordonnees du barycentre du triangle / coordinates of the barycentre of the triangle

  real(kind=kind(0.d0))   :: x12,y12,z12,x13,y13,z13
!=========================== DEBUT DU CODE EXECUTABLE ==================

  x12=x2-x1
  y12=y2-y1
  z12=z2-z1
  
  x13=x3-x1
  y13=y3-y1
  z13=z3-z1
  
  ! normale et surface du triangle 
  xn=y12*z13-z12*y13
  yn=-x12*z13+z12*x13
  zn=x12*y13-y12*x13
  
  xn=-xn
  yn=-yn
  zn=-zn
  
  surf=dsqrt(xn*xn+yn*yn+zn*zn)
  
  xm=1.d0/3.d0*(x1+x2+x3)
  ym=1.d0/3.d0*(y1+y2+y3)
  zm=1.d0/3.d0*(z1+z2+z3)
            
  IF (surf.ne.0) THEN
      xn=xn/surf
      yn=yn/surf
      zn=zn/surf
  ELSE
      xn=-1
      yn=0
      zn=0
  END IF
  
  surf=surf*0.5d0
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE TRIANGLE
!DOCFIN