
      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
      write(6,*) xmax,xmin
      write(6,*) dr
      write(6,*) 'initialize_grid',p%gnr
!      stop
      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 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,dt1

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

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

       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

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

      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

!      =============================================================
!      INTERPOLATION ROUTINES   (Mark Miller) 
!      =============================================================


!      ------------------------------------------------------------
!      interpolation driver
!      NOTE: HCODE SPECIFIC !!!
!      ------------------------------------------------------------

CCC       subroutine rginterp(child,parent,drc,drp,is)

      subroutine interpolateH3(parentptr,childptr,drc)

      type(GF3TL), pointer :: parentptr,childptr
      real, dimension(3)   :: drc   !dr of child     
      integer i,j,k,gfnum,idx
      integer, parameter :: gfcopnum = 35

      INTERFACE
      subroutine invert(flag, 
     &     in_xx,in_xy,in_xz,in_yy,in_yz,in_zz,
     &     rg, out_xx,out_xy,out_xz,out_yy,out_yz,out_zz)
c      integer nx,ny,nz
      integer flag

      real rg(:,:,:)
      real in_xx(:,:,:),in_xy(:,:,:),in_xz(:,:,:),
     &     in_yy(:,:,:),in_yz(:,:,:),in_zz(:,:,:) 

      real out_xx(:,:,:),out_xy(:,:,:),out_xz(:,:,:),
     &     out_yy(:,:,:),out_yz(:,:,:),out_zz(:,:,:) 
      end subroutine invert
      END INTERFACE 
      


!     THESE ARE THE GRID FUNCTIONS THAT NEED TO BE INTERPOLATED

      integer, dimension(gfcopnum) ::
     & gfcop = (/1,2,3,4,5,6,7,8,9,10,11,19,20,21,22,23,24,25,26,27,28,
     &          29,30,31,32,33,34,35,36,37,38,39,40,41,42/)


! NOTE: shouldn't be functions listed in gfcop used???
      do gfnum = 1,gfcopnum
      idx=gfcop(gfnum) 
!      write(6,*) 'interpolate function',idx
      do i = 1,childptr%gfcn(idx)%tlev(1)%ddim(1)
      do j = 1,childptr%gfcn(idx)%tlev(1)%ddim(2)
      do k = 1,childptr%gfcn(idx)%tlev(1)%ddim(3)
         call interpoint(parentptr,childptr,i,j,k,idx)

      enddo
      enddo
      enddo
      enddo

!  ---------------------------------------------------------  
!  THIS IS AKWARD PART OF H3 CODE
!  Initialize the grid functions simulating the coordinates
!  --------------------------------------------------------- 
       nx=childptr%gfcn(1)%tlev(1)%c(1)%n_coord
       ny=childptr%gfcn(1)%tlev(1)%c(2)%n_coord
       nz=childptr%gfcn(1)%tlev(1)%c(3)%n_coord
       dx=drc(1)
       dy=drc(2)
       dz=drc(3)
       write(6,*) dx,dy,dz
       write(6,*) nx,ny,nz

      do i=1,nx
      childptr%gfcn(xx)%tlev(1)%data(i,:,:)= (i-1-(nx-1.)/2.)*dx
      end do
      do j=1,ny
      childptr%gfcn(yy)%tlev(1)%data(:,j,:)= (j-1-(ny-1.)/2.)*dy
      end do
      do k=1,nz
      childptr%gfcn(zz)%tlev(1)%data(:,:,k)= (k-1-(nz-1.)/2.)*dz
      end do

!      r = sqrt(x**2 + y**2 + z**2);
      childptr%gfcn(rr)%tlev(1)%data= sqrt(
     & childptr%gfcn(xx)%tlev(1)%data**2 +
     & childptr%gfcn(yy)%tlev(1)%data**2 +
     & childptr%gfcn(zz)%tlev(1)%data**2)




      write(6,*) 'now covariant metric'
c     covariant metric (grid functions that aren't actually evolved)
      call invert(1,
     &      childptr%gfcn(uxx)%tlev(1)%data,
     &      childptr%gfcn(uxy)%tlev(1)%data,
     &      childptr%gfcn(uxz)%tlev(1)%data,
     &      childptr%gfcn(uyy)%tlev(1)%data,
     &      childptr%gfcn(uyz)%tlev(1)%data,
     &      childptr%gfcn(uzz)%tlev(1)%data,
     &      childptr%gfcn(rg)%tlev(1)%data,
     &      childptr%gfcn(gxx)%tlev(1)%data,
     &      childptr%gfcn(gxy)%tlev(1)%data,
     &      childptr%gfcn(gxz)%tlev(1)%data,
     &      childptr%gfcn(gyy)%tlev(1)%data,
     &      childptr%gfcn(gyz)%tlev(1)%data,
     &      childptr%gfcn(gzz)%tlev(1)%data)

      write(6,*) '------- INTERPOLATION DONE ------------'
      return
      end subroutine interpolateH3

c************************************************************************

      subroutine interpoint(parentptr,childptr,childnx,childny,childnz,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childnx,childny,childnz,childn(3),nlow(3)
      integer i,gf,tl,gfnum
      real*8 p(3),ddx(3)


      childn(1) = childnx
      childn(2) = childny
      childn(3) = childnz

      call bracket(parentptr,childptr,childn,nlow,gfnum)

!      write(6,*) '        ------> bracket'
      
      do i=1,3
!        TARGET COORDINATES
         p(i) = childptr%gfcn(gfnum)%tlev(1)%c(i)%coord(childn(i))
!        COORDINATES RELATIVE TO (nxlow,nylow,nzlow)
         ddx(i) = p(i) - parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(nlow(i))
      enddo

      if(ddx(1) .le. ddx(2)) then
         if(ddx(2) .le. ddx(3)) then
!           SIMPLEX1 : dx <= dy <= dz
            call simplex1(parentptr,childptr,childn,nlow,ddx,gfnum)
         else
            if(ddx(3) .le. ddx(1)) then
!              SIMPLEX5 : dz <= dx <= dy
               call simplex5(parentptr,childptr,childn,nlow,ddx,gfnum)
            else
!              SIMPLEX2 : dx <= dz <= dy
               call simplex2(parentptr,childptr,childn,nlow,ddx,gfnum)
            endif
         endif
      else
         if(ddx(3) .le. ddx(2)) then
!           SIMPLEX6 : dz <= dy <= dx
            call simplex6(parentptr,childptr,childn,nlow,ddx,gfnum)
         else
            if(ddx(3) .le. ddx(1)) then
!              SIMPLEX4: dy <= dz <= dx
               call simplex4(parentptr,childptr,childn,nlow,ddx,gfnum)
            else
!              SIMPLEX3: dy <= dx <= dz
               call simplex3(parentptr,childptr,childn,nlow,ddx,gfnum)
            endif
         endif
      endif

!      write(6,*) '        ------> simplex'

      return
      end subroutine interpoint

c **********************************************************************
      subroutine simplex1(parentptr,childptr,childn,nlow,ddx,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),gfnum
      real*8 ddx(3)
      integer svert(6,3)

c     X COEFF
      svert(1,1) = nlow(1) + 1
      svert(1,2) = nlow(2) + 1
      svert(1,3) = nlow(3) + 1
      svert(2,1) = nlow(1) 
      svert(2,2) = nlow(2) + 1
      svert(2,3) = nlow(3) + 1
c     Y COEFF
      svert(3,1) = nlow(1)
      svert(3,2) = nlow(2) + 1
      svert(3,3) = nlow(3) + 1
      svert(4,1) = nlow(1)
      svert(4,2) = nlow(2)
      svert(4,3) = nlow(3) + 1
c     Z COEFF
      svert(5,1) = nlow(1)
      svert(5,2) = nlow(2)
      svert(5,3) = nlow(3) + 1
      svert(6,1) = nlow(1)
      svert(6,2) = nlow(2)
      svert(6,3) = nlow(3)
      call dosimplex(parentptr,childptr,childn,nlow,ddx,svert,gfnum)
      return
      end subroutine simplex1

c **********************************************************************
      subroutine simplex2(parentptr,childptr,childn,nlow,ddx,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),gfnum
      real*8 ddx(3)
      integer svert(6,3)

c     X COEFF
      svert(1,1) = nlow(1) + 1
      svert(1,2) = nlow(2) + 1
      svert(1,3) = nlow(3) + 1
      svert(2,1) = nlow(1) 
      svert(2,2) = nlow(2) + 1
      svert(2,3) = nlow(3) + 1
c     Y COEFF
      svert(3,1) = nlow(1)
      svert(3,2) = nlow(2) + 1
      svert(3,3) = nlow(3) 
      svert(4,1) = nlow(1)
      svert(4,2) = nlow(2)
      svert(4,3) = nlow(3)
c     Z COEFF
      svert(5,1) = nlow(1)
      svert(5,2) = nlow(2) + 1
      svert(5,3) = nlow(3) + 1
      svert(6,1) = nlow(1)
      svert(6,2) = nlow(2) + 1
      svert(6,3) = nlow(3)
      call dosimplex(parentptr,childptr,childn,nlow,ddx,svert,gfnum)
      return
      end subroutine simplex2

c **********************************************************************
      subroutine simplex3(parentptr,childptr,childn,nlow,ddx,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),gfnum
      real*8 ddx(3)
      integer svert(6,3)

c     X COEFF
      svert(1,1) = nlow(1) + 1
      svert(1,2) = nlow(2) 
      svert(1,3) = nlow(3) + 1
      svert(2,1) = nlow(1) 
      svert(2,2) = nlow(2) 
      svert(2,3) = nlow(3) + 1
c     Y COEFF
      svert(3,1) = nlow(1) + 1
      svert(3,2) = nlow(2) + 1
      svert(3,3) = nlow(3) + 1
      svert(4,1) = nlow(1) + 1
      svert(4,2) = nlow(2)
      svert(4,3) = nlow(3) + 1
c     Z COEFF
      svert(5,1) = nlow(1)
      svert(5,2) = nlow(2) 
      svert(5,3) = nlow(3) + 1
      svert(6,1) = nlow(1)
      svert(6,2) = nlow(2)
      svert(6,3) = nlow(3)
      call dosimplex(parentptr,childptr,childn,nlow,ddx,svert,gfnum)
      return
      end subroutine simplex3

c **********************************************************************
      subroutine simplex4(parentptr,childptr,childn,nlow,ddx,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),gfnum
      real*8 ddx(3)
      integer svert(6,3)

c     X COEFF
      svert(1,1) = nlow(1) + 1
      svert(1,2) = nlow(2) 
      svert(1,3) = nlow(3) 
      svert(2,1) = nlow(1) 
      svert(2,2) = nlow(2) 
      svert(2,3) = nlow(3) 
c     Y COEFF
      svert(3,1) = nlow(1) + 1
      svert(3,2) = nlow(2) + 1
      svert(3,3) = nlow(3) + 1
      svert(4,1) = nlow(1) + 1
      svert(4,2) = nlow(2)
      svert(4,3) = nlow(3) + 1
c     Z COEFF
      svert(5,1) = nlow(1) + 1
      svert(5,2) = nlow(2) 
      svert(5,3) = nlow(3) + 1
      svert(6,1) = nlow(1) + 1
      svert(6,2) = nlow(2)
      svert(6,3) = nlow(3)
      call dosimplex(parentptr,childptr,childn,nlow,ddx,svert,gfnum)
      return
      end subroutine simplex4

c **********************************************************************
      subroutine simplex5(parentptr,childptr,childn,nlow,ddx,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),gfnum
      real*8 ddx(3)
      integer svert(6,3)

c     X COEFF
      svert(1,1) = nlow(1) + 1
      svert(1,2) = nlow(2) + 1
      svert(1,3) = nlow(3) 
      svert(2,1) = nlow(1) 
      svert(2,2) = nlow(2) + 1
      svert(2,3) = nlow(3) 
c     Y COEFF
      svert(3,1) = nlow(1) 
      svert(3,2) = nlow(2) + 1
      svert(3,3) = nlow(3) 
      svert(4,1) = nlow(1)
      svert(4,2) = nlow(2)
      svert(4,3) = nlow(3)
c     Z COEFF
      svert(5,1) = nlow(1) + 1
      svert(5,2) = nlow(2) + 1
      svert(5,3) = nlow(3) + 1
      svert(6,1) = nlow(1) + 1
      svert(6,2) = nlow(2) + 1
      svert(6,3) = nlow(3)
      call dosimplex(parentptr,childptr,childn,nlow,ddx,svert,gfnum)
      return
      end subroutine simplex5

c **********************************************************************
      subroutine simplex6(parentptr,childptr,childn,nlow,ddx,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),gfnum
      real*8 ddx(3)
      integer svert(6,3)

c     X COEFF
      svert(1,1) = nlow(1) + 1
      svert(1,2) = nlow(2) 
      svert(1,3) = nlow(3) 
      svert(2,1) = nlow(1) 
      svert(2,2) = nlow(2) 
      svert(2,3) = nlow(3) 
c     Y COEFF
      svert(3,1) = nlow(1) + 1
      svert(3,2) = nlow(2) + 1
      svert(3,3) = nlow(3) 
      svert(4,1) = nlow(1) + 1
      svert(4,2) = nlow(2)
      svert(4,3) = nlow(3)
c     Z COEFF
      svert(5,1) = nlow(1) + 1
      svert(5,2) = nlow(2) + 1
      svert(5,3) = nlow(3) + 1
      svert(6,1) = nlow(1) + 1
      svert(6,2) = nlow(2) + 1
      svert(6,3) = nlow(3)
      call dosimplex(parentptr,childptr,childn,nlow,ddx,svert,gfnum)
      return
      end subroutine simplex6

c **********************************************************************
      subroutine dosimplex(parentptr,childptr,childn,nlow,ddx,sv,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),sv(6,3),gfnum
      real*8 ddx(3),delx(3),coeff(3),constcoeff,tally
      integer i,j,gf,tl

      do i = 1,3
         delx(i) = parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(nlow(i)+1) -
     &          parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(nlow(i))
      enddo
      constcoeff = 
     &    parentptr%gfcn(gfnum)%tlev(1)%data(nlow(1),nlow(2),nlow(3))
      coeff(1) = 
     &    parentptr%gfcn(gfnum)%tlev(1)%data(sv(1,1),sv(1,2),sv(1,3)) -
     &    parentptr%gfcn(gfnum)%tlev(1)%data(sv(2,1),sv(2,2),sv(2,3))
      coeff(2) = 
     &    parentptr%gfcn(gfnum)%tlev(1)%data(sv(3,1),sv(3,2),sv(3,3)) -
     &    parentptr%gfcn(gfnum)%tlev(1)%data(sv(4,1),sv(4,2),sv(4,3))
      coeff(3) = 
     &    parentptr%gfcn(gfnum)%tlev(1)%data(sv(5,1),sv(5,2),sv(5,3)) -
     &    parentptr%gfcn(gfnum)%tlev(1)%data(sv(6,1),sv(6,2),sv(6,3))
      tally = constcoeff
      do j=1,3
         coeff(j) = ddx(j) * coeff(j) / delx(j)
         tally = tally + coeff(j)
      enddo
c     PUT FINAL INTERPOLATED RESULT IN CHILD TARGET
      childptr%gfcn(gfnum)%tlev(1)%data(childn(1),childn(2),childn(3)) =
     &      tally
      return
      end subroutine dosimplex

c************************************************************************

      subroutine bracket(parentptr,childptr,childn,nlow,gfnum)

      type(GF3TL), pointer :: parentptr,childptr
      integer childn(3),nlow(3),nmid,nhigh,gfnum,i,childmax,parentmax
      real*8  targetcoord,highcoord,lowcoord,midcoord

c     CHECK FOR PROPER NESTING
      do i=1,3
         if(childptr%gfcn(gfnum)%tlev(1)%c(i)%coord(1) .lt. 
     &     parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(1) ) then
            write(6,*) 'error:bracket:child coordinates out of bounds'
            stop
         endif
         childmax = childptr%gfcn(gfnum)%tlev(1)%ddim(i)
         parentmax = parentptr%gfcn(gfnum)%tlev(1)%ddim(i)
         if(childptr%gfcn(gfnum)%tlev(1)%c(i)%coord(childmax) .gt. 
     &      parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(parentmax) ) then
            write(6,*) 'error:bracket:child coordinates out of bounds'
            stop
         endif
      enddo
c     BRACKET THE TARGET COORDINATE
      do i = 1,3
         nlow(i) = 1
         nhigh = parentptr%gfcn(gfnum)%tlev(1)%ddim(i)
         lowcoord = parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(1)
         highcoord = parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(nhigh)
         targetcoord = childptr%gfcn(gfnum)%tlev(1)%c(i)%coord(childn(i))
         do while (nhigh - nlow(i) .gt. 1)
            nmid = (nhigh + nlow(i))/2
            midcoord = parentptr%gfcn(gfnum)%tlev(1)%c(i)%coord(nmid)
            if(midcoord .gt. targetcoord) then
               nhigh = nmid
               highcoord = midcoord
            else
               nlow(i) = nmid
               lowcoord = midcoord
            endif
         enddo
      enddo
       
      return     
      end subroutine bracket
      
      end module adaptive










