compute_face_area_vectors Subroutine

        subroutine compute_face_area_vectors(Ifaces, Jfaces, Kfaces, nodes)
            !< Compute face area vectors
            !<
            !< The face area vectors denote the face area both in 
            !< magnitude and direction. They are placed in the face 
            !< normal variables for further calculation.
            !<
            !< The face normal is given by d1 x d2, where d1 and d2 are
            !< the diagonals of a face
            !-----------------------------------------------------------
            
            implicit none
            type(facetype), dimension(-2:imx+3,-2:jmx+2,-2:kmx+2), intent(inout) :: Ifaces
            !< Store face quantites for I faces: normal and area
            type(facetype), dimension(-2:imx+2,-2:jmx+3,-2:kmx+2), intent(inout) :: Jfaces
            !< Store face quantites for J faces: normal and area
            type(facetype), dimension(-2:imx+2,-2:jmx+2,-2:kmx+3), intent(inout) :: Kfaces
            !< Store face quantites for K faces: normal and area
            type(nodetype), dimension(-2:imx+3,-2:jmx+3,-2:kmx+3), intent(in)  :: nodes
            !< Grid points

    
            real(wp) :: d1x, d2x, d1y, d2y, d1z, d2z
            integer :: i, j, k

            do k = -2, kmx+2
             do j = -2, jmx+2
              do i = -2, imx+3
                d1x = nodes(i, j+1, k+1)%x - nodes(i, j, k)%x
                d1y = nodes(i, j+1, k+1)%y - nodes(i, j, k)%y
                d1z = nodes(i, j+1, k+1)%z - nodes(i, j, k)%z
                d2x = nodes(i, j, k+1)%x - nodes(i, j+1, k)%x
                d2y = nodes(i, j, k+1)%y - nodes(i, j+1, k)%y
                d2z = nodes(i, j, k+1)%z - nodes(i, j+1, k)%z
                Ifaces(i, j, k)%nx = 0.5 * (d1y*d2z - d1z*d2y)
                Ifaces(i, j, k)%ny = 0.5 * (d1z*d2x - d1x*d2z)
                Ifaces(i, j, k)%nz = 0.5 * (d1x*d2y - d1y*d2x)
               end do
              end do
             end do


            do k = -2, kmx+2
             do j = -2, jmx+3
              do i = -2, imx+2
               d1x = nodes(i+1, j, k+1)%x - nodes(i, j, k)%x
               d1y = nodes(i+1, j, k+1)%y - nodes(i, j, k)%y
               d1z = nodes(i+1, j, k+1)%z - nodes(i, j, k)%z
               d2x = nodes(i+1, j, k)%x - nodes(i, j, k+1)%x
               d2y = nodes(i+1, j, k)%y - nodes(i, j, k+1)%y
               d2z = nodes(i+1, j, k)%z - nodes(i, j, k+1)%z
            
               Jfaces(i, j, k)%nx = 0.5 * (d1y*d2z - d1z*d2y)
               Jfaces(i, j, k)%ny = 0.5 * (d1z*d2x - d1x*d2z)
               Jfaces(i, j, k)%nz = 0.5 * (d1x*d2y - d1y*d2x)
              end do
             end do
            end do

            
            do k = -2, kmx+3
             do j = -2, jmx+2
              do i = -2, imx+2
               d1x = nodes(i+1, j+1, k)%x - nodes(i, j, k)%x
               d1y = nodes(i+1, j+1, k)%y - nodes(i, j, k)%y
               d1z = nodes(i+1, j+1, k)%z - nodes(i, j, k)%z
               d2x = nodes(i, j+1, k)%x - nodes(i+1, j, k)%x
               d2y = nodes(i, j+1, k)%y - nodes(i+1, j, k)%y
               d2z = nodes(i, j+1, k)%z - nodes(i+1, j, k)%z
            
               Kfaces(i, j, k)%nx = 0.5 * (d1y*d2z - d1z*d2y)
               Kfaces(i, j, k)%ny = 0.5 * (d1z*d2x - d1x*d2z)
               Kfaces(i, j, k)%nz = 0.5 * (d1x*d2y - d1y*d2x)
              end do
             end do
            end do


        end subroutine compute_face_area_vectors