module bh01     

type GF3
character(32), pointer          :: name
integer, pointer                :: n_dim
integer, pointer, dimension(:)  :: dim
real, pointer, dimension(:,:,:) :: data
real, pointer, dimension(:)     :: coord
end type GF3 

type GF3T
integer :: nt_levels
real, pointer, dimension(:) :: integration_time
type (GF3), pointer, dimension (:) :: gridfunc
end type GF3T

type GF3TL
integer :: n_func
type (GF3T), pointer, dimension (:) :: lattice
end type GF3TL 

end module bh01
!  =========== subroutine allocate ================================

subroutine allocate_grid(nf,ntlev,my_grid)
use bh01
integer, intent(in) :: nf
integer, dimension(:), intent(in)  :: ntlev 
type (GF3TL), pointer :: my_grid
character(32),target  :: myname
integer               :: ncoord
integer,target        :: ndim
integer, dimension(3), target :: dim3 = (/10,10,10/)

myname = 'any name is good'
ndim = 3

allocate(my_grid)

my_grid % n_func = nf
allocate(my_grid % lattice(my_grid % n_func))

if(size(ntlev).ne.nf) then
write(6,*) "error: size of ntlev not equal to nf"
endif

do i = 1, my_grid % n_func
my_grid % lattice(i) % nt_levels = ntlev(i)
allocate(my_grid % lattice(i) % integration_time( my_grid % lattice(i) % nt_levels))
allocate(my_grid % lattice(i) % gridfunc( my_grid % lattice(i) % nt_levels))
enddo

do i = 1, my_grid % n_func
do j= 1, my_grid % lattice(i) % nt_levels           
my_grid % lattice(i) % gridfunc(j) % name => myname  
my_grid % lattice(i) % gridfunc(j) % n_dim => ndim  
allocate(my_grid % lattice(i) % gridfunc(j) % dim(ndim))
my_grid % lattice(i) % gridfunc(j) % dim => dim3  
allocate(my_grid % lattice(i) % gridfunc(j) % data(dim3(1),dim3(2),dim3(3)))
ncoord=sum(dim3) 
allocate(my_grid % lattice(i) % gridfunc(j) % coord(ncoord)) 
enddo
enddo 

return
end
!     -------------------------------------------------------------

program NR3
use bh01
type (GF3TL), pointer  :: my_grid
integer, parameter     :: nf = 1
!      integer, dimension(nf) :: ntlev = (/3,3,2,2,1/)
integer, dimension(nf) :: ntlev

interface
subroutine allocate_grid(nf,ntlev,my_grid)
use bh01
integer, intent(in) :: nf
integer, dimension(:), intent(in)  :: ntlev 
type (GF3TL), pointer :: my_grid
character(32),target  :: myname
integer               :: ncoord
integer,target        :: ndim
end
end interface
real, parameter :: Epsilon = 0.00002
integer, parameter :: R = 10, C = 10, P = 10
real :: error = 100.0, bsum
integer :: x, y, niter, kount = 1
integer, pointer, dimension(:,:,:) :: I, J
real, pointer, dimension(:,:,:) :: Z

allocate(Z(R,C,P))
allocate(I(R,C,P))
allocate(J(R,C,P))

I = SPREAD( SPREAD([1:R],2,C), 3, P )     !Coordinates
J = SPREAD( SPREAD([1:C],1,R), 3, P )     !Coordinates

!Z is a solution of the difference equation
Z = -1.4 + 7.12*I - 9.23*J + 0.001*(I**2 - J**2)


!      for ParaSoft
ntlev = (/1/)

call allocate_grid(nf,ntlev,my_grid)

!     Install Boundaries
my_grid%lattice(1)%gridfunc(1)%data = Z

print *, "Solution range: ", minval(Z), maxval(Z)
bsum = SUM(Z(:,:,1)) + SUM(Z(:,:,P)) + SUM(Z(1,:,2:P-1)) + SUM(Z(R,:,2:P-1))+SUM(Z(2:R-1,1,2:P-1)) + SUM(Z(2:R-1,C,2:P-1))


!     Intialize interior to average on boundary
!     The "-16" prevents counting corners twice

my_grid%lattice(1)%gridfunc(1)%data(2:R-1,2:C-1,2:P-1) = bsum/(2*R*C+2*R*P+2*C*P - 16)

do while (error > Epsilon)
  kount = kount + 1
  error = update(Z,my_grid%lattice(1)%gridfunc(1)%data,R,C,P)
  if (mod(kount,20).eq. 0) print *, 'error = ', error, kount
end do

print *, 'error = ', error, kount
end program



!-----------------subroutine update-------------------------
function update(Z,grid, i, j, k)
real, dimension(i,j,k),intent(inout) :: grid
real, dimension(i,j,k),intent(in) :: Z
integer :: i, j, k, a, b, c

forall(a=2:i-1,b=2:j-1,c=2:k-1) grid(a,b,c)=(grid(a-1,b,c)+grid(a+1,b,c)+&
grid(a,b-1,c)+grid(a,b+1,c)+grid(a,b,c-1)+grid(a,b,c+1)) / 6.0

update = MAXVAL(abs(Z - grid))

return
end
