We have to consider cases. Denote $B(y,z;\rho)$ the joint bivariate normal cumulative distribution function of two correlated variables with correlation coefficient $\rho$.
The joint support is $[0,\infty] \times (-\infty, \infty)$
For $\{Y=0, Z \in (-\infty, \infty)\}$ we have
$$P(Y=0, Z\le z) = P(X\le 0 , Z\le z) = \int_{-\infty}^{z}\int_{-\infty}^0f_{XZ}(x,z)dxdz$$
$$= B(0,z;\rho_{XZ})$$
For $\{Y>0, Z \in (-\infty, \infty)\}$ we have
$$P(Y>0, Z\le z) = P(0<X\le x , Z\le z) = \int_{-\infty}^{z}\int_{0}^xf_{XZ}(x,z)dxdz$$
$$= B(x,z;\rho_{XZ}) = B(y^+,z;\rho_{XZ})$$
the last equality because for this range $Y = X$.
Bringin together using indicator functions
$$F_{YZ}(y,z) = B(0,z;\rho_{XZ})\cdot I_{\{Y=0\}} + B(y^+,z;\rho_{XZ})\cdot (1-I_{\{Y=0\}})$$
Differentiation of the two branches of the cdf will give you the corresponding densities.