From 2c0e713985e3898ccb2bec24739ca4acd2f425ec Mon Sep 17 00:00:00 2001
From: Juan ESCOBAR <juan.escobar@aero.obs-mip.fr>
Date: Fri, 26 Aug 2022 16:59:27 +0000
Subject: [PATCH] Juan 26/08/2022:ZSOLVER/*.f90 , Cray/nvidia Opt,
 mnh_expand+(LOOP)+OPENACC macro

---
 src/ZSOLVER/ppm.f90                           | 1827 +++++++++++------
 .../communication.f90                         |  206 +-
 .../datatypes.f90                             |   21 +-
 .../discretisation.f90                        |   75 +-
 .../mg_main_mnh.f90                           |    2 +-
 .../modd_mpif.f90                             |    3 +
 .../tensorproductmultigrid_Source/mode_mg.f90 |   13 +-
 .../multigrid.f90                             |   34 +-
 8 files changed, 1322 insertions(+), 859 deletions(-)
 create mode 100644 src/ZSOLVER/tensorproductmultigrid_Source/modd_mpif.f90

diff --git a/src/ZSOLVER/ppm.f90 b/src/ZSOLVER/ppm.f90
index 099e9adb2..83f0849e7 100644
--- a/src/ZSOLVER/ppm.f90
+++ b/src/ZSOLVER/ppm.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -319,14 +319,62 @@ END MODULE MODI_PPM
 !
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
+!     ########################################################################
+!!$      FUNCTION PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
+!!$               RESULT(PR)
       SUBROUTINE  PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
+!     ########################################################################
+
+  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+  USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
+
+  IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
+!
+! output source term
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+
+INTEGER :: IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG
+
+!$acc data present( PSRC, PCR, PRHO, PR )
+
+        CALL  MNH_GET_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)
+
+        CALL  PPM_01_X_D(IIU,IJU,IKU,HLBCX, KGRID, &
+                     & PSRC, PCR, PRHO, PTSTEP, PR, &
+                     & ZT3D(:,:,:,IZQL),ZT3D(:,:,:,IZQR),ZT3D(:,:,:,IZDQ),ZT3D(:,:,:,IZQ6), &
+                     & ZT3D(:,:,:,IZDMQ),ZT3D(:,:,:,IZQL0),ZT3D(:,:,:,IZQR0), ZT3D(:,:,:,IZQ60), &
+                     & ZT3D(:,:,:,IZFPOS),ZT3D(:,:,:,IZFNEG) )
+
+        CALL  MNH_REL_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)
+
+!$acc end data
+!
+CONTAINS
+!
+!     ########################################################################
+        SUBROUTINE  PPM_01_X_D(IIU,IJU,IKU,HLBCX, KGRID, &
+                     & PSRC, PCR, PRHO, PTSTEP, PR, &
+                     & ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG)
+
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
                RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_01_X - PPM_01 fully monotonic PPM advection scheme in X direction
 !!                 Colella notation
@@ -344,9 +392,6 @@ END MODULE MODI_PPM
 USE MODD_CONF
 
 USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
 use mode_mppdb
 #ifdef MNH_OPENACC
 use mode_msg
@@ -366,6 +411,9 @@ IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
 !
+#ifdef MNH_OPENACC
+INTEGER                        , INTENT(IN) :: IIU,IJU,IKU
+#endif
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
 !
 INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
@@ -379,7 +427,7 @@ REAL,                   INTENT(IN)  :: PTSTEP  ! Time step
 #ifndef MNH_OPENACC
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
 #else
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+REAL, DIMENSION(IIU,IJU,IKU), INTENT(OUT) :: PR
 #endif
 !
 !*       0.2   Declarations of local variables :
@@ -388,8 +436,8 @@ INTEGER:: IIB,IJB    ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE    ! End useful area in x,y,z directions
 !
 integer :: ji, jj, jk
-integer :: iiu, iju, iku
 #ifndef MNH_OPENACC
+integer :: iiu, iju, iku
 ! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL,ZQR
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZDQ,ZQ6
@@ -400,23 +448,32 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL0,ZQR0,ZQ60
 !
 ! advection fluxes
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
+!
+!BEG JUAN PPM_LL
+INTEGER                          :: IJS,IJN
+!END JUAN PPM_LL
 #else
 INTEGER                          :: I,J,K 
-! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZQL,ZQR
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDQ,ZQ6
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDMQ
 !
-! extra variables for the initial guess of parabolae parameters
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZQL0,ZQR0,ZQ60
+!!$!
+!!$! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
+REAL , DIMENSION(IIU,IJU,IKU) :: &
+  ZQL,ZQR, ZDQ,ZQ6, ZDMQ &
+!!$!
+!!$! extra variables for the initial guess of parabolae parameters
+  , ZQL0,ZQR0,ZQ60 &
+!!$!
+!!$! advection fluxes
+  , ZFPOS, ZFNEG
 !
-! advection fluxes
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFPOS, ZFNEG
-#endif
 INTEGER                :: IJS,IJN
+#endif
 LOGICAL                :: GWEST , GEAST
 !-------------------------------------------------------------------------------
 
+!$acc data present( PSRC, PCR, PRHO, PR , &
+!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )
+
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK(PCR, "PPM_01_X beg:PCR")
@@ -424,30 +481,6 @@ IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PSRC,"PPM_01_X beg:PSRC")
 END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZQL,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQR,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZDQ,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQ6,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZDMQ,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQL0,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQR0,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQ60,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFPOS, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFNEG, IIU, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PR , &
-!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
@@ -464,6 +497,10 @@ GEAST = LEAST_ll()
 !*              initialise & update halo & halo2 for PSRC
 !
 #ifndef MNH_OPENACC
+iiu = size( PSRC, 1 )
+iju = size( PSRC, 2 )
+iku = size( PSRC, 3 )
+
 CALL GET_HALO(PSRC, HNAME='PSRC')
 !
 PR   (:,:,:) = PSRC(:,:,:)
@@ -478,11 +515,13 @@ ZQ60 (:,:,:) = PSRC(:,:,:)
 ZFPOS(:,:,:) = PSRC(:,:,:)
 ZFNEG(:,:,:) = PSRC(:,:,:)
 #else
-CALL GET_HALO_D(PSRC,HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::PSRC')
+CALL GET_HALO_D(PSRC,HDIR="01_X", HNAME='PSRC')
 !
 !$acc kernels 
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = 1, iiu
         PR   (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQL  (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQR  (ji, jj, jk ) = PSRC(ji, jj, jk )
@@ -492,10 +531,19 @@ DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
         ZQL0 (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQR0 (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQ60 (ji, jj, jk ) = PSRC(ji, jj, jk )
-     END DO
+    end do
+  end do
+end do
 !
+#if 0
+ZFPOS(:,1:IJS,:)=PSRC(:,1:IJS,:)
+ZFNEG(:,1:IJS,:)=PSRC(:,1:IJS,:)
+ZFPOS(:,IJN:,:)=PSRC(:,IJN:,:)
+ZFNEG(:,IJN:,:)=PSRC(:,IJN:,:)
+#else
 ZFPOS(:,:,:) = PSRC(:,:,:)
 ZFNEG(:,:,:) = PSRC(:,:,:)
+#endif
 !$acc end kernels
 #endif
 !
@@ -512,7 +560,7 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCX(1) == HLBCX(2)
 #endif
 !
 ! calculate dmq
-  CALL DIF2X( PSRC, ZDMQ )
+   ZDMQ = DIF2X(PSRC)
 !
 ! monotonize the difference followinq eq. 5 in Lin94
 !
@@ -644,7 +692,14 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCX(1) == HLBCX(2)
 CASE('OPEN')
 !
 ! calculate dmq
-  CALL DIF2X( PSRC, ZDMQ )
+!
+#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
 ! for qL and qR at the boundary
@@ -652,26 +707,23 @@ CASE('OPEN')
 !  WEST BOUND
 !
   IF (GWEST) THEN
-!$acc kernels
    ZDMQ(IIB-1,IJS:IJN,:) = -ZDMQ(IIB,IJS:IJN,:)
-!$acc end kernels
   ENDIF
 !
 !  EAST BOUND
 !
   IF (GEAST) THEN 
-!$acc kernels
    ZDMQ(IIE+1,IJS:IJN,:) = -ZDMQ(IIE,IJS:IJN,:)
-!$acc end kernels
   ENDIF
 !
 ! monotonize the difference followinq eq. 5 in Lin94
 !
 !  ZDMQ(i) = Fct[ ZDMQ(i),PSRC(i),PSRC(i-1),PSRC(i+1) ]
 !
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-  DO CONCURRENT (ji = iib:iie ,  jj = ijs:ijn , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = iib, iie
         ZDMQ(ji, jj, jk ) = SIGN(                                                                  &
             MIN( ABS(ZDMQ(ji, jj, jk )),                                                           &
                  2.0 * ( PSRC(ji, jj, jk )                                                         &
@@ -679,7 +731,9 @@ CASE('OPEN')
                  2.0 * (-PSRC(ji, jj, jk )                                                         &
                         + MAX(PSRC(ji - 1, jj, jk ),PSRC(ji, jj, jk ),PSRC(ji + 1, jj, jk )) )  ), &
             ZDMQ(ji, jj, jk ) )
-  END DO
+      end do
+    end do
+  end do
 !
 !  WEST BOUND
 !
@@ -704,7 +758,7 @@ CASE('OPEN')
    CALL  GET_HALO(ZDMQ, HNAME='ZDMQ')
 #else
 !$acc end kernels
-   CALL  GET_HALO_D(ZDMQ, HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::ZDMQ')
+   CALL  GET_HALO_D(ZDMQ, HDIR="01_X", HNAME='ZDMQ')
 #endif
 !$acc kernels
 !
@@ -712,42 +766,43 @@ CASE('OPEN')
 !
 !  ZQL0(i) = Fct[ PSRC(i),PSRC(i-1),ZDMQ(i),ZDMQ(i-1) ]
 !
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT (ji = iib:iie + 1 , jj = ijs:ijn ,  jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = iib, iie + 1
         ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji-1, jj, jk ) ) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji-1, jj, jk ) ) / 6.0
-   END DO
+      end do
+    end do
+  end do
 !
 #ifndef MNH_OPENACC
    CALL  GET_HALO(ZQL0, HNAME='ZQL0')
 #else
 !$acc end kernels
-   CALL  GET_HALO_D(ZQL0,HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
+   CALL  GET_HALO_D(ZQL0,HDIR="01_X", HNAME='ZQL0')
+!$acc kernels
 #endif
 !  
 !  WEST BOUND
 !
   IF (GWEST) THEN
-!$acc kernels
    ZQL0(IIB-1,IJS:IJN,:) = ZQL0(IIB,IJS:IJN,:)
-!$acc end kernels
   ENDIF
 !
-!$acc kernels present( ZQL0, ZQR0 )
    ZQR0(IIB-1:IIE,IJS:IJN,:) = ZQL0(IIB:IIE+1,IJS:IJN,:)
-!$acc end kernels
 !
 #ifndef MNH_OPENACC
    CALL  GET_HALO(ZQR0, HNAME='ZQR0')
 #else
-   CALL  GET_HALO_D(ZQR0, HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::ZQR0')
+!$acc end kernels
+   CALL  GET_HALO_D(ZQR0, HDIR="01_X", HNAME='ZQR0')
+!$acc kernels
 #endif
 !
 !  EAST BOUND
 !
   IF (GEAST) THEN
-!$acc kernels present( ZQR0 )
    ZQR0(IIE+1,IJS:IJN,:) = ZQR0(IIE,IJS:IJN,:)
-!$acc end kernels
   ENDIF
 #ifndef MNH_OPENACC
 !
@@ -782,9 +837,11 @@ CASE('OPEN')
 !
    ZDQ(:,IJS:IJN,:) = ZQR(:,IJS:IJN,:) - ZQL(:,IJS:IJN,:)
 #else
-!$acc kernels present( ZQL0, ZQR0 )
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT (I=1:IIU , J = IJS:IJN , K=1:IKU )
+!$acc loop independent collapse(3)
+DO K=1,IKU 
+ DO J = IJS,IJN
+    ! acc loop vector(24)
+   DO I=1,IIU
 !
 ! determine initial coefficients of the parabolae
 !
@@ -816,8 +873,7 @@ DO CONCURRENT (I=1:IIU , J = IJS:IJN , K=1:IKU )
 ! recalculate coefficients of the parabolae
 !
    ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
-END DO
-!$acc end kernels
+ENDDO ; ENDDO ; ENDDO
 #endif
 !
 ! and finally calculate fluxes for the advection
@@ -828,20 +884,24 @@ END DO
 !!$   ZFPOS(IIB+1:IIE+1,:,:) = ZQR(IIB:IIE,:,:) - 0.5*PCR(IIB+1:IIE+1,:,:) * &            
 !!$        (ZDQ(IIB:IIE,:,:) - (1.0 - 2.0*PCR(IIB+1:IIE+1,:,:)/3.0)          &
 !!$        * ZQ6(IIB:IIE,:,:))
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = iib:iie + 1 , jj = ijs:ijn ,  jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = iib, iie + 1
         ZFPOS(ji, jj, jk ) = ZQR(ji - 1, jj, jk ) - 0.5 * PCR(ji, jj, jk )                  &
                             * ( ZDQ(ji - 1, jj, jk) - (1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0 ) &
                             * ZQ6(ji - 1, jj, jk) )
-END DO
+      end do
+    end do
+  end do
 !
 !
 #ifndef MNH_OPENACC
    CALL GET_HALO(ZFPOS, HNAME='ZFPOS')
 #else
 !$acc end kernels
-   CALL GET_HALO_D(ZFPOS, HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::ZFPOS')
+   CALL GET_HALO_D(ZFPOS, HDIR="01_X", HNAME='ZFPOS')
+!$acc kernels
 #endif
 !
 !
@@ -850,40 +910,40 @@ END DO
 ! advection flux at open boundary when u(IIB) > 0
 ! 
   IF (GWEST) THEN
-!$acc kernels present_cr(ZFPOS,ZQR)
    ZFPOS(IIB,IJS:IJN,:) = (PSRC(IIB-1,IJS:IJN,:) - ZQR(IIB-1,IJS:IJN,:))*PCR(IIB,IJS:IJN,:) + &
                     ZQR(IIB-1,IJS:IJN,:)
 ! PPOSX(IIB-1,:,:) is not important for the calc of advection so 
 ! we set it to 0
 !!$   ZFPOS(IIB-1,:,:) = 0.0
-!$acc end kernels
    ENDIF
 !
 !!$   ZFNEG(IIB-1:IIE,:,:) = ZQL(IIB-1:IIE,:,:) - 0.5*PCR(IIB-1:IIE,:,:) *  &            
 !!$        (ZDQ(IIB-1:IIE,:,:) + (1.0 + 2.0*PCR(IIB-1:IIE,:,:)/3.0)         &
 !!$        * ZQ6(IIB-1:IIE,:,:))
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = ijs:ijn , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = 1, iiu
    ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5*PCR(ji, jj, jk ) *      &
         ( ZDQ(ji, jj, jk ) + (1.0 + 2.0*PCR(ji, jj, jk )/3.0) * ZQ6(ji, jj, jk ) )
-END DO
+      end do
+    end do
+  end do
 !
 #ifndef MNH_OPENACC
    CALL GET_HALO(ZFNEG, HNAME='ZFNEG')
 #else
 !$acc end kernels
-   CALL GET_HALO_D(ZFNEG, HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::ZFNEG')
+   CALL GET_HALO_D(ZFNEG, HDIR="01_X", HNAME='ZFNEG')
+!$acc kernels
 #endif
 !
 !  EAST BOUND
 !
 ! advection flux at open boundary when u(IIE+1) < 0
   IF (GEAST) THEN
-!$acc kernels present_cr(ZFNEG,ZQR)
    ZFNEG(IIE+1,IJS:IJN,:) = (ZQR(IIE,IJS:IJN,:)-PSRC(IIE+1,IJS:IJN,:))*PCR(IIE+1,IJS:IJN,:) + &
                       ZQR(IIE,IJS:IJN,:)
-!$acc end kernels
   ENDIF
 !
 ! advect the actual field in X direction by U*dt
@@ -894,8 +954,9 @@ END DO
    CALL GET_HALO(PR, HNAME='PR')
 #else
    !mxm(ZQL,PRHO)
+!$acc end kernels
    CALL MXM_DEVICE(PRHO,ZQL)
-!$acc kernels present_cr(ZFPOS,ZFNEG,ZQR)
+!$acc kernels
    where ( PCR(:,:,:) > 0. )
      ZQR(:,:,:) = PCR(:,:,:) * ZQL(:,:,:) * ZFPOS(:,:,:)
    elsewhere
@@ -904,7 +965,7 @@ END DO
     !dxf(PR,ZQR)
 !$acc end kernels
    CALL DXF_DEVICE(ZQR,PR)
-   CALL GET_HALO_D(PR, HDIR="01_X", HNAME='UPDATE_HALO_ll::GET_HALO::PR')
+   CALL GET_HALO_D(PR, HDIR="01_X", HNAME='PR')
 #endif
 !
 END SELECT
@@ -918,19 +979,19 @@ END IF
 
 !$acc end data
 
-#ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-
+#ifndef MNH_OPENACC
 CONTAINS
+#else
+END SUBROUTINE PPM_01_X_D
+#endif
 !
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
 !     ########################################################################
-      SUBROUTINE DIF2X( PQ, DQ )
+      FUNCTION DIF2X(PQ) RESULT(DQ)
 !     ########################################################################
+!!
 !!****  DIF2X - leap-frog difference operator in X direction
 !!
 !!    Calculates the difference assuming periodic BC (CYCL). 
@@ -946,7 +1007,7 @@ CONTAINS
 !-------------------------------------------------------------------------------
 !
 !
-! USE MODE_ll
+USE MODE_ll
 !
 IMPLICIT NONE
 ! 
@@ -985,8 +1046,69 @@ DQ = 0.5*DQ
 
 !$acc end data
 
-END SUBROUTINE DIF2X
+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
 #else
@@ -997,14 +1119,62 @@ END FUNCTION PPM_01_X
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
+!     ########################################################################
+!!$      FUNCTION PPM_01_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP) &
+!!$               RESULT(PR)
       SUBROUTINE PPM_01_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
+!     ########################################################################
+
+        USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+        USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
+
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+REAL,                   INTENT(IN)  :: PTSTEP   ! Time step 
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT):: PSRC ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR   &  ! Courant number
+                                     , PRHO     ! density
+!
+! output source term
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+
+INTEGER :: IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG
+
+!$acc data present( PSRC, PCR, PRHO, PR )
+
+        CALL  MNH_GET_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)
+
+        CALL  PPM_01_Y_D(IIU,IJU,IKU,HLBCY, KGRID, &
+                     & PSRC, PCR, PRHO, PTSTEP, PR, &
+                     & ZT3D(:,:,:,IZQL),ZT3D(:,:,:,IZQR),ZT3D(:,:,:,IZDQ),ZT3D(:,:,:,IZQ6), &
+                     & ZT3D(:,:,:,IZDMQ),ZT3D(:,:,:,IZQL0),ZT3D(:,:,:,IZQR0), ZT3D(:,:,:,IZQ60), &
+                     & ZT3D(:,:,:,IZFPOS),ZT3D(:,:,:,IZFNEG) )
+
+        CALL  MNH_REL_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)
+
+!$acc end data
+
+CONTAINS
+!
+!     ########################################################################
+        SUBROUTINE  PPM_01_Y_D(IIU,IJU,IKU,HLBCY, KGRID, &
+                     & PSRC, PCR, PRHO, PTSTEP, PR, &
+                     & ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG)
+
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_01_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP) &
                RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_01_Y - PPM_01 fully monotonic PPM advection scheme in Y direction
 !!                 Colella notation
@@ -1023,12 +1193,9 @@ USE MODD_CONF
 
 USE MODE_ll
 #ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-use mode_mppdb
-#ifdef MNH_OPENACC
 use mode_msg
 #endif
+use mode_mppdb
 
 #ifdef MNH_BITREP
 USE MODI_BITREP
@@ -1044,6 +1211,9 @@ IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
 !
+#ifdef MNH_OPENACC
+integer, intent(in) :: iiu, iju, iku
+#endif
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
 !
 INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
@@ -1058,7 +1228,7 @@ REAL,                   INTENT(IN)  :: PTSTEP  ! Time step
 #ifndef MNH_OPENACC
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
 #else
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+REAL, DIMENSION(IIU,IJU,IKU), INTENT(OUT) :: PR
 #endif
 !
 !*       0.2   Declarations of local variables :
@@ -1066,11 +1236,11 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
 INTEGER:: IIB,IJB    ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE    ! End useful area in x,y,z directions
 !
-integer :: iiu, iju, iku
 INTEGER                          :: IIW,IIA
 !
 LOGICAL                          :: GSOUTH , GNORTH
 #ifndef MNH_OPENACC
+integer ::  iiu, iju, iku
 !
 ! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL,ZQR
@@ -1085,24 +1255,27 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 !
 #else
 ! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZQL, ZQR
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDQ, ZQ6
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDMQ
-!
+REAL, DIMENSION(IIU,IJU,IKU) :: &
+     ZQL,ZQR , ZDQ,ZQ6 , ZDMQ &
 ! extra variables for the initial guess of parabolae parameters
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZQL0,ZQR0,ZQ60
-!
+   , ZQL0,ZQR0,ZQ60 &
 ! advection fluxes
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFPOS, ZFNEG
+   , ZFPOS, ZFNEG
+
 !
+!JUAN ACC
 INTEGER                          :: I,J,K
 !
 INTEGER                          :: IKB,IKE
 INTEGER                          :: IJN,IJS
+!JUAN ACC
 #endif
 integer :: ji, jj, jk
 !-------------------------------------------------------------------------------
 
+!$acc data present( PSRC, PCR, PRHO, PR, &
+!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )
+
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK(PCR, "PPM_01_Y beg:PCR")
@@ -1110,30 +1283,6 @@ IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PSRC,"PPM_01_Y beg:PSRC")
 END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZQL,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQR,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZDQ,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQ6,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZDMQ,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQL0,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQR0,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQ60,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFPOS, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFNEG, IIU, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PR, &
-!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
@@ -1146,6 +1295,10 @@ GSOUTH=LSOUTH_ll()
 GNORTH=LNORTH_ll()
 !
 #ifndef MNH_OPENACC
+iiu = size( PSRC, 1 )
+iju = size( PSRC, 2 )
+iku = size( PSRC, 3 )
+
 CALL GET_HALO(PSRC, HNAME='PSRC')
 #else
 IJS=1
@@ -1160,15 +1313,17 @@ IKE=IKU
 !IIB=2
 !IIE=IIU-1
 !
-CALL GET_HALO_D(PSRC, HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::PSRC')
+CALL GET_HALO_D(PSRC, HDIR="01_Y", HNAME='PSRC')
 #endif
 !
 !-------------------------------------------------------------------------------
 !
 !
 !$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = 1, iiu
         PR   (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQL  (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQR  (ji, jj, jk ) = PSRC(ji, jj, jk )
@@ -1178,20 +1333,34 @@ DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
         ZQL0 (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQR0 (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQ60 (ji, jj, jk ) = PSRC(ji, jj, jk )
-END DO
-ZFPOS(:,:,:) = PSRC(:,:,:)
-ZFNEG(:,:,:) = PSRC(:,:,:)
-!$acc end kernels
-!
-SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
-!
-!*       2.1    CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
-!               ---------------------------------------------
-!
-CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
-!
-! calculate dmq
-  CALL DIF2Y( PSRC, ZDMQ )
+    end do
+  end do
+end do
+#ifndef MNH_OPENACC
+ZFPOS=PSRC
+ZFNEG=PSRC
+#else
+#if 0
+ZFPOS(1:IIW,:,IKB:IKE)=PSRC(1:IIW,:,IKB:IKE)
+ZFNEG(1:IIW,:,IKB:IKE)=PSRC(1:IIW,:,IKB:IKE)
+ZFPOS(IIA:,:,IKB:IKE)=PSRC(IIA:,:,IKB:IKE)
+ZFNEG(IIA:,:,IKB:IKE)=PSRC(IIA:,:,IKB:IKE)
+#else
+ZFPOS(:,:,:) = PSRC(:,:,:)
+ZFNEG(:,:,:) = PSRC(:,:,:)
+#endif
+#endif
+!$acc end kernels
+!
+SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
+!
+!*       2.1    CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
+!               ---------------------------------------------
+!
+CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
+!
+! calculate dmq
+   ZDMQ = DIF2Y(PSRC)
 !
 ! monotonize the difference followinq eq. 5 in Lin94
 !BEG JUAN PPM_LL01
@@ -1227,7 +1396,7 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
    CALL GET_HALO(ZDMQ, HNAME='ZDMQ')
 #else
 !$acc end kernels
-   CALL GET_HALO_D(ZDMQ,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZDMQ')
+   CALL GET_HALO_D(ZDMQ,HDIR="01_Y", HNAME='ZDMQ')
 !$acc kernels
 #endif
 !
@@ -1235,28 +1404,20 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
 !
 !    ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
 !         (ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
-#if 1
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( ji = iiw:iia ,  jj = ijb:ije + 1 , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije + 1
+      do ji = iiw, iia
         ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj-1, jk )) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj-1, jk ) ) / 6.0
-   END DO
-#else
-  DO CONCURRENT ( ji = iiw : iia, jj = ijb : ije + 1, jk = 1 : iku )
-    ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj-1, jk )) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj-1, jk ) ) / 6.0
-  END DO
-#endif
+      end do
+    end do
+  end do
 !
 #ifndef MNH_OPENACC
-CALL MPPDB_CHECK(PSRC,"PPM_01_Y: PSRC")
-CALL MPPDB_CHECK(ZDMQ,"PPM_01_Y: ZDMQ")
-CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
    CALL GET_HALO(ZQL0, HNAME='ZQL0')
 #else
 !$acc end kernels
-CALL MPPDB_CHECK(PSRC,"PPM_01_Y: PSRC")
-CALL MPPDB_CHECK(ZDMQ,"PPM_01_Y: ZDMQ")
-CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
-  CALL GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
+  CALL GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='ZQL0')
 !$acc kernels
 #endif
 !
@@ -1270,7 +1431,7 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
    CALL GET_HALO(ZQR0, HNAME='ZQR0')
 #else
 !$acc end kernels
-  CALL GET_HALO_D(ZQR0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQR0')
+  CALL GET_HALO_D(ZQR0,HDIR="01_Y", HNAME='ZQR0')
 !$acc kernels
 #endif
 !
@@ -1310,8 +1471,10 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
 !
    ZDQ(IIW:IIA,:,:) = ZQR(IIW:IIA,:,:) - ZQL(IIW:IIA,:,:)
 #else
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT( I=1:IIU , J=IJS:IJN , K=IKB:IKE )
+!$acc loop independent collapse(3)
+   DO K=IKB,IKE
+      DO J=IJS,IJN
+         DO I=1,IIU
             !
             ! determine initial coefficients of the parabolae
             !
@@ -1344,9 +1507,9 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
             !
             ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
             !
+         END DO
+      END DO
    END DO
-!$acc end kernels
-!$acc kernels
 #endif
 !
 ! and finally calculate fluxes for the advection
@@ -1361,7 +1524,7 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
   CALL GET_HALO(ZFPOS, HNAME='ZFPOS')
 #else
 !$acc end kernels
-  CALL GET_HALO_D(ZFPOS,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZFPOS')
+  CALL GET_HALO_D(ZFPOS,HDIR="01_Y", HNAME='ZFPOS')
 !$acc kernels
 #endif
 !
@@ -1378,7 +1541,7 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
   CALL GET_HALO(ZFNEG, HNAME='ZFNEG')
 #else
 !$acc end kernels
-  CALL GET_HALO_D(ZFNEG,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZFNEG')
+  CALL GET_HALO_D(ZFNEG,HDIR="01_Y", HNAME='ZFNEG')
 #endif
 !
 ! advect the actual field in Y direction by V*dt
@@ -1389,13 +1552,17 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
 #else
    CALL MYM_DEVICE(PRHO,ZQL)
 !$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT( I=1:IIU , J=IJS:IJN ,  K=IKB:IKE )
+!$acc loop independent collapse(3)
+   DO K=IKB,IKE
+      DO J=IJS,IJN
+         DO I=1,IIU
           if ( PCR(I,J,K) > 0. ) then
             ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) * ZFPOS(I,J,K)
           else
             ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) * ZFNEG(I,J,K)
           end if
+         END DO
+      END DO
    END DO
 !$acc end kernels
    CALL DYF_DEVICE(ZQR,PR)
@@ -1403,7 +1570,7 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
 #ifndef MNH_OPENACC
   CALL GET_HALO(PR, HNAME='PR')
 #else
-  CALL GET_HALO_D(PR,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::PR')
+  CALL GET_HALO_D(PR,HDIR="01_Y", HNAME='PR')
 #endif
 !
 !*       2.2    NON-CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
@@ -1412,36 +1579,40 @@ CALL MPPDB_CHECK(ZQL0,"PPM_01_Y: ZQL0")
 CASE('OPEN')
 !
 ! calculate dmq
-  CALL DIF2Y( PSRC, ZDMQ )
+#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
 !
 !  SOUTH BOUND
 !
    IF (GSOUTH) THEN
-!$acc kernels
     ZDMQ(IIW:IIA,IJB-1,:) = -ZDMQ(IIW:IIA,IJB,:)
-!$acc end kernels
    ENDIF
 !
 !  NORTH BOUND
 !
    IF (GNORTH) THEN
-!$acc kernels
     ZDMQ(IIW:IIA,IJE+1,:) = -ZDMQ(IIW:IIA,IJE,:)
-!$acc end kernels
    ENDIF
 !
 ! monotonize the difference followinq eq. 5 in Lin94
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( ji = iiw:iia ,  jj = ijb:ije , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije
+      do ji = iiw, iia
         ZDMQ(ji, jj, jk ) = SIGN( &
         MIN( ABS( ZDMQ(ji, jj, jk ) ),                                                                             &
              2.0 * (   PSRC(ji, jj, jk ) - MIN( PSRC(ji, jj-1, jk ), PSRC(ji, jj, jk ), PSRC(ji, jj+1, jk ) ) ),   &
              2.0 * ( - PSRC(ji, jj, jk ) + MAX( PSRC(ji, jj-1, jk ), PSRC(ji, jj, jk ), PSRC(ji, jj+1, jk ) ) ) ), &
         ZDMQ(ji, jj, jk ) )
-   END DO
+    end do
+  end do
+end do
 !!$   ZDMQ(:,IJB-1,:) = &
 !!$        SIGN( (MIN( ABS(ZDMQ(:,IJB-1,:)), 2.0*(PSRC(:,IJB-1,:) - &
 !!$        MIN(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:))),   &
@@ -1459,45 +1630,48 @@ CASE('OPEN')
    CALL  GET_HALO(ZDMQ, HNAME='ZDMQ')
 #else
 !$acc end kernels
-   CALL  GET_HALO_D(ZDMQ,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZDMQ')
+   CALL  GET_HALO_D(ZDMQ,HDIR="01_Y", HNAME='ZDMQ')
 !$acc kernels
 #endif
 !
 ! calculate qL and qR with the modified dmq
 !
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( ji = iiw:iia ,  jj = ijb:ije + 1 , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije + 1
+      do ji = iiw, iia
         ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj-1, jk )) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj-1, jk ) ) / 6.0
-   END DO
+      end do
+    end do
+  end do
 !
 #ifndef MNH_OPENACC
    CALL  GET_HALO(ZQL0, HNAME='ZQL0')
 #else
 !$acc end kernels
-CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
+CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='ZQL0')
+!$acc kernels
 #endif
 !
 !  SOUTH BOUND
 !
    IF (GSOUTH) THEN
-!$acc kernels
     ZQL0(IIW:IIA,IJB-1,:) = ZQL0(IIW:IIA,IJB,:)
-!$acc end kernels
    ENDIF
 !
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( ji = iiw:iia ,  jj = ijb - 1:ije ,  jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb - 1, ije
+      do ji = iiw, iia
         ZQR0(ji, jj, jk ) = ZQL0(ji, jj+1, jk )
-   END DO
-!$acc end kernels
+      end do
+    end do
+  end do
 !
 !  NORTH BOUND
 !
    IF (GNORTH) THEN
-!$acc kernels
     ZQR0(IIW:IIA,IJE+1,:) = ZQR0(IIW:IIA,IJE,:)
-!$acc end kernels
    ENDIF
 #ifndef MNH_OPENACC
 !
@@ -1533,9 +1707,10 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
    ZDQ(IIW:IIA,:,:) = ZQR(IIW:IIA,:,:) - ZQL(IIW:IIA,:,:)
 !
 #else
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( I=1:IIU , J=IJS:IJN , K=IKB:IKE )
+!$acc loop independent collapse(3)
+   DO K=IKB,IKE
+      DO J=IJS,IJN
+         DO I=1,IIU
             !
             ! determine initial coefficients of the parabolae
             !
@@ -1568,8 +1743,9 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
             !
             ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
             !
+         END DO
+      END DO
    END DO
-!$acc end kernels
 #endif
 !
 ! and finally calculate fluxes for the advection
@@ -1577,18 +1753,22 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
 !!$   ZFPOS(:,IJB+1:IJE+1,:) = ZQR(:,IJB:IJE,:) - 0.5*PCR(:,IJB+1:IJE+1,:) * &            
 !!$        (ZDQ(:,IJB:IJE,:) - (1.0 - 2.0*PCR(:,IJB+1:IJE+1,:)/3.0)        &
 !!$        * ZQ6(:,IJB:IJE,:))
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT (  ji = iiw:iia , jj = ijb:ije + 1 , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije + 1
+      do ji = iiw, iia
         ZFPOS(ji, jj, jk ) = ZQR(ji, jj-1, jk ) - 0.5 * PCR(ji, jj, jk )   &
                                       * ( ZDQ(ji, jj-1, jk ) - ( 1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0 ) * ZQ6(ji, jj-1, jk ) )
-   END DO
+      end do
+    end do
+  end do
 !
 #ifndef MNH_OPENACC
   CALL GET_HALO(ZFPOS, HNAME='ZFPOS')
 #else
 !$acc end kernels
-  CALL GET_HALO_D(ZFPOS,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZFPOS')
+  CALL GET_HALO_D(ZFPOS,HDIR="01_Y", HNAME='ZFPOS')
+!$acc kernels
 #endif
 !
 !
@@ -1597,10 +1777,8 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
 !  SOUTH BOUND
 !
    IF (GSOUTH) THEN
-!$acc kernels present_cr(ZFPOS,ZQR)
     ZFPOS(IIW:IIA,IJB,:) = (PSRC(IIW:IIA,IJB-1,:) - ZQR(IIW:IIA,IJB-1,:))*PCR(IIW:IIA,IJB,:) + &
                       ZQR(IIW:IIA,IJB-1,:)
-!$acc end kernels
    ENDIF
 !
 ! PPOSX(:,IJB-1,:) is not important for the calc of advection so 
@@ -1610,18 +1788,22 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
 !!$   ZFNEG(:,IJB-1:IJE,:) = ZQL(:,IJB-1:IJE,:) - 0.5*PCR(:,IJB-1:IJE,:) * &            
 !!$        ( ZDQ(:,IJB-1:IJE,:) + (1.0 + 2.0*PCR(:,IJB-1:IJE,:)/3.0) * &
 !!$        ZQ6(:,IJB-1:IJE,:) )
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( ji = iiw:iia , jj = 1:iju , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = iiw, iia
         ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5 * PCR(ji, jj, jk )      &
                                      * ( ZDQ(ji, jj, jk ) + ( 1.0 + 2.0 * PCR(ji, jj, jk ) / 3.0 ) * ZQ6(ji, jj, jk ) )
-   END DO
+      end do
+    end do
+  end do
 !
 #ifndef MNH_OPENACC
   CALL GET_HALO(ZFNEG, HNAME='ZFNEG')
 #else
 !$acc end kernels
-  CALL GET_HALO_D(ZFNEG,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZFNEG')
+  CALL GET_HALO_D(ZFNEG,HDIR="01_Y", HNAME='ZFNEG')
+!$acc kernels
 #endif
 !
 ! advection flux at open boundary when u(IJE+1) < 0
@@ -1629,10 +1811,8 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
 !  NORTH BOUND
 !
    IF (GNORTH) THEN
-!$acc kernels present_cr(ZFNEG,ZQR)
     ZFNEG(IIW:IIA,IJE+1,:) = (ZQR(IIW:IIA,IJE,:)-PSRC(IIW:IIA,IJE+1,:))*PCR(IIW:IIA,IJE+1,:) + &
                         ZQR(IIW:IIA,IJE,:)
-!$acc end kernels
    ENDIF
 #ifndef MNH_OPENACC
 !
@@ -1642,15 +1822,20 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 !
 #else
+!$acc end kernels
    CALL MYM_DEVICE(PRHO,ZQL)
 !$acc kernels
-!$acc_nv loop independent collapse(3)
-   DO CONCURRENT ( I=1:IIU , J=IJS:IJN , K=IKB:IKE )
+!$acc loop independent collapse(3)
+   DO K=IKB,IKE
+      DO J=IJS,IJN
+         DO I=1,IIU
           if ( PCR(I,J,K) > 0. ) then
             ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) * ZFPOS(I,J,K)
           else
             ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) * ZFNEG(I,J,K)
           end if
+         END DO
+      END DO
    END DO
 !$acc end kernels
    CALL DYF_DEVICE(ZQR,PR)
@@ -1659,7 +1844,7 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::ZQL0')
 #ifndef MNH_OPENACC
   CALL GET_HALO(PR, HNAME='PR')
 #else
-  CALL GET_HALO_D(PR,HDIR="01_Y", HNAME='UPDATE_HALO_ll::GET_HALO::PR')
+  CALL GET_HALO_D(PR,HDIR="01_Y", HNAME='PR')
 #endif
 !
 !
@@ -1674,18 +1859,80 @@ END IF
 
 !$acc end data
 
-#ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
+#ifndef MNH_OPENACC
+CONTAINS
+#else
+END SUBROUTINE PPM_01_Y_D
 #endif
+!
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+!
+!     ########################################################################
+      FUNCTION DIF2Y(PQ) RESULT(DQ)
+!     ########################################################################
+!!
+!!****  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(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: 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
+!
+!-------------------------------------------------------------------------------
 
-CONTAINS
+!$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 FUNCTION DIF2Y
 !-------------------------------------------------------------------------------
 !
 !     ########################################################################
-      SUBROUTINE DIF2Y( PQ, DQ )
+      SUBROUTINE DIF2Y_DEVICE(DQ,PQ)
 !     ########################################################################
 !!
 !!****  DIF2Y - leap-frog difference operator in Y direction
@@ -1711,8 +1958,8 @@ IMPLICIT NONE
 ! 
 !*       0.1   Declarations of dummy arguments :
 !   
-REAL, DIMENSION(:,:,:),                            INTENT(IN)  :: PQ
-REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)), INTENT(OUT) :: DQ
+REAL, DIMENSION(:,:,:), INTENT(IN)                :: PQ
+REAL, DIMENSION(:,:,:)                            :: DQ
 !
 !*       0.2   Declarations of local variables :
 !   
@@ -1744,10 +1991,12 @@ DQ = 0.5 * DQ
 
 !$acc end data
 
-END SUBROUTINE DIF2Y
+END SUBROUTINE DIF2Y_DEVICE
 ! #endif
 !
 #ifdef MNH_OPENACC
+! END SUBROUTINE PPM_01_Y_D
+
 END SUBROUTINE PPM_01_Y
 #else
 END FUNCTION PPM_01_Y
@@ -1757,13 +2006,57 @@ END FUNCTION PPM_01_Y
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
+!     ########################################################################
+!!$      FUNCTION PPM_01_Z(KGRID, PSRC, PCR, PRHO, PTSTEP) RESULT(PR)
       SUBROUTINE PPM_01_Z(KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
+!     ########################################################################
+
+  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+  USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
+
+  IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)    :: PCR &    ! Courant number
+                                      ,  PRHO  ! density
+!
+! output source term
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+
+INTEGER :: IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG
+
+!$acc data present( PSRC, PCR, PRHO, PR )
+
+        CALL  MNH_GET_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)
+
+        CALL  PPM_01_Z_D(IIU,IJU,IKU, KGRID, &
+                     & PSRC, PCR, PRHO, PTSTEP, PR, &
+                     & ZT3D(:,:,:,IZQL),ZT3D(:,:,:,IZQR),ZT3D(:,:,:,IZDQ),ZT3D(:,:,:,IZQ6), &
+                     & ZT3D(:,:,:,IZDMQ),ZT3D(:,:,:,IZQL0),ZT3D(:,:,:,IZQR0), ZT3D(:,:,:,IZQ60), &
+                     & ZT3D(:,:,:,IZFPOS),ZT3D(:,:,:,IZFNEG) )
+
+        CALL  MNH_REL_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)
+
+!$acc end data
+
+CONTAINS
+!
+!     ########################################################################
+        SUBROUTINE  PPM_01_Z_D(IIU,IJU,IKU,KGRID, &
+                     & PSRC, PCR, PRHO, PTSTEP, PR, &
+                     & ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG)
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_01_Z(KGRID, PSRC, PCR, PRHO, PTSTEP) RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_01_Z - PPM_01 fully monotonic PPM advection scheme in Z direction
 !!                 Colella notation
@@ -1775,29 +2068,30 @@ END FUNCTION PPM_01_Y
 !!
 !-------------------------------------------------------------------------------
 !
-USE MODD_CONF
-USE MODD_PARAMETERS, ONLY: JPVEXT
-
 USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,  ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-USE MODE_MPPDB
 
-#ifdef MNH_BITREP
-USE MODI_BITREP
-#endif
-USE MODI_GET_HALO
 #ifndef MNH_OPENACC
 USE MODI_SHUMAN
 #else
 USE MODI_SHUMAN_DEVICE
 #endif
+USE MODI_GET_HALO
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
+!
+USE MODD_CONF
+USE MODD_PARAMETERS
+!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODE_MPPDB
 !
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
 !
+#ifdef MNH_OPENACC
+integer,                intent(in)    :: iiu, iju, iku
+#endif
 INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
 !
 #ifndef MNH_OPENACC
@@ -1814,16 +2108,16 @@ REAL,                   INTENT(IN)  :: PTSTEP  ! Time step
 #ifndef MNH_OPENACC
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
 #else
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+REAL, DIMENSION(IIU,IJU,IKU), INTENT(OUT) :: PR
 #endif
 !
 !*       0.2   Declarations of local variables :
 !
+#ifndef MNH_OPENACC
 INTEGER :: IIU, IJU, IKU
 INTEGER:: IKB    ! Begining useful area in x,y,z directions
 INTEGER:: IKE    ! End useful area in x,y,z directions
 !
-#ifndef MNH_OPENACC
 ! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL,ZQR
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZDQ,ZQ6
@@ -1836,15 +2130,17 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL0,ZQR0,ZQ60
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 #else
 ! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZQL, ZQR
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDQ, ZQ6
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDMQ
+REAL, DIMENSION(IIU,IJU,IKU) :: &
+                                                      ZQL, ZQR, ZDQ, ZQ6, ZDMQ &
 !
 ! extra variables for the initial guess of parabolae parameters
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZQL0,ZQR0,ZQ60
+                                                     , ZQL0,ZQR0,ZQ60 &
 !
 ! advection fluxes
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFPOS, ZFNEG
+                                                     , ZFPOS, ZFNEG
+!
+INTEGER:: IKB    ! Begining useful area in x,y,z directions
+INTEGER:: IKE    ! End useful area in x,y,z directions
 !
 INTEGER                          :: I,J,K 
 #endif
@@ -1852,6 +2148,8 @@ integer                          :: ji, jj, jk
 !
 !-------------------------------------------------------------------------------
 
+!$acc data present( PSRC, PCR, PRHO, PR, &
+!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK(PCR, "PPM_01_Z beg:PCR")
@@ -1859,40 +2157,23 @@ IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PSRC,"PPM_01_Z beg:PSRC")
 END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZQL,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQR,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZDQ,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQ6,   IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZDMQ,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQL0,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQR0,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZQ60,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFPOS, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFNEG, IIU, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PR, &
-!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
 !
 IKB = 1 + JPVEXT
 IKE = SIZE(PSRC,3) - JPVEXT
+#ifndef MNH_OPENACC
+iiu = size( PSRC, 1 )
+iju = size( PSRC, 2 )
+iku = size( PSRC, 3 )
+#endif
 
 !$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = 1, iiu
         PR   (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQL  (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQR  (ji, jj, jk ) = PSRC(ji, jj, jk )
@@ -1902,7 +2183,9 @@ DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
         ZQL0 (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQR0 (ji, jj, jk ) = PSRC(ji, jj, jk )
         ZQ60 (ji, jj, jk ) = PSRC(ji, jj, jk )
-END DO
+    end do
+  end do
+end do
 #ifndef MNH_OPENACC
 ZFPOS=PSRC
 ZFNEG=PSRC
@@ -1925,51 +2208,59 @@ ZFNEG(:,:,:) = PSRC(:,:,:)
 !               --------------------------------
 ! 
 ! calculate dmq
-  call DIF2Z( PSRC, ZDMQ )
+#ifndef  MNH_OPENACC
+ZDMQ = DIF2Z(PSRC)
+#else
+CALL DIF2Z_DEVICE(ZDMQ,PSRC)
+#endif
 !$acc kernels
 !
 ! monotonize the difference followinq eq. 5 in Lin94
 ! use the periodic BC here, it doesn't matter for vertical (hopefully) 
 !
-!$acc_nv loop independent collapse(3)
-  DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = ikb:ike )
+!$acc loop independent collapse(3)
+do jk = ikb, ike
+  do jj = 1, iju
+    do ji = 1, iiu
       ZDMQ(ji, jj, jk ) = SIGN(                                                                                    &
         MIN( ABS( ZDMQ(ji, jj, jk ) ),                                                                             &
              2.0 * (   PSRC(ji, jj, jk ) - MIN( PSRC(ji, jj, jk-1 ), PSRC(ji, jj, jk ), PSRC(ji, jj, jk+1 ) ) ) ,  &
              2.0 * ( - PSRC(ji, jj, jk ) + MAX( PSRC(ji, jj, jk-1 ), PSRC(ji, jj, jk ), PSRC(ji, jj, jk+1 ) ) ) ), &
         ZDMQ(ji, jj, jk ) )
-  END DO
-!$acc end kernels
-!$acc kernels
+    end do
+  end do
+end do
 ZDMQ(:,:,IKB-1) = &
      SIGN( (MIN( ABS(ZDMQ(:,:,IKB-1)), 2.0*(PSRC(:,:,IKB-1) - &
      MIN(PSRC(:,:,IKE-1),PSRC(:,:,IKB-1),PSRC(:,:,IKB))),   &
      2.0*(MAX(PSRC(:,:,IKE-1),PSRC(:,:,IKB-1),PSRC(:,:,IKB)) - &
      PSRC(:,:,IKB-1)) )), ZDMQ(:,:,IKB-1) )
-!$acc end kernels
-!$acc kernels
 ZDMQ(:,:,IKE+1) = &
      SIGN( (MIN( ABS(ZDMQ(:,:,IKE+1)), 2.0*(PSRC(:,:,IKE+1) - &
      MIN(PSRC(:,:,IKE),PSRC(:,:,IKE+1),PSRC(:,:,IKB+1))),  &
      2.0*(MAX(PSRC(:,:,IKE),PSRC(:,:,IKE+1),PSRC(:,:,IKB+1)) - &
      PSRC(:,:,IKE+1)) )), ZDMQ(:,:,IKE+1) )
-!$acc end kernels
-!$acc kernels
 !
 ! calculate qL and qR with the modified dmq
 !
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = 1:iju ,  jk = ikb:ike + 1 )
+!$acc loop independent collapse(3)
+do jk = ikb, ike + 1
+  do jj = 1, iju
+    do ji = 1, iiu
       ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj, jk-1 ) ) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj, jk-1 ) ) / 6.0
-END DO
+    end do
+  end do
+end do
 ZQL0(:,:,IKB-1) = ZQL0(:,:,IKE)
 !
-!$acc end kernels
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT (  ji = 1:iiu , jj = 1:iju , jk = ikb-1:ike )
+!$acc loop independent collapse(3)
+do jk = ikb - 1, ike
+  do jj = 1, iju
+    do ji = 1, iiu
       ZQR0(ji, jj, jk ) = ZQL0(ji, jj, jk+1 )
-END DO
+    end do
+  end do
+end do
 ZQR0(:,:,IKE+1) = ZQR0(:,:,IKB)
 #ifndef MNH_OPENACC
 !
@@ -2021,8 +2312,10 @@ END WHERE
 !
 ZDQ = ZQR - ZQL
 #else
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( I=1:IIU , J=1:IJU , K=1:IKU )
+!$acc loop independent collapse(3)
+   DO K=1,IKU
+      DO J=1,IJU
+         DO I=1,IIU
             !
             ! determine initial coefficients of the parabolae
             !
@@ -2055,21 +2348,23 @@ DO CONCURRENT ( I=1:IIU , J=1:IJU , K=1:IKU )
             !
             ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
             !
-END DO
+         END DO
+      END DO
+   END DO
 #endif
-!$acc end kernels
 !
 ! and finally calculate fluxes for the advection
 !
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = 1:iju ,  jk = ikb + 1: ike + 1 )
+!$acc loop independent collapse(3)
+do jk = ikb + 1, ike + 1
+  do jj = 1, iju
+    do ji = 1, iiu
       ZFPOS(ji, jj, jk ) = ZQR(ji, jj, jk-1 ) - 0.5 * PCR(ji, jj, jk ) &
                                    * ( ZDQ(ji, jj, jk-1 ) - ( 1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0) * ZQ6(ji, jj, jk-1 ) )
-END DO
-!$acc end kernels
+    end do
+  end do
+end do
 !
-!$acc kernels present_cr(ZFPOS,ZQR)
 ! advection flux at open boundary when u(IKB) > 0
 ZFPOS(:,:,IKB) = (PSRC(:,:,IKB-1) - ZQR(:,:,IKB-1))*PCR(:,:,IKB) + &
                  ZQR(:,:,IKB-1)
@@ -2078,16 +2373,16 @@ ZFPOS(:,:,IKB) = (PSRC(:,:,IKB-1) - ZQR(:,:,IKB-1))*PCR(:,:,IKB) + &
 ! we set it to 0
 ZFPOS(:,:,IKB-1) = 0.0
 !
-!$acc end kernels
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = ikb - 1: ike )
+!$acc loop independent collapse(3)
+do jk = ikb - 1, ike
+  do jj = 1, iju
+    do ji = 1, iiu
       ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5 * PCR(ji, jj, jk ) &
                                    * ( ZDQ(ji, jj, jk ) + ( 1.0 + 2.0 * PCR(ji, jj, jk ) / 3.0) * ZQ6(ji, jj, jk ) )
-END DO
+    end do
+  end do
+end do
 !
-!$acc end kernels
-!$acc kernels present_cr(ZFNEG,ZQR)
 ! advection flux at open boundary when u(IKE+1) < 0
 ZFNEG(:,:,IKE+1) = (ZQR(:,:,IKE)-PSRC(:,:,IKE+1))*PCR(:,:,IKE+1) + &
                    ZQR(:,:,IKE)
@@ -2101,17 +2396,21 @@ PR = DZF( PCR*MZM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + &
 !$acc end kernels 
     CALL MZM_DEVICE(PRHO,ZQL)
 !$acc kernels 
-!$acc_nv loop independent collapse(3)
-    DO CONCURRENT ( ji = 1:iiu , jj = 1:iju , jk = 1:iku )
+!$acc loop independent collapse(3)
+do jk = 1, iku
+  do jj = 1, iju
+    do ji = 1, iiu
       if ( PCR(ji, jj, jk ) > 0. ) then
         ZQR(ji, jj, jk ) =  PCR(ji, jj, jk ) * ZQL(ji, jj, jk ) * ZFPOS(ji, jj, jk )
       else
         ZQR(ji, jj, jk ) =  PCR(ji, jj, jk ) * ZQL(ji, jj, jk ) * ZFNEG(ji, jj, jk )
       end if
-   END DO
+    end do
+  end do
+end do
     !dzf(PR,ZQR)
 !$acc end kernels
-    CALL DZF_DEVICE( ZQR, PR )
+    CALL DZF_DEVICE(ZQR,PR)
 #endif
 !
 #ifndef MNH_OPENACC
@@ -2129,20 +2428,82 @@ END IF
 
 !$acc end data
 
-#ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-
+#ifndef MNH_OPENACC
 CONTAINS
+#else
+END SUBROUTINE PPM_01_Z_D
+#endif
 !
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
 !     ########################################################################
-      SUBROUTINE DIF2Z( PQ, DQ )
+      FUNCTION DIF2Z(PQ) RESULT(DQ)
 !     ########################################################################
+!!
+!!****  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(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: 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 FUNCTION DIF2Z
+!-------------------------------------------------------------------------------
+!
+!     ########################################################################
+      SUBROUTINE DIF2Z_DEVICE(DQ,PQ)
+!     ########################################################################
+!!
 !!****  DIF2Z - leap-frog difference operator in Z direction
 !!
 !!    Calculates the difference assuming periodic BC (CYCL). 
@@ -2166,9 +2527,8 @@ IMPLICIT NONE
 ! 
 !*       0.1   Declarations of dummy arguments :
 !   
-REAL, DIMENSION(:,:,:),                            INTENT(IN)  :: PQ
-REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)), INTENT(OUT) :: DQ
-!REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+REAL, DIMENSION(:,:,:), INTENT(IN)                :: PQ
+REAL, DIMENSION(:,:,:)                            :: DQ
 !
 !*       0.2   Declarations of local variables :
 !   
@@ -2200,9 +2560,13 @@ DQ = 0.5 * DQ
 
 !$acc end data
 
-END SUBROUTINE DIF2Z
+END SUBROUTINE DIF2Z_DEVICE
 
+! #endif
+!
 #ifdef MNH_OPENACC
+! END SUBROUTINE PPM_01_Z_D
+!
 END SUBROUTINE PPM_01_Z
 #else
 END FUNCTION PPM_01_Z
@@ -2212,14 +2576,59 @@ END FUNCTION PPM_01_Z
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
-      SUBROUTINE PPM_S0_X( HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP, PR )
+!     ########################################################################
+!!$      FUNCTION PPM_S0_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
+!!$               RESULT(PR)
+SUBROUTINE PPM_S0_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP , PR)
+
+  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+  USE MODE_MNH_ZWORK, ONLY : ZPSRC_HALO2_WEST
+
+  IMPLICIT NONE
+  !
+  !*       0.1   Declarations of dummy arguments :
+  !
+  CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
+  !
+  INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+  !
+  REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
+  REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+  REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+  REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
+  !
+  ! output source term
+  REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+
+  INTEGER     :: IZFPOS,IZFNEG,IZPHAT,IZRHO_MXM,IZCR_MXM,IZCR_DXF
+
+!$acc data present( PSRC, PCR, PRHO, PR )
+
+  CALL MNH_GET_ZT3D(IZFPOS,IZFNEG,IZPHAT,IZRHO_MXM,IZCR_MXM,IZCR_DXF)
+
+  CALL PPM_S0_X_D(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP , PR,  &
+               &  ZT3D(:,:,:,IZFPOS),ZT3D(:,:,:,IZFNEG),ZT3D(:,:,:,IZPHAT),ZT3D(:,:,:,IZRHO_MXM), &
+               &  ZT3D(:,:,:,IZCR_MXM),ZT3D(:,:,:,IZCR_DXF),ZPSRC_HALO2_WEST  )
+
+  CALL MNH_REL_ZT3D(IZFPOS,IZFNEG,IZPHAT,IZRHO_MXM,IZCR_MXM,IZCR_DXF)
+
+!$acc end data
+
+CONTAINS
+!
+!     ########################################################################
+      SUBROUTINE PPM_S0_X_D(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP , PR &
+                         & ,ZFPOS,ZPHAT,ZFNEG &
+                         & ,ZRHO_MXM,ZCR_MXM,ZCR_DXF,ZPSRC_HALO2_WEST )
+
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_S0_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
                RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_S0_X - PPM  advection scheme in X direction in Skamarock 2006 
 !!                 notation - NO CONSTRAINTS
@@ -2237,10 +2646,6 @@ USE MODD_CONF
 
 USE MODE_ll
 #ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,   ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-USE MODE_MPPDB
-#ifdef MNH_OPENACC
 use mode_msg
 #endif
 
@@ -2251,6 +2656,15 @@ USE MODI_SHUMAN
 USE MODI_SHUMAN_DEVICE
 #endif
 !
+#ifdef MNH_OPENACC
+USE MODD_PARAMETERS, ONLY : JPHEXT
+!
+USE MODE_MNH_ZWORK, ONLY : IIB,IIE, IIU,IJU,IKU , IJS,IJN, GWEST,GEAST
+!
+USE MODD_IO,   ONLY : GSMONOPROC
+#endif
+USE MODE_MPPDB
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -2274,14 +2688,14 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
 !
 !*       0.2   Declarations of local variables :
 !
+#ifndef MNH_OPENACC
 INTEGER:: IIB,IJB    ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE    ! End useful area in x,y,z directions
 INTEGER                          :: IJS,IJN
-integer :: iiu, iju, iku
+INTEGER :: IIU, IJU, IKU
 !
 LOGICAL :: GWEST, GEAST
 
-#ifndef MNH_OPENACC
 ! advection fluxes
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 !
@@ -2289,25 +2703,32 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT
 !
 REAL, DIMENSION(SIZE(PCR,2),SIZE(PCR,3))             :: ZPSRC_HALO2_WEST
-TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
 #else
 ! advection fluxes
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFPOS, ZFNEG
+REAL, DIMENSION(:,:,:) :: ZFPOS, ZFNEG
 !
 ! variable at cell edges
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPHAT
+REAL, DIMENSION(:,:,:) :: ZPHAT
 !
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZCR_DXF, ZCR_MXM, ZRHO_MXM
+REAL, DIMENSION(:,:,:) :: ZRHO_MXM, ZCR_MXM  , ZCR_DXF
 INTEGER                          :: I,J,K
 !
-REAL, DIMENSION(:,:),   POINTER, CONTIGUOUS :: ZPSRC_HALO2_WEST
-REAL, DIMENSION(:,:),   POINTER, CONTIGUOUS :: ZWEST
 LOGICAL, SAVE :: GFIRST_CALL_PPM_S0_X = .TRUE.
-TYPE(HALO2LIST_ll), SAVE , POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
+REAL, DIMENSION(:,:)             :: ZPSRC_HALO2_WEST
 #endif
 
+#ifndef MNH_OPENACC
+TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
+#else
+TYPE(HALO2LIST_ll), SAVE , POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
+!
+REAL , POINTER , CONTIGUOUS , DIMENSION(:,:) :: ZWEST
+#endif
 !-------------------------------------------------------------------------------
 
+!$acc data present( PSRC, PCR, PRHO, PR , &
+!$acc &             ZFPOS, ZFNEG, ZPHAT, ZRHO_MXM, ZCR_MXM, ZCR_DXF, ZPSRC_HALO2_WEST )
+
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK(PCR, "PPM_S0_X beg:PCR")
@@ -2315,31 +2736,11 @@ IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PSRC,"PPM_S0_X beg:PSRC")
 END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZFPOS,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFNEG,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZPHAT,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZCR_DXF,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZCR_MXM,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRHO_MXM, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZPSRC_HALO2_WEST, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PR , &
-!$acc &             ZFPOS, ZFNEG, ZPHAT, ZRHO_MXM, ZCR_MXM, ZCR_DXF, ZPSRC_HALO2_WEST )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
 !
+#ifndef MNH_OPENACC
 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
 IJS=IJB
 IJN=IJE
@@ -2349,6 +2750,11 @@ IJN=IJE
 GWEST = LWEST_ll()
 GEAST = LEAST_ll()
 !
+IIU = SIZE( PSRC, 1 )
+IJU = SIZE( PSRC, 2 )
+IKU = SIZE( PSRC, 3 )
+#endif
+!
 !BEG JUAN PPM_LL
 !
 !*              initialise & update halo & halo2 for PSRC
@@ -2362,9 +2768,9 @@ IF (GFIRST_CALL_PPM_S0_X) THEN
    NULLIFY(TZ_PSRC_HALO2_ll)
    CALL INIT_HALO2_ll(TZ_PSRC_HALO2_ll,1,IIU,IJU,IKU)
 END IF
-CALL GET_HALO2_DF(PSRC, TZ_PSRC_HALO2_ll, HNAME='GET_HALO2::PSRC')
+CALL GET_HALO2_DF(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC')
 ZWEST => TZ_PSRC_HALO2_ll%HALO2%WEST
-!$acc kernels present_cr(ZPSRC_HALO2_WEST,ZWEST)
+!$acc kernels
 ZPSRC_HALO2_WEST(:,:) = ZWEST(:,:)
 !$acc end kernels
 #endif
@@ -2458,6 +2864,7 @@ CALL GET_HALO(ZFNEG, HNAME='ZFNEG') ! JUAN
    CALL GET_HALO(PR, HNAME='PR') ! JUAN
 !
 CASE ('OPEN')
+!$acc kernels
 !
 !!$   ZPHAT(IIB,:,:) = 0.5*(PSRC(IIB-1,:,:) + PSRC(IIB,:,:))
 !!$   ZPHAT(IIB-1,:,:) = ZPHAT(IIB,:,:)   ! not used
@@ -2466,25 +2873,24 @@ CASE ('OPEN')
 !  WEST BOUND 
 !
 IF (.NOT. GWEST) THEN
-!$acc kernels
    ZPHAT(IIB  ,IJS:IJN,:) = ( 7.0 * &
                       ( PSRC(IIB  ,IJS:IJN,:) + PSRC(IIB-1,IJS:IJN,:)                  ) - &
                       ( PSRC(IIB+1,IJS:IJN,:) + ZPSRC_HALO2_WEST(IJS:IJN,:)            ) ) / 12.0
 ! <=>  WEST BOUND     ( PSRC(IIB+1,IJS:IJN,:) + PSRC(IIB-2,IJS:IJN,:)                  ) ) / 12.0
-!$acc end kernels
 ENDIF
+!$acc end kernels
 !
 !   update ZPHAT HALO before next/further  utilisation 
 !
 #ifndef MNH_OPENACC
-!PW: a remettre? CALL  GET_HALO(ZPHAT, HNAME='ZPHAT')
+CALL  GET_HALO(ZPHAT, HNAME='ZPHAT')
 #else
 ! acc update self(ZPHAT)
 !CALL GET_HALO_D(ZPHAT(:,:,:), HDIR="Z0_X", HNAME='ZPHAT')
 ! acc update device(ZPHAT)
 #endif
 !
-!$acc kernels present_cr(ZFPOS,ZPHAT)
+!$acc kernels
   IF (GWEST) THEN
    ZPHAT(IIB  ,IJS:IJN,:) = 0.5*(PSRC(IIB-1,IJS:IJN,:) + PSRC(IIB,IJS:IJN,:))
    ZPHAT(IIB-1,IJS:IJN,:) = ZPHAT(IIB,IJS:IJN,:)
@@ -2512,24 +2918,22 @@ ENDIF
 !$acc end kernels
 !
 #ifndef MNH_OPENACC
-!PW: a remettre? CALL GET_HALO(ZFPOS, HNAME='ZFPOS') ! JUAN
+CALL GET_HALO(ZFPOS, HNAME='ZFPOS') ! JUAN
 #else
 ! acc update self(ZFPOS)
 !CALL GET_HALO_D(ZFPOS(:,:,:), HDIR="Z0_X", HNAME='ZFPOS') ! JUAN
 ! acc update device(ZFPOS)
 #endif
 !
+!$acc kernels
 ! positive flux on the WEST boundary
   IF (GWEST) THEN
-!$acc kernels present_cr(ZFPOS,ZPHAT)
    ZFPOS(IIB,IJS:IJN,:) = (PSRC(IIB-1,IJS:IJN,:) - ZPHAT(IIB,IJS:IJN,:))*PCR(IIB,IJS:IJN,:) + &
                      ZPHAT(IIB,IJS:IJN,:) 
 ! this is not used
    ZFPOS(IIB-1,IJS:IJN,:) = 0.0
-!$acc end kernels
   ENDIF
 !
-!$acc kernels present_cr(ZFNEG,ZPHAT)
 ! negative fluxes
 !!$   ZFNEG(IIB:IIE,:,:) = ZPHAT(IIB:IIE,:,:) + & 
 !!$        PCR(IIB:IIE,:,:)*(ZPHAT(IIB:IIE,:,:) - PSRC(IIB:IIE,:,:)) + &
@@ -2542,15 +2946,15 @@ ENDIF
 !$acc end kernels
 !
 #ifndef MNH_OPENACC
-!PW: a remettre?    CALL GET_HALO(ZFNEG, HNAME='ZFNEG') ! JUAN
+   CALL GET_HALO(ZFNEG, HNAME='ZFNEG') ! JUAN
 #else
 ! acc update self(ZFNEG)
 !CALL GET_HALO_D(ZFNEG, HDIR="Z0_X", HNAME='ZFNEG') ! JUAN
 ! acc update device(ZFNEG)
 #endif
 !
+!$acc kernels
   IF (GEAST) THEN
-!$acc kernels present_cr(ZFNEG,ZPHAT)
 !
 ! in OPEN case PCR(IIB-1) is not used, so we also set ZFNEG(IIB-1) = 0
 !
@@ -2561,7 +2965,6 @@ ENDIF
 !
    ZFNEG(IIE+1,IJS:IJN,:) = (ZPHAT(IIE+1,IJS:IJN,:) - PSRC(IIE+1,IJS:IJN,:))*PCR(IIE+1,IJS:IJN,:) + &
                        ZPHAT(IIE+1,IJS:IJN,:)
-!$acc end kernels
   ENDIF
 !
 ! calculate the advection
@@ -2571,38 +2974,35 @@ ENDIF
         DXF( PCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 #else
+!$acc end kernels
    CALL MXM_DEVICE(PRHO,ZRHO_MXM)
-!$acc kernels present_cr(ZCR_MXM,ZRHO_MXM,ZFPOS,ZFNEG)
+!$acc kernels
    ZCR_MXM =  PCR * ZRHO_MXM * ( ZFPOS*(0.5+SIGN(0.5,PCR)) + ZFNEG*(0.5-SIGN(0.5,PCR)) ) 
 !$acc end kernels
    CALL DXF_DEVICE(ZCR_MXM,ZCR_DXF)
 !$acc kernels
    PR = PSRC * PRHO - ZCR_DXF
-!$acc end kernels
 #endif
 !
 ! in OPEN case fix boundary conditions
 !
   IF (GWEST) THEN
-!$acc kernels
    WHERE ( PCR(IIB,IJS:IJN,:) <= 0. ) !  OUTFLOW condition
       PR(IIB-1,IJS:IJN,:) = 2.*PR(IIB,IJS:IJN,:) - PR(IIB+1,IJS:IJN,:)
    ELSEWHERE
       PR(IIB-1,IJS:IJN,:) = PR(IIB,IJS:IJN,:)
    END WHERE
-!$acc end kernels
   ENDIF
 !
   IF (GEAST) THEN 
-!$acc kernels
    WHERE ( PCR(IIE,IJS:IJN,:) >= 0. ) !  OUTFLOW condition
       PR(IIE+1,IJS:IJN,:) = 2.*PR(IIE,IJS:IJN,:) - PR(IIE-1,IJS:IJN,:)
    ELSEWHERE
       PR(IIE+1,IJS:IJN,:) = PR(IIE,IJS:IJN,:)
    END WHERE
-!$acc end kernels 
   ENDIF
 !
+!$acc end kernels 
 !
 !
 END SELECT
@@ -2610,8 +3010,7 @@ END SELECT
 #ifndef MNH_OPENACC
 CALL GET_HALO(PR, HNAME='PR')
 #else
-! CALL GET_HALO_D(PR, HDIR="S0_X", HNAME='PR')
-CALL GET_HALO_D(PR, HDIR="S0_X", HNAME='UPDATE_HALO_ll::GET_HALO::PR')
+CALL GET_HALO_D(PR, HDIR="S0_X", HNAME='PR')
 #endif
 !-------------------------------------------------------------------------------
 #ifndef MNH_OPENACC
@@ -2628,11 +3027,8 @@ END IF
 !$acc end data
 
 #ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
+END SUBROUTINE PPM_S0_X_D
 
-#ifdef MNH_OPENACC
 END SUBROUTINE PPM_S0_X
 #else
 END FUNCTION PPM_S0_X
@@ -2641,14 +3037,59 @@ END FUNCTION PPM_S0_X
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
+!     ########################################################################
+!!$      FUNCTION PPM_S0_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP) &
+!!$               RESULT(PR)
          SUBROUTINE PPM_S0_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
+
+  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+  USE MODE_MNH_ZWORK, ONLY : ZPSRC_HALO2_SOUTH
+
+  IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
+!
+! output source term
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
+
+  INTEGER     :: IZFPOS,IZPHAT,IZFNEG,IZRHO_MYM,IZCR_MYM,IZCR_DYF
+
+!$acc data present( PSRC, PCR, PRHO, PR )
+
+  CALL MNH_GET_ZT3D(IZFPOS,IZFNEG,IZPHAT,IZRHO_MYM,IZCR_MYM,IZCR_DYF)
+
+  CALL PPM_S0_Y_D(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP , PR,  &
+               &  ZT3D(:,:,:,IZFPOS),ZT3D(:,:,:,IZFNEG),ZT3D(:,:,:,IZPHAT),ZT3D(:,:,:,IZRHO_MYM), &
+               &  ZT3D(:,:,:,IZCR_MYM),ZT3D(:,:,:,IZCR_DYF),ZPSRC_HALO2_SOUTH  )
+
+  CALL MNH_REL_ZT3D(IZFPOS,IZFNEG,IZPHAT,IZRHO_MYM,IZCR_MYM,IZCR_DYF)
+
+!$acc end data
+
+CONTAINS
+!
+!     ########################################################################
+      SUBROUTINE PPM_S0_Y_D(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP , PR &
+                         & ,ZFPOS,ZPHAT,ZFNEG &
+                         & ,ZRHO_MYM,ZCR_MYM,ZCR_DYF,ZPSRC_HALO2_SOUTH )
+
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_S0_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP) &
                RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_S0_Y - PPM  advection scheme in Y direction in Skamarock 2006 
 !!                 notation - NO CONSTRAINTS
@@ -2664,10 +3105,6 @@ USE MODD_CONF
 
 USE MODE_ll
 #ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-USE MODE_MPPDB
-#ifdef MNH_OPENACC
 use mode_msg
 #endif
 
@@ -2678,6 +3115,16 @@ USE MODI_SHUMAN
 USE MODI_SHUMAN_DEVICE
 #endif
 !
+#ifdef MNH_OPENACC
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODD_PARAMETERS, ONLY : JPHEXT
+!
+USE MODE_MNH_ZWORK, ONLY : IJB,IJE, IIU,IJU,IKU , IIW,IIA, GSOUTH , GNORTH
+!
+USE MODD_IO,   ONLY : GSMONOPROC
+#endif
+USE MODE_MPPDB
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -2701,14 +3148,15 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
 !
 !*       0.2   Declarations of local variables :
 !
+#ifndef MNH_OPENACC
 INTEGER:: IIB,IJB    ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE    ! End useful area in x,y,z directions
+INTEGER                          :: IJS,IJN
 INTEGER                          :: IIW,IIA
 INTEGER :: IIU, IJU, IKU
 !
 LOGICAL :: GNORTH, GSOUTH
 !
-#ifndef MNH_OPENACC
 ! advection fluxes
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 !
@@ -2719,28 +3167,37 @@ TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
 TYPE(HALO2LIST_ll), POINTER      :: TZ_PHAT_HALO2_ll         ! halo2 for ZPHAT
 !
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,3))             :: ZPSRC_HALO2_SOUTH
-TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
 #else
+!
 ! advection fluxes
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFPOS, ZFNEG
+REAL, DIMENSION(:,:,:) :: ZFPOS, ZFNEG
 !
 ! variable at cell edges
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPHAT
+REAL, DIMENSION(:,:,:) :: ZPHAT
 !
+#ifndef MNH_OPENACC
+TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
+#else
 TYPE(HALO2LIST_ll), SAVE ,POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
+
 REAL , POINTER , CONTIGUOUS , DIMENSION(:,:) :: ZSOUTH
+#endif
+
 TYPE(HALO2LIST_ll), POINTER      :: TZ_PHAT_HALO2_ll         ! halo2 for ZPHAT
 !
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZCR_DYF, ZCR_MYM, ZRHO_MYM
+REAL, DIMENSION(:,:,:) :: ZRHO_MYM , ZCR_MYM , ZCR_DYF
 !
 INTEGER                          :: I,J,K
 !
 LOGICAL, SAVE :: GFIRST_CALL_PPM_S0_Y = .TRUE.
-REAL, DIMENSION(:,:),   POINTER, CONTIGUOUS :: ZPSRC_HALO2_SOUTH
+REAL, DIMENSION(:,:)             :: ZPSRC_HALO2_SOUTH
 #endif
 !
 !-------------------------------------------------------------------------------
 
+!$acc data present( PSRC, PCR, PRHO, PR , &
+!$acc &             ZFPOS, ZFNEG, ZPHAT, ZRHO_MYM, ZCR_MYM, ZCR_DYF, ZPSRC_HALO2_SOUTH )
+
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK(PCR, "PPM_S0_Y beg:PCR")
@@ -2748,31 +3205,11 @@ IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PSRC,"PPM_S0_Y beg:PSRC")
 END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZFPOS,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFNEG,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZPHAT,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZCR_DYF,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZCR_MYM,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRHO_MYM, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZPSRC_HALO2_SOUTH, IIU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PR , &
-!$acc &             ZFPOS, ZFNEG, ZPHAT, ZRHO_MYM, ZCR_MYM, ZCR_DYF, ZPSRC_HALO2_SOUTH )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
 !
+#ifndef MNH_OPENACC
 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
 IIW=IIB
 IIA=IIE
@@ -2782,6 +3219,11 @@ IIA=IIE
 GNORTH = LNORTH_ll()
 GSOUTH = LSOUTH_ll()
 !
+IIU = SIZE( PSRC, 1 )
+IJU = SIZE( PSRC, 2 )
+IKU = SIZE( PSRC, 3 )
+#endif
+!
 !-------------------------------------------------------------------------------
 !
 IF ( L2D ) THEN
@@ -2803,9 +3245,9 @@ IF (GFIRST_CALL_PPM_S0_Y) THEN
    NULLIFY(TZ_PSRC_HALO2_ll)
    CALL INIT_HALO2_ll(TZ_PSRC_HALO2_ll,1,IIU,IJU,IKU)
 END IF   
-CALL GET_HALO2_DF(PSRC, TZ_PSRC_HALO2_ll, HNAME='GET_HALO2::PSRC')
+CALL GET_HALO2_DF(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC')
 ZSOUTH => TZ_PSRC_HALO2_ll%HALO2%SOUTH(:,:)
-!$acc kernels present_cr(ZPSRC_HALO2_SOUTH,ZSOUTH)
+!$acc kernels
 ZPSRC_HALO2_SOUTH(:,:) = ZSOUTH(:,:)
 !$acc end kernels
 #endif
@@ -2825,7 +3267,7 @@ PR=PSRC
 ZPHAT(IIW:IIA,IJB+1:IJE,:) = (7.0 * &
                        (PSRC(IIW:IIA,IJB+1:IJE,:) + PSRC(IIW:IIA,IJB:IJE-1,:)) - &
                        (PSRC(IIW:IIA,IJB+2:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE-2,:))) / 12.0
-!$acc end kernels
+!$acc end kernels 
 !
 SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
 CASE ('CYCL','WALL')            ! In that case one must have HLBCY(1) == HLBCY(2)
@@ -2893,6 +3335,7 @@ CALL GET_HALO(ZFNEG, HNAME='ZFNEG') ! JUAN
 #endif
 !
 CASE ('OPEN')
+!$acc kernels 
 !
 !!$   ZPHAT(:,IJB,:) = 0.5*(PSRC(:,IJB-1,:) + PSRC(:,IJB,:))
 !!$   ZPHAT(:,IJB-1,:) = ZPHAT(:,IJB,:)   ! not used
@@ -2902,38 +3345,34 @@ CASE ('OPEN')
 !  SOUTH BOUND 
 !
   IF ( .NOT. GSOUTH) THEN
-!$acc kernels
    ZPHAT(IIW:IIA,IJB  ,:) = (7.0 * &
                       (PSRC(IIW:IIA,IJB  ,:) + PSRC(IIW:IIA,IJB-1,:)) - &
                       (PSRC(IIW:IIA,IJB+1,:) + ZPSRC_HALO2_SOUTH(IIW:IIA,:)            )) / 12.0
 !                     (PSRC(IIW:IIA,IJB+1,:) + TZ_PSRC_HALO2_ll%HALO2%SOUTH(IIW:IIA,:) )) / 12.0
 ! <=> SOUTH BOUND     (PSRC(IIW:IIA,IJB+1,:) + PSRC(IIW:IIA,IJB-2,:)                   )) / 12.0
-!$acc end kernels
   ENDIF
 !
 !TEMPO_JUAN  
+!$acc end kernels
 !
 #ifndef MNH_OPENACC
-!PW: a remettre? CALL  GET_HALO(ZPHAT, HNAME='ZPHAT')
+CALL  GET_HALO(ZPHAT, HNAME='ZPHAT')
 #else
 ! acc update self(ZPHAT)
 !CALL  GET_HALO_D(ZPHAT(:,:,:), HDIR="Z0_Y", HNAME='ZPHAT')
 ! acc update device(ZPHAT)
 #endif
 !
-  IF (GSOUTH) THEN
 !$acc kernels
+  IF (GSOUTH) THEN
    ZPHAT(IIW:IIA,IJB  ,:) = 0.5*(PSRC(IIW:IIA,IJB-1,:) + PSRC(IIW:IIA,IJB,:))
    ZPHAT(IIW:IIA,IJB-1,:) = ZPHAT(IIW:IIA,IJB,:)
-!$acc end kernels
   ENDIF
 !
 ! NORTH BOUND
 !
   IF (GNORTH) THEN
-!$acc kernels present_cr(ZPHAT)
    ZPHAT(IIW:IIA,IJE+1,:) =  0.5*(PSRC(IIW:IIA,IJE,:) + PSRC(IIW:IIA,IJE+1,:))
-!$acc end kernels
   ENDIF
 !
 !
@@ -2947,30 +3386,28 @@ CASE ('OPEN')
 !!$        PCR(:,IJB+1:IJE+1,:)*(ZPHAT(:,IJB+1:IJE+1,:) - PSRC(:,IJB:IJE,:)) - &
 !!$        PCR(:,IJB+1:IJE+1,:)*(1.0 - PCR(:,IJB+1:IJE+1,:)) * &
 !!$        (ZPHAT(:,IJB:IJE,:) - 2.0*PSRC(:,IJB:IJE,:) + ZPHAT(:,IJB+1:IJE+1,:))
-!$acc kernels present_cr(ZFPOS,ZPHAT)
-   ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZPHAT(IIW:IIA,IJB:IJE+1,:) - &
+   ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZPHAT(IIW:IIA,IJB:IJE+1,:) - & 
         PCR(IIW:IIA,IJB:IJE+1,:)*( ZPHAT(IIW:IIA,IJB:IJE+1,:) - PSRC(IIW:IIA,IJB-1:IJE  ,:) ) - &
         PCR(IIW:IIA,IJB:IJE+1,:)*( 1.0                  -  PCR(IIW:IIA,IJB  :IJE+1,:) ) * &
         (ZPHAT(IIW:IIA,IJB-1:IJE,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE,:) + ZPHAT(IIW:IIA,IJB:IJE+1,:))
 !$acc end kernels
 !
 #ifndef MNH_OPENACC
-!PW: a remettre? CALL GET_HALO(ZFPOS, HNAME='ZFPOS') ! JUAN
+CALL GET_HALO(ZFPOS, HNAME='ZFPOS') ! JUAN
 #else
 ! acc update self(ZFPOS)
 !CALL GET_HALO_D(ZFPOS(:,:,:), HDIR="Z0_Y", HNAME='ZFPOS') ! JUAN
 ! acc update device(ZFPOS)
 #endif
 !
+!$acc kernels
 ! positive flux on the SOUTH boundary
   IF (GSOUTH) THEN
-!$acc kernels present_cr(ZFPOS,ZPHAT)
    ZFPOS(IIW:IIA,IJB,:) = (PSRC(IIW:IIA,IJB-1,:) - ZPHAT(IIW:IIA,IJB,:))*PCR(IIW:IIA,IJB,:) + &
                      ZPHAT(IIW:IIA,IJB,:)
 !
 ! this is not used
    ZFPOS(IIW:IIA,IJB-1,:) = 0.0
-!$acc end kernels
   ENDIF
 ! 
 ! negative fluxes
@@ -2978,30 +3415,28 @@ CASE ('OPEN')
 !!$        PCR(:,IJB:IJE,:)*(ZPHAT(:,IJB:IJE,:) - PSRC(:,IJB:IJE,:)) + &
 !!$        PCR(:,IJB:IJE,:)*(1.0 + PCR(:,IJB:IJE,:)) * &
 !!$        (ZPHAT(:,IJB:IJE,:) - 2.0*PSRC(:,IJB:IJE,:) +ZPHAT(:,IJB+1:IJE+1,:))
-!$acc kernels present_cr(ZFNEG,ZPHAT)
-   ZFNEG(IIW:IIA,IJB-1:IJE,:) = ZPHAT(IIW:IIA,IJB-1:IJE,:) + &
+   ZFNEG(IIW:IIA,IJB-1:IJE,:) = ZPHAT(IIW:IIA,IJB-1:IJE,:) + & 
         PCR(IIW:IIA,IJB-1:IJE,:)*(ZPHAT(IIW:IIA,IJB-1:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE,:)) + &
         PCR(IIW:IIA,IJB-1:IJE,:)*(1.0 + PCR(IIW:IIA,IJB-1:IJE,:)) * &
         (ZPHAT(IIW:IIA,IJB-1:IJE,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE,:) +ZPHAT(IIW:IIA,IJB:IJE+1,:))
 !$acc end kernels
 !
 #ifndef MNH_OPENACC
-!PW: a remettre?    CALL GET_HALO(ZFNEG, HNAME='ZFNEG') ! JUAN
+   CALL GET_HALO(ZFNEG, HNAME='ZFNEG') ! JUAN
 #else
 ! acc update self(ZFNEG)
 !   CALL GET_HALO_D(ZFNEG, HDIR="Z0_Y", HNAME='ZFNEG') ! JUAN
 ! acc update device(ZFNEG)
 #endif
 !
+!$acc kernels
   IF (GNORTH) THEN
-!$acc kernels present_cr(ZFNEG,ZPHAT)
 ! this is not used
    ZFNEG(IIW:IIA,IJB-1,:) = 0.0
 !
 ! negative flux on the NORTH boundary
    ZFNEG(IIW:IIA,IJE+1,:) = (ZPHAT(IIW:IIA,IJE+1,:) - PSRC(IIW:IIA,IJE+1,:))*PCR(IIW:IIA,IJE+1,:) + &
                        ZPHAT(IIW:IIA,IJE+1,:)
-!$acc end kernels
   ENDIF
 !
 ! calculate the advection
@@ -3011,46 +3446,43 @@ CASE ('OPEN')
         DYF( PCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 #else
+!$acc end kernels
    CALL MYM_DEVICE(PRHO,ZRHO_MYM)
-!$acc kernels present_cr(ZCR_MYM,ZRHO_MYM,ZFPOS,ZFNEG)
-   ZCR_MYM =  PCR * ZRHO_MYM * ( ZFPOS(:,:,:)*(0.5+SIGN(0.5,PCR)) + ZFNEG*(0.5-SIGN(0.5,PCR)) )
+!$acc kernels
+   ZCR_MYM =  PCR* ZRHO_MYM*( ZFPOS(:,:,:)*(0.5+SIGN(0.5,PCR)) + ZFNEG*(0.5-SIGN(0.5,PCR)) ) 
 !$acc end kernels
    CALL DYF_DEVICE(ZCR_MYM,ZCR_DYF)
 !$acc kernels
    PR = PSRC * PRHO - ZCR_DYF
-!$acc end kernels
 #endif
 !
 ! in OPEN case fix boundary conditions
 !
   IF (GSOUTH) THEN
-!$acc kernels
    WHERE ( PCR(IIW:IIA,IJB,:) <= 0. ) !  OUTFLOW condition
       PR(IIW:IIA,IJB-1,:) = 1.0 * 2.*PR(IIW:IIA,IJB,:) - PR(IIW:IIA,IJB+1,:)
    ELSEWHERE
       PR(IIW:IIA,IJB-1,:) = PR(IIW:IIA,IJB,:) 
    END WHERE
-!$acc end kernels
   ENDIF
 !
   IF (GNORTH) THEN
-!$acc kernels
    WHERE ( PCR(IIW:IIA,IJE,:) >= 0. ) !  OUTFLOW condition
       PR(IIW:IIA,IJE+1,:) = 1.0 * 2.*PR(IIW:IIA,IJE,:) - PR(IIW:IIA,IJE-1,:)
    ELSEWHERE
       PR(IIW:IIA,IJE+1,:) = PR(IIW:IIA,IJE,:) 
    END WHERE
-!$acc end kernels
   ENDIF
 !
+!$acc end kernels 
+!
 !
 END SELECT
 !
 #ifndef MNH_OPENACC
 CALL GET_HALO(PR, HNAME='PR')
 #else
-! CALL GET_HALO_D(PR, HDIR="S0_Y", HNAME='PR')
-CALL GET_HALO_D(PR, HDIR="S0_Y", HNAME='UPDATE_HALO_ll::GET_HALO::PR')
+CALL GET_HALO_D(PR, HDIR="S0_Y", HNAME='PR')
 #endif
 !
 #ifndef MNH_OPENACC
@@ -3068,11 +3500,8 @@ END IF !not L2D
 !$acc end data
 
 #ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
+END SUBROUTINE PPM_S0_Y_D
 
-#ifdef MNH_OPENACC
 END SUBROUTINE PPM_S0_Y
 #else
 END FUNCTION PPM_S0_Y
@@ -3082,14 +3511,58 @@ END FUNCTION PPM_S0_Y
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
-      SUBROUTINE PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
+!     ########################################################################
+!!$      FUNCTION PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP) &
+!!$               RESULT(PR)
+SUBROUTINE PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
+
+  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+
+  IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+!
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
+!
+! output source term
+REAL, DIMENSION(:,:,:),INTENT(OUT):: PR
+
+
+  INTEGER     :: IZFPOS,IZPHAT,IZFNEG,IZRHO_MZM,IZCR_MZM,IZCR_DZF
+
+!$acc data present ( PSRC, PCR, PRHO, PR )
+
+  CALL MNH_GET_ZT3D(IZFPOS,IZFNEG,IZPHAT,IZRHO_MZM,IZCR_MZM,IZCR_DZF)
+
+  CALL PPM_S0_Z_D(KGRID, PSRC, PCR, PRHO, PTSTEP , PR,  &
+               &  ZT3D(:,:,:,IZFPOS),  ZT3D(:,:,:,IZFNEG),  ZT3D(:,:,:,IZPHAT), &
+               &  ZT3D(:,:,:,IZRHO_MZM),ZT3D(:,:,:,IZCR_MZM),ZT3D(:,:,:,IZCR_DZF)  )
+
+  CALL MNH_REL_ZT3D(IZFPOS,IZFNEG,IZPHAT,IZRHO_MZM,IZCR_MZM,IZCR_DZF)
+
+!$acc end data
+
+CONTAINS
+!
+!     ########################################################################
+SUBROUTINE PPM_S0_Z_D(KGRID, PSRC, PCR, PRHO, PTSTEP , PR  &
+                   & ,ZFPOS,ZFNEG,ZPHAT  &
+                   & ,ZRHO_MZM,ZCR_MZM,ZCR_DZF )
+
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP) &
                RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_S0_Z - PPM  advection scheme in Z direction in Skamarock 2006 
 !!                 notation - NO CONSTRAINTS
@@ -3101,22 +3574,22 @@ END FUNCTION PPM_S0_Y
 !!
 !-------------------------------------------------------------------------------
 !
-USE MODD_CONF
-USE MODD_PARAMETERS
-
 USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-USE MODE_MPPDB
-
-USE MODI_GET_HALO
 #ifndef MNH_OPENACC
 USE MODI_SHUMAN
 #else
 USE MODI_SHUMAN_DEVICE
 #endif
-
+USE MODI_GET_HALO
+!
+USE MODD_CONF
+USE MODD_PARAMETERS
+USE MODE_MPPDB
+!
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK, ONLY : IKB,IKE, IKU
+#endif
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -3138,11 +3611,10 @@ REAL, DIMENSION(:,:,:),INTENT(OUT):: PR
 !
 !*       0.2   Declarations of local variables :
 !
+#ifndef MNH_OPENACC
 INTEGER:: IKB    ! Begining useful area in x,y,z directions
 INTEGER:: IKE    ! End useful area in x,y,z directions
-INTEGER :: IIU, IJU, IKU
 !
-#ifndef MNH_OPENACC
 ! advection fluxes
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 !
@@ -3150,48 +3622,33 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT
 #else
 ! advection fluxes
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFPOS, ZFNEG
+REAL, DIMENSION(:,:,:),INTENT(OUT):: ZFPOS, ZFNEG &
 !
 ! interpolated variable at cell edges
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPHAT
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZCR_DZF, ZCR_MZM, ZRHO_MZM
+  &                                                   , ZPHAT &
+  &                                                   , ZRHO_MZM ,ZCR_MZM,ZCR_DZF
 #endif
 !
 !-------------------------------------------------------------------------------
 
-IF (MPPDB_INITIALIZED) THEN
-  !Check all IN arrays
-  CALL MPPDB_CHECK(PCR, "PPM_S0_Z beg:PCR")
-  CALL MPPDB_CHECK(PRHO,"PPM_S0_Z beg:PRHO")
-  !Check all INOUT arrays
-  CALL MPPDB_CHECK(PSRC,"PPM_S0_Z beg:PSRC")
-END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZFPOS,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFNEG,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZPHAT,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZCR_DZF,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZCR_MZM,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRHO_MZM, IIU, IJU, IKU )
-#endif
-
 !$acc data present ( PSRC, PCR, PRHO, PR , &
 !$acc &              ZFPOS, ZFNEG, ZPHAT, ZRHO_MZM, ZCR_MZM, ZCR_DZF )
 
+IF (MPPDB_INITIALIZED) THEN
+  !Check all IN arrays
+  CALL MPPDB_CHECK(PCR, "PPM_S0_Z beg:PCR")
+  CALL MPPDB_CHECK(PRHO,"PPM_S0_Z beg:PRHO")
+  !Check all INOUT arrays
+  CALL MPPDB_CHECK(PSRC,"PPM_S0_Z beg:PSRC")
+END IF
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
 !
+#ifndef MNH_OPENACC
 IKB = 1 + JPVEXT
 IKE = SIZE(PSRC,3) - JPVEXT
+#endif
 !
 !-------------------------------------------------------------------------------
 !
@@ -3200,10 +3657,12 @@ IKE = SIZE(PSRC,3) - JPVEXT
 #ifndef MNH_OPENACC
    CALL GET_HALO(PSRC, HNAME='PSRC')
 #else
-   CALL GET_HALO_D(PSRC, HNAME='UPDATE_HALO_ll::GET_HALO::PSRC')
+   CALL GET_HALO_D(PSRC, HNAME='PSRC')
 #endif
 !
-!$acc kernels present_cr(ZFPOS,ZPHAT,ZFNEG)
+#ifdef MNH_OPENACC
+!$acc kernels
+#endif
 !
 ZPHAT(:,:,IKB+1:IKE) = (7.0 * &
                        (PSRC(:,:,IKB+1:IKE) + PSRC(:,:,IKB:IKE-1)) - &
@@ -3258,11 +3717,11 @@ PR = PSRC * PRHO - &
 #else
 !$acc end kernels
    CALL MZM_DEVICE(PRHO,ZRHO_MZM)
-!$acc kernels present_cr(ZCR_MZM,ZRHO_MZM,ZFPOS,ZFNEG)
+!$acc kernels
    ZCR_MZM =  PCR* ZRHO_MZM*( ZFPOS(:,:,:)*(0.5+SIGN(0.5,PCR)) + ZFNEG*(0.5-SIGN(0.5,PCR)) ) 
    !dzf(ZCR_DZF,ZCR_MZM)
 !$acc end kernels
-   CALL DZF_DEVICE( ZCR_MZM, ZCR_DZF )
+   CALL DZF_DEVICE(ZCR_MZM,ZCR_DZF)
 !$acc kernels
    PR = PSRC * PRHO - ZCR_DZF
 #endif
@@ -3274,14 +3733,10 @@ PR = PSRC * PRHO - &
 !
 !$acc end kernels 
 !
-CALL MPPDB_CHECK(ZPHAT,"PPM_S0_Z:ZPHAT")
-CALL MPPDB_CHECK(ZFPOS,"PPM_S0_Z:ZFPOS")
-CALL MPPDB_CHECK(ZFNEG,"PPM_S0_Z:ZFNEG")
-CALL MPPDB_CHECK(PR,"PPM_S0_Z:PR")
 #ifndef MNH_OPENACC
    CALL GET_HALO(PR, HNAME='PR')
 #else
-   CALL GET_HALO_D(PR, HNAME='UPDATE_HALO_ll::GET_HALO::PR')
+   CALL GET_HALO_D(PR, HNAME='PR')
 #endif
 IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
@@ -3293,11 +3748,8 @@ END IF
 !$acc end data
 
 #ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
+END SUBROUTINE PPM_S0_Z_D
 
-#ifdef MNH_OPENACC
 END SUBROUTINE PPM_S0_Z
 #else
 END FUNCTION PPM_S0_Z
@@ -3306,14 +3758,71 @@ END FUNCTION PPM_S0_Z
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
-      SUBROUTINE PPM_S1_X(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, PR)
+!     ########################################################################
+!      FUNCTION PPM_S1_X(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
+!                        PTSTEP) RESULT(PR)
+      SUBROUTINE PPM_S1_X(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
+                          PTSTEP, PR)
+!     ########################################################################
+USE MODE_ll
+use mode_msg
+USE MODE_IO
+USE MODI_SHUMAN_DEVICE
+!
+USE MODD_CONF
+USE MODD_LUNIT
+USE MODD_PARAMETERS
+!
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOT ! density at t+dt
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step
+!
+! output source term
+REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
+
+INTEGER :: IZPHAT,IZRUT,IZFUP,IZFCOR,IZRPOS,IZRNEG
+
+!$acc data present( PSRC, PCR, PRHO, PRHOT, PR )
+
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_X', 'OpenACC: not yet implemented' )
+
+    CALL  MNH_GET_ZT3D(IZPHAT,IZRUT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
+
+    CALL PPM_S1_X_D(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, PR, &
+                    ZT3D(:,:,:,IZPHAT),ZT3D(:,:,:,IZRUT),ZT3D(:,:,:,IZFUP), &
+                    ZT3D(:,:,:,IZFCOR),ZT3D(:,:,:,IZRPOS),ZT3D(:,:,:,IZRNEG) )
+
+    CALL  MNH_REL_ZT3D(IZPHAT,IZRUT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
+
+!$acc end data
+
+  CONTAINS
+!
+!     ########################################################################
+!      FUNCTION PPM_S1_X_D(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
+!                        PTSTEP) RESULT(PR)
+      SUBROUTINE PPM_S1_X_D(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
+                          PTSTEP, PR, ZPHAT,ZRUT,ZFUP,ZFCOR,ZRPOS,ZRNEG)
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_S1_X(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
                         PTSTEP) RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_S1_X - PPM  advection scheme in X direction in Skamarock 2006 
 !!                 notation - with flux limiting for monotonicity
@@ -3325,13 +3834,7 @@ END FUNCTION PPM_S0_Z
 !!
 !-------------------------------------------------------------------------------
 !
-USE MODD_CONF
-USE MODD_PARAMETERS
-
 USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
 use mode_mppdb
 
 #ifndef MNH_OPENACC
@@ -3340,7 +3843,11 @@ USE MODI_SHUMAN
 USE MODI_SHUMAN
 USE MODI_SHUMAN_DEVICE
 #endif
-
+!
+USE MODD_CONF
+USE MODD_PARAMETERS
+!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -3370,9 +3877,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
 !
 INTEGER:: IIB,IJB,IKB    ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE,IKE    ! End useful area in x,y,z directions
-INTEGER :: IIU, IJU, IKU
 !
-#ifndef MNH_OPENACC
 ! variable at cell edges
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT, ZRUT
 !
@@ -3381,16 +3886,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFUP, ZFCOR
 !
 ! ratios for limiting the correction flux
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZRPOS, ZRNEG
-#else
-! variable at cell edges
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPHAT, ZRUT
-!
-! advection fluxes, upwind and correction
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFUP, ZFCOR
-!
-! ratios for limiting the correction flux
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRPOS, ZRNEG
-#endif
 !
 ! variables for limiting the correction flux
 REAL :: ZSRCMAX, ZSRCMIN, ZFIN, ZFOUT
@@ -3401,9 +3896,9 @@ INTEGER :: II, IJ, IK
 INTEGER                          :: IRESP             ! for prints
 !
 !-------------------------------------------------------------------------------
-#ifdef MNH_OPENACC
-call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_X', 'OpenACC: not yet implemented' )
-#endif
+
+!$acc data present( PSRC, PCR, PRHO, PRHOT, PR, &
+!$acc &             ZPHAT, ZRUT, ZFUP, ZFCOR, ZRPOS, ZRNEG )
 
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
@@ -3413,26 +3908,6 @@ IF (MPPDB_INITIALIZED) THEN
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PSRC, "PPM_S1_X beg:PSRC")
 END IF
-
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZPHAT, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRUT,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFUP,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFCOR, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRPOS, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRNEG, IIU, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PRHOT, PR, &
-!$acc &             ZPHAT, ZRUT, ZFUP, ZFCOR, ZRPOS, ZRNEG )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
@@ -3626,11 +4101,7 @@ END IF
 !$acc end data
 
 #ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-
-#ifdef MNH_OPENACC
+  END SUBROUTINE PPM_S1_X_D
 END SUBROUTINE PPM_S1_X
 #else
 END FUNCTION PPM_S1_X
@@ -3639,15 +4110,71 @@ END FUNCTION PPM_S1_X
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
-      SUBROUTINE PPM_S1_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PRHOT, &
+!     ########################################################################
+!      FUNCTION PPM_S1_Y(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
+!                        PTSTEP) RESULT(PR)
+      SUBROUTINE PPM_S1_Y(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
                           PTSTEP, PR)
+!     ########################################################################
+USE MODE_ll
+USE MODE_IO
+use mode_msg
+USE MODI_SHUMAN_DEVICE
+!
+USE MODD_CONF
+USE MODD_LUNIT
+USE MODD_PARAMETERS
+!
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOT ! density at t+dt
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step
+!
+! output source term
+REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
+
+INTEGER :: IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG
+
+!$acc data present( PSRC, PCR, PRHO, PRHOT, PR )
+
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_Y', 'OpenACC: not yet implemented' )
+
+    CALL  MNH_GET_ZT3D(IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
+
+    CALL PPM_S1_Y_D(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, PR, &
+                    ZT3D(:,:,:,IZPHAT),ZT3D(:,:,:,IZRVT),ZT3D(:,:,:,IZFUP), &
+                    ZT3D(:,:,:,IZFCOR),ZT3D(:,:,:,IZRPOS),ZT3D(:,:,:,IZRNEG) )
+
+    CALL  MNH_REL_ZT3D(IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
+
+!$acc end data
+
+  CONTAINS
+!
+!     ########################################################################
+!      FUNCTION PPM_S1_Y_D(HLBCY, KGRID, PSRC, PCR, PRHO, PRHOT, &
+!                        PTSTEP) RESULT(PR)
+      SUBROUTINE PPM_S1_Y_D(HLBCY, KGRID, PSRC, PCR, PRHO, PRHOT, &
+                          PTSTEP, PR, ZPHAT,ZRVT,ZFUP,ZFCOR,ZRPOS,ZRNEG)
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_S1_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PRHOT, &
                         PTSTEP) RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_S1_Y - PPM  advection scheme in Y direction in Skamarock 2006 
 !!                 notation - with flux limiting for monotonicity
@@ -3659,13 +4186,7 @@ END FUNCTION PPM_S1_X
 !!
 !-------------------------------------------------------------------------------
 !
-USE MODD_CONF
-USE MODD_PARAMETERS
-
 USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
 use mode_mppdb
 
 #ifndef MNH_OPENACC
@@ -3675,6 +4196,12 @@ USE MODI_SHUMAN
 USE MODI_SHUMAN_DEVICE
 #endif
 !
+USE MODD_CONF
+USE MODD_PARAMETERS
+!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+#endif
 !
 IMPLICIT NONE
 !
@@ -3705,9 +4232,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
 !
 INTEGER:: IIB,IJB,IKB   ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE,IKE   ! End useful area in x,y,z directions
-INTEGER :: IIU, IJU, IKU
 !
-#ifndef MNH_OPENACC
 ! variable at cell edges
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT, ZRVT
 !
@@ -3716,16 +4241,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFUP, ZFCOR
 !
 ! ratios for limiting the correction flux
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZRPOS, ZRNEG
-#else
-! variable at cell edges
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPHAT, ZRVT
-!
-! advection fluxes, upwind and correction
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFUP, ZFCOR
-!
-! ratios for limiting the correction flux
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRPOS, ZRNEG
-#endif
 !
 ! variables for limiting the correction flux
 REAL :: ZSRCMAX, ZSRCMIN, ZFIN, ZFOUT
@@ -3737,9 +4252,9 @@ INTEGER :: II, IJ, IK
 INTEGER  :: IRESP   ! Return code of FM-routines
 !
 !-------------------------------------------------------------------------------
-#ifdef MNH_OPENACC
-call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_Y', 'OpenACC: not yet implemented' )
-#endif
+
+!$acc data present( PSRC, PCR, PRHO, PRHOT, PR , &
+!$acc &             ZPHAT, ZRVT, ZFUP, ZFCOR, ZRPOS, ZRNEG )
 
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
@@ -3750,25 +4265,7 @@ IF (MPPDB_INITIALIZED) THEN
   CALL MPPDB_CHECK(PSRC, "PPM_S1_Y beg:PSRC")
 END IF
 
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZPHAT, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRVT,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFUP,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFCOR, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRPOS, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRNEG, IIU, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PRHOT, PR , &
-!$acc &             ZPHAT, ZRVT, ZFUP, ZFCOR, ZRPOS, ZRNEG )
-
+!
 IF ( L2D ) THEN
    PR = PSRC*PRHO
    !RETURN
@@ -3965,11 +4462,7 @@ END IF !not L2D
 !$acc end data
 
 #ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-
-#ifdef MNH_OPENACC
+  END SUBROUTINE PPM_S1_Y_D
 END SUBROUTINE PPM_S1_Y
 #else
 END FUNCTION PPM_S1_Y
@@ -3978,14 +4471,68 @@ END FUNCTION PPM_S1_Y
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !
-!     ########################################################################
 #ifdef MNH_OPENACC
-      SUBROUTINE PPM_S1_Z(KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, PR)
+!
+!     ########################################################################
+!      FUNCTION PPM_S1_Z(KGRID, PSRC, PCR, PRHO, PRHOT, &
+!                        PTSTEP) RESULT(PR)
+      SUBROUTINE PPM_S1_Z(KGRID, PSRC, PCR, PRHO, PRHOT, &
+                          PTSTEP, PR)
+!     ########################################################################
+USE MODE_ll
+USE MODE_IO
+use mode_msg
+
+USE MODI_SHUMAN_DEVICE
+!
+USE MODD_CONF
+USE MODD_LUNIT
+USE MODD_PARAMETERS
+!
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
+REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOT ! density at t+dt
+REAL,                   INTENT(IN)  :: PTSTEP  ! Time step
+!
+! output source term
+REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
+
+INTEGER :: IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG
+
+!$acc data present( PSRC, PCR, PRHO, PRHOT, PR )
+
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_Z', 'OpenACC: not yet implemented' )
+
+    CALL  MNH_GET_ZT3D(IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
+
+    CALL PPM_S1_Z_D(KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, PR, &
+                    ZT3D(:,:,:,IZPHAT),ZT3D(:,:,:,IZRVT),ZT3D(:,:,:,IZFUP), &
+                    ZT3D(:,:,:,IZFCOR),ZT3D(:,:,:,IZRPOS),ZT3D(:,:,:,IZRNEG) )
+
+    CALL  MNH_REL_ZT3D(IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
+
+!$acc end data
+
+  CONTAINS
+!     ########################################################################
+      SUBROUTINE PPM_S1_Z_D(KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, &
+                          PR, ZPHAT,ZRVT,ZFUP,ZFCOR,ZRPOS,ZRNEG)
+!     ########################################################################
 #else
+!     ########################################################################
       FUNCTION PPM_S1_Z(KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP) &
                         RESULT(PR)
-#endif
 !     ########################################################################
+#endif
 !!
 !!****  PPM_S1_Z - PPM  advection scheme in Z direction in Skamarock 2006 
 !!                 notation - with flux limiting for monotonicity
@@ -3997,13 +4544,7 @@ END FUNCTION PPM_S1_Y
 !!
 !-------------------------------------------------------------------------------
 !
-USE MODD_CONF
-USE MODD_PARAMETERS
-
 USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
 use mode_mppdb
 
 #ifndef MNH_OPENACC
@@ -4012,7 +4553,14 @@ USE MODI_SHUMAN
 USE MODI_SHUMAN
 USE MODI_SHUMAN_DEVICE
 #endif
-
+!
+USE MODD_CONF
+USE MODD_PARAMETERS
+!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+#endif
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -4040,9 +4588,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
 !
 INTEGER:: IIB,IJB,IKB   ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE,IKE   ! End useful area in x,y,z directions
-INTEGER :: IIU, IJU, IKU
 !
-#ifndef MNH_OPENACC
 ! variable at cell edges
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT, ZRVT
 !
@@ -4051,16 +4597,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFUP, ZFCOR
 !
 ! ratios for limiting the correction flux
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZRPOS, ZRNEG
-#else
-! variable at cell edges
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPHAT, ZRVT
-!
-! advection fluxes, upwind and correction
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFUP, ZFCOR
-!
-! ratios for limiting the correction flux
-REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRPOS, ZRNEG
-#endif
 !
 ! variables for limiting the correction flux
 REAL :: ZSRCMAX, ZSRCMIN, ZFIN, ZFOUT
@@ -4070,9 +4606,9 @@ REAL, PARAMETER :: ZEPS = 1.0E-16
 INTEGER :: II, IJ, IK
 !
 !-------------------------------------------------------------------------------
-#ifdef MNH_OPENACC
-call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_Z', 'OpenACC: not yet implemented' )
-#endif
+
+!$acc data present( PSRC, PCR, PRHO, PRHOT, PR, &
+!$acc &             ZPHAT, ZRVT, ZFUP, ZFCOR, ZRPOS, ZRNEG )
 
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
@@ -4083,25 +4619,6 @@ IF (MPPDB_INITIALIZED) THEN
   CALL MPPDB_CHECK(PSRC, "PPM_S1_Z beg:PSRC")
 END IF
 
-iiu = size( PSRC, 1 )
-iju = size( PSRC, 2 )
-iku = size( PSRC, 3 )
-
-#ifdef MNH_OPENACC
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZPHAT, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRVT,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFUP,  IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZFCOR, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRPOS, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZRNEG, IIU, IJU, IKU )
-#endif
-
-!$acc data present( PSRC, PCR, PRHO, PRHOT, PR, &
-!$acc &             ZPHAT, ZRVT, ZFUP, ZFCOR, ZRPOS, ZRNEG )
-
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
@@ -4359,11 +4876,7 @@ END IF
 !$acc end data
 
 #ifdef MNH_OPENACC
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-
-#ifdef MNH_OPENACC
+  END SUBROUTINE PPM_S1_Z_D
 END SUBROUTINE PPM_S1_Z
 #else
 END FUNCTION PPM_S1_Z
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/communication.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/communication.f90
index 6d5c71c75..8ee4340ab 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/communication.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/communication.f90
@@ -501,7 +501,7 @@ contains
 !  Scalar product of two fields
 !==================================================================
   subroutine scalarprod_mnh(m, a, b, s)
-#ifdef MNH
+#ifdef MNH_MPPDB
     USE MODE_MPPDB
 #endif    
     implicit none
@@ -561,9 +561,7 @@ contains
            za_st => a%st
            zb_st => b%st
            !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop collapse(3)
-#endif
+           !$acc loop collapse(3) reduction(+:local_sumt)
            do iz=0,nz+1
               do iy=icompy_min,icompy_max
                  do ix=icompx_min,icompx_max
@@ -1135,49 +1133,37 @@ contains
            ! Send to south
            if (Gneighbour_s) then
            ztab_halo_st_haloTin => tab_halo_st(level,m)%haloTin
-           !$acc kernels async(IS_SOUTH)
-#ifdef MNH_COMPILER_NVHPC            
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent ( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )
+           !$acc kernels async(IS_SOUTH) present_cr(zst,ztab_halo_st_haloTin)
+           !$mnh_do_concurrent( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )
               ztab_halo_st_haloTin(ii,ij,ik) = zst(ii,ij+a_n-halo_size,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
            ! Send to north
            if (Gneighbour_n) then
            ztab_halo_nt_haloTin => tab_halo_nt(level,m)%haloTin
-           !$acc kernels async(IS_NORTH)
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent ( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )           
+           !$acc kernels async(IS_NORTH) present_cr(zst,ztab_halo_nt_haloTin)
+           !$mnh_do_concurrent( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )           
               ztab_halo_nt_haloTin(ii,ij,ik) = zst(ii,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
            ! Send to east
            if (Gneighbour_e) then
            ztab_halo_et_haloTin => tab_halo_et(level,m)%haloTin
-           !$acc kernels async(IS_EAST)
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent ( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 ) 
+           !$acc kernels async(IS_EAST) present_cr(zst,ztab_halo_et_haloTin)
+           !$mnh_do_concurrent( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 ) 
               ztab_halo_et_haloTin(ii,ij,ik) = zst(ii+a_n-halo_size,ij-halo_size,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
            ! Send to west
            if (Gneighbour_w) then
            ztab_halo_wt_haloTin => tab_halo_wt(level,m)%haloTin
-           !$acc kernels async(IS_WEST)
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent ( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 ) 
+           !$acc kernels async(IS_WEST) present_cr(zst,ztab_halo_wt_haloTin)
+           !$mnh_do_concurrent( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 ) 
               ztab_halo_wt_haloTin(ii,ij,ik) = zst(ii,ij-halo_size,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
         end if
@@ -1386,46 +1372,34 @@ contains
         if (LUseT) then
            if (Gneighbour_n) then
            ! copy north halo for GPU managed
-           !$acc kernels async(IS_NORTH)
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif              
-           do concurrent ( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )
+           !$acc kernels async(IS_NORTH) present_cr(zst,ztab_halo_nt_haloTout)
+           !$mnh_do_concurrent( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )
               zst(ii,ij-halo_size,ik-1) = ztab_halo_nt_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
            if (Gneighbour_s) then
            ! copy south halo for GPU managed
-           !$acc kernels async(IS_SOUTH)
-#ifdef MNH_COMPILER_NVHPC              
-           !$acc loop independent collapse(3)   
-#endif              
-           do concurrent ( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )
+           !$acc kernels async(IS_SOUTH) present_cr(zst,ztab_halo_st_haloTout)
+           !$mnh_do_concurrent( ii=1:a_n,ij=1:halo_size,ik=1:nz+2 )
               zst(ii,ij+a_n,ik-1) = ztab_halo_st_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
            if (Gneighbour_w) then
            ! copy west halo for GPU managed
-           !$acc kernels async(IS_WEST)
-#ifdef MNH_COMPILER_NVHPC              
-           !$acc loop independent collapse(3)
-#endif              
-           do concurrent ( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 )
+           !$acc kernels async(IS_WEST) present_cr(zst,ztab_halo_wt_haloTout)
+           !$mnh_do_concurrent( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 )
               zst(ii-halo_size,ij-halo_size,ik-1) = ztab_halo_wt_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            end if
            if (Gneighbour_e) then
            ! copy east halo for GPU managed
-           !$acc kernels async(IS_EAST)
-#ifdef MNH_COMPILER_NVHPC
-           !$acc loop independent collapse(3)
-#endif              
-           do concurrent ( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 )
+           !$acc kernels async(IS_EAST) present_cr(zst,ztab_halo_et_haloTout)
+           !$mnh_do_concurrent( ii=1:halo_size,ij=1:a_n+2*halo_size,ik=1:nz+2 )
               zst(ii+a_n,ij-halo_size,ik-1) = ztab_halo_et_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            end if 
            ! wait for async copy of send buffer to GPU
@@ -1724,13 +1698,10 @@ contains
 #ifdef MNH_GPUDIRECT           
            zb_st => b%st
            za_st => a%st
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:a_n,ij=1:a_n,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,za_st)
+           !$mnh_do_concurrent(ii=1:a_n,ij=1:a_n,ik=1:nz+2)
               zb_st(ii,ij,ik-1) = za_st(ii,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
 #else           
            b%st(1:a_n,1:a_n,0:nz+1) = a%st(1:a_n,1:a_n,0:nz+1)
@@ -1744,32 +1715,23 @@ contains
            zb_st => b%st
            ! copy from buffer for GPU DIRECT
            ! Receive from NE
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC
-           !$acc loop independent collapse(3)
-#endif
-           do concurrent (ii=1:b_n/2,ij=1:b_n/2,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,ztab_sub_interiorT_ne_m_1_haloTout)
+           !$mnh_do_concurrent(ii=1:b_n/2,ij=1:b_n/2,ik=1:nz+2)
               zb_st(ii+b_n/2,ij,ik-1) = ztab_sub_interiorT_ne_m_1_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            ! Receive from SW
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif
-           do concurrent (ii=1:b_n/2,ij=1:b_n/2,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,ztab_sub_interiorT_sw_m_1_haloTout)
+           !$mnh_do_concurrent(ii=1:b_n/2,ij=1:b_n/2,ik=1:nz+2)
               zb_st(ii,ij+b_n/2,ik-1) = ztab_sub_interiorT_sw_m_1_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            ! Receive from SE
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-            do concurrent (ii=1:b_n/2,ij=1:b_n/2,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,ztab_sub_interiorT_se_m_1_haloTout)
+           !$mnh_do_concurrent(ii=1:b_n/2,ij=1:b_n/2,ik=1:nz+2)
                zb_st(ii+b_n/2,ij+b_n/2,ik-1) = ztab_sub_interiorT_se_m_1_haloTout(ii,ij,ik)
-            end do
-            !$acc end kernels
+           !$mnh_end_do()
+           !$acc end kernels
         end if
 #endif
       end if
@@ -1788,13 +1750,10 @@ contains
         if (LUseT) then
 #ifdef MNH_GPUDIRECT
            ztab_interiorT_ne_m_haloTin => tab_interiorT_ne(level,m)%haloTin
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif
-           do concurrent (ii=1:a_n,ij=1:a_n,ik=1:nz+2)
+           !$acc kernels present_cr(ztab_interiorT_ne_m_haloTin,za_st)
+           !$mnh_do_concurrent(ii=1:a_n,ij=1:a_n,ik=1:nz+2)
               ztab_interiorT_ne_m_haloTin(ii,ij,ik) = za_st(ii,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            !$acc host_data use_device(ztab_interiorT_ne_m_haloTin)
            call mpi_send(ztab_interiorT_ne_m_haloTin,size(ztab_interiorT_ne_m_haloTin), &
@@ -1824,13 +1783,10 @@ contains
         if (LUseT) then
 #ifdef MNH_GPUDIRECT
            ztab_interiorT_sw_m_haloTin => tab_interiorT_sw(level,m)%haloTin
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif
-           do concurrent (ii=1:a_n,ij=1:a_n,ik=1:nz+2)
+           !$acc kernels present_cr(ztab_interiorT_sw_m_haloTin,za_st)
+           !$mnh_do_concurrent(ii=1:a_n,ij=1:a_n,ik=1:nz+2)
               ztab_interiorT_sw_m_haloTin(ii,ij,ik) = za_st(ii,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            !$acc host_data use_device(ztab_interiorT_sw_m_haloTin)
            call mpi_send(ztab_interiorT_sw_m_haloTin,size(ztab_interiorT_sw_m_haloTin), &
@@ -1860,13 +1816,10 @@ contains
         if (LUseT) then
 #ifdef MNH_GPUDIRECT
            ztab_interiorT_se_m_haloTin => tab_interiorT_se(level,m)%haloTin
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif
-           do concurrent (ii=1:a_n,ij=1:a_n,ik=1:nz+2)
+           !$acc kernels present_cr(ztab_interiorT_se_m_haloTin,za_st)
+           !$mnh_do_concurrent(ii=1:a_n,ij=1:a_n,ik=1:nz+2)
               ztab_interiorT_se_m_haloTin(ii,ij,ik) = za_st(ii,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
            !$acc host_data use_device(ztab_interiorT_se_m_haloTin)
            call mpi_send(ztab_interiorT_se_m_haloTin,size(ztab_interiorT_se_m_haloTin), &
@@ -1971,13 +1924,10 @@ contains
         if (LUseT) then
 #ifdef MNH_GPUDIRECT
            ztab_sub_interiorT_ne_m_1_haloTin => tab_sub_interiorT_ne(level,m-1)%haloTin
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif
-           do concurrent (ii=1:a_n/2,ij=1:a_n/2,ik=1:nz+2)
+           !$acc kernels present_cr(ztab_sub_interiorT_ne_m_1_haloTin,za_st)
+           !$mnh_do_concurrent(ii=1:a_n/2,ij=1:a_n/2,ik=1:nz+2)
               ztab_sub_interiorT_ne_m_1_haloTin(ii,ij,ik) = za_st(ii+a_n/2,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            !$acc host_data use_device(ztab_sub_interiorT_ne_m_1_haloTin)
            call mpi_isend(ztab_sub_interiorT_ne_m_1_haloTin,size(ztab_sub_interiorT_ne_m_1_haloTin), &
@@ -2005,13 +1955,10 @@ contains
         if (LUseT) then
 #ifdef MNH_GPUDIRECT
            ztab_sub_interiorT_sw_m_1_haloTin => tab_sub_interiorT_sw(level,m-1)%haloTin
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:a_n/2,ij=1:a_n/2,ik=1:nz+2)
+           !$acc kernels present_cr(ztab_sub_interiorT_sw_m_1_haloTin,za_st)
+           !$mnh_do_concurrent(ii=1:a_n/2,ij=1:a_n/2,ik=1:nz+2)
               ztab_sub_interiorT_sw_m_1_haloTin(ii,ij,ik) = za_st(ii,ij+a_n/2,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            !$acc host_data use_device(ztab_sub_interiorT_sw_m_1_haloTin)
            call mpi_isend(ztab_sub_interiorT_sw_m_1_haloTin,size(ztab_sub_interiorT_sw_m_1_haloTin), &
@@ -2040,13 +1987,10 @@ contains
         if (LUseT) then
 #ifdef MNH_GPUDIRECT
            ztab_sub_interiorT_se_m_1_haloTin => tab_sub_interiorT_se(level,m-1)%haloTin
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:a_n/2,ij=1:a_n/2,ik=1:nz+2)
+           !$acc kernels present_cr(ztab_sub_interiorT_se_m_1_haloTin,za_st)
+           !$mnh_do_concurrent(ii=1:a_n/2,ij=1:a_n/2,ik=1:nz+2)
               ztab_sub_interiorT_se_m_1_haloTin(ii,ij,ik) = za_st(ii+a_n/2,ij+a_n/2,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
            !$acc host_data use_device(ztab_sub_interiorT_se_m_1_haloTin)
            call mpi_isend(ztab_sub_interiorT_se_m_1_haloTin,size(ztab_sub_interiorT_se_m_1_haloTin), &
@@ -2068,13 +2012,10 @@ contains
 #ifdef MNH_GPUDIRECT                     
            zb_st => b%st
            za_st => a%st
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:b_n,ij=1:b_n,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,za_st)
+           !$mnh_do_concurrent(ii=1:b_n,ij=1:b_n,ik=1:nz+2)
               zb_st(ii,ij,ik-1) = za_st(ii,ij,ik-1)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
 #else
            b%st(1:b_n,1:b_n,0:nz+1) = a%st(1:b_n,1:b_n,0:nz+1)
@@ -2104,13 +2045,10 @@ contains
            call mpi_recv(ztab_interiorT_ne_m_haloTout,size(ztab_interiorT_ne_m_haloTout), &
                 MPI_DOUBLE_PRECISION,source_rank,recv_tag,MPI_COMM_HORIZ,stat,ierr)
            !$acc end host_data
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:b_n,ij=1:b_n,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,ztab_interiorT_ne_m_haloTout)
+           !$mnh_do_concurrent(ii=1:b_n,ij=1:b_n,ik=1:nz+2)
               zb_st(ii,ij,ik-1) = ztab_interiorT_ne_m_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels           
 #else
            call mpi_recv(b%st(1,1,0),1,interiorT(level,m),source_rank,recv_tag,MPI_COMM_HORIZ,stat,ierr)           
@@ -2140,13 +2078,10 @@ contains
            call mpi_recv(ztab_interiorT_sw_m_haloTout,size(ztab_interiorT_sw_m_haloTout), &           
                 MPI_DOUBLE_PRECISION,source_rank,recv_tag,MPI_COMM_HORIZ,stat,ierr)
            !$acc end host_data
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:b_n,ij=1:b_n,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,ztab_interiorT_sw_m_haloTout)
+           !$mnh_do_concurrent(ii=1:b_n,ij=1:b_n,ik=1:nz+2)
               zb_st(ii,ij,ik-1) = ztab_interiorT_sw_m_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
 #else
            call mpi_recv(b%st(1,1,0),1,interiorT(level,m),source_rank,recv_tag,MPI_COMM_HORIZ,stat,ierr)           
@@ -2176,13 +2111,10 @@ contains
            call mpi_recv(ztab_interiorT_se_m_haloTout,size(ztab_interiorT_se_m_haloTout), &
                 MPI_DOUBLE_PRECISION,source_rank,recv_tag,MPI_COMM_HORIZ,stat,ierr)
            !$acc end host_data
-           !$acc kernels
-#ifdef MNH_COMPILER_NVHPC           
-           !$acc loop independent collapse(3)
-#endif           
-           do concurrent (ii=1:b_n,ij=1:b_n,ik=1:nz+2)
+           !$acc kernels present_cr(zb_st,ztab_interiorT_se_m_haloTout)
+           !$mnh_do_concurrent(ii=1:b_n,ij=1:b_n,ik=1:nz+2)
               zb_st(ii,ij,ik-1) = ztab_interiorT_se_m_haloTout(ii,ij,ik)
-           end do
+           !$mnh_end_do()
            !$acc end kernels
 #else
            call mpi_recv(b%st(1,1,0),1,interiorT(level,m),source_rank,recv_tag,MPI_COMM_HORIZ,stat,ierr)           
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90
index f99505180..0a99735ec 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90
@@ -360,8 +360,8 @@ private
        
        if (LUseT) then
           zphi_st => phi%st
-          !$acc kernels
-          !$acc loop collapse(3)
+          !$acc kernels present_cr(zphi_st)
+          !$acc loop collapse(3) reduction(+:tmpt) 
           do iz=1,nz
              do iy=1,nlocaly
                 do ix=1,nlocalx
@@ -477,7 +477,7 @@ private
     integer, intent(in) :: comm_horiz
     type(scalar3d), intent(in) :: phi
     character(*), intent(in)   :: filename
-    integer :: file_id = 100
+    integer :: file_id = 200
     integer :: ix,iy,iz
     integer :: nlocal
     integer :: rank, nproc, ierr
@@ -510,7 +510,7 @@ private
     write(file_id,'(" icompy_max  = ",I10)') phi%icompy_max
     write(file_id,'(" halosize    = ",I10)') phi%halo_size
 
-
+    if (LUseO) then
     do ix=1-phi%halo_size,nlocal+phi%halo_size
       do iy=1-phi%halo_size,nlocal+phi%halo_size
         do iz=0,phi%grid_param%nz+1
@@ -518,6 +518,19 @@ private
         end do
       end do
     end do
+    endif
+
+    if (LUseT) then
+       !$acc update host(phi%st)
+       do iz=0,phi%grid_param%nz+1
+          do iy=1-phi%halo_size,nlocal+phi%halo_size
+             do ix=1-phi%halo_size,nlocal+phi%halo_size
+                write(file_id,'(E24.15)') phi%st(ix,iy,iz)
+             end do
+          end do
+       end do
+    endif
+ 
     close(file_id)
   end subroutine save_scalar3d
 
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
index 471e25482..92c5800f2 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
@@ -262,7 +262,8 @@ contains
 !!$       print*,"zt1d_discretisation_allocate3d:: nt1d_discretisation_max=",nt1d_discretisation_max
 !!$       call flush(6)
     endif
-            
+
+    kindex = nt1d_discretisation_pointer_top
   end function zt1d_discretisation_allocate3d
 
   
@@ -792,18 +793,17 @@ end subroutine construct_vertical_coeff
        zc_k => vert_coeff%c
        zd_k => vert_coeff%d
     
-       !$acc kernels 
+       !$acc kernels present_cr(zb_k,zv_st)
        iz=1
-       !$acc_nv loop independent  collapse(2)
-       do concurrent(ii=iib:iie,ij=ijb:ije)
-             zv_st(ii,ij,iz)   = zd_k(iz)* ( (-zb_k(iz)-zc_k(iz))*Tij * zu_st(ii,ij,iz  )  &
-                                              +zb_k(iz)          *Tij * zu_st(ii,ij,iz+1)  )
-       end do
+       !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije )
+          zv_st(ii,ij,iz)   = zd_k(iz)* ( (-zb_k(iz)-zc_k(iz))*Tij * zu_st(ii,ij,iz  )  &
+               +zb_k(iz)          *Tij * zu_st(ii,ij,iz+1)  )
+       !$mnh_end_do()
+       
        !
 !!$       !$acc loop seq
-!!$       do iz=2,ize-1 
-          !$acc_nv loop independent  collapse(3)
-       do concurrent(ii=iib:iie,ij=ijb:ije,iz=2:ize-1)
+!!$       do iz=2,ize-1       
+       !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije , iz=2:ize-1 )
                 zv_st(ii,ij,iz) = zd_k(iz)* ( ((-zb_k(iz)-zc_k(iz))*Tij - 4.0_rl ) * zu_st(ii,ij,iz)    &
                                                 +zb_k(iz)          *Tij            * zu_st(ii,ij,iz+1)  &
                                                          +zc_k(iz) *Tij            * zu_st(ii,ij,iz-1)  &
@@ -812,15 +812,14 @@ end subroutine construct_vertical_coeff
                                            +                                         zu_st(ii,ij+1,iz) &
                                            +                                         zu_st(ii,ij-1,iz) &
                                        )
-       end do
+       !$mnh_end_do()
 !!$          end do
        !
        iz=ize
-       !$acc_nv loop independent  collapse(2)
-       do concurrent(ii=iib:iie,ij=ijb:ije)       
+       !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije )       
              zv_st(ii,ij,iz) = zd_k(iz)*  (  (-zb_k(iz)-zc_k(iz))*Tij  * zu_st(ii,ij,iz)    &
                                                        +zc_k(iz) *Tij  * zu_st(ii,ij,iz-1)  )
-       end do
+       !$mnh_end_do()
        !$acc end kernels
     endif
     
@@ -1437,6 +1436,7 @@ end subroutine construct_vertical_coeff
     integer :: nib,nie,njb,nje,nkb,nke
     integer :: ii,ij,ik
     integer :: iindex
+    integer :: n_lin 
 
     !
     !  init size , param
@@ -1444,6 +1444,8 @@ end subroutine construct_vertical_coeff
     nz = u%grid_param%nz
     nlocal = u%ix_max -u%ix_min + 1
     halo_size = u%halo_size
+    n_lin = (nlocal+2*halo_size)**2 * (nz+2)
+
  
     nib = 1-halo_size
     nie = nlocal+halo_size
@@ -1573,12 +1575,12 @@ end subroutine construct_vertical_coeff
     if (LUseO) u0(:,:,:) = u%s(:,:,:)
     if (LUseT) then
        zu_st =>  u%st
-       !$acc kernels
-       !$acc_nv loop independent collapse(3)
-       do concurrent (ii=nib:nie,ij=njb:nje,ik=nkb:nke)
-          ut0(ii,ij,ik) = zu_st(ii,ij,ik)
-       end do
-       !$acc end kernels
+!!$       !$acc kernels present_cr(ut0,zu_st)
+!!$       ! mnh_do_concurrent ( ii=nib:nie , ij=njb:nje , ik=nkb:nke )
+!!$          ut0(ii,ij,ik) = zu_st(ii,ij,ik)
+!!$       ! mnh_end_do()
+!!$       !$acc end kernels
+       call dcopy(n_lin,zu_st,1,ut0,1)   
     end if
     ! Create residual vector
     r => Tjacobi(level)%r
@@ -1592,12 +1594,12 @@ end subroutine construct_vertical_coeff
        zSut0_st => Sut0%st
        zu_st => u%st
 
-       !$acc kernels present(zSut0_st,zu_st)
-       !$acc_nv loop independent collapse(3)
-        do concurrent (ii=nib:nie,ij=njb:nje,ik=nkb:nke)
-           zSut0_st(ii,ij,ik) = zu_st(ii,ij,ik)
-        end do
-       !$acc end kernels
+!!$       !$acc kernels present(zSut0_st,zu_st)
+!!$       ! mnh_do_concurrent( ii=nib:nie , ij=njb:nje , ik=nkb:nke )
+!!$           zSut0_st(ii,ij,ik) = zu_st(ii,ij,ik)
+!!$       ! mnh_end_do()
+!!$       !$acc end kernels
+       call dcopy(n_lin,zu_st,1,zSut0_st,1)    
 
        Sutmp => Tjacobi(level)%Sutmp
        
@@ -1672,13 +1674,12 @@ end subroutine construct_vertical_coeff
        call apply_tridiag_solve_mnh_allT(iib,iie,ijb,ije,Sr,c,b,        &
                   Sut0, &
                   Sutmp,level )
-       !$acc kernels
-       !$acc_nv loop independent collapse(3)
-       do concurrent ( ix=iib:iie,iy=ijb:ije,iz=1:nz )
+       !$acc kernels present_cr(zsut0_st,zu_st)
+       !$mnh_do_concurrent( ix=iib:iie  , iy=ijb:ije , iz=1:nz )
           zu_st(ix,iy,iz) = & 
                rho*zSutmp_st(ix,iy,iz) & 
                + (1.0_rl-rho)*zSut0_st(ix,iy,iz)
-       end do ! concurrent
+       !$mnh_end_do() ! concurrent
        !$acc end kernels
     end if
 
@@ -2081,9 +2082,9 @@ end subroutine line_Jacobi_mnh
       tmp_k => Ttridiag_mnh(level)%tmp_k
       c_k => Ttridiag_mnh(level)%c_k
          
-      !$acc kernels present(zSr_st,zSu_out_st)
-      !$acc_nv loop independent collapse(2)
-      do concurrent(ii=iib:iie,ij=ijb:ije)
+      ! acc kernels present(zSr_st,zSu_out_st)
+      !$acc parallel present(zSr_st,zSu_out_st)
+      !$mnh_do_concurrent ( ii=iib:iie , ij=ijb:ije )   
          iz=1 
          zSr_st(ii,ij,iz) = zb_st(ii,ij,iz)
          !$acc loop seq
@@ -2098,7 +2099,7 @@ end subroutine line_Jacobi_mnh
          zSr_st(ii,ij,iz) = zb_st(ii,ij,iz)
       !
       ! Thomas algorithm
-      ! 
+      !
          iz = 1     
          zSu_out_st(ii,ij,iz) = zSr_st(ii,ij,iz) / (tmp_k(iz)*Tijs*zd_k(iz))
          !$acc loop seq
@@ -2116,8 +2117,10 @@ end subroutine line_Jacobi_mnh
             zSu_out_st(ii,ij,iz) = zSu_out_st(ii,ij,iz) & 
                  - c_k(iz) * zSu_out_st(ii,ij,iz+1)
          end do
-      end do
-      !$acc end kernels
+      !end do; end do
+      !$mnh_end_do()   
+      !$acc end parallel
+      ! acc end kernels
        
     end if
 
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/mg_main_mnh.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/mg_main_mnh.f90
index e20657cfe..b06568415 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/mg_main_mnh.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/mg_main_mnh.f90
@@ -175,7 +175,7 @@ end subroutine mg_main_get_u_mnh
 
     call  mg_main_mnh_init()
 
-    DO it=1,1
+    DO it=1,10
        
        call mg_main_initialise_rhs_mnh()
        call mg_main_initialise_u_mnh()
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/modd_mpif.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/modd_mpif.f90
new file mode 100644
index 000000000..7ad1e6fad
--- /dev/null
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/modd_mpif.f90
@@ -0,0 +1,3 @@
+MODULE MODD_MPIF
+  include 'mpif.h'
+END MODULE MODD_MPIF
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/mode_mg.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/mode_mg.f90
index 11ade04e6..921d0cd3e 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/mode_mg.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/mode_mg.f90
@@ -286,14 +286,15 @@ implicit none
 
   ! Save fields to disk
   if (savefields) then
-    call haloswap(mg_param%n_lev,pproc,xu_fine)
+    call haloswap_mnh(mg_param%n_lev,pproc,xu_fine)
     call save_scalar3d(MPI_COMM_HORIZ,xu_fine,"solution")
-    call volscale_scalar3d(xb_fine,1)
-    call calculate_residual(mg_param%n_lev,pproc,xb_fine,xu_fine,xr_fine)
-    call volscale_scalar3d(xb_fine,-1)
-    call volscale_scalar3d(xr_fine,-1)
-    call haloswap(mg_param%n_lev,pproc,xr_fine)
+    call volscale_scalar3d_mnh(xb_fine,1)
+    call calculate_residual_mnh(mg_param%n_lev,pproc,xb_fine,xu_fine,xr_fine)
+    call volscale_scalar3d_mnh(xb_fine,-1)
+    call volscale_scalar3d_mnh(xr_fine,-1)
+    call haloswap_mnh(mg_param%n_lev,pproc,xr_fine)
     call save_scalar3d(MPI_COMM_HORIZ,xr_fine,"residual")
+    call save_scalar3d(MPI_COMM_HORIZ,xb_fine,"righthand")
   end if
 
   if (i_am_master_mpi) then
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
index 2fd6bc0b0..db1065fe2 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
@@ -510,18 +510,14 @@ contains
       if (LUseT) then 
          zphifine_st => phifine%st
          zphicoarse_st => phicoarse%st
-         !$acc kernels
-#ifdef MNH_COMPILER_NVHPC         
-         !$acc loop independent  collapse(3)
-#endif         
-         do concurrent (ix=ix_min:ix_max, iy=iy_min:iy_max, iz=1:nz)
-                  zphicoarse_st(ix,iy,iz) =  &
+         !$acc kernels present_cr(zphifine_st,zphicoarse_st)
+         !$mnh_do_concurrent( ix=ix_min:ix_max, iy=iy_min:iy_max, iz=1:nz )
+              zphicoarse_st(ix,iy,iz) =  &
                        zphifine_st(2*ix  ,2*iy  ,iz) + &
                        zphifine_st(2*ix  ,2*iy-1,iz) + &
                        zphifine_st(2*ix-1,2*iy  ,iz) + &
                        zphifine_st(2*ix-1,2*iy-1,iz)
-
-         end do
+         !$mnh_end_do()            
          !$acc end kernels
       endif
 
@@ -847,23 +843,25 @@ contains
          zphifine_st => phifine%st
          zphicoarse_st => phicoarse%st
 
-         !$acc kernels
+         !!!!!$acc kernels
 #ifdef MNH_COMPILER_NVHPC         
-         !$acc loop independent  collapse(5)
+         !!!!!$acc loop independent  collapse(5)
 #endif         
+         !$acc parallel loop collapse(5)
          do diy = -1,0
             do dix = -1,0
-               DO CONCURRENT (ix=ix_min:ix_max, iy=iy_min:iy_max, iz=1:nz)
+               !!!!DO CONCURRENT (ix=ix_min:ix_max, iy=iy_min:iy_max, iz=1:nz)
+               do iz=1,nz ; do iy=iy_min,iy_max; do ix=ix_min,ix_max;
                   zphifine_st(2*ix+dix,2*iy+diy,iz) =      &
                        zphicoarse_st(ix,iy,iz) +                &
                        rhox*(zphicoarse_st(ix+(2*dix+1),iy,iz)  &
                        - zphicoarse_st(ix,iy,iz)) +       &
                        rhoy*(zphicoarse_st(ix,iy+(2*diy+1),iz)  &
                        - zphicoarse_st(ix,iy,iz))                  
+               end do; end do; end do
                end do
             end do
-         end do
-         !$acc end kernels
+         !!!!!$acc end kernels
       end if
 
     end subroutine loop_over_grid_linear_mnh
@@ -1205,7 +1203,7 @@ contains
 ! Multigrid V-cycle
 !==================================================================
   recursive subroutine mg_vcycle_mnh(b,u,r,finelevel,splitlevel,level,m)
-#ifdef MNH
+#ifdef MNH_MPPDB
     USE MODE_MPPDB
 #endif
     use MODE_OPENACC_SET_DEVICE
@@ -1232,7 +1230,7 @@ contains
     n_gridpoints = (nlocalx+2*halo_size) &
                  * (nlocaly+2*halo_size) &
                  * (u(level,m)%grid_param%nz+2)
-#ifdef MNH
+#ifdef MNH_MPPDB
     if ( mppdb_initialized ) then
        WRITE(ylevel, '( I3.3 )' ) level
      call Mppdb_check3D_real_mg( b(level,m)%st," mg_vcycle_mnh beg:"//ylevel//":b"       )
@@ -1354,7 +1352,7 @@ contains
       call finish_timer(t_total(level,m))
       call finish_timer(t_coarsesolve(level,m))
     end if
-#ifdef MNH
+#ifdef MNH_MPPDB
     if ( mppdb_initialized ) then
        WRITE(ylevel, '( I3.3 )' ) level
      call Mppdb_check3D_real_mg( b(level,m)%st," mg_vcycle_mnh end:"//ylevel//":b"       )
@@ -1519,7 +1517,7 @@ contains
 ! Assumes that ghosts in initial solution are up-to-date
 !==================================================================
   subroutine mg_solve_mnh(bRHS,usolution,solver_param)
-#ifdef MNH
+#ifdef MNH_MPPDB
     USE MODE_MPPDB
 #endif    
     implicit none
@@ -1739,7 +1737,7 @@ contains
       if (LUseO) call dcopy(n_gridpoints,xu_mg(level,pproc)%s,1,usolution%s,1)
       if (LUseT) call dcopy(n_gridpoints,xu_mg(level,pproc)%st,1,usolution%st,1)
     end if
-#ifdef MNH
+#ifdef MNH_MPPDB
   if ( mppdb_initialized ) then
      call Mppdb_check(      bRHS%st,"mg_solve_mnh:bRHS     "       )
      call Mppdb_check( usolution%st,"mg_solve_mnh:usolution"       )     
-- 
GitLab