diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp
index d6ef31d..365e836 100644
--- a/c_tools/zobov2/jozov2_watershed.cpp
+++ b/c_tools/zobov2/jozov2_watershed.cpp
@@ -2,6 +2,7 @@
 #include <omp.h>
 #endif
 
+#include <queue>
 #include <set>
 #include <iostream>
 #include <fstream>
@@ -13,196 +14,195 @@
 using namespace std;
 using boost::format;
 
+struct ZoneDensityPair
+{
+  int h;
+  double density;
+  double core;
+
+  bool operator<(const ZoneDensityPair& p2) const
+  {
+    return (density > p2.density) || (density==p2.density && core > p2.core);
+  }
+};
+
+typedef priority_queue<ZoneDensityPair> ZoneQueue;
+
+static void build_process_queue(ZoneQueue& q, ZONE *z, PARTICLE *p, char *inyet, int h)
+{
+  ZoneDensityPair zdp;
+  ZONE& z_h = z[h];
+  bool interior = true;
+
+  assert(inyet[h] == 1);
+  for (int za = 0; za < z_h.nadj; za++)
+   {
+     zdp.h = z_h.adj[za];
+     zdp.density = z_h.slv[za];
+     zdp.core = p[z[zdp.h].core].dens;
+     if (inyet[zdp.h] == 0)
+       {
+	 q.push(zdp);
+	 interior = false;
+       }
+   }
+  if (interior)
+    inyet[h] = 2;
+}
+
 void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol)
 {
-  int *iord;
-
-  double *sorter = new double[numZones+1];
-  /* Assign sorter by probability (could use volume instead) */
-  for (int h = 0; h < numZones; h++)
-    sorter[h] = (double)z[h].core;
-    
   /* Text output file */
 
-  printf("about to sort (pre-watershed sort) ...\n");FF;
-
-  iord = new int[numZones];
- 
-  findrtop(sorter, numZones, iord, numZones);
-  delete[] sorter;
-
 #pragma omp parallel
   {
     char *inyet, *inyet2;
     int *zonelist, *zonelist2;
     int nhl;
-    int *links = new int[NLINKS];
-    bool *done_zones;
     int prev_ii = -1;
 
     inyet = new char[numZones];
     inyet2 = new char[numZones];
     zonelist = new int[numZones];
     zonelist2 = new int[numZones];
-    done_zones = new bool[numZones];
 
     fill(inyet, inyet + numZones, 0);
     fill(inyet2, inyet2 + numZones, 0);
-    fill(done_zones, done_zones + numZones, false);
 
     nhl = 0;
 #pragma omp for schedule(dynamic,1)
     for (int h = 0; h < numZones; h++)
       {
         int nhlcount = 0;
-        float lowvol;
+        float previous_lowvol = BIGFLT, lowvol, z_cur_core_dens;
         bool beaten;
-        set<int> to_process;
+        priority_queue<ZoneDensityPair> to_process;
+	int link0;
+        int save_nhl;
+        pid_t save_npjoin;
 
         for (int hl = 0; hl < nhl; hl++)
           inyet[zonelist[hl]] = 0;
         
         zonelist[0] = h;
         inyet[h] = 1;
-        nhl = 1;
-        to_process.insert(h);
-        z[h].npjoin = z[h].np;
+        save_nhl = nhl = 1;
+        save_npjoin = z[h].npjoin = z[h].np;
+        z_cur_core_dens = p[z[h].core].dens;
+
+        build_process_queue(to_process, z, p, inyet, h);
+	
         do {
           /* Find the lowest-volume (highest-density) adjacency */
-          int nl = 0;
-          
           beaten = false;
-          lowvol = BIGFLT;
           
-          set<int>::iterator iter = to_process.begin();
-          while (iter != to_process.end())
-            {
-              int h2 = *iter;
-              ZONE& z_cur = z[h2];
-              bool interior = true, touched = false;
-              assert(inyet[h2] == 1);
-
-              for (int za = 0; za < z_cur.nadj; za++)
-                {
-                 if (inyet[z_cur.adj[za]] == 0)
-                    {
-                      interior = false;
-                      if (z_cur.slv[za] == lowvol)
-                        {
-                          links[nl] = z_cur.adj[za];
-                          touched = true;
-                          nl ++;
-                          if (nl == NLINKS) 
-                            {
-                              printf("Too many links with the same rho_sl!  Increase NLINKS from %d\n",nl);
-                              exit(0);
-                            }
-                        }
-                      else
-                        if (z_cur.slv[za] < lowvol)
-                          {
-                            lowvol = z_cur.slv[za];
-                            links[0] = z_cur.adj[za];
-                            nl = 1;
-                            touched = true;
-                          }
-                    }
-                }
-              ++iter;
-              if (interior || !touched)
-                to_process.erase(h2);
-              if (interior)
-                inyet[h2] = 2; /* No bordering exter. zones */
-            }
-          
-          if (nl == 0)
+          if (to_process.empty())
             {
               beaten = true;
               z[h].leak = maxvol;
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
               continue;
             }
-          
-          if (lowvol > voltol)
+
+	  do
+	    {
+	      lowvol = to_process.top().density;
+	      link0 = to_process.top().h; 
+	      to_process.pop();
+	    }
+	  while ((inyet[link0] != 0) && (!to_process.empty()));
+
+	  if (to_process.empty())
+	    {
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+	      beaten = true;
+	      z[h].leak = maxvol;
+	      continue;
+	    }
+
+          /* See if there's a beater */
+          if (previous_lowvol != lowvol)
+            {
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+              previous_lowvol = lowvol;
+            }
+
+	  if (lowvol > voltol)
             {
               beaten = true;
               z[h].leak = lowvol;
               continue;
             }
           
-          for (int l = 0; l < nl; l++)
-            {
-              //if (!done_zones[links[l]]) /* Equivalent to p[z[links[l]].core].dens < p[z[h].core].dens) as zones are sorted. */
-              if (p[z[links[l]].core].dens < p[z[h].core].dens)
-                beaten = true;
-            }
-          
-          if (beaten)
+          if (p[z[link0].core].dens < z_cur_core_dens)
             {
+	      beaten = true;
               z[h].leak = lowvol;
               continue;
             }
           
           /* Add everything linked to the link(s) */
           int nhl2 = 0;
-          for (int l = 0; l < nl; l++)
-            {
-              if (inyet2[links[l]] == 0)
-                {
-                  zonelist2[nhl2] = links[l];
-                  inyet2[links[l]] = 1;
-                  nhl2 ++;
-                  bool added = true;
-                  while (added && !beaten)
-                    {
-                      added = false;
-                      for (int hl = 0; (hl < nhl2) && (!beaten); hl++)
-                        {
-                          int h2 = zonelist2[hl];
-                          
-                          if (inyet2[h2] == 1) {
-                            bool interior = true; /* Guilty until proven innocent */
-                            for (int za = 0; za < z[h2].nadj; za ++) {
-                              int link2 = z[h2].adj[za];
-                              if ((inyet[link2]+inyet2[link2]) == 0) {
-                                interior = false;
-                                if (z[h2].slv[za] <= lowvol) {
-                                  //if (!done_zones[link2]) { // Equivalent to p[z[link2].core].dens < p[z[h].core].dens)
-                                  if (p[z[link2].core].dens < p[z[h].core].dens) {
-                                    beaten = true;
-                                    break;
-                                  }
-                                  zonelist2[nhl2] = link2;
-                                  inyet2[link2] = 1;
-                                  nhl2++;
-                                  added = true;
-                                }
-                              }
-                            }
-                            if (interior)
-                              inyet2[h2] = 2;
-                          }
-                        }
-                    }
-                }
-            }
-          
-          /* See if there's a beater */
-          if (beaten) {
-            z[h].leak = lowvol;
-          } else {
-            for (int h2 = 0; h2 < nhl2; h2++) {
-              int new_h = zonelist2[h2];
+	  zonelist2[0] = link0;
+	  inyet2[link0] = 1;
+	  nhl2=1;
+	  
+	  bool added = true;
+	  while (added && !beaten)
+	    {
+	      added = false;
+	      for (int hl = 0; (hl < nhl2) && (!beaten); hl++)
+		{
+		  int h2 = zonelist2[hl];
+                  
+		  if (inyet2[h2] == 1) {
+		    bool interior = true; /* Guilty until proven innocent */
+                    
+		    for (int za = 0; za < z[h2].nadj; za ++) {
+		      int link2 = z[h2].adj[za];
+		      
+		      if ((inyet[link2]+inyet2[link2]) == 0) {
+			interior = false;
+			if (z[h2].slv[za] <= lowvol) {
+			  if (p[z[link2].core].dens < z_cur_core_dens) {
+			    beaten = true;
+			    break;
+			  }
+			  zonelist2[nhl2] = link2;
+			  inyet2[link2] = 1;
+			  nhl2++;
+			  added = true;
+			}
+		      }
+		    }
+		    
+		    if (interior)
+		      inyet2[h2] = 2;
+		  }
+		}
+	    }
+	  
+	  /* See if there's a beater */
+	  if (beaten) {
+	    z[h].leak = lowvol;
+	  } else {
 
+	    for (int h2 = 0; h2 < nhl2; h2++) {
+	      int new_h = zonelist2[h2];
+	      
               zonelist[nhl] = new_h;
               assert(inyet[new_h] == 0);
+              z[h].npjoin += z[new_h].np;
               inyet[new_h] = 1;
               if (inyet2[new_h] != 2)
-                to_process.insert(new_h);
+                build_process_queue(to_process, z, p, inyet, new_h);
               nhl++;
-              z[h].npjoin += z[new_h].np;
             }
           }
-
+	  
           for (int hl = 0; hl < nhl2; hl++)
             inyet2[zonelist2[hl]] = 0;
           if (nhl/10000 > nhlcount) {
@@ -212,10 +212,18 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo
               (cout << format(" %d [%d]") % nhl % to_process.size()).flush();
             nhlcount = nhl/10000;
           }
-        } while((lowvol < BIGFLT) && (!beaten));
-        
-        z[h].denscontrast = z[h].leak/p[z[h].core].dens;
-        if (z[h].denscontrast < 1.) z[h].denscontrast = 1.;
+	}
+	while((lowvol < BIGFLT) && (!beaten));
+
+        if (!beaten)
+          { 
+            save_npjoin = z[h].npjoin;
+            save_nhl = nhl;
+          }
+	
+	z[h].denscontrast = z[h].leak/p[z[h].core].dens;
+        if (z[h].denscontrast < 1.) 
+	  z[h].denscontrast = 1.;
         
         
         /* Don't sort; want the core zone to be first */
@@ -225,27 +233,24 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo
           FF;
         }
         /* Calculate volume */
+        z[h].npjoin = save_npjoin;
         z[h].voljoin = 0.;
-        z[h].zonelist = new int[nhl];
-        z[h].numzones = nhl;
-        for (int q = 0; q < nhl; q++) {
+        z[h].zonelist = new int[save_nhl];
+        z[h].numzones = save_nhl;
+        for (int q = 0; q < save_nhl; q++) {
           z[h].voljoin += z[zonelist[q]].vol;
           z[h].zonelist[q] = zonelist[q];
         }
         
-        done_zones[h] = true;
-        z[h].nhl = nhl;
+        z[h].nhl = save_nhl;
       }
     delete[] zonelist;
     delete[] zonelist2;
-    delete[] links;
     delete[] inyet;
     delete[] inyet2;
-    delete[] done_zones;
-
+    
   }
-  delete[] iord;
-
+  
   double maxdenscontrast = 0;
 #pragma omp parallel shared(maxdenscontrast)
   {
diff --git a/c_tools/zobov2/jozov2_zones.cpp b/c_tools/zobov2/jozov2_zones.cpp
index ca59768..ebfc974 100644
--- a/c_tools/zobov2/jozov2_zones.cpp
+++ b/c_tools/zobov2/jozov2_zones.cpp
@@ -39,10 +39,13 @@ void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped,
   }
   (cout << "Post-jump ..." << endl).flush();
   
-  numZones = 0;
+  pid_t loc_NumZones = 0;
+#pragma omp parallel for schedule(static) reduction(+:loc_NumZones)
   for (pid_t i = 0; i < np; i++)
     if (numinh[i] > 0) 
-      numZones++;
+      loc_NumZones++;
+
+  numZones = loc_NumZones;
   cout << format("%d initial zones found") % numZones << endl;
 
   delete[] jumper;
diff --git a/zobov/voztie.c b/zobov/voztie.c
index 4aaaf57..5a15448 100644
--- a/zobov/voztie.c
+++ b/zobov/voztie.c
@@ -2,8 +2,21 @@
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include "voz.h"
 
+int compare_int(const void *a, const void *b)
+{
+  const int *ia = (const int *)a;
+  const int *ib = (const int *)b;
+  if ((*ia) < (*ib))
+    return -1;
+  else if ((*ia) > (*ib))
+    return 1;
+  else
+    return 0;
+}
+
 int main(int argc, char *argv[]) {
 
   FILE *part, *adj, *vol;
@@ -135,15 +148,25 @@ int main(int argc, char *argv[]) {
 	fread(&na,1,sizeof(int),part);
         pid_t pid = orig[p];
 	if (na > 0) {
-          assert(adjs[pid].nadj == 0);// || adjs[pid].nadj == na);
-	  adjs[orig[p]].nadj = na;
-	  if (adjs[pid].adj == 0)
-            adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t));
-	  if (adjs[orig[p]].adj == 0) {
+          pid_t *previous_adj = adjs[pid].adj;
+          assert(adjs[pid].nadj == 0 || adjs[pid].nadj == na);
+	  adjs[pid].nadj = na;
+          adjs[pid].adj = (pid_t *)malloc(na*sizeof(pid_t));
+	  if (adjs[pid].adj == 0) {
 	    printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
 	    exit(0);
 	  }
-	  fread(adjs[orig[p]].adj,na,sizeof(pid_t),part);
+	  fread(adjs[pid].adj,na,sizeof(pid_t),part);
+          if (previous_adj != 0)
+           {
+             qsort(previous_adj, na, sizeof(pid_t), compare_int);
+             qsort(adjs[pid].adj, na, sizeof(pid_t), compare_int);
+             if (memcmp(previous_adj, adjs[pid].adj, sizeof(pid_t)*na) != 0)
+               {
+                 printf("Inconsistent adjacencies between two divisions.\n");
+                 abort();
+               }
+           }
 	} else {
 	  printf("0"); fflush(stdout);
 	}