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

	SUBROUTINE poisson_initial
	
c routine private du module mod_cesam

c initialisation de la représentation spline de la variable Phi
c résolution de l'équation de Poisson du formalisme Matis & Zahn 2004
c par éléments finis

c les coefficients de Omega de la spline d'interpolation de la
c rotation, ie; rota(1,:), doivent être connus, ainsi que les splines
c d'interpolation de la composition chimique et des variables quasi-statiques

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

c------------------------------------------------------------------------------
	
	USE mod_donnees, ONLY : m_rot, nom_rot, nrot
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, bvald, gauss_band, intgauss,
	1 linf, no_croiss
	USE mod_variables, ONLY : dim_rot, knotr, mrot, mrott, n_rot, rota
	
	IMPLICIT NONE
			
	REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: a, b
	REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fd, fs	
	REAL(kind=dp), DIMENSION(0:1) :: ad, as
	REAL(kind=dp), DIMENSION(0:1,m_rot) :: d	
	REAL(kind=dp), DIMENSION(m_rot) :: qg, wg
	REAL(kind=dp) :: bs	
		
	INTEGER, ALLOCATABLE, DIMENSION(:) :: indpc
	INTEGER, SAVE :: bloc
	INTEGER :: col, colonne, i, id, ig, j, k, l=1, ligne, nl, rang, spi

	LOGICAL, SAVE :: init=.TRUE.
	LOGICAL :: inversible
			
c-----------------------------------------------------------------------------

2000	FORMAT(8es10.3)

c initialisations,  bloc: longueur du bloc des coeff. non idt. nuls	
	IF(init)THEN
	 init=.FALSE.
	 bloc=2*m_rot-1		!une seule variable : Phi
	ENDIF
	
c nl = rang : nombre de lignes du système linéaire	
	rang=dim_rot ; nl=rang
	
c allocations				
	ALLOCATE(a(nl,bloc),b(1,nl),indpc(nl))

c mises à 0	
	a=0.d0 ; b=0.d0 ; indpc=0

c indices de première colonne pour les produits scalaires dans gauss_band
	DO i=1,dim_rot      !pour chaque spline de la base des mrot
	 indpc(i)=MAX(1,i-m_rot+1)	!indice il n'y a qu'une variable
	ENDDO

c---------------début de la construction du sytème linéaire------------	
	    
c les produits scalaires	  
	DO k=1,n_rot-1	  
	  
c spi : indice-1 de la première spline non id.nulle  
c comme mrott(l) <= qg(ig) < mrott(l+1), l est le même pour toutes
c les abscisses de Gauss dans l'intervalle [mrott(k),mrott(k+1)]	  
	 CALL linf(mrot(k),mrott,knotr,l) ; spi=l-m_rot
	 IF(no_croiss)PRINT*,'poisson_initial, Pb. en 0'
	   
c poids, wg et abscisses, qg pour intégration Gauss d'ordre m_rot
	 CALL intgauss(mrot(k),mrot(k+1),qg,wg,m_rot)	  
	  
c pour chaque abscisse de GAUSS
	 DO ig=1,m_rot
	   
c formation des coefficients des équations
	  CALL eq_diff_poisson(0,qg(ig),ad,as,bs)	    
	   		   
c les B-splines en qg(ig), dérivées 0 à 1
	  CALL bvald(qg(ig),mrott,m_rot,l,1,d)

c contribution au système linéaire 	   
	  DO i=1,m_rot	!pour chaque spline d'indice i
	   ligne=spi+i
	   b(1,ligne)=b(1,ligne)+wg(ig)*d(0,i)*bs
		
c la matrice compressée est le jacobien 'diagonal' ie. sans les
c éléments 'non diagonaux' identiquement nuls
	   DO j=1,m_rot		!pour chaque spline j       
	    colonne=spi+j  
	    col=colonne-indpc(ligne)+1   !colonne matrice compressée
c	    PRINT*,'ligne,colonne,col',ligne,colonne,col;PAUSE'1'
	    DO id=0,1  !pour 0: NiNj, 1: NiN'j, 3: NiN"j etc...
	     a(ligne,col)=a(ligne,col)+wg(ig)*
	1    (as(id)*d(0,i)+ad(id)*d(1,i))*d(id,j)
	    ENDDO	!id
	   ENDDO	!j
	  ENDDO		!i
	 ENDDO 		!ig
	ENDDO 		!k

c quantités intégrées sur la limite externe mrot(n_rot), toujours convectif
c avec le signe -
c variables et dérivées

c formation des coefficients des équations
	CALL eq_diff_poisson(1,mrot(n_rot),ad,as,bs)	    

c les B-splines en mrot(n_rot), dérivées 0 à 1
	CALL bvald(mrot(n_rot),mrott,m_rot,l,1,d)	 
 
c contribution au système linéaire pour Phi 
	 ligne=l
	 DO j=1,m_rot		!pour chaque spline j	   
	  colonne=spi+j    
	  col=colonne-indpc(ligne)+1   !colonne matrice compressée
	  a(ligne,col)=a(ligne,col)-as(1)*d(1,j)!dPhi seul est concerné
	 ENDDO		!j
 
c----------------fin de la construction du système linéaire-------------- 

c	DO i=1,nl
c	 WRITE(*,2000)a(i,:),b(1,i)
c	ENDDO
c	PAUSE'sys lin'
c	WRITE(*,2000)b(1,:) ; PAUSE'2-nd membres'
		
c résolution du système linéaire
	 CALL gauss_band(a,b,indpc,nl,rang,bloc,1,inversible)
	 IF(.NOT.inversible)THEN
	  PRINT*,'ARRRET, matrice singulière dans poisson_initial' ; STOP
	 ENDIF
c	 WRITE(*,2000)b(1,:) ; PAUSE'solution poisson initial'
	 
c on place dans rota(5,:) les coefficients de spline pour Phi 
	 rota(5,:)=b(1,:)
	 rota(5,1)=0.d0
	 
c vérification
c	IF(.TRUE.)THEN
	IF(.FALSE.)THEN
	 ALLOCATE(fd(nrot),fs(nrot))	
	 DO i=1,n_rot	
	  CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,mrot(i),l,
	1 fs,fd)
	  WRITE(*,2000)mrot(i),fs(5),fd(5)
	 ENDDO
	 DEALLOCATE(fd,fs)
	 PAUSE'solution linéaire'
	ENDIF	 

c nettoyage		
	DEALLOCATE(a,b,indpc)

	IF(.FALSE.)CALL test_poisson
c	IF(.TRUE.)CALL test_poisson
	
	RETURN		
		
	END SUBROUTINE poisson_initial