diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index 7bbc76e6fa4e083d001d1503435f24caf5ee0c1e..cf8b664ecc9ec8ce94a4241c71eaa80c033642cb 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -61,6 +61,7 @@ CONTAINS
 !!     R.Honnert Oct 2016 : Add ZSURF and Update with AROME
 !!     Q.Rodier  01/2019 : support RM17 mixing length
 !!     R.Honnert 01/2019 : add LGZ (reduction of the mass-flux surface closure with the resolution)
+!!     S. Riette 06/2022: compute_entr_detr is inlined
 !! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -178,7 +179,7 @@ LOGICAL                          ::  GLMIX
 LOGICAL, DIMENSION(D%NIT)              :: GWORK1
 LOGICAL, DIMENSION(D%NIT,D%NKT) :: GWORK2
 
-INTEGER  :: ITEST, JLOOP
+INTEGER  :: ITEST
 
 REAL, DIMENSION(D%NIT) :: ZRC_UP, ZRI_UP, ZRV_UP,&
                                  ZRSATW, ZRSATI,&
@@ -195,6 +196,42 @@ REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
 REAL, DIMENSION(D%NIT,16) :: ZBUF
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+!                       1.3  Declaration of additional local variables for compute_entr_detr
+!
+! Variables for cloudy part
+REAL, DIMENSION(D%NIT) :: ZKIC, ZKIC_F2  ! fraction of env. mass in the muxtures
+REAL, DIMENSION(D%NIT) :: ZEPSI,ZDELTA   ! factor entrainment detrainment
+REAL                   :: ZEPSI_CLOUD    ! factor entrainment detrainment
+REAL                   :: ZCOEFFMF_CLOUD ! factor for compputing entr. detr.
+REAL, DIMENSION(D%NIT) :: ZMIXTHL,ZMIXRT ! Thetal and rt in the mixtures
+REAL, DIMENSION(D%NIT) :: ZTHMIX         ! Theta and Thetav  of mixtures
+REAL, DIMENSION(D%NIT) :: ZRVMIX,ZRCMIX,ZRIMIX ! mixing ratios in mixtures
+REAL, DIMENSION(D%NIT) :: ZTHVMIX, ZTHVMIX_F2 ! Theta and Thetav  of mixtures
+REAL, DIMENSION(D%NIT) :: ZTHV_UP_F2     ! thv_up at flux point kk+kkl
+REAL, DIMENSION(D%NIT) :: ZRSATW_ED, ZRSATI_ED ! working arrays (mixing ratio at saturation)
+REAL, DIMENSION(D%NIT) :: ZTHV           ! theta V of environment at the bottom of cloudy part  
+REAL                   :: ZKIC_INIT      !Initial value of ZKIC
+REAL                   :: ZCOTHVU              ! Variation of Thvup between bottom and top of cloudy part
+
+! Variables for dry part
+REAL                   :: ZFOESW, ZFOESI       ! saturating vapor pressure
+REAL                   :: ZDRSATODP            ! d.Rsat/dP
+REAL                   :: ZT                   ! Temperature
+REAL                   :: ZWK0D                ! Work array
+
+! Variables for dry and cloudy parts
+REAL, DIMENSION(D%NIT) :: ZCOEFF_MINUS_HALF,&  ! Variation of Thv between mass points kk-kkl and kk
+                                  ZCOEFF_PLUS_HALF     ! Variation of Thv between mass points kk and kk+kkl
+REAL, DIMENSION(D%NIT) :: ZPRE                 ! pressure at the bottom of the cloudy part
+REAL, DIMENSION(D%NIT) :: ZG_O_THVREF_ED
+REAL, DIMENSION(D%NIT) :: ZFRAC_ICE            ! fraction of ice
+REAL, DIMENSION(D%NIT) :: ZDZ_STOP,&           ! Exact Height of the LCL above flux level KK
+                          ZTHV_MINUS_HALF,&    ! Thv at flux point(kk)  
+                          ZTHV_PLUS_HALF       ! Thv at flux point(kk+kkl)
+REAL                   :: ZDZ                  ! Delta Z used in computations
+INTEGER :: JKLIM
+
 !
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT',0,ZHOOK_HANDLE)
 
@@ -241,7 +278,9 @@ IF (OENTR_DETR) THEN
   PBUO_INTEG=0.
 
   PFRAC_ICE_UP(:,:)=0.
-  PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not used
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+  PRSAT_UP(D%NIB:D%NIE,:)=PRVM(D%NIB:D%NIE,:) ! should be initialised correctly but is (normaly) not used
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
   !cloud/dry air mixture cloud content
   ZRC_MIX = 0.
@@ -259,23 +298,28 @@ CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 
 DO JSV=1,KSV
   IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
- CALL MZM_MF(D, PSVM(:,:,JSV), ZSVM_F(:,:,JSV))
+  CALL MZM_MF(D, PSVM(:,:,JSV), ZSVM_F(:,:,JSV))
 END DO
 !                     
 !          Initialisation of updraft characteristics 
-PTHL_UP(:,:)=ZTHLM_F(:,:)
-PRT_UP(:,:)=ZRTM_F(:,:)
-PU_UP(:,:)=ZUM_F(:,:)
-PV_UP(:,:)=ZVM_F(:,:)
-PSV_UP(:,:,:)=ZSVM_F(:,:,:)
-
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+PTHL_UP(D%NIB:D%NIE,:)=ZTHLM_F(D%NIB:D%NIE,:)
+PRT_UP(D%NIB:D%NIE,:)=ZRTM_F(D%NIB:D%NIE,:)
+PU_UP(D%NIB:D%NIE,:)=ZUM_F(D%NIB:D%NIE,:)
+PV_UP(D%NIB:D%NIE,:)=ZVM_F(D%NIB:D%NIE,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT,JSV=1:KSV)
+PSV_UP(D%NIB:D%NIE,:,:)=ZSVM_F(D%NIB:D%NIE,:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT,JSV=1:KSV)
 
 ! Computation or initialisation of updraft characteristics at the KKB level
 ! thetal_up,rt_up,thetaV_up, w2,Buoyancy term and mass flux (PEMF)
-
-PTHL_UP(:,D%NKB)= ZTHLM_F(:,D%NKB)+MAX(0.,MIN(ZTMAX,(PSFTH(:)/SQRT(ZTKEM_F(:,D%NKB)))*PARAMMF%XALP_PERT))
-PRT_UP(:,D%NKB) = ZRTM_F(:,D%NKB)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,D%NKB)))*PARAMMF%XALP_PERT)) 
-
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+PTHL_UP(D%NIB:D%NIE,D%NKB)= ZTHLM_F(D%NIB:D%NIE,D%NKB)+MAX(0.,MIN(ZTMAX,(PSFTH(D%NIB:D%NIE)/SQRT(ZTKEM_F(D%NIB:D%NIE,D%NKB)))* &
+                                                                       &PARAMMF%XALP_PERT))
+PRT_UP(D%NIB:D%NIE,D%NKB) = ZRTM_F(D%NIB:D%NIE,D%NKB)+MAX(0.,MIN(ZRMAX,(PSFRV(D%NIB:D%NIE)/SQRT(ZTKEM_F(D%NIB:D%NIE,D%NKB)))* &
+                                                                      &PARAMMF%XALP_PERT)) 
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 IF (OENTR_DETR) THEN
   CALL MZM_MF(D, PTHM (:,:), ZTHM_F (:,:))
@@ -283,45 +327,58 @@ IF (OENTR_DETR) THEN
   CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F (:,:))
   CALL MZM_MF(D, PRVM(:,:), ZRVM_F (:,:))
 
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
   ! thetav at mass and flux levels
-  ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
-  ZTHVM(:,:)=PTHM(:,:)*((1.+ZRVORD*PRVM(:,:))/(1.+PRTM(:,:)))
+  ZTHVM_F(D%NIB:D%NIE,:)=ZTHM_F(D%NIB:D%NIE,:)* &
+                                    &((1.+ZRVORD*ZRVM_F(D%NIB:D%NIE,:))/(1.+ZRTM_F(D%NIB:D%NIE,:)))
+  ZTHVM(D%NIB:D%NIE,:)=PTHM(D%NIB:D%NIE,:)* &
+                                    &((1.+ZRVORD*PRVM(D%NIB:D%NIE,:))/(1.+PRTM(D%NIB:D%NIE,:)))
 
-  PTHV_UP(:,:)=ZTHVM_F(:,:)
+  PTHV_UP(D%NIB:D%NIE,:)=ZTHVM_F(D%NIB:D%NIE,:)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
   ZW_UP2(:,:)=0.
-  ZW_UP2(:,D%NKB) = MAX(0.0001,(2./3.)*ZTKEM_F(:,D%NKB))
-
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  ZW_UP2(D%NIB:D%NIE,D%NKB) = MAX(0.0001,(2./3.)*ZTKEM_F(D%NIB:D%NIE,D%NKB))
 
   ! Computation of non conservative variable for the KKB level of the updraft
   ! (all or nothing ajustement)
   PRC_UP(:,D%NKB)=0.
   PRI_UP(:,D%NKB)=0.
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
   CALL TH_R_FROM_THL_RT(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
              PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
              PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE., &
              PBUF=ZBUF(:,:))
 
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
   ! compute updraft thevav and buoyancy term at KKB level
-  PTHV_UP(:,D%NKB) = ZTH_UP(:,D%NKB)*((1+ZRVORD*PRV_UP(:,D%NKB))/(1+PRT_UP(:,D%NKB)))
+  PTHV_UP(D%NIB:D%NIE,D%NKB) = ZTH_UP(D%NIB:D%NIE,D%NKB)*((1+ZRVORD*PRV_UP(D%NIB:D%NIE,D%NKB))/(1+PRT_UP(D%NIB:D%NIE,D%NKB)))
   ! compute mean rsat in updraft
-  PRSAT_UP(:,D%NKB) = ZRSATW(:)*(1-PFRAC_ICE_UP(:,D%NKB)) + ZRSATI(:)*PFRAC_ICE_UP(:,D%NKB)
-                                                            
+  PRSAT_UP(D%NIB:D%NIE,D%NKB) = ZRSATW(D%NIB:D%NIE)*(1-PFRAC_ICE_UP(D%NIB:D%NIE,D%NKB)) + &
+                              & ZRSATI(D%NIB:D%NIE)*PFRAC_ICE_UP(D%NIB:D%NIE,D%NKB)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
   ! Closure assumption for mass flux at KKB level
   !
 
-  ZG_O_THVREF(:,:)=CST%XG/ZTHVM_F(:,:)
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+  ZG_O_THVREF(D%NIB:D%NIE,:)=CST%XG/ZTHVM_F(D%NIB:D%NIE,:)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
   ! compute L_up
   GLMIX=.TRUE.
-  ZTKEM_F(:,D%NKB)=0.
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  ZTKEM_F(D%NIB:D%NIE,D%NKB)=0.
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
   !
   IF(TURB%CTURBLEN=='RM17') THEN
     CALL GZ_M_W_MF(D, PUM, PDZZ, ZWK)
     CALL MZF_MF(D, ZWK, ZDUDZ)
     CALL GZ_M_W_MF(D, PVM, PDZZ, ZWK)
     CALL MZF_MF(D, ZWK, ZDVDZ)
-    ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
+    !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+    ZSHEAR(D%NIB:D%NIE,:) = SQRT(ZDUDZ(D%NIB:D%NIE,:)**2 + ZDVDZ(D%NIB:D%NIE,:)**2)
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
   ELSE
     ZSHEAR = 0. !no shear in bl89 mixing length
   END IF
@@ -333,30 +390,35 @@ IF (OENTR_DETR) THEN
   CALL COMPUTE_BL89_ML(D, CST, CSTURB, PDZZ,ZTKEM_F(:,D%NKB),&
                       &ZG_O_THVREF(:,D%NKB),ZTHVM,D%NKB,GLMIX,.FALSE.,ZSHEAR,ZLUP)
 #endif
-  ZLUP(:)=MAX(ZLUP(:),1.E-10)
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
+  ZLUP(D%NIB:D%NIE)=MAX(ZLUP(D%NIB:D%NIE),1.E-10)
 
   ! Compute Buoyancy flux at the ground
-  ZWTHVSURF(:) = (ZTHVM_F(:,D%NKB)/ZTHM_F(:,D%NKB))*PSFTH(:)+     &
-                (0.61*ZTHM_F(:,D%NKB))*PSFRV(:)
+  ZWTHVSURF(D%NIB:D%NIE) = (ZTHVM_F(D%NIB:D%NIE,D%NKB)/ZTHM_F(D%NIB:D%NIE,D%NKB))*PSFTH(D%NIB:D%NIE)+     &
+                (0.61*ZTHM_F(D%NIB:D%NIE,D%NKB))*PSFRV(D%NIB:D%NIE)
 
   ! Mass flux at KKB level (updraft triggered if PSFTH>0.)
   IF (PARAMMF%LGZ) THEN
-    ZSURF(:)=TANH(PARAMMF%XGZ*SQRT(PDX*PDY)/ZLUP)
+    ZSURF(D%NIB:D%NIE)=TANH(PARAMMF%XGZ*SQRT(PDX*PDY)/ZLUP(D%NIB:D%NIE))
   ELSE
-    ZSURF(:)=1.
+    ZSURF(D%NIB:D%NIE)=1.
   END IF
-  WHERE (ZWTHVSURF(:)>0.)
-    PEMF(:,D%NKB) = PARAMMF%XCMF * ZSURF(:) * ZRHO_F(:,D%NKB) *  &
-            ((ZG_O_THVREF(:,D%NKB))*ZWTHVSURF*ZLUP)**(1./3.)
-    PFRAC_UP(:,D%NKB)=MIN(PEMF(:,D%NKB)/(SQRT(ZW_UP2(:,D%NKB))*ZRHO_F(:,D%NKB)),PARAMMF%XFRAC_UP_MAX)
-    ZW_UP2(:,D%NKB)=(PEMF(:,D%NKB)/(PFRAC_UP(:,D%NKB)*ZRHO_F(:,D%NKB)))**2
-    GTEST(:)=.TRUE.
+  WHERE (ZWTHVSURF(D%NIB:D%NIE)>0.)
+    PEMF(D%NIB:D%NIE,D%NKB) = PARAMMF%XCMF * ZSURF(D%NIB:D%NIE) * ZRHO_F(D%NIB:D%NIE,D%NKB) *  &
+            ((ZG_O_THVREF(D%NIB:D%NIE,D%NKB))*ZWTHVSURF(D%NIB:D%NIE)*ZLUP(D%NIB:D%NIE))**(1./3.)
+    PFRAC_UP(D%NIB:D%NIE,D%NKB)=MIN(PEMF(D%NIB:D%NIE,D%NKB)/(SQRT(ZW_UP2(D%NIB:D%NIE,D%NKB))*ZRHO_F(D%NIB:D%NIE,D%NKB)), &
+                                   &PARAMMF%XFRAC_UP_MAX)
+    ZW_UP2(D%NIB:D%NIE,D%NKB)=(PEMF(D%NIB:D%NIE,D%NKB)/(PFRAC_UP(D%NIB:D%NIE,D%NKB)*ZRHO_F(D%NIB:D%NIE,D%NKB)))**2
+    GTEST(D%NIB:D%NIE)=.TRUE.
   ELSEWHERE
-    PEMF(:,D%NKB) =0.
-    GTEST(:)=.FALSE.
+    PEMF(D%NIB:D%NIE,D%NKB) =0.
+    GTEST(D%NIB:D%NIE)=.FALSE.
   ENDWHERE
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ELSE
-  GTEST(:)=PEMF(:,D%NKB+D%NKL)>0.
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  GTEST(D%NIB:D%NIE)=PEMF(D%NIB:D%NIE,D%NKB+D%NKL)>0.
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 END IF
 
 !--------------------------------------------------------------------------
@@ -374,26 +436,29 @@ GTESTETL(:)=.FALSE.
 
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
-! IF the updraft top is reached for all column, stop the loop on levels
-  ITEST=COUNT(GTEST)
+  ! IF the updraft top is reached for all column, stop the loop on levels
+  ITEST=COUNT(GTEST(D%NIB:D%NIE))
   IF (ITEST==0) CYCLE
 
-!       Computation of entrainment and detrainment with KF90
-!       parameterization in clouds and LR01 in subcloud layer
-
+  !       Computation of entrainment and detrainment with KF90
+  !       parameterization in clouds and LR01 in subcloud layer
 
-! to find the LCL (check if JK is LCL or not)
 
-  WHERE ((PRC_UP(:,JK)+PRI_UP(:,JK)>0.).AND.(.NOT.(GTESTLCL)))
-      KKLCL(:) = JK           
-      GTESTLCL(:)=.TRUE.
+  ! to find the LCL (check if JK is LCL or not)
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
+  WHERE ((PRC_UP(D%NIB:D%NIE,JK)+PRI_UP(D%NIB:D%NIE,JK)>0.).AND.(.NOT.(GTESTLCL(D%NIB:D%NIE))))
+      KKLCL(D%NIB:D%NIE) = JK           
+      GTESTLCL(D%NIB:D%NIE)=.TRUE.
   ENDWHERE
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 
-! COMPUTE PENTR and PDETR at mass level JK
+  ! COMPUTE PENTR and PDETR at mass level JK
   IF (OENTR_DETR) THEN
     IF(JK/=D%NKB) THEN
-      ZRC_MIX(:,JK) = ZRC_MIX(:,JK-D%NKL) ! guess of Rc of mixture
-      ZRI_MIX(:,JK) = ZRI_MIX(:,JK-D%NKL) ! guess of Ri of mixture
+      !$mnh_expand_array(JI=D%NIB:D%NIE)
+      ZRC_MIX(D%NIB:D%NIE,JK) = ZRC_MIX(D%NIB:D%NIE,JK-D%NKL) ! guess of Rc of mixture
+      ZRI_MIX(D%NIB:D%NIE,JK) = ZRI_MIX(D%NIB:D%NIE,JK-D%NKL) ! guess of Ri of mixture
+      !$mnh_end_expand_array(JI=D%NIB:D%NIE)
     ENDIF
     CALL COMPUTE_ENTR_DETR(D, CST, NEB, PARAMMF, JK,D%NKB,D%NKE,D%NKL,GTEST,GTESTLCL,HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
                            PRHODREF(:,JK),ZPRES_F(:,JK),ZPRES_F(:,JK+D%NKL),&
@@ -405,187 +470,215 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
                            PENTR(:,JK),PDETR(:,JK),ZENTR_CLD(:,JK),ZDETR_CLD(:,JK),&
                            ZBUO_INTEG_DRY(:,JK), ZBUO_INTEG_CLD(:,JK), &
                            ZPART_DRY(:)   )
-    PBUO_INTEG(:,JK)=ZBUO_INTEG_DRY(:,JK)+ZBUO_INTEG_CLD(:,JK)
+    !$mnh_expand_where(JI=D%NIB:D%NIE)
+    PBUO_INTEG(D%NIB:D%NIE,JK)=ZBUO_INTEG_DRY(D%NIB:D%NIE,JK)+ZBUO_INTEG_CLD(D%NIB:D%NIE,JK)
 
     IF (JK==D%NKB) THEN
-       PDETR(:,JK)=0.
-       ZDETR_CLD(:,JK)=0.
+       PDETR(D%NIB:D%NIE,JK)=0.
+       ZDETR_CLD(D%NIB:D%NIE,JK)=0.
     ENDIF   
  
-!       Computation of updraft characteristics at level JK+KKL
-    WHERE(GTEST)
-      ZMIX1(:)=0.5*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*(PENTR(:,JK)-PDETR(:,JK))
-      PEMF(:,JK+D%NKL)=PEMF(:,JK)*EXP(2*ZMIX1(:))
+    !       Computation of updraft characteristics at level JK+KKL
+    WHERE(GTEST(D%NIB:D%NIE))
+      ZMIX1(D%NIB:D%NIE)=0.5*(PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*(PENTR(D%NIB:D%NIE,JK)-PDETR(D%NIB:D%NIE,JK))
+      PEMF(D%NIB:D%NIE,JK+D%NKL)=PEMF(D%NIB:D%NIE,JK)*EXP(2*ZMIX1(D%NIB:D%NIE))
     ENDWHERE
-  ELSE
-    GTEST(:) = (PEMF(:,JK+D%NKL)>0.)
-  END IF 
+    !$mnh_end_expand_where(JI=D%NIB:D%NIE)
+  ELSE !OENTR_DETR
+    !$mnh_expand_array(JI=D%NIB:D%NIE)
+    GTEST(D%NIB:D%NIE) = (PEMF(D%NIB:D%NIE,JK+D%NKL)>0.)
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE)
+  END IF !OENTR_DETR
   
-
-! stop the updraft if MF becomes negative
-  WHERE (GTEST.AND.(PEMF(:,JK+D%NKL)<=0.))
-    PEMF(:,JK+D%NKL)=0.
-    KKCTL(:) = JK+D%NKL
-    GTEST(:)=.FALSE.
-    PFRAC_ICE_UP(:,JK+D%NKL)=PFRAC_ICE_UP(:,JK)
-    PRSAT_UP(:,JK+D%NKL)=PRSAT_UP(:,JK)
+  ! stop the updraft if MF becomes negative
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
+  WHERE (GTEST(D%NIB:D%NIE).AND.(PEMF(D%NIB:D%NIE,JK+D%NKL)<=0.))
+    PEMF(D%NIB:D%NIE,JK+D%NKL)=0.
+    KKCTL(D%NIB:D%NIE) = JK+D%NKL
+    GTEST(D%NIB:D%NIE)=.FALSE.
+    PFRAC_ICE_UP(D%NIB:D%NIE,JK+D%NKL)=PFRAC_ICE_UP(D%NIB:D%NIE,JK)
+    PRSAT_UP(D%NIB:D%NIE,JK+D%NKL)=PRSAT_UP(D%NIB:D%NIE,JK)
   ENDWHERE
-
-
-! If the updraft did not stop, compute cons updraft characteritics at jk+KKL
-  DO JLOOP=1,D%NIT
-    IF(GTEST(JLOOP)) THEN
-      ZMIX2(JLOOP) = (PZZ(JLOOP,JK+D%NKL)-PZZ(JLOOP,JK))*PENTR(JLOOP,JK) !&
-      ZMIX3_CLD(JLOOP) = (PZZ(JLOOP,JK+D%NKL)-PZZ(JLOOP,JK))*(1.-ZPART_DRY(JLOOP))*ZDETR_CLD(JLOOP,JK) !&                   
-      ZMIX2_CLD(JLOOP) = (PZZ(JLOOP,JK+D%NKL)-PZZ(JLOOP,JK))*(1.-ZPART_DRY(JLOOP))*ZENTR_CLD(JLOOP,JK)
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
+
+  ! If the updraft did not stop, compute cons updraft characteritics at jk+KKL
+  DO JI=D%NIB,D%NIE
+    IF(GTEST(JI)) THEN
+      ZMIX2(JI) = (PZZ(JI,JK+D%NKL)-PZZ(JI,JK))*PENTR(JI,JK) !&
+      ZMIX3_CLD(JI) = (PZZ(JI,JK+D%NKL)-PZZ(JI,JK))*(1.-ZPART_DRY(JI))*ZDETR_CLD(JI,JK) !&                   
+      ZMIX2_CLD(JI) = (PZZ(JI,JK+D%NKL)-PZZ(JI,JK))*(1.-ZPART_DRY(JI))*ZENTR_CLD(JI,JK)
 #ifdef REPRO48                  
-      PTHL_UP(JLOOP,JK+D%NKL)=(PTHL_UP(JLOOP,JK)*(1.-0.5*ZMIX2(JLOOP)) + PTHLM(JLOOP,JK)*ZMIX2(JLOOP)) &
-                            /(1.+0.5*ZMIX2(JLOOP))   
-      PRT_UP(JLOOP,JK+D%NKL) =(PRT_UP (JLOOP,JK)*(1.-0.5*ZMIX2(JLOOP)) + PRTM(JLOOP,JK)*ZMIX2(JLOOP))  &
-                            /(1.+0.5*ZMIX2(JLOOP))
+      PTHL_UP(JI,JK+D%NKL)=(PTHL_UP(JI,JK)*(1.-0.5*ZMIX2(JI)) + PTHLM(JI,JK)*ZMIX2(JI)) &
+                            /(1.+0.5*ZMIX2(JI))   
+      PRT_UP(JI,JK+D%NKL) =(PRT_UP (JI,JK)*(1.-0.5*ZMIX2(JI)) + PRTM(JI,JK)*ZMIX2(JI))  &
+                            /(1.+0.5*ZMIX2(JI))
 #else
-      PTHL_UP(JLOOP,JK+D%NKL)=PTHL_UP(JLOOP,JK)*EXP(-ZMIX2(JLOOP)) + PTHLM(JLOOP,JK)*(1-EXP(-ZMIX2(JLOOP)))
-      PRT_UP(JLOOP,JK+D%NKL) =PRT_UP (JLOOP,JK)*EXP(-ZMIX2(JLOOP)) +  PRTM(JLOOP,JK)*(1-EXP(-ZMIX2(JLOOP)))
+      PTHL_UP(JI,JK+D%NKL)=PTHL_UP(JI,JK)*EXP(-ZMIX2(JI)) + PTHLM(JI,JK)*(1-EXP(-ZMIX2(JI)))
+      PRT_UP(JI,JK+D%NKL) =PRT_UP (JI,JK)*EXP(-ZMIX2(JI)) +  PRTM(JI,JK)*(1-EXP(-ZMIX2(JI)))
 #endif
     ENDIF
   ENDDO
   
-
   IF(OMIXUV) THEN
     IF(JK/=D%NKB) THEN
-      WHERE(GTEST)
-        PU_UP(:,JK+D%NKL) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
-                          0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
-                          ((PUM(:,JK+D%NKL)-PUM(:,JK))/PDZZ(:,JK+D%NKL)+&
-                           (PUM(:,JK)-PUM(:,JK-D%NKL))/PDZZ(:,JK))        )   &
-                          /(1+0.5*ZMIX2(:))
-        PV_UP(:,JK+D%NKL) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
-                          0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
-                          ((PVM(:,JK+D%NKL)-PVM(:,JK))/PDZZ(:,JK+D%NKL)+&
-                           (PVM(:,JK)-PVM(:,JK-D%NKL))/PDZZ(:,JK))    )   &
-                          /(1+0.5*ZMIX2(:))
+      !$mnh_expand_where(JI=D%NIB:D%NIE)
+      WHERE(GTEST(D%NIB:D%NIE))
+        PU_UP(D%NIB:D%NIE,JK+D%NKL) = (PU_UP(D%NIB:D%NIE,JK)*(1-0.5*ZMIX2(D%NIB:D%NIE)) + PUM(D%NIB:D%NIE,JK)*ZMIX2(D%NIB:D%NIE)+ &
+                          0.5*PARAMMF%XPRES_UV*(PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*&
+                          ((PUM(D%NIB:D%NIE,JK+D%NKL)-PUM(D%NIB:D%NIE,JK))/PDZZ(D%NIB:D%NIE,JK+D%NKL)+&
+                           (PUM(D%NIB:D%NIE,JK)-PUM(D%NIB:D%NIE,JK-D%NKL))/PDZZ(D%NIB:D%NIE,JK))        )   &
+                          /(1+0.5*ZMIX2(D%NIB:D%NIE))
+        PV_UP(D%NIB:D%NIE,JK+D%NKL) = (PV_UP(D%NIB:D%NIE,JK)*(1-0.5*ZMIX2(D%NIB:D%NIE)) + PVM(D%NIB:D%NIE,JK)*ZMIX2(D%NIB:D%NIE)+ &
+                          0.5*PARAMMF%XPRES_UV*(PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*&
+                          ((PVM(D%NIB:D%NIE,JK+D%NKL)-PVM(D%NIB:D%NIE,JK))/PDZZ(D%NIB:D%NIE,JK+D%NKL)+&
+                           (PVM(D%NIB:D%NIE,JK)-PVM(D%NIB:D%NIE,JK-D%NKL))/PDZZ(D%NIB:D%NIE,JK))    )   &
+                          /(1+0.5*ZMIX2(D%NIB:D%NIE))
       ENDWHERE
+      !$mnh_end_expand_where(JI=D%NIB:D%NIE)
     ELSE
-      WHERE(GTEST)
-        PU_UP(:,JK+D%NKL) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
-                          0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
-                          ((PUM(:,JK+D%NKL)-PUM(:,JK))/PDZZ(:,JK+D%NKL))        )   &
-                          /(1+0.5*ZMIX2(:))
-        PV_UP(:,JK+D%NKL) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
-                          0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
-                          ((PVM(:,JK+D%NKL)-PVM(:,JK))/PDZZ(:,JK+D%NKL))    )   &
-                          /(1+0.5*ZMIX2(:))
+      !$mnh_expand_where(JI=D%NIB:D%NIE)
+      WHERE(GTEST(D%NIB:D%NIE))
+        PU_UP(D%NIB:D%NIE,JK+D%NKL) = (PU_UP(D%NIB:D%NIE,JK)*(1-0.5*ZMIX2(D%NIB:D%NIE)) + PUM(D%NIB:D%NIE,JK)*ZMIX2(D%NIB:D%NIE)+ &
+                          0.5*PARAMMF%XPRES_UV*(PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*&
+                          ((PUM(D%NIB:D%NIE,JK+D%NKL)-PUM(D%NIB:D%NIE,JK))/PDZZ(D%NIB:D%NIE,JK+D%NKL))        )   &
+                          /(1+0.5*ZMIX2(D%NIB:D%NIE))
+        PV_UP(D%NIB:D%NIE,JK+D%NKL) = (PV_UP(D%NIB:D%NIE,JK)*(1-0.5*ZMIX2(D%NIB:D%NIE)) + PVM(D%NIB:D%NIE,JK)*ZMIX2(D%NIB:D%NIE)+ &
+                          0.5*PARAMMF%XPRES_UV*(PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*&
+                          ((PVM(D%NIB:D%NIE,JK+D%NKL)-PVM(D%NIB:D%NIE,JK))/PDZZ(D%NIB:D%NIE,JK+D%NKL))    )   &
+                          /(1+0.5*ZMIX2(D%NIB:D%NIE))
       ENDWHERE
-
+      !$mnh_end_expand_where(JI=D%NIB:D%NIE)
     ENDIF
-  ENDIF
+  ENDIF !OMIXUV
   DO JSV=1,KSV 
-     IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
-      WHERE(GTEST) 
-           PSV_UP(:,JK+D%NKL,JSV) = (PSV_UP (:,JK,JSV)*(1-0.5*ZMIX2(:)) + &
-                        PSVM(:,JK,JSV)*ZMIX2(:))  /(1+0.5*ZMIX2(:))
-      ENDWHERE                        
+    IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
+    !$mnh_expand_where(JI=D%NIB:D%NIE)
+    WHERE(GTEST(D%NIB:D%NIE)) 
+      PSV_UP(D%NIB:D%NIE,JK+D%NKL,JSV) = (PSV_UP(D%NIB:D%NIE,JK,JSV)*(1-0.5*ZMIX2(D%NIB:D%NIE)) + &
+                   PSVM(D%NIB:D%NIE,JK,JSV)*ZMIX2(D%NIB:D%NIE))  /(1+0.5*ZMIX2(D%NIB:D%NIE))
+    ENDWHERE
+    !$mnh_end_expand_where(JI=D%NIB:D%NIE)
   END DO  
   
- IF (OENTR_DETR) THEN
-
-! Compute non cons. var. at level JK+KKL
-  ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
-  ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
-  CALL TH_R_FROM_THL_RT(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
-          PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL),              &
-          ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE., &
-          PBUF=ZBUF(:,:))
-  WHERE(GTEST)
-    PRC_UP(:,JK+D%NKL)=ZRC_UP(:)
-    PRV_UP(:,JK+D%NKL)=ZRV_UP(:)
-    PRI_UP(:,JK+D%NKL)=ZRI_UP(:)
-    PRSAT_UP(:,JK+D%NKL) = ZRSATW(:)*(1-PFRAC_ICE_UP(:,JK+D%NKL)) + ZRSATI(:)*PFRAC_ICE_UP(:,JK+D%NKL)
-  ENDWHERE
-  
+  IF (OENTR_DETR) THEN
 
-! Compute the updraft theta_v, buoyancy and w**2 for level JK+KKL
-  WHERE(GTEST)
-    PTHV_UP(:,JK+D%NKL) = ZTH_UP(:,JK+D%NKL)*((1+ZRVORD*PRV_UP(:,JK+D%NKL))/(1+PRT_UP(:,JK+D%NKL)))
-    WHERE (ZBUO_INTEG_DRY(:,JK)>0.)
-      ZW_UP2(:,JK+D%NKL)  = ZW_UP2(:,JK) + 2.*(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)* ZBUO_INTEG_DRY(:,JK)
-    ELSEWHERE
-      ZW_UP2(:,JK+D%NKL)  = ZW_UP2(:,JK) + 2.*PARAMMF%XABUO* ZBUO_INTEG_DRY(:,JK)
+    ! Compute non cons. var. at level JK+KKL
+    !$mnh_expand_array(JI=D%NIB:D%NIE)
+    ZRC_UP(D%NIB:D%NIE)=PRC_UP(D%NIB:D%NIE,JK) ! guess = level just below
+    ZRI_UP(D%NIB:D%NIE)=PRI_UP(D%NIB:D%NIE,JK) ! guess = level just below
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE)
+    CALL TH_R_FROM_THL_RT(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
+            PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL),              &
+            ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE., &
+            PBUF=ZBUF(:,:))
+    !$mnh_expand_where(JI=D%NIB:D%NIE)
+    WHERE(GTEST(D%NIB:D%NIE))
+      PRC_UP(D%NIB:D%NIE,JK+D%NKL)=ZRC_UP(D%NIB:D%NIE)
+      PRV_UP(D%NIB:D%NIE,JK+D%NKL)=ZRV_UP(D%NIB:D%NIE)
+      PRI_UP(D%NIB:D%NIE,JK+D%NKL)=ZRI_UP(D%NIB:D%NIE)
+      PRSAT_UP(D%NIB:D%NIE,JK+D%NKL) = ZRSATW(D%NIB:D%NIE)*(1-PFRAC_ICE_UP(D%NIB:D%NIE,JK+D%NKL)) + &
+                                     & ZRSATI(D%NIB:D%NIE)*PFRAC_ICE_UP(D%NIB:D%NIE,JK+D%NKL)
     ENDWHERE
-    ZW_UP2(:,JK+D%NKL)  = ZW_UP2(:,JK+D%NKL)*(1.-(PARAMMF%XBDETR*ZMIX3_CLD(:)+PARAMMF%XBENTR*ZMIX2_CLD(:)))&
-            /(1.+(PARAMMF%XBDETR*ZMIX3_CLD(:)+PARAMMF%XBENTR*ZMIX2_CLD(:))) &
-            +2.*(PARAMMF%XABUO)*ZBUO_INTEG_CLD(:,JK)/(1.+(PARAMMF%XBDETR*ZMIX3_CLD(:)+PARAMMF%XBENTR*ZMIX2_CLD(:)))
- ENDWHERE
 
+    ! Compute the updraft theta_v, buoyancy and w**2 for level JK+KKL
+    WHERE(GTEST(D%NIB:D%NIE))
+      PTHV_UP(D%NIB:D%NIE,JK+D%NKL) = ZTH_UP(D%NIB:D%NIE,JK+D%NKL)* &
+                                    & ((1+ZRVORD*PRV_UP(D%NIB:D%NIE,JK+D%NKL))/(1+PRT_UP(D%NIB:D%NIE,JK+D%NKL)))
+      WHERE (ZBUO_INTEG_DRY(D%NIB:D%NIE,JK)>0.)
+        ZW_UP2(D%NIB:D%NIE,JK+D%NKL)  = ZW_UP2(D%NIB:D%NIE,JK) + 2.*(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)* &
+                                                                &ZBUO_INTEG_DRY(D%NIB:D%NIE,JK)
+      ELSEWHERE
+        ZW_UP2(D%NIB:D%NIE,JK+D%NKL)  = ZW_UP2(D%NIB:D%NIE,JK) + 2.*PARAMMF%XABUO* ZBUO_INTEG_DRY(D%NIB:D%NIE,JK)
+      ENDWHERE
+      ZW_UP2(D%NIB:D%NIE,JK+D%NKL)  = ZW_UP2(D%NIB:D%NIE,JK+D%NKL)*(1.-(PARAMMF%XBDETR*ZMIX3_CLD(D%NIB:D%NIE)+ &
+                                                                       &PARAMMF%XBENTR*ZMIX2_CLD(D%NIB:D%NIE)))&
+              /(1.+(PARAMMF%XBDETR*ZMIX3_CLD(D%NIB:D%NIE)+PARAMMF%XBENTR*ZMIX2_CLD(D%NIB:D%NIE))) &
+              +2.*(PARAMMF%XABUO)*ZBUO_INTEG_CLD(D%NIB:D%NIE,JK)/ &
+              &(1.+(PARAMMF%XBDETR*ZMIX3_CLD(D%NIB:D%NIE)+PARAMMF%XBENTR*ZMIX2_CLD(D%NIB:D%NIE)))
+    ENDWHERE
 
-  ! Test if the updraft has reach the ETL
-  GTESTETL(:)=.FALSE.
-  WHERE (GTEST.AND.(PBUO_INTEG(:,JK)<=0.))
-      KKETL(:) = JK+D%NKL
-      GTESTETL(:)=.TRUE.
-  ENDWHERE
+    ! Test if the updraft has reach the ETL
+    WHERE (GTEST(D%NIB:D%NIE).AND.(PBUO_INTEG(D%NIB:D%NIE,JK)<=0.))
+      KKETL(D%NIB:D%NIE) = JK+D%NKL
+      GTESTETL(D%NIB:D%NIE)=.TRUE.
+    ELSEWHERE
+      GTESTETL(D%NIB:D%NIE)=.FALSE.
+    ENDWHERE
 
-  ! Test is we have reached the top of the updraft
-  WHERE (GTEST.AND.((ZW_UP2(:,JK+D%NKL)<=0.).OR.(PEMF(:,JK+D%NKL)<=0.)))
-      ZW_UP2(:,JK+D%NKL)=0.
-      PEMF(:,JK+D%NKL)=0.
-      GTEST(:)=.FALSE.
-      PTHL_UP(:,JK+D%NKL)=ZTHLM_F(:,JK+D%NKL)
-      PRT_UP(:,JK+D%NKL)=ZRTM_F(:,JK+D%NKL)
-      PRC_UP(:,JK+D%NKL)=0.
-      PRI_UP(:,JK+D%NKL)=0.
-      PRV_UP(:,JK+D%NKL)=0.
-      PTHV_UP(:,JK+D%NKL)=ZTHVM_F(:,JK+D%NKL)
-      PFRAC_UP(:,JK+D%NKL)=0.
-      KKCTL(:)=JK+D%NKL
-  ENDWHERE
+    ! Test is we have reached the top of the updraft
+    WHERE (GTEST(D%NIB:D%NIE).AND.((ZW_UP2(D%NIB:D%NIE,JK+D%NKL)<=0.).OR.(PEMF(D%NIB:D%NIE,JK+D%NKL)<=0.)))
+        ZW_UP2(D%NIB:D%NIE,JK+D%NKL)=0.
+        PEMF(D%NIB:D%NIE,JK+D%NKL)=0.
+        GTEST(D%NIB:D%NIE)=.FALSE.
+        PTHL_UP(D%NIB:D%NIE,JK+D%NKL)=ZTHLM_F(D%NIB:D%NIE,JK+D%NKL)
+        PRT_UP(D%NIB:D%NIE,JK+D%NKL)=ZRTM_F(D%NIB:D%NIE,JK+D%NKL)
+        PRC_UP(D%NIB:D%NIE,JK+D%NKL)=0.
+        PRI_UP(D%NIB:D%NIE,JK+D%NKL)=0.
+        PRV_UP(D%NIB:D%NIE,JK+D%NKL)=0.
+        PTHV_UP(D%NIB:D%NIE,JK+D%NKL)=ZTHVM_F(D%NIB:D%NIE,JK+D%NKL)
+        PFRAC_UP(D%NIB:D%NIE,JK+D%NKL)=0.
+        KKCTL(D%NIB:D%NIE)=JK+D%NKL
+    ENDWHERE
  
-  ! compute frac_up at JK+KKL
-  WHERE (GTEST)
-    PFRAC_UP(:,JK+D%NKL)=PEMF(:,JK+D%NKL)/(SQRT(ZW_UP2(:,JK+D%NKL))*ZRHO_F(:,JK+D%NKL))
-  ENDWHERE
-
-  ! Updraft fraction must be smaller than XFRAC_UP_MAX
-  WHERE (GTEST)
-    PFRAC_UP(:,JK+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX,PFRAC_UP(:,JK+D%NKL))
-  ENDWHERE
+    ! compute frac_up at JK+KKL
+    WHERE (GTEST(D%NIB:D%NIE))
+      PFRAC_UP(D%NIB:D%NIE,JK+D%NKL)=PEMF(D%NIB:D%NIE,JK+D%NKL)/(SQRT(ZW_UP2(D%NIB:D%NIE,JK+D%NKL))*ZRHO_F(D%NIB:D%NIE,JK+D%NKL))
+    ENDWHERE
 
-  ! When cloudy and non-buoyant, updraft fraction must decrease
-  WHERE ((GTEST.AND.GTESTETL).AND.GTESTLCL)
-    PFRAC_UP(:,JK+D%NKL)=MIN(PFRAC_UP(:,JK+D%NKL),PFRAC_UP(:,JK))
-  ENDWHERE
+    ! Updraft fraction must be smaller than XFRAC_UP_MAX
+    WHERE (GTEST(D%NIB:D%NIE))
+      PFRAC_UP(D%NIB:D%NIE,JK+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX,PFRAC_UP(D%NIB:D%NIE,JK+D%NKL))
+    ENDWHERE
 
-  ! Mass flux is updated with the new updraft fraction
-  IF (OENTR_DETR) PEMF(:,JK+D%NKL)=PFRAC_UP(:,JK+D%NKL)*SQRT(ZW_UP2(:,JK+D%NKL))*ZRHO_F(:,JK+D%NKL)
+    ! When cloudy and non-buoyant, updraft fraction must decrease
+    WHERE ((GTEST(D%NIB:D%NIE).AND.GTESTETL(D%NIB:D%NIE)).AND.GTESTLCL(D%NIB:D%NIE))
+      PFRAC_UP(D%NIB:D%NIE,JK+D%NKL)=MIN(PFRAC_UP(D%NIB:D%NIE,JK+D%NKL),PFRAC_UP(D%NIB:D%NIE,JK))
+    ENDWHERE
 
- END IF
+    ! Mass flux is updated with the new updraft fraction
+    IF (OENTR_DETR) PEMF(D%NIB:D%NIE,JK+D%NKL)=PFRAC_UP(D%NIB:D%NIE,JK+D%NKL)*SQRT(ZW_UP2(D%NIB:D%NIE,JK+D%NKL))* &
+                                              &ZRHO_F(D%NIB:D%NIE,JK+D%NKL)
+    !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 
+  END IF !OENTR_DETR
 ENDDO
 
 IF(OENTR_DETR) THEN
 
-  PW_UP(:,:)=SQRT(ZW_UP2(:,:))
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+  PW_UP(D%NIB:D%NIE,:)=SQRT(ZW_UP2(D%NIB:D%NIE,:))
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
-  PEMF(:,D%NKB) =0.
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  PEMF(D%NIB:D%NIE,D%NKB) =0.
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
-! Limits the shallow convection scheme when cloud heigth is higher than 3000m.
-! To do this, mass flux is multiplied by a coefficient decreasing linearly
-! from 1 (for clouds of ZDEPTH_MAX1 m of depth) to 0 (for clouds of ZDEPTH_MAX2 m of depth).
-! This way, all MF fluxes are diminished by this amount.
-! Diagnosed cloud fraction is also multiplied by the same coefficient.
-!
-  DO JI=1,D%NIT
+  ! Limits the shallow convection scheme when cloud heigth is higher than 3000m.
+  ! To do this, mass flux is multiplied by a coefficient decreasing linearly
+  ! from 1 (for clouds of ZDEPTH_MAX1 m of depth) to 0 (for clouds of ZDEPTH_MAX2 m of depth).
+  ! This way, all MF fluxes are diminished by this amount.
+  ! Diagnosed cloud fraction is also multiplied by the same coefficient.
+  !
+  DO JI=D%NIB,D%NIE
      PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
   END DO
 
-  GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
-  GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=MAX(D%NKU,D%NKA) )
-  ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=D%NKT)
-  ZCOEF=MIN(MAX(ZCOEF,0.),1.)
-  WHERE (GWORK2) 
-    PEMF(:,:)     = PEMF(:,:)     * ZCOEF(:,:)
-    PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  GWORK1(D%NIB:D%NIE)= (GTESTLCL(D%NIB:D%NIE) .AND. (PDEPTH(D%NIB:D%NIE) > ZDEPTH_MAX1) )
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
+  DO JK=1, D%NKT
+    !$mnh_expand_array(JI=D%NIB:D%NIE)
+    GWORK2(D%NIB:D%NIE,JK) = GWORK1(D%NIB:D%NIE)
+    ZCOEF(D%NIB:D%NIE,JK) = (1.-(PDEPTH(D%NIB:D%NIE)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1))
+    ZCOEF(D%NIB:D%NIE,JK)=MIN(MAX(ZCOEF(D%NIB:D%NIE,JK),0.),1.)
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE)
+  ENDDO
+  !$mnh_expand_where(JI=D%NIB:D%NIE,JK=1:D%NKT)
+  WHERE (GWORK2(D%NIB:D%NIE,:)) 
+    PEMF(D%NIB:D%NIE,:)     = PEMF(D%NIB:D%NIE,:)     * ZCOEF(D%NIB:D%NIE,:)
+    PFRAC_UP(D%NIB:D%NIE,:) = PFRAC_UP(D%NIB:D%NIE,:) * ZCOEF(D%NIB:D%NIE,:)
   ENDWHERE
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE,JK=1:D%NKT)
 ENDIF
 
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT',1,ZHOOK_HANDLE)
@@ -662,9 +755,6 @@ USE MODD_CST,             ONLY: CST_t
 USE MODD_NEB,             ONLY: NEB_t
 USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
 !
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-
 IMPLICIT NONE
 !
 !                         
@@ -718,143 +808,114 @@ REAL, DIMENSION(D%NIT),   INTENT(OUT)    ::  PPART_DRY ! ratio of dry part at th
 !
 !                       1.2  Declaration of local variables
 !
+!     Local array declaration must be put in the compute_updraft subroutine
+!     For simplicity all local variables (including scalars) are moved in the compute_updraft subroutine
 !
 
-! Variables for cloudy part
-REAL, DIMENSION(D%NIT) :: ZKIC, ZKIC_F2  ! fraction of env. mass in the muxtures
-REAL, DIMENSION(D%NIT) :: ZEPSI,ZDELTA   ! factor entrainment detrainment
-REAL                   :: ZEPSI_CLOUD    ! factor entrainment detrainment
-REAL                   :: ZCOEFFMF_CLOUD ! factor for compputing entr. detr.
-REAL, DIMENSION(D%NIT) :: ZMIXTHL,ZMIXRT ! Thetal and rt in the mixtures
-REAL, DIMENSION(D%NIT) :: ZTHMIX         ! Theta and Thetav  of mixtures
-REAL, DIMENSION(D%NIT) :: ZRVMIX,ZRCMIX,ZRIMIX ! mixing ratios in mixtures
-REAL, DIMENSION(D%NIT) :: ZTHVMIX, ZTHVMIX_F2 ! Theta and Thetav  of mixtures
-REAL, DIMENSION(D%NIT) :: ZTHV_UP_F2     ! thv_up at flux point kk+kkl
-REAL, DIMENSION(D%NIT) :: ZRSATW, ZRSATI ! working arrays (mixing ratio at saturation)
-REAL, DIMENSION(D%NIT) :: ZTHV           ! theta V of environment at the bottom of cloudy part  
-REAL                   :: ZKIC_INIT      !Initial value of ZKIC
-REAL                   :: ZCOTHVU              ! Variation of Thvup between bottom and top of cloudy part
-
-! Variables for dry part
-REAL                   :: ZFOESW, ZFOESI       ! saturating vapor pressure
-REAL                   :: ZDRSATODP            ! d.Rsat/dP
-REAL                   :: ZT                   ! Temperature
-REAL                   :: ZWK                  ! Work array
-
-! Variables for dry and cloudy parts
-REAL, DIMENSION(D%NIT) :: ZCOEFF_MINUS_HALF,&  ! Variation of Thv between mass points kk-kkl and kk
-                                  ZCOEFF_PLUS_HALF     ! Variation of Thv between mass points kk and kk+kkl
-REAL, DIMENSION(D%NIT) :: ZPRE                 ! pressure at the bottom of the cloudy part
-REAL, DIMENSION(D%NIT) :: ZG_O_THVREF
-REAL, DIMENSION(D%NIT) :: ZFRAC_ICE            ! fraction of ice
-REAL                   :: ZRVORD               ! RV/RD
-REAL, DIMENSION(D%NIT) :: ZDZ_STOP,&           ! Exact Height of the LCL above flux level KK
-                          ZTHV_MINUS_HALF,&    ! Thv at flux point(kk)  
-                          ZTHV_PLUS_HALF       ! Thv at flux point(kk+kkl)
-REAL                   :: ZDZ                  ! Delta Z used in computations
-INTEGER :: JI, JLOOP
-REAL, DIMENSION(D%NIT, 16) :: ZBUF
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !----------------------------------------------------------------------------------
                         
 !                1.3 Initialisation
 !                ------------------
-  IF (LHOOK) CALL DR_HOOK('COMPUTE_ENTR_DETR',0,ZHOOK_HANDLE)
-  
-  ZRVORD   = CST%XRV / CST%XRD   !=1.607
-  ZG_O_THVREF(:)=CST%XG/PTHVM(:,KK)
-  ZCOEFFMF_CLOUD=PARAMMF%XENTR_MF * CST%XG / PARAMMF%XCRAD_MF
   
-  ZFRAC_ICE(:)=PFRAC_ICE(:) ! to not modify fraction of ice
- 
-  ZPRE(:)=PPRE_MINUS_HALF(:)
+ZCOEFFMF_CLOUD=PARAMMF%XENTR_MF * CST%XG / PARAMMF%XCRAD_MF
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+ZG_O_THVREF_ED(D%NIB:D%NIE)=CST%XG/PTHVM(D%NIB:D%NIE,KK)
+
+ZFRAC_ICE(D%NIB:D%NIE)=PFRAC_ICE(D%NIB:D%NIE) ! to not modify fraction of ice
+
+ZPRE(D%NIB:D%NIE)=PPRE_MINUS_HALF(D%NIB:D%NIE)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 !                1.4 Estimation of PPART_DRY
-  DO JLOOP=1,SIZE(OTEST)
-    IF(OTEST(JLOOP) .AND. OTESTLCL(JLOOP)) THEN
-      !No dry part when condensation level is reached
-      PPART_DRY(JLOOP)=0.
-      ZDZ_STOP(JLOOP)=0.
-      ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)
-    ELSE IF (OTEST(JLOOP) .AND. .NOT. OTESTLCL(JLOOP)) THEN
-      !Temperature at flux level KK
-      ZT=PTH_UP(JLOOP)*(PPRE_MINUS_HALF(JLOOP)/CST%XP00) ** (CST%XRD/CST%XCPD)
-      !Saturating vapor pressure at flux level KK
-      ZFOESW = MIN(EXP( CST%XALPW - CST%XBETAW/ZT - CST%XGAMW*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JLOOP))
-      ZFOESI = MIN(EXP( CST%XALPI - CST%XBETAI/ZT - CST%XGAMI*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JLOOP))
-      !Computation of d.Rsat / dP (partial derivations with respect to P and T
-      !and use of T=Theta*(P/P0)**(R/Cp) to transform dT into dP with theta_up
-      !constant at the vertical)
-      ZDRSATODP=(CST%XBETAW/ZT-CST%XGAMW)*(1-ZFRAC_ICE(JLOOP))+(CST%XBETAI/ZT-CST%XGAMI)*ZFRAC_ICE(JLOOP)
-      ZDRSATODP=((CST%XRD/CST%XCPD)*ZDRSATODP-1.)*PRSAT_UP(JLOOP)/ &
-                  &(PPRE_MINUS_HALF(JLOOP)-(ZFOESW*(1-ZFRAC_ICE(JLOOP)) + ZFOESI*ZFRAC_ICE(JLOOP)))
-      !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
-      !where Rsat is equal to PRT_UP
-      ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)+(PRT_UP(JLOOP)-PRSAT_UP(JLOOP))/ZDRSATODP
-      !Fraction of dry part (computed with pressure and used with heights, no
-      !impact found when using log function here and for pressure on flux levels
-      !computation)
-      PPART_DRY(JLOOP)=MAX(0., MIN(1., (PPRE_MINUS_HALF(JLOOP)-ZPRE(JLOOP))/(PPRE_MINUS_HALF(JLOOP)-PPRE_PLUS_HALF(JLOOP))))
-      !Height above flux level KK of the cloudy part
-      ZDZ_STOP(JLOOP) = (PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*PPART_DRY(JLOOP)
-    ELSE
-      PPART_DRY(JLOOP)=0. ! value does not matter, here
-    END IF
-  END DO
+DO JI=D%NIB,D%NIE
+  IF(OTEST(JI) .AND. OTESTLCL(JI)) THEN
+    !No dry part when condensation level is reached
+    PPART_DRY(JI)=0.
+    ZDZ_STOP(JI)=0.
+    ZPRE(JI)=PPRE_MINUS_HALF(JI)
+  ELSE IF (OTEST(JI) .AND. .NOT. OTESTLCL(JI)) THEN
+    !Temperature at flux level KK
+    ZT=PTH_UP(JI)*(PPRE_MINUS_HALF(JI)/CST%XP00) ** (CST%XRD/CST%XCPD)
+    !Saturating vapor pressure at flux level KK
+    ZFOESW = MIN(EXP( CST%XALPW - CST%XBETAW/ZT - CST%XGAMW*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JI))
+    ZFOESI = MIN(EXP( CST%XALPI - CST%XBETAI/ZT - CST%XGAMI*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JI))
+    !Computation of d.Rsat / dP (partial derivations with respect to P and T
+    !and use of T=Theta*(P/P0)**(R/Cp) to transform dT into dP with theta_up
+    !constant at the vertical)
+    ZDRSATODP=(CST%XBETAW/ZT-CST%XGAMW)*(1-ZFRAC_ICE(JI))+(CST%XBETAI/ZT-CST%XGAMI)*ZFRAC_ICE(JI)
+    ZDRSATODP=((CST%XRD/CST%XCPD)*ZDRSATODP-1.)*PRSAT_UP(JI)/ &
+                &(PPRE_MINUS_HALF(JI)-(ZFOESW*(1-ZFRAC_ICE(JI)) + ZFOESI*ZFRAC_ICE(JI)))
+    !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
+    !where Rsat is equal to PRT_UP
+    ZPRE(JI)=PPRE_MINUS_HALF(JI)+(PRT_UP(JI)-PRSAT_UP(JI))/ZDRSATODP
+    !Fraction of dry part (computed with pressure and used with heights, no
+    !impact found when using log function here and for pressure on flux levels
+    !computation)
+    PPART_DRY(JI)=MAX(0., MIN(1., (PPRE_MINUS_HALF(JI)-ZPRE(JI))/(PPRE_MINUS_HALF(JI)-PPRE_PLUS_HALF(JI))))
+    !Height above flux level KK of the cloudy part
+    ZDZ_STOP(JI) = (PZZ(JI,KK+KKL)-PZZ(JI,KK))*PPART_DRY(JI)
+  ELSE
+    PPART_DRY(JI)=0. ! value does not matter, here
+  END IF
+END DO
 
 !               1.5 Gradient and flux values of thetav
-  IF(KK/=KKB)THEN
-    ZCOEFF_MINUS_HALF(:)=((PTHVM(:,KK)-PTHVM(:,KK-KKL))/PDZZ(:,KK))
-    ZTHV_MINUS_HALF(:) = PTHVM(:,KK) - ZCOEFF_MINUS_HALF(:)*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))
-  ELSE
-    ZCOEFF_MINUS_HALF(:)=0.
-    ZTHV_MINUS_HALF(:) = PTHVM(:,KK)
-  ENDIF
-  ZCOEFF_PLUS_HALF(:)  = ((PTHVM(:,KK+KKL)-PTHVM(:,KK))/PDZZ(:,KK+KKL))
-  ZTHV_PLUS_HALF(:)  = PTHVM(:,KK) + ZCOEFF_PLUS_HALF(:)*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+IF(KK/=KKB)THEN
+  ZCOEFF_MINUS_HALF(D%NIB:D%NIE)=((PTHVM(D%NIB:D%NIE,KK)-PTHVM(D%NIB:D%NIE,KK-KKL))/PDZZ(D%NIB:D%NIE,KK))
+  ZTHV_MINUS_HALF(D%NIB:D%NIE) = PTHVM(D%NIB:D%NIE,KK) - &
+                               & ZCOEFF_MINUS_HALF(D%NIB:D%NIE)*0.5*(PZZ(D%NIB:D%NIE,KK+KKL)-PZZ(D%NIB:D%NIE,KK))
+ELSE
+  ZCOEFF_MINUS_HALF(D%NIB:D%NIE)=0.
+  ZTHV_MINUS_HALF(D%NIB:D%NIE) = PTHVM(D%NIB:D%NIE,KK)
+ENDIF
+ZCOEFF_PLUS_HALF(D%NIB:D%NIE)  = ((PTHVM(D%NIB:D%NIE,KK+KKL)-PTHVM(D%NIB:D%NIE,KK))/PDZZ(D%NIB:D%NIE,KK+KKL))
+ZTHV_PLUS_HALF(D%NIB:D%NIE)  = PTHVM(D%NIB:D%NIE,KK) + &
+                             & ZCOEFF_PLUS_HALF(D%NIB:D%NIE)*0.5*(PZZ(D%NIB:D%NIE,KK+KKL)-PZZ(D%NIB:D%NIE,KK))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 !               2  Dry part computation:
 !                  Integral buoyancy and computation of PENTR and PDETR for dry part
 !               --------------------------------------------------------------------
 
-DO JLOOP=1,SIZE(OTEST)
-  IF (OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.) THEN
+DO JI=D%NIB,D%NIE
+  IF (OTEST(JI) .AND. PPART_DRY(JI)>0.) THEN
     !Buoyancy computation in two parts to use change of gradient of theta v of environment
     !Between flux level KK and min(mass level, bottom of cloudy part)
-    ZDZ=MIN(ZDZ_STOP(JLOOP),(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*0.5)
-    PBUO_INTEG_DRY(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ*&
-                (0.5 * (  - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ  &
-                  - ZTHV_MINUS_HALF(JLOOP) + PTHV_UP(JLOOP) )
+    ZDZ=MIN(ZDZ_STOP(JI),(PZZ(JI,KK+KKL)-PZZ(JI,KK))*0.5)
+    PBUO_INTEG_DRY(JI) = ZG_O_THVREF_ED(JI)*ZDZ*&
+                (0.5 * (  - ZCOEFF_MINUS_HALF(JI))*ZDZ  &
+                  - ZTHV_MINUS_HALF(JI) + PTHV_UP(JI) )
 
     !Between mass flux KK and bottom of cloudy part (if above mass flux)
-    ZDZ=MAX(0., ZDZ_STOP(JLOOP)-(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*0.5)
-    PBUO_INTEG_DRY(JLOOP) = PBUO_INTEG_DRY(JLOOP) + ZG_O_THVREF(JLOOP)*ZDZ*&
-                (0.5 * (  - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ &
-                  - PTHVM(JLOOP,KK) + PTHV_UP(JLOOP) )
+    ZDZ=MAX(0., ZDZ_STOP(JI)-(PZZ(JI,KK+KKL)-PZZ(JI,KK))*0.5)
+    PBUO_INTEG_DRY(JI) = PBUO_INTEG_DRY(JI) + ZG_O_THVREF_ED(JI)*ZDZ*&
+                (0.5 * (  - ZCOEFF_PLUS_HALF(JI))*ZDZ &
+                  - PTHVM(JI,KK) + PTHV_UP(JI) )
 
     !Entr//Detr. computation
-    IF (PBUO_INTEG_DRY(JLOOP)>=0.) THEN
-      PENTR(JLOOP) = 0.5/(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)*&
-                 LOG(1.+ (2.*(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)/PW_UP2(JLOOP,KK))* &
-                 PBUO_INTEG_DRY(JLOOP))
-      PDETR(JLOOP) = 0.
+    IF (PBUO_INTEG_DRY(JI)>=0.) THEN
+      PENTR(JI) = 0.5/(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)*&
+                 LOG(1.+ (2.*(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)/PW_UP2(JI,KK))* &
+                 PBUO_INTEG_DRY(JI))
+      PDETR(JI) = 0.
     ELSE
-      PENTR(JLOOP) = 0.
-      PDETR(JLOOP) = 0.5/(PARAMMF%XABUO)*&
-                 LOG(1.+ (2.*(PARAMMF%XABUO)/PW_UP2(JLOOP,KK))* &
-                 (-PBUO_INTEG_DRY(JLOOP)))
+      PENTR(JI) = 0.
+      PDETR(JI) = 0.5/(PARAMMF%XABUO)*&
+                 LOG(1.+ (2.*(PARAMMF%XABUO)/PW_UP2(JI,KK))* &
+                 (-PBUO_INTEG_DRY(JI)))
     ENDIF
-    PENTR(JLOOP) = PARAMMF%XENTR_DRY*PENTR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))    
-    PDETR(JLOOP) = PARAMMF%XDETR_DRY*PDETR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
+    PENTR(JI) = PARAMMF%XENTR_DRY*PENTR(JI)/(PZZ(JI,KK+KKL)-PZZ(JI,KK))    
+    PDETR(JI) = PARAMMF%XDETR_DRY*PDETR(JI)/(PZZ(JI,KK+KKL)-PZZ(JI,KK))
     !Minimum value of detrainment
-    ZWK=PLUP(JLOOP)-0.5*(PZZ(JLOOP,KK)+PZZ(JLOOP,KK+KKL))
-    ZWK=SIGN(MAX(1., ABS(ZWK)), ZWK) ! ZWK must not be zero
-    PDETR(JLOOP) = MAX(PPART_DRY(JLOOP)*PARAMMF%XDETR_LUP/ZWK, PDETR(JLOOP))
+    ZWK0D=PLUP(JI)-0.5*(PZZ(JI,KK)+PZZ(JI,KK+KKL))
+    ZWK0D=SIGN(MAX(1., ABS(ZWK0D)), ZWK0D) ! ZWK0D must not be zero
+    PDETR(JI) = MAX(PPART_DRY(JI)*PARAMMF%XDETR_LUP/ZWK0D, PDETR(JI))
   ELSE
     !No dry part, condensation reached (OTESTLCL)
-    PBUO_INTEG_DRY(JLOOP) = 0.
-    PENTR(JLOOP)=0.
-    PDETR(JLOOP)=0.
+    PBUO_INTEG_DRY(JI) = 0.
+    PENTR(JI)=0.
+    PDETR(JI)=0.
   ENDIF
 ENDDO
 
@@ -863,44 +924,48 @@ ENDDO
 
 !               3.1 Integral buoyancy for cloudy part
 
-  ! Compute theta_v of updraft at flux level KK+KKL                   
-  !MIX variables are used to avoid declaring new variables
-  !but we are dealing with updraft and not mixture
-  ZRCMIX(:)=PRC_UP(:)
-  ZRIMIX(:)=PRI_UP(:)
-  CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
-               PPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
-               ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
-               ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
-               PBUF=ZBUF)
-  ZTHV_UP_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
-
-  ! Integral buoyancy for cloudy part
-  DO JLOOP=1,SIZE(OTEST)
-    IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)<1.) THEN
-      !Gradient of Theta V updraft over the cloudy part, assuming that thetaV updraft don't change
-      !between flux level KK and bottom of cloudy part
-      ZCOTHVU=(ZTHV_UP_F2(JLOOP)-PTHV_UP(JLOOP))/((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*(1-PPART_DRY(JLOOP)))
-
-      !Computation in two parts to use change of gradient of theta v of environment
-      !Between bottom of cloudy part (if under mass level) and mass level KK
-      ZDZ=MAX(0., 0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP))
-      PBUO_INTEG_CLD(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ*&
-              (0.5*( ZCOTHVU - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ &
-                - (PTHVM(JLOOP,KK)-ZDZ*ZCOEFF_MINUS_HALF(JLOOP)) + PTHV_UP(JLOOP) )
-
-      !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
-      ZDZ=(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-MAX(ZDZ_STOP(JLOOP),0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))
-      PBUO_INTEG_CLD(JLOOP) = PBUO_INTEG_CLD(JLOOP)+ZG_O_THVREF(JLOOP)*ZDZ*&
-                        (0.5*( ZCOTHVU - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ&
-                - (PTHVM(JLOOP,KK)+(0.5*((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))-ZDZ)*ZCOEFF_PLUS_HALF(JLOOP)) +&
-                PTHV_UP(JLOOP) )
+! Compute theta_v of updraft at flux level KK+KKL                   
+!MIX variables are used to avoid declaring new variables
+!but we are dealing with updraft and not mixture
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+ZRCMIX(D%NIB:D%NIE)=PRC_UP(D%NIB:D%NIE)
+ZRIMIX(D%NIB:D%NIE)=PRI_UP(D%NIB:D%NIE)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
+CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
+             PPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
+             ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
+             ZRSATW_ED, ZRSATI_ED,OOCEAN=.FALSE.,&
+             PBUF=ZBUF)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+ZTHV_UP_F2(D%NIB:D%NIE) = ZTHMIX(D%NIB:D%NIE)*(1.+ZRVORD*ZRVMIX(D%NIB:D%NIE))/(1.+PRT_UP(D%NIB:D%NIE))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
+
+! Integral buoyancy for cloudy part
+DO JI=D%NIB,D%NIE
+  IF(OTEST(JI) .AND. PPART_DRY(JI)<1.) THEN
+    !Gradient of Theta V updraft over the cloudy part, assuming that thetaV updraft don't change
+    !between flux level KK and bottom of cloudy part
+    ZCOTHVU=(ZTHV_UP_F2(JI)-PTHV_UP(JI))/((PZZ(JI,KK+KKL)-PZZ(JI,KK))*(1-PPART_DRY(JI)))
+
+    !Computation in two parts to use change of gradient of theta v of environment
+    !Between bottom of cloudy part (if under mass level) and mass level KK
+    ZDZ=MAX(0., 0.5*(PZZ(JI,KK+KKL)-PZZ(JI,KK))-ZDZ_STOP(JI))
+    PBUO_INTEG_CLD(JI) = ZG_O_THVREF_ED(JI)*ZDZ*&
+            (0.5*( ZCOTHVU - ZCOEFF_MINUS_HALF(JI))*ZDZ &
+              - (PTHVM(JI,KK)-ZDZ*ZCOEFF_MINUS_HALF(JI)) + PTHV_UP(JI) )
+
+    !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
+    ZDZ=(PZZ(JI,KK+KKL)-PZZ(JI,KK))-MAX(ZDZ_STOP(JI),0.5*(PZZ(JI,KK+KKL)-PZZ(JI,KK)))
+    PBUO_INTEG_CLD(JI) = PBUO_INTEG_CLD(JI)+ZG_O_THVREF_ED(JI)*ZDZ*&
+                      (0.5*( ZCOTHVU - ZCOEFF_PLUS_HALF(JI))*ZDZ&
+              - (PTHVM(JI,KK)+(0.5*((PZZ(JI,KK+KKL)-PZZ(JI,KK)))-ZDZ)*ZCOEFF_PLUS_HALF(JI)) +&
+              PTHV_UP(JI) )
 
-    ELSE
-      !No cloudy part
-      PBUO_INTEG_CLD(JLOOP)=0.
-    END IF
-  END DO
+  ELSE
+    !No cloudy part
+    PBUO_INTEG_CLD(JI)=0.
+  END IF
+END DO
 
 !               3.2 Critical mixed fraction for KK+KKL flux level (ZKIC_F2) and
 !                   for bottom of cloudy part (ZKIC), then a mean for the cloudy part
@@ -911,78 +976,82 @@ ENDDO
 !                   We determine the zero crossing of the linear curve
 !                   evaluating the derivative using ZMIXF=0.1
                 
-  ZKIC_INIT=0.1  ! starting value for critical mixed fraction for CLoudy Part
-
-  !  Compute thetaV of environment at the bottom of cloudy part
-  !    and cons then non cons. var. of mixture at the bottom of cloudy part
-
-  !   JI computed to avoid KKL(KK-KKL) being < KKL*KKB
-JI=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
-DO JLOOP=1,SIZE(OTEST)
-  IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.5) THEN
-    ZDZ=ZDZ_STOP(JLOOP)-0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
-    ZTHV(JLOOP)= PTHVM(JLOOP,KK)+ZCOEFF_PLUS_HALF(JLOOP)*ZDZ
-    ZMIXTHL(JLOOP) = ZKIC_INIT * &
-               (PTHLM(JLOOP,KK)+ZDZ*(PTHLM(JLOOP,KK+KKL)-PTHLM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) + &
-               (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
-    ZMIXRT(JLOOP)  = ZKIC_INIT * &
-               (PRTM(JLOOP,KK)+ZDZ*(PRTM(JLOOP,KK+KKL)-PRTM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) +   &
-               (1. - ZKIC_INIT)*PRT_UP(JLOOP)
-  ELSEIF(OTEST(JLOOP)) THEN
-    ZDZ=0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP)
-    ZTHV(JLOOP)= PTHVM(JLOOP,KK)-ZCOEFF_MINUS_HALF(JLOOP)*ZDZ
-    ZMIXTHL(JLOOP) = ZKIC_INIT * &
-               (PTHLM(JLOOP,KK)-ZDZ*(PTHLM(JLOOP,KK)-PTHLM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
-               (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
-    ZMIXRT(JLOOP)  = ZKIC_INIT * &
-               (PRTM(JLOOP,KK)-ZDZ*(PRTM(JLOOP,KK)-PRTM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
-               (1. - ZKIC_INIT)*PRT_UP(JLOOP)
+ZKIC_INIT=0.1  ! starting value for critical mixed fraction for CLoudy Part
+
+!  Compute thetaV of environment at the bottom of cloudy part
+!    and cons then non cons. var. of mixture at the bottom of cloudy part
+
+!   JKLIM computed to avoid KKL(KK-KKL) being < KKL*KKB
+JKLIM=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
+DO JI=D%NIB,D%NIE
+  IF(OTEST(JI) .AND. PPART_DRY(JI)>0.5) THEN
+    ZDZ=ZDZ_STOP(JI)-0.5*(PZZ(JI,KK+KKL)-PZZ(JI,KK))
+    ZTHV(JI)= PTHVM(JI,KK)+ZCOEFF_PLUS_HALF(JI)*ZDZ
+    ZMIXTHL(JI) = ZKIC_INIT * &
+               (PTHLM(JI,KK)+ZDZ*(PTHLM(JI,KK+KKL)-PTHLM(JI,KK))/PDZZ(JI,KK+KKL)) + &
+               (1. - ZKIC_INIT)*PTHL_UP(JI)
+    ZMIXRT(JI)  = ZKIC_INIT * &
+               (PRTM(JI,KK)+ZDZ*(PRTM(JI,KK+KKL)-PRTM(JI,KK))/PDZZ(JI,KK+KKL)) +   &
+               (1. - ZKIC_INIT)*PRT_UP(JI)
+  ELSEIF(OTEST(JI)) THEN
+    ZDZ=0.5*(PZZ(JI,KK+KKL)-PZZ(JI,KK))-ZDZ_STOP(JI)
+    ZTHV(JI)= PTHVM(JI,KK)-ZCOEFF_MINUS_HALF(JI)*ZDZ
+    ZMIXTHL(JI) = ZKIC_INIT * &
+               (PTHLM(JI,KK)-ZDZ*(PTHLM(JI,KK)-PTHLM(JI,JKLIM))/PDZZ(JI,KK)) + &
+               (1. - ZKIC_INIT)*PTHL_UP(JI)
+    ZMIXRT(JI)  = ZKIC_INIT * &
+               (PRTM(JI,KK)-ZDZ*(PRTM(JI,KK)-PRTM(JI,JKLIM))/PDZZ(JI,KK)) + &
+               (1. - ZKIC_INIT)*PRT_UP(JI)
   ELSE
 #ifdef REPRO55
-    ZMIXTHL(JLOOP) = 0.1
+    ZMIXTHL(JI) = 0.1
 #else
-    ZMIXTHL(JLOOP) = 300.
+    ZMIXTHL(JI) = 300.
 #endif
-    ZMIXRT(JLOOP) = 0.1
+    ZMIXRT(JI) = 0.1
   ENDIF
 ENDDO
-  CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
-               ZPRE,ZMIXTHL,ZMIXRT,&
-               ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
-               ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
-               PBUF=ZBUF)
-  ZTHVMIX(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
-
-  !  Compute cons then non cons. var. of mixture at the flux level KK+KKL  with initial ZKIC
-  ZMIXTHL(:) = ZKIC_INIT * 0.5*(PTHLM(:,KK)+PTHLM(:,KK+KKL))+(1. - ZKIC_INIT)*PTHL_UP(:)
-  ZMIXRT(:)  = ZKIC_INIT * 0.5*(PRTM(:,KK)+PRTM(:,KK+KKL))+(1. - ZKIC_INIT)*PRT_UP(:)
-  CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
-               PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
-               ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
-               ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
-               PBUF=ZBUF)
-  ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
-
-  !Computation of mean ZKIC over the cloudy part
-DO JLOOP=1,SIZE(OTEST)
-  IF (OTEST(JLOOP)) THEN
+CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
+             ZPRE,ZMIXTHL,ZMIXRT,&
+             ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
+             ZRSATW_ED, ZRSATI_ED,OOCEAN=.FALSE.,&
+             PBUF=ZBUF)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+ZTHVMIX(D%NIB:D%NIE) = ZTHMIX(D%NIB:D%NIE)*(1.+ZRVORD*ZRVMIX(D%NIB:D%NIE))/(1.+ZMIXRT(D%NIB:D%NIE))
+
+!  Compute cons then non cons. var. of mixture at the flux level KK+KKL  with initial ZKIC
+ZMIXTHL(D%NIB:D%NIE) = ZKIC_INIT * 0.5*(PTHLM(D%NIB:D%NIE,KK)+PTHLM(D%NIB:D%NIE,KK+KKL))+(1. - ZKIC_INIT)*PTHL_UP(D%NIB:D%NIE)
+ZMIXRT(D%NIB:D%NIE)  = ZKIC_INIT * 0.5*(PRTM(D%NIB:D%NIE,KK)+PRTM(D%NIB:D%NIE,KK+KKL))+(1. - ZKIC_INIT)*PRT_UP(D%NIB:D%NIE)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
+CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
+             PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
+             ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
+             ZRSATW_ED, ZRSATI_ED,OOCEAN=.FALSE.,&
+             PBUF=ZBUF)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
+ZTHVMIX_F2(D%NIB:D%NIE) = ZTHMIX(D%NIB:D%NIE)*(1.+ZRVORD*ZRVMIX(D%NIB:D%NIE))/(1.+ZMIXRT(D%NIB:D%NIE))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
+
+!Computation of mean ZKIC over the cloudy part
+DO JI=D%NIB,D%NIE
+  IF (OTEST(JI)) THEN
     ! Compute ZKIC at the bottom of cloudy part
     ! Thetav_up at bottom is equal to Thetav_up at flux level KK
-    IF (ABS(PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))<1.E-10) THEN
-      ZKIC(JLOOP)=1.
+    IF (ABS(PTHV_UP(JI)-ZTHVMIX(JI))<1.E-10) THEN
+      ZKIC(JI)=1.
     ELSE
-      ZKIC(JLOOP) = MAX(0.,PTHV_UP(JLOOP)-ZTHV(JLOOP))*ZKIC_INIT /  &  
-                 (PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))
+      ZKIC(JI) = MAX(0.,PTHV_UP(JI)-ZTHV(JI))*ZKIC_INIT /  &  
+                 (PTHV_UP(JI)-ZTHVMIX(JI))
     END IF
     ! Compute ZKIC_F2 at flux level KK+KKL
-    IF (ABS(ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))<1.E-10) THEN
-      ZKIC_F2(JLOOP)=1.
+    IF (ABS(ZTHV_UP_F2(JI)-ZTHVMIX_F2(JI))<1.E-10) THEN
+      ZKIC_F2(JI)=1.
     ELSE
-      ZKIC_F2(JLOOP) = MAX(0.,ZTHV_UP_F2(JLOOP)-ZTHV_PLUS_HALF(JLOOP))*ZKIC_INIT /  &  
-                 (ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))
+      ZKIC_F2(JI) = MAX(0.,ZTHV_UP_F2(JI)-ZTHV_PLUS_HALF(JI))*ZKIC_INIT /  &  
+                 (ZTHV_UP_F2(JI)-ZTHVMIX_F2(JI))
     END IF
     !Mean ZKIC over the cloudy part
-    ZKIC(JLOOP)=MAX(MIN(0.5*(ZKIC(JLOOP)+ZKIC_F2(JLOOP)),1.),0.)
+    ZKIC(JI)=MAX(MIN(0.5*(ZKIC(JI)+ZKIC_F2(JI)),1.),0.)
   END IF
 END DO
 
@@ -991,45 +1060,44 @@ END DO
 !                   in eq. (7) and (8) using eq. (5). Here we compute the ratio
 !                   of integrals without computing delta Me
 
-  !Constant PDF
-  !For this PDF, eq. (5) is delta Me=0.5*delta Mt
-DO JLOOP=1,SIZE(OTEST)
-  IF(OTEST(JLOOP)) THEN
-    ZEPSI(JLOOP) = ZKIC(JLOOP)**2. !integration multiplied by 2
-    ZDELTA(JLOOP) = (1.-ZKIC(JLOOP))**2. !idem
+!Constant PDF
+!For this PDF, eq. (5) is delta Me=0.5*delta Mt
+DO JI=D%NIB,D%NIE
+  IF(OTEST(JI)) THEN
+    ZEPSI(JI) = ZKIC(JI)**2. !integration multiplied by 2
+    ZDELTA(JI) = (1.-ZKIC(JI))**2. !idem
   ENDIF
 ENDDO
 
-  !Triangular PDF
-  !Calculus must be verified before activating this part, but in this state,
-  !results on ARM case are almost identical
-  !For this PDF, eq. (5) is also delta Me=0.5*delta Mt
-  !WHERE(OTEST)
-  !  !Integration multiplied by 2
-  !  WHERE(ZKIC<0.5)
-  !    ZEPSI(:)=8.*ZKIC(:)**3/3.
-  !    ZDELTA(:)=1.-4.*ZKIC(:)**2+8.*ZKIC(:)**3/3.
-  !  ELSEWHERE
-  !    ZEPSI(:)=5./3.-4*ZKIC(:)**2+8.*ZKIC(:)**3/3.
-  !    ZDELTA(:)=8.*(1.-ZKIC(:))**3/3.
-  !  ENDWHERE
-  !ENDWHERE
+!Triangular PDF
+!Calculus must be verified before activating this part, but in this state,
+!results on ARM case are almost identical
+!For this PDF, eq. (5) is also delta Me=0.5*delta Mt
+!WHERE(OTEST(D%NIB:D%NIE))
+!  !Integration multiplied by 2
+!  WHERE(ZKIC<0.5)
+!    ZEPSI(D%NIB:D%NIE)=8.*ZKIC(D%NIB:D%NIE)**3/3.
+!    ZDELTA(D%NIB:D%NIE)=1.-4.*ZKIC(D%NIB:D%NIE)**2+8.*ZKIC(D%NIB:D%NIE)**3/3.
+!  ELSEWHERE
+!    ZEPSI(D%NIB:D%NIE)=5./3.-4*ZKIC(D%NIB:D%NIE)**2+8.*ZKIC(D%NIB:D%NIE)**3/3.
+!    ZDELTA(D%NIB:D%NIE)=8.*(1.-ZKIC(D%NIB:D%NIE))**3/3.
+!  ENDWHERE
+!ENDWHERE
 
 !               3.4 Computation of PENTR and PDETR
-DO JLOOP=1,SIZE(OTEST)
-  IF(OTEST(JLOOP)) THEN
-    ZEPSI_CLOUD=MIN(ZDELTA(JLOOP), ZEPSI(JLOOP))
-    PENTR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZEPSI_CLOUD
-    PDETR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZDELTA(JLOOP)
-    PENTR(JLOOP) = PENTR(JLOOP)+PENTR_CLD(JLOOP)
-    PDETR(JLOOP) = PDETR(JLOOP)+PDETR_CLD(JLOOP)
+DO JI=D%NIB,D%NIE
+  IF(OTEST(JI)) THEN
+    ZEPSI_CLOUD=MIN(ZDELTA(JI), ZEPSI(JI))
+    PENTR_CLD(JI) = (1.-PPART_DRY(JI))*ZCOEFFMF_CLOUD*PRHODREF(JI)*ZEPSI_CLOUD
+    PDETR_CLD(JI) = (1.-PPART_DRY(JI))*ZCOEFFMF_CLOUD*PRHODREF(JI)*ZDELTA(JI)
+    PENTR(JI) = PENTR(JI)+PENTR_CLD(JI)
+    PDETR(JI) = PDETR(JI)+PDETR_CLD(JI)
   ELSE
-    PENTR_CLD(JLOOP) = 0.
-    PDETR_CLD(JLOOP) = 0.
+    PENTR_CLD(JI) = 0.
+    PDETR_CLD(JI) = 0.
   ENDIF
 ENDDO
 
-IF (LHOOK) CALL DR_HOOK('COMPUTE_ENTR_DETR',1,ZHOOK_HANDLE)
 END SUBROUTINE COMPUTE_ENTR_DETR
 END SUBROUTINE COMPUTE_UPDRAFT
 END MODULE MODE_COMPUTE_UPDRAFT
diff --git a/src/common/turb/mode_compute_updraft_raha.F90 b/src/common/turb/mode_compute_updraft_raha.F90
index f6818168e1722e9383e6aa1469a6e2c25f2f6a08..a17f1deade8c192ef983e472001547a2c2497504 100644
--- a/src/common/turb/mode_compute_updraft_raha.F90
+++ b/src/common/turb/mode_compute_updraft_raha.F90
@@ -230,13 +230,15 @@ ZTH_UP(:,:)=0.
 PFRAC_UP(:,:)=0.
 PTHV_UP(:,:)=0.
 
-PBUO_INTEG=0.
-ZBUO      =0.
+PBUO_INTEG(:,:)=0.
+ZBUO(:,:)      =0.
 
 !no ice cloud coded yet 
 PRI_UP(:,:)=0.
 PFRAC_ICE_UP(:,:)=0.
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not used
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 ! Initialisation of environment variables at t-dt
 
@@ -254,10 +256,12 @@ CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 !END DO  
 
 !          Initialisation of updraft characteristics 
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PTHL_UP(:,:)=ZTHLM_F(:,:)
 PRT_UP(:,:)=ZRTM_F(:,:)
 PU_UP(:,:)=ZUM_F(:,:)
 PV_UP(:,:)=ZVM_F(:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PSV_UP(:,:,:)=0.
 !IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) then
 !    PSV_UP(:,:,:)=ZSVM_F(:,:,:)
@@ -266,48 +270,58 @@ PSV_UP(:,:,:)=0.
 ! Computation or initialisation of updraft characteristics at the KKB level
 ! thetal_up,rt_up,thetaV_up, w�,Buoyancy term and mass flux (PEMF)
 
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 PTHL_UP(:,D%NKB)= ZTHLM_F(:,D%NKB)+MAX(0.,MIN(ZTMAX,(PSFTH(:)/SQRT(ZTKEM_F(:,D%NKB)))*PARAMMF%XALP_PERT))
 PRT_UP(:,D%NKB) = ZRTM_F(:,D%NKB)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,D%NKB)))*PARAMMF%XALP_PERT)) 
 
 ZQT_UP(:) = PRT_UP(:,D%NKB)/(1.+PRT_UP(:,D%NKB))
 ZTHS_UP(:,D%NKB)=PTHL_UP(:,D%NKB)*(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(:))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 CALL MZM_MF(D, PTHM (:,:), ZTHM_F(:,:))
 CALL MZM_MF(D, PPABSM(:,:), ZPRES_F(:,:))
 CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F(:,:))
 CALL MZM_MF(D, PRVM(:,:), ZRVM_F(:,:))
 
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ! thetav at mass and flux levels 
 ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
 ZTHVM(:,:)=PTHM(:,:)*((1.+ZRVORD*PRVM(:,:))/(1.+PRTM(:,:)))
 
 PTHV_UP(:,:)= ZTHVM_F(:,:)
-PRV_UP (:,:)= ZRVM_F (:,:)
+PRV_UP(:,:) = ZRVM_F(:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 ZW_UP2(:,:)=ZEPS
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 ZW_UP2(:,D%NKB) = MAX(0.0001,(1./6.)*ZTKEM_F(:,D%NKB))
-GTEST = (ZW_UP2(:,D%NKB) > ZEPS)  
+GTEST(:) = (ZW_UP2(:,D%NKB) > ZEPS)  
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 ! Computation of non conservative variable for the KKB level of the updraft
 ! (all or nothing ajustement)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 PRC_UP(:,D%NKB)=0.
 PRI_UP(:,D%NKB)=0.
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
              PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
              PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
              PBUF=ZBUF)
 
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 ! compute updraft thevav and buoyancy term at KKB level             
 PTHV_UP(:,D%NKB) = ZTH_UP(:,D%NKB)*((1+ZRVORD*PRV_UP(:,D%NKB))/(1+PRT_UP(:,D%NKB))) 
 ! compute mean rsat in updraft
 PRSAT_UP(:,D%NKB) = ZRSATW(:)*(1-PFRAC_ICE_UP(:,D%NKB)) + ZRSATI(:)*PFRAC_ICE_UP(:,D%NKB)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 !Tout est commente pour tester dans un premier temps la s�paration en deux de la 
 !  boucle verticale, une pour w et une pour PEMF                                                            
-
-ZG_O_THVREF=CST%XG/ZTHVM_F
-
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+ZG_O_THVREF(:,:)=CST%XG/ZTHVM_F(:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 !  Definition de l'alimentation au sens de la fermeture de Hourdin et al
 
@@ -316,6 +330,7 @@ ZALIM_STAR_TOT(:) = 0.    ! <== Normalization of ZALIM_STAR
 IALIM(:)          = D%NKB   ! <== Top level of the alimentation layer
 
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL   !  Vertical loop
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
   ZZDZ(:,JK)   = MAX(ZEPS,PZZ(:,JK+D%NKL)-PZZ(:,JK))       ! <== Delta Z between two flux level
   ZZZ(:,JK)    = MAX(0.,0.5*(PZZ(:,JK+D%NKL)+PZZ(:,JK)) )  ! <== Hight of mass levels
   ZDTHETASDZ(:,JK) = (ZTHVM_F(:,JK)-ZTHVM_F(:,JK+D%NKL))   ! <== Delta theta_v
@@ -325,13 +340,16 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL   !  Vertical loop
      ZALIM_STAR_TOT(:) = ZALIM_STAR_TOT(:)+ZALIM_STAR(:,JK)*ZZDZ(:,JK)
      IALIM(:)          = JK
   ENDWHERE
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO
 
 ! Normalization of ZALIM_STAR
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL   !  Vertical loop   
-   WHERE (ZALIM_STAR_TOT > ZEPS)
+   !$mnh_expand_where(JI=D%NIB:D%NIE)
+   WHERE (ZALIM_STAR_TOT(:) > ZEPS)
       ZALIM_STAR(:,JK)  = ZALIM_STAR(:,JK)/ZALIM_STAR_TOT(:)
    ENDWHERE
+   !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO
 ZALIM_STAR_TOT(:) = 0.
 
@@ -358,7 +376,7 @@ ZPHI(:) = 0.
 
 
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL
-
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
 ! IF the updraft top is reached for all column, stop the loop on levels
 
 !  ITEST=COUNT(GTEST)
@@ -369,12 +387,10 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
 
 ! to find the LCL (check if JK is LCL or not)
-
-  WHERE ((PRC_UP(:,JK)+PRI_UP(:,JK)>0.).AND.(.NOT.(GTESTLCL)))
+    WHERE ((PRC_UP(:,JK)+PRI_UP(:,JK)>0.).AND.(.NOT.(GTESTLCL(:))))
       KKLCL(:) = JK           
       GTESTLCL(:)=.TRUE.
-  ENDWHERE
-
+    ENDWHERE
 
 ! COMPUTE PENTR and PDETR at mass level JK
 
@@ -386,7 +402,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
     ZRC_UP(:)  = PRC_UP(:,JK)
     ZRI_UP(:)  = PRI_UP(:,JK) ! guess 
     ZRV_UP(:)  = PRV_UP(:,JK)
-    ZBUO      (:,JK) = ZG_O_THVREF(:,JK)*(PTHV_UP(:,JK) - ZTHVM_F(:,JK))
+    ZBUO(:,JK) = ZG_O_THVREF(:,JK)*(PTHV_UP(:,JK) - ZTHVM_F(:,JK))
     PBUO_INTEG(:,JK) = ZBUO(:,JK)*(PZZ(:,JK+D%NKL)-PZZ(:,JK))
     
     ZDZ(:)   = MAX(ZEPS,PZZ(:,JK+D%NKL)-PZZ(:,JK))
@@ -416,7 +432,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
    
 ! If the updraft did not stop, compute cons updraft characteritics at jk+1
-  WHERE(GTEST)     
+  WHERE(GTEST(:))     
     ZZTOP(:) = MAX(ZZTOP(:),PZZ(:,JK+D%NKL))
     ZMIX2(:) = (PZZ(:,JK+D%NKL)-PZZ(:,JK))*PENTR(:,JK) !&
     ZMIX3(:) = (PZZ(:,JK+D%NKL)-PZZ(:,JK))*PDETR(:,JK) !&                   
@@ -425,7 +441,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
     ZTHSM(:,JK) = PTHLM(:,JK)*(1.+PARAMMF%XLAMBDA_MF*ZQTM(:))
     ZTHS_UP(:,JK+D%NKL)=(ZTHS_UP(:,JK)*(1.-0.5*ZMIX2(:)) + ZTHSM(:,JK)*ZMIX2(:)) &
                           /(1.+0.5*ZMIX2(:))   
-    PRT_UP(:,JK+D%NKL)=(PRT_UP (:,JK)*(1.-0.5*ZMIX2(:)) + PRTM(:,JK)*ZMIX2(:))  &
+    PRT_UP(:,JK+D%NKL)=(PRT_UP(:,JK)*(1.-0.5*ZMIX2(:)) + PRTM(:,JK)*ZMIX2(:))  &
                           /(1.+0.5*ZMIX2(:))
     ZQT_UP(:) = PRT_UP(:,JK+D%NKL)/(1.+PRT_UP(:,JK+D%NKL))
     PTHL_UP(:,JK+D%NKL)=ZTHS_UP(:,JK+D%NKL)/(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(:))                      
@@ -434,25 +450,25 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
   IF(OMIXUV) THEN
     IF(JK/=D%NKB) THEN
-      WHERE(GTEST)
-        PU_UP(:,JK+D%NKL) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
+      WHERE(GTEST(:))
+        PU_UP(:,JK+D%NKL) = (PU_UP(:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
                           0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
                           ((PUM(:,JK+D%NKL)-PUM(:,JK))/PDZZ(:,JK+D%NKL)+&
                            (PUM(:,JK)-PUM(:,JK-D%NKL))/PDZZ(:,JK))        )   &
                           /(1+0.5*ZMIX2(:))
-        PV_UP(:,JK+D%NKL) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
+        PV_UP(:,JK+D%NKL) = (PV_UP(:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
                           0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
                           ((PVM(:,JK+D%NKL)-PVM(:,JK))/PDZZ(:,JK+D%NKL)+&
                            (PVM(:,JK)-PVM(:,JK-D%NKL))/PDZZ(:,JK))    )   &
                           /(1+0.5*ZMIX2(:))
       ENDWHERE
     ELSE
-      WHERE(GTEST)
-        PU_UP(:,JK+D%NKL) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
+      WHERE(GTEST(:))
+        PU_UP(:,JK+D%NKL) = (PU_UP(:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
                           0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
                           ((PUM(:,JK+D%NKL)-PUM(:,JK))/PDZZ(:,JK+D%NKL))        )   &
                           /(1+0.5*ZMIX2(:))
-        PV_UP(:,JK+D%NKL) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
+        PV_UP(:,JK+D%NKL) = (PV_UP(:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
                           0.5*PARAMMF%XPRES_UV*(PZZ(:,JK+D%NKL)-PZZ(:,JK))*&
                           ((PVM(:,JK+D%NKL)-PVM(:,JK))/PDZZ(:,JK+D%NKL))    )   &
                           /(1+0.5*ZMIX2(:))
@@ -463,7 +479,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 !  DO JSV=1,ISV 
 !     IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
 !      WHERE(GTEST) 
-!           PSV_UP(:,JK+KKL,JSV) = (PSV_UP (:,JK,JSV)*(1-0.5*ZMIX2(:)) + &
+!           PSV_UP(:,JK+KKL,JSV) = (PSV_UP(:,JK,JSV)*(1-0.5*ZMIX2(:)) + &
 !                        PSVM(:,JK,JSV)*ZMIX2(:))  /(1+0.5*ZMIX2(:))
 !      ENDWHERE                        
 !  ENDDO  
@@ -473,11 +489,13 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
   ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
   ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
   ZRV_UP(:)=PRV_UP(:,JK)
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
   CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
           PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL),              &
           ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
           PBUF=ZBUF)
-  WHERE(GTEST)
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
+  WHERE(GTEST(:))
     ZT_UP(:) = ZTH_UP(:,JK+D%NKL)*PEXNM(:,JK+D%NKL)
     ZCP(:) = CST%XCPD + CST%XCL * ZRC_UP(:)
     ZLVOCPEXN(:)=(CST%XLVTT + (CST%XCPV-CST%XCL) *  (ZT_UP(:)-CST%XTT) ) / ZCP(:) / PEXNM(:,JK+D%NKL)
@@ -491,22 +509,22 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
   
 
 ! Compute the updraft theta_v, buoyancy and w**2 for level JK+1   
- WHERE(GTEST)
+  WHERE(GTEST(:))
 !      PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL)))
       PTHV_UP(:,JK+D%NKL) = ZTH_UP(:,JK+D%NKL)*(1.+0.608*PRV_UP(:,JK+D%NKL) - PRC_UP(:,JK+D%NKL))
- ENDWHERE
+  ENDWHERE
 
 
 ! Test if the updraft has reach the ETL
   GTESTETL(:)=.FALSE.
-  WHERE (GTEST.AND.(PBUO_INTEG(:,JK)<=0.))
+  WHERE (GTEST(:).AND.(PBUO_INTEG(:,JK)<=0.))
       KKETL(:) = JK+D%NKL
       GTESTETL(:)=.TRUE.
   ENDWHERE
 
 ! Test is we have reached the top of the updraft
 
-  WHERE (GTEST.AND.((ZW_UP2(:,JK+D%NKL)<=ZEPS)))
+  WHERE (GTEST(:).AND.((ZW_UP2(:,JK+D%NKL)<=ZEPS)))
       ZW_UP2(:,JK+D%NKL)=ZEPS
       GTEST(:)=.FALSE.
       PTHL_UP(:,JK+D%NKL)=ZTHLM_F(:,JK+D%NKL)
@@ -518,24 +536,29 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
       PFRAC_UP(:,JK+D%NKL)=0.
       KKCTL(:)=JK+D%NKL
   ENDWHERE
-
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO 
 
 ! Closure assumption for mass flux at KKB+1 level (Mass flux is supposed to be 0 at KKB level !)                                                 
 !   Hourdin et al 2002 formulation
 
 
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 ZZTOP(:) = MAX(ZZTOP(:),ZEPS)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 DO JK=D%NKB+D%NKL,D%NKE-D%NKL,D%NKL !  Vertical loop
-   WHERE(JK<=IALIM)
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
+   WHERE(JK<=IALIM(:))
      ZALIM_STAR_TOT(:) = ZALIM_STAR_TOT(:) + ZALIM_STAR(:,JK)*ZALIM_STAR(:,JK)*ZZDZ(:,JK)/PRHODREF(:,JK)
-   ENDWHERE  
+   ENDWHERE
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO   
 
-WHERE (ZALIM_STAR_TOT*ZZTOP > ZEPS)
+!$mnh_expand_where(JI=D%NIB:D%NIE)
+WHERE (ZALIM_STAR_TOT(:)*ZZTOP(:) > ZEPS)
    ZPHI(:) =  ZW_MAX(:)/(PARAMMF%XR*ZZTOP(:)*ZALIM_STAR_TOT(:))
-ENDWHERE      
+ENDWHERE
 
 GTEST(:) = .TRUE.
 PEMF(:,D%NKB+D%NKL) = ZPHI(:)*ZZDZ(:,D%NKB)*ZALIM_STAR(:,D%NKB)
@@ -543,13 +566,15 @@ PEMF(:,D%NKB+D%NKL) = ZPHI(:)*ZZDZ(:,D%NKB)*ZALIM_STAR(:,D%NKB)
 PFRAC_UP(:,D%NKB+D%NKL)=PEMF(:,D%NKB+D%NKL)/(SQRT(ZW_UP2(:,D%NKB+D%NKL))*ZRHO_F(:,D%NKB+D%NKL))
 PFRAC_UP(:,D%NKB+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX,PFRAC_UP(:,D%NKB+D%NKL))
 PEMF(:,D%NKB+D%NKL) = ZRHO_F(:,D%NKB+D%NKL)*PFRAC_UP(:,D%NKB+D%NKL)*SQRT(ZW_UP2(:,D%NKB+D%NKL))
+!$mnh_end_expand_where(JI=D%NIB:D%NIE)
 
 DO JK=D%NKB+D%NKL,D%NKE-D%NKL,D%NKL !  Vertical loop
+  !$mnh_expand_where(JI=D%NIB:D%NIE)
      
-   GTEST = (ZW_UP2(:,JK) > ZEPS)  
+  GTEST(:) = (ZW_UP2(:,JK) > ZEPS)  
 
-  WHERE (GTEST)
-    WHERE(JK<IALIM)
+  WHERE (GTEST(:))
+    WHERE(JK<IALIM(:))
        PEMF(:,JK+D%NKL) = MAX(0.,PEMF(:,JK) + ZPHI(:)*ZZDZ(:,JK)*(PENTR(:,JK) - PDETR(:,JK)))
     ELSEWHERE
        ZMIX1(:)=ZZDZ(:,JK)*(PENTR(:,JK)-PDETR(:,JK))
@@ -561,11 +586,15 @@ DO JK=D%NKB+D%NKL,D%NKE-D%NKL,D%NKL !  Vertical loop
     PFRAC_UP(:,JK+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX,PFRAC_UP(:,JK+D%NKL))
     PEMF(:,JK+D%NKL) = ZRHO_F(:,JK+D%NKL)*PFRAC_UP(:,JK+D%NKL)*SQRT(ZW_UP2(:,JK+D%NKL))
   ENDWHERE
-
+  !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO
 
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PW_UP(:,:)=SQRT(ZW_UP2(:,:))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 PEMF(:,D%NKB) =0.
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 ! Limits the shallow convection scheme when cloud heigth is higher than 3000m.
 ! To do this, mass flux is multiplied by a coefficient decreasing linearly
@@ -573,18 +602,26 @@ PEMF(:,D%NKB) =0.
 ! This way, all MF fluxes are diminished by this amount.
 ! Diagnosed cloud fraction is also multiplied by the same coefficient.
 !
-DO JI=1,SIZE(PTHM,1) 
+DO JI=D%NIB,D%NIE 
    PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
 END DO
 
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
-GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=D%NKT )
-ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=D%NKT)
-ZCOEF=MIN(MAX(ZCOEF,0.),1.)
-WHERE (GWORK2) 
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
+DO JK=D%NKTB,D%NKTE
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  GWORK2(:,JK) = GWORK1(:)
+  ZCOEF(:,JK) = (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1))
+  ZCOEF(:,JK)=MIN(MAX(ZCOEF(:,JK),0.),1.)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
+ENDDO
+!$mnh_expand_where(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+WHERE (GWORK2(:,:)) 
    PEMF(:,:)     = PEMF(:,:)     * ZCOEF(:,:)
    PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
 ENDWHERE
+!$mnh_end_expand_where(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAF_RAHA',1,ZHOOK_HANDLE)
 !
diff --git a/src/common/turb/mode_compute_updraft_rhcj10.F90 b/src/common/turb/mode_compute_updraft_rhcj10.F90
index b5bfe8cf2f6965333331f8e898efee07a5c25f52..9e97f76bbae74c6f10a46cd438daece65a8cb100 100644
--- a/src/common/turb/mode_compute_updraft_rhcj10.F90
+++ b/src/common/turb/mode_compute_updraft_rhcj10.F90
@@ -242,7 +242,9 @@ ZBUO      =0.
 !no ice cloud coded yet
 PRI_UP(:,:)=0.
 PFRAC_ICE_UP(:,:)=0.
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not used
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 ! Initialisation of environment variables at t-dt
 
@@ -264,10 +266,12 @@ CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 !END DO
 
 !          Initialisation of updraft characteristics 
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PTHL_UP(:,:)=ZTHLM_F(:,:)
 PRT_UP(:,:)=ZRTM_F(:,:)
 PU_UP(:,:)=ZUM_F(:,:)
 PV_UP(:,:)=ZVM_F(:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PSV_UP(:,:,:)=0.
 ! This updraft is not yet ready to use scalar variables
 !IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) then
@@ -277,7 +281,7 @@ PSV_UP(:,:,:)=0.
 ! Computation or initialisation of updraft characteristics at the KKB level
 ! thetal_up,rt_up,thetaV_up, w,Buoyancy term and mass flux (PEMF)
 
-DO JI=1,D%NIT
+DO JI=D%NIB,D%NIE
   !PTHL_UP(JI,KKB)= ZTHLM_F(JI,KKB)+MAX(0.,MIN(ZTMAX,(PSFTH(JI)/SQRT(ZTKEM_F(JI,KKB)))*XALP_PERT))
   !PRT_UP(JI,KKB) = ZRTM_F(JI,KKB)+MAX(0.,MIN(ZRMAX,(PSFRV(JI)/SQRT(ZTKEM_F(JI,KKB)))*XALP_PERT)) 
   PTHL_UP(JI,D%NKB)= ZTHLM_F(JI,D%NKB)
@@ -292,30 +296,36 @@ CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F(:,:))
 CALL MZM_MF(D, PRVM(:,:), ZRVM_F(:,:))
 
 ! thetav at mass and flux levels 
-DO JK=1,D%NKT
-  DO JI=1,D%NIT
+DO JK=D%NKTB,D%NKTE
+  DO JI=d%NIB,D%NIE
     ZTHVM_F(JI,JK)=ZTHM_F(JI,JK)*((1.+ZRVORD*ZRVM_F(JI,JK))/(1.+ZRTM_F(JI,JK)))
   ENDDO
 ENDDO
 
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PTHV_UP(:,:)= ZTHVM_F(:,:)
-PRV_UP (:,:)= ZRVM_F (:,:)
+PRV_UP(:,:)= ZRVM_F(:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 ZW_UP2(:,:)=ZEPS
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 !ZW_UP2(:,KKB) = MAX(0.0001,(3./6.)*ZTKEM_F(:,KKB))
 ZW_UP2(:,D%NKB) = MAX(0.0001,(2./3.)*ZTKEM_F(:,D%NKB))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 ! Computation of non conservative variable for the KKB level of the updraft
 ! (all or nothing ajustement)
 
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 PRC_UP(:,D%NKB)=0.
 PRI_UP(:,D%NKB)=0.
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
              PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
              PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
              PBUF=ZBUF)
 
-DO JI=1,D%NIT
+DO JI=D%NIB,D%NIE
   ! compute updraft thevav and buoyancy term at KKB level             
   PTHV_UP(JI,D%NKB) = ZTH_UP(JI,D%NKB)*((1+ZRVORD*PRV_UP(JI,D%NKB))/(1+PRT_UP(JI,D%NKB))) 
   ! compute mean rsat in updraft
@@ -325,35 +335,43 @@ ENDDO
 !Tout est commente pour tester dans un premier temps la separation en deux de la 
 !  boucle verticale, une pour w et une pour PEMF                                                            
 
-ZG_O_THVREF=CST%XG/ZTHVM_F
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+ZG_O_THVREF(:,:)=CST%XG/ZTHVM_F(:,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 ! Calcul de la fermeture de Julien Pergaut comme limite max de PHY
 
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL   !  Vertical loop
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIE
     ZZDZ(JI,JK)   = MAX(ZEPS,PZZ(JI,JK+D%NKL)-PZZ(JI,JK))  ! <== Delta Z between two flux level
   ENDDO
 ENDDO
 
 ! compute L_up
 GLMIX=.TRUE.
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 ZTKEM_F(:,D%NKB)=0.
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 !
 IF(TURB%CTURBLEN=='RM17') THEN
   CALL GZ_M_W_MF(D, PUM, PDZZ, ZWK)
   CALL MZF_MF(D, ZWK, ZDUDZ)
   CALL GZ_M_W_MF(D, PVM, PDZZ, ZWK)
   CALL MZF_MF(D, ZWK, ZDVDZ)
-  ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+  ZSHEAR(:,:) = SQRT(ZDUDZ(:,:)**2 + ZDVDZ(:,:)**2)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ELSE
-  ZSHEAR = 0. !no shear in bl89 mixing length
+  ZSHEAR(:,:) = 0. !no shear in bl89 mixing length
 END IF
 !
 CALL COMPUTE_BL89_ML(D, CST, CSTURB, PDZZ,ZTKEM_F(:,D%NKB),ZG_O_THVREF(:,D%NKB), &
                        ZTHVM_F,D%NKB,GLMIX,.TRUE.,ZSHEAR,ZLUP)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 ZLUP(:)=MAX(ZLUP(:),1.E-10)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
-DO JI=1,D%NIT
+DO JI=D%NIB,D%NIE
   ! Compute Buoyancy flux at the ground
   ZWTHVSURF = (ZTHVM_F(JI,D%NKB)/ZTHM_F(JI,D%NKB))*PSFTH(JI)+     &
               (0.61*ZTHM_F(JI,D%NKB))*PSFRV(JI)
@@ -401,7 +419,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
  
 ! to find the LCL (check if JK is LCL or not)
 
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIE
     IF ((PRC_UP(JI,JK)+PRI_UP(JI,JK)>0.).AND.(.NOT.(GTESTLCL(JI)))) THEN
       KKLCL(JI) = JK           
       GTESTLCL(JI)=.TRUE.
@@ -416,18 +434,20 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
   ! Compute theta_v of updraft at flux level JK    
     
+    !$mnh_expand_array(JI=D%NIB:D%NIE)
     ZRC_UP(:)   =PRC_UP(:,JK) ! guess
     ZRI_UP(:)   =PRI_UP(:,JK) ! guess 
     ZRV_UP(:)   =PRV_UP(:,JK)
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE)
     CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
                PPABSM(:,JK),PTHL_UP(:,JK),PRT_UP(:,JK),&
                ZTH_UP(:,JK),ZRV_UP,ZRC_UP,ZRI_UP,ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
                PBUF=ZBUF)
     
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIE
     IF (GTEST(JI)) THEN
-      PTHV_UP   (JI,JK) = ZTH_UP(JI,JK)*(1.+ZRVORD*ZRV_UP(JI))/(1.+PRT_UP(JI,JK))
-      ZBUO      (JI,JK) = ZG_O_THVREF(JI,JK)*(PTHV_UP(JI,JK) - ZTHVM_F(JI,JK))    
+      PTHV_UP(JI,JK)    = ZTH_UP(JI,JK)*(1.+ZRVORD*ZRV_UP(JI))/(1.+PRT_UP(JI,JK))
+      ZBUO(JI,JK)       = ZG_O_THVREF(JI,JK)*(PTHV_UP(JI,JK) - ZTHVM_F(JI,JK))    
       PBUO_INTEG(JI,JK) = ZBUO(JI,JK)*(PZZ(JI,JK+D%NKL)-PZZ(JI,JK))
       
       ZDZ(JI)   = MAX(ZEPS,PZZ(JI,JK+D%NKL)-PZZ(JI,JK))
@@ -476,7 +496,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
   IF(OMIXUV) THEN
     IF(JK/=D%NKB) THEN
-      DO JI=1,D%NIT
+      DO JI=D%NIB,D%NIE
         IF(GTEST(JI)) THEN
           PU_UP(JI,JK+D%NKL) = (PU_UP (JI,JK)*(1-0.5*ZMIX2(JI)) + PUM(JI,JK)*ZMIX2(JI)+ &
                             0.5*PARAMMF%XPRES_UV*(PZZ(JI,JK+D%NKL)-PZZ(JI,JK))*&
@@ -491,7 +511,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
         ENDIF
       ENDDO
     ELSE
-      DO JI=1,D%NIT
+      DO JI=D%NIB,D%NIE
         IF(GTEST(JI)) THEN
           PU_UP(JI,JK+D%NKL) = (PU_UP (JI,JK)*(1-0.5*ZMIX2(JI)) + PUM(JI,JK)*ZMIX2(JI)+ &
                             0.5*PARAMMF%XPRES_UV*(PZZ(JI,JK+D%NKL)-PZZ(JI,JK))*&
@@ -516,16 +536,18 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 !  ENDDO  
   
   
-! Compute non cons. var. at level JK+KKL
+  ! Compute non cons. var. at level JK+KKL
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
   ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
   ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
   ZRV_UP(:)=PRV_UP(:,JK)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
   CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
           PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL),              &
           ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
           PBUF=ZBUF)
 
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIE
     IF(GTEST(JI)) THEN
       !ZT_UP(JI) = ZTH_UP(JI,JK+KKL)*PEXNM(JI,JK+KKL)
       !ZCP(JI) = XCPD + XCL * ZRC_UP(JI)
@@ -548,13 +570,13 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
     ENDIF
   ENDDO
 
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIB
     IF(GTEST(JI)) THEN
       PEMF(JI,JK+D%NKL)=PEMF(JI,JK)*EXP(ZMIX1(JI))
     ENDIF
   ENDDO
 
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIE
     IF(GTEST(JI)) THEN
       ! Updraft fraction must be smaller than XFRAC_UP_MAX
       PFRAC_UP(JI,JK+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX, &
@@ -564,7 +586,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
   ENDDO
 
 ! Test if the updraft has reach the ETL
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIE
     IF (GTEST(JI) .AND. (PBUO_INTEG(JI,JK)<=0.)) THEN
       KKETL(JI) = JK+D%NKL
     ENDIF
@@ -572,7 +594,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
 
 ! Test is we have reached the top of the updraft
-  DO JI=1,D%NIT
+  DO JI=D%NIB,D%NIB
     IF (GTEST(JI) .AND. ((ZW_UP2(JI,JK+D%NKL)<=ZEPS).OR.(PEMF(JI,JK+D%NKL)<=ZEPS))) THEN
       ZW_UP2   (JI,JK+D%NKL)=ZEPS
       PEMF     (JI,JK+D%NKL)=0.
@@ -590,8 +612,12 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 
 ENDDO   ! Fin de la boucle verticale 
 
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 PW_UP(:,:)=SQRT(ZW_UP2(:,:))
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 PEMF(:,D%NKB) =0.
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 ! Limits the shallow convection scheme when cloud heigth is higher than 3000m.
 ! To do this, mass flux is multiplied by a coefficient decreasing linearly
@@ -599,16 +625,22 @@ PEMF(:,D%NKB) =0.
 ! This way, all MF fluxes are diminished by this amount.
 ! Diagnosed cloud fraction is also multiplied by the same coefficient.
 !
-DO JI=1,D%NIT
+DO JI=D%NIB,D%NIE
    PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
 ENDDO
 
+!$mnh_expand_array(JI=D%NIB:D%NIE)
 GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
-GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=D%NKT )
-ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=D%NKT)
-ZCOEF(:,:)=MIN(MAX(ZCOEF(:,:),0.),1.)
-DO JK=1, D%NKT
-  DO JI=1,D%NIT
+!$mnh_end_expand_array(JI=D%NIB:D%NIE)
+DO JK=D%NKTB,D%NKTE
+  !$mnh_expand_array(JI=D%NIB:D%NIE)
+  GWORK2(:,JK) = GWORK1(:)
+  ZCOEF(:,JK) = (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1))
+  ZCOEF(:,JK)=MIN(MAX(ZCOEF(:,JK),0.),1.)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE)
+ENDDO
+DO JK=D%NKTB, D%NKTE
+  DO JI=D%NIB,D%NIE
     IF (GWORK2(JI,JK)) THEN
       PEMF(JI,JK)     = PEMF(JI,JK)     * ZCOEF(JI,JK) 
       PFRAC_UP(JI,JK) = PFRAC_UP(JI,JK) * ZCOEF(JI,JK)