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

	SUBROUTINE initialise_poisson
	
c routine private du module mod_evol

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, xli	
		
	INTEGER, ALLOCATABLE, DIMENSION(:) :: indpc
	INTEGER, SAVE :: bloc
	INTEGER :: col, colonne, i, id, ig, ik, 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~~~~~~~~~~~~~~au centre~~~~~~~~~~~~~~~~~~~~~~~~~

c	ligne=1 ; b(1,ligne)=0.d0 ; a(ligne,:)=0.d0

c Phi=0
c	a(ligne,1)=1.d0
		
c dPhi / dnu = 0

c localisation
	CALL linf(mrot(1),mrott,knotr,l)
	
c les B-splines		
	CALL bvald(mrot(1),mrott,m_rot,l,1,d)

c avec le système linéaire, le second membre est nul
c au premier membre les dérivées des B-splines en mrot(1)
	DO i=1,m_rot
	 a(ligne,i)=d(1,i)
	ENDDO	

c~~~~~~~~~~~limite externe~~~~~~~~~~~~~~~~~~~

c quantité intégrée sur la limite externe mrot(n_rot), toujours convectif
c avec le signe -

c localisation
	CALL linf(mrot(n_rot),mrott,knotr,l)

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 pour le système linéaire, le second membre est inchangé, on retire au
c premier membre la quantité intégrée	 
	 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~~~~~~~~~~~~~

!!!!!!!!!!!!!!!
c Il semble que la moins mauvaise solution soit obtenue sans se soucier de
c continuité ou de quantités intégrées
!!!!!!!!!!!!!!!

	If(.TRUE.)Then	!poisson continu
c	If(.FALSE.)Then	
	
	 DO ik=1,ndis
	  CALL linf(mrot(idis(ik)),mrott,knotr,l) ; spi=l-m_rot

c convectif à droite de la limite, utilisation de la spline spi+1
	  IF(lmix(mrot(idis(ik))))THEN	  
	   ligne=spi+1
	   col=m_rot-1 ; a(ligne,:)=0.d0
	   a(ligne,col)=1.d0 ; a(ligne,col+1)=-1.d0
	   
c convectif à gauche de la limite, utilisation de la spline spi
	  ELSE  
	     ligne=spi    
	     col=m_rot ; a(ligne,:)=0.d0
	     a(ligne,col)=1.d0 ; a(ligne,col+1)=-1.d0   
	  ENDIF		!lmix
	 ENDDO		!ik
	 	
	Else		!poisson discontinu
	 
c quantités intégrées sur les faces radiatives des limites ZR/ZC
	 DO ik=1,ndis

c radiatif à gauche de la limite affecte la spline spi avec le signe +
	  IF(lmix(mrot(idis(ik))))THEN
	  
c on se déplace un peu à gauche de la limite	  
	   xli=mrot(idis(ik))*un_eps+mrot(idis(ik)-1)*eps

c localisation
	   CALL linf(xli,mrott,knotr,l) ; spi=l-m_rot

c formation des coefficients des équations
	   CALL eq_diff_poisson(2,xli,ad,as,bs)	    

c les B-splines en mrot(idis(ik)), dérivées 0 à 1
c avec le l de xli, mrot(idis(ik)) pour avoir des valeurs exactes
	   CALL bvald(mrot(idis(ik)),mrott,m_rot,l,1,d)	 
	    
c contribution au système linéaire
	   ligne=nrot*(l-1)+1
	   DO j=1,m_rot		!pour chaque spline j	   
	    colonne=nrot*(spi+j-1)+1    
	    col=colonne-indpc(ligne)+1   !colonne matrice compressée
	    DO id=0,1
	     a(ligne,col)=a(ligne,col)+as(id)*d(id,j)
	    ENDDO	!id
	   ENDDO	!j

c radiatif à droite de la limite, affecte la spline spi+1 avec le signe -
	  ELSE

c localisation
	   CALL linf(mrot(idis(ik)),mrott,knotr,l) ; spi=l-m_rot	  
	  
c formation des coefficients des équations
	   CALL eq_diff_poisson(2,mrot(idis(ik)),ad,as,bs)	    

c les B-splines en mrot(idis(ik)), dérivées 0 à 1
	   CALL bvald(mrot(idis(ik)),mrott,m_rot,l,1,d)	 

c contribution au système linéaire, Omega, Uf et Phi sont concernés
	   ligne=nrot*spi+1
	   DO j=1,m_rot		!pour chaque spline j	   
	    colonne=nrot*(spi+j-1)+1    
	    col=colonne-indpc(ligne)+1   !colonne matrice compressée
	    DO id=0,1
	     a(ligne,col)=a(ligne,col)-as(id)*d(id,j)
	    ENDDO	!id
	   ENDDO	!j
	  ENDIF		!lmix
	 ENDDO		!ik
	 
	 Endif	!suppression de la continuité
 
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)) ; ik=1	
	 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)
	
	RETURN		
		
	END SUBROUTINE initialise_poisson