read_faces Subroutine

   subroutine read_faces(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, i, id
      integer :: owner, neighbor, patch
      real(rk) :: nx, ny, nz, area, cx, cy, cz

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('mesh_io', 'could not open '//trim(filename))

      read(unit_id, *, iostat=ios) m%nfaces
      if (ios /= 0 .or. m%nfaces <= 0) call fatal_error('mesh_io', 'invalid faces header')

      allocate(m%faces(m%nfaces))
      do i = 1, m%nfaces
         read(unit_id, *, iostat=ios) id, owner, neighbor, patch, nx, ny, nz, area, cx, cy, cz
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading face')
         if (id < 1 .or. id > m%nfaces) call fatal_error('mesh_io', 'face id out of range')
         if (owner < 1 .or. owner > m%ncells) call fatal_error('mesh_io', 'face owner out of range')
         if (neighbor < 0 .or. neighbor > m%ncells) call fatal_error('mesh_io', 'face neighbor out of range')
         if (area <= 0.0_rk) call fatal_error('mesh_io', 'face area must be positive')
         m%faces(id)%id = id
         m%faces(id)%owner = owner
         m%faces(id)%neighbor = neighbor
         m%faces(id)%patch = patch
         m%faces(id)%normal = [nx, ny, nz]
         m%faces(id)%area = area
         m%faces(id)%center = [cx, cy, cz]
      end do

      close(unit_id)
   end subroutine read_faces