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

	SUBROUTINE resout_chim(dt,ok)

c routine private 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 microscopique  par itération Newton-Raphson

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 : Krot, langue, m_ch, nchim,
	1 nom_elem, nucleo, precic, precix
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, bvald, gauss_band, intgauss,
	1 linf, no_croiss
	USE mod_variables, ONLY : chim, dim_ch, knotc, mc, mct, n_ch
	
	IMPLICIT NONE

	REAL (kind=dp), INTENT(in) :: dt
	LOGICAL, INTENT(out) :: ok
		
	REAL(kind=dp), DIMENSION(nchim,nchim,0:1) :: ad, as		
	REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: a, b, temp
	REAL(kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: qg, wg	
	REAL(kind=dp), DIMENSION(nchim,0:1) :: y
	REAL(kind=dp), DIMENSION(0:1,m_ch) :: d	
	REAL(kind=dp), DIMENSION(nchim) :: bd, bs
	
	REAL(kind=dp), SAVE :: dtp=-HUGE(dt), err_lim, err_max, preci	
	REAL(kind=dp) :: corr, er, errg
		
	INTEGER, PARAMETER :: compt_max=15
	INTEGER, ALLOCATABLE, DIMENSION(:) :: indpc, indpc0
	INTEGER, SAVE :: bloc, mg
	INTEGER :: col, colonne, compt, i, id, ie, ig, ipe, iv, ive, j, k,
	1 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
	
c bloc: longueur du bloc des coeff. non idt. nuls	 
	 bloc=nchim*(2*m_ch-1)
	 err_max=10.d0*precic
	 
c nombre de points pour l'intégration de Gauss
	 mg=NINT(REAL(m_ch))/2.d0+1
	 ALLOCATE(qg(mg),wg(mg))	 

c gestion de la précision	 
	 err_lim=10.d0*err_max	 
	 init=.FALSE.
	 SELECT CASE(Krot)
	 CASE(3,4)
	  preci=precic*10.d0
	 CASE DEFAULT
	  preci=precic	 
	 END SELECT
	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=nchim*dim_ch ; nl=rang
		
	ALLOCATE(a(nl,bloc),b(1,nl),indpc(nl),indpc0(nl),temp(nchim,dim_ch))

c indices de première colonne pour les produits scalaires
c ils seront changés dans gauss_band
	indpc0=0 ; ligne=0
	DO i=1,dim_ch      !pour chaque spline de la base des mc
	 k=MAX(1,i-m_ch+1)      !indice s'il n'y avait qu'une variable
	 DO ie=1,nchim         !avec toutes les variables
	  ligne=ligne+1 ; indpc0(ligne)=nchim*(k-1)+1	  
	 ENDDO
	ENDDO
	 
c itérations Newton Raphson: boucle infinie B1
	compt=0
	B1: DO compt=1,compt_max
	 indpc=indpc0

c initialisations
	 a=0.d0 ; b=0.d0 ; ligne=0
	 
c---------------début de la construction du sytème linéaire------------	
	    
c les produits scalaires, pour chaque intervalle de m_ch à dim_ch
	 B3: DO k=m_ch,dim_ch

c on évite les points multiples
	  IF(mct(k) >= mct(k+1))CYCLE B3	 
	  	 
c spi : indice-1 de la première spline non id.nulle  
c comme mct(l) <= qg(ig) < mct(l+1), l est le même pour toutes
c les abscisses de Gauss dans l'intervalle [mct(k),mct(k+1)]
	  CALL linf(mct(k),mct,knotc,l) ; spi=l-m_ch
	  
c poids, wg et abscisses, qg pour intégration Gauss d'ordre m_ch
	  CALL intgauss(mct(k),mct(k+1),qg,wg,mg)

c pour chaque abscisse de GAUSS
	  DO ig=1,mg

c variables et dérivées
	   CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,qg(ig),l,bs,bd)
	   IF(no_croiss)PRINT*,'resout_chim, Pb. en 1'
	   y(:,0)=bs ; y(:,1)=bd
	   
c formation des coefficients des équations
	   CALL eq_diff_chim(qg(ig),y,dt,ad,as,bd,bs)
	  	   		   
c les B-splines en qg(ig), dérivées 0 à 1
	   CALL bvald(qg(ig),mct,m_ch,l,1,d)   
	   	   
c contribution au système linéaire
	   DO ie=1,nchim		!pour chaque équation
	    DO i=1,m_ch		!pour chaque spline d'indice i
	     ligne=nchim*(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_ch		!pour chaque spline j
	      DO iv=1,nchim	!pour chaque variable	       
	       colonne=nchim*(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)*d(id,j)+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 B3	!k
	 
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	  PRINT*,i,indpc(i)
c	  WRITE(*,2000)a(i,:),b(1,i) ; IF(mod(i,200) == 0)PAUSE'a200'
c	  IF(MAXVAL(ABS(a(i,:))) == 0.d0)PAUSE'a=0'
c	 ENDDO
c	 PAUSE'a'
c	ENDIF	
	
c	PRINT*,nl
c	DO i=1,nl
c	 IF(MAXVAL(ABS(a(i,:))) == 0.d0)THEN
c	  PRINT*,'ligne ',i,' nulle'
c	  WRITE(*,2000)a(i,1:7),b(1,i)
c	  PAUSE'ligne nulle'
c	 ENDIF
c	ENDDO
 	
c résolution du système linéaire
	 CALL gauss_band(a,b,indpc,nl,rang,bloc,1,inversible)
	 IF(.NOT.inversible)THEN
	  PRINT*,'ARRET, matrice singulière dans resout_chim' ; STOP
	 ENDIF
	 
c	 WRITE(*,2000)b(1,:) ; PAUSE'solution'
	
c vecteur des corrections
	 temp=RESHAPE(b,SHAPE(temp))
	 
c	 DO i=1,dim_ch!,100
c	  PRINT*,i
c	  WRITE(*,2000)temp(1,i),chim(1,i),temp(3,i),chim(3,i)
c	  WRITE(*,2000)chim(:,i)
c	  IF(mod(i,50) == 0)PAUSE'solution'
c	 ENDDO
c	 PAUSE'corrections au pas de 100'

c estimation de la précision NR/gramme
	 errg=0.d0 ; corr=1.d0
	 DO spl=1,dim_ch
	  DO iv=1,nchim
	   er=ABS(temp(iv,spl))*nucleo(iv)	   
	   IF(er > errg)THEN
	    sple=spl ; ive=iv ; errg=er
c	    PRINT*,iv,spl,nom_elem(iv)
	   ENDIF
	   
c limitation des corrections	   
	   IF(corr*ABS(temp(iv,spl)) > 0.6d0*ABS(chim(iv,spl)))corr=corr/2.d0
	  ENDDO		!chim
	 ENDDO		!spl
	 
c	 PRINT*,sple,ive ; WRITE(*,2000)er,err ; PAUSE

	corr=1.d0

c nouvelle solution
	 chim=chim-temp*corr
	 
c écritures
	 ipe=sple/m_ch+1	!indice de couche
	 WRITE(usl_evol,100)compt,errg,nom_elem(ive),corr,
	1 sple,ipe,mc(ipe)**(2.d0/3.d0)	
100	 FORMAT('resout_chim, itération:',i3,', corr/gramme max:',es8.1,
	1 ', sur ',a,/,'corr:',es8.1,', B-spline:',i4,', couche:',i4,
	2 ', m=',es10.3)
	 
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 on accepte une mauvaise convergence pour un nouveau pas temporel
	 ok= errg <= preci
	 IF(compt <= 1)THEN
	  CYCLE B1
	 ELSEIF(ok)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_chim, number of iterations :',compt
	  CASE DEFAULT
	   PRINT*,'conv. forcée dans resout_chim, 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_chim, number of iterations :',compt
	  CASE DEFAULT
	   PRINT*,'pas de conv. dans resout_chim, nb. itérations :',compt
	  END SELECT
	  EXIT B1
	 ENDIF
	ENDDO B1
	dtp=dt

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

c	PAUSE'sortie'
	
	RETURN
	
	END SUBROUTINE resout_chim