subroutine octree_interpolate_three(nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
                                    field1,f1,field2,f2,field3,f3)

! octree_interpolate_many specialized for 3
! optional parameters removed.
! Author: Douglas Guptill
! Date: 2009-07-19

! This function returns the value of three fields (fieldi) known at the nodes
! of an octree by trilinear interpolation

! nf is the number of fields being interpolate (must be 3)
! icon is the connectivity matrix
! nleaves is the number of leaves in the octree
! fieldi are the arrays of dimension nfield containing the fields
!     known at the nodes of the octree and to be interpolated
! x,y,z are the location of the point where the fields are to be interpolated
! fi are the resulting interpolated fields

implicit none

integer noctree,octree(noctree),nleaves,icon(8,nleaves)
integer nfield,nf
double precision field1(nfield),field2(nfield),field3(nfield)
double precision f1,f2,f3
double precision x,y,z,x0,y0,z0,dxyz,r,s,t,h(8),phi,xt,yt,zt
integer leaf,level,loc,k,iii,jjj,kkk,ii

! function modified by JEAN BRAUN on September 26 2005
! to correct for an error in the logics that led to an interpolation
! from an octree to another identical octree with differences in the
! interpolated function. The reason for this problem was related to
! bad faces or hanging nodes. Indeed, for a hanging node it was very likely
! that the leaf that was detected as the leaf in which the node resides
! was in fact a leaf where the node was a hanging node (ie not one of the
! 4 corner nodes). This meant that the interpolated value was not equal
! to the "constrained" value imposed by the linear constraint at the
! hanging node. To correct for this we first check if the node can
! be interpolated with r,s,t values that are equal to 1 or -1. If this is
! true than this value is chosen as this would correspond to a nodal value

xt=x
yt=y
zt=z

if (xt.lt.-1.e-11 .or. xt.gt.1.d0+1.d-11) return
if (yt.lt.-1.e-11 .or. yt.gt.1.d0+1.d-11) return
if (zt.lt.-1.e-11 .or. zt.gt.1.d0+1.d-11) return

if (x.lt.1.e-11) xt=1.e-11
if (x.gt.1.d0-1.d-11) xt=1.d0-1.d-11
if (y.lt.1.e-11) yt=1.e-11
if (y.gt.1.d0-1.d-11) yt=1.d0-1.d-11
if (z.lt.1.e-11) zt=1.e-11
if (z.gt.1.d0-1.d-11) zt=1.d0-1.d-11

if (nf.ne.3) write(*,*) 'error: octree_interpolate_three, nf = ', nf

do kkk=-1,1,2
   do jjj=-1,1,2
      do iii=-1,1,2

         xt=x+iii*1.d-10
         yt=y+jjj*1.d-10
         zt=z+kkk*1.d-10

         if  (xt*(xt-1.d0).ge.0d0 .or. &
              yt*(yt-1.d0).ge.0d0 .or. &
              zt*(zt-1.d0).ge.0d0) goto 111

         call octree_find_leaf (octree,noctree,xt,yt,zt,leaf,level,loc,x0,y0,z0,dxyz)

         r=(x-x0)/dxyz*2.d0-1.d0
         s=(y-y0)/dxyz*2.d0-1.d0
         t=(z-z0)/dxyz*2.d0-1.d0
         h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
         h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
         h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
         h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
         h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
         h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
         h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
         h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
         phi=0.d0
         do k=1,8
            phi=phi+h(k)*field1(icon(k,leaf))
         enddo
         f1=phi
         if (nf.eq.1) goto 222
         phi=0.d0
         do k=1,8
            phi=phi+h(k)*field2(icon(k,leaf))
         enddo
         f2=phi
         if (nf.eq.2) goto 222
         phi=0.d0
         do k=1,8
            phi=phi+h(k)*field3(icon(k,leaf))
         enddo
         f3=phi
         if (nf.eq.3) goto 222

222      continue
         !if (abs(r)-1.d0+abs(s)-1.d0+abs(t)-1.d0.lt.1.d-10) return
         if  (abs(abs(r)-1.d0).lt.1.d-10 .and. &
              abs(abs(s)-1.d0).lt.1.d-10 .and. &
              abs(abs(t)-1.d0).lt.1.d-10) return
         
111      continue
         
      enddo
   enddo
enddo

return
end