diff --git a/src/grid.c b/src/grid.c
new file mode 100644
index 0000000000000000000000000000000000000000..3e60bbcc85b4bb840a00f3c4e68327af9d50deb5
--- /dev/null
+++ b/src/grid.c
@@ -0,0 +1,295 @@
+/** @file grid.c
+ *  @brief General grid functions common for cdr and poisson grids.
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "grid.h"
+#include "parameters.h"
+#include "proto.h"
+#include "species.h"
+
+/** @brief Gets the maximum depth of a tree (the level of its deepest node) */
+int
+grid_max_depth_r (grid_t *grid)
+{
+  int l, lc;
+  grid_t *child;
+
+  l = grid->level;
+  iter_childs (grid, child) {
+    lc = grid_max_depth_r (child);
+    l = (l > lc)? l: lc;
+  }
+  return l;
+}
+
+/** @brief Gets the minimum r of a grid and all its descendants */
+double
+grid_rmin_r (grid_t *grid)
+{
+  int lev;
+  double minr, cminr;
+  grid_t *child;
+
+  lev = grid->level;
+  minr = r_at (grid->r0, lev);
+
+  iter_childs (grid, child) {
+    cminr = grid_rmin_r (child);
+    minr = (minr > cminr)? cminr: minr;
+  }
+  return minr;
+}
+
+/** @brief Recursively prints one grid and all its descendants. */
+void
+grid_print_r (grid_t *grid, int indent)
+{
+  int i;
+  grid_t *child;
+
+  for (i = 0; i < indent; i++) fputs ("  ", stdout);
+  printf ("grid (" grid_printf_str") {\n", grid_printf_args(grid));
+
+  iter_childs (grid, child) {
+    grid_print_r (child, indent + 1);
+  }
+  for (i = 0; i < indent; i++) fputs ("  ", stdout);
+  fputs ("}\n", stdout);
+}
+
+/** @brief Checks whether a point is in the interior or on a given boundary
+ * of a grid.
+ *
+ * Does the grid contain the point (i, j) (in its own coordinates)?\n
+ * @a check has the following meaning:\n
+ *  GRID_INSIDE : Checks if the point is in the interior of the grid\n
+ *  BND_MASK_*  : Checks if it is on the given boundary (BND_MASK_ALL checks for
+ *               any boundary).\n
+ */
+int
+grid_contains (grid_t *grid, int i, int j, int check)
+{
+  int r0, r1, z0, z1;
+  int r;
+
+  r0 = grid->r0;
+  z0 = grid->z0;
+  r1 = grid->r1 - 1;
+  z1 = grid->z1 - 1;
+
+  r = FALSE;
+
+  if (check & BND_MASK_LEFT) {
+    r0 -= 1;
+  }
+
+  if (check & BND_MASK_RIGHT) {
+    r1 += 1;
+  }
+
+  if (check & BND_MASK_BOTTOM) {
+    z0 -= 1;
+  }
+
+  if (check & BND_MASK_TOP) {
+    z1 += 1;
+  }
+
+  if (r0 < i && i < r1 && z0 < j && j < z1) 
+    r = TRUE;
+
+  if (check & GRID_INSIDE) {
+    return r;
+  } else {
+    return (r && !grid_contains (grid, i, j, GRID_INSIDE));
+  }
+}
+
+/** @brief Checks whether two grids overlap and return the rectangle that they
+ * in common.
+ *
+ * Checks whether two grids overlap and return the rectangle that they 
+ * have in common in the coordinates of the finest grid in 
+ * (*left, *bottom) - (*right, *top)\n
+ *
+ * Note that *left, *bottom... can be modified even if the two grids do
+ * not overlap.
+ * @a buff1 and @a buff2 are integers that specify the numbers of cells
+ * that are taken as buffer in the boundaries for @a grid1 and @a grid2.
+ */
+int
+grid_overlap (grid_t *grid1, grid_t *grid2, int buff1, int buff2, 
+	      int *left, int *bottom, int *right, int *top, int *level_diff)
+{
+  return grid_overlap_with_shifts (grid1, grid2, buff1, buff2, 
+				  left, bottom, right, top, level_diff,
+				  0, 0);
+}
+
+/** @brief Checks the overlap between two grids.
+ *
+ * Checks the overlap between two grids before it shifts the
+ * coarser grid by (shift_r << (level_diff - 1), shift_z << (level_diff - 1)
+ * (in finer grid coordinates).
+ *
+ *  Q: WTF! WHY!!!? \n
+ *  A: We use this to implement the mapping between electric fields.  Since
+ *     The electric field is evaluated not in the cell centers but in the
+ *     borders, when we interpolate we have to take a shift into account.
+ *     See mapper.c and mapper.h for more info.
+ */
+int
+grid_overlap_with_shifts (grid_t *grid1, grid_t *grid2, int buff1, int buff2, 
+			 int *left, int *bottom, int *right, int *top, 
+			 int *level_diff, int shift_r, int shift_z)
+{
+  int sign;
+
+  debug (3, "grid_overlap(" grid_printf_str", " grid_printf_str 
+	 ", buff1 = %d, buff2 = %d)\n",
+	 grid_printf_args(grid1), grid_printf_args(grid2),
+	 buff1, buff2);
+
+  sign = 1;
+
+  /* If level_diff < 0, grid2 is finer */
+  *level_diff = grid1->level - grid2->level;
+
+  if (*level_diff < 0) {
+    sign = -1;
+    *level_diff = - *level_diff;
+    XCHG (grid1, grid2);
+    XCHG (buff1, buff2);
+  }
+
+  /* At this point, level_diff >= 0 and grid1 is finer 
+     or equivalent to grid2 */
+
+  *left = MAX_AT_LEVEL_WITH_SHIFT (grid1->r0 - buff1, grid2->r0 - buff2, 
+				   *level_diff, shift_r);
+  *right = MIN_AT_LEVEL_WITH_SHIFT (grid1->r1 + buff1, grid2->r1 + buff2, 
+				    *level_diff, shift_r);
+
+  /* Note: *right is not contained in the grid.  Hence, if *left == *right
+     the two grids do not overlap.  Same for *bottom and *top, below. */
+  debug (3, "*left = %d, *right = %d\n", *left, *right);
+  if (*left >= *right) {
+    *level_diff *= sign;
+    return FALSE;
+  }
+
+  *bottom = MAX_AT_LEVEL_WITH_SHIFT (grid1->z0 - buff1, grid2->z0 - buff2, 
+				     *level_diff, shift_z);
+  *top = MIN_AT_LEVEL_WITH_SHIFT (grid1->z1 + buff1, grid2->z1 + buff2, 
+				  *level_diff, shift_z);
+
+  debug (3, "*bottom = %d, *top = %d\n", *bottom, *top);
+  if (*bottom >= *top) {
+    *level_diff *= sign;
+    return FALSE;
+  }
+  debug (3, "pois_overlap = TRUE\n");
+
+  *level_diff *= sign;
+  return TRUE;
+}
+
+/** @brief Makes a grid inherit the ext_bound field of its parent,
+ *
+ * with the appropiate modifications.
+ */
+void
+grid_inherit_ext_bound (grid_t *grid)
+{
+  grid_t *parent = grid->parent;
+
+  debug (3, "grid_inherit_ext_bound (...)\n");
+
+  grid->ext_bound = BND_NONE;
+
+  /* If the parent does not have any external boundary, the child will also
+     not have */
+  if (BND_NONE == parent->ext_bound) {
+    return;
+  }
+
+  if (parent->ext_bound & BND_MASK_LEFT) {
+    if ((parent->r0 << 1) == grid->r0) grid->ext_bound |= BND_MASK_LEFT;
+  }
+
+  if (parent->ext_bound & BND_MASK_RIGHT) {
+    if ((parent->r1 << 1) == grid->r1) grid->ext_bound |= BND_MASK_RIGHT;
+  }
+
+  if (parent->ext_bound & BND_MASK_BOTTOM) {
+    if ((parent->z0 << 1) == grid->z0) grid->ext_bound |= BND_MASK_BOTTOM;
+  }
+
+  if (parent->ext_bound & BND_MASK_TOP) {
+    if ((parent->z1 << 1) == grid->z1) grid->ext_bound |= BND_MASK_TOP;
+  }
+
+}
+
+
+/** @brief Finds the finest descendant of grid that contains the point given by
+ * the (@a r, @a z) coordinates.
+ *
+ * Or NULL if the point is not contained in grid.
+ * We assume that only ONE sub-grid can contain a given point.
+ */
+grid_t *
+grid_finest_containing_r (grid_t *grid, double r, double z)
+{
+  grid_t *leaf, *ret;
+  
+  if (r < er_r_at(grid->r0 - 1, grid->level) || 
+      r > er_r_at(grid->r1 - 1, grid->level))
+    return NULL;
+
+  if (z < ez_z_at(grid->z0 - 1, grid->level) || 
+      z > ez_z_at(grid->z1 - 1, grid->level))
+    return NULL;
+
+  iter_childs (grid, leaf) {
+    ret = grid_finest_containing_r (leaf, r, z);
+    if (NULL != ret)
+      return ret;
+  }
+  return grid;
+}
+
+/** @brief How many children does one grid have? */
+int
+grid_howmany_children (grid_t *grid)
+{
+  grid_t *leaf;
+  int count;
+
+  count = 0;
+  iter_childs (grid, leaf) {
+    count ++;
+  }
+  return count;
+}
+
+/** @brief Finds child number @a n of the given grid.
+ *
+ * If it doesn't have enough children just returns NULL
+ */
+grid_t*
+grid_get_child (grid_t *grid, int n)
+{
+  grid_t *leaf;
+  int count;
+
+  count = 0;
+  iter_childs (grid, leaf) {
+    if (count == n) return leaf;
+    count ++;
+  }
+  return NULL;
+}