Here I use a latent variable approach. This readily extends to the continuous/categorical case.
The idea is to treat a continuous variable (the latent variable) as laying behind the ordered categories that are actually observed (by splitting up the continuous variable at breakpoints).
So for the two variables that are independent, we define breakpoints that give the desired proportions in each category. Then the third continuous variable, correlated with the other two, is also split up in similar fashion. It's common to use standardized normal variables for the latent variables, but other distributions could be used.
The example below is in R but I have annotated it to help conversion to other platforms.
set.seed(10345) # just to make sure if you run this we have the same results
xu=rnorm(50) # draw 50 observations from continuous latent variables
yu=rnorm(50) #
zu= 0.8*xu+0.6*yu # the latent variables have correlations 0 between x and y,
# 0.8 between x and z, and 0.6 between y and z
cor(cbind(xu,yu,zu)) # sample correlations will be similar to those population values
px=c(.3,.2,.5) # our selected population proportions in the marginal categories
py=c(.1,.2,.4,.3)
pz=c(.1,.2,.4,.2,.1)
xc=cut(xu,qnorm(cumsum(c(0,px))),labels=c("AI","AII","AIII")) # convert to ord. categ.
yc=cut(yu,qnorm(cumsum(c(0,py))),labels=LETTERS[1:4])
zc=cut(zu,qnorm(cumsum(c(0,pz))),labels=letters[1:5])
Now let's see the relationships between variables:
table(xc,yc) #examine the resulting data. xc,yc populations are independent
yc
xc A B C D
AI 1 7 9 2
AII 0 4 11 7
AIII 2 5 18 14
> table(xc,zc) #xc,zc dependent
zc
xc a b c d e
AI 4 11 4 0 0
AII 0 2 19 1 0
AIII 0 1 18 12 8
> table(yc,zc) #yc,zc dependent
zc
yc a b c d e
A 1 1 1 0 0
B 2 7 5 1 1
C 1 5 27 5 0
D 0 1 8 7 7
How correlations between the latent variables work.
I chose $X_u$ and $Y_u$ ($u$ for "underlying"; I'd have put $l$ for "latent", but it tends to look like a "1") to be two independent standard normal variates. You can make them correlated with a third variate, $Z_u$, by making $Z_u$ a linear combination of $X_u$, $Y_u$, and an independent noise variate $\epsilon$, which we'll also take to be standard normal here.
If we write $Z^*=aX_u+bY_u+c\epsilon$ then $Z^*$ is normal, but not standard normal.
$\text{Cov}(Z^*,X_u)=\text{Cov}(aX_u+bY_u+c\epsilon,X)=a\,\sigma^2_X=a$
Similarly $\text{Cov}(Z^*,Y_u)=b$ and $\text{Cov}(Z^*,\epsilon)=c$.
$\text{Var}(Z^*)=a^2+b^2+c^2$
So $\text{Cor}(Z^*,X_u)=\frac{a}{\sqrt{a^2+b^2+c^2}}$ and So $\text{Cor}(Z^*,Y_u)=\frac{b}{\sqrt{a^2+b^2+c^2}}$.
But I want $Z_u$ to have variance $1$, so if we define $Z_u=\frac{Z^*}{\sqrt{a^2+b^2+c^2}}$ then
$\text{Var}(Z_u)=\frac{a^2+b^2+c^2}{a^2+b^2+c^2}=1$
In the example, I chose $a=0.8,b=0.6,c=0$, which has $a^2+b^2+c^2=1$ and in that case $Z_u=Z^*$, and we have $\text{Cor}(Z_u,X_u)=a=0.8$ and $\text{Cor}(Z_u,Y_u)=b$.
If you choose to have $\text{Cor}(Z_u,X_u)=\rho\,,$ then $-\sqrt{1-\rho^2}\leq\text{Cor}(Z_u,Y_u)\leq \sqrt{1-\rho^2}$ (with the limits being achieved when $c=0$).
Note that these are population correlations, not sample correlations.
In the example you mention in comments, $a=b=\frac{1}{2}$, and $c=0$ which gives $\text{Cor}(Z^*,X_u)=\frac{a}{\sqrt{a^2+b^2+c^2}}=\frac{1/2}{\sqrt{(1/2)^2+(1/2)^2}}=\sqrt{\frac{1}{2}}\approx 0.7071$
-- but now to make $Z_u$ standard normal we need to divide through by
$\sqrt{a^2+b^2+c^2}=\sqrt{\frac{1}{2}}$, i.e.
$Z_u=Z^*/\sqrt{\frac{1}{2}}=\sqrt{2}Z^*$.