diff --git a/src/LIB/SURCOUCHE/src/modd_confz.f90 b/src/LIB/SURCOUCHE/src/modd_confz.f90
index 9712fe2a7bd2dd3c1f1d35c47e02e16a0a106362..0942f667d93d57608fc7bb7290a65940b830ddfe 100644
--- a/src/LIB/SURCOUCHE/src/modd_confz.f90
+++ b/src/LIB/SURCOUCHE/src/modd_confz.f90
@@ -54,6 +54,8 @@ LOGICAL,SAVE      :: LMNH_MPI_ALLTOALLV_REMAP = .FALSE. ! default remap with sen
 INTEGER,SAVE      :: NZ_SPLITTING   = 10 ! /!\ setting of NZ_SPLITTING by namelist for 'EXPERT' use only for DEBUG 
                                          ! 'STANDARD' user use LMNH_MPI_ALLTOALLV_REMAP=T/F only !!!
                                          !  IZ=1=flat_inv;  IZ=2=flat_invz ;  IZ=1+2=the two ; IZ=4 alltoall ; +8=P1/P2 splitting
-!JUAN
+
+INTEGER,SAVE      :: NPMAX_T1DFLAT_R = 130 ! Used to determine max size of buffer ZT1DFLAT in mode_mnh_zwork.f90
+                                           ! (3D size of the mesh * NPMAX_T1DFLAT_R)
 !
 END MODULE MODD_CONFZ
diff --git a/src/LIB/SURCOUCHE/src/modn_confz.f90 b/src/LIB/SURCOUCHE/src/modn_confz.f90
index 3938154448c85d8a002d5af0c1485a3323c99b10..232f4e99d273ef4e56564cf978de131f154a31b8 100644
--- a/src/LIB/SURCOUCHE/src/modn_confz.f90
+++ b/src/LIB/SURCOUCHE/src/modn_confz.f90
@@ -43,6 +43,7 @@ USE MODD_CONFZ
 IMPLICIT NONE
 !
 NAMELIST/NAM_CONFZ/ NZ_VERB,NZ_PROC,NB_PROCIO_R,NB_PROCIO_W,MPI_BUFFER_SIZE,LMNH_MPI_BSEND & 
-                   ,LMNH_MPI_ALLTOALLV_REMAP,NZ_SPLITTING !JUAN Z_SPLITTING
+                   ,LMNH_MPI_ALLTOALLV_REMAP,NZ_SPLITTING &
+                   ,NPMAX_T1DFLAT_R 
 !
 END MODULE MODN_CONFZ
diff --git a/src/MNH/mode_mnh_zwork.f90 b/src/MNH/mode_mnh_zwork.f90
index 38b5e1c721bd97668f198c6f18be58f1d0eeb505..73fb5b2f0b54d8f0acbde36d9b155878822c0eb8 100644
--- a/src/MNH/mode_mnh_zwork.f90
+++ b/src/MNH/mode_mnh_zwork.f90
@@ -17,6 +17,8 @@ MODULE MODE_MNH_ZWORK
 
   use mode_msg
 
+  USE  MODD_CONFZ, ONLY : NPMAX_T1DFLAT_R
+
   IMPLICIT NONE
 
   INTEGER, SAVE :: IIB,IJB,IKB  ! Begining useful area in x,y,z directions
@@ -40,8 +42,6 @@ MODULE MODE_MNH_ZWORK
   INTEGER , ALLOCATABLE, DIMENSION (:)      :: NT3D_POOL
   INTEGER                                   :: NT3D_TOP , NT3D_TOP_MAX = 0
   INTEGER                                   :: NT3D_TOP_CURRENT(JPMAX_T3D+1) ,  NT3D_TOP_CURRENT_INDEX = 0
-  !REAL    , ALLOCATABLE, DIMENSION(:,:,:,:) , TARGET :: ZT3D_A1,ZT3D_A2,ZT3D_A3,ZT3D_A4
-  !REAL    , POINTER    , DIMENSION(:,:,:,:)          :: ZT3D
   REAL,SAVE    , ALLOCATABLE, TARGET , DIMENSION(:,:,:,:)          :: ZT3D
 
   REAL,SAVE    , ALLOCATABLE, TARGET , DIMENSION(:)                :: ZT1D_OSIZE
@@ -149,8 +149,8 @@ MODULE MODE_MNH_ZWORK
 
 
 !------ Real 1DFLAT pool
-  INTEGER, PARAMETER                                 :: JPMAX_T1DFLAT_R = 300      !Used to determine max size of buffer ZT1DFLAT
-                                                                                   !(3D size of the mesh * JPMAX_T1DFLAT_R)
+!!$  INTEGER,                SAVE                       :: NPMAX_T1DFLAT_R = 300      !Used to determine max size of buffer ZT1DFLAT
+!!$                                                                                   !(3D size of the mesh * NPMAX_T1DFLAT_R)
   INTEGER,                SAVE                       :: NPMAX_POOL_T1DFLAT_R = 250 !Maximum size of the pool (max number of arrays)
   INTEGER(KIND=MNHINT64), ALLOCATABLE, DIMENSION (:) :: NT1DFLAT_POOL_R   !Position in ZT1DFLAT of the beginning of each array
   INTEGER(KIND=MNHINT64), ALLOCATABLE, DIMENSION (:) :: NT1DFLAT_SIZE_R   !Size of each array
@@ -347,7 +347,7 @@ CONTAINS
 
 !------ Real 1DFLAT pool
 
-       NT1DFLAT_MAXSIZE_R = INT( IIU, KIND=MNHINT64 ) * IJU * IKU * JPMAX_T1DFLAT_R
+       NT1DFLAT_MAXSIZE_R = INT( IIU, KIND=MNHINT64 ) * IJU * IKU * NPMAX_T1DFLAT_R
        ALLOCATE( ZT1DFLAT(NT1DFLAT_MAXSIZE_R) )
        !$acc enter data create( ZT1DFLAT )
 
@@ -1539,7 +1539,9 @@ CONTAINS
     NTOT_GETSIZE_ZT1DFLAT  = NTOT_GETSIZE_ZT1DFLAT + KSIZE
 
     IF ( NT1DFLAT_POS_R + KSIZE > NT1DFLAT_MAXSIZE_R ) THEN
-      print*,"MNH_GET_ZT1DFLAT ZT1DFLAT too small, JPMAX_T1DFLAT_R =" , JPMAX_T1DFLAT_R
+      print*,"MNH_GET_ZT1DFLAT:: ZT1DFLAT too small, NPMAX_T1DFLAT_R =" , NPMAX_T1DFLAT_R
+      print*,"MNH_GET_ZT1DFLAT:: augmente NPMAX_T1DFLAT_R in NAM_CONFZ"
+      flush(6)
       WRITE( YSIZE,  '( I0 )' ) KSIZE
       WRITE( YAVAIL, '( I0 )' ) NT1DFLAT_MAXSIZE_R - NT1DFLAT_POS_R
       WRITE( YMAX,   '( I0 )' ) NT1DFLAT_MAXSIZE_R
diff --git a/src/MNH/slow_terms.f90 b/src/MNH/slow_terms.f90
index c962c71122cd656eb0f75ee0cbd7d151e941897e..d26ffcbdf857de3655b1cf3445201cba5e3ba5c2 100644
--- a/src/MNH/slow_terms.f90
+++ b/src/MNH/slow_terms.f90
@@ -166,9 +166,16 @@ USE MODD_PARAMETERS, only: JPVEXT
 use mode_budget,     only: Budget_store_init, Budget_store_end
 use mode_mppdb
 
+#ifdef MNH_OPENACC
+  USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
+#endif
 #if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
 use modi_bitrep
 #endif
+#if defined(MNH_COMPILER_CCE) && defined(MNH_BITREP_OMP)
+!$mnh_undef(LOOP)
+!$mnh_undef(OPENACC)
+#endif
 
 IMPLICIT NONE
 !
@@ -213,19 +220,21 @@ INTEGER :: IKE           !  the microphysical sources have to be computed
 !
 REAL    :: ZTSPLITR      ! Small time step for rain sedimentation
 !
-REAL,    DIMENSION(:,:,:), ALLOCATABLE :: ZT,ZW,ZW1,ZW2,ZW3,ZEXNT,ZDZZ  ! Work arrays
-LOGICAL, DIMENSION(:,:,:), ALLOCATABLE :: G3D
+REAL,    DIMENSION(:,:,:), POINTER,CONTIGUOUS :: ZT,ZW,ZW1,ZW2,ZW3,ZEXNT,ZDZZ  ! Work arrays
+LOGICAL, DIMENSION(:,:,:), POINTER,CONTIGUOUS :: G3D
 !
 INTEGER                              :: JI,JJ,IC,JL ! loop control for packed array
+INTEGER                              :: JIU, JJU, JKU
 
 !-------------------------------------------------------------------------------
 
+JIU = size(pzz,1)
+JJU = size(pzz,2)
+JKU = size(pzz,3)
+
 !$acc data present( PZZ, PRHODJ, PRHODREF, PCLDFR, PTHT, PRVT, PRCT, PRRT,      &
 !$acc &             PPABST, PTHS, PRVS, PRCS, PRRS, PINPRR, PINPRR3D, PEVAP3D )
 
-! !$acc &     copyin( XC1RC, XC2RC, XC1RE, XC2RE, XCEXRA, XCEXRE, XCEXRS, XCEXVT, XCRA, XCRS, XDIVA, XTHCO, &
-! !$acc &             XALPW, XBETAW, XGAMW, XCL, XCPD, XCPV, XLVTT, XMD, XMV, XP00, XRD, XRHOLW, XRV, XTT )
-
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK(PZZ,     "SLOW_TERMS beg:PZZ")
@@ -247,16 +256,29 @@ IF (MPPDB_INITIALIZED) THEN
   CALL MPPDB_CHECK(PEVAP3D, "SLOW_TERMS beg:PEVAP3D")
 END IF
 
-allocate( zt   ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( zw   ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( zw1  ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( zw2  ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( zw3  ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( zexnt( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( zdzz ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
-allocate( g3d   ( size(prhodj,1), size(prhodj,2), size(prhodj,3) ) )
+#ifndef MNH_OPENACC
+allocate( zt   ,JIU,JJU,JKU )
+allocate( zw   ,JIU,JJU,JKU )
+allocate( zw1  ,JIU,JJU,JKU )
+allocate( zw2  ,JIU,JJU,JKU )
+allocate( zw3  ,JIU,JJU,JKU )
+allocate( zexnt,JIU,JJU,JKU )
+allocate( zdzz ,JIU,JJU,JKU )
+allocate( g3d   )
+#else
+!Pin positions in the pools of MNH memory
+CALL MNH_MEM_POSITION_PIN()
+CALL MNH_MEM_GET( zt   ,JIU,JJU,JKU )
+CALL MNH_MEM_GET( zw   ,JIU,JJU,JKU )
+CALL MNH_MEM_GET( zw1  ,JIU,JJU,JKU )
+CALL MNH_MEM_GET( zw2  ,JIU,JJU,JKU )
+CALL MNH_MEM_GET( zw3  ,JIU,JJU,JKU )
+CALL MNH_MEM_GET( zexnt,JIU,JJU,JKU )
+CALL MNH_MEM_GET( zdzz ,JIU,JJU,JKU )
+CALL MNH_MEM_GET( g3d  ,JIU,JJU,JKU )
+#endif
 
-!$acc data create( zt, zw, zw1, zw2, zw3, zexnt, zdzz, g3d )
+!$acc data present( zt, zw, zw1, zw2, zw3, zexnt, zdzz, g3d )
 !
 !*       1.     COMPUTE THE LOOP BOUNDS AND EXNER FUNCTION
 !   	        ------------------------------------------
@@ -279,33 +301,24 @@ if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'SEDI', prrs(:,
 !
 !*       2.1    time splitting loop initialization        
 !
-!$acc kernels
 ZTSPLITR = PTSTEP / REAL(KSPLITR)       ! Small time step
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!$acc kernels
 ! 
 ZW1(:,:,:) = PRRS(:,:,:) * PTSTEP
 ZW2(:,:,:) = 0.
 ZW3(:,:,:) = 0.
 !
-#ifndef MNH_OPENACC
 G3D(:,:,:)=.FALSE.
 DO JK=IKE,1,-1
+  !$mnh_expand_where(JI=1:JIU,JJ=1:JJU) 
   WHERE( ZW1(:,:,JK)>0. .OR. G3D(:,:,JK+1) )
     G3D(:,:,JK)=.TRUE.
   END WHERE
+  !$mnh_end_expand_where()
 END DO
-#else
-g3d(:, :, : ) = .false.
-do jk = ike, 1, -1
-  do jj = 1, size(pzz,2)
-    do ji = 1, size(pzz,1)
-      if ( zw1(ji, jj, jk ) >0. .or. g3d(ji, jj, jk+1 ) ) then
-        g3d(ji, jj, jk ) = .true.
-      end if
-    end do
-  end do
-end do
-#endif
 !
+!$mnh_expand_where(JI=1:JIU,JJ=1:JJU,JK=1:JKU)
 WHERE (G3D(:,:,:))
 #if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
   ZW3(:,:,:) = PRHODREF(:,:,:) ** (XCEXRS-XCEXVT)
@@ -313,7 +326,9 @@ WHERE (G3D(:,:,:))
   ZW3(:,:,:) = BR_POW( PRHODREF(:,:,:), XCEXRS - XCEXVT )
 #endif
 END WHERE
+!$mnh_end_expand_where()
 !
+!$mnh_expand_where(JI=1:JIU,JJ=1:JJU)
 WHERE (ZW1(:,:,IKE+1)>0.)
 #if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
   ZW2(:,:,IKE+1) =   XCRS                         &
@@ -325,6 +340,7 @@ WHERE (ZW1(:,:,IKE+1)>0.)
                   * BR_POW( PRHODREF(:,:,IKE+1), XCEXRS - XCEXVT )
 #endif
 END WHERE
+!$mnh_end_expand_where()
 !
 !
 !*       2.2    small time step integration
@@ -337,9 +353,7 @@ DO JN=1,KSPLITR
 !
 
   DO JK=IKE,IKB,-1
-!$acc loop collapse(2) independent
-    DO JJ = 1,SIZE( ZW1,2)
-      DO JI = 1,SIZE( ZW1,1)
+    !$mnh_do_concurrent(JI=1:JIU,JJ=1:JJU)
         IF ( ( ZW1(JI,JJ,JK+1)>0. ) .AND. ( ZW1(JI,JJ,JK)>0. ) )  THEN
 !
 !*       2.2.2    compute the rain flux
@@ -360,9 +374,8 @@ DO JN=1,KSPLITR
                     ( ZW2(JI,JJ,JK+1)-ZW2(JI,JJ,JK) ) /    &
                     ( ZDZZ(JI,JJ,JK) * PRHODREF(JI,JJ,JK) )
         END IF
-      END DO
-    ENDDO
-  ENDDO
+    !$mnh_end_do()
+  END DO     
 !
 !*       2.2.4    compute the explicit accumulated precipitations
 !                 -----------------------------------------------
@@ -378,6 +391,7 @@ END DO
 !
 PRRS(:,:,:) = ZW1(:,:,:) / PTSTEP
 !$acc end kernels
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
 !*       2.5     budget storage
 !
@@ -396,6 +410,7 @@ if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'ACCR', prrs(:,
 !
 !$acc kernels
 G3D(:,:,:) = PRCT(:,:,:)>0.0 .AND. PRRT(:,:,:)>0.0 .AND. PRCS(:,:,:)>0.0
+!$mnh_expand_where(JI=1:JIU,JJ=1:JJU,JK=1:JKU)
 WHERE ( G3D(:,:,:) )
 #if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
   ZW(:,:,:) = XCRA * PRCT(:,:,:)                        &
@@ -411,6 +426,7 @@ WHERE ( G3D(:,:,:) )
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
 END WHERE
+!$mnh_end_expand_where()
 !$acc end kernels
 !
 !*       3.2     budget storage
@@ -432,14 +448,17 @@ if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'AUTO', prrs(:,
 !$acc kernels
 IF ( HSUBG_AUCV == 'CLFR' ) THEN
  G3D(:,:,:) = PRCT(:,:,:)>0.0 .AND. PRCS(:,:,:)>0.0 .AND. PCLDFR(:,:,:)>0.0
+ !$mnh_expand_where(JI=1:JIU,JJ=1:JJU,JK=1:JKU)
  WHERE ( G3D(:,:,:) )
   ZW(:,:,:) = XC1RC * MAX(PRCT(:,:,:) /(PCLDFR(:,:,:)) - XC2RC / PRHODREF(:,:,:),0.)
   ZW(:,:,:) = MIN ( (ZW(:,:,:)* PCLDFR(:,:,:)), PRCS(:,:,:) )
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
  END WHERE
+ !$mnh_end_expand_where()
 ELSE
  G3D(:,:,:) = PRCT(:,:,:)>0.0 .AND. PRCS(:,:,:)>0.0
+ !$mnh_expand_where(JI=1:JIU,JJ=1:JJU,JK=1:JKU)
  WHERE ( G3D(:,:,:) )
   ZW(:,:,:) = XC1RC * MAX ( PRCT(:,:,:) - XC2RC / PRHODREF(:,:,:), 0. )
   ZW(:,:,:) = MIN ( ZW(:,:,:) , PRCS(:,:,:) )
@@ -447,6 +466,7 @@ ELSE
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
  END WHERE
+ !$mnh_end_expand_where()
 END IF
 !$acc end kernels
 !
@@ -467,6 +487,7 @@ if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'REVA', prrs(:,
 !$acc kernels
 PEVAP3D(:,:,:)=0.
 G3D(:,:,:) = PRRT(:,:,:)>0.0 .AND. PRCT(:,:,:)==0.0
+!$mnh_expand_where(JI=1:JIU,JJ=1:JJU,JK=1:JKU)
 WHERE ( G3D(:,:,:) )
 !
 !*       5.1    compute the Exner function
@@ -531,6 +552,7 @@ WHERE ( G3D(:,:,:) )
                       + XCL * (PRCT(:,:,:) + PRRT(:,:,:)) ) )
   PEVAP3D(:,:,:)=ZW(:,:,:)
 END WHERE
+!$mnh_end_expand_where()
 !$acc end kernels
 !
 !*       5.8     budget storage
@@ -565,6 +587,13 @@ END IF
 
 !$acc end data
 
+#ifndef MNH_OPENACC
+deallocate (zt, zw, zw1, zw2, zw3, zexnt, zdzz, g3d )
+#else
+!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
+CALL MNH_MEM_RELEASE()
+#endif
+
 !$acc end data
 
 END SUBROUTINE SLOW_TERMS
diff --git a/src/Makefile b/src/Makefile
index ef8e3f15572713f0af1aa4f09e74053e70f58b58..ee69febd6d0ad21fe7075d4e0278e34ea5c4dfa6 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -17,14 +17,6 @@ obj_flist        = $(notdir $(shell find $(1) -follow -type f \
                -name 'spll_*.f*' | sed -e 's/\(.*\)\(\.\).*/\1.o/g' )) 
 #
 ##########################################################
-#                                                        #
-# COMPILER & ARCHITECTURE CONFIGURATION                  #
-#                                                        #
-##########################################################
-#
-include Rules.$(ARCH)$(F).mk
-#
-##########################################################
 ##########################################################
 ##########################################################
 #                                                        #
@@ -115,10 +107,17 @@ endif
 endif
 ##########################################################
 #                                                        #
-# VPATH ADJUSTEMENT :                                    #
+# COMPILER & ARCHITECTURE CONFIGURATION                  #
 #                                                        #
 ##########################################################
 #
+include Rules.$(ARCH)$(F).mk
+#
+##########################################################
+#                                                        #
+# VPATH ADJUSTEMENT :                                    #
+#                                                        #
+##########################################################
 #
 #VPATH := $(filter-out $(VPATH_EXCLUDE),$(VPATH))
 #
diff --git a/src/Rules.LXcray.mk b/src/Rules.LXcray.mk
index d21f8318df6236687ea32d9e8aec8a5f072bf2d1..5fdd9a8d4f93384bbcfd84e11b4122cb243e807d 100644
--- a/src/Rules.LXcray.mk
+++ b/src/Rules.LXcray.mk
@@ -168,7 +168,8 @@ FX90FLAGS =  $(OPT)
 # -132 
 #
 #LDFLAGS    =  -Wl,-noinhibit-exec  -Wl,-warn-once $(PAR)
-LDFLAGS    =   -Wl,-warn-once $(PAR) $(OPT_BASE)
+LDFLAGS    =  -Wl,--allow-multiple-definition -Wl,-warn-once $(PAR) $(OPT_BASE)
+# -hsystem_alloc
 #
 # preprocessing flags 
 #
diff --git a/src/Rules.LXnvhpc2202.mk b/src/Rules.LXnvhpc2202.mk
index 15ef6ced160560b3f1ca44055813aef81adb69cc..341f0d2789b401a44806d27d1c60a949f31ba440 100644
--- a/src/Rules.LXnvhpc2202.mk
+++ b/src/Rules.LXnvhpc2202.mk
@@ -18,17 +18,18 @@
 #TP= -tp=sandybridge
 #TP= -tp=skylake
 #TP= -tp=p7-64
-TP= -tp=px
+#TP= -tp=px
+TP ?= -tp=native
 #
 #Version of CUDA
 #(8.0 at least if compute capability >= 6.0)
-CUDALEVEL=cuda11.0
+CUDALEVEL ?= cuda11.0
 #
 #Compute capability of GPU
 #
 #OPT_CPTCAP=cc35,cc50,cc70
 #OPT_CPTCAP=cc35,cc50,cc60,cc70
-OPT_CPTCAP=cc70
+OPT_CPTCAP ?= cc70
 #Aeropc45: cc50
 #Nuwa: cc35
 #Ouessant Firestone K80: cc35
@@ -90,7 +91,7 @@ OPT_NOCB  = $(OPT_BASE) $(OPT_PERF2)
 #
 # List of Files with compilation problem in O2
 #
-OBJS_REPROD= spll_mode_sum_ll.o mode_sum_ll.mod
+OBJS_REPROD= spll_mode_sum_ll.o mode_sum_ll.mod  spll_mode_repro_sum.o mode_repro_sum.mod
 
 OBJS_O1_OPENACC= spll_ice4_tendencies.o spll_turb_ver_thermo_flux.o spll_mode_device.o mode_device.mod spll_mppdb_check3d_real_mg.o \
                  spll_mppdb_check0d_real_mg.o spll_mode_turb.o mode_turb.mod spll_modd_les_n.o modd_les_n.mod
@@ -289,7 +290,6 @@ $(OBJS_O1) : OPT = $(OPT_BASE) $(OPT_PERF1)
 
 OBJS_O0= spll_mode_mppdb.o mode_mppdb.mod \
          spll_fft55.o spll_fft.o spll_flat_invz.o \
-         spll_mode_repro_sum.o mode_repro_sum.mod \
          spll_modd_les_n.o modd_les_n.mod \
          spll_default_desfm_n.o \
          spll_modd_pack_gr_field_n.o modd_pack_gr_field_n.mod
diff --git a/src/ZSOLVER/dyn_sources.f90 b/src/ZSOLVER/dyn_sources.f90
new file mode 100644
index 0000000000000000000000000000000000000000..aa8feb2fc99f4836c19fd04492f442826b0457a1
--- /dev/null
+++ b/src/ZSOLVER/dyn_sources.f90
@@ -0,0 +1,607 @@
+!MNH_LIC Copyright 1994-2021 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.
+!-----------------------------------------------------------------
+!     #######################
+      MODULE MODI_DYN_SOURCES
+!     #######################
+!
+INTERFACE
+!
+      SUBROUTINE DYN_SOURCES( KRR,KRRL, KRRI,                              &
+                              PUT, PVT, PWT, PTHT, PRT,                    &
+                              PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY,   &
+                              PRHODJ, PZZ, PTHVREF, PEXNREF,               &
+                              PRUS, PRVS, PRWS, PRTHS                      )
+!
+INTEGER,                  INTENT(IN)    :: KRR  ! Total number of water var.
+INTEGER,                  INTENT(IN)    :: KRRL ! Number of liquid water var.
+INTEGER,                  INTENT(IN)    :: KRRI ! Number of ice water var.
+!
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT, PVT, PWT   ! variables
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT            !     at
+REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT             !  time t
+!
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCORIOX   !    2D horizontal
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCORIOY   !      Coriolis 
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCORIOZ   !       factors
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCURVX    !    2D horizontal
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCURVY    !  curvature factors
+!
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ    ! Density * Jacobian
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ       ! Height (z)
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHVREF   ! Virtual Temperature
+                                          ! of the reference state
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PEXNREF   ! Exner function
+                                          ! of the reference state
+!
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRUS, PRVS, PRWS ! Sources of Momentum
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS            ! Sources of theta
+!
+END SUBROUTINE DYN_SOURCES
+!
+END INTERFACE
+!
+END MODULE MODI_DYN_SOURCES 
+!     ######################################################################
+      SUBROUTINE DYN_SOURCES( KRR,KRRL, KRRI,                              &
+                              PUT, PVT, PWT, PTHT, PRT,                    &
+                              PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY,   &
+                              PRHODJ, PZZ, PTHVREF, PEXNREF,               &
+                              PRUS, PRVS, PRWS, PRTHS                      )
+!     ######################################################################
+!
+!!****  *DYN_SOURCES * - routine to compute the curvature, coriolis and gravity terms
+!!
+!!    PURPOSE
+!!    -------
+!!      The purpose of this routine is to compute the dynamical sources
+!!    (curvature, coriolis and gravity terms) for each component of the
+!!    momentum. The curvature terms arise in case of non-cartesian geometry
+!!    only. For both curvature and coriolis terms, the "thin shell" 
+!!    approximation can be assumed or not. Gravity affects only the vertical
+!!    component of the momentum.
+!!      This routine also adds the Lagrangien derivative of the pressure as a
+!!    source for the theta evolution.
+!!
+!!     
+!!**  METHOD
+!!    ------
+!!      The horizontal curvature and coriolis field arrays are expanded in
+!!    the vertical direction with the SPREAD intrinsics to match with the
+!!    Shuman operators. The source term for theta due to the advection of the
+!!    Exner function is treated in an advective form (or non-conservative form).
+!!    Values of the calorific capacity of the melting and of the potential
+!!    virtual temperature are recovered taking care of the number of water
+!!    variables. The different sources terms are stored for the budget
+!!    computations.
+!!
+!!    EXTERNAL
+!!    --------
+!!      MXM,MYM,MZM : Shuman functions (mean operators)
+!!      MXF,MYF,MZF : Shuman functions (mean operators)
+!!      GZ_M_W      : projection along the vertical direction of the gradient
+!!                    vector. It acts on a field localized in mass point and
+!!                    the result is localized in w point.
+!!      BUDGET      : Stores the different budget components
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_CONF     : contains configuration variables for all models
+!!           LTHINSHELL  = .TRUE.  if the THINSHELL approximation is made
+!!           LCARTESIAN  = .TRUE.  if the CARTESIAN approximation is made
+!!      Module MODD_CST      : 
+!!           XG           gravity acceleration
+!!           XRADIUS      Earth radius
+!!           XRD,XRV      Gas constant for dry air and wator vapor
+!!           XCPD,XCPV    Cp for dry air and wator vapor
+!!           XCL,XCI      C (calorific capacity) for liquid and solid water
+!!      Module MODD_DYN      : contains the parameters for the dynamics
+!!           LCORIO      = .FALSE. if the earth rotation is neglected
+!!    
+!!      Module MODD_BUDGET:
+!!         NBUMOD       : model in which budget is calculated
+!!         CBUTYPE      : type of desired budget
+!!                          'CART' for cartesian box configuration
+!!                          'MASK' for budget zone defined by a mask 
+!!                          'NONE'  ' for no budget
+!!         LBU_RTH      : logical for budget of RTH (potential temperature)
+!!                        .TRUE. = budget of RTH        
+!!                        .FALSE. = no budget of RTH
+!!         LBU_RU       : logical for budget of RU (wind component along x)
+!!                        .TRUE. = budget of RU         
+!!                        .FALSE. = no budget of RU 
+!!         LBU_RV       : logical for budget of RV (wind component along y)
+!!                        .TRUE. = budget of RV         
+!!                        .FALSE. = no budget of RV 
+!!         LBU_RW        : logical for budget of RW (wind component along z)
+!!                        .TRUE. = budget of RW         
+!!                        .FALSE. = no budget of RW 
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation ( routine DYN_SOURCE )
+!!
+!!    AUTHOR
+!!    ------
+!!	J.-P. Pinty      * Laboratoire d'Aerologie*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    06/06/94 
+!!	Corrections 06/08/94 (J.-P. Lafore) 
+!!	Corrections 17/10/94 (Stein) For LCORIO
+!!	Corrections 22/12/94 (Stein) add the pressure term in the theta evolution
+!!  Corrections 30/12/94 (J.P. Lafore) bug corrections for the pressure term 
+!!  Corrections 16/03/95 (Stein) remove R from the historical variables and
+!!                                   correction of the pressure term   
+!!  Corrections 14/04/95 (Masson Stein) bugs in the curvatureand Coriolis
+!!                                   terms for PRUS,PRVS, PRWS         
+!!  Corrections 01/04/95 (Ph. Hereil J. Nicolau) add the budget computation
+!!  Corrections 16/10/95 (J. Stein)     change the budget calls 
+!!  Corrections 16/02/96 (J.-P. Pinty)  Introduce L1D switches
+!!  Corrections 19/12/96 (J.-P. Pinty)  Update the CALL BUDGET
+!!  Corrections 03/12/02  (P. Jabouille)  add no thinshell condition
+!!  Correction     06/10  (C.Lac) Exclude L1D for Coriolis term 
+!!  Modification   03/11  (C.Lac) Split the gravity term due to buoyancy
+!!   J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!  P. Wautelet    02/2020: use the new data structures and subroutines for budgets
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+use modd_budget,      only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, &
+                            NBUDGET_U, NBUDGET_V, NBUDGET_W, NBUDGET_TH, &
+                            tbudgets
+USE MODD_CONF
+USE MODD_CST
+USE MODD_DYN
+USE MODD_DYN_n,      ONLY: LOCEAN
+!
+use mode_budget,     only: Budget_store_init, Budget_store_end
+USE MODE_MPPDB
+
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN
+#ifdef MNH_OPENACC
+USE MODI_SHUMAN_DEVICE
+#endif
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
+#endif
+
+
+IMPLICIT NONE
+!  
+!*       0.1   Declarations of dummy arguments :
+!
+INTEGER,                  INTENT(IN)    :: KRR  ! Total number of water var.
+INTEGER,                  INTENT(IN)    :: KRRL ! Number of liquid water var.
+INTEGER,                  INTENT(IN)    :: KRRI ! Number of ice water var.
+!
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT, PVT, PWT   ! variables
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT            !     at
+REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT             !  time t
+!
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCORIOX   !    2D horizontal
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCORIOY   !      Coriolis 
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCORIOZ   !       factors
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCURVX    !    2D horizontal
+REAL, DIMENSION(:,:),     INTENT(IN)    :: PCURVY    !  curvature factors
+!
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ    ! dry Density * Jacobian
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ       ! Height (z)
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHVREF   ! Virtual Temperature
+                                          ! of the reference state
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PEXNREF   ! Exner function
+                                          ! of the reference state
+!
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRUS, PRVS, PRWS ! Sources of Momentum
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS            ! Sources of theta
+!
+!
+!*       0.2   Declarations of local variables :
+!
+REAL       ::  ZCPD_OV_RD   ! = CPD / RD
+REAL       ::  ZG_OV_CPD    ! =-XG / XCPD
+INTEGER    ::  JWATER       ! loop index on the different types of water
+INTEGER    ::  JIU, JJU, JKU     
+INTEGER    ::  JI,JJ,JK     
+INTEGER    ::  IKU          ! array size along the k direction 
+REAL       ::  ZD1          ! DELTA1 (switch 0/1) for thinshell approximation
+#ifndef MNH_OPENACC
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) ::           &
+                              ZWORK1, ZWORK2, ZWORK3, ZRUT, ZRVT
+#else
+REAL, DIMENSION(:,:,:), pointer , contiguous ::ZWORK1, ZWORK2, ZWORK3, ZRUT, ZRVT
+#endif
+#ifdef MNH_OPENACC
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE,ZTMP2_DEVICE,ZTMP3_DEVICE,ZTMP4_DEVICE,ZRVTT,ZRUTT
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP5_DEVICE,ZTMP6_DEVICE,ZTMP7_DEVICE,ZTMP8_DEVICE
+!!!REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP9_DEVICE,ZTMP10_DEVICE,ZTMP11_DEVICE,ZTMP12_DEVICE
+#endif                      
+!
+!-------------------------------------------------------------------------------
+!
+
+JIU = size( PUT, 1 )
+JJU = size( PUT, 2 )
+JKU = size( PUT, 3 )
+
+#ifdef MNH_OPENACC
+CALL MNH_MEM_POSITION_PIN()
+
+CALL MNH_MEM_GET( ztmp1_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp2_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp3_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp4_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp5_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp6_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp7_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ztmp8_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp9_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp10_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp11_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp12_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRVTT, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRUTT, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZWORK1, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZWORK2, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZWORK3, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRUT, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRVT, JIU, JJU, JKU )
+
+#endif
+
+!
+!*       1.     COMPUTES THE TRUE VELOCITY COMPONENTS
+!	        -------------------------------------
+!
+
+!$acc data copyin             ( PUT, PVT, PWT,PTHT, PRT,                     &
+!$acc &                        PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY,   &
+!$acc &                        PZZ, PTHVREF, PEXNREF)               &       
+!$acc &   present               (PRUS,PRVS, PRWS, PRTHS,PRHODJ)        
+!$acc update device(PRUS, PRVS, PRWS,PRTHS,PRHODJ)
+!$acc update device  ( PUT, PVT, PWT,PTHT, PRT,                     &
+!$acc &                        PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY,   &
+!$acc &                        PZZ, PTHVREF, PEXNREF)
+
+
+#ifndef MNH_OPENACC
+ZRUT(:,:,:) = PUT(:,:,:) * MXM(PRHODJ(:,:,:))
+ZRVT(:,:,:) = PVT(:,:,:) * MYM(PRHODJ(:,:,:))
+
+print *,"PUT(1,1,1)",PUT(1,1,1)
+#else
+CALL MXM_DEVICE(PRHODJ(:,:,:), ZTMP1_DEVICE (:,:,:))
+CALL MYM_DEVICE(PRHODJ(:,:,:), ZTMP2_DEVICE (:,:,:))
+!$acc kernels present(ZRUT,ZRVT,PUT,PVT,ZTMP1_DEVICE,ZTMP2_DEVICE)
+ZRUT(:,:,:) = PUT(:,:,:) * ZTMP1_DEVICE(:,:,:)
+ZRVT(:,:,:) = PVT(:,:,:) * ZTMP2_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+!
+IKU = SIZE(PUT,3)
+!
+!-------------------------------------------------------------------------------
+!
+!*       2.     COMPUTES THE CURVATURE TERMS
+!    	        ----------------------------
+!
+! Only when earth rotation is considered but not in 1D and CARTESIAN cases
+!
+IF ((.NOT.L1D).AND.(.NOT.LCARTESIAN) )  THEN
+  if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'CURV', prus(:, :, :) )
+  if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'CURV', prvs(:, :, :) )
+
+  IF ( LTHINSHELL ) THEN           !  THINSHELL approximation
+#ifndef  MNH_OPENACC
+    ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=JKU ) / XRADIUS
+    ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=JKU ) / XRADIUS
+#else
+ !$acc kernels present_cr(ZWORK1,ZWORK2)
+  !$mnh_expand_array(JI=1:JIU  , JJ=1:JJU , JK=1:JKU)   
+    ZWORK1(:,:,:) = PCURVX(:,:) / XRADIUS
+    ZWORK2(:,:,:) = PCURVY(:,:) / XRADIUS
+  !$mnh_end_expand_array()
+!$acc end kernels
+#endif    
+#ifndef MNH_OPENACC
+    PRUS(:,:,:) = PRUS                                   &
+    + ZRUT * MXM(  MYF(PVT) * ZWORK1 )                   &
+    + MXM( MYF(ZRVT*PVT) * ZWORK2 )
+!
+    PRVS(:,:,:) = PRVS                                   &
+    - MYM( MXF(ZRUT*PUT) * ZWORK1 )                      &
+    - ZRVT * MYM( MXF(PUT) * ZWORK2 ) 
+!
+#else
+!$acc kernels
+ZRVTT(:,:,:)=ZRVT(:,:,:)*PVT(:,:,:)
+ZRUTT(:,:,:)=ZRUT(:,:,:)*PUT(:,:,:)
+!$acc end kernels
+CALL MYF_DEVICE(PVT,ZTMP1_DEVICE )   !MYF(PVT)
+CALL MYF_DEVICE(ZRVTT,ZTMP2_DEVICE ) !MYF(ZRVT*PVT)
+CALL MXF_DEVICE(ZRUTT,ZTMP3_DEVICE ) !MXF(ZRUT*PUT)
+CALL MXF_DEVICE(PUT,ZTMP4_DEVICE )   !MXF(PUT)
+
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP1_DEVICE,ZTMP1_DEVICE,ZTMP1_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZTMP1_DEVICE(:,:,:)* ZWORK1(:,:,:)   !MYF(PVT) * ZWORK1
+ZTMP2_DEVICE(:,:,:)=ZTMP2_DEVICE(:,:,:)* ZWORK2(:,:,:)   !MYF(ZRVT*PVT) * ZWORK2
+ZTMP3_DEVICE(:,:,:)=ZTMP3_DEVICE(:,:,:)* ZWORK1(:,:,:)   !MXF(ZRUT*PUT) * ZWORK1
+ZTMP4_DEVICE(:,:,:)=ZTMP4_DEVICE(:,:,:) * ZWORK2(:,:,:)  !MXF(PUT) * ZWORK2
+!$acc end kernels
+
+CALL MXM_DEVICE(ZTMP1_DEVICE,ZTMP5_DEVICE)   !MXM(  MYF(PVT) * ZWORK1 )
+CALL MXM_DEVICE(ZTMP2_DEVICE,ZTMP6_DEVICE)   !MXM( MYF(ZRVT*PVT) * ZWORK2 )
+
+CALL MYM_DEVICE(ZTMP3_DEVICE,ZTMP7_DEVICE)  !MYM( MXF(ZRUT*PUT) * ZWORK1 )
+CALL MYM_DEVICE(ZTMP4_DEVICE,ZTMP8_DEVICE) !MYM( MXF(PUT) * ZWORK2 )
+
+!$acc kernels present_cr(PRUS,PRVS)
+PRUS(:,:,:) = PRUS(:,:,:)                                  &
++ ZRUT(:,:,:) * ZTMP5_DEVICE(:,:,:)                        &
+    + ZTMP6_DEVICE(:,:,:)
+
+PRVS(:,:,:) = PRVS(:,:,:)                                   &
+    - ZTMP7_DEVICE(:,:,:)                                  &
+    - ZRVT(:,:,:) * ZTMP8_DEVICE(:,:,:)
+
+!$acc end kernels
+#endif
+
+  ELSE                           !  NO THINSHELL approximation
+    if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'CURV', prws(:, :, :) )
+#ifndef MNH_OPENACC
+ ZWORK3(:,:,:) = 1.0 / ( XRADIUS + MZF(PZZ(:,:,:)) )
+ ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=JKU )
+ ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=JKU )
+#else
+CALL MZF_DEVICE(PZZ,ZWORK3) !MZF(PZZ(:,:,:))
+!$acc kernels present_cr(ZWORK3,ZWORK1,ZWORK2)
+    ZWORK3(:,:,:) = 1.0 / ( XRADIUS + ZWORK3(:,:,:) )
+    ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=JKU )
+    ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=JKU )
+!$acc end kernels
+#endif
+    IF (MPPDB_INITIALIZED) THEN
+    CALL MPPDB_CHECK3DM("DYN_SOURCES:ZWORK3,ZWORK1,ZWORK2",PRECISION,&
+                      & ZWORK3,ZWORK1,ZWORK2,&
+                      & MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 ) , &
+                      & MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 ) ,&
+                      & MYF(PVT) * ZWORK1 - MZF(PWT) , &
+                      & MYF(PVT) , MZF(PWT) , MXM(PWT) , MYM(PWT) )
+    CALL MPPDB_CHECK3DM("DYN_SOOURCES:SUITE",PRECISION,&
+         &  MXM(ZRVT),MXM(PVT),MXM(PWT),MXM(ZWORK1),MXM(ZWORK2),MXM(ZWORK3),&
+         &  ZRUT,ZRVT,PRUS,PRVS,PRWS )
+    ENDIF
+!
+#ifndef MNH_OPENACC
+    PRUS(:,:,:) = PRUS                                              &
+    + MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 )                        &
+    + ZRUT * MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 )
+!
+    PRVS(:,:,:) = PRVS                                              &
+    - MYM( MXF(ZRUT*PUT) * ZWORK1 * ZWORK3 )                        &
+    - ZRVT * MYM( (MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3 )
+!
+    PRWS(:,:,:) = PRWS                                              &
+    +MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+#else
+!$acc kernels present_cr(ZRVTT,ZRUTT)
+ZRVTT(:,:,:)=ZRVT(:,:,:)*PVT(:,:,:)
+ZRUTT(:,:,:)=ZRUT(:,:,:)*PUT(:,:,:)
+
+!PRUS(:,:,:) = PRUS                                              &
+!   + MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 )                        &
+!   + ZRUT * MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 )
+!$acc end kernels
+CALL MYF_DEVICE(ZRVTT,ZTMP1_DEVICE )  !MYF(ZRVT*PVT)
+CALL MYF_DEVICE(PVT,ZTMP2_DEVICE )    !MYF(PVT)
+CALL MZF_DEVICE(PWT,ZTMP3_DEVICE )    !MZF(PWT)
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP4_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZTMP1_DEVICE(:,:,:)* ZWORK2(:,:,:) * ZWORK3(:,:,:)  ! MYF(ZRVT*PVT) * ZWORK2 * ZWORK3
+ZTMP4_DEVICE(:,:,:)=(ZTMP2_DEVICE(:,:,:)* ZWORK1(:,:,:) - ZTMP3_DEVICE(:,:,:)) * ZWORK3(:,:,:)  !( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3
+!$acc end kernels
+CALL MXM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) ! MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 ) 
+CALL MXM_DEVICE(ZTMP4_DEVICE,ZTMP5_DEVICE) ! MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 )
+!$acc kernels present_cr(PRUS)
+    PRUS(:,:,:) = PRUS(:,:,:)                                &
+    +ZTMP2_DEVICE(:,:,:)                                     &
+    + ZRUT(:,:,:) * ZTMP5_DEVICE(:,:,:)
+!$acc end kernels
+
+!    PRVS(:,:,:) = PRVS                                              &
+!    - MYM( MXF(ZRUT*PUT) * ZWORK1 * ZWORK3 )                        &
+!    - ZRVT * MYM( (MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3 )
+CALL MXF_DEVICE(ZRUTT,ZTMP1_DEVICE ) !MXF(ZRUT*PUT)
+CALL MXF_DEVICE(PUT,ZTMP2_DEVICE )   !MXF(PUT)
+!CALL MZF_DEVICE(PWT,ZTMP3_DEVICE )   !MZF(PWT)
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP4_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZTMP1_DEVICE(:,:,:)* ZWORK1(:,:,:) * ZWORK3(:,:,:) !MXF(ZRUT*PUT) * ZWORK1 * ZWORK3
+ZTMP4_DEVICE(:,:,:)=(ZTMP2_DEVICE(:,:,:)* ZWORK2(:,:,:) + ZTMP3_DEVICE(:,:,:)) * ZWORK3(:,:,:) !(MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3
+!$acc end kernels
+CALL MYM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !MYM( MXF(ZRUT*PUT) * ZWORK1 * ZWORK3 )
+CALL MYM_DEVICE(ZTMP4_DEVICE,ZTMP5_DEVICE) !MYM( (MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3 )
+!$acc kernels present_cr(PRVS)
+    PRVS(:,:,:) = PRVS(:,:,:)                                &
+    - ZTMP2_DEVICE(:,:,:)                                     &
+    - ZRVT(:,:,:) * ZTMP5_DEVICE(:,:,:)
+!$acc end kernels
+
+!    PRWS(:,:,:) = PRWS                                              &
+!    +MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+CALL MYF_DEVICE(ZRVTT,ZTMP1_DEVICE ) !MXF(ZRUT*PUT)
+CALL MXF_DEVICE(ZRUTT,ZTMP2_DEVICE ) !MYF(ZRVT*PVT)
+!$acc kernels present_cr(ZTMP3_DEVICE)
+ZTMP3_DEVICE(:,:,:)=(ZTMP2_DEVICE(:,:,:)+ZTMP1_DEVICE(:,:,:)) * ZWORK3(:,:,:) ! (MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3
+!$acc end kernels
+CALL MZM_DEVICE(ZTMP3_DEVICE,ZTMP1_DEVICE) ! MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+!$acc kernels present_cr(PRWS)
+PRWS(:,:,:) = PRWS(:,:,:)+ZTMP1_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+
+
+    if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'CURV', prws(:, :, :) )
+  END IF
+
+  if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'CURV', prus(:, :, :) )
+  if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'CURV', prvs(:, :, :) )
+END IF
+!
+!
+!-------------------------------------------------------------------------------
+!
+!*       3.     COMPUTES THE CORIOLIS TERMS
+!	        ---------------------------
+!
+IF (LCORIO)   THEN 
+  if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'COR', prus(:, :, :) )
+  if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'COR', prvs(:, :, :) )
+!$acc kernels present_cr(ZWORK3)
+  ZWORK3(:,:,:) = SPREAD( PCORIOZ(:,:),DIM=3,NCOPIES=JKU ) * PRHODJ(:,:,:)
+!$acc end kernels
+#ifndef MNH_OPENACC
+PRUS(:,:,:) = PRUS + MXM( ZWORK3 * MYF(PVT) )
+PRVS(:,:,:) = PRVS - MYM( ZWORK3 * MXF(PUT) )
+#else
+CALL MYF_DEVICE(PVT,ZTMP1_DEVICE) !MYF(PVT)
+CALL MXF_DEVICE(PUT,ZTMP2_DEVICE) !MXF(PUT)
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP2_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZWORK3(:,:,:) * ZTMP1_DEVICE(:,:,:)
+ZTMP2_DEVICE(:,:,:)=ZWORK3(:,:,:) * ZTMP2_DEVICE(:,:,:)
+!$acc end kernels
+CALL MXM_DEVICE(ZTMP1_DEVICE, ZTMP3_DEVICE) ! MXM( ZWORK3 * MYF(PVT)
+CALL MYM_DEVICE(ZTMP2_DEVICE, ZTMP4_DEVICE) ! MYM( ZWORK3 * MXF(PUT)
+!$acc kernels present_cr(PRUS,PRVS)
+PRUS(:,:,:) = PRUS(:,:,:) + ZTMP3_DEVICE(:,:,:)
+PRVS(:,:,:) = PRVS(:,:,:) - ZTMP4_DEVICE(:,:,:)
+!$acc end kernels
+#endif  
+!
+
+IF ((.NOT.LTHINSHELL) .AND. (.NOT.L1D)) THEN
+if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'COR', prws(:, :, :) )
+!$acc kernels present_cr(ZWORK1,ZWORK2)
+ZWORK1(:,:,:) = SPREAD( PCORIOX(:,:),DIM=3,NCOPIES=JKU) * PRHODJ(:,:,:) 
+ZWORK2(:,:,:) = SPREAD( PCORIOY(:,:),DIM=3,NCOPIES=JKU) * PRHODJ(:,:,:)
+!$acc end kernels
+
+CALL MPPDB_CHECK3DM("DYN_SOOURCES:CORIOLIS",PRECISION,&
+ & ZWORK1,ZWORK2,ZWORK3 )
+!
+#ifndef MNH_OPENACC
+    PRUS(:,:,:) = PRUS - MXM( ZWORK2 * MZF(PWT) )
+!
+    PRVS(:,:,:) = PRVS - MYM( ZWORK1 * MZF(PWT) )
+!
+    PRWS(:,:,:) = PRWS + MZM( ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT) )
+#else
+!***PRUS(:,:,:) = PRUS - MXM( ZWORK2 * MZF(PWT) )
+CALL MZF_DEVICE(PWT,ZTMP1_DEVICE) !MZF(PWT)
+!$acc kernels present_cr(ZTMP2_DEVICE)
+ZTMP2_DEVICE(:,:,:)=ZWORK2(:,:,:) * ZTMP1_DEVICE(:,:,:)    !ZWORK2 * MZF(PWT)
+!$acc end kernels
+CALL MXM_DEVICE(ZTMP2_DEVICE,ZTMP3_DEVICE) !MXM( ZWORK2 * MZF(PWT) )
+!$acc kernels present_cr(PRUS)
+    PRUS(:,:,:) = PRUS(:,:,:) - ZTMP3_DEVICE(:,:,:)
+!$acc end kernels
+
+!*** PRVS(:,:,:) = PRVS - MYM( ZWORK1 * MZF(PWT) )
+!$acc kernels present_cr(ZTMP2_DEVICE)
+ZTMP2_DEVICE(:,:,:)=ZWORK1(:,:,:) * ZTMP1_DEVICE(:,:,:)    !ZWORK1 * MZF(PWT)
+!$acc end kernels
+CALL MYM_DEVICE(ZTMP2_DEVICE,ZTMP3_DEVICE) !MYM( ZWORK1 * MZF(PWT) )
+!$acc kernels present_cr(PRVS)
+    PRVS(:,:,:) = PRVS(:,:,:) - ZTMP2_DEVICE(:,:,:)
+!$acc end kernels
+
+!*** PRWS(:,:,:) = PRWS + MZM( ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT) )
+CALL MXF_DEVICE(PUT,ZTMP1_DEVICE) !MXF(PUT)
+CALL MYF_DEVICE(PVT,ZTMP2_DEVICE) !MYF(PVT)
+!$acc kernels present_cr(ZTMP3_DEVICE)
+ZTMP3_DEVICE(:,:,:)= ZWORK2(:,:,:) * ZTMP1_DEVICE(:,:,:)+ ZWORK1(:,:,:) * ZTMP2_DEVICE(:,:,:) !ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT)
+!$acc end kernels
+CALL MZM_DEVICE(ZTMP3_DEVICE,ZTMP1_DEVICE) !MZM( ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT) )
+!$acc kernels present_cr(PRWS)
+    PRWS(:,:,:) = PRWS(:,:,:) + ZTMP1_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+
+    if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'COR', prws(:, :, :) )
+  END IF
+
+  if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'COR', prus(:, :, :) )
+  if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'COR', prvs(:, :, :) )
+END IF
+!!!!$acc end data
+!
+!-------------------------------------------------------------------------------
+!
+!*       4.     COMPUTES THE THETA SOURCE TERM DUE TO THE REFERENCE PRESSURE
+!	        ------------------------------------------------------------
+!
+IF (LCARTESIAN .OR. LTHINSHELL) THEN
+  ZD1=0.
+ELSE
+  ZD1=1.
+ENDIF
+!
+IF( .NOT.L1D ) THEN
+ IF (.NOT. LOCEAN) THEN
+!
+  IF (KRR > 0) THEN
+    if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'PREF', prths(:, :, :) )
+
+
+    ZCPD_OV_RD = XCPD / XRD
+    ZG_OV_CPD  = -XG / XCPD
+!
+! stores the specific heat capacity (Cph) in ZWORK1
+
+!$acc kernels present_cr( ZWORK1)
+    ZWORK1(:,:,:) = XCPD + XCPV * PRT(:,:,:,1)  ! gas mixing
+!$acc loop seq
+    DO JWATER = 2,1+KRRL             !  loop on the liquid components  
+      ZWORK1(:,:,:) = ZWORK1(:,:,:) + XCL * PRT(:,:,:,JWATER)
+    END DO
+!$acc loop seq
+    DO JWATER = 2+KRRL,1+KRRL+KRRI   !  loop on the solid components   
+      ZWORK1(:,:,:) = ZWORK1(:,:,:) + XCI * PRT(:,:,:,JWATER)
+    END DO
+!$acc end kernels
+! computes the source term
+!
+#ifndef MNH_OPENACC
+    PRTHS(:,:,:) = PRTHS(:,:,:) +  PRHODJ(:,:,:)                             &
+      * ( ( XRD + XRV * PRT(:,:,:,1) ) * ZCPD_OV_RD / ZWORK1(:,:,:)  - 1. )  &
+      * PTHT(:,:,:)/PEXNREF(:,:,:)*MZF(PWT(:,:,:))*(ZG_OV_CPD/PTHVREF(:,:,:) &
+      -ZD1*4./7.*PEXNREF(:,:,:)/( XRADIUS+MZF(PZZ(:,:,:)) ))
+#else
+CALL MZF_DEVICE(PWT,ZTMP1_DEVICE) !MZF(PWT(:,:,:))
+CALL MZF_DEVICE(PZZ(:,:,:),ZTMP2_DEVICE(:,:,:))  !MZF(PZZ(:,:,:))
+!$acc kernels present_cr(PRTHS)
+PRTHS(:,:,:) = PRTHS(:,:,:) +  PRHODJ(:,:,:)                             &
+      * ( ( XRD + XRV * PRT(:,:,:,1) ) * ZCPD_OV_RD / ZWORK1(:,:,:)  - 1. )  &
+      * PTHT(:,:,:)/PEXNREF(:,:,:)*ZTMP1_DEVICE*(ZG_OV_CPD/PTHVREF(:,:,:) &
+      -ZD1*4./7.*PEXNREF(:,:,:)/( XRADIUS+ZTMP2_DEVICE(:,:,:) ))
+!$acc end kernels
+#endif
+    if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'PREF', prths(:, :, :) )
+  END IF
+ END IF
+!
+END IF
+!
+!-------------------------------------------------------------------------------
+#ifdef MNH_OPENACC
+CALL MNH_MEM_RELEASE()
+#endif
+!$acc end data
+!$acc update self(PRUS, PRVS, PRWS,PRTHS)
+END SUBROUTINE DYN_SOURCES 
diff --git a/src/ZSOLVER/get_halo.f90 b/src/ZSOLVER/get_halo.f90
index e53338d1f9e3db6a9798375bf8c17f4289b2ab20..a8665d2059ebde26c45b582cbbce6d343be75aba 100644
--- a/src/ZSOLVER/get_halo.f90
+++ b/src/ZSOLVER/get_halo.f90
@@ -88,23 +88,49 @@ INTERFACE
 END INTERFACE
 !
 INTERFACE
-   SUBROUTINE GET_HALO_START_D(PSRC,KNB_REQ,KREQ,HDIR)
+   SUBROUTINE GET_HALO_START_D(PSRC,KNB_REQ,KREQ,&
+     PZNORTH_IN , PZSOUTH_IN , PZWEST_IN , PZEAST_IN , &
+     PZNORTH_OUT, PZSOUTH_OUT, PZWEST_OUT, PZEAST_OUT, &
+     KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,KIHALO_1,&
+     HDIR)
+     
      IMPLICIT NONE
      !
-     REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+     REAL, DIMENSION(KIIU,KIJU,KIKU), INTENT(INOUT) :: PSRC    ! variable at t
      ! acc declare present (PSRC)
      INTEGER                               :: KNB_REQ ,  KREQ(8)
+     REAL            :: PZSOUTH_IN ( KIIB:KIIE   , KIJB:KIJB+KIHALO_1   , KIKU ) ,&
+                        PZNORTH_IN ( KIIB:KIIE   , KIJE-KIHALO_1:KIJE   , KIKU ) ,&
+                        PZWEST_IN  ( KIIB:KIIB+KIHALO_1   , KIJB:KIJE   , KIKU ) ,&
+                        PZEAST_IN  ( KIIE-KIHALO_1:KIIE   , KIJB:KIJE   , KIKU ) ,&
+       !
+                        PZSOUTH_OUT (   KIIB:KIIE   , 1:KIJB-1 , KIKU ) ,&
+                        PZNORTH_OUT (   KIIB:KIIE   , KIJE+1:KIJU , KIKU ) ,&
+                        PZWEST_OUT  ( 1:KIIB-1 ,   KIJB:KIJE   , KIKU ) ,&
+                        PZEAST_OUT  ( KIIE+1:KIIU ,   KIJB:KIJE   , KIKU ) 
+     INTEGER         :: KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,KIHALO_1
      CHARACTER(len=4), OPTIONAL :: HDIR ! to send only halo on X or Y direction
      !
    END SUBROUTINE GET_HALO_START_D
 END INTERFACE
 INTERFACE
-   SUBROUTINE GET_HALO_STOP_D(PSRC,KNB_REQ,KREQ,HDIR)
+   SUBROUTINE GET_HALO_STOP_D(PSRC,KNB_REQ,KREQ,&
+     PZNORTH_OUT, PZSOUTH_OUT, PZWEST_OUT, PZEAST_OUT, &
+     KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,&
+     HDIR)
+     
      IMPLICIT NONE
-     REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+     !
+     REAL, DIMENSION(KIIU,KIJU,KIKU), INTENT(INOUT) :: PSRC    ! variable at t
      ! acc declare present (PSRC)
-     INTEGER                               :: KNB_REQ , KREQ(8)
+     INTEGER                               :: KNB_REQ ,  KREQ(8)
+     REAL            :: PZSOUTH_OUT (   KIIB:KIIE   , 1:KIJB-1 , KIKU ) ,&
+                        PZNORTH_OUT (   KIIB:KIIE   , KIJE+1:KIJU , KIKU ) ,&
+                        PZWEST_OUT  ( 1:KIIB-1 ,   KIJB:KIJE   , KIKU ) ,&
+                        PZEAST_OUT  ( KIIE+1:KIIU ,   KIJB:KIJE   , KIKU ) 
+     INTEGER         :: KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU
      CHARACTER(len=4), OPTIONAL :: HDIR ! to send only halo on X or Y direction
+     !
    END SUBROUTINE GET_HALO_STOP_D
 END INTERFACE
 INTERFACE
@@ -409,8 +435,9 @@ END MODULE MODD_HALO_D
 USE MODD_HALO_D
 
 !USE MODE_MNH_ZWORK, ONLY : GWEST , GEAST, GSOUTH , GNORTH
-!USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
-!USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE 
+USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
+USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE
+
 !!
 !USE MODE_DEVICE
 USE MODE_MPPDB
@@ -427,22 +454,30 @@ INTEGER            :: INB_REQ , IREQ(8)
 !
 CALL MPPDB_CHECK(PSRC,"GET_HALO_D big:PSRC")
 !
-CALL GET_HALO_START_D(PSRC,INB_REQ,IREQ,HDIR)
-CALL GET_HALO_STOP_D(PSRC,INB_REQ,IREQ,HDIR)
+CALL GET_HALO_START_D(PSRC,INB_REQ,IREQ,&
+     ZNORTH_IN , ZSOUTH_IN , ZWEST_IN , ZEAST_IN , &
+     ZNORTH_OUT, ZSOUTH_OUT, ZWEST_OUT, ZEAST_OUT, &
+     IIB,IIE,IJB,IJE,IIU,IJU,IKU,IHALO_1,HDIR)
+
+CALL GET_HALO_STOP_D(PSRC,INB_REQ,IREQ,&
+     ZNORTH_OUT, ZSOUTH_OUT, ZWEST_OUT, ZEAST_OUT, &
+     IIB,IIE,IJB,IJE,IIU,IJU,IKU,HDIR)
 !
 CALL MPPDB_CHECK(PSRC,"GET_HALO_D end:PSRC")
 !
 END SUBROUTINE GET_HALO_D
 !     #########################
-      SUBROUTINE GET_HALO_START_D(PSRC,KNB_REQ,KREQ,HDIR)
+SUBROUTINE GET_HALO_START_D(PSRC,KNB_REQ,KREQ,&
+     PZNORTH_IN , PZSOUTH_IN , PZWEST_IN , PZEAST_IN , &
+     PZNORTH_OUT, PZSOUTH_OUT, PZWEST_OUT, PZEAST_OUT, &
+     KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,KIHALO_1,&
+     HDIR)
 !     #########################
 #define MNH_GPUDIRECT
 !
 USE MODD_HALO_D
 
 USE MODE_MNH_ZWORK, ONLY : GWEST , GEAST, GSOUTH , GNORTH
-USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
-USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE 
 !
 USE MODD_VAR_ll, ONLY : NMNH_COMM_WORLD
 USE MODD_MPIF,    ONLY : MPI_STATUSES_IGNORE
@@ -453,9 +488,20 @@ USE MODE_MPPDB
 !
 IMPLICIT NONE
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+REAL, DIMENSION(KIIU,KIJU,KIKU), INTENT(INOUT) :: PSRC    ! variable at t
 ! acc declare present (PSRC)
 INTEGER                               :: KNB_REQ ,  KREQ(8)
+REAL            :: PZSOUTH_IN ( KIIB:KIIE   , KIJB:KIJB+KIHALO_1   , KIKU ) ,&
+                   PZNORTH_IN ( KIIB:KIIE   , KIJE-KIHALO_1:KIJE   , KIKU ) ,&
+                   PZWEST_IN  ( KIIB:KIIB+KIHALO_1   , KIJB:KIJE   , KIKU ) ,&
+                   PZEAST_IN  ( KIIE-KIHALO_1:KIIE   , KIJB:KIJE   , KIKU ) ,&
+       !
+                   PZSOUTH_OUT (   KIIB:KIIE   , 1:KIJB-1 , KIKU ) ,&
+                   PZNORTH_OUT (   KIIB:KIIE   , KIJE+1:KIJU , KIKU ) ,&
+                   PZWEST_OUT  ( 1:KIIB-1 ,   KIJB:KIJE   , KIKU ) ,&
+                   PZEAST_OUT  ( KIIE+1:KIIU ,   KIJB:KIJE   , KIKU ) 
+INTEGER         :: KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,KIHALO_1
+
 CHARACTER(len=4), OPTIONAL :: HDIR ! to send only halo on X or Y direction
 !
 INTEGER                          :: IERROR                 ! error return code 
@@ -464,19 +510,15 @@ INTEGER,PARAMETER :: IS_WEST=1 , IS_EAST=2, IS_SOUTH=3, IS_NORTH=4
 LOGICAL      :: LX , LY
 INTEGER      :: NB_REQ, IERR
 !
-INTEGER :: JI,JJ,JK, JIU,JJU,JKU
-
-JIU = SIZE(PSRC,1)
-JJU = SIZE(PSRC,2)
-JKU = SIZE(PSRC,3)
+INTEGER :: JI,JJ,JK 
 
 CALL INIT_HALO_D()
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 !$acc data present (PSRC) &
-!$acc present (ZNORTH_IN, ZSOUTH_IN, ZWEST_IN, ZEAST_IN) &
-!$acc present (ZNORTH_OUT, ZSOUTH_OUT, ZWEST_OUT, ZEAST_OUT) 
+!$acc present (PZNORTH_IN, PZSOUTH_IN, PZWEST_IN, PZEAST_IN) &
+!$acc present (PZNORTH_OUT, PZSOUTH_OUT, PZWEST_OUT, PZEAST_OUT) 
 
 
 LX = .FALSE.
@@ -490,7 +532,7 @@ ELSE
 !!$LY = ( HDIR == "01_Y"  .OR. HDIR == "S0_Y" )
 LX = ( HDIR == "01_X" .OR. HDIR == "S0_X" .OR. HDIR == "S0_Y" )
 LY = ( HDIR == "01_Y" .OR. HDIR == "S0_Y" .OR. HDIR == "S0_X" )
-!!$print *,"IIB=",IIB," HDIR=",HDIR," LX=",LX," LY=",LY ; call flush(6)
+!!$print *,"KIIB=",KIIB," HDIR=",HDIR," LX=",LX," LY=",LY ; call flush(6)
 END IF
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -504,10 +546,10 @@ NB_REQ = 0
 IF (LX) THEN
    IF (.NOT. GWEST) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZWEST_OUT)
+      !$acc host_data use_device(PZWEST_OUT)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_IRECV(ZWEST_OUT,SIZE(ZWEST_OUT),MNHREAL_MPI,NP_WEST-1,1000+IS_EAST,&
+      CALL MPI_IRECV(PZWEST_OUT,SIZE(PZWEST_OUT),MNHREAL_MPI,NP_WEST-1,1000+IS_EAST,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -515,10 +557,10 @@ IF (LX) THEN
    END IF
    IF (.NOT.GEAST) THEN 
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZEAST_OUT)
+      !$acc host_data use_device(PZEAST_OUT)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_IRECV(ZEAST_OUT,SIZE(ZEAST_OUT),MNHREAL_MPI,NP_EAST-1,1000+IS_WEST,&
+      CALL MPI_IRECV(PZEAST_OUT,SIZE(PZEAST_OUT),MNHREAL_MPI,NP_EAST-1,1000+IS_WEST,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -529,10 +571,10 @@ END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZSOUTH_OUT)
+      !$acc host_data use_device(PZSOUTH_OUT)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_IRECV(ZSOUTH_OUT,SIZE(ZSOUTH_OUT),MNHREAL_MPI,NP_SOUTH-1,1000+IS_NORTH,&
+      CALL MPI_IRECV(PZSOUTH_OUT,SIZE(PZSOUTH_OUT),MNHREAL_MPI,NP_SOUTH-1,1000+IS_NORTH,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -540,10 +582,10 @@ IF (LY) THEN
    ENDIF
    IF (.NOT.GNORTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZNORTH_OUT)
+      !$acc host_data use_device(PZNORTH_OUT)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_IRECV(ZNORTH_OUT,SIZE(ZNORTH_OUT),MNHREAL_MPI,NP_NORTH-1,1000+IS_SOUTH,&
+      CALL MPI_IRECV(PZNORTH_OUT,SIZE(PZNORTH_OUT),MNHREAL_MPI,NP_NORTH-1,1000+IS_SOUTH,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -558,15 +600,15 @@ END IF
 IF (LX) THEN
    IF (.NOT. GWEST) THEN
       !$acc kernels async(IS_WEST)
-      !$mnh_expand_array(JI=IIB:IIB+IHALO_1 , JJ=IJB:IJE , JK=1:JKU )
-           ZWEST_IN ( IIB:IIB+IHALO_1  ,    IJB:IJE  , : )  = PSRC( IIB:IIB+IHALO_1  ,  IJB:IJE  , : )
+      !$mnh_expand_array(JI=KIIB:KIIB+KIHALO_1 , JJ=KIJB:KIJE , JK=1:KIKU )
+           PZWEST_IN ( KIIB:KIIB+KIHALO_1  ,    KIJB:KIJE  , : )  = PSRC( KIIB:KIIB+KIHALO_1  ,  KIJB:KIJE  , : )
       !$mnh_end_expand_array()
       !$acc end kernels
    END IF
    IF (.NOT.GEAST) THEN
       !$acc kernels async(IS_EAST)
-      !$mnh_expand_array(JI=IIE-IHALO_1:IIE , JJ=IJB:IJE , JK=1:JKU)
-           ZEAST_IN ( IIE-IHALO_1:IIE  ,    IJB:IJE  , : )  = PSRC( IIE-IHALO_1:IIE  ,  IJB:IJE  , : )
+      !$mnh_expand_array(JI=KIIE-KIHALO_1:KIIE , JJ=KIJB:KIJE , JK=1:KIKU)
+           PZEAST_IN ( KIIE-KIHALO_1:KIIE  ,    KIJB:KIJE  , : )  = PSRC( KIIE-KIHALO_1:KIIE  ,  KIJB:KIJE  , : )
       !$mnh_end_expand_array()
       !$acc end kernels
    ENDIF
@@ -575,15 +617,15 @@ END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
       !$acc kernels async(IS_SOUTH)
-      !$mnh_expand_array(JI=IIB:IIE , JJ=IJB:IJB+IHALO_1 , JK=1:JKU )
-           ZSOUTH_IN ( IIB:IIE  ,    IJB:IJB+IHALO_1  , : ) = PSRC( IIB:IIE  ,    IJB:IJB+IHALO_1  , : )
+      !$mnh_expand_array(JI=KIIB:KIIE , JJ=KIJB:KIJB+KIHALO_1 , JK=1:KIKU )
+           PZSOUTH_IN ( KIIB:KIIE  ,    KIJB:KIJB+KIHALO_1  , : ) = PSRC( KIIB:KIIE  ,    KIJB:KIJB+KIHALO_1  , : )
       !$mnh_end_expand_array()
       !$acc end kernels
    ENDIF
    IF (.NOT.GNORTH) THEN
       !$acc kernels async(IS_NORTH)
-      !$mnh_expand_array(JI=IIB:IIE , JJ=IJE-IHALO_1:IJE , JK=1:JKU )
-           ZNORTH_IN ( IIB:IIE  ,    IJE-IHALO_1:IJE  , : ) = PSRC( IIB:IIE  ,    IJE-IHALO_1:IJE  , : )
+      !$mnh_expand_array(JI=KIIB:KIIE , JJ=KIJE-KIHALO_1:KIJE , JK=1:KIKU )
+           PZNORTH_IN ( KIIB:KIIE  ,    KIJE-KIHALO_1:KIJE  , : ) = PSRC( KIIB:KIIE  ,    KIJE-KIHALO_1:KIJE  , : )
       !$mnh_end_expand_array()
       !$acc end kernels
    ENDIF   
@@ -597,12 +639,12 @@ ENDIF
 IF (LX) THEN
    IF (.NOT. GWEST) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZWEST_IN)
+      !$acc host_data use_device(PZWEST_IN)
 #else
-      !$acc update host(ZWEST_IN)
+      !$acc update host(PZWEST_IN)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_ISEND(ZWEST_IN,SIZE(ZWEST_IN)  ,MNHREAL_MPI,NP_WEST-1,1000+IS_WEST,&
+      CALL MPI_ISEND(PZWEST_IN,SIZE(PZWEST_IN)  ,MNHREAL_MPI,NP_WEST-1,1000+IS_WEST,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -610,12 +652,12 @@ IF (LX) THEN
    END IF
    IF (.NOT.GEAST) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZEAST_IN)
+      !$acc host_data use_device(PZEAST_IN)
 #else
-      !$acc update host(ZEAST_IN)
+      !$acc update host(PZEAST_IN)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_ISEND(ZEAST_IN,SIZE(ZEAST_IN)  ,MNHREAL_MPI,NP_EAST-1,1000+IS_EAST,&
+      CALL MPI_ISEND(PZEAST_IN,SIZE(PZEAST_IN)  ,MNHREAL_MPI,NP_EAST-1,1000+IS_EAST,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -626,12 +668,12 @@ END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZSOUTH_IN)
+      !$acc host_data use_device(PZSOUTH_IN)
 #else
-      !$acc update host(ZSOUTH_IN)
+      !$acc update host(PZSOUTH_IN)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_ISEND(ZSOUTH_IN,SIZE(ZSOUTH_IN)  ,MNHREAL_MPI,NP_SOUTH-1,1000+IS_SOUTH,&
+      CALL MPI_ISEND(PZSOUTH_IN,SIZE(PZSOUTH_IN)  ,MNHREAL_MPI,NP_SOUTH-1,1000+IS_SOUTH,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -639,12 +681,12 @@ IF (LY) THEN
    ENDIF
    IF (.NOT.GNORTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZNORTH_IN)
+      !$acc host_data use_device(PZNORTH_IN)
 #else
-      !$acc update host(ZNORTH_IN)
+      !$acc update host(PZNORTH_IN)
 #endif
       NB_REQ = NB_REQ + 1
-      CALL MPI_ISEND(ZNORTH_IN,SIZE(ZNORTH_IN)  ,MNHREAL_MPI,NP_NORTH-1,1000+IS_NORTH,&
+      CALL MPI_ISEND(PZNORTH_IN,SIZE(PZNORTH_IN)  ,MNHREAL_MPI,NP_NORTH-1,1000+IS_NORTH,&
                      NMNH_COMM_WORLD,KREQ(NB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -659,15 +701,18 @@ KNB_REQ = NB_REQ
 END SUBROUTINE GET_HALO_START_D
 !
 !     #########################
-      SUBROUTINE GET_HALO_STOP_D(PSRC,KNB_REQ,KREQ,HDIR)
+SUBROUTINE GET_HALO_STOP_D(PSRC,KNB_REQ,KREQ,&
+     PZNORTH_OUT, PZSOUTH_OUT, PZWEST_OUT, PZEAST_OUT, &
+     KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,&
+     HDIR)
 !     #########################
 #define MNH_GPUDIRECT
 !
 USE MODD_HALO_D
 
 USE MODE_MNH_ZWORK, ONLY : GWEST , GEAST, GSOUTH , GNORTH
-USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
-USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE 
+!!$USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
+!!$USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE 
 !
 USE MODD_VAR_ll, ONLY : NMNH_COMM_WORLD
 USE MODD_MPIF,    ONLY : MPI_STATUSES_IGNORE
@@ -678,9 +723,15 @@ USE MODE_MPPDB
 !
 IMPLICIT NONE
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+REAL, DIMENSION(KIIU,KIJU,KIKU), INTENT(INOUT) :: PSRC    ! variable at t
 ! acc declare present (PSRC)
 INTEGER                               :: KNB_REQ , KREQ(8)
+REAL            :: PZSOUTH_OUT (   KIIB:KIIE   , 1:KIJB-1 , KIKU ) ,&
+                   PZNORTH_OUT (   KIIB:KIIE   , KIJE+1:KIJU , KIKU ) ,&
+                   PZWEST_OUT  ( 1:KIIB-1 ,   KIJB:KIJE   , KIKU ) ,&
+                   PZEAST_OUT  ( KIIE+1:KIIU ,   KIJB:KIJE   , KIKU ) 
+INTEGER         :: KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU
+
 CHARACTER(len=4), OPTIONAL :: HDIR ! to send only halo on X or Y direction
 !
 INTEGER                          :: IERROR                 ! error return code 
@@ -689,19 +740,14 @@ INTEGER,PARAMETER :: IS_WEST=1 , IS_EAST=2, IS_SOUTH=3, IS_NORTH=4
 LOGICAL      :: LX , LY
 INTEGER      :: NB_REQ, IERR
 !
-INTEGER :: JI,JJ,JK, JIU,JJU,JKU
-
-JIU = SIZE(PSRC,1)
-JJU = SIZE(PSRC,2)
-JKU = SIZE(PSRC,3)
+INTEGER :: JI,JJ,JK
 
 CALL INIT_HALO_D()
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 !$acc data  present (PSRC) &
-!$acc present (ZNORTH_IN, ZSOUTH_IN, ZWEST_IN, ZEAST_IN) &
-!$acc present (ZNORTH_OUT, ZSOUTH_OUT, ZWEST_OUT, ZEAST_OUT)
+!$acc & present (PZNORTH_OUT,PZSOUTH_OUT,PZWEST_OUT,PZEAST_OUT)
 
 LX = .FALSE.
 LY = .FALSE. 
@@ -714,7 +760,7 @@ ELSE
 !!$LY = ( HDIR == "01_Y"  .OR. HDIR == "S0_Y" )
 LX = ( HDIR == "01_X" .OR. HDIR == "S0_X" .OR. HDIR == "S0_Y" )
 LY = ( HDIR == "01_Y" .OR. HDIR == "S0_Y" .OR. HDIR == "S0_X" )
-!!$print *,"IIB=",IIB," HDIR=",HDIR," LX=",LX," LY=",LY ; call flush(6)
+!!$print *,"KIIB=",KIIB," HDIR=",HDIR," LX=",LX," LY=",LY ; call flush(6)
 END IF
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -730,21 +776,21 @@ CALL MPI_WAITALL(NB_REQ,KREQ,MPI_STATUSES_IGNORE,IERR)
 IF (LX) THEN
    IF (.NOT.GWEST) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZWEST_OUT) async(IS_WEST)
+   !$acc update device(PZWEST_OUT) async(IS_WEST)
 #endif
    !$acc kernels async(IS_WEST)
-   !$mnh_expand_array(JI=1:IIB-1 , JJ=IJB:IJE , JK=1:JKU )
-        PSRC( 1:IIB-1  ,      IJB:IJE      , : ) = ZWEST_OUT( 1:IIB-1  ,   IJB:IJE    , : )
+   !$mnh_expand_array(JI=1:KIIB-1 , JJ=KIJB:KIJE , JK=1:KIKU )
+        PSRC( 1:KIIB-1  ,      KIJB:KIJE      , : ) = PZWEST_OUT( 1:KIIB-1  ,   KIJB:KIJE    , : )
    !$mnh_end_expand_array()
    !$acc end kernels
    ENDIF
    IF (.NOT.GEAST) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZEAST_OUT) async(IS_EAST)
+   !$acc update device(PZEAST_OUT) async(IS_EAST)
 #endif
    !$acc kernels async(IS_EAST)
-   !$mnh_expand_array(JI=IIE+1:IIU , JJ=IJB:IJE , JK=1:JKU )
-        PSRC( IIE+1:IIU  ,      IJB:IJE      , : ) = ZEAST_OUT( IIE+1:IIU  ,   IJB:IJE    , : )  
+   !$mnh_expand_array(JI=KIIE+1:KIIU , JJ=KIJB:KIJE , JK=1:KIKU )
+        PSRC( KIIE+1:KIIU  ,      KIJB:KIJE      , : ) = PZEAST_OUT( KIIE+1:KIIU  ,   KIJB:KIJE    , : )  
    !$mnh_end_expand_array()
    !$acc end kernels
    ENDIF
@@ -752,21 +798,21 @@ END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZSOUTH_OUT) async(IS_SOUTH)
+   !$acc update device(PZSOUTH_OUT) async(IS_SOUTH)
 #endif
    !$acc kernels async(IS_SOUTH)
-   !$mnh_expand_array(JI=IIB:IIE , JJ=1:IJB-1 , JK=1:JKU )
-        PSRC(      IIB:IIE       ,  1:IJB-1 , : ) = ZSOUTH_OUT(  IIB:IIE     , 1:IJB-1  , : )
+   !$mnh_expand_array(JI=KIIB:KIIE , JJ=1:KIJB-1 , JK=1:KIKU )
+        PSRC(      KIIB:KIIE       ,  1:KIJB-1 , : ) = PZSOUTH_OUT(  KIIB:KIIE     , 1:KIJB-1  , : )
    !$mnh_end_expand_array()
    !$acc end kernels
    ENDIF
    IF (.NOT.GNORTH) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZNORTH_OUT) async(IS_NORTH)
+   !$acc update device(PZNORTH_OUT) async(IS_NORTH)
 #endif
    !$acc kernels async(IS_NORTH)
-   !$mnh_expand_array(JI=IIB:IIE , JJ=IJE+1:IJU , JK=1:JKU )
-        PSRC(      IIB:IIE       , IJE+1:IJU , : ) = ZNORTH_OUT (  IIB:IIE     , IJE+1:IJU  , : )
+   !$mnh_expand_array(JI=KIIB:KIIE , JJ=KIJE+1:KIJU , JK=1:KIKU )
+        PSRC(      KIIB:KIIE       , KIJE+1:KIJU , : ) = PZNORTH_OUT (  KIIB:KIIE     , KIJE+1:KIJU  , : )
    !$mnh_end_expand_array()
    !$acc end kernels
    ENDIF
@@ -2256,33 +2302,74 @@ CALL MPPDB_CHECK(PSRC,"GET_HALO2_DD end:PSRC")
 END SUBROUTINE GET_HALO2_DD
 !-------------------------------------------------------------------------------
 !     ########################################
-      SUBROUTINE GET_HALO2_DF(PSRC, TP_PSRC_HALO2F_ll, HNAME)
+SUBROUTINE GET_HALO2_DF(PSRC, TP_PSRC_HALO2F_ll, HNAME)
 !     ########################################
 #define MNH_GPUDIRECT
+  !
+USE MODD_HALO_D  
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
+USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE  
+IMPLICIT NONE
 !
-USE MODD_HALO_D
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
+TYPE(HALO2LIST_ll), POINTER  :: TP_PSRC_HALO2F_ll  ! halo2 for SRC
+character(len=*), optional, intent(in) :: HNAME ! Name of the field to be added
+!
+!  local var
+!
+REAL , DIMENSION(:,:) , POINTER , CONTIGUOUS :: ZH2F_EAST,ZH2F_WEST,ZH2F_NORTH,ZH2F_SOUTH
+!
+ZH2F_EAST => TP_PSRC_HALO2F_ll%HALO2%EAST
+ZH2F_WEST => TP_PSRC_HALO2F_ll%HALO2%WEST
+ZH2F_NORTH => TP_PSRC_HALO2F_ll%HALO2%NORTH
+ZH2F_SOUTH => TP_PSRC_HALO2F_ll%HALO2%SOUTH
+!
+CALL  GET_HALO2_DF_DIM(PSRC,&
+      ZSOUTH2F_IN, ZNORTH2F_IN, ZWEST2F_IN, ZEAST2F_IN,&
+      ZSOUTH2F_OUT, ZNORTH2F_OUT, ZWEST2F_OUT, ZEAST2F_OUT,&
+      ZH2F_EAST,ZH2F_WEST,ZH2F_NORTH,ZH2F_SOUTH,&
+      IIB,IIE,IJB,IJE,IIU,IJU,IKU,IHALO2_1,IHALO2,&
+      HNAME)
+!
+CONTAINS
+!
+SUBROUTINE GET_HALO2_DF_DIM(PSRC,&
+      PZSOUTH2F_IN, PZNORTH2F_IN, PZWEST2F_IN, PZEAST2F_IN,&
+      PZSOUTH2F_OUT, PZNORTH2F_OUT, PZWEST2F_OUT, PZEAST2F_OUT,&
+      PZH2F_EAST,PZH2F_WEST,PZH2F_NORTH,PZH2F_SOUTH,&
+      KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,KIHALO2_1,KIHALO2,&
+      HNAME)  
+!        
 USE MODE_ll
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll
-USE MODD_PARAMETERS, ONLY : JPHEXT
 !
-USE MODD_IO,        ONLY : GSMONOPROC
 USE MODE_MNH_ZWORK, ONLY : GWEST , GEAST, GSOUTH , GNORTH
-USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU
-USE MODE_MNH_ZWORK, ONLY : IIB,IJB ,IIE,IJE
-!
-USE MODD_CONF, ONLY : NHALO
-USE MODE_DEVICE
+
 USE MODE_MPPDB
 !
-USE MODD_VAR_ll,     ONLY : IP,NPROC,NP1,NP2
+USE MODD_VAR_ll,     ONLY : NPROC
 USE MODD_VAR_ll,     ONLY : NMNH_COMM_WORLD
 USE MODD_MPIF,       ONLY : MPI_STATUSES_IGNORE
 USE MODD_PRECISION,  ONLY : MNHREAL_MPI
 !
 IMPLICIT NONE
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
-TYPE(HALO2LIST_ll), POINTER  :: TP_PSRC_HALO2F_ll  ! halo2 for SRC
+REAL, DIMENSION(KIIU,KIJU,KIKU), INTENT(INOUT) :: PSRC    ! variable at t
+REAL :: PZSOUTH2F_IN ( KIIB:KIIE   , KIJB:KIJB+KIHALO2_1   , KIKU ) ,&
+        PZNORTH2F_IN ( KIIB:KIIE   , KIJE-KIHALO2_1:KIJE   , KIKU ) ,&
+        PZWEST2F_IN  ( KIIB:KIIB+KIHALO2_1   , KIJB:KIJE   , KIKU ) ,&
+        PZEAST2F_IN  ( KIIE-KIHALO2_1:KIIE   , KIJB:KIJE   , KIKU ) ,&
+        !
+        PZSOUTH2F_OUT (   KIIB:KIIE   , KIJB-KIHALO2:KIJB-1 , KIKU ) ,&
+        PZNORTH2F_OUT (   KIIB:KIIE   , KIJE+1:KIJE+KIHALO2 , KIKU ) ,&
+        PZWEST2F_OUT  ( KIIB-KIHALO2:KIIB-1 ,   KIJB:KIJE   , KIKU ) ,&
+        PZEAST2F_OUT  ( KIIE+1:KIIE+KIHALO2 ,   KIJB:KIJE   , KIKU ) ,&
+        !
+        PZH2F_EAST ( KIJU , KIKU ) ,&
+        PZH2F_WEST ( KIJU , KIKU ) ,&
+        PZH2F_NORTH( KIIU , KIKU ) ,&
+        PZH2F_SOUTH( KIIU , KIKU )
+INTEGER ::  KIIB,KIIE,KIJB,KIJE,KIIU,KIJU,KIKU,KIHALO2_1,KIHALO2
 character(len=*), optional, intent(in) :: HNAME ! Name of the field to be added
 !
 character(len=:), allocatable    :: yname
@@ -2295,7 +2382,7 @@ LOGICAL      :: LX , LY
 INTEGER      :: INB_REQ , IREQ(8)
 INTEGER      :: IERR
 
-REAL , DIMENSION(:,:) , POINTER , CONTIGUOUS :: ZH2F_EAST,ZH2F_WEST,ZH2F_NORTH,ZH2F_SOUTH
+!!$REAL , DIMENSION(:,:) , POINTER , CONTIGUOUS :: ZH2F_EAST,ZH2F_WEST,ZH2F_NORTH,ZH2F_SOUTH
 
 CALL MPPDB_CHECK(PSRC,"GET_HALO2_DF beg:PSRC")
 
@@ -2307,8 +2394,9 @@ end if
 CALL INIT_HALO_D()
 
 !$acc data present ( PSRC ) &
-!$acc present (ZNORTH2F_IN, ZSOUTH2F_IN, ZWEST2F_IN, ZEAST2F_IN) &
-!$acc present (ZNORTH2F_OUT, ZSOUTH2F_OUT, ZWEST2F_OUT, ZEAST2F_OUT)
+!$acc & present (PZNORTH2F_IN, PZSOUTH2F_IN, PZWEST2F_IN, PZEAST2F_IN) &
+!$acc & present (PZNORTH2F_OUT, PZSOUTH2F_OUT, PZWEST2F_OUT, PZEAST2F_OUT) &
+!$acc & present (PZH2F_EAST,PZH2F_WEST,PZH2F_NORTH,PZH2F_SOUTH)
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -2326,10 +2414,10 @@ INB_REQ = 0
 IF (LX) THEN
    IF (.NOT. GWEST) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZWEST2F_OUT)
+      !$acc host_data use_device(PZWEST2F_OUT)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_IRECV(ZWEST2F_OUT,SIZE(ZWEST2F_OUT),MNHREAL_MPI,NP_WEST-1,1000+IS_EAST,&
+      CALL MPI_IRECV(PZWEST2F_OUT,SIZE(PZWEST2F_OUT),MNHREAL_MPI,NP_WEST-1,1000+IS_EAST,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2337,10 +2425,10 @@ IF (LX) THEN
    END IF
    IF (.NOT.GEAST) THEN 
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZEAST2F_OUT)
+      !$acc host_data use_device(PZEAST2F_OUT)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_IRECV(ZEAST2F_OUT,SIZE(ZEAST2F_OUT),MNHREAL_MPI,NP_EAST-1,1000+IS_WEST,&
+      CALL MPI_IRECV(PZEAST2F_OUT,SIZE(PZEAST2F_OUT),MNHREAL_MPI,NP_EAST-1,1000+IS_WEST,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2351,10 +2439,10 @@ END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZSOUTH2F_OUT)
+      !$acc host_data use_device(PZSOUTH2F_OUT)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_IRECV(ZSOUTH2F_OUT,SIZE(ZSOUTH2F_OUT),MNHREAL_MPI,NP_SOUTH-1,1000+IS_NORTH,&
+      CALL MPI_IRECV(PZSOUTH2F_OUT,SIZE(PZSOUTH2F_OUT),MNHREAL_MPI,NP_SOUTH-1,1000+IS_NORTH,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2362,10 +2450,10 @@ IF (LY) THEN
    ENDIF
    IF (.NOT.GNORTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZNORTH2F_OUT)
+      !$acc host_data use_device(PZNORTH2F_OUT)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_IRECV(ZNORTH2F_OUT,SIZE(ZNORTH2F_OUT),MNHREAL_MPI,NP_NORTH-1,1000+IS_SOUTH,&
+      CALL MPI_IRECV(PZNORTH2F_OUT,SIZE(PZNORTH2F_OUT),MNHREAL_MPI,NP_NORTH-1,1000+IS_SOUTH,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2380,28 +2468,28 @@ END IF
 IF (LX) THEN
    IF (.NOT. GWEST) THEN
    !$acc kernels async(IS_WEST)
-      ZWEST2F_IN ( IIB:IIB+IHALO2_1  ,    IJB:IJE  , : )  = PSRC( IIB:IIB+IHALO2_1  ,  IJB:IJE  , : )
-!!$      ZWEST2F_IN ( : , : )  = PSRC( IIB+1  , : , : )
+      PZWEST2F_IN ( KIIB:KIIB+KIHALO2_1  ,    KIJB:KIJE  , : )  = PSRC( KIIB:KIIB+KIHALO2_1  ,  KIJB:KIJE  , : )
+!!$      PZWEST2F_IN ( : , : )  = PSRC( KIIB+1  , : , : )
    !$acc end kernels
       END IF
    IF (.NOT.GEAST) THEN
    !$acc kernels async(IS_EAST)
-      ZEAST2F_IN ( IIE-IHALO2_1:IIE  ,    IJB:IJE  , : )  = PSRC( IIE-IHALO2_1:IIE  ,  IJB:IJE  , : )
-!!$      ZEAST2F_IN ( : , : )  = PSRC( IIE-1 ,  :  , : )
+      PZEAST2F_IN ( KIIE-KIHALO2_1:KIIE  ,    KIJB:KIJE  , : )  = PSRC( KIIE-KIHALO2_1:KIIE  ,  KIJB:KIJE  , : )
+!!$      PZEAST2F_IN ( : , : )  = PSRC( KIIE-1 ,  :  , : )
    !$acc end kernels
       ENDIF
 END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
    !$acc kernels async(IS_SOUTH)
-   ZSOUTH2F_IN ( IIB:IIE  ,    IJB:IJB+IHALO2_1  , : ) = PSRC( IIB:IIE  ,    IJB:IJB+IHALO2_1  , : )
-!!$      ZSOUTH2F_IN ( : , : ) = PSRC( : , IJB+1 , : )
+   PZSOUTH2F_IN ( KIIB:KIIE  ,    KIJB:KIJB+KIHALO2_1  , : ) = PSRC( KIIB:KIIE  ,    KIJB:KIJB+KIHALO2_1  , : )
+!!$      PZSOUTH2F_IN ( : , : ) = PSRC( : , KIJB+1 , : )
    !$acc end kernels
       ENDIF
    IF (.NOT.GNORTH) THEN
    !$acc kernels async(IS_NORTH)
-      ZNORTH2F_IN ( IIB:IIE  ,    IJE-IHALO2_1:IJE  , : ) = PSRC( IIB:IIE  ,    IJE-IHALO2_1:IJE  , : )
-!!$      ZNORTH2F_IN ( : , : ) = PSRC( : , IJE-1  , : )      
+      PZNORTH2F_IN ( KIIB:KIIE  ,    KIJE-KIHALO2_1:KIJE  , : ) = PSRC( KIIB:KIIE  ,    KIJE-KIHALO2_1:KIJE  , : )
+!!$      PZNORTH2F_IN ( : , : ) = PSRC( : , KIJE-1  , : )      
    !$acc end kernels
    ENDIF
 ENDIF
@@ -2414,12 +2502,12 @@ ENDIF
 IF (LX) THEN
    IF (.NOT. GWEST) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZWEST2F_IN)
+      !$acc host_data use_device(PZWEST2F_IN)
 #else
-      !$acc update host(ZWEST2F_IN)
+      !$acc update host(PZWEST2F_IN)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_ISEND(ZWEST2F_IN,SIZE(ZWEST2F_IN)  ,MNHREAL_MPI,NP_WEST-1,1000+IS_WEST,&
+      CALL MPI_ISEND(PZWEST2F_IN,SIZE(PZWEST2F_IN)  ,MNHREAL_MPI,NP_WEST-1,1000+IS_WEST,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2427,12 +2515,12 @@ IF (LX) THEN
    END IF
    IF (.NOT.GEAST) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZEAST2F_IN)
+      !$acc host_data use_device(PZEAST2F_IN)
 #else
-      !$acc update host(ZEAST2F_IN)
+      !$acc update host(PZEAST2F_IN)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_ISEND(ZEAST2F_IN,SIZE(ZEAST2F_IN)  ,MNHREAL_MPI,NP_EAST-1,1000+IS_EAST,&
+      CALL MPI_ISEND(PZEAST2F_IN,SIZE(PZEAST2F_IN)  ,MNHREAL_MPI,NP_EAST-1,1000+IS_EAST,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2443,12 +2531,12 @@ END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZSOUTH2F_IN)
+      !$acc host_data use_device(PZSOUTH2F_IN)
 #else
-      !$acc update host(ZSOUTH2F_IN)
+      !$acc update host(PZSOUTH2F_IN)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_ISEND(ZSOUTH2F_IN,SIZE(ZSOUTH2F_IN)  ,MNHREAL_MPI,NP_SOUTH-1,1000+IS_SOUTH,&
+      CALL MPI_ISEND(PZSOUTH2F_IN,SIZE(PZSOUTH2F_IN)  ,MNHREAL_MPI,NP_SOUTH-1,1000+IS_SOUTH,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2456,12 +2544,12 @@ IF (LY) THEN
    ENDIF
    IF (.NOT.GNORTH) THEN
 #ifdef MNH_GPUDIRECT
-      !$acc host_data use_device(ZNORTH2F_IN)
+      !$acc host_data use_device(PZNORTH2F_IN)
 #else
-      !$acc update host(ZNORTH2F_IN)
+      !$acc update host(PZNORTH2F_IN)
 #endif
       INB_REQ = INB_REQ + 1
-      CALL MPI_ISEND(ZNORTH2F_IN,SIZE(ZNORTH2F_IN)  ,MNHREAL_MPI,NP_NORTH-1,1000+IS_NORTH,&
+      CALL MPI_ISEND(PZNORTH2F_IN,SIZE(PZNORTH2F_IN)  ,MNHREAL_MPI,NP_NORTH-1,1000+IS_NORTH,&
                      NMNH_COMM_WORLD,IREQ(INB_REQ),IERR)
 #ifdef MNH_GPUDIRECT
       !$acc end host_data
@@ -2485,44 +2573,40 @@ END IF
 IF (LX) THEN
    IF (.NOT.GWEST) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZWEST2F_OUT) async(IS_WEST)
+   !$acc update device(PZWEST2F_OUT) async(IS_WEST)
 #endif
-   ZH2F_WEST => TP_PSRC_HALO2F_ll%HALO2%WEST  
    !$acc kernels async(IS_WEST)
-   PSRC( 1:IIB-1  ,      IJB:IJE      , : ) = ZWEST2F_OUT( 1:IIB-1  ,   IJB:IJE    , : )
-   ZH2F_WEST( IJB:IJE , : ) = ZWEST2F_OUT( IIB-2, IJB:IJE , : )      
+   PSRC( 1:KIIB-1  ,      KIJB:KIJE      , : ) = PZWEST2F_OUT( 1:KIIB-1  ,   KIJB:KIJE    , : )
+   PZH2F_WEST( KIJB:KIJE , : ) = PZWEST2F_OUT( KIIB-2, KIJB:KIJE , : )      
    !$acc end kernels
    ENDIF
    IF (.NOT.GEAST) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZEAST2F_OUT) async(IS_EAST)
+   !$acc update device(PZEAST2F_OUT) async(IS_EAST)
 #endif
-   ZH2F_EAST => TP_PSRC_HALO2F_ll%HALO2%EAST
    !$acc kernels async(IS_EAST)
-   PSRC( IIE+1:IIU  ,      IJB:IJE      , : ) = ZEAST2F_OUT( IIE+1:IIU  ,   IJB:IJE    , : )
-   ZH2F_EAST( IJB:IJE , : ) = ZEAST2F_OUT( IIE+2 , IJB:IJE , : )
+   PSRC( KIIE+1:KIIU  ,      KIJB:KIJE      , : ) = PZEAST2F_OUT( KIIE+1:KIIU  ,   KIJB:KIJE    , : )
+   PZH2F_EAST( KIJB:KIJE , : ) = PZEAST2F_OUT( KIIE+2 , KIJB:KIJE , : )
    !$acc end kernels
    ENDIF
 END IF
 IF (LY) THEN
    IF (.NOT.GSOUTH) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZSOUTH2F_OUT) async(IS_SOUTH)
+   !$acc update device(PZSOUTH2F_OUT) async(IS_SOUTH)
 #endif
-   ZH2F_SOUTH => TP_PSRC_HALO2F_ll%HALO2%SOUTH    
    !$acc kernels async(IS_SOUTH)
-   PSRC(      IIB:IIE       ,  1:IJB-1 , : ) = ZSOUTH2F_OUT(  IIB:IIE     , 1:IJB-1  , : )
-   ZH2F_SOUTH( IIB:IIE , : ) = ZSOUTH2F_OUT( IIB:IIE , IJB-2 , : )  
+   PSRC(      KIIB:KIIE       ,  1:KIJB-1 , : ) = PZSOUTH2F_OUT(  KIIB:KIIE     , 1:KIJB-1  , : )
+   PZH2F_SOUTH( KIIB:KIIE , : ) = PZSOUTH2F_OUT( KIIB:KIIE , KIJB-2 , : )  
    !$acc end kernels
    ENDIF
    IF (.NOT.GNORTH) THEN
 #ifndef MNH_GPUDIRECT
-   !$acc update device(ZNORTH2F_OUT) async(IS_NORTH)
+   !$acc update device(PZNORTH2F_OUT) async(IS_NORTH)
 #endif
-   ZH2F_NORTH => TP_PSRC_HALO2F_ll%HALO2%NORTH
    !$acc kernels async(IS_NORTH)
-   PSRC(      IIB:IIE       , IJE+1:IJU , : ) = ZNORTH2F_OUT (  IIB:IIE     , IJE+1:IJU  , : )
-   ZH2F_NORTH( IIB:IIE , : ) = ZNORTH2F_OUT ( IIB:IIE , IJE+2 , : )      
+   PSRC(      KIIB:KIIE       , KIJE+1:KIJU , : ) = PZNORTH2F_OUT (  KIIB:KIIE     , KIJE+1:KIJU  , : )
+   PZH2F_NORTH( KIIB:KIIE , : ) = PZNORTH2F_OUT ( KIIB:KIIE , KIJE+2 , : )      
    !$acc end kernels
    ENDIF
 END IF
@@ -2532,6 +2616,8 @@ END IF
 
 CALL MPPDB_CHECK(PSRC,"GET_HALO2_DF end:PSRC")
 
+END SUBROUTINE GET_HALO2_DF_DIM
+
 END SUBROUTINE GET_HALO2_DF
 !
 !     ###################################################
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
index e45406f0d75a57fb559a342a5e47f04a9a50300c..afc9d71dcf2276959a97e60486eb2012d92e86ae 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
@@ -2003,7 +2003,7 @@ end subroutine line_Jacobi_mnh
 
     real, dimension(:,:,:) , pointer, contiguous :: zSr_st , zb_st , zSu_in_st , zSu_out_st
     real, dimension(:)     , pointer, contiguous :: za_k, zb_k, zc_k, zd_k 
-    integer :: ii,ij
+!!$    integer :: ii,ij
 
     integer , parameter :: max_lev = 128
     logical , save                                 :: Lfirst_call_tridiag_mnhallT = .true.
@@ -2020,6 +2020,7 @@ end subroutine line_Jacobi_mnh
 
     real(kind=rl) :: Tijs
 
+    integer :: nib,nie,njb,nje,nzb,nze
     
     if (LUseT ) then
 
@@ -2036,6 +2037,10 @@ end subroutine line_Jacobi_mnh
       zc_k  => vert_coeff%c
       zd_k  => vert_coeff%d
 
+      nib = Lbound(zb_st,1) ; nie = Ubound(zb_st,1)
+      njb = Lbound(zb_st,2) ; nje = Ubound(zb_st,2)
+      nzb = Lbound(zb_st,3) ; nze = Ubound(zb_st,3)
+
       if ( Lfirst_call_level_tridiag_mnhallT(level) ) then
          Lfirst_call_level_tridiag_mnhallT(level) = .false.
 
@@ -2081,48 +2086,82 @@ 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 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
-         do iz=2,nz-1
-            zSr_st(ii,ij,iz) = zb_st(ii,ij,iz) - zd_k(iz) * ( &
-                 zSu_in_st(ii+1,ij,iz) + &
-                 zSu_in_st(ii-1,ij,iz) + &
-                 zSu_in_st(ii,ij+1,iz) + &
-                 zSu_in_st(ii,ij-1,iz) )
-         end do
-         iz=nz
-         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
-         do iz=2,nz-1
-            zSu_out_st(ii,ij,iz) = (zSr_st(ii,ij,iz) / (Tijs*zd_k(iz)) & 
-                 - zSu_out_st(ii,ij,iz-1)*zc_k(iz)) / tmp_k(iz)
-         end do
-         iz=nz
-         zSu_out_st(ii,ij,iz) = (zSr_st(ii,ij,iz) / (Tijs*zd_k(iz)) &
-              -  zSu_out_st(ii,ij,iz-1)*zc_k(iz)) / tmp_k(iz)         
-      !     
-      ! STEP 2: back substitution
-         !$acc loop seq
-         do iz=nz-1,1,-1             
-            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; end do
-      !$mnh_end_do()   
-      !$acc end parallel
-      ! acc end kernels
+
+!!$      print*,"apply_tridiag::level,Sr%ix_min/%ix_max/%icompx_min/%icompx_max",&
+!!$           Sr%ix_min,Sr%ix_max,Sr%icompx_min,Sr%icompx_max
+!!$      print*,"apply_tridiag::level,iib,iie,ijb,ije,nz",level,iib,iie,ijb,ije,nz
+!!$      print*,"apply_tridiag::level,nib,nie,njb,nje,nzb,nze",level,&
+!!$                                   nib,nie,njb,nje,nzb,nze
+!!$      print*,"apply_tridiag::level, zSr_st=L/U/S",level,Lbound(zSr_st),Ubound(zSr_st),Shape(zSr_st)
+!!$      print*,"apply_tridiag::lvel, zc_k  =L/U/S",level,Lbound(zc_k),Ubound(zc_k),Shape(zc_k)
+
+      call  apply_tridiag_solve_mnh_allT_dim(&
+           zSr_st,zSu_out_st,zb_st,zSu_in_st,&
+           zc_k,zd_k,tmp_k,c_k &
+      )
        
-    end if
+      end if
+
+    contains
+      
+      subroutine  apply_tridiag_solve_mnh_allT_dim(&
+           pzSr_st,pzSu_out_st,pzb_st,pzSu_in_st,&
+           pzc_k,pzd_k,ptmp_k,pc_k &
+        )
+
+        implicit none
+
+        real :: pzSr_st    (nib:nie,njb:nje,nzb:nze) ,&
+                pzSu_out_st(nib:nie,njb:nje,nzb:nze) ,&
+                pzb_st     (nib:nie,njb:nje,nzb:nze) ,&
+                pzSu_in_st (nib:nie,njb:nje,nzb:nze) ,&
+                pzc_k(nz),pzd_k(nz),ptmp_k(nz),pc_k(nz)
+        !
+        ! local var
+        !
+        integer :: ii,ij,iz
+                   
+        ! acc kernels present(pzSr_st,pzSu_out_st,pzb_st,pzSu_in_st)
+        !$acc parallel present(pzSr_st,pzSu_out_st,pzb_st,pzSu_in_st)
+        !$acc &        present(pzc_k,pzd_k,ptmp_k,pc_k)
+        !$mnh_do_concurrent ( ii=iib:iie , ij=ijb:ije )   
+        iz=1 
+        pzSr_st(ii,ij,iz) = pzb_st(ii,ij,iz)
+        !$acc loop seq
+        do iz=2,nz-1
+           pzSr_st(ii,ij,iz) = pzb_st(ii,ij,iz) - pzd_k(iz) * ( &
+                pzSu_in_st(ii+1,ij,iz) + &
+                pzSu_in_st(ii-1,ij,iz) + &
+                pzSu_in_st(ii,ij+1,iz) + &
+                pzSu_in_st(ii,ij-1,iz) )
+        end do
+        iz=nz
+        pzSr_st(ii,ij,iz) = pzb_st(ii,ij,iz)
+        !
+        ! Thomas algorithm
+        !
+        iz = 1     
+        pzSu_out_st(ii,ij,iz) = pzSr_st(ii,ij,iz) / (ptmp_k(iz)*Tijs*pzd_k(iz))
+        !$acc loop seq
+        do iz=2,nz-1
+           pzSu_out_st(ii,ij,iz) = (pzSr_st(ii,ij,iz) / (Tijs*pzd_k(iz)) & 
+                - pzSu_out_st(ii,ij,iz-1)*pzc_k(iz)) / ptmp_k(iz)
+        end do
+        iz=nz
+        pzSu_out_st(ii,ij,iz) = (pzSr_st(ii,ij,iz) / (Tijs*pzd_k(iz)) &
+             -  pzSu_out_st(ii,ij,iz-1)*pzc_k(iz)) / ptmp_k(iz)         
+        !     
+        ! STEP 2: back substitution
+        !$acc loop seq
+        do iz=nz-1,1,-1             
+           pzSu_out_st(ii,ij,iz) = pzSu_out_st(ii,ij,iz) & 
+                - pc_k(iz) * pzSu_out_st(ii,ij,iz+1)
+        end do
+        !$mnh_end_do()   
+        !$acc end parallel
+        ! acc end kernels
+         
+      end subroutine apply_tridiag_solve_mnh_allT_dim
 
   end subroutine apply_tridiag_solve_mnh_allT
   !==================================================================