I think it may be a little quicker to compute the mean differences as a matrix multiplication rather than using the mean subscript operator twice. It involves a lot more computation, but the advantage is not having to split the large matrix up.

y = j(N1, 1, 1/N1) // j(N-N1, 1, -1/(N-N1));
nulldist = sample(alldata, N//NRepl, "WOR") * y;

Does casting the problem as a matrix multiplication also have potential benefits from the multi-threading?