Frozen-density embedding (FDE) represents a versatile embedding scheme to describe the environmental effect on electron dynamics in molecular systems. The extension of the general theory of FDE to the real-time time-dependent Kohn-Sham method has previously been presented and implemented in plane waves and periodic boundary conditions [Pavanello, M.; J. Chem. Phys. 2015, 142, 154116]. In the current paper, we extend our recent formulation of the real-time time-dependent Kohn-Sham method based on localized basis set functions and developed within the Psi4NumPy framework to the FDE scheme. The latter has been implemented in its "uncoupled" flavor (in which the time evolution is only carried out for the active subsystem, while the environment subsystems remain at their ground state), using and adapting the FDE implementation already available in the PyEmbed module of the scripting framework PyADF. The implementation was facilitated by the fact that both Psi4NumPy and PyADF, being native Python API, provided an ideal framework of development using the Python advantages in terms of code readability and reusability. We employed this new implementation to investigate the stability of the time-propagation procedure, which is based on an efficient predictor/corrector second-order midpoint Magnus propagator employing an exact diagonalization, in combination with the FDE scheme. We demonstrate that the inclusion of the FDE potential does not introduce any numerical instability in time propagation of the density matrix of the active subsystem, and in the limit of the weak external field, the numerical results for low-lying transition energies are consistent with those obtained using the reference FDE calculations based on the linear-response TDDFT. The method is found to give stable numerical results also in the presence of a strong external field inducing nonlinear effects. Preliminary results are reported for high harmonic generation (HHG) of a water molecule embedded in a small water cluster. The effect of the embedding potential is evident in the HHG spectrum reducing the number of the well-resolved high harmonics at high energy with respect to the free water. This is consistent with a shift toward lower ionization energy passing from an isolated water molecule to a small water cluster. The computational burden for the propagation step increases approximately linearly with the size of the surrounding frozen environment. Furthermore, we have also shown that the updating frequency of the embedding potential may be significantly reduced, much less than one per time step, without jeopardizing the accuracy of the transition energies.