!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              UPDATE_STRAIN_ON_CLOUD    Feb. 2007                             |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine update_cloud_fields (cl,ov,osolve,vo,params)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine is used to update the total strain value stored
! on the particles of the 3D cloud from the strainrate obtained
! by interpolation and differentiation of the velocity field
! from the octreev structure

! cl is the cloud object
! ov is the velocity octree
! dt is time step

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

use definitions
use invariants

implicit none
      
type (octreev) ov
type (cloud) cl
type (octreesolve) osolve
type (void) vo
type (parameters) params

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

integer iproc,nproc,ierr,i,k,ic,leaf,level,loc,err,mpe,j,inode
double precision exx,eyy,ezz,exy,eyz,ezx,e2d
double precision jcb(3,3),jcbi(3,3)
double precision dhdr(8),dhds(8),dhdt(8),x(8),y(8),z(8),vx(8),vy(8),vz(8)
double precision dhdx(8),dhdy(8),dhdz(8)
double precision x0,y0,z0,dxyz,r,s,t,volume
double precision avrg_x,avrg_y,avrg_z
logical is_in_void(8)
integer nv
double precision,dimension(:),allocatable::strain,lsf0,temp,press,lsf
double precision,dimension(:),allocatable::nodal_pressure
integer, dimension(:),allocatable :: countnode
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

INCLUDE 'mpif.h'

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

mpe=8

allocate (strain(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%strain in update_cloud_fields$')

strain=0.d0

do i=1+iproc,cl%np,nproc
   call octree_find_leaf (ov%octree,ov%noctree, &
                          cl%x(i),cl%y(i),cl%z(i), &
                          leaf,level,loc,x0,y0,z0,dxyz)

!if(ov%whole_leaf_in_fluid(leaf)) then !MODIF NAZE PR RITSKE

   is_in_void=.false.
   do k=1,mpe
      inode=ov%icon(k,leaf)
      x(k)=ov%x(inode)
      y(k)=ov%y(inode)
      z(k)=ov%z(inode)
      vx(k)=ov%unode(inode)
      vy(k)=ov%vnode(inode)
      vz(k)=ov%wnode(inode)
      if ((vo%node(inode))==1) is_in_void(k)=.true.
   end do

   if (count(is_in_void)/=0 .and. count(is_in_void)<8) then
      nv=count(is_in_void)
      avrg_x=sum(vx,.not.is_in_void)/dble(8-nv)
      avrg_y=sum(vy,.not.is_in_void)/dble(8-nv)
      avrg_z=sum(vz,.not.is_in_void)/dble(8-nv)
      do k=1,mpe
         if (is_in_void(k)) then
            vx(k)=avrg_x
            vy(k)=avrg_y
            vz(k)=avrg_z
         end if
      end do
   end if



!   do k=1,mpe
!      x(k)=ov%x(ov%icon(k,leaf))
!      y(k)=ov%y(ov%icon(k,leaf))
!      z(k)=ov%z(ov%icon(k,leaf))
!      vx(k)=ov%unode(ov%icon(k,leaf))
!      vy(k)=ov%vnode(ov%icon(k,leaf))
!      vz(k)=ov%wnode(ov%icon(k,leaf))
!   end do

   r=(cl%x(i)-x0)/dxyz*2.d0-1.d0
   s=(cl%y(i)-y0)/dxyz*2.d0-1.d0
   t=(cl%z(i)-z0)/dxyz*2.d0-1.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

   e2d=sqrt(second_invariant(exx,eyy,ezz,exy,eyz,ezx))
   strain(i)=cl%strain(i)+e2d*params%dt

!end if ! END MODIF

enddo


cl%strain=0.d0

call mpi_allreduce (strain,cl%strain,cl%np,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

if (iproc.eq.0.d0) then
   write (8,*) 'Strain on cloud ', minval(cl%strain),maxval(cl%strain)
end if

deallocate (strain)


allocate (lsf(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsf in update_cloud_fields$')
lsf=0.d0
if (osolve%nlsf.ne.0) lsf=osolve%lsf(1:osolve%nnode,1)

!allocate(nodal_pressure(ov%nnode))
!allocate(countnode(ov%nnode))
!nodal_pressure=0.d0
!do i=1,ov%nleaves
!   do j=1,8
!      k=osolve%icon(j,i)
!      nodal_pressure(k)  = nodal_pressure(k) + ov%pressure(i)  
!      countnode(k)       = countnode(k)      + 1 
!   end do
!end do
!do i=1,ov%nnode
!   if (countnode(i).gt.0.d0) then
!      nodal_pressure(i)  = nodal_pressure(i) / countnode(i)
!   end if
!end do

allocate (lsf0(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsf0 in update_cloud_fields$')
allocate (temp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc temp in update_cloud_fields$')
allocate (press(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc press in update_cloud_fields$')

lsf0=0.d0
temp=0.d0
press=0.d0

do k=1+iproc,cl%np,nproc
   call octree_interpolate_three (3,ov%octree,ov%noctree,ov%icon, &
                                 ov%nleaves,ov%nnode, &
                                 cl%x(k),cl%y(k),cl%z(k), &
                                 lsf,lsf0(k), &
                                 ov%temp,temp(k), &
                                 ov%temporary_nodal_pressure,press(k))
!   call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,lsf                        ,ov%x,ov%y,ov%z,ov%nnode,cl%x(k),cl%y(k),cl%z(k),lsf0(k))
!   call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%temp                    ,ov%x,ov%y,ov%z,ov%nnode,cl%x(k),cl%y(k),cl%z(k),temp(k))
!   call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%temporary_nodal_pressure,ov%x,ov%y,ov%z,ov%nnode,cl%x(k),cl%y(k),cl%z(k),press(k))
enddo

cl%lsf0=0.d0
cl%temp=0.d0
cl%press=0.d0

call mpi_allreduce (lsf0, cl%lsf0, cl%np, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (temp, cl%temp, cl%np, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (press,cl%press,cl%np, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

deallocate (lsf0)
deallocate (temp)
deallocate (press)
deallocate (lsf)
!deallocate (nodal_pressure)
!deallocate (countnode)

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