!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              MAKE_PRESSURE    Apr. 2007                                      |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine make_pressure (params,icon,xg,yg,zg,viscosity0,penalty,expon,unode, &
                         vnode,wnode,temp,pressure,nnode,r0,s0,t0,rst,         &
                         eviscosity)

use definitions
use invariants
!use mpi

implicit none

include 'mpif.h'

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine is a stripped down version of make_matrix but is used in the
! calculation of the pressure

! icon is connectivity array for the current element
! xg,yg,zg are the global coordinate arrays of length nnode
! viscosity0 is viscosity
! penalty is incompressibility penalty factor
! expon is nonlinear viscosity exponent
! unode,vnode,wnode is the velocity array (obtained from previous time step
! or at least containing the proper velocity at the fixed dofs)
! temp and pressure anre temperature and pressure
! nnode is number of nodes
! r0,s0,t0 are the bottom left back corner of the part of the element
!   we are now integrating (in local r,s,t coordinates)
! rst is the size of the part of the element we are integrating

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

type(parameters) params
integer icon(params%mpe)
double precision xg(nnode),yg(nnode),zg(nnode)
double precision viscosity0
double precision penalty
double precision expon
double precision unode(nnode),vnode(nnode),wnode(nnode)
double precision temp(nnode)
double precision pressure
integer nnode
double precision r0,s0,t0,rst
double precision eviscosity

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
!integer nint
integer iproc,nproc,ierr,ic,k,iint,err
double precision viscosity
double precision exx,eyy,ezz,exy,eyz,ezx,E2d,E3d
double precision f1(3,3),defgrad(3,3)
double precision detf1,detf,e2dref,viscomean
double precision velox,veloy,veloz,tau,uvwnorm,xmin,xmax,ymin,ymax
double precision zmin,zmax,dx,dy,dz,alpha,div
double precision eps,r,s,t,w,volume
double precision jcb(3,3),jcbi(3,3),jcbp(3,3),jcbip(3,3)
double precision,dimension(:),allocatable :: x,y,z
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
logical is_plastic_dummy,flag_vrm_pb

double precision :: rr,ss,tt,ww

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

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

eps=tiny(eps)

allocate (x(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc x in make_pressure$')
allocate (y(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc y in make_pressure$')
allocate (z(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc z in make_pressure$')
allocate (dhdx(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdx in make_pressure$')
allocate (dhdy(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdy in make_pressure$')
allocate (dhdz(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdz in make_pressure$')

rr=0.d0
ss=0.d0
tt=0.d0
ww=8.d0

do k=1,params%mpe
   x(k)=xg(icon(k))
   y(k)=yg(icon(k))
   z(k)=zg(icon(k))*params%vex
enddo

! integrate over nint integration points
! here nint has been chosen to 1 because incompresibility
! is imposed in the FE matrix with a single integration point

   r=r0+rst*(rr+1.d0)/2.d0
   s=s0+rst*(ss+1.d0)/2.d0
   t=t0+rst*(tt+1.d0)/2.d0
   w=ww

   ! added by Jean 13/2/08
   r=rr
   s=ss
   t=tt
   w=ww
   !end addition

   call compute_dhdx_dhdy_dhdz (params%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,params%mpe
      ic=icon(k)
      exx=exx+dhdx(k)*unode(ic)
      eyy=eyy+dhdy(k)*vnode(ic)
      ezz=ezz+dhdz(k)*wnode(ic)
   enddo

! Change by Jean 13/2/08
! divergence should not be scaled by volume or weight
!   div=(exx+eyy+ezz)*volume*w
   div=exx+eyy+ezz
! end of change

if (params%bulkvisc) then
  pressure=div*penalty
else
  pressure=div*penalty*eviscosity
endif

deallocate (x)
deallocate (y)
deallocate (z)
deallocate (dhdx)
deallocate (dhdy)
deallocate (dhdz)

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