C FLHLLC    SOURCE    BECC      10/05/05    21:15:10     6674
      SUBROUTINE FLHLLC(NESP,
     &     GAMG,ROG,PG,UNG,UTG,
     &     GAMD,ROD,PD,UND,UTD,
     &     YG,YD,FLU1,
     &     CELLT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FLHLLC
C
C DESCRIPTION       :  Formulation Volumes Finis  pour les Equations
C                      d'Euler Multi-Especes relatives à un melange
C                      de gaz ideals.
C
C                      Calcul du flux aux interfaces avec la methode
C                      de HLLC (avec gamma = const)
C                      Voir: * Batten et al. SIAM J. Sci. Comp.
C                              Vol. 18 pp 1553-6, 1997
C                            * Batten et al. JCP, 1997
C                              Vol. 137 pp 38--
C
C                      N.B.: SPLIT DE FRACTIONS MASSIQUES ET UNT A LA
C                            LARROUTUROU
C
C LANGUAGE          :  FORTRAN 77
C
C AUTEUR            :  A. BECCANTINI DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C APPELES
C
C************************************************************************
C
C**** Entrées:
C
C     NESP            =  nombre d'especes considérées dans les Equations
C                        d'Euler
C
C     GAMG, GAMD      =  les "gamma" du gaz (gauche et droite)
C
C     ROG, ROD        =  les densités
C
C     PG, PD          =  les pressions
C
C     UNG, UND        =  vitesses normales
C
C     UTG, UTD        =  vitesses tangentielles
C
C     YG, YD          =  tables des fractiones massiques
C
C**** Sorties:
C
C     FLU1            =  table du flux a l'interface dans le repaire
C                        (n,t), i.e.
C                        (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
C                         rho*un*y1, ...)
C
C     CELLT           =  condition de stabilité, i.e.
C
C                        dT/diamax < cellt (temps/longueur)
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créé le 11.09.2000
C
C************************************************************************
C
C N.B.: Toutes les variables sont DECLAREES
C
C     It exactly captures:
C
C     Stationary shocks (test for instance
C
C       GAMG=1.4D0
C       GAMD=1.4D0
C       ROG=1.0D0
C       ROD=4.5584471D0
C       PG=1.0D5
C       PD=0.18279373D7
C       UNG=0.14877919D4
C       UND=0.32638130D3
C       UTG=0.0D0
C       UTD=0.0D0
C     ) and stationary contacts
C
      IMPLICIT INTEGER(I-N)
      INTEGER NESP
     &     , IESP
      REAL*8  GAMG,ROG,PG,UNG,UTG
     &     ,GAMD,ROD,PD,UND,UTD
     &     ,YG(*),YD(*),FLU1(*),CELLT
     &     ,HALF
      PARAMETER (HALF=0.5D0)
      REAL*8 USGGM1,USGDM1,HTG,HTD,REING,REIND,SQRG,SQRD,USOSQR
     &     ,HTAVE,UNAVE,REIAVE,USGAVE,GAMAVE,ROAVE,PAVE,ASONAV
     &     ,ASONG,ASOND,SG,SD,SM,PSTAR,RETG,RETD,OMEGA
     &     ,ROGS,RUNGS,RETGS,RODS,RUNDS,RETDS
C
      USGGM1=1.0D0/(GAMG-1.0D0)
      USGDM1=1.0D0/(GAMD-1.0D0)
C
C     Evaluation of the Roe average state to compute the non-linear wave
C     speeds (According to Shuye, JCP 142, pag. 217--, 1998
C
      REING=PG*USGGM1
      REIND=PD*USGDM1
      HTG=(HALF*UNG*UNG)+(GAMG*REING/ROG)
      HTD=(HALF*UND*UND)+(GAMD*REIND/ROD)
      SQRG=ROG**0.5D0
      SQRD=ROD**0.5D0
      USOSQR=1.0D0/(SQRG+SQRD)
      HTAVE=((SQRG*HTG)+(SQRD*HTD))*USOSQR
      UNAVE=((SQRG*UNG)+(SQRD*UND))*USOSQR
      USGAVE=((SQRG*USGGM1)+(SQRD*USGDM1))*USOSQR
      GAMAVE=(1.0D0/USGAVE)+1.0D0
      REIAVE=((SQRG*REING)+(SQRD*REIND))*USOSQR
C
      ROAVE=REIAVE*GAMAVE/(HTAVE-(HALF*UNAVE*UNAVE))
      PAVE=REIAVE*(GAMAVE-1.0D0)
C
      ASONAV=(GAMAVE*PAVE/ROAVE)**0.5D0
C
C**** Evaluation of the speed of the Genuinely-Non-Linear Waves
C
      ASONG=(GAMG*PG/ROG)**0.5D0
      ASOND=(GAMD*PD/ROD)**0.5D0
C
      SG=UNG-ASONG
      SG=MIN(SG,(UNAVE-ASONAV))
      SD=UND+ASOND
      SD=MAX(SD,(UNAVE+ASONAV))
C
C**** Speed of the contact discontinuity
C
      SM=ROD*UND*(SD-UND)-ROG*UNG*(SG-UNG)+PG-PD
      SM=SM/((ROD*(SD-UND))-(ROG*(SG-UNG)))
C
C**** Interfacial pressure
C
      PSTAR=ROG*(UNG-SG)*(UNG-SM)+PG
      RETG=(USGGM1*PG)+(HALF*ROG*((UNG*UNG)+(UTG*UTG)))
      RETD=(USGDM1*PD)+(HALF*ROD*((UND*UND)+(UTD*UTD)))
C
      IF(SG .GT. 0)THEN
         FLU1(1)=ROG*UNG
         FLU1(2)=ROG*UNG*UNG+PG
         FLU1(4)=UNG*(RETG+PG)
      ELSEIF(SM .GT. 0)THEN
         OMEGA=1.0D0/(SG-SM)
         ROGS=ROG*(SG-UNG)
         ROGS=ROGS*OMEGA
         RUNGS=((SG-UNG)*(ROG*UNG))+(PSTAR - PG)
         RUNGS=RUNGS*OMEGA
         RETGS=((SG-UNG)*RETG)+(PSTAR*SM)-(PG*UNG)
         RETGS=RETGS*OMEGA
C
         FLU1(1)=ROGS*SM
         FLU1(2)=RUNGS*SM+PSTAR
         FLU1(4)=(RETGS+PSTAR)*SM
C
      ELSEIF(SD .GT. 0)THEN
         OMEGA=1.0D0/(SD-SM)
         RODS=ROD*(SD-UND)
         RODS=RODS*OMEGA
         RUNDS=((SD-UND)*(ROD*UND))+(PSTAR - PD)
         RUNDS=RUNDS*OMEGA
         RETDS=((SD-UND)*RETD)+(PSTAR*SM)-(PD*UND)
         RETDS=RETDS*OMEGA
C
         FLU1(1)=RODS*SM
         FLU1(2)=RUNDS*SM+PSTAR
         FLU1(4)=(RETDS+PSTAR)*SM
      ELSE
         FLU1(1)=ROD*UND
         FLU1(2)=ROD*UND*UND+PD
         FLU1(4)=UND*(RETD+PD)
      ENDIF
C
      FLU1(3)=(FLU1(1)*(UTG+UTD))+(ABS(FLU1(1))*(UTG-UTD))
      FLU1(3)=HALF*FLU1(3)
C
      DO IESP=1,NESP,1
         FLU1(4+IESP)=(FLU1(1)*(YG(IESP)+YD(IESP)))+
     &        (ABS(FLU1(1))*(YG(IESP)-YD(IESP)))
         FLU1(4+IESP)=HALF*FLU1(4+IESP)
      ENDDO
C
      CELLT=1.0D0/(MAX(SG,SD))
C
      RETURN
      END