!  Z_read2_saha.f90
!
!  FUNCTIONS:
!  Z_read2_saha - read SAHA-S3 data tables eos_saha_***.dat, prepared for CESAM;
!                 the result is written in xz1(mx,nt,nqs,nclm) and xz2(mx,nt,nqs,nclm).
!               - read SAHA_S3 tables of Bspline coefficients eos_saha_Bspl_***.dat;
!                 the result is written in Beos1(iqs,it,imv,ix) and Beos2(iqs,it,imv,ix).
!
!  Files eos_saha_***.dat and eos_saha_Bspl_***.dat are saved in the directory CESAM-DATA.
!  Asterisks *** in the names of files denote metallicity Z for the tables:
!  000 Z=0; 010  Z=0.010; 015  Z=0.015; 020  Z=0.020.
!
!********************************************************************************************
!
!  SUBROUTINE: Z_read2_saha
!
!  sense of variables in a line of table eos_saha_***.dat
!   P, dyne/cm^2                   u, egr/g          du/drho, CGS             {du/dT}_P=cp, CGS
!   {dlgP/dlgrho}_T=chi_rho        delta             Gamma1                   grad_ad
!   Gamma3-1                       rho, g/cm^3       {dlgP/dlgX}_T,rho        {du/dx}_T,rho, CGS
!
!********************************************************************************************

	SUBROUTINE Z_read2_saha

	USE mod_donnees, ONLY : Z0
	USE mod_numerique, ONLY : entre
	USE mod_variables, ONLY : sortie

	IMPLICIT NONE

! Variables

	INTEGER ::i, idep, ix, it, iqs, imv                 ! loop indices
	LOGICAL :: ok

	CHARACTER (len=4) :: ch3='.dat'
	CHARACTER (len=5) :: ch2='Bspl_'
	CHARACTER (len=9) :: ch1='eos_saha_'
	CHARACTER (len=3), DIMENSION(4) :: ch=(/ '000', '010', '015', '020' /)
	CHARACTER (len=30) :: lst
	CHARACTER (len=50) :: chain

c---------------------------------------------------------------------------


! Metallicity Ztab for the two SAHA-S3 tables used for Z-interpolation
! and for two tables of Bspline coefficients
! Ztab can be 0.000, 0.010, 0.015, 0.020

	IF(entre(0.00d0,0.01d0,Z0))THEN
	 Ztab1=0.00d0
	 Ztab2=0.010d0
	 idep=1
	ELSEIF(entre(0.01d0,0.015d0,Z0))THEN
	 Ztab1=0.010d0
	 Ztab2=0.015d0
	 idep=2
	ELSEIF(entre(0.015d0,0.020d0,Z0))THEN
	 Ztab1=0.015d0
	 Ztab2=0.020d0
	 idep=3
	ELSE
	 WRITE(*,5)Z0
5	 FORMAT('etat_saha, SORTIE: métallicité Z0=',es10.3,
	1 ', en dehors de [0.00,0.02]')
	 CALL sortie('read2_saha 1')
	ENDIF

c vérification de l'existence et décompression des tables d'EOS dans SUN_STAR_DATA
c lecture et exploitation
	WRITE(*,6)Z0
6	FORMAT('Recherche des données pour la métallicité Z0=',es10.3)
	DO i=idep,idep+1
	 lst=ch1//ch(i)//ch3
	 chain=TRIM(nom_chemin)//TRIM(lst)
	 INQUIRE(file=TRIM(chain)//'.gz',exist=ok)
	 IF(ok)THEN
	  WRITE(*,11)TRIM(chain)//'.gz'
11	  FORMAT("décompression du fichier d'EOS: ",a)
	  CALL SYSTEM('gunzip '//TRIM(chain)//'.gz')
	 ELSE
	  INQUIRE(file=TRIM(chain),exist=ok)
	  IF(.NOT.ok)THEN
	   WRITE(*,12)TRIM(chain)
12	   FORMAT("Arrêt, fichier d'EOS non trouvé:",a)
	   CALL sortie('read2_saha 2')
          ENDIF
         ENDIF

! Reading files with SAHA-S3 data for given metallicities
	 OPEN(unit=60,form='formatted',status='old',action='read',file=TRIM(chain))
	 do ix=1,mx
	  do it=1,nt
	   do iqs=1,nqs
	    IF(i == idep)THEN
	     read (60,1005,ERR=2,END=3) (xz1(ix,imv,it,iqs), imv=1,mv)
	    ELSE
	     read (60,1005,ERR=2,END=3) (xz2(ix,imv,it,iqs), imv=1,mv)
	    ENDIF
1005	    format(1x,12e18.10)
	   enddo
	  enddo
	 enddo
	 close(unit=60)

c recompression du fichier d'EOS
	 CALL SYSTEM('gzip '//TRIM(chain))
	 WRITE(*,13)TRIM(chain)
13	 FORMAT("recompression du fichier d'EOS: ",a)

c vérification de l'existence et décompression des tables d'EOS BSPL dans SUN_STAR_DATA
c lecture et exploitation
	 lst=ch1//ch2//ch(i)//ch3
	 chain=TRIM(nom_chemin)//TRIM(lst)
	 INQUIRE(file=TRIM(chain)//'.gz',exist=ok)
	 IF(ok)THEN
	  WRITE(*,11)TRIM(chain)//'.gz'
	  CALL SYSTEM('gunzip '//TRIM(chain)//'.gz')
	 ELSE
	  INQUIRE(file=TRIM(chain),exist=ok)
	  IF(.NOT.ok)THEN
	   WRITE(*,12)TRIM(chain)
	   CALL sortie('read2_saha 3')
	  ENDIF
	 ENDIF

! Reading files with tables of Bspline coefficients for given metallicities Ztab1 and Ztab2
	 OPEN(unit=60,form='formatted',status='old',action='read',file=TRIM(chain))
	 do ix=1,mx
	  do it=1,nt
	   do iqs=1,nqs
	    IF(i == idep)THEN
	     read (60,1006,ERR=2,END=3) (Beos1(iqs,it,imv,ix), imv=1,mv)
	    ELSE
	     read (60,1006,ERR=2,END=3) (Beos2(iqs,it,imv,ix), imv=1,mv)
	    ENDIF
1006	    format(12e25.16)
	   enddo
	  enddo
	 enddo
	 close(unit=60)

c compression du fichier d'EOS BSPL
	 CALL SYSTEM('gzip '//TRIM(chain))
	 WRITE(*,13)TRIM(chain)
	ENDDO

! initialization B-spline knots
	call aug_knot(Xkn,Xknot)
	call aug_knot(lgTkn,Tknot)
	call aug_knot(lgQskn,Qsknot)

	RETURN

2	PRINT*,'erreur de lecture dans le fichier: ',TRIM(chain)
	CALL sortie('read2_saha 4')

3	PRINT*,'fin de fichier dans le fichier: ',TRIM(chain)
	CALL sortie('read2_saha 5')

	END SUBROUTINE Z_read2_saha