  We report a study of the decay $D^0 \to K^0_S K^0_S$ using 921~fb$^{-1}$ of
data collected at or near the $\Upsilon(4S)$ and $\Upsilon(5S)$ resonances with
the Belle detector at the KEKB asymmetric energy $e^+e^-$ collider. The
measured time-integrated $CP$ asymmetry is $ A_{CP}(D^0 \to K^0_S K^0_S) =
(-0.02 \pm 1.53 \pm 0.02 \pm 0.17) \%$, and the branching fraction is
$\mathcal{B} (D^{0}\rightarrow K_{S}^{0}K_{S}^{0})$ = (1.321 $\pm$ 0.023 $\pm$
0.036 $\pm$ 0.044) $\times$ 10$^{-4}$, where the first uncertainty is
statistical, the second is systematic, and the third is due to the
normalization mode ($D^0 \to K_S^0 \pi^0$). These results are significantly
more precise than previous measurements available for this mode. The $A_{CP}$
measurement is consistent with the standard model expectation.
