INCLUDE '../SOURCE/mod_opa.f'

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

	PROGRAM test_opa
	
c	test pour les routines d'opacité
c	le nom de la routine d'opacité et du fichier à utiliser sont lus
c	dans le fichier mon_modele.don qu'il convient d'adapter

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

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

	USE mod_donnees, ONLY: ab_min, ihe4, lit_nl, nchim,
	1 nom_elem, nom_fich2, nom_xheavy, nucleo, print_ctes, w_rot,
	2 x0, y0, z0	
	USE mod_kind
	USE mod_nuc, ONLY: nuc
	USE mod_opa, ONLY: opa
		
	IMPLICIT NONE

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: jac
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: xchim, dxchim, ex
	REAL (kind=dp), DIMENSION(5) :: epsilo

	REAL (kind=dp), PARAMETER :: dd=1.d-6, unpdd=1.d0+dd
	REAL (kind=dp) :: be7, b8, et, ero, f17, hh, n13, o15
	REAL (kind=dp) :: t, ro, kappa, dkapdt, dkapdr, dkapdx
	REAL (kind=dp) :: kappa0, dkapdt0, dkapdr0, dkapdx0
	
	REAL (kind=dp) :: dstor,stor,stor0
	INTEGER :: ktest	

	LOGICAL :: deriv
		
c---------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

c lecture des données
	nom_fich2='test'
c	nom_fich2='5msol'	
	CALL lit_nl(w_rot)

c appel d'initialisation pour tabulation des réactions nucléaires
c allocations fictives
	ALLOCATE(xchim(0),dxchim(0),jac(0,0),ex(0))
	CALL nuc(1.5d+07,1.5d+02,xchim,dxchim,jac,.FALSE.,0,
	1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17)

c détermination des abondances initiales
	DEALLOCATE(xchim,dxchim,ex)
	ALLOCATE(xchim(nchim),dxchim(nchim),ex(nchim))
	CALL nuc(1.5d+07,1.5d+02,xchim,dxchim,jac,.FALSE.,1,
	1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17)
	DEALLOCATE(dxchim,ex,jac)
	
c pour	ppcno3a12Ne	
	t= 1.551E+08 ; ro= 8.274E+03	
	xchim=(/ 0.d0, 0.d0, 4.507d-02, 4.377d-02, 3.190d-06,
	1 9.329d-09, 2.972d-08, 1.679d-02, 1.153d-07, 1.040d-04, 6.783d-04,
	2 2.198d-07 /)		

c par gramme	
	xchim=xchim*nucleo
	ktest=6
	B2: DO
	 SELECT CASE(ktest)
	 CASE(1)
	  t=6.554E+03 ; ro=4.149E-07 ; WRITE(*,1)t,ro
1	  FORMAT('t=',es10.3,', ro=',es10.3)	  
	 CASE(2)	  
	  t=4.d3 ; ro=1.d-10 ; WRITE(*,1)t,ro
	 CASE(3)	  
	  t=9.0380E+03 ; ro=2.1197E-08 ; WRITE(*,1)t,ro
	 CASE(4)	  
	  t=8.180E+03 ; ro=4.364E-07 ; WRITE(*,1)t,ro
	 CASE(5)	  
	  t=2.119E+07 ; ro=9.021E+02 ; WRITE(*,1)t,ro
	 CASE(0)	  
	  t=3.026E+06 ; ro=6.272E-01 ; WRITE(*,1)t,ro	!opa=1.452E+01	  	  
	 CASE(6)	  
	  t=7.937E+08 ; ro=7.381E+05 ; WRITE(*,1)t,ro
	  
	 CASE(7)	  
c pour	ppcno3a12Ne	
	  t= 1.9E+08 ; ro= 8.274E+03	
	  xchim=(/ 0.d0, 0.d0, 4.507d-02, 4.377d-02, 3.190d-06,
	1 9.329d-09, 2.972d-08, 1.679d-02, 1.153d-07, 1.040d-04, 6.783d-04,
	2 2.198d-07 /)		

c par gramme	
	  xchim=xchim*nucleo
	  WRITE(*,2000)t,ro ; WRITE(*,2000)xchim
	  PAUSE'composition chimique modifiée' 	  

c	 CASE()	  
c	  t= ; ro= ; WRITE(*,1)t,ro
	 CASE DEFAULT	 
	  EXIT B2
	 END SELECT
	 ktest=ktest+1 ; deriv=.TRUE.	
	 CALL opa(xchim,t,ro,kappa0,dkapdt0,dkapdr0,dkapdx0)
	 PRINT*,'kappa,dkapdt,dkapdr,dkapdx'	
	 WRITE(*,2000)kappa0,dkapdt0,dkapdr0,dkapdx0
	 deriv=.FALSE.	

	 WRITE(*,*)
	 stor0=t ; stor=stor0*unpdd
	 IF(stor == 0.d0)stor=dd
	 dstor=stor-stor0
	 t=stor
	 CALL opa(xchim,t,ro,kappa,dkapdt,dkapdr,dkapdx)
	 t=stor0
	 WRITE(*,*)'dkapdt'
	 WRITE(*,2000)dkapdt0,(kappa-kappa0)/dstor
	 PAUSE'dérivée par rapport à T'
	 
	 WRITE(*,*)
	 stor0=ro ; stor=stor0*unpdd
	 IF(stor == 0.d0)stor=dd
	 dstor=stor-stor0
	 ro=stor
	 CALL opa(xchim,t,ro,kappa,dkapdt,dkapdr,dkapdx)
	 ro=stor0
	 WRITE(*,*)'dkapdr'
	 WRITE(*,2000)dkapdr0,(kappa-kappa0)/dstor
	 PAUSE'dérivée par rapport à ro'

c pour le test d kappa / dX, mettre dans la routine d'opacité z=z0	 
	 WRITE(*,*)
	 stor0=xchim(1) ; stor=stor0*unpdd
	 IF(stor == 0.d0)stor=dd
	 dstor=stor-stor0
	 xchim(1)=stor
	 CALL opa(xchim,t,ro,kappa,dkapdt,dkapdr,dkapdx)
	 xchim(1)=stor0
	 WRITE(*,*)'dkapdx'		!,kappa,kappa0,dstor
	 WRITE(*,2000)dkapdx0,(kappa-kappa0)/dstor
	 PAUSE'dérivée par rapport à X'	  	 
	ENDDO B2
	 
	STOP
	
	END PROGRAM test_opa