You can try the -hc option in g_clustsize which gives you the distribution of clusters (monomers, dimers, trimers etc) in your trajectory as a whole.
To get the number of clusters over time you do not want to look at your trajectory in its entirety but instead divide it into time windows of appropriate size and run g_clustsize on each of them. The -hc (histo-clust.time1.xvg, histo-clust.time2.xvg ...) outputs gives you the information you need.
Oh, and keep in mind that for molecular clustering the index file option is not available, meaning g_clustsize will consider every molecule in your trajectory which might not exactly be what you want.
To focus on a subset of specific molecules you have to first generate a new tpr file containing only the molecules of interest and then run g_clustsize on that tpr.
For the extraction of the total number of clusters over time I followed Christian's approach. More analytically, I divided the trajectory in time windows, where each time window commences from 0 ns.
Instead of g_clustsize, I employed g_cluster. Because my study system is the unfolding of a biomolecule, which remains intact during the whole simulation. The clustering applies to the different conformations of the biomolecule during MD (clustered based on their RMSd).
For this purpose I wrote a "raw" python script, which is attached to this message.
I have to say that I feel a bit sorry, because I have not managed to find the "optimum path" using the output of the -clid (from g_cluster).
You probably want g_clustsize rather than g_cluster. g_cluster is for doing RMSD-based clustering of structures, not computing properties of multiple molecules that form clusters. You're likely getting the error due to differences in RMSD never meeting the threshold.