(1) Modified the algorithm to search for continuous regions, used matlab's function bwconncomp. The use of this function reduces the code for finding connected regions to one line in case of a non-periodic domain. Unfortunately it requires some additional code in the case of a periodic domain, in which case the regions sliced apart by the meridional domain boundary have to be glued together again.

(2) Modified the algorithm to construct the sparse matrix for the least-square problem.

It seems that the  modifications speed up the computation by 10 to 20 percent, depending on resolution, and I think this is mostly due to (2).