C GMAT2D    SOURCE    CHAT      05/01/13    00:18:34     5004
       SUBROUTINE GMAT2D(GMT,WVECT,GA,V_INF,DIAM)
C
C      Voir kon171.eso
C
c-----------------------------------------------------
c  EQUATIONS:  2D Euler Equations of Gas Dynamics
c
c  REFERENCE:  AIAA Journal, Sept. 1998
c              "Low-Diffusion Flux-Splitting Methods
c               for Flows at All Speeds"
c-----------------------------------------------------
c  This subroutine provides the preconditioning matrix
c  which is the product of the two matrices:
c     a)  the 'original' preconditioning matrix
c          given in the Reference and
c     b)  the jacobian matrix which is
c
c                  dW/dU
c        with W=(p,u,v,T) (T-temperature) and
c             U=(rho,rho*u,rho*v,rho*e).
c-----------------------------------------------------
c  INPUT:
c
c        wvect   --  vector of the primitive variables
c                    (rho, u, v, p);
c
c        ga      --  ratio of the specific heats;
c
c        v_inf   --  reference velocity
c
c        diam    --  diameter of the cell
c
c  OUTPUT:
c
c        gmt     -- preconditioning matrix 4 by 4.
c-----------------------------------------------------
      IMPLICIT INTEGER(I-N)
       real*8 wvect(4),gmt(4,4)
       real*8 ga,gm1,a
       real*8 rho,p,u,v, v_inf
       real*8 rho_p,ur,theta
       real*8 epsil,h,qq,m,t1,t2
       real*8 diam,lambda,mref,mref2,mach,mach2,dtaum1
       parameter(epsil=1.0d0)
c--------------------------------------------------
c  Input
c--------------------------------------------------
       rho=wvect(1)
       u=wvect(2)
       v=wvect(3)
       p=wvect(4)
c--------------------------------------------------
       gm1=ga-1.0d0
       a=sqrt(ga*p/rho)
c--------------------------------------------------
       rho_p=ga/(a*a)
       h=0.5d0*(u*u+v*v)+ga*p/(gm1*rho)
c--------------------------------------------------
c  Reference velocity
c--------------------------------------------------
       qq=sqrt(u*u+v*v)
       if(qq .lt. (epsil*v_inf)) then
          ur = epsil*v_inf
       else
          ur = qq
       endif
       if(ur .gt. a)then
          ur=a
       endif
c----------------------------------------------------
       theta=1.0d0/(ur*ur) - (1.0d0-ga)*rho_p/ga
c----------------------------------------------------
       m=(1.0d0-ga)*(theta-rho_p)
       t1=1.0d0+m*(h-qq*qq)+(theta-rho_p)*ga/rho_p
       t2=h*(theta-rho_p)*((1.0d0-ga)*(h-qq*qq)+ga/rho_p)
c----- We compute u+a------------
       mref=ur/a
       mref2=mref*mref
       mach=qq/a
       mach2=mach*mach
       lambda=(1.0D0-mref2)**2
       lambda=(lambda*mach2)+(4.0D0*mref2)
       lambda=sqrt(lambda)*a
       lambda=((1.0D0+mref2)*qq)+lambda
       lambda=lambda*0.5D0
       dtaum1=lambda/diam
c--------------------------------
       gmt(1,1)=t1 * dtaum1
       gmt(1,2)=u*m * dtaum1
       gmt(1,3)=v*m * dtaum1
       gmt(1,4)=-1.0D0 * m * dtaum1
c-------------------------------
       gmt(2,1)=u*(t1-1.0d0) * dtaum1
       gmt(2,2)=(1.0d0+u*u*m) * dtaum1
       gmt(2,3)=u*v*m * dtaum1
       gmt(2,4)=-1.0D0 * m*u * dtaum1
c-------------------------------
       gmt(3,1)=v*(t1-1.0d0) * dtaum1
       gmt(3,2)=u*v*m * dtaum1
       gmt(3,3)=(1.0d0+v*v*m) * dtaum1
       gmt(3,4)=-1.0D0 * m*v * dtaum1
c-------------------------------
       gmt(4,1)=t2 * dtaum1
       gmt(4,2)=u*h*m * dtaum1
       gmt(4,3)=v*h*m * dtaum1
       gmt(4,4)=(-m*h+1.0d0) * dtaum1
c-------------------------------
       return
       end