
C  ************************************************
C   Subroutine : wolff_update
C  ************************************************

      SUBROUTINE wolff_update
      INCLUDE  "exp_wolff.h"
      INTEGER  LogicByte,d1
      LOGICAL  ANDemptyQ, emptyQ
      EXTERNAL ANDemptyQ, emptyQ

C---Choose a source cell
      CALL set_up 
      
      CALL clearQ
      
      cycle = 0
      IF( procnum  .EQ.  source_pnum ) THEN
         CALL enQ(source_local)
         size=size+1
         spin(source_local(row),source_local(col)) =
     *                   - spin(source_local(row),source_local(col))
         cluster(source_local(row),source_local(col)) = .TRUE.
         new_counter = 1
      ELSE
         new_counter = 0
      ENDIF

      all_emptyQ = .FALSE.
       
      DO WHILE( .NOT. all_emptyQ )
         cycle=cycle+1
         old_counter = new_counter
         new_counter = 0

         CALL next_cycle

         all_emptyQ = emptyQ(Qsize)

         LogicByte = 1
         d1=KXCOMB(all_emptyQ,ANDemptyQ,LogicByte,1,IALNOD,0,ANDTYPE)
      ENDDO

      RETURN
      END



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

      LOGICAL FUNCTION ANDemptyQ(A1,A2,size)
      LOGICAL*1  A1,A2
      INTEGER    size
      
      A1= A1 .AND. A2
c      ANDemptyQ= .TRUE.

      RETURN
      END




C *************************************************************************
C 
C    Subroutine : set_up                                                   
C    Task     : Choose a source cell and destination cell generating random
C               numbers. All processors use same random seed.              
C
C **************************************************************************

      SUBROUTINE  set_up
      INCLUDE  "exp_wolff.h"
      INTEGER  buf(row:col)
      INTEGER  i, j, origin, d1
      INTEGER  temp

      IF( procnum .EQ. 0 ) THEN
        buf(row) = int_ran()
        buf(col) = int_ran()
      ENDIF
      
      origin = 0
      d1 = KXBROD(buf,origin,8,IALNOD,0,RANTYPE)   

      source(row) = buf(row)
      source(col) = buf(col)
   
      temp = MOD(source(col), twidth)
      source_pnum = temp / width
      source_tnum = source(col) / twidth

      source_local(row)  = source(row)
      source_local(col)  = MOD(temp, width) + width * source_tnum

      size = 0

      DO i=0, Nrow-1
         DO j=0, Ncol-1 
            cluster(i,j) = .FALSE.
         ENDDO
      ENDDO

      RETURN
      END




C ***********************************************************
C
C    Subroutine : next_cycle
C
C ***********************************************************

      SUBROUTINE next_cycle
      INCLUDE  "exp_wolff.h"
      INTEGER  i

      CALL clearB

      DO i=0, old_counter-1
         CALL local_expansion
      ENDDO

      CALL communication   

C--- Debug
C      call print_spin

      RETURN
      END
   



C ***************************************************************
C                                                                
C      Subroutine : local_expansion
C                                                               
C ***************************************************************

      SUBROUTINE local_expansion
      INCLUDE  "exp_wolff.h"
      INTEGER dad(row:col), son(row:col)
      INTEGER c(row:val)   
      LOGICAL  lbound, rbound
      EXTERNAL lbound, rbound


      CALL  deQ(dad)

C--- Up
      IF( dad(row) .EQ. 0)  THEN
        son(row) = Nrow-1
      ELSE
        son(row) = dad(row)-1
      ENDIF
      son(col) = dad(col)
      IF( .NOT. cluster(son(row),son(col))
     *    .AND. ( -spin(dad(row),dad(col)) .EQ. 
     *             spin(son(row),son(col)) )
     *    .AND. ( ranmar() .LT. p ) )  THEN
        cluster(son(row),son(col)) = .TRUE.
        spin(son(row),son(col)) = - spin(son(row),son(col)) 
        CALL  enQ(son)
        size=size+1
        new_counter=new_counter+1
      ENDIF

C--- Down
      IF( dad(row) .EQ. Nrow-1 )  THEN
        son(row) = 0
      ELSE
        son(row) = dad(row)+1
      ENDIF
      son(col) = dad(col)
      IF( .NOT. cluster(son(row),son(col)) 
     *    .AND. ( -spin(dad(row),dad(col)) .EQ. 
     *             spin(son(row),son(col)) ) 
     *    .AND. ( ranmar() .LT. p ) )  THEN
          cluster(son(row),son(col)) = .TRUE.
          spin(son(row),son(col)) = - spin(son(row),son(col)) 
          CALL  enQ(son)
          size=size+1
          new_counter=new_counter+1
      ENDIF

C--- Left
      IF( .NOT. lbound(dad(col)) ) THEN
         son(row) = dad(row)
         son(col) = dad(col) - 1
         IF( .NOT. cluster(son(row),son(col))
     *       .AND. ( -spin(dad(row),dad(col)) .EQ.
     *                spin(son(row),son(col)) )
     *       .AND. ( ranmar() .LT. p ) )  THEN
             cluster(son(row),son(col)) = .TRUE.
             spin(son(row),son(col)) = - spin(son(row),son(col)) 
             CALL  enQ(son)
             size=size+1
             new_counter=new_counter+1
	 ENDIF
      ELSE 
         IF( procnum .NE. 0 )  THEN
           c(row) = dad(row)
           c(col) = dad(col) + width - 1
           c(val) = spin(dad(row),dad(col))
           CALL putB(send_Ldata, c, send_Lptr)
         ELSE 
           c(row) = dad(row)
           IF( dad(col) .GT. 0 )  THEN
              c(col) = dad(col)-1
           ELSE
              c(col) = Ncol-1
           ENDIF
           c(val) = spin(dad(row),dad(col))
           CALL putB(send_Ldata, c, send_Lptr)
         ENDIF
      ENDIF 
            
C--- Right    
      IF( .NOT. rbound(dad(col)) ) THEN
         son(row) = dad(row)
         son(col) = dad(col) + 1
         IF( .NOT. cluster(son(row),son(col)) 
     *       .AND. ( -spin(dad(row),dad(col)) .EQ. 
     *                spin(son(row),son(col)) ) 
     *       .AND. ( ranmar() .LT. p ) )  THEN
             cluster(son(row),son(col)) = .TRUE.
             spin(son(row),son(col)) = - spin(son(row),son(col)) 
             CALL enQ(son)
             size=size+1
             new_counter=new_counter+1
	 ENDIF
      ELSE  
         IF( procnum .NE. proc_order(nproc) )  THEN  ! proc_order(npoc) has the # of last processor.
            c(row) = dad(row)
            c(col) = dad(col) - width + 1
            c(val) = spin(dad(row),dad(col))
            CALL putB(send_Rdata, c, send_Rptr)
         ELSE 
            c(row) = dad(row)
            IF( dad(col) .LT. Ncol-1 )  THEN
              c(col) = dad(col) +1
            ELSE
              c(col) = 0
            ENDIF
            c(val) = spin(dad(row),dad(col))
            CALL putB(send_Rdata, c, send_Rptr)
         ENDIF
      ENDIF

      RETURN
      END



C ***************************************************************************
C
C     Subroutine : communication
C 
C ***************************************************************************

      SUBROUTINE  communication
      INCLUDE  "exp_wolff.h"
      INTEGER  sizeB
      EXTERNAL sizeB
      INTEGER  d1,d2,d3,d4
      INTEGER  type1,type2,type3,type4
      INTEGER  sLlength,rRlength,sRlength,rLlength

      send_Lsize = sizeB(send_Lptr)
      send_Rsize = sizeB(send_Rptr)

      type1=EDGE+1
      d1 = KXCHAN( recv_Rsize, intsize, rpnum, type1,  
     &             send_Lsize, intsize, lpnum, type1 )

      sLlength=send_Lsize*3*intsize
      rRlength=recv_Rsize*3*intsize
      type2=EDGE+2    
      d2 = KXCHAN( recv_Rdata, rRlength, rpnum, type2,
     &             send_Ldata, sLlength, lpnum, type2 )


      type3=EDGE+3
      d3 = KXCHAN( recv_Lsize, intsize, lpnum, type3,
     &             send_Rsize, intsize, rpnum, type3 )

      sRlength=send_Rsize*3*intsize
      rLlength=recv_Lsize*3*intsize
      type4=EDGE+4
      d4 = KXCHAN(recv_Ldata, rLlength, lpnum, type4,
     &            send_Rdata, sRlength, rpnum, type4 )
      

      CALL  process_recv_lB
      CALL  process_recv_rB

      RETURN
      END


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

      SUBROUTINE   process_recv_lB
      INCLUDE  "exp_wolff.h"
      INTEGER  c(row:val)
      INTEGER  son(row:col)
      LOGICAL  emptyB
      EXTERNAL emptyB

      CALL set_ptrB(recv_Lptr, recv_lsize-1)

      DO WHILE( .NOT. emptyB(recv_Lptr) ) 
         CALL getB(recv_Ldata,recv_Lptr,c)

         IF( ( .NOT. cluster(c(row),c(col)) ) .AND.
     *       ( -c(val) .EQ. spin(c(row),c(col)) ) .AND.
     *       ( ranmar() .LT. p ) )  THEN
             cluster(c(row),c(col)) = .TRUE.
             size=size+1
             spin(c(row),c(col)) = - spin(c(row),c(col)) 
             son(row) = c(row)
             son(col) = c(col)
             CALL  enQ(son)
             new_counter=new_counter+1
          ENDIF
       ENDDO
 
       RETURN
       END
       


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

      SUBROUTINE   process_recv_rB
      INCLUDE  "exp_wolff.h"
      INTEGER  c(row:val)
      INTEGER  son(row:col)
      LOGICAL  emptyB
      EXTERNAL emptyB

      CALL set_ptrB(recv_Rptr, recv_rsize-1)

      DO WHILE( .NOT. emptyB(recv_Rptr) ) 
         CALL getB(recv_Rdata,recv_Rptr,c)

         IF( ( .NOT. cluster(c(row),c(col)) ) .AND.
     *       ( -c(val) .EQ. spin(c(row),c(col)) ) .AND.
     *       ( ranmar() .LT. p ) )  THEN
             cluster(c(row),c(col)) = .TRUE.
             size=size+1
             spin(c(row),c(col)) = - spin(c(row),c(col)) 
             son(row) = c(row)
             son(col) = c(col)
             CALL  enQ(son)
             new_counter=new_counter+1
          ENDIF
       ENDDO
 
       RETURN
       END






