random_number_init Subroutine

  subroutine random_number_init(shift)
    integer,optional                 :: shift
    integer                          :: i, n, clock
    integer,dimension(:),allocatable :: seed
    call RANDOM_SEED(size = n)
    allocate(seed(n))
    call SYSTEM_CLOCK(COUNT=clock)
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
    if(present(shift))seed=seed+shift
    call RANDOM_SEED(PUT = seed)
    deallocate(seed)
  end subroutine random_number_init