c	INCLUDE '../SOURCE/mod_donnees.f'
c	INCLUDE '../SOURCE/mod_nuc.f'
c	INCLUDE '../SOURCE/mod_cesam.f'	

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

	PROGRAM test_jacobien_reac_nuc

c	test pour le jacobien des réactions nucléaires

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

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

	USE mod_donnees, ONLY : ab_min, lit_nl, nchim,
	1 nom_elem, nom_fich2, nucleo, rot_solid, secon6, w_rot
	USE mod_kind
	USE mod_nuc, ONLY : nuc
		
	IMPLICIT NONE

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: jac, jac0    
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: anuc, comp, dcomp, ex, ex0,
	1 dcomp0
	REAL (kind=dp), DIMENSION(5) :: epsilon, epsilon0
	REAL (kind=dp), PARAMETER :: dd=1.d-4, unpdd=1.d0+dd 
	REAL (kind=dp) :: t, ro, et, et0, ero, ero0, hhe, be7e, b8e, n13e,
	1 o15e, f17e, stor, stor0, dstor

	INTEGER :: fait, i, j

	LOGICAL :: deriv
	
	CHARACTER (len=5), ALLOCATABLE, DIMENSION(:) :: nom

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

2000	FORMAT(8es10.3)

c lecture du fichier de données
c	nom_fich2='test_jacob'
	nom_fich2='5msol'
	CALL lit_nl(w_rot)

c	test sur la routine de réactions nucléaires

c	t=2.593d+07 ; ro=5.077d+00
c	t=1.2E+08 ; ro=2.351E+02	
c	t=1.858E+08 ; ro=2.351E+02
c	t=1.491E+07 ; ro=5.881E+01	
	t= 1.589E+08 ; ro= 7.850E+03
	
	ALLOCATE(jac(0,0),comp(0),dcomp(0),ex(0))

	fait=0
	CALL nuc(t,ro,comp,dcomp,jac,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)
	
	DEALLOCATE(jac,comp,dcomp,ex)
	ALLOCATE(anuc(nchim),jac(nchim,nchim),jac0(nchim,nchim),comp(nchim),
	1 dcomp(nchim),dcomp0(nchim),ex(nchim),ex0(nchim))

c calcul des abondances initiales
	fait=1
	CALL nuc(t,ro,comp,dcomp,jac,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)

c	PRINT*,'composition chimique par mole' 
c	WRITE(*,2000)comp(1:nchim)
c	IF(nom_nuc(1:6) /= 'ppcno9')THEN
c	 comp(2)=1.d-18 ; comp(3)=1.d-6
c	 PRINT*,'composition chimique modifiée'   
c	 WRITE(*,2000)comp(1:nchim)    
c	ENDIF

c pour ppcno12 nouvelle composition chimique et t, ro 
c	t= 1.439E+07 ; ro= 5.562E+01
c	comp=(/	6.988E-01, 2.871E-18, 1.287E-05, 6.881E-02, 1.379E-16,
c	1 6.077E-13, 1.444E-04, 4.098E-05, 1.791E-04, 7.429E-09, 6.028E-04,
c	2 2.402E-07, 2.080E-04 /)

c pour	ppcno3a12Ne
	t= 1.589d8 ; ro= 7.850d3
	comp=(/	2.937d-38, 1.000d-50, 3.781d-02, 3.596d-02, 3.250d-06,
	1 4.664d-11, 2.962d-08, 2.446d-02, 1.065d-07, 1.040d-04, 6.789d-04,
	2 6.160d-08 /)
	
c pour	ppcno3a12Ne	
	t= 1.551E+08 ; ro= 8.274E+03	
	comp=(/	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 /)		
	WRITE(*,2000)t,ro ; WRITE(*,2000)comp ;	PAUSE'composition chimique' 

c les tests
	fait=2 ; deriv=.TRUE.
c	fait=3 ; deriv=.TRUE.
	CALL nuc(t,ro,comp,dcomp0,jac0,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)

c	WRITE(*,2000)epsilon(1:4)/secon6
c	WRITE(*,2000)t,ro,epsilon(1:4),et ; PAUSE'epsilon'
	WRITE(*,2000)dcomp0 ; PAUSE'dcomp0'

c conservation des baryons
	PRINT*,'conservation des baryons'	
	anuc=ANINT(nucleo)	
	WRITE(*,2000)DOT_PRODUCT(dcomp0,anuc) ; PAUSE'vérif' 

c	PRINT*,'jac'
c	DO i=1,nchim
c	 PRINT*,i
c	 WRITE(*,2000)jac0(i,1:nchim)
c	ENDDO
c	PAUSE'jac0'

	ALLOCATE(nom(nchim))
	DO i=1,nchim
	 nom(i)=nom_elem(i)//' '
	ENDDO
 	PRINT*,nom ; PAUSE
	
	deriv=.FALSE.
	B1: DO i=1,nchim
	 IF(comp(i) < ab_min(i))CYCLE B1
	 stor0=comp(i) ; stor=stor0*unpdd
	 IF(stor <= 0.d0)stor=dd
	 dstor=stor-stor0 ; comp(i)=stor
	 CALL nuc(t,ro,comp,dcomp,jac,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)
	 PRINT*,(nom_elem(j),' ',j=1,nchim) 
	 PRINT*,'dérivée de dcomp / ',nom_elem(i)
	 WRITE(*,2000)(jac0(j,i),(dcomp(j)-dcomp0(j))/dstor,j=1,nchim)
	 comp(i)=stor0
	 PAUSE 
	ENDDO B1
	
c vérification des dérivées de epsilon/ T, ro, Xi
	fait=3 ; deriv=.TRUE.
	CALL nuc(t,ro,comp,dcomp0,jac,deriv,fait,
	1 epsilon0,et0,ero0,ex0,hhe,be7e,b8e,n13e,o15e,f17e)

	stor0=t ; deriv=.FALSE.     !dérivée/T
	stor=stor0*unpdd ; dstor=stor-stor0 ; t=stor
	CALL nuc(t,ro,comp,dcomp0,jac,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)
	PRINT* ; WRITE(*,2000)et0,(epsilon(1)-epsilon0(1))/dstor
	PAUSE'dérivée epsilon/ T' ; t=stor0

	stor0=ro    !dérivée/ro
	stor=stor0*unpdd ; dstor=stor-stor0 ; ro=stor
	CALL nuc(t,ro,comp,dcomp0,jac,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)
	PRINT* ; WRITE(*,2000)ero0,(epsilon(1)-epsilon0(1))/dstor
	PAUSE'dérivée epsilon /ro' ; ro=stor0

	DO i=1,nchim
	 stor0=comp(i)    !dérivée/Xi
	 stor=stor0*unpdd  
	 IF(stor <= 0.d0)stor=dd
	 dstor=stor-stor0 ; comp(i)=stor
	 CALL nuc(t,ro,comp,dcomp0,jac,deriv,fait,
	1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)
	 PRINT*,'dérivée / ',nom_elem(i) ; comp(i)=stor0  
	 WRITE(*,2000)ex0(i),(epsilon(1)-epsilon0(1))/dstor
	 PAUSE'dérivée epsilon /Xi'
	ENDDO
	
	STOP
    
	END PROGRAM test_jacobien_reac_nuc