flood Subroutine

   subroutine flood(masque, taille, nx, ny, niv)
   !================================================================================================
   !< @note
   !<
   !< Perform some kind of flood fill or connected component labeling on a grid (masque),
   !< starting from an initial '1' element found and spreading out to adjacent '1' elements,
   !< updating them to a specified value or zero if no value (niv) is specified.
   !<
   !< @endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in   )                        :: nx, ny          !! *Input dimensions*
   integer(kind=I4), intent(in   ), optional              :: niv             !! *Optional input level*
   integer(kind=I4), intent(inout), dimension(1:nx, 1:ny) :: masque          !! *Input/output matrix*
   integer(kind=I4), intent(  out)                        :: taille          !! *Output scalar*

   integer(kind=I4) :: i, j, p, q, k, ind, ival
   integer(kind=I4) :: nb_traite

   ! Arrays to hold position indices and a tracking matrix
   integer(kind=I4), dimension(nx*ny)  :: liste_x, liste_y
   integer(kind=I4), dimension(nx, ny) :: deja_fait

      ! Initialize arrays to zero
      liste_x = 0
      liste_y = 0
      deja_fait = 0

      ! Initial position
      i = 0
      j = 0

      ! Find the first occurrence of the element "1" in masque
o:    do i = 1, nx
         do j = 1, ny
            if ( masque(i, j) == 1 ) exit o
         enddo
      enddo o

      ! If no element is found, set output size to zero and exit subroutine
      if ( i == nx + 1 .and. j == ny + 1 ) then
         taille = 0
         return
      endif

      ! Initialize processing of found element
      ind = 1
      liste_x(ind) = i
      liste_y(ind) = j

      ! Mark element as processed
      deja_fait(i, j) = 1

      ! Start processing of elements
      nb_traite = 1

      ! Explore neighbors in 8 directions (up, down, left, right and combinations)
      do
         i = liste_x(nb_traite)
         j = liste_y(nb_traite)

         ! Break the loop if outside bounds
         if ( i*j == 0 ) exit

         ! Check south direction
         k = 1
         do

            p = i
            q = j - k

            if ( q < 1 ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check south-west direction
         k = 1
         do

            p = i - k
            q = j - k

            if ( p < 1 ) exit
            if ( q < 1 ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check south-east direction
         k = 1
         do

            p = i + k
            q = j - k

            if ( p > nx ) exit
            if ( q < 1  ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check north direction
         k = 1
         do

            p = i
            q = j + k

            if ( q > nx  ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check north-west direction
         k = 1
         do

            p = i - k
            q = j + k

            if ( p < 1  ) exit
            if ( q > ny ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check north-east direction
         k = 1
         do

            p = i + k
            q = j + k

            if ( p > nx  ) exit
            if ( q > ny ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check west direction
         k = 1
         do

            p = i - k
            q = j

            if ( p < 1  ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Check east direction
         k = 1
         do

            p = i + k
            q = j

            if ( p > nx  ) exit
            if (    masque(p, q) == 0 ) exit
            if ( deja_fait(p, q) == 0 ) then
               ind = ind + 1
               liste_x(ind) = p
               liste_y(ind) = q
               deja_fait(p, q) = 1
            endif

            k = k + 1

         enddo

         ! Increment processed count
         nb_traite = nb_traite + 1

      enddo

      ! Set output size to number of processed elements
      taille = ind

      ! Update masque values based on presence of niv and processed nodes
      ival = 0
      if ( present( niv ) ) ival = niv
      where( deja_fait == 1) masque = ival ! Fill the connected region in 'masque'

   return
   endsubroutine flood