Sampling correlated categorical variables

I am looking for a way to sample correlated categorical (non-binary) variables, and in particular I am interested in the category counts:

I have a set of $n$ correlated categorical random variables $X_i$ taking values $A$, $B$ or $C$, with marginal probabilities $p_{iA}$, $p_{iB}$, and $p_{iC}=1-p_{iA}-p_{iB}$. If I represent each $X_i$ with dummy variables $A_i$, $B_i$ and $C_i$, I can represent the correlation between $X_i$ and $X_j$ using 9 parameters $<cov(U_i,V_j)|U,V\in \{A,B,C\}>$. Note that if $i=j$ this covariance is fixed at $p_{iU}-p_{iU}^2$ or $-p_{iU}p_{iV}$ because they are dummy variables. I am interested in drawing samples from the distribution of counts of A, B and C (i.e. sample all $X_i$ and count up the number of $A$s, $B$s and $C$s, to give one sample of counts).

Some simplifications I would be okay with (if it helps):

1. $cov(U_i,V_j)=cov(V_i,U_j)$
2. $\forall_{i,U} p_{iU}=p_U$
3. $\forall_{i,U} p_{iU}=p_U$ and $cov(U_i,V_j)=f(i,j)\times cov(U_i,V_i)$

Some almost solutions:

1. If $\forall_{i,U} p_{iU}=p_U$ and $\forall_{i\ne j,U,V}cov(U_i,V_j)=0$, this would be a multinomial distribution.
2. If $n$ were small, I could specify the relative probabilities of all the outcomes and use the integral method, but $n$ is large enough that $3^n$ probabilities can't be reasonably stored.
3. If I had only two values, I could use methods described by Leisch et al. [Friedrich Leisch, Andreas Weingessel, and Kurt Hornik. On the generation of correlated artificial binary data. Technical Report 13, SFB Adaptive Information Systems and Modelling in Economics and Management Science, Wirtschaftsuniversit¨at Wien, Augasse 2-6, A-1090 Wien, Austria, 1998.].
4. I suspect that the right copula would allow me to extend Leisch et al. to the 3-value problem, but I would need a good way (analytic or fast numerical way) to specify the copula parameters in terms of $cov(U_i,V_j)$.

Any thoughts on how to sample either the categorical variables or their counts?