From 3615ca4adc682603d9e2e60bda0a53ec10bed31e Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Fri, 16 Dec 2022 11:43:04 +0100
Subject: [PATCH] Quentin 16/12/2022: bugfix uniformize use of structure D% by
 using local integers

---
 src/common/aux/gradient_m_phy.F90             |   9 +-
 src/common/aux/gradient_u_phy.F90             |   2 +-
 .../micro/mode_ice4_sedimentation_split.F90   |   8 +
 .../micro/mode_ice4_sedimentation_stat.F90    |   4 +-
 src/common/turb/mode_compute_bl89_ml.F90      |   6 +-
 .../turb/mode_compute_mf_cloud_direct.F90     |   6 +-
 src/common/turb/mode_compute_updraft.F90      | 228 ++++++++--------
 src/common/turb/mode_compute_updraft_raha.F90 |   4 +-
 .../turb/mode_compute_updraft_rhcj10.F90      | 256 +++++++++---------
 src/common/turb/mode_mf_turb_expl.F90         |  18 +-
 src/common/turb/mode_prandtl.F90              |   9 +-
 src/common/turb/mode_turb_ver_thermo_corr.F90 |   6 +-
 src/common/turb/shuman_mf.F90                 |  12 +-
 13 files changed, 292 insertions(+), 276 deletions(-)

diff --git a/src/common/aux/gradient_m_phy.F90 b/src/common/aux/gradient_m_phy.F90
index 7cd6eb087..2b1636b63 100644
--- a/src/common/aux/gradient_m_phy.F90
+++ b/src/common/aux/gradient_m_phy.F90
@@ -69,7 +69,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)                :: PY       ! vari
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGZ_M_W  ! result at flux
                                                               ! side
 !
-INTEGER :: IKT,IKTB,IKTE,IIB,IJB,IIE,IJE,IKA,IKU
+INTEGER :: IKT,IKTB,IKTE,IIB,IJB,IIE,IJE,IKA,IKU,IKL
 INTEGER :: JI,JJ,JK
 !-------------------------------------------------------------------------------
 !
@@ -86,6 +86,7 @@ IJB=D%NJBC
 IKT=D%NKT
 IKA=D%NKA
 IKU=D%NKU
+IKL=D%NKL
 !
 DO JK=IKTB,IKTE 
   DO JJ=IJB,IJE 
@@ -318,6 +319,7 @@ IIE=D%NIEC
 IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
+IKT=D%NKT
 !
 !----------------------------------------------------------------------------
 !
@@ -435,7 +437,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGX_M_U  ! result at flux
                                                               ! side
 REAL, DIMENSION(D%NIT*D%NJT*D%NKT) :: ZGX_M_U
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT):: ZY, ZDXX,ZDZZ,ZDZX
-INTEGER  IIU,IKU,JI,JK,IKL, IKA, IKU
+INTEGER  IIU,IKU,JI,JK,IKL, IKA
 !
 INTEGER :: JJK,IJU
 INTEGER :: JIJK,JIJKOR,JIJKEND
@@ -454,7 +456,6 @@ IJU=D%NJT
 IKU=D%NKT
 IKL=D%NKL
 IKA=D%NKA
-IKU=D%NKU
 IF (.NOT. OFLAT) THEN
   JIJKOR  = 1 + JPHEXT + IIU*IJU*(JPVEXT_TURB+1 - 1)
   JIJKEND = IIU*IJU*(IKU-JPVEXT_TURB)
@@ -587,7 +588,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT),INTENT(OUT) :: PGY_M_V  ! result at flux
                                                               ! side
 !REAL, DIMENSION(D%NIT*D%NJT*D%NKT) :: ZGY_M_V
 !REAL, DIMENSION(D%NIT,D%NJT,D%NKT):: ZY, ZDYY,ZDZZ,ZDZY
-INTEGER  IJU,IKU,JI,JJ,JK,IKL, IKA, IKU
+INTEGER  IJU,IKU,JI,JJ,JK,IKL, IKA
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/common/aux/gradient_u_phy.F90 b/src/common/aux/gradient_u_phy.F90
index 75ce7dac0..ff685a0c1 100644
--- a/src/common/aux/gradient_u_phy.F90
+++ b/src/common/aux/gradient_u_phy.F90
@@ -69,7 +69,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDZZ    ! metric coefficient
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGZ_U_UW ! result UW point
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: PA_WORK, PDZZ_WORK
 !
-INTEGER :: JI,JJ,JK, IIB, IIE, IJB, IJE
+INTEGER :: JI,JJ,JK, IIB, IIE, IJB, IJE,IKT
 
 !
 !*       0.2   declaration of local variables
diff --git a/src/common/micro/mode_ice4_sedimentation_split.F90 b/src/common/micro/mode_ice4_sedimentation_split.F90
index 2c078a26e..a7cb83129 100644
--- a/src/common/micro/mode_ice4_sedimentation_split.F90
+++ b/src/common/micro/mode_ice4_sedimentation_split.F90
@@ -89,6 +89,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), OPTIONAL, INTENT(OUT)   :: PFPR    ! upp
 !
 !
 INTEGER                                                             :: JI, JJ, JK
+INTEGER :: IKTB, IKTE, IKB, IKL, IIE, IIB, IJB, IJE
 INTEGER                                                             :: IRR !Workaround of PGI bug with OpenACC (at least up to 18.10 version)
 LOGICAL                                                             :: GSEDIC !Workaround of PGI bug with OpenACC (at least up to 18.10 version)
 LOGICAL                                                             :: GPRESENT_PFPR, GPRESENT_PSEA
@@ -116,6 +117,13 @@ IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_SPLIT', 0, ZHOOK_HANDLE)
 GSEDIC = OSEDIC
 IRR    = KRR
 !
+IKTB=D%NKTB
+IKTE=D%NKTE
+IIB=D%NIB
+IIE=D%NIE
+IJB=D%NJB
+IJE=D%NJE
+!
 IF (PRESENT(PFPR)) THEN
   GPRESENT_PFPR = .TRUE.
 ELSE
diff --git a/src/common/micro/mode_ice4_sedimentation_stat.F90 b/src/common/micro/mode_ice4_sedimentation_stat.F90
index 8ce021921..33e4e6e01 100644
--- a/src/common/micro/mode_ice4_sedimentation_stat.F90
+++ b/src/common/micro/mode_ice4_sedimentation_stat.F90
@@ -89,7 +89,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), OPTIONAL, INTENT(OUT)   :: PFPR    ! upp
 !*       0.2  declaration of local variables
 !
 LOGICAL :: LLSEA_AND_TOWN
-INTEGER :: JRR, JI, JJ, JK, IKB, IKE,IKL, IIE, IIB, IJB, IJE
+INTEGER :: JRR, JI, JJ, JK, IKB, IKE,IKL, IIE, IIB, IJB, IJE, IKTB, IKTE
 INTEGER :: ISHIFT, IK, IKPLUS
 REAL :: ZQP, ZP1, ZINVTSTEP, ZGAC, ZGC, ZGAC2, ZGC2, ZRAYDEFO
 REAL, DIMENSION(D%NIT)     :: ZWSEDW1, ZWSEDW2 ! sedimentation speed
@@ -119,6 +119,8 @@ IIB=D%NIB
 IIE=D%NIE
 IJB=D%NJB
 IJE=D%NJE
+IKTB=D%NKTB
+IKTE=D%NKTE
 !
 IF ( PRESENT( PFPR ) ) THEN
  !Set to 0. to avoid undefined values (in files)
diff --git a/src/common/turb/mode_compute_bl89_ml.F90 b/src/common/turb/mode_compute_bl89_ml.F90
index ce14d0b56..36008959d 100644
--- a/src/common/turb/mode_compute_bl89_ml.F90
+++ b/src/common/turb/mode_compute_bl89_ml.F90
@@ -133,7 +133,7 @@ IF (OUPORDN.EQV..TRUE.) THEN
    ZVPT_DEP(IIJB:IIJE)=ZHLVPT(IIJB:IIJE,KK) ! departure point is on flux level
    !$mnh_end_expand_array(JIJ=IIJB:IIJE)
    !We must compute what happens between flux level KK and mass level KK
-   DO J1D=D%NIJB,D%NIJE
+   DO J1D=IIJB,IIJE
      ZTEST0=0.5+SIGN(0.5,ZINTE(J1D)) ! test if there's energy to consume
      ! Energy consumed if parcel cross the entire layer
      ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D)      *      &
@@ -169,7 +169,7 @@ IF (OUPORDN.EQV..TRUE.) THEN
  DO JKK=KK+IKL,IKE,IKL
     IF(ZTESTM > 0.) THEN
       ZTESTM=0
-      DO J1D=D%NIJB,D%NIJE
+      DO J1D=IIJB,IIJE
         ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
         ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D)      *      &
             (ZHLVPT(J1D,JKK) - ZVPT_DEP(J1D))   &
@@ -210,7 +210,7 @@ IF (OUPORDN.EQV..FALSE.) THEN
  DO JKK=KK,IKB,-IKL
     IF(ZTESTM > 0.) THEN
       ZTESTM=0
-      DO J1D=D%NIJB,D%NIJE
+      DO J1D=IIJB,IIJE
         ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
          ZPOTE(J1D) = ZTEST0*(-PG_O_THVREF(J1D)      *      &
             (ZHLVPT(J1D,JKK) - PVPT(J1D,KK)) &
diff --git a/src/common/turb/mode_compute_mf_cloud_direct.F90 b/src/common/turb/mode_compute_mf_cloud_direct.F90
index 0a2f7def9..2323b77f5 100644
--- a/src/common/turb/mode_compute_mf_cloud_direct.F90
+++ b/src/common/turb/mode_compute_mf_cloud_direct.F90
@@ -72,7 +72,7 @@ REAL, DIMENSION(D%NIJT,D%NKT),   INTENT(OUT)  :: PCF_MF         ! and cloud frac
 !
 !*                    0.1  Declaration of local variables
 !
-INTEGER  :: JI,JK, JK0, IKB,IKE,IKL
+INTEGER  :: JI,JK, JK0, IKB,IKE,IKL,IIJB,IIJE
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
 !*                    0.2 Initialisation
@@ -82,6 +82,8 @@ IF (LHOOK) CALL DR_HOOK('COMPUTE_MF_CLOUD_DIRECT',0,ZHOOK_HANDLE)
 IKB=D%NKB
 IKE=D%NKE
 IKL=D%NKL
+IIJB=D%NIJB
+IIJE=D%NIJE
 !*      1. COMPUTATION OF SUBGRID CLOUD
 !          ----------------------------
 
@@ -92,7 +94,7 @@ PRC_MF(:,:)=0.
 PRI_MF(:,:)=0.
 PCF_MF(:,:)=0.
 
-DO JI=D%NIJB,D%NIJE
+DO JI=IIJB,IIJE
 #ifdef REPRO48
   JK0=KKLCL(JI)-IKL ! first mass level with cloud
   JK0=MAX(JK0, MIN(IKB,IKE)) !protection if KKL=1
diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index adbb47b8d..9bafea7b1 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -516,19 +516,19 @@ DO JK=IKB,IKE-IKL,IKL
   !$mnh_end_expand_where(JIJ=IIJB:IIJE)
 
   ! If the updraft did not stop, compute cons updraft characteritics at jk+KKL
-  DO JI=D%NIJB,D%NIJE
-    IF(GTEST(JI)) THEN
-      ZMIX2(JI) = (PZZ(JI,JK+IKL)-PZZ(JI,JK))*PENTR(JI,JK) !&
-      ZMIX3_CLD(JI) = (PZZ(JI,JK+IKL)-PZZ(JI,JK))*(1.-ZPART_DRY(JI))*ZDETR_CLD(JI,JK) !&                   
-      ZMIX2_CLD(JI) = (PZZ(JI,JK+IKL)-PZZ(JI,JK))*(1.-ZPART_DRY(JI))*ZENTR_CLD(JI,JK)
+  DO JIJ=IIJB,IIJE
+    IF(GTEST(JIJ)) THEN
+      ZMIX2(JIJ) = (PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*PENTR(JIJ,JK) !&
+      ZMIX3_CLD(JIJ) = (PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*(1.-ZPART_DRY(JIJ))*ZDETR_CLD(JIJ,JK) !&                   
+      ZMIX2_CLD(JIJ) = (PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*(1.-ZPART_DRY(JIJ))*ZENTR_CLD(JIJ,JK)
 #ifdef REPRO48                  
-      PTHL_UP(JI,JK+IKL)=(PTHL_UP(JI,JK)*(1.-0.5*ZMIX2(JI)) + PTHLM(JI,JK)*ZMIX2(JI)) &
-                            /(1.+0.5*ZMIX2(JI))   
-      PRT_UP(JI,JK+IKL) =(PRT_UP (JI,JK)*(1.-0.5*ZMIX2(JI)) + PRTM(JI,JK)*ZMIX2(JI))  &
-                            /(1.+0.5*ZMIX2(JI))
+      PTHL_UP(JIJ,JK+IKL)=(PTHL_UP(JIJ,JK)*(1.-0.5*ZMIX2(JIJ)) + PTHLM(JIJ,JK)*ZMIX2(JIJ)) &
+                            /(1.+0.5*ZMIX2(JIJ))   
+      PRT_UP(JIJ,JK+IKL) =(PRT_UP (JIJ,JK)*(1.-0.5*ZMIX2(JIJ)) + PRTM(JIJ,JK)*ZMIX2(JIJ))  &
+                            /(1.+0.5*ZMIX2(JIJ))
 #else
-      PTHL_UP(JI,JK+IKL)=PTHL_UP(JI,JK)*EXP(-ZMIX2(JI)) + PTHLM(JI,JK)*(1-EXP(-ZMIX2(JI)))
-      PRT_UP(JI,JK+IKL) =PRT_UP (JI,JK)*EXP(-ZMIX2(JI)) +  PRTM(JI,JK)*(1-EXP(-ZMIX2(JI)))
+      PTHL_UP(JIJ,JK+IKL)=PTHL_UP(JIJ,JK)*EXP(-ZMIX2(JIJ)) + PTHLM(JIJ,JK)*(1-EXP(-ZMIX2(JIJ)))
+      PRT_UP(JIJ,JK+IKL) =PRT_UP (JIJ,JK)*EXP(-ZMIX2(JIJ)) +  PRTM(JIJ,JK)*(1-EXP(-ZMIX2(JIJ)))
 #endif
     ENDIF
   ENDDO
@@ -678,8 +678,8 @@ IF(OENTR_DETR) THEN
   ! This way, all MF fluxes are diminished by this amount.
   ! Diagnosed cloud fraction is also multiplied by the same coefficient.
   !
-  DO JI=D%NIJB,D%NIJE
-     PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
+  DO JIJ=IIJB,IIJE
+     PDEPTH(JIJ) = MAX(0., PZZ(JIJ,KKCTL(JIJ)) -  PZZ(JIJ,KKLCL(JIJ)) )
   END DO
 
   !$mnh_expand_array(JIJ=IIJB:IIJE)
@@ -846,35 +846,35 @@ ZPRE(IIJB:IIJE)=PPRE_MINUS_HALF(IIJB:IIJE)
 !$mnh_end_expand_array(JIJ=IIJB:IIJE)
 
 !                1.4 Estimation of PPART_DRY
-DO JI=D%NIJB,D%NIJE
-  IF(OTEST(JI) .AND. OTESTLCL(JI)) THEN
+DO JIJ=IIJB,IIJE
+  IF(OTEST(JIJ) .AND. OTESTLCL(JIJ)) 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
+    PPART_DRY(JIJ)=0.
+    ZDZ_STOP(JIJ)=0.
+    ZPRE(JIJ)=PPRE_MINUS_HALF(JIJ)
+  ELSE IF (OTEST(JIJ) .AND. .NOT. OTESTLCL(JIJ)) THEN
     !Temperature at flux level KK
-    ZT=PTH_UP(JI)*(PPRE_MINUS_HALF(JI)/CST%XP00) ** (CST%XRD/CST%XCPD)
+    ZT=PTH_UP(JIJ)*(PPRE_MINUS_HALF(JIJ)/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))
+    ZFOESW = MIN(EXP( CST%XALPW - CST%XBETAW/ZT - CST%XGAMW*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JIJ))
+    ZFOESI = MIN(EXP( CST%XALPI - CST%XBETAI/ZT - CST%XGAMI*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JIJ))
     !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)))
+    ZDRSATODP=(CST%XBETAW/ZT-CST%XGAMW)*(1-ZFRAC_ICE(JIJ))+(CST%XBETAI/ZT-CST%XGAMI)*ZFRAC_ICE(JIJ)
+    ZDRSATODP=((CST%XRD/CST%XCPD)*ZDRSATODP-1.)*PRSAT_UP(JIJ)/ &
+                &(PPRE_MINUS_HALF(JIJ)-(ZFOESW*(1-ZFRAC_ICE(JIJ)) + ZFOESI*ZFRAC_ICE(JIJ)))
     !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
+    ZPRE(JIJ)=PPRE_MINUS_HALF(JIJ)+(PRT_UP(JIJ)-PRSAT_UP(JIJ))/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))))
+    PPART_DRY(JIJ)=MAX(0., MIN(1., (PPRE_MINUS_HALF(JIJ)-ZPRE(JIJ))/(PPRE_MINUS_HALF(JIJ)-PPRE_PLUS_HALF(JIJ))))
     !Height above flux level KK of the cloudy part
-    ZDZ_STOP(JI) = (PZZ(JI,KK+KKL)-PZZ(JI,KK))*PPART_DRY(JI)
+    ZDZ_STOP(JIJ) = (PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))*PPART_DRY(JIJ)
   ELSE
-    PPART_DRY(JI)=0. ! value does not matter, here
+    PPART_DRY(JIJ)=0. ! value does not matter, here
   END IF
 END DO
 
@@ -897,44 +897,44 @@ ZTHV_PLUS_HALF(IIJB:IIJE)  = PTHVM(IIJB:IIJE,KK) + &
 !                  Integral buoyancy and computation of PENTR and PDETR for dry part
 !               --------------------------------------------------------------------
 
-DO JI=D%NIJB,D%NIJE
-  IF (OTEST(JI) .AND. PPART_DRY(JI)>0.) THEN
+DO JIJ=IIJB,IIJE
+  IF (OTEST(JIJ) .AND. PPART_DRY(JIJ)>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(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) )
+    ZDZ=MIN(ZDZ_STOP(JIJ),(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))*0.5)
+    PBUO_INTEG_DRY(JIJ) = ZG_O_THVREF_ED(JIJ)*ZDZ*&
+                (0.5 * (  - ZCOEFF_MINUS_HALF(JIJ))*ZDZ  &
+                  - ZTHV_MINUS_HALF(JIJ) + PTHV_UP(JIJ) )
 
     !Between mass flux KK and bottom of cloudy part (if above mass flux)
-    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) )
+    ZDZ=MAX(0., ZDZ_STOP(JIJ)-(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))*0.5)
+    PBUO_INTEG_DRY(JIJ) = PBUO_INTEG_DRY(JIJ) + ZG_O_THVREF_ED(JIJ)*ZDZ*&
+                (0.5 * (  - ZCOEFF_PLUS_HALF(JIJ))*ZDZ &
+                  - PTHVM(JIJ,KK) + PTHV_UP(JIJ) )
 
     !Entr//Detr. computation
-    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.
+    IF (PBUO_INTEG_DRY(JIJ)>=0.) THEN
+      PENTR(JIJ) = 0.5/(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)*&
+                 LOG(1.+ (2.*(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)/PW_UP2(JIJ,KK))* &
+                 PBUO_INTEG_DRY(JIJ))
+      PDETR(JIJ) = 0.
     ELSE
-      PENTR(JI) = 0.
-      PDETR(JI) = 0.5/(PARAMMF%XABUO)*&
-                 LOG(1.+ (2.*(PARAMMF%XABUO)/PW_UP2(JI,KK))* &
-                 (-PBUO_INTEG_DRY(JI)))
+      PENTR(JIJ) = 0.
+      PDETR(JIJ) = 0.5/(PARAMMF%XABUO)*&
+                 LOG(1.+ (2.*(PARAMMF%XABUO)/PW_UP2(JIJ,KK))* &
+                 (-PBUO_INTEG_DRY(JIJ)))
     ENDIF
-    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))
+    PENTR(JIJ) = PARAMMF%XENTR_DRY*PENTR(JIJ)/(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))    
+    PDETR(JIJ) = PARAMMF%XDETR_DRY*PDETR(JIJ)/(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))
     !Minimum value of detrainment
-    ZWK0D=PLUP(JI)-0.5*(PZZ(JI,KK)+PZZ(JI,KK+KKL))
+    ZWK0D=PLUP(JIJ)-0.5*(PZZ(JIJ,KK)+PZZ(JIJ,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))
+    PDETR(JIJ) = MAX(PPART_DRY(JIJ)*PARAMMF%XDETR_LUP/ZWK0D, PDETR(JIJ))
   ELSE
     !No dry part, condensation reached (OTESTLCL)
-    PBUO_INTEG_DRY(JI) = 0.
-    PENTR(JI)=0.
-    PDETR(JI)=0.
+    PBUO_INTEG_DRY(JIJ) = 0.
+    PENTR(JIJ)=0.
+    PDETR(JIJ)=0.
   ENDIF
 ENDDO
 
@@ -960,29 +960,29 @@ ZTHV_UP_F2(IIJB:IIJE) = ZTHMIX(IIJB:IIJE)*(1.+ZRVORD*ZRVMIX(IIJB:IIJE))/(1.+PRT_
 !$mnh_end_expand_array(JIJ=IIJB:IIJE)
 
 ! Integral buoyancy for cloudy part
-DO JI=D%NIJB,D%NIJE
-  IF(OTEST(JI) .AND. PPART_DRY(JI)<1.) THEN
+DO JIJ=IIJB,IIJE
+  IF(OTEST(JIJ) .AND. PPART_DRY(JIJ)<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)))
+    ZCOTHVU=(ZTHV_UP_F2(JIJ)-PTHV_UP(JIJ))/((PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))*(1-PPART_DRY(JIJ)))
 
     !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) )
+    ZDZ=MAX(0., 0.5*(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))-ZDZ_STOP(JIJ))
+    PBUO_INTEG_CLD(JIJ) = ZG_O_THVREF_ED(JIJ)*ZDZ*&
+            (0.5*( ZCOTHVU - ZCOEFF_MINUS_HALF(JIJ))*ZDZ &
+              - (PTHVM(JIJ,KK)-ZDZ*ZCOEFF_MINUS_HALF(JIJ)) + PTHV_UP(JIJ) )
 
     !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) )
+    ZDZ=(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))-MAX(ZDZ_STOP(JIJ),0.5*(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK)))
+    PBUO_INTEG_CLD(JIJ) = PBUO_INTEG_CLD(JIJ)+ZG_O_THVREF_ED(JIJ)*ZDZ*&
+                      (0.5*( ZCOTHVU - ZCOEFF_PLUS_HALF(JIJ))*ZDZ&
+              - (PTHVM(JIJ,KK)+(0.5*((PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK)))-ZDZ)*ZCOEFF_PLUS_HALF(JIJ)) +&
+              PTHV_UP(JIJ) )
 
   ELSE
     !No cloudy part
-    PBUO_INTEG_CLD(JI)=0.
+    PBUO_INTEG_CLD(JIJ)=0.
   END IF
 END DO
 
@@ -1002,32 +1002,32 @@ ZKIC_INIT=0.1  ! starting value for critical mixed fraction for CLoudy Part
 
 !   JKLIM computed to avoid KKL(KK-KKL) being < KKL*KKB
 JKLIM=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
-DO JI=D%NIJB,D%NIJE
-  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)
+DO JIJ=IIJB,IIJE
+  IF(OTEST(JIJ) .AND. PPART_DRY(JIJ)>0.5) THEN
+    ZDZ=ZDZ_STOP(JIJ)-0.5*(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))
+    ZTHV(JIJ)= PTHVM(JIJ,KK)+ZCOEFF_PLUS_HALF(JIJ)*ZDZ
+    ZMIXTHL(JIJ) = ZKIC_INIT * &
+               (PTHLM(JIJ,KK)+ZDZ*(PTHLM(JIJ,KK+KKL)-PTHLM(JIJ,KK))/PDZZ(JIJ,KK+KKL)) + &
+               (1. - ZKIC_INIT)*PTHL_UP(JIJ)
+    ZMIXRT(JIJ)  = ZKIC_INIT * &
+               (PRTM(JIJ,KK)+ZDZ*(PRTM(JIJ,KK+KKL)-PRTM(JIJ,KK))/PDZZ(JIJ,KK+KKL)) +   &
+               (1. - ZKIC_INIT)*PRT_UP(JIJ)
+  ELSEIF(OTEST(JIJ)) THEN
+    ZDZ=0.5*(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))-ZDZ_STOP(JIJ)
+    ZTHV(JIJ)= PTHVM(JIJ,KK)-ZCOEFF_MINUS_HALF(JIJ)*ZDZ
+    ZMIXTHL(JIJ) = ZKIC_INIT * &
+               (PTHLM(JIJ,KK)-ZDZ*(PTHLM(JIJ,KK)-PTHLM(JIJ,JKLIM))/PDZZ(JIJ,KK)) + &
+               (1. - ZKIC_INIT)*PTHL_UP(JIJ)
+    ZMIXRT(JIJ)  = ZKIC_INIT * &
+               (PRTM(JIJ,KK)-ZDZ*(PRTM(JIJ,KK)-PRTM(JIJ,JKLIM))/PDZZ(JIJ,KK)) + &
+               (1. - ZKIC_INIT)*PRT_UP(JIJ)
   ELSE
 #ifdef REPRO55
-    ZMIXTHL(JI) = 0.1
+    ZMIXTHL(JIJ) = 0.1
 #else
-    ZMIXTHL(JI) = 300.
+    ZMIXTHL(JIJ) = 300.
 #endif
-    ZMIXRT(JI) = 0.1
+    ZMIXRT(JIJ) = 0.1
   ENDIF
 ENDDO
 CALL TH_R_FROM_THL_RT(CST,NEB,D%NIJT,HFRAC_ICE,ZFRAC_ICE,&
@@ -1054,25 +1054,25 @@ ZTHVMIX_F2(IIJB:IIJE) = ZTHMIX(IIJB:IIJE)*(1.+ZRVORD*ZRVMIX(IIJB:IIJE))/(1.+ZMIX
 !$mnh_end_expand_array(JIJ=IIJB:IIJE)
 
 !Computation of mean ZKIC over the cloudy part
-DO JI=D%NIJB,D%NIJE
-  IF (OTEST(JI)) THEN
+DO JIJ=IIJB,IIJE
+  IF (OTEST(JIJ)) 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(JI)-ZTHVMIX(JI))<1.E-10) THEN
-      ZKIC(JI)=1.
+    IF (ABS(PTHV_UP(JIJ)-ZTHVMIX(JIJ))<1.E-10) THEN
+      ZKIC(JIJ)=1.
     ELSE
-      ZKIC(JI) = MAX(0.,PTHV_UP(JI)-ZTHV(JI))*ZKIC_INIT /  &  
-                 (PTHV_UP(JI)-ZTHVMIX(JI))
+      ZKIC(JIJ) = MAX(0.,PTHV_UP(JIJ)-ZTHV(JIJ))*ZKIC_INIT /  &  
+                 (PTHV_UP(JIJ)-ZTHVMIX(JIJ))
     END IF
     ! Compute ZKIC_F2 at flux level KK+KKL
-    IF (ABS(ZTHV_UP_F2(JI)-ZTHVMIX_F2(JI))<1.E-10) THEN
-      ZKIC_F2(JI)=1.
+    IF (ABS(ZTHV_UP_F2(JIJ)-ZTHVMIX_F2(JIJ))<1.E-10) THEN
+      ZKIC_F2(JIJ)=1.
     ELSE
-      ZKIC_F2(JI) = MAX(0.,ZTHV_UP_F2(JI)-ZTHV_PLUS_HALF(JI))*ZKIC_INIT /  &  
-                 (ZTHV_UP_F2(JI)-ZTHVMIX_F2(JI))
+      ZKIC_F2(JIJ) = MAX(0.,ZTHV_UP_F2(JIJ)-ZTHV_PLUS_HALF(JIJ))*ZKIC_INIT /  &  
+                 (ZTHV_UP_F2(JIJ)-ZTHVMIX_F2(JIJ))
     END IF
     !Mean ZKIC over the cloudy part
-    ZKIC(JI)=MAX(MIN(0.5*(ZKIC(JI)+ZKIC_F2(JI)),1.),0.)
+    ZKIC(JIJ)=MAX(MIN(0.5*(ZKIC(JIJ)+ZKIC_F2(JIJ)),1.),0.)
   END IF
 END DO
 
@@ -1083,10 +1083,10 @@ END DO
 
 !Constant PDF
 !For this PDF, eq. (5) is delta Me=0.5*delta Mt
-DO JI=D%NIJB,D%NIJE
-  IF(OTEST(JI)) THEN
-    ZEPSI(JI) = ZKIC(JI)**2. !integration multiplied by 2
-    ZDELTA(JI) = (1.-ZKIC(JI))**2. !idem
+DO JIJ=IIJB,IIJE
+  IF(OTEST(JIJ)) THEN
+    ZEPSI(JIJ) = ZKIC(JIJ)**2. !integration multiplied by 2
+    ZDELTA(JIJ) = (1.-ZKIC(JIJ))**2. !idem
   ENDIF
 ENDDO
 
@@ -1106,16 +1106,16 @@ ENDDO
 !ENDWHERE
 
 !               3.4 Computation of PENTR and PDETR
-DO JI=D%NIJB,D%NIJE
-  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)
+DO JIJ=IIJB,IIJE
+  IF(OTEST(JIJ)) THEN
+    ZEPSI_CLOUD=MIN(ZDELTA(JIJ), ZEPSI(JIJ))
+    PENTR_CLD(JIJ) = (1.-PPART_DRY(JIJ))*ZCOEFFMF_CLOUD*PRHODREF(JIJ)*ZEPSI_CLOUD
+    PDETR_CLD(JIJ) = (1.-PPART_DRY(JIJ))*ZCOEFFMF_CLOUD*PRHODREF(JIJ)*ZDELTA(JIJ)
+    PENTR(JIJ) = PENTR(JIJ)+PENTR_CLD(JIJ)
+    PDETR(JIJ) = PDETR(JIJ)+PDETR_CLD(JIJ)
   ELSE
-    PENTR_CLD(JI) = 0.
-    PDETR_CLD(JI) = 0.
+    PENTR_CLD(JIJ) = 0.
+    PDETR_CLD(JIJ) = 0.
   ENDIF
 ENDDO
 
diff --git a/src/common/turb/mode_compute_updraft_raha.F90 b/src/common/turb/mode_compute_updraft_raha.F90
index 069068e89..fa7e5a371 100644
--- a/src/common/turb/mode_compute_updraft_raha.F90
+++ b/src/common/turb/mode_compute_updraft_raha.F90
@@ -634,8 +634,8 @@ PEMF(IIJB:IIJE,IKB) =0.
 ! This way, all MF fluxes are diminished by this amount.
 ! Diagnosed cloud fraction is also multiplied by the same coefficient.
 !
-DO JI=D%NIJB,D%NIJE 
-  PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
+DO JIJ=IIJB,IIJE 
+  PDEPTH(JIJ) = MAX(0., PZZ(JIJ,KKCTL(JIJ)) -  PZZ(JIJ,KKLCL(JIJ)) )
 END DO
 
 !$mnh_expand_array(JIJ=IIJB:IIJE)
diff --git a/src/common/turb/mode_compute_updraft_rhcj10.F90 b/src/common/turb/mode_compute_updraft_rhcj10.F90
index c5f8ce0a8..4aabe9609 100644
--- a/src/common/turb/mode_compute_updraft_rhcj10.F90
+++ b/src/common/turb/mode_compute_updraft_rhcj10.F90
@@ -289,13 +289,13 @@ PSV_UP(IIJB:IIJE,1:IKT,:)=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=D%NIJB,D%NIJE
-  !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,IKB)= ZTHLM_F(JI,IKB)
-  PRT_UP(JI,IKB) = ZRTM_F(JI,IKB)
-  !ZQT_UP(JI) = PRT_UP(JI,KKB)/(1.+PRT_UP(JI,KKB))
-  !ZTHS_UP(JI,KKB)=PTHL_UP(JI,KKB)*(1.+XLAMBDA_MF*ZQT_UP(JI))
+DO JIJ=IIJB,IIJE
+  !PTHL_UP(JIJ,KKB)= ZTHLM_F(JIJ,KKB)+MAX(0.,MIN(ZTMAX,(PSFTH(JIJ)/SQRT(ZTKEM_F(JIJ,KKB)))*XALP_PERT))
+  !PRT_UP(JIJ,KKB) = ZRTM_F(JIJ,KKB)+MAX(0.,MIN(ZRMAX,(PSFRV(JIJ)/SQRT(ZTKEM_F(JIJ,KKB)))*XALP_PERT)) 
+  PTHL_UP(JIJ,IKB)= ZTHLM_F(JIJ,IKB)
+  PRT_UP(JIJ,IKB) = ZRTM_F(JIJ,IKB)
+  !ZQT_UP(JIJ) = PRT_UP(JIJ,KKB)/(1.+PRT_UP(JIJ,KKB))
+  !ZTHS_UP(JIJ,KKB)=PTHL_UP(JIJ,KKB)*(1.+XLAMBDA_MF*ZQT_UP(JIJ))
 ENDDO
 
 CALL MZM_MF(D, PTHM (:,:), ZTHM_F(:,:))
@@ -305,8 +305,8 @@ CALL MZM_MF(D, PRVM(:,:), ZRVM_F(:,:))
 
 ! thetav at mass and flux levels 
 DO JK=1,IKT
-  DO JI=D%NIB,D%NIJE
-    ZTHVM_F(JI,JK)=ZTHM_F(JI,JK)*((1.+ZRVORD*ZRVM_F(JI,JK))/(1.+ZRTM_F(JI,JK)))
+  DO JIJ=D%NIB,D%NIJE
+    ZTHVM_F(JIJ,JK)=ZTHM_F(JIJ,JK)*((1.+ZRVORD*ZRVM_F(JIJ,JK))/(1.+ZRTM_F(JIJ,JK)))
   ENDDO
 ENDDO
 
@@ -333,11 +333,11 @@ CALL TH_R_FROM_THL_RT(CST,NEB,D%NIJT,HFRAC_ICE,PFRAC_ICE_UP(:,IKB),ZPRES_F(:,IKB
              PRV_UP(:,IKB),PRC_UP(:,IKB),PRI_UP(:,IKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
              PBUF=ZBUF, KB=D%NIJB, KE=D%NIJE)
 
-DO JI=D%NIJB,D%NIJE
+DO JIJ=IIJB,IIJE
   ! compute updraft thevav and buoyancy term at KKB level             
-  PTHV_UP(JI,IKB) = ZTH_UP(JI,IKB)*((1+ZRVORD*PRV_UP(JI,IKB))/(1+PRT_UP(JI,IKB))) 
+  PTHV_UP(JIJ,IKB) = ZTH_UP(JIJ,IKB)*((1+ZRVORD*PRV_UP(JIJ,IKB))/(1+PRT_UP(JIJ,IKB))) 
   ! compute mean rsat in updraft
-  PRSAT_UP(JI,IKB) = ZRSATW(JI)*(1-PFRAC_ICE_UP(JI,IKB)) + ZRSATI(JI)*PFRAC_ICE_UP(JI,IKB)
+  PRSAT_UP(JIJ,IKB) = ZRSATW(JIJ)*(1-PFRAC_ICE_UP(JIJ,IKB)) + ZRSATI(JIJ)*PFRAC_ICE_UP(JIJ,IKB)
 ENDDO
 
 !Tout est commente pour tester dans un premier temps la separation en deux de la 
@@ -350,8 +350,8 @@ ZG_O_THVREF(IIJB:IIJE,1:IKT)=CST%XG/ZTHVM_F(IIJB:IIJE,1:IKT)
 ! Calcul de la fermeture de Julien Pergaut comme limite max de PHY
 
 DO JK=IKB,IKE-IKL,IKL   !  Vertical loop
-  DO JI=D%NIJB,D%NIJE
-    ZZDZ(JI,JK)   = MAX(ZEPS,PZZ(JI,JK+IKL)-PZZ(JI,JK))  ! <== Delta Z between two flux level
+  DO JIJ=IIJB,IIJE
+    ZZDZ(JIJ,JK)   = MAX(ZEPS,PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))  ! <== Delta Z between two flux level
   ENDDO
 ENDDO
 
@@ -379,21 +379,21 @@ CALL COMPUTE_BL89_ML(D, CST, CSTURB, PDZZ,ZTKEM_F(:,IKB),ZG_O_THVREF(:,IKB), &
 ZLUP(IIJB:IIJE)=MAX(ZLUP(IIJB:IIJE),1.E-10)
 !$mnh_end_expand_array(JIJ=IIJB:IIJE)
 
-DO JI=D%NIJB,D%NIJE
+DO JIJ=IIJB,IIJE
   ! Compute Buoyancy flux at the ground
-  ZWTHVSURF = (ZTHVM_F(JI,IKB)/ZTHM_F(JI,IKB))*PSFTH(JI)+     &
-              (0.61*ZTHM_F(JI,IKB))*PSFRV(JI)
+  ZWTHVSURF = (ZTHVM_F(JIJ,IKB)/ZTHM_F(JIJ,IKB))*PSFTH(JIJ)+     &
+              (0.61*ZTHM_F(JIJ,IKB))*PSFRV(JIJ)
 
   ! Mass flux at KKB level (updraft triggered if PSFTH>0.)
   IF (ZWTHVSURF>0.010) THEN ! <==  Not 0 Important to have stratocumulus !!!!!
-    PEMF(JI,IKB) = PARAMMF%XCMF * ZRHO_F(JI,IKB) * ((ZG_O_THVREF(JI,IKB))*ZWTHVSURF*ZLUP(JI))**(1./3.)
-    PFRAC_UP(JI,IKB)=MIN(PEMF(JI,IKB)/(SQRT(ZW_UP2(JI,IKB))*ZRHO_F(JI,IKB)),PARAMMF%XFRAC_UP_MAX)
-    !PEMF(JI,KKB) = ZRHO_F(JI,KKB)*PFRAC_UP(JI,KKB)*SQRT(ZW_UP2(JI,KKB))
-    ZW_UP2(JI,IKB)=(PEMF(JI,IKB)/(PFRAC_UP(JI,IKB)*ZRHO_F(JI,IKB)))**2
-    GTEST(JI)=.TRUE.
+    PEMF(JIJ,IKB) = PARAMMF%XCMF * ZRHO_F(JIJ,IKB) * ((ZG_O_THVREF(JIJ,IKB))*ZWTHVSURF*ZLUP(JIJ))**(1./3.)
+    PFRAC_UP(JIJ,IKB)=MIN(PEMF(JIJ,IKB)/(SQRT(ZW_UP2(JIJ,IKB))*ZRHO_F(JIJ,IKB)),PARAMMF%XFRAC_UP_MAX)
+    !PEMF(JIJ,KKB) = ZRHO_F(JIJ,KKB)*PFRAC_UP(JIJ,KKB)*SQRT(ZW_UP2(JIJ,KKB))
+    ZW_UP2(JIJ,IKB)=(PEMF(JIJ,IKB)/(PFRAC_UP(JIJ,IKB)*ZRHO_F(JIJ,IKB)))**2
+    GTEST(JIJ)=.TRUE.
   ELSE
-    PEMF(JI,IKB) =0.
-    GTEST(JI)=.FALSE.
+    PEMF(JIJ,IKB) =0.
+    GTEST(JIJ)=.FALSE.
   ENDIF
 ENDDO
 
@@ -427,10 +427,10 @@ DO JK=IKB,IKE-IKL,IKL
  
 ! to find the LCL (check if JK is LCL or not)
 
-  DO JI=D%NIJB,D%NIJE
-    IF ((PRC_UP(JI,JK)+PRI_UP(JI,JK)>0.).AND.(.NOT.(GTESTLCL(JI)))) THEN
-      KKLCL(JI) = JK           
-      GTESTLCL(JI)=.TRUE.
+  DO JIJ=IIJB,IIJE
+    IF ((PRC_UP(JIJ,JK)+PRI_UP(JIJ,JK)>0.).AND.(.NOT.(GTESTLCL(JIJ)))) THEN
+      KKLCL(JIJ) = JK           
+      GTESTLCL(JIJ)=.TRUE.
     ENDIF
   ENDDO
 
@@ -452,83 +452,83 @@ DO JK=IKB,IKE-IKL,IKL
                ZTH_UP(:,JK),ZRV_UP,ZRC_UP,ZRI_UP,ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
                PBUF=ZBUF, KB=D%NIJB, KE=D%NIJE)
     
-  DO JI=D%NIJB,D%NIJE
-    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))    
-      PBUO_INTEG(JI,JK) = ZBUO(JI,JK)*(PZZ(JI,JK+IKL)-PZZ(JI,JK))
+  DO JIJ=IIJB,IIJE
+    IF (GTEST(JIJ)) THEN
+      PTHV_UP(JIJ,JK)    = ZTH_UP(JIJ,JK)*(1.+ZRVORD*ZRV_UP(JIJ))/(1.+PRT_UP(JIJ,JK))
+      ZBUO(JIJ,JK)       = ZG_O_THVREF(JIJ,JK)*(PTHV_UP(JIJ,JK) - ZTHVM_F(JIJ,JK))    
+      PBUO_INTEG(JIJ,JK) = ZBUO(JIJ,JK)*(PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))
       
-      ZDZ(JI)   = MAX(ZEPS,PZZ(JI,JK+IKL)-PZZ(JI,JK))
-      ZTEST(JI) = PARAMMF%XA1*ZBUO(JI,JK) -  PARAMMF%XB*ZW_UP2(JI,JK)
+      ZDZ(JIJ)   = MAX(ZEPS,PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))
+      ZTEST(JIJ) = PARAMMF%XA1*ZBUO(JIJ,JK) -  PARAMMF%XB*ZW_UP2(JIJ,JK)
 
       !  Ancien calcul de la vitesse
-      ZCOE(JI)      = ZDZ(JI)
-      IF (ZTEST(JI)>0.) THEN
-        ZCOE(JI)    = ZDZ(JI)/(1.+ PARAMMF%XBETA1)
+      ZCOE(JIJ)      = ZDZ(JIJ)
+      IF (ZTEST(JIJ)>0.) THEN
+        ZCOE(JIJ)    = ZDZ(JIJ)/(1.+ PARAMMF%XBETA1)
       ENDIF
 
       !  Convective Vertical speed computation
-      ZWCOE(JI)         = (1.-PARAMMF%XB*ZCOE(JI))/(1.+PARAMMF%XB*ZCOE(JI))
-      ZBUCOE(JI)        =  2.*ZCOE(JI)/(1.+PARAMMF%XB*ZCOE(JI))
+      ZWCOE(JIJ)         = (1.-PARAMMF%XB*ZCOE(JIJ))/(1.+PARAMMF%XB*ZCOE(JIJ))
+      ZBUCOE(JIJ)        =  2.*ZCOE(JIJ)/(1.+PARAMMF%XB*ZCOE(JIJ))
 
       ! Second Rachel bug correction (XA1 has been forgotten)
-      ZW_UP2(JI,JK+IKL) = MAX(ZEPS,ZW_UP2(JI,JK)*ZWCOE(JI) + PARAMMF%XA1*ZBUO(JI,JK)*ZBUCOE(JI) )  
-      ZW_MAX(JI) = MAX(ZW_MAX(JI), SQRT(ZW_UP2(JI,JK+IKL)))
-      ZWUP_MEAN(JI)     = MAX(ZEPS,0.5*(ZW_UP2(JI,JK+IKL)+ZW_UP2(JI,JK)))
+      ZW_UP2(JIJ,JK+IKL) = MAX(ZEPS,ZW_UP2(JIJ,JK)*ZWCOE(JIJ) + PARAMMF%XA1*ZBUO(JIJ,JK)*ZBUCOE(JIJ) )  
+      ZW_MAX(JIJ) = MAX(ZW_MAX(JIJ), SQRT(ZW_UP2(JIJ,JK+IKL)))
+      ZWUP_MEAN(JIJ)     = MAX(ZEPS,0.5*(ZW_UP2(JIJ,JK+IKL)+ZW_UP2(JIJ,JK)))
  
       !  Entrainement and detrainement
 
       ! First Rachel bug correction (Parenthesis around 1+beta1 ==> impact is small)   
-      PENTR(JI,JK)  = MAX(0.,(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*(PARAMMF%XA1*ZBUO(JI,JK)/ZWUP_MEAN(JI)-PARAMMF%XB))
-      ZDETR_BUO(JI) = MAX(0., -(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*PARAMMF%XA1*ZBUO(JI,JK)/ZWUP_MEAN(JI))
-      ZDETR_RT(JI)  = PARAMMF%XC*SQRT(MAX(0.,(PRT_UP(JI,JK) - ZRTM_F(JI,JK))) / MAX(ZEPS,ZRTM_F(JI,JK)) / ZWUP_MEAN(JI))
-      PDETR(JI,JK)  = ZDETR_RT(JI)+ZDETR_BUO(JI)
+      PENTR(JIJ,JK)  = MAX(0.,(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*(PARAMMF%XA1*ZBUO(JIJ,JK)/ZWUP_MEAN(JIJ)-PARAMMF%XB))
+      ZDETR_BUO(JIJ) = MAX(0., -(PARAMMF%XBETA1/(1.+PARAMMF%XBETA1))*PARAMMF%XA1*ZBUO(JIJ,JK)/ZWUP_MEAN(JIJ))
+      ZDETR_RT(JIJ)  = PARAMMF%XC*SQRT(MAX(0.,(PRT_UP(JIJ,JK) - ZRTM_F(JIJ,JK))) / MAX(ZEPS,ZRTM_F(JIJ,JK)) / ZWUP_MEAN(JIJ))
+      PDETR(JIJ,JK)  = ZDETR_RT(JIJ)+ZDETR_BUO(JIJ)
    
       ! If the updraft did not stop, compute cons updraft characteritics at jk+1
-      ZZTOP(JI) = MAX(ZZTOP(JI),PZZ(JI,JK+IKL))
-      ZMIX2(JI) = (PZZ(JI,JK+IKL)-PZZ(JI,JK))*PENTR(JI,JK) !&
-
-      !ZQTM(JI) = PRTM(JI,JK)/(1.+PRTM(JI,JK))            
-      !ZTHSM(JI,JK) = PTHLM(JI,JK)*(1.+XLAMBDA_MF*ZQTM(JI))
-      !ZTHS_UP(JI,JK+KKL)=(ZTHS_UP(JI,JK)*(1.-0.5*ZMIX2(JI)) + ZTHSM(JI,JK)*ZMIX2(JI)) &
-      !                     /(1.+0.5*ZMIX2(JI))
-      PRT_UP(JI,JK+IKL) =(PRT_UP (JI,JK)*(1.-0.5*ZMIX2(JI)) + PRTM(JI,JK)*ZMIX2(JI))  &
-                           /(1.+0.5*ZMIX2(JI))
-      !ZQT_UP(JI) = PRT_UP(JI,JK+KKL)/(1.+PRT_UP(JI,JK+KKL))
-      !PTHL_UP(JI,JK+KKL)=ZTHS_UP(JI,JK+KKL)/(1.+XLAMBDA_MF*ZQT_UP(JI))
-      PTHL_UP(JI,JK+IKL)=(PTHL_UP(JI,JK)*(1.-0.5*ZMIX2(JI)) + PTHLM(JI,JK)*ZMIX2(JI)) &
-                           /(1.+0.5*ZMIX2(JI))
+      ZZTOP(JIJ) = MAX(ZZTOP(JIJ),PZZ(JIJ,JK+IKL))
+      ZMIX2(JIJ) = (PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*PENTR(JIJ,JK) !&
+
+      !ZQTM(JIJ) = PRTM(JIJ,JK)/(1.+PRTM(JIJ,JK))            
+      !ZTHSM(JIJ,JK) = PTHLM(JIJ,JK)*(1.+XLAMBDA_MF*ZQTM(JIJ))
+      !ZTHS_UP(JIJ,JK+KKL)=(ZTHS_UP(JIJ,JK)*(1.-0.5*ZMIX2(JIJ)) + ZTHSM(JIJ,JK)*ZMIX2(JIJ)) &
+      !                     /(1.+0.5*ZMIX2(JIJ))
+      PRT_UP(JIJ,JK+IKL) =(PRT_UP (JIJ,JK)*(1.-0.5*ZMIX2(JIJ)) + PRTM(JIJ,JK)*ZMIX2(JIJ))  &
+                           /(1.+0.5*ZMIX2(JIJ))
+      !ZQT_UP(JIJ) = PRT_UP(JIJ,JK+KKL)/(1.+PRT_UP(JIJ,JK+KKL))
+      !PTHL_UP(JIJ,JK+KKL)=ZTHS_UP(JIJ,JK+KKL)/(1.+XLAMBDA_MF*ZQT_UP(JIJ))
+      PTHL_UP(JIJ,JK+IKL)=(PTHL_UP(JIJ,JK)*(1.-0.5*ZMIX2(JIJ)) + PTHLM(JIJ,JK)*ZMIX2(JIJ)) &
+                           /(1.+0.5*ZMIX2(JIJ))
     ENDIF  ! GTEST
   ENDDO
   
 
   IF(OMIXUV) THEN
     IF(JK/=IKB) THEN
-      DO JI=D%NIJB,D%NIJE
-        IF(GTEST(JI)) THEN
-          PU_UP(JI,JK+IKL) = (PU_UP (JI,JK)*(1-0.5*ZMIX2(JI)) + PUM(JI,JK)*ZMIX2(JI)+ &
-                            0.5*PARAMMF%XPRES_UV*(PZZ(JI,JK+IKL)-PZZ(JI,JK))*&
-                            ((PUM(JI,JK+IKL)-PUM(JI,JK))/PDZZ(JI,JK+IKL)+&
-                             (PUM(JI,JK)-PUM(JI,JK-IKL))/PDZZ(JI,JK))        )   &
-                            /(1+0.5*ZMIX2(JI))
-          PV_UP(JI,JK+IKL) = (PV_UP (JI,JK)*(1-0.5*ZMIX2(JI)) + PVM(JI,JK)*ZMIX2(JI)+ &
-                            0.5*PARAMMF%XPRES_UV*(PZZ(JI,JK+IKL)-PZZ(JI,JK))*&
-                            ((PVM(JI,JK+IKL)-PVM(JI,JK))/PDZZ(JI,JK+IKL)+&
-                             (PVM(JI,JK)-PVM(JI,JK-IKL))/PDZZ(JI,JK))    )   &
-                            /(1+0.5*ZMIX2(JI))
+      DO JIJ=IIJB,IIJE
+        IF(GTEST(JIJ)) THEN
+          PU_UP(JIJ,JK+IKL) = (PU_UP (JIJ,JK)*(1-0.5*ZMIX2(JIJ)) + PUM(JIJ,JK)*ZMIX2(JIJ)+ &
+                            0.5*PARAMMF%XPRES_UV*(PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*&
+                            ((PUM(JIJ,JK+IKL)-PUM(JIJ,JK))/PDZZ(JIJ,JK+IKL)+&
+                             (PUM(JIJ,JK)-PUM(JIJ,JK-IKL))/PDZZ(JIJ,JK))        )   &
+                            /(1+0.5*ZMIX2(JIJ))
+          PV_UP(JIJ,JK+IKL) = (PV_UP (JIJ,JK)*(1-0.5*ZMIX2(JIJ)) + PVM(JIJ,JK)*ZMIX2(JIJ)+ &
+                            0.5*PARAMMF%XPRES_UV*(PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*&
+                            ((PVM(JIJ,JK+IKL)-PVM(JIJ,JK))/PDZZ(JIJ,JK+IKL)+&
+                             (PVM(JIJ,JK)-PVM(JIJ,JK-IKL))/PDZZ(JIJ,JK))    )   &
+                            /(1+0.5*ZMIX2(JIJ))
         ENDIF
       ENDDO
     ELSE
-      DO JI=D%NIJB,D%NIJE
-        IF(GTEST(JI)) THEN
-          PU_UP(JI,JK+IKL) = (PU_UP (JI,JK)*(1-0.5*ZMIX2(JI)) + PUM(JI,JK)*ZMIX2(JI)+ &
-                            0.5*PARAMMF%XPRES_UV*(PZZ(JI,JK+IKL)-PZZ(JI,JK))*&
-                            ((PUM(JI,JK+IKL)-PUM(JI,JK))/PDZZ(JI,JK+IKL))        ) &
-                            /(1+0.5*ZMIX2(JI))
-          PV_UP(JI,JK+IKL) = (PV_UP (JI,JK)*(1-0.5*ZMIX2(JI)) + PVM(JI,JK)*ZMIX2(JI)+ &
-                            0.5*PARAMMF%XPRES_UV*(PZZ(JI,JK+IKL)-PZZ(JI,JK))*&
-                            ((PVM(JI,JK+IKL)-PVM(JI,JK))/PDZZ(JI,JK+IKL))    )   &
-                            /(1+0.5*ZMIX2(JI))
+      DO JIJ=IIJB,IIJE
+        IF(GTEST(JIJ)) THEN
+          PU_UP(JIJ,JK+IKL) = (PU_UP (JIJ,JK)*(1-0.5*ZMIX2(JIJ)) + PUM(JIJ,JK)*ZMIX2(JIJ)+ &
+                            0.5*PARAMMF%XPRES_UV*(PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*&
+                            ((PUM(JIJ,JK+IKL)-PUM(JIJ,JK))/PDZZ(JIJ,JK+IKL))        ) &
+                            /(1+0.5*ZMIX2(JIJ))
+          PV_UP(JIJ,JK+IKL) = (PV_UP (JIJ,JK)*(1-0.5*ZMIX2(JIJ)) + PVM(JIJ,JK)*ZMIX2(JIJ)+ &
+                            0.5*PARAMMF%XPRES_UV*(PZZ(JIJ,JK+IKL)-PZZ(JIJ,JK))*&
+                            ((PVM(JIJ,JK+IKL)-PVM(JIJ,JK))/PDZZ(JIJ,JK+IKL))    )   &
+                            /(1+0.5*ZMIX2(JIJ))
         ENDIF
       ENDDO
     ENDIF
@@ -555,66 +555,66 @@ DO JK=IKB,IKE-IKL,IKL
           ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
           PBUF=ZBUF, KB=D%NIJB, KE=D%NIJE)
 
-  DO JI=D%NIJB,D%NIJE
-    IF(GTEST(JI)) THEN
-      !ZT_UP(JI) = ZTH_UP(JI,JK+KKL)*PEXNM(JI,JK+KKL)
-      !ZCP(JI) = XCPD + XCL * ZRC_UP(JI)
-      !ZLVOCPEXN(JI)=(XLVTT + (XCPV-XCL) *  (ZT_UP(JI)-XTT) ) / ZCP(JI) / PEXNM(JI,JK+KKL)
-      !PRC_UP(JI,JK+KKL)=MIN(0.5E-3,ZRC_UP(JI))  ! On ne peut depasser 0.5 g/kg (autoconversion donc elimination !)
-      !PTHL_UP(JI,JK+KKL) = PTHL_UP(JI,JK+KKL)+ZLVOCPEXN(JI)*(ZRC_UP(JI)-PRC_UP(JI,JK+KKL))
-      PRC_UP(JI,JK+IKL)=ZRC_UP(JI)
-      PRV_UP(JI,JK+IKL)=ZRV_UP(JI)
-      PRI_UP(JI,JK+IKL)=ZRI_UP(JI)
-      !PRT_UP(JI,JK+KKL)  = PRC_UP(JI,JK+KKL) + PRV_UP(JI,JK+KKL)
-      PRSAT_UP(JI,JK+IKL) = ZRSATW(JI)*(1-PFRAC_ICE_UP(JI,JK+IKL)) + ZRSATI(JI)*PFRAC_ICE_UP(JI,JK+IKL)
+  DO JIJ=IIJB,IIJE
+    IF(GTEST(JIJ)) THEN
+      !ZT_UP(JIJ) = ZTH_UP(JIJ,JK+KKL)*PEXNM(JIJ,JK+KKL)
+      !ZCP(JIJ) = XCPD + XCL * ZRC_UP(JIJ)
+      !ZLVOCPEXN(JIJ)=(XLVTT + (XCPV-XCL) *  (ZT_UP(JIJ)-XTT) ) / ZCP(JIJ) / PEXNM(JIJ,JK+KKL)
+      !PRC_UP(JIJ,JK+KKL)=MIN(0.5E-3,ZRC_UP(JIJ))  ! On ne peut depasser 0.5 g/kg (autoconversion donc elimination !)
+      !PTHL_UP(JIJ,JK+KKL) = PTHL_UP(JIJ,JK+KKL)+ZLVOCPEXN(JIJ)*(ZRC_UP(JIJ)-PRC_UP(JIJ,JK+KKL))
+      PRC_UP(JIJ,JK+IKL)=ZRC_UP(JIJ)
+      PRV_UP(JIJ,JK+IKL)=ZRV_UP(JIJ)
+      PRI_UP(JIJ,JK+IKL)=ZRI_UP(JIJ)
+      !PRT_UP(JIJ,JK+KKL)  = PRC_UP(JIJ,JK+KKL) + PRV_UP(JIJ,JK+KKL)
+      PRSAT_UP(JIJ,JK+IKL) = ZRSATW(JIJ)*(1-PFRAC_ICE_UP(JIJ,JK+IKL)) + ZRSATI(JIJ)*PFRAC_ICE_UP(JIJ,JK+IKL)
 
       ! Compute the updraft theta_v, buoyancy and w**2 for level JK+1   
       !PTHV_UP(IIJB:IIJE,JK+KKL) = PTH_UP(IIJB:IIJE,JK+KKL)*((1+ZRVORD*PRV_UP(IIJB:IIJE,JK+KKL))/(1+PRT_UP(IIJB:IIJE,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))
+      !PTHV_UP(JIJ,JK+KKL) = ZTH_UP(JIJ,JK+KKL)*(1.+0.608*PRV_UP(JIJ,JK+KKL) - PRC_UP(JIJ,JK+KKL))
       !! A corriger pour utiliser q et non r !!!!      
-      !ZMIX1(JI)=ZZDZ(JI,JK)*(PENTR(JI,JK)-PDETR(JI,JK))
-      PTHV_UP(JI,JK+IKL) = ZTH_UP(JI,JK+IKL)*((1+ZRVORD*PRV_UP(JI,JK+IKL))/(1+PRT_UP(JI,JK+IKL)))
-      ZMIX1(JI)=ZZDZ(JI,JK)*(PENTR(JI,JK)-PDETR(JI,JK))
+      !ZMIX1(JIJ)=ZZDZ(JIJ,JK)*(PENTR(JIJ,JK)-PDETR(JIJ,JK))
+      PTHV_UP(JIJ,JK+IKL) = ZTH_UP(JIJ,JK+IKL)*((1+ZRVORD*PRV_UP(JIJ,JK+IKL))/(1+PRT_UP(JIJ,JK+IKL)))
+      ZMIX1(JIJ)=ZZDZ(JIJ,JK)*(PENTR(JIJ,JK)-PDETR(JIJ,JK))
     ENDIF
   ENDDO
 
-  DO JI=D%NIJB,D%NIJE
-    IF(GTEST(JI)) THEN
-      PEMF(JI,JK+IKL)=PEMF(JI,JK)*EXP(ZMIX1(JI))
+  DO JIJ=IIJB,IIJE
+    IF(GTEST(JIJ)) THEN
+      PEMF(JIJ,JK+IKL)=PEMF(JIJ,JK)*EXP(ZMIX1(JIJ))
     ENDIF
   ENDDO
 
-  DO JI=D%NIJB,D%NIJE
-    IF(GTEST(JI)) THEN
+  DO JIJ=IIJB,IIJE
+    IF(GTEST(JIJ)) THEN
       ! Updraft fraction must be smaller than XFRAC_UP_MAX
-      PFRAC_UP(JI,JK+IKL)=MIN(PARAMMF%XFRAC_UP_MAX, &
-                             &PEMF(JI,JK+IKL)/(SQRT(ZW_UP2(JI,JK+IKL))*ZRHO_F(JI,JK+IKL)))
-      !PEMF(JI,JK+KKL) = ZRHO_F(JI,JK+KKL)*PFRAC_UP(JI,JK+KKL)*SQRT(ZW_UP2(JI,JK+KKL))
+      PFRAC_UP(JIJ,JK+IKL)=MIN(PARAMMF%XFRAC_UP_MAX, &
+                             &PEMF(JIJ,JK+IKL)/(SQRT(ZW_UP2(JIJ,JK+IKL))*ZRHO_F(JIJ,JK+IKL)))
+      !PEMF(JIJ,JK+KKL) = ZRHO_F(JIJ,JK+KKL)*PFRAC_UP(JIJ,JK+KKL)*SQRT(ZW_UP2(JIJ,JK+KKL))
     ENDIF
   ENDDO
 
 ! Test if the updraft has reach the ETL
-  DO JI=D%NIJB,D%NIJE
-    IF (GTEST(JI) .AND. (PBUO_INTEG(JI,JK)<=0.)) THEN
-      KKETL(JI) = JK+IKL
+  DO JIJ=IIJB,IIJE
+    IF (GTEST(JIJ) .AND. (PBUO_INTEG(JIJ,JK)<=0.)) THEN
+      KKETL(JIJ) = JK+IKL
     ENDIF
   ENDDO
 
 
 ! Test is we have reached the top of the updraft
-  DO JI=D%NIJB,D%NIJE
-    IF (GTEST(JI) .AND. ((ZW_UP2(JI,JK+IKL)<=ZEPS).OR.(PEMF(JI,JK+IKL)<=ZEPS))) THEN
-      ZW_UP2   (JI,JK+IKL)=ZEPS
-      PEMF     (JI,JK+IKL)=0.
-      GTEST    (JI)       =.FALSE.
-      PTHL_UP  (JI,JK+IKL)=ZTHLM_F(JI,JK+IKL)
-      PRT_UP   (JI,JK+IKL)=ZRTM_F(JI,JK+IKL)
-      PRC_UP   (JI,JK+IKL)=0.
-      PRI_UP   (JI,JK+IKL)=0.
-      PRV_UP   (JI,JK+IKL)=ZRVM_F (JI,JK+IKL)
-      PTHV_UP  (JI,JK+IKL)=ZTHVM_F(JI,JK+IKL)
-      PFRAC_UP (JI,JK+IKL)=0.
-      KKCTL    (JI)       =JK+IKL
+  DO JIJ=IIJB,IIJE
+    IF (GTEST(JIJ) .AND. ((ZW_UP2(JIJ,JK+IKL)<=ZEPS).OR.(PEMF(JIJ,JK+IKL)<=ZEPS))) THEN
+      ZW_UP2   (JIJ,JK+IKL)=ZEPS
+      PEMF     (JIJ,JK+IKL)=0.
+      GTEST    (JIJ)       =.FALSE.
+      PTHL_UP  (JIJ,JK+IKL)=ZTHLM_F(JIJ,JK+IKL)
+      PRT_UP   (JIJ,JK+IKL)=ZRTM_F(JIJ,JK+IKL)
+      PRC_UP   (JIJ,JK+IKL)=0.
+      PRI_UP   (JIJ,JK+IKL)=0.
+      PRV_UP   (JIJ,JK+IKL)=ZRVM_F (JIJ,JK+IKL)
+      PTHV_UP  (JIJ,JK+IKL)=ZTHVM_F(JIJ,JK+IKL)
+      PFRAC_UP (JIJ,JK+IKL)=0.
+      KKCTL    (JIJ)       =JK+IKL
     ENDIF
   ENDDO
 
@@ -633,8 +633,8 @@ PEMF(IIJB:IIJE,IKB) =0.
 ! This way, all MF fluxes are diminished by this amount.
 ! Diagnosed cloud fraction is also multiplied by the same coefficient.
 !
-DO JI=D%NIJB,D%NIJE
-   PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
+DO JIJ=IIJB,IIJE
+   PDEPTH(JIJ) = MAX(0., PZZ(JIJ,KKCTL(JIJ)) -  PZZ(JIJ,KKLCL(JIJ)) )
 ENDDO
 
 !$mnh_expand_array(JIJ=IIJB:IIJE)
@@ -648,10 +648,10 @@ DO JK=1,IKT
   !$mnh_end_expand_array(JIJ=IIJB:IIJE)
 ENDDO
 DO JK=1,IKT
-  DO JI=D%NIJB,D%NIJE
-    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) 
+  DO JIJ=IIJB,IIJE
+    IF (GWORK2(JIJ,JK)) THEN
+      PEMF(JIJ,JK)     = PEMF(JIJ,JK)     * ZCOEF(JIJ,JK) 
+      PFRAC_UP(JIJ,JK) = PFRAC_UP(JIJ,JK) * ZCOEF(JIJ,JK) 
     ENDIF
   ENDDO
 ENDDO
diff --git a/src/common/turb/mode_mf_turb_expl.F90 b/src/common/turb/mode_mf_turb_expl.F90
index 4c65ba81b..d7b3eb76a 100644
--- a/src/common/turb/mode_mf_turb_expl.F90
+++ b/src/common/turb/mode_mf_turb_expl.F90
@@ -177,20 +177,20 @@ ENDIF
 !          --------------------------------------------
 
 DO JK=IKB,IKE-IKL,IKL
-  DO JI=D%NIJB,D%NIJE
-    !PTHLDT(JI,JK) = (PFLXZTHLMF(JI,JK  ) - PFLXZTHLMF(JI,JK+IKL)) / PRHODJ(JI,JK)
-    PRTDT(JI,JK) = (PFLXZRMF(JI,JK) - PFLXZRMF(JI,JK+IKL)) / PRHODJ(JI,JK)
-    ZQTDT(JI,JK) = PRTDT(JI,JK)/(1.+ ZRTM_F(JI,JK)*ZRTM_F(JI,JK))
-    ZTHSDT(JI,JK)= (ZFLXZTHSMF(JI,JK) - ZFLXZTHSMF(JI,JK+IKL)) / PRHODJ(JI,JK)
-    PTHLDT(JI,JK) = ZTHSDT(JI,JK)/(1.+PARAMMF%XLAMBDA_MF*ZQTM(JI,JK)) - ZTHLM_F(JI,JK)*PARAMMF%XLAMBDA_MF*ZQTDT(JI,JK)
+  DO JIJ=IIJB,IIJE
+    !PTHLDT(JIJ,JK) = (PFLXZTHLMF(JIJ,JK  ) - PFLXZTHLMF(JIJ,JK+IKL)) / PRHODJ(JIJ,JK)
+    PRTDT(JIJ,JK) = (PFLXZRMF(JIJ,JK) - PFLXZRMF(JIJ,JK+IKL)) / PRHODJ(JIJ,JK)
+    ZQTDT(JIJ,JK) = PRTDT(JIJ,JK)/(1.+ ZRTM_F(JIJ,JK)*ZRTM_F(JIJ,JK))
+    ZTHSDT(JIJ,JK)= (ZFLXZTHSMF(JIJ,JK) - ZFLXZTHSMF(JIJ,JK+IKL)) / PRHODJ(JIJ,JK)
+    PTHLDT(JIJ,JK) = ZTHSDT(JIJ,JK)/(1.+PARAMMF%XLAMBDA_MF*ZQTM(JIJ,JK)) - ZTHLM_F(JIJ,JK)*PARAMMF%XLAMBDA_MF*ZQTDT(JIJ,JK)
   ENDDO
 END DO
 
 IF (OMIXUV) THEN
   DO JK=IKB,IKE-IKL,IKL
-    DO JI=D%NIJB,D%NIJE
-      PUDT(JI,JK) = (PFLXZUMF(JI,JK) - PFLXZUMF(JI,JK+IKL)) / PRHODJ(JI,JK)
-      PVDT(JI,JK) = (PFLXZVMF(JI,JK) - PFLXZVMF(JI,JK+IKL)) / PRHODJ(JI,JK)
+    DO JIJ=IIJB,IIJE
+      PUDT(JIJ,JK) = (PFLXZUMF(JIJ,JK) - PFLXZUMF(JIJ,JK+IKL)) / PRHODJ(JIJ,JK)
+      PVDT(JIJ,JK) = (PFLXZVMF(JIJ,JK) - PFLXZVMF(JIJ,JK+IKL)) / PRHODJ(JIJ,JK)
     ENDDO
   END DO
 ENDIF  
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index ae2612035..120b784a5 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -738,7 +738,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN)    :: PF_LIM  ! Value of F when Phi3 i
 REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PF      ! function F to smooth
 !
 REAL, DIMENSION(D%NIJT,D%NKT) :: ZCOEF
-INTEGER :: JIJ,JK, IIJB,IIJE
+INTEGER :: JIJ,JK, IIJB,IIJE,IKT
 !
 !* adds a artificial correction to smooth the function near the discontinuity
 !  point at Phi3 = Phi_lim
@@ -747,8 +747,7 @@ INTEGER :: JIJ,JK, IIJB,IIJE
 !
 IIJE=D%NIJE
 IIJB=D%NIJB
-
-
+IKT=D%NKT
 !
 !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)    
 ZCOEF(IIJB:IIJE,1:IKT) = MAX(MIN((  10.*(1.-PPHI3(IIJB:IIJE,1:IKT)/CSTURB%XPHI_LIM)) ,1.), 0.) 
@@ -845,7 +844,7 @@ SUBROUTINE PSI_SV(D,CSTURB,KSV,PREDTH1,PREDR1,PREDS1,PRED2THS,PRED2RS,PPHI3,PPSI
   REAL, DIMENSION(D%NIJT,D%NKT),   INTENT(IN) :: PPSI3
   REAL, DIMENSION(D%NIJT,D%NKT,KSV),INTENT(OUT) :: PPSI_SV
 !
-  INTEGER :: IKB, IKE, IIJB,IIJE
+  INTEGER :: IKB, IKE, IIJB,IIJE, IKT
   INTEGER :: JSV,JIJ,JK
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -956,7 +955,7 @@ IF (HTURBDIM=='3DIM') THEN
   END IF
 ELSE
         !* 1DIM case
-DO JK=1,D%IKT
+DO JK=1,IKT
     DO JIJ=IIJB,IIJE
       IF ( ABS(PPHI3(JIJ,JK)-CSTURB%XPHI_LIM) < 1.E-12 ) THEN
          PD_PHI3DTDZ_O_DDTDZ(JIJ,JK)=PPHI3(JIJ,JK)*&
diff --git a/src/common/turb/mode_turb_ver_thermo_corr.F90 b/src/common/turb/mode_turb_ver_thermo_corr.F90
index cdf672b16..6d71d9196 100644
--- a/src/common/turb/mode_turb_ver_thermo_corr.F90
+++ b/src/common/turb/mode_turb_ver_thermo_corr.F90
@@ -326,13 +326,13 @@ REAL, DIMENSION(D%NIJT,D%NKT)  ::  &
        ZWKPHIPSI1,ZWKPHIPSI2,&
        ZWKPHIPSI3,ZWKPHIPSI4       ! working var. for shuman operators (array syntax)
 
-INTEGER             :: IIJB, IIJE, IKB,IKE,IKT,IKA,IKU ! index value for the mass points of the domain 
+INTEGER             :: IIJB, IIJE, IKB,IKE,IKT,IKA ! index value for the mass points of the domain 
 INTEGER             :: IKU  ! array sizes
 INTEGER             :: IKL
 INTEGER             :: JIJ, JK ! loop indexes 
 
-REAL, DIMENSION(D%NIJT,MIN(IKA+JPVEXT_TURB*IKL,IKA+JPVEXT_TURB*IKL+2*IKL):&
-                            MAX(IKA+JPVEXT_TURB*IKL,IKA+JPVEXT_TURB*IKL+2*IKL))&
+REAL, DIMENSION(D%NIJT,MIN(D%NKA+JPVEXT_TURB*D%NKL,D%NKA+JPVEXT_TURB*D%NKL+2*D%NKL):&
+                            MAX(D%NKA+JPVEXT_TURB*D%NKL,D%NKA+JPVEXT_TURB*D%NKL+2*D%NKL))&
                     :: ZCOEFF
                                     ! coefficients for the uncentred gradient
                                     ! computation near the ground, defined in
diff --git a/src/common/turb/shuman_mf.F90 b/src/common/turb/shuman_mf.F90
index b2c36ff65..2c6dabccd 100644
--- a/src/common/turb/shuman_mf.F90
+++ b/src/common/turb/shuman_mf.F90
@@ -210,7 +210,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: PMZM   ! result at flux localizati
 !
 INTEGER :: JK, JIJ
 INTEGER :: IIJB,IIJE ! physical horizontal domain indices
-INTEGER :: IKA,IKU,IKT
+INTEGER :: IKA,IKU,IKT,IKL
 !
 !
 IIJE=D%NIJE
@@ -218,6 +218,7 @@ IIJB=D%NIJB
 IKA=D%NKA
 IKU=D%NKU
 IKT=D%NKT
+IKL=D%NKL
 !
 !-------------------------------------------------------------------------------
 !
@@ -300,7 +301,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: PDZF   ! result at mass
 !
 INTEGER :: JK, JIJ
 INTEGER :: IIJB,IIJE ! physical horizontal domain indices
-INTEGER :: IKA,IKU,IKT
+INTEGER :: IKA,IKU,IKT,IKL
 !
 !-------------------------------------------------------------------------------
 !
@@ -309,6 +310,7 @@ IIJB=D%NIJB
 IKA=D%NKA
 IKU=D%NKU
 IKT=D%NKT
+IKL=D%NKL
 !
 !*       1.    DEFINITION OF DZF
 !              ------------------
@@ -389,7 +391,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: PDZM   ! result at flux
 !
 INTEGER :: JK, JIJ
 INTEGER :: IIJB,IIJE ! physical horizontal domain indices
-INTEGER :: IKA,IKU,IKT
+INTEGER :: IKA,IKU,IKT,IKL
 !
 !-------------------------------------------------------------------------------
 !
@@ -398,6 +400,7 @@ IIJB=D%NIJB
 IKA=D%NKA
 IKU=D%NKU
 IKT=D%NKT
+IKL=D%NKL
 !
 !*       1.    DEFINITION OF DZM
 !              ------------------
@@ -479,7 +482,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: PGZ_M_W  ! result at flux side
 !
 INTEGER  JK, JIJ
 INTEGER :: IIJB,IIJE ! physical horizontal domain indices
-INTEGER :: IKA,IKU,IKT
+INTEGER :: IKA,IKU,IKT,IKL
 !-------------------------------------------------------------------------------
 !
 IIJE=D%NIJE
@@ -487,6 +490,7 @@ IIJB=D%NIJB
 IKA=D%NKA
 IKU=D%NKU
 IKT=D%NKT
+IKL=D%NKL
 !
 !*       1.    COMPUTE THE GRADIENT ALONG Z
 !              -----------------------------
-- 
GitLab