module parameters
integer,parameter  :: PX=4,PY=1       !number of processors
integer, parameter :: rho=2           !grid ratio
integer, parameter :: maxlevel=5      !maximum allowed levels
real, parameter    :: dtbydr=1.0    
real, parameter    :: dx=1.0,dy=1.0,dz=1.0
end module parameters

module adaptive
use parameters
integer, parameter  :: ngf=2
integer, parameter  :: tlevels=2
integer, parameter  :: n1id=1,n2id=2,n3id=2,n4id=2,n5id=2
integer, parameter  :: g1x=10,g2x=20,g3x=40,g4x=80,g5x=160
integer, parameter  :: g1y=10,g2y=20,g3y=40,g4y=80,g5y=160
integer, parameter  :: g1z=1,g2z=1,g3z=1,g4z=1,g5z=1

real, dimension(n1id,ngf,tlevels,g1x,g1y,g1z) :: glev1
real, dimension(n2id,ngf,tlevels,g2x,g2y,g2z) :: glev2
real, dimension(n3id,ngf,tlevels,g3x,g3y,g3z) :: glev3
real, dimension(n4id,ngf,tlevels,g4x,g4y,g4z) :: glev4
real, dimension(n5id,ngf,tlevels,g5x,g5y,g5z) :: glev5

!HPF$ PROCESSORS PCS(PX,PY)
!HPF$ DISTRIBUTE (*,*,*,BLOCK,BLOCK,*) ONTO PCS :: glev1,glev2,glev3
!HPF$ DISTRIBUTE (*,*,*,BLOCK,BLOCK,*) ONTO PCS :: glev4,glev5

type ADAPTIVE_MESH
integer :: level
integer :: grid_id
type(ADAPTIVE_MESH), pointer :: neighbor
type(ADAPTIVE_MESH), pointer :: sibling
type(ADAPTIVE_MESH), pointer :: parent
type(ADAPTIVE_MESH), pointer :: child
real, dimension(3)           :: rmin
real, dimension(3)           :: rmax
integer, dimension(3)        :: nr 
real t
end type ADAPTIVE_MESH

type LEVEL_INFO
type(ADAPTIVE_MESH),pointer :: head
integer                     :: curmax
real                        :: threshold
real,dimension(0:3)         :: dr
end type LEVEL_INFO

type(LEVEL_INFO), dimension(maxlevel), SAVE :: lp

contains

subroutine allocate_LEVEL_INFO

do i=1,maxlevel
nullify(lp(i)%head)
lp(i)%curmax=0
lp(i)%dr(1)=dx/(rho**(i-1))
lp(i)%dr(2)=dy/(rho**(i-1))
lp(i)%dr(3)=dz/(rho**(i-1))
lp(i)%dr(0)=lp(i)%dr(1)*dtbydr
write(6,*) i,lp(i)%dr
enddo

lp(1)%threshold=0.25
lp(2)%threshold=0.45
lp(3)%threshold=0.65
!lp(4)%threshold=100.0
!lp(5)%threshold=100.0

return
end subroutine allocate_LEVEL_INFO


subroutine allocate_ADAPTIVE_MESH(level,parent,mesh,t,rmin,rmax,nr)
integer level
type(ADAPTIVE_MESH), pointer  :: parent, mesh
real, dimension(3) :: rmin,rmax
integer, dimension(3) :: nr
real, dimension(0:3) :: dr

type(ADAPTIVE_MESH), pointer :: work

nullify(work)

allocate(mesh)
mesh%parent=>parent
nullify(mesh%child)
nullify(mesh%neighbor)
nullify(mesh%sibling)
mesh%t=t
mesh%rmin=rmin
mesh%rmax=rmax
mesh%nr=nr
mesh%level=level 

if(associated(parent)) then
  if(associated(parent%child)) then 
     work=>parent%child

     DO WHILE (associated(work))
     if(.not.associated(work%sibling)) then
         work%sibling=>mesh
         mesh%grid_id=work%grid_id+1
         work=>work%sibling%sibling
     else
         work=>work%sibling
     endif
     ENDDO
  else   
     parent%child=>mesh
     mesh%grid_id=1
  endif
else    
  nullify(mesh%parent)
  mesh%grid_id=1
  if(level <> 1) write(6,*) 'WARNING: (1) orfant grid'
endif
lp(level)%curmax=mesh%grid_id
write(6,*) 'end of allocate_ADAPTIVE_MESH:',lp%curmax
if(.not.associated(lp(level)%head)) lp(level)%head=>mesh
return 
end subroutine allocate_ADAPTIVE_MESH

end module adaptive

program driver
use parameters
use adaptive      
real, dimension(:,:,:), pointer :: gridpointer

type(ADAPTIVE_MESH),pointer     :: mesh,parent
integer, dimension(3)           :: nr
real, dimension(3)              :: gmin,gmax,gd
character(256)                  :: grid_name
real t
integer level

INTERFACE
subroutine initialize_grid_functions(t)
use adaptive
real                                     :: t 
end subroutine initialize_grid_functions

subroutine update(anew,aold,t,dt,dr,nx,ny)
use parameters
real, dimension(:,:) :: anew,aold
real                 :: t,dt
integer              :: nx,ny
real, dimension(2)   :: dr
!HPF$ PROCESSORS PCS(PX,PY)
!HPF$ DISTRIBUTE *(BLOCK,BLOCK) ONTO *PCS :: anew,aold
end subroutine update
END INTERFACE


nullify(mesh)
nullify(parent)
call allocate_LEVEL_INFO

!initialize base grid
grid_name='test1'
gmin=1.0
gmax(1)=10.0
gmax(2)=10.0
gmax(3)=1.0
gd=1.0
level=1
t=0
nr(1)=int((gmax(1)-gmin(1)+gd(1))/gd(1) + 0.5)
nr(2)=int((gmax(2)-gmin(2)+gd(2))/gd(2) + 0.5)
nr(3)=int((gmax(3)-gmin(3)+gd(3))/gd(3) + 0.5)
write(6,*) nr

call allocate_ADAPTIVE_MESH(level,parent,mesh,t,gmin,gmax,nr)
write(6,*) mesh%grid_id

call initialize_grid_functions(t)

write(6,*) glev1(1,1,1,:,:,:)

grid_name='test2'
gmin(1)=2.0
gmin(2)=2.0
gmin(3)=1.0
gmax(1)=4.0
gmax(2)=6.0
gmax(3)=1.0
gd  =0.5
level=2

nr(1)=int((gmax(1)-gmin(1)+gd(1))/gd(1) + 0.5)
nr(2)=int((gmax(2)-gmin(2)+gd(2))/gd(2) + 0.5)
nr(3)=int((gmax(3)-gmin(3)+gd(3))/gd(3) + 0.5)

parent=>mesh
call allocate_ADAPTIVE_MESH(level,parent,mesh,t,gmin,gmax,nr)
write(6,*) mesh%grid_id

grid_name='test2'
gmin(1)=5.0
gmin(2)=7.0
gmin(3)=1.0
gmax(1)=9.0
gmax(2)=9.0
gmax(3)=1.0
gd  =0.5
level=2

nr(1)=int((gmax(1)-gmin(1)+gd(1))/gd(1) + 0.5)
nr(2)=int((gmax(2)-gmin(2)+gd(2))/gd(2) + 0.5)
nr(3)=int((gmax(3)-gmin(3)+gd(3))/gd(3) + 0.5)
call allocate_ADAPTIVE_MESH(level,parent,mesh,t,gmin,gmax,nr)
write(6,*) mesh%grid_id

do i=1,maxlevel
if(associated(lp(i)%head)) then
   mesh=>lp(i)%head
   DO WHILE (associated(mesh))
    write(6,*) mesh%rmin
    mesh=>mesh%sibling
   ENDDO
endif
enddo

call update(glev1(1,1,1,:,:,1),glev1(1,1,2,:,:,1),t,lp(1)%dr(0), &
lp(1)%dr(1:2),lp(1)%head%nr(1),lp(1)%head%nr(2))


stop
end


subroutine initialize_grid_functions(t)
!-------------------------------------------------------------------------
! this is a wrapper routine for user supplied subroutine u_init_gf
! that initialize grid functions at initial time
!
! the current version supports initialization on level 1 only
! and consequently initializes functions defined on the base grid
! pointed by lp(1)%ml
!
! AUTHOR: Tom Haupt, NPAC at Syracuse University, haupt@npac.syr.edu
!-------------------------------------------------------------------------
use adaptive
real                                     :: t 

type(ADAPTIVE_MESH),pointer              :: grid
real                                     :: dt
real, dimension(3)                       :: dr
integer                                  :: id

      INTERFACE
      subroutine u_init_gf(id,t,dt,dr)
      use adaptive
      integer             :: id 
      real                :: t,dt
      real, dimension(3)  :: dr
      end subroutine u_init_gf
      END INTERFACE


!initialize base grid, i.e. level=1, id=1
id=1
call u_init_gf(id,t,lp(1)%dr(0),lp(1)%dr(1:3))

return 
end subroutine initialize_grid_functions

subroutine u_init_gf(id,t,dt,dr)
!---------------------------------------------------------------------------
! At initial time (t) set of grid functions defined in GF3TL structure
! is initialized. Some of the function may have to be initialized on 
! more than one time level
!
! User is responsible for providing this routine.
! --------------------------------------------------------------------------

use adaptive
integer             :: id
real                :: t,dt
real, dimension(3)  :: dr

!VERSION for TRANSPORT demo
!set to zero all grid functions at all time levels

glev1(id,:,:,:,:,:)=0

return  
end subroutine u_init_gf

subroutine update(anew,aold,t,dt,dr,nx,ny)
use parameters
!--------------------------------------------------------------------
!user function called from user supplied update driver: gridstep
!VERSION for TRANSPORT demo
!--------------------------------------------------------------------
real, dimension(:,:)   :: anew,aold
real                   :: t,dt
integer                :: nx,ny
real, dimension(2)     :: dr

!HPF$ PROCESSORS PCS(PX,PY)
!HPF$ DISTRIBUTE *(BLOCK,BLOCK) ONTO *PCS :: anew,aold

real, parameter :: cx=1.0,cy=1.0
real  sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1
write(6,*) 'start update'

 sx=cx*dt/dr(1)
 sy=cy*dt/dr(2)
 sx2=sx**2
 sy2=sy**2
 sxy=sx*sy
 sx1=1.0/(1.0+sx)
 sy1=1.0/(1.0+sy)
 sxy1=1.0/(1.0+sx+sy)

write(6,*) dr,dt,t
write(6,*) nx,ny
write(6,*) 'AOLD'
write(6,*) aold(5,:), '*'

anew(2:nx-1,2:ny-1)=aold(2:nx-1,2:ny-1) &
  -(aold(3:nx,2:ny-1)-aold(1:nx-2,2:ny-1))*sx/2.0 &
  -(aold(2:nx-1,3:ny)-aold(2:nx-1,1:ny-2))*sy/2.0 &
  +(aold(3:nx,2:ny-1)-2.0*aold(2:nx-1,2:ny-1)     &
  +aold(1:nx-2,2:ny-1))*sx2/2.0                   &
  +(aold(2:nx-1,3:ny)-2.0*aold(2:nx-1,2:ny-1)+aold(2:nx-1,1:ny-2))*sy2/2.0 &
  +(aold(3:nx,3:ny)-aold(3:nx,1:ny-2)-aold(1:nx-2,3:ny) &
  +aold(1:nx-2,1:ny-2))*sxy/4.0 
   
                                                           


!Forall(i=2:nx-1,j=2:ny-1) anew(i,j)=                                        &
!       aold(i,j)                                                            &
!       -(aold(i+1,j)-aold(i-1,j))*sx/2.0-(aold(i,j+1)-aold(i,j-1))*sy/2.0   &
!       +(aold(i+1,j)-2.0*aold(i,j)+aold(i-1,j))*sx2/2.0                     &
!       +(aold(i,j+1)-2.0*aold(i,j)+aold(i,j-1))*sy2/2.0                     &
!       +(aold(i+1,j+1)-aold(i+1,j-1)-aold(i-1,j+1)+aold(i-1,j-1))*sxy/4.0 

write(6,*) 'ANEW'
write(6,*) anew(5,:), '*'

return
end









