
      module myinterpolate
      use bh02


      contains 

      subroutine interpolateH3(parentptr,childptr)

      type(GF3TL), pointer :: parentptr,childptr
      integer i,j,k,gfnum
      integer, parameter :: gfcopnum = 35

c     THESE ARE THE GRID FUNCTIONS THAT NEED TO BE INTERPOLATED

      integer, dimension(35) ::
     & 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/)

      do gfnum = 1,gfcopnum
      do i = 1,childptr%gfcn(gfnum)%tlev(1)%ddim(1)
      do j = 1,childptr%gfcn(gfnum)%tlev(1)%ddim(2)
      do k = 1,childptr%gfcn(gfnum)%tlev(1)%ddim(3)
         call interpoint(parentptr,childptr,i,j,k,gfnum)
      enddo
      enddo
      enddo
      enddo
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)

      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)

      do i=1,3
c        TARGET COORDINATES
         p(i) = childptr%gfcn(gfnum)%tlev(1)%c(i)%coord(childn(i))
c        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
c           SIMPLEX1 : dx <= dy <= dz
            call simplex1(parentptr,childptr,childn,nlow,ddx,gfnum)
         else
            if(ddx(3) .le. ddx(1)) then
c              SIMPLEX5 : dz <= dx <= dy
               call simplex5(parentptr,childptr,childn,nlow,ddx,gfnum)
            else
c              SIMPLEX2 : dx <= dz <= dy
               call simplex2(parentptr,childptr,childn,nlow,ddx,gfnum)
            endif
         endif
      else
         if(ddx(3) .le. ddx(2)) then
c           SIMPLEX6 : dz <= dy <= dx
            call simplex6(parentptr,childptr,childn,nlow,ddx,gfnum)
         else
            if(ddx(3) .le. ddx(1)) then
c              SIMPLEX4: dy <= dz <= dx
               call simplex4(parentptr,childptr,childn,nlow,ddx,gfnum)
            else
c              SIMPLEX3: dy <= dx <= dz
               call simplex3(parentptr,childptr,childn,nlow,ddx,gfnum)
            endif
         endif
      endif
      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 myinterpolate
