      module adaptive
      use types
    
      use mesh

      type ADAPTIVE_MESH

      integer(I4B) :: g
      integer(I4B) :: level
      integer(I4B), dimension(3)  :: gnr   ! grid's size
      integer(I4B), dimension(3)  :: gflb  ! grid's left boundary type
      integer(I4B), dimension(3)  :: gfrb  ! grid's right bounary type
      real(SP), dimension(3)   :: grmin ! coordinate of grid's left boundary
      real(SP), dimension(3)   :: grmax ! coordinate of grid's right boundary
      type (ADAPTIVE_MESH), pointer :: child
      type (ADAPTIVE_MESH), pointer :: sibiling
      type (ADAPTIVE_MESH), pointer :: parent
      type (ADAPTIVE_MESH), pointer :: neighbor
      type (GF3TL), pointer :: grid

      end type ADAPTIVE_MESH

      contains

      subroutine initialize_grid(p,level,grid_name,xmin,ymin,zmin, &
                                 xmax,ymax,zmax,dr)
      type(ADAPTIVE_MESH), pointer  :: p
      integer(I4B)                  :: level
      character(200)                :: grid_name 
      real(SP)                      :: xmin,ymin,zmin
      real(SP)                      :: xmax,ymax,zmax
      real(SP),dimension(3)         :: dr

      character(256), target :: name 
      character(1) :: a,b,c

      allocate(p,stat=ierr)
      if(ierr>0) then 
          write(6,*) 'INITIALIZE GRID(1): not enough memory'
          stop
      endif 

      p%g=get_g(0)
      p%level=level
      p%gflb=1
      p%gfrb=1
      p%grmin=(/xmin,ymin,zmin/)
      p%grmax=(/xmax,ymax,zmax/)

      p%gnr=1.0
      p%gnr(1)=int((xmax-xmin+dr(1))/dr(1) + 0.5)
      if(maxdim>1) p%gnr(2)=int((ymax-ymin+dr(2))/dr(2) + 0.5)
      if(maxdim>2) p%gnr(3)=int((zmax-zmin+dr(3))/dr(3) + 0.5)

write(71,*) 'p%gnr',p%gnr
write(71,*) ymax,ymin,dr(2)

      nullify (p%child)
      nullify (p%sibiling)
      nullify (p%parent)
      nullify (p%neighbor)
      allocate(p%grid,stat=ierr)
      if(ierr>0) then
          write(6,*) 'INITIALIZE GRID(2): not enough memory'
          stop
      endif 

       a = CHAR(p%g/100+48)        
       b = CHAR(MOD(p%g,100)/10+48)
       c = CHAR(MOD(p%g,10)+48)
       name = 'grid_' // a // b // c  // ' :: ' //grid_name
!       write(6,*) trim(name)
 
      call allocate_grid(p%grid,name,p%grmin,p%grmax,dr)

      return
      end subroutine initialize_grid

subroutine coarse_to_fine(grid,ngf,rho)
!----------------------------------------------------------------------
! AUTHOR: Tom Haupt, NPAC at Syracuse University, haupt@npac.syr.edu
!----------------------------------------------------------------------
type(ADAPTIVE_MESH), pointer :: grid
integer(I4B) ngf,rho

integer(I4B) i,j,k,i1,j1,k1,ni,nj,nk
type(ADAPTIVE_MESH), pointer :: coarse


!find parent
coarse=>grid%parent
IF(.not.associated(coarse)) THEN
    write(6,*) 'align_meshes: orfant grid detected!'
    stop
ENDIF

i=(grid%grmin(1)-coarse%grmin(1))/rho + 1
j=(grid%grmin(2)-coarse%grmin(2))/rho + 1
k=(grid%grmin(3)-coarse%grmin(3))/rho + 1

i0=1
j0=1
k0=1

i1=grid%gnr(1)
j1=grid%gnr(2)
k1=grid%gnr(3)

ni=(i1+1)/rho
nj=(j1+1)/rho
nk=(k1+1)/rho

write(71,*) 'COARSE'
write(71,*)  coarse%grid%gfcn(ngf)%tlev(1)%data(5,:,1)
write(71,*) i,j,k
write(71,*) i1,j1,k1
write(71,*) ni,nj,nk
write(71,*) rho


    grid%grid%gfcn(ngf)%tlev(1)%data=0.0
    grid%grid%gfcn(ngf)%tlev(1)%data(i0:i1:rho,j0:j1:rho,k0:k1:rho)= &
              coarse%grid%gfcn(ngf)%tlev(1)%data(i:ni,j:nj,k:nk)

write(71,*) 'AND FINE'
write(71,*)  grid%grid%gfcn(ngf)%tlev(1)%data

return
end subroutine coarse_to_fine

subroutine align_grids(grid,idim,istride,ibeg,iend)
!----------------------------------------------------------------------
! To transfer data form grids on differen refinement levels, one 
! have to align them first. This utility routine provides this
! functionality.
!
! Given a grid on a fine level and a refinement ratio (stride)
! this routine returns indices of its parent that correspond to
! the grid boundary (ibeg,iend) along dimension idim
!
! This routine is called by inject 
!
! AUTHOR: Tom Haupt, NPAC at Syracuse University, haupt@npac.syr.edu
!----------------------------------------------------------------------

type(ADAPTIVE_MESH), pointer :: grid
integer(I4B) idim,istride,ibeg,iend

type(ADAPTIVE_MESH), pointer :: coarse

!find parent
coarse=>grid%parent
IF(.not.associated(coarse)) THEN
    write(6,*) 'align_meshes: orfant grid detected!'
    stop
ENDIF

!     returns indices on the coarse grids    

!      write(6,*) idim,coarse%gnr(idim)
!      write(6,*) coarse%grmax(idim),coarse%grmin(idim)
      IF(abs(coarse%grmax(idim)-coarse%grmin(idim))<1e-10) THEN
!       one- or two-dimensionl grids
        ibeg=1
        iend=1
        return
      ELSE         
      dr=(coarse%grmax(idim)-coarse%grmin(idim))/(coarse%gnr(idim)-1)
      ibeg=(grid%grmin(idim)-coarse%grmin(idim))/dr+1
      iend=(grid%gnr(idim)-1)/istride + ibeg
!      write(6,*) 'ALIGN:',ibeg,iend
      return
      ENDIF 
      
      end subroutine align_grids

subroutine linject(head,curlev,is)
      type(ADAPTIVE_MESH),pointer :: head
      integer(I4B) curlev,is

      integer(I4B) ix,iy,iz,jx,jy,jz
      integer(I4B) nx,ny,nz
      type (ADAPTIVE_MESH), pointer :: coarse,fine

      fine=>head

      nx=fine%gnr(1)
      ny=fine%gnr(2)
      nz=fine%gnr(3)

      DO WHILE(associated(fine))
        coarse => fine%parent
        call align_grids(fine,1,is,ix,jx)
        call align_grids(fine,2,is,iy,jy)
if(maxdim>2) then
  call align_grids(fine,3,is,iz,jz)
else
 iz=1
 jz=1
endif
  
i0=1
j0=1
z0=1

write(69,*) 'LINJECT'
write(69,*) nx,ny,nz
write(69,*) ix,jx,iy,jy
write(69,*) is, 'igf',coarse%grid%n_gfcn
write(69,*) fine%grid%gfcn(1)%tlev(1)%data(1:nx:is,1:ny:is,1:nz:is)
write(69,*) 'all'
write(69,*) fine%grid%gfcn(1)%tlev(1)%data

!do not transfer edges!!!
ix=ix+1
jx=jx-1
iy=iy+1
jy=jy-1
i0=i0+is
nx=nx-is
j0=j0+is
ny=ny-is

write(69,*) 'before'
write(69,*) coarse%grid%gfcn(1)%tlev(1)%data(5,:,1)

!this loop is unsafe!!!

        DO igf=1,coarse%grid%n_gfcn-1
         DO it=1,coarse%grid%gfcn(igf)%n_t
           coarse%grid%gfcn(igf)%tlev(it)%data(ix:jx,iy:jy,iz:jz) = &
           fine%grid%gfcn(igf)%tlev(it)%data(i0:nx:is,j0:ny:is,z0:nz:is)
         ENDDO
        ENDDO
        fine=>fine%neighbor
      ENDDO

write(69,*) 'after'
write(69,*) coarse%grid%gfcn(1)%tlev(1)%data(5,:,1)

      return
      end subroutine linject


      subroutine find_parent(first,parent,clmin,clmax)
      type(ADAPTIVE_MESH),pointer :: first,parent
      real(SP), dimension(3) :: clmin,clmax
   
      type(ADAPTIVE_MESH),pointer :: p

      nullify(p)
      p=>first

write(6,*) 'CLMIN,CLMAX',clmin,clmax

      DO WHILE(associated(p))
write(6,*) 'PARENT min ',p%grmin
write(6,*) 'PARENT max ',p%grmax

       IF(clmin(1).ge.p%grmin(1).and.clmax(1).le.p%grmax(1).and. &  
          clmin(2).ge.p%grmin(2).and.clmax(2).le.p%grmax(2).and. &   
          clmin(3).ge.p%grmin(3).and.clmax(3).le.p%grmax(3)) THEN   
       EXIT
       ENDIF
       p=>p%neighbor
      ENDDO

      parent=>p
      return
      end subroutine find_parent

subroutine initialize_data(child,parent,old_head,drc,drp,dt,is)
!--------------------------------------------------------------------------
! User and/or application dependent code is burried inside routines
! init_gf 
!
! TO BE CLEANED UP !!!
!
! AUTHOR: Tom Haupt, NPAC at Syracuse University, haupt@npac.syr.edu
!--------------------------------------------------------------------------

type(ADAPTIVE_MESH),pointer    :: child,parent,old_head
real(SP), dimension(3)         :: drc,drp
real(SP)                       :: dt 
integer(I4B) is

real(SP) dt0,time,dt1

      INTERFACE
       subroutine u_init_gf(grid,t,dt,dr)
       use amr
       type(GF3TL),pointer :: grid
       real(SP)              :: t,dt
       real(SP),dimension(3) :: dr
       end subroutine u_init_gf

       subroutine interpolateH3(parentptr,childptr,drc)
       use types
       use amr
       type(GF3TL), pointer :: parentptr,childptr
       real(SP), dimension(3)   :: drc   !dr of child     
       end subroutine interpolateH3

      END INTERFACE

!      write(6,*) 'start initialize_data',dt
      dt1=dt   
      do i=1,parent%grid%n_gfcn
       dt0=0
        do j=1,parent%grid%gfcn(i)%n_t 
        child%grid%gfcn(i)%t(1)=parent%grid%gfcn(i)%t(1)-dt1
        dt1=dt0-dt1
        enddo 
      enddo

      time=child%grid%gfcn(1)%t(1)

!c     the following is specific for AD code
!c          do j=1,child%grid%gfcn(R2)%tlev(1)%c(1)%n_coord
!c          child%grid%gfcn(R2)%tlev(1)%data(j,1,1)=
!c     &    child%grid%gfcn(R2)%tlev(1)%c(1)%coord(j)**2      
!c          child%grid%gfcn(R3)%tlev(1)%data(j,1,1)=
!c     &    child%grid%gfcn(R2)%tlev(1)%c(1)%coord(j)*
!c     &    child%grid%gfcn(R2)%tlev(1)%data(j,1,1)       
!c          enddo
         
!c      --------------------------------------------------------
!c      IS THIS THE INITIAL TIME?
!c      --------------------------------------------------------
       if(abs(child%grid%gfcn(1)%t(1)-t_init).lt.1.0d-10) then
!        force proper boundary conditions 
! Note: Do we need it?        call rginterp(child,parent,drc,drp,is)
         call u_init_gf(child%grid,time,dt,drc) !from first principles     
       else
!! used to be (AD)  call rginterp(child,parent,drc,drp,is)
         write(6,*) 'call interpolateH3'
!         call interpolateH3(parent%grid,child%grid,drc)
!         write(6,*) 'return from interpolateH3'
!         do i=35,47
!          write(6,*) 'interp',i,child%grid%gfcn(i)%tlev(1)%data(3,3,3),
!     &    parent%grid%gfcn(i)%tlev(1)%data(3,3,3)
!         enddo
!         write(6,*) 'call rgtran' 
!        move data from old level if exist
! NOTE:  to ba done later         call rgtran(child,old_head,drc)
         write(6,*) 'WARNING: rgtran commented out'
       endif

      return
      end subroutine initialize_data

!     ---------------------------------------------------------
!       GRID ID SERVER
!     ---------------------------------------------------------

      integer(I4B) function get_g(iflag)
      integer(I4B) iflag
      integer(I4B),dimension(1000),SAVE :: g_table

      if(iflag.eq. -1) then
        g_table=-1
        get_g=0
        return
      endif   
      if(iflag.eq.0) then
        get_g=1 
        DO WHILE(g_table(get_g).gt.0)
          get_g=get_g+1
        enddo
        g_table(get_g)=get_g     
        return
      endif   
      if(iflag.gt.0) then
        g_table(iflag)=-1
        get_g=-1
      endif

      return
      end function get_g
      
      end module adaptive










