include 'mod_chim'

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

	PROGRAM test_rk_imps

c	test pour la routines rk_imps d'integration temporelle

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

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

	USE mod_chim, ONLY : iw, nbelem, nchim, nuc, nucleo, ordre, rk_imps
	USE mod_donnees, ONLY : dtmin, lit_nl, ini_ctes_phys, nom_fich2,
	1	print_ctes
			
	USE mod_kind

	IMPLICIT NONE

	INTEGER, PARAMETER :: pn=40
	REAL (kind=dp), ALLOCATABLE, dimension(:,:) :: jac, compx 
	REAL (kind=dp), ALLOCATABLE, dimension(:) :: comp0, dcomp, esti, ex,
	1	compy
	REAL (kind=dp), dimension(pn) :: t, ro, dm, t_t, ro_t
	REAL (kind=dp), dimension(5) :: e	
	REAL (kind=dp) :: som, som0, et, ero, hh, be7, b8, n13, o15, f17, 
	1	dtt, dt_t, dro, dro_t, dt, vt, age, age_max
	
	INTEGER, ALLOCATABLE, dimension(:)  :: anuc
	INTEGER :: i, kk, j

	logical :: ok

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

2000	FORMAT(8es10.3)

c	initialisations

	dt=100.d0

c	lecture de test.don

	nom_fich2='mon_modele' ; CALL lit_nl ; CALL ini_ctes_phys
	CALL print_ctes(6)	

c	initialisations
		
	CALL nuc(vt,vt,comp0,dcomp,jac,.false.,0,
	1	e,et,ero,ex,hh,be7,b8,n13,o15,f17)
				
	kk=1	!nombre de couches dans la ZC
	kk=10

c	ordre=1	!ordre du R_K, Euler, ordre 1
	ordre=2	!ordre du R_K, Lobatto IIIc, ordre 2
c	ordre=6	!ordre du R_K, Lobatto IIIc, ordre 6	
		
	ALLOCATE(esti(nbelem),comp0(nbelem),compy(nbelem),compx(nbelem,kk),
	1	anuc(nchim))

c	initialisation de la composition chimique comp0
	
	CALL nuc(vt,vt,comp0,dcomp,jac,.false.,1,
	1	e,et,ero,ex,hh,be7,b8,n13,o15,f17)
	
	IF(iw > 1)comp0(iw)=5.d-6
	
	anuc=float(nint(nucleo(1:nchim)))	!nombre de nucleons
	
	DO i=1,kk
	 compx(:,i)=comp0
	ENDDO

	som0=SUM(anuc*comp0(1:nchim))	

c	dm, t et ro dans la ZM

	dm(1:kk)=1.d0/kk
	t(kk)=1.42d7 ; t_t(kk)=1.35d7
	ro(kk)=113 ; ro_t(kk)=109

	t(1)=1.561d7 ; t_t(1)=1.547d7
	ro(1)=154.d0 ; ro_t(1)=149.d0

	dtt=(   t(kk)-  t(1))/kk ; dt_t=(t_t(kk)-t_t(1))/kk
	dro=    (ro(kk)-ro(1))  /kk ; dro_t=(ro_t(kk)-ro_t(1))/kk
	DO i=2,kk-1
	 t(i)=    t(i-1)+dtt ; t_t(i)=t_t(i-1)+dt_t
	 ro(i)=    ro(i-1)+dro ; ro_t(i)=ro_t(i-1)+dro_t
	ENDDO

	PRINT*,'t, t_t, ro, ro_t, dm'
	DO i=1,kk
	 WRITE(*,2000)t(i),t_t(i),ro(i),ro_t(i),dm(i)
	ENDDO

c	integration
	
	CALL rk_imps(t_t,ro_t,compx,t,ro,compy,dt,esti,ok,kk,dm)
	PRINT* ; PRINT*,'X(t), X(t+dt)'
	WRITE(*,2000)comp0 ; WRITE(*,2000)compy ; PAUSE
	
c	conservation des baryons

	som=SUM(anuc*compy(1:nchim))
	PRINT* ; PRINT*,'conservation des baryons'	
	WRITE(*,"('som0=',es10.3,', som=',es10.3,', dsom=',es10.3)")som0,som,
	1	(som-som0)/som0
c	PAUSE

	PRINT*
	age=0.d0 ; age_max=4750.d0
	b1: DO WHILE(age < age_max)
	 CALL rk_imps(t_t,ro_t,compx,t,ro,compy,dt,esti,ok,kk,dm)	
	 age=age+dt
	 WRITE(*,"('age=',es10.3,', dt=',es10.3)")age,dt
	 WRITE(*,2000)compy ; PAUSE
	 IF(age >= age_max)EXIT b1
	 DO i=1,kk
	  compx(:,i)=compy
	 ENDDO
	 dt=min(dt*1.2d0,250.d0,age_max-age)
	ENDDO b1

c	conservation des baryons

	som=sum(anuc*compy(1:nchim))
	PRINT* ; PRINT*,'conservation des baryons'	
	WRITE(*,"('som0=',es10.3,', som=',es10.3,', dsom=',es10.3)")som0,som,
	1	(som-som0)/som0	
	STOP
		 	 
	END PROGRAM test_rk_imps