text:
      PROGRAM SU2COR
****************************************************************
*
* Quantum Physics Program: Computation of masses of elementary
* particles in the framework of the Quark-Gluon theory,
* using a monte carlo method
*
* SPEC Benchmark Program 103.su2cor, adopted from the program "su2cor"
* (SNI-internal name: ASW01) written by B. Bunk, University of
* Wuppertal, Germany, 1987
* Previous SPEC version (in suite CFP92): 089.su2cor
*
* Contact: Dr. Wilfried Stehling, Siemens Nixdorf, D 552
*          81730 Muenchen, Germany
*
* Modifications for the SPEC version:
*
* - Internal timing statements removed or commented out
* - Detailed output supressed
*
* Modifications March 1991 by
*   Aenne Scharbert,
*   Marcus Schwankl,
*   Reinhold Weicker,
*   Andrea Wittmann,
*       Siemens Nixdorf Information Systems, STM OS 323,
*       Otto-Hahn-Ring 6, W-8000 Muenchen 83, Germany
*
* Modifications Nov. 1993 by
*   Wilfried Stehling, Siemens Nixdorf Informationssysteme AG,
*   D 552, 81730 Muenchen, Germany
* Increased problem size for longer execution time:
* Grid size MAXVEC was 4**3*8, is now 8**3*8.
* File su2cor.data.f and input values have been modified accordingly.
*
* Modifications Sept. 1994 by
*   Reinhold Weicker, Siemens Nixdorf
* Decreased number of iterations
* Removed subroutines that are not necessary for SPEC version
*
* Modifications Nov. 1994 by
*   Reinhold Weicker, Siemens Nixdorf
* Increased number of iterations
* Set grid size to 8**3*8+7
*
****************************************************************
C     PROGRAM SU2COR
C***************************************************************
C****    B. BUNK, UMIVERSITAET WUPPERTAL 1987
C****    BERECHNUNG VON GLUEBALL-MASSEN IM RAHMEN DER
C****    SU2 EICHTHEORIE.
C****    MONTE CARLO RECHNUNG, ERZEUGUNG DER KONFIGURATIONEN
C****    MIT DER WAERMEBAD-METHODE.
C****    GITTERGROESSE: 8**3*16 , SPEICHERBEDARF 5 MBYTES
C****    ANZAHL DER ERZEUGTEN KONFIGURATIONEN : 20
C****    DIE ANZAHL WIRD IN DER EINGABEDATEI UEBER DIE UNIT 5
C****    EINGELESEN. ES IST DIE DRITTE ZAHL DES LETZTEN RECORDS.
C***************************************************************
*
*   English translation of (original) comment:
*
*       Computation of glueball masses in the framework of the
*       SU2 calibration theory.
*       Monte Carlo computation, generation of the configuration
*       with the "warm bath" method.
*       Grid size: 8**3*16, Memory size needed: 5 MByte
*       Number of configurations generated: 20
*       The number is read from the input file via unit 5;
*       it is the third number of the last record
*
****************************************************************

      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXATT=5)
      PARAMETER(ERROR=1.E-8)
      PARAMETER(ERROR4=1.E-7,NDMAX4=8)
      PARAMETER(ERROR2=1.E-7,NDMAX2=10)

      REAL*8 U,ACPT,ENERGY,CD0,CD1,CD2,CD3,RD,SD,BCOEFF,BETINV

      REAL*8 W11SUM(20,3),W11INT(20,3)
      REAL*8 W22SUM(20,3),W22INT(20,3)

      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
      COMMON/PAR3/BETA,GAMMA,ITER,IFREQ,MFREQ
      COMMON/GEOM1/NEWLAT(2**4,-4:4)
      COMMON/GEOM2/ICODE(2**4,-4:4)
      COMMON/GEOM3/INDPBC(MAXVEC,2*4)
      COMMON/SW/ACPT,ENERGY
      COMMON/BESS1/CD0(0:200,2:20),CD1(0:200,2:20)
      COMMON/BESS2/CD2(0:200,2:20),CD3(0:200,2:20)
      COMMON/BESS3/RD(20),SD(20)
      COMMON/BESS4/BCOEFF
      COMMON/INTBLK/BETINV,NERR4,NERR2
      COMMON/CONF/U(MAXVEC,0:3,2**4,4)
      EXTERNAL TRINIT

      WRITE (6,*) 'SPEC benchmark 103.su2cor'

      WRITE (6,*) ' '
C
C     GET START INSTRUCTION
C
*     For historical reasons, we keep the variable ISTART.
*     In the original version, it used to be either -1
*     or the channel number of the channel from which input
*     values are read.

      ISTART = -1

C     GET BASIC PARAMETERS

      READ(5,*)NDIM,(LSIZE(N),N=1,NDIM)

      IF(NDIM.NE.4 .OR. NDIM.GT.4)STOP

      WRITE(6,1000)
1000  FORMAT(1X,75('+')/' PROGRAM SU2V1COR'
     X   ,' -- MODIFIED HEAT BATH ALGORITHM'
     X   /'    -- VECTORIZED PROGRAM IN STANDARD FORTRAN'
     X   /'    -- GLUEBALL CORRELATIONS FOR 1X1 AND 2X2 LOOP OPERATORS'
     X   /'    -- INTEGRAL EXPRESSIONS FOR LOOP OPERATORS USED')
      WRITE(6,1005)NDIM,(LSIZE(N),N=1,NDIM)
1005  FORMAT(//' DIMENSION =',I3
     X        /' SIZE      =',4I3)
C
C     COMPUTE BASIC QUANTITIES
C
      CALL GEOM
C
      WRITE(6,1015)LVEC,MAXATT,200,ERROR,NDMAX4,ERROR4,NDMAX2,ERROR2
1015  FORMAT(//' VECTOR LENGTH: ',I8
     X        /' IN SWEEP:  MAXATT =',I4
     X        /' IN BESTAB: NSTEPS =',I4,'      ERROR  =',E11.3
     X        /' IN INT4:   NDMAX4 =',I4,'      ERROR4 =',E11.3
     X        /' IN INT2:   NDMAX2 =',I4,'      ERROR2 =',E11.3)
      IF(LVEC.GT.MAXVEC) THEN
	WRITE (6,*) 'LVEC=',LVEC
	STOP
	ENDIF
      IF(MOD(LVEC,2).NE.0) THEN
	WRITE (6,*) 'LVEC,mod 2=',LVEC
	STOP
	ENDIF
      IF(LSIZE(NDIM).GT.20) THEN
	WRITE (6,*) 'LSIZE(NDIM)=',LSIZE(NDIM)
	STOP
	ENDIF
      IF(20.LT.NDMAX4 .OR. 20.LT.NDMAX2) THEN
	WRITE (6,*) 'NDMAX.LT.NDMAX4 .OR. NDMAX.LT.NDMAX2'
	STOP
	ENDIF

C     INITIALIZE THE CONFIGURATION

      CALL INIT(ISTART)
      OPEN (11,FILE='SU2COR.MODEL',STATUS='OLD')
C     DO 90 I4=1,4
C     DO 90 I3=1,16
C     DO 90 I2=0,3
C     DO 90 I1=1,LVEC
      READ(11,111) ((((U(I1,I2,I3,I4),I1=1,LVEC),I2=0,3),I3=1,16),
     1I4=1,4)
  90  CONTINUE
 111  FORMAT(D16.10)

C     GET RUN PARAMETERS

200   CONTINUE
      READ(5,*,END=999)BETA,NTHERM,NMEAS,NPRINT,MFREQ,NSAVE,ISAVE
      IF(BETA.LE.0.)GOTO 999
      GAMMA=0.
      WRITE(6,1017)BETA
1017  FORMAT(//1X,75('-')/' BETA  =',F8.4)
C
      CALL BESTAB(BETA,NDIM)
C
      IF(NTHERM.GE.0)THEN
      WRITE(6,1018)NTHERM
1018  FORMAT(/' CLEAR COUNTERS'
     X      //' THERMALISATION:',I6,' ITERATIONS')
      ITER=0
      IFREQ=0
      NERR2=0
      NERR4=0
      CALL CLEAR2(2*NDIM+9+8*(LSIZE(NDIM)/2+1),10)
      CALL CLEAR4(8,2+LSIZE(NDIM)/2)
      ELSE
      WRITE(6,1019)ITER
1019  FORMAT(/' STATISTICS CONTINUED'
     X      //' OLD ITERATIONS:',I6)
      ENDIF
      WRITE(6,1020)NMEAS,MFREQ
1020  FORMAT( ' MEASUREMENT   :',I6,' ITERATIONS'
     X   //' FREQUENCY OF CORRELATION MEASUREMENTS:',I5,' ITERATION(S)')
C
      IF(ISAVE.GT.0)THEN
      WRITE(6,1030)NSAVE,ISAVE
1030  FORMAT(/' SAVE FREQUENCY:',I6,' ITERATIONS, UNIT =',I3)
      ELSE
      WRITE(6,'(/'' CONFIGURATION NOT SAVED'')')
      ENDIF
C
C     THERMALISATION
C
      IF(NTHERM.GT.0)THEN
      IPR=0
      ISA=0
      WRITE(6,1040)(I,I=1,NDIM)
1040  FORMAT(/29X,'E',15X,'FUNDAMENTAL POLYAKOV LINE AVERAGES'
     X     /'   ITER   ACPT      CURRENT     AVERAGE',4(7X,I2,1X))
      DO 50 ITER=1,NTHERM
      CALL SWEEP
C+++++      CALL POLYA
      IPR=IPR+1
      IF(IPR.EQ.NPRINT .OR. ITER.EQ.NTHERM)THEN
      WRITE(6,1050)ITER,ACPT,ENERGY,AVER2(1),(AVER2(I),I=2,2*NDIM,2)
1050  FORMAT(I7,F8.4,2F12.6,2X,5F10.6)
      IPR=0
      ENDIF
      ISA=ISA+1
      IF(ISA.GE.NSAVE .AND. ISAVE.GT.0)THEN
C     CALL SAVE(ISAVE)
      ISA=0
      ENDIF
50    CONTINUE
      ITER=0
      IFREQ=0
      NERR2=0
      NERR4=0
      CALL CLEAR2(2*NDIM+9+8*(LSIZE(NDIM)/2+1),10)
      CALL CLEAR4(8,2+LSIZE(NDIM)/2)
C     IF(ISAVE.GT.0) CALL SAVE(ISAVE)
      ENDIF
C
C     MEASUREMENT
C
      IPR=0
      ISA=0
      WRITE(6,1040)(I,I=1,NDIM)
      DO 60 IT=1,NMEAS
      ITER=ITER+1
      CALL SWEEP
C+++++      CALL POLYA
      IPR=IPR+1
      IF(IPR.EQ.NPRINT .OR. IT.EQ.NMEAS)THEN
      WRITE(6,1050)ITER,ACPT,ENERGY,AVER2(1),(AVER2(I),I=2,2*NDIM,2)
      IPR=0
      ENDIF
      IFREQ=IFREQ+1
      IF(IFREQ.GE.MFREQ)THEN
      CALL CORR
      IFREQ=0
      ENDIF
      ISA=ISA+1
      IF((ISA.GE.NSAVE .OR. IT.EQ.NMEAS) .AND. ISAVE.GT.0)THEN
C     CALL SAVE(ISAVE)
      ISA=0
      CALL LOOPS(W11SUM,W11INT,W22SUM,W22INT)
      SUG = 0.
      DO 70 IPL = 1, 3
      SU = 0.
      DO 71 ITIT = 1, LSIZE(NDIM)
          SU = SU + W11SUM(ITIT,IPL)
   71 CONTINUE
      SU = SU / REAL(LSIZE(NDIM))
      SUG = SUG + SU
   70 CONTINUE
      SUG = SUG/3
      ENDIF
60    CONTINUE
C
C     RESULTS AND ERROR ANALYSIS
C
      IF(ITER.GT.0) CALL EVAL
C
      GOTO 200
C PROGRAM NEEDS THE LABEL
999   CONTINUE
      STOP
      END
C*************************EVAL*****************************
C
C     PRINT RESULTS AND ERROR ANALYSIS
C
      SUBROUTINE EVAL
*     ***************
C
      PARAMETER(MAXVEC=8**3*8+7)
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
C
      REAL*8 ED,PF,PA
      REAL*8 S110,S112,S110I,S112I,S220,S222,S220I,S222I
      REAL*8 C110,C112,C110I,C112I,C220,C222,C220I,C222I
      REAL*8 COVAR4
C
      CALL STAT2(1,ED,SIGE,RE,0)
C
*      WRITE(6,1080)
*1080  FORMAT(/' POLYAKOV LINE AVERAGES:'
*     X   /'  I',12X,'FUNDAMENTAL       CORR.FAC.'
*     X   ,'            ADJOINT         CORR.FAC.')
      DO 70 I=1,NDIM
      CALL STAT2(2*I,PF,SIGPF,RPF,0)
      CALL STAT2(2*I+1,PA,SIGPA,RPA,0)
70    CONTINUE
C
*     X      /11X,'W11(0+)',25X,'W22(0+)'
*     X      ,25X,'W11(2+)',25X,'W22(2+)')
      CALL STAT2(10,S110,ES110,RS110,0)
      CALL STAT2(11,S112,ES112,RS112,0)
      CALL STAT2(12,S110I,ES110I,RS110I,0)
      CALL STAT2(13,S112I,ES112I,RS112I,0)
      CALL STAT2(14,S220,ES220,RS220,0)
      CALL STAT2(15,S222,ES222,RS222,0)
      CALL STAT2(16,S220I,ES220I,RS220I,0)
      CALL STAT2(17,S222I,ES222I,RS222I,0)
      MAXDST=LSIZE(NDIM)/2
      NDIST=MAXDST+1
      DO 80 IDIST=0,MAXDST
      IND=IDIST+1
      IND1=IND+1
C
      CALL STAT2(17+IND,C110,EC110,RC110,0)
      C110=(C110-S110**2) /NVOL
      EC110=COVAR4(1,IND1,IND1)
     X      + 4.*S110*( S110*COVAR4(1,1,1) - COVAR4(1,1,IND1) )
      EC110=SQRT(MAX(0.,EC110))*RC110 /NVOL
C
      CALL STAT2(17+NDIST+IND,C112,EC112,RC112,0)
      C112=C112 /NVOL
      EC112=EC112 /NVOL
C
      CALL STAT2(17+4*NDIST+IND,C220,EC220,RC220,0)
      C220=(C220-S220**2) /NVOL
      EC220=COVAR4(5,IND1,IND1)
     X      + 4.*S220*( S220*COVAR4(5,1,1) - COVAR4(5,1,IND1) )
      EC220=SQRT(MAX(0.,EC220))*RC220 /NVOL
C
      CALL STAT2(17+5*NDIST+IND,C222,EC222,RC222,0)
      C222=C222 /NVOL
      EC222=EC222 /NVOL
C
      CALL STAT2(17+2*NDIST+IND,C110I,EC110I,RC110I,0)
      C110I=(C110I-S110I**2) /NVOL
      EC110I=COVAR4(3,IND1,IND1)
     X       + 4.*S110I*( S110I*COVAR4(3,1,1) - COVAR4(3,1,IND1) )
      EC110I=SQRT(MAX(0.,EC110I))*RC110I /NVOL
C
      CALL STAT2(17+3*NDIST+IND,C112I,EC112I,RC112I,0)
      C112I=C112I /NVOL
      EC112I=EC112I /NVOL
C
      CALL STAT2(17+6*NDIST+IND,C220I,EC220I,RC220I,0)
      C220I=(C220I-S220I**2) /NVOL
      EC220I=COVAR4(7,IND1,IND1)
     X       + 4.*S220I*( S220I*COVAR4(7,1,1) - COVAR4(7,1,IND1) )
      EC220I=SQRT(MAX(0.,EC220I))*RC220I /NVOL
C
      CALL STAT2(17+7*NDIST+IND,C222I,EC222I,RC222I,0)
      C222I=C222I /NVOL
      EC222I=EC222I /NVOL
C
80    CONTINUE
C
      END
C*************************SWEEP**************************
C
C     ONE VECTORISED SWEEP USING A MODIFIED HEAT BATH ALGORITHM
C
      SUBROUTINE SWEEP
*     ****************
C
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXATT=5)
C
      REAL*8 U,ACPT,ENERGY
      REAL*8 BETINV,PI,ZERO,ONE,TWO,EPS
      REAL*8 B,W1,W2,BNORM,BINV,WT,W0,C1,C2,R1,R2
      PARAMETER(PI=3.141592654D0,ZERO=0.D0,ONE=1.D0,TWO=2.D0,EPS=1.E-20)
C
      REAL*8 EDOUB
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
      COMMON/PAR3/BETA,GAMMA,ITER,IFREQ,MFREQ
      COMMON/SW/ACPT,ENERGY
      COMMON/CONF/U(MAXVEC,0:3,2**4,4)
C
      DIMENSION B(MAXVEC,0:3),W1(MAXVEC,0:3),W2(MAXVEC,0:3)
      DIMENSION BNORM(MAXVEC),BINV(MAXVEC),WT(MAXVEC),W0(MAXVEC)
      DIMENSION C1(MAXVEC),C2(MAXVEC),R1(MAXVEC),R2(MAXVEC)
C+++++EQUIVALENCE (BNORM(1),W2(1,0)),(BINV(1),W2(1,1))
C+++++EQUIVALENCE (WT(1),W2(1,2)),(W0(1),W2(1,3))
      LOGICAL LACPT(MAXVEC)
C
      EDOUB=0.D0
      ENERGY=ZERO
      IACPT=0
      BETINV=1.D0/BETA
C
C     LOOP OVER SUB-LATTICES  ILAT  AND LINK DIRECTIONS  I1
C
      DO 200 ILAT=1,2**NDIM
      DO 200 I1=1,NDIM
C
C     ACCUMULATE THE INTERACTING PLAQUETTES IN  B(...)
C
      DO 210 I=1,LVEC
      B(I,0)=ZERO
      B(I,1)=ZERO
      B(I,2)=ZERO
      B(I,3)=ZERO
      LACPT(I)=.FALSE.
210   CONTINUE
C
      DO 220 I2=1,NDIM
      IF(I1.EQ.I2)GOTO 220
      CALL INTACT(ILAT,I1, I2,W1)
      CALL INTACT(ILAT,I1,-I2,W2)
      DO 230 I=1,LVEC
      B(I,0)=B(I,0)+W1(I,0)+W2(I,0)
      B(I,1)=B(I,1)+W1(I,1)+W2(I,1)
      B(I,2)=B(I,2)+W1(I,2)+W2(I,2)
      B(I,3)=B(I,3)+W1(I,3)+W2(I,3)
230   CONTINUE
220   CONTINUE
C
C     HEAT BATH STEP
C
      DO 300 I=1,LVEC
      BNORM(I)=SQRT( B(I,0)**2 + B(I,1)**2 + B(I,2)**2 + B(I,3)**2 )
      BINV(I)=ONE/BNORM(I)
      C1(I)=BETINV*BINV(I)
300   CONTINUE
C
      DO 310 IATT=1,MAXATT
      CALL TRNGV(R1,LVEC,ONE)
      CALL TRNGV(R2,LVEC,ONE)
      DO 311 I=1,LVEC
      WT(I)=ONE + C1(I)*LOG(EPS+R1(I))
      IF( R2(I)**2 .LE. (ONE-WT(I)**2) )THEN
      W0(I)=WT(I)
      LACPT(I)=.TRUE.
      ENDIF
311   CONTINUE
310   CONTINUE
C
      CALL TRNGV(R1,LVEC,TWO)
      CALL TRNGV(R2,LVEC,PI)
C
      DO 320 I=1,LVEC
      IF(LACPT(I))THEN
C
      IACPT=IACPT+1
      EDOUB=EDOUB + BNORM(I)*W0(I)
      W1(I,0)=BINV(I) * W0(I)
C
      BINV(I)=BINV(I) * SQRT(ONE-W0(I)**2)
      R1(I)=R1(I) - ONE
      W1(I,3)=BINV(I) * R1(I)
C
      C1(I)=TAN(R2(I))
      C2(I)=C1(I)**2
      BINV(I)=BINV(I) * SQRT(ONE-R1(I)**2) / (ONE+C2(I))
      W1(I,1)=BINV(I) * (ONE-C2(I))
      W1(I,2)=BINV(I) * TWO*C1(I)
C
      ELSE
C
      W1(I,0)= U(I,0,ILAT,I1)
      W1(I,1)=-U(I,1,ILAT,I1)
      W1(I,2)=-U(I,2,ILAT,I1)
      W1(I,3)=-U(I,3,ILAT,I1)
C
      ENDIF
320   CONTINUE
C
C     COMPUTE AND STORE THE NEW LINK VARIABLES
C
      CALL MATMAT(W1,B,W2,LVEC,MAXVEC)
C
      DO 400 I=1,LVEC
      IF(LACPT(I))THEN
      U(I,0,ILAT,I1)=W2(I,0)
      U(I,1,ILAT,I1)=W2(I,1)
      U(I,2,ILAT,I1)=W2(I,2)
      U(I,3,ILAT,I1)=W2(I,3)
      ELSE
      EDOUB=EDOUB + W2(I,0)
      ENDIF
400   CONTINUE
C
200   CONTINUE
C
C     END OF LOOP OVER SUB-LATTICES AND LINK DIRECTIONS
C
      ENERGY=EDOUB/(4*NPLAQS)
      CALL ACCUM2(1,DBLE(ENERGY))
      ACPT=FLOAT(IACPT)/NLINKS
      END
C*************************GEOM****************************
C
C     COMPUTE GEOMETRY VECTORS
C
      SUBROUTINE GEOM
*     ***************
C
      PARAMETER(MAXVEC=8**3*8+7)
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
      COMMON/GEOM1/NEWLAT(2**4,-4:4)
      COMMON/GEOM2/ICODE(2**4,-4:4)
      COMMON/GEOM3/INDPBC(MAXVEC,2*4)
C
      DIMENSION LHALF(4),IBASE(4),JUMP(4)
      DIMENSION IX(4),IUP(4),IDN(4),ISTEP(4)
      LOGICAL LZ(4)
C
C     COMPUTE SOME BASIC PARAMETERS
C
      NDSP=NDIM-1
C
      NPTS=1
      LVEC=1
      DO 10 I=1,NDIM
      IF(MOD(LSIZE(I),2) .NE. 0) STOP
      NVOL=NPTS
      NPTS=NPTS*LSIZE(I)
      LHALF(I)=LSIZE(I)/2
      LVEC=LVEC*LHALF(I)
10    CONTINUE
C
      NLINKS=NPTS*NDIM
      NPLPP=(NDIM*(NDIM-1))/2
      NPLAQS=NPTS*NPLPP
C
C     COMPUTE THE VECTORS  NEWLAT(..)  AND  ICODE(..)
C
C        NEWLAT(ILAT,I1)  DETERMINES THE INDEX OF THE SUB-LATTICE
C              REACHED FROM  ILAT  IN DIRECTION  I1
C        ICODE(ILAT,I1)  CONTAINS THE INDEX OF THE PERMUTATION VECTOR
C              INDPBC  NEEDED ON SUB-LATTICE  ILAT  FOR A STEP  I1
C              ICODE=0 IF NO PERMUTATION IS NEEDED
C        I1 = +/-1,...,+/-NDIM
C
      DO 20 I=1,NDIM
      ISTEP(I)=2**(I-1)
      LZ(I)=.FALSE.
20    CONTINUE
C
      ILAT=0
100   CONTINUE
      ILAT=ILAT+1
C
      NEWLAT(ILAT,0)=-1
      ICODE(ILAT,0)=-1
C
      DO 110 I1=1,NDIM
C
      NEWLAT(ILAT, I1)=ILAT+ISTEP(I1)
      NEWLAT(ILAT,-I1)=ILAT+ISTEP(I1)
C
      IF(LZ(I1))THEN
      ICODE(ILAT, I1)=I1
      ICODE(ILAT,-I1)=0
      ELSE
      ICODE(ILAT, I1)=0
      ICODE(ILAT,-I1)=NDIM+I1
      ENDIF
C
110   CONTINUE
C
      DO 190 I=1,NDIM
      LZ(I)=.NOT.LZ(I)
      ISTEP(I)=-ISTEP(I)
      IF(LZ(I))GOTO 100
190   CONTINUE
C
C     COMPUTE THE INDEX VECTORS  INDPBC(...)  FOR PERMUTATIONS OF
C              VECTORS SUBJECT TO PERIODIC BOUNDARY CONDITIONS
C         INDPBC(IPT,IND)  CONTAINS THE LOCATION, WITHIN ITS VECTOR,
C              OF THE NEIGHBOUR OF THE (BLOCKED) POINT  IPT .
C              THE SHIFT  +/-I1  IS CODED IN  IND  AS DEFINED
C              DURING THE CONSTRUCTION OF THE ARRAY  ICODE(...) .
C
      IBASE(1)=1
      DO 210 I=2,NDIM
      IBASE(I)=IBASE(I-1)*LHALF(I-1)
      JUMP(I-1)=IBASE(I)-IBASE(I-1)
210   CONTINUE
      JUMP(NDIM)=IBASE(NDIM)*(LHALF(NDIM)-1)
C
      DO 220 I=1,NDIM
      IUP(I)=IBASE(I)
      IDN(I)=JUMP(I)
      IX(I)=1
220   CONTINUE
      IPT=0
C
300   CONTINUE
      IPT=IPT+1
C
      DO 310 I1=1,NDIM
      INDPBC(IPT,I1)=IPT+IUP(I1)
      INDPBC(IPT,NDIM+I1)=IPT+IDN(I1)
310   CONTINUE
C
      DO 390 I=1,NDIM
      IX(I)=IX(I)+1
      IDN(I)=-IBASE(I)
      IF(IX(I).LT.LHALF(I))GOTO 300
      IF(IX(I).EQ.LHALF(I))THEN
      IUP(I)=-JUMP(I)
      GOTO 300
      ENDIF
      IX(I)=1
      IUP(I)=IBASE(I)
      IDN(I)=JUMP(I)
390   CONTINUE
C
      END
C****************MATRIX MULTIPLICATIONS*********************
C
C     VECTORIZED MULTIPLICATION OF SU(2) MATRICES
C
C          LVEC: VECTOR LENGTH
C          LDIM: FIRST DIMENSIONAL PARAMETER OF U,V,W
C                    IN THE CALLING PROGRAM
C
      SUBROUTINE MATMAT(U,V,W,LVEC,LDIM)
*     *****************
C
      REAL*8 U,V,W
C
      DIMENSION U(LDIM,0:3),V(LDIM,0:3),W(LDIM,0:3)
C
C     W = U * V
C
      DO 10 I=1,LVEC
      W(I,0)=U(I,0)*V(I,0)-U(I,1)*V(I,1)-U(I,2)*V(I,2)-U(I,3)*V(I,3)
      W(I,1)=U(I,0)*V(I,1)+U(I,1)*V(I,0)+U(I,2)*V(I,3)-U(I,3)*V(I,2)
      W(I,2)=U(I,0)*V(I,2)-U(I,1)*V(I,3)+U(I,2)*V(I,0)+U(I,3)*V(I,1)
      W(I,3)=U(I,0)*V(I,3)+U(I,1)*V(I,2)-U(I,2)*V(I,1)+U(I,3)*V(I,0)
10    CONTINUE
      RETURN
C
      ENTRY ADDMM(U,V,W,LVEC,LDIM)
C
C     W = W + U * V
C
      DO 20 I=1,LVEC
      W(I,0)=W(I,0)
     X       +U(I,0)*V(I,0)-U(I,1)*V(I,1)-U(I,2)*V(I,2)-U(I,3)*V(I,3)
      W(I,1)=W(I,1)
     X       +U(I,0)*V(I,1)+U(I,1)*V(I,0)+U(I,2)*V(I,3)-U(I,3)*V(I,2)
      W(I,2)=W(I,2)
     X       +U(I,0)*V(I,2)-U(I,1)*V(I,3)+U(I,2)*V(I,0)+U(I,3)*V(I,1)
      W(I,3)=W(I,3)
     X       +U(I,0)*V(I,3)+U(I,1)*V(I,2)-U(I,2)*V(I,1)+U(I,3)*V(I,0)
20    CONTINUE
      RETURN
C
      ENTRY ADJMAT(U,V,W,LVEC,LDIM)
C
C     W = ADJ(U) * V
C
      DO 30 I=1,LVEC
      W(I,0)=U(I,0)*V(I,0)+U(I,1)*V(I,1)+U(I,2)*V(I,2)+U(I,3)*V(I,3)
      W(I,1)=U(I,0)*V(I,1)-U(I,1)*V(I,0)-U(I,2)*V(I,3)+U(I,3)*V(I,2)
      W(I,2)=U(I,0)*V(I,2)+U(I,1)*V(I,3)-U(I,2)*V(I,0)-U(I,3)*V(I,1)
      W(I,3)=U(I,0)*V(I,3)-U(I,1)*V(I,2)+U(I,2)*V(I,1)-U(I,3)*V(I,0)
30    CONTINUE
      RETURN
C
      ENTRY ADDAM(U,V,W,LVEC,LDIM)
C
C     W = W + ADJ(U) * V
C
      DO 40 I=1,LVEC
      W(I,0)=W(I,0)
     X       +U(I,0)*V(I,0)+U(I,1)*V(I,1)+U(I,2)*V(I,2)+U(I,3)*V(I,3)
      W(I,1)=W(I,1)
     X       +U(I,0)*V(I,1)-U(I,1)*V(I,0)-U(I,2)*V(I,3)+U(I,3)*V(I,2)
      W(I,2)=W(I,2)
     X       +U(I,0)*V(I,2)+U(I,1)*V(I,3)-U(I,2)*V(I,0)-U(I,3)*V(I,1)
      W(I,3)=W(I,3)
     X       +U(I,0)*V(I,3)-U(I,1)*V(I,2)+U(I,2)*V(I,1)-U(I,3)*V(I,0)
40    CONTINUE
      RETURN
C
      ENTRY MATADJ(U,V,W,LVEC,LDIM)
C
C     W = U * ADJ(V)
C
      DO 50 I=1,LVEC
      W(I,0)=+U(I,0)*V(I,0)+U(I,1)*V(I,1)+U(I,2)*V(I,2)+U(I,3)*V(I,3)
      W(I,1)=-U(I,0)*V(I,1)+U(I,1)*V(I,0)-U(I,2)*V(I,3)+U(I,3)*V(I,2)
      W(I,2)=-U(I,0)*V(I,2)+U(I,1)*V(I,3)+U(I,2)*V(I,0)-U(I,3)*V(I,1)
      W(I,3)=-U(I,0)*V(I,3)-U(I,1)*V(I,2)+U(I,2)*V(I,1)+U(I,3)*V(I,0)
50    CONTINUE
      RETURN
C
      ENTRY ADDMA(U,V,W,LVEC,LDIM)
C
C     W = W + U * ADJ(V)
C
      DO 60 I=1,LVEC
      W(I,0)=W(I,0)
     X       +U(I,0)*V(I,0)+U(I,1)*V(I,1)+U(I,2)*V(I,2)+U(I,3)*V(I,3)
      W(I,1)=W(I,1)
     X       -U(I,0)*V(I,1)+U(I,1)*V(I,0)-U(I,2)*V(I,3)+U(I,3)*V(I,2)
      W(I,2)=W(I,2)
     X       -U(I,0)*V(I,2)+U(I,1)*V(I,3)+U(I,2)*V(I,0)-U(I,3)*V(I,1)
      W(I,3)=W(I,3)
     X       -U(I,0)*V(I,3)-U(I,1)*V(I,2)+U(I,2)*V(I,1)+U(I,3)*V(I,0)
60    CONTINUE
      RETURN
C
      ENTRY ADJADJ(U,V,W,LVEC,LDIM)
C
C     W = ADJ(U) * ADJ(V)
C
      DO 70 I=1,LVEC
      W(I,0)=+U(I,0)*V(I,0)-U(I,1)*V(I,1)-U(I,2)*V(I,2)-U(I,3)*V(I,3)
      W(I,1)=-U(I,0)*V(I,1)-U(I,1)*V(I,0)+U(I,2)*V(I,3)-U(I,3)*V(I,2)
      W(I,2)=-U(I,0)*V(I,2)-U(I,1)*V(I,3)-U(I,2)*V(I,0)+U(I,3)*V(I,1)
      W(I,3)=-U(I,0)*V(I,3)+U(I,1)*V(I,2)-U(I,2)*V(I,1)-U(I,3)*V(I,0)
70    CONTINUE
      RETURN
C
      ENTRY ADDAA(U,V,W,LVEC,LDIM)
C
C     W = W + ADJ(U) * ADJ(V)
C
      DO 80 I=1,LVEC
      W(I,0)=W(I,0)
     X       +U(I,0)*V(I,0)-U(I,1)*V(I,1)-U(I,2)*V(I,2)-U(I,3)*V(I,3)
      W(I,1)=W(I,1)
     X       -U(I,0)*V(I,1)-U(I,1)*V(I,0)+U(I,2)*V(I,3)-U(I,3)*V(I,2)
      W(I,2)=W(I,2)
     X       -U(I,0)*V(I,2)-U(I,1)*V(I,3)-U(I,2)*V(I,0)+U(I,3)*V(I,1)
      W(I,3)=W(I,3)
     X       -U(I,0)*V(I,3)+U(I,1)*V(I,2)-U(I,2)*V(I,1)-U(I,3)*V(I,0)
80    CONTINUE
      END
C*************************PERM**************************************
C
C     PERMUTE A VECTOR OF SU(2) MATRICES  U(..,..)  ACCORDING TO THE
C        INDEX VECTOR  INDPBC(..,IC)   AND RETURN THE RESULT IN
C        V(..,..)
C
C        LVEC: VECTOR LENGTH
C        LDIM: FIRST DIMENSIONAL PARAMETER OF U,V,W
C                 IN THE CALLING PROGRAM
C
      SUBROUTINE PERM(U,IC,V,LVEC)
*     ***************
C
      PARAMETER(MAXVEC=8**3*8+7)
C
      REAL*8 U,V
C
      COMMON/GEOM3/INDPBC(MAXVEC,2*4)
C
      DIMENSION U(MAXVEC,0:3),V(MAXVEC,0:3)
C
      DO 10 I=2,LVEC,2
      V(I-1,0)=U(INDPBC(I-1,IC),0)
      V(I-1,1)=U(INDPBC(I-1,IC),1)
      V(I-1,2)=U(INDPBC(I-1,IC),2)
      V(I-1,3)=U(INDPBC(I-1,IC),3)
      V(I,0)=U(INDPBC(I,IC),0)
      V(I,1)=U(INDPBC(I,IC),1)
      V(I,2)=U(INDPBC(I,IC),2)
      V(I,3)=U(INDPBC(I,IC),3)
10    CONTINUE
      END
C*************************INTACT************************************
C
C     COMPUTE THE (VECTORISED) PRODUCT OF THREE INTERACTING NEIGHBOURS
C          OF A LINK  (ILAT2,I1)  AROUND THE PLAQUETTE  (I1,I2)
C       I1=1,..,NDIM
C       I2=+/-1,..,+/-NDIM
C
      SUBROUTINE INTACT(ILAT2,I1,I2,PROD)
*     *****************
C
      PARAMETER(MAXVEC=8**3*8+7)
C
      REAL*8 U
      REAL*8 PROD,W1,W2
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
      COMMON/PAR3/BETA,GAMMA,ITER,IFREQ,MFREQ
      COMMON/GEOM1/NEWLAT(2**4,-4:4)
      COMMON/GEOM2/ICODE(2**4,-4:4)
      COMMON/CONF/U(MAXVEC,0:3,2**4,4)
C
      DIMENSION PROD(MAXVEC,0:3),W1(MAXVEC,0:3),W2(MAXVEC,0:3)
C
C     THE SUB-LATTICE POINTS ARE LABELLED AS FOLLOWS:
C
C               I1
C                .
C                .    4   5   6
C                .    1   2   3
C                . . . . . . . . . . ABS(I2)
C
C     ILAT2  REFERS TO POINT '2'
C     '2-5' IS THE CENTRAL LINK  (ILAT2,I1)
C
C     THE PRODUCT  '2-3-6-5'  (FOR POSITIVE  I2)
C              OR  '2-1-4-5'  (FOR NEGATIVE  I2)
C        IS RETURNED IN  PROD(...)
C
      IF(I2.GT.0)THEN
C
C     CASE  I1 POSITIVE , I2 POSITIVE
C
      ILAT3=NEWLAT(ILAT2,I2)
      ILAT5=NEWLAT(ILAT2,I1)
      IC23=ICODE(ILAT2,I2)
      IC25=ICODE(ILAT2,I1)
C
      IF(IC23.EQ.0)THEN
      CALL MATMAT(U(1,0,ILAT3,I1),U(1,0,ILAT2,I2),W1,
     X                        LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,ILAT3,I1),IC23,W2,LVEC)
      CALL MATMAT(W2,U(1,0,ILAT2,I2),W1,LVEC,MAXVEC)
      ENDIF
C
      IF(IC25.EQ.0)THEN
      CALL ADJMAT(U(1,0,ILAT5,I2),W1,PROD,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,ILAT5,I2),IC25,W2,LVEC)
      CALL ADJMAT(W2,W1,PROD,LVEC,MAXVEC)
      ENDIF
C
      ELSE
C
C     CASE  I1 POSITIVE , I2 NEGATIVE
C
      ILAT1=NEWLAT(ILAT2,I2)
      ILAT4=NEWLAT(ILAT1,I1)
      IC21=ICODE(ILAT2,I2)
      IC14=ICODE(ILAT1,I1)
C
      CALL MATADJ(U(1,0,ILAT1,I1),U(1,0,ILAT1,-I2),W1,
     X                    LVEC,MAXVEC)
      IF(IC14.EQ.0)THEN
      IF(IC21.EQ.0)THEN
      CALL MATMAT(U(1,0,ILAT4,-I2),W1,PROD,LVEC,MAXVEC)
      ELSE
      CALL MATMAT(U(1,0,ILAT4,-I2),W1,W2,LVEC,MAXVEC)
      CALL PERM(W2,IC21,PROD,LVEC)
      ENDIF
      ELSE
      IF(IC21.EQ.0)THEN
      CALL PERM(U(1,0,ILAT4,-I2),IC14,W2,LVEC)
      CALL MATMAT(W2,W1,PROD,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,ILAT4,-I2),IC14,PROD,LVEC)
      CALL MATMAT(PROD,W1,W2,LVEC,MAXVEC)
      CALL PERM(W2,IC21,PROD,LVEC)
      ENDIF
      ENDIF
C
      ENDIF
C
      END
C*************************BESTABNN*****************************
C
      SUBROUTINE BESTAB(BETA,NDIM)
*     *****************
C
      PARAMETER(ERROR=1.E-8)
C
      REAL*8 CD0,CD1,CD2,CD3,RD,SD,BCOEFF
C
      COMMON/BESS1/CD0(0:200,2:20),CD1(0:200,2:20)
      COMMON/BESS2/CD2(0:200,2:20),CD3(0:200,2:20)
      COMMON/BESS3/RD(20),SD(20)
      COMMON/BESS4/BCOEFF
C
      DIMENSION CD(-1:200+2,2:20)
C
      BCOEFF=FLOAT(200)/(2*(NDIM-1))**2
C
      STEP=((NDIM-1)*BETA)**2 /200
      BHALF=.5D0*BETA
C
C     LOOP OVER POINTS
C
      DO 100 IBSQ=-1,200+2
      U=IBSQ*STEP
C
C     COMPUTE THE NUMERATOR AND DENOMINATOR SUMS FOR  CD(BETA,B)
C
      DENOM=1.D0
      CD(IBSQ,2)=.5D0
      DO 10 ND=3,20
      CD(IBSQ,ND)=CD(IBSQ,ND-1) /ND
10    CONTINUE
      S=1.
      K=0
C
20    CONTINUE
      K=K+1
      S=S*U/(K*(K+1))
      DENOM=DENOM+S
      S1=S
      DO 30 ND=2,20
      S1=S1/(K+ND)
      CD(IBSQ,ND)=CD(IBSQ,ND)+S1
30    CONTINUE
      IF(ABS(S).GT.DENOM*ERROR)GOTO 20
C
C     COMPUTE  CD(BETA,B)
C
      F=1./DENOM
      DO 40 ND=2,20
      F=F*BHALF
      CD(IBSQ,ND)=F*CD(IBSQ,ND)
40    CONTINUE
C
100   CONTINUE
C
C     COMPUTE THE COEFFICIENTS OF THE INTERPOLATING POLYNOMIALS
C
      DO 200 IBSQ=0,200
      DO 200 ND=2,20
      CD0(IBSQ,ND)=CD(IBSQ,ND)
      CD1(IBSQ,ND)=-CD(IBSQ-1,ND)/3.D0 - .5D0*CD(IBSQ,ND)
     X                    + CD(IBSQ+1,ND) - CD(IBSQ+2,ND)/6.D0
      CD2(IBSQ,ND)=.5D0*( CD(IBSQ-1,ND) + CD(IBSQ+1,ND) ) - CD(IBSQ,ND)
      CD3(IBSQ,ND)=( CD(IBSQ+2,ND) - CD(IBSQ-1,ND) )/6.D0
     X                   + .5D0*( CD(IBSQ,ND) - CD(IBSQ+1,ND) )
200   CONTINUE
C
C     COMPUTE THE SUMS FOR  RD(BETA)
C
      RD(1)=1.D0
      DO 50 ND=2,20
      RD(ND)=RD(ND-1)/ND
50    CONTINUE
      U=BHALF**2
      S=1.
      K=0
C
60    CONTINUE
      K=K+1
      S=S*U/K**2
      S1=S
      DO 70 ND=1,20
      S1=S1/(K+ND)
      RD(ND)=RD(ND)+S1
70    CONTINUE
      IF(S.GT.RD(1)*ERROR)GOTO 60
C
C     COMPUTE  RD(BETA)  AND  SD(BETA)
C
      F=1./RD(1)
      RD(1)=1.
      DO 80 ND=2,20
      D=ND
      F=F*BHALF
      RD(ND)=D*F*RD(ND)
      SD(ND)=D/(D-1.D0)*RD(ND-1) - (D+1.D0)/BETA*RD(ND)
80    CONTINUE
      SD(1)=.5D0*RD(2)
      END
C*************************BESPOL****************************
C
C     COMPUTE THE INTERPOLATING POLYNOMIAL FOR TABULATED
C          RATIOS OF BESSEL FUNCTIONS
C
      SUBROUTINE BESPOL(ND,IND,DIFF,POL,LVEC)
*     *****************
C
      PARAMETER(ERROR=1.E-8)
C
      REAL*8 CD0,CD1,CD2,CD3
      REAL*8 DIFF,POL
C
      COMMON/BESS1/CD0(0:200,2:20),CD1(0:200,2:20)
      COMMON/BESS2/CD2(0:200,2:20),CD3(0:200,2:20)
C
      DIMENSION IND(LVEC),DIFF(LVEC),POL(LVEC)
C
      DO 10 I=2,LVEC,2
      POL(I-1)=
     X     ( CD0(IND(I-1),ND)+DIFF(I-1)*( CD1(IND(I-1),ND)+DIFF(I-1)*
     X     ( CD2(IND(I-1),ND)+DIFF(I-1)*  CD3(IND(I-1),ND) )))
      POL(I)=( CD0(IND(I),ND)+DIFF(I)*( CD1(IND(I),ND)+DIFF(I)*
     X       ( CD2(IND(I),ND)+DIFF(I)*  CD3(IND(I),ND) )))
10    CONTINUE
      END
C*************************INT4V*****************************
C
C     COMPUTE THE 4-LINK INTEGRAL EXPRESSIONS FOR THE
C          PLAQUETTE AVERAGE (VECTORISED VERSION)
C
      SUBROUTINE INT4V(B1,B2,B3,B4,PLINT)
*     ****************
C
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(ERROR=1.E-8)
      PARAMETER(ERROR4=1.E-7,NDMAX4=8)
C
      REAL*8 BCOEFF,BETINV,RD,SD
      REAL*8 B1,B2,B3,B4,PLINT,W1,W2,TRACE,BSQ,D1,D2,D3,D4
      REAL*8 SUM1,SUM2,CHIDM1,CHID,BES1,BES2,BES3,BES4,F6
      REAL*8 B1SQ,B2SQ,B3SQ,B4SQ,CHIOLD,F5,ERR4,ZERO,ONE,TWO
      PARAMETER(ERR4=ERROR4,ZERO=0.D0,ONE=1.D0,TWO=2.D0)
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/INTBLK/BETINV,NERR4,NERR2
      COMMON/BESS3/RD(20),SD(20)
      COMMON/BESS4/BCOEFF
C
      DIMENSION B1(MAXVEC,0:3),B2(MAXVEC,0:3)
      DIMENSION B3(MAXVEC,0:3),B4(MAXVEC,0:3),PLINT(MAXVEC)
C
      DIMENSION W1(MAXVEC,0:3),W2(MAXVEC,0:3),TRACE(MAXVEC),BSQ(MAXVEC)
      DIMENSION D1(MAXVEC),D2(MAXVEC),D3(MAXVEC),D4(MAXVEC)
      DIMENSION I1(MAXVEC),I2(MAXVEC),I3(MAXVEC),I4(MAXVEC)
      DIMENSION SUM1(MAXVEC),SUM2(MAXVEC),CHIDM1(MAXVEC),CHID(MAXVEC)
      DIMENSION BES1(MAXVEC),BES2(MAXVEC),BES3(MAXVEC),BES4(MAXVEC)
      DIMENSION F6(MAXVEC)
C
C+++++EQUIVALENCE (BES1,W1(1,0)),(BES2,W1(1,1))
C+++++EQUIVALENCE (BES3,W1(1,2)),(BES4,W1(1,3))
C+++++EQUIVALENCE (D1,W2(1,0)),(D2,W2(1,1))
C+++++EQUIVALENCE (D3,W2(1,2)),(D4,W2(1,3))
C
C     PRODUCT B2*B1
C
      CALL MATMAT(B2,B1,W1,LVEC,MAXVEC)
C
C     PRODUCT B3*B4
C
      CALL MATMAT(B3,B4,W2,LVEC,MAXVEC)
C
      DO 50 I=1,LVEC
C
      TRACE(I)=TWO*( W1(I,0)*W2(I,0) + W1(I,1)*W2(I,1)
     X             + W1(I,2)*W2(I,2) + W1(I,3)*W2(I,3) )
C
      B1SQ=B1(I,0)**2 + B1(I,1)**2 + B1(I,2)**2 + B1(I,3)**2
      B2SQ=B2(I,0)**2 + B2(I,1)**2 + B2(I,2)**2 + B2(I,3)**2
      B3SQ=B3(I,0)**2 + B3(I,1)**2 + B3(I,2)**2 + B3(I,3)**2
      B4SQ=B4(I,0)**2 + B4(I,1)**2 + B4(I,2)**2 + B4(I,3)**2
C
      BSQ(I)=B1SQ*B2SQ*B3SQ*B4SQ
C
      D1(I)=BCOEFF*B1SQ
      I1(I)=D1(I)
      D1(I)=D1(I)-I1(I)
C
      D2(I)=BCOEFF*B2SQ
      I2(I)=D2(I)
      D2(I)=D2(I)-I2(I)
C
      D3(I)=BCOEFF*B3SQ
      I3(I)=D3(I)
      D3(I)=D3(I)-I3(I)
C
      D4(I)=BCOEFF*B4SQ
      I4(I)=D4(I)
      D4(I)=D4(I)-I4(I)
C
C     INITIALISE THE SUMS AND RECURSION RELATIONS FOR D=1
C
      SUM1(I)=SD(1)
      SUM2(I)=ONE
C
      CHIDM1(I)=ZERO
      CHID(I)=ONE
C
50    CONTINUE
C
C     LOOP OVER D=2,3,4,...,NDMAX4
C
      DO 100 ND=2,NDMAX4
C
      CALL BESPOL(ND,I1,D1,BES1,LVEC)
      CALL BESPOL(ND,I2,D2,BES2,LVEC)
      CALL BESPOL(ND,I3,D3,BES3,LVEC)
      CALL BESPOL(ND,I4,D4,BES4,LVEC)
C
      DO 150 I=1,LVEC
C
C     RECURSION RELATION FOR  CHI
C
      CHIOLD=CHIDM1(I)
      CHIDM1(I)=CHID(I)
      CHID(I)=TRACE(I)*CHID(I) - BSQ(I)*CHIOLD
C
C     SUMMATION
C
      F5=CHID(I)*BES1(I)*BES2(I)*BES3(I)*BES4(I)
      F6(I)=RD(ND)*F5
      SUM1(I)=SUM1(I) + SD(ND)*F5
      SUM2(I)=SUM2(I) + F6(I)
150   CONTINUE
C
100   CONTINUE
C
C     END OF LOOP OVER D
C
      DO 200 I=1,LVEC
      PLINT(I)=SUM1(I)/SUM2(I)
      IF(ABS(F6(I)) .GT. SUM2(I)*ERR4) NERR4=NERR4+1
200   CONTINUE
      END
C*************************INT2V*****************************
C
C     COMPUTE THE 2-LINK INTEGRAL EXPRESSION FOR THE CORNER
C          (VECTORISED VERSION)
C
      SUBROUTINE INT2V(B1,B2,W,AV)
*     ****************
C
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(ERROR=1.E-8)
      PARAMETER(ERROR2=1.E-7,NDMAX2=10)
C
      REAL*8 BCOEFF,BETINV,RD,SD
      REAL*8 B1,B2,W,AV,TRACE,BSQ,D1,D2,SUM1,SUM2,PHID,CHID
      REAL*8 C1,C2,F4,DENOM
      REAL*8 B1SQ,B2SQ,D,PHIOLD,F2,F3,ERR2,ZERO,ONE
      PARAMETER(ERR2=ERROR2,ZERO=0.D0,ONE=1.D0)
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/INTBLK/BETINV,NERR4,NERR2
      COMMON/BESS3/RD(20),SD(20)
      COMMON/BESS4/BCOEFF
C
      DIMENSION B1(MAXVEC,0:3),B2(MAXVEC,0:3)
      DIMENSION W(MAXVEC,0:3),AV(MAXVEC,0:3)
C
      DIMENSION TRACE(MAXVEC),BSQ(MAXVEC)
      DIMENSION I1(MAXVEC),I2(MAXVEC),D1(MAXVEC),D2(MAXVEC)
      DIMENSION SUM1(MAXVEC),SUM2(MAXVEC),PHID(MAXVEC),CHID(MAXVEC)
      DIMENSION C1(MAXVEC),C2(MAXVEC)
      DIMENSION F4(MAXVEC),DENOM(MAXVEC)
C
C     PRODUCT B2*B1
C
      CALL MATMAT(B2,B1,AV,LVEC,MAXVEC)
C
      DO 50 I=1,LVEC
C
      TRACE(I)= AV(I,0)*W(I,0) + AV(I,1)*W(I,1)
     X        + AV(I,2)*W(I,2) + AV(I,3)*W(I,3)
C
      B1SQ=B1(I,0)**2 + B1(I,1)**2 + B1(I,2)**2 + B1(I,3)**2
      B2SQ=B2(I,0)**2 + B2(I,1)**2 + B2(I,2)**2 + B2(I,3)**2
C
      BSQ(I)=B1SQ*B2SQ
C
      D1(I)=BCOEFF*B1SQ
      I1(I)=D1(I)
      D1(I)=D1(I)-I1(I)
C
      D2(I)=BCOEFF*B2SQ
      I2(I)=D2(I)
      D2(I)=D2(I)-I2(I)
C
C     INITIALISE THE SUMS AND RECURSION RELATIONS FOR D=1
C
      SUM1(I)=SD(1)
      SUM2(I)=ZERO
      DENOM(I)=ONE
C
      CHID(I)=ONE
      PHID(I)=ZERO
C
50    CONTINUE
C
C     LOOP OVER D=2,3,4,...,NDMAX2
C
      DO 100 ND=2,NDMAX2
      D=ND
C
      CALL BESPOL(ND,I1,D1,C1,LVEC)
      CALL BESPOL(ND,I2,D2,C2,LVEC)
C
      DO 150 I=1,LVEC
C
C     RECURSION RELATION FOR  CHI  AND  PHI
C
      PHIOLD=PHID(I)
      PHID(I)=D*CHID(I) + TRACE(I)*PHID(I)
      CHID(I)=(TRACE(I)*PHID(I) - BSQ(I)*PHIOLD) /(D-ONE)
C
C     SUMMATION
C
      F2=C1(I)*C2(I)
      F3=CHID(I)*F2
      F4(I)=RD(ND)*F3
      SUM1(I)=SUM1(I) + SD(ND)*F3
      SUM2(I)=SUM2(I) + RD(ND)*PHID(I)*F2
      DENOM(I)=DENOM(I) + F4(I)
C
150   CONTINUE
C
100   CONTINUE
C
C     END OF LOOP OVER D
C
C     COMPUTE THE CORNER AVERAGE
C
      DO 200 I=1,LVEC
C
      C2(I)=BETINV*SUM2(I)/DENOM(I)
      C1(I)=SUM1(I)/DENOM(I) - TRACE(I)*C2(I)
C
      AV(I,0)=C1(I)*W(I,0) + C2(I)*AV(I,0)
      AV(I,1)=C1(I)*W(I,1) + C2(I)*AV(I,1)
      AV(I,2)=C1(I)*W(I,2) + C2(I)*AV(I,2)
      AV(I,3)=C1(I)*W(I,3) + C2(I)*AV(I,3)
C
      IF(ABS(F4(I)) .GT. DENOM(I)*ERR2) NERR2=NERR2+1
C
200   CONTINUE
      END
C*************************CORR***************************
C
C     ACCUMULATE THE GLUEBALL CORRELATIONS
C       (SUBROUTINE SPECIALISED TO 4 DIMENSIONS)
C
      SUBROUTINE CORR
*     ***************
C
      PARAMETER(MAXVEC=8**3*8+7)
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
C
      REAL*8 S110,S112,S110I,S112I
      REAL*8 S220,S222,S220I,S222I
      REAL*8 W11SUM(20,3),W11INT(20,3)
      REAL*8 W22SUM(20,3),W22INT(20,3)
      REAL*8 OP110(20),OP112(20),OP110I(20),OP112I(20)
      REAL*8 OP220(20),OP222(20),OP220I(20),OP222I(20)
      REAL*8 PR110(0:20/2),PR112(0:20/2)
      REAL*8 PR110I(0:20/2),PR112I(0:20/2)
      REAL*8 PR220(0:20/2),PR222(0:20/2)
      REAL*8 PR220I(0:20/2),PR222I(0:20/2)
      REAL*8 VECT(2+20/2,8)
C
      NT=LSIZE(NDIM)
      MAXDST=NT/2
      NDIST=MAXDST+1
C
      S110=0.D0
      S112=0.D0
      S110I=0.D0
      S112I=0.D0
C
      S220=0.D0
      S222=0.D0
      S220I=0.D0
      S222I=0.D0
C
      CALL LOOPS(W11SUM,W11INT,W22SUM,W22INT)
C
      DO 10 IT=1,NT
C
      OP110(IT)=W11SUM(IT,3)+W11SUM(IT,2)+W11SUM(IT,1)
      OP112(IT)=W11SUM(IT,3)+W11SUM(IT,2)-2*W11SUM(IT,1)
      OP110I(IT)=W11INT(IT,3)+W11INT(IT,2)+W11INT(IT,1)
      OP112I(IT)=W11INT(IT,3)+W11INT(IT,2)-2*W11INT(IT,1)
C
      OP220(IT)=W22SUM(IT,3)+W22SUM(IT,2)+W22SUM(IT,1)
      OP222(IT)=W22SUM(IT,3)+W22SUM(IT,2)-2*W22SUM(IT,1)
      OP220I(IT)=W22INT(IT,3)+W22INT(IT,2)+W22INT(IT,1)
      OP222I(IT)=W22INT(IT,3)+W22INT(IT,2)-2*W22INT(IT,1)
C
      S110=S110 + OP110(IT)
      S112=S112 + OP112(IT)
      S110I=S110I + OP110I(IT)
      S112I=S112I + OP112I(IT)
C
      S220=S220 + OP220(IT)
      S222=S222 + OP222(IT)
      S220I=S220I + OP220I(IT)
      S222I=S222I + OP222I(IT)
10    CONTINUE
C
      CALL ACCUM2(10,S110/NT)
      CALL ACCUM2(11,S112/NT)
      CALL ACCUM2(12,S110I/NT)
      CALL ACCUM2(13,S112I/NT)
C
      CALL ACCUM2(14,S220/NT)
      CALL ACCUM2(15,S222/NT)
      CALL ACCUM2(16,S220I/NT)
      CALL ACCUM2(17,S222I/NT)
C
      VECT(1,1)=S110/NT
      VECT(1,2)=S112/NT
      VECT(1,3)=S110I/NT
      VECT(1,4)=S112I/NT
      VECT(1,5)=S220/NT
      VECT(1,6)=S222/NT
      VECT(1,7)=S220I/NT
      VECT(1,8)=S222I/NT
C
      DO 20 IDIST=0,MAXDST
      PR110(IDIST)=0.D0
      PR112(IDIST)=0.D0
      PR110I(IDIST)=0.D0
      PR112I(IDIST)=0.D0
C
      PR220(IDIST)=0.D0
      PR222(IDIST)=0.D0
      PR220I(IDIST)=0.D0
      PR222I(IDIST)=0.D0
20    CONTINUE
C
      DO 30 IT1=1,NT
      DO 30 IT2=IT1,NT
      IDIST=MIN(IT2-IT1,NT+IT1-IT2)
      PR110(IDIST)=PR110(IDIST) + OP110(IT1)*OP110(IT2)
      PR112(IDIST)=PR112(IDIST) + OP112(IT1)*OP112(IT2)
      PR110I(IDIST)=PR110I(IDIST) + OP110I(IT1)*OP110I(IT2)
      PR112I(IDIST)=PR112I(IDIST) + OP112I(IT1)*OP112I(IT2)
C
      PR220(IDIST)=PR220(IDIST) + OP220(IT1)*OP220(IT2)
      PR222(IDIST)=PR222(IDIST) + OP222(IT1)*OP222(IT2)
      PR220I(IDIST)=PR220I(IDIST) + OP220I(IT1)*OP220I(IT2)
      PR222I(IDIST)=PR222I(IDIST) + OP222I(IT1)*OP222I(IT2)
30    CONTINUE
C
      DO 40 IDIST=0,MAXDST
      IND=IDIST+1
      IND1=IND+1
      FACT=1.D0/NT
      IF(2*IDIST.EQ.NT) FACT=2.D0/NT
C
      CALL ACCUM2(17+IND,PR110(IDIST)*FACT)
      CALL ACCUM2(17+NDIST+IND,PR112(IDIST)*FACT)
      CALL ACCUM2(17+2*NDIST+IND,PR110I(IDIST)*FACT)
      CALL ACCUM2(17+3*NDIST+IND,PR112I(IDIST)*FACT)
C
      CALL ACCUM2(17+4*NDIST+IND,PR220(IDIST)*FACT)
      CALL ACCUM2(17+5*NDIST+IND,PR222(IDIST)*FACT)
      CALL ACCUM2(17+6*NDIST+IND,PR220I(IDIST)*FACT)
      CALL ACCUM2(17+7*NDIST+IND,PR222I(IDIST)*FACT)
C
      VECT(IND1,1)=PR110(IDIST)*FACT
      VECT(IND1,2)=PR112(IDIST)*FACT
      VECT(IND1,3)=PR110I(IDIST)*FACT
      VECT(IND1,4)=PR112I(IDIST)*FACT
      VECT(IND1,5)=PR220(IDIST)*FACT
      VECT(IND1,6)=PR222(IDIST)*FACT
      VECT(IND1,7)=PR220I(IDIST)*FACT
      VECT(IND1,8)=PR222I(IDIST)*FACT
40    CONTINUE
C
      DO 50 IBL=1,8
      CALL ACCUM4(IBL,VECT(1,IBL))
50    CONTINUE
      END
C*************************LOOPS**************************
C
C     MEASUREMENT OF SUMS OF LOOPS (VECTORISED)
C          W11SUM(IT,IPLANE) : NAIVE PLAQUETTE EXPRESSION
C          W11INT(IT,IPLANE) : 4-LINK INTEGRAL
C          W22SUM(IT,IPLANE) : NAIVE 2X2 LOOP
C          W22INT(IT,IPLANE) : 2X2 LOOP INTEGRAL (4 CORNERS)
C        IT: TIME SLICE
C        IPLANE=1,2,...  FOR SPATIAL PLANES (1,2),(1,3),...
C
      SUBROUTINE LOOPS(W11SUM,W11INT,W22SUM,W22INT)
*     ****************
C
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(NST=(4-1)*(4-2),MAXSL=2**(4-1))
C
      REAL*8 U,BETINV
      REAL*8 R,S,W1,W2,WP,WR,WL,W,WLOOP1,WLOOP2
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
      COMMON/PAR3/BETA,GAMMA,ITER,IFREQ,MFREQ
      COMMON/GEOM1/NEWLAT(2**4,-4:4)
      COMMON/GEOM2/ICODE(2**4,-4:4)
      COMMON/INTBLK/BETINV,NERR4,NERR2
      COMMON/CONF/U(MAXVEC,0:3,2**4,4)
C
      DIMENSION R(MAXVEC,0:3,MAXSL,NST),S(MAXVEC,0:3,MAXSL,NST)
      DIMENSION W1(MAXVEC,0:3),W2(MAXVEC,0:3),WP(MAXVEC,0:3)
      DIMENSION WR(MAXVEC,0:3),WL(MAXVEC,0:3),W(MAXVEC,0:3)
      DIMENSION WLOOP1(MAXVEC),WLOOP2(MAXVEC)
C
      DIMENSION MPL(4,4),MPLANE(4,4)
C
      REAL*8 W11SUM(20,*),W11INT(20,*)
      REAL*8 W22SUM(20,*),W22INT(20,*)
C
      NT=LSIZE(NDIM)
      NSPLPP=(NDSP*(NDSP-1))/2
      LSPVEC=NVOL/2**NDSP
C
      BETINV=1.D0/BETA
C
      IPLANE=0
      DO 10 I1=1,NDSP-1
      DO 10 I2=I1+1,NDSP
      IPLANE=IPLANE+1
      MPLANE(I1,I2)=IPLANE
      MPLANE(I2,I1)=IPLANE
10    CONTINUE
C
      I1I2=0
      DO 20 I1=1,NDSP
      DO 20 I2=1,NDSP
      IF(I1.EQ.I2)GOTO 20
      I1I2=I1I2+1
      MPL(I1,I2)=I1I2
20    CONTINUE
C
      DO 40 IPLANE=1,NSPLPP
      DO 40 IT=1,NT
      W11SUM(IT,IPLANE)=0.D0
      W11INT(IT,IPLANE)=0.D0
      W22SUM(IT,IPLANE)=0.D0
      W22INT(IT,IPLANE)=0.D0
40    CONTINUE
C
C     LOOP OVER SUB-LATTICES WITH ODD AND EVEN TIME-SLICES
C
      DO 900 IFIRST=1,2
C
      IADD=0
      IF(IFIRST.EQ.2)IADD=2**NDSP
C.....................................................
C     COMPUTE THE SPACE-LIKE PRODUCTS R AND S
C
C     LOOP OVER SPATIAL PLAQUETTES (ILAT2,I1,I2)
C
      DO 100 ILAT2=1,2**NDSP
      DO 100 I1=1,NDSP-1
      DO 100 I2=I1+1,NDSP
C
C               I1
C                .
C                .       5   6
C                .       2   3
C                . . . . . . . . . . I2
C
C     ILAT2  REFERS TO POINT '2'
C
      I1I2=MPL(I1,I2)
      I2I1=MPL(I2,I1)
C
      ILAT3=NEWLAT(ILAT2,I2)
      ILAT5=NEWLAT(ILAT2,I1)
C
      IU2=ILAT2+IADD
      IU3=ILAT3+IADD
      IU5=ILAT5+IADD
C
      IC23=ICODE(ILAT2,I2)
      IC25=ICODE(ILAT2,I1)
      IC32=ICODE(ILAT3,-I2)
      IC52=ICODE(ILAT5,-I1)
C
C     PRODUCT '2-3-6' IN W1
C
      IF(IC23.EQ.0)THEN
      CALL MATMAT(U(1,0,IU3,I1),U(1,0,IU2,I2),W1,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU3,I1),IC23,WP,LVEC)
      CALL MATMAT(WP,U(1,0,IU2,I2),W1,LVEC,MAXVEC)
      ENDIF
C
C     PRODUCT '2-3-6-5' IN R
C
      IF(IC25.EQ.0)THEN
      CALL ADJMAT(U(1,0,IU5,I2),W1,R(1,0,ILAT2,I1I2),LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU5,I2),IC25,WP,LVEC)
      CALL ADJMAT(WP,W1,R(1,0,ILAT2,I1I2),LVEC,MAXVEC)
      ENDIF
C
C     PRODUCT '5-2-3-6' IN S
C
      IF(IC52.EQ.0)THEN
      CALL MATADJ(W1,U(1,0,IU2,I1),S(1,0,ILAT5,I2I1),LVEC,MAXVEC)
      ELSE
      CALL MATADJ(W1,U(1,0,IU2,I1),WP,LVEC,MAXVEC)
      CALL PERM(WP,IC52,S(1,0,ILAT5,I2I1),LVEC)
      ENDIF
C
C     PRODUCT '2-5-6' IN W2
C
      IF(IC25.EQ.0)THEN
      CALL MATMAT(U(1,0,IU5,I2),U(1,0,IU2,I1),W2,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU5,I2),IC25,WP,LVEC)
      CALL MATMAT(WP,U(1,0,IU2,I1),W2,LVEC,MAXVEC)
      ENDIF
C
C     PRODUCT '2-5-6-3' IN R
C
      IF(IC23.EQ.0)THEN
      CALL ADJMAT(U(1,0,IU3,I1),W2,R(1,0,ILAT2,I2I1),LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU3,I1),IC23,WP,LVEC)
      CALL ADJMAT(WP,W2,R(1,0,ILAT2,I2I1),LVEC,MAXVEC)
      ENDIF
C
C     PRODUCT '3-2-5-6' IN S
C
      IF(IC32.EQ.0)THEN
      CALL MATADJ(W2,U(1,0,IU2,I2),S(1,0,ILAT3,I1I2),LVEC,MAXVEC)
      ELSE
      CALL MATADJ(W2,U(1,0,IU2,I2),WP,LVEC,MAXVEC)
      CALL PERM(WP,IC32,S(1,0,ILAT3,I1I2),LVEC)
      ENDIF
C
C     COMPUTE AND ACCUMULATE THE SPATIAL PLAQUETTES
C
      DO 110 I=1,LVEC
      WLOOP1(I)= W1(I,0)*W2(I,0) + W1(I,1)*W2(I,1)
     X         + W1(I,2)*W2(I,2) + W1(I,3)*W2(I,3)
110   CONTINUE
C
      IPLANE=MPLANE(I1,I2)
      N=0
      DO 120 IT=IFIRST,NT,2
      DO 130 I=1,LSPVEC
      W11SUM(IT,IPLANE)=W11SUM(IT,IPLANE) + DBLE(WLOOP1(N+I))
130   CONTINUE
      N=N+LSPVEC
120   CONTINUE
C
100   CONTINUE
C
C     END OF LOOP OVER SPATIAL PLAQUETTES
C.....................................................
C     COMPUTE THE SPACE-LIKE 2X2 LOOPS
C
C     LOOP OVER SPATIAL PLANES (ILAT2,I1,I2)
C
      DO 200 ILAT2=1,2**NDSP
      DO 200 I1=1,NDSP-1
      DO 200 I2=I1+1,NDSP
C
C               I1
C                .    7   8   9
C                .    4   5   6
C                .    1   2   3
C                . . . . . . . . . . I2
C
C     ILAT2  REFERS TO POINT '2'
C
      I1I2=MPL(I1,I2)
C
      ILAT5=NEWLAT(ILAT2,I1)
      IC25=ICODE(ILAT2,I1)
C
C     PRODUCT '5-4-1-2-3-6-5' IN W1
C
      CALL MATADJ(R(1,0,ILAT2,I1I2),S(1,0,ILAT2,I1I2),W1,LVEC,MAXVEC)
C
C     PRODUCT '5-4-7-8-9-6-5' IN W2
C
      IF(IC25.EQ.0)THEN
      CALL ADJMAT(R(1,0,ILAT5,I1I2),S(1,0,ILAT5,I1I2),W2,LVEC,MAXVEC)
      ELSE
      CALL ADJMAT(R(1,0,ILAT5,I1I2),S(1,0,ILAT5,I1I2),WP,LVEC,MAXVEC)
      CALL PERM(WP,IC25,W2,LVEC)
      ENDIF
C
C     COMPUTE AND ACCUMULATE THE SPATIAL 2X2 LOOPS
C
      DO 210 I=1,LVEC
      WLOOP2(I)= W1(I,0)*W2(I,0) + W1(I,1)*W2(I,1)
     X         + W1(I,2)*W2(I,2) + W1(I,3)*W2(I,3)
210   CONTINUE
C
      IPLANE=MPLANE(I1,I2)
      N=0
      DO 220 IT=IFIRST,NT,2
      DO 230 I=1,LSPVEC
      W22SUM(IT,IPLANE)=W22SUM(IT,IPLANE) + DBLE(WLOOP2(N+I))
230   CONTINUE
      N=N+LSPVEC
220   CONTINUE
C
200   CONTINUE
C
C     END OF LOOP OVER SPATIAL PLANES
C..........................................................
C     COMPUTE THE TIME-LIKE PRODUCTS
C       AND REPLACE  R  AND  S  BY THE COMPLEMENTARY
C       BACKGROUND EXPRESSIONS  (SUM - R)  AND  (SUM - S)
C
C     LOOP OVER SPATIAL LINKS (ILAT,I1)
C
      DO 300 ILAT=1,2**NDSP
      DO 300 I1=1,NDSP
C
      CALL INTACT(ILAT+IADD,I1, NDIM,W1)
      CALL INTACT(ILAT+IADD,I1,-NDIM,W2)
      DO 310 I=1,LVEC
      W1(I,0)=W1(I,0) + W2(I,0)
      W1(I,1)=W1(I,1) + W2(I,1)
      W1(I,2)=W1(I,2) + W2(I,2)
      W1(I,3)=W1(I,3) + W2(I,3)
310   CONTINUE
C
      DO 320 I2=1,NDSP
      IF(I1.EQ.I2)GOTO 320
      I1I2=MPL(I1,I2)
      DO 321 I=1,LVEC
      W1(I,0)=W1(I,0) + R(I,0,ILAT,I1I2) + S(I,0,ILAT,I1I2)
      W1(I,1)=W1(I,1) + R(I,1,ILAT,I1I2) + S(I,1,ILAT,I1I2)
      W1(I,2)=W1(I,2) + R(I,2,ILAT,I1I2) + S(I,2,ILAT,I1I2)
      W1(I,3)=W1(I,3) + R(I,3,ILAT,I1I2) + S(I,3,ILAT,I1I2)
321   CONTINUE
320   CONTINUE
C
      DO 330 I2=1,NDSP
      IF(I1.EQ.I2)GOTO 330
      I1I2=MPL(I1,I2)
      DO 331 I=1,LVEC
      R(I,0,ILAT,I1I2)=W1(I,0) - R(I,0,ILAT,I1I2)
      R(I,1,ILAT,I1I2)=W1(I,1) - R(I,1,ILAT,I1I2)
      R(I,2,ILAT,I1I2)=W1(I,2) - R(I,2,ILAT,I1I2)
      R(I,3,ILAT,I1I2)=W1(I,3) - R(I,3,ILAT,I1I2)
      S(I,0,ILAT,I1I2)=W1(I,0) - S(I,0,ILAT,I1I2)
      S(I,1,ILAT,I1I2)=W1(I,1) - S(I,1,ILAT,I1I2)
      S(I,2,ILAT,I1I2)=W1(I,2) - S(I,2,ILAT,I1I2)
      S(I,3,ILAT,I1I2)=W1(I,3) - S(I,3,ILAT,I1I2)
331   CONTINUE
330   CONTINUE
C
300   CONTINUE
C
C     END OF LOOP OVER SPATIAL LINKS
C..........................................................
C     COMPUTE THE SPACE-LIKE 1X1 AND 2X2 INTEGRALS
C
C     LOOP OVER SPATIAL PLANES (ILAT1,I1,I2)
C
      DO 400 ILAT1=1,2**NDSP
      DO 400 I1=1,NDSP-1
      DO 400 I2=I1+1,NDSP
C
C               I1
C                .    7   8   9
C                .    4   5   6
C                .    1   2   3
C                . . . . . . . . . . I2
C
C     ILAT1  REFERS TO POINT '1'
C
      I1I2=MPL(I1,I2)
      I2I1=MPL(I2,I1)
C
      ILAT2=NEWLAT(ILAT1,I2)
      ILAT3=NEWLAT(ILAT2,I2)
      ILAT4=NEWLAT(ILAT1,I1)
      ILAT5=NEWLAT(ILAT2,I1)
      ILAT6=NEWLAT(ILAT3,I1)
      ILAT7=NEWLAT(ILAT4,I1)
      ILAT8=NEWLAT(ILAT5,I1)
C
      IU2=ILAT2+IADD
      IU4=ILAT4+IADD
      IU5=ILAT5+IADD
C
      IC12=ICODE(ILAT1,I2)
      IC14=ICODE(ILAT1,I1)
      IC23=ICODE(ILAT2,I2)
      IC25=ICODE(ILAT2,I1)
      IC45=ICODE(ILAT4,I2)
      IC47=ICODE(ILAT4,I1)
      IC56=ICODE(ILAT5,I2)
      IC58=ICODE(ILAT5,I1)
C
C     COMPUTE THE PLAQUETTE INTEGRAL FOR   4  5
C                                          1  2
C
      IF(IC12.EQ.0)THEN
      IF(IC14.EQ.0)THEN
      CALL INT4V(R(1,0,ILAT1,I2I1),S(1,0,ILAT2,I1I2),
     X           S(1,0,ILAT4,I2I1),R(1,0,ILAT1,I1I2),WLOOP1)
      ELSE
      CALL PERM(S(1,0,ILAT4,I2I1),IC14,WP,LVEC)
      CALL INT4V(R(1,0,ILAT1,I2I1),S(1,0,ILAT2,I1I2),
     X           WP               ,R(1,0,ILAT1,I1I2),WLOOP1)
      ENDIF
      ELSE
      CALL PERM(S(1,0,ILAT2,I1I2),IC12,W1,LVEC)
      IF(IC14.EQ.0)THEN
      CALL INT4V(R(1,0,ILAT1,I2I1),W1,
     X           S(1,0,ILAT4,I2I1),R(1,0,ILAT1,I1I2),WLOOP1)
      ELSE
      CALL PERM(S(1,0,ILAT4,I2I1),IC14,WP,LVEC)
      CALL INT4V(R(1,0,ILAT1,I2I1),W1,
     X           WP               ,R(1,0,ILAT1,I1I2),WLOOP1)
      ENDIF
      ENDIF
C
C     COMPUTE THE 2X2 LOOP INTEGRALS
C
C     CORNER '2-3-6' IN  W1  (REFERENCE POINT '2')
C
      IF(IC25.EQ.0)THEN
      CALL MATMAT(U(1,0,IU5,I2),U(1,0,IU2,I1),W,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU5,I2),IC25,WP,LVEC)
      CALL MATMAT(WP,U(1,0,IU2,I1),W,LVEC,MAXVEC)
      ENDIF
C
      IF(IC23.EQ.0)THEN
      CALL INT2V(R(1,0,ILAT2,I2I1),S(1,0,ILAT3,I1I2),W,W1)
      ELSE
      CALL PERM(S(1,0,ILAT3,I1I2),IC23,WP,LVEC)
      CALL INT2V(R(1,0,ILAT2,I2I1),WP,W,W1)
      ENDIF
C
C     CORNER '6-9-8' IN  W2  (REFERENCE POINT '5')
C
      CALL MATADJ(U(1,0,IU5,I1),U(1,0,IU5,I2),W,LVEC,MAXVEC)
C
      DO 450 I=1,LVEC
      WR(I,0)= S(I,0,ILAT8,I2I1)
      WR(I,1)=-S(I,1,ILAT8,I2I1)
      WR(I,2)=-S(I,2,ILAT8,I2I1)
      WR(I,3)=-S(I,3,ILAT8,I2I1)
450   CONTINUE
      IF(IC58.EQ.0)THEN
      IF(IC56.EQ.0)THEN
      CALL INT2V(S(1,0,ILAT6,I1I2),WR,W,W2)
      ELSE
      CALL PERM(S(1,0,ILAT6,I1I2),IC56,WP,LVEC)
      CALL INT2V(WP,WR,W,W2)
      ENDIF
      ELSE
      CALL PERM(WR,IC58,WL,LVEC)
      IF(IC56.EQ.0)THEN
      CALL INT2V(S(1,0,ILAT6,I1I2),WL,W,W2)
      ELSE
      CALL PERM(S(1,0,ILAT6,I1I2),IC56,WP,LVEC)
      CALL INT2V(WP,WL,W,W2)
      ENDIF
      ENDIF
C
C     PRODUCT '2-3-6-9-8' IN  WR  (REFERENCE POINT '2')
C
      IF(IC25.EQ.0)THEN
      CALL MATMAT(W2,W1,WR,LVEC,MAXVEC)
      ELSE
      CALL PERM(W2,IC25,WP,LVEC)
      CALL MATMAT(WP,W1,WR,LVEC,MAXVEC)
      ENDIF
C
C     CORNER '2-1-4' IN  W1  (REFERENCE POINT '1')
C
      IF(IC12.EQ.0)THEN
      IF(IC14.EQ.0)THEN
      CALL ADJMAT(U(1,0,IU4,I2),U(1,0,IU2,I1),W,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU4,I2),IC14,WP,LVEC)
      CALL ADJMAT(WP,U(1,0,IU2,I1),W,LVEC,MAXVEC)
      ENDIF
      ELSE
      CALL PERM(U(1,0,IU2,I1),IC12,W2,LVEC)
      IF(IC14.EQ.0)THEN
      CALL ADJMAT(U(1,0,IU4,I2),W2,W,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU4,I2),IC14,WP,LVEC)
      CALL ADJMAT(WP,W2,W,LVEC,MAXVEC)
      ENDIF
      ENDIF
C
      DO 460 I=1,LVEC
      WL(I,0)= R(I,0,ILAT1,I2I1)
      WL(I,1)=-R(I,1,ILAT1,I2I1)
      WL(I,2)=-R(I,2,ILAT1,I2I1)
      WL(I,3)=-R(I,3,ILAT1,I2I1)
460   CONTINUE
      CALL INT2V(WL,R(1,0,ILAT1,I1I2),W,W1)
C
C     CORNER '4-7-8' IN  W2  (REFERENCE POINT '4')
C
      IF(IC45.EQ.0)THEN
      CALL MATMAT(U(1,0,IU5,I1),U(1,0,IU4,I2),W,LVEC,MAXVEC)
      ELSE
      CALL PERM(U(1,0,IU5,I1),IC45,WP,LVEC)
      CALL MATMAT(WP,U(1,0,IU4,I2),W,LVEC,MAXVEC)
      ENDIF
C
      IF(IC47.EQ.0)THEN
      CALL INT2V(R(1,0,ILAT4,I1I2),S(1,0,ILAT7,I2I1),W,W2)
      ELSE
      CALL PERM(S(1,0,ILAT7,I2I1),IC47,WP,LVEC)
      CALL INT2V(R(1,0,ILAT4,I1I2),WP,W,W2)
      ENDIF
C
C     PRODUCT '2-1-4-7-8' IN  WL  (REFERENCE POINT '1')
C
      IF(IC14.EQ.0)THEN
      CALL MATMAT(W2,W1,WL,LVEC,MAXVEC)
      ELSE
      CALL PERM(W2,IC14,WP,LVEC)
      CALL MATMAT(WP,W1,WL,LVEC,MAXVEC)
      ENDIF
C
C     COMPUTE THE TRACE (REFERENCE POINT '1')
C
      IF(IC12.EQ.0)THEN
      DO 470 I=1,LVEC
      WLOOP2(I)= WR(I,0)*WL(I,0) + WR(I,1)*WL(I,1)
     X         + WR(I,2)*WL(I,2) + WR(I,3)*WL(I,3)
470   CONTINUE
      ELSE
      CALL PERM(WR,IC12,WP,LVEC)
      DO 480 I=1,LVEC
      WLOOP2(I)= WP(I,0)*WL(I,0) + WP(I,1)*WL(I,1)
     X         + WP(I,2)*WL(I,2) + WP(I,3)*WL(I,3)
480   CONTINUE
      ENDIF
C
C     ACCUMULATE THE 1X1 AND 2X2 LOOP INTEGRALS
C
      IPLANE=MPLANE(I1,I2)
      N=0
      DO 490 IT=IFIRST,NT,2
      DO 495 I=1,LSPVEC
      W11INT(IT,IPLANE)=W11INT(IT,IPLANE) + DBLE(WLOOP1(N+I))
      W22INT(IT,IPLANE)=W22INT(IT,IPLANE) + DBLE(WLOOP2(N+I))
495   CONTINUE
      N=N+LSPVEC
490   CONTINUE
C
400   CONTINUE
C
C     END OF LOOP OVER SPATIAL PLANES
C...........................................................
C
900   CONTINUE
C
C     END OF LOOP OVER ODD/EVEN SUB-LATTICES
C
      END
C*************************STAT2NN***************************
C                                               B.BUNK 1985
C     RELAXATION ANALYSIS
C
C     FOR INTERNAL COUNTERS IN REAL*8
C        REMOVE  C+++++
C
C     CLEAR COUNTERS
C
      SUBROUTINE CLEAR2(NVAR1,NCORR1)
*     *****************
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXVAR=25+4*20,MAXCOR=10)
* new common order against missalignments          a.w.
      COMMON/ST2/NVAR,NCORR,SAV(MAXVAR),F(MAXVAR,MAXCOR)
      COMMON/ST2/SUM(MAXVAR,0:MAXCOR),H(MAXVAR,MAXCOR),NDAT(MAXVAR)
      REAL*8 SAV,F,SUM,H
      SAVE/ST2/
      IF(NVAR1.GT.MAXVAR)THEN
      WRITE(*,'('' ERROR IN STAT2: NVAR ='',I5)')NVAR1
      STOP
      ENDIF
      IF(NCORR1.GT.MAXCOR)THEN
      WRITE(*,'('' ERROR IN STAT2: NCORR ='',I5)')NCORR1
      STOP
      ENDIF
      NVAR=NVAR1
      NCORR=NCORR1
      DO 10 IVAR=1,NVAR
      NDAT(IVAR)=0
      SAV(IVAR)=0.D0
      SUM(IVAR,0)=0.D0
      DO 10 ICORR=1,NCORR
      F(IVAR,ICORR)=0.D0
      SUM(IVAR,ICORR)=0.D0
10    CONTINUE
      END
C
C     ACCUMULATE DATA
C
      SUBROUTINE ACCUM2(IVAR,DFN)
*     *****************
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXVAR=25+4*20,MAXCOR=10)
* new common order against missalignments          a.w.
      COMMON/ST2/NVAR,NCORR,SAV(MAXVAR),F(MAXVAR,MAXCOR)
      COMMON/ST2/SUM(MAXVAR,0:MAXCOR),H(MAXVAR,MAXCOR),NDAT(MAXVAR)
      REAL*8 SAV,F,SUM,H,DFN
      SAVE/ST2/
      IF(IVAR.LT.1 .OR. IVAR.GT.NVAR)RETURN
      NDAT(IVAR)=NDAT(IVAR)+1
C+++++      DFN=FN
      SAV(IVAR)=SAV(IVAR)+DFN
      SUM(IVAR,0)=SUM(IVAR,0)+DFN*DFN
      DO 10 ICORR=1,NCORR
10    SUM(IVAR,ICORR)=SUM(IVAR,ICORR)+DFN*F(IVAR,ICORR)
      DO 20 ICORR=NCORR,2,-1
20    F(IVAR,ICORR)=F(IVAR,ICORR-1)
      F(IVAR,1)=DFN
      IF(NDAT(IVAR).LE.NCORR) H(IVAR,NDAT(IVAR))=DFN
      END
C
C     COMPUTE AVERAGES
C
      FUNCTION AVER2(IVAR)
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXVAR=25+4*20,MAXCOR=10)
* new common order against missalignments          a.w.
      COMMON/ST2/NVAR,NCORR,SAV(MAXVAR),F(MAXVAR,MAXCOR)
      COMMON/ST2/SUM(MAXVAR,0:MAXCOR),H(MAXVAR,MAXCOR),NDAT(MAXVAR)
      REAL*8 SAV,F,SUM,H
      SAVE/ST2/
      AVER2=0.D0
      ND=NDAT(IVAR)
      IF(IVAR.LT.1 .OR. IVAR.GT.NVAR .OR. ND.LT.1)RETURN
      AVER2=SAV(IVAR)/ND
      END
C
C     STATISTICAL ANALYSIS, INCLUDING RELAXATION EFFECTS
C
      SUBROUTINE STAT2(IVAR,AV,SIG,RATIO,IUNIT)
*     ****************
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXVAR=25+4*20,MAXCOR=10)
* new common order against missalignments          a.w.
      COMMON/ST2/NVAR,NCORR,SAV(MAXVAR),F(MAXVAR,MAXCOR)
      COMMON/ST2/SUM(MAXVAR,0:MAXCOR),H(MAXVAR,MAXCOR),NDAT(MAXVAR)
      REAL*8 SAV,F,SUM,H,AV,S1,S2
      SAVE/ST2/
      DIMENSION GAM(0:MAXCOR),R(MAXCOR),T(MAXCOR)
      AV=0.D0
      SIG=0.D0
      RATIO=0.D0
      IF(IVAR.LT.1 .OR. IVAR.GT.NVAR)RETURN
      ND=NDAT(IVAR)
      NMAX=MIN(ND-1,NCORR)
      IF(IUNIT.GT.0)THEN
      WRITE(IUNIT,1000) ND,NMAX
1000  FORMAT(' RELAXATION ANALYSIS WITH NDAT =',I5,'   ,NCORR =',I5)
      WRITE(IUNIT,'('' INTERNAL REAL*8'')')
      ENDIF
      IF(ND.LT.1)RETURN
      AV=SAV(IVAR)/ND
      IF(IUNIT.GT.0) WRITE(IUNIT,'('' AVERAGE ='',F12.6)') AV
      GAM(0)=SUM(IVAR,0)/ND-AV*AV
      SIG0=SQRT(MAX( GAM(0)/ND ,0.))
      IF(NMAX.LT.1)GOTO 90
C
      N1=MAX(NMAX/5,1)
      N2=NMAX/2
      SIG=0.D0
      IFLAG=NMAX+1
      U=GAM(0)
      S1=0.D0
      S2=0.D0
C
      DO 10 N=1,NMAX
      S1=S1+H(IVAR,N)
      S2=S2+F(IVAR,N)
      GAM(N)=( SUM(IVAR,N)+(S1+S2-(ND+N)*AV)*AV )/(ND-N)
C
      R(N)=0.D0
      IF(U.NE.0.) R(N)=GAM(N)/U
      U=GAM(N)
C
      T(N)=0.D0
      IF(IFLAG.LT.N)GOTO 10
      IF(R(N).GT.0.D0 .AND. R(N).LT.1.D0)THEN
      T(N)=-1.D0/LOG(R(N))
      SIG=SIG + GAM(N)*(ND-N)
      ELSE
      IFLAG=N
      ENDIF
10    CONTINUE
C
      IF(N2.GT.N1 .AND. IFLAG.GT.N2)THEN
      TAU=(GAM(N1)/GAM(N2))**(1.D0/(N2-N1))
      SIG=GAM(0) + 2.D0*GAM(N2)/(TAU-1.D0)
      DO 20 N=1,N2
      SIG=SIG + 2.D0*GAM(N)
20    CONTINUE
      TAU=1.D0/LOG(TAU)
      ELSE
      SIG=GAM(0) + 2.D0/ND*SIG
      TAU=0.D0
      ENDIF
      SIG=SQRT(MAX( SIG/ND ,0.))
      RATIO=0.D0
      IF(SIG0.GT.0.D0) RATIO=SIG/SIG0
C
      IF(IUNIT.LE.0)RETURN
C
      WRITE(IUNIT,1010) SIG,(GAM(N),N=0,NMAX)
1010  FORMAT(' SIGMA   =',F12.6/' GAMMA(0)=',F12.6/' GAMMA(N)='/
     X   (10F12.6))
      WRITE(IUNIT,1020) (R(N),N=1,NMAX)
1020  FORMAT(' GAMMA(N)/GAMMA(N-1) ='/(10F12.6))
      WRITE(IUNIT,1030)(T(N),N=1,NMAX)
1030  FORMAT(' RELAXATION TIME T(N) ='/(10F12.4))
      WRITE(IUNIT,1035)N1,N2,TAU
1035  FORMAT(' RELAX. TIME T(',I2,',',I2,') =',F10.4)
      WRITE(IUNIT,1040) RATIO
1040  FORMAT(' SIGMA/SIGMA(UNCORR.) =',F10.4)
90    CONTINUE
      IF(IUNIT.GT.0) WRITE(IUNIT,1050) SIG0
1050  FORMAT(' SIGMA(UNCORRELATED)  =',F12.6)
      END

C********************STAT4*********************************
C                                               B.BUNK 1986
C     STATISTICAL ANALYSIS FOR  NBL  BLOCKS
C        OF  NVAR  CORRELATED VARIABLES EACH
C
C     FOR INTERNAL COUNTERS IN REAL*8:
C        REMOVE C+++++
C
C     CLEAR COUNTERS
C
      SUBROUTINE CLEAR4(NBL1,NVAR1)
*     *****************
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXBL=8,MAXVAR=2+20/2)
      COMMON/ST4/NBL,NVAR,NDAT(MAXBL)
      COMMON/ST4/SAV(MAXBL,MAXVAR),SUM(MAXBL,MAXVAR,MAXVAR)
      SAVE/ST4/
      REAL*8 SAV,SUM
      IF(NBL1.GT.MAXBL)THEN
      WRITE(*,'('' ERROR IN STAT4: NBL ='',I5)')NBL1
      STOP
      ENDIF
      IF(NVAR1.GT.MAXVAR)THEN
      WRITE(*,'('' ERROR IN STAT4: NVAR ='',I5)')NVAR1
      STOP
      ENDIF
      NBL=NBL1
      NVAR=NVAR1
      DO 10 IBL=1,NBL
      NDAT(IBL)=0
      DO 10 I=1,NVAR
      SAV(IBL,I)=0.D0
      DO 10 J=1,NVAR
10    SUM(IBL,I,J)=0.D0
      END
C
C     ACCUMULATE DATA
C
      SUBROUTINE ACCUM4(IBL,F)
*     *****************
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXBL=8,MAXVAR=2+20/2)
      COMMON/ST4/NBL,NVAR,NDAT(MAXBL)
      COMMON/ST4/SAV(MAXBL,MAXVAR),SUM(MAXBL,MAXVAR,MAXVAR)
      SAVE/ST4/
      REAL*8 SAV,SUM,DF
      REAL*8 F(*)
C+++++      DIMENSION F(*)
      IF(IBL.LT.1 .OR. IBL.GT.NBL)RETURN
      NDAT(IBL)=NDAT(IBL)+1
      DO 10 I=1,NVAR
      DF=F(I)
      SAV(IBL,I)=SAV(IBL,I)+DF
      DO 10 J=I,NVAR
10    SUM(IBL,I,J)=SUM(IBL,I,J)+DF*F(J)
      END
C
C     COMPUTE AVERAGES
C
      REAL*8 FUNCTION AVER4(IBL,IVAR)
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXBL=8,MAXVAR=2+20/2)
      COMMON/ST4/NBL,NVAR,NDAT(MAXBL)
      COMMON/ST4/SAV(MAXBL,MAXVAR),SUM(MAXBL,MAXVAR,MAXVAR)
      SAVE/ST4/
      REAL*8 SAV,SUM
      AVER4=0.D0
      IF(IBL.LT.1 .OR. IBL.GT.NBL .OR. IVAR.LT.1 .OR. IVAR.GT.NVAR
     X                            .OR. NDAT(IBL).LT.1)RETURN
      AVER4=SAV(IBL,IVAR)/NDAT(IBL)
      END
C
C     COMPUTE ERRORS
C
      REAL*8 FUNCTION SIGMA4(IBL,IVAR)
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXBL=8,MAXVAR=2+20/2)
      COMMON/ST4/NBL,NVAR,NDAT(MAXBL)
      COMMON/ST4/SAV(MAXBL,MAXVAR),SUM(MAXBL,MAXVAR,MAXVAR)
      SAVE/ST4/
      REAL*8 SAV,SUM,ZERO
      PARAMETER(ZERO=0.D0)
      SIGMA4=0.D0
      IF(IBL.LT.1 .OR. IBL.GT.NBL .OR. IVAR.LT.1 .OR. IVAR.GT.NVAR
     X                            .OR. NDAT(IBL).LT.1)RETURN
      SIGMA4=SQRT(MAX(SUM(IBL,IVAR,IVAR)-SAV(IBL,IVAR)**2/NDAT(IBL)
     X               ,ZERO))/NDAT(IBL)
      END
C
C     COMPUTE COVARIANCE MATRIX
C
      REAL*8 FUNCTION COVAR4(IBL,I,J)
      PARAMETER(MAXVEC=8**3*8+7)
      PARAMETER(MAXBL=8,MAXVAR=2+20/2)
      COMMON/ST4/NBL,NVAR,NDAT(MAXBL)
      COMMON/ST4/SAV(MAXBL,MAXVAR),SUM(MAXBL,MAXVAR,MAXVAR)
      SAVE/ST4/
      REAL*8 SAV,SUM
      COVAR4=0.D0
      IMIN=MIN(I,J)
      IMAX=MAX(I,J)
      ND=NDAT(IBL)
      IF(IBL.LT.1 .OR. IBL.GT.NBL .OR. IMIN.LT.1 .OR. IMAX.GT.NVAR
     X                            .OR. ND.LT.1)RETURN
      COVAR4=((SUM(IBL,IMIN,IMAX)-SAV(IBL,I)*SAV(IBL,J)/ND)/ND)/ND
      END

C****************************TRNGV*************************************
C
C     TRNGV PROVIDES A VECTOR OF RANDOM NUMBERS
C       WITH UNIFORM DISTRIBUTION IN  (0.,RMAX)
C
      SUBROUTINE TRNGV(RANVEC,LVEC,RMAX)
*     ****************
C
      PARAMETER(MAXVEC=8**3*8+7)
      REAL*8 RANVEC,RMAX,R1,RNORM
      PARAMETER(R1=2.D0**(-31))
C
      COMMON/TRNBLK/IREG(0:249+MAXVEC)
      SAVE/TRNBLK/
C
      DIMENSION RANVEC(LVEC)
      DIMENSION ISEED(0:249)
C
      N103=LVEC/103
      LEFT=LVEC - 103*N103
C
C SOME LINES ARE NEW
      IFIRST=0
      IMAX=2147483647
      DO 100 N=1,N103
      DO 10 I=IFIRST,IFIRST+102
C     IREG(I+250)=IEOR(IREG(I+147),IREG(I))
      IREG(I+250)=IREG(I+147)-IREG(I)
      IF (IREG(I+250).LE.0) IREG(I+250)=IREG(I+250)+IMAX
10    CONTINUE
      IFIRST=IFIRST+103
100   CONTINUE
C
C SOME LINES ARE NEW
      DO 20 I=IFIRST,IFIRST+LEFT-1
C     IREG(I+250)=IEOR(IREG(I+147),IREG(I))
      IREG(I+250)=IREG(I+147)-IREG(I)
      IF (IREG(I+250).LE.0) IREG(I+250)=IREG(I+250)+IMAX
20    CONTINUE
C
      RNORM=R1*RMAX
      DO 30 I=1,LVEC
      RANVEC(I)=RNORM*IREG(I+249)
30    CONTINUE
C
      DO 40 I=0,249
      IREG(I)=IREG(I+LVEC)
40    CONTINUE
C
      RETURN
C
C
      ENTRY TRNGET(ISEED)
C
      DO 50 I=0,249
      ISEED(I)=IREG(I)
50    CONTINUE
      RETURN
C
      ENTRY TRNSET(ISEED)
C
      DO 60 I=0,249
      IREG(I)=ISEED(I)
60    CONTINUE
      END
C
      SUBROUTINE INIT(ISTART)
*     ***************
C
      PARAMETER(MAXVEC=8**3*8+7)
C
      REAL*8 U,BETINV
      REAL*8 ZERO,ONE,TWO
      PARAMETER(ZERO=0.,ONE=1.,TWO=2.)
C
      COMMON/PAR1/NDIM,LSIZE(4),LVEC,NDSP
      COMMON/PAR2/NPTS,NLINKS,NPLAQS,NPLPP,NVOL
      COMMON/PAR3/BETA,GAMMA,ITER,IFREQ,MFREQ
      COMMON/INTBLK/BETINV,NERR4,NERR2
      COMMON/CONF/U(MAXVEC,0:3,2**4,4)
C
      ITER=0
      IFREQ=0
      NERR2=0
      NERR4=0
      CALL CLEAR2(2*NDIM+9+8*(LSIZE(NDIM)/2+1),10)
      CALL CLEAR4(8,2+LSIZE(NDIM)/2)
      WRITE(6,'(//'' RANDOM START'')')
      DO 30 I1=1,NDIM
      DO 30 ILAT=1,2**NDIM
      DO 30 I30=1,4
      CALL TRNGV1(LVEC)
      GOTO 30
30    CONTINUE
C
      END
C
      SUBROUTINE TRNGV1(LVEC)
*     ****************
C
      PARAMETER(MAXVEC=8**3*8+7)
      REAL*8 R1
      PARAMETER(R1=2.D0**(-31))
C
      COMMON/TRNBLK/IREG(0:249+MAXVEC)
      SAVE/TRNBLK/
C
      N103=LVEC/103
      LEFT=LVEC - 103*N103
C
C SOME LINES ARE NEW
      IFIRST=0
      IMAX=2147483647
      DO 100 N=1,N103
      DO 10 I=IFIRST,IFIRST+102
C     IREG(I+250)=IEOR(IREG(I+147),IREG(I))
      IREG(I+250)=IREG(I+147)-IREG(I)
      IF (IREG(I+250).LE.0) IREG(I+250)=IREG(I+250)+IMAX
10    CONTINUE
      IFIRST=IFIRST+103
100   CONTINUE
C
C SOME LINES ARE NEW
      DO 20 I=IFIRST,IFIRST+LEFT-1
C     IREG(I+250)=IEOR(IREG(I+147),IREG(I))
      IREG(I+250)=IREG(I+147)-IREG(I)
      IF (IREG(I+250).LE.0) IREG(I+250)=IREG(I+250)+IMAX
20    CONTINUE
C
      DO 40 I=0,249
      IREG(I)=IREG(I+LVEC)
40    CONTINUE
C
      RETURN
      END


program demo_M_escape
          use M_escape, only : esc, esc_mode, update
             write(*,'(a)') esc('TEST MANNER=RAW:')
             call esc_mode(manner='raw')
             write(*,'(a)') esc('<r>RED</r>,<g>GREEN</g>,BLUE')
             write(*,'(a)') esc('<c>CYAN</c>,<m>MAGENTA</g>,<y>YELLOW</y>')
             write(*,'(a)') esc('<w>WHITE</w> and <e>EBONY</e>')

             write(*,'(a)') esc('Adding <bo>bold</bo>')
             write(*,'(a)') esc('<bo><r>RED</r>,<g>GREEN</g>,BLUE</bo>')
             write(*,'(a)') esc('<bo><c>CYAN</c>,<m>MAGENTA</g>,<y>YELLOW</y></bo>')
             write(*,'(a)') esc('<bo><w>WHITE</w> and <e>EBONY</e></bo>')

             write(*,'(a)') esc('Adding
underline
')
             write(*,'(a)') esc('<bo>
<r>RED</r>,<g>GREEN</g>,BLUE
</bo>')
             write(*,'(a)') esc('<bo>
<c>CYAN</c>,<m>MAGENTA</g>,<y>YELLOW</y>
</bo>')
             write(*,'(a)') esc('<bo>
<w>WHITE</w> and <e>EBONY</e>
</bo>')

             write(*,'(a)') esc('Adding
italic
')
             write(*,'(a)') esc('<bo>
<it><r>RED</r>,<g>GREEN</g>,BLUE</it>
</bo>')
             write(*,'(a)') esc('<bo>
<it><c>CYAN</c>,<m>MAGENTA</g>,<y>YELLOW</it></y>
</bo>')
             write(*,'(a)') esc('<bo>
<it><w>WHITE</w> and <e>EBONY</e>
</bo>')

             write(*,'(a)') esc('Adding <in>inverse</in>')
             write(*,'(a)') esc('<in><bo>
<it><r>RED</r>,<g>GREEN</g>,BLUE</it>
</bo></in>')
             write(*,'(a)') esc('<in><bo>
<it><c>CYAN</c>,<m>MAGENTA</g>,<y>YELLOW</it></y>
</bo></in>')
             write(*,'(a)') esc('<in><bo>
<it><w>WHITE</w> and <e>EBONY</e>
</bo></in>')
       end program demo_M_escape


text:
!======================================================================!
  SUBROUTINE ComSav()
!----------------------------------------------------------------------!
! Writes: *.com, convert.scr and CONVERT.scr                           !
! ~~~~~~~                                                              !
!------------------------------[Modules]-------------------------------!
  USE all_mod
  USE gen_mod
  USE par_mod
!----------------------------------------------------------------------!
  IMPLICIT NONE
!-------------------------------[Locals]-------------------------------!
  CHARACTER :: namOut*80
  INTEGER   :: sub
!--------------------------------[CVS]---------------------------------!
!  $Id: ComSav.f90,v 1.4 2002/10/30 16:29:18 niceno Exp $
!  $Source: /home/muhamed/.CVSROOT/T-Rex/Divide/ComSav.f90,v $
!======================================================================!

!<<<<<<<<<<<<<<<<<<<<<<<<<!
!                         !
!     create com file     !
!                         !
!<<<<<<<<<<<<<<<<<<<<<<<<<!
  do sub=1,Nsub
    call NamFil(sub, namOut, '.com', len_trim('.com'))
    open(9, FILE=namOut)
    write(6, *) 'Now creating the file:', namOut
    call NamFil(sub, namOut, '.com', 0)
    write(9,*) namOut
    close(9)
  end do

  open(9, FILE = 'convert.scr')
  do sub=1,Nsub
    call NamFil(sub, namOut, '.com', len_trim('.com'))
    write(9,*) './B2A < ', namOut
  end do
  close(9)

  open(9, FILE = 'CONVERT.SCR')
  do sub=1,Nsub
    call NamFil(sub, namOut, '.com', len_trim('.com'))
    write(9,*) './A2B < ', namOut
  end do
  close(9)

  END SUBROUTINE ComSav

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
PROGRAM KPP_ROOT_Driver

  USE KPP_ROOT_Model
  USE KPP_ROOT_Initialize, ONLY: Initialize

      KPP_REAL :: T, DVAL(NSPEC)
      KPP_REAL :: RSTATE(20)
      INTEGER :: i

!~~~> TIME VARIABLES

      STEPMIN = 0.0d0
      STEPMAX = 0.0d0

      DO i=1,NVAR
        RTOL(i) = 1.0d-4
        ATOL(i) = 1.0d-3
      END DO

      CALL Initialize()
      CALL InitSaveData()

!~~~> Time loop
      T = TSTART
kron: DO WHILE (T < TEND)

        CALL GetMass( C, DVAL )
        WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,       &
                   ( TRIM(SPC_NAMES(MONITOR(i))),           &
                     C(MONITOR(i))/CFACTOR, i=1,NMONITOR ), &
                   ( TRIM(SMASS(i)), DVAL(i)/CFACTOR, i=1,NMASS )
        TIME = T
        CALL SaveData()
        CALL Update_SUN()
        CALL Update_RCONST()

        CALL INTEGRATE( TIN = T, TOUT = T+DT, RSTATUS_U = RSTATE, &
        ICNTRL_U = (/ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 /) )
        T = RSTATE(1)

      END DO kron
!~~~> End Time loop

      CALL GetMass( C, DVAL )
      WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,     &
               ( TRIM(SPC_NAMES(MONITOR(i))),           &
                 C(MONITOR(i))/CFACTOR, i=1,NMONITOR ), &
               ( TRIM(SMASS(i)), DVAL(i)/CFACTOR, i=1,NMASS )
      TIME = T
      CALL SaveData()
      CALL CloseSaveData()

990   FORMAT('Done[%]. Time ',20(4X,A12))
991   FORMAT(F6.1,'%. T=',E9.3,2X,200(A,'=',E11.4,'; '))

END PROGRAM KPP_ROOT_Driver

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
subroutine tara1dHydSer(arg)

use cfgio_mod, only: cfg_t, parse_cfg

implicit none

type(cfg_t):: cfg

real (kind=8), parameter :: pi = 3.14159265358979323846d0

include "fftw3.f"

integer (kind=4) N,Nh

real (kind = 8), dimension (:), allocatable :: x, ux, dux, ddux

complex (kind=8), dimension(:), allocatable :: uk,duk,dduk

integer*8 :: plan
real*8 :: L

integer :: i,t
real*8 dx,k,time,time_min,time_max,dt,nu
real (kind = 8), dimension (:), allocatable :: force_u_n, force_u_o

character (len=90) :: filename
character (len=32) :: arg

cfg = parse_cfg(arg)

call cfg%get("grid","Nx",N)

Nh = (N/2) + 1

allocate(x(N))
allocate(ux(N))
allocate(dux(N))
allocate(ddux(N))

allocate(force_u_n(N))
allocate(force_u_o(N))

allocate(uk(Nh))
allocate(duk(Nh))
allocate(dduk(Nh))

call cfg%get("length","Lx",L)

call cfg%get("resistivity","nu",nu)

call cfg%get("time","time_min",time_min)
call cfg%get("time","time_max",time_max)
call cfg%get("time","dt",dt)

L  = L*pi
dx = L/dfloat(N)

do i = 1,N
  x(i) = dfloat(i-1)*dx
  ux(i) = dsin(x(i))
enddo

do time = time_min, time_max, dt

t = nint(time/dt) - int(time_min/dt)

call dfftw_plan_dft_r2c_1d(plan, N, ux, uk, FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(plan, ux, uk)
call dfftw_destroy_plan(plan)

do i = 1,Nh
  k = 2.0d0*pi*dfloat(i-1)/L
  duk(i) = (0.0d0,1.0d0) * k * uk(i)
  dduk(i) = - k * k * uk(i)
enddo

call dfftw_plan_dft_c2r_1d(plan, N, duk, dux, FFTW_ESTIMATE)
call dfftw_execute_dft_c2r(plan, duk, dux)
call dfftw_destroy_plan(plan)

call dfftw_plan_dft_c2r_1d(plan, N, dduk, ddux, FFTW_ESTIMATE)
call dfftw_execute_dft_c2r(plan, dduk, ddux)
call dfftw_destroy_plan(plan)

do i = 1,N
dux(i) = dux(i)/dfloat(N)
ddux(i) = ddux(i)/dfloat(N)
  if (t .ge. 1000 .and. mod(t,500) == 0) then
    write(filename, '("output/fort.",I8.8)') t
    open(unit=t,file=filename,status='unknown')
    write(t,*) t, x(i), ux(i)
  endif
enddo

do i = 1,N
  force_u_n(i) = -ux(i)*dux(i) + nu*ddux(i)
    if (t==0) then
    ux(i) = ux(i) + dt * force_u_n(i)
    else
    ux(i) = ux(i) + dt* ( (3.0/2.0) * force_u_n(i) - (1.0/2.0) * force_u_o(i))
    endif
 force_u_o(i) = force_u_n(i)
enddo

enddo

end subroutine tara1dHydSer

max_stars_repo_name: RupakMukherjee/Navier-Stokes

max_issues_repo_name: RupakMukherjee/TARA

max_forks_repo_name: RupakMukherjee/TARA


text:
C
C	$Id: pwt.f,v 1.6 2008-07-27 00:59:05 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE PWT(IOS,STATUS)
C
C  START UP THE CHARACTER WRITING SEQUENCE TO THE DEVICE
C
C  OUTPUT
C       IOS-IO STATUS VALID ONLY IF STATUS EQUALS AN I/O ERROR
C       STATUS-THE ERROR STATUS DEFINED BY COMMON TREROR
C
      COMMON/TREROR/ ALLOK, MFRCHK, MTOPER, METRDC, REDERR, TYPCHG
     1             ,INVTYP, MINVLD, TYPERR, FRMEND, ENCINT, IVDCDT
     2             ,GCOERR, GCRERR, GCCERR, FCOERR, FCRERR, FCCERR
     3             ,PLIDXG, PMIDXG, TXIDXG, PGIDXG, INVLMT, CELERR
     4             ,COIERR, COLNRM, UNKNOW, UNKOPC, ENDMTF, VNEROR
     5             ,BADRSZ, DEVOUT, NOVERS, BADFNT, PGMERR, FASERR
     6             ,HINERR, VDWERR, RDWERR, RIXLIM
      INTEGER        ALLOK, MFRCHK, MTOPER, METRDC, REDERR, TYPCHG
     1             ,INVTYP, MINVLD, TYPERR, FRMEND, ENCINT, IVDCDT
     2             ,GCOERR, GCRERR, GCCERR, FCOERR, FCRERR, FCCERR
     3             ,PLIDXG, PMIDXG, TXIDXG, PGIDXG, INVLMT, CELERR
     4             ,COIERR, COLNRM, UNKNOW, UNKOPC, ENDMTF, VNEROR
     5             ,BADRSZ, DEVOUT, NOVERS, BADFNT, PGMERR, FASERR
     6             ,HINERR, VDWERR, RDWERR, RIXLIM
      COMMON /PGKSCOM/ CLPDAT, CLPFLG, POLIDX, LINTYP, LINWTH, LINCOL,
     1                LINRGB, MARIDX, MARSIZ, MARCOL, MARRGB, TXTIDX,
     2                INTSTL, PATIDX, FILCOL, FILRGB, MARTYP, HORIZ ,
     3                VERT  , PATH  , CHIGHT, XU    , YU    , XB    ,
     4                YB    , TXTCOL, FINDEX, CEXPN , CSPACE, CURIDX,
     5                CURIMP, CURINT, COLMOD, FILIDX, TXTRGB, PROPN ,
     6                FIRSTX, FIRSTY, LASTX , LASTY , TRATIO, CBV   ,
     7                CUV   , CHP   , CVP   , CCSPAC, CHHORZ, CDV   ,
     8                PMRSIZ, CLPX  , CLPY  , CLPP  , DEFREP, DEFLEN,
     9                VDCECR, TRANVN, TXTPRE, HATIDX, FILRPT, ASFSRF,
     A                ASFSDF, MAPMOD, VERSOK, PCBFST, CPGLEN, CLPNUL,
     B                MTDLEN
      INTEGER         CMPMAX, ASFMAX
      PARAMETER      (CMPMAX=256    , ASFMAX=18)
      LOGICAL         CLPFLG, PROPN , CHHORZ, DEFREP, CLPNUL, MAPMOD,
     1                VERSOK, PCBFST
      REAL            LINWTH, MARSIZ, CEXPN , CSPACE, TRATIO, CBV(2),
     1                CUV(2), CDV(2), CCSPAC
      INTEGER         CLPDAT(4)     , POLIDX, LINTYP, LINCOL, LINRGB(3),
     1                MARIDX, MARCOL, MARRGB(3)     , TXTIDX, INTSTL,
     2                PATIDX, FILCOL, FILRGB(3)     , MARTYP, HORIZ ,
     3                VERT  , PATH  , XU    , YU    , XB    , YB    ,
     4                TXTCOL, FINDEX, CURIMP(CMPMAX), CURINT(CMPMAX*3),
     5                CHIGHT, COLMOD, FILIDX, TXTRGB(3)     , CURIDX,
     6                FIRSTX, FIRSTY, LASTX , LASTY , CHP   , CVP   ,
     7                PMRSIZ, CLPX  , CLPY  , CLPP  , DEFLEN, TRANVN,
     8                TXTPRE, HATIDX, FILRPT(2)     , VDCECR(4)     ,
     9                ASFSRF(ASFMAX), ASFSDF(ASFMAX), CPGLEN, MTDLEN
      COMMON /GKSCHR/ MTDESC
      CHARACTER*80    MTDESC
      COMMON /TRMTYP/ FOURBT, MULTBT, TWOBT, PTSFLG, CURCOL, IPWRIT,
     1                IFRAME, ISETOP, IRANG, IPTS, ICAS, IINT, IORNT,
     2                ISIZE, IFONT, IPAT, ISSIZ, ICENT, ICOLR,
     3                MCENT, MSIZE, MOR, OLDX, OLDY, MCASE, MINTEN,
     4                MFONT, MDASH, MSPOT, MCOLOR, MLIMIT, PTSYM,
     5                SAVX, SAVY
      INTEGER FOURBT, MULTBT, TWOBT, PTSFLG, CURCOL, IPWRIT,
     1                IFRAME, ISETOP, IRANG, IPTS, ICAS, IINT, IORNT,
     2                ISIZE, IFONT, IPAT, ISSIZ, ICENT, ICOLR,
     3                MCENT, MSIZE, MOR, OLDX, OLDY, MCASE, MINTEN,
     4                MFONT, MDASH, MSPOT, MCOLOR(3), MLIMIT(4), PTSYM,
     5                SAVX, SAVY
      COMMON /TRTYPE/ METMIN, METEXT, METPRT, METHED, MIOPCL, MXOPCL,
     1                MIOPID, MXOPID, ASCDEC, ASCHEX, ASCOCT, BINARY,
     2                GAHNOR, GALEFT, GACENT, GARITE, GAVNOR, GATOP ,
     3                GACAP , GAHALF, GABASE, GABOTT, GRIGHT, GLEFT ,
     4                GUP   ,
     5                GDOWN , CINDEX, CDIRCT, PENUP , PENDN , FLAGVL,
     6                SEPRVL, ASCTEK, TYPRGB, TYPBGR,
     7                TYPHLS, TYPMON, SOLID , DASHED, INVSBL, SOLPAT,
     8                ASFIND, ASFBND, ASFLNT, ASFLNW, ASFLNC, ASFMRT,
     9                ASFMRS, ASFMRC, ASFFIS, ASFFHI, ASFFPI, ASFFCO,
     A                ASFFPT, ASFFPW, ASFFPC, ASFTFI, ASFTPR, ASFTCX,
     B                ASFTCS, ASFTCO, ASCFLT, PHOLLO, PSOLID, PPATTR,
     C                PHATCH, PEMPTY, HHORIZ, HVERTI, HPOSLP, HNESLP,
     D                HHOAVE, HPOANE
      INTEGER         METMIN, METEXT, METPRT, METHED, MIOPCL, MXOPCL,
     1                MIOPID, MXOPID, ASCDEC, ASCHEX, ASCOCT, BINARY,
     2                GAHNOR, GALEFT, GACENT, GARITE, GAVNOR, GATOP ,
     3                GACAP , GAHALF, GABASE, GABOTT, GRIGHT, GLEFT ,
     4                GUP   ,
     5                GDOWN , CINDEX, CDIRCT, PENUP , PENDN , FLAGVL,
     6                SEPRVL, ASCTEK, TYPRGB, TYPBGR,
     7                TYPHLS, TYPMON, SOLID , DASHED, INVSBL, SOLPAT,
     8                ASFIND, ASFBND, ASFLNT, ASFLNW, ASFLNC, ASFMRT,
     9                ASFMRS, ASFMRC, ASFFIS, ASFFHI, ASFFPI, ASFFCO,
     A                ASFFPT, ASFFPW, ASFFPC, ASFTFI, ASFTPR, ASFTCX,
     B                ASFTCS, ASFTCO, ASCFLT, PHOLLO, PSOLID, PPATTR,
     C                PHATCH, PEMPTY, HHORIZ, HVERTI, HPOSLP, HNESLP,
     D                HHOAVE, HPOANE
      COMMON /TRINST/ OPCL, OPID, LEN, CONT
      INTEGER OPCL, OPID, LEN
      LOGICAL CONT
      COMMON /CAPFNT/ CHRSTR, CHREND, CHRWDT, CHRSCH, FNTRHT,
     1                FNTTOP, FNTCAP, FNTHLF, FNTBAS, FNTBOT,
     2                FNTTYP, FNTSCH, CORXST, CORXLN, CORYST,
     3                CORYLN, CORPST, CORPLN, PBEGST, PBEGLN,
     4                PENDST, PENDLN, CORSCH, NEWCLS, CPNTRS,
     5                CPNLST, CSTRKS
      INTEGER         CHRSM1, CHRSM2
      PARAMETER       (CHRSM1=128, CHRSM2=5121)
      INTEGER         CHRSTR, CHREND, CHRWDT,
     1                CHRSCH(11)    , FNTRHT, FNTTOP,
     2                FNTCAP, FNTHLF, FNTBAS, FNTBOT, FNTTYP,
     3                FNTSCH(9)     , CORXST, CORXLN,
     4                CORYST, CORYLN, CORPST, CORPLN,
     5                PBEGST, PBEGLN, PENDST, PENDLN,
     6                CORSCH(10)    ,
     7                NEWCLS(300)   ,
     8                CPNTRS(CHRSM1), CPNLST        ,
     9                CSTRKS(CHRSM2)
C
C  XLCN   --  Bit offset to extract X coordinates from fontcap.
C  YLCN   --  Bit offset to extract Y coordinates from fontcap.
C  PLCN   --  Bit offset to extract PEN control bit from fontcap.
C  FNTMAX --  Number of fonts obtained from FONT LIST (defaults to 1).
C  FNTS   --  An array of pointers associating indices obtained from
C             FONT INDEX elements with the internally supported fonts.
C             Consult the documentation in BLOCKDATA TRNDATX for a
C             list of supported fonts.
C  FLEFT  --  An array of character left values for proportionally
C             spaced fonts.
C  FRIGHT --  An array of character right values for proportionally
C             spaced fonts.
C  MXFLST --  The maximum number of strings accepted in a FONT LIST
C             element.
C  NSFNTS --  The number of distinct fonts supported internally.
C
      PARAMETER (NSFNTS=20,MXFLST=300)
      COMMON /FNTCOM/ XLCN  , YLCN  , PLCN  , FNTMAX,
     1                FNTS  , FLEFT , FRIGHT, SFLENS
      INTEGER         FLEFT(CHRSM1) , FRIGHT(CHRSM1), FNTS(MXFLST)  ,
     1                XLCN, YLCN    , PLCN  , FNTMAX
      CHARACTER*30 SFONTS(NSFNTS)
      COMMON /SUPFNT/ SFONTS
C
C  Size of common CAPFNT in words.
C
      INTEGER WRLEN
      PARAMETER (WRLEN=5600)
      COMMON /CAPLIN/ PLAVBL, LDSSTR, LDSSIZ, LDTSTR, LDTSIZ,
     1                LMSSTR, LMSSIZ, LMTSTR, LMTSIZ, LCSSTR,
     2                LCSSIZ, LCTSTR, LCTSIZ, LINFIN, LINFMT,
     3                LWSSTR, LWSSIZ, LWTSTR, LWTSIZ, LWTFIN,
     4                LWTFMT, LWTRNG, LWTSCF, LBSSTR, LBSSIZ,
     5                LBTSTR, LBTSIZ, LPSSTR, LPSSIZ, LPTSTR,
     6                LPTSIZ
      INTEGER         LDSMAX, LDTMAX, LMSMAX, LMTMAX, LCSMAX,
     1                LCTMAX, LVCFMX, LWSMAX, LWTMAX, LWTFMX,
     2                LBSMAX, LBTMAX, LPSMAX, LPTMAX
      PARAMETER   (LDSMAX=30, LDTMAX=15, LMSMAX=30, LMTMAX=15,
     1             LCSMAX=30, LCTMAX=15, LVCFMX=8,  LWSMAX=30,
     2             LWTMAX=15, LWTFMX=8 , LBSMAX=30, LBTMAX=15,
     3             LPSMAX=20, LPTMAX=20)
      INTEGER         LDSSTR(LDSMAX), LDSSIZ, LDTSTR(LDTMAX),
     1                LDTSIZ, LMSSTR(LMSMAX), LMSSIZ,
     2                LMTSTR(LMTMAX), LMTSIZ, LCSSTR(LCSMAX),
     3                LCSSIZ, LCTSTR(LCTMAX), LCTSIZ, LINFIN(8),
     4                LINFMT(LVCFMX,4)      , LWSSTR(LWSMAX),
     5                LWSSIZ, LWTSTR(LWTMAX), LWTSIZ, LWTFIN(8),
     6                LWTFMT(LWTFMX,4)      , LWTRNG(2)     ,
     7                LBSSTR(LBSMAX), LBSSIZ, LBTSTR(LBTMAX),
     8                LBTSIZ, LPSSTR(LPSMAX), LPSSIZ,
     9                LPTSTR(LPTMAX), LPTSIZ
      LOGICAL         PLAVBL
      REAL            LWTSCF,LINRIN(8),LWTRIN(8)
      INTEGER         LENLIN
      PARAMETER   (LENLIN=1+LDSMAX+1+LDTMAX+1+LMSMAX+1+
     1             LMTMAX+1+LCSMAX+1+LCTMAX+1+8+LVCFMX*4+
     2             LWSMAX+1+LWTMAX+1+8+LWTFMX*4+2+1+LBSMAX+
     3             1+LBTMAX+1+LPSMAX+1+LPTMAX+1)
      EQUIVALENCE (LINFIN,LINRIN), (LWTFIN,LWTRIN)
C
      INTEGER IOS, STATUS
      INTEGER K, CCOUNT, CBYTE, NULL, ICNT, STRX, STRY
      CHARACTER*300 TCHAR
      REAL DEGRAD,METCOR
C
C  256 CHARACTERS MAX ALLOWED
C
      INTEGER MXCH
      PARAMETER (MXCH=300)
      INTEGER ACHAR(MXCH)
      DATA NULL /0/
C
C  DEGREES TO RADIANS CONVERSION AND COORDINATE EXPANDER
C
      DATA DEGRAD,METCOR/.01745329,32767.0/
C
C  DEFINE INSTRUCTION COUNT DATA AND CHARACTER BYTE LENGTH
C
      DATA ICNT,CBYTE/8,8/
C
      STATUS = ALLOK
C
C  GET THE CHARACTER COUNT (NO ERROR CHECKING REQUIRED)
C
      CALL MNINST(ICNT,IOS,STATUS)
      CCOUNT = OPCL
      IF (CCOUNT.GT.MXCH) CCOUNT = MXCH
C
C  GET THE CHARACTER STRING
C
      DO 100 K=1,CCOUNT
        CALL MNINST(CBYTE,IOS,STATUS)
        ACHAR(K) = OPCL
        IF (STATUS.NE.ALLOK) RETURN
C
 100  CONTINUE
C
C  TEST IF ANOTHER READ NEEDED TO GET TO A 16 BIT BOUNDARY
C
      IF (0 .NE. MOD(CCOUNT,2)) CALL MNINST(CBYTE,IOS,STATUS)
C
C  GET THE COORDINATE LOCATION FOR THE STRING
C       DO NOT PASS THE COMMON VALUES BECAUSE THEY ARE RESET BY PLTIT
C
      STRX = OLDX
      STRY = OLDY
C
      CALL GSTXFP(MFONT,2)
      TCHAR = ' '
      DO 130 I=1,CCOUNT
      TCHAR(I:I) = CHAR(ACHAR(I))
  130 CONTINUE
      CALL PLOTIT(0,0,2)
      CALL WTSTR(FLOAT(STRX)/32767.,FLOAT(STRY)/32767.,
     -           TCHAR(1:CCOUNT),MSIZE/32,MOR,MCENT-1)
C
      RETURN
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE IFP1S (LIST,ISTOR,NLIST)
C
C     THIS ROUTINE FINDS ANY OVERLAPPING INTERVALS IN A SET LIST.
C     IT WILL ALSO CHECK SINGLES
C
      INTEGER         OTPE
      DIMENSION       LIST(1),ISTOR(1)
      CHARACTER       UFM*23,UWM*25
      COMMON /XMSSG / UFM,UWM
      COMMON /SYSTEM/ SYSBUF,OTPE,INX(6),NLPP,INX1(2),LINE
C
      IPAIR = 0
      DO 100 I = 1,NLIST
      IF (LIST(I) .GT. 0) GO TO 100
   10 IF (IPAIR   .NE. 0) GO TO 40
   20 L = 2*IPAIR + 1
      ISTOR(L  ) = LIST(I-1)
      ISTOR(L+1) = LIST(I  )
      IPAIR = IPAIR + 1
   30 LIST(I  ) = 0
      LIST(I-1) = 0
      GO TO 100
C
C     PAIR FOUND - CHECK FOR OVERLAP
C
   40 L  = 1
      IN = IABS(LIST(I-1))
      IOUT = IABS(LIST(I  ))
      K  = 2*L - 1
   50 IF (IN.GE.ISTOR(K)   .AND. IN.LE.IABS(ISTOR(K+1))  ) GO TO 60
      IF (IOUT.GE.ISTOR(K) .AND. IOUT.LE.IABS(ISTOR(K+1))) GO TO 60
      L  = L + 1
      IF (L .LE. IPAIR) GO TO 50
C
C     STORE NEW PAIR IN LIST
C
      GO TO 20
C
C     ERROR IN INTERVAL
C
   60 LIST(I-1) = MIN0(IN,ISTOR(K))
      LIST(I  ) =-MAX0(IOUT,IABS(ISTOR(K+1)))
      IF (LIST(I-1).EQ.ISTOR(K) .AND. LIST(I).EQ.ISTOR(K+1)) GO TO 30
      IX = IABS(ISTOR(K+1))
      WRITE  (OTPE,70) UWM,IN,IOUT,ISTOR(K),IX
   70 FORMAT (A25,' 621, INTERVAL',I8,' THRU',I8,' OVERLAPS INTERVAL',
     1       I8,' THRU', I8,'. THE MAXIMUM INTERVAL WILL BE USED.')
      LINE = LINE + 3
      IF (LINE .GE. NLPP) CALL PAGE
C
C     REMOVE PAIR L FROM ISTOR
C
   80 IF (L .GE. IPAIR) GO TO 90
      M = 2*L + 1
      K = 2*L - 1
      ISTOR(K  ) = ISTOR(M  )
      ISTOR(K+1) = ISTOR(M+1)
      L = L + 1
      GO TO 80
   90 IPAIR = IPAIR - 1
      GO TO 10
  100 CONTINUE
C
C     ALL PAIRS PROCESSED - TRY SINGLES
C
      ISING = 0
      M  = 2*IPAIR
      DO 180 I = 1,NLIST
      IN = LIST(I)
      IF (IPAIR   .EQ. 0) GO TO 140
      IF (LIST(I) .EQ. 0) GO TO 180
C
C     CHECK EACH PAIR
C
      L = 1
  110 K = 2*L - 1
      IF (IN.GE.ISTOR(K) .AND. IN.LE.IABS(ISTOR(K+1))) GO TO 120
      L = L + 1
      IF (L-IPAIR) 110,110,140
C
C     ERROR -- PAIR CONTAINS SINGLE
C
  120 IN1 = IABS(ISTOR(K+1))
      WRITE  (OTPE,130) UWM,IN,ISTOR(K),IN1
  130 FORMAT (A25,' 619, SET MEMBER',I8,' BELONGS TO',I8,' THRU',I8)
      LINE = LINE + 3
      IF (LINE .GE. NLPP) CALL PAGE
      GO TO 180
C
C     CHECK FOR DUPLICATE SINGLES
C
  140 IF (ISING .EQ. 0) GO TO 170
      DO 160 K = 1,ISING
      L = 2*IPAIR + K
      IF (IN .NE. ISTOR(L)) GO TO 160
      WRITE  (OTPE,150) UWM,IN
  150 FORMAT (A25,' 620, SET MEMBER',I8,' IS DUPLICATED IN SET LIST.')
      LINE = LINE + 3
      IF (LINE .GE. NLPP) CALL PAGE
      GO TO 180
  160 CONTINUE
  170 M = M + 1
      ISING = ISING + 1
      ISTOR(M) = IN
  180 CONTINUE
C
C     COPY GOOD STUFF INTO LIST
C
      DO 190 I = 1,M
  190 LIST(I) = ISTOR(I)
      NLIST = M
C
C     SORT LIST
C
      N1 = M - 1
      DO 230 I = 1,N1
      N2 = I + 1
      DO 220 K = N2,M
      IF (IABS(LIST(I))-IABS(LIST(K))) 220,220,210
C
C     SWITCH
C
  210 IN = LIST(I)
      LIST(I) = LIST(K)
      LIST(K) = IN
  220 CONTINUE
  230 CONTINUE
C
      RETURN
      END


text:
module test
contains

end module test

program main
    use plantFEM
    use test
    implicit none

    real(real64),parameter :: duration=10.0d0
    integer(int32),parameter :: sampling = 4096
    type(Random_) :: random
    type(IO_) :: f,z
    type(Math_) :: Math
    integer(int32) :: i
    complex(complex64 ),allocatable :: wave(:),spectre(:)
    real(real64),allocatable :: omega(:), t(:)
    real(real64) :: dt,set_omega,set_omega_sub

    set_omega = (100.0d0/2.0d0/math%PI)
    set_omega_sub = (1000.0d0/2.0d0/math%PI)

    print *, "w1= ", set_omega
    print *, "w2= ", set_omega_sub

    create_input_wave: block
        t    = linspace([0.0d0, duration], sampling)
        dt   = duration/dble(sampling)
        wave = zeros(4096)
        do i=1, size(wave)
            wave(i) = random%gauss(mu=0.0d0, sigma=1.0d0)
            wave(i) = sin( set_omega*t(i) )+sin( set_omega_sub*t(i) )+random%gauss(mu=0.0d0, sigma=1.0d0)
        enddo
        wave = taper(wave,margin=0.050d0)
    end block create_input_wave

    show_input_wave: block
        call f%open("input_wave.txt")
        call f%write( t,real(wave) )
        call f%close()
        call f%plot(option="with lines")
    end block show_input_wave

    run_FFT: block
        omega   = linspace([0.0d0, dble(sampling)/duration*2.0d0*math%PI ], sampling)
        spectre = zeros(sampling)
        spectre =  FFT(wave)
        call z%open("spectre.txt")
        do i=1,size(spectre)/2
            call z%write(omega(i),abs( spectre(i)))
        enddo
        call z%close()
        call z%plot(option="with lines",logscale=true)
    end block run_FFT

end program main

max_stars_repo_name: kazulagi/faba


text:
module external_force_module

  use bl_types
  use bl_constants_module
  use multifab_module
  use ml_layout_module
  use multifab_fill_ghost_module
  use probin_multispecies_module, only: Dbar, sigma, beta, alpha1, delta
  use probin_common_module, only: prob_lo, prob_hi, prob_type, molmass

  implicit none

  private

  public :: external_source

contains

  subroutine external_source(mla,rho,fluxdiv,dx,time)

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(in   ) :: rho(:)
    type(multifab) , intent(inout) :: fluxdiv(:)
    real(kind=dp_t), intent(in   ) :: dx(:,:),time

    ! local variables
    integer :: i,n,ng,dm,nlevs
    integer :: lo(mla%dim),hi(mla%dim)

    ! pointer for rho,fluxdiv
    real(kind=dp_t), pointer :: dp(:,:,:,:)
    real(kind=dp_t), pointer :: dp1(:,:,:,:)

    type(bl_prof_timer), save :: bpt

    call build(bpt,"external_source")

    dm    = mla%dim         ! dimensionality
    ng    = fluxdiv(1)%ng   ! number of ghost cells
    nlevs = mla%nlevel      ! number of levels

    do n=1,nlevs
       do i = 1, nfabs(rho(n))
          dp  => dataptr(rho(n),i)
          dp1 => dataptr(fluxdiv(n),i)
          lo = lwb(get_box(rho(n),i))
          hi = upb(get_box(rho(n),i))
          select case (dm)
          case (2)
             call external_source_2d(dp(:,:,1,:),dp1(:,:,1,:),&
                                     ng,lo,hi,dx(n,:),time)
          case (3)
             call external_source_3d(dp(:,:,:,:),dp1(:,:,:,:),&
                                     ng,lo,hi,dx(n,:),time)
          end select
       end do
    end do

    call destroy(bpt)

  end subroutine external_source

  subroutine external_source_2d(rho,fluxdiv,ng,lo,hi,dx,time)

    integer          :: lo(:),hi(:),ng
    real(kind=dp_t)  :: rho(lo(1)-ng:,lo(2)-ng:,:)
    real(kind=dp_t)  :: fluxdiv(lo(1)-ng:,lo(2)-ng:,:)
    real(kind=dp_t)  :: dx(:),time

    ! local variables
    integer          :: i,j
    real(kind=dp_t)  :: x,y,r,L(2),D12,D23,r_temp,r_temp1

    L(1:2) = prob_hi(1:2)-prob_lo(1:2) ! Domain length
    D12    = Dbar(1)
    D23    = Dbar(3)

    select case (abs(prob_type))

    !==== for m1 = m2 = m and also for m1 != m2 != m, 2 species ====!
    case(4)

     ! for specific box, now start loops over alloted cells
     do j=lo(2), hi(2)
        y = prob_lo(2) + (dble(j)+0.5d0) * dx(2) - 0.5d0
        do i=lo(1), hi(1)
           x = prob_lo(1) + (dble(i)+0.5d0) * dx(1) - 0.5d0

             r = sqrt((x-L(1)*0.5d0)**2 + (y-L(2)*0.5d0)**2)

             r_temp = -(dexp(-(beta*time) - (r**2*(1.0d0 + time))/(4.0d0*D12*time))*(16.0d0*beta*D12**2*&
                      dexp(r**2/(4.0d0*D12))*M_PI*time + alpha1*(r**2 + 4.0d0*beta*D12*time)))/(64.0d0*&
                      D12**3*M_PI**2*time**2)

             fluxdiv(i,j,1) = fluxdiv(i,j,1) + r_temp
             fluxdiv(i,j,2) = fluxdiv(i,j,2) - r_temp

        end do
     end do

    !==== for m2=m3, D12=D13 and Grad(w2)=0, 3 species ====!
    case(5)

     ! for specific box, now start loops over alloted cells
     do j=lo(2), hi(2)
        y = prob_lo(2) + (dble(j)+0.5d0) * dx(2) - 0.5d0
        do i=lo(1), hi(1)
           x = prob_lo(1) + (dble(i)+0.5d0) * dx(1) - 0.5d0

             r = sqrt((x-L(1)*0.5d0)**2 + (y-L(2)*0.5d0)**2)
             r_temp  = (alpha1*(2.0d0*alpha1*(D12*molmass(1) - D23*molmass(2))*(r - sigma)*(r +&
                       sigma) - D12*dexp(r**2/(2.0d0*sigma**2))*molmass(1)*(r**2 - 2.0d0*sigma**2)))/&
                       (dexp(r**2/sigma**2)*molmass(1)*sigma**4)

             r_temp1 = (dexp(-r**2/(2.0d0*sigma**2) - beta*time)*(-(beta*D12*delta*dexp(r**2/(2.0d0*&
                       sigma**2))*molmass(1)*sigma**4) + alpha1*delta*(D12*(D12 - D23)*molmass(1)*&
                       r**2 + 2.0d0*D12*(-D12 + D23)*molmass(1)*sigma**2 + beta*(D12*molmass(1) -&
                       D23*molmass(2))*sigma**4)))/(D12*molmass(1)*sigma**4)

             fluxdiv(i,j,1) = fluxdiv(i,j,1) + r_temp
             fluxdiv(i,j,2) = fluxdiv(i,j,2) + r_temp1
             fluxdiv(i,j,3) = fluxdiv(i,j,3) - (r_temp + r_temp1)

             !== to benchmark eqn1 ==!
             !r_temp = (alpha1*D12*(-(dexp(r**2/(2.0d0*sigma**2))*(r**2 - 2.0d0*sigma**2)) +&
             !         2.0d0*beta*(-r**2 + sigma**2)))/(dexp(r**2/sigma**2)*sigma**4)
             !fluxdiv(i,j,1) = fluxdiv(i,j,1) + r_temp
             !fluxdiv(i,j,2) = fluxdiv(i,j,2) - r_temp
             !=======================!

        end do
     end do

    end select

  end subroutine external_source_2d

  subroutine external_source_3d(rho,fluxdiv,ng,lo,hi,dx,time)

    integer          :: lo(:),hi(:),ng
    real(kind=dp_t)  :: rho(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:,:)
    real(kind=dp_t)  :: fluxdiv(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:,:)
    real(kind=dp_t)  :: dx(:),time

    ! local variables
    integer          :: i,j,k
    real(kind=dp_t)  :: x,y,z,r,L(3),D12,D23,r_temp,r_temp1

    L(1:3) = prob_hi(1:3)-prob_lo(1:3) ! Domain length
    D12    = Dbar(1)
    D23    = Dbar(3)

    select case (abs(prob_type))

    !==== for m1 = m2 = m and m1 != m2 != m, 2 species ====!
    case(4)

     ! for specific box, now start loops over alloted cells
     do k=lo(3), hi(3)
        z = prob_lo(3) + (dble(k)+0.5d0) * dx(3) - 0.5d0
        do j=lo(2), hi(2)
           y = prob_lo(2) + (dble(j)+0.5d0) * dx(2) - 0.5d0
           do i=lo(1), hi(1)
              x = prob_lo(1) + (dble(i)+0.5d0) * dx(1) - 0.5d0

              r = sqrt((x-L(1)*0.5d0)**2 + (y-L(2)*0.5d0)**2 + (z-L(3)*0.5d0)**2)

              r_temp = (dexp(-(beta*time) - (r**2*(1.0d0 + time))/(4.0d0*D12*time))*&
                       (0.02244839026564582d0*beta*D12**3.0d0*dexp(r**2/(4.0d0*D12))*(D12*&
                       time)**3.5d0 + alpha1*D12**1.5d0*(0.0001259825563796855d0*r**2*(D12*time)**2.5d0 +&
                       0.000503930225518742d0*beta*(D12*time)**3.5d0)))/(D12**3.0d0*(D12*time)**5.0d0)

              fluxdiv(i,j,k,1) = fluxdiv(i,j,k,1) - r_temp
              fluxdiv(i,j,k,2) = fluxdiv(i,j,k,2) + r_temp

           end do
        end do
     end do

    !==== for m2=m3, D12=D13 and Grad(w2)=0, 3 species ====!
    case(5)

     do k=lo(3), hi(3)
        z = prob_lo(3) + (dble(k)+0.5d0) * dx(3) - 0.5d0
        do j=lo(2), hi(2)
           y = prob_lo(2) + (dble(j)+0.5d0) * dx(2) - 0.5d0
           do i=lo(1), hi(1)
              x = prob_lo(1) + (dble(i)+0.5d0) * dx(1) - 0.5d0

              r = sqrt((x-L(1)*0.5d0)**2 + (y-L(2)*0.5d0)**2 + (z-L(3)*0.5d0)**2)

               r_temp = (alpha1*(-(D12*dexp(r**2/(2.0d0*sigma**2))*molmass(1)*(r**2 - 3.0d0*sigma**2)) +&
                        alpha1*(D12*molmass(1) - D23*molmass(2))*(2.0d0*r**2 - 3.0d0*sigma**2)))/&
                        (dexp(r**2/sigma**2)*molmass(1)*sigma**4)

              r_temp1 = (dexp(-r**2/(2.0d0*sigma**2) - beta*time)*(-(beta*D12*delta*dexp(r**2/(2.0d0*&
                        sigma**2))*molmass(1)*sigma**4) + alpha1*delta*(D12*(D12 - D23)*molmass(1)*&
                        r**2 + 3.0d0*D12*(-D12 + D23)*molmass(1)*sigma**2 + beta*(D12*molmass(1) -&
                        D23*molmass(2))*sigma**4)))/(D12*molmass(1)*sigma**4)

              fluxdiv(i,j,k,1) = fluxdiv(i,j,k,1) + r_temp
              fluxdiv(i,j,k,2) = fluxdiv(i,j,k,2) + r_temp1
              fluxdiv(i,j,k,3) = fluxdiv(i,j,k,3) - (r_temp + r_temp1)

           end do
        end do
     end do

    end select

  end subroutine external_source_3d

end module external_force_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures