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

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

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

	USE mod_donnees, ONLY: aradia, clight, f_eos, ihe4, lit_nl, m_ch,
	1 nchim, ne, nom_elem, nom_etat,  nom_opa, ord_qs, re_nu, 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_opa, ONLY : opa	
	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 xchim0, 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) :: cte4, p, t, stor0, stor, dstor,
	1 ro, drop, drot, drox, u, dup, dut, dux, dt,
	2 delta, deltap, deltat, deltax, cp, nu, nu_rad, nu_rad0,
	3 dcpp, dcpt, dcpx, gradad, dgradadp, dgradadt,
	4 dgradadx, alfa, beta, gamma1, xs, kappa, dkapdt,dkapdr,dkapdx,
	5 kappa0, dkapdt0, dkapdr0, dkapdx0, dnu_radx
	
	INTEGER :: i, idep, j, k, l=1
		
	LOGICAL :: deriv
	
	CHARACTER (len=80) :: chaine
		
c-----------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

c lecture du fichier	 
	CALL lit_binaire(chaine,dt)
	
c constante pour re_nu
	cte4=re_nu*aradia/clight*4.d0/15.d0
	
c pour initialisation des constantes		
	CALL lit_nl(w_rot)
	
c nom de l'équation d'état à utiliser, 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 nom de l'opacité à tester, s'il diffère de celui du fichier *.don	
c	nom_opa='opa_yveline'
c	nom_opa='opa_yveline_lisse'	
c	nom_opa='opa_gong'
	nom_opa='opa_houdek9'	
	
c allocations et initialisation des nucleo
	ALLOCATE(dfqs(ne),dxchim(nchim), ex(nchim), fqs(ne), xchim(nchim),
	1 xchim0(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) ; xchim0=xchim

c équation d'état	 
	 CALL etat(p,t,xchim0,.TRUE.,
	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)	 

c valeurs initialses
	 deriv=.TRUE.
	 CALL opa(xchim0,t,ro,kappa0,dkapdt0,dkapdr0,dkapdx0)
	 PRINT*,'kappa,dkapdt,dkapdr,dkapdx'	
	 WRITE(*,2000)kappa0,dkapdt0,dkapdr0,dkapdx0

c pour nu_rad, dérivée par GRAMME	 
	 nu_rad0=cte4*t**4/kappa0/ro**2
	 dnu_radx=-nu_rad0*((dkapdx0+dkapdr0*drox)/kappa0+2.d0*drox/ro)
	 	  	 	 
	 deriv=.FALSE.	

	 WRITE(*,*)'------------------------'
	 stor0=t ; stor=stor0*unpdx
	 dstor=stor-stor0
	 t=stor
	 CALL opa(xchim0,t,ro,kappa,dkapdt,dkapdr,dkapdx)
	 t=stor0
	 WRITE(*,*)'dkapdt'
	 WRITE(*,2000)dkapdt0,(kappa-kappa0)/dstor
	 PAUSE'dérivée par rapport à T'
	 
	 stor0=ro ; stor=stor0*unpdx
	 dstor=stor-stor0
	 ro=stor
	 CALL opa(xchim0,t,ro,kappa,dkapdt,dkapdr,dkapdx)
	 ro=stor0
	 WRITE(*,*)'dkapdro'
	 WRITE(*,2000)dkapdr0,(kappa-kappa0)/dstor
	 PAUSE'dérivée par rapport à ro'

c dérivées par GRAMME
c pour le test d kappa / dX, mettre dans la routine d'opacité z=z0
	 xchim=xchim0
	 DO j=1,nchim
	  stor0=xchim(j) ; stor=stor0*unpdx
	  dstor=stor-stor0
	  xchim(j)=stor ; xchim(ihe4)=xchim(ihe4)-dstor
	  CALL opa(xchim,t,ro,kappa,dkapdt,dkapdr,dkapdx)
	  WRITE(*,*)'dkapdx',nom_elem(j)		!,kappa,kappa0,dstor
	  WRITE(*,2000)dkapdx0,(kappa-kappa0)/dstor
	 ENDDO
	 PAUSE'dérivées par rapport à xchim'	 

c dérivée de re_nu /X
c pour le test d kappa / dX, mettre dans la routine d'opacité z=z0
	 xchim=xchim0	
	 stor0=xchim(1) ; stor=stor0*unpdx
	 dstor=stor-stor0
	 xchim(1)=stor ; xchim(ihe4)=xchim(ihe4)-dstor

c équation d'état	 
	 CALL etat(p,t,xchim,.TRUE.,
	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)	 
c opacité
	  CALL opa(xchim,t,ro,kappa,dkapdt,dkapdr,dkapdx)	 
	 WRITE(*,*)'dnu_radX'

c pour nu_rad	 	 
	 nu_rad=cte4*t**4/kappa/ro**2	 
	 WRITE(*,2000)dnu_radx,(nu_rad-nu_rad0)/dstor
	 	 
	ENDDO
	
	STOP

	END PROGRAM test_opa2