c balot.f: version 5 : avec optimisations
c compilation optimale avec intel fortran 8.0 :
c ifort -o balot balot.f -O3 -axN -ftz -ip -ipo -openmp

c compilation gnu : g77 -o balot balot.f  -O3 -ffast-math

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=0,xrecif=2.d0,hrecif=0.263d0,irecif=1,
     & ignuplot=1,itecplot=0,
     & nbaret=n*(m+1)+(n+1)*m,tmax=.1d0,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.and.ignuplot.eq.1) 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

               call calcpentes(n,m,wn,imiroir,dx,dy,
     $              dxr,dyr,
     $              dxu,dyu,dxv,dyv,
     $              dxp,dyp,dxyv,dyyv)
               
            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
c$$$               fx(1)=r*u
c$$$               fx(2)=r*u*u+p
c$$$               fx(3)=r*u*v
c$$$               fx(4)=((p+gam*pinf)/(gam-1.)+.5*r*(u*u+v*v)+p)*u
c$$$
c$$$c     calcul du flux suivant y
c$$$               fy(1)=r*v
c$$$               fy(2)=r*u*v
c$$$               fy(3)=r*v*v+p
c$$$               fy(4)=((p+gam*pinf)/(gam-1.)+.5*r*(u*u+v*v)+p)*v
c$$$
c$$$c     calcul du flux dans la direction normale
c$$$               do ii=1,4
c$$$                  flux(ii)=xnorm(ia)*fx(ii)+ynorm(ia)*fy(ii)
c$$$               end do
               
               flux(1)=r*un
               flux(2)=r*u*un+p*xnorm(ia)
               flux(3)=r*v*un+p*ynorm(ia)
               flux(4)=((p+gam*pinf)/(gam-1.)+.5*r*(u*u+v*v)+p)*un

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"

      if (ignuplot.eq.1) then

         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

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

      endif

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)
            
      if (itecplot.eq.1) then

         open(53,file='tecplot.dat')

         write(53,*)  ' TITLE= "t=',t,'"'
         write(53,*) 'VARIABLES ="X","Y","u","v","rho","yv","p" '
         write(53,1001) nbvol*4 , nbvol
 1001    format('ZONE N=',i8,' E=',i8,'F=FEPOINT , ET=quadrilateral') 

         if (iordre.eq.2) then

            call calcpentes(n,m,wnp1,imiroir,dx,dy,
     $           dxr,dyr,
     $           dxu,dyu,dxv,dyv,
     $           dxp,dyp,dxyv,dyyv)
         
         endif

         do iv=1,nbvol

            call primitive(wnp1(1,iv),r,u,v,p,yv)
            xmil=0
            ymil=0
            do ii=1,4
               xmil=xmil+xs(ivol(ii,iv))
               ymil=ymil+ys(ivol(ii,iv))
            enddo
            xmil=.25d0*xmil
            ymil=.25d0*ymil

       

            do ii=1,4
               xx=xs(ivol(ii,iv))
               yy=ys(ivol(ii,iv))
	       if(idur(iv).ne.-1)then
               write(53,'(7(x,e12.5))') xx,yy,
     $              u+(dxu(iv)*(xx-xmil)+dyu(iv)*(yy-ymil)),
     $              v+(dxv(iv)*(xx-xmil)+dyv(iv)*(yy-ymil)),
     $              r+(dxr(iv)*(xx-xmil)+dyr(iv)*(yy-ymil)),
     $              yv+(dxyv(iv)*(xx-xmil)+dyyv(iv)*(yy-ymil)),
     $              p+(dxp(iv)*(xx-xmil)+dyp(iv)*(yy-ymil))
              else
               write(53,'(7(x,e12.5))') xx,yy,
     $              0.,
     $              0.,
     $              1500.,
     $              2.,
     $              1.d5
              endif
c               write(*,*) u,dxu(iv),dyu(iv),iv,xmil,ymil,xx,yy
            enddo

         enddo

         ie=1
         do iv=1,nbvol
            write(53,'(4i8)') ie,ie+1,ie+2,ie+3
            ie=ie+4
         enddo

         close(53)


      endif
      

      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)
      parameter(
     & 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) 

      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 solveur VFROE
c ne marche pas

c$$$      ub=0.5d0*(ug+ud)
c$$$      rb=0.5d0*(rg+rd)
c$$$      pb=0.5d0*(pg+pd)
c$$$      psib=0.5d0*(psig+psid)
c$$$
c$$$      call eos(gamb,pinfb,psib)
c$$$
c$$$      cb=dsqrt(gamb*(pb+pinfb)/rb)
c$$$
c$$$      a1=.5d0/cb**2*(-rb*cb*(ud-ug)+pd-pg)
c$$$      a3=.5d0/cb**2*(rb*cb*(ud-ug)+pd-pg)
c$$$
c$$$      if (ub.gt.0) then
c$$$
c$$$         r=rg+a1
c$$$         u=ug-a1*cb/rb
c$$$         p=pg+a1*cb*cb
c$$$         psi=psig
c$$$
c$$$      else
c$$$
c$$$         r=rd-a3
c$$$         u=ud-a3*cb/rb
c$$$         p=pd-a3*cb*cb
c$$$         psi=psid
c$$$
c$$$      endif
c$$$
c$$$      return

c initialisation

      err=1.d0 
      eps=1.d-10
      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 débogage
c      write(*,*) 'crit1',pd,eps*pd,pg,eps*pd
c      crit=ug-ud-xhia(pinfd,fd+1.d0,1.d0/rd,pd,pmin+eps*pd)-
c     $     xhia(pinfg,fg+1.d0,1.d0/rg,pg,pmin+eps*pg)

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

c valeur initiale newton
      pn=0.5d0*(pg+pd)
c nombre max d'itérations
      newton=100
c nombre d'itérations pour le schéma rapide
      irap=1

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

c une seule itération de riemann
c ne marche pas !!!!!!!!!!!!!
c      do while((err.gt.eps).and.(it.lt.irap)) 

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

         ff=ug-ud
         df=0.d0


         xhig=xhia(pinfg,gg, 1.d0/rg, pg, pn)

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

c     terme de droite (voir Rouy)
         ff=ff-xhia(pinfd,gd, 1.d0/rd, pd, pn)
         df=df-dxhia(pinfd,gd, 1.d0/rd, pd, pn)

         dp=ff/df

         pn=pn-dp

c erreur relative
         err = dabs(dp/pn)

      enddo

c historique des convergences      
      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
    
c vitesse de la discontinuité de contact
      um = ug-xhig

c simplification faible mach
c à supprimer pour le solveur exact
c-----------------------------------
      u=um
      p=pm

      if (u.gt.0) then 
         r=rg
         psi=psig
      else
         if (u.lt.0) then
            r=rd
            psi=psid
         else
            r=0.5d0*(rg+rd)
            psi=0.5d0*(psig+psid)
         endif
      endif

      return
c----------------------------------
c fin simplification

      r2=1.d0/ha(pinfd,gd,1.d0/rd,pd,pn)
      r1=1.d0/ha(pinfg,gg,1.d0/rg,pg,pn)

      um1=um
      um2 = ud+xhia(pinfd,gd, 1.d0/rd, pd, pn)
      


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

      subroutine calcpentes(n,m,wn,imiroir,dx,dy,
     $     dxr,dyr,
     $     dxu,dyu,dxv,dyv,
     $     dxp,dyp,dxyv,dyyv)

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

      dimension wn(5,*),dxr(*),dyr(*),
     $     dxu(*),dyu(*),dxv(*),dyv(*),
     $     dxp(*),dyp(*),dxyv(*),dyyv(*)

      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
c               pn=p-9.81d0*rn*dy
               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
c               ps=p+9.81d0*rs*dy
               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
      
      return
      end