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

	SUBROUTINE resout_rota4(dt,ok)

c routine public du module mod_evol
	
c résolution par éléments finis Galerkin du système d'équa. diff.
c ord. non linéaires de la diffusion du moment cinétique
c par itération Newton-Raphson, formalisme de Matis & Zahn 2004, Krot=4

c	fait=0 : on forme les conditions limites
c	fait=1:
c	on forme les produits scalaire pour toutes les splines de la base
c	continue - sauf la première te la dernière - par intégration de Gauss
c	(algorithme inspiré de celui de Schumaker, 5.22 p. 203),
c	les quantités intégrées sont nulles

c	 as : dérivées/ Xi, Xi' des coeff. de la spline du prod. scalaire
c	 ad : dérivées/ Xi, Xi' des coeff. de la der. du produit scalaire

c	 as(i,j,k)=coefficient de la k-ième dérivée (k=0,1)
c	 de la j-ième variable de la i-ième équation < . Ni >

c	 ad(i,j,k)=coefficient de la k-ième dérivée (k=0,1)
c	 de la j-ième variable de la i-ième équation < . dNi/dx > 

c	 a(<spline 1 . équation 1: variable 1,2,..,ne>
c	   <spline 1 . équation 2: variable 1,2,..,ne>
c		 :     :     :     :     :     
c	   <spline 1 . équation ne: variable 1,2,..,ne>
c	   <spline 2 . équation ne: variable 1,2,..,ne>
c	   <spline 2 . équation ne: variable 1,2,..,ne>
c           :     :     :     :     :
c	   <spline 2 . équation ne: variable 1,2,..,ne>
c	   <spline 3 . équation ne: variable 1,2,..,ne>
c           :     :     :     :     :
c	   <spline dim . équation ne: variable 1,2,..,ne>)

c	 bs(i)=second membre de la i-ième équation < . Ni > 
c	 bd(i)=second membre de la i-ième équation < . d/dx Ni >

c	 b(<spline 1 . équation 1, 2,.., ne>
c	  <spline 2 . équation 1, 2,.., ne>
c           :     :     :     :
c	  <spline dims . équation 1, 2,.., nvar>)

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

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

	USE mod_donnees, ONLY : langue, lim_jpz, m_rot, nom_rot, nrot,
	1 precix, rot_min
	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), INTENT(in) :: dt
	LOGICAL, INTENT(out) :: ok
		
	REAL(kind=dp), DIMENSION(nrot,nrot,0:1) :: ad, as		
	REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: a, b, temp
	REAL(kind=dp), DIMENSION(nrot,0:1) :: y
	REAL(kind=dp), DIMENSION(0:1,m_rot) :: d	
	REAL(kind=dp), DIMENSION(nrot) :: bd, bs
	REAL(kind=dp), DIMENSION(m_rot) :: qg, wg
	REAL(kind=dp), PARAMETER :: explose=1.d2
	REAL(kind=dp), SAVE :: dtp=-HUGE(dt), err_lim, err_max
	REAL(kind=dp) :: er, err, errg, Ulim, xli
		
	INTEGER, PARAMETER, DIMENSION(5) :: iter_rot=(/ 1, 0, 0, 0, 0 /)
	INTEGER, ALLOCATABLE, DIMENSION(:) :: indpc, indpc0
	INTEGER, PARAMETER :: compt_max=10 , itest2=5
	INTEGER, SAVE :: bloc
	INTEGER :: col, colonne, compt, i, id, ie, ik, ig, ipe,
	1 iv, ive, j, k, l=1, ligne, nl, rang, spi, spl, sple

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

2000	FORMAT(8es10.3)
2001	FORMAT(10es8.1)

c initialisations
	IF(init)THEN
	 init=.FALSE.
	 	
c bloc: longueur du bloc des coeff. non idt. nuls 	
	 bloc=nrot*(2*m_rot-1)
	 
c les niveaux de précision	 
	 err_max=10.d0*precix
	 err_lim=10.d0*err_max

	ENDIF

c NOTATIONS (hélas incohérentes) pour les développements sur B-splines
c n_ch : nombre VARIABLE de points élément de mod_variables
c nch : nombre FIXE de fonctions élément de mod_donnees
c m_ch : ordre FIXE des splines élément de mod_donnees 
c mc(n_ch) : abscisses VARIABLES élément de mod_variables

c nl = rang : nombre de lignes du système linéaire	
	rang=nrot*dim_rot ; nl=rang
	
c allocations				
	ALLOCATE(a(nl,bloc),b(1,nl),indpc(nl),indpc0(nl),temp(nrot,dim_rot))

c indices de première colonne pour les produits scalaires
c dans gauss_band
	indpc0=0 ; ligne=0
	DO i=1,dim_rot      !pour chaque spline de la base des mrot
	 k=MAX(1,i-m_rot+1)      !indice s'il n'y avait qu'une variable
	 DO ie=1,nrot         !avec toutes les variables
	  ligne=ligne+1 ; indpc0(ligne)=nrot*(k-1)+1
	 ENDDO
	ENDDO
	 
c itérations Newton Raphson: boucle infinie B1
	compt=0
	B1: DO
	 
c initialisations
	 a=0.d0 ; b=0.d0 ; ligne=0 ; indpc=indpc0
	 
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*,'resout_rota4, 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 variables et dérivées
	   CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,
	1  qg(ig),l,bs,bd)
	   IF(no_croiss)PRINT*,'resout_rota4, Pb. en 1'
	   y(:,0)=bs ; y(:,1)=bd
	   
c formation des coefficients des équations
	   CALL eq_diff_rota4(0,qg(ig),y,dt,ad,as,bd,bs,compt)	    
	   		   
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 ie=1,nrot		!pour chaque équation
	    DO i=1,m_rot	!pour chaque spline d'indice i
	     ligne=nrot*(spi+i-1)+ie
	     b(1,ligne)=b(1,ligne)+wg(ig)*(d(0,i)*bs(ie)+d(1,i)*bd(ie))

c	     WRITE(*,2003)ligne,d(0,i),d(1,i),bs(ie),bd(ie),y(1,0),
c	1    y(1,1),qg(ig),b(1,ligne) ; PAUSE'b'
2003	     FORMAT(i3,9es8.1)
		
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
	      DO iv=1,nrot	!pour chaque variable	       
	       colonne=nrot*(spi+j-1)+iv    
	       col=colonne-indpc(ligne)+1   !colonne matrice compressée
c	       PRINT*,'ligne,colonne,col,iv',ligne,colonne,col,iv;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(ie,iv,id)*d(0,i)+ad(ie,iv,id)*d(1,i))*d(id,j)
	       ENDDO    !id
	      ENDDO	!iv variable
	     ENDDO	!j
	    ENDDO	!i
	   ENDDO	!ie équation	   
	  ENDDO 	!ig
	 ENDDO 		!k

c~~~~~~~~~~~~~~continuités~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c	IF(.FALSE.)THEN	!suppression des discontinuités

c continuités suivant le formalisme de Matis & Zahn 2004
c on assure la continuité de Omega, U et Phi aux limites ZR/ZC
c on égale les coefficients de spline à droite et à gauche de la discontinuité
c il faut se placer coté convectif

c continuité de Omega et Phi et conditions limites de JPZh pour U 
	IJPZ: IF(lim_jpz)THEN	
	 DO ik=1,ndis
	  CALL linf(mrot(idis(ik)),mrott,knotr,l) ; spi=l-m_rot
	  DO ie=1,nrot				!pour Omega, -U et Phi
	   SELECT CASE(ie)
	   CASE(1,2,5)		!phi continu    	   

c convectif à droite de la limite, utilisation de la spline spi+1
	    IF(lmix(mrot(idis(ik))))THEN	  
	     ligne=nrot*spi+ie
	     b(1,ligne)=rota(ie,spi)+(-1.d0)**ie*rota(ie,spi+1)
	     col=nrot*(m_rot-2)+ie ; a(ligne,:)=0.d0
	     a(ligne,col)=1.d0 ; a(ligne,col+nrot)=(-1.d0)**ie
	   
c convectif à gauche de la limite, utilisation de la spline spi
	    ELSE  	   
	     ligne=nrot*(spi-1)+ie
	     b(1,ligne)=rota(ie,spi)-rota(ie,spi+1)	    
	     col=nrot*(m_rot-1)+ie ; a(ligne,:)=0.d0
	     a(ligne,col)=1.d0 ; a(ligne,col+nrot)=-1.d0
	    ENDIF		!lmix	     
	   END SELECT	    
	  ENDDO	!ie
	 ENDDO		!ik

c conditions limites de PM
	ELSE Ijpz
	 
c pour Omega(ie=1) et Phi(ie=5)--------------------
	 DO ik=1,ndis
	  CALL linf(mrot(idis(ik)),mrott,knotr,l) ; spi=l-m_rot
	  
c pour Omega(ie=1) et Phi(ie=5)
	  DO ie=1,nrot
	   SELECT CASE(ie)
	   CASE(1,5)

c convectif à droite de la limite, utilisation de la spline spi+1
	    IF(lmix(mrot(idis(ik))))THEN	  
	     ligne=nrot*spi+ie
	     b(1,ligne)=rota(ie,spi)-rota(ie,spi+1)
	     col=nrot*(m_rot-2)+ie ; a(ligne,:)=0.d0
	     a(ligne,col)=1.d0 ; a(ligne,col+nrot)=-1.d0
	   
c convectif à gauche de la limite, utilisation de la spline spi
	    ELSE  
	     ligne=nrot*(spi-1)+ie
	     b(1,ligne)=rota(ie,spi)-rota(ie,spi+1)	    
	     col=nrot*(m_rot-1)+ie ; a(ligne,:)=0.d0
	     a(ligne,col)=1.d0 ; a(ligne,col+nrot)=-1.d0
	    ENDIF	!lmix
	   END SELECT
	  ENDDO
	 ENDDO		!ik
	
c pour le flux--------------------
c la continuité du flux repose sur des dérivées à droite et à gauche des limites
c donc sur m_rot splines de part et d'autre, le bloc devrait faire 2m_rot et
c non pas 2mrot-1. On conserve des blocs de 2mrot-1 en utilisant, dans les ZC
c la place laissée libre par U=0, pour conserver Dc dOmega/dnu.
c Avec convection à droite de la limite, le coté droit nécessite m_rot splines
	 DO ik=1,ndis
	  CALL linf(mrot(idis(ik)),mrott,knotr,l) ; spi=l-m_rot

c convectif à droite de la limite (ZC externe), utilisation de la spline spi
c à gauche de la limite, coté ZR : Omega U / 5 + Dv dOmega/dnu
c à droite de la limite, coté ZC : Uf=rota(2,spi+1)=Dc dOmega/dnu
c b(ligne)=Omega (U+Ulim) / 5 + Dv dOmega/dnu-Uf, Ulim=d Rzc /dt
c convectif à droite de la limite on passe ici pour ik= 1, 2, 4, .. ndis
c utilisation des splines 1, 2, m_rot à gauche de la limite :
	  IF(lmix(mrot(idis(ik))))THEN	  

c on se place un peu à gauche de la limite, coté ZR	  
	   xli=mrot(idis(ik))*un_eps+mrot(idis(ik)-1)*eps

c variables et dérivées
	   CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,
	1  xli,l,bs,bd)
	   IF(no_croiss)PRINT*,'resout_rota4, Pb. en 1'
	   y(:,0)=bs ; y(:,1)=bd
	   
c contribution de d/dt Rzc
c convectif à droite perte dans la ZC, gain dans la ZR
	   Ulim=+dr_zc(ik)

c formation des coefficients des équations
	   CALL eq_diff_rota4(3,xli,y,dt,ad,as,bd,bs,compt,Ulim)
	   
c les B-splines en xli, dérivées 0 à 1
	   CALL bvald(xli,mrott,m_rot,l,1,d)	   

c contribution au système linéaire, ligne de la spline spi ie. coté radiatif
	   ligne=nrot*(spi-1)+2		!spline spi, 2 pour U
	   b(1,ligne)=bs(2)-rota(2,spi+1) ; a(ligne,:)=0.d0
	   	   
	   DO j=1,m_rot		!pour chaque spline j	   
	    DO iv=1,2		!pour chaque variable
	     col=nrot*(j-1)+iv	!splines 1, 2, à gauche
	     DO id=0,1	     
	      a(ligne,col)=a(ligne,col)+as(2,iv,id)*d(id,j)	      	      
	     ENDDO	!id
	    ENDDO	!iv
	   ENDDO 	!j
	   
c contribution de Uf, spline m_rot (1, 6, 11, 16, 21)
	   col=nrot*m_rot+2 ; a(ligne,col)=-1.d0
	   
c	   PRINT*,'col pour +Uf spi+1; {b(ligne),rota(2)}',col,lmix(xli)
c	   WRITE(*,2000)bs(2),rota(2,spi-1),rota(2,spi),rota(2,spi+1)
c	   PRINT*,'ik, lmix',ik,lmix(mrot(idis(ik)))
c	   PAUSE'lmix'

	  ELSE
c convectif à gauche de la limite (noyau CV), utilisation de la spline spi+1
c à droite de la limite, coté ZR : Omega U / 5 + Dv dOmega/dnu
c à gauche de la limite, coté ZC : Uf=rota(2,spi)=Dc dOmega/dnu
c b(ligne)=Omega (U-Ulim) / 5 + Dv dOmega/dnu-Uf, Ulim=d Rzc /dt
c convectif à gauche de la limite on passe ici pour ik= (0,) 1, 3, .. ndis-1
c utilisation des splines m_rot, m_rot+1... m_rot+m_rot-1 à droite de la limite
c utilisation de la spline m_rot-1 à gauche de la limite
	   xli=mrot(idis(ik))
	   CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,
	1  xli,l,bs,bd)
	   y(:,0)=bs ; y(:,1)=bd	
	   	
c contribution de d/dt Rzc
c convectif à gauche gain dans la Zc, perte dans ZR
	   Ulim=-dr_zc(ik)
	   	  
c formation des coefficients des équations
	   CALL eq_diff_rota4(3,xli,y,dt,ad,as,bd,bs,compt,Ulim)
	   
c les B-splines en xli, dérivées 0 à 1
	   CALL bvald(xli,mrott,m_rot,l,1,d)	 

c contribution au système linéaire
	   ligne=nrot*spi+2
	   b(1,ligne)=bs(2)-rota(2,spi) ; a(ligne,:)=0.d0	   
	   	   
	   DO j=1,m_rot		!pour chaque spline j	   
	    DO iv=1,2		!pour chaque variable	       
	     col=nrot*(m_rot+j-2)+iv	!splines m_rot, m_rot+1	...     
	     DO id=0,1	         	     
	      a(ligne,col)=a(ligne,col)+as(2,iv,id)*d(id,j)      
	     ENDDO	!id
	    ENDDO	!iv
	   ENDDO	!j
	   
c contribution de Uf, spline m_rot-1	(1, 6, 11, 16, 21)
	   col=nrot*(m_rot-2)+iv ; a(ligne,col)=-1.d0
c	   PRINT*,'col pour -Uf spi; {b(ligne),rota(2)}',col,lmix(xli)
c	   WRITE(*,2000)bs(2),rota(2,spi-1),rota(2,spi),rota(2,spi+1)   
c	   PRINT*,'ik, lmix',ik,lmix(mrot(idis(ik)))	   

	  ENDIF		!lmix
	 ENDDO		!ik
	 	
	ENDIF Ijpz
	
c	ENDIF	!suppression des continuités

c~~~~~~~~~~~~~~~~~~~quantités intégrées~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c quantités intégrées sur la limite externe mrot(n_rot), toujours convectif
	CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,
	1 mrot(n_rot),l,bs,bd)
	spi=l-m_rot
	IF(no_croiss)PRINT*,'resout_rota4, Pb. en 1'
	y(:,0)=bs ; y(:,1)=bd

c formation des coefficients des équations
	CALL eq_diff_rota4(1,mrot(n_rot),y,dt,ad,as,bd,bs,compt)	    

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 avec les conditions limites de JPZ contribution de Uf, pour lim_jpz=.TRUE.
	IF(lim_jpz)THEN		!pour Uf
	 ligne=nrot*(l-1)+2
	 b(1,ligne)=b(1,ligne)-bs(2)
	 DO j=1,m_rot		!pour chaque spline j	   
	  DO iv=1,2	!pour chaque variable
	   colonne=nrot*(spi+j-1)+iv    
	   col=colonne-indpc(ligne)+1   !colonne matrice compressée
	   DO id=0,1
	    a(ligne,col)=a(ligne,col)-as(2,iv,id)*d(id,j)
	   ENDDO	!id
	  ENDDO	!iv variable
	 ENDDO	!j
	ENDIF	 

c contribution au système linéaire pour Phi
	 ligne=nrot*(l-1)+5
	 b(1,ligne)=b(1,ligne)-bs(5)
	 DO j=1,m_rot				!pour chaque spline j	   
	  colonne=nrot*(spi+j-1)+5		!dPhi seul est concerné   
	  col=colonne-indpc(ligne)+1   		!colonne matrice compressée
	  a(ligne,col)=a(ligne,col)-as(5,5,1)*d(1,j)!dPhi seul est concerné
c	  Print*,'ligne,col,{b,bs,a,d,as}',ligne,col
c	  Write(*,2000)b(1,ligne),a(ligne,col),d(1,j),as(5,5,1),bs(5)
	 ENDDO		!j
	 
c~~~~~~~~~~~~~

c	IF(.FALSE.)THEN	!sans discontinuités
	 
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 variables et dérivées
	   CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,
	1  xli,l,bs,bd)
	   spi=l-m_rot
	   IF(no_croiss)PRINT*,'resout_rota4, Pb. en 1'
	   y(:,0)=bs ; y(:,1)=bd

c formation des coefficients des équations
	   CALL eq_diff_rota4(2,xli,y,dt,ad,as,bd,bs,compt)	    

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 Omega, U et Phi sont concernés	   
	   DO ie=1,nrot	!pour chaque équation
	    SELECT CASE(ie)
	    CASE(1,2,5)
	     ligne=nrot*(l-1)+ie
	     b(1,ligne)=b(1,ligne)+bs(ie)
	     DO j=1,m_rot		!pour chaque spline j	   
	      DO iv=1,nrot	!pour chaque variable
	       colonne=nrot*(spi+j-1)+iv    
	       col=colonne-indpc(ligne)+1   !colonne matrice compressée
	       DO id=0,1
	        a(ligne,col)=a(ligne,col)+as(ie,iv,id)*d(id,j)
	       ENDDO	!id
	      ENDDO	!iv variable
	     ENDDO	!j
	    END SELECT	      
	   ENDDO 	!ie	  

c radiatif à droite de la limite, affecte la spline spi+1 avec le signe -
	  ELSE
	  
c variables et dérivées
	   CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,
	1  mrot(idis(ik)),l,bs,bd)
	   spi=l-m_rot
	   IF(no_croiss)PRINT*,'resout_rota4, Pb. en 1'
	   y(:,0)=bs ; y(:,1)=bd

c formation des coefficients des équations
	   CALL eq_diff_rota4(2,mrot(idis(ik)),y,dt,ad,as,bd,bs,compt)	    

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     
	   DO ie=1,nrot		!pour chaque équation
	    SELECT CASE(ie)
	    CASE(1,2,5)   	       
	     ligne=nrot*spi+ie
	     b(1,ligne)=b(1,ligne)-bs(ie)
	     DO j=1,m_rot		!pour chaque spline j	   
	      DO iv=1,nrot	!pour chaque variable
	       colonne=nrot*(spi+j-1)+iv    
	       col=colonne-indpc(ligne)+1   !colonne matrice compressée
	       DO id=0,1
	        a(ligne,col)=a(ligne,col)-as(ie,iv,id)*d(id,j)
	       ENDDO	!id
	      ENDDO	!iv variable
	     ENDDO	!j
	    END SELECT
	   ENDDO	!ie	   	   	   
	  ENDIF		!lmix
	 ENDDO		!ik
	 
c	 ENDIF	!sans discontinuité
	 	 	 
c----------------fin de la construction du système linéaire-------------- 

c	 PRINT*,ligne,nl,bloc ; PAUSE'ligne,nl,bloc avant a'
c	 DO i=1,nl
c	 DO i=500,510
c	 DO i=5,nl,5	 
c	  PRINT*,i,indpc(i)
c	  WRITE(*,2000)a(i,:),b(1,i)
c	  IF(mod(i,200) == 0)PAUSE'a200'
c	  IF(MAXVAL(ABS(a(i,:))) == 0.d0)PAUSE'a=0'
c	 ENDDO
c	 PAUSE'a'
c	 Write(*,2000)b(1,5:nl:5) ; PAUSE'2-nd membres'	
c	 PRINT*,nl
c	 DO i=1,nl
c	  IF(MAXVAL(ABS(a(i,:))) == 0.d0)PRINT*,i
c	 ENDDO
c	 PAUSE'lignes nulles'
		
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 resout_rota4' ; STOP
	 ENDIF
	 
c	 WRITE(*,2000)b(1,:) ; PAUSE'solution'
	
c construction du vecteur temporaire des corrections
	 temp=RESHAPE(b,SHAPE(temp))
	 
c	  B5: DO i=1,dim_rot
c	  PRINT*,i
c	  WRITE(*,2000)temp(1,i),rota(1,i),temp(3,i),rota(3,i)
c	  WRITE(*,2000)temp(:,i)
c	  WRITE(*,2000)rota(:,i)
c	  IF(MINVAL(ABS(rota(:,i))) == 0.d0)CYCLE B5
c	  WRITE(*,2000)(temp(j,i)/rota(j,i),j=1,nrot)	  
c	  WRITE(*,2000)rota(:,i)
c	  IF(mod(i,200) == 0)PAUSE'corrections'
c	 ENDDO B5
c	 PAUSE'corrections'
	 	
c estimation de la précision
	 err=0.d0 ; errg=0.d0
	 DO spl=1,dim_rot
	  B4: DO iv=1,nrot
	   IF(ABS(rota(iv,spl)) < rot_min(iv))CYCLE B4
	   er=ABS(temp(iv,spl)/rota(iv,spl))
	   IF(er > err)THEN
	    sple=spl ; ive=iv ; err=er
c	    PRINT*,iv,spl,nom_rot(iv)
c	    WRITE(*,2000)er,err,temp(iv,spl),rota(iv,spl),rot_min(iv)
	   ENDIF
	   IF(COUNT(iv == iter_rot) > 0)errg=MAX(errg,er)
	  ENDDO B4	!var
	 ENDDO		!spl	 
c	 PRINT*,sple,ive ; WRITE(*,2000)er,err ; PAUSE
	 
c écritures
	 ipe=sple/m_rot+1 ; compt=compt+1
	 WRITE(*,100)compt,errg,nom_rot(1),nom_rot(itest2),err,
	1 nom_rot(ive),sple,ipe,mrot(ipe)**(2.d0/3.d0)
100	 FORMAT('resout_rota4, itération:',i3,', erreur max:',es8.1,
	1 ', sur ',a,' ou ',a,/,'err. max. globale:',es8.1,
	2 ', sur ',a,', B-spline:',i4,', couche:',i4,', m=',es10.3)
	 WRITE(*,101)mrot(ipe),nom_rot,rota(:,sple),temp(:,sple)
101	 FORMAT(/,'variables et corrections au maximum d''erreur, nu=',es10.3,/,
	1 5(3x,a4,3x),/,5es10.3,/,5es10.3,/)	
	 
c	 DO i=1,n_rot
c	  CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,mrot(i),
c	1  l,bs,bd)
c	  WRITE(*,2000)SQRT(mrot(i))**3,bs
c	 ENDDO	 
c	 PAUSE'S'

c	 PAUSE'avant écritures'
c	 DO i=1,dim_rot
c	  PRINT*,i
c	  WRITE(*,2000)temp(1,i),rota(1,i),temp(3,i),rota(3,i)
c	  WRITE(*,2000)temp(:,i)
c	  WRITE(*,2000)rota(:,i)
c	  IF(mod(i,200) == 0)PAUSE'corrections'
c	 ENDDO
c	 PAUSE'écritures'

c limitation des corrections
	 DO spl=1,dim_rot
	  B5: DO iv=1,nrot	  
	   IF(ABS(rota(iv,spl)-temp(iv,spl)) < rot_min(iv))CYCLE B5
	   temp(iv,spl)=SIGN(MIN(ABS(temp(iv,spl)),0.6d0*ABS(rota(iv,spl))),
	1  temp(iv,spl))
	  ENDDO B5	!iv
	 ENDDO		!spl

c corrections
	 rota=rota-temp

c mises à 0 en dessous de 1.d-20
	WHERE(ABS(rota) < 1.d-20)rota=0.d0
	 
c Si le centre est convectif U=0, mais on y garde Dc d/dnu Omega
c Psi et Phi sont toujours nuls au centre
	 rota(3,1)=0.d0 ; rota(5,1)=0.d0	 

c gestion de la précision	 
c on poursuit si l'erreur est < err_lim après compt_max 
c il y a convergence forcée si l'erreur est < err_max après 2*compt_max
c pas de test de précision pour le premier pas temporel
c sortie avec une erreur > explose
c	 PRINT*,compt,compt_max
c	 WRITE(*,2000)errg,precix,err_max,err_lim	 
	 	 	 
	 IF(compt <= 1)CYCLE B1	 
	 ok= errg <= precix	 	  
	 IF(ok .OR. errg > explose)THEN
	  EXIT B1	  
	 ELSEIF(compt < compt_max)THEN
	  CYCLE B1	 
	 ELSEIF(errg < err_max)THEN
	  SELECT CASE(langue)
	  CASE('english')
	   PRINT*,'bad conv. in resout_rota4, number of iterations :',compt
	  CASE DEFAULT
	   PRINT*,'conv. forcée dans resout_rota4, nb. itérations :',compt
	  END SELECT
	  ok=.TRUE. ; EXIT B1
	 ELSEIF(errg < err_lim .AND. compt < 2*compt_max)THEN
	  CYCLE B1
	 ELSEIF(dt /= dtp)THEN
	  ok=.TRUE. ; EXIT B1 	  
	 ELSE  
	  SELECT CASE(langue)
	  CASE('english')
	   PRINT*,'no conv. in resout_rota4, number of iterations :',compt  
	  CASE DEFAULT
	   PRINT*,'pas de conv. dans resout_rota4, nb. itérations :',compt
	  END SELECT
	  EXIT B1
	 ENDIF
	ENDDO B1
	dtp=dt

c	DO i=1,n_rot
c	 CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE.,mrot(i),
c	1 l,bs,bd)
c	 WRITE(*,2000)mrot(i),bs
c	ENDDO
c	PAUSE'solution'	

c nettoyage		
	DEALLOCATE(a,b,indpc,indpc0,temp)

c	PAUSE'sortie resout_rota4'
	
	RETURN
	
	END SUBROUTINE resout_rota4