include'mod_chim'

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

	program test_acc_rad

c	lecture du fichier pour G. Alecian pour test de calcul des forces radiatives
c	Auteur: P. Morel, Departement J.D. Cassini, O.C.A., Observatoire de Nice
c	CESAM
c       Adaptation pour tests des modules g_rad: G.Alecian (7/9/99->...)

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

	use mod_cesam
	use mod_chim
	use mod_donnees
	use mod_kind

	implicit none

	integer	itot,i,j,k,nb
	
	real (kind=dp), dimension(:,:), allocatable :: dg_rad, ioni
	real (kind=dp), dimension(:), allocatable :: g_rad, g_rad0, xi
	real (kind=dp), dimension(pnchim) :: ychim
	real (kind=dp), parameter :: dx=1.d-5, unpdx=1.d0+dx, dkapx=0.d0
	real (kind=dp) :: rtot, ltot, teff, mass, r, ray, t,ro, kap, grad,
	1	lum, nel, mstar, grav, dstor, stor, stor0

	character*1 :: bidon

*====================================================== fin
*** locales reservees a ce main.f de test
	integer iunit
	real*4 logT, logPop(0:pzi)
	character*3	massest
	character*64 NFIa,NFib
*****************************************

2000	format(1x,8es10.3)
2001	format(1x,15es10.3)
2002	format(1x,i5,1p,30('	',es10.3))
2003	format(1x,'i	ychim	logT',30('	',a5,i3))
2004	format(1x,8es15.8)
2005	format(6x,1p,30('	',es10.3))

c	print*, 'Version du 24/11/99'
c	print*, 'Quelle masse pour MOREL_Fe56_ (1Msol=100):'
c	read(5,'(a3)') massest
	massest='100'
	NFib='~/CESAM4/F_RAD/MOREL_Fe56_'//massest//'.data'
	open(unit=1,form='formatted',status='old',file=NFib)

	OPEN(UNIT=10,FILE='grcesam.job',STATUS='unknown',CARRIAGECONTROL='LIST')

	read(1,'(a1)') bidon
	read(1,90)nchim,(nom_elem(i),i=1,nchim)
90	format(i3,14(1x,a4))	
	write(*,*) 'liste des elements traites:'
	write(*,90)nchim,(nom_elem(i),i=1,nchim)

	if(nchim>pnchim) then
	 print*, 'nchim>pnchim =>stop'
	 stop
	endif

	read(1,'(a1)') bidon
	read(1,2000) (nucleo(i),i=1,nchim)
	write(*,*) 'nombre de nucleons:'
	write(*,2001)(nucleo(i),i=1,nchim)
	read(1,'(a1)') bidon
	read(1,2000)(zi(i),i=1,nchim)
	write(*,*) 'charge des noyaux:'
	write(*,2001)(zi(i),i=1,nchim)
	read(1,'(a1)') bidon
	write(*,*) 'mtot,rtot,ltot,teff:'
	read(1,2004)mtot,rtot,ltot,teff
	write(*,2000)mtot,rtot,ltot,teff
	read(1,'(a1)') bidon	
	read(1,'(a1)') bidon
	read(1,'(a1)') bidon
c	write(11,*)'mass	ray	t	ro	kap	grad	lum	nel'

	write(*,*)
	write(*,*) 'Accelerations:'
	write(*,91)nchim,(nom_elem(i),i=1,nchim)
91	format(i3,14('  ',a4))
	
c	boucle sur les couches pour un pas de temps

	ltot=ltot/lsol

	mstar=mtot/msol
	nb=nchim+1
	
	allocate(dg_rad(nb,nb),g_rad(nb),g_rad0(nb),xi(nb),
	1	ioni(0:nint(maxval(zi)),nchim))
	
	f_rad=.true.		

	do i=1,2000,100
	 read(1,2004,end=10)mass,ray,t,ro,kap,grad,lum,nel
	 do j=1,nchim
	  read(1,2001)ychim(j),(ioni(k,j),k=0,nint(zi(j)))
	 enddo
	 grav=-g*mtot*(1.d0-mass)/ray**2
	 r=ray/rsol
c	 print*,'mstar,ltot,ray,mtot,mass,r,grav'
c	 write(*,2000)mstar,ltot,ray,mtot,mass,r,grav
c	 pause

	 xi(1:nchim)=ychim(1:nchim)	
	 call acc_rad(mstar,ltot,r,t,kap,dkapx,nel,xi,ioni,grav,
	1	g_rad0,dg_rad)
	
	 print*,i
	 write(*,2000)mass,ray,t,ro,grav
	 write(*,2000)g_rad0
c	 call pause('apres premier appel')
	 	   
	 if(t > 5.d4)then
c	 if(t > 5.d10)then	 
	  write(*,2000)grav,(g_rad0(j),j=1,nchim)	 
	  do j=1,nchim 	 
	   stor0=xi(j)
	   stor=stor0*unpdx
	   dstor=stor-stor0
	   xi(j)=stor
	   call acc_rad(mstar,ltot,r,t,kap,dkapx,nel,xi,ioni,grav,
	1	g_rad,dg_rad)
	   xi(j)=stor0
	   print*,'derivee / ',nom_elem(j)
	   write(*,2000)((g_rad(k)-g_rad0(k))/dstor,dg_rad(k,j),k=1,nchim)
	   pause
	  enddo
	 endif
	enddo
10	close(unit=1)
	itot=i-1
	print*,'nombre de couches ',itot

	stop

	end program test_acc_rad