We propose a modification to the nonadiabatic surface hopping calculation method formulated in a paper by Yu et al. [Phys. Chem. Chem. Phys. 16, 25883 (2014)], which is a multidimensional extension of the Zhu-Nakamura theory with a practical diabatic gradient estimation algorithm. In our modification, their diabatic gradient estimation algorithm, which is based on a simple interpolation of the adiabatic potential energy surfaces, is replaced by an algorithm using the numerical derivatives of the adiabatic gradients. We then apply the algorithm to several models of nonadiabatic dynamics, both analytic and ab initio models, to numerically demonstrate that our method indeed widens the applicability and robustness of their method. We also discuss the validity and limitations of our new nonadiabatic surface hopping method while considering in mind potential applications to excited-state dynamics of biomolecules or unconventional nonadiabatic dynamics such as radiation decay processes in ultraintense X-ray fields.