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

	SUBROUTINE initialise_rota
	
c	routine private du module mod_cesam
c	initialisation de la rotation rota(Omega,U,theta,Lambda,Psi)

c	pour un modèle de ZAMS ou de PMS initial ou encore pour un modèle
c	repris en ASCII les masses de rota et de bp coincident

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_rot, ne,
	1 nrot, ord_qs, pi, rsol, w_rot
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, no_croiss
	USE mod_variables, ONLY : bp, dim_ch, dim_rot, knot, knotc, knotr,
	1 mc, mct, mrot, mrott, n_ch, n_qs, n_rot, q, qt, rota

	IMPLICIT NONE
	
	REAL (kind=dp), DIMENSION(ne) :: dfqs, fqs
	REAL (kind=dp), DIMENSION(nrot) :: dfrot, frot		
	REAL (kind=dp), SAVE :: cte1_0	
	REAL (kind=dp) :: r, rstar

	INTEGER :: i, l=1
				
c--------------------------------------------------------------------------------

2000	FORMAT(8es10.3)

	cte1_0=16.d0*pi*rsol**6/9.d0/g/msol**2
	
c	initialisation de l'interpolation de la rotation
c	pour plus de clarté, certaines définitions sont volontairement
c	naïves

	n_rot=n_ch ; knotr=n_rot+m_rot ; dim_rot=knotr-m_rot	
	ALLOCATE(mrot(n_rot),mrott(knotr),rota(nrot,dim_rot))
	mrot=mc ; rota=0.d0

c	détermination de rstar
	CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(n_qs),l,fqs,dfqs)
	IF(no_croiss)PRINT*,'Pb. en 1 dans initialise_rota'
	IF(en_masse)THEN
	 rstar=SQRT(ABS(fqs(3)))
	ELSE
	 rstar=ABS(fqs(3))
	ENDIF
	
c	initialisation de rota(1,i)=Omega et U=rota(1,i)	
	DO i=1,n_qs
	 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),l,fqs,dfqs)
	 IF(no_croiss)PRINT*,'Pb. en 2 dans initialise_rota'
	 IF(en_masse)THEN
	  r=SQRT(ABS(fqs(3)))
	 ELSE
	  r=ABS(fqs(3))
	 ENDIF
	 rota(1,i)=w_initial(r,rstar) ; rota(2,i)=u_initial(r,rstar) 
	ENDDO
	
c	pour le calcul de Theta initialisation provisoire de rota

	CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.FALSE.,mrot(1),
	1 l,frot,dfrot)
	IF(no_croiss)PRINT*,'Pb. en 3 dans initialise_rota'
			
c	initialisation de Theta et Psi=-Theta (n_ch et n_qs coïncident)
	DO i=1,n_rot
	 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),l,fqs,dfqs)
	 IF(no_croiss)PRINT*,'Pb. en  4 dans initialise_rota'
	 IF(en_masse)THEN
	  r=SQRT(ABS(fqs(3)))
	 ELSE
	  r=ABS(fqs(3))
	 ENDIF
	 CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,mrot(i),
	1  l,frot,dfrot)
	 IF(no_croiss)PRINT*,'Pb. en 5 dans initialise_rota'
	 IF(mrot(i) <= 0.d0)THEN
	  rota(3,i)=0.d0
	 ELSE 
	  rota(3,i)=cte1_0*r**6*EXP(fqs(7))/mrot(i)**2*frot(1)*dfrot(1)
	 ENDIF
	 rota(5,i)=-rota(3,i)
	ENDDO
c	DO i=1,n_rot
c	 WRITE(*,2000)mrot(i),rota(:,i)
c	ENDDO
c	PAUSE'rota_ini'
	
	RETURN	
	
	END SUBROUTINE initialise_rota