Quantile of a conditional cumulative bivariate normal distribution and R Implementation

I would like to compute a quantile Z of a conditional cumulative bivariate normal distribution. I have two questions: 1) Do I formulate the problem correctly? 2) How would I compute this quantile Z with R?

1) Problem Formulation: The joint distribution of X,Y is given by: \begin{align} \left(X,Y \right) \sim N\left(0,\Sigma \right) \end{align} The conditioning set for Y is that Y is below its q% quantile $\{Y < \phi^{-1}(q)\sigma_{22}\}$. The conditonal distribution should equal p. Im interested in the p$\%$ quantile Z. Using the general law of conditional probabilities I get:

\begin{align*} p= \Pr \left(X < Z |Y < \phi^{-1}(q)\sigma_{22}\right) = \frac{\int_{-\infty}^{Z} \int_{-\infty}^{\phi^{-1}\sigma_{22}} f_{X,Y}(u,v)\, du \,dx}{q}. \end{align*}

Here $f_{X,Y}$ is the joint normal density which is known. This is how far I get. Can I summarize this expression any further?

2) Given $\Sigma,q,p$, how can I compute Z? I know that I could just simulate the distribution and form the empirical conditional distribution function. Are there any other, maybe more elegant ways that use the nice properties of the multivariate normal distribution?