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

	PROGRAM test_etat2
	
	
c test des dérivées de l'équation d'état utilisant un modèle en binaire
c Auteur: P. Morel, Département J.D. Cassini, O.C.A.

c EFF ne donne pas toujours des valeurs exactes pour
c drox, dux, drotx et dutx car X et Y sont reliés par X+Y+Z=1
c ce dont on ne semble pas tenir compte 

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

	USE mod_donnees, ONLY: f_eos, ihe4, lit_nl, m_ch, nchim, ne, nom_elem,
	1 nom_etat, ord_qs, w_rot
	USE mod_etat, ONLY: etat
	USE mod_exploit, ONLY : lit_binaire
	USE mod_kind
	USE mod_nuc, ONLY: nuc
	USE mod_numerique, ONLY : bsp1dn	
	USE mod_variables, ONLY : bp, chim, chim_gram, knot, knotc, mc, mct,
	1 n_ch, n_qs, q, qt
			
	IMPLICIT NONE
	
	REAL (kind=dp), DIMENSION(0,0) :: jac
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: dfqs, dxchim, fqs, xchim,
	1 ex
	REAL (kind=dp), DIMENSION(5) :: eps

	REAL (kind=dp), PARAMETER :: dx=1.d-6, unpdx=1.d0+dx
	REAL (kind=dp) :: be7, b8, et, ero, f17, hh, n13, o15

	REAL (kind=dp) :: p, t, stor0, stor, dstor,
	1 ro, drop, drot, drox, u, dup, dut, dux, dt,
	2 delta, deltap, deltat, deltax, cp, nu,
	3 dcpp, dcpt, dcpx, gradad, dgradadp, dgradadt,
	4 dgradadx, alfa, beta, gamma1, xs
	REAL (kind=dp) :: ro0, drop0, drot0, drox0, u0, dup0, dut0, dux0,
	2 delta0, deltap0, deltat0, deltax0, cp0,
	3 dcpp0, dcpt0, dcpx0, gradad0, dgradadp0, dgradadt0,
	4 dgradadx0, alfa0, beta0, gamma10
	
	INTEGER :: i, idep, k, l=1
		
	LOGICAL :: deriv
	
	CHARACTER (len=80) :: chaine
		
c-----------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

c lecture du fichier	 
	CALL lit_binaire(chaine,dt)
	
c pour initialisation des constantes		
	CALL lit_nl(w_rot)
	
c nom de l'équation d'état à tester, s'il diffère de celui du fichier *.don
c	nom_etat='etat_gong2'
	nom_etat='etat_eff'
c	nom_etat='etat_ceff'
c	nom_etat='etat_opal'
c	nom_etat='etat_opalX'		
	
c allocations et initialisation des nucleo
	ALLOCATE(dfqs(ne),dxchim(nchim), ex(nchim), fqs(ne), xchim(nchim))
	IF(ALLOCATED(nom_elem))DEALLOCATE(nom_elem)
	CALL nuc(2.593d+07,5.077d+00,xchim,dxchim,jac,.FALSE.,0,
	1 eps,et,ero,ex,hh,be7,b8,n13,o15,f17)		
	CALL nuc(2.593d+07,5.077d+00,xchim,dxchim,jac,.FALSE.,1,
	1 eps,et,ero,ex,hh,be7,b8,n13,o15,f17)

c couche à couche
	idep=300 ; k=0
	DO i=idep, idep+k
	 
c variables quasi-statiques	 
	 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),l,fqs,dfqs)
	 p=EXP(fqs(1)) ; t=EXP(fqs(2)) ; nu=fqs(5)

c la composition chimique	 
	 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,nu,l,xchim,
	1 dxchim)
	 
c composition chimique par gramme
	 CALL chim_gram(xchim,dxchim)

c valeurs initialses
	 deriv=.TRUE.
	 CALL etat(p,t,xchim,deriv,
	1 ro0,drop0,drot0,drox0,u0,dup0,dut0,dux0,
	2 delta0,deltap0,deltat0,deltax0,cp0,dcpp0,dcpt0,dcpx0,
	3 gradad0,dgradadp0,dgradadt0,dgradadx0,alfa0,beta0,gamma10)
	 WRITE(*,*)'p, t, X'
	 WRITE(*,2000)p,t,xchim(1)
	 WRITE(*,*)'ro, u, delta, cp, gradad, alfa, beta, gamma1'	 
	 WRITE(*,2000)ro0,u0,delta0,cp0,gradad0,alfa0,beta0,gamma10

	 deriv=.FALSE. ; WRITE(*,*)'--------------------'	
	 WRITE(*,*)
	 WRITE(*,3)p
3	 FORMAT('dérivées par rapport à P=',es10.3)
	 stor0=p ; stor=stor0*unpdx
	 IF(stor == 0.d0)stor=dx
	 dstor=stor-stor0
	 p=stor	
	 CALL etat(p,t,xchim,deriv,
	1 ro,drop,drot,drox,u,dup,dut,dux,
	2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
	 p=stor0
	 WRITE(*,*)'drop,deltap,dgradadp'
	 WRITE(*,2000)drop0,(ro-ro0)/dstor,deltap0,(delta-delta0)/dstor,
	1 dgradadp0,(gradad-gradad0)/dstor
	 WRITE(*,*)'dup,dcpp'
	 WRITE(*,2000)dup0,(u-u0)/dstor,dcpp0,(cp-cp0)/dstor
	 PAUSE'fin dérivées par rapport à P'
	 WRITE(*,4)t
4	 FORMAT('dérivées par rapport à T=',es10.3)
	 stor0=t ; stor=stor0*unpdx
	 IF(stor == 0.d0)stor=dx
	 dstor=stor-stor0
	 t=stor
	 CALL etat(p,t,xchim,deriv,
	1 ro,drop,drot,drox,u,dup,dut,dux,
	2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
	 t=stor0
	 WRITE(*,*)'drot,deltat,dgradadt'
	 WRITE(*,2000)drot0,(ro-ro0)/dstor,deltat0,(delta-delta0)/dstor,
	1 dgradadt0,(gradad-gradad0)/dstor
	 WRITE(*,*)'dut,dcpt'
	 WRITE(*,2000)dut0,(u-u0)/dstor,dcpt0,(cp-cp0)/dstor
	 PAUSE'fin dérivées par rapport à T'
	 	 
	 xs=xchim(1)
	 WRITE(*,5)xchim(1)
5	 FORMAT('dérivées par rapport à X=',es10.3)
	 stor0=xchim(1)
	 stor=stor0*unpdx
	 IF(stor == 0.d0)stor=dx
	 dstor=stor-stor0	 
	 xchim(1)=stor ; xchim(ihe4)=xchim(ihe4)-dstor
	 CALL etat(p,t,xchim,deriv,
	1 ro,drop,drot,drox,u,dup,dut,dux,
	2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
	 xchim(1)=stor
	 WRITE(*,*)'drox,deltax,dgradadx'
	 WRITE(*,2000)drox0,(ro-ro0)/dstor,deltax0,(delta-delta0)/dstor,
	1 dgradadx0,(gradad-gradad0)/dstor
	 WRITE(*,*)'dux,dcpx'
	 WRITE(*,2000)dux0,(u-u0)/dstor,dcpx0,(cp-cp0)/dstor
	 PAUSE'fin dérivées par rapport à X'		
	ENDDO
	
	STOP

	END PROGRAM test_etat2