This is a partial answer but too long for a comment.
Let us first prove the property stated i.e.
$$
\mathbf{A} \mbox{ unitary }\Longrightarrow \mbox{ all roots of }p(z) \mbox{ are of modulus } 1
$$
Firsly $|p(0)|=|det(-\mathbf{A})|=1$,
so $p$ has not zero as a root.
Second, if $z\not=0$, one can write
$$
det(\mathbf{D}(z)-\mathbf{A})
=det(\mathbf{D}(z))det(I-\mathbf{D}(z)^{-1}\mathbf{A})
$$
so $p(z)=0$ is equivalent to $1\in sp(\mathbf{D}(z)^{-1}\mathbf{A})$ and to the existence of $\mathbf{v}\not=0$ such that
$$
\mathbf{D}(z)^{-1}\mathbf{A}\mathbf{v}=\mathbf{v}
$$
for such $\mathbf{v}$ one has
$
\mathbf{A}\mathbf{v}=\mathbf{D}(z)\mathbf{v}
$
and, writing $z=\rho e^{it}$ one gets
$
\mathbf{A}\mathbf{v}=\mathbf{D}(\rho)\mathbf{D}(e^{it})\mathbf{v}\ .
$
Now, $\mathbf{D}(e^{it})$ being unitary, one gets finally
$$
||\mathbf{v}||_2=||\mathbf{D}(e^{-it})\mathbf{A}\mathbf{v}||_2=||\mathbf{D}(\rho)\mathbf{v}||_2\qquad \mbox{(*).}
$$
As it is easy to check that for all $\mathbf{v}$ and $\rho>0$,
$$
||\mathbf{D}(\rho)\mathbf{v}||_2\geq \rho ||\mathbf{v}||_2 \mbox{ if } \rho>1\ ;\ ||\mathbf{D}(\rho)\mathbf{v}||_2\leq \rho ||\mathbf{v}||_2 \mbox{ if } \rho<1\ ,
$$
the result follows from (*).
For $\alpha=(m_1,\cdots ,m_N)\in (\mathbb{N}_{\geq 1})^N$, let
$$
\mathbf{D}(\alpha,z)=diag(z^{m_1}, \dots, z^{m_N})
$$
for $z$ fixed, one has
$$
\mathbf{D}(\alpha,z)\mathbf{D}(\beta,z)=\mathbf{D}(\alpha+\beta,z)
$$
these matrices form a semigroup.
For $\alpha$ fixed
$$
\mathbf{D}(\alpha,z_1)\mathbf{D}(\alpha,z_2)=\mathbf{D}(\alpha,z_1z_2)
$$
for $z\not=0$ these matrices form a group. The set of all these matrices is normalised by the monomial matrices, i.e. the semi-direct product of the (Weyl) group of permutation matrices and of diagonal matrices as, if $W=W(\sigma)$ is a permutation matrix and if $D$ is a diagonal (regular) matrix, one has
$$
WD\mathbf{D}(\alpha,z)D^{-1}W^{-1}=\mathbf{D}(\alpha_\sigma,z)
$$
As was remarked by Victor, matrices such as unitary (see above) and upper (or lower) triangular with unitary diagonal possess the property.
Their conjugates through the monomial group possess also the property.
How to test (algorithmically) that a matrix is conjugated of a unitary matrix through the monomial group ?
Firstly, as was remarked as the permutation matrices and as the diagonal ones with unitary spectrum are unitary, to be such is equivalent of being conjugated of a unitary matrix through the diagonal group of matrices with strictly positive eigenvalues i.e. for a matrix $B$ test whether it exists a unitary matrix and $R=diag(r_1,\cdots ,r_N)$ such that
$$
B=R^{-1}AR
$$
(one can even restrict to the special group of them, but we will not use this)
(analysis) suppose it were the case, then $B^*R^2B=R^2$
(synthesis)
- find all the diagonal matrices $D$ which fulfil $B^*DB=D$ (it is a linear system with $N$ variables, i.e. diagonal eigenvectors for the eigenvalue $1$ of the linear transformation $D\rightarrow B^*DB$).
- among them select, if possible, a $D$ with strictly positive spectrum and set $R=\sqrt{D}$ then $A=R^{-1}AR$ is unitary.
Remark It can happen that the transformation $D\rightarrow B^*DB$ admit diagonal eigenvectors for the eigenvalue $1$ none of which is strictly positive. As the procedure provides a necessary and sufficient condition, the corresponding matrix is not diagonally conjugate to a unitary matrix.
(Counter)-example Set
$$
B=
\begin{pmatrix}
\sqrt{2} & 1\cr
1 & \sqrt{2}
\end{pmatrix}
$$
then, solving
$$
B^*\begin{pmatrix}
x & 0\cr
0 & y
\end{pmatrix}
B
=
\begin{pmatrix}
x & 0\cr
0 & y
\end{pmatrix}
$$
yields $x=-y$ so there are eigenvectors as
$$
\begin{pmatrix}
1 & 0\cr
0 & -1
\end{pmatrix}
$$
but none of them is strictly positive.