You can also do this "topologically".  The idea is to separately calculate the (l-adic etale) Euler characteristics of the stack Y_0(p) in char. 0 and in char. p, then see how they must relate.  Here we go:

Over C, and hence over Q_p-bar as well, the Euler characteristic of Y_0(p) is (p+1)*(-1/12), since Y_0(p) is a (p+1)-fold cover of Y=M_{ell}.

On the other hand, over F_p-bar, up to "homeomorphism" Y_0(p) is two copies of Y glued at the supersingular points, so the Euler characteristic is 2*(-1/12) - S (where S is the "number" of supersingular elliptic curves).

However, since Y_0(p) has semi-stable reduction over Z_p, the generic fiber is gotten ("topologially") from the special fiber by fattening up the nodal points into closed cylinders then removing the open centers from these; thus our first (char 0) Euler characteristic is equal to our second (char p) Euler characteristic minus S.  Comparing and solving for S gives the formula pretty quickly.