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

	program test_thermo_atm

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

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

	use mod_chim
	use mod_donnees
	use mod_etat
	use mod_kind		
	use mod_opacite		
	use mod_quasi_static

	implicit none

	real (kind=dp), allocatable, dimension(:) :: xchim
		
	real (kind=dp) :: pt, p, t, m, l, r, dlpp, tb,
	1	ro, ro0, drop, drot, teff, grav, dtsdtau, dtsdteff, dtsdg,
	2	ro_ext, dro_grav, dro_teff, f_tau, df_tau, d2f_tau,
	3	grad, grad0, dgradpt, dgradp, dgradt, dgradr, dgradrs,
	4	dgradm, dgradlpp, dgradtau, nh1, nhe1, nhe2, deltap, deltat,
	5	kap, kap0, dkapp, dkapt, tau, rstar, delta, cp,
	6	hp, hp0, dhppt, dhpp, dhpt, dhpr, dhpm, alfa, beta, gamma1,
	7	gradad, gradad0, dgradadp, dgradadt,
	8	gradrad, stor, stor0, dd, unpdd, dstor,
	9	gam, gam0, dgampt, dgamp, dgamt, dgamr, dgamrs, dgamm,
	1	dgamtau, dgamlpp
	
	integer :: ktest, i
		
	logical :: radiatif

c--------------------------------------------------------------------
	
2000	format(1x,8es10.3)
2002	format(1x,9es10.3)
2001	format(1x,3es21.14)

	dd=1.d-6
	unpdd=1.d0+dd
	
	nom_fich2='mon_modele' ; call lit_nl ; call ini_ctes_phys ; call print_ctes(6)

	mstar=mtot

	l=1.d0 ; r=1.d0  ; ihe4=3
	nchim=3 ; precix=1.d-4
	iw=-100 ; n_atm=100
	m_qs=2
	allocate(nom_elem(nchim),xchim(nchim))
	nom_elem='    '

	r_zc=0.d0 ; r_ov=0.d0 ; lim=0 ; wrot=0.d0
	
	xchim(1)=x0 ; xchim(2)=0.7d-5 ; xchim(ihe4)=y0-xchim(2)	
	
	m=1.d0
	l=7.000D-01
	r=8.831D-01
	rstar=8.832D-01
	tau_max=20.d0
	
	pt=1.754D+05
	p=1.5D+05
	dlpp=8.463D-01
	pturb=.true.
	tau=1.736D+00
	teff=5.623D+03
	grav=1.d4
	t=6.554D+03
	
	pt=2.403D+05	
	p=2.403D+05
	dlpp=1.d0
	pturb=.false.
	t=9.339D+03
	tau=9.267D+00
	teff=5.624D+03
	xchim(1)=7.130D-01
	l=7.000D-01
	grav=3.513D+04
	r=8.831D-01
	rstar=8.831D-01	
	ktest=1
	
	pt=2000	
	p=2000
	dlpp=1.d0
	pturb=.false.
	t=4400.d0
	tau=0.004d0
	teff=5.624D+03
	xchim(1)=7.130D-01
	l=1.d0
	grav=3.513D+04
	r=.9999d0
	rstar=1.d0	
	ktest=1
	
10	call tdetau(tau,teff,grav,tb,dtsdtau,dtsdteff,dtsdg,
	1	ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau)
	
	print*,'apres t(tau)'
	
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro0,drop,drot,kap0,dkapp,dkapt,
	3	gradad0,dgradadp,dgradadt,
	4	grad0,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam0,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp0,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.false.)
	

	if(radiatif)then
	 write(*,*)'RADIATIF, pturb=',pturb
	else
	 write(*,*)'CONVECTIF, pturb=',pturb
	endif
	write(*,*)'pt,p,t,xchim(1),m,l,r,dlpp,tau,rstar'	
	write(*,2000)pt,p,t,xchim(1),m,l,r,dlpp,tau,rstar
	write(*,*)'ro0,kap0,gradad0,grad0,alfa,beta,gamma1'	
	write(*,2000)ro0,kap0,gradad0,grad0,alfa,beta,gamma1
	write(*,*)'gam0,hp0'	
	write(*,2000)gam0,hp0
	
	write(*,*)
	write(*,*)'derivees par rapport a pt'
	stor0=pt
	stor=stor0*unpdd
	if(stor == 0.)stor=dd
	dstor=stor-stor0
	pt=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,	
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	pt=stor0
	write(*,*)'dgradpt dgampt dhppt radiatif',radiatif
	write(*,2000)dgradpt,(grad-grad0)/dstor,
	4	dgampt,(gam-gam0)/dstor,
	5	dhppt,(hp-hp0)/dstor
c	pause'derivees par rapport a Pt'

	write(*,*)
	write(*,*)'derivees par rapport a p'
	stor0=p
	stor=stor0*unpdd
	if(stor == 0.d0)stor=dd
	dstor=stor-stor0
	p=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	p=stor0
	write(*,*)'drop dgradp dkapp dgradadp dgamp dhpp radiatif',radiatif	
	write(*,2000)drop,(ro-ro0)/dstor,
	1	dgradp,(grad-grad0)/dstor,
	2	dkapp,(kap-kap0)/dstor,
	4	dgradadp,(gradad-gradad0)/dstor,dgamp,(gam-gam0)/dstor,
	5	dhpp,(hp-hp0)/dstor	
c	pause'derivees par rapport a p'

	write(*,*)
	write(*,*)'derivees par rapport a t'
	stor0=t
	stor=stor0*unpdd
	if(stor == 0.d0)stor=dd
	dstor=stor-stor0
	t=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	t=stor0
	write(*,*)'drot dgradt dkapt dgradadt dgamt dhpt radiatif',radiatif
	write(*,2000)drot,(ro-ro0)/dstor,
	1	dgradt,(grad-grad0)/dstor,
	2	dkapt,(kap-kap0)/dstor,
	4	dgradadt,(gradad-gradad0)/dstor,dgamt,(gam-gam0)/dstor,
	5	dhpt,(hp-hp0)/dstor	
c	pause'derivees par rapport a T'
	
	write(*,*)
	write(*,*)'derivees par rapport a m'
	stor0=m
	stor=stor0*unpdd
	if(stor == 0.)stor=dd
	dstor=stor-stor0
	m=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	m=stor0
	write(*,*)'dgradm dhpm dgamm radiatif',radiatif
	write(*,2000)dgradm,(grad-grad0)/dstor,dgamm,(gam-gam0)/dstor,
	1	dhpm,(hp-hp0)/dstor
c	pause'derivees par rapport a m'
	
	write(*,*)
	write(*,*)'derivees par rapport a r'
	stor0=r
	stor=stor0*unpdd
	if(stor == 0.)stor=dd
	dstor=stor-stor0
	r=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	r=stor0
	write(*,*)'dgradr dgamr dhpr radiatif',radiatif
	write(*,2000)dgradr,(grad-grad0)/dstor,dgamr,(gam-gam0)/dstor,
	1	dhpr,(hp-hp0)/dstor
c	pause'derivees par rapport a r'

	write(*,*)
	write(*,*)'derivees par rapport a rstar'
	stor0=rstar
	stor=stor0*unpdd
	if(stor == 0.)stor=dd
	dstor=stor-stor0
	rstar=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	rstar=stor0
	write(*,*)'dgradrs dgamrs radiatif',radiatif
	write(*,2000)dgradrs,(grad-grad0)/dstor,dgamrs,(gam-gam0)/dstor
c	pause'derivees par rapport a rstar'
	
	write(*,*)	
	write(*,*)'derivees par rapport a dlpp'
	write(*,*)
	stor0=dlpp
	stor=stor0*unpdd
	if(stor == 0.)stor=dd
	dstor=stor-stor0
	dlpp=stor
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	dlpp=stor0
	write(*,*)'dgradlpp, dgamlpp, radiatif',radiatif
	write(*,2000)dgradlpp,(grad-grad0)/dstor,dgamlpp,(gam-gam0)/dstor
c	pause'derivees par rapport a dlpp'

	write(*,*)
	write(*,*)'derivees par rapport a tau'
	stor0=tau
	stor=stor0*unpdd
	if(stor == 0.)stor=dd
	dstor=stor-stor0
	tau=stor	
	call tdetau(tau,teff,grav,tb,dtsdtau,dtsdteff,dtsdg,
	1	ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau)
	call thermo_atm(pt,p,t,xchim,m,l,r,dlpp,
	1	tau,df_tau,d2f_tau,rstar,	
	2	ro,drop,drot,kap,dkapp,dkapt,
	3	gradad,dgradadp,dgradadt,
	4	grad,dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm,dgradtau,dgradlpp,
	5	gam,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm,dgamtau,dgamlpp,	
	6	hp,dhppt,dhpp,dhpt,dhpr,dhpm,delta,deltap,deltat,
	7	cp,gradrad,alfa,beta,gamma1,radiatif,.true.)
	tau=stor0
	write(*,*)'dgradtau dgamtau radiatif',radiatif
	write(*,2000)dgradtau,(grad-grad0)/dstor,dgamtau,(gam-gam0)/dstor
c	pause'derivees par rapport a tau'
	
	print*,'pt, p, t, X, v, dlpp, Pturb, radiatif ',pturb,radiatif
	write(*,2000)pt,p,t,xchim(1),dlpp
	pause'nouveau test'	
	
	goto(11,12,13,14,100),ktest

11	pt=2.D+05
	t=9680.d0
	tau=10.d0
	p=2.D+05
	dlpp=1.d0
	pturb=.false.
	m=1.d0
	l=1.d0
	r=1.d0
	rstar=1.d0
	ktest=ktest+1
	goto10
	
12	pt=1.2529D+05
	p=.6D+05
	t=9680.d0
	tau=1.86d0
	dlpp=4.d0
	pturb=.true.
	m=1.d0
	l=1.d0
	r=1.d0
	rstar=1.d0
	ktest=ktest+1
	goto10
	
13	pt=1.9D+05
	p=1.D+05
	tau=10.d0
	t=8620.d0	
	dlpp=4.d0
	pturb=.true.
	m=1.d0
	l=1.d0
	r=1.d0
	rstar=1.d0
	ktest=ktest+1
	goto10
	
14	pt=2.402D+05
	p=2.D+05
	t=8620.d0
	teff=5596.d0
	tau=9.267d0
	dlpp=.5572d0
	pturb=.true.
	m=1.d0
	l=.638d0
	r=.8514d0
	rstar=.8515d0
	ktest=ktest+1
	goto10

100	stop

	end program test_thermo_atm