random_order Subroutine

  SUBROUTINE random_order(order, n)
    !     generate a random ordering of the integers 1 ... n.
    integer, intent(in)  :: n
    integer, intent(out) :: order(n)
    !     local variables
    integer :: i, j, k
    real(8) :: wk
    do i = 1, n
       order(i) = i
    end do
    !     starting at the end, swap the current last indicator with one
    !     randomly chosen from those preceeding it.
    do i = n, 2, -1
       call random_number(wk)
       j = 1 + i * wk
       if (j < i) then
          k = order(i)
          order(i) = order(j)
          order(j) = k
       end if
    end do
    return
  end subroutine random_order