$z_1,\ldots,z_n,x_n$ together form a $n+1$-dimensional multivariate normal distribution. Proof: as $z$s are independent normal, linear combinations of them are normal. Furthermore, as $x_n$ is defined as a linear combination of the $z$s, any linear combination of the $z$s and $x_n$ is a linear combination of the $z$s and thus normal.
Mean of the multinormal distribution in question is $\mathbf{0}$. The covariance is such that the covariance of any $z_i,z_j$ is 0 if $i\neq j$ and 1 if $i=j$. Covariance of $z_i,x_n$ is $c_k$ and the variance of $x_n$ is $\sum_{k=1}^n c_n^2$. Now, to obtain the conditional distribution of $(z_1,\ldots,z_n)$ conditional on $x_n = 0$, we may directly apply the known equations for conditional distributions of multivariate normal distribution, with $\mathbf{x}_1 = (z_1,\ldots,z_n),~\mathbf{x}_2=x_n,\mathbf{a}=0$.
As the observed value 0 equals the mean of $x_n$, the conditional mean equals the unconditional mean, i.e.,
\begin{equation}
E((z_1,\ldots,z_n) \mid x_n = 0 ) = (0,\ldots,0).
\end{equation}
To obtain the covariance, we need to evaluate the Schur complement,
\begin{equation}
\bar{\Sigma} = \Sigma_{11} - \Sigma_{12}\Sigma^{-1}_{22}\Sigma_{21}.
\end{equation}
Plugging in the values of the distribution in question,
\begin{equation}
= I_n - \begin{pmatrix} c_1 & \ldots & c_n \end{pmatrix}^T\left(\sum_{k=1}^n c_k^2\right)^{-1}\begin{pmatrix}c_1 & \ldots & c_n \end{pmatrix}.
\end{equation}
Evaluating this element by element, we get
\begin{array}{lll}
Cov(z_i,z_i) & = & 1 - \frac{c_i^2}{\sum_{k=1}^n c_k^2},& i \in \{1,\ldots,n\} \\
Cov(z_i,z_j) & = & - \frac{c_ic_j}{\sum_{k=1}^n c_k^2}, & i\neq j;~ i,j \in \{1,\ldots,n\}.
\end{array}
Since the conditional distribution is multivariate normal, you can sample from it using standard methods, e.g., the Cholesky decomposition approach.
