C FRUSB3    SOURCE    KLOCZKO   05/06/14    21:15:14     5111
C FRUSB3    SOURCE    CHAT      05/01/13    00:09:18     5004
      SUBROUTINE FRUSB3(NESP,
     &                  GAMG,ROG,PG,UNG,UTG,UVG,
     &                  GAMD,ROD,PD,UND,UTD,UVD,
     &                  YG,YD,VINF,FLU,
     &                  CELLT)
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FRUSB3
C
C DESCRIPTION       :  Formulation Volumes Finis  pour les Equations
C                      d'Euler Multi-Especes relatives à un melange
C                      de gaz parfaits.
C
C                      Calcul du flux aux interfaces avec la methode
C                      de Rusanov preconditionée
C
C LANGUAGE          :  FORTRAN 77
C
C AUTEUR            :  T. KLOCZKO DM2S/SFME/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     UVG, UVD        =  vitesses tangentielles
C
C     VINF            = vitesse de cut-off
C
C**** Sorties:
C
C     FLU             =  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*ut,
C                         rho*un*ht, rho*un*y1, ...)
C
C     CELLT           =  condition de stabilité, i.e.
C
C                        dT/diamax < cellt
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créé le 02.05.05
C
C
C
C************************************************************************
C
C N.B.: Toutes les variables sont DECLAREES
C
C
      IMPLICIT INTEGER(I-N)
      INTEGER NESP,I1,I2
      REAL*8  GAMG,ROG,PG,UNG,UTG,UVG,HG,RETG,CG
     &       ,GAMD,ROD,PD,UND,UTD,UVD,HD,RETD,CD
     &       ,YG(*),YD(*),VINF
     &       ,FLU(*),CELLT, FLUM
     &       ,UNF,UTF,UVF,ROF,PF,GAMF,HF,VPMAX,PHIM
     &       ,AMAT(5,5),XP(5),FC(5),DS(5)
     &       ,GAM1,XU,CF,Q2,COEF
C
C**** Etat moyen
C
      UNF  = 0.5D0 * (UNG  + UND)
      UTF  = 0.5D0 * (UTG  + UTD)
      UVF  = 0.5D0 * (UVG  + UVD)
      ROF  = 0.5D0 * (ROG  + ROD)
      PF   = 0.5D0 * (PG   + PD)
      GAMF = 0.5D0 * (GAMG + GAMD)
      GAM1 = GAMF - 1.D0
      CF    = (GAMF * PF / ROF)**0.5D0
C
C**** We include in the cut-off UNG,UND,UTG,UTD in order to
C     avoid low diffusivity in stagnation regions
C
      VINF = MAX(VINF,((UNG**2 + UTG**2 + UVG**2)**0.5D0))
      VINF = MAX(VINF,((UND**2 + UTD**2 + UVD**2)**0.5D0))
C
C*** Calcul de la condition de stabilité
C
      CELLT = 1.0D0 / (CF + ABS(UNF))
C
C**** DELTA
C
      RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG)
     &     + (0.5D0 * ROG * (UNG**2 + UTG**2 + UVG**2))
      RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD)
     &     + (0.5D0 * ROD * (UND**2 + UTD**2 + UVD**2))
C
      HG   = (1.D0 / ROG) * (RETG + PG)
      HD   = (1.D0 / ROD) * (RETD + PD)
      HF   = 0.5D0 * (HG + HD)
C
      XP(1) = ROG - ROD
      XP(2) = (ROG * UNG) - (ROD * UND)
      XP(3) = (ROG * UTG) - (ROD * UTD)
      XP(4) = (ROG * UVG) - (ROD * UVD)
      XP(5) = RETG - RETD
C
C*** Calcul du nombre de Mach de référence PHIM
C
      CG = SQRT((GAMG * PG) / ROG)
      CD = SQRT((GAMD * PD) / ROD)
      IF(VINF.GE.MAX(CD,CG))THEN
      PHIM = 1.D0
      ELSE
      PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD))
      ENDIF
C
C*** Calcul des valeurs propres
C
      XU    = ((1.D0 - PHIM**2) * UNF)**2 + 4.D0 * (PHIM * CF)**2
C
      VPMAX = 0.5D0 * ((1.D0 + PHIM**2) * ABS(UNF) + SQRT(XU))
C
C
C*** Calcul de la matrice de préconditionnement
C
      Q2    = 0.5D0 * (UNF**2 + UTF**2 + UVF**2)
      COEF  = (1.0D0 / PHIM**2 - 1.0D0) * GAM1 / CF**2
C
      AMAT(1,1) = COEF * 1.D0 * Q2    + 1.D0
      AMAT(1,2) = COEF * 1.D0 * (-UNF)
      AMAT(1,3) = COEF * 1.D0 * (-UTF)
      AMAT(1,4) = COEF * 1.D0 * (-UVF)
      AMAT(1,5) = COEF * 1.D0 * 1.D0
C
      AMAT(2,1) = COEF * UNF  * Q2
      AMAT(2,2) = COEF * UNF  * (-UNF) + 1.D0
      AMAT(2,3) = COEF * UNF  * (-UTF)
      AMAT(2,4) = COEF * UNF  * (-UVF)
      AMAT(2,5) = COEF * UNF  * 1.D0
C
      AMAT(3,1) = COEF * UTF  * Q2
      AMAT(3,2) = COEF * UTF  * (-UNF)
      AMAT(3,3) = COEF * UTF  * (-UTF) + 1.D0
      AMAT(3,4) = COEF * UTF  * (-UVF)
      AMAT(3,5) = COEF * UTF  * 1.D0
C
      AMAT(4,1) = COEF * UVF  * Q2
      AMAT(4,2) = COEF * UVF  * (-UNF)
      AMAT(4,3) = COEF * UVF  * (-UTF)
      AMAT(4,4) = COEF * UVF  * (-UVF) + 1.D0
      AMAT(4,5) = COEF * UVF  * 1.D0
C
      AMAT(5,1) = COEF * HF   * Q2
      AMAT(5,2) = COEF * HF   * (-UNF)
      AMAT(5,3) = COEF * HF   * (-UTF)
      AMAT(5,4) = COEF * HF   * (-UVF)
      AMAT(5,5) = COEF * HF   * 1.D0  + 1.D0
C
C*** Calcul de la  dissipation
C
      DO I1 = 1,5
         DS(I1) = 0.D0
      DO I2 = 1,5
          DS(I1) =  DS(I1) + 0.5D0 * VPMAX * AMAT(I1,I2) * XP(I2)
      ENDDO
      ENDDO
C
C*** Calcul du flux convectif
C
      FC(1) = 0.5D0 * ((ROG * UNG)            + (ROD * UND))
      FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD))
      FC(3) = 0.5D0 * ((ROG * UNG * UTG)      + (ROD * UND * UTD))
      FC(4) = 0.5D0 * ((ROG * UNG * UVG)      + (ROD * UND * UVD))
      FC(5) = 0.5D0 * ((UNG * (RETG + PG))    + (UND * (RETD + PD)))
C
C*** Calcul du flux numérique
C
      DO I1 = 1,5
         FLU(I1) = FC(I1) + DS(I1)
      ENDDO
C
C*** Partie multi-espèces
C
      FLUM = FLU(1)
      IF(FLUM .GT. 0)THEN
         DO I1 = 1, NESP, 1
            FLU(5+I1)=FLUM * YG(I1)
         ENDDO
      ELSE
         DO I1 = 1, NESP, 1
            FLU(5+I1)=FLUM * YD(I1)
         ENDDO
      ENDIF
C
      RETURN
      END