!************************************************************************!
!Resolution de l'equation de la Chaleur sur le domaine [0,1]x[0,1]x[0,1]
!     discretisation spatiale : differences finies 2nd ordre
!     discretisation temporelle : Crank-Nicholson, 2nd ordre
!     avec comme solveur le Gradient Conjugue.
!     d(u)/dt - \lambda * Delta  u =  f(x,y,z)
!     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 Sept 11 11:02:41 2001
!
!     Nature           : parallelisation par OpenMP
!                        version 1 : insertion classique de directives
!
!     Auteur           : Guy Moebs (CRIHAN - France)
!                        <Guy.Moebs@crihan.fr>
!
!***********************************************************************!
!
      PROGRAM Chaleur
!
      USE numerics
!
      IMPLICIT   NONE
!
INTERFACE
!
      INCLUDE     'i_domaine.h'
!
      INCLUDE     'i_fctfx.h'
!
      INCLUDE     'i_scdmb.h'
!
      INCLUDE     'i_gradconj.h'
!
      INCLUDE     'i_nrmerr.h'
!
END INTERFACE
!
!---Declaration  des  variables
!
!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
!
!Nombre  iterations  en  temps
      INTEGER   ::  itg
!
!Nombre  d'iterations  maximum  gradient  conjugue
      INTEGER   ::  it_max
!
!Nombre  d'iterations  en  temps
      INTEGER   ::  nbiter
!
!Nombre  de  threads
      INTEGER   ::  nbthd
!
!$    EXTERNAL      OMP_GET_NUM_THREADS
!$    INTEGER   ::  OMP_GET_NUM_THREADS
!
!Code  erreur
      INTEGER   ::  ierr = 0
      INTEGER   ::  i
      LOGICAL   ::  fichex
!
!Scalaires  gradient  conjugue
      REAL(rp)  ::  rnorm_k
      INTEGER   ::  itgc
      INTEGER   ::  nbitgc
!
!Scalaire  Crank  Nicholson
      REAL(rp)  ::  rnudt2
!
!Normes  diverses
      REAL(rp),  DIMENSION(2)  ::  gnrm2
      REAL(rp),  DIMENSION(2)  ::  gnrmoo
!
!Temps  courant
      REAL(rp)  ::  tps
!
!Mesure  du  temps
      REAL(dp)  ::  t1, t2
      REAL(dp)  ::  timef
!
!Frequence  des  sorties
      INTEGER   ::  ifreq
!
!Consommation  memoire
      REAL(rp)  ::  gmem
!
!Convergence  gradient  conjugue
      REAL(rp)  ::  prec
!
!
!***********************************************************************!
!
!Lecture des parametres
      machep = MAX( EPSILON(machep), 1E-15_rp )
!
      INQUIRE (file='C3D.dat', exist=fichex)
      IF ( .NOT. fichex ) THEN
         WRITE (6,1020) 'STOP STOP STOP STOP STOP'
         WRITE (6,1020) 'Le fichier d entree blas.dat n est pas present'
         WRITE (6,1020) 'Le programme s arrete'
 1020 FORMAT (A)
         STOP
      END IF
!
      OPEN(19,  FILE='C3D.dat',  STATUS='OLD',  FORM='FORMATTED')
      READ(19, 101) nbthd
      READ(19, 101) ntx
      READ(19, 101) nty
      READ(19, 101) ntz
      READ(19, 101) nbiter
      READ(19, 102) dt
      READ(19, 102) rlambda
      READ(19, 102) prec
      READ(19, 101) it_max
      READ(19, 101) ifreq
      CLOSE(19)
!
  101 FORMAT ( I5 )
  102 FORMAT ( E14.4 )
!
      WRITE (6,201) 'Equation de la chaleur 3D'
      WRITE (6,202) 'machep = ', machep
!$    WRITE (6,201) 'Nombre de threads OpenMP : ', nbthd
      WRITE (6,201) 'Nombre de points en X : ', ntx
      WRITE (6,201) 'Nombre de points en Y : ', nty
      WRITE (6,201) 'Nombre de points en Z : ', ntz
!
      WRITE (6,201) 'Nombre d iterations en temps : ', nbiter
      WRITE (6,202) 'Pas de temps : ', dt
      WRITE (6,202) 'Viscosite cinematique : ', rlambda
!
      WRITE (6,202) 'Precision gradient    : ', prec
      WRITE (6,201) 'Nbre max d iterations : ', it_max
      WRITE (6,201) 'Frequence des sorties : ', ifreq
!
  201 FORMAT ( A, I5 )
  202 FORMAT ( A, E12.4 )
!
!Determiner  les  indices  du  domaine
      CALL domaine(  )
!
! Consommation memoire
      gmem = 0.0625_rp * REAL(ex-sx+3, rp) * REAL(ey-sy+3, rp) * REAL(ez-sz+3, rp)
      WRITE (6,1030) gmem
 1030 FORMAT ( 'Volume memoire de (u,b,p,q,r) : ', F14.0, ' ko', /)
!
!Allocation  dynamique  des  tableaux  u, u_exact, f,  b,  p,  q,  r
      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), &
                 u_exact(sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), STAT = ierr )
      IF ( ierr /= 0 ) THEN
         WRITE (6,201) 'ierr = ', ierr
         WRITE (6,201) 'Probleme allocation'
         STOP
      END IF
!
!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
!
!Mesure  du  temps  en  seconde  dans  la  boucle  en  temps
      t1  = timef ()
!
!************************************************************************!
!
!Debut de la region parallele
!
      CALL OMP_SET_NUM_THREADS ( nbthd )
!
!$OMP PARALLEL DEFAULT( NONE ) &
!$OMP SHARED( u, b, u_exact, f1, f2, p, q, r ) &
!$OMP SHARED( tps, dt, rnorm_k, itgc, nbitgc ) &
!$OMP SHARED( gnrm2,  gnrmoo ) &
!$OMP SHARED( nbthd, nbiter, prec, machep ) &
!$OMP SHARED( ifreq, it_max, rnudt2, rlambda ) &
!$OMP PRIVATE( itg )
!
!$OMP MASTER
         nbthd = 1
!$       nbthd = OMP_GET_NUM_THREADS ()
!$       WRITE (6,1040) 'Nombre de threads dans la region parallele : ', nbthd
!$OMP END MASTER
 1040 FORMAT (A, I5, /)
!
!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
!=================================
!
      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
!
!Calcul  second  membre  Crank  Nicholson
         CALL  scdmb(  b,  f1,  f2,  u,  tps,  rnudt2  )
!
!Avancement en temps
!$OMP MASTER
         tps = tps + dt
!$OMP END MASTER
!
!Resolution  gradient  conjugue
         CALL  gradconj(  u,  b,  p,  q,  r,  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
!
!************************************************************************!
!
!Fin   de la region parallele
!$OMP END PARALLEL
!
!Mesure  du  temps  a  la  sortie  de  la  boucle
      t2  =  ( timef () - t1 ) * 0.001_dp
!
      nbitgc = nbitgc + itgc
!
!Performances  du  code
      t1 = REAL(ntx,rp) * REAL(nty,rp) * REAL(ntz,rp) * ( &
           23.0E-06_rp  * REAL(nbitgc,rp) &                  ! gradconj
         + 30.0E-06_rp  * REAL(nbiter)    &                  ! scdmb
         + 12.0E-06_rp  * REAL(nbiter/ifreq + 1, rp) ) / t2  ! nrmerr
      WRITE (6,202) 'Temps ecoule : ', t2
      WRITE (6,203) 'Performances cumulees       : ', t1, ' MFlops'
      WRITE (6,203) 'Performances par processeur : ', t1 / REAL(nbthd,rp), ' MFlops'
!
  203 FORMAT ( A, F8.2, A )
!
!Desallocation  des  tableaux  u,  u_exact,  b,  p,  q,  r  et  f
      DEALLOCATE(  u,  u_exact,  b,  f1,  f2,  q,  p,  r  )
!
      STOP
      END  PROGRAM  Chaleur
