From ac24b04b39e1996beb903c309d3d269b53d21230 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Mon, 3 Jun 2019 16:40:23 +0200
Subject: [PATCH] Philippe 03/06/2019: remove PACK/UNPACK intrinsics (to get
 more performance and better OpenACC support) in rain_ice(_red) subroutines

---
 src/MNH/ice4_fast_rg.f90            |  50 ++++--
 src/MNH/ice4_fast_rh.f90            |  74 +++++++--
 src/MNH/ice4_fast_rs.f90            |  67 ++++++--
 src/MNH/ice4_nucleation_wrapper.f90 |   8 +-
 src/MNH/ice4_rsrimcg_old.f90        |  23 ++-
 src/MNH/ice4_tendencies.f90         |  17 +-
 src/MNH/rain_ice.f90                |  54 +++---
 src/MNH/rain_ice_fast_rg.f90        |  50 ++++--
 src/MNH/rain_ice_fast_rh.f90        |  60 +++++--
 src/MNH/rain_ice_fast_rs.f90        |  68 ++++++--
 src/MNH/rain_ice_nucleation.f90     |  12 +-
 src/MNH/rain_ice_red.f90            | 245 ++++++++++++++++++++--------
 src/MNH/rain_ice_warm.f90           |  14 +-
 13 files changed, 545 insertions(+), 197 deletions(-)

diff --git a/src/MNH/ice4_fast_rg.f90 b/src/MNH/ice4_fast_rg.f90
index a4c0d17e1..8a52f858e 100644
--- a/src/MNH/ice4_fast_rg.f90
+++ b/src/MNH/ice4_fast_rg.f90
@@ -93,6 +93,7 @@ SUBROUTINE ICE4_FAST_RG(KSIZE, LDSOFT, LDCOMPUTE, KRR, &
 !!    -------------
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !
 !*      0. DECLARATIONS
@@ -167,6 +168,7 @@ INTEGER, PARAMETER :: IRCDRYG=1, IRIDRYG=2, IRIWETG=3, IRSDRYG=4, IRSWETG=5, IRR
 LOGICAL, DIMENSION(SIZE(PRHODREF)) :: GDRY, LLDRYG, GMASK
 INTEGER :: IGDRY
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZVEC1, ZVEC2, ZVEC3
+INTEGER, DIMENSION(SIZE(PRHODREF)) :: I1
 INTEGER, DIMENSION(SIZE(PRHODREF)) :: IVEC1, IVEC2
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZZW, &
                                    ZRDRYG_INIT, & !Initial dry growth rate of the graupeln
@@ -247,7 +249,17 @@ ELSE
 ENDIF
 
 ! Wet and dry collection of rs on graupel (6.2.1)
-GDRY(:)=PRST(:)>XRTMIN(5) .AND. PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:)
+IGDRY = 0
+DO JJ = 1, SIZE(GDRY)
+  IF (PRST(JJ)>XRTMIN(5) .AND. PRGT(JJ)>XRTMIN(6) .AND. LDCOMPUTE(JJ)) THEN
+    IGDRY = IGDRY + 1
+    I1(IGDRY) = JJ
+    GDRY(JJ) = .TRUE.
+  ELSE
+    GDRY(JJ) = .FALSE.
+  END IF
+END DO
+
 IF(LDSOFT) THEN
   WHERE(.NOT. GDRY(:))
     PRG_TEND(:, IRSDRYG)=0.
@@ -256,13 +268,14 @@ IF(LDSOFT) THEN
 ELSE
   PRG_TEND(:, IRSDRYG)=0.
   PRG_TEND(:, IRSWETG)=0.
-  IGDRY=COUNT(GDRY(:))
   IF(IGDRY>0)THEN
     !
     !*       6.2.3  select the (PLBDAG,PLBDAS) couplet
     !
-    ZVEC1(1:IGDRY)=PACK(PLBDAG(:), MASK=GDRY(:))
-    ZVEC2(1:IGDRY)=PACK(PLBDAS(:), MASK=GDRY(:))
+    DO JJ = 1, IGDRY
+      ZVEC1(JJ) = PLBDAG(I1(JJ))
+      ZVEC2(JJ) = PLBDAS(I1(JJ))
+    END DO
     !
     !*       6.2.4  find the next lower indice for the PLBDAG and for the PLBDAS
     !               in the geometrical set of (Lbda_g,Lbda_s) couplet use to
@@ -289,7 +302,10 @@ ELSE
                     - XKER_SDRYG(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                          *(ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:)=UNPACK(VECTOR=ZVEC3(1:IGDRY), MASK=GDRY(:), FIELD=0.0)
+    ZZW(:) = 0.
+    DO JJ = 1, IGDRY
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     WHERE(GDRY(:))
       PRG_TEND(:, IRSWETG)=XFSDRYG*ZZW(:)                         & ! RSDRYG
@@ -306,21 +322,32 @@ ENDIF
 !
 !*       6.2.6  accretion of raindrops on the graupeln
 !
-GDRY(:)=PRRT(:)>XRTMIN(3) .AND. PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:)
+IGDRY = 0
+DO JJ = 1, SIZE(GDRY)
+  IF (PRRT(JJ)>XRTMIN(3) .AND. PRGT(JJ)>XRTMIN(6) .AND. LDCOMPUTE(JJ)) THEN
+    IGDRY = IGDRY + 1
+    I1(IGDRY) = JJ
+    GDRY(JJ) = .TRUE.
+  ELSE
+    GDRY(JJ) = .FALSE.
+  END IF
+END DO
+
 IF(LDSOFT) THEN
   WHERE(.NOT. GDRY(:))
     PRG_TEND(:, IRRDRYG)=0.
   END WHERE
 ELSE
   PRG_TEND(:, IRRDRYG)=0.
-  IGDRY=COUNT(GDRY(:))
   !
   IF(IGDRY>0) THEN
     !
     !*       6.2.8  select the (PLBDAG,PLBDAR) couplet
     !
-    ZVEC1(1:IGDRY)=PACK(PLBDAG(:), MASK=GDRY(:))
-    ZVEC2(1:IGDRY)=PACK(PLBDAR(:), MASK=GDRY(:))
+    DO JJ = 1, IGDRY
+      ZVEC1(JJ) = PLBDAG(I1(JJ))
+      ZVEC2(JJ) = PLBDAR(I1(JJ))
+    END DO
     !
     !*       6.2.9  find the next lower indice for the PLBDAG and for the PLBDAR
     !               in the geometrical set of (Lbda_g,Lbda_r) couplet use to
@@ -347,7 +374,10 @@ ELSE
                     - XKER_RDRYG(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                          *(ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:)=UNPACK(VECTOR=ZVEC3(1:IGDRY), MASK=GDRY, FIELD=0.)
+    ZZW(:) = 0.
+    DO JJ = 1, IGDRY
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     WHERE(GDRY(:))
       PRG_TEND(:, IRRDRYG) = XFRDRYG*ZZW(:)                    & ! RRDRYG
diff --git a/src/MNH/ice4_fast_rh.f90 b/src/MNH/ice4_fast_rh.f90
index e590f79a6..9678dceb1 100644
--- a/src/MNH/ice4_fast_rh.f90
+++ b/src/MNH/ice4_fast_rh.f90
@@ -83,6 +83,7 @@ SUBROUTINE ICE4_FAST_RH(KSIZE, LDSOFT, LDCOMPUTE, LDWETG, &
 !!    -------------
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !
 !*      0. DECLARATIONS
@@ -153,6 +154,7 @@ INTEGER, PARAMETER :: IRCWETH=1, IRRWETH=2, IRIDRYH=3, IRIWETH=4, IRSDRYH=5, IRS
 LOGICAL, DIMENSION(SIZE(PRHODREF)) :: GHAIL, GWET, GMASK, LLWETH, LLDRYH
 INTEGER :: IHAIL, IGWET
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZVEC1, ZVEC2, ZVEC3
+INTEGER, DIMENSION(SIZE(PRHODREF)) :: I1
 INTEGER, DIMENSION(SIZE(PRHODREF)) :: IVEC1, IVEC2
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZZW, &
                                    ZRDRYH_INIT, ZRWETH_INIT, &
@@ -195,7 +197,17 @@ ENDIF
 !
 !*       7.2.1  accretion of aggregates on the hailstones
 !
-GWET(:) = PRHT(:)>XRTMIN(7) .AND. PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:)
+IGWET = 0
+DO JJ = 1, SIZE(GWET)
+  IF (PRHT(JJ)>XRTMIN(7) .AND. PRST(JJ)>XRTMIN(5) .AND. LDCOMPUTE(JJ)) THEN
+    IGWET = IGWET + 1
+    I1(IGWET) = JJ
+    GWET(JJ) = .TRUE.
+  ELSE
+    GWET(JJ) = .FALSE.
+  END IF
+END DO
+
 IF(LDSOFT) THEN
   WHERE(.NOT. GWET(:))
     PRH_TEND(:, IRSWETH)=0.
@@ -204,13 +216,14 @@ IF(LDSOFT) THEN
 ELSE
   PRH_TEND(:, IRSWETH)=0.
   PRH_TEND(:, IRSDRYH)=0.
-  IGWET=COUNT(GWET(:))
   IF(IGWET>0)THEN
     !
     !*       7.2.3  select the (PLBDAH,PLBDAS) couplet
     !
-    ZVEC1(1:IGWET) = PACK( PLBDAH(:),MASK=GWET(:) )
-    ZVEC2(1:IGWET) = PACK( PLBDAS(:),MASK=GWET(:) )
+    DO JJ = 1, IGWET
+      ZVEC1(JJ) = PLBDAH(I1(JJ))
+      ZVEC2(JJ) = PLBDAS(I1(JJ))
+    END DO
     !
     !*       7.2.4  find the next lower indice for the PLBDAG and for the PLBDAS
     !               in the geometrical set of (Lbda_h,Lbda_s) couplet use to
@@ -237,7 +250,10 @@ ELSE
                   - XKER_SWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                          * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(1:IGWET),MASK=GWET,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGWET
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     WHERE(GWET(:))
       PRH_TEND(:, IRSWETH)=XFSWETH*ZZW(:)                       & ! RSWETH
@@ -253,7 +269,17 @@ ENDIF
 !
 !*       7.2.6  accretion of graupeln on the hailstones
 !
-GWET(:) = PRHT(:)>XRTMIN(7) .AND. PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:)
+IGWET = 0
+DO JJ = 1, SIZE(GWET)
+  IF (PRHT(JJ)>XRTMIN(7) .AND. PRGT(JJ)>XRTMIN(6) .AND. LDCOMPUTE(JJ)) THEN
+    IGWET = IGWET + 1
+    I1(IGWET) = JJ
+    GWET(JJ) = .TRUE.
+  ELSE
+    GWET(JJ) = .FALSE.
+  END IF
+END DO
+
 IF(LDSOFT) THEN
   WHERE(.NOT. GWET(:))
     PRH_TEND(:, IRGWETH)=0.
@@ -262,13 +288,14 @@ IF(LDSOFT) THEN
 ELSE
   PRH_TEND(:, IRGWETH)=0.
   PRH_TEND(:, IRGDRYH)=0.
-  IGWET=COUNT(GWET(:))
   IF(IGWET>0)THEN
     !
     !*       7.2.8  select the (PLBDAH,PLBDAG) couplet
     !
-    ZVEC1(1:IGWET) = PACK( PLBDAH(:),MASK=GWET(:) )
-    ZVEC2(1:IGWET) = PACK( PLBDAG(:),MASK=GWET(:) )
+    DO JJ = 1, IGWET
+      ZVEC1(JJ) = PLBDAH(I1(JJ))
+      ZVEC2(JJ) = PLBDAG(I1(JJ))
+    END DO
     !
     !*       7.2.9  find the next lower indice for the PLBDAH and for the PLBDAG
     !               in the geometrical set of (Lbda_h,Lbda_g) couplet use to
@@ -295,7 +322,10 @@ ELSE
                    - XKER_GWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                         * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(1:IGWET),MASK=GWET,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGWET
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     WHERE(GWET(:))
       PRH_TEND(:, IRGWETH)=XFGWETH*ZZW(:)                       & ! RGWETH
@@ -315,20 +345,31 @@ ENDIF
 !
 !*       7.2.11  accretion of raindrops on the hailstones
 !
-GWET(:) = PRHT(:)>XRTMIN(7) .AND. PRRT(:)>XRTMIN(3) .AND. LDCOMPUTE(:)
+IGWET = 0
+DO JJ = 1, SIZE(GWET)
+  IF (PRHT(JJ)>XRTMIN(7) .AND. PRRT(JJ)>XRTMIN(3) .AND. LDCOMPUTE(JJ)) THEN
+    IGWET = IGWET + 1
+    I1(IGWET) = JJ
+    GWET(JJ) = .TRUE.
+  ELSE
+    GWET(JJ) = .FALSE.
+  END IF
+END DO
+
 IF(LDSOFT) THEN
   WHERE(.NOT. GWET(:))
     PRH_TEND(:, IRRWETH)=0.
   END WHERE
 ELSE
   PRH_TEND(:, IRRWETH)=0.
-  IGWET=COUNT(GWET(:))
   IF(IGWET>0)THEN
     !
     !*       7.2.12  select the (PLBDAH,PLBDAR) couplet
     !
-    ZVEC1(1:IGWET)=PACK(PLBDAH(:), MASK=GWET(:))
-    ZVEC2(1:IGWET)=PACK(PLBDAR(:), MASK=GWET(:))
+    DO JJ = 1, IGWET
+      ZVEC1(JJ) = PLBDAH(I1(JJ))
+      ZVEC2(JJ) = PLBDAR(I1(JJ))
+    END DO
     !
     !*       7.2.13 find the next lower indice for the PLBDAH and for the PLBDAR
     !               in the geometrical set of (Lbda_h,Lbda_r) couplet use to
@@ -355,7 +396,10 @@ ELSE
                     - XKER_RWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                          *(ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:)=UNPACK(VECTOR=ZVEC3(1:IGWET), MASK=GWET, FIELD=0.)
+    ZZW(:) = 0.
+    DO JJ = 1, IGWET
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     WHERE(GWET(:))
       PRH_TEND(:, IRRWETH) = XFRWETH*ZZW(:)                    & ! RRWETH
diff --git a/src/MNH/ice4_fast_rs.f90 b/src/MNH/ice4_fast_rs.f90
index ca78eb9e1..b4cf9bd3b 100644
--- a/src/MNH/ice4_fast_rs.f90
+++ b/src/MNH/ice4_fast_rs.f90
@@ -76,6 +76,7 @@ SUBROUTINE ICE4_FAST_RS(KSIZE, LDSOFT, LDCOMPUTE, &
 !!    -------------
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !
 !*      0. DECLARATIONS
@@ -134,6 +135,7 @@ INTEGER, PARAMETER :: IRCRIMS=1, IRCRIMSS=2, IRSRIMCG=3, IRRACCS=4, IRRACCSS=5,
 LOGICAL, DIMENSION(SIZE(PRHODREF)) :: GRIM, GACC, GMASK
 INTEGER :: IGRIM, IGACC
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZVEC1, ZVEC2, ZVEC3
+INTEGER, DIMENSION(SIZE(PRHODREF)) :: I1
 INTEGER, DIMENSION(SIZE(PRHODREF)) :: IVEC1, IVEC2
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZZW, ZZW2, ZZW6, ZFREEZ_RATE
 INTEGER :: JJ
@@ -169,7 +171,16 @@ END WHERE
 !
 !*       5.1    cloud droplet riming of the aggregates
 !
-GRIM(:) = PRCT(:)>XRTMIN(2) .AND. PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:)
+IGRIM = 0
+DO JJ = 1, SIZE(GRIM)
+  IF (PRCT(JJ)>XRTMIN(2) .AND. PRST(JJ)>XRTMIN(5) .AND. LDCOMPUTE(JJ)) THEN
+    IGRIM = IGRIM + 1
+    I1(IGRIM) = JJ
+    GRIM(JJ) = .TRUE.
+  ELSE
+    GRIM(JJ) = .FALSE.
+  END IF
+END DO
 !
 ! Collection of cloud droplets by snow: this rate is used for riming (T<0) and for conversion/melting (T>0)
 IF(LDSOFT) THEN
@@ -182,13 +193,14 @@ ELSE
   PRS_TEND(:, IRCRIMS)=0.
   PRS_TEND(:, IRCRIMSS)=0.
   PRS_TEND(:, IRSRIMCG)=0.
-  IGRIM = COUNT(GRIM(:))
   !
   IF(IGRIM>0) THEN
     !
     !        5.1.1  select the PLBDAS
     !
-    ZVEC1(1:IGRIM) = PACK( PLBDAS(:),MASK=GRIM(:) )
+    DO JJ = 1, IGRIM
+      ZVEC1(JJ) = PLBDAS(I1(JJ))
+    END DO
     !
     !        5.1.2  find the next lower indice for the PLBDAS in the geometrical
     !               set of Lbda_s used to tabulate some moments of the incomplete
@@ -204,7 +216,10 @@ ELSE
     !
     ZVEC1(1:IGRIM) =   XGAMINC_RIM1( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                      - XGAMINC_RIM1( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
-    ZZW(:) = UNPACK( VECTOR=ZVEC1(1:IGRIM),MASK=GRIM,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGRIM
+      ZZW(I1(JJ)) = ZVEC1(JJ)
+    END DO
     !
     !        5.1.4  riming of the small sized aggregates
     !
@@ -220,11 +235,17 @@ ELSE
     !
     ZVEC1(1:IGRIM) =  XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                     - XGAMINC_RIM2( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
-    ZZW(:) = UNPACK( VECTOR=ZVEC1(1:IGRIM),MASK=GRIM,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGRIM
+      ZZW(I1(JJ)) = ZVEC1(JJ)
+    END DO
 
     ZVEC1(1:IGRIM) =  XGAMINC_RIM4( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                     - XGAMINC_RIM4( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
-    ZZW2(:) = UNPACK( VECTOR=ZVEC1(1:IGRIM),MASK=GRIM,FIELD=0.0)
+    ZZW2(:) = 0.
+    DO JJ = 1, IGRIM
+      ZZW2(I1(JJ)) = ZVEC1(JJ)
+    END DO
     !
     !        5.1.6  riming-conversion of the large sized aggregates into graupeln
     !
@@ -277,7 +298,17 @@ PA_TH(:) = PA_TH(:) + PRCRIMSG(:)*(PLSFACT(:)-PLVFACT(:))
 !
 !*       5.2    rain accretion onto the aggregates
 !
-GACC(:) = PRRT(:)>XRTMIN(3) .AND. PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:)
+IGACC = 0
+DO JJ = 1, SIZE(GACC)
+  IF (PRRT(JJ)>XRTMIN(3) .AND. PRST(JJ)>XRTMIN(5) .AND. LDCOMPUTE(JJ)) THEN
+    IGACC = IGACC + 1
+    I1(IGACC) = JJ
+    GACC(JJ) = .TRUE.
+  ELSE
+    GACC(JJ) = .FALSE.
+  END IF
+END DO
+
 IF(LDSOFT) THEN
   WHERE(.NOT. GACC(:))
     PRS_TEND(:, IRRACCS)=0.
@@ -288,14 +319,15 @@ ELSE
   PRS_TEND(:, IRRACCS)=0.
   PRS_TEND(:, IRRACCSS)=0.
   PRS_TEND(:, IRSACCRG)=0.
-  IGACC = COUNT(GACC(:))
   IF(IGACC>0)THEN
     !
     !
     !        5.2.1  select the (PLBDAS,PLBDAR) couplet
     !
-    ZVEC1(1:IGACC) = PACK( PLBDAS(:),MASK=GACC(:) )
-    ZVEC2(1:IGACC) = PACK( PLBDAR(:),MASK=GACC(:) )
+    DO JJ = 1, IGACC
+      ZVEC1(JJ) = PLBDAS(I1(JJ))
+      ZVEC2(JJ) = PLBDAR(I1(JJ))
+    END DO
     !
     !        5.2.2  find the next lower indice for the PLBDAS and for the PLBDAR
     !               in the geometrical set of (Lbda_s,Lbda_r) couplet use to
@@ -322,7 +354,10 @@ ELSE
                     - XKER_RACCSS(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                           * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(1:IGACC),MASK=GACC,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGACC
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     !        5.2.4  raindrop accretion on the small sized aggregates
     !
@@ -346,7 +381,10 @@ ELSE
                     -  XKER_RACCS(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                            * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(1:IGACC),MASK=GACC(:),FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGACC
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     WHERE(GACC(:))
       PRS_TEND(:, IRRACCS) = ZZW(:)*ZZW6(:)
     END WHERE
@@ -361,7 +399,10 @@ ELSE
                       - XKER_SACCRG(IVEC2(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
                                                             * (ZVEC2(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(1:IGACC),MASK=GACC,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGACC
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
     !
     !        5.2.6  raindrop accretion-conversion of the large sized aggregates
     !               into graupeln
diff --git a/src/MNH/ice4_nucleation_wrapper.f90 b/src/MNH/ice4_nucleation_wrapper.f90
index 5a4c97f80..6bd37bd55 100644
--- a/src/MNH/ice4_nucleation_wrapper.f90
+++ b/src/MNH/ice4_nucleation_wrapper.f90
@@ -40,6 +40,7 @@ SUBROUTINE ICE4_NUCLEATION_WRAPPER(KIT, KJT, KKT, LDMASK, &
 !!    MODIFICATIONS
 !!    -------------
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !
 !
@@ -135,8 +136,11 @@ IF(INEGT>0) THEN
                        ZTHT, ZPRES, ZRHODREF, ZEXN, ZLSFACT, ZZT, &
                        ZRVT, &
                        ZCIT, ZRVHENI_MR, ZB_TH, ZB_RV, ZB_RI)
-  PRVHENI_MR(:,:,:)=UNPACK(ZRVHENI_MR(:), MASK=GNEGT(:,:,:), FIELD=0.0)
-  PCIT(:,:,:)      =UNPACK(ZCIT(:),       MASK=GNEGT(:,:,:), FIELD=PCIT(:,:,:))
+  PRVHENI_MR(:,:,:)= 0.0
+  DO JL=1, INEGT
+    PRVHENI_MR(I1(JL), I2(JL), I3(JL)) = ZRVHENI_MR(JL)
+    PCIT      (I1(JL), I2(JL), I3(JL)) = ZCIT      (JL)
+  END DO
 END IF
 !
 DEALLOCATE(GLDCOMPUTE)
diff --git a/src/MNH/ice4_rsrimcg_old.f90 b/src/MNH/ice4_rsrimcg_old.f90
index 5d2471f7b..01f8f0310 100644
--- a/src/MNH/ice4_rsrimcg_old.f90
+++ b/src/MNH/ice4_rsrimcg_old.f90
@@ -43,6 +43,7 @@ SUBROUTINE ICE4_RSRIMCG_OLD(KSIZE, ODSOFT, ODCOMPUTE, &
 !!    -------------
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !
 !*      0. DECLARATIONS
@@ -76,7 +77,7 @@ INTEGER :: IGRIM, IGACC
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZVEC1, ZVEC2, ZVEC3
 INTEGER, DIMENSION(SIZE(PRHODREF)) :: IVEC1, IVEC2
 REAL, DIMENSION(SIZE(PRHODREF)) :: ZZW, ZZW2, ZZW6
-INTEGER :: JJ
+INTEGER :: JL
 !-------------------------------------------------------------------------------
 !
 !
@@ -87,14 +88,23 @@ INTEGER :: JJ
 PRSRIMCG_MR(:)=0.
 !
 IF(.NOT. ODSOFT) THEN
-  GRIM(:) = PRCT(:)>XRTMIN(2) .AND. PRST(:)>XRTMIN(5) .AND. ODCOMPUTE(:) .AND. PT(:)<XTT
-  IGRIM = COUNT(GRIM(:))
+  IGRIM = 0
+  GRIM(:) = .FALSE.
+  DO JL = 1, SIZE(GRIM)
+    IF ( PRCT(JL)>XRTMIN(2) .AND. PRST(JL)>XRTMIN(5) .AND. ODCOMPUTE(JL) .AND. PT(JL)<XTT ) THEN
+      IGRIM = IGRIM + 1
+      IVEC1(IGRIM) = Jl
+      GRIM(JL) = .TRUE.
+    END IF
+  END DO
   !
   IF(IGRIM>0 .AND. CSNOWRIMING=='OLD ') THEN
     !
     !        5.1.1  select the PLBDAS
     !
-    ZVEC1(1:IGRIM) = PACK( PLBDAS(:),MASK=GRIM(:) )
+    DO JL = 1, IGRIM
+      ZVEC1(JL) = PLBDAS(IVEC1(JL))
+    END DO
     !
     !        5.1.2  find the next lower indice for the PLBDAS in the geometrical
     !               set of Lbda_s used to tabulate some moments of the incomplete
@@ -111,7 +121,10 @@ IF(.NOT. ODSOFT) THEN
     !
     ZVEC1(1:IGRIM) =  XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                     - XGAMINC_RIM2( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
-    ZZW(:) = UNPACK( VECTOR=ZVEC1(1:IGRIM),MASK=GRIM,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JL = 1, IGRIM
+      ZZW(IVEC1(JL)) = ZVEC1(JL)
+    END DO
 
     !
     !        5.1.6  riming-conversion of the large sized aggregates into graupeln
diff --git a/src/MNH/ice4_tendencies.f90 b/src/MNH/ice4_tendencies.f90
index 5e536be5a..3259c42fe 100644
--- a/src/MNH/ice4_tendencies.f90
+++ b/src/MNH/ice4_tendencies.f90
@@ -8,7 +8,7 @@ INTERFACE
 SUBROUTINE ICE4_TENDENCIES(KSIZE, KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKL, &
                           &KRR, ODSOFT, ODCOMPUTE, &
                           &OWARM, HSUBG_RC_RR_ACCR, HSUBG_RR_EVAP, HSUBG_AUCV_RC, HSUBG_PR_PDF, &
-                          &PEXN, PRHODREF, PLVFACT, PLSFACT, LDMICRO, K1, K2, K3, &
+                          &PEXN, PRHODREF, PLVFACT, PLSFACT, K1, K2, K3, &
                           &PPRES, PCF, PSIGMA_RC, &
                           &PCIT, &
                           &PT, PTHT, &
@@ -40,7 +40,6 @@ REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PEXN
 REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRHODREF
 REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLVFACT
 REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLSFACT
-LOGICAL, DIMENSION(KIT,KJT,KKT), INTENT(IN) :: LDMICRO
 INTEGER, DIMENSION(KSIZE),    INTENT(IN)    :: K1
 INTEGER, DIMENSION(KSIZE),    INTENT(IN)    :: K2
 INTEGER, DIMENSION(KSIZE),    INTENT(IN)    :: K3
@@ -135,7 +134,7 @@ END MODULE MODI_ICE4_TENDENCIES
 SUBROUTINE ICE4_TENDENCIES(KSIZE, KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKL, &
                           &KRR, ODSOFT, ODCOMPUTE, &
                           &OWARM, HSUBG_RC_RR_ACCR, HSUBG_RR_EVAP, HSUBG_AUCV_RC, HSUBG_PR_PDF, &
-                          &PEXN, PRHODREF, PLVFACT, PLSFACT, LDMICRO, K1, K2, K3, &
+                          &PEXN, PRHODREF, PLVFACT, PLSFACT, K1, K2, K3, &
                           &PPRES, PCF, PSIGMA_RC, &
                           &PCIT, &
                           &PT, PTHT, &
@@ -164,7 +163,8 @@ SUBROUTINE ICE4_TENDENCIES(KSIZE, KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, K
 !!
 !!    MODIFICATIONS
 !!    -------------
-!!
+!
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !
 !*      0. DECLARATIONS
@@ -204,7 +204,6 @@ REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PEXN
 REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRHODREF
 REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLVFACT
 REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLSFACT
-LOGICAL, DIMENSION(KIT,KJT,KKT), INTENT(IN) :: LDMICRO
 INTEGER, DIMENSION(KSIZE),    INTENT(IN)    :: K1
 INTEGER, DIMENSION(KSIZE),    INTENT(IN)    :: K2
 INTEGER, DIMENSION(KSIZE),    INTENT(IN)    :: K3
@@ -401,8 +400,12 @@ IF(KSIZE>0) THEN
                         PRHODREF, ZRCT, PCF, PSIGMA_RC,&
                         PHLC_HCF, PHLC_LCF, PHLC_HRC, PHLC_LRC, ZRF)
   !Diagnostic of precipitation fraction
-  PRAINFR(:,:,:)=UNPACK(ZRF(:), MASK=LDMICRO(:,:,:), FIELD=0.)
-  ZRRT3D(:,:,:)=PRRT3D(:,:,:)-UNPACK(PRRHONG_MR(:), MASK=LDMICRO(:,:,:), FIELD=0.)
+  PRAINFR(:,:,:) = 0.
+  ZRRT3D (:,:,:) = PRRT3D(:,:,:)
+  DO JL=1,KSIZE
+    PRAINFR(K1(JL), K2(JL), K3(JL)) = ZRF(JL)
+    ZRRT3D (K1(JL), K2(JL), K3(JL)) = ZRRT3D(K1(JL), K2(JL), K3(JL)) - PRRHONG_MR(JL)
+  END DO
   CALL ICE4_RAINFR_VERT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKL, PRAINFR(:,:,:), ZRRT3D(:,:,:))
   DO JL=1,KSIZE
     ZRF(JL)=PRAINFR(K1(JL), K2(JL), K3(JL))
diff --git a/src/MNH/rain_ice.f90 b/src/MNH/rain_ice.f90
index b449c0a0c..04c6953c1 100644
--- a/src/MNH/rain_ice.f90
+++ b/src/MNH/rain_ice.f90
@@ -227,7 +227,6 @@ END MODULE MODI_RAIN_ICE
 !!                      land, sea and urban areas in the cloud sedimentation.
 !!      (D. Degrauwe), 2013-11: Export upper-air precipitation fluxes PFPR.
 !!      (S. Riette) Nov 2013 Protection against null sigma
-!!      Juan 24/09/2012: for BUG Pgi rewrite PACK function on mode_pack_pgi
 !!      (C. Lac) FIT temporal scheme : instant M removed
 !!      (JP Pinty), 01-2014 : ICE4 : partial reconversion of hail to graupel
 !!              July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
@@ -243,6 +242,7 @@ END MODULE MODI_RAIN_ICE
 !  P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !*       0.    DECLARATIONS
 !              ------------
@@ -258,9 +258,6 @@ use MODD_RAIN_ICE_DESCR, only: XLBEXR, XLBR, XRTMIN
 use MODD_RAIN_ICE_PARAM, only: XCRIAUTC
 
 use MODE_MSG
-#ifdef MNH_PGI
-USE MODE_PACK_PGI
-#endif
 use MODE_RAIN_ICE_FAST_RG,             only: RAIN_ICE_FAST_RG
 use MODE_RAIN_ICE_FAST_RH,             only: RAIN_ICE_FAST_RH
 use MODE_RAIN_ICE_FAST_RI,             only: RAIN_ICE_FAST_RI
@@ -580,7 +577,9 @@ IF( IMICRO >= 0 ) THEN
 !
   IF (LBU_ENABLE .OR. LLES_CALL) THEN
     ALLOCATE(ZRHODJ(IMICRO))
-    ZRHODJ(:) = PACK( PRHODJ(:,:,:),MASK=GMICRO(:,:,:) )
+    DO JL=1,IMICRO
+      ZRHODJ(JL) = PRHODJ(I1(JL),I2(JL),I3(JL))
+    END DO
   ELSE
     ALLOCATE(ZRHODJ(0))
   END IF
@@ -773,8 +772,10 @@ IF( IMICRO >= 0 ) THEN
   ENDIF
 
   !Diagnostic of precipitation fraction
-  ZW(:,:,:) = 0.
-  PRAINFR(:,:,:) = UNPACK( ZRF(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
+  PRAINFR(:,:,:) = 0.
+  DO JL=1,IMICRO
+    PRAINFR(I1(JL),I2(JL),I3(JL)) = ZRF(JL)
+  END DO
   CALL ICE4_RAINFR_VERT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKT, KKL, PRAINFR, PRRT(:,:,:))
   DO JL=1,IMICRO
     ZRF(JL)=PRAINFR(I1(JL),I2(JL),I3(JL))
@@ -809,7 +810,8 @@ IF( IMICRO >= 0 ) THEN
   IF( OWARM ) THEN    !  Check if the formation of the raindrops by the slow
                       !  warm processes is allowed
     PEVAP3D(:,:,:)= 0.
-    CALL RAIN_ICE_WARM(GMICRO, ZRHODREF, ZRVT, ZRCT, ZRRT, ZHLC_HCF, ZHLC_LCF, ZHLC_HRC, ZHLC_LRC,           &
+    CALL RAIN_ICE_WARM(GMICRO, IMICRO, I1, I2, I3,                                                           &
+                       ZRHODREF, ZRVT, ZRCT, ZRRT, ZHLC_HCF, ZHLC_LCF, ZHLC_HRC, ZHLC_LRC,                   &
                        ZRHODJ, ZPRES, ZZT, ZLBDAR, ZLBDAR_RF, ZLVFACT, ZCJ, ZKA, ZDV, ZRF, ZCF, ZTHT, ZTHLT, &
                        PRHODJ, PTHS, PRVS, ZRVS, ZRCS, ZRRS, ZTHS, ZUSW, PEVAP3D)
   END IF
@@ -861,29 +863,23 @@ IF( IMICRO >= 0 ) THEN
 !
 !
 !
-  ZW(:,:,:) = PRVS(:,:,:)
-  PRVS(:,:,:) = UNPACK( ZRVS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-  ZW(:,:,:) = PRCS(:,:,:)
-  PRCS(:,:,:) = UNPACK( ZRCS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-  ZW(:,:,:) = PRRS(:,:,:)
-  PRRS(:,:,:) = UNPACK( ZRRS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-  ZW(:,:,:) = PRIS(:,:,:)
-  PRIS(:,:,:) = UNPACK( ZRIS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-  ZW(:,:,:) = PRSS(:,:,:)
-  PRSS(:,:,:) = UNPACK( ZRSS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-  ZW(:,:,:) = PRGS(:,:,:)
-  PRGS(:,:,:) = UNPACK( ZRGS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
+  DO JL=1,IMICRO
+    PRVS(I1(JL),I2(JL),I3(JL)) = ZRVS(JL)
+    PRCS(I1(JL),I2(JL),I3(JL)) = ZRCS(JL)
+    PRRS(I1(JL),I2(JL),I3(JL)) = ZRRS(JL)
+    PRIS(I1(JL),I2(JL),I3(JL)) = ZRIS(JL)
+    PRSS(I1(JL),I2(JL),I3(JL)) = ZRSS(JL)
+    PRGS(I1(JL),I2(JL),I3(JL)) = ZRGS(JL)
+    PTHS(I1(JL),I2(JL),I3(JL)) = ZTHS(JL)
+    PCIT(I1(JL),I2(JL),I3(JL)) = ZCIT(JL)
+    !
+    PRAINFR(I1(JL),I2(JL),I3(JL)) = ZRF(JL)
+  END DO
   IF ( KRR == 7 ) THEN
-    ZW(:,:,:) = PRHS(:,:,:)
-    PRHS(:,:,:) = UNPACK( ZRHS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
+    DO JL=1,IMICRO
+      PRHS(I1(JL),I2(JL),I3(JL)) = ZRHS(JL)
+    END DO
   END IF
-  ZW(:,:,:) = PTHS(:,:,:)
-  PTHS(:,:,:) = UNPACK( ZTHS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-  ZW(:,:,:) = PCIT(:,:,:)
-  PCIT(:,:,:) = UNPACK( ZCIT(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!
-  ZW(:,:,:) = PRAINFR(:,:,:)
-  PRAINFR(:,:,:) = UNPACK( ZRF(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
 !
 !
 !
diff --git a/src/MNH/rain_ice_fast_rg.f90 b/src/MNH/rain_ice_fast_rg.f90
index af8226a36..41ce09cbf 100644
--- a/src/MNH/rain_ice_fast_rg.f90
+++ b/src/MNH/rain_ice_fast_rg.f90
@@ -6,6 +6,7 @@
 ! Modifications:
 !  P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !-----------------------------------------------------------------
 MODULE MODE_RAIN_ICE_FAST_RG
 
@@ -78,6 +79,7 @@ REAL,     DIMENSION(:),     intent(out)   :: PRWETG   ! Wet growth rate of the g
 !
 INTEGER                              :: IGDRY
 INTEGER                              :: JJ
+INTEGER, DIMENSION(size(PRHODREF))   :: I1
 INTEGER, DIMENSION(:), ALLOCATABLE   :: IVEC1, IVEC2      ! Vectors of indices for interpolations
 LOGICAL, DIMENSION(size(PRHODREF))   :: GDRY              ! Test where to compute dry growth
 REAL,    DIMENSION(size(PRHODREF))   :: ZZW               ! Work array
@@ -130,9 +132,17 @@ REAL,    DIMENSION(size(PRHODREF),7) :: ZZW1              ! Work arrays
 !
 !*       6.2.1  accretion of aggregates on the graupeln
 !
-  GDRY(:) = (PRST(:)>XRTMIN(5)) .AND. (PRGT(:)>XRTMIN(6)) .AND. (PRSS(:)>0.0)
-  IGDRY = COUNT( GDRY(:) )
-!
+  IGDRY = 0
+  DO JJ = 1, SIZE(GDRY)
+    IF ( PRST(JJ)>XRTMIN(5) .AND. PRGT(JJ)>XRTMIN(6) .AND. PRSS(JJ)>0.0 ) THEN
+      IGDRY = IGDRY + 1
+      I1(IGDRY) = JJ
+      GDRY(JJ) = .TRUE.
+    ELSE
+      GDRY(JJ) = .FALSE.
+    END IF
+  END DO
+
   IF( IGDRY>0 ) THEN
 !
 !*       6.2.2  allocations
@@ -145,8 +155,10 @@ REAL,    DIMENSION(size(PRHODREF),7) :: ZZW1              ! Work arrays
 !
 !*       6.2.3  select the (PLBDAG,PLBDAS) couplet
 !
-    ZVEC1(:) = PACK( PLBDAG(:),MASK=GDRY(:) )
-    ZVEC2(:) = PACK( PLBDAS(:),MASK=GDRY(:) )
+    DO JJ = 1, IGDRY
+      ZVEC1(JJ) = PLBDAG(I1(JJ))
+      ZVEC2(JJ) = PLBDAS(I1(JJ))
+    END DO
 !
 !*       6.2.4  find the next lower indice for the PLBDAG and for the PLBDAS
 !               in the geometrical set of (Lbda_g,Lbda_s) couplet use to
@@ -173,7 +185,10 @@ REAL,    DIMENSION(size(PRHODREF),7) :: ZZW1              ! Work arrays
                     - XKER_SDRYG(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                          * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGDRY
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
 !
     WHERE( GDRY(:) )
       ZZW1(:,3) = MIN( PRSS(:),XFSDRYG*ZZW(:)                         & ! RSDRYG
@@ -193,8 +208,16 @@ REAL,    DIMENSION(size(PRHODREF),7) :: ZZW1              ! Work arrays
 !
 !*       6.2.6  accretion of raindrops on the graupeln
 !
-  GDRY(:) = (PRRT(:)>XRTMIN(3)) .AND. (PRGT(:)>XRTMIN(6)) .AND. (PRRS(:)>0.0)
-  IGDRY = COUNT( GDRY(:) )
+  IGDRY = 0
+  DO JJ = 1, SIZE(GDRY)
+    IF ( PRRT(JJ)>XRTMIN(3) .AND. PRGT(JJ)>XRTMIN(6) .AND. PRRS(JJ)>0.0 ) THEN
+      IGDRY = IGDRY + 1
+      I1(IGDRY) = JJ
+      GDRY(JJ) = .TRUE.
+    ELSE
+      GDRY(JJ) = .FALSE.
+    END IF
+  END DO
 !
   IF( IGDRY>0 ) THEN
 !
@@ -208,8 +231,10 @@ REAL,    DIMENSION(size(PRHODREF),7) :: ZZW1              ! Work arrays
 !
 !*       6.2.8  select the (PLBDAG,PLBDAR) couplet
 !
-    ZVEC1(:) = PACK( PLBDAG(:),MASK=GDRY(:) )
-    ZVEC2(:) = PACK( PLBDAR(:),MASK=GDRY(:) )
+    DO JJ = 1, IGDRY
+      ZVEC1(JJ) = PLBDAG(I1(JJ))
+      ZVEC2(JJ) = PLBDAR(I1(JJ))
+    END DO
 !
 !*       6.2.9  find the next lower indice for the PLBDAG and for the PLBDAR
 !               in the geometrical set of (Lbda_g,Lbda_r) couplet use to
@@ -236,7 +261,10 @@ REAL,    DIMENSION(size(PRHODREF),7) :: ZZW1              ! Work arrays
                     - XKER_RDRYG(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                          * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGDRY
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
 !
     WHERE( GDRY(:) )
       ZZW1(:,4) = MIN( PRRS(:),XFRDRYG*ZZW(:)                    & ! RRDRYG
diff --git a/src/MNH/rain_ice_fast_rh.f90 b/src/MNH/rain_ice_fast_rh.f90
index 92be5f71c..d5f7813dd 100644
--- a/src/MNH/rain_ice_fast_rh.f90
+++ b/src/MNH/rain_ice_fast_rh.f90
@@ -6,6 +6,7 @@
 ! Modifications:
 !  P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !-----------------------------------------------------------------
 MODULE MODE_RAIN_ICE_FAST_RH
 
@@ -72,6 +73,7 @@ REAL,     DIMENSION(:),     intent(inout) :: PUSW     ! Undersaturation over wat
 !
 INTEGER                              :: IHAIL, IGWET
 INTEGER                              :: JJ
+INTEGER, DIMENSION(size(PRHODREF))   :: I1
 INTEGER, DIMENSION(:), ALLOCATABLE   :: IVEC1, IVEC2      ! Vectors of indices for interpolations
 LOGICAL, DIMENSION(size(PRHODREF))   :: GWET              ! Test where to compute wet growth
 LOGICAL, DIMENSION(size(PRHODREF))   :: GHAIL             ! Test where to compute hail growth
@@ -81,8 +83,16 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
 !
 !-------------------------------------------------------------------------------
 !
-  GHAIL(:) = PRHT(:)>XRTMIN(7)
-  IHAIL = COUNT(GHAIL(:))
+  IHAIL = 0
+  DO JJ = 1, SIZE(GHAIL)
+    IF ( PRHT(JJ)>XRTMIN(7) ) THEN
+      IHAIL = IHAIL + 1
+      I1(IHAIL) = JJ
+      GHAIL(JJ) = .TRUE.
+    ELSE
+      GHAIL(JJ) = .FALSE.
+    END IF
+  END DO
 !
   IF( IHAIL>0 ) THEN
 !
@@ -104,8 +114,16 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
 !
 !*       7.2.1  accretion of aggregates on the hailstones
 !
-    GWET(:) = GHAIL(:) .AND. (PRST(:)>XRTMIN(5) .AND. PRSS(:)>0.0)
-    IGWET = COUNT( GWET(:) )
+    IGWET = 0
+    DO JJ = 1, SIZE(GWET)
+      IF ( GHAIL(JJ) .AND. PRST(JJ)>XRTMIN(5) .AND. PRSS(JJ)>0.0 ) THEN
+        IGWET = IGWET + 1
+        I1(IGWET) = JJ
+        GWET(JJ) = .TRUE.
+      ELSE
+        GWET(JJ) = .FALSE.
+      END IF
+    END DO
 !
     IF( IGWET>0 ) THEN
 !
@@ -119,8 +137,10 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
 !
 !*       7.2.3  select the (PLBDAH,PLBDAS) couplet
 !
-      ZVEC1(:) = PACK( PLBDAH(:),MASK=GWET(:) )
-      ZVEC2(:) = PACK( PLBDAS(:),MASK=GWET(:) )
+      DO JJ = 1, IGWET
+        ZVEC1(JJ) = PLBDAH(I1(JJ))
+        ZVEC2(JJ) = PLBDAS(I1(JJ))
+      END DO
 !
 !*       7.2.4  find the next lower indice for the PLBDAG and for the PLBDAS
 !               in the geometrical set of (Lbda_h,Lbda_s) couplet use to
@@ -147,7 +167,10 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
                      - XKER_SWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                           * (ZVEC1(JJ) - 1.0)
       END DO
-      ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
+      ZZW(:) = 0.
+      DO JJ = 1, IGWET
+        ZZW(I1(JJ)) = ZVEC3(JJ)
+      END DO
 !
       WHERE( GWET(:) )
         ZZW1(:,3) = MIN( PRSS(:),XFSWETH*ZZW(:)                       & ! RSWETH
@@ -166,8 +189,16 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
 !
 !*       7.2.6  accretion of graupeln on the hailstones
 !
-    GWET(:) = GHAIL(:) .AND. (PRGT(:)>XRTMIN(6) .AND. PRGS(:)>0.0)
-    IGWET = COUNT( GWET(:) )
+    IGWET = 0
+    DO JJ = 1, SIZE(GWET)
+      IF ( GHAIL(JJ) .AND. PRGT(JJ)>XRTMIN(6) .AND. PRGS(JJ)>0.0 ) THEN
+        IGWET = IGWET + 1
+        I1(IGWET) = JJ
+        GWET(JJ) = .TRUE.
+      ELSE
+        GWET(JJ) = .FALSE.
+      END IF
+    END DO
 !
     IF( IGWET>0 ) THEN
 !
@@ -181,8 +212,10 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
 !
 !*       7.2.8  select the (PLBDAH,PLBDAG) couplet
 !
-      ZVEC1(:) = PACK( PLBDAH(:),MASK=GWET(:) )
-      ZVEC2(:) = PACK( PLBDAG(:),MASK=GWET(:) )
+      DO JJ = 1, IGWET
+        ZVEC1(JJ) = PLBDAH(I1(JJ))
+        ZVEC2(JJ) = PLBDAG(I1(JJ))
+      END DO
 !
 !*       7.2.9  find the next lower indice for the PLBDAH and for the PLBDAG
 !               in the geometrical set of (Lbda_h,Lbda_g) couplet use to
@@ -209,7 +242,10 @@ REAL,    DIMENSION(size(PRHODREF),6) :: ZZW1              ! Work arrays
                      - XKER_GWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                           * (ZVEC1(JJ) - 1.0)
       END DO
-      ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
+      ZZW(:) = 0.
+      DO JJ = 1, IGWET
+        ZZW(I1(JJ)) = ZVEC3(JJ)
+      END DO
 !
       WHERE( GWET(:) )
         ZZW1(:,5) = MAX(MIN( PRGS(:),XFGWETH*ZZW(:)                       & ! RGWETH
diff --git a/src/MNH/rain_ice_fast_rs.f90 b/src/MNH/rain_ice_fast_rs.f90
index 72d7c02bc..7f90da8c0 100644
--- a/src/MNH/rain_ice_fast_rs.f90
+++ b/src/MNH/rain_ice_fast_rs.f90
@@ -6,6 +6,7 @@
 ! Modifications:
 !  P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !-----------------------------------------------------------------
 MODULE MODE_RAIN_ICE_FAST_RS
 
@@ -69,6 +70,7 @@ REAL,     DIMENSION(:),     INTENT(INOUT) :: PTHS     ! Theta source
 !
 INTEGER                              :: IGRIM, IGACC
 INTEGER                              :: JJ
+INTEGER, DIMENSION(size(PRHODREF))   :: I1
 INTEGER, DIMENSION(:), ALLOCATABLE   :: IVEC1, IVEC2      ! Vectors of indices for interpolations
 LOGICAL, DIMENSION(size(PRHODREF))   :: GRIM              ! Test where to compute riming
 LOGICAL, DIMENSION(size(PRHODREF))   :: GACC              ! Test where to compute accretion
@@ -81,11 +83,17 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
 !
   ZZW1(:,:) = 0.0
 !
-! GRIM(:) = (PRCT(:)>0.0) .AND. (PRST(:)>0.0) .AND.            &
-  GRIM(:) = (PRCT(:)>XRTMIN(2)) .AND. (PRST(:)>XRTMIN(5)) .AND.            &
-                                (PRCS(:)>0.0) .AND. (PZT(:)<XTT)
-  IGRIM = COUNT( GRIM(:) )
-!
+  IGRIM = 0
+  DO JJ = 1, SIZE(GRIM)
+    IF (PRCT(JJ)>XRTMIN(2) .AND. PRST(JJ)>XRTMIN(5) .AND. PRCS(JJ)>0.0 .AND. PZT(JJ)<XTT ) THEN
+      IGRIM = IGRIM + 1
+      I1(IGRIM) = JJ
+      GRIM(JJ) = .TRUE.
+    ELSE
+      GRIM(JJ) = .FALSE.
+    END IF
+  END DO
+  !
   IF( IGRIM>0 ) THEN
 !
 !        5.1.0  allocations
@@ -96,7 +104,9 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
 !
 !        5.1.1  select the PLBDAS
 !
-    ZVEC1(:) = PACK( PLBDAS(:),MASK=GRIM(:) )
+    DO JJ = 1, IGRIM
+      ZVEC1(JJ) = PLBDAS(I1(JJ))
+    END DO
 !
 !        5.1.2  find the next lower indice for the PLBDAS in the geometrical
 !               set of Lbda_s used to tabulate some moments of the incomplete
@@ -112,7 +122,10 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
 !
     ZVEC1(1:IGRIM) =   XGAMINC_RIM1( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                      - XGAMINC_RIM1( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
-    ZZW(:) = UNPACK( VECTOR=ZVEC1(:),MASK=GRIM,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGRIM
+      ZZW(I1(JJ)) = ZVEC1(JJ)
+    END DO
 !
 !        5.1.4  riming of the small sized aggregates
 !
@@ -131,7 +144,10 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
 !
     ZVEC1(1:IGRIM) =  XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                     - XGAMINC_RIM2( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
-    ZZW(:) = UNPACK( VECTOR=ZVEC1(:),MASK=GRIM,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGRIM
+      ZZW(I1(JJ)) = ZVEC1(JJ)
+    END DO
 !
 !        5.1.6  riming-conversion of the large sized aggregates into graupeln
 !
@@ -170,10 +186,18 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
 !*       5.2    rain accretion onto the aggregates
 !
   ZZW1(:,2:3) = 0.0
-   GACC(:) = (PRRT(:)>XRTMIN(3)) .AND. (PRST(:)>XRTMIN(5)) .AND.            &
-                            (PRRS(:)>0.0) .AND. (PZT(:)<XTT)
-  IGACC = COUNT( GACC(:) )
-!
+
+  IGACC = 0
+  DO JJ = 1, SIZE(GACC)
+    IF (PRRT(JJ)>XRTMIN(3) .AND. PRST(JJ)>XRTMIN(5) .AND. PRRS(JJ)>0.0 .AND. PZT(JJ)<XTT ) THEN
+      IGACC = IGACC + 1
+      I1(IGACC) = JJ
+      GACC(JJ) = .TRUE.
+    ELSE
+      GACC(JJ) = .FALSE.
+    END IF
+  END DO
+  !
   IF( IGACC>0 ) THEN
 !
 !        5.2.0  allocations
@@ -186,8 +210,10 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
 !
 !        5.2.1  select the (PLBDAS,PLBDAR) couplet
 !
-    ZVEC1(:) = PACK( PLBDAS(:),MASK=GACC(:) )
-    ZVEC2(:) = PACK( PLBDAR(:),MASK=GACC(:) )
+    DO JJ = 1, IGACC
+      ZVEC1(JJ) = PLBDAS(I1(JJ))
+      ZVEC2(JJ) = PLBDAR(I1(JJ))
+    END DO
 !
 !        5.2.2  find the next lower indice for the PLBDAS and for the PLBDAR
 !               in the geometrical set of (Lbda_s,Lbda_r) couplet use to
@@ -214,7 +240,10 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
                     - XKER_RACCSS(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                           * (ZVEC1(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGACC
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
 !
 !        5.2.4  raindrop accretion on the small sized aggregates
 !
@@ -241,7 +270,9 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
                     -  XKER_RACCS(IVEC2(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
                                                            * (ZVEC2(JJ) - 1.0)
     END DO
-    ZZW1(:,2) = ZZW1(:,2)*UNPACK( VECTOR=ZVEC3(:),MASK=GACC(:),FIELD=0.0 )
+    DO JJ = 1, IGACC
+      ZZW1(I1(JJ), 2) =  ZZW1(I1(JJ), 2 ) * ZVEC3(JJ)
+    END DO
                                                                        !! RRACCS!
 !        5.2.5  perform the bilinear interpolation of the normalized
 !               SACCRG-kernel
@@ -254,7 +285,10 @@ REAL,    DIMENSION(size(PRHODREF),4) :: ZZW1              ! Work arrays
                     - XKER_SACCRG(IVEC2(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
                                                           * (ZVEC2(JJ) - 1.0)
     END DO
-    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
+    ZZW(:) = 0.
+    DO JJ = 1, IGACC
+      ZZW(I1(JJ)) = ZVEC3(JJ)
+    END DO
 !
 !        5.2.6  raindrop accretion-conversion of the large sized aggregates
 !               into graupeln
diff --git a/src/MNH/rain_ice_nucleation.f90 b/src/MNH/rain_ice_nucleation.f90
index 08b7a107e..959c7bb0d 100644
--- a/src/MNH/rain_ice_nucleation.f90
+++ b/src/MNH/rain_ice_nucleation.f90
@@ -6,6 +6,7 @@
 ! Modifications:
 !  P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !-----------------------------------------------------------------
 MODULE MODE_RAIN_ICE_NUCLEATION
 
@@ -131,7 +132,10 @@ IF( INEGT >= 1 ) THEN
 !*       3.1.2   update the r_i and r_v mixing ratios
 !
     ZZW(:) = MIN( ZZW(:),50.E3 ) ! limitation provisoire a 50 l^-1
-    ZW(:,:,:) = UNPACK( ZZW(:),MASK=GNEGT(:,:,:),FIELD=0.0 )
+    ZW(:,:,:) = 0.0
+    DO JL=1, INEGT
+      ZW(I1(JL), I2(JL), I3(JL)) = ZZW( JL )
+    END DO
     ZW(:,:,:) = MAX( ZW(:,:,:) ,0.0 ) *XMNU0/(PRHODREF(:,:,:)*PTSTEP)
     PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:)
     PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:)
@@ -146,8 +150,10 @@ IF( INEGT >= 1 ) THEN
     END IF
                                  ! f(L_s*(RVHENI))
     ZZW(:) = MAX( ZZW(:)+ZCIT(:),ZCIT(:) )
-    PCIT(:,:,:) = MAX( UNPACK( ZZW(:),MASK=GNEGT(:,:,:),FIELD=0.0 ) , &
-                       PCIT(:,:,:) )
+    PCIT(:,:,:) = MAX( PCIT(:,:,:), 0.0 )
+    DO JL=1, INEGT
+      PCIT(I1(JL), I2(JL), I3(JL)) = MAX( ZZW( JL ), PCIT(I1(JL), I2(JL), I3(JL)), 0.0 )
+    END DO
   END IF
   DEALLOCATE(ZSSI)
   DEALLOCATE(ZUSW)
diff --git a/src/MNH/rain_ice_red.f90 b/src/MNH/rain_ice_red.f90
index 50b6ef1ff..321d59013 100644
--- a/src/MNH/rain_ice_red.f90
+++ b/src/MNH/rain_ice_red.f90
@@ -228,7 +228,6 @@ END MODULE MODI_RAIN_ICE_RED
 !!                      land, sea and urban areas in the cloud sedimentation.
 !!      (D. Degrauwe), 2013-11: Export upper-air precipitation fluxes PFPR.
 !!      (S. Riette) Nov 2013 Protection against null sigma
-!!      Juan 24/09/2012: for BUG Pgi rewrite PACK function on mode_pack_pgi
 !!      (C. Lac) FIT temporal scheme : instant M removed
 !!      (JP Pinty), 01-2014 : ICE4 : partial reconversion of hail to graupel
 !!              July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
@@ -242,6 +241,7 @@ END MODULE MODI_RAIN_ICE_RED
 !!                  02/2019 C.Lac add rain fraction as an output field
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
+!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !
 !*       0.    DECLARATIONS
 !              ------------
@@ -256,9 +256,6 @@ USE MODD_VAR_ll,         ONLY: IP
 
 USE MODE_ll
 USE MODE_MSG
-#ifdef MNH_PGI
-USE MODE_PACK_PGI
-#endif
 use mode_tools,          only: Countjv
 
 USE MODI_BUDGET
@@ -767,7 +764,7 @@ DO WHILE(ANY(ZTIME(:)<PTSTEP)) ! Loop to *really* compute tendencies
     CALL ICE4_TENDENCIES(IMICRO, IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKT, KKL, &
                         &KRR, LSOFT, LLCOMPUTE, &
                         &OWARM, CSUBG_RC_RR_ACCR, CSUBG_RR_EVAP, HSUBG_AUCV_RC, CSUBG_PR_PDF, &
-                        &ZEXN, ZRHODREF, ZLVFACT, ZLSFACT, ODMICRO, I1, I2, I3, &
+                        &ZEXN, ZRHODREF, ZLVFACT, ZLSFACT, I1, I2, I3, &
                         &ZPRES, ZCF, ZSIGMA_RC, &
                         &ZCIT, &
                         &ZZT, ZTHT, &
@@ -1019,15 +1016,17 @@ ENDDO
 !               ---------------------
 !
 IF(IMICRO>0) THEN
-  ZW(:,:,:) = 0.
-  ZHLC_HCF3D(:,:,:) = UNPACK(ZHLC_HCF(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))
-  ZW(:,:,:) = 0.
-  ZHLC_LCF3D(:,:,:) = UNPACK(ZHLC_LCF(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))
-  ZW(:,:,:) = 0.
-  ZHLC_HRC3D(:,:,:) = UNPACK(ZHLC_HRC(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))
-  ZW(:,:,:) = 0.
-  ZHLC_LRC3D(:,:,:) = UNPACK(ZHLC_LRC(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))
-  PCIT(:,:,:) = UNPACK(ZCIT(:), MASK=ODMICRO(:,:,:), FIELD=PCIT(:,:,:))
+  ZHLC_HCF3D(:,:,:)=0.
+  ZHLC_LCF3D(:,:,:)=0.
+  ZHLC_HRC3D(:,:,:)=0.
+  ZHLC_LRC3D(:,:,:)=0.
+  DO JL=1,IMICRO
+    ZHLC_HCF3D(I1(JL), I2(JL), I3(JL)) = ZHLC_HCF(JL)
+    ZHLC_LCF3D(I1(JL), I2(JL), I3(JL)) = ZHLC_LCF(JL)
+    ZHLC_HRC3D(I1(JL), I2(JL), I3(JL)) = ZHLC_HRC(JL)
+    ZHLC_LRC3D(I1(JL), I2(JL), I3(JL)) = ZHLC_LRC(JL)
+    PCIT(I1(JL), I2(JL), I3(JL)) = ZCIT(JL)
+  END DO
 ELSE
   PRAINFR(:,:,:)=0.
   ZHLC_HCF3D(:,:,:)=0.
@@ -1037,7 +1036,10 @@ ELSE
   PCIT(:,:,:) = 0.
 ENDIF
 IF(OWARM) THEN
-  PEVAP3D(:,:,:)=UNPACK(ZRREVAV(:), MASK=ODMICRO(:,:,:), FIELD=0.)
+  PEVAP3D(:,:,:) = 0.
+  DO JL=1,IMICRO
+    PEVAP3D(I1(JL), I2(JL), I3(JL)) = ZRREVAV(JL)
+  END DO
 ENDIF
 !
 !
@@ -1076,17 +1078,26 @@ IF(GEXT_TEND) THEN
   ZTHT(:) = ZTHT(:) - ZEXT_TH(:) * PTSTEP
 ENDIF
 !Tendencies computed from difference between old state and new state (can be negative)
-ZW_RVS(:,:,:) = (UNPACK(ZRVT(:), MASK=ODMICRO(:,:,:), FIELD=PRVT(:,:,:)) - PRVT(:,:,:))*ZINV_TSTEP
-ZW_RCS(:,:,:) = (UNPACK(ZRCT(:), MASK=ODMICRO(:,:,:), FIELD=PRCT(:,:,:)) - PRCT(:,:,:))*ZINV_TSTEP
-ZW_RRS(:,:,:) = (UNPACK(ZRRT(:), MASK=ODMICRO(:,:,:), FIELD=PRRT(:,:,:)) - PRRT(:,:,:))*ZINV_TSTEP
-ZW_RIS(:,:,:) = (UNPACK(ZRIT(:), MASK=ODMICRO(:,:,:), FIELD=PRIT(:,:,:)) - PRIT(:,:,:))*ZINV_TSTEP
-ZW_RSS(:,:,:) = (UNPACK(ZRST(:), MASK=ODMICRO(:,:,:), FIELD=PRST(:,:,:)) - PRST(:,:,:))*ZINV_TSTEP
-ZW_RGS(:,:,:) = (UNPACK(ZRGT(:), MASK=ODMICRO(:,:,:), FIELD=PRGT(:,:,:)) - PRGT(:,:,:))*ZINV_TSTEP
-IF(KRR==7) THEN
-  ZW_RHS(:,:,:) = (UNPACK(ZRHT(:), MASK=ODMICRO(:,:,:), FIELD=PRHT(:,:,:)) - PRHT(:,:,:))*ZINV_TSTEP
-ELSE
+  ZW_RVS(:,:,:) = 0.
+  ZW_RCS(:,:,:) = 0.
+  ZW_RRS(:,:,:) = 0.
+  ZW_RIS(:,:,:) = 0.
+  ZW_RSS(:,:,:) = 0.
+  ZW_RGS(:,:,:) = 0.
   ZW_RHS(:,:,:) = 0.
-ENDIF
+  DO JL=1,IMICRO
+    ZW_RVS(I1(JL), I2(JL), I3(JL)) = ( ZRVT(JL) - PRVT(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+    ZW_RCS(I1(JL), I2(JL), I3(JL)) = ( ZRCT(JL) - PRCT(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+    ZW_RRS(I1(JL), I2(JL), I3(JL)) = ( ZRRT(JL) - PRRT(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+    ZW_RIS(I1(JL), I2(JL), I3(JL)) = ( ZRIT(JL) - PRIT(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+    ZW_RSS(I1(JL), I2(JL), I3(JL)) = ( ZRST(JL) - PRST(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+    ZW_RGS(I1(JL), I2(JL), I3(JL)) = ( ZRGT(JL) - PRGT(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+  END DO
+  IF(KRR==7) THEN
+  DO JL=1,IMICRO
+    ZW_RHS(I1(JL), I2(JL), I3(JL)) = ( ZRHT(JL) - PRHT(I1(JL), I2(JL), I3(JL)) ) * ZINV_TSTEP
+  END DO
+END IF
 ZW_THS(:,:,:) = (ZW_RCS(:,:,:)+ZW_RRS(:,:,:)                            )*ZZ_LVFACT(:,:,:) + &
               & (ZW_RIS(:,:,:)+ZW_RSS(:,:,:)+ZW_RGS(:,:,:)+ZW_RHS(:,:,:))*ZZ_LSFACT(:,:,:)
 !We apply these tendencies to the S variables
@@ -1107,7 +1118,9 @@ CALL CORRECT_NEGATIVITIES(KRR, ZW_RVS, ZW_RCS, ZW_RRS, &
 !
 IF(LBU_ENABLE) THEN
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RVHENI(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RVHENI(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:)
   PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*ZZ_LSFACT(:,:,:)
@@ -1116,7 +1129,9 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'HENU_BU_RRI')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCHONI(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCHONI(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:)
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
@@ -1125,7 +1140,9 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'HON_BU_RRI')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RRHONG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRHONG(JL) * ZINV_TSTEP
+  END DO
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
@@ -1134,7 +1151,9 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'SFR_BU_RRG')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RVDEPS(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RVDEPS(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:)
   PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*ZZ_LSFACT(:,:,:)
@@ -1143,21 +1162,27 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'DEPS_BU_RRS')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RIAGGS(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIAGGS(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:)
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'AGGS_BU_RRI')
   IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'AGGS_BU_RRS')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RIAUTS(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIAUTS(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:)
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'AUTS_BU_RRI')
   IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'AUTS_BU_RRS')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RVDEPG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RVDEPG(JL) * ZINV_TSTEP
+  END DO
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*ZZ_LSFACT(:,:,:)
@@ -1167,21 +1192,27 @@ IF(LBU_ENABLE) THEN
 
   IF(OWARM) THEN
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RCAUTR(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCAUTR(JL) * ZINV_TSTEP
+    END DO
     PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
     PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
     IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'AUTO_BU_RRC')
     IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'AUTO_BU_RRR')
 
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RCACCR(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCACCR(JL) * ZINV_TSTEP
+    END DO
     PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
     PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
     IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'ACCR_BU_RRC')
     IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'ACCR_BU_RRR')
 
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RREVAV(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RREVAV(JL) * ZINV_TSTEP
+    END DO
     PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
     PRVS(:,:,:) = PRVS(:,:,:) + ZW(:,:,:)
     PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*ZZ_LVFACT(:,:,:)
@@ -1191,17 +1222,23 @@ IF(LBU_ENABLE) THEN
   ENDIF
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCRIMSS(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCRIMSS(JL) * ZINV_TSTEP
+  END DO
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCRIMSG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCRIMSG(JL) * ZINV_TSTEP
+  END DO
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RSRIMCG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSRIMCG(JL) * ZINV_TSTEP
+  END DO
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
   IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'RIM_BU_RTH')
@@ -1210,17 +1247,23 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'RIM_BU_RRG')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RRACCSS(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRACCSS(JL) * ZINV_TSTEP
+  END DO
   PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
   PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RRACCSG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRACCSG(JL) * ZINV_TSTEP
+  END DO
   PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RSACCRG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSACCRG(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'ACC_BU_RTH')
@@ -1229,11 +1272,15 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'ACC_BU_RRG')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RSMLTG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSMLTG(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCMLTSR, MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCMLTSR(JL) * ZINV_TSTEP
+  END DO
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
   IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'CMEL_BU_RRS')
@@ -1242,16 +1289,22 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'CMEL_BU_RRR')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RICFRRG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RICFRRG(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RRCFRIG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRCFRIG(JL) * ZINV_TSTEP
+  END DO
   PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RICFRR(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RICFRR(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
   IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'CFRZ_BU_RTH')
@@ -1260,21 +1313,29 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'CFRZ_BU_RRG')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCWETG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCWETG(JL) * ZINV_TSTEP
+  END DO
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RRWETG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRWETG(JL) * ZINV_TSTEP
+  END DO
   PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RIWETG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIWETG(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RSWETG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSWETG(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'WETG_BU_RTH')
@@ -1286,7 +1347,9 @@ IF(LBU_ENABLE) THEN
 
   IF(KRR==7) THEN
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RWETGH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RWETGH(JL) * ZINV_TSTEP
+    END DO
     PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'GHCV_BU_RRG')
@@ -1294,21 +1357,29 @@ IF(LBU_ENABLE) THEN
   END IF
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCDRYG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCDRYG(JL) * ZINV_TSTEP
+  END DO
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RRDRYG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRDRYG(JL) * ZINV_TSTEP
+  END DO
   PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RIDRYG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIDRYG(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RSDRYG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSDRYG(JL) * ZINV_TSTEP
+  END DO
   PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
   IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'DRYG_BU_RTH')
@@ -1319,7 +1390,9 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'DRYG_BU_RRG')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RGMLTR(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RGMLTR(JL) * ZINV_TSTEP
+  END DO
   PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
   PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
@@ -1329,25 +1402,35 @@ IF(LBU_ENABLE) THEN
 
   IF(KRR==7) THEN
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RCWETH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCWETH(JL) * ZINV_TSTEP
+    END DO
     PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RRWETH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRWETH(JL) * ZINV_TSTEP
+    END DO
     PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RIWETH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIWETH(JL) * ZINV_TSTEP
+    END DO
     PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RSWETH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSWETH(JL) * ZINV_TSTEP
+    END DO
     PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RGWETH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RGWETH(JL) * ZINV_TSTEP
+    END DO
     PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'WETH_BU_RTH')
@@ -1358,36 +1441,50 @@ IF(LBU_ENABLE) THEN
     IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'WETH_BU_RRH')
 
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RGWETH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RGWETH(JL) * ZINV_TSTEP
+    END DO
     PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'HGCV_BU_RRG')
     IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'HGCV_BU_RRH')
 
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RCDRYH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCDRYH(JL) * ZINV_TSTEP
+    END DO
     PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RRDRYH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RRDRYH(JL) * ZINV_TSTEP
+    END DO
     PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RIDRYH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIDRYH(JL) * ZINV_TSTEP
+    END DO
     PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RSDRYH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RSDRYH(JL) * ZINV_TSTEP
+    END DO
     PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RGDRYH(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RGDRYH(JL) * ZINV_TSTEP
+    END DO
     PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:)
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RDRYHG(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RDRYHG(JL) * ZINV_TSTEP
+    END DO
     PRHS(:,:,:) = PRHS(:,:,:) - ZW(:,:,:)
     PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:)
     IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'DRYH_BU_RTH')
@@ -1399,7 +1496,9 @@ IF(LBU_ENABLE) THEN
     IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'DRYH_BU_RRH')
 
     ZW(:,:,:) = 0.
-    ZW(:,:,:)=UNPACK(ZTOT_RHMLTR(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+    DO JL=1,IMICRO
+      ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RHMLTR(JL) * ZINV_TSTEP
+    END DO
     PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
     PRHS(:,:,:) = PRHS(:,:,:) - ZW(:,:,:)
     PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
@@ -1409,7 +1508,9 @@ IF(LBU_ENABLE) THEN
   ENDIF
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RIMLTC(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RIMLTC(JL) * ZINV_TSTEP
+  END DO
   PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:)
   PRCS(:,:,:) = PRCS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
@@ -1418,7 +1519,9 @@ IF(LBU_ENABLE) THEN
   IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'IMLT_BU_RRI')
 
   ZW(:,:,:) = 0.
-  ZW(:,:,:)=UNPACK(ZTOT_RCBERI(:), MASK=ODMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP
+  DO JL=1,IMICRO
+    ZW(I1(JL), I2(JL), I3(JL)) = ZTOT_RCBERI(JL) * ZINV_TSTEP
+  END DO
   PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
   PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:)
   PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:))
diff --git a/src/MNH/rain_ice_warm.f90 b/src/MNH/rain_ice_warm.f90
index d83fad363..54a8c315c 100644
--- a/src/MNH/rain_ice_warm.f90
+++ b/src/MNH/rain_ice_warm.f90
@@ -5,6 +5,7 @@
 !-----------------------------------------------------------------
 ! Modifications:
 !  P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
+!  P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
 !-----------------------------------------------------------------
 MODULE MODE_RAIN_ICE_WARM
 
@@ -16,7 +17,8 @@ MODULE MODE_RAIN_ICE_WARM
 
 CONTAINS
 
-SUBROUTINE RAIN_ICE_WARM(OMICRO, PRHODREF, PRVT, PRCT, PRRT, PHLC_HCF, PHLC_LCF, PHLC_HRC, PHLC_LRC,           &
+SUBROUTINE RAIN_ICE_WARM(OMICRO, KMICRO, K1, K2, K3,                                                           &
+                         PRHODREF, PRVT, PRCT, PRRT, PHLC_HCF, PHLC_LCF, PHLC_HRC, PHLC_LRC,                   &
                          PRHODJ, PPRES, PZT, PLBDAR, PLBDAR_RF, PLVFACT, PCJ, PKA, PDV, PRF, PCF, PTHT, PTHLT, &
                          PRHODJ3D, PTHS3D, PRVS3D, PRVS, PRCS, PRRS, PTHS, PUSW, PEVAP3D)
 !
@@ -38,6 +40,10 @@ IMPLICIT NONE
 !*       0.1   Declarations of dummy arguments :
 !
 LOGICAL,  DIMENSION(:,:,:), intent(in)    :: OMICRO   ! Test where to compute all processes
+INTEGER,                    intent(in)    :: KMICRO
+INTEGER,  DIMENSION(:),     intent(in)    :: K1
+INTEGER,  DIMENSION(:),     intent(in)    :: K2
+INTEGER,  DIMENSION(:),     intent(in)    :: K3
 REAL,     DIMENSION(:),     intent(in)    :: PRHODREF ! RHO Dry REFerence
 REAL,     DIMENSION(:),     intent(in)    :: PRVT     ! Water vapor m.r. at t
 REAL,     DIMENSION(:),     intent(in)    :: PRCT     ! Cloud water m.r. at t
@@ -73,6 +79,7 @@ REAL,     DIMENSION(:,:,:), INTENT(INOUT) :: PEVAP3D  ! Rain evap profile
 !
 !*       0.2  declaration of local variables
 !
+INTEGER                         :: JL
 REAL, DIMENSION(size(PRHODREF)) :: ZZW  ! Work array
 REAL, DIMENSION(size(PRHODREF)) :: ZZW2 ! Work array
 REAL, DIMENSION(size(PRHODREF)) :: ZZW3 ! Work array
@@ -230,7 +237,10 @@ REAL, DIMENSION(size(PRHODREF)) :: ZZW4 ! Work array
     IF (LBUDGET_RR) CALL BUDGET (                                               &
                      UNPACK(PRRS(:)*PRHODJ(:),MASK=OMICRO(:,:,:),FIELD=0.0),    &
                                                               8,'REVA_BU_RRR')
-    PEVAP3D(:,:,:)=UNPACK(ZZW(:),MASK=OMICRO(:,:,:),FIELD=PEVAP3D(:,:,:))
+
+    DO JL = 1, KMICRO
+      PEVAP3D(K1(JL), K2(JL), K3(JL)) = ZZW( JL )
+    END DO
 !
   END SUBROUTINE RAIN_ICE_WARM
 
-- 
GitLab