      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.001
      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

!I, J arrays are coordinates arrays only for temporay purpose.
!Z array contains exact solutions. This is used for error calculation.
      allocate(Z(R,C,P))
      allocate(I(R,C,P))
      allocate(J(R,C,P))

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

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

! Free memory of I and J arrays
      deallocate(I)
      deallocate(J)

!      for ParaSoft
       ntlev = (/1/)

      call allocate_grid(nf,ntlev,my_grid)
      
!     initialize time step 1
      do i=1,nf
        my_grid % lattice(i) % gridfunc(1) % data = 1.0
       enddo

!     for those guys who has more than 1 step defined: update

      do i=1,nf
      if (my_grid % lattice(i) % nt_levels .gt.1) then
        my_grid % lattice(i) % gridfunc(2) % data = 2.0
       endif
       enddo

!     print sample results

      do i=1,nf
       if (my_grid % lattice(i) % nt_levels .gt.1) then
        write(6,*) i,
     &   my_grid % lattice(i) % gridfunc(2) % data(1,2,3)
        else
        write(6,*) i,
     &   my_grid % lattice(i) % gridfunc(1) % data(1,2,3)
        endif 
       enddo

! Free memory of Z array
      deallocate(Z)

      end 
