      !***********************************************************************!
      SUBROUTINE  initialisation(  w,  p,  r,  b,  w_exact, rnormk  )
      !***********************************************************************!
!
      USE numerics
!
      IMPLICIT  NONE
!
INTERFACE
!
      INCLUDE 'i_saxpy.h'
!
      INCLUDE 'i_prodscal.h'
!
      INCLUDE 'i_pmv.h'
!
END INTERFACE
!
!--Declaration  des  arguments
!
!Solution  w
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  w
!
!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
!
!Solution  exacte
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  w_exact
!
! norme  du  residu  initial
      REAL(rp),  INTENT(OUT)  ::  rnormk
!
!Compteurs
      INTEGER   ::  i,  j,  k
!
!Coordonnees  globales  suivant  x,  y,  z
      REAL(rp)  ::  x,   y,   z
      REAL(rp)  ::  ax,  ay,  az
!
!
      !***********************************************************************!
!
!Initialisation  des  matrices
!$OMP DO
      DO  k  =  sz - 1,  ez + 1
          DO  j  =  sy - 1,   ey + 1
              DO  i  =  sx - 1,  ex + 1
                  w(i,j,k) = zero
                  p(i,j,k) = zero
                  w_exact(i,j,k) = zero
               END DO
            END DO
         END DO
!$OMP END DO
!
!Initialisation  du  second  membre  et  calcul  de  la  solution  exacte
!$OMP DO
      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)
!
                  w(i,j,k) = 0.2_rp                 
                  b(i,j,k)  =  ax * ay * az &
                            -  2.0_rp * rlambda * ( ( ax + ay ) * az + ax * ay )
                  w_exact(i,j,k)  =  ax * ay * az
!
              END  DO
          END  DO
      END  DO
!$OMP END DO
!
!=================
! pmv   : r = Aw
!=================
      CALL  pmv(  r,  w  )
!
!===================
! saxpy : r = b - Aw
!===================
      CALL  saxpy(  r,  -one,  b,  one  )
!
!=======
! p <- r
!=======
!$OMP DO
      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
!$OMP END DO
!
!==========
! || r ||^2
!==========
      CALL prodscal ( rnormk, r, r )
!
      RETURN
      END  SUBROUTINE  initialisation
