
      module adaptive
    
      use mesh

      type ADAPTIVE_MESH

      integer :: g
      integer :: level
      integer, dimension(3)  :: gnr   ! grid's size
      integer, dimension(3)  :: gflb  ! grid's left boundary type
      integer, dimension(3)  :: gfrb  ! grid's right bounary type
      real, dimension(3)   :: grmin ! coordinate of grid's left boundary
      real, 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                      :: level
      character(200)               :: grid_name 
      real                       :: xmin,ymin,zmin
      real                       :: xmax,ymax,zmax
      real,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)=(xmax-xmin)/dr(1) + 1.5
      p%gnr(2)=(ymax-ymin)/dr(2) + 1.5
      p%gnr(3)=(zmax-zmin)/dr(3) + 1.5
      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
c       write(6,*) trim(name)
 
      call allocate_grid(p%grid,name,p%grmin,p%grmax,dr)

      return
      end subroutine initialize_grid

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

      type(ADAPTIVE_MESH), pointer :: coarse

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

!     returns indices on the coarse grids    

c      write(6,*) idim,coarse%gnr(idim)
c      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
c      write(6,*) 'ALIGN:',ibeg,iend
      return
      ENDIF 
      
      end subroutine align_grids

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

      integer ix,iy,iz,jx,jy,jz
      integer 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)
        call align_grids(fine,3,is,iz,jz)

        DO igf=1,6
c        DO igf=1,coarse%grid%n_gfcn-2
           coarse%grid%gfcn(igf)%tlev(1)%data(ix:jx,iy:jy,iz:jz) =
     &     fine%grid%gfcn(igf)%tlev(1)%data(1:nx:is,1:ny:is,1:nz:is)
        ENDDO
        fine=>fine%neighbor
      ENDDO

      return
      end subroutine linject


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

      write(6,*) 'x',clmin(1),clmax(1)
c      write(6,*) 'y',clmin(2),clmax(2)
c      write(6,*) 'z',clmin(3),clmax(3)
      nullify(p)
      p=>first
      DO WHILE(associated(p))
       write(6,*) p%grmin(1),p%grmax(1) 
       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)
      type(ADAPTIVE_MESH),pointer  :: child,parent,old_head
      real, dimension(3)         :: drc,drp
      real                       :: dt 
      integer is

      real dt0,time

      INTERFACE
       subroutine init_gf(grid,t,dt,dr)
       use bh02
       type(GF3TL),pointer :: grid
       real              :: t,dt,dr
       end subroutine init_gf

      END INTERFACE
   
      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)-dt0
        dt=dt0-dt
        enddo 
      enddo

      time=child%grid%gfcn(1)%t(1)
          do j=1,child%grid%gfcn(R2)%tlev(1)%c(1)%n_coord
          child%grid%gfcn(R2)%tlev(1)%data(j,1,1)=
     &    child%grid%gfcn(R2)%tlev(1)%c(1)%coord(j)**2      
          child%grid%gfcn(R3)%tlev(1)%data(j,1,1)=
     &    child%grid%gfcn(R2)%tlev(1)%c(1)%coord(j)*
     &    child%grid%gfcn(R2)%tlev(1)%data(j,1,1)       
          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 
         call rginterp(child,parent,drc,drp,is)
         call init_gf(child%grid,time,dt,drc(1)) !from first principles     
       else
         call rginterp(child,parent,drc,drp,is)
!        move data from old level if exist
         call rgtran(child,old_head,drc)
       endif

      return
      end subroutine initialize_data

       subroutine rgtran(new,old,ldr)  
       type(ADAPTIVE_MESH), pointer :: new
       type(ADAPTIVE_MESH), pointer :: old
       real, dimension(3) :: ldr

       type(ADAPTIVE_MESH), pointer :: o
       real dr
c       integer align_m,lev
       integer lev
                
       nullify(o)

c       if(new%level.ne.old%level) then
c         write(6,*) 'rgtran: levels does not match!'
c         return
c       endif

       o=>old
       lev=new%level
       dr=ldr(1)
       nr=new%gnr(1)
 

c      ------------------------------------------------------------
c      FOR EACH GRID IN OLD LIST ...
c      ------------------------------------------------------------
        
       DO WHILE(associated(o))
          
c        ----------------------------------------------------------
c        ... IF THERE IS ANY OVERLAP WITH GRID G ...
c        ----------------------------------------------------------
         if(o%grmax(1).ge.new%grmax(1). and.
     &      o%grmin(1).le.new%grmin(1)) then  

c          --------------------------------------------------------
c          ... COMPUTE OVERLAP INDICES ...
c          --------------------------------------------------------
           if(new%grmin(1).ge.o%grmin(1)) then
             ig1 = 1
             itg1= align_m(new%grmin(1),o%grmin(1),dr)
           else
             ig1= align_m(o%grmin(1),new%grmin(1),dr)
             igt1=1
           endif
        
           if(new%grmax(1).ge.o%grmax(1)) then
             ig2=nr+1-align_m(new%grmax(1),o%grmax(1),dr) 
           else
             ig2=nr
           endif
           ig2=ig2+1

c          --------------------------------------------------------
c          ... AND COPY ALL FCNS, WHICH NEED TO BE COPIED
c          --------------------------------------------------------

          ntran=6 !all grid functions but ESTRE,R2 and R3

          do i=1,ntran
           do it=1,new%grid%gfcn(i)%n_t 
           new%grid%gfcn(i)%tlev(it)%data(itg1:ig2,1,1) 
     &     =o%grid%gfcn(i)%tlev(it)%data(ig1:ig2,1,1) 
           enddo
          enddo

          endif 

          o=>o%neighbor               
        

       ENDDO

       return
       end subroutine rgtran
c ---------------------------------------------------------------------
 

       integer function align_m(a1,a2,ainc)
       real a1,a2,ainc
       real, parameter :: tol=1.0d-06
    
          align_m=int((a1-a2)/ainc) + 1
          if(a1-(a2+align_m*ainc).gt.tol*ainc) then
            write(6,*) 'align: bad parameters'
            stop
          endif
      return
      end function align_m

      integer function get_g(iflag)
      integer iflag
      integer,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

       subroutine rginterp(child,parent,drc,drp,is)
       type(ADAPTIVE_MESH),pointer :: child,parent
       real,dimension(3)         :: drc,drp  
       integer is

       real                      :: dr,dr2
       integer nr,nr2,irho,nint,it,tlev

       integer fp 
       
       interface
         subroutine dvi2q1(v1,v2,n1,irho)
         implicit       none
         integer        n1,          irho
         real         v1(:),      v2(:)
         real         rhom1 
         end subroutine dvi2q1

         subroutine dvsvsa(v1,v2,v3,s1,s2,n)
         real      v1(:), v2(:), v3(:)
         real     s1, s2
         integer     n
         end subroutine dvsvsa

       end interface

!      reduction to 1-d       

       dr = drc(1) 
       dr2= drp(1)
       irho=int(dr2/dr)
       nr=child%gnr(1)
       nr2=parent%gnr(1)



       nint=6 !PPN,PIN,PHN,MASS,MOVR,DMDR; skip ESTRE,R2 and R3

       do i=1,nint
       it=child%grid%gfcn(i)%n_t

c      --------------------------------------------------------
c      SPATIALLY INTERPOLATE MOST ADVANCED DATA 
c      --------------------------------------------------------

c          (linear interpolation)

          call align_grids(child,1,is,ibeg,iend)

          call dvi2q1(
     &      parent%grid%gfcn(i)%tlev(1)%data(ibeg:iend,1,1),
     &      child%grid%gfcn(i)%tlev(1)%data(:,1,1),  
     &      iend-ibeg+1,irho)
c        write(6,*) 'time level 1'

c      --------------------------------------------------------
c      SPATIALLY INTERPOLATE RETARDED DATA
c      --------------------------------------------------------
       if(it>1) then 
!      only two time levels assumed          
c          do ilev=it,2,-2
          ilev=2
c         ospint is set to 2: linear interpolation
c       write(6,*)   parent%grid%gfcn(i)%tlev(ilev)%data(ibeg:iend,1,1) 
           call dvi2q1(
     &      parent%grid%gfcn(i)%tlev(ilev)%data(ibeg:iend,1,1),
     &      child%grid%gfcn(i)%tlev(ilev)%data(:,1,1),             
     &      iend-ibeg+1,irho)
c          enddo 

c      --------------------------------------------------------
c       INTERPOLATE TEMPORALLY
c      --------------------------------------------------------

c       write(6,*) child%grid%gfcn(i)%tlev(2)%data(:,1,1) 

           child%grid%gfcn(i)%tlev(1)%data(:,1,1)=
     & (1.0d0 - dr/dr2)*child%grid%gfcn(i)%tlev(2)%data(:,1,1) 
     &       + (dr/dr2)*child%grid%gfcn(i)%tlev(1)%data(:,1,1) 
       endif
       enddo 

c       write(6,*) 'return from rginterp'
       return
       end subroutine rginterp

      end module adaptive










