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

	SUBROUTINE evol(compt,dt,dts,reprend)

c subroutine public du module mod_evol

c gestion de l'évolution temporelle de la composition chimique
c les points d'intégration en comp. chim. sont les points de raccord
c les m_zc, en m/Mstar sont déterminés dans lim_zc

c sans diffusion, les ZC sont mélangées, dans ZR on conserve W
c avec diffusion, W est diffusé

c on ne tient pas compte de la différence Ptot/Pgaz négligeable
c dans les régions ou les réactions nucléaires sont actives

c on reconstruit le vecteur de composition chimique avec
c discontinuité des xchim aux limites ZR/ZC

c avec pression turbulente 8 inconnues
c sans pression turbulente 7 inconnues, Ptot=Pgaz

c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k

c 23 09 96 : suppression du mouvement des ZR/ZC avec diffusion
c 09 10 96 : modif de la gestion de somme dm/somme r2 dm dans ZC
c 26 06 97 : remplacement du moment angulaire par la v. angulaire
c 25 08 97 : mise en place des variables eulériennes
c 20 10 99 : ajout des variables de structure au temps t dans diffus
c 19 11 99 : suppression de nh1, nhe1, nhe2, lamb
c 18 04 00 ; coeff_diff ---> diffm, difft, age dans diffus
c 30 07 00 : introduction F95

c entrées
c	compt=0: compteur des itérations globales, compt=0 pour
c	la première intération globale de chaque nouveau pas temporel
c	dt: pas temporel (fixé dans UPDATE)

c sorties
c	reprend: il faut réduire dt variation de TdS ou non CV
c	dts : estimation du pas temporel à utiliser pour le pas suivant

c modèle totalement convectif: lim=1,jlim(1)=n,lconv(1)=.FALSE.
c modèle totalement radiatif: lim=0,jlim(i)=-100,lconv(i)=.FALSE.

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

	USE mod_donnees, ONLY : ab_ini, agemax, baratine, cephe,
	1 diffusion, dtmax, dtmin, en_m23, g, grille_fixe,
	2 ihe4, Ipg, Kdes_rot, Krot, langue, lisse,
	3 l_pertm, l_pertmt, mdot, msol, mtot, m_qs, m_rot, m_ch, m_ptm,
	4 nchim, ne, nom_diffm, nom_elem, nom_fich2, no_discon,
	5 npt_lisse, nrot, nucleo, n_min_ZC, pi, pnzc, pturb, ordre,
	6 ord_qs, ord_rot, precit, rsol, r_qs, secon6, thw, typ, nvth
	USE mod_etat, ONLY : dege_elec, etat
	USE mod_kind
	USE mod_nuc, ONLY : mzc_ext, nuzc_ext
	USE mod_numerique, ONLY : bsp1dn, bsp_dis, coll, entre, inside, linf,
	1 lisse_sumd, newspl_gal, noeud, no_croiss, shell
	USE mod_variables, ONLY : age, bp, bp_t, chim, chim_t,
	1 dim_ch, dim_rot, inter, jlim, jlim_t, knot, knotc, knotc_t,
	2 knot_ptm, knot_t, knotr, knotr_t, lconv, lconv_t, lim,
	3 lim_t, log10_teff, log10_teff_t, mc, mct, mct_t,
	4 mc_fixe, mc_t, mrot, mrot_t, mrott, mrott_t, mstar,
	5 m_zc, m_zc_t, m_zc23, m_stat, m_stat_t, nc_fixe, n_ch,
	6 n_ch_t, n_ptm, n_qs, n_qs_t, n_rot, n_rot_t, old_ptm, q, qt,
	7 qt_t, q_t, r_zc, r_zc_t, r_stat, r_stat_t, rota, rota_t,
	8 sortie, tot_conv, tot_rad, xt_ptm, x_ptm, vth

	IMPLICIT NONE

	TYPE d_abon
	 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: masse
	 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: abon
	 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dabon
	 INTEGER, ALLOCATABLE, DIMENSION(:) :: couche
	END TYPE d_abon
	TYPE (d_abon), SAVE :: err_max

	REAL (kind=dp), INTENT(in) :: dt
	INTEGER, INTENT(in) :: compt
	REAL (kind=dp), INTENT(out) :: dts
	LOGICAL, INTENT(out) :: reprend

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: compx,
	1 chim_zc, chim_zr, rota_tmp
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: dm, mc_tmp,
	1 mrott_tmp, ro, ro_t, t, t_t
	REAL (kind=dp), DIMENSION(nchim) :: ab_max, bidd, comp, comp1, comp2,
	1 compy, dfdx, esti, estim, f, mc_maxi
	REAL (kind=dp), DIMENSION(ne) :: dfe, fe
	REAL (kind=dp), DIMENSION(nrot) :: omega_t

	REAL (kind=dp), SAVE :: ln0t8, ln1t8, ln2t8
	REAL (kind=dp) :: a, alfa, b, beta, bid, cp,  dcpp, dcpt, dcpx,
	1 dege, dege_max, delta, deltap, deltat, deltax, dgradadp, dgradadt,
	2 dgradadx, drop, drot, drox, dtn, dtnew, dup, dut, dux,
	3 est, gamma1, gradad, mass_zc, mk, mk32, mk32_t,
	4 mk_t, mmax, norm, p, p_t, r, r_ext, r_t, u, w_max

	INTEGER, DIMENSION(nchim) :: kmaxi
	INTEGER, DIMENSION(1) :: ich_max
	INTEGER, SAVE :: ll, kmax
	INTEGER :: i, ij, ich_abon_max, izc, j, k, knotr_tmp,
	1 nadd, nc_tmp, n_pt_zc

	LOGICAL, SAVE :: init=.TRUE.
	LOGICAL :: loop, l_nzc, ok

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

2000	FORMAT(8es10.3)
2001	FORMAT(es10.3,3es22.15)
2002	FORMAT(13es8.1)

c	PRINT*,'n,nchim,nchim,n_ch,knotc',n_qs,nchim,nchim,n_ch,knotc
c	PRINT*,'dt,t_inf,lim,jlim,l_conv',dt,t_inf,lim,
c	1 (jlim(i),i=1,lim),(lconv(i),i=1,lim)
c	PRINT*,ALLOCATED(x_ptm) ; PAUSE'entrée evol'

c-----------initialisations début------------------------------
	IF(init)THEN
	 init=.FALSE.

	 IF(pturb)THEN
	  WRITE(*,11) ; WRITE(2,11)
11	  FORMAT(/,'Evolution temporelle de la composition chimique',/,
	1 'on ne tient pas compte de la différence Ptot-Pgaz',/,
	2 'qui est négligeable dans les régions où les',/,
	3 'réactions thermonucléaires sont actives',/)
	 ENDIF

	 SELECT CASE(Krot)
	 CASE(3,4)
	  SELECT CASE(langue)
	  CASE('english')
           WRITE(*,1022)TRIM(thw(Krot)) ; WRITE(2,1022)TRIM(thw(Krot))
1022       FORMAT(/,'Diffusion of ang. mom. with collocation',/
	1  'Formalism : ',a)
           WRITE(*,1023)m_rot ; WRITE(2,1023)m_rot
1023       FORMAT('order of B-splines:',i3)
	  CASE DEFAULT
           WRITE(*,22)TRIM(thw(Krot)) ; WRITE(2,22)TRIM(thw(Krot))
22         FORMAT(/,'Diffusion du moment cinétique par collocation.',/
	1  'Formalisme : ',a)
           WRITE(*,23)ord_rot ; WRITE(2,23)ord_rot
23         FORMAT('ordre des B-splines:',i3)
	  END SELECT
	  IF(.NOT.ALLOCATED(xcoll_rot))ALLOCATE(xcoll_rot(0))	!initialisation
	 END SELECT
	 kmax=MAX(nchim,ne)

c baratine=.FALSE. permet de dérouter sur le fichier mon_modele_103 les
c informations concernant le déroulement des calculs de l'évolution de la
c composition chimique et de la vitesse angulaire
	 IF(baratine)THEN
	  usl_evol=6
	 ELSE
	  usl_evol=103
	  OPEN(unit=103,form='formatted',status='unknown',!access='append',
	1 file=TRIM(nom_fich2)//'_evol')
	 ENDIF

c allocations
	 ALLOCATE(err_max%masse(nchim),err_max%abon(nchim),
	1 err_max%dabon(nchim),err_max%couche(nchim))
c err_max%masse: masse au max de variation
c err_max%abon: isotope présentant le max de variation
c err_max%dabon: variation de l'abondance
c err_max%couche: indice de la couche du max

c limites pour ne pas rater les blue loops
	 ln1t8=LOG(0.9d8) ; ln2t8=LOG(2.d8) ; ln0t8=LOG(0.8d8)

c si on tient compte explicitement des discontinuités, no_discon=.FALSE.
	 PRINT*
	 IF(no_discon)THEN
	  PRINT*,'On ne tiendra pas compte des discontinuités de composition chimique'
	  IF(lisse)PRINT*,'on effectuera un léger lissage avec contour'
	 ELSE
	  ALLOCATE(chimd(nchim,2*pnzc))
	  PRINT*,'On tiendra compte des discontinuités de composition chimique'
	  IF(lisse)PRINT*,'qui seront ensuite étalées sur ',2*npt_lisse+1,' couches'
	 ENDIF

	ENDIF	!init
c-----------initialisations fin------------------------------

c	PRINT*,'m_zc23(1:lim),m_zc(1:lim),dt'
c	WRITE(*,2000)m_zc23(1:lim),m_zc(1:lim),dt
c	PAUSE'entrée evol'

c initialisations des tableaux de précision
	est=-100.d0 ; estim=-100.d0 ; ab_max=-100.d0
	mc_maxi=-100.d0 ; kmaxi=-100

	IF(.NOT. ALLOCATED(err_max%masse))ALLOCATE(err_max%masse(nchim),
	1 err_max%abon(nchim),err_max%dabon(nchim),err_max%couche(nchim))
	err_max%masse=0.d0 ; err_max%abon=0.d0 ; err_max%dabon=0.d0
	err_max%couche=0

c sortie si dt < dtmin
	IF(dt < dtmin)THEN
	 SELECT CASE(langue)
	 CASE('english')
	  WRITE(*,1010)dt,dtmin ; WRITE(2,1010)dt,dtmin
1010	  FORMAT('STOP : beginning of evol, dt=',es10.3,' < dtmin=',es10.3)
	 CASE DEFAULT
	  WRITE(*,10)dt,dtmin ; WRITE(2,10)dt,dtmin
10	  FORMAT('ARRET : Entrée de evol, dt=',es10.3,' < dtmin=',es10.3)
	 END SELECT
	 CALL sortie('evol 1')
	ENDIF

c détermination de r_ext
	CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(n_qs),ll,fe,dfe)
	IF(no_croiss)PRINT*,'Pb. at 1 in evol'
	IF(en_m23)THEN
	 r_ext=SQRT(fe(3))
	ELSE
	 r_ext=fe(3)
	ENDIF
c	WRITE(*,2000)f(1:ne) ; WRITE(*,2000)r_ext ; PAUSE'r_ext'

c modèle totalement convectif: lim=1, jlim(1)=n_qs, lconv(1)=.FALSE.
c nzc=1, convd(1)=1, convf(1)=nc_tmp
c modèle totalement radiatif: lim=0, jlim(i)=-100, lconv(i)=.FALSE.
c nzc=0, convd(1)=nc_tmp, convf(0)=1
	tot_conv=lim == 1 .AND. (jlim(1) == n_qs) .AND. .NOT.lconv(1)
	tot_rad=lim == 0
c	PRINT*,lim,lconv(1:lim),jlim(1:lim)
c	PRINT*,m_zc(1:lim)
c	PAUSE'evol 0,lim,lconv(1:lim),jlim(1:lim)'

c adaptation de la limite de la zone externe TOUJOURS convective
	I0: IF(.NOT.tot_conv .AND. .NOT.tot_rad)THEN
	 I2: IF(.NOT.lconv(lim))THEN	!la dernière limite est fin de ZC
	  I1: IF(lim == 1)THEN		!si c'est l'unique limite
	   SELECT CASE(langue)	!que l'on déplace à la fin du modèle
	   CASE('english')	!qui devient totalement convectif
	    WRITE(usl_evol,1008)
1008	    FORMAT(/,'fully mixed model')
	   CASE DEFAULT
	    WRITE(usl_evol,8)
8	    FORMAT(/,'Modèle complètement mélangé')
	   END SELECT
	   tot_conv=.TRUE.	!rappel de la définition de tot_conv
	   lim=1 ; jlim(1)=n_qs ; lconv(1)=.FALSE.

c la dernière limite, fin de ZC est supprimée et lim > 1
	  ELSE I1
	   SELECT CASE(langue)
	   CASE('english')
	    WRITE(usl_evol,1005)lim
1005	    FORMAT(/,'extend to the surface of the outer ZC, limit #',i3)
	   CASE DEFAULT
	    WRITE(usl_evol,5)lim
5	    FORMAT(/,'prolongement à la surface de la ZC externe, limite',i3)
	   END SELECT
	   lim=lim-1
	  ENDIF I1
	 ENDIF I2
	ENDIF I0
c	PRINT*,'lim,lconv(1:lim),jlim(1:lim)',lim,lconv(1:lim),jlim(1:lim)
c	PRINT*,'m_zc(1:lim)',m_zc(1:lim)
c	PAUSE'pause 00'

c limites ZR/ZC effectives
	SELECT CASE(langue)
	CASE('english')
	 WRITE(usl_evol,1006)
1006	 FORMAT(/,'limits RZ/CZ used for the convective mixing :')
	CASE DEFAULT
	 WRITE(usl_evol,6)
6	 FORMAT(/,'limites ZR/ZC utilisées pour le mélange convectif :')
	END SELECT

c cas totalement radiatif, masse fictive de la ZC externe
	IF(tot_rad)THEN
	 SELECT CASE(langue)
	 CASE('english')
	  WRITE(usl_evol,1007)
1007	  FORMAT('Unusual condition : fully radiative model')
	 CASE DEFAULT
	  WRITE(usl_evol,7)
7	  FORMAT(/,'Situation anormale : modèle complètement radiatif')
	 END SELECT
	 lim=0 ; jlim=-100 ; lconv=.FALSE.
	 mzc_ext=1.d-20 ;  nuzc_ext=mstar**(2.d0/3.d0)

c cas totalement convectif, masse de la ZC = mstar
	ELSEIF(tot_conv)THEN
	 SELECT CASE(langue)
	 CASE('english')
	  WRITE(usl_evol,1009)
1009	  FORMAT(/'fully convective model')
	 CASE DEFAULT
	  WRITE(usl_evol,9)
9	  FORMAT(/'modèle complètement convectif')
	  mzc_ext=mstar ; nuzc_ext=0.d0
	 END SELECT

	ELSE
	 DO i=1,lim
	  IF(lconv(i))THEN
	   SELECT CASE(langue)
	   CASE('english')
	    WRITE(usl_evol,1012)jlim(i),m_zc(i),1.d0-m_zc(i)/mstar,r_zc(i),
	1   r_zc(i)/r_ext
1012	    FORMAT('beginning of CZ shell : ',i4,/,'mass=',es10.3,
	1   ', 1-M/Mstar=',es10.3,', radius=',es10.3,
	2   ', radius/Rstar=',es10.3)
	   CASE DEFAULT
	    WRITE(usl_evol,12)jlim(i),m_zc(i),1.d0-m_zc(i)/mstar,r_zc(i),
	1   r_zc(i)/r_ext
12	    FORMAT('début de ZC couche : ',i4,/,'masse=',es10.3,
	1  ', 1-m/Mstar=',es10.3,', rayon=',es10.3,', rayon/Rstar=',es10.3)
	   END SELECT
	  ELSE
	   SELECT CASE(langue)
	   CASE('english')
	    WRITE(usl_evol,1017)jlim(i),m_zc(i),1.d0-m_zc(i)/mstar,r_zc(i),
	1   r_zc(i)/r_ext
1017	    FORMAT('end of CZ shell : ',i4,/,'mass=',es10.3,
	1   ', 1-m/Mstar=',es10.3,', radius=',es10.3,
	2   ', radius/Rstar=',es10.3)
	   CASE DEFAULT
	    WRITE(usl_evol,17)jlim(i),m_zc(i),1.d0-m_zc(i)/mstar,r_zc(i),
	1   r_zc(i)/r_ext
17	    FORMAT('fin de ZC couche : ',i4,/,'masse=',es10.3,
	1   ', 1-m/Mstar=',es10.3,', rayon=',es10.3,', rayon/Rstar=',es10.3)
	   END SELECT
	  ENDIF
	 ENDDO
	ENDIF

c mc_tmp : masses (m^2/3) temporaires pour tabul. de la comp.chim.
c nc_temp : nombre temporaire de masses
	IF(ALLOCATED(mc_tmp))DEALLOCATE(mc_tmp)
	ALLOCATE(mc_tmp(n_qs+n_qs_t+2*pnzc))	!allocation généreuse

c on prend la répartition déduite de bp ou la grille fixe
	IF(grille_fixe)THEN
	 mc_tmp(1:nc_fixe)=mc_fixe ; nc_tmp=nc_fixe
	ELSE
	 DO i=1,n_qs
	  CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),ll,fe,dfe)
	  IF(no_croiss)PRINT*,'Pb. at 2 in evol'
	  IF(en_m23)THEN
	   mc_tmp(i)=fe(5)	!m^2/3
	  ELSE
	   mc_tmp(i)=fe(5)**2	!(m^1/3)**2
	  ENDIF
	 ENDDO
	 nc_tmp=n_qs
	ENDIF
c	PRINT*,nc_tmp ; WRITE(*,2000)mc_tmp(1:nc_tmp) ; PAUSE'1'

c couches par ordre croissant
	CALL shell(nc_tmp,mc_tmp)
c	PRINT*,mc_tmp(nc_tmp-4:nc_tmp) ; PAUSE'21'

c suppression des multiplicités, actualisation de nc_tmp
	k=1
	DO i=2,nc_tmp
	 IF(mc_tmp(k) /= mc_tmp(i))THEN
	  k=k+1 ; mc_tmp(k)=mc_tmp(i)
	 ENDIF
	ENDDO
	nc_tmp=k

c écritures pour tests
c	PRINT*,'lim,lconv(1:lim),jlim(1:lim),nc_tmp',
c	1 lim,lconv(1:lim),jlim(1:lim),nc_tmp
c	PRINT*,'m_zc(1:lim),SQRT(mc_tmp(nc_tmp))**3',
c	1 m_zc(1:lim),SQRT(mc_tmp(nc_tmp))**3
c	DO i=1,lim
c	 PRINT*,m_zc23(i),m_zc(i)
c	ENDDO
c	PAUSE'm_zc23(i),m_zc(i)'

c test de dépassement
	IF(MAX(MAXVAL(m_zc23(1:lim)),mc_tmp(nc_tmp)) /= mc_tmp(nc_tmp))THEN
	 PRINT*,'EVOL',nc_tmp,lim,mc_tmp(nc_tmp)
	 PRINT*,m_zc23(1:lim)
	 PRINT*,'dépassement peut-être dû à une erreur de troncature?'
c	 m_zc23(lim)=MIN(m_zc23(lim),mc_tmp(nc_tmp))
	ENDIF

c on introduit les limites ZR/ZC en remplacement du point le plus proche
c délimitation des zones convectives
c il y a nzc zones convectives chacune entre convd(.) et convf(.)
	convd=-100 ; convf=-100
	IF(tot_conv)THEN
	 nzc=1 ; convd(1)=1 ; convf(1)=nc_tmp
	ELSEIF(tot_rad)THEN
	 nzc=0 ; convd(1)=nc_tmp ; convf(0)=1
	ELSE
	 nzc=0		!nombre de ZC
	 j=1		!indice des limites
	 B4: DO i=2,nc_tmp
	  b=mc_tmp(i)-m_zc23(j) ; a=m_zc23(j)-mc_tmp(i-1)
c	  IF(j >= 4)THEN
c	   PRINT*,'i,j',i,j,m_zc23(j),mc_tmp(i)
c	   WRITE(*,2000)a,b,a*b,mc_tmp(i-1),m_zc23(j),mc_tmp(i)
c	  ENDIF
	  IF(a*b < 0.d0)CYCLE B4	!a*b >0 m_zc23 dans [i-1, i+1]
	  IF(b < a)THEN			!point le plus proche
	   mc_tmp(i)=m_zc23(j) ; ij=i
	  ELSE
	   mc_tmp(i-1)=m_zc23(j) ; ij=i-1	!si égalité coté gauche
	  ENDIF
	  IF(lconv(j))THEN		!début de ZC
	   nzc=nzc+1			!une ZC en plus
	   convd(nzc)=ij
c	   PRINT*,'ZCd+,j,i,nzc,convd(nzc)',j,i,nzc,convd(nzc)
	  ELSE				!fin de ZC
	   IF(nzc == 0)THEN	!début de la ZC en 1 car pas encore identifié
	    nzc=1 ; convd(nzc)=1	!la zone externe est ZC
c	    PRINT*,'ZCd1,j,i,nzc,convd(nzc)',j,i,nzc,convd(nzc)
	   ENDIF
	   convf(nzc)=ij
c	   PRINT*,'ZCf,j,i,nzc,convf(nzc)',j,i,nzc,convf(nzc)
	  ENDIF
	  j=j+1
c	  PRINT*,'4,j,lim',j,lim
	  IF(j > lim)EXIT B4
	 ENDDO B4

c si l'atmosphère descend bas, mc_tmp(nc_tmp) < m_zc23(j), on sort avec j=lim et
c non pas j=lim+1; la ZC externe est oubliée; on impose nzc=nzc+1 et
c convd(nzc)=nc_tmp-n_min_ZC et on ramème m_zc23(lim)=mc_tmp(convd(nzc))
	 IF(j == lim)THEN
	  nzc=nzc+1
	  convd(nzc)=nc_tmp-n_min_ZC
	  m_zc23(lim)=mc_tmp(convd(nzc))
	 ENDIF
	 IF(nzc >= 1 .AND. convf(nzc) < 0)convf(nzc)=nc_tmp
	ENDIF

c la ZC externe doit avoir au moins n_min_ZC couches
c on recule la limite ou, si besoin, on étend la ZC précédente à l'extérieur
	nadd=n_min_ZC+1+convd(nzc)-convf(nzc)
c	PRINT*,'nzc,nadd,convd(nzc),convf(nzc),nc_tmp',nzc,nadd,convd(nzc),
c	1 convf(nzc),nc_tmp
	IF(nadd > 0)THEN
	 convd(nzc)=nc_tmp-n_min_ZC
	 IF(nzc > 1 .AND. convf(nzc-1) >= convd(nzc))THEN
	  convf(nzc-1)=convf(nzc) ; nzc=nzc-1
	 ENDIF
c	 PRINT*,'nzc,nadd,convd(nzc),convf(nzc)',nzc,nadd,convd(nzc),convf(nzc)
c	 PAUSE'nadd'
	ENDIF

c masse de la ZC externe en Msol
	IF(.NOT.(tot_conv .OR. tot_rad))THEN
	 nuzc_ext=mc_tmp(convd(nzc)) ; mzc_ext=mstar-m_zc(lim)
c	 WRITE(*,2000)mzc_ext,nuzc_ext,m_zc(lim)**(2.d0/3.d0) ; PAUSE'mzc_ext'
	ENDIF
c	PRINT*,'lim= ',lim,', nzc=',nzc
c	DO i=1,nzc
c	 PRINT*,convd(i),convf(i)
c	 WRITE(*,2000)SQRT(mc_tmp(convd(i)))**3,SQRT(mc_tmp(convf(i)))**3
c	ENDDO
c	PAUSE'verif'

c------------------------Pour vérifications-------------------------------

	IF(.FALSE.)THEN
c	IF(.TRUE.)THEN
	 PRINT*,'vérif,lim,nzc,nc_tmp,mc_tmp(nc_tmp)',lim,nzc,nc_tmp,mc_tmp(nc_tmp)
	 PRINT*,'convd(1:nzc),convf(1:nzc),mc_tmp(convd(1:nzc)),mc_tmp(convf(1:nzc))'
	 PRINT*,convd(1:nzc),convf(1:nzc),mc_tmp(convd(1:nzc)),mc_tmp(convf(1:nzc))
	 PAUSE'convd,convf'
	ELSEIF(.FALSE.)THEN
c	ELSEIF(.TRUE.)THEN
	 j=1
	 DO i=1,nc_tmp
c	  ok=lmix(mc_tmp(i))
	  IF(i == convd(j))THEN
	   IF(i /= 1)PRINT*,mc_tmp(i),i!,ok
	   PRINT*,'début de ZC',j,i ; PRINT*,mc_tmp(i),i!,ok
	  ELSEIF(i == convf(j))THEN
	   PRINT*,mc_tmp(i),i!,ok
	   PRINT*,'fin de ZC',j,i
	   IF(i /= nc_tmp)PRINT*,mc_tmp(i),i!,ok
	   j=MIN(j+1,nzc)
	  ELSE
	   PRINT*,mc_tmp(i),i!,ok
	  ENDIF
	 ENDDO
	 PAUSE'avant ajustement des ZR et ZC'
	ENDIF

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

c limites fictives pour faciliter les algorithmes
	IF(convd(1) /= 1)THEN
	 convf(0)=1
	ELSE
	 convf(0)=10000
	ENDIF
	IF(convf(nzc) /= nc_tmp)THEN
	 convd(nzc+1)=nc_tmp
	ELSE
	 convd(nzc+1)=-10000
	ENDIF
	SELECT CASE(langue)
	CASE('english')
	 WRITE(usl_evol,1002)nc_tmp
1002	 FORMAT(/,'number of shells for the interpolation of chemicals:',i5)
	CASE DEFAULT
	 WRITE(usl_evol,2)nc_tmp
2	 FORMAT('nombre de couches utilisées pour la comp. chim.:',i5)
	END SELECT

c on fixe le centre
	mc_tmp(1)=0.d0

c le vecteur de mélange pour la composition chimique: xmix, mix
c le mélange convectif doit affecter la couche contenant une limite ZR/ZC
c avec rotation on fixe la grille à celle de la 1-ère itération
c sans rotation on adapte la grille
	IF(convd(1) == 1)THEN
	 n_mix=2*nzc ; j=0
	ELSE
	 n_mix=2*nzc+1 ; j=1
	 x_mix(1)=mc_tmp(1) ; mix(1)=.FALSE.
	ENDIF
	DO i=1,nzc
	 j=j+1
	 x_mix(j)=mc_tmp(convd(i)) ; mix(j)=.TRUE.
	 j=j+1
	 x_mix(j)=mc_tmp(convf(i)) ; mix(j)=.FALSE.
	ENDDO

c initialisation éventuelle de 	x_mix_t
	IF(n_mix_t <= 0)THEN
	 n_mix_t=n_mix ; x_mix_t=x_mix
	ENDIF

c estimation de ndis et idis : nombre et tableau des discontinuités
c	PRINT*,'nzc,n_mix',nzc,n_mix
	idis=-HUGE(1)
	IF(tot_conv)THEN
	 ndis=0
	ELSE
	 ndis=2*nzc-1		!ZC externe pas de limite discon externe
	 IF(convd(1) == 1)THEN
	  ndis=ndis-1 ; idis(0)=1 !coeur CV pas de limite discon interne
	 ENDIF
	 idis(ndis+1)=nc_tmp		!à l'extérieur
	 j=0
	 DO i=1,nzc
	  IF(convd(i) /= 1)THEN
	   j=j+1 ; idis(j)=convd(i)
c	   PRINT*,'1,i,j,idis(j)',i,j,idis(j)
	  ENDIF
	  j=j+1 ; idis(j)=convf(i)
c	  PRINT*,'2,i,j,idis(j)',i,j,idis(j)
	 ENDDO
	ENDIF

c on recale sur les discontinuités idis, m_zc et r_zc qui peuvent avoir
c été légèrement décalés
c	PRINT*,'nzc,ndis,idis(0:ndis+1)',nzc,lim,ndis,idis(0:ndis+1)
c	PRINT*,'mc_tmp(idis(1:ndis)),m_zc23(1:lim)'
c	WRITE(*,2000)mc_tmp(idis(1:ndis)),m_zc23(1:lim)
	DO i=1,ndis
	 bid=mc_tmp(idis(i))	!m^2/3
	 IF(en_m23)THEN
	  CALL inter('m',bp,q,qt,n_qs,knot,bid,fe,dfe,r_stat,m_stat)
	  fe(3)=SQRT(fe(3))	!r
	 ELSE
	  CALL inter('m',bp,q,qt,n_qs,knot,SQRT(bid),fe,dfe,	!m^1/3
	1 r_stat,m_stat)
	 ENDIF

c recalage de m_zc et r_zc
	 m_zc(i)=SQRT(bid)**3 ; r_zc(i)=fe(3)
c	 PRINT*,'idis=',i
c         WRITE(*,2000)SQRT(bid)**3,m_zc(i),fe(3),r_zc(i)
c	 WRITE(*,2000)m_zc(i),r_zc(i)
c	 PAUSE'ndis'
	ENDDO

c---------------------Pour vérifications-----------------------------

c	IF(.TRUE.)THEN
	IF(.FALSE.)THEN
	 PRINT*,'nzc,nc_tmp',nzc,nc_tmp
	 PRINT*,'convd',convd(1:nzc+1)
	 PRINT*,'convf',convf(0:nzc)
	 DO i=1,nzc
	  PRINT*,'i,convd(i),convf(i),mc_tmp(convd(i)),mc_tmp(convf(i))'
	  PRINT*,i,convd(i),convf(i),mc_tmp(convd(i)),mc_tmp(convf(i))
	 ENDDO
	 PRINT* ; PAUSE'en tête'
	 j=1
	 DO i=1,nc_tmp
	  ok=lmix(mc_tmp(i))
	  IF(i == convd(j))THEN
	   IF(i /= 1)PRINT*,mc_tmp(i),i,mc_tmp(i)**(3.d0/2.d0),ok
	   PRINT*,'début de ZC',j,i
	   PRINT*,mc_tmp(i),i,mc_tmp(i)**(3.d0/2.d0),ok ; PAUSE'début de ZC'
	  ELSEIF(i == convf(j))THEN
	   PRINT*,mc_tmp(i),i,mc_tmp(i)**(3.d0/2.d0),ok
	   PRINT*,'fin de ZC',j,i ; PAUSE'fin de ZC'
	   IF(i /= nc_tmp)PRINT*,mc_tmp(i),i,mc_tmp(i)**(3.d0/2.d0),ok
	   j=MIN(j+1,nzc)
	  ELSE
	   PRINT*,mc_tmp(i),i,mc_tmp(i)**(3.d0/2.d0),ok
	  ENDIF
	 ENDDO
	 PAUSE'après ajustement des ZR et des ZC'
	ENDIF

c mises à 0 pour la précision sur l'intégration de la comp. chim.
	 err_max%masse=0.d0 ; err_max%abon=0.d0 ; err_max%dabon=0.d0
	 err_max%couche=0

c-----------------------  diffusion des éléments chimiques-----------

	ID: IF(diffusion)THEN

c intégration
	 CALL diffus(ok,dt,mc_tmp,nc_tmp)
	 reprend=.NOT.ok
	 IF(reprend)THEN
	  SELECT CASE(langue)
	  CASE('english')
	   WRITE(*,1018) ; WRITE(2,1018)
1018	   FORMAT('Reinitialize and dt --> dt/2')
	  CASE DEFAULT
	   WRITE(*,18) ; WRITE(2,18)
18	   FORMAT('Réinitialisation et dt --> dt/2')
	  END SELECT
	  IF(dt > dtmin)THEN
	   RETURN
	  ELSE
	   SELECT CASE(langue)
	   CASE('english')
	    WRITE(*,1019)dt,dtmin ; WRITE(2,1019)dt,dtmin
1019	    FORMAT('STOP, dt =,'es10.3,' <',es10.3,' = dtmin')
	   CASE DEFAULT
	    WRITE(*,19)dt,dtmin ; WRITE(2,19)dt,dtmin
19	    FORMAT('ARRET, dt =,'es10.3,' <',es10.3,' = dtmin')
	   END SELECT
	   CALL sortie('evol 2')
	  ENDIF
	 ENDIF

c recherche des variations maximales des isotopes
	 i=1
	 B1: DO
	  IF(i >= convf(nzc))EXIT B1
	  CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(i),
	1 ll,comp1,bidd)
	  IF(no_croiss)PRINT*,'Pb. at 2 in evol'

c s'il y a perte de masse x_ptm en m^2/3 ou m^1/3
	  IF(l_pertm .OR. l_pertmt)THEN
	   IF(en_m23)THEN
	    bid=MIN(mc(i),x_ptm(n_ptm))		!m^2/3
	   ELSE
	    bid=MIN(SQRT(ABS(mc(i))),x_ptm(n_ptm))	!car mc(i)en m^2/3
	   ENDIF
	   CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,
	1  .TRUE.,bid,ll,f,dfdx)
	   IF(no_croiss)PRINT*,'Pb. at 3 in evol'

c mk_t: masse en m^1/3 ou m^2/3 au temps t mk_t
	   IF(.NOT.en_m23)f(1)=ABS(f(1))**(2.d0/3.d0)
	   mk_t=MIN(MAX(mc_t(1),f(1)),mc_t(n_ch_t))
	  ELSE
	   mk_t=MIN(MAX(mc_t(1),mc(i)),mc_t(n_ch_t))
	  ENDIF

c abondances au temps t
	  CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch,
	1 knotc_t,.TRUE.,mk_t,ll,comp2,bidd)
	  IF(no_croiss)PRINT*,'Pb. at 4 in evol'

c variation des abondances par gramme
	  esti=ABS(comp1-comp2)*nucleo

c si on est dans une ZC
	  IF(lmix(mc(i)))THEN
	   B2: DO izc=1,nzc
	    IF(entre(convd(izc),convf(izc),i))THEN
	     k=(convd(izc)+convf(izc))/2
	     DO j=1,nchim
	      IF(esti(j) > err_max%dabon(j))THEN
	       err_max%masse(j)=mc(k)	!mc(i)
	       err_max%abon(j)=comp2(j)*nucleo(j)
	       err_max%dabon(j)=esti(j)
	       err_max%couche(j)=1111*izc
	      ENDIF
	     ENDDO
	     i=convf(izc)+1
	     CYCLE B1 	!saut par dessus la ZC
	    ENDIF
	   ENDDO B2

c si on est dans une ZR
	  ELSE
c saut ZR <--> ZC entre t et t+dt, alors esti est ignoré
	   DO j=1,nchim
	    IF(esti(j) > err_max%dabon(j))THEN
	     err_max%masse(j)=mc(i)
	     err_max%abon(j)=comp2(j)*nucleo(j)
	     err_max%dabon(j)=esti(j)
	     err_max%couche(j)=i
	    ENDIF
	   ENDDO
	   i=i+1
	  ENDIF
	 ENDDO B1
c	 PAUSE'après diffusion'

c----------------sans diffusion mélange classique, début-------------

	ELSE ID
c	 PRINT*,'knotc',knotc,nchim ; WRITE(*,2000)dt
c	 PRINT*,'avant intégr.: n_ch_t,knotc_t,compt',n_ch_t,knotc_t,compt
c	 PRINT*,'mc_t',mc_t(n_ch_t) ; WRITE(*,2000)(mc_t(i),i=1,nc_tmp)
c	 PRINT*,'mct_t',mct_t(knotc_t)
c	 WRITE(*,2000)(mct_t(i),i=1,knotc_t) ; PRINT*,'xchim'
c	 DO i=1,nc_t
c	  CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch,
c	1 knotc_t,.TRUE.,mc_t(i),ll,compx,compy)
c	  WRITE(*,2000)mc_t(i),(compx(j)*nucleo(j),j=1,MIN(7,nchim))
c	  WRITE(*,2002)mc_t(i),(compx(j)*nucleo(j),j=1,nchim)
c	 ENDDO
c	 PAUSE'avant intégration'

c on intègre par rapport à t en mélangeant les ZM
c	 d'abord les ZC puis les ZR
c	 on garde dans chim_zc(nchim,nzc) les valeurs des ZM
c	 la variable en masse est mc=m**2/3
c	 convf(0)=1 si convd(1) /= 1, sinon convf(0)=10000
c	 convd(nzc+1)=nc_tmp si convf(nzc) /= nc_tmp,
c	 sinon convd(nzc+1)=-10000
c	 dans ZM on mélange MW, dans ZR on conserve MW
c	 on reconstruit le vecteur de composition chimique avec
c	 discontinuité des xchim aux LMR.

c les ZC d'abord
	 ALLOCATE(chim_zc(nchim,nzc))
	 Bzc: DO izc=1,nzc
	  mass_zc=0.d0

c n_pt_zc : nb. de points dans la ZC, mk=m_milieu^2/3, mk32=m_milieu^1/3
	  n_pt_zc=convf(izc)-convd(izc)
	  ALLOCATE(compx(nchim,n_pt_zc),ro(n_pt_zc),ro_t(n_pt_zc),
	1 t(n_pt_zc),t_t(n_pt_zc),dm(n_pt_zc))
	  DO i=convd(izc),convf(izc)-1
	   k=i-convd(izc)+1
	   mk=(mc_tmp(i)+mc_tmp(i+1))/2.d0	!point milieu m^2/3
	   mk32=SQRT(ABS(mk))			!en m^1/3=mk^3/2
	   dm(k)=(mc_tmp(i+1)-mc_tmp(i))*mk32
	   mass_zc=mass_zc+dm(k)		!masse de la ZC
c	   WRITE(*,2000)mass_zc,dm(k),mc_tmp(i),mc_tmp(i+1),SQRT(mk),mk,mk32

c composition chimique en mk_t au temps t, et en mk au temps t+dt
c on tient compte de la perte/gain de masse
	   IF(l_pertm .OR. l_pertmt)THEN
	    IF(en_m23)THEN
	     bid=MAX(x_ptm(1),MIN(mk,x_ptm(n_ptm)))
	    ELSE
	     bid=MAX(x_ptm(1),MIN(mk32,x_ptm(n_ptm)))
	    ENDIF
	    CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,
	1  .TRUE.,bid,ll,f,dfdx)	!masse au temps t mk_t
	    IF(no_croiss)PRINT*,'Pb. at 5 in evol'
	    IF(en_m23)THEN
	     bid=f(1)		!m^2/3
	    ELSE
	     bid=f(1)**2	!(m^1/3)**2
	    ENDIF
	    mk_t=MIN(MAX(mc_t(1),bid),mc_t(n_ch_t))	!m^2/3
	   ELSE
	    mk_t=mk			!m^2/3
	   ENDIF
	   mk32_t=SQRT(ABS(mk_t))	!m^1/3

c abondances au temps t
	   CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch,
	1  knotc_t,.TRUE.,MAX(mct_t(1),MIN(mk_t,mct_t(n_ch_t))),ll,
	2  comp1,bidd)
	   IF(no_croiss)PRINT*,'Pb. at 6 in evol'
c	   PRINT*,'après bsp1dn 1, izc,i,k',izc,i,k
c	   WRITE(*,2000)mk,(comp1(j),j=1,MIN(7,nchim))
c	   PAUSE'après bsp1dn 1'

c comp1, compx : abondances au temps t
	   WHERE(comp1 < 1.d-40)comp1=1.d-40 ; compx(:,k)=comp1(:)

c comp2, compy : abondances au temps t+dt
c compx, compy : abondances par mole, comp1, comp2 par gramme

	   IF(compt == 0)THEN		!comp. chim. en t+dt
	    comp2=comp1
	   ELSE
	    CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
	1   knotc,.TRUE.,MIN(mk,mct(knotc)),ll,comp2,bidd)
	    IF(no_croiss)PRINT*,'Pb. at 7 in evol'
c	    WRITE(*,2000)mk,(compy(j),j=1,MIN(7,nchim))
c	    PAUSE'après bsp1dn 2'
	    WHERE(comp2<1.d-40)comp2=1.d-40
	   ENDIF		!new

c comp1 et comp2 par gramme pour appel à l'EOS
	   comp1=comp1*nucleo ; comp2=comp2*nucleo

c P et T au temps t+dt par interpolation inverse en m_stat par inter
c	   PRINT*,'avant, izc,i,k',izc,i,k
	   IF(en_m23)THEN
	    CALL inter('m',bp,q,qt,n_qs,knot,MIN(mk,m_stat(n_qs)),fe,dfe,
	1   r_stat,m_stat)
	   ELSE
	    CALL inter('m',bp,q,qt,n_qs,knot,MIN(mk32,m_stat(n_qs)),
	1   fe,dfe,r_stat,m_stat)
	   ENDIF
	   IF(pturb)THEN	!avec pression turbulente
	    p=EXP(fe(Ipg))
	   ELSE			!sans pression turbulente
	    p=EXP(fe(1))
	   ENDIF
	   t(k)=EXP(fe(2))

c les ro au temps t+dt
	   CALL etat(p,t(k),comp2,.FALSE.,			!les ro
	1  ro(k),drop,drot,drox,u,dup,dut,dux,
	2  delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3  gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
c	   WRITE(*,2000)p,t(k),ro(k),mk

c P et T au temps t par interpolation inverse en m^2/3 ou m^1/3
	   IF(en_m23)THEN
	    CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,
	1   MIN(mk_t,m_stat_t(n_qs_t)),fe,dfe,r_stat_t,m_stat_t)
	   ELSE
	    CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,
	1   MIN(mk32_t,m_stat_t(n_qs_t)),fe,dfe,r_stat_t,m_stat_t)
	   ENDIF
	   IF(pturb)THEN	!avec pression turbulente
	    p_t=EXP(fe(Ipg))
	   ELSE			!sans pression turbulente
	    p_t=EXP(fe(1))
	   ENDIF
	   t_t(k)=EXP(fe(2))

c les ro au temps t
	   CALL etat(p_t,t_t(k),comp1,.FALSE.,
	1  ro_t(k),drop,drot,drox,u,dup,dut,dux,
	2  delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3  gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
c	   WRITE(*,2000)p_t,t_t(k),ro_t(k),mk

c l_nzc : on est dans la ZC externe
	   l_nzc= izc == nzc
	  ENDDO		!i dans chaque ZC
c	  PRINT*,k,n_pt_zc ; PAUSE'evol: n_pt_zc'
c	  WRITE(*,2000)compx(1,1:k)
c	  WRITE(*,2000)dm ; WRITE(*,2000)t
c	  WRITE(*,2000)ro ; WRITE(*,2000)mass_zc ; PAUSE'après les dm(k)'

c normalisation des dm
	  dm=dm/mass_zc

c intégration temporelle
	  CALL rk_imps(t_t,ro_t,compx,t,ro,compy,dt,esti,ok,n_pt_zc,
	1 dm,l_nzc)

c  on garde, provisoirement, les abondances dans chim_zc
	  IF(ok)THEN		!il y a convergence
	   chim_zc(:,izc)=compy(:) !chim_zc: comp. chim. de la ZC

c recherche des variations maximales des isotopes
	   DO j=1,nchim
	    IF(esti(j) > err_max%dabon(j))THEN
	     err_max%masse(j)=mc_tmp(i)
	     err_max%abon(j)=compy(j)*nucleo(j)
	     err_max%dabon(j)=esti(j)
	     err_max%couche(j)=1111*izc
	    ENDIF
	   ENDDO

c difficulté il n'y a pas eu de convergence NR dans rk_imps
	  ELSE
c	   PRINT* ; PRINT*,'evol sortie 1 ok:',ok
	   reprend=.TRUE. ; RETURN
	  ENDIF
c	  PRINT*,'zone mélangee',izc
c	  WRITE(*,2000)mc_tmp(convd(izc)),mc_tmp(convf(izc))
c	  WRITE(*,2002)(chim_zc(1:nchim,izc) ; PAUSE'ZM'
	  DEALLOCATE(compx,ro,ro_t,t,t_t,dm)
	 ENDDO	Bzc!izc

c ensuite les ZR
	 ALLOCATE(compx(nchim,1),ro(1),ro_t(1),t(1),t_t(1),dm(1),
	1 chim_zr(nchim,nc_tmp))
	 dm=1.d0 ; l_nzc=.FALSE.
	 Bzr: DO izc=1,nzc+1
c	  PRINT*,'ZR',izc,convf(izc-1),convd(izc)
	  DO i=convf(izc-1),convd(izc)

c mk_t = m^2/3 au temps t correspondant à mc_tmp(i) au temps t+dt
c mk32_t = m^1/3 au temps t (correspondant à mc_tmp(i) au temps t+dt)
	   IF(l_pertm .OR. l_pertmt)THEN
	    IF(en_m23)THEN
	     bid=mc_tmp(i)		!m^2/3
	    ELSE
	     bid=SQRT(ABS(mc_tmp(i)))	!m^1/3
	    ENDIF
	    bid=MIN(bid,x_ptm(n_ptm))	!limite en cas de croissance de M
	    CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,
	1   .TRUE.,bid,ll,f,dfdx)!masse au temps t = mc_tmp(i)
	    IF(no_croiss)PRINT*,'Pb. at 8 in evol'
	    IF(.NOT.en_m23)f(1)=f(1)**2		!m^2/3
	    mk_t=inside(mc_t(1),mc_t(n_ch_t),f(1))	!m^2/3
	   ELSE
	    mk_t=inside(mc_t(1),mc_t(n_ch_t),mc_tmp(i)) !m^2/3
	   ENDIF
	   mk32_t=SQRT(mk_t)		!m^1/3

c comp1: comp.chim. au temps t, comp2: comp.chim. au temps t+dt
c pour l'initialisation, comp2=comp1
c comp1 et comp2 serviront au calcul de ro et ro_t aux temps t+dt et t
	   CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch,
	1  knotc_t,.TRUE.,mk_t,ll,comp1,bidd)
	   IF(no_croiss)PRINT*,'Pb. at 9 in evol'
c	   PRINT*,'après bsp1dn 1, izc,i',izc,i
c	   WRITE(*,2000)mc_tmp(i),(comp1(j),j=1,MIN(7,nchim))
c	   PAUSE'après bsp1dn 1 radiatif'

	   comp1=MAX(comp1,1.d-100) ; compx(:,1)=comp1
	   IF(compt == 0)THEN		!comp. chim. en t+dt
	    comp2=comp1
	   ELSE
	    CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
	1   knotc,.TRUE.,MIN(mc_tmp(i),mct(knotc)),ll,comp2,bidd)
	    IF(no_croiss)PRINT*,'Pb. at 10 in evol'
c	    WRITE(*,2000)compy ; PAUSE'après bsp1dn 2'
	    comp2=MAX(comp2,1.d-100)
	   ENDIF	!compt=0
	   comp1=comp1*nucleo ; comp2=comp2*nucleo

c P et T au temps t+dt, par interpolation inverse en m^2/3 ou m^1/3
c	   PRINT*,'avant, izc,i',izc,i
	   IF(en_m23)THEN
	    bid=mc_tmp(i)		!m^2/3
	   ELSE
	    bid=SQRT(ABS(mc_tmp(i)))	!m^1/3
	   ENDIF
	   CALL inter('m',bp,q,qt,n_qs,knot,MIN(bid,m_stat(n_qs)),fe,dfe,
	1  r_stat,m_stat)
	   IF(pturb)THEN	!avec pression turbulente
	    p=EXP(fe(Ipg))
	   ELSE			!sans pression turbulente
	    p=EXP(fe(1))
	   ENDIF
	   t(1)=EXP(fe(2))

c ro au temps t+dt
	   CALL etat(p,t(1),comp2,.FALSE.,			!les ro
	1  ro(1),drop,drot,drox,u,dup,dut,dux,
	2  delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3  gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
c	   WRITE(*,2000)p,t(1),ro(1),mc_tmp(i)

c s'il y a perte de masse
	   IF(l_pertm .OR. l_pertmt)THEN
	    CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,
	1   .TRUE.,MIN(bid,x_ptm(n_ptm)),ll,f,dfdx)
	    IF(no_croiss)PRINT*,'Pb. at 11 in evol'
	    mk_t=f(1)		!m^2/3 ou m^1/3
	   ELSE
	    mk_t=bid
	   ENDIF

c P et T au temps t, par interpolation inverse en bid=m^1/3 ou m^2/3
	   CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,
	1  MIN(mk_t,m_stat_t(n_qs_t)),fe,dfe,r_stat_t,m_stat_t)
	   IF(pturb)THEN	!avec pression turbulente
	    p_t=EXP(fe(Ipg))
	   ELSE			!sans pression turbulente
	    p_t=EXP(fe(1))
	   ENDIF
	   t_t(1)=EXP(fe(2))

c ro au temps t
	   CALL etat(p_t,t_t(1),comp1,.FALSE.,
	1  ro_t(1),drop,drot,drox,u,dup,dut,dux,
	2  delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3  gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
c	   WRITE(*,2000)p_t,t_t(1),ro_t(1),mc_tmp(i)

c intégration temporelle, esti=ABS(X(t+dt)-X(t)) ie. variation par gramme
	   CALL rk_imps(t_t,ro_t,compx,t,ro,compy,dt,esti,ok,1,dm,l_nzc)
c	   WRITE(*,2000)mc_tmp(i),t,ro,compy(1),compx(1,1)
c	   PRINT*,'dans ZR',i

c recherche des variations maximales des isotopes
	   DO j=1,nchim
	    IF(esti(j) > err_max%dabon(j))THEN
	     err_max%masse(j)=mc_tmp(i)
	     err_max%abon(j)=compy(j)*nucleo(j)
	     err_max%dabon(j)=esti(j)
	     err_max%couche(j)=i
	    ENDIF
	   ENDDO
c	   PRINT*,err_max%couche ; WRITE(*,2000)err_max%dabon
c	   WRITE(*,2000)err_max%abon ; WRITE(*,2000)SQRT(mc_tmp(i))**3
c	   WRITE(*,2000)compy ; WRITE(*,2000)esti
c	   PAUSE'err_max'

c on place les abondances dans chim_zr
	   IF(ok)THEN		!il y a eu convergence
	    chim_zr(:,i)=compy

c difficulté
c err_max%masse: masse au max de variation
c err_max%abon: isotope présentant le max de variation
c err_max%dabon: variation de l'abondance
c err_max%couche: indice de la couche du max
	   ELSE
c	    PRINT* ; PRINT*,'evol sortie 2, ok=',ok
	    reprend=.TRUE.
	    ich_max=MAXLOC(err_max%dabon)
	    ich_abon_max=ich_max(1)
	    est=err_max%dabon(ich_abon_max)
	    mc_maxi=err_max%masse
	    mmax=SQRT(ABS(mc_maxi(ich_abon_max)))**3
c	    PRINT*,ok,err_max%couche,ich_abon_max,nom_elem(ich_abon_max)
c	    WRITE(*,2000)mmax,est
c	    PAUSE'mmax1'
	    SELECT CASE(langue)
	    CASE('english')
	     WRITE(*,2021)mmax,est,nom_elem(ich_abon_max),
	1    err_max%couche(ich_abon_max)
	     WRITE(2,2021)mmax,est,nom_elem(ich_abon_max),
	1    err_max%couche(ich_abon_max)
2021	     FORMAT(/,'EVOL/rk_imps, large change chem. in ZR, M/Msun=',
	1    es10.3,/,'var=',es10.3,', chemical:',a4,', shell:',i4,
	2    ', resumption dt--> dt/2',/)
	    CASE DEFAULT
	     WRITE(*,21)mmax,est,nom_elem(ich_abon_max),
	1    err_max%couche(ich_abon_max)
	     WRITE(2,21)mmax,est,nom_elem(ich_abon_max),
	1    err_max%couche(ich_abon_max)
21	     FORMAT(/,'EVOL/rk_imps, var. comp.chim. dans ZR, M/Msol=',
	1    es10.3,/,'var=',es10.3,', élément:',a4,', couche:',i4,
	2    ', reprise dt--> dt/2',/)
	    END SELECT
c	    PAUSE'mmax'
	    RETURN
	   ENDIF
	  ENDDO	!i
	 ENDDO	Bzr	!izc
c	 PAUSE'après intégration'

c on modifie le vecteur de composition chimique
c	  convf(0)=1 si convd(1) /= 1, sinon convf(0)=10000
c	  convd(nzc+1)=nc_tmp si convf(nzc) /= nc_tmp,
c	  sinon convd(nzc+1)=-10000
c	  PRINT*,nzc;  PRINT*,(convf(i),i=0,nzc)
c	  PRINT*,convd(1:nzc+1) ; PAUSE'nzc'

c redéfinition du vecteur de composition chimique nc_tmp-->n_ch
c chimd : abondances aux discontinuités, idis(ndis): indices cont.
c la seconde dim. du vecteur de comp.chim. chim(nchim,n_ch+ndis)
c tient compte des discontinuités, la dimension du vecteur nodal
c correspondant sera knotc=n_ch+m_ch+ndis, dim_ch=n_ch+ndis
c transfert des masses temporaires
c la dimension du vecteur dual mct est knotc=n_ch+m_ch+ndis
	 n_ch=nc_tmp ; DEALLOCATE(chim,mc,mct) ; ALLOCATE(mc(n_ch))
	 mc(1:n_ch)=mc_tmp(1:n_ch)

c avec no_discon=.TRUE. on ignore les discontinuités
	 Ind: IF(no_discon)THEN
	  ALLOCATE(chim(nchim,n_ch),mct(n_ch+m_ch))

c j est l'indice de la discontinuité, à la fin de la séquence j=ndis
c coté gauche on place dans chim, suivant le cas, chim_zr ou chim_zc
c coté droit on place dans chimd, suivant le cas chim_zc ou chim_zr
	  j=0
	  DO izc=1,nzc+1
c	   PRINT*,'1,izc,convf(izc-1),convd(izc)',izc,convf(izc-1),convd(izc)
	   DO i=convf(izc-1),convd(izc)	!zone radiative  ZR
	    chim(:,i)=chim_zr(:,i)
	   ENDDO

	   IF(izc <= nzc)THEN		!pour ZM
c	    PRINT*,'2,izc,convd(izc),convf(izc)',izc,convd(izc),convf(izc)
	    DO i=convd(izc),convf(izc)
	     chim(:,i)=chim_zc(:,izc)
	    ENDDO	!i=convd , convf
	   ENDIF	!izc <= nzc
	  ENDDO		!nzc

c interpolation sans discontinuité, si lisse=.TRUE. utilisation de la base duale
	  CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.FALSE.,
	1 mc(1),ll,comp2,bidd,lisse)

c dimension de l'espace d'interpolation
	  dim_ch=knotc-m_ch

c on tient compte des discontinuités, qui seront lissées si lisse=.TRUE.
	 ELSE Ind
	  ALLOCATE(chim(nchim,n_ch+ndis),mct(n_ch+m_ch+ndis))

c j est l'indice de la discontinuité, à la fin de la séquence j=ndis
c coté gauche on place dans chim, suivant le cas, chim_zr ou chim_zc
c coté droit on place dans chimd, suivant le cas chim_zc ou chim_zr
	  j=0
	  DO izc=1,nzc+1
c	   PRINT*,'izc,convf(izc-1),convd(izc),j',izc,convf(izc-1),convd(izc),j
	   DO i=convf(izc-1),convd(izc)	!zone radiative  ZR
	    IF(i == convf(izc-1) .AND. i > 1)THEN
	     IF(i == idis(j+1))THEN
	      j=j+1 ; chimd(:,j)=chim_zr(:,i)
c	      PRINT*,'ZR droite de ZM j,i,izc',j,i,izc
c	      PRINT*,'M_zc=',SQRT(ABS(mc(i)))**3/mstar
c	      PAUSE'début de ZR à droite de ZM'
	     ENDIF
	    ELSE
	     chim(:,i)=chim_zr(:,i)
	    ENDIF
	   ENDDO

	   IF(izc <= nzc)THEN		!pour ZM
c	    PRINT*,'izc,convd(izc),convf(izc),j',izc,convd(izc),convf(izc),j
	    DO i=convd(izc),convf(izc)
	     IF(i == convd(izc) .AND. i > 1)THEN
	      IF(i == idis(j+1))THEN
	       j=j+1 ; chimd(:,j)=chim_zc(:,izc)
c	       PRINT*,'ZM droite de ZR j,i,izc',j,i,izc
c	       WRITE(*,2000)chim_zc(1,izc),chimd(1,j)
c	       PAUSE'début de ZM à droite de ZR'
	      ENDIF
	     ELSE
	      chim(:,i)=chim_zc(:,izc)
	     ENDIF	!i == convd
	    ENDDO	!i=convd , convf
	   ENDIF	!izc <= nzc
	  ENDDO		!nzc

c vérification: on doit avoir j=ndis
	  IF(j /= ndis)THEN
	   PRINT*,'evol, on doit avoir: j=',j,' /= ndis=',ndis
	   PRINT*,'ARRET'
	   CALL sortie('evol 3')
	  ENDIF

c point double sur chaque discontinuité (eps: parameter défini dans mod_evol)
	  CALL bsp_dis(nchim,mc,chim,ndis,idis,chimd,eps,n_ch,m_ch,mct,knotc)
	  IF(no_croiss)THEN
	   PRINT*,'Arrêt bsp_dis dans evol' ; CALL sortie('evol 4')
	  ENDIF

c dimension de l'espace d'interpolation
	  dim_ch=knotc-m_ch

c lissage éventuel par lisse_sumd
	  IF(lisse)THEN
	   CALL lisse_sumd(npt_lisse,idis,ndis,nchim,chim,mc,mct,m_ch,
	1  n_ch,knotc)
	   IF(no_croiss)THEN
	   PRINT*,'Arrêt lisse_sumd dans evol' ; CALL sortie('evol 5')
	   ENDIF
	   dim_ch=knotc-m_ch
	  ENDIF
	 ENDIF Ind
c	 PRINT*,'après intégration : n_ch,knotc',n_ch,knotc
c	 PRINT*,'mc',mc(n_ch) ; WRITE(*,2000)mc
c	 PRINT*,'mct',mct(knotc),mc(n_ch)
c	 WRITE(*,2000)mct ; PRINT*,'xchim'
c	 DO i=1,nc
c	  CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
c	1 knotc,.TRUE.,mc(i),ll,compx,compy)
c	  WRITE(*,2000)mc(i),(compx(j)*nucleo(j),j=1,MIN(7,nchim))
c	  WRITE(*,2002)mc(i),(compx(j)*nucleo(j),j=1,nchim)
c	  WRITE(*,2002)mc(i),(compx(j),j=1,nchim)
c	 ENDDO
c	 PAUSE'après intégration'

c	 DO i=1,nc
c	  CALL inter('m',bp,q,qt,n_qs,knot,MIN(m_stat(i),m_stat(n_qs)),
c	1 fe,dfe,r_stat,m_stat)
c	  CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
c	1 knotc,.TRUE.,m_stat(i),ll,compx,compy)
c	  IF(r_stat(i) > 0.d0)PRINT*,m_stat(i),r_stat(i),,compx(1)
c	  IF(fe(3) > 0.d0)THEN
c	   PRINT*,m_stat(i),fe(3)
c	  ELSE
c	   PRINT*,m_stat(i),fe(3)
c	  ENDIF
c	 ENDDO
c	 PAUSE'après test'

c nettoyage
	 DEALLOCATE(compx,chim_zc,chim_zr,dm,ro,ro_t,t,t_t)
	ENDIF ID

c---------diffusion ou intégration et mélange classiques fin------------

c normalisation: afin d'obtenir Somme Xi=1
c on compense le fait que les nucleo i.e. les masses des éléments
c chimiques ne sont pas des nombres entiers
c pour chaque élément, dans la base des splines, on forme
c norn = somme des coefficients X nucleo, norm est proche de 1
c on divise ensuite les coefficients par norm
c ce qui est valable parce que la somme des B-splines normalisées = 1
	IF(nchim > 1)THEN
	 WHERE(chim < 1.d-50)chim=1.d-50
	 DO i=1,dim_ch	!dim_ch=knotc-m_ch: dimension de la base
	  norm=SUM(chim(1:nchim,i)*nucleo(1:nchim))
	  chim(1:nchim,i)=chim(1:nchim,i)/norm
	 ENDDO
c	 PAUSE
	ENDIF

c~~~~le vecteur nodal pour les équations différentielles de la rotation~~~~~~

c formation des points de collocation, qui serviront pour rota et poisson
c pour un ordre d'interpolation m_rot+r_qs=ord_rot, il y a m_rot points de
c collocation dans chaque intervalle, la multiplicité des noeuds étant mrot
c il y a ncoll_rot=(n_rot-1)*m_rot points de collocation et r_qs points limite
c r_qs ordre des équations différentielles
c ord_rot=m_rot+r_qs ordre d'interpolation
c avec nd discontinuités il convient d'ajouter nzc*r_qs points limite
c soit (nd+1)*r_qs points limite
c sans discontinuité interne dimension du vecteur nodal knotr=dim_rot+m_rot+r_qs
c avec discontinuités internes knotr=dim_rot+m_rot+r_qs(nzc+1)
c dim_rot, knotr,  éléments de mod_variables, définis dans base_rota
c ord_rot est initialisé dans ini_rota,
c m_rot, r_qs sont des éléments de mod_donnees
c ce n'est qu'au second passage que les modifications de masse totale
c (planétoïdes, perte de masse...) sont pris en compte

c-------------------- diffusion du moment cinétique---------------------

	SELECT CASE(Krot)
	CASE(3,4)

c interpolation de variables thermodynamiques (calcul de dro/dr)
	 CALL tab_vth

c mrot --> mrot_tmp
	 ALLOCATE(mrott_tmp(SIZE(mrott)),rota_tmp(nrot,SIZE(rota,2)))
	 rota_tmp=rota ; mrott_tmp=mrott ; knotr_tmp=knotr ; n_rot=nc_tmp

c mc --> mrot
	 DEALLOCATE(mrot) ; ALLOCATE(mrot(n_rot))
	 mrot=mc_tmp(1:n_rot) ; mrot(1)=0.d0

c génération de la base de de Boor pour la rotation
	 CALL base_rota

c spline sur le nouveau vecteur nodal mrott
	 DEALLOCATE(rota) ;  ALLOCATE(rota(nrot,dim_rot))

c	 CALL newspl(nrot,mrot,mrott_tmp,knotr_tmp,ord_rot,mrott,knotr,ord_rot,
c	1 rota_tmp,rota,.TRUE.)
	 CALL newspl_gal(nrot,mrot,mrott_tmp,knotr_tmp,ord_rot,rota_tmp,
	1 mrott,knotr,ord_rot,rota)
	 IF(no_croiss)PRINT*,'Pb. en 2 dans evol'

c les points de collocation
	  ncoll_rot=(n_rot-1)*m_rot
	  DEALLOCATE(xcoll_rot) ; ALLOCATE(xcoll_rot(ncoll_rot))
	  CALL coll(mrot,n_rot,m_rot,xcoll_rot)

c nettoyage
	  DEALLOCATE(mrott_tmp,rota_tmp)

c intégration
	 PRINT*
	 SELECT CASE(langue)
	 CASE('english')
	  WRITE(usl_evol,*)'---- diffusion of angular momentun (begin) ----'
	 CASE DEFAULT
	  WRITE(usl_evol,*)'---- diffusion du moment cinétique (début) ----'
	 END SELECT

c résolution des équations
	CALL resout_rota(dt,ok)
	reprend=.NOT.ok ; IF(reprend)RETURN

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

c écriture des variables de la rotation off line
	 IF(Kdes_rot == 3)CALL ecrit_rota(dt)
	 SELECT CASE(langue)
	 CASE('english')
	  WRITE(usl_evol,*)'---- diffusion of angular momentun (end) ----'
	 CASE DEFAULT
	  WRITE(usl_evol,*)'---- diffusion du moment cinétique (fin) ----'
	 END SELECT

c ~~~~~~~~~~~conservation locale du moment cinétique~~~~~~~~~~~~~~~~

	CASE(5)

c valeur moyenne sur la couche: omega=r omega_t/(3r -2r_t)
c on se place aux masses (m^2/3) de la composition chimique
	 DEALLOCATE(mrot,mrott,rota) ; n_rot=n_ch
	 ALLOCATE(mrot(n_rot),mrott(n_rot+ord_rot),rota(nrot,n_rot))
	 mrot=mc

c rayon au temps t+dt
	 DO i=2,n_rot
	  IF(en_m23)THEN
	   bid=mrot(i)			!m^2/3
	  ELSE
	   bid=SQRT(ABS(mrot(i)))	!m^1/3
	  ENDIF
	  CALL inter('m',bp,q,qt,n_qs,knot,MIN(bid,m_stat(n_qs)),fe,dfe,
	1 r_stat,m_stat)
	  IF(en_m23)THEN
	   r=SQRT(fe(3))
	  ELSE
	   r=fe(3)
	  ENDIF

c rayon au temps t, masse au temps t s'il y a perte de masse
	  IF(l_pertm .OR. l_pertmt)THEN
	   bid=MIN(bid,x_ptm(n_ptm))
	   CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,
	1  .TRUE.,bid,ll,f,dfdx)		!masse au temps t
	   IF(no_croiss)PRINT*,'Pb. at 14 in evol'
	   IF(.NOT.en_m23)f(1)=ABS(f(1))**2	!m^2/3
	  ELSE
	   f(1)=mrot(i)				!m^2/3
	  ENDIF
	  mk=SQRT(f(1))**3  			!m

c omega_t au temps t
	  mk_t=MIN(MAX(mrot_t(1),f(1)),mrot_t(n_ch_t))
	  CALL bsp1dn(nrot,rota_t,mrot_t,mrott_t,n_rot_t,ord_rot,knotr_t,
	1 .TRUE.,mk_t,ll,fe,dfe)
	  IF(no_croiss)PRINT*,'Pb. at 15 in evol'
	  omega_t(1)=fe(1)

c r_t au temps t
	  CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,
	1 MIN(bid,m_stat_t(n_qs_t)),fe,dfe,r_stat_t,m_stat_t)
	  IF(en_m23)THEN
	   r_t=SQRT(fe(3))
	  ELSE
	   r_t=fe(3)
	  ENDIF

c omega(t+dt)
	  rota(1,i)=r*Omega_t(1)/(3.d0*r-2.d0*r_t)
	  r=r*rsol ; mk=mk*msol
	  w_max=SQRT(g*mk/r**3)/2.d0
	  IF(rota(1,i) > w_max)THEN
	   WRITE(*,2000)rota(1,i), w_max,r_t
	   PAUSE'Omega > w_max'
	  ENDIF
	 ENDDO
	 rota(1,1)=rota(1,2)

c tabulation nouvel omega
	 CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.FALSE.,
	1 mrot(1),ll,fe,dfe)
	 dim_rot=knotr-ord_rot
	END SELECT	!Krot

c mc_tmp est désormais inutile
	DEALLOCATE(mc_tmp)

c---------------------estimation du pas temporel suivant--------------

c extractions
	ich_max=MAXLOC(err_max%dabon)
	ich_abon_max=ich_max(1)
	est=err_max%dabon(ich_abon_max)
	kmax=err_max%couche(ich_abon_max)
	mc_maxi=err_max%masse
	mc_maxi=SQRT(ABS(mc_maxi))**3
	ab_max=err_max%abon
	estim=err_max%dabon
	kmaxi=err_max%couche
c	WRITE(*,2000)err_max%dabon ; WRITE(*,2000)mc_maxi ; PRINT*,kmaxi
c	PAUSE'extractions'

c nouveau dt
c	WRITE(*,2000)dt,precit,est,precit/est
	IF(est /= 0.d0)THEN
	 IF(diffusion)THEN
	  dtnew=dt*precit/est
	 ELSE
	  dtnew=dt*(precit/est)**(1.d0/(ordre+1))
	 ENDIF
	ELSE
	 dtnew=1.2d0*dt
	ENDIF

c limitation de la variation temporelle de la dégénérescence
c variation relative maximale
	dege_max=-100.d0
	B3: DO i=1,n_qs
	 p=EXP(bp(1,m_qs*(i-1)+1))
	 p_t=EXP(bp(2,m_qs*(i-1)+1)) !p_t:vT
	 mk=bp(5,m_qs*(i-1)+1)	!m^1/3 ou m^2/3
	 IF(.NOT.en_m23)mk=mk**2	!m^2/3

c composition chimique au temps t+dt
	 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
	1 knotc,.TRUE.,mk,ll,comp,comp2)
	 comp1=comp*nucleo	!par gramme

c ro au temps t+dt
	 CALL etat(p,p_t,comp1,.FALSE.,	!ro: a:vT
	1 a,drop,drot,drox,u,dup,dut,dux,
	2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
c	 WRITE(*,2000)p,t(k),ro(k),mk

c dégénérescence au temps t+dt
	 dege=dege_elec(comp,p_t,a)	!dege_elec(X,T,ro) p_t=T :vT, ro=a:vT
	 dege_max=MAX(dege,dege_max)
	 IF(dege <= 0.d0)EXIT B3
	ENDDO B3
c	PRINT*,'dege' ; WRITE(*,2000)dege

c pas temporel estimé, entre 0.8 dt et 1.2 dt pour éviter trop de variations
	SELECT CASE(Krot)
	CASE(3,4)
	 dtn=MAX(0.8d0*dt,MIN(dtnew,1.1d0*dt))
	CASE DEFAULT
	 IF(cephe)THEN
c	  dtn=MAX(0.5d0*dt,MIN(dtnew,1.2d0*dt))
	  dtn=MIN(dtnew,1.2d0*dt)
	 ELSE
	  dtn=MAX(0.8d0*dt,MIN(dtnew,1.2d0*dt))
	 ENDIF
	END SELECT
c	PRINT*,'dt,dtn,dtnew' ; WRITE(*,2000)dt,dtn,dtnew

c limitation du pas temporel en cas de dégénérescence moyenne
	IF(dege_max > 8.d0 .AND.
	1 entre(ln0t8,ln1t8,bp(2,1)))dtnew=MIN(5.d-4,dtnew)

c limitation du pas temporel au voisinage de l'amorçage du 3 alpha
	loop = ihe4 > 0.d0 .AND. entre(ln1t8,ln2t8,bp(2,1))

c ou pour ne pas rater les blue loops
	IF(loop)loop = chim(ihe4,1)*nucleo(ihe4) > 0.1d0
        IF(loop)THEN
         dtn=MIN(5.d-2,dtn)

c en cas de dégénérescence forte
c	 IF(dege_max > 10.d0)dtn=MIN(dtn,1.d-6)
        ENDIF

c écritures
	WRITE(usl_evol,111)kmax,est,nom_elem(ich_abon_max),dege_max,dt,dtn
111	FORMAT(/,'EVOL couche ',i4,', var. max.:',es10.3,' pour ',a4,
	1 ', dégé. élec. max:',es10.3,/,'dt=',es10.3,', dt optimal=',es10.3,
	2 ', masse(Msol)/abondance/variation/couche:')
	WRITE(usl_evol,114)nom_elem(1:MIN(10,nchim))
114	FORMAT(10(2x,a4,2x))
	WRITE(usl_evol,112)mc_maxi(1:MIN(10,nchim))
	WRITE(usl_evol,112)(ab_max(i),i=1,MIN(10,nchim))
	WRITE(usl_evol,112)estim(1:MIN(10,nchim))
112	FORMAT(10es8.1)
	WRITE(usl_evol,113)kmaxi(1:MIN(10,nchim))
113	FORMAT(10(2x,i4,2x))

c pas temporel estimé
	dts=MIN(dtmax,dtn,agemax-age-dt) ; reprend=.FALSE.
c	PRINT*,'fin evol dt,dts,dtn,dtmax,agemax,age,agemax-age-dt,Krot=',Krot
c	WRITE(*,2000)dt,dts,dtn,dtmax,agemax,age,agemax-age-dt
c	PAUSE'fin de evol'

	RETURN

	CONTAINS
	 INCLUDE 'base_rota.f'

	END SUBROUTINE evol