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

      subroutine interpolateH3(parentptr,childptr,drc)

      type(GF3TL), pointer :: parentptr,childptr
      real(SP), dimension(3)   :: drc   !dr of child     
      integer(I4B) i,j,k,gfnum,idx
      integer(I4B), 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) 
      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(I4B), 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%c(1)%n_coord
       ny=childptr%c(2)%n_coord
       nz=childptr%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'
!     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

!************************************************************************

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

      type(GF3TL), pointer :: parentptr,childptr
      integer(I4B) childnx,childny,childnz,childn(3),nlow(3)
      integer(I4B) i,gf,tl,gfnum
      real(DP) 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%c(i)%coord(childn(i))
!        COORDINATES RELATIVE TO (nxlow,nylow,nzlow)
         ddx(i) = p(i) - parentptr%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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

      do i = 1,3
         delx(i) = parentptr%c(i)%coord(nlow(i)+1) - &
               parentptr%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
!     PUT FINAL INTERPOLATED RESULT IN CHILD TARGET
      childptr%gfcn(gfnum)%tlev(1)%data(childn(1),childn(2),childn(3)) = &
           tally
      return
      end subroutine dosimplex

!************************************************************************

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

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

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