CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        PROGRAM ISING
C     Monte Carlo simulation of the two-dimensional Ising Model
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL INIT         !get input from screen
      CALL ARCHON       !calculate the thermodynamic quantities
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ARCHON  
C calculates the thermodynamic quantities:  energy, magnetization
C susceptibility, and specific heat at constant field and intrctn strngth
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Local variables:
      INTEGER SPIN(MAXX,MAXY)    !spin configuration
C     all of the thermodynamic quant have 2 indices 
C     (not all array elements are used, e.g. CHI(sweep,value))
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL MAG(3,3)              !magnetization
      REAL ENERGY(3,3)           !energy
      REAL CB(3,3)               !specific heat
      REAL CHI(3,3)              !susceptibility
      INTEGER ITHERM             !thermalization index
      INTEGER ITER               !iteration index
      REAL ACCPT                 !acceptance ratio
      INTEGER IX,IY              !horiz and vert indices
      INTEGER MORE,IGRP          !how many more groups, group index
      INTEGER ISWEEP             !sweep index
      INTEGER SWEEP,GROUP,TOTAL  !which level of calculation
C Functions:
      REAL RANNOS                !generates a random number
      DATA SWEEP,GROUP,TOTAL/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      NUNIT=12
      OPEN(NUNIT,FILE='EX8.OUT',STATUS='UNKNOWN')
      CALL PRMOUT(NUNIT)         !write out parameter summary
C
      DO 5 IX=1,NX              !random initial spin configuration
         DO 6 IY=1,NY
              IF (RANNOS(DSEED) .GT. .5) THEN
                 SPIN(IX,IY)=1
              ELSE
                 SPIN (IX,IY)=-1
              END IF
6        CONTINUE
5     CONTINUE
C
      DO 10 ITHERM=1,NTHERM     !thermalize init config
         CALL METROP(SPIN,ACCPT)
         WRITE (NUNIT,7) ITHERM,ACCPT
7        FORMAT (5X,' thermalization sweep ',I4,
     +           ', acceptance ratio =',F5.3) 
10    CONTINUE            
C
      CALL ZERO(TOTAL,MAG,ENERGY,CHI,CB)  !zero total averages
      MORE=NGROUP
C
15    CONTINUE                                !allow for more groups
        DO 20 IGRP=NGROUP-MORE+1,NGROUP       !loop over groups
           CALL ZERO(GROUP,MAG,ENERGY,CHI,CB) !zero group averages
C
           CALL TITLES(NUNIT)
C
           DO 30 ITER=1,NFREQ*NSIZE
              CALL METROP(SPIN,ACCPT)     !make a sweep of the lattice
C                               
              IF (MOD(ITER,NFREQ) .EQ. 0) THEN  !include in averages 
                  ISWEEP=ITER/NFREQ               !which sweep is it
                  CALL ZERO(SWEEP,MAG,ENERGY,CHI,CB) !zero sweep totals
                  CALL SUM(SPIN,MAG,ENERGY)!sweep totals, add to group
                  CALL  
     +              SWPOUT(NUNIT,MAG,ENERGY,ACCPT,IGRP,ISWEEP)
              END IF
30         CONTINUE  
           CALL AVERAG(MAG,ENERGY,CHI,CB,IGRP,NUNIT) !calc total averages
20      CONTINUE
C
      PRINT*, ' How many more groups? [0 - stop]'
      READ(*,*)MORE
C
      IF(MORE.EQ.0) THEN
        WRITE(*,*)' *** See output in EX8.OUT ***'
        WRITE(NUNIT,*)' '
        WRITE(NUNIT,*)
     &    ' *********   Spin lattice  (-1=down, 1=up)  *********'
        WRITE(NUNIT,*)' '
        DO IX=1,NX
          WRITE(NUNIT,123)(SPIN(IX,IY),IY=1,NY)
123       FORMAT(70(I2))
        ENDDO
      ENDIF
C
      IF (MORE .GT. 0) THEN
           NGROUP=NGROUP+MORE
           GOTO 15
         ELSE
           STOP
      END IF
C     
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE METROP(SPIN,ACCPT)
C make one sweep of the lattic using the Metropolis algorithm 
C to generate a new configuration
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Input/Output variables:
      INTEGER SPIN(MAXX,MAXY)     !spin configuration (I/O)
      REAL ACCPT                  !acceptance ratio (output)
C Local variables:
      INTEGER IX,IY               !horiz and vert indices
      INTEGER IXM1,IXP1,IYM1,IYP1 !indices of nearest neighbors
      INTEGER NNSUM               !sum of nearest neighbors
C Function:
      REAL RANNOS                 !generates a random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      ACCPT=0.                          !zero acceptance ratio
      DO 10 IX=1,NX              
        IXP1=IX+1                       !nearest neighbors
        IF (IX .EQ. NX) IXP1=1          !with periodic b.c.
        IXM1=IX-1
        IF (IX .EQ. 1) IXM1=NX
C     
        DO 20 IY=1,NY
          IYP1=IY+1                     !nearest neighbors
          IF (IY .EQ. NY) IYP1=1        !with periodic b.c.
          IYM1=IY-1
          IF (IY .EQ. 1) IYM1=NY
C                                       !term to weight new config
          NNSUM=SPIN(IX,IYP1)+SPIN(IX,IYM1)+SPIN(IXP1,IY)+SPIN(IXM1,IY)
C
          IF (RANNOS(DSEED) .LT. RATIO(NNSUM,SPIN(IX,IY))) THEN
             SPIN(IX,IY)=-SPIN(IX,IY)   !flip the spin
             ACCPT=ACCPT+1              !update accept count
          END IF
C
20      CONTINUE
10    CONTINUE
      ACCPT=ACCPT/NSPIN                 !make it a ratio
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SUM(SPIN,MAG,ENERGY)
C calculate magnetization and energy for this sweep
C add these values to the group averages
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Input/output variables:
      INTEGER SPIN(MAXX,MAXY)    !spin configuration (input)
C     all of the thermodynamic quant have 2 indices
C     (not all array elements are used, e.g. CHI(sweep,value))
C     first index is the level: sweep, group, or total
C     second index is the quantity: value, square, or sigma**2
      REAL MAG(3,3)              !magnetization (I/O)
      REAL ENERGY(3,3)           !energy (I/O)
C Local variables:
      INTEGER PAIRS              !interaction sum
      INTEGER SWEEP,GROUP        !which level of calculation
      INTEGER VALUE,SQUARE       !which quantity
      INTEGER IX,IY              !horiz and vert indices
      INTEGER IXM1,IYM1          !neighbor indices
      DATA SWEEP,GROUP /1,2/
      DATA VALUE,SQUARE/1,2/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PAIRS=0                      !zero pair sum
      DO 10 IY=1,NY
         IYM1=IY-1                 !neighbor just below
         IF (IY .EQ. 1) IYM1=NY    !periodic b.c.
         DO 20 IX=1,NX
           IXM1=IX-1               !neighbor to the left
           IF (IX .EQ. 1) IXM1=NX  !periodic b.c.
C
C          this method of summing pairs does not count twice
           PAIRS=PAIRS+SPIN(IX,IY)*(SPIN(IX,IYM1)+SPIN(IXM1,IY))
C          magnetization is the sum of the spins  (Eq. 8.21a)
           MAG(SWEEP,VALUE)=MAG(SWEEP,VALUE)+SPIN(IX,IY)
C
20       CONTINUE
10    CONTINUE
C
      MAG(SWEEP,SQUARE)=MAG(SWEEP,VALUE)**2
      ENERGY(SWEEP,VALUE)=-J*PAIRS-B*MAG(SWEEP,VALUE)   !Eq 8.18
      ENERGY(SWEEP,SQUARE)=ENERGY(SWEEP,VALUE)**2
C
C     add sweep contributions to group sums
      MAG(GROUP,VALUE)=MAG(GROUP,VALUE)+MAG(SWEEP,VALUE)
      MAG(GROUP,SQUARE)=MAG(GROUP,SQUARE)+MAG(SWEEP,SQUARE)
      ENERGY(GROUP,VALUE)=ENERGY(GROUP,VALUE)+ENERGY(SWEEP,VALUE)
      ENERGY(GROUP,SQUARE)=ENERGY(GROUP,SQUARE)+ENERGY(SWEEP,SQUARE)
C
      RETURN
      END      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE AVERAG(MAG,ENERGY,CHI,CB,IGROUP,NUNIT) 
C find group averages from group sums and add these to total averages;
C calculate uncertainties and display results
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Input/Output variables:
C     all of the thermodynamic quant have 2 indicies
C     (not all array elements are used, e.g. CHI(sweep,value))
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL MAG(3,3)              !magnetization
      REAL ENERGY(3,3)           !energy
      REAL CB(3,3)               !specific heat
      REAL CHI(3,3)              !susceptibility
      INTEGER IGROUP             !group index (input)
C Local variables:
      REAL M,MSIG1,MSIG2         !magnetization and uncertainties
      REAL E,ESIG1,ESIG2         !energy and uncertainties
      REAL SUS,SUSSIG            !susceptibility and uncertainty
      REAL C,CSIG                !specific heat and uncertainty
      INTEGER IQUANT             !index the quantity
      INTEGER GROUP,TOTAL        !which level of calculation
      INTEGER VALUE,SQUARE,SIGSQ !which quantity
      DATA GROUP,TOTAL/2,3/
      DATA VALUE,SQUARE,SIGSQ/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     calculate group averages and uncertainties from group sums
      DO 10 IQUANT=VALUE,SQUARE
         MAG(GROUP,IQUANT)=MAG(GROUP,IQUANT)/NSIZE
         ENERGY(GROUP,IQUANT)=ENERGY(GROUP,IQUANT)/NSIZE
10    CONTINUE
      CHI(GROUP,VALUE)=MAG(GROUP,SQUARE)-MAG(GROUP,VALUE)**2
      MAG(GROUP,SIGSQ)=CHI(GROUP,VALUE)/NSIZE
      IF (MAG(GROUP,SIGSQ) .LT. 0.) MAG(GROUP,SIGSQ)=0.
      CB(GROUP,VALUE)=ENERGY(GROUP,SQUARE)-ENERGY(GROUP,VALUE)**2
      ENERGY(GROUP,SIGSQ)=CB(GROUP,VALUE)/NSIZE
      IF (ENERGY(GROUP,SIGSQ) .LT. 0.) ENERGY(GROUP,SIGSQ)=0.
      CHI(GROUP,SQUARE)=CHI(GROUP,VALUE)**2
      CB(GROUP,SQUARE)=CB(GROUP,VALUE)**2
C
C     add group averages to total sums
      DO 20 IQUANT=VALUE,SIGSQ
         MAG(TOTAL,IQUANT)=MAG(TOTAL,IQUANT)+MAG(GROUP,IQUANT)
         ENERGY(TOTAL,IQUANT)=ENERGY(TOTAL,IQUANT)+ENERGY(GROUP,IQUANT)
         CHI(TOTAL,IQUANT)=CHI(TOTAL,IQUANT)+CHI(GROUP,IQUANT)
         CB(TOTAL,IQUANT)=CB(TOTAL,IQUANT)+CB(GROUP,IQUANT)
20    CONTINUE
C
C     find total averages using total sums accumulated so far
      M=MAG(TOTAL,VALUE)/IGROUP
      MSIG1=(MAG(TOTAL,SQUARE)/IGROUP-M**2)/IGROUP/NSIZE
      IF (MSIG1 .LT. 0) MSIG1=0.
      MSIG1=SQRT(MSIG1)
      MSIG2=SQRT(MAG(TOTAL,SIGSQ))/IGROUP
C
      E=ENERGY(TOTAL,VALUE)/IGROUP
      ESIG1=(ENERGY(TOTAL,SQUARE)/IGROUP-E**2)/IGROUP/NSIZE
      IF (ESIG1 .LT. 0) ESIG1=0.
      ESIG1=SQRT(ESIG1)
      ESIG2=SQRT(ENERGY(TOTAL,SIGSQ))/IGROUP
C
      SUS=CHI(TOTAL,VALUE)/IGROUP
      SUSSIG=(CHI(TOTAL,SQUARE)/IGROUP-SUS**2)/IGROUP
      IF (SUSSIG .LT. 0.) SUSSIG=0.
      SUSSIG=SQRT(SUSSIG)
C
      C=CB(TOTAL,VALUE)/IGROUP
      CSIG=(CB(TOTAL,SQUARE)/IGROUP-C**2)/IGROUP
      IF (CSIG .LT. 0.) CSIG=0.
      CSIG=SQRT(CSIG)
C
C     write out summary
      CALL TXTOUT(MAG,ENERGY,CB,CHI,E,ESIG1,ESIG2,
     +    M,MSIG1,MSIG2,SUS,SUSSIG,C,CSIG,IGROUP,NUNIT)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ZERO(ILEVEL,MAG,ENERGY,CHI,CB)
C zero sums for ILEVEL thermodynamic values
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Input/Output variables:
      INTEGER ILEVEL             !which level to zero  (input)
C     all of the thermodynamic quant have 2 indices
C     (not all array elements are used, e.g. CHI(sweep,value))
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL MAG(3,3)              !magnetization (output)
      REAL ENERGY(3,3)           !energy (output)
      REAL CB(3,3)               !specific heat (output)
      REAL CHI(3,3)              !susceptibility (output)
C Local variable:                                                              
      INTEGER IQUANT             !which quantity
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 10 IQUANT=1,3 
           MAG(ILEVEL,IQUANT)=0.
           ENERGY(ILEVEL,IQUANT)=0.
           CHI(ILEVEL,IQUANT)=0.
           CB(ILEVEL,IQUANT)=0.
10    CONTINUE
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INIT
C intializes constants, displays header screen,
C gets parameters from screen
C intializes menu arrays for input parameters
C calculates all derivative parameters
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Local variables:
      INTEGER IF              !possible values for sum of neighb. spins  
      INTEGER IOPT            !option for screen input
      INTEGER IDSEED          !input initial seed from screen
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C     display header screen     
      PRINT*, 'EXAMPLE 8'
      PRINT*, 'Monte Carlo simulation of the 2-D Ising Model'
      PRINT*, 'using the Metropolis algorithm'
C 
      B=0.
      J=3.
      NX=20
      NY=20
      DSEED=DBLE(54767)
      NTHERM=20
      NFREQ=5
      NSIZE=10
      NGROUP=10
      PRINT*, ' Want "default" run? [0=yes, 1=no]'
      READ(*,*)IOPT
      IF(IOPT.EQ.0) GO TO 100
C
      PRINT*,'Enter value for magnetic field (units of kT)'
      PRINT*,' [-20 to 20]'  
      READ(*,*)B
C
      PRINT*,'Enter value for interaction strength (units of kT)'
      PRINT*,' [-20 to 20]'  
      READ(*,*)J
C 
      PRINT*,'Enter number of X lattice points [2-79]'
      READ(*,*)NX
C 
      PRINT*, 'Enter number of Y lattice points [2-20]'
      READ(*,*)NY
C 
      PRINT*, 'Integer random number seed for init fluctuations'
      PRINT*,' [1000 to 99999]'
      READ(*,*)IDSEED
      DSEED=DBLE(IDSEED)
C 
      PRINT*, 'Number of thermalization sweeps [0 - 1000]'
      READ(*,*)NTHERM
C 
      PRINT*,'Enter sampling frequency (to avoid correlations)'
      PRINT*,' [1 - 100]'
      READ(*,*)NFREQ
C 
      PRINT*, 'Number of samples in a group [1 - 1000]'
      READ(*,*)NSIZE
C                   
      PRINT*, 'Enter number of groups [1 - 1000]'
      READ(*,*)NGROUP
C
C     calculate derivative parameters
100   NSPIN=NX*NY
      DO 10 IF=-4,4,2    !ratio of prob.; not all matrix elem are used
            RATIO(IF,-1)=EXP(2*(J*IF+B))
            RATIO(IF,1) =1./RATIO(IF,-1)
10    CONTINUE
C 
      RETURN
      END                            
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PRMOUT(NUNIT)
C write out parameter summary to NUNIT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
      INTEGER NUNIT                   !output unit number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      WRITE (NUNIT,5)
      WRITE (NUNIT,6)
      WRITE (NUNIT,7) B
      WRITE (NUNIT,8) J
      WRITE (NUNIT,10) NX,NY
      WRITE (NUNIT,15) NTHERM
      WRITE (NUNIT,20) NFREQ,NSIZE,NGROUP
      WRITE (NUNIT,*) '**********************************************'
C      
5     FORMAT (' Output from example 8:')
6     FORMAT (' Monte Carlo Simulation of the',/,
     +        ' 2-D Ising Model using the Metropolis Algorithm',/,
     +        ' **********************************************')
7     FORMAT (' Magnetic field (units of kT) =',1PE12.5)
8     FORMAT (' Interaction strength (units of kT) =',1PE12.5)
10    FORMAT (' NX =',I3,5X,' NY =',I3)
15    FORMAT (' number of thermalization sweeps =',I4)
20    FORMAT (' sweep frequency = ',I4,' group size =',I4,/,
     +        ' number of groups =', I4)
C                               
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SWPOUT(NUNIT,MAG,ENERGY,ACCPT,IGROUP,ISWEEP)
C and write out data for this sweep
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Input variables:
C     all of the thermodynamic quant have 2 indicies
C     (not all array elements are used, e.g. CHI(sweep,value))
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL MAG(3,3)              !magnetization
      REAL ENERGY(3,3)           !energy
      REAL ACCPT                 !acceptance ratio
      INTEGER ISWEEP,IGROUP      !sweep and group index
      INTEGER NUNIT              !output unit number
C Local variables:
      INTEGER SWEEP              !which level of calculation
      INTEGER VALUE              !which quantity
      DATA SWEEP /1/
      DATA VALUE /1/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      WRITE (NUNIT,11) IGROUP,ISWEEP,NSIZE,ACCPT,
     +       ENERGY(SWEEP,VALUE)/NSPIN,MAG(SWEEP,VALUE)/NSPIN
11    FORMAT (7X,I3,7X,I3,'/',I3,7X,F5.3,7X,F9.5,7X,F9.3)
      RETURN
      END  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TXTOUT(MAG,ENERGY,CB,CHI,E,ESIG1,ESIG2,
     +    M,MSIG1,MSIG2,SUS,SUSSIG,C,CSIG,IGROUP,NUNIT)
C write out averages and uncertaintes to NUNIT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E8'
C Input variables:
C     all of the thermodynamic quant have 2 indicies
C     (not all array elements are used, e.g. CHI(sweep,value))
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL MAG(3,3)              !magnetization
      REAL ENERGY(3,3)           !energy
      REAL CB(3,3)               !specific heat
      REAL CHI(3,3)              !susceptibility
      INTEGER IGROUP             !group index
      REAL M,MSIG1,MSIG2         !magnetization and uncertainties
      REAL E,ESIG1,ESIG2         !energy and uncertainties
      REAL SUS,SUSSIG            !susceptibility and uncertainty
      REAL C,CSIG                !specific heat and uncertainty
      INTEGER NUNIT              !output unit number
C Local variables:
      INTEGER GROUP              !which level of calculation
      INTEGER VALUE,SIGSQ        !which quantity
      DATA GROUP /2/
      DATA VALUE,SIGSQ/1,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      WRITE (NUNIT,30) IGROUP,NGROUP
      WRITE (NUNIT,32)
      WRITE (NUNIT,33)
      WRITE (NUNIT,35) 
     +  ENERGY(GROUP,VALUE)/NSPIN,SQRT(ENERGY(GROUP,SIGSQ))/NSPIN,
     +  MAG(GROUP,VALUE)/NSPIN,SQRT(MAG(GROUP,SIGSQ))/NSPIN,
     +  CHI(GROUP,VALUE)/NSPIN,CB(GROUP,VALUE)/NSPIN
      WRITE (NUNIT,40) E/NSPIN,ESIG1/NSPIN,ESIG2/NSPIN,
     +  M/NSPIN,MSIG1/NSPIN,MSIG2/NSPIN,
     +  SUS/NSPIN,SUSSIG/NSPIN,C/NSPIN,CSIG/NSPIN
C
      WRITE (*,*) ' '
      WRITE (*,30) IGROUP,NGROUP
      WRITE (*,32)
      WRITE (*,33)
      WRITE (*,35) 
     +  ENERGY(GROUP,VALUE)/NSPIN,SQRT(ENERGY(GROUP,SIGSQ))/NSPIN,
     +  MAG(GROUP,VALUE)/NSPIN,SQRT(MAG(GROUP,SIGSQ))/NSPIN,
     +  CHI(GROUP,VALUE)/NSPIN,CB(GROUP,VALUE)/NSPIN
      WRITE (*,40) E/NSPIN,ESIG1/NSPIN,ESIG2/NSPIN,
     +  M/NSPIN,MSIG1/NSPIN,MSIG2/NSPIN,
     +  SUS/NSPIN,SUSSIG/NSPIN,C/NSPIN,CSIG/NSPIN
      WRITE (*,*) ' '
C
30    FORMAT ('  Group ',I3,' (out of ',I4,') averages')
32    FORMAT (14X,'Energy',13X,'Magnetization',5X,'Susceptibility',
     +        2X,'Specific Heat')
33    FORMAT (14X,'------',13X,'-------------',5X,'--------------',
     +        2X,'-------------')
35    FORMAT ('  group ',2(1X,F7.3,'+-',F5.3,6X),
     +        2(2X,F6.3,7X))
40    FORMAT ('  total ',2(1X,F7.3,'+-',F5.3,'/',F5.3),
     +        2(1X,F6.3,'+-',F6.3))
C
      RETURN
      END                                   
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TITLES(NUNIT)
C write out text data titles
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Input/Output variables:
      INTEGER NUNIT                   !output unit number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      WRITE (NUNIT,10)
      WRITE (NUNIT,11)
10    FORMAT (6X,'Group',3X,'Sweep/Out of',5X,'Accpt',9X,'Energy',
     +       6X,'Magnetization')      
11    FORMAT (6X,'-----',3X,'------------',5X,'-----',9X,'------',
     +       6X,'-------------')      
C
      RETURN                                                            
      END
                                                                    
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL FUNCTION RANNOS(DSEED)
C returns a uniformly distributed random number between 0 and 1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DOUBLE PRECISION         DSEED
      DOUBLE PRECISION         D2P31M,D2P31
      DATA             D2P31M/2147483647.D0/
      DATA             D2P31 /2147483711.D0/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DSEED = MOD(16807.D0*DSEED,D2P31M)
      RANNOS = DSEED / D2P31
      RETURN
      END