One possibility would be to use blavaan for this, which uses Stan under the hood. You might be able to use mplus2lavaan() from the lavaan package to translate your model to lavaan, then pass the lavaan model to blavaan.

Or, if you wanted to do it directly in blavaan: assuming that A1-A3 are either binary or ordinal (as opposed to >2 unordered categories, which are not supported), you could do something like

m1 <- ' A =~ A1 + A2 + A3
B =~ B1 + B2 + B3
B1 ~~ A2
B2 ~~ A3 '
fit <- bcfa(m1, data = mydata, ordered = paste0("A", 1:3))

Possibly with an extra line to specify that the correlation between A and B is 0 (if desired, or it is possibly required to identify the model).

The underlying Stan code is general and not tailored to your model, so not best for learning Stan. It would be worth looking at brms as well, see this thread.

Thank you for the suggestion. I have attempted to use blavaan for this problem before, however, my parameter estimates were biased. I wanted to have a Stan code of my own to have more freedom with attempts to reparameterize the model.

I would be interested to see the evidence of parameter bias in blavaan.

For coding this model yourself, you might look at the multivariate probit example at the bottom of this page. It does the Albert-Chib thing, where you end up with z parameters that are the continuous responses underlying ordinal variables. Once you have those z parameters, you could model them jointly with the continuous data in the usual way.

That example involves binary data, and it is not trivial to extend to variables with >2 ordered categories. I could add more detail there if you need it.