c@a
c@versb
C-----------------------------------------------------------------------
C
CVERS                  Code_Saturne version 1.3
C                      ------------------------
C
C     This file is part of the Code_Saturne Kernel, element of the
C     Code_Saturne CFD tool.
C
C     Copyright (C) 1998-2008 EDF S.A., France
C
C     contact: saturne-support@edf.fr
C
C     The Code_Saturne Kernel is free software; you can redistribute it
C     and/or modify it under the terms of the GNU General Public License
C     as published by the Free Software Foundation; either version 2 of
C     the License, or (at your option) any later version.
C
C     The Code_Saturne Kernel is distributed in the hope that it will be
C     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
C     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with the Code_Saturne Kernel; if not, write to the
C     Free Software Foundation, Inc.,
C     51 Franklin St, Fifth Floor,
C     Boston, MA  02110-1301  USA
C
C-----------------------------------------------------------------------
c@verse
                        SUBROUTINE USCFTH
C                       *****************
C     --------------------------------------------------------------
     & ( IDBIA0 , IDBRA0 ,
     &   NDIM   , NCELET , NCEL   , NFAC   , NFABOR , NFML   , NPRFML ,
     &   NNOD   , LNDFAC , LNDFBR , NCELBR ,
     &   NVAR   , NSCAL  , NPHAS  ,
     &   ICCFTH , IMODIF , IPHAS  ,
     &   NIDEVE , NRDEVE , NITUSE , NRTUSE ,
     &   IFACEL , IFABOR , IFMFBR , IFMCEL , IPRFML ,
     &   IPNFAC , NODFAC , IPNFBR , NODFBR ,
     &   IDEVEL , ITUSER , IA     ,
     &   XYZCEN , SURFAC , SURFBO , CDGFAC , CDGFBO , XYZNOD , VOLUME ,
     &   DT     , RTP    , RTPA   , PROPCE , PROPFA , PROPFB ,
     &   COEFA  , COEFB  ,
     &   SORTI1 , SORTI2 , GAMAGR , XMASMR ,
     &   RDEVEL , RTUSER , RA     )
C     --------------------------------------------------------------
C***********************************************************************
C FONCTION :
C ---------
c@foncb
CFONC
CFONC ROUTINE UTILISATEUR POUR ENTREE DES PARAMETRES THERMODYNAMIQUES
CFONC
CFONC
CFONC     CE SOUS PROGRAMME UTILISATEUR EST OBLIGATOIRE
CFONC     =============================================
CFONC
CFONC
CFONC
CFONC CALCUL DE LA PROPRIETE THERMODYNAMIQUE VOULUE
CFONC EN FONCTION DES PARAMETRES DONNES ET DE LA THERMO CHOISIE
CFONC
CFONC
CFONC LOIS THERMODYNAMIQUES IMPLEMENTEES :
CFONC ====================================
CFONC
CFONC   1. GAZ PARFAIT (fournir la masse molaire XMASML)
CFONC
CFONC   2. GAZ PARFAIT A GAMMA VARIABLE (à adapter)
CFONC
CFONC   3. VAN DER WAALS (pas encore implemente)
CFONC
CFONC
CFONC CALCULS IMPLEMENTES :
CFONC =====================
CFONC
CFONC   VARIABLE     P   rho   T    e    h    s  epsilon-CvT  C.L.
CFONC       CODE     1    2    3    4    5    6      7         9
CFONC
CFONC
CFONC   CODE DU CALCUL : *VARIABLE i : i
CFONC
CFONC                                       di|
CFONC                    *DERIVEE PARTIELLE --| : ijk
CFONC                                       dj|k
CFONC
CFONC                    *CALCUL DES VARIABLES A PARTIR DE i ET j : ij
CFONC
CFONC   Et pour un reperage automatique, on utilise les nombres premiers
CFONC            P   rho   T   e   s
CFONC            2     3   5   7   13  (* 10000 pour les calculs aux cellules)
CFONC
CFONC
CFONC
CFONC   -> OPTIONS DE CALCUL                           : ICCFTH = -1
CFONC
CFONC   -> INITIALISATIONS PAR DEFAUT                  : ICCFTH = 0
CFONC
CFONC   -> CALCUL DE GAMMA                             : ICCFTH = 1
CFONC
CFONC   -> VERIFICATION DE rho                         : ICCFTH = -2
CFONC
CFONC   -> VERIFICATION DE E                           : ICCFTH = -4
CFONC
CFONC   -> CALCUL DE T ET e EN FONCTION DE P ET rho    : ICCFTH = 12 ou  60000
CFONC
CFONC   -> CALCUL DE rho ET e EN FONCTION DE P ET T    : ICCFTH = 13 ou 100000
CFONC
CFONC   -> CALCUL DE rho ET T EN FONCTION DE P ET e    : ICCFTH = 14 ou 140000
CFONC
CFONC   -> CALCUL DE P ET e EN FONCTION DE rho ET T    : ICCFTH = 23 ou 150000
CFONC
CFONC   -> CALCUL DE P ET T EN FONCTION DE rho ET e    : ICCFTH = 24 ou 210000
CFONC
CFONC                 2    dP |
CFONC   -> CALCUL DE c  = ----|                        : ICCFTH = 126
CFONC                     drho|s
CFONC
CFONC                       dP|
CFONC   -> CALCUL DE beta = --|                        : ICCFTH = 162
CFONC                       ds|rho
CFONC
CFONC                     de|
CFONC   -> CALCUL DE Cv = --|                          : ICCFTH = 432
CFONC                     dT|rho
CFONC
CFONC   -> CALCUL DE L'ENTROPIE                        : ICCFTH = 6
CFONC
CFONC
CFONC   -> CALCUL DE epsilon - Cv.T                    : ICCFTH = 7
CFONC
CFONC
CFONC   -> CALCUL DES CONDITIONS AUX LIMITES
CFONC               - SYMETRIE                         : ICCFTH = 90
CFONC               - PAROI                            : ICCFTH = 91
CFONC               - ENTREE                           : ICCFTH = 92
CFONC               - SORTIE                           : ICCFTH = 93
CFONC               - T ET e EN FONCTION DE P ET rho   : ICCFTH = 912 ou  60900
CFONC               - rho ET e EN FONCTION DE P ET T   : ICCFTH = 913 ou 100900
CFONC               - rho ET T EN FONCTION DE P ET e   : ICCFTH = 914 ou 140900
CFONC               - P ET e EN FONCTION DE rho ET T   : ICCFTH = 923 ou 150900
CFONC               - P ET T EN FONCTION DE rho ET e   : ICCFTH = 924 ou 210900
CFONC
CFONC
c@fonce
C-----------------------------------------------------------------------
C                                 ARGUMENTS
c@argub
CARGU .______________.____._____.______________________________________.
CARGU !    NOM       !TYPE!MODE !                   ROLE               !
CARGU !______________!____!_____!______________________________________!
CARGU ! IDBIA0       ! E  !  -> ! NUMERO DE LA 1ERE CASE LIBRE DANS IA !
CARGU ! IDBRA0       ! E  !  -> ! NUMERO DE LA 1ERE CASE LIBRE DANS RA !
CARGU ! NDIM         ! E  !  -> ! DIMENSION DE L'ESPACE                !
CARGU ! NCELET       ! E  !  -> ! NOMBRE D'ELEMENTS HALO COMPRIS       !
CARGU ! NCEL         ! E  !  -> ! NOMBRE D'ELEMENTS ACTIFS             !
CARGU ! NFAC         ! E  !  -> ! NOMBRE DE FACES INTERNES             !
CARGU ! NFABOR       ! E  !  -> ! NOMBRE DE FACES DE BORD              !
CARGU ! NFML         ! E  !  -> ! NOMBRE DE FAMILLES D ENTITES         !
CARGU ! NPRFML       ! E  !  -> ! NOMBRE DE PROPRIETESE DES FAMILLES   !
CARGU ! NNOD         ! E  !  -> ! NOMBRE DE SOMMETS                    !
CARGU ! LNDFAC       ! E  !  -> ! LONGUEUR DU TABLEAU NODFAC (OPTIONNEL!
CARGU ! LNDFBR       ! E  !  -> ! LONGUEUR DU TABLEAU NODFBR (OPTIONNEL!
CARGU ! NCELBR       ! E  !  -> ! NOMBRE D'ELEMENTS AYANT AU MOINS UNE !
CARGU !              !    !     ! FACE DE BORD                         !
CARGU ! NVAR         ! E  !  -> ! NOMBRE TOTAL DE VARIABLES            !
CARGU ! NSCAL        ! E  !  -> ! NOMBRE TOTAL DE SCALAIRES            !
CARGU ! NPHAS        ! E  !  -> ! NOMBRE DE PHASES                     !
CARGU ! ICCFTH       ! E  !  -> ! NUMERO DU TYPE DE CALCUL DEMANDE     !
CARGU ! IMODIF       ! E  !  -> ! SI CALCUL AUX CELLULES : MODIF DE RTP!
CARGU !              !    !     !   QUAND IMODIF > 0 (EN PARTICULIER   !
CARGU !              !    !     !   AUTOMATISE L'INIT (VOIR USCFXI)    !
CARGU !              !    !     ! SI CALCUL AUX FACES : NUMERO DE FACE !
CARGU !              !    !     !   AUTOMATISE L'INIT (VOIR USCFXI)    !
CARGU ! IPHAS        ! E  !  -> ! NUMERO DE LA PHASE COURANTE          !
CARGU ! NIDEVE NRDEVE! E  !  -> ! LONGUEUR DE IDEVEL RDEVEL            !
CARGU ! NITUSE NRTUSE! E  !  -> ! LONGUEUR DE ITUSER RTUSER            !
CARGU ! IFACEL       ! TE !  -> ! ELEMENTS VOISINS D'UNE FACE INTERNE  !
CARGU ! (2, NFAC)    !    !     !                                      !
CARGU ! IFABOR       ! TE !  -> ! ELEMENT  VOISIN  D'UNE FACE DE BORD  !
CARGU ! (NFABOR)     !    !     !                                      !
CARGU ! IFMFBR       ! TE !  -> ! NUMERO DE FAMILLE D'UNE FACE DE BORD !
CARGU ! (NFABOR)     !    !     !                                      !
CARGU ! IFMCEL       ! TE !  -> ! NUMERO DE FAMILLE D'UNE CELLULE      !
CARGU ! (NCELET)     !    !     !                                      !
CARGU ! IPRFML       ! TE !  -> ! PROPRIETES D'UNE FAMILLE             !
CARGU ! NFML  ,NPRFML!    !     !                                      !
CARGU ! IPNFAC       ! TE !  -> ! POSITION DU PREMIER NOEUD DE CHAQUE  !
CARGU !   (LNDFAC)   !    !     !  FACE INTERNE DANS NODFAC (OPTIONNEL)!
CARGU ! NODFAC       ! TE !  -> ! CONNECTIVITE FACES INTERNES/NOEUDS   !
CARGU !   (NFAC+1)   !    !     !  (OPTIONNEL)                         !
CARGU ! IPNFBR       ! TE !  -> ! POSITION DU PREMIER NOEUD DE CHAQUE  !
CARGU !   (LNDFBR)   !    !     !  FACE DE BORD DANS NODFBR (OPTIONNEL)!
CARGU ! NODFBR       ! TE !  -> ! CONNECTIVITE FACES DE BORD/NOEUDS    !
CARGU !   (NFABOR+1) !    !     !  (OPTIONNEL)                         !
CARGU ! IDEVEL(NIDEVE! TE ! <-> ! TAB ENTIER COMPLEMENTAIRE DEVELOPEMT !
CARGU ! ITUSER(NITUSE! TE ! <-> ! TAB ENTIER COMPLEMENTAIRE UTILISATEUR!
CARGU ! IA(*)        ! TR !  -  ! MACRO TABLEAU ENTIER                 !
CARGU ! XYZCEN       ! TR !  -> ! POINT ASSOCIES AUX VOLUMES DE CONTROL!
CARGU ! (NDIM,NCELET !    !     !                                      !
CARGU ! SURFAC       ! TR !  -> ! VECTEUR SURFACE DES FACES INTERNES   !
CARGU ! (NDIM,NFAC)  !    !     !                                      !
CARGU ! SURFBO       ! TR !  -> ! VECTEUR SURFACE DES FACES DE BORD    !
CARGU ! (NDIM,NFABOR)!    !     !                                      !
CARGU ! CDGFAC       ! TR !  -> ! CENTRE DE GRAVITE DES FACES INTERNES !
CARGU ! (NDIM,NFAC)  !    !     !                                      !
CARGU ! CDGFBO       ! TR !  -> ! CENTRE DE GRAVITE DES FACES DE BORD  !
CARGU ! (NDIM,NFABOR)!    !     !                                      !
CARGU ! XYZNOD       ! TR !  -> ! COORDONNES DES NOEUDS (OPTIONNEL)    !
CARGU ! (NDIM,NNOD)  !    !     !                                      !
CARGU ! VOLUME       ! TR !  -> ! VOLUME D'UN DES NCELET ELEMENTS      !
CARGU ! (NCELET      !    !     !                                      !
CARGU ! DT(NCELET)   ! TR !  -> ! VALEUR DU PAS DE TEMPS               !
CARGU ! RTP,RTPA     ! TR ! <-> ! VARIABLES DE CALCUL AU CENTRE DES    !
CARGU ! (NCELET,*)   !    !     !    CELLULES                          !
CARGU ! PROPCE       ! TR ! <-> ! PROPRIETES PHYSIQUES AU CENTRE DES   !
CARGU ! (NCELET,*)   !    !     !    CELLULES                          !
CARGU ! PROPFA       ! TR ! <-> ! PROPRIETES PHYSIQUES AU CENTRE DES   !
CARGU !  (NFAC,*)    !    !     !    FACES INTERNES                    !
CARGU ! PROPFB       ! TR ! <-> ! PROPRIETES PHYSIQUES AU CENTRE DES   !
CARGU !  (NFABOR,*)  !    !     !    FACES DE BORD                     !
CARGU ! COEFA COEFB  ! TR ! <-> ! CONDITIONS AUX LIMITES AUX           !
CARGU !  (NFABOR,*)  !    !     !    FACES DE BORD                     !
CARGU ! SORTI1,2(*)  ! TR ! <-  ! VARIABLES DE SORTIE                  !
CARGU !              !    !     !   (INUTILISE SI ICCFTH.LT.0)         !
CARGU ! GAMAGR(*)    ! TR ! <-  ! CONSTANTE GAMMA EQUIVALENT DU GAZ    !
CARGU !              !    !     !   (INUTILISE SI ICCFTH.LT.0)         !
CARGU !              !    !     !   (PREMIERE CASE UTILISEE EN G.P.)   !
CARGU ! XMASMR(*)    ! TR ! <-  ! MASSE MOLAIRE DES CONSTITUANTS DU GAZ!
CARGU !              !    !     !   (INUTILISE SI ICCFTH.LT.0)         !
CARGU ! RDEVEL(NRDEVE! TR ! <-> ! TAB REEL COMPLEMENTAIRE DEVELOPEMT   !
CARGU ! RTUSER(NRTUSE! TR ! <-> ! TAB REEL COMPLEMENTAIRE UTILISATEUR  !
CARGU ! RA(*)        ! TR !  -  ! MACRO TABLEAU REEL                   !
CARGU !______________!____!_____!______________________________________!
c@argue
C
c@commb
CCOMM                             COMMONS
CCOMM .______________.____._____.______________________________________.
CCOMM !    NOM       !TYPE!MODE !                   ROLE               !
CCOMM !______________!____!_____!______________________________________!
CCOMM !______________!____!_____!______________________________________!
c@comme
C
C     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
C            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
C     MODE : -> DONNEE, <- RESULTAT, <-> DONNEE MODIFIEE,
C            - TABLEAU DE TRAVAIL
C***********************************************************************
C
      IMPLICIT NONE
C
C***********************************************************************
C     DONNEES EN COMMON
C***********************************************************************
C
      INCLUDE "paramx.h"
      INCLUDE "numvar.h"
      INCLUDE "optcal.h"
      INCLUDE "cstphy.h"
      INCLUDE "cstnum.h"
      INCLUDE "pointe.h"
      INCLUDE "entsor.h"
      INCLUDE "parall.h"
      INCLUDE "ppppar.h"
      INCLUDE "ppthch.h"
      INCLUDE "ppincl.h"
C
C***********************************************************************
C
C ARGUMENTS
C
      INTEGER          IDBIA0 , IDBRA0
      INTEGER          NDIM   , NCELET , NCEL   , NFAC   , NFABOR
      INTEGER          NFML   , NPRFML
      INTEGER          NNOD   , LNDFAC , LNDFBR , NCELBR
      INTEGER          NVAR   , NSCAL  , NPHAS
      INTEGER          ICCFTH   , IMODIF , IPHAS
      INTEGER          NIDEVE , NRDEVE , NITUSE , NRTUSE
C
      INTEGER          IFACEL(2,NFAC) , IFABOR(NFABOR)
      INTEGER          IFMFBR(NFABOR) , IFMCEL(NCELET)
      INTEGER          IPRFML(NFML,NPRFML)
      INTEGER          IPNFAC(NFAC+1), NODFAC(LNDFAC)
      INTEGER          IPNFBR(NFABOR+1), NODFBR(LNDFBR)
      INTEGER          IDEVEL(NIDEVE), ITUSER(NITUSE), IA(*)
C
      DOUBLE PRECISION XYZCEN(NDIM,NCELET)
      DOUBLE PRECISION SURFAC(NDIM,NFAC), SURFBO(NDIM,NFABOR)
      DOUBLE PRECISION CDGFAC(NDIM,NFAC), CDGFBO(NDIM,NFABOR)
      DOUBLE PRECISION XYZNOD(NDIM,NNOD), VOLUME(NCELET)
      DOUBLE PRECISION DT(NCELET), RTP(NCELET,*), RTPA(NCELET,*)
      DOUBLE PRECISION PROPCE(NCELET,*),PROPFA(NFAC,*),PROPFB(NFABOR,*)
      DOUBLE PRECISION COEFA(NFABOR,*), COEFB(NFABOR,*)
C
      DOUBLE PRECISION SORTI1(*), SORTI2(*), GAMAGR(*), XMASMR(*)
      DOUBLE PRECISION RDEVEL(NRDEVE), RTUSER(NRTUSE), RA(*)
C
C VARIABLES LOCALES
C
      INTEGER          IDEBIA, IDEBRA
      INTEGER          IIPH   , IFAC0
      INTEGER          IERR
      INTEGER          IEL    , IFAC   , IVAR
      INTEGER          IPRIPH , IRHIPH , ITKIPH , IENIPH
      INTEGER          IUIPH  , IVIPH  , IWIPH
      INTEGER          ICLP   , ICLR   , ICLT   , ICLE
      INTEGER          ICLU   , ICLV   , ICLW
      INTEGER          IUTILE
      DOUBLE PRECISION GAMAGP , XMASML , ENINT
      DOUBLE PRECISION XMACH  , XMACHI , XMACHE , DXMACH
C
      INTEGER          NPMAX
      PARAMETER (NPMAX = 1000)
      DOUBLE PRECISION CSTGR(NPMAX)
C
C***********************************************************************
C
C TEST_A_ENLEVER_POUR_UTILISER_LE_SOUS_PROGRAMME_DEBUT
C=======================================================================
C
      IF(1.EQ.1) RETURN
C
C=======================================================================
C TEST_A_ENLEVER_POUR_UTILISER_LE_SOUS_PROGRAMME_FIN
C
C=======================================================================
C 0. INITIALISATION
C     Pas d'intervention utilisateur requise.
C=======================================================================
C
C     Pointeurs mémoire
      IDEBIA = IDBIA0
      IDEBRA = IDBRA0
C
C     Indicateur d'erreur (stop si non nul)
      IERR   = 0
C
C     Position des variables
      IF(ICCFTH.GE.0.OR.ICCFTH.LE.-2) THEN
        IPRIPH = IPR(IPHAS)
        IRHIPH = ISCA(IRHO  (IPHAS))
        ITKIPH = ISCA(ITEMPK(IPHAS))
        IENIPH = ISCA(IENERG(IPHAS))
        IUIPH = IU(IPHAS)
        IVIPH = IV(IPHAS)
        IWIPH = IW(IPHAS)
        ICLP = ICLRTP(IPRIPH,ICOEF)
        ICLR = ICLRTP(IRHIPH,ICOEF)
        ICLT = ICLRTP(ITKIPH,ICOEF)
        ICLE = ICLRTP(IENIPH,ICOEF)
        ICLU = ICLRTP(IUIPH,ICOEF)
        ICLV = ICLRTP(IVIPH,ICOEF)
        ICLW = ICLRTP(IWIPH,ICOEF)
      ENDIF
C
C
C     Si calcul aux cellules,
C       indique, si > 0, que RTP doit etre modifié
C     Si calcul aux faces,
C       indique le numéro de la face considérée
      IFAC0 = IMODIF
C
C=======================================================================
C 1. CHOIX DE LA THERMODYNAMIQUE POUR CHAQUE PHASE
C     Intervention utilisateur requise.
C=======================================================================
C
C --> IEOS = 1 : Gaz Parfait avec Gamma constant
C --> IEOS = 2 : Gaz Parfait avec Gamma variable
C --> IEOS = 3 : Equation d'etat de Van Der Waals
C
      DO IIPH = 1, NPHAS
        IEOS(IIPH) = 1
      ENDDO
C
C
C            NE PAS OUBLIER DE RENSEIGNER LES PARAMETRES
C                  DANS LA THERMODYNAMIQUE CHOISIE
C            ===========================================
C
C
C=======================================================================
C 2.  GAZ PARFAIT
C=======================================================================
C
      IF(IEOS(IPHAS).EQ.1) THEN
C
C=======================================================================
C 2.1. PARAMETRES A RENSEIGNER PAR L'UTILISATEUR
C     Intervention utilisateur requise.
C=======================================================================
C
C     Pour chaque phase
C
C     Masse molaire du gaz (en kg/mol)
C
        IF(ICCFTH.GE.0) THEN
C
          IF(IPHAS.EQ.1) THEN
C
            XMASML = 28.8D-3
C
          ENDIF
C
        ENDIF
C
C=======================================================================
C 2.2. LOIS IMPLEMENTEES PAR DEFAUT
C     Pas d'intervention utilisateur requise.
C=======================================================================
C
C     CALCUL DE LA CONSTANTE GAMAGP
C     =============================
C
C     On la suppose supérieure ou égale à 1.
C     On la calcule à chaque appel, même si cela peut paraître cher, dans
C       par cohérence avec le cas gamma variable, pour lequel, on n'a pas
C       prévu de la conserver. Avec un save et un test, on pourrait
C       s'affranchir du calcul si nécessaire.
C
        IF(ICCFTH.GT.0) THEN
C
          GAMAGP = 1.D0 + RR/(XMASML*CP0(IPHAS)-RR)
C
          IF(GAMAGP.LT.1.D0) THEN
            WRITE(NFECRA,1010) GAMAGP
            CALL CSEXIT (1)
          ENDIF
C
C     On renvoie gamma si demande
          IF(ICCFTH.EQ.1) THEN
            GAMAGR(1) = GAMAGP
          ENDIF
C
        ENDIF
C
C
C     OPTIONS DE CALCUL  : Cp ET Cv CONSTANTS (hypothèse gaz parfait)
C     =================
C
C     La valeur de CP0 est à fournir dans usini1. La valeur de CV0
C       est calculée plus bas (à partir de CP0).
C
C
        IF(ICCFTH.EQ.-1) THEN
C
          ICP(IPHAS) = 0
          ICV(IPHAS) = 0
C
C
C     INITIALISATIONS PAR DEFAUT (avant uscfxi)
C     ==========================
C     T0 est forcément positif (vérifié dans verini)
C
        ELSEIF(ICCFTH.EQ.0) THEN
C
          CV0(IPHAS) = CP0(IPHAS) - RR/XMASML
C
          IF ( ISUITE .EQ. 0 ) THEN
            DO IEL = 1, NCEL
              RTP(IEL,IRHIPH) = P0(IPHAS)*XMASML/(RR*T0(IPHAS))
              RTP(IEL,IENIPH) = CV0(IPHAS)*T0(IPHAS)
            ENDDO
          ENDIF
C
C
C     VERIFICATION DE rho
C     ===================
C
        ELSEIF(ICCFTH.EQ.-2) THEN
C
C
C --- Clipping si rho est exactement nul + write + stop
C       En effet, si c'est le cas, on va nécessairement se planter
C       dans la thermo (bon, d'accord, pas dans tous les cas, mais
C       pour le moment, je présume que les fluides à traiter seront
C       assez classiques, avec des pbs à rho = 0)
C       Appelé en fin de resolution de la masse volumique (apres
C       clipping classique et avant communication parallele).
C
          IERR = 0
          DO IEL = 1, NCEL
            IF(RTP(IEL,IRHIPH).LE.0.D0) THEN
              RTP(IEL,IRHIPH) = EPZERO
              IERR = IERR + 1
            ENDIF
          ENDDO
          IF(IRANGP.GE.0) THEN
            CALL PARCPT (IERR)
          ENDIF
          IF(IERR.GT.0) THEN
            NTMABS = NTCABS
            WRITE(NFECRA,8000)IERR, EPZERO
          ENDIF
C
C
C     VERIFICATION DE e
C     ===================
C
        ELSEIF(ICCFTH.EQ.-4) THEN
C
C
C --- Clipping si e est exactement nul + write + stop
C       En effet, si c'est le cas, on va nécessairement se planter
C       dans la thermo (bon, d'accord, pas dans tous les cas, mais
C       pour le moment, je présume que les fluides à traiter seront
C       assez classiques, avec des pbs à e = 0)
C
          IERR = 0
          DO IEL = 1, NCEL
            ENINT = RTP(IEL,IENIPH)
     &               - 0.5D0*( RTP(IEL,IUIPH)**2
     &                       + RTP(IEL,IVIPH)**2
     &                       + RTP(IEL,IWIPH)**2 )
            IF(ENINT.LE.0.D0) THEN
              RTP(IEL,IENIPH) = EPZERO
     &               + 0.5D0*( RTP(IEL,IUIPH)**2
     &                       + RTP(IEL,IVIPH)**2
     &                       + RTP(IEL,IWIPH)**2 )
              IERR = IERR + 1
            ENDIF
          ENDDO
          IF(IRANGP.GE.0) THEN
            CALL PARCPT (IERR)
          ENDIF
          IF(IERR.GT.0) THEN
            NTMABS = NTCABS
            WRITE(NFECRA,8100)IERR, EPZERO
          ENDIF
C
C
C     CALCUL DE T ET e EN FONCTION DE P ET rho :
C     ========================================
C
        ELSEIF(ICCFTH.EQ.12.OR.ICCFTH.EQ.60000) THEN
C
C     Test des valeurs de rho
          IERR = 0
          DO IEL = 1, NCEL
            IF(RTP(IEL,IRHIPH).LE.0.D0) THEN
              WRITE(NFECRA,3010)RTP(IEL,IRHIPH),IEL
            ENDIF
          ENDDO
C     Si pb on s'arrete, car rho a été fourni par l'utilisateur
C       (initialisation erronée sans doute)
          IF(IERR.EQ.1) THEN
            CALL CSEXIT (1)
          ENDIF
C
          DO IEL = 1, NCEL
C     Temperature
            SORTI1(IEL) = XMASML*RTP(IEL,IPRIPH)/(RR*RTP(IEL,IRHIPH))
C     Energie totale
            SORTI2(IEL) = CV0(IPHAS)*SORTI1(IEL)
     &           + 0.5D0*( RTP(IEL,IUIPH)**2 + RTP(IEL,IVIPH)**2
     &                                       + RTP(IEL,IWIPH)**2 )
          ENDDO
C
C     Affectation a RTP
          IF(IMODIF.GT.0) THEN
            DO IEL = 1, NCEL
              RTP(IEL,ITKIPH) = SORTI1(IEL)
              RTP(IEL,IENIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C     CALCUL DE rho ET e EN FONCTION DE P ET T :
C     ========================================
C
        ELSEIF(ICCFTH.EQ.13.OR.ICCFTH.EQ.100000) THEN
C
C     Test des valeurs de T
          IERR = 0
          DO IEL = 1, NCEL
            IF(RTP(IEL,ITKIPH).LE.0.D0) THEN
              WRITE(NFECRA,2010)RTP(IEL,ITKIPH),IEL
            ENDIF
          ENDDO
C     Si pb on s'arrete, car T a été fourni par l'utilisateur
C       (initialisation erronée par exemple, non en K ?)
          IF(IERR.EQ.1) THEN
            CALL CSEXIT (1)
          ENDIF
C
          DO IEL = 1, NCEL
C     Masse volumique
            SORTI1(IEL) = XMASML*RTP(IEL,IPRIPH)/(RR*RTP(IEL,ITKIPH))
C     Energie totale
            SORTI2(IEL) = CV0(IPHAS)*RTP(IEL,ITKIPH)
     &           + 0.5D0*( RTP(IEL,IUIPH)**2 + RTP(IEL,IVIPH)**2
     &                                       + RTP(IEL,IWIPH)**2 )
          ENDDO
C
C     Affectation a RTP
          IF(IMODIF.GT.0) THEN
            DO IEL = 1, NCEL
              RTP(IEL,IRHIPH) = SORTI1(IEL)
              RTP(IEL,IENIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C     CALCUL DE rho ET T EN FONCTION DE P ET e :
C     ========================================
C
        ELSEIF(ICCFTH.EQ.14.OR.ICCFTH.EQ.140000) THEN
C
          DO IEL = 1, NCEL
C     Energie interne (évite de diviser par T)
            ENINT = RTP(IEL,IENIPH)
     &               - 0.5D0*( RTP(IEL,IUIPH)**2
     &                       + RTP(IEL,IVIPH)**2
     &                       + RTP(IEL,IWIPH)**2 )
C     Masse volumique
            SORTI1(IEL) = RTP(IEL,IPRIPH) / ( (GAMAGP-1.D0) * ENINT )
C     Temperature
            SORTI2(IEL) = XMASML * (GAMAGP-1.D0) * ENINT / RR
          ENDDO
C
C     Affectation a RTP
          IF(IMODIF.GT.0) THEN
            DO IEL = 1, NCEL
              RTP(IEL,IRHIPH) = SORTI1(IEL)
              RTP(IEL,ITKIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C     CALCUL DE P ET e EN FONCTION DE rho ET T :
C     ========================================
C
        ELSEIF(ICCFTH.EQ.23.OR.ICCFTH.EQ.150000) THEN
C
          DO IEL = 1, NCEL
C     Pression
            SORTI1(IEL) = RTP(IEL,IRHIPH)*RTP(IEL,ITKIPH)*RR/XMASML
C     Energie totale
            SORTI2(IEL) = CV0(IPHAS)*RTP(IEL,ITKIPH)
     &           + 0.5D0*( RTP(IEL,IUIPH)**2 + RTP(IEL,IVIPH)**2
     &                                       + RTP(IEL,IWIPH)**2 )
          ENDDO
C
C     Affectation a RTP
          IF(IMODIF.GT.0) THEN
            DO IEL = 1, NCEL
              RTP(IEL,IPRIPH) = SORTI1(IEL)
              RTP(IEL,IENIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C     CALCUL DE P ET T EN FONCTION DE rho ET e :
C     ========================================
C
        ELSEIF(ICCFTH.EQ.24.OR.ICCFTH.EQ.210000) THEN
C
          DO IEL = 1, NCEL
C     Energie interne (évite de diviser par T)
            ENINT = RTP(IEL,IENIPH)
     &               - 0.5D0*( RTP(IEL,IUIPH)**2
     &                       + RTP(IEL,IVIPH)**2
     &                       + RTP(IEL,IWIPH)**2 )
C     Pression
            SORTI1(IEL) = (GAMAGP-1.D0) * RTP(IEL,IRHIPH) * ENINT
C     Temperature
            SORTI2(IEL) = XMASML * (GAMAGP-1.D0) * ENINT / RR
          ENDDO
C
C       Affectation a RTP
          IF(IMODIF.GT.0) THEN
            DO IEL = 1, NCEL
              RTP(IEL,IPRIPH) = SORTI1(IEL)
              RTP(IEL,ITKIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C                2                            2          P
C     CALCUL DE c  EN FONCTION DE P ET rho : c  = gamma*---
C     ====================================              rho
C
        ELSEIF(ICCFTH.EQ.126) THEN
C
C     Test des valeurs de rho (on peut éliminer ce test pour optimiser)
          IUTILE = 0
          IF(IUTILE.EQ.1) THEN
            IERR = 0
            DO IEL = 1, NCEL
              IF(RTP(IEL,IRHIPH).LE.0.D0) THEN
                WRITE(NFECRA,4010)RTP(IEL,IRHIPH),IEL
              ENDIF
            ENDDO
            IF(IERR.EQ.1) THEN
              CALL CSEXIT (1)
            ENDIF
          ENDIF
C
          DO IEL = 1, NCEL
            SORTI1(IEL) = GAMAGP * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH)
          ENDDO
C
C                                                    gamma
C CALCUL DE beta EN FONCTION DE P ET rho : beta = rho
C ======================================
C
        ELSEIF(ICCFTH.EQ.162) THEN
C
C     Test des valeurs de rho (on peut éliminer ce test pour optimiser)
          IUTILE = 0
          IF(IUTILE.EQ.1) THEN
            IERR = 0
            DO IEL = 1, NCEL
              IF(RTP(IEL,IRHIPH).LT.0.D0) THEN
                WRITE(NFECRA,4020)RTP(IEL,IRHIPH),IEL
              ENDIF
            ENDDO
            IF(IERR.EQ.1) THEN
              CALL CSEXIT (1)
            ENDIF
          ENDIF
C
          DO IEL = 1, NCEL
            SORTI1(IEL) = RTP(IEL,IRHIPH)**GAMAGP
          ENDDO
C
C
C     CALCUL DE LA CHALEUR MASSIQUE A VOLUME CONSTANT
C     ===============================================
C
C     elle est constante !
C
C                                                           P
C     CALCUL DE L'ENTROPIE EN FONCTION DE P ET rho : s = --------
C     ============================================          gamma
C                                                        rho
        ELSEIF(ICCFTH.EQ.6) THEN
C
C     Test des valeurs de rho (on peut éliminer ce test pour optimiser)
          IERR = 0
          DO IEL = 1, NCEL
            IF(RTP(IEL,IRHIPH).LE.0.D0) THEN
              WRITE(NFECRA,4030)RTP(IEL,IRHIPH),IEL
            ENDIF
          ENDDO
          IF(IERR.EQ.1) THEN
            CALL CSEXIT (1)
          ENDIF
C
          DO IEL = 1, NCEL
            SORTI1(IEL) = RTP(IEL,IPRIPH) / (RTP(IEL,IRHIPH)**GAMAGP)
          ENDDO
C
C
C     CALCUL DE epsilon - Cv.T : epsilon - Cv.T = 0
C     ========================
C
        ELSEIF(ICCFTH.EQ.7) THEN
C
C --- A l'interieur du domaine
C
          DO IEL = 1, NCEL
            SORTI1(IEL) = 0.D0
          ENDDO
C
C --- Sur les bords
C
          DO IFAC = 1, NFABOR
            SORTI2(IFAC) = 0.D0
          ENDDO
C
C
C     CALCUL DES CONDITIONS AUX LIMITES : sur la face IFAC = IFAC0
C     =================================
C
C     PAROI
C     -----
C
        ELSEIF(ICCFTH.EQ.91) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- Calcul du nombre de Mach normal a la frontiere
C
          XMACH =
     &         ( RTP(IEL,IUIPH)*SURFBO(1,IFAC)
     &         + RTP(IEL,IVIPH)*SURFBO(2,IFAC)
     &         + RTP(IEL,IWIPH)*SURFBO(3,IFAC) ) / RA(ISRFBN+IFAC-1)
     &         / SQRT( GAMAGP * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH) )
C
C --- Pression
C
C     On utilise un Neumann en pression. Bien que cela ne permette
C       pas d'utiliser Rusanov, on espere que cela stabilise.
C     On ajoute un test pour eviter de basculer
C       de detente a choc d'un pas de temps à l'autre
C
C     Detente
          IF(XMACH.LT.0.D0.AND.COEFB(IFAC,ICLP).LE.1.D0) THEN
C
            IF(XMACH.GT.2.D0/(1.D0-GAMAGP)) THEN
              COEFB(IFAC,ICLP) = (1.D0 + (GAMAGP-1.D0)/2.D0 * XMACH)
     &             ** (2.D0*GAMAGP/(GAMAGP-1.D0))
C         Si on detend trop fort : Dirichlet nul en pression
C           (la valeur de COEFB ici est utilisee comme indicateur
C            et sera retraitée plus tard dans cfxtcl)
            ELSE
              COEFB(IFAC,ICLP) = RINFIN
            ENDIF
C
C     Choc
          ELSEIF(XMACH.GT.0.D0.AND.COEFB(IFAC,ICLP).GE.1.D0) THEN
C
            COEFB(IFAC,ICLP) = 1.D0 + GAMAGP*XMACH
     &            *( (GAMAGP+1.D0)/4.D0*XMACH
     &                + SQRT(1.D0 + (GAMAGP+1.D0)**2/16.D0*XMACH**2) )
C
C     Oscillation ou Mach nul
          ELSE
            COEFB(IFAC,ICLP) = 1.D0
          ENDIF

C
C     SYMETRIE
C     --------
C
        ELSEIF(ICCFTH.EQ.90) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- On impose un Neumann Homogene en Pression par defaut
C         Pas d'intervention requise
C
C
C --- ENTREE SUBSONIQUE, RHO ET U DONNES (subsonique par hypothèse)
C     ----------------------------------
C
C     Contrairement à la conception initiale, on impose un Dirichlet
C       explicite sur la pression et non pas un flux (on se base cpdt
C       sur la meme valeur physique).
C     Ceci a l'avantage de permettre l'utilisation de Rusanov pour
C       lisser les conditions d'entree.
C     En outre, ceci permet de ne pas avoir à gérer le remplissage de
C       coefb ici (c'est secondaire, car il faut le faire pour la paroi).
C     Si on a des oscillations en temps, on pourra essayer d'eviter
C       le basculement de détente à choc d'un pas de temps à l'autre.
C     La pertinence de cette démarche reste à démontrer.
C
        ELSEIF(ICCFTH.EQ.92) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- Calcul du nombre de Mach entrant normalement a la frontiere
C
          XMACHI =
     &         ( RTP(IEL,IUIPH)*SURFBO(1,IFAC)
     &         + RTP(IEL,IVIPH)*SURFBO(2,IFAC)
     &         + RTP(IEL,IWIPH)*SURFBO(3,IFAC) ) / RA(ISRFBN+IFAC-1)
     &         / SQRT( GAMAGP * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH) )
          XMACHE =
     &         ( COEFA(IFAC,ICLU)*SURFBO(1,IFAC)
     &         + COEFA(IFAC,ICLV)*SURFBO(2,IFAC)
     &         + COEFA(IFAC,ICLW)*SURFBO(3,IFAC) ) /RA(ISRFBN+IFAC-1)
     &         / SQRT( GAMAGP * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH) )
          DXMACH = XMACHI - XMACHE
C
C --- Pression : Detente
          IF(DXMACH.LE.0.D0) THEN
C
            IF(DXMACH.GT.2.D0/(1.D0-GAMAGP)) THEN
              COEFA(IFAC,ICLP) = RTP(IEL,IPRIPH)*
     &             ( (1.D0 + (GAMAGP-1.D0)*0.50D0*DXMACH)
     &               ** (2.D0*GAMAGP/(GAMAGP-1.D0))    )
            ELSEIF(DXMACH.LE.2.D0/(1.D0-GAMAGP) ) THEN
              COEFA(IFAC,ICLP) = 0.D0
            ENDIF
C
C --- Pression : Choc
          ELSE
            COEFA(IFAC,ICLP) = RTP(IEL,IPRIPH)*
     &           (  1.D0 + GAMAGP*DXMACH
     &           *( (GAMAGP+1.D0)*0.25D0*DXMACH
     &           + SQRT(1.D0 + (GAMAGP+1.D0)**2/16.D0*DXMACH**2) )  )
          ENDIF
C
C --- Energie totale
          COEFA(IFAC,ICLP) = RTP(IEL,IPRIPH)
          COEFA(IFAC,ICLE) =
     &         COEFA(IFAC,ICLP)/((GAMAGP-1.D0)*COEFA(IFAC,ICLR))
     &         + 0.5D0*(COEFA(IFAC,ICLU)**2
     &                + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C
C --- ENTREE SUBSONIQUE, RHO*U ET RHO*U*H DONNES (subsonique par hypothèse)
C     ------------------------------------------
C
C     Ceci reste a completer :
        ELSEIF(ICCFTH.EQ.94) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
          WRITE(NFECRA,7000)
C
          CALL CSEXIT (1)
C         ===========
C
C     P   : a calculer avec un newton
C     u et rho en fonction de P
C     E en fonction de H
C     (je l'ai ecrit, reste à le coder ...)
C
C
C     SORTIE SUBSONIQUE (subsonique par hypothèse)
C     -----------------
C
        ELSEIF(ICCFTH.EQ.93) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- Detente :
          IF(COEFA(IFAC,ICLP).LE.RTP(IEL,IPRIPH)) THEN
C
C       Masse volumique
            COEFA(IFAC,ICLR) = RTP(IEL,IRHIPH)
     &           * (COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH))**(1.D0/GAMAGP)
C
C       Vitesse
            COEFA(IFAC,ICLU) = RTP(IEL,IUIPH)
     &           + 2.D0/(GAMAGP-1.D0)
     &           * SQRT(GAMAGP*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH))
     &           * (1.D0-(COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH)
     &                        )**((GAMAGP-1.D0)/(2.D0*GAMAGP)))
     &           * SURFBO(1,IFAC)/RA(ISRFBN+IFAC-1)
C
            COEFA(IFAC,ICLV) = RTP(IEL,IVIPH)
     &           + 2.D0/(GAMAGP-1.D0)
     &           * SQRT( GAMAGP*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH))
     &           * (1.D0-(COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH)
     &                        )**((GAMAGP-1.D0)/(2.D0*GAMAGP)))
     &           * SURFBO(2,IFAC)/RA(ISRFBN+IFAC-1)
C
            COEFA(IFAC,ICLW) = RTP(IEL,IWIPH)
     &           + 2.D0/(GAMAGP-1.D0)
     &           * SQRT( GAMAGP*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH))
     &           * (1.D0-(COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH)
     &                        )**((GAMAGP-1.D0)/(2.D0/GAMAGP)))
     &           * SURFBO(3,IFAC)/RA(ISRFBN+IFAC-1)
C
C       Energie totale
            COEFA(IFAC,ICLE) =
     &           COEFA(IFAC,ICLP)/((GAMAGP-1.D0)*COEFA(IFAC,ICLR))
     &           + 0.5D0*(COEFA(IFAC,ICLU)**2
     &                  + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C --- Choc :
          ELSE
C
C       Masse volumique
            COEFA(IFAC,ICLR) = RTP(IEL,IRHIPH)
     &           * ( (GAMAGP+1.D0)*COEFA(IFAC,ICLP)
     &             + (GAMAGP-1.D0)*RTP(IEL,IPRIPH) )
     &           / ( (GAMAGP-1.D0)*COEFA(IFAC,ICLP)
     &             + (GAMAGP+1.D0)*RTP(IEL,IPRIPH) )
C
C       Vitesse
            COEFA(IFAC,ICLU) = RTP(IEL,IUIPH)
     &           - (COEFA(IFAC,ICLP)-RTP(IEL,IPRIPH))
     &           * SQRT(2.D0/
     &                  (RTP(IEL,IRHIPH)
     &                   *((GAMAGP+1.D0)*COEFA(IFAC,ICLP)
     &                    +(GAMAGP-1.D0)*RTP(IEL,IPRIPH) )))
     &           * SURFBO(1,IFAC)/RA(ISRFBN+IFAC-1)
C
            COEFA(IFAC,ICLV) = RTP(IEL,IVIPH)
     &           - (COEFA(IFAC,ICLP)-RTP(IEL,IPRIPH))
     &           * SQRT(2.D0/
     &                  (RTP(IEL,IRHIPH)
     &                   *((GAMAGP+1.D0)*COEFA(IFAC,ICLP)
     &                    +(GAMAGP-1.D0)*RTP(IEL,IPRIPH) )))
     &           * SURFBO(2,IFAC)/RA(ISRFBN+IFAC-1)
C
            COEFA(IFAC,ICLW) = RTP(IEL,IWIPH)
     &           - (COEFA(IFAC,ICLP)-RTP(IEL,IPRIPH))
     &           * SQRT(2.D0/
     &                  (RTP(IEL,IRHIPH)
     &                   *((GAMAGP+1.D0)*COEFA(IFAC,ICLP)
     &                    +(GAMAGP-1.D0)*RTP(IEL,IPRIPH) )))
     &           * SURFBO(3,IFAC)/RA(ISRFBN+IFAC-1)
C
C       Energie totale
            COEFA(IFAC,ICLE) =
     &           COEFA(IFAC,ICLP)/((GAMAGP-1.D0)*COEFA(IFAC,ICLR))
     &           + 0.5D0*(COEFA(IFAC,ICLU)**2
     &                  + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
          ENDIF
C
C
C     T ET e EN FONCTION DE P ET rho
C     ------------------------------
C      on a donné des valeurs positives strictement (hypothèse)
C
        ELSEIF(ICCFTH.EQ.912.OR.ICCFTH.EQ.60900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Temperature
          COEFA(IFAC,ICLT) =
     &         XMASML*COEFA(IFAC,ICLP)/(RR*COEFA(IFAC,ICLR))
C
C     Energie totale
          COEFA(IFAC,ICLE) =
     &         CV0(IPHAS)*COEFA(IFAC,ICLT)
     &         + 0.5D0*( COEFA(IFAC,ICLU)**2
     &                 + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2 )
C
C
C     rho ET e EN FONCTION DE P ET T
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.913.OR.ICCFTH.EQ.100900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Masse volumique
          COEFA(IFAC,ICLR) =
     &         XMASML*COEFA(IFAC,ICLP)/(RR*COEFA(IFAC,ICLT))
C
C     Energie totale
          COEFA(IFAC,ICLE) =
     &         CV0(IPHAS)*COEFA(IFAC,ICLT)
     &         + 0.5D0*( COEFA(IFAC,ICLU)**2
     &                 + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2 )
C
C
C     rho ET T EN FONCTION DE P ET e
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.914.OR.ICCFTH.EQ.140900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Masse volumique
          COEFA(IFAC,ICLR) = COEFA(IFAC,ICLP)/( (GAMAGP-1.D0)*
     &         (COEFA(IFAC,ICLE)
     &         - 0.5D0*( COEFA(IFAC,ICLU)**2
     &                 + COEFA(IFAC,ICLV)**2
     &                 + COEFA(IFAC,ICLW)**2 ) ) )
C
C     Temperature
          COEFA(IFAC,ICLT)=
     &         XMASML*COEFA(IFAC,ICLP)/(RR*COEFA(IFAC,ICLR))
C
C
C     P ET e EN FONCTION DE rho ET T
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.923.OR.ICCFTH.EQ.150900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Pression
          COEFA(IFAC,ICLP) = COEFA(IFAC,ICLR)*RR/XMASML
     &                                       *COEFA(IFAC,ICLT)
C
C     Energie totale
          COEFA(IFAC,ICLE) = CV0(IPHAS) * COEFA(IFAC,ICLT)
     &         + 0.5D0*( COEFA(IFAC,ICLU)**2
     &                 + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2 )
C
C
C     P ET T EN FONCTION DE rho ET e
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.924.OR.ICCFTH.EQ.210900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Pression
          COEFA(IFAC,ICLP) = (GAMAGP-1.D0)*COEFA(IFAC,ICLR)
     &          *( COEFA(IFAC,ICLE)
     &            - 0.5D0*( COEFA(IFAC,ICLU)**2
     &                    + COEFA(IFAC,ICLV)**2
     &                    + COEFA(IFAC,ICLW)**2 ) )

C
C     Temperature
          COEFA(IFAC,ICLT)=
     &         XMASML*COEFA(IFAC,ICLP)/(RR*COEFA(IFAC,ICLR))
C
C
C --- Fin de test sur GAZ PARFAIT
        ENDIF
C
C
C=======================================================================
C 2.  GAZ PARFAIT AVEC GAMMA VARIABLE
C=======================================================================
C
      ELSEIF(IEOS(IPHAS).EQ.2) THEN
C
C***********************************************************************
C
C     PARAMETRES
C     ==========
C
C-----------------------------------------------------------------------
C
C     EXEMPLES  (A recopier et a adapter apres 'PARAMETRES')
        IF(0.EQ.1) THEN
C-----------------------------------------------------------------------
C
C --> EXEMPLE 1 : GAZ PARFAIT A 3 CONSTITUANTS
C
C     Phase
        IF(IPHAS.EQ.1) THEN
C
C     Masses molaires des constituants (en kg/mol)
          CSTGR(1)  = 18.D-3
          CSTGR(2)  = 32.D-3
          CSTGR(3)  = 28.D-3
C
C     Calcul de la masse molaire du melange au centre des cellules
          IF(ICCFTH.GT.0) THEN
            DO IEL = 1, NCEL
              XMASMR(IEL) = 1.D0 / ( RTP(IEL,ISCA(1))/CSTGR(1)
     &                             + RTP(IEL,ISCA(2))/CSTGR(2)
     &                             + RTP(IEL,ISCA(3))/CSTGR(3) )
            ENDDO

C
C     Calcul du GAMMA equivalent au centre des cellules
            DO IEL = 1, NCEL
              GAMAGR(IEL) = PROPCE(IEL,IPPROC(ICP(IPHAS)))
     &           / ( PROPCE(IEL,IPPROC(ICP(IPHAS))) - RR/XMASMR(IEL) )
            ENDDO
          ENDIF
C
        ENDIF
C
C-----------------------------------------------------------------------
C
C     FIN DES EXEMPLES
        ENDIF
C-----------------------------------------------------------------------
C
C***********************************************************************
C
C
C TEST DE VALIDITE DE GAMAGR : GAMAGR >= 1
C ==========================
C
        IERR = 0
C
        DO IEL = 1, NCEL
          IF(ICCFTH.GT.0 .AND. GAMAGR(IEL).LT.1.D0) THEN
            IERR = 1
            WRITE(NFECRA,1020) IEL, GAMAGR(IEL)
          ENDIF
        ENDDO
C
        IF(IERR.EQ.1) THEN
          CALL CSEXIT (1)
        ENDIF
C
C
C OPTIONS DE CALCUL  : Cp ET Cv VARIABLES
C =================
C
        IF(ICCFTH.EQ.-1) THEN
C
          ICP(IPHAS) = 1
          CP0(IPHAS) = EPZERO
          ICV(IPHAS) = 1
          CV0(IPHAS) = EPZERO
C
C
C INITIALISATIONS PAR DEFAUT :
C ==========================
C
        ELSEIF(ICCFTH.EQ.0) THEN
C
          DO IEL = 1, NCEL
            PROPCE(IEL,IPPROC(ICP(IPHAS))) = CP0(IPHAS)
            PROPCE(IEL,IPPROC(ICV(IPHAS))) =
     &           CP0(IPHAS) - RR/XMASMR(IEL)
            RTP(IEL,IRHIPH) = P0(IPHAS)*XMASMR(IEL)/RR/T0(IPHAS)
            RTP(IEL,IENIPH) =
     &           PROPCE(IEL,IPPROC(ICV(IPHAS)))*T0(IPHAS)
          ENDDO
C
C
C CALCUL DE T ET e EN FONCTION DE P ET rho :
C ========================================
C
        ELSEIF(ICCFTH.EQ.12) THEN
C
          DO IEL = 1, NCEL
C
C     Temperature
            SORTI1(IEL) =
     &           XMASMR(IEL)/RR*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH)
C
C     Energie totale
            SORTI2(IEL) = PROPCE(IEL,IPPROC(ICV(IPHAS)))*SORTI1(IEL)
     &    + 0.5D0*( RTP(IEL,IUIPH)**2
     &           + RTP(IEL,IVIPH)**2 + RTP(IEL,IWIPH)**2 )
C
          ENDDO
C
C
          IF(IMODIF.GT.0) THEN
C       Affectation directe a RTP
            DO IEL = 1, NCEL
              RTP(IEL,ITKIPH) = SORTI1(IEL)
              RTP(IEL,IENIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C CALCUL DE rho ET e EN FONCTION DE P ET T :
C ========================================
C
        ELSEIF(ICCFTH.EQ.13) THEN
C
          DO IEL = 1, NCEL
C
C     Masse volumique
            SORTI1(IEL) =
     &           XMASMR(IEL)/RR*RTP(IEL,IPRIPH)/RTP(IEL,ITKIPH)
C
C     Energie totale
            SORTI2(IEL) =
     &           PROPCE(IEL,IPPROC(ICV(IPHAS)))*RTP(IEL,ITKIPH)
     &    + 0.5D0*( RTP(IEL,IUIPH)**2
     &           + RTP(IEL,IVIPH)**2 + RTP(IEL,IWIPH)**2 )
C
          ENDDO
C
C
          IF(IMODIF.GT.0) THEN
C       Affectation directe a RTP
            DO IEL = 1, NCEL
              RTP(IEL,IRHIPH) = SORTI1(IEL)
              RTP(IEL,IENIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C CALCUL DE rho ET T EN FONCTION DE P ET e :
C ========================================
C
        ELSEIF(ICCFTH.EQ.14) THEN
C
          DO IEL = 1, NCEL
C
C     Masse volumique
            SORTI1(IEL) =
     &           RTP(IEL,IPRIPH)/(GAMAGR(IEL)-1.D0)/( RTP(IEL,IENIPH)
     &  - 0.5D0*( RTP(IEL,IUIPH)**2
     &           + RTP(IEL,IVIPH)**2 + RTP(IEL,IWIPH)**2 ) )
C
C     Temperature
            SORTI2(IEL) = XMASMR(IEL)/RR*RTP(IEL,IPRIPH)/SORTI1(IEL)
C
          ENDDO
C
C
          IF(IMODIF.GT.0) THEN
C       Affectation directe a RTP
            DO IEL = 1, NCEL
              RTP(IEL,IRHIPH) = SORTI1(IEL)
              RTP(IEL,ITKIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C CALCUL DE P ET e EN FONCTION DE rho ET T :
C ========================================
C
        ELSEIF(ICCFTH.EQ.23) THEN
C
          DO IEL = 1, NCEL
C
C     Pression
            SORTI1(IEL) =
     &           RTP(IEL,IRHIPH)*RR/XMASMR(IEL)*RTP(IEL,ITKIPH)
C
C     Energie totale
            SORTI2(IEL) =
     &           PROPCE(IEL,IPPROC(ICV(IPHAS)))*RTP(IEL,ITKIPH)
     &    + 0.5D0*( RTP(IEL,IUIPH)**2
     &           + RTP(IEL,IVIPH)**2 + RTP(IEL,IWIPH)**2 )
C
          ENDDO
C
C
          IF(IMODIF.GT.0) THEN
C       Affectation directe a RTP
            DO IEL = 1, NCEL
              RTP(IEL,IPRIPH) = SORTI1(IEL)
              RTP(IEL,IENIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C
C CALCUL DE P ET T EN FONCTION DE rho ET e :
C ========================================
C
        ELSEIF(ICCFTH.EQ.24) THEN
C
          DO IEL = 1, NCEL
C
C     Pression
            SORTI1(IEL) =
     &           (GAMAGR(IEL)-1.D0)*RTP(IEL,IRHIPH)*( RTP(IEL,IENIPH)
     &  - 0.5D0*( RTP(IEL,IUIPH)**2
     &           + RTP(IEL,IVIPH)**2 + RTP(IEL,IWIPH)**2 ) )
C
C     Temperature
            SORTI2(IEL) = XMASMR(IEL)/RR*SORTI1(IEL)/RTP(IEL,IRHIPH)
C
          ENDDO
C
          IF(IMODIF.GT.0) THEN
C       Affectation directe a RTP
            DO IEL = 1, NCEL
              RTP(IEL,IPRIPH) = SORTI1(IEL)
              RTP(IEL,ITKIPH) = SORTI2(IEL)
            ENDDO
          ENDIF
C
C            2                            2          P
C CALCUL DE c  EN FONCTION DE P ET rho : c  = gamma*---
C ====================================              rho
C
        ELSEIF(ICCFTH.EQ.126) THEN
C
          DO IEL = 1, NCEL
C
C --- Test de positivite de P
            IF(RTP(IEL,IPRIPH).LT.0.D0) THEN
              WRITE(NFECRA,1110) IEL , RTP(IEL,IPRIPH)
              IERR = 1
C
C --- Test de positivite de rho
            ELSEIF(RTP(IEL,IRHIPH).LE.0.D0) THEN
              WRITE(NFECRA,1120) IEL , RTP(IEL,IRHIPH)
              IERR = 1
C
            ELSE
C
C --- Calcul
              SORTI1(IEL) =
     &             GAMAGR(IEL) * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH)
C
            ENDIF
C
          ENDDO
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C                                                    gamma
C CALCUL DE beta EN FONCTION DE P ET rho : beta = rho
C ======================================
C
        ELSEIF(ICCFTH.EQ.162) THEN
C
          DO IEL = 1, NCEL
C
C --- Test de positivite de rho
            IF(RTP(IEL,IRHIPH).LT.0.D0) THEN
              WRITE(NFECRA,1220) IEL , RTP(IEL,IRHIPH)
              IERR = 1
C
            ELSE
C
C --- Calcul
              SORTI1(IEL) = RTP(IEL,IRHIPH)**GAMAGR(IEL)
C
            ENDIF
C
          ENDDO
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C                           R
C CALCUL DE Cv : Cv = Cp - ---
C ============              M
C
        ELSEIF(ICCFTH.EQ.432) THEN
C
          DO IEL = 1, NCEL
C
            SORTI1(IEL) = PROPCE(IEL,IPPROC(ICP(IPHAS)))-RR/XMASMR(IEL)
C
          ENDDO
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C                                                       P
C CALCUL DE L'ENTROPIE EN FONCTION DE P ET rho : s = --------
C ============================================          gamma
C                                                    rho
        ELSEIF(ICCFTH.EQ.6) THEN
C
          DO IEL = 1, NCEL
C
C --- Test de positivite de P
            IF(RTP(IEL,IPRIPH).LT.0.D0) THEN
              WRITE(NFECRA,1310) IEL , RTP(IEL,IPRIPH)
              IERR = 1
C
C --- Test de positivite de rho
            ELSEIF(RTP(IEL,IRHIPH).LE.0.D0) THEN
              WRITE(NFECRA,1320) IEL , RTP(IEL,IRHIPH)
              IERR = 1
C
            ELSE
C
C --- Calcul
              SORTI1(IEL) =
     &             RTP(IEL,IPRIPH) / (RTP(IEL,IRHIPH)**GAMAGR(IEL))
C
            ENDIF
C
          ENDDO
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C
C CALCUL DE epsilon - Cv.T : epsilon - Cv.T = 0
C ========================
C
        ELSEIF(ICCFTH.EQ.7) THEN
C
C --- A l'interieur du domaine
C
          DO IEL = 1, NCEL
C
            SORTI1(IEL) = 0.D0
C
          ENDDO
C
C --- Sur les bords
C
          DO IFAC = 1, NFABOR
C
            SORTI2(IFAC) = 0.D0
C
          ENDDO
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C
C CALCUL DES CONDITIONS AUX LIMITES :
C =================================
C
C     PAROI/SYMETRIE
C     --------------
        ELSEIF(ICCFTH.EQ.91) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- Calcul du nombre de Mach normal a la frontiere
C
          XMACH = ( RTP(IEL,IUIPH)*SURFBO(1,IFAC)
     &           + RTP(IEL,IVIPH)*SURFBO(2,IFAC)
     &           + RTP(IEL,IWIPH)*SURFBO(3,IFAC) ) / RA(ISRFBN+IFAC-1)
     &         / SQRT( GAMAGR(IEL)*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH) )
C
C --- Pression / Entropie : Detente
C
          COEFA(IFAC,ICLP) = 0.D0
C
          IF(XMACH.LE.0.D0 .AND. XMACH.GT.2.D0/(1.D0-GAMAGR(IEL))) THEN
            COEFB(IFAC,ICLP) = (1.D0 + (GAMAGR(IEL)-1.D0)/2.D0 * XMACH)
     &           ** (2.D0*GAMAGR(IEL)/(GAMAGR(IEL)-1.D0))
            COEFB(IFAC,ICLT) = 1.D0
C
          ELSEIF(XMACH.LE.2.D0/(1.D0-GAMAGR(IEL)) ) THEN
            COEFB(IFAC,ICLP) = 0.D0
            COEFB(IFAC,ICLT) = 1.D0
C
C --- Pression / Entropie : Choc
C
          ELSE
            COEFB(IFAC,ICLP) = 1.D0 + GAMAGR(IEL)*XMACH
     &            *( (GAMAGR(IEL)+1.D0)/4.D0*XMACH
     &           + SQRT(1.D0 + (GAMAGR(IEL)+1.D0)**2/16.D0*XMACH**2) )
            COEFB(IFAC,ICLT) = COEFB(IFAC,ICLP)/(1.D0-COEFB(IFAC,ICLP))
     &          / RTP(IEL,IPRIPH) * ( RTP(IEL,IRHIPH)
     &              * (RTP(IEL,IUIPH)**2
     &                +RTP(IEL,IVIPH)**2+RTP(IEL,IWIPH)**2)
     &              + RTP(IEL,IPRIPH) *(1.D0-COEFB(IFAC,ICLP)) )
          ENDIF
C
C --- Energie totale : Epsilon sup
C
          COEFA(IFAC,ICLE) = 0.D0
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C     ENTREE
C     ------
        ELSEIF(ICCFTH.EQ.92) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- Calcul du nombre de Mach entrant normalement a la frontiere
C
          XMACHI = ( RTP(IEL,IUIPH)*SURFBO(1,IFAC)
     &         + RTP(IEL,IVIPH)*SURFBO(2,IFAC)
     &         + RTP(IEL,IWIPH)*SURFBO(3,IFAC) )/RA(ISRFBN+IFAC-1)
     &         / SQRT(GAMAGR(IEL)*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH))
          XMACHE = ( COEFA(IFAC,ICLU)*SURFBO(1,IFAC)
     &         + COEFA(IFAC,ICLV)*SURFBO(2,IFAC)
     &         + COEFA(IFAC,ICLW)*SURFBO(3,IFAC) )/RA(ISRFBN+IFAC-1)
     &         / SQRT(GAMAGR(IEL)*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH))
          DXMACH = XMACHI - XMACHE
C
C --- Pression : Detente
          IF(DXMACH.LE.0.D0) THEN
C
            IF(DXMACH.GT.2.D0/(1.D0-GAMAGR(IEL))) THEN
              COEFA(IFAC,ICLP) = RTP(IEL,IPRIPH)*
     &             ( (1.D0 + (GAMAGR(IEL)-1.D0)*0.50D0*DXMACH)
     &               ** (2.D0*GAMAGR(IEL)/(GAMAGR(IEL)-1.D0))  )
            ELSEIF(DXMACH.LE.2.D0/(1.D0-GAMAGR(IEL)) ) THEN
              COEFA(IFAC,ICLP) = 0.D0
            ENDIF
C
C --- Pression : Choc
          ELSE
            COEFA(IFAC,ICLP) = RTP(IEL,IPRIPH)*
     &           (  1.D0 + GAMAGR(IEL)*DXMACH
     &           *( (GAMAGR(IEL)+1.D0)*0.25D0*DXMACH
     &           + SQRT(1.D0 + (GAMAGR(IEL)+1.D0)**2/16.D0
     &                                           *DXMACH**2) )  )
          ENDIF
C
C --- Energie totale
C
          COEFA(IFAC,ICLP) = RTP(IEL,IPRIPH)
C
          COEFA(IFAC,ICLE) =
     &         COEFA(IFAC,ICLP)/((GAMAGR(IEL)-1.D0)*COEFA(IFAC,ICLR))
     &         + 0.5D0*(COEFA(IFAC,ICLU)**2
     &                + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C     SORTIE
C     ------
        ELSEIF(ICCFTH.EQ.93) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C --- Calcul du nombre de Mach sortant normalement a la frontiere
C
          XMACH = ( RTP(IEL,IUIPH)*SURFBO(1,IFAC)
     &           + RTP(IEL,IVIPH)*SURFBO(2,IFAC)
     &           + RTP(IEL,IWIPH)*SURFBO(3,IFAC) ) / RA(ISRFBN+IFAC-1)
     &         / SQRT(GAMAGR(IEL)*RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH))
C
C --- Sortie supersonique : Dirichlet sur toutes les variables
C     -------------------
          IF(XMACH.GE.1.D0) THEN
            DO IVAR = 1, NVAR
              COEFA(IFAC,ICLRTP(IVAR,ICOEF)) = RTP(IEL,IVAR)
            ENDDO
C
C           Entropie
            COEFA(IFAC,ICLT) =
     &           RTP(IEL,IPRIPH)/RTP(IEL,IRHIPH)**GAMAGR(IEL)
C
C --- Sortie subsonique
C     -----------------
          ELSEIF(XMACH.LT.1.D0 .AND. XMACH.GE.0.D0) THEN
C
C --- Detente :
            IF(COEFA(IFAC,ICLP).LE.RTP(IEL,IPRIPH)) THEN
C
C       Masse volumique
              COEFA(IFAC,ICLR) = RTP(IEL,IRHIPH)
     &             * (COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH))
     &                **(1.D0/GAMAGR(IEL))
C
C       Vitesse
              COEFA(IFAC,ICLU) = RTP(IEL,IUIPH)
     &             + 2.D0/(GAMAGR(IEL)-1.D0)
     & * SQRT( GAMAGR(IEL) * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH) )
     &             * ( 1.D0
     & - (COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH))
     &               **((GAMAGR(IEL)-1.D0)/2.D0/GAMAGR(IEL)) )
     & * SURFBO(1,IFAC) / RA(ISRFBN+IFAC-1)
C
              COEFA(IFAC,ICLV) = RTP(IEL,IVIPH)
     &             + 2.D0/(GAMAGR(IEL)-1.D0)
     & * SQRT( GAMAGR(IEL) * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH) )
     &             * ( 1.D0
     & - (COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH))
     &               **((GAMAGR(IEL)-1.D0)/2.D0/GAMAGR(IEL)) )
     & * SURFBO(2,IFAC) / RA(ISRFBN+IFAC-1)
C
              COEFA(IFAC,ICLW) = RTP(IEL,IWIPH)
     &             + 2.D0/(GAMAGR(IEL)-1.D0)
     & * SQRT( GAMAGR(IEL) * RTP(IEL,IPRIPH) / RTP(IEL,IRHIPH) )
     &             * ( 1.D0
     & - (COEFA(IFAC,ICLP)/RTP(IEL,IPRIPH))
     &               **((GAMAGR(IEL)-1.D0)/2.D0/GAMAGR(IEL)) )
     & * SURFBO(3,IFAC) / RA(ISRFBN+IFAC-1)
C
C       Energie totale
              COEFA(IFAC,ICLE) = COEFA(IFAC,ICLP)
     & /( (GAMAGR(IEL)-1.D0)*COEFA(IFAC,ICLR) )
     &              + 0.5D0*(COEFA(IFAC,ICLU)**2
     &                     + COEFA(IFAC,ICLV)**2
     &                     + COEFA(IFAC,ICLW)**2)
C
C       Entropie
              COEFA(IFAC,ICLT) = COEFA(IFAC,ICLP)
     &                             /COEFA(IFAC,ICLR)**GAMAGR(IEL)
C
C --- Choc :
            ELSE
C
C       Masse volumique
              COEFA(IFAC,ICLR) = RTP(IEL,IRHIPH)
     & * ( (GAMAGR(IEL)+1.D0)*COEFA(IFAC,ICLP)
     &   + (GAMAGR(IEL)-1.D0)*RTP(IEL,IPRIPH) )
     & / ( (GAMAGR(IEL)-1.D0)*COEFA(IFAC,ICLP)
     &   + (GAMAGR(IEL)+1.D0)*RTP(IEL,IPRIPH) )
C
C       Vitesse
              COEFA(IFAC,ICLU) = RTP(IEL,IUIPH)
     & - (COEFA(IFAC,ICLP)-RTP(IEL,IPRIPH))*SQRT(2.D0/RTP(IEL,IRHIPH)
     & / ( (GAMAGR(IEL)+1.D0)*COEFA(IFAC,ICLP)
     &   + (GAMAGR(IEL)-1.D0)*RTP(IEL,IPRIPH) ))
     & * SURFBO(1,IFAC) / RA(ISRFBN+IFAC-1)
C
              COEFA(IFAC,ICLV) = RTP(IEL,IVIPH)
     & - (COEFA(IFAC,ICLP)-RTP(IEL,IPRIPH))*SQRT(2.D0/RTP(IEL,IRHIPH)
     & / ( (GAMAGR(IEL)+1.D0)*COEFA(IFAC,ICLP)
     &   + (GAMAGR(IEL)-1.D0)*RTP(IEL,IPRIPH) ))
     & * SURFBO(2,IFAC) / RA(ISRFBN+IFAC-1)
C
              COEFA(IFAC,ICLW) = RTP(IEL,IWIPH)
     & - (COEFA(IFAC,ICLP)-RTP(IEL,IPRIPH))*SQRT(2.D0/RTP(IEL,IRHIPH)
     & / ( (GAMAGR(IEL)+1.D0)*COEFA(IFAC,ICLP)
     &   + (GAMAGR(IEL)-1.D0)*RTP(IEL,IPRIPH) ))
     & * SURFBO(3,IFAC) / RA(ISRFBN+IFAC-1)
C
C       Energie totale
              COEFA(IFAC,ICLE) = COEFA(IFAC,ICLP)
     & /( (GAMAGR(IEL)-1.D0)*COEFA(IFAC,ICLR) )
     &     + 0.5D0*(COEFA(IFAC,ICLU)**2
     &            + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C       Entropie
              COEFA(IFAC,ICLT) = COEFA(IFAC,ICLP)
     &                             /COEFA(IFAC,ICLR)**GAMAGR(IEL)
C
            ENDIF
C
          ELSE
            WRITE(NFECRA,*) 'ICCFTH = ',ICCFTH,'  MACH = ',XMACH
            IERR = 1
          ENDIF
C
          IF(IERR.EQ.1) CALL CSEXIT (1)
C
C
C     T ET e EN FONCTION DE P ET rho
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.912.OR.ICCFTH.EQ.60900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Temperature
          COEFA(IFAC,ICLT) = XMASMR(IEL)/RR*COEFA(IFAC,ICLP)
     &                                        /COEFA(IFAC,ICLR)
C
C     Energie totale
          COEFA(IFAC,ICLE) = PROPCE(IEL,IPPROC(ICV(IPHAS)))
     &               * COEFA(IFAC,ICLT) + 0.5D0*( COEFA(IFAC,ICLU)**2
     &                    + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C
C     rho ET e EN FONCTION DE P ET T
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.913.OR.ICCFTH.EQ.100900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Masse volumique
          COEFA(IFAC,ICLR) = XMASMR(IEL)/RR*COEFA(IFAC,ICLP)
     &                                       /COEFA(IFAC,ICLT)
C
C     Energie totale
          COEFA(IFAC,ICLE) = PROPCE(IEL,IPPROC(ICV(IPHAS)))
     &               * COEFA(IFAC,ICLT) + 0.5D0*( COEFA(IFAC,ICLU)**2
     &                    + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C
C     rho ET T EN FONCTION DE P ET e
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.914.OR.ICCFTH.EQ.140900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Masse volumique
          COEFA(IFAC,ICLR) = COEFA(IFAC,ICLP)/(GAMAGR(IEL)-1.D0)
     &           / (COEFA(IFAC,ICLE) - 0.5D0*( COEFA(IFAC,ICLU)**2
     &                 + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2 ))
C
C     Temperature
          COEFA(IFAC,ICLT)= XMASMR(IEL)/RR*COEFA(IFAC,ICLP)
     &                                       /COEFA(IFAC,ICLR)
C
C
C     P ET e EN FONCTION DE rho ET T
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.923.OR.ICCFTH.EQ.150900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Pression
          COEFA(IFAC,ICLP) = COEFA(IFAC,ICLR)*RR/XMASMR(IEL)
     &                                       *COEFA(IFAC,ICLT)
C
C     Energie totale
          COEFA(IFAC,ICLE) = PROPCE(IEL,IPPROC(ICV(IPHAS)))
     &               * COEFA(IFAC,ICLT) + 0.5D0*( COEFA(IFAC,ICLU)**2
     &                    + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2)
C
C
C     P ET T EN FONCTION DE rho ET e
C     ------------------------------
C
        ELSEIF(ICCFTH.EQ.924.OR.ICCFTH.EQ.210900) THEN
C
          IFAC = IFAC0
          IEL  = IFABOR(IFAC)
C
C     Pression
          COEFA(IFAC,ICLP) = (GAMAGR(IEL)-1.D0)*COEFA(IFAC,ICLR)
     &          *( COEFA(IFAC,ICLE) - 0.5D0*( COEFA(IFAC,ICLU)**2
     &                + COEFA(IFAC,ICLV)**2 + COEFA(IFAC,ICLW)**2 ) )

C
C     Temperature
          COEFA(IFAC,ICLT)= XMASMR(IEL)/RR*COEFA(IFAC,ICLP)
     &                                       /COEFA(IFAC,ICLR)
C
C
C --- Fin de test sur GAZ PARFAIT A GAMMA VARIABLE
        ENDIF
C
C --- Fin de test sur les thermos
      ENDIF
C
C
C--------
C FORMATS
C--------
C
 1010 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********                                               ',/,
     &'@                                                            ',/,
     &'@    GAMMA DOIT ETRE UN REEL SUPERIEUR OU EGAL A 1           ',/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 1020 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********                                               ',/,
     &'@                                                            ',/,
     &'@    GAMMA DOIT ETRE UN REEL SUPERIEUR OU EGAL A 1           ',/,
     &'@    IL EST NEGATIF OU NUL DANS LA CELLULE ',I10              ,/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
 2010 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********   CALCUL DE LA MASSE VOLUMIQUE                ',/,
     &'@                                                            ',/,
     &'@    LA TEMPERATURE DOIT ETRE UN REEL POSITIF STRICTEMENT    ',/,
     &'@    ELLE VAUT ',E12.4   ,' DANS LA CELLULE ',I10             ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 3010 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********   CALCUL DE LA TEMPERATURE                    ',/,
     &'@                                                            ',/,
     &'@    LA MASSE VOLUMIQUE DOIT ETRE UN REEL POSITIF STRICTEMENT',/,
     &'@    ELLE VAUT ',E12.4   ,' DANS LA CELLULE ',I10             ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
 4010 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMODYNAMIQUE ''C2''',/,
     &'@                                                            ',/,
     &'@    LA MASSE VOLUMIQUE DOIT ETRE UN REEL POSITIF STRICTEMENT',/,
     &'@    ELLE VAUT ',E12.4   ,' DANS LA CELLULE ',I10             ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 4020 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMO        ''BETA''',/,
     &'@                                                            ',/,
     &'@    LA MASSE VOLUMIQUE DOIT ETRE UN REEL POSITIF            ',/,
     &'@    ELLE VAUT ',E12.4   ,' DANS LA CELLULE ',I10             ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 4030 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USCFTH                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMODYNAMIQUE ''S'' ',/,
     &'@                                                            ',/,
     &'@    LA MASSE VOLUMIQUE DOIT ETRE UN REEL POSITIF STRICTEMENT',/,
     &'@    ELLE VAUT ',E12.4   ,' DANS LA CELLULE ',I10             ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
 7000 FORMAT (
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : CONDITION A LA LIMITE NON DISPONIBLE USCFTH ',/,
     &'@    *********                                               ',/,
     &'@                                                            ',/,
     &'@    La condition a la limite de type debit et debit         ',/,
     &'@      enthalpique imposes n''est pas disponible dans la     ',/,
     &'@      version courante.                                     ',/,
     &'@                                                            ',/,
     &'@    Le calcul ne peut pas etre execute.                     ',/,
     &'@                                                            ',/,
     &'@    Modifier uscfcl.                                        ',/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 8000 FORMAT (
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : MASSE VOLUMIQUE NEGATIVE OU NULLE (USCFTH)  ',/,
     &'@    *********                                               ',/,
     &'@                                                            ',/,
     &'@    rho est negatif ou nul en ',I10   ,' cellules           ',/,
     &'@    On le limite a ',E12.4                                   ,/,
     &'@    et on met fin au calcul.                                ',/,
     &'@                                                            ',/,
     &'@    Si on souhaite que le calcul se poursuive, on peut      ',/,
     &'@      forcer un clipping standard par valeur inferieure en  ',/,
     &'@      renseignant la variable SCAMIN associee au scalaire   ',/,
     &'@      representant la masse volumique : ISCA(IRHO(IPHAS))   ',/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 8100 FORMAT (
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : ENERGIE INTERNE NEGATIVE OU NULLE (USCFTH)  ',/,
     &'@    *********                                               ',/,
     &'@                                                            ',/,
     &'@    e-0.5*u*u est negatif ou nul en ',I10   ,' cellules     ',/,
     &'@    On le limite a ',E12.4                                   ,/,
     &'@    et on met fin au calcul.                                ',/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
C
C
C formats a eliminer apres mise au point de gama variable
C
C

 1110 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USTHMO                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMODYNAMIQUE ''C2''',/,
     &'@                                                            ',/,
     &'@    P DOIT ETRE UN REEL POSITIF                             ',/,
     &'@    IL EST NEGATIF DANS LA CELLULE ',I10                     ,/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
 1120 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USTHMO                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMODYNAMIQUE ''C2''',/,
     &'@                                                            ',/,
     &'@    RHO DOIT ETRE UN REEL POSITIF STRICTEMENT               ',/,
     &'@    IL EST NEGATIF OU NUL DANS LA CELLULE ',I10              ,/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
 1220 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USTHMO                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMO ''BETA''       ',/,
     &'@                                                            ',/,
     &'@    RHO DOIT ETRE UN REEL POSITIF                           ',/,
     &'@    IL EST NEGATIF DANS LA CELLULE ',I10                     ,/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
 1310 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USTHMO                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMODYNAMIQUE ''S'' ',/,
     &'@                                                            ',/,
     &'@    P DOIT ETRE UN REEL POSITIF                             ',/,
     &'@    IL EST NEGATIF DANS LA CELLULE ',I10                     ,/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)
C
 1320 FORMAT(
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/,
     &'@ @@ ATTENTION : PROBLEME DANS USTHMO                        ',/,
     &'@    *********   CALCUL DE LA VARIABLE THERMODYNAMIQUE ''S'' ',/,
     &'@                                                            ',/,
     &'@    RHO DOIT ETRE UN REEL POSITIF STRICTEMENT               ',/,
     &'@    IL EST NEGATIF OU NUL DANS LA CELLULE ',I10              ,/,
     &'@    IL A POUR VALEUR ',E12.4                                 ,/,
     &'@                                                            ',/,
     &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,
     &'@                                                            ',/)

C
C----
C FIN
C----
C
      RETURN
C
      END
c@z
