diff --git a/sources/mesh.F90 b/sources/mesh.F90
index c6068aa..692e1dc 100644
--- a/sources/mesh.F90
+++ b/sources/mesh.F90
@@ -592,7 +592,7 @@ module mesh
     use blocks  , only : append_datablock, remove_datablock, link_blocks
     use helpers , only : print_message
     use mpitools, only : check_status
-    use mpitools, only : npmax, nproc, nodes, lprocs
+    use mpitools, only : nproc, nodes, lprocs
     use mpitools, only : send_array, receive_array
 #endif /* MPI */
 
@@ -601,16 +601,15 @@ module mesh
     integer, intent(out) :: status
 
 #ifdef MPI
-    logical         :: flag
-    integer(kind=4) :: np, nl, nc, nr
+    integer(kind=4) :: l, m, n, p
 
     type(block_meta), pointer :: pmeta
     type(block_data), pointer :: pdata
 
     integer(kind=4) :: tag1, tag2
 
-    integer(kind=4), dimension(nodes)        :: nb
-    integer(kind=4), dimension(lprocs,nodes) :: lb
+    integer(kind=4), dimension(nodes)        :: leafs_per_node
+    integer(kind=4), dimension(lprocs,nodes) :: leafs_per_process
 
     character(len=*), parameter :: loc = 'MESH::redistribute_blocks()'
 #endif /* MPI */
@@ -620,50 +619,53 @@ module mesh
     status = 0
 
 #ifdef MPI
-! calculate the new data block distribution between processes
+! calculate the new data block distribution between nodes and processes on each
+! node
 !
-    nl       = mod(get_nleafs(), nodes)
-    nb(:)    = get_nleafs() / nodes
-    nb(1:nl) = nb(1:nl) + 1
-    do nc = 1, nodes
-      nl          = mod(nb(nc), lprocs)
-      lb( :  ,nc) = nb(nc) / lprocs
-      lb(1:nl,nc) = lb(1:nl,nc) + 1
+    n                   = mod(get_nleafs(), nodes)
+    leafs_per_node( : ) = get_nleafs() / nodes
+    leafs_per_node(1:n) = leafs_per_node(1:n) + 1
+    do n = 1, nodes
+      m                        = mod(leafs_per_node(n), lprocs)
+      leafs_per_process( : ,n) = leafs_per_node(n) / lprocs
+      leafs_per_process(1:m,n) = leafs_per_process(1:m,n) + 1
     end do
 
-    nc = 1
-    nr = 1
-    np = 0
-    nl = 0
+! update the %process field of leafs with respect to the new division,
+! and redistribute data blocks at the same time
+!
+    p = 0 ! process count
+    l = 0 ! leaf per process count
+    m = 1 ! node process count
+    n = 1 ! node count
 
     tag2 = 0
+
     pmeta => list_meta
     do while (associated(pmeta))
 
-! exchange data blocks between processes
-!
-      if (pmeta%process /= np) then
-        if (pmeta%leaf) then
+      if (pmeta%leaf) then
+        l = l + 1
+        if (pmeta%process /= p) then
           changed = .true.
 
           tag1 = tag2 + 1
           tag2 = tag1 + 1
 
-          if (nproc == pmeta%process) then
+          if (nproc == pmeta%process .and. associated(pmeta%data)) then
+            pdata => pmeta%data
+            call send_array(p, tag1, pdata%uu)
+            call send_array(p, tag2, pdata%q)
 
-            call send_array(np, tag1, pmeta%data%uu)
-            call send_array(np, tag2, pmeta%data%q)
-
-            call remove_datablock(pmeta%data, status)
+            call remove_datablock(pdata, status)
             if (status /= 0) then
               call print_message(loc, "Could not remove data block!")
               go to 1000
             end if
-
+            nullify(pmeta%data)
           end if
 
-          if (nproc == np) then
-
+          if (nproc == p .and. .not. associated(pmeta%data)) then
             call append_datablock(pdata, status)
             if (status /= 0) then
               call print_message(loc, "Could not append new data block!")
@@ -673,43 +675,20 @@ module mesh
 
             call receive_array(pmeta%process, tag1, pmeta%data%uu)
             call receive_array(pmeta%process, tag2, pmeta%data%q)
-
           end if
+
+          pmeta%process = p
         end if
 
-        pmeta%process = np
-
-      end if
-
-! determine the new block distribution
-!
-      if (pmeta%leaf) then
-        nl = nl + 1
-
-        if (nl >= lb(nr,nc)) then
-
-          np = min(npmax, np + 1)
-          nl = 0
-          nr = nr + 1
-          if (nr > lprocs) then
-            nr = 1
-            nc = nc + 1
+        if (l >= leafs_per_process(m,n)) then
+          l = 0
+          p = p + 1
+          m = m + 1
+          if (m > lprocs) then
+            m = 1
+            n = n + 1
           end if
-          flag = nc <= nodes
-          if (flag) flag = lb(nr,nc) == 0
-          do while(flag)
-            np = min(npmax, np + 1)
-            nr = nr + 1
-            if (nr > lprocs) then
-              nr = 1
-              nc = nc + 1
-            end if
-            flag = nc <= nodes
-            if (flag) flag = lb(nr,nc) == 0
-          end do
-
         end if
-
       end if
 
       pmeta => pmeta%next