subroutine octree_interpolate_three_derivative &
     (nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
     field1,f1,df1x,df1y,df1z, &
     field2,f2,df2x,df2y,df2z, &
     field3,f3,df3x,df3y,df3z)

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

! This function returns the value of several fields (fieldi) known at the nodes
! of an octree by trilinear interpolation as well as their 3 spatial derivatives

! 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 df1x,df2x,df3x
double precision df1y,df2y,df3y
double precision df1z,df2z,df3z
double precision x,y,z,x0,y0,z0,dxyz,r,s,t,h(8),phi,xt,yt,zt
double precision phix,phiy,phiz
integer leaf,level,loc,k,iii,jjj,kkk,ii,jj,kk,ic,ijk
double precision dhdr(8),dhds(8),dhdt(8),xg(8),yg(8),zg(8),jcb(3,3),jcbi(3,3),volume
double precision dhdx(8),dhdy(8),dhdz(8)

! 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 loeaf in which the node resides
! was in fact a leave 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

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
dhdr(1)=-(1.d0-s)*(1.d0-t)/8.d0
dhdr(2)=(1.d0-s)*(1.d0-t)/8.d0
dhdr(3)=-(1.d0+s)*(1.d0-t)/8.d0
dhdr(4)=(1.d0+s)*(1.d0-t)/8.d0
dhdr(5)=-(1.d0-s)*(1.d0+t)/8.d0
dhdr(6)=(1.d0-s)*(1.d0+t)/8.d0
dhdr(7)=-(1.d0+s)*(1.d0+t)/8.d0
dhdr(8)=(1.d0+s)*(1.d0+t)/8.d0
dhds(1)=-(1.d0-r)*(1.d0-t)/8.d0
dhds(2)=-(1.d0+r)*(1.d0-t)/8.d0
dhds(3)=(1.d0-r)*(1.d0-t)/8.d0
dhds(4)=(1.d0+r)*(1.d0-t)/8.d0
dhds(5)=-(1.d0-r)*(1.d0+t)/8.d0
dhds(6)=-(1.d0+r)*(1.d0+t)/8.d0
dhds(7)=(1.d0-r)*(1.d0+t)/8.d0
dhds(8)=(1.d0+r)*(1.d0+t)/8.d0
dhdt(1)=-(1.d0-r)*(1.d0-s)/8.d0
dhdt(2)=-(1.d0+r)*(1.d0-s)/8.d0
dhdt(3)=-(1.d0-r)*(1.d0+s)/8.d0
dhdt(4)=-(1.d0+r)*(1.d0+s)/8.d0
dhdt(5)=(1.d0-r)*(1.d0-s)/8.d0
dhdt(6)=(1.d0+r)*(1.d0-s)/8.d0
dhdt(7)=(1.d0-r)*(1.d0+s)/8.d0
dhdt(8)=(1.d0+r)*(1.d0+s)/8.d0
ijk=0
do kk=0,1
do jj=0,1
do ii=0,1
ijk=ijk+1
xg(ijk)=x0+ii*dxyz
yg(ijk)=y0+jj*dxyz
zg(ijk)=z0+kk*dxyz
enddo
enddo
enddo
jcb=0.
do k=1,8
jcb(1,1)=jcb(1,1)+dhdr(k)*xg(k)
jcb(1,2)=jcb(1,2)+dhdr(k)*yg(k)
jcb(1,3)=jcb(1,3)+dhdr(k)*zg(k)
jcb(2,1)=jcb(2,1)+dhds(k)*xg(k)
jcb(2,2)=jcb(2,2)+dhds(k)*yg(k)
jcb(2,3)=jcb(2,3)+dhds(k)*zg(k)
jcb(3,1)=jcb(3,1)+dhdt(k)*xg(k)
jcb(3,2)=jcb(3,2)+dhdt(k)*yg(k)
jcb(3,3)=jcb(3,3)+dhdt(k)*zg(k)
enddo

volume=jcb(1,1)*jcb(2,2)*jcb(3,3)+jcb(1,2)*jcb(2,3)*jcb(3,1) &
+jcb(2,1)*jcb(3,2)*jcb(1,3) &
-jcb(1,3)*jcb(2,2)*jcb(3,1)-jcb(1,2)*jcb(2,1)*jcb(3,3) &
-jcb(2,3)*jcb(3,2)*jcb(1,1)

jcbi(1,1)=(jcb(2,2)*jcb(3,3)-jcb(2,3)*jcb(3,2))/volume
jcbi(2,1)=(jcb(2,3)*jcb(3,1)-jcb(2,1)*jcb(3,3))/volume
jcbi(3,1)=(jcb(2,1)*jcb(3,2)-jcb(2,2)*jcb(3,1))/volume
jcbi(1,2)=(jcb(1,3)*jcb(3,2)-jcb(1,2)*jcb(3,3))/volume
jcbi(2,2)=(jcb(1,1)*jcb(3,3)-jcb(1,3)*jcb(3,1))/volume
jcbi(3,2)=(jcb(1,2)*jcb(3,1)-jcb(1,1)*jcb(3,2))/volume
jcbi(1,3)=(jcb(1,2)*jcb(2,3)-jcb(1,3)*jcb(2,2))/volume
jcbi(2,3)=(jcb(1,3)*jcb(2,1)-jcb(1,1)*jcb(2,3))/volume
jcbi(3,3)=(jcb(1,1)*jcb(2,2)-jcb(1,2)*jcb(2,1))/volume

do k=1,8
dhdx(k)=jcbi(1,1)*dhdr(k)+jcbi(1,2)*dhds(k)+jcbi(1,3)*dhdt(k)
dhdy(k)=jcbi(2,1)*dhdr(k)+jcbi(2,2)*dhds(k)+jcbi(2,3)*dhdt(k)
dhdz(k)=jcbi(3,1)*dhdr(k)+jcbi(3,2)*dhds(k)+jcbi(3,3)*dhdt(k)
enddo

  phi=0.d0
  phix=0.d0
  phiy=0.d0
  phiz=0.d0
    do k=1,8
    ic=icon(k,leaf)
    phi=phi+h(k)*field1(ic)
    phix=phix+dhdx(k)*field1(ic)
    phiy=phiy+dhdy(k)*field1(ic)
    phiz=phiz+dhdz(k)*field1(ic)
    enddo
  f1=phi
  df1x=phix
  df1y=phiy
  df1z=phiz
if (nf.eq.1) goto 222
  phi=0.d0
  phix=0.d0
  phiy=0.d0
  phiz=0.d0
    do k=1,8
    ic=icon(k,leaf)
    phi=phi+h(k)*field2(ic)
    phix=phix+dhdx(k)*field2(ic)
    phiy=phiy+dhdy(k)*field2(ic)
    phiz=phiz+dhdz(k)*field2(ic)
    enddo
  f2=phi
  df2x=phix
  df2y=phiy
  df2z=phiz
if (nf.eq.2) goto 222
  phi=0.d0
  phix=0.d0
  phiy=0.d0
  phiz=0.d0
    do k=1,8
    ic=icon(k,leaf)
    phi=phi+h(k)*field3(ic)
    phix=phix+dhdx(k)*field3(ic)
    phiy=phiy+dhdy(k)*field3(ic)
    phiz=phiz+dhdz(k)*field3(ic)
    enddo
  f3=phi
  df3x=phix
  df3y=phiy
  df3z=phiz


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

