diff --git a/src/mesonh/ext/phys_paramn.f90 b/src/mesonh/ext/phys_paramn.f90
index c335a4652eaa503c75ca1612af8df6a555b385ce..6a3a0f32ef89d83d894532c82a2be176470c0ab7 100644
--- a/src/mesonh/ext/phys_paramn.f90
+++ b/src/mesonh/ext/phys_paramn.f90
@@ -1496,8 +1496,8 @@ END IF
               XCEI, XCEI_MIN, XCEI_MAX, XCOEF_AMPL_SAT,                              &
               XTHT, XRT,                                                             &
               XRUS, XRVS, XRWS, XRTHS, XRRS, XRSVS, XRTKES, XRTKEMS, XSIGS, XWTHVMF, &
-              XTHW_FLUX, XRCW_FLUX, XSVW_FLUX,XDYP, XTHP, XTR, XDISS, XLEM,          &
-              TBUDGETS, KBUDGETS=SIZE(TBUDGETS)      )
+              XTHW_FLUX, XRCW_FLUX, XSVW_FLUX,XDYP, XTHP, XTR, XDISS,          &
+              TBUDGETS, KBUDGETS=SIZE(TBUDGETS),PLEM=XLEM      )
 !
 IF (LRMC01) THEN
   CALL ADD2DFIELD_ll( TZFIELDS_ll, XSBL_DEPTH, 'PHYS_PARAM_n::XSBL_DEPTH' )
diff --git a/src/mesonh/turb/bl_depth_diag.f90 b/src/mesonh/turb/bl_depth_diag.f90
deleted file mode 100644
index 2e7fb121cdda00386511bb0254c18249bec192bb..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/bl_depth_diag.f90
+++ /dev/null
@@ -1,200 +0,0 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 turb 2006/05/18 13:07:25
-!-----------------------------------------------------------------
-!    ################ 
-     MODULE MODI_BL_DEPTH_DIAG  
-!    ################ 
-!
-INTERFACE BL_DEPTH_DIAG  
-!
-!
-      FUNCTION BL_DEPTH_DIAG_3D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
-
-INTEGER,                INTENT(IN)           :: KKB          ! bottom point
-INTEGER,                INTENT(IN)           :: KKE          ! top point
-REAL, DIMENSION(:,:),   INTENT(IN)           :: PSURF        ! surface flux
-REAL, DIMENSION(:,:),   INTENT(IN)           :: PZS          ! orography
-REAL, DIMENSION(:,:,:), INTENT(IN)           :: PFLUX        ! flux
-REAL, DIMENSION(:,:,:), INTENT(IN)           :: PZZ          ! altitude of flux points
-REAL,                   INTENT(IN)           :: PFTOP_O_FSURF! Flux at BL top / Surface flux
-REAL, DIMENSION(SIZE(PSURF,1),SIZE(PSURF,2)) :: BL_DEPTH_DIAG_3D
-!
-END FUNCTION BL_DEPTH_DIAG_3D
-!
-!
-      FUNCTION BL_DEPTH_DIAG_1D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
-INTEGER,                INTENT(IN)           :: KKB          ! bottom point
-INTEGER,                INTENT(IN)           :: KKE          ! top point
-REAL,                   INTENT(IN)           :: PSURF        ! surface flux
-REAL,                   INTENT(IN)           :: PZS          ! orography
-REAL, DIMENSION(:),     INTENT(IN)           :: PFLUX        ! flux
-REAL, DIMENSION(:),     INTENT(IN)           :: PZZ          ! altitude of flux points
-REAL,                   INTENT(IN)           :: PFTOP_O_FSURF! Flux at BL top / Surface flux
-REAL                                         :: BL_DEPTH_DIAG_1D
-!
-END FUNCTION BL_DEPTH_DIAG_1D
-!
-END INTERFACE
-!
-END MODULE MODI_BL_DEPTH_DIAG
-!
-!-------------------------------------------------------------------------------
-!
-!    ################ 
-     MODULE MODI_BL_DEPTH_DIAG_3D  
-!    ################ 
-!
-!
-INTERFACE
-!
-!
-      FUNCTION BL_DEPTH_DIAG_3D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
-INTEGER,                INTENT(IN)           :: KKB          ! bottom point
-INTEGER,                INTENT(IN)           :: KKE          ! top point
-REAL, DIMENSION(:,:),   INTENT(IN)           :: PSURF        ! surface flux
-REAL, DIMENSION(:,:),   INTENT(IN)           :: PZS          ! orography
-REAL, DIMENSION(:,:,:), INTENT(IN)           :: PFLUX        ! flux
-REAL, DIMENSION(:,:,:), INTENT(IN)           :: PZZ          ! altitude of flux points
-REAL,                   INTENT(IN)           :: PFTOP_O_FSURF! Flux at BL top / Surface flux
-REAL, DIMENSION(SIZE(PSURF,1),SIZE(PSURF,2)) :: BL_DEPTH_DIAG_3D
-!
-END FUNCTION BL_DEPTH_DIAG_3D
-!
-!
-END INTERFACE
-!
-END MODULE MODI_BL_DEPTH_DIAG_3D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BL_DEPTH_DIAG_3D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
-!
-!
-!!****  *SBL_DEPTH* - computes SBL depth
-!!
-!!    PURPOSE
-!!    -------
-!
-!!**  METHOD
-!!    ------
-!!
-!!    SBL is defined as the layer where momentum flux is equal to XSBL_FRAC of its surface value
-!!    
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!      V. Masson * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original         nov. 2005
-!!
-!! --------------------------------------------------------------------------
-!       
-!*      0. DECLARATIONS
-!          ------------
-!
-!*      0.1  declarations of arguments
-!
-IMPLICIT NONE
-!
-INTEGER,                INTENT(IN)           :: KKB          ! bottom point
-INTEGER,                INTENT(IN)           :: KKE          ! top point
-REAL, DIMENSION(:,:),   INTENT(IN)           :: PSURF        ! surface flux
-REAL, DIMENSION(:,:),   INTENT(IN)           :: PZS          ! orography
-REAL, DIMENSION(:,:,:), INTENT(IN)           :: PFLUX        ! flux
-REAL, DIMENSION(:,:,:), INTENT(IN)           :: PZZ          ! altitude of flux points
-REAL,                   INTENT(IN)           :: PFTOP_O_FSURF! Flux at BL top / Surface flux
-REAL, DIMENSION(SIZE(PSURF,1),SIZE(PSURF,2)) :: BL_DEPTH_DIAG_3D
-!
-!
-!       0.2  declaration of local variables
-!
-INTEGER :: JI,JJ,JK ! loop counters
-INTEGER :: IKL      ! +1 : MesoNH levels -1: Arome
-REAL    :: ZFLX     ! flux at top of BL
-!
-!----------------------------------------------------------------------------
-!
-IF (KKB < KKE) THEN
-  IKL=1
-ELSE
-  IKL=-1
-ENDIF
-
-BL_DEPTH_DIAG_3D(:,:) = 0.
-!
-
-DO JJ=1,SIZE(PSURF,2)
-  DO JI=1,SIZE(PSURF,1)
-    IF (PSURF(JI,JJ)==0.) CYCLE
-    DO JK=KKB,KKE,IKL
-      IF (PZZ(JI,JJ,JK-IKL)<=PZS(JI,JJ)) CYCLE
-      ZFLX = PSURF(JI,JJ) * PFTOP_O_FSURF
-      IF ( (PFLUX(JI,JJ,JK)-ZFLX)*(PFLUX(JI,JJ,JK-IKL)-ZFLX) <= 0. ) THEN
-        BL_DEPTH_DIAG_3D(JI,JJ) = (PZZ  (JI,JJ,JK-IKL) - PZS(JI,JJ))     &
-                         + (PZZ  (JI,JJ,JK) - PZZ  (JI,JJ,JK-IKL))    &
-                         * (ZFLX            - PFLUX(JI,JJ,JK-IKL)  )  &
-                         / (PFLUX(JI,JJ,JK) - PFLUX(JI,JJ,JK-IKL)   )
-        EXIT
-      END IF
-    END DO
-  END DO
-END DO
-!
-BL_DEPTH_DIAG_3D(:,:) = BL_DEPTH_DIAG_3D(:,:) / (1. - PFTOP_O_FSURF)
-!
-END FUNCTION BL_DEPTH_DIAG_3D
-!
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-FUNCTION BL_DEPTH_DIAG_1D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
-!
-USE MODI_BL_DEPTH_DIAG_3D
-IMPLICIT NONE
-!
-INTEGER,                INTENT(IN)           :: KKB          ! bottom point
-INTEGER,                INTENT(IN)           :: KKE          ! top point
-REAL,                   INTENT(IN)           :: PSURF        ! surface flux
-REAL,                   INTENT(IN)           :: PZS          ! orography
-REAL, DIMENSION(:),     INTENT(IN)           :: PFLUX        ! flux
-REAL, DIMENSION(:),     INTENT(IN)           :: PZZ          ! altitude of flux points
-REAL,                   INTENT(IN)           :: PFTOP_O_FSURF! Flux at BL top / Surface flux
-REAL                                         :: BL_DEPTH_DIAG_1D
-!
-REAL, DIMENSION(1,1)             :: ZSURF
-REAL, DIMENSION(1,1)             :: ZZS
-REAL, DIMENSION(1,1,SIZE(PFLUX)) :: ZFLUX
-REAL, DIMENSION(1,1,SIZE(PZZ))   :: ZZZ
-REAL, DIMENSION(1,1)             :: ZBL_DEPTH_DIAG
-!
-ZSURF        = PSURF
-ZZS          = PZS
-ZFLUX(1,1,:) = PFLUX(:)
-ZZZ  (1,1,:) = PZZ  (:)
-!
-ZBL_DEPTH_DIAG = BL_DEPTH_DIAG_3D(KKB,KKE,ZSURF,ZZS,ZFLUX,ZZZ,PFTOP_O_FSURF)
-!
-BL_DEPTH_DIAG_1D = ZBL_DEPTH_DIAG(1,1)
-!
-!-------------------------------------------------------------------------------
-!
-END FUNCTION BL_DEPTH_DIAG_1D
diff --git a/src/mesonh/turb/mode_turb_ver.f90 b/src/mesonh/turb/mode_turb_ver.f90
index abec958b5d8d80efb56d6bca107d0fbf9f16ec8f..bdd4f3299de591714a22075dc82677c8be97b684 100644
--- a/src/mesonh/turb/mode_turb_ver.f90
+++ b/src/mesonh/turb/mode_turb_ver.f90
@@ -219,18 +219,18 @@ USE MODD_LES
 USE MODD_NSV,            ONLY: NSV
 !
 !USE MODE_PRANDTL, ONLY: PRANDTL
-USE MODI_EMOIST
-USE MODI_ETHETA
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
 USE MODI_GRADIENT_M
 USE MODI_GRADIENT_W
-USE MODI_TURB
+!USE MODI_TURB
 USE MODE_TURB_VER_THERMO_FLUX, ONLY: TURB_VER_THERMO_FLUX
 USE MODE_TURB_VER_THERMO_CORR, ONLY: TURB_VER_THERMO_CORR
 USE MODE_TURB_VER_DYN_FLUX, ONLY: TURB_VER_DYN_FLUX
 USE MODE_TURB_VER_SV_FLUX, ONLY: TURB_VER_SV_FLUX
 USE MODE_TURB_VER_SV_CORR, ONLY: TURB_VER_SV_CORR
 USE MODI_LES_MEAN_SUBGRID
-USE MODI_SBL_DEPTH
+USE MODE_SBL_DEPTH, ONLY: SBL_DEPTH
 USE MODI_SECOND_MNH
 !
 USE MODE_IO_FIELD_WRITE, only: IO_Field_write
diff --git a/src/mesonh/turb/mode_turb_ver_sv_corr.f90 b/src/mesonh/turb/mode_turb_ver_sv_corr.f90
index 0544e22ace27dfbccdb862de87be1179e910b30a..4deb47ce7f0a66cff11aacf03717b547bbd5103f 100644
--- a/src/mesonh/turb/mode_turb_ver_sv_corr.f90
+++ b/src/mesonh/turb/mode_turb_ver_sv_corr.f90
@@ -69,8 +69,8 @@ USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : MZF
-USE MODI_EMOIST
-USE MODI_ETHETA
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
 USE MODI_LES_MEAN_SUBGRID
 !
 USE MODI_SECOND_MNH
diff --git a/src/mesonh/turb/mode_turb_ver_sv_flux.f90 b/src/mesonh/turb/mode_turb_ver_sv_flux.f90
index 624f9b12baa16700eac53587cd1ef4f829f0164f..271ee734e6aabd3d2cda09b7eba343113fa6e402 100644
--- a/src/mesonh/turb/mode_turb_ver_sv_flux.f90
+++ b/src/mesonh/turb/mode_turb_ver_sv_flux.f90
@@ -227,8 +227,8 @@ USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
 USE MODE_TRIDIAG, ONLY: TRIDIAG
-USE MODI_EMOIST
-USE MODI_ETHETA
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
 USE MODI_LES_MEAN_SUBGRID
 !
 USE MODI_SECOND_MNH
diff --git a/src/mesonh/turb/modi_turb.f90 b/src/mesonh/turb/modi_turb.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a87b97854a7cf790abcf810fb5f56d9c3bc2b354
--- /dev/null
+++ b/src/mesonh/turb/modi_turb.f90
@@ -0,0 +1,143 @@
+!     ######spl
+     MODULE MODI_TURB  
+!    ################ 
+!
+INTERFACE
+!
+!
+      SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,      &
+              & KSPLIT,KMODEL_CL,                                     &
+              & OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,               &
+              & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HCLOUD,PIMPL,      &
+              & PTSTEP,TPFILE,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,           &
+              & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE,    &
+              & PRHODJ,PTHVREF,                                       &
+              & PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
+              & PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT,                  &
+              & PLENGTHM,PLENGTHH,MFMOIST,                            &
+              & PBL_DEPTH, PSBL_DEPTH,                                &
+              & PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,                &
+              & PTHLT,PRT,                                            &
+              & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES,PRTKEMS,PSIGS,&
+              & PFLXZTHVMF,PWTH,PWRC,PWSV,PDYP,PTHP,PTDIFF,PTDISS,    &
+              & TBUDGETS, KBUDGETS                                   ,&
+              & PTR,PDISS,PEDR,PLEM                            )
+
+USE MODD_BUDGET, ONLY : TBUDGETDATA
+USE MODD_IO, ONLY : TFILEDATA
+!
+!
+!
+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  
+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.
+CHARACTER(LEN=*),DIMENSION(:),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
+INTEGER,                INTENT(IN)   :: KSPLIT        ! number of time-splitting
+INTEGER,                INTENT(IN)   :: KMODEL_CL     ! model number for cloud mixing length
+LOGICAL,                INTENT(IN)   ::  OTURB_FLX    ! switch to write the
+                                 ! turbulent fluxes in the syncronous FM-file
+LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
+                                 ! diagnostic fields in the syncronous FM-file
+LOGICAL,                INTENT(IN)   ::  OSUBG_COND   ! switch for SUBGrid 
+                                 ! CONDensation
+LOGICAL,                INTENT(IN)   ::  ORMC01       ! switch for RMC01 lengths in SBL
+CHARACTER(LEN=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
+                                                      ! turbulence scheme
+CHARACTER(LEN=4),       INTENT(IN)   ::  HTURBLEN     ! kind of mixing length
+CHARACTER(LEN=4),       INTENT(IN)   ::  HTOM         ! kind of Third Order Moment
+CHARACTER(LEN=4),       INTENT(IN)   ::  HTURBLEN_CL  ! kind of cloud mixing length
+REAL,                   INTENT(IN)   ::  PIMPL        ! degree of implicitness
+CHARACTER (LEN=4),      INTENT(IN)   ::  HCLOUD       ! Kind of microphysical scheme
+REAL,                   INTENT(IN)   ::  PTSTEP       ! timestep 
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDXX,PDYY,PDZZ,PDZX,PDZY
+                                        ! metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN)   :: PZZ       !  physical distance 
+! between 2 succesive grid points along the K direction
+REAL, DIMENSION(:,:),   INTENT(IN)      ::  PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
+! Director Cosinus along x, y and z directions at surface w-point
+REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCOSSLOPE       ! cosinus of the angle
+                                 ! between i and the slope vector
+REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE       ! sinus of the angle
+                                 ! between i and the slope vector
+REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODJ    ! dry density * Grid size
+REAL, DIMENSION(:,:,:), INTENT(IN)      ::  MFMOIST ! moist mass flux dual scheme
+REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PTHVREF   ! Virtual Potential
+                                        ! Temperature of the reference state
+!
+REAL, DIMENSION(:,:),   INTENT(IN)      ::  PSFTH,PSFRV,   &
+! normal surface fluxes of theta and Rv 
+                                            PSFU,PSFV
+! normal surface fluxes of (u,v) parallel to the orography
+REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PSFSV
+! normal surface fluxes of Scalar var. 
+!
+!    prognostic variables at t- deltat
+REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABST      ! Pressure at time t
+REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUT,PVT,PWT ! wind components
+REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKET       ! TKE
+REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVT        ! passive scal. var.
+REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PSRCT       ! Second-order flux
+                      ! s'rc'/2Sigma_s2 at time t-1 multiplied by Lambda_3
+REAL, DIMENSION(:,:),     INTENT(INOUT) :: PBL_DEPTH  ! BL height for TOMS
+REAL, DIMENSION(:,:),     INTENT(INOUT) :: PSBL_DEPTH ! SBL depth for RMC01
+!
+!    variables for cloud mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PCEI ! Cloud Entrainment instability
+                                                 ! index to emphasize localy 
+                                                 ! turbulent fluxes
+REAL, INTENT(IN)      ::  PCEI_MIN ! minimum threshold for the instability index CEI
+REAL, INTENT(IN)      ::  PCEI_MAX ! maximum threshold for the instability index CEI
+REAL, INTENT(IN)      ::  PCOEF_AMPL_SAT ! saturation of the amplification coefficient
+!
+!   thermodynamical variables which are transformed in conservative var.
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PTHLT       ! conservative pot. temp.
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRT         ! water var.  where 
+                             ! PRT(:,:,:,1) is the conservative mixing ratio        
+!
+! sources of momentum, conservative potential temperature, Turb. Kin. Energy, 
+! TKE dissipation
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PRUS,PRVS,PRWS,PRTHLS,PRTKES
+! Source terms for all water kinds, PRRS(:,:,:,1) is used for the conservative
+! mixing ratio
+REAL, DIMENSION(:,:,:),   INTENT(IN)    ::  PRTKEMS
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRRS 
+! Source terms for all passive scalar variables
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRSVS
+! Sigma_s at time t+1 : square root of the variance of the deviation to the 
+! saturation
+REAL, DIMENSION(:,:,:), INTENT(OUT)     ::  PSIGS
+REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PFLXZTHVMF 
+!                                           MF contribution for vert. turb. transport
+!                                           used in the buoy. prod. of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWTH       ! heat flux
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWRC       ! cloud water flux
+REAL, DIMENSION(:,:,:,:),INTENT(OUT) :: PWSV       ! scalar flux
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDYP  ! Dynamical production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHP  ! Thermal production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTDIFF   ! Transport production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTDISS ! Dissipation of TKE
+!
+TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
+INTEGER, INTENT(IN) :: KBUDGETS
+!
+! length scale from vdfexcu
+REAL, DIMENSION(:,:,:), INTENT(IN)    :: PLENGTHM, PLENGTHH
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PTR   ! Transport production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PDISS ! Dissipation of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PEDR  ! EDR
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PLEM  ! Mixing length
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE TURB
+!
+END INTERFACE
+!
+END MODULE MODI_TURB
diff --git a/src/mesonh/turb/th_r_from_thl_rt_1d.f90 b/src/mesonh/turb/th_r_from_thl_rt_1d.f90
deleted file mode 100644
index 8c58fe8d0f9b784c89c6f3f0455e2fe425b6b11d..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/th_r_from_thl_rt_1d.f90
+++ /dev/null
@@ -1,203 +0,0 @@
-!MNH_LIC Copyright 2006-2021 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!     ######spl
-      MODULE MODI_TH_R_FROM_THL_RT_1D
-!     ###############################
-!
-      INTERFACE
-      SUBROUTINE TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI,         &
-                                  PRSATW, PRSATI, PRR, PRS, PRG, PRH     )
-CHARACTER(len=1),   INTENT(IN)    :: HFRAC_ICE
-REAL, DIMENSION(:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:), INTENT(IN) :: PP     ! Pressure
-REAL, DIMENSION(:), INTENT(IN) :: PTHL   ! Liquid pot. temp.
-REAL, DIMENSION(:), INTENT(IN) :: PRT    ! Total mixing ratios 
-REAL, DIMENSION(:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:), INTENT(OUT):: PTH    ! Potential temp.
-REAL, DIMENSION(:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(INOUT):: PRL  ! cloud mixing ratio
-REAL, DIMENSION(:), INTENT(INOUT):: PRI  ! ice   mixing ratio
-REAL, DIMENSION(:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
-
-      END SUBROUTINE TH_R_FROM_THL_RT_1D
-      END INTERFACE
-      END MODULE MODI_TH_R_FROM_THL_RT_1D
-!     ######spl
-      SUBROUTINE TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI,         &
-                                  PRSATW, PRSATI, PRR, PRS, PRG, PRH     )
-!     #################################################################
-!
-!
-!!****  *TH_R_FROM_THL_RT_1D* - computes the non-conservative variables
-!!                          from conservative variables
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!**  METHOD
-!!    ------
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!      Julien PERGAUD      * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original         13/03/06
-!!      S. Riette April 2011 : ice added, allow ZRLTEMP to be negative
-!!                             we use dQsat/dT to help convergence
-!!                             use of optional PRR, PRS, PRG, PRH
-!!      S. Riette Nov 2016: support for HFRAC_ICE='S'
-!!
-!! --------------------------------------------------------------------------
-!
-!*      0. DECLARATIONS
-!          ------------
-!
-USE MODE_COMPUTE_FRAC_ICE, ONLY : COMPUTE_FRAC_ICE
-USE MODD_CST
-USE MODD_DYN_n, ONLY : LOCEAN
-USE MODE_THERMO
-!
-IMPLICIT NONE
-!
-!
-!*      0.1  declarations of arguments
-!
-CHARACTER(LEN=1),   INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:), INTENT(IN) :: PP          ! Pressure
-REAL, DIMENSION(:), INTENT(IN) :: PTHL    ! thetal to transform into th
-REAL, DIMENSION(:),INTENT(IN)  :: PRT    ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:), INTENT(OUT):: PTH    ! th
-REAL, DIMENSION(:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(INOUT):: PRL    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(INOUT):: PRI    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
-!
-!-------------------------------------------------------------------------------
-!
-!       0.2  declaration of local variables
-INTEGER                       :: II ! Loop control
-INTEGER                       :: JITER ! number of iterations
-REAL, DIMENSION(SIZE(PTHL,1))   :: ZEXN
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZRVSAT,ZCPH,ZRLTEMP,ZCPH2
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZT,ZLVOCPEXN,ZLSOCPEXN
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZDRSATODT,ZDRSATODTW,ZDRSATODTI
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZFOESW, ZFOESI
-INTEGER, DIMENSION(SIZE(PTHL,1)) :: IERR
-!----------------------------------------------------------------------------
-!
-!*      1 Initialisation
-!         --------------
-!
-!
-!
-!Number of iterations
-JITER=2
-!
-!Computation of ZCPH2 depending on dummy arguments received
-ZCPH2(:)=0
-IF(PRESENT(PRR)) ZCPH2(:)=ZCPH2(:) + XCL*PRR(:)
-IF(PRESENT(PRS)) ZCPH2(:)=ZCPH2(:) + XCI*PRS(:)
-IF(PRESENT(PRG)) ZCPH2(:)=ZCPH2(:) + XCI*PRG(:)
-IF(PRESENT(PRH)) ZCPH2(:)=ZCPH2(:) + XCI*PRH(:)
-!
-!Computation of an approximate state thanks to PRL and PRI guess
-ZEXN(:)=(PP(:)/XP00) ** (XRD/XCPD)
-ZT(:)=PTHL(:)*ZEXN(:)
-PRV(:)=PRT(:)-PRL(:)-PRI(:)
-ZCPH(:)=XCPD+ XCPV * PRV(:)+ XCL * PRL(:) + XCI * PRI(:) + ZCPH2(:)
-ZLVOCPEXN(:) = (XLVTT + (XCPV-XCL) * (ZT(:)-XTT)) &
-                        /(ZCPH(:)*ZEXN(:))
-ZLSOCPEXN(:) = (XLSTT + (XCPV-XCI) * (ZT(:)-XTT)) &
-                        /(ZCPH(:)*ZEXN(:))
-PTH(:)=PTHL(:)+ZLVOCPEXN(:)*PRL(:)+ZLSOCPEXN(:)*PRI(:)
-!
-!
-!       2 Iteration
-!         ---------
-
-DO II=1,JITER
-  IF (LOCEAN) THEN
-    ZT=PTH                  
-  ELSE
-    ZT(:)=PTH(:)*ZEXN(:)
-  END IF
-  !Computation of liquid/ice fractions
-  PFRAC_ICE(:) = 0.
-  WHERE(PRL(:)+PRI(:) > 1.E-20)
-    PFRAC_ICE(:) = PRI(:) / (PRL(:)+PRI(:))
-  ENDWHERE
-  CALL COMPUTE_FRAC_ICE(HFRAC_ICE,PFRAC_ICE(:),ZT(:), IERR(:))
-  !Computation of Rvsat and dRsat/dT
-  !In this version QSAT, QSATI, DQSAT and DQASATI functions are not used
-  !due to performance issue
-  ZFOESW(:) = MIN(EXP( XALPW - XBETAW/ZT(:) - XGAMW*LOG(ZT(:))  ), 0.99*PP(:))
-  ZFOESI(:) = MIN(EXP( XALPI - XBETAI/ZT(:) - XGAMI*LOG(ZT(:))  ), 0.99*PP(:))
-  PRSATW(:) = XRD/XRV*ZFOESW(:)/PP(:) / (1.+(XRD/XRV-1.)*ZFOESW(:)/PP(:))
-  PRSATI(:) = XRD/XRV*ZFOESI(:)/PP(:) / (1.+(XRD/XRV-1.)*ZFOESI(:)/PP(:))
-  ZDRSATODTW(:) = PRSATW(:) / (1.+(XRD/XRV-1.)*ZFOESW(:)/PP(:) ) &
-                   * (XBETAW/ZT(:)**2 - XGAMW/ZT(:))*(1+PRT(:))
-  ZDRSATODTI(:) = PRSATI(:) / (1.+(XRD/XRV-1.)*ZFOESI(:)/PP(:) ) &
-                   * (XBETAI/ZT(:)**2 - XGAMI/ZT(:))*(1+PRT(:))
-  !PRSATW(:) =  QSAT(ZT(:),PP(:)) !qsatw
-  !PRSATI(:) = QSATI(ZT(:),PP(:)) !qsati
-  !ZDRSATODTW(:) =  DQSAT(ZT(:),PP(:),PRSATW(:))*(1+PRT(:))
-  !ZDRSATODTI(:) = DQSATI(ZT(:),PP(:),PRSATI(:))*(1+PRT(:))
-  PRSATW(:) = PRSATW(:)*(1+PRT(:))
-  PRSATI(:) = PRSATI(:)*(1+PRT(:))
-  ZRVSAT(:) = PRSATW(:)*(1-PFRAC_ICE(:)) + PRSATI(:)*PFRAC_ICE(:)
-  ZDRSATODT(:) = (ZDRSATODTW(:)*(1-PFRAC_ICE(:))+ &
-            & ZDRSATODTI(:)*PFRAC_ICE(:))
-
-  !Computation of new PRL, PRI and PRV
-  !Correction term applied to (PRV(:)-ZRVSAT(:)) is computed assuming that
-  !ZLVOCPEXN, ZLSOCPEXN and ZCPH don't vary to much with T. It takes into account
-  !the variation (estimated linear) of Qsat with T
-  ZRLTEMP(:)=(PRV(:)-ZRVSAT(:))/ &
-                &(1 + ZDRSATODT(:)*ZEXN(:)* &
-                &     (ZLVOCPEXN(:)*(1-PFRAC_ICE(:))+ZLSOCPEXN(:)*PFRAC_ICE(:)))
-  ZRLTEMP(:)=MIN(MAX(-PRL(:)-PRI(:), ZRLTEMP(:)),PRV(:))
-  PRV(:)=PRV(:)-ZRLTEMP(:)
-  PRL(:)=PRL(:)+PRI(:)+ZRLTEMP(:)
-  PRI(:)=PFRAC_ICE(:)     * (PRL(:))
-  PRL(:)=(1-PFRAC_ICE(:)) * (PRT(:) - PRV(:))
-
-  !Computation of Cph (as defined in Meso-NH doc, equation 2.2, to be used with mixing ratios)
-  ZCPH(:)=XCPD+ XCPV * PRV(:)+ XCL * PRL(:) + XCI * PRI(:) + ZCPH2(:)
-
-  !Computation of L/Cph/EXN, then new PTH
-  ZLVOCPEXN(:) = (XLVTT + (XCPV-XCL) * (ZT(:)-XTT)) &
-                    /(ZCPH(:)*ZEXN(:))
-  ZLSOCPEXN(:) = (XLSTT + (XCPV-XCI) * (ZT(:)-XTT)) &
-                    /(ZCPH(:)*ZEXN(:))
-  PTH(:)=PTHL(:)+ZLVOCPEXN(:)*PRL(:)+ZLSOCPEXN(:)*PRI(:)
-
-  !Computation of estimated mixing ration at saturation
-  !To compute the adjustement a first order development was used
-  PRSATW(:)=PRSATW(:) + ZDRSATODTW(:)*(PTH(:)*ZEXN(:)-ZT(:))
-  PRSATI(:)=PRSATI(:) + ZDRSATODTI(:)*(PTH(:)*ZEXN(:)-ZT(:))
-ENDDO
-
-
-END SUBROUTINE TH_R_FROM_THL_RT_1D
diff --git a/src/mesonh/turb/th_r_from_thl_rt_2d.f90 b/src/mesonh/turb/th_r_from_thl_rt_2d.f90
deleted file mode 100644
index d610a72bda073171f6154bdc8e6a6f1e6763300f..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/th_r_from_thl_rt_2d.f90
+++ /dev/null
@@ -1,128 +0,0 @@
-!MNH_LIC Copyright 2006-2019 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!     ######spl
-      MODULE MODI_TH_R_FROM_THL_RT_2D
-      INTERFACE
-      SUBROUTINE TH_R_FROM_THL_RT_2D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI,         &
-                                  PRSATW, PRSATI, PRR, PRS, PRG, PRH     )
-CHARACTER(len=1),     INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:,:), INTENT(IN) :: PP          ! Pressure
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHL    ! thetal to transform into th
-REAL, DIMENSION(:,:),INTENT(IN)  :: PRT    ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:,:), INTENT(OUT):: PTH    ! th
-REAL, DIMENSION(:,:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRL    ! vapor mixing ratio
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRI    ! vapor mixing ratio
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
-
-      END SUBROUTINE TH_R_FROM_THL_RT_2D
-
-      END INTERFACE
-
-      END MODULE MODI_TH_R_FROM_THL_RT_2D
-
-!     ######spl
-      SUBROUTINE TH_R_FROM_THL_RT_2D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI,         &
-                                  PRSATW, PRSATI, PRR, PRS, PRG, PRH     )
-!     #################################################################
-!
-!
-!!****  *TH_R_FROM_THL_RT_2D* - computes the non-conservative variables
-!!                          from conservative variables
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!**  METHOD
-!!    ------
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!      Julien PERGAUD      * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original         13/03/06
-!!      Sébastien Riette April 2011: code moved in th_r_from_thl_rt_3D
-!!
-!! --------------------------------------------------------------------------
-!       
-!*      0. DECLARATIONS
-!          ------------
-!
-!
-USE MODI_TH_R_FROM_THL_RT_3D
-
-IMPLICIT NONE
-!
-!
-!*      0.1  declarations of arguments
-!
-CHARACTER(LEN=1),     INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:,:), INTENT(IN) :: PP     ! Pressure
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHL   ! Liquid pot. temp.
-REAL, DIMENSION(:,:), INTENT(IN) :: PRT    ! Total mixing ratios
-REAL, DIMENSION(:,:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:,:), INTENT(OUT):: PTH    ! Potential temp.
-REAL, DIMENSION(:,:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRL  ! cloud mixing ratio
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRI  ! ice   mixing ratio
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
-
-!
-!-------------------------------------------------------------------------------
-!
-!       0.2  declaration of local variables
-!
-!----------------------------------------------------------------------------
-!
-REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2)) :: ZRR, ZRS, ZRG, ZRH
-INTEGER :: JK
-!----------------------------------------------------------------------------
-!
-!*      1 Initialisation
-!         --------------
-!
-ZRR(:,:)=0.
-ZRS(:,:)=0.
-ZRG(:,:)=0.
-ZRH(:,:)=0.
-IF(PRESENT(PRR)) ZRR(:,:)=PRR(:,:)
-IF(PRESENT(PRS)) ZRS(:,:)=PRS(:,:)
-IF(PRESENT(PRG)) ZRG(:,:)=PRG(:,:)
-IF(PRESENT(PRH)) ZRH(:,:)=PRH(:,:)
-!
-!
-!       2 Call of 1d version
-!         ------------------
-DO JK=1, SIZE(PTHL,2)
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE(:,JK),PP(:,JK),             &
-                                PTHL(:,JK), PRT(:,JK), PTH(:,JK),       &
-                                PRV(:,JK), PRL(:,JK), PRI(:,JK),        &
-                                PRSATW(:,JK), PRSATI(:,JK),                &
-                                ZRR(:,JK), ZRS(:,JK), ZRG(:,JK), ZRH(:,JK))
-ENDDO
-
-
-END SUBROUTINE TH_R_FROM_THL_RT_2D
diff --git a/src/mesonh/turb/th_r_from_thl_rt_3d.f90 b/src/mesonh/turb/th_r_from_thl_rt_3d.f90
deleted file mode 100644
index 0591be50c12cee5be78a7522530ddfa5617a35f1..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/th_r_from_thl_rt_3d.f90
+++ /dev/null
@@ -1,126 +0,0 @@
-!MNH_LIC Copyright 2011-2019 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!     ######spl
-      MODULE MODI_TH_R_FROM_THL_RT_3D
-!     ###############################
-INTERFACE
-!
-      SUBROUTINE TH_R_FROM_THL_RT_3D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI, &
-                                  PRSATW, PRSATI, PRR, PRS, PRG, PRH       )
-
-CHARACTER(len=1),       INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PP          ! Pressure
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHL    ! thetal to transform into th
-REAL, DIMENSION(:,:,:),INTENT(IN)  :: PRT    ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(:,:,:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PTH    ! th
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRL    ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRI    ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
-!
-END SUBROUTINE TH_R_FROM_THL_RT_3D
-!
-END INTERFACE
-!
-END MODULE MODI_TH_R_FROM_THL_RT_3D
-!     ######spl
-      SUBROUTINE TH_R_FROM_THL_RT_3D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI, &
-                                  PRSATW, PRSATI, PRR, PRS, PRG, PRH       )
-!     #################################################################
-!
-!
-!!****  *TH_R_FROM_THL_RT_3D* - computes the non-conservative variables
-!!                          from conservative variables
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!**  METHOD
-!!    ------
-!!
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!      Julien PERGAUD      * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original         15/03/2011
-!!      S. Riette April 2011 : code moved in th_r_from_thl_rt_1d
-!!
-!! --------------------------------------------------------------------------
-!
-!*      0. DECLARATIONS
-!          ------------
-!
-USE MODI_TH_R_FROM_THL_RT_1D
-!
-IMPLICIT NONE
-!
-!
-!*      0.1  declarations of arguments
-!
-CHARACTER(LEN=1),       INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PP          ! Pressure
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHL    ! thetal to transform into th
-REAL, DIMENSION(:,:,:),INTENT(IN)  :: PRT    ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(:,:,:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PTH    ! th
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRL    ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRI    ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
-!
-!-------------------------------------------------------------------------------
-!
-!       0.2  declaration of local variables
-REAL, DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2),SIZE(PTHL,3)) :: ZRR, ZRS, ZRG, ZRH
-INTEGER :: JJ, JK
-!----------------------------------------------------------------------------
-!
-!*      1 Initialisation
-!         --------------
-!
-ZRR(:,:,:)=0.
-ZRS(:,:,:)=0.
-ZRG(:,:,:)=0.
-ZRH(:,:,:)=0.
-IF(PRESENT(PRR)) ZRR(:,:,:)=PRR(:,:,:)
-IF(PRESENT(PRS)) ZRS(:,:,:)=PRS(:,:,:)
-IF(PRESENT(PRG)) ZRG(:,:,:)=PRG(:,:,:)
-IF(PRESENT(PRH)) ZRH(:,:,:)=PRH(:,:,:)
-!
-!
-!       2 Call of 1d version
-!         ------------------
-DO JK=1, SIZE(PTHL,3)
-  DO JJ=1, SIZE(PTHL,2)
-    CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE(:,JJ,JK),PP(:,JJ,JK),             &
-                                  PTHL(:,JJ,JK), PRT(:,JJ,JK), PTH(:,JJ,JK),       &
-                                  PRV(:,JJ,JK), PRL(:,JJ,JK), PRI(:,JJ,JK),        &
-                                  PRSATW(:,JJ,JK), PRSATI(:,JJ,JK),                &
-                                  ZRR(:,JJ,JK), ZRS(:,JJ,JK), ZRG(:,JJ,JK), ZRH(:,JJ,JK))
-  ENDDO
-ENDDO
-
-END SUBROUTINE TH_R_FROM_THL_RT_3D
diff --git a/src/mesonh/turb/tke_eps_sources.f90 b/src/mesonh/turb/tke_eps_sources.f90
deleted file mode 100644
index 4efe246beff2b829a835448166f9ad33bd6ea51f..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/tke_eps_sources.f90
+++ /dev/null
@@ -1,485 +0,0 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!     ###########################
-      MODULE MODI_TKE_EPS_SOURCES
-!     ###########################
-INTERFACE
-!
-      SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP,PTRH,  &
-                      PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,                  &
-                      PTSTEP,PIMPL,PEXPL,                                   &
-                      HTURBLEN,HTURBDIM,                                    &
-                      TPFILE,OTURB_DIAG,                                    &
-                      PTP,PRTKES,PRTKESM, PRTHLS,PCOEF_DISS,PTR,PDISS       )
-!
-USE MODD_IO, ONLY: TFILEDATA
-!
-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(:,:,:),  INTENT(IN)   ::  PTKEM        ! TKE at t-deltat
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PLM          ! mixing length         
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PLEPS        ! dissipative length
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRHODJ       ! density * grid volume
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                       ! metric coefficients
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PZZ          ! physical height w-pt
-REAL,                    INTENT(IN)   ::  PTSTEP       ! Time step 
-REAL,                    INTENT(IN)   ::  PEXPL, PIMPL ! Coef. temporal. disc.
-CHARACTER(len=4),        INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                       ! turbulence scheme
-CHARACTER(len=4),        INTENT(IN)   ::  HTURBLEN     ! kind of mixing length
-TYPE(TFILEDATA),         INTENT(IN)   ::  TPFILE       ! Output file
-LOGICAL,                 INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                  ! diagnostic fields in the syncronous FM-file
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PDP          ! Dyn. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTRH
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTP          ! Ther. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTKES       ! RHOD * Jacobian *
-                                                       ! TKE at t+deltat
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRTKESM      ! Advection source 
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTHLS       ! Source of Theta_l
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PCOEF_DISS   ! 1/(Cph*Exner)
-REAL, DIMENSION(:,:,:),  INTENT(OUT)  ::  PTR          ! Transport prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(OUT)  ::  PDISS        ! Dissipati prod. of TKE
-!
-!
-!
-END SUBROUTINE TKE_EPS_SOURCES
-!
-END INTERFACE
-!
-END MODULE MODI_TKE_EPS_SOURCES
-!
-!     ##################################################################
-      SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP,  &
-                      PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,        &
-                      PTSTEP,PIMPL,PEXPL,                              &
-                      HTURBLEN,HTURBDIM,                               &
-                      TPFILE,OTURB_DIAG,                               &
-                      PTP,PRTKES,PRTKESM, PRTHLS,PCOEF_DISS,PTR,PDISS  )
-!     ##################################################################
-!
-!
-!!****  *TKE_EPS_SOURCES* - routine to compute the sources of the turbulent 
-!!      evolutive variables: TKE and its dissipation when it is taken into 
-!!      account. The contribution to the heating of tke dissipation is computed.
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this routine is to compute the sources necessary for
-!     the evolution of the turbulent kinetic energy and its dissipation 
-!     if necessary.
-!
-!!**  METHOD
-!!    ------
-!!      The vertical turbulent flux is computed in an off-centered 
-!!    implicit scheme (a Crank-Nicholson type with coefficients different 
-!!    than 0.5), which allows to vary the degree of implicitness of the 
-!!    formulation.
-!!      In high resolution, the horizontal transport terms are also
-!!    calculated, but explicitly. 
-!!      The evolution of the dissipation as a variable is made if 
-!!    the parameter HTURBLEN is set equal to KEPS. The same reasoning 
-!!    made for TKE applies.
-!!
-!!    EXTERNAL
-!!    --------
-!!      GX_U_M,GY_V_M,GZ_W_M
-!!      GX_M_U,GY_M_V          :  Cartesian vertical gradient operators
-!!
-!!      MXF,MXM.MYF,MYM,MZF,MZM:  Shuman functions (mean operators)
-!!      DZF                    :  Shuman functions (difference operators)     
-!!
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      Module MODD_CST : contains physical constants
-!!
-!!           XG         : gravity constant
-!!
-!!      Module MODD_CTURB: contains the set of constants for
-!!                        the turbulence scheme
-!!
-!!           XCET,XCED  : transport and dissipation cts. for the TKE
-!!           XCDP,XCDD,XCDT: constants from the parameterization of
-!!                        the K-epsilon equation
-!!           XTKEMIN,XEPSMIN : minimum values for the TKE and its
-!!                        dissipation
-!!
-!!      Module MODD_PARAMETERS: 
-!!
-!!           JPVEXT
-!!      Module MODD_BUDGET:
-!!         NBUMOD       : model in which budget is calculated
-!!         CBUTYPE      : type of desired budget
-!!                          'CART' for cartesian box configuration
-!!                          'MASK' for budget zone defined by a mask 
-!!                          'NONE'  ' for no budget
-!!         LBU_RTKE     : logical for budget of RTKE (turbulent kinetic energy)
-!!                        .TRUE. = budget of RTKE       
-!!                        .FALSE. = no budget of RTKE
-!!
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book 2 of documentation (routine TKE_EPS_SOURCES)
-!!      Book 1 of documentation (Chapter: Turbulence)
-!!
-!!    AUTHOR
-!!    ------
-!!      Joan Cuxart             * INM and Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original       August 23, 1994
-!!      Modifications: Feb 14, 1995 (J.Cuxart and J.Stein)
-!!                                  Doctorization and Optimization
-!!                     June 29, 1995 (J.Stein) TKE budget
-!!                     June 28, 1995 (J.Cuxart) Add LES tools
-!!      Modifications: February 29, 1996 (J. Stein) optimization
-!!      Modifications: May 6, 1996 (N. Wood) Extend some loops over
-!!                                              the outer points
-!!      Modifications: August 30, 1996 (P. Jabouille)  calcul ZFLX at the
-!!                                                      IKU level
-!!                     October 10, 1996 (J.Stein)  set Keff at t-deltat
-!!                     Oct 8, 1996 (Cuxart,Sanchez) Var.LES: XETR_TF,XDISS_TF
-!!                     December 20, 1996 (J.-P. Pinty) update the CALL BUDGET
-!!                     November 24, 1997 (V. Masson) bug in <v'e>
-!!                                                   removes the DO loops
-!!                     Augu. 9, 1999 (J.Stein) TKE budget correction
-!!                     Mar 07  2001 (V. Masson and J. Stein) remove the horizontal 
-!!                                         turbulent transports of Tke computation
-!!                     Nov 06, 2002 (V. Masson) LES budgets
-!!                     July 20, 2003 (J.-P. Pinty P Jabouille) add the dissipative heating
-!!                     May   2006    Remove KEPS
-!!                     October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
-!!                                              change of YCOMMENT
-!!                     2012-02 Y. Seity,  add possibility to run with reversed 
-!!                                    vertical levels
-!!                     2015-01 (J. Escobar) missing get_halo(ZRES) for JPHEXT<> 1 
-!!     J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
-!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!  P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
-!  P. Wautelet    02/2020: use the new data structures and subroutines for budgets
-! --------------------------------------------------------------------------
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-USE MODD_ARGSLIST_ll,    ONLY: LIST_ll
-use modd_budget,         only: lbudget_tke, lbudget_th, NBUDGET_TKE, NBUDGET_TH, tbudgets
-USE MODD_CONF
-USE MODD_CST
-USE MODD_CTURB
-USE MODD_DIAG_IN_RUN,    ONLY: LDIAG_IN_RUN, XCURRENT_TKE_DISS
-use modd_field,          only: tfielddata, TYPEREAL
-USE MODD_IO,             ONLY: TFILEDATA
-USE MODD_LES
-USE MODD_PARAMETERS
-!
-use mode_budget,         only: Budget_store_add, Budget_store_end, Budget_store_init
-USE MODE_IO_FIELD_WRITE, only: IO_Field_write
-USE MODE_ll
-!
-USE MODI_GET_HALO
-USE MODI_GRADIENT_M
-USE MODI_GRADIENT_U
-USE MODI_GRADIENT_V
-USE MODI_GRADIENT_W
-USE MODI_LES_MEAN_SUBGRID
-USE MODI_SHUMAN
-USE MODI_TRIDIAG_TKE
-!
-IMPLICIT NONE
-!
-!
-!*       0.1  declarations of arguments
-!
-!
-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(:,:,:),  INTENT(IN)   ::  PTKEM        ! TKE at t-deltat
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PLM          ! mixing length         
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PLEPS        ! dissipative length
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRHODJ       ! density * grid volume
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                       ! metric coefficients
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PZZ          ! physical height w-pt
-REAL,                    INTENT(IN)   ::  PTSTEP       ! Time step 
-REAL,                    INTENT(IN)   ::  PEXPL, PIMPL ! Coef. temporal. disc.
-CHARACTER(len=4),        INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                       ! turbulence scheme
-CHARACTER(len=4),        INTENT(IN)   ::  HTURBLEN     ! kind of mixing length
-TYPE(TFILEDATA),         INTENT(IN)   ::  TPFILE       ! Output file
-LOGICAL,                 INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                  ! diagnostic fields in the syncronous FM-file
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PDP          ! Dyn. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTRH
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTP          ! Ther. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTKES       ! RHOD * Jacobian *
-                                                       ! TKE at t+deltat
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTHLS       ! Source of Theta_l
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PCOEF_DISS   ! 1/(Cph*Exner)
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRTKESM      ! Advection source 
-REAL, DIMENSION(:,:,:),  INTENT(OUT)  ::  PTR          ! Transport prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(OUT)  ::  PDISS        ! Dissipati prod. of TKE
-!
-!
-!
-!*       0.2  declaration of local variables
-!
-REAL, DIMENSION(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3))::         &
-       ZA,       & ! under diagonal elements of the tri-diagonal matrix involved
-                   ! in the temporal implicit scheme
-       ZRES,     & ! treated variable at t+ deltat when the turbu-
-                   ! lence is the only source of evolution added to the ones
-                   ! considered in ZSOURCE. This variable is also used to
-                   ! temporarily store some diagnostics stored in FM file
-       ZFLX,     & ! horizontal or vertical flux of the treated variable
-       ZSOURCE,  & ! source of evolution for the treated variable
-       ZKEFF       ! effectif diffusion coeff = LT * SQRT( TKE )
-!LOGICAL,DIMENSION(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3)) :: GTKENEG
-!                   ! 3D mask .T. if TKE < XTKEMIN
-INTEGER             :: IIB,IIE,IJB,IJE,IKB,IKE
-                                    ! Index values for the Beginning and End
-                                    ! mass points of the domain 
-INTEGER             :: IIU,IJU,IKU  ! array size in the 3 dimensions 
-!
-TYPE(LIST_ll), POINTER :: TZFIELDDISS_ll ! list of fields to exchange
-INTEGER                :: IINFO_ll       ! return code of parallel routine
-TYPE(TFIELDDATA) :: TZFIELD
-!
-!----------------------------------------------------------------------------
-NULLIFY(TZFIELDDISS_ll)
-!
-!*       1.   PRELIMINARY COMPUTATIONS
-!             ------------------------
-!
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
-IIU=SIZE(PTKEM,1)
-IJU=SIZE(PTKEM,2)
-IKB=KKA+JPVEXT_TURB*KKL
-IKE=KKU-JPVEXT_TURB*KKL
-!
-! compute the effective diffusion coefficient at the mass point
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) 
-
-if (lbudget_th)  call Budget_store_init( tbudgets(NBUDGET_TH),  'DISSH', prthls(:, :, :) )
-
-!----------------------------------------------------------------------------
-!
-!*       2.   TKE EQUATION  
-!             ------------
-!
-!*       2.1  Horizontal turbulent explicit transport
-!
-!
-! Complete the sources of TKE with the horizontal turbulent explicit transport
-!
-IF (HTURBDIM=='3DIM') THEN
-  PTR=PTRH
-ELSE
-  PTR=0.
-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))
-!
-! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
-! + (Dynamical Production) + (Thermal Production) - (dissipation) 
-ZFLX(:,:,:) = XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:)
-ZSOURCE(:,:,:) = ( PRTKES(:,:,:) +  PRTKESM(:,:,:) ) / PRHODJ(:,:,:) &
-   - PTKEM(:,:,:) / PTSTEP &
-   + PDP(:,:,:) + PTP(:,:,:) + PTR(:,:,:) - PEXPL * ZFLX(:,:,:) * PTKEM(:,:,:)
-!
-!*       2.2  implicit vertical TKE transport
-!
-!
-! Compute the vector giving the elements just under the diagonal for the 
-! matrix inverted in TRIDIAG 
-!
-ZA(:,:,:)     = - PTSTEP * XCET * &
-                MZM(ZKEFF) * MZM(PRHODJ) / PDZZ**2
-!
-! Compute TKE at time t+deltat: ( stored in ZRES )
-!
-CALL TRIDIAG_TKE(KKA,KKU,KKL,PTKEM,ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,&
-            & ZSOURCE,PTSTEP*ZFLX,ZRES)
-CALL GET_HALO(ZRES)
-!
-!* diagnose the dissipation
-!
-IF (LDIAG_IN_RUN) THEN
-  XCURRENT_TKE_DISS = ZFLX(:,:,:) * PTKEM(:,:,:) &
-                                  *(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:))
-  CALL ADD3DFIELD_ll( TZFIELDDISS_ll, XCURRENT_TKE_DISS, 'TKE_EPS_SOURCES::XCURRENT_TKE_DISS' )
-  CALL UPDATE_HALO_ll(TZFIELDDISS_ll,IINFO_ll)
-  CALL CLEANLIST_ll(TZFIELDDISS_ll)
-ENDIF
-!
-! TKE must be greater than its minimum value
-!
-! CL : Now done at the end of the time step in ADVECTION_METSV
-!GTKENEG =  ZRES <= XTKEMIN 
-!WHERE ( GTKENEG ) 
-!  ZRES = XTKEMIN
-!END WHERE
-!
-IF ( LLES_CALL .OR.                         &
-     (OTURB_DIAG .AND. tpfile%lopened)  ) THEN
-!
-! Compute the cartesian vertical flux of TKE in ZFLX
-!
-
-  ZFLX(:,:,:)   = - XCET * MZM(ZKEFF) *   &
-                  DZM(PIMPL * ZRES + PEXPL * PTKEM ) / PDZZ
-!
-  ZFLX(:,:,IKB) = 0.
-  ZFLX(:,:,KKA) = 0.
-!
-! Compute the whole turbulent TRansport of TKE:
-!
-  PTR(:,:,:)= PTR - DZF( MZM(PRHODJ) * ZFLX / PDZZ ) /PRHODJ
-!
-! Storage in the LES configuration
-!
-  IF (LLES_CALL) THEN
-    CALL LES_MEAN_SUBGRID( MZF(ZFLX), X_LES_SUBGRID_WTke )
-    CALL LES_MEAN_SUBGRID( -PTR, X_LES_SUBGRID_ddz_WTke )
-  END IF
-!
-END IF
-!
-!*       2.4  stores the explicit sources for budget purposes
-!
-if (lbudget_tke) then
-  ! Dynamical production
-  call Budget_store_add( tbudgets(NBUDGET_TKE), 'DP', pdp(:, :, :) * prhodj(:, :, :) )
-  ! Thermal production
-  call Budget_store_add( tbudgets(NBUDGET_TKE), 'TP', ptp(:, :, :) * prhodj(:, :, :) )
-  ! Dissipation
-  call Budget_store_add( tbudgets(NBUDGET_TKE), 'DISS', -xced * sqrt( ptkem(:, :, :) ) / pleps(:, :, :) &
-                         * ( pexpl * ptkem(:, :, :) + pimpl * zres(:, :, :) ) * prhodj(:, :, :) )
-end if
-!
-!*       2.5  computes the final RTKE and stores the whole turbulent transport
-!              with the removal of the advection part 
-
-if (lbudget_tke) then
-  !Store the previous source terms in prtkes before initializing the next one
-  PRTKES(:,:,:) = PRTKES(:,:,:) + PRHODJ(:,:,:) *                                                           &
-                  ( PDP(:,:,:) + PTP(:,:,:)                                                                 &
-                    - XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * ( PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:) ) )
-
-  call Budget_store_init( tbudgets(NBUDGET_TKE), 'TR', prtkes(:, :, :) )
-end if
-
-PRTKES(:,:,:) = ZRES(:,:,:) * PRHODJ(:,:,:) / PTSTEP -  PRTKESM(:,:,:)
-!
-! stores the whole turbulent transport
-!
-if (lbudget_tke) call Budget_store_end( tbudgets(NBUDGET_TKE), 'TR', prtkes(:, :, :) )
-
-!----------------------------------------------------------------------------
-!
-!*       3.   COMPUTE THE DISSIPATIVE HEATING
-!             -------------------------------
-!
-PRTHLS(:,:,:) = PRTHLS(:,:,:) + XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * &
-                (PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) * PRHODJ(:,:,:) * PCOEF_DISS(:,:,:)
-
-if (lbudget_th) call Budget_store_end( tbudgets(NBUDGET_TH), 'DISSH', prthls(:, :, :) )
-
-!----------------------------------------------------------------------------
-!
-!*       4.   STORES SOME DIAGNOSTICS
-!             -----------------------
-!
-PDISS(:,:,:) =  -XCED * (PTKEM(:,:,:)**1.5) / PLEPS(:,:,:)
-!
-IF ( OTURB_DIAG .AND. tpfile%lopened ) THEN
-!
-! stores the dynamic production 
-!
-  TZFIELD%CMNHNAME   = 'DP'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'DP'
-  TZFIELD%CUNITS     = 'm2 s-3'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_DP'
-  TZFIELD%NGRID      = 1
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PDP)
-!
-! stores the thermal production 
-!
-  TZFIELD%CMNHNAME   = 'TP'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'TP'
-  TZFIELD%CUNITS     = 'm2 s-3'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_TP'
-  TZFIELD%NGRID      = 1
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PTP)
-!
-! stores the whole turbulent transport
-!
-  TZFIELD%CMNHNAME   = 'TR'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'TR'
-  TZFIELD%CUNITS     = 'm2 s-3'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_TR'
-  TZFIELD%NGRID      = 1
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PTR)
-!
-! stores the dissipation of TKE 
-!
-  TZFIELD%CMNHNAME   = 'DISS'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'DISS'
-  TZFIELD%CUNITS     = 'm2 s-3'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_DISS'
-  TZFIELD%NGRID      = 1
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PDISS)
-END IF
-!
-! Storage in the LES configuration of the Dynamic Production of TKE and
-! the dissipation of TKE 
-! 
-IF (LLES_CALL ) THEN
-  CALL LES_MEAN_SUBGRID( PDISS, X_LES_SUBGRID_DISS_Tke )
-END IF
-!
-!----------------------------------------------------------------------------
-! 
-!
-!----------------------------------------------------------------------------
-!
-END SUBROUTINE TKE_EPS_SOURCES
diff --git a/src/mesonh/turb/turb.f90 b/src/mesonh/turb/turb.f90
index fff7f9fed9f377f94378cedf8275d8bf8aab960c..dc38dc779f8a6973af95ac7c77b62c83966d7cd4 100644
--- a/src/mesonh/turb/turb.f90
+++ b/src/mesonh/turb/turb.f90
@@ -4,11 +4,6 @@
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !    ################ 
-     MODULE MODI_TURB  
-!    ################ 
-!
-INTERFACE
-!
       SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,      &
               & KSPLIT,KMODEL_CL,                                     &
               & OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,               &
@@ -23,139 +18,9 @@ INTERFACE
               & PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,                &
               & PTHLT,PRT,                                            &
               & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES,PRTKEMS,PSIGS,&
-              & PFLXZTHVMF,PWTH,PWRC,PWSV,PDYP,PTHP,PTDIFF,PTDISS,PLEM,&
-              & TBUDGETS, KBUDGETS                                   )
-
-!
-USE MODD_IO, ONLY: TFILEDATA
-USE MODD_BUDGET, ONLY: TBUDGETDATA
-!
-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=AR
-INTEGER,                INTENT(IN)   :: KMI           ! model index number  
-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.
-CHARACTER(LEN=*),DIMENSION(:),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
-INTEGER,                INTENT(IN)   :: KSPLIT        ! number of time-splitting
-INTEGER,                INTENT(IN)   :: KMODEL_CL     ! model number for cloud mixing length
-LOGICAL,                INTENT(IN)   ::  OTURB_FLX    ! switch to write the
-                                 ! turbulent fluxes in the syncronous FM-file
-LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                 ! diagnostic fields in the syncronous FM-file
-LOGICAL,                INTENT(IN)   ::  OSUBG_COND   ! switch for SUBGrid 
-                                 ! CONDensation
-LOGICAL,                INTENT(IN)   ::  ORMC01       ! switch for RMC01 lengths in SBL
-CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                      ! turbulence scheme
-CHARACTER(len=4),       INTENT(IN)   ::  HTURBLEN     ! kind of mixing length
-CHARACTER(len=4),       INTENT(IN)   ::  HTOM         ! kind of Third Order Moment
-CHARACTER(len=4),       INTENT(IN)   ::  HTURBLEN_CL  ! kind of cloud mixing length
-                                                      ! surface friction flux
-REAL,                   INTENT(IN)   ::  PIMPL        ! degree of implicitness
-CHARACTER (LEN=4),      INTENT(IN)   ::  HCLOUD       ! Kind of microphysical scheme
-REAL,                   INTENT(IN)   ::  PTSTEP       ! timestep 
-TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDXX,PDYY,PDZZ,PDZX,PDZY
-                                        ! metric coefficients
-REAL, DIMENSION(:,:,:), INTENT(IN)   :: PZZ       !  physical distance 
-! between 2 succesive grid points along the K direction
-REAL, DIMENSION(:,:),   INTENT(IN)      ::  PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
-! Director Cosinus along x, y and z directions at surface w-point
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCOSSLOPE       ! cosinus of the angle
-                                 ! between i and the slope vector
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE       ! sinus of the angle
-                                 ! between i and the slope vector
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODJ    ! dry density * Grid size
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  MFMOIST ! moist mass flux dual scheme
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PTHVREF   ! Virtual Potential
-                                        ! Temperature of the reference state
-!
-REAL, DIMENSION(:,:),   INTENT(IN)      ::  PSFTH,PSFRV,   &
-! normal surface fluxes of theta and Rv 
-                                            PSFU,PSFV
-! normal surface fluxes of (u,v) parallel to the orography 
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PSFSV
-! normal surface fluxes of Scalar var. 
-!
-!    prognostic variables at t- deltat
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABST      ! Pressure at time t
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUT,PVT,PWT ! wind components
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKET       ! TKE
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVT        ! passive scal. var.
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PSRCT       ! Second-order flux
-                      ! s'rc'/2Sigma_s2 at time t-1 multiplied by Lambda_3
-REAL, DIMENSION(:,:),     INTENT(INOUT) :: PBL_DEPTH  ! BL depth for TOMS
-REAL, DIMENSION(:,:),     INTENT(INOUT) :: PSBL_DEPTH ! SBL depth for RMC01
-!
-!
-!    variables for cloud mixing length
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PCEI ! Cloud Entrainment instability
-                                                 ! index to emphasize localy 
-                                                 ! turbulent fluxes
-REAL, INTENT(IN)      ::  PCEI_MIN ! minimum threshold for the instability index CEI
-REAL, INTENT(IN)      ::  PCEI_MAX ! maximum threshold for the instability index CEI
-REAL, INTENT(IN)      ::  PCOEF_AMPL_SAT ! saturation of the amplification coefficient
-!   thermodynamical variables which are transformed in conservative var.
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PTHLT       ! conservative pot. temp.
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRT         ! water var.  where 
-                             ! PRT(:,:,:,1) is the conservative mixing ratio        
-!
-! sources of momentum, conservative potential temperature, Turb. Kin. Energy, 
-! TKE dissipation
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PRUS,PRVS,PRWS,PRTHLS,PRTKES
-! Source terms for all water kinds, PRRS(:,:,:,1) is used for the conservative
-! mixing ratio
-REAL, DIMENSION(:,:,:),   INTENT(IN)    ::  PRTKEMS
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRRS 
-! Source terms for all passive scalar variables
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRSVS
-! Sigma_s at time t+1 : square root of the variance of the deviation to the 
-! saturation 
-REAL, DIMENSION(:,:,:), INTENT(OUT)     ::  PSIGS
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PFLXZTHVMF 
-!                                           MF contribution for vert. turb. transport
-!                                           used in the buoy. prod. of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWTH       ! heat flux
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWRC       ! cloud water flux
-REAL, DIMENSION(:,:,:,:),INTENT(OUT) :: PWSV       ! scalar flux
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDYP  ! Dynamical production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHP  ! Thermal production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTDIFF   ! Transport production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTDISS ! Dissipation of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PLEM  ! Mixing length
-TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
-INTEGER, INTENT(IN) :: KBUDGETS
-! length scale from vdfexcu
-REAL, DIMENSION(:,:,:), INTENT(IN)    :: PLENGTHM, PLENGTHH
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE TURB
-!
-END INTERFACE
-!
-END MODULE MODI_TURB
-!
-!     #################################################################
-      SUBROUTINE TURB(KKA, KKU, KKL, KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,   &
-                KSPLIT,KMODEL_CL,                                     &
-                OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,               &
-                HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HCLOUD,PIMPL,      &
-                PTSTEP,TPFILE,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,           &
-                PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE,    &
-                PRHODJ,PTHVREF,                                       &
-                PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
-                PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT,                  &
-                PLENGTHM,PLENGTHH,MFMOIST,                            &
-                PBL_DEPTH, PSBL_DEPTH,                                &
-                PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,                &
-                PTHLT,PRT,                                            &
-                PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES,PRTKEMS,PSIGS,&
-                PFLXZTHVMF,PWTH,PWRC,PWSV,PDYP,PTHP,PTDIFF,PTDISS,PLEM,&
-                TBUDGETS, KBUDGETS                                    )
+              & PFLXZTHVMF,PWTH,PWRC,PWSV,PDYP,PTHP,PTDIFF,PTDISS,&
+              & TBUDGETS, KBUDGETS                                   ,&
+              & PTR,PDISS,PEDR,PLEM                            )
 !     #################################################################
 !
 !
@@ -366,7 +231,7 @@ USE MODD_BUDGET, ONLY: LBUDGET_U,  LBUDGET_V,  LBUDGET_W,  LBUDGET_TH, LBUDGET_R
                             LBUDGET_RR, LBUDGET_RI, LBUDGET_RS, LBUDGET_RG, LBUDGET_RH, LBUDGET_SV,  &
                             NBUDGET_U,  NBUDGET_V,  NBUDGET_W,  NBUDGET_TH, NBUDGET_RV, NBUDGET_RC,  &
                             NBUDGET_RR, NBUDGET_RI, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1, &
-                            TBUDGETS
+                            TBUDGETDATA
 USE MODD_CONF
 USE MODD_CST
 USE MODD_CTURB
@@ -379,17 +244,17 @@ USE MODD_PARAMETERS, ONLY: JPVEXT_TURB
 USE MODD_PARAM_LIMA
 USE MODD_TURB_n, ONLY: XCADAP
 !
+USE MODE_BL89, ONLY: BL89
 USE MODI_GRADIENT_M
 USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
-USE MODI_BL89
 USE MODE_TURB_VER, ONLY : TURB_VER
 USE MODI_ROTATE_WIND
 USE MODI_TURB_HOR_SPLT 
 USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
 USE MODI_SHUMAN
 USE MODI_LES_MEAN_SUBGRID
-USE MODI_RMC01
+USE MODE_RMC01, ONLY: RMC01
 USE MODI_GRADIENT_W
 USE MODE_TM06, ONLY: TM06
 USE MODI_UPDATE_LM
@@ -398,10 +263,10 @@ USE MODI_GET_HALO
 USE MODE_BUDGET,         ONLY: BUDGET_STORE_INIT, BUDGET_STORE_END
 USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 USE MODE_SBL
-USE MODE_SOURES_NEG_CORRECT, ONLY: SOURCES_NEG_CORRECT
+USE MODE_SOURCES_NEG_CORRECT, ONLY: SOURCES_NEG_CORRECT
 !
-USE MODI_EMOIST
-USE MODI_ETHETA
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
 !
 USE MODI_SECOND_MNH
 !
@@ -509,7 +374,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDYP  ! Dynamical production of TKE
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHP  ! Thermal production of TKE
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTDIFF   ! Transport production of TKE
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTDISS ! Dissipation of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PLEM  ! Mixing length
 !
 TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
 INTEGER, INTENT(IN) :: KBUDGETS
@@ -517,16 +381,22 @@ INTEGER, INTENT(IN) :: KBUDGETS
 ! length scale from vdfexcu
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PLENGTHM, PLENGTHH
 !
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PTR   ! Transport production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PDISS ! Dissipation of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PEDR  ! EDR
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PLEM  ! Mixing length
+!
+!
 !-------------------------------------------------------------------------------
 !
 !       0.2  declaration of local variables
 !
-REAL, ALLOCATABLE, DIMENSION(:,:,:) ::&
+REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) ::     &
           ZCP,                        &  ! Cp at t-1
           ZEXN,                       &  ! EXN at t-1
           ZT,                         &  ! T at t-1
           ZLOCPEXNM,                  &  ! Lv/Cp/EXNREF at t-1
-          ZLMW,                       &  ! Turbulent mixing length (work array)
+          ZLM,ZLMW,                   &  ! Turbulent mixing length (+ work array)
           ZLEPS,                      &  ! Dissipative length
           ZTRH,                       &  ! Dynamic and Thermal Production of TKE
           ZATHETA,ZAMOIST,            &  ! coefficients for s = f (Thetal,Rnp)
@@ -535,9 +405,9 @@ REAL, ALLOCATABLE, DIMENSION(:,:,:) ::&
           ZMWTH,ZMWR,ZMTH2,ZMR2,ZMTHR,&  ! 3rd order moments
           ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,&  ! opposite of verticale derivate of 3rd order moments
           ZTHLM, ZTR, ZDISS              ! initial potential temp.
-REAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::     &
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),SIZE(PRT,4)) ::     &
           ZRM                            ! initial mixing ratio 
-REAL, ALLOCATABLE, DIMENSION(:,:) ::  ZTAU11M,ZTAU12M,  &
+REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2)) ::  ZTAU11M,ZTAU12M,  &
                                                  ZTAU22M,ZTAU33M,  &
             ! tangential surface fluxes in the axes following the orography
                                                  ZUSLOPE,ZVSLOPE,  &
@@ -571,58 +441,25 @@ REAL                :: ZALPHA       ! work coefficient :
 !                                   !   BL89 mixing length near the surface
 !
 REAL :: ZTIME1, ZTIME2
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)):: ZTT,ZEXNE,ZLV,ZLS,ZCPH,ZCOR
 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3))::  ZSHEAR, ZDUDZ, ZDVDZ
 TYPE(TFIELDDATA) :: TZFIELD
 !
-!------------------------------------------------------------------------------------------
-ALLOCATE (                                                        &
-          ZCP(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),         &  
-          ZEXN(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        &  
-          ZT(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),          &  
-          ZLOCPEXNM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),   & 
-          ZLMW(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        & 
-          ZLEPS(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &  
-          ZTRH(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        &  
-          ZATHETA(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),     &
-          ZAMOIST(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),     & 
-          ZCOEF_DISS(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),  &  
-          ZFRAC_ICE(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),   &  
-          ZMWTH(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &
-          ZMWR(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        &
-          ZMTH2(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &
-          ZMR2(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        &
-          ZMTHR(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &  
-          ZFWTH(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &
-          ZFWR(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        &
-          ZFTH2(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &
-          ZFR2(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),        &
-          ZFTHR(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)),       &  
-          ZTHLM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3))        )                        
-
-ALLOCATE ( ZRM(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),SIZE(PRT,4))   )
-
-ALLOCATE ( &
-           ZTAU11M(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZTAU12M(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZTAU22M(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZTAU33M(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZUSLOPE(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZVSLOPE(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZCDUEFF(SIZE(PTHLT,1),SIZE(PTHLT,2)),  &
-           ZUSTAR(SIZE(PTHLT,1),SIZE(PTHLT,2)),   &
-           ZLMO(SIZE(PTHLT,1),SIZE(PTHLT,2)),     &
-           ZRVM(SIZE(PTHLT,1),SIZE(PTHLT,2)),     &
-           ZSFRV(SIZE(PTHLT,1),SIZE(PTHLT,2))     )
-
-!------------------------------------------------------------------------------------------
-!
 !*      1.PRELIMINARIES
 !         -------------
 !
 !*      1.1 Set the internal domains, ZEXPL 
 !
 !
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB',0,ZHOOK_HANDLE)
+!
+IF (LHARAT .AND. HTURBDIM /= '1DIM') THEN
+  CALL ABOR1('LHARATU only implemented for option HTURBDIM=1DIM!')
+ENDIF
+IF (LHARAT .AND. LLES_CALL) THEN
+  CALL ABOR1('LHARATU not implemented for option LLES_CALL')
+ENDIF
+
 IKT=SIZE(PTHLT,3)          
 IKTB=1+JPVEXT_TURB              
 IKTE=IKT-JPVEXT_TURB
@@ -725,7 +562,7 @@ IF (KRRL >=1) THEN
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 3
     TZFIELD%LTIMEDEP   = .TRUE.
-    CALL IO_Field_write(TPFILE,TZFIELD,ZATHETA)
+    CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZATHETA)
 ! 
     TZFIELD%CMNHNAME   = 'AMOIST'
     TZFIELD%CSTDNAME   = ''
@@ -737,7 +574,7 @@ IF (KRRL >=1) THEN
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 3
     TZFIELD%LTIMEDEP   = .TRUE.
-    CALL IO_Field_write(TPFILE,TZFIELD,ZAMOIST)
+    CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZAMOIST)
   END IF
 !
 ELSE
@@ -772,6 +609,8 @@ END IF
 !          -----------------------------------------
 !
 !
+IF (.NOT. LHARAT) THEN
+
 SELECT CASE (HTURBLEN)
 !
 !*      3.1 BL89 mixing length
@@ -779,7 +618,7 @@ SELECT CASE (HTURBLEN)
 
   CASE ('BL89')
     ZSHEAR=0.
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM)
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM)
 !
 !*      3.2 RM17 mixing length
 !           ------------------
@@ -788,7 +627,7 @@ SELECT CASE (HTURBLEN)
     ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ)))
     ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ)))
     ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM)
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM)
 !
 !*      3.3 Grey-zone combined RM17 & Deardorff mixing lengths 
 !           --------------------------------------------------
@@ -797,7 +636,7 @@ SELECT CASE (HTURBLEN)
     ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ)))
     ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ)))
     ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM)
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM)
 
     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.
@@ -805,37 +644,37 @@ SELECT CASE (HTURBLEN)
     ! For LES grid meshes, this is equivalent to Deardorff : the base mixing lentgh is the horizontal grid mesh, 
     !                      and it is limited by a stability-based length (RM17), as was done in Deardorff length (but taking into account shear as well)
     ! For grid meshes in the grey zone, then this is the smaller of the two.
-    PLEM = MIN(PLEM,XCADAP*ZLMW)
+    ZLM = MIN(ZLM,XCADAP*ZLMW)
 !
 !*      3.4 Delta mixing length
 !           -------------------
 !
   CASE ('DELT')
-    CALL DELT(PLEM,ODZ=.TRUE.)
+    CALL DELT(ZLM,ODZ=.TRUE.)
 !
 !*      3.5 Deardorff mixing length
 !           -----------------------
 !
   CASE ('DEAR')
-    CALL DEAR(PLEM)
+    CALL DEAR(ZLM)
 !
 !*      3.6 Blackadar mixing length
 !           -----------------------
 !
   CASE ('BLKR')
    ZL0 = 100.
-   PLEM(:,:,:) = ZL0
+   ZLM(:,:,:) = ZL0
 
    ZALPHA=0.5**(-1.5)
    !
    DO JK=IKTB,IKTE
-     PLEM(:,:,JK) = ( 0.5*(PZZ(:,:,JK)+PZZ(:,:,JK+KKL)) - &
+     ZLM(:,:,JK) = ( 0.5*(PZZ(:,:,JK)+PZZ(:,:,JK+KKL)) - &
      & PZZ(:,:,KKA+JPVEXT_TURB*KKL) ) * PDIRCOSZW(:,:)
-     PLEM(:,:,JK) = ZALPHA  * PLEM(:,:,JK) * ZL0 / ( ZL0 + ZALPHA*PLEM(:,:,JK) )
+     ZLM(:,:,JK) = ZALPHA  * ZLM(:,:,JK) * ZL0 / ( ZL0 + ZALPHA*ZLM(:,:,JK) )
    END DO
 !
-   PLEM(:,:,IKTB-1) = PLEM(:,:,IKTB)
-   PLEM(:,:,IKTE+1) = PLEM(:,:,IKTE)
+   ZLM(:,:,IKTB-1) = ZLM(:,:,IKTB)
+   ZLM(:,:,IKTE+1) = ZLM(:,:,IKTE)
 !
 !
 !
@@ -844,12 +683,17 @@ END SELECT
 !*      3.5 Mixing length modification for cloud
 !           -----------------------
 IF (KMODEL_CL==KMI .AND. HTURBLEN_CL/='NONE') CALL CLOUD_MODIF_LM
+ENDIF  ! end LHARRAT
 
 !
 !*      3.6 Dissipative length
 !           ------------------
-!
-ZLEPS(:,:,:)=PLEM(:,:,:)
+
+IF (LHARAT) THEN
+  ZLEPS=PLENGTHM*(3.75**2.)
+ELSE
+  ZLEPS=ZLM
+ENDIF
 !
 !*      3.7 Correction in the Surface Boundary Layer (Redelsperger 2001)
 !           ----------------------------------------
@@ -864,7 +708,7 @@ IF (ORMC01) THEN
     ZSFRV=0.
     ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRVM,PSFTH,ZSFRV)
   END IF
-  CALL RMC01(HTURBLEN,KKA,KKU,KKL,PZZ,PDXX,PDYY,PDZZ,PDIRCOSZW,PSBL_DEPTH,ZLMO,PLEM,ZLEPS)
+  CALL RMC01(HTURBLEN,KKA,KKU,KKL,PZZ,PDXX,PDYY,PDZZ,PDIRCOSZW,PSBL_DEPTH,ZLMO,ZLM,ZLEPS)
 END IF
 !
 !RMC01 is only applied on RM17 in ADAP
@@ -874,14 +718,14 @@ IF (HTURBLEN=='ADAP') ZLEPS = MIN(ZLEPS,ZLMW*XCADAP)
 !           ----------------------------------------------------------
 !
 IF (HTURBDIM=="3DIM") THEN
-  CALL UPDATE_LM(HLBCX,HLBCY,PLEM,ZLEPS)
+  CALL UPDATE_LM(HLBCX,HLBCY,ZLM,ZLEPS)
 END IF
 !
 !*      3.9 Mixing length correction if immersed walls 
 !           ------------------------------------------
 !
 IF (LIBM) THEN
-   CALL IBM_MIXINGLENGTH(PLEM,ZLEPS,XIBM_XMUT,XIBM_LS(:,:,:,1),PTKET)
+   CALL IBM_MIXINGLENGTH(ZLM,ZLEPS,XIBM_XMUT,XIBM_LS(:,:,:,1),PTKET)
 ENDIF
 !----------------------------------------------------------------------------
 !
@@ -968,39 +812,40 @@ ENDIF
 !*      5. TURBULENT SOURCES
 !          -----------------
 !
-if ( lbudget_u )  call Budget_store_init( tbudgets(NBUDGET_U ), 'VTURB', prus  (:, :, :)    )
-if ( lbudget_v )  call Budget_store_init( tbudgets(NBUDGET_V ), 'VTURB', prvs  (:, :, :)    )
-if ( lbudget_w )  call Budget_store_init( tbudgets(NBUDGET_W ), 'VTURB', prws  (:, :, :)    )
+IF( LBUDGET_U )  CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_U ), 'VTURB', PRUS(:, :, :)    )
+IF( LBUDGET_V )  CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_V ), 'VTURB', PRVS(:, :, :)    )
+IF( LBUDGET_W )  CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_W ), 'VTURB', PRWS(:, :, :)    )
 
-if ( lbudget_th ) then
-  if ( krri >= 1 .and. krrl >= 1 ) then
-    call Budget_store_init( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) &
-                                                                          + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) )
-  else if ( krrl >= 1 ) then
-    call Budget_store_init( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) )
-  else
-    call Budget_store_init( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) )
-  end if
-end if
+IF( LBUDGET_TH ) THEN
+  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) + ZLVOCPEXNM(:, :, :) * PRRS(:, :, :, 2) &
+                                                                          + ZLSOCPEXNM(:, :, :) * PRRS(:, :, :, 4) )
+  ELSE IF( KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) + ZLOCPEXNM(:, :, :) * PRRS(:, :, :, 2) )
+  ELSE
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) )
+  END IF
+END IF
+
+IF( LBUDGET_RV ) THEN
+  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) - PRRS(:, :, :, 4) )
+  ELSE IF( KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) )
+  ELSE
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:, :, :, 1) )
+  END IF
+END IF
 
-if ( lbudget_rv ) then
-  if ( krri >= 1 .and. krrl >= 1 ) then
-    call Budget_store_init( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) )
-  else if ( krrl >= 1 ) then
-    call Budget_store_init( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) )
-  else
-    call Budget_store_init( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) )
-  end if
-end if
+IF( LBUDGET_RC ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RC), 'VTURB', PRRS  (:, :, :, 2) )
+IF( LBUDGET_RI ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RI), 'VTURB', PRRS  (:, :, :, 4) )
 
-if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), 'VTURB', prrs  (:, :, :, 2) )
-if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI), 'VTURB', prrs  (:, :, :, 4) )
+IF( LBUDGET_SV ) THEN
+  DO JSV = 1, NSV
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'VTURB', PRSVS(:, :, :, JSV) )
+  END DO
+END IF
 
-if ( lbudget_sv ) then
-  do jsv = 1, nsv
-    call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + jsv), 'VTURB', prsvs(:, :, :, jsv) )
-  end do
-end if
 CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI,               &
           OTURB_FLX,                                     &
           HTURBDIM,HTOM,PIMPL,ZEXPL,                     &
@@ -1011,81 +856,84 @@ CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI,               &
           PSFTH,PSFRV,PSFSV,PSFTH,PSFRV,PSFSV,           &
           ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU33M,               &
           PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,    &
-          PTKET,PLEM,PLENGTHM,PLENGTHH,ZLEPS,MFMOIST,    &
+          PTKET,ZLM,PLENGTHM,PLENGTHH,ZLEPS,MFMOIST,     &
           ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE,     &
           ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,PBL_DEPTH,         &
           PSBL_DEPTH,ZLMO,                               &
           PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,              &
           PDYP,PTHP,PSIGS,PWTH,PWRC,PWSV                 )
 
-if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'VTURB', prus(:, :, :) )
-if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'VTURB', prvs(:, :, :) )
-if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'VTURB', prws(:, :, :) )
+IF( LBUDGET_U ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_U), 'VTURB', PRUS(:, :, :) )
+IF( LBUDGET_V ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_V), 'VTURB', PRVS(:, :, :) )
+IF( LBUDGET_W ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_W), 'VTURB', PRWS(:, :, :) )
 
-if ( lbudget_th ) then
-  if ( krri >= 1 .and. krrl >= 1 ) then
-    call Budget_store_end( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) &
-                                                                          + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) )
-  else if ( krrl >= 1 ) then
-    call Budget_store_end( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) )
-  else
-    call Budget_store_end( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) )
-  end if
-end if
+IF( LBUDGET_TH ) THEN
+  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) + ZLVOCPEXNM(:, :, :) * PRRS(:, :, :, 2) &
+                                                                          + ZLSOCPEXNM(:, :, :) * PRRS(:, :, :, 4) )
+  ELSE IF( KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) + ZLOCPEXNM(:, :, :) * PRRS(:, :, :, 2) )
+  ELSE
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) )
+  END IF
+END IF
 
-if ( lbudget_rv ) then
-  if ( krri >= 1 .and. krrl >= 1 ) then
-    call Budget_store_end( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) )
-  else if ( krrl >= 1 ) then
-    call Budget_store_end( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) )
-  else
-    call Budget_store_end( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) )
-  end if
-end if
+IF( LBUDGET_RV ) THEN
+  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) - PRRS(:, :, :, 4) )
+  ELSE IF( KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) )
+  ELSE
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:, :, :, 1) )
+  END IF
+END IF
 
-if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'VTURB', prrs(:, :, :, 2) )
-if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI), 'VTURB', prrs(:, :, :, 4) )
+IF( LBUDGET_RC ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RC), 'VTURB', PRRS(:, :, :, 2) )
+IF( LBUDGET_RI ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RI), 'VTURB', PRRS(:, :, :, 4) )
 
-if ( lbudget_sv )  then
-  do jsv = 1, nsv
-    call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + jsv), 'VTURB', prsvs(:, :, :, jsv) )
-  end do
-end if
-!
-if ( hturbdim == '3DIM' ) then
-  if ( lbudget_u  ) call Budget_store_init( tbudgets(NBUDGET_U ), 'HTURB', prus  (:, :, :) )
-  if ( lbudget_v  ) call Budget_store_init( tbudgets(NBUDGET_V ), 'HTURB', prvs  (:, :, :) )
-  if ( lbudget_w  ) call Budget_store_init( tbudgets(NBUDGET_W ), 'HTURB', prws  (:, :, :) )
+IF( LBUDGET_SV )  THEN
+  DO JSV = 1, NSV
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'VTURB', PRSVS(:, :, :, JSV) )
+  END DO
+END IF
+!
+#ifdef REPRO48
+#else
+IF( HTURBDIM == '3DIM' ) THEN
+#endif
+  IF( LBUDGET_U  ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_U ), 'HTURB', PRUS  (:, :, :) )
+  IF( LBUDGET_V  ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_V ), 'HTURB', PRVS  (:, :, :) )
+  IF( LBUDGET_W  ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_W ), 'HTURB', PRWS  (:, :, :) )
 
-  if (lbudget_th)  then
-    if ( krri >= 1 .and. krrl >= 1 ) then
-      call Budget_store_init( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) &
-                                                                             + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) )
-    else if ( krrl >= 1 ) then
-      call Budget_store_init( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) )
-    else
-      call Budget_store_init( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) )
-    end if
-  end if
+  IF(LBUDGET_TH)  THEN
+    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:, :, :) + ZLVOCPEXNM(:, :, :) * PRRS(:, :, :, 2) &
+                                                                             + ZLSOCPEXNM(:, :, :) * PRRS(:, :, :, 4) )
+    ELSE IF( KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:, :, :) + ZLOCPEXNM(:, :, :) * PRRS(:, :, :, 2) )
+    ELSE
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:, :, :) )
+    END IF
+  END IF
 
-  if ( lbudget_rv ) then
-    if ( krri >= 1 .and. krrl >= 1 ) then
-      call Budget_store_init( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) )
-    else if ( krrl >= 1 ) then
-      call Budget_store_init( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) )
-    else
-      call Budget_store_init( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) )
-    end if
-  end if
+  IF( LBUDGET_RV ) THEN
+    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) - PRRS(:, :, :, 4) )
+    ELSE IF( KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) )
+    ELSE
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:, :, :, 1) )
+    END IF
+  END IF
 
-  if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), 'HTURB', prrs(:, :, :, 2) )
-  if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI), 'HTURB', prrs(:, :, :, 4) )
+  IF( LBUDGET_RC ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:, :, :, 2) )
+  IF( LBUDGET_RI ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:, :, :, 4) )
 
-  if ( lbudget_sv )  then
-    do jsv = 1, nsv
-      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + jsv), 'HTURB', prsvs(:, :, :, jsv) )
-    end do
-  end if
+  IF( LBUDGET_SV )  THEN
+    DO JSV = 1, NSV
+      CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:, :, :, JSV) )
+    END DO
+  END IF
 
     CALL TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP,        &
           HLBCX,HLBCY,OTURB_FLX,OSUBG_COND,                    &
@@ -1097,46 +945,49 @@ if ( hturbdim == '3DIM' ) then
           PSFTH,PSFRV,PSFSV,                                   &
           ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU22M,ZTAU33M,             &
           PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,          &
-          PTKET,PLEM,ZLEPS,                                    &
+          PTKET,ZLM,ZLEPS,                                    &
           ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE,           &
           PDYP,PTHP,PSIGS,                                     &
           ZTRH,                                                &
           PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS                     )
 
-  if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'HTURB', prus(:, :, :) )
-  if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'HTURB', prvs(:, :, :) )
-  if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'HTURB', prws(:, :, :) )
+  IF( LBUDGET_U ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_U), 'HTURB', PRUS(:, :, :) )
+  IF( LBUDGET_V ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_V), 'HTURB', PRVS(:, :, :) )
+  IF( LBUDGET_W ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_W), 'HTURB', PRWS(:, :, :) )
 
-  if ( lbudget_th ) then
-    if ( krri >= 1 .and. krrl >= 1 ) then
-      call Budget_store_end( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) &
-                                                                            + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) )
-    else if ( krrl >= 1 ) then
-      call Budget_store_end( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) )
-    else
-      call Budget_store_end( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) )
-    end if
-  end if
+  IF( LBUDGET_TH ) THEN
+    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:, :, :) + ZLVOCPEXNM(:, :, :) * PRRS(:, :, :, 2) &
+                                                                            + ZLSOCPEXNM(:, :, :) * PRRS(:, :, :, 4) )
+    ELSE IF( KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:, :, :) + ZLOCPEXNM(:, :, :) * PRRS(:, :, :, 2) )
+    ELSE
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:, :, :) )
+    END IF
+  END IF
 
-  if ( lbudget_rv ) then
-    if ( krri >= 1 .and. krrl >= 1 ) then
-      call Budget_store_end( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) )
-    else if ( krrl >= 1 ) then
-      call Budget_store_end( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) )
-    else
-      call Budget_store_end( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) )
-    end if
-  end if
+  IF( LBUDGET_RV ) THEN
+    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) - PRRS(:, :, :, 4) )
+    ELSE IF( KRRL >= 1 ) THEN
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:, :, :, 1) - PRRS(:, :, :, 2) )
+    ELSE
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:, :, :, 1) )
+    END IF
+  END IF
 
-  if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'HTURB', prrs(:, :, :, 2) )
-  if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI), 'HTURB', prrs(:, :, :, 4) )
+  IF( LBUDGET_RC ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:, :, :, 2) )
+  IF( LBUDGET_RI ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:, :, :, 4) )
 
-  if ( lbudget_sv )  then
-    do jsv = 1, nsv
-      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + jsv), 'HTURB', prsvs(:, :, :, jsv) )
-    end do
-  end if
-end if
+  IF( LBUDGET_SV )  THEN
+    DO JSV = 1, NSV
+      CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:, :, :, JSV) )
+    END DO
+  END IF
+#ifdef REPRO48
+#else
+END IF
+#endif
 !----------------------------------------------------------------------------
 !
 !*      6. EVOLUTION OF THE TKE AND ITS DISSIPATION 
@@ -1149,7 +1000,19 @@ end if
 
 !  6.2 TKE evolution equation
 
-CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,PLEM,ZLEPS,PDYP,ZTRH,     &
+IF (.NOT. LHARAT) THEN
+!
+IF (LBUDGET_TH)  THEN
+  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
+                                                          & + ZLSOCPEXNM * PRRS(:,:,:,4) )
+  ELSE IF ( KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
+  ELSE
+    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
+  END IF
+END IF
+CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,ZLM,ZLEPS,PDYP,ZTRH,     &
                      PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,            &
                      PTSTEP,PIMPL,ZEXPL,                             &
                      HTURBLEN,HTURBDIM,                              &
@@ -1157,6 +1020,19 @@ CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,PLEM,ZLEPS,PDYP,ZTRH,     &
                      PTHP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,&
                      TBUDGETS,KBUDGETS,PRTKESM=PRTKEMS)
 
+IF (LBUDGET_TH)  THEN
+  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
+                                                          & + ZLSOCPEXNM * PRRS(:,:,:,4) )
+  ELSE IF ( KRRL >= 1 ) THEN
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
+  ELSE
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
+  END IF
+END IF
+!
+ENDIF
+!
 !----------------------------------------------------------------------------
 !
 !*      7. STORES SOME INFORMATIONS RELATED TO THE TURBULENCE SCHEME
@@ -1176,7 +1052,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   TZFIELD%NTYPE      = TYPEREAL
   TZFIELD%NDIMS      = 3
   TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PLEM)
+  CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM)
 !
   IF (KRR /= 0) THEN
 !
@@ -1192,7 +1068,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 3
     TZFIELD%LTIMEDEP   = .TRUE.
-    CALL IO_Field_write(TPFILE,TZFIELD,PTHLT)
+    CALL IO_FIELD_WRITE(TPFILE,TZFIELD,PTHLT)
 !
 ! stores the conservative mixing ratio
 !
@@ -1206,7 +1082,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 3
     TZFIELD%LTIMEDEP   = .TRUE.
-    CALL IO_Field_write(TPFILE,TZFIELD,PRT(:,:,:,1))
+    CALL IO_FIELD_WRITE(TPFILE,TZFIELD,PRT(:,:,:,1))
    END IF
 END IF
 !
@@ -1290,7 +1166,7 @@ IF (LLES_CALL) THEN
 !*     12. LES mixing end dissipative lengths, presso-correlations
 !          -------------------------------------------------------
 !
-  CALL LES_MEAN_SUBGRID(PLEM,X_LES_SUBGRID_LMix)
+  CALL LES_MEAN_SUBGRID(ZLM,X_LES_SUBGRID_LMix)
   CALL LES_MEAN_SUBGRID(ZLEPS,X_LES_SUBGRID_LDiss)
 !
 !* presso-correlations for subgrid Tke are equal to zero.
@@ -1302,8 +1178,10 @@ IF (LLES_CALL) THEN
   XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
 END IF
 !
+IF(PRESENT(PLEM)) PLEM = ZLM
 !----------------------------------------------------------------------------
 !
+IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
 CONTAINS
 !
 !
@@ -1503,6 +1381,8 @@ REAL                :: ZD           ! distance to the surface
 !
 !-------------------------------------------------------------------------------
 !
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB:DELT',0,ZHOOK_HANDLE)
 IF (ODZ) THEN
   ! Dz is take into account in the computation
   DO JK = IKTB,IKTE ! 1D turbulence scheme
@@ -1564,6 +1444,7 @@ END IF
 PLM(:,:,KKA) = PLM(:,:,IKB  )
 PLM(:,:,KKU  ) = PLM(:,:,IKE)
 !
+IF (LHOOK) CALL DR_HOOK('TURB:DELT',1,ZHOOK_HANDLE)
 END SUBROUTINE DELT
 !
 !     ####################
@@ -1607,6 +1488,8 @@ REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) ::     &
 !-------------------------------------------------------------------------------
 !
 !   initialize the mixing length with the mesh grid
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB:DEAR',0,ZHOOK_HANDLE)
 ! 1D turbulence scheme
 PLM(:,:,IKTB:IKTE) = PZZ(:,:,IKTB+KKL:IKTE+KKL) - PZZ(:,:,IKTB:IKTE)
 PLM(:,:,KKU) = PLM(:,:,IKE)
@@ -1720,6 +1603,7 @@ PLM(:,:,KKA) = PLM(:,:,IKB  )
 PLM(:,:,IKE  ) = PLM(:,:,IKE-KKL)
 PLM(:,:,KKU  ) = PLM(:,:,KKU-KKL)
 !
+IF (LHOOK) CALL DR_HOOK('TURB:DEAR',1,ZHOOK_HANDLE)
 END SUBROUTINE DEAR
 !
 !     #########################
@@ -1787,6 +1671,8 @@ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZLM_CLOUD
 !*       1.    INITIALISATION
 !              --------------
 !
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',0,ZHOOK_HANDLE)
 ZPENTE = ( PCOEF_AMPL_SAT - 1. ) / ( PCEI_MAX - PCEI_MIN ) 
 ZCOEF_AMPL_CEI_NUL = 1. - ZPENTE * PCEI_MIN
 !
@@ -1811,7 +1697,7 @@ WHERE ( PCEI(:,:,:) <  PCEI_MAX .AND.                                        &
 !              ------------------------------------------
 !
 IF (HTURBLEN_CL == HTURBLEN) THEN
-  ZLM_CLOUD(:,:,:) = PLEM(:,:,:)
+  ZLM_CLOUD(:,:,:) = ZLM(:,:,:)
 ELSE
   SELECT CASE (HTURBLEN_CL)
 !
@@ -1849,16 +1735,16 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   TZFIELD%NTYPE      = TYPEREAL
   TZFIELD%NDIMS      = 3
   TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PLEM)
+  CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM)
 ENDIF
 !
 ! Amplification of the mixing length when the criteria are verified
 !
-WHERE (ZCOEF_AMPL(:,:,:) /= 1.) PLEM(:,:,:) = ZCOEF_AMPL(:,:,:)*ZLM_CLOUD(:,:,:)
+WHERE (ZCOEF_AMPL(:,:,:) /= 1.) ZLM(:,:,:) = ZCOEF_AMPL(:,:,:)*ZLM_CLOUD(:,:,:)
 !
 ! Cloud mixing length in the clouds at the points which do not verified the CEI
 !
-WHERE (PCEI(:,:,:) == -1.) PLEM(:,:,:) = ZLM_CLOUD(:,:,:)
+WHERE (PCEI(:,:,:) == -1.) ZLM(:,:,:) = ZLM_CLOUD(:,:,:)
 !
 !
 !*       5.    IMPRESSION
@@ -1875,7 +1761,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   TZFIELD%NTYPE      = TYPEREAL
   TZFIELD%NDIMS      = 3
   TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,ZCOEF_AMPL)
+  CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZCOEF_AMPL)
   !
   TZFIELD%CMNHNAME   = 'LM_CLOUD'
   TZFIELD%CSTDNAME   = ''
@@ -1886,10 +1772,11 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   TZFIELD%NGRID      = 1
   TZFIELD%NTYPE      = TYPEREAL
   TZFIELD%NDIMS      = 3
-  CALL IO_Field_write(TPFILE,TZFIELD,ZLM_CLOUD)
+  CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM_CLOUD)
   !
 ENDIF
 !
+IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',1,ZHOOK_HANDLE)
 END SUBROUTINE CLOUD_MODIF_LM
 !
 END SUBROUTINE TURB