Algorithms for creating a directed network with a given 3-node motifs distribution
As I understand, you are interested in an approximate solution. I have two heuristic ideas, maybe they can serve as a starting point.

However, if you want more (e.g. randomness in the Haar-measure sense or fixing other properties (graph diameter, spectrum, number of connected components, ...)) they are not sufficient. 

Rewiring

Especially when you want to have distributions of network motifs close to a 'random' one.


Generate a graph with desired distribution of in-degrees and out-degrees.
For each pair of nodes $(x,y)$ calculate the number of each motif $m_i(x,y)$
Now, in each step:

Randomly draw two edges.
Perform rewiring (i.e. $(a\rightarrow b, c \rightarrow d) \mapsto (a\rightarrow d, c \rightarrow b) $) if it is possible and if it results in distribution of motifs closer to the desired one (e.g. in the sense of the Kullback-Leibler divergence or perhaps the Euclidian distance). Modify $4\times13$ affected $m_i(x,y)$.



As variance is of order $\sqrt{n}$, each rewiring takes a fixed time, you need roughly $n^3/\sqrt{n}$ steps, so I guess that it should scale as $O(n^{2.5})$.

Addition

Alternatively, you can try to add edges to an empty graph in such way it produces the desired distribution. Also, one need to hold $m_i(x,y)$.

For each step one randomly draws a pair of nodes and put an edge if it results in a desired distribution at this step (perhaps with some randomness).

Then the time cost is $O(n^2)$ (or rather $O(\#\text{edges})$ ), so better by $\sqrt{n}$ than the rewiring.

The tricky part is to calculate the desired distribution at each step (or rather - at each edge count). For example, even if one wants to have only $(\bullet\leftrightarrows \bullet\quad \bullet)$ and $(\bullet\quad \bullet\quad \bullet)$ 3-motifs, then when starting one needs to create mostly $(\bullet\rightarrow \bullet\quad \bullet)$ 3-motifs.

I would guess that the desired distribution at a given edge count $E$ is the one resulting from the final motif distribution with removed edges with the probability $1-\frac{E}{E_{final}}$.