View difference between Paste ID: 9mXhcjka and K8GMEfMx

 PROGRAM diffusion_1D
 ! *******************************************************************
 !     Auteurs:
 !     Date de creation:
 !     Date de modification:
 !
 !     Ce programme permet de resoudre une equation de diffusion 1D:
 !
 !     dC/dt = K d2C/dx2
 !
 !     C(x,t): la concentration en kg/m3
 !     K: la diffusivite en m2/s
 !     t: le temps en s
 !     x: la variable d'espace en m
 !
 ! *******************************************************************
 !
 ! ==========================
 ! Declaration des variables
 ! ==========================
 !
   IMPLICIT NONE                      ! Toutes les variables doivent être declarees
 !
   INTEGER :: i                       ! indice de la variable spatiale
   INTEGER :: n                       ! indice de la variable temporelle
   INTEGER :: NX                      ! nombre d'intervales suivant Ox
   INTEGER :: NT,NTU                  ! nombre d'iterations totales et demandées
 !
   PARAMETER(NX = 100)                ! Fixe le nombre de points suivant Ox
   PARAMETER(NT = 2000)               ! Fixe le nombre d'iterations max
 !
   REAl :: M0                         ! la masse initiale de polluant
   REAl :: c0                         ! la concentration initiale c(x,0)=c0
   REAL :: K                          ! la diffusivite
   REAL :: U0                          ! la speed of wind
   REAL :: tfin                       ! le temps max a atteindre
   REAL :: dt                         ! le pas de temps
   REAL :: dx                         ! le pas d'espace
   REAL :: L                          ! longueur du domaine de calcul
   REAL :: Mt                         ! Masse totale de colorant
   REAL :: per
 
 !
   REAL, DIMENSION(0:NX, 0:NT)      :: ca   !  Solution analytique
   REAL, DIMENSION(0:NX)      :: x    ! la variable d'espace x
   REAL, DIMENSION(0:NT)      :: t    ! le temps
   REAL, DIMENSION(0:NX,0:NT) :: c    ! la concentration c(x,t)
 !
   INTEGER           :: LENGTH        ! Longeur de la chaine de caracteres du nom des fichiers
   CHARACTER(LEN=4)  :: INT_TO_STRING ! Pour conversion entier -> chaine de caractères
   CHARACTER(LEN=35) :: NOMFICHANIM   ! Nom du fichier de resultats
   INTEGER :: LENGTH2
   CHARACTER(LEN=35) :: NOMFICHANIM2
 !
 ! ===================================================
 ! Affectation des grandeurs physiques et géométriques
 ! ===================================================
 !
   K = .001
   U0 = 0.5
   M0 = 1.
   L = 10.
   dx = L / NX
   c0 = M0 / dx
   per = ACOS(-1.) ! per, c'est pi!
 
 !
 ! =================================
 ! Creation et affichage du maillage
 ! =================================
 !
   WRITE(*,*) '--------------------------------'
   WRITE(*,*) 'Affichage des noeuds du maillage'
   WRITE(*,*) '--------------------------------'
   DO i=0,NX
      x(i) = i*dx
      WRITE(*,*) 'X(',i,')=', x(i)
   ENDDO
 !
 ! ===================================
 ! Definition des paramètres temporels
 ! ===================================
 !
 !
 ! Definition du pas de temps
 !
   WRITE(*,*) 'Donnez votre pas de temps, aussi connu sous le nom de dt :'
   READ(*,*) dt
 !
 ! Calcul du nombre d'iterations a effectuer NTU
 !
   WRITE(*,*) 'Quand est-ce qu''on s''arrete ?'
   READ(*,*) tfin
   NTU = tfin / dt
   IF (NTU.GT.NT) THEN
      WRITE(*,*) 'Désolé garcon, le temps max est: ', NT*dt
      NTU = NT
      WRITE(*,*) 'Le nombre d''iterations sera donc : ', NTU, ', désolé mec'
   ENDIF
 !
 ! Definition du tableau temps
 !
   DO n=0,NTU
      t(n) = n*dt
   ENDDO
 !
 ! =====================================
 ! Mise à zéro du champ de concentration
 ! =====================================
 !
   DO i=0,NX
      if (i<25) then
         c(i, 0) = MAX(0, i-15)
      else
         c(i, 0) = MAX(0, 35-i)
     end if
   ENDDO
 !
 ! ==================
 ! Condition initiale
 ! ==================
 !
   !c(NX/2, 0) = c0
 !
 ! ================================
 ! Calcul du champ de concentration
 ! ================================
 !
   WRITE(*,*) 'Calcul en cours, va chercher un café...'
 !
 ! -----------------------------
 ! Debut des iterations en temps
 ! -----------------------------
   DO n=0,NTU
 !
      WRITE(*,*) '-------------------------------'
      WRITE(*,*) 'Iteration n=',n,' avec t=',t(n)
      WRITE(*,*) '-------------------------------'
 !
 ! Calcul de la solution au pas de temps n+1
 !
     DO i=1,NX-1
         c(i, n+1) = U0*dt/(dx*2) * ( c(i-1, n) - c(i+1, n)) + c(i, n)
         ca(i, n+1) = c0 * (i*dx - U0*n*dt)
     ENDDO
 
 !
 ! Conditions aux limites
 !
     c(0, n+1) = 0
 
 ! ---------------------------
 ! Fin des iterations en temps
 ! ---------------------------
 !
   ENDDO
 
   ! Vérification de la masse totale
     Mt = 0
     DO i=1,NX-1
         Mt = Mt + c(i, NTU)
     ENDDO
     Mt = Mt + .5 * (c(0, NTU) + c(NX, NTU))
     WRITE(*,*) 'M0 = ', M0, ', Mt = ', Mt*dx
 !
   WRITE(*,*) 'Victoire ! Le calcul est termine avec succes \o/'
 !
 !
 ! ==========================================================================
 ! Ecriture du champ de concentration dans des fichiers pour chaque iteration
 ! ==========================================================================
 !
   WRITE(*,*) 'Ecriture des resultats dans des fichiers ...'
 !
   DO n=0,NTU
 !
 ! Definition du nom du fichier
 !
 !
      IF (n.LT.10)                    WRITE (int_to_string, '(I1)') n
      IF ((n.LT.100).AND.(n.GE.10))   WRITE (int_to_string, '(I2)') n
      IF ((n.LT.1000).AND.(n.GE.100)) WRITE (int_to_string, '(I3)') n
      IF ((n.LT.10000).AND.(n.GE.1000)) WRITE (int_to_string, '(I4)') n
 !
      NOMFICHANIM='diffusion_1D_'//int_to_string
      LENGTH=INDEX(NOMFICHANIM,' ')-1
      IF (LENGTH.LT.0) LENGTH=LEN(NOMFICHANIM)
      NOMFICHANIM=NOMFICHANIM(1:LENGTH)//'.vtk'
 
 
      NOMFICHANIM2='diffusion_analytique_'//int_to_string
      LENGTH2=INDEX(NOMFICHANIM2,' ')-1
      IF (LENGTH2.LT.0) LENGTH2=LEN(NOMFICHANIM2)
      NOMFICHANIM2=NOMFICHANIM2(1:LENGTH2)//'.vtk'
 
 !
 ! Ouverture du fichier
 !
      OPEN(                      &
           UNIT   = 9,           &
           FILE   = NOMFICHANIM, &
           STATUS = 'UNKNOWN',   &
           FORM   = 'FORMATTED'  &
           )
 
     OPEN(                      &
           UNIT   = 10,           &
           FILE   = NOMFICHANIM2, &
           STATUS = 'UNKNOWN',   &
           FORM   = 'FORMATTED'  &
           )
 !
 ! Ecriture dans le fichier
 !
      WRITE(9,'(''# vtk DataFile Version 2.0'')')
      WRITE(9,'(''2D scalar data'')')
      WRITE(9,'(''           '')')
 !
      WRITE(9,'(''ASCII'')')
 !
      WRITE(9,'(''DATASET STRUCTURED_GRID'')')
      WRITE(9,'(''DIMENSIONS'',3i7)') NX+1,1,1
      WRITE(9,'(''POINTS '',i7,'' float'')') (NX +1)*1*1
      DO i=0,NX
         WRITE(9,*) x(i),0,0
      ENDDO
 !
      WRITE(9,'(''POINT_DATA '',i7)') (NX+1)*1*1
 !
      WRITE(9,'(''SCALARS concentration float'')')
      WRITE(9,'(''LOOKUP_TABLE default'')')
      DO i=0,NX
         WRITE(9,*) c(i,n)
      ENDDO
 
 
      ! Solution analytique
      WRITE(10,'(''# vtk DataFile Version 2.0'')')
      WRITE(10,'(''2D scalar data'')')
      WRITE(10,'(''           '')')
 !
      WRITE(10,'(''ASCII'')')
 !
      WRITE(10,'(''DATASET STRUCTURED_GRID'')')
      WRITE(10,'(''DIMENSIONS'',3i7)') NX+1,1,1
      WRITE(10,'(''POINTS '',i7,'' float'')') (NX +1)*1*1
      DO i=0,NX
         WRITE(10,*) x(i),0,0
      ENDDO
 !
      WRITE(10,'(''POINT_DATA '',i7)') (NX+1)*1*1
 !
      WRITE(10,'(''SCALARS concentration float'')')
      WRITE(10,'(''LOOKUP_TABLE default'')')
      DO i=0,NX
         WRITE(10,*) ca(i, n)
      ENDDO
 !
   ENDDO
 !
 !     Fermeture du fichier
 !
   CLOSE(9)
   CLOSE(10)
 !
 !     FIN DU PROGRAMME
 !----------------------------------------------------------------------
 END PROGRAM diffusion_1D