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

	PROGRAM test_mu_mol
	
c test des dérivées de la routine mu_mol

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

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

	USE mod_donnees, ONLY : lit_nl, mu_saha, m_ch, nchim, ne, nom_elem,
	1 ord_qs, rsol, w_rot
	USE mod_etat, ONLY :  etat, saha	!, mu_mol
	USE mod_exploit, ONLY : lit_binaire
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn
	USE mod_nuc, ONLY : nuc	
	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(:) :: dfqs, dxchim, dxchimg,
	1 dxchim0, dxchimg0, dlnmu_x, fqs, dgrad_mux, dmu_x, xchim, xchimg,
	2 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-5, unpdx=1.d0+dx
	REAL (kind=dp) :: alfa, beta, be7e, b8e, cp, cte1, dcpp, dcpt, dcpx,
	1 ddx, delta, 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,
	4 hp, mu, mu0, nu, nel, n13e, o15e, p, r, ro, ro0, t, u, Zbar
	
	INTEGER :: i, idep, iplus, j, k, 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 constantes
	 cte1=rsol/2.d0
	
c allocations	 
	ALLOCATE(dfqs(ne), dlnmu_x(nchim), dgrad_mux(nchim), dmu_x(nchim),
	1 dxchim(nchim), dxchimg(nchim), dxchim0(nchim), dxchimg0(nchim),
	2 fqs(ne), xchim(nchim), xchimg(nchim), xchim0(nchim), xchimg0(nchim))		

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
	idep=3 ; iplus=0
	B1: DO i=idep,idep+iplus
	 PRINT*
		
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)) ; nu=fqs(5) 

c échelle de hauteur de pression : -dR / d ln P	
	 hp=-cte1*dfqs(3)/dfqs(1)/r

c	 PRINT*,'p,t,r,nu,hp,dfqs(3),dfqs(1)'	 
c	 WRITE(*,2000)p,t,r,nu,hp,dfqs(3),dfqs(1)	

c la composition chimique	 
	 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,nu,l,xchim0,
	1 dxchim0)
c	 PRINT*,'xchim,dxchim' ; WRITE(*,2000)xchim0 ; WRITE(*,2000)dxchim0


	DO k=1,nchim
	 DO j=1,2
	  IF(j == 1)THEN
	   xchim=xchim0 ; dxchim=dxchim0
	  ELSE
	   xchim(k)=xchim0(k)*unpdx ; dxchim(k)=dxchim0(k)!*unpdx
	   ddx=xchim(k)-xchim0(k)   
	  ENDIF
		
c composition chimique par gramme, équation d'état
	  xchimg=ABS(xchim) ; dxchimg=dxchim
	  CALL chim_gram(xchimg,dxchimg)
	  CALL etat(p,t,xchimg,.FALSE.,ro,drop,drot,drox,u,dup,dut,dux,
	1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
	
c ro*hp est inchangé par un dX	
	  IF(j == 1)ro0=ro

c	  PRINT*,'ro' ;  WRITE(*,2000)ro	

c appel à mu_mol
	  CALL mu_mol(dxchim,hp,nu,r,ro0,t,xchim,dlnmu,dlnmu_x,grad_mu,
	1 dgrad_mux,mu,dmu_x,nel,Zbar)
	
c	  PRINT*,'dlnmu,dlnmu_x,grad_mu,dgrad_mux,mu,dmu_x,nel,Zbar'
c	  WRITE(*,2000)dlnmu,dlnmu_x,grad_mu,dgrad_mux,mu,dmu_x,nel,Zbar
	  
	  IF(j == 1)THEN
	   dlnmu0=dlnmu ; grad_mu0=grad_mu ; mu0=mu
	  ELSE
	   PRINT*,'dérivées/',nom_elem(k),', dlnmu, grad_mu, mu'
	   WRITE(*,2000)(dlnmu-dlnmu0)/ddx,dlnmu_x(k),
	1  (grad_mu-grad_mu0)/ddx,dgrad_mux(k), (mu-mu0)/ddx,dmu_x(k)
	  ENDIF
	 ENDDO
	 
	ENDDO
	ENDDO B1
			
	STOP
	
	CONTAINS
	
	INCLUDE '../SOURCE/mu_mol.f'
	
	END PROGRAM test_mu_mol