!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              SMOOTH_PRESSURES    Jun. 2007                                   |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine smooth_pressures (osolve,ov,params)

use definitions

implicit none

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

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|
type (octreesolve) osolve
type (octreev) ov
type (parameters) params

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer ip,jp,i,ii,j,k,ie,nleaves,iv
integer ierr,iproc,nproc,nb
integer leaf,level,loc
integer, dimension(:),allocatable :: levs
integer, dimension(:), allocatable :: nb_nb
double precision leaf_size,alpha,eps,vol,coords(4),kh_ip
double precision :: x0,y0,z0,dxyz,rhoo,sumpress,s,kh,dist_ij,x,y,z
double precision, dimension(:),allocatable :: xc,yc,zc,xfield,yfield,zfield
double precision, dimension(:),allocatable :: field,newfield 
double precision, dimension(:),allocatable :: rho,correc 
double precision, dimension(:),allocatable :: countnode
character*72  :: shift

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

INCLUDE 'mpif.h'

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

shift=' '

select case (params%smoothing_type)

case(0)
!==============================================================================
!=====[No smoothing]===========================================================
!==============================================================================

osolve%spressure=osolve%pressure

case(1)
!==============================================================================
!=====[Center -> Nodes -> Center]==============================================
!==============================================================================

allocate (countnode(osolve%nnode)) 
allocate (field(osolve%nnode))

field=0
countnode=0

! A few comments for reference by dwhipp - 12/09 
! This loop goes through each element, loops through the associated nodes and 
! adds the elemental pressure to field (the nodal sum). countnode tracks the 
! number of times a pressure value is contributed 
do i=1,osolve%nleaves
   do j=1,8
      k=osolve%icon(j,i)
      field(k)           = field(k)      + osolve%pressure(i)     
      countnode(k)       = countnode(k)  + 1 
   end do
end do

! This loop goes through each node and divides the field value from above by 
! countnode, to yield a nodal value that averages the elemental pressures from 
! each element that shares that node 
do i=1,osolve%nnode
   if (countnode(i).gt.0.d0) then
      field(i)  = field(i) / countnode(i)
   end if
end do

! This stores the averaged nodal pressures
ov%temporary_nodal_pressure=field

do ie=1,osolve%nleaves
   osolve%spressure(ie)=sum(field(osolve%icon(1:8,ie)))/8.d0
enddo

deallocate(countnode)
deallocate(field)

if (iproc.eq.0 .and. params%debug>=1) then
   write(*,'(a,2e13.5)') shift//'smooth pressures :',minval(osolve%spressure), &
                                                     maxval(osolve%spressure)
end if

case(2)
!==============================================================================
!=====[Center -> Nodes -> Center, weighted by element volumes]=================
!==============================================================================

allocate (countnode(osolve%nnode)) 
allocate (field(osolve%nnode))

field=0
countnode=0

do i=1,osolve%nleaves
   vol=(abs( osolve%x(osolve%icon(1,i))-osolve%x(osolve%icon(2,i)) ) )**3
   if (vol.le.0.d0) print *,'pb with vol in pressures'
   do j=1,8
      k=osolve%icon(j,i)
      field(k)           = field(k)      +osolve%pressure(i)       *vol
      countnode(k)       = countnode(k)  +                         vol
   end do
end do

do i=1,osolve%nnode
   if (countnode(i).gt.0.d0) then
      field(i)  = field(i) / countnode(i)
   end if
end do

do ie=1,osolve%nleaves
   osolve%spressure(ie)=sum(field(osolve%icon(1:8,ie)))/8.d0
enddo

deallocate(countnode)
deallocate(field)

if (iproc.eq.0 .and. params%debug>=1) then
   write(*,'(a,2e13.5)') shift//'smooth pressures :',minval(osolve%spressure), &
                                                     maxval(osolve%spressure)
end if

!==============================================================================
!=====[error]==================================================================
!==============================================================================

case default
   call stop_run('pb with smoothing_type in spressures')

end select


return 
end

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