SUBROUTINE rhoof_PTX(lgT,lgP,X,eos,deos)
! give RHO as function of P, T, X according to SAHA-S3 EOS.
! suppose to be a part of module mod_etat_saha
! equivalent to rhoofp from CESAM etat_opalZ code

	IMPLICIT NONE
	REAL*8, INTENT(in) :: lgT,lgP,X

	REAL*8, INTENT(out), DIMENSION(mv) :: eos
	REAL*8, DIMENSION(mv,3) :: deos

c	REAL*8, DIMENSION(mv) :: vars
	REAL*8, DIMENSION(nqs) :: rho_line, P_line
	REAL*8 lgQs
	REAL*8 lgrho, lgrho_i, lgP_i
	INTEGER iter

!  save all available pairs of P and rho for given T and X
	CALL int2_onlgQs(X,lgT,P_line,rho_line)

!  check this line for universality
	IF ((lgP-P_line(1)*(1.d0-SIGN(tol_extr,P_line(1))))*
	1 (lgP-P_line(nqs)*(1.d0+SIGN(tol_extr,P_line(nqs))))<=0) THEN
   !determine rho using iterations
	 CALL interp1(P_line,rho_line,lgP,lgrho)
	 lgQs=lgrho-2.25d0*(lgT-6.0d0)

	 CALL Bspline2f(x,lgT,lgQs,eos)   ! computes only eos(1) and eos(5)
	 lgP_i=log10(eos(1))
	 iter=1

	 DO WHILE (ABS(lgP-lgP_i)>1d-10)
	  lgrho_i=lgrho + (lgP-lgP_i)/eos(5)
	  lgQs=lgrho_i-2.25d0*(lgT-6.0d0)
	  CALL Bspline2f(x,lgT,lgQs,eos) ! computes only eos(1) and eos(5)
	  lgP_i=log10(eos(1))
	  lgrho=lgrho_i
	  iter=iter+1
	  IF (iter.gt.20) THEN
	   PRINT*, 'too many iter', lgP,lgP_i
	   PAUSE
	  END IF
	 ENDDO

	ELSE
	 PRINT*, 'SAHA_S3 \ rhoof_PTX: out of lgP boundary'
	 PRINT*, 'P_line(1)=', P_line(1), 'P_line(nqs)=', P_line(nqs), 'lgP=',lgP
	 STOP
	END IF

	lgQs=lgrho-2.25d0*(lgT-6.0d0)
	CALL Bspline(x,lgT,lgQs,eos,deos) ! interpolate all functions eos and derivatives deos

	RETURN

	END SUBROUTINE rhoof_PTX