C FAUUPT    SOURCE    BECC      10/06/11    21:15:06     6690
      SUBROUTINE FAUUPT(nesp,nsca,idim,
     &     gaml,rhol,pl,unl,utl,uvl,etherl,
     &     gamr,rhor,pr,unr,utr,uvr,etherr,
     &     yl,yr,sl,sr,v_inf,f,
     &     cellt)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FAUUPT ('FLUX for Bas Mach' and thermally perfect)
C
C DESCRIPTION       :
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. Beccantini, DM2S/SFME/LTMF
C
C************************************************************************
C
c----------------------------------------------------------------------
c GENERAL DESCRIPTION:
c   This subroutine provides the numerical flux function
c   defined at the cell interface; this flux is given in
c   the NORMAL DIRECTION (nvect)
c   The low-mach number corrections are made for the flux functions
c
c EQUATIONS: 2D Euler equations of gas dynamics
c
c
c REFERENCE: 1) JCP, 129, 364-382 (1996)
c               " A Sequel to AUSM: AUSM+ ";
c               M.S.Liou
c            2) AIAA Journal, Sept. 1998
c               "Low-Diffusion Flux-Splitting Methods for Flows at
C               All Speeds"
c               J.R.Edwards and M.S.Liou
C            3) JCP, 214, 137-170 (2006)
C               A Sequel to AUSM, Part 2 AUSM+-up for all speeds
c               M.S.Liou
C
c----------------------------------------------------------------------
c INPUT:
c
      implicit integer(i-n)
      integer nesp, nsca, idim,iesp
      real*8 gaml, rhol, pl, unl, utl, uvl, etherl
     &     , gamr, rhor, pr, unr, utr, uvr, etherr
     &     , yl(*), yr(*), v_inf, f(*)
     &     , cellt, sl(*), sr(*)
C     Some functions and parameters
      real*8 fau_m4, fau_p5, kp, ku, sigma
C
      real*8 htotl, htotr, asonl, asonr, qq
     &     , uco_l, uco_r, mco0, asonm, ml, mr, mbar2
     &     , mzero, fson, alpha, rhom, mhalf, uhalf
     &     , fp5p, fp5m, phalf

      PARAMETER(kp=0.25d0, ku=0.75d0, sigma=1.0D0)
C      PARAMETER(kp=0.25d0, ku=0.0d0, sigma=1.0D0)
c------------------------------------------------------------------
c Computation of the specific total enthalpies.
c------------------------------------------------------------------
      htotl = ((unl*unl) + (utl*utl) + (uvl*uvl)) / 2.0d0
      htotl = htotl + etherl
      htotl = htotl + (pl / rhol)
      htotr = ((unr*unr) + (utr*utr) + (uvr*uvr)) / 2.0d0
      htotr = htotr + etherr
      htotr = htotr + (pr / rhor)
*      write(*,*) htotl, htotr
c------------------------------------------------------------------
c   Computing reference velocity
c   see Eq.(70) of the Ref.3).
c------------------------------------------------------------------
      asonl=sqrt(gaml*pl/rhol)
      asonr=sqrt(gamr*pr/rhor)
      qq=0.5D0*sqrt((unl*unl)+(utl*utl)+(uvl*uvl))
c------------------------------------------------------------------
      if(qq .lt. v_inf) then
         uco_l = v_inf
      else
         uco_l = qq
      endif
c-----------------------------
      if(uco_l .ge. asonl) then
         uco_l = asonl
      endif
c-------------------------------------------------------------------
      qq=0.5D0*sqrt((unr*unr)+(utr*utr)+(uvr*uvr))
      if(qq .lt. v_inf) then
         uco_r = v_inf
      else
         uco_r = qq
      endif
c-------------------------
      if(uco_r .ge. asonr) then
         uco_r = asonr
      endif
c-------------------------------------------------------------------
c  mco0 is the equivalent of M_{\inf} in Ref 3.
c-------------------------------------------------------------------
      mco0 = max((uco_r / asonr), (uco_l / asonl))
c      write(*,*) 'mco0 = ', mco0
c-------------------------------------------------------------------
c Computation of the speed of sound;
c numerical speed of sound at the interface is taken as an average
c of the speeds of sounds of the neighbouring cells
c-------------------------------------------------------------------
      asonm=0.5d0*(asonl + asonr)
c-------------------------------------------------------------------
c Computing numerical Mach number
c ml, mr, mbar2 (used in the term for Mu, which does not involves
c the cut-off speed, see equation (70), ref 3.)
c Then, even if mc0 = 1, the term involving kp acts!
c-------------------------------------------------------------------
      ml = unl / asonm
      mr = unr / asonm
      mbar2 = ((unl*unl)+(unr*unr))/(2.0d0*asonm*asonm)
*      write(*,*) asonm
*      write(*,*) ml, mr, mbar2
c-------------------------------
c-------------------------------------------------------------------
c Scaling function for the speed of sound
c We redefine mzero (see (71), ref 3).
c-------------------------------------------------------------------
C
      mzero = max((mco0 * mco0), mbar2)
      mzero = min(1.0D0, mzero)
      mzero = mzero ** 0.5D0
C
      fson = mzero * (2.0D0 - mzero)
      alpha = (3.0D0 / 16.0D0) * (-4.0D0 + (5.0D0 * fson * fson))
C      write(*,*) 'alpha = ', alpha
C-------------------------------------------------------------------
C mhalf (see (73), ref 3)
C First, contribution of Mu; then contribution of ml and mr.
C-------------------------------------------------------------------
C
      rhom = 0.5D0 * (rhol + rhor)
      mhalf = max(1.0D0 - (sigma * mbar2), 0.0D0)
      mhalf = -1.0D0 * mhalf * kp * (pr - pl) /
     &     (fson * rhom * asonm * asonm)
C
      mhalf = mhalf + fau_m4(ml,1.0D0) + fau_m4(mr,-1.0D0)
c      write(*,*) 'mhalf = ', mhalf
      uhalf = mhalf * asonm
c      write(*,*) 'unl, unr ' , unl, unr
c      write(*,*) 'mhalf, uhalf ', mhalf, uhalf
C-------------------------------------------------------------------
C phalf
C-------------------------------------------------------------------
      fp5p = fau_p5(ml,1.0D0,alpha)
      fp5m = fau_p5(mr,-1.0D0,alpha)
      phalf = -1.0D0 * ku * fp5m * fp5p * (rhol + rhor) * fson * asonm *
     &     (unr - unl)
      phalf = phalf + (fp5p * pl) + (fp5m * pr)
c      write(*,*) 'phalf = ', phalf
c---------------------------------------------------------------------
c Computing numerical fluxes
c---------------------------------------------------------------------
      if(uhalf .ge. 0.0d0) then
         f(1) = uhalf * rhol
         f(2) = (f(1) * unl) + phalf
         f(3) = f(1) * utl
         if (idim .eq. 3) then
            f(4) = f(1) * uvl
         endif
         f(2 + idim) = f(1) * htotl
         do iesp = 1, nesp, 1
            f(2 + idim + iesp) = f(1) * yl(iesp)
         enddo
         do iesp = 1, nsca, 1
            f(2 + idim + nesp + iesp) = f(1) * sl(iesp)
         enddo
      else
         f(1) = uhalf * rhor
         f(2) = (f(1) * unr) + phalf
         f(3) = f(1) * utr
         if (idim .eq. 3) then
            f(4) = f(1) * uvr
         endif
         f(idim + 2) = f(1) * htotr
         do iesp = 1, nesp, 1
            f(2 + idim + iesp) = f(1) * yr(iesp)
         enddo
         do iesp = 1, nsca, 1
            f(2 + idim + nesp + iesp) = f(1) * sr(iesp)
         enddo
      endif
c----------------------------------------------------------------------
      cellt = max(abs(unl),abs(unr)) + asonm
c       write(*,*) 'cellt  = ', cellt
c----------------------------------------------------------------------
      return
      end