!************************************************************************!
!Resolution de l'equation de Poisson en utilisant Jacobi
!     sur le domaine [0,1]x[0,1]x[0,1] par une methode aux differences finies,
!     avec comme solveur Jacobi.
!     - Delta  u  =  f(x,y,z)
!                 =  - 2 * { z*(z-1) * [ x*(x-1) + y*(y-1) ]
!                          + x*(x-1)*y*(y-1) }
!     u sur les bords vaut 0
!     La solution exacte est u =  x*y*z*(x-1)*(y-1)*(z-1)
!
!     Dans cette version, on se donne
!     le nombre de points interieurs total suivant x (ntx),
!     le nombre de points inferieurs total suivant y (nty) et 
!     le nombre de points inferieurs total suivant z (ntz).
!     ntx =? nty =? ntz
!     hx represente le pas suivant x, hx = 1 / (ntx+1),
!     hy represente le pas suivant y, hy = 1 / (nty+1),
!     hz represente le pas suivant z, hz = 1 / (ntz+1).
!
!!!
!!!
!
!     Cree le          : Tue May 15 18:02:41 2001
!
!     Nature           : parallelisation par OpenMP
!
!     Auteur           : Guy Moebs (CRIHAN - France)
!                        <Guy.Moebs@crihan.fr>
!
!***********************************************************************!
!
      PROGRAM Poisson
!
      USE numerics
!
      IMPLICIT   NONE
!
INTERFACE
!
      INCLUDE 'i_domaine.h'
!
      INCLUDE 'i_initialisation.h'
!
      INCLUDE 'i_calcul.h'
!
      INCLUDE 'i_converge.h'
!
      INCLUDE 'i_erreur.h'
!
END INTERFACE
!
!---Declaration  des  variables
!
!Solution  u  a  l'iteration  n  et  n+1
      REAL(rp),  DIMENSION(:,:,:,:),  ALLOCATABLE  ::  u
!
!Solution exacte
      REAL(rp),  DIMENSION(:,:,:),    ALLOCATABLE  ::  u_exact
!
!Second membre
      REAL(rp),  DIMENSION(:,:,:),    ALLOCATABLE  ::  f
!
!Pointeurs  de  tableau
      INTEGER   ::  iold = 0
      INTEGER   ::  inew = 1
!
!Nombre  iterations  en  temps
      INTEGER   ::  it, itg
!
!Nombre  iterations  maximum  en  temps
      INTEGER   ::  it_max
!
!Nombre de threads
      INTEGER   ::  nbproc
!$    EXTERNAL      OMP_GET_THREAD_NUM
!$    INTEGER   ::  OMP_GET_THREAD_NUM
!
!$    EXTERNAL      OMP_GET_NUM_THREADS
!$    INTEGER   ::  OMP_GET_NUM_THREADS
!
!Code erreur
      INTEGER   ::  ierr = 0
!
!Convergence  et  erreur
      REAL(rp)  ::  gdiff
      REAL(rp)  ::  gerr
!
!Mesure  du  temps
      REAL(dp)  ::  t1, t2
      REAL(dp)  ::  timef
!
!Frequence des sorties
      INTEGER   ::  ifreq
!
!Compteurs
      INTEGER   ::  j
!
!Consommation memoire
      REAL(rp)  ::  gmem
!
!Test  de  convergence
      LOGICAL   ::  convergence
      REAL(rp)  ::  prec
!
!Precision machine
      REAL(rp)  ::  machep
!
!***********************************************************************!
!
!Critere de convergence
      machep = EPSILON(machep)
      prec = MAX( machep, 1E-10_rp )
!
!Connaitre  le  nombre  de  points  interieurs  total
      READ  (5,*) ntx, nty, ntz, it_max, ifreq
!
! Affichage des parametres de l'execution
      WRITE (6,1010) 'Equation de Poisson 3D'
      WRITE (6,1020) 'Critere d arret, prec = ', prec
      WRITE (6,1020) 'Precision machine, machep = ', machep
      WRITE (6,1030) 'ntx    : nombre de points en x = ', ntx
      WRITE (6,1030) 'nty    : nombre de points en y = ', nty
      WRITE (6,1030) 'ntz    : nombre de points en z = ', ntz
      WRITE (6,1030) 'it_max : nombre maximal d iterations pour la convergence = ', it_max
      WRITE (6,1030) 'ifreq  : frequence d affichage des residus = ', ifreq
!
 1010 FORMAT ( A, I5 )
 1020 FORMAT ( A, E12.4, / )
 1030 FORMAT (A, I6, /)
!
!Determiner  les  indices  de  calcul
      CALL domaine(  )
!
! Consommation memoire
      gmem = 0.0234375_rp * REAL(ex-sx+3, rp) * REAL(ey-sy+3, rp) * REAL(ez-sz+3, rp)
      WRITE (6,1040) gmem
 1040 FORMAT ( 'Volume memoire de (u,f) : ', F14.0, ' ko', /)
!
!Allocation  dynamique  des  tableaux
      ALLOCATE(  u      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1,  0:1), &
                 f      (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 u_exact(sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), STAT = ierr )
      IF ( ierr /= 0 ) THEN
         WRITE (6,1010) 'ierr = ', ierr
         WRITE (6,1010) 'Probleme allocation'
         STOP
      END IF
!
!Mesure  du  temps  en  seconde  dans  la  boucle  iterative
      t1  = timef ()
!
!Schema  iteratif  en  temps
      it  =  0
      gdiff = zero
      gerr = zero
      convergence  =  .FALSE.
!
!************************************************************************!
!
!Debut de la region parallele
!
!$OMP PARALLEL DEFAULT( NONE ) &
!$OMP SHARED( u, f, u_exact ) &
!$OMP SHARED( itg, nbproc ) &
!$OMP SHARED( prec, it_max, ifreq, gdiff, gerr ) &
!$OMP PRIVATE( j ) &
!$OMP FIRSTPRIVATE( iold, inew, it, convergence )
!
!$OMP MASTER
      nbproc = 1
!$    nbproc = OMP_GET_NUM_THREADS ()
      WRITE (6,1030) 'Nombre de threads : ', nbproc
!
!Residus / erreurs
      WRITE (6,200)
      WRITE (6,201)
      WRITE (6,200)
  200 FORMAT ( '+-----------+------------+')
  201 FORMAT ( '| Iteration |   Ecart    |')
!$OMP END MASTER
!
!Initialisation  du  second  membre  et  de  la  solution  u  sur  les  bords
      CALL  initialisation(  u,  f,  u_exact  )
!
      DO  WHILE  ( (.NOT. convergence) .AND. (it < it_max) .AND. (gdiff < 50.0_rp) )
!
         it  =  it  +  1
         j = iold
         iold = inew
         inew = j
!
!$OMP BARRIER
!
!$OMP MASTER
         gdiff = zero
!$OMP END MASTER
!
!Calcul  de  u  a  l'iteration  it+1
         CALL  calcul(  u,  f,  inew  )
!
!Calcul  de  l'ecart  sur  le  domaine
         CALL  converge(  u,  gdiff  )
!
!Arret  du  programe  si  on  a  atteint  la  precision  souhaitee
         convergence  =  ( gdiff  <  prec )
!
!Synchronisation pour gdiff
!$OMP BARRIER
!
!Affichage  pour  le  processeur  0  de  la  difference
!$OMP MASTER
         itg = it
         IF  ( MOD(it,ifreq)  ==  0 ) WRITE (6,202) itg, gdiff
  202 FORMAT ( '|  ', I6, '   |', E11.4, ' |')
!$OMP END MASTER
!
      END  DO
!
!************************************************************************!
!
!Validation avec une solution analytique
      CALL  erreur(  u,  u_exact,  inew,  gerr  )
!
!Fin   de la region parallele
!$OMP END PARALLEL
!
!Mesure  du  temps  a  la  sortie  de  la  boucle
      t2  = ( timef () - t1 ) * 0.001_dp
!
!Affichage  du  temps  de  convergence  par  le  processeur  0
      IF ( MOD(itg, ifreq)  /=  0 )  THEN
         WRITE (6,202) itg, gdiff
      END IF
      WRITE (6,200)
!
      t1 = 12.0E-06_rp * REAL(ntx,rp) * REAL(nty,rp) * REAL(ntz,rp) * REAL(itg,rp) / t2
      WRITE (6,203) 'Erreur maximale : ', gerr
      WRITE (6,203) 'Temps ecoule : ', t2
      WRITE (6,204) 'Performances cumulees       : ', t1, ' MFlops'
      WRITE (6,204) 'Performances par processeur : ', t1 / REAL(nbproc,rp), ' MFlops'
!
  203 FORMAT ( A, E12.4 )
  204 FORMAT ( A, F8.2, A )
!
!Desallocation  des  tableaux  u,  u_exact  et  f
      DEALLOCATE(  u,  u_exact,  f)
!
      STOP
      END  PROGRAM  Poisson

