diff --git a/include/rz_array.h b/include/rz_array.h
new file mode 100644
index 0000000000000000000000000000000000000000..5cd4b08809e2566805fedd09cb3c0e7a349334ba
--- /dev/null
+++ b/include/rz_array.h
@@ -0,0 +1,86 @@
+/** @file rz_array.h
+ *  @brief Definitions for 2D/3D arrays
+ */
+
+#ifndef _RZ_ARRAY_H_
+
+#include "cstream.h"
+
+typedef struct rz_array_t rz_array_t;
+
+struct rz_array_t {
+  REAL *data;
+ 
+  int strides[3];
+  int r0, z0, theta0, nr, nz;
+  int ntheta;
+  int dim;
+
+  int len;
+
+  rz_array_t *host;
+};
+
+#define R_INDX 0
+#define Z_INDX 1
+#define THETA_INDX 2
+
+/* The 3D macros */
+#define RZTP_OK(_A, _R, _Z, _T)						\
+  ((_R) >= (_A)->r0 && (_R) < (_A)->nr + (_A)->r0 + 4			\
+   && (_Z) >= (_A)->z0 && (_Z) < (_A)->nz + (_A)->z0 + 4		\
+   && (_T) >= (_A)->theta0 && (_T) < (_A)->ntheta + (_A)->theta0 + 4)	\
+
+#define __RZTP(_A, _R, _Z, _T) ((_A)->data				\
+			      +	((_R) - (_A)->r0) * (_A)->strides[R_INDX] \
+			      + ((_Z) - (_A)->z0) * (_A)->strides[Z_INDX] \
+			      + (((_T) - (_A)->theta0) * (_A)->strides[THETA_INDX]))
+
+#if defined (DEBUG_LEVEL) && DEBUG_LEVEL > 4
+/* Look, mum! I am also doing array bound checking.
+ * But, alas, since in some parts of the code we perform direct pointer
+ * arithmetic, this check will not detect all possible out-of-indexes
+ * (but we catch most of them).
+ */
+# define RZTP(_A, _R, _Z, _T) (RZTP_OK (_A, _R, _Z, _T)?		\
+			       __RZTP (_A, _R, _Z, _T):			\
+			       (fprintf (stderr, \
+					 "%s:%d: Out of bounds ir = %d, iz = %d, itheta = %d\n", \
+					 __FILE__, __LINE__, _R, _Z, _T), \
+				fprintf (stderr, \
+					 "->r0 = %d ->z0 = %d ->theta0 = %d ->nr = %d ->nz = %d ->ntheta = %d\n", \
+					 (_A)->r0, (_A)->z0, (_A)->theta0, \
+					 (_A)->nr, (_A)->nz, (_A)->ntheta), \
+				exit(-1), (double*) NULL))
+#else
+# define RZTP(_A, _R, _Z, _T) __RZTP(_A, _R, _Z, _T)
+#endif
+
+#define RZT(_A,_R,_Z,_T) (*RZTP(_A, _R, _Z, _T))
+#define RZTm(_A,_R,_Z,_T) ((fprintf (stderr, \
+					 "%s:%d: CHECK ir = %d, iz = %d, itheta = %d\n", \
+					 __FILE__, __LINE__, _R, _Z, _T), \
+				fprintf (stderr, \
+					 "->r0 = %d ->z0 = %d ->theta0 = %d ->nr = %d ->nz = %d ->ntheta = %d\n", \
+					 (_A)->r0, (_A)->z0, (_A)->theta0, \
+					 (_A)->nr, (_A)->nz, (_A)->ntheta), \
+				fprintf (stderr, \
+					 "->data = %g\n", \
+					 (_A)->data), \
+				exit(-1), (double*) NULL))
+
+
+/* These are valid for 2D arrays.  When applied to a 3D array, they
+ * return the value with theta=0.
+ */
+#define RZP(_A,_R,_Z) RZTP(_A, _R, _Z, (_A)->theta0)
+#define RZ(_A,_R,_Z) RZT(_A, _R, _Z, (_A)->theta0)
+
+#define BND_CND_HNEUMANN   1
+#define BND_CND_HDIRICHLET -1
+
+#define BND_INWARD -1
+#define BND_OUTWARD 1
+
+#define _RZ_ARRAY_H_
+#endif