diff --git a/example/ALL_test.cpp b/example/ALL_test.cpp
index 5827b7df32e17cd763300179e6f0ea636870d77f..9fec20ce0890b73da6fa6f6cb76e955fc4da858f 100644
--- a/example/ALL_test.cpp
+++ b/example/ALL_test.cpp
@@ -417,6 +417,7 @@ int main(int argc, char **argv) {
     // print out chosen method
     if (localRank == 0) {
       switch (chosen_method) {
+      case ALL::LB_t::TENSOR_MAX:
       case ALL::LB_t::TENSOR:
         std::cout << "chosen method: TENSOR" << std::endl;
         break;
@@ -457,6 +458,7 @@ int main(int argc, char **argv) {
     int nvertices = 2;
 
     switch (chosen_method) {
+    case ALL::LB_t::TENSOR_MAX:
     case ALL::LB_t::TENSOR:
       nvertices = 2;
       break;
@@ -488,6 +490,7 @@ int main(int argc, char **argv) {
     std::vector<ALL::Point<double>> old_vertices(nvertices, lp);
 
     switch (chosen_method) {
+    case ALL::LB_t::TENSOR_MAX:
     case ALL::LB_t::TENSOR:
       vertices.at(0) = lp;
       vertices.at(1) = up;
@@ -1113,6 +1116,7 @@ int main(int argc, char **argv) {
         std::cout << "beginning particles transfer..." << std::endl;
 #endif
       switch (chosen_method) {
+      case ALL::LB_t::TENSOR_MAX:
       case ALL::LB_t::TENSOR: {
         int n_transfer[2];
         int n_recv[2];
diff --git a/example/create_logo.cpp b/example/create_logo.cpp
index be1ec4588336f53c6924ae1f4d5baea5edbc149a..ef8674d22ec4e494aa8d6e02fd935b90e28e4573 100644
--- a/example/create_logo.cpp
+++ b/example/create_logo.cpp
@@ -12,7 +12,7 @@
 int main(int argv, char** argc)
 {
     // if applicable read in size in x
-    const double Lx = (argv >= 3)?atof(argc[1]):LX_DEFAULT;
+    const double Lx = (argv >= 2)?atof(argc[1]):LX_DEFAULT;
     // if applicable read in size in y
     const double Ly = (argv >= 3)?atof(argc[2]):LY_DEFAULT;
     // if applicable read in number of particles in z
diff --git a/include/ALL.hpp b/include/ALL.hpp
index d4a3a2e2afb51ca0b5795609878f124ebd9dd0b4..0665cff43fb64dd70b3532baf21d24783e05cb96 100644
--- a/include/ALL.hpp
+++ b/include/ALL.hpp
@@ -38,9 +38,9 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #ifdef ALL_VORONOI_ACTIVE
 #include "ALL_Voronoi.hpp"
 #endif
+#include "ALL_ForceBased.hpp"
 #include "ALL_Histogram.hpp"
 #include "ALL_Staggered.hpp"
-#include "ALL_ForceBased.hpp"
 #include <algorithm>
 #include <iomanip>
 #include <memory>
@@ -49,8 +49,8 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <string>
 #include <vector>
 
-#include <sys/stat.h>
 #include <cerrno>
+#include <sys/stat.h>
 
 #ifdef ALL_VTK_OUTPUT
 #include <limits>
@@ -90,6 +90,8 @@ enum LB_t : int {
 #endif
   /// histogram based load balancing
   HISTOGRAM = 4,
+  /// tensor based load balancing using maximum values
+  TENSOR_MAX = 5,
   /// Unimplemented load balancing NEVER SUPPORTED
   UNIMPLEMENTED = 128
 };
@@ -102,12 +104,12 @@ public:
   /// default constructor, that sets up internal data structures and initializes
   /// internal values shared between all balancing methods
   ALL()
-      :
-        loadbalancing_step(0)
+      : loadbalancing_step(0)
 #ifdef ALL_VTK_OUTPUT
-        , vtk_init(false)
+        ,
+        vtk_init(false)
 #endif
-        {
+  {
     // determine correct MPI data type for template T
     if (std::is_same<T, double>::value)
       MPIDataTypeT = MPI_DOUBLE;
@@ -157,10 +159,12 @@ public:
     case LB_t::HISTOGRAM:
       balancer.reset(new Histogram_LB<T, W>(d, std::vector<W>(10), g));
       break;
+    case LB_t::TENSOR_MAX:
+      balancer.reset(new TensorMax_LB<T, W>(d, (W)0, g));
+      break;
     default:
-      throw InvalidArgumentException(
-          __FILE__, __func__, __LINE__,
-          "Unknown type of loadbalancing passed.");
+      throw InvalidArgumentException(__FILE__, __func__, __LINE__,
+                                     "Unknown type of loadbalancing passed.");
     }
     balancer->setDimension(d);
     balancer->setGamma(g);
@@ -176,8 +180,7 @@ public:
   /// balancing step
   /// @param[in] g the value used for the correction value, if required by the
   /// chosen balancing method
-  ALL(const LB_t m, const int d, const std::vector<Point<T>> &inp,
-      const T g)
+  ALL(const LB_t m, const int d, const std::vector<Point<T>> &inp, const T g)
       : ALL(m, d, g) {
     balancer->setVertices(inp);
     calculate_outline();
@@ -208,7 +211,7 @@ public:
   void setCommunicator(const MPI_Comm comm) {
     int comm_type;
     MPI_Topo_test(comm, &comm_type);
-    if ( comm_type == MPI_CART) {
+    if (comm_type == MPI_CART) {
       communicator = comm;
       balancer->setCommunicator(communicator);
       // set procGridLoc and procGridDim
@@ -217,8 +220,8 @@ public:
       int size[dimension];
       int periodicity[dimension];
       MPI_Cart_get(communicator, dimension, size, periodicity, location);
-      procGridLoc.assign(location, location+dimension);
-      procGridDim.assign(size, size+dimension);
+      procGridLoc.assign(location, location + dimension);
+      procGridDim.assign(size, size + dimension);
       procGridSet = true;
     } else {
       int rank;
@@ -286,9 +289,7 @@ public:
 
   /// method to call the setup of the chosen balancing method (not all methods
   /// require a setup)
-  void setup() {
-    balancer->setup();
-  }
+  void setup() { balancer->setup(); }
 
   /// method the trigger the balancing step, that updates the vertices according
   /// to the previously provided work and chosen method
@@ -341,7 +342,8 @@ public:
       d_max = (double)(max_work) * (double)n_ranks / (double)sum_work - 1.0;
 
 
-      // internal needs to be readded to argument list const bool internal = false
+      // internal needs to be readded to argument list
+      const bool internal = false
       for (int i = 0; i < estimation_limit && d_max > 0.1 && !internal; ++i)
       {
           balancer->setWork(tmp_work);
@@ -390,10 +392,14 @@ public:
       balancer->balance(loadbalancing_step);
       calculate_outline();
       break;
+    case LB_t::TENSOR_MAX:
+      balancer->balance(loadbalancing_step);
+      calculate_outline();
+      break;
+
     default:
-      throw InvalidArgumentException(
-          __FILE__, __func__, __LINE__,
-          "Unknown type of loadbalancing passed.");
+      throw InvalidArgumentException(__FILE__, __func__, __LINE__,
+                                     "Unknown type of loadbalancing passed.");
     }
     loadbalancing_step++;
     return balancer->getVertices();
@@ -414,9 +420,7 @@ public:
   /// performed
   /// @result std::vector<ALL::Point<T>> containing the resulting vertices after
   /// the balancing step
-  std::vector<Point<T>> &getVertices() {
-    return balancer->getVertices();
-  };
+  std::vector<Point<T>> &getVertices() { return balancer->getVertices(); };
 
   /// get the dimension of the ALL::Point<T> objects used to describe the
   /// vertices of the domains
@@ -435,16 +439,14 @@ public:
   /// method to provide a list of the ranks of the neighbors the local domain
   /// has in all directions
   /// @result vector if neighboring ranks
-  std::vector<int> &getNeighbors() {
-    return balancer->getNeighbors();
-  }
+  std::vector<int> &getNeighbors() { return balancer->getNeighbors(); }
 
   /// method to provide a list of neighboring vertices, e.g. required for
   /// VORONOI
   /// @result vector of neighboring vertices
   /// neighboring vertices are stored in
   std::vector<T> &getNeighborVertices() {
-      return balancer->getNeighborVertices();
+    return balancer->getNeighborVertices();
   }
 
   /// method to set the size of the system, e.g. required for the HISTOGRAM
@@ -724,7 +726,6 @@ public:
       cell->SetNumberOfTuples(n_ranks);
       cell->SetName("cell id");
 
-
       for (int n = 0; n < n_ranks; ++n) {
 
 #ifdef VTK_CELL_ARRAY_V2
@@ -732,20 +733,22 @@ public:
         vtkIdType pointIds[8] = {8 * n + 0, 8 * n + 1, 8 * n + 2, 8 * n + 3,
                                  8 * n + 4, 8 * n + 5, 8 * n + 6, 8 * n + 7};
 
-        vtkIdType faces[48] = { 3, 8 * n + 0, 8 * n + 2, 8 * n + 1,
-                                3, 8 * n + 1, 8 * n + 2, 8 * n + 3,
-                                3, 8 * n + 0, 8 * n + 4, 8 * n + 2,
-                                3, 8 * n + 2, 8 * n + 4, 8 * n + 6,
-                                3, 8 * n + 2, 8 * n + 6, 8 * n + 3,
-                                3, 8 * n + 3, 8 * n + 6, 8 * n + 7,
-                                3, 8 * n + 1, 8 * n + 5, 8 * n + 3,
-                                3, 8 * n + 3, 8 * n + 5, 8 * n + 7,
-                                3, 8 * n + 0, 8 * n + 4, 8 * n + 1,
-                                3, 8 * n + 1, 8 * n + 4, 8 * n + 5,
-                                3, 8 * n + 4, 8 * n + 6, 8 * n + 5,
-                                3, 8 * n + 5, 8 * n + 6, 8 * n + 7};
-
-        unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12, faces);
+        vtkIdType faces[48] = {
+            3,         8 * n + 0, 8 * n + 2, 8 * n + 1, 
+            3,         8 * n + 1, 8 * n + 2, 8 * n + 3, 
+            3,         8 * n + 0, 8 * n + 4, 8 * n + 2, 
+            3,         8 * n + 2, 8 * n + 4, 8 * n + 6, 
+            3,         8 * n + 2, 8 * n + 6, 8 * n + 3, 
+            3,         8 * n + 3, 8 * n + 6, 8 * n + 7, 
+            3,         8 * n + 1, 8 * n + 5, 8 * n + 3, 
+            3,         8 * n + 3, 8 * n + 5, 8 * n + 7, 
+            3,         8 * n + 0, 8 * n + 4, 8 * n + 1, 
+            3,         8 * n + 1, 8 * n + 4, 8 * n + 5, 
+            3,         8 * n + 4, 8 * n + 6, 8 * n + 5, 
+            3,         8 * n + 5, 8 * n + 6, 8 * n + 7};
+
+        unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
+                                         faces);
 #else
         // define grid points, i.e. vertices of local domain
         vtkIdType pointIds[8] = {8 * n + 0, 8 * n + 1, 8 * n + 2, 8 * n + 3,
@@ -788,9 +791,8 @@ public:
                                          faces->GetPointer());
 #endif
         work->SetValue(
-            n,
-            global_vertices[n * (nVertices * balancer->getDimension() + 1) +
-                            8 * balancer->getDimension()]);
+            n, global_vertices[n * (nVertices * balancer->getDimension() + 1) +
+                               8 * balancer->getDimension()]);
         cell->SetValue(n, (T)n);
       }
       unstructuredGrid->GetCellData()->AddArray(work);
@@ -807,7 +809,7 @@ public:
       // writer->SetDataModeToBinary();
       writer->Write();
 
-      //delete[] global_vertices;
+      // delete[] global_vertices;
     }
   }
 #endif
@@ -866,25 +868,24 @@ private:
   void createDirectory(const char *filename) {
     errno = 0;
     mode_t mode = S_IRWXU | S_IRWXG; // rwx for user and group
-    if(mkdir(filename, mode)) {
-      if(errno != EEXIST) {
-        throw FilesystemErrorException(
-            __FILE__, __func__, __LINE__,
-            "Could not create output directory.");
+    if (mkdir(filename, mode)) {
+      if (errno != EEXIST) {
+        throw FilesystemErrorException(__FILE__, __func__, __LINE__,
+                                       "Could not create output directory.");
       }
     }
     // check permission in case directory already existed, but had wrong
     // permissions
     struct stat attr;
     stat(filename, &attr);
-    if( (attr.st_mode & S_IRWXU) != S_IRWXU) {
+    if ((attr.st_mode & S_IRWXU) != S_IRWXU) {
       throw FilesystemErrorException(
           __FILE__, __func__, __LINE__,
-          "Possibly already existing directory does not have correct permissions (rwx) for user");
+          "Possibly already existing directory does not have correct "
+          "permissions (rwx) for user");
     }
   }
 
-
   /// limit for the estimation procedure, which is an experimental way to use
   /// non-updated work-loads to provide better results for the final vertex
   /// positions
@@ -920,7 +921,7 @@ private:
 #endif
 };
 
-}//namespace ALL
+} // namespace ALL
 
 #endif
 
diff --git a/include/ALL_Tensor.hpp b/include/ALL_Tensor.hpp
index f43be8cbc97b0c038ca216ae25108c38cfa2ce1a..258f25dc57f036b9bda55b8bbd6f5088cae3e923 100644
--- a/include/ALL_Tensor.hpp
+++ b/include/ALL_Tensor.hpp
@@ -32,9 +32,9 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define ALL_TENSOR_HEADER_INCLUDED
 
 #include "ALL_CustomExceptions.hpp"
+#include "ALL_Defines.h"
 #include "ALL_Functions.hpp"
 #include "ALL_LB.hpp"
-#include "ALL_Defines.h"
 #include <exception>
 #include <mpi.h>
 
@@ -78,7 +78,7 @@ public:
 
   /// method to execute a load-balancing step
   /// @param[in] step number of the load-balancing step
-  virtual void balance(int step) override;
+  virtual void balance(int step) override {balance(step, MPI_SUM);}
 
   // getter for variables (passed by reference to avoid
   // direct access to private members of the object)
@@ -92,6 +92,9 @@ public:
   /// @param[in] data pointer to the data structure
   virtual void setAdditionalData(known_unused const void *data) override {}
 
+  /// method to execute a load-balancing step
+  /// @param[in] step number of the load-balancing step
+  virtual void balance(int step, MPI_Op reductionMode);
 private:
   // type for MPI communication
   MPI_Datatype MPIDataTypeT;
@@ -104,6 +107,10 @@ private:
   // list of neighbors
   std::vector<int> neighbors;
   std::vector<int> nNeighbors;
+  
+protected:
+
+
 };
 
 // setup routine for the tensor-based load-balancing scheme
@@ -173,17 +180,20 @@ template <class T, class W> void Tensor_LB<T, W>::setup() {
   }
 }
 
-template <class T, class W> void Tensor_LB<T, W>::balance(int) {
+template <class T, class W> void Tensor_LB<T, W>::balance(int step, MPI_Op reductionMode) {
   int dim = this->getDimension();
   std::vector<Point<T>> newVertices = this->vertices;
   this->prevVertices = this->vertices;
 
+  bool localIsSum = reductionMode == MPI_SUM;
+  bool localIsMax = reductionMode == MPI_MAX;
+
   // loop over all available dimensions
   for (int i = 0; i < dim; ++i) {
     W work_local_plane;
     // collect work from all processes in the same plane
     MPI_Allreduce(this->getWork().data(), &work_local_plane, 1, MPIDataTypeW,
-                  MPI_SUM, communicators.at(i));
+                  reductionMode, communicators.at(i));
 
     // correct right border:
 
@@ -255,11 +265,28 @@ template <class T, class W> void Tensor_LB<T, W>::balance(int) {
 }
 
 // provide list of neighbors
-template <class T, class W>
-std::vector<int> &Tensor_LB<T, W>::getNeighbors() {
+template <class T, class W> std::vector<int> &Tensor_LB<T, W>::getNeighbors() {
   return neighbors;
 }
 
-}//namespace ALL
+/// Load-balancing scheme based on a independent local scheme, where the
+/// workloads of domains are reduced over their cartesian position, on reduction
+/// for each of the dimensions. These reduced values are then compared against
+/// the values gathered for the neighbor domains. Then the boundaries between
+/// the process groups are adjusted according to the ratio of cummulated work
+/// loads on each of them. The border is shifted in the direction of the domain
+/// with more cummulated work. This is independently done for each of the
+/// dimensions. For this method, the reduction is a maximum over the respective
+/// dimension, instead of a sum.
+/// @tparam T data for vertices and related data
+/// @tparam W data for work and related data
+template <class T, class W> class TensorMax_LB : public Tensor_LB<T, W> {
+public:
+  TensorMax_LB() {}
+  TensorMax_LB(int d, W w, T g) : Tensor_LB<T, W>(d, w, g) {}
+  virtual void balance(int step) override {Tensor_LB<T,W>::balance(step, MPI_MAX);}
+};
+
+} // namespace ALL
 
 #endif
diff --git a/src/ALL_module.F90 b/src/ALL_module.F90
index ad9dcfd14937b5065ece390acbc6b176960de49f..44e12bcd9373b0a2f404c6b3786d6b524b7807cc 100644
--- a/src/ALL_module.F90
+++ b/src/ALL_module.F90
@@ -35,510 +35,511 @@
 !! available there. For a simple example on how to use this see ``examples/``
 !! directory.
 module ALL_module
-    use iso_c_binding
-    use iso_fortran_env
-    implicit none
-    private
-    ! The different loadbalancing methods. Must be kept in sync with C side.
-    integer(c_int), public, parameter :: ALL_STAGGERED = 0
-    integer(c_int), public, parameter :: ALL_TENSOR = 1
-    integer(c_int), public, parameter :: ALL_FORCEBASED = 2
+   use iso_c_binding
+   use iso_fortran_env
+   implicit none
+   private
+   ! The different loadbalancing methods. Must be kept in sync with C side.
+   integer(c_int), public, parameter :: ALL_STAGGERED = 0
+   integer(c_int), public, parameter :: ALL_TENSOR = 1
+   integer(c_int), public, parameter :: ALL_FORCEBASED = 2
 #ifdef ALL_VORONOI_ACTIVE
-    integer(c_int), public, parameter :: ALL_VORONOI = 3
+   integer(c_int), public, parameter :: ALL_VORONOI = 3
 #endif
-    integer(c_int), public, parameter :: ALL_HISTOGRAM = 4
-    integer(c_int), public, parameter :: ALL_UNIMPLEMENTED = 128
-    ! The different error IDs. Must be kept in sync with CustomException class
-    integer(c_int), public, parameter :: ALL_ERROR_GENERIC = 1
-    integer(c_int), public, parameter :: ALL_ERROR_POINTDIMENSIONMISSMATCH = 2
-    integer(c_int), public, parameter :: ALL_ERROR_INVALIDCOMMTYPE = 3
-    integer(c_int), public, parameter :: ALL_ERROR_INVALIDARGUMENT = 4
-    integer(c_int), public, parameter :: ALL_ERROR_OUTOFBOUNDS = 5
-    integer(c_int), public, parameter :: ALL_ERROR_INTERNAL = 6
-    integer(c_int), public, parameter :: ALL_ERROR_FILESYSTEM = 6
-    ! Maximum character length of error descriptions
-    integer, public, parameter :: ALL_ERROR_LENGTH = 1024
-    ! Direct interface with C wrapper
-    ! TODO add intent() to dummy arguments
-    interface
-        function all_init_c(method, dim, gamma) result(this) bind(C)
-            use iso_c_binding
-            integer(c_int), value :: method
-            integer(c_int), value :: dim
-            real(c_double), value :: gamma
-            type(c_ptr) :: this
-        end function
-        subroutine all_finalize_c(obj) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-        end subroutine
-        subroutine all_set_gamma_c(obj, gamma) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            real(c_double), value :: gamma
-        end subroutine
-        subroutine all_set_proc_grid_params_c(obj, nloc, loc, nsize, size) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: nloc
-            integer(c_int), dimension(nloc) :: loc
-            integer(c_int), value :: nsize
-            integer(c_int), dimension(nsize) :: size
-        end subroutine
-        subroutine all_set_proc_tag_c(obj, tag) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: tag
-        end subroutine
-        subroutine all_set_min_domain_size_c(obj, dim, domain_size) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: dim
-            real(c_double), dimension(dim) :: domain_size
-        end subroutine
-        subroutine all_set_work_c(obj, work) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            real(c_double), value :: work
-        end subroutine
-        subroutine all_set_work_multi_c(obj, work, dim) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: dim
-            real(c_double), dimension(dim) :: work
-        end subroutine
-        subroutine all_set_vertices_c(obj, n, dim, vertices) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: n
-            integer(c_int), value :: dim
-            real(c_double), dimension(n*dim) :: vertices
-        end subroutine
-        subroutine all_set_communicator_c(obj, comm) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer, value :: comm
-        end subroutine
-        subroutine all_set_sys_size_c(obj, size, dim) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: dim
-            real(c_double), dimension(dim) :: size
-        end subroutine
-        subroutine all_set_method_data_histogram_c(obj, nbins) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), dimension(3) :: nbins
-        end subroutine
-        subroutine all_setup_c(obj) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-        end subroutine
-        subroutine all_balance_c(obj) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-        end subroutine
-        subroutine all_get_gamma_c(obj, gamma) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            real(c_double) :: gamma
-        end subroutine
-        subroutine all_get_number_of_vertices_c(obj, n) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int) :: n
-        end subroutine
-        subroutine all_get_vertices_c(obj, n, vertices) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: n
-            real(c_double), dimension(*) :: vertices
-        end subroutine
-        subroutine all_get_prev_vertices_c(obj, n, vertices) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: n
-            real(c_double), dimension(*) :: vertices
-        end subroutine
-        subroutine all_get_dimension_c(obj,dim) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int) :: dim
-        end subroutine
-        subroutine all_get_length_of_work_c(obj, length) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int) :: length
-        end subroutine
-        subroutine all_get_work_c(obj, work) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            real(c_double) :: work
-        end subroutine
-        subroutine all_get_work_array_c(obj, work, length) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: length
-            real(c_double), dimension(length) :: work
-        end subroutine
-        subroutine all_get_number_of_neighbors_c(obj, count) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int) :: count
-        end subroutine
-        subroutine all_get_neighbors_c(obj, neighbors, count) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: count
-            integer(c_int), dimension(count) :: neighbors
-        end subroutine
+   integer(c_int), public, parameter :: ALL_HISTOGRAM = 4
+   integer(c_int), public, parameter :: ALL_TENSOR_MAX = 5
+   integer(c_int), public, parameter :: ALL_UNIMPLEMENTED = 128
+   ! The different error IDs. Must be kept in sync with CustomException class
+   integer(c_int), public, parameter :: ALL_ERROR_GENERIC = 1
+   integer(c_int), public, parameter :: ALL_ERROR_POINTDIMENSIONMISSMATCH = 2
+   integer(c_int), public, parameter :: ALL_ERROR_INVALIDCOMMTYPE = 3
+   integer(c_int), public, parameter :: ALL_ERROR_INVALIDARGUMENT = 4
+   integer(c_int), public, parameter :: ALL_ERROR_OUTOFBOUNDS = 5
+   integer(c_int), public, parameter :: ALL_ERROR_INTERNAL = 6
+   integer(c_int), public, parameter :: ALL_ERROR_FILESYSTEM = 6
+   ! Maximum character length of error descriptions
+   integer, public, parameter :: ALL_ERROR_LENGTH = 1024
+   ! Direct interface with C wrapper
+   ! TODO add intent() to dummy arguments
+   interface
+      function all_init_c(method, dim, gamma) result(this) bind(C)
+         use iso_c_binding
+         integer(c_int), value :: method
+         integer(c_int), value :: dim
+         real(c_double), value :: gamma
+         type(c_ptr) :: this
+      end function
+      subroutine all_finalize_c(obj) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+      end subroutine
+      subroutine all_set_gamma_c(obj, gamma) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         real(c_double), value :: gamma
+      end subroutine
+      subroutine all_set_proc_grid_params_c(obj, nloc, loc, nsize, size) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: nloc
+         integer(c_int), dimension(nloc) :: loc
+         integer(c_int), value :: nsize
+         integer(c_int), dimension(nsize) :: size
+      end subroutine
+      subroutine all_set_proc_tag_c(obj, tag) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: tag
+      end subroutine
+      subroutine all_set_min_domain_size_c(obj, dim, domain_size) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: dim
+         real(c_double), dimension(dim) :: domain_size
+      end subroutine
+      subroutine all_set_work_c(obj, work) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         real(c_double), value :: work
+      end subroutine
+      subroutine all_set_work_multi_c(obj, work, dim) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: dim
+         real(c_double), dimension(dim) :: work
+      end subroutine
+      subroutine all_set_vertices_c(obj, n, dim, vertices) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: n
+         integer(c_int), value :: dim
+         real(c_double), dimension(n*dim) :: vertices
+      end subroutine
+      subroutine all_set_communicator_c(obj, comm) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer, value :: comm
+      end subroutine
+      subroutine all_set_sys_size_c(obj, size, dim) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: dim
+         real(c_double), dimension(dim) :: size
+      end subroutine
+      subroutine all_set_method_data_histogram_c(obj, nbins) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), dimension(3) :: nbins
+      end subroutine
+      subroutine all_setup_c(obj) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+      end subroutine
+      subroutine all_balance_c(obj) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+      end subroutine
+      subroutine all_get_gamma_c(obj, gamma) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         real(c_double) :: gamma
+      end subroutine
+      subroutine all_get_number_of_vertices_c(obj, n) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int) :: n
+      end subroutine
+      subroutine all_get_vertices_c(obj, n, vertices) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: n
+         real(c_double), dimension(*) :: vertices
+      end subroutine
+      subroutine all_get_prev_vertices_c(obj, n, vertices) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: n
+         real(c_double), dimension(*) :: vertices
+      end subroutine
+      subroutine all_get_dimension_c(obj, dim) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int) :: dim
+      end subroutine
+      subroutine all_get_length_of_work_c(obj, length) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int) :: length
+      end subroutine
+      subroutine all_get_work_c(obj, work) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         real(c_double) :: work
+      end subroutine
+      subroutine all_get_work_array_c(obj, work, length) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: length
+         real(c_double), dimension(length) :: work
+      end subroutine
+      subroutine all_get_number_of_neighbors_c(obj, count) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int) :: count
+      end subroutine
+      subroutine all_get_neighbors_c(obj, neighbors, count) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: count
+         integer(c_int), dimension(count) :: neighbors
+      end subroutine
 #ifdef ALL_VTK_OUTPUT
-        subroutine all_print_vtk_outlines_c(obj, step) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: step
-        end subroutine
-        subroutine all_print_vtk_vertices_c(obj, step) bind(C)
-            use iso_c_binding
-            type(c_ptr), value :: obj
-            integer(c_int), value :: step
-        end subroutine
+      subroutine all_print_vtk_outlines_c(obj, step) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: step
+      end subroutine
+      subroutine all_print_vtk_vertices_c(obj, step) bind(C)
+         use iso_c_binding
+         type(c_ptr), value :: obj
+         integer(c_int), value :: step
+      end subroutine
 #endif
-        function all_errno_c() result(errno) bind(C)
-            use iso_c_binding
-            integer(c_int) :: errno
-        end function
-        subroutine all_reset_errno_c() bind(C)
-        end subroutine
-        subroutine all_errdesc_c(text, length) bind(C)
-            use iso_c_binding
-            character(len=1,kind=c_char) :: text(*)
-            integer(c_size_t), value :: length
-        end subroutine
-    end interface
+      function all_errno_c() result(errno) bind(C)
+         use iso_c_binding
+         integer(c_int) :: errno
+      end function
+      subroutine all_reset_errno_c() bind(C)
+      end subroutine
+      subroutine all_errdesc_c(text, length) bind(C)
+         use iso_c_binding
+         character(len=1, kind=c_char) :: text(*)
+         integer(c_size_t), value :: length
+      end subroutine
+   end interface
 
-    !> The object oriented API is this object. It contains all relevant functions
-    type, public :: ALL_t
-        type(c_ptr), private :: object = C_NULL_PTR !< pointer to C++ object used on C side
-        integer, private :: dim = 0 !< dimensionality of system, set during init
-    contains
-        procedure :: init => ALL_init
-        procedure :: finalize => ALL_finalize
-        procedure :: set_gamma => ALL_set_gamma
-        procedure :: set_proc_grid_params => ALL_set_proc_grid_params
-        procedure :: set_proc_tag => ALL_set_proc_tag
-        procedure :: set_min_domain_size => ALL_set_min_domain_size
-        procedure :: set_work => ALL_set_work
-        procedure :: set_work_multi => ALL_set_work_multi
-        procedure :: set_vertices => ALL_set_vertices
+   !> The object oriented API is this object. It contains all relevant functions
+   type, public :: ALL_t
+      type(c_ptr), private :: object = C_NULL_PTR !< pointer to C++ object used on C side
+      integer, private :: dim = 0 !< dimensionality of system, set during init
+   contains
+      procedure :: init => ALL_init
+      procedure :: finalize => ALL_finalize
+      procedure :: set_gamma => ALL_set_gamma
+      procedure :: set_proc_grid_params => ALL_set_proc_grid_params
+      procedure :: set_proc_tag => ALL_set_proc_tag
+      procedure :: set_min_domain_size => ALL_set_min_domain_size
+      procedure :: set_work => ALL_set_work
+      procedure :: set_work_multi => ALL_set_work_multi
+      procedure :: set_vertices => ALL_set_vertices
 #ifdef ALL_USE_F08
-        procedure :: set_communicator => ALL_set_communicator_f08
+      procedure :: set_communicator => ALL_set_communicator_f08
 #else
-        procedure :: set_communicator => ALL_set_communicator_int
+      procedure :: set_communicator => ALL_set_communicator_int
 #endif
-        procedure :: set_sys_size => ALL_set_sys_size
-        procedure :: set_method_data_histgram => ALL_set_method_data_histgram
-        procedure :: setup => ALL_setup
-        procedure :: balance => ALL_balance
-        procedure :: get_gamma => ALL_get_gamma
-        procedure :: get_number_of_vertices => ALL_get_number_of_vertices
-        procedure :: get_vertices => ALL_get_vertices
-        procedure :: get_vertices_alloc => ALL_get_vertices_alloc
-        procedure :: get_prev_vertices => ALL_get_prev_vertices
-        procedure :: get_dimension => ALL_get_dimension
-        procedure :: get_length_of_work => ALL_get_length_of_work
-        procedure :: get_work => ALL_get_work
-        procedure :: get_work_array => ALL_get_work_array
-        procedure :: get_number_of_neighbors => ALL_get_number_of_neighbors
-        procedure :: get_neighbors => ALL_get_neighbors
+      procedure :: set_sys_size => ALL_set_sys_size
+      procedure :: set_method_data_histgram => ALL_set_method_data_histgram
+      procedure :: setup => ALL_setup
+      procedure :: balance => ALL_balance
+      procedure :: get_gamma => ALL_get_gamma
+      procedure :: get_number_of_vertices => ALL_get_number_of_vertices
+      procedure :: get_vertices => ALL_get_vertices
+      procedure :: get_vertices_alloc => ALL_get_vertices_alloc
+      procedure :: get_prev_vertices => ALL_get_prev_vertices
+      procedure :: get_dimension => ALL_get_dimension
+      procedure :: get_length_of_work => ALL_get_length_of_work
+      procedure :: get_work => ALL_get_work
+      procedure :: get_work_array => ALL_get_work_array
+      procedure :: get_number_of_neighbors => ALL_get_number_of_neighbors
+      procedure :: get_neighbors => ALL_get_neighbors
 #ifdef ALL_VTK_OUTPUT
-        procedure :: print_vtk_outlines => ALL_print_vtk_outlines
-        procedure :: print_vtk_vertices => ALL_print_vtk_vertices
+      procedure :: print_vtk_outlines => ALL_print_vtk_outlines
+      procedure :: print_vtk_vertices => ALL_print_vtk_vertices
 #endif
-    end type
+   end type
 
-    ! This is the old, but still supported API of separate functions
-    public :: ALL_init
-    public :: ALL_finalize
-    public :: ALL_set_gamma
-    public :: ALL_set_proc_grid_params
-    public :: ALL_set_proc_tag
-    public :: ALL_set_min_domain_size
-    public :: ALL_set_work
-    public :: ALL_set_work_multi
-    public :: ALL_set_vertices
-    public :: ALL_set_communicator
-    public :: ALL_set_sys_size
-    public :: ALL_set_method_data_histgram
-    public :: ALL_setup
-    public :: ALL_balance
-    public :: ALL_get_gamma
-    public :: ALL_get_number_of_vertices
-    public :: ALL_get_vertices
-    public :: ALL_get_vertices_alloc
-    public :: ALL_get_prev_vertices
-    public :: ALL_get_dimension
-    public :: ALL_get_length_of_work
-    public :: ALL_get_work
-    public :: ALL_get_work_array
-    public :: ALL_get_number_of_neighbors
-    public :: ALL_get_neighbors
+   ! This is the old, but still supported API of separate functions
+   public :: ALL_init
+   public :: ALL_finalize
+   public :: ALL_set_gamma
+   public :: ALL_set_proc_grid_params
+   public :: ALL_set_proc_tag
+   public :: ALL_set_min_domain_size
+   public :: ALL_set_work
+   public :: ALL_set_work_multi
+   public :: ALL_set_vertices
+   public :: ALL_set_communicator
+   public :: ALL_set_sys_size
+   public :: ALL_set_method_data_histgram
+   public :: ALL_setup
+   public :: ALL_balance
+   public :: ALL_get_gamma
+   public :: ALL_get_number_of_vertices
+   public :: ALL_get_vertices
+   public :: ALL_get_vertices_alloc
+   public :: ALL_get_prev_vertices
+   public :: ALL_get_dimension
+   public :: ALL_get_length_of_work
+   public :: ALL_get_work
+   public :: ALL_get_work_array
+   public :: ALL_get_number_of_neighbors
+   public :: ALL_get_neighbors
 #ifdef ALL_VTK_OUTPUT
-    public :: ALL_print_vtk_outlines
-    public :: ALL_print_vtk_vertices
+   public :: ALL_print_vtk_outlines
+   public :: ALL_print_vtk_vertices
 #endif
-    public :: ALL_error
-    public :: ALL_reset_error
-    public :: ALL_error_description
+   public :: ALL_error
+   public :: ALL_reset_error
+   public :: ALL_error_description
 
-    interface ALL_set_communicator
-        module procedure ALL_set_communicator_int
+   interface ALL_set_communicator
+      module procedure ALL_set_communicator_int
 #ifdef ALL_USE_F08
-        module procedure ALL_set_communicator_f08
+      module procedure ALL_set_communicator_f08
 #endif
-    end interface
+   end interface
 contains
-    !> Initialises the loadbalancer
-    subroutine ALL_init(this, method, dim, gamma)
-        class(ALL_t), intent(out) :: this !< teh ALL object is returned
-        integer, intent(in) :: method !< Must be one of the ALL_* method values
-        integer, intent(in) :: dim !< dimensionality of system
-        real(c_double), intent(in) :: gamma !< gamma value for load balancer (ignored for TENSOR and STAGGERED)
-        this%object = all_init_c(method, dim, gamma)
-        this%dim = dim
-    end subroutine
-    !> Delete the loadbalance object
-    subroutine ALL_finalize(this)
-        class(ALL_t), intent(in) :: this
-        call all_finalize_c(this%object)
-    end subroutine
-    !> Set gamma value for load balancer
-    subroutine ALL_set_gamma(this, gamma)
-        class(ALL_t), intent(in) :: this
-        real(c_double), intent(in) :: gamma
-        call all_set_gamma_c(this%object, gamma)
-    end subroutine
-    !> Set the grid parameters for this process
-    subroutine ALL_set_proc_grid_params(this, loc, ranks)
-        class(ALL_t), intent(in) :: this
-        integer, dimension(this%dim), intent(in) :: loc !< index of domain in `dim` directions (0-indexed!)
-        integer, dimension(this%dim), intent(in) :: ranks !< total number of domains in `dim` directions
-        call all_set_proc_grid_params_c(this%object, this%dim, loc, this%dim, ranks)
-    end subroutine
-    !> Set process identifying tag for output
-    subroutine ALL_set_proc_tag(this, tag)
-        class(ALL_t), intent(in) :: this
-        integer, intent(in) :: tag !< tag of local process, only output in VTK outlines
-        call all_set_proc_tag_c(this%object, tag)
-    end subroutine
-    !> Set a minimum domain size
-    subroutine ALL_set_min_domain_size(this, domain_size)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(this%dim), intent(in) :: domain_size !< minimum domain size
-        call all_set_min_domain_size_c(this%object, this%dim, domain_size)
-    end subroutine
-    !> Set the work of this process
-    subroutine ALL_set_work(this, work)
-        class(ALL_t), intent(in) :: this
-        real(c_double), intent(in) :: work !< work of this domain
-        call all_set_work_c(this%object, work)
-    end subroutine
-    !> Set multi dimensional work of this process
-    subroutine ALL_set_work_multi(this, work)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:), intent(in) :: work !< multi dimensional work of this domain
+   !> Initialises the loadbalancer
+   subroutine ALL_init(this, method, dim, gamma)
+      class(ALL_t), intent(out) :: this !< teh ALL object is returned
+      integer, intent(in) :: method !< Must be one of the ALL_* method values
+      integer, intent(in) :: dim !< dimensionality of system
+      real(c_double), intent(in) :: gamma !< gamma value for load balancer (ignored for TENSOR and STAGGERED)
+      this%object = all_init_c(method, dim, gamma)
+      this%dim = dim
+   end subroutine
+   !> Delete the loadbalance object
+   subroutine ALL_finalize(this)
+      class(ALL_t), intent(in) :: this
+      call all_finalize_c(this%object)
+   end subroutine
+   !> Set gamma value for load balancer
+   subroutine ALL_set_gamma(this, gamma)
+      class(ALL_t), intent(in) :: this
+      real(c_double), intent(in) :: gamma
+      call all_set_gamma_c(this%object, gamma)
+   end subroutine
+   !> Set the grid parameters for this process
+   subroutine ALL_set_proc_grid_params(this, loc, ranks)
+      class(ALL_t), intent(in) :: this
+      integer, dimension(this%dim), intent(in) :: loc !< index of domain in `dim` directions (0-indexed!)
+      integer, dimension(this%dim), intent(in) :: ranks !< total number of domains in `dim` directions
+      call all_set_proc_grid_params_c(this%object, this%dim, loc, this%dim, ranks)
+   end subroutine
+   !> Set process identifying tag for output
+   subroutine ALL_set_proc_tag(this, tag)
+      class(ALL_t), intent(in) :: this
+      integer, intent(in) :: tag !< tag of local process, only output in VTK outlines
+      call all_set_proc_tag_c(this%object, tag)
+   end subroutine
+   !> Set a minimum domain size
+   subroutine ALL_set_min_domain_size(this, domain_size)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(this%dim), intent(in) :: domain_size !< minimum domain size
+      call all_set_min_domain_size_c(this%object, this%dim, domain_size)
+   end subroutine
+   !> Set the work of this process
+   subroutine ALL_set_work(this, work)
+      class(ALL_t), intent(in) :: this
+      real(c_double), intent(in) :: work !< work of this domain
+      call all_set_work_c(this%object, work)
+   end subroutine
+   !> Set multi dimensional work of this process
+   subroutine ALL_set_work_multi(this, work)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:), intent(in) :: work !< multi dimensional work of this domain
 #ifndef NDEBUG
-        if(size(work) /= this%dim) then
-            write(error_unit, '(a)') "ALL_set_work_multi: Wrong dimension for work"
-            stop
-        endif
+      if (size(work) /= this%dim) then
+         write (error_unit, '(a)') "ALL_set_work_multi: Wrong dimension for work"
+         stop
+      end if
 #endif
-        call all_set_work_multi_c(this%object, work, size(work))
-    end subroutine
-    !> Set new vertices
+      call all_set_work_multi_c(this%object, work, size(work))
+   end subroutine
+   !> Set new vertices
     !!
     !! @todo write interface for single rank array of vertices
-    subroutine ALL_set_vertices(this, vertices)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:,:), intent(in) :: vertices !< vertices of domain, for `n` domains must have shape: (dim,n)
+   subroutine ALL_set_vertices(this, vertices)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:, :), intent(in) :: vertices !< vertices of domain, for `n` domains must have shape: (dim,n)
 #ifndef NDEBUG
-        if(size(vertices,1) /= this%dim) then
-            write(error_unit,'(a)') "ALL_set_vertices: Wrong dimension for vertices"
-            stop
-        endif
+      if (size(vertices, 1) /= this%dim) then
+         write (error_unit, '(a)') "ALL_set_vertices: Wrong dimension for vertices"
+         stop
+      end if
 #endif
-        call all_set_vertices_c(this%object, size(vertices,2), this%dim, vertices)
-    end subroutine
-    !> Set the MPI communicator with an mpi_f08 object
+      call all_set_vertices_c(this%object, size(vertices, 2), this%dim, vertices)
+   end subroutine
+   !> Set the MPI communicator with an mpi_f08 object
 #ifdef ALL_USE_F08
-    subroutine ALL_set_communicator_f08(this, comm)
-        use mpi_f08
-        class(ALL_t), intent(in) :: this
-        type(MPI_Comm), intent(in) :: comm !< MPI Communicator
-        call all_set_communicator_c(this%object, comm%mpi_val)
-    end subroutine
+   subroutine ALL_set_communicator_f08(this, comm)
+      use mpi_f08
+      class(ALL_t), intent(in) :: this
+      type(MPI_Comm), intent(in) :: comm !< MPI Communicator
+      call all_set_communicator_c(this%object, comm%mpi_val)
+   end subroutine
 #endif
-    !> Set the MPI communicator with an ``mpi`` oder ``mpif.h`` communicator
-    subroutine ALL_set_communicator_int(this, comm)
-        class(ALL_t), intent(in) :: this
-        integer, intent(in) :: comm !< MPI Communicator, not type checked!
-        call all_set_communicator_c(this%object, comm)
-    end subroutine
-    !> Set size of system, which is required for the histogram method
-    subroutine ALL_set_sys_size(this, syssize)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:), intent(in) :: syssize !< system size
+   !> Set the MPI communicator with an ``mpi`` oder ``mpif.h`` communicator
+   subroutine ALL_set_communicator_int(this, comm)
+      class(ALL_t), intent(in) :: this
+      integer, intent(in) :: comm !< MPI Communicator, not type checked!
+      call all_set_communicator_c(this%object, comm)
+   end subroutine
+   !> Set size of system, which is required for the histogram method
+   subroutine ALL_set_sys_size(this, syssize)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:), intent(in) :: syssize !< system size
 #ifndef NDEBUG
-        if(size(syssize) /= this%dim) then
-            write(error_unit, '(a)') "ALL_set_size: Wrong dimension for size"
-            stop
-        endif
+      if (size(syssize) /= this%dim) then
+         write (error_unit, '(a)') "ALL_set_size: Wrong dimension for size"
+         stop
+      end if
 #endif
-        call all_set_sys_size_c(this%object, syssize, size(syssize))
-    end subroutine
-    !> Set number of bins for histogram method
-    subroutine ALL_set_method_data_histgram(this, nbins)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), dimension(3), intent(in) :: nbins !< Number of bins per dimension
-        call all_set_method_data_histogram_c(this%object, nbins)
-    end subroutine
-    !> Set up the loadbalancer
-    subroutine ALL_setup(this)
-        class(ALL_t), intent(in) :: this
-        call all_setup_c(this%object)
-    end subroutine
-    !> Run loadbalancer and calculate new vertices
-    subroutine ALL_balance(this)
-        class(ALL_t), intent(in) :: this
-        call all_balance_c(this%object)
-    end subroutine
-    !> Retrieve currently set gamma value of balancer
-    subroutine ALL_get_gamma(this, gamma)
-        class(ALL_t), intent(in) :: this
-        real(c_double), intent(out) :: gamma
-        call all_get_gamma_c(this%object, gamma)
-    end subroutine
-    !> Retrieve number of vertices for current vertices
-    subroutine ALL_get_number_of_vertices(this, n)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), intent(out) :: n !< set to number of new vertices
-        call all_get_number_of_vertices_c(this%object, n)
-    end subroutine
-    !> Retrieve current vertices
-    subroutine ALL_get_vertices(this, vertices)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:,:), intent(out) :: vertices !< set to new vertices, must be exact size (dim,n)
-        integer(c_int) :: n_verts
-        call this%get_number_of_vertices(n_verts)
+      call all_set_sys_size_c(this%object, syssize, size(syssize))
+   end subroutine
+   !> Set number of bins for histogram method
+   subroutine ALL_set_method_data_histgram(this, nbins)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), dimension(3), intent(in) :: nbins !< Number of bins per dimension
+      call all_set_method_data_histogram_c(this%object, nbins)
+   end subroutine
+   !> Set up the loadbalancer
+   subroutine ALL_setup(this)
+      class(ALL_t), intent(in) :: this
+      call all_setup_c(this%object)
+   end subroutine
+   !> Run loadbalancer and calculate new vertices
+   subroutine ALL_balance(this)
+      class(ALL_t), intent(in) :: this
+      call all_balance_c(this%object)
+   end subroutine
+   !> Retrieve currently set gamma value of balancer
+   subroutine ALL_get_gamma(this, gamma)
+      class(ALL_t), intent(in) :: this
+      real(c_double), intent(out) :: gamma
+      call all_get_gamma_c(this%object, gamma)
+   end subroutine
+   !> Retrieve number of vertices for current vertices
+   subroutine ALL_get_number_of_vertices(this, n)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), intent(out) :: n !< set to number of new vertices
+      call all_get_number_of_vertices_c(this%object, n)
+   end subroutine
+   !> Retrieve current vertices
+   subroutine ALL_get_vertices(this, vertices)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:, :), intent(out) :: vertices !< set to new vertices, must be exact size (dim,n)
+      integer(c_int) :: n_verts
+      call this%get_number_of_vertices(n_verts)
 #ifndef NDEBUG
-        if(size(vertices,1) /= this%dim .or. size(vertices,2) < n_verts) then
-            write(error_unit,'(a)') "ALL_get_vertices: vertices array not large enough!"
-            stop
-        endif
+      if (size(vertices, 1) /= this%dim .or. size(vertices, 2) < n_verts) then
+         write (error_unit, '(a)') "ALL_get_vertices: vertices array not large enough!"
+         stop
+      end if
 #endif
-        call all_get_vertices_c(this%object, n_verts, vertices)
-    end subroutine
-    !> Same function as get_vertices, but takes an allocatable array, and resizes automatically
-    subroutine ALL_get_vertices_alloc(this, vertices)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:,:), allocatable, intent(inout) :: vertices !< set to new vertices, may be reallocated to fit
-        integer(c_int) :: n_verts
-        call this%get_number_of_vertices(n_verts)
-        if(allocated(vertices)) then
-            if(size(vertices,1) /= this%dim .or. size(vertices,2) < n_verts) then
-                deallocate(vertices)
-            endif
-        endif
-        if(.not.allocated(vertices)) allocate(vertices(this%dim, n_verts))
-        call all_get_vertices_c(this%object, n_verts, vertices)
-    end subroutine
-    !> Retrieve vertices from before loadbalancing
-    subroutine ALL_get_prev_vertices(this, vertices)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:,:), intent(out) :: vertices !< set to prev vertices, must be exact size (dim,n), unchecked!
+      call all_get_vertices_c(this%object, n_verts, vertices)
+   end subroutine
+   !> Same function as get_vertices, but takes an allocatable array, and resizes automatically
+   subroutine ALL_get_vertices_alloc(this, vertices)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:, :), allocatable, intent(inout) :: vertices !< set to new vertices, may be reallocated to fit
+      integer(c_int) :: n_verts
+      call this%get_number_of_vertices(n_verts)
+      if (allocated(vertices)) then
+         if (size(vertices, 1) /= this%dim .or. size(vertices, 2) < n_verts) then
+            deallocate (vertices)
+         end if
+      end if
+      if (.not. allocated(vertices)) allocate (vertices(this%dim, n_verts))
+      call all_get_vertices_c(this%object, n_verts, vertices)
+   end subroutine
+   !> Retrieve vertices from before loadbalancing
+   subroutine ALL_get_prev_vertices(this, vertices)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:, :), intent(out) :: vertices !< set to prev vertices, must be exact size (dim,n), unchecked!
 #ifndef NDEBUG
-        ! TODO Check against number of vertices not only dimensionality
-        if(size(vertices,1) /= this%dim) then
-            write(error_unit,'(a)') "ALL_get_prev_vertices: vertices array not large enough!"
-            stop
-        endif
+      ! TODO Check against number of vertices not only dimensionality
+      if (size(vertices, 1) /= this%dim) then
+         write (error_unit, '(a)') "ALL_get_prev_vertices: vertices array not large enough!"
+         stop
+      end if
 #endif
-        call all_get_prev_vertices_c(this%object, size(vertices,2), vertices)
-    end subroutine
-    !> Get current dimension from loadbalancer
-    subroutine ALL_get_dimension(this, dim)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), intent(out) :: dim
-        call all_get_dimension_c(this%object,dim)
-    end subroutine
-    !> Retrieve length of work array
-    subroutine ALL_get_length_of_work(this, length)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), intent(out) :: length
-        call all_get_length_of_work_c(this%object, length)
-    end subroutine
-    !> Retrieve first element of work array
-    subroutine ALL_get_work(this, work)
-        class(ALL_t), intent(in) :: this
-        real(c_double), intent(out) :: work
-        call all_get_work_c(this%object, work)
-    end subroutine
-    !> Retrieve work array, which must already be the correct size
-    subroutine ALL_get_work_array(this, work)
-        class(ALL_t), intent(in) :: this
-        real(c_double), dimension(:), intent(out) :: work
-        integer :: length
-        call ALL_get_length_of_work(this, length)
-        if(size(work) /= length) then
-            write(error_unit,'(a)') "ALL_get_work_array: work has wrong length!"
-            stop
-        endif
-        call all_get_work_array_c(this%object, work, size(work))
-    end subroutine
-    !> Retrieve number of neighbors (i.e. length of neighbors list)
-    subroutine ALL_get_number_of_neighbors(this, count)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), intent(out) :: count
-        call all_get_number_of_vertices_c(this%object, count)
-    end subroutine
-    !> Retrieve list of neighboring ranks (must be correct size already)
-    subroutine ALL_get_neighbors(this, neighbors)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), dimension(:), intent(out) :: neighbors
-        integer :: count
-        call ALL_get_number_of_neighbors(this, count)
-        if(size(neighbors) /= count) then
-            write(error_unit,'(a)') "ALL_get_neighbors: neighbors has wrong length!"
-            stop
-        endif
-        call all_get_neighbors_c(this%object, neighbors, size(neighbors))
-    end subroutine
+      call all_get_prev_vertices_c(this%object, size(vertices, 2), vertices)
+   end subroutine
+   !> Get current dimension from loadbalancer
+   subroutine ALL_get_dimension(this, dim)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), intent(out) :: dim
+      call all_get_dimension_c(this%object, dim)
+   end subroutine
+   !> Retrieve length of work array
+   subroutine ALL_get_length_of_work(this, length)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), intent(out) :: length
+      call all_get_length_of_work_c(this%object, length)
+   end subroutine
+   !> Retrieve first element of work array
+   subroutine ALL_get_work(this, work)
+      class(ALL_t), intent(in) :: this
+      real(c_double), intent(out) :: work
+      call all_get_work_c(this%object, work)
+   end subroutine
+   !> Retrieve work array, which must already be the correct size
+   subroutine ALL_get_work_array(this, work)
+      class(ALL_t), intent(in) :: this
+      real(c_double), dimension(:), intent(out) :: work
+      integer :: length
+      call ALL_get_length_of_work(this, length)
+      if (size(work) /= length) then
+         write (error_unit, '(a)') "ALL_get_work_array: work has wrong length!"
+         stop
+      end if
+      call all_get_work_array_c(this%object, work, size(work))
+   end subroutine
+   !> Retrieve number of neighbors (i.e. length of neighbors list)
+   subroutine ALL_get_number_of_neighbors(this, count)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), intent(out) :: count
+      call all_get_number_of_vertices_c(this%object, count)
+   end subroutine
+   !> Retrieve list of neighboring ranks (must be correct size already)
+   subroutine ALL_get_neighbors(this, neighbors)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), dimension(:), intent(out) :: neighbors
+      integer :: count
+      call ALL_get_number_of_neighbors(this, count)
+      if (size(neighbors) /= count) then
+         write (error_unit, '(a)') "ALL_get_neighbors: neighbors has wrong length!"
+         stop
+      end if
+      call all_get_neighbors_c(this%object, neighbors, size(neighbors))
+   end subroutine
 
 #ifdef ALL_VTK_OUTPUT
-    !> Print VTK outlines (must be enabled in build step)
-    subroutine ALL_print_vtk_outlines(this, step)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), intent(in) :: step !< current step, used for filename
-        call all_print_vtk_outlines_c(this%object, step)
-    end subroutine
-    !> Print VTK domain vertices (must be enabled in build step)
-    subroutine ALL_print_vtk_vertices(this, step)
-        class(ALL_t), intent(in) :: this
-        integer(c_int), intent(in) :: step !< current step, used for filename
-        call all_print_vtk_vertices_c(this%object, step)
-    end subroutine
+   !> Print VTK outlines (must be enabled in build step)
+   subroutine ALL_print_vtk_outlines(this, step)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), intent(in) :: step !< current step, used for filename
+      call all_print_vtk_outlines_c(this%object, step)
+   end subroutine
+   !> Print VTK domain vertices (must be enabled in build step)
+   subroutine ALL_print_vtk_vertices(this, step)
+      class(ALL_t), intent(in) :: this
+      integer(c_int), intent(in) :: step !< current step, used for filename
+      call all_print_vtk_vertices_c(this%object, step)
+   end subroutine
 #endif
-    function ALL_error() result(error)
-        integer :: error
-        error = all_errno_c()
-    end function
-    subroutine ALL_reset_error()
-        call all_reset_errno_c()
-    end subroutine
-    function ALL_error_description() result(desc)
-        character(len=ALL_ERROR_LENGTH) :: desc
-        call all_errdesc_c(desc, int(ALL_ERROR_LENGTH,c_size_t))
-    end function
+   function ALL_error() result(error)
+      integer :: error
+      error = all_errno_c()
+   end function
+   subroutine ALL_reset_error()
+      call all_reset_errno_c()
+   end subroutine
+   function ALL_error_description() result(desc)
+      character(len=ALL_ERROR_LENGTH) :: desc
+      call all_errdesc_c(desc, int(ALL_ERROR_LENGTH, c_size_t))
+   end function
 end module
 
 ! VIM modeline for indentation