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

	SUBROUTINE pertm(dt)

c routine générique pour le calcul de la perte/gain de masse
c routine PRIVATE du module mod_static

c la perte de masse solaire due uniquement à la transformation par les réactions
c nucléaires de H en He4, est de l'ordre de 7d-14Msol / an soit 3.5d-4 sur 5Gy
c quantité négligeable que si on recherche l'effet d'une perte/gain de masse
c de valeur absolue supérieure à 1d-12 Msol/an. Cette perte de masse affecte
c essentiellement le noyau où les réactions sont actives: soit environ le
c 1/3 plus interne. Avec l'allocation automatique des points de masse
c on ne connait pas la masse m(t) qui donnera m(t+dt) qui, seule, est connue.
c Pour la connaissance des quantités au temps t permettant le calcul des dérivées
c temporelles (TdS/dt et dchim/dt) et pour celle de la composition chimique
c on doit faire le calcul à rebour interpolation m(t+dt)--->m(t) en tenant compte
c des variations de masse dues a E=mc**2 et à une éventuelle autre
c perte/gain de masse (mdot > 0 : gain de masse, mdot < 0 : perte de masse)
c la perte de masse externe (c.a.d. classique) est supposée concentrée
c dans la couche [n_ptm-1,n_ptm], elle sera prise en compte dans l'atmosphère
c par le traitement de la condition limite externe, la masse totale mstar,
c tenant compte de la perte de masse restant celle au rayon bolométrique

c on entre avec mstar=masse au temps t
c on sort  avec mstar=masse au temps t+dt

c on intègre les variations dues a E=mc**2
c d'où la variation de masse totale et la nouvelle masse totale mstar
c calcul de la masse qui correspondait, au temps t, à la masse en chaque point
c tabule l'ancienne masse en fonction de la nouvelle (en m**2/3 ou m^1/3)
c on néglige la variation de masse de l'atmosphère entre les instants t+dt et t

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

c 25 08 97 : mise en place des variables eulériennes
c 30 07 00 : introduction F95
c 02 08 12 : introduction systématique de la perte de masse due à e=mc^2
c	     rassemblement de tous les cas de perte de masse externe

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

c entrées
c	bp,q,n_qs,qt,knot,chim,mc,nc,mct,knotc : modele au temps t
c	dt : pas temporel

c entrées/sorties
c	mstar : masse totale avec perte de masse
c	en entrée au temps t, en sortie au temps t+dt

c sorties
c	old_ptm,x_ptm,xt_ptm,n_ptm,knot_ptm : interpolation de
c	l'ancienne masse en fonction de la nouvelle masse en m^2/3 ou m^1/3
C
c----------------------------------------------------------------

	USE mod_donnees, ONLY : clight, en_m23, eve, l_pertm, l_pertmt, msol,
	1 nchim, mdot, m_ch, m_ptm, m_qs, nchim, ne, nom_pertm, pturb,
	2 ord_qs, secon6, t_inf
	USE mod_etat, ONLY : etat
	USE mod_kind
	USE mod_nuc, ONLY : l_planet, nuc, planetoides
	USE mod_numerique, ONLY : bsp1dn, no_croiss, sum_n
	USE mod_variables, ONLY : bp, chim, chim_gram, knot, knotc,
	1 knot_ptm, log10_teff, lstar, mc, mct, mstar, mstar_t, n_ch,
	2 n_ptm, n_qs, old_ptm, q, qt, rstar, sortie, x_ptm, xt_ptm

	IMPLICIT NONE

	REAL (kind=dp), INTENT(in) :: dt

	REAL (kind=dp), DIMENSION(nchim,nchim) :: jac
	REAL (kind=dp), DIMENSION(1,n_qs) :: dm
	REAL (kind=dp), DIMENSION(ne) :: df, f
	REAL (kind=dp), DIMENSION(nchim) :: dxchim, depsx, xchim, xchimm
	REAL (kind=dp), DIMENSION(n_qs+6) :: mt
	REAL (kind=dp), DIMENSION(n_qs) :: m, vt1, vt2
	REAL (kind=dp), SAVE, DIMENSION(5) :: epsilon

	REAL (kind=dp), SAVE :: cte1, ehh, ebe7, eb8, ef17, en13, eo15, mev
	REAL (kind=dp) :: alfa, beta, be7, bid, b8, cp, cte2, c1, c2, c3, dcpp,
	1 dcpt, dcpx, def, delta, deltap, deltat, deltax, depsro, depst,
	2 dmass, dplanet, drop, drot, drox, dup, dut, dux, gamma1, gradad, dgradadp,
	3 dgradadt, dgradadx, hh, m_new, m_planet, n13, o15, f17, p, ro, t, u

	INTEGER :: i, l=1

	LOGICAL, SAVE :: init=.TRUE.

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

2000	FORMAT(8es10.3)

c	PRINT*,ALLOCATED(x_ptm) ; PAUSE'entrée pertm'

c initialisations
	dmass=0.d0 ; def=0.d0 ; dplanet=0.d0

c si on tient compte des pertes de masse dues aux réactions thermonucléaires
	IF(l_pertmt)THEN
	 IF(init)THEN
	  init=.FALSE.
	  cte1=secon6/clight/clight
	  mev=1.d6*eve	!énergie des neutrinos selon Clayton p.380 et 392
	  ehh=mev*0.263d0 ; ebe7=mev*0.80d0 ; eb8=mev*7.2d0
	  en13=mev*0.71d0 ; eo15=mev*1.d0 ; ef17=mev*0.94d0
	  WRITE(*,1) ; WRITE(2,1)
1	  FORMAT(/,'on tient compte de la perte de Masse due à E=mc**2',/)
	 ENDIF

c e=mc^2, d masse (Msol) = énergie rayonnée/c^2 pendant dt (My), dm=e/c^2 * dt / Msol
	 cte2=cte1*dt

c extraction des masses vt2 : m^1/3 ou m^2/3, m en Msol
c calcul des réactions thermonucléaires puis des pertes par neutrinos
c m(i)=delta m dû à e=mc2 pendant dt en fraction de Mstar x Msol
	 dm=0.d0 ; m=0.d0
	 DO i=2,n_qs
	  CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),l,f,df)
	  IF(no_croiss)PRINT*,'Pb. at 1 in pertm'
c	  WRITE(*,2000)f
	  IF(pturb)THEN		!avec pression turbulente 7 inconnues
	   p=EXP(f(7))
	  ELSE			!sans pression turbulente 6 inconnues
	   p=EXP(f(1))
	  ENDIF
	  t=EXP(f(2))
	  vt2(i)=f(5)	!masses en m^2/3 ou m^1/3
	  IF(en_m23)THEN	!bid=m^2/3 pour mc
	   bid=ABS(f(5))
	  ELSE
	   bid=f(5)**2
	  ENDIF
	  m(i)=SQRT(bid)**3	!masse en Msol au temps t+dt
c	  WRITE(*,2000)m(i),f(5)
	  CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
	1 knotc,.TRUE.,MIN(bid,mc(n_ch)),l,xchim,dxchim)
	  IF(no_croiss)PRINT*,'Pb. at 2 in pertm'
	  IF(t > t_inf)THEN
	   xchimm=xchim
	   CALL chim_gram(xchimm,dxchim)
	   CALL etat(p,t,xchimm,.FALSE.,
	1  ro,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 perte de masse locale en grammes due aux réactions
	   CALL nuc(t,ro,xchim,dxchim,jac,.FALSE.,3,	!3 : epsilon
	1  epsilon,depst,depsro,depsx,hh,be7,b8,n13,o15,f17)

c hh,be7,b8,n13,o15,f17 : nombre de neutrinos par unité de masse
c et de temps pour les diverses réactions
c perte de masse locale due neutrinos en m(i)
	   CALL nuc(t,ro,xchim,dxchim,jac,.FALSE.,4,	!4 : neutrinos
	1  epsilon,depst,depsro,depsx,hh,be7,b8,n13,o15,f17)

c perte de masse locale totale (Msol) pendant dt (My) au point de masse m (Msol)
	   dm(1,i)=(epsilon(1)+hh*ehh+be7*ebe7+b8*eb8+n13*en13+
	1  o15*eo15+f17*ef17)*cte2
	  ENDIF
c	  WRITE(*,2000)m(i),f(5),dm(1,i)
	 ENDDO
c	 WRITE(*,2000)dt,cte2,mstar_t ; WRITE(*,2000)dm(1,:)
c	 WRITE(*,2000)m

c suppression des inversions en masse, elles sont anormales avec perte
c de masse mais pouvent se produire avec gain de masse
	 n_ptm=n_qs ; i=n_ptm-1
	 DO WHILE(i > 1)
	  IF(m(i) >= m(i+1))THEN
	   n_ptm=n_ptm-1
c	   PRINT*,'perte_tot, inversion en',i,n_ptm,m(i),m(i+1) ; PAUSE
	   DO l=i,n_ptm
	    m(l)=m(l+1) ; dm(1,l)=dm(1,l+1) ; vt2(l)=vt2(l+1)
	   ENDDO
	   i=MIN(i,n_ptm-1)
	  ELSE
	   i=i-1
	  ENDIF
	 ENDDO

c spline des pertes dm (Msol) en fonction de m (Msol) pour intégration
	 CALL bsp1dn(1,dm,m,mt,n_ptm,m_ptm,knot_ptm,.FALSE.,m(1),l,xchim,
	1 dxchim)
         IF(no_croiss)THEN
          CALL sortie('pertm 1')
         ENDIF

c initialisation de l'intégrale sum_n, des pertes de masse aux points de masse m
	 CALL sum_n(1,dm,mt,m_ptm,knot_ptm,.FALSE.,m(1),m(n_ptm),f)

c calcul de l'intégrale def, des pertes de masse en tous les points de masse m
c vt1 = masse au temps t avant perte de masse vt1-->old_ptm masse au temps t+dt
	 DO i=1,n_ptm
	  CALL sum_n(1,dm,mt,m_ptm,knot_ptm,.TRUE.,m(1),m(i),f)
	  def=f(1)
	  vt1(i)=m(i)-def
c	  PRINT*,i ; WRITE(*,2000)def,f(1)
	 ENDDO
c         cte2=secon6/clight/clight/msol*dt

c si n_ptm < n_qs, la boucle suivante est effective
c	 DO i=n_ptm+1,n_qs
c	  vt1(i)=m(i)-def
c	 ENDDO

c vt1 en m**1/3 ou m^2/3, avec perte de masse due à e=mc^2
	 vt1=vt1**(1.d0/3.d0)
	 IF(en_m23)vt1=vt1**2

c suppression des masses négatives anormales vt2 en m**2/3 ou m^1/3
	 vt1(1)=0.d0 ; vt2(1)=0.d0 ; i=2
	 DO WHILE(i <= n_ptm)
	  IF(vt1(i) <= 0.d0 .OR. vt2(i) <= 0.d0)THEN
	   n_ptm=n_ptm-1
	   DO l=i,n_ptm
	    vt1(l)=vt1(l+1) ; vt2(l)=vt2(l+1)
	   ENDDO
	  ELSE
	   i=i+1
	  ENDIF
	 ENDDO

c suppression des inversions, elles sont anormales avec
c perte de masse mais peuvent se produire avec gain de masse
	 i=n_ptm-1
	 DO WHILE(i > 1)
	  IF(vt1(i) >= vt1(i+1) .OR. vt2(i) >= vt2(i+1))THEN
c	   PRINT*,'perte_ext, inversion en',i,n_ptm,
c	1  vt1(i),vt1(i+1),vt2(i),vt2(i+1)
c	   PAUSE
	   n_ptm=n_ptm-1
	   DO l=i,n_ptm
	    vt1(l)=vt1(l+1) ; vt2(l)=vt2(l+1)
	   ENDDO
	   i=MIN(i,n_ptm-1)
	  ELSE
	   i=i-1
	  ENDIF
	 ENDDO

c allocations
	 IF(ALLOCATED(old_ptm))DEALLOCATE(old_ptm,x_ptm,xt_ptm)
	 ALLOCATE(old_ptm(1,n_ptm),x_ptm(n_ptm),xt_ptm(n_ptm+m_ptm))
	 old_ptm(1,1:n_ptm)=vt1(1:n_ptm) ; x_ptm(1:n_ptm)=vt2(1:n_ptm)

c old_ptm=vt1 est la masse (en m^1/3 ou m^2/3) au temps t qui deviendra
c bp(5,:)=vt2 (en m^1/3 ou m^2/3) au temps t+dt compte tenu de la perte de masse
c la composition chimique au temps t est celle en vt1

c tabulation de old_ptm en fonction de x_ptm
	 CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,.FALSE.,
	1 x_ptm(1),l,f,df)
	 IF(no_croiss)THEN
	  CALL sortie('pertm 2')
	 ENDIF
c	 PRINT*,mstar,dt ; PAUSE'sortie pertm'

c détection d'une inversion de masse dans pertm
	 IF(n_ptm /= n_qs)THEN
	  WRITE(*,2)n_qs,n_ptm
2	  FORMAT("inversion de masse dans pertm n_qs=",i4,'/= n_ptm=',i4)
	 ENDIF

c on ne tient pas compte des pertes de masse dues à E=mc^2
	ELSEIF(l_pertm)THEN

c extraction des masses, suppression des inversions au temps t+dt en m^2/3 ou m^1/3
	 DO i=1,n_qs
	  vt1(i)=bp(5,m_qs*(i-1)+1)
	 ENDDO

c on s'assure de la stricte croissance
	 n_ptm=1 ; vt2(n_ptm)=vt1(1)
	 DO i=2,n_qs
	  IF(vt1(i) > vt2(n_ptm))THEN
	   n_ptm=n_ptm+1 ; vt2(n_ptm)=vt1(i)
	  ENDIF
	 ENDDO

c allocations
	 IF(ALLOCATED(old_ptm))DEALLOCATE(old_ptm,x_ptm,xt_ptm)
	 ALLOCATE(old_ptm(1,n_ptm),x_ptm(n_ptm),xt_ptm(n_ptm+m_ptm))
	 x_ptm=vt2(1:n_ptm) ; old_ptm(1,1:n_ptm)=vt2(1:n_ptm)

c tabulation de la répartition de masse sans inversion
	 CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,.FALSE.,
	1 x_ptm(1),l,f,df)
	 IF(no_croiss)THEN
	  CALL sortie('pertm 3')
	 ENDIF
	ENDIF

c avec perte/gain de masse externe
	IF(ABS(mdot) /= 0.d0)THEN

c cas prévus de perte de masse externe
	 SELECT CASE(nom_pertm)

c avec du vent et no_pertm on peut avoir mdot /= 0 par erreur
	 CASE ('no_pertm')
	  dmass=0.d0

c perte de masse "stationary dust-driven winds" Dominik & al. 1990 A&A, 240, 365
	 CASE ('pertm_dominik')
	  c3=LOG10(lstar)-10.d0**LOG10_teff/600.d0
	  IF(c3 > 0.d0)THEN
	   c1=-3.7d0-0.93d0*mstar
	   c2=1.71d0+2.59d0*mstar
	   dmass=10.d0**c1*c3**c2*1.d6*dt
	  ELSE
	   dmass=0.d0
	  ENDIF

c perte/gain de masse externe, selon signe de mdot en Msol/an
	 CASE ('pertm_ext')
	  dmass=mdot*1.d6*dt

c compte tenu des pertes/gains externes, la masse totale ne pourra pas franchir Msol
	 CASE ('pertm_msol')
	  IF(mstar_t == 1.d0)THEN
	   m_new=1.d0
	  ELSE
	   m_new=mstar_t+mdot*1.d6*dt-def
	   IF(mdot < 0.d0 .AND. m_new < 1.d0)THEN
	    m_new=1.d0
	   ELSEIF(mdot > 0.d0 .AND. m_new > 1.d0)THEN
	    m_new=1.d0
	   ENDIF
	  ENDIF

c nouvelle masse
          mstar=m_new
          RETURN

c perte de masse solaire empirique de Reimers
c cité dans l'article Rasio & al.ApJ 470, 1187, 1996 (taux=0.6d-13 Msol/y)
	 CASE ('pertm_reimers')
	  dmass=4.d6*mdot*lstar*rstar/mstar*dt

c perte de masse empirique de Reimers si > perte de masse de Dominik
	 CASE ('pertm_reim-domi')
	  dmass=4.d6*mdot*lstar*rstar/mstar*dt
	  c3=LOG10(lstar)-10.d0**LOG10_teff/600.d0
	  IF(c3 > 0.d0)THEN
	   c1=-3.7d0-0.93d0*mstar
	   c2=1.71d0+2.59d0*mstar
	   dmass=MAX(dmass,10.d0**c1*c3**c2*1.d6*dt)
	  ENDIF

c perte de masse empirique de Waldron
c cité dans l'article A.Ap 229 (1990) 469-474
	 CASE ('pertm_waldron')
	  dmass=10.d0**(1.07d0*LOG10(lstar)+1.77d0*LOG10(rstar)-14.3d0)*1.d6*dt

	 CASE DEFAULT
	  PRINT*,'routine de perte/gain de masse inconnue: ',nom_pertm
	  PRINT*,'routines connues: pertm_ext, pertm_msol'
	  PRINT*,'pertm_reimers, pertm_waldron'
	  PRINT*,'arrêt' ; CALL sortie('pertm 4')
	 ENDSELECT
 	ENDIF

c contribution des planétoïdes, m_planet en Msol/My
	IF(l_planet)THEN
	 CALL planetoides(m_planet=m_planet)
	 dplanet=m_planet*dt
	ENDIF

c bilan total mstar(t+dt)=mstar(t)+m_dot*dt-defaut de masse + planet
	mstar=mstar_t+dmass-def+dplanet
c	PRINT*,'mstar,mstar_t',mstar,mstar_t,MAXVAL(old_ptm),MAXVAL(x_ptm)
c	WRITE(*,2000)m_new,mstar,mstar_t,mdot,dmass,def
c	PAUSE'pertm'

	RETURN

	END SUBROUTINE pertm