INCLUDE 'mod_opa'

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

	PROGRAM test_opa_opal2

c	test pour les routines d'opacité opal type2

c	Auteur: P.Morel, Département J.D. Cassini, O.C.A.

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

	USE mod_donnees, ONLY : ihe4, nchim, nom_elem
	USE mod_kind

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: xchim	
	REAL (kind=dp) :: t, ro, kappa, dkapdt, dkapdr, dkapdx

	REAL (kind=sp), DIMENSION(19) :: xmet=0.
	REAL (kind=sp) :: opact,dopact,dopacr,dopactd,
	1 fedge,ftredge,fzedge,slr, slt,z_opal,x_opal,xc_opal,xo_opal,
	1 lt_opal,slr_opal,fcn_opal,fcon_opal,fcnone_opal,fu_opal,x,y,z

	INTEGER, PARAMETER :: nmet=19		
	INTEGER :: i
			
c---------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

	CALL set_opal_dir('/home/bilou/SUN_STAR_DATA/')
	CALL read_opal_dump(23,'/home/bilou/OPAL/z14sdtst_GN93hz_withCNO.bin')

	nchim=10 ; ALLOCATE(nom_elem(nchim),xchim(nchim))	
	xchim(1) =3.560E-01 ; nom_elem(1) =' H1 '
	xchim(2) =8.718E-06 ; nom_elem(2) ='He3 ' 
	xchim(3) =6.249E-01 ; nom_elem(3) ='He4 ' 
	xchim(4) =1.804E-05 ; nom_elem(4) ='C12 '
	xchim(5) =4.933E-06 ; nom_elem(5) ='C13 '
	xchim(6) =4.683E-03 ; nom_elem(6) ='N14 '
	xchim(7) =2.067E-07 ; nom_elem(7) ='N15 '
	xchim(8) =8.425E-03 ; nom_elem(8) ='O16 ' 
	xchim(9) =5.373E-04 ; nom_elem(9) ='O17 '
	xchim(10)=5.389E-03 ; nom_elem(10)='Si28'

	ihe4=3 ; y=SUM(xchim(ihe4-1:ihe4)) ; z=SUM(xchim(ihe4+1:nchim))
	x=1.d0-y-z ; WRITE(*,2000)x,y,z
	
	CALL opa_opal2(xchi,t,ro,kap,dkapt,dkapro,dkapx)

c détermination des abondances/gr des éléments de OPAL
	DO i=1,nchim	
	 IF(nom_elem(i)(1:2)=='C1')THEN
	  xmet(1)=xmet(1)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='N1')THEN
	  xmet(2)=xmet(2)+xchim(i)
	 ELSEIF(nom_elem(i)(1:1)=='O')THEN
	  xmet(3)=xmet(3)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Ne')THEN
	  xmet(4)=xmet(4)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Na')THEN
	  xmet(5)=xmet(5)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Mg')THEN
	  xmet(6)=xmet(6)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Al')THEN
	  xmet(7)=xmet(7)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Si')THEN
	  xmet(8)=xmet(7)+xchim(i)
	 ELSEIF(nom_elem(i)(1:1)=='P')THEN
	  xmet(9)=xmet(9)+xchim(i)
	 ELSEIF(nom_elem(i)(1:1)=='S')THEN
	  xmet(10)=xmet(10)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Cl')THEN
	  xmet(11)=xmet(11)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Ar')THEN
	  xmet(12)=xmet(12)+xchim(i)
	 ELSEIF(nom_elem(i)(1:1)=='K')THEN
	  xmet(13)=xmet(13)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Ca')THEN
	  xmet(14)=xmet(14)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Ti')THEN
	  xmet(15)=xmet(15)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Cr')THEN
	  xmet(16)=xmet(16)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Mn')THEN
	  xmet(17)=xmet(17)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Fe')THEN
	  xmet(18)=xmet(18)+xchim(i)
	 ELSEIF(nom_elem(i)(1:2)=='Ni')THEN
	  xmet(19)=xmet(19)+xchim(i)
	 ENDIF
	ENDDO
		  
	WRITE(*,2000) xmet

	t=1.553E+07 ; ro=1.494E+02	
	slt=LOG10(t)-6. ; slr=LOG10(ro)-3.*slt
	CALL opal_x_cno_fu(x,slt,slr,xmet,nmet,1.)	
	CALL ask_last_opac(opact,dopact,dopacr,dopactd,
	1 fedge,ftredge,fzedge)
	CALL ask_last_xcnou( z_opal, x_opal, xc_opal, xo_opal, slt_opal,
	1 slr_opal, fcn_opal, fcon_opal, fcnone_opal, fu_opal )
	
1	FORMAT(' INIT: ztot=',es10.3,' slt,slr=',2es10.3,' nmet=',i4)
2	FORMAT(' xmet(C,N,O,Ne,...)=',4es10.3,3(/6x,5es10.3))
3	FORMAT('opact=log10(Kappa)=',es10.3,/,
	1 'DERIVATIVES dopact,dopacr,dopactd=',3es10.3,/,
	2 'EDGE FACTORS fedge,ftredge,fzedge=',3es10.3)
4	FORMAT('Z=',es10.3,', X=',es10.3,', exC,O=',2es10.3 ,/,'slt,slr=',
	1 2es10.3,/,'fCN,fCON,fCNONe,fu=',4es10.3)
	
	WRITE(*,3) opact,dopact,dopacr,dopactd,fedge,ftredge,fzedge
	WRITE(*,4) z_opal, x_opal, xc_opal, xo_opal, slt_opal,
	1 slr_opal, fcn_opal, fcon_opal, fcnone_opal, fu_opal
	PRINT*,' [note the bad Z value]'
	
	STOP

	END PROGRAM test_opa_opal2