!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              COMPUTE_DHDX_DHDY_DHDZ    Nov. 2006                             |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine compute_dhdx_dhdy_dhdz (mpe,r,s,t,x,y,z,dhdx,dhdy,dhdz,volume)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

implicit none

integer mpe
double precision r,s,t
double precision x(mpe),y(mpe),z(mpe)
double precision dhdx(mpe),dhdy(mpe),dhdz(mpe)
double precision volume

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision,dimension(:),allocatable :: h,dhdr,dhds,dhdt
double precision jcb(3,3),jcbi(3,3)
double precision eps
integer k,err
integer nint

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


nint=8

eps=tiny(eps)

allocate (h(mpe)   ,stat=err) ; if (err.ne.0) call stop_run ('Error alloc h    in compute_dhdx_dhdy_dhdz$')
allocate (dhdr(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdr in compute_dhdx_dhdy_dhdz$')
allocate (dhds(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhds in compute_dhdx_dhdy_dhdz$')
allocate (dhdt(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdt in compute_dhdx_dhdy_dhdz$')
   
   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

   jcb=0.
   do k=1,mpe
      jcb(1,1)=jcb(1,1)+dhdr(k)*x(k)
      jcb(1,2)=jcb(1,2)+dhdr(k)*y(k)
      jcb(1,3)=jcb(1,3)+dhdr(k)*z(k)
      jcb(2,1)=jcb(2,1)+dhds(k)*x(k)
      jcb(2,2)=jcb(2,2)+dhds(k)*y(k)
      jcb(2,3)=jcb(2,3)+dhds(k)*z(k)
      jcb(3,1)=jcb(3,1)+dhdt(k)*x(k)
      jcb(3,2)=jcb(3,2)+dhdt(k)*y(k)
      jcb(3,3)=jcb(3,3)+dhdt(k)*z(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)

   if (volume.le.eps) then
      call stop_run ('Element bow-tied or collapsed in make_matrix$')
   endif

   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,mpe
      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


deallocate (h)
deallocate (dhdr)
deallocate (dhds)
deallocate (dhdt)


return

end

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