      !***********************************************************************!
      SUBROUTINE  gradconj(  u,  b,  p,  q,  r,  it_max,  prec,  itgc,  &
                             rnorm_k,  rnudt2,  ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d  )
      !***********************************************************************!
!
      USE numerics
!
      IMPLICIT  NONE
!
!--Interface
      INTERFACE
!
      include 'i_pmv.h'
!
      include 'i_saxpy.h'
!
      include 'i_prodscal.h'
!
      END INTERFACE
!
!--Declaration  des  arguments
!
!Solution  u
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(INOUT)  ::  u
!
!Second  membre
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(INOUT)  ::  b
!
!Vecteurs  intermediaires
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(INOUT)  ::  p
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(INOUT)  ::  q
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(INOUT)  ::  r
!
!Critere  de  convergence
      REAL(rp),  INTENT(IN)   ::  prec
!
!Scalaires  gradient  conjugue
      REAL(rp),  INTENT(OUT)  ::  rnorm_k
!
!Scalire Crank  Nicholson
      REAL(rp),  INTENT(IN)   ::  rnudt2
!
!Nombre  d'iterations  maximum  du gradient
      INTEGER,   INTENT(IN)   ::  it_max
!
!Nombre  iterations  en  temps
      INTEGER,   INTENT(OUT)  ::  itgc
!
!Tableau  contenant  les  voisins  du  domaine  courant
      INTEGER,   PARAMETER  ::  NB_VOISINS  =  6
      INTEGER,   DIMENSION(NB_VOISINS),  INTENT(IN)  ::  voisin
!
!Constantes MPI
      INTEGER,   INTENT(IN)  ::  comm3d
      INTEGER,   INTENT(IN)  ::  ip_xy,  ip_xz,  ip_yz
!
!--Declaration des variales locales
!
!Scalaires  gradient  conjugue
      REAL(rp)  ::  alpha_k
      REAL(rp)  ::  beta_k
!
!Compteurs
      INTEGER   ::  i,  j,  k
!
!Test  de  convergence
      LOGICAL   ::  convergence
!
!***********************************************************************!
!
!=================
! pmv   : r = Au
!=================
      CALL  pmv(  r,  u,  - rnudt2  )
!
!===================
! saxpy : r = b - Au
!===================
      CALL  saxpy(  r,  -one,  b,  one  )
!
!=======
! p <- r
!=======
      DO  k  =  sz,  ez
          DO  j  =  sy,  ey
              DO  i  =  sx,  ex
                  p(i,j,k) = r(i,j,k)
              END  DO
          END  DO
      END  DO
!
!==========
! || r ||^2
!==========
      CALL prodscal(  rnorm_k,  r,  r,  comm3d  )
!
      itgc = 0
      convergence = ( SQRT(rnorm_k) < prec )
!
      DO  WHILE  ( (.NOT. convergence)  .AND.  (itgc  <  it_max) )
!
          itgc  =  itgc  +  1
!
!===================
! calcul de \alpha_k
!===================
!
!Echange  des  points  aux  interfaces  pour  p
         CALL  communication(  p,  ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d  )
!
! pmv   :  q_k = A.p_k
         CALL  pmv(  q,  p,  - rnudt2  )
!
! (q_k,p_k)
         CALL  prodscal( alpha_k,  q,  p,  comm3d  )
!
         alpha_k = rnorm_k / alpha_k
         beta_k  = rnorm_k
!
!================
! calcul de u_k+1
!================
         CALL  saxpy(  u,  one,  p,    alpha_k  )
!
!================
! calcul de r_k+1
!================
         CALL  saxpy(  r,  one,  q,  - alpha_k  )
!
!==================
! calcul de \beta_k
!==================
!
! || r_k+1 ||^2
         CALL  prodscal(  rnorm_k,  r,  r,  comm3d  )
!
! \beta_k
         beta_k  = rnorm_k / beta_k
!
!================
! calcul de p_k+1
!================
         CALL  saxpy(  p,  beta_k,  r,  one  )
!
!Arret  du  programe  si  on  a  atteint  la  precision  souhaitee
         convergence  =  ( SQRT(rnorm_k)  <  prec )
!
      END  DO
!
      RETURN
      END  SUBROUTINE  gradconj
