PROGRAM Shielding_a_Nuclear_Reactor
!-----------------------------------------------------------------------
! This program uses a random number generator to simulate neutrons 
! entering a shield and to determine what percentage reaches the
! outside.  The neutrons are assumed to move forward, backward, left
! and right with equal likelihood and to die within the shield if a
! certain number of changes of direction have occurred.  Identifiers 
! used are:
!     Thickness            : thickness of shield
!     DirectionChangeLimit : limit on # of direction changes before
!                            energy dissipated
!     RandomReal           : a random real number in the range 0 to 1
!     NewDirection         : a random integer 1, 2, 3, or 4 
!                            representing direction
!     OldDirection         : previous direction of neutron
!     NumDirectionChanges  : number of changes of direction
!     Forward              : net units forward traveled
!     NumNeutrons          : number of neutrons simulated
!     NumEscaped           : number of neutrons reaching outside of
!                            shield
!     I                    : loop control variable
!
! Input:   Thickness, DirectionChangeLimit, and NumNeutrons
! Output:  Percentage of neutrons that reach the outside
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: RandomReal
  INTEGER :: Thickness, DirectionChangeLimit, NewDirection, &
             OldDirection, NumDirectionChanges, Forward, &
             NumNeutrons, NumEscaped, I

  PRINT *, "Enter thickness of shield, limit on # of direction"
  PRINT *, "changes, and the number of neutrons to simulate:"
  READ *, Thickness, DirectionChangeLimit, NumNeutrons

  NumEscaped = 0

  ! Begin the simulation

  CALL RANDOM_SEED
  DO I = 1, NumNeutrons
     Forward = 0
     OldDirection = 0
     NumDirectionChanges = 0

    ! Repeat the following until neutron reaches outside of
    ! shield, returns inside reactor, or dies within shield

    DO
       CALL RANDOM_NUMBER(RandomReal)
       NewDirection = 1 + INT(4 * RandomReal)
       IF (NewDirection /= OldDirection) THEN
          NumDirectionChanges = NumDirectionChanges + 1
          OldDirection = NewDirection
       END IF
       IF (NewDirection == 1) THEN
          Forward = Forward + 1
       ELSE IF (NewDirection == 2) THEN
          Forward = Forward - 1
       END IF
       IF ((Forward >= Thickness) .OR. (Forward <= 0) .OR.  &
           (NumDirectionChanges >= DirectionChangeLimit)) EXIT
    END DO

    IF (Forward == Thickness)  NumEscaped = NumEscaped + 1

  END DO

  PRINT '(1X, F5.2, "% of the neutrons escaped")', &
        100 * REAL(NumEscaped) / REAL(NumNeutrons)

END PROGRAM Shielding_a_Nuclear_Reactor