C FSTIFF    SOURCE    BECC      11/05/26    21:15:28     6981
      SUBROUTINE FSTIFF(NDIM, NPAR,
     &     PC_L, GAM_L,
     &     PC_R, GAM_R,
     &     PMIN,
     &     RHO_L, P_L, U_L, UT1_L, UT2_L, W_L,
     &     RHO_R, P_R, U_R, UT1_R, UT2_R, W_R,
     &     RHO_B, P_B, U_B, D_B,
     &     RHO_A, P_A, U_A, D_A,
     &     Z,
     &     UINT,
     &     F_LAG, CELLT,
     &     LOGAN, LOGNC, LOGVAC)
*
*************************************************************************
*
* project           :  CAST3M, Europlexus...
*
* name              :  fstiff
*
* description       :  stiffened gas mixture.
*
* language          :  fortran 77
*
* author            :  a. beccantini den/dm2s/sfme/ltmf
*
*************************************************************************
*
* called sub        :
*
* called by         :
*
*************************************************************************
*
***** input
*
*     ndim                =  dimension (1, 2 or 3)
*
*     npar                =  number of parameters involved in the
*                            equations
*
*     pc_, gam_           =  parameters involved in the eos
*
*     rho_, p_, u_, ut1_, ut2_, w_,
*                         =  density, pressure, normal and tangential
*                            velocities, parameter vector
*                            nb dim(w_) = npar
*
*     z                   =  value of the characteristic variable x/t in
*                            which we want to evaluate the flux f_lag
*
***** output
*
*     uint                =  speed of the contact discontinuity
*
*     f_lag               =  ale interfacial flux in (n,t1,t2), i.e.
*                            rho*(un - x/t)         mass
*                            rho*(un - x/t)*un + p  momentum along n
*                            rho*(un - x/t)*ut1     momentum along t1
*                            rho*(un - x/t)*ut2     momentum along t2
*                            rho*(un - x/t)*ht      energy
*                            rho*(un - x/t)*w(1)    parameter 1
*                            ...
*
*     nb
*     according to nkonga, comput methods appl. mech engnr 190, 2000
*     \dep{u}{x} + \dep{f(u)}{x} = 0
*     z = speed of the surface
*     f_z = f - z u
*
*
*     cellt               =  stability condition, i.e.
*
*                            dt/dx < cellt (dimension 1/velocity)
*
*     logan               =  if true, anomaly detected
*
*     lognc               =  if true, the Newton for the Riemann problem
*                            solution did not converge
*
*************************************************************************
*
*     02/12/2009      created
*     25/05/2011      evolution in CAST3M
*
*************************************************************************
*
* n.b.: all variables are declared
*
C      IMPLICIT NONE
      IMPLICIT INTEGER(I-N)
      INTEGER NDIM, NPAR, IPAR
      REAL*8 PMIN
      REAL*8 PC_L, GAM_L, PC_R, GAM_R, Z
     &     , RHO_L, P_L, U_L, UT1_L, UT2_L, W_L(NPAR)
     &     , RHO_R, P_R, U_R, UT1_R, UT2_R, W_R(NPAR)
     &     , F_LAG(NDIM + 3 + NPAR)
     &     , CELLT
     &     , D_L, D_B, D_A, D_R, RHO_B, P_B, U_B
     &     , RHO_A, P_A, U_A
     &     , FLURHO, FLURU, FLURT1, FLURT2, FLURET, COEF
     &     , UINT
      LOGICAL LOGAN, LOGDEB, LOGNC, LOGVAC
      PARAMETER (LOGDEB = .FALSE.)
*
*************************************************************************
******** field-by-field decomposition ***********************************
*************************************************************************
*
      CALL RACSTI(
     &     PC_L, GAM_L,
     &     PC_R, GAM_R,
     &     PMIN,
     &     RHO_L, P_L, U_L,
     &     RHO_R, P_R, U_R,
     &     D_L, D_B, D_A, D_R,
     &     RHO_B, P_B, U_B,
     &     RHO_A, P_A, U_A,
     &     LOGDEB, LOGAN, LOGNC, LOGVAC)
*
      IF (LOGDEB .OR. LOGNC) THEN
*      IF (LOGDEB .OR. LOGNC .OR. LOGVAC) THEN
         WRITE(*,*)
         WRITE(*,*) 'RHO_L, P_L, U_L, D_L'
         WRITE(*,*) RHO_L, P_L, U_L, D_L
         WRITE(*,*) 'RHO_B, P_B, U_B, D_B'
         WRITE(*,*) RHO_B, P_B, U_B, D_B
         WRITE(*,*) 'RHO_A, P_A, U_A, D_A'
         WRITE(*,*) RHO_A, P_A, U_A, D_A
         WRITE(*,*) 'RHO_R, P_R, U_R, D_R'
         WRITE(*,*) RHO_R, P_R, U_R, D_R
         WRITE(*,*)
      ENDIF
      IF (LOGAN) THEN
         WRITE(*,*) 'SUBROUTINE FSTIFF.F'
         WRITE(*,*) 'ANOMALY DETECTED'
         GOTO 9999
      ENDIF
*
*************************************************************************
******** flux ***********************************************************
*************************************************************************
*
      CALL FLUSTI(
     &     PC_L, GAM_L,
     &     PC_R, GAM_R,
     &     Z,
     &     D_L, D_B, D_A, D_R,
     &     RHO_L, P_L, U_L, UT1_L, UT2_L,
     &     RHO_B, P_B, U_B,
     &     RHO_A, P_A, U_A,
     &     RHO_R, P_R, U_R, UT1_R, UT2_R,
     &     UINT,
     &     FLURHO, FLURU, FLURT1, FLURT2, FLURET,
     &     LOGDEB, LOGAN)
      IF (LOGDEB) THEN
C      IF (LOGVAC) THEN
         WRITE(*,*)
         WRITE(*,*) 'FLURHO, FLURU, FLURT1, FLURT2, FLURET'
         WRITE(*,*) FLURHO, FLURU, FLURT1, FLURT2, FLURET
         WRITE(*,*)
         stop
      ENDIF
      IF (LOGAN) THEN
         WRITE(*,*) 'SUBROUTINE FSTIFF.F'
         WRITE(*,*) 'ANOMALY DETECTED'
         GOTO 9999
      ENDIF
*
***** condition for the cfl
*
      CELLT = MAX (ABS(D_L), ABS(D_R))
*
      F_LAG(1) = FLURHO
      F_LAG(2) = FLURU
      F_LAG(3) = FLURT1
      IF (NDIM .EQ. 3) THEN
*        nb if ndim = 1 and npar = 0, index 4 does not
*           exist
         F_LAG(4) = FLURT2
      ENDIF
      F_LAG(2 + NDIM) = FLURET
      COEF = 1.0D0
      IF (FLURHO .LT. 0) THEN
         COEF = -1.0D0
      ENDIF
*
      DO IPAR = 1, NPAR, 1
         F_LAG(2 + NDIM + IPAR) = 0.5D0 * FLURHO * (
     &        ((COEF + 1.0D0) * W_L(IPAR)) +
     &        ((1.0D0 - COEF) * W_R(IPAR))
     &        )
      ENDDO
 9999 CONTINUE
      RETURN
      END