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

	SUBROUTINE test_poisson
	
c routine private du module mod_cesam	
	
c test de l'intégration de l'équation de Poisson
c on forme la relation 3.23 de Zahn 1992

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

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

	USE mod_donnees, ONLY : en_masse, g, msol, m_ch, m_qs, m_rot, nchim,
	1 ne, nrot, ord_qs, pi, rsol
	USE mod_etat, ONLY: etat
	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), SAVE, ALLOCATABLE, DIMENSION(:) :: fd, fdq, fq, fs
	REAL (kind=dp), DIMENSION(nchim) :: dxchim, xchim		
	REAL (kind=dp), SAVE, DIMENSION(4) :: cte_0
	REAL (kind=dp), DIMENSION(4) :: f
	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,
	2 drop, drot, drox,dup, dut,  dux, gradad, gamma1,
	3 grav, mr, nu, 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.

c constantes	 
	 cte3_0=g*msol/rsol**2	
	 
	 cts_0(1)=pi*rsol**3/msol
	 cts_0(2)=2.d0*pi*rsol**3/msol 
	 
	 cte_0(1)=4.d0/3.d0*rsol	 
	 cte_0(2)=2.d0/rsol
	 cte_0(3)=8.d0*pi*rsol**2/3.d0/msol
	 cte_0(4)=4.d0/3.d0*rsol**3/g/msol
	 	 	 	 	 
	 ALLOCATE(fd(nrot),fdq(ne),fq(ne),fs(nrot))	  
	ENDIF

	DO i=2,n_qs
	 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),l,fq,fdq)	
	 nu=fq(5)
	 p=EXP(fq(1)) ; t=EXP(fq(2))
	 IF(en_masse)THEN
	  r=SQRT(ABS(fq(3))) ; nu12=SQRT(fq(5))
	 ELSE
	  r=ABS(fq(3)) ; nu12=fq(5)**(1/3.d0)
	 ENDIF
	 nu=nu12**2 ; nu32=nu12**3
	 
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 coeff_rota4'	 
	
c composition chimique par gramme, équation d'état
	 CALL chim_gram(xchim,dxchim)
	 CALL etat(p,t,xchim,.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 la rotation
	 CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,nu,l,fs,fd)

	 cts(1)=r**3*ro/nu32
	 cts(2)=r**3*ro/nu32	
	 cts=cts*cts_0

	 f(1)=r/grav*(1.d0-cts(1))*fs(1)**2
	 f(2)=1.d0/r/grav*(cts(2)-1.d0)*fs(5)
	 f(3)=r**2*ro/grav/nu12*fd(5)
	 f(4)=r**3/nu32*fs(1)**2
	 f=f*cte_0
	 
c écritures 
	 WRITE(*,2000)nu, f(1:3)
	 WRITE(*,2000)f(1)+f(2)-f(3),f(4) ; PAUSE'test'	
	ENDDO

	RETURN	
	
	END SUBROUTINE test_poisson