indice_fractal Subroutine

   subroutine indice_fractal(tab_in, long, larg, indf)
   !================================================================================================
   !! Function that returns the fractal dimension with the box counting method
   implicit none
   integer(kind=I4), intent(in )                            :: long   !! *surface array length*
   integer(kind=I4), intent(in )                            :: larg   !! *surface array width*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in !! *surface array*
   real   (kind=R8), intent(out), dimension(3)              :: indf   !! *result: ordinate at origin, slope, R2*

      integer(kind=I4) :: i, j, k, ib, nbhmax, nni, nb, dec_i, dec_j, ri, rj
      real   (kind=R8) :: lboite, hhmax, hhmin, ddh, t1, t2, t3, t4, hmax, hmin
      real   (kind=R8) :: yibarr, yi_m_yichap, yi_m_yibarr, r2adj

      real(kind=R8), dimension(1:2) :: var

      real   (kind=R8), allocatable, dimension(:,:) :: tab, Jf
      real   (kind=R8), allocatable, dimension(:)   :: tab_nni

      integer(kind=I4), allocatable, dimension(:,:) :: pas_i, pas_j
      integer(kind=I4), allocatable, dimension(:)   :: nnb

      type(moment_stat)                             :: mom

      call calc_moments(   tab = reshape( tab_in(1:long, 1:larg), [long * larg] ),  &  !
                            mx = mom,                                               &  !
                        nb_mom = 2 )                                                   !

      allocate(tab(1:long, 1:larg))
      tab(1:long, 1:larg) = (tab_in(1:long, 1:larg) -mom%mu)/mom%si

      hhmin = minval(tab(1:long, 1:larg))     ! hauteur min de la surface
      tab(1:long, 1:larg) = tab(1:long, 1:larg) +hhmin +1._R8

      hhmax = maxval(tab(1:long, 1:larg))     ! hauteur max de la surface
      hhmin = minval(tab(1:long, 1:larg))     ! hauteur min de la surface
      ddh   = hhmax -hhmin                    ! amplitude

      nbhmax = nint( log( UN*min(long-1,larg-1) )/log(2.) )      ! nbre max de comptages, attention il y a n-1 intervalles pour n points
      if (2**nbhmax > min(long-1, larg-1)) nbhmax = nbhmax -1

      allocate(      Jf(1:nbhmax, 1:2) )
      allocate( tab_nni(1:nbhmax) )

      nb = 2**nbhmax ! nombre de boîtes dans une direction
      allocate( pas_i(1:nb, 1:nbhmax), pas_j(1:nb, 1:nbhmax) )
      allocate(   nnb(1:nb) )

      pas_i = 0
      pas_j = 0

      pas_i(1:nb, nbhmax) = (long -1)/nb ! longueur moyenne d'UN pas selon x (le pas n'est pas forcément constant : nbpts /= 2**n +1, par exemple)
      pas_j(1:nb, nbhmax) = (larg -1)/nb ! id selon y

      ri = mod(long -1, nb) ! si la division au-dessus ne tombe pas juste, c'est le reste selon x
      rj = mod(larg -1, nb) ! id selon y

      if ( ri>0 ) then ! s'il y a un résidu, répartition régulière de ce résidu sur le découpage le plus fin
         do i = 1, ri
            pas_i( (i-1)*(nb/ri) +1, nbhmax ) = pas_i( (i-1)*(nb/ri) +1, nbhmax ) +1
         enddo
      endif
      if ( rj>0 ) then
         do j = 1, rj
            pas_j( (j-1)*(nb/rj) +1, nbhmax ) = pas_j( (j-1)*(nb/rj) +1, nbhmax ) +1
         enddo
      endif

      do ib = nbhmax-1, 1, -1 ! agglomération des pas 2 à 2 pour former le pas de la boîte englobante
         do k = 1, 2**ib
            pas_i(k, ib) = pas_i(2*(k-1)+1, ib+1) +pas_i(2*k, ib+1)
            pas_j(k, ib) = pas_j(2*(k-1)+1, ib+1) +pas_j(2*k, ib+1)
         enddo
      enddo

      do ib = 1, nbhmax               ! niveau de découpage
         nb = 2**ib                   ! nombre de boîtes dans une direction pour ce découpage
         lboite = (ddh -2*EPS_R8)/nb  ! taille z de la boîte

         nni   = 0
         dec_i = 1
         do i = 1, nb ! numéro i de la colonne de boîtes

            dec_j = 1
            do j = 1, nb ! numéro j de la colonne de boîtes
               ! on considère le plan résultant de la coupe de la colonne par la surface
               t1 = tab(dec_i              , dec_j              )
               t2 = tab(dec_i +pas_i(i, ib), dec_j              )
               t3 = tab(dec_i +pas_i(i, ib), dec_j +pas_j(j, ib))
               t4 = tab(dec_i              , dec_j +pas_j(j, ib))
               ! ce plan traverse plusieurs boîtes, il suffit de les compter
               hmax  = max(t1, t2, t3, t4)
               hmin  = min(t1, t2, t3, t4)
               nni   = nni + ceiling( hmax/lboite ) &
                           -   floor( hmin/lboite )

               dec_j = dec_j +pas_j(j, ib)
            enddo

            dec_i = dec_i +pas_i(i, ib)
         enddo

         tab_nni(ib) = log(UN*nni)
         Jf(ib, 1)   = UN
         Jf(ib, 2)   = log(UN/lboite)

      enddo

      call moindres_carres_lineaire( nb_var = 2,                   & !
                                     nb_pts = nbhmax,              & !
                                        hij = tab_nni(1:nbhmax),   & !
                                       beta = var(1:2),            & !
                                         Jf = Jf(1:nbhmax, 1:2) )    !

      yibarr = sum( tab_nni(1:nbhmax) )/(nbhmax)
      yi_m_yichap = 0
      yi_m_yibarr = 0
      do i = 1, nbhmax
         yi_m_yichap = yi_m_yichap +( tab_nni(i) -( Jf(i, 1)*var(1) +Jf(i, 2)*var(2) ) )**2
         yi_m_yibarr = yi_m_yibarr +( tab_nni(i) -yibarr )**2
      enddo
      r2adj = UN -(yi_m_yichap/(nbhmax -2))/(yi_m_yibarr/(nbhmax -1))

      indf = [var(2), var(1), r2adj] ! fractal index first

      deallocate( tab, tab_nni, Jf, pas_i, pas_j, nnb )

   return
   endsubroutine indice_fractal