      !***********************************************************************!
      SUBROUTINE  fctfx(  u,  f1,  f2,  p, u_exact,  tps  )
      !***********************************************************************!
!
      USE numerics
!
      IMPLICIT  NONE
!
!--Declaration  des  arguments
!
!Solution  u
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  u
!
!Vecteur  intermediaire
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  p
!
!Composante  spatiale,  force  exterieure
      REAL(rp),  DIMENSION(sx  :ex  ,sy  :ey  ,sz  :ez  ),  INTENT(OUT)  ::  f1
      REAL(rp),  DIMENSION(sx  :ex  ,sy  :ey  ,sz  :ez  ),  INTENT(OUT)  ::  f2
!
!Solution  exacte
      REAL(rp),  DIMENSION(sx-1:ex+1,sy-1:ey+1,sz-1:ez+1),  INTENT(OUT)  ::  u_exact
!
!Instant courant
      REAL(rp),  INTENT(IN)  ::  tps
!
!--Declaration des variales locales
!
!Compteurs
      INTEGER   ::  i,  j,  k
!
!Coordonnees  globales  suivant  x,  y,  z
      REAL(rp)  ::  x,   y,   z
      REAL(rp)  ::  ax,  ay,  az
!
!Valeur initiale
      REAL(rp)  ::  ut
!
!
      !***********************************************************************!
!
!Instant initial
      CALL fctut( tps, ut )
!
!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) = zero
                  p(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)
!
                  u (i,j,k) =  ax * ay * az * ut
                  f1(i,j,k) =  ax * ay * az
                  f2(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  fctfx
