compute_faces_metrics Subroutine

   elemental subroutine compute_faces_metrics(self)
   !< Compute block faces metrics.
   class(block_object), intent(inout) :: self         !< Block.
   type(vector)                       :: triplet(1:4) !< Dummy vectors.
   real(R8P)                          :: signi        !< Sign of direction of normals along I coordinate.
   real(R8P)                          :: signj        !< Sign of direction of normals along J coordinate.
   real(R8P)                          :: signk        !< Sign of direction of normals along K coordinate.
   integer(I4P)                       :: i            !< Counter.
   integer(I4P)                       :: j            !< Counter.
   integer(I4P)                       :: k            !< Counter.

   associate(node=>self%node, ni=>self%signature%ni, nj=>self%signature%nj, nk=>self%signature%nk)
      i = max(1, ni)
      j = max(1, nj)
      k = max(1, nk)
      triplet(1) = node(i, j,   k)%vertex - node(i, j-1, k-1)%vertex
      triplet(2) = node(i, j-1, k)%vertex - node(i, j,   k-1)%vertex
      triplet(3) = triplet(1).cross.triplet(2)
      triplet(4) = (0.25_R8P * (node(i,   j, k  )%vertex + node(i,   j-1, k  )%vertex   + &
                                node(i,   j, k-1)%vertex + node(i,   j-1, k-1)%vertex)) - &
                   (0.25_R8P * (node(i-1, j, k  )%vertex + node(i-1, j-1, k  )%vertex   + &
                                node(i-1, j, k-1)%vertex + node(i-1, j-1, k-1)%vertex))
      signi = sign(1._R8P, (triplet(3).dot.triplet(4)))

      triplet(1) = node(i, j, k  )%vertex - node(i-1, j, k-1)%vertex
      triplet(2) = node(i, j, k-1)%vertex - node(i-1, j, k  )%vertex
      triplet(3) = triplet(1).cross.triplet(2)
      triplet(4) = (0.25_R8P * (node(i, j,   k  )%vertex + node(i-1, j,   k  )%vertex   + &
                                node(i, j,   k-1)%vertex + node(i-1, j,   k-1)%vertex)) - &
                   (0.25_R8P * (node(i, j-1, k  )%vertex + node(i-1, j-1, k  )%vertex   + &
                                node(i, j-1, k-1)%vertex + node(i-1, j-1, k-1)%vertex))
      signj = sign(1._R8P, (triplet(3).dot.triplet(4)))

      triplet(1) = node(i,   j, k)%vertex - node(i-1, j-1, k)%vertex
      triplet(2) = node(i-1, j, k)%vertex - node(i,   j-1, k)%vertex
      triplet(3) = triplet(1).cross.triplet(2)
      triplet(4) = (0.25_R8P * (node(i, j, k    )%vertex + node(i-1, j, k    )%vertex   + &
                                node(i, j-1, k  )%vertex + node(i-1, j-1, k  )%vertex)) - &
                   (0.25_R8P * (node(i, j, k-1  )%vertex + node(i-1, j, k-1  )%vertex   + &
                                node(i, j-1, k-1)%vertex + node(i-1, j-1, k-1)%vertex))
      signk = sign(1._R8P, (triplet(3).dot.triplet(4)))

      do k=1, nk
         do j=1, nj
            do i=0, ni
               call self%face_i(i, j, k)%compute_metrics(pt1 = node(i, j-1, k-1)%vertex, &
                                                         pt2 = node(i, j  , k-1)%vertex, &
                                                         pt3 = node(i, j  , k  )%vertex, &
                                                         pt4 = node(i, j-1, k  )%vertex, signd=signi)
            enddo
         enddo
      enddo

      do k=1, nk
         do j=0, nj
            do i=1, ni
               call self%face_j(i, j, k)%compute_metrics(pt1 = node(i-1, j, k-1)%vertex, &
                                                         pt2 = node(i-1, j, k  )%vertex, &
                                                         pt3 = node(i  , j, k  )%vertex, &
                                                         pt4 = node(i  , j, k-1)%vertex, signd=signj)
            enddo
         enddo
      enddo

      do k=0, nk
         do j=1, nj
            do i=1, ni
               call self%face_k(i, j, k)%compute_metrics(pt1 = node(i-1, j-1, k)%vertex, &
                                                         pt2 = node(i  , j-1, k)%vertex, &
                                                         pt3 = node(i  , j  , k)%vertex, &
                                                         pt4 = node(i-1, j  , k)%vertex, signd=signk)
            enddo
         enddo
      enddo
   endassociate
   endsubroutine compute_faces_metrics