      !***********************************************************************!
      SUBROUTINE  initialisation(  u,  f,  u_exact  )
      !***********************************************************************!
!
      USE numerics
!
      IMPLICIT  NONE
!
!--Declaration  des  arguments
!
!Solution  u  a  l'iteration  n  et  n+1
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1,0:1),  INTENT(OUT)  ::  u
!
!Second  membre
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),      INTENT(OUT)  ::  f
!
!Solution exacte
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),      INTENT(OUT)  ::  u_exact
!
!Compteurs
      INTEGER   ::  i,  j,  k
!
!Coordonnees  globales  suivant  x  et  y
      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
                  u(i,j,k,0) = zero
                  u(i,j,k,1) = zero
                  f(i,j,k)   = zero
                  u_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)
!
                  f(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
!$OMP END DO
!
      RETURN
      END  SUBROUTINE  initialisation
