icpu =-1;
  int iatoms=0;
  for(int i=0;i<NCPUS;i++) ci[i] = i;

  int starti,startj;
  int ii,jj;
  for(ii=0;ii<Dimension;ii++){
    for(jj=0;jj<Dimension;jj++){  
      icpu++;
      if(icpu > NCPUS-1){
	break;
      }
      CellCpus[ii][jj][0] = icpu;
      cc[icpu] += CellCpus[ii][jj][3];
      isum++;
    }
    if(icpu > NCPUS-1) {
      break;
    }
  }
  starti=ii;
  startj=jj;
  // Now quick sort the cc and ci arrays so that I can put cells in the lowest

  int Limit = NumParticles/(NCPUS); // target size
  while(isum < TotalCells){
    sort(cc,ci,0,NCPUS);
    for(int i=0;i<NCPUS;i++)
    icpu = -1;
    for(ii=starti; ii<Dimension;ii++){
      for(jj=startj;jj<Dimension;jj++){
	icpu++;
	while(cc[icpu]+CellCpus[ii][jj][3] > Limit && icpu !=NCPUS-1) {
	  icpu++;
	}
	if(icpu > NCPUS-1 || isum >= TotalCells){
	  break;
	}
	CellCpus[ii][jj][0] = ci[icpu];
	cc[icpu] += CellCpus[ii][jj][3];
	isum ++;
      }
      if(icpu > NCPUS-1 || isum >= TotalCells){
	break;
      }
      startj =0;
    }
    starti=ii;
    startj=jj;
  }
  free(cc);free(ci);

  int LocalCPUnCells = 0;
  vector <int> LocalCPUCells;
  
  for(int i=0;i<Dimension;i++){
    for(int j=0;j<Dimension;j++){
      if(MY_PE == CellCpus[i][j][0]){
	LocalCPUnCells++;
	LocalCPUCells.push_back(i);
	LocalCPUCells.push_back(j);
      }
    }
  }
  
void swap(int *a, int *b){
  int t=*a; 
  *a=*b; 
  *b=t;
}

void sort(int arr[], int arr2[], int beg, int end){
  if (end > beg + 1){
    int piv = arr[beg], l = beg + 1, r = end;
    while (l < r){
      if (arr[l] <= piv)
        l++;
      else{
	r--;
        swap(&arr[l], &arr[r]);
	swap(&arr2[l], &arr2[r]);
      }
    }
    l--;
    swap(&arr[l], &arr[beg]);
    swap(&arr2[l],&arr2[beg]);
    sort(arr, arr2, beg, l);
    sort(arr, arr2, r, end);
  }