      !***********************************************************************!
      SUBROUTINE  boucle(  nbthd, nbiter, ifreq,  &
                           it_max,  prec, nbitgc, t2  )
      !***********************************************************************!
!
      USE numerics
!
      IMPLICIT  NONE
!
      INCLUDE 'comthd.h'
!
      INCLUDE 'p_comptr.h'
!
      INCLUDE 'u_comptr.h'
!
!--Interface
      INTERFACE
!
      include 'i_gm_dims_create.h'
!
      include 'i_gm_coords_create.h'
!
      include 'i_domaine.h'
!
      include 'i_voisinage.h'
!
      include 'i_assign_ptr.h'
!
      include 'i_fctfx.h'
!
      include 'i_nrmerr.h'
!
      include 'i_u_comm.h'
!
      include 'i_scdmb.h'
!
      include 'i_gradconj.h'
!
      END INTERFACE
!
!--Declaration  des  arguments
!
!Nombre  de  threads
      INTEGER,  INTENT(INOUT)  ::  nbthd
!
!Nombre  d'iterations  en  temps
      INTEGER,  INTENT(IN)  ::  nbiter
!
!Nombre  d'iterations  maximum  gradient  conjugue
      INTEGER,  INTENT(IN)  ::  it_max
!
!Frequence  des  sorties
      INTEGER,  INTENT(IN)  ::  ifreq
!
!Precision convergence GRadient Conjugue
      REAL(rp),  INTENT(IN)  ::  prec 
!
!Nombre total d'iteration Grasdient Conjugue
      INTEGER,   INTENT(OUT)  ::  nbitgc
!
!Temps ecoule
      REAL(dp),  INTENT(OUT)  ::  t2
!
!--Declaration des variales locales
!
!Solution  u
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  u
!
!Solution  exacte
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  u_exact
!
!Second  membre
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  b
!
!Forcage  exterieur
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  f1
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  f2
!
!Vecteurs  de  travail
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  p
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  q
      REAL(rp),  DIMENSION(:,:,:),  ALLOCATABLE  ::  r
!
!Tableau contenant toutes les adresses memoire
      INTEGER(ip),  DIMENSION(:,:),  ALLOCATABLE ::  pglobal_tab
      INTEGER(ip),  DIMENSION(:,:),  ALLOCATABLE ::  uglobal_tab
!
!Buffers locaux
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  ptabN
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  ptabE
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  ptabS
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  ptabW
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  ptabAV
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  ptabAR
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  utabN
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  utabE
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  utabS
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  utabW
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  utabAV
      REAL(rp),  DIMENSION(:,:),  ALLOCATABLE  ::  utabAR
!
!Normes  diverses
      REAL(rp),  DIMENSION(2)  ::  gnrm2
      REAL(rp),  DIMENSION(2)  ::  gnrmoo
!
!Temps  courant
      REAL(rp)  ::  tps
!
!Mesure  du  temps
      REAL(dp)  ::  t1
      REAL(dp)  ::  timef
!
!Scalaires  gradient  conjugue
      REAL(rp)  ::  rnorm_k
      INTEGER   ::  itgc
      INTEGER   ::  itg
!
!Scalaire  Crank  Nicholson
      REAL(rp)  ::  rnudt2
!
      REAL(rp)  ::  gmem
!
!$    EXTERNAL      OMP_GET_THREAD_NUM
!$    INTEGER   ::  OMP_GET_THREAD_NUM
!
!$    EXTERNAL      OMP_GET_NUM_THREADS
!$    INTEGER   ::  OMP_GET_NUM_THREADS
!
!Tableau  contenant  les  voisins  du  domaine  courant
      INTEGER,  PARAMETER                  ::  NB_VOISINS  =  6
!
      INTEGER,  DIMENSION(NB_VOISINS)      ::  voisin
      INTEGER,  PARAMETER                  ::  N=1,  E=2,   S=3
      INTEGER,  PARAMETER                  ::  W=4,  AV=5,  AR=6
!
!Constantes  decomposition de domaine
      INTEGER,  PARAMETER  ::  nbdims = 3
      INTEGER,  DIMENSION(1:nbdims)  ::  dims
      LOGICAL,  DIMENSION(1:nbdims)  ::  periods
      INTEGER,  DIMENSION(1:nbdims)  ::  coords
!
      INTEGER  ::  ierr, i
!
!***********************************************************************!
!
!Initialisation des scalaires communs
      rnudt2 = 0.5_rp * rlambda * dt
!
      itgc   = 0
      nbitgc = 0
!
      rnorm_k = zero
      tps     = zero
!
      DO  i  =  1, 2
         gnrm2(i)  = zero
         gnrmoo(i) = zero
      END DO
!
!Debut de la region parallele
!
!$    CALL OMP_SET_NUM_THREADS ( nbthd )
!
      ALLOCATE ( uglobal_tab (6, 0:nbthd-1),  pglobal_tab (6, 0:nbthd-1), STAT = ierr )
      IF ( ierr /= 0 ) THEN
         WRITE (6,201) 'ierr = ', ierr
         WRITE (6,201) 'Probleme allocation'
  201 FORMAT ( A, I5 )
         STOP
      END IF
!
!$OMP PARALLEL DEFAULT( NONE ) &
!$OMP SHARED( uglobal_tab, pglobal_tab) &
!$OMP SHARED( tps, dt, rnorm_k, itgc, nbitgc ) &
!$OMP SHARED( nbthd, nbiter, prec, machep ) &
!$OMP SHARED( ifreq, it_max, rnudt2, rlambda ) &
!$OMP SHARED( ntx, nty, ntz, t1, t2, gnrm2, gnrmoo ) &
!$OMP PRIVATE( itg, ierr, voisin, dims, periods, coords ) &
!$OMP PRIVATE( p, q, r, b, u, u_exact, f1, f2, gmem, i ) &
!$OMP PRIVATE( utabN, utabE, utabS, utabW, utabAV, utabAR ) &
!$OMP PRIVATE( ptabN, ptabE, ptabS, ptabW, ptabAV, ptabAR )
!
         ierr = 0
!
         mythd = 0
!$       mythd = OMP_GET_THREAD_NUM ()
!
!$OMP MASTER
         nbthd = 1
!$       nbthd = OMP_GET_NUM_THREADS ()
         WRITE (6,'(A,I5,/)') 'Nombre de threads OpenMP : ', nbthd
!$OMP END MASTER
!
!Connaitre
!  le  nombre  de  processeurs  selon  x
!  le  nombre  de  processeurs  selon  y
!  le  nombre  de  processeurs  selon  z  en  fonction  du  nombre  total  de  processeurs
      dims(1:nbdims)  =  0
      IF ( ntx == 1 ) dims(1) = 1
      IF ( nty == 1 ) dims(2) = 1
      IF ( ntz == 1 ) dims(3) = 1
!
      CALL gm_dims_create( nbdims, dims, nbthd, ierr )
!
!Creation  de  la  grille  de  processeurs  3D  sans  periodicite
      periods(:)  =  .FALSE.
!
      CALL gm_coords_create( nbdims, dims, coords )
!
!Recherche  de  ses  6  voisins  pour  chaque  processeur
      CALL  voisinage( nbdims, dims, coords, periods, voisin )
!
!Determiner  les  indices  du  domaine
      CALL domaine( nbdims, dims, coords )
!
!Allocation  dynamique  des  tableaux  b, u_exact, f1, f2, q, et r
      ierr = 0
      ALLOCATE(  u  (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 b  (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 f1 (sx  :ex  , sy  :ey  , sz  :ez  ), &
                 f2 (sx  :ex  , sy  :ey  , sz  :ez  ), &
                 p  (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 q  (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 r  (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), &
                 utabE  (sy:ey,sz:ez), &
                 utabW  (sy:ey,sz:ez), &
                 utabS  (sx:ex,sz:ez), &
                 utabN  (sx:ex,sz:ez), &
                 utabAV (sx:ex,sy:ey), &
                 utabAR (sx:ex,sy:ey), &
                 ptabE  (sy:ey,sz:ez), &
                 ptabW  (sy:ey,sz:ez), &
                 ptabS  (sx:ex,sz:ez), &
                 ptabN  (sx:ex,sz:ez), &
                 ptabAV (sx:ex,sy:ey), &
                 ptabAR (sx:ex,sy:ey), &
                 u_exact(sx  :ex  ,sy  :ey  ,sz  :ez  ), STAT = ierr )
      IF ( ierr /= 0 ) THEN
         WRITE (6,201) 'ierr = ', ierr
         WRITE (6,201) 'Probleme allocation'
         STOP
      END IF
!
!Assignation des adresses memoires des buffers des voisins
      CALL assign_ptr ( pglobal_tab,  uglobal_tab,  voisin, nbthd,  &
                        u, utabN, utabE, utabS, utabW, utabAV, utabAR,  &
                        p, ptabN, ptabE, ptabS, ptabW, ptabAV, ptabAR )
!
! Consommation memoire
!$OMP MASTER
      gmem = 0.0625_rp * REAL(ex-sx+3, rp) * REAL(ey-sy+3, rp) * REAL(ez-sz+3, rp)
      WRITE (6,'(A, F14.0, A, /)') 'Volume memoire de (u,b,p,q,r) : ', gmem, ' ko'
!$OMP END MASTER
!
!Initialisation  des matrices
      CALL  fctfx(  u, f1, f2, p, u_exact, tps  )
!
!Erreurs initiales
      CALL  nrmerr(  gnrm2,  gnrmoo,  u,  u_exact,  tps  )
!
!$OMP MASTER
      WRITE (6, 301) 'TPS', 'uL2', 'uLoo', 'eL2', 'eLoo', 'itGC', 'res'
      WRITE (6, 302) tps,  gnrm2(1),  gnrmoo(1),  gnrm2(2),  gnrmoo(2)
  301 FORMAT ( 4x, A3, 7x, A3, 8x, A4, 9x, A3, 8x, A4, 7x, A4, 5x, A3 )
  302 FORMAT ( F8.4, 4E12.4, 4x, I4, E12.4 )
!$OMP END MASTER
!
!=================================
! Debut  de  la  boucle  en  temps
!=================================
!
!Mesure  du  temps  en  seconde  dans  la  boucle  en  temps
!$OMP MASTER
      t1  = timef ()
!$OMP END MASTER
!
      DO itg = 1, nbiter
!
!$OMP MASTER
         DO  i  =  1, 2
            gnrm2(i)  = zero
            gnrmoo(i) = zero
         END DO
!
         nbitgc = nbitgc + itgc
         itgc = 0
         rnorm_k = zero
!$OMP END MASTER
!
!Echange  des  points  aux  interfaces  pour  u
         CALL u_comm( u, utabN, utabE, utabS, utabW, utabAV, utabAR, voisin )
!
!Calcul  second  membre  Crank  Nicholson
         CALL  scdmb(  u,  b,  f1,  f2,  tps,  rnudt2  )
!
!Avancement en temps
!$OMP MASTER
         tps = tps + dt
!$OMP END MASTER
!
!Resolution  gradient  conjugue
         CALL  gradconj( u,  b,  p,  q,  r,  voisin,  &
                         ptabN,  ptabE,  ptabS,  ptabW,  ptabAV,  ptabAR,  &
                         it_max,  prec,  itgc,  rnorm_k,  rnudt2 )
!
!Calcul des normes
         IF  ( MOD(itg,ifreq)  ==  0 ) THEN
!
            CALL  nrmerr(  gnrm2,  gnrmoo,  u,  u_exact,  tps  )
!
!$OMP MASTER
            WRITE (6, 302) tps,  gnrm2(1),  gnrmoo(1),  gnrm2(2),  gnrmoo(2),  itgc,  SQRT(rnorm_k)
!$OMP END MASTER
!
         END IF
!
!=================================
! Fin    de  la  boucle  en  temps
!=================================
!
      END  DO
!
!Mesure  du  temps  a  la  sortie  de  la  boucle
!$OMP MASTER
      t2  =  ( timef () - t1 ) * 0.001_dp
!$OMP END MASTER
!
!Desallocation  des  tableaux  u_exact, b, f1, f2, q et r
      DEALLOCATE(  u_exact, u, b, f1, f2, p, q, r, &
                   utabN, utabE, utabS, utabW, utabAV, utabAR,  &
                   ptabN, ptabE, ptabS, ptabW, ptabAV, ptabAR )
!
!************************************************************************!
!
!Fin   de la region parallele
!$OMP END PARALLEL
!
      DEALLOCATE ( uglobal_tab, pglobal_tab )
!
      nbitgc = nbitgc + itgc
!
      RETURN
      END  SUBROUTINE  boucle
