!-------------------------------------------------------------------------------------
!					MODULES
!--------------------------------------------------------------------------------------
!    Constantes du programme
	module constante
	   integer,parameter			:: nd1=20
	   integer,parameter			:: nd2=2
	   integer,parameter			:: ncm=4	! Nb options calcul
	   integer,parameter			:: nca=122	! Longueur chaine car.
	end module constante
!    ----------------------------------------------------------------------------------
!    Parametres de calcul
	module parametre
	   use constante,only: ncm,nca
	   integer				:: n,m,niter,niterp,nopt,nbloc,c,b,nv,kx,ky
	   real(8)				:: beta,rmu,wmin,wmax,u,remp
	   real(8)				:: alpha1,alpha2,ski
	   real(8),parameter			:: eps=1.d-8
	   character(len=nca)			:: ficpara='mesh.dat'
	   character(len=nca)			:: ficin,ficout,ficres
	   character(len=nca)			:: ficcov,repweb
	   character(len=1),dimension(ncm)	:: ynmes
	   integer,dimension(500,2)		:: v
	end module parametre
!--------------------------------------------------------------------------------------
!				   PROGRAMME PRINCIPAL
!--------------------------------------------------------------------------------------
	Program meshmem
!--------------------------------------------------------------------------------------
! Nom	     :	meshmem	[MAIN]							0
!	
! Version    : 	4_4.0 Beta (29 juin 1999)
!
! Fonction   :	Methode de l'entropie maximum selon l'algorithme
!		de Meshkov et Berkov
!
! Auteur(s)  : 	H. Touchette selon la version de S.Moukouri
!
! Revision   :  29 juillet 1999 D. Poulin
!
! Adresse    : 	Centre de recherche en physique du solide (CRPS)
!		Departement de physique, Universite de Sherbrooke
!		Sherbrooke, Quebec, Canada, J1K 2R1
!
! References :  [1] Meshkov et Berkov, Int. J. Mod. Phys. C 5, 987 (1994)
!		[2] H.Touchette, Aspects numeriques des simulations Monte Carlo du
!		    modele de Hubbard, rapport technique du CRPS, ete 1997.
!		[3] Jarrell et Gubernatis, Physics Reports 269 (1996)
!		[4] D. Poulin et H. Touchette, Aspects numériques des simulations Monte
!		    Carlo quantique du modele de Hubbard, Rapport technique du CRPS, 1999.
!
! Notes: 	- Version avec le moyennage sur le parametre de regularisation (alpha)
!		- Traitement des donnees et du kernel avec la matrice de covariance
!		  (optionnel)
!		- Les donnees qui s'etalent de 0 a TRAN (TRAN+1 valeurs) sont stockees
!		  dans g(k), k=1..TRAN+1  (TRAN = n)
!		- Options de calcul
!		  YNMES(1) : Matrice de covariance
!		  YNMES(2) : Imposition de la symétrie particule-trou
!		  YNMES(3) : Doubler le nombre de mesure par symétrie élémentaire
!                 YNMES(4) : Prolongement d'une susceptibilité
!
!--------------------------------------------------------------------------------------
	use constante,only: nd1,nd2
	use parametre
	implicit none
!     Variables scalaires
	real(8)				:: ti,tf,td
	real(8)				:: chi2,dw,entr,dtau
	real(8)				:: sm
	real(8)				:: alpha,dalpha
	real(8)				:: da,ppr
	real(8)				:: prmax,almax,smax,chi2max
	integer				:: i,j,k,l
	integer				:: iter,iterp,iterv
!     Fonction
	character(8)			:: date
	integer				:: mclock
	real(8)				:: rtsafe,delta
!     Tableaux
	integer,allocatable		:: it(:)
	real(8),allocatable		:: vgam(:),vgam2(:),a(:),w(:)
	real(8),allocatable		:: d(:),e(:),ap(:),as(:)
	real(8),allocatable		:: veca(:),phi(:),phib(:),dm(:)
	real(8),allocatable		:: t(:),g(:),err(:)
	real(8),allocatable		:: gami(:,:),selfm(:,:),ek(:,:)
	real(8),allocatable		:: gama(:,:),dmata(:,:)
	real(8),allocatable		:: covm(:,:),covd(:),cove(:)
	real(8),allocatable		:: aa(:,:),am(:),pr(:),lpr(:),al(:)
	real(8),allocatable		:: ss0(:),ss1(:),ss2(:)
!--------------------------------------------------------------------------------------
!
	ti=mclock()
!
!Lecture parametres de calcul
	call lecture_para(n,g,err)
!
!Ouverture des fichiers d'ecriture
	open(unit=7,file=ficout)
	open(unit=9,file=ficres)
	rewind(unit=7)
	rewind(unit=9)
	
	do 222 iterv=1,nv
	kx=v(iterv,1)
	ky=v(iterv,2)
	write(7,'(A)')'#-------------------------------'
	write(7,fmt='(2(A,I3))')'#kx=',kx,'  ky=',ky
	write(9,'(A)')'#-------------------------------'
	write(9,fmt='(2(A,I3))')'#kx=',kx,'  ky=',ky
!
!Declaration dynamique des tableaux
	allocate (t(n),it(n),g(n),err(n))
	allocate (vgam(m),vgam2(m),a(m),w(m),d(m),e(m),ap(m),as(m))
	allocate (veca(m),phi(m),phib(m),dm(m))
	allocate (ek(n,m),gami(m,m),selfm(m,m),gama(m,m),dmata(m,m))
	allocate (covm(n,n),covd(n),cove(n))
	allocate (aa(niterp,m),am(m),pr(niterp),lpr(niterp),al(niterp))
	allocate (ss0(niterp),ss1(niterp),ss2(niterp))
!
!Lecture donnees
	call lecture_donnees(it,g,err,iterv)
!
!Traitement des donnees par la matrice de covariance
	if(ynmes(1).eq.'Y'.or.ynmes(1).eq.'y') then
	   call lecture_cova(n,covm)
	   call tred(covm,n,n,covd,cove)
	   call tqli(covd,cove,n,n,covm)
	   call invmg(n,n,covm)
           if(covd(n)<1D-10) covd(n)=1D-10
           do i=n,1,-1
              if(covd(i)<1D-10) covd(i)=covd(i+1)
           end do
           g=matmul(covm,g)
           do i=1,n
              err(i)=sqrt(covd(i))
           end do
           write(9,*)'--valeurs propres--'
           do i=1,n
              write(9,*)g(i),'		',err(i)
           end do
           write(9,*) '--------------------------'
	endif
!
!Paramètres de calcul
	dw=wmax-wmin
	dw=dw/m
	dtau=beta/(n)
	write(*,*) 'dtau=',dtau
!
!controle de l'erreur
	do k=1,n
           t(k)=it(k)*dtau
	end do
	do i=1,m
	   w(i)=wmin+(i-1)*dw-rmu
	   dm(i)=1.d0/((m-1))
	end do
!
	alpha=alpha1
	dalpha=(alpha2-alpha1)/niterp
	do iterp=1,niterp
	   alpha=alpha+dalpha
	   al(iterp)=alpha
           ppr=0
!
!modèle par défaut
	   do i=1,m
	      a(i)=dm(i)
              ppr=ppr+a(i)
	   end do
           do i=1,m
              a(i)=a(i)/ppr
              ap(i)=a(i)
           end do
!
!Noyau de la transformation
	   do i=1,m
	      do k=1,n
	         ek(k,i)=dexp(-w(i)*t(k))/(1.d0+dexp(-beta*w(i)))
	      end do
	   end do
!
	   if(ynmes(1).eq.'Y'.or.ynmes(1).eq.'y') ek=matmul(covm,ek)
	   do iter=1,niter
	      do i=1,m
	         as(i)=ski*a(i)+(1.d0-ski)*ap(i)+eps
	      end do
	      do i=1,m
	         sm=0.d0
		 do k=1,n
		    sm=sm+g(k)*dw*ek(k,i)/err(k)**2
		 end do
		 veca(i)=sm-alpha*dw*(dlog(as(i)/dm(i))-1.d0)   !F'(A(n))
	      end do
	      do j=1,m
		 do i=1,m
		    sm=0.d0
		    do l=1,n
		       sm=sm+dw*dw*ek(l,i)*ek(l,j)/err(l)**2
		    end do
		    dmata(i,j)=sm+alpha*dw*delta(i,j)/as(j)     !F''(A(n))
		    selfm(i,j)=sm
		 end do
	      end do
	      call tred(dmata,m,m,d,e)
	      call tqli(d,e,m,m,dmata)
	      do j=1,m
		 sm=0.d0
		 do i=1,m
		    sm=sm+dmata(i,j)*veca(i)
		 end do
		 phib(j)=sm/dsqrt(d(j))
	      end do
	      do j=1,m
		 do i=1,m
		    gama(i,j)=dmata(i,j)/dsqrt(d(j))
		    gami(i,j)=dsqrt(d(i))*dmata(j,i)
		 end do
	      end do
	      do i=1,m
	         sm=0.d0
		 do j=1,m
		    sm=sm-gama(i,j)*phib(j)
		 end do
		 vgam(i)=sm
	      end do
	      do i=1,m
		 sm=0.d0
		 do j=1,m
		    sm=sm-gami(i,j)*vgam(j)
		 end do
		 vgam2(i)=sm
	      end do
	      call nnls(gami,m,m,m,vgam2,phi)
	      ss0(iterp)=0.d0
	      ss1(iterp)=0.d0
	      ss2(iterp)=0.d0
	      do i=1,m
		 a(i)=ap(i)
		 ap(i)=phi(i)
		 ss0(iterp)=ss0(iterp)+ap(i)*dw
		 ss1(iterp)=ss1(iterp)+ap(i)*w(i)*dw
		 ss2(iterp)=ss2(iterp)+ap(i)*w(i)*w(i)*dw
	      end do
	   end do
!
	   entr=0.d0
	   do i=1,m
	      entr=entr+(dm(i)-as(i)+as(i)*dlog(as(i)/dm(i)))*dw	!#
	   end do
 	   chi2=0.d0
	   do k=1,n
	      sm=0.d0
	      do i=1,m
		 sm=sm+dw*ek(k,i)*phi(i)
	      end do
	      chi2=chi2+.5d0*((g(k)-sm)/err(k))**2
	   end do
	   do i=1,m
	      do j=1,m
		 selfm(i,j)=selfm(i,j)*dsqrt(as(i)*as(j))*dw		!#
	      end do
	   end do
	   call tred(selfm,m,m,d,e)
	   call tqli(d,e,m,m,selfm)
	   
	   da=1.d0
	   do i=1,m
	      aa(iterp,i)=as(i)
	      da=da*dsqrt(alpha/(alpha+d(i)))
	   end do
	   lpr(iterp)=-(chi2+alpha*entr)+log(da)-log(alpha)
	   if (iterp==1) then
	      prmax=lpr(iterp)
	      almax=alpha
              chi2max=chi2*2
	      smax=-entr*alpha
	   else
	      if (lpr(iterp)>prmax) then
	         prmax=lpr(iterp)
	         almax=alpha
		 chi2max=chi2*2
		 smax=-entr*alpha
	      end if
	   end if
	end do
	write(9,*)'almax                   chi2 max              S max* almax'
	write(9,*)almax,'	',chi2max,'	',smax
	write(9,*)'---------------------------------------------------'
!
	am=0.d0
	do iterp=1,niterp
	   pr(iterp)=0.d0
	   do i=1,niterp
	      pr(iterp)=pr(iterp)+exp(lpr(i)-lpr(iterp))
	   end do
	   pr(iterp)=1.d0/pr(iterp)
	   do i=1,m
	      am(i)=am(i)+pr(iterp)*aa(iterp,i)
	   end do
	end do
	call sortie(w,am,al,pr,ss0,ss1,ss2)
!
	tf=mclock()
	td=(tf-ti)/6000
	write(*,*) td,'Temps de calcul (min CPU)'
!
!	Deallocation dynamique des tableaux
	deallocate (t,it,ek,g,err)
	deallocate (vgam,vgam2,a,w,d,e,ap,as)
	deallocate (veca,phi,phib,dm)
	deallocate (gami,selfm,gama,dmata)
	deallocate (covm,covd,cove)
	deallocate (aa,am,pr,al)
	
222	end do

	close(unit=7)
	close(unit=9)
!
33	format(4(1x,f10.5))
	end program meshmem
!--------------------------------------------------------------------------------------
!					ROUTINES
!--------------------------------------------------------------------------------------
! Nom		: lecture_para							1
! Fonction	: Lecture des parametres de calcul
! Auteur(s)	: S. Allen
! Revision	: 12 mars 1998 D. Poulin
! Notes		:
!--------------------------------------------------------------------------------------
SUBROUTINE LECTURE_PARA
        USE CONSTANTE
        USE PARAMETRE
	IMPLICIT NONE
	INTEGER				:: I,I1,I2,I3,I4,I5,I6,I7,I8,I9,I10
	CHARACTER(LEN=NCA)		:: TP1,TP2,TP3,TP4,TP5,TP6,TP7,TP10
	CHARACTER(LEN=NCA)		:: PREFIX='mes',NANUM,NANUM2,BUFFER
	CHARACTER(LEN=10)		:: NUMTOASCII
        OPEN(UNIT=4,FILE=FICPARA,STATUS='OLD',ACCESS='SEQUENTIAL',POSITION=&
            &'REWIND',ERR=110)
        READ(4,*) TP2
        READ(4,*) BETA
        READ(4,*) N
        READ(4,*) U
        READ(4,*) RMU
        READ(4,*) REMP
	READ(4,*) C
	READ(4,*) B
	READ(4,*) M
	READ(4,*) WMIN
	READ(4,*) WMAX
	READ(4,*) NITERP
	READ(4,*) ALPHA1
	READ(4,*) ALPHA2
	READ(4,*) NITER
	READ(4,*) SKI
	READ(4,*) NV
	READ(4,*) NOPT
	DO I=1,NOPT
	    READ(4,*) YNMES(I)
	END DO
	READ(4,*) BUFFER
	READ(4,*) BUFFER
	DO I=1,NV
           READ(4,*) V(I,1),V(I,2)
        END DO
        CLOSE(UNIT=4)
!
TP1=NUMTOASCII(INT(ABS(U)))
	TP3=NUMTOASCII(C)
	TP4=NUMTOASCII(B)
	TP5=NUMTOASCII(INT(BETA))
	TP6=NUMTOASCII(KX)
	TP7=NUMTOASCII(KY)
	TP10=NUMTOASCII(INT(100*REMP))
	I1=INDEX(TP1,' ')
	I2=INDEX(TP2,' ')
	I3=INDEX(TP3,' ')
	I4=INDEX(TP4,' ')
	I5=INDEX(TP5,' ')
	I6=INDEX(PREFIX,' ')
	I8=INDEX(TP6,' ')-1
	I9=INDEX(TP7,' ')-1
	I10=INDEX(TP10,' ')
	NANUM=TP2(1:I2-1)//'.'//PREFIX(1:I6-1)//TP1(1:I1-1)//TP10(1:I10-1)//&
	    &TP3(1:I3-1)//TP4(1:I4-1)//TP5(1:I5-1)//' '
	NANUM2=TP2(1:I2-1)//'.'//'sae'//TP1(1:I1-1)//TP10(1:I10-1)//&
	    &TP3(1:I3-1)//TP4(1:I4-1)//TP5(1:I5-1)//' '
	I7=INDEX(NANUM,' ')-1
	IF(YNMES(4)=='Y'.OR.YNMES(4)=='y') THEN
	   FICIN=NANUM2(1:I7)//'re.dat'
	ELSE
	   FICIN=NANUM2(1:I7)//'gi.dat'
	END IF
	FICOUT=NANUM(1:I7)//'ps'//'.dat'
	FICCOV=NANUM2(1:I7)//'cv.dat'
	FICRES=NANUM(1:I7)//'al.dat'
!
	OPEN(UNIT=33,FILE=FICRES,POSITION='REWIND')
	   WRITE(33,*) TP2,'	germe'
	   WRITE(33,*) BETA,'	beta'
           WRITE(33,*) N,'	ntran'
           WRITE(33,*) U,'	U'
           WRITE(33,*) RMU,'	mu'
	   WRITE(33,*) NBLOC,'	nbloc'
	   WRITE(33,*) C,'	dimension C'
	   WRITE(33,*) B,'	dimension B'
	   WRITE(33,*) M,'	nb w'
	   WRITE(33,*) WMIN,'	w min'
	   WRITE(33,*) WMAX,'	w max'
	   WRITE(33,*) NITERP,'	niterp'
	   WRITE(33,*) ALPHA1,'	alpha1'
	   WRITE(33,*) ALPHA2,'	alpha2'
	   WRITE(33,*) NITER,'	niter'
	   WRITE(33,*) SKI,'	ski'
	CLOSE(UNIT=33)
	
        RETURN
110     WRITE(*,*) 'Ne peut ouvrir le fichier FICIN'; STOP
        END SUBROUTINE LECTURE_PARA

!--------------------------------------------------------------------------------------
! Nom		: lecture_donnees	[sub]					2
! Fonction	: Lecture des donnees
! Auteur(s)	: D. Poulin
! Notes		:
!--------------------------------------------------------------------------------------


	SUBROUTINE LECTURE_DONNEES(IT,G,ERR,ITERV)
	USE CONSTANTE
	USE PARAMETRE
	IMPLICIT NONE
	INTEGER					:: iterv,NB,KP
	INTEGER					:: VX,VY,T
	INTEGER,DIMENSION(4)			:: QX,QY
	INTEGER,INTENT(OUT),DIMENSION(N+1)	:: IT
	REAL(8)					:: GL,EL
	REAL(8),DIMENSION(N+1,4)		:: GR,ERRT
	REAL(8),DIMENSION(N)			:: G,ERR
	REAL(8),DIMENSION(N+1)			:: GSYM
	CHARACTER(LEN=NCA)			:: LIGNE
	CHARACTER(LEN=1),DIMENSION(NCM)		:: YNMS
	CHARACTER(LEN=9)			:: BUFFER
	INTEGER					:: NUMY,L,I,J,K,NO,MESC,BL,P
	REAL(8)					:: TEMP,SUM,SUM2

	QX(1)=KX
	QY(1)=KY
	NB=1
	IF(YNMES(2)=='y'.OR.YNMES(2)=='Y') THEN
	   IF(YNMES(3)=='y'.OR.YNMES(3)=='Y') THEN
	      QX(2)=C-KX
	      QY(2)=B-KY
	      IF(QX(2)==0) QX(2)=C
	      IF(QY(2)==0) QY(2)=B
	      QX(3)=MOD(KX+C/2,C)
	      QY(3)=MOD(KY+B/2,B)
	      IF(QX(3)==0) QX(3)=C
	      IF(QY(3)==0) QY(3)=B
	      IF(C==1) QX(3)=1
	      IF(B==1) QY(3)=1
	      QX(4)=C-QX(3)
	      QY(4)=B-QY(3)
	      IF(QX(4)==0) QX(4)=C
	      IF(QY(4)==0) QY(4)=B
	      NB=4
	   ELSE
	      QX(2)=MOD(KX+C/2,C)
	      QY(2)=MOD(KY+B/2,B)
	      IF(QX(2)==0) QX(2)=C
	      IF(QY(2)==0) QY(2)=B
	      IF(C==1) QX(2)=1
	      IF(B==1) QY(2)=1
	      NB=2
	   END IF
	ELSE
	   IF(YNMES(3)=='y'.OR.YNMES(3)=='Y') THEN
	      QX(2)=C-KX
	      QY(2)=B-KY
	      IF(QX(2)==0) QX(2)=C
	      IF(QY(2)==0) QY(2)=B
	      NB=2
	   END IF	   
	END IF

	DO I=1,NB
	   OPEN(UNIT=11,FILE=FICIN,STATUS='OLD',ACCESS='SEQUENTIAL',POSITION=&
	       &'REWIND',ERR=100)
	   IF(YNMES(4)=='Y'.OR.YNMES(4)=='y') THEN
	      DO
	         READ(11,*) BUFFER
	         IF(BUFFER(1:6)=='SUSMAG') EXIT
	      END DO
	      READ(11,*) BUFFER
	      READ(11,*) BUFFER
	   ELSE
	      READ(11,*) LIGNE
	      READ(11,*) NO
	      READ(11,*) MESC
	      DO J=1,6
	          READ(11,*) BUFFER
	      END DO
	   END IF
	   DO
	      READ(11,*) VY,VX,T,GL,EL
	      IF(VX==QX(I).AND.VY==QY(I)) EXIT
	   END DO
	   GR(1,I)=GL
	   ERRT(1,I)=EL
	   IT(1)=T
	   DO K=2,N+1
	      READ(11,*) NUMY,NUMY,IT(K),GR(K,I),ERRT(K,I)
	   END DO
	   CLOSE(UNIT=11)
	END DO
        
        IF(YNMES(3)=='Y'.OR.YNMES(3)=='y') THEN
        IF(YNMES(2)=='Y'.OR.YNMES(2)=='y') THEN
	   DO K=1,N
	      KP=N+2-K
	      G(K)=(GR(K,1)+GR(K,2)+GR(KP,3)+GR(KP,4))/4
	      SUM=(GR(K,1)**2+GR(K,2)**2+GR(KP,3)**2+GR(KP,4)**2)/4
	      SUM2=(ERRT(K,1)**2+ERRT(K,2)**2+ERRT(KP,3)**2+ERRT(KP,4)**2)/4
	      ERR(K)=DSQRT(SUM+SUM2-G(K)**2)
	   END DO
	ELSE
           DO K=1,N
              G(K)=(GR(K,1)+GR(K,2))/2
              ERR(K)=(ERRT(K,1)**2+ERRT(K,2)**2+GR(K,1)**2+GR(K,2)**2)/2-G(K)**2
              ERR(K)=SQRT(ERR(K))
           END DO
        END IF
        ELSE
        IF(YNMES(2)=='Y'.OR.YNMES(2)=='y') THEN
	   DO K=1,N
	      KP=N+2-K
	      G(K)=(GR(K,1)+GR(KP,2))/2
	      SUM=(GR(K,1)**2+GR(KP,2)**2)/2
	      SUM2=(ERRT(K,1)**2+ERRT(KP,2)**2)/2
	      ERR(K)=DSQRT(SUM+SUM2-G(K)**2)
	   END DO
	ELSE
           DO K=1,N
              G(K)=GR(K,1)
              ERR(K)=SQRT(ERRT(K,1))
           END DO
        END IF
        END IF
!
	IF(ITERV==1) THEN
	   WRITE(7,fmt='(A)')'#-------------------------------------------------------------'
	   WRITE(7,fmt='(A,I5,A)') '#',NO,'    No. simulation'
	   WRITE(7,fmt='(A,I8,A)') '#',MESC,'      Nb mesures effectuees'
	   WRITE(7,fmt='(A,A,A)') '# ',YNMES(1),'	MATRICE DE COVARIANCE'
	   WRITE(7,fmt='(A,A,A)') '# ',YNMES(2),'	SYMÉTRIE DE G'
	   WRITE(7,fmt='(A,A,A)') '# ',YNMES(3),'	DOUBLER LE NOMBRE DE MESURES'
	END IF
	CLOSE(UNIT=11)
	RETURN
!    Messages d'erreurs
100	WRITE(*,*) 'Erreur : Ne peut ouvrir le fichier ',FICIN; STOP
	END SUBROUTINE LECTURE_DONNEES
!
!--------------------------------------------------------------------------------------
! Nom		: lecture_cova		[sub]					3
! Fonction	: Lecture de la matrice de covariance
! Auteur(s)	: H.Touchette
! Revision	: 4 aout 1997 (Touchette)
! Notes		: - UNIT 3
!--------------------------------------------------------------------------------------
	subroutine lecture_cova(n,covm)
	use parametre,only: ficcov,kx,ky,c,b
	implicit none
	integer					:: i,j,k1,k2,l1,l2,p,lim
	integer,intent(in)			:: n
	real(8),dimension (n,n),intent(out)	:: covm
	character				:: buffer
!
	open(unit=3,file=ficcov)
	rewind(unit=3)
	do i=1,5
	   read(3,*) buffer
	end do
	lim=(kx-1)*b
	do p=1,lim
	  do i=1,n
	    do j=1,n
	      read(3,*) buffer
	    end do
	  end do
	end do
	do p=1,ky-1
	  do i=1,n
	    do j=1,n
	      read(3,*) buffer
	    end do
	  end do
	end do
	do i=1,n
	   do j=1,n
	      read(3,*)  k1,k2,l1,l2,covm(i,j)
	   end do
	end do
	end subroutine lecture_cova
!--------------------------------------------------------------------------------------
! Nom		: sortie		[sub]					4
! Fonction	: Impression des resultats
! Auteur(s)	: H.Touchette
! Revision	: 4 aout 1997 (Touchette)
!		  13 mars 1998 D. Poulin
! Notes		:
!--------------------------------------------------------------------------------------
	subroutine sortie(w,am,al,pr,ss0,ss1,ss2)
	use parametre,only: m,niterp
	implicit none
	integer					:: i
	real(8),dimension (m),intent(in)	:: w,am
	real(8),dimension (niterp),intent(in)	:: al,pr
	real(8),dimension(niterp)		:: ss0,ss1,ss2
!
	do i=1,m
	   write(7,34) w(i),am(i)
	end do
	do i=1,niterp
	   write(9,35) al(i),pr(i),ss0(i),ss1(i),ss2(i)
	end do
	
34	format(2(f10.5,2x))
35	format(5(f10.5,2x))
	end subroutine sortie
!--------------------------------------------------------------------------------------
! Nom		: funcd 		[sub]					6
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	subroutine funcd(m,d,alpha,entr,falpha,df)
	implicit real*8(a-h,o-z)
	dimension d(m)
	falpha=2.d0*alpha*entr
	df=2.d0*entr
	do i=1,m
	   falpha=falpha-d(i)/(alpha+d(i))
	   df=df+d(i)/(alpha+d(i))**2
	end do
	end subroutine funcd
!--------------------------------------------------------------------------------------
! Nom		: TQLI			[sub]					8
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE TQLI(D,E,N,NP,Z)
      implicit real*8(a-h,o-z)
      DIMENSION D(NP),E(NP),Z(NP,NP)
      DIMENSION WKSP(Np),IWKSP(Np)
      IF (N.GT.1) THEN
      DO 11 I=2,N
      E(I-1)=E(I)
11    CONTINUE
      E(N)=0.d0
      DO 15 L=1,N
      ITER=0
1     DO 12 M=L,N-1
      DD=dABS(D(M))+dABS(D(M+1))
      IF (dABS(E(M))+DD.EQ.DD) GO TO 2
12    CONTINUE
      M=N
2     IF(M.NE.L)THEN
      IF(ITER.EQ.1000)PAUSE 'too many iterations'
      ITER=ITER+1
      G=(D(L+1)-D(L))/(2.d0*E(L))
      R=dSQRT(G**2+1.d0)
      G=D(M)-D(L)+E(L)/(G+dSIGN(R,G))
      S=1.d0
      C=1.d0
      P=0.d0
      DO 14 I=M-1,L,-1
      F=S*E(I)
      B=C*E(I)
      IF(dABS(F).GE.dABS(G))THEN
      C=G/F
      R=dSQRT(C**2+1.d0)
      E(I+1)=F*R
      S=1.d0/R
      C=C*S
      ELSE
      S=F/G
      R=dSQRT(S**2+1.d0)
      E(I+1)=G*R
      C=1.d0/R
      S=S*C
      ENDIF
      G=D(I+1)-P
      R=(D(I)-G)*S+2.d0*C*B
      P=S*R
      D(I+1)=G+P
      G=C*R-B
      DO 13 K=1,N
      F=Z(K,I+1)
      Z(K,I+1)=S*Z(K,I)+C*F
      Z(K,I)=C*Z(K,I)-S*F
13    CONTINUE
14    CONTINUE
      D(L)=D(L)-P
      E(L)=G
      E(M)=0.d0
      GO TO 1
      ENDIF
15    CONTINUE
      ENDIF
      call SORT3(N,np,D,Z,WKSP,IWKSP)
      RETURN
      END
!--------------------------------------------------------------------------------------
! Nom		: TRED			[sub]					9
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE TRED(A,N,NP,D,E)
      implicit real*8(a-h,o-z)
      DIMENSION A(NP,NP),D(NP),E(NP)
      IF(N.GT.1)THEN
      DO 18 I=N,2,-1
      L=I-1
      H=0.d0
      SCALE=0.d0
      IF(L.GT.1)THEN
      DO 11 K=1,L
      SCALE=SCALE+dABS(A(I,K))
11    CONTINUE
      IF(SCALE.EQ.0.d0)THEN
      E(I)=A(I,L)
      ELSE
      DO 12 K=1,L
      A(I,K)=A(I,K)/SCALE
      H=H+A(I,K)**2
12    CONTINUE
      F=A(I,L)
      G=-dSIGN(dSQRT(H),F)
      E(I)=SCALE*G
      H=H-F*G
      A(I,L)=F-G
      F=0.d0
      DO 15 J=1,L
      A(J,I)=A(I,J)/H
      G=0.d0
      DO 13 K=1,J
      G=G+A(J,K)*A(I,K)
13    CONTINUE
      IF(L.GT.J)THEN
      DO 14 K=J+1,L
      G=G+A(K,J)*A(I,K)
14    CONTINUE
      ENDIF
      E(J)=G/H
      F=F+E(J)*A(I,J)
15    CONTINUE
      HH=F/(H+H)
      DO 17 J=1,L
      F=A(I,J)
      G=E(J)-HH*F
      E(J)=G
      DO 16 K=1,J
      A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16    CONTINUE
17    CONTINUE
      ENDIF
      ELSE
      E(I)=A(I,L)
      ENDIF
      D(I)=H
18    CONTINUE
      ENDIF
      D(1)=0.d0
      E(1)=0.d0
      DO 23 I=1,N
      L=I-1
      IF(D(I).NE.0.d0)THEN
      DO 21 J=1,L
      G=0.d0
      DO 19 K=1,L
      G=G+A(I,K)*A(K,J)
19    CONTINUE
      DO 20 K=1,L
      A(K,J)=A(K,J)-G*A(K,I)
20    CONTINUE
21    CONTINUE
      ENDIF
      D(I)=A(I,I)
      A(I,I)=1.d0
      IF(L.GE.1)THEN
      DO 22 J=1,L
      A(I,J)=0.d0
      A(J,I)=0.d0
22    CONTINUE
      ENDIF
23    CONTINUE
      RETURN
      END
!--------------------------------------------------------------------------------------
! Nom		: SORT3			[sub]					10
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE SORT3(N,np,RA,RB,WKSP,IWKSP)
      implicit real*8(a-h,o-z)
      DIMENSION RA(Np),RB(Np,Np),WKSP(Np),IWKSP(Np)
      CALL INDEXX(N,np,RA,IWKSP)
      DO 11 J=1,N
      WKSP(J)=RA(J)
11    CONTINUE
      DO 12 J=1,N
      RA(J)=WKSP(IWKSP(J))
12    CONTINUE
      do i=1,n
      DO 13 J=1,N
      WKSP(J)=RB(i,J)
13    CONTINUE
      DO 14 J=1,N
      RB(i,J)=WKSP(IWKSP(J))
14    CONTINUE
      end do
      RETURN
      END
!--------------------------------------------------------------------------------------
! Nom		: INDEXX		[sub]					11
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE INDEXX(N,np,ARRIN,INDX)
      implicit real*8(a-h,o-z)
      DIMENSION ARRIN(Np),INDX(Np)
      DO 11 J=1,N
      INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
      IF(L.GT.1)THEN
      L=L-1
      INDXT=INDX(L)
      Q=ARRIN(INDXT)
      ELSE
      INDXT=INDX(IR)
      Q=ARRIN(INDXT)
      INDX(IR)=INDX(1)
      IR=IR-1
      IF(IR.EQ.1)THEN
      INDX(1)=INDXT
      RETURN
      ENDIF
      ENDIF
      I=L
      J=L+L
20    IF(J.LE.IR)THEN
      IF(J.LT.IR)THEN
      IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
      ENDIF
      IF(Q.LT.ARRIN(INDX(J)))THEN
      INDX(I)=INDX(J)
      I=J
      J=J+J
      ELSE
      J=IR+1
      ENDIF
      GO TO 20
      ENDIF
      INDX(I)=INDXT
      GO TO 10
      return
      END
!--------------------------------------------------------------------------------------
! Nom		: nnls	 		[sub]					12
! Fonction	: Minimisation positive des moindres carres
!		  (non-negative least-squares)
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      Subroutine nnls(a,mda,m,n,b,x)
      implicit real*8(a-h,o-z)
      dimension a(mda,n),b(m),x(n),w(n),zz(m)
      dimension index(n),dzt(n)
      zero=0.
      one=1.
      two=2.
      factor=0.01
!
      mode=1
      if(m.gt.0.and.n.gt.0) go to 10
      mode=2
      return
10    iter=0
      itmax=3*n
      do 20 i=1,n
      x(i)=zero
20    index(i)=i
      iz2=n
      iz1=1
      nstep=0
      npp1=1
30    continue
      if(iz1.gt.iz2.or.nstep.ge.m) go to 350
      do 50 iz=iz1,iz2
      j=index(iz)
      sm=zero
      do 40 l=npp1,m
40    sm=sm+a(l,j)*b(l)
50    w(j)=sm
60    wmax=zero
      do 70 iz=iz1,iz2
      j=index(iz)
      if(w(j).le.wmax) go to 70
      wmax=w(j)
      izmax=iz
70    continue
      if(wmax) 350,350,80
80    iz=izmax
      j=index(iz)
      asave=a(npp1,j)
      call h12(1,npp1,npp1+1,m,a(1,j),1,up,dzt,1,1,0)
!     call h12(1,npp1,npp1+1,m,a(1,j),1,up,dummy,1,1,0)
      unorm=zero
      if(nstep.eq.0) go to 100
      do 90 l=1,nstep
90    unorm=unorm+a(l,j)**2
100   unorm=dsqrt(unorm)
      if(diff(unorm+dabs(a(npp1,j))*factor,unorm))130,130,110
110   do 120 l=1,m
120   zz(l)=b(l)
      call h12(2,npp1,npp1+1,m,a(1,j),1,up,zz,1,1,1)
      ztest=zz(npp1)/a(npp1,j)
      if(ztest)130,130,140
130   a(npp1,j)=asave
      w(j)=zero
      go to 60
140   do 150 l=1,m
150   b(l)=zz(l)
      index(iz)=index(iz1)
      index(iz1)=j
      iz1=iz1+1
      nstep=npp1
      npp1=npp1+1
      if(iz1.gt.iz2) go to 170
      do 160 jz=iz1,iz2
      jj=index(jz)
160   call h12(2,nstep,npp1,m,a(1,j),1,up,a(1,jj),1,mda,1)
170   continue
      if(nstep.eq.m) go to 190
      do 180 l=npp1,m
180   a(l,j)=zero
190   continue
      w(j)=zero
      assign 200 to next
      go to 400
200   continue
210   iter=iter+1
      if(iter.le.itmax) go to 220
      mode=3
      go to 350
220   continue
      alpha=two
      do 240 ip=1,nstep
      l=index(ip)
      if(zz(ip))230,230,240
230   t=-x(l)/(zz(ip)-x(l))
      if(alpha.le.t) go to 240
      alpha=t
      jj=ip
240   continue
      if(alpha.eq.two) go to 330
      do 250 ip=1,nstep
      l=index(ip)
250   x(l)=x(l)+alpha*(zz(ip)-x(l))
      i=index(jj)
260   x(i)=zero
      if(jj.eq.nstep) go to 290
      jj=jj+1
      do 280 j=jj,nstep
      ii=index(j)
      index(j-1)=ii
      call g1(a(j-1,ii),a(j,ii),cc,ss,a(j-1,ii))
      a(j,ii)=zero
      do 270 l=1,n
      if(l.ne.ii) call g2(cc,ss,a(j-1,l),a(j,l))
270   continue
280   call g2(cc,ss,b(j-1),b(j))
290   npp1=nstep
      nstep=nstep-1
      iz1=iz1-1
      index(iz1)=i
      do 300 jj=1,nstep
      i=index(jj)
      if(x(i)) 260,260,300
300   continue
      do 310 i=1,m
310   zz(i)=b(i)
      assign 320 to next
      go to 400
320   continue
      go to 210
330   do 340 ip=1,nstep
      i=index(ip)
340   x(i)=zz(ip)
      go to 30
350   sm=zero
      if(npp1.gt.m) go to 370
      do 360 i=npp1,m
360   sm=sm+b(i)**2
      go to 390
370   do 380 j=1,n
380   w(j)=zero
390   rnorm=dsqrt(sm)
      return
400   do 430 l=1,nstep
      ip=nstep+1-l
      if(l.eq.1) go to 420
      do 410 ii=1,ip
410   zz(ii)=zz(ii)-a(ii,jj)*zz(ip+1)
420   jj=index(ip)
430   zz(ip)=zz(ip)/a(ip,jj)
      go to next, (200,320)
      end
!--------------------------------------------------------------------------------------
! Nom		: h12			[sub]					13
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      subroutine h12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
      implicit real*8(a-h,o-z)
      dimension u(1,m),c(m)
      one=1.
      if(0.ge.lpivot.or.lpivot.ge.l1.or.l1.gt.m) return
      cl=dabs(u(1,lpivot))
      if(mode.eq.2) go to 60
      do 10 j=l1,m
10    cl=dmax1(dabs(u(1,j)),cl)
      if(cl) 130,130,20
20    clinv=one/cl
      sm=(u(1,lpivot)*clinv)**2
      do 30 j=l1,m
30    sm=sm+(u(1,j)*clinv)**2
      sm1=sm
      cl=cl*dsqrt(sm1)
      if(u(1,lpivot)) 50,50,40
40    cl=-cl
50    up=u(1,lpivot)-cl
      u(1,lpivot)=cl
      go to 70
60    if(cl) 130,130,70
70    if(ncv.le.0) return
      b=up*u(1,lpivot)
      if(b) 80,130,130
80    b=one/b
      i2=1-icv+ice*(lpivot-1)
      incr=ice*(l1-lpivot)
      do 120 j=1,ncv
      i2=i2+icv
      i3=i2+incr
      i4=i3
      sm=c(i2)*up
      do 90 i=l1,m
      sm=sm+c(i3)*u(1,i)
90    i3=i3+ice
      if(sm) 100,120,100
100   sm=sm*b
      c(i2)=c(i2)+sm*up
      do 110 i=l1,m
      c(i4)=c(i4)+sm*u(1,i)
110   i4=i4+ice
120   continue
130   return
      end
!--------------------------------------------------------------------------------------
! Nom		: g1 			[sub]					14
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	subroutine g1 (a,b,cs,si,sig)
	implicit real*8(a-h,o-z)
	zero=0.
	one=1.
	if(dabs(a).le.dabs(b)) go to 10
	xr=b/a
	yr=dsqrt(one+xr**2)
	cs=dsign(one/yr,a)
	si=cs*xr
	sig=dabs(a)*yr
	return
10	if(b) 20,30,20
20	xr=a/b
	yr=dsqrt(one+xr**2)
	si=dsign(one/yr,b)
	cs=si*xr
	sig=dabs(b)*yr
	return
30	sig=zero
	cs=zero
	si=one
	return
	end subroutine g1
!--------------------------------------------------------------------------------------
! Nom		: g2 			[sub]					15
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	subroutine g2(cs,si,x,y)
	implicit real*8(a-h,o-z)
	xr=cs*x+si*y
	y=-si*x+cs*y
	x=xr
	end subroutine g2
!--------------------------------------------------------------------------------------
! Nom		: diff 			[function]				16
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	function diff(x,y)
	real*8 diff,x,y
	diff=x-y
 	end function diff
!--------------------------------------------------------------------------------------
! Nom		: delta 		[function]				17
! Fonction	: Delta de Kronecker
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	function delta(i,j)
	real*8 delta
	integer i,j
	delta=0.d0
	if(i.eq.j) delta=1.d0
	end function delta
!--------------------------------------------------------------------------------------
! Nom		: INVMG 	[SUB]						20-(1)
! Fonction	: Inversion d'une matrice carree
! Auteurs	: H.Touchette
! Revision	: 31 juillet 97 (Touchette)
! Notes		: A -> A^(-1)=AI
!		  Provient de Numerical Recipes, p.40.
! Dependances	: LUDCMP,LUBKSB (externes)
!--------------------------------------------------------------------------------------
	SUBROUTINE INVMG(NP,NM,A)
	IMPLICIT NONE
	REAL(8)					:: D		! Var. locale
	INTEGER					:: I,J		! var. locale
	INTEGER,INTENT(IN)			:: NP,NM
	INTEGER,DIMENSION (NP)			:: INDX   	! Var. locale
	REAL(8),DIMENSION (NP,NP)		:: Y		! Var. locale
	REAL(8),DIMENSION (NP,NP),INTENT(INOUT)	:: A
!
	DO I=1,NM
	   DO J=1,NM
	      Y(I,J)=0.D0
	   END DO
	   Y(I,I)=1.D0
	END DO
	CALL LUDCMP(A,NM,NP,INDX,D)
	DO J=1,NM
	   CALL LUBKSB(A,NM,NP,INDX,Y(1,J))
	END DO
	A=Y
	END SUBROUTINE INVMG
!--------------------------------------------------------------------------------------
! Nom		: LUDCMP 	[SUB]						19
! Fonction	: Decomposition gaussienne A = LU
! Auteurs	:
! Revision	:
! Notes		: Provient entierement de Numerical Recipes
!--------------------------------------------------------------------------------------
      	SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      	IMPLICIT REAL(8) (A-H,O-Z)
      	PARAMETER (NMAX=200,TINY=1.0E-20)
      	DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      	D=1.D0
      	DO 12 I=1,N
           AAMAX=0.D0
           DO 11 J=1,N
              IF(ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11         CONTINUE
           IF(AAMAX.EQ.0.) PAUSE '--- Erreur : Singular matrix LUDCMP'
           VV(I)=1.D0/AAMAX
12     	CONTINUE
      	DO 19 J=1,N
           IF(J.GT.1) THEN
              DO 14 I=1,J-1
                 SUM=A(I,J)
                 IF(I.GT.1) THEN
                    DO 13 K=1,I-1
                       SUM=SUM-A(I,K)*A(K,J)
13                  CONTINUE
                    A(I,J)=SUM
                 ENDIF
14            CONTINUE
           ENDIF
           AAMAX=0.D0
           DO 16 I=J,N
              SUM=A(I,J)
              IF(J.GT.1) THEN
                 DO 15 K=1,J-1
                    SUM=SUM-A(I,K)*A(K,J)
15               CONTINUE
                 A(I,J)=SUM
              ENDIF
              DUM=VV(I)*ABS(SUM)
              IF(DUM.GE.AAMAX) THEN
                 IMAX=I
                 AAMAX=DUM
              ENDIF
16         CONTINUE
           IF(J.NE.IMAX) THEN
              DO 17 K=1,N
                 DUM=A(IMAX,K)
                 A(IMAX,K)=A(J,K)
                 A(J,K)=DUM
17            CONTINUE
              D=-D
              VV(IMAX)=VV(J)
           ENDIF
           INDX(J)=IMAX
           IF(J.NE.N) THEN
             IF(A(J,J).EQ.0.) A(J,J)=TINY
             DUM=1./A(J,J)
             DO 18 I=J+1,N
                A(I,J)=A(I,J)*DUM
18           CONTINUE
           ENDIF
19    	CONTINUE
      	IF(A(N,N).EQ.0.) A(N,N)=TINY
      	RETURN
      	END
!--------------------------------------------------------------------------------------
! Nom		: LUBKSB 	[SUB]						20-(2)
! Fonction	: Backsubstitution (pour l'inversion)
! Auteurs	:
! Revision	:
! Notes		: Provient de Numerical Recipes
!--------------------------------------------------------------------------------------
	SUBROUTINE LUBKSB(A,N,NP,INDX,B)
	INTEGER N,NP,INDX(N)
	REAL(8) A(NP,NP),B(N)
	INTEGER I,II,J,LL
	REAL(8) SUM
	II=0
	DO I=1,N
	   LL=INDX(I)
	   SUM=B(LL)
	   B(LL)=B(I)
	   IF(II.NE.0) THEN
	      DO J=II,I-1
	         SUM=SUM-A(I,J)*B(J)
	      END DO
	   ELSE IF(SUM.NE.0) THEN
	      II=I
	   END IF
	   B(I)=SUM
	END DO
	DO I=N,1,-1
	   SUM=B(I)
	   DO J=I+1,N
	      SUM=SUM-A(I,J)*B(J)
	   END DO
	   B(I)=SUM/A(I,I)
	END DO
	RETURN
	END
!--------------------------------------------------------------------------------------
! Nom		:NUMTOASCII
! Fonction	:Transforme en CHARACTER
! Auteur	:C. Foisi
!--------------------------------------------------------------------------------------
	CHARACTER(LEN=10) FUNCTION NUMTOASCII(I)
	IMPLICIT NONE
	CHARACTER(LEN=10)		:: TMP
	INTEGER,INTENT(IN)		:: I
	INTEGER				:: J,K
	J=I
	K=1
	NUMTOASCII(1:10)='          '
	TMP=NUMTOASCII
	DO WHILE (J.GT.0.OR.K.EQ.1)
	    TMP(K:K)=CHAR(ICHAR('0')+MOD(J,10))
	    K=K+1
	    J=J/10
	END DO
	J=INDEX(TMP,' ')
	DO K=1,J-1
	    NUMTOASCII(K:K)=TMP(J-K:J-K)
	END DO
	END FUNCTION NUMTOASCII
