diff --git a/src/ZSOLVER/ppm.f90 b/src/ZSOLVER/ppm.f90
index 11f508d76bfa111481bd58ae5b49fd3485a52b04..150ef055a45211c6956fbda8f9cf9bd40688b480 100644
--- a/src/ZSOLVER/ppm.f90
+++ b/src/ZSOLVER/ppm.f90
@@ -693,7 +693,12 @@ CASE('OPEN')
 !
 ! calculate dmq
 !
+#ifndef  MNH_OPENACC  
    ZDMQ = DIF2X(PSRC)
+#else
+   CALL DIF2X_DEVICE(ZDMQ,PSRC)
+#endif
+   
 !$acc kernels
 !
 ! overwrite the values on the boundary to get second order difference
@@ -1042,6 +1047,67 @@ DQ = 0.5*DQ
 !$acc end data
 
 END FUNCTION DIF2X
+!-------------------------------------------------------------------------------
+!
+!     ########################################################################
+      SUBROUTINE DIF2X_DEVICE(DQ,PQ)
+!     ########################################################################
+!!
+!!****  DIF2X - leap-frog difference operator in X direction
+!!
+!!    Calculates the difference assuming periodic BC (CYCL). 
+!!
+!!    DQ(I) = 0.5 * (PQ(I+1) - PQ(I-1))
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!    
+!!    18.3.2006.  T. Maric - original version
+!!    07/2010     J.Escobar : Correction for reproducility
+!!    04/2017     J.Escobar : initialize realistic value in all HALO pts
+!-------------------------------------------------------------------------------
+!
+!
+USE MODE_ll
+!
+IMPLICIT NONE
+! 
+!*       0.1   Declarations of dummy arguments :
+!   
+REAL, DIMENSION(:,:,:), INTENT(IN)                :: PQ
+REAL, DIMENSION(:,:,:)                            :: DQ
+!
+!*       0.2   Declarations of local variables :
+!   
+INTEGER :: IIB,IJB        ! Begining useful area in x,y directions
+INTEGER :: IIE,IJE        ! End useful area in x,y directions
+!
+!-------------------------------------------------------------------------------
+
+!$acc data present( PQ, DQ )
+!
+!*       1.0.     COMPUTE THE DOMAIN DIMENSIONS
+!                 -----------------------------
+!
+!!$CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+IIB=2 ; IIE = SIZE(PQ,1) -1
+IJB=2 ; IJE = SIZE(PQ,2) -1
+!
+!-------------------------------------------------------------------------------
+!
+!*       2.0.     COMPUTE THE DIFFERENCE
+!                 ----------------------
+!
+!$acc kernels
+DQ(IIB:IIE,:,:) = PQ(IIB+1:IIE+1,:,:) - PQ(IIB-1:IIE-1,:,:)
+DQ(IIB-1,:,:) = PQ(IIB,:,:) - PQ(IIE-1,:,:)
+DQ(IIE+1,:,:) = PQ(IIB+1,:,:) - PQ(IIE,:,:)  
+DQ = 0.5*DQ
+!$acc end kernels
+
+!$acc end data
+
+END SUBROUTINE DIF2X_DEVICE
 !
 #ifdef MNH_OPENACC
 END SUBROUTINE PPM_01_X
@@ -1515,7 +1581,11 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
 CASE('OPEN')
 !
 ! calculate dmq
+#ifndef  MNH_OPENACC   
    ZDMQ = DIF2Y(PSRC)
+#else
+   CALL DIF2Y_DEVICE(ZDMQ,PSRC)
+#endif   
 !$acc kernels
 ! overwrite the values on the boundary to get second order difference
 ! for qL and qR at the boundary
@@ -1861,6 +1931,69 @@ DQ = 0.5 * DQ
 !$acc end data
 
 END FUNCTION DIF2Y
+!-------------------------------------------------------------------------------
+!
+!     ########################################################################
+      SUBROUTINE DIF2Y_DEVICE(DQ,PQ)
+!     ########################################################################
+!!
+!!****  DIF2Y - leap-frog difference operator in Y direction
+!!
+!!    Calculates the difference assuming periodic BC (CYCL). 
+!!
+!!    DQ(J) = 0.5 * (PQ(J+1) - PQ(J-1))
+!!
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!    
+!!    18.3.2006.  T. Maric - original version, works only for periodic boundary
+!!                           conditions and on one domain
+!!    04/2017     J.Escobar : initialize realistic value in all HALO pts
+!!
+!-------------------------------------------------------------------------------
+!
+!
+USE MODE_ll
+!
+IMPLICIT NONE
+! 
+!*       0.1   Declarations of dummy arguments :
+!   
+REAL, DIMENSION(:,:,:), INTENT(IN)                :: PQ
+REAL, DIMENSION(:,:,:)                            :: DQ
+!
+!*       0.2   Declarations of local variables :
+!   
+INTEGER :: IIB,IJB        ! Begining useful area in x,y directions
+INTEGER :: IIE,IJE        ! End useful area in x,y directions
+!
+!-------------------------------------------------------------------------------
+
+!$acc data present(PQ, DQ)
+!
+!*       1.0.     COMPUTE THE DOMAIN DIMENSIONS
+!                 -----------------------------
+!
+!!$CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+IIB=2 ; IIE = SIZE(PQ,1) -1
+IJB=2 ; IJE = SIZE(PQ,2) -1
+!
+!-------------------------------------------------------------------------------
+!
+!*       2.0.     COMPUTE THE DIFFERENCE
+!                 ----------------------
+!
+!$acc kernels
+DQ(:,IJB:IJE,:) = PQ(:,IJB+1:IJE+1,:) - PQ(:,IJB-1:IJE-1,:)
+DQ(:,IJB-1,:) = PQ(:,IJB,:) - PQ(:,IJE-1,:)
+DQ(:,IJE+1,:) = PQ(:,IJB+1,:) - PQ(:,IJE,:) 
+DQ = 0.5 * DQ
+!$acc end kernels
+
+!$acc end data
+
+END SUBROUTINE DIF2Y_DEVICE
 ! #endif
 !
 #ifdef MNH_OPENACC
@@ -2077,7 +2210,11 @@ ZFNEG(:,:,:) = PSRC(:,:,:)
 !               --------------------------------
 ! 
 ! calculate dmq
+#ifndef  MNH_OPENACC
 ZDMQ = DIF2Z(PSRC)
+#else
+CALL DIF2Z_DEVICE(ZDMQ,PSRC)
+#endif
 !$acc kernels
 !
 ! monotonize the difference followinq eq. 5 in Lin94
@@ -2404,6 +2541,70 @@ DQ = 0.5 * DQ
 !$acc end data
 
 END FUNCTION DIF2Z
+!-------------------------------------------------------------------------------
+!
+!     ########################################################################
+      SUBROUTINE DIF2Z_DEVICE(DQ,PQ)
+!     ########################################################################
+!!
+!!****  DIF2Z - leap-frog difference operator in Z direction
+!!
+!!    Calculates the difference assuming periodic BC (CYCL). 
+!!
+!!    DQ(K) = 0.5 * (PQ(K+1) - PQ(K-1))
+!!
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!    
+!!    18.3.2006.  T. Maric - original version
+!!
+!-------------------------------------------------------------------------------
+!
+!
+USE MODE_ll
+USE MODD_CONF
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+! 
+!*       0.1   Declarations of dummy arguments :
+!   
+REAL, DIMENSION(:,:,:), INTENT(IN)                :: PQ
+REAL, DIMENSION(:,:,:)                            :: DQ
+!
+!*       0.2   Declarations of local variables :
+!   
+INTEGER :: IKB    ! Begining useful area in z directions
+INTEGER :: IKE    ! End useful area in z directions
+!
+!-------------------------------------------------------------------------------
+
+!$acc data present( PQ, DQ )
+!
+!*       1.0.     COMPUTE THE DOMAIN DIMENSIONS
+!                 -----------------------------
+!
+!CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+IKB = 1 + JPVEXT
+IKE = SIZE(PQ,3) - JPVEXT
+!
+!-------------------------------------------------------------------------------
+!
+!*       2.0.     COMPUTE THE DIFFERENCE
+!                 ----------------------
+!
+!$acc kernels
+DQ(:,:,IKB:IKE) = PQ(:,:,IKB+1:IKE+1) - PQ(:,:,IKB-1:IKE-1)
+DQ(:,:,IKB-1) = -DQ(:,:,IKB)
+DQ(:,:,IKE+1) = -DQ(:,:,IKE)
+DQ = 0.5 * DQ
+!$acc end kernels
+
+!$acc end data
+
+END SUBROUTINE DIF2Z_DEVICE
+
 ! #endif
 !
 #ifdef MNH_OPENACC