      program Semidif
c     resolution numerique de l equation  SEMI DIFferentielle
c            D 1/2 u = thf*f(u(t)) + thg*g(t)
c                                   francois dubois et stephanie mengue
c=======================================================================
c
c     Version 1 : cf rapport interne IAT n° 334/2000
c       7 1 2K    23 01 2K     06-02-2K     22-04-2000   06 mai 2000
c                                        Tourne sur HP   26 mai 2000
c=======================================================================
c     Version 2  : introduction du cinquieme scheme
c                  methode mixte par collocation
c                  26 juin 2000     18 aout 2000
c=======================================================================
c     Version 3  : Condition initiale non nulle
c                  20 octobre 2000
c=======================================================================
c     Version 4  : Suite de la methode mixte
c                  u appartient a P1 et f a P1 aussi
c                  08  novembre  2000, 30 nov 2k, 21 dec 2000.
c                  version asci   
c                  Marne La Vallee le 12 janvier 2001.
c                  21 mai 2001, 25 juin 2001,    version aci.
c
c                  voir un rapport de recherche a venir
c=======================================================================
c     fd et sm,   12 janvier 2001, 21, 22  mai 2001.
c=======================================================================
c
c     parametres de commande
c
c        ischema   choix du schema numerique
c               1  racine carree du schema a deux points
c               2  racine  carree du schema a trois points
c               3  schema de msallam
c               4  elements finis P1
c               5  methode mixte u dans P1, f(u) dans P0.
c               6  methode mixte u dans P1, f(u) dans P1.
c        joulin 0  Ces general a voir plus bas
c               1  cas test de Joulin
c        idrap     choix de la fonction f(u)
c               0  f = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine  carree de u
c               3                        u puissance deux tiers
c               4                                    trois quarts
c               5                                    quatre cinquiemes
c               6      proportionnelle a u
c
c        ifonc     choix de la fonction g(t)
c               0  g = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine carree de t
c               3      fonction lineaire
c               4      proportionnelle a t multiplie par racine de t
c               5      quadratique
c               6      harmonique               ref = Lubich 86
c               7      sinus de racine de t     ref = Lubich 86
c               8      polynome de racine de t  ref = Diethelm
c
c        thf       parametre multiplicatif de la fonction f
c        thg       parametre multiplicatif de la fonction g
c        uinit     valeur initiale
c        np        nombre de pas de temps
c        tfin      temps final
c        dt        valeur du pas de temps
c        cc        coefficients pour le calcul du second membre
c        ccdd      idem pour le tableau des gradients des flux
c        uu        vecteur des solutions discretes
c        phi       tableau des flux moyens dans un intervalle de temps
c        dphi      idem pour leur derivee en temps
c

      parameter (nmax=2600000)
      implicit double precision (a-h,o-z)
      dimension cc(nmax) , ccdd(nmax) ,
     &          uu(nmax) , phi(nmax) , dphi(nmax)
      common / joul1 / joulin, beta, eee
c              beta     parametre de la fonction du test de Joulin
c              eee      idem
c=================================================================
c     initialisation des parametres
c
      idata=30
      open (unit=idata, file = 'input4' ,form='formatted')
      call init ( idata , ischema ,  idrap , ifonc , thf , thg ,
     &                  uinit ,  np , nmail , tfin)
c
      npmax = nmax
      if ((joulin.ne.0).and.(ischema.ne.5)) then
      print * , ' Pour le cas test de Joulin, il faut utiliser',
     &          ' le schema numero 5 '
      print * , ' Tchao ! '
      stop
      end if
c                                              Boucle sur les maillages
c
      open (unit=11, file = 'log-erreurs' ,form='formatted')
      open (unit=12, file = 'erreurs' ,form='formatted')
      do 1 ,  j = 1 , nmail
      dt        = tfin / np
c
      print * , 'numero du maillage ', j ,
     &          '          nb de points en temps ' , np
c
c     pretraitement :  coefficients de l'integrale de convolution
      call coeff (ischema , npmax , dt , czero , cc , ccdd )
c
c     boucle en temps
      call  temps (ischema , idrap , ifonc , thf , thg , uinit ,
     &             np , dt ,czero , cc , ccdd , uu , phi , dphi )
c
c                              Calcul de l erreur et sauvegarde l'icelle
      errlinf = 0.d0
      errl2   = 0.d0
      if (ifonc.eq.0)  goto 2
      if (idrap.eq.0)  goto 2
      if (idrap.ne.ifonc) then
c      print * , ' Je ne garantis rien pour le calul des erreurs '
c      print * , ' Etes vous bien sur de la coherence de vos donneess ?'
      end if
 2    do 10 k = 1 , np
      if (ifonc.eq.0) then
      if (idrap.eq.1)    uex = sqrt(k*dt)
      if (idrap.eq.2)    uex = k*dt
      if (idrap.eq.3)    uex = k*dt * sqrt(k*dt)
      if (idrap.eq.4)    uex = (k*dt) **2
      if (idrap.eq.5)    uex = k*dt * k*dt* sqrt(k*dt)
      end if
      if (idrap.eq.0) then
      if (ifonc.eq.1)    uex = sqrt(k*dt)
      if (ifonc.eq.2)    uex = (k*dt)
      if (ifonc.eq.3)    uex = k*dt * sqrt(k*dt)
      if (ifonc.eq.4)    uex = (k*dt) **2
      if (ifonc.eq.5)    uex = k*dt * k*dt* sqrt(k*dt)
      end if
      if (idrap.eq.ifonc) then
      if (idrap.eq.1)    uex = sqrt(k*dt)
      if (idrap.eq.2)    uex = k*dt
      if (idrap.eq.3)    uex = k*dt * sqrt(k*dt)
      if (idrap.eq.4)    uex = (k*dt) **2
      if (idrap.eq.5)    uex = k*dt * k*dt* sqrt(k*dt)
      end if
      errlinf = max ( errlinf , abs(uu(k)-uex) )
c
      if (k.lt.np) then
      errl2   = errl2 + ( abs(uu(k)-uex) )**2
      else
      errl2   = errl2 + 0.5d0 * ( abs(uu(k)-uex) )**2
      end if
  10  continue
      errl2   = sqrt(dt) * sqrt(errl2)
c
      err1    = log(errlinf)/log(2.d0)
      err2    = log(errl2)/log(2.d0)
      write(12,*) np, errlinf, errl2
      write(11,*) log(np+0.d0) / log(2.d0) , err1, err2
c
c                        sauvegarde des derniers resultats sur fichier
      if (j.eq.nmail) then
      open (unit=10, file = 'resultat' ,form='formatted')
      write(10,998) 0.d0, uinit
      if (np.le.1000) then
      do 100  jj = 1 , np
      write(10,998) jj*dt , uu(jj)
 100  continue
      else
      jkk       = np / 1000
      print * , ' pas d echantillonnage pour ecriture ' , jkk
c
c                                              modif fd     21 mai 01
      write(10,998) dt , uu(1)
      do 110  l = 1 , 1000
      pro       = dfloat(l) / 1000.d0
      tapp      = pro * tfin
      jj        = 1 + tapp / dt
c      print * , jj
      write(10,998) jj*dt , uu(jj)
c                                          fin modif fd     21 mai 01
 110  continue
      end if
c
      end if
c
c      print 997 , ' np = ' , np , '              uu(np) = ' , uu(np)
      np        = 2 * np
 1    continue
c
      stop
 998  format (e25.18, 1X, e25.18)
 997  format (a , i8 , a ,  e25.18)
      end
c
c
c==========================================================================
c
      subroutine init ( idata , ischema ,  idrap , ifonc , thf , thg ,
     &                  uinit ,  np , nmail , tfin)
c
c     lecture des donnees sur fichier
c
c        ischema   choix du schema numerique
c        joulin    vous voudriez bien le savoir, hein ?
c        idrap     choix de la fonction f(u)
c        ifonc     choix de la fonction g(t)
c        thf       parametre multiplicatif de la fonction f
c        thg       parametre multiplicatif de la fonction g
c        uinit     valeur initiale
c        np        nombre de pas de temps
c        nmail     nombre de maillages
c        tfin      temps final
c        beta      parametre de la fonction du test de Joulin
c        eee       idem
c=================================================================
c     fd et sm     21 decembre 2k.
c=================================================================
c
      implicit double precision (a-h,o-z)
      common / joul1 / joulin, beta, eee
      character*72 phrase
c
      print 999 , '              SEMIDIF4        25 juin  2001'
      read (idata,996) phrase
      read (idata,*) ischema
      print 998 ,  phrase(1:40) ,  ischema
      read (idata,996) phrase
      read (idata,*) joulin
      print 998 ,  phrase(1:40) ,  joulin
      read (idata,996) phrase
      read (idata,*) idrap
      print 998 ,  phrase(1:40) ,  idrap
      read (idata,996) phrase
      read (idata,*) ifonc
      print 998 ,  phrase(1:40) ,  ifonc
      read (idata,996) phrase
      read (idata,*) thf
      print 997 ,  phrase(1:40) ,  thf
      read (idata,996) phrase
      read (idata,*) thg
      print 997 ,  phrase(1:40) ,  thg
      read (idata,996) phrase
      read (idata,*) uinit
      print 997 ,  phrase(1:40) ,  uinit
      read (idata,996) phrase
      read (idata,*) np
      print 998 ,  phrase(1:40) ,  np
      read (idata,996) phrase
      read (idata,*) nmail
      print 998 ,  phrase(1:40) ,  nmail
      read (idata,996) phrase
      read (idata,*) tfin
      print 997 ,  phrase(1:40) ,  tfin
      read (idata,996) phrase
      read (idata,*) beta
      print 997 ,  phrase(1:40) ,  beta
      read (idata,996) phrase
      read (idata,*) eee
      print 997 ,  phrase(1:40) ,  eee
c
      return
 999  format (1x,a)
 998  format (1x,a,1x,i8)
 997  format (1x,a,1x,e13.6)
 996  format (a)
      end
c
c
c==========================================================================
c
      subroutine coeff ( ischema , npmax , dt , czero , cc , ccdd )
c     calcul des coefficients pour l'integrale de convolution
c     parametres d entree
c        ischema   choix du schema numerique
c               1  racine carree du schema a deux points
c               2  racine  carree du schema a trois points
c               3  schema de msallam
c               4  elements finis P1
c               5  methode mixte u dans P1, f(u) dans P0.
c               6  methode mixte u dans P1, f(u) dans P1.
c        npmax     nombre maximal de pas de temps
c        dt        valeur du pas de temps
c     parametre de sortie
c        czero     premier des coefficients cc pour l indice zero
c        cc        coefficients pour le calcul du second membre
c        ccdd      idem pour le gradient du flux
c
c=======================================================================
c     fd et sm     7 1 2K    23 01 2K   06-02-2K     26-06-2000
c               08 nov 2000 - 30 ans depuis le bal tragique a Colombey
c                 30 nov 2k.
c=======================================================================
c
      implicit double precision (a-h,o-z)
      dimension cc(*) , ccdd(*)
      pi=3.141592653589793238462643383D+00
c                                racine carree du schema a deux points
      if (ischema.eq.1) then
      czero    = 1.d0
      unsaa    = 1.d0 / czero
      cc(1)    = -unsaa/2
      cc(2)    = unsaa * (0.d0 - cc(1)*cc(1)) / 2.d0
      do 11  j = 3 , npmax
      jp       = j - 2*(j/2)
c                                  j est impair
      if (jp.eq.1) then
      i        = (j-1)/2
      c        = 0.d0
      do 12  k = 1 , i
      c        = c - cc(k) * cc(j-k)
 12   continue
      cc(j)    = unsaa * c
      else
c                                  j est pair
      i        = j / 2
      c        = - 0.5d0 * cc(i) * cc(i)
      do 13  k = 1 , i-1
      c        = c - cc(k) * cc(j-k)
 13   continue
      cc(j)    = unsaa * c
      end if
c      write (6,*) j  , 'cc  initial' , cc(j)
 11   continue
c
      do 14  j = 1 , npmax
      cc(j)    = -cc(j)
c      print * , '  cc  ' , j , cc(j)
 14   continue
      end if
c
c                                racine carree du schema a trois points
      if (ischema.eq.2) then
      czero    = sqrt(3.d0/2.d0)
      aa       = sqrt(3.d0/2.d0)
      unsaa    = 1.d0 / aa
      cc(1)    = -unsaa
      cc(2)    = unsaa * (0.5d0 - cc(1)*cc(1)) / 2.d0
c      write (6,*) 0  , 'cc  initial' , aa
c      write (6,*) 1  , 'cc  initial' , cc(1)
c      write (6,*) 2  , 'cc  initial' , cc(2)
      do 21  j = 3 , npmax
      jp       = j - 2*(j/2)
c
c     print * ,  j  , jp , j/2
c
c                                  j est impair
      if (jp.eq.1) then
      i        = (j-1)/2
      c        = 0.d0
      do 22  k = 1 , i
      c        = c - cc(k) * cc(j-k)
 22   continue
      cc(j)    = unsaa * c
      else
c                                  j est pair
      i        = j / 2
      c        = - 0.5d0 * cc(i) * cc(i)
      do 23  k = 1 , i-1
      c        = c - cc(k) * cc(j-k)
 23   continue
      cc(j)    = unsaa * c
      end if
c      write (6,*) j  , 'cc  initial' , cc(j)
 21   continue
c
      do 24  j = 1 , npmax
      cc(j)    = -cc(j)
c      write (6,*) j  , 'cc  ' , cc(j)
 24   continue
      end if
c
c         schema de Regis Msallam : calcul simplifie de la convolution
c
      if (ischema.eq.3) then
      czero    = 1.d0
      do 31  j = 1 , npmax
      cc(j)    = sqrt(dt/pi) * (sqrt(j+1.d0)-sqrt(j-1.d0))
  31  continue
      end if
c                        schema par elements finis conformes P1
      if (ischema.eq.4) then
      czero    = 1.d0
      do 41  j = 1 , npmax
      cc(j)    = 4.d0 * sqrt(dt/pi) * ((j+1.d0)*sqrt(j+1.d0)
     &                             - 2.d0 * (j+0.d0)*sqrt(j+0.d0)
     &                            + (j-1.d0)*sqrt(j-1.d0)) / 3.d0
c      print 999 , j , cc(j)
  41  continue
      end if
c                               methode mixte u dans P1, f(u) dans P0.
      if (ischema.eq.5) then
      czero    = 1.d0
      do 51  j = 1 , npmax
      cc(j)    = 2.d0 * sqrt(dt/pi) * (sqrt(j+1.d0)-sqrt(j+0.d0))
c      print 999 , j , cc(j)
  51  continue
      end if
c                               methode mixte u dans P1, f(u) dans P1.
      if (ischema.eq.6) then
      czero    = 1.d0
      do 61  j = 1 , npmax
      cc(j)    = 2.d0 * sqrt(dt/pi) * (sqrt(j+1.d0)-sqrt(j+0.d0))
  61  continue
      do 62  j = 1 , npmax
      ccdd(j)  =  sqrt(dt/pi)*(
     &            (2.d0*j+1.d0)*(sqrt(j+1.d0)-sqrt(j+0.d0))
     &          - 2.d0 *((j+1d0)*sqrt(j+1.d0) - (j+0.d0)*sqrt(j+0.d0))
     &                    / 3.d0  )
  62  continue
      end if
c
      return
 999  format (1x , i6 , 3X , e25.18)
      end
c
c=======================================================================
c
      subroutine temps (ischema , idrap , ifonc , thf , thg , uinit ,
     &                  np , dt ,czero , cc , ccdd , uu , phi , dphi )
c     boucle en temps du schema aux differences
c     parametres d entree
c        ischema   choix du schema numerique
c               1  racine carree du schema a deux points
c               2  racine  carree du schema a trois points
c               3  schema de msallam
c               4  elements finis P1
c               5  methode mixte u dans P1, f(u) dans P0
c               6  methode mixte u dans P1, f(u) dans P1.
c        idrap     choix de la fonction f(u)
c               0  f = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine  carree de u
c               3                        u puissance deux tiers
c               4                                    trois quarts
c               5                                    quatre cinquiemes
c               6      proportionnel a u
c        ifonc     choix de la fonction g(t)
c               0  g = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine carree de t
c               3      fonction lineaire
c               4      proportionnelle a t multiplie par racine de t
c               5      quadratique
c               6      harmonique               ref = Lubich 86
c               7      sinus de racine de t     ref = Lubich 86
c               8      polynome de racine de t  ref = Diethelm
c
c        thf       parametre multiplicatif de la fonction f
c        thg       parametre multiplicatif de la fonction g
c        uinit     condition initiale
c        np        nombre de pas de temps
c        dt        valeur du pas de temps
c        czero     premier des coefficients cc
c        cc        coefficients qui permettent le calcul du second membre
c     parametre de sortie
c        uu        vecteur des solutions discretes
c        phi       tableau des flux moyens dans un intervalle de temps
c        dphi      idem pour leur derivee en temps
c
c=======================================================================
c           fd et sm     7 1 2K   23-01-2K   06-02-2K   24 avril 2000
c                       26 06 2K   18 aout 2k    20 octobre 2000
c=======================================================================
c
      implicit double precision (a-h,o-z)
      common / joul1 / joulin , beta, eee
      dimension cc(*) , ccdd(*) , uu(*) , phi(*) , dphi(*)
      pi=3.141592653589793238462643383D+00
c
c     calcul du coefficient a en fonction du schema numerique
      if (ischema.eq.1) a = sqrt(dt)
      if (ischema.eq.2) a = sqrt(dt)
      if (ischema.eq.3) a = sqrt(dt/pi)
      if (ischema.eq.4) a = 4.d0 * sqrt(dt/pi) / 3.d0
      if (ischema.eq.5) a = sqrt(dt/pi)
      if (ischema.eq.6) a = sqrt(dt/pi) * (1.+1.d0/3.d0)
c
      azero   = czero
c
c     boucle en temps
      do 1  k = 1 , np
c                                                Calcul du premier point
      if (k.eq.1) then
c      print * , ' toto   ischema  = ' , ischema
      if (joulin.eq.0) then
      call  premier ( ischema , idrap , ifonc ,   thf , thg , uinit ,
     &                dt ,czero , uu , phi , dphi )
c      print * , ' apres premier ' , uu(1)
      else
      call newzer ( dt , u , flux )
      uu(1)   =  u
c      print * , ' u(1) = ' , u
      phi(1)  =  flux
      end if
      else
c               calcul du second membre ie de l'integrale de convolution
      b       = 0.d0
c                         Introduction du drift g(t) dependant du temps
c                                         ==========
c                         et aussi d autres termes de convolution ???
c
      if (ischema.le.2) then
      tk      = k * dt
      call gg ( ifonc, dt , tk , gk , zmug , zmugpr )
      gk      = thg * gk
      b       = b + sqrt(dt) * gk
      end if
c
      if ((ischema.eq.3).or.(ischema.eq.4)) then
      tzero   = 0.d0
      call gg ( ifonc, dt , tzero , gzero , zmug , zmugpr )
      gzero   = thg * gzero
      if (ischema.eq.3) then
      b       = b + sqrt(dt/pi) * (sqrt(k+0.d0)-sqrt(k-1.d0)) * gzero
      end if
      if (ischema.eq.4) then
      b       = b + 2.d0 * sqrt(dt/pi) *((3.d0 -2.d0*k)*sqrt(k+0.d0)
     &                   + 2.d0 *(k-1.d0)*sqrt(k-1.d0)) * gzero / 3.d0
      end if
      do 20  j = 1 , k-1
      tj       = j * dt
      call gg ( ifonc, dt , tj , gj , zmug , zmugpr )
      gj       = thg * gj
      b        = b + cc(k-j) * gj
  20  continue
      tk       = k * dt
      call gg ( ifonc, dt , tk , gk , zmug , zmugpr )
      gk       = thg * gk
      if (ischema.eq.3) b = b + sqrt(dt/pi) * gk
      if (ischema.eq.4) b = b + 4.d0 * sqrt(dt/pi) * gk / 3.d0
      end if
c
      if (ischema.eq.5) then
      if (joulin.eq.0) then
      do 50 j = 1 , k-1
c     print * , j , phi(j)
      b       = b + cc(k-j) * phi(j)
  50  continue
c      print * , ' b = ' , b
      tkm1    = (k-1.d0) * dt
      call gg ( ifonc, dt , tkm1 , gkm1 , zmug , zmugpr )
      b       = b + uinit + 2.d0 * sqrt(dt/pi) * thg * zmug
c      print * , ' b = ' , b
      else
c                                                   cas Joulin
      sm      = 0.d0
      do 51 j = 1 , k-1
c     print * , j , phi(j)
      sm      = sm + cc(k-j) * phi(j)
  51  continue
c      print * , ' second membre  ' , sm
      um1     =   uu(k-1)
      kk      =   k
      end if
      end if
c
      if (ischema.eq.6) then
      do 60 j = 1 , k-1
c     print * , j , phi(j)
      b       = b + cc(k-j) * phi(j) + ccdd(k-j) * dt * dphi(j)
  60  continue
c      print * , ' b = ' , b
      tkm1    = (k-1.d0) * dt
      call gg ( ifonc, dt , tkm1 , gkm1 , zmug , zmugpr )
      b       = b + uinit + 2.d0 * sqrt(dt/pi) * thg *
     &          ( zmug + dt*zmugpr / 3.d0 )
      end if
c                                        fin du terme de drift
c
c
c          debut du terme dynamique au sein du second membre
c                   ===============
c
      do 2  j = 1 , k-1
c     print *, 'k,uu(j)',k,uu(j)
      if (ischema.le.2)then
      b       = b + cc(k-j) * uu(j)
      end if
      if ((ischema.eq.3).or.(ischema.eq.4)) then
      u       = uu(j)
      du      = 1.d0
      call ff(idrap,du,u,f,fd,zmuf)
      f       = thf * f
      b       = b + cc(k-j) * f
      end if
c      print * , 'k,j,cc,f,uu(j)',k,j,cc(k-j),f,uu(j)
 2    continue
c
      if (ischema.eq.3) then
      u       = 0.d0
      du      = 1.d0
      call ff(idrap,du,u,f,fd,zmuf)
      f       = thf * f
      fd      = thf * fd
      b       = b + sqrt(dt/pi) * (sqrt(k+0.d0)-sqrt(k-1.d0)) * f
      end if
c
      if (ischema.eq.4) then
      u       = 0.d0
      du      = 1.d0
      call ff(idrap,du,u,f,fd,zmuf)
      f       = thf * f
      b       = b + 2.d0 * sqrt(dt/pi) *((3.d0 -2.d0*k)*sqrt(k+0.d0)
     &                   + 2.d0 *(k-1.d0)*sqrt(k-1.d0)) * f / 3.d0
      end if
c
      if (ischema.eq.5)then
c                formule de quadrature tres simple pour le calcul
c                de la moyenne de f  entre u(k-1) et u(k).
      if (joulin.eq.0) then
      ukm1    = uu(k-1)
      du      = 1.d0
      call ff(idrap,du,ukm1,fkm1,fd,zmuf)
c      print * , ' fkm1 = ' , fkm1
c      print * , ' db = ' , sqrt(dt/pi) * thf * fkm1
      b       = b + sqrt(dt/pi) * thf * fkm1
      end if
      end if
c
      if (ischema.eq.6)then
      ukm1    = uu(k-1)
      du      = 1.d0
      call ff(idrap,du,ukm1,fkm1,fd,zmuf)
      b       = b + sqrt(dt/pi) * thf * (1-1.d0/3.d0) * fkm1
      end if
c
c      print * ,'b=',b
c
c
c====================================================================
c     resolution de l equation       (1) azero  u - a thf f(u) = b
c
      uzero  = uu(k-1)
      if (joulin.eq.0) then
      call newton (idrap, thf , azero , a , b , uzero , u)
      uu(k)  = u
c      print * ,' u = ' , u
      else
      um1     =   uu(k-1)
      kk      =   k
      call newjoulin ( sm , dt , kk , um1 , u , flux)
      uu(k)   =  u
      phi(k)  =  flux
      end if
c
c                             incrementation du tableau phi(*)
      if (ischema.eq.5) then
      if (joulin.eq.0) then
      tkm1    = (k-1) * dt
      call gg ( ifonc, dt , tkm1 , gkm1 , zmugk , zmugkp )
      ukm1    = uu(k-1)
      du      = uu(k) - uu(k-1)
      call ff( idrap,du, ukm1, fkm1 ,fd,zmufk )
      uk      = uu(k)
      call ff( idrap,du, uk, fk,fd,zmufk )
      phi(k) = thf * 0.5d0 * (fkm1 + fk) + thg * zmugk
      end if
      end if
c                             incrementation de phi(*) et dphi(*)
      if (ischema.eq.6) then
      tkm1    = (k-1) * dt
      call gg ( ifonc, dt , tkm1 , gkm1 , zmugk , zmugkp )
      ukm1    = uu(k-1)
      du      = uu(k) - uu(k-1)
      call ff( idrap,du, ukm1, fkm1 ,fd,zmufk )
      uk      = uu(k)
      call ff( idrap,du, uk, fk,fd,zmufk )
      phi(k) = thf * 0.5d0 * (fkm1 + fk) + thg * zmugk
      dphi(k) = thf * 2.d0 * (-fkm1 + fk) + thg * zmugkp
c      print * , k , ' dphi(k) ' , dphi(k)
      end if
c
      end if
c
c      write(6,999)  'temps ',k, '  solution ' , uu(k)
c
c
  1   continue
      return
 999  format (1x , a , 1x , i7 , 1x, a ,  e13.6 )
      end
c
c=======================================================================
c
      subroutine premier ( ischema , idrap , ifonc ,  thf , thg ,
     &                     uinit ,  dt , czero , uu , phi , dphi )
c     calcul du premier point du schema aux differences
c     parametres d entree
c        ischema   choix du schema numerique
c               1  racine carree du schema a deux points
c               2  racine  carree du schema a trois points
c               3  schema de msallam
c               4  elements finis P1
c               5  methode mixte u dans P1, f(u) dans P0
c               6  methode mixte u dans P1, f(u) dans P1.
c        idrap     choix de la fonction f(u)
c               0  f = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine  carree de u
c               3                        u puissance deux tiers
c               4                                    trois quarts
c               5                                    quatre cinquiemes
c               6      proportionnel a u
c        ifonc     choix de la fonction g(t)
c               0  g = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine carree de t
c               3      fonction lineaire
c               4      proportionnelle a t multiplie par racine de t
c               5      quadratique
c               6      harmonique               ref = Lubich 86
c               7      sinus de racine de t     ref = Lubich 86
c               8      polynome de racine de t  ref = Diethelm
c        thf       parametre multiplicatif de la fonction f
c        thg       parametre multiplicatif de la fonction g
c        uinit     condition initiale
c        dt        valeur du pas de temps
c        czero     premier des coefficients cc
c    parametre de sortie
c        uu        vecteur des solutions discretes
c        phi       tableau des flux moyens dans un intervalle de temps
c        dphi      idem pour leur derivee en temps
c
c=======================================================================
c      fd et sm      23-01-2K    06-02-2K    27 fevrier 2000  24 04 2K
c                    06 mai 2000, 26 juin 2K  18 aout 2000   20 oct 2K
c=======================================================================
c
      implicit double precision (a-h,o-z)
      dimension   uu (*), phi(*) , dphi(*)
      pi=3.141592653589793238462643383D+00
c
c     calcul du coefficient a en fonction du schema numerique
      if (ischema.eq.1) a = sqrt(dt)
      if (ischema.eq.2) a = sqrt(dt)
      if (ischema.eq.3) a = sqrt(dt/pi)
      if (ischema.eq.4) a = 4.d0 * sqrt(dt/pi) / 3.d0
      if (ischema.eq.5) a =  sqrt(dt/pi)
      if (ischema.eq.6) a =  sqrt(dt/pi) * (1.+1.d0/3.d0)
c
      azero   = czero
c
c     calcul du second membre ie de l'integrale de convolution
c               =============
c
c     calcul des premieres valeurs de la fonction g(t) du temps
c                             Ici la fonction du temps imposee
      tzero   = 0.d0
      tun     = dt
      call gg ( ifonc, dt , tzero , gzero , zmug0 , zmugp0 )
      gzero   = thg * gzero
c      print * , ifonc , tzero , gzero
      call gg ( ifonc, dt , tun , gun , zmug1 , zmugp1 )
      gun     = thg * gun
c
      if (ischema.le.2) then
      b       = sqrt(dt) * gun
      end if
      if (ischema.eq.3) then
      b       = sqrt(dt/pi) * ( gzero + gun )
      end if
      if (ischema.eq.4) then
      b       = sqrt(dt/pi) * ( 2.d0 * gzero  + 4.d0 * gun )/3.d0
      end if
      if (ischema.eq.5) then
      b       = 2.d0 * sqrt(dt/pi) * thg * zmug0
      end if
      if (ischema.eq.6) then
      b       = sqrt(dt/pi) * thg * ( 2.d0 * zmug0 + dt*zmugp0 / 3.d0 )
      end if
c                                      Ici commence la dynamique
      if (ischema.eq.3) then
      u       = 0.d0
      du      = 1.d0
      call ff(idrap,du,u,f,fd,zmuf)
      f       = thf * f
      b       = b + sqrt(dt/pi) * f
      end if
c
      if (ischema.eq.4) then
      u       = 0.d0
      du      = 1.d0
      call ff(idrap,du,u,f,fd,zmuf)
      f       = thf * f
      k       = 1
      k       = 1.d0
      b       = b + 2.*sqrt(dt/pi) *((3.d0-2.d0*k)*sqrt(k+0.d0)
     &              + 2.d0 * (k-1.d0)*sqrt(k-1.d0)) * f / 3.d0
      end if
c
      if (ischema.eq.5) then
      uzero   = uinit
      du      = 1.d0
      call ff(idrap,du,uzero,fzero,fd,zmuf)
      b       = b +  uinit + sqrt(dt/pi) * thf * fzero
      end if
c
      if (ischema.eq.6) then
      uzero   = uinit
      du      = 1.d0
      call ff(idrap,du,uzero,fzero,fd,zmuf)
      b  = b +  uinit + sqrt(dt/pi) * thf * (1-1.d0/3.d0) * fzero
      end if
c      print * , ' premier    b  =   ' , b

c
c====================================================================
c                             Inialisation de l algorithme de newton
      if (idrap.eq.1)  then
c                             sois zen !  on resout l equation ici...
      u      = 0.d0
      du     = 1.d0
      call ff(idrap,du,u,f,fd,zmuf)
      f      =  thf * f
      uzero  = (b+a*f)/azero
      end if
c
      if (idrap.eq.2) then
c                    sois zen !  f est de la forme aaa * racine (u)
c                    donc aaa = f(1) zen sois !
      u      = 1.d0
      du     = 1.d0
      call ff(idrap,du,u, aaa ,fd,zmuf)
      uuzz1  =  thf * thf  * (aaa/azero) * (aaa/azero)
      uuzz2  =  b / azero
      uzero  =  0.5 * ( uuzz1 + uuzz2 )
c      uzero  = uuzz1
      end if
c
      if (idrap.eq.3) then
c                    sois zen !  idem plus haut mais cette fois
c                    f est de la forme aaa * u puissance deux tiers
      u      = 1.d0
      du     = 1.d0
      call ff(idrap,du,u, aaa ,fd,zmuf)
      uuzz1  =  thf * thf  *thf * (aaa / azero) ** 3
      uuzz2  =  b / azero
      uzero  =  0.5 * ( uuzz1 + uuzz2 )
      end if
c
      if (idrap.eq.4) then
c                    f est de la forme aaa * u puissance trois quarts.
      u      = 1.d0
      du     = 1.d0
      call ff(idrap,du,u, aaa ,fd,zmuf)
      uuzz1  = (thf * aaa / azero) ** 4
      uuzz2  =  b / azero
      uzero  =  0.5 * ( uuzz1 + uuzz2 )
      end if
c
      if (idrap.eq.5) then
c               f est de la forme aaa * u puissance quatre cinquiemes
      u      = 1.d0
      du     = 1.d0
      call ff(idrap,du,u, aaa ,fd,zmuf)
      uuzz1  = (thf * aaa / azero) ** 5
      uuzz2  =  b / azero
      uzero  =  0.5 * ( uuzz1 + uuzz2 )
      end if
c
      if (idrap.eq.6) then
c                                          f est de la forme aaa * u
      uzero  =  b / (azero - a*thf)
      end if
c---------------------------------------------------------------------
      call   newton ( idrap, thf , azero , a , b , uzero , u )
c
c---------------------------------------------------------------------
      uu(1)  = u
c
      if (ischema.eq.5) then
      tzero   = 0.d0
      call gg ( ifonc, dt , tzero , gzero , zmug0 , zmugp0 )
      uzero   = uinit
      du      = uu(1) - uzero
      call ff(idrap,du,uzero, fzero ,fd,zmuf0)
      uun     = uu(1)
      call ff(idrap,du,uun, fun ,fd,zmuf0)
c
c           attention : calcul des moyennes de deux facons differentes...
c
      phi(1)  = thf * 0.5d0 * (fzero + fun)  + thg * zmug0
      end if
c
      if (ischema.eq.6) then
      tzero   = 0.d0
      call gg ( ifonc, dt , tzero , gzero , zmug0 , zmugp0 )
      uzero   = uinit
      du      = uu(1) - uzero
      call ff(idrap,du,uzero, fzero ,fd,zmuf0)
      uun     = uu(1)
      call ff(idrap,du,uun, fun ,fd,zmuf0)
      phi(1)  = thf * 0.5d0 * (fzero + fun)  + thg * zmug0
      phi(1)  = thf * 0.5d0 * (fzero + fun)  + thg * zmug0
      end if
c
      return
      end
c
c=======================================================================
c
      subroutine newzer ( dt , u , flux )
c                                                    Cas test de Joulin
c     resolution de l equation  (1) au premier pas de temps
c
c          I    u = 2 sqrt(dt/pi) * ( log(u) - 1 +
c     (1)  I              + eee *dt^beta * (beta + 1 - beta * dt )/
c          I                           (beta * (beta+1) * u ) )
c
c     parametres d entree
c        beta      exposant de la dynamique en temps
c        eee       constante de Audounet et Roquejoffre
c        dt        valeur du pas de temps
c     parametre de sortie
c        u         solution de l equation (1)
c        flux      flux de la solution
c
c=======================================================================
c        fd et sm     7 01  2K     23 1 2K     06-02-2k     24-04-2000
c                     6 mai 2000    6 juin 2K   14 septembre 2000
c                     21 dec 2000.
c=======================================================================
c
      implicit double precision (a-h,o-z)
      common / joul1 / joulin, beta, eee
      pi=3.141592653589793238462643383D+00
      nie    =  100
c                                                   initialisation
c      print 999 , ' beta ' , beta
c      print 999 , ' eee  ' , eee
c      print 999 , '   dt ' , dt

      a      =  2.d0 * sqrt(dt/pi)
      b      =  eee * dt**beta * ( beta + 1.d0 - beta*dt ) /
     &                           (beta * ( beta + 1.d0 ) )
c
c      print 999 , '    a  ' , a
c      print 999 , '    b  ' , b
      u      =  sqrt ( b )
c      print *,' uzero  =  ',  u
c
      do 1 i = 1 , nie
      f      = u - a * ( log(u)-1.d0 ) - a * b / u
      fd     = 1.d0  - a / u + a*b / (u*u)
      du     = - f / fd
c      write(6,*)  '   numero   u  f   fd  du  '  ,  i , u , f , fd, du
      if (abs(du).ge.abs(u)/10.d0) then
      v      = u + 0.1d0 * du
c      write(6,*)  '  hello '
      else
      v      = u +  du
      end if
c
      if (u*v.le.0) then
      print * , ' Danger a la ' , i , '  ieme iteration du newton'
      du = 0.1 * du
      end if
      u      = v
c      write (6,*) 'Newton n°',i , '  du =  ',du  ,'  u =', u
 1    continue
c
      flux   = log(u)-1.d0 + b/u
      return
 999  format (1x,a,1x,e20.13)
      end
c
c=======================================================================
c
      subroutine newjoulin (sm , dt , k , um1 , u , flux)
c                                                  Newton pour Joulin
c     resolution de l equation   pas de temps
c     parametres d entree
c        sm        morceau de l integrale de convolution
c        beta      exposant de la dynamique en temps
c        eee       constante de Audounet et Roquejoffre
c        dt        valeur du pas de temps
c        k         numero du pas de temps
c        um1       solution au pas de temps d avant
c     parametre de sortie
c        u         solution de l equation (1)
c        flux      flux de la solution
c
c=======================================================================
c        fd et sm     7 01  2K     23 1 2K     06-02-2k     24-04-2000
c                     6 mai 2000    6 juin 2K   14 septembre 2000,
c                     04 octobre 2000, 21 dec 2000.
c=======================================================================
c
      implicit double precision (a-h,o-z)
      common / joul1 / joulin, beta, eee
      pi=3.141592653589793238462643383D+00
c      print * , sm , beta , eee , dt , k , um1
      nie    =  100
c                                                   initialisation
      a      =  2.d0 * sqrt(dt/pi)
      u      =  um1
c
      do 1 i = 1 , nie
      if (i.eq.1) then
      phi1   = log(u)
      dphi1  = 1.d0 / u
      else
      phi1   = log(u) - 1.d0 + um1 * ( log(u) - log(um1) ) / (u - um1)
      dphi1  = 1.d0 / u + um1 / ( u * (u - um1) ) -
     &             um1 * ( log(u) - log(um1) ) / ((u - um1)*(u - um1))
      end if
c
      tt = (k-0.d0)*dt
      if (tt.le.1.d0) then
      phi2   = eee *  dt**beta *
     &         ( (k-1.d0)**beta * (1.d0 -(k-1.d0)*dt) / um1   +
     &         8.d0 *(k-0.5d0)**beta * (1.d0 -(k-0.5d0)*dt)/(um1+u) +
     &          ( (k-0.d0)**beta * (1.d0 -(k-0.d0)*dt) ) / u  )/6.d0
      dphi2  = -  eee *  dt**beta * (
     &  8.d0 *(k-0.5d0)**beta * (1.d0 -(k-0.5d0)*dt)/((um1+u)*(um1+u))
     &   +  ((k-0.d0)**beta ) *  (1.d0 -(k-0.d0)*dt) / (u*u)   )/6.d0
      else
      phi2   = 0.d0
      dphi2  = -0.d0
      end if
c
      f      = u - a * ( phi1 + phi2 ) - sm
      df     = 1.d0 - a *  ( dphi1 + dphi2 )
      du     = - f / df
c      write(6,*)  'numero   u    f   df  du  '  ,  i , u , f , df, du
      if (abs(du).ge.abs(u)/10.d0) then
      v      = u + 0.1d0 * du
c      write(6,*)  '  hello '
      else
      v      = u +  du
      end if
c
      if (u*v.le.0) then
      print * , ' Danger a la ' , i , '  ieme iteration du newton'
      du = 0.1 * du
      end if
      u      = v
c      write (6,*) 'Newton n°',i , '  du =  ',du  ,'  u =', u
 1    continue
c
      flux   = phi1 + phi2
      return
      end
c
c=======================================================================
c
      subroutine newton ( idrap, thf , azero , a , b , uzero , u )
c     resolution de l equation     (1)       azero u - a thf f(u) = b
c     parametres d entree
c        idrap     choix de la fonction f
c               0  f = zero
c               1      constante
c               2      fonction racine carree
c               3      proportionnelle a u puissance trois quarts
c               4                                    deux tiers
c               5                                    quatre cinquiemes
c               6      proportionnel a u
c        thf       parametre multiplicatif de la fonction f
c        azero     coefficient dans l equation ci-dessus
c        a         coefficient dans l equation ci-dessus
c        b         second membre
c        uzero     condition initiale
c     parametre de sortie
c        u         solution de l equation (1)
c
c=======================================================================
c        fd et sm     7 01  2K     23 1 2K     06-02-2k     24-04-2000
c                     6 mai 2000
c=======================================================================
c
      implicit double precision (a-h,o-z)
      nie    =  300
      u      = uzero
c      print *,'uzero  =  ',uzero
      do 1 i = 1 , nie
      dduu   = 1.
      call ff ( idrap ,  dduu , u , f , fd , zmuf )
      f      = thf * f
      fd     = thf * fd
      du     = (b-azero*u+a*f)/(azero-a*fd)
c      write(6,*)  'u,azero,a,b,f,fd',u,azero,a,b,f,fd
c                     rajouter ici un test sur la valeur absolue de du
      v      = u + du
      if (u*v.le.0) then
      print * , ' Danger a la ' , i , '  ieme iteration du newton'
      du = 0.1 * du
      end if
      u      = u + du
c     write (6,*) 'Newton n°',i , '  du =  ',du  ,'  u =', u
 1    continue
      return
      end
c
c=======================================================================
c
      subroutine ff ( idrap, du , u , f , fd , zmuf )
c     calcul de f(u) et df sur dx (u) pour des fonctions f parametrees
c                                                           par idrap
c     parametres d entree
c        idrap     choix de la fonction f(u)
c               0  f = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine  carree de u
c               3                        u puissance deux tiers
c               4                                    trois quarts
c               5                                    quatre cinquiemes
c               6      proportionnel a u
c        du        pas de u donne pour faire la moyenne
c        u         argument
c     parametres de sortie
c        f         f(u)
c        fd        df sur dx (u)
c        zmuf      moyenne de f entre u et u+du
c
c
c=======================================================================
c         fd et sm     7 1 2K   06-02-2K  27 fevrier 2000  24-04-2K
c                      06 mai 2000     26 juin 2K
c=======================================================================
c
      implicit double precision (a-h,o-z)
      pi=3.141592653589793238462643383D+00
c
      if (idrap.eq.0) then
            f  =  0.d0
            fd =  0.d0
          zmuf =  0.d0
      end if
c
      if (idrap.eq.1) then
            f  =  0.5d0 * sqrt (pi)
            fd =  0.d0
          zmuf =  0.5d0 * sqrt (pi)
      end if
c
      if (idrap.eq.2) then
            f  =  2.d0 * sqrt (u/pi)
            fd =  1.d0 / sqrt(u*pi)
          zmuf =  2.d0 * sqrt (1.d0/pi) *((u+du)**1.5d0 - u**1.5d0) /
     &                                         du
      end if
c
      if (idrap.eq.3) then
      deutier  = 2.d0/3.d0
      cinqtie  = 5.d0 / 3.d0
            f  =  3.d0 * sqrt(pi) * ( u**deutier ) / 4.d0
            fd =  0.5d0 * sqrt(pi) /  u **(1.d0/3.d0)
          zmuf =  9.d0 * sqrt(pi) * ((u+du)**cinqtie - u**cinqtie)
     &                                 / (20.d0 *du)
      end if
c
      if (idrap.eq.4) then
      troiqu  = 3.d0 / 4.d0
      septqu  = 7.d0 / 4.d0
            f  =  8.d0 * ( u**troiqu ) / (3.d0 * sqrt(pi) )
            fd =  2.d0 / ( sqrt(pi) * (  u **(1.d0/4.d0)))
          zmuf =  32.d0 * ((u+du)**septqu - u**septqu )
     &                      / (21.d0 * sqrt(pi) * du )
      end if
c
      if (idrap.eq.5) then
      quacin   = 4.d0 / 5.d0
      neucin   = 9.d0 / 5.d0
            f  =  15.d0 * sqrt(pi) * ( u**quacin ) / 16.d0
            fd =  0.75d0 * sqrt(pi) /  u **(1.d0/5.d0)
          zmuf =  5.d0 * 15.d0 * sqrt(pi) * ((u+du)**neucin
     &           - u**neucin)          / (16.d0 * 7.d0 * du)
      end if
c
      if (idrap.eq.6) then
            f  =  u
            fd =  1.d0
      end if
c      print *, 'idrap,u, f et fd',idrap, u,f,fd
c
      return
      end

c
c=======================================================================
c
      subroutine gg (ifonc, dt , t , g , zmug , zmugpr )
c     calcul de  g(t) pour des fonctions g parametrees
c     parametres d entree
c        ifonc     choix de la fonction g(t)
c               0  g = zero
c               1      constante egale a 1/2 de racine de pi
c               2      proportionnelle a la racine carree de t
c               3      fonction lineaire
c               4      proportionnelle a t multiplie par racine de t
c               5      quadratique
c               6      harmonique               ref = Lubich 86
c               7      sinus de racine de t     ref = Lubich 86
c               8      polynome de racine de t  ref = Diethelm
c        dt        intervalle de temps sur lequel on moyenne  g
c        t         argument
c     parametre  de sortie
c        g         g(t)
c        zmug      moyenne de g entre g(t) et g(d+dt)
c        zmugpr    gradient de g au point milieu de t a t+dt
c
c=======================================================================
c             fd et sm    06-02-2K  27 fevrier 2000  24 avril 2000
c                         26 juin 2000   18 aout 2000  20 oct 2000
c                         08 nov 2K
c=======================================================================
c
      implicit double precision (a-h,o-z)
      pi=3.141592653589793238462643383D+00
c
      if (ifonc.eq.0) then
            g  = 0.
d0         zmug  = 0.d0
       zmugpr  = 0.d0
      end if
c
      if (ifonc.eq.1) then
            g  =  0.5d0 * sqrt (pi)
         zmug  =  0.5d0 * sqrt (pi)
       zmugpr  =  0.d0
      end if
c
      if (ifonc.eq.2) then
      g     =  2.d0 * sqrt (t/pi)
      zmug  =  4.d0 *((t+dt)**1.5d0 - t**1.5d0) /(3.d0 * sqrt(pi) * dt)
      zmugpr    =  sqrt((t+0.5d0*dt)*pi)
      end if
c
      if (ifonc.eq.3) then
            g  =  3.d0 * sqrt (pi) * t / 4.d0
         zmug  =  3.d0 * sqrt (pi) *((t+dt)**2 - t**2)/(8.d0 * dt)
       zmugpr  =  3.d0 * sqrt (pi)  / 4.d0
      end if
c
      if (ifonc.eq.4) then
            g  =  8.d0 * t *sqrt(t) / ( 3.d0 * sqrt(pi) )
         zmug  =  16.d0 *((t+dt)**2.5d0 - t**2.5d0 )
     &                                 /(15.d0 * sqrt(pi) * dt)
       zmugpr  =  4.d0 * sqrt(t) / ( sqrt(pi) )
      end if
c
      if (ifonc.eq.5) then
            g  =  15.d0 * sqrt(pi) * t*t / 16.d0
         zmug  =  5.d0 * sqrt (pi) *( (t+dt)** 3 - t**3 )
     &                                 /(16.d0 *  dt)
       zmugpr  =  15.d0 * sqrt(pi) * t / 8.d0
      end if
c
      if (ifonc.eq.6) then
            g  =  t / (1.d0 + t)
         zmug  =  1.d0 - log((1.d0+t+dt)/(1.d0+t)) / dt
       zmugpr  =  1.d0 / (1.d0 + t + 0.5d0*dt)**2
      end if
c
      if (ifonc.eq.7) then
            g  =  sin(sqrt(t)) / sqrt(pi)
c                        calcul approche pour ce cas delicat
         zmug  =  ( sin(sqrt(t)) + 4*sin(sqrt(t+0.5*dt)) +
     &              sin(sqrt(t+dt)) ) / (6.d0 *  sqrt(pi) )
       zmugpr  =  cos (sqrt(t+0.5d0*dt)) /(2.d0 * sqrt((t+0.5d0*dt)*pi))
      end if
c
      if (ifonc.eq.8) then
            g  = t*t +  8.d0 * t *sqrt(t) / ( 3.d0 * sqrt(pi) )
c                                                        modif 25 juin
         zmug  = ( (t+dt)** 3 - t**3 ) /(3.d0 *  dt)  
c                                                        modif 25 juin
     &           + 16.d0 *((t+dt)**2.5d0 - t**2.5d0 )
     &                                 /(15.d0 * sqrt(pi) * dt)
       zmugpr  = 2*(t+0.5d0*dt) + 4.d0 * sqrt(t) / ( sqrt(pi) )
      end if
c      print *, 'ifonc  ' ,  t , g, zmug , zmugpr c
      return
      end

   
