diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_exyz/kahan_sum.F src/kahan_sum.F
*** src_exyz/kahan_sum.F	2012-06-04 15:07:26.312132812 +0100
--- src/kahan_sum.F	2012-06-04 15:09:45.346622453 +0100
***************
*** 44,50 ****
  
    IMPLICIT NONE
    PRIVATE
!   PUBLIC :: accurate_sum
    INTEGER, PARAMETER :: sp=KIND(0.0),dp=KIND(0.0D0)
    REAL(KIND=sp),    PARAMETER :: szero=0.0_sp
    REAL(KIND=dp),    PARAMETER :: dzero=0.0_dp
--- 44,50 ----
  
    IMPLICIT NONE
    PRIVATE
!   PUBLIC :: accurate_sum, pairwise_sum
    INTEGER, PARAMETER :: sp=KIND(0.0),dp=KIND(0.0D0)
    REAL(KIND=sp),    PARAMETER :: szero=0.0_sp
    REAL(KIND=dp),    PARAMETER :: dzero=0.0_dp
***************
*** 60,65 ****
--- 60,69 ----
       kahan_sum_s6,      kahan_sum_d6,      kahan_sum_c6,      kahan_sum_z6,&
       kahan_sum_s7,      kahan_sum_d7,      kahan_sum_c7,      kahan_sum_z7
    END INTERFACE accurate_sum
+   INTERFACE pairwise_sum
+     MODULE PROCEDURE & 
+      pairwise_sum_d1
+   END INTERFACE pairwise_sum
   CONTAINS
  ! *****************************************************************************
    FUNCTION kahan_sum_s1(array,mask) RESULT(ks)
***************
*** 123,128 ****
--- 127,164 ----
      ENDIF
    END FUNCTION kahan_sum_d1
  
+   RECURSIVE FUNCTION pairwise_sum_d1(array,mask) RESULT(ps)
+     REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: array
+     LOGICAL, DIMENSION(:), INTENT(IN), &
+       OPTIONAL                               :: mask
+     REAL(KIND=dp)                            :: ps
+ 
+     INTEGER :: N, Nh
+ 
+     N = size(array)
+     Nh = size(array)/2
+ 
+     if (N == 0) then ! terminate
+       ps = 0.0_dp
+     else if (N == 1) then ! terminate
+       if (present(mask)) then
+ 	if (mask(1)) then
+ 	  ps = array(1)
+ 	else
+ 	  ps = 0.0_dp
+ 	endif
+       else
+ 	ps = array(1)
+       endif
+     else ! N >= 2
+       if (present(mask)) then
+ 	ps = pairwise_sum_d1(array(1:Nh),mask(1:Nh)) + pairwise_sum_d1(array(Nh+1:N),mask(Nh+1:N))
+       else
+ 	ps = pairwise_sum_d1(array(1:Nh)) + pairwise_sum_d1(array(Nh+1:N))
+       endif
+     endif
+   END FUNCTION pairwise_sum_d1
+ 
  ! *****************************************************************************
    FUNCTION kahan_sum_c1(array,mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:), &
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_exyz/pw_grids.F src/pw_grids.F
*** src_exyz/pw_grids.F	2012-06-04 15:06:02.795630124 +0100
--- src/pw_grids.F	2012-06-04 15:10:01.975642426 +0100
***************
*** 878,884 ****
            i2 = SIZE ( yz_mask, 2 )
            ALLOCATE ( yz_index(2,i1*i2), STAT = ierr )
            CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure)
!           CALL order_mask ( yz_mask, yz_index, error )
            DO i = 1, i1*i2
               lo(1) = yz_index(1,i)
               lo(2) = yz_index(2,i)
--- 878,884 ----
            i2 = SIZE ( yz_mask, 2 )
            ALLOCATE ( yz_index(2,i1*i2), STAT = ierr )
            CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure)
!           CALL order_mask ( yz_mask, yz_index, error , pw_grid%para%my_pos)
            DO i = 1, i1*i2
               lo(1) = yz_index(1,i)
               lo(2) = yz_index(2,i)
***************
*** 1099,1109 ****
    !------------------------------------------------------------------------------
  
  ! *****************************************************************************
!   SUBROUTINE order_mask ( yz_mask, yz_index, error )
  
      INTEGER, DIMENSION(:, :), INTENT(IN)     :: yz_mask
      INTEGER, DIMENSION(:, :), INTENT(OUT)    :: yz_index
      TYPE(cp_error_type), INTENT(inout)       :: error
  
      CHARACTER(len=*), PARAMETER :: routineN = 'order_mask', &
        routineP = moduleN//':'//routineN
--- 1099,1110 ----
    !------------------------------------------------------------------------------
  
  ! *****************************************************************************
!   SUBROUTINE order_mask ( yz_mask, yz_index, error , my_pos)
  
      INTEGER, DIMENSION(:, :), INTENT(IN)     :: yz_mask
      INTEGER, DIMENSION(:, :), INTENT(OUT)    :: yz_index
      TYPE(cp_error_type), INTENT(inout)       :: error
+ integer, intent(in) :: my_pos
  
      CHARACTER(len=*), PARAMETER :: routineN = 'order_mask', &
        routineP = moduleN//':'//routineN
***************
*** 1112,1120 ****
--- 1113,1125 ----
                                                  j1, j2
      INTEGER, ALLOCATABLE, DIMENSION(:)       :: cindex, icol, irow, rindex
      LOGICAL                                  :: failure
+ !NB load balance
+ integer :: im, ic, jc, ii, jj, iit, jjt, ii_lb, ii_ub, jj_lb, jj_ub
  
  !------------------------------------------------------------------------------
  
+ 
+ #ifdef ORIGINAL_RAY_DISTRIBUTION
      failure = .FALSE.
  
      i1 = SIZE ( yz_mask, 1 )
***************
*** 1155,1160 ****
--- 1160,1323 ----
      DEALLOCATE ( irow, rindex, icol, cindex, STAT = ierr )
      CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure)
  
+ #else
+ 
+ #ifdef RAY_SPIRAL_DISTRIBUTION_SLOW
+ 
+    icount = 1
+    ic=i1/2
+    jc=i2/2
+    ii=1
+    jj=1
+    if (ii > 0 .and. ii <= i1 .and. jj > 0 .and. jj <= i2) then
+       if (yz_mask(ii,jj) /= 0) then
+ 	 yz_index(1,icount) = ii
+ 	 yz_index(2,icount) = jj
+ 	 icount = icount + 1
+       endif
+    endif
+    do im=1, max(ic,jc)
+ 
+       if (mod(i1,2) == 0) then
+ 	 ii_lb = max(1-im,1-ic+1)
+       else
+ 	 ii_lb = max(1-im,1-ic)
+       endif
+       ii_ub = min(1+im,1+ic)
+ 
+       if (mod(i2,2) == 0) then
+ 	 jj_lb = max(1-im,1-jc+1)
+       else
+ 	 jj_lb = max(1-im,1-jc)
+       endif
+       jj_ub = min(1+im,1+jc)
+ 
+       if (im < ic .or. (im == ic .and. mod(i1,2) /= 0)) then
+ 	 ii = 1-im+i1
+ 	 do jj=jj_lb, jj_ub
+ 	    jjt = jj
+ 	    if (jjt <= 0) jjt = jjt + i2
+ 	    if (ii > 0 .and. ii <= i1 .and. jjt > 0 .and. jjt <= i2) then
+ 	       if (yz_mask(ii,jjt) /= 0) then
+ 		  yz_index(1,icount) = ii
+ 		  yz_index(2,icount) = jjt
+ 		  icount = icount + 1
+ 	       endif
+ 	    endif
+ 	 end do
+       endif
+ 
+       if (im <= ic) then
+ 	 ii = 1+im
+ 	 do jj=jj_lb, jj_ub
+ 	    jjt = jj
+ 	    if (jjt <= 0) jjt = jjt + i2
+ 	    if (ii > 0 .and. ii <= i1 .and. jjt > 0 .and. jjt <= i2) then
+ 	       if (yz_mask(ii,jjt) /= 0) then
+ 		  yz_index(1,icount) = ii
+ 		  yz_index(2,icount) = jjt
+ 		  icount = icount + 1
+ 	       endif
+ 	    endif
+ 	 end do
+       endif
+ 
+       if (im < jc .or. (im == jc .and. mod(i2,2) /= 0)) then
+ 	 jj = 1-im+i2
+ 	 do ii=1-im+1,1+im-1
+ 	    iit = ii
+ 	    if (iit <= 0) iit = iit + i1
+ 	    if (iit > 0 .and. iit <= i1 .and. jj > 0 .and. jj <= i2) then
+ 	       if (yz_mask(iit,jj) /= 0) then
+ 		  yz_index(1,icount) = iit
+ 		  yz_index(2,icount) = jj
+ 		  icount = icount + 1
+ 	       endif
+ 	    endif
+ 	 end do
+       endif
+ 
+       if (im <= jc) then
+ 	 jj = 1+im
+ 	 do ii=1-im+1,1+im-1
+ 	    iit = ii
+ 	    if (iit <= 0) iit = iit + i1
+ 	    if (iit > 0 .and. iit <= i1 .and. jj > 0 .and. jj <= i2) then
+ 	       if (yz_mask(iit,jj) /= 0) then
+ 		  yz_index(1,icount) = iit
+ 		  yz_index(2,icount) = jj
+ 		  icount = icount + 1
+ 	       endif
+ 	    endif
+ 	 end do
+       endif
+ 
+    end do
+ 
+ #else
+ 
+     i1 = SIZE ( yz_mask, 1 )
+     i2 = SIZE ( yz_mask, 2 )
+     yz_index = 0
+ 
+    icount = 1
+    ic=i1/2
+    jc=i2/2
+    ii=ic
+    jj=jc
+    if (ii > 0 .and. ii <= i1 .and. jj > 0 .and. jj <= i2) then
+       if (yz_mask(ii,jj) /= 0) then
+ 	 yz_index(1,icount) = ii
+ 	 yz_index(2,icount) = jj
+ 	 icount = icount + 1
+       endif
+    endif
+    do im=1, max(ic+1,jc+1)
+       ii = ic-im
+       do jj=jc-im,jc+im
+ 	 if (ii > 0 .and. ii <= i1 .and. jj > 0 .and. jj <= i2) then
+ 	    if (yz_mask(ii,jj) /= 0) then
+ 	       yz_index(1,icount) = ii
+ 	       yz_index(2,icount) = jj
+ 	       icount = icount + 1
+ 	    endif
+ 	 endif
+       end do
+       ii = ic+im
+       do jj=jc-im,jc+im
+ 	 if (ii > 0 .and. ii <= i1 .and. jj > 0 .and. jj <= i2) then
+ 	    if (yz_mask(ii,jj) /= 0) then
+ 	       yz_index(1,icount) = ii
+ 	       yz_index(2,icount) = jj
+ 	       icount = icount + 1
+ 	    endif
+ 	 endif
+       end do
+       jj = jc-im
+       do ii=ic-im+1,ic+im-1
+ 	 if (ii > 0 .and. ii <= i1 .and. jj > 0 .and. jj <= i2) then
+ 	    if (yz_mask(ii,jj) /= 0) then
+ 	       yz_index(1,icount) = ii
+ 	       yz_index(2,icount) = jj
+ 	       icount = icount + 1
+ 	    endif
+ 	 endif
+       end do
+       jj = jc+im
+       do ii=ic-im+1,ic+im-1
+ 	 if (ii > 0 .and. ii <= i1 .and. jj > 0 .and. jj <= i2) then
+ 	    if (yz_mask(ii,jj) /= 0) then
+ 	       yz_index(1,icount) = ii
+ 	       yz_index(2,icount) = jj
+ 	       icount = icount + 1
+ 	    endif
+ 	 endif
+       end do
+    end do
+ 
+ #endif
+ #endif
+ 
    END SUBROUTINE order_mask
  
  ! *****************************************************************************
