We perform a numerical study of black hole formation from the spherically symmetric collapse of a massless scalar field. The calculations are done in Painleve-Gullstrand (PG) coordinates that extend across apparent horizons and allow the numerical evolution to proceed until the onset of singularity formation. We generate spacetime maps of the collapse and illustrate the evolution of apparent horizons and trapping surfaces for various initial data. We also study the critical behaviour and find the expected Choptuik scaling with universal values for the critical exponent and echoing period consistent with the accepted values of $gamma=0.374$ and $Delta = 3.44$, respectively. The subcritical curvature scaling exhibits the expected oscillatory behaviour but the form of the periodic oscillations in the supercritical mass scaling relation, while universal with respect to initial PG data, is non-standard: it is non-sinusoidal with large amplitude cusps.