!--------------------------------------------------------------------------------------
!				     MODULES EXTERIEURS
!--------------------------------------------------------------------------------------
!     Constantes de dimensionnement a priori (constantes proteges)
	MODULE CONSTANTES	   
	   INTEGER,PARAMETER		:: ND1=10		! Dimension max N 
	   INTEGER,PARAMETER		:: ND2=10		! Dimension max M	
	   INTEGER,PARAMETER		:: NP=ND1*ND2		! N*M
	   INTEGER,PARAMETER		:: NL=51		! TRAN (NL > TRAN+1)
	   INTEGER,PARAMETER		:: NCOR=50		! Nb blocs 
	   INTEGER,PARAMETER		:: NZT=7		! Nb options ZETA
	   INTEGER,PARAMETER		:: NSC=1		! Nb res. scalaire ext
	   INTEGER,PARAMETER		:: NRC=2		! Nb res. scalaire int
	   INTEGER,PARAMETER		:: NTE=5		! Nb res. temps egal
	   INTEGER,PARAMETER		:: NTI=5		! Nb res. temps inegal
	   INTEGER,PARAMETER		:: NTG=3		! Nb res. fct. Green
	   INTEGER,PARAMETER		:: NXC=NSC+NRC  	! Nb res. scalaire tot
	   INTEGER,PARAMETER		:: NCM=NSC+NTE+NTI+NTG	! Nb res. total (YNMES) 
	   INTEGER,PARAMETER		:: NCA=60		! Longueur caracteres
	END MODULE CONSTANTES		
!     ---------------------------------------------------------------------------------	   		     		
!     Fichiers externes (parametres ou variables remanents)
	MODULE FICHIERS
	   USE CONSTANTES,ONLY: NCA
	   CHARACTER(LEN=*),PARAMETER		:: FICREF='tampin.dat'	! Source (||)
	   CHARACTER(LEN=NCA),SAVE		:: FICIN		! Entree
	   CHARACTER(LEN=NCA),SAVE		:: FICSED		! Germes	
	   CHARACTER(LEN=NCA),SAVE		:: FICHSIN		! HS archivee
	   CHARACTER(LEN=NCA),SAVE		:: FICHSOUT		! HS sauvegarde
	   CHARACTER(LEN=NCA),SAVE		:: FICGIOUT		! Sortie Green
	   CHARACTER(LEN=NCA),SAVE		:: FICPAOUT		! Sortie para
	   CHARACTER(LEN=NCA),SAVE		:: FICCVOUT		! Sortie covar
	   CHARACTER(LEN=NCA),SAVE		:: FICWEB		! Sortie web
	   CHARACTER(LEN=NCA),SAVE		:: REPWEB		! Repert web
	END MODULE FICHIERS 
!     ---------------------------------------------------------------------------------
!     Parametres de simulation (variables non-remanentes,declaration du type seulement)	
	MODULE PARA_SIM
	   USE CONSTANTES,ONLY: NCM,NZT,NP
	   TYPE PARAMETRE					
	     SEQUENCE
	     REAL(KIND(1.D0))			:: BETA,U,MU,REMPLD,TPAR,TPER,TPRM
	     INTEGER				:: TRAN,N,M,MES,NWARM,NMODUL
	     INTEGER				:: NR1,NR2,NR3,NR4
	     INTEGER				:: NBLOC,NIMPRI,NOPT
	     INTEGER				:: NVEC
	     INTEGER,DIMENSION(NP,2)		:: VEC
	     CHARACTER(LEN=1),DIMENSION(NCM)	:: YNMES	
	     CHARACTER(LEN=1),DIMENSION(NZT)	:: ZETA
	   END TYPE PARAMETRE
	END MODULE PARA_SIM
!     ---------------------------------------------------------------------------------
!     Forme factorisee de la matrice de hopping (variables remanentes et constantes)
	MODULE HOPPING
	   USE CONSTANTES,ONLY: NP
	   REAL(KIND(1.D0)),DIMENSION (NP)	:: DV		! Val. propres de H 
	   REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: ZZ		! Vect. propres de H
	   REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: H		! exp(K)
	   REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: HM		! exp(-K)
	   REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: H2P		! exp(K/2)
	   REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: H2M		! exp(-K/2)
	   SAVE							! Remanence des valeurs
	END MODULE HOPPING
!     ---------------------------------------------------------------------------------
!     Localisation des chaines et des colonnes (variables remanentes et constantes) 
	MODULE UTILITAIRE
	   USE CONSTANTES,ONLY: NP,NL,ND1,ND2
	   INTEGER,DIMENSION (NP)			:: CHI		! Localisation chaines
	   INTEGER,DIMENSION (NP)			:: COLI		! Localisation colonnes
	   REAL(KIND(1.D0)),DIMENSION (0:NL)		:: W		! Coefficient integrale 
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,NP,NP)	:: COSMAG	! Cosinus precalcules
	   SAVE
	END MODULE UTILITAIRE
!     ---------------------------------------------------------------------------------
!     Identification (variables remanentes et constantes)
	MODULE IDENT
	   USE CONSTANTES,ONLY: NCA
	   CHARACTER(LEN=NCA)			:: NANUM	! Nom fichier sortie
	   CHARACTER(LEN=NCA)			:: PREFIX	! Prefixe sortie 
	   CHARACTER(LEN=NCA)			:: SUFFIX	! Suffixe sortie
	   CHARACTER(LEN=8)			:: DTE		! Date
	   INTEGER				:: TI		! Heure initiale
	   INTEGER				:: IDNUM	! Germe initial
	   SAVE							! Remanence des valeurs
	END MODULE IDENT
!     ---------------------------------------------------------------------------------
!     Module d'interfacage (cf. WRITEOUT)
	MODULE INTER
	   INTERFACE 
	      SUBROUTINE WRITEOUT(PARA,NU,E,VE,F,VF)
	      USE CONSTANTES,ONLY: ND1,ND2,NL
	      USE PARA_SIM
	      INTEGER,INTENT(IN)						:: NU
	      REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(IN),OPTIONAL		:: E,VE
	      REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL),INTENT(IN),OPTIONAL	:: F,VF
	      TYPE(PARAMETRE),INTENT(IN)					:: PARA
	      END SUBROUTINE WRITEOUT
	   END INTERFACE
	END MODULE INTER
!--------------------------------------------------------------------------------------
!				PROGRAMME PRINCIPAL
!--------------------------------------------------------------------------------------
	PROGRAM HUBB
!--------------------------------------------------------------------------------------
! Nom	     :	HUBB	[MAIN]							0
!	
! Version    : 	4_2  (10 aout 1999) (version anterieure: hubb4_0)
!
! Fonction   :	Simulation Monte Carlo du modele de Hubbard basee sur la
!		methode du determinant (BSS)
!
! Auteur     : 	H.Touchette selon la version trnb de A.-M.Tremblay
!
! Adresse    : 	Centre de recherche en physique du solide (CRPS)
!		Departement de physique, Universite de Sherbrooke
!		Sherbrooke, Quebec, Canada, J1K 2R1
!		hugo.touchette@physique.usherb.ca  (H. Touchette)
!		tremblay@physique.usherb.ca	   (A.-M. Tremblay)
!
! References :  [1] Programme trnb3, prog10#, et references citees	
!		[2] Notes personnelles de Hugo Touchette, dept. de physique,
!		    Universite de Sherbrooke, 1996
!		[3] H.Touchette, Rapports techniques des projets MCQ, CRPS, 1996-1997  
!		    http://www.physique.usherb.ca/~tremblay/hugo/
!		[4] D. Poulin et H. Touchette, Aspects numériques des simulations Monte
!		    Carlo quantiques du modele de Hubbard, Rapport technique du CRPS, 1999.
!		[5] A.-M. Tremblay, Méthode Monte Carlo pour électrons sur réseau, 1993.
!
! Notes: 	
! - Version Fortran90 de trnb ou prog10 (compilation xlf90)
! - Traitement du cas U > 0 et U < 0
! - Tient compte de la symetrie GUP-GDWN pour U < 0
! - Moyennage optionnel des fcts de Green a U > 0 (MESURE) 
! - Nouvelles options de calcul ZETA(1 a 6)
! - Stabilisation du calcul des fonctions de Green a temps egal et inegal
! - Nouveau generateur de nombres pseudo-aleatoires GGUBFS (Numerical Recipes)
! - Remplacement des COMMON par des MODULE 	
! - Nouveau MODULE d'accumulation pour le calcul des statistiques
! - Changements pour la version parallele notes (||)
! - Nouveau format de sortie et automatisation de l'archivage sur le web
! - Sortie optionnelle des configuration HS. Forme: NM colonnes par TRAN lignes
! - Calcul optionnel de la matrice de covariance (COVARIANCE)
! - Calcul optionnel de la self-energie (SELF)
! 		
!--------------------------------------------------------------------------------------
	USE CONSTANTES,ONLY: NP,NL,NCA
	USE FICHIERS
	USE PARA_SIM
	USE IDENT
	IMPLICIT NONE
!     Variables scalaires
	CHARACTER(LEN=7)		:: FORMATHS
	CHARACTER(LEN=3)		:: TP
	REAL(KIND(1.D0))		:: SIG
	REAL(KIND(1.D0))		:: DTAU,CTE2
	REAL(KIND(1.D0))		:: VU,VD,VARUP,VARDWN
	REAL(KIND(1.D0))		:: RUP,RDWN,R,P,ALEA
	REAL(KIND(1.D0))		:: DETU,DETD,DETUN,DETDN
        INTEGER                         :: LENFIC
	INTEGER				:: NM,MESC,NCHN,NHS,I,J
	INTEGER				:: NCG,NST
	INTEGER				:: NSWEEP,L,NS		! Boucles principales
	INTEGER(SELECTED_INT_KIND(9))	:: DSEED,DSEEDI		! Germes	
!     Fonctions	
	CHARACTER(3)			:: NUMTOASCII		! conversion(int->c)
	REAL(KIND(1.D0)) 		:: GGUBFS		! Gen nb aleatoires
!     Vecteurs, matrices	
	REAL(KIND(1.D0)),DIMENSION (2)		:: COMPT		! Compteurs	
	REAL(KIND(1.D0)),DIMENSION (NP,NL)	:: VUP,VDWN		! Interaction 	
	REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: GUP,GDWN		! Fonctions de Green
	REAL(KIND(1.D0)),DIMENSION (NP,NP)	:: GUPT,GDWNT		! Fonctions de Green
	INTEGER,DIMENSION (NP,NL)		:: SPIN			! Spins ±1 champ HS
!     Type derive
	TYPE(PARAMETRE)				:: PARA			! Parametres simulation
!     ---------------------------------------------------------------------------------
	CALL SYSTEM_CLOCK(TI)						! Temps depart
	CALL DATE_AND_TIME(DTE)						! Date depart
!
!	Lecture noms de fichiers. Pour la version parallele UNIT=13 -> UNIT='*' (||)
	OPEN(UNIT=13,FILE=FICREF,STATUS='OLD',POSITION='REWIND',ERR=111)
	READ(13,*) SUFFIX
	READ(13,*) PREFIX
	READ(13,*) FICSED
	READ(13,*) FICIN
	READ(13,*) FICHSIN
	READ(13,*) FICWEB
	READ(13,*) REPWEB
	CLOSE(UNIT=13)
!
!	Initialisation du germe
	LENFIC=LEN_TRIM(FICSED)
	OPEN(UNIT=4,FILE=FICSED(1:LENFIC),STATUS='OLD',ERR=111)
	REWIND(UNIT=4)
	READ(4,*,ERR=112) DSEED
	DSEEDI=-1*ABS(DSEED)	
	CLOSE(UNIT=4)
!     ----------------------------BLOC D'INITIALISATION--------------------------------
	CALL ENTREE(PARA)					! Parametres (module)
	NM=PARA%N*PARA%M					! Nb de sites reels
	TP=NUMTOASCII(NM)
	I=LEN_TRIM(TP)
	FORMATHS(1:4+I)='('//TP(1:I)//'I4)'
	CALL IDENTIFICATION(PARA,DSEEDI)			! Id fichiers
	CALL INIT_LOC(PARA%N,PARA%M)				! Vec. localisation	
	CALL DEFCOSMAG(PARA%N,PARA%M)				! Cosinus pre-calcules
	CALL INITW(PARA%TRAN)					! Coef. d'integration
	CALL INITH3(NM,PARA)					! H dense
	CALL INITV(PARA,DSEED,VU,VD,SPIN,VUP,VDWN)		! Matrices interaction 		
!
!	Ouverture des fichiers d'ecriture 
	OPEN(UNIT=12,FILE=FICPAOUT,DELIM='APOSTROPHE')
	IF(PARA%ZETA(3).EQ.'Y'.OR.PARA%ZETA(3).EQ.'y')THEN
	    OPEN(UNIT=10,FILE=FICHSOUT)
	    REWIND(UNIT=10)
	ENDIF
	OPEN(UNIT=11,FILE=FICGIOUT,DELIM='APOSTROPHE')
	REWIND(UNIT=12)
	REWIND(UNIT=11)
!     ---------------------------------SIMULATION--------------------------------------
!	Initialisation des constantes de depart
	SIG=1.D0						! Signe initial
	MESC=0							! Mesures effectuees
	COMPT(1:2)=0.D0						! Flips acceptes
	DTAU=PARA%BETA/DBLE(PARA%TRAN)
	CTE2=2*DLOG(VU)						! U < 0	(pour R)
!
	! Boucles principales
  	DO 100 NSWEEP=1,PARA%MES+PARA%NWARM			! Mesures+rechauf
	   DO 120 L=1,PARA%TRAN					! Tranches temps
!	         
	      ! Calcul des fcts de Green a temps egal pour la tranche L
	      ! Calcul complet a tous les NMODUL
	      NCHN=MOD(L,PARA%NMODUL)
	      IF(NCHN.EQ.1.OR.L.EQ.1) THEN			! Calcul complet
	         CALL GCALS(NM,L,PARA,VUP,GUP,DETU)
	         SELECT CASE(PARA%U>0)				! Symetrie GUP,GDWN
	            CASE(.TRUE.)
	               CALL GCALS(NM,L,PARA,VDWN,GDWN,DETD)
	            CASE(.FALSE.)
	               CALL MAC10(NP,NM,GUP,GDWN)
	               DETD=DETU
	         END SELECT
	         IF(NSWEEP.GT.PARA%NWARM.AND.L.EQ.1) THEN
		    MESC=MESC+1
	            CALL MESURE(GUP,GDWN,VUP,VDWN,MESC,COMPT,SIG,PARA)
	            IF(MESC.EQ.PARA%MES) GOTO 101				! Fin
	         ENDIF
	      ELSE								! Wrap
	         CALL MAC4251(NM,2,GUP,VUP(:,L-1))
	         SELECT CASE(PARA%U>0)				! Symetrie GUP,GDWN
	            CASE(.TRUE.);  CALL MAC4251(NM,2,GDWN,VDWN(:,L-1)) 
	            CASE(.FALSE.); CALL MAC10(NP,NM,GUP,GDWN)
	         END SELECT	
	      ENDIF
!
		 ! Mouvements globaux 
		 IF((PARA%ZETA(4).EQ.'Y'.OR.PARA%ZETA(4).EQ.'y').AND.L.EQ.1) THEN
	   	    NCG=INT(GGUBFS(DSEED)*NM)+1	
	   	    NST=SUM(SPIN(NCG,1:PARA%TRAN))		
		    CALL GLOBAL_MOVE(PARA%TRAN,VU,VD,PARA%U,SPIN,VUP,VDWN,NCG)
		    CALL GCALS(NM,L,PARA,VUP,GUPT,DETUN)
		    SELECT CASE(PARA%U>0)
		       CASE(.TRUE.)
		          CALL GCALS(NM,L,PARA,VDWN,GDWNT,DETDN)
		       CASE(.FALSE.)
		          CALL MAC10(NP,NM,GUPT,GDWNT)
		          DETDN=DETUN
		    END SELECT
		    RUP=DETUN/DETU 
		    RDWN=DETDN/DETD 
                    SELECT CASE(PARA%U>0)
                       CASE(.TRUE.);  R=RUP*RDWN				! U > 0
                       CASE(.FALSE.); R=RUP*RDWN*DEXP(CTE2*DBLE(NST))		! U < 0 
                    END SELECT
		    P=DABS(R)/(DABS(R)+1.D0)				
		    ALEA=GGUBFS(DSEED)
		    IF(ALEA.LE.P) THEN
		       COMPT(2)=COMPT(2)+1.D0
		       SIG=SIG*DSIGN(1.D0,R)				
		       CALL MAC10(NP,NM,GUPT,GUP)
		       CALL MAC10(NP,NM,GDWNT,GDWN)
		    ELSE
		       CALL GLOBAL_MOVE(PARA%TRAN,VU,VD,PARA%U,SPIN,VUP,VDWN,NCG)
		    ENDIF
		 ENDIF
!
	      ! Balayage des sites d'une tranche de temps
	      DO 130 NS=1,NM
!
		 ! Definition des constantes selon U et la valeur du champ HS	            
	         IF(PARA%U.GT.0) THEN
                    IF(SPIN(NS,L).EQ.1) THEN			
		       VARUP=VD/VU-1.D0				! SPIN = 1,U > 0
                       VARDWN=VU/VD-1.D0
		    ELSE
		       VARUP=VU/VD-1.D0				! SPIN =-1,U > 0
                       VARDWN=VD/VU-1.D0
		    ENDIF
                 ELSE
                    IF(SPIN(NS,L).EQ.1) THEN
                       VARUP=VD/VU-1.D0				! SPIN = 1,U < 0
                       VARDWN=VD/VU-1.D0
                    ELSE
                       VARUP=VU/VD-1.D0				! SPIN =-1,U < 0
                       VARDWN=VU/VD-1.D0
                    ENDIF
                 ENDIF
!
                 ! Calcul des determinants pour la probabilite de flip 
                 RUP=1.D0+(1.D0-GUP(NS,NS))*VARUP
                 RDWN=1.D0+(1.D0-GDWN(NS,NS))*VARDWN
                 IF(PARA%U.GT.0) THEN
        	    R=RUP*RDWN					! U > 0
		 ELSE
		    R=RUP*RDWN*DEXP(CTE2*DBLE(SPIN(NS,L)))	! U < 0
		 END IF
		 P=DABS(R)/(DABS(R)+1.D0)			! Probabilite
!
		 ! Changement de la valeur du champ HS selon la valeur de ALEA
		 ALEA=GGUBFS(DSEED)
		 IF(ALEA.LE.P) THEN				
	            COMPT(1)=COMPT(1)+1.D0			! Nb spins flippes
	            SIG=SIG*DSIGN(1.D0,R)			! Signe du det.
	            IF(PARA%U.GT.0) THEN
                       IF(SPIN(NS,L).EQ.1) THEN
		          VUP(NS,L)=VD				! SPIN = 1,U > 0
                          VDWN(NS,L)=VU
                          SPIN(NS,L)=-1
		       ELSE
		          VUP(NS,L)=VU				! SPIN =-1,U > 0
                          VDWN(NS,L)=VD
                          SPIN(NS,L)=1
                       ENDIF
                    ELSE
                       IF(SPIN(NS,L).EQ.1) THEN
	 	          VUP(NS,L)=VD				! SPIN = 1,U < 0
                          VDWN(NS,L)=VD
                          SPIN(NS,L)=-1
		       ELSE
		          VUP(NS,L)=VU				! SPIN =-1,U < 0
                          VDWN(NS,L)=VU
                          SPIN(NS,L)=1
		       ENDIF
                    ENDIF
                    CALL UPDTG(NP,NM,NS,GUP,VARUP)		! Mise a jour de G
                    CALL UPDTG(NP,NM,NS,GDWN,VARDWN)
                 ENDIF		    
130           END DO					
120	   END DO 
!
	   ! Sauvegarde de la configuration HS a tous les NIMPRI
	   NHS=MOD(NSWEEP,PARA%NIMPRI)
	   IF(NHS.EQ.0.AND.(PARA%ZETA(3).EQ.'Y'.OR.PARA%ZETA(3).EQ.'y')) THEN
	      WRITE(10,FMT='(I8,A)') DSEED,' Germe'
	      WRITE(10,FMT='(I8,A)') NSWEEP+1,' Mesures effectuees'
	      WRITE(10,FMT=FORMATHS)((SPIN(J,I),J=1,NM),I=1,PARA%TRAN)
	   ENDIF
100	END DO
!     -----------------------------FIN DE LA SIMULATION--------------------------------
101	CONTINUE
!	Remplacement du germe
	OPEN(UNIT=4,FILE=FICSED(1:LENFIC),STATUS='OLD',ERR=113)
	REWIND(UNIT=4)
	WRITE(4,*) DSEED
	CLOSE(UNIT=4)
!
!	Fermeture des fichiers 
	CLOSE(UNIT=12)
	IF(PARA%ZETA(3).EQ.'Y'.OR.PARA%ZETA(3).EQ.'y')THEN	
	    CLOSE(UNIT=10)
	ENDIF
	CLOSE(UNIT=11)
	STOP
111	WRITE(*,*) '--- Erreur : Fichier inexistant (MAIN) ---'
	STOP
112	WRITE(*,*) '--- Erreur : incapable de lire le SEED ---'
	STOP
113     WRITE(*,*) '--- Erreur : Incapable de ouvrir le fichier SEED ---'
        STOP
!
	END PROGRAM HUBB
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
!				  PROCEDURES EXTERNES
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! Nom		: ACCUMULATION_STATISTIQUE	[MODULE]			1
! Fonction	: Accumulation des resultats des mesures et controle de la sortie
! Auteurs	: H.Touchette
! Revision	: 17 novembre 1996
! References	: [1] H.Touchette, Algorithme MCQ: version Fortran 90, rapport 
!		      technique du CRPS, automne 1996.
! Notes		: - Les variables d'accumulation sont rangees dans un MODULE et doivent 
!		    etre explicitement declarees avec l'attribut SAVE pour qu'elles
!		    demeurent remanentes.
!		  - Pour la signification des etages (6) d'accumulation, cf. MESURE2
!		  - NW3 controle la remise a zera de l'etage 2 a tous les blocs
!
! Dependances	: 		       |- Declarations des variables d'accumulation 
!		  MODULE accumulation -|- ACCUMULATOR [SUB] 
!--------------------------------------------------------------------------------------	            
	MODULE ACCUMULATION_STATISTIQUE
!	   ----------------------------------------------------------------------------          
!          Variables d'accumulation pour les mesures (variables remanentes)
	   USE CONSTANTES,ONLY: ND1,ND2,NL,NTE,NTI,NTG,NXC,NCOR
	   REAL(KIND(1.D0))					:: SIGMOY,SBLOC,C1
	   REAL(KIND(1.D0)),DIMENSION (NCOR)			:: AUTOCORR	! Correlation 
	   REAL(KIND(1.D0)),DIMENSION (6,0:NXC)			:: XM		! Res. courants
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,6,NTE)		:: XSTRU	! Fact. struct
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,6,0:NL,NTI)	:: XSTRUI	! Fact. structI
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,6,0:NL,NTG)	:: XGIUP	! Inegal G
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,6,0:NL,NTG)	:: XGIDN	! Inegal G
	   REAL(KIND(1.D0)),DIMENSION(:,:,:,:),ALLOCATABLE	:: COVGU,COVXU,COVJU
	   REAL(KIND(1.D0)),DIMENSION(:,:,:,:),ALLOCATABLE	:: COVMU,COVXM,COVJM
	   SAVE							! Remanence des valeurs
!	   ----------------------------------------------------------------------------
	CONTAINS						
!	   ----------------------------------------------------------------------------
	   SUBROUTINE ACCUMUL(COMPT,MESC,SIG,PARA,DMU,XSTP,FACST,FACSTI,FGIUP,FGIDN)
	   USE PARA_SIM
	   IMPLICIT NONE
	   INTERFACE
		SUBROUTINE SORTIEC(PARA,MESC,C1,AUTOCOR,COVM,COVX,COVJ)
		USE CONSTANTES,ONLY: NL,ND1,ND2,NCOR
		USE FICHIERS,ONLY: FICCVOUT
		USE IDENT,ONLY: IDNUM
		USE PARA_SIM
		IMPLICIT NONE
		REAL(KIND(1.D0)),INTENT(IN)                              	:: C1
		INTEGER,INTENT(IN)                             			:: MESC
		REAL(KIND(1.D0)),DIMENSION(NCOR),INTENT(IN)              	:: AUTOCOR
		REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL,0:NL),INTENT(IN),OPTIONAL:: COVM,COVX,COVJ
		TYPE(PARAMETRE),INTENT(IN)                              	:: PARA
		END SUBROUTINE SORTIEC
	   END INTERFACE
	   REAL(KIND(1.D0))						:: DERIV,XMUPROP 
	   REAL(KIND(1.D0))						:: FLUCT,ENERG
	   INTEGER							:: NW1,NW2,ICQ,J
	   INTEGER							:: I,NM,NCOMES,NR	
	   REAL(KIND(1.D0)),INTENT(IN)					:: DMU,SIG
	   INTEGER,INTENT(IN)						:: MESC	
	   REAL(KIND(1.D0)),DIMENSION (2)				:: COMPT,XCOMP		
	   REAL(KIND(1.D0)),DIMENSION (0:NXC),INTENT(IN)		:: XSTP		
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,NTE),INTENT(IN)		:: FACST
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,NTI),INTENT(IN)	:: FACSTI
	   REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,NTG),INTENT(IN)	:: FGIUP,FGIDN	
	   TYPE(PARAMETRE),INTENT(IN)					:: PARA
!
	   NM=PARA%N*PARA%M
!
!	   Initialisation des variables a la premiere mesure
	   IF(MESC.EQ.1) THEN
	      SIGMOY=0.D0
	      SBLOC=0.D0 
	      C1=0.D0
	      XM(1:6,0:NXC)=0.D0
	      XSTRU(1:PARA%N,1:PARA%M,1:6,1:NTE)=0.D0
	      XGIUP(1:PARA%N,1:PARA%M,1:6,0:PARA%TRAN,1:NTG)=0.D0
	      XGIDN(1:PARA%N,1:PARA%M,1:6,0:PARA%TRAN,1:NTG)=0.D0
	      XSTRUI(1:PARA%N,1:PARA%M,1:6,0:PARA%TRAN,1:NTI)=0.D0
	      AUTOCORR(1:NCOR)=0.D0
	      IF(PARA%ZETA(5).EQ.'Y'.OR.PARA%ZETA(5).EQ.'y')THEN
		ALLOCATE (COVGU(ND1,ND2,0:NL,NCOR),COVMU(ND1,ND2,0:NL,0:NL))
		COVGU(1:PARA%N,1:PARA%M,0:PARA%TRAN,1:NCOR)=0.D0
	      ENDIF
	      IF(PARA%ZETA(6).EQ.'Y'.OR.PARA%ZETA(6).EQ.'y')THEN
		ALLOCATE (COVXU(ND1,ND2,0:NL,NCOR),COVXM(ND1,ND2,0:NL,0:NL))
		COVXU(1:PARA%N,1:PARA%M,0:PARA%TRAN,1:NCOR)=0.d0
	      ENDIF
              IF(PARA%ZETA(7).EQ.'Y'.OR.PARA%ZETA(7).EQ.'y')THEN
                ALLOCATE (COVJU(ND1,ND2,0:NL,NCOR),COVJM(ND1,ND2,0:NL,0:NL))
                COVJU(1:PARA%N,1:PARA%M,0:PARA%TRAN,1:NCOR)=0.D0
              ENDIF
	   ENDIF
!	
!	   Incrementation du signe moyen et SBLOC
	   SIGMOY=SIGMOY+SIG
	   SBLOC=SBLOC+SIG
!
!	   Initialisation des variables pour les statistiques 
	   NW1=0
	   NW2=0
	   IF(MOD(MESC,PARA%NBLOC).EQ.0) THEN			! Calcul par blocs
	      NW1=1
	      C1=C1+1.D0
	   ENDIF
	   NCOMES=MOD(MESC,PARA%NIMPRI)
	   IF(NCOMES.EQ.0.OR.MESC.EQ.PARA%MES) NW2=1		! Calcul impression
!
!	   Choix des statistiques
!	   ----------------------------------------------------------------------------
!	   Statistiques sur les resultats courants (non-optionnel)		! 1
	   CALL MESURE2(NXC,0,NXC,SIG,SBLOC,SIGMOY,C1,NW1,NW2,0,XSTP,XM)
!
!	   Statistique sur les quantites a temps egal (optionnel)
	   NR=PARA%NR1
           DO J=1,PARA%NR2
              NR=NR+1	
	      IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	         CALL MESURE3(ND1,ND2,1,PARA%N,PARA%M,SIG,SBLOC,SIGMOY,C1,&
		             & NW1,NW2,1,FACST(:,:,J),XSTRU(:,:,:,J))
	      ENDIF
	   END DO
!
!	   Statistique sur les quantites a temps inegal (optionnel)
	   DO J=1,PARA%NR3
              NR=NR+1
	      IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	         DO I=0,PARA%TRAN
	            CALL MESURE3(ND1,ND2,1,PARA%N,PARA%M,SIG,SBLOC,SIGMOY,C1,&
	      	                & NW1,NW2,0,FACSTI(:,:,I,J),XSTRUI(:,:,:,I,J))
	         END DO
	      ENDIF
           END DO
!	
!	   Statistiques sur les fcts de Green et les quantites derivees
           DO J=1,PARA%NR4
	      NR=NR+1
	      IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	         DO I=0,PARA%TRAN
	            CALL MESURE3(ND1,ND2,1,PARA%N,PARA%M,SIG,SBLOC,SIGMOY,C1,&
	      		        & NW1,NW2,0,FGIUP(:,:,I,J),XGIUP(:,:,:,I,J))
	            CALL MESURE3(ND1,ND2,1,PARA%N,PARA%M,SIG,SBLOC,SIGMOY,C1,&
	      	                & NW1,NW2,1,FGIDN(:,:,I,J),XGIDN(:,:,:,I,J))
	         END DO
	      ENDIF
           END DO
!
!	   Controle des blocs
	   IF(MOD(MESC,PARA%NBLOC).EQ.0) THEN
	      AUTOCORR(nint(C1))=XM(2,2)/SBLOC
	      IF(PARA%ZETA(5).EQ.'y'.OR.PARA%ZETA(5).EQ.'Y')THEN
		COVGU(1:PARA%N,1:PARA%M,0:PARA%TRAN,nint(C1))=&
		& XGIUP(1:PARA%N,1:PARA%M,2,0:PARA%TRAN,1)/SBLOC
	      ENDIF
	      IF(PARA%ZETA(6).EQ.'y'.OR.PARA%ZETA(6).EQ.'Y')THEN
		COVXU(1:PARA%N,1:PARA%M,0:PARA%TRAN,nint(C1))=&
		& XSTRUI(1:PARA%N,1:PARA%M,2,0:PARA%TRAN,1)/SBLOC
	      ENDIF
              IF(PARA%ZETA(7).EQ.'y'.OR.PARA%ZETA(7).EQ.'Y')THEN
                COVJU(1:PARA%N,1:PARA%M,0:PARA%TRAN,nint(C1))=&
                & XSTRUI(1:PARA%N,1:PARA%M,2,0:PARA%TRAN,1)/SBLOC
              ENDIF
	      XM(2,0:NXC)=0.D0					! Remise a 0 de l'etage 2
	      XGIUP(1:PARA%N,1:PARA%M,2,0:PARA%TRAN,1:PARA%NR4)=0.D0	! Remise a 0 de l'etage 2
	      XSTRUI(1:PARA%N,1:PARA%M,2,0:PARA%TRAN,1:PARA%NR3)=0.D0
	      SBLOC=0.D0
	   ENDIF
!
!	   Controle de la sortie	
	   NCOMES=MOD(MESC,PARA%NIMPRI)
	   IF(NCOMES.EQ.0.OR.MESC.EQ.PARA%MES) THEN
!
	      ! Calcul mu propose et fluct. nb particules
	      DERIV=-DMU/SIGMOY-PARA%BETA*NM*XM(5,0)**2			
	      XMUPROP=PARA%MU+(PARA%REMPLD-XM(5,0))/DERIV		
	      XCOMP(1)=COMPT(1)/(MESC+PARA%NWARM)/(NM*PARA%TRAN)
	      XCOMP(2)=COMPT(2)/(MESC+PARA%NWARM)		
	      FLUCT=DSQRT(ABS(DERIV/(PARA%BETA*NM)))			
!
	      ! Calcul de l'energie des electrons libres
	      CALL ENE_TOT(PARA,ENERG)
!
	      ! Sortie des resultats
	      CALL SORTIE1(PARA,MESC,XCOMP,SIGMOY,FLUCT,XMUPROP,ENERG,XM(5,:),XM(6,:))
	      CALL SORTIE2(PARA,XSTRU(:,:,5,:),XSTRU(:,:,6,:),&
	      		  &XSTRUI(:,:,5,:,:),XSTRUI(:,:,6,:,:))	
	      CALL SORTIE3(PARA,MESC,XGIUP(:,:,5,:,:),XGIDN(:,:,5,:,:),&
	         	  &XGIUP(:,:,6,:,:),XGIDN(:,:,6,:,:))
!
	      IF(MESC.EQ.PARA%MES) THEN
	         IF(PARA%ZETA(1).EQ.'y'.OR.PARA%ZETA(1).EQ.'Y') THEN
	            CALL AUTO_CORRELATION(C1,AUTOCORR)
	         ENDIF
		 ICQ=0
		 IF(PARA%ZETA(5).EQ.'y'.OR.PARA%ZETA(5).EQ.'Y') THEN
                    ICQ=ICQ+1
	            CALL COVARIANCE(PARA,C1,COVGU,XGIUP(:,:,5,:,1),COVMU)
	         ENDIF
		 IF(PARA%ZETA(6).EQ.'y'.OR.PARA%ZETA(6).EQ.'Y') THEN
                    ICQ=ICQ+2
		    CALL COVARIANCE(PARA,C1,COVXU,XSTRUI(:,:,5,:,1),COVXM)
		 ENDIF
                 IF (PARA%ZETA(7).EQ.'Y'.OR.PARA%ZETA(7).EQ.'y')THEN
                    ICQ=ICQ+4
                    CALL COVARIANCE(PARA,C1,COVJU,XSTRUI(:,:,5,:,4),COVJM)
                 ENDIF
                 SELECT CASE (ICQ)
                 CASE (7)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVM=COVMU,COVX=COVXM,&
                                & COVJ=COVJM)
                 CASE (6)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVX=COVXM,COVJ=COVJM)
                 CASE (5)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVM=COVMU,COVJ=COVJM)
                 CASE (4)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVJ=COVJM)
                 CASE (3)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVM=COVMU,COVX=COVXM)
                 CASE (2)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVX=COVXM)
                 CASE (1)
                    CALL SORTIEC(PARA,MESC,C1,AUTOCORR,COVM=COVMU)
                 CASE (0)
	            IF(PARA%ZETA(1).EQ.'y'.OR.PARA%ZETA(1).EQ.'Y') THEN
		       CALL SORTIEC(PARA,MESC,C1,AUTOCORR)
		    ENDIF
                 END SELECT
	         CALL SORTIE4(PARA,XMUPROP,XM(5,0))
		 IF(PARA%ZETA(5).EQ.'Y'.OR.PARA%ZETA(5).EQ.'y')THEN
		    DEALLOCATE (COVMU,COVGU)
		 ENDIF
                 IF(PARA%ZETA(6).EQ.'Y'.OR.PARA%ZETA(6).EQ.'y')THEN
                    DEALLOCATE (COVXM,COVXU)
                 ENDIF
                 IF(PARA%ZETA(7).EQ.'Y'.OR.PARA%ZETA(7).EQ.'y')THEN
                    DEALLOCATE (COVJM,COVJU)
                 ENDIF
	      ENDIF
	   ENDIF
!	
	   END SUBROUTINE ACCUMUL
!	   ----------------------------------------------------------------------------
	END MODULE ACCUMULATION_STATISTIQUE
!--------------------------------------------------------------------------------------
! Nom		: MESURE 	[SUB]						2
! Fonction	: Appel des routines de mesure 
! Auteurs	: AM.Tremblay,S.Caille,H.Touchette
! Revision	: 20 fevrier 1996 (Touchette)
! Notes		: - L'extension de cette routine pour de nouvelles mesures est decrite
!		    dans les rapports des projets MCQ,H.Touchette,CPRS,Automne 1996
!		  - FGIUP(2,:,:,:)=Re{G(iwn,Qx,Qy)}
!		  - FGIDN(2,:,:,:)=Im{G(iwn,Qx,Qy)}
!		  - FGIUP(3,:,:,:)=Re{self-energie}
!		  - FGIDN(3,:,:,:)=Im{self-energie}
!		  - Le moyennage des fcts de Green up et down a U > 0 est optionnel.
!		    Enlevez les commentaires ci-dessous pour activer cette option.
!--------------------------------------------------------------------------------------	
	SUBROUTINE MESURE(GUP,GDN,VUP,VDN,MESC,COMPT,SIG,PARA)
	USE CONSTANTES,ONLY: NP,ND1,ND2,NL,NTE,NTI,NTG,NXC
	USE PARA_SIM
	USE ACCUMULATION_STATISTIQUE
	IMPLICIT NONE
	INTEGER						:: NM,L,I,NR,ICQ,NCHN,TR ! Var locale
	REAL(KIND(1.D0))				:: DETU,DTAU
	REAL(KIND(1.D0)),SAVE				:: DMU		! Remanent
	REAL(KIND(1.D0)),INTENT(IN)			:: SIG
	INTEGER,INTENT(IN)				:: MESC
	REAL(KIND(1.D0)),DIMENSION (2)			:: COMPT 
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDN
	REAL(KIND(1.D0)),DIMENSION (NP,NL),INTENT(IN)	:: VUP,VDN
	REAL(KIND(1.D0)),DIMENSION (0:NXC)		:: XSTP		! Res. courants
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,NTE)	:: FACST	! Fact. struct
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,NTI)	:: FACSTI	! Fact. structI	
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: GIUP,GIDN	! Green
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: GIPUP,GIPDN
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: GEUP,GEDN
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,NTG)	:: FGIUP,FGIDN ! Green,Self
	TYPE(PARAMETRE),INTENT(IN)			:: PARA	
!
	NM=PARA%N*PARA%M
	IF(MESC.EQ.1) DMU=0.D0				
!
!	Calcul des fonctions de Green a temps inegal
	ICQ=0
        NR=PARA%NR1+PARA%NR2
        DO I=1,3
           NR=NR+1
           IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y')THEN
	       ICQ=2
	   ENDIF
        END DO
        DO I=1,PARA%NR3-3
           NR=NR+1
           IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y')THEN
	       ICQ=3
           ENDIF
        END DO
	NR=NR+1
	IF((PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y').AND.(ICQ.EQ.0))THEN
	   ICQ=1
	ENDIF
	
	DO 555, TR=0,PARA%TRAN
	
	IF (TR.LE.PARA%TRAN/2) THEN
	   L=TR
	ELSE
	   L=PARA%TRAN*3/2-TR+1
	END IF
	
	IF(ICQ.GT.0)THEN
	   CALL GICALS(NM,PARA,GUP,VUP,GIUP,L,1)
	   IF(PARA%U.GT.0.)THEN 				! Green inegal
	      CALL GICALS(NM,PARA,GDN,VDN,GIDN,L,2)
	   ELSE
	      GIDN=GIUP
	   END IF
	   IF(ICQ.GT.1) THEN
	      IF(ICQ.GT.2)THEN				! Green inegal
	      NCHN=MOD(L,PARA%NMODUL)
	      IF(L==0.OR.L==PARA%TRAN)THEN
                  GEUP(1:NM,1:NM)=GUP(1:NM,1:NM)
                  GEDN(1:NM,1:NM)=GDN(1:NM,1:NM)
              ELSE IF(NCHN==0) THEN
                 CALL GCALS(NM,L+1,PARA,VUP,GEUP(:,:),DETU)
                 IF(PARA%U.GT.0.) THEN
                    CALL GCALS(NM,L+1,PARA,VDN,GEDN(:,:),DETU)
                 ELSE
                    GEDN=GEUP
                 END IF
              ELSE IF(L.LE.PARA%TRAN/2) THEN
                 CALL MAC4251(NM,2,GEUP,VUP(:,L))
                 IF(PARA%U.GT.0.) THEN
                    CALL MAC4251(NM,2,GEDN,VDN(:,L))
                 ELSE
                    GEDN=GEUP
                 END IF
              ELSE
                 CALL MAC4251(NM,5,GEUP,VUP(:,L+1))
                 IF(PARA%U.GT.0.) THEN
                    CALL MAC4251(NM,5,GEDN,VDN(:,L+1))
                 ELSE
                    GEDN=GEUP
                 END IF
              END IF
	      ENDIF
              CALL GIPCALS(NM,PARA,GUP,VUP,GIPUP,L,1)
              IF(PARA%U.GT.0.) THEN
                 CALL GIPCALS(NM,PARA,GDN,VDN,GIPDN,L,2)
              ELSE
                 GIPDN=GIPUP
              END IF
	   ENDIF
	ENDIF
!
!	Moyennage des fcts de Green up et down pour U > 0
!	GIUP(1:NM,1:NM,0:PARA%TRAN)=&
!	& (GIUP(1:NM,1:NM,0:PARA%TRAN)+GIDN(1:NM,1:NM,0:PARA%TRAN))/2
!	GIDN(1:NM,1:NM,0:PARA%TRAN)=GIUP(1:NM,1:NM,0:PARA%TRAN)	
!
! 	Mesures : choix et appel des mesures a effectuer		
!	-------------------------------------------------------------------------------
!
	NR=PARA%NR1+PARA%NR2+PARA%NR3+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   CALL FTQ(PARA%N,PARA%M,GIUP,FGIUP(:,:,L,1),PARA%NVEC,PARA%VEC(:,:))	
	   CALL FTQ(PARA%N,PARA%M,GIDN,FGIDN(:,:,L,1),PARA%NVEC,PARA%VEC(:,:))	
	ENDIF
	NR=NR+2
	IF(TR==PARA%TRAN) THEN
	   IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	      CALL MATSUBARA(PARA%TRAN,PARA%BETA,&
	   		 &0,FGIUP(:,:,:,1),FGIUP(:,:,:,2),PARA%NVEC,PARA%VEC(:,:))	
	      CALL MATSUBARA(PARA%TRAN,PARA%BETA,&
	   		 &1,FGIUP(:,:,:,1),FGIDN(:,:,:,2),PARA%NVEC,PARA%VEC(:,:))	
	   ENDIF
	   IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	      CALL SELF(PARA,FGIUP(:,:,:,2),FGIDN(:,:,:,2),FGIUP(:,:,:,3),FGIDN(:,:,:,3))
	   ENDIF
	END IF
!
        NR=PARA%NR1+PARA%NR2+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL SUSCEPI(PARA%N,PARA%M,GIUP,GIDN,GIPUP,GIPDN,1,FACSTI(:,:,L,1)&
                        &,PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL SUSCEPI(PARA%N,PARA%M,GIUP,GIDN,GIPUP,GIPDN,2,FACSTI(:,:,L,2)&
                        &,PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL SUSCEPI(PARA%N,PARA%M,GIUP,GIDN,GIPUP,GIPDN,3,FACSTI(:,:,L,3)&
                        &,PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(10).EQ.'y') THEN
           CALL COURANT(PARA%N,PARA%M,GIUP,GIDN,GIPUP,GIPDN,GEUP,GEDN&
			&,GUP,GDN,FACSTI(:,:,L,4),PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	    CALL SUSMAG(PARA%N,PARA%M,GIUP,GIDN,GIPUP,GIPDN,GEUP,GEDN&
			&,GUP,GDN,FACSTI(:,:,L,5),PARA%NVEC,PARA%VEC(:,:))
	ENDIF
	
555	END DO
!
        NR=PARA%NR1+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL MAGST(PARA%N,PARA%M,GUP(:,:),GDN(:,:),FACST(:,:,1),PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL DENSITE(PARA%N,PARA%M,GUP(:,:),GDN(:,:),FACST(:,:,2),PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL RHOE(PARA%N,PARA%M,GUP(:,:),GDN(:,:),FACST(:,:,3),PARA%NVEC,PARA%VEC(:,:))
        ENDIF
        NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
           CALL PSCOR(PARA%N,PARA%M,GUP(:,:),GDN(:,:),FACST(:,:,4),PARA%NVEC,PARA%VEC(:,:))
        ENDIF
	NR=NR+1
	DTAU=PARA%BETA/PARA%TRAN
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'Y')THEN
	    CALL SUMSUSCEP(PARA%TRAN,DTAU,FACSTI(:,:,:,1),&
			   &FACST(:,:,5),PARA%NVEC,PARA%VEC(:,:))
	ENDIF
!
        IF(PARA%YNMES(1).EQ.'Y'.OR.PARA%YNMES(1).EQ.'y') THEN
           CALL DSITE(PARA,GUP(:,:),GDN(:,:),XSTP(1),XSTP(2))
        ENDIF
!
!	Calcul du nombre moyen d'occupation et mu propose
	CALL OCC_REMP(NP,NM,SIG,GUP(:,:),GDN(:,:),DMU,XSTP(0),PARA%BETA) 
!
!	Exportation des resultats de mesures pour l'accumulation 
	CALL ACCUMUL(COMPT,MESC,SIG,PARA,DMU,XSTP,FACST,FACSTI,FGIUP,FGIDN)
			
!
	END SUBROUTINE MESURE
!--------------------------------------------------------------------------------------
! Nom		: MESURE2 	[SUB]						3-(1)
! Fonction	: Calcul des statistiques sur les mesures (resultats courants)
! Auteurs	: S.Caille,AM.Tremblay,H.Touchette
! Revision	: Decembre 1990
! Reference 	: [1] H.Touchette, Algorithme MCQ: version Fortran 90, rapport 
!		      technique du CRPS, automne 96.
! Notes		: 
! Variables	: XME(i,j) -> i: etage d'accumulation
!	      		      j: indice du vecteur contenant les resultats
!    
! Etages 		Signification
!   1		Accumulation en vue de la moyenne totale
!   2		Accumulation sur un bloc en vue de la variance
!   3		Moyenne des blocs
!   4		Moyenne des blocs au carre	
!   5		Moyenne totale
!   6		Variance et ecart-type
!
! SIG 	  	: valeur du signe
! SBLOC   	: nombre de mesures par bloc
! SIGMOY  	: valeur moyenne du signe
! C1	  	: nombre de blocs
! NW1,NW2,NW3 	: indicateurs de conditions:
! Si NW1=1, Le calcul est fait pour le nombre de blocs; le compteur est mis a zero.
! Si NW2=1, Le calcul est fait pour en prevision de l'impression.
! Si NW3=1, L'etage 2 n'est pas remis a zero dans la routine
!--------------------------------------------------------------------------------------
	SUBROUTINE MESURE2(NP,N1,M1,SIG,SBLOC,SIGMOY,C1,NW1,NW2,NW3,VME,XME)
	IMPLICIT NONE
	INTEGER							:: I		! Var. locale
	REAL(KIND(1.D0)),INTENT(IN)				:: SBLOC,SIGMOY,C1,SIG
	INTEGER,INTENT(IN)					:: NP,N1,M1,NW1,NW2,NW3
	REAL(KIND(1.D0)),DIMENSION (N1:NP),INTENT(IN) 		:: VME
	REAL(KIND(1.D0)),DIMENSION(6,N1:NP),INTENT(INOUT)	:: XME	
!
	DO I=N1,M1
	   XME(1,I)=XME(1,I)+VME(I)*SIG
	   XME(2,I)=XME(2,I)+VME(I)*SIG
	END DO
	IF(NW1.EQ.1) THEN
	   DO I=N1,M1
	      XME(3,I)=XME(3,I)+XME(2,I)/SBLOC
	      XME(4,I)=XME(4,I)+(XME(2,I)/SBLOC)**2
	   END DO
	ENDIF
	IF(NW1.EQ.1.AND.NW3.EQ.1) XME(2,N1:M1)=0.D0
	IF(NW2.EQ.1) THEN
	   DO I=N1,M1
	      XME(5,I)=XME(1,I)/SIGMOY
	      XME(6,I)=XME(4,I)/C1-(XME(3,I)/C1)**2
	      XME(6,I)=DSQRT(XME(6,I)/C1)
	   END DO
	ENDIF
	END SUBROUTINE MESURE2
!--------------------------------------------------------------------------------------
! Nom		: MESURE3 	[SUB]						3-(2)
! Fonction	: Calcul des statistiques sur les facteurs de structures
! Auteurs	: H.Touchette
! Revision	: 7 novembre 1996
! Notes		: Voir commentaires de MESURE2 
!--------------------------------------------------------------------------------------    
	SUBROUTINE MESURE3(NP1,NP2,N1,N,M,SIG,SBLOC,SIGMOY,C1,NW1,NW2,NW3,VME,XME)
	IMPLICIT NONE
	INTEGER								:: I,J	! Var. locales 
	REAL(KIND(1.D0)),INTENT(IN)					:: SBLOC,SIGMOY,C1,SIG
	INTEGER,INTENT(IN)						:: NP1,NP2,N1,N,M
	INTEGER,INTENT(IN)						:: NW1,NW2,NW3
	REAL(KIND(1.D0)),DIMENSION (N1:NP1,N1:NP2),INTENT(IN)		:: VME
	REAL(KIND(1.D0)),DIMENSION(N1:NP1,N1:NP2,6),INTENT(INOUT)	:: XME
!
	DO J=N1,M
	   DO I=N1,N
	      XME(I,J,1)=XME(I,J,1)+VME(I,J)*SIG
	      XME(I,J,2)=XME(I,J,2)+VME(I,J)*SIG
	   END DO
	END DO	
	IF(NW1.EQ.1) THEN
	   DO J=N1,M
	      DO I=N1,N
	         XME(I,J,3)=XME(I,J,3)+XME(I,J,2)/SBLOC
	         XME(I,J,4)=XME(I,J,4)+(XME(I,J,2)/SBLOC)**2
	      END DO
	   END DO	
	ENDIF
	IF(NW1.EQ.1.AND.NW3.EQ.1) XME(N1:N,N1:M,2)=0.D0
	IF(NW2.EQ.1) THEN
	   DO J=N1,M
	      DO I=N1,N
	         XME(I,J,5)=XME(I,J,1)/SIGMOY
	         XME(I,J,6)=XME(I,J,4)/C1-(XME(I,J,3)/C1)**2
	         XME(I,J,6)=DSQRT(XME(I,J,6)/C1)
	      END DO
	   END DO	
	ENDIF
	END SUBROUTINE MESURE3
!--------------------------------------------------------------------------------------
! Nom		: ENTREE 	[SUB]						4
! Fonction	: Lecture des parametres et verification de compatibilite
! Auteurs	: H.Touchette,AM.Tremblay,A.Reid,H.Nelisse,C.Bourbonnais,C.Boily
! Revision	: 15 juillet 1997 (Touchette)
! Notes		: - Le nombre d'option (ZETA) = NOPT
!		  - Le module CONSTANTES fixe toutes les bornes des parametres 
!		  - Voir le fichier d'entree pour une description des parametres
!--------------------------------------------------------------------------------------	
	SUBROUTINE ENTREE(PARA)
	USE FICHIERS,ONLY: FICIN
	USE CONSTANTES
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER					:: I,J,K,FLAG,NR	! Var. locales
	TYPE(PARAMETRE),INTENT(OUT) 		:: PARA
	CHARACTER				:: BUFFER
!
	J=LEN_TRIM(FICIN)
	OPEN(UNIT=3,FILE=FICIN(1:J),STATUS='OLD',ACCESS='SEQUENTIAL',ERR=666)
	REWIND(UNIT=3)
!
	READ(3,*,ERR=650) PARA%BETA
	READ(3,*,ERR=651) PARA%TRAN
        READ(3,*,ERR=652) PARA%U
        READ(3,*,ERR=653) PARA%MU
        READ(3,*,ERR=654) PARA%REMPLD
        READ(3,*) PARA%TPAR
        READ(3,*) PARA%TPER
        READ(3,*) PARA%TPRM
        READ(3,*) PARA%N
        READ(3,*) PARA%M
        READ(3,*) PARA%NVEC
	READ(3,*) BUFFER
	IF(PARA%NVEC.NE.0) THEN
	   DO I=1,PARA%NVEC
	      READ(3,*) PARA%VEC(I,1),PARA%VEC(I,2)
	   END DO
	END IF
        READ(3,*) PARA%MES
        READ(3,*) PARA%NWARM
        READ(3,*) PARA%NMODUL
        READ(3,*) PARA%NBLOC
	READ(3,*) PARA%NIMPRI
        READ(3,*) PARA%NOPT
	DO I=1,PARA%NOPT
	   READ(3,*) PARA%ZETA(I)
	END DO
	READ(3,*) PARA%NR1
	READ(3,*) PARA%NR2
	READ(3,*) PARA%NR3
	READ(3,*) PARA%NR4
	NR=PARA%NR1+PARA%NR2+PARA%NR3+PARA%NR4
	DO I=1,NR
	   READ(3,*) PARA%YNMES(I)
	END DO
	
	IF(PARA%NVEC==0) THEN
	   PARA%NVEC=PARA%N*PARA%M
	   DO I=1,PARA%N
	      DO J=1,PARA%M
	         K=PARA%N*(J-1)+I
	         PARA%VEC(K,1)=I
	         PARA%VEC(K,2)=J
	      END DO
	   END DO
	END IF
!
!	Confirmation de l'entree a l'unite '*'
        WRITE(*,6) PARA%BETA,		'BETA'
        WRITE(*,5) PARA%TRAN,		'TRAN'
        WRITE(*,6) PARA%U,		'U'
        WRITE(*,6) PARA%MU,		'MU'
        WRITE(*,6) PARA%REMPLD,		'REMPLD'
        WRITE(*,6) PARA%TPAR,		'TPAR'
        WRITE(*,6) PARA%TPER,		'TPER'
        WRITE(*,6) PARA%TPRM,		'TPRM'
        WRITE(*,5) PARA%N,		'N'
        WRITE(*,5) PARA%M,		'M'
        WRITE(*,5) PARA%NVEC,	'Nb de vecteurs'
	WRITE(*,*) 'KX        KY'
	DO I=1,PARA%NVEC
	   WRITE(*,*) PARA%VEC(I,1),'        ',PARA%VEC(I,2)
	END DO
        WRITE(*,5) PARA%MES,		'MES'
        WRITE(*,5) PARA%NWARM,		'NWARM'
        WRITE(*,5) PARA%NMODUL,		'NMODUL'
        WRITE(*,5) PARA%NBLOC,		'NBLOC'
        WRITE(*,5) PARA%NIMPRI,		'NIMPRI'
        WRITE(*,5) PARA%NOPT,		'NOPT'
        DO I=1,PARA%NOPT
           WRITE(*,*) PARA%ZETA(I)
        END DO
        WRITE(*,5) PARA%NR1,		'NR1'
        WRITE(*,5) PARA%NR2,		'NR2'
        WRITE(*,5) PARA%NR3,		'NR3'
        WRITE(*,5) PARA%NR4,		'NR4'
        DO I=1,NR
           WRITE(*,*) PARA%YNMES(I)
	END DO
	
!
!	Verifications
	FLAG=0
166	IF(MOD(PARA%TRAN,PARA%NMODUL).NE.0)   				GOTO 667
168	IF(MOD(PARA%MES,PARA%NBLOC).NE.0)				GOTO 668
!169	IF(PARA%NBLOC.GT.PARA%MES) 					GOTO 669
170	IF(PARA%N.GT.ND1) 						GOTO 670
171	IF(PARA%M.GT.ND2) 						GOTO 671
172	IF(PARA%TRAN.GT.NL) 						GOTO 672
173	IF(PARA%MES/PARA%NBLOC.GT.NCOR) 				GOTO 673
174	IF(PARA%NR1.GT.NSC)						GOTO 674
175	IF(PARA%NR2.GT.NTE)						GOTO 675
176	IF(PARA%NR3.GT.NTI)						GOTO 676
177	IF(PARA%NR4.GT.NTG)						GOTO 677
178	NR=PARA%NR1+PARA%NR2+PARA%NR3+1
	IF((PARA%YNMES(NR+1).EQ.'Y'.OR.PARA%YNMES(NR+1).EQ.'y').AND.&
	  &(PARA%YNMES(NR).NE.'Y'.AND.PARA%YNMES(NR).NE.'y'))		GOTO 678
179	IF((PARA%YNMES(NR+2).EQ.'Y'.OR.PARA%YNMES(NR+2).EQ.'y').AND.&
	  &(PARA%YNMES(NR+1).NE.'Y'.AND.PARA%YNMES(NR+1).NE.'y'))	GOTO 679
180	IF((PARA%ZETA(5).EQ.'Y'.OR.PARA%ZETA(5).EQ.'y').AND.&
	  &(PARA%YNMES(NR).NE.'Y'.AND.PARA%YNMES(NR).NE.'y'))		GOTO 680
181	NR=PARA%NR1+PARA%NR2+1
	IF((PARA%YNMES(PARA%NR1+5).EQ.'Y'.OR.PARA%YNMES(PARA%NR1+5).EQ.'y')&
	  &.AND.(PARA%YNMES(NR).NE.'Y'.AND.PARA%YNMES(NR).NE.'y'))	GOTO 681
182	IF((PARA%ZETA(6).EQ.'Y'.OR.PARA%ZETA(6).EQ.'y').AND.&
	  &(PARA%YNMES(NR).NE.'Y'.AND.PARA%YNMES(NR).NE.'y'))		GOTO 682
183     IF((PARA%ZETA(7).EQ.'Y'.OR.PARA%ZETA(7).EQ.'y').AND.&
          &(PARA%YNMES(NR+3).NE.'Y'.AND.PARA%YNMES(NR+3).NE.'y'))	GOTO 683
184	IF(FLAG.EQ.1) STOP
!
	CLOSE(UNIT=3)
!
5	FORMAT(I10,A10)
6	FORMAT(F10.5,A10)
	RETURN
!	Messages d'erreurs
650	WRITE(*,*) '--- Erreur : Incapable de lire BETA ---';STOP
651	WRITE(*,*) '--- Erreur : Incapable de lire TRAN ---';STOP
652	WRITE(*,*) '--- Erreur : Incapable de lire U ---';STOP
653	WRITE(*,*) '--- Erreur : incapable de lire MU ---';STOP
654	WRITE(*,*) '--- Erreur : Incapable de lire REMPLD ---';STOP
666	WRITE(*,*) '--- Erreur : Fichier inexistant (ENTREE) ---'; STOP
667	WRITE(*,*) '--- Erreur : TRAN pas multiple de NMODUL ---'; 	FLAG=1;GOTO 168
668	WRITE(*,*) '--- Erreur : MES pas multiple NBLOC ---'; 	  	FLAG=1;GOTO 170
!669	WRITE(*,*) '--- Erreur : NBLOC > MES ---'; 		   	FLAG=1;GOTO 170
670	WRITE(*,*) '--- Erreur : N > ND1 ---'; 			   	FLAG=1;GOTO 171
671	WRITE(*,*) '--- Erreur : M > ND2 ---'; 			   	FLAG=1;GOTO 172
672	WRITE(*,*) '--- Erreur : TRAN > NL ---'; 			FLAG=1;GOTO 173
673	WRITE(*,*) '--- Erreur : Nb de blocs > NCOR ---'; 		FLAG=1;GOTO 174 
674	WRITE(*,*) '--- Erreur : NR1 > NSC ---'; 			FLAG=1;GOTO 175
675	WRITE(*,*) '--- Erreur : NR2 > NTE ---'; 			FLAG=1;GOTO 176
676	WRITE(*,*) '--- Erreur : NR3 > NTI ---'; 			FLAG=1;GOTO 177
677	WRITE(*,*) '--- Erreur : NR4 > NTG ---'; 			FLAG=1;GOTO 178
678	WRITE(*,*) '--- Erreur : MATSUBARA sans fct de Green ---'; 	FLAG=1;GOTO 179
679	WRITE(*,*) '--- Erreur : SELF sans MATSUBARA ---'; 		FLAG=1;GOTO 180
680	WRITE(*,*) '--- Erreur : COVARIANCE sans fcts de GREEN ---';	FLAG=1;GOTO 181
681	WRITE(*,*) '--- Erreur : Integrale sur tau de SUSCEPI sans SUSCEPI ---';FLAG=1;GOTO 182
682	WRITE(*,*) '--- Erreur : Covariance sans suscep de paire ---';	FLAG=1;GOTO 183
683     WRITE(*,*) '--- Erreur : Covariance sans Fct de corr de courant ---';FLAG=1;GOTO 184
	END SUBROUTINE
!--------------------------------------------------------------------------------------
! Nom		: SORTIE1 	[SUB]						5-(1)
! Fonction	: Impressions des resultats
! Auteurs	: H.Touchette
! Revision	: 15 juillet 1997 (Touchette)
! Notes		: - UNIT=12 
!		    
!--------------------------------------------------------------------------------------
	SUBROUTINE SORTIE1(PARA,MESC,COMP,SIGMOY,FLUCT,XMUPROP,ENERG,OX,VX)	
	USE CONSTANTES,ONLY: NXC
	USE PARA_SIM
	USE IDENT
	IMPLICIT NONE
	INTEGER						:: NM,MESU,I
	REAL(KIND(1.D0)),INTENT(IN)			:: SIGMOY
	REAL(KIND(1.D0)),INTENT(IN)			:: FLUCT,XMUPROP,ENERG
	INTEGER,INTENT(IN)				:: MESC	
	REAL(KIND(1.D0)),DIMENSION (2)			:: COMP	
	REAL(KIND(1.D0)),DIMENSION (0:NXC),INTENT(IN)	:: OX,VX
	TYPE(PARAMETRE),INTENT(IN) 			:: PARA
!
	NM=PARA%N*PARA%M
!
!	Sortie des parametres et resultats au fichier (UNIT=12)
	WRITE(12,*) '------------------------------------------------------------'
	WRITE(12,*)  DTE,		'Date mm/dd/yy'
	WRITE(12,*) IDNUM,		'No. simulation'
	WRITE(12,*)  NANUM,		'Prefixe fichier'
	WRITE(12,*) ' '
	WRITE(12,*) 'Parametres de simulation'
	WRITE(12,*) ' '
	WRITE(12,17) PARA%BETA,		'BETA'
	WRITE(12,15) PARA%TRAN,		'TRAN'
	WRITE(12,17) PARA%U,		'U'
	WRITE(12,17) PARA%MU,		'MU'
	WRITE(12,17) PARA%REMPLD,	'Remplissage desire'
	WRITE(12,17) PARA%TPAR,		'TPAR'
	WRITE(12,17) PARA%TPER,		'TPER'
	WRITE(12,17) PARA%TPRM,		'TPRM'
	WRITE(12,15) PARA%N,		'Nombre de sites par chaine N'
	WRITE(12,15) PARA%M,		'Nombre de chaines M'
	WRITE(12,15) PARA%NVEC,	'Nb de vecteurs'
	WRITE(12,*) 'KX        KY'
	DO I=1,PARA%NVEC
	   WRITE(12,16) PARA%VEC(I,1),PARA%VEC(I,2)
	END DO
	WRITE(12,15) PARA%MES,		'Nombre total de mesures'
	WRITE(12,15) PARA%NWARM,		'Nombre de rechauffement'
	WRITE(12,15) PARA%NMODUL,	'NMODUL : Nb de tranche de stabilisation'
	WRITE(12,15) PARA%NBLOC,		'Nombre de mesures dans un bloc'
	WRITE(12,15) PARA%NIMPRI,	'Borne impression'
	WRITE(12,15) PARA%NOPT,		'Nombre option de calcul'  
	WRITE(12,13,ADVANCE='NO') PARA%ZETA(1)
	WRITE(12,*)		'Fct autocorrelation'
	WRITE(12,13,ADVANCE='NO') PARA%ZETA(2)
	WRITE(12,*)		'Lecture configuration HS initiale'
	WRITE(12,13,ADVANCE='NO') PARA%ZETA(3)
	WRITE(12,*)		'Sauvegarde des configuration HS'
	WRITE(12,13,ADVANCE='NO') PARA%ZETA(4)
	WRITE(12,*)		'Mouvements globaux'
	WRITE(12,13,ADVANCE='NO') PARA%ZETA(5)
	WRITE(12,*)		'Matrice de covariance'
	WRITE(12,13,ADVANCE='NO') PARA%ZETA(6)
	WRITE(12,*)		'Matrice de cov de suscep'
        WRITE(12,13,ADVANCE='NO') PARA%ZETA(7)
        WRITE(12,*)              'Matrice de cov de courant'
	WRITE(12,15) PARA%NR1,		'Nombre de routines scalaire'
	WRITE(12,15) PARA%NR2,		'Nombre de routines temps egal'
	WRITE(12,15) PARA%NR3,		'Nombre de routines temps inegal'
	WRITE(12,15) PARA%NR4,		'Nombre de routines fct de Green'
	
	DO I=1,PARA%NR1+PARA%NR2+PARA%NR3+PARA%NR4
	    IF(PARA%YNMES(I).EQ.'Y'.OR.PARA%YNMES(I).EQ.'y')THEN
		MESU=MESC
	    ELSE
		MESU=0
	    ENDIF
	    WRITE(12,14,ADVANCE='NO') PARA%YNMES(I),MESU
	    SELECT CASE (I)
	    CASE (1)
		WRITE(12,*)	'DSITE'
	    CASE (2)
		WRITE(12,*)	'MAGST'
	    CASE (3)
		WRITE(12,*)	'DENSITE'
	    CASE (4)
		WRITE(12,*)	'RHOE'
	    CASE (5)
		WRITE(12,*)	'PSCOR'
	    CASE (6)
		WRITE(12,*)	'Integrale sur tau de SUSCEPI'
	    CASE (7)
		WRITE(12,*)	'SUSCEPI'
	    CASE (8)
		WRITE(12,*)	'SUSCEPIe'
	    CASE (9)
		WRITE(12,*)	'SUSCEPIh'
            CASE (10)
                WRITE(12,*)	'COURANT (tau,Qx,Qy)'
	    CASE (11)
		WRITE(12,*)	'SUSMAG'
	    CASE (12)
		WRITE(12,*)	'Fct de Green (qx,qy)'
	    CASE (13)
		WRITE(12,*)	'MATSUBARA G(iwn,Qx,Qy)'
	    CASE (14)
		WRITE(12,*)	'SELF'
	    END SELECT
	END DO
	WRITE(12,*) ' '
	WRITE(12,*) 'Parametres de la bande'
	WRITE(12,*) ' '
	WRITE(12,20) OX(0),VX(0),	'Remplissage mesure'
	WRITE(12,21) XMUPROP,		'Mu propose'
	WRITE(12,21) COMP(1),		'Proportion de spin-flip accepte'
	WRITE(12,21) COMP(2),		'Proportion de global move accepte'
	WRITE(12,21) SIGMOY/DBLE(MESC),	'Signe moyen du produit des det.'
	WRITE(12,*) ' '
	WRITE(12,20) OX(1)/NM,VX(1)/NM,	'Energie cinetique'
	WRITE(12,20) OX(2)/NM,VX(2)/NM,	'Energie totale'
	IF(PARA%YNMES(1).EQ.'Y'.OR.PARA%YNMES(1).EQ.'y')THEN
	   WRITE(12,21) ENERG,		'Energie des electrons libres'
	ENDIF
	WRITE(12,21) FLUCT,		'Fluctuation du nombre de particules'
!
!	Format d'impression
13	FORMAT(A2,2x)
14	FORMAT(A2,2X,I6,2X)
15	FORMAT(I10,2x,A)
16	FORMAT(I3,7X,I3)
17	FORMAT(F10.6,2x,A)
20	FORMAT(F10.7,1x,F10.7,2x,A)
21	FORMAT(F10.7,2x,A)
!
	END SUBROUTINE SORTIE1
!--------------------------------------------------------------------------------------
! Nom		: SORTIE2 	[SUB]						5-(2)
! Fonction	: Impressions des facteurs de structures a temps egal et inegal
! Auteurs	: H.Touchette
! Revision	: 20 juillet 1997 (Touchette)
! Notes		: - Utilise WRITEOUT par l'intermediare de l'interface INTER
!                 - UNIT=12
!--------------------------------------------------------------------------------------
	SUBROUTINE SORTIE2(PARA,X,VX,XI,VXI)
	USE CONSTANTES,ONLY: NL,ND1,ND2,NTE,NTI
	USE PARA_SIM
	USE INTER
	IMPLICIT NONE
	INTEGER						:: NR
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,NTE)	:: X,VX
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,1:NTI)	:: XI,VXI
	TYPE(PARAMETRE),INTENT(IN)			:: PARA 
!
	NR=2
	WRITE(12,*) ' ' 
	WRITE(12,*) 'Quantites a temps egal'
	WRITE(12,*) ' '
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Susceptibilite magnetique (MAGST)'
	   WRITE(12,*) 'qx    qy      x         err(x)'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,E=X(:,:,1),VE=VX(:,:,1))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Densite impulsion (DENSITE)'
	   WRITE(12,*) 'qx    qy      k         err(k)'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,E=X(:,:,2),VE=VX(:,:,2))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Densite de charge (RHOE)'
	   WRITE(12,*) 'qx    qy      x         err(x)'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,E=X(:,:,3),VE=VX(:,:,3))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Susceptibilite de paires (PSCOR)'
	   WRITE(12,*) 'qx    qy      x         err(x)'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,E=X(:,:,4),VE=VX(:,:,4))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Integrale sur tau de la susceptibilite de paire dynamique'
	   WRITE(12,*) 'qx    qy      x         err(x)'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,E=X(:,:,5),VE=VX(:,:,5))
	ENDIF
	NR=PARA%NR1+PARA%NR2+1
	WRITE(12,*) 'Quantites a temps inegal'
	WRITE(12,*) ' '
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Susceptibilite de paires dynamique (SUSCEPI)'
	   WRITE(12,*) ' qx  qy  tau     X            err'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,F=XI(:,:,:,1),VF=VXI(:,:,:,1))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Susceptibilite de paires electrons (SUSCEPIe)'
	   WRITE(12,*) ' qx  qy  tau     X            err'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,F=XI(:,:,:,2),VF=VXI(:,:,:,2))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(12,*) 'Susceptibilite de paires de trous (SUSCEPIh)'
	   WRITE(12,*) ' qx  qy  tau     X            err'
	   WRITE(12,*) ' '
	   CALL WRITEOUT(PARA,12,F=XI(:,:,:,3),VF=VXI(:,:,:,3))
	ENDIF
	NR=NR+1
        IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y')THEN
           WRITE(12,*) 'Fonction de correlation courant-courant'
           WRITE(12,*) ' qx  qy  tau     X            err'
           WRITE(12,*) ' '
           CALL WRITEOUT(PARA,12,F=XI(:,:,:,4),VF=VXI(:,:,:,4))
        ENDIF
        NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y')THEN
	    WRITE(12,*) 'Susceptibilité magnétique'
	    WRITE(12,*) ' qx  qy  tau    X            err'
	    WRITE(12,*) ' '
	    CALL WRITEOUT(PARA,12,F=XI(:,:,:,5),VF=VXI(:,:,:,5))
	ENDIF
!
	END SUBROUTINE SORTIE2
!--------------------------------------------------------------------------------------
! Nom		: SORTIE3 	[SUB]						5-(3)
! Fonction	: Impressions des fonctions de Green
! Auteurs	: H.Touchette
! Revision	: 20 juillet 1997 (Touchette)
! Notes		: - Utilise WRITEOUT par l'intermediare de l'interface INTER
!		  - UNIT=11
!--------------------------------------------------------------------------------------
	SUBROUTINE SORTIE3(PARA,MESC,GU,GD,VGU,VGD)
	USE CONSTANTES,ONLY: NL,ND1,ND2,NTG
	USE PARA_SIM
	USE INTER
	USE IDENT,ONLY: IDNUM
	IMPLICIT NONE
	INTEGER								:: NR
	INTEGER,INTENT(IN)						:: MESC
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,NTG),INTENT(IN)	:: GU,GD 
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL,NTG),INTENT(IN)	:: VGU,VGD
	TYPE(PARAMETRE),INTENT(IN)					:: PARA 
!
	WRITE(11,*) '------------------------------------------------------------'
	WRITE(11,97) IDNUM,		'No. simulation'
	WRITE(11,97) MESC,		'Nb mesures effectuees'
	WRITE(11,*) ' '
	NR=PARA%NR1+PARA%NR2+PARA%NR3+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Fonctions de Green, espace reciproque'
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Spin-up'
	   WRITE(11,*) ' qx  qy  tau     G            err'
	   CALL WRITEOUT(PARA,11,F=GU(:,:,:,1),VF=VGU(:,:,:,1))
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Spin-down'
	   WRITE(11,*) ' qx  qy  tau     G            err'
	   CALL WRITEOUT(PARA,11,F=GD(:,:,:,1),VF=VGD(:,:,:,1))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Fonctions de Green de Matsubara, espace reciproque'
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Partie reelle'
	   WRITE(11,*) ' qx  qy  iwn     G            err'
	   CALL WRITEOUT(PARA,11,F=GU(:,:,:,2),VF=VGU(:,:,:,2))
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Partie imaginaire'
	   WRITE(11,*) ' qx  qy  iwn     G            err'
	   CALL WRITEOUT(PARA,11,F=GD(:,:,:,2),VF=VGD(:,:,:,2))
	ENDIF
	NR=NR+1
	IF(PARA%YNMES(NR).EQ.'Y'.OR.PARA%YNMES(NR).EQ.'y') THEN
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Self-energie, espace reciproque'
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Partie reelle'
	   WRITE(11,*) ' qx  qy  iwn     S            err'
	   CALL WRITEOUT(PARA,11,F=GU(:,:,:,3),VF=VGU(:,:,:,3))
	   WRITE(11,*) ' '
	   WRITE(11,*) 'Partie imaginaire'
	   WRITE(11,*) ' qx  qy  iwn     S            err'
	   CALL WRITEOUT(PARA,11,F=GD(:,:,:,3),VF=VGD(:,:,:,3))
	ENDIF
	WRITE(11,*) ' '
97	FORMAT(I9,12X,A)
98 	FORMAT(2(I4),2(I6),3x,F10.7)
!
	END SUBROUTINE SORTIE3
!--------------------------------------------------------------------------------------
! Nom		: SORTIEC 	[SUB]						5-(4)
! Fonction	: Impressions de la fonction d'autocorrelation et de la matrice 
!		  de covariance
! Auteurs	: H.Touchette
! Revision	: 6 aout 1997 (Touchette)
! Notes		: - UNIT 15
!		  - Les valeurs de la matrice de covariance pour la derniere tranche
!		    de temps imaginaire ne sont pas imprimees du fait qu'elles 
!		    derivent exactement des valeurs a tau = 0. 
!--------------------------------------------------------------------------------------
	SUBROUTINE SORTIEC(PARA,MESC,C1,AUTOCOR,COVM,COVX,COVJ)
	USE CONSTANTES,ONLY: NL,ND1,ND2,NCOR 
	USE FICHIERS,ONLY: FICCVOUT
	USE IDENT,ONLY: IDNUM
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER								:: I,J,L1,L2,K
	INTEGER								:: NC1
	REAL(KIND(1.D0)),INTENT(IN)					:: C1
	INTEGER,INTENT(IN)						:: MESC
	REAL(KIND(1.D0)),DIMENSION (NCOR),INTENT(IN)			:: AUTOCOR
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:Nl,0:NL),INTENT(IN),OPTIONAL:: COVM,COVX,COVJ
	TYPE(PARAMETRE),INTENT(IN)					:: PARA 
!
	NC1=INT(C1)
	OPEN(UNIT=15,FILE=FICCVOUT)
	REWIND(UNIT=15)
	WRITE(15,77) IDNUM,		'No. simulation'
	WRITE(15,77) MESC,		'Nb mesures effectuees'
	WRITE(15,*) ' '
	IF(PARA%ZETA(5).EQ.'Y'.OR.PARA%ZETA(5).EQ.'y') THEN
	   WRITE(15,*) 'Matrice de covariance spin-up'
	   WRITE(15,*) ' '
	   WRITE(15,77) NC1,		'Nb blocs'
	   WRITE(15,*) ' '
	   WRITE(15,*) ' qx   qy    tau1   tau2         C'
	   DO K=1,PARA%NVEC
	      J=PARA%VEC(K,1)
	      I=PARA%VEC(K,2)
	      DO L1=0,PARA%TRAN-1
	         DO L2=0,PARA%TRAN-1
	            WRITE(15,78) I,J,L1,L2,COVM(J,I,L2,L1)
	         END DO
	      END DO
	   END DO
	ENDIF
	IF(PARA%ZETA(6).EQ.'y'.OR.PARA%ZETA(6).EQ.'Y')THEN
	   WRITE(15,*)' '
	   WRITE(15,*)'Matrice de covariance de la susceptibilite de paires'
	   WRITE(15,*)' '
	   WRITE(15,77) NC1,		'Nb blocs'
	   WRITE(15,*)' '
	   WRITE(15,*)' qx   qy   tau1   tau2        C'
	   DO K=1,PARA%NVEC
	      J=PARA%VEC(K,1)
	      I=PARA%VEC(K,2)
	      DO L1=0,PARA%TRAN-1
	        DO L2=0,PARA%TRAN-1
	           WRITE(15,78) I,J,L1,L2,COVX(J,I,L2,L1)
	        END DO
  	      END DO
	   END DO
	ENDIF
        IF(PARA%ZETA(7).EQ.'y'.OR.PARA%ZETA(7).EQ.'Y')THEN
           WRITE(15,*)' '
           WRITE(15,*)'Matrice de covariance de la Fct de corr de courant'
           WRITE(15,*)' '
           WRITE(15,77) NC1,            'Nb blocs'
           WRITE(15,*)' '
           WRITE(15,*)' qx   qy   tau1   tau2        C'
           DO K=1,PARA%NVEC
	      J=PARA%VEC(K,1)
	      I=PARA%VEC(K,2)
              DO L1=0,PARA%TRAN-1
                 DO L2=0,PARA%TRAN-1
                    WRITE(15,78) I,J,L1,L2,COVJ(J,I,L2,L1)
                 END DO
              END DO
           END DO
        ENDIF
	IF(PARA%ZETA(1).EQ.'Y'.OR.PARA%ZETA(1).EQ.'y') THEN
	   WRITE(15,*) 'Fonction autocorrelation'
	   WRITE(15,*) ' '
	   WRITE(15,77) NC1,		'Nb blocs'
	   WRITE(15,*) ' '
	   WRITE(15,*) ' n         C(n)'
	   DO I=1,NC1
	      WRITE(15,79) I-1,AUTOCOR(I)
	   END DO	               
	ENDIF
	CLOSE(UNIT=15)
!
77	FORMAT(I9,12X,A)
78	FORMAT(2(I4,1X),2(I6,1X),E16.4)
79	FORMAT(I4,1X,E16.4)
	END SUBROUTINE SORTIEC
!--------------------------------------------------------------------------------------
! Nom		: SORTIE4 	[SUB]						5-(5)
! Fonction	: Temps de calcul et archivage des simulations pour les pages web
! Auteurs	: H.Touchette
! Revision	: 22 juillet 1997 (Touchette)
! Notes		: UNIT=12  : Temps de calcul
!		  UNIT=13 : Archivage
!--------------------------------------------------------------------------------------
	SUBROUTINE SORTIE4(PARA,XMU,REMP)
	USE FICHIERS,ONLY: FICWEB,FICPAOUT,FICHSOUT,FICGIOUT,REPWEB
	USE PARA_SIM
	USE IDENT
	IMPLICIT NONE
	REAL(KIND(1.D0)),PARAMETER	:: spj=8.64d+4,mpj=1.44d+3
	REAL(KIND(1.D0))		:: TD
	REAL(KIND(1.D0)),INTENT(IN)	:: XMU,REMP
	INTEGER				:: IO,IH,IG,IR		! Var locales
	INTEGER				:: TF,CR,CM,IDTE,FDTE
	INTEGER				:: ASCIITONUM
	CHARACTER(LEN=8)		:: DTEF
	TYPE(PARAMETRE),INTENT(IN)	:: PARA
!
	CALL SYSTEM_CLOCK(TF,CR,CM)				! Heure arrivee
	IF(TF.LT.TI)THEN
	   TF=TF+CM
	ENDIF
	TD=REAL(TF-TI,KIND(1.D0))/(6.D1*CR)			! Diff en minute
	CALL DATE_AND_TIME(DTEF)
	IDTE=ASCIITONUM(DTE)
	FDTE=ASCIITONUM(DTEF)
	IF((FDTE-IDTE)*spj.gt.real(CM,kind(1.d0))/real(CR,kind(1.d0)))THEN
	   TD=mod(TD,mpj)+mpj*(FDTE-IDTE)
	ENDIF
	WRITE(12,*) ' '
	WRITE(12,*) TD,'  = Temps de calcul (minutes CPU)'
!
	IO=INDEX(FICPAOUT,' ')
	IH=INDEX(FICHSOUT,' ')
	IG=INDEX(FICGIOUT,' ')
	IR=INDEX(REPWEB,' ')
	OPEN(UNIT=17,FILE=FICWEB,POSITION="APPEND")
	WRITE(17,*) '<PRE>'
	WRITE(17,102) DTE,IDNUM,PARA%N,PARA%M,PARA%U,PARA%BETA,PARA%MU,XMU,REMP,&
	              &PARA%MES,PARA%NWARM,PARA%NBLOC,TD,REPWEB
	WRITE(17,103) '<A HREF="',REPWEB(1:IR-1),FICPAOUT(1:IO-1),'">',&
		&FICPAOUT(1:IO-1),'</A>'
	WRITE(17,103) '<A HREF="',REPWEB(1:IR-1),FICHSOUT(1:IH-1),'">',&
		&FICHSOUT(1:IH-1),'</A>'
	WRITE(17,103) '<A HREF="',REPWEB(1:IR-1),FICGIOUT(1:IG-1),'">',&
		&FICGIOUT(1:IG-1),'</A>'
	WRITE(17,*) '</PRE>'
	CLOSE(UNIT=17)
102	FORMAT(A8,1X,I8,2(1X,I2),5(1X,F6.2),3(1X,I8),1X,E16.2,6X,A)
103	FORMAT(106X,A9,A,A,A2,A,A4)
	END SUBROUTINE SORTIE4
!--------------------------------------------------------------------------------------
! Nom		: WRITEOUT		[SUB]					6
! Fonction	: Ecriture des resultats
! Auteurs	: H.Touchette
! Revision	: 20 juin 1997 (Touchette)
! Notes		: - Cette routine a ete construite de facon a illustrer l'utilisation
!		    d'arguments optionnels a l'aide d'une interface explicite
!		  - Utilitaires de SORTIE2
!		  - E : Resultats a temps egal UNIT=12 
!		  - F : Resultats a temps inegal UNIT 12 ou 11 (NU)
!--------------------------------------------------------------------------------------
	SUBROUTINE WRITEOUT(PARA,NU,E,VE,F,VF)
	USE CONSTANTES,ONLY: ND1,ND2,NL
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER								:: I,J,L,K
	INTEGER,INTENT(IN)						:: NU
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(IN),OPTIONAL	:: E,VE
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL),INTENT(IN),OPTIONAL	:: F,VF
	TYPE(PARAMETRE),INTENT(IN)					:: PARA
!
	IF(PRESENT(E)) THEN
	   DO K=1,PARA%NVEC
	      J=PARA%VEC(K,1)
	      I=PARA%VEC(K,2)
	      WRITE(NU,324) I,J,E(J,I),VE(J,I)
	   END DO
	ELSE IF(PRESENT(F)) THEN
	   DO K=1,PARA%NVEC
	      J=PARA%VEC(K,1)
	      I=PARA%VEC(K,2)
	      DO L=0,PARA%TRAN
	         WRITE(NU,325) I,J,L,F(J,I,L),VF(J,I,L)
	      END DO
	   END DO
	ELSE
	   WRITE(*,*) '--- Erreur : WRITEOUT ---'
	ENDIF
324	FORMAT(2(1X,I3),2(2X,F14.6))
325	FORMAT(3(1X,I3),2(1X,F14.6))
	END SUBROUTINE WRITEOUT
!--------------------------------------------------------------------------------------
! Nom		: GRASCHP	[SUB]						7
! Fonction	: Decomposition UDV par la methode de Gram-Schmidt modifiee + pivot
! Auteurs	: 
! Revision	: 
! References	: S.Moukouri
! Notes		: 
! Dependances	: ORTHOGP,DONORM (internes)
!--------------------------------------------------------------------------------------
	SUBROUTINE GRASCHP(NP,NM,AMAT,DVEC,RMAT)
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: TEMP		! Variable locale
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER						:: I,J		! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP),INTENT(OUT)	:: DVEC
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT)	:: RMAT
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: AMAT
!
 	CALL ORTHOGP(NP,NM,AMAT,RMAT)
 	DO I=1,NM
 	   DVEC(I)=RMAT(I,I)
 	   IF(DVEC(I).EQ.0.0) THEN
 	      RMAT(I,I)=1.D0
 	   ELSE IF(DVEC(I).NE.1.0) THEN
 	      TEMP=1.D0/DVEC(I)
 	      DO J=1,NM
 	         RMAT(I,J)=RMAT(I,J)*TEMP
 	      END DO
 	   END IF
 	END DO
 	RETURN
 	END
!
!	Routines internes
!
 	SUBROUTINE ORTHOGP(NP,NM,AMAT,RMAT)
 	IMPLICIT NONE
 	REAL(KIND(1.D0)) 				:: AMAXSIZE,COLSIZE	! Loc.
 	INTEGER,INTENT(IN)				:: NP,NM
 	INTEGER						:: I,J,K,IB	! Variables locales
 	REAL(KIND(1.D0)),DIMENSION (NP)			:: IDONE	! Variable locale	
 	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: AMAT
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT)	:: RMAT
  	DO I=1,NM
 	   DO J=1,NM
 	      RMAT(I,J)=0.D0
 	   END DO
 	   RMAT(I,I)=1.D0
 	END DO
  	DO I=1,NM
 	   IDONE(I)=0
 	END DO
 	DO I=1,NM
 	   AMAXSIZE=-10.0
 	   IB=-1
 	   DO J=1,NM
 	      IF(IDONE(J).EQ.0) THEN
 	         COLSIZE=0.0
 	         DO K=1,NM
 	            COLSIZE=COLSIZE+AMAT(K,J)**2
 	         END DO
 	         IF(COLSIZE.GT.AMAXSIZE) THEN
 	            AMAXSIZE=COLSIZE
 	            IB=J
 	         END IF
 	      END IF
 	   END DO
 	   IDONE(IB)=1
 	   CALL DONORM(NP,NM,AMAT(:,IB),RMAT(IB,IB))
 	   DO J=1,NM
 	      IF(IDONE(J).EQ.0) THEN
 	         DO K=1,NM
 	            RMAT(IB,J)=RMAT(IB,J)+AMAT(K,IB)*AMAT(K,J)
 	         END DO
 	         DO K=1,NM
 	            AMAT(K,J)=AMAT(K,J)-RMAT(IB,J)*AMAT(K,IB)
 	         END DO
 	      END IF
 	   END DO
 	END DO
 	END SUBROUTINE ORTHOGP
!
 	SUBROUTINE DONORM(NP,NM,VEC,VNORM)
 	IMPLICIT NONE
 	REAL(KIND(1.D0)) 				:: VNORM,TEMP	! Variables locales
 	INTEGER						:: I		! Variable locale
 	INTEGER,INTENT(IN)				:: NP,NM
 	REAL(KIND(1.D0)),DIMENSION (NP),INTENT(INOUT)	:: VEC
  	TEMP=0.D0
 	DO I=1,NM
 	   TEMP=TEMP+VEC(I)*VEC(I)
 	END DO
 	VNORM=DSQRT(TEMP)
 	IF(VNORM.NE.0.0) THEN
 	   TEMP=1.D0/VNORM
 	   DO I=1,NM
 	      VEC(I)=VEC(I)*TEMP
 	   END DO
 	END IF
 	END SUBROUTINE DONORM
!--------------------------------------------------------------------------------------
! Nom		: MULD		[SUB]						8
! Fonction	: Multiplication par une matrice diagonale
! Auteurs	: AM.Tremblay
! Revision	: 16 aout 1989 (Tremblay)
! Notes		: L=0 : multiplication par la gauche A = D * A (D:diag)
!		  L=2 : multiplication par la droite A = A * D
!--------------------------------------------------------------------------------------
	SUBROUTINE MULD(NP,NM,L,A,D)
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM,L
	INTEGER						:: I,J
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: A
	REAL(KIND(1.D0)),DIMENSION (NP),INTENT(IN)	:: D
	IF(L.EQ.0) THEN						! L = 0						
	   DO I=1,NM
	      DO J=1,NM
	         A(I,J)=D(I)*A(I,J)
	      END DO
	   END DO
	ELSE IF(L.EQ.2) THEN					! L = 2
	   DO J=1,NM
	      DO I=1,NM
	         A(I,J)=A(I,J)*D(J)
	      END DO
	   END DO
	ELSE							! Erreur
	   WRITE(*,*) '--- Erreur : MULD ---'
	ENDIF
	END SUBROUTINE MULD
!--------------------------------------------------------------------------------------
! Nom		: MULTD		[SUB]						9
! Fonction	: Multiplication par une matrice triangulaire droite
! Auteurs	: AM.Tremblay
! Revision	: 16 aout 1989 (Tremblay)
! Notes		: L=0 : multiplication par la gauche R = A(transpose) * T
!		  L=2 : multiplication par la droite R = T * A(transpose)
!			T = triangulaire droite NM*NM
!			A = Matrice source (par laquelle on fait la multiplication)
!			R = Resultat de la multiplication
!--------------------------------------------------------------------------------------
	SUBROUTINE MULTD(NP,NM,L,A,T,R)
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: S		! Variable locale
	INTEGER,INTENT(IN)				:: NP,NM,L
	INTEGER						:: I,J,K	! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: A
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: T
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT)	:: R
	IF(L.EQ.2) THEN						! L = 2
	   DO I=1,NM
	      DO J=1,NM
	         S=0.D0
	         DO K=I,NM
	            S=S+T(I,K)*A(J,K)
	         END DO
	         R(I,J)=S
	      END DO
	   END DO
	ELSE IF(L.EQ.0) THEN					! L = 0
	   DO I=1,NM
	      DO J=1,NM
	         S=0.D0
	         DO K=1,J
	            S=S+A(K,I)*T(K,J)
	         END DO
	         R(I,J)=S
	      END DO
	   END DO
	ELSE							! Erreur
	   WRITE(*,*) '--- Erreur : MULTD ---'
	ENDIF
	END SUBROUTINE MULTD
!--------------------------------------------------------------------------------------
! Nom		: MULTD2	[SUB]						10
! Fonction	: Multiplication de deux matrices triangulaires droites
! Auteurs	: AM.tremblay
! Revision	: 16 aout 1989 (Tremblay)
! Notes		: T2 = T1 * T2 (T1,T2: triangulaires droites NM*NM)
!--------------------------------------------------------------------------------------
	SUBROUTINE MULTD2(NP,NM,T1,T2)
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER						:: I,J,K,L	! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP)			:: R		! Variable locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: T1
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: T2
	DO I=1,NM
	   DO J=1,I
	      R(J)=0.D0
	      DO K=J,I
	         R(J)=R(J)+T1(J,K)*T2(K,I)
	      END DO
	   END DO
	   DO L=1,I
	      T2(L,I)=R(L)
	   END DO
	END DO
	END SUBROUTINE MULTD2
!--------------------------------------------------------------------------------------
! Nom		: MAC1		[SUB]						11-(1)
! Fonction	: A * B => B
! Auteurs	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC1(NP,NM,A,B)
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER						:: I,J,K,L	! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP)			:: C		! Variable locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: A
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: B
	DO I=1,NM
	   DO J=1,NM
	      C(J)=0.D0
	      DO K=1,NM
	         C(J)=C(J)+A(J,K)*B(K,I)
	      END DO
	   END DO
	   DO L=1,NM
	      B(L,I)=C(L)
	   END DO
	END DO
	END SUBROUTINE MAC1
!--------------------------------------------------------------------------------------
! Nom		: MAC2		[SUB]						11-(2)
! Fonction	: A * B => A
! Auteurs	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC2(NP,NM,A,B)
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER						:: I,J,K,L	! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP)			:: C		! Variable locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: A
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: B
	DO I=1,NM
	   DO J=1,NM
	      C(J)=0.D0
	      DO K=1,NM
	         C(J)=C(J)+A(I,K)*B(K,J)
	      END DO
	   END DO
	   DO L=1,NM
	      A(I,L)=C(L)
	   END DO
	END DO
	END SUBROUTINE MAC2
!--------------------------------------------------------------------------------------
! Nom		: MAC8		[SUB]						11-(8)
! Fonction	: CTE* A => B
! Auteurs	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC8(NP,NM,A,B,CTE)
	IMPLICIT NONE
	REAL(KIND(1.D0)),INTENT(IN)			:: CTE
	INTEGER,INTENT(IN)				:: NP,NM	
	INTEGER						:: I,J		! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: A
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT)	:: B
	DO J=1,NM
	   DO I=1,NM
	      B(I,J)=CTE*A(I,J)
	   END DO
	END DO
	END SUBROUTINE MAC8
!--------------------------------------------------------------------------------------
! Nom		: MAC9		[SUB]						11-(9)
! Fonction	: A + I => A
! Auteurs	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC9(NP,NM,A)
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER						:: I		! Variable locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: A
	DO I=1,NM
	   A(I,I)=A(I,I)+1.D0
	END DO
	END SUBROUTINE MAC9
!--------------------------------------------------------------------------------------
! Nom		: MAC10		[SUB]						11-(10)
! Fonction	: A => B
! Auteurs	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC10(NP,NM,A,B)
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER						:: I,J		! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: A
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT)	:: B
	DO J=1,NM
	   DO I=1,NM
	      B(I,J)=A(I,J)
	   END DO
	END DO
	END SUBROUTINE MAC10
!--------------------------------------------------------------------------------------
! Nom		: MAC11		[SUB]						11-(11)
! Fonction	: TRANSPOSE(A) => A
! Auteurs	:
! Revision	: 30 juillet 97 (Touchette)
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC11(NP,NM,A)
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: C 		! Var. locale
	INTEGER						:: I,J		! Var. locales
	INTEGER,INTENT(IN)				:: NP,NM
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: A			
	DO I=2,NM
	   DO J=1,I-1
	      C=A(I,J)
	      A(I,J)=A(J,I)
	      A(J,I)=C
	   END DO
	END DO
	END SUBROUTINE MAC11
!--------------------------------------------------------------------------------------
! Nom		: INVTD		[SUB]						12
! Fonction	: Inversion d'une matrice triangulaire droite
! Auteurs	: AM.Tremblay
! Revision	: 16 aout 1989 (Tremblay)
! Notes		: - A: triangulaire droite NM*NM
!		  - Ne sert pas dans le programme
!--------------------------------------------------------------------------------------
	SUBROUTINE INVTD(NP,NM,A)
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: S		! Variable locale
	INTEGER						:: I,J,K	! Variables locales
	INTEGER,INTENT(IN)				:: NP,NM
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: A			
        REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: R		! Variable locale 
        DO J=1,NM	
           DO I=1,NM
              R(I,J)=0.D0
           END DO
           R(I,I)=1.D0
        END DO
	DO K=NM,2,-1
	   DO I=K-1,1,-1
	      S=0.D0
	      DO J=I+1,K
	         S=S-A(I,J)*R(J,K)
	      END DO
	      R(I,K)=S
	   END DO
	END DO
	CALL MAC10(NP,NM,R,A)
	END SUBROUTINE INVTD
!--------------------------------------------------------------------------------------
! Nom		: UPDTG		[SUB]						13
! Fonction	: Mise a jour de la fonction de Green si un changement de spin HS
!		  a ete accepte.
! Auteurs	: AM.Tremblay
! Revision	:
! Reference	: Blankenbecler et al. (1981)
! Notes		:
!--------------------------------------------------------------------------------------	
	SUBROUTINE UPDTG(NP,NM,NS,G,VAR)
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: FACT		! Variables locales
	REAL(KIND(1.D0)),INTENT(IN)			:: VAR
	INTEGER,INTENT(IN)				:: NP,NM,NS
	INTEGER						:: I,J		! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP)			:: VV,VH	! Variables locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: G
!
	FACT=1.D0+(1.D0-G(NS,NS))*VAR
	FACT=VAR/FACT
	DO I=1,NM
	   VV(I)=G(I,NS)
	   VH(I)=G(NS,I)
	END DO
	DO I=1,NM
	   DO J=1,NM
	      G(I,J)=G(I,J)+VV(I)*VH(J)*FACT
	   END DO
	END DO
	DO J=1,NM
	   G(NS,J)=G(NS,J)-VH(J)*FACT
	END DO
	END SUBROUTINE UPDTG
!--------------------------------------------------------------------------------------
! Nom		: INITH3	[SUB]						14
! Fonction	: Initialisation de la matrice de hopping (forme dense diagonalisee)
! Auteurs	: Liang Chen,H.Touchette
! Revision	: 21 novembre 1996 (Touchette)
! Notes		: H est factorisee sous la forme H = ZZ(transpose) * DV * ZZ
!		  ZZ : matrice orthogonale des vecteurs propres de H
!		  DV : matrice diagonale des valeurs propres de H 
!		
!--------------------------------------------------------------------------------------
	SUBROUTINE INITH3(NM,PARA)
	USE CONSTANTES,ONLY: NP
	USE HOPPING
	USE PARA_SIM
	IMPLICIT NONE
	REAL(KIND(1.D0))		:: TX,TY,DNORM,DNORM1,PI,DTAU	! Var. locales
	REAL(KIND(1.D0))		:: DKX,DKY,PHASE,QX,QY,ENGY	! Var. locales
	INTEGER				:: II,I,J,IK,KX,KY,KXRF,KYRF	! Var. locales
	INTEGER				:: IKRF,IX,IR,IY,ICOUNT		! Var. locales
	REAL(KIND(1.D0)),DIMENSION (NP)	:: DVN				! Var. locale
	INTEGER,INTENT(IN)		:: NM
	INTEGER,DIMENSION (NM)		:: I_X,I_Y			! Var. locales
	TYPE(PARAMETRE),INTENT(IN)	:: PARA
!
	DTAU=PARA%BETA/DBLE(PARA%TRAN)
	PI=DACOS(-1.D0)
	TX=1.D0
	IF(PARA%N==2) TX=0.5D0
	TY=1.D0
	IF(PARA%M==1) TY=0.D0
	IF(PARA%M==2) TY=0.5D0
	DNORM=DSQRT(2.D0)/DSQRT(DBLE(NM))
	DNORM1=1.D0/DSQRT(DBLE(NM))
	DKX=2.D0*PI/DBLE(PARA%N)
	DKY=2.D0*PI/DBLE(PARA%M)
	ICOUNT=0
!
!	Reperage cartesien
	II=0
	DO J=1,PARA%M
	   DO I=1,PARA%N
	      II=II+1
	      I_X(II)=I
	      I_Y(II)=J
	   END DO
	END DO
!	
!	Definition de la matrice de transformation ZZ
	DO IK=1,NM
	   KX=I_X(IK)
	   KY=I_Y(IK)
	   KXRF=PARA%N-KX
	   IF(KXRF==0) KXRF=PARA%N
	   KYRF=PARA%M-KY
	   IF(KYRF==0) KYRF=PARA%M
	   IKRF=(KYRF-1)*PARA%N+KXRF
	   IF(IKRF>IK) THEN
	      ICOUNT=ICOUNT+2
	      DO IR=1,NM
	         IX=I_X(IR)
	         IY=I_Y(IR)
	         PHASE=DKX*DBLE(IX*KX)+DKY*DBLE(IY*KY)
	         ZZ(ICOUNT,IR)=DNORM*DCOS(PHASE)
	         ZZ(ICOUNT-1,IR)=DNORM*DSIN(PHASE)
	      END DO
	      QX=DKX*DBLE(KX)
	      QY=DKY*DBLE(KY)
	      ENGY=2.D0*PARA%TPAR*TX*DCOS(QX)+2.D0*PARA%TPER*TY*DCOS(QY)+&
	      	  &4.D0*PARA%TPRM*TX*TY*DCOS(QX)*DCOS(QY)
	      DV(ICOUNT)=DEXP((ENGY+PARA%MU)*DTAU)
	      DV(ICOUNT-1)=DEXP((ENGY+PARA%MU)*DTAU)
	   ENDIF
	   IF(IKRF==IK) THEN
	      ICOUNT=ICOUNT+1
	      DO IR=1,NM
	         IX=I_X(IR)
	         IY=I_Y(IR)
	         PHASE=DKX*DBLE(IX*KX)+DKY*DBLE(IY*KY)
	         ZZ(ICOUNT,IR)=DNORM1*DCOS(PHASE)
	      END DO
	      QX=DKX*DBLE(KX)
	      QY=DKY*DBLE(KY)
	      ENGY=2.D0*PARA%TPAR*TX*DCOS(QX)+2.D0*PARA%TPER*TY*DCOS(QY)+&
	          &4.D0*PARA%TPRM*TX*TY*DCOS(QX)*DCOS(QY)
	      DV(ICOUNT)=DEXP((ENGY+PARA%MU)*DTAU)
	   ENDIF
	END DO
!
!	Calcul de H = exp(K) = ZZT*DV*ZZ
	CALL MAC10(NP,NM,ZZ,H)					! H<=ZZ
	CALL MAC11(NP,NM,H)					! H<=TRANSPOSE(ZZ)
	CALL MULD(NP,NM,2,H,DV)					! H<=H*DV
	CALL MAC2(NP,NM,H,ZZ)					! H<=H*ZZ	
!	
!	Calcul de HM = exp(-K) = ZZT*DV^-1*ZZ
	DO I=1,NM
	   DVN(I)=1/DV(I)
	END DO
	CALL MAC10(NP,NM,ZZ,HM)
	CALL MAC11(NP,NM,HM)
	CALL MULD(NP,NM,2,HM,DVN)
	CALL MAC2(NP,NM,HM,ZZ)	
!
!	Calcul de H2P = exp(K/2) = ZZT*SQRT(DV)*ZZ
	DO I=1,NM
	   DVN(I)=DSQRT(DV(I))
	END DO
	CALL MAC10(NP,NM,ZZ,H2P)
	CALL MAC11(NP,NM,H2P)
	CALL MULD(NP,NM,2,H2P,DVN)
	CALL MAC2(NP,NM,H2P,ZZ)
!
!	Calcul de H2M = exp(-K/2) = ZZT*SQRT(DV)^-1*ZZ	
	DO I=1,NM
	   DVN(I)=1/DSQRT(DV(I))
	END DO
	CALL MAC10(NP,NM,ZZ,H2M)
	CALL MAC11(NP,NM,H2M)
	CALL MULD(NP,NM,2,H2M,DVN)
	CALL MAC2(NP,NM,H2M,ZZ)
!	
	END SUBROUTINE INITH3
!--------------------------------------------------------------------------------------
! Nom		: INITV	[SUB]							15
! Fonction	: Initialisation des registres VUP,VDWN et SPIN (variables du champ HS)
! Auteurs	: AM.Tremblay,A.Reid,H.Nelisse,C.Bourbonnais,C.Boily,L.Chen,H.Touchette
! Revision	: 22 juillet 1997 (Touchette)
! Notes		: - ZETA(2)='y' ou 'Y' : Initialisation a partir d'une configuration
!		    archivee dans FICHSIN sinon, initialisation aleatoire.
!--------------------------------------------------------------------------------------
	SUBROUTINE INITV(PARA,DSEED,CTEU,CTED,SPIN,VUP,VDN)
	USE CONSTANTES,ONLY: NP,NL,NCA
	USE FICHIERS,ONLY: FICHSIN
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER						:: I,J,NM,INUM	! Var. locales
	REAL(KIND(1.D0))				:: DTAU,DEX	! Var. locales
	REAL(KIND(1.D0))				:: ASQT,XLAMBDA	! Var. locales
	REAL(KIND(1.D0))				:: GGUBFS
	REAL(KIND(1.D0)),INTENT(OUT)			:: CTEU,CTED
	REAL(KIND(1.D0)),DIMENSION (NP,NL),INTENT(OUT)	:: VUP,VDN
	INTEGER,DIMENSION (NP,NL),INTENT(OUT)		:: SPIN
	INTEGER(SELECTED_INT_KIND(9)),INTENT(INOUT)	:: DSEED
	TYPE(PARAMETRE),INTENT(IN)			:: PARA
!
!	Initialisation de la configuration
	NM=PARA%N*PARA%M
	SELECT CASE(PARA%ZETA(2))
	   CASE('y','Y')				! Configuration archivee
	      J=LEN_TRIM(FICHSIN)
	      OPEN(UNIT=16,FILE=FICHSIN(1:J),ERR=65)
	      REWIND(UNIT=16)
	      READ(16,*) INUM
	      READ(16,*) INUM
	      DO I=1,PARA%TRAN					
	         READ(16,*)(SPIN(J,I),J=1,NM)			
	      END DO
	      CLOSE(UNIT=16)
	   CASE DEFAULT					! Configuration aleatoire
	      DO I=1,PARA%TRAN
	         DO J=1,NM
	            IF(GGUBFS(DSEED).GE.0.5) THEN
	               SPIN(J,I)=1
	            ELSE
	               SPIN(J,I)=-1
	            ENDIF
	         END DO
	      END DO
	END SELECT
!
!	Initialisation des registres VUP et VDN
	DTAU=PARA%BETA/DBLE(PARA%TRAN)
	DEX=DEXP(0.5D0*DTAU*DABS(PARA%U))
	ASQT=DSQRT(DEX*DEX-1.D0)
	XLAMBDA=DLOG(DEX+ASQT)
	CTEU=DEXP(XLAMBDA)
	CTED=DEXP(-XLAMBDA)
	SELECT CASE(PARA%U.GT.0)
	   CASE(.TRUE.)						! U > 0						
	      DO I=1,PARA%TRAN
	         DO J=1,NM
	            IF(SPIN(J,I).EQ.1) THEN			! SPIN = 1
	               VUP(J,I)=CTEU
	               VDN(J,I)=CTED
	            ELSE					! SPIN =-1
	               VUP(J,I)=CTED
	               VDN(J,I)=CTEU
	            ENDIF
	         END DO
	      END DO
	   CASE(.FALSE.)					! U < 0
	      DO I=1,PARA%TRAN
	         DO J=1,NM
	            IF(SPIN(J,I).EQ.1) THEN			! SPIN = 1
	               VUP(J,I)=CTEU
	               VDN(J,I)=CTEU
	            ELSE					! SPIN =-1
	               VUP(J,I)=CTED
	               VDN(J,I)=CTED
	            ENDIF
	         END DO
	      END DO
	END SELECT
	RETURN
65	WRITE(*,*) '--- Erreur : Fichier inexistant (INITV) ---'
	STOP
	END SUBROUTINE INITV	
!--------------------------------------------------------------------------------------
! Nom		: MAC4251	[SUB]						16
! Fonction	: Multiplication par la gauche, la droite ou les deux (enroulement
!	 	  de la fonction de Green) au 1er ordre de la decomposition de Trotter. 
! Auteurs	: H.Touchette
! Revision	: 27 octobre 1996 (Touchette)
! Notes		: Utilise la forme dense de H sous sa forme factorisee H = ZZT DV ZZ
!		  (forme de Liang Chen)
!--------------------------------------------------------------------------------------
	SUBROUTINE MAC4251(NM,IST,G,V)
	USE CONSTANTES,ONLY: NP
	USE HOPPING
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NM,IST
	INTEGER						:: I
	REAL(KIND(1.D0)),DIMENSION (NP),INTENT(IN)	:: V
	REAL(KIND(1.D0)),DIMENSION (NP)			:: VN		! Variable locale	
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(INOUT):: G	
!
!	Produit par la droite - IST = 2, 3 ou 4 : G'= G*(H*V)^-1
	IF(IST.GT.1) THEN
	   DO I=1,NM
	      VN(I)=1/V(I)					! VN=exp(-V)
	   END DO
	   IF(IST.GT.3)THEN
	      CALL MAC1(NP,NM,HM,G)				! G<=exp(-K)*G
	      CALL MULD(NP,NM,0,G,VN)				! G<=VN*G
	   ELSE
	      CALL MULD(NP,NM,2,G,VN)				! G<=G*VN
	      CALL MAC2(NP,NM,G,HM)				! G<=G*exp(-K)
	   ENDIF
        ENDIF
!
!	Produit par la gauche - IST = 0, 1 ou 2 : G'= (H*V)*G
        IF(IST.LT.3.OR.IST.GT.4) THEN
	   IF(IST.LT.1.OR.IST.GT.4)THEN
	      CALL MAC2(NP,NM,G,H)				! G<=G*exp(K)
	      CALL MULD(NP,NM,2,G,V)				! G<=G*V
	   ELSE
	      CALL MULD(NP,NM,0,G,V)				! G<=V*G
	      CALL MAC1(NP,NM,H,G)				! G<=exp(K)*G
	   ENDIF
        ENDIF
!		
	END SUBROUTINE MAC4251
!--------------------------------------------------------------------------------------
! Nom		: INITW 	[SUB]						17
! Fonction	: Calcul des coefficients pour le l'integration 
! Auteurs	:
! Revision	: 1 decembre 1996 (Touchette)
! Notes		: - Simpson (integrale a 3 points) : ! S
!		  - Trapeze (integrale a 2 points) : ! T 
!		  - Noter la borne TRAN+1
!--------------------------------------------------------------------------------------	
	SUBROUTINE INITW(NTRAN)
	USE UTILITAIRE,ONLY: W
	IMPLICIT NONE
	INTEGER,INTENT(IN)		:: NTRAN
	INTEGER				:: I			! Variable locale
	W(0)=1.D0/3.D0						! S
!	W(0)=1.D0/2.D0						! T
	DO I=1,NTRAN-1
	   W(I)=4.D0/3.D0-2.D0/3.D0*MOD(I+1,2)			! S
!	   W(I)=1.D0						! T
	END DO
	W(NTRAN)=1.D0/3.D0					! S
!	W(NTRAN)=1.D0/2.D0					! T					
	END SUBROUTINE INITW
!--------------------------------------------------------------------------------------
! Nom		: LUDCMP 	[SUB]						18
! Fonction	: Decomposition gaussienne A = LU
! Auteurs	:
! Revision	:
! Notes		: Provient entierement de Numerical Recipes
!--------------------------------------------------------------------------------------	
      	SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      	IMPLICIT REAL(KIND(1.D0)) (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.) THEN
	      WRITE(*,*) '--- Erreur : Singular matrix LUDCMP'
	      STOP
	   ENDIF
           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		: INVMG 	[SUB]						19-(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(KIND(1.D0))				:: D		! Var. locale
	INTEGER						:: I,J		! var. locale
	INTEGER,INTENT(IN)				:: NP,NM
	INTEGER,DIMENSION (NP)				:: INDX   	! Var. locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: Y		! Var. locale
	REAL(KIND(1.D0)),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
	CALL MAC10(NP,NM,Y,A)						
	END SUBROUTINE INVMG
!--------------------------------------------------------------------------------------
! Nom		: LUBKSB 	[SUB]						19-(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(KIND(1.D0)) A(NP,NP),B(N)
	INTEGER I,II,J,LL
	REAL(KIND(1.D0)) 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		: GGUBFS 	[FUNCTION]					20
! Fonction	: Generateur de nombres pseudo-aleatoires
! Auteurs	: 
! Revision	: 15 decembre 1996
! References	: [1] Numerical Recipes, Fortran 90, p.1142.
! Notes		: Remplace le generateur de P. L'Ecuyer.
!
!--------------------------------------------------------------------------------------
	FUNCTION GGUBFS(IDUM)
	IMPLICIT NONE
	INTEGER,PARAMETER 		:: K4B=SELECTED_INT_KIND(9)
	INTEGER(K4B),INTENT(INOUT) 	:: IDUM
	REAL(KIND(1.D0)) 		:: GGUBFS
	INTEGER(K4B),PARAMETER 		:: IA=16807,IM=2147483647,IQ=127773,IR=2836
	REAL(KIND(1.D0)),SAVE 		:: AM
	INTEGER(K4B),SAVE 		:: IX=-1,IY=-1,K
!
	IF(IDUM<=0.OR.IY<0) THEN
	   AM=NEAREST(1.0,-1.0)/IM
	   IY=IOR(IEOR(888889999,ABS(IDUM)),1)
	   IX=IEOR(777755555,ABS(IDUM))
	   IDUM=ABS(IDUM)+1
	ENDIF
	IX=IEOR(IX,ISHFT(IX,13))
	IX=IEOR(IX,ISHFT(IX,-17))
	IX=IEOR(IX,ISHFT(IX,5))
	K=IY/IQ
	IY=IA*(IY-K*IQ)-IR*K
	IF(IY<0) IY=IY+IM
	GGUBFS=AM*IOR(IAND(IM,IEOR(IX,IY)),1)
	END FUNCTION GGUBFS
!--------------------------------------------------------------------------------------
! Nom		: INIT_LOC 	[SUB]						21
! Fonction	: Initialisation des vecteurs de localisation des sites
! Auteurs	:
! Revision	: 22 octobre 1996 (Touchette)
! Notes		: Le site I se trouve sur CHI(I) et COLI(I)
!--------------------------------------------------------------------------------------
	SUBROUTINE INIT_LOC(N,M)
	USE UTILITAIRE,ONLY: CHI,COLI
	IMPLICIT NONE	
	INTEGER,INTENT(IN)		:: N,M	
	INTEGER				:: I,NM			! Variables locales
	NM=N*M
	DO I=1,NM
	   CHI(I)=(I+N-1)/N
	   COLI(I)=I-(CHI(I)-1)*N
	END DO
	END SUBROUTINE INIT_LOC
!--------------------------------------------------------------------------------------
! Nom		: DEFCOSMAG	[SUB]						22
! Fonction	: Construction de la matrice des cosinus pre-calcules
! Auteurs	: A.Veilleux,H.Touchette
! Revision	: 20 fevrier 1996 (Touchette)
! Notes		: Utilise par MAGST,PSCOR. Capacite fixee par ND1 et ND2
!--------------------------------------------------------------------------------------
	SUBROUTINE DEFCOSMAG(N,M)
	USE UTILITAIRE,ONLY: CHI,COLI,COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))		:: PHI,PI		! Variables locales
	INTEGER,INTENT(IN)		:: N,M
	INTEGER				:: KX,KY,I,J,NX,NY	! Variables locales
	INTEGER				:: NM,NCHI,NCHJ		! Variables locales
	INTEGER				:: NCOLJ,NCOLI		! Variables locales
!
	NM=N*M
	COSMAG(1:N,1:M,1:NM,1:NM)=0.D0				! Boucle implicite
	PI=DACOS(-1.D0)
	DO KX=1,N
	   DO KY=1,M
	     DO I=1,NM
	        NCHI=CHI(I)
	        NCOLI=COLI(I)
	        DO J=1,NM
	           NCHJ=CHI(J)
	           NCOLJ=COLI(J)
	           NX=NCOLI-NCOLJ
	           NY=NCHI-NCHJ
	           PHI=DCOS(2.D0*PI*KX/N*NX+2.D0*PI*KY/M*NY)
	           COSMAG(KX,KY,I,J)=PHI
	        END DO
	     END DO
	   END DO
	END DO
	END SUBROUTINE DEFCOSMAG
!--------------------------------------------------------------------------------------
! Nom		: AUTO_CORRELATION	[SUB]					23
! Fonction	: Calcul de la fonction d'autocorrelation
! Auteurs	: H.Touchette
! Revision	: 27 novembre 1996 (Touchette)
! References	: [1] Negele & Orland, Quantum Many-Particle Systems, Addison-Wesley,
!		      p.445.
!		  [2] Lacaze et al., cond-mat/9610193 28 octobre 1996.
!		  [3] Koonin & Meredith, Computational Physics, Addison-Wesley, 1989.
! Notes		: - Calcul seulement a la fin d'une simulation
!		  - Les donnees dans AUTOCORR sont remplacees a la sortie par la
!		    fonction d'autocorrelation 
!--------------------------------------------------------------------------------------
	SUBROUTINE AUTO_CORRELATION(C1,AUTOCORR)
	USE CONSTANTES,ONLY: NCOR
	IMPLICIT NONE	
	INTEGER						:: MC1,I,K		! Var locales
	REAL(KIND(1.D0))				:: SUM,MOY1,MOY2	! var locales
	REAL(KIND(1.D0)),INTENT(IN)			:: C1
	REAL(KIND(1.D0)),DIMENSION (NCOR)		:: TEMP			! Var locale
	REAL(KIND(1.D0)),DIMENSION (NCOR),INTENT(INOUT)	:: AUTOCORR
!
	MC1=INT(C1)
	MOY1=0.D0
	MOY2=0.D0
	DO I=1,MC1
	   MOY1=MOY1+AUTOCORR(I)
	   MOY2=MOY2+AUTOCORR(I)*AUTOCORR(I)
	END DO
	MOY1=(MOY1/C1)**2
	MOY2=(MOY2/C1)-MOY1					
	DO K=0,MC1-1
	   SUM=0.D0
	   DO I=1,MC1-K
	      SUM=SUM+AUTOCORR(I)*AUTOCORR(I+K)
	   END DO
	   SUM=SUM/DBLE(MC1-K)						
	   TEMP(K+1)=(SUM-MOY1)/MOY2
	END DO	
	AUTOCORR(1:MC1)=TEMP(1:MC1)
	END SUBROUTINE AUTO_CORRELATION
!--------------------------------------------------------------------------------------
! Nom		: IDENTIFICATION	[SUB]					24
! Fonction	: Initialisation des parametres d'identification des fichiers
! Auteurs	: H.Touchette
! Revision	: 15 juillet 1997 (Touchette)
! References	: [1] AIX XL Fortran Compiler/6000, Language Reference, p.33-
! Notes		: 
!--------------------------------------------------------------------------------------
	SUBROUTINE IDENTIFICATION(PARA,DSEEDI)
	USE FICHIERS,ONLY: FICPAOUT,FICHSOUT,FICGIOUT,FICCVOUT,REPWEB
	USE IDENT
	USE PARA_SIM
	IMPLICIT NONE
	CHARACTER(LEN=5)				:: NUMTOASCII
	CHARACTER(LEN=NCA)				:: TP1,TP2,TP3,TP4,TP5,TPI	! Loc.
	INTEGER						:: I1,I2,I3,I4,I5,I6,I7,II	! Loc.
	INTEGER(SELECTED_INT_KIND(9)),INTENT(IN)	:: DSEEDI
	TYPE(PARAMETRE),INTENT(IN)			:: PARA
!
	IDNUM=ABS(DSEEDI)
	TP1=NUMTOASCII(INT(DABS(PARA%U)))
	TP2=NUMTOASCII(INT(ANINT(100*(PARA%REMPLD))))
	TP3=NUMTOASCII(PARA%N)
	TP4=NUMTOASCII(PARA%M)
	TP5=NUMTOASCII(INT(ANINT(PARA%BETA)))
	TPI=NUMTOASCII(IDNUM)
	I1=INDEX(TP1,' ')
	I2=INDEX(TP2,' ')
	I3=INDEX(TP3,' ')
	I4=INDEX(TP4,' ')
	I5=INDEX(TP5,' ')
	II=INDEX(TPI,' ')
	I6=INDEX(PREFIX,' ')
	I7=INDEX(SUFFIX,' ')
	NANUM=PREFIX(1:I6-1)//TP1(1:I1-1)//TP2(1:I2-1)//TP3(1:I3-1)//TP4(1:I4-1)//TP5(1:I5-1)
	I6=INDEX(NANUM,' ')
	FICPAOUT=TPI(1:II-1)//'.'//NANUM(1:I6-1)//'re'//'.'//SUFFIX(1:I7-1)//' '
	FICHSOUT=TPI(1:II-1)//'.'//NANUM(1:I6-1)//'hs'//'.'//SUFFIX(1:I7-1)//' '
	FICGIOUT=TPI(1:II-1)//'.'//NANUM(1:I6-1)//'gi'//'.'//SUFFIX(1:I7-1)//' '
	FICCVOUT=TPI(1:II-1)//'.'//NANUM(1:I6-1)//'cv'//'.'//SUFFIX(1:I7-1)//' '
	END SUBROUTINE IDENTIFICATION
!--------------------------------------------------------------------------------------
! Nom		: GLOBAL_MOVE	[SUB]						26
! Fonction	: Basculement des spins d'un site sur toutes les tranches de temps
! Auteurs	: H.Touchette
! Revision	: 6 mai 1997 (Touchette)
! References	: [1] Scalettar et al.,PRB,44,10502,1991.
! Notes		: 
!--------------------------------------------------------------------------------------
	SUBROUTINE GLOBAL_MOVE(NTRAN,CUP,CDN,DU,SPN,VUP,VDN,NCH)
	USE CONSTANTES,ONLY: NP,NL
	INTEGER						:: I		! Var. locale
	INTEGER,INTENT(IN)				:: NTRAN,NCH		
	REAL(KIND(1.D0)),INTENT(IN)			:: CUP,CDN,DU
	INTEGER,DIMENSION(NP,NL),INTENT(INOUT)		:: SPN
	REAL(KIND(1.D0)),DIMENSION(NP,NL),INTENT(INOUT)	:: VUP,VDN
!
	IF(DU.GT.0) THEN
	   DO I=1,NTRAN 
	      IF(SPN(NCH,I).EQ.1) THEN
	         VUP(NCH,I)=CDN
	         VDN(NCH,I)=CUP
	         SPN(NCH,I)=-1
	      ELSE
	         VUP(NCH,I)=CUP
	         VDN(NCH,I)=CDN
	         SPN(NCH,I)=1
	      ENDIF
	   END DO
	ELSE
	   DO I=1,NTRAN 
	      IF(SPN(NCH,I).EQ.1) THEN
	         VUP(NCH,I)=CDN
	         VDN(NCH,I)=CDN
	         SPN(NCH,I)=-1
	      ELSE
	         VUP(NCH,I)=CUP
	         VDN(NCH,I)=CUP
	         SPN(NCH,I)=1
	      ENDIF
	   END DO
	ENDIF
	END SUBROUTINE GLOBAL_MOVE
!--------------------------------------------------------------------------------------
!			    	  ROUTINES - TEMPS EGAL
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! Nom		: GCALS 	[SUB]						ME1
! Fonction	: Calcul de la fonction de Green a temps egal
! Auteurs	: AM.Tremblay,A.Reid,H.Nelisse,C.Bourbonnais,C.Boily,H.Touchette
! Revision	: 28 juillet 1997 (Touchette)
! Reference 	: H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		  hiver 96 - ete 97.
! Notes		: - Avec procedure de stabilisation (orthogonalisation + pivot)
!		  - Calcul optionnel du determinant de G^-1 
!--------------------------------------------------------------------------------------
	SUBROUTINE GCALS(NM,L,PARA,V,G,DET)
	USE CONSTANTES,ONLY: NP,NL
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER						:: I,J,K	! Var. locales
	REAL(KIND(1.D0))				:: DETD,DETR 	! Var. locales
	INTEGER,INTENT(IN)				:: NM,L
	REAL(KIND(1.D0)),INTENT(OUT)			:: DET
	INTEGER,DIMENSION (NP)				:: INDX		! Var. locale
	REAL(KIND(1.D0)),DIMENSION (NP)			:: D		! Var. locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: U,R,R2,R3	! Var. locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT)	:: G
	REAL(KIND(1.D0)),DIMENSION (NP,NL),INTENT(IN)	:: V
	TYPE(PARAMETRE),INTENT(IN)			:: PARA
!
!	Initialisation des matrices
	DO J=1,NM
	   DO I=1,NM
	      U(I,J)=0.D0							
	      R(I,J)=0.D0							
	   END DO
	   U(J,J)=1.D0
	   R(J,J)=1.D0
	   D(J)=1.D0
	END DO
!
!	Boucle de L a TRAN
	K=L
	DO WHILE(K.LT.PARA%TRAN+1)
	   DO I=1,PARA%NMODUL					! U=B(K+NMODUL)...B(K)
	      CALL MAC4251(NM,1,U,V(:,K))						
	      K=K+1
	   END DO
	   CALL MULD(NP,NM,2,U,D)				! U = U*D
	   CALL GRASCHP(NP,NM,U,D,R2)				! UDV (V=R2) 
	   CALL MAC1(NP,NM,R2,R)				! V finale = R
	END DO
!
!	Boucle de 1 a L-1 
	K=1
	DO WHILE(K.LT.L)			
	   DO I=1,PARA%NMODUL
	      CALL MAC4251(NM,1,U,V(:,K))
	      K=K+1
	   END DO
	   CALL MULD(NP,NM,2,U,D)
	   CALL GRASCHP(NP,NM,U,D,R2)
	   CALL MAC1(NP,NM,R2,R)
	END DO
!
!	Decomposition UDV du produit des matrices B 
!	                 -1  -1 -1  -1 -1
!	Calcul de   G = V  (U  V +D)  U
!
!	           -1 -1
!	Calcul de U  V  +D
!
	CALL MAC10(NP,NM,R,R3)					! V													
	CALL INVMG(NP,NM,R)					! V^-1					
 	CALL MAC10(NP,NM,U,G)													         
	CALL MAC11(NP,NM,U)					! U <= U^T = U^-1					
	CALL MAC2(NP,NM,U,R)					! U <= U^-1 V^-1	
	DO I=1,NM					
	   U(I,I)=U(I,I)+D(I)					! U <= U^-1 V^-1 +D
	END DO
!
!	Decomposition UDV avant inversion				      -1 -1
	CALL GRASCHP(NP,NM,U,D,R2)				! U'D'V' <= (U  V  +D)
!
!	Calcul optionnel du determinant 
	IF((PARA%ZETA(4).EQ.'Y'.OR.PARA%ZETA(4).EQ.'y').AND.L.EQ.1) THEN
	   DETD=1.D0
	   DO I=1,NM
	      DETD=DETD*D(I)
	   END DO
	   CALL MAC1(NP,NM,R2,R3)				
	   CALL LUDCMP(R3,NM,NP,INDX,DETR)
	   DO I=1,NM
	      DETR=DETR*R3(I,I)
	   END DO
	   DET=DETR*DETD
	ENDIF
!
!                  -1  -1 -1      -1
!	Inversion V  R2  D  (U3 U)						   
	CALL INVMG(NP,NM,R2)					! V'^-1
	CALL MAC1(NP,NM,R,R2)					! V^-1 * V'^-1						
	DO J=1,NM						!        								
	   DO I=1,NM						! (V^-1 V'^-1)*D^-1						
	      R2(I,J)=R2(I,J)/D(J)
	   END DO
	END DO
	CALL MAC2(NP,NM,G,U)					! U U'									   
	CALL MAC11(NP,NM,G)					! G <= (U U')^T
	CALL MAC1(NP,NM,R2,G)					! V2*G						
	END SUBROUTINE
!--------------------------------------------------------------------------------------
! Nom		: OCC_REMP		[SUB]					ME2
! Fonction	: Calcul des parametres de la bande (remplissage, mu propose)
! Auteurs	: 
! Revision	: 6 novembre 1996 (Touchette)
! Notes		: 
!--------------------------------------------------------------------------------------
	SUBROUTINE OCC_REMP(NP,NM,SIG,GUP,GDWN,DMU,SOCC,XBETA)
	IMPLICIT NONE
	INTEGER						:: I,J			! Var. locales
	REAL(KIND(1.D0))				:: SOM,SI,SJ,SUP,SDW	! Var. locales
	INTEGER,INTENT(IN)				:: NP,NM
	REAL(KIND(1.D0)),INTENT(INOUT)			:: DMU
	REAL(KIND(1.D0)),INTENT(OUT)			:: SOCC
	REAL(KIND(1.D0)),INTENT(IN)			:: XBETA,SIG
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDWN
!
	SOM=0.D0
	DO I=1,NM
	   SOM=SOM+2.D0-GUP(I,I)-GDWN(I,I)
	END DO
	SOCC=SOM/NM
	SOM=0.D0
	DO I=1,NM
	   SI=2.D0-GUP(I,I)-GDWN(I,I)
	   DO J=1,NM
	      SJ=2.D0-GUP(J,J)-GDWN(J,J)
	      SUP=GUP(J,I)*GUP(I,J)
	      SDW=GDWN(J,I)*GDWN(I,J)
	      SOM=SOM+SI*SJ-SUP-SDW
	   END DO
	   SOM=SOM+GUP(I,I)+GDWN(I,I)
	END DO
	DMU=DMU-SIG*XBETA*SOM/NM
	END SUBROUTINE OCC_REMP
!--------------------------------------------------------------------------------------
! Nom		: ENETOT	[SUB]						ME3
! Fonction	: Calcul de l'energie totale des electrons libres ayant le meme
!		  potentiel chimique mu
! Auteurs	: AM.Tremblay
! Revision	: 27 octobre 1996 (Touchette)
! Notes		:
!--------------------------------------------------------------------------------------
	SUBROUTINE ENE_TOT(PARA,ENERG)
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER					:: I,J,NM		! Var. locales
	REAL(KIND(1.D0)),INTENT(OUT)		:: ENERG		
	REAL(KIND(1.D0)) 			:: KX,KY		! Var. locales
	REAL(KIND(1.D0))			:: PI,EPS,DEN,ENE
	TYPE(PARAMETRE),INTENT(IN)		:: PARA
!
	PI=DACOS(-1.D0)
	NM=PARA%N*PARA%M
	ENERG=0.D0
	   DO I=0,PARA%N-1
	      KX=2.D0*PI*DBLE(I)/DBLE(PARA%N)
	      DO J=0,PARA%M-1
	         KY=2.D0*PI*DBLE(J)/DBLE(PARA%M)
	         EPS=ENE(PARA%N,PARA%M,PARA%TPAR,PARA%TPER,PARA%TPRM,KX,KY)
	         DEN=DEXP(PARA%BETA*(EPS-PARA%MU))+1
	         ENERG=ENERG+EPS/DEN
	      END DO
	   END DO
	ENERG=ENERG*2.D0/DBLE(NM)
	END SUBROUTINE ENE_TOT
!--------------------------------------------------------------------------------------
! Nom		: ENE		[FUNCTION]					ME4
! Fonction	: Relation de dispersion des electrons libres
! Auteurs	: H.Touchette
! Revision	: 15 juillet 1997 (Touchette)
! Notes		: - Tient compte des cas speciaux N,M = 2, de la periodicite du reseau,
!		    de meme que des seconds voisins.
!
!--------------------------------------------------------------------------------------	
	FUNCTION ENE(N,M,TPAR,TPER,TPRM,KX,KY)
	REAL(KIND(1.D0))			:: ENE
	INTEGER,INTENT(IN)			:: N,M
	REAL(KIND(1.D0)),INTENT(IN)		:: TPAR,TPER,TPRM,KX
	REAL(KIND(1.D0)),INTENT(IN)		:: KY
	REAL(KIND(1.D0))			:: TPA,TPE,TP2		! Var. locales
!
	TPA=TPAR
	TPE=TPER
	TP2=TPRM
	IF(N.EQ.2) THEN
	   TPA=TPA/2
	   TP2=TP2/2
	ENDIF
	IF(M.EQ.2) THEN
	   TPE=TPE/2
	   TP2=TP2/2
	ENDIF
	IF(M.EQ.1) THEN
	   ENE=-2.D0*TPA*DCOS(KX)
	ELSE
	   ENE=-2.D0*(TPA*DCOS(KX)+TPE*DCOS(KY))-4*TP2*DCOS(KX)*DCOS(KY)
	ENDIF
	END FUNCTION ENE	 
!--------------------------------------------------------------------------------------
! Nom		: DSITE		[SUB]						ME5
! Fonction	: Calcul de l'energie 
! Auteurs	: AM.Tremblay,A.Reid,C.Boily
! Revision	: 21 juillet 1997 (Touchette)
! Reference	: [1] AM.Tremblay, Code Hirsch,Departement de physique de l'Universite
!		      de Sherbrooke.
! Notes		: - A cause des conditions periodiques, le cas 2 sites est special
!		  - Pas d'extension pour t' !$
!--------------------------------------------------------------------------------------
	SUBROUTINE DSITE(PARA,GUP,GDN,CIN,TOT)
	USE CONSTANTES,ONLY: NP
	USE UTILITAIRE,ONLY: CHI,COLI
	USE PARA_SIM
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: CTE,POT,CIN1,CIN2	! Var. locales
	REAL(KIND(1.D0))				:: TPA,TPE,TP2		! Var. locales
	INTEGER						:: NM,I,II,II2		! Var. locales
	INTEGER						:: NCHI,NCOLI		! Var. locales
	REAL(KIND(1.D0)),INTENT(OUT)			:: CIN,TOT
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDN	
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: BUP,BDN		! Var. locales
	TYPE(PARAMETRE),INTENT(IN)			:: PARA
!
	CTE=-1.D0
	NM=PARA%M*PARA%N
	CALL MAC8(NP,NM,GUP,BUP,CTE)
	CALL MAC8(NP,NM,GDN,BDN,CTE)
	CALL MAC9(NP,NM,BUP)					! BUP = 1 - GUP
	CALL MAC9(NP,NM,BDN)					! BDN = 1 - GDN
!
	CIN=0.D0
	POT=0.D0
	TOT=0.D0
	TPA=PARA%TPAR
	TPE=PARA%TPER
	TP2=PARA%TPRM
	IF(PARA%N.EQ.2) TPA=TPA/2.D0
	IF(PARA%M.EQ.2) TPE=TPE/2.D0
	DO I=1,NM
	   NCHI=CHI(I)
	   NCOLI=COLI(I)
	   II=I+1
	   II2=I-PARA%N
	   IF(NCOLI.EQ.PARA%N) II=I-PARA%N+1
	   IF(NCHI.EQ.1) II2=NM-PARA%N+I
	   CIN1=BUP(II,I)+BUP(I,II)+BDN(II,I)+BDN(I,II)
	   CIN2=BUP(II2,I)+BUP(I,II2)+BDN(II2,I)+BDN(I,II2)
	   IF(PARA%M.EQ.1) CIN2=0.D0
	   CIN=CIN+TPA*CIN1+TPE*CIN2
	   POT=POT+PARA%U*BUP(I,I)*BDN(I,I)
        END DO
	CIN=-CIN
	TOT=POT+CIN
	END SUBROUTINE DSITE
!--------------------------------------------------------------------------------------
! Nom		: DENSITE		[SUB]					ME6
! Fonction	: Calcul de l'impulsion 
! Auteurs	: AM.Tremblay,A.Reid,C.Boily
! Revision	: 21 juillet 1997 (Touchette)
! Reference	: [1] AM.Tremblay, Code Hirsch,Departement de physique de l'Universite
!		      de Sherbrooke.
! Notes		: 
!--------------------------------------------------------------------------------------
	SUBROUTINE DENSITE(N,M,GUP,GDN,X,NV,V)
	USE CONSTANTES,ONLY: NP,ND1,ND2
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: SUM,CTE		! Var. locales
	REAL(KIND(1.D0))				:: GSU,GSD		! Var. locales
	INTEGER						:: K,I,J,KX,KY,NM		! Var. locales
	INTEGER,INTENT(IN)				:: N,M,NV
	INTEGER,DIMENSION(NP,2)				:: V
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDN
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: BUP,BDN		! Var. locales
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(OUT):: X
!	
	NM=N*M
	CTE=-1.D0
	CALL MAC8(NP,NM,GUP,BUP,CTE)
	CALL MAC8(NP,NM,GDN,BDN,CTE)
	CALL MAC9(NP,NM,BUP)					! BUP = 1 - GUP
	CALL MAC9(NP,NM,BDN)					! BDN = 1 - GDN
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   SUM=0.D0
	   DO I=1,NM
	      DO J=1,NM
	         GSU=BUP(I,J)+BUP(J,I)
	         GSD=BDN(I,J)+BDN(J,I)
	         SUM=SUM+COSMAG(KX,KY,I,J)*(GSU+GSD)
	      END DO
	   END DO
	   X(KX,KY)=SUM/NM
	END DO   
	END SUBROUTINE DENSITE
!--------------------------------------------------------------------------------------
! Nom		: PSCOR		[SUB]						ME7
! Fonction	: Calcul de Ps(Qx,Qy) : s-wave pair-field correlation function 
!		  (facteur de structure de paire)
! Auteurs	: AM.Tremblay,H.Touchette (1996)
! Revision	: Fevrier 1996 (Touchette)
! Reference	: [1] A.Moreo,D.J.Scalapino,PRB 45,7544,1992.
!		  [2] H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		      hiver 96 - ete 97.
! Notes		: 
!--------------------------------------------------------------------------------------
	SUBROUTINE PSCOR(N,M,GUP,GDN,X,NV,V)
	USE CONSTANTES,ONLY: NP,ND1,ND2
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: GIJ,SUM		! Var. locales
	REAL(KIND(1.D0))				:: CTE,PHI		! Var. locales
	INTEGER						:: K,KX,KY,I,J		! Var. locales
	INTEGER						:: NM			! Var. locales
	INTEGER,INTENT(IN)				:: N,M,NV
	INTEGER,DIMENSION(NP,2)				:: V
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDN
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: AUP,ADW		! Var. locales
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(OUT)	:: X
!
	NM=N*M
	CTE=-1.D0
	CALL MAC8(NP,NM,GUP,AUP,CTE)
	CALL MAC8(NP,NM,GDN,ADW,CTE)
	CALL MAC9(NP,NM,AUP)					! AUP = 1 - GUP
	CALL MAC9(NP,NM,ADW)					! ADW = 1 - GDWN
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   SUM=0.D0
	   DO I=1,NM
	      DO J=1,NM
	         PHI=COSMAG(KX,KY,I,J)
	         GIJ=AUP(J,I)*ADW(J,I)+GUP(I,J)*GDN(I,J)
	         SUM=SUM+PHI*GIJ
	      END DO
	   END DO
	   X(KX,KY)=SUM/NM
	END DO
	END SUBROUTINE PSCOR
!--------------------------------------------------------------------------------------
! Nom		: MAGST		[SUB]						ME8
! Fonction	: Calcul du facteur de structure magnetique 
! Auteurs	: AM.Tremblay,A.Reid,H.Nelisse,C.Bourbonnais,C.Boily
! Revision	: 6 fevrier 1996 (Touchette)
! References	: [1] J.E.Hirsch,PRB 31,4409,1985.
!		  [2] Loh,Gubernatis pour le meilleur estimateur.
! Notes		: - Choix de l'estimateur [S+,S-] et [Sz], On active l'estimateur
! 	 	    voulu en enlevant ses marques de commentaires et en placant
! 	 	    des marques de commentaires pour l'autre estimateur.
!
!--------------------------------------------------------------------------------------
	SUBROUTINE MAGST(N,M,GUP,GDN,X,NV,V)
	USE CONSTANTES,ONLY: NP,ND1,ND2
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: PHI,CTE,SUM	! Variables locales
	REAL(KIND(1.D0))				:: GI,GJ,GIJ	! Variables locales
	INTEGER						:: L,I,J,KX,KY,NM	! Variables locales
	INTEGER,INTENT(IN)				:: N,M,NV
	INTEGER,DIMENSION(NP,2)				:: V
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDN
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: AUP,ADW 	! Variables locales
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(OUT):: X	
!
	NM=N*M
	CTE=-1.D0
        CALL MAC8(NP,NM,GUP,AUP,CTE)
	CALL MAC8(NP,NM,GDN,ADW,CTE)
        CALL MAC9(NP,NM,AUP)					! AUP = 1 - GUP	
        CALL MAC9(NP,NM,ADW) 					! ADW = 1 - GDWN  
!
	DO L=1,NV
	   KX=V(L,1)
	   KY=V(L,2)
	   SUM=0.D0
	   DO I=1,NM
	      ! Utilise pour l'estimateur Sz (et non pour S+,S-)
              GI=AUP(I,I)-ADW(I,I)
	      DO J=1,NM
	         ! Calcul de l'estimateur Sz (les 2 lignes qui suivent)
                 GJ=AUP(J,J)-ADW(J,J)
                 GIJ=GI*GJ+AUP(J,I)*GUP(I,J)+ADW(J,I)*GDN(I,J)
	         ! Calcul de l'estimateur S+,S-
	         ! GIJ=AUP(J,I)*GDN(I,J)+ADW(J,I)*GUP(I,J)
		 PHI=COSMAG(KX,KY,I,J)
		 SUM=SUM+PHI*GIJ
	      END DO
	   END DO
	   X(KX,KY)=SUM/NM
	END DO
	END SUBROUTINE MAGST
!--------------------------------------------------------------------------------------
! Nom		: RHOE		[SUB]						ME9
! Fonction	: Susceptibilite de charge a temps egal
! Auteurs	: H.Touchette
! Revision	: 16 juillet 1997 (Touchette)
! References	: H.Touchette, Aspects numeriques des simulation Monte Carlo quantique 
!		  du modele de Hubbard, rapport technique du CRPS, departement de 
!		  physique, UdeS, ete 1997.
! Notes		: - C(Qx,Qy)
!--------------------------------------------------------------------------------------
	SUBROUTINE RHOE(N,M,GUP,GDN,FST,NV,V)
	USE CONSTANTES,ONLY: NP,NL,ND1,ND2
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	INTEGER						:: K,I,J,KX,KY,NM		! Var. locales
	REAL(KIND(1.D0))				:: G1,G2,CTE		! Var. locales
	REAL(KIND(1.D0))				:: PROD,PHI,ACC		! Var. locales
	INTEGER,INTENT(IN)				:: N,M,NV
	INTEGER,DIMENSION(NP,2)				:: V
	REAL(KIND(1.D0)),DIMENSION (NP,NP)		:: AUP,ADN		! Var. locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: GUP,GDN
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(OUT):: FST
!
	NM=N*M
	CTE=-1.D0
        CALL MAC8(NP,NM,GUP,AUP,CTE)
	CALL MAC8(NP,NM,GDN,ADN,CTE)
        CALL MAC9(NP,NM,AUP)					! AUP = 1 - GUP	
        CALL MAC9(NP,NM,ADN) 					! ADW = 1 - GDWN 
! 	
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   ACC=0.D0
	   DO I=1,NM
	      G1=AUP(I,I)+ADN(I,I)
	      DO J=1,NM
	         G2=AUP(J,J)+ADN(J,J)
	         PROD=G1*G2+AUP(J,I)*GUP(I,J)+ADN(J,I)*GDN(I,J)
	         PHI=COSMAG(KX,KY,I,J)
	         ACC=ACC+PROD*PHI
	      END DO
	   END DO
	   FST(KX,KY)=ACC/NM
	END DO
	END SUBROUTINE RHOE
!--------------------------------------------------------------------------------------
!		 	   	ROUTINES - TEMPS INEGAL
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! Nom		: GICALS	[SUB]						MI1
! Fonction	: Calcul des fonctions de Green a temps inegal au second ordre
!		  avec procedure de stabilisation
! Auteurs	: H.Touchette (1996)
! Revision	: 30 juillet 1997 (Touchette)
! Reference	: [1] H.Touchette, Aspects numeriques des simulations Monte Carlo  
!		      quantique du modele de Hubbard, rapport technique du CRPS,  
!		      departement de physique, UdeS, ete 1997.
!		  [2] H.G.Matuttis,Ph.D. Thesis,University of Regensburg,1995.
! Notes		: - En entree, on recoit G(l=1,l=1) au 1er ordre mais on calcule
!		    G(tau) au 2e ordre par MACGD
!		  
!--------------------------------------------------------------------------------------
	SUBROUTINE GICALS(NM,PARA,G,V,GI,L,SPIN)
	USE CONSTANTES,ONLY: NP,NL
	USE PARA_SIM
	IMPLICIT NONE
	REAL(KIND(1.D0)),SAVE				:: CTE
	INTEGER						:: J,K,NTOP	! Var. locales
	INTEGER,INTENT(IN)				:: NM
	REAL(KIND(1.D0)),DIMENSION (NP,2),SAVE		:: D		! Var. locale
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: G
	REAL(KIND(1.D0)),DIMENSION (NP,NP)              :: GP       ! Var. locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP,2),SAVE	:: AP
	REAL(KIND(1.D0)),DIMENSION (NP,NP,2),SAVE	:: U,R,R2	! Var. locales 
	REAL(KIND(1.D0)),DIMENSION (NP,NL),INTENT(IN)	:: V
	REAL(KIND(1.D0)),DIMENSION(NP,NP),INTENT(OUT)	:: GI
	TYPE(PARAMETRE),INTENT(IN)			:: PARA
	INTEGER						:: L,SPIN
!
	IF(L==0.OR.L==PARA%TRAN) THEN
	   CTE=-1.D0
	   IF(L==0) THEN
	      CALL MACGD(NP,NM,G,GI(:,:))
	   ELSE
	      CALL MAC8(NP,NM,G,AP(:,:,SPIN),CTE)				! AP<=-G
	      CALL MAC9(NP,NM,AP(:,:,SPIN))					! AP<=1-G
	      CALL MACGD(NP,NM,AP(:,:,SPIN),GI(:,:))
	   END IF
!
	   DO J=1,NM
	      DO K=1,NM
	         U(K,J,SPIN)=0.D0
	         R(K,J,SPIN)=0.D0
	      END DO
	      U(J,J,SPIN)=1.D0
	      R(J,J,SPIN)=1.D0
	      D(J,SPIN)=1.D0
	   END DO
	ELSE
!
	   IF(L.LE.PARA%TRAN/2) THEN
	      CALL MAC4251(NM,1,U(:,:,SPIN),V(:,L))
	      NTOP=MOD(L,PARA%NMODUL)
	      IF(NTOP.EQ.0) THEN
	         CALL MULD(NP,NM,2,U(:,:,SPIN),D(:,SPIN))
	         CALL GRASCHP(NP,NM,U(:,:,SPIN),D(:,SPIN),R2(:,:,SPIN))
	         CALL MAC1(NP,NM,R2(:,:,SPIN),R(:,:,SPIN))
	      ENDIF
!             
	      CALL MAC10(NP,NM,G,GP)	
	      CALL MAC1(NP,NM,R(:,:,SPIN),GP)
	      CALL MULD(NP,NM,0,GP,D(:,SPIN))
	      CALL MAC1(NP,NM,U(:,:,SPIN),GP)
	      CALL MACGD(NP,NM,GP,GI(:,:))
	   ELSE
!
              CALL MAC4251(NM,4,U(:,:,SPIN),V(:,L+1))
              NTOP=MOD(L,PARA%NMODUL)
              IF(NTOP.EQ.0) THEN
                 CALL MULD(NP,NM,2,U(:,:,SPIN),D(:,SPIN))
                 CALL GRASCHP(NP,NM,U(:,:,SPIN),D(:,SPIN),R2(:,:,SPIN))
                 CALL MAC1(NP,NM,R2(:,:,SPIN),R(:,:,SPIN))
              ENDIF
!
	      CALL MAC10(NP,NM,AP(:,:,SPIN),GP)
	      CALL MAC1(NP,NM,R(:,:,SPIN),GP)
	      CALL MULD(NP,NM,0,GP,D(:,SPIN))
	      CALL MAC1(NP,NM,U(:,:,SPIN),GP)
	      CALL MACGD(NP,NM,GP,GI(:,:)) 			
   	   END IF
        END IF
!
!
	END SUBROUTINE GICALS
!--------------------------------------------------------------------------------------
! Nom           : GIPCALS        [SUB]
! Fonction      : Calcul des fonctions de Green a temps inegal au second ordre
!                 avec procedure de stabilisation.
!		  Cette sous-routine évalue <c^+(l)c(0)>
! Auteurs       : S.Allen (janvier 1999)
! Revision      : 24 janvier 1999 (Allen)
!--------------------------------------------------------------------------------------
        SUBROUTINE GIPCALS(NM,PARA,G,V,GI,L,SPIN)
        USE CONSTANTES,ONLY: NP,NL
        USE PARA_SIM
        IMPLICIT NONE
        REAL(KIND(1.D0))                                :: CTE
        INTEGER                                         :: J,K,NTOP   ! Var. locales
        INTEGER,INTENT(IN)                              :: NM
        REAL(KIND(1.D0)),DIMENSION (NP,2),SAVE          :: D         	! Var. locale
        REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)   :: G
        REAL(KIND(1.D0)),DIMENSION (NP,NP)              :: GP       ! Var. locales
	REAL(KIND(1.D0)),DIMENSION (NP,NP,2),SAVE	:: AP
        REAL(KIND(1.D0)),DIMENSION (NP,NP,2),SAVE       :: U,R,R2       ! Var. locales
        REAL(KIND(1.D0)),DIMENSION (NP,NL),INTENT(IN)   :: V
        REAL(KIND(1.D0)),DIMENSION(NP,NP),INTENT(OUT)	:: GI
        TYPE(PARAMETRE),INTENT(IN)                      :: PARA
        INTEGER						:: L,SPIN
!
        CTE=-1.D0
        IF(L==0.OR.L==PARA%TRAN) THEN
           IF(L==PARA%TRAN) THEN
              CALL MACGD(NP,NM,G,GI(:,:))
           ELSE
              CALL MAC8(NP,NM,G,AP(:,:,SPIN),CTE)                               ! AP<=-G
              CALL MAC9(NP,NM,AP(:,:,SPIN))                                     ! AP<=1-G
              CALL MACGD(NP,NM,AP(:,:,SPIN),GI(:,:))
           END IF
!
           DO J=1,NM
              DO K=1,NM
                 U(K,J,SPIN)=0.D0
                 R(K,J,SPIN)=0.D0
              END DO
              U(J,J,SPIN)=1.D0
              R(J,J,SPIN)=1.D0
              D(J,SPIN)=1.D0
           END DO
        ELSE
!
           IF(L.LE.PARA%TRAN/2) THEN
              CALL MAC4251(NM,3,U(:,:,SPIN),V(:,L))
              NTOP=MOD(L,PARA%NMODUL)
              IF(NTOP.EQ.0) THEN
                 CALL MULD(NP,NM,0,U(:,:,SPIN),D(:,SPIN))
                 CALL GRASCHP(NP,NM,U(:,:,SPIN),D(:,SPIN),R2(:,:,SPIN))
                 CALL MAC2(NP,NM,R(:,:,SPIN),U(:,:,SPIN))
	         CALL MAC10(NP,NM,R2(:,:,SPIN),U(:,:,SPIN))
              ENDIF
!
              CALL MAC10(NP,NM,AP(:,:,SPIN),GP)
              CALL MAC2(NP,NM,GP,R(:,:,SPIN))
              CALL MULD(NP,NM,2,GP,D(:,SPIN))
              CALL MAC2(NP,NM,GP,U(:,:,SPIN))
              CALL MACGD(NP,NM,GP,GI(:,:))

           ELSE 
        
              CALL MAC4251(NM,0,U(:,:,SPIN),V(:,L+1))
              NTOP=MOD(L,PARA%NMODUL)
              IF(NTOP.EQ.0) THEN
                 CALL MULD(NP,NM,0,U(:,:,SPIN),D(:,SPIN))
                 CALL GRASCHP(NP,NM,U(:,:,SPIN),D(:,SPIN),R2(:,:,SPIN))
                 CALL MAC2(NP,NM,R(:,:,SPIN),U(:,:,SPIN))
	         CALL MAC10(NP,NM,R2(:,:,SPIN),U(:,:,SPIN))
              ENDIF
!
              CALL MAC10(NP,NM,G,GP)
              CALL MAC2(NP,NM,GP,R(:,:,SPIN))
              CALL MULD(NP,NM,2,GP,D(:,SPIN))
              CALL MAC2(NP,NM,GP,U(:,:,SPIN))
              CALL MACGD(NP,NM,GP,GI(:,:))
           END IF
        END IF
!

	CALL MAC11(NP,NM,GI(:,:))
!
        END SUBROUTINE GIPCALS

!--------------------------------------------------------------------------------------
! Nom		: MACGD		[SUB]						MI2
! Fonction	: Calcul des fonctions de Green au 2e ordre a partir du 1er ordre
! Auteurs	: H.Touchette (1996)
! Revision	: 20 novembre 1996 (Touchette)
! Reference	: [1] H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		      hiver 96 - ete 97.
!		  [2] H.G.Matuttis,Ph.D. Thesis,University of Regensburg,1995.
! Notes		: - G(2)=exp(-K/2)*G(1)*exp(K/2) (Multiplication a gauche et a droite)
!		  - Pour laisser les fcts de Green au 1er ordre, mettre les lignes 
!		    MAC1 et MAC2 en commentaires.
!		  
!--------------------------------------------------------------------------------------
	SUBROUTINE MACGD(NP,NM,G,GP)
	USE HOPPING,ONLY: H2P,H2M
	IMPLICIT NONE
	INTEGER,INTENT(IN)				:: NP,NM
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)	:: G
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(OUT) 	:: GP
	CALL MAC10(NP,NM,G,GP)					! GP<=G
!	CALL MAC2(NP,NM,GP,H2P)					! GP<=GP*H2P
!	CALL MAC1(NP,NM,H2M,GP)					! GP<=H2M*GP=H2M*G*H2P
	END SUBROUTINE MACGD
!--------------------------------------------------------------------------------------
! Nom		: FTQ		[SUB]	(Fourier Transform Q-space)		MI3
! Fonction	: Transformee de Fourier des fonctions de Green a temps inegal
! Auteurs	: H.Touchette (1996)
! Revision	: 20 novembre 1996 (Touchette)
! Reference	: [1] H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		      hiver 96 - ete 97.
! Notes		: 
!--------------------------------------------------------------------------------------
	SUBROUTINE FTQ(N,M,G,FG,NV,V)
	USE CONSTANTES,ONLY: NP,ND1,ND2,NL
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))					:: SUM
	INTEGER,INTENT(IN)					:: N,M,NV
	INTEGER,DIMENSION(NP,2)					:: V
	INTEGER							:: K,NM,I,J	! Var locales
	INTEGER							:: KX,KY	! Var locales
	REAL(KIND(1.D0)),DIMENSION(NP,NP),INTENT(IN)		:: G
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2),INTENT(OUT)		:: FG
!
	NM=N*M
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   SUM=0.D0
	   DO I=1,NM
	      DO J=1,NM
	         SUM=SUM+COSMAG(KX,KY,I,J)*G(I,J)
	      END DO
	   END DO
	   FG(KX,KY)=SUM/NM
	END DO
	END SUBROUTINE FTQ
!--------------------------------------------------------------------------------------
! Nom		: SUSCEPI	[SUB]						MI4
! Fonction	: Calcul du facteur de structure de paire a temps inegal et calcul  
!		  de la susceptibilite thermodynamique
! Auteurs	: H.Touchette (1996)
! Revision	: 20 novembre 1996 (Touchette)
! Reference	: [1] H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		      hiver 96 - ete 97.
! Notes		: - Integrale sur tau de X1 en option 
!		  - X2 = SUM(0,TRAN) X(tau) (susceptibilite thermodynamique)
!--------------------------------------------------------------------------------------
	SUBROUTINE SUSCEPI(N,M,GIUP,GIDWN,GIPUP,GIPDN,IQ,X1,NV,V)
	USE CONSTANTES,ONLY: NP,NL,ND1,ND2
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))					:: SUM,CI1,CI2	! loc.
	REAL(KIND(1.D0))					:: GPS1,GPS2	! loc.
	INTEGER							:: K,NM,I,J,KX,KY! loc.
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(OUT)	:: X1
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)		:: GIUP,GIDWN,GIPUP,GIPDN
	INTEGER,INTENT(IN)					:: N,M,IQ,NV
	INTEGER,DIMENSION(NP,2)					:: V
!
	CI1=1.d0
	CI2=1.d0
	IF(IQ.EQ.2)THEN
	   CI2=0.d0
	ELSE IF(IQ.EQ.3)THEN
	   CI1=0.d0
	ENDIF
	NM=N*M
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   SUM=0.D0
	   DO I=1,NM
	      DO J=1,NM
	         GPS1=GIUP(I,J)*GIDWN(I,J)
	         GPS2=GIPUP(J,I)*GIPDN(J,I)
	         SUM=SUM+(GPS1*CI2+GPS2*CI1)*COSMAG(KX,KY,I,J)
	      END DO
	   END DO
	   X1(KX,KY)=SUM/NM
	END DO
	END SUBROUTINE SUSCEPI
!--------------------------------------------------------------------------------------
! Nom		: SUMSUSCEP      [SUB]
! Fonction	: Calcul de l'integrale sur tau de la susceptibilite de paires
! Auteur	: S. Allen
! Revision	: 23 janvier 1999 (Allen)
! Note		: Integration par la methode de Filon
!--------------------------------------------------------------------------------------
	SUBROUTINE SUMSUSCEP(TRAN,DTAU,X1,X2,NV,V)
	USE CONSTANTES,ONLY: NL,ND1,ND2,NP
	USE UTILITAIRE,ONLY: W
	IMPLICIT NONE
	REAL(KIND(1.D0))					:: SUM		! loc
	INTEGER							:: K,KX,KY,L	! loc
	REAL(KIND(1.D0)),INTENT(IN)				:: DTAU
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL),INTENT(IN)	:: X1
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2),INTENT(OUT)		:: X2
	INTEGER,INTENT(IN)					:: TRAN,NV
	INTEGER,DIMENSION(NP,2)					:: V
!
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   SUM=0.D0
	   DO L=0,TRAN
	      SUM=SUM+W(L)*X1(KX,KY,L)
	   END DO
	   X2(KX,KY)=SUM*DTAU
	END DO
	END SUBROUTINE SUMSUSCEP
!--------------------------------------------------------------------------------------
! Nom           : COURANT      [SUB]
! Fonction      : Calcul de la fonction de correlation courant-courant
! Auteur        : S. Allen
! Revision      : 23 janvier 1999 (Allen)
! Note          : Pour evaluer la densité de superfluide (Scalapino et al. PRB 47 p.7995)
!--------------------------------------------------------------------------------------
	SUBROUTINE COURANT(N,M,GIUP,GIDN,GIPUP,GIPDN,GEUP,GEDN,GUP,GDN,X,NV,V)
	USE CONSTANTES,ONLY: NL,NP,ND1,ND2
	USE UTILITAIRE,ONLY: COSMAG
	IMPLICIT NONE
	REAL(KIND(1.D0))					:: SUM,G1,G2
	INTEGER							:: K,I,J,KX,KY,I2,J2,NM
	REAL(KIND(1.D0)),DIMENSION(NP,NP),INTENT(IN)		:: GIUP,GIDN,GIPUP,GIPDN
	REAL(KIND(1.D0)),DIMENSION(NP,NP),INTENT(IN)		:: GEUP,GEDN,GUP,GDN
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2),INTENT(OUT)		:: X
	INTEGER,INTENT(IN)					:: N,M,NV
	INTEGER,DIMENSION(NP,2)					:: V
!
	NM=N*M
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   SUM=0.d0
	   DO I=1,NM
	      I2=I+1
	      IF(MOD(I,N).EQ.0)THEN
		 I2=I2-N
	      ENDIF
	      G1=GEUP(I2,I)-GEUP(I,I2)+GEDN(I2,I)-GEDN(I,I2)
	      DO J=1,NM
	         J2=J+1
		 IF(MOD(J,N).EQ.0)THEN
		    J2=J2-N
		 ENDIF
		 G2=GUP(J2,J)-GUP(J,J2)+GDN(J2,J)-GDN(J,J2)
		 SUM=SUM+(G1*G2-GIUP(I,J2)*GIPUP(I2,J)+GIUP(I2&
		    &,J2)*GIPUP(I,J)+GIUP(I,J)*GIPUP(I2,J2)-&
		    &GIUP(I2,J)*GIPUP(I,J2)-GIDN(I,J2)*GIPDN(&
		    &I2,J)+GIDN(I2,J2)*GIPDN(I,J)+GIDN(I,J)*&
		    &GIPDN(I2,J2)-GIDN(I2,J)*GIPDN(I,J2))*&
		    &COSMAG(KX,KY,J,I)
	      END DO
	      X(KY,KX)=SUM/NM
	   END DO
	END DO
	END SUBROUTINE COURANT
!--------------------------------------------------------------------------------------
! Nom           : SUSMAG       [SUB]
! Fonction      : Calcul de la susceptibilite magnetique
! Auteurs       : S. Allen
!                     25 juin 1999.
!--------------------------------------------------------------------------------------
        SUBROUTINE SUSMAG(N,M,GIUP,GIDWN,GIPUP,GIPDN,GEUP,GEDN,GUP,GDN,X1,NV,V)
        USE CONSTANTES,ONLY: NP,NL,ND1,ND2
        USE UTILITAIRE,ONLY: COSMAG
        IMPLICIT NONE
        REAL(KIND(1.D0))                                        :: SUM,GPS1,GPS2! loc.
        INTEGER                                                 :: K,NM,I,J,KX,KY! loc.
        REAL(KIND(1.D0)),DIMENSION (ND1,ND2),INTENT(OUT)	:: X1
        REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)      	:: GIUP,GIDWN,GIPUP,GIPDN
	REAL(KIND(1.D0)),DIMENSION (NP,NP),INTENT(IN)		:: GEUP,GEDN,GUP,GDN
        INTEGER,INTENT(IN)                                      :: N,M,NV
        INTEGER,DIMENSION(NP,2)					:: V
!
       NM=N*M
       DO K=1,NV
           KX=V(K,1)
           KY=V(K,2)
           SUM=0.D0
           DO I=1,NM
              DO J=1,NM
                 GPS1=GEUP(I,I)-GEDN(I,I)
                 GPS2=GUP(J,J)-GDN(J,J)
                 SUM=SUM+(GPS1*GPS2+GIUP(I,J)*GIPUP(I,J)+GIDWN(I,J)&!
	             &*GIPDN(I,J))*COSMAG(KX,KY,I,J)
              END DO
           END DO
           X1(KX,KY)=SUM/NM
        END DO

        END SUBROUTINE SUSMAG
!--------------------------------------------------------------------------------------
! Nom		: MATSUBARA	[SUB]						MI5
! Fonction	: Calcul des fonctions de Green en frequences de Matsubara
! Auteurs	: H.Touchette (1996)
! Revision	: 20 novembre 1996 (Touchette)
! Reference	: [1] H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		      hiver 96 - ete 97.
! Notes		: - G(iwn,Qx,Qy) = Re{G(iwn,Qx,Qy)} + Im{G(iwn,Qx,Qy)}
!		    Re{G(iwn,Qx,Qy)} - IST=0
!		    Im{G(iwn,Qx,Qy)} - IST=1
!		  - Integration selon la regle de Filon (cf ref)
!--------------------------------------------------------------------------------------
	SUBROUTINE MATSUBARA(NTRAN,BETA,IST,GI,GIM,NV,V)
	USE CONSTANTES,ONLY: NL,ND1,ND2,NP
	IMPLICIT NONE
	REAL(KIND(1.D0))					:: PI,WMAT,DTAU	! Var locales
	REAL(KIND(1.D0)),INTENT(IN)				:: BETA
	INTEGER,INTENT(IN)					:: NTRAN,IST,NV
	INTEGER,DIMENSION(NP,2)					:: V
	INTEGER							:: K,KX,KY,LW	! Var locales
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL),INTENT(IN)	:: GI
	REAL(KIND(1.D0)),DIMENSION (ND1,ND2,0:NL),INTENT(OUT)	:: GIM
!
	PI=DACOS(-1.D0)
	DTAU=BETA/DBLE(NTRAN)
	DO K=1,NV
	   KX=V(K,1)
	   KY=V(K,2)
	   DO LW=0,NTRAN
	      WMAT=PI*DBLE(2*LW+1)/BETA
	      SELECT CASE(IST)
	         CASE(0)						
	            CALL FILON(0.D0,BETA,WMAT,NTRAN,0,GI(KX,KY,:),GIM(LW,KX,KY))
	         CASE(1)						 
	            CALL FILON(0.D0,BETA,WMAT,NTRAN,1,GI(KX,KY,:),GIM(LW,KX,KY))
	         CASE DEFAULT
	            WRITE(*,*) '--- Erreur : MATSUBARA ---'
	       END SELECT 
	   END DO
	END DO
	END SUBROUTINE MATSUBARA
!--------------------------------------------------------------------------------------
! Nom		: FILON		[SUB]						MI6
! Fonction	: Integration selon la regle de Filon
! Auteurs	: Liang Chen (1996)
! Revision	: 4 decembre 1996 (Touchette)
! Reference	: [1] H.Touchette, Rapports des projets Monte Carlo quantique, CRPS,
!		      hiver 96 - ete 97.
!		  [2] Davis & Rabinowitz, Methods of Numerical Integration, Academic
!		      Press, 1984. QA 299.3 D28 1984
!		  [3] Davis & Rabinowitz, Numerical Integration, Blaisdell, 1967.
!		      QA 299.3 D3N 1967
! Notes		: - IST=0 : Integrale de f(x)*cos(t*x)
!		  - IST=1 : Integrale de f(x)*sin(t*x)
!		  - Utilise un nombre impair de points: de 0 a N_EVEN 
!
!--------------------------------------------------------------------------------------
	SUBROUTINE FILON(A,B,K,N_EVEN,IST,Y,RESULT)
	USE CONSTANTES,ONLY: NL
	IMPLICIT NONE
	REAL(KIND(1.D0))				:: H,ALPHA,BETA,GAMMA	! Var locales
	REAL(KIND(1.D0))				:: C2N,C2N_1,KH,KA,KB	! Var locales
	REAL(KIND(1.D0))				:: SKH,CKH,KH2,KH3	! Var locales
	REAL(KIND(1.D0)),INTENT(IN)			:: A,B,K
	REAL(KIND(1.D0)),INTENT(OUT)			:: RESULT
	INTEGER,INTENT(IN)				:: N_EVEN,IST
	INTEGER						:: I
	REAL(KIND(1.D0)),DIMENSION (0:NL),INTENT(IN)	:: Y
!
	H=(B-A)/DBLE(N_EVEN)
	KH=K*H
	KH2=KH*KH
	KH3=KH2*KH
	KA=K*A
	KB=K*B
	SKH=DSIN(KH)
	CKH=DCOS(KH)
	ALPHA=(KH2+KH*SKH*CKH-2.D0*SKH*SKH)/KH3
	BETA=2.D0*(KH*(1.D0+CKH*CKH)-2.D0*SKH*CKH)/KH3
	GAMMA=4.D0*(SKH-KH*CKH)/KH3
	SELECT CASE(IST)
	   CASE(0)
	      C2N=0.5D0*Y(0)*DCOS(KA)+0.5D0*Y(N_EVEN)*DCOS(KB)
	      DO I=2,N_EVEN-1,2
	         C2N=C2N+Y(I)*DCOS(KA+DBLE(I)*KH)
	      END DO
	      C2N_1=0.D0
	      DO I=1,N_EVEN,2
	         C2N_1=C2N_1+Y(I)*DCOS(KA+DBLE(I)*KH)
	      END DO
	      RESULT=H*(ALPHA*(Y(N_EVEN)*DSIN(KB)-Y(0)*DSIN(KA))+BETA*C2N+GAMMA*C2N_1)
	   CASE(1)
	      C2N=0.5D0*Y(0)*DSIN(KA)+0.5D0*Y(N_EVEN)*DSIN(KB)
	      DO I=2,N_EVEN-1,2
	         C2N=C2N+Y(I)*DSIN(KA+DBLE(I)*KH)
	      END DO
	      C2N_1=0.D0
	      DO I=1,N_EVEN,2
	         C2N_1=C2N_1+Y(I)*DSIN(KA+DBLE(I)*KH)
	      END DO
	      RESULT=H*(ALPHA*(-Y(N_EVEN)*DCOS(KB)+Y(0)*DCOS(KA))+BETA*C2N+GAMMA*C2N_1) 
	   CASE DEFAULT
	      WRITE(*,*) '--- Erreur : FILON ---'
	END SELECT
	END SUBROUTINE FILON
!--------------------------------------------------------------------------------------
! Nom		: SELF		[SUB]						MI7
! Fonction	: Calcul de la self-energie en frequences de Matsubara
! Auteurs	: H.Touchette
! Revision	: 20 juin 1997 (Touchette)
! References	: H.Touchette, Aspects numeriques des simulations Monte Carlo quantique 
!		  du modele de Hubbard, rapport technique du CRPS, departement de 
!		  physique, UdeS, ete 1997.
! Notes		: - SR = Re{sigma(iwn,Qx,Qy)}
!		  - SI = Im{sigma(iwn,Qx,Qy)}
!--------------------------------------------------------------------------------------
	SUBROUTINE SELF(PARA,GR,GI,SR,SI)
	USE CONSTANTES,ONLY: NL,ND1,ND2
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER							:: K,I,J,L
	REAL(KIND(1.D0))					:: WMAT,PI,DEN,ENE
	REAL(KIND(1.D0))					:: QX,QY,ZETAK
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL),INTENT(IN)	:: GR,GI
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL),INTENT(OUT) 	:: SR,SI
	TYPE(PARAMETRE),INTENT(IN)				:: PARA
!
	PI=DACOS(-1.D0)
	DO K=1,PARA%NVEC
	   I=PARA%VEC(K,1)
	   J=PARA%VEC(K,2)
	   QX=2.D0*PI*DBLE(I)/DBLE(PARA%N)
	   QY=2.D0*PI*DBLE(J)/DBLE(PARA%M)
	   ZETAK=ENE(PARA%N,PARA%M,PARA%TPAR,PARA%TPER,PARA%TPRM,QX,QY)-PARA%MU
	   DO L=0,PARA%TRAN
	      WMAT=PI*DBLE(2*L+1)/PARA%BETA
	      DEN=GR(I,J,L)*GR(I,J,L)+GI(I,J,L)*GI(I,J,L)
	      SR(I,J,L)=-ZETAK-GR(I,J,L)/DEN
	      SI(I,J,L)=WMAT+GI(I,J,L)/DEN
	   END DO
	END DO
	END SUBROUTINE SELF
!--------------------------------------------------------------------------------------
! Nom		: COVARIANCE		[SUB]					MI8
! Fonction	: Calcul de la covariance pour la methode d'entropie maximum
! Auteurs	: H.Touchette
! References	: [1] H.Touchette, Aspects numeriques des simulations Monte Carlo  
!		      quantique du modele de Hubbard, rapport technique du CRPS,  
!		      departement de physique, UdeS, ete 1997.
!		  [2] Jarrel et al., Phys.Rep.,269,1996,p.145.
!		  [3] Gubernatis et al., PRB,4,12,1991,p.6011.
! Notes		: - Matrice de covariance COVM 
!		  - On calcule cette matrice seulement pour une composante de spin
!		    et a la fin d'une simulation
!--------------------------------------------------------------------------------------
	SUBROUTINE COVARIANCE(PARA,C1,COVG,XG,COVM)
	USE CONSTANTES,ONLY: NL,ND1,ND2,NCOR
	USE PARA_SIM
	IMPLICIT NONE
	INTEGER								:: L1,L2	! Loc.
	INTEGER								:: I,J,K,NMES	! Loc.
	REAL(KIND(1.D0))						:: SUM		! Loc.
	REAL(KIND(1.D0))						:: D1,D2	! Loc.
	REAL(KIND(1.D0)),INTENT(IN)					:: C1
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL),INTENT(IN)		:: XG
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL,NCOR),INTENT(IN)	:: COVG
	REAL(KIND(1.D0)),DIMENSION(ND1,ND2,0:NL,0:NL),INTENT(OUT)	:: COVM
	TYPE(PARAMETRE),INTENT(IN)					:: PARA
!
	NMES=INT(C1)
	DO I=1,PARA%N
	   DO J=1,PARA%M
	      DO L1=0,PARA%TRAN
	         DO L2=0,PARA%TRAN
	            SUM=0.D0
	            DO K=1,NMES
	               D1=(XG(I,J,L1)-COVG(I,J,L1,K))
	               D2=(XG(I,J,L2)-COVG(I,J,L2,K))
	               SUM=SUM+D1*D2
	            END DO
	            COVM(I,J,L2,L1)=SUM/(C1*(C1-1.D0))
	         END DO
	      END DO
	   END DO
	END DO
	END SUBROUTINE COVARIANCE
!--------------------------------------------------------------------------------------
!       ----------       FONCTION NUMTOASCII        ----------
!       Fonction servant a creer des suffixes numeraux a
!       partir de la liste ascii de 0 a 9.
!       Cf. Christian Foisy (1996), (S. Allen version f90 1998)

        CHARACTER(*) FUNCTION NUMTOASCII(I)
        IMPLICIT NONE
        INTEGER                                         :: J,K
        INTEGER,INTENT(IN)                              :: I
        J=I
        K=LEN(NUMTOASCII)
	NUMTOASCII(1:K-1)=' '
	NUMTOASCII(K:K)='0'
        DO WHILE (J.GT.0.AND.K.GT.0)
           NUMTOASCII(K:K)=CHAR(ICHAR('0')+MOD(J,10))
           K=K-1
           J=J/10
        END DO
        NUMTOASCII=ADJUSTL(NUMTOASCII)
        END FUNCTION NUMTOASCII
!
        INTEGER FUNCTION ASCIITONUM(CARACT)
        IMPLICIT NONE
        CHARACTER(LEN=*),INTENT(IN)                     :: CARACT
        INTEGER                                         :: I0,J,L,K,NUM
	I0=IACHAR('0')
        K=0
        NUM=0
        DO
            K=K+1
            IF(CARACT(K:K).NE.' ') EXIT
        END DO
        J=INDEX(CARACT(K:),' ')
        IF(J.EQ.0)THEN
            J=LEN(CARACT)+1
        ENDIF
        DO L=K,J-1
            NUM=NUM+10**(J-L-1)*(IACHAR(CARACT(L:L))-I0)
        END DO
	ASCIITONUM=NUM
        END FUNCTION ASCIITONUM
