Generating the merged gg->4l+0,1,2j

The class of loop-induced diagrams (with the vector bosons attached to a closed fermion loop) is a gauge-invariant class of contribution that you can indeed treat separately with the conventional LO merging procedures.

You are however right that start from two additional jets and beyond, you must do some extra work to select only this loop-induced contribution and disentangle it from the existing tree-level contributions and their loop corrections when the jets are quarks.
Two things must be done:

a) Make sure existing tree-level diagrams are ignored and the loops considered are squared against themselves.
-> This can be done by using [noborn=QCD] instead of [QCD] in the process definition.

b) Make sure that only the 'loop-induced class of diagrams are kept'. This cannot be done simply from the interactive interface by modifying the process definition. You will have to modify the MG5aMC code in a specific user-dedicated function to do this.
-> Modify the function 'user_filter' in 'madgraph/loop/loop_diagram_generation.py' as follows:
  b.1) Specify that you want to edit the loop_filter manually by changing:
              edit_filter_manually = False
           to
              edit_filter_manually = True
  b.2) Specify that you want to remove all diagrams which are not closed fermion loops with at least one vector boson attached to them, by adding the following code within the for-loop over the diagrams (i.e. somewhere below the "for diag in self['loop_diagrams']" next to the commented examples):
"""
            if any(abs(pdg) not in range(1,7) for pdg in diag.get_loop_lines_pdgs()) or \
               [abs(pdg) in [22,23,24,25] for pdg in diag.get_pdgs_attached_to_loop(structs)].count(True)<1:
                valid_diag = False
"""
This code will make sure that only closed fermion loops with at least one vector boson attached to it are considered.
(you can look at the diagrams selected with the command 'display diagrams')
Remember that with your suggested process definition you will include all non-resonant, singy-resonant and doubly-resonant diagrams, not only the doubly resonant ones which, depending on your analysis, you might really be interested in.

When present, you can (and should) include the tree-level contributions separately as well, because they are likely to contribute more than the loop-induced counterpart.
You can even include NLO corrections for those without fear of double-counting because even though the same loops could appear their, they will be squared against the tree-level diagrams, and not themselves. This is basically because the loop-induced contribution for these processes is basically a part of the NNLO correction, not NLO.

Finally, regarding the complexity of the computation, the 4 leptons come from the decay of two vector bosons, so they should really count as two (although the fact that you also includes non-resonant and singly resonant diagrams increases significantly the number of diagrams and integration channels).
So, if you have a cluster of reasonable size, you should be able to do 4l + 2j. To go beyond that, you will need considerably more resources, but there is no show-stopper from a technical standpoint in the code.