C     ************
      PROGRAM QCD1
C     ************
C
C
C     THIS PROGRAM PERFORMS A MONTE-CARLO-SIMULATION OF THE (3+1)-DIMENSIONAL
C     SU3 LATTICE GAUGE THEORIE. THE  GAUGE VARIABLES ARE UPDATED  USING THE
C     METROPOLIS ALGORITHM.
C
C
C     LATEST REVISION NOV. 1989; E.KEHL
C
C
C     NAMING OF PARAMETERS:
C
C
C     NX    : NUMBER OF LATTICE POINTS I X-DIMENSION
C     NY    : NUMBER OF LATTICE POINTS I X-DIMENSION
C     NZ    : NUMBER OF LATTICE POINTS I Z-DIMENSION
C     NT    : NUMBER OF LATTICE POINTS I T-DIMENSION
C     NV    : TOTAL NUMBER OF LATTICE POINTS OF THE INTERIOR LATTICE
C     NVH   : HALF THE NUMBER OF LATTICE POINTS (ONE CHESS BOARD COLOUR)
C     NXYZ  : NUMBER OF LATTICE POINTS OF THE 3-D HYPERPLANE T=C
C     NXYT  : NUMBER OF LATTICE POINTS OF THE 3-D HYPERPLANE Z=C
C     NXZT  : NUMBER OF LATTICE POINTS OF THE 3-D HYPERPLANE Y=C
C     NYZT  : NUMBER OF LATTICE POINTS OF THE 3-D HYPERPLANE X=C
C     NXY   : NUMBER OF LATTICE POINTS OF THE 2-D HYPERPLANE Z=C, T=C
C     NXZ   : NUMBER OF LATTICE POINTS OF THE 2-D HYPERPLANE Y=C, T=C
C     NXT   : NUMBER OF LATTICE POINTS OF THE 2-D HYPERPLANE Y=C, Z=C
C     NYZ   : NUMBER OF LATTICE POINTS OF THE 2-D HYPERPLANE X=C, T=C
C     NYT   : NUMBER OF LATTICE POINTS OF THE 2-D HYPERPLANE X=C, Z=C
C     NZT   : NUMBER OF LATTICE POINTS OF THE 2-D HYPERPLANE X=C, Y=C
C     NXE   : NUMBER OF LATTICE POINTS IN X-DIMENSION PLUS TWO
C     NYE   : NUMBER OF LATTICE POINTS IN Y-DIMENSION PLUS TWO
C     NZE   : NUMBER OF LATTICE POINTS IN Z-DIMENSION PLUS TWO
C     NTE   : NUMBER OF LATTICE POINTS IN T-DIMENSION PLUS TWO
C     NXYE  : NUMBER OF LATTICE POINTS OF THE 2-HYP. Z=C,T=C EXTENDED
C     NXYZE : NUMBER OF LATTICE POINTS OF THE 3-HYP.T=C EXTENDED
C     NVE   : NUMBER OF LATTICE POINTS OF THE EXTENDED LATTICE
C
C
C     NAMING OF LOOP VARIABLES:
C
C
C     ISD   : SEND DIRECTION             ;  1,2,3,4 (X,Y,Z,T)
C     IBF   : SEND BACKWARD OR FORWARD   ;  1,2     (B,F)
C     ILD   : LINK DIRECTION             ;  1,2,3,4 (X,Y,Z,T)
C     ICC   : CHESS BOARD COLOUR         ;  1,2     (BLACK,WHITE)
C     ISM   : COMPONENT OF SU(3) MATRIX  ;  1,...,18
C
C
C     NAMING OF ARRAYS AND VARIABLES
C
C
C     LINKS     : SU3 LINK VARIABLES OF THE LATTICE
C     WKI       : WORK ARRAYS FOR TEMPORARY RESULTS
C     TBL       : TABLE OF SU3 MATRICES CLOSE TO THE UNIT MATRIX
C     GCB       : INDEX LIST FOR GATHERING ONE CHESS  BOARD COLOUR
C     IRAND     : ARRAY OF INTEGER RANDOM NUMBERS
C
C
      PARAMETER( NR1   =   10                         )
      PARAMETER( NW1   =   6                         )
C
      INCLUDE 'qcd1.inc'
C
      PARAMETER( NX    =   NS                         )
      PARAMETER( NY    =   NS                         )
      PARAMETER( NZ    =   NS                         )
      PARAMETER( NV    =   NX  * NY * NZ * NT         )
      PARAMETER( NXYZ  =   NX  * NY * NZ              )
      PARAMETER( NVH   = ( NV            + 1 ) / 2    )
      PARAMETER( NXYZH = ( NX  * NY * NZ + 1 ) / 2    )
      PARAMETER( NXYTH = ( NX  * NY * NT + 1 ) / 2    )
      PARAMETER( NXZTH = ( NX  * NZ * NT + 1 ) / 2    )
      PARAMETER( NYZTH = ( NY  * NZ * NT + 1 ) / 2    )
      PARAMETER( NXYH  = ( NX  * NY      + 1 ) / 2    )
      PARAMETER( NXZH  = ( NX  * NZ      + 1 ) / 2    )
      PARAMETER( NXTH  = ( NX  * NT      + 1 ) / 2    )
      PARAMETER( NYZH  = ( NY  * NZ      + 1 ) / 2    )
      PARAMETER( NYTH  = ( NY  * NT      + 1 ) / 2    )
      PARAMETER( NZTH  = ( NZ  * NT      + 1 ) / 2    )
      PARAMETER( NXE   =   NX  + 2                    )
      PARAMETER( NYE   =   NY  + 2                    )
      PARAMETER( NZE   =   NZ  + 2                    )
      PARAMETER( NTE   =   NT  + 2                    )
      PARAMETER( NXYE  =   NXE * NYE                  )
      PARAMETER( NXYZE =   NXE * NYE * NZE            )
      PARAMETER( NVE   =   NXE * NYE * NZE * NTE      )
C
C
      DOUBLE PRECISION  LINKS( 1:NVE, 1:18, 1:4 )
C
      DOUBLE PRECISION  WK1( 1:NVH, 1:18 )
      DOUBLE PRECISION  WK2( 1:NVH, 1:18 )
      DOUBLE PRECISION  WK3( 1:NVH, 1:18 )
      DOUBLE PRECISION  WK4( 1:NVH, 1:18 )
      DOUBLE PRECISION  WK5( 1:NVH, 1:18 )
C
      DOUBLE PRECISION  TBL  ( 1:NVH, 1:2, 1:18 )
C
      DOUBLE PRECISION  BETA,SIGMA,
     >                  WLR,WLI,
     >                  WLSUMR,WLSUMI,
     >                  ZERO,RTIME,TSTART,TEND
C
C
      INTEGER  GCB  ( 1:NVH, 1:2  )
      INTEGER  GTS  ( 1:NXYZ      )
      INTEGER  IRAND( 1:NVH       )
C
      INTEGER  GX( 1:NYZTH, 1:2, 1:4 )
      INTEGER  GY( 1:NXZTH, 1:2, 1:4 )
      INTEGER  GZ( 1:NXYTH, 1:2, 1:4 )
      INTEGER  GT( 1:NXYZH, 1:2, 1:4 )
C
      INTEGER  GXY( 1:NZTH, 1:2, 1:4 )
      INTEGER  GXZ( 1:NYTH, 1:2, 1:4 )
      INTEGER  GXT( 1:NYZH, 1:2, 1:4 )
      INTEGER  GYZ( 1:NXTH, 1:2, 1:4 )
      INTEGER  GYT( 1:NXZH, 1:2, 1:4 )
      INTEGER  GZT( 1:NXYH, 1:2, 1:4 )
C
C
      INTEGER  LX(2,4),LY(2,4),
     >         LZ(2,4),LT(2,4)
C
      INTEGER  LXY(2,4),LXZ(2,4),
     >         LXT(2,4),LYZ(2,4),
     >         LYT(2,4),LZT(2,4)
C
      INTEGER  LCB(2)
      INTEGER  SHT(4)
      INTEGER  SEED,START,ORDEXP
C
C
      LOGICAL     RESULT
C
C
      PARAMETER( ZERO = 0.D0 )
C
      COMMON / PERF / RTIME(15),IPRO
C
      DATA ISU3 / 10 /
C
C
C     READ AND CHECK INPUT PARAMETERS
C
C
      iccos = 0
c
      OPEN( NR1, FILE = 'qcd1.dat' )
c--      OPEN( NW1, FILE = 'result' )
      REWIND( NR1 )
c--      REWIND( NW1 )
C
C
      READ( NR1, 1000 ) START,IAVG,IAVGDD,SEED,ITEMAX,
     >                  METHIT,ITBLUP,ORDEXP,BETA,SIGMA,
     >                  IPRO
 1000 FORMAT(7(1X/),8(57X,I6/),2(57X,F13.6/),57X,I6)
C
C
C     **********
      CALL CHKIN(
C     **********
     >           NS,NT,NW1,
     >           START,IAVG,IAVGDD,SEED,ITEMAX,
     >           METHIT,ITBLUP,ORDEXP,BETA,SIGMA)
C
C
C     **********
      CALL HEAD(
C     **********
     >          NS,NT,NW1,
     >          START,IAVG,IAVGDD,SEED,ITEMAX,
     >          METHIT,ITBLUP,ORDEXP,BETA,SIGMA )
C
C
      SHT(1) = 1
      SHT(2) = NXE
      SHT(3) = NXYE
      SHT(4) = NXYZE
C
C
C     INITIALIZATION
C
C
C     *********
      CALL INIT(
C     *********
     >          NX,NY,NZ,NT,
     >          NXYZH,NXYTH,NXZTH,NYZTH,
     >          NXYH,NXZH,NXTH,NYZH,NYTH,NZTH,
     >          NXE,NXYE,NXYZE,
     >          NVH,NVE,NXYZ,
     >          SEED,START,
     >          LINKS,WK1,
     >          GCB,GTS,LCB,
     >          GX,GY,GZ,GT,
     >          GXY,GXZ,GXT,GYZ,GYT,GZT,
     >          LX,LY,LZ,LT,
     >          LXY,LXZ,LXT,LYZ,LYT,LZT)
C
C
C     **********
      CALL MKTBL(NVH,ORDEXP,SIGMA,TBL,WK1,WK2,WK3,WK4)
C     **********
C
C
C     EXCHANGING SURFACE VARIABLES
C
C
      DO 20 JCC = 1,2
C
         ICC = JCC + ICCOS
         IF( ICC .EQ. 3 ) ICC = 1
         ICCOS = ICC - 1
C
      DO 20 ILD = 1,4
C
C
C     ***********
      CALL SUREXC(
C     ***********
     >            NXYZH,NXYTH,NXZTH,NYZTH,
     >            NXYH,NXZH,NXTH,NYZH,NYTH,NZTH,
     >            ICC,ILD,
     >            LINKS,WK1,NVE,NVH,
     >            GX,GY,GZ,GT,
     >            GXY,GXZ,GXT,GYZ,GYT,GZT,
     >            LX,LY,LZ,LT,
     >            LXY,LXZ,LXT,LYZ,LYT,LZT)
C
C
   20 CONTINUE
C
C
C     INITIALIZATION
C
C
      ITSAVG = 0
      WLSUMR = ZERO
      WLSUMI = ZERO
      ISU3M  = ISU3
      ISU3R  = ISU3
C
C
C     THE ITERATION SWEEP
C
C
      DO 15 I1 = 1,15
         RTIME(I1) = 0.0D0
 15   CONTINUE
C
C
      CALL TIMER(RTIME(1))
C
C
      DO 30 ITS = 1,ITEMAX
C
      IACPT = 0
C
C
C     UPDATEING  THE TABLE OF RANDOM SU3 MATRICES
C
C
      IF( MOD(ITS,ITBLUP) .EQ. 0 ) THEN
C
          CALL TIMER(TSTART)
C
C
C         **********
          CALL MKTBL(NVH,ORDEXP,SIGMA,TBL,WK1,WK2,WK3,WK4)
C         **********
C
C
          CALL TIMER(TEND)
          RTIME(2) = RTIME(2) + TEND - TSTART
C
      ENDIF
C
C
C     UPDATE OF LINK VARIABLES
C
C
C     LOOP OVER TWO CHESS-BOARD COLOURS
C
C
      DO 40 JCC = 1,2
C
C
      ICC = JCC + ICCOS
      IF( ICC .EQ. 3 ) ICC = 1
      ICCOS = ICC - 1
      NVL = LCB( ICC )
C
C
      DO 50 ILD = 1,4
C
C
      CALL TIMER(TSTART)
C     ***********
      CALL STAPLE(
C     ***********
     >            NVH,NVL,NVE,ICC,ILD,
     >            SHT,GCB,LINKS,
     >            WK1,WK2,WK3,WK4 )
      CALL TIMER(TEND)
      RTIME(3) = RTIME(3) + TEND - TSTART
C
C
C     ***********
      CALL GATHER( NVH,NVL,WK1,NVE,LINKS(1,1,ILD),GCB(1,ICC),0 )
C     ***********
C
C
      CALL TIMER(TSTART)
C     **********
      CALL METRO(
C     **********
     >           NVH,NVL,NV,
     >           METHIT,IACPT,BETA,
     >           TBL,
     >           WK1,WK2,WK3,WK4,WK5,
     >           IRAND)
      CALL TIMER(TEND)
      RTIME(4) = RTIME(4) + TEND - TSTART
C
C
C     CHECK LINK VARIABLES ON SU3 CHARACTERISTICS
C
C
      IF( MOD( ITS, ISU3R ) .EQ. 0 ) THEN
C
C
          CALL TIMER(TSTART)
C         ***********
          CALL SU3TST( NVH,NVL,WK1,RESULT,ITS )
C         ***********
C
C
          IF( RESULT ) THEN
C
C
C             ***********
              CALL RENORM( NVH,WK1 )
C             ***********
C
C
              ISU3M = MAX( 1, ISU3M - ISU3 )
          ELSE
              ISU3M = ISU3M + ISU3
C
          ENDIF
          ISU3R = ISU3R + ISU3M
          CALL TIMER(TEND)
          RTIME(5) = RTIME(5) + TEND - TSTART
C
C
      ENDIF
C
C
C     ***********
      CALL SCATER( NVE,NVL,LINKS(1,1,ILD),NVH,WK1,GCB(1,ICC),0 )
C     ***********
C
C
      CALL TIMER(TSTART)
C     ***********
      CALL SUREXC(
C     ***********
     >            NXYZH,NXYTH,NXZTH,NYZTH,
     >            NXYH,NXZH,NXTH,NYZH,NYTH,NZTH,
     >            ICC,ILD,
     >            LINKS,WK1,NVE,NVH,
     >            GX,GY,GZ,GT,
     >            GXY,GXZ,GXT,GYZ,GYT,GZT,
     >            LX,LY,LZ,LT,
     >            LXY,LXZ,LXT,LYZ,LYT,LZT)
      CALL TIMER(TEND)
      RTIME(6) = RTIME(6) + TEND - TSTART
C
C
   50 CONTINUE
   40 CONTINUE
C
C
C     CALCULATION OF EXPECTION VALUES
C
C
      CALL TIMER(TSTART)
      IF( MOD(ITS,IAVG) .EQ. 0 .AND.
     >        ITS .GT. IAVGDD ) THEN
C
C
C         **********
          CALL WLINE(
C         **********
     >               NVE,NXYZ,NXYZE,NT,
     >               WLR,WLI,
     >               LINKS,
     >               GTS,
     >               WK1,WK2,WK3)
C
C
          ITSAVG = ITSAVG + 1
C
C         ***********
          CALL EXPVAL(
C         ***********
     >                NV,NXYZ,METHIT,
     >                ITS,ITSAVG,NW1,
     >                WLR,WLI,IACPT,
     >                WLSUMR,WLSUMI)
C
      ENDIF
      CALL TIMER(TEND)
      RTIME(7) = RTIME(7) + TEND - TSTART
C
C
   30 CONTINUE
C
C
      CALL TIMER(TEND)
      RTIME(1) = TEND - RTIME(1)
C
C
C     PROVIDE PROFILING INFORMATION
C
C
C     ***********
      CALL PRFILE(
C     ***********
     >             NV,NW1,ISU3M,
     >             ITEMAX,METHIT)
C
C
      CLOSE( NR1 )
c--      CLOSE( NW1 )
      STOP
C
C
      END
