calcul_normales Subroutine

   subroutine calcul_normales(tab_in, long, larg, scale_xyz, cone_angle, hori, print_mask)
   !================================================================================================
   !! Function that returns the fraction of surface nearly horizontal (normal less than 5 degrees
   !!        from a vertical line)
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                           :: long         !! *surface array length*
   integer(kind=I4), intent(in )                           :: larg         !! *surface array width*
   real   (kind=R8), intent(in )                           :: cone_angle   !! *cone angle*
   real   (kind=R8), intent(out)                           :: hori         !! *fraction of facets nearly horizontal*
   logical(kind=I4), intent(in ), optional                 :: print_mask   !! *mask output ?*
   real   (kind=R8), intent(in ), dimension(1:long,1:larg) :: tab_in       !! *surface array*
   real   (kind=R8), intent(in ), dimension(1:3)           :: scale_xyz    !! *lag along x, y and scale z*

      integer(kind=I4) :: i, j, ino, nb_no, ua
      real   (kind=R8) :: z1, z2, z3, z4
      real   (kind=R8) :: n1x, n1y, n1z, n2x, n2y, n2z, nx, ny, nz, norme
      real   (kind=R8) :: phi, angle_moy, ech_x, ech_y, ech_z, ech_r, d13, d24, d42

      character(len=16) :: str

      real(kind=R8), allocatable, dimension(:)     :: angle
      real(kind=R8), allocatable, dimension(:,:,:) :: vec

      nb_no = long*larg !  nombre de noeuds auquels on associera une normale

      allocate( angle(1:nb_no) )                ! recevra les angles à la verticale pour chaque noeud
      allocate(   vec(1:long, 1:larg, 1:3) )    ! recevra pour chaque noeud la contribution des facettes triangulaires connectées

      ech_x = scale_xyz(1) !SCALE_IMG%dx * unit2IUf(SCALE_IMG%dx_unit)
      ech_y = scale_xyz(2) !SCALE_IMG%dy * unit2IUf(SCALE_IMG%dy_unit)
      ech_z = scale_xyz(3)

      ech_r = ech_x/ech_y

      d13 = (2./3.)*sqrt(  ech_x**2    +  ech_y**2    )/2
      d42 = (2./3.)*sqrt( (ech_x/2)**2 +  ech_y**2    )
      d24 = (2./3.)*sqrt(  ech_x**2    + (ech_y/2)**2 )

      ! Raisonnement sur chaque carré du domaine
      vec(1:long, 1:larg, 1:3) = 0
      do j = 1, larg -1
      do i = 1, long -1

         z1 = tab_in(i   , j   ) * ech_z
         z4 = tab_in(i   , j +1) * ech_z
         z2 = tab_in(i +1, j   ) * ech_z
         z3 = tab_in(i +1, j +1) * ech_z

         ! Triangle 1
         n1x = (z1 -z2)/ech_x
         n1y = (z1 -z4)/ech_y
         n1z = 1

         norme = sqrt( n1x**2 +n1y**2 +n1z**2)
         n1x = n1x/norme
         n1y = n1y/norme
         n1z = n1z/norme


         ! Triangle 2
         n2x = (z4 -z3)/ech_x
         n2y = (z2 -z3)/ech_y
         n2z = 1

         norme = sqrt( n2x**2 +n2y**2 +n2z**2)
         n2x = n2x/norme
         n2y = n2y/norme
         n2z = n2z/norme

         ! triangle 1
         vec(i +0, j +0, 1:3) = vec(i +0, j +0, 1:3) +[n1x, n1y, n1z]/d13 ! node 1
         vec(i +0, j +1, 1:3) = vec(i +0, j +1, 1:3) +[n1x, n1y, n1z]/d24 ! node 2
         vec(i +1, j +0, 1:3) = vec(i +1, j +0, 1:3) +[n1x, n1y, n1z]/d42 ! node 4

         ! triangle 2
         vec(i +1, j +1, 1:3) = vec(i +1, j +1, 1:3) +[n2x, n2y, n2z]/d13 ! node 3
         vec(i +0, j +1, 1:3) = vec(i +0, j +1, 1:3) +[n2x, n2y, n2z]/d42 ! node 2
         vec(i +1, j +0, 1:3) = vec(i +1, j +0, 1:3) +[n2x, n2y, n2z]/d24 ! node 4


      enddo
      enddo

      ino = 0
      do j = 1, larg
      do i = 1, long

         nx = vec(i, j, 1)
         ny = vec(i, j, 2)
         nz = vec(i, j, 3)

         norme = sqrt( nx**2 +ny**2 +nz**2)

         nx = nx/norme
         ny = ny/norme
         nz = nz/norme

         phi = acos( nz )*180./PI_R8
         ino = ino +1
         angle(ino) = phi

      enddo
      enddo

      ! angle moyen de la surface
      angle_moy = 0!sum( angle(1:nb_no) )/nb_no

      ! tout ce qui dépasse l'angle moyen (à cone_angle° près) est masqué à 0
      where( abs(angle - angle_moy) < cone_angle )
         angle = 1.
      elsewhere
         angle = 0.
      endwhere
      ! pourcentage de facettes quasi horizontales
      hori = 100*sum( angle(1:nb_no) )/nb_no

      if ( present(print_mask) ) then

         str = repeat( ' ', len(str) )
         write( str, '(I12.12)' ) ino
         call get_unit(ua)
         open( unit = ua, file = "out/mask_angle.txt")
            write( ua, '('//trim(str)//'I2.1)' ) ( int( angle(i) ), i = 1, long * larg )
         close( ua )

      endif

      deallocate( angle, vec )

   return
   endsubroutine calcul_normales