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

	SUBROUTINE rq_reac(comp,t,ro,r,drt,dro,drx,q,dqt,dqo,dqx,mue,dmuex,ecranw)

c routine private du module mod_nuc

c calcul des R et Q des réactions thermonucléaires
c et de l'effet d'écran, on utilisera la tabulation tabul_nuc
c au premier appel nchim est connu

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

c	effet d'écran faible ou fort de Mitler ApJ 212, 513, 1977
c	mitler=.TRUE. : écrantage de Mitler, sinon écran faible
c	exploitation F95 18.10.95

c entrées :
c	comp : abondances en mole
c	t : température
c	ro : densité

c sorties :
c	r(i) : taux de réaction g /s /mole/mole (/mole)
c	drt, dro, drx : dérivées
c	q(i): énergie ( erg /sec /gr )
c	dqt, dqo, dqx, dmuex : dérivées
c	mue : 1./masse moléculaire moyenne par électron
c	ecranw : écrantage, sortie OPTIONAL pour tests

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

	USE mod_donnees, ONLY : amu, echarg, eve, i_ex, kbol, hpl, me, mitler,
	1 nchim, nucleo, pi, t_inf, zi
	USE mod_numerique, ONLY : bsp1dn
	USE mod_kind
	USE mod_variables, ONLY : sortie

	IMPLICIT NONE

	REAL (kind=dp), INTENT(in), DIMENSION(:) :: comp
	REAL (kind=dp), INTENT(in) :: t, ro
	REAL (kind=dp), INTENT(out), DIMENSION(:,:) :: drx, dqx
	REAL (kind=dp), INTENT(out), OPTIONAL, DIMENSION(nreac) :: ecranw
	REAL (kind=dp), INTENT(out), DIMENSION(:) :: r, drt, dro, q,
	1 dqt, dqo, dmuex
	REAL (kind=dp), INTENT(out) :: mue

	REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: decrx, dzetax

	REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: decro, decrt, dr,
	1 ecran, lambda, mz1z2, z1z2, zz53

	REAL (kind=dp), DIMENSION(nchim) :: abon, datx, dbidx, dc1x, dfx, df12x,
	1 dkx, dnex, dzstarx

	REAL (kind=dp), SAVE :: p53=5.d0/3.d0
	REAL (kind=dp), SAVE, DIMENSION(2) :: zeta, dzetat, dzetaro
	REAL (kind=dp), SAVE :: cte1, cte2, cte3, cte4, cte5, cte6, cte7,
	1 cte8, cte9, cte10
	REAL (kind=dp) :: zstar, ne, c1, f12, at, k, bid, bid1, df,
	1 e0, eps, beta, dnero, dc1ro, df12t, df12ro, datt, datro, dkt,
	2 dkro, dbidt, dbidro, dft, dfro

	INTEGER, SAVE :: l=1
	INTEGER :: i, j

	LOGICAL, SAVE :: fort=.TRUE., init=.TRUE.

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

2000	FORMAT(8es10.3)
2001	FORMAT(10es8.1)
2002	FORMAT(20i4)

	IF(init)THEN	!en cas de modifications assurer la cohérence
	 init=.FALSE.	!avec tabul_nuc

c	 PRINT*,nchim ; PAUSE'rq_reac entrée'

c tabulation des réactions nucléaires
	 CALL tabul_nuc

c écritures diverses
	 IF(mitler)THEN
	  WRITE(2,4) ; WRITE(*,4)
4         FORMAT(/,'écrantage de Mitler')
	 ELSE
	  WRITE(2,5) ; WRITE(*,5)
5         FORMAT(/,'écrantage faible ou fort')
	 ENDIF
c	 WRITE(*,2000)q0 ; PAUSE'les q'

	 WRITE(2,1)t_inf ; WRITE(*,1)t_inf
1	 FORMAT('Température d''amorçage des réactions nucléaires:',
	1 es10.3)
 	 WRITE(2,2)t_sup ; WRITE(*,2)t_sup
2	 FORMAT('Température maximale pour les réactions nucléaires:',
	1 es10.3)

c quelques constantes
	 cte1=2.d0/3.d0*pi*echarg ; cte2=4.d0*pi*echarg**2/kbol
	 cte3=hpl**3/4.d0/pi/SQRT(2.d0*me*kbol)**3*3.d0/2.d0
	 cte4=3.d0/4.d0/pi ; cte5=-18.d0/5.d0
	 cte6=1.2204d0*(1.d-6)**(2.d0/3.d0)*eve*1.d3	!3', Lang 4-44 A en MeV
	 cte7=5.d0/3.d0 ; cte8=0.001037d0
	 cte9=7.d0*5.47d2/9.6d2 ; cte10=5.d0/8.d0

c quantités à conserver
	 ALLOCATE(z1z2(nreac),lambda(nreac),mz1z2(nreac),zz53(nreac))
 	 DO i=1,nreac
	  z1z2(i)=izz(i,1)*izz(i,2)
	  zz53(i)=(izz(i,1)+izz(i,2))**p53-izz(i,1)**p53-izz(i,2)**p53
	  mz1z2(i)=cte6*(z1z2(i)**2*ar(i))**(1.d0/3.d0)	!Lang 4-44
	  lambda(i)=z1z2(i)*echarg**2		!eq. 30'
	 ENDDO

c autres allocations
	 ALLOCATE(decro(nreac),decrt(nreac),decrx(nreac,nchim),
	1 dr(nreac),dzetax(2,nchim),ecran(nreac))
c	 PAUSE'sortie init'

	 RETURN	!le premier appel est un appel d'initialisation

	ENDIF

	IF(t < t_inf)THEN
	 r=0.d0 ; q=0.d0 ; drt=0.d0 ; dro=0.d0 ; dqt=0.d0 ; dqo=0.d0
	 RETURN
	ELSEIF(t > t_sup)THEN
	 WRITE(*,3)t,t_sup
3	 FORMAT('T=',es10.3,' > T_sup='es10.3,
	1', hors de la tabulation des réactions nucléaires')
	 WRITE(2,3)t,t_sup ; CALL sortie('rq_reac')	!STOP
	ENDIF

c comp --> abon
	abon=comp ; WHERE(abon < 0.d0)abon=0.d0

c effet d'écran Mitler
	Im: IF(mitler)THEN
	 mue=DOT_PRODUCT(zi,abon)	!1/poids moléculaire par électron
	 dmuex=zi
         dzstarx=dmuex**2
         zstar=DOT_PRODUCT(dzstarx,abon)
	 zstar=zstar/mue		!eq. 6"
	 dnero=mue/amu
	 ne=ro*dnero	!nb. électrons, mue=1/poids mol. par électron
	 c1=cte1*ne				!eq. 17'
	 dc1ro=cte1*dnero
	 DO j=1,nchim
	  dzstarx(j)=dzstarx(j)/mue-zstar/mue*dmuex(j)
	 ENDDO
	 dnex=ne/mue*dmuex
	 dc1x=cte1*dnex
	 f12=cte3*ne/SQRT(t)**3		!Clayton eq. 2.59 X 3/2, Kipp 15.39
	 df12t=-3.d0/2.d0*f12/t
	 df12ro=f12/ne*dnero
	 df12x=f12*dnex/ne
	 at=4.d0/9.d0*f12**(4.d0/3.d0)		!eq. 22"
	 datt= 4.d0/3.d0*at/f12*df12t
	 datro=4.d0/3.d0*at/f12*df12ro
	 datx=4.d0/3.d0*at/f12*df12x
	 at=1.d0/SQRT(at+1.d0)			!eq. 22"
	 datt=  -datt/2.d0*at**3		!eq. 22"
	 datro=-datro/2.d0*at**3		!eq. 22"
         datx=-datx/2.d0/at**3
	 k=cte2*ne/t*(zstar+at)			!eq 22'
	 dkt= k*(datt /(zstar+at)-1.d0/t)
	 dkro=k*(datro/(zstar+at)+dnero/ne)
	 DO j=1,nchim
	  dkx(j)=k*((dzstarx(j)+datx(j))/(zstar+at)+dnex(j)/ne)
	 ENDDO
	 k=SQRT(k)				!eq. 22'
	 dkt=dkt/2.d0/k ; dkro=dkro/2.d0/k
	 dkx=dkx/2.d0/k

c	 PRINT*,'zstar,mue,ne,c1,f12,at,k'
c	 WRITE(*,2000)zstar,mue,ne,c1,f12,at,k

	 B1: DO i=1,nreac
	  IF(z1z2(i) /= 0.d0)THEN	!pour passer les captures e-
	   DO j=1,2
	    zeta(j)=cte4*k**3/ne*izz(i,j)	!eq. 24
	    dzetat(j)= zeta(j) *3./k*dkt
	    dzetaro(j)=zeta(j)*(3./k*dkro-dnero/ne)
	    DO l=1,nchim
	     dzetax(j,l)=zeta(j)*(3./k*dkx(l)-dnex(l)/ne)
	    ENDDO
	   ENDDO
	   bid1=(zeta(1)+zeta(2)+1.d0)**cte7		!eq. 25 partie
	   bid=bid1+1.d0
	   dbidt= bid1*cte7/(zeta(1)+zeta(2)+1.d0)*(dzetat(1) +dzetat(2))
	   dbidro=bid1*cte7/(zeta(1)+zeta(2)+1.d0)*(dzetaro(1)+dzetaro(2))
	   DO l=1,nchim
	    dbidx(l)=bid1*cte7/(zeta(1)+zeta(2)+1.d0)*(dzetax(1,l)
	1   +dzetax(2,l))
	   ENDDO
	   bid1=-(zeta(1)+1.d0)**cte7 		!eq. 25 partie
	   bid=bid+bid1
	   dbidt= dbidt +bid1*cte7/(zeta(1)+1.d0)*dzetat(1)
	   dbidro=dbidro+bid1*cte7/(zeta(1)+1.d0)*dzetaro(1)
	   DO l=1,nchim
	    dbidx(l)=dbidx(l)+bid1*cte7/(zeta(1)+1.d0)*dzetax(1,l)
	   ENDDO
	   bid1=-(zeta(2)+1.d0)**cte7 		!eq. 25 partie
	   bid=bid+bid1
	   dbidt= dbidt +bid1*cte7/(zeta(2)+1.d0)*dzetat(2)
	   dbidro=dbidro+bid1*cte7/(zeta(2)+1.d0)*dzetaro(2)
	   DO l=1,nchim
	    dbidx(l)=dbidx(l)+bid1*cte7/(zeta(2)+1.d0)*dzetax(2,l)
	   ENDDO
	   df=cte5*c1**2/k**5*bid		!eq. 25
	   dft= df*(-5.d0*dkt /k+dbidt /bid)
	   dfro=df*(-5.d0*dkro/k+dbidro/bid+2.d0*dc1ro/c1)
	   DO l=1,nchim
	    dfx(l)=df*(-5.d0*dkx(l)/k+dbidx(l)/bid+2.d0*dc1x(l)/c1)
	   ENDDO

	   ecran(i)=ABS(df) ; decrt(i)=SIGN(dft,df)
	   decro(i)=SIGN(dfro,df)
	   DO l=1,nchim
	    decrx(i,l)=SIGN(dfx(l),df)
	   ENDDO

	   IF(.FALSE.)THEN	!suppression du terme d'ordre supérieur
c	   IF(.TRUE.)THEN
	    e0=t**(2.d0/3.d0)*mz1z2(i)	!eq. 3' i.e. Lang 4-44 E0 Gamov
	    bid1=(zeta(1)+zeta(2)+1.d0)**cte7-cte7*(zeta(1)+zeta(2))-1.d0!eq.32'
	    beta=cte8*df/lambda(i)/k*bid1/bid  !eq. 32'
	    eps=beta*(lambda(i)*k/(e0+df))**3	!eq. 38' 39 + remarque eq. 44
	    ecran(i)=ecran(i)-cte10*e0*eps*(1.d0+cte9*eps)	!eq. 58
	   ENDIF
	   ecran(i)=ecran(i)/kbol/t			!eq. 58
	   decrt(i)=decrt(i)/kbol/t-ecran(i)/t
	   decro(i)=decro(i)/kbol/t
	   DO l=1,nchim
	    decrx(i,l)=decrx(i,l)/kbol/t
	   ENDDO

c exp(f)
c          WRITE(*,7)ecran(i)
7	   FORMAT('mitler',es10.3)
	   ecran(i)=EXP(ecran(i))			!eq. 58
	   decrt(i)=ecran(i)*decrt(i) ; decro(i)=ecran(i)*decro(i)
	   DO l=1,nchim
	    decrx(i,l)=ecran(i)*decrx(i,l)
	   ENDDO

c	   WRITE(*,2000)df,cte10*e0*eps*(1.d0+cte9*eps)
c	   PRINT*,'zeta(1),zeta(2),bid,df,e0,bid1,beta,eps'
c	   WRITE(*,2000)zeta(1),zeta(2),bid,df,e0,bid1,beta,eps
c	   PRINT*,'t,ecran(i),izz(i,1),izz(i,2)',izz(i,1),izz(i,2)
c	   WRITE(*,2000)t,ecran(i),EXP(.188*SQRT(ro*(abon(1)+3.)/2./
c	1  (t*1.d-6)**3))
c	   WRITE(*,2000)EXP(.188*SQRT(ro*(abon(1)+3.)/2./(t*1.d-6)**3))
c	   PAUSE
	  ELSE
	   ecran(i)=1.d0 !si z1 ou z2=0 par exemple pour Be7(e-,nu g)Li7
	   decrt(i)=0.d0 ; decro(i)=0.d0 ; decrx(i,:)=0.d0
	  ENDIF
	 ENDDO	B1

c écran Salpeter Clayton 4-221, Kippenhahn 18-57
	ELSE Im
	 mue=DOT_PRODUCT(zi,abon)	!1/mue
	 dmuex=zi
	 zstar=0.d0
	 DO j=1,nchim
	  zstar=zstar+zi(j)*(zi(j)+1.d0)*abon(j)
	  dzstarx(j)=zi(j)*(zi(j)+1.d0)
	 ENDDO

	 B2: DO i=1,nreac
c	  IF(izz(i,1) > 0)THEN
	  I12: IF(z1z2(i) > 0)THEN
	   ecran(i)=1.88d8*z1z2(i)*SQRT(zstar*ro/t**3)	!Clayton 4-221, Kippenhahn 18-57
	   decrt(i)=-3.d0/2.d0*ecran(i)/t
	   decro(i)=ecran(i)/ro/2.d0
	   DO j=1,nchim
	    decrx(i,j)=ecran(i)/2.d0/zstar*dzstarx(j)
	   ENDDO
c	   WRITE(*,8)ecran(i)
8	   FORMAT('faible',es10.3)

c écran fort Kippenhahn 18-58, Clayton 4-226 !!mue=1/mue!!
	   IF(ecran(i) > 1.d0)THEN
c	   IF(.FALSE.)THEN
c	   IF(.TRUE.)THEN
            bid=ecran(i)
	    bid1=2.05d5*zz53(i)*(ro*mue)**(1.d0/3.d0)/t
c	    WRITE(*,9)bid1 ; PAUSE'bid'
9	    FORMAT('fort',es10.3)
            IF(bid1 > bid)THEN
	     ecran(i)=bid1
	     decrt(i)=-ecran(i)/t
	     decro(i)=ecran(i)/ro/3.d0
	     DO j=1,nchim
	      decrx(i,j)=ecran(i)/3.d0/mue*dmuex(i)
	     ENDDO
	     IF(fort)THEN
              fort=.FALSE.
              WRITE(*,6)bid,bid1,TRIM(nom_reac(i)) ; WRITE(2,6)bid,bid1,TRIM(nom_reac(i))
6	      FORMAT("Au moins une fois, écrantage fort:",es10.3,
	1     ' --->', es10.3,' réaction:',a)
	     ENDIF
	    ENDIF
	   ENDIF

c EXP(ecran), ecran >= 1
	   ecran(i)=EXP(ecran(i))
	   decrt(i)=ecran(i)*decrt(i)
	   decro(i)=ecran(i)*decro(i)
	   DO j=1,nchim
	    decrx(i,j)=ecran(i)*decrx(i,j)
	   ENDDO

c photodissociations ou capture électroniques
	  ELSE I12
	   ecran(i)=1.d0 !si z1 ou z2=0 par exemple pour Be7(e-,nu g)Li7
	   decrt(i)=0.d0 ; decro(i)=0.d0 ; decrx(i,:)=0.d0
	  ENDIF I12
	 ENDDO B2
c	 PRINT*,'ecran(6)',ecran(6)
	ENDIF	Im	!mitler

c interpolation en ln t
c	PRINT*,nreac,n_temp,m_temp,knot_temp,l
c	WRITE(*,2000)temp(1),ttemp(1),LOG(t)
c	WRITE(*,2000)taux_reac(1,1)
	CALL bsp1dn(nreac,taux_reac,temp,ttemp,n_temp,m_temp,knot_temp,
	1 .TRUE.,LOG(t),l,r,dr)
c	PRINT*,'dr(6)',dr(6)

c limitation de la réaction 6 : Be7(e-,nu g)Li7 (Caughlan & Fowler)
c toujours utilisée dans le PP
	IF(t < 1.d6)THEN
	 bid=LOG(1.51d-7/2.d0/ro/(1.d0+abon(1)*nucleo(1)))
	 r(6)=MIN(r(6),bid)
	ENDIF

c taux des réactions Rij (écran > = 1)
	DO i=1,nreac
	 IF(i == i3al)THEN
	  r(i)=EXP(r(i))*ro**2*ecran(i)		!réaction 3 alpha
	  drt(i)=r(i)*(dr(i)/t+decrt(i)/ecran(i))
	  dro(i)=r(i)*(2.d0/ro+decro(i)/ecran(i))
	  DO j=1,nchim
	   drx(i,j)=r(i)/ecran(i)*decrx(i,j)
	  ENDDO
	 ELSE		!les autres réactions
	  r(i)=EXP(r(i))*ro*ecran(i)
	  drt(i)=r(i)*(dr(i)/t+decrt(i)/ecran(i))
	  dro(i)=r(i)*(1.d0/ro+decro(i)/ecran(i))
	  DO j=1,nchim
	   drx(i,j)=r(i)/ecran(i)*decrx(i,j)
	  ENDDO
	 ENDIF

c énergie libérée par les réactions Q, q0, du module mod_nuc,  est formé dans tabul_nuc
	 q(i)=r(i)*q0(i) ; dqt(i)=drt(i)*q0(i) ; dqo(i)=dro(i)*q0(i)
	 DO j=1,nchim
	  dqx(i,j)=drx(i,j)*q0(i)
	 ENDDO
	ENDDO	!i

c sortie OPTIONAL pour tests
	IF(PRESENT(ecranw))ecranw=ecran

c pour test de l'écrantage
c	DO i=1,nreac
c	 PRINT*,'reac:',i,', ',nom_reac(i)
c	 WRITE(*,2000)LOG(ecran(i)),ecran(i)
c	ENDDO

	RETURN

	END SUBROUTINE rq_reac