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 get_grid_point_num(Npoints)
 call get_grid_origin(origin)
 call get_grid_opposite(opposite)
 call get_grid_step_size(step_size)

END_PROVIDER

!subroutine write_grid_eplf
! implicit none
! integer :: i
! integer :: l
! integer :: ix, iy, iz
! if (.not.mpi_master) then
!  return
! endif
! open(unit=99,file=grid_cube_filename,status='UNKNOWN',action='WRITE')
! write (99,*) 'Cube File'
! write (99,*) 'Analytical EPLF grid'
! write (99,10) nucl_num,(grid_eplf_origin(i), i=1,3)
! write (99,10) grid_eplf_x_num, grid_eplf_step(1), 0., 0.
! write (99,10) grid_eplf_y_num, 0., grid_eplf_step(2), 0.
! write (99,10) grid_eplf_z_num, 0., 0., grid_eplf_step(3)
! do i=1,nucl_num
!  write (99,11) int(nucl_charge(i)), nucl_charge(i), (nucl_coord(i,l),l=1,3)
! enddo
! do ix = 1, grid_eplf_x_num
!  do iy = 1, grid_eplf_y_num
!    write (99,20) (grid_eplf(ix,iy,iz), iz=1, grid_eplf_z_num)
!  enddo
! enddo
! 10   format (2X,I3,3(2X,F10.6))
! 11   format (2X,I3,4(2X,F10.6))
! 20   format (6(E13.5))
! close(99)
!end    

BEGIN_SHELL [ /usr/bin/python ]

grids = [ \
 "eplf",
 "elf",
 "density", 
]

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

 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 = 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)
    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)
   enddo
 IRP_ENDIF

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

 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

END_PROVIDER

subroutine set_grid_data_$X(buffer)
 real   :: buffer(grid_x_num,grid_y_num,grid_z_num)
 if (mpi_master) then
   call ezfio_set_grid_data_$X(buffer)
 endif
end

subroutine set_grid_data_$X_grad(buffer)
 real   :: buffer(grid_x_num,grid_y_num,grid_z_num,4)
 if (mpi_master) then
   call ezfio_set_grid_data_$X_grad(buffer)
 endif
end

subroutine set_grid_data_$X_lapl(buffer)
 real   :: buffer(grid_x_num,grid_y_num,grid_z_num)
 if (mpi_master) then
   call ezfio_set_grid_data_$X_lapl(buffer)
 endif
end
"""

for grid in grids: 
  print template.replace("$X",grid)

END_SHELL