Let H _n = 1 + ½ + … + be the n-th partial sum of the harmonic series. A classical result of Wolstenholme states that, if p > 3 is prime, the numerator of H _p–l is divisible by p ². Here we consider, for a given prime p, the set J _p of n for which p divides the numerator of H _n. This set J _p had been previously determined for p = 2,3,5,7. One of our results is that J ₁₁ contains exactly 638 integers, the largestof which is a number of 31 decimal digits. We determine J _p for all p < 550 with three exceptions: 83, 127 and 397.The computation is based on a new p-adically convergent formula for the quantity H _pn – H _n/p. We describe a probabilistic model for the sets J _p, based on branching processes. The model predicts that |J _p| = O(p ²(log log p)^2+∊), and that there are infinitely many p with |J _p| ≥ p ²(log log p)². This strengthens an earlier conjecture of Eswarathasan and Levine that |J _p| is finite for all p. Another prediction of the model is that there will be infin itely many pairs (n,p) for which p ³ divides the numerator of H _n, but only finitely many for which p ⁴ divides H _n.It has been conjectured that there are infinitely many p for which |J_p | = 3. We give a probabilistic argument that suggests that such primes have a density 1/e in the set of all primes, and experimentally confirm this by a determination of all such p ≤ 10⁵.