!-------------------------------------------------------------------------------

!                      Code_Saturne version 2.0.4n!                      --------------------------

!     This file is part of the Code_Saturne Kernel, element of the
!     Code_Saturne CFD tool.

!     Copyright (C) 1998-2009 EDF S.A., France

!     contact: saturne-support@edf.fr

!     The Code_Saturne Kernel is free software; you can redistribute it
!     and/or modify it under the terms of the GNU General Public License
!     as published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.

!     The Code_Saturne Kernel is distributed in the hope that it will be
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with the Code_Saturne Kernel; if not, write to the
!     Free Software Foundation, Inc.,
!     51 Franklin St, Fifth Floor,
!     Boston, MA  02110-1301  USA

!-------------------------------------------------------------------------------

subroutine uscfth &
!================

 ( 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     )

!===============================================================================
! Purpose:
! -------

!    User subroutine.

!    Define thermodynamic laws (especially for the compressible flow scheme).

!    This user subroutine is mandatory for the compressible flow scheme.


! Introduction
! ============

! This user subroutine allows to define all physical properties and
! variables, through the implementation of thermodynamic laws.


! Avalable thermodynamic laws
! ===========================

!  1. Perfect gas (the molar mass 'xmasml' must be provided)
!  2. Perfect gas with non constant Gamma (example to be adapted)
!  3. Van Der Waals (not yet implemented)


! Implemented calculations
! ========================

! This user subroutine implements the computation of several quantities.
! Each calculation has to be explicitly implemented in the appropriate
! section below (already done for perfect gas).


! Selection of the quantity to return
! ===================================

! When calling the user subroutine, the integer 'iccfth' specifies which
! calculation has to be performed (and which quantity has to be returned).
! The values for 'iccfth' for each case are provided below.
! For some configurations, two systems of references are used for 'iccfth'
! (this is useful to make tests easier to implement in the calling
! subroutines): both systems are explained hereafter for information.

! First system:

!   the variables are referred to using an index i:
!     Variable  P  rho  T   e   h   s  'internal energy - CvT'
!        Index  1   2   3   4   5   6              7

!   iccfth is as follows, depending on which quantity needs to be computed:
!     - compute all variables at cell centers from variable i
!                                              and variable j (i<j):
!               => iccfth = 10*i+j
!     - compute all variables at boundary faces from variable i
!                                                and variable j (i<j):
!               => iccfth = 10*i+j+900

! Second system:

!   the variables are referred to using a different index i:
!     Variable  P  rho  T  e  s
!        Index  2   3   5  7  13

!   iccfth is as follows, depending on which quantity needs to be computed:
!     - compute all variables at cell centers from variable i
!                                              and variable j (i<j):
!               => iccfth = i*j*10000
!     - compute all variables at boundary faces from variable i
!                                                and variable j (i<j):
!               => iccfth = i*j*10000+900

! Other quantities:

!   the variables are referred to using the index of the first system.
!   iccfth is defined as follows:
!     - compute variable i at cell centers (for s and 'internal energy-CvT')
!               => iccfth = i
!                                   \partial(variable i)|
!     - compute partial derivative  --------------------|
!                                   \partial(variable j)|variable k
!               => iccfth = 100*i+10*j+k
!     - compute boundary conditions, resp. symmetry, wall, inlet, outlet:
!               => iccfth = 91, 92, 93, 94


! Values of iccfth
! ================

! To summarize, the values for iccfth are as follows:

!   Values at the cell centers:

!   -> set calculation options (cst/variable cp)   : iccfth = -1
!   -> set default initialization                  : iccfth =  0
!   -> calculate gamma                             : iccfth =  1
!   -> verification of the density                 : iccfth = -2
!   -> verification of the energy                  : iccfth = -4
!   -> calculation of temperature and energy
!                     from pressure and density    : iccfth =  12 or  60000
!   -> calculation of density and energy
!                     from pressure and temperature: iccfth =  13 or 100000
!   -> calculation of density and temperature
!                     from pressure and energy     : iccfth =  14 or 140000
!   -> calculation of pressure and energy
!                     from density and temperature : iccfth =  23 or 150000
!   -> calculation of pressure and temperature
!                     from density and energy      : iccfth =  24 or 210000
!
!                      2    dP |
!   -> calculation of c  = ----|                   : iccfth = 126
!                          drho|s
!
!                            dP|
!   -> calculation of beta = --|                   : iccfth = 162
!                            ds|rho
!
!                          de|
!   -> calculation of Cv = --|                     : iccfth = 432
!                          dT|rho
!
!   -> calculation of entropie                     : iccfth =   6
!
!
!   Values at the boundary faces
!
!   -> calculation of the boundary conditions:
!     - symmetry                                   : iccfth =  90
!     - wall                                       : iccfth =  91
!     - inlet                                      : iccfth =  92
!     - outlet                                     : iccfth =  93
!     - different outlet,not implemented yet       : iccfth =  94
!
!   -> calculation of the variables at the faces for boundary conditions:
!     - temperature and energy
!         from pressure and density                : iccfth = 912 ou  60900
!     - density and energy
!         from pressure and temperature            : iccfth = 913 ou 100900
!     - density and temperature
!         from pressure and energy                 : iccfth = 914 ou 140900
!     - pressure and energy
!         from density and temperature             : iccfth = 923 ou 150900
!     - pressure and temperature
!         from density and energy                  : iccfth = 924 ou 210900


!   Values at the cell centers and at the boundary faces

!   -> calculation of 'internal energy - Cv.T'     : iccfth =   7

!-------------------------------------------------------------------------------
! Arguments
!__________________.____._____.________________________________________________.
!    nom           !type!mode !                   role                         !
!__________________!____!_____!________________________________________________!
! idbia0           ! i  ! <-- ! number of first free position in ia            !
! idbra0           ! i  ! <-- ! number of first free position in ra            !
! ndim             ! i  ! <-- ! spatial dimension                              !
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
! ncel             ! i  ! <-- ! number of cells                                !
! nfac             ! i  ! <-- ! number of interior faces                       !
! nfabor           ! i  ! <-- ! number of boundary faces                       !
! nfml             ! i  ! <-- ! number of families (group classes)             !
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
! nnod             ! i  ! <-- ! number of vertices                             !
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
! nvar             ! i  ! <-- ! total number of variables                      !
! nscal            ! i  ! <-- ! total number of scalars                        !
! nphas            ! i  ! <-- ! number of phases                               !
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
! iprfml           ! ia ! <-- ! property numbers per family                    !
!  (nfml, nprfml)  !    !     !                                                !
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
! nodfac(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary developpement !
! ituser(nituse    ! ia ! <-> ! user-reserved integer work array               !
! ia(*)            ! ia ! --- ! main integer work array                        !
! xyzcen           ! ra ! <-- ! cell centers                                   !
!  (ndim, ncelet)  !    !     !                                                !
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
!  (ndim, nfac)    !    !     !                                                !
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
!  (ndim, nfavor)  !    !     !                                                !
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
!  (ndim, nfac)    !    !     !                                                !
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
!  (ndim, nfabor)  !    !     !                                                !
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
!  (ndim, nnod)    !    !     !                                                !
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
!  (ncelet, *)     !    !     !  (at current and preceding time steps)         !
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
!  (nfabor, *)     !    !     !                                                !
! sorti1,2(*)      ! ra ! --> ! output variable (unused if iccfth.lt.0)        !
! gamagr(*)        ! ra ! --> ! equivalent "gamma" constant of the gas         !
!                  !    !     !   (unused if iccfth.lt.0)                      !
!                  !    !     !   (first value only used for perfect gas)      !
! xmasmr(*)        ! ra ! --> ! molar mass of the components of the gas        !
!                  !    !     !   (unused if iccfth.lt.0)                      !
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary developpement    !
! rtuser(nituse    ! ra ! <-> ! user-reserved real work array                  !
! ra(*)            ! ra ! --- ! main real work array                           !
!__________________!____!_____!________________________________________________!

!     Type: i (integer), r (real), s (string), a (array), l (logical),
!           and composite types (ex: ra real array)
!     mode: <-- input, --> output, <-> modifies data, --- work array
!===============================================================================

implicit none

!===============================================================================
! Common blocks
!===============================================================================

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"

!===============================================================================

! Arguments

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

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(*)

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,*)

double precision sorti1(*), sorti2(*), gamagr(*), xmasmr(*)
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)

! Local variables

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

integer          npmax
parameter (npmax = 1000)
double precision cstgr(npmax)

!===============================================================================


!===============================================================================
! 0. Initialization.
!    No user input required.
!===============================================================================

! Memory pointers
idebia = idbia0
idebra = idbra0

! Error indicator (stop if non zero)
ierr   = 0

! Rank of the variables in their associated arrays
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

! For calculation of values at the cell centers,
!   ifac0 > indicates that the array rtp must be modified
! For calculation of values at the cell faces,
!   ifac0 is the number of the current face
ifac0 = imodif

!===============================================================================
! 1. Thermodynamic law choice (for each phase)
!    User input required.
!===============================================================================

! --> ieos = 1: Perfect gas with constant Gamma
! --> ieos = 2: Perfect gas with variable Gamma
! --> ieos = 3: Van Der Waals

do iiph = 1, nphas
  ieos(iiph) = 1
enddo


! Warning: once the thermodynamic law has been chosen,
! =======  the remainder of the user subroutine must be modified


!===============================================================================
! 2. Perfect gas
!===============================================================================

if(ieos(iphas).eq.1) then

!===============================================================================
! 2.1. Parameters to be completed by the user
!===============================================================================

! For each phase

! --- Molar mass of the gas (kg/mol)

  if(iccfth.ge.0) then

    if(iphas.eq.1) then

      xmasml = 28.8d-3

    endif

  endif

!===============================================================================
! 2.2. Default laws
!      No user input required.
!===============================================================================

! --- Calculation of the constant gamagp

  if(iccfth.gt.0) then

    ! Gamagp is supposed to be superior or equal to 1.
    ! It is computed at each call, even if this may seem costly,
    !   to be coherent with the "constant gamma" case for which this
    !   constant is not saved. A ''save'' instruction and a test would
    !   be sufficient to avoid computing gamagp at each call if necessary.

    gamagp = 1.d0 + rr/(xmasml*cp0(iphas)-rr)

    if(gamagp.lt.1.d0) then
      write(nfecra,1010) gamagp
      call csexit (1)
    endif

    ! Gamma is returned if required

    if(iccfth.eq.1) then
      gamagr(1) = gamagp
    endif

  endif


! --- Calculation options: constant Cp and Cv (perfect gas)

  if(iccfth.eq.-1) then

    ! The value for the isobaric specific heat Cp0 must be provided in
    !   the user subroutine ''usini1''. The value for the isochoric
    !   specific heat Cv0 is calculated in a subsequent section (from Cp0)

    icp(iphas) = 0
    icv(iphas) = 0


! --- Default initializations (before uscfxi)

!     T0 is positive (this assumption has been checked in
!       the user programme 'verini')

  elseif(iccfth.eq.0) then

    cv0(iphas) = cp0(iphas) - rr/xmasml

    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


! --- Verification of the density

  elseif(iccfth.eq.-2) then

    ! If the density is lower or equal to zero: clipping, write and stop.
    !   Indeed, if this is the case, the thermodynamic computations will
    !   most probably fail.
    ! This call is done at the end of the density calculation (after
    !   a classical clipping and before parallel communications).

    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


! --- Verification of the energy

  elseif(iccfth.eq.-4) then

    ! If the total energy <= zero: clipping, write and stop
    !   Indeed, if this is the case, the thermodynamic computations will
    !   most probably fail.

    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


! --- Calculation of temperature and energy from pressure and density

  elseif(iccfth.eq.12.or.iccfth.eq.60000) then

    ! Verification of the values of the density
    ierr = 0
    do iel = 1, ncel
      if(rtp(iel,irhiph).le.0.d0) then
        write(nfecra,3010)rtp(iel,irhiph),iel
      endif
    enddo
    ! Stop if a negative value is detected (since the density has been
    ! provided by the user, one potential cause is a wrong user
    ! initialization)
    if(ierr.eq.1) then
      call csexit (1)
    endif

    do iel = 1, ncel
      ! Temperature
      sorti1(iel) = xmasml*rtp(iel,ipriph)/(rr*rtp(iel,irhiph))
      ! Total energy
      sorti2(iel) = cv0(iphas)*sorti1(iel)                        &
           + 0.5d0*( rtp(iel,iuiph)**2 + rtp(iel,iviph)**2        &
                                       + rtp(iel,iwiph)**2 )
    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,itkiph) = sorti1(iel)
        rtp(iel,ieniph) = sorti2(iel)
      enddo
    endif


! --- Calculation of density and energy from pressure and temperature:

  elseif(iccfth.eq.13.or.iccfth.eq.100000) then

    ! Verification of the values of the temperature
    ierr = 0
    do iel = 1, ncel
      if(rtp(iel,itkiph).le.0.d0) then
        write(nfecra,2010)rtp(iel,itkiph),iel
      endif
    enddo
    ! Stop if a negative value is detected (since the temperature has been
    ! provided by the user, one potential cause is a wrong user
    ! initialization: a value not provided in Kelvin for example)
    if(ierr.eq.1) then
      call csexit (1)
    endif

    do iel = 1, ncel
      ! Density
      sorti1(iel) = xmasml*rtp(iel,ipriph)/(rr*rtp(iel,itkiph))
      ! Total energy
      sorti2(iel) = cv0(iphas)*rtp(iel,itkiph)                    &
           + 0.5d0*( rtp(iel,iuiph)**2 + rtp(iel,iviph)**2        &
                                       + rtp(iel,iwiph)**2 )
    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,irhiph) = sorti1(iel)
        rtp(iel,ieniph) = sorti2(iel)
      enddo
    endif


! --- Calculation of density and temperature from pressure and energy

  elseif(iccfth.eq.14.or.iccfth.eq.140000) then

    do iel = 1, ncel
      ! Internal energy (to avoid the need to divide by the temperature
      ! to compute density)
      enint = rtp(iel,ieniph)                                     &
               - 0.5d0*( rtp(iel,iuiph)**2                        &
                       + rtp(iel,iviph)**2                        &
                       + rtp(iel,iwiph)**2 )
      ! Density
      sorti1(iel) = rtp(iel,ipriph) / ( (gamagp-1.d0) * enint )
      ! Temperature
      sorti2(iel) = xmasml * (gamagp-1.d0) * enint / rr
    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,irhiph) = sorti1(iel)
        rtp(iel,itkiph) = sorti2(iel)
      enddo
    endif


! --- Calculation of pressure and energy from density and temperature

  elseif(iccfth.eq.23.or.iccfth.eq.150000) then

    do iel = 1, ncel
      ! Pressure
      sorti1(iel) = rtp(iel,irhiph)*rtp(iel,itkiph)*rr/xmasml
      ! Total energy
      sorti2(iel) = cv0(iphas)*rtp(iel,itkiph)                    &
           + 0.5d0*( rtp(iel,iuiph)**2 + rtp(iel,iviph)**2        &
                                       + rtp(iel,iwiph)**2 )
    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,ipriph) = sorti1(iel)
        rtp(iel,ieniph) = sorti2(iel)
      enddo
    endif


! --- Calculation of pressure and temperature from density and energy

  elseif(iccfth.eq.24.or.iccfth.eq.210000) then

    do iel = 1, ncel
      ! Internal energy (to avoid the need to divide by the temperature
      ! to compute density)
      enint = rtp(iel,ieniph)                                     &
               - 0.5d0*( rtp(iel,iuiph)**2                        &
                       + rtp(iel,iviph)**2                        &
                       + rtp(iel,iwiph)**2 )
      ! Pressure
      sorti1(iel) = (gamagp-1.d0) * rtp(iel,irhiph) * enint
      ! Temperature
      sorti2(iel) = xmasml * (gamagp-1.d0) * enint / rr
    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,ipriph) = sorti1(iel)
        rtp(iel,itkiph) = sorti2(iel)
      enddo
    endif


!                     2                            2         P
! --- Calculation of c from pressure and density: c = gamma*---
!                                                           rho

  elseif(iccfth.eq.126) then

    ! Verification of the values of the density
    !   This test can be discarded to reduce the CPU time (if
    !     density is <= 0, the calculation will simply fail)
    !   It is discarded here with iutile = 0
    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

    do iel = 1, ncel
      sorti1(iel) = gamagp * rtp(iel,ipriph) / rtp(iel,irhiph)
    enddo


!                                                              gamma
! --- Calculation of beta from pressure and density: beta = rho

  elseif(iccfth.eq.162) then

    ! Verification of the values of the density
    !   This test can be discarded to reduce the CPU time (if
    !     density is <= 0, the calculation will simply fail)
    !   It is discarded here with iutile = 0
    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

    do iel = 1, ncel
      sorti1(iel) = rtp(iel,irhiph)**gamagp
    enddo


! --- Calculation of the isochoric specific heat

    ! It is a constant: nothing to do


!                                                                  P
! --- Calculation of the entropy from pressure and density: s = --------
!                                                                  gamma
!                                                               rho

  elseif(iccfth.eq.6) then

    ! Verification of the values of the density
    !   This test can be discarded to reduce the CPU time (if
    !     density is <= 0, the calculation will simply fail)
    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

    do iel = 1, ncel
      sorti1(iel) = rtp(iel,ipriph) / (rtp(iel,irhiph)**gamagp)
    enddo


! --- Calculation of 'internal energy - Cv.T'

  elseif(iccfth.eq.7) then

    ! It is zero for a perfect gas

    !   At the cell centers
    do iel = 1, ncel
      sorti1(iel) = 0.d0
    enddo

    !   On the boundary faces
    do ifac = 1, nfabor
      sorti2(ifac) = 0.d0
    enddo


! --- Calculation of the boundary conditions on the face ifac = ifac0

!  -- Wall

  elseif(iccfth.eq.91) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Calculation of the Mach number at the boundary face, using the
    !   cell center velocity projected on the vector normal to the boundary
    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) )

    ! Pressure

    !   A Neumann boundary condition is used. This does not allow to use
    !     the Rusanov scheme, but some stabilization effect is expected.
    !     A test based on the value of coefb at the previous time step
    !     is implemented to avoid oscillating between a rarefaction
    !     situation and a shock configuration from one time step to the
    !     next.

    !   Rarefaction
    if(xmach.lt.0.d0.and.coefb(ifac,iclp).le.1.d0) then

      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))
      else
        ! In case the rarefaction is too strong, a zero Dirichlet value
        !   is used for pressure (the value of coefb is used here as an
        !   indicator and will be modified later in cfxtcl)
        coefb(ifac,iclp) = rinfin
      endif

      !  Shock
    elseif(xmach.gt.0.d0.and.coefb(ifac,iclp).ge.1.d0) then

      coefb(ifac,iclp) = 1.d0 + gamagp*xmach                      &
            *( (gamagp+1.d0)/4.d0*xmach                           &
                + sqrt(1.d0 + (gamagp+1.d0)**2/16.d0*xmach**2) )

      !  Oscillation between rarefaction and shock or zero Mach number
    else
      coefb(ifac,iclp) = 1.d0
    endif


!  -- Symmetry

  elseif(iccfth.eq.90) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! A zero flux condition (homogeneous Neumann condition) is
    !   prescribed by default.
    ! No user input required


!  -- Subsonic inlet with prescribed density and velocity

    ! The subsonic nature of the inlet is postulated.

    ! Further testing may be required here. Contrary to the initial
    !   development, an explicit Dirichlet condition is prescribed for
    !   pressure instead of a Neumann condition (however, the same
    !   physical value for pressure is used).
    ! The advantage of this approach is to allow the use of the Rusanov
    !   scheme to stabilize the user defined inlet conditions.
    ! Moreover, with this approach, coefb does not have to be filled in
    !   here (it is not a major point, since coefb has to be filled in
    !   for the wall boundary condition anyway)
    ! Shall an oscillatory behavior (in time) be observed, it might be
    !   worth trying to add a test to avoid switching between
    !   rarefaction and shock from one time step to the other (just as
    !   for the wall boundary condition).
    ! The relevance of this approach remains to be demonstrated.

  elseif(iccfth.eq.92) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Calculation of the Mach number at the boundary face, using the
    !   cell center velocity projected on the vector normal to the boundary
    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

    ! Pressure: rarefaction wave (Rusanov)
    if(dxmach.le.0.d0) then

      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

      ! Pressure: shock (Rusanov)
    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

    ! This choice overrides the previous Rusanov choice
    coefa(ifac,iclp) = rtp(iel,ipriph)

    ! Total energy
    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)


!  -- Subsonic inlet with prescribed mass and enthalpy flow rates
    ! The quantities prescribed are rho*u and rho*u*h

    ! The subsonic nature of the inlet is postulated.

    ! This section remains to be implemented: stop for the moment

    ! One may proceed as follows:
    !   Pressure computed with a Newton method
    !   Velocity and density computed from pressure
    !   Total energy computed from enthalpy
    !   (written on paper, to be implemented: contact the user support)

  elseif(iccfth.eq.94) then

    ifac = ifac0
    iel  = ifabor(ifac)

    write(nfecra,7000)

    call csexit (1)
    !==========


!  -- Subsonic outlet

    ! The subsonic nature of the inlet is postulated.

  elseif(iccfth.eq.93) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Rarefaction case
    if(coefa(ifac,iclp).le.rtp(iel,ipriph)) then

      ! Density
      coefa(ifac,iclr) = rtp(iel,irhiph)                          &
           * (coefa(ifac,iclp)/rtp(iel,ipriph))**(1.d0/gamagp)

      ! Velocity
      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)

      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)

      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)

      ! Total energy
      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)

    ! Shock
    else

      ! Density
      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) )

      ! Velocity
      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)

      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)

      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)

      ! Total energy
      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)

    endif


! --- Calculation of temperature and energy from pressure and density

    ! It is postulated that the pressure and density values are
    !   strictly positive

  elseif(iccfth.eq.912.or.iccfth.eq.60900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Temperature
    coefa(ifac,iclt) =                                            &
         xmasml*coefa(ifac,iclp)/(rr*coefa(ifac,iclr))

    ! Energie totale
    coefa(ifac,icle) =                                            &
         cv0(iphas)*coefa(ifac,iclt)                              &
         + 0.5d0*( coefa(ifac,iclu)**2                            &
                 + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2 )


! --- Calculation of density and energy from pressure and temperature

  elseif(iccfth.eq.913.or.iccfth.eq.100900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Density
    coefa(ifac,iclr) =                                            &
         xmasml*coefa(ifac,iclp)/(rr*coefa(ifac,iclt))

    ! Total energy
    coefa(ifac,icle) =                                            &
         cv0(iphas)*coefa(ifac,iclt)                              &
         + 0.5d0*( coefa(ifac,iclu)**2                            &
                 + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2 )


! --- Calculation of density and temperature from pressure and total energy

  elseif(iccfth.eq.914.or.iccfth.eq.140900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Density
    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 ) ) )

    ! Temperature
    coefa(ifac,iclt)=                                             &
         xmasml*coefa(ifac,iclp)/(rr*coefa(ifac,iclr))


! --- Calculation of pressure and energy from density and temperature

  elseif(iccfth.eq.923.or.iccfth.eq.150900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Pressure
    coefa(ifac,iclp) = coefa(ifac,iclr)*rr/xmasml                 &
                                       *coefa(ifac,iclt)

    ! Total energy
    coefa(ifac,icle) = cv0(iphas) * coefa(ifac,iclt)              &
         + 0.5d0*( coefa(ifac,iclu)**2                            &
                 + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2 )


! --- Calculation of pressure and temperature from density and energy

  elseif(iccfth.eq.924.or.iccfth.eq.210900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Pressure
    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 ) )


    ! Temperature
    coefa(ifac,iclt)=                                             &
         xmasml*coefa(ifac,iclp)/(rr*coefa(ifac,iclr))


! --- End of the treatment of the perfect gas
  endif


!===============================================================================
! 3. Perfect gas with variable gamma
!===============================================================================

! This section requires further checking and testing

elseif(ieos(iphas).eq.2) then

!===============================================================================

!===============================================================================
! 3.1. Parameters to be completed by the user
!===============================================================================



! --- Examples (to be copied and adapted in section ''3.1. Parameters ...''

!-------------------------------------------------------------------------------
! This test allows the user to ensure that the version of this subroutine
!   used is that from his case definition, and not that from the library.

  if(0.eq.1) then

! --- Ex. 1: Perfect gas containing 3 components
!     Molar mass, gamma

    ! Phase
  if(iphas.eq.1) then

    ! Molar mass of the components (kg/mol)
    cstgr(1)  = 18.d-3
    cstgr(2)  = 32.d-3
    cstgr(3)  = 28.d-3

    if(iccfth.gt.0) then

      ! Calculation of the molar mass of the mixture at cell centers
      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

      ! Calculation of the equivalent gamma of the mixture at cell centers
      do iel = 1, ncel
        gamagr(iel) = propce(iel,ipproc(icp(iphas)))              &
           / ( propce(iel,ipproc(icp(iphas))) - rr/xmasmr(iel) )
      enddo

    endif

  endif

  endif

!-------------------------------------------------------------------------------

! End of the examples


! Verification of the values of gamagr: gamagr >= 1., otherwise stop

  ierr = 0

  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

  if(ierr.eq.1) then
    call csexit (1)
  endif


! --- Calculation options: variable Cp and Cv
!     (isobaric and isochoric specific heat)

  if(iccfth.eq.-1) then

    icp(iphas) = 1
    cp0(iphas) = epzero
    icv(iphas) = 1
    cv0(iphas) = epzero


! Default initializations

  elseif(iccfth.eq.0) then

    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


! --- Calculation of temperature and energy from pressure and density

  elseif(iccfth.eq.12) then

    do iel = 1, ncel

      ! Temperature
      sorti1(iel) =                                               &
           xmasmr(iel)/rr*rtp(iel,ipriph)/rtp(iel,irhiph)

      ! Total energy
      sorti2(iel) = propce(iel,ipproc(icv(iphas)))*sorti1(iel)    &
    + 0.5d0*( rtp(iel,iuiph)**2                                   &
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 )

    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,itkiph) = sorti1(iel)
        rtp(iel,ieniph) = sorti2(iel)
      enddo
    endif


! --- Calculation of density and energy from pressure and temperature:

  elseif(iccfth.eq.13) then

    do iel = 1, ncel

      ! Density
      sorti1(iel) =                                               &
           xmasmr(iel)/rr*rtp(iel,ipriph)/rtp(iel,itkiph)

      ! Total energy
      sorti2(iel) =                                               &
           propce(iel,ipproc(icv(iphas)))*rtp(iel,itkiph)         &
    + 0.5d0*( rtp(iel,iuiph)**2                                   &
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 )

    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,irhiph) = sorti1(iel)
        rtp(iel,ieniph) = sorti2(iel)
      enddo
    endif


! --- Calculation of density and temperature from pressure and energy

  elseif(iccfth.eq.14) then

    do iel = 1, ncel

      ! Density
      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 ) )

      ! Temperature
      sorti2(iel) = xmasmr(iel)/rr*rtp(iel,ipriph)/sorti1(iel)

    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,irhiph) = sorti1(iel)
        rtp(iel,itkiph) = sorti2(iel)
      enddo
    endif


! --- Calculation of pressure and energy from density and temperature

  elseif(iccfth.eq.23) then

    do iel = 1, ncel

      ! Pressure
      sorti1(iel) =                                               &
           rtp(iel,irhiph)*rr/xmasmr(iel)*rtp(iel,itkiph)

      ! Total energy
      sorti2(iel) =                                               &
           propce(iel,ipproc(icv(iphas)))*rtp(iel,itkiph)         &
    + 0.5d0*( rtp(iel,iuiph)**2                                   &
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 )

    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,ipriph) = sorti1(iel)
        rtp(iel,ieniph) = sorti2(iel)
      enddo
    endif


! --- Calculation of pressure and temperature from density and energy

  elseif(iccfth.eq.24) then

    do iel = 1, ncel

      ! Pressure
      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 ) )

      ! Temperature
      sorti2(iel) = xmasmr(iel)/rr*sorti1(iel)/rtp(iel,irhiph)

    enddo

    ! Transfer to the array rtp
    if(imodif.gt.0) then
      do iel = 1, ncel
        rtp(iel,ipriph) = sorti1(iel)
        rtp(iel,itkiph) = sorti2(iel)
      enddo
    endif

!                     2                            2         P
! --- Calculation of c from pressure and density: c = gamma*---
!                                                           rho

  elseif(iccfth.eq.126) then

    do iel = 1, ncel

      ! Verification of the positivity of the pressure
      if(rtp(iel,ipriph).lt.0.d0) then
        write(nfecra,1110) iel , rtp(iel,ipriph)
        ierr = 1

      ! Verification of the positivity of the density
      elseif(rtp(iel,irhiph).le.0.d0) then
        write(nfecra,1120) iel , rtp(iel,irhiph)
        ierr = 1

      else

        ! Computation
        sorti1(iel) =                                             &
             gamagr(iel) * rtp(iel,ipriph) / rtp(iel,irhiph)

      endif

    enddo

    ! Stop if error detected
    if(ierr.eq.1) call csexit (1)


!                                                              gamma
! --- Calculation of beta from pressure and density: beta = rho

  elseif(iccfth.eq.162) then

    do iel = 1, ncel

      ! Verification of the positivity of the density
      if(rtp(iel,irhiph).lt.0.d0) then
        write(nfecra,1220) iel , rtp(iel,irhiph)
        ierr = 1

      else

        ! Computation
        sorti1(iel) = rtp(iel,irhiph)**gamagr(iel)

      endif

    enddo

    ! Stop if error detected
    if(ierr.eq.1) call csexit (1)


! --- Calculation of the isochoric specific heat: Cv = Cp - R/M

  elseif(iccfth.eq.432) then

    do iel = 1, ncel

      sorti1(iel) = propce(iel,ipproc(icp(iphas)))-rr/xmasmr(iel)

    enddo

    ! Stop if error detected (kept by consistance with other sections)
    if(ierr.eq.1) call csexit (1)

!                                                                  P
! --- Calculation of the entropy from pressure and density: s = --------
!                                                                  gamma
!                                                               rho

  elseif(iccfth.eq.6) then

    do iel = 1, ncel

      ! Verification of the positivity of the pressure
      if(rtp(iel,ipriph).lt.0.d0) then
        write(nfecra,1310) iel , rtp(iel,ipriph)
        ierr = 1

      ! Verification of the positivity of the density
      elseif(rtp(iel,irhiph).le.0.d0) then
        write(nfecra,1320) iel , rtp(iel,irhiph)
        ierr = 1

      else

        ! Computation
        sorti1(iel) =                                             &
             rtp(iel,ipriph) / (rtp(iel,irhiph)**gamagr(iel))

      endif

    enddo

    ! Stop if error detected
    if(ierr.eq.1) call csexit (1)


! --- Calculation of 'internal energy - Cv.T'

  elseif(iccfth.eq.7) then

    ! It is zero for a perfect gas

    !   At the cell centers
    do iel = 1, ncel
      sorti1(iel) = 0.d0
    enddo

    !   On the boundary faces
    do ifac = 1, nfabor
      sorti2(ifac) = 0.d0
    enddo

    ! Stop if error detected (kept by consistance with other sections)
    if(ierr.eq.1) call csexit (1)


! --- Calculation of the boundary conditions on the face ifac = ifac0

!  -- Wall/symmetry

  elseif(iccfth.eq.91) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Calculation of the Mach number at the boundary face, using the
    !   cell center velocity projected on the vector normal to the boundary
    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) )

    coefa(ifac,iclp) = 0.d0

    ! Pression and entropy: rarefaction

    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

    elseif(xmach.le.2.d0/(1.d0-gamagr(iel)) ) then
      coefb(ifac,iclp) = 0.d0
      coefb(ifac,iclt) = 1.d0

      ! Pressure and entropy: shock

    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

    ! Total energy: 'internal energy - Cv T'

    coefa(ifac,icle) = 0.d0

    ! Stop if error detected
    if(ierr.eq.1) call csexit (1)


!  -- Inlet

  elseif(iccfth.eq.92) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Calculation of the Mach number at the boundary face, using the
    !   cell center velocity projected on the vector normal to the boundary
    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

    ! Pressure: rarefaction wave
    if(dxmach.le.0.d0) then

      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

    ! Pressure: shock
    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

    ! This choice overrides the previous Rusanov choice
    coefa(ifac,iclp) = rtp(iel,ipriph)

    ! Total energy
    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)

!  -- Outlet

  elseif(iccfth.eq.93) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Calculation of the Mach number at the boundary face, using the
    !   cell center velocity projected on the vector normal to the boundary
    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))

    ! Supersonic outlet: Dirichlet for all variables
    if(xmach.ge.1.d0) then
      do ivar = 1, nvar
        coefa(ifac,iclrtp(ivar,icoef)) = rtp(iel,ivar)
      enddo

      ! Entropy
      coefa(ifac,iclt) =                                          &
           rtp(iel,ipriph)/rtp(iel,irhiph)**gamagr(iel)

    ! Subsonic outlet
    elseif(xmach.lt.1.d0 .and. xmach.ge.0.d0) then

      ! Rarefaction:
      if(coefa(ifac,iclp).le.rtp(iel,ipriph)) then

        ! Density
        coefa(ifac,iclr) = rtp(iel,irhiph)                        &
             * (coefa(ifac,iclp)/rtp(iel,ipriph))                 &
                **(1.d0/gamagr(iel))

        ! Velocity
        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)

        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)

        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)

        ! Total energy
        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)

        ! Entropy
        coefa(ifac,iclt) = coefa(ifac,iclp)                       &
                             /coefa(ifac,iclr)**gamagr(iel)

      ! Shock:
      else

        ! Density
        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) )

        ! Velocity
        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)

        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)

        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)

        ! Total energy
        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)

        ! Entropy
        coefa(ifac,iclt) = coefa(ifac,iclp)                       &
                             /coefa(ifac,iclr)**gamagr(iel)

      endif

    else
      write(nfecra,*) 'iccfth = ',iccfth,'  Mach = ',xmach
      ierr = 1
    endif

    if(ierr.eq.1) call csexit (1)


! --- Calculation of temperature and energy from pressure and density

  elseif(iccfth.eq.912.or.iccfth.eq.60900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Temperature
    coefa(ifac,iclt) = xmasmr(iel)/rr*coefa(ifac,iclp)            &
                                        /coefa(ifac,iclr)

    ! Total energy
    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)


! --- Calculation of density and energy from pressure and temperature

  elseif(iccfth.eq.913.or.iccfth.eq.100900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Density
    coefa(ifac,iclr) = xmasmr(iel)/rr*coefa(ifac,iclp)            &
                                       /coefa(ifac,iclt)

    ! Total energy
    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)


! --- Calculation of density and temperature from pressure and total energy

  elseif(iccfth.eq.914.or.iccfth.eq.140900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Density
    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 ))

    ! Temperature
    coefa(ifac,iclt)= xmasmr(iel)/rr*coefa(ifac,iclp)             &
                                       /coefa(ifac,iclr)


! --- Calculation of pressure and energy from density and temperature

  elseif(iccfth.eq.923.or.iccfth.eq.150900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Pressure
    coefa(ifac,iclp) = coefa(ifac,iclr)*rr/xmasmr(iel)            &
                                       *coefa(ifac,iclt)

    ! Total energy
    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)


! --- Calculation of pressure and temperature from density and energy

  elseif(iccfth.eq.924.or.iccfth.eq.210900) then

    ifac = ifac0
    iel  = ifabor(ifac)

    ! Pressure
    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 ) )


    ! Temperature
    coefa(ifac,iclt)= xmasmr(iel)/rr*coefa(ifac,iclp)             &
                                       /coefa(ifac,iclr)


! --- End of perfect gas with variable gamma
  endif

! --- End of test on the thermodynamic laws
endif


!--------
! Formats
!--------

 1010 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     Gamma = ',e12.4   ,/,                                      &
'@     Gamma must be a real number greater or equal to 1.',/,     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 1020 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     In cell ',i10   ,', Gamma = ',e12.4   ,/,                  &
'@     Gamma must be a real number greater or equal to 1.',/,     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 2010 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     The computation of density failed.',/,                     &
'@',/,                                                            &
'@     Temperature = ',e12.4   ,' in cell ',i10  ,/,              &
'@     Temperature must be strictly positive.',/,                 &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 3010 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     The computation of temperature failed.',/,                 &
'@',/,                                                            &
'@     Density = ',e12.4   ,' in cell ',i10  ,/,                  &
'@     Density must be strictly positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 4010 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the squared speed of sound failed.',/,  &
'@',/,                                                            &
'@     Density = ',e12.4   ,' in cell ',i10  ,/,                  &
'@     Density must be strictly positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 4020 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the variable beta failed.',/,           &
'@',/,                                                            &
'@     Density = ',e12.4   ,' in cell ',i10  ,/,                  &
'@     Density must be strictly positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 4030 format(                                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the entropy failed.',/,                 &
'@',/,                                                            &
'@     Density = ',e12.4   ,' in cell ',i10  ,/,                  &
'@     Density must be strictly positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 7000 format (                                                    &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     The boundary condition of the type ''prescribed mass',/,   &
'@     and enthalpy flow rates '' is not available in the ',/,    &
'@     current release.',/,                                       &
'@',/,                                                            &
'@     Modify the user subroutine ''uscfth''.',/,                 &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 8000 format (                                                    &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     Negative values of the density were encountered ',/,       &
'@     in ',i10   ,' cells.',/,                                   &
'@     The density was clipped at ',e12.4  ,/                     &
'@     The run was stopped.',/,                                   &
'@',/,                                                            &
'@     If it is desired to continue the run in spite of this ',/, &
'@     behavior, it is possible to force a standard clipping ',/, &
'@     by setting a minimum value for the density variable in',/, &
'@     the GUI or in the user subroutine ''usini1'' (set the ',/, &
'@     scamin value associated to the variable ',/,               &
'@     isca(irho(iphas)).',/,                                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 8100 format (                                                          &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with constant gamma.',/,                 &
'@',/,                                                            &
'@     Negative values of the internal energy were encountered',/,&
'@     in ',i10   ,' cells.',/,                                   &
'@     The internal energy  was clipped at ',e12.4  ,/            &
'@     The run was stopped.',/,                                   &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)


! The following formats may be discarded if or when the
! gamma variable option will have been fixed


 1110 format(                                                           &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with variable gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the squared speed of sound failed.',/,  &
'@',/,                                                            &
'@     In cell ',i10   ,' Pressure = ',e12.4   ,/,                &
'@     Pressure must be positive.',/,                             &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 1120 format(                                                           &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with variable gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the squared speed of sound failed.',/,  &
'@',/,                                                            &
'@     In cell ',i10   ,' Density = ',e12.4   ,/,                 &
'@     Density must be strictly positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 1220 format(                                                           &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with variable gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the variable beta failed.',/,           &
'@',/,                                                            &
'@     In cell ',i10   ,' Density = ',e12.4   ,/,                 &
'@     Density must be strictly positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 1310 format(                                                           &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with variable gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the entropy failed.',/,                 &
'@',/,                                                            &
'@     In cell ',i10   ,' Pressure = ',e12.4   ,/,                &
'@     Pressure must be positive.',/,                             &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)
 1320 format(                                                           &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/,                                                            &
'@ @@ WARNING:    stop in thermodynamics computations',/,         &
'@    =======',/,                                                 &
'@     Error encountered in the user subroutine ''uscfth'', ',/,  &
'@       for perfect gas with variable gamma.',/,                 &
'@',/,                                                            &
'@     The computation of the entropy failed.',/,                 &
'@',/,                                                            &
'@     In cell ',i10   ,' Density = ',e12.4   ,/,                 &
'@     Density must be striclty positive.',/,                     &
'@',/,                                                            &
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
'@',/)


!----
! End
!----

return
end subroutine
