diff --git a/README_new.md b/README_new.md
index 92afa7fd847394fbe19e4ba4119ad9a89d0b1c52..74686b26721e16bd69f65c9710859b15d2e1de9c 100644
--- a/README_new.md
+++ b/README_new.md
@@ -123,6 +123,10 @@ from `ALL.hpp`. All classes are encapsulated in the `ALL` name space.
 Simple (and not so simple) examples are available in the `/examples` sub
 directory. The simple examples are also documented at length.
 
+Errors are treated as exceptions that are thrown of a (sub) class of
+`ALL::CustomException`. Likely candidates that may throw are `balance`,
+`setup`, `printVTKoutlines`.
+
 ### ALL object
 The ALL object can be constructed with
 
@@ -239,6 +243,12 @@ becomes
     call ALL_get_work(work) !or
     !call all%get_work(work)
 
+Since there is no exception handling in Fortran, the error handling from
+the libc is borrowed. So there are additional function that retrieve and
+reset an error number as well as an additional error message describing
+the error. These must be queried and handled by the user. An example of
+this usage is shown in the `ALL_Staggered_f` example when printing the
+VTK outlines.
 
 ### Details of the balancing methods
 #### Tensor based
diff --git a/include/ALL.hpp b/include/ALL.hpp
index 9635dfff53840c9266405a8027a4857e92269ed0..903b5b9a8905c7996d2a361683afabef0ae0ef8d 100644
--- a/include/ALL.hpp
+++ b/include/ALL.hpp
@@ -283,12 +283,7 @@ public:
   /// method to call the setup of the chosen balancing method (not all methods
   /// require a setup)
   void setup() {
-    try {
-      balancer->setup();
-    } catch (CustomException &e) {
-      std::cout << e.what() << std::endl;
-      MPI_Abort(MPI_COMM_WORLD, -1);
-    }
+    balancer->setup();
   }
 
   /// method the trigger the balancing step, that updates the vertices according
@@ -297,107 +292,102 @@ public:
   /// for some recursive calls of the routine by some methods
   /// @result std::vector<ALL::Point<T>> containing the shifted set of vertices
   std::vector<Point<T>> &balance() {
-    try {
-      nVertices = balancer->getNVertices();
-      switch (method) {
-      case LB_t::TENSOR:
-        balancer->balance(loadbalancing_step);
-        break;
-      case LB_t::STAGGERED:
-
-        /* ToDo: Check if repetition is useful at all and
-         *       change it to be included for all methods,
-         *       not only STAGGERED */
-        /*
-        // estimation to improve speed of convergence
-        // changing vertices during estimation
-        std::vector<Point<T>> tmp_vertices = balancer->getVertices();
-        // saved original vertices
-        std::vector<Point<T>> old_vertices = balancer->getVertices();
-        // old temporary vertices
-        std::vector<Point<T>> old_tmp_vertices = balancer->getVertices();
-        // result temporary vertices
-        std::vector<Point<T>> result_vertices = balancer->getVertices();
-        // temporary work
-        W tmp_work = balancer->getWork().at(0);
-        // old temporary work
-        W old_tmp_work = balancer->getWork().at(0);
-
-        W max_work;
-        W sum_work;
-        double d_max;
-        int n_ranks;
-
-        // compute work density on local domain
-        double V = 1.0;
-        for (int i = 0; i < balancer->getDimension(); ++i)
-            V *= ( outline.at(1)[i] - outline.at(0)[i] );
-        double rho = workArray.at(0) / V;
-
-        // collect maximum work in system
-        MPI_Allreduce(workArray.data(), &max_work, 1, MPIDataTypeW, MPI_MAX,
-        communicator); MPI_Allreduce(workArray.data(), &sum_work, 1,
-        MPIDataTypeW, MPI_SUM, communicator);
-        MPI_Comm_size(communicator,&n_ranks);
-        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
-        for (int i = 0; i < estimation_limit && d_max > 0.1 && !internal; ++i)
-        {
-            balancer->setWork(tmp_work);
-            balancer->setup();
-            balancer->balance(true);
-            bool sane = true;
-            // check if the estimated boundaries are not too deformed
-            for (int j = 0; j < balancer->getDimension(); ++j)
-                sane = sane && (result_vertices.at(0)[j] <=
-        old_vertices.at(1)[j]);
-            MPI_Allreduce(MPI_IN_PLACE,&sane,1,MPI_CXX_BOOL,MPI_LAND,communicator);
-            if (sane)
-            {
-                old_tmp_vertices = tmp_vertices;
-                tmp_vertices = result_vertices;
-                V = 1.0;
-                for (int i = 0; i < balancer->getDimension(); ++i)
-                    V *= ( tmp_vertices.at(1)[i] - tmp_vertices.at(0)[i] );
-                old_tmp_work = tmp_work;
-                tmp_work = rho * V;
-            }
-            else if (!sane || i == estimation_limit - 1)
-            {
-                balancer->getVertices() = old_tmp_vertices;
-                workArray->at(0) = old_tmp_work;
-                calculate_outline();
-                i = estimation_limit;
-            }
-        }
+    nVertices = balancer->getNVertices();
+    switch (method) {
+    case LB_t::TENSOR:
+      balancer->balance(loadbalancing_step);
+      break;
+    case LB_t::STAGGERED:
+
+      /* ToDo: Check if repetition is useful at all and
+       *       change it to be included for all methods,
+       *       not only STAGGERED */
+      /*
+      // estimation to improve speed of convergence
+      // changing vertices during estimation
+      std::vector<Point<T>> tmp_vertices = balancer->getVertices();
+      // saved original vertices
+      std::vector<Point<T>> old_vertices = balancer->getVertices();
+      // old temporary vertices
+      std::vector<Point<T>> old_tmp_vertices = balancer->getVertices();
+      // result temporary vertices
+      std::vector<Point<T>> result_vertices = balancer->getVertices();
+      // temporary work
+      W tmp_work = balancer->getWork().at(0);
+      // old temporary work
+      W old_tmp_work = balancer->getWork().at(0);
+
+      W max_work;
+      W sum_work;
+      double d_max;
+      int n_ranks;
+
+      // compute work density on local domain
+      double V = 1.0;
+      for (int i = 0; i < balancer->getDimension(); ++i)
+          V *= ( outline.at(1)[i] - outline.at(0)[i] );
+      double rho = workArray.at(0) / V;
+
+      // collect maximum work in system
+      MPI_Allreduce(workArray.data(), &max_work, 1, MPIDataTypeW, MPI_MAX,
+      communicator); MPI_Allreduce(workArray.data(), &sum_work, 1,
+      MPIDataTypeW, MPI_SUM, communicator);
+      MPI_Comm_size(communicator,&n_ranks);
+      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
+      for (int i = 0; i < estimation_limit && d_max > 0.1 && !internal; ++i)
+      {
+          balancer->setWork(tmp_work);
+          balancer->setup();
+          balancer->balance(true);
+          bool sane = true;
+          // check if the estimated boundaries are not too deformed
+          for (int j = 0; j < balancer->getDimension(); ++j)
+              sane = sane && (result_vertices.at(0)[j] <=
+      old_vertices.at(1)[j]);
+          MPI_Allreduce(MPI_IN_PLACE,&sane,1,MPI_CXX_BOOL,MPI_LAND,communicator);
+          if (sane)
+          {
+              old_tmp_vertices = tmp_vertices;
+              tmp_vertices = result_vertices;
+              V = 1.0;
+              for (int i = 0; i < balancer->getDimension(); ++i)
+                  V *= ( tmp_vertices.at(1)[i] - tmp_vertices.at(0)[i] );
+              old_tmp_work = tmp_work;
+              tmp_work = rho * V;
+          }
+          else if (!sane || i == estimation_limit - 1)
+          {
+              balancer->getVertices() = old_tmp_vertices;
+              workArray->at(0) = old_tmp_work;
+              calculate_outline();
+              i = estimation_limit;
+          }
+      }
 
-        ((Staggered_LB<T,W>*)balancer.get())->setVertices(outline->data());
-        */
-        balancer->balance(loadbalancing_step);
-        break;
-      case LB_t::UNSTRUCTURED:
-        balancer->balance(loadbalancing_step);
-        break;
+      ((Staggered_LB<T,W>*)balancer.get())->setVertices(outline->data());
+      */
+      balancer->balance(loadbalancing_step);
+      break;
+    case LB_t::UNSTRUCTURED:
+      balancer->balance(loadbalancing_step);
+      break;
 #ifdef ALL_VORONOI_ACTIVE
-      case LB_t::VORONOI:
-        balancer->balance(loadbalancing_step);
-        break;
+    case LB_t::VORONOI:
+      balancer->balance(loadbalancing_step);
+      break;
 #endif
-      case LB_t::HISTOGRAM:
-        balancer->balance(loadbalancing_step);
-        break;
-      default:
-        throw InvalidArgumentException(
-            __FILE__, __func__, __LINE__,
-            "Unknown type of loadbalancing passed.");
-      }
-      loadbalancing_step++;
-    } catch (CustomException &e) {
-      std::cout << e.what() << std::endl;
-      MPI_Abort(MPI_COMM_WORLD, -1);
+    case LB_t::HISTOGRAM:
+      balancer->balance(loadbalancing_step);
+      break;
+    default:
+      throw InvalidArgumentException(
+          __FILE__, __func__, __LINE__,
+          "Unknown type of loadbalancing passed.");
     }
+    loadbalancing_step++;
     return balancer->getVertices();
   }
 
@@ -499,12 +489,7 @@ public:
   /// next-near neighbors instead of only with next neighbors
   /// @param[in] minSize the minimum size of a domain in each dimension
   void setMinDomainSize(const std::vector<T> &minSize) {
-    try {
-      (balancer.get())->setMinDomainSize(minSize);
-    } catch (CustomException &e) {
-      std::cout << e.what() << std::endl;
-      MPI_Abort(MPI_COMM_WORLD, -1);
-    }
+    (balancer.get())->setMinDomainSize(minSize);
   }
 
 #ifdef ALL_VTK_OUTPUT
@@ -666,132 +651,126 @@ public:
   /// @param[in] step the number of the loadbalancing step used for numbering
   /// the output files
   void printVTKvertices(const int step) {
-    try {
-      int n_ranks;
-      int local_rank;
+    int n_ranks;
+    int local_rank;
 
-      MPI_Comm_rank(communicator, &local_rank);
-      MPI_Comm_size(communicator, &n_ranks);
+    MPI_Comm_rank(communicator, &local_rank);
+    MPI_Comm_size(communicator, &n_ranks);
 
-      // local vertices
-      // (vertices + work)
-      T local_vertices[nVertices * balancer->getDimension() + 1];
+    // local vertices
+    // (vertices + work)
+    T local_vertices[nVertices * balancer->getDimension() + 1];
 
-      for (int v = 0; v < nVertices; ++v) {
-        for (int d = 0; d < balancer->getDimension(); ++d) {
-          local_vertices[v * balancer->getDimension() + d] =
-              balancer->getVertices().at(v)[d];
-        }
+    for (int v = 0; v < nVertices; ++v) {
+      for (int d = 0; d < balancer->getDimension(); ++d) {
+        local_vertices[v * balancer->getDimension() + d] =
+            balancer->getVertices().at(v)[d];
       }
-      local_vertices[nVertices * balancer->getDimension()] =
-          (T)balancer->getWork().at(0);
+    }
+    local_vertices[nVertices * balancer->getDimension()] =
+        (T)balancer->getWork().at(0);
 
-      T *global_vertices;
-      if (local_rank == 0) {
-        global_vertices =
-            new T[n_ranks * (nVertices * balancer->getDimension() + 1)];
-      }
+    T *global_vertices;
+    if (local_rank == 0) {
+      global_vertices =
+          new T[n_ranks * (nVertices * balancer->getDimension() + 1)];
+    }
 
-      // collect all works and vertices on a single process
-      MPI_Gather(local_vertices, nVertices * balancer->getDimension() + 1,
-                 MPIDataTypeT, global_vertices,
-                 nVertices * balancer->getDimension() + 1, MPIDataTypeT, 0,
-                 communicator);
-
-      if (local_rank == 0) {
-        auto vtkpoints = vtkSmartPointer<vtkPoints>::New();
-        auto unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
-
-        // enter vertices into unstructured grid
-        for (int i = 0; i < n_ranks; ++i) {
-          for (int v = 0; v < nVertices; ++v) {
-            vtkpoints->InsertNextPoint(
-                global_vertices[i * (nVertices * balancer->getDimension() + 1) +
-                                v * balancer->getDimension() + 0],
-                global_vertices[i * (nVertices * balancer->getDimension() + 1) +
-                                v * balancer->getDimension() + 1],
-                global_vertices[i * (nVertices * balancer->getDimension() + 1) +
-                                v * balancer->getDimension() + 2]);
-          }
-        }
-        unstructuredGrid->SetPoints(vtkpoints);
-
-        // data arrays for work and cell id
-        auto work = vtkSmartPointer<vtkFloatArray>::New();
-        auto cell = vtkSmartPointer<vtkFloatArray>::New();
-        work->SetNumberOfComponents(1);
-        work->SetNumberOfTuples(n_ranks);
-        work->SetName("work");
-        cell->SetNumberOfComponents(1);
-        cell->SetNumberOfTuples(n_ranks);
-        cell->SetName("cell id");
-
-        for (int n = 0; n < n_ranks; ++n) {
-          // define grid points, i.e. vertices of local domain
-          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};
-
-          auto faces = vtkSmartPointer<vtkCellArray>::New();
-          // setup faces of polyhedron
-          vtkIdType f0[3] = {8 * n + 0, 8 * n + 2, 8 * n + 1};
-          vtkIdType f1[3] = {8 * n + 1, 8 * n + 2, 8 * n + 3};
-
-          vtkIdType f2[3] = {8 * n + 0, 8 * n + 4, 8 * n + 2};
-          vtkIdType f3[3] = {8 * n + 2, 8 * n + 4, 8 * n + 6};
-
-          vtkIdType f4[3] = {8 * n + 2, 8 * n + 6, 8 * n + 3};
-          vtkIdType f5[3] = {8 * n + 3, 8 * n + 6, 8 * n + 7};
-
-          vtkIdType f6[3] = {8 * n + 1, 8 * n + 5, 8 * n + 3};
-          vtkIdType f7[3] = {8 * n + 3, 8 * n + 5, 8 * n + 7};
-
-          vtkIdType f8[3] = {8 * n + 0, 8 * n + 4, 8 * n + 1};
-          vtkIdType f9[3] = {8 * n + 1, 8 * n + 4, 8 * n + 5};
-
-          vtkIdType fa[3] = {8 * n + 4, 8 * n + 6, 8 * n + 5};
-          vtkIdType fb[3] = {8 * n + 5, 8 * n + 6, 8 * n + 7};
-
-          faces->InsertNextCell(3, f0);
-          faces->InsertNextCell(3, f1);
-          faces->InsertNextCell(3, f2);
-          faces->InsertNextCell(3, f3);
-          faces->InsertNextCell(3, f4);
-          faces->InsertNextCell(3, f5);
-          faces->InsertNextCell(3, f6);
-          faces->InsertNextCell(3, f7);
-          faces->InsertNextCell(3, f8);
-          faces->InsertNextCell(3, f9);
-          faces->InsertNextCell(3, fa);
-          faces->InsertNextCell(3, fb);
-
-          unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
-                                           faces->GetPointer());
-          work->SetValue(
-              n,
-              global_vertices[n * (nVertices * balancer->getDimension() + 1) +
-                              8 * balancer->getDimension()]);
-          cell->SetValue(n, (T)n);
+    // collect all works and vertices on a single process
+    MPI_Gather(local_vertices, nVertices * balancer->getDimension() + 1,
+               MPIDataTypeT, global_vertices,
+               nVertices * balancer->getDimension() + 1, MPIDataTypeT, 0,
+               communicator);
+
+    if (local_rank == 0) {
+      auto vtkpoints = vtkSmartPointer<vtkPoints>::New();
+      auto unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
+
+      // enter vertices into unstructured grid
+      for (int i = 0; i < n_ranks; ++i) {
+        for (int v = 0; v < nVertices; ++v) {
+          vtkpoints->InsertNextPoint(
+              global_vertices[i * (nVertices * balancer->getDimension() + 1) +
+                              v * balancer->getDimension() + 0],
+              global_vertices[i * (nVertices * balancer->getDimension() + 1) +
+                              v * balancer->getDimension() + 1],
+              global_vertices[i * (nVertices * balancer->getDimension() + 1) +
+                              v * balancer->getDimension() + 2]);
         }
-        unstructuredGrid->GetCellData()->AddArray(work);
-        unstructuredGrid->GetCellData()->AddArray(cell);
-
-        createDirectory("vtk_vertices");
-        std::ostringstream filename;
-        filename << "vtk_vertices/ALL_vtk_vertices_" << std::setw(7)
-                 << std::setfill('0') << step << ".vtu";
-        auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
-        writer->SetInputData(unstructuredGrid);
-        writer->SetFileName(filename.str().c_str());
-        writer->SetDataModeToAscii();
-        // writer->SetDataModeToBinary();
-        writer->Write();
-
-        delete[] global_vertices;
       }
-
-    } catch (CustomException &e) {
-      std::cout << e.what() << std::endl;
-      MPI_Abort(MPI_COMM_WORLD, -1);
+      unstructuredGrid->SetPoints(vtkpoints);
+
+      // data arrays for work and cell id
+      auto work = vtkSmartPointer<vtkFloatArray>::New();
+      auto cell = vtkSmartPointer<vtkFloatArray>::New();
+      work->SetNumberOfComponents(1);
+      work->SetNumberOfTuples(n_ranks);
+      work->SetName("work");
+      cell->SetNumberOfComponents(1);
+      cell->SetNumberOfTuples(n_ranks);
+      cell->SetName("cell id");
+
+      for (int n = 0; n < n_ranks; ++n) {
+        // define grid points, i.e. vertices of local domain
+        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};
+
+        auto faces = vtkSmartPointer<vtkCellArray>::New();
+        // setup faces of polyhedron
+        vtkIdType f0[3] = {8 * n + 0, 8 * n + 2, 8 * n + 1};
+        vtkIdType f1[3] = {8 * n + 1, 8 * n + 2, 8 * n + 3};
+
+        vtkIdType f2[3] = {8 * n + 0, 8 * n + 4, 8 * n + 2};
+        vtkIdType f3[3] = {8 * n + 2, 8 * n + 4, 8 * n + 6};
+
+        vtkIdType f4[3] = {8 * n + 2, 8 * n + 6, 8 * n + 3};
+        vtkIdType f5[3] = {8 * n + 3, 8 * n + 6, 8 * n + 7};
+
+        vtkIdType f6[3] = {8 * n + 1, 8 * n + 5, 8 * n + 3};
+        vtkIdType f7[3] = {8 * n + 3, 8 * n + 5, 8 * n + 7};
+
+        vtkIdType f8[3] = {8 * n + 0, 8 * n + 4, 8 * n + 1};
+        vtkIdType f9[3] = {8 * n + 1, 8 * n + 4, 8 * n + 5};
+
+        vtkIdType fa[3] = {8 * n + 4, 8 * n + 6, 8 * n + 5};
+        vtkIdType fb[3] = {8 * n + 5, 8 * n + 6, 8 * n + 7};
+
+        faces->InsertNextCell(3, f0);
+        faces->InsertNextCell(3, f1);
+        faces->InsertNextCell(3, f2);
+        faces->InsertNextCell(3, f3);
+        faces->InsertNextCell(3, f4);
+        faces->InsertNextCell(3, f5);
+        faces->InsertNextCell(3, f6);
+        faces->InsertNextCell(3, f7);
+        faces->InsertNextCell(3, f8);
+        faces->InsertNextCell(3, f9);
+        faces->InsertNextCell(3, fa);
+        faces->InsertNextCell(3, fb);
+
+        unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
+                                         faces->GetPointer());
+        work->SetValue(
+            n,
+            global_vertices[n * (nVertices * balancer->getDimension() + 1) +
+                            8 * balancer->getDimension()]);
+        cell->SetValue(n, (T)n);
+      }
+      unstructuredGrid->GetCellData()->AddArray(work);
+      unstructuredGrid->GetCellData()->AddArray(cell);
+
+      createDirectory("vtk_vertices");
+      std::ostringstream filename;
+      filename << "vtk_vertices/ALL_vtk_vertices_" << std::setw(7)
+               << std::setfill('0') << step << ".vtu";
+      auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+      writer->SetInputData(unstructuredGrid);
+      writer->SetFileName(filename.str().c_str());
+      writer->SetDataModeToAscii();
+      // writer->SetDataModeToBinary();
+      writer->Write();
+
+      delete[] global_vertices;
     }
   }
 #endif