To enhance the scalability and performance of the traditional finite-difference time-domain (FDTD) methods, a three-dimensional summation-by-parts simultaneous approximation term (SBP-SAT) FDTD method is developed to solve complex electromagnetic problems. It is theoretically stable and can be further used for multiple mesh blocks with different mesh sizes. This paper mainly focuses on the fundamental theoretical aspects upon its three-dimensional implementation, the SAT for various boundary conditions, and the numerical dispersion properties and the comparison with the FDTD method. The proposed SBP-SAT FDTD method inherits all the merits of the FDTD method, which is matrix-free, easy to implement, and has the same level of accuracy with a negligible overhead of runtime (0.13\%) and memory usage (1.2\%). Four numerical examples are carried out to validate the effectiveness of the proposed method.