PROGRAM sandpile
*     simple one-dimensional sandpile model
      IMPLICIT NONE
      INTEGER N(0:51), grain, height(51)
      INTEGER IC1_, IC2_, IR_, L, move(100), topple, unstable
      LOGICAL TBkeyinput, Ltmp1_
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL initial(height,N,L,IC1_,IC2_)
      grain = 0
      DO WHILE(.NOT.TBkeyinput())
         height(1) = height(1) + 1          ! add grain to first site
         CALL GWselvp(IR_, IC1_)
         CALL GWsetogn(IR_, 0, L + height(1))
         CALL GWsrect(IR_, 1.1, height(1) - 0.4, 
     +                     1.9, height(1) + 0.4, 16)
         topple = 0                         ! number of grains that topple
         grain = grain + 1                  ! number of grains added
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            CALL check(height,L,move,unstable)
            IF(unstable .EQ. 1) THEN
               CALL slide(height,L,move,IC1_)
               topple = topple + 1
            END IF
            Ltmp1_ = unstable .NE. 0
         END DO
         N(topple) = N(topple) + 1
         CALL show_distribution(N,L,grain,IC2_)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(height,N,L,IC1_,IC2_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      INTEGER height(51), N(0:51)
      REAL GWaspect
      INTEGER IR_, i, topple, L, IC1_, IC2_
      CALL GWsetogn(IR_, 0, 1)
      CALL GWvport(IR_, 0.0, 0.0, 0.7*GWaspect(1), 1.0)
      CALL GWinput(IR_, 'lattice size L = ', Stmp1_)
      READ(Stmp1_,*) L                      ! suggest L = 20
      DO i = 1, L
         height(i) = 0
      END DO
      height(L+1) = 0                       ! boundary condition
      CALL GWindow(IR_, -1.0, -2.0, REAL(L+1), 40.0)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWline(IR_, 1.0, -0.05, REAL(L+1), -0.05)
      CALL GWsavevp(IR_, IC1_)
*     initialize array
      DO topple = 0, L
         N(topple) = 0
      END DO
      CALL GWvport(IR_, 0.72*GWaspect(1), 0.4, GWaspect(1), 1.0)
      CALL GWindow(IR_, -0.1, -0.05, REAL(L+1), 1.01)
      CALL GWline(IR_, 0.0, -0.01, L+0.2, -0.01)
      CALL GWline(IR_, -0.05, -0.01, -0.05, 1.0)
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetogn(IR_, 1, 0)
      END
C
      SUBROUTINE check(height,L,move,unstable)
      IMPLICIT NONE
      INTEGER height(51)
      INTEGER i, L, move(100), unstable
      unstable = 0
      DO i = 1, L
         IF(height(i) - height(i+1) .GT. 1) THEN
            move(i) = 1
            unstable = 1
         ELSE
            move(i) = 0
         END IF
      END DO
      END
C
      SUBROUTINE slide(height,L,move,IC1_)
      IMPLICIT NONE
      INTEGER height(51)
      INTEGER IR_, i, L, move(100), IC1_
      CALL GWselvp(IR_, IC1_)
      DO i = 1, L
         IF(move(i) .EQ. 1) THEN
            CALL GWerase(IR_, i*L+height(i), -1)
            height(i) = height(i) - 1       ! sand topples
            IF(i .LT. L) THEN
*              next site receives grain
               height(i+1) = height(i+1) + 1
               CALL GWsetogn(IR_, 0, (i + 1)*L+height(i + 1))
               CALL GWsrect(IR_, i + 0.1, height(i + 1) - 0.4, 
     +                           i + 0.9, height(i + 1) + 0.4, 16)
            END IF
         END IF
      END DO
      END
C
      SUBROUTINE show_distribution(N,L,grain,IC2_)
      IMPLICIT NONE
      INTEGER N(0:51), grain
      INTEGER IR_, topple, L, IC2_
*     erase previous graph
      CALL GWselvp(IR_, IC2_)
      DO topple = 0, L
         IF(N(topple) .GT. 0) THEN
            CALL GWerase(IR_, topple+10000, -1)
            CALL GWsetogn(IR_, 0, topple+10000)
            CALL GWsrect(IR_, REAL(topple), 0.0, 
     +                        topple+0.2, REAL(N(topple))/grain, 0)
         END IF
      END DO
      END