      !***********************************************************************!
      SUBROUTINE  gradconj(  u,  b,  p,  q,  r,  voisin,  &
                             ptabN,  ptabE,  ptabS,  ptabW,  ptabAV,  ptabAR,  &
                             it_max,  prec,  itgc,  rnorm_k,  rnudt2 )
      !***********************************************************************!
!
      USE numerics
!
      IMPLICIT  NONE
!
      INCLUDE 'comthd.h'
!
!--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
!
!Buffers locaux de donnees
      REAL(rp),  DIMENSION(sy:ey,sz:ez),  INTENT(OUT)  ::  ptabE
      REAL(rp),  DIMENSION(sy:ey,sz:ez),  INTENT(OUT)  ::  ptabW
      REAL(rp),  DIMENSION(sx:ex,sz:ez),  INTENT(OUT)  ::  ptabS
      REAL(rp),  DIMENSION(sx:ex,sz:ez),  INTENT(OUT)  ::  ptabN
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1),  INTENT(OUT)  ::  ptabAV
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1),  INTENT(OUT)  ::  ptabAR
!
!Tableau  contenant  les  voisins  du  domaine  courant
      INTEGER,  PARAMETER  ::  NB_VOISINS  =  6
      INTEGER,  PARAMETER  ::  N=1,  E=2,  S=3,  W=4,  AV=5,  AR=6
      INTEGER,  DIMENSION(NB_VOISINS),  INTENT(IN)  ::  voisin
!
!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
!
!--Declaration des variales locales
!
!Scalaires  gradient  conjugue
      REAL(rp)  ::  alpha_k = zero
      REAL(rp)  ::  beta_k  = zero
!
!Compteurs
      INTEGER   ::  i,  j,  k
!
!Nombre  iterations  en  temps
      INTEGER   ::  it
!
!Test  de  convergence
      LOGICAL   ::  convergence
!
!***********************************************************************!
!
!Initialisations
!
!=================
! 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  )
!
      it = 0
      convergence = ( SQRT(rnorm_k) < prec )
!
      DO  WHILE  ( (.NOT. convergence)  .AND.  (it  <  it_max) )
!
          it  =  it  +  1
!
!===================
! calcul de \alpha_k
!===================
!
!Echange  des  points  aux  interfaces  pour  p
         CALL  p_comm( p, ptabN, ptabE, ptabS, ptabW, ptabAV, ptabAR, voisin )
!
! pmv   :  q_k = A.p_k
         CALL  pmv(  q,  p,  - rnudt2  )
!
! (q_k,p_k)
         CALL prodscal ( alpha_k,  q,  p  )
!
!$OMP MASTER
         itgc = it
         alpha_k = rnorm_k / alpha_k
         beta_k  = rnorm_k
         rnorm_k = zero
!$OMP END MASTER
!
!$OMP BARRIER
!
!================
! 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  )
!
! \beta_k
!$OMP MASTER
         alpha_k = zero
         beta_k  = rnorm_k / beta_k
!$OMP END MASTER
!
!$OMP BARRIER
!
!================
! 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

