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

	PROGRAM test_difft_smc
	
c test des dérivées de la routine difft_smc

c Auteur: P.Morel, Département Cassiopée, O.C.A., CESAM2k

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

	USE mod_donnees, ONLY : d_turb, ihe4, lit_nl, mu_saha, m_ch, nchim, ne,
	1 nom_elem, nom_etat, nom_opa, nucleo, ord_qs, rsol, w_rot
	USE mod_etat, ONLY :  etat, saha
	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), ALLOCATABLE, DIMENSION(:,:,:) :: dd, dd0
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: 	d, d0, y, y0	
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: dfqs, dxchim, dxchimg,
	1 dxchim0, dxchimg0, fqs, xchim, xchimg, xchim0, xchimg0
	
	REAL (kind=dp), DIMENSION(0,0) :: jac(0,0)
	REAL (kind=dp), DIMENSION(0) :: eps(0), ex(0)	

	REAL (kind=dp), PARAMETER :: dx=1.d-6, unpdx=1.d0+dx
	REAL (kind=dp) :: alfa, beta, be7e, b8e, cp, cte1, dcpp, dcpt, dcpx,
	1 ddx, delta, deff, deltap, deltat, deltax, dlnmu, dlnmu0,
	2 drop, drot, drox, dt, dup, dut, dux, et, ero, f17e, gradad,
	3 dgradadp, dgradadt, dgradadx, grad_mu, grad_mu0, gamma1, hhe, lum,
	4 hp, m, mu, mu0, nu, nel, n13e, o15e, p, r, ro, ro0, t, u, Zbar
	
	INTEGER :: i, idep, k, kp, l=1
	
	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 nom de l'opacité à tester, s'il diffère de celui du fichier *.don	
	nom_opa='opa_yveline'
c	nom_opa='opa_yveline_lisse'	
c	nom_opa='opa_gong'
c	nom_opa='opa_houdek9'	

c initialisations
	deff=0.d0
	
c allocations	 
	ALLOCATE(d(nchim,nchim), dd(nchim,nchim,nchim), dd0(nchim,nchim,nchim),
	1 dfqs(ne),
	2 dxchim(nchim), dxchimg(nchim), dxchim0(nchim), dxchimg0(nchim),
	3 d0(nchim,nchim), fqs(ne), xchim(nchim), xchimg(nchim), xchim0(nchim),
	4 xchimg0(nchim), y(nchim,0:1),y0(nchim,0:1))		

c initialisation des nucleo
	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,hhe,be7e,b8e,n13e,o15e,f17e)		
	CALL nuc(2.593d+07,5.077d+00,xchim,dxchim,jac,.FALSE.,1,
	1 eps,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e)
	
	mu_saha=.FALSE.
c	mu_saha=.TRUE.

c couche à couche
c	idep=2 ; kp=15
	idep=163 ; kp=15	
c	idep=950 ; kp=15	
	
	DO i=idep, idep+kp
	 PRINT*,'i=',i
	 
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)) ; r=SQRT(fqs(3))
	 lum=SQRT(fqs(4))**3 ; nu=fqs(5) ; m=SQRT(nu)**3 

c la composition chimique	 
	 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,nu,l,xchim,
	1 dxchim)
c	 WRITE(*,2000)xchim(1:4)
	 y0(:,0)=xchim ; y0(:,1)=dxchim ; d0=0.d0 ; dd0=0.d0
	 CALL difft_smc(.TRUE.,deff,lum,m,p,r,t,y0,d0,dd0)
	 	 	
	 B1: DO k=1,1!nchim
	 
c	  EXIT B1
	  
	  d=0.d0 ; dd=0.d0
	  y=y0 ; y(k,0)=y(k,0)*unpdx ; ddx=y(k,0)-y0(k,0)
	  SELECT CASE(k)
	  CASE(1)
	   y(ihe4,0)=y(ihe4,0)-ddx*nucleo(1)/nucleo(ihe4)
	  CASE DEFAULT
	   y(1,0)=y(1,0)-ddx*nucleo(k)/nucleo(1)
	  END SELECT	   
	  CALL difft_smc(.TRUE.,deff,lum,m,p,r,t,y,d,dd)	  
	  PRINT*,'dérivées/',nom_elem(k)
	  WRITE(*,2000)(d(k,k)-d0(k,k))/ddx,dd0(k,k,1:3)
	 ENDDO B1
	ENDDO
	  		
	STOP	
	
	CONTAINS
	
	INCLUDE '../SOURCE/difft_smc.f'	
	

	END PROGRAM test_difft_smc