C FDV32F    SOURCE    CB215821  16/04/21    21:16:48     8920
      SUBROUTINE fDV32F(INDMET,
     &                  alphaL, uvnL, uvtL, uvvL, ulnL, ultL, ulvL,
     &                  pL, TvL, TlL, rvL, rlL,
     &                  alphaR, uvnR, uvtR, uvvR, ulnR, ultR, ulvR,
     &                  pR, TvR, TlR, rvR, rlR,
     &                  F)
C***********************************************************************
C NOM         : fDV32F
C DESCRIPTION : Calculation of the Two phase AUSMDV flux
C               (3D two fluid model)
C
C LANGAGE     : ESOPE
C AUTEUR      : José R. Garcia Cascales (CEA/DEN/DM2S/SFME/LTMF)
C               mél : fd1@semt2.smts.cea.fr
C***********************************************************************
C APPELES          :
C APPELES (E/S)    :
C     inputs       : INDMET, variable to say if we want
C                    preconditioned version or not
C                    primitive variable at left side
C                    AL, UVNL, UVTL, UVVL, ULNL, ULTL, ULVL,
C                    PL, TVL, TLL, RVG, RLG,
C                    primitive variables at right side
C                    AR, UVNR, UVTR, UVVR, ULNR, ULTR, ULVR,
C                    PR, TVR, TLR, RVR, RLR
C
C     output       : F, 3D AUSMDV or "preconditioned" AUSMDV flux
C
C APPELES (BLAS)   :
C APPELES (CALCUL) :
C APPELE PAR       :
C***********************************************************************
C SYNTAXE GIBIANE    :
C ENTREES            :
C ENTREES/SORTIES    :
C SORTIES            :
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1, 13/03/2002, version initiale
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
*
      real*8 F(10)
      real*8 alphaL, alphaR, uvnL, uvnR, uvtL, uvtR,uvvL, uvvR,
     &       ulnL, ulnR, ultL, ultR, ulvL,ulvR,
     &       pL, pR, evL, evR, elL,
     &       elR, hvL, hvR, hlL, hlR, TvL, TvR, TlL, TlR
      integer INDMET
C variables for AUSMDV scheme
      real*8 rvL, rvR, rlL, rlR, avL, avR, alL, alR
      real*8 alm, avm, MlL, MlR, MvL, MvR, avpm, alpm,
     &       dotmv, dotml, Mvplus, Mvmin, Mlplus, Mlmin,
     &       Pvplus, Pvmin, Plplus, Plmin,
     &       fvL, fvR, flL, flR, wvplus, wvmin, wlplus, wlmin

      real*8 gamv, cpv, gaml, cpl, pil
      PARAMETER(gamv = 1.4D0, cpv = 1008.D0,
     &          gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8)

      evL = cpv*TvL/gamv
      evR = cpv*TvR/gamv
      elL = cpl*TlL/gaml + pil/rlL
      elR = cpl*TlR/gaml + pil/rlR

      hvL = evL + pL/rvL
      hvR = evR + pR/rvR
      hlL = elL + pL/rlL
      hlR = elR + pR/rlR

C Speed of sound at left and rigth sides
      avL = SQRT((gamv - 1.D0)*Cpv*TvL)
      avR = SQRT((gamv - 1.D0)*Cpv*TvR)
      avm = SQRT(avL*avR)
      alL = SQRT((gaml - 1.D0)*Cpl*TlL)
      alR = SQRT((gaml - 1.D0)*Cpl*TlR)
      alm = SQRT(alL*alR)

C Match number at both sides
      MlL = ulnL/alm
      MlR = ulnR/alm
      MvL = uvnL/avm
      MvR = uvnR/avm

C AUSMDV scheme
C      fvL = pL/(alphaL*rvL)
C      fvR = pR/(alphaR*rvR)
C      flL = pL/((1.D0 - alphaL)*rlL)
C      flR = pR/((1.D0 - alphaR)*rlR)

      fvL = 1.D0/(alphaL*rvL)
      fvR = 1.D0/(alphaR*rvR)
      flL = 1.D0/((1.D0 - alphaL)*rlL)
      flR = 1.D0/((1.D0 - alphaR)*rlR)

C     fvL = pL/(rvL)
C     fvR = pR/(rvR)
C     flL = pL/(rlL)
C     flR = pR/(rlR)

C     fvL = 1.D0/(rvL)
C     fvR = 1.D0/(rvR)
C     flL = 1.D0/(rlL)
C     flR = 1.D0/(rlR)

      wvplus = 2.D0*fvL/(fvL + fvR)
      wvmin = 2.D0*fvR/(fvL + fvR)
      wlplus = 2.D0*flL/(flL + flR)
      wlmin = 2.D0*flR/(flL + flR)

      call mpADV(wvplus, wvmin, MvL, MvR, Mvplus, Mvmin, Pvplus, Pvmin)
      call mpADV(wlplus, wlmin, MlL, MlR, Mlplus, Mlmin, Plplus, Plmin)

C AUSM schemes mass fluxes
      dotmv = avm*(alphaL*rvL*Mvplus + alphaR*rvR*Mvmin)
      dotml = alm*((1.D0 - alphaL)*rlL*Mlplus + (1.D0 - alphaR)*rlR
     $     *Mlmin)
      if(INDMET .EQ. 4)then
         dotml = dotml + 0.01d0*(pL - pR)
      end if
C AUSM schemes pressure fluxes
      avpm = Pvplus*alphaL*pL + Pvmin*alphaR*pR
      alpm = Plplus*(1.D0 - alphaL)*pL + Plmin*(1.D0 - alphaR)
     $     *pR

C Definition of convective part of the flux

      if (dotmv .GT. 0.D0) then

         F(1) = dotmv*1.D0
         F(3) = dotmv*uvnL
         F(5) = dotmv*uvtL
         F(7) = dotmv*uvvL
         F(9) = dotmv*(hvL + (uvnL**2 + uvtL**2 + uvvL**2)/2.D0)
      else
         F(1) = dotmv*1.D0
         F(3) = dotmv*uvnR
         F(5) = dotmv*uvtR
         F(7) = dotmv*uvvR
         F(9) = dotmv*(hvR + (uvnR**2 + uvtR**2 + uvvR**2)/2.D0)
      end if

      if (dotml .GT. 0.D0) then
         F(2) = dotml*1.D0
         F(4) = dotml*ulnL
         F(6) = dotml*ultL
         F(8) = dotml*ulvL
         F(10) = dotml*(hlL + (ulnL**2 + ultL**2 + ulvL**2)/2.D0)
      else
         F(2) = dotml*1.D0
         F(4) = dotml*ulnR
         F(6) = dotml*ultR
         F(8) = dotml*ulvR
         F(10) = dotml*(hlR + (ulnR**2 + ultR**2 + ulvR**2)/2.D0)
      end if

C Addition of the pressure term

      F(3) = F(3) + avpm
      F(4) = F(4) + alpm

      END