Monte Carlo Pi with R + Fortran + MPI
Here are a few pointers. I should first say that I am more familiar with R;  I know next to nothing about Fortran but I could read your code and guess about a few things that will likely improve your code: both in simplicity and speed.

1) Why the need for two inputs DARTS and ROUNDS? You know that in the end you will have a total of DARTS * ROUNDS simulations: could there just be a single input num_simulations that will be split evenly among your numprocs? This should get you rid of the loop inside snowpi so it could only call dboard once.

2) You should make your slave (snowpi) as simple as possible. IMO, it should not know about things like proc_num or numprocs if you just tell it how many simulations it needs to do. Let R do the math upstream and pass it down to the slave. I don't know much about random generators but maybe the same can apply with the random seeds: if you could let R compute a seed for each slave and pass it down as an (optional) argument to the slave, that would be best. (However, if knowing the proc_num and numprocs is crucial for seed selection and can only happen on the Fortran side, then disregard my comment...)

3) This one I feel is important: do not compute averages too early. They are costly and can make you lose precision. This line in your code (avepi = ((avepi*i) + pi_est)/(i + 1)) is particularly problematic: you are doing so many unnecessary multiplications and divisions... Instead, you should just be counting the number of darts, i.e. incrementing a counter. Your slave should just be returning an integer. Let R add all these integers together and do a single division by num_simulations at the end.

To summarize 1, 2, and 3, it would be great (if possible) if you could just tell your slave: here is your seed number n; now throw x darts and tell me how many ended up in the circle.

4) Now looking at your R code. Yes, the use of &lt;&lt;- in your code is not very elegant. That's because you are trying to use clusterEvalQ which may not be the best option from the snow package. Instead, clusterApply or clusterMap are better suited for functional programming. Your code could look like this, I let you figure out the rest:

FSnowPi &lt;- function(num_simulations, seed_num = 123456L) {
  retvals &lt;- .Fortran("snowpi",
                      num_inside      = 1L,  # an integer
                      num_simulations = as.integer(num_simulations),
                      seed_num = as.integer(seed_num))
  return(retvals$num_inside)
}

SnowPi  &lt;- function(num_simulations, num_procs) {
  # input checks
  stopifnot(is.integer(num_simulations),
            is.integer(num_procs),
            num_simulations &gt; 0L,
            num_procs &gt; 0L)
  # initialize cluster
  cl &lt;- snow::makeCluster(num_procs, type = "MPI")
  # balance the number of simulations per node
  node_num_simu &lt;- balance_evenly(num_simulations, num_procs)
  # compute random seeds
  node_seeds &lt;- sample.int(1e6, num_procs)   # please review...
  # run in parallel
  num_inside &lt;- snow::clusterMap(cl, MyPi::FSnowPi, node_num_simu, node_seeds)
  # close the cluster
  snow::stopCluster(cl)

  pi_estimate &lt;- 4 * Reduce(`+`, num_inside) / num_simulations
  return(pi_estimate)
}

balance_evenly &lt;- function(total, among) {
   min_load &lt;- floor(total / among)
   max_load &lt;- min_load + 1L
   rest &lt;- total - min_load * among
   load &lt;- c(rep(max_load, rest),
             rep(min_load, among - rest))
   return(load)
}


I hope some of it is helpful. Please let me know. And thanks for sharing your Fortran code; I know at least one other R package that uses it:  fields whose rdist() function for computing euclidean distances is extremely fast. I want to give a try one day.