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

 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

 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

IRP_IF caca
BEGIN_PROVIDER [ integer, grid_partition_$X, (grid_x_num,grid_y_num,grid_z_num) ]
 implicit none
 BEGIN_DOC  
! Create the topological partition of $X
 END_DOC
 
END_PROVIDER
IRP_ENDIF