OctreeBitPlus.f90

                   abs(nodez(inode)-z)<1.d-10 ) then
                  f=field(inode)
                  return
               end if
            end do
         end if
      end do
   end do
end do
!==============


xt=x
yt=y
zt=z

do kkk=-1,1,2
   zt=z+kkk*1.d-10
   do jjj=-1,1,2
      yt=y+jjj*1.d-10
      do iii=-1,1,2
         xt=x+iii*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)*field(icon(k,leaf))
         enddo
         f=phi
         if (abs(r)-1.d0+abs(s)-1.d0+abs(t)-1.d0.lt.1.d-10) return
         111 continue
      enddo
   enddo
enddo

return
end














!===================================!
!=====[OCTREE_INTERPOLATE_MANY]=====!
!===================================!

subroutine octree_interpolate_many (nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
                                    field1,f1,field2,f2,field3,f3, &
                                    field4,f4,field5,f5,field6,f6, &
                                    field7,f7,field8,f8,field9,f9, &
                                    field10,f10,field11,f11,field12,f12, &
                                    field13,f13,field14,f14,field15,f15)
! This function returns the value of several fields (fieldi) known at the nodes
! of an octree by trilinear interpolation

! nf is the number of fields being interpolate (must be comprised between 1 and 15)
! 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

! Note that the number of arguments to this routine depends on the number of
! fields to be interpolated (nf). This is why some of the arguments are declared
! as optional

implicit none

optional :: field2,f2,field3,f3,field4,f4,field5,f5,field6,f6
optional :: field7,f7,field8,f8,field9,f9,field10,f10
optional :: field11,f11,field12,f12,field13,f13,field14,f14,field15,f15

integer noctree,octree(noctree),nleaves,icon(8,nleaves)
integer nfield,nf
double precision field1(nfield),field2(nfield),field3(nfield),field4(nfield), &
                 field5(nfield),field6(nfield),field7(nfield),field8(nfield)
double precision field9(nfield),field10(nfield),field11(nfield),field12(nfield), &
                 field13(nfield),field14(nfield),field15(nfield)
double precision f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15
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 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
  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
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field4(icon(k,leaf))
    enddo
  f4=phi
if (nf.eq.4) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field5(icon(k,leaf))
    enddo
  f5=phi
if (nf.eq.5) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field6(icon(k,leaf))
    enddo
  f6=phi
if (nf.eq.6) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field7(icon(k,leaf))
    enddo
  f7=phi
if (nf.eq.7) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field8(icon(k,leaf))
    enddo
  f8=phi
if (nf.eq.8) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field9(icon(k,leaf))
    enddo
  f9=phi
if (nf.eq.9) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field10(icon(k,leaf))
    enddo
  f10=phi
if (nf.eq.10) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field11(icon(k,leaf))
    enddo
  f11=phi
if (nf.eq.11) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field12(icon(k,leaf))
    enddo
  f12=phi
if (nf.eq.12) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field13(icon(k,leaf))
    enddo
  f13=phi
if (nf.eq.13) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field14(icon(k,leaf))
    enddo
  f14=phi
if (nf.eq.14) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field15(icon(k,leaf))
    enddo
  f15=phi
if (nf.gt.15) stop 'too many field'

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

!==============================================!
!=====[OCTREE_INTERPOLATE_MANY_DERIVATIVE]=====!
!==============================================!

subroutine octree_interpolate_many_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, &
                                    field4,f4,df4x,df4y,df4z, &
                                    field5,f5,df5x,df5y,df5z, &
                                    field6,f6,df6x,df6y,df6z, &
                                    field7,f7,df7x,df7y,df7z, &
                                    field8,f8,df8x,df8y,df8z, &
                                    field9,f9,df9x,df9y,df9z, &
                                    field10,f10,df10x,df10y,df10z, &
                                    field11,f11,df11x,df11y,df11z, &
                                    field12,f12,df12x,df12y,df12z, &
                                    field13,f13,df13x,df13y,df13z, &
                                    field14,f14,df14x,df14y,df14z, &
                                    field15,f15,df15x,df15y,df15z)
! 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 comprised between 1 and 15)
! 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

! Note that the number of oarguments to this routine depends on the number of
! fields to be interpolated (nf). This is why some of the arguments are declared
! as optional

implicit none

optional :: field2,f2,field3,f3,field4,f4,field5,f5,field6,f6
optional :: field7,f7,field8,f8,field9,f9,field10,f10
optional :: field11,f11,field12,f12,field13,f13,field14,f14,field15,f15
optional :: df2x,df3x,df4x,df5x,df6x,df7x,df8x,df9x
optional :: df10x,df11x,df12x,df13x,df14x,df15x
optional :: df2y,df3y,df4y,df5y,df6y,df7y,df8y,df9y
optional :: df10y,df11y,df12y,df13y,df14y,df15y
optional :: df2z,df3z,df4z,df5z,df6z,df7z,df8z,df9z
optional :: df10z,df11z,df12z,df13z,df14z,df15z


integer noctree,octree(noctree),nleaves,icon(8,nleaves)
integer nfield,nf
double precision field1(nfield),field2(nfield),field3(nfield),field4(nfield), &
                 field5(nfield),field6(nfield),field7(nfield),field8(nfield)
double precision field9(nfield),field10(nfield),field11(nfield),field12(nfield), &
                 field13(nfield),field14(nfield),field15(nfield)
double precision f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15
double precision df1x,df2x,df3x,df4x,df5x,df6x,df7x,df8x,df9x
double precision df10x,df11x,df12x,df13x,df14x,df15x
double precision df1y,df2y,df3y,df4y,df5y,df6y,df7y,df8y,df9y
double precision df10y,df11y,df12y,df13y,df14y,df15y
double precision df1z,df2z,df3z,df4z,df5z,df6z,df7z,df8z,df9z
double precision df10z,df11z,df12z,df13z,df14z,df15z
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
if (nf.eq.3) 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)*field4(ic)
    phix=phix+dhdx(k)*field4(ic)
    phiy=phiy+dhdy(k)*field4(ic)
    phiz=phiz+dhdz(k)*field4(ic)
    enddo
  f4=phi
  df4x=phix
  df4y=phiy
  df4z=phiz
if (nf.eq.4) 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)*field5(ic)
    phix=phix+dhdx(k)*field5(ic)
    phiy=phiy+dhdy(k)*field5(ic)
    phiz=phiz+dhdz(k)*field5(ic)
    enddo
  f5=phi
  df5x=phix
  df5y=phiy
  df5z=phiz
if (nf.eq.5) 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)*field6(ic)
    phix=phix+dhdx(k)*field6(ic)
    phiy=phiy+dhdy(k)*field6(ic)
    phiz=phiz+dhdz(k)*field6(ic)
    enddo
  f6=phi
  df6x=phix
  df6y=phiy
  df6z=phiz
if (nf.eq.6) 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)*field7(ic)
    phix=phix+dhdx(k)*field7(ic)
    phiy=phiy+dhdy(k)*field7(ic)
    phiz=phiz+dhdz(k)*field7(ic)
    enddo
  f7=phi
  df7x=phix
  df7y=phiy
  df7z=phiz
if (nf.eq.7) 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)*field8(ic)
    phix=phix+dhdx(k)*field8(ic)
    phiy=phiy+dhdy(k)*field8(ic)
    phiz=phiz+dhdz(k)*field8(ic)
    enddo
  f8=phi
  df8x=phix
  df8y=phiy
  df8z=phiz
if (nf.eq.8) 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)*field9(ic)
    phix=phix+dhdx(k)*field9(ic)
    phiy=phiy+dhdy(k)*field9(ic)
    phiz=phiz+dhdz(k)*field9(ic)
    enddo
  f9=phi
  df9x=phix
  df9y=phiy
  df9z=phiz
if (nf.eq.9) 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)*field10(ic)
    phix=phix+dhdx(k)*field10(ic)
    phiy=phiy+dhdy(k)*field10(ic)
    phiz=phiz+dhdz(k)*field10(ic)
    enddo
  f10=phi
  df10x=phix
  df10y=phiy
  df10z=phiz
if (nf.eq.10) 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)*field11(ic)
    phix=phix+dhdx(k)*field11(ic)
    phiy=phiy+dhdy(k)*field11(ic)
    phiz=phiz+dhdz(k)*field11(ic)
    enddo
  f11=phi
  df11x=phix
  df11y=phiy
  df11z=phiz
if (nf.eq.11) 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)*field12(ic)
    phix=phix+dhdx(k)*field12(ic)
    phiy=phiy+dhdy(k)*field12(ic)
    phiz=phiz+dhdz(k)*field12(ic)
    enddo
  f12=phi
  df12x=phix
  df12y=phiy
  df12z=phiz
if (nf.eq.12) 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)*field13(ic)
    phix=phix+dhdx(k)*field13(ic)
    phiy=phiy+dhdy(k)*field13(ic)
    phiz=phiz+dhdz(k)*field13(ic)
    enddo
  f13=phi
  df13x=phix
  df13y=phiy
  df13z=phiz
if (nf.eq.13) 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)*field14(ic)
    phix=phix+dhdx(k)*field14(ic)
    phiy=phiy+dhdy(k)*field14(ic)
    phiz=phiz+dhdz(k)*field14(ic)
    enddo
  f14=phi
  df14x=phix
  df14y=phiy
  df14z=phiz
if (nf.eq.14) 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)*field15(ic)
    phix=phix+dhdx(k)*field15(ic)
    phiy=phiy+dhdy(k)*field15(ic)
    phiz=phiz+dhdz(k)*field15(ic)
    enddo
  f15=phi
  df15x=phix
  df15y=phiy
  df15z=phiz
if (nf.gt.15) stop 'too many field'

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


!-------------------------------------

Subroutine mrgrnk (XDONT, IRNGT, np)
!Subroutine mrgrnk (XDONT, IRNGT)
! __________________________________________________________
!   MRGRNK = Merge-sort ranking of an array
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! __________________________________________________________
      Integer XDONT(np),IRNGT(np)
!      Integer, Dimension (:), Intent (In)  :: XDONT
!      Integer, Dimension (:), Intent (Out) :: IRNGT
! __________________________________________________________
      Integer :: XVALA, XVALB
!
      Integer, Dimension(:),allocatable :: JWRKT
!      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
      Select Case (NVAL)
      Case (:0)
         Return
      Case (1)
         IRNGT (1) = 1
         Return
      Case Default
         Continue
      End Select

allocate (JWRKT(np))

!
!  Fill-in the index array, creating ordered couples
!
      Do IIND = 2, NVAL, 2
         If (XDONT(IIND-1) <= XDONT(IIND)) Then
            IRNGT (IIND-1) = IIND - 1
            IRNGT (IIND) = IIND
         Else
            IRNGT (IIND-1) = IIND
            IRNGT (IIND) = IIND - 1
         End If
      End Do
      If (Modulo(NVAL, 2) /= 0) Then
         IRNGT (NVAL) = NVAL
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      LMTNA = 2
      LMTNC = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (NVAL <= 2) Exit
!
!   Loop on merges of A and B into C
!
         Do IWRKD = 0, NVAL - 1, 4
            If ((IWRKD+4) > NVAL) Then
               If ((IWRKD+2) >= NVAL) Exit
!
!   1 2 3
!
               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
!
!   1 3 2
!
               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
                  IRNG2 = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNG2
!
!   3 1 2
!
               Else
                  IRNG1 = IRNGT (IWRKD+1)
                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNG1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
!
!   1 3 x x
!
            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   1 3 2 4
                  IRNGT (IWRKD+3) = IRNG2
               Else
!   1 3 4 2
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+4) = IRNG2
               End If
!
!   3 x x x
!
            Else
               IRNG1 = IRNGT (IWRKD+1)
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
                  IRNGT (IWRKD+2) = IRNG1
                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   3 1 2 4
                     IRNGT (IWRKD+3) = IRNG2
                  Else
!   3 1 4 2
                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                     IRNGT (IWRKD+4) = IRNG2
                  End If
               Else
!   3 4 1 2
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+3) = IRNG1
                  IRNGT (IWRKD+4) = IRNG2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (LMTNA >= NVAL) Exit
         IWRKF = 0
         LMTNC = 2 * LMTNC
!
!   Loop on merges of A and B into C
!
         Do
            IWRK = IWRKF
            IWRKD = IWRKF + 1
            JINDA = IWRKF + LMTNA
            IWRKF = IWRKF + LMTNC
            If (IWRKF >= NVAL) Then
               If (JINDA >= NVAL) Exit
               IWRKF = NVAL
            End If
            IINDA = 1
            IINDB = JINDA + 1
!
!   Shortcut for the case when the max of A is smaller
!   than the min of B. This line may be activated when the
!   initial set is already close to sorted.
!
!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
!
!  One steps in the C subset, that we build in the final rank array
!
!  Make a copy of the rank array for the merge iteration
!
            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
!
            XVALA = XDONT (JWRKT(IINDA))
            XVALB = XDONT (IRNGT(IINDB))
!
            Do
               IWRK = IWRK + 1
!
!  We still have unprocessed values in both A and B
!
               If (XVALA > XVALB) Then
                  IRNGT (IWRK) = IRNGT (IINDB)
                  IINDB = IINDB + 1
                  If (IINDB > IWRKF) Then
!  Only A still with unprocessed values
                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
                     Exit
                  End If
                  XVALB = XDONT (IRNGT(IINDB))
               Else
                  IRNGT (IWRK) = JWRKT (IINDA)
                  IINDA = IINDA + 1
                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
                  XVALA = XDONT (JWRKT(IINDA))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 2 * LMTNA
      End Do
!

deallocate (JWRKT)

      Return
!
End Subroutine mrgrnk