PROGRAM latgas
*     simulate lattice gas model with periodic boundary conditions
*     must be modifed slightly to include barrier
*     graphics routines not listed
      IMPLICIT NONE
      INTEGER lat(10000),nn(10000,0:5)
      INTEGER Lx,Ly,nstep,nplot,istep
      INTEGER rule(0:1024)
      CALL initial(lat,Lx,Ly,nstep,nplot)
*     CALL initial_graphics(Lx,Ly)
      CALL nntable(nn,Lx,Ly)
      CALL ruletable(rule)
*     CALL plot(lat,Lx,Ly)
      DO istep = 1,nstep
         CALL update(lat,rule,nn,Lx,Ly)
*        IF (MOD(istep,nplot) .EQ. 0) THEN
*             CALL plot(lat,Lx,Ly)
*        END IF
      END DO
*     CALL end_graphics
      END

      SUBROUTINE initial(lat,Lx,Ly,nstep,nplot)
      IMPLICIT NONE
      INTEGER lat(10000)
      INTEGER Lx,Ly,nstep,nplot
      INTEGER i,j,n
      Lx = 50
      Ly = 20
      nstep = 100
      nplot = 1
*     begin with no particles
      DO n = 1,Lx*Ly
         lat(n) = 0
      END DO
*     fill block in center of lattice with 6 particles per site
      DO j = Ly/2-1,Ly/2+1
         DO i = Lx/2-1,Lx/2+1
            lat((j-1)*Lx + i) = 63
         END DO
      END DO
      END

      SUBROUTINE update(lat,rule,nn,Lx,Ly)
      IMPLICIT NONE
      INTEGER Lx,Ly,i,j,n,dir,barrier
      PARAMETER (barrier = 7)
      INTEGER latn(10000),lat(10000),nv(0:5)
      INTEGER rule(0:1024),nn(10000,0:5)
*     bounce back boundary conditions v goes to -v
      DATA nv/3,4,5,0,1,2/
      DATA latn/10000*0/
      save latn
*     move particles
      DO j = 1,Ly
         DO i = 1,Lx
            n = (j-1)*Lx + i
               DO dir = 0,5
               IF (BTEST(lat(nn(n,dir)),dir)) THEN
                  IF (BTEST(lat(n),barrier)) THEN
*                    reflection
                     latn(nn(n,dir)) = IBSET(latn(nn(n,dir)),nv(dir))
                  else
*                    particle moves from nearest neighbor
                     latn(n) = IBSET(latn(n),dir)
                  END IF
               END IF
               END DO
*              WRITE(9,*)n,lat(n),latn(n)
            END DO
      END DO
*     collisions
      DO i = 1,Lx
         DO j = 1,Ly
            n = (j-1)*Lx + i
            IF (.NOT.BTEST(lat(n),barrier)) THEN
               lat(n) = rule(latn(n))
               latn(n) = 0
            END IF
         END DO
       END DO
       END

      SUBROUTINE ruletable(rule)
      IMPLICIT NONE
*     6-bit saturated deterministic rule
      INTEGER rule(0:1024)
      INTEGER i
      DO i = 0,256
          rule(i) = i
      END DO
      rule(21) = 42
      rule(42) = 21
      rule(9) = 36
      rule(18) = 9
      rule(36) = 18
      rule(27) = 45
      rule(45) = 54
      rule(54) = 27
      rule(19) = 37
      rule(37) = 19
      rule(50) = 41
      rule(41) = 50
      rule(22) = 13
      rule(13) = 22
      rule(26) = 44
      rule(44) = 26
      rule(11) = 38
      rule(38) = 11
      rule(25) = 52
      rule(52) = 25
      END

      SUBROUTINE nntable(nn,Lx,Ly)
*     nn(n,dir) gives the neighbor whose particle moving in
*     direction dir will move to site n.
      IMPLICIT NONE
      INTEGER nn(10000,0:5)
      INTEGER Lx,Ly,i,j,ip,im,jp,jm,n
*     odd values of y
      DO j = 1,Ly,2
         DO i = 1,Lx
            n = (j-1)*Lx + i
            ip = i + 1
            IF (ip .GT. Lx) ip = 1
            im = i - 1
            IF (im .LT. 1) im = Lx
            jp = j + 1
            IF (jp .GT. Ly) jp = 1
            jm = j - 1
            IF (jm .LT. 1) jm = Ly
            nn(n,0) = (j - 1)*Lx + im
            nn(n,1) = (jp - 1)*Lx + im
            nn(n,2) = (jp - 1)*Lx + i
            nn(n,3) = (j - 1)*Lx + ip
            nn(n,4) = (jm - 1)*Lx + i
            nn(n,5) = (jm - 1)*Lx + im
         END DO
      END DO
*     even values of y
      DO j = 2,Ly,2
         DO i = 1,Lx
            n = (j - 1)*Lx + i
            ip = i + 1
            IF (ip .GT. Lx) ip = 1
            im = i - 1
            IF (im .LT. 1) im = Lx
            jp = j + 1
            IF (jp .GT. Ly) jp = 1
            jm = j - 1
            IF (jm .LT. 1) jm = Ly
            nn(n,0) = (j - 1)*Lx + im
            nn(n,1) = (jp - 1)*Lx + i
            nn(n,2) = (jp - 1)*Lx + ip
            nn(n,3) = (j - 1)*Lx + ip
            nn(n,4) = (jm - 1)*Lx + ip
            nn(n,5) = (jm - 1)*Lx + i
         END DO
      END DO
      END