diff --git a/sample/scale.cpp b/sample/scale.cpp
new file mode 100644
index 0000000..80f2bb3
--- /dev/null
+++ b/sample/scale.cpp
@@ -0,0 +1,20 @@
+#include "yorick.hpp"
+
+using namespace CosmoTool;
+
+int main(int argc, char **argv)
+{
+  uint32_t *dimList;
+  uint32_t dimRank;
+  ProgressiveInput<float> input = ProgressiveInput<float>::loadArrayProgressive("displacement.nc", dimList, dimRank);
+  ProgressiveOutput<float> output = ProgressiveOutput<float>::saveArrayProgressive("scaledDispl.nc", dimList, dimRank);
+
+  uint32_t N = 3 * dimList[1] * dimList[2] * dimList[3];
+
+  for (uint32_t i = 0; i < N; i++)
+   {
+     output.put(input.read() * 500./.65);
+   } 
+
+  return 0;
+}
diff --git a/sample/testBQueue.cpp b/sample/testBQueue.cpp
new file mode 100644
index 0000000..9187e24
--- /dev/null
+++ b/sample/testBQueue.cpp
@@ -0,0 +1,33 @@
+#include <iostream>
+#include "bqueue.hpp"
+
+using namespace std;
+
+int main(int argc, char **argv)
+{
+  CosmoTool::BoundedQueue<int,int,100000> bq(4);
+  
+  for (int i = 10; i >= 0; i--)
+    {
+      bq.push(i, i);
+      
+      int *prio = bq.getPriorities();
+      for (int j = 0; j < 4; j++)
+	cout << prio[j] << " ";
+
+      cout << endl;
+    }
+
+  for (int i = 1; i >= -2; i--)
+    {
+      bq.push(i, i);
+      
+      int *prio = bq.getPriorities();
+      for (int j = 0; j < 4; j++)
+	cout << prio[j] << " ";
+
+      cout << endl;
+    }
+
+  return 0;
+}
diff --git a/sample/testkd.cpp b/sample/testkd.cpp
new file mode 100644
index 0000000..a3c76e7
--- /dev/null
+++ b/sample/testkd.cpp
@@ -0,0 +1,94 @@
+#include <ctime>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include "mykdtree.hpp"
+
+#define NTRY 5000000
+#define ND 2
+
+using namespace std;
+using namespace CosmoTool;
+
+typedef KDTree<ND,char> MyTree;
+typedef KDCell<ND,char> MyCell;
+
+MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells)
+{
+    MyCell *near2 = 0;
+    double R2 = INFINITY;
+    for (int i = 0; i < Ncells; i++)
+      {
+	double d2 = 0;
+	for (int j = 0; j < ND; j++)
+	  {
+	    double delta = xc[j]-cells[i].coord[j];
+	    d2 += delta*delta;
+	  }
+	if (d2 < R2)
+	  {
+	    near2 = &cells[i];
+	    R2 = d2;
+	  }
+      }
+    return near2;
+}
+
+int main()
+{
+  uint32_t Ncells = 100000;
+  MyCell *cells = new MyCell[Ncells];
+  
+  for (int i = 0; i < Ncells; i++)
+    {
+      cells[i].active = true;
+      for (int l = 0; l < ND; l++)
+	cells[i].coord[l] = drand48();
+    }
+  
+  MyTree tree(cells, Ncells);
+  
+  MyTree::coords *xc = new MyTree::coords[NTRY];
+
+  cout << "Generating seeds..." << endl;
+  for (int k = 0; k < NTRY; k++)
+    {
+      for (int l = 0; l < ND; l++)
+	xc[k][l] = drand48();
+    }
+
+  // Check consistency
+  cout << "Check consistency..." << endl;
+#if 0
+  for (int k = 0; k < NTRY; k++) {
+    MyCell *near = tree.getNearestNeighbour(xc[k]);
+    MyCell *near2 = findNearest(xc[k], cells, Ncells);
+    assert(near == near2);
+  }
+#endif  
+  cout << "Check timing..." << endl;
+  // Check timing
+  clock_t startTimer = clock();
+  
+  for (int k = 0; k < NTRY; k++) {
+    MyCell *near = tree.getNearestNeighbour(xc[k]);
+  }
+
+  clock_t endTimer = clock();
+  clock_t delta = endTimer-startTimer;
+  double myTime = delta*1.0/CLOCKS_PER_SEC * 1000000.0;
+
+  cout << "KDTree search/sec = " << myTime/NTRY << " us" << endl;
+  
+  startTimer = clock();
+  for (int k = 0; k < NTRY; k++) {
+    MyCell *near = findNearest(xc[k], cells, Ncells);
+  }
+  endTimer = clock();
+  delta = endTimer-startTimer;
+  myTime = delta*1.0/CLOCKS_PER_SEC * 1000000.0;
+
+  cout << "Direct search/sec = " << myTime/NTRY << " us" << endl;
+  
+  return 0;
+}
diff --git a/sample/testkd2.cpp b/sample/testkd2.cpp
new file mode 100644
index 0000000..0688200
--- /dev/null
+++ b/sample/testkd2.cpp
@@ -0,0 +1,78 @@
+#include <ctime>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include "mykdtree.hpp"
+
+#define NTRY 3
+#define ND 2
+
+using namespace std;
+using namespace CosmoTool;
+
+typedef KDTree<ND,char> MyTree;
+typedef KDCell<ND,char> MyCell;
+
+MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells)
+{
+    MyCell *near2 = 0;
+    double R2 = INFINITY;
+    for (int i = 0; i < Ncells; i++)
+      {
+	double d2 = 0;
+	for (int j = 0; j < ND; j++)
+	  {
+	    double delta = xc[j]-cells[i].coord[j];
+	    d2 += delta*delta;
+	  }
+	if (d2 < R2)
+	  {
+	    near2 = &cells[i];
+	    R2 = d2;
+	  }
+      }
+    return near2;
+}
+
+int main()
+{
+  uint32_t Ncells = 100000;
+  MyCell *cells = new MyCell[Ncells];
+  
+  for (int i = 0; i < Ncells; i++)
+    {
+      cells[i].active = true;
+      for (int l = 0; l < ND; l++)
+	cells[i].coord[l] = drand48();
+    }
+  
+  MyTree tree(cells, Ncells);
+  
+  MyTree::coords *xc = new MyTree::coords[NTRY];
+
+  cout << "Generating seeds..." << endl;
+  for (int k = 0; k < NTRY; k++)
+    {
+      for (int l = 0; l < ND; l++)
+	xc[k][l] = drand48();
+    }
+
+  // Check consistency
+  cout << "Check consistency..." << endl;
+  MyCell **ngb = new MyCell *[12];
+  
+  for (int k = 0; k < NTRY; k++) {
+    cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl;
+    tree.getNearestNeighbours(xc[k], 12, ngb);
+
+    for (uint32_t i = 0; i < 12; i++)
+      {
+	double d2 = 0;
+	for (int l = 0; l < 3; l++)
+	  d2 += ({double delta = xc[k][l] - ngb[i]->coord[l]; delta*delta;});
+	cout << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << endl;
+      }
+  }
+  
+  return 0;
+}