c balot.f: version 3

c Ce programme simule par volumes finis un écoulement
c d'air et d'eau en présence de pesanteur.
c Les deux phases obéissent à la loi des
c gaz raides. L'équation d'énergie est prise en compte.
c Le domaine de calcul est rectangulaire. On peut retirer
c certains volumes pour simuler un récif.

c234567
      program balot

      implicit double precision (a-h,o-z)

c constantes du calcul:
c n: nombre de vf horizontaux
c m: nombre de vf verticaux

c init = 1 -->  on continue un calcul stocké dans "sauv"
c init = 0 -->  on utilise la condition
c initiale de "plotuv.dat"
c Dans ce dernier cas, il faut faire:
c "rm plotuv.dat"
c "g77 -o balot balot.f ; balot"
c "g77 -o soliton soliton.f ; soliton"
c "balot"

c nbvol: nombre de vf
c nbsom: nombre de sommets
c nbaret: nombre d'arêtes
c xmin,xmax:dimensions en x
c ymin,ymin:dimensions en y

c tmax: durée de la simulation
c l'instant initial est lu dans la 
c sauvegarde si init = 1

c cfl: cfl
c g: accélération de la pesanteur
c iordre: ordre de la méthode
c imiroir = -1 --> obstacle
c imiroir =  1 --> on recopie (debugage)
c raf : fréquence d'affichage (en secondes)
c irecif = 1 --> calcul avec un recif
c xrecif: abscisse du récif
c hrecif: hauteur du récif

      parameter (n=100,m=20,nbvol=n*m,nbsom=(n+1)*(m+1),init=0,
     & raf=10,xrecif=3.d0,hrecif=0.263d0,irecif=0,
     & nbaret=n*(m+1)+(n+1)*m,tmax=0.25d0,iordre=2,imiroir=-1,
     & cfl=2.d0,g=9.81d0,xmin=-2.d0,xmax=4.5d0,ymin=-.31d0,ymax=.5d0)

c tableaux
c xs,ys: coordonnées des sommets
c iaret: tableaux des arêtes
c iaret(1,ia) et le n° du premier volume voisin de ia
c iaret(2,ia) et le n° du second volume voisin de ia
c si iaret(2,ia)=0, ia est une arête de bord
c iext: permet de construire le maillage (inutile après)
c along, xnorm, ynorm: tableaux des longueurs et coordonnées
c des  normales aux arêtes
c ivol(ii,iv): n° du sommet ii (ii=1..4) du volume iv
c vol(iv): volume du vf n° iv
c perim(iv): périmètre du vf n° iv
c wn et wnp1: stockage de w à l'étape n et n+1
c fx, fy, flux: stockage des flux
c idur : repérage des volumes dans le récif

      dimension xs(nbsom),ys(nbsom),iaret(2,nbaret),
     &  ivol(4,nbvol),iext(2,nbaret),along(nbaret),
     & xnorm(nbaret),ynorm(nbaret),vol(nbvol),
     & wn(5,nbvol),wnp1(5,nbvol),perim(nbvol),
     & fx(5),fy(5),flux(5),wdemi(5,nbvol)
      real*4 tarray(2)

c tableaux utiles pour la montée en ordre
c dérivées suivant x et y de rho, u, v, p, yv

      dimension dxr(nbvol),dyr(nbvol),
     & dxu(nbvol),dyu(nbvol),
     & dxv(nbvol),dyv(nbvol),
     &     dxp(nbvol),dyp(nbvol),
     &     dxyv(nbvol),dyyv(nbvol),
     &     idur(0:nbvol)
      
      common/iter/iterriem/err/errriem
      
      iterriem=0
      errriem=0.d0

c     variables de base

      write(*,*) 'nb de volumes=',nbvol
      write(*,*) 'nb de sommets=',nbsom
      write(*,*) 'nb d''arêtes=',nbaret

c     pas en x et en y

      dx=(xmax-xmin)/n
      dy=(ymax-ymin)/m

c     remplissage des tableaux de connectivité

      write(*,*) 'construction du maillage'

c     les sommets

      do i=0,n
         do j=0,m
            is=i+j*(n+1)+1
            xs(is)=i*dx+xmin
            ys(is)=j*dy+ymin
         end do
      end do


c     les volumes

      do i=0,n-1
         do j=0,m-1
            iv=i+j*n+1
            ivol(1,iv)=i+j*(n+1)+1
            ivol(2,iv)=i+1+j*(n+1)+1
            ivol(3,iv)=i+1+(j+1)*(n+1)+1
            ivol(4,iv)=i+(j+1)*(n+1)+1
            vol(iv)=dx*dy
            perim(iv)=2*dx+2*dy
         end do
      end do

c     les arêtes horizontales

      icomptar=0

      do i=0,n-1
         do j=0,m
            ia=i+j*n+1
            iaret(1,ia)=i+j*n+1
            iaret(2,ia)=i+(j-1)*n+1
            iext(1,ia)=1
            iext(2,ia)=3

c     correction en haut et en bas pour que
c     le vecteur normal soit toujours sortant
            if (j.eq.m) then
               iaret(1,ia)=i+(j-1)*n+1
               iext(1,ia)=3
               iext(2,ia)=0
               iaret(2,ia)=0
            endif
            if (j.eq.0) then
               iaret(2,ia)=0
               iext(2,ia)=0
            endif 

            icomptar=icomptar+1
         end do
      end do

c     les arêtes verticales

      do i=0,n
         do j=0,m-1
            ia=i+j*(n+1)+1+icomptar
            iaret(1,ia)=i-1+j*n+1
            iaret(2,ia)=i+j*n+1
            iext(1,ia)=2
            iext(2,ia)=4
c     correction à gauche et à droite pour que
c     le vecteur normal soit toujours sortant
            if (i.eq.0) then
               iaret(1,ia)=i+j*n+1
               iext(1,ia)=4
               iext(2,ia)=0
               iaret(2,ia)=0
            endif
            if (i.eq.n) then
               iaret(2,ia)=0
               iext(2,ia)=0
            endif 
         end do
      end do


c     calcul des variables d'arêtes:
c     longueur, vecteurs normaux
      do ia=1,nbaret
         iv=iaret(1,ia)
         iloc2=iext(1,ia)
         iloc1=iloc2+1
         if(iloc1.eq.5) iloc1=1
         is1=ivol(iloc1,iv)
         is2=ivol(iloc2,iv)
         x1=xs(is1)
         y1=ys(is1)
         x2=xs(is2)
         y2=ys(is2)
         xmil=0.5*(x1+x2)
         ymil=0.5*(y1+y2)
         dlong=dsqrt((x1-x2)**2+(y1-y2)**2)
         along(ia)=dlong
         dxx=(x2-x1)/dlong
         dyy=(y2-y1)/dlong
         xn=-dyy
         yn=dxx
         xnorm(ia)=xn
         ynorm(ia)=yn
      enddo

c repérage des volumes dont le centre est dans
c le récif ( x > xrecif et y-ymin < hrecif )

      do iv=1,nbvol
         idur(iv)=1
      enddo
      idur(0)=-1

      if (irecif.eq.1) then
         do iv=1,nbvol
            xmil=0
            ymil=0
            do is=1,4
               xmil=xmil+xs(ivol(is,iv))
               ymil=ymil+ys(ivol(is,iv))
            enddo
            xmil=.25d0*xmil
            ymil=.25d0*ymil
            if ((xmil.ge.xrecif).and.((ymil-ymin).le.hrecif)) then
               idur(iv)=-1
            endif
         enddo
      endif
      
c les arêtes à l'interface entre l'eau et le récif
c sont marquées comme "miroir"

      if (irecif.eq.1) then
         do ia=1,nbaret
            if (idur(iaret(1,ia))*idur(iaret(2,ia)).lt.0) then
               iaret(2,ia)=0
            endif
         enddo
      endif

c     dessin du maillage
c     pour visualiser le maillage et les arêtes dans gnuplot
c     taper: plot 'mail' with line
      if (init.eq.0) then
         open(1,file='mail')
         do iv=1,nbvol
            write(1,*) xs(ivol(1,iv)),ys(ivol(1,iv))
            write(1,*) xs(ivol(2,iv)),ys(ivol(2,iv))
            write(1,*) xs(ivol(3,iv)),ys(ivol(3,iv))
            write(1,*) xs(ivol(4,iv)),ys(ivol(4,iv))
            write(1,*) xs(ivol(1,iv)),ys(ivol(1,iv))
            write(1,*)
         end do
         do ia=1,nbaret
            iv=iaret(1,ia)
            iloc2=iext(1,ia)
            iloc1=iloc2+1
            if(iloc1.eq.5) iloc1=1
            is1=ivol(iloc1,iv)
            is2=ivol(iloc2,iv)
            x1=xs(is1)
            y1=ys(is1)
            x2=xs(is2)
            y2=ys(is2)
            xmil=0.5*(x1+x2)
            ymil=0.5*(y1+y2)
            coef=.25*dsqrt(dx*dx+dy*dy)*.5
            write(1,*) xmil,ymil
            write(1,*) xmil+xnorm(ia)*coef,ymil+ynorm(ia)*coef
            write(1,*)
         enddo
         close(1)
      endif

c     initialisation des variables
      write(*,*) 'initialisations'

c     écriture de la liste des centres des volumes
c     fichier d'entrée du programme "soliton.f"
      if (init.eq.0) then
         write(*,*) 'écriture du maillage'
         write(*,*) 'il faut peut-être relancer "soliton.f"'
         write(*,*) 'pour actualiser la condition  initiale'
         open(88,file='pcal.dat')
         write(88,*) nbvol
c     calcul du centre du volume
         do iv=1,nbvol
            xmil=0
            ymil=0
            do ii=1,4
               xmil=xmil+xs(ivol(ii,iv))
               ymil=ymil+ys(ivol(ii,iv))
            enddo
            xmil=xmil/4
            ymil=ymil/4
            
            write(88,*) iv,xmil,ymil
            
         enddo
         close(88)
      endif
      
c     lecture des valeurs du fichier plotuv.dat
c     si init = 0
      
      if (init.eq.0) then
         write(*,*) 'lecture de "plotuv.dat"'
         open(89,file='plotuv.dat')
         do iv=1,nbvol
            read(89,*) iu,xmil,ymil,r,u,v,xp
            yv=0.d0
            if (iu.ne.iv) stop
            if (r.gt.500) yv=1.d0
            call conservative(wnp1(1,iv),r,u,v,xp,yv)
         enddo
         tini=0
         close(89)
      else
c     si init = 1, lecture en vrac dans "sauv"
         write(*,*) 'lecture de la sauvegarde binaire'
         open(97,form='unformatted',file='sauv')
         read(97) ni,mi
         if (ni.ne.n.or.mi.ne.m) then
            write(*,*) 'erreur de dimension'
            stop
         endif
         read(97) tini
         read(97) is
         read(97)ir
         if (is.ne.iordre) then
            write(*,*) 'les ordres ne correspondent pas'
            stop
         endif
         if (ir.ne.irecif) then
            write(*,*) 'problème de récif'
            stop
         endif
         read(97) ((wnp1(ii,i),ii=1,5),i=1,nbvol)
         close(97)
      endif
      
c     initialisation des pentes à 0

      do iv=1,nbvol
         dxr(iv)=0.d0
         dyr(iv)=0.d0
         dxu(iv)=0.d0
         dyu(iv)=0.d0
         dxv(iv)=0.d0
         dyv(iv)=0.d0
         dxp(iv)=0.d0
         dyp(iv)=0.d0
         dxyv(iv)=0.d0
         dyyv(iv)=0.d0
      end do

      if (irecif.eq.1) then
         write(*,*) 'récif présent'
      else
         write(*,*) 'pas de récif'
      endif

      write(*,*) 'calcul à l''ordre',iordre

c     boucle en temps
      write(*,*) 'instant initial=',tini
      write(*,*) 'instant final prévu=',tini+tmax
c     pas de temps
      dt=0
c     nbre d'itération
      it=0

      write(*,*) 'itérations en temps'
      t=tini
      cpu0=etime(tarray)
      cpu1=cpu0
      cpu2=cpu0

      do while(t-tini.lt.tmax)
         it=it+1

c     initialisation
c     on copie w(n+1) dans w(n)
c     et dans wdemi
         do iv=1,nbvol
            do ii=1,5
               wn(ii,iv)=wnp1(ii,iv)
               wdemi(ii,iv)=wnp1(ii,iv)
            end do
         end do

c     predicteur : ic=2
c     correcteur : ic=1
c pour l'ordre 2 en temps, on fait:

c wnp1 += dt/2 * f(wn)
c wn = wnp1
c wnp1= wdemi
c wnp1 += dt * f(wn)

c pour l'ordre 1, on fait juste:
c wnp1 += dt * f(wn)

         do ic=iordre,1,-1

c     ----------ordre 2--------------------------------------------
c     calcul des pentes par muscl
c     si méthode d'ordre 2

            if (iordre.eq.2) then

c     calcul des pentes en x et y
               do i=0,n-1
                  do j=0,m-1
c     volume central
                     iv=i+j*n+1
                     call primitive(wn(1,iv),r,u,v,p,yv)
c     volume ouest
                     if (i.ne.0) then
                        ivo=(i-1)+j*n+1
                        call primitive(wn(1,ivo),ro,uo,vo,po,yvo)
                     else
                        ro=r
                        uo=u*imiroir
                        vo=v
                        po=p
                        yvo=yv
                     endif
c     volume est
                     if (i.ne.n-1) then
                        ive=(i+1)+j*n+1
                        call primitive(wn(1,ive),re,ue,ve,pe,yve)
                     else
                        re=r
                        ue=u*imiroir
                        ve=v
                        pe=p
                        yve=yv
                     endif

c     volume nord
                     if (j.ne.m-1) then
                        ivn=i+(j+1)*n+1
                        call primitive(wn(1,ivn),rn,un,vn,pn,yvn)
                     else
                        rn=r
                        un=u
                        vn=v*imiroir
                        pn=p
                        yvn=yv
                     endif

c     volume sud
                     if (j.ne.0) then
                        ivs=i+(j-1)*n+1
                        call primitive(wn(1,ivs),rs,us,vs,ps,yvs)
                     else
                        rs=r
                        us=u
                        vs=v*imiroir
                        ps=p
                        yvs=yv
                     endif

c     calcul des pentes en x
                     call minmod(dx, ro, r, re, dxr(iv))
                     call minmod(dx, uo, u, ue, dxu(iv))
                     call minmod(dx, vo, v, ve, dxv(iv))
                     call minmod(dx, po, p, pe, dxp(iv))
                     call minmod(dx,yvo,yv,yve, dxyv(iv))

c     calcul des pentes en y
                     call minmod(dy, rs, r, rn, dyr(iv))
                     call minmod(dy, us, u, un, dyu(iv))
                     call minmod(dy, vs, v, vn, dyv(iv))
                     call minmod(dy, ps, p, pn, dyp(iv))
                     call minmod(dy,yvs,yv,yvn, dyyv(iv))

                  enddo
               enddo

            endif
c     ----------fin ordre 2--------------------------------------------


c     calcul du pas de temps associé à la cfl
c     nécessaire pour la stabilité du schéma
            if (ic.eq.iordre) then
               do iv=1,nbvol
                  call primitive(wn(1,iv),r,u,v,p,yv)
c     cette fonction récupère les coefs de la loi des gaz
c     raides connaissant yv
                  call eos(gam,pinf,yv)
                  cloc=dsqrt(gam*(p+pinf)/r)
                  vloc=cloc+dsqrt(u*u+v*v)
                  dtloc=vol(iv)/perim(iv)*cfl/vloc
                  if (iv.eq.1) then 
                     dt = dtloc
                     cmax=cloc
                     cmin=cloc
                     vmax=vloc
                  endif
                  if (dtloc.lt.dt) dt = dtloc
                  if (cloc.gt.cmax) cmax=cloc
                  if (vloc.gt.vmax) vmax=vloc
                  if (cloc.lt.cmin) cmin=cloc
               end do
               
               t=t+dt
               cpu2=etime(tarray)
               if (cpu2-cpu1.gt.raf) then
                  cpu1=etime(tarray)
                  cpu2=cpu1
                  write(*,1000) t,dt,vmax
 1000             format('t=',e12.5,' dt=',e12.5,' vmax=',e12.5)
               endif
            endif

c     termes source de pesanteur
            do iv=1,nbvol
               if (idur(iv).ne.-1) then
                  wnp1(3,iv)=wnp1(3,iv)-wn(1,iv)*g*dt/ic
                  wnp1(4,iv)=wnp1(4,iv)-wn(3,iv)*g*dt/ic
c     ordre deux sur la fraction de masse
      wnp1(5,iv)=wnp1(5,iv)-dt/ic*
     $              (wn(2,iv)/wn(1,iv)*dxyv(iv)
     $              +wn(3,iv)/wn(1,iv)*dyyv(iv))
               endif
            end do
            
c     boucle sur les aretes 

            do ia=1,nbaret
c     on récupère les n° des deux vf voisins 
c     "gauche" et "droit" : par convention la normale
c     à l'arête est orientée de "gauche" à "droite"
               ivg=iaret(1,ia)
               ivd=iaret(2,ia)
c     variables primitives du vf (g)
               call primitive(wn(1,ivg),rg,ug,vg,pg,yg)

c     ----------ordre 2--------------------------------------------
c     calcul des variables primitives
c     corrigées par les pentes
               if (iordre.eq.2) then
                  rg=rg+0.5d0*(xnorm(ia)*dx*dxr(ivg)
     $                 +ynorm(ia)*dy*dyr(ivg))
                  ug=ug+0.5d0*(xnorm(ia)*dx*dxu(ivg)
     $                 +ynorm(ia)*dy*dyu(ivg))
                  vg=vg+0.5d0*(xnorm(ia)*dx*dxv(ivg)
     $                 +ynorm(ia)*dy*dyv(ivg))
                  pg=pg+0.5d0*(xnorm(ia)*dx*dxp(ivg)
     $                 +ynorm(ia)*dy*dyp(ivg))
                  yg=yg+0.5d0*(xnorm(ia)*dx*dxyv(ivg)
     $                 +ynorm(ia)*dy*dyyv(ivg))
               endif
c     ----------fin ordre 2--------------------------------------------

c     on se place dans le repère lié à la normale à l'arête
c     calcul de la vitesse normale et tangentielle pour le vf (g)
               ung=ug*xnorm(ia)+vg*ynorm(ia)
               utg=-ug*ynorm(ia)+vg*xnorm(ia)
c     idem pour le vf (d) si l'arête n'est pas frontière
               if (ivd.ne.0) then
                  call primitive(wn(1,ivd),rd,ud,vd,pd,yd)

c     ----------ordre 2--------------------------------------------
c     calcul des variables primitives
c     corrigées par les pentes
                  if (iordre.eq.2) then
                     rd=rd-0.5d0*(xnorm(ia)*dx*dxr(ivd)
     $                    +ynorm(ia)*dy*dyr(ivd))
                     ud=ud-0.5d0*(xnorm(ia)*dx*dxu(ivd)
     $                    +ynorm(ia)*dy*dyu(ivd))
                     vd=vd-0.5d0*(xnorm(ia)*dx*dxv(ivd)
     $                    +ynorm(ia)*dy*dyv(ivd))
                     pd=pd-0.5d0*(xnorm(ia)*dx*dxp(ivd)
     $                    +ynorm(ia)*dy*dyp(ivd))
                     yd=yd-0.5d0*(xnorm(ia)*dx*dxyv(ivd)
     $                    +ynorm(ia)*dy*dyyv(ivd))
                  endif
c     ----------fin ordre 2--------------------------------------------

                  und=ud*xnorm(ia)+vd*ynorm(ia)
                  utd=-ud*ynorm(ia)+vd*xnorm(ia)
               else
c     si le vf (d) n'existe pas (arête frontière),
c     on recopie dans (d) l'état (g) ou bien on calcule
c     l'état miroir de (g)
                  rd=rg
                  pd=pg
                  yd=yg
c     condition imposée (test)
c     und=ung
c     condition miroir
                  und=ung*imiroir
                  utd=utg
               endif
c     résolution du problème de riemann dans la direction
c     normale à l'arête
               xi=0.d0
               call riemann(rg,ung,pg,yg,rd,und,pd,yd,xi,r,un,p,y)
C     calcul de unmoins et unplus pour appliquer la correction
c     non-conservative d'Abgrall-Saurel
               if (un.gt.0) then
                  unmoins=0
                  unplus=un
                  ut = utg
               else
                  unmoins=un
                  unplus=0
                  ut = utd
               endif

c     retour au repère (x,y)
               u=un*xnorm(ia)-ut*ynorm(ia)
               v=un*ynorm(ia)+ut*xnorm(ia)

               call eos(gam,pinf,y)
c     calcul du flux suivant x
               fx(1)=r*u
               fx(2)=r*u*u+p
               fx(3)=r*u*v
               fx(4)=((p+gam*pinf)/(gam-1.)+.5*r*(u*u+v*v)+p)*u

c     calcul du flux suivant y
               fy(1)=r*v
               fy(2)=r*u*v
               fy(3)=r*v*v+p
               fy(4)=((p+gam*pinf)/(gam-1.)+.5*r*(u*u+v*v)+p)*v

c     calcul du flux dans la direction normale
               do ii=1,4
                  flux(ii)=xnorm(ia)*fx(ii)+ynorm(ia)*fy(ii)
               end do

c     correction non-conservative pour le vf (d)
c     et incrément (s'il existe)
               flux(5)=-unplus*(yd-yg)

               if (ivd.ne.0.and.idur(ivd).ne.-1) then
                  do ii=1,5
                     wnp1(ii,ivd)=wnp1(ii,ivd)+dt/ic*along(ia)/vol(ivd)
     &                    *flux(ii)
                  enddo
               endif

c     correction non-conservative pour le vf (g)
c     et incrément
               flux(5)=unmoins*(yd-yg)

               if (idur(ivg).ne.-1) then
                  do ii=1,5
                     wnp1(ii,ivg)=wnp1(ii,ivg)-dt/ic*along(ia)/vol(ivg)
     &                    *flux(ii)
                  end do
               endif
            end do
 
c préparation du correcteur
            if (ic.eq.2) then
               do iv=1,nbvol
                  do ii=1,5
                     wn(ii,iv)=wnp1(ii,iv)
                     wnp1(ii,iv)=wdemi(ii,iv)
                  end do
               end do
            endif

         end do

      end do
      write(*,*) 'fin du calcul'
      write(*,*) 'itérations Riemann=',iterriem
      write(*,*) 'erreur Riemann=',errriem
      write(*,*) 'instant initial=',tini
      write(*,*) 'instant final  =',t
      cpu2=etime(tarray)
      write(*,*) 'temps CPU (secondes)',cpu2-cpu0

c tracé et sauvegarde
c voici des exemples de commandes
c gnuplot pour obtenir de jolis tracés !
c pour plus de détails consulter l'aide de gnuplot
c en tapant : help

c$$$ set view 0 , 0
c$$$ set contour
c$$$ set nosurface
c$$$ splot 'rho' w l
c$$$ pause -1 "densite"
c$$$ splot 'u' w l
c$$$ pause -1 "vitesse x"
c$$$ splot 'v' w l
c$$$ pause -1 "vitesse y"
c$$$ splot 'p' w l
c$$$ pause -1 "pression"
c$$$ splot 'yv' w l
c$$$ pause -1 "fraction liq"
c$$$ set cntrparam levels incremental 0,0.5,1
c$$$ splot 'yv' w l
c$$$ pause -1 "surface libre"

      open(2,file='rho')
      open(3,file='u')
      open(4,file='v')
      open(15,file='p')
      open(16,file='yv')
      open(7,file='vit')

c écriture du fichier de commande gnuplot
c s'utilise en faisant "gnuplot plotcom"
c nécessite une version de gnuplot > 3.8

      open(77,file='plotcom')

      write(77,*) 'set view 0,0,1,1'
      write(77,*) 'unset surface'
      write(77,*) 'set pm3d'
      write(77,*) 'splot ''rho'''
      write(77,*) 'pause -1'
      write(77,*) 'splot ''yv'''
      write(77,*) 'pause -1'
      close(77)

      write(*,*) 'écriture des fichiers de visu'

c écriture des variables primitives
      do j=0,m-1
         do i=0,n-1
            call primitive(wnp1(1,i+j*n+1),rho,u,v,xp,yv)            
            if (idur(i+j*n+1).eq.-1) then
               yv=2
               rho=1500
               u=0
               v=0
               xp=1e5
            endif
            write(2,*) rho
            write(3,*) u
            write(4,*) v
            write(15,*) xp
            write(16,*) yv
            write(7,*) dsqrt(u*u+v*v)
         enddo
         write(2,*)
         write(3,*)
         write(4,*)
         write(15,*)
         write(16,*)
         write(7,*)
      enddo

c sauvegarde en vrac des résultats

      write(*,*) 'sauvegarde binaire'

      open(98,form='unformatted',file='sauv')
      write(98) n,m
      write(98) t
      write(98) iordre
      write(98) irecif
      write(98) ((wnp1(ii,i),ii=1,5),i=1,nbvol)
      close(98)
            

      close(2)
      close(3)
      close(4)
      close(15)
      close(16)
      close(7)

      

      end


c calcul des 3 pentes à partir de d1, d2 et d3
c dans a,b,c et calcul du minmod de a,b et c
c renvoie le résultat dans x

      subroutine minmod(h,d1,d2,d3,x)

      implicit double precision(a-h,o-z)

      a=(d2-d1)/h
      b=(d3-d2)/h
      c=(d3-d1)/h*0.5d0

      isign=1
      if (a.lt.0) isign=-1
      x=dabs(a)

      if (dabs(b).lt.x) x=dabs(b)
      if (dabs(c).lt.x) x=dabs(c)
      if (a*b.lt.0.or.a*c.lt.0) isign=0

      x=x*isign
      

      return
      end

c ce sous-programme permet de calculer les variables
c conservatives w(1-5) connaissant les variables dites
c "primitives": densité, vitesse, pression, fraction de masse.
      subroutine conservative(w,r,u,v,p,yv)

      implicit double precision(a-h,o-z)
      dimension w(*)

      w(1)=r
      w(2)=r*u
      w(3)=r*v

      call eos(gam,pinf,yv)

      w(4)=.5*r*(u*u+v*v)+(p+gam*pinf)/(gam-1.)
c pour la correction saurel-abgrall, la 5ème variable
c n'est pas tout à fait conservative
      w(5)=yv

      return 
      end

c ce sous-programme permet de calculer les variables
c primitives  connaissant les variables conservatives
c c'est la fonction inverse du sous-programme précédent
      subroutine primitive(w,r,u,v,p,yv)

      implicit double precision(a-h,o-z)
      dimension w(*)  

      r=w(1)
      u=w(2)/r
      v=w(3)/r
      yv=w(5)

      call eos(gam,pinf,yv)

      p=(gam-1.d0)*(w(4)-.5d0*r*(u*u+v*v))-gam*pinf

      return 
      end


c ce sous-programme renvoie les coefficients de la loi
c d'état en fonction de la fraction yv.
      subroutine eos(gam,pinf,yv)

      implicit double precision(a-h,o-z)

      gam1=1.1d0
      gam2=1.1d0
      pinf1=263636.d0
      pinf2=-99636.d0

      t0=yv/(gam1-1.d0)+(1.d0-yv)/(gam2-1.d0)
      gam=1.d0+1.d0/t0
      t0=yv*gam1*pinf1/(gam1-1.d0)+(1.d0-yv)*gam2*pinf2/(gam2-1.d0)
      pinf=(gam-1.d0)*t0/gam

      return
      end


c la suite concerne la résolution du problème de Riemann
c ...

c234567
c$$$      program euler
c$$$      implicit double precision (a-h,o-z)
c$$$
c$$$      nmax=1000
c$$$      a=-6
c$$$      b=6
c$$$      h=(b-a)/nmax
c$$$
c$$$      rg=1000
c$$$      ug=0.
c$$$      pg=1.d5
c$$$      yg=1
c$$$
c$$$      rd=1
c$$$      ud=0
c$$$      pd=1.d5
c$$$      yd=0
c$$$
c$$$      time=1
c$$$
c$$$      open(1,file='rho')
c$$$      open(2,file='u')
c$$$      open(3,file='p')
c$$$      open(4,file='y')
c$$$      do i=0,nmax
c$$$         xi=(a+i*h)/time
c$$$         call riemann(rg,ug,pg,yg,rd,ud,pd,yd,xi,r,u,p,y)
c$$$         write(1,*) a+i*h,r
c$$$         write(2,*) a+i*h,u
c$$$         write(3,*) a+i*h,p
c$$$         write(4,*) a+i*h,y
c$$$      enddo
c$$$
c$$$      close(1)
c$$$      close(2)
c$$$      close(3)
c$$$      close(4)
c$$$
c$$$      return
c$$$      end

c courbe de Hugoniot pour les chocs
c234567
      function phia(pinfa,ga,ta,pa,pl)

      implicit double precision (a-h,o-z)

      t0=dsqrt(dabs((pl-pa)*(ta-ha(pinfa,ga,ta,pa,pl))))
      if (pl.le.pa) t0=-t0
  
      phia=t0
      
      return
      end

c dérivée par rapport à pl de la fonction précédente
c234567
      function dphia(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)
      pi=pl+pinfa
      pi0=pa+pinfa

      t0=.5d0*dsqrt(2.d0*ta/((ga-1.d0)*pi0+(ga+1.d0)*pi))
     $  -.5d0*dsqrt(((ga-1.d0)*pi0+(ga+1.d0)*pi)/2.d0/ta)
     &  *dha(pinfa,ga,ta,pa,pl) 

      if (pl.le.pa) t0=-t0
      dphia= t0

      return
      end

c calcul de tau=1/rho du côté (l) d'une 1 ou 3 onde en fonction de pl
      function ha(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)

      pi=pl+pinfa
      pi0=pa+pinfa
c cas du choc
      if (pl.gt.pa)  then
         t = ta*((ga+1.d0)*pi0+(ga-1.d0)*pi)/
     &    ((ga+1.d0)*pi+(ga-1.d0)*pi0)
c cas de la détente
      else 
         t = (pi0/pi)**(1.d0/ga)*ta
      endif
      ha=t

      return
      end

c dérivée par rapport à pl de la fonction précédente
      function dha(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)

      pi=pl+pinfa
      pi0=pa+pinfa
c choc
      if (pl.gt.pa) then
         t = -4.d0*ga*ta*pi0/(pi*(ga+1.d0)+pi0*(ga-1.d0))**2
c     détente (ne sert qu'au déboguage)
      else
         t = -ta*(pi0)**(1.d0/ga)/ga*(pi)**(-(ga+1.d0)/ga)
      endif
      dha=t

      return
      end

c fonction puissance
      function powa(a,b)
      implicit double precision (a-h,o-z)

      powa = a**b

      return
      end
     


c courbes isentropiques pour les détentes
      function psia(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)
  
      c0=dsqrt(ga*(pa+pinfa)*ta)
c      write(*,*) 'powa',pl,pa,pinfa,ga
      t1=2.d0*c0/(ga-1.d0)*(powa((pl+pinfa)/
     & (pa+pinfa),(ga-1.d0)/2.d0/ga)-1.d0)
      psia=t1

      return
      end


c dérivée par rapport à pl de la fonction précédente
      function dpsia(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)

      c0=dsqrt(ga*(pa+pinfa)*ta)
      t0=c0/ga*powa(pa+pinfa,(1.d0-ga)/2.d0/ga)*
     &  powa(pl+pinfa,-(ga+1.d0)/2.d0/ga)
      dpsia = t0

      return
      end

c calcul de la pression en fonction de la vitesse dans une détente
      function pp(pinfa,ga,ta,pa,Psi)
      implicit double precision (a-h,o-z)

      c0=dsqrt(ga*(pa+pinfa)*ta)
      c=(ga-1.d0)/(ga+1.d0)*(Psi+2.d0/(ga-1.d0)*c0)
      t0=c*c/ga/ta/powa(pa+pinfa,1.d0/ga)
      t0=powa(t0,ga/(ga-1.d0))-pinfa
      pp = t0

      return
      end
  
c courbes mixtes chocs-détentes
      function xhia(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)

      if (pl.gt.pa) then
         t0=phia(pinfa,ga,ta,pa,pl)
      else
         t0=psia(pinfa,ga,ta,pa,pl)
      endif
      xhia = t0
      return
      end

c dérivée par rapport à pl de la fonction précédente
      function dxhia(pinfa,ga,ta,pa,pl)
      implicit double precision (a-h,o-z)

      if (pl-pa.gt.0) then
c         write(*,*) 'choc 1',pl,pa,pl-pa
         t0=dphia(pinfa,ga,ta,pa,pl)
c         write(*,*) 'choc 2'
      else 
c         write(*,*) 'detente1'
         t0=dpsia(pinfa,ga,ta,pa,pl)
c         write(*,*) 'detente2'
      endif
      dxhia = t0

      return
      end


c fonction de résolution du problème de Riemann
      subroutine riemann(rg,ug,pg,psig,rd,ud,pd,psid,xi,r,u,p,psi)

      implicit double precision (a-h,o-z)
      dimension vit(5)
      
      common/iter/iterriem/err/errriem

c initialisation

      call eos(gamg,pinfg,psig)
      call eos(gamd,pinfd,psid)

      fg=gamg-1.d0
      fd=gamd-1.d0

c initialisation
      pn=0.d0
      err=1.d0
      eps=1.d-14
      it=0
      pmin=-pinfg
      if (pmin.lt.-pinfd) pmin=-pinfd
  
      gg=fg+1.d0
      gd=fd+1.d0

      
c critère d'apparition du vide
c      write(*,*) 'crit1',pd,eps*pd,pg,eps*pd
      crit=ug-ud-xhia(pinfd,fd+1.d0,1.d0/rd,pd,pmin+eps*pd)-
     $     xhia(pinfg,fg+1.d0,1.d0/rg,pg,pmin+eps*pg)

c      write(*,*) 'crit2'
      if (crit.le.0.d0) then
         write(*,*) 'apparition du vide'
         stop
      endif

c valeur initiale newton
c      pn=pg
c      if (pn.lt.pd) pn=pd
      pn=.5d0*(pg+pd) 
c nombre max d'itérations
      newton=100


      do while((err.gt.eps).and.(it.lt.newton)) 

         it=it+1
        
c     terme de gauche (voir Rouy)

         ff=ug-ud
         df=0.d0

         ff=ff-xhia(pinfg,gg, 1.d0/rg, pg, pn)
         df=df-dxhia(pinfg,gg, 1.d0/rg, pg, pn)

c     terme de droite (voir Rouy)
c         write(*,*) 'coucou1'
         ff=ff-xhia(pinfd,gd, 1.d0/rd, pd, pn)
c         write(*,*) 'coucou2',rd,pd,pn
         df=df-dxhia(pinfd,gd, 1.d0/rd, pd, pn)
c         write(*,*) 'coucou3'
         dp=ff/df

         pn=pn-dp

c         if (dabs(dp).lt.dabs(ff)) then
c            err = dabs(ff)
c         else
            err = dabs(dp/pn)
c         endif
         
      enddo
      
      if (it.gt.iterriem) iterriem=it
      if (err.gt.errriem) errriem=err

      if (it.eq.newton) then
         write(*,*) 'non convergence '
	 write(*,*) 'ff,dp,pn',ff,dp,pn
	 write(*,*) 'rg,ug,pg,rd,ud,pd,psid,psig',
     &    rg,ug,pg,rd,ud,pd,psid,psig
         stop
      endif
    
      pm=pn
      r2=1.d0/ha(pinfd,gd,1.d0/rd,pd,pn)
      r1=1.d0/ha(pinfg,gg,1.d0/rg,pg,pn)
    
c vitesse de la discontinuité de contact
      um1 = ug-xhia(pinfg,gg, 1.d0/rg, pg, pn)
      um2 = ud+xhia(pinfd,gd, 1.d0/rd, pd, pn)
      um=.5d0*(um1+um2)

c     vitesses caractéristiques
c     1-détente
      if (pm.le.pg) then
         vit(1)=ug-1.d0/dpsia(pinfg,gg,1.d0/rg,pg,pg)/rg
         vit(2)=um1-1.d0/dpsia(pinfg,gg,1.d0/rg,pg,pm)/r1
c     1-choc
      else
         vit(1)=ug-dsqrt( ((gg+1.d0)*(pm+pinfg)+
     & (gg-1.d0)*(pg+pinfg))*.5d0/rg)
         vit(2)=vit(1)
      endif

c     contact
      vit(3)=um

c     3-détente
      if (pm.le.pd) then
       vit(4)=um2+1.d0/dpsia(pinfd,gd,1.d0/rd,pd,pm)/r2
       vit(5)=ud+1.d0/dpsia(pinfd,gd,1.d0/rd,pd,pd)/rd
    
c     3-choc
      else      
         vit(4)=ud+dsqrt( ((gd+1.d0)*(pm+pinfd)+
     & (gd-1.d0)*(pd+pinfd))*.5d0/rd)
         vit(5)=vit(4)
      endif
      

c cette fonction calcule les variables primitives
c de la solution du problème de Riemann en xi=x/t

c état gauche
      if (xi.lt.vit(1)) then
         r=rg
         p=pg
         u=ug
         f=fg
         pinf=pinfg
         psi=psig
      endif


c     calcul dans la 1-détente
      if (xi.ge.vit(1).and.xi.lt.vit(2)) then
         p=pp(pinfg,fg+1.d0,1.d0/rg,pg,ug-xi)
         r=1.d0/ha(pinfg,fg+1.d0,1.d0/rg,pg,p)
         f=fg
         pinf=pinfg
         psi=psig
         u=ug-psia(pinfg,fg+1.d0,1.d0/rg,pg,p)
      endif
  
c     état (I) (voir Rouy)
      if (xi.ge.vit(2).and.xi.lt.vit(3)) then
         r=r1
         p=pm
         u=um1
         f=fg
         pinf=pinfg
         psi=psig
      endif

c     état (II) (voir Rouy)
      if (xi.ge.vit(3).and.xi.le.vit(4)) then
         r=r2
         p=pm
         u=um2
         f=fd
         pinf=pinfd
         psi=psid
      endif

c calcul dans la 3-détente
      if (xi.gt.vit(4).and.xi.lt.vit(5)) then
         p=pp(pinfd,fd+1.d0,1.d0/rd,pd,xi-ud)
         r=1.d0/ha(pinfd,fd+1.d0,1.d0/rd,pd,p)
         f=fd
         pinf=pinfd
         psi=psid
         u=ud+psia(pinfd,fd+1.d0,1.d0/rd,pd,p)
      endif
  
c     état droit
      if (xi.ge.vit(5)) then
         r=rd
         p=pd
         u=ud
         f=fd
         pinf=pinfd
         psi=psid
      endif

      return
      end
