      !***********************************************************************!
      SUBROUTINE  initialisation(  u,  b,  p,  r,  u_exact, rnorm_k, &
                                   ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d  )
      !***********************************************************************!
!
      USE MPI
!
      USE numerics
!
      IMPLICIT  NONE
!
INTERFACE
!
      INCLUDE 'i_communication.h'
!
      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(OUT)  ::  u
!
!Solution exacte
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  u_exact
!
!Second  membre
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  b
!
!Vecteurs intermediaires
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  p
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  r
!
!Norme initiale
      REAL(rp),  INTENT(OUT)  ::  rnorm_k
!
!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
!
!Compteurs
      INTEGER   ::  i,  j,  k
!
!Coordonnees  globales  suivant  x  et  y
      REAL(rp)  ::  x,  y,  z
      REAL(rp)  ::  ax, ay, az
!
!
      !***********************************************************************!
!
!Initialisation  des  matrices
      DO  k  =  sz - 1,  ez + 1
          DO  j  =  sy - 1,  ey + 1
              DO  i  =  sx - 1,  ex + 1
                  u(i,j,k) = zero
                  p(i,j,k) = zero
              END DO
          END DO
      END DO
!
!Pas d'espace selon x, y et z
      hx  =  1.0_rp  /  real( ntx + 1, rp)
      hy  =  1.0_rp  /  real( nty + 1, rp)
      hz  =  1.0_rp  /  real( ntz + 1, rp)
!
!Initialisation  du  second  membre
      DO  k  =  sz,  ez
          z  =  REAL(k, rp) * hz
          az =  z * (z - 1.0_rp)
!
          DO  j  =  sy,  ey
              y  =  REAL(j, rp) * hy
              ay =  y * (y - 1.0_rp)
!
              DO  i  =  sx,  ex
                  x  =  REAL(i, rp) * hx
                  ax =  x * (x - 1.0_rp)
!
                  u(i,j,k) = 0.2_rp
                  b(i,j,k) = 2.0_rp * ( ( ax + ay ) * az + ax * ay )
                  u_exact(i,j,k)   =  ax * ay * az
!
              END  DO
          END  DO
      END  DO
!
!Echange  des  points  aux  interfaces  pour  u
      CALL communication(  u,  ip_xy,  ip_xz,  ip_yz,  voisin,  comm3d  )
!
!=================
! pmv   : r = Au
!=================
      CALL  pmv(  r,  u  )
!
!===================
! 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  )
!
      RETURN
      END  SUBROUTINE  initialisation
