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

	SUBROUTINE eq_diff_poisson(fait,nu,ad,as,bs)	
	
c routine private du module mod_evol

c formation des équations à résoudre dans poisson_initial
c pour la rotation par la méthode des éléments finis Galerkin

c Auteur: P.Morel, Département Cassiopée, O.C.A.
c CESAM2k

c entrées
c	fait=0 : point courant, fait=1 : extérieur
c	nu: point de calcul en m**2/3
c	dt : pas temporel

c sorties	
c	as, ad : coefficients de la spline/dérivée
c	bs : second membre

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

	USE mod_donnees, ONLY : en_masse, g, msol, mu_saha, m_ch, m_qs,
	1 m_rot, nchim, ne, nrot, phi, pi, rsol, zi
	USE mod_etat, ONLY: etat, saha
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, no_croiss
	USE mod_variables, ONLY : bp, chim, chim_gram,
	1 inter, knot, knotc, knotr, mc, mct, mrot, mrott, m23, n_ch, n_qs, q,
	2 qt, n_rot, rota, r2
				
	IMPLICIT NONE

	REAL (kind=dp), INTENT(in) :: nu
	INTEGER, INTENT(in) :: fait
	
	REAL (kind=dp), INTENT(out), DIMENSION(0:) :: ad, as
	REAL (kind=dp), INTENT(out) :: bs
	
	
	REAL (kind=dp), DIMENSION(0:NINT(MAXVAL(zi)),nchim) :: ioni
	REAL (kind=dp), SAVE, ALLOCATABLE, DIMENSION(:) :: fd, fs, z_bar	
	REAL (kind=dp), DIMENSION(nchim) :: dxchim, dxchimg, xchim, xchimg	
	REAL (kind=dp), DIMENSION(ne) :: dfqs, fqs
	REAL (kind=dp), SAVE, DIMENSION(4) :: cte_0
	REAL (kind=dp), DIMENSION(4) :: frl
	REAL (kind=dp), SAVE, DIMENSION(2) ::cts_0
	REAL (kind=dp), DIMENSION(2) :: cts
	
	REAL (kind=dp), SAVE ::	cte3_0	
	REAL (kind=dp) :: alfa, beta, cp, dcpp, dcpt, dcpx, delta, deltap,
	1 deltat,deltax, dgradadp, dgradadt, dgradadx, dlnp,
	2 dlnmu, dlnro_nu, drop, drot, drox, dup, dut,  dux, eta, grad,
	3 gradad, gamma1, grav, mr, mu, nel, p, nu12, nu32, r, ro, t, u

	INTEGER, SAVE :: l=1
	INTEGER :: i
		
	LOGICAL, SAVE :: init=.TRUE.
		
c----------------------------------------------------------------------

2000	FORMAT(8es10.3)

c initialisations
	IF(init)THEN
	 init=.FALSE.
	 
	 cte3_0=g*msol/rsol**2	 
	 
	 cte_0(1)=8.d0*pi*rsol**3/3.d0/msol	 
	 cte_0(2)=1.d0
	 cte_0(3)=4.d0*pi*g*rsol**3/3.d0
	 cte_0(4)=8.d0*pi*g*rsol**3/3.d0	 
	 	 
	 cts_0(1)=9.d0*msol/4.d0/pi/rsol**3
	 cts_0(2)=4.d0*pi*g*rsol
	 
	 ALLOCATE(fd(nrot),fs(nrot))
	 
c r2 et m23 ne sont peut être pas encore alloués	
	 IF(.NOT.ALLOCATED(r2))THEN
	  ALLOCATE(r2(n_qs),m23(n_qs))
	  DO i=1,n_qs	   
	   r2(i)=bp(3,m_qs*(i-1)+1) ; m23(i)=bp(5,m_qs*(i-1)+1)
	  ENDDO
	 ENDIF
	 
c charges moyennes dans l'hypothèse totalement ionisé
	 ALLOCATE(z_bar(nchim))	 
	 IF(.NOT.mu_saha)z_bar=zi	 
	 	  
	ENDIF

c~~~~~~~~~~~~~~~~fin des initialisations~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c les coefficients au temps t+dt	
	
c nu12 etc...
	nu12=SQRT(nu) ; nu32=nu12**3
			
c Omega et dOmega/dnu		
	 CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,nu,l,fs,fd)
		
c ATTENTION, dérivées / m^2/3 ou m effectuées dans inter		 
	CALL inter('m23',bp,q,qt,n_qs,knot,nu,fqs,dfqs,r2,m23)
	IF(no_croiss)PRINT*,'Pb en 0 dans eq_diff_poisson'
	p=EXP(fqs(1)) ; t=EXP(fqs(2))
	IF(en_masse)THEN
	 r=SQRT(ABS(fqs(3))) ; dlnp=dfqs(1)
	ELSE
	 r=ABS(fqs(3)) ; dlnp=3.d0/2.d0*dfqs(1)*nu12
	ENDIF
	
c grad = d ln T / d ln P
	grad=dfqs(2)/dfqs(1)	
	
c gravité
	grav=cte3_0*nu32/r**2
	
c composition chimique
	mr=MAX(mc(1),MIN(nu,mc(n_ch)))
	CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mr,l,xchim,dxchim)
	IF(no_croiss)PRINT*,'Pb en 1 dans eq_diff_poisson'	 
	
c composition chimique par gramme, équation d'état
	xchimg=ABS(xchim) ; dxchimg=dxchim
	CALL chim_gram(xchimg,dxchimg)
	CALL etat(p,t,xchimg,.FALSE.,ro,drop,drot,drox,u,dup,dut,dux,
	1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)
				
c poids moléculaire moyen	
c appel à saha pour les charges moyennes
	IF(mu_saha)CALL saha(xchim,t,ro,ioni,z_bar,nel,eta)
	mu=1.d0/DOT_PRODUCT((1.d0+z_bar),xchim)
c	mu=16.d0/(23.d0*xchimg(1)+3.d0*xchimg(ihe4)+9.d0)	
	dlnmu=-mu*DOT_PRODUCT((1.d0+z_bar),dxchim)
c	dlnmu=-mu/16.d0*(23.d0*dxchimg(1)+3.d0*dxchimg(ihe4))

c d ln ro / d nu
	dlnro_nu=(alfa-delta*grad)*dlnp+phi*dlnmu
		
c les coefficients frl(1)=F25, frl(2)=F27, frl(3)=F29, frl(4)=F28
	cts(1)=nu12/r**3/ro
	cts(2)=r*ro/grav*dlnro_nu
	cts=cts*cts_0

	frl(1)=r**3*ro/nu12
	frl(2)=cts(1)+cts(2)
	frl(3)=r**3*ro/grav*dlnro_nu
	frl(4)=r**3*ro/grav
	frl=frl*cte_0
		
c les équations

c mises à 0
	ad=0.d0 ; as=0.d0 ; bs=0.d0	

	SELECT CASE(fait)
	CASE(0)
c Phi
	 bs=-frl(3)*fs(1)**2-frl(4)*fs(1)*fd(1)
	 as(0)=frl(2)
	 ad(0)=1.d0
	 ad(1)=frl(1)
	 
c	IF(nu < mrot(3))WRITE(*,2000)nu,frl,fs(1),fd(1)

c condition limite à l'extérieur
c pour le système linéaire, le second membre est inchangé, on retire au
c premier membre la quantité intégrée	 
	CASE(1)
c	 as(1)=frl(1)*2.d0/3.d0
	 as(1)=frl(1)/2.d0
	 	
c parties intégrées, au signe près sur la face radiative des limites ZR/ZC
	CASE(2)
	 as(0)=1.d0 ; as(1)=frl(1)	 
	 
	END SELECT
	
c	WRITE(*,2000)nu,nu12,r,ro,grav,dlnro_nu
c	WRITE(*,2000)cts,cts_0 ; WRITE(*,2000)frl,cte_0
c	WRITE(*,2000)bs,as,ad ; PAUSE'eq_diff_poisson'

	RETURN

	END SUBROUTINE eq_diff_poisson