From 85a538e554505f767862e3edca5603ae9b3fdd9e Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Tue, 14 Feb 2017 09:44:40 +0100
Subject: [PATCH] Philippe 10-20/02/2017: OpenACC: work on rain_ice.f90

---
 src/MNH/rain_ice.f90 | 705 +++++++++++++++++++++++++++++--------------
 1 file changed, 486 insertions(+), 219 deletions(-)

diff --git a/src/MNH/rain_ice.f90 b/src/MNH/rain_ice.f90
index 5a0c4b83c..9e01f3d9b 100644
--- a/src/MNH/rain_ice.f90
+++ b/src/MNH/rain_ice.f90
@@ -332,14 +332,14 @@ LOGICAL,                OPTIONAL, INTENT(IN)    :: OCONVHG ! Switch for conversi
 !
 ! IN variables
 !
-!$acc declare present(PDZZ,PRHODJ,PRHODREF,PEXNREF,PPABST, &
-!$acc &               PCLDFR,PTHT,PRVT,PRCT,PRRT,PRIT,PRST,PRGT, &
-!$acc &               PSIGS)
+!$acc declare present(PDZZ,PRHODJ,PRHODREF,PEXNREF,PPABST,               &
+!$acc &               PCLDFR,PTHT,PRVT,PRCT,PRRT,PRIT,PRST,PRGT,         &
+!$acc &               PSIGS)                                             &
 !
 ! INOUT variables
 !
-!$acc declare present(PCIT,PTHS,PRVS,PRCS,PRRS,PRIS,PRSS,PRGS,&
-!$acc &               PINPRC,PINPRR,PINPRR3D,PEVAP3D,PINPRS,PINPRG,PINDEP)
+!$acc &       present(PCIT,PTHS,PRVS,PRCS,PRRS,PRIS,PRSS,PRGS,           &
+!$acc &               PINPRC,PINDEP,PINPRR,PINPRR3D,PEVAP3D,PINPRS,PINPRG)
 
 !
 ! OUT variables
@@ -379,6 +379,7 @@ LOGICAL, DIMENSION(:), ALLOCATABLE :: GACC ! Test where to compute accretion
 LOGICAL, DIMENSION(:), ALLOCATABLE :: GDRY ! Test where to compute dry growth
 LOGICAL, DIMENSION(:), ALLOCATABLE :: GWET  ! Test where to compute wet growth
 LOGICAL, DIMENSION(:), ALLOCATABLE :: GHAIL ! Test where to compute hail growth
+LOGICAL, DIMENSION(:), ALLOCATABLE :: GWORK
 LOGICAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2)):: GDEP
 INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC1,IVEC2       ! Vectors of indices for
         ! interpolations
@@ -465,10 +466,19 @@ REAL, DIMENSION(SIZE(XRTMIN))     :: ZRTMIN
 !
 INTEGER , DIMENSION(SIZE(GMICRO)) :: I1,I2,I3 ! Used to replace the COUNT
 INTEGER                           :: JL       ! and PACK intrinsics
-!$acc declare create(GNEGT,GMICRO,&
-!$acc &              ZW,ZT, &
-!$acc &              ZZW) &
-!$acc &       device_resident(I1,I2,I3)
+!$acc declare create(GSEDIMR,GSEDIMC,GSEDIMI,GSEDIMS,GSEDIMG,GSEDIMH,             &
+!$acc &              GNEGT,GMICRO,GRIM,GACC,GDRY,GWET,GHAIL,GWORK,GDEP,           &
+!$acc &              IVEC1,IVEC2,ZVEC1,ZVEC2,ZVEC3,ZW,                            &
+!$acc &              ZPRCS,ZPRRS,ZPRSS,ZPRGS,ZPRHS,                               &
+!$acc &              ZWSED,ZWSEDW1,ZWSEDW2,ZCONC_TMP,ZT,ZRAY,ZLBC,ZFSEDC,         &
+!$acc &              ZRVT,ZRCT,ZRRT,ZRIT,ZRST,ZRGT,ZRHT,ZCIT,                     &
+!$acc &              ZRVS,ZRCS,ZRRS,ZRIS,ZRSS,ZRGS,ZRHS,ZTHS,                     &
+!$acc &              ZRHODREF,ZRHODREFR,ZRHODREFI,ZRHODREFS,ZRHODREFG,ZRHODREFH,  &
+!$acc &              ZRHODJ,ZZT,ZPRES,ZZW,ZLSFACT,ZLVFACT,                        &
+!$acc &              ZUSW,ZSSI,ZLBDAR,ZLBDAS,ZLBDAG,ZLBDAH,ZRDRYG,ZRWETG,         &
+!$acc &              ZAI,ZCJ,ZKA,ZDV,ZZW1,ZRTMIN)                                 &
+!$acc &       device_resident(ZCRIAUTI,ZEXNREF,ZSIGMA_RC,ZCF,ZCC,ZFSEDC1D,ZWLBDC, &
+!$acc &                       ZCONC,ZRAY1D,ZWLBDA,I1,I2,I3)
 !
 #ifdef _OPENACC
 PRINT *,'OPENACC: RAIN_ICE being implemented'
@@ -528,9 +538,9 @@ GMICRO(:,:,:) = .FALSE.
 !$acc update self(GMICRO)
 
 #ifndef _OPENACC
-IMICRO = COUNTJV( GMICRO(:,:,:),I1(:),I2(:),I3(:))
+IMICRO = COUNTJV3D( GMICRO(:,:,:),I1(:),I2(:),I3(:))
 #else
-CALL COUNTJV_DEVICE(GMICRO(:,:,:),I1(:),I2(:),I3(:),IMICRO)
+CALL COUNTJV3D_DEVICE(GMICRO(:,:,:),I1(:),I2(:),I3(:),IMICRO)
 #endif
 
 IF( IMICRO >= 0 ) THEN
@@ -574,6 +584,7 @@ IF( IMICRO >= 0 ) THEN
   ALLOCATE(ZCJ(IMICRO))
   ALLOCATE(ZKA(IMICRO))
   ALLOCATE(ZDV(IMICRO))
+  ALLOCATE(GWORK(IMICRO))
 !
   IF ( KRR == 7 ) THEN
     ALLOCATE(ZZW1(IMICRO,7))
@@ -627,23 +638,30 @@ IF( IMICRO >= 0 ) THEN
   ZLSFACT(:) = (XLSTT+(XCPV-XCI)*(ZZT(:)-XTT))/ZZW(:) ! L_s/(Pi_ref*C_ph)
   ZLVFACT(:) = (XLVTT+(XCPV-XCL)*(ZZT(:)-XTT))/ZZW(:) ! L_v/(Pi_ref*C_ph)
 
-!$acc end kernels
-!$acc update self(ZZW)
-!acc kernels
 #ifndef MNH_BITREP
   ZZW(:) = EXP( XALPI - XBETAI/ZZT(:) - XGAMI*ALOG(ZZT(:) ) )
 #else
   ZZW(:) = BR_EXP( XALPI - XBETAI/ZZT(:) - XGAMI*BR_LOG(ZZT(:) ) )
 #endif
-!$acc update device(ZZW)
-!$acc kernels
   ZSSI(:) = ZRVT(:)*( ZPRES(:)-ZZW(:) ) / ( (XMV/XMD) * ZZW(:) ) - 1.0
                                                     ! Supersaturation over ice
 !$acc end kernels
-!
+!$acc update self(ZRHODREF,ZZT,ZPRES,ZLSFACT,ZLVFACT,ZSSI,ZRCT,ZRRT,ZRIT,ZRST,ZRGT,ZCIT,&
+!$acc &           ZRCS,ZRRS,ZRIS,ZRSS,ZRGS,ZTHS)
+#ifdef _OPENACC
+  IF (KRR == 7 ) THEN
+!$acc update self(ZRHT,ZRHS)
+  END IF
+#endif
+  !
   IF (LBU_ENABLE .OR. LLES_CALL .OR. LCHECK) THEN
     ALLOCATE(ZRHODJ(IMICRO))
-    ZRHODJ(:) = PACK( PRHODJ(:,:,:),MASK=GMICRO(:,:,:) )
+!$acc kernels
+    DO JL=1,IMICRO
+      ZRHODJ(JL) = PRHODJ(I1(JL),I2(JL),I3(JL))
+    END DO
+!$acc end kernels
+!$acc update self(ZRHODJ)
   END IF
 !
   CALL RAIN_ICE_SLOW
@@ -656,6 +674,7 @@ IF( IMICRO >= 0 ) THEN
 !
 !*       3.1    compute the slope parameter Lbda_r
 !
+!$acc kernels
   WHERE( ZRRT(:)>0.0 )
 #ifndef MNH_BITREP
     ZLBDAR(:)  = XLBR*( ZRHODREF(:)*MAX( ZRRT(:),XRTMIN(3) ) )**XLBEXR
@@ -663,10 +682,11 @@ IF( IMICRO >= 0 ) THEN
     ZLBDAR(:)  = XLBR*BR_POW( ZRHODREF(:)*MAX( ZRRT(:),XRTMIN(3) ) ,XLBEXR)
 #endif
   END WHERE
+!$acc end kernels
+!$acc update self(ZLBDAR)
 !
   IF( OWARM ) THEN    !  Check if the formation of the raindrops by the slow
                       !  warm processes is allowed
-    PEVAP3D(:,:,:)= 0.
     CALL RAIN_ICE_WARM
   END IF
 !
@@ -709,39 +729,58 @@ IF( IMICRO >= 0 ) THEN
 !
 !
 !
+#ifndef _OPENACC
   ZW(:,:,:) = PRVS(:,:,:)
   PRVS(:,:,:) = UNPACK( ZRVS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRVS)
   ZW(:,:,:) = PRCS(:,:,:)
   PRCS(:,:,:) = UNPACK( ZRCS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRCS)
   ZW(:,:,:) = PRRS(:,:,:)
   PRRS(:,:,:) = UNPACK( ZRRS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRRS)
   ZW(:,:,:) = PRIS(:,:,:)
   PRIS(:,:,:) = UNPACK( ZRIS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRIS)
   ZW(:,:,:) = PRSS(:,:,:)
   PRSS(:,:,:) = UNPACK( ZRSS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRSS)
   ZW(:,:,:) = PRGS(:,:,:)
   PRGS(:,:,:) = UNPACK( ZRGS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRGS)
   IF ( KRR == 7 ) THEN
     ZW(:,:,:) = PRHS(:,:,:)
     PRHS(:,:,:) = UNPACK( ZRHS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PRHS)
   END IF
   ZW(:,:,:) = PTHS(:,:,:)
   PTHS(:,:,:) = UNPACK( ZTHS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PTHS)
   ZW(:,:,:) = PCIT(:,:,:)
   PCIT(:,:,:) = UNPACK( ZCIT(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-!$acc update device(PCIT)
+#else
+!$acc update device(ZRCS,ZRRS,ZRIS,ZRRS,ZRSS,ZRGS,ZTHS)
+!$acc kernels
+  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)
+  END DO
+!$acc end kernels
+!$acc update self(PRVS,PRCS,PRRS,PRIS,PRSS,PRGS,PTHS,PCIT)
+  IF ( KRR == 7 ) THEN
+!$acc update device(ZRHS)
+!$acc kernels
+    DO JL=1,IMICRO
+      PRHS(I1(JL),I2(JL),I3(JL)) = ZRHS(JL)
+    END DO
+!$acc end kernels
+!$acc update self(PRHS)
+  END IF
+#endif  
+print *,'PW: RAIN_ICE 12'
 !
 !
 !
   DEALLOCATE(ZZW1)
+  DEALLOCATE(GWORK)
   DEALLOCATE(ZDV)
   DEALLOCATE(ZCJ)
   DEALLOCATE(ZRDRYG)
@@ -969,7 +1008,7 @@ ENDIF
 !
    IF (OSEDIC) THEN
 #ifdef _OPENACC
-    PRINT *,'OPENACC: RAIN_ICE_SEDIMENTATION_SPLIT::OSDIC=.T. not yet implemented'
+    PRINT *,'OPENACC: RAIN_ICE_SEDIMENTATION_SPLIT::OSEDIC=.T. not yet implemented'
     CALL ABORT
 #endif
     ZRAY(:,:,:)   = 0.
@@ -1017,6 +1056,15 @@ ILENALLOCS = 0
 ILENALLOCG = 0
 IF ( KRR == 7 ) ILENALLOCH = 0
 !$acc end kernels
+!$acc update self(GSEDIMR,GSEDIMI,GSEDIMS,GSEDIMG,ZRTMIN)
+#ifdef _OPENACC
+IF (OSEDIC) THEN
+!$acc update self(GSEDIMC)
+END IF
+IF (KRR==7) THEN
+!$acc update self(GSEDIMH)
+END IF
+#endif
 !
 ! ZPiS = Specie i source creating during the current time step
 ! PRiS = Source of the previous time step
@@ -1074,86 +1122,88 @@ DO JN = 1 , KSPLITR
  IF ( KRR == 7 ) GSEDIMH(IIB:IIE,IJB:IJE,IKTB:IKTE) =            &
                   PRHS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(7)
 !
- IF (OSEDIC) ISEDIMC = COUNTJV( GSEDIMC(:,:,:),IC1(:),IC2(:),IC3(:))
- ISEDIMR = COUNTJV( GSEDIMR(:,:,:),IR1(:),IR2(:),IR3(:))
- ISEDIMI = COUNTJV( GSEDIMI(:,:,:),II1(:),II2(:),II3(:))
- ISEDIMS = COUNTJV( GSEDIMS(:,:,:),IS1(:),IS2(:),IS3(:))
- ISEDIMG = COUNTJV( GSEDIMG(:,:,:),IG1(:),IG2(:),IG3(:))
- IF ( KRR == 7 ) ISEDIMH = COUNTJV( GSEDIMH(:,:,:),IH1(:),IH2(:),IH3(:))
+ IF (OSEDIC) ISEDIMC = COUNTJV3D( GSEDIMC(:,:,:),IC1(:),IC2(:),IC3(:))
+ ISEDIMR = COUNTJV3D( GSEDIMR(:,:,:),IR1(:),IR2(:),IR3(:))
+ ISEDIMI = COUNTJV3D( GSEDIMI(:,:,:),II1(:),II2(:),II3(:))
+ ISEDIMS = COUNTJV3D( GSEDIMS(:,:,:),IS1(:),IS2(:),IS3(:))
+ ISEDIMG = COUNTJV3D( GSEDIMG(:,:,:),IG1(:),IG2(:),IG3(:))
+ IF ( KRR == 7 ) ISEDIMH = COUNTJV3D( GSEDIMH(:,:,:),IH1(:),IH2(:),IH3(:))
 !
 !*       2.1   for cloud
 !
- IF (OSEDIC) THEN
-  ZWSED(:,:,:) = 0.
-  IF( JN==1 ) PRCS(:,:,:) = PRCS(:,:,:) * PTSTEP
-  IF( ISEDIMC >= 1 ) THEN
-    IF ( ISEDIMC .GT. ILENALLOCC ) THEN
-      IF ( ILENALLOCC .GT. 0 ) THEN
-        DEALLOCATE (ZRCS, ZRHODREFC, ILISTC,ZWLBDC,ZCONC,ZRCT,  &
-                    ZZT,ZPRES,ZRAY1D,ZFSEDC1D,ZWLBDA,ZCC )
+  IF (OSEDIC) THEN
+    ZWSED(:,:,:) = 0.
+    IF( JN==1 ) PRCS(:,:,:) = PRCS(:,:,:) * PTSTEP
+    IF( ISEDIMC >= 1 ) THEN
+      IF ( ISEDIMC .GT. ILENALLOCC ) THEN
+        IF ( ILENALLOCC .GT. 0 ) THEN
+          DEALLOCATE (ZRCS, ZRHODREFC, ILISTC,ZWLBDC,ZCONC,ZRCT,  &
+                      ZZT,ZPRES,ZRAY1D,ZFSEDC1D,ZWLBDA,ZCC )
+        END IF
+        ILENALLOCC = MAX (IOLDALLOCC, 2*ISEDIMC )
+        IOLDALLOCC = ILENALLOCC
+        ALLOCATE(ZRCS(ILENALLOCC), ZRHODREFC(ILENALLOCC), ILISTC(ILENALLOCC), &
+          ZWLBDC(ILENALLOCC), ZCONC(ILENALLOCC), ZRCT(ILENALLOCC), ZZT(ILENALLOCC), &
+          ZPRES(ILENALLOCC), ZRAY1D(ILENALLOCC), ZFSEDC1D(ILENALLOCC), &
+          ZWLBDA(ILENALLOCC), ZCC(ILENALLOCC)  )
       END IF
-      ILENALLOCC = MAX (IOLDALLOCC, 2*ISEDIMC )
-      IOLDALLOCC = ILENALLOCC
-      ALLOCATE(ZRCS(ILENALLOCC), ZRHODREFC(ILENALLOCC), ILISTC(ILENALLOCC), &
-        ZWLBDC(ILENALLOCC), ZCONC(ILENALLOCC), ZRCT(ILENALLOCC), ZZT(ILENALLOCC), &
-        ZPRES(ILENALLOCC), ZRAY1D(ILENALLOCC), ZFSEDC1D(ILENALLOCC), &
-        ZWLBDA(ILENALLOCC), ZCC(ILENALLOCC)  )
-    END IF
 !
-    DO JL=1,ISEDIMC
-      ZRCS(JL) = PRCS(IC1(JL),IC2(JL),IC3(JL))
-      ZRHODREFC(JL) =  PRHODREF(IC1(JL),IC2(JL),IC3(JL))
-      ZWLBDC(JL) = ZLBC(IC1(JL),IC2(JL),IC3(JL))
-      ZCONC(JL) = ZCONC3D(IC1(JL),IC2(JL),IC3(JL))
-      ZRCT(JL) = PRCT(IC1(JL),IC2(JL),IC3(JL))
-      ZZT(JL) = PTHT(IC1(JL),IC2(JL),IC3(JL))
-      ZPRES(JL) = PPABST(IC1(JL),IC2(JL),IC3(JL))
-      ZRAY1D(JL) = ZRAY(IC1(JL),IC2(JL),IC3(JL))
-      ZFSEDC1D(JL) = ZFSEDC(IC1(JL),IC2(JL),IC3(JL))
-    END DO
+      DO JL=1,ISEDIMC
+        ZRCS(JL) = PRCS(IC1(JL),IC2(JL),IC3(JL))
+        ZRHODREFC(JL) =  PRHODREF(IC1(JL),IC2(JL),IC3(JL))
+        ZWLBDC(JL) = ZLBC(IC1(JL),IC2(JL),IC3(JL))
+        ZCONC(JL) = ZCONC3D(IC1(JL),IC2(JL),IC3(JL))
+        ZRCT(JL) = PRCT(IC1(JL),IC2(JL),IC3(JL))
+        ZZT(JL) = PTHT(IC1(JL),IC2(JL),IC3(JL))
+        ZPRES(JL) = PPABST(IC1(JL),IC2(JL),IC3(JL))
+        ZRAY1D(JL) = ZRAY(IC1(JL),IC2(JL),IC3(JL))
+        ZFSEDC1D(JL) = ZFSEDC(IC1(JL),IC2(JL),IC3(JL))
+      END DO
 !
-    ILISTLENC = 0
-    DO JL=1,ISEDIMC
-     IF( ZRCS(JL) .GT. ZRTMIN(2) ) THEN
-       ILISTLENC = ILISTLENC + 1
-       ILISTC(ILISTLENC) = JL
-     END IF
-    END DO
-       DO JJ = 1, ILISTLENC
-          JL = ILISTC(JJ)
-          IF (ZRCS(JL) .GT. ZRTMIN(2) .AND. ZRCT(JL) .GT. XRTMIN(2)) THEN
-            ZWLBDC(JL) = ZWLBDC(JL) * ZCONC(JL) / (ZRHODREFC(JL) * ZRCT(JL))
+      ILISTLENC = 0
+      DO JL=1,ISEDIMC
+        IF( ZRCS(JL) .GT. ZRTMIN(2) ) THEN
+          ILISTLENC = ILISTLENC + 1
+          ILISTC(ILISTLENC) = JL
+        END IF
+      END DO
+!
+      DO JJ = 1, ILISTLENC
+        JL = ILISTC(JJ)
+        IF (ZRCS(JL) .GT. ZRTMIN(2) .AND. ZRCT(JL) .GT. XRTMIN(2)) THEN
+          ZWLBDC(JL) = ZWLBDC(JL) * ZCONC(JL) / (ZRHODREFC(JL) * ZRCT(JL))
 #ifndef MNH_BITREP
-            ZWLBDC(JL) = ZWLBDC(JL)**XLBEXC
+          ZWLBDC(JL) = ZWLBDC(JL)**XLBEXC
 #else
-            ZWLBDC(JL) = BR_POW(ZWLBDC(JL),XLBEXC)
+          ZWLBDC(JL) = BR_POW(ZWLBDC(JL),XLBEXC)
 #endif
-            ZRAY1D(JL) = ZRAY1D(JL) / ZWLBDC(JL) !! ZRAY : mean diameter=M(1)/2
+          ZRAY1D(JL) = ZRAY1D(JL) / ZWLBDC(JL) !! ZRAY : mean diameter=M(1)/2
 #ifndef MNH_BITREP
-            ZZT(JL)    = ZZT(JL) * (ZPRES(JL)/XP00)**(XRD/XCPD)
+          ZZT(JL)    = ZZT(JL) * (ZPRES(JL)/XP00)**(XRD/XCPD)
 #else
-            ZZT(JL)    = ZZT(JL) * BR_POW(ZPRES(JL)/XP00,XRD/XCPD)
+          ZZT(JL)    = ZZT(JL) * BR_POW(ZPRES(JL)/XP00,XRD/XCPD)
 #endif
-            ZWLBDA(JL) = 6.6E-8*(101325./ZPRES(JL))*(ZZT(JL)/293.15)
-            ZCC(JL)    = XCC*(1.+1.26*ZWLBDA(JL)/ZRAY1D(JL)) !! XCC modified for cloud
+          ZWLBDA(JL) = 6.6E-8*(101325./ZPRES(JL))*(ZZT(JL)/293.15)
+          ZCC(JL)    = XCC*(1.+1.26*ZWLBDA(JL)/ZRAY1D(JL)) !! XCC modified for cloud
 #ifndef MNH_BITREP
-            ZWSED (IC1(JL),IC2(JL),IC3(JL))= ZRHODREFC(JL)**(-XCEXVT +1 ) *   &
+          ZWSED (IC1(JL),IC2(JL),IC3(JL))= ZRHODREFC(JL)**(-XCEXVT +1 ) *   &
               ZWLBDC(JL)**(-XDC)*ZCC(JL)*ZFSEDC1D(JL) * ZRCS(JL)
 #else
-            ZWSED (IC1(JL),IC2(JL),IC3(JL))= BR_POW(ZRHODREFC(JL),-XCEXVT +1) *   &
+          ZWSED (IC1(JL),IC2(JL),IC3(JL))= BR_POW(ZRHODREFC(JL),-XCEXVT +1) *   &
               BR_POW(ZWLBDC(JL),-XDC)*ZCC(JL)*ZFSEDC1D(JL) * ZRCS(JL)
 #endif
-          END IF
-       END DO
+        END IF
+      END DO
+!$acc update device(ZZT,ZPRES)
+    END IF
+    DO JK = IKTB , IKTE
+      PRCS(:,:,JK) = PRCS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
+    END DO
+    PINPRC(:,:) = PINPRC(:,:) + ZWSED(:,:,IKB) / XRHOLW / KSPLITR 
+    IF( JN==KSPLITR ) THEN
+      PRCS(:,:,:) = PRCS(:,:,:) * ZINVTSTEP
+    END IF
   END IF
-       DO JK = IKTB , IKTE
-         PRCS(:,:,JK) = PRCS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
-       END DO
-      PINPRC(:,:) = PINPRC(:,:) + ZWSED(:,:,IKB) / XRHOLW / KSPLITR 
-      IF( JN==KSPLITR ) THEN
-        PRCS(:,:,:) = PRCS(:,:,:) * ZINVTSTEP
-      END IF
- END IF
 !
 !*       2.2   for rain
 !
@@ -1169,6 +1219,8 @@ DO JN = 1 , KSPLITR
       ALLOCATE(ZRRS(ILENALLOCR), ZRHODREFR(ILENALLOCR), ILISTR(ILENALLOCR))
     END IF
 !
+!$acc update device(PRRS)
+!$acc kernels
     DO JL=1,ISEDIMR
       ZRRS(JL) = PRRS(IR1(JL),IR2(JL),IR3(JL))
       ZRHODREFR(JL) =  PRHODREF(IR1(JL),IR2(JL),IR3(JL))
@@ -1176,30 +1228,34 @@ DO JN = 1 , KSPLITR
 !
     ILISTLENR = 0
     DO JL=1,ISEDIMR
-     IF( ZRRS(JL) .GT. ZRTMIN(3) ) THEN
-       ILISTLENR = ILISTLENR + 1
-       ILISTR(ILISTLENR) = JL
-     END IF
+      IF( ZRRS(JL) .GT. ZRTMIN(3) ) THEN
+        ILISTLENR = ILISTLENR + 1
+        ILISTR(ILISTLENR) = JL
+      END IF
     END DO
-       DO JJ = 1, ILISTLENR
-          JL = ILISTR(JJ)
+!$acc end kernels
+!$acc update self(ZRHODREFR,ZRRS)
+!acc kernels
+    DO JJ = 1, ILISTLENR
+      JL = ILISTR(JJ)
 #ifndef MNH_BITREP
-           ZWSED (IR1(JL),IR2(JL),IR3(JL))= XFSEDR  * ZRRS(JL)**XEXSEDR *   &
-                                        ZRHODREFR(JL)**(XEXSEDR-XCEXVT)
+      ZWSED (IR1(JL),IR2(JL),IR3(JL))= XFSEDR  * ZRRS(JL)**XEXSEDR *   &
+                                       ZRHODREFR(JL)**(XEXSEDR-XCEXVT)
 #else
-           ZWSED (IR1(JL),IR2(JL),IR3(JL))= XFSEDR  * BR_POW(ZRRS(JL),XEXSEDR) *   &
-                                        BR_POW(ZRHODREFR(JL),XEXSEDR-XCEXVT)
+      ZWSED (IR1(JL),IR2(JL),IR3(JL))= XFSEDR  * BR_POW(ZRRS(JL),XEXSEDR) *   &
+                                       BR_POW(ZRHODREFR(JL),XEXSEDR-XCEXVT)
 #endif
-       END DO
+    END DO
+!acc end kernels
+  END IF
+  DO JK = IKTB , IKTE
+    PRRS(:,:,JK) = PRRS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
+  END DO
+  PINPRR(:,:) = PINPRR(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
+  PINPRR3D(:,:,:) = PINPRR3D(:,:,:) + ZWSED(:,:,1:IKT)/XRHOLW/KSPLITR 
+  IF ( JN==KSPLITR ) THEN
+    PRRS(:,:,:) = PRRS(:,:,:) * ZINVTSTEP
   END IF
-       DO JK = IKTB , IKTE
-         PRRS(:,:,JK) = PRRS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
-       END DO
-       PINPRR(:,:) = PINPRR(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
-       PINPRR3D(:,:,:) = PINPRR3D(:,:,:) + ZWSED(:,:,1:IKT)/XRHOLW/KSPLITR 
-      IF( JN==KSPLITR ) THEN
-        PRRS(:,:,:) = PRRS(:,:,:) * ZINVTSTEP
-      END IF
 !
 !*       2.3   for pristine ice
 !
@@ -1491,6 +1547,7 @@ ENDIF
       
 
 ZRTMIN(:)    = XRTMIN(:) * ZINVTSTEP
+!$acc update device(ZRTMIN)
 !
 IF (OSEDIC) THEN
   ZPRCS(:,:,:) = 0.0
@@ -1534,7 +1591,7 @@ END DO
        !estimation of q' taking into account incomming ZWSED
        ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)
 
-       JCOUNT=COUNTJV2((PRCS(:,:,JK) > ZRTMIN(2) .AND. PRCT(:,:,JK) > ZRTMIN(2)) .OR. &
+       JCOUNT=COUNTJV2D((PRCS(:,:,JK) > ZRTMIN(2) .AND. PRCT(:,:,JK) > ZRTMIN(2)) .OR. &
                        (ZQP(:,:) > ZRTMIN(2)),I1(:),I2(:))
        DO JL=1, JCOUNT
          JI=I1(JL)
@@ -1620,7 +1677,7 @@ END DO
      !estimation of q' taking into account incomming ZWSED
      ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)
 
-     JCOUNT=COUNTJV2((PRRS(:,:,JK) > ZRTMIN(3)) .OR. &
+     JCOUNT=COUNTJV2D((PRRS(:,:,JK) > ZRTMIN(3)) .OR. &
                      (ZQP(:,:) > ZRTMIN(3)),I1(:),I2(:))
      DO JL=1, JCOUNT
        JI=I1(JL)
@@ -1682,7 +1739,7 @@ END DO
      !estimation of q' taking into account incomming ZWSED
      ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)
 
-     JCOUNT=COUNTJV2((PRIS(:,:,JK) > MAX(ZRTMIN(4),1.0E-7 )) .OR. &
+     JCOUNT=COUNTJV2D((PRIS(:,:,JK) > MAX(ZRTMIN(4),1.0E-7 )) .OR. &
                      (ZQP(:,:) > MAX(ZRTMIN(4),1.0E-7 )),I1(:),I2(:))
      DO JL=1, JCOUNT
        JI=I1(JL)
@@ -1753,7 +1810,7 @@ END DO
      !estimation of q' taking into account incomming ZWSED
      ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)
 
-     JCOUNT=COUNTJV2((PRSS(:,:,JK) > ZRTMIN(5)) .OR. &
+     JCOUNT=COUNTJV2D((PRSS(:,:,JK) > ZRTMIN(5)) .OR. &
                      (ZQP(:,:) > ZRTMIN(5)),I1(:),I2(:))
      DO JL=1, JCOUNT
        JI=I1(JL)
@@ -1816,7 +1873,7 @@ END DO
      !estimation of q' taking into account incomming ZWSED
      ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)
 
-     JCOUNT=COUNTJV2((PRGS(:,:,JK) > ZRTMIN(6)) .OR. &
+     JCOUNT=COUNTJV2D((PRGS(:,:,JK) > ZRTMIN(6)) .OR. &
                      (ZQP(:,:) > ZRTMIN(6)),I1(:),I2(:))
      DO JL=1, JCOUNT
        JI=I1(JL)
@@ -1877,7 +1934,7 @@ END DO
      !estimation of q' taking into account incomming ZWSED
      ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)
 
-     JCOUNT=COUNTJV2((PRHS(:,:,JK)+ZQP(JI,JJ) > ZRTMIN(7)) .OR. &
+     JCOUNT=COUNTJV2D((PRHS(:,:,JK)+ZQP(JI,JJ) > ZRTMIN(7)) .OR. &
                      (ZQP(:,:) > ZRTMIN(7)),I1(:),I2(:))
      DO JL=1, JCOUNT
        JI=I1(JL)
@@ -1998,9 +2055,9 @@ GNEGT(:,:,:) = .FALSE.
 GNEGT(IIB:IIE,IJB:IJE,IKTB:IKTE) = ZT(IIB:IIE,IJB:IJE,IKTB:IKTE)<XTT
 !$acc end kernels
 #ifndef _OPENACC
-INEGT = COUNTJV( GNEGT(:,:,:),I1(:),I2(:),I3(:))
+INEGT = COUNTJV3D( GNEGT(:,:,:),I1(:),I2(:),I3(:))
 #else
-CALL COUNTJV_DEVICE(GNEGT(:,:,:),I1(:),I2(:),I3(:),INEGT)
+CALL COUNTJV3D_DEVICE(GNEGT(:,:,:),I1(:),I2(:),I3(:),INEGT)
 #endif
 
 IF( INEGT >= 1 ) THEN
@@ -2060,19 +2117,17 @@ IF( INEGT >= 1 ) THEN
   ZZW(:) = ZZW(:) - ZCIT(:)
   ZZWMAX = MAXVAL(ZZW(:))
 !$acc end kernels
+!$acc update self(ZUSW,ZSSI)
   IF( ZZWMAX > 0.0 ) THEN
 !
 !*       3.1.2   update the r_i and r_v mixing ratios
 !
 !$acc kernels
-    ZZW(:) = MIN( ZZW(:),50.E3 ) ! limitation provisoire a 50 l^-1
-!$acc end kernels
-#ifndef _OPENACC
-    ZW(:,:,:) = UNPACK( ZZW(:),MASK=GNEGT(:,:,:),FIELD=0.0 )
-#else
-    CALL UNPACK_DEVICE (ZZW,GNEGT,0.0,ZW)
-#endif
-!$acc kernels
+  ZZW(:) = MIN( ZZW(:),50.E3 ) ! limitation provisoire a 50 l^-1
+  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(:,:,:)
@@ -2094,7 +2149,12 @@ IF( INEGT >= 1 ) THEN
     PCIT(:,:,:) = MAX( UNPACK( ZZW(:),MASK=GNEGT(:,:,:),FIELD=0.0 ) , &
                        PCIT(:,:,:) )
 #else
-    CALL UNPACK_DEVICE(ZZW,GNEGT,0.0,ZW)
+!$acc kernels
+  ZW(:,:,:) = 0.0
+  DO JL=1,INEGT
+    ZW(I1(JL),I2(JL),I3(JL)) = ZZW(JL)
+  END DO
+!$acc end kernels
 !$acc kernels
     PCIT(:,:,:) = MAX( ZW(:,:,:), PCIT(:,:,:) )
 !$acc end kernels
@@ -2129,14 +2189,11 @@ IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'HENU_BU_RRI')
 !
 IMPLICIT NONE
 !
-LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
-!$acc declare create(GWORK)
 !-------------------------------------------------------------------------------
 !
 !
 !*       3.2     compute the homogeneous nucleation source: RCHONI
 !
-  ALLOCATE(GWORK(IMICRO))
 !
 !$acc kernels
   ZZW(:) = 0.0
@@ -2153,6 +2210,7 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
     ZTHS(:) = ZTHS(:) + ZZW(:)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RCHONI))
   ENDWHERE
 !$acc end kernels
+!$acc update self(ZTHS,ZRCS,ZRIS)
 !
   IF (LBUDGET_TH) CALL BUDGET (                                                &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),  &
@@ -2176,6 +2234,7 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
     ZTHS(:) = ZTHS(:) + ZZW(:)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RRHONG))
   ENDWHERE
 !$acc end kernels
+!$acc update self(ZTHS,ZRRS,ZRGS)
 !
   IF (LBUDGET_TH) CALL BUDGET (                                                &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),  &
@@ -2189,14 +2248,13 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
 !
 !*       3.4    compute the deposition, aggregation and autoconversion sources
 !
-!acc kernels
+!$acc kernels
   ZKA(:) = 2.38E-2 + 0.0071E-2 * ( ZZT(:) - XTT )          ! k_a
 #ifndef MNH_BITREP
   ZDV(:) = 0.211E-4 * (ZZT(:)/XTT)**1.94 * (XP00/ZPRES(:)) ! D_v
 #else
   ZDV(:) = 0.211E-4 * BR_POW(ZZT(:)/XTT,1.94) * (XP00/ZPRES(:)) ! D_v
 #endif
-!acc end kernels
 !
 !*       3.4.1  compute the thermodynamical function A_i(T,P)
 !*              and the c^prime_j (in the ventilation factor)
@@ -2237,7 +2295,8 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
 #endif
   END WHERE
   ZZW(:) = 0.0
-  WHERE ( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>0.0) )
+  GWORK = (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>0.0)
+  WHERE ( GWORK(:) )
 #ifndef MNH_BITREP
     ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) *                               &
              ( X0DEPS*ZLBDAS(:)**XEX0DEPS + X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS )
@@ -2251,6 +2310,9 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
     ZRVS(:) = ZRVS(:) - ZZW(:)
     ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
   END WHERE
+!$acc end kernels
+!$acc update self(ZAI,ZCJ,ZKA,ZDV)
+!$acc update self(ZLBDAS,ZTHS,ZRVS,ZRSS)
   IF (LBUDGET_TH) CALL BUDGET (                                                &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),  &
                                                               4,'DEPS_BU_RTH')
@@ -2263,8 +2325,10 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
 !
 !*       3.4.4  compute the aggregation on r_s: RIAGGS
 !
+!$acc kernels present(ZLBDAS,ZRHODREF)
   ZZW(:) = 0.0
-  WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>0.0) )
+  GWORK(:) = (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>0.0)
+  WHERE ( GWORK(:) )
 #ifndef MNH_BITREP
     ZZW(:) = MIN( ZRIS(:),XFIAGGS * EXP( XCOLEXIS*(ZZT(:)-XTT) ) &
                    * ZRIT(:)                      &
@@ -2279,6 +2343,8 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
     ZRSS(:)  = ZRSS(:)  + ZZW(:)
     ZRIS(:)  = ZRIS(:)  - ZZW(:)
   END WHERE
+!$acc end kernels
+!$acc update self(ZRIS,ZRSS)
   IF (LBUDGET_RI) CALL BUDGET (                                                 &
                      UNPACK(ZRIS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
                                                               9,'AGGS_BU_RRI')
@@ -2289,13 +2355,15 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
 !*       3.4.5  compute the autoconversion of r_i for r_s production: RIAUTS
 !
   ALLOCATE(ZCRIAUTI(IMICRO))
+!$acc kernels
 #ifndef MNH_BITREP
   ZCRIAUTI(:)=MIN(XCRIAUTI,10**(XACRIAUTI*(ZZT(:)-XTT)+XBCRIAUTI))
 #else
   ZCRIAUTI(:)=MIN(XCRIAUTI, BR_POW(10.,XACRIAUTI*(ZZT(:)-XTT)+XBCRIAUTI) )
 #endif
   ZZW(:) = 0.0
-  WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRIS(:)>0.0) )
+  GWORK(:) = (ZRIT(:)>XRTMIN(4)) .AND. (ZRIS(:)>0.0)
+  WHERE ( GWORK(:) )
 #ifndef MNH_BITREP
     ZZW(:) = MIN( ZRIS(:),XTIMAUTI * EXP( XTEXAUTI*(ZZT(:)-XTT) ) &
                             * MAX( ZRIT(:)-ZCRIAUTI(:),0.0 ) )
@@ -2306,6 +2374,8 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
     ZRSS(:)  = ZRSS(:)  + ZZW(:)
     ZRIS(:)  = ZRIS(:)  - ZZW(:)
   END WHERE
+!$acc end kernels
+!$acc update self(ZRIS,ZRSS)
   DEALLOCATE(ZCRIAUTI)
   IF (LBUDGET_RI) CALL BUDGET (                                                 &
                      UNPACK(ZRIS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
@@ -2317,6 +2387,7 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
 !*       3.4.6  compute the deposition on r_g: RVDEPG
 !
 !
+!$acc kernels
   WHERE ( ZRGT(:)>0.0 )
 #ifndef MNH_BITREP
     ZLBDAG(:)  = XLBG*( ZRHODREF(:)*MAX( ZRGT(:),XRTMIN(6) ) )**XLBEXG
@@ -2325,7 +2396,8 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
 #endif
   END WHERE
   ZZW(:) = 0.0
-  WHERE ( (ZRGT(:)>XRTMIN(6)) .AND. (ZRGS(:)>0.0) )
+  GWORK(:) = (ZRGT(:)>XRTMIN(6)) .AND. (ZRGS(:)>0.0)
+  WHERE ( GWORK(:) )
 #ifndef MNH_BITREP
     ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) *                               &
              ( X0DEPG*ZLBDAG(:)**XEX0DEPG + X1DEPG*ZCJ(:)*ZLBDAG(:)**XEX1DEPG )
@@ -2339,6 +2411,8 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
     ZRVS(:) = ZRVS(:) - ZZW(:)
     ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
   END WHERE
+!$acc end kernels
+!$acc update self(ZLBDAG,ZTHS,ZRVS,ZRGS)
   IF (LBUDGET_TH) CALL BUDGET (                                                 &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),   &
                                                               4,'DEPG_BU_RTH')
@@ -2348,8 +2422,6 @@ LOGICAL,DIMENSION(:),ALLOCATABLE :: GWORK
   IF (LBUDGET_RG) CALL BUDGET (                                                 &
                      UNPACK(ZRGS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
                                                              11,'DEPG_BU_RRG')
-!
-  DEALLOCATE(GWORK)
 !
   END SUBROUTINE RAIN_ICE_SLOW
 !
@@ -2369,10 +2441,12 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
 !
 !*       4.2    compute the autoconversion of r_c for r_r production: RCAUTR
 !
+!$acc kernels
     ZZW(:) = 0.0
 !
     IF ( HSUBG_AUCV == 'CLFR' ) THEN
-      WHERE( (ZRCT(:)>0.0) .AND. (ZRCS(:)>0.0) .AND. (ZCF(:)>0.0) )
+      GWORK(:) = (ZRCT(:)>0.0) .AND. (ZRCS(:)>0.0) .AND. (ZCF(:)>0.0)
+      WHERE( GWORK(:) )
        ZZW(:) = XTIMAUTC*MAX( ZRCT(:)/(ZCF(:)) -XCRIAUTC/ZRHODREF(:),0.0)
        ZZW(:) = MIN( ZRCS(:),(ZCF(:))*ZZW(:))
        ZRCS(:) = ZRCS(:) - ZZW(:)
@@ -2399,12 +2473,15 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
         ENDIF
       END DO
     ELSE
-       WHERE( (ZRCT(:)>XRTMIN(2)) .AND. (ZRCS(:)>0.0) )
+       GWORK(:) = (ZRCT(:)>XRTMIN(2)) .AND. (ZRCS(:)>0.0)
+       WHERE( GWORK(:) )
         ZZW(:) = MIN( ZRCS(:),XTIMAUTC*MAX( ZRCT(:)-XCRIAUTC/ZRHODREF(:),0.0 ) )
         ZRCS(:) = ZRCS(:) - ZZW(:)
         ZRRS(:) = ZRRS(:) + ZZW(:)
       END WHERE
     END IF
+!$acc end kernels
+!$acc update self(ZRCS,ZRRS)
 !
       IF (LBUDGET_RC) CALL BUDGET (                                               &
                        UNPACK(ZRCS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
@@ -2415,8 +2492,10 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
 !
 !*       4.3    compute the accretion of r_c for r_r production: RCACCR
 !
+!$acc kernels copyin(XEXCACCR,XRTMIN(2:3)) present(GWORK,ZZW,ZRHODREF,ZLBDAR,ZRCT,ZRRT,ZRCS,ZRRS) default(none)
     ZZW(:) = 0.0
-    WHERE( (ZRCT(:)>XRTMIN(2)) .AND. (ZRRT(:)>XRTMIN(3)) .AND. (ZRCS(:)>0.0) )
+    GWORK(:) = (ZRCT(:)>XRTMIN(2)) .AND. (ZRRT(:)>XRTMIN(3)) .AND. (ZRCS(:)>0.0)
+    WHERE( GWORK(:) )
 #ifndef MNH_BITREP
       ZZW(:) = MIN( ZRCS(:),XFCACCR * ZRCT(:)                &
                                     * ZLBDAR(:)**XEXCACCR    &
@@ -2429,6 +2508,8 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
       ZRCS(:) = ZRCS(:) - ZZW(:)
       ZRRS(:) = ZRRS(:) + ZZW(:)
     END WHERE
+!$acc end kernels
+!$acc update self(ZRCS,ZRRS)
     IF (LBUDGET_RC) CALL BUDGET (                                               &
                      UNPACK(ZRCS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
                                                               7,'ACCR_BU_RRC')
@@ -2438,8 +2519,10 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
 !
 !*       4.4    compute the evaporation of r_r: RREVAV
 !
+!$acc kernels present(ZLBDAR)
     ZZW(:) = 0.0
-    WHERE( (ZRRT(:)>XRTMIN(3)) .AND. (ZRCT(:)<=XRTMIN(2)) )
+    GWORK(:) = (ZRRT(:)>XRTMIN(3)) .AND. (ZRCT(:)<=XRTMIN(2))
+    WHERE( GWORK(:) )
 #ifndef MNH_BITREP
       ZZW(:)  = EXP( XALPW - XBETAW/ZZT(:) - XGAMW*ALOG(ZZT(:) ) ) ! es_w
 #else
@@ -2462,6 +2545,9 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
       ZRVS(:) = ZRVS(:) + ZZW(:)
       ZTHS(:) = ZTHS(:) - ZZW(:)*ZLVFACT(:)
     END WHERE
+!$acc end kernels
+!$acc update self(ZTHS,ZRVS,ZRRS)
+print *,'PW: RAIN_ICE_WARM 04'
     IF (LBUDGET_TH) CALL BUDGET (                                               &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),   &
                                                               4,'REVA_BU_RTH')
@@ -2471,9 +2557,19 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
     IF (LBUDGET_RR) CALL BUDGET (                                               &
                      UNPACK(ZRRS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
                                                               8,'REVA_BU_RRR')
-    ZW(:,:,:)=PEVAP3D(:,:,:)
-    PEVAP3D(:,:,:)=UNPACK(ZZW(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:))
+
+#ifndef _OPENACC
+    PEVAP3D(:,:,:)=UNPACK(ZZW(:),MASK=GMICRO(:,:,:),FIELD=0.0)
+#else
+!$acc kernels
+  DO JL=1,IMICRO
+    PEVAP3D(I1(JL),I2(JL),I3(JL)) = ZZW(JL)
+  END DO
+!$acc end kernels
+!$acc update self(PEVAP3D)
+#endif
 !
+print *,'PW: RAIN_ICE_WARM 05'
   END SUBROUTINE RAIN_ICE_WARM
 !
 !-------------------------------------------------------------------------------
@@ -2486,17 +2582,31 @@ REAL :: ZCRIAUTC             ! Critical cloud mixing ratio
 !
 IMPLICIT NONE
 !
+INTEGER , DIMENSION(:),ALLOCATABLE :: I1 ! Used to replace the COUNT
+INTEGER                            :: JL ! and PACK intrinsics
+!acc declare device_resident(I1)
+!$acc declare create(I1)
+!
 !-------------------------------------------------------------------------------
 !
 !*       5.1    cloud droplet riming of the aggregates
 !
+!$acc kernels
   ZZW1(:,:) = 0.0
+!$acc end kernels
 !
   ALLOCATE(GRIM(IMICRO))
+  ALLOCATE(I1(IMICRO)) !I1 is bigger than necessary but it easier to do it now (instead of computing IGRIM before allocating I1)
+!$acc kernels
 ! GRIM(:) = (ZRCT(:)>0.0) .AND. (ZRST(:)>0.0) .AND.            &
   GRIM(:) = (ZRCT(:)>XRTMIN(2)) .AND. (ZRST(:)>XRTMIN(5)) .AND.            &
                                 (ZRCS(:)>0.0) .AND. (ZZT(:)<XTT)
+!$acc end kernels
+#ifndef _OPENACC
   IGRIM = COUNT( GRIM(:) )
+#else
+  CALL COUNTJV1D_DEVICE(GRIM(:),I1(:),IGRIM)
+#endif
 !
   IF( IGRIM>0 ) THEN
 !
@@ -2509,30 +2619,52 @@ IMPLICIT NONE
 !
 !        5.1.1  select the ZLBDAS
 !
+#ifndef _OPENACC
     ZVEC1(:) = PACK( ZLBDAS(:),MASK=GRIM(:) )
+#else
+!$acc kernels
+    DO JL=1,IGRIM
+      ZVEC1(JL)  = ZLBDAS(I1(JL))
+    END DO
+!$acc end kernels
+#endif
 !
 !        5.1.2  find the next lower indice for the ZLBDAS in the geometrical
 !               set of Lbda_s used to tabulate some moments of the incomplete 
 !               gamma function
 !
-    ZVEC2(1:IGRIM) = MAX( 1.00001, MIN( FLOAT(NGAMINC)-0.00001,           &
+!$acc kernels
+    ZVEC2(:) = MAX( 1.00001, MIN( FLOAT(NGAMINC)-0.00001,           &
 #ifndef MNH_BITREP
-                          XRIMINTP1 * LOG( ZVEC1(1:IGRIM) ) + XRIMINTP2 ) )
+                          XRIMINTP1 * LOG( ZVEC1(:) ) + XRIMINTP2 ) )
 #else
-                          XRIMINTP1 * BR_LOG( ZVEC1(1:IGRIM) ) + XRIMINTP2 ) )
+                          XRIMINTP1 * BR_LOG( ZVEC1(:) ) + XRIMINTP2 ) )
 #endif
-    IVEC2(1:IGRIM) = INT( ZVEC2(1:IGRIM) )
-    ZVEC2(1:IGRIM) = ZVEC2(1:IGRIM) - FLOAT( IVEC2(1:IGRIM) )
+    IVEC2(:) = INT( ZVEC2(:) )
+    ZVEC2(:) = ZVEC2(:) - FLOAT( IVEC2(:) )
 !
 !        5.1.3  perform the linear interpolation of the normalized
 !               "2+XDS"-moment of the incomplete gamma function
 !
-    ZVEC1(1:IGRIM) =   XGAMINC_RIM1( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
-                     - XGAMINC_RIM1( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
+    ZVEC1(:) =   XGAMINC_RIM1( IVEC2(:)+1 )* ZVEC2(:)      &
+                     - XGAMINC_RIM1( IVEC2(:)   )*(ZVEC2(:) - 1.0)
+!$acc end kernels
+#ifndef _OPENACC
     ZZW(:) = UNPACK( VECTOR=ZVEC1(:),MASK=GRIM,FIELD=0.0 )
+#else
+!$acc kernels
+    ZZW(:) = 0.0
+    DO JL=1,IGRIM
+      ZZW(I1(JL)) = ZVEC1(JL)
+    END DO
+!$acc end kernels
+#endif
 !
 !        5.1.4  riming of the small sized aggregates
 !
+!$acc kernels default(none) &
+!$acc &       copyin(XEXCRIMSS,XGAMINC_RIM2) &
+!$acc &       present(GRIM,ZLBDAS,ZRHODREF,ZLSFACT,ZLVFACT,ZZW1,IVEC2,ZVEC1,ZVEC2,ZRCT,ZRCS,ZRSS,ZTHS)
     WHERE ( GRIM(:) )
 #ifndef MNH_BITREP
       ZZW1(:,1) = MIN( ZRCS(:),                                &
@@ -2553,14 +2685,27 @@ IMPLICIT NONE
 !        5.1.5  perform the linear interpolation of the normalized
 !               "XBS"-moment of the incomplete gamma function
 !
-    ZVEC1(1:IGRIM) =  XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
-                - XGAMINC_RIM2( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
+    ZVEC1(:) =  XGAMINC_RIM2( IVEC2(:)+1 )* ZVEC2(:)      &
+                - XGAMINC_RIM2( IVEC2(:)   )*(ZVEC2(:) - 1.0)
+!$acc end kernels
+#ifndef _OPENACC
     ZZW(:) = UNPACK( VECTOR=ZVEC1(:),MASK=GRIM,FIELD=0.0 )
+#else
+!$acc kernels
+    ZZW(:) = 0.0
+    DO JL=1,IGRIM
+      ZZW(I1(JL)) = ZVEC1(JL)
+    END DO
+!$acc end kernels
+#endif
 !
 !        5.1.6  riming-conversion of the large sized aggregates into graupeln
 !
 !
-    WHERE ( GRIM(:) .AND. (ZRSS(:)>0.0) )
+!$acc kernels default(none) copyin(XEXCRIMSG,XEXSRIMCG) &
+!$acc & present(GRIM,GWORK,ZLBDAS,ZRHODREF,ZLSFACT,ZLVFACT,ZZW1,ZRCT,ZRCS,ZRSS,ZRGS,ZTHS)
+    GWORK(:) = GRIM(:) .AND. (ZRSS(:)>0.0)
+    WHERE ( GWORK(:) )
 #ifndef MNH_BITREP
       ZZW1(:,2) = MIN( ZRCS(:),                     &
                XCRIMSG * ZRCT(:)                & ! RCRIMSG
@@ -2585,6 +2730,8 @@ IMPLICIT NONE
       ZRGS(:) = ZRGS(:) + ZZW1(:,2)+ZZW1(:,3)
       ZTHS(:) = ZTHS(:) + ZZW1(:,2)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RCRIMSG))
     END WHERE
+!$acc end kernels
+!$acc update self(ZTHS,ZRCS,ZRSS,ZRGS)
     DEALLOCATE(IVEC2)
     DEALLOCATE(IVEC1)
     DEALLOCATE(ZVEC2)
@@ -2606,11 +2753,17 @@ IMPLICIT NONE
 !
 !*       5.2    rain accretion onto the aggregates
 !
-  ZZW1(:,2:3) = 0.0
   ALLOCATE(GACC(IMICRO))
-   GACC(:) = (ZRRT(:)>XRTMIN(3)) .AND. (ZRST(:)>XRTMIN(5)) .AND.            &
+!$acc kernels
+  ZZW1(:,2:3) = 0.0
+  GACC(:) = (ZRRT(:)>XRTMIN(3)) .AND. (ZRST(:)>XRTMIN(5)) .AND.            &
                           (ZRRS(:)>0.0) .AND. (ZZT(:)<XTT)
+!$acc end kernels
+#ifndef _OPENACC
   IGACC = COUNT( GACC(:) )
+#else
+   CALL COUNTJV1D_DEVICE(GACC(:),I1(:),IGACC)
+#endif
 !
   IF( IGACC>0 ) THEN
 !
@@ -2624,30 +2777,40 @@ IMPLICIT NONE
 !
 !        5.2.1  select the (ZLBDAS,ZLBDAR) couplet
 !
+#ifndef _OPENACC
     ZVEC1(:) = PACK( ZLBDAS(:),MASK=GACC(:) )
     ZVEC2(:) = PACK( ZLBDAR(:),MASK=GACC(:) )
+#else
+!$acc kernels
+    DO JL=1,IGACC
+      ZVEC1(JL)  = ZLBDAS(I1(JL))
+      ZVEC2(JL)  = ZLBDAR(I1(JL))
+    END DO
+!$acc end kernels
+#endif
 !
 !        5.2.2  find the next lower indice for the ZLBDAS and for the ZLBDAR
 !               in the geometrical set of (Lbda_s,Lbda_r) couplet use to
 !               tabulate the RACCSS-kernel
 !
-    ZVEC1(1:IGACC) = MAX( 1.00001, MIN( FLOAT(NACCLBDAS)-0.00001,           &
+!$acc kernels copyin(XKER_RACCSS) present(IVEC1,IVEC2,ZVEC1,ZVEC2,ZVEC3)
+    ZVEC1(:) = MAX( 1.00001, MIN( FLOAT(NACCLBDAS)-0.00001,           &
 #ifndef MNH_BITREP
-                          XACCINTP1S * LOG( ZVEC1(1:IGACC) ) + XACCINTP2S ) )
+                          XACCINTP1S * LOG( ZVEC1(:) ) + XACCINTP2S ) )
 #else
-                          XACCINTP1S * BR_LOG( ZVEC1(1:IGACC) ) + XACCINTP2S ) )
+                          XACCINTP1S * BR_LOG( ZVEC1(:) ) + XACCINTP2S ) )
 #endif
-    IVEC1(1:IGACC) = INT( ZVEC1(1:IGACC) )
-    ZVEC1(1:IGACC) = ZVEC1(1:IGACC) - FLOAT( IVEC1(1:IGACC) )
+    IVEC1(:) = INT( ZVEC1(:) )
+    ZVEC1(:) = ZVEC1(:) - FLOAT( IVEC1(:) )
 !
-    ZVEC2(1:IGACC) = MAX( 1.00001, MIN( FLOAT(NACCLBDAR)-0.00001,           &
+    ZVEC2(:) = MAX( 1.00001, MIN( FLOAT(NACCLBDAR)-0.00001,           &
 #ifndef MNH_BITREP
-                          XACCINTP1R * LOG( ZVEC2(1:IGACC) ) + XACCINTP2R ) )
+                          XACCINTP1R * LOG( ZVEC2(:) ) + XACCINTP2R ) )
 #else
-                          XACCINTP1R * BR_LOG( ZVEC2(1:IGACC) ) + XACCINTP2R ) )
+                          XACCINTP1R * BR_LOG( ZVEC2(:) ) + XACCINTP2R ) )
 #endif
-    IVEC2(1:IGACC) = INT( ZVEC2(1:IGACC) )
-    ZVEC2(1:IGACC) = ZVEC2(1:IGACC) - FLOAT( IVEC2(1:IGACC) )
+    IVEC2(:) = INT( ZVEC2(:) )
+    ZVEC2(:) = ZVEC2(:) - FLOAT( IVEC2(:) )
 !
 !        5.2.3  perform the bilinear interpolation of the normalized
 !               RACCSS-kernel
@@ -2660,10 +2823,21 @@ IMPLICIT NONE
                     - XKER_RACCSS(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                    * (ZVEC1(JJ) - 1.0)
     END DO
+!$acc end kernels
+#ifndef _OPENACC
     ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
+#else
+!$acc kernels
+    ZZW(:) = 0.0
+    DO JL=1,IGACC
+      ZZW(I1(JL)) = ZVEC3(JL)
+    END DO
+!$acc end kernels
+#endif
 !
 !        5.2.4  raindrop accretion on the small sized aggregates
 !
+!$acc kernels present(GACC,ZLBDAS,ZRHODREF,ZZW,ZLSFACT,ZLVFACT,ZLBDAR,ZZW1,ZRRS,ZRSS,ZTHS)
     WHERE ( GACC(:) )
 #ifndef MNH_BITREP
       ZZW1(:,2) =                                            & !! coef of RRACCS
@@ -2683,10 +2857,12 @@ IMPLICIT NONE
       ZRSS(:) = ZRSS(:) + ZZW1(:,4)
       ZTHS(:) = ZTHS(:) + ZZW1(:,4)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RRACCSS))
     END WHERE
+!$acc end kernels
 !
 !        5.2.4b perform the bilinear interpolation of the normalized
 !               RACCS-kernel
 !
+!$acc kernels copyin(XKER_RACCS) present(IVEC1,IVEC2,ZVEC1,ZVEC2,ZVEC3)
     DO JJ = 1,IGACC
       ZVEC3(JJ) =  (   XKER_RACCS(IVEC2(JJ)+1,IVEC1(JJ)+1)* ZVEC1(JJ)          &
                     -  XKER_RACCS(IVEC2(JJ)+1,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
@@ -2695,11 +2871,27 @@ IMPLICIT NONE
                     -  XKER_RACCS(IVEC2(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
                                                            * (ZVEC2(JJ) - 1.0)
     END DO
+!$acc end kernels
+#ifndef _OPENACC
+!$acc update self(GACC)
     ZZW1(:,2) = ZZW1(:,2)*UNPACK( VECTOR=ZVEC3(:),MASK=GACC(:),FIELD=0.0 ) 
+#else
+!$acc update self(I1)
+!$acc kernels present(I1,ZZW1,ZVEC3)
+    ZZW1(1:I1(1)-1,2) = 0.0
+    ZZW1(I1(1),2) = ZZW1(I1(1),2)*ZVEC3(1)
+    DO JL=2,IGACC
+      ZZW1(I1(JL-1)+1:I1(JL)-1,2) = 0.0
+      ZZW1(I1(JL),2) = ZZW1(I1(JL),2)*ZVEC3(JL)
+    END DO
+    ZZW1(I1(IGACC)+1:,2) = 0.0
+!$acc end kernels
+#endif
                                                                        !! RRACCS!
 !        5.2.5  perform the bilinear interpolation of the normalized
 !               SACCRG-kernel
 !
+!$acc kernels copyin(XKER_SACCRG) present(IVEC1,IVEC2,ZVEC1,ZVEC2,ZVEC3)
     DO JJ = 1,IGACC
       ZVEC3(JJ) =  (  XKER_SACCRG(IVEC2(JJ)+1,IVEC1(JJ)+1)* ZVEC1(JJ)          &
                     - XKER_SACCRG(IVEC2(JJ)+1,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
@@ -2708,15 +2900,31 @@ IMPLICIT NONE
                     - XKER_SACCRG(IVEC2(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
                                      * (ZVEC2(JJ) - 1.0)
     END DO
+!$acc end kernels
+! #ifndef _OPENACC
+!$acc update self(GACC,ZVEC3)
     ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
+! #else
+! !acc kernels
+!     ZZW(:) = 0.0
+!     DO JL=1,IGACC
+!       ZZW(I1(JL)) = ZVEC3(JL)
+!     END DO
+! !acc end kernels
+! #endif
 !
 !        5.2.6  raindrop accretion-conversion of the large sized aggregates
 !               into graupeln
 !
+!$acc kernels
     WHERE ( GACC(:) .AND. (ZRSS(:)>0.0) )
       ZZW1(:,2) = MAX( MIN( ZRRS(:),ZZW1(:,2)-ZZW1(:,4) ),0.0 )       ! RRACCSG
     END WHERE
-    WHERE ( GACC(:) .AND. (ZRSS(:)>0.0) .AND. ZZW1(:,2)>0.0 )
+!$acc end kernels
+!$acc update self(ZZW1,ZRRS,ZRSS,ZTHS)
+!acc kernels copyin(ZZW) present(GWORK,GACC,ZLBDAS,ZRHODREF,ZLSFACT,ZLVFACT,ZLBDAR,ZZW1,ZRRS,ZRSS,ZRGS,ZTHS) default(none)
+    GWORK(:) = GACC(:) .AND. (ZRSS(:)>0.0) .AND. ZZW1(:,2)>0.0
+    WHERE ( GWORK(:) )
 #ifndef MNH_BITREP
       ZZW1(:,3) = MIN( ZRSS(:),XFSACCRG*ZZW(:)*                     & ! RSACCRG
             ( ZLBDAS(:)**(XCXS-XBS) )*( ZRHODREF(:)**(-XCEXVT-1.) ) &
@@ -2736,6 +2944,7 @@ IMPLICIT NONE
       ZTHS(:) = ZTHS(:) + ZZW1(:,2)*(ZLSFACT(:)-ZLVFACT(:)) !
                  ! f(L_f*(RRACCSG))
     END WHERE
+!acc end kernels
     DEALLOCATE(IVEC2)
     DEALLOCATE(IVEC1)
     DEALLOCATE(ZVEC3)
@@ -2758,8 +2967,11 @@ IMPLICIT NONE
 !
 !*       5.3    Conversion-Melting of the aggregates
 !
+!$acc update self(ZRVT)
+!acc kernels
   ZZW(:) = 0.0
-  WHERE( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>0.0) .AND. (ZZT(:)>XTT) )
+  GWORK(:) = (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>0.0) .AND. (ZZT(:)>XTT)
+  WHERE( GWORK(:) )
     ZZW(:) = ZRVT(:)*ZPRES(:)/((XMV/XMD)+ZRVT(:)) ! Vapor pressure
     ZZW(:) =  ZKA(:)*(XTT-ZZT(:)) +                                 &
            ( ZDV(:)*(XLVTT + ( XCPV - XCL ) * ( ZZT(:) - XTT )) &
@@ -2789,6 +3001,7 @@ IMPLICIT NONE
     ZRSS(:) = ZRSS(:) - ZZW(:)
     ZRGS(:) = ZRGS(:) + ZZW(:)
   END WHERE
+!acc end kernels
   IF (LBUDGET_RS) CALL BUDGET (                                                 &
                      UNPACK(ZRSS(:)*ZRHODJ(:),MASK=GMICRO(:,:,:),FIELD=0.0),    &
                                                              10,'CMEL_BU_RRS')
@@ -2812,6 +3025,7 @@ IMPLICIT NONE
 !
 !*       6.1    rain contact freezing
 !
+!acc kernels
   ZZW1(:,3:4) = 0.0
   WHERE( (ZRIT(:)>XRTMIN(4)) .AND. (ZRRT(:)>XRTMIN(3)) .AND.  &
                              (ZRIS(:)>0.0) .AND. (ZRRS(:)>0.0) )
@@ -2835,6 +3049,7 @@ IMPLICIT NONE
     ZRGS(:) = ZRGS(:) + ZZW1(:,3)+ZZW1(:,4)
     ZTHS(:) = ZTHS(:) + ZZW1(:,4)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*RRCFRIG)
   END WHERE
+!acc end kernels
   IF (LBUDGET_TH) CALL BUDGET (                                                 &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),   &
                                                               4,'CFRZ_BU_RTH')
@@ -3033,6 +3248,7 @@ IMPLICIT NONE
   END IF
 !
   ZRDRYG(:) = ZZW1(:,1) + ZZW1(:,2) + ZZW1(:,3) + ZZW1(:,4)
+!$acc update device(ZRDRYG)
   DEALLOCATE(GDRY)
 !
 !*       6.3    compute the Wet growth case
@@ -3075,6 +3291,7 @@ IMPLICIT NONE
                           ( ZRHODREF(:)*(XLMTT-XCL*(XTT-ZZT(:))) )   )
 #endif
   END WHERE
+!$acc update device(ZRWETG)
 !
 !*       6.4    Select Wet or Dry case
 !
@@ -3237,6 +3454,10 @@ IMPLICIT NONE
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef _OPENACC
+PRINT *,'OPENACC: RAIN_ICE_FAST_RH not yet implemented'
+CALL ABORT
+#endif
   ALLOCATE(GHAIL(IMICRO))
   GHAIL(:) = ZRHT(:)>XRTMIN(7)
   IHAIL = COUNT(GHAIL(:))
@@ -3560,6 +3781,7 @@ IMPLICIT NONE
     ZRIS(:) = 0.0
     ZCIT(:) = 0.0
   END WHERE
+!$acc update device(ZCIT)
   IF (LBUDGET_TH) CALL BUDGET (                                                 &
                  UNPACK(ZTHS(:),MASK=GMICRO(:,:,:),FIELD=PTHS)*PRHODJ(:,:,:),   &
                                                               4,'IMLT_BU_RTH')
@@ -3602,8 +3824,7 @@ IMPLICIT NONE
 !
 !-------------------------------------------------------------------------------
 !
-!
-  FUNCTION COUNTJV(LTAB,I1,I2,I3) RESULT(IC)
+  FUNCTION COUNTJV3D(LTAB,I1,I2,I3) RESULT(IC)
 !
 !*      0. DECLARATIONS
 !          ------------
@@ -3613,8 +3834,8 @@ IMPLICIT NONE
 !*       0.2  declaration of local variables
 !
 !
-LOGICAL, DIMENSION(:,:,:) :: LTAB ! Mask
-INTEGER, DIMENSION(:) :: I1,I2,I3 ! Used to replace the COUNT and PACK
+LOGICAL, DIMENSION(:,:,:), INTENT(IN)  :: LTAB     ! Mask
+INTEGER, DIMENSION(:),     INTENT(OUT) :: I1,I2,I3 ! Used to replace the COUNT and PACK
 INTEGER :: JI,JJ,JK,IC
 !
 !-------------------------------------------------------------------------------
@@ -3633,12 +3854,12 @@ DO JK = 1,SIZE(LTAB,3)
   END DO
 END DO
 !
-END FUNCTION COUNTJV
+END FUNCTION COUNTJV3D
 !
 !-------------------------------------------------------------------------------
 !
 #ifdef _OPENACC
-  SUBROUTINE COUNTJV_DEVICE(LTAB,I1,I2,I3,IC)
+  SUBROUTINE COUNTJV3D_DEVICE(LTAB,I1,I2,I3,IC)
 !
 !*      0. DECLARATIONS
 !          ------------
@@ -3673,12 +3894,12 @@ DO JK = 1,SIZE(LTAB,3)
 END DO
 !$acc end kernels
 !
-END SUBROUTINE COUNTJV_DEVICE
+END SUBROUTINE COUNTJV3D_DEVICE
 #endif
 !
 !-------------------------------------------------------------------------------
 !
-  FUNCTION COUNTJV2(LTAB,I1,I2) RESULT(IC)
+  FUNCTION COUNTJV2D(LTAB,I1,I2) RESULT(IC)
 !
 !*      0. DECLARATIONS
 !          ------------
@@ -3688,8 +3909,8 @@ IMPLICIT NONE
 !*       0.2  declaration of local variables
 !
 !
-LOGICAL, DIMENSION(:,:) :: LTAB ! Mask
-INTEGER, DIMENSION(:) :: I1,I2 ! Used to replace the COUNT and PACK
+LOGICAL, DIMENSION(:,:), INTENT(IN)  :: LTAB  ! Mask
+INTEGER, DIMENSION(:),   INTENT(OUT) :: I1,I2 ! Used to replace the COUNT and PACK
 INTEGER :: JI,JJ,IC
 !
 !-------------------------------------------------------------------------------
@@ -3705,62 +3926,108 @@ DO JJ = 1,SIZE(LTAB,2)
   END DO
 END DO
 !
-END FUNCTION COUNTJV2
+END FUNCTION COUNTJV2D
 !
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE UNPACK_DEVICE (PVEC,OMASK,PFIELD,PMAT)
+#ifdef _OPENACC
+  SUBROUTINE COUNTJV2D_DEVICE(LTAB,I1,I2,IC)
+!
+!*      0. DECLARATIONS
+!          ------------
 !
 IMPLICIT NONE
 !
-REAL,DIMENSION(:),       INTENT(IN)  :: PVEC
-LOGICAL,DIMENSION(:,:,:),INTENT(IN)  :: OMASK
-REAL,                    INTENT(IN)  :: PFIELD
-REAL,DIMENSION(:,:,:),   INTENT(OUT) :: PMAT
-!$acc declare present(PVEC,OMASK,PMAT)
+!*       0.2  declaration of local variables
 !
-INTEGER :: IDX, JI, JJ, JK
-INTEGER :: JIMAX, JJMAX, JKMAX
 !
-!PW: TODO: opti: store idx in mask array (logical->integer)
+LOGICAL, DIMENSION(:,:), INTENT(IN)  :: LTAB  ! Mask
+INTEGER, DIMENSION(:),   INTENT(OUT) :: I1,I2 ! Used to replace the COUNT and PACK
+INTEGER,                 INTENT(OUT) :: IC    ! Count
+!$acc declare present(LTAB,I1,I2)
 !
-IF ( SIZE(OMASK,1)/=SIZE(PMAT,1) .OR. SIZE(OMASK,2)/=SIZE(PMAT,2) .OR. SIZE(OMASK,3)/=SIZE(PMAT,3) ) THEN
-  PRINT *,'FATAL: UNPACK_DEVICE: arrays are not conformant'
-  CALL ABORT
-END IF
+INTEGER :: JI,JJ,IC
 !
-IDX = 1
-JIMAX = SIZE(PMAT,1)
-JJMAX = SIZE(PMAT,2)
-JKMAX = SIZE(PMAT,3)
+!-------------------------------------------------------------------------------
 !
-!$acc update self(OMASK,PVEC)
-!acc kernels
-DO JK = 1,JKMAX
-  DO JJ = 1,JJMAX
-    DO JI = 1,JIMAX
-      IF (OMASK(JI,JJ,JK)) THEN
-        PMAT(JI,JJ,JK) = PVEC(IDX)
-        IDX = IDX + 1
-      ELSE
-        PMAT(JI,JJ,JK) = PFIELD
-      END IF
-    END DO
+!$acc kernels
+IC = 0
+DO JJ = 1,SIZE(LTAB,2)
+  DO JI = 1,SIZE(LTAB,1)
+    IF( LTAB(JI,JJ) ) THEN
+      IC = IC +1
+      I1(IC) = JI
+      I2(IC) = JJ
+    END IF
   END DO
 END DO
-!acc end kernels
-!$acc update device(PMAT)
+!$acc end kernels
 !
-IF (IDX-1 > SIZE(PVEC)) THEN
-  PRINT *,'FATAL: UNPACK_DEVICE: PVEC is too small'
-  CALL ABORT
-ELSE IF (IDX-1 < SIZE(PVEC)) THEN
-  PRINT *,'WARNING: UNPACK_DEVICE: some elements of PVEC were not used (',IDX-1,'/',SIZE(PVEC),')'
-!ELSE
-!  PRINT *,'INFO: UNPACK_DEVICE: used all elements of PVEC (',IDX-1,'/',SIZE(PVEC),')'
-END IF
+END SUBROUTINE COUNTJV2D_DEVICE
+#endif
+!
+!-------------------------------------------------------------------------------
+!
+  FUNCTION COUNTJV1D(LTAB,I1) RESULT(IC)
+!
+!*      0. DECLARATIONS
+!          ------------
+!
+IMPLICIT NONE
+!
+!*       0.2  declaration of local variables
+!
+!
+LOGICAL, DIMENSION(:), INTENT(IN)  :: LTAB ! Mask
+INTEGER, DIMENSION(:), INTENT(OUT) :: I1   ! Used to replace the COUNT and PACK
+INTEGER :: JI,IC
 !
-END SUBROUTINE UNPACK_DEVICE
+!-------------------------------------------------------------------------------
+!
+IC = 0
+DO JI = 1,SIZE(LTAB,1)
+  IF( LTAB(JI) ) THEN
+    IC = IC +1
+    I1(IC) = JI
+  END IF
+END DO
+!
+END FUNCTION COUNTJV1D
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef _OPENACC
+  SUBROUTINE COUNTJV1D_DEVICE(LTAB,I1,IC)
+!
+!*      0. DECLARATIONS
+!          ------------
+!
+IMPLICIT NONE
+!
+!*       0.2  declaration of local variables
+!
+!
+LOGICAL, DIMENSION(:), INTENT(IN)  :: LTAB     ! Mask
+INTEGER, DIMENSION(:), INTENT(OUT) :: I1       ! Used to replace the COUNT and PACK
+INTEGER,               INTENT(OUT) :: IC       ! Count
+!$acc declare present(LTAB,I1)
+!
+INTEGER :: JI,IC
+!
+!-------------------------------------------------------------------------------
+!
+!$acc kernels
+IC = 0
+DO JI = 1,SIZE(LTAB,1)
+  IF( LTAB(JI) ) THEN
+    IC = IC +1
+    I1(IC) = JI
+  END IF
+END DO
+!$acc end kernels
+!
+END SUBROUTINE COUNTJV1D_DEVICE
+#endif
 !
 !-------------------------------------------------------------------------------
 !
-- 
GitLab