   
C
C  Take measurements of energy and magnetization.
C

C--------------------------------------------------------
      SUBROUTINE init_measure
      INCLUDE "exp_wolff.h"

C--- Initialize sums and sum of squares.
      sum_energy = 0.0
      sum_spins  = 0.0
      sum_enersq = 0.0
      sum_spinsq = 0.0
      sum_size   = 0.0
      sum_sizesq = 0.0

      RETURN
      END


C---------------------------------------------------------
       SUBROUTINE  measure
       INCLUDE "exp_wolff.h"
       integer nt,Nr
       INTEGER sum_sp, sum_ene
       REAL    ene, mag
       INTEGER i,j,k
       INTEGER d1,d2,d3,type,bytes
       LOGICAL  rbound
       EXTERNAL rbound
       INTEGER  length
       INTEGER inbuf(0:intsize*Nrow_max*Ncol_max)
       INTEGER outbuf(0:intsize*Nrow_max*Ncol_max)


       length = intsize*Nrow_max*Ncol_max
       
       k = 0
       IF( procnum .GT. 0 ) THEN  
        DO j=0, ntile-1
             DO i=0, Nrow-1 
                outbuf(k) = spin(i,j*width)
                k=k+1
             ENDDO
          ENDDO
       ELSE
          DO j=1, ntile-1
             DO i=0, Nrow-1
                outbuf(k) = spin(i,j*width)
                k=k+1
             ENDDO
          ENDDO
          DO i=0, Nrow-1
             outbuf(k) = spin(i,0)
             k=k+1
          ENDDO
       ENDIF
    
       type=70
       d1 = KXCHAN(inbuf, length, rpnum, type,
     *             outbuf,length, lpnum, type)

       sum_ene = 0

       DO j=0, Ncol-1
          IF( .NOT. rbound(j) )  THEN
             DO i=0, Nrow-2
                sum_ene = sum_ene + (-spin(i,j)) 
     &                   * ( spin(i,j+1) + spin(i+1,j) )
             ENDDO
             sum_ene = sum_ene + (-spin(Nrow-1,j))
     &             * ( spin(Nrow-1,j+1) + spin(0,j) )
          ELSE
             tnum = j / width
             DO i=0, Nrow-2
                sum_ene = sum_ene + (-spin(i,j))
     &                   * ( inbuf(i+tnum*Nrow) + spin(i+1,j) )
             ENDDO
             sum_ene = sum_ene + (-spin(Nrow-1,j))
     &             * ( inbuf(Nrow-1+tnum*Nrow) + spin(0,j) )
          ENDIF
       ENDDO


       d2 = KXCMOP(sum_ene, 1, 6, ADD_OP)
       ene = FLOAT(sum_ene) / vol
       sum_energy = sum_energy + ene
       sum_enersq = sum_enersq + (ene * ene)


       sum_sp = 0
       DO i=0, Nrow-1
          DO j=0, Ncol-1
             sum_sp = sum_sp + spin(i,j)
          ENDDO 
       ENDDO
       d3 = KXCMOP(sum_sp, 1, 6, ADD_OP)
       mag = FLOAT(sum_sp) / vol
       IF(mag .LT. 0.0) THEN
          mag = -mag
       ENDIF

       sum_spins = sum_spins + mag
       sum_spinsq = sum_spinsq + (mag * mag)


       RETURN
       END


C--------------------------------------------------------------

C
C        Take averages, find errors and variances, 
C        and print Out the results
C                                                      

       SUBROUTINE  output
       INCLUDE "exp_wolff.h"
       REAL*8  specific_heat, susceptibility
       REAL*8  err_ener, err_spin, err_size, var_size
       REAL*8  measure

       measure = FLOAT(measurements) * FLOAT(sweeps)

C  ========================  Take averages ==================
       sum_energy = sum_energy / FLOAT(measurements)
       sum_enersq = sum_enersq / FLOAT(measurements)
       sum_spins  = sum_spins  / FLOAT(measurements)
       sum_spinsq = sum_spinsq / FLOAT(measurements)
       sum_size   = sum_size   / measure
       sum_sizesq = sum_sizesq / measure

C  ================  Work out errors and variances.  =====================
       specific_heat  = sum_enersq - (sum_energy * sum_energy)
       susceptibility = sum_spinsq - (sum_spins  * sum_spins)
       var_size       = sum_sizesq - (sum_size   * sum_size)
       err_ener       = sqrt(abs(specific_heat  / FLOAT(measurements)))
       err_spin       = sqrt(abs(susceptibility / FLOAT(measurements)))
       err_size       = sqrt(abs(var_size       / measure ))
       specific_heat  = specific_heat  * vol
       susceptibility = susceptibility * vol

C ====================  Print out results.  =============================
       WRITE(6,10) '*************************************************'
       WRITE(6,*) '      Statistics(Measurements=', measurements, ')'
       WRITE(6,*) '*************************************************'
       WRITE(6,20) ' Average energy = ' ,sum_energy, '+/-' , err_ener 
       WRITE(6,20) ' Average magnet = ' ,sum_spins, '+/-', err_spin 
       WRITE(6,30) ' Specific heat  = ' ,specific_heat        
       WRITE(6,30) ' Susceptibility = ' ,susceptibility 
       WRITE(6,*)  ' Average size   = ' , sum_size, '+/-' ,err_size 
       WRITE(6,*) '*************************************************'
 10    FORMAT(//,2X,A)
 20    FORMAT(4X,A,F10.5,1X,A,F10.5)
 30    FORMAT(4X,A,F10.5)

       RETURN
       END


