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

	SUBROUTINE eq_ini_rota4(fait,nu,as,bs)	
	
c routine PRIVATE du module mod_evol

c formation des équations à résoudre pour la résolution de de l'équation de
c Poisson. Omega et dOmega sont fixés, le système est alors linéaire

c Auteur: P.Morel, Département Cassiopée, O.C.A., 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 : g, en_m23, msol, m_ch, nrot, nvth, ord_rot, pi,
	1 rsol, w_rot
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, no_croiss
	USE mod_variables, ONLY : knotc, knotr, mc, mct, mrot, mrott,
	1 n_ch, n_rot, rota, vth
				
	IMPLICIT NONE

	REAL (kind=dp), INTENT(in) :: nu
	INTEGER, INTENT(in) :: fait
	
	REAL (kind=dp), INTENT(out), DIMENSION(:,:,0:) :: as
	REAL (kind=dp), INTENT(out), DIMENSION(:) :: bs
	REAL (kind=dp), DIMENSION(nvth) :: dfvth, fvth
	REAL (kind=dp), SAVE, DIMENSION(3) :: cte_0
	REAL (kind=dp), DIMENSION(3) :: fr
	REAL (kind=dp), SAVE, DIMENSION(2) ::cts_0
	REAL (kind=dp), DIMENSION(2) :: cts
	
	REAL (kind=dp), SAVE ::	cte1_0	
	REAL (kind=dp) :: dlnro, grav, nu12, nu32, Omega, r, ro

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

2000	FORMAT(8es10.3)

c initialisations
	IF(init)THEN
	 init=.FALSE.
	 
c pour la gravité
	 cte1_0=g*msol/rsol**2
	 	 
c pour les fr (F)
	 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
	 	 
c pour les F* 	 
	 cts_0(1)=9.d0*msol/4.d0/pi/rsol**3
	 cts_0(2)=4.d0*pi*g*rsol
	 
	ENDIF

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

c les coefficients au temps t+dt	
	
c nu12 etc...
	nu12=SQRT(nu) ; nu32=nu12**3

c utilisation de la tabulation vth	
c fvth(1)=lnP, fvth(2)=lnT, fvth(5)=ln ro, 
c fvth(6)=cp, fvth(7)=delta, fvth(8)=gamma1, fvth(9)=ln µ,
c fvth(10)=ln kap, fvth(11)=degene
c en m^2/3 : fvth(3)=r**2, fvth(4)=l**2/3
c en m^1/3 : fvth(3)=r, fvth(4)=l	
	CALL bsp1dn(nvth,vth,mc,mct,n_ch,m_ch,knotc,.TRUE.,
	1 MAX(mc(1),MIN(nu,mc(n_ch))),l,fvth,dfvth)
	IF(no_croiss)PRINT*,'Pb. at 2 in eq_ini_rota4'	
	ro=EXP(fvth(5)) ; dlnro=dfvth(5) ; r=SQRT(ABS(fvth(3)))
	grav=cte1_0*nu32/r**2
	
	IF(.NOT.en_m23)THEN
	 PRINT*,'en m13 modifier fvth(3)=r, fvth(4)=l' ; STOP
	ENDIF
				
c les cts*
	cts(1)=nu12/r**3/ro
	cts(2)=r*ro/grav*dlnro
	cts=cts*cts_0

c les coefficients fr(1)=F25, fr(2)=F27, fr(3)=F29
	fr(1)=r**3*ro/nu12
	fr(2)=cts(1)+cts(2)
	fr(3)=r**3*ro/grav*dlnro
	fr=fr*cte_0

	Omega=ABS(w_rot)
		
c mises à 0
	as=0.d0 ; bs=0.d0	

	SELECT CASE(fait)
	CASE(0)

c Omega	
	 bs(1)=Omega ; as(1,1,0)=1.d0
	 
c U=Omega/100 (Th. Corbard)	 
c	 bs(2)=Omega/100.d0 ; as(2,2,0)=1.d0
	 
c U, Psi, Lambda, Tau, Upsilon = 0.d0		
	 bs(2:6)=0.d0
	 DO i=2,6
	  as(i,i,0)=1.d0
	 ENDDO

c dPi - F27 Phi = F29 Omega^2 (+ F28 Omega*dOmega = 0)
	 bs(7)=fr(3)*Omega**2 ; as(7,7,0)=-fr(2) ; as(7,8,1)=1.d0
	 
c Pi - Phi - F25 dPhi=0
	 bs(8)=0.d0 ; as(8,7,0)=-1.d0 ; as(8,7,1)=-fr(1) ; as(8,8,0)=1.d0

c condition limite au centre
	CASE(1)
	
c Omega	
	 bs(1)=Omega ; as(1,1,0)=1.d0	 

c U, Psi =0
	 bs(2:3)=0.d0 ; as(2,2,0)=1.d0 ; as(3,3,0)=1.d0
	 
c Phi = 0	
	 bs(4)=0.d0 ; as(4,7,0)=1.d0
	 	
c condition limite en surface
	CASE(2)
	
c Lambda, Tau, Upsilon =0.d0
	 bs(1:3)=0 ; as(1,4,0)=1.d0 ; as(2,5,0)=1.d0 ; as(3,6,0)=1.d0

c Phi + F25/3 dPhi=0		
	 bs(4)=0.d0 ; as(4,7,0)=1.d0 ; as(4,7,1)=fr(1)/3.d0	 
	 
	END SELECT

	RETURN

	END SUBROUTINE eq_ini_rota4