How do I make a single legend for many subplots?

There is also a nice function get_legend_handles_labels() you can call on the last axis (if you iterate over them) that would collect everything you need from label= arguments:

handles, labels = ax.get_legend_handles_labels()
fig.legend(handles, labels, loc='upper center')

If the pyplot interface is being used instead of the Axes interface, use:

handles, labels = plt.gca().get_legend_handles_labels()

To remove legends from subplots, see Remove the legend on a matplotlib figure.

To merge txinx legends, see Secondary axis with twinx(): how to add to legend.

figlegend may be what you're looking for: matplotlib.pyplot.figlegend

An example is at Figure legend demo.

Another example:

plt.figlegend(lines, labels, loc = 'lower center', ncol=5, labelspacing=0.)

Or:

fig.legend(lines, labels, loc = (0.5, 0), ncol=5)

For the automatic positioning of a single legend in a figure with many axes, like those obtained with subplots(), the following solution works really well:

plt.legend(lines, labels, loc = 'lower center', bbox_to_anchor = (0, -0.1, 1, 1),
           bbox_transform = plt.gcf().transFigure)

With bbox_to_anchor and bbox_transform=plt.gcf().transFigure, you are defining a new bounding box of the size of your figureto be a reference for loc. Using (0, -0.1, 1, 1) moves this bounding box slightly downwards to prevent the legend to be placed over other artists.

OBS: Use this solution after you use fig.set_size_inches() and before you use fig.tight_layout()

You just have to ask for the legend once, outside of your loop.

For example, in this case I have 4 subplots, with the same lines, and a single legend.

from matplotlib.pyplot import *

ficheiros = ['120318.nc', '120319.nc', '120320.nc', '120321.nc']

fig = figure()
fig.suptitle('concentration profile analysis')

for a in range(len(ficheiros)):
    # dados is here defined
    level = dados.variables['level'][:]

    ax = fig.add_subplot(2,2,a+1)
    xticks(range(8), ['0h','3h','6h','9h','12h','15h','18h','21h']) 
    ax.set_xlabel('time (hours)')
    ax.set_ylabel('CONC ($\mu g. m^{-3}$)')

    for index in range(len(level)):
        conc = dados.variables['CONC'][4:12,index] * 1e9
        ax.plot(conc,label=str(level[index])+'m')

    dados.close()

ax.legend(bbox_to_anchor=(1.05, 0), loc='lower left', borderaxespad=0.)
         # it will place the legend on the outer right-hand side of the last axes

show()