  A multiscalar and nonrenormalizable $B-L$ extension of the standard model
(SM) with $S_4$ symmetry which successfully explains the recent observed
neutrino oscillation data is proposed. The tiny neutrino masses and their
hierarchies are generated via the type-I seesaw mechanism. The model reproduces
the recent experiments of neutrino mixing angles and Dirac CP violating phase
in which the atmospheric angle $(\theta_{23})$ and the reactor angle
$(\theta_{13})$ get the best-fit values while the solar angle $(\theta_{12})$
and Dirac CP violating phase ($\delta $) belong to $3\, \si $ range of the
best-fit value for normal hierarchy (NH). For inverted hierarchy (IH),
$\theta_{13}$ gets the best-fit value and $\theta_{23}$ together with $\de $
belongs to $1\, \si $ range while $\theta_{12}$ belongs to $3\, \si $ range of
the best-fit value. The effective neutrino masses are predicted to be $\langle
m_{ee}\rangle=6.81 \,\, \mbox{meV}$ for NH and $\langle m_{ee}\rangle=48.48\,\,
\mbox{meV}$ for IH being in good agreement with the most recent experimental
data.
