The three lowest full three-dimensional adiabatic and three diabatic global potential energy surfaces are reported for the title system. The accurate ab initio method (MCSCF/MRCI) with larger basis sets (aug-cc-pVQZ) is used to reduce the adiabatic potential energies, and the global adiabatic potential energy surfaces are deduced by a three-dimensional B-spline fitting method. The conical intersections and the mixing angles between the lowest three adiabatic potential energy surfaces are precisely studied. The most possible nonadiabatic reaction pathways are predicted, i.e., N(D) + H(X∑) → NH(2A') → CI (1A'-2A') → NH(1A') → CI (1A″-1A') → NH(1A″) → NH(X∑) + H(S). The products of the first excited state (NH(aΔ) + H(S)) and the second excited state (NH(b∑) + H(S)) can be generated in these nonadiabatic reaction pathways too.