roll_smooth Subroutine

   subroutine roll_smooth(tabin, tabou, long, larg, scale_xyz, sgn, ray, omp, nb_div)
   !================================================================================================
   !< @note
   !<
   !< A ball of radius "ray" rolls on / below the surface, hence defining a closing or an opening enveloppe.
   !<
   !< @endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                            :: long        !! *2D array length*
   integer(kind=I4), intent(in )                            :: larg        !! *2D array width*
   integer(kind=I4), intent(in )                            :: sgn         !! *+ 1: dilation, -1:erosion*
   integer(kind=I4), intent(in )                            :: nb_div      !! *number of macro elements along an axis*
   logical(kind=I4), intent(in )                            :: omp         !! *if multithreading*
   real   (kind=R8), intent(in )                            :: ray         !! *roll radius*
   real   (kind=R8), intent(in ), dimension(1:3)            :: scale_xyz   !! *lag along x, y and scale z*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tabin       !! *2D array in*
   real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: tabou       !! *2D array out*

      integer(kind=I4) :: i, j, ii, jj, nb_th
      integer(kind=I4) :: hw

      integer(kind=I4) :: ik, jk, idiv, jdiv, ista, iend, jsta, jend

      real(kind=R8) :: h1, h2, ht, tmp, delta_h_max
      real(kind=R8) :: ech_x, ech_y

      real(kind=R8), allocatable, dimension(:,:) :: elem, tab_tmp

      integer(kind=I4), allocatable, dimension(:,:) :: thw

      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)

      h1 = minval( tabin(1:long, 1:larg) )
      h2 = maxval( tabin(1:long, 1:larg) )

      delta_h_max = h2 - h1

      ht = min( -abs(h1), -abs(h2) )
      h2 = max( +abs(h1), +abs(h2) )
      h1 = ht

      ! the normal width of the ball is : hw = int( ray / ech_x )
      ! However the ball curvature makes the ball height sometimes higher than the surfaces heights
      hw = int( sqrt( 2*delta_h_max*ray - delta_h_max**2 ) / ech_x ) !+ 1

      allocate( elem(   -hw    :hw,        -hw    :hw       ) )
      allocate( tab_tmp(-hw + 1:hw + long, -hw + 1:hw + larg) ) ! surface is extended

      tab_tmp(1:long, 1:larg) = tabin(1:long, 1:larg)           ! original surface

      do i = 1, long

         tab_tmp(i,  -hw + 1:         0) = tab_tmp(i,    1)
         tab_tmp(i, larg + 1: larg + hw) = tab_tmp(i, larg)

      enddo

      do j = -hw + 1, larg + hw

         tab_tmp( -hw + 1:        0, j) = tab_tmp(1,    j)
         tab_tmp(long + 1:long + hw, j) = tab_tmp(long, j)

      enddo

      nb_th = 1
      if (omp) then
         nb_th = omp_get_num_procs()
      endif

      allocate( thw(1:nb_div, 1:nb_div) )    ! number of macro squares on which the ball active width is determined

      idiv = long / nb_div
      jdiv = larg / nb_div

      !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th) IF (omp)
      !$OMP DO  SCHEDULE (STATIC, max(nb_div/nb_th, 1)) PRIVATE(jk, ik, delta_h_max, ista, iend, jsta, jend)

      do jk = 1, nb_div

      if (jk==nb_div) then ; jend = larg ; else ; jend = jk * jdiv ; endif ; jsta = 1 + jdiv * (jk-1)

      do ik = 1, nb_div

      if (ik==nb_div) then ; iend = long ; else ; iend = ik * idiv ; endif ; ista = 1 + idiv * (ik-1)

            delta_h_max = maxval( tab_tmp(-hw + ista:hw + iend, -hw + jsta:hw + jend) ) - & !
                          minval( tab_tmp(-hw + ista:hw + iend, -hw + jsta:hw + jend) )

            ! beyond the present width, the ball height is greater than the surface height
            thw(ik, jk) = int( sqrt( 2*delta_h_max*ray - delta_h_max**2 ) / ech_x )

      enddo

      enddo

      !$OMP END DO
      !$OMP END PARALLEL

      if (sgn == +1) then

         do jj = -hw, hw
         do ii = -hw, hw
            tmp = ray**2 - (ii * ech_x)**2 - (jj * ech_y)**2
            if ( tmp < 0. ) then
               elem(ii, jj) = ray               + 1.1 * abs(h2)
            else
               ! the ball location is a little above the surface
               elem(ii, jj) = ray - sqrt( tmp )
            endif
         enddo
         enddo

      else

         do jj = -hw, hw
         do ii = -hw, hw
            tmp = ray**2 - (ii * ech_x)**2 - (jj * ech_y)**2
            if ( tmp < 0. ) then
               elem(ii, jj) = -ray               - 1.1 * abs(h1)
            else
               ! the ball location is a little below the surface
               elem(ii, jj) = -ray + sqrt( tmp )
            endif
         enddo
         enddo

      endif

      !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th) IF (omp)
      !$OMP DO SCHEDULE (STATIC,larg/nb_th) PRIVATE(i, j, jk, ik, hw)

      do j = 1, larg ; jk = 1 + j/(jdiv + 1)
      do i = 1, long ; ik = 1 + i/(idiv + 1)

         hw = thw(ik, jk)

         tabou(i, j) = - minval( sgn * ( elem(-hw    :hw    , -hw    :hw    ) - & !
                                      tab_tmp(-hw + i:hw + i, -hw + j:hw + j) ) )

      enddo
      enddo

      !$OMP END DO
      !$OMP END PARALLEL

      if (sgn == +1) then

         tabou(1:long, 1:larg) = + tabou(1:long, 1:larg) + ray

      else

         tabou(1:long, 1:larg) = - tabou(1:long, 1:larg) - ray

      endif

      deallocate( elem, tab_tmp, thw )

   return
   endsubroutine roll_smooth