BEGIN_PROVIDER [ integer, grid_x_num ]
&BEGIN_PROVIDER [ integer, grid_y_num ]
&BEGIN_PROVIDER [ integer, grid_z_num ]
&BEGIN_PROVIDER [ real   , grid_step  , (3) ]
&BEGIN_PROVIDER [ real   , grid_origin, (3) ]

 real, parameter :: UNDEFINED=123456789.123456789

 BEGIN_DOC  
! Number of grid points in x, y, z directions
 END_DOC

 integer :: Npoints(3)
 real :: step_size(3)
 real :: origin(3)
 real :: opposite(3)

 Npoints  (:)  = 80
 origin   (:)  = UNDEFINED
 opposite (:)  = UNDEFINED
 step_size(:)  = UNDEFINED

 call get_grid_point_num(Npoints)
 call get_grid_origin(origin)
 call get_grid_opposite(opposite)
 call get_grid_step_size(step_size)
 
 if (origin(1) == UNDEFINED) then
   integer :: i,l
   do l=1,3
    origin(l) = nucl_coord(1,l)
    do i=2,nucl_num
     origin(l) = min(origin(l),nucl_coord(i,l))
    enddo
    origin(l) = origin(l) - 4.
   enddo
 endif

 if (opposite(1) == UNDEFINED) then
   do l=1,3
    opposite(l) = nucl_coord(1,l)
    do i=2,nucl_num
     opposite(l) = max(opposite(l),nucl_coord(i,l))
    enddo
    opposite(l) = opposite(l) + 4.
   enddo
 endif

 if (step_size(1) == UNDEFINED) then
   do l=1,3
     step_size(l) = (opposite(l) - origin(l))/float(Npoints(l))
   enddo
 endif

 do l=1,3
   grid_origin(l) = origin(l)
   grid_step(l)   = step_size(l)
 enddo
 grid_x_num = Npoints(1)
 grid_y_num = Npoints(2)
 grid_z_num = Npoints(3)

 call set_grid_point_num(Npoints)
 call set_grid_origin(origin)
 call set_grid_opposite(opposite)
 call set_grid_step_size(step_size)

END_PROVIDER

BEGIN_TEMPLATE

BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
 implicit none
 BEGIN_DOC
! $X on a grid
 END_DOC

 IRP_IF MPI
   include 'mpif.h'
 IRP_ENDIF

 integer :: ix, iy, iz
 integer :: ibegin, iend
 real, parameter :: UNDEFINED=123456789.123456789

 grid_$X(1,1,1) = UNDEFINED
 call get_grid_data_$X(grid_$X)
 if (grid_$X(1,1,1) /= UNDEFINED) then
   return
 endif

 do iz=1,grid_z_num 
  do iy=1,grid_y_num 
   do ix=1,grid_x_num 
    grid_$X(ix,iy,iz) = 0.
   enddo
  enddo
 enddo

 if (mpi_master) then
   print *, 'Computing $X' 
 endif

 integer :: icount
 icount = mpi_size
 do iz=1,grid_z_num
  if (mpi_master) then
    print *,  int(100*dble(iz)/dble(grid_z_num)), '%'
  endif
  point(3) = grid_origin(3)+(iz-1)*grid_step(3)
  do iy=1,grid_y_num 
   point(2) = grid_origin(2)+(iy-1)*grid_step(2)
   do ix=1,grid_x_num 
    icount -= 1
    if (icount == mpi_rank) then
     point(1) = grid_origin(1)+(ix-1)*grid_step(1)
     TOUCH point
     grid_$X(ix,iy,iz) = $X_value_p
    endif
    if (icount == 0) then
      icount = mpi_size
    endif
   enddo
  enddo
 enddo

 IRP_IF MPI
   integer :: dim, ierr
   do iz=1,grid_z_num 
    real   :: buffer(grid_x_num*grid_y_num+1)
    icount = 0
    do iy=1,grid_y_num 
     do ix=1,grid_x_num 
      buffer(icount+ix) = grid_$X(ix,iy,iz)
     enddo
     icount = icount + grid_x_num
    enddo
    dim = grid_x_num * grid_y_num
    call MPI_REDUCE(buffer,grid_$X(1,1,iz),dim,mpi_real, &
       mpi_sum,0,MPI_COMM_WORLD,ierr)
    if (ierr /= MPI_SUCCESS) then
      call abrt(irp_here,'Unable to fetch buffer')
    endif
    call barrier
   enddo

 IRP_ENDIF
 call set_grid_data_$X(grid_$X)

END_PROVIDER


 BEGIN_PROVIDER [ real, grid_$X_grad, (grid_x_num,grid_y_num,grid_z_num,4) ]
&BEGIN_PROVIDER [ real, grid_$X_lapl, (grid_x_num,grid_y_num,grid_z_num) ]
 BEGIN_DOC
! Laplacian of $X on a grid
 END_DOC
 implicit none
 BEGIN_DOC
! Gradient and lapacian of $X on a grid. 4th dimension of the grad is its norm.
 END_DOC

 IRP_IF MPI
   include 'mpif.h'
 IRP_ENDIF

 integer :: ix, iy, iz, it
 integer :: ibegin, iend
 real, parameter :: UNDEFINED=123456789.123456789

 grid_$X_lapl(1,1,1) = UNDEFINED
 call get_grid_data_$X_lapl(grid_$X)
 if (grid_$X_lapl(1,1,1) /= UNDEFINED) then
   return
 endif

 do iz=1,grid_z_num 
  do iy=1,grid_y_num 
   do ix=1,grid_x_num 
    do it=1,4
     grid_$X_grad(ix,iy,iz,it) = 0.
    enddo
    grid_$X_lapl(ix,iy,iz) = 0.
   enddo
  enddo
 enddo

 if (mpi_master) then
   print *, 'Computing $X' 
 endif

 integer :: icount
 icount = mpi_size
 do iz=1,grid_z_num
  if (mpi_master) then
    print *,  int(100*dble(iz)/dble(grid_z_num)), '%'
  endif
  point(3) = grid_origin(3)+(iz-1)*grid_step(3)
  do iy=1,grid_y_num 
   point(2) = grid_origin(2)+(iy-1)*grid_step(2)
   do ix=1,grid_x_num 
    icount = icount-1
    if (icount == mpi_rank) then
     point(1) = grid_origin(1)+(ix-1)*grid_step(1)
     TOUCH point
     do it=1,3
      grid_$X_grad(ix,iy,iz,it) = $X_grad_p(it)
     enddo
     grid_$X_grad(ix,iy,iz,4) = $X_grad_p(1)**2 + $X_grad_p(2)**2 + $X_grad_p(3)**2
     grid_$X_lapl(ix,iy,iz) = $X_lapl_p
     grid_$X_grad(ix,iy,iz,4) = sqrt(grid_$X_grad(ix,iy,iz,4))
    endif
    if (icount == 0) then
      icount = mpi_size
    endif
   enddo
  enddo
 enddo

 IRP_IF MPI
  integer :: dim, ierr
  do it=1,4 
   do iz=1,grid_z_num 
    real   :: buffer(grid_x_num*grid_y_num)
    icount = 0
    do iy=1,grid_y_num 
     do ix=1,grid_x_num 
      buffer(icount+ix) = grid_$X_grad(ix,iy,iz,it)
     enddo
     icount = icount + grid_x_num
    enddo
    dim = grid_x_num * grid_y_num
    call MPI_REDUCE(buffer,grid_$X_grad(1,1,iz,it),dim,mpi_real, &
       mpi_sum,0,MPI_COMM_WORLD,ierr)
    icount = 0
    do iy=1,grid_y_num 
     do ix=1,grid_x_num 
      buffer(icount+ix) = grid_$X_lapl(ix,iy,iz)
     enddo
     icount = icount + grid_x_num
    enddo
    dim = grid_x_num * grid_y_num
    call MPI_REDUCE(buffer,grid_$X_lapl(1,1,iz),dim,mpi_real, &
       mpi_sum,0,MPI_COMM_WORLD,ierr)
   enddo
  enddo
 IRP_ENDIF
 call set_grid_data_$X_grad(grid_$X_grad)
 call set_grid_data_$X_lapl(grid_$X_lapl)

END_PROVIDER

SUBST [ X ]
 eplf;;
 elf;;
 density;; 

END_TEMPLATE

BEGIN_TEMPLATE

BEGIN_PROVIDER [ real, grid_$X_partition, (grid_x_num,grid_y_num,grid_z_num) ]
 implicit none
 BEGIN_DOC  
! Create the topological partition of $X
 END_DOC
 
 integer :: ix, iy, iz
 real, parameter :: UNDEFINED=123456789.123456789

 if (mpi_master) then
   print *, 'Attribution : $X' 

   grid_$X_partition(1,1,1) = UNDEFINED
   call get_grid_data_$X_partition(grid_$X_partition)
   if (grid_$X_partition(1,1,1) /= UNDEFINED) then
     return
   endif

   do iz=1,grid_z_num 
    do iy=1,grid_y_num 
     do ix=1,grid_x_num 
      grid_$X_partition(ix,iy,iz) = 0.
     enddo
    enddo
   enddo

   integer :: icount, imax
   imax = 100
   do iz=1,grid_z_num 
    do iy=1,grid_y_num 
     do ix=1,grid_x_num 
      call find_attribution(ix,iy,iz,grid_$X,grid_$X_partition,imax)
     enddo
    enddo
   enddo

   call set_grid_data_$X_partition(grid_$X_partition)
   print *,  int(current_partition_index), ' basins found'
 endif

END_PROVIDER

SUBST [ X ]
 elf;;
 eplf;;
 density;;

END_TEMPLATE

BEGIN_PROVIDER [ real, current_partition_index ]
 implicit none
 BEGIN_DOC  
! Current number for the space partitioning
 END_DOC
 current_partition_index = 0.
END_PROVIDER

recursive subroutine find_attribution(ix,iy,iz,g,a,depth)
 implicit none

 integer, intent(in) :: ix, iy, iz
 real, intent(in) :: g(grid_x_num,grid_y_num,grid_z_num)
 real, intent(inout) :: a(grid_x_num,grid_y_num,grid_z_num)
 integer, intent(in) :: depth

 real :: buffer(-1:1,-1:1,-1:1)

 integer :: i,j,k

 if (depth == 0) then
   return
 endif
 if (a(ix,iy,iz) /= 0.) then
  return
 else if (ix == 1) then
  call find_attribution(2,iy,iz,g,a,depth-1)
  a(ix,iy,iz) = a(2,iy,iz)
 else if (ix == grid_x_num) then
  call find_attribution(ix-1,iy,iz,g,a,depth-1)
  a(ix,iy,iz) = a(ix-1,iy,iz)
 else if (iy == 1) then
  call find_attribution(ix,2,iz,g,a,depth-1)
  a(ix,iy,iz) = a(ix,2,iz)
 else if (iy == grid_y_num) then
  call find_attribution(ix,iy-1,iz,g,a,depth-1)
  a(ix,iy,iz) = a(ix,iy-1,iz)
 else if (iz == 1) then
  call find_attribution(ix,iy,2,g,a,depth-1)
  a(ix,iy,iz) = a(ix,iy,2)
 else if (iz == grid_z_num) then
  call find_attribution(ix,iy,iz-1,g,a,depth-1)
  a(ix,iy,iz) = a(ix,iy,iz-1)

 else

   ! Copy the environment of the current point
   do i=-1,1
    do j=-1,1
     buffer(-1,j,i) = g(ix-1, iy+j, iz+i)
     buffer( 0,j,i) = g(ix  , iy+j, iz+i)
     buffer( 1,j,i) = g(ix+1, iy+j, iz+i)
    enddo
   enddo 

   ! Find the index of the maximum value
   real :: m
   integer :: im(3)
   m  = buffer(0,0,0)
   real :: average, variance
   im(1) = 0
   im(2) = 0
   im(3) = 0
   average = 0.
   variance = 0.
   do i=-1,1
    do j=-1,1
     do k=-1,1
      if (buffer(k,j,i) > m) then
       m = buffer(k,j,i)
       im(1) = k
       im(2) = j
       im(3) = i
      endif
      average += buffer(k,j,i)
      variance += buffer(k,j,i)**2
     enddo
    enddo
   enddo 
   variance = variance / 27.
   average = average / 27.
   variance = (variance - average**2)

   ! Avoid strict equality between points
   buffer(0,0,0) += sqrt(variance/26.)


   if ( (im(1) == 0).and.(im(2) == 0).and.(im(3) == 0) ) then
     ! If the maximum is at the center, a new maximum is found
     if (a(ix,iy,iz) == 0.) then
         current_partition_index += 1.
         a(ix,iy,iz) = current_partition_index
         real :: x,y,z
         call int_to_cart(ix,iy,iz,x,y,z)
         print '(3(2X,F10.7))', x,y,z
     endif
   else
     ! Otherwise, get the partition index of the max value
     call find_attribution(ix+im(1), iy+im(2), iz+im(3),g,a,depth-1)
     a(ix,iy,iz) = a(ix+im(1), iy+im(2), iz+im(3))
   endif

 endif

end

subroutine int_to_cart(ix,iy,iz,x,y,z)
 implicit none
 integer, intent(in)  :: ix, iy, iz
 real, intent(out) :: x, y, z
  x = grid_origin(1)+(ix-1)*grid_step(1)
  y = grid_origin(2)+(iy-1)*grid_step(2)
  z = grid_origin(3)+(iz-1)*grid_step(3)
end