diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..a64c2f8df121e44db5c0f6dfa98bb4e6ad08b669
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+.*.swp
+.*.swo
diff --git a/docs/TODO b/docs/TODO
index 3c75ded6b079d8541595e7be7a74e073224f446f..dcc73620ea66c018ef54eb8421e94e2a4e68cf5b 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -23,6 +23,7 @@ Ecrire doc sur marche à suivre pour intégrer un nouveau développement:
   - dev dans MNH à faire en array-syntax
   - dev dans AROME à faire en boucles do
   - intégration dans PHYEX: en array-syntax avec directives mnh_expand
+  - vérifier les directives mnh_expand à l'aide du script verify_mnh_expand.py
   - les 3 tests suivants doivent donner les mêmes résultats (au bit près) dans chacun des deux modèles:
     - compilation directe sans activer mnh_expand
     - compilation en activant mnh_expand
@@ -52,7 +53,6 @@ Pb identifiés à corriger plus tard:
   - th_r_from_thl_rt appelée partout, il faudrait limiter à OTEST
   - doute sur le codage de MODD_PRECISION
   - appel à abort à travers print_msg non testé
-  - indentation inorrecte dans les blocs mnh_expand
   - sedimentation momentum non branchée (et à trasformer comme sedim_stat)
   - si possible, modifier ice4_sedimentation_split* dans le même esprit que stat
   - il faudrait nettoyer les interfaces pour supprimer les clés passées directement
@@ -80,3 +80,5 @@ Regarder s'il ne serait pas possible/souhaitable de supprimer modd_lunit de PHYE
 Nettoyage des répertoires aux nécessaire
 
 Initialiser dans AROME la variable ldiag_in_run de MODD_DIAG_IN_RUN pour pouvoir phaser le modd
+
+Dans shallow_mf (et toutes les routines appelées en-dessous) il faut remplacer l'utilisation des D%NIB, D%NIE, D%NIT par ce qui sera utilisé dans la turbulence
diff --git a/src/common/micro/mode_ice4_compute_pdf.F90 b/src/common/micro/mode_ice4_compute_pdf.F90
index b1d379c55a85bb75c1d6b7cae2d40e87caca3a10..7ccb88c1274867edbcff476743ebf920ea8dc2cc 100644
--- a/src/common/micro/mode_ice4_compute_pdf.F90
+++ b/src/common/micro/mode_ice4_compute_pdf.F90
@@ -152,8 +152,7 @@ ELSEIF(HSUBG_AUCV_RC=='PDF ') THEN
       PHLC_LCF(:)=0.
       PHLC_HRC(:)=PRCT(:)
       PHLC_LRC(:)=0.
-    ELSEWHERE(PRCT(:)> (ZRCRAUTC(:)-PSIGMA_RC(:)) .AND. &
-            & PRCT(:)<=(ZRCRAUTC(:)+PSIGMA_RC(:))       )
+    ELSEWHERE(PRCT(:)> (ZRCRAUTC(:)-PSIGMA_RC(:)) .AND. PRCT(:)<=(ZRCRAUTC(:)+PSIGMA_RC(:))       )
       PHLC_HCF(:)=(PRCT(:)+PSIGMA_RC(:)-ZRCRAUTC(:))/ &
                   &(2.*PSIGMA_RC(:))
       PHLC_LCF(:)=MAX(0., PCF(:)-PHLC_HCF(:))
diff --git a/src/common/turb/mode_compute_mf_cloud_bigaus.F90 b/src/common/turb/mode_compute_mf_cloud_bigaus.F90
index 767fc9718ce52b3cc10232fee184e0553ebd48af..180532c20cab2bf81bc79b47634932894c890f6f 100644
--- a/src/common/turb/mode_compute_mf_cloud_bigaus.F90
+++ b/src/common/turb/mode_compute_mf_cloud_bigaus.F90
@@ -124,44 +124,44 @@ ZOMEGA_UP_M(:)=0.
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL
   !$mnh_expand_array(JI=D%NIB:D%NIE)
   !Vertical integration over the entire column but only buoyant points are used
-  !ZOMEGA_UP_M(:)=ZOMEGA_UP_M(:) + &
-  !                ZEMF_M(:,JK) * &
-  !                MAX(0.,(ZTHV_UP_M(:,JK)-PTHVM(:,JK))) * &
-  !                (PZZ(:,JK+KKL)-PZZ(:,JK)) / &
-  !                (PTHM(:,JK) * PRHODREF(:,JK))
+  !ZOMEGA_UP_M(D%NIB:D%NIE)=ZOMEGA_UP_M(D%NIB:D%NIE) + &
+  !                ZEMF_M(D%NIB:D%NIE,JK) * &
+  !                MAX(0.,(ZTHV_UP_M(D%NIB:D%NIE,JK)-PTHVM(D%NIB:D%NIE,JK))) * &
+  !                (PZZ(D%NIB:D%NIE,JK+KKL)-PZZ(D%NIB:D%NIE,JK)) / &
+  !                (PTHM(D%NIB:D%NIE,JK) * PRHODREF(D%NIB:D%NIE,JK))
 
   !Vertical integration over the entire column
-  ZOMEGA_UP_M(:)=ZOMEGA_UP_M(:) + &
-                 ZEMF_M(:,JK) * &
-                 (ZTHV_UP_M(:,JK)-PTHVM(:,JK)) * &
-                 (PZZ(:,JK+D%NKL)-PZZ(:,JK)) / &
-                 (PTHM(:,JK) * PRHODREF(:,JK))
+  ZOMEGA_UP_M(D%NIB:D%NIE)=ZOMEGA_UP_M(D%NIB:D%NIE) + &
+                 ZEMF_M(D%NIB:D%NIE,JK) * &
+                 (ZTHV_UP_M(D%NIB:D%NIE,JK)-PTHVM(D%NIB:D%NIE,JK)) * &
+                 (PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK)) / &
+                 (PTHM(D%NIB:D%NIE,JK) * PRHODREF(D%NIB:D%NIE,JK))
   !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 ENDDO
 !$mnh_expand_array(JI=D%NIB:D%NIE)
-ZOMEGA_UP_M(:)=MAX(ZOMEGA_UP_M(:), 1.E-20)
-ZOMEGA_UP_M(:)=(CST%XG*ZOMEGA_UP_M(:))**(1./3.)
+ZOMEGA_UP_M(D%NIB:D%NIE)=MAX(ZOMEGA_UP_M(D%NIB:D%NIE), 1.E-20)
+ZOMEGA_UP_M(D%NIB:D%NIE)=(CST%XG*ZOMEGA_UP_M(D%NIB:D%NIE))**(1./3.)
 !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 !computation of alpha up
 DO JK=D%NKA,D%NKU,D%NKL
   !$mnh_expand_array(JI=D%NIB:D%NIE)
-  ZALPHA_UP_M(:,JK)=ZEMF_M(:,JK)/(PARAMMF%XALPHA_MF*PRHODREF(:,JK)*ZOMEGA_UP_M(:))
-  ZALPHA_UP_M(:,JK)=MAX(0., MIN(ZALPHA_UP_M(:,JK), 1.))
+  ZALPHA_UP_M(D%NIB:D%NIE,JK)=ZEMF_M(D%NIB:D%NIE,JK)/(PARAMMF%XALPHA_MF*PRHODREF(D%NIB:D%NIE,JK)*ZOMEGA_UP_M(D%NIB:D%NIE))
+  ZALPHA_UP_M(D%NIB:D%NIE,JK)=MAX(0., MIN(ZALPHA_UP_M(D%NIB:D%NIE,JK), 1.))
   !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 ENDDO
 
 !computation of sigma of the distribution
 DO JK=D%NKA,D%NKU,D%NKL
   !$mnh_expand_array(JI=D%NIB:D%NIE)
-  ZSIGMF(:,JK)=ZEMF_M(:,JK) * &
-               (ZRT_UP_M(:,JK) - PRTM(:,JK)) * &
-               PDEPTH(:) * ZGRAD_Z_RT(:,JK) / &
-               (PARAMMF%XSIGMA_MF * ZOMEGA_UP_M(:) * PRHODREF(:,JK))
+  ZSIGMF(D%NIB:D%NIE,JK)=ZEMF_M(D%NIB:D%NIE,JK) * &
+               (ZRT_UP_M(D%NIB:D%NIE,JK) - PRTM(D%NIB:D%NIE,JK)) * &
+               PDEPTH(D%NIB:D%NIE) * ZGRAD_Z_RT(D%NIB:D%NIE,JK) / &
+               (PARAMMF%XSIGMA_MF * ZOMEGA_UP_M(D%NIB:D%NIE) * PRHODREF(D%NIB:D%NIE,JK))
   !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 ENDDO
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-ZSIGMF(:,:)=SQRT(MAX(ABS(ZSIGMF(:,:)), 1.E-40))
+ZSIGMF(D%NIB:D%NIE,:)=SQRT(MAX(ABS(ZSIGMF(D%NIB:D%NIE,:)), 1.E-40))
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 !
 !*      2. PDF integration
@@ -170,19 +170,20 @@ ZSIGMF(:,:)=SQRT(MAX(ABS(ZSIGMF(:,:)), 1.E-40))
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 !The mean of the distribution is ZRT_UP
 !Computation of ZA and ZGAM (=efrc(ZA)) coefficient
-ZA(:,:)=(ZRSAT_UP_M(:,:)-ZRT_UP_M(:,:))/(sqrt(2.)*ZSIGMF(:,:))
+ZA(D%NIB:D%NIE,:)=(ZRSAT_UP_M(D%NIB:D%NIE,:)-ZRT_UP_M(D%NIB:D%NIE,:))/(sqrt(2.)*ZSIGMF(D%NIB:D%NIE,:))
 
 !Approximation of erf function
-ZGAM(:,:)=1-SIGN(1., ZA(:,:))*SQRT(1-EXP(-4*ZA(:,:)**2/CST%XPI))
+ZGAM(D%NIB:D%NIE,:)=1-SIGN(1., ZA(D%NIB:D%NIE,:))*SQRT(1-EXP(-4*ZA(D%NIB:D%NIE,:)**2/CST%XPI))
 
 !computation of cloud fraction
-PCF_MF(:,:)=MAX( 0., MIN(1.,0.5*ZGAM(:,:) * ZALPHA_UP_M(:,:)))
+PCF_MF(D%NIB:D%NIE,:)=MAX( 0., MIN(1.,0.5*ZGAM(D%NIB:D%NIE,:) * ZALPHA_UP_M(D%NIB:D%NIE,:)))
 
 !computation of condensate, then PRC and PRI
-ZCOND(:,:)=(EXP(-ZA(:,:)**2)-ZA(:,:)*SQRT(CST%XPI)*ZGAM(:,:))*ZSIGMF(:,:)/SQRT(2.*CST%XPI) * ZALPHA_UP_M(:,:)
-ZCOND(:,:)=MAX(ZCOND(:,:), 0.) !due to approximation of ZGAM value, ZCOND could be slightly negative
-PRC_MF(:,:)=(1.-ZFRAC_ICE_UP_M(:,:)) * ZCOND(:,:)
-PRI_MF(:,:)=(   ZFRAC_ICE_UP_M(:,:)) * ZCOND(:,:)
+ZCOND(D%NIB:D%NIE,:)=(EXP(-ZA(D%NIB:D%NIE,:)**2)-ZA(D%NIB:D%NIE,:)*SQRT(CST%XPI)*ZGAM(D%NIB:D%NIE,:))* &
+                    &ZSIGMF(D%NIB:D%NIE,:)/SQRT(2.*CST%XPI) * ZALPHA_UP_M(D%NIB:D%NIE,:)
+ZCOND(D%NIB:D%NIE,:)=MAX(ZCOND(D%NIB:D%NIE,:), 0.) !due to approximation of ZGAM value, ZCOND could be slightly negative
+PRC_MF(D%NIB:D%NIE,:)=(1.-ZFRAC_ICE_UP_M(D%NIB:D%NIE,:)) * ZCOND(D%NIB:D%NIE,:)
+PRI_MF(D%NIB:D%NIE,:)=(   ZFRAC_ICE_UP_M(D%NIB:D%NIE,:)) * ZCOND(D%NIB:D%NIE,:)
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 !
 IF (LHOOK) CALL DR_HOOK('COMPUTE_MF_CLOUD_BIGAUS',1,ZHOOK_HANDLE)
diff --git a/src/common/turb/mode_compute_mf_cloud_stat.F90 b/src/common/turb/mode_compute_mf_cloud_stat.F90
index 4b1ca91a2f14206fab9254d662ea447622bea4a9..28032ab6218ed7a02c471b81b9c68bf7034e9698 100644
--- a/src/common/turb/mode_compute_mf_cloud_stat.F90
+++ b/src/common/turb/mode_compute_mf_cloud_stat.F90
@@ -117,15 +117,16 @@ IF (KRRL > 0)  THEN
     CALL MZM_MF(D, PTHLM(:,:), ZFLXZ(:,:))
     CALL GZ_M_W_MF(D, PTHLM(:,:), PDZZ(:,:), ZWK(:,:))
     !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-    ZFLXZ(:,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(:,:)*(PTHL_UP(:,:)-ZFLXZ(:,:)) * ZWK(:,:)
+    ZFLXZ(D%NIB:D%NIE,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(D%NIB:D%NIE,:)* &
+                         & (PTHL_UP(D%NIB:D%NIE,:)-ZFLXZ(D%NIB:D%NIE,:)) * ZWK(D%NIB:D%NIE,:)
     !
     !   Avoid negative values
-    ZFLXZ(:,:) = MAX(0.,ZFLXZ(:,:))
+    ZFLXZ(D%NIB:D%NIE,:) = MAX(0.,ZFLXZ(D%NIB:D%NIE,:))
     !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
     CALL MZF_MF(D, ZFLXZ(:,:), PSIGMF(:,:))
     !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-    PSIGMF(:,:) = PSIGMF(:,:) * ZATHETA(:,:)**2
+    PSIGMF(D%NIB:D%NIE,:) = PSIGMF(D%NIB:D%NIE,:) * ZATHETA(D%NIB:D%NIE,:)**2
     !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 !
@@ -138,21 +139,22 @@ IF (KRRL > 0)  THEN
     CALL MZM_MF(D, PRTM(:,:), ZFLXZ(:,:))
     CALL GZ_M_W_MF(D, PRTM(:,:), PDZZ(:,:), ZWK(:,:))
     !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-    ZFLXZ(:,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(:,:)*(PRT_UP(:,:)-ZFLXZ(:,:)) * ZWK(:,:)
+    ZFLXZ(D%NIB:D%NIE,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(D%NIB:D%NIE,:)* &
+                         & (PRT_UP(D%NIB:D%NIE,:)-ZFLXZ(D%NIB:D%NIE,:)) * ZWK(D%NIB:D%NIE,:)
     !
     !   Avoid negative values
-    ZFLXZ(:,:) = MAX(0.,ZFLXZ(:,:))
+    ZFLXZ(D%NIB:D%NIE,:) = MAX(0.,ZFLXZ(D%NIB:D%NIE,:))
     !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
     CALL MZF_MF(D, ZFLXZ(:,:), ZWK(:,:))
     !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-    PSIGMF(:,:) = PSIGMF(:,:) + ZAMOIST(:,:) **2 * ZWK(:,:)
+    PSIGMF(D%NIB:D%NIE,:) = PSIGMF(D%NIB:D%NIE,:) + ZAMOIST(D%NIB:D%NIE,:) **2 * ZWK(D%NIB:D%NIE,:)
     !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 !
 !        1.3  Vertical part of Sigma_s
 !
     !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-    PSIGMF(:,:) =  SQRT( MAX (PSIGMF(:,:) , 0.) )
+    PSIGMF(D%NIB:D%NIE,:) =  SQRT( MAX (PSIGMF(D%NIB:D%NIE,:) , 0.) )
     !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ELSE
   PSIGMF(:,:) = 0.
diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index cf8b664ecc9ec8ce94a4241c71eaa80c033642cb..6e3a533dc072a25b2bf35e0e91a744e24d64e745 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -349,7 +349,7 @@ IF (OENTR_DETR) THEN
   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(:,:))
+             PBUF=ZBUF(:,:), KB=D%NIB, KE=D%NIE)
 
   !$mnh_expand_array(JI=D%NIB:D%NIE)
   ! compute updraft thevav and buoyancy term at KKB level
@@ -570,7 +570,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
     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(:,:))
+            PBUF=ZBUF(:,:), KB=D%NIB, KE=D%NIE)
     !$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)
@@ -935,7 +935,7 @@ 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)
+             PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
 !$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)
@@ -1015,7 +1015,7 @@ 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)
+             PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
 !$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))
 
@@ -1027,7 +1027,7 @@ 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)
+             PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
 !$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)
diff --git a/src/common/turb/mode_compute_updraft_raha.F90 b/src/common/turb/mode_compute_updraft_raha.F90
index a17f1deade8c192ef983e472001547a2c2497504..bb1dc222569a9a1e3d2c3154c15ee5f00c17ea88 100644
--- a/src/common/turb/mode_compute_updraft_raha.F90
+++ b/src/common/turb/mode_compute_updraft_raha.F90
@@ -119,7 +119,7 @@ REAL, DIMENSION(D%NIT,D%NKT),   INTENT(INOUT)::  PEMF,PDETR,PENTR ! Mass_flux,
                                                           ! detrainment,entrainment
 REAL, DIMENSION(D%NIT,D%NKT),   INTENT(INOUT) :: PBUO_INTEG       ! Integrated Buoyancy 
 INTEGER, DIMENSION(D%NIT),  INTENT(INOUT)::  KKLCL,KKETL,KKCTL! LCL, ETL, CTL                                     
-REAL, DIMENSION(:),     INTENT(OUT)   :: PDEPTH           ! Deepness of cloud
+REAL, DIMENSION(D%NIT),     INTENT(OUT)   :: PDEPTH           ! Deepness of cloud
 !                       1.2  Declaration of local variables
 !
 !
@@ -236,9 +236,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)
+!$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)
 
 ! Initialisation of environment variables at t-dt
 
@@ -251,31 +251,33 @@ CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 
 !DO JSV=1,ISV 
 ! IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
-!   ZSVM_F(:,KKB:IKU,JSV) = 0.5*(PSVM(:,KKB:IKU,JSV)+PSVM(:,1:IKU-1,JSV))
-!   ZSVM_F(:,1,JSV)       = ZSVM_F(:,KKB,JSV) 
+!   ZSVM_F(D%NIB:D%NIE,KKB:IKU,JSV) = 0.5*(PSVM(D%NIB:D%NIE,KKB:IKU,JSV)+PSVM(D%NIB:D%NIE,1:IKU-1,JSV))
+!   ZSVM_F(D%NIB:D%NIE,1,JSV)       = ZSVM_F(D%NIB:D%NIE,KKB,JSV) 
 !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)
+!$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)
 PSV_UP(:,:,:)=0.
 !IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) then
-!    PSV_UP(:,:,:)=ZSVM_F(:,:,:)
+!    PSV_UP(D%NIB:D%NIE,:,:)=ZSVM_F(D%NIB:D%NIE,:,:)
 !ENDIF
 
 ! 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)) 
+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)) 
 
-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(:))
+ZQT_UP(D%NIB:D%NIE) = PRT_UP(D%NIB:D%NIE,D%NKB)/(1.+PRT_UP(D%NIB:D%NIE,D%NKB))
+ZTHS_UP(D%NIB:D%NIE,D%NKB)=PTHL_UP(D%NIB:D%NIE,D%NKB)*(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(D%NIB:D%NIE))
 !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 CALL MZM_MF(D, PTHM (:,:), ZTHM_F(:,:))
@@ -283,45 +285,46 @@ 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)
+!$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(:,:)
-PRV_UP(:,:) = ZRVM_F(:,:)
-!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+PTHV_UP(D%NIB:D%NIE,:)= ZTHVM_F(D%NIB:D%NIE,:)
+PRV_UP(D%NIB:D%NIE,:) = ZRVM_F(D%NIB:D%NIE,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
 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)  
+ZW_UP2(D%NIB:D%NIE,D%NKB) = MAX(0.0001,(1./6.)*ZTKEM_F(D%NIB:D%NIE,D%NKB))
+GTEST(D%NIB:D%NIE) = (ZW_UP2(D%NIB:D%NIE,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.
+PRC_UP(D%NIB:D%NIE,D%NKB)=0.
+PRI_UP(D%NIB:D%NIE,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), &
+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)
+             PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
 
 !$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)
 
 !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                                                            
-!$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)
+!$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)
 
 !  Definition de l'alimentation au sens de la fermeture de Hourdin et al
 
@@ -331,14 +334,14 @@ 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
+  ZZDZ(D%NIB:D%NIE,JK)   = MAX(ZEPS,PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))       ! <== Delta Z between two flux level
+  ZZZ(D%NIB:D%NIE,JK)    = MAX(0.,0.5*(PZZ(D%NIB:D%NIE,JK+D%NKL)+PZZ(D%NIB:D%NIE,JK)) )  ! <== Hight of mass levels
+  ZDTHETASDZ(D%NIB:D%NIE,JK) = (ZTHVM_F(D%NIB:D%NIE,JK)-ZTHVM_F(D%NIB:D%NIE,JK+D%NKL))   ! <== Delta theta_v
   
-  WHERE ((ZTHVM_F(:,JK+D%NKL)<ZTHVM_F(:,JK)) .AND. (ZTHVM_F(:,D%NKB)>=ZTHVM_F(:,JK)))
-     ZALIM_STAR(:,JK)  = SQRT(ZZZ(:,JK))*ZDTHETASDZ(:,JK)/ZZDZ(:,JK)
-     ZALIM_STAR_TOT(:) = ZALIM_STAR_TOT(:)+ZALIM_STAR(:,JK)*ZZDZ(:,JK)
-     IALIM(:)          = JK
+  WHERE ((ZTHVM_F(D%NIB:D%NIE,JK+D%NKL)<ZTHVM_F(D%NIB:D%NIE,JK)) .AND. (ZTHVM_F(D%NIB:D%NIE,D%NKB)>=ZTHVM_F(D%NIB:D%NIE,JK)))
+     ZALIM_STAR(D%NIB:D%NIE,JK)  = SQRT(ZZZ(D%NIB:D%NIE,JK))*ZDTHETASDZ(D%NIB:D%NIE,JK)/ZZDZ(D%NIB:D%NIE,JK)
+     ZALIM_STAR_TOT(D%NIB:D%NIE) = ZALIM_STAR_TOT(D%NIB:D%NIE)+ZALIM_STAR(D%NIB:D%NIE,JK)*ZZDZ(D%NIB:D%NIE,JK)
+     IALIM(D%NIB:D%NIE)          = JK
   ENDWHERE
   !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO
@@ -346,8 +349,8 @@ ENDDO
 ! Normalization of ZALIM_STAR
 DO JK=D%NKB,D%NKE-D%NKL,D%NKL   !  Vertical loop   
    !$mnh_expand_where(JI=D%NIB:D%NIE)
-   WHERE (ZALIM_STAR_TOT(:) > ZEPS)
-      ZALIM_STAR(:,JK)  = ZALIM_STAR(:,JK)/ZALIM_STAR_TOT(:)
+   WHERE (ZALIM_STAR_TOT(D%NIB:D%NIE) > ZEPS)
+      ZALIM_STAR(D%NIB:D%NIE,JK)  = ZALIM_STAR(D%NIB:D%NIE,JK)/ZALIM_STAR_TOT(D%NIB:D%NIE)
    ENDWHERE
    !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO
@@ -377,164 +380,170 @@ 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
+  ! IF the updraft top is reached for all column, stop the loop on levels
 
-!  ITEST=COUNT(GTEST)
-!  IF (ITEST==0) CYCLE
+  !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(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
 
-! 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.
-    ENDWHERE
-
-! COMPUTE PENTR and PDETR at mass level JK
+  ! COMPUTE PENTR and PDETR at mass level JK
 
     
-!  Buoyancy is computed on "flux" levels where updraft variables are known
+  !  Buoyancy is computed on "flux" levels where updraft variables are known
 
   ! Compute theta_v of updraft at flux level JK    
     
-    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))
-    PBUO_INTEG(:,JK) = ZBUO(:,JK)*(PZZ(:,JK+D%NKL)-PZZ(:,JK))
-    
-    ZDZ(:)   = MAX(ZEPS,PZZ(:,JK+D%NKL)-PZZ(:,JK))
-    ZTEST(:) = PARAMMF%XA1*ZBUO(:,JK) -  PARAMMF%XB*ZW_UP2(:,JK)
+  ZRC_UP(D%NIB:D%NIE)  = PRC_UP(D%NIB:D%NIE,JK)
+  ZRI_UP(D%NIB:D%NIE)  = PRI_UP(D%NIB:D%NIE,JK) ! guess 
+  ZRV_UP(D%NIB:D%NIE)  = PRV_UP(D%NIB:D%NIE,JK)
+  ZBUO(D%NIB:D%NIE,JK) = ZG_O_THVREF(D%NIB:D%NIE,JK)*(PTHV_UP(D%NIB:D%NIE,JK) - ZTHVM_F(D%NIB:D%NIE,JK))
+  PBUO_INTEG(D%NIB:D%NIE,JK) = ZBUO(D%NIB:D%NIE,JK)*(PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))
+  
+  ZDZ(D%NIB:D%NIE)   = MAX(ZEPS,PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))
+  ZTEST(D%NIB:D%NIE) = PARAMMF%XA1*ZBUO(D%NIB:D%NIE,JK) -  PARAMMF%XB*ZW_UP2(D%NIB:D%NIE,JK)
 
-    ZCOE(:)      = ZDZ(:)
-    WHERE (ZTEST(:)>0.)
-      ZCOE(:)    = ZDZ(:)/(1.+ PARAMMF%XBETA1)
-    ENDWHERE   
+  ZCOE(D%NIB:D%NIE)      = ZDZ(D%NIB:D%NIE)
+  WHERE (ZTEST(D%NIB:D%NIE)>0.)
+    ZCOE(D%NIB:D%NIE)    = ZDZ(D%NIB:D%NIE)/(1.+ PARAMMF%XBETA1)
+  ENDWHERE   
 
-!  Calcul de la vitesse
+  !  Calcul de la vitesse
 
-    ZWCOE(:)         = (1.-PARAMMF%XB*ZCOE(:))/(1.+PARAMMF%XB*ZCOE(:))
-    ZBUCOE(:)        =  2.*ZCOE(:)/(1.+PARAMMF%XB*ZCOE(:))
-    
-    ZW_UP2(:,JK+D%NKL) = MAX(ZEPS,ZW_UP2(:,JK)*ZWCOE(:) + PARAMMF%XA1*ZBUO(:,JK)*ZBUCOE(:) )  
-    ZW_MAX(:) = MAX(ZW_MAX(:), SQRT(ZW_UP2(:,JK+D%NKL)))
-    ZWUP_MEAN(:)     = MAX(ZEPS,0.5*(ZW_UP2(:,JK+D%NKL)+ZW_UP2(:,JK)))
- 
-!  Entrainement et detrainement
+  ZWCOE(D%NIB:D%NIE)         = (1.-PARAMMF%XB*ZCOE(D%NIB:D%NIE))/(1.+PARAMMF%XB*ZCOE(D%NIB:D%NIE))
+  ZBUCOE(D%NIB:D%NIE)        =  2.*ZCOE(D%NIB:D%NIE)/(1.+PARAMMF%XB*ZCOE(D%NIB:D%NIE))
 
-   PENTR(:,JK)  = MAX(0.,(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*(PARAMMF%XA1*ZBUO(:,JK)/ZWUP_MEAN(:)-PARAMMF%XB))
-   
-   ZDETR_BUO(:) = MAX(0., -(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*PARAMMF%XA1*ZBUO(:,JK)/ZWUP_MEAN(:))
-   ZDETR_RT(:)  = PARAMMF%XC*SQRT(MAX(0.,(PRT_UP(:,JK) - ZRTM_F(:,JK))) / MAX(ZEPS,ZRTM_F(:,JK)) / ZWUP_MEAN(:))
-   PDETR(:,JK)  = ZDETR_RT(:)+ZDETR_BUO(:)
+  ZW_UP2(D%NIB:D%NIE,JK+D%NKL) = MAX(ZEPS,ZW_UP2(D%NIB:D%NIE,JK)*ZWCOE(D%NIB:D%NIE) + &
+                                         &PARAMMF%XA1*ZBUO(D%NIB:D%NIE,JK)*ZBUCOE(D%NIB:D%NIE))
+  ZW_MAX(D%NIB:D%NIE) = MAX(ZW_MAX(D%NIB:D%NIE), SQRT(ZW_UP2(D%NIB:D%NIE,JK+D%NKL)))
+  ZWUP_MEAN(D%NIB:D%NIE)     = MAX(ZEPS,0.5*(ZW_UP2(D%NIB:D%NIE,JK+D%NKL)+ZW_UP2(D%NIB:D%NIE,JK)))
+
+  !  Entrainement et detrainement
+
+  PENTR(D%NIB:D%NIE,JK)  = MAX(0.,(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))* &
+                                 &(PARAMMF%XA1*ZBUO(D%NIB:D%NIE,JK)/ZWUP_MEAN(D%NIB:D%NIE)-PARAMMF%XB))
+  
+  ZDETR_BUO(D%NIB:D%NIE) = MAX(0., -(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*PARAMMF%XA1*ZBUO(D%NIB:D%NIE,JK)/ &
+                                   &ZWUP_MEAN(D%NIB:D%NIE))
+  ZDETR_RT(D%NIB:D%NIE)  = PARAMMF%XC*SQRT(MAX(0.,(PRT_UP(D%NIB:D%NIE,JK) - ZRTM_F(D%NIB:D%NIE,JK))) / &
+                                          &MAX(ZEPS,ZRTM_F(D%NIB:D%NIE,JK)) / ZWUP_MEAN(D%NIB:D%NIE))
+  PDETR(D%NIB:D%NIE,JK)  = ZDETR_RT(D%NIB:D%NIE)+ZDETR_BUO(D%NIB:D%NIE)
 
    
-! If the updraft did not stop, compute cons updraft characteritics at jk+1
-  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) !&                   
+  ! If the updraft did not stop, compute cons updraft characteritics at jk+1
+  WHERE(GTEST(D%NIB:D%NIE))     
+    ZZTOP(D%NIB:D%NIE) = MAX(ZZTOP(D%NIB:D%NIE),PZZ(D%NIB:D%NIE,JK+D%NKL))
+    ZMIX2(D%NIB:D%NIE) = (PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*PENTR(D%NIB:D%NIE,JK) !&
+    ZMIX3(D%NIB:D%NIE) = (PZZ(D%NIB:D%NIE,JK+D%NKL)-PZZ(D%NIB:D%NIE,JK))*PDETR(D%NIB:D%NIE,JK) !&                   
            
-    ZQTM(:) = PRTM(:,JK)/(1.+PRTM(:,JK))            
-    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(:))  &
-                          /(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(:))                      
+    ZQTM(D%NIB:D%NIE) = PRTM(D%NIB:D%NIE,JK)/(1.+PRTM(D%NIB:D%NIE,JK))            
+    ZTHSM(D%NIB:D%NIE,JK) = PTHLM(D%NIB:D%NIE,JK)*(1.+PARAMMF%XLAMBDA_MF*ZQTM(D%NIB:D%NIE))
+    ZTHS_UP(D%NIB:D%NIE,JK+D%NKL)=(ZTHS_UP(D%NIB:D%NIE,JK)*(1.-0.5*ZMIX2(D%NIB:D%NIE)) + ZTHSM(D%NIB:D%NIE,JK)*ZMIX2(D%NIB:D%NIE))&
+                          /(1.+0.5*ZMIX2(D%NIB:D%NIE))   
+    PRT_UP(D%NIB:D%NIE,JK+D%NKL)=(PRT_UP(D%NIB:D%NIE,JK)*(1.-0.5*ZMIX2(D%NIB:D%NIE)) + PRTM(D%NIB:D%NIE,JK)*ZMIX2(D%NIB:D%NIE))  &
+                          /(1.+0.5*ZMIX2(D%NIB:D%NIE))
+    ZQT_UP(D%NIB:D%NIE) = PRT_UP(D%NIB:D%NIE,JK+D%NKL)/(1.+PRT_UP(D%NIB:D%NIE,JK+D%NKL))
+    PTHL_UP(D%NIB:D%NIE,JK+D%NKL)=ZTHS_UP(D%NIB:D%NIE,JK+D%NKL)/(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(D%NIB:D%NIE))                      
   ENDWHERE
   
 
   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(:))
+      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
     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(:))
+      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
 
     ENDIF
   ENDIF
-!  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(:)) + &
-!                        PSVM(:,JK,JSV)*ZMIX2(:))  /(1+0.5*ZMIX2(:))
-!      ENDWHERE                        
-!  ENDDO  
+  !DO JSV=1,ISV 
+  !   IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
+  !    WHERE(GTEST(D%NIB:D%NIE)) 
+  !         PSV_UP(D%NIB:D%NIE,JK+KKL,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                        
+  !ENDDO  
   
   
-! 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
-  ZRV_UP(:)=PRV_UP(:,JK)
+  ! Compute non cons. var. at level JK+KKL
+  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
+  ZRV_UP(D%NIB:D%NIE)=PRV_UP(D%NIB:D%NIE,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)
+          PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
   !$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)
-    PRC_UP(:,JK+D%NKL)=MIN(0.5E-3,ZRC_UP(:))  ! On ne peut depasser 0.5 g/kg (autoconversion donc elimination !)
-    PTHL_UP(:,JK+D%NKL) = PTHL_UP(:,JK+D%NKL)+ZLVOCPEXN(:)*(ZRC_UP(:)-PRC_UP(:,JK+D%NKL))
-    PRV_UP(:,JK+D%NKL)=ZRV_UP(:)
-    PRI_UP(:,JK+D%NKL)=ZRI_UP(:)
-    PRT_UP(:,JK+D%NKL)  = PRC_UP(:,JK+D%NKL) + PRV_UP(:,JK+D%NKL)
-    PRSAT_UP(:,JK+D%NKL) = ZRSATW(:)*(1-PFRAC_ICE_UP(:,JK+D%NKL)) + ZRSATI(:)*PFRAC_ICE_UP(:,JK+D%NKL)
+  WHERE(GTEST(D%NIB:D%NIE))
+    ZT_UP(D%NIB:D%NIE) = ZTH_UP(D%NIB:D%NIE,JK+D%NKL)*PEXNM(D%NIB:D%NIE,JK+D%NKL)
+    ZCP(D%NIB:D%NIE) = CST%XCPD + CST%XCL * ZRC_UP(D%NIB:D%NIE)
+    ZLVOCPEXN(D%NIB:D%NIE)=(CST%XLVTT + (CST%XCPV-CST%XCL) *  (ZT_UP(D%NIB:D%NIE)-CST%XTT) ) / &
+                          &ZCP(D%NIB:D%NIE) / PEXNM(D%NIB:D%NIE,JK+D%NKL)
+    PRC_UP(D%NIB:D%NIE,JK+D%NKL)=MIN(0.5E-3,ZRC_UP(D%NIB:D%NIE))  ! On ne peut depasser 0.5 g/kg (autoconversion donc elimination !)
+    PTHL_UP(D%NIB:D%NIE,JK+D%NKL) = PTHL_UP(D%NIB:D%NIE,JK+D%NKL)+ &
+                                  & ZLVOCPEXN(D%NIB:D%NIE)*(ZRC_UP(D%NIB:D%NIE)-PRC_UP(D%NIB:D%NIE,JK+D%NKL))
+    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)
+    PRT_UP(D%NIB:D%NIE,JK+D%NKL)  = PRC_UP(D%NIB:D%NIE,JK+D%NKL) + PRV_UP(D%NIB:D%NIE,JK+D%NKL)
+    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
   
 
-! Compute the updraft theta_v, buoyancy and w**2 for level JK+1   
-  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))
+  ! Compute the updraft theta_v, buoyancy and w**2 for level JK+1   
+  WHERE(GTEST(D%NIB:D%NIE))
+    !PTHV_UP(D%NIB:D%NIE,JK+KKL) = ZTH_UP(D%NIB:D%NIE,JK+KKL)*((1+ZRVORD*PRV_UP(D%NIB:D%NIE,JK+KKL))/(1+PRT_UP(D%NIB:D%NIE,JK+KKL)))
+    PTHV_UP(D%NIB:D%NIE,JK+D%NKL) = ZTH_UP(D%NIB:D%NIE,JK+D%NKL)* &
+                                  & (1.+0.608*PRV_UP(D%NIB:D%NIE,JK+D%NKL) - PRC_UP(D%NIB:D%NIE,JK+D%NKL))
   ENDWHERE
 
 
-! Test if the updraft has reach the ETL
-  GTESTETL(:)=.FALSE.
-  WHERE (GTEST(:).AND.(PBUO_INTEG(:,JK)<=0.))
-      KKETL(:) = JK+D%NKL
-      GTESTETL(:)=.TRUE.
+  ! Test if the updraft has reach the ETL
+  GTESTETL(D%NIB:D%NIE)=.FALSE.
+  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.
   ENDWHERE
 
-! Test is we have reached the top of the updraft
-
-  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)
-      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
+  ! 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)<=ZEPS)))
+    ZW_UP2(D%NIB:D%NIE,JK+D%NKL)=ZEPS
+    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
   !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO 
@@ -544,56 +553,60 @@ ENDDO
 
 
 !$mnh_expand_array(JI=D%NIB:D%NIE)
-ZZTOP(:) = MAX(ZZTOP(:),ZEPS)
+ZZTOP(D%NIB:D%NIE) = MAX(ZZTOP(D%NIB:D%NIE),ZEPS)
 !$mnh_end_expand_array(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)
-   WHERE(JK<=IALIM(:))
-     ZALIM_STAR_TOT(:) = ZALIM_STAR_TOT(:) + ZALIM_STAR(:,JK)*ZALIM_STAR(:,JK)*ZZDZ(:,JK)/PRHODREF(:,JK)
+   WHERE(JK<=IALIM(D%NIB:D%NIE))
+     ZALIM_STAR_TOT(D%NIB:D%NIE) = ZALIM_STAR_TOT(D%NIB:D%NIE) + ZALIM_STAR(D%NIB:D%NIE,JK)**2* &
+                                                               & ZZDZ(D%NIB:D%NIE,JK)/PRHODREF(D%NIB:D%NIE,JK)
    ENDWHERE
   !$mnh_end_expand_where(JI=D%NIB:D%NIE)
 ENDDO   
 
 !$mnh_expand_where(JI=D%NIB:D%NIE)
-WHERE (ZALIM_STAR_TOT(:)*ZZTOP(:) > ZEPS)
-   ZPHI(:) =  ZW_MAX(:)/(PARAMMF%XR*ZZTOP(:)*ZALIM_STAR_TOT(:))
+WHERE (ZALIM_STAR_TOT(D%NIB:D%NIE)*ZZTOP(D%NIB:D%NIE) > ZEPS)
+ ZPHI(D%NIB:D%NIE) =  ZW_MAX(D%NIB:D%NIE)/(PARAMMF%XR*ZZTOP(D%NIB:D%NIE)*ZALIM_STAR_TOT(D%NIB:D%NIE))
 ENDWHERE
 
-GTEST(:) = .TRUE.
-PEMF(:,D%NKB+D%NKL) = ZPHI(:)*ZZDZ(:,D%NKB)*ZALIM_STAR(:,D%NKB)
+GTEST(D%NIB:D%NIE) = .TRUE.
+PEMF(D%NIB:D%NIE,D%NKB+D%NKL) = ZPHI(D%NIB:D%NIE)*ZZDZ(D%NIB:D%NIE,D%NKB)*ZALIM_STAR(D%NIB:D%NIE,D%NKB)
 ! Updraft fraction must be smaller than XFRAC_UP_MAX
-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))
+PFRAC_UP(D%NIB:D%NIE,D%NKB+D%NKL)=PEMF(D%NIB:D%NIE,D%NKB+D%NKL)/ &
+                                 &(SQRT(ZW_UP2(D%NIB:D%NIE,D%NKB+D%NKL))*ZRHO_F(D%NIB:D%NIE,D%NKB+D%NKL))
+PFRAC_UP(D%NIB:D%NIE,D%NKB+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX,PFRAC_UP(D%NIB:D%NIE,D%NKB+D%NKL))
+PEMF(D%NIB:D%NIE,D%NKB+D%NKL) = ZRHO_F(D%NIB:D%NIE,D%NKB+D%NKL)*PFRAC_UP(D%NIB:D%NIE,D%NKB+D%NKL)* &
+                              & SQRT(ZW_UP2(D%NIB:D%NIE,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(D%NIB:D%NIE) = (ZW_UP2(D%NIB:D%NIE,JK) > ZEPS)  
 
-  WHERE (GTEST(:))
-    WHERE(JK<IALIM(:))
-       PEMF(:,JK+D%NKL) = MAX(0.,PEMF(:,JK) + ZPHI(:)*ZZDZ(:,JK)*(PENTR(:,JK) - PDETR(:,JK)))
+  WHERE (GTEST(D%NIB:D%NIE))
+    WHERE(JK<IALIM(D%NIB:D%NIE))
+      PEMF(D%NIB:D%NIE,JK+D%NKL) = MAX(0.,PEMF(D%NIB:D%NIE,JK) + ZPHI(D%NIB:D%NIE)*ZZDZ(D%NIB:D%NIE,JK)* &
+                                                                & (PENTR(D%NIB:D%NIE,JK) - PDETR(D%NIB:D%NIE,JK)))
     ELSEWHERE
-       ZMIX1(:)=ZZDZ(:,JK)*(PENTR(:,JK)-PDETR(:,JK))
-       PEMF(:,JK+D%NKL)=PEMF(:,JK)*EXP(ZMIX1(:))
+      ZMIX1(D%NIB:D%NIE)=ZZDZ(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(ZMIX1(D%NIB:D%NIE))
     ENDWHERE
 
 ! Updraft fraction must be smaller than XFRAC_UP_MAX
-    PFRAC_UP(:,JK+D%NKL)=PEMF(:,JK+D%NKL)/(SQRT(ZW_UP2(:,JK+D%NKL))*ZRHO_F(:,JK+D%NKL))
-    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))
+    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))
+    PFRAC_UP(D%NIB:D%NIE,JK+D%NKL)=MIN(PARAMMF%XFRAC_UP_MAX,PFRAC_UP(D%NIB:D%NIE,JK+D%NKL))
+    PEMF(D%NIB:D%NIE,JK+D%NKL) = ZRHO_F(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))
   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,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)
 !$mnh_expand_array(JI=D%NIB:D%NIE)
-PEMF(:,D%NKB) =0.
+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.
@@ -603,25 +616,25 @@ PEMF(:,D%NKB) =0.
 ! 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)) )
+  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) )
+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=D%NKTB,D%NKTE
+DO JK=1,D%NKT
   !$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.)
+  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=D%NKTB:D%NKTE)
-WHERE (GWORK2(:,:)) 
-   PEMF(:,:)     = PEMF(:,:)     * ZCOEF(:,:)
-   PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
+!$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=D%NKTB:D%NKTE)
+!$mnh_end_expand_where(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
 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 9e97f76bbae74c6f10a46cd438daece65a8cb100..c3fcc7fc5a5ae05b519be3ec2b443de64c3c28ab 100644
--- a/src/common/turb/mode_compute_updraft_rhcj10.F90
+++ b/src/common/turb/mode_compute_updraft_rhcj10.F90
@@ -242,9 +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)
+!$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)
 
 ! Initialisation of environment variables at t-dt
 
@@ -259,23 +259,23 @@ CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 !DO JSV=1,ISV
 !  IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
 ! *** SR merge AROME/Meso-nh: following two lines come from the AROME version
-!   ZSVM_F(:,KKB:IKU,JSV) = 0.5*(PSVM(:,KKB:IKU,JSV)+PSVM(:,1:IKU-1,JSV))
-!   ZSVM_F(:,1,JSV)       = ZSVM_F(:,KKB,JSV) 
+!   ZSVM_F(D%NIB:D%NIE,KKB:IKU,JSV) = 0.5*(PSVM(D%NIB:D%NIE,KKB:IKU,JSV)+PSVM(D%NIB:D%NIE,1:IKU-1,JSV))
+!   ZSVM_F(D%NIB:D%NIE,1,JSV)       = ZSVM_F(D%NIB:D%NIE,KKB,JSV) 
 ! *** the following single line comes from the Meso-NH version
-!  ZSVM_F(:,:,JSV) = MZM_MF(KKA,KKU,KKL,PSVM(:,:,JSV))
+!  ZSVM_F(D%NIB:D%NIE,:,JSV) = MZM_MF(KKA,KKU,KKL,PSVM(D%NIB:D%NIE,:,JSV))
 !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.
+!$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)
+PSV_UP(D%NIB:D%NIE,:,:)=0.
 ! This updraft is not yet ready to use scalar variables
 !IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) then
-!    PSV_UP(:,:,:)=ZSVM_F(:,:,:)
+!    PSV_UP(D%NIB:D%NIE,:,:)=ZSVM_F(D%NIB:D%NIE,:,:)
 !ENDIF
 
 ! Computation or initialisation of updraft characteristics at the KKB level
@@ -296,34 +296,34 @@ CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F(:,:))
 CALL MZM_MF(D, PRVM(:,:), ZRVM_F(:,:))
 
 ! thetav at mass and flux levels 
-DO JK=D%NKTB,D%NKTE
+DO JK=1,D%NKT
   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(:,:)
-!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+PTHV_UP(D%NIB:D%NIE,:)= ZTHVM_F(D%NIB:D%NIE,:)
+PRV_UP(D%NIB:D%NIE,:)= ZRVM_F(D%NIB:D%NIE,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 
 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))
+!ZW_UP2(D%NIB:D%NIE,KKB) = MAX(0.0001,(3./6.)*ZTKEM_F(D%NIB:D%NIE,KKB))
+ZW_UP2(D%NIB:D%NIE,D%NKB) = MAX(0.0001,(2./3.)*ZTKEM_F(D%NIB:D%NIE,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.
+PRC_UP(D%NIB:D%NIE,D%NKB)=0.
+PRI_UP(D%NIB:D%NIE,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)
+             PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
 
 DO JI=D%NIB,D%NIE
   ! compute updraft thevav and buoyancy term at KKB level             
@@ -335,9 +335,9 @@ 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                                                            
 
-!$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)
+!$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)
 
 ! Calcul de la fermeture de Julien Pergaut comme limite max de PHY
 
@@ -350,7 +350,7 @@ ENDDO
 ! compute L_up
 GLMIX=.TRUE.
 !$mnh_expand_array(JI=D%NIB:D%NIE)
-ZTKEM_F(:,D%NKB)=0.
+ZTKEM_F(D%NIB:D%NIE,D%NKB)=0.
 !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 !
 IF(TURB%CTURBLEN=='RM17') THEN
@@ -358,17 +358,17 @@ IF(TURB%CTURBLEN=='RM17') THEN
   CALL MZF_MF(D, ZWK, ZDUDZ)
   CALL GZ_M_W_MF(D, PVM, PDZZ, ZWK)
   CALL MZF_MF(D, ZWK, 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)
+  !$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
+  ZSHEAR(D%NIB:D%NIE,:) = 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)
+ZLUP(D%NIB:D%NIE)=MAX(ZLUP(D%NIB:D%NIE),1.E-10)
 !$mnh_end_expand_array(JI=D%NIB:D%NIE)
 
 DO JI=D%NIB,D%NIE
@@ -435,14 +435,14 @@ 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)
+    ZRC_UP(D%NIB:D%NIE)   =PRC_UP(D%NIB:D%NIE,JK) ! guess
+    ZRI_UP(D%NIB:D%NIE)   =PRI_UP(D%NIB:D%NIE,JK) ! guess 
+    ZRV_UP(D%NIB:D%NIE)   =PRV_UP(D%NIB:D%NIE,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)
+               PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
     
   DO JI=D%NIB,D%NIE
     IF (GTEST(JI)) THEN
@@ -530,22 +530,22 @@ 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(:)) + &
-!                        PSVM(:,JK,JSV)*ZMIX2(:))  /(1+0.5*ZMIX2(:))
+!           PSV_UP(D%NIB:D%NIE,JK+KKL,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                        
 !  ENDDO  
   
   
   ! 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)
+  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
+  ZRV_UP(D%NIB:D%NIE)=PRV_UP(D%NIB:D%NIE,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)
+          PBUF=ZBUF, KB=D%NIB, KE=D%NIE)
 
   DO JI=D%NIB,D%NIE
     IF(GTEST(JI)) THEN
@@ -561,7 +561,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
       PRSAT_UP(JI,JK+D%NKL) = ZRSATW(JI)*(1-PFRAC_ICE_UP(JI,JK+D%NKL)) + ZRSATI(JI)*PFRAC_ICE_UP(JI,JK+D%NKL)
 
       ! Compute the updraft theta_v, buoyancy and w**2 for level JK+1   
-      !PTHV_UP(:,JK+KKL) = PTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL)))
+      !PTHV_UP(D%NIB:D%NIE,JK+KKL) = PTH_UP(D%NIB:D%NIE,JK+KKL)*((1+ZRVORD*PRV_UP(D%NIB:D%NIE,JK+KKL))/(1+PRT_UP(D%NIB:D%NIE,JK+KKL)))
       !PTHV_UP(JI,JK+KKL) = ZTH_UP(JI,JK+KKL)*(1.+0.608*PRV_UP(JI,JK+KKL) - PRC_UP(JI,JK+KKL))
       !! A corriger pour utiliser q et non r !!!!      
       !ZMIX1(JI)=ZZDZ(JI,JK)*(PENTR(JI,JK)-PDETR(JI,JK))
@@ -612,11 +612,11 @@ 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,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)
 !$mnh_expand_array(JI=D%NIB:D%NIE)
-PEMF(:,D%NKB) =0.
+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.
@@ -630,16 +630,16 @@ DO JI=D%NIB,D%NIE
 ENDDO
 
 !$mnh_expand_array(JI=D%NIB:D%NIE)
-GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
+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=D%NKTB,D%NKTE
+DO JK=1,D%NKT
   !$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.)
+  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
-DO JK=D%NKTB, D%NKTE
+DO JK=1, D%NKT
   DO JI=D%NIB,D%NIE
     IF (GWORK2(JI,JK)) THEN
       PEMF(JI,JK)     = PEMF(JI,JK)     * ZCOEF(JI,JK) 
diff --git a/src/common/turb/mode_mf_turb.F90 b/src/common/turb/mode_mf_turb.F90
index 761a73faa5b97e103018b21ea7e24c299fca6022..537e16c14db54119da557d425da69494ccd0b387 100644
--- a/src/common/turb/mode_mf_turb.F90
+++ b/src/common/turb/mode_mf_turb.F90
@@ -156,17 +156,17 @@ CALL MZM_MF(D, PRTM(:,:), PFLXZRMF(:,:))
 CALL MZM_MF(D, PTHVM(:,:), PFLXZTHVMF(:,:))
 
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-PFLXZTHMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-PFLXZTHMF(:,:))
-PFLXZRMF(:,:) =  PEMF(:,:)*(PRT_UP(:,:)-PFLXZRMF(:,:))
-PFLXZTHVMF(:,:) = PEMF(:,:)*(PTHV_UP(:,:)-PFLXZTHVMF(:,:))
+PFLXZTHMF(D%NIB:D%NIE,:) = PEMF(D%NIB:D%NIE,:)*(PTHL_UP(D%NIB:D%NIE,:)-PFLXZTHMF(D%NIB:D%NIE,:))
+PFLXZRMF(D%NIB:D%NIE,:) =  PEMF(D%NIB:D%NIE,:)*(PRT_UP(D%NIB:D%NIE,:)-PFLXZRMF(D%NIB:D%NIE,:))
+PFLXZTHVMF(D%NIB:D%NIE,:) = PEMF(D%NIB:D%NIE,:)*(PTHV_UP(D%NIB:D%NIE,:)-PFLXZTHVMF(D%NIB:D%NIE,:))
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 IF (OMIXUV) THEN
   CALL MZM_MF(D, PUM(:,:), PFLXZUMF(:,:))
   CALL MZM_MF(D, PVM(:,:), PFLXZVMF(:,:))
   !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  PFLXZUMF(:,:) =  PEMF(:,:)*(PU_UP(:,:)-PFLXZUMF(:,:))
-  PFLXZVMF(:,:) =  PEMF(:,:)*(PV_UP(:,:)-PFLXZVMF(:,:))
+  PFLXZUMF(D%NIB:D%NIE,:) =  PEMF(D%NIB:D%NIE,:)*(PU_UP(D%NIB:D%NIE,:)-PFLXZUMF(D%NIB:D%NIE,:))
+  PFLXZVMF(D%NIB:D%NIE,:) =  PEMF(D%NIB:D%NIE,:)*(PV_UP(D%NIB:D%NIE,:)-PFLXZVMF(D%NIB:D%NIE,:))
   !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ELSE
   PFLXZUMF(:,:) = 0.
@@ -191,8 +191,8 @@ CALL TRIDIAG_MASSFLUX(D,PTHLM,PFLXZTHMF,-PEMF,PTSTEP,PIMPL,  &
 ! compute new flux and THL tendency
 CALL MZM_MF(D, ZVARS(:,:), PFLXZTHMF(:,:))
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-PFLXZTHMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-PFLXZTHMF(:,:))
-PTHLDT(:,:)= (ZVARS(:,:)-PTHLM(:,:))/PTSTEP
+PFLXZTHMF(D%NIB:D%NIE,:) = PEMF(D%NIB:D%NIE,:)*(PTHL_UP(D%NIB:D%NIE,:)-PFLXZTHMF(D%NIB:D%NIE,:))
+PTHLDT(D%NIB:D%NIE,:)= (ZVARS(D%NIB:D%NIE,:)-PTHLM(D%NIB:D%NIE,:))/PTSTEP
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 !
@@ -203,8 +203,8 @@ CALL TRIDIAG_MASSFLUX(D,PRTM(:,:),PFLXZRMF,-PEMF,PTSTEP,PIMPL,  &
 ! compute new flux and RT tendency
 CALL MZM_MF(D, ZVARS(:,:), PFLXZRMF(:,:))
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-PFLXZRMF(:,:) =  PEMF(:,:)*(PRT_UP(:,:)-PFLXZRMF(:,:))
-PRTDT(:,:) = (ZVARS(:,:)-PRTM(:,:))/PTSTEP
+PFLXZRMF(D%NIB:D%NIE,:) =  PEMF(D%NIB:D%NIE,:)*(PRT_UP(D%NIB:D%NIE,:)-PFLXZRMF(D%NIB:D%NIE,:))
+PRTDT(D%NIB:D%NIE,:) = (ZVARS(D%NIB:D%NIE,:)-PRTM(D%NIB:D%NIE,:))/PTSTEP
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 !
 
@@ -219,8 +219,8 @@ IF (OMIXUV) THEN
   ! compute new flux and U tendency
   CALL MZM_MF(D, ZVARS(:,:), PFLXZUMF(:,:))
   !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  PFLXZUMF(:,:) = PEMF(:,:)*(PU_UP(:,:)-PFLXZUMF(:,:))
-  PUDT(:,:)= (ZVARS(:,:)-PUM(:,:))/PTSTEP
+  PFLXZUMF(D%NIB:D%NIE,:) = PEMF(D%NIB:D%NIE,:)*(PU_UP(D%NIB:D%NIE,:)-PFLXZUMF(D%NIB:D%NIE,:))
+  PUDT(D%NIB:D%NIE,:)= (ZVARS(D%NIB:D%NIE,:)-PUM(D%NIB:D%NIE,:))/PTSTEP
   !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
   !
   !
@@ -233,8 +233,8 @@ IF (OMIXUV) THEN
   ! compute new flux and V tendency
   CALL MZM_MF(D, ZVARS(:,:), PFLXZVMF(:,:))
   !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  PFLXZVMF(:,:) = PEMF(:,:)*(PV_UP(:,:)-PFLXZVMF(:,:))
-  PVDT(:,:)= (ZVARS(:,:)-PVM(:,:))/PTSTEP
+  PFLXZVMF(D%NIB:D%NIE,:) = PEMF(D%NIB:D%NIE,:)*(PV_UP(D%NIB:D%NIE,:)-PFLXZVMF(D%NIB:D%NIE,:))
+  PVDT(D%NIB:D%NIE,:)= (ZVARS(D%NIB:D%NIE,:)-PVM(D%NIB:D%NIE,:))/PTSTEP
   !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ELSE
   PUDT(:,:)=0.
@@ -250,7 +250,7 @@ DO JSV=1,KSV
 
   CALL MZM_MF(D, PSVM(:,:,JSV), PFLXZSVMF(:,:,JSV))
   !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  PFLXZSVMF(:,:,JSV) = PEMF(:,:)*(PSV_UP(:,:,JSV)-PFLXZSVMF(:,:,JSV))
+  PFLXZSVMF(D%NIB:D%NIE,:,JSV) = PEMF(D%NIB:D%NIE,:)*(PSV_UP(D%NIB:D%NIE,:,JSV)-PFLXZSVMF(D%NIB:D%NIE,:,JSV))
   !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
   !
   ! 3.5 Compute the tendency for scalar variables
@@ -261,8 +261,8 @@ DO JSV=1,KSV
   ! compute new flux and Sv tendency
   CALL MZM_MF(D, ZVARS, PFLXZSVMF(:,:,JSV))
   !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  PFLXZSVMF(:,:,JSV) = PEMF(:,:)*(PSV_UP(:,:,JSV)-PFLXZSVMF(:,:,JSV))
-  PSVDT(:,:,JSV)= (ZVARS(:,:)-PSVM(:,:,JSV))/PTSTEP
+  PFLXZSVMF(D%NIB:D%NIE,:,JSV) = PEMF(D%NIB:D%NIE,:)*(PSV_UP(D%NIB:D%NIE,:,JSV)-PFLXZSVMF(D%NIB:D%NIE,:,JSV))
+  PSVDT(D%NIB:D%NIE,:,JSV)= (ZVARS(D%NIB:D%NIE,:)-PSVM(D%NIB:D%NIE,:,JSV))/PTSTEP
   !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 ENDDO
diff --git a/src/common/turb/mode_mf_turb_expl.F90 b/src/common/turb/mode_mf_turb_expl.F90
index 7332c477ae257a847b9cab26d651d05adc6eef7b..cace0eb3ac7af50ab043dd2e25d2a56ba4063c0b 100644
--- a/src/common/turb/mode_mf_turb_expl.F90
+++ b/src/common/turb/mode_mf_turb_expl.F90
@@ -131,29 +131,29 @@ PVDT   = 0.
 CALL MZM_MF(D, PRTM (:,:), ZRTM_F(:,:))
 CALL MZM_MF(D, PTHLM(:,:), ZTHLM_F(:,:))
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-ZQTM(:,:)   = ZRTM_F(:,:)/(1.+ZRTM_F(:,:))
-ZQT_UP(:,:) = PRT_UP(:,:)/(1.+PRT_UP(:,:))
-ZTHS_UP(:,:)= PTHL_UP(:,:)*(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(:,:))
-ZTHSM(:,:)  = ZTHLM_F(:,:)*(1.+PARAMMF%XLAMBDA_MF*ZQTM(:,:))
+ZQTM(D%NIB:D%NIE,:)   = ZRTM_F(D%NIB:D%NIE,:)/(1.+ZRTM_F(D%NIB:D%NIE,:))
+ZQT_UP(D%NIB:D%NIE,:) = PRT_UP(D%NIB:D%NIE,:)/(1.+PRT_UP(D%NIB:D%NIE,:))
+ZTHS_UP(D%NIB:D%NIE,:)= PTHL_UP(D%NIB:D%NIE,:)*(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(D%NIB:D%NIE,:))
+ZTHSM(D%NIB:D%NIE,:)  = ZTHLM_F(D%NIB:D%NIE,:)*(1.+PARAMMF%XLAMBDA_MF*ZQTM(D%NIB:D%NIE,:))
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 CALL MZM_MF(D, PTHLM(:,:), PFLXZTHLMF(:,:))
 CALL MZM_MF(D, PRTM(:,:), PFLXZRMF(:,:))
 CALL MZM_MF(D, PTHVM(:,:), PFLXZTHVMF(:,:))
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-PFLXZTHLMF(:,:)  = PEMF(:,:)*(PTHL_UP(:,:)-PFLXZTHLMF(:,:))  ! ThetaL
-PFLXZRMF(:,:)    = PEMF(:,:)*(PRT_UP(:,:)-PFLXZRMF(:,:))  ! Rt
-PFLXZTHVMF(:,:)  = PEMF(:,:)*(PTHV_UP(:,:)-PFLXZTHVMF(:,:))  ! ThetaV
+PFLXZTHLMF(D%NIB:D%NIE,:)  = PEMF(D%NIB:D%NIE,:)*(PTHL_UP(D%NIB:D%NIE,:)-PFLXZTHLMF(D%NIB:D%NIE,:))  ! ThetaL
+PFLXZRMF(D%NIB:D%NIE,:)    = PEMF(D%NIB:D%NIE,:)*(PRT_UP(D%NIB:D%NIE,:)-PFLXZRMF(D%NIB:D%NIE,:))  ! Rt
+PFLXZTHVMF(D%NIB:D%NIE,:)  = PEMF(D%NIB:D%NIE,:)*(PTHV_UP(D%NIB:D%NIE,:)-PFLXZTHVMF(D%NIB:D%NIE,:))  ! ThetaV
 
-ZFLXZTHSMF(:,:)  = PEMF(:,:)*(ZTHS_UP(:,:)-ZTHSM(:,:))    ! Theta S flux
+ZFLXZTHSMF(D%NIB:D%NIE,:)  = PEMF(D%NIB:D%NIE,:)*(ZTHS_UP(D%NIB:D%NIE,:)-ZTHSM(D%NIB:D%NIE,:))    ! Theta S flux
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 IF (OMIXUV) THEN
   CALL MZM_MF(D, PUM(:,:), PFLXZUMF(:,:))
   CALL MZM_MF(D, PVM(:,:), PFLXZVMF(:,:))
   !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  PFLXZUMF(:,:) =  PEMF(:,:)*(PU_UP(:,:)-PFLXZUMF(:,:))  ! U
-  PFLXZVMF(:,:) =  PEMF(:,:)*(PV_UP(:,:)-PFLXZVMF(:,:))  ! V
+  PFLXZUMF(D%NIB:D%NIE,:) =  PEMF(D%NIB:D%NIE,:)*(PU_UP(D%NIB:D%NIE,:)-PFLXZUMF(D%NIB:D%NIE,:))  ! U
+  PFLXZVMF(D%NIB:D%NIE,:) =  PEMF(D%NIB:D%NIE,:)*(PV_UP(D%NIB:D%NIE,:)-PFLXZVMF(D%NIB:D%NIE,:))  ! V
   !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ELSE
   PFLXZUMF(:,:) = 0.
diff --git a/src/common/turb/mode_thl_rt_from_th_r_mf.F90 b/src/common/turb/mode_thl_rt_from_th_r_mf.F90
index c2fb615c87fd9b367a3c512db174ba31f4d5ccf8..c83a0ff578f714f84542d789806b55590b454cce 100644
--- a/src/common/turb/mode_thl_rt_from_th_r_mf.F90
+++ b/src/common/turb/mode_thl_rt_from_th_r_mf.F90
@@ -78,48 +78,63 @@ REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT)  :: PRT      ! total non precip. water
 !----------------------------------------------------------------------------
 REAL, DIMENSION(D%NIT,D%NKT) :: ZCP, ZT
 REAL, DIMENSION(D%NIT,D%NKT) :: ZLVOCPEXN, ZLSOCPEXN
-INTEGER :: JRR
+INTEGER :: JRR, JI, JK
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !----------------------------------------------------------------------------
 !
 !
 IF (LHOOK) CALL DR_HOOK('THL_RT_FRM_TH_R_MF',0,ZHOOK_HANDLE)
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 !temperature
-ZT(:,:) = PTH(:,:) * PEXN(:,:)
+ZT(D%NIB:D%NIE,:) = PTH(D%NIB:D%NIE,:) * PEXN(D%NIB:D%NIE,:)
 
 !Cp
-ZCP=CST%XCPD
-IF (KRR > 0) ZCP(:,:) = ZCP(:,:) + CST%XCPV * PR(:,:,1)
+ZCP(D%NIB:D%NIE,:)=CST%XCPD
+IF (KRR > 0) ZCP(D%NIB:D%NIE,:) = ZCP(D%NIB:D%NIE,:) + CST%XCPV * PR(D%NIB:D%NIE,:,1)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 DO JRR = 2,1+KRRL  ! loop on the liquid components
-  ZCP(:,:)  = ZCP(:,:) + CST%XCL * PR(:,:,JRR)
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+  ZCP(D%NIB:D%NIE,:)  = ZCP(D%NIB:D%NIE,:) + CST%XCL * PR(D%NIB:D%NIE,:,JRR)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 END DO
 DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
-  ZCP(:,:)  = ZCP(:,:)  + CST%XCI * PR(:,:,JRR)
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
+  ZCP(D%NIB:D%NIE,:)  = ZCP(D%NIB:D%NIE,:)  + CST%XCI * PR(D%NIB:D%NIE,:,JRR)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 END DO
 
 IF ( KRRL >= 1 ) THEN
   IF ( KRRI >= 1 ) THEN
+    !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
     !ZLVOCPEXN and ZLSOCPEXN
-    ZLVOCPEXN(:,:)=(CST%XLVTT + (CST%XCPV-CST%XCL) *  (ZT(:,:)-CST%XTT) ) / ZCP(:,:) / PEXN(:,:)
-    ZLSOCPEXN(:,:)=(CST%XLSTT + (CST%XCPV-CST%XCI) *  (ZT(:,:)-CST%XTT) ) / ZCP(:,:) / PEXN(:,:)
+    ZLVOCPEXN(D%NIB:D%NIE,:)=(CST%XLVTT + (CST%XCPV-CST%XCL) *  (ZT(D%NIB:D%NIE,:)-CST%XTT) ) & 
+                            &/ ZCP(D%NIB:D%NIE,:) / PEXN(D%NIB:D%NIE,:)
+    ZLSOCPEXN(D%NIB:D%NIE,:)=(CST%XLSTT + (CST%XCPV-CST%XCI) *  (ZT(D%NIB:D%NIE,:)-CST%XTT) ) &
+                            &/ ZCP(D%NIB:D%NIE,:) / PEXN(D%NIB:D%NIE,:)
     ! Rnp 
-    PRT(:,:)  = PR(:,:,1)  + PR(:,:,2)  + PR(:,:,4)
+    PRT(D%NIB:D%NIE,:)  = PR(D%NIB:D%NIE,:,1)  + PR(D%NIB:D%NIE,:,2)  + PR(D%NIB:D%NIE,:,4)
     ! Theta_l 
-    PTHL(:,:)  = PTH(:,:)  - ZLVOCPEXN(:,:) * PR(:,:,2) &
-                           - ZLSOCPEXN(:,:) * PR(:,:,4)
+    PTHL(D%NIB:D%NIE,:)  = PTH(D%NIB:D%NIE,:)  - ZLVOCPEXN(D%NIB:D%NIE,:) * PR(D%NIB:D%NIE,:,2) &
+                           - ZLSOCPEXN(D%NIB:D%NIE,:) * PR(D%NIB:D%NIE,:,4)
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
   ELSE
+    !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
     !ZLVOCPEXN
-    ZLVOCPEXN(:,:)=(CST%XLVTT + (CST%XCPV-CST%XCL) *  (ZT(:,:)-CST%XTT) ) / ZCP(:,:) / PEXN(:,:)
+    ZLVOCPEXN(D%NIB:D%NIE,:)=(CST%XLVTT + (CST%XCPV-CST%XCL) *  (ZT(D%NIB:D%NIE,:)-CST%XTT) ) &
+                            &/ ZCP(D%NIB:D%NIE,:) / PEXN(D%NIB:D%NIE,:)
     ! Rnp
-    PRT(:,:)  = PR(:,:,1)  + PR(:,:,2) 
+    PRT(D%NIB:D%NIE,:)  = PR(D%NIB:D%NIE,:,1)  + PR(D%NIB:D%NIE,:,2) 
     ! Theta_l
-    PTHL(:,:) = PTH(:,:)  - ZLVOCPEXN(:,:) * PR(:,:,2)
+    PTHL(D%NIB:D%NIE,:) = PTH(D%NIB:D%NIE,:)  - ZLVOCPEXN(D%NIB:D%NIE,:) * PR(D%NIB:D%NIE,:,2)
+    !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
   END IF
 ELSE
+  !$mnh_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
   ! Rnp = rv
-  PRT(:,:)  = PR(:,:,1)
+  PRT(D%NIB:D%NIE,:)  = PR(D%NIB:D%NIE,:,1)
   ! Theta_l = Theta
-  PTHL(:,:) = PTH(:,:)
+  PTHL(D%NIB:D%NIE,:) = PTH(D%NIB:D%NIE,:)
+  !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=1:D%NKT)
 END IF
 IF (LHOOK) CALL DR_HOOK('THL_RT_FRM_TH_R_MF',1,ZHOOK_HANDLE)
 END SUBROUTINE THL_RT_FROM_TH_R_MF
diff --git a/src/common/turb/shallow_mf.F90 b/src/common/turb/shallow_mf.F90
index f3ff7848d42cf7a1c9c49d3b7974ece4dd5bc868..c59028930612593fdcec07a62d02b97dc047fcec 100644
--- a/src/common/turb/shallow_mf.F90
+++ b/src/common/turb/shallow_mf.F90
@@ -171,6 +171,7 @@ REAL, DIMENSION(D%NIT,D%NKT) ::     &
           ZEMF_O_RHODREF,                         & ! entrainment/detrainment
           ZBUO_INTEG                                ! integrated buoyancy
 REAL, DIMENSION(D%NIT,D%NKT) :: ZFRAC_ICE
+REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
 
 REAL, DIMENSION(D%NIT,D%NKT,KSV) ::  &
                                           ZSV_UP,&  ! updraft scalar var.
@@ -201,12 +202,15 @@ ENDIF
 ZFRAC_ICE(:,:) = 0.
 IF (KRR.GE.4) THEN
   !$mnh_expand_where(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-  WHERE(PRM(:,:,2)+PRM(:,:,4) > 1.E-20)
-    ZFRAC_ICE(:,:) = PRM(:,:,4) / (PRM(:,:,2)+PRM(:,:,4))
+  WHERE(PRM(D%NIB:D%NIE,:,2)+PRM(D%NIB:D%NIE,:,4) > 1.E-20)
+    ZFRAC_ICE(D%NIB:D%NIE,:) = PRM(D%NIB:D%NIE,:,4) / (PRM(D%NIB:D%NIE,:,2)+PRM(D%NIB:D%NIE,:,4))
   ENDWHERE
   !$mnh_end_expand_where(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ENDIF
-CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEB,ZFRAC_ICE(:,:),PTHM(:,:)*PEXNM(:,:), IERR(:,:))
+!$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+ZWK(D%NIB:D%NIE,:)=PTHM(D%NIB:D%NIE,:)*PEXNM(D%NIB:D%NIE,:)
+!$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
+CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEB,ZFRAC_ICE(:,:),ZWK(:,:), IERR(:,:))
 
 ! Conservative variables at t-dt
 CALL THL_RT_FROM_TH_R_MF(D, CST, KRR,KRRL,KRRI,    &
@@ -215,7 +219,7 @@ CALL THL_RT_FROM_TH_R_MF(D, CST, KRR,KRRL,KRRI,    &
 
 ! Virtual potential temperature at t-dt
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-ZTHVM(:,:) = PTHM(:,:)*((1.+CST%XRV / CST%XRD *PRM(:,:,1))/(1.+ZRTM(:,:))) 
+ZTHVM(D%NIB:D%NIE,:) = PTHM(D%NIB:D%NIE,:)*((1.+CST%XRV / CST%XRD *PRM(D%NIB:D%NIE,:,1))/(1.+ZRTM(D%NIB:D%NIE,:))) 
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 ! 
 !!! 2. Compute updraft
@@ -288,7 +292,7 @@ CALL COMPUTE_MF_CLOUD(D, CST, PARAMMF, KRR, KRRL, KRRI, &
 !!!    ------------------------------------------------------------------------
 !
 !$mnh_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
-ZEMF_O_RHODREF(:,:)=PEMF(:,:)/PRHODREF(:,:)
+ZEMF_O_RHODREF(D%NIB:D%NIE,:)=PEMF(D%NIB:D%NIE,:)/PRHODREF(D%NIB:D%NIE,:)
 !$mnh_end_expand_array(JI=D%NIB:D%NIE,JK=D%NKTB:D%NKTE)
 
 IF ( PIMPL_MF > 1.E-10 ) THEN  
diff --git a/src/common/turb/th_r_from_thl_rt.func.h b/src/common/turb/th_r_from_thl_rt.func.h
index 62cbcd5c79fd7af3ac78ffaf24ed9c13783a754c..5a571b688d6ca46c05818bd43cc87947491f043a 100644
--- a/src/common/turb/th_r_from_thl_rt.func.h
+++ b/src/common/turb/th_r_from_thl_rt.func.h
@@ -5,7 +5,7 @@
       SUBROUTINE TH_R_FROM_THL_RT(CST, NEB, KT, HFRAC_ICE,PFRAC_ICE,PP, &
                                   PTHL, PRT, PTH, PRV, PRL, PRI,           &
                                   PRSATW, PRSATI, PRR, PRS, PRG, PRH, OOCEAN,&
-                                  PBUF)
+                                  PBUF, KB, KE)
 ! ******* TO BE INCLUDED IN THE *CONTAINS* OF A SUBROUTINE, IN ORDER TO EASE AUTOMATIC INLINING ******
 ! => Don't use drHook !!!
 ! "compute_frac_ice.func.h" must be included at the same time
@@ -74,13 +74,15 @@ REAL, DIMENSION(KT), INTENT(INOUT):: PRI    ! vapor mixing ratio
 REAL, DIMENSION(KT), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
 REAL, DIMENSION(KT), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
 REAL, DIMENSION(KT, 16), INTENT(OUT) :: PBUF ! buffer to replace automatic arrays
+INTEGER, OPTIONAL :: KB !first index to deal with (default is 1)
+INTEGER, OPTIONAL :: KE !last index to deal with (default if KT)
 !
 !-------------------------------------------------------------------------------
 !
 !       0.2  declaration of local variables
 INTEGER                       :: II ! Loop control
 INTEGER                       :: JITER ! number of iterations
-INTEGER                       :: J
+INTEGER                       :: J, IB, IE
 INTEGER, PARAMETER :: IEXN=1, IRVSAT=2, ICPH=3, IRLTEMP=4, ICPH2=5, IT=6, ILVOCPEXN=7, ILSOCPEXN=8, &
                     & IDRSATODT=9, IDRSATODTW=10, IDRSATODTI=11, IFOESW=12, IFOESI=13, &
                     & ILOGT=14, I99PP=15, I1PRT=16
@@ -93,20 +95,22 @@ REAL :: ZVAR1, ZVAR2, ZTPOW2, ZDELT
 !
 !
 !
+IB=MERGE(KB, 1, PRESENT(KB))
+IE=MERGE(KE, KT, PRESENT(KE))
 !Number of iterations
 JITER=2
 !
-!Computation of PBUF(:, ICPH2) depending on dummy arguments received
-PBUF(:, ICPH2)=0
-IF(PRESENT(PRR)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCL*PRR(:)
-IF(PRESENT(PRS)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCI*PRS(:)
-IF(PRESENT(PRG)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCI*PRG(:)
-IF(PRESENT(PRH)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCI*PRH(:)
+!Computation of PBUF(IB:IE, ICPH2) depending on dummy arguments received
+PBUF(IB:IE, ICPH2)=0
+IF(PRESENT(PRR)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCL*PRR(IB:IE)
+IF(PRESENT(PRS)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCI*PRS(IB:IE)
+IF(PRESENT(PRG)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCI*PRG(IB:IE)
+IF(PRESENT(PRH)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCI*PRH(IB:IE)
 !
 !Computation of an approximate state thanks to PRL and PRI guess
-PBUF(:, IEXN)=(PP(:)/CST%XP00) ** CST%RDSCPD
+PBUF(IB:IE, IEXN)=(PP(IB:IE)/CST%XP00) ** CST%RDSCPD
 
-DO J=1,KT
+DO J=IB,IE
   PBUF(J, I99PP)=0.99*PP(J)
   PRV(J)=PRT(J)-PRL(J)-PRI(J)
   PBUF(J, ICPH)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + PBUF(J, ICPH2)
@@ -124,27 +128,27 @@ ENDDO
 
 DO II=1,JITER
   IF (OOCEAN) THEN
-    PBUF(:, IT)=PTH(:)
+    PBUF(IB:IE, IT)=PTH(IB:IE)
   ELSE
-    PBUF(:, IT)=PTH(:)*PBUF(:, IEXN)
+    PBUF(IB:IE, IT)=PTH(IB:IE)*PBUF(IB:IE, IEXN)
   END IF
   !Computation of liquid/ice fractions
-  PFRAC_ICE(:) = 0.
-  DO J=1, KT
+  PFRAC_ICE(IB:IE) = 0.
+  DO J=IB, IE
     IF(PRL(J)+PRI(J) > 1.E-20) THEN
       PFRAC_ICE(J) = PRI(J) / (PRL(J)+PRI(J))
     ENDIF
   ENDDO
-  CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEB,PFRAC_ICE(:),PBUF(:, IT))
+  CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEB,PFRAC_ICE(IB:IE),PBUF(IB:IE, IT))
 
   !Computation of Rvsat and dRsat/dT
   !In this version QSAT, QSATI, DQSAT and DQASATI functions are not used
   !due to performance issue
 
   ! Log does not vectorize on all compilers:
-  PBUF(:, ILOGT)=LOG(PBUF(:, IT))
+  PBUF(IB:IE, ILOGT)=LOG(PBUF(IB:IE, IT))
 
-  DO J=1,KT
+  DO J=IB, IE
     PBUF(J, IFOESW) = MIN(EXP( CST%XALPW - CST%XBETAW/PBUF(J, IT) - CST%XGAMW*PBUF(J, ILOGT)  ), PBUF(J, I99PP))
     PBUF(J, IFOESI) = MIN(EXP( CST%XALPI - CST%XBETAI/PBUF(J, IT) - CST%XGAMI*PBUF(J, ILOGT)  ), PBUF(J, I99PP))
     PRSATW(J) = CST%XRD/CST%XRV*PBUF(J, IFOESW)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESW)/PP(J))