c*************************************************************************

	PROGRAM test_saha

c	test pour saha
c	Auteur: P. Morel, Departement J.D. Cassini, O.C.A.
c	CESAM2k

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

	use mod_donnees, ONLY : abe7, ac12, ac13, ah, ah2, ahe3, ahe4,
	1 ali7, an14, an15, ao16, ao17, afe56, lit_nl, nbelem, nchim,
	2 nom_elem, nom_fich2, nucleo, w_rot, zi
	use mod_etat, ONLY : saha
	use mod_exploit, ONLY : read_ascii, var	
	use mod_kind

	IMPLICIT NONE

	INTEGER, PARAMETER :: pnosc=4500, pnelem=15

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: xchim, ioni
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: t, ro, ychim, z_bar
	REAL (kind=dp) :: eta, nel

	INTEGER :: i, iconst, itot, j, k, nglob, nvar
	
	CHARACTER (len=4), ALLOCATABLE, DIMENSION(:) :: nom_elem_lu
	CHARACTER (len=40) :: nom_fich
	CHARACTER (len=80), DIMENSION(4) :: abid	
	CHARACTER (len=80) :: modele
	
c------------------------------------------------------------------------
	
2000	FORMAT(8es10.3)
	
c	lecture de test.don

	nom_fich2='test' ; CALL lit_nl(w_rot)	
		
c	PRINT*,'entrer le nom du fichier sans l''extension .osc'
c	PRINT*,'Exemple: evry, Sortie: zzzz'
c	READ*,nom_fich
	nom_fich='test.ascii'
	
	CALL read_ascii(nom_fich,itot,nglob,nvar,abid)
	
c	initialisation de nucleo: masses des noyaux

	ALLOCATE(t(itot),ro(itot),xchim(nchim,itot),nucleo(nbelem),
	1 zi(nbelem))
	DO i=1,itot
	 t(i)=var(3,i) ; ro(i)=var(5,i)
	 DO j=1,nchim
	  xchim(j,i)=var(nvar+j,i)	!max(var(ivar+j,i),1.d-30)	
	 ENDDO
	ENDDO
c	PRINT*,itot
c	WRITE(*,2000)t
c	pause
	
	DO i=1,nbelem
	 IF(nom_elem(i) == ' H1 ' )THEN
	  nucleo(i)=ah ; zi(i)=1
	 ELSEIF(nom_elem(i) ==  ' H2 ')THEN
	  nucleo(i)=ah2 ; zi(i)=1
	 ELSEIF(nom_elem(i) == 'He3 ')THEN
	  nucleo(i)=ahe3 ; zi(i)=2
	 ELSEIF(nom_elem(i) == 'He4 ')THEN
	  nucleo(i)=ahe4 ; zi(i)=2
	 ELSEIF(nom_elem(i) == 'Li7 ')THEN
	  nucleo(i)=ali7 ; zi(i)=3
	 ELSEIF(nom_elem(i) == 'Be7 ')THEN
	  nucleo(i)=abe7 ; zi(i)=4
	 ELSEIF(nom_elem(i) == 'C12 ')THEN
	  nucleo(i)=ac12 ; zi(i)=6   	  
	 ELSEIF(nom_elem(i) == 'C13 ')THEN
	  nucleo(i)=ac13 ; zi(i)=6	  
	 ELSEIF(nom_elem(i) == 'N14 ')THEN
	  nucleo(i)=an14 ; zi(i)=7
	 ELSEIF(nom_elem(i) == 'N15 ')THEN
	  nucleo(i)=an15 ; zi(i)=7
	 ELSEIF(nom_elem(i) == 'O16 ')THEN
	  nucleo(i)=ao16 ; zi(i)=8
	 ELSEIF(nom_elem(i) == 'O17 ')THEN
	  nucleo(i)=ao17 ; zi(i)=8
	 ELSEIF(nom_elem(i) == 'Fe56')THEN
	  nucleo(i)=afe56 ; zi(i)=26
	 ELSEIF(nom_elem(i) == 'Fe55')THEN
	  nucleo(i)=55 ; zi(i)=26
	 ELSEIF(nom_elem(i) == 'Si28')THEN
	  nucleo(i)=28 ; zi(i)=14
	 ELSEIF(nom_elem(i) == ' Ex ')THEN
	  nucleo(i)=28 ; zi(i)=13
	 ELSEIF(nom_elem(i) == ' B11')THEN
	  nucleo(i)=11 ; zi(i)=5	    
	 ENDIF
	ENDDO
	
	i=0
	DO j=1,nchim
	 i=max(i,nint(zi(j)))
	ENDDO
	ALLOCATE(ychim(nchim),ioni(0:i,nchim),z_bar(nchim))

	DO i=1,itot,200
	 ychim(:)=xchim(:,i)/nucleo(:)
	 CALL saha(ychim,t(i),ro(i),ioni,z_bar,nel,eta)
	 PRINT*,'t,ro,nel,eta'
	 WRITE(*,2000)t(i),ro(i),nel,eta
	 PRINT*,'charges moyennes'
	 WRITE(*,"(13a4)")(nom_elem(j),j=1,nbelem)
	 WRITE(*,2000)z_bar	  
	 DO k=1,nchim
	  PRINT*,nom_elem(k)
	  WRITE(*,2000)(ioni(j,k),j=0,nint(zi(k)))
	 ENDDO

	 PAUSE
	ENDDO
	
	STOP
			 				
	END PROGRAM test_saha