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

	program test_diffusion
	
c	test pour les coefficients de diffusion
c	CESAM95

c	Auteur: P. Morel, Departement J.D. Cassini, O.C.A.
	
c---------------------------------------------------------------------

	use mod_chim
	use mod_donnees
	use mod_exploit
	use mod_kind

	implicit none
	
	integer, parameter :: pnosc=800, pnelem=15
	
	real (kind=dp), allocatable, dimension(:,:,:) :: dd
	real (kind=dp), allocatable, dimension(:,:) :: var, y, d, d0, dv	
	real (kind=dp), allocatable, dimension(:) :: r, p, t, ro, l, m,
	1	xchim, v, v0, glob, kap, gradrad, gradad, gamma1, cp, delta	
	real (kind=dp), dimension(3) :: m_zc=0.d0, r_zc=0.d0
	real (kind=dp), parameter :: dx=1.d-5, unpdx=1.d0+dx
	real (kind=dp) :: mstar, stor, stor0, dstor, age, cte8,
	1	drox=0.d0, dkapx=0.d0, dgradradx=0.d0, dgradadx=0.d0

	real (kind=sp), allocatable, dimension(:) :: xdes, dh, dhe4, dli7,
	1	vh, vhe4, vli7	
	real (kind=sp) :: xmax, xmin, ymax, ymin
	
	integer :: iconst, ivar, i, itot, j, k, ll, ides, ili7=-1, lim=0	
	
	logical :: ok, dessin, ps, melange, deriv
	
	character (len=4), allocatable, dimension(:) :: nom_elem_lu		
	character (len=9) :: device, titre='diffusion'	
	character (len=80) :: modele,nom_fich
	
c---------------------------------------------------------------------

c	itot: nombre de points
c	iconst: nombre de constantes global
c	ivar: nombre de variables du tableau var sans les elements chimiques
c	modele: nom du modele
c	nom_elem: nom des elements utilises
c	nbelem: nombre d'elements utilises
c	nchim: nombre d'elements chimiques utilises	
c	glob: variables globales
c		glob(1)=mstar*msol
c		glob(2)=rstar*rsol
c		glob(3)=ltot*lsol
c		glob(4)=z
c		glob(5)=x0
c		glob(6)=alpha
c		glob(7)=9./4.
c		glob(8)=1./162.
c		glob(9)=1.
c		glob(10)=1.
c		glob(11)=d2p
c		glob(12)=d2ro
c		glob(13)=age
c		glob(14)=vsal
c		glob(15)=0.
c	var: variables
c		var(1,i)=r*rsol
c	 	var(2,i)=log(m) -1.d38 au centre
c		var(3,i)=t
c		var(4,i)=p
c		var(5,i)=ro
c		var(6,i)=gradient
c		var(7,i)=l
c		var(8,i)=kap
c		var(9,i)=energie thermo+gravifique
c		var(10,i)=grand Gamma1
c		var(11,i)=gradient adiabatique
c		var(12,i)=delta
c		var(13,i)=cp
c		var(14,i)=1 / (mu elec.)
c		var(15,i)=vaissala, 0 au centre
c	 	var(16,i)=w(i)
c	 	var(17,i)=d ln kappa / d ln T
c	 	var(18,i)=d ln kappa / d ln ro
c	 	var(19,i)=d epsilon(nuc) / d ln T
c	 	var(20,i)=d epsilon(nuc) / d ln ro
c		var(21,i)=Ptot/Pgas
c		var(22,i)=gradrad
c	  	var(22+j,i)=xchim(j), j=1,nbelem
		
2000	format(1x,8es10.3)
2001	format(1x,10es8.1)
	
c	lecture de test.don

	nom_fich2='test' ; call lit_nl ; call ini_ctes_phys ; call print_ctes(6)	

	cte8=3.d0/16.d0/pi/aradia/clight/g*lsol/msol
	
	nom_fich='evry'
	
	allocate(glob(15),var(22+pnelem,pnosc),nom_elem_lu(pnelem))
	
	call read_osc(nom_fich,itot,iconst,ivar,modele,glob,var,
	1	nbelem,nchim,iw,nom_elem_lu,ok)	
		
	mstar=glob(1)/msol
	allocate(r(itot),p(itot),t(itot),ro(itot),l(itot),m(itot),
	1	kap(itot),y(nbelem,itot),gradrad(itot),gradad(itot),
	2	xdes(itot),dh(itot),dhe4(itot),dli7(itot),vh(itot),
	3	vhe4(itot),vli7(itot),gamma1(itot),cp(itot),delta(itot),
	4	nom_elem(nbelem),nucleo(nbelem),zi(nbelem))
	m=0.d0
	if(.not.ok)stop
	do i=1,itot	 
	 if(var(2,i) > -1.d30)m(i)=exp(var(2,i))
	 p(i)=var(4,i)
	 t(i)=var(3,i)	 
 	 ro(i)=var(5,i)
	 l(i)=var(7,i)/lsol
	 r(i)=var(1,i)/rsol
	 kap(i)=var(8,i)
	 gradrad(i)=var(22,i)
	 gradad(i)=var(11,i)
	 gamma1(i)=var(10,i)
	 cp(i)=var(13,i)
	 delta(i)=var(12,i)
	 do j=1,nbelem
	  y(j,i)=max(var(ivar+j,i),1.d-30)	 
	 enddo
	enddo
		
	nom_elem=nom_elem_lu(1:nbelem)	
	deallocate(var,glob,nom_elem_lu)
	
	do i=1,nbelem
	 if(nom_elem(i) == ' H1 ')then
	  nucleo(i)=ah
	  zi(i)=1
	 elseif(nom_elem(i) == ' H2 ')then
	  nucleo(i)=ah2
	  zi(i)=1
	 elseif(nom_elem(i) == 'He3 ')then
	  nucleo(i)=ahe3
	  zi(i)=2
	 elseif(nom_elem(i) == 'He4 ')then
	  nucleo(i)=ahe4
	  zi(i)=2
	  ihe4=i
	 elseif(nom_elem(i) == 'Li7 ')then
	  nucleo(i)=ali7
	  zi(i)=3
	  ili7=i
	 elseif(nom_elem(i) == 'Be7 ')then
	  nucleo(i)=abe7
	  zi(i)=4
	 elseif(nom_elem(i) == 'C12 ')then
	  nucleo(i)=ac12
	  zi(i)=6
	 elseif(nom_elem(i) == 'C13 ')then
	  nucleo(i)=ac13
	  zi(i)=6
	 elseif(nom_elem(i) == 'N14 ')then
	  nucleo(i)=an14
	  zi(i)=7
	 elseif(nom_elem(i) == 'N15 ')then
	  nucleo(i)=an15
	  zi(i)=7
	 elseif(nom_elem(i) == 'O16 ')then
	  nucleo(i)=ao16
	  zi(i)=8
	 elseif(nom_elem(i) == 'O17 ')then
	  nucleo(i)=ao17
	  zi(i)=8
	 elseif(nom_elem(i) == ' Ex ')then
	  nucleo(i)=28
	  zi(i)=13
	 else
	  print*,'erreur de nom d''element',nom_elem(i)
	  stop
	 endif
	enddo
	
	allocate(xchim(nbelem),v0(nbelem),dv(nbelem,nbelem),d0(nbelem,nbelem),
	1	dd(nbelem,nbelem,nbelem),v(nbelem),d(nbelem,nbelem))
	
	age=0.d0
	ides=0

c	le test

	b1: do ll=1,itot	
c	 if(t(ll) < 3.872E+06)cycle b1
c	 if(t(ll) > 3.933E+06)cycle b1
c	 if(t(ll) < 1.d5)cycle b1
	 if(t(ll) < 1.d6)cycle b1
	 melange=gradrad(ll) > gradad(ll)
	 if(melange)cycle b1
	 if(m(ll)*r(ll) <= 0.d0)cycle b1
	 xchim(:)=y(:,ll)/nucleo(:)	!comp.chim. par mole
	 
	 write(*,2000)t(ll),r(ll) ; write(*,2001)xchim	 
	 
	 v0=0.d0 ; d0=0.d0 ; dv=0.d0 ; dd=0.d0	 
	 call diffm(p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox,kap(ll),dkapx,
	1	gradrad(ll),dgradradx,gradad(ll),dgradadx,mstar,xchim,d0,dd,v0,dv)
	 print*,'V ',nom_diffm
	 write(*,2001)(v0(i),i=1,nbelem)
	 print*,'D ',nom_diffm
	 do i=1,nbelem
	  write(*,2001)(d0(i,j),j=1,nbelem)
	 enddo
	 
c	 v0=0.d0 ; d0=0.d0 ; dv=0.d0 ; dd=0.d0	 
	 call difft(melange,p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox,
	1	kap(ll),dkapx,age,mstar,gradrad(ll),dgradradx,gradad(ll),
	2	dgradadx,m_zc,r_zc,lim,gamma1(ll),cp(ll),delta(ll),xchim,
	3	d0,dd,v0,dv)
	 print*,'melange=',melange
	 print*,'V ',nom_difft
	 write(*,2001)(v0(i),i=1,nbelem)
	 print*,'D ',nom_difft
	 do i=1,nbelem
	  write(*,2001)(d0(i,j),j=1,nbelem)
	 enddo
	 	 
	 pause'dans coeff_diff'
	  
c	 do k=1,nbelem
c	  print*,'dD / Xk, k=',k	  
c	  do i=1,nbelem
c	   write(*,2000)(dd(i,j,k),j=1,nbelem)
c	  enddo
c	 enddo
c	 pause
	
	 deriv=.false.
c	 deriv=.true.
	 if(deriv)then
	  write(*,2000)r(ll),t(ll)
	  do i=1,nchim	  
	   v=0.d0 ; d=0.d0 ; dv=0.d0 ; dd=0.d0
	   stor0=xchim(i) ; stor=stor0*unpdx ; dstor=stor-stor0 ; xchim(i)=stor
	 
c	   call diffm(p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox,kap(ll),dkapx,
c	1	gradrad(ll),dgradradx,gradad(ll),dgradadx,mstar,xchim,d,dd,v,dv)

	   call difft(melange,p(ll),t(ll),r(ll),l(ll),m(ll),ro(ll),drox,
	1	kap(ll),dkapx,age,mstar,gradrad(ll),dgradradx,gradad(ll),
	2	dgradadx,m_zc,r_zc,lim,gamma1(ll),cp(ll),delta(ll),xchim,
	3	d,dd,v,dv)
	  
	   xchim(i)=stor0
	   write(*,"('derivee / ',a4)")nom_elem(i)
c	   write(*,2000)v0 ; write(*,2000)v
	   write(*,2000)((v(j)-v0(j))/dstor,dv(j,i),j=1,nchim)
	  enddo
	  pause
	 endif	  	
	
	 ides=ides+1
	 xdes(ides)=r(ll)
	 dh(ides)=d0(1,1) ; dhe4(ides)=d0(ihe4,ihe4)
	 if(ili7 > 0)dli7(ides)=d0(ili7,ili7)	 
	 vh(ides)=v0(1)   ; vhe4(ides)=abs(v0(ihe4))
	 if(ili7 > 0)vli7(ides)=abs(v0(ili7))	 
	 	 	 
	enddo b1

	dessin=.false. ; dessin=.true.
	ps=.false.
	if(dessin)then
	
	 xmax=maxval(xdes)*1.05
	 xmin=0.
	 
	 if(ili7 > 0)then
	  ymax=max(maxval(dh),maxval(dhe4),maxval(dli7))*1.05
	 else
	  ymax=max(maxval(dh),maxval(dhe4))*1.05
	 endif
	 ymin=-ymax/10.
	
	 if(.not.ps)then
	  device='/xw'
	 else
	  device='/ps'
	 endif	

c	 write(*,*)'device ? /xw, /PS, /CPS'
c	 read(5,'(a)')device
	 call pgbegin(0,device,1,2)
	 call pgenv(xmin,xmax,ymin,ymax,0,0)
	 call pglabel('R/Rsol','Diff',titre)

	 call pgsls(1)
	 call pgline(ides,xdes,dh)
	 call pgsls(2)
	 call pgline(ides,xdes,dhe4)
	 call pgsls(3)
	 if(ili7 > 0)call pgline(ides,xdes,dli7)
	 
	 call pgsls(1)
	 if(ili7 > 0)then
	  ymax=max(maxval(vh),maxval(vhe4),maxval(vli7))*1.05
	 else
	  ymax=max(maxval(vh),maxval(vhe4))*1.05
	 endif
	 ymin=-ymax/10.
	 call pgenv(xmin,xmax,ymin,ymax,0,0)
	 call pglabel('R/Rsol','Vdiff',titre)
	 call pgsls(1)
	 call pgline(ides,xdes,vh)
	 call pgsls(2)
	 call pgline(ides,xdes,vhe4)
	 call pgsls(3)
	 if(ili7 > 0)call pgline(ides,xdes,vli7) 	 
	
	 call pgend
	 
	endif
	
	stop
	
	end program test_diffusion