!************************************************************************!
!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 3 : decomposition de domaines
!                        gestion des interfaces avec des pointeurs
!                        et des buffers locaux
!
!     Auteur           : Guy Moebs (CRIHAN - France)
!                        <Guy.Moebs@crihan.fr>
!
!***********************************************************************!
!
      PROGRAM Chaleur
!
      USE numerics
!
      IMPLICIT   NONE
!
      INCLUDE 'comthd.h'
!
INTERFACE
!
      include  'i_boucle.h'
!
END INTERFACE
!
!---Declaration  des  variables
!
!Scalaires  gradient  conjugue
      INTEGER   ::  nbitgc
!
!Nombre  d'iterations  maximum  gradient  conjugue
      INTEGER   ::  it_max
!
!Nombre  d'iterations  en  temps
      INTEGER   ::  nbiter
!
!Nombre  de  threads
      INTEGER   ::  nbthd
!
!Temps  ecoule
      REAL(dp)  ::  tps_ecoule, t
!
!
!Frequence  des  sorties
      INTEGER   ::  ifreq
!
!Convergence  gradient  conjugue
      REAL(rp)  ::  prec
!
!***********************************************************************!
!
!Lecture des parametres
      machep = MAX( EPSILON(machep), 1E-15_rp )
!
      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,202) '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 )
!
!************************************************************************!
!
      CALL boucle ( nbthd, nbiter, ifreq, &
                    it_max, prec, nbitgc, tps_ecoule  )
!
!************************************************************************!
!
!Performances  du  code
      t = 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)  & ! nrmerr
          ) / tps_ecoule
      WRITE (6,202) 'Temps ecoule : ', tps_ecoule
      WRITE (6,203) 'Performances cumulees       : ', t, ' MFlops'
      WRITE (6,203) 'Performances par processeur : ', t / REAL(nbthd,rp), ' MFlops'
!
  203 FORMAT ( A, F8.2, A )
!
      STOP
      END  PROGRAM  Chaleur
