PROGRAM lattice_gas
*     simulation of particle diffusion in a lattice gas
      IMPLICIT NONE
      REAL*8 L_2, site(40,40), t, tmax, x0(1200), y0(1200)
      INTEGER L, N, im1, IC1_, im2, IC2_, IR_, x(1200), y(1200)
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL initial(L,L_2,N,tmax,im2,IC2_)
      CALL lattice(site,x,y,x0,y0,L,N,im1,IC1_)
      DO t = 1, tmax
         CALL move(site,x,y,x0,y0,L,N,im1,IC1_)
         CALL mydata(x,y,x0,y0,L,L_2,N,t,im2,IC2_)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(L,L_2,N,tmax,im2,IC2_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 L_2, tmax
      REAL GWaspect
      REAL dmy, rnd
      INTEGER IR_, L, N, im2, IC2_
      dmy = rnd(-1)
      L = 40                                ! linear dimension of lattice
      L_2 = 0.5D0*L
      N = 500                               ! number of particles
      tmax = 50                     ! # Monte Carlo steps per particle
      CALL GWcolor(IR_, 0, 2)
      CALL GWvport(IR_, 0.52*GWaspect(1), 0.0, GWaspect(1), 1.0)
      CALL GWindow(IR_, -1.0, -1.0, REAL(tmax + 1), REAL(L_2 + 1))
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, N + 1)
      CALL GWcolor(IR_, 19, 7)
      CALL GWsetpen(IR_, 19, -1, -1, 4)
      WRITE(Stmp1_,*) 'Plot of <R^2(t)> versus t'
      CALL GWputtxt(IR_, 1.0, REAL(L_2), Stmp1_)
      CALL GWrect(IR_, 0.0, 0.0, REAL(tmax), REAL(L_2))
      CALL GWcmbmrk(im2, 0, 1, REAL(L_2/tmax/4), 19, -1, 4)
      END
C
      SUBROUTINE lattice(site,x,y,x0,y0,L,N,im1,IC1_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 site(40,40), x0(1200), y0(1200)
      REAL r, xwin, ywin, GWaspect, rnd
      INTEGER IR_, L, N, i, ix, iy, xadd, yadd, x(1200), y(1200),
     +       im1, IC1_
      LOGICAL Ltmp1_
      CALL GWvport(IR_, 0.0, 0.2, 0.5*GWaspect(1), 1.0)
      r = 0.5                               ! size of box
      CALL compute_aspect_ratio(REAL(L+r),xwin,ywin)
      CALL GWindow(IR_, -0.1*xwin, -0.1*ywin, xwin, ywin)
      CALL GWsavevp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, N + 1)
      CALL GWcolor(IR_, 19, 7)
      WRITE(Stmp1_,*) 'density = ', DBLE(N)/L/L
      CALL GWputtxt(IR_, 0.0, ywin*0.93, Stmp1_)
      CALL GWsetpen(IR_, 13, -1, -1, 4)
      CALL GWrect(IR_, 1-r, 1-r, L+r, L+r)
      DO iy = 1, L                          ! draw lattice sites 
         DO ix = 1, L
            site(ix,iy) = 0
            CALL GWsetpxl(IR_, REAL(ix), REAL(iy), 13)
         END DO
      END DO
      i = 0
      CALL GWanchor(IR_, 1)
      CALL GWsetpen(IR_, 16, -1, -1, 4)
      CALL GWcmbmrk(im1, 0, 6, 1.0, 16, -1, 4)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         xadd = int(L*rnd(0)) + 1
         yadd = int(L*rnd(0)) + 1
         IF(site(xadd,yadd) .EQ. 0) THEN
            i = i + 1                       ! number of particles added
            site(xadd,yadd) = -1            ! site occupied
            CALL GWsetogn(IR_, 0, i)
            CALL GWputcmb(IR_, im1, REAL(xadd), REAL(yadd), 0.0, 0.0, 0)
            x(i) = xadd
            y(i) = yadd
            x0(i) = x(i)                    ! x-coordinate at t = 0
            y0(i) = y(i)
         END IF
         Ltmp1_ = i .NE. N
      END DO
      END
C
      SUBROUTINE move(site,x,y,x0,y0,L,N,im1,IC1_)
      IMPLICIT NONE
      REAL*8 site(40,40), x0(1200), y0(1200)
      REAL rnd
      INTEGER IR_, particle, i, xtrial, ytrial, x(1200), y(1200), L, N,
     +       im1, IC1_
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, N + 1)
      DO particle = 1, N
         i = int(N*rnd(0)) + 1
         xtrial = x(i)
         ytrial = y(i)
         CALL direction(xtrial,ytrial,L)
         IF(site(xtrial,ytrial) .GE. 0) THEN    ! unoccupied site
            site(x(i),y(i)) = 0
            CALL GWsrect(IR_, x(i)-0.5,y(i)-0.5, x(i)+0.5,y(i)+0.5, -1)
            CALL GWerase(IR_, i, 0)
            CALL GWsetogn(IR_, 0, i)
            x(i) = xtrial
            y(i) = ytrial
            site(x(i),y(i)) = -1            ! new site occupied
            CALL GWputcmb(IR_, im1, REAL(x(i)), REAL(y(i)), 0.0, 0.0, 0)
         END IF
      END DO
      END
C
      SUBROUTINE direction(xtemp,ytemp,L)
*     choose random direction and use periodic boundary conditions
      IMPLICIT NONE
      REAL rnd
      INTEGER L, dir, xtemp, ytemp
      dir = int(4*rnd(0)) + 1
      IF(dir .EQ. 1) THEN
         xtemp = xtemp + 1
         IF(xtemp .GT. L) xtemp = 1
      ELSE IF(dir .EQ. 2) THEN
         xtemp = xtemp - 1
         IF(xtemp .LT. 1) xtemp = L
      ELSE IF(dir .EQ. 3) THEN
         ytemp = ytemp + 1
         IF(ytemp .GT. L) ytemp = 1
      ELSE IF(dir .EQ. 4) THEN
         ytemp = ytemp - 1
         IF(ytemp .LT. 1) ytemp = L
      ENDIF
      END
C
      SUBROUTINE mydata(x,y,x0,y0,L,L_2,N,t,im2,IC2_)
      IMPLICIT NONE
      REAL*8 L_2, R2, R2bar, dx, dy, x0(1200), y0(1200), t
      INTEGER IR_, i, x(1200), y(1200), L, N, im2, IC2_
      R2 = 0
      DO i = 1, N
         dx = x(i) - x0(i)
         dy = y(i) - y0(i)
         CALL separation(dx,dy,L,L_2)
         R2 = R2 + dx*dx + dy*dy
      END DO
      R2bar = R2/N
      CALL GWselvp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, N + 1)
      CALL GWputcmb(IR_, im2, REAL(t), REAL(R2bar), 0.0, 0.0, 0)
      END
C
      SUBROUTINE separation(dx,dy,L,L_2)
*     use periodic boundary conditions to determine separation
      IMPLICIT NONE
      REAL*8 L_2, dx, dy
      INTEGER L
      IF(abs(dx) .GT. L_2) dx = dx - sign(1.0D0,dx)*L
      IF(abs(dy) .GT. L_2) dy = dy - sign(1.0D0,dy)*L
      END