diff --git a/PH/add_zstar_ue_us.f90 b/PH/add_zstar_ue_us.f90
index 6b58b70f6..0cb5b453f 100644
--- a/PH/add_zstar_ue_us.f90
+++ b/PH/add_zstar_ue_us.f90
@@ -23,6 +23,9 @@ subroutine add_zstar_ue_us(imode0,npe)
   USE wavefunctions_module,    ONLY : evc
   USE io_files, ONLY: iunigk
   USE phcom
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
 
   integer, intent(in) :: imode0, npe
@@ -62,7 +65,7 @@ subroutine add_zstar_ue_us(imode0,npe)
            pdsp = (0.d0,0.d0)
            call psidspsi (ik, u (1, mode), pdsp,npw)
 #ifdef __PARA
-           call reduce(2*nbnd*nbnd,pdsp)
+           call mp_sum(pdsp, intra_pool_comm )
 #endif
            !
            ! add the term of the double summation
diff --git a/PH/addnlcc.f90 b/PH/addnlcc.f90
index bec145ce3..591f4a356 100644
--- a/PH/addnlcc.f90
+++ b/PH/addnlcc.f90
@@ -18,6 +18,9 @@ subroutine addnlcc (imode0, drhoscf, npe)
   use scf, only : rho, rho_core
   USE kinds, only : DP
   use phcom
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
 
   integer :: imode0, npe
@@ -111,7 +114,7 @@ subroutine addnlcc (imode0, drhoscf, npe)
   !
   ! collect contributions from all pools (sum over k-points)
   !
-  call reduce (18 * nat * nat, dyn1)
+  call mp_sum ( dyn1, intra_pool_comm )
 #endif
   dyn (:,:) = dyn(:,:) + dyn1(:,:) 
   deallocate (dvaux)
diff --git a/PH/cch_psi_all.f90 b/PH/cch_psi_all.f90
index a913b6e3d..192c4d627 100644
--- a/PH/cch_psi_all.f90
+++ b/PH/cch_psi_all.f90
@@ -19,6 +19,9 @@ subroutine cch_psi_all (n, h, ah, e, ik, m)
   use becmod
   USE kinds, only : DP
   use phcom
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
 
   integer :: n, m, ik
@@ -79,7 +82,7 @@ subroutine cch_psi_all (n, h, ah, e, ik, m)
        npwx, spsi, npwx, (0.d0, 0.d0) , ps, nbnd)
   ps (:,:) = ps(:,:) * alpha_pv
 #ifdef __PARA
-  call reduce (2 * nbnd * m, ps)
+  call mp_sum (ps, intra_pool_comm)
 #endif
 
   hpsi (:,:) = (0.d0, 0.d0)
diff --git a/PH/cgsolve_all.f90 b/PH/cgsolve_all.f90
index 902677376..c2fa9ea73 100644
--- a/PH/cgsolve_all.f90
+++ b/PH/cgsolve_all.f90
@@ -57,6 +57,9 @@ subroutine cgsolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, &
   !   revised (to reduce memory) 29 May 2004 by S. de Gironcoli
   !
   USE kinds, only : DP
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
   !
   !   first the I/O variables
@@ -155,7 +158,7 @@ subroutine cgsolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, &
      enddo
      kter_eff = kter_eff + DBLE (lbnd) / DBLE (nbnd)
 #ifdef __PARA
-     call reduce (lbnd, rho )
+     call mp_sum(  rho(1:lbnd) , intra_pool_comm )
 #endif
      do ibnd = nbnd, 1, -1
         if (conv(ibnd).eq.0) then
@@ -212,8 +215,8 @@ subroutine cgsolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, &
         end if
      end do
 #ifdef __PARA
-     call reduce (lbnd, a)
-     call reduce (lbnd, c)
+     call mp_sum(  a(1:lbnd), intra_pool_comm )
+     call mp_sum(  c(1:lbnd), intra_pool_comm )
 #endif
      lbnd=0
      do ibnd = 1, nbnd
diff --git a/PH/ch_psi_all.f90 b/PH/ch_psi_all.f90
index f9e048e7d..74f65048f 100644
--- a/PH/ch_psi_all.f90
+++ b/PH/ch_psi_all.f90
@@ -20,6 +20,9 @@ subroutine ch_psi_all (n, h, ah, e, ik, m)
   USE noncollin_module, ONLY : noncolin, npol
   USE kinds, only : DP
   use phcom
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
 
   integer :: n, m, ik
@@ -93,7 +96,7 @@ subroutine ch_psi_all (n, h, ah, e, ik, m)
   ENDIF
   ps (:,:) = ps(:,:) * alpha_pv
 #ifdef __PARA
-  call reduce (2 * nbnd * m, ps)
+  call mp_sum ( ps, intra_pool_comm )
 #endif
 
   hpsi (:,:) = (0.d0, 0.d0)
diff --git a/PH/compute_nldyn.f90 b/PH/compute_nldyn.f90
index 856eee323..a5d1b639e 100644
--- a/PH/compute_nldyn.f90
+++ b/PH/compute_nldyn.f90
@@ -23,6 +23,9 @@ subroutine compute_nldyn (wdyn, wgg, becq, alpq)
   USE kinds, only : DP
   USE uspp_param, ONLY: nh
   use phcom
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
 
   complex(DP) :: becq (nkb, npol, nbnd, nksq), alpq(nkb, npol, nbnd, 3, nksq), &
@@ -353,7 +356,7 @@ subroutine compute_nldyn (wdyn, wgg, becq, alpq)
      enddo
   enddo
 #ifdef __PARA
-  call reduce (2 * 3 * nat * 3 * nat, dynwrk)
+  call mp_sum ( dynwrk, intra_pool_comm )
 #endif
   do nu_i = 1, 3 * nat
      do nu_j = 1, 3 * nat
diff --git a/PH/d2ionq.f90 b/PH/d2ionq.f90
index e9cdfc18e..57dccfce8 100644
--- a/PH/d2ionq.f90
+++ b/PH/d2ionq.f90
@@ -24,6 +24,9 @@ subroutine d2ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
   USE io_global,  ONLY : stdout
   USE kinds, only : DP
   USE constants, ONLY: e2, tpi, fpi
+  USE mp_global, ONLY: intra_pool_comm
+  USE mp,        ONLY: mp_sum
+
   implicit none
   !
   !  first the dummy variables
@@ -223,7 +226,7 @@ subroutine d2ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
   enddo
 #ifdef __PARA
 100 continue
-  call reduce (18 * nat * nat, dy3)
+  call mp_sum ( dy3, intra_pool_comm )
 #endif
   !
   !   The dynamical matrix was computed in cartesian axis and now we put