SUBROUTINE etat_saha(pp,tt,xchim,deriv,
	1 ro,drop,drot,drox,u,dup,dut,dux,
	2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)

! routine private du module mod_etat

! Interface of the equation of state SAHA-S3 with CESAM2k.
! Computations are performed for current metallicity Z(r,t).
!
! input :
!	p : pressure
!	t : temperature
!	xchim : chemical composition
!	deriv=.FALSE. : avoide computations of some derivatives

! output :
!	ro : density and derivatives
!	u : internal energy and derivatives
!   ...

      IMPLICIT NONE

	REAL*8, INTENT(in), DIMENSION(:) :: xchim
	REAL*8, INTENT(in) :: pp, tt
	LOGICAL, INTENT(in) :: deriv

	REAL*8, INTENT(out) :: ro,drop,drot,drox,u,dup,dut,dux
	REAL*8, INTENT(out) :: delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx
	REAL*8, INTENT(out) :: gradad,dgradadp,dgradadt,dgradadx,alfa,beta, gamma1

	REAL*8 ro1,drop1,drot1,drox1,u1,dup1,dut1,dux1
	REAL*8 delta1,deltap1,deltat1,deltax1,cp1,dcpp1,dcpt1,dcpx1
	REAL*8 gradad1,dgradadp1,dgradadt1,dgradadx1,alfa1,beta1, gamma1_1

	REAL*8 ro2,drop2,drot2,drox2,u2,dup2,dut2,dux2
	REAL*8 delta2,deltap2,deltat2,deltax2,cp2,dcpp2,dcpt2,dcpx2
	REAL*8 gradad2,dgradadp2,dgradadt2,dgradadx2,alfa2,beta2, gamma1_2

	real*8 tx2
	REAL*8, DIMENSION(23) :: Y_var,Y_var1,Y_var2

	REAL*8  x, z

	LOGICAL, SAVE :: init=.TRUE.

! ________________________________________________________________________________

	x=xchim(1)

! Computing current metallicity  z
	IF(nchim > 1)THEN
	 z=1.d0-SUM(xchim(1:ihe4))
	ELSE
	 z=z0
	ENDIF

! Range checking
	IF (Z < 0.0d0 .OR. Z > 0.1d0) THEN
	 PRINT*, 'Metallicity Z is out of the range of the SAHA-S3 tables'
	 PRINT*, 'Z = ', z
	 PRINT*, 'Ztab1 = ', Ztab1, ',   Ztab2 = ', Ztab2
	 STOP
	ENDIF

! Reading the SAHA_S3 tables for two z's
	IF(init)THEN
	 init=.FALSE.
	 PRINT*, '------ Equation of state SAHA-S3 ------'
	 CALL Z_read2_saha
	ENDIF

! Computing derivatives for Z=Ztab1
	xz=>xz1
	Beos=>Beos1
	CALL derivatives(pp,tt,x,deriv,
	1 ro1,drop1,drot1,drox1,u1,dup1,dut1,dux1,
	2 delta1,deltap1,deltat1,deltax1,cp1,dcpp1,dcpt1,dcpx1,
	3 gradad1,dgradadp1,dgradadt1,dgradadx1,alfa1,beta1,gamma1_1)

	Y_var1=(/ro1,drop1,drot1,drox1,u1,dup1,dut1,dux1,
	1 delta1,deltap1,deltat1,deltax1,cp1,dcpp1,dcpt1,dcpx1,
	2 gradad1,dgradadp1,dgradadt1,dgradadx1,alfa1,beta1,gamma1_1/)

! Computing derivatives for Z=Ztab2
	xz=>xz2
	Beos=>Beos2
	CALL derivatives(pp,tt,x,deriv,
	1 ro2,drop2,drot2,drox2,u2,dup2,dut2,dux2,
	2 delta2,deltap2,deltat2,deltax2,cp2,dcpp2,dcpt2,dcpx2,
	3 gradad2,dgradadp2,dgradadt2,dgradadx2,alfa2,beta2,gamma1_2)

	Y_var2=(/ro2,drop2,drot2,drox2,u2,dup2,dut2,dux2,
	1 delta2,deltap2,deltat2,deltax2,cp2,dcpp2,dcpt2,dcpx2,
	2 gradad2,dgradadp2,dgradadt2,dgradadx2,alfa2,beta2,gamma1_2/)

! Interpolating derivatives for given metallicity Z
	tx2=(Z-Ztab1)/(Ztab2-Ztab1)
	Y_var=(1-tx2)*Y_var1+tx2*Y_var2

	ro=Y_var(1)
	drop=Y_var(2)
	drot=Y_var(3)
	drox=Y_var(4)

	u=Y_var(5)
	dup=Y_var(6)
	dut=Y_var(7)
	dux=Y_var(8)

	delta=Y_var(9)
	deltap=Y_var(10)
	deltat=Y_var(11)
	deltax=Y_var(12)

	cp=Y_var(13)
	dcpp=Y_var(14)
	dcpt=Y_var(15)
	dcpx=Y_var(16)

	gradad=Y_var(17)
	dgradadp=Y_var(18)
	dgradadt=Y_var(19)
	dgradadx=Y_var(20)

	alfa=Y_var(21)
	beta=Y_var(22)
	gamma1=Y_var(23)

	RETURN

	END SUBROUTINE etat_saha