! mettre ici les variables globales de votre programme
! par exemple: b, gamma, etc.
module global
  implicit none
  integer,parameter               :: nvar=2,nmax=300   ! nombre de variables, nb de cellules
  real(kind=kind(0.d0)),parameter :: g=9.81d0   ! acceleration pesanteur (pour st venant)
end module global


! ce programme resout les equations de st venant par le schema de vfroe
! les equations de st venant ou de l'eau peu profonde sont
! decrites dans le polycopie systemes hyperboliques
! sur mon site web. elles permettent de modélider une rupture de barrage par
! exemple. Les inconnues sont la hauteur d'eau h et la vitesse horizontale de la 
! colonne d'eau u.
! vous devez modifier le programme pour prendre en compte les 
! equations de la mhd ideale (7 equations au lieu de 2)
program vfroe_stvenant

  use global
  implicit none

  ! variables conservatives en tn et tn+1
  real(kind=kind(0.d0)),dimension(nvar,0:nmax+1) :: wn,wnp1
  ! milieux des cellules
  real(kind=kind(0.d0)),dimension(0:nmax+1) :: xmil(0:nmax+1)
  ! tableaux utiles (flux, etats primitifs)
  real(kind=kind(0.d0)),dimension(nvar) :: flux(nvar),yL,yR

  ! declarations
  real(kind=kind(0.d0)) :: a1,a2,hLini,hRini,uLini,uRini,tmax,dx,cfl,t,dt
  real(kind=kind(0.d0)) :: vmax,h,u,c,vloc
  integer ::    i,ii

  ! initialisations
  ! bornes en espace
  a1=-2.d0
  a2=2.d0

  ! données initiales d'un pb de riemann
  !hauteur d'eau et vitesse a gauche
  hLini=2.d0
  uLini=0
  ! hauteur d'eau et vitesse a droite
  hRini=1.d0
  uRini=0

  ! initialisations des variables primitives
  ! (a modifier: initialiser les 7 variables rho ,u ,p ,v2, v3 B2 B3)
  yL(1)=hLini
  yL(2)=uLini
  yR(1)=hRini
  yR(2)=uRini

  ! borne en temps
  !write(*,*) 'entrer le temps final'
  !read(*,*) tmax
  tmax=0.1d0
  ! pas d'espace
  dx=(a2-a1)/nmax
  ! coefficient de CFL
  cfl=0.8d0

  ! vecteur initiaux conservatifs dans les cellules
  ! on initialise aussi les deux cellules fictives
  ! i=0 et i=nmax+1

  ! conditions initiales 

  do i=0,nmax+1
     xmil(i)=a1+i*dx-0.5d0*dx
     if (xmil(i).lt.0) then
        call conservatives(yL,wn(1,i))
     else
        call conservatives(yR,wn(1,i))
     endif
     do ii=1,nvar
        wnp1(ii,i)=wn(ii,i)
     end do
  enddo

  ! debut de la boucle en temps
  t=0.d0
  do while(t.lt.tmax)
     ! recherche de la vitesse max
     ! pour calculer la cfl
     vmax=0.d0
     do i=0,nmax+1
        call cvmax(wn(1,i),vloc)
        if (vloc.gt.vmax) vmax=vloc
     enddo
     dt=cfl*dx/2.d0/vmax
     ! boucle sur les cellules internes i de 1 à nmax
     do i=1,nmax
        ! flux à droite
        call fluxnum(wn(1,i),wn(1,i+1),flux)
        ! ajout dans la cellule
        do ii=1,nvar
           wnp1(ii,i)=wnp1(ii,i)-dt/dx*flux(ii)
        end do
        ! flux à gauche
        call fluxnum(wn(1,i-1),wn(1,i),flux)
        do ii=1,nvar
           wnp1(ii,i)=wnp1(ii,i)+dt/dx*flux(ii)
        end do
     enddo
     ! actualisations
     do i=1,nmax
        do ii=1,nvar
           wn(ii,i)=wnp1(ii,i)
        end do
     enddo
     ! décommentez les lignes suivantes si vous voulez changer
     ! de condition limite
     !wn(1,0)=wnp1(1,1)
     !wn(2,0)=-wnp1(2,1)
     !wn(1,nmax+1)=wnp1(1,nmax)
     !wn(2,nmax+1)=-wnp1(2,nmax)

     t=t+dt
     write(*,*) 't=',t
  enddo

  ! écriture de la solution exacte et de la
  ! solution numérique dans 4 fichiers
  ! pour tracé dans "gnuplot"
  ! taper "gnuplot" puis, dans gnuplot
  ! "  plot 'hnum', 'hriem'  " par exemple
  ! a modifier pour la mhd bien sur...

  open(1,file='hnum')
  open(2,file='hriem')
  open(3,file='unum')
  open(4,file='uriem')

  !      do i=1,nmax
  !         write(1,*) xmil(i),wn(i)
  !         call solexacte(xmil(i),t,w)
  !         write(2,*) xmil(i),w
  !      enddo

  do i=1,nmax 
     write(1,*) xmil(i) , wn(1,i)
     write(3,*) xmil(i), wn(2,i)/wn(1,i)
  enddo

! tracé de la solution exacte dans le cas de st venant
! (pour la mhd, il faudra utiliser maple pour fabriquer les fichiers
! de solutions exactes)
  do i=1,nmax
     call riemann(hLini,uLini,hRini,uRini,xmil(i)/t,h,u)
     write(2,*) xmil(i) , h
     write(4,*) xmil(i) , u
  enddo


  close(1)
  close(2)
  close(3)
  close(4)

end program vfroe_stvenant



! calcul du flux numerique
! en utilisant un solveur de vfroe
subroutine fluxnum(wL,wR,flux)

  use global
  implicit none
  real(kind=kind(0.d0)),dimension(nvar),intent(in) :: wL,wR
  real(kind=kind(0.d0)),dimension(nvar),intent(out) :: flux
  real(kind=kind(0.d0)) :: xi
  real(kind=kind(0.d0)),dimension(nvar) :: yL,yR,y,w

  call primitives(wL,yL)
  call primitives(wR,yR)

  ! les deux lignes suivantes correspondent
  ! au schema de godunov
  !call riemann(yL(1),yL(2),yR(1),yR(2),0.d0,y(1),y(2))
  !call conservatives(y,w)
  xi=0.d0
  !resolution du pb de riemann linearise
  call riemann_vfroe(wL,wR,xi,w)
  ! flux de l'etat intermediaire
  call cflux(w,flux)

end subroutine fluxnum

! passage des variables primitives aux variables conservatives
! (a modifier)
subroutine conservatives(y,w)

  use global
  implicit none
  real(kind=kind(0.d0)),dimension(nvar),intent(in) :: y
  real(kind=kind(0.d0)),dimension(nvar),intent(out) :: w

  w(1)=y(1)          ! h
  w(2)=y(1)*y(2)     ! h*u

end subroutine conservatives

! passage des variables conservatives aux
! variables primitives
subroutine primitives(w,y)

  use global
  implicit none
  real(kind=kind(0.d0)),dimension(nvar),intent(in) :: w
  real(kind=kind(0.d0)),dimension(nvar),intent(out) :: y

  y(1)=w(1)           ! h    : hauteur d'eau
  y(2)=w(2)/w(1)      ! u    : vitesse horizontale de la colonne d'eau

end subroutine primitives

! calcul du flux connaissant un etat conservatif w
subroutine cflux(w,flux)

  use global
  implicit none
  real(kind=kind(0.d0)),dimension(nvar),intent(in) :: w
  real(kind=kind(0.d0)),dimension(nvar),intent(out) :: flux

  real(kind=kind(0.d0)),dimension(nvar)  :: y

  call primitives(w,y)

  flux(1)=w(2)                          ! h*u
  flux(2)=y(1)*y(2)**2+g/2*y(1)**2      ! h*u**2 + 1/2*g*h*u**2
                                        ! g = acceleration peasanteur



end subroutine cflux

! calcul de la plus grande vp en valeur absolue de la matrice jacobienne
! ici les vp sont u +/- sqrt(g*h)
! à modifier pour la mhd...
subroutine cvmax(w,vmax)

  use global
  implicit none
  real(kind=kind(0.d0)),dimension(nvar),intent(in) :: w
  real(kind=kind(0.d0)),intent(out) :: vmax

  real(kind=kind(0.d0)) :: h,u,c

  h=w(1)
  u=w(2)/h
  c=dsqrt(g*h)
  
  vmax=dabs(u)+c

end subroutine cvmax


! solveur de riemann linearise pour st venant
subroutine riemann_vfroe(wL,wR,xi,wstar)
  !=========================== DOCUMENTATION ============================= 
  !    NOM
  ! riemann_vfroe(wL,wR,vn,xi,w)
  !    DESCRIPTION  
  ! fonction de résolution du problème de Riemann linéarisé pour le schéma de vfroe
  ! 
  !    ENTREES 
  !  wL  : état gauche
  !  wR  : état droit
  !  xi  : x/t 
  !
  !    SORTIES 
  !  wstar : état solution en x/t
  !
  !    SOURCE
  !=========================== DEBUT DES DECLARATIONS ==================== 
  !.1-----  Implicit, Use
  use global
  implicit none
  !.2-----  Declaration
  real(kind=kind(0.d0)),dimension(nvar), intent(in) :: wL , wR 
  real(kind=kind(0.d0)),intent(in):: xi
  real(kind=kind(0.d0)),dimension(nvar),intent(out):: wstar

  real(kind=kind(0.d0)),dimension(nvar) :: yL,yR,ystar,ymil

  real(kind=kind(0.d0)),dimension(:),allocatable   :: du,adu,aadu
  real(kind=kind(0.d0)),dimension(:,:),allocatable :: A
  real(kind=kind(0.d0)) :: h,u,c,a0,a1,a2,signum
  integer :: i
  !=========================== DEBUT DU CODE EXECUTABLE ==================
  !
  !-------  Initialisations
  !
  ! matrice jacobienne
  allocate(A(nvar,nvar))
  A=0.d0
  ! vecteurs de travail
  allocate(du(nvar))
  allocate(adu(nvar))
  allocate(aadu(nvar))

  ! passage en variables primitives
  call primitives(wL,yL)
  call primitives(wR,yR)

  ! calcul de l'etat moyen autour duquel
  ! on linearise les equations
  ymil=0.5d0*(yL+yR)

  h=ymil(1)
  u=ymil(2)
  c=sqrt(g*h)

  ! jacobienne de l'état moyen
  A(1,1) = u
  A(1,2) = h
  A(2,1) = g
  A(2,2) = u

  ! du = yR-yL
  do i=1,nvar      
     du(i) =yR(i)-yL(i)
  end do
  ! vecteur A*du
  adu =matmul(A,du)
  ! vecteur A**2 du
  aadu=matmul(A,adu)

  !coefficients des polynômes de décentrage (de degré 2 ici car j'ai interpolé aussi 
  ! en x=u bien que u ne soit pas valeur propre...)
  ! le calcul n'est valable que pour x/t=0
  a2 = (signum(u+c)-2.D0*signum(u)+signum(u-c))/c**2/2.D0
  a1 = (-2.D0*signum(u+c)*u+signum(u+c)*c-2.D0*signum(u-c)*u+ & 
       4.D0*signum(u)*u-c*signum(u-c))/c**2/2.D0
  a0 = (2.D0*c**2*signum(u)+signum(u+c)*u**2-signum(u+c)*u*c- &
       2.D0*signum(u)*u**2+signum(u-c)*u*c+signum(u-c)*u**2)/c**2/2.D0

  ! pour la mhd il faudra optimiser le calcul par l'algorithme de
  ! horner car le polynôme de matrice est de degré 7 !!!

  ! calcul de l'état décentré
  do i=1,nvar      
     ystar(i)=ymil(i)-0.5d0*(a2*aadu(i)+a1*adu(i)+a0*du(i))
  end do

  ! retour aux variables conservatives
  call conservatives(ystar,wstar)

  !===========================   FIN DE LA ROUTINE    ====================
end subroutine riemann_vfroe
!DOCFIN 

 function signum(x)

  implicit none
  real(kind=kind(0.d0))::x, signum

  if (x.gt.0.d0) then 
     signum=1.d0
  else if (x.lt.0.d0) then
     signum=-1.d0
  else
     signum=0.0d0
  end if


end function signum




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! fonction Z(h1,h2)

function z(h1,h2)
  implicit double precision (a-h,o-z)

  g=9.81

  if (h1.gt.h2) then
     zsol=sqrt(2.d0)*sqrt(g*(h1+h2)/h1/h2)/2.d0
  endif

  if (h1.le.h2) then
     zsol=2.d0*dsqrt(g)/(sqrt(h1)+sqrt(h2))
  endif

  z=zsol

end function z
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! derivee de z par rapport a h1

function dzdh1(h1,h2)
  implicit double precision (a-h,o-z)

  g=9.81

  if (h1.gt.h2) then
     dzsol= sqrt(2.d0)/sqrt(g*(h1+h2)/h1/h2)*   &
          (g/h1/h2-g*(h1+h2)/h1**2/h2)/4.d0
  endif

  if (h1.le.h2) then
     dzsol=-2.d0*sqrt(g)/(sqrt(h1)+sqrt(h2))**2/sqrt(h1)/2
  endif

  dzdh1=dzsol


end function dzdh1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! resolution du pb de riemann

subroutine riemann(hg,ug,hd,ud,xi,h,u)
  implicit double precision (a-h,o-z)

  ! accel pesanteur
  g=9.81d0

  ! valeur initiale de hstar pour newton
  hstar=.5*(hg+hd)

  ! debut des iterations
1 continue

  ! calcul du residu (on veut que res=0)
  res= ud-ug+(hstar-hd)*z(hstar,hd)  &
       + (hstar-hg)*z(hstar,hg)

  ! calcul de l'accroissement
  dh=z(hstar,hd) +(hstar-hd)*dzdh1(hstar,hd)  &
       + z(hstar,hg) +(hstar-hg)*dzdh1(hstar,hg) 

  dh=-res/dh

  ! mise a jour
  hstar=hstar+dh

  ! test d'arret
  if ((dabs(dh).gt.1.d-10).and.(dabs(res).gt.1.d-10)) then 
     goto 1
  endif

  !      write(*,*) hstar

  ustar=ud+(hstar-hd)*z(hstar,hd)
  !      write(*,*) ustar,ug-(hstar-hg)*z(hstar,hg),res

  ! celerite du "son" a gauche droite et milieu
  cg=dsqrt(g*hg)
  cd=dsqrt(g*hd)
  cstar=dsqrt(g*hstar)

  ! calcul de lambda(i,plus) et lambda(i,moins)
  ! onde 1
  if (hstar.gt.hg) then
     ! choc
     xlambda1m=(hg*ug-hstar*ustar)/(hg-hstar)
     xlambda1p=xlambda1m
  else
     ! detente
     xlambda1m=ug-cg
     xlambda1p=ustar-cstar
  endif

  ! onde 2
  if (hstar.gt.hd) then
     ! choc
     xlambda2m=(hd*ud-hstar*ustar)/(hd-hstar)
     xlambda2p=xlambda2m
  else
     ! detente
     xlambda2p=ud+cd
     xlambda2m=ustar+cstar
  endif

  !      write(*,*) xlambda1m,xlambda1p,xlambda2m,xlambda2p

  ! en fonction de xi on evalue h et u

  if (xi.lt.xlambda1m) then
     h=hg
     u=ug
  endif

  if ((xi.ge.xlambda1m).and.(xi.lt.xlambda1p)) then
     ! 1-detente
     h=(2.d0*dsqrt(hg)/3.d0+(ug-xi)/3.d0/dsqrt(g))**2;
     u=ug-(h-hg)*z(h,hg)
  endif

  if ((xi.ge.xlambda1p).and.(xi.lt.xlambda2m)) then
     ! etat milieu
     h=hstar
     u=ustar
  endif

  if ((xi.ge.xlambda2m).and.(xi.lt.xlambda2p)) then
     ! 2-detente
     h=(2.d0*dsqrt(hd)/3.d0+(xi-ud)/3.d0/dsqrt(g))**2;
     u=ud+(h-hd)*z(h,hd)
  endif

  if (xi.ge.xlambda2p) then
     h=hd
     u=ud
  endif

  return

end subroutine riemann


