Merge lp:~nickpapior/libgridxc/reshape into lp:libgridxc

[Nick Papior]

Changes all reshape calls in mesh3d to an intrinsic copy.

The test4 has been tested with mpirun -n 4/5/6 and results are the same.

I am not sure whether all array_* calls are executed due to various present-statements in the routines.

[Alberto Garcia (albertog)]

OK. I always wondered whether reshape was safe to use...

review: Approve

[Nick Papior (nickpapior)]

Great, I will commence fixing this in SIESTA-4.0.

=== added file 'src/array.F90'
--- src/array.F90	1970-01-01 00:00:00 +0000
+++ src/array.F90	2017-02-16 10:29:06 +0000
@@ -0,0 +1,974 @@
+!!@LICENSE
+!
+!******************************************************************************
+! MODULE array
+! Provides utilities for copying data from one array shape to another
+! array shape.
+! Written by Nick R. Papior, Feb.2017
+!******************************************************************************
+!
+!   PUBLIC procedures available from this module:
+! array_copy    : Copies an array from one dimension to another array with another dimension
+! array_add     : Adds an array from one dimension to another array with another dimension
+!
+!   PUBLIC parameters, types, and variables available from this module:
+! none
+!
+! The intent of this module is to bypass the use of the intrinsic reshape
+! function.
+! Sadly many compilers implement the reshape with a guard on pointers such
+! that a temporary array is created upon reshaping. For very large segments
+! this turns out to be a substantial memory consumption on the heap which
+! may easily be circumvented by doing explicit copying.
+! These routines are extremely simple and does nothing but copy data from
+! one shape to another.
+!
+! Currently only these copies and adds are implemented:
+!
+!   1D -> 2D [integer, real(sp), real(dp)]
+!   1D -> 3D [integer, real(sp), real(dp)]
+!   2D -> 1D [integer, real(sp), real(dp)]
+!   3D -> 1D [integer, real(sp), real(dp)]
+
+module m_array
+
+  use precision, only: sp, dp
+  use sys,       only: die       ! Termination routine
+
+  private
+
+  public :: array_copy
+  public :: array_add
+
+  interface array_copy
+     module procedure ac_1d_2d_ip, ac_1d_2d_sp, ac_1d_2d_dp
+     module procedure ac_1d_3d_ip, ac_1d_3d_sp, ac_1d_3d_dp
+     module procedure ac_2d_1d_ip, ac_2d_1d_sp, ac_2d_1d_dp
+     module procedure ac_3d_1d_ip, ac_3d_1d_sp, ac_3d_1d_dp
+     module procedure ac_4d_1d_ip, ac_4d_1d_sp, ac_4d_1d_dp
+  end interface array_copy
+
+  interface array_add
+     module procedure aa_1d_2d_ip, aa_1d_2d_sp, aa_1d_2d_dp
+     module procedure aa_1d_3d_ip, aa_1d_3d_sp, aa_1d_3d_dp
+     module procedure aa_1d_4d_ip, aa_1d_4d_sp, aa_1d_4d_dp
+     module procedure aa_2d_1d_ip, aa_2d_1d_sp, aa_2d_1d_dp
+     module procedure aa_3d_1d_ip, aa_3d_1d_sp, aa_3d_1d_dp
+     module procedure aa_4d_1d_ip, aa_4d_1d_sp, aa_4d_1d_dp
+  end interface array_add
+
+contains
+
+  ! Copies a 1D array to a 2D array
+  subroutine ac_1d_2d_ip(from_1D, to_1D, in1D, from_2D, to_2D, out2D)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(in) :: in1D(:)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    integer, intent(inout) :: out2D(:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i2D, j2D
+
+    i2D = from_2D(1)
+    j2D = from_2D(2)
+    do i1D = from_1D, to_1D
+       out2D(i2D, j2D) = in1D(i1D)
+       i2D = i2D + 1
+       if ( i2D > to_2D(1) ) then
+          i2D = from_2D(1)
+          j2D = j2D + 1
+       end if
+    end do
+
+    if ( i2D /= from_2D(1) ) &
+         call die("integer: 1D->2D failed (i)")
+    if ( j2D <= to_2D(2) ) &
+         call die("integer: 1D->2D failed (j)")
+
+  end subroutine ac_1d_2d_ip
+  subroutine ac_1d_2d_sp(from_1D, to_1D, in1D, from_2D, to_2D, out2D)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(in) :: in1D(:)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(sp), intent(inout) :: out2D(:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i2D, j2D
+
+    i2D = from_2D(1)
+    j2D = from_2D(2)
+    do i1D = from_1D, to_1D
+       out2D(i2D, j2D) = in1D(i1D)
+       i2D = i2D + 1
+       if ( i2D > to_2D(1) ) then
+          i2D = from_2D(1)
+          j2D = j2D + 1
+       end if
+    end do
+
+    if ( i2D /= from_2D(1) ) &
+         call die("real: 1D->2D failed (i)")
+    if ( j2D <= to_2D(2) ) &
+         call die("real: 1D->2D failed (j)")
+
+  end subroutine ac_1d_2d_sp
+  subroutine ac_1d_2d_dp(from_1D, to_1D, in1D, from_2D, to_2D, out2D)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(in) :: in1D(:)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(dp), intent(inout) :: out2D(:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i2D, j2D
+
+    i2D = from_2D(1)
+    j2D = from_2D(2)
+    do i1D = from_1D, to_1D
+       out2D(i2D, j2D) = in1D(i1D)
+       i2D = i2D + 1
+       if ( i2D > to_2D(1) ) then
+          i2D = from_2D(1)
+          j2D = j2D + 1
+       end if
+    end do
+
+    if ( i2D /= from_2D(1) ) &
+         call die("double: 1D->2D failed (i)")
+    if ( j2D <= to_2D(2) ) &
+         call die("double: 1D->2D failed (j)")
+
+  end subroutine ac_1d_2d_dp
+
+  ! Copies a 1D array to a 3D array
+  subroutine ac_1d_3d_ip(from_1D, to_1D, in1D, from_3D, to_3D, out3D)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(in) :: in1D(:)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    integer, intent(inout) :: out3D(:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i3D, j3D, k3D
+
+    i3D = from_3D(1)
+    j3D = from_3D(2)
+    k3D = from_3D(3)
+    do i1D = from_1D, to_1D
+       out3D(i3D, j3D, k3D) = in1D(i1D)
+       i3D = i3D + 1
+       if ( i3D > to_3D(1) ) then
+          i3D = from_3D(1)
+          j3D = j3D + 1
+       end if
+       if ( j3D > to_3D(2) ) then
+          j3D = from_3D(2)
+          k3D = k3D + 1
+       end if
+    end do
+
+    if ( i3D /= from_3D(1) ) &
+         call die("integer: 1D->3D failed (i)")
+    if ( j3D /= from_3D(2) ) &
+         call die("integer: 1D->3D failed (j)")
+    if ( k3D <= to_3D(3) ) &
+         call die("integer: 1D->3D failed (k)")
+
+  end subroutine ac_1d_3d_ip
+  subroutine ac_1d_3d_sp(from_1D, to_1D, in1D, from_3D, to_3D, out3D)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(in) :: in1D(:)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(sp), intent(inout) :: out3D(:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i3D, j3D, k3D
+
+    i3D = from_3D(1)
+    j3D = from_3D(2)
+    k3D = from_3D(3)
+    do i1D = from_1D, to_1D
+       out3D(i3D, j3D, k3D) = in1D(i1D)
+       i3D = i3D + 1
+       if ( i3D > to_3D(1) ) then
+          i3D = from_3D(1)
+          j3D = j3D + 1
+       end if
+       if ( j3D > to_3D(2) ) then
+          j3D = from_3D(2)
+          k3D = k3D + 1
+       end if
+    end do
+
+    if ( i3D /= from_3D(1) ) &
+         call die("real: 1D->3D failed (i)")
+    if ( j3D /= from_3D(2) ) &
+         call die("real: 1D->3D failed (j)")
+    if ( k3D <= to_3D(3) ) &
+         call die("real: 1D->3D failed (k)")
+
+  end subroutine ac_1d_3d_sp
+  subroutine ac_1d_3d_dp(from_1D, to_1D, in1D, from_3D, to_3D, out3D)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(in) :: in1D(:)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(dp), intent(inout) :: out3D(:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i3D, j3D, k3D
+
+    i3D = from_3D(1)
+    j3D = from_3D(2)
+    k3D = from_3D(3)
+    do i1D = from_1D, to_1D
+       out3D(i3D, j3D, k3D) = in1D(i1D)
+       i3D = i3D + 1
+       if ( i3D > to_3D(1) ) then
+          i3D = from_3D(1)
+          j3D = j3D + 1
+       end if
+       if ( j3D > to_3D(2) ) then
+          j3D = from_3D(2)
+          k3D = k3D + 1
+       end if
+    end do
+
+    if ( i3D /= from_3D(1) ) &
+         call die("double: 1D->3D failed (i)")
+    if ( j3D /= from_3D(2) ) &
+         call die("double: 1D->3D failed (j)")
+    if ( k3D <= to_3D(3) ) &
+         call die("double: 1D->3D failed (k)")
+
+  end subroutine ac_1d_3d_dp
+
+  ! Copies a 2D array to a 1D array
+  subroutine ac_2d_1d_ip(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    integer, intent(in) :: in2D(:,:)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i2D, j2D, i1D
+
+    i1D = from_1D
+    do j2D = from_2D(2), to_2D(2)
+       do i2D = from_2D(1), to_2D(1)
+          out1D(i1D) = in2D(i2D, j2D)
+          i1D = i1D + 1
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("integer: 2D->1D failed")
+
+  end subroutine ac_2d_1d_ip
+  subroutine ac_2d_1d_sp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(sp), intent(in) :: in2D(:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i2D, j2D, i1D
+
+    i1D = from_1D
+    do j2D = from_2D(2), to_2D(2)
+       do i2D = from_2D(1), to_2D(1)
+          out1D(i1D) = in2D(i2D, j2D)
+          i1D = i1D + 1
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("real: 2D->1D failed")
+
+  end subroutine ac_2d_1d_sp
+  subroutine ac_2d_1d_dp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(dp), intent(in) :: in2D(:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i2D, j2D, i1D
+
+    i1D = from_1D
+    do j2D = from_2D(2), to_2D(2)
+       do i2D = from_2D(1), to_2D(1)
+          out1D(i1D) = in2D(i2D, j2D)
+          i1D = i1D + 1
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("double: 2D->1D failed")
+
+  end subroutine ac_2d_1d_dp
+  
+  ! Copies a 3D array to a 1D array
+  subroutine ac_3d_1d_ip(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    integer, intent(in) :: in3D(:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i3D, j3D, k3D, i1D
+
+    i1D = from_1D
+    do k3D = from_3D(3), to_3D(3)
+       do j3D = from_3D(2), to_3D(2)
+          do i3D = from_3D(1), to_3D(1)
+             out1D(i1D) = in3D(i3D, j3D, k3D)
+             i1D = i1D + 1
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("integer: 3D->1D failed")
+
+  end subroutine ac_3d_1d_ip
+  subroutine ac_3d_1d_sp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(sp), intent(in) :: in3D(:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i3D, j3D, k3D, i1D
+
+    i1D = from_1D
+    do k3D = from_3D(3), to_3D(3)
+       do j3D = from_3D(2), to_3D(2)
+          do i3D = from_3D(1), to_3D(1)
+             out1D(i1D) = in3D(i3D, j3D, k3D)
+             i1D = i1D + 1
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("real: 3D->1D failed")
+
+  end subroutine ac_3d_1d_sp
+  subroutine ac_3d_1d_dp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(dp), intent(in) :: in3D(:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i3D, j3D, k3D, i1D
+
+    i1D = from_1D
+    do k3D = from_3D(3), to_3D(3)
+       do j3D = from_3D(2), to_3D(2)
+          do i3D = from_3D(1), to_3D(1)
+             out1D(i1D) = in3D(i3D, j3D, k3D)
+             i1D = i1D + 1
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("double: 3D->1D failed")
+
+  end subroutine ac_3d_1d_dp
+
+  ! Copies a 4D array to a 1D array
+  subroutine ac_4d_1d_ip(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    integer, intent(in) :: in4D(:,:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i4D, j4D, k4D, m4D, i1D
+
+    i1D = from_1D
+    do m4D = from_4D(4), to_4D(4)
+       do k4D = from_4D(3), to_4D(3)
+          do j4D = from_4D(2), to_4D(2)
+             do i4D = from_4D(1), to_4D(1)
+                out1D(i1D) = in4D(i4D, j4D, k4D, m4D)
+                i1D = i1D + 1
+             end do
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("integer: 4D->1D failed")
+
+  end subroutine ac_4d_1d_ip
+  subroutine ac_4d_1d_sp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    real(sp), intent(in) :: in4D(:,:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i4D, j4D, k4D, m4D, i1D
+
+    i1D = from_1D
+    do m4D = from_4D(4), to_4D(4)
+       do k4D = from_4D(3), to_4D(3)
+          do j4D = from_4D(2), to_4D(2)
+             do i4D = from_4D(1), to_4D(1)
+                out1D(i1D) = in4D(i4D, j4D, k4D, m4D)
+                i1D = i1D + 1
+             end do
+          end do
+       end do
+    end do
+    
+    if ( i1D <= to_1D ) &
+         call die("real: 4D->1D failed")
+
+  end subroutine ac_4d_1d_sp
+  subroutine ac_4d_1d_dp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    real(dp), intent(in) :: in4D(:,:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i4D, j4D, k4D, m4D, i1D
+
+    i1D = from_1D
+    do m4D = from_4D(4), to_4D(4)
+       do k4D = from_4D(3), to_4D(3)
+          do j4D = from_4D(2), to_4D(2)
+             do i4D = from_4D(1), to_4D(1)
+                out1D(i1D) = in4D(i4D, j4D, k4D, m4D)
+                i1D = i1D + 1
+             end do
+          end do
+       end do
+    end do
+    
+    if ( i1D <= to_1D ) &
+         call die("double: 4D->1D failed")
+
+  end subroutine ac_4d_1d_dp
+  
+
+  ! Adds a 1D array to a 2D array
+  subroutine aa_1d_2d_ip(from_1D, to_1D, a1D, from_2D, to_2D, out2D)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(in) :: a1D(:)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    integer, intent(inout) :: out2D(:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i2D, j2D
+
+    i2D = from_2D(1)
+    j2D = from_2D(2)
+    do i1D = from_1D, to_1D
+       out2D(i2D, j2D) = out2D(i2D, j2D) + a1D(i1D)
+       i2D = i2D + 1
+       if ( i2D > to_2D(1) ) then
+          i2D = from_2D(1)
+          j2D = j2D + 1
+       end if
+    end do
+
+    if ( i2D /= from_2D(1) ) &
+         call die("integer: 1D+>2D failed (i)")
+    if ( j2D <= to_2D(2) ) &
+         call die("integer: 1D+>2D failed (j)")
+
+  end subroutine aa_1d_2d_ip
+  subroutine aa_1d_2d_sp(from_1D, to_1D, a1D, from_2D, to_2D, out2D)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(in) :: a1D(:)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(sp), intent(inout) :: out2D(:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i2D, j2D
+
+    i2D = from_2D(1)
+    j2D = from_2D(2)
+    do i1D = from_1D, to_1D
+       out2D(i2D, j2D) = out2D(i2D, j2D) + a1D(i1D)
+       i2D = i2D + 1
+       if ( i2D > to_2D(1) ) then
+          i2D = from_2D(1)
+          j2D = j2D + 1
+       end if
+    end do
+
+    if ( i2D /= from_2D(1) ) &
+         call die("real: 1D+>2D failed (i)")
+    if ( j2D <= to_2D(2) ) &
+         call die("real: 1D+>2D failed (j)")
+
+  end subroutine aa_1d_2d_sp
+  subroutine aa_1d_2d_dp(from_1D, to_1D, a1D, from_2D, to_2D, out2D)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(in) :: a1D(:)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(dp), intent(inout) :: out2D(:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i2D, j2D
+
+    i2D = from_2D(1)
+    j2D = from_2D(2)
+    do i1D = from_1D, to_1D
+       out2D(i2D, j2D) = out2D(i2D, j2D) + a1D(i1D)
+       i2D = i2D + 1
+       if ( i2D > to_2D(1) ) then
+          i2D = from_2D(1)
+          j2D = j2D + 1
+       end if
+    end do
+
+    if ( i2D /= from_2D(1) ) &
+         call die("double: 1D+>2D failed (i)")
+    if ( j2D <= to_2D(2) ) &
+         call die("double: 1D+>2D failed (j)")
+
+  end subroutine aa_1d_2d_dp
+
+  ! Adds a 1D array to a 3D array
+  subroutine aa_1d_3d_ip(from_1D, to_1D, a1D, from_3D, to_3D, out3D)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(in) :: a1D(:)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    integer, intent(inout) :: out3D(:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i3D, j3D, k3D
+
+    i3D = from_3D(1)
+    j3D = from_3D(2)
+    k3D = from_3D(3)
+    do i1D = from_1D, to_1D
+       out3D(i3D, j3D, k3D) = out3D(i3D, j3D, k3D) + a1D(i1D)
+       i3D = i3D + 1
+       if ( i3D > to_3D(1) ) then
+          i3D = from_3D(1)
+          j3D = j3D + 1
+       end if
+       if ( j3D > to_3D(2) ) then
+          j3D = from_3D(2)
+          k3D = k3D + 1
+       end if
+    end do
+
+    if ( i3D /= from_3D(1) ) &
+         call die("integer: 1D+>3D failed (i)")
+    if ( j3D /= from_3D(2) ) &
+         call die("integer: 1D+>3D failed (j)")
+    if ( k3D <= to_3D(3) ) &
+         call die("integer: 1D+>3D failed (k)")
+
+  end subroutine aa_1d_3d_ip
+  subroutine aa_1d_3d_sp(from_1D, to_1D, a1D, from_3D, to_3D, out3D)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(in) :: a1D(:)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(sp), intent(inout) :: out3D(:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i3D, j3D, k3D
+
+    i3D = from_3D(1)
+    j3D = from_3D(2)
+    k3D = from_3D(3)
+    do i1D = from_1D, to_1D
+       out3D(i3D, j3D, k3D) = out3D(i3D, j3D, k3D) + a1D(i1D)
+       i3D = i3D + 1
+       if ( i3D > to_3D(1) ) then
+          i3D = from_3D(1)
+          j3D = j3D + 1
+       end if
+       if ( j3D > to_3D(2) ) then
+          j3D = from_3D(2)
+          k3D = k3D + 1
+       end if
+    end do
+
+    if ( i3D /= from_3D(1) ) &
+         call die("real: 1D+>3D failed (i)")
+    if ( j3D /= from_3D(2) ) &
+         call die("real: 1D+>3D failed (j)")
+    if ( k3D <= to_3D(3) ) &
+         call die("real: 1D+>3D failed (k)")
+
+  end subroutine aa_1d_3d_sp
+  subroutine aa_1d_3d_dp(from_1D, to_1D, a1D, from_3D, to_3D, out3D)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(in) :: a1D(:)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(dp), intent(inout) :: out3D(:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i3D, j3D, k3D
+
+    i3D = from_3D(1)
+    j3D = from_3D(2)
+    k3D = from_3D(3)
+    do i1D = from_1D, to_1D
+       out3D(i3D, j3D, k3D) = out3D(i3D, j3D, k3D) + a1D(i1D)
+       i3D = i3D + 1
+       if ( i3D > to_3D(1) ) then
+          i3D = from_3D(1)
+          j3D = j3D + 1
+       end if
+       if ( j3D > to_3D(2) ) then
+          j3D = from_3D(2)
+          k3D = k3D + 1
+       end if
+    end do
+
+    if ( i3D /= from_3D(1) ) &
+         call die("double: 1D+>3D failed (i)")
+    if ( j3D /= from_3D(2) ) &
+         call die("double: 1D+>3D failed (j)")
+    if ( k3D <= to_3D(3) ) &
+         call die("double: 1D+>3D failed (k)")
+
+  end subroutine aa_1d_3d_dp
+
+  ! Adds a 1D array to a 4D array
+  subroutine aa_1d_4d_ip(from_1D, to_1D, a1D, from_4D, to_4D, out4D)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(in) :: a1D(:)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    integer, intent(inout) :: out4D(:,:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i4D, j4D, k4D, m4D
+
+    i4D = from_4D(1)
+    j4D = from_4D(2)
+    k4D = from_4D(3)
+    m4D = from_4D(4)
+    do i1D = from_1D, to_1D
+       out4D(i4D, j4D, k4D, m4D) = out4D(i4D, j4D, k4D, m4D) + a1D(i1D)
+       i4D = i4D + 1
+       if ( i4D > to_4D(1) ) then
+          i4D = from_4D(1)
+          j4D = j4D + 1
+       end if
+       if ( j4D > to_4D(2) ) then
+          j4D = from_4D(2)
+          k4D = k4D + 1
+       end if
+       if ( k4D > to_4D(3) ) then
+          k4D = from_4D(3)
+          m4D = m4D + 1
+       end if
+    end do
+
+    if ( i4D /= from_4D(1) ) &
+         call die("integer: 1D+>4D failed (i)")
+    if ( j4D /= from_4D(2) ) &
+         call die("integer: 1D+>4D failed (j)")
+    if ( k4D /= from_4D(3) ) &
+         call die("integer: 1D+>4D failed (k)")
+    if ( m4D <= to_4D(4) ) &
+         call die("integer: 1D+>4D failed (m)")
+
+  end subroutine aa_1d_4d_ip
+  subroutine aa_1d_4d_sp(from_1D, to_1D, a1D, from_4D, to_4D, out4D)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(in) :: a1D(:)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    real(sp), intent(inout) :: out4D(:,:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i4D, j4D, k4D, m4D
+
+    i4D = from_4D(1)
+    j4D = from_4D(2)
+    k4D = from_4D(3)
+    m4D = from_4D(4)
+    do i1D = from_1D, to_1D
+       out4D(i4D, j4D, k4D, m4D) = out4D(i4D, j4D, k4D, m4D) + a1D(i1D)
+       i4D = i4D + 1
+       if ( i4D > to_4D(1) ) then
+          i4D = from_4D(1)
+          j4D = j4D + 1
+       end if
+       if ( j4D > to_4D(2) ) then
+          j4D = from_4D(2)
+          k4D = k4D + 1
+       end if
+       if ( k4D > to_4D(3) ) then
+          k4D = from_4D(3)
+          m4D = m4D + 1
+       end if
+    end do
+
+    if ( i4D /= from_4D(1) ) &
+         call die("real: 1D+>4D failed (i)")
+    if ( j4D /= from_4D(2) ) &
+         call die("real: 1D+>4D failed (j)")
+    if ( k4D /= from_4D(3) ) &
+         call die("real: 1D+>4D failed (k)")
+    if ( m4D <= to_4D(4) ) &
+         call die("real: 1D+>4D failed (m)")
+
+  end subroutine aa_1d_4d_sp
+  subroutine aa_1d_4d_dp(from_1D, to_1D, a1D, from_4D, to_4D, out4D)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(in) :: a1D(:)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    real(dp), intent(inout) :: out4D(:,:,:,:)
+
+    ! Local variables for copying
+    integer :: i1D, i4D, j4D, k4D, m4D
+
+    i4D = from_4D(1)
+    j4D = from_4D(2)
+    k4D = from_4D(3)
+    m4D = from_4D(4)
+    do i1D = from_1D, to_1D
+       out4D(i4D, j4D, k4D, m4D) = out4D(i4D, j4D, k4D, m4D) + a1D(i1D)
+       i4D = i4D + 1
+       if ( i4D > to_4D(1) ) then
+          i4D = from_4D(1)
+          j4D = j4D + 1
+       end if
+       if ( j4D > to_4D(2) ) then
+          j4D = from_4D(2)
+          k4D = k4D + 1
+       end if
+       if ( k4D > to_4D(3) ) then
+          k4D = from_4D(3)
+          m4D = m4D + 1
+       end if
+    end do
+
+    if ( i4D /= from_4D(1) ) &
+         call die("double: 1D+>4D failed (i)")
+    if ( j4D /= from_4D(2) ) &
+         call die("double: 1D+>4D failed (j)")
+    if ( k4D /= from_4D(3) ) &
+         call die("double: 1D+>4D failed (k)")
+    if ( m4D <= to_4D(4) ) &
+         call die("double: 1D+>4D failed (m)")
+
+  end subroutine aa_1d_4d_dp
+  
+  ! Adds a 2D array to a 1D array
+  subroutine aa_2d_1d_ip(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    integer, intent(in) :: in2D(:,:)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i2D, j2D, i1D
+
+    i1D = from_1D
+    do j2D = from_2D(2), to_2D(2)
+       do i2D = from_2D(1), to_2D(1)
+          out1D(i1D) = out1D(i1D) + in2D(i2D, j2D)
+          i1D = i1D + 1
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("integer: 2D+>1D failed")
+
+  end subroutine aa_2d_1d_ip
+  subroutine aa_2d_1d_sp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(sp), intent(in) :: in2D(:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i2D, j2D, i1D
+
+    i1D = from_1D
+    do j2D = from_2D(2), to_2D(2)
+       do i2D = from_2D(1), to_2D(1)
+          out1D(i1D) = out1D(i1D) + in2D(i2D, j2D)
+          i1D = i1D + 1
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("real: 2D+>1D failed")
+
+  end subroutine aa_2d_1d_sp
+  subroutine aa_2d_1d_dp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_2D(2), to_2D(2)
+    real(dp), intent(in) :: in2D(:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i2D, j2D, i1D
+
+    i1D = from_1D
+    do j2D = from_2D(2), to_2D(2)
+       do i2D = from_2D(1), to_2D(1)
+          out1D(i1D) = out1D(i1D) + in2D(i2D, j2D)
+          i1D = i1D + 1
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("double: 2D+>1D failed")
+
+  end subroutine aa_2d_1d_dp
+  
+  ! Adds a 3D array to a 1D array
+  subroutine aa_3d_1d_ip(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    integer, intent(in) :: in3D(:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i3D, j3D, k3D, i1D
+
+    i1D = from_1D
+    do k3D = from_3D(3), to_3D(3)
+       do j3D = from_3D(2), to_3D(2)
+          do i3D = from_3D(1), to_3D(1)
+             out1D(i1D) = out1D(i1D) + in3D(i3D, j3D, k3D)
+             i1D = i1D + 1
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("integer: 3D+>1D failed")
+
+  end subroutine aa_3d_1d_ip
+  subroutine aa_3d_1d_sp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(sp), intent(in) :: in3D(:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i3D, j3D, k3D, i1D
+
+    i1D = from_1D
+    do k3D = from_3D(3), to_3D(3)
+       do j3D = from_3D(2), to_3D(2)
+          do i3D = from_3D(1), to_3D(1)
+             out1D(i1D) = out1D(i1D) + in3D(i3D, j3D, k3D)
+             i1D = i1D + 1
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("real: 3D+>1D failed")
+
+  end subroutine aa_3d_1d_sp
+  subroutine aa_3d_1d_dp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_3D(3), to_3D(3)
+    real(dp), intent(in) :: in3D(:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i3D, j3D, k3D, i1D
+
+    i1D = from_1D
+    do k3D = from_3D(3), to_3D(3)
+       do j3D = from_3D(2), to_3D(2)
+          do i3D = from_3D(1), to_3D(1)
+             out1D(i1D) = out1D(i1D) + in3D(i3D, j3D, k3D)
+             i1D = i1D + 1
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("double: 3D+>1D failed")
+
+  end subroutine aa_3d_1d_dp
+
+  ! Adds a 4D array to a 1D array
+  subroutine aa_4d_1d_ip(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    integer, intent(in) :: in4D(:,:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    integer, intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i4D, j4D, k4D, m4D, i1D
+
+    i1D = from_1D
+    do m4D = from_4D(4), to_4D(4)
+       do k4D = from_4D(3), to_4D(3)
+          do j4D = from_4D(2), to_4D(2)
+             do i4D = from_4D(1), to_4D(1)
+                out1D(i1D) = out1D(i1D) + in4D(i4D, j4D, k4D, m4D)
+                i1D = i1D + 1
+             end do
+          end do
+       end do
+    end do
+    
+    if ( i1D <= to_1D ) &
+         call die("integer: 4D+>1D failed")
+
+  end subroutine aa_4d_1d_ip
+  subroutine aa_4d_1d_sp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    real(sp), intent(in) :: in4D(:,:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(sp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i4D, j4D, k4D, m4D, i1D
+
+    i1D = from_1D
+    do m4D = from_4D(4), to_4D(4)
+       do k4D = from_4D(3), to_4D(3)
+          do j4D = from_4D(2), to_4D(2)
+             do i4D = from_4D(1), to_4D(1)
+                out1D(i1D) = out1D(i1D) + in4D(i4D, j4D, k4D, m4D)
+                i1D = i1D + 1
+             end do
+          end do
+       end do
+    end do
+    
+    if ( i1D <= to_1D ) &
+         call die("real: 4D+>1D failed")
+
+  end subroutine aa_4d_1d_sp
+  subroutine aa_4d_1d_dp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
+    integer, intent(in) :: from_4D(4), to_4D(4)
+    real(dp), intent(in) :: in4D(:,:,:,:)
+    integer, intent(in) :: from_1D, to_1D
+    real(dp), intent(inout) :: out1D(:)
+
+    ! Local variables for copying
+    integer :: i4D, j4D, k4D, m4D, i1D
+
+    i1D = from_1D
+    do m4D = from_4D(4), to_4D(4)
+       do k4D = from_4D(3), to_4D(3)
+          do j4D = from_4D(2), to_4D(2)
+             do i4D = from_4D(1), to_4D(1)
+                out1D(i1D) = out1D(i1D) + in4D(i4D, j4D, k4D, m4D)
+                i1D = i1D + 1
+             end do
+          end do
+       end do
+    end do
+
+    if ( i1D <= to_1D ) &
+         call die("double: 4D+>1D failed")
+
+  end subroutine aa_4d_1d_dp
+
+end module m_array
+    

=== modified file 'src/gridxc.F90'
--- src/gridxc.F90	2017-01-20 11:55:29 +0000
+++ src/gridxc.F90	2017-02-16 10:29:06 +0000
@@ -605,7 +605,9 @@
   USE m_cellXC, only: gridxc_cellXC => cellXC  ! XC for a periodic unit cell
   USE xcmod,    only: gridxc_getXC => getXC   ! Returns XC functional(s)
   USE xcmod,    only: gridxc_setXC => setXC   ! Sets XC functional(s)
+#ifdef LIBXC
   USE xcmod,    only: gridxc_setXC_libxc => setXC_libxc_ids   ! Sets XC functional(s)
+#endif
 
 ! Secondary entry points for testers and lower-level programming
   USE m_ldaxc,  only: gridxc_ldaxc => ldaxc    ! LDA-XC functionals

=== modified file 'src/makefile'
--- src/makefile	2017-01-20 11:55:29 +0000
+++ src/makefile	2017-02-16 10:29:06 +0000
@@ -29,35 +29,46 @@
 #
 ARCH_MAKE_DEFAULT=$(MAIN_OBJDIR)/fortran.mk
 ARCH_MAKE?=$(ARCH_MAKE_DEFAULT)
-include $(ARCH_MAKE)
+
+# Ensure that we can ALWAYS run make dep in the MAIN directory
+ifeq ($(VPATH),.)
+ -include $(ARCH_MAKE)
+else
+ include $(ARCH_MAKE)
+endif
+
+# Default library commands
+AR ?= ar
 RANLIB ?= ranlib
-#
+
 ifdef WITH_LIBXC
- #
+
  ifndef LIBXC_ROOT
    $(error you need to define LIBXC_ROOT in your arch.make)
  endif
- #
+
  LIBXC_MK_FILE_DIR ?= $(LIBXC_ROOT)
  FPPFLAGS += -DLIBXC
  include $(LIBXC_MK_FILE_DIR)/libxc.mk
+
 endif
-#
+
 module:  libGridXC.a
-#
+
 dep:
 	sfmakedepend --depend=obj  --modext=o *.f *.f90 *.F *.F90
-#
+
 GRID_SRCS= alloc.F90 bessph.f cellsubs.f \
              chkgmx.f  interpolation.f90  \
              minvec.f m_io.f moreParallelSubs.F90 \
              m_fft_gpfa.F  m_walltime.f90 \
-	     precision.F radfft.f sorting.f sys.F
+	     precision.F radfft.f sorting.f sys.F \
+             array.F90
 GRID_OBJS:= $(GRID_SRCS:.f=.o)
 GRID_OBJS:= $(GRID_OBJS:.F=.o)
 GRID_OBJS:= $(GRID_OBJS:.f90=.o)
 GRID_OBJS:= $(GRID_OBJS:.F90=.o)
-#
+
 LOCAL_SRCS= atomxc.F90 cellxc.F90 debugxc.f90 ggaxc.F ldaxc.F \
             gridxc.F90 vdwxc.F90 xcmod.f90  fft3d.F90 fftr.F90 \
             mesh1d.f90 mesh3d.F90 am05.f90 vv_vdwxc.F90 \
@@ -66,12 +77,12 @@
 LOCAL_OBJS:= $(LOCAL_OBJS:.F=.o)
 LOCAL_OBJS:= $(LOCAL_OBJS:.f90=.o)
 LOCAL_OBJS:= $(LOCAL_OBJS:.F90=.o)
-#
+
 ALL_OBJS=$(GRID_OBJS) $(LOCAL_OBJS)
-#
-#
+
+
 INCFLAGS:= $(LIBXC_INCFLAGS) $(INCFLAGS)
-#
+
 ifdef WITH_MPI
   MPI_INTERFACE=libmpi_f90.a
   FPPFLAGS:= -DMPI $(FPPFLAGS)
@@ -79,21 +90,21 @@
 else
   MPI_INTERFACE=
 endif
-#
+
 libmpi_f90.a: 
 	@(cd MPI_instr ; $(MAKE) \
                         "VPATH=$(VPATH)/MPI_instr" )
-#
+
 libGridXC.a: $(MPI_INTERFACE) $(ALL_OBJS)
-	ar $(ARFLAGS_EXTRA) cru libGridXC.a $(ALL_OBJS) 
+	$(AR) $(ARFLAGS_EXTRA) cru libGridXC.a $(ALL_OBJS) 
 	@if [ ! -z "$(MPI_INTERFACE)" ] ; then \
-            ar $(ARFLAGS_EXTRA) cru libGridXC.a MPI_instr/*.o ; fi
+            $(AR) $(ARFLAGS_EXTRA) cru libGridXC.a MPI_instr/*.o ; fi
 	-$(RANLIB) libGridXC.a
-#
+
 MODULES_TO_INSTALL=gridxc.mod gridxc_config.mod m_atomxc.mod \
                    m_cellxc.mod xcmod.mod m_ldaxc.mod m_ggaxc.mod \
                    m_fft_gpfa.mod debugxc.mod mesh3d.mod
-#
+
 ifdef WITH_LIBXC
  ifdef WITH_MPI
 setup_mk_file:
@@ -112,7 +123,7 @@
 	cp gridxc.mk.in gridxc.mk
  endif
 endif
-#
+
 module: setup_mk_file libGridXC.a
 	mkdir -p $(MAIN_OBJDIR)/lib 
 	mkdir -p $(MAIN_OBJDIR)/include
@@ -129,17 +140,20 @@
 	rm -f $(GRIDXC_PREFIX)/gridxc.mk
 
 	cp -p $(MAIN_OBJDIR)/include/* $(GRIDXC_PREFIX)/include
-	cp -p $(MAIN_OBJDIR)/lib/libGridXC.a  $(GRIDXC_PREFIX)/lib
+	cp -p $(MAIN_OBJDIR)/lib/libGridXC.a $(GRIDXC_PREFIX)/lib
 	cp -p gridxc.mk $(GRIDXC_PREFIX)
-#
+
 
 clean:
 	rm -f *.o *.*d *.a gridxc.mk
 	(cd MPI_instr ; $(MAKE) \
                         "VPATH=$(VPATH)/MPI_instr" clean)
+
+
 #
 # DO NOT DELETE THIS LINE - used by make depend
 am05.o: precision.o sys.o
+array.o: precision.o sys.o
 atomxc.o: alloc.o debugxc.o ggaxc.o gridxc_config.o ldaxc.o mesh1d.o
 atomxc.o: precision.o radfft.o sys.o vdwxc.o xcmod.o
 bessph.o: precision.o sys.o
@@ -157,15 +171,16 @@
 ldaxc.o: precision.o sys.o
 m_io.o: sys.o
 mesh1d.o: interpolation.o precision.o
-mesh3d.o: alloc.o debugxc.o gridxc_config.o precision.o sorting.o sys.o
+mesh3d.o: alloc.o array.o debugxc.o gridxc_config.o precision.o sorting.o sys.o
 minvec.o: cellsubs.o precision.o sorting.o sys.o
-moreParallelSubs.o: alloc.o gridxc_config.o m_io.o precision.o sys.o
+moreParallelSubs.o: alloc.o array.o gridxc_config.o m_io.o precision.o sys.o
 radfft.o: alloc.o bessph.o m_fft_gpfa.o precision.o
 vdwxc.o: alloc.o debugxc.o ggaxc.o interpolation.o ldaxc.o mesh1d.o precision.o
 vdwxc.o: radfft.o sys.o vv_vdwxc.o
 vv_vdwxc.o: alloc.o debugxc.o interpolation.o mesh1d.o precision.o radfft.o
 vv_vdwxc.o: sys.o
 xcmod.o: precision.o sys.o vdwxc.o
+m_array.o: array.o
 m_atomxc.o: atomxc.o
 m_bessph.o: bessph.o
 m_cellxc.o: cellxc.o

=== modified file 'src/mesh3d.F90'
--- src/mesh3d.F90	2016-05-09 12:15:45 +0000
+++ src/mesh3d.F90	2017-02-16 10:29:06 +0000
@@ -29,6 +29,8 @@
 ! use alloc,      only: re_alloc        ! Re-allocation routine
 ! use sys,        only: die             ! Termination routine
 ! use sorting,    only: ordix           ! Finds the sorting index of a vector
+! use m_array, only: array_copy         ! Copies flattened arrays (instead of reshaping)
+! use m_array, only: array_add          ! Adds flattened arrays (instead of reshaping)
 !
 !   USED module parameters:
 ! use precision,  only: dp              ! Real double precision type
@@ -648,6 +650,9 @@
   use debugXC,   only: udebug          ! Output file unit for debug info
 #endif /* DEBUG_XC */
 
+  ! Copy and add for limiting temporary arrays
+  use m_array, only: array_copy, array_add
+
 ! Used MPI procedures and types
 #ifdef MPI
   use mpi_instr, only: MPI_AllGather
@@ -1624,17 +1629,10 @@
   integer,intent(in) :: box(2,3)                 ! My node's mesh box
   integer,intent(out):: boxes(2,3,0:totNodes-1)  ! All node's boxes
 
-  integer,allocatable:: buffer(:)
-  integer :: box1d(6)    ! To avoid the creation of a temporary array
-
 #ifdef MPI
 ! Gather the boxes of all nodes
-  allocate( buffer(6*totNodes) )
-  box1d(1:6) = reshape(box,(/6/))
-  call MPI_AllGather( box1d , 6, MPI_Integer, &
-                      buffer, 6, MPI_Integer, COMM, MPIerror )
-  boxes = reshape( buffer, (/2,3,totNodes/) )
-  deallocate( buffer )
+  call MPI_AllGather( box(1,1) , 6, MPI_Integer, &
+                      boxes(1,1,0), 6, MPI_Integer, COMM, MPIerror )
 #else
 ! Copy the only node's box
   boxes(:,:,0) = box(:,:)
@@ -2018,7 +2016,7 @@
   integer:: orderNode1(nNodes), orderNode2(2*nNodes)
   logical:: found
   real(dp):: nTrsfNode(2*nNodes)
-  integer,pointer:: buffer(:)=>null(), inTrsf(:,:)=>null(), &
+  integer, pointer:: inTrsf(:,:)=>null(), &
                     myTrsf(:)=>null(), outTrsf(:,:)=>null()
 
 ! In serial execution, do nothing
@@ -2047,11 +2045,8 @@
   myTrsf(1:nTrsf) = trsfDir(:)*(trsfNode(:)+1)
 
 ! Gather initial transfer sequence of all nodes
-  call re_alloc( buffer, 1,maxTrsf*nNodes, name=myName//'buffer' )
   call MPI_AllGather( myTrsf, maxTrsf, MPI_Integer, &
-                      buffer, maxTrsf, MPI_Integer, COMM, MPIerror )
-  inTrsf = reshape( buffer, (/maxTrsf,nNodes/) )
-  call de_alloc( buffer, name=myName//'buffer' )
+                      inTrsf(1,1), maxTrsf, MPI_Integer, COMM, MPIerror )
 
 ! Find the number of transfers of each node
   nTrsfNode(1:nNodes) = count( inTrsf(:,1:nNodes)/=0, dim=1 )
@@ -2399,7 +2394,7 @@
   character(len=*),parameter:: errHead = myName//'ERROR: '
   integer,allocatable:: dstBoxes(:,:,:), srcBoxes(:,:,:), &
                         trsfDir(:), trsfNode(:)
-  integer :: comBox(2,3), comShape(4), &
+  integer :: comBox(2,3), arrayS(4,2), &
              dstComBox(2,3,maxParts), dstNode, dstSize, &
              i, ic, ix, iNode, iPart, iTask, iTrsf, &
              maxMesh, nData, nParts, nTrsf, partSize, &
@@ -2573,7 +2568,7 @@
   if (present(prjData)) trsfSize = trsfSize + 3*maxMesh*maxParts
 
 ! Allocate buffer to transfer data
-  call re_alloc( trsfBuff, 1,trsfSize, name=myName//'trsfBuff' )
+  call re_alloc( trsfBuff, 1, trsfSize, name=myName//'trsfBuff' )
 
 ! Loop on transfer sequence
   trueTrsf = 0
@@ -2624,12 +2619,18 @@
       ! Copy data from each part of common box to a transfer buffer
       if (present(dstData)) then
         do iPart = 1,nParts
-          partSize = sizeSum(iPart) - sizeSum(iPart-1)
-          trsfBuff(sizeSum(iPart-1)+1:sizeSum(iPart)) =                    &
-            reshape( srcData(srcComBox(1,1,iPart):srcComBox(2,1,iPart),    &
-                             srcComBox(1,2,iPart):srcComBox(2,2,iPart),    &
-                             srcComBox(1,3,iPart):srcComBox(2,3,iPart),:), &
-                     (/partSize/) )
+          !partSize = sizeSum(iPart) - sizeSum(iPart-1)
+          arrayS(1:3,1) = srcComBox(1,:,iPart) + 1
+          arrayS(4,1) = 1
+          arrayS(1:3,2) = srcComBox(2,:,iPart) + 1
+          arrayS(4,2) = nData
+          call array_copy(arrayS(:,1), arrayS(:,2), srcData(:,:,:,:), &
+               sizeSum(iPart-1)+1, sizeSum(iPart), trsfBuff(:) )
+          !trsfBuff(sizeSum(iPart-1)+1:sizeSum(iPart)) =                    &
+          !  reshape( srcData(srcComBox(1,1,iPart):srcComBox(2,1,iPart),    &
+          !                   srcComBox(1,2,iPart):srcComBox(2,2,iPart),    &
+          !                   srcComBox(1,3,iPart):srcComBox(2,3,iPart),:), &
+          !           (/partSize/) )
         end do ! iPart
         trsfSize = sizeSum(nParts)
       else ! (.not.present(dstData))
@@ -2645,9 +2646,11 @@
           ! Find projection of srcData in comBox
           call projection( nMesh, srcBox, srcData(:,:,:,1), comBox, prjBuff )
           ! Copy projection data onto a transfer buffer
+          call array_copy([1, 1], [3, maxMesh], prjBuff(:,:), &
+               trsfSize+1, trsfSize+3*maxMesh, trsfBuff(:))
           ! AG: Work around gfortran bug: specify 0 lbound in reshape
-          trsfBuff(trsfSize+1:trsfSize+3*maxMesh) = &
-             reshape( prjBuff(:,0:), (/3*maxMesh/) )
+          !trsfBuff(trsfSize+1:trsfSize+3*maxMesh) = &
+          !   reshape( prjBuff(:,0:), (/3*maxMesh/) )
           trsfSize = trsfSize + 3*maxMesh
         end do ! iPart
       end if ! (present(prjData))
@@ -2721,15 +2724,21 @@
       ! Copy buffer data for each part of common box to destination array
       if (present(dstData)) then
         do iPart = 1,nParts
-          comShape(1:3) = dstComBox(2,:,iPart) - dstComBox(1,:,iPart) + 1
-          comShape(4) = nData
-          dstData(dstComBox(1,1,iPart):dstComBox(2,1,iPart),     &
-                  dstComBox(1,2,iPart):dstComBox(2,2,iPart),     &
-                  dstComBox(1,3,iPart):dstComBox(2,3,iPart),:) = &
-          dstData(dstComBox(1,1,iPart):dstComBox(2,1,iPart),     &
-                  dstComBox(1,2,iPart):dstComBox(2,2,iPart),     &
-                  dstComBox(1,3,iPart):dstComBox(2,3,iPart),:) + &
-            reshape( trsfBuff(sizeSum(iPart-1)+1:sizeSum(iPart)), comShape )
+          arrayS(1:3,1) = dstComBox(1,:,iPart) + 1
+          arrayS(4,1) = 1
+          arrayS(1:3,2) = dstComBox(2,:,iPart) + 1
+          arrayS(4,2) = nData
+          call array_add(sizeSum(iPart-1)+1, sizeSum(iPart), trsfBuff(:), &
+               arrayS(:,1), arrayS(:,2), dstData(:,:,:,:))
+          !comShape(1:3) = dstComBox(2,:,iPart) - dstComBox(1,:,iPart) + 1
+          !comShape(4) = nData
+          !dstData(dstComBox(1,1,iPart):dstComBox(2,1,iPart),     &
+          !        dstComBox(1,2,iPart):dstComBox(2,2,iPart),     &
+          !        dstComBox(1,3,iPart):dstComBox(2,3,iPart),:) = &
+          !dstData(dstComBox(1,1,iPart):dstComBox(2,1,iPart),     &
+          !        dstComBox(1,2,iPart):dstComBox(2,2,iPart),     &
+          !        dstComBox(1,3,iPart):dstComBox(2,3,iPart),:) + &
+          !  reshape( trsfBuff(sizeSum(iPart-1)+1:sizeSum(iPart)), comShape )
         end do ! iPart
         trsfSize = sizeSum(nParts)
       else ! (.not.present(dstData))
@@ -2743,8 +2752,10 @@
           comBox(1,:) = srcComBox(1,:,iPart) + srcBox(1,:)
           comBox(2,:) = srcComBox(2,:,iPart) + srcBox(1,:)
           ! Copy transfer buffer onto a projection buffer
-          prjBuff(:,:) = &
-            reshape( trsfBuff(trsfSize+1:trsfSize+3*maxMesh), (/3,maxMesh/) )
+          call array_copy(trsfSize+1, trsfSize+3*maxMesh, trsfBuff(:), &
+               [1, 1], [3, maxMesh], prjBuff(:,:))
+          !prjBuff(:,:) = &
+          !  reshape( trsfBuff(trsfSize+1:trsfSize+3*maxMesh), (/3,maxMesh/) )
           ! Add projection of common box to output array
           do ix = 1,3
             prjData(ix,dstComBox(1,ix,iPart):dstComBox(2,ix,iPart)) = &

=== modified file 'src/moreParallelSubs.F90'
--- src/moreParallelSubs.F90	2016-05-09 12:15:45 +0000
+++ src/moreParallelSubs.F90	2017-02-16 10:29:06 +0000
@@ -453,6 +453,7 @@
                                  a1, b1, c1, &
                                  a2, b2, &
                                  a3 )
+  use m_array, only: array_copy
   implicit none
   character(len=*),intent(in)    :: op  ! ('sum'|'prod'|'max'|'min')
   integer,optional,intent(inout)                 :: a0, b0, c0, d0, e0, f0
@@ -530,17 +531,23 @@
   end if
   if (present(a2)) then
     m = size(a2)
-    sendBuff(n+1:n+m) = reshape( a2, (/m/) )
+    call array_copy(lbound(a2), ubound(a2), a2(:,:), &
+         n+1, n+m, sendBuff(:))
+    !sendBuff(n+1:n+m) = reshape( a2, (/m/) )
     n = n + m
   end if
   if (present(b2)) then
     m = size(b2)
-    sendBuff(n+1:n+m) = reshape( b2, (/m/) )
+    call array_copy(lbound(b2), ubound(b2), b2(:,:), &
+         n+1, n+m, sendBuff(:))
+    !sendBuff(n+1:n+m) = reshape( b2, (/m/) )
     n = n + m
   end if
   if (present(a3)) then
     m = size(a3)
-    sendBuff(n+1:n+m) = reshape( a3, (/m/) )
+    call array_copy(lbound(a3), ubound(a3), a3(:,:,:), &
+         n+1, n+m, sendBuff(:))
+    !sendBuff(n+1:n+m) = reshape( a3, (/m/) )
     n = n + m
   end if
 
@@ -604,17 +611,23 @@
   end if
   if (present(a2)) then
     m = size(a2)
-    a2 = reshape( recvBuff(n+1:n+m), shape(a2) )
+    call array_copy(n+1, n+m, recvBuff(:), &
+         lbound(a2), ubound(a2), a2(:,:))
+    !a2 = reshape( recvBuff(n+1:n+m), shape(a2) )
     n = n + m
   end if
   if (present(b2)) then
     m = size(b2)
-    b2 = reshape( recvBuff(n+1:n+m), shape(b2) )
+    call array_copy(n+1, n+m, recvBuff(:), &
+         lbound(b2), ubound(b2), b2(:,:))
+    !b2 = reshape( recvBuff(n+1:n+m), shape(b2) )
     n = n + m
   end if
   if (present(a3)) then
     m = size(a3)
-    a3 = reshape( recvBuff(n+1:n+m), shape(a3) )
+    call array_copy(n+1, n+m, recvBuff(:), &
+         lbound(a3), ubound(a3), a3(:,:,:))
+    !a3 = reshape( recvBuff(n+1:n+m), shape(a3) )
     n = n + m
   end if
 
@@ -631,7 +644,7 @@
                                     a1, b1, c1, &
                                     a2, b2, &
                                     a3 )
-
+  use m_array, only: array_copy
   implicit none
   character(len=*),intent(in)    :: op  ! ('sum'|'prod'|'max'|'min')
   real(dp),optional,intent(inout)                 :: a0, b0, c0, d0, e0, f0
@@ -709,17 +722,23 @@
   end if
   if (present(a2)) then
     m = size(a2)
-    sendBuff(n+1:n+m) = reshape( a2, (/m/) )
+    call array_copy(lbound(a2), ubound(a2), a2(:,:), &
+         n+1, n+m, sendBuff(:))
+    !sendBuff(n+1:n+m) = reshape( a2, (/m/) )
     n = n + m
   end if
   if (present(b2)) then
     m = size(b2)
-    sendBuff(n+1:n+m) = reshape( b2, (/m/) )
+    call array_copy(lbound(b2), ubound(b2), b2(:,:), &
+         n+1, n+m, sendBuff(:))
+    !sendBuff(n+1:n+m) = reshape( b2, (/m/) )
     n = n + m
   end if
   if (present(a3)) then
     m = size(a3)
-    sendBuff(n+1:n+m) = reshape( a3, (/m/) )
+    call array_copy(lbound(a3), ubound(a3), a3(:,:,:), &
+         n+1, n+m, sendBuff(:))
+    !sendBuff(n+1:n+m) = reshape( a3, (/m/) )
     n = n + m
   end if
 
@@ -783,17 +802,23 @@
   end if
   if (present(a2)) then
     m = size(a2)
-    a2 = reshape( recvBuff(n+1:n+m), shape(a2) )
+    call array_copy(n+1, n+m, recvBuff(:), &
+         lbound(a2), ubound(a2), a2(:,:))
+    !a2 = reshape( recvBuff(n+1:n+m), shape(a2) )
     n = n + m
   end if
   if (present(b2)) then
     m = size(b2)
-    b2 = reshape( recvBuff(n+1:n+m), shape(b2) )
+    call array_copy(n+1, n+m, recvBuff(:), &
+         lbound(b2), ubound(b2), b2(:,:))
+    !b2 = reshape( recvBuff(n+1:n+m), shape(b2) )
     n = n + m
   end if
   if (present(a3)) then
     m = size(a3)
-    a3 = reshape( recvBuff(n+1:n+m), shape(a3) )
+    call array_copy(n+1, n+m, recvBuff(:), &
+         lbound(a3), ubound(a3), a3(:,:,:))
+    !a3 = reshape( recvBuff(n+1:n+m), shape(a3) )
     n = n + m
   end if
 

=== modified file 'src/sorting.f'
--- src/sorting.f	2016-05-09 12:15:45 +0000
+++ src/sorting.f	2017-02-16 10:29:06 +0000
@@ -103,7 +103,7 @@
       integer, intent(in) :: m, n     ! Dimensions of array x
       integer, intent(out):: indx(n)  ! Increasing order of x(1,:)
 
-      integer:: child, child2, k, nFamily, parent
+      integer:: k, nFamily, parent
       real(dp):: ageTol
 
       ! Construct the heap (family tree)

=== modified file 'src/xcmod.F90'
--- src/xcmod.F90	2017-01-25 21:34:36 +0000
+++ src/xcmod.F90	2017-02-16 10:29:06 +0000
@@ -189,7 +189,7 @@
     character(len=*), intent(inout) ::  auth
 
     call die("Libxc not compiled in. Cannot handle " //  &
-              trim(XCfunc(i)) // " " // trim(XCauth(i)))
+              trim(func) // " " // trim(auth))
   end subroutine process_libxc_spec
   
 #else

=== modified file 'version.info'
--- version.info	2016-10-31 16:23:50 +0000
+++ version.info	2017-02-16 10:29:06 +0000
@@ -1,2 +1,2 @@
-libGridXC-devel-29
+libGridXC-devel-30--reshape-2