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

	subroutine AT3PopIas59tdetau(tau,teff,grav,t,dtsdtau,dtsdteff,dtsdg,
	1	ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau)

c	loi t(tau) des modeles Atlas12 avec CGM, alpha=0.5

c	Auteur: P. Morel, Departement J.D. Cassini, O.C.A., Observatoire de Nice
c       Adaptations: L. Piau, SAp-Saclay
c	CESAM2k

c entree :
c	tau : profondeur optique Rosseland
c	teff : temperature effective
c	grav : gravite

c sortie :
c	t : temperature
c	dtsd* : derivees t/ tau, teff, grav
c	dtsd* : derivees t/ tau, teff, grav
c	dro_** : derivees ro_ext/ teff, grav
c	f_tau, df_tau, df_tau2 : f, d f / d tau, d2 f / d2 tau

c       difference AT1PopIas59tdetau : appel a gr3interpatmas59cgm


 	USE mod_numerique, only :  linf, bsp1ddn
        USE mod_variables, only :  age
	USE mod_kind

	implicit none

	integer, parameter :: itauP=72, itauT=72, pmm=4

	real (kind=dp) :: tau,teff,grav,t,dtsdtau,dtsdteff,dtsdg,teffc,cte1,
	1	taus(itauT),ts(1,itauT),taut(itauT+2*pmm),fx(1,pmm),f_tau,
	2	df_tau,d2f_tau,ro_ext0,ro_ext,dro_grav,dro_teff, lograv,
	2	tauPe(itauP),Ptaue(itauP)

	real (kind=dp), SAVE :: agecT=0.d0

	real (kind=dp),dimension(25) :: yy1=(/
     3 4.0000D+03, 4.1000D+03, 4.2000D+03, 4.3000D+03, 4.4000D+03,
     4 4.5000D+03, 4.6000D+03, 4.7000D+03, 4.8000D+03, 4.9000D+03,
     5 5.0000D+03, 5.1000D+03, 5.2000D+03, 5.3000D+03, 5.4000D+03,
     5 5.5000D+03, 5.6000D+03, 5.7000D+03, 5.8000D+03, 5.9000D+03,
     5 6.0000D+03, 6.1000D+03, 6.2000D+03, 6.3000D+03, 6.4000D+03/)
	real (kind=dp),dimension(11) :: yy2=(/
     1 0.5000D+00, 1.0000D+00, 1.5000D+00, 2.0000D+00, 2.5000D+00,
     2 3.0000D+00, 3.5000D+00, 4.0000D+00, 4.5000D+00, 5.0000D+00,
     2 5.5000D+00/)

	integer, SAVE :: knot,l,nx,mx,i,l1,l1o,l2,l2o

	logical, SAVE :: init=.true.


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

	if(init)then
	 init=.false.
c	 agec=age
c	 write(*,*)'Grav',lograv
c	 l1o=-1
c	 l2o=-1
	 lograv=log10(grav)
	 write(*,*)'Grav',lograv
c	 call gr3interpatmas59cgm(teff,lograv,taus,ts,tauPe,Ptaue,
c	1	tau_min,ro_ext0,teffc)
	 call piau3(teff,lograv,taus,ts,tauPe,Ptaue,tau_min,ro_ext0,teffc)
	 rad=.false.
	 l=1
	 mx=pmm
	 nx=itauT

	 cte1=(3./4.)**.25
	 do i=1,itauT
	  ts(1,i)=(ts(1,i)/teffc)**4*4.d0/3.d0
	 enddo
	 write(2,1)teffc,tau_min,ro_ext0
1	 format(/,1x,'loi t(tau,teff,grav), ATLAS 12 non purement radiative',/,
	1 1x,'teff=',1pd10.3,' tau_min=',1pd10.3,' ro_ext=',1pd10.3,/)
	 write(6,1)teffc,tau_min,ro_ext0
c	 write(6,*)'ts'
c	 write(6,*)ts(1,:)
c	 write(6,*)'taus'
c	 write(6,*)taus
	 call bsp1ddn(1,ts,taus,taut,nx,mx,knot,.false.,taus(1),l,fx)
c	 write(*,*)'ts',teffc*(ts(1,:)*3.d0/4.d0)**.25
c	 write(*,*)'taut',taut
c	 tau_min=1.d-4
c	 ro_ext0=3.55d-9
         call linf(teff,yy1,25,l1)
         call linf(lograv,yy2,11,l2)
         if(abs(teff-yy1(l1)) .gt. abs(teff-yy1(l1+1))) l1=l1+1
         if(abs(lograv-yy2(l2)) .gt. abs(lograv-yy2(l2+1))) l2=l2+1
	 write(6,*)'Chargement de la table atmospherique T(tau)'
	 write(6,2000)l1,l1o,l2,l2o
	 write(6,2001)age,yy1(l1),teff,yy2(l2),lograv
c	 write(6,*)age,yy1(l1),teff,yy2(l2),lograv
	 write(2,*)'Chargement de la table atmospherique T(tau)'
	 write(2,2000)l1,l1o,l2,l2o
	 write(2,2001)age,yy1(l1),teff,yy2(l2),lograv
 2000	 format(4i3)
 2001	 format('Age=',es10.3' Teff tableau=',f6.1,' Teff modele=',
	1	f6.1,' logg tableau=',f5.2,' logg modele=',f5.2)

	 l1o=l1
	 l2o=l2
	endif

c	write(*,*)'ts',teffc*(ts(1,:)*3.d0/4.d0)**.25
c	write(*,*)'teffc',teffc
c	write(*,*)'logg',lograv
c	write(*,*)'taus',taus

c	write(*,*)'tau',tau
c	lograv=log10(grav)
c	call interpatm1(teff,lograv,taus,ts,tauPe,Ptaue,
c	1	tau_min,ro_ext0,teffc)
c	do i=1,itauT
c	  ts(1,i)=(ts(1,i)/teffc)**4*4.d0/3.d0
c	enddo
	lograv=log10(grav)
c	write(*,*)'age,agecT',age,agecT
	if (age .gt. agecT) then
	 agecT=age
         call linf(teff,yy1,25,l1)
         call linf(lograv,yy2,11,l2)
c         if(abs(teff-yy1(l1)) .gt.
c	1	abs(teff-yy1(l1+1))) l1=l1+1
c         if(abs(lograv-yy2(l2)) .gt.
c	1	abs(lograv-yy2(l2+1))) l2=l2+1
c	 write(*,*)'m',teff,lograv
c	 if(l1 .ne. l1o .or. l2 .ne. l2o) then
c	   call gr3interpatmas59cgm(teff,lograv,taus,ts,tauPe,Ptaue,
c	1	tau_min,ro_ext0,teffc)
	   call piau3(teff,lograv,taus,ts,tauPe,Ptaue,tau_min,ro_ext0,teffc)

c	   write(6,*)'ts'
c	   write(6,*)ts(1,:)
c	   write(6,*)'taus'
c	   write(6,*)taus

	   do i=1,itauT
	    ts(1,i)=(ts(1,i)/teffc)**4*4.d0/3.d0
	   enddo

	   call bsp1ddn(1,ts,taus,taut,nx,mx,knot,.false.,taus(1),l,fx)
	   write(2,*)'Changement de table atmospherique T(tau)'
	   write(2,2000)l1,l1o,l2,l2o
	   write(2,2001)age,yy1(l1),teff,yy2(l2),lograv
	   write(6,*)'Changement de table atmospherique T(tau)'
	   write(6,2000)l1,l1o,l2,l2o
	   write(6,2001)age,yy1(l1),teff,yy2(l2),lograv
	   l1o=l1
	   l2o=l2
c	   pause
c	 endif
	endif
c	write(6,*)l1,l1o,l2,l2o,teff,lograv,teffc

	call bsp1ddn(1,ts,taus,taut,nx,mx,knot,.true.,tau,l,fx)
	f_tau=fx(1,1)
	df_tau=fx(1,2)
	d2f_tau=fx(1,3)
	t=teff*cte1*f_tau**.25
	dtsdtau=t/4./f_tau*df_tau
	dtsdteff=t/teff
	dtsdg=0.
	ro_ext=ro_ext0
	dro_grav=0.
	dro_teff=0.

	return

	end subroutine AT3PopIas59tdetau