!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              COMPUTE_DIVERGENCE    Oct. 2007                                 |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine compute_divergence (params,osolve,ov,vo,istep,iter)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine calculates the divergence inside each element 
! by using the velocity information at the nodes.

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
use gauss
!use mpi

implicit none

include 'mpif.h'

type(parameters) params
type (octreev) ov
type (octreesolve) osolve
type (void) vo
integer istep
integer iter


!------------------------------------------------------------------------------|
!(((((((((((((((( 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,dimension(:),allocatable :: divergence,divergencetemp
double precision :: div,volume
integer err,ierr,nproc,iproc,k,i,iint,nint,mpe

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

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

nint=8

mpe=params%mpe

allocate (divergence(osolve%nleaves),stat=err)
if (err.ne.0) call stop_run ('Error alloc divergence in compute_divergence$')
allocate (divergencetemp(osolve%nleaves),stat=err)
if (err.ne.0) call stop_run ('Error alloc divergencetemp in compute_divergence$')
allocate (x(mpe),y(mpe),z(mpe),stat=err)
if (err.ne.0) call stop_run ('Error alloc x,y,z in compute_divergence$')
allocate (vx(mpe),vy(mpe),vz(mpe),stat=err)
if (err.ne.0) call stop_run ('Error alloc vx,vy,vz in compute_divergence$')
allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=err)
if (err.ne.0) call stop_run ('Error alloc dhdx,dhdy,dhdz in compute_divergence$')


divergence=0.d0
divergencetemp=0.d0

do i=1+iproc,osolve%nleaves,nproc
   div=0.d0
   if (vo%leaf(i).eq.0) then 
      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
      do iint=1,nint
         call compute_dhdx_dhdy_dhdz (mpe,rr(iint),ss(iint),tt(iint),x,y,z,dhdx,dhdy,dhdz,volume)
         do k=1,mpe
            div=div+(dhdx(k)*vx(k)+dhdy(k)*vy(k)+dhdz(k)*vz(k))*volume*ww(iint)
         enddo
      enddo
   end if
   divergencetemp(i)=div
end do

call mpi_reduce (divergencetemp, divergence, osolve%nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)

if (iproc.eq.0) write (8,*) 'divergence ',minval(divergence),maxval(divergence)

call DoRuRe_div_stats (params%doDoRuRe,istep,iter,osolve%nleaves,divergence)



deallocate (divergence)
deallocate (divergencetemp)
deallocate (x,y,z)
deallocate (vx,vy,vz)
deallocate (dhdx,dhdy,dhdz)

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|