Numerical reproducibility and stability of large scale scientific simulations, especially climate modeling, on distributed memory parallel computers are becoming critical issues. In particular, global summation of distributed arrays is most susceptible to rounding errors, and their propagation and accumulation cause uncertainty in final simulation results. We analyzed several accurate summation methods and found that two methods are particularly effective to improve (ensure) reproducibility and stability -- Kahan's self-compensated summation and Bailey's double-double precision summation. We provide an MPI operator MPI_SUMMDD to work with MPI collective operations to ensure a scalable implementation on a large number of processors. The final methods are particularly simple to adopt in practical codes.