PROGRAM RIEMANN 
C C     This program computes the solution of a 1D  
c     relativistic Riemann problem with 
C     initial data UL if X<0.5 and UR if X>0.5 
C     in the whole spatial domain [0, 1] 
C
CG    GLOBAL DATA (COMMON BLOCKS) REFERENCED:
CG    /GAMMA/, /STATES/, /LS/, /RS/
C
CF    FILE ACCESS:
CF    solution.dat
C
CM    MODULES CALLED:
CM    GEP, RAREF
C
      PROGRAM RIEMAN

      IMPLICIT NONE

C     -------
C     COMMON BLOCKS
C     -------

      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL,
     & RHOR, PR, UR, HR, CSR, VELR, WR
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, WL,
     & RHOR, PR, UR, HR, CSR, VELR, WR

      DOUBLE PRECISION RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL
      COMMON /LS/      RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL

      DOUBLE PRECISION RHORS, URS, HRS, CSRS, VELRS, VSHOCKR
      COMMON /RS/      RHORS, URS, HRS, CSRS, VELRS, VSHOCKR

      DOUBLE PRECISION GAMMA
      COMMON /GAMMA/   GAMMA

C     ---------
C     INTERNAL VARIABLES
C     ---------

      INTEGER         MN, N, I, ILOOP
      PARAMETER      (MN = 400)

      DOUBLE PRECISION TOL, PMIN, PMAX, DVEL1, DVEL2, CHECK

      DOUBLE PRECISION PS, VELS

      DOUBLE PRECISION RHOA(MN), PA(MN), VELA(MN), UA(MN)

      DOUBLE PRECISION XI

      DOUBLE PRECISION RAD(MN), X1, X2, X3, X4, X5, T

C     -------
C     INITIAL STATES
C     -------

      WRITE(*,*) ' ADIABATIC INDEX OF THE GAS: '
      READ (*,*)   GAMMA

      WRITE(*,*) ' TIME FOR THE SOLUTION: '
      READ (*,*)   T

C     ----- C     LEFT STATE C     -----

      WRITE(*,*) ' -- LEFT STATE -- '
      WRITE(*,*) '      PRESSURE     : '
      READ (*,*) PL
      WRITE(*,*) '      DENSITY      : '
      READ (*,*) RHOL
      WRITE(*,*) '      FLOW VELOCITY: '
      READ (*,*) VELL

C     ------
C     RIGHT STATE
C     ------

      WRITE(*,*) ' -- RIGHT STATE -- '
      WRITE(*,*) '      PRESSURE     : '
      READ (*,*) PR
      WRITE(*,*) '      DENSITY      : '
      READ (*,*) RHOR
      WRITE(*,*) '      FLOW VELOCITY: '
      READ (*,*) VELR

C     ------------------------------
C     SPECIFIC INTERNAL ENERGY, SPECIFIC ENTHALPY, SOUND SPEED AND 
C     FLOW LORENTZ FACTORS IN THE INITIAL STATES
C     ------------------------------ 

      UL  = PL/(GAMMA-1.D0)/RHOL
      UR  = PR/(GAMMA-1.D0)/RHOR

      HL  = 1.D0+UL+PL/RHOL
      HR  = 1.D0+UR+PR/RHOR

      CSL = DSQRT(GAMMA*PL/RHOL/HL)
      CSR = DSQRT(GAMMA*PR/RHOR/HR)

      WL  = 1.D0/DSQRT(1.D0-VELL**2)
      WR  = 1.D0/DSQRT(1.D0-VELR**2)

C     --------
C     NUMBER OF POINTS
C     --------

      N   = 400

C     -------------
C     TOLERANCE FOR THE SOLUTION
C     -------------

      TOL = 0.D0

C

      ILOOP = 0

      PMIN  = (PL + PR)/2.D0
      PMAX  = PMIN

 5    ILOOP = ILOOP + 1

      PMIN  = 0.5D0*MAX(PMIN,0.D0)
      PMAX  = 2.D0*PMAX

      CALL GETDVEL(PMIN, DVEL1)

      CALL GETDVEL(PMAX, DVEL2)

      CHECK = DVEL1*DVEL2
      IF (CHECK.GT.0.D0) GOTO 5

C     ---------------------------
C     PRESSURE AND FLOW VELOCITY IN THE INTERMEDIATE STATES
C     ---------------------------

      CALL GETP(PMIN, PMAX, TOL, PS)

      VELS = 0.5D0*(VELLS + VELRS)

C     ---------------
C     SOLUTION ON THE NUMERICAL MESH
C     ---------------

C     -----------
C     POSITIONS OF THE WAVES
C     -----------

      IF (PL.GE.PS) THEN

        X1 = 0.5D0 + (VELL - CSL )/(1.D0 - VELL*CSL )*T
        X2 = 0.5D0 + (VELS - CSLS)/(1.D0 - VELS*CSLS)*T

      ELSE

        X1 = 0.5D0 + VSHOCKL*T
        X2 = X1

      END IF

      X3 = 0.5D0 + VELS*T

      IF (PR.GE.PS) THEN

        X4 = 0.5D0 + (VELS + CSRS)/(1.D0 + VELS*CSRS)*T
        X5 = 0.5D0 + (VELR + CSR )/(1.D0 + VELR*CSR )*T

      ELSE

        X4 = 0.5D0 + VSHOCKR*T
        X5 = X4

      END IF

C     ----------
C     SOLUTION ON THE MESH
C     ----------

      DO 100 I=1,N

        RAD(I) = DFLOAT(I)/DFLOAT(N)

 100  CONTINUE

      DO 120 I=1,N

        IF (RAD(I).LE.X1) THEN

          PA(I)   = PL
          RHOA(I) = RHOL
          VELA(I) = VELL
          UA(I)   = UL

        ELSE IF (RAD(I).LE.X2) THEN

          XI = (RAD(I) - 0.5D0)/T

          CALL RAREF(XI, RHOL,    PL,    UL,    CSL,  VELL,  'L',
     & RHOA(I), PA(I), UA(I),       VELA(I))

        ELSE IF (RAD(I).LE.X3) THEN

          PA(I)   = PS
          RHOA(I) = RHOLS
          VELA(I) = VELS
          UA(I)   = ULS

        ELSE IF (RAD(I).LE.X4) THEN

          PA(I)   = PS
          RHOA(I) = RHORS
          VELA(I) = VELS
          UA(I)   = URS

        ELSE IF (RAD(I).LE.X5) THEN

          XI = (RAD(I) - 0.5D0)/T

          CALL RAREF(XI, RHOR,    PR,    UR,    CSR,  VELR,  'R',
     & RHOA(I), PA(I), UA(I),       VELA(I))

        ELSE

          PA(I)   = PR
          RHOA(I) = RHOR
          VELA(I) = VELR
          UA(I)   = UR

        END IF

 120  CONTINUE

      OPEN (3,FILE='solution.dat',FORM='FORMATTED',STATUS='NEW')

      WRITE(3,150) N, T
 150  FORMAT(I5,1X,F10.5)

      DO 60 I=1,N
        WRITE(3,200) RAD(I),PA(I),RHOA(I),VELA(I),UA(I)
 60   CONTINUE

 200  FORMAT(5(E15.8,1X))

      CLOSE(3)

      STOP
      END

C     ----------
CN    NAME: G E T D V E L
C     ----------

CP    PURPOSE:
CP    COMPUTE THE DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE
CP    STATES FOR GIVEN LEFT AND RIGHT STATES AND PRESSURE
C
CU    USAGE:
CU    CALL GETDVEL( P, DVEL )
C
CA    ARGUMENTS:
CA    NAME / I,IO,O  / TYPE(DIMENSION)
CA       DESCRIPTION
CA    P    / I       / DOUBLEPRECISION
CA       PRESSURE IN THE INTERMEDIATE STATE
CA    DVEL   / O       / DOUBLEPRECISION
CA       DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE SATES
C
CG    GLOBAL DATA (COMMON BLOCKS) REFERENCED:
CG    /GAMMA/, /STATES/
C
CF    FILE ACCESS:
CF    NONE
C
CM    MODULES CALLED:
CM    GETVEL
C
      SUBROUTINE GETDVEL( P, DVEL )

      IMPLICIT NONE

C     -----
C     ARGUMENTS
C     -----

      DOUBLEPRECISION P, DVEL

C     -------
C     COMMON BLOCKS
C     -------

      DOUBLE PRECISION RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL
      COMMON /LS/      RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL

      DOUBLE PRECISION RHORS,URS,HRS,CSRS,VELRS,VSHOCKR
      COMMON /RS/      RHORS,URS,HRS,CSRS,VELRS,VSHOCKR

      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL,
     & RHOR, PR, UR, HR, CSR, VELR, WR
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, WL,
     & RHOR, PR, UR, HR, CSR, VELR, WR

      DOUBLE PRECISION GAMMA
      COMMON /GAMMA/   GAMMA

C     -----
C     LEFT WAVE
C     -----

      CALL GETVEL(P, RHOL, PL, UL,  HL,  CSL,  VELL,  WL, 'L',
     & RHOLS,    ULS, HLS, CSLS, VELLS, VSHOCKL )

C     -----
C     RIGHT WAVE
C     -----

      CALL GETVEL(P, RHOR, PR, UR,  HR,  CSR,  VELR,  WR, 'R',
     & RHORS,    URS, HRS, CSRS, VELRS, VSHOCKR )

      DVEL = VELLS - VELRS

      RETURN
      END

C     -------
CN    NAME: G E T P
C     -------

CP    PURPOSE:
CP    FIND THE PRESSURE IN THE INTERMEDIATE STATE OF A RIEMANN PROBLEM IN
CP    RELATIVISTIC HYDRODYNAMICS
C
CD    DESCRIPTION:
CD    THIS ROUTINE USES A COMBINATION OF INTERVAL BISECTION AND INVERSE
CD    QUADRATIC INTERPOLATION TO FIND THE ROOT IN A SPECIFIED INTERVAL.
C
CU    USAGE:
CU    CALL GETP( PMIN, PMAX, TOL, PS )
C
CA    ARGUMENTS:
CA    NAME  / I,IO,O  / TYPE(DIMENSION)
CA       DESCRIPTION
CA    PMIN    / I       / DOUBLEPRECISION
CA       THE LEFT ENDPOINT OF THE INTERVAL
CA    PMAX    / I       / DOUBLEPRECISION
CA       THE RIGHT ENDPOINT OF THE INTERVAL
CA    TOL   / I       / DOUBLEPRECISION
CA       TOLERANCE:  DESIRED LENGTH OF THE INTERVAL OF UNCERTAINTY
CA       OF THE RESULT ( .GE. 0)
CA    PS    / O       / DOUBLEPRECISION
CA       PRESSURE IN THE INTERMEDIATE STATE WITHIN THE SPECIFIED TOLERANCE 
CA       (OR TO MACHINE PRECISION IF THE TOLERANCE IS ZERO).
C
CG    GLOBAL DATA (COMMON BLOCKS) REFERENCED:
CG    /GAMMA/, /STATES/
C
CF    FILE ACCESS:
CF    NONE
C
CM    MODULES CALLED:
CM    GETDVEL
C
CC    COMMENTS:
CC    IT IS ASSUMED THAT DVEL(PMIN) AND DVEL(PMAX) HAVE OPPOSITE SIGNS WITHOUT
CC    A CHECK.
CC    THIS ROUTINE IS FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATION",
CC    BY G. E. FORSYTHE, M. A. MALCOLM, AND C. B. MOLER,
CC    PRENTICE-HALL, ENGLEWOOD CLIFFS N.J.
C       SUBROUTINE GETP( PMIN, PMAX, TOL, PS )

      IMPLICIT NONE

C     -----
C     ARGUMENTS
C     -----

      DOUBLEPRECISION PMIN, PMAX, TOL, PS

C     -------
C     COMMON BLOCKS
C     -------

      DOUBLEPRECISION GAMMA
      COMMON /GAMMA/  GAMMA

      DOUBLEPRECISION RHOL, PL, UL, HL, CSL, VELL, WL,
     & RHOR, PR, UR, HR, CSR, VELR, WR
      COMMON /STATES/ RHOL, PL, UL, HL, CSL, VELL, WL,
     & RHOR, PR, UR, HR, CSR, VELR, WR

C     ---------
C     INTERNAL VARIABLES
C     ---------

      DOUBLEPRECISION A, B, C, D, E, EPS, FA, FB, FC, TOL1,
     & XM, P, Q, R, S

C     -------------
C     COMPUTE MACHINE PRECISION
C     -------------

      EPS  = 1.D0
10    EPS  = EPS/2.D0
      TOL1 = 1.D0 + EPS
      IF( TOL1 .GT. 1.D0 ) GO TO 10

C     -------
C     INITIALIZATION
C     -------

      A  = PMIN
      B  = PMAX
      CALL GETDVEL(A,FA)
      CALL GETDVEL(B,FB)

C     -----
C     BEGIN STEP
C     -----

20    C  = A
      FC = FA
      D  = B - A
      E  = D
30    IF( DABS(FC) .GE. DABS(FB) )GO TO 40
      A  = B
      B  = C
      C  = A
      FA = FB
      FB = FC
      FC = FA

C     --------
C     CONVERGENCE TEST
C     --------

40    TOL1 = 2.D0*EPS*DABS(B) + 0.5D0*TOL
      XM   = 0.5D0*(C - B)
      IF( DABS(XM) .LE. TOL1 ) GO TO 90
      IF( FB .EQ. 0.D0 ) GO TO 90

C     ------------
C     IS BISECTION NECESSARY?
C     ------------

      IF( DABS(E) .LT. TOL1 ) GO TO 70
      IF( DABS(FA) .LE. DABS(FB) ) GO TO 70

C     ------------------
C     IS QUADRATIC INTERPOLATION POSSIBLE?
C     ------------------

      IF( A .NE. C ) GO TO 50

C     ----------
C     LINEAR INTERPOLATION
C     ----------

      S = FB/FA
      P = 2.D0*XM*S
      Q = 1.D0 - S
      GO TO 60

C     ----------------
C     INVERSE QUADRATIC INTERPOLATION
C     ----------------

50    Q = FA/FC
      R = FB/FC
      S = FB/FA
      P = S*(2.D0*XM*Q*(Q - R) - (B - A)*(R - 1.D0))
      Q = (Q - 1.D0)*(R - 1.D0)*(S - 1.D0)

C     ------
C     ADJUST SIGNS
C     ------

60    IF( P .GT. 0.D0 ) Q = -Q
      P = DABS(P)

C     --------------
C     IS INTERPOLATION ACCEPTABLE?
C     --------------

      IF( (2.D0*P) .GE. (3.D0*XM*Q-DABS(TOL1*Q)) ) GO TO 70
      IF( P .GE. DABS(0.5D0*E*Q) ) GO TO 70
      E = D
      D = P/Q
      GO TO 80

C     -----
C     BISECTION
C     -----

70    D = XM
      E = D

C     -------
C     COMPLETE STEP
C     -------

80    A  = B
      FA = FB
      IF( DABS(D) .GT. TOL1 ) B = B+D
      IF( DABS(D) .LE. TOL1 ) B = B+DSIGN(TOL1,XM)
      CALL GETDVEL(B,FB)
      IF( (FB*(FC/DABS(FC))) .GT. 0.D0) GO TO 20
      GO TO 30

C     --
C     DONE
C     --

90    PS = B

      RETURN
      END

C     ---------
CN    NAME: G E T V E L
C     ---------

CP    PURPOSE:
CP    COMPUTE THE FLOW VELOCITY BEHIND A RAREFACTION OR SHOCK IN TERMS OF THE
CP    POST-WAVE PRESSURE FOR A GIVEN STATE AHEAD THE WAVE IN A RELATIVISTIC
CP    FLOW
C
CD    DESCRIPTION:
C
CU    USAGE:
CU    CALL GETVEL( P, RHOA, PA, UA, HA, CSA, VELA, WA, S, VEL )
C
CA    ARGUMENTS:
CA    NAME / I,IO,O  / TYPE(DIMENSION)
CA       DESCRIPTION
CA    P      / I       / DOUBLEPRECISION
CA       THE POST-WAVE PRESSURE
CA    RHOA   / I       / DOUBLEPRECISION
CA       THE DENSITY AHEAD THE WAVE
CA    PA     / I       / DOUBLEPRECISION
CA       THE PRESSURE AHEAD THE WAVE
CA    UA     / I       / DOUBLEPRECISION
CA       THE SPECIFIC INTERNAL ENERGY AHEAD THE WAVE
CA    HA     / I       / DOUBLEPRECISION
CA       THE SPECIFIC ENTHALPY AHEAD THE WAVE
CA    CSA    / I       / DOUBLEPRECISION
CA       THE LOCAL SOUND SPEED AHEAD THE WAVE
CA    VELA   / I       / DOUBLEPRECISION
CA       THE FLOW VELOCITY AHEAD THE WAVE
CA    WA     / I       / DOUBLEPRECISION
CA       THE FLOW LORENTZ FACTOR AHEAD THE WAVE
CA    S      / I       / CHARACTER
CA       THE DIRECTION OF PROPAGATION OF THE WAVE (LEFT OR RIGHT)
CA    RHO    / O       / DOUBLEPRECISION
CA       THE DENSITY IN THE POST-WAVE STATE
CA    U      / O       / DOUBLEPRECISION
CA       THE SPECIFIC INTERNAL ENERGY IN THE POST-WAVE STATE
CA    H      / O       / DOUBLEPRECISION
CA       THE SPECIFIC ENTHALPY IN THE POST-WAVE STATE
CA    CS     / O       / DOUBLEPRECISION
CA       THE LOCAL SOUND SPEED IN THE POST-WAVE STATE
CA    VEL    / O       / DOUBLEPRECISION
CA       THE FLOW VELOCITY IN THE POST-WAVE STATE
CA    VSHOCK / O       / DOUBLEPRECISION
CA       THE SHOCK VELOCITY
C
CG    GLOBAL DATA (COMMON BLOCKS) REFERENCED:
CG    /GAMMA/
C
CF    FILE ACCESS:
CF    NONE
C
CM    MODULES CALLED:
CM    NONE
C
CC    COMMENTS:
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER,
CC    J. FLUID MECH., (1994)

      SUBROUTINE GETVEL( P, RHOA, PA, UA, HA, CSA, VELA, WA, S,
     & RHO,      U,  H,  CS,  VEL,  VSHOCK )

      IMPLICIT NONE

C     -----
C     ARGUMENTS
C     -----

      DOUBLE PRECISION P, RHOA, PA, UA, HA, CSA, VELA, WA 
      CHARACTER*1      S
      DOUBLE PRECISION RHO, U, H, CS, VEL, VSHOCK

C     -------
C     COMMON BLOCKS
C     -------

      DOUBLE PRECISION GAMMA
      COMMON /GAMMA/   GAMMA

C     ---------
C     INTERNAL VARIABLES
C     ---------

      DOUBLE PRECISION A, B, C, SIGN
      DOUBLE PRECISION J, WSHOCK
      DOUBLE PRECISION K, SQGL1

C     ---------------
C     LEFT OR RIGHT PROPAGATING WAVE
C     ---------------

      IF (S.EQ.'L') SIGN = -1.D0

      IF (S.EQ.'R') SIGN =  1.D0

C

      IF (P.GT.PA) THEN

C       ---
C       SHOCK
C       ---

        A  = 1.D0+(GAMMA-1.D0)*(PA-P)/GAMMA/P
        B  = 1.D0-A
        C  = HA*(PA-P)/RHOA-HA**2

C       ----------------
C       CHECK FOR UNPHYSICAL ENTHALPIES
C       ----------------

        IF (C.GT.(B**2/4.D0/A)) STOP
     & 'GETVEL: UNPHYSICAL SPECIFIC ENTHALPY IN INTERMEDIATE STATE'

C       -----------------------------
C       SPECIFIC ENTHALPY IN THE POST-WAVE STATE
C       (FROM THE EQUATION OF STATE AND THE TAUB ADIABAT,
C       EQ.(74), MM94)
C       -----------------------------

        H = (-B+DSQRT(B**2-4.D0*A*C))/2.D0/A

C       ---------------
C       DENSITY IN THE POST-WAVE STATE
C       (FROM EQ.(73), MM94)
C       ---------------

        RHO = GAMMA*P/(GAMMA-1.D0)/(H-1.D0)

C       ------------------------
C       SPECIFIC INTERNAL ENERGY IN THE POST-WAVE STATE
C       (FROM THE EQUATION OF STATE)
C       ------------------------

        U = P/(GAMMA-1.D0)/RHO

C       --------------------------
C       MASS FLUX ACROSS THE WAVE 
C       (FROM THE RANKINE-HUGONIOT RELATIONS, EQ.(71), MM94)
C       --------------------------

        J = SIGN*DSQRT((P-PA)/(HA/RHOA-H/RHO))

C       ----------
C       SHOCK VELOCITY
C       (FROM EQ.(86), MM94
C       ----------

        A      = J**2+(RHOA*WA)**2
        B      = -VELA*RHOA**2*WA**2
        VSHOCK = (-B+SIGN*J**2*DSQRT(1.D0+RHOA**2/J**2))/A
        WSHOCK = 1.D0/DSQRT(1.D0-VSHOCK**2)

C       -------------------
C       FLOW VELOCITY IN THE POST-SHOCK STATE
C       (FROM EQ.(67), MM94)
C       -------------------

        A = WSHOCK*(P-PA)/J+HA*WA*VELA
        B = HA*WA+(P-PA)*(WSHOCK*VELA/J+1.D0/RHOA/WA)

        VEL = A/B

C       ---------------------
C       LOCAL SOUND SPEED IN THE POST-SHOCK STATE
C       (FROM THE EQUATION OF STATE)
C       ---------------------

        CS = DSQRT(GAMMA*P/RHO/H)

      ELSE

C       ------
C       RAREFACTION
C       ------

C       ---------------------------
C       POLITROPIC CONSTANT OF THE GAS ACROSS THE RAREFACTION
C       ---------------------------

        K = PA/RHOA**GAMMA

C       ---------------
C       DENSITY BEHIND THE RAREFACTION
C       ---------------

        RHO = (P/K)**(1.D0/GAMMA)

C       ------------------------
C       SPECIFIC INTERNAL ENERGY BEHIND THE RAREFACTION
C       (FROM THE EQUATION OF STATE)
C       ------------------------

        U = P/(GAMMA-1.D0)/RHO

C       --------------------
C       LOCAL SOUND SPEED BEHIND THE RAREFACTION
C       (FROM THE EQUATION OF STATE)
C       --------------------

        CS = DSQRT(GAMMA*P/(RHO+GAMMA*P/(GAMMA-1.D0)))

C       ------------------
C       FLOW VELOCITY BEHIND THE RAREFACTION
C       ------------------

        SQGL1 = DSQRT(GAMMA-1.D0)
        A   = (1.D0+VELA)/(1.D0-VELA)*
     & ((SQGL1+CSA)/(SQGL1-CSA)*
     & (SQGL1-CS )/(SQGL1+CS ))**(-SIGN*2.D0/SQGL1)

        VEL = (A-1.D0)/(A+1.D0)

      END IF

      END

C     --------
CN    NAME: R A R E F
C     --------

CP    PURPOSE:
CP    COMPUTE THE FLOW STATE IN A RAREFACTION FOR GIVEN PRE-WAVE STATE
C CD    DESCRIPTION:
C CU    USAGE:
CU    CALL RAREF( XI, RHOA, PA, UA, CSA, VELA, S, RHO, P, U, VEL )
C CA    ARGUMENTS:
CA    NAME / I,IO,O  / TYPE(DIMENSION)
CA       DESCRIPTION
CA    XI   / I       / DOUBLEPRECISION
CA       SIMILARITY VARIABLE
CA    RHOA   / I       / DOUBLEPRECISION
CA       THE DENSITY AHEAD THE RAREFACTION
CA    PA     / I       / DOUBLEPRECISION
CA       THE PRESSURE AHEAD THE RAREFACTION
CA    UA     / I       / DOUBLEPRECISION
CA       THE SPECIFIC INTERNAL ENERGY AHEAD THE RAREFACTION
CA    CSA    / I       / DOUBLEPRECISION
CA       THE LOCAL SOUND SPEED AHEAD THE RAREFACTION
CA    VELA   / I       / DOUBLEPRECISION
CA       THE FLOW VELOCITY AHEAD THE RAREFACTION
CA    S      / I       / CHARACTER
CA       THE DIRECTION OF PROPAGATION OF THE WAVE (LEFT OR RIGHT)
CA    RHO    / O       / DOUBLEPRECISION
CA       THE DENSITY IN A POINT WITHIN THE RAREFACTION
CA    P      / O       / DOUBLEPRECISION
CA       THE PRESSURE IN A POINT WITHIN THE RAREFACTION
CA    U      / O       / DOUBLEPRECISION
CA       THE SPECIFIC INTERNAL ENERGY IN A POINT WITHIN THE RAREFACTION
CA    CS     / O       / DOUBLEPRECISION
CA       THE LOCAL SOUND SPEED IN IN A POINT WITHIN THE RAREFACTION
CA    VEL    / O       / DOUBLEPRECISION
CA       THE FLOW VELOCITY IN IN A POINT WITHIN THE RAREFACTION
C
CG    GLOBAL DATA (COMMON BLOCKS) REFERENCED:
CG    /GAMMA/
C
CF    FILE ACCESS:
CF    NONE
C
CM    MODULES CALLED:
CM    NONE
C
CC    COMMENTS:
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER,
CC    J. FLUID MECH., (1994)

      SUBROUTINE RAREF( XI, RHOA, PA, UA, CSA, VELA, S, RHO, P, U, VEL )

      IMPLICIT NONE

C     -----
C     ARGUMENTS
C     -----

      DOUBLE PRECISION XI

      DOUBLE PRECISION RHOA, PA, UA, CSA, VELA

      CHARACTER        S

      DOUBLE PRECISION RHO, P, U, VEL

C     -------
C     COMMON BLOCKS
C     -------

      DOUBLE PRECISION GAMMA
      COMMON /GAMMA/   GAMMA

C     ---------
C     INTERNAL VARIABLES
C     ---------

      DOUBLE PRECISION B, C, D, K, L, V, OCS2, FCS2, DFDCS2, CS2, SIGN

C     ---------------
C     LEFT OR RIGHT PROPAGATING WAVE
C     ---------------

      IF (S.EQ.'L') SIGN =  1.D0

      IF (S.EQ.'R') SIGN = -1.D0

      B    = DSQRT(GAMMA - 1.D0)
      C    = (B + CSA)/(B - CSA)
      D    = -SIGN*B/2.D0
      K    = (1.D0 + XI)/(1.D0 - XI)
      L    = C*K**D
      V    = ((1.D0 - VELA)/(1.D0 + VELA))**D

      OCS2 = CSA

25    FCS2   = L*V*(1.D0 + SIGN*OCS2)**D*(OCS2 - B) +
     & (1.D0 - SIGN*OCS2)**D*(OCS2 + B)

      DFDCS2 = L*V*(1.D0 + SIGN*OCS2)**D*
     & (1.D0 + SIGN*D*(OCS2 - B)/(1.D0 + SIGN*OCS2)) +
     & (1.D0 - SIGN*OCS2)**D*      & (1.D0 - SIGN*D*(OCS2 + B)/(1.D0 - SIGN*OCS2))

      CS2 = OCS2 - FCS2/DFDCS2

      IF (ABS(CS2 - OCS2)/OCS2.GT.5.E-7)THEN
        OCS2 = CS2
        GOTO 25
      END IF

      VEL = (XI + SIGN*CS2)/(1.D0 + SIGN*XI*CS2)

      RHO = RHOA*((CS2**2*(GAMMA - 1.D0 - CSA**2))/
     & (CSA**2*(GAMMA - 1.D0 - CS2**2)))
     & **(1.D0/(GAMMA - 1.D0))

      P   = CS2**2*(GAMMA - 1.D0)*RHO/(GAMMA - 1.D0 - CS2**2)/GAMMA

      U   = P/(GAMMA - 1.D0)/RHO

      RETURN 
      END