We report in this paper our numerical analysis of energy level spacing statistics for the one-dimensional spin-$1/2$ XXZ model in random on-site longitudinal magnetic fields $B_i$ ($-hleq B_ileq h$)). We concentrate on the strong disorder limit $J_{perp}<<J_z,h)$ where $J_z$ and $J_{perp}$ are the (nearest neighbor) spin interaction strength in $z$- and planar ($xy$)- directions, respectively. The system is expected to be in a many-body localized (MBL) state in this parameter regime. By analyzing the energy-level spacing statistics as a function of strength of random magnetic field $h$, energy of the many-body state $E$, the number of spin-$uparrow$ particles in the system $M=sum_i(s_i^z+{1over2})$ and the spin interaction strengths $J_z$ and $J_{perp}$, we show that there exists a small parameter region $J_zsim h$ where ergodic behaviour emerges at the middle of the many-body energy spectrum when $Msim{Nover2}$ ($N=$ length of spin chain). The emerging ergodic phase shows qualitatively different behaviour compared with the usual ergodic phase that exists in the weak-disorder limit.