From f7e86c47085bf226ca40dfe3ad483fd5a5ddb6bb Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 29 Mar 2022 15:17:09 +0200
Subject: [PATCH] Quentin 29/03/2022: clean, remove KKA,KKU,KKL args in turb.
 Replaced by the type D%NKA/NKU/NKL

---
 src/common/turb/mode_bl89.F90                 |   5 +-
 src/common/turb/mode_prandtl.F90              |  12 +-
 src/common/turb/mode_tke_eps_sources.F90      |  28 +-
 src/common/turb/mode_turb_ver.F90             |  21 +-
 src/common/turb/mode_turb_ver_dyn_flux.F90    | 199 +++++++-------
 src/common/turb/mode_turb_ver_sv_corr.F90     |  30 +--
 src/common/turb/mode_turb_ver_sv_flux.F90     |  41 ++-
 src/common/turb/mode_turb_ver_thermo_corr.F90 | 245 +++++++++---------
 src/common/turb/mode_turb_ver_thermo_flux.F90 | 206 ++++++++-------
 src/common/turb/turb.F90                      |  15 +-
 10 files changed, 380 insertions(+), 422 deletions(-)

diff --git a/src/common/turb/mode_bl89.F90 b/src/common/turb/mode_bl89.F90
index a97742d72..71f0eab13 100644
--- a/src/common/turb/mode_bl89.F90
+++ b/src/common/turb/mode_bl89.F90
@@ -6,7 +6,7 @@ MODULE MODE_BL89
 IMPLICIT NONE
 CONTAINS
 !     ######spl
-      SUBROUTINE BL89(D,CST,CSTURB,KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM,OOCEAN,HPROGRAM)
+      SUBROUTINE BL89(D,CST,CSTURB,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM,OOCEAN,HPROGRAM)
       USE PARKIND1, ONLY : JPRB
       USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 !     #########################################################
@@ -71,9 +71,6 @@ IMPLICIT NONE
 TYPE(DIMPHYEX_t),         INTENT(IN)  :: D
 TYPE(CST_t),              INTENT(IN)  :: CST
 TYPE(CSTURB_t),           INTENT(IN)  :: CSTURB
-INTEGER,                  INTENT(IN)  :: KKA      !near ground array index
-INTEGER,                  INTENT(IN)  :: KKU      !uppest atmosphere array index
-INTEGER,                  INTENT(IN)  :: KKL      !vert. levels type 1=MNH -1=ARO
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT),   INTENT(IN)  :: PZZ
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT),   INTENT(IN)  :: PDZZ
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT),   INTENT(IN)  :: PTHVREF
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index c8492b523..4a6199e29 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -162,8 +162,8 @@ IMPLICIT NONE
 !*      0.1  declarations of arguments
 !
 TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
-TYPE(CST_t),                  INTENT(IN)    :: CST
-TYPE(CSTURB_t),                  INTENT(IN)    :: CSTURB
+TYPE(CST_t),            INTENT(IN)    :: CST
+TYPE(CSTURB_t),         INTENT(IN)    :: CSTURB
 INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
 INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
 INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
@@ -224,7 +224,6 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::  &
 !                                                     
 INTEGER :: IKB      ! vertical index value for the first inner mass point
 INTEGER :: IKE      ! vertical index value for the last inner mass point
-INTEGER::  ISV                      ! number of scalar variables       
 INTEGER::  JSV                      ! loop index for the scalar variables  
 
 INTEGER :: JLOOP
@@ -252,7 +251,6 @@ ENDIF
 !
 IKB = D%NKB
 IKE = D%NKE 
-ISV  =SIZE(PSVM,4)
 !
 PETHETA(:,:,:) = MZM(ETHETA(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM,OOCEAN,OCOMPUTE_SRC), KKA, KKU, KKL)
 PEMOIST(:,:,:) = MZM(EMOIST(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM,OOCEAN), KKA, KKU, KKL)
@@ -321,11 +319,11 @@ END IF
 !---------------------------------------------------------------------------
 !
 !          For the scalar variables
-DO JSV=1,ISV
+DO JSV=1,KSV
   PREDS1(:,:,:,JSV)=CSTURB%XCTV*PBLL_O_E(:,:,:)*GZ_M_W(KKA, KKU, KKL,PSVM(:,:,:,JSV),PDZZ)
 END DO
 !
-DO JSV=1,ISV
+DO JSV=1,KSV
   ZW2=SIGN(1.,PREDS1(:,:,:,JSV))
   PREDS1(:,:,:,JSV)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDS1(:,:,:,JSV))
 END DO
@@ -412,7 +410,7 @@ END IF   ! end of the if structure on the turbulence dimensionnality
 !
 !           5. Prandtl numbers for scalars
 !              ---------------------------
-DO JSV=1,ISV
+DO JSV=1,KSV
 !
   IF(HTURBDIM=='1DIM') THEN
 !        1D case
diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index 9f981d80e..0478fea7f 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -6,7 +6,7 @@ MODULE MODE_TKE_EPS_SOURCES
 IMPLICIT NONE
 CONTAINS
       SUBROUTINE TKE_EPS_SOURCES(D,CST,CSTURB,BUCONF,HPROGRAM,         &
-                    & KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP,             &
+                    & KMI,PTKEM,PLM,PLEPS,PDP,                         &
                     & PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,        &
                     & PTSTEP,PIMPL,PEXPL,                              &
                     & HTURBLEN,HTURBDIM,                               &
@@ -165,10 +165,6 @@ TYPE(DIMPHYEX_t),        INTENT(IN)   :: D
 TYPE(CST_t),             INTENT(IN)   :: CST
 TYPE(CSTURB_t),          INTENT(IN)   :: CSTURB
 TYPE(TBUDGETCONF_t),     INTENT(IN)   :: BUCONF
-INTEGER,                 INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                 INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                 INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
-
 INTEGER,                 INTENT(IN)   ::  KMI          ! model index number  
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)   ::  PTKEM        ! TKE at t-deltat
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)   ::  PLM          ! mixing length         
@@ -234,7 +230,6 @@ NULLIFY(TZFIELDDISS_ll)
 !*       1.   PRELIMINARY COMPUTATIONS
 !             ------------------------
 !
-
 IF (LHOOK) CALL DR_HOOK('TKE_EPS_SOURCES',0,ZHOOK_HANDLE)
 !
 IKB=D%NKB
@@ -259,12 +254,11 @@ ELSE
 END IF
 !
 !
-!
 !*       2.2  Explicit TKE sources except horizontal turbulent transport 
 !
 ! extrapolate the dynamic production with a 1/Z law from its value at the 
 ! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB)
-PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+KKL)/PDZZ(:,:,IKB))
+PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+D%NKL)/PDZZ(:,:,IKB))
 !
 ! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
 ! + (Dynamical Production) + (Thermal Production) - (dissipation) 
@@ -280,11 +274,11 @@ ZSOURCE(:,:,:) = ( PRTKES(:,:,:) +  PRTKEMS(:,:,:) ) / PRHODJ(:,:,:) &
 ! matrix inverted in TRIDIAG 
 !
 ZA(:,:,:)     = - PTSTEP * CSTURB%XCET * &
-                MZM(ZKEFF, D%NKA, KKU, KKL) * MZM(PRHODJ, D%NKA, KKU, KKL) / PDZZ**2
+                MZM(ZKEFF, D%NKA, D%NKU, D%NKL) * MZM(PRHODJ, D%NKA, D%NKU, D%NKL) / PDZZ**2
 !
 ! Compute TKE at time t+deltat: ( stored in ZRES )
 !
-CALL TRIDIAG_TKE(D,D%NKA,KKU,KKL,PTKEM,ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,&
+CALL TRIDIAG_TKE(D,D%NKA,D%NKU,D%NKL,PTKEM,ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,&
             & ZSOURCE,PTSTEP*ZFLX,ZRES)
 CALL GET_HALO(ZRES)
 !
@@ -313,21 +307,20 @@ IF ( OLES_CALL .OR.                         &
 !
 ! Compute the cartesian vertical flux of TKE in ZFLX
 !
-
-  ZFLX(:,:,:)   = - CSTURB%XCET * MZM(ZKEFF, D%NKA, KKU, KKL) *   &
-                  DZM(PIMPL * ZRES + PEXPL * PTKEM, D%NKA, KKU, KKL) / PDZZ
+  ZFLX(:,:,:)   = - CSTURB%XCET * MZM(ZKEFF, D%NKA, D%NKU, D%NKL) *   &
+                  DZM(PIMPL * ZRES + PEXPL * PTKEM, D%NKA, D%NKU, D%NKL) / PDZZ
 !
   ZFLX(:,:,IKB) = 0.
   ZFLX(:,:,D%NKA) = 0.
 !
 ! Compute the whole turbulent TRansport of TKE:
 !
-  ZTR(:,:,:)= ZTR - DZF(MZM(PRHODJ, D%NKA, KKU, KKL) * ZFLX / PDZZ, D%NKA, KKU, KKL) /PRHODJ
+  ZTR(:,:,:)= ZTR - DZF(MZM(PRHODJ, D%NKA, D%NKU, D%NKL) * ZFLX / PDZZ, D%NKA, D%NKU, D%NKL) /PRHODJ
 !
 ! Storage in the LES configuration
 !
   IF (OLES_CALL) THEN
-    CALL LES_MEAN_SUBGRID(MZF(ZFLX, D%NKA, KKU, KKL), X_LES_SUBGRID_WTke )
+    CALL LES_MEAN_SUBGRID(MZF(ZFLX, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WTke )
     CALL LES_MEAN_SUBGRID(-ZTR, X_LES_SUBGRID_ddz_WTke )
   END IF
 !
@@ -364,7 +357,7 @@ PRTKES(:,:,:) = ZRES(:,:,:) * PRHODJ(:,:,:) / PTSTEP -  PRTKEMS(:,:,:)
 ! stores the whole turbulent transport
 !
 IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TKE), 'TR', PRTKES(:, :, :) )
-
+!
 !----------------------------------------------------------------------------
 !
 !*       3.   COMPUTE THE DISSIPATIVE HEATING
@@ -448,9 +441,6 @@ IF (OLES_CALL ) THEN
 END IF
 !
 !----------------------------------------------------------------------------
-! 
-!
-!----------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('TKE_EPS_SOURCES',1,ZHOOK_HANDLE)
 END SUBROUTINE TKE_EPS_SOURCES
diff --git a/src/common/turb/mode_turb_ver.F90 b/src/common/turb/mode_turb_ver.F90
index 19b8f2219..d1eb3ae42 100644
--- a/src/common/turb/mode_turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -5,7 +5,7 @@
 MODULE MODE_TURB_VER
 IMPLICIT NONE
 CONTAINS
-SUBROUTINE TURB_VER(D,CST,CSTURB,TURBN,KKA,KKU,KKL,KRR,KRRL,KRRI,   &
+SUBROUTINE TURB_VER(D,CST,CSTURB,TURBN,KRR,KRRL,KRRI,   &
                       OTURB_FLX,OOCEAN,ODEEPOC,OHARAT,OCOMPUTE_SRC, &
                       KSV,KSV_LGBEG,KSV_LGEND,                      &
                       HTURBDIM,HTOM,PIMPL,PEXPL,                    &
@@ -247,9 +247,6 @@ TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
 TYPE(CST_t),            INTENT(IN)   :: CST
 TYPE(CSTURB_t),         INTENT(IN)   :: CSTURB
 TYPE(TURB_t),           INTENT(IN)   :: TURBN
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
 INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
 INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
 INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.
@@ -412,7 +409,7 @@ IKE=D%NKE
 ! 3D Redelsperger numbers
 !
 !
-CALL PRANDTL(D,CST,CSTURB,KKA,KKU,KKL,KRR,KSV,KRRI,OTURB_FLX, &
+CALL PRANDTL(D,CST,CSTURB,D%NKA,D%NKU,D%NKL,KRR,KSV,KRRI,OTURB_FLX, &
              HTURBDIM,OOCEAN,OHARAT,O2D,OCOMPUTE_SRC,&
              TPFILE,                               &
              PDXX,PDYY,PDZZ,PDZX,PDZY,             &
@@ -438,9 +435,9 @@ ZSQRT_TKE = SQRT(PTKEM)
 !
 ! gradients of mean quantities at previous time-step
 !
-ZDTH_DZ = GZ_M_W(KKA, KKU, KKL,PTHLM(:,:,:),PDZZ)
+ZDTH_DZ = GZ_M_W(D%NKA, D%NKU, D%NKL,PTHLM(:,:,:),PDZZ)
 ZDR_DZ  = 0.
-IF (KRR>0) ZDR_DZ  = GZ_M_W(KKA, KKU, KKL,PRM(:,:,:,1),PDZZ)
+IF (KRR>0) ZDR_DZ  = GZ_M_W(D%NKA, D%NKU, D%NKL,PRM(:,:,:,1),PDZZ)
 !
 !
 ! Denominator factor in 3rd order terms
@@ -496,7 +493,7 @@ ELSE
 ENDIF
 !
   CALL  TURB_VER_THERMO_FLUX(D,CST,CSTURB,TURBN,                      &
-                        KKA,KKU,KKL,KRR,KRRL,KRRI,KSV,                &
+                        KRR,KRRL,KRRI,KSV,                            &
                         OTURB_FLX,HTURBDIM,HTOM,OOCEAN,ODEEPOC,OHARAT,&
                         OCOUPLES,OLES_CALL,OCOMPUTE_SRC,              &
                         PIMPL,PEXPL,PTSTEP,HPROGRAM,TPFILE,           &
@@ -514,7 +511,7 @@ ENDIF
                         PRTHLS,PRRS,ZTHLP,ZRP,PTP,PWTH,PWRC )
 !
   CALL  TURB_VER_THERMO_CORR(D,CST,CSTURB,                            &
-                        KKA,KKU,KKL,KRR,KRRL,KRRI,KSV,                &
+                        KRR,KRRL,KRRI,KSV,                            &
                         OTURB_FLX,HTURBDIM,HTOM, OHARAT,OCOMPUTE_SRC, &
                         OCOUPLES,OLES_CALL,                           &                        
                         PIMPL,PEXPL,TPFILE,                           &
@@ -546,7 +543,7 @@ ENDIF
 !
 IF (OHARAT) ZLM=PLENGTHM
 !
-CALL  TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KKA,KKU,KKL,KSV,O2D,OFLAT,  &
+CALL  TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KSV,O2D,OFLAT,           &
                       OTURB_FLX,KRR,OOCEAN,OHARAT,OCOUPLES,OLES_CALL,&
                       HTURBDIM,PIMPL,PEXPL,PTSTEP,TPFILE,           &
                       PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
@@ -567,7 +564,7 @@ CALL  TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KKA,KKU,KKL,KSV,O2D,OFLAT,  &
 IF (OHARAT) ZLM=PLENGTHH
 !
 IF (SIZE(PSVM,4)>0)                                                 &
-CALL  TURB_VER_SV_FLUX(D,CST,CSTURB,KKA,KKU,KKL,ONOMIXLG,           &
+CALL  TURB_VER_SV_FLUX(D,CST,CSTURB,ONOMIXLG,                       &
                       KSV,KSV_LGBEG,KSV_LGEND,                      &
                       OTURB_FLX,HTURBDIM,OHARAT,OBLOWSNOW,OLES_CALL,&
                       PIMPL,PEXPL,PTSTEP,                           &
@@ -581,7 +578,7 @@ CALL  TURB_VER_SV_FLUX(D,CST,CSTURB,KKA,KKU,KKL,ONOMIXLG,           &
 !
 !
 IF (SIZE(PSVM,4)>0 .AND. OLES_CALL)                                 &
-CALL  TURB_VER_SV_CORR(D,CST,CSTURB,KKA,KKU,KKL,KRR,KRRL,KRRI,OOCEAN,&
+CALL  TURB_VER_SV_CORR(D,CST,CSTURB,KRR,KRRL,KRRI,OOCEAN,           &
                       PDZZ,KSV,KSV_LGBEG,KSV_LGEND,ONOMIXLG,        &
                       OBLOWSNOW,OLES_CALL,OCOMPUTE_SRC,             &
                       PTHLM,PRM,PTHVREF,                            &
diff --git a/src/common/turb/mode_turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
index f4aa457c5..8766729ad 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -5,7 +5,7 @@
 MODULE MODE_TURB_VER_DYN_FLUX
 IMPLICIT NONE
 CONTAINS
-SUBROUTINE TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KKA,KKU,KKL,KSV,O2D,OFLAT,&
+SUBROUTINE TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KSV,O2D,OFLAT,&
                       OTURB_FLX,KRR, OOCEAN,OHARAT,OCOUPLES,OLES_CALL,&
                       HTURBDIM,PIMPL,PEXPL,                         &
                       PTSTEP,                                       &
@@ -240,9 +240,6 @@ TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
 TYPE(CST_t),            INTENT(IN)   :: CST
 TYPE(CSTURB_t),         INTENT(IN)   :: CSTURB
 TYPE(TURB_t),           INTENT(IN)   :: TURBN
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
 INTEGER,                INTENT(IN)   :: KSV           ! number of scalar variables
 LOGICAL,                INTENT(IN)   ::  OTURB_FLX    ! switch to write the
                                  ! turbulent fluxes in the syncronous FM-file
@@ -362,7 +359,7 @@ ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
 IF (OHARAT) THEN
   ZKEFF(:,:,:) =  PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
 ELSE 
-  ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, KKU, KKL)
+  ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, D%NKU, D%NKL)
 ENDIF
 
 !
@@ -380,7 +377,7 @@ ZVSLOPEM(:,:,1)=PVSLOPEM(:,:)
 ! Preparation of the arguments for TRIDIAG_WIND 
 !
 ZA(:,:,:)    = -PTSTEP * ZCMFS *                              &
-              MXM( ZKEFF ) * MXM(MZM(PRHODJ, D%NKA, KKU, KKL)) / &
+              MXM( ZKEFF ) * MXM(MZM(PRHODJ, D%NKA, D%NKU, D%NKL)) / &
               MXM( PDZZ )**2
 !
 !
@@ -405,11 +402,11 @@ IF (OOCEAN) THEN  ! OCEAN MODEL ONLY
   ! Sfx flux assumed to be in SI & at vorticity point
   IF (OCOUPLES) THEN  
     ZSOURCE(:,:,IKE:IKE) = TURBN%XSSUFL_C(:,:,1:1)/PDZZ(:,:,IKE:IKE) &
-         *0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE))) 
+         *0.5 * ( 1. + MXM(PRHODJ(:,:,D%NKU:D%NKU)) / MXM(PRHODJ(:,:,IKE:IKE))) 
   ELSE
     ZSOURCE(:,:,IKE)     = XSSUFL(:,:)
     ZSOURCE(:,:,IKE:IKE) = ZSOURCE (:,:,IKE:IKE) /PDZZ(:,:,IKE:IKE) &
-        *0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
+        *0.5 * ( 1. + MXM(PRHODJ(:,:,D%NKU:D%NKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
   ENDIF
   !No flux at the ocean domain bottom
    ZSOURCE(:,:,IKB)           = 0.
@@ -444,7 +441,7 @@ ENDIF !end ocean or atmosphere cases
 !
 ! Obtention of the split U at t+ deltat 
 !
-CALL TRIDIAG_WIND(D%NKA,KKU,KKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL,   &
+CALL TRIDIAG_WIND(D%NKA,D%NKU,D%NKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL,   &
                   MXM(PRHODJ),ZSOURCE,ZRES)
 ! 
 !  Compute the equivalent tendency for the U wind component
@@ -457,7 +454,7 @@ PRUS(:,:,:)=PRUS(:,:,:)+MXM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PUM(:,:,:))/PTSTEP
 ! vertical flux of the U wind component
 !
 ZFLXZ(:,:,:)     = -ZCMFS * MXM(ZKEFF) * &
-                  DZM(PIMPL*ZRES + PEXPL*PUM, D%NKA, KKU, KKL) / MXM(PDZZ)
+                  DZM(PIMPL*ZRES + PEXPL*PUM, D%NKA, D%NKU, D%NKL) / MXM(PDZZ)
 !
 ! surface flux 
 ZFLXZ(:,:,IKB:IKB)   =   MXM(PDZZ(:,:,IKB:IKB))  *                &
@@ -470,8 +467,8 @@ ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB)
 IF (OOCEAN) THEN !ocean model at phys sfc (ocean domain top)
   ZFLXZ(:,:,IKE:IKE)   =   MXM(PDZZ(:,:,IKE:IKE))  *                &
                            ZSOURCE(:,:,IKE:IKE)                     &
-                           / 0.5 / ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
-  ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE) 
+                           / 0.5 / ( 1. + MXM(PRHODJ(:,:,D%NKU:D%NKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
+  ZFLXZ(:,:,D%NKU) = ZFLXZ(:,:,IKE) 
 END IF
 !
 IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -496,19 +493,19 @@ PWU(:,:,:) = ZFLXZ(:,:,:)
 ! Contribution to the dynamic production of TKE
 ! compute the dynamic production at the mass point
 !
-PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, D%NKA, KKU, KKL)), D%NKA, KKU, KKL)
+PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, D%NKA, D%NKU, D%NKL)), D%NKA, D%NKU, D%NKL)
 !
-! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
 PDP(:,:,IKB:IKB) = - MXF (                                                      &
-  ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PUM(:,:,IKB+KKL:IKB+KKL)-PUM(:,:,IKB:IKB))  &
-                         / MXM(PDZZ(:,:,IKB+KKL:IKB+KKL))                   &
+  ZFLXZ(:,:,IKB+D%NKL:IKB+D%NKL) * (PUM(:,:,IKB+D%NKL:IKB+D%NKL)-PUM(:,:,IKB:IKB))  &
+                         / MXM(PDZZ(:,:,IKB+D%NKL:IKB+D%NKL))                   &
                          )
 !
 IF (OOCEAN) THEN
-  ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+  ! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
   PDP(:,:,IKE:IKE) = - MXF (                                                      &
-    ZFLXZ(:,:,IKE-KKL:IKE-KKL) * (PUM(:,:,IKE:IKE)-PUM(:,:,IKE-KKL:IKE-KKL))  &
-                           / MXM(PDZZ(:,:,IKE-KKL:IKE-KKL))                   &
+    ZFLXZ(:,:,IKE-D%NKL:IKE-D%NKL) * (PUM(:,:,IKE:IKE)-PUM(:,:,IKE-D%NKL:IKE-D%NKL))  &
+                           / MXM(PDZZ(:,:,IKE-D%NKL:IKE-D%NKL))                   &
                            ) 
 END IF
 !
@@ -516,9 +513,9 @@ END IF
 ! 
 IF (OLES_CALL) THEN
   CALL SECOND_MNH(ZTIME1)
-  CALL LES_MEAN_SUBGRID(MZF(MXF(ZFLXZ), D%NKA, KKU, KKL), X_LES_SUBGRID_WU ) 
-  CALL LES_MEAN_SUBGRID(MZF(MXF(GZ_U_UW(PUM,PDZZ, D%NKA, KKU, KKL) &
-                       & *ZFLXZ), D%NKA, KKU, KKL), X_LES_RES_ddxa_U_SBG_UaU )
+  CALL LES_MEAN_SUBGRID(MZF(MXF(ZFLXZ), D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WU ) 
+  CALL LES_MEAN_SUBGRID(MZF(MXF(GZ_U_UW(PUM,PDZZ, D%NKA, D%NKU, D%NKL) &
+                       & *ZFLXZ), D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_U_SBG_UaU )
   CALL LES_MEAN_SUBGRID( ZCMFS * ZKEFF, X_LES_SUBGRID_Km )
   CALL SECOND_MNH(ZTIME2)
   XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
@@ -530,52 +527,52 @@ END IF
 IF(HTURBDIM=='3DIM') THEN
   ! Compute the source for the W wind component
                 ! used to compute the W source at the ground
-  ZFLXZ(:,:,D%NKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation 
+  ZFLXZ(:,:,D%NKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+D%NKL) ! extrapolation 
  IF (OOCEAN) THEN
-   ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation 
+   ZFLXZ(:,:,D%NKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-D%NKL) ! extrapolation 
  END IF     
      
   !
   IF (.NOT. OFLAT) THEN
     PRWS(:,:,:)= PRWS                                      &
-                -DXF( MZM(MXM(PRHODJ) /PDXX, D%NKA, KKU, KKL)  * ZFLXZ )  &
-                +DZM(PRHODJ / MZF(PDZZ, D%NKA, KKU, KKL) *                &
-                      MXF(MZF(ZFLXZ*PDZX, D%NKA, KKU, KKL) / PDXX ),      &
-                     D%NKA, KKU, KKL)
+                -DXF( MZM(MXM(PRHODJ) /PDXX, D%NKA, D%NKU, D%NKL)  * ZFLXZ )  &
+                +DZM(PRHODJ / MZF(PDZZ, D%NKA, D%NKU, D%NKL) *                &
+                      MXF(MZF(ZFLXZ*PDZX, D%NKA, D%NKU, D%NKL) / PDXX ),      &
+                     D%NKA, D%NKU, D%NKL)
   ELSE
-    PRWS(:,:,:)= PRWS -DXF(MZM(MXM(PRHODJ) /PDXX, D%NKA, KKU, KKL)  * ZFLXZ )
+    PRWS(:,:,:)= PRWS -DXF(MZM(MXM(PRHODJ) /PDXX, D%NKA, D%NKU, D%NKL)  * ZFLXZ )
   END IF
   !
   ! Complete the Dynamical production with the W wind component 
   !
-  ZA(:,:,:)=-MZF(MXF(ZFLXZ * GX_W_UW(PWM,PDXX,PDZZ,PDZX, D%NKA, KKU, KKL)), D%NKA, KKU, KKL)
+  ZA(:,:,:)=-MZF(MXF(ZFLXZ * GX_W_UW(PWM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)), D%NKA, D%NKU, D%NKL)
   !
   !
-  ! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+  ! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
   ZA(:,:,IKB:IKB) = - MXF (                                                  &
-   ZFLXZ(:,:,IKB+KKL:IKB+KKL) *                                              &
-     ( DXM( PWM(:,:,IKB+KKL:IKB+KKL) )                                       &
-      -MXM(  (PWM(:,:,IKB+2*KKL:IKB+2*KKL   )-PWM(:,:,IKB+KKL:IKB+KKL))      &
-              /(PDZZ(:,:,IKB+2*KKL:IKB+2*KKL)+PDZZ(:,:,IKB+KKL:IKB+KKL))     &
-            +(PWM(:,:,IKB+KKL:IKB+KKL)-PWM(:,:,IKB:IKB  ))                   &
-              /(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB  ))               &
+   ZFLXZ(:,:,IKB+D%NKL:IKB+D%NKL) *                                              &
+     ( DXM( PWM(:,:,IKB+D%NKL:IKB+D%NKL) )                                       &
+      -MXM(  (PWM(:,:,IKB+2*D%NKL:IKB+2*D%NKL   )-PWM(:,:,IKB+D%NKL:IKB+D%NKL))      &
+              /(PDZZ(:,:,IKB+2*D%NKL:IKB+2*D%NKL)+PDZZ(:,:,IKB+D%NKL:IKB+D%NKL))     &
+            +(PWM(:,:,IKB+D%NKL:IKB+D%NKL)-PWM(:,:,IKB:IKB  ))                   &
+              /(PDZZ(:,:,IKB+D%NKL:IKB+D%NKL)+PDZZ(:,:,IKB:IKB  ))               &
           )                                                                  &
-        * PDZX(:,:,IKB+KKL:IKB+KKL)                                          &
-     ) / (0.5*(PDXX(:,:,IKB+KKL:IKB+KKL)+PDXX(:,:,IKB:IKB)))                 &
+        * PDZX(:,:,IKB+D%NKL:IKB+D%NKL)                                          &
+     ) / (0.5*(PDXX(:,:,IKB+D%NKL:IKB+D%NKL)+PDXX(:,:,IKB:IKB)))                 &
                           )
   !
 IF (OOCEAN) THEN
-  ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+  ! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
   ZA(:,:,IKE:IKE) = - MXF (                                                  &
-   ZFLXZ(:,:,IKE-KKL:IKE-KKL) *                                              &
-     ( DXM( PWM(:,:,IKE-KKL:IKE-KKL) )                                       &
-      -MXM(  (PWM(:,:,IKE-2*KKL:IKE-2*KKL   )-PWM(:,:,IKE-KKL:IKE-KKL))      &
-              /(PDZZ(:,:,IKE-2*KKL:IKE-2*KKL)+PDZZ(:,:,IKE-KKL:IKE-KKL))     &
-            +(PWM(:,:,IKE-KKL:IKE-KKL)-PWM(:,:,IKE:IKE  ))                   &
-              /(PDZZ(:,:,IKE-KKL:IKE-KKL)+PDZZ(:,:,IKE:IKE  ))               &
+   ZFLXZ(:,:,IKE-D%NKL:IKE-D%NKL) *                                              &
+     ( DXM( PWM(:,:,IKE-D%NKL:IKE-D%NKL) )                                       &
+      -MXM(  (PWM(:,:,IKE-2*D%NKL:IKE-2*D%NKL   )-PWM(:,:,IKE-D%NKL:IKE-D%NKL))      &
+              /(PDZZ(:,:,IKE-2*D%NKL:IKE-2*D%NKL)+PDZZ(:,:,IKE-D%NKL:IKE-D%NKL))     &
+            +(PWM(:,:,IKE-D%NKL:IKE-D%NKL)-PWM(:,:,IKE:IKE  ))                   &
+              /(PDZZ(:,:,IKE-D%NKL:IKE-D%NKL)+PDZZ(:,:,IKE:IKE  ))               &
           )                                                                  &
-         * PDZX(:,:,IKE-KKL:IKE-KKL)                                         &
-     ) / (0.5*(PDXX(:,:,IKE-KKL:IKE-KKL)+PDXX(:,:,IKE:IKE)))                 &
+         * PDZX(:,:,IKE-D%NKL:IKE-D%NKL)                                         &
+     ) / (0.5*(PDXX(:,:,IKE-D%NKL:IKE-D%NKL)+PDXX(:,:,IKE:IKE)))                 &
                           )
 END IF
   !
@@ -586,16 +583,16 @@ END IF
   IF (OLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
     CALL LES_MEAN_SUBGRID(MZF(MXF(GX_W_UW(PWM,PDXX,&
-      PDZZ,PDZX, D%NKA, KKU, KKL)*ZFLXZ), D%NKA, KKU, KKL), X_LES_RES_ddxa_W_SBG_UaW )
-    CALL LES_MEAN_SUBGRID(MXF(GX_M_U(D%NKA, KKU, KKL,PTHLM,PDXX,PDZZ,PDZX)&
-      * MZF(ZFLXZ, D%NKA, KKU, KKL)), X_LES_RES_ddxa_Thl_SBG_UaW )
+      PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*ZFLXZ), D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_W_SBG_UaW )
+    CALL LES_MEAN_SUBGRID(MXF(GX_M_U(D%NKA, D%NKU, D%NKL,PTHLM,PDXX,PDZZ,PDZX)&
+      * MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)), X_LES_RES_ddxa_Thl_SBG_UaW )
     IF (KRR>=1) THEN
-      CALL LES_MEAN_SUBGRID(MXF(GX_U_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, KKU, KKL)&
-      *MZF(ZFLXZ, D%NKA, KKU, KKL)),X_LES_RES_ddxa_Rt_SBG_UaW )
+      CALL LES_MEAN_SUBGRID(MXF(GX_U_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)&
+      *MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)),X_LES_RES_ddxa_Rt_SBG_UaW )
     END IF
     DO JSV=1,KSV
       CALL LES_MEAN_SUBGRID( MXF(GX_U_M(PSVM(:,:,:,JSV),PDXX,PDZZ,&
-      PDZX, D%NKA, KKU, KKL)*MZF(ZFLXZ, D%NKA, KKU, KKL)),X_LES_RES_ddxa_Sv_SBG_UaW(:,:,:,JSV) )
+      PDZX, D%NKA, D%NKU, D%NKL)*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)),X_LES_RES_ddxa_Sv_SBG_UaW(:,:,:,JSV) )
     END DO
     CALL SECOND_MNH(ZTIME2)
     XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
@@ -613,7 +610,7 @@ END IF
 ! Preparation of the arguments for TRIDIAG_WIND
 !!
 ZA(:,:,:)    = - PTSTEP * ZCMFS *                              &
-              MYM( ZKEFF ) * MYM(MZM(PRHODJ, D%NKA, KKU, KKL)) / &
+              MYM( ZKEFF ) * MYM(MZM(PRHODJ, D%NKA, D%NKU, D%NKL)) / &
               MYM( PDZZ )**2
 !
 !
@@ -635,11 +632,11 @@ ZCOEFS(:,:,1:1)=MYM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
 IF (OOCEAN) THEN ! Ocean case
   IF (OCOUPLES) THEN
     ZSOURCE(:,:,IKE:IKE) =  TURBN%XSSVFL_C(:,:,1:1)/PDZZ(:,:,IKE:IKE) &
-        *0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) ) 
+        *0.5 * ( 1. + MYM(PRHODJ(:,:,D%NKU:D%NKU)) / MYM(PRHODJ(:,:,IKE:IKE)) ) 
   ELSE 
     ZSOURCE(:,:,IKE) = XSSVFL(:,:)
     ZSOURCE(:,:,IKE:IKE) = ZSOURCE(:,:,IKE:IKE)/PDZZ(:,:,IKE:IKE) &
-        *0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+        *0.5 * ( 1. + MYM(PRHODJ(:,:,D%NKU:D%NKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
   END IF
   !No flux at the ocean domain bottom
   ZSOURCE(:,:,IKB) = 0.
@@ -672,7 +669,7 @@ ENDIF ! End of Ocean or Atmospher Cases
 ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
 ! 
 !  Obtention of the split V at t+ deltat 
-CALL TRIDIAG_WIND(D%NKA,KKU,KKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL,  &
+CALL TRIDIAG_WIND(D%NKA,D%NKU,D%NKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL,  &
                   MYM(PRHODJ),ZSOURCE,ZRES)
 !
 ! Compute the equivalent tendency for the V wind component
@@ -685,7 +682,7 @@ PRVS(:,:,:)=PRVS(:,:,:)+MYM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PVM(:,:,:))/PTSTEP
 !  vertical flux of the V wind component
 !
 ZFLXZ(:,:,:)   = -ZCMFS * MYM(ZKEFF) * &
-                DZM(PIMPL*ZRES + PEXPL*PVM, D%NKA, KKU, KKL) / MYM(PDZZ)
+                DZM(PIMPL*ZRES + PEXPL*PVM, D%NKA, D%NKU, D%NKL) / MYM(PDZZ)
 !
 ZFLXZ(:,:,IKB:IKB)   =   MYM(PDZZ(:,:,IKB:IKB))  *                       &
   ( ZSOURCE(:,:,IKB:IKB)                                                 &
@@ -698,8 +695,8 @@ ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB)
 IF (OOCEAN) THEN
   ZFLXZ(:,:,IKE:IKE)   =   MYM(PDZZ(:,:,IKE:IKE))  *                &
       ZSOURCE(:,:,IKE:IKE)                                          &
-      / 0.5 / ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
-  ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE) 
+      / 0.5 / ( 1. + MYM(PRHODJ(:,:,D%NKU:D%NKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+  ZFLXZ(:,:,D%NKU) = ZFLXZ(:,:,IKE) 
 END IF
 !
 IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -724,20 +721,20 @@ PWV(:,:,:) = ZFLXZ(:,:,:)
 !  Contribution to the dynamic production of TKE
 ! compute the dynamic production contribution at the mass point
 !
-ZA(:,:,:) = - MZF(MYF(ZFLXZ * GZ_V_VW(PVM,PDZZ, D%NKA, KKU, KKL)), D%NKA, KKU, KKL)
+ZA(:,:,:) = - MZF(MYF(ZFLXZ * GZ_V_VW(PVM,PDZZ, D%NKA, D%NKU, D%NKL)), D%NKA, D%NKU, D%NKL)
 !
-! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
 ZA(:,:,IKB:IKB)  =                                                 &
                  - MYF (                                          &
-ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PVM(:,:,IKB+KKL:IKB+KKL)-PVM(:,:,IKB:IKB))  &
-                       / MYM(PDZZ(:,:,IKB+KKL:IKB+KKL))               &
+ZFLXZ(:,:,IKB+D%NKL:IKB+D%NKL) * (PVM(:,:,IKB+D%NKL:IKB+D%NKL)-PVM(:,:,IKB:IKB))  &
+                       / MYM(PDZZ(:,:,IKB+D%NKL:IKB+D%NKL))               &
                        )
 !
 IF (OOCEAN) THEN
-  ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+  ! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
   ZA(:,:,IKE:IKE) = - MYF (                                                  &
-   ZFLXZ(:,:,IKE-KKL:IKE-KKL) * (PVM(:,:,IKE:IKE)-PVM(:,:,IKE-KKL:IKE-KKL))  &
-                          / MYM(PDZZ(:,:,IKE-KKL:IKE-KKL))                   &
+   ZFLXZ(:,:,IKE-D%NKL:IKE-D%NKL) * (PVM(:,:,IKE:IKE)-PVM(:,:,IKE-D%NKL:IKE-D%NKL))  &
+                          / MYM(PDZZ(:,:,IKE-D%NKL:IKE-D%NKL))                   &
                           )
 END IF
 !
@@ -747,9 +744,9 @@ PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
 !
 IF (OLES_CALL) THEN
   CALL SECOND_MNH(ZTIME1)
-  CALL LES_MEAN_SUBGRID(MZF(MYF(ZFLXZ), D%NKA, KKU, KKL), X_LES_SUBGRID_WV ) 
-  CALL LES_MEAN_SUBGRID(MZF(MYF(GZ_V_VW(PVM,PDZZ, D%NKA, KKU, KKL)*&
-                    & ZFLXZ), D%NKA, KKU, KKL), X_LES_RES_ddxa_V_SBG_UaV )
+  CALL LES_MEAN_SUBGRID(MZF(MYF(ZFLXZ), D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WV ) 
+  CALL LES_MEAN_SUBGRID(MZF(MYF(GZ_V_VW(PVM,PDZZ, D%NKA, D%NKU, D%NKL)*&
+                    & ZFLXZ), D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_V_SBG_UaV )
   CALL SECOND_MNH(ZTIME2)
   XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
 END IF
@@ -759,51 +756,51 @@ END IF
 !
 IF(HTURBDIM=='3DIM') THEN
   ! Compute the source for the W wind component
-  ZFLXZ(:,:,D%NKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation 
+  ZFLXZ(:,:,D%NKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+D%NKL) ! extrapolation 
   IF (OOCEAN) THEN
-    ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation 
+    ZFLXZ(:,:,D%NKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-D%NKL) ! extrapolation 
   END IF
   !
   IF (.NOT. O2D) THEN 
     IF (.NOT. OFLAT) THEN
       PRWS(:,:,:)= PRWS(:,:,:)                               &
-                  -DYF( MZM(MYM(PRHODJ) /PDYY, D%NKA, KKU, KKL) * ZFLXZ )   &
-                  +DZM(PRHODJ / MZF(PDZZ, D%NKA, KKU, KKL) *                &
-                        MYF(MZF(ZFLXZ*PDZY, D%NKA, KKU, KKL) / PDYY ),      &
-                       D%NKA, KKU, KKL)
+                  -DYF( MZM(MYM(PRHODJ) /PDYY, D%NKA, D%NKU, D%NKL) * ZFLXZ )   &
+                  +DZM(PRHODJ / MZF(PDZZ, D%NKA, D%NKU, D%NKL) *                &
+                        MYF(MZF(ZFLXZ*PDZY, D%NKA, D%NKU, D%NKL) / PDYY ),      &
+                       D%NKA, D%NKU, D%NKL)
     ELSE
-      PRWS(:,:,:)= PRWS(:,:,:) -DYF(MZM(MYM(PRHODJ) /PDYY, D%NKA, KKU, KKL) * ZFLXZ )
+      PRWS(:,:,:)= PRWS(:,:,:) -DYF(MZM(MYM(PRHODJ) /PDYY, D%NKA, D%NKU, D%NKL) * ZFLXZ )
     END IF
   END IF
   ! 
   ! Complete the Dynamical production with the W wind component 
   IF (.NOT. O2D) THEN
-    ZA(:,:,:) = - MZF(MYF(ZFLXZ * GY_W_VW(PWM,PDYY,PDZZ,PDZY, D%NKA, KKU, KKL)), D%NKA, KKU, KKL)
+    ZA(:,:,:) = - MZF(MYF(ZFLXZ * GY_W_VW(PWM,PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL)), D%NKA, D%NKU, D%NKL)
   !
-  ! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+  ! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
     ZA(:,:,IKB:IKB) = - MYF (                                              &
-     ZFLXZ(:,:,IKB+KKL:IKB+KKL) *                                          &
-       ( DYM( PWM(:,:,IKB+KKL:IKB+KKL) )                                   &
-        -MYM(  (PWM(:,:,IKB+2*KKL:IKB+2*KKL)-PWM(:,:,IKB+KKL:IKB+KKL))     &
-                /(PDZZ(:,:,IKB+2*KKL:IKB+2*KKL)+PDZZ(:,:,IKB+KKL:IKB+KKL)) &
-              +(PWM(:,:,IKB+KKL:IKB+KKL)-PWM(:,:,IKB:IKB  ))               &
-                /(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB  ))           &
+     ZFLXZ(:,:,IKB+D%NKL:IKB+D%NKL) *                                          &
+       ( DYM( PWM(:,:,IKB+D%NKL:IKB+D%NKL) )                                   &
+        -MYM(  (PWM(:,:,IKB+2*D%NKL:IKB+2*D%NKL)-PWM(:,:,IKB+D%NKL:IKB+D%NKL))     &
+                /(PDZZ(:,:,IKB+2*D%NKL:IKB+2*D%NKL)+PDZZ(:,:,IKB+D%NKL:IKB+D%NKL)) &
+              +(PWM(:,:,IKB+D%NKL:IKB+D%NKL)-PWM(:,:,IKB:IKB  ))               &
+                /(PDZZ(:,:,IKB+D%NKL:IKB+D%NKL)+PDZZ(:,:,IKB:IKB  ))           &
             )                                                              &
-          * PDZY(:,:,IKB+KKL:IKB+KKL)                                      &
-       ) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB)))             &
+          * PDZY(:,:,IKB+D%NKL:IKB+D%NKL)                                      &
+       ) / (0.5*(PDYY(:,:,IKB+D%NKL:IKB+D%NKL)+PDYY(:,:,IKB:IKB)))             &
                             )
   !
     IF (OOCEAN) THEN
      ZA(:,:,IKE:IKE) = - MYF (                                              &
-      ZFLXZ(:,:,IKE-KKL:IKE-KKL) *                                          &
-        ( DYM( PWM(:,:,IKE-KKL:IKE-KKL) )                                   &
-         -MYM(  (PWM(:,:,IKE-2*KKL:IKE-2*KKL)-PWM(:,:,IKE-KKL:IKE-KKL))     &
-                 /(PDZZ(:,:,IKE-2*KKL:IKE-2*KKL)+PDZZ(:,:,IKE-KKL:IKE-KKL)) &
-               +(PWM(:,:,IKE-KKL:IKE-KKL)-PWM(:,:,IKE:IKE  ))               &
-                 /(PDZZ(:,:,IKE-KKL:IKE-KKL)+PDZZ(:,:,IKE:IKE  ))           &
+      ZFLXZ(:,:,IKE-D%NKL:IKE-D%NKL) *                                          &
+        ( DYM( PWM(:,:,IKE-D%NKL:IKE-D%NKL) )                                   &
+         -MYM(  (PWM(:,:,IKE-2*D%NKL:IKE-2*D%NKL)-PWM(:,:,IKE-D%NKL:IKE-D%NKL))     &
+                 /(PDZZ(:,:,IKE-2*D%NKL:IKE-2*D%NKL)+PDZZ(:,:,IKE-D%NKL:IKE-D%NKL)) &
+               +(PWM(:,:,IKE-D%NKL:IKE-D%NKL)-PWM(:,:,IKE:IKE  ))               &
+                 /(PDZZ(:,:,IKE-D%NKL:IKE-D%NKL)+PDZZ(:,:,IKE:IKE  ))           &
              )                                                              &
-           * PDZY(:,:,IKE-KKL:IKE-KKL)                                      &
-        ) / (0.5*(PDYY(:,:,IKE-KKL:IKE-KKL)+PDYY(:,:,IKE:IKE)))             &
+           * PDZY(:,:,IKE-D%NKL:IKE-D%NKL)                                      &
+        ) / (0.5*(PDYY(:,:,IKE-D%NKL:IKE-D%NKL)+PDYY(:,:,IKE:IKE)))             &
                             )
     END IF
 !    
@@ -816,14 +813,14 @@ IF(HTURBDIM=='3DIM') THEN
   IF (OLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
     CALL LES_MEAN_SUBGRID(MZF(MYF(GY_W_VW(PWM,PDYY,&
-                         &PDZZ,PDZY, D%NKA, KKU, KKL)*ZFLXZ), D%NKA, KKU, KKL), &
+                         &PDZZ,PDZY, D%NKA, D%NKU, D%NKL)*ZFLXZ), D%NKA, D%NKU, D%NKL), &
                          &X_LES_RES_ddxa_W_SBG_UaW , .TRUE. )
-    CALL LES_MEAN_SUBGRID(MYF(GY_M_V(D%NKA, KKU, KKL,PTHLM,PDYY,PDZZ,PDZY)*&
-                         &MZF(ZFLXZ, D%NKA, KKU, KKL)), &
+    CALL LES_MEAN_SUBGRID(MYF(GY_M_V(D%NKA, D%NKU, D%NKL,PTHLM,PDYY,PDZZ,PDZY)*&
+                         &MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)), &
                          &X_LES_RES_ddxa_Thl_SBG_UaW , .TRUE. )
     IF (KRR>=1) THEN
       CALL LES_MEAN_SUBGRID(MYF(GY_V_M(PRM(:,:,:,1),PDYY,PDZZ,&
-                           &PDZY, D%NKA, KKU, KKL)*MZF(ZFLXZ, D%NKA, KKU, KKL)),&
+                           &PDZY, D%NKA, D%NKU, D%NKL)*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)),&
                            &X_LES_RES_ddxa_Rt_SBG_UaW , .TRUE. )
     END IF
     CALL SECOND_MNH(ZTIME2)
@@ -840,7 +837,7 @@ END IF
 !
 IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. HTURBDIM == '1DIM') THEN
   ZFLXZ(:,:,:)= (2./3.) * PTKEM(:,:,:)                     &
-     -ZCMFS*PLM(:,:,:)*SQRT(PTKEM(:,:,:))*GZ_W_M(PWM,PDZZ, D%NKA, KKU, KKL)
+     -ZCMFS*PLM(:,:,:)*SQRT(PTKEM(:,:,:))*GZ_W_M(PWM,PDZZ, D%NKA, D%NKU, D%NKL)
   ! to be tested &
   !   +XCMFB*(4./3.)*PLM(:,:,:)/SQRT(PTKEM(:,:,:))*PTP(:,:,:) 
   ! stores the W variance
diff --git a/src/common/turb/mode_turb_ver_sv_corr.F90 b/src/common/turb/mode_turb_ver_sv_corr.F90
index d77d74449..55a477d10 100644
--- a/src/common/turb/mode_turb_ver_sv_corr.F90
+++ b/src/common/turb/mode_turb_ver_sv_corr.F90
@@ -5,7 +5,7 @@
 MODULE MODE_TURB_VER_SV_CORR
 IMPLICIT NONE
 CONTAINS
-SUBROUTINE TURB_VER_SV_CORR(D,CST,CSTURB,KKA,KKU,KKL,KRR,KRRL,KRRI,OOCEAN,&
+SUBROUTINE TURB_VER_SV_CORR(D,CST,CSTURB,KRR,KRRL,KRRI,OOCEAN,&
                       PDZZ,KSV,KSV_LGBEG,KSV_LGEND,ONOMIXLG,        &
                       OBLOWSNOW,OLES_CALL,OCOMPUTE_SRC,             &
                       PTHLM,PRM,PTHVREF,                            &
@@ -81,12 +81,10 @@ IMPLICIT NONE
 !
 !
 !
-TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
-TYPE(CST_t),            INTENT(IN)    :: CST
-TYPE(CSTURB_t),         INTENT(IN)    :: CSTURB
-INTEGER,                INTENT(IN)  :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,                INTENT(IN)  :: KKL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-INTEGER,                INTENT(IN)   :: KSV, KSV_LGBEG, KSV_LGEND ! number of scalar variables
+TYPE(DIMPHYEX_t),       INTENT(IN)   ::  D
+TYPE(CST_t),            INTENT(IN)   ::  CST
+TYPE(CSTURB_t),         INTENT(IN)   ::  CSTURB
+INTEGER,                INTENT(IN)   ::  KSV, KSV_LGBEG, KSV_LGEND ! number of scalar variables
 LOGICAL,                INTENT(IN)   ::  OOCEAN       ! switch for Ocean model version
 LOGICAL,                INTENT(IN)   ::  ONOMIXLG     ! to use turbulence for lagrangian variables (modd_conf)
 LOGICAL,                INTENT(IN)   ::  OLES_CALL    ! compute the LES diagnostics at current time-step
@@ -156,10 +154,10 @@ DO JSV=1,KSV
   !
   IF (OLES_CALL) THEN
     ! approximation: diagnosed explicitely (without implicit term)
-    ZFLXZ(:,:,:) =  PPSI_SV(:,:,:,JSV)*GZ_M_W(D%NKA, KKU, KKL,PSVM(:,:,:,JSV),PDZZ)**2
-    ZFLXZ(:,:,:) = ZCSV / ZCSVD * PLM * PLEPS * MZF(ZFLXZ(:,:,:), D%NKA, KKU, KKL)
+    ZFLXZ(:,:,:) =  PPSI_SV(:,:,:,JSV)*GZ_M_W(D%NKA, D%NKU, D%NKL,PSVM(:,:,:,JSV),PDZZ)**2
+    ZFLXZ(:,:,:) = ZCSV / ZCSVD * PLM * PLEPS * MZF(ZFLXZ(:,:,:), D%NKA, D%NKU, D%NKL)
     CALL LES_MEAN_SUBGRID(-2.*ZCSVD*SQRT(PTKEM)*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV) )
-    CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Sv2(:,:,:,JSV) )
+    CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, D%NKU, D%NKL)*ZFLXZ, X_LES_RES_W_SBG_Sv2(:,:,:,JSV) )
   END IF
   !
   ! covariance ThvSv
@@ -168,18 +166,18 @@ DO JSV=1,KSV
     ! approximation: diagnosed explicitely (without implicit term)
     ZA(:,:,:)   =  ETHETA(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM,OOCEAN,OCOMPUTE_SRC)
     ZFLXZ(:,:,:)= ( CSTURB%XCSHF * PPHI3 + ZCSV * PPSI_SV(:,:,:,JSV) )              &
-                  *  GZ_M_W(D%NKA, KKU, KKL,PTHLM,PDZZ)                          &
-                  *  GZ_M_W(D%NKA, KKU, KKL,PSVM(:,:,:,JSV),PDZZ)
-    ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCTSVD) * MZF(ZFLXZ, D%NKA, KKU, KKL)
+                  *  GZ_M_W(D%NKA, D%NKU, D%NKL,PTHLM,PDZZ)                          &
+                  *  GZ_M_W(D%NKA, D%NKU, D%NKL,PSVM(:,:,:,JSV),PDZZ)
+    ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCTSVD) * MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)
     CALL LES_MEAN_SUBGRID( ZA*ZFLXZ, X_LES_SUBGRID_SvThv(:,:,:,JSV) ) 
     CALL LES_MEAN_SUBGRID( -CST%XG/PTHVREF/3.*ZA*ZFLXZ, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE.)
     !
     IF (KRR>=1) THEN
       ZA(:,:,:)   =  EMOIST(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM,OOCEAN)
       ZFLXZ(:,:,:)= ( ZCSV * PPSI3 + ZCSV * PPSI_SV(:,:,:,JSV) )             &
-                    *  GZ_M_W(D%NKA, KKU, KKL,PRM(:,:,:,1),PDZZ)                 &
-                    *  GZ_M_W(D%NKA, KKU, KKL,PSVM(:,:,:,JSV),PDZZ)
-      ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCQSVD) * MZF(ZFLXZ, D%NKA, KKU, KKL)
+                    *  GZ_M_W(D%NKA, D%NKU, D%NKL,PRM(:,:,:,1),PDZZ)                 &
+                    *  GZ_M_W(D%NKA, D%NKU, D%NKL,PSVM(:,:,:,JSV),PDZZ)
+      ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCQSVD) * MZF(ZFLXZ, D%NKA, D%NKU, D%NKL)
       CALL LES_MEAN_SUBGRID( ZA*ZFLXZ, X_LES_SUBGRID_SvThv(:,:,:,JSV) , .TRUE.)
       CALL LES_MEAN_SUBGRID( -CST%XG/PTHVREF/3.*ZA*ZFLXZ, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE.)
     END IF
diff --git a/src/common/turb/mode_turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
index 52487054d..892a5c360 100644
--- a/src/common/turb/mode_turb_ver_sv_flux.F90
+++ b/src/common/turb/mode_turb_ver_sv_flux.F90
@@ -5,7 +5,7 @@
 MODULE MODE_TURB_VER_SV_FLUX
 IMPLICIT NONE
 CONTAINS
-SUBROUTINE TURB_VER_SV_FLUX(D,CST,CSTURB,KKA,KKU,KKL,ONOMIXLG,      &
+SUBROUTINE TURB_VER_SV_FLUX(D,CST,CSTURB,ONOMIXLG,                  &
                       KSV,KSV_LGBEG,KSV_LGEND,                      &
                       OTURB_FLX,HTURBDIM,OHARAT,OBLOWSNOW,OLES_CALL,&
                       PIMPL,PEXPL,                                  &
@@ -241,9 +241,6 @@ IMPLICIT NONE
 TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
 TYPE(CST_t),            INTENT(IN)   :: CST
 TYPE(CSTURB_t),         INTENT(IN)   :: CSTURB
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
 INTEGER,                INTENT(IN)   :: KSV, &
                                         KSV_LGBEG, KSV_LGEND ! number of scalar variables
 LOGICAL,                INTENT(IN)   ::  OTURB_FLX    ! switch to write the
@@ -282,9 +279,6 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KSV), INTENT(INOUT) ::  PRSVS
                             ! cumulated sources for the prognostic variables
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KSV), INTENT(OUT)  :: PWSV        ! scalar flux
 !
-!
-!
-!
 !*       0.2  declaration of local variables
 !
 !
@@ -304,7 +298,6 @@ INTEGER             :: IKT          ! array size in k direction
 INTEGER             :: IKTB,IKTE    ! start, end of k loops in physical domain 
 INTEGER             :: JSV          ! loop counters
 INTEGER             :: JK           ! loop
-INTEGER             :: ISV          ! number of scalar var.
 !
 REAL :: ZTIME1, ZTIME2
 
@@ -326,12 +319,10 @@ IKTE=D%NKTE
 IKB=D%NKB
 IKE=D%NKE             
 !
-ISV=SIZE(PSVM,4)
-!
 IF (OHARAT) THEN
   ZKEFF(:,:,:) =  PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50.*MFMOIST(:,:,:)
 ELSE
-  ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, KKU, KKL)
+  ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, D%NKU, D%NKL)
 ENDIF
 !
 IF(OBLOWSNOW) THEN
@@ -345,18 +336,18 @@ ENDIF
 !*       8.   SOURCES OF PASSIVE SCALAR VARIABLES
 !             -----------------------------------
 !
-DO JSV=1,ISV
+DO JSV=1,KSV
 !
   IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
 !
 ! Preparation of the arguments for TRIDIAG 
     IF (OHARAT) THEN
       ZA(:,:,:)    = -PTSTEP*   &
-                   ZKEFF * MZM(PRHODJ, D%NKA, KKU, KKL) /   &
+                   ZKEFF * MZM(PRHODJ, D%NKA, D%NKU, D%NKL) /   &
                    PDZZ**2
     ELSE
       ZA(:,:,:)    = -PTSTEP*ZCSV*PPSI_SV(:,:,:,JSV) *   &
-                   ZKEFF * MZM(PRHODJ, D%NKA, KKU, KKL) /   &
+                   ZKEFF * MZM(PRHODJ, D%NKA, D%NKU, D%NKL) /   &
                    PDZZ**2
     ENDIF
   ZSOURCE(:,:,:) = 0.
@@ -381,7 +372,7 @@ DO JSV=1,ISV
   ZSOURCE(:,:,IKE) = 0.
 !
 ! Obtention of the split JSV scalar variable at t+ deltat  
-  CALL TRIDIAG(D,D%NKA,KKU,KKL,PSVM(:,:,:,JSV),ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,ZSOURCE,ZRES)
+  CALL TRIDIAG(D,D%NKA,D%NKU,D%NKL,PSVM(:,:,:,JSV),ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,ZSOURCE,ZRES)
 !
 !  Compute the equivalent tendency for the JSV scalar variable
   PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV)+    &
@@ -390,8 +381,8 @@ DO JSV=1,ISV
   IF ( (OTURB_FLX .AND. TPFILE%LOPENED) .OR. OLES_CALL ) THEN
     ! Diagnostic of the cartesian vertical flux
     !
-    ZFLXZ(:,:,:) = -ZCSV * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM), D%NKA, KKU, KKL) / PDZZ * &
-                  DZM(PIMPL*ZRES(:,:,:) + PEXPL*PSVM(:,:,:,JSV), D%NKA, KKU, KKL)
+    ZFLXZ(:,:,:) = -ZCSV * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM), D%NKA, D%NKU, D%NKL) / PDZZ * &
+                  DZM(PIMPL*ZRES(:,:,:) + PEXPL*PSVM(:,:,:,JSV), D%NKA, D%NKU, D%NKL)
     ! surface flux
     !* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
     ! (in presence of slopes)
@@ -408,10 +399,10 @@ DO JSV=1,ISV
     ! the IKB flux gives the ground value
     ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB)
     DO JK=IKTB+1,IKTE-1
-      PWSV(:,:,JK,JSV)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+      PWSV(:,:,JK,JSV)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+D%NKL))
     END DO
-    PWSV(:,:,IKB,JSV)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
-    PWSV(:,:,IKE,JSV)=PWSV(:,:,IKE-KKL,JSV)
+    PWSV(:,:,IKB,JSV)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+D%NKL))
+    PWSV(:,:,IKE,JSV)=PWSV(:,:,IKE-D%NKL,JSV)
  END IF
   !
   IF (OTURB_FLX .AND. TPFILE%LOPENED) THEN
@@ -435,13 +426,13 @@ DO JSV=1,ISV
   !
   IF (OLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
-    CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_WSv(:,:,:,JSV) )
-    CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, D%NKA, KKU, KKL)*MZF(ZFLXZ, D%NKA, KKU, KKL), &
+    CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WSv(:,:,:,JSV) )
+    CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, D%NKA, D%NKU, D%NKL)*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), &
                           X_LES_RES_ddxa_W_SBG_UaSv(:,:,:,JSV) )
-    CALL LES_MEAN_SUBGRID(MZF(GZ_M_W(D%NKA, KKU, KKL,PSVM(:,:,:,JSV),PDZZ)*ZFLXZ, D%NKA, KKU, KKL), &
+    CALL LES_MEAN_SUBGRID(MZF(GZ_M_W(D%NKA, D%NKU, D%NKL,PSVM(:,:,:,JSV),PDZZ)*ZFLXZ, D%NKA, D%NKU, D%NKL), &
                           X_LES_RES_ddxa_Sv_SBG_UaSv(:,:,:,JSV) )
-    CALL LES_MEAN_SUBGRID(-ZCSVP*SQRT(PTKEM)/PLM*MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_SvPz(:,:,:,JSV) )
-    CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_W_SBG_WSv(:,:,:,JSV) )
+    CALL LES_MEAN_SUBGRID(-ZCSVP*SQRT(PTKEM)/PLM*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_SvPz(:,:,:,JSV) )
+    CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_W_SBG_WSv(:,:,:,JSV) )
     CALL SECOND_MNH(ZTIME2)
     XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
   END IF
diff --git a/src/common/turb/mode_turb_ver_thermo_corr.F90 b/src/common/turb/mode_turb_ver_thermo_corr.F90
index 66b704cb5..80849004c 100644
--- a/src/common/turb/mode_turb_ver_thermo_corr.F90
+++ b/src/common/turb/mode_turb_ver_thermo_corr.F90
@@ -6,7 +6,7 @@ MODULE MODE_TURB_VER_THERMO_CORR
 IMPLICIT NONE
 CONTAINS      
 SUBROUTINE TURB_VER_THERMO_CORR(D,CST,CSTURB,                       &
-                      KKA,KKU,KKL,KRR,KRRL,KRRI,KSV,                &
+                      KRR,KRRL,KRRI,KSV,                            &
                       OTURB_FLX,HTURBDIM,HTOM,OHARAT,OCOMPUTE_SRC,  &
                       OCOUPLES,OLES_CALL,                           &
                       PIMPL,PEXPL,TPFILE,                           &
@@ -235,9 +235,6 @@ IMPLICIT NONE
 TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
 TYPE(CST_t),            INTENT(IN)   :: CST
 TYPE(CSTURB_t),         INTENT(IN)   :: CSTURB
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
 INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
 INTEGER,                INTENT(IN)   :: KSV           ! number of scalar var.
 INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
@@ -339,8 +336,8 @@ INTEGER             :: ILENCH       ! Length of comment string in LFIFM file
 INTEGER             :: IKB,IKE      ! I index values for the Beginning and End
 INTEGER             :: IKU  ! array sizes
 
-REAL, DIMENSION(D%NIT,D%NJT,MIN(D%NKA+JPVEXT_TURB*KKL,D%NKA+JPVEXT_TURB*KKL+2*KKL):&
-                            MAX(D%NKA+JPVEXT_TURB*KKL,D%NKA+JPVEXT_TURB*KKL+2*KKL))&
+REAL, DIMENSION(D%NIT,D%NJT,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
@@ -370,16 +367,16 @@ GUSERV = (KRR/=0)
 !
 !  compute the coefficients for the uncentred gradient computation near the 
 !  ground
-ZCOEFF(:,:,IKB+2*KKL)= - PDZZ(:,:,IKB+KKL) /      &
-       ( (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) * PDZZ(:,:,IKB+2*KKL) )
-ZCOEFF(:,:,IKB+KKL)=   (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) /      &
-       ( PDZZ(:,:,IKB+KKL) * PDZZ(:,:,IKB+2*KKL) )
-ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2*KKL)+2.*PDZZ(:,:,IKB+KKL)) /      &
-       ( (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) * PDZZ(:,:,IKB+KKL) )
+ZCOEFF(:,:,IKB+2*D%NKL)= - PDZZ(:,:,IKB+D%NKL) /      &
+       ( (PDZZ(:,:,IKB+2*D%NKL)+PDZZ(:,:,IKB+D%NKL)) * PDZZ(:,:,IKB+2*D%NKL) )
+ZCOEFF(:,:,IKB+D%NKL)=   (PDZZ(:,:,IKB+2*D%NKL)+PDZZ(:,:,IKB+D%NKL)) /      &
+       ( PDZZ(:,:,IKB+D%NKL) * PDZZ(:,:,IKB+2*D%NKL) )
+ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2*D%NKL)+2.*PDZZ(:,:,IKB+D%NKL)) /      &
+       ( (PDZZ(:,:,IKB+2*D%NKL)+PDZZ(:,:,IKB+D%NKL)) * PDZZ(:,:,IKB+D%NKL) )
 !
 !
 IF (OHARAT) THEN
-PLMF=MZF(PLM, D%NKA, KKU, KKL)
+PLMF=MZF(PLM, D%NKA, D%NKU, D%NKL)
 PLEPSF=PLMF
 !  function MZF produces -999 for level IKU (82 for 80 levels)
 !  so put these to normal value as this level (82) is indeed calculated
@@ -387,7 +384,7 @@ PLMF(:,:,D%NKT)=0.001
 PLEPSF(:,:,D%NKT)=0.001
 ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
 ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, KKU, KKL)
+ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, D%NKU, D%NKL)
 ENDIF
 !
 
@@ -419,9 +416,9 @@ END IF
 ! Compute the turbulent variance F and F' at time t-dt.
 !
 IF (OHARAT) THEN
-  ZF      (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ**2, D%NKA, KKU, KKL)
+  ZF      (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ**2, D%NKA, D%NKU, D%NKL)
 ELSE
-  ZF      (:,:,:) = CSTURB%XCTV*PLM*PLEPS*MZF(PPHI3*PDTH_DZ**2, D%NKA, KKU, KKL)
+  ZF      (:,:,:) = CSTURB%XCTV*PLM*PLEPS*MZF(PPHI3*PDTH_DZ**2, D%NKA, D%NKU, D%NKL)
 ENDIF
   ZDFDDTDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
   !
@@ -429,42 +426,42 @@ ENDIF
   !
   ! d(w'th'2)/dz
   IF (GFTH2) THEN
-    ZF       = ZF       + M3_TH2_WTH2(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,PLEPS,&
+    ZF       = ZF       + M3_TH2_WTH2(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,PLEPS,&
      & PSQRT_TKE) * PFTH2
-    ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WTH2_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+    ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WTH2_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
      & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA) * PFTH2
   END IF
   !
   ! d(w'2th')/dz
   IF (GFWTH) THEN
-    ZF       = ZF       + M3_TH2_W2TH(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,PDTH_DZ,&
-     & PLM,PLEPS,PTKEM) * MZF(PFWTH, D%NKA, KKU, KKL)
-    ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_W2TH_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,&
-     & PLM,PLEPS,PTKEM,GUSERV) * MZF(PFWTH, D%NKA, KKU, KKL)
+    ZF       = ZF       + M3_TH2_W2TH(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,PDTH_DZ,&
+     & PLM,PLEPS,PTKEM) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
+    ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_W2TH_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,&
+     & PLM,PLEPS,PTKEM,GUSERV) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
   END IF
   !
   IF (KRR/=0) THEN
     ! d(w'r'2)/dz
     IF (GFR2) THEN
-      ZF       = ZF       + M3_TH2_WR2(D,CSTURB,D%NKA,KKU,KKL,PD,PLEPS,PSQRT_TKE,PBLL_O_E,&
+      ZF       = ZF       + M3_TH2_WR2(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,PLEPS,PSQRT_TKE,PBLL_O_E,&
        & PEMOIST,PDTH_DZ) * PFR2
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WR2_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,&
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WR2_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,&
        & PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTH_DZ) * PFR2
     END IF
     !
     ! d(w'2r')/dz
     IF (GFWR) THEN
-      ZF       = ZF       + M3_TH2_W2R(D,CSTURB,D%NKA,KKU,KKL,PD,PLM,PLEPS,PTKEM,PBLL_O_E,&
-       & PEMOIST,PDTH_DZ) * MZF(PFWR, D%NKA, KKU, KKL)
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_W2R_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,&
-       & PLM,PLEPS,PTKEM,PBLL_O_E,PEMOIST,PDTH_DZ) * MZF(PFWR, D%NKA, KKU, KKL)
+      ZF       = ZF       + M3_TH2_W2R(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,PLM,PLEPS,PTKEM,PBLL_O_E,&
+       & PEMOIST,PDTH_DZ) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_W2R_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,&
+       & PLM,PLEPS,PTKEM,PBLL_O_E,PEMOIST,PDTH_DZ) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
     END IF
     !
     ! d(w'th'r')/dz
     IF (GFTHR) THEN
-      ZF       = ZF       + M3_TH2_WTHR(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PD,PLEPS,PSQRT_TKE,&
+      ZF       = ZF       + M3_TH2_WTHR(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PD,PLEPS,PSQRT_TKE,&
        & PBLL_O_E,PEMOIST,PDTH_DZ) * PFTHR
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WTHR_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WTHR_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
        & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTH_DZ) * PFTHR
     END IF
 
@@ -473,32 +470,32 @@ ENDIF
   ZFLXZ(:,:,:)   = ZF                                                              &
   !     + PIMPL * CSTURB%XCTV*PLM*PLEPS                                                   &
   !        *MZF(D_PHI3DTDZ2_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,PDTH_DZ,HTURBDIM,GUSERV)   &
-  !             *DZM(PTHLP - PTHLM, D%NKA, KKU, KKL) / PDZZ                                        ) &
-        + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM, D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)
+  !             *DZM(PTHLP - PTHLM, D%NKA, D%NKU, D%NKL) / PDZZ                                        ) &
+        + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM, D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)
   !
   ! special case near the ground ( uncentred gradient )
   IF (OHARAT) THEN
   ZFLXZ(:,:,IKB) =  PLMF(:,:,IKB)   &
      * PLEPSF(:,:,IKB)                                         &
   *( PEXPL *                                                  &
-     ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL)             &
-      +ZCOEFF(:,:,IKB+KKL  )*PTHLM(:,:,IKB+KKL  )             & 
+     ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLM(:,:,IKB+2*D%NKL)             &
+      +ZCOEFF(:,:,IKB+D%NKL  )*PTHLM(:,:,IKB+D%NKL  )             & 
       +ZCOEFF(:,:,IKB      )*PTHLM(:,:,IKB  )   )**2          &
     +PIMPL *                                                  &
-     ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL)             &
-      +ZCOEFF(:,:,IKB+KKL  )*PTHLP(:,:,IKB+KKL  )             &
+     ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLP(:,:,IKB+2*D%NKL)             &
+      +ZCOEFF(:,:,IKB+D%NKL  )*PTHLP(:,:,IKB+D%NKL  )             &
       +ZCOEFF(:,:,IKB      )*PTHLP(:,:,IKB  )   )**2          &
    ) 
    ELSE
-  ZFLXZ(:,:,IKB) = CSTURB%XCTV * PPHI3(:,:,IKB+KKL) * PLM(:,:,IKB)   &
+  ZFLXZ(:,:,IKB) = CSTURB%XCTV * PPHI3(:,:,IKB+D%NKL) * PLM(:,:,IKB)   &
      * PLEPS(:,:,IKB)                                         &
   *( PEXPL *                                                  &
-     ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL)             &
-      +ZCOEFF(:,:,IKB+KKL  )*PTHLM(:,:,IKB+KKL  )             & 
+     ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLM(:,:,IKB+2*D%NKL)             &
+      +ZCOEFF(:,:,IKB+D%NKL  )*PTHLM(:,:,IKB+D%NKL  )             & 
       +ZCOEFF(:,:,IKB      )*PTHLM(:,:,IKB  )   )**2          &
     +PIMPL *                                                  &
-     ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL)             &
-      +ZCOEFF(:,:,IKB+KKL  )*PTHLP(:,:,IKB+KKL  )             &
+     ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLP(:,:,IKB+2*D%NKL)             &
+      +ZCOEFF(:,:,IKB+D%NKL  )*PTHLP(:,:,IKB+D%NKL  )             &
       +ZCOEFF(:,:,IKB      )*PTHLP(:,:,IKB  )   )**2          &
    ) 
    ENDIF
@@ -532,7 +529,7 @@ ENDIF
   IF (OLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
     CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Thl2 ) 
-    CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Thl2 )
+    CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, D%NKU, D%NKL)*ZFLXZ, X_LES_RES_W_SBG_Thl2 )
     CALL LES_MEAN_SUBGRID(-2.*CSTURB%XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Thl2 ) 
     CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_ThlThv ) 
     CALL LES_MEAN_SUBGRID(-CSTURB%XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. ) 
@@ -547,9 +544,9 @@ ENDIF
 !
     ! Compute the turbulent variance F and F' at time t-dt.
 IF (OHARAT) THEN
-    ZF      (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ*PDR_DZ, D%NKA, KKU, KKL)
+    ZF      (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ*PDR_DZ, D%NKA, D%NKU, D%NKL)
 ELSE
-    ZF      (:,:,:) = CSTURB%XCTV*PLM*PLEPS*MZF(0.5*(PPHI3+PPSI3)*PDTH_DZ*PDR_DZ, D%NKA, KKU, KKL)
+    ZF      (:,:,:) = CSTURB%XCTV*PLM*PLEPS*MZF(0.5*(PPHI3+PPSI3)*PDTH_DZ*PDR_DZ, D%NKA, D%NKU, D%NKL)
 ENDIF
     ZDFDDTDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
     ZDFDDRDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
@@ -558,51 +555,51 @@ ENDIF
     !
     ! d(w'th'2)/dz
     IF (GFTH2) THEN
-      ZF       = ZF       + M3_THR_WTH2(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PD,PLEPS,PSQRT_TKE,&
+      ZF       = ZF       + M3_THR_WTH2(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PD,PLEPS,PSQRT_TKE,&
        & PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WTH2_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WTH2_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
        & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WTH2_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WTH2_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
        & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA) * PFTH2
     END IF
     !
     ! d(w'2th')/dz
     IF (GFWTH) THEN
-      ZF       = ZF       + M3_THR_W2TH(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PD,PLM,PLEPS,PTKEM,&
-       & PDR_DZ) * MZF(PFWTH, D%NKA, KKU, KKL)
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_W2TH_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
-       & PD,PLM,PLEPS,PTKEM,PBLL_O_E,PDR_DZ,PETHETA) * MZF(PFWTH, D%NKA, KKU, KKL)
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_W2TH_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
-       & PD,PLM,PLEPS,PTKEM) * MZF(PFWTH, D%NKA, KKU, KKL)
+      ZF       = ZF       + M3_THR_W2TH(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PD,PLM,PLEPS,PTKEM,&
+       & PDR_DZ) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_W2TH_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
+       & PD,PLM,PLEPS,PTKEM,PBLL_O_E,PDR_DZ,PETHETA) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_W2TH_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
+       & PD,PLM,PLEPS,PTKEM) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
     END IF
     !
     ! d(w'r'2)/dz
     IF (GFR2) THEN
-      ZF       = ZF       + M3_THR_WR2(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PD,PLEPS,PSQRT_TKE,&
+      ZF       = ZF       + M3_THR_WR2(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PD,PLEPS,PSQRT_TKE,&
        & PBLL_O_E,PEMOIST,PDTH_DZ) * PFR2
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WR2_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WR2_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,&
        & PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST) * PFR2
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WR2_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WR2_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,&
        & PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTH_DZ) * PFR2
     END IF
     !
       ! d(w'2r')/dz
     IF (GFWR) THEN
-      ZF       = ZF       + M3_THR_W2R(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PD,PLM,PLEPS,PTKEM,&
-      & PDTH_DZ) * MZF(PFWR, D%NKA, KKU, KKL)
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_W2R_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,&
-      & PLM,PLEPS,PTKEM) * MZF(PFWR, D%NKA, KKU, KKL)
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_W2R_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,&
-      & PLM,PLEPS,PTKEM,PBLL_O_E,PDTH_DZ,PEMOIST) * MZF(PFWR, D%NKA, KKU, KKL)
+      ZF       = ZF       + M3_THR_W2R(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PD,PLM,PLEPS,PTKEM,&
+      & PDTH_DZ) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_W2R_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,&
+      & PLM,PLEPS,PTKEM) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_W2R_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,&
+      & PLM,PLEPS,PTKEM,PBLL_O_E,PDTH_DZ,PEMOIST) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
     END IF
     !
     ! d(w'th'r')/dz
     IF (GFTHR) THEN
-      ZF       = ZF       + M3_THR_WTHR(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,PLEPS,&
+      ZF       = ZF       + M3_THR_WTHR(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,PLEPS,&
       & PSQRT_TKE) * PFTHR
-      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WTHR_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+      ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WTHR_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
       & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA) * PFTHR
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WTHR_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,&
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WTHR_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,&
       & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST) * PFTHR
     END IF
     !
@@ -610,24 +607,24 @@ ENDIF
     ZFLXZ(:,:,:)   = ZF                                                     &
         + PIMPL * PLMF*PLEPSF*0.5                                        &
           * MZF(( 2.  & 
-                 ) *PDR_DZ  *DZM(PTHLP - PTHLM, D%NKA, KKU, KKL    ) / PDZZ               &
+                 ) *PDR_DZ  *DZM(PTHLP - PTHLM, D%NKA, D%NKU, D%NKL    ) / PDZZ               &
                 +( 2.                                                    &
-                 ) *PDTH_DZ *DZM(PRP   - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ               &
-               , D%NKA, KKU, KKL)                                                            &
-        + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)         &
-        + PIMPL * ZDFDDRDZ * MZF(DZM(PRP   - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)
+                 ) *PDTH_DZ *DZM(PRP   - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ               &
+               , D%NKA, D%NKU, D%NKL)                                                            &
+        + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)         &
+        + PIMPL * ZDFDDRDZ * MZF(DZM(PRP   - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)
     ELSE
     ZFLXZ(:,:,:)   = ZF                                                     &
         + PIMPL * CSTURB%XCTV*PLM*PLEPS*0.5                                        &
           * MZF(( D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*dthdz)/ddthdz term
                   +D_PSI3DTDZ_O_DDTDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*dthdz)/ddthdz term
-                 ) *PDR_DZ  *DZM(PTHLP - PTHLM, D%NKA, KKU, KKL) / PDZZ               &
+                 ) *PDR_DZ  *DZM(PTHLP - PTHLM, D%NKA, D%NKU, D%NKL) / PDZZ               &
                 +( D_PHI3DRDZ_O_DDRDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*drdz )/ddrdz term
                   +D_PSI3DRDZ_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*drdz )/ddrdz term
-                 ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ               &
-               , D%NKA, KKU, KKL)                                                            &
-        + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)         &
-        + PIMPL * ZDFDDRDZ * MZF(DZM(PRP   - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)
+                 ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ               &
+               , D%NKA, D%NKU, D%NKL)                                                            &
+        + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)         &
+        + PIMPL * ZDFDDRDZ * MZF(DZM(PRP   - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)
     ENDIF
     !
     ! special case near the ground ( uncentred gradient )
@@ -635,36 +632,36 @@ ENDIF
     ZFLXZ(:,:,IKB) =                                            & 
     (1. )   &
     *( PEXPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PTHLM(:,:,IKB+KKL  )             & 
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLM(:,:,IKB+2*D%NKL)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PTHLM(:,:,IKB+D%NKL  )             & 
         +ZCOEFF(:,:,IKB      )*PTHLM(:,:,IKB      ))            &
-      *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PRM(:,:,IKB+KKL,1  )             & 
+      *( ZCOEFF(:,:,IKB+2*D%NKL)*PRM(:,:,IKB+2*D%NKL,1)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRM(:,:,IKB+D%NKL,1  )             & 
         +ZCOEFF(:,:,IKB      )*PRM(:,:,IKB  ,1    ))            &
       +PIMPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PTHLP(:,:,IKB+KKL  )             &
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLP(:,:,IKB+2*D%NKL)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PTHLP(:,:,IKB+D%NKL  )             &
         +ZCOEFF(:,:,IKB      )*PTHLP(:,:,IKB      ))            &
-      *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL  )             &
-        +ZCOEFF(:,:,IKB+KKL  )*PRP(:,:,IKB+KKL    )             & 
+      *( ZCOEFF(:,:,IKB+2*D%NKL)*PRP(:,:,IKB+2*D%NKL  )             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRP(:,:,IKB+D%NKL    )             & 
         +ZCOEFF(:,:,IKB      )*PRP(:,:,IKB        ))            &
      ) 
     ELSE 
     ZFLXZ(:,:,IKB) =                                            & 
-    (CSTURB%XCHT1 * PPHI3(:,:,IKB+KKL) + CSTURB%XCHT2 * PPSI3(:,:,IKB+KKL))   &
+    (CSTURB%XCHT1 * PPHI3(:,:,IKB+D%NKL) + CSTURB%XCHT2 * PPSI3(:,:,IKB+D%NKL))   &
     *( PEXPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PTHLM(:,:,IKB+KKL  )             & 
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLM(:,:,IKB+2*D%NKL)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PTHLM(:,:,IKB+D%NKL  )             & 
         +ZCOEFF(:,:,IKB      )*PTHLM(:,:,IKB      ))            &
-      *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PRM(:,:,IKB+KKL,1  )             & 
+      *( ZCOEFF(:,:,IKB+2*D%NKL)*PRM(:,:,IKB+2*D%NKL,1)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRM(:,:,IKB+D%NKL,1  )             & 
         +ZCOEFF(:,:,IKB      )*PRM(:,:,IKB  ,1    ))            &
       +PIMPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PTHLP(:,:,IKB+KKL  )             &
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PTHLP(:,:,IKB+2*D%NKL)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PTHLP(:,:,IKB+D%NKL  )             &
         +ZCOEFF(:,:,IKB      )*PTHLP(:,:,IKB      ))            &
-      *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL  )             &
-        +ZCOEFF(:,:,IKB+KKL  )*PRP(:,:,IKB+KKL    )             & 
+      *( ZCOEFF(:,:,IKB+2*D%NKL)*PRP(:,:,IKB+2*D%NKL  )             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRP(:,:,IKB+D%NKL    )             & 
         +ZCOEFF(:,:,IKB      )*PRP(:,:,IKB        ))            &
      ) 
     ENDIF
@@ -699,7 +696,7 @@ ENDIF
 IF (OLES_CALL) THEN
       CALL SECOND_MNH(ZTIME1)
       CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_THlRt ) 
-      CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_ThlRt )
+      CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, D%NKU, D%NKL)*ZFLXZ, X_LES_RES_W_SBG_ThlRt )
       CALL LES_MEAN_SUBGRID(-2.*CSTURB%XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_ThlRt ) 
       CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_RtThv ) 
       CALL LES_MEAN_SUBGRID(-CSTURB%XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. ) 
@@ -715,9 +712,9 @@ END IF
 !
     ! Compute the turbulent variance F and F' at time t-dt.
 IF (OHARAT) THEN
-    ZF      (:,:,:) = PLMF*PLEPSF*MZF(PDR_DZ**2, D%NKA, KKU, KKL)
+    ZF      (:,:,:) = PLMF*PLEPSF*MZF(PDR_DZ**2, D%NKA, D%NKU, D%NKL)
   ELSE
-    ZF      (:,:,:) = CSTURB%XCTV*PLM*PLEPS*MZF(PPSI3*PDR_DZ**2, D%NKA, KKU, KKL)
+    ZF      (:,:,:) = CSTURB%XCTV*PLM*PLEPS*MZF(PPSI3*PDR_DZ**2, D%NKA, D%NKU, D%NKL)
 ENDIF
     ZDFDDRDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
     !
@@ -725,42 +722,42 @@ ENDIF
     !
     ! d(w'r'2)/dz
     IF (GFR2) THEN
-      ZF       = ZF       + M3_R2_WR2(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,PLEPS,&
+      ZF       = ZF       + M3_R2_WR2(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,PLEPS,&
       & PSQRT_TKE) * PFR2
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WR2_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,&
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WR2_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,&
       & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST) * PFR2
     END IF
     !
     ! d(w'2r')/dz
     IF (GFWR) THEN
-      ZF       = ZF       + M3_R2_W2R(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,PDR_DZ,&
-      & PLM,PLEPS,PTKEM) * MZF(PFWR, D%NKA, KKU, KKL)
-      ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_W2R_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,&
-      & PD,PLM,PLEPS,PTKEM,GUSERV) * MZF(PFWR, D%NKA, KKU, KKL)
+      ZF       = ZF       + M3_R2_W2R(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,PDR_DZ,&
+      & PLM,PLEPS,PTKEM) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
+      ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_W2R_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,&
+      & PD,PLM,PLEPS,PTKEM,GUSERV) * MZF(PFWR, D%NKA, D%NKU, D%NKL)
     END IF
     !
     IF (KRR/=0) THEN
       ! d(w'r'2)/dz
       IF (GFTH2) THEN
-        ZF       = ZF       + M3_R2_WTH2(D,CSTURB,D%NKA,KKU,KKL,PD,PLEPS,PSQRT_TKE,&
+        ZF       = ZF       + M3_R2_WTH2(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,PLEPS,PSQRT_TKE,&
         & PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
-        ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTH2_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,&
+        ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTH2_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,&
         & PREDTH1,PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
       END IF
       !
       ! d(w'2r')/dz
       IF (GFWTH) THEN
-        ZF       = ZF       + M3_R2_W2TH(D,CSTURB,D%NKA,KKU,KKL,PD,PLM,PLEPS,PTKEM,&
-        & PBLL_O_E,PETHETA,PDR_DZ) * MZF(PFWTH, D%NKA, KKU, KKL)
-        ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_W2TH_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,&
-        & PD,PLM,PLEPS,PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * MZF(PFWTH, D%NKA, KKU, KKL)
+        ZF       = ZF       + M3_R2_W2TH(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,PLM,PLEPS,PTKEM,&
+        & PBLL_O_E,PETHETA,PDR_DZ) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
+        ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_W2TH_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,&
+        & PD,PLM,PLEPS,PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * MZF(PFWTH, D%NKA, D%NKU, D%NKL)
       END IF
       !
       ! d(w'th'r')/dz
       IF (GFTHR) THEN
-        ZF       = ZF       + M3_R2_WTHR(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PD,PLEPS,&
+        ZF       = ZF       + M3_R2_WTHR(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PD,PLEPS,&
         & PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTHR
-        ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTHR_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,&
+        ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTHR_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,&
         & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTHR
       END IF
   
@@ -770,14 +767,14 @@ ENDIF
     ZFLXZ(:,:,:)   = ZF                                                              &
           + PIMPL * PLMF*PLEPSF                                                  &
             *MZF(2.*PDR_DZ*    &
-                 DZM(PRP - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL) &
-          + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)
+                 DZM(PRP - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL) &
+          + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)
    ELSE
     ZFLXZ(:,:,:)   = ZF                                                              &
           + PIMPL * CSTURB%XCTV*PLM*PLEPS                                                  &
             *MZF(D_PSI3DRDZ2_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,PDR_DZ,HTURBDIM,GUSERV)    &
-                 *DZM(PRP - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL) &
-          + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ, D%NKA, KKU, KKL)
+                 *DZM(PRP - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL) &
+          + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ, D%NKA, D%NKU, D%NKL)
   ENDIF
     !
     ! special case near the ground ( uncentred gradient )
@@ -785,24 +782,24 @@ ENDIF
     ZFLXZ(:,:,IKB) =  PLMF(:,:,IKB)   &
         * PLEPSF(:,:,IKB)                                        &
     *( PEXPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PRM(:,:,IKB+KKL,1  )             & 
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PRM(:,:,IKB+2*D%NKL,1)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRM(:,:,IKB+D%NKL,1  )             & 
         +ZCOEFF(:,:,IKB      )*PRM(:,:,IKB  ,1    ))**2         &
       +PIMPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL)               &
-        +ZCOEFF(:,:,IKB+KKL  )*PRP(:,:,IKB+KKL  )               &
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PRP(:,:,IKB+2*D%NKL)               &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRP(:,:,IKB+D%NKL  )               &
         +ZCOEFF(:,:,IKB      )*PRP(:,:,IKB      ))**2           &
      ) 
    ELSE
-    ZFLXZ(:,:,IKB) = CSTURB%XCHV * PPSI3(:,:,IKB+KKL) * PLM(:,:,IKB)   &
+    ZFLXZ(:,:,IKB) = CSTURB%XCHV * PPSI3(:,:,IKB+D%NKL) * PLM(:,:,IKB)   &
         * PLEPS(:,:,IKB)                                        &
     *( PEXPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1)             &
-        +ZCOEFF(:,:,IKB+KKL  )*PRM(:,:,IKB+KKL,1  )             & 
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PRM(:,:,IKB+2*D%NKL,1)             &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRM(:,:,IKB+D%NKL,1  )             & 
         +ZCOEFF(:,:,IKB      )*PRM(:,:,IKB  ,1    ))**2         &
       +PIMPL *                                                  &
-       ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL)               &
-        +ZCOEFF(:,:,IKB+KKL  )*PRP(:,:,IKB+KKL  )               &
+       ( ZCOEFF(:,:,IKB+2*D%NKL)*PRP(:,:,IKB+2*D%NKL)               &
+        +ZCOEFF(:,:,IKB+D%NKL  )*PRP(:,:,IKB+D%NKL  )               &
         +ZCOEFF(:,:,IKB      )*PRP(:,:,IKB      ))**2           &
      ) 
   ENDIF
@@ -832,7 +829,7 @@ ENDIF
     IF (OLES_CALL) THEN
       CALL SECOND_MNH(ZTIME1)
       CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Rt2 ) 
-      CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Rt2 )
+      CALL LES_MEAN_SUBGRID(MZF(PWM, D%NKA, D%NKU, D%NKL)*ZFLXZ, X_LES_RES_W_SBG_Rt2 )
       CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_RtThv , .TRUE. ) 
       CALL LES_MEAN_SUBGRID(-CSTURB%XA3*PBETA*PEMOIST*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. )
       CALL LES_MEAN_SUBGRID(-2.*CSTURB%XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Rt2 ) 
@@ -848,7 +845,7 @@ ENDIF
   IF ( KRRL > 0 ) THEN
     ! Extrapolate PSIGS at the ground and at the top
     PSIGS(:,:,D%NKA) = PSIGS(:,:,IKB)
-    PSIGS(:,:,KKU) = PSIGS(:,:,IKE)
+    PSIGS(:,:,D%NKU) = PSIGS(:,:,IKE)
 #ifdef REPRO48
     PSIGS(:,:,:) =  MAX (PSIGS(:,:,:) , 0.)
     PSIGS(:,:,:) =  SQRT(PSIGS(:,:,:))
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index 5f8e070d8..b11ec118d 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -7,7 +7,7 @@ IMPLICIT NONE
 CONTAINS
       
 SUBROUTINE TURB_VER_THERMO_FLUX(D,CST,CSTURB,TURBN,                 &
-                      KKA,KKU,KKL,KRR,KRRL,KRRI,KSV,                &
+                      KRR,KRRL,KRRI,KSV,                            &
                       OTURB_FLX,HTURBDIM,HTOM,OOCEAN,ODEEPOC,OHARAT,&
                       OCOUPLES,OLES_CALL, OCOMPUTE_SRC,             &
                       PIMPL,PEXPL,PTSTEP,HPROGRAM,                  &
@@ -269,9 +269,6 @@ TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
 TYPE(CST_t),            INTENT(IN)   :: CST
 TYPE(CSTURB_t),         INTENT(IN)   :: CSTURB
 TYPE(TURB_t),           INTENT(IN)   :: TURBN
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
 INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
 INTEGER,                INTENT(IN)   :: KSV           ! number of scalar var.
 INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
@@ -394,8 +391,7 @@ INTEGER                    :: IIB,IJB       ! Lower bounds of the physical
 INTEGER                    :: IIE,IJE       ! Upper bounds of the physical
                                             ! sub-domain in x and y directions
 !
-! TODO TO BE REMOVED OUTSIDE OF TURB ? : SHOULD BE NI/JMAX_ll+2*JPHEXT
-! WARNING WRONG DIMENSION, NIMAX_ll not initialized with AROME
+! NIMPORTE QUOI : TODO TO BE REMOVED OUTSIDE OF TURB ? :
 REAL, DIMENSION(1)  :: ZXHAT_ll  !  Position x in the conformal
                                                  ! plane (array on the complete domain)
 REAL, DIMENSION(1)  :: ZYHAT_ll  !   Position y in the conformal
@@ -475,7 +471,7 @@ IF (OHARAT) THEN
 ! OHARAT so TKE and length scales at half levels!
   ZKEFF(:,:,:) =  PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
 ELSE
-  ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, KKU, KKL)
+  ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, D%NKU, D%NKL)
 ENDIF
 !
 ! Define a cloud mask with ri and rc (used after with a threshold) for Leonard terms
@@ -516,10 +512,10 @@ END IF
 ! Compute the turbulent flux F and F' at time t-dt.
 !
 IF (OHARAT) THEN
-  ZF      (:,:,:) = -ZKEFF*DZM(PTHLM, D%NKA, KKU, KKL)/PDZZ
+  ZF      (:,:,:) = -ZKEFF*DZM(PTHLM, D%NKA, D%NKU, D%NKL)/PDZZ
   ZDFDDTDZ(:,:,:) = -ZKEFF
 ELSE
-  ZF      (:,:,:) = -CSTURB%XCSHF*PPHI3*ZKEFF*DZM(PTHLM, D%NKA, KKU, KKL)/PDZZ
+  ZF      (:,:,:) = -CSTURB%XCSHF*PPHI3*ZKEFF*DZM(PTHLM, D%NKA, D%NKU, D%NKL)/PDZZ
   ZDFDDTDZ(:,:,:) = -CSTURB%XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
 END IF
 !
@@ -527,20 +523,20 @@ IF (TURBN%LHGRAD) THEN
  ! Compute the Leonard terms for thl
  ZDELTAX= XXHAT(3) - XXHAT(2)
  ZF_LEONARD (:,:,:)= TURBN%XCOEFHGRADTHL*ZDELTAX*ZDELTAX/12.0*(      &
-                 MXF(GX_W_UW(PWM(:,:,:), XDXX,XDZZ,XDZX,D%NKA,KKU,KKL))&
-                *MZM(GX_M_M(PTHLM(:,:,:),XDXX,XDZZ,XDZX,D%NKA, KKU, KKL), D%NKA, KKU, KKL)  &
-              +  MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY,D%NKA,KKU,KKL))  &
-                *MZM(GY_M_M(PTHLM(:,:,:),XDYY,XDZZ,XDZY,D%NKA, KKU, KKL), D%NKA, KKU, KKL) )
+                 MXF(GX_W_UW(PWM(:,:,:), XDXX,XDZZ,XDZX,D%NKA,D%NKU,D%NKL))&
+                *MZM(GX_M_M(PTHLM(:,:,:),XDXX,XDZZ,XDZX,D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)  &
+              +  MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY,D%NKA,D%NKU,D%NKL))  &
+                *MZM(GY_M_M(PTHLM(:,:,:),XDYY,XDZZ,XDZY,D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL) )
 END IF
 !
 ! Effect of 3rd order terms in temperature flux (at flux point)
 !
 ! d(w'2th')/dz
 IF (GFWTH) THEN
-  Z3RDMOMENT= M3_WTH_W2TH(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,ZKEFF,PTKEM)
+  Z3RDMOMENT= M3_WTH_W2TH(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,ZKEFF,PTKEM)
 !
   ZF       = ZF       + Z3RDMOMENT * PFWTH
-  ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2TH_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+  ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2TH_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
    & PD,PBLL_O_E,PETHETA,ZKEFF,PTKEM) * PFWTH
 END IF
 !
@@ -548,41 +544,41 @@ END IF
 IF (GFTH2) THEN
   Z3RDMOMENT= M3_WTH_WTH2(D,CSTURB,PREDTH1,PREDR1,PD,PBLL_O_E,PETHETA)
 !
-  ZF       = ZF       + Z3RDMOMENT * MZM(PFTH2, D%NKA, KKU, KKL)
+  ZF       = ZF       + Z3RDMOMENT * MZM(PFTH2, D%NKA, D%NKU, D%NKL)
   ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTH2_O_DDTDZ(D,CSTURB,Z3RDMOMENT,PREDTH1,PREDR1,&
-    & PD,PBLL_O_E,PETHETA) * MZM(PFTH2, D%NKA, KKU, KKL)
+    & PD,PBLL_O_E,PETHETA) * MZM(PFTH2, D%NKA, D%NKU, D%NKL)
 END IF
 !
 ! d(w'2r')/dz
 IF (GFWR) THEN
-  ZF       = ZF       + M3_WTH_W2R(D,CSTURB,D%NKA,KKU,KKL,PD,ZKEFF,&
+  ZF       = ZF       + M3_WTH_W2R(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,ZKEFF,&
     & PTKEM,PBLL_O_E,PEMOIST,PDTH_DZ) * PFWR
-  ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2R_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,&
+  ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2R_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,&
     & PD,ZKEFF,PTKEM,PBLL_O_E,PEMOIST) * PFWR
 END IF
 !
 ! d(w'r'2)/dz
 IF (GFR2) THEN
-  ZF       = ZF       + M3_WTH_WR2(D,CSTURB,D%NKA,KKU,KKL,PD,ZKEFF,PTKEM,&
-    & PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST,PDTH_DZ) * MZM(PFR2, D%NKA, KKU, KKL)
-  ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WR2_O_DDTDZ(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PREDR1,PD,&
-    & ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST) * MZM(PFR2, D%NKA, KKU, KKL)
+  ZF       = ZF       + M3_WTH_WR2(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,ZKEFF,PTKEM,&
+    & PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST,PDTH_DZ) * MZM(PFR2, D%NKA, D%NKU, D%NKL)
+  ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WR2_O_DDTDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PREDR1,PD,&
+    & ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST) * MZM(PFR2, D%NKA, D%NKU, D%NKL)
 END IF
 !
 ! d(w'th'r')/dz
 IF (GFTHR) THEN
-  Z3RDMOMENT= M3_WTH_WTHR(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PD,ZKEFF,PTKEM,PSQRT_TKE,PBETA,&
+  Z3RDMOMENT= M3_WTH_WTHR(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PD,ZKEFF,PTKEM,PSQRT_TKE,PBETA,&
     & PLEPS,PEMOIST)
 !
-  ZF       = ZF       + Z3RDMOMENT * MZM(PFTHR, D%NKA, KKU, KKL)
+  ZF       = ZF       + Z3RDMOMENT * MZM(PFTHR, D%NKA, D%NKU, D%NKL)
   ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTHR_O_DDTDZ(D,CSTURB,Z3RDMOMENT,PREDTH1,&
-    & PREDR1,PD,PBLL_O_E,PETHETA) * MZM(PFTHR, D%NKA, KKU, KKL)
+    & PREDR1,PD,PBLL_O_E,PETHETA) * MZM(PFTHR, D%NKA, D%NKU, D%NKL)
 END IF
 ! compute interface flux
 IF (OCOUPLES) THEN   ! Autocoupling O-A LES
   IF (OOCEAN) THEN    ! ocean model in coupled case 
     ZF(:,:,IKE) =  (TURBN%XSSTFL_C(:,:,1)+TURBN%XSSRFL_C(:,:,1)) &
-                  *0.5* ( 1. + PRHODJ(:,:,KKU)/PRHODJ(:,:,IKE) )
+                  *0.5* ( 1. + PRHODJ(:,:,D%NKU)/PRHODJ(:,:,IKE) )
   ELSE                ! atmosph model in coupled case
     ZF(:,:,IKB) =  TURBN%XSSTFL_C(:,:,1) &
                   *0.5* ( 1. + PRHODJ(:,:,D%NKA)/PRHODJ(:,:,IKB) )
@@ -604,7 +600,7 @@ ELSE  ! No coupling O and A cases
   END IF
 !
   IF (OOCEAN) THEN
-    ZF(:,:,IKE) = XSSTFL(:,:) *0.5*(1. + PRHODJ(:,:,KKU) / PRHODJ(:,:,IKE))
+    ZF(:,:,IKE) = XSSTFL(:,:) *0.5*(1. + PRHODJ(:,:,D%NKU) / PRHODJ(:,:,IKE))
   ELSE !end ocean case (in nocoupled case)
     ! atmos top
 #ifdef REPRO48
@@ -615,7 +611,7 @@ ELSE  ! No coupling O and A cases
 END IF !end no coupled cases
 !
 ! Compute the split conservative potential temperature at t+deltat
-CALL TRIDIAG_THERMO(D,D%NKA,KKU,KKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
+CALL TRIDIAG_THERMO(D,D%NKA,D%NKU,D%NKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
                     PRHODJ,PTHLP)
 !
 ! Compute the equivalent tendency for the conservative potential temperature
@@ -623,12 +619,12 @@ CALL TRIDIAG_THERMO(D,D%NKA,KKU,KKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
 ZRWTHL(:,:,:)= PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
 ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
 IF (TURBN%LHGRAD) THEN
- DO JK=1,KKU
+ DO JK=1,D%NKU
   ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
  END DO
  WHERE ( (ZCLD_THOLD(:,:,:) >= TURBN%XCLDTHOLD) .AND. ( ZALT(:,:,:) >= TURBN%XALTHGRAD) )
-  ZRWTHL(:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:), D%NKA, KKU, KKL)*ZF_LEONARD(:,:,:),XDZZ,&
-                   D%NKA, KKU, KKL)
+  ZRWTHL(:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:), D%NKA, D%NKU, D%NKL)*ZF_LEONARD(:,:,:),XDZZ,&
+                   D%NKA, D%NKU, D%NKL)
  END WHERE
 END IF
 !
@@ -639,7 +635,7 @@ PRTHLS(:,:,:)= PRTHLS(:,:,:)  + ZRWTHL(:,:,:)
 !  Conservative potential temperature flux : 
 !
 ZFLXZ(:,:,:)   = ZF                                                &
-               + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, D%NKA, KKU, KKL) / PDZZ
+               + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, D%NKA, D%NKU, D%NKL) / PDZZ
 ! replace the flux by the Leonard terms
 IF (TURBN%LHGRAD) THEN
  WHERE ( (ZCLD_THOLD(:,:,:) >= TURBN%XCLDTHOLD) .AND. ( ZALT(:,:,:) >= TURBN%XALTHGRAD) )
@@ -649,22 +645,22 @@ END IF
 !
 ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB) 
 IF (OOCEAN) THEN
-  ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+  ZFLXZ(:,:,D%NKU) = ZFLXZ(:,:,IKE)
 END IF
 !  
 DO JK=IKTB+1,IKTE-1
-  PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+  PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+D%NKL))
 END DO
 !
-PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL)) 
+PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+D%NKL)) 
 !
 IF (OOCEAN) THEN
-  PWTH(:,:,IKE)=0.5*(ZFLXZ(:,:,IKE)+ZFLXZ(:,:,IKE+KKL))
+  PWTH(:,:,IKE)=0.5*(ZFLXZ(:,:,IKE)+ZFLXZ(:,:,IKE+D%NKL))
   PWTH(:,:,D%NKA)=0. 
-  PWTH(:,:,KKU)=ZFLXZ(:,:,KKU)
+  PWTH(:,:,D%NKU)=ZFLXZ(:,:,D%NKU)
 ELSE
-  PWTH(:,:,D%NKA)=0.5*(ZFLXZ(:,:,D%NKA)+ZFLXZ(:,:,D%NKA+KKL))
-  PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
+  PWTH(:,:,D%NKA)=0.5*(ZFLXZ(:,:,D%NKA)+ZFLXZ(:,:,D%NKA+D%NKL))
+  PWTH(:,:,IKE)=PWTH(:,:,IKE-D%NKL)
 END IF
 !
 IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -684,20 +680,20 @@ END IF
 !
 ! Contribution of the conservative temperature flux to the buoyancy flux
 IF (OOCEAN) THEN
-  PTP(:,:,:)= CST%XG*CST%XALPHAOC * MZF(ZFLXZ,D%NKA, KKU, KKL )
+  PTP(:,:,:)= CST%XG*CST%XALPHAOC * MZF(ZFLXZ,D%NKA, D%NKU, D%NKL )
 ELSE
   IF (KRR /= 0) THEN
-    PTP(:,:,:)  =  PBETA * MZF( MZM(PETHETA,D%NKA, KKU, KKL) * ZFLXZ,D%NKA, KKU, KKL )
+    PTP(:,:,:)  =  PBETA * MZF( MZM(PETHETA,D%NKA, D%NKU, D%NKL) * ZFLXZ,D%NKA, D%NKU, D%NKL )
     PTP(:,:,IKB)=  PBETA(:,:,IKB) * PETHETA(:,:,IKB) *   &
-                   0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+                   0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+D%NKL) )
   ELSE
-    PTP(:,:,:)=  PBETA * MZF( ZFLXZ,D%NKA, KKU, KKL )
+    PTP(:,:,:)=  PBETA * MZF( ZFLXZ,D%NKA, D%NKU, D%NKL )
   END IF
 END IF 
 !
 ! Buoyancy flux at flux points
 ! 
-PWTHV = MZM(PETHETA, D%NKA, KKU, KKL) * ZFLXZ
+PWTHV = MZM(PETHETA, D%NKA, D%NKU, D%NKL) * ZFLXZ
 PWTHV(:,:,IKB) = PETHETA(:,:,IKB) * ZFLXZ(:,:,IKB)
 !
 IF (OOCEAN) THEN
@@ -710,14 +706,14 @@ IF(HPROGRAM/='AROME  ') THEN
  IF ( KRRL >= 1 ) THEN
    IF ( KRRI >= 1 ) THEN
      PRRS(:,:,:,2) = PRRS(:,:,:,2) -                                        &
-                     PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, D%NKA, KKU, KKL) &
+                     PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, D%NKA, D%NKU, D%NKL) &
                      *(1.0-PFRAC_ICE(:,:,:))
      PRRS(:,:,:,4) = PRRS(:,:,:,4) -                                        &
-                     PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, D%NKA, KKU, KKL) &
+                     PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, D%NKA, D%NKU, D%NKL) &
                      *PFRAC_ICE(:,:,:)
    ELSE
      PRRS(:,:,:,2) = PRRS(:,:,:,2) -                                        &
-                     PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, D%NKA, KKU, KKL)
+                     PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, D%NKA, D%NKU, D%NKL)
    END IF
  END IF
 END IF
@@ -726,22 +722,22 @@ END IF
 ! 
 IF (OLES_CALL) THEN
   CALL SECOND_MNH(ZTIME1)
-  CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_WThl ) 
-  CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_W_SBG_WThl )
-  CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, D%NKA, KKU, KKL)*MZF(ZFLXZ, D%NKA, KKU, KKL),&
+  CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WThl ) 
+  CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_W_SBG_WThl )
+  CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, D%NKA, D%NKU, D%NKL)*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL),&
       & X_LES_RES_ddxa_W_SBG_UaThl )
-  CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_ddxa_Thl_SBG_UaThl )
-  CALL LES_MEAN_SUBGRID(-CSTURB%XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_ThlPz ) 
-  CALL LES_MEAN_SUBGRID(MZF(MZM(PETHETA, D%NKA, KKU, KKL)*ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_WThv ) 
+  CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_Thl_SBG_UaThl )
+  CALL LES_MEAN_SUBGRID(-CSTURB%XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_ThlPz ) 
+  CALL LES_MEAN_SUBGRID(MZF(MZM(PETHETA, D%NKA, D%NKU, D%NKL)*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WThv ) 
   IF (KRR>=1) THEN
-    CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_ddxa_Rt_SBG_UaThl )
+    CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_Rt_SBG_UaThl )
   END IF
   !* diagnostic of mixing coefficient for heat
-  ZA = DZM(PTHLP, D%NKA, KKU, KKL)
+  ZA = DZM(PTHLP, D%NKA, D%NKU, D%NKL)
   WHERE (ZA==0.) ZA=1.E-6
   ZA = - ZFLXZ / ZA * PDZZ
   ZA(:,:,IKB) = CSTURB%XCSHF*PPHI3(:,:,IKB)*ZKEFF(:,:,IKB)
-  ZA = MZF(ZA, D%NKA, KKU, KKL)
+  ZA = MZF(ZA, D%NKA, D%NKU, D%NKL)
   ZA = MIN(MAX(ZA,-1000.),1000.)
   CALL LES_MEAN_SUBGRID( ZA, X_LES_SUBGRID_Kh   ) 
   !
@@ -767,10 +763,10 @@ IF (KRR /= 0) THEN
   ! Compute the turbulent flux F and F' at time t-dt.
   !
  IF (OHARAT) THEN
-  ZF      (:,:,:) = -ZKEFF*DZM(PRM(:,:,:,1), D%NKA, KKU, KKL)/PDZZ
+  ZF      (:,:,:) = -ZKEFF*DZM(PRM(:,:,:,1), D%NKA, D%NKU, D%NKL)/PDZZ
   ZDFDDRDZ(:,:,:) = -ZKEFF
  ELSE
-  ZF      (:,:,:) = -CSTURB%XCSHF*PPSI3*ZKEFF*DZM(PRM(:,:,:,1), D%NKA, KKU, KKL)/PDZZ
+  ZF      (:,:,:) = -CSTURB%XCSHF*PPSI3*ZKEFF*DZM(PRM(:,:,:,1), D%NKA, D%NKU, D%NKL)/PDZZ
   ZDFDDRDZ(:,:,:) = -CSTURB%XCSHF*ZKEFF*D_PSI3DRDZ_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)
  ENDIF
   !
@@ -778,20 +774,20 @@ IF (KRR /= 0) THEN
   IF (TURBN%LHGRAD) THEN
     ZDELTAX= XXHAT(3) - XXHAT(2)
     ZF_LEONARD (:,:,:)= TURBN%XCOEFHGRADRM*ZDELTAX*ZDELTAX/12.0*(        &
-                MXF(GX_W_UW(PWM(:,:,:),  XDXX,XDZZ,XDZX,D%NKA,KKU,KKL))       &
-                *MZM(GX_M_M(PRM(:,:,:,1),XDXX,XDZZ,XDZX,D%NKA,KKU,KKL),D%NKA,KKU,KKL) &
-                +MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY,D%NKA,KKU,KKL))        &
-                *MZM(GY_M_M(PRM(:,:,:,1),XDYY,XDZZ,XDZY,D%NKA,KKU,KKL),D%NKA,KKU,KKL) )
+                MXF(GX_W_UW(PWM(:,:,:),  XDXX,XDZZ,XDZX,D%NKA,D%NKU,D%NKL))       &
+                *MZM(GX_M_M(PRM(:,:,:,1),XDXX,XDZZ,XDZX,D%NKA,D%NKU,D%NKL),D%NKA,D%NKU,D%NKL) &
+                +MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY,D%NKA,D%NKU,D%NKL))        &
+                *MZM(GY_M_M(PRM(:,:,:,1),XDYY,XDZZ,XDZY,D%NKA,D%NKU,D%NKL),D%NKA,D%NKU,D%NKL) )
    END IF
   !
   ! Effect of 3rd order terms in temperature flux (at flux point)
   !
   ! d(w'2r')/dz
   IF (GFWR) THEN
-    Z3RDMOMENT= M3_WR_W2R(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,ZKEFF,PTKEM)
+    Z3RDMOMENT= M3_WR_W2R(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,ZKEFF,PTKEM)
   !
     ZF       = ZF       + Z3RDMOMENT * PFWR
-    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2R_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2R_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,&
      & PBLL_O_E,PEMOIST,ZKEFF,PTKEM) * PFWR
   END IF
   !
@@ -799,35 +795,35 @@ IF (KRR /= 0) THEN
   IF (GFR2) THEN
     Z3RDMOMENT= M3_WR_WR2(D,CSTURB,PREDR1,PREDTH1,PD,PBLL_O_E,PEMOIST)
   !
-    ZF       = ZF       + Z3RDMOMENT * MZM(PFR2, D%NKA, KKU, KKL)
+    ZF       = ZF       + Z3RDMOMENT * MZM(PFR2, D%NKA, D%NKU, D%NKL)
     ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WR2_O_DDRDZ(D,CSTURB,Z3RDMOMENT,PREDR1,&
-     & PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFR2, D%NKA, KKU, KKL)
+     & PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFR2, D%NKA, D%NKU, D%NKL)
   END IF
   !
   ! d(w'2th')/dz
   IF (GFWTH) THEN
-    ZF       = ZF       + M3_WR_W2TH(D,CSTURB,D%NKA,KKU,KKL,PD,ZKEFF,&
+    ZF       = ZF       + M3_WR_W2TH(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,ZKEFF,&
      & PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * PFWTH
-    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2TH_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,& 
+    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2TH_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,& 
      & PD,ZKEFF,PTKEM,PBLL_O_E,PETHETA) * PFWTH
   END IF
   !
   ! d(w'th'2)/dz
   IF (GFTH2) THEN
-    ZF       = ZF       + M3_WR_WTH2(D,CSTURB,D%NKA,KKU,KKL,PD,ZKEFF,PTKEM,&
-    & PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDR_DZ) * MZM(PFTH2, D%NKA, KKU, KKL)
-    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTH2_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,PREDR1,PREDTH1,PD,&
-     &ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA) * MZM(PFTH2, D%NKA, KKU, KKL)
+    ZF       = ZF       + M3_WR_WTH2(D,CSTURB,D%NKA,D%NKU,D%NKL,PD,ZKEFF,PTKEM,&
+    & PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDR_DZ) * MZM(PFTH2, D%NKA, D%NKU, D%NKL)
+    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTH2_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDR1,PREDTH1,PD,&
+     &ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA) * MZM(PFTH2, D%NKA, D%NKU, D%NKL)
   END IF
   !
   ! d(w'th'r')/dz
   IF (GFTHR) THEN
-    Z3RDMOMENT= M3_WR_WTHR(D,CSTURB,D%NKA,KKU,KKL,PREDTH1,PD,ZKEFF,PTKEM,PSQRT_TKE,PBETA,&
+    Z3RDMOMENT= M3_WR_WTHR(D,CSTURB,D%NKA,D%NKU,D%NKL,PREDTH1,PD,ZKEFF,PTKEM,PSQRT_TKE,PBETA,&
      & PLEPS,PETHETA)
   !
-    ZF       = ZF       + Z3RDMOMENT * MZM(PFTHR, D%NKA, KKU, KKL)
-    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTHR_O_DDRDZ(D,CSTURB,D%NKA,KKU,KKL,Z3RDMOMENT,PREDR1, &
-     & PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR, D%NKA, KKU, KKL)
+    ZF       = ZF       + Z3RDMOMENT * MZM(PFTHR, D%NKA, D%NKU, D%NKL)
+    ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTHR_O_DDRDZ(D,CSTURB,D%NKA,D%NKU,D%NKL,Z3RDMOMENT,PREDR1, &
+     & PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR, D%NKA, D%NKU, D%NKL)
   END IF
   !
   ! compute interface flux
@@ -870,7 +866,7 @@ IF (KRR /= 0) THEN
     END IF
   END IF!end no coupled cases
   ! Compute the split conservative potential temperature at t+deltat
-  CALL TRIDIAG_THERMO(D,D%NKA,KKU,KKL,PRM(:,:,:,1),ZF,ZDFDDRDZ,PTSTEP,PIMPL,&
+  CALL TRIDIAG_THERMO(D,D%NKA,D%NKU,D%NKL,PRM(:,:,:,1),ZF,ZDFDDRDZ,PTSTEP,PIMPL,&
                       PDZZ,PRHODJ,PRP)
   !
   ! Compute the equivalent tendency for the conservative mixing ratio
@@ -879,11 +875,11 @@ IF (KRR /= 0) THEN
   !
   ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
   IF (TURBN%LHGRAD) THEN
-   DO JK=1,KKU
+   DO JK=1,D%NKU
     ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
    END DO
    WHERE ( (ZCLD_THOLD(:,:,:) >= TURBN%XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= TURBN%XALTHGRAD ) )
-    ZRWRNP (:,:,:) =  -GZ_W_M(MZM(PRHODJ(:,:,:),D%NKA,KKU,KKL)*ZF_LEONARD(:,:,:),XDZZ,D%NKA,KKU,KKL)
+    ZRWRNP (:,:,:) =  -GZ_W_M(MZM(PRHODJ(:,:,:),D%NKA,D%NKU,D%NKL)*ZF_LEONARD(:,:,:),XDZZ,D%NKA,D%NKU,D%NKL)
    END WHERE
   END IF
   !
@@ -894,7 +890,7 @@ IF (KRR /= 0) THEN
   ! cons. mixing ratio flux :
   !
   ZFLXZ(:,:,:)   = ZF                                                &
-                 + PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), D%NKA, KKU, KKL) / PDZZ 
+                 + PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), D%NKA, D%NKU, D%NKL) / PDZZ 
   !
   ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
   IF (TURBN%LHGRAD) THEN
@@ -906,11 +902,11 @@ IF (KRR /= 0) THEN
   ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB) 
   !
   DO JK=IKTB+1,IKTE-1
-   PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+   PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+D%NKL))
   END DO
-  PWRC(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
-  PWRC(:,:,D%NKA)=0.5*(ZFLXZ(:,:,D%NKA)+ZFLXZ(:,:,D%NKA+KKL))
-  PWRC(:,:,IKE)=PWRC(:,:,IKE-KKL)
+  PWRC(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+D%NKL))
+  PWRC(:,:,D%NKA)=0.5*(ZFLXZ(:,:,D%NKA)+ZFLXZ(:,:,D%NKA+D%NKL))
+  PWRC(:,:,IKE)=PWRC(:,:,IKE-D%NKL)
   !
   !
   IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -930,17 +926,17 @@ IF (KRR /= 0) THEN
   !
   ! Contribution of the conservative water flux to the Buoyancy flux
   IF (OOCEAN) THEN
-     ZA(:,:,:)=  -CST%XG*CST%XBETAOC  * MZF(ZFLXZ, D%NKA, KKU, KKL )
+     ZA(:,:,:)=  -CST%XG*CST%XBETAOC  * MZF(ZFLXZ, D%NKA, D%NKU, D%NKL )
   ELSE
-    ZA(:,:,:)   =  PBETA * MZF( MZM(PEMOIST, D%NKA, KKU, KKL) * ZFLXZ,D%NKA,KKU,KKL )
+    ZA(:,:,:)   =  PBETA * MZF( MZM(PEMOIST, D%NKA, D%NKU, D%NKL) * ZFLXZ,D%NKA,D%NKU,D%NKL )
     ZA(:,:,IKB) =  PBETA(:,:,IKB) * PEMOIST(:,:,IKB) *   &
-                   0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+                   0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+D%NKL) )
     PTP(:,:,:) = PTP(:,:,:) + ZA(:,:,:)
   END IF
   !
   ! Buoyancy flux at flux points
   ! 
-  PWTHV          = PWTHV          + MZM(PEMOIST, D%NKA, KKU, KKL) * ZFLXZ
+  PWTHV          = PWTHV          + MZM(PEMOIST, D%NKA, D%NKU, D%NKL) * ZFLXZ
   PWTHV(:,:,IKB) = PWTHV(:,:,IKB) + PEMOIST(:,:,IKB) * ZFLXZ(:,:,IKB)
   IF (OOCEAN) THEN
     PWTHV(:,:,IKE) = PWTHV(:,:,IKE) + PEMOIST(:,:,IKE)* ZFLXZ(:,:,IKE)
@@ -952,14 +948,14 @@ IF(HPROGRAM/='AROME  ') THEN
    IF ( KRRL >= 1 ) THEN
      IF ( KRRI >= 1 ) THEN
        PRRS(:,:,:,2) = PRRS(:,:,:,2) -                                        &
-                       PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,D%NKA,KKU,KKL )       &
+                       PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,D%NKA,D%NKU,D%NKL )       &
                        *(1.0-PFRAC_ICE(:,:,:))
        PRRS(:,:,:,4) = PRRS(:,:,:,4) -                                        &
-                       PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,D%NKA,KKU,KKL )       &
+                       PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,D%NKA,D%NKU,D%NKL )       &
                        *PFRAC_ICE(:,:,:)
      ELSE
        PRRS(:,:,:,2) = PRRS(:,:,:,2) -                                        &
-                       PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,D%NKA,KKU,KKL )
+                       PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,D%NKA,D%NKU,D%NKL )
      END IF
    END IF
 END IF
@@ -968,14 +964,14 @@ END IF
 ! 
   IF (OLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
-    CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_WRt ) 
-    CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_W_SBG_WRt )
-    CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, D%NKA, KKU, KKL)*MZF(ZFLXZ, D%NKA, KKU, KKL),&
+    CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WRt ) 
+    CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_W_SBG_WRt )
+    CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, D%NKA, D%NKU, D%NKL)*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL),&
     & X_LES_RES_ddxa_W_SBG_UaRt )
-    CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_ddxa_Thl_SBG_UaRt )
-    CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, D%NKA, KKU, KKL), X_LES_RES_ddxa_Rt_SBG_UaRt )
-    CALL LES_MEAN_SUBGRID(MZF(MZM(PEMOIST, D%NKA, KKU, KKL)*ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_WThv , .TRUE. ) 
-    CALL LES_MEAN_SUBGRID(-CSTURB%XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_RtPz ) 
+    CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_Thl_SBG_UaRt )
+    CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_RES_ddxa_Rt_SBG_UaRt )
+    CALL LES_MEAN_SUBGRID(MZF(MZM(PEMOIST, D%NKA, D%NKU, D%NKL)*ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WThv , .TRUE. ) 
+    CALL LES_MEAN_SUBGRID(-CSTURB%XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_RtPz ) 
     CALL SECOND_MNH(ZTIME2)
     XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
   END IF
@@ -997,18 +993,18 @@ IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. OLES_CALL) .AND. (KRRL > 0) ) THEN
 ! With OHARAT is true tke and length scales at half levels
 ! yet modify to use length scale and tke at half levels from vdfexcuhl
  IF (OHARAT) THEN
-  ZA(:,:,:)   = DZM(PIMPL * PTHLP + PEXPL * PTHLM, D%NKA, KKU, KKL) / PDZZ *       &
+  ZA(:,:,:)   = DZM(PIMPL * PTHLP + PEXPL * PTHLM, D%NKA, D%NKU, D%NKL) / PDZZ *       &
                   (-PLM*PSQRT_TKE) 
  ELSE
-  ZA(:,:,:)   = DZM(PIMPL * PTHLP + PEXPL * PTHLM, D%NKA, KKU, KKL) / PDZZ *       &
-                  (-PPHI3*MZM(PLM*PSQRT_TKE, D%NKA, KKU, KKL)) * CSTURB%XCSHF 
+  ZA(:,:,:)   = DZM(PIMPL * PTHLP + PEXPL * PTHLM, D%NKA, D%NKU, D%NKL) / PDZZ *       &
+                  (-PPHI3*MZM(PLM*PSQRT_TKE, D%NKA, D%NKU, D%NKL)) * CSTURB%XCSHF 
  ENDIF
   ZA(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
                * PDIRCOSZW(:,:)
   !  
   ! compute <w Rc>
-  ZFLXZ(:,:,:) = MZM(PAMOIST * 2.* PSRCM, D%NKA, KKU, KKL) * ZFLXZ(:,:,:) + &
-                 MZM(PATHETA * 2.* PSRCM, D%NKA, KKU, KKL) * ZA(:,:,:)
+  ZFLXZ(:,:,:) = MZM(PAMOIST * 2.* PSRCM, D%NKA, D%NKU, D%NKL) * ZFLXZ(:,:,:) + &
+                 MZM(PATHETA * 2.* PSRCM, D%NKA, D%NKU, D%NKL) * ZA(:,:,:)
   ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB) 
   !                 
   ! store the liquid water mixing ratio vertical flux
@@ -1030,7 +1026,7 @@ IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. OLES_CALL) .AND. (KRRL > 0) ) THEN
 !
   IF (OLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
-    CALL LES_MEAN_SUBGRID( MZF(ZFLXZ, D%NKA, KKU, KKL), X_LES_SUBGRID_WRc ) 
+    CALL LES_MEAN_SUBGRID( MZF(ZFLXZ, D%NKA, D%NKU, D%NKL), X_LES_SUBGRID_WRc ) 
     CALL SECOND_MNH(ZTIME2)
     XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
   END IF
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 46dbec082..01d939b6e 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -500,15 +500,12 @@ IKE=D%NKE
 ZEXPL = 1.- PIMPL
 ZRVORD= CST%XRV / CST%XRD
 !
-!
 !Copy data into ZTHLM and ZRM only if needed
 IF (HTURBLEN=='BL89' .OR. HTURBLEN=='RM17' .OR. ORMC01) THEN
   ZTHLM(:,:,:) = PTHLT(:,:,:)
   ZRM(:,:,:,:) = PRT(:,:,:,:)
 END IF
 !
-!
-!
 !----------------------------------------------------------------------------
 !
 !*      2. COMPUTE CONSERVATIVE VARIABLES AND RELATED QUANTITIES
@@ -649,7 +646,7 @@ SELECT CASE (HTURBLEN)
 
   CASE ('BL89')
     ZSHEAR=0.
-    CALL BL89(D,CST,CSTURB,D%NKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN,HPROGRAM)
+    CALL BL89(D,CST,CSTURB,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN,HPROGRAM)
 !
 !*      3.2 RM17 mixing length
 !           ------------------
@@ -658,7 +655,7 @@ SELECT CASE (HTURBLEN)
     ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ,D%NKA,KKU,KKL),D%NKA,KKU,KKL))
     ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ,D%NKA,KKU,KKL),D%NKA,KKU,KKL))
     ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
-    CALL BL89(D,CST,CSTURB,D%NKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN,HPROGRAM)
+    CALL BL89(D,CST,CSTURB,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN,HPROGRAM)
 !
 !*      3.3 Grey-zone combined RM17 & Deardorff mixing lengths 
 !           --------------------------------------------------
@@ -667,7 +664,7 @@ SELECT CASE (HTURBLEN)
     ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ,D%NKA,KKU,KKL),D%NKA,KKU,KKL))
     ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ,D%NKA,KKU,KKL),D%NKA,KKU,KKL))
     ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
-    CALL BL89(D,CST,CSTURB,D%NKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN,HPROGRAM)
+    CALL BL89(D,CST,CSTURB,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN,HPROGRAM)
 
     CALL DELT(ZLMW,ODZ=.FALSE.)
     ! The minimum mixing length is chosen between Horizontal grid mesh (not taking into account the vertical grid mesh) and RM17.
@@ -881,7 +878,7 @@ IF( BUCONF%LBUDGET_SV ) THEN
   END DO
 END IF
 
-CALL TURB_VER(D, CST,CSTURB,TURBN,D%NKA,KKU,KKL,KRR, KRRL, KRRI,&
+CALL TURB_VER(D, CST,CSTURB,TURBN,KRR, KRRL, KRRI,       &
           OTURB_FLX, OOCEAN, ODEEPOC, OHARAT,OCOMPUTE_SRC,&
           KSV,KSV_LGBEG,KSV_LGEND,                       &
           HTURBDIM,HTOM,PIMPL,ZEXPL,                     &
@@ -1064,7 +1061,7 @@ ELSE
 END IF
 !
 CALL TKE_EPS_SOURCES(D,CST,CSTURB,BUCONF,HPROGRAM,&
-                   &D%NKA,KKU,KKL,KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH,       &
+                   & KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH,       &
                    & PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,            &
                    & PTSTEP,PIMPL,ZEXPL,                         &
                    & HTURBLEN,HTURBDIM,                              &
@@ -1760,7 +1757,7 @@ ELSE
 !           ------------------
   CASE ('BL89','RM17','ADAP')
     ZSHEAR=0.
-    CALL BL89(D,CST,CSTURB,D%NKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD,OOCEAN,HPROGRAM)
+    CALL BL89(D,CST,CSTURB,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD,OOCEAN,HPROGRAM)
 !
 !*         3.2 Delta mixing length
 !           -------------------
-- 
GitLab