CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                       
C                THIS SOFTWARE COPYRIGHTED BY CALTECH, 1983, 1986, 1987 
C                                                                       
C                PURE GAUGE SU(3) MONTE CARLO                           
C                THIS CODE HAS A LONG HISTORY:                          
C                    EARLY VERSIONS OF THIS DEVELOPED FOR THE MARK I    
C                    HYPERCUBE BY STEVE OTTO AND PAUL STOLORZ (19       
C                    JON FLOWER AND STEVE OTTO EXTENDED AND IMPROVED THE
C                    CODE FOR THE MARK II MACHINES. (1984-86)           
C                    THE CURRENT INCARNATION IS A "PORT" TO THE NCUBE   
C                    DONE BY CLIVE BAILLIE, ED FELTEN, STEVE OTTO,      
C                    AND DAVID WALKER  (1987)                           
C                                                                       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                       
C               PURE GAUGE SU(3) THEORY IN 3+1 DIMENSIONS               
C               LINK VARIABLE UPDATING DONE VIA CABIBBO-MARINARI HEAT BA
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE          
C     AND DAVID WALKER, SUMMER 1987                                     
C                                                                       
      PROGRAM QCD2                                                      
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      PARAMETER ( MBAR=18*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
      PARAMETER ( LISMAX=(L1*L2*L3*L4)/2+1 )                            
      COMMON /COLCOM/ LISRED(LISMAX),LISBLK(LISMAX),NRED,NBLACK         
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      REAL UBAR(MBAR),SBAR(MBAR)                                        
      COMMON /COMUBR/ UBAR,SBAR                                         
C                                                                       
      INTEGER NUMSTG                                                    
      COMMON /COMFIL/ NUMSTG                                            
C                                                                       
      INTEGER I,J,K,NUMOFF,AVGNUM,FLG                                   
      INTEGER R,T,NUMBAR,DIR,RMIN,RMAX,TMIN(13),TMAX(13)                
      INTEGER R1MAX,R2MAX,R1,R2,R3,TMX,TMN                              
      CHARACTER C                                                       
      INTEGER SWEEPS,PINVAL                                             
C                                                                       
      REAL FC(2)                                                        
      REAL TIME(6)                                                      
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
C                                                                       
      OPEN(5,FILE='LGI5',STATUS='OLD')                                  
      OPEN(6,FILE='LGO6',STATUS='NEW')                                  
      OPEN(8,FILE='LGV',STATUS='NEW')                                   
C                                                                       
      CALL CPUTIM(TIME(3))                                              
      CALL ELAPSE(TIME(4))                                              
C                                                                       
      CALL PRNSET(42)                                                   
C                                                                       
      WRITE(8,201)                                                      
201   FORMAT(1X,
     * 'OUTPUT FOR PERFECT CLUB BENCHMARK: QCD INSTRUMENTED VERSION'/
     *       1X,'$Revision: 1.A.1 $ $Author: psharma $'/)                   
      WRITE(8,202)                                                      
202   FORMAT(1X,'VALIDATION PARAMETERS:'/)                              
C                                                                       
      WRITE(6,*) 'PURE GAUGE QCD USING PSEUDO HEAT BATH'                
      WRITE(6,*) '-------------------------------------'                
      WRITE(6,*) ' '                                                    
      WRITE(6,*) 'ENTER COMMAND (? FOR HELP)'                           
      WRITE(6,*) ' '                                                    
   10 READ (5,200) C                                                    
C                                                                       
      IF (C.EQ.'R') THEN                                                
        CALL SETSED()                                                   
      ELSEIF (C.EQ.'D') THEN                                            
        CALL DMPSED()                                                   
      ELSEIF (C.EQ.'S') THEN                                            
        CALL SETPAR()                                                   
      ELSEIF (C.EQ.'P') THEN                                            
        CALL PROJEC()                                                   
      ELSEIF (C.EQ.'U') THEN                                            
        WRITE(6,*) 'ENTER NUMBER OF TEMPORAL LINK HITS TO AVERAGE OVER' 
        READ (5,*) AVGNUM                                               
        WRITE(6,*) AVGNUM                                               
C       CLOCK = SECNDS(0)                                               
        CALL CPUTIM(TIME(1))                                            
        CALL LINKBR(AVGNUM,3,UBAR)                                      
        CALL CPUTIM(TIME(2))                                            
C       WRITE(6,101) SECNDS(0)-CLOCK                                    
        WRITE(6,101) TIME(2)-TIME(1)                                    
        WRITE(6,*) 'UBARS MADE.'                                        
      ELSEIF (C.EQ.'M') THEN                                            
        WRITE(6,*) 'ENTER DIR (0<->X,1<->Y,2<->Z), RMIN, RMAX'          
        READ (5,*) DIR,RMIN,RMAX                                        
        WRITE(6,*) DIR,RMIN,RMAX                                        
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
        DO 1 R=RMIN,RMAX                                                
          WRITE(6,*) 'FOR R = ',R,' ENTER TMIN, TMAX '                  
          READ (5,*) TMIN(R+1),TMAX(R+1)                                
          WRITE(6,*) TMIN(R+1),TMAX(R+1)                                
    1   CONTINUE                                                        
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
C                                                                       
        DO 2 R=1,13                                                     
          DO 2 T=1,13                                                   
            WHICH(R,T) = 0                                              
    2   CONTINUE                                                        
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
C                                                                       
        DO 3 R=RMIN,RMAX                                                
          DO 3 T=1,13                                                   
            IF((T.GE.TMIN(R+1)).AND.(T.LE.TMAX(R+1))) THEN              
              WHICH(R+1,T+1) = 1                                        
            ENDIF                                                       
    3   CONTINUE                                                        
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
C                                                                       
        WRITE(6,*) 'ENTER NUMBER OF SPATIAL LINK HITS TO AVERAGE OVER'  
        READ (5,*) AVGNUM                                               
        WRITE(6,*) AVGNUM                                               
        WRITE(6,*) 'MEASURING WILSON (ON-AXIS) LOOPS IN DIRECTION ',DIR 
C       CLOCK = SECNDS(0)                                               
        CALL CPUTIM(TIME(1))                                            
        CALL MEASUR(DIR,RMIN,RMAX,AVGNUM)                               
        CALL CPUTIM(TIME(2))                                            
        WRITE(6,102) TIME(2)-TIME(1)                                    
C       WRITE(6,102) SECNDS(0)-CLOCK                                    
        FC(1) = NUMSTG                                                  
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
        DO 4 R=RMIN,RMAX                                                
          DO 4 T=TMIN(R+1),TMAX(R+1)                                    
            FC(2) = WLOOP(R+1,T+1)                                      
            FC(2) = FC(2)/FC(1)                                         
            WRITE(6,*) 'W',R,'X',T,':',FC(2)                            
            WRITE(8,*) 'W',R,'X',T,':',FC(2)                            
C                                                                       
    4   CONTINUE                                                        
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
C                                                                       
        IVAL=0                                                          
        IF (ABS((FC(2)-.633259)/.633259).GE.1E-2) IVAL=1                
C                                                                       
      ELSEIF (C.EQ.'I') THEN                                            
C        NUMOFF = 0                                                     
      ELSEIF (C.EQ.'A') THEN                                            
        WRITE(6,*) 'ENTER R1,R2,R3,TMN,TMX'                             
        READ (5,*) R1,R2,R3,TMN,TMX                                     
        WRITE(6,*) R1,R2,R3,TMN,TMX                                     
        WRITE(6,*) 'MEASURING BARYON LOOPS ...'                         
C       CLOCK = SECNDS(0)                                               
        CALL CPUTIM(TIME(1))                                            
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
        DO 5 T=TMN,TMX                                                  
          CALL QQQMEA(R1,R2,R3,T)                                       
    5   CONTINUE                                                        
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
        CALL CPUTIM(TIME(2))                                            
        WRITE(6,103) TIME(2)-TIME(1)                                    
C       WRITE(6,103) SECNDS(0)-CLOCK                                    
C                                                                       
         IF(IVAL.EQ.0) THEN                                             
           WRITE(8,203)                                                 
203        FORMAT(//1X,'RESULTS FOR THIS RUN ARE:   VALID')             
         ELSE                                                           
           WRITE(8,204)                                                 
204        FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')             
         ENDIF                                                          
C                                                                       
      ELSEIF (C.EQ.'E') THEN                                            
        WRITE(6,*) 'ENTER DIR (0<->X,1<->Y,2<->Z)'                      
        READ (5,*) DIR                                                  
        WRITE(6,*) DIR                                                  
        WRITE(6,*) 'ENTER NUMBER OF SPATIAL LINK HITS TO AVERAGE OVER'  
        READ (5,*) NUMBAR                                               
        WRITE(6,*) NUMBAR                                               
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
        DO 6 I=1,13                                                     
          DO 6 J=1,13                                                   
            DO 6 K=1,6                                                  
    6          RES(I,J,K)=0.0                                           
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
        WRITE(6,*) 'ENTER R1MAX AND R2MAX :'                            
        READ (5,*) R1MAX,R2MAX                                          
        WRITE(6,*) R1MAX,R2MAX                                          
        WRITE(6,*) 'NOW ENTER THE QUADRANT TABLE OF TMINS ;'            
        WRITE(6,*) R1MAX,'ROWS &',R2MAX,'COLS (0 ==> OMIT)'             
        READ(5,*)((WHICH(I,J), J=1,R2MAX),I=1,R1MAX)                    
        WRITE(6,*) 'MEASURING LOOPS IN DIRECTION ',DIR                  
C                                                                       
C       CLOCK = SECNDS(0)                                               
        CALL CPUTIM(TIME(1))                                            
        CALL ROTMEA(DIR,NUMBAR,R1MAX,R2MAX)                             
        CALL CPUTIM(TIME(2))                                            
        WRITE(6,105) TIME(2)-TIME(1)                                    
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
C       WRITE(6,105) SECNDS(0)-CLOCK                                    
        DO 7 I=1,R1MAX                                                  
          DO 7 J=1,R2MAX                                                
            IF(WHICH(I,J).NE.0) THEN                                    
              TMX = MIN(12,WHICH(I,J)+5)                                
              DO 8 T=WHICH(I,J),TMX-1                                   
                WRITE(6,*) 'R (',I+1,',',J+1,') X',T,':',               
     .                   RES(I,J,T-WHICH(I,J)+1)/REAL(NUMSTG)           
                WRITE(8,*) 'R (',I+1,',',J+1,') X',T,':',               
     .                   RES(I,J,T-WHICH(I,J)+1)/REAL(NUMSTG)           
    8           CONTINUE                                                
            ENDIF                                                       
    7     CONTINUE                                                      
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
C                                                                       
        IF (ABS((RES(1,1,1-WHICH(1,1)+1)/REAL(NUMSTG)                   
     .          -.434777)/.434777).GE.1E-2) IVAL=1                      
        IF (ABS((RES(1,1,2-WHICH(1,1)+1)/REAL(NUMSTG)                   
     .          -.299319)/.299319).GE.1E-2) IVAL=1                      
        IF (ABS((RES(1,1,3-WHICH(1,1)+1)/REAL(NUMSTG)                   
     .          -.205511)/.205511).GE.1E-2) IVAL=1                      
        IF (ABS((RES(1,1,4-WHICH(1,1)+1)/REAL(NUMSTG)                   
     .         -.136702)/.136702).GE.1E-2) IVAL=1                       
        IF (ABS((RES(1,1,5-WHICH(1,1)+1)/REAL(NUMSTG)                   
     .      -9.14717E-02)/9.14717E-02).GE.1E-2) IVAL=1                  
C                                                                       
      ELSEIF (C.EQ.'C') THEN                                            
        WRITE(6,*) '0 <-> COLD START, 1 <-> RESTART'                    
        READ (5,*) FLG                                                  
        WRITE(6,*) FLG                                                  
        CALL INIT(FLG)                                                  
        WRITE(6,*) 'INITIALIZED SUCCESSFULLY.'                          
      ELSEIF (C.EQ.'B') THEN                                            
        CALL LATTOT()                                                   
        WRITE(6,*) 'BACKED UP SUCCESSFULLY'                             
      ELSEIF (C.EQ.'F') THEN                                            
        WRITE(6,*) 'INPUT # OF SWEEPS AND PROJECT INTERVAL'             
        READ (5,*) SWEEPS,PINVAL                                        
        WRITE(6,*) SWEEPS,PINVAL                                        
        WRITE(6,*) 'RUNNING...'                                         
C       CLOCK = SECNDS(0)                                               
        CALL CPUTIM(TIME(1))                                            
        CALL UPDATE(SWEEPS,PINVAL)                                      
        CALL CPUTIM(TIME(2))                                            
        WRITE(6,104) TIME(2)-TIME(1)                                    
C       WRITE(6,104) SECNDS(0)-CLOCK                                    
        WRITE(6,*) 'LATTICE UPDATED'                                    
      ELSEIF (C.EQ.'?') THEN                                            
        WRITE(6,*) 'COMMANDS (AND THEIR PARAMETERS) ARE:'               
        WRITE(6,*) '===================================='               
        WRITE(6,*) 'C :             INITIALIZE (COLD OR READ)'          
        WRITE(6,*) 'D :             UPDATES (SWEEPS, PINVAL)'           
        WRITE(6,*) 'B :             BACK UP LATTICE (FILENAME)'         
        WRITE(6,*) 'S :             SET PARAMETERS'                     
        WRITE(6,*) 'R :             RESET RANDOM NUMBER SEEDS'          
        WRITE(6,*) 'D :             DUMP CURRENT SEED TO FILE'          
        WRITE(6,*) 'P :             PROJECT MATRICES ONTO SU(3)'        
        WRITE(6,*) 'M :             MEASURE ON AXIS LOOPS'              
        WRITE(6,*) 'E :             MEASURE OFF AXIS LOOP'              
        WRITE(6,*) 'A :             MEASURE BARYON LOOPS'               
        WRITE(6,*) 'U :             MAKE U BARS'                        
        WRITE(6,*) 'Q :             QUIT'                               
      ELSEIF (C.EQ.'Q') THEN                                            
        WRITE(6,*) 'FINISHED.'                                          
C                                                                       
        CALL CPUTIM(TIME(5))                                            
        CALL ELAPSE(TIME(6))                                            
C                                                                       
      WRITE(8,*)                                                        
      WRITE(8,*) 'ELAPSED CPU TIME IN SECONDS: ',TIME(5)-TIME(3)        
C MFLOPS = MFLOP 1 CPU CRAY X-MP/CPU TIME                               
      WRITE(8,*) 'MFLOPS RATE: ', 259.550/(TIME(5)-TIME(3))             
      WRITE(8,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',                 
     *            MOD(TIME(6)-TIME(4)+1000000.0,1000000.0)              
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
C                                                                       
        STOP                                                            
      ELSE                                                              
        WRITE(6,*) 'UNKNOWN COMMAND (? FOR HELP)'                       
      ENDIF                                                             
      GOTO 10                                                           
C                                                                       
  101 FORMAT(/5X,'Time for LINKBR = ',F11.4,' seconds')                 
  102 FORMAT(/5X,'Time for MEASUR = ',F11.4,' seconds')                 
  103 FORMAT(/5X,'Time for QQQMEA = ',F11.4,' seconds')                 
  104 FORMAT(/5X,'Time for UPDATE  = ',F11.4,' seconds')                
  105 FORMAT(/5X,'Time for ROTMEA = ',F11.4,' seconds')                 
  200 FORMAT(A1)                                                        
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
      END                                                               
                                                                        
                                                                        
                                                                        
                                                                        
C                                                                       
      SUBROUTINE CFILL(A,B)                                             
      REAL A(*),B(*)                                                    
CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (3)                                                  
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
C                                                                       
      DO 1 I=1,12                                                       
        A(I) = B(I)                                                     
    1   CONTINUE                                                        
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
      A(13) =  A(3)*A(11)-A(4)*A(12)-A(5)*A(9) +A(6)*A(10)              
      A(14) = -A(4)*A(11)-A(3)*A(12)+A(6)*A(9) +A(5)*A(10)              
      A(15) =  A(5)*A(7) -A(6)*A(8) -A(1)*A(11)+A(2)*A(12)              
      A(16) = -A(6)*A(7) -A(5)*A(8) +A(2)*A(11)+A(1)*A(12)              
      A(17) =  A(1)*A(9) -A(2)*A(10)-A(3)*A(7) +A(4)*A(8)               
      A(18) = -A(2)*A(9) -A(1)*A(10)+A(4)*A(7) +A(3)*A(8)               
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C     HEAT BATH - USES CABIBBO-MARINARI ALGORITHM INSTEAD OF            
C                 METROPOLIS                                            
C                 3 SUBGROUP VERSION                                    
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE          
C     AND DAVID WALKER, SUMMER 1987                                     
C                                                                       
      SUBROUTINE CHOOS(SST,PU)                                          
      REAL SST(*)                                                       
      REAL PU(*)                                                        
C                                                                       
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL R0,R1,R2,R3,R,A0,A1,A2,A3,XX,EFAC,RFAC,NORM                  
      REAL PART(18),PLAQ(18),MAT(18),DIS,B0,B1,B2,B3,TEMP(18)           
      INTEGER I,J                                                       
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
C                                                                       
      DO 1 I=1,18                                                       
        PART(I) = 0.0                                                   
        MAT(I)  = PU(I)                                                 
    1   CONTINUE                                                        
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
C                                                                       
      DO 2 J=1,6                                                        
        DO 2 I=1,18                                                     
          PART(I) = PART(I)+SST(I+18*(J-1))                             
    2     CONTINUE                                                      
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
C                                                                       
      DO 3 I=1,3                                                        
        CALL MULT(MAT,PART,PLAQ)                                        
        IF(I.EQ.1) THEN                                                 
          R0 = 0.5*(PLAQ(1)+PLAQ(9))                                    
          R1 = 0.5*(PLAQ(2)-PLAQ(10))                                   
          R2 = 0.5*(PLAQ(3)-PLAQ(7))                                    
          R3 = 0.5*(PLAQ(4)+PLAQ(8))                                    
        ELSEIF(I.EQ.2) THEN                                             
          R0 = 0.5*(PLAQ(9)+PLAQ(17))                                   
          R1 = 0.5*(PLAQ(10)-PLAQ(18))                                  
          R2 = 0.5*(PLAQ(11)-PLAQ(15))                                  
          R3 = 0.5*(PLAQ(12)+PLAQ(16))                                  
        ELSE                                                            
          R0 = 0.5*(PLAQ(1)+PLAQ(17))                                   
          R1 = 0.5*(PLAQ(2)-PLAQ(18))                                   
          R2 = 0.5*(PLAQ(5)-PLAQ(13))                                   
          R3 = 0.5*(PLAQ(6)+PLAQ(14))                                   
        ENDIF                                                           
        R    = R0*R0+R1*R1+R2*R2+R3*R3                                  
        R    = SQRT(R)                                                  
        R0   = R0/R                                                     
        R1   = R1/R                                                     
        R2   = R2/R                                                     
        R3   = R3/R                                                     
        XX   = BETA*R*2.0/3.0                                           
        EFAC = EXP(-2.0*XX)                                             
        RFAC = 1.0-EFAC                                                 
C                                                                       
    4   A0   = 1.0+LOG(EFAC+RFAC*PRANF(0))/XX                           
        XTEMP= PRANF(0)                                                 
        IF((XTEMP*XTEMP).GT.(1.0-A0*A0)) GOTO 4                         
        NORM = SQRT(1.0-A0*A0)                                          
C                                                                       
    5   A1   = 2.0*PRANF(0)-1.0                                         
        A2   = 2.0*PRANF(0)-1.0                                         
        A3   = 2.0*PRANF(0)-1.0                                         
        DIS  = A1*A1+A2*A2+A3*A3                                        
        IF(DIS.GT.1.0) GOTO 5                                           
C                                                                       
        NORM = NORM/SQRT(DIS)                                           
        A1   = A1*NORM                                                  
        A2   = A2*NORM                                                  
        A3   = A3*NORM                                                  
        B0   = A0*R0+A1*R1+A2*R2+A3*R3                                  
        B1   =  -A0*R1+A1*R0+A2*R3-A3*R2                                
        B2   = A2*R0+A3*R1-A0*R2-A1*R3                                  
        B3   =  -A2*R1+A3*R0-A0*R3+A1*R2                                
        IF(I.EQ.1) THEN                                                 
          TEMP(1)  = B0*MAT(1)-B1*MAT(2)-B2*MAT(3)-B3*MAT(4)            
          TEMP(2)  = B1*MAT(1)+B0*MAT(2)+B3*MAT(3)-B2*MAT(4)            
          TEMP(3)  = B2*MAT(1)-B3*MAT(2)+B0*MAT(3)+B1*MAT(4)            
          TEMP(4)  = B3*MAT(1)+B2*MAT(2)-B1*MAT(3)+B0*MAT(4)            
          TEMP(7)  = B0*MAT(7)-B1*MAT(8)-B2*MAT(9)-B3*MAT(10)           
          TEMP(8)  = B1*MAT(7)+B0*MAT(8)+B3*MAT(9)-B2*MAT(10)           
          TEMP(9)  = B2*MAT(7)-B3*MAT(8)+B0*MAT(9)+B1*MAT(10)           
          TEMP(10) = B3*MAT(7)+B2*MAT(8)-B1*MAT(9)+B0*MAT(10)           
          TEMP(13) = B0*MAT(13)-B1*MAT(14)-B2*MAT(15)-B3*MAT(16)        
          TEMP(14) = B1*MAT(13)+B0*MAT(14)+B3*MAT(15)-B2*MAT(16)        
          TEMP(15) = B2*MAT(13)-B3*MAT(14)+B0*MAT(15)+B1*MAT(16)        
          TEMP(16) = B3*MAT(13)+B2*MAT(14)-B1*MAT(15)+B0*MAT(16)        
          MAT(1)   = TEMP(1)                                            
          MAT(2)   = TEMP(2)                                            
          MAT(3)   = TEMP(3)                                            
          MAT(4)   = TEMP(4)                                            
          MAT(7)   = TEMP(7)                                            
          MAT(8)   = TEMP(8)                                            
          MAT(9)   = TEMP(9)                                            
          MAT(10)  = TEMP(10)                                           
          MAT(13)  = TEMP(13)                                           
          MAT(14)  = TEMP(14)                                           
          MAT(15)  = TEMP(15)                                           
          MAT(16)  = TEMP(16)                                           
        ELSEIF(I.EQ.2) THEN                                             
          TEMP(3)  = B0*MAT(3)-B1*MAT(4)-B2*MAT(5)-B3*MAT(6)            
          TEMP(4)  = B1*MAT(3)+B0*MAT(4)+B3*MAT(5)-B2*MAT(6)            
          TEMP(5)  = B2*MAT(3)-B3*MAT(4)+B0*MAT(5)+B1*MAT(6)            
          TEMP(6)  = B3*MAT(3)+B2*MAT(4)-B1*MAT(5)+B0*MAT(6)            
          TEMP(9)  = B0*MAT(9)-B1*MAT(10)-B2*MAT(11)-B3*MAT(12)         
          TEMP(10) = B1*MAT(9)+B0*MAT(10)+B3*MAT(11)-B2*MAT(12)         
          TEMP(11) = B2*MAT(9)-B3*MAT(10)+B0*MAT(11)+B1*MAT(12)         
          TEMP(12) = B3*MAT(9)+B2*MAT(10)-B1*MAT(11)+B0*MAT(12)         
          TEMP(15) = B0*MAT(15)-B1*MAT(16)-B2*MAT(17)-B3*MAT(18)        
          TEMP(16) = B1*MAT(15)+B0*MAT(16)+B3*MAT(17)-B2*MAT(18)        
          TEMP(17) = B2*MAT(15)-B3*MAT(16)+B0*MAT(17)+B1*MAT(18)        
          TEMP(18) = B3*MAT(15)+B2*MAT(16)-B1*MAT(17)+B0*MAT(18)        
          MAT(3)   = TEMP(3)                                            
          MAT(4)   = TEMP(4)                                            
          MAT(5)   = TEMP(5)                                            
          MAT(6)   = TEMP(6)                                            
          MAT(9)   = TEMP(9)                                            
          MAT(10)  = TEMP(10)                                           
          MAT(11)  = TEMP(11)                                           
          MAT(12)  = TEMP(12)                                           
          MAT(15)  = TEMP(15)                                           
          MAT(16)  = TEMP(16)                                           
          MAT(17)  = TEMP(17)                                           
          MAT(18)  = TEMP(18)                                           
        ELSE                                                            
          TEMP(1)  = B0*MAT(1)-B1*MAT(2)-B2*MAT(5)-B3*MAT(6)            
          TEMP(2)  = B1*MAT(1)+B0*MAT(2)+B3*MAT(5)-B2*MAT(6)            
          TEMP(5)  = B2*MAT(1)-B3*MAT(2)+B0*MAT(5)+B1*MAT(6)            
          TEMP(6)  = B3*MAT(1)+B2*MAT(2)-B1*MAT(5)+B0*MAT(6)            
          TEMP(7)  = B0*MAT(7)-B1*MAT(8)-B2*MAT(11)-B3*MAT(12)          
          TEMP(8)  = B1*MAT(7)+B0*MAT(8)+B3*MAT(11)-B2*MAT(12)          
          TEMP(11) = B2*MAT(7)-B3*MAT(8)+B0*MAT(11)+B1*MAT(12)          
          TEMP(12) = B3*MAT(7)+B2*MAT(8)-B1*MAT(11)+B0*MAT(12)          
          MAT(1)   = TEMP(1)                                            
          MAT(2)   = TEMP(2)                                            
          MAT(5)   = TEMP(5)                                            
          MAT(6)   = TEMP(6)                                            
          MAT(7)   = TEMP(7)                                            
          MAT(8)   = TEMP(8)                                            
          MAT(11)  = TEMP(11)                                           
          MAT(12)  = TEMP(12)                                           
        ENDIF                                                           
    3   CONTINUE                                                        
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
C                                                                       
      CALL CFILL(PU,MAT)                                                
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE CPYMAT(OUTMAT,INMAT,N)                                 
      REAL INMAT(*), OUTMAT(*)                                          
      INTEGER N                                                         
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
C                                                                       
      DO 1 I=1,N                                                        
    1   OUTMAT(I) = INMAT(I)                                            
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE ICPYMT(OUTMAT,INMAT,N)                                 
      INTEGER INMAT(*), OUTMAT(*), N                                    
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
                                                                        
C     ADDED BY LPOINTER 4-3-89                                          
                                                                        
      DO 1 I=1,N                                                        
   1     OUTMAT(I) = INMAT(I)                                           
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
                                                                        
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE DMPSED()                                               
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
C                                                                       
      OPEN(UNIT=1,FILE='SEED',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')     
      REWIND(1)                                                         
      INO = IGETSD(0)                                                   
      WRITE (1,10) INO                                                  
   10 FORMAT(I12)                                                       
      WRITE(6,*) 'O.K. ...... DONE IT !'                                
      CLOSE(1)                                                          
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
      END                                                               
C                                                                       
CCCCCCCCCCCCCCCCCCCCCCCCC END OF ROUTINE DMPSED CCCCCCCCCCCCCCCCCCCC    
C                                                                       
                                                                        
                                                                        
                                                                        
                                                                        
       SUBROUTINE GETPOS(LATT,COORD,INDEX)                              
       INTEGER COORD(4), LATT(4), INDEX                                 
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
C                                                                       
       J = INDEX                                                        
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
       DO 1 I=1,4                                                       
          COORD(I) = MOD(J,LATT(I))                                     
    1     J = J/LATT(I)                                                 
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
C                                                                       
       RETURN                                                           
       END                                                              
                                                                        
      SUBROUTINE HIT(COORD,ORIGIN,DIR)                                  
      INTEGER COORD(4),DIR,ORIGIN                                       
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MBAR=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MBAR)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      INTEGER OTHEND,SITE,CNT                                           
      INTEGER PU                                                        
      REAL ST(18,6)                                                     
      INTEGER PTEMP(24)                                                 
C                                                                       
      INTEGER PSTR(24,4)                                                
      SAVE PSTR                                                         
C                                                                       
      DATA PSTR/ 2,5,6,14,6,5,2,14,3,5,7,14,7,5,3,14,4,5,8,14,8,5,4,14, 
     .           1,6,5,14,5,6,1,14,3,6,7,14,7,6,3,14,4,6,8,14,8,6,4,14, 
     .           1,7,5,14,5,7,1,14,2,7,6,14,6,7,2,14,4,7,8,14,8,7,4,14, 
     .           1,8,5,14,5,8,1,14,2,8,6,14,6,8,2,14,3,8,7,14,7,8,3,14 /
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
C                                                                       
C      CALL CPYMAT(PTEMP,PSTR(1,DIR+1),24) LPOINTER 4-3-89              
      CALL ICPYMT(PTEMP,PSTR(1,DIR+1),24)                               
C                                                                       
      IF(COORD(DIR+1).EQ.LATT1(DIR+1)) THEN                             
        OTHEND = 0                                                      
      ELSE                                                              
        OTHEND = COORD(DIR+1)+1                                         
      ENDIF                                                             
C                                                                       
      COORD(DIR+1) = OTHEND                                             
      SITE = COORD(1)*MOV(1)+COORD(2)*MOV(2)+COORD(3)*MOV(3)+           
     .           COORD(4)*MOV(4)                                        
      CNT = 1                                                           
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
      DO 1 J=1,6                                                        
        COORD(DIR+1) = OTHEND                                           
        CALL SYSLOP(COORD,SITE,PTEMP(CNT),ST(1,J),NCHAR)                
        CNT = NCHAR+CNT                                                 
    1 CONTINUE                                                          
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
C                                                                       
      PU = ORIGIN+ROT(DIR+1)                                            
      CALL CHOOS(ST(1,1),U1(PU+1))                                      
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
                                                                        
                                                                        
      SUBROUTINE INIT(FLAG)                                             
C                                                                       
      INTEGER NUMSTG                                                    
      COMMON /COMFIL/ NUMSTG                                            
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      INTEGER UMAX                                                      
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3, UMAX=72*L1*L2*L3*L4 )               
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      REAL FACFI,FACIF                                                  
      COMMON /COMFAC/ FACFI,FACIF                                       
C                                                                       
      INTEGER FLAG,BIGNUM,MOVX,BUFF,REM,BUFNUM                          
      INTEGER PU                                                        
      CHARACTER*5 BNAME                                                 
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
C                                                                       
      MOVX = 72                                                         
      BUFF = 256                                                        
      BIGNUM = 4*LATT(1)*LATT(2)*LATT(3)*LATT(4)                        
C                                                                       
      IF (FLAG.EQ.1) THEN                                               
        WRITE(6,*) 'BACKUP FILE NAME ? : '                              
        READ(5,100) BNAME                                               
        OPEN(UNIT=1,FILE = BNAME,ACCESS = 'SEQUENTIAL',STATUS='UNKNOWN',
     .       FORM='UNFORMATTED')                                        
        REWIND(1)                                                       
        BIGNUM = MOVX*NUMSTG                                            
        BUFNUM = BIGNUM/BUFF                                            
        REM = MOD(BIGNUM, BUFF)                                         
        WRITE(6,*) 'BACKUP FILE NAME = ',BNAME,'BUFNUM = ',BUFNUM,      
     .  'REM = ',REM                                                    
        PU=1                                                            
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
        DO 1 I=0,BUFNUM-1                                               
          READ (1) (U1(J),J=PU,PU+BUFF)                                 
          PU = PU + BUFF                                                
          WRITE(6,*) 'B'                                                
    1     CONTINUE                                                      
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
        WRITE(6,*) 'E'                                                  
        CLOSE(1)                                                        
      ELSE                                                              
        PU=1                                                            
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
        DO 2 I=0,BIGNUM-1                                               
          U1(PU)    = 1.0                                               
          U1(PU+1)  = 0.0                                               
          U1(PU+2)  = 0.0                                               
          U1(PU+3)  = 0.0                                               
          U1(PU+4)  = 0.0                                               
          U1(PU+5)  = 0.0                                               
          U1(PU+6)  = 0.0                                               
          U1(PU+7)  = 0.0                                               
          U1(PU+8)  = 1.0                                               
          U1(PU+9)  = 0.0                                               
          U1(PU+10) = 0.0                                               
          U1(PU+11) = 0.0                                               
          U1(PU+12) = 0.0                                               
          U1(PU+13) = 0.0                                               
          U1(PU+14) = 0.0                                               
          U1(PU+15) = 0.0                                               
          U1(PU+16) = 1.0                                               
          U1(PU+17) = 0.0                                               
          PU=PU+ROT(2)                                                  
    2     CONTINUE                                                      
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
      ENDIF                                                             
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
C                                                                       
      RETURN                                                            
  100 FORMAT(A)                                                         
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      END                                                               
                                                                        
C                                                                       
C     THE ROUTINES USED TO INITIALIZE AND DUMP THE LATTICE.             
C                                                                       
C     CONVERTED FORM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, SUMMER 1987                                     
C                                                                       
      SUBROUTINE LATTOT()                                               
C                                                                       
      INTEGER UMAX                                                      
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      INTEGER NUMSTG                                                    
      COMMON /COMFIL/ NUMSTG                                            
C                                                                       
      INTEGER BIGNUM                                                    
      INTEGER BUFNUM,REM                                                
      INTEGER PU                                                        
      INTEGER BUFF,MOVX                                                 
      CHARACTER*5 BNAME                                                 
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
C                                                                       
      WRITE(6,*) 'ENTER 5-LETTER BACKUP FILE NAME ? : '                 
      READ(5,100) BNAME                                                 
      WRITE(6,*) BNAME                                                  
      OPEN(UNIT=1,FILE = BNAME,ACCESS = 'SEQUENTIAL',STATUS='UNKNOWN')  
C     OPEN(UNIT=1,FILE = BNAME,ACCESS = 'SEQUENTIAL',STATUS='UNKNOWN',  
C    .FORM='UNFORMATTED')                                               
      REWIND(1)                                                         
      MOVX = 72                                                         
      BUFF = 256                                                        
      BIGNUM = MOVX*NUMSTG                                              
      BUFNUM = BIGNUM/BUFF                                              
C      REM = MOD(BIGNUM, BUFF)                                          
      WRITE(6,*) 'BACKUP FILE NAME = ',BNAME,'BUFNUM = ',BUFNUM,'REM = '
      PU = 1                                                            
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
      DO 2 I=0,BUFNUM-1                                                 
C       WRITE (1) (U1(J),J=PU,PU+BUFF)                                  
        WRITE (1,101) (U1(J),J=PU,PU+BUFF)                              
  101   FORMAT(9F8.3)                                                   
        PU = PU + BUFF                                                  
        WRITE(6,*) 'B'                                                  
    2   CONTINUE                                                        
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
      WRITE(6,*) 'E'                                                    
      CLOSE(1)                                                          
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
      RETURN                                                            
  100 FORMAT(A)                                                         
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
      END                                                               
                                                                        
C                                                                       
C     THIS ROUTINE MAKES BARRED LEGS ONLY IN THE ARRAY LEGS             
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE LEGMAK(T,TT)                                           
      INTEGER TT,T                                                      
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      INTEGER STRB(26)                                                  
C                                                                       
      INTEGER SITE,COORD(4),PLEGS                                       
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
C                                                                       
      DO 1 I=1,T                                                        
    1   STRB(I) = 18                                                    
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
C                                                                       
      DO 2 I=T+1,2*T+1                                                  
    2   STRB(I) = 30                                                    
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
C                                                                       
      STRB(2*T+1) = 14                                                  
C                                                                       
      COORD(4) = TT                                                     
      SITE  = TT*MOV(4)                                                 
      PLEGS = 1                                                         
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
      DO 3 I=0,LATT(3)-1                                                
        COORD(3) = I                                                    
        DO 3 J=0,LATT(2)-1                                              
          COORD(2) = J                                                  
          DO 3 K=0,LATT(1)-1                                            
            COORD(1) = K                                                
            CALL SYSLOP(COORD,SITE,STRB,LEGS(PLEGS),NCHAR)              
            PLEGS = PLEGS+18                                            
            SITE  = SITE+MOV(1)                                         
    3       CONTINUE                                                    
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C     LINKBR MULTI-HITS AND AVERAGES U MATRICES IN ANY DIRECTION        
C     AND PUTS THE RESULTS IN THE ARRAY POINTED TO BY "WHERE".          
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE LINKBR(N,DIR,WHERE)                                    
      REAL WHERE(*)                                                     
      INTEGER N,DIR                                                     
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      INTEGER COORD(4),OTHEND,SITE,ENDSIT,PT,CNT                        
      INTEGER PU                                                        
      INTEGER PSTR(24)                                                  
      REAL AVMAT(18),MAT2(18),ST(18,6)                                  
C                                                                       
      INTEGER ENV(24,4)                                                 
      SAVE ENV                                                          
C                                                                       
      DATA ENV/ 2,5,6,14,6,5,2,14,3,5,7,14,7,5,3,14,4,5,8,14,8,5,4,14,  
     .          1,6,5,14,5,6,1,14,3,6,7,14,7,6,3,14,4,6,8,14,8,6,4,14,  
     .          1,7,5,14,5,7,1,14,2,7,6,14,6,7,2,14,4,7,8,14,8,7,4,14,  
     .          1,8,5,14,5,8,1,14,2,8,6,14,6,8,2,14,3,8,7,14,7,8,3,14 / 
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
C                                                                       
      PT = 1                                                            
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
      DO 1 ICORD4=0,LATT(4)-1                                           
        COORD(4) = ICORD4                                               
        DO 1 ICORD3=0,LATT(3)-1                                         
          COORD(3) = ICORD3                                             
          DO 1 ICORD2=0,LATT(2)-1                                       
            COORD(2) = ICORD2                                           
            DO 1 ICORD1=0,LATT(1)-1                                     
              COORD(1) = ICORD1                                         
              CNT = 1                                                   
              SITE = ICORD1*MOV(1)+ICORD2*MOV(2)+ICORD3*MOV(3)          
     .                                            +ICORD4*MOV(4)        
              PU = SITE+ROT(DIR+1)                                      
              CALL CPYMAT(AVMAT,U1(PU+1),18)                            
              IF(N.GT.1) THEN                                           
                DO 2 J=1,18                                             
    2             MAT2(J) = U1(PU+J)                                    
                IF(COORD(DIR+1).EQ.LATT1(DIR+1)) THEN                   
                  OTHEND = 0                                            
                  ENDSIT = SITE - LATT1(DIR+1)*MOV(DIR+1)               
                ELSE                                                    
                  OTHEND = COORD(DIR+1) +1                              
                  ENDSIT = SITE + MOV(DIR+1)                            
                ENDIF                                                   
C                CALL CPYMAT(PSTR,ENV(1,DIR+1),24) LPOINTER 4-3-89      
                CALL ICPYMT(PSTR,ENV(1,DIR+1),24)                       
C                                                                       
                DO 3 J=1,6                                              
                  COORD(DIR+1) = OTHEND                                 
                  CALL SYSLOP(COORD,ENDSIT,PSTR(CNT),ST(1,J),NCHAR)     
    3             CNT = CNT+NCHAR                                       
                DO 4 J=1,N-1                                            
                  CALL CHOOS(ST(1,1),MAT2)                              
                  DO 5 I=1,18                                           
    5               AVMAT(I) = AVMAT(I)+MAT2(I)                         
    4             CONTINUE                                              
                DO 6 I=1,18                                             
    6             AVMAT(I) = AVMAT(I)/N                                 
              ENDIF                                                     
                DO 7 K=1,18                                             
    7            WHERE(PT+K-1) = AVMAT(K)                               
              PT = PT + 18                                              
    1         CONTINUE                                                  
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
                                                                        
                                                                        
                                                                        
C                                                                       
C     S. W. OTTO AND J. W. FLOWER, SUMMER 1985.                         
C                                                                       
C     CONSTRUCTS THE STRINGS FOR MEASURING PLANAR RxT WILSON LOOPS.     
C     IF R=1 THEN INSTEAD OF FETCHING A BARRED LEG FOR THE T DIRECTION W
C     CONSTRUCT IT EXPLICITLY USING THE REGULAR LINKS.                  
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE LPSTRM(PSTRT,PSTRU,PSTRB,R,T,DIR)                      
      INTEGER PSTRT(60),PSTRU(60),PSTRB(60)                             
      INTEGER R,T,DIR                                                   
C                                                                       
      INTEGER RR,TT,CNT,TEMP                                            
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
C                                                                       
C  FIRST GET THE TOP OF THE LOOP BARRING AS MANY LINKS AS POSSIBLE      
C                                                                       
      CNT = 1                                                           
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
      DO 1 TT=0,T-1                                                     
        PSTRT(CNT) = 26                                                 
    1   CNT = CNT+1                                                     
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
C                                                                       
      TEMP = 23+DIR                                                     
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
      DO 2 RR=0,R-1                                                     
        PSTRT(CNT) = TEMP                                               
    2   CNT = CNT+1                                                     
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
C                                                                       
      TEMP = 5+DIR                                                      
      PSTRT(CNT) = TEMP                                                 
      CNT = CNT+1                                                       
C                                                                       
      TEMP = 19+DIR                                                     
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
      DO 3 RR=1,R-2                                                     
        PSTRT(CNT) = TEMP                                               
    3   CNT = CNT+1                                                     
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
C                                                                       
      IF(R.GE.2) THEN                                                   
        TEMP = 5+DIR                                                    
        PSTRT(CNT) = TEMP                                               
        CNT = CNT+1                                                     
      ENDIF                                                             
CX10551 CALL DOINX( 551 )                                               
CY10551 CALL DOINY( 551 )                                               
C                                                                       
      DO 4 TT=0,T-1                                                     
        PSTRT(CNT) = 30                                                 
    4   CNT = CNT+1                                                     
CX20552 CALL DOOUTX( 552 )                                              
CY20552 CALL DOOUTY( 552 )                                              
      PSTRT(CNT) = 14                                                   
C                                                                       
C  NOW CONSTRUCT THE TIMELIKE LEG. IF R=1 MAKE EXPLICITLY               
C                                                                       
      CNT = 1                                                           
      TEMP = 23+DIR                                                     
CX10553 CALL DOINX( 553 )                                               
CY10553 CALL DOINY( 553 )                                               
      DO 5 RR=0,R-1                                                     
        PSTRU(CNT) = TEMP                                               
    5   CNT = CNT+1                                                     
CX20554 CALL DOOUTX( 554 )                                              
CY20554 CALL DOOUTY( 554 )                                              
C                                                                       
      IF(R.GT.1) THEN                                                   
        PSTRU(CNT) = 36                                                 
        CNT = CNT+1                                                     
      ELSE                                                              
CX10555 CALL DOINX( 555 )                                               
CY10555 CALL DOINY( 555 )                                               
        DO 6 TT=0,T-1                                                   
          PSTRU(CNT) = 4                                                
    6     CNT = CNT+1                                                   
CX20556 CALL DOOUTX( 556 )                                              
CY20556 CALL DOOUTY( 556 )                                              
CX10557 CALL DOINX( 557 )                                               
CY10557 CALL DOINY( 557 )                                               
        DO 7 TT=0,T-1                                                   
          PSTRU(CNT) = 30                                               
    7     CNT = CNT+1                                                   
CX20558 CALL DOOUTX( 558 )                                              
CY20558 CALL DOOUTY( 558 )                                              
      ENDIF                                                             
C                                                                       
      TEMP = 27+DIR                                                     
CX10559 CALL DOINX( 559 )                                               
CY10559 CALL DOINY( 559 )                                               
      DO 8 RR=0,R-1                                                     
        PSTRU(CNT) = TEMP                                               
    8   CNT = CNT+1                                                     
CX20560 CALL DOOUTX( 560 )                                              
CY20560 CALL DOOUTY( 560 )                                              
C                                                                       
      PSTRU(CNT) = 14                                                   
      CNT = CNT+1                                                       
C                                                                       
C  FINALLY MAKE THE BOTTOM OF THE LOOP                                  
C                                                                       
      CNT = 1                                                           
      TEMP = 1+DIR                                                      
      PSTRB(CNT) = TEMP                                                 
C                                                                       
      CNT = CNT+1                                                       
      TEMP = 15+DIR                                                     
CX10561 CALL DOINX( 561 )                                               
CY10561 CALL DOINY( 561 )                                               
      DO 9 RR=1,R-2                                                     
        PSTRB(CNT) = TEMP                                               
    9   CNT = CNT+1                                                     
CX20562 CALL DOOUTX( 562 )                                              
CY20562 CALL DOOUTY( 562 )                                              
C                                                                       
      IF(R.GE.2) THEN                                                   
        TEMP = 1+DIR                                                    
        PSTRB(CNT) = TEMP                                               
        CNT = CNT+1                                                     
      ENDIF                                                             
C                                                                       
      TEMP = 27+DIR                                                     
CX10563 CALL DOINX( 563 )                                               
CY10563 CALL DOINY( 563 )                                               
      DO 10 RR=0,R-1                                                    
        PSTRB(CNT) = TEMP                                               
   10   CNT = CNT+1                                                     
CX20564 CALL DOOUTX( 564 )                                              
CY20564 CALL DOOUTY( 564 )                                              
C                                                                       
      PSTRB(CNT) = 14                                                   
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
                                                                        
                                                                        
                                                                        
C                                                                       
C     STEVE OTTO (AND JON FLOWER) JULY 1984                             
C                                                                       
C     SETS UP INITIAL MOVEMENT TABLES.                                  
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, SUMMER 1987                                     
C                                                                       
      SUBROUTINE MAKE                                                   
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
CX10565 CALL DOINX( 565 )                                               
CY10565 CALL DOINY( 565 )                                               
C                                                                       
      DO 1 I=1,4                                                        
    1   ROT(I) = 18*(I-1)                                               
CX20566 CALL DOOUTX( 566 )                                              
CY20566 CALL DOOUTY( 566 )                                              
C                                                                       
      MOV(1) = 72                                                       
      MOV(2) = MOV(1)*LATT(1)                                           
      MOV(3) = MOV(2)*LATT(2)                                           
      MOV(4) = MOV(3)*LATT(3)                                           
CX10567 CALL DOINX( 567 )                                               
CY10567 CALL DOINY( 567 )                                               
C                                                                       
      DO 2 I=1,4                                                        
    2   LATT1(I) = LATT(I)-1                                            
CX20568 CALL DOOUTX( 568 )                                              
CY20568 CALL DOOUTY( 568 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C                                                                       
C     JON FLOWER, STEVE OTTO   MARCH 1985.  MOD: JUNE 85                
C                                                                       
C     MEASURES WILSON LOOPS UP TO RxT BARRING EVERYTHING IN SIGHT....   
C     ALSO MAKES LEGS.                                                  
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE MEASUR(DIR,RMIN,RMAX,NUMHIT)                           
      INTEGER DIR,RMIN,RMAX,NUMHIT                                      
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      PARAMETER ( MBAR=18*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL UBAR(MBAR),SBAR(MBAR)                                        
      COMMON /COMUBR/ UBAR,SBAR                                         
C                                                                       
      INTEGER R,NCHAR,COORD(4),TFLAG,SITE,T,PLEGS                       
      INTEGER STRT(60),STRU(60),STRB(60)                                
      REAL UP(18),DOWN(18),TOP(18),BOTTOM(18),STEMP(18,2)               
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
CX10569 CALL DOINX( 569 )                                               
CY10569 CALL DOINY( 569 )                                               
C                                                                       
      DO 1 R=0,12                                                       
        DO 1 T=0,12                                                     
          WLOOP(R+1,T+1) = 0.0                                          
    1     CONTINUE                                                      
CX20570 CALL DOOUTX( 570 )                                              
CY20570 CALL DOOUTY( 570 )                                              
C                                                                       
      CALL LINKBR(NUMHIT,DIR,SBAR)                                      
CX10571 CALL DOINX( 571 )                                               
CY10571 CALL DOINY( 571 )                                               
C                                                                       
C     ONE TIME SLICE AT A TIME                                          
C                                                                       
      DO 3 ICORD4=0,LATT(4)-1                                           
        COORD(4) = ICORD4                                               
        DO 4 T=1,12                                                     
          TFLAG = 0                                                     
          DO 5 R=RMIN,RMAX                                              
            IF(WHICH(R+1,T+1).EQ.1) TFLAG=1                             
    5       CONTINUE                                                    
C                                                                       
C    TFLAG=1 MEANS WE NEED LEGS OF THIS T                               
C                                                                       
          IF(TFLAG.EQ.1) CALL LEGMAK(T,ICORD4)                          
C                                                                       
          DO 6 R=RMIN,RMAX                                              
            IF((R.LE.T).AND.(WHICH(R+1,T+1).EQ.1)) THEN                 
              CALL LPSTRM(STRT,STRU,STRB,R,T,DIR)                       
              SITE = ICORD4*MOV(4)                                      
              PLEGS = 1                                                 
              DO 8 ICORD3=0,LATT(3)-1                                   
                COORD(3) = ICORD3                                       
                DO 8 ICORD2=0,LATT(2)-1                                 
                  COORD(2) = ICORD2                                     
                  DO 8 ICORD1=0,LATT(1)-1                               
                     COORD(1) = ICORD1                                  
C                                                                       
C                    FETCH AND DAGGER THE DOWN LEG                      
C                                                                       
                     CALL UDAG(LEGS(PLEGS),DOWN)                        
C                                                                       
C                    FIND THE TOP                                       
C                                                                       
                     CALL SYSLOP(COORD,SITE,STRT,TOP,NCHAR)             
C                                                                       
C                    FIND THE UP                                        
C                                                                       
                     CALL SYSLOP(COORD,SITE,STRU,UP,NCHAR)              
C                                                                       
C                    FIND THE BOTTOM                                    
C                                                                       
                     CALL SYSLOP(COORD,SITE,STRB,BOTTOM,NCHAR)          
C                                                                       
C                    PUT THEM ALL TOGETHER                              
C                                                                       
                     CALL MULT(BOTTOM,UP,STEMP(1,1))                    
                     CALL MULT(STEMP(1,1),TOP,STEMP(1,2))               
                     CALL MULT(STEMP(1,2),DOWN,STEMP(1,1))              
                     WLOOP(R+1,T+1) = WLOOP(R+1,T+1)+(STEMP(1,1)+       
     .                                  STEMP(9,1)+STEMP(17,1))/3.0     
                     PLEGS = PLEGS+18                                   
                     SITE = SITE+MOV(1)                                 
    8                CONTINUE                                           
            ENDIF                                                       
    6       CONTINUE                                                    
    4     CONTINUE                                                      
    3   CONTINUE                                                        
CX20572 CALL DOOUTX( 572 )                                              
CY20572 CALL DOOUTY( 572 )                                              
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE MULT(A,B,C)                                            
      REAL A(18),B(18),C(18)                                            
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
C                                                                       
      C(1)  = A(1)*B(1)-A(2)*B(2)+A(7)*B(3)-A(8)*B(4)+A(13)*B(5)        
     .       -A(14)*B(6)                                                
      C(2)  = A(1)*B(2)+A(2)*B(1)+A(7)*B(4)+A(8)*B(3)+A(13)*B(6)        
     .       +A(14)*B(5)                                                
      C(3)  = A(3)*B(1)-A(4)*B(2)+A(9)*B(3)-A(10)*B(4)+A(15)*B(5)       
     .       -A(16)*B(6)                                                
      C(4)  = A(3)*B(2)+A(4)*B(1)+A(9)*B(4)+A(10)*B(3)+A(15)*B(6)       
     .       +A(16)*B(5)                                                
      C(5)  = A(5)*B(1)-A(6)*B(2)+A(11)*B(3)-A(12)*B(4)+A(17)*B(5)      
     .       -A(18)*B(6)                                                
      C(6)  = A(5)*B(2)+A(6)*B(1)+A(11)*B(4)+A(12)*B(3)+A(17)*B(6)      
     .       +A(18)*B(5)                                                
      C(7)  = A(1)*B(7)-A(2)*B(8)+A(7)*B(9)-A(8)*B(10)+A(13)*B(11)      
     .       -A(14)*B(12)                                               
      C(8)  = A(1)*B(8)+A(2)*B(7)+A(7)*B(10)+A(8)*B(9)+A(13)*B(12)      
     .       +A(14)*B(11)                                               
      C(9)  = A(3)*B(7)-A(4)*B(8)+A(9)*B(9)-A(10)*B(10)+A(15)*B(11)     
     .       -A(16)*B(12)                                               
      C(10) = A(3)*B(8)+A(4)*B(7)+A(9)*B(10)+A(10)*B(9)+A(15)*B(12)     
     .       +A(16)*B(11)                                               
      C(11) = A(5)*B(7)-A(6)*B(8)+A(11)*B(9)-A(12)*B(10)+A(17)*B(11)    
     .       -A(18)*B(12)                                               
      C(12) = A(5)*B(8)+A(6)*B(7)+A(11)*B(10)+A(12)*B(9)+A(17)*B(12)    
     .       +A(18)*B(11)                                               
      C(13) = A(1)*B(13)-A(2)*B(14)+A(7)*B(15)-A(8)*B(16)+A(13)*B(17)   
     .       -A(14)*B(18)                                               
      C(14) = A(1)*B(14)+A(2)*B(13)+A(7)*B(16)+A(8)*B(15)+A(13)*B(18)   
     .       +A(14)*B(17)                                               
      C(15) = A(3)*B(13)-A(4)*B(14)+A(9)*B(15)-A(10)*B(16)+A(15)*B(17)  
     .       -A(16)*B(18)                                               
      C(16) = A(3)*B(14)+A(4)*B(13)+A(9)*B(16)+A(10)*B(15)+A(15)*B(18)  
     .       +A(16)*B(17)                                               
      C(17) = A(5)*B(13)-A(6)*B(14)+A(11)*B(15)-A(12)*B(16)+A(17)*B(17) 
     .       -A(18)*B(18)                                               
      C(18) = A(5)*B(14)+A(6)*B(13)+A(11)*B(16)+A(12)*B(15)+A(17)*B(18) 
     .       +A(18)*B(17)                                               
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PRNSE2(ISEED)                                          
      PARAMETER (MULT = 1103515245, IADD=12345)                         
      INTEGER AAA,BBB,RANDX                                             
      COMMON /COMPRA/ AAA,BBB,RANDX                                     
      SAVE /COMPRA/                                                     
CX30033 CALL SBINX (33)                                                 
CY30033 CALL SBINY (33)                                                 
      AAA = MULT                                                        
      BBB = IADD                                                        
      RANDX = ISEED                                                     
CX40034 CALL SBOUTX (34)                                                
CY40034 CALL SBOUTY (34)                                                
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C    J. W. FLOWER   AUG 1985                                            
C                                                                       
C     CONSTRUCTS THE THREE QUARK OBSERVABLE :                           
C                                                                       
C                            (1)  (2)  (3)                              
C       EPSILON   EPSILON    U    U    U                                
C              IJK       PQR  IP   JQ   KR                              
C                                                                       
C    THIS CODE IS CERTAINLY NON-OPTIMAL BUT SACRIFICES                  
C            SPEED FOR CODING SIMPLICITY                                
C                                                                       
C    CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,          
C    AND DAVID WALKER, 1987                                             
C                                                                       
      SUBROUTINE OBSERV(TOT,U1,U2,U3)                                   
      REAL TOT(2),U1(2,9),U2(2,9),U3(2,9)                               
C                                                                       
      INTEGER P,Q,R,EPSILO(3,3,3)                                       
      REAL FAC                                                          
CX30035 CALL SBINX (35)                                                 
CY30035 CALL SBINY (35)                                                 
CX10573 CALL DOINX( 573 )                                               
CY10573 CALL DOINY( 573 )                                               
C                                                                       
      DO 1 I=1,3                                                        
        DO 1 J=1,3                                                      
          DO 1 K=1,3                                                    
    1       EPSILO(I,J,K) = 0                                           
CX20574 CALL DOOUTX( 574 )                                              
CY20574 CALL DOOUTY( 574 )                                              
C                                                                       
      EPSILO(1,2,3) = 1                                                 
      EPSILO(2,3,1) = 1                                                 
      EPSILO(3,1,2) = 1                                                 
      EPSILO(2,1,3) = -1                                                
      EPSILO(1,3,2) = -1                                                
      EPSILO(3,2,1) = -1                                                
CX10575 CALL DOINX( 575 )                                               
CY10575 CALL DOINY( 575 )                                               
C                                                                       
      DO 2 I=0,2                                                        
        DO 2 P=0,2                                                      
          DO 2 J=0,2                                                    
            DO 2 Q=0,2                                                  
              DO 2 K=0,2                                                
                IF(EPSILO(I+1,J+1,K+1).NE.0) THEN                       
                DO 3 R=0,2                                              
                  IF(EPSILO(P+1,Q+1,R+1).NE.0) THEN                     
                    FAC = EPSILO(I+1,J+1,K+1)*EPSILO(P+1,Q+1,R+1)       
                    TOT(1) = TOT(1)+ FAC*U1(1,3*I+P+1)*U2(1,3*J+Q+1)    
     .                       *U3(1,3*K+R+1)                             
                    TOT(1) = TOT(1)- FAC*U1(2,3*I+P+1)*U2(2,3*J+Q+1)    
     .                       *U3(1,3*K+R+1)                             
                    TOT(1) = TOT(1)- FAC*U1(1,3*I+P+1)*U2(2,3*J+Q+1)    
     .                       *U3(2,3*K+R+1)                             
                    TOT(1) = TOT(1)- FAC*U1(2,3*I+P+1)*U2(1,3*J+Q+1)    
     .                       *U3(2,3*K+R+1)                             
                    TOT(2) = TOT(2)+ FAC*U1(1,3*I+P+1)*U2(1,3*J+Q+1)    
     .                       *U3(2,3*K+R+1)                             
                    TOT(2) = TOT(2)+ FAC*U1(1,3*I+P+1)*U2(2,3*J+Q+1)    
     .                       *U3(1,3*K+R+1)                             
                    TOT(2) = TOT(2)+ FAC*U1(2,3*I+P+1)*U2(1,3*J+Q+1)    
     .                       *U3(1,3*K+R+1)                             
                    TOT(2) = TOT(2)- FAC*U1(2,3*I+P+1)*U2(2,3*J+Q+1)    
     .                       *U3(2,3*K+R+1)                             
                  ENDIF                                                 
    3             CONTINUE                                              
                ENDIF                                                   
    2           CONTINUE                                                
CX20576 CALL DOOUTX( 576 )                                              
CY20576 CALL DOOUTY( 576 )                                              
CX40036 CALL SBOUTX (36)                                                
CY40036 CALL SBOUTY (36)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
        SUBROUTINE PRNSET(ISEED)                                        
        COMMON /RAN/IRAND0,IRAND1,IAB0,IAB1,ICB0,ICB1                   
        SAVE /RAN/                                                      
C                                                                       
        PARAMETER (MULT0=20077,MULT1=16838,IADD0=12345,IADD1=0)         
        PARAMETER (IPOW16=2**16)                                        
CX30037 CALL SBINX (37)                                                 
CY30037 CALL SBINY (37)                                                 
C                                                                       
        ISEED0=MOD(ISEED,IPOW16)                                        
        ISEED1=ISEED/IPOW16                                             
        IAB0=MULT0                                                      
        IAB1=MULT1                                                      
        ICB0=IADD0                                                      
        ICB1=IADD1                                                      
        CALL LMULT(J0,J1,ISEED0,ISEED1,IAB0,IAB1)                       
        CALL LADD(IRAND0,IRAND1,J0,J1,ICB0,ICB1)                        
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
        RETURN                                                          
        END                                                             
C                                                                       
CCCCCCCCCCCCCCCCCCCCC END OF ROUTINE PRNSET CCCCCCCCCCCCCCCCCCCCCCC C   
C                                                                       
        FUNCTION PRANF(IDUMMY)                                          
        COMMON /RAN/IRAND0,IRAND1,IAB0,IAB1,ICB0,ICB1                   
        SAVE /RAN/                                                      
C                                                                       
        DOUBLE PRECISION RANDB,DIVFAC,POW16                             
        PARAMETER (DIVFAC=2.147483648D9,POW16=6.5536D4)                 
CX30039 CALL SBINX (39)                                                 
CY30039 CALL SBINY (39)                                                 
C                                                                       
        RANDB=(DBLE(IRAND0)+POW16*DBLE(IRAND1))/DIVFAC                  
        PRANF=REAL(RANDB)                                               
        CALL LMULT(J0,J1,IRAND0,IRAND1,IAB0,IAB1)                       
        CALL LADD(IRAND0,IRAND1,J0,J1,ICB0,ICB1)                        
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C                                                                       
CCCCCCCCCCCCCCCCCCCCC END OF ROUTINE PRANF CCCCCCCCCCCCCCCCCCCCCCC C    
C                                                                       
        SUBROUTINE LADD(I0,I1,J0,J1,K0,K1)                              
C                                                                       
        PARAMETER (IPOW16=2**16,IPOW15=2**15)                           
CX30041 CALL SBINX (41)                                                 
CY30041 CALL SBINY (41)                                                 
C                                                                       
        I0=MOD(K0+J0,IPOW16)                                            
        I1=MOD((K0+J0)/IPOW16+K1+J1,IPOW15)                             
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C                                                                       
CCCCCCCCCCCCCCCCCCCCC END OF ROUTINE LADD CCCCCCCCCCCCCCCCCCCCCCC C     
C                                                                       
        SUBROUTINE LMULT(I0,I1,J0,J1,K0,K1)                             
C                                                                       
        PARAMETER (IPOW30=2**30,IPOW16=2**16,IPOW15=2**15)              
CX30043 CALL SBINX (43)                                                 
CY30043 CALL SBINY (43)                                                 
C                                                                       
        K = K0*J0                                                       
        IF(K.LT.0) K = (K+IPOW30)+IPOW30                                
        I0=MOD(K,IPOW16)                                                
        L = K0*J1+K1*J0                                                 
        IF(L.LT.0) L = (L+IPOW30)+IPOW30                                
        K = K/IPOW16+L                                                  
        IF(K.LT.0) K = (K+IPOW30)+IPOW30                                
        I1=MOD(K,IPOW15)                                                
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C                                                                       
CCCCCCCCCCCCCCCCCCCCCCCCCCC END OF ROUTINE LMULT CCCCCCCCCCCCCCCCCCCCCCC
C                                                                       
        FUNCTION IGETSD(IDUMMY)                                         
C                                                                       
        COMMON /RAN/IRAND0,IRAND1,IAB0,IAB1,ICB0,ICB1                   
        SAVE /RAN/                                                      
CX30045 CALL SBINX (45)                                                 
CY30045 CALL SBINY (45)                                                 
C                                                                       
        IGETSD = IRAND0 + IRAND1*2**16                                  
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C                                                                       
CCCCCCCCCCCCCCCCCCCCCCCCCCC  END OF ROUTINE IGETSD  CCCCCCCCCCCCCCCCCCCC
C                                                                       
                                                                        
C                                                                       
C     STEVE OTTO  AND JON FLOWER, JULY 1984.                            
C                                                                       
C     DUE TO THE STRANGE SHAPE OF THE SU(3) GROUP MANIFOLD (WHATEVER IT 
C     THE LINK MATRICES TEND TO DRIFT OFF AND MUST BE PROJECTED BACK    
C     PERIODICALLY. THE SITUATION IS OBVIOUSLY NOT HELPED BY THE LIMITED
C     ACCURACY OF THE FIXED POINT FORMAT.                             */
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE PROJEC                                                 
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      REAL X(12),Z(6),FAC                                               
      INTEGER SITE                                                      
      INTEGER PU                                                        
CX30047 CALL SBINX (47)                                                 
CY30047 CALL SBINY (47)                                                 
C                                                                       
      NUMLIN = 4*LATT(1)*LATT(2)*LATT(3)*LATT(4)                        
CX10577 CALL DOINX( 577 )                                               
CY10577 CALL DOINY( 577 )                                               
C                                                                       
      DO 1 I=0,NUMLIN-1                                                 
        SITE = ROT(2)*I                                                 
        PU = SITE                                                       
        DO 10 J=1,12                                                    
   10     X(J) = U1(PU+J)                                               
        Z(1) = X(3)*X(11) - X(4)*X(12) - X(5)*X(9) + X(6)*X(10)         
        Z(2) = -X(4)*X(11) - X(3)*X(12) + X(6)*X(9) + X(5)*X(10)        
        Z(3) = X(5)*X(7) - X(6)*X(8) - X(1)*X(11) + X(2)*X(12)          
        Z(4) = -X(6)*X(7) - X(5)*X(8) + X(2)*X(11) + X(1)*X(12)         
        Z(5) = X(1)*X(9) - X(2)*X(10) - X(3)*X(7) + X(4)*X(8)           
        Z(6) = -X(2)*X(9) - X(1)*X(10) + X(4)*X(7) + X(3)*X(8)          
        FAC  = 0.0                                                      
        DO 2 J=1,6                                                      
    2     FAC = FAC+X(J)*X(J)                                           
        FAC = 1.0/SQRT(FAC)                                             
        DO 3 J=1,6                                                      
    3     X(J) = X(J)*FAC                                               
        X(7)  = Z(3)*X(5) - Z(4)*X(6) - Z(5)*X(3) + Z(6)*X(4)           
        X(8)  = -Z(3)*X(6) - Z(4)*X(5) + Z(5)*X(4) + Z(6)*X(3)          
        X(9)  = Z(5)*X(1) - Z(6)*X(2) - Z(1)*X(5) + Z(2)*X(6)           
        X(10) = -Z(5)*X(2) - Z(6)*X(1) + Z(1)*X(6) + Z(2)*X(5)          
        X(11) = Z(1)*X(3) - Z(2)*X(4) - Z(3)*X(1) + Z(4)*X(2)           
        X(12) = -Z(1)*X(4) - Z(2)*X(3) + Z(3)*X(2) + Z(4)*X(1)          
        FAC=0.0                                                         
        DO 4 J=1,6                                                      
    4     FAC = FAC+X(6+J)*X(6+J)                                       
        FAC = 1.0/SQRT(FAC)                                             
        DO 5 J=1,6                                                      
    5     X(6+J) = X(6+J)*FAC                                           
        CALL CFILL(U1(PU+1),X)                                          
    1 CONTINUE                                                          
CX20578 CALL DOOUTX( 578 )                                              
CY20578 CALL DOOUTY( 578 )                                              
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C                                                                       
C     J. W. FLOWER & STEVE OTTO, SEPTEMBER 1984.                        
C                                                                       
C     USED TO ACTUALLY CALCULATE THE THREE QUARK OBSERVABLES.       */  
C                                                                       
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE QQQLPC(R1,R2,R3,T,LPSTR,TOT)                           
      INTEGER R1,R2,R3,T                                                
      INTEGER LPSTR(512)                                                
      REAL TOT(2)                                                       
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      INTEGER COORD(4),ODIR,SITE,CNT                                    
      INTEGER PSTR(512)                                                 
      REAL U1(18),U2(18),U3(18),BITS(18,3),TEMP(18)                     
CX30049 CALL SBINX (49)                                                 
CY30049 CALL SBINY (49)                                                 
C                                                                       
      SITE = T*MOV(4)                                                   
C                                                                       
      COORD(4) = T                                                      
CX10579 CALL DOINX( 579 )                                               
CY10579 CALL DOINY( 579 )                                               
      DO 1 ICORD3=0,LATT(3)-1                                           
        COORD(3) = ICORD3                                               
        DO 1 ICORD2=0,LATT(2)-1                                         
          COORD(2) = ICORD2                                             
          DO 1 ICORD1=0,LATT(1)-1                                       
            COORD(1) = ICORD1                                           
C                                                                       
C     FIRST GET THE MATRIX CORRESPONDING TO DIRECTION R1.               
C                                                                       
C            CALL CPYMAT(PSTR,LPSTR,512) LPOINTER 4-3-89                
            CALL ICPYMT(PSTR,LPSTR,512)                                 
            CNT = 1                                                     
            DO 2 I=1,3                                                  
              CALL SYSLOP(COORD,SITE,PSTR(CNT),BITS(1,I),NCHAR)         
    2         CNT = NCHAR+CNT                                           
            IF(R1.GT.0) THEN                                            
              CALL MULT(BITS(1,1),BITS(1,2),TEMP)                       
              CALL MULT(TEMP,BITS(1,3),U1)                              
            ELSE                                                        
              DO 3 I=1,18                                               
    3           U1(I) = BITS(I,2)                                       
             ENDIF                                                      
C                                                                       
C      NEXT GET THE MATRIX CORRESPONDING TO DIRECTION R2.               
C                                                                       
            DO 4 I=1,3                                                  
              CALL SYSLOP(COORD,SITE,PSTR(CNT),BITS(1,I),NCHAR)         
    4         CNT = CNT+NCHAR                                           
            IF(R2.GT.0) THEN                                            
              CALL MULT(BITS(1,1),BITS(1,2),TEMP)                       
              CALL MULT(TEMP,BITS(1,3),U2)                              
            ELSE                                                        
              DO 5 I=1,18                                               
    5           U2(I) = BITS(I,2)                                       
            ENDIF                                                       
C                                                                       
C       NOW LOOP OVER THE TWO (CURRENTLY) DIRECTIONS FOR THE R3 LEG     
C       CALCULATING THE OBSERVABLES AS WE GO.                           
C                                                                       
            DO 6 ODIR=1,2                                               
              DO 7 I=1,3                                                
                CALL SYSLOP(COORD,SITE,PSTR(CNT),BITS(1,I),NCHAR)       
    7           CNT = CNT+NCHAR                                         
              IF(R3.GT.0) THEN                                          
                CALL MULT(BITS(1,1),BITS(1,2),TEMP)                     
                CALL MULT(TEMP,BITS(1,3),U3)                            
              ELSE                                                      
                DO 8 I=1,18                                             
    8             U3(I) = BITS(I,2)                                     
              ENDIF                                                     
C                                                                       
C        CALCULATE OBSERVABLES                                          
C                                                                       
              CALL OBSERV(TOT,U1,U2,U3)                                 
C                                                                       
    6         CONTINUE                                                  
C                                                                       
C        MOVE TO NEXT SITE IN LATTICE.                                  
C                                                                       
            SITE = SITE+MOV(1)                                          
    1     CONTINUE                                                      
CX20580 CALL DOOUTX( 580 )                                              
CY20580 CALL DOOUTY( 580 )                                              
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C                                                                       
C     J .W .FLOWER  AUG 1985                                            
C                                                                       
C                                                                       
C     FROM ABOVE THE THREE QUARK ORIENTATION IS SPECIFIED AS FOLLOWS :  
C                                                                       
C                                                                       
C                                   * ^                                 
C                                   |                                   
C                                   | R3                                
C                                   |                                   
C                     * ------------|-------------------- *             
C                     <---- R1 ---->|<-------- R2 -------->             
C                                                                       
C     THE LIMB WITH R1 & R2 IS TAKEN TO LIE IN DIRECTION "DIR" AND THE  
C     OTHER LIMB IS TAKEN IN THE TWO TRANSVERSE DIRECTIONS SUCCESSIVELY.
C     NOTE : SINCE WE HAVE ACCESS TO ONLY BARRED 'LEGS' WE NEED TO MAKE 
C     SOME TIMELIKE PIECES ON THE FLY.                                  
C                                                                     */
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,         
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE QQQLPS(R1,R2,R3,T,DIR,STR)                             
      INTEGER R1,R2,R3,T,DIR                                            
      INTEGER STR(*)                                                    
      INTEGER CTR,ODIR,TEMP                                             
CX30051 CALL SBINX (51)                                                 
CY30051 CALL SBINY (51)                                                 
C                                                                       
C  DO R1 LIMB - GO BACK IN "DIR" FIRST.                                 
C                                                                       
      CTR=1                                                             
      TEMP=5+DIR                                                        
CX10581 CALL DOINX( 581 )                                               
CY10581 CALL DOINY( 581 )                                               
      DO 1 I=0,R1-1                                                     
        STR(CTR) = TEMP                                                 
    1   CTR=CTR+1                                                       
CX20582 CALL DOOUTX( 582 )                                              
CY20582 CALL DOOUTY( 582 )                                              
C                                                                       
      TEMP=23+DIR                                                       
CX10583 CALL DOINX( 583 )                                               
CY10583 CALL DOINY( 583 )                                               
      DO 2 I=0,R1-1                                                     
        STR(CTR) = TEMP                                                 
    2   CTR=CTR+1                                                       
CX20584 CALL DOOUTX( 584 )                                              
CY20584 CALL DOOUTY( 584 )                                              
C                                                                       
      STR(CTR) = 14                                                     
      CTR=CTR+1                                                         
C                                                                       
      TEMP=27+DIR                                                       
CX10585 CALL DOINX( 585 )                                               
CY10585 CALL DOINY( 585 )                                               
      DO 3 I=0,R1-1                                                     
        STR(CTR) = TEMP                                                 
    3   CTR=CTR+1                                                       
CX20586 CALL DOOUTX( 586 )                                              
CY20586 CALL DOOUTY( 586 )                                              
C                                                                       
      IF(R1.NE.0) THEN                                                  
        STR(CTR) = 36                                                   
        CTR=CTR+1                                                       
      ELSE                                                              
CX10587 CALL DOINX( 587 )                                               
CY10587 CALL DOINY( 587 )                                               
        DO 4 I=0,T-1                                                    
          STR(CTR) = 4                                                  
    4     CTR=CTR+1                                                     
CX20588 CALL DOOUTX( 588 )                                              
CY20588 CALL DOOUTY( 588 )                                              
CX10589 CALL DOINX( 589 )                                               
CY10589 CALL DOINY( 589 )                                               
        DO 5 I=0,T-1                                                    
          STR(CTR) = 30                                                 
    5     CTR=CTR+1                                                     
CX20590 CALL DOOUTX( 590 )                                              
CY20590 CALL DOOUTY( 590 )                                              
      ENDIF                                                             
C                                                                       
      TEMP=23+DIR                                                       
CX10591 CALL DOINX( 591 )                                               
CY10591 CALL DOINY( 591 )                                               
      DO 6 I=0,R1-1                                                     
        STR(CTR) = TEMP                                                 
    6   CTR=CTR+1                                                       
CX20592 CALL DOOUTX( 592 )                                              
CY20592 CALL DOOUTY( 592 )                                              
C                                                                       
      STR(CTR) = 14                                                     
      CTR=CTR+1                                                         
C                                                                       
      TEMP=27+DIR                                                       
CX10593 CALL DOINX( 593 )                                               
CY10593 CALL DOINY( 593 )                                               
      DO 7 I=0,R1-1                                                     
        STR(CTR) = TEMP                                                 
    7   CTR=CTR+1                                                       
CX20594 CALL DOOUTX( 594 )                                              
CY20594 CALL DOOUTY( 594 )                                              
CX10595 CALL DOINX( 595 )                                               
CY10595 CALL DOINY( 595 )                                               
C                                                                       
      DO 8 I=0,T-1                                                      
        STR(CTR) = 26                                                   
    8   CTR=CTR+1                                                       
CX20596 CALL DOOUTX( 596 )                                              
CY20596 CALL DOOUTY( 596 )                                              
C                                                                       
      TEMP=1+DIR                                                        
CX10597 CALL DOINX( 597 )                                               
CY10597 CALL DOINY( 597 )                                               
      DO 9 I=0,R1-1                                                     
        STR(CTR) = TEMP                                                 
    9   CTR=CTR+1                                                       
CX20598 CALL DOOUTX( 598 )                                              
CY20598 CALL DOOUTY( 598 )                                              
CX10599 CALL DOINX( 599 )                                               
CY10599 CALL DOINY( 599 )                                               
C                                                                       
      DO 10 I=0,T-1                                                     
        STR(CTR) = 30                                                   
   10   CTR=CTR+1                                                       
CX20600 CALL DOOUTX( 600 )                                              
CY20600 CALL DOOUTY( 600 )                                              
C                                                                       
      STR(CTR) = 14                                                     
      CTR=CTR+1                                                         
C                                                                       
C   DO R2 LIMB - GO FORWARD IN "DIR" FIRST.                             
C                                                                       
      TEMP=1+DIR                                                        
CX10601 CALL DOINX( 601 )                                               
CY10601 CALL DOINY( 601 )                                               
      DO 11 I=0,R2-1                                                    
        STR(CTR) = TEMP                                                 
   11   CTR=CTR+1                                                       
CX20602 CALL DOOUTX( 602 )                                              
CY20602 CALL DOOUTY( 602 )                                              
C                                                                       
      TEMP=27+DIR                                                       
CX10603 CALL DOINX( 603 )                                               
CY10603 CALL DOINY( 603 )                                               
      DO 12 I=0,R2-1                                                    
        STR(CTR) = TEMP                                                 
   12   CTR=CTR+1                                                       
CX20604 CALL DOOUTX( 604 )                                              
CY20604 CALL DOOUTY( 604 )                                              
C                                                                       
      STR(CTR) = 14                                                     
      CTR=CTR+1                                                         
C                                                                       
      TEMP=23+DIR                                                       
CX10605 CALL DOINX( 605 )                                               
CY10605 CALL DOINY( 605 )                                               
      DO 13 I=0,R2-1                                                    
        STR(CTR) = TEMP                                                 
   13   CTR=CTR+1                                                       
CX20606 CALL DOOUTX( 606 )                                              
CY20606 CALL DOOUTY( 606 )                                              
C                                                                       
      IF(R2.NE.0) THEN                                                  
        STR(CTR) = 36                                                   
        CTR=CTR+1                                                       
      ELSE                                                              
CX10607 CALL DOINX( 607 )                                               
CY10607 CALL DOINY( 607 )                                               
        DO 14 I=0,T-1                                                   
          STR(CTR) = 4                                                  
   14     CTR=CTR+1                                                     
CX20608 CALL DOOUTX( 608 )                                              
CY20608 CALL DOOUTY( 608 )                                              
CX10609 CALL DOINX( 609 )                                               
CY10609 CALL DOINY( 609 )                                               
        DO 15 I=0,T-1                                                   
          STR(CTR) = 30                                                 
   15     CTR=CTR+1                                                     
CX20610 CALL DOOUTX( 610 )                                              
CY20610 CALL DOOUTY( 610 )                                              
      ENDIF                                                             
C                                                                       
      TEMP=27+DIR                                                       
CX10611 CALL DOINX( 611 )                                               
CY10611 CALL DOINY( 611 )                                               
      DO 16 I=0,R2-1                                                    
        STR(CTR) = TEMP                                                 
   16   CTR=CTR+1                                                       
CX20612 CALL DOOUTX( 612 )                                              
CY20612 CALL DOOUTY( 612 )                                              
C                                                                       
      STR(CTR) = 14                                                     
      CTR=CTR+1                                                         
C                                                                       
      TEMP=23+DIR                                                       
CX10613 CALL DOINX( 613 )                                               
CY10613 CALL DOINY( 613 )                                               
      DO 17 I=0,R2-1                                                    
        STR(CTR) = TEMP                                                 
   17   CTR=CTR+1                                                       
CX20614 CALL DOOUTX( 614 )                                              
CY20614 CALL DOOUTY( 614 )                                              
CX10615 CALL DOINX( 615 )                                               
CY10615 CALL DOINY( 615 )                                               
C                                                                       
      DO 18 I=0,T-1                                                     
        STR(CTR) = 26                                                   
   18   CTR=CTR+1                                                       
CX20616 CALL DOOUTX( 616 )                                              
CY20616 CALL DOOUTY( 616 )                                              
C                                                                       
      TEMP=5+DIR                                                        
CX10617 CALL DOINX( 617 )                                               
CY10617 CALL DOINY( 617 )                                               
      DO 19 I=0,R2-1                                                    
        STR(CTR) = TEMP                                                 
   19   CTR=CTR+1                                                       
CX20618 CALL DOOUTX( 618 )                                              
CY20618 CALL DOOUTY( 618 )                                              
CX10619 CALL DOINX( 619 )                                               
CY10619 CALL DOINY( 619 )                                               
C                                                                       
      DO 20 I=0,T-1                                                     
        STR(CTR) = 30                                                   
   20   CTR=CTR+1                                                       
CX20620 CALL DOOUTX( 620 )                                              
CY20620 CALL DOOUTY( 620 )                                              
C                                                                       
      STR(CTR) = 14                                                     
      CTR=CTR+1                                                         
CX10621 CALL DOINX( 621 )                                               
CY10621 CALL DOINY( 621 )                                               
C                                                                       
C  NOW GET OTHER LEG ... LOOP OVER TRANSVERSE DIRECTIONS                
C   .... AND ALSO GO FORWARD AND (EVENTUALLY) BACK TO R3.               
C                                                                       
      DO 21 NTRANS=1,2                                                  
        ODIR = MOD((DIR + NTRANS),3)                                    
        TEMP=1+ODIR                                                     
        DO 22 I=0,R3-1                                                  
          STR(CTR) = TEMP                                               
   22     CTR=CTR+1                                                     
        TEMP=27+ODIR                                                    
        DO 23 I=0,R3-1                                                  
          STR(CTR) = TEMP                                               
   23     CTR=CTR+1                                                     
        STR(CTR) = 14                                                   
        CTR=CTR+1                                                       
C                                                                       
        TEMP=23+ODIR                                                    
        DO 24 I=0,R3-1                                                  
          STR(CTR) = TEMP                                               
   24     CTR=CTR+1                                                     
        IF(R3.NE.0) THEN                                                
          STR(CTR) = 36                                                 
          CTR=CTR+1                                                     
        ELSE                                                            
          DO 25 I=0,T-1                                                 
            STR(CTR) = 4                                                
   25       CTR=CTR+1                                                   
          DO 26 I=0,T-1                                                 
            STR(CTR) = 30                                               
   26       CTR=CTR+1                                                   
        ENDIF                                                           
        TEMP=27+ODIR                                                    
        DO 27 I=0,R3-1                                                  
          STR(CTR) = TEMP                                               
   27     CTR=CTR+1                                                     
        STR(CTR) = 14                                                   
        CTR=CTR+1                                                       
C                                                                       
        TEMP=23+ODIR                                                    
        DO 28 I=0,R3-1                                                  
          STR(CTR) = TEMP                                               
   28     CTR=CTR+1                                                     
        DO 29 I=0,T-1                                                   
          STR(CTR) = 26                                                 
   29     CTR=CTR+1                                                     
        TEMP=5+ODIR                                                     
        DO 30 I=0,R3-1                                                  
          STR(CTR) = TEMP                                               
   30     CTR=CTR+1                                                     
        DO 31 I=0,T-1                                                   
          STR(CTR) = 30                                                 
   31     CTR=CTR+1                                                     
        STR(CTR) = 14                                                   
C                                                                       
        CTR=CTR+1                                                       
   21   CONTINUE                                                        
CX20622 CALL DOOUTX( 622 )                                              
CY20622 CALL DOOUTY( 622 )                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
                                                                        
                                                                        
                                                                        
C                                                                       
C   J. W. FLOWER & STEVE OTTO  SEPTEMBER 1984                           
C                                                                       
C   MEASURES THREE QUARK OBSERVABLES USING UBARS.                       
C                                                                       
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,           
C   AND DAVID WALKER, 1987                                              
C                                                                       
      SUBROUTINE QQQMEA(R1,R2,R3,TT)                                    
      INTEGER R1,R2,R3,TT                                               
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      INTEGER NUMSTG                                                    
      COMMON /COMFIL/ NUMSTG                                            
C                                                                       
      INTEGER DIR,T                                                     
CMWB 6.22.88      CHARACTER*512 LPSTR                                   
      INTEGER LPSTR(512)                                                
      REAL TOT(2),VALS(2,15)                                            
CX30053 CALL SBINX (53)                                                 
CY30053 CALL SBINY (53)                                                 
C                                                                       
      VALS(1,TT) = 0.0                                                  
      VALS(2,TT) = 0.0                                                  
CX10623 CALL DOINX( 623 )                                               
CY10623 CALL DOINY( 623 )                                               
C                                                                       
      DO 1 T=0,LATT(4)-1                                                
        CALL LEGMAK(TT,T)                                               
        TOT(1) = 0.0                                                    
        TOT(2) = 0.0                                                    
        DO 2 DIR=0,2                                                    
          CALL QQQLPS(R1,R2,R3,TT,DIR,LPSTR)                            
          CALL QQQLPC(R1,R2,R3,T,LPSTR,TOT)                             
    2     CONTINUE                                                      
        VALS(1,TT) = VALS(1,TT) + TOT(1)                                
        VALS(2,TT) = VALS(2,TT) + TOT(2)                                
    1   CONTINUE                                                        
CX20624 CALL DOOUTX( 624 )                                              
CY20624 CALL DOOUTY( 624 )                                              
C                                                                       
      VALS(1,TT) = VALS(1,TT) / REAL(NUMSTG*3*2)                        
      VALS(2,TT) = VALS(2,TT) / REAL(NUMSTG*3*2)                        
      WRITE(6,*) 'Q (',R1,R2,R3,')','  X ',TT,' : ',                    
     *           VALS(1,TT),VALS(2,TT)                                  
CX40054 CALL SBOUTX (54)                                                
CY40054 CALL SBOUTY (54)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C   MAKE THE STRINGS USED TO MEASURE OFF-AXIS WILSON LOOPS.             
C   STRING COMES IN THREE PARTS.                                        
C    (i)  OUTWARD MOVES AND FAKE BACK                                   
C    (ii) FETCH LEG AND FAKE BACK.                                      
C    (iii GET TOP PART AND FAKE BACK.                                   
C   BARS ALL T LINKS EXCEPT IF R1+R2=1.                                 
C                                                                       
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,           
C   AND DAVID WALKER, 1987                                              
C                                                                       
      SUBROUTINE ROTLPS(PSTR,R1,R2,T,DIR1,DIR2)                         
      INTEGER PSTR(256)                                                 
      INTEGER R1,R2,T,DIR1,DIR2                                         
C                                                                       
      INTEGER R,DIR,CNT,TEMP                                            
CX30055 CALL SBINX (55)                                                 
CY30055 CALL SBINY (55)                                                 
C                                                                       
C   GO OUT IN SPATIAL DIRECTION.                                        
C                                                                       
      NUM1=0                                                            
      NUM2=0                                                            
      CNT = 1                                                           
    1 CONTINUE                                                          
      IF((NUM1.LT.R1) .AND. (NUM2.LT.R2)) THEN                          
        TEMP=1+DIR2                                                     
        PSTR(CNT) = TEMP                                                
        CNT=CNT+1                                                       
        NUM2=NUM2+1                                                     
        IF((NUM1+NUM2).LT.(R1+R2-1)) THEN                               
          TEMP=15+DIR1                                                  
        ELSE                                                            
          TEMP=1+DIR1                                                   
        ENDIF                                                           
        PSTR(CNT) = TEMP                                                
        CNT=CNT+1                                                       
        NUM1=NUM1+1                                                     
        GOTO 1                                                          
      ENDIF                                                             
C                                                                       
      IF(NUM1.LT.R1) THEN                                               
        DIR=DIR1                                                        
      ELSE                                                              
        DIR=DIR2                                                        
      ENDIF                                                             
CX10625 CALL DOINX( 625 )                                               
CY10625 CALL DOINY( 625 )                                               
C                                                                       
      DO 2 R=(NUM1+NUM2),(R1+R2-1)                                      
        IF((DIR.EQ.DIR1).AND.(R.LT.R1+R2-1).AND.(R.GT.0)) THEN          
          TEMP=15+DIR                                                   
        ELSE                                                            
          TEMP=1+DIR                                                    
        ENDIF                                                           
        PSTR(CNT) = TEMP                                                
        CNT=CNT+1                                                       
    2   CONTINUE                                                        
CX20626 CALL DOOUTX( 626 )                                              
CY20626 CALL DOOUTY( 626 )                                              
C                                                                       
C   NOW FAKE BACK... THIS PART IS EASY !                                
C                                                                       
      TEMP=27+DIR1                                                      
CX10627 CALL DOINX( 627 )                                               
CY10627 CALL DOINY( 627 )                                               
      DO 3 I=0,R1-1                                                     
        PSTR(CNT) = TEMP                                                
    3   CNT=CNT+1                                                       
CX20628 CALL DOOUTX( 628 )                                              
CY20628 CALL DOOUTY( 628 )                                              
      TEMP=27+DIR2                                                      
CX10629 CALL DOINX( 629 )                                               
CY10629 CALL DOINY( 629 )                                               
      DO 4 I=0,R2-1                                                     
        PSTR(CNT) = TEMP                                                
    4   CNT=CNT+1                                                       
CX20630 CALL DOOUTX( 630 )                                              
CY20630 CALL DOOUTY( 630 )                                              
      PSTR(CNT) = 14                                                    
      CNT=CNT+1                                                         
C                                                                       
C   NOW GO FETCH THE LEG. IF R1+R2 = 1 THEN WE HAVE TO DO               
C   THE LEG ON THE FLY                                                  
C                                                                       
      TEMP=23+DIR1                                                      
CX10631 CALL DOINX( 631 )                                               
CY10631 CALL DOINY( 631 )                                               
      DO 5 I=0,R1-1                                                     
        PSTR(CNT) = TEMP                                                
    5   CNT=CNT+1                                                       
CX20632 CALL DOOUTX( 632 )                                              
CY20632 CALL DOOUTY( 632 )                                              
      TEMP=23+DIR2                                                      
CX10633 CALL DOINX( 633 )                                               
CY10633 CALL DOINY( 633 )                                               
      DO 6 I=0,R2-1                                                     
        PSTR(CNT) = TEMP                                                
    6   CNT=CNT+1                                                       
CX20634 CALL DOOUTX( 634 )                                              
CY20634 CALL DOOUTY( 634 )                                              
      IF((R1+R2).NE.1) THEN                                             
        PSTR(CNT) = 36                                                  
        CNT=CNT+1                                                       
      ELSE                                                              
CX10635 CALL DOINX( 635 )                                               
CY10635 CALL DOINY( 635 )                                               
        DO 7 I=0,T-1                                                    
          PSTR(CNT) = 4                                                 
    7     CNT=CNT+1                                                     
CX20636 CALL DOOUTX( 636 )                                              
CY20636 CALL DOOUTY( 636 )                                              
CX10637 CALL DOINX( 637 )                                               
CY10637 CALL DOINY( 637 )                                               
        DO 8 I=0,T-1                                                    
          PSTR(CNT) = 30                                                
    8     CNT=CNT+1                                                     
CX20638 CALL DOOUTX( 638 )                                              
CY20638 CALL DOOUTY( 638 )                                              
      ENDIF                                                             
      TEMP=27+DIR1                                                      
CX10639 CALL DOINX( 639 )                                               
CY10639 CALL DOINY( 639 )                                               
      DO 9 I=0,R1-1                                                     
        PSTR(CNT) = TEMP                                                
    9   CNT=CNT+1                                                       
CX20640 CALL DOOUTX( 640 )                                              
CY20640 CALL DOOUTY( 640 )                                              
      TEMP=27+DIR2                                                      
CX10641 CALL DOINX( 641 )                                               
CY10641 CALL DOINY( 641 )                                               
      DO 10 I=0,R2-1                                                    
        PSTR(CNT) = TEMP                                                
   10   CNT=CNT+1                                                       
CX20642 CALL DOOUTX( 642 )                                              
CY20642 CALL DOOUTY( 642 )                                              
      PSTR(CNT) = 14                                                    
      CNT=CNT+1                                                         
C                                                                       
C   NOW FETCH THE TOP PART OF THE LOOP.                                 
C                                                                       
      TEMP=23+DIR1                                                      
CX10643 CALL DOINX( 643 )                                               
CY10643 CALL DOINY( 643 )                                               
      DO 11 I=0,R1-1                                                    
        PSTR(CNT) = TEMP                                                
   11   CNT=CNT+1                                                       
CX20644 CALL DOOUTX( 644 )                                              
CY20644 CALL DOOUTY( 644 )                                              
      TEMP=23+DIR2                                                      
CX10645 CALL DOINX( 645 )                                               
CY10645 CALL DOINY( 645 )                                               
      DO 12 I=0,R2-1                                                    
        PSTR(CNT) = TEMP                                                
   12   CNT=CNT+1                                                       
CX20646 CALL DOOUTX( 646 )                                              
CY20646 CALL DOOUTY( 646 )                                              
CX10647 CALL DOINX( 647 )                                               
CY10647 CALL DOINY( 647 )                                               
      DO 13 I=0,T-1                                                     
        PSTR(CNT) = 26                                                  
   13   CNT=CNT+1                                                       
CX20648 CALL DOOUTX( 648 )                                              
CY20648 CALL DOOUTY( 648 )                                              
C                                                                       
      NUM1=0                                                            
      NUM2=0                                                            
   14 CONTINUE                                                          
      IF((NUM1.LT.R1).AND.(NUM2.LT.R2)) THEN                            
        TEMP=5+DIR2                                                     
        PSTR(CNT) = TEMP                                                
        CNT=CNT+1                                                       
        NUM2=NUM2+1                                                     
        IF((NUM1+NUM2).LT.(R1+R2-1)) THEN                               
          TEMP=19+DIR1                                                  
        ELSE                                                            
          TEMP=5+DIR1                                                   
        ENDIF                                                           
        PSTR(CNT) = TEMP                                                
        CNT=CNT+1                                                       
        NUM1=NUM1+1                                                     
        GOTO 14                                                         
      ENDIF                                                             
C                                                                       
      IF(NUM1.LT.R1) THEN                                               
        DIR=DIR1                                                        
      ELSE                                                              
        DIR=DIR2                                                        
      ENDIF                                                             
CX10649 CALL DOINX( 649 )                                               
CY10649 CALL DOINY( 649 )                                               
      DO 15 R=(NUM1+NUM2),(R1+R2-1)                                     
        IF((DIR.EQ.DIR1).AND.(R.LT.(R1+R2-1)).AND.(R.GT.0)) THEN        
          TEMP=19+DIR                                                   
        ELSE                                                            
          TEMP=5+DIR                                                    
        ENDIF                                                           
        PSTR(CNT) = TEMP                                                
   15   CNT=CNT+1                                                       
CX20650 CALL DOOUTX( 650 )                                              
CY20650 CALL DOOUTY( 650 )                                              
CX10651 CALL DOINX( 651 )                                               
CY10651 CALL DOINY( 651 )                                               
C                                                                       
      DO 16 I=0,T-1                                                     
        PSTR(CNT) = 30                                                  
   16   CNT=CNT+1                                                       
CX20652 CALL DOOUTX( 652 )                                              
CY20652 CALL DOOUTY( 652 )                                              
C                                                                       
      PSTR(CNT) = 14                                                    
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
                                                                        
                                                                        
                                                                        
C                                                                       
C   J. W. FLOWER  OCT 1985                                              
C                                                                       
C   THIS CODE MEASURES THE OFF AXIS WILSON LOOPS.                       
C                                                                       
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,           
C   AND DAVID WALKER, 1987                                              
C                                                                       
      SUBROUTINE ROTMEA(DIR,NUMHIT,R1MAX,R2MAX)                         
      INTEGER DIR,NUMHIT,R1MAX,R2MAX                                    
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      PARAMETER ( MBAR=18*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL UBAR(MBAR),SBAR(MBAR)                                        
      COMMON /COMUBR/ UBAR,SBAR                                         
C                                                                       
      INTEGER TMAX,T,PLEGS,CNT                                          
      INTEGER COORD(4),NCHAR,TFLAG,SITE,TMP                             
      REAL BOT(18),TOP(18),UP(18),DOWN(18),STEMP(18,2)                  
      REAL TEMP                                                         
      INTEGER PSTR(256),LPSTR(256)                                      
CX30057 CALL SBINX (57)                                                 
CY30057 CALL SBINY (57)                                                 
C                                                                       
      CALL LINKBR(NUMHIT,DIR,SBAR)                                      
CX10653 CALL DOINX( 653 )                                               
CY10653 CALL DOINY( 653 )                                               
C                                                                       
      DO 2 ICORD4=0,LATT(4)-1                                           
        COORD(4) = ICORD4                                               
        DO 2 T=1,12                                                     
          TFLAG=0                                                       
          DO 3 I=1,R1MAX                                                
            DO 3 J=1,R2MAX                                              
              IF(WHICH(I,J).NE.0) THEN                                  
                TMAX = MIN(12,WHICH(I,J)+5)                             
                IF((WHICH(I,J).LE.T).AND.(T.LE.TMAX)) TFLAG = 1         
              ENDIF                                                     
    3         CONTINUE                                                  
          IF(TFLAG.NE.1) GOTO2                                          
          CALL LEGMAK(T,ICORD4)                                         
C                                                                       
C    NOW LOOP THROUGH THE REQUIRED LOOPS.                               
C                                                                       
          DO 4 I=1,R1MAX                                                
            DO 4 J=1,R2MAX                                              
              IF(WHICH(I,J).NE.0) THEN                                  
                WLOOP(I,J) = 0.                                         
                TMAX = MIN(12,WHICH(I,J)+5)                             
                IF((T.GE.WHICH(I,J)).AND.(T.LE.TMAX)) THEN              
                  CALL ROTLPS(LPSTR,I,J,T,DIR,MOD((DIR+1),3))           
                  SITE = COORD(4)*MOV(4)                                
                  PLEGS=1                                               
                  DO 5 ICORD3=0,LATT(3)-1                               
                    COORD(3) = ICORD3                                   
                    DO 5 ICORD2=0,LATT(2)-1                             
                      COORD(2) = ICORD2                                 
                      DO 5 ICORD1=0,LATT(1)-1                           
                        COORD(1) = ICORD1                               
                        CALL UDAG(LEGS(PLEGS),DOWN)                     
C                        CALL CPYMAT(PSTR,LPSTR,256) LPOINTER 4-3-89    
                        CALL ICPYMT(PSTR,LPSTR,256)                     
                        CNT=1                                           
                        CALL SYSLOP(COORD,SITE,PSTR(CNT),BOT,NCHAR)     
                        CNT=CNT+NCHAR                                   
                        CALL SYSLOP(COORD,SITE,PSTR(CNT),UP,NCHAR)      
                        CNT=CNT+NCHAR                                   
                        CALL SYSLOP(COORD,SITE,PSTR(CNT),TOP,NCHAR)     
                        CALL MULT(BOT,UP,STEMP(1,1))                    
                        CALL MULT(STEMP(1,1),TOP,STEMP(1,2))            
                        CALL MULT(STEMP(1,2),DOWN,STEMP(1,1))           
                        WLOOP(I,J)=WLOOP(I,J)+(STEMP(1,1)+STEMP(9,1)    
     .                             +STEMP(17,1))/3.0                    
                        PLEGS=PLEGS+18                                  
                        SITE = SITE + MOV(1)                            
    5                   CONTINUE                                        
                  ENDIF                                                 
                ENDIF                                                   
    4           CONTINUE                                                
C                                                                       
            DO 6 I=1,R1MAX                                              
              DO 6 J=1,R2MAX                                            
                IF(WHICH(I,J).NE.0) THEN                                
                  TEMP = WLOOP(I,J)                                     
                  TMP  = T-WHICH(I,J)+1                                 
                  RES(I,J,TMP) = RES(I,J,TMP) + TEMP                    
                ENDIF                                                   
    6           CONTINUE                                                
    2     CONTINUE                                                      
CX20654 CALL DOOUTX( 654 )                                              
CY20654 CALL DOOUTY( 654 )                                              
CX40058 CALL SBOUTX (58)                                                
CY40058 CALL SBOUTY (58)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE SETCOL(LATT)                                           
      INTEGER LATT(4)                                                   
C                                                                       
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( LISMAX=(L1*L2*L3*L4)/2+1 )                            
      COMMON /COLCOM/ LISRED(LISMAX),LISBLK(LISMAX),NRED,NBLACK         
CX30059 CALL SBINX (59)                                                 
CY30059 CALL SBINY (59)                                                 
C                                                                       
      NRED   = 0                                                        
      NBLACK = 0                                                        
      IPTR   = 0                                                        
      INDEX4 = 0                                                        
CX10655 CALL DOINX( 655 )                                               
CY10655 CALL DOINY( 655 )                                               
      DO 4 I=1,LATT(4)                                                  
         INDEX3 = INDEX4                                                
         DO 3 J=1,LATT(3)                                               
            INDEX2 = INDEX3                                             
            DO 2 K=1,LATT(2)                                            
               INDEX1 = INDEX2                                          
               DO 1 L=1,LATT(1)                                         
                  IF (INDEX1.EQ.0) THEN                                 
                     NRED = NRED+1                                      
                     LISRED(NRED) = IPTR                                
                  ELSE                                                  
                     NBLACK = NBLACK+1                                  
                     LISBLK(NBLACK) = IPTR                              
                  ENDIF                                                 
                  IPTR   = IPTR+72                                      
    1             INDEX1 = 1-INDEX1                                     
    2          INDEX2 = 1-INDEX2                                        
    3       INDEX3 = 1-INDEX3                                           
    4    INDEX4 = 1-INDEX4                                              
CX20656 CALL DOOUTX( 656 )                                              
CY20656 CALL DOOUTY( 656 )                                              
CX40060 CALL SBOUTX (60)                                                
CY40060 CALL SBOUTY (60)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C                                                                       
C   STEVE OTTO AND JON FLOWER   JULY 1984.                              
C                                                                       
C   RECEIVES  PARAMETER INPUT FOR THE SU3 ROUTINES.                     
C                                                                       
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,           
C   AND DAVID WALKER, 1987                                              
C                                                                       
      SUBROUTINE SETPAR                                                 
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      INTEGER NUMSTG                                                    
      COMMON /COMFIL/ NUMSTG                                            
CX30061 CALL SBINX (61)                                                 
CY30061 CALL SBINY (61)                                                 
C                                                                       
    1 WRITE(6,*) 'ENTER X,Y,Z,T LATTICE SIZE ?'                         
      READ (5,*) LATT                                                   
      WRITE(6,*) LATT                                                   
      NUM = LATT(1)*LATT(2)*LATT(3)*LATT(4)                             
C                                                                       
      IF(NUM*72.GT.UMAX) THEN                                           
        WRITE(6,*) 'NOT ENOUGH ROOM FOR GAUGE FIELD'                    
        WRITE(6,*) '*******************************'                    
        WRITE(6,*) 'TRY AGAIN'                                          
        GOTO 1                                                          
      ENDIF                                                             
C                                                                       
C      NUMSIT = LATT(1)*LATT(2)*LATT(3)*LATT(4)                         
      NUMSTG = LATT(1)*LATT(2)*LATT(3)*LATT(4)                          
      WRITE(6,*) 'ENTER VALUE OF BETA '                                 
      READ (5,*) BETA                                                   
      WRITE(6,*) BETA                                                   
C                                                                       
      CALL MAKE()                                                       
C                                                                       
      CALL SETCOL(LATT)                                                 
C                                                                       
      WRITE(6,*) 'PARAMETERS NOW SET.'                                  
      WRITE(6,*) 'NOW NEED TO INITIALIZE..........'                     
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE SETSED()                                               
C                                                                       
      INTEGER SEED,MODE                                                 
CX30063 CALL SBINX (63)                                                 
CY30063 CALL SBINY (63)                                                 
C                                                                       
      WRITE(6,*) 'SELECT MODE OF SEED SETTING: '                        
      WRITE(6,*) '0  <->  BY HAND'                                      
      WRITE(6,*) '1  <->  READ IN FROM STANDARD FILE, SEED'             
      READ (5,*) MODE                                                   
      IF( MODE.EQ.0) THEN                                               
        WRITE(6,*) 'INPUT SEED'                                         
        READ (5,*) SEED                                                 
        CALL PRNSET(SEED)                                               
      ELSE                                                              
        OPEN(UNIT=1,FILE='SEED',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')   
        REWIND(1)                                                       
        READ(1,10) SEED                                                 
   10   FORMAT(I12)                                                     
        CALL PRNSET(SEED)                                               
        CLOSE(1)                                                        
      ENDIF                                                             
      WRITE(6,*) 'SEED SET'                                             
CX40064 CALL SBOUTX (64)                                                
CY40064 CALL SBOUTY (64)                                                
      RETURN                                                            
      END                                                               
C                                                                       
CCCCCCCCCCCCCCCCCCCCCCCCC END OF ROUTINE SETSED CCCCCCCCCCCCCCCCCCCC    
C                                                                       
                                                                        
                                                                        
                                                                        
                                                                        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C                                                                     C 
C          THIS SOFTWARE COPYRIGHTED BY CALTECH, 1983, 1986, 1987     C 
C                                                                     C 
C          PURE GAUGE SU(3) MONTE CARLO                               C 
C                                                                     C 
C          THIS CODE HAS A LONG HISTORY:                              C 
C              EARLY VERSIONS OF THIS DEVELOPED FOR THE MARK I        C 
C              HYPERCUBE BY STEVE OTTO AND PAUL STOLORZ (1981-83)     C 
C              JON FLOWER AND STEVE OTTO EXTENDED AND IMPROVED THE    C 
C              CODE FOR THE MARK II MACHINES. (1984-86)               C 
C              THE CURRENT INCARNATION IS A "PORT" TO THE NCUBE       C 
C              DONE BY CLIVE BAILLIE, ED FELTEN, STEVE OTTO,          C 
C              AND DAVID WALKER  (1987)                               C 
C                                                                     C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C                                                                       
C       SYSLOP SYSTEMATICALLY CALCULATES WILSON LOOPS                   
C       FOR SU3 THEORY IN 3+1 DIM WITH THE INT STRINGS                  
C       AS INPUT, ALSO KEEPS TRACK OF THE LENGTH OF THE STRING.         
C                                                                       
C       CODING OF LATTICE PATHS:                                        
C                1, 2, 3, 4 +X,+Y,+Z,+T                                 
C                5, 6, 7, 8 -X,-Y,-Z,-T                                 
C               15,16,17,18  X,+Y,+Z,+T USING BARS.                     
C               19,20,21,22 -X,-Y,-Z,-T USING BARS.                     
C               23,24,25,26 +X,+Y,+Z,+T FAKE STEPS                      
C               27,28,29,30 -X,-Y,-Z,-T FAKE STEPS                      
C               36          RETURNS LEG                                 
C                                                                       
C       CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,       
C       AND DAVID WALKER, 1987                                          
C                                                                       
      SUBROUTINE SYSLOP(COORD,SITE,PTR,PLOOP,NN)                        
      INTEGER COORD(4),NN,SITE,SITES                                    
      INTEGER PTR(*)                                                    
      REAL PLOOP(18)                                                    
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      PARAMETER ( MBAR=18*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      REAL UBAR(MBAR),SBAR(MBAR)                                        
      COMMON /COMUBR/ UBAR,SBAR                                         
C                                                                       
      REAL FTEMP(18,2),FILMAT(18),EXTRA(18)                             
      INTEGER PU,PBAR                                                   
      INTEGER PLEGS,TINDEX,SETFLG                                       
CX30065 CALL SBINX (65)                                                 
CY30065 CALL SBINY (65)                                                 
C                                                                       
      SITES  = SITE                                                     
      INDEX  = 0                                                        
      SETFLG = 0                                                        
      NN     = 0                                                        
C                                                                       
    1 NN = NN+1                                                         
      IC = PTR(NN)                                                      
      IF(IC.EQ.14) GOTO 2                                               
C                                                                       
      IF(IC.LE.4) THEN                                                  
        IND = IC                                                        
        PU  = SITES+ROT(IND)                                            
        IF(SETFLG.EQ.0) THEN                                            
          CALL CPYMAT(FTEMP(1,INDEX+1),U1(PU+1),18)                     
          SETFLG = 1                                                    
        ELSE                                                            
          CALL CPYMAT(FILMAT,U1(PU+1),18)                               
          TINDEX = 1-INDEX                                              
          CALL MULT(FTEMP(1,INDEX+1),FILMAT,FTEMP(1,TINDEX+1))          
          INDEX  = TINDEX                                               
        ENDIF                                                           
        COORD(IND) = MOD( COORD(IND)+1, LATT1(IND)+1 )                  
        IF(COORD(IND).EQ.0) THEN                                        
          SITES = SITES-MOV(IND)*LATT1(IND)                             
        ELSE                                                            
          SITES = SITES+MOV(IND)                                        
        ENDIF                                                           
C                                                                       
      ELSEIF(IC.LE.8) THEN                                              
        IND = IC - 4                                                    
        IF(COORD(IND).EQ.0) THEN                                        
          COORD(IND) = LATT1(IND)                                       
          SITES = SITES+MOV(IND)*LATT1(IND)                             
        ELSE                                                            
          COORD(IND) = COORD(IND)-1                                     
          SITES = SITES-MOV(IND)                                        
        ENDIF                                                           
        PU = SITES+ROT(IND)                                             
        IF(SETFLG.EQ.0) THEN                                            
          CALL CPYMAT(FILMAT,U1(PU+1),18)                               
          CALL UDAG(FILMAT,FTEMP(1,INDEX+1))                            
          SETFLG = 1                                                    
        ELSE                                                            
          CALL CPYMAT(EXTRA,U1(PU+1),18)                                
          CALL UDAG(EXTRA,FILMAT)                                       
          TINDEX = 1-INDEX                                              
          CALL MULT(FTEMP(1,INDEX+1),FILMAT,FTEMP(1,TINDEX+1))          
          INDEX  = TINDEX                                               
        ENDIF                                                           
C                                                                       
      ELSEIF(IC.EQ.36) THEN                                             
        LSITE = SITES - COORD(4)*MOV(4)                                 
        PLEGS = (LSITE/72)*18                                           
CX10657 CALL DOINX( 657 )                                               
CY10657 CALL DOINY( 657 )                                               
        DO 4 I=1,18                                                     
    4     FTEMP(I,INDEX+1) = LEGS(PLEGS+I)                              
CX20658 CALL DOOUTX( 658 )                                              
CY20658 CALL DOOUTY( 658 )                                              
        SETFLG = 1                                                      
C                                                                       
      ELSEIF(IC.LE.18) THEN                                             
        IND = IC - 14                                                   
        PBAR = (SITES/72)*18                                            
        IF(IND.EQ.4) THEN                                               
          CALL CPYMAT(FILMAT,UBAR(PBAR+1),18)                           
        ELSE                                                            
          CALL CPYMAT(FILMAT,SBAR(PBAR+1),18)                           
        ENDIF                                                           
        IF(SETFLG.EQ.0) THEN                                            
CX10659 CALL DOINX( 659 )                                               
CY10659 CALL DOINY( 659 )                                               
          DO 3 J=1,18                                                   
    3       FTEMP(J,INDEX+1) = FILMAT(J)                                
CX20660 CALL DOOUTX( 660 )                                              
CY20660 CALL DOOUTY( 660 )                                              
          SETFLG = 1                                                    
        ELSE                                                            
          TINDEX = 1-INDEX                                              
          CALL MULT(FTEMP(1,INDEX+1),FILMAT,FTEMP(1,TINDEX+1))          
          INDEX  = TINDEX                                               
        ENDIF                                                           
        COORD(IND) = MOD( COORD(IND)+1, LATT1(IND)+1 )                  
        IF(COORD(IND).EQ.0) THEN                                        
          SITES = SITES-MOV(IND)*LATT1(IND)                             
        ELSE                                                            
          SITES = SITES+MOV(IND)                                        
        ENDIF                                                           
C                                                                       
      ELSEIF(IC.LE.22) THEN                                             
        IND = IC - 18                                                   
        IF(COORD(IND).EQ.0) THEN                                        
          COORD(IND) = LATT1(IND)                                       
          SITES = SITES+ MOV(IND)*LATT1(IND)                            
        ELSE                                                            
          COORD(IND) = COORD(IND)-1                                     
          SITES = SITES-MOV(IND)                                        
        ENDIF                                                           
        PBAR = (SITES/72)*18                                            
        IF(IND.EQ.4) THEN                                               
          CALL CPYMAT(EXTRA,UBAR(PBAR+1),18)                            
        ELSE                                                            
          CALL CPYMAT(EXTRA,SBAR(PBAR+1),18)                            
        ENDIF                                                           
        IF(SETFLG.EQ.0) THEN                                            
          CALL UDAG(EXTRA,FTEMP(1,INDEX+1))                             
          SETFLG = 1                                                    
        ELSE                                                            
          CALL UDAG(EXTRA,FILMAT)                                       
          TINDEX = 1-INDEX                                              
          CALL MULT(FTEMP(1,INDEX+1),FILMAT,FTEMP(1,TINDEX+1))          
          INDEX  = TINDEX                                               
        ENDIF                                                           
C                                                                       
      ELSEIF(IC.LE.26) THEN                                             
        IND = IC - 22                                                   
        COORD(IND) = MOD( COORD(IND)+1, LATT1(IND)+1 )                  
        IF(COORD(IND).EQ.0) THEN                                        
          SITES = SITES-MOV(IND)*LATT1(IND)                             
        ELSE                                                            
          SITES = SITES+MOV(IND)                                        
        ENDIF                                                           
C                                                                       
      ELSEIF(IC.LE.30) THEN                                             
        IND = IC - 26                                                   
        IF(COORD(IND).EQ.0) THEN                                        
          COORD(IND) = LATT1(IND)                                       
          SITES = SITES+MOV(IND)*LATT1(IND)                             
        ELSE                                                            
          COORD(IND) = COORD(IND)-1                                     
          SITES = SITES-MOV(IND)                                        
        ENDIF                                                           
C                                                                       
      ENDIF                                                             
      GOTO 1                                                            
C                                                                       
    2 CONTINUE                                                          
CX10661 CALL DOINX( 661 )                                               
CY10661 CALL DOINY( 661 )                                               
C                                                                       
      DO 5 I=1,18                                                       
    5   PLOOP(I) = FTEMP(I,INDEX+1)                                     
CX20662 CALL DOOUTX( 662 )                                              
CY20662 CALL DOOUTY( 662 )                                              
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C                                                                       
C   UDAG EVALUATES THE COMPLEX CONJUGATE OF THE INPUT 3x3 MATRIX        
C                                                                       
      SUBROUTINE UDAG(PIN,POUT)                                         
      REAL PIN(18),POUT(18)                                             
CX30067 CALL SBINX (67)                                                 
CY30067 CALL SBINY (67)                                                 
C                                                                       
      POUT(1)  =  PIN(1)                                                
      POUT(7)  =  PIN(3)                                                
      POUT(13) =  PIN(5)                                                
      POUT(2)  = -PIN(2)                                                
      POUT(8)  = -PIN(4)                                                
      POUT(14) = -PIN(6)                                                
      POUT(3)  =  PIN(7)                                                
      POUT(9)  =  PIN(9)                                                
      POUT(15) =  PIN(11)                                               
      POUT(4)  = -PIN(8)                                                
      POUT(10) = -PIN(10)                                               
      POUT(16) = -PIN(12)                                               
      POUT(5)  =  PIN(13)                                               
      POUT(11) =  PIN(15)                                               
      POUT(17) =  PIN(17)                                               
      POUT(6)  = -PIN(14)                                               
      POUT(12) = -PIN(16)                                               
      POUT(18) = -PIN(18)                                               
CX40068 CALL SBOUTX (68)                                                
CY40068 CALL SBOUTY (68)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
                                                                        
                                                                        
C                                                                       
C     J. W. FLOWER AND STEVE OTTO   JULY 1984.                          
C                                                                       
C     UPDATES THE LATICE "SWEEPS" TIMES WITH PROJECTIONS EVERY "PINVAL" 
C                                                                       
C     CONVERRTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,        
C     AND DAVID WALKER, 1987                                            
C                                                                       
      SUBROUTINE UPDATE(SWEEPS,PINVAL)                                  
      INTEGER SWEEPS,PINVAL                                             
C                                                                       
      INTEGER DIMS                                                      
      PARAMETER (DIMS=4)                                                
      PARAMETER ( L1=12, L2=12, L3=12, L4=12 )                          
      PARAMETER ( MLEG=18*L1*L2*L3 )                                    
      INTEGER UMAX                                                      
      PARAMETER ( UMAX=72*L1*L2*L3*L4 )                                 
      INTEGER LATT(DIMS),LATT1(DIMS),MOV(DIMS),ROT(DIMS),WHICH(13,13)   
      REAL WLOOP(13,13),RES(13,13,6),BETA,LEGS(MLEG)                    
      COMMON /COMSYS/ BETA,WLOOP,RES,LEGS,LATT,LATT1,MOV,ROT,WHICH      
      PARAMETER ( LISMAX=(L1*L2*L3*L4)/2+1 )                            
      COMMON /COLCOM/ LISRED(LISMAX),LISBLK(LISMAX),NRED,NBLACK         
C                                                                       
C                                                                       
      REAL U1(UMAX)                                                     
      COMMON /COMU1/ U1                                                 
C                                                                       
      INTEGER COORD(4),DIR,ORIGIN                                       
CX30069 CALL SBINX (69)                                                 
CY30069 CALL SBINY (69)                                                 
CX10663 CALL DOINX( 663 )                                               
CY10663 CALL DOINY( 663 )                                               
C                                                                       
      DO 1000 LOOP=0,SWEEPS-1                                           
C                                                                       
C  FIRST DO RED POINTS                                                  
C                                                                       
CX10665 CALL DOINX( 665 )                                               
CY10665 CALL DOINY( 665 )                                               
        DO 1 I=1,NRED                                                   
           ORIGIN = LISRED(I)                                           
           INDEX  = ORIGIN/72                                           
           CALL GETPOS(LATT,COORD,INDEX)                                
           DO 10 DIR=0,3                                                
              CALL HIT(COORD,ORIGIN,DIR)                                
   10         CONTINUE                                                  
    1      CONTINUE                                                     
CX10666 CALL DOOUTX( 666 )                                               
CY10666 CALL DOOUTY( 666 )                                               
C                                                                       
C  NOW DO BLACK POINTS                                                  
C                                                                       
CX10667 CALL DOINX( 667 )                                               
CY10667 CALL DOINY( 667 )                                               
        DO 2 I=1,NBLACK                                                 
           ORIGIN = LISBLK(I)                                           
           INDEX  = ORIGIN/72                                           
           CALL GETPOS(LATT,COORD,INDEX)                                
           DO 20 DIR=0,3                                                
              CALL HIT(COORD,ORIGIN,DIR)                                
   20         CONTINUE                                                  
    2      CONTINUE                                                     
CX10668 CALL DOOUTX( 668 )                                               
CY10668 CALL DOOUTY( 668 )                                               
C                                                                       
        IF(MOD(LOOP,PINVAL).EQ.(PINVAL-1)) CALL PROJEC()                
 1000   CONTINUE                                                        
CX20664 CALL DOOUTX( 664 )                                              
CY20664 CALL DOOUTY( 664 )                                              
CX40070 CALL SBOUTX (70)                                                
CY40070 CALL SBOUTY (70)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
