!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              DO_LEAF_MEASUREMENTS    Nov. 2006                               |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
         
subroutine do_leaf_measurements (osolve,ov,params,istep,iter,iter_nl) 
      
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine calculates the second invariant of the strain rate, 
! and a value for the criterion used by the improve_octree subroutine,
! inside each element by using the velocity information at the nodes.
! mpe is number of nodes per element (8)
! icon is connectivity array 
! xnode,ynode,znode are the global coordinate arrays of length nnode
! unode,vnode,wnode are the velocity arrays of length nnode (obtained 
! from previous time step or at least containing the proper velocity 
! at the fixed dofs)
! nnode is number of nodes
! nleaves is the number of leaves/elements
! e2d is the second invariant of the strain rate
! crit is the criterion

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

use constants
use definitions
use gauss
use invariants
!use mpi

implicit none

include 'mpif.h'

type (octreev) ov
type (octreesolve) osolve
type (parameters) params
integer istep,iter,iter_nl

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

double precision,dimension(:),allocatable :: x,y,z,vx,vy,vz
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
double precision :: r,s,t,volume
double precision :: exx,eyy,ezz,exy,eyz,ezx,J2d,J3d
integer err,ierr,nproc,iproc,k,i,nleaves,mpe
double precision, dimension(:), allocatable :: e2dtemp,e3dtemp,crittemp,lodetemp
double precision, dimension(:), allocatable :: qtemp,dilatrtemp
character*72  :: shift

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

mpe=params%mpe

nleaves=osolve%nleaves

allocate (x(mpe),y(mpe),z(mpe),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc x,y,z in compute_e2d_crit$')
allocate (vx(mpe),vy(mpe),vz(mpe),stat=err)       ; if (err.ne.0) call stop_run ('Error alloc vx,vy,vz in compute_e2d_crit$')
allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdx,dhdy,dhdz in compute_e2d_crit$')
allocate (e2dtemp(nleaves),stat=err)              ; if (err.ne.0) call stop_run ('Error alloc e2dtemp in compute_e2d_crit$')
allocate (e3dtemp(nleaves),stat=err)              ; if (err.ne.0) call stop_run ('Error alloc e2dtemp in compute_e3d_crit$')
allocate (crittemp(nleaves),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc crittemp in compute_e2d_crit$')
allocate (lodetemp(nleaves),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc lodetemp in compute_e2d_crit$')
allocate (qtemp(nleaves),stat=err)                ; if (err.ne.0) call stop_run ('Error alloc qtemp in compute_e2d_crit$')
allocate (dilatrtemp(nleaves),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc dilatrtemp in compute_e2d_crit$')

e2dtemp=0.d0
e3dtemp=0.d0
crittemp=0.d0
lodetemp=0.d0
qtemp=0.d0
dilatrtemp=0.d0

osolve%e2d=0.d0
osolve%e3d=0.d0
osolve%crit=0.d0
osolve%lode=0.d0
osolve%q=0.d0
osolve%dilatr=0.d0

do i=1+iproc,nleaves,nproc

   do k=1,mpe
      x(k) =osolve%x(osolve%icon(k,i))
      y(k) =osolve%y(osolve%icon(k,i))
      z(k) =osolve%z(osolve%icon(k,i))
      vx(k)=ov%unode(osolve%icon(k,i))
      vy(k)=ov%vnode(osolve%icon(k,i))
      vz(k)=ov%wnode(osolve%icon(k,i))
   enddo

   r=0.d0
   s=0.d0
   t=0.d0

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

   exx=0.d0
   eyy=0.d0
   ezz=0.d0
   exy=0.d0
   eyz=0.d0
   ezx=0.d0
   do k=1,mpe
      exx=exx+dhdx(k)*vx(k)
      eyy=eyy+dhdy(k)*vy(k)
      ezz=ezz+dhdz(k)*vz(k)
      exy=exy+(dhdx(k)*vy(k)+dhdy(k)*vx(k))/2.d0
      eyz=eyz+(dhdy(k)*vz(k)+dhdz(k)*vy(k))/2.d0
      ezx=ezx+(dhdz(k)*vx(k)+dhdx(k)*vz(k))/2.d0
   enddo

   ! Dilatation rate
   dilatrtemp(i)=(exx+eyy+ezz)

   call deviatoric (params%invariants_2d,exx,eyy,ezz,exy,eyz,ezx)

!   if (.not.osolve%is_plastic(i)) then
!      qtemp(i)=0.d0
!   else
      qtemp(i)=2.d0*osolve%eviscosity(i)*sqrt(second_invariant(exx,eyy,ezz,exy,eyz,ezx))
!   end if

   J2d=second_invariant(exx,eyy,ezz,exy,eyz,ezx)
   J3d=third_invariant (exx,eyy,ezz,exy,eyz,ezx)

   e2dtemp(i)=sqrt(J2d)

   if (J2d/=0.d0) lodetemp(i)=lode_angle (J2d,J3d)/pi

   if (J3d<0.d0) then
      e3dtemp(i)=-abs(J3d)**(1./3.d0)
   else
      e3dtemp(i)=J3d**(1./3.d0)
   end if

   select case (params%refine_criterion)
      case(1) 
         crittemp(i)=e2dtemp(i)
      case(2) 
         crittemp(i)=sqrt(exx**2+eyy**2+ezz**2)
      case(3)
         crittemp(i)=sqrt(J2d)*abs(x(1)-x(2)) 
      case default
         crittemp(i)=0.d0
   end select

end do

call mpi_reduce (e2dtemp, osolve%e2d, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (e3dtemp, osolve%e3d, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (crittemp,osolve%crit,nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (lodetemp,osolve%lode,nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (qtemp, osolve%q, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (dilatrtemp, osolve%dilatr, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)

deallocate (x,y,z)
deallocate (vx,vy,vz)
deallocate (dhdx,dhdy,dhdz)
deallocate (crittemp)
deallocate (e2dtemp)
deallocate (e3dtemp)
deallocate (lodetemp) 
deallocate (qtemp) 
deallocate (dilatrtemp)

if (iproc.eq.0 .and. params%debug>=1 .and. params%compute_qpgram) then
   write(*,'(a,f11.5)')  shift//'min(q) =', minval(osolve%q)
   write(*,'(a,f11.5)')  shift//'max(q) =', maxval(osolve%q)
end if     

call DoRuRe_vel_stats  (params%doDoRuRe,istep,iter,iter_nl,ov%nnode,ov%unode,ov%vnode,ov%wnode)
call DoRuRe_leaf_stats (params%doDoRuRe,osolve%nleaves,osolve%e2d,osolve%e3d,osolve%pressure,osolve%q,osolve%dilatr,istep,iter,iter_nl)

return
end

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|