C  TM-GAS with walls, after Tiffoli and Margolus.
C  The universe is represented as an integer array
C  of 0's and 1's ( 1 = `present') using the
C  Margolus neighborhoods.  See notes for the
C  details.  This version has terminal output for
C  initial testing.
C
C  File:  tmgas.f
C  Headers:  tmgas.hp  tmgas.hc  tmgas.hf
C
C  Parameters are 
C
C       R, C - dimensions of the universe
C
C  Common information in /TMGASFE/, /TMGASPE/
C
C    I, J     - Array coordinates
C    EM,  OM  - Position for even, odd phases
C    MHE, MVE - Horiz., vert. flags for even phase
C    MHO, MVO - Similarly for odd phase
C    WLE, WLO - True in cells blocked by walls
C    PH       - logical, .FALSE. on even phase
C
C  Phase conventions ( X = O or E ):
C
C              MHX
C	      F   T
C	    +---+---+
C	  F | 0 | 1 |  Entry in box
C     MVX   +---+---+  is XM
C	  T | 2 | 3 |
C	    +---+---+


C  Run Nc cycles of TM-GAS on the universe U,
C  leaving final state in U.

      subroutine tmgas(U, Nc)
      INCLUDE 'tmgas.hc'
      integer U(R, C)
CMF$  LAYOUT U(,)
CMPF  ONDPU U
      integer Nc

      integer, array(R, C) :: UH, UV, UW
      logical, array(R, C) :: MHX, MVX, WLX, CLL
      integer kc

C   UH = state of opposite horiz cell
C   UV = 2*(opposite vert) + (opp diag)
C   MHX, MVX are flags for current phase
C   CLL = Collision flag  

D     print *, U(1, 1:5)
D     print *, U(2, 1:5)
D     print *, U(3, 1:5)
D     print *, U(4, 1:5)
D     print *, U(5, 1:5)
      do, k = 1, Nc
	if ( PH ) then
	  MHX = MHO
          MVX = MVO
	  WLX = WLO
	else
	  MHX = MHE
	  MVX = MVE
	  WLX = WLE
	end if
C   Get block state to every cell
	where ( MHX )
	  UH = CSHIFT(U, DIM=2, SHIFT=-1)
	elsewhere
	  UH = CSHIFT(U, DIM=2, SHIFT= 1)
	end where
	UW = ishft(U, 1) + UH      ! 2*x + opp
	where ( MVX )
	  UV = CSHIFT(UW, DIM=1, SHIFT=-1)
	elsewhere
	  UV = CSHIFT(UW, DIM=1, SHIFT= 1)
	end where
C   Identify blocks with collisions or walls
	CLL =     ( UW == 1 .and. UV == 2 )
     *       .or. ( UW == 2 .and. UV == 1 )
     *       .or. WLX
C   Calculate new states
	if ( PH ) then          ! Odd phase
	  where ( MVX .eqv. MHX )
	    UW = UH
	  elsewhere
	    UW = ishft(UV, -1)
	  end where
	else                    ! Even phase
	  where ( MVX .eqv. MHX )
    	    UW = ishft(UV, -1)
	  elsewhere
	    UW = UH
	  end where
	end if
	where ( .not. CLL ) U = UW
	PH = .not. PH           ! Phase complete
      end do
      end subroutine tmgas

C  Print upper left corner (rp, cp) on terminal

      subroutine prgas(U, rp, cp, cy)
      INCLUDE 'tmgas.hp'
      integer  U(R, C), rp, cp, cy
CMF$  LAYOUT U(,)
CMPF  ONDPU U

      character*100 buff
      integer j, k
      print 21, cy
 21   FORMAT (/"Time =", I6/)
      do, j = 1, rp
	do, k = 1, cp
	  if ( U(j, k) == 0 ) then
	    buff(k:k+1) = " "
	  else
	    buff(k:k+1) = "H"
	  end if
	end do ! k
	print '(A)', buff(1:cp)
      end do ! j
      print *              ! Blank line
      end subroutine prgas

C  Initialize system with no walls

      subroutine tminit()
      INCLUDE 'tmgas.hc'

      I = SPREAD([1:R], 2, C)
      J = SPREAD([1:C], 1, R)
      where ( btest(I, 0) )     ! I is odd
	EM = 1
      elsewhere
	EM = 0
      endwhere
      where ( btest(J, 0) ) EM = EM + 2
      OM = CSHIFT(CSHIFT(EM, DIM=1, SHIFT=1),
     *            DIM=2, SHIFT=2)
      MHE = btest(EM, 0)
      MVE = btest(EM, 1)
      MHO = btest(OM, 0)
      MHO = btest(OM, 1)
      WLE = .FALSE.
      WLO = .FALSE.
      PH  = .FALSE.
      end subroutine tminit

C  Construct matrix of random integers  0 <= IR <= 100
C  This is needed only for mpfortran

      subroutine rmat(IR, r, c)
      integer r, c, IR(r, c)
CMPF  ONDPU IR

CMPF  mpl MPF_RANDOM, MPF_SRANDOM
      real, parameter :: Two31 = 2147483648.0
      integer RR(r), k

      call mpf_srandom([1:R] + 12131992)
      do, k = 1, r
        call mpf_random(RR)
        IR(:, k) = 100.0*(RR/Two31) + 0.5
      end do
      end subroutine rmat

C  Interactive program to manage the universe

      program Zeus  ! (also Jupiter)
      INCLUDE 'tmgas.hp'
      integer pr, cy, ty, Ux(R, C)
      integer eltm     ! *** mpp ***

      print 20, R, C
 20   FORMAT (" Universe is a ", I4, " by ", I4, " torus.")
      call tminit()
C  Later we will construct wall here.
      print '(" Probability of particle (0-100) > " $)'
      read *, pr
C  Random initial state of the universe
      print *, "Generating random universe ...."
      call rmat(Ux, R, C)
C      call CMF_random(Ux, 101)     ! *** CM5 ***
      where ( Ux <= pr )
	Ux = 1
      elsewhere
	Ux = 0
      end where
      ty = 0		! Total cycles so far
      call prgas(Ux, 16, 64, ty)
      print *, "Specify 0 or fewer cycles to quit."
      do while (.true.)
	print '(" Number of cycles > " $)'
	read *, cy
	if ( cy <= 0 ) go to 99    ! Finished
        call mpTimerStart()     ! *** mpp ***
C	call CM_timer_clear(1)  ! *** CM5 ***
C	call CM_timer_start(1)  ! *** CM5 ***
	call tmgas(Ux, cy)
        eltm = mpTimerElapsed() ! *** mpp ***
        print 21, eltm          ! *** mpp ***
 21     FORMAT ("Elapsed time = ", I5, " msec.")
C	call CM_timer_stop(1)   ! *** CM5 ***
C	call CM_timer_print(1)  ! *** CM5 ***
	ty = ty + cy
	call prgas(Ux, 16, 64, ty)
      end do
 99   continue  ! For exit
      end program Zeus

